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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26