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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26