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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1986 - (show annotations)
Thu Nov 6 23:33:14 2008 UTC (10 years, 10 months ago) by jfenwick
File MIME type: text/plain
File size: 11954 byte(s)
Warning removal
1
2 /*******************************************************
3 *
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
14
15 /**************************************************************/
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 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 {
37 #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 #ifdef Finley_TRACE
47 double time0=Finley_timer();
48 #endif
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 /* set up the global dimensions of the mesh */
58
59 NE0=MAX(1,numElements[0]);
60 NE1=MAX(1,numElements[1]);
61 N0=N_PER_E*NE0+1;
62 N1=N_PER_E*NE1+1;
63
64 /* allocate mesh: */
65 sprintf(name,"Rectangular %d x %d mesh",N0,N1);
66 out=Finley_Mesh_alloc(name,DIM,order, reduced_order, mpi_info);
67 if (! Finley_noError()) {
68 Paso_MPIInfo_free( mpi_info );
69 return NULL;
70 }
71
72 Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Rec4,out->order,out->reduced_order,mpi_info));
73 if (useElementsOnFace) {
74 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 } else {
83 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 }
92 Finley_Mesh_setPoints(out,Finley_ElementFile_alloc(Point1,
93 out->order,
94 out->reduced_order,
95 mpi_info));
96 if (! Finley_noError()) {
97 Paso_MPIInfo_free( mpi_info );
98 Finley_Mesh_free(out);
99 return NULL;
100 }
101
102 /* 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 } else {
110 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 }
116 offset0=e_offset0*N_PER_E;
117 offset1=e_offset1*N_PER_E;
118 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
121 /* get the number of surface elements */
122
123 NFaceElements=0;
124 if (!periodic[0] && (local_NE0>0)) {
125 NDOF0=N0;
126 if (e_offset0 == 0) NFaceElements+=local_NE1;
127 if (local_NE0+e_offset0 == NE0) NFaceElements+=local_NE1;
128 } else {
129 NDOF0=N0-1;
130 }
131 if (!periodic[1] && (local_NE1>0)) {
132 NDOF1=N1;
133 if (e_offset1 == 0) NFaceElements+=local_NE0;
134 if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0;
135 } else {
136 NDOF1=N1-1;
137 }
138
139 /* allocate tables: */
140
141 Finley_NodeFile_allocTable(out->Nodes,local_N0*local_N1);
142 Finley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1);
143 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
144
145 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
170 out->Elements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1);
171 out->Elements->Tag[k]=0;
172 out->Elements->Owner[k]=myRank;
173
174 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 if (!periodic[0] && (local_NE0>0)) {
185 /* ** 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
194 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
217 out->FaceElements->Id[k]=(i1+e_offset1)+totalNECount;
218 out->FaceElements->Tag[k]=2;
219 out->FaceElements->Owner[k]=myRank;
220
221 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 if (!periodic[1] && (local_NE1>0)) {
236 /* ** 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 }
298
299 if (!Finley_noError()) {
300 Finley_Mesh_free(out);
301 }
302 /* free up memory */
303 Paso_MPIInfo_free( mpi_info );
304 #ifdef Finley_TRACE
305 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
306 #endif
307
308 return out;
309 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26