/[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 1811 - (show 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
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 double time0=Finley_timer();
47
48 /* get MPI information */
49 mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
50 if (! Finley_noError()) {
51 return NULL;
52 }
53 myRank=mpi_info->rank;
54
55 /* set up the global dimensions of the mesh */
56
57 NE0=MAX(1,numElements[0]);
58 NE1=MAX(1,numElements[1]);
59 N0=N_PER_E*NE0+1;
60 N1=N_PER_E*NE1+1;
61
62 /* allocate mesh: */
63 sprintf(name,"Rectangular %d x %d mesh",N0,N1);
64 out=Finley_Mesh_alloc(name,DIM,order, reduced_order, mpi_info);
65 if (! Finley_noError()) {
66 Paso_MPIInfo_free( mpi_info );
67 return NULL;
68 }
69
70 Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Rec4,out->order,out->reduced_order,mpi_info));
71 if (useElementsOnFace) {
72 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 } else {
81 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 }
90 Finley_Mesh_setPoints(out,Finley_ElementFile_alloc(Point1,
91 out->order,
92 out->reduced_order,
93 mpi_info));
94 if (! Finley_noError()) {
95 Paso_MPIInfo_free( mpi_info );
96 Finley_Mesh_free(out);
97 return NULL;
98 }
99
100 /* 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 } else {
108 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 }
114 offset0=e_offset0*N_PER_E;
115 offset1=e_offset1*N_PER_E;
116 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
119 /* get the number of surface elements */
120
121 NFaceElements=0;
122 if (!periodic[0] && (local_NE0>0)) {
123 NDOF0=N0;
124 if (e_offset0 == 0) NFaceElements+=local_NE1;
125 if (local_NE0+e_offset0 == NE0) NFaceElements+=local_NE1;
126 } else {
127 NDOF0=N0-1;
128 }
129 if (!periodic[1] && (local_NE1>0)) {
130 NDOF1=N1;
131 if (e_offset1 == 0) NFaceElements+=local_NE0;
132 if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0;
133 } else {
134 NDOF1=N1-1;
135 }
136
137 /* allocate tables: */
138
139 Finley_NodeFile_allocTable(out->Nodes,local_N0*local_N1);
140 Finley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1);
141 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
142
143 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
168 out->Elements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1);
169 out->Elements->Tag[k]=0;
170 out->Elements->Owner[k]=myRank;
171
172 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 if (!periodic[0] && (local_NE0>0)) {
183 /* ** 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
192 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
215 out->FaceElements->Id[k]=(i1+e_offset1)+totalNECount;
216 out->FaceElements->Tag[k]=2;
217 out->FaceElements->Owner[k]=myRank;
218
219 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 if (!periodic[1] && (local_NE1>0)) {
234 /* ** 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 }
296
297 if (!Finley_noError()) {
298 Finley_Mesh_free(out);
299 }
300 /* free up memory */
301 Paso_MPIInfo_free( mpi_info );
302 #ifdef Finley_TRACE
303 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
304 #endif
305
306 return out;
307 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26