/[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 1733 - (show annotations)
Thu Aug 28 01:46:07 2008 UTC (11 years, 3 months ago) by gross
File MIME type: text/plain
File size: 34952 byte(s)
bug fix: if the numer of elements to be generated suppose to be zero
still face elements and nodes are generated. this is fixed now.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26