/[escript]/branches/domexper/dudley/src/Mesh_tet4.c
ViewVC logotype

Contents of /branches/domexper/dudley/src/Mesh_tet4.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 787 - (show annotations)
Wed Jul 26 01:46:45 2006 UTC (12 years, 10 months ago) by bcumming
Original Path: trunk/finley/src/Mesh_hex8.c
File MIME type: text/plain
File size: 36667 byte(s)
MPI update
Each element (normal elements, faceElements, ContactElements and points)
are now assigned a unique global id to streamline per-element
calculations and file IO of element data.



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 #endif
116
117 #ifdef PASO_MPI
118 Finley_Mesh* Finley_RectangularMesh_Hex8_singleCPU(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace,Paso_MPIInfo *mpi_info)
119 #else
120 Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
121 #endif
122 {
123 dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES,M0,M1,M2;
124 index_t node0;
125 Finley_Mesh* out;
126 char name[50];
127 double time0=Finley_timer();
128 NE0=MAX(1,numElements[0]);
129 NE1=MAX(1,numElements[1]);
130 NE2=MAX(1,numElements[2]);
131 N0=NE0+1;
132 N1=NE1+1;
133 N2=NE2+1;
134
135 if (N0<=MIN(N1,N2)) {
136 if (N1 <= N2) {
137 M0=1;
138 M1=N0;
139 M2=N0*N1;
140 } else {
141 M0=1;
142 M2=N0;
143 M1=N0*N2;
144 }
145 } else if (N1<=MIN(N2,N0)) {
146 if (N2 <= N0) {
147 M1=1;
148 M2=N1;
149 M0=N2*N1;
150 } else {
151 M1=1;
152 M0=N1;
153 M2=N1*N0;
154 }
155 } else {
156 if (N0 <= N1) {
157 M2=1;
158 M0=N2;
159 M1=N2*N0;
160 } else {
161 M2=1;
162 M1=N2;
163 M0=N1*N2;
164 }
165 }
166
167
168 NFaceElements=0;
169 if (!periodic[0]) {
170 NDOF0=N0;
171 NFaceElements+=2*NE1*NE2;
172 } else {
173 NDOF0=N0-1;
174 }
175 if (!periodic[1]) {
176 NDOF1=N1;
177 NFaceElements+=2*NE0*NE2;
178 } else {
179 NDOF1=N1-1;
180 }
181 if (!periodic[2]) {
182 NDOF2=N2;
183 NFaceElements+=2*NE0*NE1;
184 } else {
185 NDOF2=N2-1;
186 }
187
188 /* allocate mesh: */
189
190 sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
191
192 #ifndef PASO_MPI
193 out=Finley_Mesh_alloc(name,3,order);
194 #else
195 out=Finley_Mesh_alloc(name,3,order,mpi_info);
196 #endif
197 if (! Finley_noError()) return NULL;
198
199 #ifdef PASO_MPI
200 out->Elements=Finley_ElementFile_alloc(Hex8,out->order,mpi_info);
201 if (useElementsOnFace) {
202 out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,mpi_info);
203 out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order,mpi_info);
204 } else {
205 out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order,mpi_info);
206 out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order,mpi_info);
207 }
208 out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
209 #else
210 out->Elements=Finley_ElementFile_alloc(Hex8,out->order);
211 if (useElementsOnFace) {
212 out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order);
213 out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order);
214 } else {
215 out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order);
216 out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order);
217 }
218 out->Points=Finley_ElementFile_alloc(Point1,out->order);
219 #endif
220 if (! Finley_noError()) {
221 Finley_Mesh_dealloc(out);
222 return NULL;
223 }
224
225
226 /* allocate tables: */
227 Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);
228 #ifdef PASO_MPI
229 Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0*NDOF1*NDOF2, 0, 0 );
230 #endif
231 Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);
232 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
233 if (! Finley_noError()) {
234 Finley_Mesh_dealloc(out);
235 return NULL;
236 }
237
238 /* set nodes: */
239
240 #pragma omp parallel for private(i0,i1,i2,k)
241 for (i2=0;i2<N2;i2++) {
242 for (i1=0;i1<N1;i1++) {
243 for (i0=0;i0<N0;i0++) {
244 k=M0*i0+M1*i1+M2*i2;
245 out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0)/DBLE(N0-1)*Length[0];
246 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
247 out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
248 out->Nodes->Id[k]=i0+N0*i1+N0*N1*i2;
249 out->Nodes->Tag[k]=0;
250 out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1) +M2*(i2%NDOF2);
251 #ifdef PASO_MPI
252 out->Nodes->Dom[k]=NODE_INTERNAL;
253 #endif
254 }
255 }
256 }
257 /* tags for the faces: */
258 /* tags for the faces: */
259 if (!periodic[2]) {
260 for (i1=0;i1<N1;i1++) {
261 for (i0=0;i0<N0;i0++) {
262 out->Nodes->Tag[M0*i0+M1*i1+M2*0]+=100;
263 out->Nodes->Tag[M0*i0+M1*i1+M2*(N2-1)]+=200;
264 }
265 }
266 }
267 if (!periodic[1]) {
268 for (i2=0;i2<N2;i2++) {
269 for (i0=0;i0<N0;i0++) {
270 out->Nodes->Tag[M0*i0+M1*0+M2*i2]+=10;
271 out->Nodes->Tag[M0*i0+M1*(N1-1)+M2*i2]+=20;
272 }
273 }
274 }
275 if (!periodic[0]) {
276 for (i2=0;i2<N2;i2++) {
277 for (i1=0;i1<N1;i1++) {
278 out->Nodes->Tag[M0*0+M1*i1+M2*i2]+=1;
279 out->Nodes->Tag[M0*(N0-1)+M1*i1+M2*i2]+=2;
280 }
281 }
282 }
283
284 /* set the elements: */
285
286 #pragma omp parallel for private(i0,i1,i2,k,node0)
287 for (i2=0;i2<NE2;i2++) {
288 for (i1=0;i1<NE1;i1++) {
289 for (i0=0;i0<NE0;i0++) {
290 k=i0+NE0*i1+NE0*NE1*i2;
291 node0=i0+i1*N0+N0*N1*i2;
292
293 out->Elements->Id[k]=k;
294 out->Elements->Tag[k]=0;
295 out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
296 #ifdef PASO_MPI
297 out->Elements->Dom[k]=ELEMENT_INTERNAL;
298 #endif
299
300 out->Elements->Nodes[INDEX2(0,k,8)]=node0;
301 out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
302 out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0+1;
303 out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0;
304 out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0*N1;
305 out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0*N1+1;
306 out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0*N1+N0+1;
307 out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0*N1+N0;
308
309 }
310 }
311 }
312 out->Elements->minColor=0;
313 out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
314
315 /* face elements: */
316
317 if (useElementsOnFace) {
318 NUMNODES=8;
319 } else {
320 NUMNODES=4;
321 }
322 totalNECount=NE0*NE1*NE2;
323 faceNECount=0;
324
325 /* these are the quadrilateral elements on boundary 1 (x3=0): */
326 if (!periodic[2]) {
327 /* ** elements on boundary 100 (x3=0): */
328
329 #pragma omp parallel for private(i0,i1,k,node0)
330 for (i1=0;i1<NE1;i1++) {
331 for (i0=0;i0<NE0;i0++) {
332 k=i0+NE0*i1+faceNECount;
333 node0=i0+i1*N0;
334
335 out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
336 out->FaceElements->Tag[k]=100;
337 out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
338 #ifdef PASO_MPI
339 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
340 #endif
341
342 if (useElementsOnFace) {
343 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
344 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
345 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
346 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
347 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0*N1;
348 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+N0;
349 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
350 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1+1;
351 } else {
352 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
353 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
354 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
355 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
356 }
357 }
358 }
359 totalNECount+=NE1*NE0;
360 faceNECount+=NE1*NE0;
361
362 /* ** elements on boundary 200 (x3=1) */
363
364 #pragma omp parallel for private(i0,i1,k,node0)
365 for (i1=0;i1<NE1;i1++) {
366 for (i0=0;i0<NE0;i0++) {
367 k=i0+NE0*i1+faceNECount;
368 node0=i0+i1*N0+N0*N1*(NE2-1);
369
370 out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
371 out->FaceElements->Tag[k]=200;
372 out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
373 #ifdef PASO_MPI
374 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
375 #endif
376
377 if (useElementsOnFace) {
378 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
379 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
380 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
381 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
382 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
383 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
384 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+1;
385 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;
386 } else {
387 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
388 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
389 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
390 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
391 }
392
393
394 }
395 }
396 totalNECount+=NE1*NE0;
397 faceNECount+=NE1*NE0;
398 }
399 if (!periodic[0]) {
400 /* ** elements on boundary 001 (x1=0): */
401
402 #pragma omp parallel for private(i1,i2,k,node0)
403 for (i2=0;i2<NE2;i2++) {
404 for (i1=0;i1<NE1;i1++) {
405 k=i1+NE1*i2+faceNECount;
406 node0=i1*N0+N0*N1*i2;
407
408 out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
409 out->FaceElements->Tag[k]=1;
410 out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
411 #ifdef PASO_MPI
412 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
413 #endif
414
415 if (useElementsOnFace) {
416 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
417 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
418 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
419 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
420 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
421 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+1;
422 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
423 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0+1;
424 } else {
425 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
426 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
427 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
428 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
429 }
430 }
431 }
432 totalNECount+=NE1*NE2;
433 faceNECount+=NE1*NE2;
434
435 /* ** elements on boundary 002 (x1=1): */
436
437 #pragma omp parallel for private(i1,i2,k,node0)
438 for (i2=0;i2<NE2;i2++) {
439 for (i1=0;i1<NE1;i1++) {
440 k=i1+NE1*i2+faceNECount;
441 node0=(NE0-1)+i1*N0+N0*N1*i2 ;
442
443 out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
444 out->FaceElements->Tag[k]=2;
445 out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
446 #ifdef PASO_MPI
447 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
448 #endif
449
450 if (useElementsOnFace) {
451 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
452 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
453 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
454 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
455 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
456 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0;
457 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0;
458 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1;
459 } else {
460 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
461 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
462 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
463 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
464 }
465 }
466 }
467 totalNECount+=NE1*NE2;
468 faceNECount+=NE1*NE2;
469 }
470 if (!periodic[1]) {
471 /* ** elements on boundary 010 (x2=0): */
472
473 #pragma omp parallel for private(i0,i2,k,node0)
474 for (i2=0;i2<NE2;i2++) {
475 for (i0=0;i0<NE0;i0++) {
476 k=i0+NE0*i2+faceNECount;
477 node0=i0+N0*N1*i2;
478
479 out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
480 out->FaceElements->Tag[k]=10;
481 out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;
482 #ifdef PASO_MPI
483 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
484 #endif
485
486 if (useElementsOnFace) {
487 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
488 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
489 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
490 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
491 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;
492 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0+1;
493 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1*N0+N0+1;
494 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*N0+N0;
495 } else {
496 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
497 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
498 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
499 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
500 }
501 }
502 }
503 totalNECount+=NE0*NE2;
504 faceNECount+=NE0*NE2;
505
506 /* ** elements on boundary 020 (x2=1): */
507
508 #pragma omp parallel for private(i0,i2,k,node0)
509 for (i2=0;i2<NE2;i2++) {
510 for (i0=0;i0<NE0;i0++) {
511 k=i0+NE0*i2+faceNECount;
512 node0=i0+(NE1-1)*N0+N0*N1*i2;
513
514 out->FaceElements->Tag[k]=20;
515 out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
516 out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
517 #ifdef PASO_MPI
518 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
519 #endif
520
521 if (useElementsOnFace) {
522 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
523 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
524 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
525 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
526 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
527 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1;
528 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+1;
529 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
530 } else {
531 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
532 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
533 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
534 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
535 }
536
537 }
538 }
539 totalNECount+=NE0*NE2;
540 faceNECount+=NE0*NE2;
541 }
542 out->FaceElements->minColor=0;
543 out->FaceElements->maxColor=23;
544
545 #ifdef PASO_MPI
546 Finley_ElementFile_setDomainFlags( out->Elements );
547 Finley_ElementFile_setDomainFlags( out->FaceElements );
548 Finley_ElementFile_setDomainFlags( out->ContactElements );
549 Finley_ElementFile_setDomainFlags( out->Points );
550
551 /* reorder the degrees of freedom */
552 Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
553 #endif
554
555 /* condense the nodes: */
556 Finley_Mesh_resolveNodeIds(out);
557
558 /* prepare mesh for further calculatuions:*/
559 Finley_Mesh_prepare(out) ;
560
561 #ifdef Finley_TRACE
562 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
563 #endif
564
565 if (! Finley_noError()) {
566 Finley_Mesh_dealloc(out);
567 return NULL;
568 }
569 return out;
570 }
571
572 #ifdef PASO_MPI
573 Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
574 {
575 dim_t N0,N1,N2,N0t,NDOF0t,NE0,NE1,NE2,i0,i1,i2,kk,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES;//,M0,M1,M2;
576 dim_t idCount, NE0_local, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal;
577 bool_t dom_left, dom_right, dom_internal;
578 index_t firstNode=0, DOFcount=0, node0, node1, node2;
579 index_t targetDomain=-1, firstNodeConstruct, j;
580 bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE, boundaryLeft=FALSE, boundaryRight=FALSE;
581 index_t *indexBackward=NULL, *indexForward=NULL,*facePerm=NULL, *forwardDOF=NULL, *backwardDOF=NULL;
582 Finley_Mesh* out;
583
584 char name[50];
585 Paso_MPIInfo *mpi_info = NULL;
586 double time0=Finley_timer();
587
588 index_t face0[] = {3, 0, 4, 7, 2, 1, 5, 6};
589 index_t face1[] = {1, 2, 6, 5, 0, 3, 7, 4};
590 index_t face2[] = {0, 3, 2, 1, 4, 7, 6, 5};
591 index_t face3[] = {4, 5, 6, 7, 0, 1, 2, 3};
592 index_t face4[] = {0, 1, 5, 4, 3, 2, 6, 7};
593 index_t face5[] = {3, 7, 6, 2, 0, 4, 5, 1};
594 NE0=MAX(1,numElements[0]);
595 NE1=MAX(1,numElements[1]);
596 NE2=MAX(1,numElements[2]);
597 N0=NE0+1;
598 N1=NE1+1;
599 N2=NE2+1;
600
601
602 /* get MPI information */
603 mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
604 if (! Finley_noError())
605 return NULL;
606
607 /* use the serial version to generate the mesh for the 1-CPU case */
608 if( mpi_info->size==1 )
609 {
610 out = Finley_RectangularMesh_Hex8_singleCPU( numElements, Length, periodic, order, useElementsOnFace, mpi_info );
611 return out;
612 }
613
614 if( mpi_info->rank==0 )
615 domLeft = TRUE;
616 if( mpi_info->rank==mpi_info->size-1 )
617 domRight = TRUE;
618 if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )
619 domInternal = TRUE;
620
621 /* dimensions of the local subdomain */
622 domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal );
623
624 /* count Degrees of Freedom along each axis */
625 NDOF0 = N0 - periodic[0];
626 NDOF1 = N1 - periodic[1];
627 NDOF2 = N2 - periodic[2];
628
629 /* count face elements */
630 /* internal face elements */
631 NFaceElements = 0;
632 if( !periodic[0] )
633 NFaceElements += (domLeft+domRight)*NE1*NE2;
634 if( !periodic[1] )
635 NFaceElements += 2*numElementsLocal*NE2;
636 if( !periodic[2] )
637 NFaceElements += 2*numElementsLocal*NE1;
638
639 boundaryLeft = !domLeft || periodicLocal[0];
640 boundaryRight = !domRight || periodicLocal[1];
641 N0t = numNodesLocal + boundaryRight + boundaryLeft;
642 NDOF0t = numDOFLocal + boundaryRight + boundaryLeft;
643 firstNodeConstruct = firstNode - boundaryLeft;
644 firstNodeConstruct = firstNodeConstruct<0 ? N0-2 : firstNodeConstruct;
645
646 /* allocate mesh: */
647 sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
648
649 out=Finley_Mesh_alloc(name,3,order,mpi_info);
650 if (! Finley_noError()) return NULL;
651
652 out->Elements=Finley_ElementFile_alloc(Hex8,out->order,mpi_info);
653 if (useElementsOnFace) {
654 out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,mpi_info);
655 out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order,mpi_info);
656 } else {
657 out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order,mpi_info);
658 out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order,mpi_info);
659 }
660 out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
661
662 if (! Finley_noError()) {
663 Finley_Mesh_dealloc(out);
664 return NULL;
665 }
666
667 /* allocate tables: */
668 Finley_NodeFile_allocTable(out->Nodes,N0t*N1*N2);
669 Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal*NDOF1*NDOF2, NDOF1*NDOF2*2, 0 );
670 Finley_ElementFile_allocTable(out->Elements,(numElementsLocal)*NE1*NE2);
671 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
672 if (! Finley_noError()) {
673 Finley_Mesh_dealloc(out);
674 return NULL;
675 }
676
677 /* set nodes: */
678 /* INTERNAL/BOUNDARY NODES */
679 k=0;
680 #pragma omp parallel for private(i0,i1,i2,k)
681 for (i2=0;i2<N2;i2++) {
682 for (i1=0;i1<N1;i1++) {
683 for (i0=0;i0<N0t;i0++,k++) {
684 out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE((i0+firstNodeConstruct) % N0)/DBLE(N0-1)*Length[0];
685 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
686 out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
687 out->Nodes->Id[k]=k;
688 out->Nodes->Tag[k]=0;
689 out->Nodes->degreeOfFreedom[k]=i0 + (i1%NDOF1)*N0t + (i2%NDOF2)*N0t*N1;
690 out->Nodes->Dom[k]=NODE_INTERNAL;
691 }
692 }
693 }
694
695 /* mark the nodes that reference external and boundary DOF as such */
696 if( boundaryLeft ){
697 for( i1=0; i1<N1; i1++ )
698 for( i2=0; i2<N2; i2++ ) {
699 out->Nodes->Dom[N1*N0t*i2+N0t*i1] = NODE_EXTERNAL;
700 out->Nodes->Dom[N1*N0t*i2+N0t*i1+1] = NODE_BOUNDARY;
701 }
702 }
703 if( boundaryRight ){
704 for( i1=0; i1<N1; i1++ )
705 for( i2=0; i2<N2; i2++ ) {
706 out->Nodes->Dom[N1*N0t*i2+N0t*(i1+1)-1] = NODE_EXTERNAL;
707 out->Nodes->Dom[N1*N0t*i2+N0t*(i1+1)-2] = NODE_BOUNDARY;
708 }
709 }
710 if( periodicLocal[0] ){
711 for( i1=0; i1<N1; i1++ )
712 for( i2=0; i2<N2; i2++ ) {
713 out->Nodes->degreeOfFreedom[N1*N0t*i2+i1*N0t+3] = out->Nodes->degreeOfFreedom[N1*N0t*i2+i1*N0t+2];
714 out->Nodes->Dom[N1*N0t*i2+N0t*i1+3] = NODE_BOUNDARY;
715 }
716 }
717
718 /* tag Nodes that are referenced by face elements */
719 if (!periodic[2]) {
720 for (i1=0;i1<N1;i1++) {
721 for (i0=0;i0<N0t;i0++) {
722 out->Nodes->Tag[i0 + N0t*i1]+=100;
723 out->Nodes->Tag[i0 + N0t*i1 + N0t*N1*(N2-1)]+=200;
724 }
725 }
726 }
727 if (!periodic[1]) {
728 for (i2=0;i2<N2;i2++) {
729 for (i0=0;i0<N0t;i0++) {
730 out->Nodes->Tag[i0 + i2*N1*N0t]+=10;
731 out->Nodes->Tag[i0 + (i2+1)*N1*N0t-N0t]+=20;
732 }
733 }
734 }
735 if (!periodic[0] && !domInternal ) {
736 for (i2=0;i2<N2;i2++) {
737 for (i1=0;i1<N1;i1++) {
738 if( domLeft )
739 out->Nodes->Tag[i1*N0t + i2*N0t*N1]+=1;
740 if( domRight )
741 out->Nodes->Tag[(i1+1)*N0t-1 + i2*N0t*N1]+=2;
742 }
743 }
744 }
745
746 /* form the boudary communication information */
747 forwardDOF = MEMALLOC(NDOF1*NDOF2,index_t);
748 backwardDOF = MEMALLOC(NDOF1*NDOF2,index_t);
749 if( !(mpi_info->size==2 && periodicLocal[0])){
750 if( boundaryLeft ) {
751 targetDomain = mpi_info->rank-1 < 0 ? mpi_info->size-1 : mpi_info->rank-1;
752 for( i2=0; i2<NDOF2; i2++ ){
753 for( i1=0; i1<NDOF1; i1++ ){
754 forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+1];
755 backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t];
756 }
757 }
758 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
759 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
760 }
761 if( boundaryRight ) {
762 targetDomain = mpi_info->rank+1 > mpi_info->size-1 ? 0 : mpi_info->rank+1;
763 for( i2=0; i2<NDOF2; i2++ ){
764 for( i1=0; i1<NDOF1; i1++ ){
765 forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-2];
766 backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-1];
767 }
768 }
769 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
770 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
771 }
772 } else{
773 /* periodic boundary conditions with 2 domains, need to change the order in which domain 0 passes boundary data */
774 targetDomain = 1;
775
776 for( i2=0; i2<NDOF2; i2++ ){
777 for( i1=0; i1<NDOF1; i1++ ){
778 forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-2];
779 backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-1];
780 }
781 }
782 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
783 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
784
785 for( i2=0; i2<NDOF2; i2++ ){
786 for( i1=0; i1<NDOF1; i1++ ){
787 forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+1];
788 backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t];
789 }
790 }
791 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
792 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
793 }
794 MEMFREE( forwardDOF );
795 MEMFREE( backwardDOF );
796 /* set the elements: */
797
798 /* INTERNAL elements */
799 k = 0;
800 #pragma omp parallel for private(i0,i1,i2,k,node0)
801 for (i2=0;i2<NE2;i2++) {
802 for (i1=0;i1<NE1;i1++) {
803 for (i0=0;i0<numElementsLocal;i0++,k++) {
804 node0 = (periodicLocal[0] && !i0) ? i1*N0t + i2*N1*N0t : i1*N0t + i2*N1*N0t + i0 + periodicLocal[0];
805
806 out->Elements->Id[k]=((firstNodeConstruct+i0)%NE0)*NE1*NE2 + NE1*i2 + i1;
807 out->Elements->Tag[k]=0;
808 out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
809 out->Elements->Dom[k]=ELEMENT_INTERNAL;
810
811 out->Elements->Nodes[INDEX2(0,k,8)]=node0;
812 out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
813 out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0t+1;
814 out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0t;
815 out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0t*N1;
816 out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0t*N1+1;
817 out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0t*N1+N0t+1;
818 out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0t*N1+N0t;
819
820 }
821 }
822 }
823 out->Elements->minColor=0;
824 out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
825 if( boundaryLeft )
826 for( i2=0; i2<NE2; i2++ )
827 for( i1=0; i1<NE1; i1++ )
828 out->Elements->Dom[i2*NE1*numElementsLocal+i1*numElementsLocal]=ELEMENT_BOUNDARY;
829 if( boundaryRight )
830 for( i2=0; i2<NE2; i2++ )
831 for( i1=0; i1<NE1; i1++ )
832 out->Elements->Dom[i2*NE1*numElementsLocal+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
833
834 Finley_ElementFile_setDomainFlags( out->Elements );
835
836 /* face elements: */
837 if (useElementsOnFace) {
838 NUMNODES=8;
839 } else {
840 NUMNODES=4;
841 }
842 totalNECount=out->Elements->numElements;
843 faceNECount=0;
844 idCount = totalNECount;
845
846 /* these are the quadrilateral elements on boundary 1 (x3=0): */
847 numElementsInternal = numElementsLocal-nodesExternal[0]-nodesExternal[1];
848 if (!periodic[2]) {
849 /* elements on boundary 100 (x3=0): */
850
851 #pragma omp parallel for private(i0,i1,k)
852 for (i1=0;i1<NE1;i1++) {
853 for (i0=0; i0<numElementsLocal; i0++) {
854 k=i0+numElementsLocal*i1+faceNECount;
855 kk=i0 + i1*numElementsLocal;
856 facePerm = face2;
857
858 out->FaceElements->Id[k]=idCount++;
859 out->FaceElements->Tag[k]=100;
860 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
861 out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
862
863 for( j=0; j<NUMNODES; j++ )
864 out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
865 }
866 }
867 if( boundaryLeft ){
868 for( i1=0; i1<NE1; i1++ )
869 out->FaceElements->Dom[faceNECount+i1*numElementsLocal]=ELEMENT_BOUNDARY;
870 if( periodicLocal[0] )
871 for( i1=0; i1<NE1; i1++ )
872 out->FaceElements->Dom[faceNECount+i1*numElementsLocal+1]=ELEMENT_BOUNDARY;
873 }
874 if( boundaryRight )
875 for( i1=0; i1<NE1; i1++ )
876 out->FaceElements->Dom[faceNECount+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
877 totalNECount+=NE1*numElementsLocal;
878 faceNECount+=NE1*numElementsLocal;
879
880 /* ** elements on boundary 200 (x3=1) */
881
882 #pragma omp parallel for private(i0,i1,k)
883 for (i1=0;i1<NE1;i1++) {
884 for (i0=0;i0<numElementsLocal;i0++) {
885 k=i0+numElementsLocal*i1+faceNECount;
886 kk=i0+i1*numElementsLocal+numElementsLocal*NE1*(NE2-1);
887 facePerm = face3;
888
889 out->FaceElements->Id[k]=idCount++;
890 out->FaceElements->Tag[k]=200;
891 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
892 out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
893
894 for( j=0; j<NUMNODES; j++ )
895 out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
896 }
897 }
898 if( boundaryLeft ){
899 for( i1=0; i1<NE1; i1++ )
900 out->FaceElements->Dom[faceNECount+i1*numElementsLocal]=ELEMENT_BOUNDARY;
901 if( periodicLocal[0] )
902 for( i1=0; i1<NE1; i1++ )
903 out->FaceElements->Dom[faceNECount+i1*numElementsLocal+1]=ELEMENT_BOUNDARY;
904 }
905 if( boundaryRight )
906 for( i1=0; i1<NE1; i1++ )
907 out->FaceElements->Dom[faceNECount+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
908 totalNECount+=NE1*numElementsLocal;
909 faceNECount+=NE1*numElementsLocal;
910 }
911 if (!periodic[0] && !domInternal) {
912 /* ** elements on boundary 001 (x1=0): */
913 if( domLeft ){
914 #pragma omp parallel for private(i1,i2,k)
915 for (i2=0;i2<NE2;i2++) {
916 for (i1=0;i1<NE1;i1++) {
917 k=i1+NE1*i2+faceNECount;
918 kk=i1*numElementsLocal + i2*numElementsLocal*NE1;
919 facePerm = face0;
920
921 out->FaceElements->Id[k]=idCount++;
922 out->FaceElements->Tag[k]=1;
923 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
924 out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
925
926 for( j=0; j<NUMNODES; j++ )
927 out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
928 }
929 }
930 totalNECount+=NE1*NE2;
931 faceNECount+=NE1*NE2;
932 }
933 /* ** elements on boundary 002 (x1=1): */
934 if( domRight ) {
935 #pragma omp parallel for private(i1,i2,k)
936 for (i2=0;i2<NE2;i2++) {
937 for (i1=0;i1<NE1;i1++) {
938 k=i1+NE1*i2+faceNECount;
939 kk=(i1+1)*numElementsLocal + i2*numElementsLocal*NE1 - 1;
940 facePerm = face1;
941
942 out->FaceElements->Id[k]=idCount++;
943 out->FaceElements->Tag[k]=2;
944 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
945 out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
946
947 for( j=0; j<NUMNODES; j++ )
948 out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
949 }
950 }
951 totalNECount+=NE1*NE2;
952 faceNECount+=NE1*NE2;
953 }
954 }
955 if (!periodic[1]) {
956 /* ** elements on boundary 010 (x2=0): */
957
958 #pragma omp parallel for private(i0,i2,k)
959 for (i2=0;i2<NE2;i2++) {
960 for (i0=0;i0<numElementsLocal;i0++) {
961 k=i0+numElementsLocal*i2+faceNECount;
962 kk=i0+numElementsLocal*NE1*i2;
963 facePerm = face4;
964
965 out->FaceElements->Id[k]=idCount++;
966 out->FaceElements->Tag[k]=10;
967 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
968 out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;
969
970 for( j=0; j<NUMNODES; j++ )
971 out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
972 }
973 }
974 if( boundaryLeft ){
975 for( i2=0; i2<NE2; i2++ )
976 out->FaceElements->Dom[faceNECount+i2*numElementsLocal]=ELEMENT_BOUNDARY;
977 if( periodicLocal[0] )
978 for( i2=0; i2<NE2; i2++ )
979 out->FaceElements->Dom[faceNECount+i2*numElementsLocal+1]=ELEMENT_BOUNDARY;
980 }
981 if( boundaryRight )
982 for( i2=0; i2<NE2; i2++ )
983 out->FaceElements->Dom[faceNECount+(i2+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
984 totalNECount+=numElementsLocal*NE2;
985 faceNECount+=numElementsLocal*NE2;
986
987 /* ** elements on boundary 020 (x2=1): */
988
989 #pragma omp parallel for private(i0,i2,k)
990 for (i2=0;i2<NE2;i2++) {
991 for (i0=0;i0<numElementsLocal;i0++) {
992 k=i0+numElementsLocal*i2+faceNECount;
993 kk=i0+numElementsLocal*NE1*(i2+1)-numElementsLocal;
994 facePerm = face5;
995
996 out->FaceElements->Tag[k]=20;
997 out->FaceElements->Id[k]=idCount++;
998 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
999 out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
1000
1001 for( j=0; j<NUMNODES; j++ )
1002 out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
1003 }
1004 }
1005 if( boundaryLeft ){
1006 for( i2=0; i2<NE2; i2++ )
1007 out->FaceElements->Dom[faceNECount+i2*numElementsLocal]=ELEMENT_BOUNDARY;
1008 if( periodicLocal[0] )
1009 for( i2=0; i2<NE2; i2++ )
1010 out->FaceElements->Dom[faceNECount+i2*numElementsLocal+1]=ELEMENT_BOUNDARY;
1011 }
1012 if( boundaryRight )
1013 for( i2=0; i2<NE2; i2++ )
1014 out->FaceElements->Dom[faceNECount+(i2+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
1015 totalNECount+=numElementsLocal*NE2;
1016 faceNECount+=numElementsLocal*NE2;
1017 }
1018 out->FaceElements->elementDistribution->numInternal = faceNECount;
1019
1020 out->FaceElements->minColor=0;
1021 out->FaceElements->maxColor=23;
1022 out->FaceElements->numElements=faceNECount;
1023
1024 Finley_ElementFile_setDomainFlags( out->FaceElements );
1025
1026 /* setup distribution info for other elements */
1027 Finley_ElementFile_setDomainFlags( out->ContactElements );
1028 Finley_ElementFile_setDomainFlags( out->Points );
1029
1030 /* reorder the degrees of freedom */
1031 Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
1032
1033 /* condense the nodes: */
1034 Finley_Mesh_resolveNodeIds(out);
1035 if( !Finley_MPI_noError(mpi_info) )
1036 {
1037 Paso_MPIInfo_dealloc( mpi_info );
1038 Finley_Mesh_dealloc(out);
1039 return NULL;
1040 }
1041
1042 /* prepare mesh for further calculatuions:*/
1043 Finley_Mesh_prepare(out);
1044 if( !Finley_MPI_noError(mpi_info) )
1045 {
1046 Paso_MPIInfo_dealloc( mpi_info );
1047 Finley_Mesh_dealloc(out);
1048 return NULL;
1049 }
1050
1051 /* free up memory */
1052 Paso_MPIInfo_dealloc( mpi_info );
1053
1054 #ifdef Finley_TRACE
1055 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
1056 #endif
1057
1058 return out;
1059 }
1060 #endif
1061

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26