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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 751 - (show annotations)
Mon Jun 26 01:46:34 2006 UTC (13 years, 4 months ago) by bcumming
File MIME type: text/plain
File size: 64654 byte(s)
Changes relating to the MPI version of escript
The standard OpenMP version of escript is unchanged

- updated data types (Finley_Mesh, Finley_NodeFile, etc) to store meshes
  over multiple MPI processes.
- added CommBuffer code in Paso for communication of Data associated
  with distributed meshes
- updates in Finley and Escript to support distributed data and operations
  on distributed data (such as interpolation).
- construction of RHS in MPI, so that simple explicit schemes (such as
  /docs/examples/wave.py without IO and the Locator) can run in MPI.
- updated mesh generation for first order line, rectangle and brick
  meshes and second order line meshes in MPI.        
- small changes to trunk/SConstruct and trunk/scons/ess_options.py to
  build the MPI version, these changes are turned off by default.

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 first order elements (Hex8) 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 /* Author: gross@access.edu.au */
24 /* Version: $Id$ */
25
26 /**************************************************************/
27
28 #include "RectangularMesh.h"
29
30 /**************************************************************/
31
32 #ifdef PASO_MPI
33 /* get the number of nodes/elements for domain with rank=rank, of size processors
34 where n is the total number of nodes/elements in the global domain */
35 static index_t domain_MODdim( index_t rank, index_t size, index_t n )
36 {
37 rank = size-rank-1;
38
39 if( rank < n%size )
40 return (index_t)floor(n/size)+1;
41 return (index_t)floor(n/size);
42 }
43
44
45 /* Determines the number of nodes/elements etc along an axis which is numElementsGlobal long for domain rank */
46 /* A bit messy, but it only has to be done once... */
47 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 )
48 {
49 index_t i0;
50 dim_t numNodesGlobal = numElementsGlobal+1;
51
52 (*numNodesLocal) = domain_MODdim( rank, size, numNodesGlobal );
53
54 numElementsLocal[0] = numNodesLocal[0]+1;
55 periodicLocal[0] = periodicLocal[1] = FALSE;
56 nodesExternal[0] = nodesExternal[1] = 1;
57 if( periodic )
58 {
59 if( size==1 )
60 {
61 numElementsLocal[0] = numElementsGlobal;
62 nodesExternal[0] = nodesExternal[1] = 0;
63 periodicLocal[0] = periodicLocal[1] = TRUE;
64 }
65 else
66 {
67 if( rank==0 )
68 {
69 periodicLocal[0] = TRUE;
70 numNodesLocal[0]++;
71 }
72 if( rank==(size-1) )
73 {
74 periodicLocal[1] = TRUE;
75 numNodesLocal[0]--;
76 numElementsLocal[0]--;
77 }
78 }
79 }
80 else if( !periodic )
81 {
82 if( rank==0 ){
83 nodesExternal[0]--;
84 numElementsLocal[0]--;
85 }
86 if( rank==(size-1) )
87 {
88 nodesExternal[1]--;
89 numElementsLocal[0]--;
90 }
91 }
92 numNodesExternal[0] = nodesExternal[0]+nodesExternal[1];
93
94 numElementsInternal[0] = numElementsLocal[0];
95 if( (rank==0) && (rank==size-1) );
96
97 else if( !periodic && ( (rank==0) ^ (rank==size-1) ) )
98 numElementsInternal[0] -= 1;
99 else
100 numElementsInternal[0] -= 2;
101
102 firstNode[0] = 0;
103 for( i0=0; i0<rank; i0++ )
104 firstNode[0] += domain_MODdim( i0, size, numNodesGlobal );
105
106 numDOFLocal[0] = numNodesLocal[0];
107 if( periodicLocal[0] )
108 {
109 numDOFLocal[0]--;
110 }
111 DOFExternal[0] = nodesExternal[0];
112 DOFExternal[1] = nodesExternal[1];
113 }
114
115 void print_mesh_statistics( Finley_Mesh *out )
116 {
117 index_t i, j, N;
118
119 printf( "\nNodes\n=====\n\n" );
120 printf( "\t%d internal DOF\n\t%d boundary DOF\n\t%d local DOF\n\t%d external DOF\n", out->Nodes->degreeOfFreedomDistribution->numInternal, out->Nodes->degreeOfFreedomDistribution->numBoundary, out->Nodes->degreeOfFreedomDistribution->numLocal, out->Nodes->degreeOfFreedomDistribution->numExternal);
121 for( i=0; i<out->Nodes->numNodes; i++ )
122 printf( "node %d\t: id %d \tDOF %d \t: tag %d \t: coordinates [%3g, %3g, %3g]\n", i, out->Nodes->Id[i], out->Nodes->degreeOfFreedom[i], out->Nodes->Tag[i], out->Nodes->Coordinates[INDEX2(0,i,3)], out->Nodes->Coordinates[INDEX2(1,i,3)], out->Nodes->Coordinates[INDEX2(2,i,3)] );
123
124 printf( "Elements\n========\n\n" );
125 printf( "\t%d internal\n\t%d boundary\n\t%d local\n", out->Elements->elementDistribution->numInternal, out->Elements->elementDistribution->numBoundary, out->Elements->elementDistribution->numLocal );
126 N = out->Elements->ReferenceElement->Type->numNodes;
127 for( i=0; i<out->Elements->numElements; i++ ){
128 printf( "element %d \t: id %d \t: nodes [ %3d, %3d, %3d, %3d, %3d, %3d, %3d, %3d ]", i, out->Elements->Id[i], out->Elements->Nodes[INDEX2(0,i,8)], out->Elements->Nodes[INDEX2(1,i,8)], out->Elements->Nodes[INDEX2(2,i,8)], out->Elements->Nodes[INDEX2(3,i,8)], out->Elements->Nodes[INDEX2(4,i,8)], out->Elements->Nodes[INDEX2(5,i,8)], out->Elements->Nodes[INDEX2(6,i,8)], out->Elements->Nodes[INDEX2(7,i,8)] );
129 printf( " DOF [ %3d", out->Nodes->degreeOfFreedom[out->Nodes->Id[out->Elements->Nodes[INDEX2(0,i,N)]]] );
130 for( j=1; j<N; j++ )
131 printf( ", %3d", out->Nodes->degreeOfFreedom[out->Nodes->Id[out->Elements->Nodes[INDEX2(j,i,N)]]] );
132 printf( " ]\n" );
133 }
134
135 printf( "\nFace Elements\n==============\n\n" );
136 printf( "\t%d internal\n\t%d boundary\n\t%d local\n", out->FaceElements->elementDistribution->numInternal, out->FaceElements->elementDistribution->numBoundary, out->FaceElements->elementDistribution->numLocal );
137 N = out->FaceElements->ReferenceElement->Type->numNodes;
138 for( i=0; i<out->FaceElements->numElements; i++ ){
139 printf( "face element %d \t: id %d \t: nodes [ %3d", i, out->FaceElements->Id[i], out->FaceElements->Nodes[INDEX2(0,i,N)] );
140 for( j=1; j<N; j++ )
141 printf( ", %3d", out->FaceElements->Nodes[INDEX2(j,i,N)] );
142 printf( " ] DOF [ %3d", out->Nodes->degreeOfFreedom[out->Nodes->Id[out->FaceElements->Nodes[INDEX2(0,i,N)]]] );
143 for( j=1; j<N; j++ )
144 printf( ", %3d", out->Nodes->degreeOfFreedom[out->Nodes->Id[out->FaceElements->Nodes[INDEX2(j,i,N)]]] );
145 printf( " ]\n" );
146 }
147 }
148
149 #endif
150
151 #ifndef PASO_MPI
152 Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
153 #else
154 Finley_Mesh* Finley_RectangularMesh_Hex8_singleCPU(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
155 #endif
156 {
157 dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES,M0,M1,M2;
158 index_t node0;
159 Finley_Mesh* out;
160 char name[50];
161 double time0=Finley_timer();
162 NE0=MAX(1,numElements[0]);
163 NE1=MAX(1,numElements[1]);
164 NE2=MAX(1,numElements[2]);
165 N0=NE0+1;
166 N1=NE1+1;
167 N2=NE2+1;
168 #ifdef PASO_MPI
169 Paso_MPIInfo *mpi_info = NULL;
170
171 /* get MPI information */
172 mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
173 if (! Finley_noError())
174 return NULL;
175 #endif
176
177 if (N0<=MIN(N1,N2)) {
178 if (N1 <= N2) {
179 M0=1;
180 M1=N0;
181 M2=N0*N1;
182 } else {
183 M0=1;
184 M2=N0;
185 M1=N0*N2;
186 }
187 } else if (N1<=MIN(N2,N0)) {
188 if (N2 <= N0) {
189 M1=1;
190 M2=N1;
191 M0=N2*N1;
192 } else {
193 M1=1;
194 M0=N1;
195 M2=N1*N0;
196 }
197 } else {
198 if (N0 <= N1) {
199 M2=1;
200 M0=N2;
201 M1=N2*N0;
202 } else {
203 M2=1;
204 M1=N2;
205 M0=N1*N2;
206 }
207 }
208
209
210 NFaceElements=0;
211 if (!periodic[0]) {
212 NDOF0=N0;
213 NFaceElements+=2*NE1*NE2;
214 } else {
215 NDOF0=N0-1;
216 }
217 if (!periodic[1]) {
218 NDOF1=N1;
219 NFaceElements+=2*NE0*NE2;
220 } else {
221 NDOF1=N1-1;
222 }
223 if (!periodic[2]) {
224 NDOF2=N2;
225 NFaceElements+=2*NE0*NE1;
226 } else {
227 NDOF2=N2-1;
228 }
229
230 /* allocate mesh: */
231
232 sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
233
234 #ifndef PASO_MPI
235 out=Finley_Mesh_alloc(name,3,order);
236 #else
237 out=Finley_Mesh_alloc(name,3,order,mpi_info);
238 #endif
239 if (! Finley_noError()) return NULL;
240
241 #ifdef PASO_MPI
242 out->Elements=Finley_ElementFile_alloc(Hex8,out->order,mpi_info);
243 if (useElementsOnFace) {
244 out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,mpi_info);
245 out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order,mpi_info);
246 } else {
247 out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order,mpi_info);
248 out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order,mpi_info);
249 }
250 out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
251 #else
252 out->Elements=Finley_ElementFile_alloc(Hex8,out->order);
253 if (useElementsOnFace) {
254 out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order);
255 out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order);
256 } else {
257 out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order);
258 out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order);
259 }
260 out->Points=Finley_ElementFile_alloc(Point1,out->order);
261 #endif
262 if (! Finley_noError()) {
263 Finley_Mesh_dealloc(out);
264 return NULL;
265 }
266
267
268 /* allocate tables: */
269 Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);
270 #ifdef PASO_MPI
271 Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0*NDOF1*NDOF2, 0, 0 );
272 #endif
273 Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);
274 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
275 if (! Finley_noError()) {
276 Finley_Mesh_dealloc(out);
277 return NULL;
278 }
279
280 /* set nodes: */
281
282 #pragma omp parallel for private(i0,i1,i2,k)
283 for (i2=0;i2<N2;i2++) {
284 for (i1=0;i1<N1;i1++) {
285 for (i0=0;i0<N0;i0++) {
286 k=M0*i0+M1*i1+M2*i2;
287 out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0)/DBLE(N0-1)*Length[0];
288 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
289 out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
290 out->Nodes->Id[k]=i0+N0*i1+N0*N1*i2;
291 out->Nodes->Tag[k]=0;
292 out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1) +M2*(i2%NDOF2);
293 }
294 }
295 }
296 /* tags for the faces: */
297 /* tags for the faces: */
298 if (!periodic[2]) {
299 for (i1=0;i1<N1;i1++) {
300 for (i0=0;i0<N0;i0++) {
301 out->Nodes->Tag[M0*i0+M1*i1+M2*0]+=100;
302 out->Nodes->Tag[M0*i0+M1*i1+M2*(N2-1)]+=200;
303 }
304 }
305 }
306 if (!periodic[1]) {
307 for (i2=0;i2<N2;i2++) {
308 for (i0=0;i0<N0;i0++) {
309 out->Nodes->Tag[M0*i0+M1*0+M2*i2]+=10;
310 out->Nodes->Tag[M0*i0+M1*(N1-1)+M2*i2]+=20;
311 }
312 }
313 }
314 if (!periodic[0]) {
315 for (i2=0;i2<N2;i2++) {
316 for (i1=0;i1<N1;i1++) {
317 out->Nodes->Tag[M0*0+M1*i1+M2*i2]+=1;
318 out->Nodes->Tag[M0*(N0-1)+M1*i1+M2*i2]+=2;
319 }
320 }
321 }
322
323 /* set the elements: */
324
325 #pragma omp parallel for private(i0,i1,i2,k,node0)
326 for (i2=0;i2<NE2;i2++) {
327 for (i1=0;i1<NE1;i1++) {
328 for (i0=0;i0<NE0;i0++) {
329 k=i0+NE0*i1+NE0*NE1*i2;
330 node0=i0+i1*N0+N0*N1*i2;
331
332 out->Elements->Id[k]=k;
333 out->Elements->Tag[k]=0;
334 out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
335
336 out->Elements->Nodes[INDEX2(0,k,8)]=node0;
337 out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
338 out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0+1;
339 out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0;
340 out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0*N1;
341 out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0*N1+1;
342 out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0*N1+N0+1;
343 out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0*N1+N0;
344
345 }
346 }
347 }
348 out->Elements->minColor=0;
349 out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
350
351 /* face elements: */
352
353 if (useElementsOnFace) {
354 NUMNODES=8;
355 } else {
356 NUMNODES=4;
357 }
358 totalNECount=NE0*NE1*NE2;
359 faceNECount=0;
360
361 /* these are the quadrilateral elements on boundary 1 (x3=0): */
362 if (!periodic[2]) {
363 /* ** elements on boundary 100 (x3=0): */
364
365 #pragma omp parallel for private(i0,i1,k,node0)
366 for (i1=0;i1<NE1;i1++) {
367 for (i0=0;i0<NE0;i0++) {
368 k=i0+NE0*i1+faceNECount;
369 node0=i0+i1*N0;
370
371 out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
372 out->FaceElements->Tag[k]=100;
373 out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
374
375 if (useElementsOnFace) {
376 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
377 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
378 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
379 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
380 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0*N1;
381 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+N0;
382 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
383 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1+1;
384 } else {
385 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
386 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
387 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
388 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
389 }
390 }
391 }
392 totalNECount+=NE1*NE0;
393 faceNECount+=NE1*NE0;
394
395 /* ** elements on boundary 200 (x3=1) */
396
397 #pragma omp parallel for private(i0,i1,k,node0)
398 for (i1=0;i1<NE1;i1++) {
399 for (i0=0;i0<NE0;i0++) {
400 k=i0+NE0*i1+faceNECount;
401 node0=i0+i1*N0+N0*N1*(NE2-1);
402
403 out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
404 out->FaceElements->Tag[k]=200;
405 out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
406
407 if (useElementsOnFace) {
408 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
409 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
410 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
411 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
412 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
413 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
414 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+1;
415 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;
416 } else {
417 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
418 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
419 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
420 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
421 }
422
423
424 }
425 }
426 totalNECount+=NE1*NE0;
427 faceNECount+=NE1*NE0;
428 }
429 if (!periodic[0]) {
430 /* ** elements on boundary 001 (x1=0): */
431
432 #pragma omp parallel for private(i1,i2,k,node0)
433 for (i2=0;i2<NE2;i2++) {
434 for (i1=0;i1<NE1;i1++) {
435 k=i1+NE1*i2+faceNECount;
436 node0=i1*N0+N0*N1*i2;
437
438 out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
439 out->FaceElements->Tag[k]=1;
440 out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
441
442 if (useElementsOnFace) {
443 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
444 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
445 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
446 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
447 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
448 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+1;
449 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
450 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0+1;
451 } else {
452 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
453 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
454 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
455 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
456 }
457 }
458 }
459 totalNECount+=NE1*NE2;
460 faceNECount+=NE1*NE2;
461
462 /* ** elements on boundary 002 (x1=1): */
463
464 #pragma omp parallel for private(i1,i2,k,node0)
465 for (i2=0;i2<NE2;i2++) {
466 for (i1=0;i1<NE1;i1++) {
467 k=i1+NE1*i2+faceNECount;
468 node0=(NE0-1)+i1*N0+N0*N1*i2 ;
469
470 out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
471 out->FaceElements->Tag[k]=2;
472 out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
473
474 if (useElementsOnFace) {
475 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
476 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
477 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
478 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
479 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
480 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0;
481 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0;
482 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1;
483 } else {
484 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
485 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
486 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
487 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
488 }
489 }
490 }
491 totalNECount+=NE1*NE2;
492 faceNECount+=NE1*NE2;
493 }
494 if (!periodic[1]) {
495 /* ** elements on boundary 010 (x2=0): */
496
497 #pragma omp parallel for private(i0,i2,k,node0)
498 for (i2=0;i2<NE2;i2++) {
499 for (i0=0;i0<NE0;i0++) {
500 k=i0+NE0*i2+faceNECount;
501 node0=i0+N0*N1*i2;
502
503 out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
504 out->FaceElements->Tag[k]=10;
505 out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;
506
507 if (useElementsOnFace) {
508 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
509 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
510 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
511 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
512 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;
513 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0+1;
514 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1*N0+N0+1;
515 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*N0+N0;
516 } else {
517 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
518 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
519 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
520 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
521 }
522 }
523 }
524 totalNECount+=NE0*NE2;
525 faceNECount+=NE0*NE2;
526
527 /* ** elements on boundary 020 (x2=1): */
528
529 #pragma omp parallel for private(i0,i2,k,node0)
530 for (i2=0;i2<NE2;i2++) {
531 for (i0=0;i0<NE0;i0++) {
532 k=i0+NE0*i2+faceNECount;
533 node0=i0+(NE1-1)*N0+N0*N1*i2;
534
535 out->FaceElements->Tag[k]=20;
536 out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
537 out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
538
539 if (useElementsOnFace) {
540 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
541 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
542 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
543 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
544 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
545 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1;
546 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+1;
547 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
548 } else {
549 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
550 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
551 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
552 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
553 }
554
555 }
556 }
557 totalNECount+=NE0*NE2;
558 faceNECount+=NE0*NE2;
559 }
560 out->FaceElements->minColor=0;
561 out->FaceElements->maxColor=23;
562
563 /* face elements done: */
564
565 #ifdef PASO_MPI
566 /* make sure that the trivial distribution data is correct */
567 out->FaceElements->elementDistribution->numBoundary = 0;
568 out->FaceElements->elementDistribution->numLocal = out->FaceElements->elementDistribution->numInternal = faceNECount;
569 out->Elements->elementDistribution->numBoundary = 0;
570 out->Elements->elementDistribution->numLocal = out->Elements->elementDistribution->numInternal = out->Elements->numElements;
571 out->ContactElements->elementDistribution->numLocal = out->ContactElements->elementDistribution->numInternal = out->ContactElements->elementDistribution->numInternal = 0;
572 out->Points->elementDistribution->numLocal = out->Points->elementDistribution->numInternal = out->Points->elementDistribution->numInternal = 0;
573
574 out->Nodes->degreeOfFreedomDistribution->numInternal = out->Nodes->degreeOfFreedomDistribution->numLocal;
575 out->Nodes->degreeOfFreedomDistribution->numBoundary = 0;
576 #endif
577 /* condense the nodes: */
578 Finley_Mesh_resolveNodeIds(out);
579
580 #ifdef PASO_MPI
581 /* setup the CommBuffer */
582 Finley_NodeDistribution_formCommBuffer( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );
583 if ( !Finley_MPI_noError( mpi_info )) {
584 if( Finley_noError() )
585 Finley_setError( PASO_MPI_ERROR, "Error on another MPI process" );
586 Paso_MPIInfo_dealloc( mpi_info );
587 Finley_Mesh_dealloc(out);
588 return NULL;
589 }
590 #endif
591
592 /* prepare mesh for further calculatuions:*/
593 Finley_Mesh_prepare(out) ;
594
595 #ifdef Finley_TRACE
596 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
597 #endif
598
599 if (! Finley_noError()) {
600 Finley_Mesh_dealloc(out);
601 return NULL;
602 }
603 return out;
604 }
605
606 #ifdef PASO_MPI
607 Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
608 {
609 dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES;//,M0,M1,M2;
610 dim_t idCount, NE0_local, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal;
611 bool_t dom_left, dom_right, dom_internal;
612
613 index_t N0dom;
614 index_t firstNode=0, DOFcount=0, forwardDOF[2], backwardDOF[2], node0, node1, node2;
615 index_t targetDomain=-1;
616 bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE;
617 index_t *indexForward=NULL;
618 Finley_Mesh* out;
619
620 char name[50];
621 Paso_MPIInfo *mpi_info = NULL;
622 double time0=Finley_timer();
623
624 NE0=MAX(1,numElements[0]);
625 NE1=MAX(1,numElements[1]);
626 NE2=MAX(1,numElements[2]);
627 N0=NE0+1;
628 N1=NE1+1;
629 N2=NE2+1;
630
631
632 /* get MPI information */
633 mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
634 if (! Finley_noError())
635 return NULL;
636
637 /* use the serial version to generate the mesh for the 1-CPU case */
638 if( mpi_info->size==1 )
639 {
640 Paso_MPIInfo_dealloc( mpi_info );
641 out = Finley_RectangularMesh_Hex8_singleCPU( numElements, Length, periodic, order, useElementsOnFace );
642 //print_mesh_statistics( out );
643 return out;
644 }
645
646 if( mpi_info->rank==0 )
647 domLeft = TRUE;
648 if( mpi_info->rank==mpi_info->size-1 )
649 domRight = TRUE;
650 if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )
651 domInternal = TRUE;
652
653 /* dimensions of the local subdomain */
654 domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal );
655
656 /* count Degrees of Freedom along each axis */
657 NDOF0 = N0 - periodic[0];
658 NDOF1 = N1 - periodic[1];
659 NDOF2 = N2 - periodic[2];
660
661 /* count face elements */
662 /* internal face elements */
663 NFaceElements = 0;
664 if( !periodic[0] )
665 NFaceElements += (domLeft+domRight)*NE1*NE2;
666 if( !periodic[1] )
667 NFaceElements += 2*numElementsInternal*NE2;
668 if( !periodic[2] )
669 NFaceElements += 2*numElementsInternal*NE1;
670 /* boundary face elements */
671 /* this is looks nasty, but it beats a bunch of nested if/then/else carry-on */
672 NFaceElements += 2*( 2 - (domLeft + domRight)*(!periodic[0]) )*( (!periodic[1])*NE2 + (!periodic[2])*NE1 );
673
674
675 /* allocate mesh: */
676 sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
677
678 out=Finley_Mesh_alloc(name,3,order,mpi_info);
679 if (! Finley_noError()) return NULL;
680
681 out->Elements=Finley_ElementFile_alloc(Hex8,out->order,mpi_info);
682 if (useElementsOnFace) {
683 out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,mpi_info);
684 out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order,mpi_info);
685 } else {
686 out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order,mpi_info);
687 out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order,mpi_info);
688 }
689 out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
690
691 if (! Finley_noError()) {
692 Finley_Mesh_dealloc(out);
693 return NULL;
694 }
695
696
697 /* allocate tables: */
698 Finley_NodeFile_allocTable(out->Nodes,(numNodesLocal+2-!periodic[0]*(domLeft+domRight))*N1*N2);
699 Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal*NDOF1*NDOF2, (DOFExternal[0]+DOFExternal[1])*NDOF1*NDOF2, 0 );
700 Finley_ElementFile_allocTable(out->Elements,(numElementsLocal)*NE1*NE2);
701 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
702 if (! Finley_noError()) {
703 Finley_Mesh_dealloc(out);
704 return NULL;
705 }
706
707 /* set nodes: */
708 /* INTERNAL/BOUNDARY NODES */
709 #pragma omp parallel for private(i0,i1,i2,k)
710 for (k=0,i2=0;i2<N2;i2++) {
711 for (i1=0;i1<N1;i1++) {
712 for (i0=0;i0<numNodesLocal-domLeft*periodic[0];i0++,k++) {
713 out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0+firstNode)/DBLE(N0-1)*Length[0];
714 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
715 out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
716 out->Nodes->Id[k]=k;
717 out->Nodes->Tag[k]=0;
718 out->Nodes->degreeOfFreedom[k]=k;
719 }
720 }
721 }
722 if( domLeft && periodic[0] ) {
723 for (i2=0;i2<N2;i2++) {
724 for (i1=0;i1<N1;i1++, k++) {
725 out->Nodes->Coordinates[INDEX2(0,k,3)]=Length[0];
726 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
727 out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
728 out->Nodes->Id[k]=k;
729 out->Nodes->Tag[k]=0;
730 out->Nodes->degreeOfFreedom[k]=i1*numNodesLocal + i2*numNodesLocal*N1;
731 }
732 }
733 /* tag the faces for this special case */
734 if( !periodic[1] )
735 {
736 for( i2=0; i2<N2; i2++ ){
737 out->Nodes->Tag[k + (i2-N2)*N1 ] += 10;
738 out->Nodes->Tag[k + (i2+1-N2)*N1 -1] += 20;
739 }
740 }
741 if( !periodic[2] )
742 {
743 for( i1=0; i1<N1; i1++ ){
744 out->Nodes->Tag[k -N1*N2 +i1] += 100;
745 out->Nodes->Tag[k -N1 +i1] += 200;
746 }
747 }
748 }
749 /* tags for the faces: */
750 N0dom = (numNodesLocal-periodicLocal[0]);
751 if (!periodic[2]) {
752 for (i1=0;i1<N1;i1++) {
753 for (i0=0;i0<N0dom;i0++) {
754 out->Nodes->Tag[i0 + N0dom*i1]+=100;
755 out->Nodes->Tag[i0 + N0dom*i1 + N0dom*N1*(N2-1)]+=200;
756 }
757 }
758 }
759 if (!periodic[1]) {
760 for (i2=0;i2<N2;i2++) {
761 for (i0=0;i0<N0dom;i0++) {
762 out->Nodes->Tag[i0 + i2*N1*N0dom]+=10;
763 out->Nodes->Tag[i0 + (i2+1)*N1*N0dom-N0dom]+=20;
764 }
765 }
766 }
767 if (!periodic[0] && !domInternal ) {
768 for (i2=0;i2<N2;i2++) {
769 for (i1=0;i1<N1;i1++) {
770 if( domLeft )
771 out->Nodes->Tag[i1*N0dom + i2*N0dom*N1]+=1;
772 if( domRight )
773 out->Nodes->Tag[(i1+1)*N0dom-1 + i2*N0dom*N1]+=2;
774 }
775 }
776 }
777 /* setup the forward communication data for the boundary nodes that we have just defined */
778 /* the case where there are 2 subdomains and periodic[0]=true has to be treated
779 as a special case to because the two domains have two interface boundaries to one-another */
780 indexForward = TMPMEMALLOC( NDOF1*NDOF2, index_t );
781 if( mpi_info->size>2 || !periodic[0] ){
782 if( domInternal || domRight || periodicLocal[0] )
783 {
784 for( int i=0; i<NDOF2; i++ )
785 for( int j=0; j<NDOF1; j++ )
786 indexForward[j+i*NDOF1] = numDOFLocal*j+NDOF1*numDOFLocal*i;
787 targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;
788 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
789 }
790 if( domInternal || domLeft || periodicLocal[1] )
791 {
792 for( int i=0; i<NDOF2; i++ )
793 for( int j=0; j<NDOF1; j++ )
794 indexForward[j+i*NDOF1] = numDOFLocal*(j+1)-1+NDOF1*numDOFLocal*i;
795 targetDomain = (mpi_info->rank+1) % mpi_info->size;
796 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
797 }
798 }
799 else {
800 for( int i=0; i<NDOF2; i++ )
801 for( int j=0; j<NDOF1; j++ )
802 indexForward[j+i*NDOF1] = numDOFLocal*(j+1)-1+NDOF1*numDOFLocal*i;
803 targetDomain = (mpi_info->rank+1) % mpi_info->size;
804 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
805
806 for( int i=0; i<NDOF2; i++ )
807 for( int j=0; j<NDOF1; j++ )
808 indexForward[j+i*NDOF1] = numDOFLocal*j+NDOF1*numDOFLocal*i;
809 targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;
810 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
811 }
812
813 /* EXTERNAL NODES */
814 /* left hand boundary */
815 DOFcount = NDOF1*NDOF2*numDOFLocal;
816 if( (domLeft && periodic[0]) || !domLeft ) {
817 if( (domLeft && periodic[0]) )
818 for (i2=0;i2<N2;i2++) {
819 for (i1=0;i1<N1;i1++, k++) {
820 out->Nodes->Coordinates[INDEX2(0,k,3)]=(1.-DBLE(1)/DBLE(N0-1))*Length[0];
821 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
822 out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
823 out->Nodes->Id[k]=k;
824 out->Nodes->Tag[k]=0;
825 out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;
826 }
827 }
828 else
829 for (i2=0;i2<N2;i2++) {
830 for (i1=0;i1<N1;i1++, k++) {
831 out->Nodes->Coordinates[INDEX2(0,k,3)]=(DBLE(firstNode-1)/DBLE(N0-1))*Length[0];
832 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
833 out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
834 out->Nodes->Id[k]=k;
835 out->Nodes->Tag[k]=0;
836 out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;
837 }
838 }
839 DOFcount += NDOF1*NDOF2;
840 targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;
841 if( !periodic[1] ){
842 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, out->Nodes->degreeOfFreedom + (k-NDOF1*NDOF2) );
843 }
844 else {
845 for( int i=0; i<NDOF2; i++ )
846 for( int j=0; j<NDOF1; j++ )
847 indexForward[j+i*NDOF1] = DOFcount - NDOF1*NDOF2 + j + i*NDOF1;
848 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
849 }
850
851 /* tag the faces for this special case */
852 if( !periodic[1] )
853 {
854 for( i1=0; i1<N1; i1++ ){
855 out->Nodes->Tag[k -N1*N2 +i1] += 10;
856 out->Nodes->Tag[k -N1 +i1] += 20;
857 }
858 }
859 if( periodic[2] )
860 {
861 for( i2=0; i2<N2; i2++ ){
862 out->Nodes->Tag[k +(i2-N2)*N1 ] += 100;
863 out->Nodes->Tag[k +(i2-N2+1)*N1 -1] += 200;
864 }
865 }
866 }
867 if( (domRight && periodic[0]) || !domRight )
868 {
869 if( domRight && periodic[0] )
870 for (i2=0;i2<N2;i2++) {
871 for (i1=0;i1<N1;i1++, k++) {
872 out->Nodes->Coordinates[INDEX2(0,k,3)]=Length[0];
873 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
874 out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
875 out->Nodes->Id[k]=k;
876 out->Nodes->Tag[k]=0;
877 out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;
878 }
879 }
880 else
881 for (i2=0;i2<N2;i2++) {
882 for (i1=0;i1<N1;i1++, k++) {
883 out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(firstNode+numNodesLocal-periodicLocal[0])/DBLE(N0-1)*Length[0];
884 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
885 out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
886 out->Nodes->Id[k]=k;
887 out->Nodes->Tag[k]=0;
888 out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;
889 }
890 }
891 DOFcount += NDOF1*NDOF2;
892
893 targetDomain = mpi_info->rank+1 < mpi_info->size? mpi_info->rank+1 : 0;
894 if( !periodic[1] ){
895 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, out->Nodes->degreeOfFreedom + (k-NDOF1*NDOF2) );
896 }
897 else {
898 for( int i=0; i<NDOF2; i++ )
899 for( int j=0; j<NDOF1; j++ )
900 indexForward[j+i*NDOF1] = DOFcount - NDOF1*NDOF2 + j + i*NDOF1;
901 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
902 }
903
904 /* tag the faces for this special case */
905 if( !periodic[1] )
906 {
907 for( i1=0; i1<N1; i1++ ){
908 out->Nodes->Tag[k -N1*N2 +i1] += 10;
909 out->Nodes->Tag[k -N1 +i1] += 20;
910 }
911 }
912 if( !periodic[2] )
913 {
914 for( i2=0; i2<N2; i2++ ){
915 out->Nodes->Tag[k +(i2-N2)*N1 ] += 100;
916 out->Nodes->Tag[k +(i2-N2+1)*N1 -1] += 200;
917 }
918 }
919 }
920 out->Nodes->degreeOfFreedomDistribution->numInternal = NDOF1*NDOF2*(numDOFLocal - 2 + domRight*(!periodic[0]) + domLeft*(!periodic[0]));
921 out->Nodes->degreeOfFreedomDistribution->numBoundary = out->Nodes->degreeOfFreedomDistribution->numLocal - out->Nodes->degreeOfFreedomDistribution->numInternal;
922
923 TMPMEMFREE( indexForward );
924 /* set the elements: */
925
926 /* INTERNAL elements */
927 N0dom = (numNodesLocal-periodicLocal[0]);
928 k = 0;
929 #pragma omp parallel for private(i0,i1,i2,k,node0)
930 for (i2=0;i2<NE2;i2++) {
931 for (i1=0;i1<NE1;i1++) {
932 for (i0=0;i0<numElementsInternal;i0++,k++) {
933 node0=i0+i1*N0dom+N0dom*N1*i2;
934
935 out->Elements->Id[k]=k;
936
937 out->Elements->Tag[k]=0;
938 out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
939
940 out->Elements->Nodes[INDEX2(0,k,8)]=node0;
941 out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
942 out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0dom+1;
943 out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0dom;
944 out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0dom*N1;
945 out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0dom*N1+1;
946 out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0dom*N1+N0dom+1;
947 out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0dom*N1+N0dom;
948
949 }
950 }
951 }
952 out->Elements->elementDistribution->numInternal = NE1*NE2*numElementsInternal;
953 out->Elements->elementDistribution->numBoundary = 0;
954
955 /* BOUNDARY Elements */
956 /* left boundary */
957 if( !domLeft )
958 {
959 for (i2=0;i2<NE2;i2++) {
960 node0 = numNodesLocal*N1*N2 + i2*N1;
961 for (i1=0;i1<NE1;i1++,node0++,k++) {
962 out->Elements->Id[k]=k;
963 out->Elements->Tag[k]=0;
964 out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);
965
966 out->Elements->Nodes[INDEX2(0,k,8)]=node0;
967 out->Elements->Nodes[INDEX2(1,k,8)]=i1*numNodesLocal + i2*numNodesLocal*N1;
968 out->Elements->Nodes[INDEX2(2,k,8)]=(i1+1)*numNodesLocal + i2*numNodesLocal*N1;
969 out->Elements->Nodes[INDEX2(3,k,8)]=node0+1;
970 out->Elements->Nodes[INDEX2(4,k,8)]=node0+N1;
971 out->Elements->Nodes[INDEX2(5,k,8)]=i1*numNodesLocal + (i2+1)*numNodesLocal*N1;
972 out->Elements->Nodes[INDEX2(6,k,8)]=(i1+1)*numNodesLocal + (i2+1)*numNodesLocal*N1;
973 out->Elements->Nodes[INDEX2(7,k,8)]=node0+N1+1;
974 }
975 }
976 out->Elements->elementDistribution->numBoundary += NE1*NE2;
977 }
978 /* the left periodic boundary is done a little differently to a left internal boundary */
979 else if( (domLeft && periodic[0]) )
980 {
981 for (i2=0;i2<NE2;i2++) {
982 node0 = numDOFLocal*N1*N2 + i2*N1;
983 node1 = node0 + N1*N2;
984 for (i1=0;i1<NE1;i1++,k++,node0++,node1++) {
985
986 out->Elements->Id[k]=k;
987 out->Elements->Tag[k]=0;
988 out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);
989
990 out->Elements->Nodes[INDEX2(0,k,8)]=node1;
991 out->Elements->Nodes[INDEX2(1,k,8)]=node0;
992 out->Elements->Nodes[INDEX2(2,k,8)]=node0+1;
993 out->Elements->Nodes[INDEX2(3,k,8)]=node1+1;
994 out->Elements->Nodes[INDEX2(4,k,8)]=node1+N1;
995 out->Elements->Nodes[INDEX2(5,k,8)]=node0+N1;
996 out->Elements->Nodes[INDEX2(6,k,8)]=node0+N1+1;
997 out->Elements->Nodes[INDEX2(7,k,8)]=node1+N1+1;
998 }
999 }
1000 out->Elements->elementDistribution->numBoundary += NE1*NE2;
1001 }
1002 /* right boundary */
1003 if( !domRight || (domRight && periodic[0]) ){
1004 for (i2=0;i2<NE2;i2++) {
1005 for (i1=0;i1<NE1;i1++,node0++,node1+=numDOFLocal,k++) {
1006 node1 = numDOFLocal -1 + numDOFLocal*i1 + N1*numDOFLocal*i2;
1007 node0 = (numNodesLocal + domInternal + periodicLocal[0])*N1*N2 + i2*N1 + i1;
1008
1009 out->Elements->Id[k]=k;
1010 out->Elements->Tag[k]=0;
1011 out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
1012
1013 out->Elements->Nodes[INDEX2(0,k,8)]=node1;
1014 out->Elements->Nodes[INDEX2(1,k,8)]=node0;
1015 out->Elements->Nodes[INDEX2(2,k,8)]=node0+1;
1016 out->Elements->Nodes[INDEX2(3,k,8)]=node1+N0dom;
1017 out->Elements->Nodes[INDEX2(4,k,8)]=node1+N0dom*N1;
1018 out->Elements->Nodes[INDEX2(5,k,8)]=node0+N1;
1019 out->Elements->Nodes[INDEX2(6,k,8)]=node0+N1+1;
1020 out->Elements->Nodes[INDEX2(7,k,8)]=node1+N0dom*N1+N0dom;
1021 }
1022 }
1023 out->Elements->elementDistribution->numBoundary += NE1*NE2;
1024 }
1025
1026 out->Elements->minColor=0;
1027 out->Elements->maxColor=0;//COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
1028 out->Elements->elementDistribution->numLocal = out->Elements->elementDistribution->numInternal + out->Elements->elementDistribution->numBoundary;
1029
1030 /* face elements: */
1031
1032 if (useElementsOnFace) {
1033 NUMNODES=8;
1034 } else {
1035 NUMNODES=4;
1036 }
1037 totalNECount=k;
1038 faceNECount=0;
1039 idCount = totalNECount;
1040
1041 /* Do internal face elements for each boundary face first */
1042 /* these are the quadrilateral elements on boundary 1 (x3=0): */
1043 numElementsInternal = numElementsLocal-nodesExternal[0]-nodesExternal[1];
1044 if (!periodic[2]) {
1045 /* elements on boundary 100 (x3=0): */
1046
1047 #pragma omp parallel for private(i0,i1,k,node0)
1048 for (i1=0;i1<NE1;i1++) {
1049 for (i0=0; i0<numElementsInternal; i0++) {
1050 k=i0+numElementsInternal*i1+faceNECount;
1051 node0=i0+i1*numDOFLocal;
1052
1053 out->FaceElements->Id[k]=idCount++;
1054 out->FaceElements->Tag[k]=100;
1055 out->FaceElements->Color[k]=0;//(i0%2)+2*(i1%2);
1056
1057 if (useElementsOnFace) {
1058 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1059 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;
1060 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal+1;
1061 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
1062 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal*N1;
1063 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1064 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1065 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1066 } else {
1067 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1068 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;
1069 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal+1;
1070 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
1071 }
1072 }
1073 }
1074 totalNECount+=NE1*numElementsInternal;
1075 faceNECount+=NE1*numElementsInternal;
1076
1077 /* ** elements on boundary 200 (x3=1) */
1078
1079 #pragma omp parallel for private(i0,i1,k,node0)
1080 for (i1=0;i1<NE1;i1++) {
1081 for (i0=0;i0<numElementsInternal;i0++) {
1082 k=i0+numElementsInternal*i1+faceNECount;
1083 node0=i0+i1*numDOFLocal+numDOFLocal*N1*(NE2-1);
1084
1085 out->FaceElements->Id[k]=idCount++;
1086 out->FaceElements->Tag[k]=200;
1087 out->FaceElements->Color[k]=0;//(i0%2)+2*(i1%2)+4;
1088
1089 if (useElementsOnFace) {
1090 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ numDOFLocal * N1;
1091 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ numDOFLocal * N1+1;
1092 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal+1;
1093 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal;
1094 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1095 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
1096 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal+1;
1097 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal;
1098 } else {
1099 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ numDOFLocal * N1;
1100 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ numDOFLocal * N1+1;
1101 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal+1;
1102 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal;
1103 }
1104 }
1105 }
1106 totalNECount+=NE1*numElementsInternal;
1107 faceNECount+=NE1*numElementsInternal;
1108 }
1109 if (!periodic[0] && !domInternal) {
1110 /* ** elements on boundary 001 (x1=0): */
1111 if( domLeft ){
1112 #pragma omp parallel for private(i1,i2,k,node0)
1113 for (i2=0;i2<NE2;i2++) {
1114 for (i1=0;i1<NE1;i1++) {
1115 k=i1+NE1*i2+faceNECount;
1116 node0=i1*numDOFLocal+numDOFLocal*N1*i2;
1117
1118 out->FaceElements->Id[k]=idCount++;
1119 out->FaceElements->Tag[k]=1;
1120 out->FaceElements->Color[k]=0; //(i2%2)+2*(i1%2)+8;
1121
1122 if (useElementsOnFace) {
1123 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1124 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1;
1125 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1126 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal;
1127 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
1128 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1129 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1130 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal+1;
1131 } else {
1132 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1133 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1;
1134 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1135 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal;
1136 }
1137 }
1138 }
1139 totalNECount+=NE1*NE2;
1140 faceNECount+=NE1*NE2;
1141 }
1142 /* ** elements on boundary 002 (x1=1): */
1143 if( domRight ) {
1144 #pragma omp parallel for private(i1,i2,k,node0)
1145 for (i2=0;i2<NE2;i2++) {
1146 for (i1=0;i1<NE1;i1++) {
1147 k=i1+NE1*i2+faceNECount;
1148 node0=(numDOFLocal-2)+i1*numDOFLocal+numDOFLocal*N1*i2 ;
1149
1150 out->FaceElements->Id[k]=idCount++;
1151 out->FaceElements->Tag[k]=2;
1152 out->FaceElements->Color[k]=0;//(i2%2)+2*(i1%2)+12;
1153
1154 if (useElementsOnFace) {
1155 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
1156 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal+1;
1157 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1158 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1159 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1160 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal;
1161 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1162 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal*N1;
1163 } else {
1164 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
1165 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal+1;
1166 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1167 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1168 }
1169 }
1170 }
1171 totalNECount+=NE1*NE2;
1172 faceNECount+=NE1*NE2;
1173 }
1174 }
1175 if (!periodic[1]) {
1176 /* ** elements on boundary 010 (x2=0): */
1177
1178 #pragma omp parallel for private(i0,i2,k,node0)
1179 for (i2=0;i2<NE2;i2++) {
1180 for (i0=0;i0<numElementsInternal;i0++) {
1181 k=i0+numElementsInternal*i2+faceNECount;
1182 node0=i0+numDOFLocal*N1*i2;
1183
1184 out->FaceElements->Id[k]=idCount++;
1185 out->FaceElements->Tag[k]=10;
1186 out->FaceElements->Color[k]=0;//(i0%2)+2*(i2%2)+16;
1187
1188 if (useElementsOnFace) {
1189 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1190 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
1191 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*numDOFLocal+1;
1192 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;
1193 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal;
1194 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal+1;
1195 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal+1;
1196 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;
1197 } else {
1198 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1199 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
1200 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*numDOFLocal+1;
1201 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;
1202 }
1203 }
1204 }
1205 totalNECount+=numElementsInternal*NE2;
1206 faceNECount+=numElementsInternal*NE2;
1207
1208 /* ** elements on boundary 020 (x2=1): */
1209
1210 #pragma omp parallel for private(i0,i2,k,node0)
1211 for (i2=0;i2<NE2;i2++) {
1212 for (i0=0;i0<numElementsInternal;i0++) {
1213 k=i0+numElementsInternal*i2+faceNECount;
1214 node0=i0+(NE1-1)*numDOFLocal+numDOFLocal*N1*i2;
1215
1216 out->FaceElements->Tag[k]=20;
1217 out->FaceElements->Id[k]=idCount++;
1218 out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
1219
1220 if (useElementsOnFace) {
1221 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;
1222 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1223 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1224 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal+1;
1225 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1226 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1;
1227 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1228 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
1229 } else {
1230 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;
1231 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1232 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1233 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal+1;
1234 }
1235
1236 }
1237 }
1238 totalNECount+=numElementsInternal*NE2;
1239 faceNECount+=numElementsInternal*NE2;
1240 }
1241 out->FaceElements->elementDistribution->numInternal = faceNECount;
1242
1243 /* now do the boundary face elements */
1244 /* LHS */
1245 if( !domLeft )
1246 {
1247 if( !periodic[2] ) {
1248
1249 /* x3=0 */
1250 for( i1=0; i1<NE1; i1++ )
1251 {
1252 k = i1+faceNECount;
1253 node0 = i1*numNodesLocal;
1254 node1 = numNodesLocal*N1*N2 + i1;
1255
1256 out->FaceElements->Tag[k]=200;
1257 out->FaceElements->Id[k]=idCount++;
1258 out->FaceElements->Color[k]=0;
1259
1260 if( useElementsOnFace )
1261 {
1262 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1263 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;
1264 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal;
1265 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
1266 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+N1;
1267 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1+1;
1268 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1269 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numNodesLocal*N1;
1270 } else {
1271 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1272 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;
1273 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal;
1274 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
1275 }
1276 }
1277 faceNECount += NE1;
1278 totalNECount += NE1;
1279
1280 /* x3=1 */
1281 for( i1=0; i1<NE1; i1++ )
1282 {
1283 k = i1+faceNECount;
1284 node0 = numNodesLocal*N1*(NE2-1) + i1*numNodesLocal;
1285 node1 = numNodesLocal*N1*N2 + i1 + (NE2-1)*N1;
1286
1287 out->FaceElements->Tag[k]=200;
1288 out->FaceElements->Id[k]=idCount++;
1289 out->FaceElements->Color[k]=0;
1290
1291 if( useElementsOnFace )
1292 {
1293 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;
1294 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal*N1;
1295 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1296 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;
1297 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;
1298 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0;
1299 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal;
1300 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+1;
1301 } else {
1302 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;
1303 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal*N1;
1304 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1305 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;
1306 }
1307 }
1308 faceNECount += NE1;
1309 totalNECount += NE1;
1310 }
1311
1312 if( !periodic[1] ) {
1313 /* x2=0 */
1314 for( i2=0; i2<NE2; i2++ )
1315 {
1316 k = i2+faceNECount;
1317 node0 = i2*numNodesLocal*N1;
1318 node1 = numNodesLocal*N1*N2 + i2*N1;
1319
1320 out->FaceElements->Tag[k]=20;
1321 out->FaceElements->Id[k]=idCount++;
1322 out->FaceElements->Color[k]=0;
1323
1324 if( useElementsOnFace )
1325 {
1326 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1327 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
1328 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1;
1329 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;
1330 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+1;
1331 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numNodesLocal;
1332 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1333 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+N1+1;
1334 } else {
1335 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1336 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
1337 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1;
1338 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;
1339 }
1340 }
1341 faceNECount += NE2;
1342 totalNECount += NE2;
1343
1344 /* x2=1 */
1345 for( i2=0; i2<NE2; i2++ )
1346 {
1347 k = i2+faceNECount;
1348 node0 = i2*numNodesLocal*N1 + numNodesLocal*(NE1-1);
1349 node1 = numNodesLocal*N1*N2 + i2*N1 + (NE1-1);
1350
1351 out->FaceElements->Tag[k]=20;
1352 out->FaceElements->Id[k]=idCount++;
1353 out->FaceElements->Color[k]=0;
1354
1355 if( useElementsOnFace )
1356 {
1357 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;
1358 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0;
1359 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal*N1;
1360 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1;
1361 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;
1362 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numNodesLocal;
1363 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1364 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;
1365 } else {
1366 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;
1367 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numNodesLocal;
1368 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1369 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;
1370 }
1371 }
1372 faceNECount += NE2;
1373 totalNECount += NE2;
1374 }
1375 }
1376
1377 /* RHS */
1378 if( !domRight || periodicLocal[1] )
1379 {
1380 /* the case of left hand boundary domain and periodic boundary condition on its left hand boundary */
1381 if( domLeft && periodic[0] ){
1382 if( !periodic[2] ) {
1383
1384 /* x3=0 */
1385 for( i1=0; i1<NE1; i1++ )
1386 {
1387 k = i1+faceNECount;
1388 node0 = numDOFLocal*N1*N2 + i1;
1389 node1 = numNodesLocal*N1*N2 + i1;
1390
1391 out->FaceElements->Tag[k]=200;
1392 out->FaceElements->Id[k]=idCount++;
1393 out->FaceElements->Color[k]=0;
1394
1395 if( useElementsOnFace )
1396 {
1397 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1398 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;
1399 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
1400 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
1401 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+N1;
1402 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1+1;
1403 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1+1;
1404 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1;
1405 } else {
1406 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1407 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;
1408 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
1409 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
1410 }
1411 }
1412 faceNECount += NE1;
1413 totalNECount += NE1;
1414
1415 /* x3=1 */
1416 for( i1=0; i1<NE1; i1++ )
1417 {
1418 k = i1+faceNECount;
1419 node0 = numDOFLocal*N1*N2 + i1 + (NE2-1)*N1;
1420 node1 = numNodesLocal*N1*N2 + i1 + (NE2-1)*N1;
1421
1422 out->FaceElements->Tag[k]=200;
1423 out->FaceElements->Id[k]=idCount++;
1424 out->FaceElements->Color[k]=0;
1425
1426 if( useElementsOnFace )
1427 {
1428 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;
1429 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+1;
1430 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+1;
1431 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0;
1432 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;
1433 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;
1434 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;
1435 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1;
1436 } else {
1437 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;
1438 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;
1439 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;
1440 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1;
1441 }
1442 }
1443 faceNECount += NE1;
1444 totalNECount += NE1;
1445 }
1446
1447 if( !periodic[1] ) {
1448 /* x2=0 */
1449 for( i2=0; i2<NE2; i2++ )
1450 {
1451 k = i2+faceNECount;
1452 node0 = numDOFLocal*N1*N2 + i2*N1;
1453 node1 = numNodesLocal*N1*N2 + i2*N1;
1454
1455 out->FaceElements->Tag[k]=20;
1456 out->FaceElements->Id[k]=idCount++;
1457 out->FaceElements->Color[k]=0;
1458
1459 if( useElementsOnFace )
1460 {
1461 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1462 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
1463 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1;
1464 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;
1465 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+1;
1466 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
1467 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1+1;
1468 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+N1+1;
1469 } else {
1470 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1471 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
1472 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1;
1473 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;
1474 }
1475 }
1476 faceNECount += NE2;
1477 totalNECount += NE2;
1478
1479 /* x2=1 */
1480 for( i2=0; i2<NE2; i2++ )
1481 {
1482 k = i2+faceNECount;
1483 node0 = numDOFLocal*N1*N2 + i2*N1 + N1-2;
1484 node1 = numNodesLocal*N1*N2 + i2*N1 + N1-2;
1485
1486 out->FaceElements->Tag[k]=20;
1487 out->FaceElements->Id[k]=idCount++;
1488 out->FaceElements->Color[k]=0;
1489
1490 if( useElementsOnFace )
1491 {
1492 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;
1493 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0;
1494 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1;
1495 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1;
1496 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;
1497 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
1498 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;
1499 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;
1500 } else {
1501 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;
1502 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
1503 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;
1504 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;
1505 }
1506 }
1507 faceNECount += NE2;
1508 totalNECount += NE2;
1509 }
1510
1511 }
1512 if( !periodic[2] ) {
1513 /* x3=0 */
1514 for( i1=0; i1<NE1; i1++ )
1515 {
1516 k = i1+faceNECount;
1517 node0 = numDOFLocal*(i1+1) - 1;
1518 node1 = (numNodesLocal+periodicLocal[0]+periodicLocal[1])*N1*N2 + i1 + domInternal*N1*N2;
1519
1520 out->FaceElements->Tag[k]=200;
1521 out->FaceElements->Id[k]=idCount++;
1522 out->FaceElements->Color[k]=0;
1523
1524 if( useElementsOnFace )
1525 {
1526 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1527 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;
1528 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+1;
1529 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1;
1530 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal*N1;
1531 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1532 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+N1+1;
1533 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+N1;
1534 } else {
1535 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1536 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;
1537 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+1;
1538 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1;
1539 }
1540 }
1541 faceNECount += NE1;
1542 totalNECount += NE1;
1543
1544 /* x3=1 */
1545 for( i1=0; i1<NE1; i1++ )
1546 {
1547 k = i1+faceNECount;
1548 node0 = numDOFLocal*N1*(NE2-1) + (i1+1)*numDOFLocal - 1;
1549 node1 = numNodesLocal*N1*N2 + i1 + (NE2-1)*N1 + (domInternal+periodicLocal[1]+periodicLocal[0])*N1*N2;
1550
1551 out->FaceElements->Tag[k]=200;
1552 out->FaceElements->Id[k]=idCount++;
1553 out->FaceElements->Color[k]=0;
1554
1555 if( useElementsOnFace )
1556 {
1557 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1558 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal;
1559 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+1;
1560 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1;
1561 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal*N1;
1562 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1563 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;
1564 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1;
1565 } else {
1566 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal*N1;
1567 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1568 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;
1569 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1;
1570 }
1571 }
1572 faceNECount += NE1;
1573 totalNECount += NE1;
1574 }
1575 if( !periodic[1] ) {
1576 /* x2=0 */
1577 for( i2=0; i2<NE2; i2++ )
1578 {
1579 k = i2+faceNECount;
1580 node0 = N1*numDOFLocal*i2 + numDOFLocal - 1;
1581 node1 = numNodesLocal*N1*N2 + i2*N1 + (domInternal+periodicLocal[0]+periodicLocal[1])*N1*N2;
1582
1583 out->FaceElements->Tag[k]=20;
1584 out->FaceElements->Id[k]=idCount++;
1585 out->FaceElements->Color[k]=0;
1586
1587 if( useElementsOnFace )
1588 {
1589 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1590 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1;
1591 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1;
1592 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;
1593 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal;
1594 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+1;
1595 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+N1+1;
1596 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;
1597 } else {
1598 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1599 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1;
1600 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1;
1601 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;
1602 }
1603 }
1604 faceNECount += NE2;
1605 totalNECount += NE2;
1606
1607 /* x2=1 */
1608 for( i2=0; i2<NE2; i2++ )
1609 {
1610 k = i2+faceNECount;
1611 node0 = numDOFLocal*N1*i2 + NE1*numDOFLocal - 1;
1612 node1 = numNodesLocal*N1*N2 + i2*N1 + (NE1-1) + (domInternal+periodicLocal[0]+periodicLocal[1])*N1*N2;
1613
1614 out->FaceElements->Tag[k]=20;
1615 out->FaceElements->Id[k]=idCount++;
1616 out->FaceElements->Color[k]=0;
1617
1618 if( useElementsOnFace ){
1619 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1620 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1;
1621 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+N1;
1622 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N1*numDOFLocal;
1623 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;
1624 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+1;
1625 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;
1626 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;
1627 } else {
1628 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;
1629 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+1;
1630 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;
1631 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;
1632 }
1633 }
1634 faceNECount += NE2;
1635 totalNECount += NE2;
1636 }
1637 }
1638 out->FaceElements->minColor=0;
1639 out->FaceElements->maxColor=0;//23;
1640
1641 out->FaceElements->elementDistribution->numBoundary = faceNECount - out->FaceElements->elementDistribution->numInternal;
1642 out->FaceElements->elementDistribution->numLocal = faceNECount;
1643
1644
1645 /* setup distribution info for other elements */
1646 out->ContactElements->elementDistribution->numLocal = out->ContactElements->elementDistribution->numInternal = out->ContactElements->elementDistribution->numInternal = 0;
1647 out->Points->elementDistribution->numLocal = out->Points->elementDistribution->numInternal = out->Points->elementDistribution->numInternal = 0;
1648
1649 /* condense the nodes: */
1650 Finley_Mesh_resolveNodeIds( out );
1651
1652 /* setup the CommBuffer */
1653 Finley_NodeDistribution_formCommBuffer( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );
1654 if ( !Finley_MPI_noError( mpi_info )) {
1655 if( Finley_noError() )
1656 Finley_setError( PASO_MPI_ERROR, "Error on another MPI process" );
1657 Paso_MPIInfo_dealloc( mpi_info );
1658 Finley_Mesh_dealloc(out);
1659 return NULL;
1660 }
1661
1662 Finley_NodeDistribution_calculateIndexExternal( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );
1663
1664 /* prepare mesh for further calculatuions:*/
1665 Finley_Mesh_prepare(out) ;
1666
1667 // print_mesh_statistics( out );
1668
1669 #ifdef Finley_TRACE
1670 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
1671 #endif
1672
1673 if (! Finley_noError()) {
1674 Finley_Mesh_dealloc(out);
1675 return NULL;
1676 }
1677
1678 return out;
1679 }
1680 #endif
1681

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26