/[escript]/trunk/finley/src/Mesh_hex20.c
ViewVC logotype

Contents of /trunk/finley/src/Mesh_hex20.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1062 - (show annotations)
Mon Mar 26 06:17:53 2007 UTC (12 years, 7 months ago) by gross
File MIME type: text/plain
File size: 48608 byte(s)
reduced integration schemes are implemented now for grad, integrate, etc. Tests still to be added.
1 /*
2 ************************************************************
3 * Copyright 2006 by ACcESS MNRF *
4 * *
5 * http://www.access.edu.au *
6 * Primary Business: Queensland, Australia *
7 * Licensed under the Open Software License version 3.0 *
8 * http://www.opensource.org/licenses/osl-3.0.php *
9 * *
10 ************************************************************
11 */
12
13 /**************************************************************/
14
15 /* Finley: generates rectangular meshes */
16
17 /* Generates a numElements[0] x numElements[1] x numElements[2] mesh with second order elements (Hex20) in the brick */
18 /* [0,Length[0]] x [0,Length[1]] x [0,Length[2]]. order is the desired accuracy of the */
19 /* integration scheme. */
20
21
22 /**************************************************************/
23
24 /* Author: gross@access.edu.au */
25 /* Version: $Id$
26
27 /**************************************************************/
28
29 #include "RectangularMesh.h"
30
31 #ifdef PASO_MPI
32 /* get the number of nodes/elements for domain with rank=rank, of size processors
33 where n is the total number of nodes/elements in the global domain */
34 static index_t domain_MODdim( index_t rank, index_t size, index_t n )
35 {
36 rank = size-rank-1;
37
38 if( rank < n%size )
39 return (index_t)floor(n/size)+1;
40 return (index_t)floor(n/size);
41 }
42
43
44 /* Determines the number of nodes/elements etc along an axis which is numElementsGlobal long for domain rank */
45 /* A bit messy, but it only has to be done once... */
46 static void domain_calculateDimension( index_t rank, dim_t size, dim_t numElementsGlobal, bool_t periodic, dim_t *numNodesLocal, dim_t *numDOFLocal, dim_t *numElementsLocal, dim_t *numElementsInternal, dim_t *firstNode, dim_t *nodesExternal, dim_t *DOFExternal, dim_t *numNodesExternal, bool_t *periodicLocal, dim_t *DOFBoundary )
47 {
48 index_t i0;
49 dim_t numNodesGlobal = numElementsGlobal+1;
50
51 (*numNodesLocal) = 2*domain_MODdim( rank, size, numNodesGlobal ) - 1;
52 if( rank<size-1 ) // add on node for right hand boundary
53 (*numNodesLocal) += 1;
54
55 numElementsLocal[0] = domain_MODdim( rank, size, numNodesGlobal )+1;
56 periodicLocal[0] = periodicLocal[1] = FALSE;
57 nodesExternal[0] = nodesExternal[1] = 1;
58 if( periodic )
59 {
60 if( size==1 )
61 {
62 numElementsLocal[0] = numElementsGlobal;
63 nodesExternal[0] = nodesExternal[1] = 0;
64 periodicLocal[0] = periodicLocal[1] = TRUE;
65 }
66 else
67 {
68 if( rank==0 )
69 {
70 periodicLocal[0] = TRUE;
71 numNodesLocal[0]++;
72 }
73 if( rank==(size-1) )
74 {
75 periodicLocal[1] = TRUE;
76 numNodesLocal[0]--;
77 numElementsLocal[0]--;
78 }
79 }
80 }
81 else if( !periodic )
82 {
83 if( rank==0 ){
84 nodesExternal[0]--;
85 numElementsLocal[0]--;
86 }
87 if( rank==(size-1) )
88 {
89 nodesExternal[1]--;
90 numElementsLocal[0]--;
91 }
92 }
93 nodesExternal[0]*=2;
94 numNodesExternal[0] = nodesExternal[0]+nodesExternal[1];
95
96 numElementsInternal[0] = numElementsLocal[0];
97 if( (rank==0) && (rank==size-1) );
98
99 else if( !periodic && ( (rank==0) ^ (rank==size-1) ) )
100 numElementsInternal[0] -= 1;
101 else
102 numElementsInternal[0] -= 2;
103
104 firstNode[0] = 0;
105 for( i0=0; i0<rank; i0++ )
106 firstNode[0] += domain_MODdim( i0, size, numNodesGlobal );
107 firstNode[0] *= 2;
108
109 numDOFLocal[0] = numNodesLocal[0];
110 if( periodicLocal[0] )
111 numDOFLocal[0]--;
112
113 DOFExternal[0] = nodesExternal[0];
114 DOFExternal[1] = nodesExternal[1];
115
116 if( size>1 )
117 {
118 DOFBoundary[0] = periodicLocal[0]*1 + (rank>0)*1;
119 DOFBoundary[1] = periodicLocal[1]*2 + (rank<(size-1))*2;
120 }
121 else
122 {
123 DOFBoundary[0] = DOFBoundary[1] = 0;
124 }
125 }
126 #endif
127 /**************************************************************/
128 #ifdef PASO_MPI
129 Finley_Mesh* Finley_RectangularMesh_Hex20_singleCPU(dim_t* numElements,double* Length,bool_t* periodic,index_t order, index_t reduced_order, bool_t useElementsOnFace, Paso_MPIInfo *mpi_info) {
130 #else
131 Finley_Mesh* Finley_RectangularMesh_Hex20(dim_t* numElements,double* Length,bool_t* periodic,index_t order, index_t reduced_order, bool_t useElementsOnFace) {
132 #endif
133 dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES,M0,M1,M2;
134 index_t node0;
135 Finley_Mesh* out;
136 char name[50];
137 double time0=Finley_timer();
138
139 NE0=MAX(1,numElements[0]);
140 NE1=MAX(1,numElements[1]);
141 NE2=MAX(1,numElements[2]);
142 N0=2*NE0+1;
143 N1=2*NE1+1;
144 N2=2*NE2+1;
145
146 if (N0<=MIN(N1,N2)) {
147 if (N1 <= N2) {
148 M0=1;
149 M1=N0;
150 M2=N0*N1;
151 } else {
152 M0=1;
153 M2=N0;
154 M1=N0*N2;
155 }
156 } else if (N1<=MIN(N2,N0)) {
157 if (N2 <= N0) {
158 M1=1;
159 M2=N1;
160 M0=N2*N1;
161 } else {
162 M1=1;
163 M0=N1;
164 M2=N1*N0;
165 }
166 } else {
167 if (N0 <= N1) {
168 M2=1;
169 M0=N2;
170 M1=N2*N0;
171 } else {
172 M2=1;
173 M1=N2;
174 M0=N1*N2;
175 }
176 }
177
178 NFaceElements=0;
179 if (!periodic[0]) {
180 NDOF0=N0;
181 NFaceElements+=2*NE1*NE2;
182 } else {
183 NDOF0=N0-1;
184 }
185 if (!periodic[1]) {
186 NDOF1=N1;
187 NFaceElements+=2*NE0*NE2;
188 } else {
189 NDOF1=N1-1;
190 }
191 if (!periodic[2]) {
192 NDOF2=N2;
193 NFaceElements+=2*NE0*NE1;
194 } else {
195 NDOF2=N2-1;
196 }
197
198 /* allocate mesh: */
199
200 sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
201 #ifdef PASO_MPI
202 out=Finley_Mesh_alloc(name,3,order, reduced_order, mpi_info);
203
204 if (! Finley_noError()) return NULL;
205
206 out->Elements=Finley_ElementFile_alloc(Hex20,out->order,out->reduced_order,mpi_info);
207 if (useElementsOnFace) {
208 out->FaceElements=Finley_ElementFile_alloc(Hex20Face,out->order,out->reduced_order,mpi_info);
209 out->ContactElements=Finley_ElementFile_alloc(Hex20Face_Contact,out->order,out->reduced_order,mpi_info);
210 } else {
211 out->FaceElements=Finley_ElementFile_alloc(Rec8,out->order,out->reduced_order,mpi_info);
212 out->ContactElements=Finley_ElementFile_alloc(Rec8_Contact,out->order,out->reduced_order,mpi_info);
213 }
214 out->Points=Finley_ElementFile_alloc(Point1,out->order,out->reduced_order,mpi_info);
215 if (! Finley_noError()) {
216 Finley_Mesh_dealloc(out);
217 return NULL;
218 }
219
220
221 /* allocate tables: */
222 Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);
223 Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);
224 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
225 Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0*NDOF1*NDOF2, 0, 0 );
226 Finley_ElementDistribution_allocTable( out->Elements->elementDistribution, NE0*NE1*NE2, NE0*NE1*NE2);
227 Finley_ElementDistribution_allocTable( out->FaceElements->elementDistribution, NFaceElements, NFaceElements );
228 #else
229 out=Finley_Mesh_alloc(name,3,order, reduced_order);
230
231 if (! Finley_noError()) return NULL;
232
233 out->Elements=Finley_ElementFile_alloc(Hex20,out->order,out->reduced_order);
234 if (useElementsOnFace) {
235 out->FaceElements=Finley_ElementFile_alloc(Hex20Face,out->order,out->reduced_order);
236 out->ContactElements=Finley_ElementFile_alloc(Hex20Face_Contact,out->order,out->reduced_order);
237 } else {
238 out->FaceElements=Finley_ElementFile_alloc(Rec8,out->order,out->reduced_order);
239 out->ContactElements=Finley_ElementFile_alloc(Rec8_Contact,out->order,out->reduced_order);
240 }
241 out->Points=Finley_ElementFile_alloc(Point1,out->order,out->reduced_order);
242 if (! Finley_noError()) {
243 Finley_Mesh_dealloc(out);
244 return NULL;
245 }
246
247
248 /* allocate tables: */
249 Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);
250 Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);
251 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
252 #endif
253 if (! Finley_noError()) {
254 Finley_Mesh_dealloc(out);
255 return NULL;
256 }
257
258 /* create nodes */
259 #pragma omp parallel for private(i0,i1,i2,k)
260 for (i2=0;i2<N2;i2++) {
261 for (i1=0;i1<N1;i1++) {
262 for (i0=0;i0<N0;i0++) {
263 k=M0*i0+M1*i1+M2*i2;
264 out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0)/DBLE(N0-1)*Length[0];
265 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
266 out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
267 out->Nodes->Id[k]=i0+N0*i1+N0*N1*i2;
268 out->Nodes->Tag[k]=0;
269 out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1) +M2*(i2%NDOF2);
270 #ifdef PASO_MPI
271 out->Nodes->Dom[k]=NODE_INTERNAL;
272 #endif
273 }
274 }
275 }
276
277
278 /* tags for the faces: */
279 if (!periodic[2]) {
280 for (i1=0;i1<N1;i1++) {
281 for (i0=0;i0<N0;i0++) {
282 out->Nodes->Tag[M0*i0+M1*i1+M2*0]+=100;
283 out->Nodes->Tag[M0*i0+M1*i1+M2*(N2-1)]+=200;
284 }
285 }
286 }
287 if (!periodic[1]) {
288 for (i2=0;i2<N2;i2++) {
289 for (i0=0;i0<N0;i0++) {
290 out->Nodes->Tag[M0*i0+M1*0+M2*i2]+=10;
291 out->Nodes->Tag[M0*i0+M1*(N1-1)+M2*i2]+=20;
292 }
293 }
294 }
295 if (!periodic[0]) {
296 for (i2=0;i2<N2;i2++) {
297 for (i1=0;i1<N1;i1++) {
298 out->Nodes->Tag[M0*0+M1*i1+M2*i2]+=1;
299 out->Nodes->Tag[M0*(N0-1)+M1*i1+M2*i2]+=2;
300 }
301 }
302 }
303
304 /* set the elements: */
305
306 #pragma omp parallel for private(i0,i1,i2,k,node0)
307 for (i2=0;i2<NE2;i2++) {
308 for (i1=0;i1<NE1;i1++) {
309 for (i0=0;i0<NE0;i0++) {
310 k=i0+NE0*i1+NE0*NE1*i2;
311 node0=2*i0+2*i1*N0+2*N0*N1*i2;
312
313 out->Elements->Id[k]=k;
314 out->Elements->Tag[k]=0;
315 out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);
316 #ifdef PASO_MPI
317 out->Elements->Dom[k]=ELEMENT_INTERNAL;
318 #endif
319
320 out->Elements->Nodes[INDEX2(0,k,20)]=node0;
321 out->Elements->Nodes[INDEX2(1,k,20)]=node0+2;
322 out->Elements->Nodes[INDEX2(2,k,20)]=node0+2*N0+2;
323 out->Elements->Nodes[INDEX2(3,k,20)]=node0+2*N0;
324 out->Elements->Nodes[INDEX2(4,k,20)]=node0+2*N0*N1;
325 out->Elements->Nodes[INDEX2(5,k,20)]=node0+2*N0*N1+2;
326 out->Elements->Nodes[INDEX2(6,k,20)]=node0+2*N0*N1+2*N0+2;
327 out->Elements->Nodes[INDEX2(7,k,20)]=node0+2*N0*N1+2*N0;
328 out->Elements->Nodes[INDEX2(8,k,20)]=node0+1;
329 out->Elements->Nodes[INDEX2(9,k,20)]=node0+N0+2;
330 out->Elements->Nodes[INDEX2(10,k,20)]=node0+2*N0+1;
331 out->Elements->Nodes[INDEX2(11,k,20)]=node0+N0;
332 out->Elements->Nodes[INDEX2(12,k,20)]=node0+N0*N1;
333 out->Elements->Nodes[INDEX2(13,k,20)]=node0+N0*N1+2;
334 out->Elements->Nodes[INDEX2(14,k,20)]=node0+N0*N1+2*N0+2;
335 out->Elements->Nodes[INDEX2(15,k,20)]=node0+N0*N1+2*N0;
336 out->Elements->Nodes[INDEX2(16,k,20)]=node0+2*N0*N1+1;
337 out->Elements->Nodes[INDEX2(17,k,20)]=node0+2*N0*N1+N0+2;
338 out->Elements->Nodes[INDEX2(18,k,20)]=node0+2*N0*N1+2*N0+1;
339 out->Elements->Nodes[INDEX2(19,k,20)]=node0+2*N0*N1+N0;
340 }
341 }
342 }
343 out->Elements->minColor=0;
344 out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
345
346 /* face elements: */
347
348 if (useElementsOnFace) {
349 NUMNODES=20;
350 } else {
351 NUMNODES=8;
352 }
353 totalNECount=NE0*NE1*NE2;
354 faceNECount=0;
355
356 /* these are the quadrilateral elements on boundary 1 (x3=0): */
357
358 if (!periodic[2]) {
359 /* ** elements on boundary 100 (x3=0): */
360 #pragma omp parallel for private(i0,i1,k,node0)
361 for (i1=0;i1<NE1;i1++) {
362 for (i0=0;i0<NE0;i0++) {
363 k=i0+NE0*i1+faceNECount;
364 node0=2*i0+2*i1*N0;
365
366 out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
367 out->FaceElements->Tag[k]=100;
368 out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
369 #ifdef PASO_MPI
370 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
371 #endif
372
373 if (useElementsOnFace) {
374 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
375 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;
376 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0+2;
377 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2;
378 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2*N0*N1;
379 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1+2*N0;
380 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
381 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0*N1+2;
382 out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+N0;
383 out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+2*N0+1;
384 out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+N0+2;
385 out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+1;
386 out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+N0*N1;
387 out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+N0*N1+2*N0;
388 out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+N0*N1+2*N0+2;
389 out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+N0*N1+2;
390 out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+2*N0*N1+N0;
391 out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;
392 out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+2*N0*N1+N0+2;
393 out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+2*N0*N1+1;
394 } else {
395 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
396 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;
397 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0+2;
398 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2;
399 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;
400 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0+1;
401 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+2;
402 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
403 }
404 }
405 }
406 totalNECount+=NE1*NE0;
407 faceNECount+=NE1*NE0;
408
409 /* ** elements on boundary 200 (x3=1) */
410 #pragma omp parallel for private(i0,i1,k,node0)
411 for (i1=0;i1<NE1;i1++) {
412 for (i0=0;i0<NE0;i0++) {
413 k=i0+NE0*i1+faceNECount;
414 node0=2*i0+2*i1*N0+2*N0*N1*(NE2-1);
415
416 out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
417 out->FaceElements->Tag[k]=200;
418 out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
419 #ifdef PASO_MPI
420 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
421 #endif
422
423 if (useElementsOnFace) {
424 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0*N1;
425 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1+2;
426 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
427 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0*N1+2*N0;
428
429 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
430 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2;
431 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0+2;
432 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0;
433
434 out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+2*N0*N1+1;
435 out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+2*N0*N1+N0+2;
436 out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;
437 out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+2*N0*N1+N0;
438
439 out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+N0*N1;
440 out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+N0*N1+2;
441 out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+N0*N1+2*N0+2;
442 out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+N0*N1+2*N0;
443
444 out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+1;
445 out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+N0+2;
446 out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+2*N0+1;
447 out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+N0;
448
449 } else {
450 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0*N1;
451 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1+2;
452 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
453 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0*N1+2*N0;
454 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2*N0*N1+1;
455 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1+N0+2;
456 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;
457 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0*N1+N0;
458 }
459 }
460 }
461 totalNECount+=NE1*NE0;
462 faceNECount+=NE1*NE0;
463 }
464 if (!periodic[0]) {
465 /* ** elements on boundary 001 (x1=0): */
466
467 #pragma omp parallel for private(i1,i2,k,node0)
468 for (i2=0;i2<NE2;i2++) {
469 for (i1=0;i1<NE1;i1++) {
470 k=i1+NE1*i2+faceNECount;
471 node0=2*i1*N0+2*N0*N1*i2;
472
473 out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
474 out->FaceElements->Tag[k]=1;
475 out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
476 #ifdef PASO_MPI
477 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
478 #endif
479
480 if (useElementsOnFace) {
481 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
482 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1;
483 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0;
484 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0;
485
486 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2;
487 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1+2;
488 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
489 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0+2;
490
491 out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+N0*N1;
492 out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+2*N0*N1+N0;
493 out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+N0*N1+2*N0;
494 out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+N0;
495
496 out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+1;
497 out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+2*N0*N1+1;
498 out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;
499 out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+2*N0+1;
500
501 out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+N0*N1+2;
502 out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+2*N0*N1+N0+2;
503 out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+N0*N1+2*N0+2;
504 out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+N0+2;
505
506 } else {
507 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
508 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1;
509 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0;
510 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0;
511 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0*N1;
512 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1+N0;
513 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+2*N0;
514 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;
515 }
516 }
517 }
518 totalNECount+=NE1*NE2;
519 faceNECount+=NE1*NE2;
520
521 /* ** elements on boundary 002 (x1=1): */
522
523 #pragma omp parallel for private(i1,i2,k,node0)
524 for (i2=0;i2<NE2;i2++) {
525 for (i1=0;i1<NE1;i1++) {
526 k=i1+NE1*i2+faceNECount;
527 node0=2*(NE0-1)+2*i1*N0+2*N0*N1*i2 ;
528
529 out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
530 out->FaceElements->Tag[k]=2;
531 out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
532 #ifdef PASO_MPI
533 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
534 #endif
535
536 if (useElementsOnFace) {
537 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;
538 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0+2;
539 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
540 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0*N1+2;
541
542 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
543 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0;
544 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2*N0;
545 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0*N1;
546
547 out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+N0+2;
548 out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+N0*N1+2*N0+2;
549 out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+2*N0*N1+N0+2;
550 out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+N0*N1+2;
551
552 out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+1;
553 out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+2*N0+1;
554 out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;
555 out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+2*N0*N1+1;
556
557 out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+N0;
558 out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+N0*N1+2*N0;
559 out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+2*N0*N1+N0;
560 out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+N0*N1;
561
562 } else {
563 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;
564 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0+2;
565 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
566 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0*N1+2;
567 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0+2;
568 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+2*N0+2;
569 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+N0+2;
570 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1+2;
571 }
572
573 }
574 }
575 totalNECount+=NE1*NE2;
576 faceNECount+=NE1*NE2;
577 }
578 if (!periodic[1]) {
579 /* ** elements on boundary 010 (x2=0): */
580
581 #pragma omp parallel for private(i0,i2,k,node0)
582 for (i2=0;i2<NE2;i2++) {
583 for (i0=0;i0<NE0;i0++) {
584 k=i0+NE0*i2+faceNECount;
585 node0=2*i0+2*N0*N1*i2;
586
587 out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
588 out->FaceElements->Tag[k]=10;
589 out->FaceElements->Color[k]=(i2%2)+2*(i0%2)+16;
590 #ifdef PASO_MPI
591 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
592 #endif
593
594 if (useElementsOnFace) {
595 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
596 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2;
597 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N1*N0+2;
598 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N1*N0;
599
600 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2*N0;
601 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0+2;
602 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N1*N0+2*N0+2;
603 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N1*N0+2*N0;
604
605 out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+1;
606 out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+N0*N1+2;
607 out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+2*N1*N0+1;
608 out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+N1*N0;
609
610 out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+N0;
611 out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+N0+2;
612 out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+2*N1*N0+N0+2;
613 out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+2*N1*N0+N0;
614
615 out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+2*N0+1;
616 out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+N0*N1+2*N0+2;
617 out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+2*N1*N0+2*N0+1;
618 out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+N1*N0+2*N0;
619
620 } else {
621 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
622 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2;
623 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N1*N0+2;
624 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N1*N0;
625 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
626 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+2;
627 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N1*N0+1;
628 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*N0;
629 }
630 }
631 }
632 totalNECount+=NE0*NE2;
633 faceNECount+=NE0*NE2;
634
635 /* ** elements on boundary 020 (x2=1): */
636
637 #pragma omp parallel for private(i0,i2,k,node0)
638 for (i2=0;i2<NE2;i2++) {
639 for (i0=0;i0<NE0;i0++) {
640 k=i0+NE0*i2+faceNECount;
641 node0=2*i0+2*(NE1-1)*N0+2*N0*N1*i2;
642
643 out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
644 out->FaceElements->Tag[k]=20;
645 out->FaceElements->Color[k]=(i2%2)+2*(i0%2)+20;
646 #ifdef PASO_MPI
647 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
648 #endif
649
650 if (useElementsOnFace) {
651 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;
652 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1+2*N0;
653 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
654 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0+2;
655
656 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
657 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1;
658 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2;
659 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2;
660
661 out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+N1*N0+2*N0;
662 out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+2*N1*N0+2*N0+1;
663 out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+N0*N1+2*N0+2;
664 out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+2*N0+1;
665
666 out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+N0;
667 out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+2*N0*N1+N0;
668 out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+2*N0*N1+N0+2;
669 out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+N0+2;
670
671 out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+N1*N0;
672 out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+2*N1*N0+1;
673 out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+N0*N1+2;
674 out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+1;
675 } else {
676 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;
677 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1+2*N0;
678 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
679 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0+2;
680 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N1*N0+2*N0;
681 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N1*N0+2*N0+1;
682 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+2*N0+2;
683 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0+1;
684 }
685 }
686 }
687 totalNECount+=NE0*NE2;
688 faceNECount+=NE0*NE2;
689 }
690 out->FaceElements->minColor=0;
691 out->FaceElements->maxColor=24;
692
693 #ifdef PASO_MPI
694 Finley_ElementFile_setDomainFlags( out->Elements );
695 Finley_ElementFile_setDomainFlags( out->FaceElements );
696 Finley_ElementFile_setDomainFlags( out->ContactElements );
697 Finley_ElementFile_setDomainFlags( out->Points );
698
699 /* reorder the degrees of freedom */
700 Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
701 #endif
702
703 /* condense the nodes: */
704 Finley_Mesh_resolveNodeIds(out);
705
706 /* prepare mesh for further calculatuions:*/
707 Finley_Mesh_prepare(out) ;
708
709 #ifdef Finley_TRACE
710 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
711 #endif
712 if (Finley_noError()) {
713 if (! Finley_Mesh_isPrepared(out) ) {
714 Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");
715 }
716 } else {
717 Finley_Mesh_dealloc(out);
718 }
719 return out;
720 }
721 #ifdef PASO_MPI
722 Finley_Mesh* Finley_RectangularMesh_Hex20(dim_t* numElements,double* Length,bool_t* periodic,index_t order, index_t reduced_order, bool_t useElementsOnFace) {
723 dim_t N0,N0t,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,totalNECount,faceNECount,NDOF0,NDOF0t,NDOF1,NDOF2,NFaceElements,NUMNODES,M0,M1,M2;
724 dim_t kk,iI, NE0_local, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal, DOFBoundary[2];
725 bool_t dom_left, dom_right, dom_internal;
726 index_t firstNode=0, DOFcount=0, node0, node1, node2, idCount;
727 index_t targetDomain=-1, firstNodeConstruct, j;
728 bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE, boundaryLeft=FALSE, boundaryRight=FALSE;
729 index_t *indexBackward=NULL, *indexForward=NULL,*facePerm=NULL, *forwardDOF=NULL, *backwardDOF=NULL;
730 Finley_Mesh* out;
731 char name[50];
732 double time0=Finley_timer();
733 Paso_MPIInfo *mpi_info = NULL;
734
735 NE0=MAX(1,numElements[0]);
736 NE1=MAX(1,numElements[1]);
737 NE2=MAX(1,numElements[2]);
738 N0=2*NE0+1;
739 N1=2*NE1+1;
740 N2=2*NE2+1;
741
742 index_t face0[] = {0,4,7,3,1,5,6,2,12,19,15,11,8,16,18,10,13,17,14,9};
743 index_t face1[] = {1,2,6,5,0,3,7,4,9,14,17,13,8,10,18,16,11,15,19,12};
744 index_t face2[] = {0,1,5,4,3,2,6,7,8,13,16,12,11,9,17,19,14,18,15,10};
745 index_t face3[] = {3,7,6,2,0,4,5,1,15,18,14,10,11,19,17,9,12,16,13,8};
746 index_t face4[] = {0,3,2,1,4,7,6,5,11,10,9,8,15,14,13,12,19,18,17,16};
747 index_t face5[] = {4,5,6,7,0,1,2,3,16,17,18,19,12,13,14,15,8,9,10,11};
748
749 index_t face0R[] = {0,4,7,3,12,19,15,11};
750 index_t face1R[] = {1,2,6,5,9,14,17,13};
751 index_t face2R[] = {0,1,5,4,8,13,16,12};
752 index_t face3R[] = {3,7,6,2,15,18,14,10};
753 index_t face4R[] = {0,3,2,1,11,10,9,8};
754 index_t face5R[] = {4,5,6,7,16,17,18,19};
755
756 /* get MPI information */
757 mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
758 if (! Finley_noError())
759 return NULL;
760
761 /* use the serial version to generate the mesh for the 1-CPU case */
762 if( mpi_info->size==1 )
763 {
764 out = Finley_RectangularMesh_Hex20_singleCPU( numElements, Length, periodic, order, reduced_order, useElementsOnFace, mpi_info);
765 return out;
766 }
767
768 if( mpi_info->rank==0 )
769 domLeft = TRUE;
770 if( mpi_info->rank==mpi_info->size-1 )
771 domRight = TRUE;
772 if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )
773 domInternal = TRUE;
774
775 /* dimensions of the local subdomain */
776 domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal, DOFBoundary );
777
778 NFaceElements=0;
779 if (!periodic[0]) {
780 NDOF0=N0;
781 NFaceElements+=(domRight+domLeft)*NE1*NE2;
782 } else {
783 NDOF0=N0-1;
784 }
785 if (!periodic[1]) {
786 NDOF1=N1;
787 NFaceElements+=2*numElementsLocal*NE2;
788 } else {
789 NDOF1=N1-1;
790 }
791 if (!periodic[2]) {
792 NDOF2=N2;
793 NFaceElements+=2*numElementsLocal*NE1;
794 } else {
795 NDOF2=N2-1;
796 }
797
798 boundaryLeft = !domLeft || periodicLocal[0];
799 boundaryRight = !domRight || periodicLocal[1];
800 N0t = numNodesLocal + boundaryRight + boundaryLeft*2;
801 NDOF0t = numDOFLocal + boundaryRight + boundaryLeft*2;
802 firstNodeConstruct = firstNode - 2*boundaryLeft;
803 firstNodeConstruct = firstNodeConstruct<0 ? N0+firstNodeConstruct-1 : firstNodeConstruct;
804
805 /* allocate mesh: */
806
807 sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
808 out=Finley_Mesh_alloc(name,3,order, reduced_order, mpi_info);
809
810 if (! Finley_noError()) return NULL;
811
812 out->Elements=Finley_ElementFile_alloc(Hex20,out->order,out->reduced_order, mpi_info);
813 if (useElementsOnFace) {
814 out->FaceElements=Finley_ElementFile_alloc(Hex20Face,out->order,out->reduced_order, mpi_info);
815 out->ContactElements=Finley_ElementFile_alloc(Hex20Face_Contact,out->order,out->reduced_order, mpi_info);
816 } else {
817 out->FaceElements=Finley_ElementFile_alloc(Rec8,out->order,out->reduced_order,mpi_info);
818 out->ContactElements=Finley_ElementFile_alloc(Rec8_Contact,out->order,out->reduced_order,mpi_info);
819 }
820 out->Points=Finley_ElementFile_alloc(Point1,out->order,out->reduced_order,mpi_info);
821 if (! Finley_noError()) {
822 Finley_Mesh_dealloc(out);
823 return NULL;
824 }
825
826 /* allocate tables: */
827 Finley_NodeFile_allocTable(out->Nodes,N0t*N1*N2);
828 Finley_ElementFile_allocTable(out->Elements,numElementsLocal*NE1*NE2);
829 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
830
831 Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal*NDOF1*NDOF2, NDOF1*NDOF2*3, 0 );
832 Finley_ElementDistribution_allocTable( out->Elements->elementDistribution, numElementsLocal*NE1*NE2, NE1*NE2*(numElementsLocal-boundaryRight*(!periodic[1])) );
833 Finley_ElementDistribution_allocTable( out->FaceElements->elementDistribution, NFaceElements, NFaceElements-2*boundaryRight*(NE2*(!periodic[1])+NE1*(!periodic[2])) );
834 if (! Finley_noError()) {
835 Finley_Mesh_dealloc(out);
836 return NULL;
837 }
838
839 k=0;
840 #pragma omp parallel for private(i0,i1,i2,k)
841 for (i2=0;i2<N2;i2++) {
842 for (i1=0;i1<N1;i1++) {
843 for (i0=0;i0<N0t;i0++,k++) {
844 out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE((i0+firstNodeConstruct) % N0)/DBLE(N0-1)*Length[0];
845 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
846 out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
847 out->Nodes->Id[k]=k;
848 out->Nodes->Tag[k]=0;
849 out->Nodes->degreeOfFreedom[k]=i0 + (i1%NDOF1)*N0t + (i2%NDOF2)*N0t*N1;
850 out->Nodes->Dom[k]=NODE_INTERNAL;
851 }
852 }
853 }
854
855 /* mark the nodes that reference external and boundary DOF as such */
856 if( boundaryLeft ){
857 for( i1=0; i1<N1; i1++ )
858 for( i2=0; i2<N2; i2++ ) {
859 out->Nodes->Dom[N1*N0t*i2+N0t*i1] = NODE_EXTERNAL;
860 out->Nodes->Dom[N1*N0t*i2+N0t*i1+1] = NODE_EXTERNAL;
861 out->Nodes->Dom[N1*N0t*i2+N0t*i1+2] = NODE_BOUNDARY;
862 }
863 }
864 if( boundaryRight ){
865 for( i1=0; i1<N1; i1++ )
866 for( i2=0; i2<N2; i2++ ) {
867 out->Nodes->Dom[N1*N0t*i2+N0t*(i1+1)-1] = NODE_EXTERNAL;
868 out->Nodes->Dom[N1*N0t*i2+N0t*(i1+1)-2] = NODE_BOUNDARY;
869 out->Nodes->Dom[N1*N0t*i2+N0t*(i1+1)-3] = NODE_BOUNDARY;
870 }
871 }
872 if( periodicLocal[0] ){
873 for( i1=0; i1<N1; i1++ )
874 for( i2=0; i2<N2; i2++ ) {
875 out->Nodes->degreeOfFreedom[N1*N0t*i2+i1*N0t+3] = out->Nodes->degreeOfFreedom[i1*N0t+2];
876 out->Nodes->Dom[N1*N0t*i2+N0t*i1+3] = NODE_BOUNDARY;
877 }
878 }
879
880 /* tag Nodes that are referenced by face elements */
881 if (!periodic[2]) {
882 for (i1=0;i1<N1;i1++) {
883 for (i0=0;i0<N0t;i0++) {
884 out->Nodes->Tag[i0 + N0t*i1]+=100;
885 out->Nodes->Tag[i0 + N0t*i1 + N0t*N1*(N2-1)]+=200;
886 }
887 }
888 }
889 if (!periodic[1]) {
890 for (i2=0;i2<N2;i2++) {
891 for (i0=0;i0<N0t;i0++) {
892 out->Nodes->Tag[i0 + i2*N1*N0t]+=10;
893 out->Nodes->Tag[i0 + (i2+1)*N1*N0t-N0t]+=20;
894 }
895 }
896 }
897 if (!periodic[0] && !domInternal ) {
898 for (i2=0;i2<N2;i2++) {
899 for (i1=0;i1<N1;i1++) {
900 if( domLeft )
901 out->Nodes->Tag[i1*N0t + i2*N0t*N1]+=1;
902 if( domRight )
903 out->Nodes->Tag[(i1+1)*N0t-1 + i2*N0t*N1]+=2;
904 }
905 }
906 }
907
908 /* form the boudary communication information */
909 forwardDOF = MEMALLOC(NDOF1*NDOF2*2,index_t);
910 backwardDOF = MEMALLOC(NDOF1*NDOF2*2,index_t);
911 if( !(mpi_info->size==2 && periodicLocal[0])){
912 if( boundaryLeft ) {
913 targetDomain = mpi_info->rank-1 < 0 ? mpi_info->size-1 : mpi_info->rank-1;
914 for( iI=0, i2=0; i2<NDOF2; i2++ ){
915 for( i1=0; i1<NDOF1; i1++, iI+=2 ){
916 forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+2];
917 backwardDOF[iI] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t];
918 backwardDOF[iI+1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+1];
919 }
920 }
921 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
922 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2*2, backwardDOF );
923 }
924 if( boundaryRight ) {
925 targetDomain = mpi_info->rank+1 > mpi_info->size-1 ? 0 : mpi_info->rank+1;
926 for( iI=0,i2=0; i2<NDOF2; i2++ ){
927 for( i1=0; i1<NDOF1; i1++, iI+=2 ){
928 forwardDOF[iI] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-3];
929 forwardDOF[iI+1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-2];
930 backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-1];
931 }
932 }
933 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2*2, forwardDOF );
934 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
935 }
936 } else{
937 /* periodic boundary conditions with 2 domains, need to change the order in which domain 0 passes boundary data */
938 targetDomain = 1;
939
940 for( iI=0,i2=0; i2<NDOF2; i2++ ){
941 for( i1=0; i1<NDOF1; i1++, iI+=2 ){
942 forwardDOF[iI] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-3];
943 forwardDOF[iI+1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-2];
944 backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-1];
945 }
946 }
947 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2*2, forwardDOF );
948 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
949
950 for( iI=0, i2=0; i2<NDOF2; i2++ ){
951 for( i1=0; i1<NDOF1; i1++, iI+=2 ){
952 forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+2];
953 backwardDOF[iI] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t];
954 backwardDOF[iI+1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+1];
955 }
956 }
957 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
958 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2*2, backwardDOF );
959 }
960 MEMFREE( forwardDOF );
961 MEMFREE( backwardDOF );
962 /* set the elements: */
963
964 k=0;
965 #pragma omp parallel for private(i0,i1,i2,k,node0)
966 for (i2=0;i2<NE2;i2++) {
967 for (i1=0;i1<NE1;i1++) {
968 for (i0=0;i0<numElementsLocal;i0++,k++) {
969 node0 = (periodicLocal[0] && !i0) ? 2*(i1*N0t + i2*N1*N0t) : 2*(i1*N0t + i2*N1*N0t + i0) + periodicLocal[0];
970
971 out->Elements->Id[k]=((firstNodeConstruct/2+i0)%NE0)*NE1*NE2 + NE1*i2 + i1;
972 out->Elements->Tag[k]=0;
973 out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);
974 out->Elements->Dom[k]=ELEMENT_INTERNAL;
975
976 out->Elements->Nodes[INDEX2(0,k,20)]=node0;
977 out->Elements->Nodes[INDEX2(1,k,20)]=node0+2;
978 out->Elements->Nodes[INDEX2(2,k,20)]=node0+2*N0t+2;
979 out->Elements->Nodes[INDEX2(3,k,20)]=node0+2*N0t;
980 out->Elements->Nodes[INDEX2(4,k,20)]=node0+2*N0t*N1;
981 out->Elements->Nodes[INDEX2(5,k,20)]=node0+2*N0t*N1+2;
982 out->Elements->Nodes[INDEX2(6,k,20)]=node0+2*N0t*N1+2*N0t+2;
983 out->Elements->Nodes[INDEX2(7,k,20)]=node0+2*N0t*N1+2*N0t;
984 out->Elements->Nodes[INDEX2(8,k,20)]=node0+1;
985 out->Elements->Nodes[INDEX2(9,k,20)]=node0+N0t+2;
986 out->Elements->Nodes[INDEX2(10,k,20)]=node0+2*N0t+1;
987 out->Elements->Nodes[INDEX2(11,k,20)]=node0+N0t;
988 out->Elements->Nodes[INDEX2(12,k,20)]=node0+N0t*N1;
989 out->Elements->Nodes[INDEX2(13,k,20)]=node0+N0t*N1+2;
990 out->Elements->Nodes[INDEX2(14,k,20)]=node0+N0t*N1+2*N0t+2;
991 out->Elements->Nodes[INDEX2(15,k,20)]=node0+N0t*N1+2*N0t;
992 out->Elements->Nodes[INDEX2(16,k,20)]=node0+2*N0t*N1+1;
993 out->Elements->Nodes[INDEX2(17,k,20)]=node0+2*N0t*N1+N0t+2;
994 out->Elements->Nodes[INDEX2(18,k,20)]=node0+2*N0t*N1+2*N0t+1;
995 out->Elements->Nodes[INDEX2(19,k,20)]=node0+2*N0t*N1+N0t;
996 }
997 }
998 }
999 out->Elements->minColor=0;
1000 out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
1001
1002 if( boundaryLeft )
1003 for( i2=0; i2<NE2; i2++ )
1004 for( i1=0; i1<NE1; i1++ )
1005 out->Elements->Dom[i2*NE1*numElementsLocal+i1*numElementsLocal]=ELEMENT_BOUNDARY;
1006 if( boundaryRight )
1007 for( i2=0; i2<NE2; i2++ )
1008 for( i1=0; i1<NE1; i1++ )
1009 out->Elements->Dom[i2*NE1*numElementsLocal+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
1010
1011 out->Elements->numElements = numElementsLocal*NE1*NE2;
1012 Finley_ElementFile_setDomainFlags( out->Elements );
1013
1014 /* face elements: */
1015
1016 if (useElementsOnFace) {
1017 NUMNODES=20;
1018 } else {
1019 NUMNODES=8;
1020 }
1021 totalNECount=out->Elements->numElements;
1022 faceNECount=0;
1023 idCount = totalNECount;
1024
1025 /* these are the quadrilateral elements on boundary 1 (x3=0): */
1026 numElementsInternal = numElementsLocal-nodesExternal[0]-nodesExternal[1];
1027 if (!periodic[0] && !domInternal) {
1028 /* ** elements on boundary 001 (x1=0): */
1029 if( domLeft ){
1030 #pragma omp parallel for private(i1,i2,k)
1031 for (i2=0;i2<NE2;i2++) {
1032 for (i1=0;i1<NE1;i1++) {
1033 k=i1+NE1*i2+faceNECount;
1034 kk=i1*numElementsLocal + i2*numElementsLocal*NE1;
1035 facePerm =!useElementsOnFace ? face0R : face0;
1036
1037 out->FaceElements->Id[k]=idCount++;
1038 out->FaceElements->Tag[k]=1;
1039 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
1040 out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
1041
1042 for( j=0; j<NUMNODES; j++ )
1043 out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,20)];
1044 }
1045 }
1046 totalNECount+=NE1*NE2;
1047 faceNECount+=NE1*NE2;
1048 }
1049 /* ** elements on boundary 002 (x1=1): */
1050 if( domRight ) {
1051 #pragma omp parallel for private(i1,i2,k)
1052 for (i2=0;i2<NE2;i2++) {
1053 for (i1=0;i1<NE1;i1++) {
1054 k=i1+NE1*i2+faceNECount;
1055 kk=(i1+1)*numElementsLocal + i2*numElementsLocal*NE1 - 1;
1056 facePerm =!useElementsOnFace ? face1R : face1;
1057
1058 out->FaceElements->Id[k]=idCount++;
1059 out->FaceElements->Tag[k]=2;
1060 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
1061 out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
1062
1063 for( j=0; j<NUMNODES; j++ )
1064 out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,20)];
1065 }
1066 }
1067 totalNECount+=NE1*NE2;
1068 faceNECount+=NE1*NE2;
1069 }
1070 }
1071 if (!periodic[1]) {
1072 /* ** elements on boundary 010 (x2=0): */
1073
1074 #pragma omp parallel for private(i0,i2,k)
1075 for (i2=0;i2<NE2;i2++) {
1076 for (i0=0;i0<numElementsLocal;i0++) {
1077 k=i0+numElementsLocal*i2+faceNECount;
1078 kk=i0+numElementsLocal*NE1*i2;
1079 facePerm =!useElementsOnFace ? face2R : face2;
1080
1081 out->FaceElements->Id[k]=idCount++;
1082 out->FaceElements->Tag[k]=10;
1083 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
1084 out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;
1085
1086 for( j=0; j<NUMNODES; j++ )
1087 out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,20)];
1088 }
1089 }
1090 if( boundaryLeft ){
1091 for( i2=0; i2<NE2; i2++ )
1092 out->FaceElements->Dom[faceNECount+i2*numElementsLocal]=ELEMENT_BOUNDARY;
1093 if( periodicLocal[0] )
1094 for( i2=0; i2<NE2; i2++ )
1095 out->FaceElements->Dom[faceNECount+i2*numElementsLocal+1]=ELEMENT_BOUNDARY;
1096 }
1097 if( boundaryRight )
1098 for( i2=0; i2<NE2; i2++ )
1099 out->FaceElements->Dom[faceNECount+(i2+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
1100 totalNECount+=numElementsLocal*NE2;
1101 faceNECount+=numElementsLocal*NE2;
1102
1103 /* ** elements on boundary 020 (x2=1): */
1104
1105 #pragma omp parallel for private(i0,i2,k)
1106 for (i2=0;i2<NE2;i2++) {
1107 for (i0=0;i0<numElementsLocal;i0++) {
1108 k=i0+numElementsLocal*i2+faceNECount;
1109 kk=i0+numElementsLocal*NE1*(i2+1)-numElementsLocal;
1110 facePerm =!useElementsOnFace ? face3R : face3;
1111
1112 out->FaceElements->Tag[k]=20;
1113 out->FaceElements->Id[k]=idCount++;
1114 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
1115 out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
1116
1117 for( j=0; j<NUMNODES; j++ )
1118 out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,20)];
1119 }
1120 }
1121 if( boundaryLeft ){
1122 for( i2=0; i2<NE2; i2++ )
1123 out->FaceElements->Dom[faceNECount+i2*numElementsLocal]=ELEMENT_BOUNDARY;
1124 if( periodicLocal[0] )
1125 for( i2=0; i2<NE2; i2++ )
1126 out->FaceElements->Dom[faceNECount+i2*numElementsLocal+1]=ELEMENT_BOUNDARY;
1127 }
1128 if( boundaryRight )
1129 for( i2=0; i2<NE2; i2++ )
1130 out->FaceElements->Dom[faceNECount+(i2+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
1131 totalNECount+=numElementsLocal*NE2;
1132 faceNECount+=numElementsLocal*NE2;
1133 }
1134 if (!periodic[2]) {
1135 /* elements on boundary 100 (x3=0): */
1136
1137 #pragma omp parallel for private(i0,i1,k)
1138 for (i1=0;i1<NE1;i1++) {
1139 for (i0=0; i0<numElementsLocal; i0++) {
1140 k=i0+numElementsLocal*i1+faceNECount;
1141 kk=i0 + i1*numElementsLocal;
1142 facePerm =!useElementsOnFace ? face4R : face4;
1143
1144 out->FaceElements->Id[k]=idCount++;
1145 out->FaceElements->Tag[k]=100;
1146 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
1147 out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
1148
1149 for( j=0; j<NUMNODES; j++ )
1150 out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,20)];
1151 }
1152 }
1153 if( boundaryLeft ){
1154 for( i1=0; i1<NE1; i1++ )
1155 out->FaceElements->Dom[faceNECount+i1*numElementsLocal]=ELEMENT_BOUNDARY;
1156 if( periodicLocal[0] )
1157 for( i1=0; i1<NE1; i1++ )
1158 out->FaceElements->Dom[faceNECount+i1*numElementsLocal+1]=ELEMENT_BOUNDARY;
1159 }
1160 if( boundaryRight )
1161 for( i1=0; i1<NE1; i1++ )
1162 out->FaceElements->Dom[faceNECount+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
1163 totalNECount+=NE1*numElementsLocal;
1164 faceNECount+=NE1*numElementsLocal;
1165
1166 /* ** elements on boundary 200 (x3=1) */
1167
1168 #pragma omp parallel for private(i0,i1,k)
1169 for (i1=0;i1<NE1;i1++) {
1170 for (i0=0;i0<numElementsLocal;i0++) {
1171 k=i0+numElementsLocal*i1+faceNECount;
1172 kk=i0+i1*numElementsLocal+numElementsLocal*NE1*(NE2-1);
1173 facePerm = !useElementsOnFace ? face5R : face5;
1174
1175 out->FaceElements->Id[k]=idCount++;
1176 out->FaceElements->Tag[k]=200;
1177 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
1178 out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
1179
1180 for( j=0; j<NUMNODES; j++ )
1181 out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,20)];
1182 }
1183 }
1184 if( boundaryLeft ){
1185 for( i1=0; i1<NE1; i1++ )
1186 out->FaceElements->Dom[faceNECount+i1*numElementsLocal]=ELEMENT_BOUNDARY;
1187 if( periodicLocal[0] )
1188 for( i1=0; i1<NE1; i1++ )
1189 out->FaceElements->Dom[faceNECount+i1*numElementsLocal+1]=ELEMENT_BOUNDARY;
1190 }
1191 if( boundaryRight )
1192 for( i1=0; i1<NE1; i1++ )
1193 out->FaceElements->Dom[faceNECount+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
1194 totalNECount+=NE1*numElementsLocal;
1195 faceNECount+=NE1*numElementsLocal;
1196 }
1197 out->FaceElements->elementDistribution->numInternal = faceNECount;
1198
1199 out->FaceElements->minColor=0;
1200 out->FaceElements->maxColor=23;
1201 out->FaceElements->numElements=faceNECount;
1202 Finley_ElementFile_setDomainFlags( out->FaceElements );
1203
1204 /* setup distribution info for other elements */
1205 Finley_ElementFile_setDomainFlags( out->ContactElements );
1206 Finley_ElementFile_setDomainFlags( out->Points );
1207
1208 /* reorder the degrees of freedom */
1209 Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
1210
1211 /* condense the nodes: */
1212 Finley_Mesh_resolveNodeIds(out);
1213 if( !Finley_MPI_noError(mpi_info) )
1214 {
1215 Paso_MPIInfo_dealloc( mpi_info );
1216 Finley_Mesh_dealloc(out);
1217 return NULL;
1218 }
1219
1220 /* prepare mesh for further calculatuions:*/
1221 Finley_Mesh_prepare(out);
1222 if( !Finley_MPI_noError(mpi_info) )
1223 {
1224 Paso_MPIInfo_dealloc( mpi_info );
1225 Finley_Mesh_dealloc(out);
1226 return NULL;
1227 }
1228
1229 /* free up memory */
1230 Paso_MPIInfo_dealloc( mpi_info );
1231
1232 //print_mesh_statistics( out, FALSE );
1233
1234 #ifdef Finley_TRACE
1235 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
1236 #endif
1237 if (Finley_noError()) {
1238 if (!Finley_Mesh_isPrepared(out) ) {
1239 Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");
1240 }
1241 }
1242 return out;
1243 }
1244 #endif
1245

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26