/[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 1811 - (hide annotations)
Thu Sep 25 23:11:13 2008 UTC (11 years ago) by ksteube
File MIME type: text/plain
File size: 11923 byte(s)
Copyright updated in all files

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26