/[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 2748 - (show annotations)
Tue Nov 17 07:32:59 2009 UTC (9 years, 8 months ago) by gross
File MIME type: text/plain
File size: 11326 byte(s)
Macro elements are implemented now. VTK writer for macro elements still needs testing.
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2009 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=0,Nstride1=0, local_NE0, local_NE1, local_N0=0, local_N1=0, global_i0, global_i1;
40 index_t offset0=0, offset1=0, e_offset0=0, e_offset1=0;
41 dim_t totalNECount,faceNECount,NDOF0=0,NDOF1=0,NFaceElements,NN;
42 Finley_ReferenceElementSet *refPoints=NULL, *refContactElements=NULL, *refFaceElements=NULL, *refElements=NULL;
43 index_t node0, myRank;
44 Finley_Mesh* out;
45 Paso_MPIInfo *mpi_info = NULL;
46 char name[50];
47 #ifdef Finley_TRACE
48 double time0=Finley_timer();
49 #endif
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 /* set up the global dimensions of the mesh */
59
60 NE0=MAX(1,numElements[0]);
61 NE1=MAX(1,numElements[1]);
62 N0=N_PER_E*NE0+1;
63 N1=N_PER_E*NE1+1;
64
65 /* allocate mesh: */
66 sprintf(name,"Rectangular %d x %d mesh",N0,N1);
67 out=Finley_Mesh_alloc(name,DIM,order, reduced_order, mpi_info);
68 if (! Finley_noError()) {
69 Paso_MPIInfo_free( mpi_info );
70 return NULL;
71 }
72 refElements= Finley_ReferenceElementSet_alloc(Rec4,out->order,out->reduced_order);
73 if (useElementsOnFace) {
74 refFaceElements=Finley_ReferenceElementSet_alloc(Rec4Face, out->order, out->reduced_order);
75 refContactElements=Finley_ReferenceElementSet_alloc(Rec4Face_Contact, out->order, out->reduced_order);
76 } else {
77 refFaceElements=Finley_ReferenceElementSet_alloc(Line2, out->order, out->reduced_order);
78 refContactElements=Finley_ReferenceElementSet_alloc(Line2_Contact, out->order, out->reduced_order);
79 }
80 refPoints=Finley_ReferenceElementSet_alloc(Point1, out->order, out->reduced_order);
81
82 if ( Finley_noError()) {
83
84 Finley_Mesh_setPoints(out,Finley_ElementFile_alloc(refPoints, mpi_info));
85 Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(refContactElements, mpi_info));
86 Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(refFaceElements, mpi_info));
87 Finley_Mesh_setElements(out,Finley_ElementFile_alloc(refElements, mpi_info));
88
89 /* work out the largest dimension */
90 if (N1==MAX(N0,N1)) {
91 Nstride0=1;
92 Nstride1=N0;
93 local_NE0=NE0;
94 e_offset0=0;
95 Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1);
96 } else {
97 Nstride0=N1;
98 Nstride1=1;
99 Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0);
100 local_NE1=NE1;
101 e_offset1=0;
102 }
103 offset0=e_offset0*N_PER_E;
104 offset1=e_offset1*N_PER_E;
105 local_N0=local_NE0>0 ? local_NE0*N_PER_E+1 : 0;
106 local_N1=local_NE1>0 ? local_NE1*N_PER_E+1 : 0;
107
108 /* get the number of surface elements */
109
110 NFaceElements=0;
111 if (!periodic[0] && (local_NE0>0)) {
112 NDOF0=N0;
113 if (e_offset0 == 0) NFaceElements+=local_NE1;
114 if (local_NE0+e_offset0 == NE0) NFaceElements+=local_NE1;
115 } else {
116 NDOF0=N0-1;
117 }
118 if (!periodic[1] && (local_NE1>0)) {
119 NDOF1=N1;
120 if (e_offset1 == 0) NFaceElements+=local_NE0;
121 if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0;
122 } else {
123 NDOF1=N1-1;
124 }
125
126 /* allocate tables: */
127
128
129 Finley_NodeFile_allocTable(out->Nodes,local_N0*local_N1);
130 Finley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1);
131 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
132
133 }
134 if (Finley_noError()) {
135 /* create nodes */
136 #pragma omp parallel for private(i0,i1,k,global_i0,global_i1)
137 for (i1=0;i1<local_N1;i1++) {
138 for (i0=0;i0<local_N0;i0++) {
139 k=i0+local_N0*i1;
140 global_i0=i0+offset0;
141 global_i1=i1+offset1;
142 out->Nodes->Coordinates[INDEX2(0,k,DIM)]=DBLE(global_i0)/DBLE(N0-1)*Length[0];
143 out->Nodes->Coordinates[INDEX2(1,k,DIM)]=DBLE(global_i1)/DBLE(N1-1)*Length[1];
144 out->Nodes->Id[k]=Nstride0*global_i0+Nstride1*global_i1;
145 out->Nodes->Tag[k]=0;
146 out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
147 +Nstride1*(global_i1%NDOF1);
148 }
149 }
150 /* set the elements: */
151 NN=out->Elements->numNodes;
152 #pragma omp parallel for private(i0,i1,k,node0)
153 for (i1=0;i1<local_NE1;i1++) {
154 for (i0=0;i0<local_NE0;i0++) {
155
156 k=i0+local_NE0*i1;
157 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
158
159 out->Elements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1);
160 out->Elements->Tag[k]=0;
161 out->Elements->Owner[k]=myRank;
162
163 out->Elements->Nodes[INDEX2(0,k,NN)]=node0;
164 out->Elements->Nodes[INDEX2(1,k,NN)]=node0+Nstride0;
165 out->Elements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1+Nstride0;
166 out->Elements->Nodes[INDEX2(3,k,NN)]=node0+Nstride1;
167 }
168 }
169 /* face elements */
170 NN=out->FaceElements->numNodes;
171 totalNECount=NE0*NE1;
172 faceNECount=0;
173 if (!periodic[0] && (local_NE0>0)) {
174 /* ** elements on boundary 001 (x1=0): */
175
176 if (e_offset0 == 0) {
177 #pragma omp parallel for private(i1,k,node0)
178 for (i1=0;i1<local_NE1;i1++) {
179
180 k=i1+faceNECount;
181 node0=Nstride1*N_PER_E*(i1+e_offset1);
182
183 out->FaceElements->Id[k]=i1+e_offset1+totalNECount;
184 out->FaceElements->Tag[k]=1;
185 out->FaceElements->Owner[k]=myRank;
186 if (useElementsOnFace) {
187 out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1;
188 out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
189 out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride0;
190 out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride1+Nstride0;
191 } else {
192 out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1;
193 out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
194 }
195 }
196 faceNECount+=local_NE1;
197 }
198 totalNECount+=NE1;
199 /* ** elements on boundary 002 (x1=1): */
200 if (local_NE0+e_offset0 == NE0) {
201 #pragma omp parallel for private(i1,k,node0)
202 for (i1=0;i1<local_NE1;i1++) {
203 k=i1+faceNECount;
204 node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1);
205
206 out->FaceElements->Id[k]=(i1+e_offset1)+totalNECount;
207 out->FaceElements->Tag[k]=2;
208 out->FaceElements->Owner[k]=myRank;
209
210 if (useElementsOnFace) {
211 out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride0;
212 out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1+Nstride0;
213 out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1;
214 out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0;
215 } else {
216 out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride0;
217 out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1+Nstride0;
218 }
219 }
220 faceNECount+=local_NE1;
221 }
222 totalNECount+=NE1;
223 }
224 if (!periodic[1] && (local_NE1>0)) {
225 /* ** elements on boundary 010 (x2=0): */
226 if (e_offset1 == 0) {
227 #pragma omp parallel for private(i0,k,node0)
228 for (i0=0;i0<local_NE0;i0++) {
229 k=i0+faceNECount;
230 node0=Nstride0*N_PER_E*(i0+e_offset0);
231
232 out->FaceElements->Id[k]=e_offset0+i0+totalNECount;
233 out->FaceElements->Tag[k]=10;
234 out->FaceElements->Owner[k]=myRank;
235
236 if (useElementsOnFace) {
237 out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0;
238 out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride0;
239 out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1+Nstride0;
240 out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride1;
241 } else {
242 out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0;
243 out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride0;
244 }
245 }
246 faceNECount+=local_NE0;
247 }
248 totalNECount+=NE0;
249 /* ** elements on boundary 020 (x2=1): */
250 if (local_NE1+e_offset1 == NE1) {
251 #pragma omp parallel for private(i0,k,node0)
252 for (i0=0;i0<local_NE0;i0++) {
253 k=i0+faceNECount;
254 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1);
255
256 out->FaceElements->Id[k]=i0+e_offset0+totalNECount;
257 out->FaceElements->Tag[k]=20;
258 out->FaceElements->Owner[k]=myRank;
259 if (useElementsOnFace) {
260 out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1+Nstride0;
261 out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1;
262 out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0;
263 out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride0;
264 } else {
265 out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1+Nstride0;
266 out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1;
267 }
268 }
269 faceNECount+=local_NE0;
270 }
271 totalNECount+=NE0;
272 }
273 }
274 if (Finley_noError()) {
275 /* add tag names */
276 Finley_Mesh_addTagMap(out,"top", 20);
277 Finley_Mesh_addTagMap(out,"bottom", 10);
278 Finley_Mesh_addTagMap(out,"left", 1);
279 Finley_Mesh_addTagMap(out,"right", 2);
280 }
281 /* prepare mesh for further calculatuions:*/
282 if (Finley_noError()) {
283 Finley_Mesh_resolveNodeIds(out);
284 }
285 if (Finley_noError()) {
286 Finley_Mesh_prepare(out, optimize);
287 }
288
289 /* free up memory */
290 Finley_ReferenceElementSet_dealloc(refPoints);
291 Finley_ReferenceElementSet_dealloc(refContactElements);
292 Finley_ReferenceElementSet_dealloc(refFaceElements);
293 Finley_ReferenceElementSet_dealloc(refElements);
294 Paso_MPIInfo_free( mpi_info );
295
296 return out;
297 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26