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

Contents of /trunk-mpi-branch/finley/src/Mesh_hex20.c

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26