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

Annotation of /trunk/finley/src/Mesh_rec8.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2856 - (hide annotations)
Mon Jan 18 04:14:37 2010 UTC (9 years, 8 months ago) by gross
File MIME type: text/plain
File size: 14396 byte(s)
FunctionSpaces provide now some information about their approximation order.
1 jgs 150
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2548 * Copyright (c) 2003-2009 by University of Queensland
5 ksteube 1811 * 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 ksteube 1312
14 ksteube 1811
15 jgs 82 /**************************************************************/
16    
17     /* Finley: generates rectangular meshes */
18    
19     /* Generates a numElements[0] x numElements[1] mesh with second order elements (Rec8) in the rectangle */
20     /* [0,Length[0]] x [0,Length[1]]. order is the desired accuracy of the integration scheme. */
21    
22     /**************************************************************/
23    
24     #include "RectangularMesh.h"
25    
26 bcumming 782
27 ksteube 1312 Finley_Mesh* Finley_RectangularMesh_Rec8(dim_t* numElements,
28     double* Length,
29     bool_t* periodic,
30     index_t order,
31     index_t reduced_order,
32     bool_t useElementsOnFace,
33     bool_t useFullElementOrder,
34 gross 2722 bool_t useMacroElements,
35 ksteube 1312 bool_t optimize)
36 bcumming 782 {
37 ksteube 1312 #define N_PER_E 2
38     #define DIM 2
39 gross 2748 dim_t N0,N1,NE0,NE1,i0,i1,k,Nstride0=0,Nstride1=0;
40     dim_t totalNECount,faceNECount,NDOF0=0,NDOF1=0,NFaceElements,NN, local_NE0, local_NE1, local_N0=0, local_N1=0;
41     index_t e_offset1, e_offset0, offset0=0, offset1=0, global_i0, global_i1;
42     Finley_ReferenceElementSet *refPoints=NULL, *refContactElements=NULL, *refFaceElements=NULL, *refElements=NULL;
43 ksteube 1312 index_t node0, myRank;
44     Finley_Mesh* out;
45     Paso_MPIInfo *mpi_info = NULL;
46     char name[50];
47 gross 2722 bool_t generateAllNodes= useFullElementOrder || useMacroElements;
48 jfenwick 1986 #ifdef Finley_TRACE
49 ksteube 1312 double time0=Finley_timer();
50 jfenwick 1986 #endif
51 bcumming 782
52 ksteube 1312 /* get MPI information */
53     mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
54     if (! Finley_noError()) {
55     return NULL;
56 bcumming 782 }
57 ksteube 1312 myRank=mpi_info->rank;
58 bcumming 782
59 ksteube 1312 /* set up the global dimensions of the mesh */
60 bcumming 782
61 jgs 82 NE0=MAX(1,numElements[0]);
62     NE1=MAX(1,numElements[1]);
63 ksteube 1312 N0=N_PER_E*NE0+1;
64     N1=N_PER_E*NE1+1;
65 jgs 82
66 ksteube 1312 /* allocate mesh: */
67 jgs 82 sprintf(name,"Rectangular %d x %d mesh",N0,N1);
68 gross 2856 out=Finley_Mesh_alloc(name,DIM, mpi_info);
69 ksteube 1312 if (! Finley_noError()) {
70     Paso_MPIInfo_free( mpi_info );
71 jgs 82 return NULL;
72     }
73 gross 2722 if (generateAllNodes) {
74 gross 2748 /* Finley_setError(SYSTEM_ERROR,"full element order for Hex elements is not supported yet."); */
75 gross 2722 if (useMacroElements) {
76 gross 2856 refElements= Finley_ReferenceElementSet_alloc(Rec9Macro,order,reduced_order);
77 gross 2722 } else {
78 gross 2856 refElements=Finley_ReferenceElementSet_alloc(Rec9, order,reduced_order);
79 gross 2722 }
80 ksteube 1312 if (useElementsOnFace) {
81     Finley_setError(SYSTEM_ERROR,"rich elements for Rec9 elements is not supported yet.");
82     } else {
83 gross 2748 if (useMacroElements) {
84 gross 2856 refFaceElements=Finley_ReferenceElementSet_alloc(Line3Macro, order, reduced_order);
85 gross 2722 } else {
86 gross 2856 refFaceElements=Finley_ReferenceElementSet_alloc(Line3, order, reduced_order);
87 gross 2722 }
88 gross 2856 refContactElements=Finley_ReferenceElementSet_alloc(Line3_Contact, order, reduced_order);
89 jgs 82 }
90 gross 2748
91 ksteube 1312 } else {
92 gross 2856 refElements= Finley_ReferenceElementSet_alloc(Rec8,order,reduced_order);
93 ksteube 1312 if (useElementsOnFace) {
94 gross 2856 refFaceElements= Finley_ReferenceElementSet_alloc(Rec8Face ,order,reduced_order);
95     refContactElements=Finley_ReferenceElementSet_alloc(Rec8Face_Contact, order, reduced_order);
96 gross 2748
97 ksteube 1312 } else {
98 gross 2856 refFaceElements= Finley_ReferenceElementSet_alloc(Line3 ,order,reduced_order);
99     refContactElements=Finley_ReferenceElementSet_alloc(Line3_Contact, order, reduced_order);
100 gross 2748
101 jgs 82 }
102     }
103 gross 2856 refPoints=Finley_ReferenceElementSet_alloc(Point1, order, reduced_order);
104 jgs 82
105 bcumming 782
106 gross 2748 if ( Finley_noError()) {
107    
108     Finley_Mesh_setPoints(out,Finley_ElementFile_alloc(refPoints, mpi_info));
109     Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(refContactElements, mpi_info));
110     Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(refFaceElements, mpi_info));
111     Finley_Mesh_setElements(out,Finley_ElementFile_alloc(refElements, mpi_info));
112 bcumming 782
113 gross 2748 /* work out the largest dimension */
114     if (N1==MAX(N0,N1)) {
115     Nstride0=1;
116     Nstride1=N0;
117     local_NE0=NE0;
118     e_offset0=0;
119     Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1);
120     } else {
121     Nstride0=N1;
122     Nstride1=1;
123     Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0);
124     local_NE1=NE1;
125     e_offset1=0;
126     }
127     offset0=e_offset0*N_PER_E;
128     offset1=e_offset1*N_PER_E;
129     local_N0=local_NE0>0 ? local_NE0*N_PER_E+1 : 0;
130     local_N1=local_NE1>0 ? local_NE1*N_PER_E+1 : 0;
131 bcumming 782
132 gross 2748 /* get the number of surface elements */
133 bcumming 782
134 gross 2748 NFaceElements=0;
135     if (!periodic[0] && (local_NE0>0)) {
136     NDOF0=N0;
137     if (e_offset0 == 0) NFaceElements+=local_NE1;
138     if (local_NE0+e_offset0 == NE0) NFaceElements+=local_NE1;
139     } else {
140     NDOF0=N0-1;
141     }
142     if (!periodic[1] && (local_NE1>0)) {
143     NDOF1=N1;
144     if (e_offset1 == 0) NFaceElements+=local_NE0;
145     if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0;
146     } else {
147     NDOF1=N1-1;
148     }
149    
150     /* allocate tables: */
151 bcumming 782
152 gross 2748 Finley_NodeFile_allocTable(out->Nodes,local_N0*local_N1);
153     Finley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1);
154     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
155     }
156    
157 ksteube 1312 if (Finley_noError()) {
158     /* create nodes */
159    
160     #pragma omp parallel for private(i0,i1,k,global_i0,global_i1)
161     for (i1=0;i1<local_N1;i1++) {
162     for (i0=0;i0<local_N0;i0++) {
163     k=i0+local_N0*i1;
164     global_i0=i0+offset0;
165     global_i1=i1+offset1;
166     out->Nodes->Coordinates[INDEX2(0,k,DIM)]=DBLE(global_i0)/DBLE(N0-1)*Length[0];
167     out->Nodes->Coordinates[INDEX2(1,k,DIM)]=DBLE(global_i1)/DBLE(N1-1)*Length[1];
168     out->Nodes->Id[k]=Nstride0*global_i0+Nstride1*global_i1;
169     out->Nodes->Tag[k]=0;
170     out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
171     +Nstride1*(global_i1%NDOF1);
172     }
173     }
174     /* set the elements: */
175     NN=out->Elements->numNodes;
176     #pragma omp parallel for private(i0,i1,k,node0)
177     for (i1=0;i1<local_NE1;i1++) {
178     for (i0=0;i0<local_NE0;i0++) {
179    
180     k=i0+local_NE0*i1;
181     node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
182 bcumming 782
183 ksteube 1312 out->Elements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1);
184     out->Elements->Tag[k]=0;
185     out->Elements->Owner[k]=myRank;
186 bcumming 782
187 ksteube 1312 out->Elements->Nodes[INDEX2(0,k,NN)]=node0;
188     out->Elements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride0;
189     out->Elements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride1+2*Nstride0;
190     out->Elements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride1;
191     out->Elements->Nodes[INDEX2(4,k,NN)]=node0+1*Nstride0;
192     out->Elements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1+2*Nstride0;
193     out->Elements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride1+1*Nstride0;
194     out->Elements->Nodes[INDEX2(7,k,NN)]=node0+Nstride1;
195 gross 2722 if (generateAllNodes) {
196 ksteube 1312 out->Elements->Nodes[INDEX2(8,k,NN)]=node0+1*Nstride1+1*Nstride0;
197     }
198     }
199 bcumming 782 }
200 ksteube 1312 /* face elements */
201     NN=out->FaceElements->numNodes;
202     totalNECount=NE0*NE1;
203     faceNECount=0;
204 gross 1818 if (!periodic[0] && (local_NE0>0)) {
205 ksteube 1312 /* ** elements on boundary 001 (x1=0): */
206    
207     if (e_offset0 == 0) {
208     #pragma omp parallel for private(i1,k,node0)
209     for (i1=0;i1<local_NE1;i1++) {
210    
211     k=i1+faceNECount;
212     node0=Nstride1*N_PER_E*(i1+e_offset1);
213 bcumming 782
214 ksteube 1312 out->FaceElements->Id[k]=i1+e_offset1+totalNECount;
215     out->FaceElements->Tag[k]=1;
216     out->FaceElements->Owner[k]=myRank;
217     if (useElementsOnFace) {
218     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride1;
219     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
220     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride0;
221     out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride1+2*Nstride0;
222     out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+Nstride1;
223     out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+1*Nstride0;
224     out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride1+2*Nstride0;
225     out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+2*Nstride1+1*Nstride0;
226     } else {
227     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride1;
228     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
229     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1;
230     }
231     }
232     faceNECount+=local_NE1;
233     }
234     totalNECount+=NE1;
235     /* ** elements on boundary 002 (x1=1): */
236     if (local_NE0+e_offset0 == NE0) {
237     #pragma omp parallel for private(i1,k,node0)
238     for (i1=0;i1<local_NE1;i1++) {
239     k=i1+faceNECount;
240     node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1);
241 bcumming 782
242 ksteube 1312 out->FaceElements->Id[k]=(i1+e_offset1)+totalNECount;
243     out->FaceElements->Tag[k]=2;
244     out->FaceElements->Owner[k]=myRank;
245 bcumming 782
246 ksteube 1312 if (useElementsOnFace) {
247     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride0;
248     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1+2*Nstride0;
249     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride1;
250     out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0;
251     out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+Nstride1+2*Nstride0;
252     out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+2*Nstride1+1*Nstride0;
253     out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride1;
254     out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+1*Nstride0;
255     } else {
256     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride0;
257     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1+2*Nstride0;
258     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1+2*Nstride0;
259     }
260     }
261     faceNECount+=local_NE1;
262     }
263     totalNECount+=NE1;
264 bcumming 782 }
265 gross 1818 if (!periodic[1] && (local_NE1>0)) {
266 ksteube 1312 /* ** elements on boundary 010 (x2=0): */
267     if (e_offset1 == 0) {
268     #pragma omp parallel for private(i0,k,node0)
269     for (i0=0;i0<local_NE0;i0++) {
270     k=i0+faceNECount;
271     node0=Nstride0*N_PER_E*(i0+e_offset0);
272    
273     out->FaceElements->Id[k]=e_offset0+i0+totalNECount;
274     out->FaceElements->Tag[k]=10;
275     out->FaceElements->Owner[k]=myRank;
276    
277     if (useElementsOnFace) {
278     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0;
279     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride0;
280     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride1+2*Nstride0;
281     out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride1;
282     out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+1*Nstride0;
283     out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1+2*Nstride0;
284     out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride1+1*Nstride0;
285     out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+Nstride1;
286     } else {
287     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0;
288     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride0;
289     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+1*Nstride0;
290     }
291     }
292     faceNECount+=local_NE0;
293     }
294     totalNECount+=NE0;
295     /* ** elements on boundary 020 (x2=1): */
296     if (local_NE1+e_offset1 == NE1) {
297     #pragma omp parallel for private(i0,k,node0)
298     for (i0=0;i0<local_NE0;i0++) {
299     k=i0+faceNECount;
300     node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1);
301 bcumming 782
302 ksteube 1312 out->FaceElements->Id[k]=i0+e_offset0+totalNECount;
303     out->FaceElements->Tag[k]=20;
304     out->FaceElements->Owner[k]=myRank;
305     if (useElementsOnFace) {
306     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride1+2*Nstride0;
307     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1;
308     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0;
309     out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride0;
310     out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+2*Nstride1+1*Nstride0;
311     out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1;
312     out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+1*Nstride0;
313     out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+Nstride1+2*Nstride0;
314     } else {
315     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride1+2*Nstride0;
316     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1;
317     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride1+1*Nstride0;
318     }
319     }
320     faceNECount+=local_NE0;
321     }
322     totalNECount+=NE0;
323     }
324 gross 2748 }
325     if (Finley_noError()) {
326 ksteube 1312 /* add tag names */
327     Finley_Mesh_addTagMap(out,"top", 20);
328     Finley_Mesh_addTagMap(out,"bottom", 10);
329     Finley_Mesh_addTagMap(out,"left", 1);
330     Finley_Mesh_addTagMap(out,"right", 2);
331 gross 2748 }
332     /* prepare mesh for further calculatuions:*/
333     if (Finley_noError()) {
334 ksteube 1312 Finley_Mesh_resolveNodeIds(out);
335 gross 2748 }
336     if (Finley_noError()) {
337 ksteube 1312 Finley_Mesh_prepare(out, optimize);
338 gross 2748 }
339     if (!Finley_noError()) {
340 ksteube 1312 Finley_Mesh_free(out);
341 gross 2748 }
342     /* free up memory */
343     Finley_ReferenceElementSet_dealloc(refPoints);
344     Finley_ReferenceElementSet_dealloc(refContactElements);
345     Finley_ReferenceElementSet_dealloc(refFaceElements);
346     Finley_ReferenceElementSet_dealloc(refElements);
347     Paso_MPIInfo_free( mpi_info );
348 bcumming 782
349 gross 2748 return out;
350 bcumming 782 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26