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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1981 - (show annotations)
Thu Nov 6 05:27:33 2008 UTC (10 years, 8 months ago) by jfenwick
File MIME type: text/plain
File size: 34889 byte(s)
More warning removal.

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 /**************************************************************/
16
17 /* Finley: generates rectangular meshes */
18
19 /* Generates a numElements[0] x numElements[1] x numElements[2] mesh with second order elements (Hex20) in the brick */
20 /* [0,Length[0]] x [0,Length[1]] x [0,Length[2]]. order is the desired accuracy of the */
21 /* integration scheme. */
22
23
24 /**************************************************************/
25
26 #include "RectangularMesh.h"
27
28 Finley_Mesh* Finley_RectangularMesh_Hex20(dim_t* numElements,
29 double* Length,
30 bool_t* periodic,
31 index_t order,
32 index_t reduced_order,
33 bool_t useElementsOnFace,
34 bool_t useFullElementOrder,
35 bool_t optimize)
36 {
37 #define N_PER_E 2
38 #define DIM 3
39 dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,Nstride0,Nstride1,Nstride2, local_NE0, local_NE1, local_NE2;
40 dim_t totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements, local_N0, local_N1, local_N2, NN;
41 index_t node0, myRank, e_offset0, e_offset1, e_offset2, offset0, offset1, offset2, global_i0, global_i1, global_i2;
42 Finley_Mesh* out;
43 Paso_MPIInfo *mpi_info = NULL;
44 char name[50];
45 #ifdef Finley_TRACE
46 double time0=Finley_timer();
47 #endif
48
49 /* get MPI information */
50 mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
51 if (! Finley_noError()) {
52 return NULL;
53 }
54 myRank=mpi_info->rank;
55
56 /* set up the global dimensions of the mesh */
57
58 NE0=MAX(1,numElements[0]);
59 NE1=MAX(1,numElements[1]);
60 NE2=MAX(1,numElements[2]);
61 N0=N_PER_E*NE0+1;
62 N1=N_PER_E*NE1+1;
63 N2=N_PER_E*NE2+1;
64
65 /* allocate mesh: */
66 sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
67 out=Finley_Mesh_alloc(name,DIM,order, reduced_order, mpi_info);
68 if (! Finley_noError()) {
69 Paso_MPIInfo_free( mpi_info );
70 return NULL;
71 }
72
73 if (useFullElementOrder) {
74 /* Finley_setError(SYSTEM_ERROR,"full element order for Hex elements is not supported yet."); */
75 Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Hex27,
76 out->order,
77 out->reduced_order,
78 mpi_info));
79 if (useElementsOnFace) {
80 Finley_setError(SYSTEM_ERROR,"rich elements for Hex27 elements is not supported yet.");
81 } else {
82 Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Rec9,
83 out->order,
84 out->reduced_order,
85 mpi_info));
86 Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Rec9_Contact,
87 out->order,
88 out->reduced_order,
89 mpi_info));
90 }
91
92 } else {
93 Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Hex20,out->order,out->reduced_order,mpi_info));
94 if (useElementsOnFace) {
95 Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Hex20Face,
96 out->order,
97 out->reduced_order,
98 mpi_info));
99 Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Hex20Face_Contact,
100 out->order,
101 out->reduced_order,
102 mpi_info));
103 } else {
104 Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Rec8,
105 out->order,
106 out->reduced_order,
107 mpi_info));
108 Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Rec8_Contact,
109 out->order,
110 out->reduced_order,
111 mpi_info));
112 }
113 }
114 Finley_Mesh_setPoints(out,Finley_ElementFile_alloc(Point1,
115 out->order,
116 out->reduced_order,
117 mpi_info));
118 if (! Finley_noError()) {
119 Paso_MPIInfo_free( mpi_info );
120 Finley_Mesh_free(out);
121 return NULL;
122 }
123
124 /* work out the largest dimension */
125 if (N2==MAX3(N0,N1,N2)) {
126 Nstride0=1;
127 Nstride1=N0;
128 Nstride2=N0*N1;
129 local_NE0=NE0;
130 e_offset0=0;
131 local_NE1=NE1;
132 e_offset1=0;
133 Paso_MPIInfo_Split(mpi_info,NE2,&local_NE2,&e_offset2);
134 } else if (N1==MAX3(N0,N1,N2)) {
135 Nstride0=N2;
136 Nstride1=N0*N2;
137 Nstride2=1;
138 local_NE0=NE0;
139 e_offset0=0;
140 Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1);
141 local_NE2=NE2;
142 e_offset2=0;
143 } else {
144 Nstride0=N1*N2;
145 Nstride1=1;
146 Nstride2=N1;
147 Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0);
148 local_NE1=NE1;
149 e_offset1=0;
150 local_NE2=NE2;
151 e_offset2=0;
152 }
153 offset0=e_offset0*N_PER_E;
154 offset1=e_offset1*N_PER_E;
155 offset2=e_offset2*N_PER_E;
156 local_N0=local_NE0>0 ? local_NE0*N_PER_E+1 : 0;
157 local_N1=local_NE1>0 ? local_NE1*N_PER_E+1 : 0;
158 local_N2=local_NE0>0 ? local_NE2*N_PER_E+1 : 0;
159
160 /* get the number of surface elements */
161
162 NFaceElements=0;
163 if (!periodic[2] && (local_NE2>0) ) {
164 NDOF2=N2;
165 if (offset2==0) NFaceElements+=local_NE1*local_NE0;
166 if (local_NE2+e_offset2 == NE2) NFaceElements+=local_NE1*local_NE0;
167 } else {
168 NDOF2=N2-1;
169 }
170
171 if (!periodic[0] && (local_NE0>0) ) {
172 NDOF0=N0;
173 if (e_offset0 == 0) NFaceElements+=local_NE1*local_NE2;
174 if (local_NE0+e_offset0 == NE0) NFaceElements+=local_NE1*local_NE2;
175 } else {
176 NDOF0=N0-1;
177 }
178 if (!periodic[1] && (local_NE1>0) ) {
179 NDOF1=N1;
180 if (e_offset1 == 0) NFaceElements+=local_NE0*local_NE2;
181 if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0*local_NE2;
182 } else {
183 NDOF1=N1-1;
184 }
185 /* allocate tables: */
186
187 Finley_NodeFile_allocTable(out->Nodes,local_N0*local_N1*local_N2);
188 Finley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1*local_NE2);
189 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
190
191 if (Finley_noError()) {
192 /* create nodes */
193
194 #pragma omp parallel for private(i0,i1,i2,k,global_i0,global_i1,global_i2)
195 for (i2=0;i2<local_N2;i2++) {
196 for (i1=0;i1<local_N1;i1++) {
197 for (i0=0;i0<local_N0;i0++) {
198 k=i0+local_N0*i1+local_N0*local_N1*i2;
199 global_i0=i0+offset0;
200 global_i1=i1+offset1;
201 global_i2=i2+offset2;
202 out->Nodes->Coordinates[INDEX2(0,k,DIM)]=DBLE(global_i0)/DBLE(N0-1)*Length[0];
203 out->Nodes->Coordinates[INDEX2(1,k,DIM)]=DBLE(global_i1)/DBLE(N1-1)*Length[1];
204 out->Nodes->Coordinates[INDEX2(2,k,DIM)]=DBLE(global_i2)/DBLE(N2-1)*Length[2];
205 out->Nodes->Id[k]=Nstride0*global_i0+Nstride1*global_i1+Nstride2*global_i2;
206 out->Nodes->Tag[k]=0;
207 out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
208 +Nstride1*(global_i1%NDOF1)
209 +Nstride2*(global_i2%NDOF2);
210 }
211 }
212 }
213 /* set the elements: */
214 NN=out->Elements->numNodes;
215 #pragma omp parallel for private(i0,i1,i2,k,node0)
216 for (i2=0;i2<local_NE2;i2++) {
217 for (i1=0;i1<local_NE1;i1++) {
218 for (i0=0;i0<local_NE0;i0++) {
219
220 k=i0+local_NE0*i1+local_NE0*local_NE1*i2;
221 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
222
223 out->Elements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+NE0*NE1*(i2+e_offset2);
224 out->Elements->Tag[k]=0;
225 out->Elements->Owner[k]=myRank;
226
227 out->Elements->Nodes[INDEX2(0,k,NN)] =node0 ;
228 out->Elements->Nodes[INDEX2(1,k,NN)] =node0+ 2*Nstride0;
229 out->Elements->Nodes[INDEX2(2,k,NN)] =node0+ 2*Nstride1+2*Nstride0;
230 out->Elements->Nodes[INDEX2(3,k,NN)] =node0+ 2*Nstride1;
231 out->Elements->Nodes[INDEX2(4,k,NN)] =node0+2*Nstride2 ;
232 out->Elements->Nodes[INDEX2(5,k,NN)] =node0+2*Nstride2 +2*Nstride0;
233 out->Elements->Nodes[INDEX2(6,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
234 out->Elements->Nodes[INDEX2(7,k,NN)] =node0+2*Nstride2+2*Nstride1 ;
235 out->Elements->Nodes[INDEX2(8,k,NN)] =node0+ 1*Nstride0;
236 out->Elements->Nodes[INDEX2(9,k,NN)] =node0+ 1*Nstride1+2*Nstride0;
237 out->Elements->Nodes[INDEX2(10,k,NN)]=node0+ 2*Nstride1+1*Nstride0;
238 out->Elements->Nodes[INDEX2(11,k,NN)]=node0+ 1*Nstride1 ;
239 out->Elements->Nodes[INDEX2(12,k,NN)]=node0+1*Nstride2 ;
240 out->Elements->Nodes[INDEX2(13,k,NN)]=node0+1*Nstride2 +2*Nstride0;
241 out->Elements->Nodes[INDEX2(14,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
242 out->Elements->Nodes[INDEX2(15,k,NN)]=node0+1*Nstride2+2*Nstride1 ;
243 out->Elements->Nodes[INDEX2(16,k,NN)]=node0+2*Nstride2 +1*Nstride0;
244 out->Elements->Nodes[INDEX2(17,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
245 out->Elements->Nodes[INDEX2(18,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
246 out->Elements->Nodes[INDEX2(19,k,NN)]=node0+2*Nstride2+1*Nstride1 ;
247 if (useFullElementOrder) {
248 out->Elements->Nodes[INDEX2(20,k,NN)]=node0+ 1*Nstride1+1*Nstride0;
249 out->Elements->Nodes[INDEX2(21,k,NN)]=node0+1*Nstride2 +1*Nstride0;
250 out->Elements->Nodes[INDEX2(22,k,NN)]=node0+1*Nstride2+1*Nstride1+2*Nstride0;
251 out->Elements->Nodes[INDEX2(23,k,NN)]=node0+1*Nstride2+2*Nstride1+1*Nstride0;
252 out->Elements->Nodes[INDEX2(24,k,NN)]=node0+1*Nstride2+1*Nstride1 ;
253 out->Elements->Nodes[INDEX2(25,k,NN)]=node0+2*Nstride2+1*Nstride1+1*Nstride0;
254 out->Elements->Nodes[INDEX2(26,k,NN)]=node0+1*Nstride2+1*Nstride1+1*Nstride0;
255 }
256 }
257 }
258 }
259 /* face elements */
260 NN=out->FaceElements->numNodes;
261 totalNECount=NE0*NE1*NE2;
262 faceNECount=0;
263 /* these are the quadrilateral elements on boundary 1 (x3=0): */
264 if (!periodic[2] && (local_NE2>0)) {
265 /* ** elements on boundary 100 (x3=0): */
266 if (offset2==0) {
267 #pragma omp parallel for private(i0,i1,k,node0)
268 for (i1=0;i1<local_NE1;i1++) {
269 for (i0=0;i0<local_NE0;i0++) {
270
271 k=i0+local_NE0*i1+faceNECount;
272 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
273
274 out->FaceElements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+totalNECount;
275 out->FaceElements->Tag[k]=100;
276 out->FaceElements->Owner[k]=myRank;
277
278 if (useElementsOnFace) {
279 out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0 ;
280 out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0 +2*Nstride1 ;
281 out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0 +2*Nstride1+2*Nstride0;
282 out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+ 2*Nstride0 ;
283 out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+2*Nstride2 ;
284 out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+2*Nstride2+2*Nstride1 ;
285 out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
286 out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+2*Nstride2 +2*Nstride0;
287 out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+ 1*Nstride1 ;
288 out->FaceElements->Nodes[INDEX2(9,k,NN)] =node0+ 2*Nstride1+1*Nstride0;
289 out->FaceElements->Nodes[INDEX2(10,k,NN)]=node0+ 1*Nstride1+2*Nstride0;
290 out->FaceElements->Nodes[INDEX2(11,k,NN)]=node0+ 1*Nstride0;
291 out->FaceElements->Nodes[INDEX2(12,k,NN)]=node0+1*Nstride2 ;
292 out->FaceElements->Nodes[INDEX2(13,k,NN)]=node0+1*Nstride2+2*Nstride1 ;
293 out->FaceElements->Nodes[INDEX2(14,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
294 out->FaceElements->Nodes[INDEX2(15,k,NN)]=node0+1*Nstride2 +2*Nstride0;
295 out->FaceElements->Nodes[INDEX2(16,k,NN)]=node0+2*Nstride2+1*Nstride1;
296 out->FaceElements->Nodes[INDEX2(17,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
297 out->FaceElements->Nodes[INDEX2(18,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
298 out->FaceElements->Nodes[INDEX2(19,k,NN)]=node0+2*Nstride2 +1*Nstride0;
299 } else {
300 out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0 ;
301 out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+ 2*Nstride1 ;
302 out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+ 2*Nstride1+2*Nstride0;
303 out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+ 2*Nstride0;
304 out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+ 1*Nstride1 ;
305 out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+ 2*Nstride1+1*Nstride0;
306 out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+ 1*Nstride1+2*Nstride0;
307 out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+ 1*Nstride0;
308 if (useFullElementOrder){
309 out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+ 1*Nstride1+1*Nstride0;
310 }
311 }
312 }
313 }
314 faceNECount+=local_NE1*local_NE0;
315 }
316 totalNECount+=NE1*NE0;
317 /* ** elements on boundary 200 (x3=1) */
318 if (local_NE2+e_offset2 == NE2) {
319 #pragma omp parallel for private(i0,i1,k,node0)
320 for (i1=0;i1<local_NE1;i1++) {
321 for (i0=0;i0<local_NE0;i0++) {
322
323 k=i0+local_NE0*i1+faceNECount;
324 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(NE2-1);
325
326 out->FaceElements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+totalNECount;
327 out->FaceElements->Tag[k]=200;
328 out->FaceElements->Owner[k]=myRank;
329 if (useElementsOnFace) {
330 out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0+2*Nstride2 ;
331 out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+2*Nstride2+ 2*Nstride0;
332 out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
333 out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+2*Nstride2+2*Nstride1 ;
334
335 out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0 ;
336 out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+2*Nstride0 ;
337 out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+ 2*Nstride1+2*Nstride0;
338 out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+ 2*Nstride1;
339
340 out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+2*Nstride2+ 1*Nstride0;
341 out->FaceElements->Nodes[INDEX2(9,k,NN)] =node0+2*Nstride2+1*Nstride1+2*Nstride0;
342 out->FaceElements->Nodes[INDEX2(10,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
343 out->FaceElements->Nodes[INDEX2(11,k,NN)]=node0+2*Nstride2+1*Nstride1 ;
344
345 out->FaceElements->Nodes[INDEX2(12,k,NN)]=node0+1*Nstride2;
346 out->FaceElements->Nodes[INDEX2(13,k,NN)]=node0+1*Nstride2 +2*Nstride0;
347 out->FaceElements->Nodes[INDEX2(14,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
348 out->FaceElements->Nodes[INDEX2(15,k,NN)]=node0+1*Nstride2+2*Nstride1 ;
349
350 out->FaceElements->Nodes[INDEX2(16,k,NN)]=node0+ 1*Nstride0;
351 out->FaceElements->Nodes[INDEX2(17,k,NN)]=node0+ 1*Nstride1+2*Nstride0;
352 out->FaceElements->Nodes[INDEX2(18,k,NN)]=node0+ 2*Nstride1+1*Nstride0;
353 out->FaceElements->Nodes[INDEX2(19,k,NN)]=node0+ 1*Nstride1 ;
354
355 } else {
356 out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0+2*Nstride2 ;
357 out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+2*Nstride2 +2*Nstride0;
358 out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
359 out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+2*Nstride2+2*Nstride1 ;
360 out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+2*Nstride2 +1*Nstride0;
361 out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+2*Nstride2+1*Nstride1+2*Nstride0;
362 out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+2*Nstride2+2*Nstride1+1*Nstride0;
363 out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+2*Nstride2+1*Nstride1 ;
364 if (useFullElementOrder){
365 out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+2*Nstride2+1*Nstride1+1*Nstride0;
366 }
367 }
368 }
369 }
370 faceNECount+=local_NE1*local_NE0;
371 }
372 totalNECount+=NE1*NE0;
373 }
374 if (!periodic[0] && (local_NE0>0)) {
375 /* ** elements on boundary 001 (x1=0): */
376
377 if (e_offset0 == 0) {
378 #pragma omp parallel for private(i1,i2,k,node0)
379 for (i2=0;i2<local_NE2;i2++) {
380 for (i1=0;i1<local_NE1;i1++) {
381
382 k=i1+local_NE1*i2+faceNECount;
383 node0=Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
384 out->FaceElements->Id[k]=(i1+e_offset1)+NE1*(i2+e_offset2)+totalNECount;
385 out->FaceElements->Tag[k]=1;
386 out->FaceElements->Owner[k]=myRank;
387
388 if (useElementsOnFace) {
389 out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0 ;
390 out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+2*Nstride2 ;
391 out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+2*Nstride2+2*Nstride1 ;
392 out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+2*Nstride1 ;
393
394 out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+2*Nstride0 ;
395 out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+2*Nstride2+2*Nstride0 ;
396 out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
397 out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+2*Nstride1+2*Nstride0 ;
398
399 out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+1*Nstride2 ;
400 out->FaceElements->Nodes[INDEX2(9,k,NN)] =node0+2*Nstride2+1*Nstride1 ;
401 out->FaceElements->Nodes[INDEX2(10,k,NN)]=node0+1*Nstride2+2*Nstride1 ;
402 out->FaceElements->Nodes[INDEX2(11,k,NN)]=node0+ 1*Nstride1 ;
403
404 out->FaceElements->Nodes[INDEX2(12,k,NN)]=node0+ 1*Nstride0;
405 out->FaceElements->Nodes[INDEX2(13,k,NN)]=node0+2*Nstride2 +1*Nstride0;
406 out->FaceElements->Nodes[INDEX2(14,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
407 out->FaceElements->Nodes[INDEX2(15,k,NN)]=node0+2*Nstride1+ 1*Nstride0;
408
409 out->FaceElements->Nodes[INDEX2(16,k,NN)]=node0+1*Nstride2+ 2*Nstride0;
410 out->FaceElements->Nodes[INDEX2(17,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
411 out->FaceElements->Nodes[INDEX2(18,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
412 out->FaceElements->Nodes[INDEX2(19,k,NN)]=node0+1*Nstride1+ 2*Nstride0;
413
414 } else {
415 out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0 ;
416 out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+2*Nstride2 ;
417 out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+2*Nstride2+2*Nstride1 ;
418 out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+ 2*Nstride1 ;
419 out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+1*Nstride2 ;
420 out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+2*Nstride2+1*Nstride1 ;
421 out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+1*Nstride2+2*Nstride1 ;
422 out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+ 1*Nstride1 ;
423 if (useFullElementOrder){
424 out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+1*Nstride2+1*Nstride1 ;
425 }
426 }
427 }
428 }
429 faceNECount+=local_NE1*local_NE2;
430 }
431 totalNECount+=NE1*NE2;
432
433 /* ** elements on boundary 002 (x1=1): */
434 if (local_NE0+e_offset0 == NE0) {
435 #pragma omp parallel for private(i1,i2,k,node0)
436 for (i2=0;i2<local_NE2;i2++) {
437 for (i1=0;i1<local_NE1;i1++) {
438 k=i1+local_NE1*i2+faceNECount;
439 node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
440 out->FaceElements->Id[k]=(i1+e_offset1)+NE1*(i2+e_offset2)+totalNECount;
441 out->FaceElements->Tag[k]=2;
442 out->FaceElements->Owner[k]=myRank;
443
444 if (useElementsOnFace) {
445 out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+ 2*Nstride0;
446 out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+ 2*Nstride1+2*Nstride0;
447 out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
448 out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride2+ 2*Nstride0;
449
450 out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0 ;
451 out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+ 2*Nstride1 ;
452 out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride2+2*Nstride1 ;
453 out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+2*Nstride2 ;
454
455 out->FaceElements->Nodes[INDEX2(8,k,NN)]=node0+ 1*Nstride1+2*Nstride0;
456 out->FaceElements->Nodes[INDEX2(9,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
457 out->FaceElements->Nodes[INDEX2(10,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
458 out->FaceElements->Nodes[INDEX2(11,k,NN)]=node0+1*Nstride2+ 2*Nstride0;
459
460 out->FaceElements->Nodes[INDEX2(12,k,NN)]=node0+ 1*Nstride0;
461 out->FaceElements->Nodes[INDEX2(13,k,NN)]=node0+ 2*Nstride1+1*Nstride0;
462 out->FaceElements->Nodes[INDEX2(14,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
463 out->FaceElements->Nodes[INDEX2(15,k,NN)]=node0+2*Nstride2+ 1*Nstride0;
464
465 out->FaceElements->Nodes[INDEX2(16,k,NN)]=node0+ 1*Nstride1 ;
466 out->FaceElements->Nodes[INDEX2(17,k,NN)]=node0+1*Nstride2+2*Nstride1 ;
467 out->FaceElements->Nodes[INDEX2(18,k,NN)]=node0+2*Nstride2+1*Nstride1 ;
468 out->FaceElements->Nodes[INDEX2(19,k,NN)]=node0+1*Nstride2 ;
469
470 } else {
471 out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0 +2*Nstride0;
472 out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+ 2*Nstride1+2*Nstride0;
473 out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
474 out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride2+ 2*Nstride0;
475 out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+ 1*Nstride1+2*Nstride0;
476 out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
477 out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
478 out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+1*Nstride2 +2*Nstride0;
479 if (useFullElementOrder){
480 out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+1*Nstride2+1*Nstride1+2*Nstride0;
481 }
482 }
483
484 }
485 }
486 faceNECount+=local_NE1*local_NE2;
487 }
488 totalNECount+=NE1*NE2;
489 }
490 if (!periodic[1] && (local_NE1>0)) {
491 /* ** elements on boundary 010 (x2=0): */
492 if (e_offset1 == 0) {
493 #pragma omp parallel for private(i0,i2,k,node0)
494 for (i2=0;i2<local_NE2;i2++) {
495 for (i0=0;i0<local_NE0;i0++) {
496 k=i0+local_NE0*i2+faceNECount;
497 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride2*N_PER_E*(i2+e_offset2);
498
499 out->FaceElements->Id[k]=(i2+e_offset2)+NE2*(e_offset0+i0)+totalNECount;
500 out->FaceElements->Tag[k]=10;
501 out->FaceElements->Owner[k]=myRank;
502 if (useElementsOnFace) {
503 out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0 ;
504 out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+ 2*Nstride0;
505 out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride2 +2*Nstride0;
506 out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride2 ;
507
508 out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+ 2*Nstride1 ;
509 out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+2*Nstride1+ 2*Nstride0;
510 out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
511 out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+2*Nstride2+2*Nstride1 ;
512
513 out->FaceElements->Nodes[INDEX2(8,k,NN)]=node0+ 1*Nstride0;
514 out->FaceElements->Nodes[INDEX2(9,k,NN)]=node0+1*Nstride2+ 2*Nstride0;
515 out->FaceElements->Nodes[INDEX2(10,k,NN)]=node0+2*Nstride2+ 1*Nstride0;
516 out->FaceElements->Nodes[INDEX2(11,k,NN)]=node0+1*Nstride2 ;
517
518 out->FaceElements->Nodes[INDEX2(12,k,NN)]=node0+ 1*Nstride1 ;
519 out->FaceElements->Nodes[INDEX2(13,k,NN)]=node0+ 1*Nstride1+2*Nstride0;
520 out->FaceElements->Nodes[INDEX2(14,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
521 out->FaceElements->Nodes[INDEX2(15,k,NN)]=node0+2*Nstride2+1*Nstride1 ;
522
523 out->FaceElements->Nodes[INDEX2(16,k,NN)]=node0+ 2*Nstride1+1*Nstride0;
524 out->FaceElements->Nodes[INDEX2(17,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
525 out->FaceElements->Nodes[INDEX2(18,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
526 out->FaceElements->Nodes[INDEX2(19,k,NN)]=node0+1*Nstride2+2*Nstride1 ;
527
528 } else {
529 out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0 ;
530 out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+ 2*Nstride0;
531 out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride2+ 2*Nstride0;
532 out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride2 ;
533 out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+ 1*Nstride0;
534 out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+1*Nstride2+ 2*Nstride0;
535 out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride2+ 1*Nstride0;
536 out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+1*Nstride2 ;
537 if (useFullElementOrder){
538 out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+1*Nstride2+ 1*Nstride0;
539 }
540 }
541 }
542 }
543 faceNECount+=local_NE0*local_NE2;
544 }
545 totalNECount+=NE0*NE2;
546 /* ** elements on boundary 020 (x2=1): */
547 if (local_NE1+e_offset1 == NE1) {
548 #pragma omp parallel for private(i0,i2,k,node0)
549 for (i2=0;i2<local_NE2;i2++) {
550 for (i0=0;i0<local_NE0;i0++) {
551 k=i0+local_NE0*i2+faceNECount;
552 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1)+Nstride2*N_PER_E*(i2+e_offset2);
553
554 out->FaceElements->Id[k]=(i2+e_offset2)+NE2*(i0+e_offset0)+totalNECount;
555 out->FaceElements->Tag[k]=20;
556 out->FaceElements->Owner[k]=myRank;
557
558 if (useElementsOnFace) {
559 out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+ 2*Nstride1 ;
560 out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride2+2*Nstride1 ;
561 out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
562 out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride1+2*Nstride0 ;
563
564 out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0 ;
565 out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+2*Nstride2 ;
566 out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride2+ 2*Nstride0;
567 out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+ 2*Nstride0;
568
569 out->FaceElements->Nodes[INDEX2(8,k,NN)]=node0+1*Nstride2+2*Nstride1 ;
570 out->FaceElements->Nodes[INDEX2(9,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
571 out->FaceElements->Nodes[INDEX2(10,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
572 out->FaceElements->Nodes[INDEX2(11,k,NN)]=node0+ 2*Nstride1+1*Nstride0;
573
574 out->FaceElements->Nodes[INDEX2(12,k,NN)]=node0+ 1*Nstride1 ;
575 out->FaceElements->Nodes[INDEX2(13,k,NN)]=node0+2*Nstride2+1*Nstride1 ;
576 out->FaceElements->Nodes[INDEX2(14,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
577 out->FaceElements->Nodes[INDEX2(15,k,NN)]=node0+ 1*Nstride1+2*Nstride0;
578
579 out->FaceElements->Nodes[INDEX2(16,k,NN)]=node0+1*Nstride2 ;
580 out->FaceElements->Nodes[INDEX2(17,k,NN)]=node0+2*Nstride2 +1*Nstride0;
581 out->FaceElements->Nodes[INDEX2(18,k,NN)]=node0+1*Nstride2 +2*Nstride0;
582 out->FaceElements->Nodes[INDEX2(19,k,NN)]=node0+ 1*Nstride0;
583 } else {
584 out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+ 2*Nstride1 ;
585 out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride2+2*Nstride1 ;
586 out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
587 out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+ 2*Nstride1+2*Nstride0;
588 out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+1*Nstride2+2*Nstride1 ;
589 out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
590 out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
591 out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+ 2*Nstride1+1*Nstride0;
592 if (useFullElementOrder){
593 out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+1*Nstride2+2*Nstride1+1*Nstride0;
594 }
595 }
596 }
597 }
598 faceNECount+=local_NE0*local_NE2;
599 }
600 totalNECount+=NE0*NE2;
601 }
602 /* add tag names */
603 Finley_Mesh_addTagMap(out,"top", 200);
604 Finley_Mesh_addTagMap(out,"bottom", 100);
605 Finley_Mesh_addTagMap(out,"left", 1);
606 Finley_Mesh_addTagMap(out,"right", 2);
607 Finley_Mesh_addTagMap(out,"front", 10);
608 Finley_Mesh_addTagMap(out,"back", 20);
609
610 /* prepare mesh for further calculatuions:*/
611 if (Finley_noError()) {
612 Finley_Mesh_resolveNodeIds(out);
613 }
614 if (Finley_noError()) {
615 Finley_Mesh_prepare(out, optimize);
616 }
617 }
618
619 if (!Finley_noError()) {
620 Finley_Mesh_free(out);
621 }
622 /* free up memory */
623 Paso_MPIInfo_free( mpi_info );
624 #ifdef Finley_TRACE
625 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
626 #endif
627
628 return out;
629 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26