/[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 759 - (show annotations)
Thu Jun 29 01:53:23 2006 UTC (16 years, 7 months ago) by bcumming
File MIME type: text/plain
File size: 64656 byte(s) 
- added directory pythonMPI to the source tree. this directory contains
  the c++ wrapper that is used to run python scripts in parallel for the
  MPI version of escript/finley
- updated the SConstruct and ./scons/ess_options.py for conditional MPI
  compilation. To compile the MPI version on ESS uncomment the #define
  PASO_MPI in ./paso/src/Paso.h and add the command line option
  useMPI=yes when running scons.
- fixed a compile time error in the MPI build in  
  finley/src/CPPAdapter/MeshAdapterFactory.cpp
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 k=0;
710 #pragma omp parallel for private(i0,i1,i2,k)
711 for (i2=0;i2<N2;i2++) {
712 for (i1=0;i1<N1;i1++) {
713 for (i0=0;i0<numNodesLocal-domLeft*periodic[0];i0++,k++) {
714 out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0+firstNode)/DBLE(N0-1)*Length[0];
715 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
716 out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
717 out->Nodes->Id[k]=k;
718 out->Nodes->Tag[k]=0;
719 out->Nodes->degreeOfFreedom[k]=k;
720 }
721 }
722 }
723 if( domLeft && periodic[0] ) {
724 for (i2=0;i2<N2;i2++) {
725 for (i1=0;i1<N1;i1++, k++) {
726 out->Nodes->Coordinates[INDEX2(0,k,3)]=Length[0];
727 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
728 out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
729 out->Nodes->Id[k]=k;
730 out->Nodes->Tag[k]=0;
731 out->Nodes->degreeOfFreedom[k]=i1*numNodesLocal + i2*numNodesLocal*N1;
732 }
733 }
734 /* tag the faces for this special case */
735 if( !periodic[1] )
736 {
737 for( i2=0; i2<N2; i2++ ){
738 out->Nodes->Tag[k + (i2-N2)*N1 ] += 10;
739 out->Nodes->Tag[k + (i2+1-N2)*N1 -1] += 20;
740 }
741 }
742 if( !periodic[2] )
743 {
744 for( i1=0; i1<N1; i1++ ){
745 out->Nodes->Tag[k -N1*N2 +i1] += 100;
746 out->Nodes->Tag[k -N1 +i1] += 200;
747 }
748 }
749 }
750 /* tags for the faces: */
751 N0dom = (numNodesLocal-periodicLocal[0]);
752 if (!periodic[2]) {
753 for (i1=0;i1<N1;i1++) {
754 for (i0=0;i0<N0dom;i0++) {
755 out->Nodes->Tag[i0 + N0dom*i1]+=100;
756 out->Nodes->Tag[i0 + N0dom*i1 + N0dom*N1*(N2-1)]+=200;
757 }
758 }
759 }
760 if (!periodic[1]) {
761 for (i2=0;i2<N2;i2++) {
762 for (i0=0;i0<N0dom;i0++) {
763 out->Nodes->Tag[i0 + i2*N1*N0dom]+=10;
764 out->Nodes->Tag[i0 + (i2+1)*N1*N0dom-N0dom]+=20;
765 }
766 }
767 }
768 if (!periodic[0] && !domInternal ) {
769 for (i2=0;i2<N2;i2++) {
770 for (i1=0;i1<N1;i1++) {
771 if( domLeft )
772 out->Nodes->Tag[i1*N0dom + i2*N0dom*N1]+=1;
773 if( domRight )
774 out->Nodes->Tag[(i1+1)*N0dom-1 + i2*N0dom*N1]+=2;
775 }
776 }
777 }
778 /* setup the forward communication data for the boundary nodes that we have just defined */
779 /* the case where there are 2 subdomains and periodic[0]=true has to be treated
780 as a special case to because the two domains have two interface boundaries to one-another */
781 indexForward = TMPMEMALLOC( NDOF1*NDOF2, index_t );
782 if( mpi_info->size>2 || !periodic[0] ){
783 if( domInternal || domRight || periodicLocal[0] )
784 {
785 for( int i=0; i<NDOF2; i++ )
786 for( int j=0; j<NDOF1; j++ )
787 indexForward[j+i*NDOF1] = numDOFLocal*j+NDOF1*numDOFLocal*i;
788 targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;
789 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
790 }
791 if( domInternal || domLeft || periodicLocal[1] )
792 {
793 for( int i=0; i<NDOF2; i++ )
794 for( int j=0; j<NDOF1; j++ )
795 indexForward[j+i*NDOF1] = numDOFLocal*(j+1)-1+NDOF1*numDOFLocal*i;
796 targetDomain = (mpi_info->rank+1) % mpi_info->size;
797 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
798 }
799 }
800 else {
801 for( int i=0; i<NDOF2; i++ )
802 for( int j=0; j<NDOF1; j++ )
803 indexForward[j+i*NDOF1] = numDOFLocal*(j+1)-1+NDOF1*numDOFLocal*i;
804 targetDomain = (mpi_info->rank+1) % mpi_info->size;
805 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
806
807 for( int i=0; i<NDOF2; i++ )
808 for( int j=0; j<NDOF1; j++ )
809 indexForward[j+i*NDOF1] = numDOFLocal*j+NDOF1*numDOFLocal*i;
810 targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;
811 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
812 }
813
814 /* EXTERNAL NODES */
815 /* left hand boundary */
816 DOFcount = NDOF1*NDOF2*numDOFLocal;
817 if( (domLeft && periodic[0]) || !domLeft ) {
818 if( (domLeft && periodic[0]) )
819 for (i2=0;i2<N2;i2++) {
820 for (i1=0;i1<N1;i1++, k++) {
821 out->Nodes->Coordinates[INDEX2(0,k,3)]=(1.-DBLE(1)/DBLE(N0-1))*Length[0];
822 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
823 out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
824 out->Nodes->Id[k]=k;
825 out->Nodes->Tag[k]=0;
826 out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;
827 }
828 }
829 else
830 for (i2=0;i2<N2;i2++) {
831 for (i1=0;i1<N1;i1++, k++) {
832 out->Nodes->Coordinates[INDEX2(0,k,3)]=(DBLE(firstNode-1)/DBLE(N0-1))*Length[0];
833 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
834 out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
835 out->Nodes->Id[k]=k;
836 out->Nodes->Tag[k]=0;
837 out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;
838 }
839 }
840 DOFcount += NDOF1*NDOF2;
841 targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;
842 if( !periodic[1] ){
843 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, out->Nodes->degreeOfFreedom + (k-NDOF1*NDOF2) );
844 }
845 else {
846 for( int i=0; i<NDOF2; i++ )
847 for( int j=0; j<NDOF1; j++ )
848 indexForward[j+i*NDOF1] = DOFcount - NDOF1*NDOF2 + j + i*NDOF1;
849 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
850 }
851
852 /* tag the faces for this special case */
853 if( !periodic[1] )
854 {
855 for( i1=0; i1<N1; i1++ ){
856 out->Nodes->Tag[k -N1*N2 +i1] += 10;
857 out->Nodes->Tag[k -N1 +i1] += 20;
858 }
859 }
860 if( periodic[2] )
861 {
862 for( i2=0; i2<N2; i2++ ){
863 out->Nodes->Tag[k +(i2-N2)*N1 ] += 100;
864 out->Nodes->Tag[k +(i2-N2+1)*N1 -1] += 200;
865 }
866 }
867 }
868 if( (domRight && periodic[0]) || !domRight )
869 {
870 if( domRight && periodic[0] )
871 for (i2=0;i2<N2;i2++) {
872 for (i1=0;i1<N1;i1++, k++) {
873 out->Nodes->Coordinates[INDEX2(0,k,3)]=Length[0];
874 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
875 out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
876 out->Nodes->Id[k]=k;
877 out->Nodes->Tag[k]=0;
878 out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;
879 }
880 }
881 else
882 for (i2=0;i2<N2;i2++) {
883 for (i1=0;i1<N1;i1++, k++) {
884 out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(firstNode+numNodesLocal-periodicLocal[0])/DBLE(N0-1)*Length[0];
885 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
886 out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
887 out->Nodes->Id[k]=k;
888 out->Nodes->Tag[k]=0;
889 out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;
890 }
891 }
892 DOFcount += NDOF1*NDOF2;
893
894 targetDomain = mpi_info->rank+1 < mpi_info->size? mpi_info->rank+1 : 0;
895 if( !periodic[1] ){
896 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, out->Nodes->degreeOfFreedom + (k-NDOF1*NDOF2) );
897 }
898 else {
899 for( int i=0; i<NDOF2; i++ )
900 for( int j=0; j<NDOF1; j++ )
901 indexForward[j+i*NDOF1] = DOFcount - NDOF1*NDOF2 + j + i*NDOF1;
902 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
903 }
904
905 /* tag the faces for this special case */
906 if( !periodic[1] )
907 {
908 for( i1=0; i1<N1; i1++ ){
909 out->Nodes->Tag[k -N1*N2 +i1] += 10;
910 out->Nodes->Tag[k -N1 +i1] += 20;
911 }
912 }
913 if( !periodic[2] )
914 {
915 for( i2=0; i2<N2; i2++ ){
916 out->Nodes->Tag[k +(i2-N2)*N1 ] += 100;
917 out->Nodes->Tag[k +(i2-N2+1)*N1 -1] += 200;
918 }
919 }
920 }
921 out->Nodes->degreeOfFreedomDistribution->numInternal = NDOF1*NDOF2*(numDOFLocal - 2 + domRight*(!periodic[0]) + domLeft*(!periodic[0]));
922 out->Nodes->degreeOfFreedomDistribution->numBoundary = out->Nodes->degreeOfFreedomDistribution->numLocal - out->Nodes->degreeOfFreedomDistribution->numInternal;
923
924 TMPMEMFREE( indexForward );
925 /* set the elements: */
926
927 /* INTERNAL elements */
928 N0dom = (numNodesLocal-periodicLocal[0]);
929 k = 0;
930 #pragma omp parallel for private(i0,i1,i2,k,node0)
931 for (i2=0;i2<NE2;i2++) {
932 for (i1=0;i1<NE1;i1++) {
933 for (i0=0;i0<numElementsInternal;i0++,k++) {
934 node0=i0+i1*N0dom+N0dom*N1*i2;
935
936 out->Elements->Id[k]=k;
937
938 out->Elements->Tag[k]=0;
939 out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
940
941 out->Elements->Nodes[INDEX2(0,k,8)]=node0;
942 out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
943 out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0dom+1;
944 out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0dom;
945 out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0dom*N1;
946 out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0dom*N1+1;
947 out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0dom*N1+N0dom+1;
948 out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0dom*N1+N0dom;
949
950 }
951 }
952 }
953 out->Elements->elementDistribution->numInternal = NE1*NE2*numElementsInternal;
954 out->Elements->elementDistribution->numBoundary = 0;
955
956 /* BOUNDARY Elements */
957 /* left boundary */
958 if( !domLeft )
959 {
960 for (i2=0;i2<NE2;i2++) {
961 node0 = numNodesLocal*N1*N2 + i2*N1;
962 for (i1=0;i1<NE1;i1++,node0++,k++) {
963 out->Elements->Id[k]=k;
964 out->Elements->Tag[k]=0;
965 out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);
966
967 out->Elements->Nodes[INDEX2(0,k,8)]=node0;
968 out->Elements->Nodes[INDEX2(1,k,8)]=i1*numNodesLocal + i2*numNodesLocal*N1;
969 out->Elements->Nodes[INDEX2(2,k,8)]=(i1+1)*numNodesLocal + i2*numNodesLocal*N1;
970 out->Elements->Nodes[INDEX2(3,k,8)]=node0+1;
971 out->Elements->Nodes[INDEX2(4,k,8)]=node0+N1;
972 out->Elements->Nodes[INDEX2(5,k,8)]=i1*numNodesLocal + (i2+1)*numNodesLocal*N1;
973 out->Elements->Nodes[INDEX2(6,k,8)]=(i1+1)*numNodesLocal + (i2+1)*numNodesLocal*N1;
974 out->Elements->Nodes[INDEX2(7,k,8)]=node0+N1+1;
975 }
976 }
977 out->Elements->elementDistribution->numBoundary += NE1*NE2;
978 }
979 /* the left periodic boundary is done a little differently to a left internal boundary */
980 else if( (domLeft && periodic[0]) )
981 {
982 for (i2=0;i2<NE2;i2++) {
983 node0 = numDOFLocal*N1*N2 + i2*N1;
984 node1 = node0 + N1*N2;
985 for (i1=0;i1<NE1;i1++,k++,node0++,node1++) {
986
987 out->Elements->Id[k]=k;
988 out->Elements->Tag[k]=0;
989 out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);
990
991 out->Elements->Nodes[INDEX2(0,k,8)]=node1;
992 out->Elements->Nodes[INDEX2(1,k,8)]=node0;
993 out->Elements->Nodes[INDEX2(2,k,8)]=node0+1;
994 out->Elements->Nodes[INDEX2(3,k,8)]=node1+1;
995 out->Elements->Nodes[INDEX2(4,k,8)]=node1+N1;
996 out->Elements->Nodes[INDEX2(5,k,8)]=node0+N1;
997 out->Elements->Nodes[INDEX2(6,k,8)]=node0+N1+1;
998 out->Elements->Nodes[INDEX2(7,k,8)]=node1+N1+1;
999 }
1000 }
1001 out->Elements->elementDistribution->numBoundary += NE1*NE2;
1002 }
1003 /* right boundary */
1004 if( !domRight || (domRight && periodic[0]) ){
1005 for (i2=0;i2<NE2;i2++) {
1006 for (i1=0;i1<NE1;i1++,node0++,node1+=numDOFLocal,k++) {
1007 node1 = numDOFLocal -1 + numDOFLocal*i1 + N1*numDOFLocal*i2;
1008 node0 = (numNodesLocal + domInternal + periodicLocal[0])*N1*N2 + i2*N1 + i1;
1009
1010 out->Elements->Id[k]=k;
1011 out->Elements->Tag[k]=0;
1012 out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
1013
1014 out->Elements->Nodes[INDEX2(0,k,8)]=node1;
1015 out->Elements->Nodes[INDEX2(1,k,8)]=node0;
1016 out->Elements->Nodes[INDEX2(2,k,8)]=node0+1;
1017 out->Elements->Nodes[INDEX2(3,k,8)]=node1+N0dom;
1018 out->Elements->Nodes[INDEX2(4,k,8)]=node1+N0dom*N1;
1019 out->Elements->Nodes[INDEX2(5,k,8)]=node0+N1;
1020 out->Elements->Nodes[INDEX2(6,k,8)]=node0+N1+1;
1021 out->Elements->Nodes[INDEX2(7,k,8)]=node1+N0dom*N1+N0dom;
1022 }
1023 }
1024 out->Elements->elementDistribution->numBoundary += NE1*NE2;
1025 }
1026
1027 out->Elements->minColor=0;
1028 out->Elements->maxColor=0;//COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
1029 out->Elements->elementDistribution->numLocal = out->Elements->elementDistribution->numInternal + out->Elements->elementDistribution->numBoundary;
1030
1031 /* face elements: */
1032
1033 if (useElementsOnFace) {
1034 NUMNODES=8;
1035 } else {
1036 NUMNODES=4;
1037 }
1038 totalNECount=k;
1039 faceNECount=0;
1040 idCount = totalNECount;
1041
1042 /* Do internal face elements for each boundary face first */
1043 /* these are the quadrilateral elements on boundary 1 (x3=0): */
1044 numElementsInternal = numElementsLocal-nodesExternal[0]-nodesExternal[1];
1045 if (!periodic[2]) {
1046 /* elements on boundary 100 (x3=0): */
1047
1048 #pragma omp parallel for private(i0,i1,k,node0)
1049 for (i1=0;i1<NE1;i1++) {
1050 for (i0=0; i0<numElementsInternal; i0++) {
1051 k=i0+numElementsInternal*i1+faceNECount;
1052 node0=i0+i1*numDOFLocal;
1053
1054 out->FaceElements->Id[k]=idCount++;
1055 out->FaceElements->Tag[k]=100;
1056 out->FaceElements->Color[k]=0;//(i0%2)+2*(i1%2);
1057
1058 if (useElementsOnFace) {
1059 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1060 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;
1061 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal+1;
1062 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
1063 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal*N1;
1064 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1065 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1066 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1067 } else {
1068 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1069 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;
1070 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal+1;
1071 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
1072 }
1073 }
1074 }
1075 totalNECount+=NE1*numElementsInternal;
1076 faceNECount+=NE1*numElementsInternal;
1077
1078 /* ** elements on boundary 200 (x3=1) */
1079
1080 #pragma omp parallel for private(i0,i1,k,node0)
1081 for (i1=0;i1<NE1;i1++) {
1082 for (i0=0;i0<numElementsInternal;i0++) {
1083 k=i0+numElementsInternal*i1+faceNECount;
1084 node0=i0+i1*numDOFLocal+numDOFLocal*N1*(NE2-1);
1085
1086 out->FaceElements->Id[k]=idCount++;
1087 out->FaceElements->Tag[k]=200;
1088 out->FaceElements->Color[k]=0;//(i0%2)+2*(i1%2)+4;
1089
1090 if (useElementsOnFace) {
1091 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ numDOFLocal * N1;
1092 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ numDOFLocal * N1+1;
1093 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal+1;
1094 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal;
1095 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1096 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
1097 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal+1;
1098 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal;
1099 } else {
1100 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ numDOFLocal * N1;
1101 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ numDOFLocal * N1+1;
1102 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal+1;
1103 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal;
1104 }
1105 }
1106 }
1107 totalNECount+=NE1*numElementsInternal;
1108 faceNECount+=NE1*numElementsInternal;
1109 }
1110 if (!periodic[0] && !domInternal) {
1111 /* ** elements on boundary 001 (x1=0): */
1112 if( domLeft ){
1113 #pragma omp parallel for private(i1,i2,k,node0)
1114 for (i2=0;i2<NE2;i2++) {
1115 for (i1=0;i1<NE1;i1++) {
1116 k=i1+NE1*i2+faceNECount;
1117 node0=i1*numDOFLocal+numDOFLocal*N1*i2;
1118
1119 out->FaceElements->Id[k]=idCount++;
1120 out->FaceElements->Tag[k]=1;
1121 out->FaceElements->Color[k]=0; //(i2%2)+2*(i1%2)+8;
1122
1123 if (useElementsOnFace) {
1124 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1125 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1;
1126 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1127 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal;
1128 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
1129 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1130 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1131 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal+1;
1132 } else {
1133 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1134 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1;
1135 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1136 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal;
1137 }
1138 }
1139 }
1140 totalNECount+=NE1*NE2;
1141 faceNECount+=NE1*NE2;
1142 }
1143 /* ** elements on boundary 002 (x1=1): */
1144 if( domRight ) {
1145 #pragma omp parallel for private(i1,i2,k,node0)
1146 for (i2=0;i2<NE2;i2++) {
1147 for (i1=0;i1<NE1;i1++) {
1148 k=i1+NE1*i2+faceNECount;
1149 node0=(numDOFLocal-2)+i1*numDOFLocal+numDOFLocal*N1*i2 ;
1150
1151 out->FaceElements->Id[k]=idCount++;
1152 out->FaceElements->Tag[k]=2;
1153 out->FaceElements->Color[k]=0;//(i2%2)+2*(i1%2)+12;
1154
1155 if (useElementsOnFace) {
1156 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
1157 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal+1;
1158 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1159 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1160 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1161 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal;
1162 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1163 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal*N1;
1164 } else {
1165 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
1166 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal+1;
1167 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1168 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1169 }
1170 }
1171 }
1172 totalNECount+=NE1*NE2;
1173 faceNECount+=NE1*NE2;
1174 }
1175 }
1176 if (!periodic[1]) {
1177 /* ** elements on boundary 010 (x2=0): */
1178
1179 #pragma omp parallel for private(i0,i2,k,node0)
1180 for (i2=0;i2<NE2;i2++) {
1181 for (i0=0;i0<numElementsInternal;i0++) {
1182 k=i0+numElementsInternal*i2+faceNECount;
1183 node0=i0+numDOFLocal*N1*i2;
1184
1185 out->FaceElements->Id[k]=idCount++;
1186 out->FaceElements->Tag[k]=10;
1187 out->FaceElements->Color[k]=0;//(i0%2)+2*(i2%2)+16;
1188
1189 if (useElementsOnFace) {
1190 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1191 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
1192 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*numDOFLocal+1;
1193 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;
1194 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal;
1195 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal+1;
1196 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal+1;
1197 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;
1198 } else {
1199 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1200 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
1201 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*numDOFLocal+1;
1202 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;
1203 }
1204 }
1205 }
1206 totalNECount+=numElementsInternal*NE2;
1207 faceNECount+=numElementsInternal*NE2;
1208
1209 /* ** elements on boundary 020 (x2=1): */
1210
1211 #pragma omp parallel for private(i0,i2,k,node0)
1212 for (i2=0;i2<NE2;i2++) {
1213 for (i0=0;i0<numElementsInternal;i0++) {
1214 k=i0+numElementsInternal*i2+faceNECount;
1215 node0=i0+(NE1-1)*numDOFLocal+numDOFLocal*N1*i2;
1216
1217 out->FaceElements->Tag[k]=20;
1218 out->FaceElements->Id[k]=idCount++;
1219 out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
1220
1221 if (useElementsOnFace) {
1222 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;
1223 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1224 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1225 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal+1;
1226 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1227 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1;
1228 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1229 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
1230 } else {
1231 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;
1232 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1233 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1234 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal+1;
1235 }
1236
1237 }
1238 }
1239 totalNECount+=numElementsInternal*NE2;
1240 faceNECount+=numElementsInternal*NE2;
1241 }
1242 out->FaceElements->elementDistribution->numInternal = faceNECount;
1243
1244 /* now do the boundary face elements */
1245 /* LHS */
1246 if( !domLeft )
1247 {
1248 if( !periodic[2] ) {
1249
1250 /* x3=0 */
1251 for( i1=0; i1<NE1; i1++ )
1252 {
1253 k = i1+faceNECount;
1254 node0 = i1*numNodesLocal;
1255 node1 = numNodesLocal*N1*N2 + i1;
1256
1257 out->FaceElements->Tag[k]=200;
1258 out->FaceElements->Id[k]=idCount++;
1259 out->FaceElements->Color[k]=0;
1260
1261 if( useElementsOnFace )
1262 {
1263 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1264 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;
1265 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal;
1266 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
1267 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+N1;
1268 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1+1;
1269 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1270 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numNodesLocal*N1;
1271 } else {
1272 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1273 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;
1274 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal;
1275 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
1276 }
1277 }
1278 faceNECount += NE1;
1279 totalNECount += NE1;
1280
1281 /* x3=1 */
1282 for( i1=0; i1<NE1; i1++ )
1283 {
1284 k = i1+faceNECount;
1285 node0 = numNodesLocal*N1*(NE2-1) + i1*numNodesLocal;
1286 node1 = numNodesLocal*N1*N2 + i1 + (NE2-1)*N1;
1287
1288 out->FaceElements->Tag[k]=200;
1289 out->FaceElements->Id[k]=idCount++;
1290 out->FaceElements->Color[k]=0;
1291
1292 if( useElementsOnFace )
1293 {
1294 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;
1295 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal*N1;
1296 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1297 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;
1298 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;
1299 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0;
1300 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal;
1301 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+1;
1302 } else {
1303 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;
1304 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal*N1;
1305 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1306 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;
1307 }
1308 }
1309 faceNECount += NE1;
1310 totalNECount += NE1;
1311 }
1312
1313 if( !periodic[1] ) {
1314 /* x2=0 */
1315 for( i2=0; i2<NE2; i2++ )
1316 {
1317 k = i2+faceNECount;
1318 node0 = i2*numNodesLocal*N1;
1319 node1 = numNodesLocal*N1*N2 + i2*N1;
1320
1321 out->FaceElements->Tag[k]=20;
1322 out->FaceElements->Id[k]=idCount++;
1323 out->FaceElements->Color[k]=0;
1324
1325 if( useElementsOnFace )
1326 {
1327 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1328 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
1329 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1;
1330 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;
1331 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+1;
1332 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numNodesLocal;
1333 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1334 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+N1+1;
1335 } else {
1336 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1337 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
1338 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1;
1339 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;
1340 }
1341 }
1342 faceNECount += NE2;
1343 totalNECount += NE2;
1344
1345 /* x2=1 */
1346 for( i2=0; i2<NE2; i2++ )
1347 {
1348 k = i2+faceNECount;
1349 node0 = i2*numNodesLocal*N1 + numNodesLocal*(NE1-1);
1350 node1 = numNodesLocal*N1*N2 + i2*N1 + (NE1-1);
1351
1352 out->FaceElements->Tag[k]=20;
1353 out->FaceElements->Id[k]=idCount++;
1354 out->FaceElements->Color[k]=0;
1355
1356 if( useElementsOnFace )
1357 {
1358 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;
1359 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0;
1360 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal*N1;
1361 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1;
1362 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;
1363 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numNodesLocal;
1364 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1365 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;
1366 } else {
1367 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;
1368 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numNodesLocal;
1369 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1370 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;
1371 }
1372 }
1373 faceNECount += NE2;
1374 totalNECount += NE2;
1375 }
1376 }
1377
1378 /* RHS */
1379 if( !domRight || periodicLocal[1] )
1380 {
1381 /* the case of left hand boundary domain and periodic boundary condition on its left hand boundary */
1382 if( domLeft && periodic[0] ){
1383 if( !periodic[2] ) {
1384
1385 /* x3=0 */
1386 for( i1=0; i1<NE1; i1++ )
1387 {
1388 k = i1+faceNECount;
1389 node0 = numDOFLocal*N1*N2 + i1;
1390 node1 = numNodesLocal*N1*N2 + i1;
1391
1392 out->FaceElements->Tag[k]=200;
1393 out->FaceElements->Id[k]=idCount++;
1394 out->FaceElements->Color[k]=0;
1395
1396 if( useElementsOnFace )
1397 {
1398 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1399 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;
1400 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
1401 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
1402 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+N1;
1403 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1+1;
1404 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1+1;
1405 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1;
1406 } else {
1407 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1408 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;
1409 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
1410 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
1411 }
1412 }
1413 faceNECount += NE1;
1414 totalNECount += NE1;
1415
1416 /* x3=1 */
1417 for( i1=0; i1<NE1; i1++ )
1418 {
1419 k = i1+faceNECount;
1420 node0 = numDOFLocal*N1*N2 + i1 + (NE2-1)*N1;
1421 node1 = numNodesLocal*N1*N2 + i1 + (NE2-1)*N1;
1422
1423 out->FaceElements->Tag[k]=200;
1424 out->FaceElements->Id[k]=idCount++;
1425 out->FaceElements->Color[k]=0;
1426
1427 if( useElementsOnFace )
1428 {
1429 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;
1430 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+1;
1431 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+1;
1432 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0;
1433 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;
1434 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;
1435 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;
1436 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1;
1437 } else {
1438 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;
1439 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;
1440 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;
1441 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1;
1442 }
1443 }
1444 faceNECount += NE1;
1445 totalNECount += NE1;
1446 }
1447
1448 if( !periodic[1] ) {
1449 /* x2=0 */
1450 for( i2=0; i2<NE2; i2++ )
1451 {
1452 k = i2+faceNECount;
1453 node0 = numDOFLocal*N1*N2 + i2*N1;
1454 node1 = numNodesLocal*N1*N2 + i2*N1;
1455
1456 out->FaceElements->Tag[k]=20;
1457 out->FaceElements->Id[k]=idCount++;
1458 out->FaceElements->Color[k]=0;
1459
1460 if( useElementsOnFace )
1461 {
1462 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1463 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
1464 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1;
1465 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;
1466 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+1;
1467 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
1468 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1+1;
1469 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+N1+1;
1470 } else {
1471 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1472 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
1473 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1;
1474 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;
1475 }
1476 }
1477 faceNECount += NE2;
1478 totalNECount += NE2;
1479
1480 /* x2=1 */
1481 for( i2=0; i2<NE2; i2++ )
1482 {
1483 k = i2+faceNECount;
1484 node0 = numDOFLocal*N1*N2 + i2*N1 + N1-2;
1485 node1 = numNodesLocal*N1*N2 + i2*N1 + N1-2;
1486
1487 out->FaceElements->Tag[k]=20;
1488 out->FaceElements->Id[k]=idCount++;
1489 out->FaceElements->Color[k]=0;
1490
1491 if( useElementsOnFace )
1492 {
1493 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;
1494 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0;
1495 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1;
1496 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1;
1497 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;
1498 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
1499 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;
1500 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;
1501 } else {
1502 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;
1503 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
1504 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;
1505 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;
1506 }
1507 }
1508 faceNECount += NE2;
1509 totalNECount += NE2;
1510 }
1511
1512 }
1513 if( !periodic[2] ) {
1514 /* x3=0 */
1515 for( i1=0; i1<NE1; i1++ )
1516 {
1517 k = i1+faceNECount;
1518 node0 = numDOFLocal*(i1+1) - 1;
1519 node1 = (numNodesLocal+periodicLocal[0]+periodicLocal[1])*N1*N2 + i1 + domInternal*N1*N2;
1520
1521 out->FaceElements->Tag[k]=200;
1522 out->FaceElements->Id[k]=idCount++;
1523 out->FaceElements->Color[k]=0;
1524
1525 if( useElementsOnFace )
1526 {
1527 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1528 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;
1529 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+1;
1530 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1;
1531 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal*N1;
1532 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1533 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+N1+1;
1534 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+N1;
1535 } else {
1536 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1537 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;
1538 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+1;
1539 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1;
1540 }
1541 }
1542 faceNECount += NE1;
1543 totalNECount += NE1;
1544
1545 /* x3=1 */
1546 for( i1=0; i1<NE1; i1++ )
1547 {
1548 k = i1+faceNECount;
1549 node0 = numDOFLocal*N1*(NE2-1) + (i1+1)*numDOFLocal - 1;
1550 node1 = numNodesLocal*N1*N2 + i1 + (NE2-1)*N1 + (domInternal+periodicLocal[1]+periodicLocal[0])*N1*N2;
1551
1552 out->FaceElements->Tag[k]=200;
1553 out->FaceElements->Id[k]=idCount++;
1554 out->FaceElements->Color[k]=0;
1555
1556 if( useElementsOnFace )
1557 {
1558 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1559 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal;
1560 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+1;
1561 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1;
1562 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal*N1;
1563 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1564 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;
1565 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1;
1566 } else {
1567 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal*N1;
1568 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1569 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;
1570 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1;
1571 }
1572 }
1573 faceNECount += NE1;
1574 totalNECount += NE1;
1575 }
1576 if( !periodic[1] ) {
1577 /* x2=0 */
1578 for( i2=0; i2<NE2; i2++ )
1579 {
1580 k = i2+faceNECount;
1581 node0 = N1*numDOFLocal*i2 + numDOFLocal - 1;
1582 node1 = numNodesLocal*N1*N2 + i2*N1 + (domInternal+periodicLocal[0]+periodicLocal[1])*N1*N2;
1583
1584 out->FaceElements->Tag[k]=20;
1585 out->FaceElements->Id[k]=idCount++;
1586 out->FaceElements->Color[k]=0;
1587
1588 if( useElementsOnFace )
1589 {
1590 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1591 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1;
1592 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1;
1593 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;
1594 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal;
1595 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+1;
1596 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+N1+1;
1597 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;
1598 } else {
1599 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1600 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1;
1601 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1;
1602 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;
1603 }
1604 }
1605 faceNECount += NE2;
1606 totalNECount += NE2;
1607
1608 /* x2=1 */
1609 for( i2=0; i2<NE2; i2++ )
1610 {
1611 k = i2+faceNECount;
1612 node0 = numDOFLocal*N1*i2 + NE1*numDOFLocal - 1;
1613 node1 = numNodesLocal*N1*N2 + i2*N1 + (NE1-1) + (domInternal+periodicLocal[0]+periodicLocal[1])*N1*N2;
1614
1615 out->FaceElements->Tag[k]=20;
1616 out->FaceElements->Id[k]=idCount++;
1617 out->FaceElements->Color[k]=0;
1618
1619 if( useElementsOnFace ){
1620 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1621 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1;
1622 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+N1;
1623 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N1*numDOFLocal;
1624 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;
1625 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+1;
1626 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;
1627 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;
1628 } else {
1629 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;
1630 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+1;
1631 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;
1632 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;
1633 }
1634 }
1635 faceNECount += NE2;
1636 totalNECount += NE2;
1637 }
1638 }
1639 out->FaceElements->minColor=0;
1640 out->FaceElements->maxColor=0;//23;
1641
1642 out->FaceElements->elementDistribution->numBoundary = faceNECount - out->FaceElements->elementDistribution->numInternal;
1643 out->FaceElements->elementDistribution->numLocal = faceNECount;
1644
1645
1646 /* setup distribution info for other elements */
1647 out->ContactElements->elementDistribution->numLocal = out->ContactElements->elementDistribution->numInternal = out->ContactElements->elementDistribution->numInternal = 0;
1648 out->Points->elementDistribution->numLocal = out->Points->elementDistribution->numInternal = out->Points->elementDistribution->numInternal = 0;
1649
1650 /* condense the nodes: */
1651 Finley_Mesh_resolveNodeIds( out );
1652
1653 /* setup the CommBuffer */
1654 Finley_NodeDistribution_formCommBuffer( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );
1655 if ( !Finley_MPI_noError( mpi_info )) {
1656 if( Finley_noError() )
1657 Finley_setError( PASO_MPI_ERROR, "Error on another MPI process" );
1658 Paso_MPIInfo_dealloc( mpi_info );
1659 Finley_Mesh_dealloc(out);
1660 return NULL;
1661 }
1662
1663 Finley_NodeDistribution_calculateIndexExternal( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );
1664
1665 /* prepare mesh for further calculatuions:*/
1666 Finley_Mesh_prepare(out) ;
1667
1668 // print_mesh_statistics( out );
1669
1670 #ifdef Finley_TRACE
1671 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
1672 #endif
1673
1674 if (! Finley_noError()) {
1675 Finley_Mesh_dealloc(out);
1676 return NULL;
1677 }
1678
1679 return out;
1680 }
1681 #endif
1682

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
  ViewVC Help
Powered by ViewVC 1.1.26