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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1986 - (hide annotations)
Thu Nov 6 23:33:14 2008 UTC (11 years ago) by jfenwick
File MIME type: text/plain
File size: 11954 byte(s)
Warning removal
1 jgs 150
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4     * Copyright (c) 2003-2008 by University of Queensland
5     * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13 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 first order elements (Rec4) in the rectangle */
20     /* [0,Length[0]] x [0,Length[1]]. order is the desired accuracy of the integration scheme. */
21    
22    
23     /**************************************************************/
24    
25     #include "RectangularMesh.h"
26    
27    
28 ksteube 1312 Finley_Mesh* Finley_RectangularMesh_Rec4(dim_t* numElements,
29     double* Length,
30     bool_t* periodic,
31     index_t order,
32     index_t reduced_order,
33     bool_t useElementsOnFace,
34     bool_t useFullElementOrder,
35     bool_t optimize)
36 bcumming 730 {
37 ksteube 1312 #define N_PER_E 1
38     #define DIM 2
39     dim_t N0,N1,NE0,NE1,i0,i1,k,Nstride0,Nstride1, local_NE0, local_NE1, local_N0, local_N1, global_i0, global_i1;
40     index_t offset0, offset1, e_offset0, e_offset1;
41     dim_t totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NN;
42     index_t node0, myRank;
43     Finley_Mesh* out;
44     Paso_MPIInfo *mpi_info = NULL;
45     char name[50];
46 jfenwick 1986 #ifdef Finley_TRACE
47 ksteube 1312 double time0=Finley_timer();
48 jfenwick 1986 #endif
49 bcumming 730
50 ksteube 1312 /* get MPI information */
51     mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
52     if (! Finley_noError()) {
53     return NULL;
54 bcumming 730 }
55 ksteube 1312 myRank=mpi_info->rank;
56 bcumming 730
57 ksteube 1312 /* set up the global dimensions of the mesh */
58 bcumming 730
59 jgs 82 NE0=MAX(1,numElements[0]);
60     NE1=MAX(1,numElements[1]);
61 ksteube 1312 N0=N_PER_E*NE0+1;
62     N1=N_PER_E*NE1+1;
63 jgs 82
64 ksteube 1312 /* allocate mesh: */
65 jgs 82 sprintf(name,"Rectangular %d x %d mesh",N0,N1);
66 ksteube 1312 out=Finley_Mesh_alloc(name,DIM,order, reduced_order, mpi_info);
67     if (! Finley_noError()) {
68     Paso_MPIInfo_free( mpi_info );
69 bcumming 782 return NULL;
70     }
71 bcumming 787
72 ksteube 1312 Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Rec4,out->order,out->reduced_order,mpi_info));
73 jgs 82 if (useElementsOnFace) {
74 ksteube 1312 Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Rec4Face,
75     out->order,
76     out->reduced_order,
77     mpi_info));
78     Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Rec4Face_Contact,
79     out->order,
80     out->reduced_order,
81     mpi_info));
82 jgs 82 } else {
83 ksteube 1312 Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Line2,
84     out->order,
85     out->reduced_order,
86     mpi_info));
87     Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Line2_Contact,
88     out->order,
89     out->reduced_order,
90     mpi_info));
91 jgs 82 }
92 ksteube 1312 Finley_Mesh_setPoints(out,Finley_ElementFile_alloc(Point1,
93     out->order,
94     out->reduced_order,
95     mpi_info));
96 jgs 150 if (! Finley_noError()) {
97 ksteube 1312 Paso_MPIInfo_free( mpi_info );
98     Finley_Mesh_free(out);
99 jgs 82 return NULL;
100     }
101 jgs 153
102 ksteube 1312 /* work out the largest dimension */
103     if (N1==MAX(N0,N1)) {
104     Nstride0=1;
105     Nstride1=N0;
106     local_NE0=NE0;
107     e_offset0=0;
108     Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1);
109 jgs 82 } else {
110 ksteube 1312 Nstride0=N1;
111     Nstride1=1;
112     Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0);
113     local_NE1=NE1;
114     e_offset1=0;
115 jgs 82 }
116 ksteube 1312 offset0=e_offset0*N_PER_E;
117     offset1=e_offset1*N_PER_E;
118 gross 1733 local_N0=local_NE0>0 ? local_NE0*N_PER_E+1 : 0;
119     local_N1=local_NE1>0 ? local_NE1*N_PER_E+1 : 0;
120 jgs 82
121 ksteube 1312 /* get the number of surface elements */
122 jgs 82
123 ksteube 1312 NFaceElements=0;
124 gross 1733 if (!periodic[0] && (local_NE0>0)) {
125 ksteube 1312 NDOF0=N0;
126     if (e_offset0 == 0) NFaceElements+=local_NE1;
127     if (local_NE0+e_offset0 == NE0) NFaceElements+=local_NE1;
128 gross 964 } else {
129 ksteube 1312 NDOF0=N0-1;
130 jgs 82 }
131 gross 1733 if (!periodic[1] && (local_NE1>0)) {
132 ksteube 1312 NDOF1=N1;
133     if (e_offset1 == 0) NFaceElements+=local_NE0;
134     if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0;
135 bcumming 730 } else {
136     NDOF1=N1-1;
137     }
138    
139 ksteube 1312 /* allocate tables: */
140 bcumming 730
141 ksteube 1312 Finley_NodeFile_allocTable(out->Nodes,local_N0*local_N1);
142     Finley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1);
143 bcumming 730 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
144    
145 ksteube 1312 if (Finley_noError()) {
146     /* create nodes */
147     #pragma omp parallel for private(i0,i1,k,global_i0,global_i1)
148     for (i1=0;i1<local_N1;i1++) {
149     for (i0=0;i0<local_N0;i0++) {
150     k=i0+local_N0*i1;
151     global_i0=i0+offset0;
152     global_i1=i1+offset1;
153     out->Nodes->Coordinates[INDEX2(0,k,DIM)]=DBLE(global_i0)/DBLE(N0-1)*Length[0];
154     out->Nodes->Coordinates[INDEX2(1,k,DIM)]=DBLE(global_i1)/DBLE(N1-1)*Length[1];
155     out->Nodes->Id[k]=Nstride0*global_i0+Nstride1*global_i1;
156     out->Nodes->Tag[k]=0;
157     out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
158     +Nstride1*(global_i1%NDOF1);
159     }
160     }
161     /* set the elements: */
162     NN=out->Elements->numNodes;
163     #pragma omp parallel for private(i0,i1,k,node0)
164     for (i1=0;i1<local_NE1;i1++) {
165     for (i0=0;i0<local_NE0;i0++) {
166    
167     k=i0+local_NE0*i1;
168     node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
169 bcumming 787
170 ksteube 1312 out->Elements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1);
171     out->Elements->Tag[k]=0;
172     out->Elements->Owner[k]=myRank;
173 bcumming 782
174 ksteube 1312 out->Elements->Nodes[INDEX2(0,k,NN)]=node0;
175     out->Elements->Nodes[INDEX2(1,k,NN)]=node0+Nstride0;
176     out->Elements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1+Nstride0;
177     out->Elements->Nodes[INDEX2(3,k,NN)]=node0+Nstride1;
178     }
179     }
180     /* face elements */
181     NN=out->FaceElements->numNodes;
182     totalNECount=NE0*NE1;
183     faceNECount=0;
184 gross 1733 if (!periodic[0] && (local_NE0>0)) {
185 ksteube 1312 /* ** elements on boundary 001 (x1=0): */
186    
187     if (e_offset0 == 0) {
188     #pragma omp parallel for private(i1,k,node0)
189     for (i1=0;i1<local_NE1;i1++) {
190    
191     k=i1+faceNECount;
192     node0=Nstride1*N_PER_E*(i1+e_offset1);
193 bcumming 730
194 ksteube 1312 out->FaceElements->Id[k]=i1+e_offset1+totalNECount;
195     out->FaceElements->Tag[k]=1;
196     out->FaceElements->Owner[k]=myRank;
197     if (useElementsOnFace) {
198     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1;
199     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
200     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride0;
201     out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride1+Nstride0;
202     } else {
203     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1;
204     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
205     }
206     }
207     faceNECount+=local_NE1;
208     }
209     totalNECount+=NE1;
210     /* ** elements on boundary 002 (x1=1): */
211     if (local_NE0+e_offset0 == NE0) {
212     #pragma omp parallel for private(i1,k,node0)
213     for (i1=0;i1<local_NE1;i1++) {
214     k=i1+faceNECount;
215     node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1);
216 bcumming 730
217 ksteube 1312 out->FaceElements->Id[k]=(i1+e_offset1)+totalNECount;
218     out->FaceElements->Tag[k]=2;
219     out->FaceElements->Owner[k]=myRank;
220 bcumming 730
221 ksteube 1312 if (useElementsOnFace) {
222     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride0;
223     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1+Nstride0;
224     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1;
225     out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0;
226     } else {
227     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride0;
228     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1+Nstride0;
229     }
230     }
231     faceNECount+=local_NE1;
232     }
233     totalNECount+=NE1;
234     }
235 gross 1733 if (!periodic[1] && (local_NE1>0)) {
236 ksteube 1312 /* ** elements on boundary 010 (x2=0): */
237     if (e_offset1 == 0) {
238     #pragma omp parallel for private(i0,k,node0)
239     for (i0=0;i0<local_NE0;i0++) {
240     k=i0+faceNECount;
241     node0=Nstride0*N_PER_E*(i0+e_offset0);
242    
243     out->FaceElements->Id[k]=e_offset0+i0+totalNECount;
244     out->FaceElements->Tag[k]=10;
245     out->FaceElements->Owner[k]=myRank;
246    
247     if (useElementsOnFace) {
248     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0;
249     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride0;
250     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1+Nstride0;
251     out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride1;
252     } else {
253     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0;
254     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride0;
255     }
256     }
257     faceNECount+=local_NE0;
258     }
259     totalNECount+=NE0;
260     /* ** elements on boundary 020 (x2=1): */
261     if (local_NE1+e_offset1 == NE1) {
262     #pragma omp parallel for private(i0,k,node0)
263     for (i0=0;i0<local_NE0;i0++) {
264     k=i0+faceNECount;
265     node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1);
266    
267     out->FaceElements->Id[k]=i0+e_offset0+totalNECount;
268     out->FaceElements->Tag[k]=20;
269     out->FaceElements->Owner[k]=myRank;
270     if (useElementsOnFace) {
271     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1+Nstride0;
272     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1;
273     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0;
274     out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride0;
275     } else {
276     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1+Nstride0;
277     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1;
278     }
279     }
280     faceNECount+=local_NE0;
281     }
282     totalNECount+=NE0;
283     }
284     /* add tag names */
285     Finley_Mesh_addTagMap(out,"top", 20);
286     Finley_Mesh_addTagMap(out,"bottom", 10);
287     Finley_Mesh_addTagMap(out,"left", 1);
288     Finley_Mesh_addTagMap(out,"right", 2);
289    
290     /* prepare mesh for further calculatuions:*/
291     if (Finley_noError()) {
292     Finley_Mesh_resolveNodeIds(out);
293     }
294     if (Finley_noError()) {
295     Finley_Mesh_prepare(out, optimize);
296     }
297 bcumming 730 }
298    
299 ksteube 1312 if (!Finley_noError()) {
300     Finley_Mesh_free(out);
301 bcumming 730 }
302 bcumming 782 /* free up memory */
303 ksteube 1312 Paso_MPIInfo_free( mpi_info );
304 bcumming 782 #ifdef Finley_TRACE
305     printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
306     #endif
307 ksteube 1312
308     return out;
309 bcumming 730 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26