/[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 2722 - (show annotations)
Fri Oct 16 06:45:01 2009 UTC (10 years, 5 months ago) by gross
File MIME type: text/plain
File size: 35702 byte(s)
First steps towards macro elements. at the monent a  macro element is identical to its counterpart of 2nd order.




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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26