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

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26