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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26