/[escript]/trunk/dudley/src/Mesh_tri3.c
ViewVC logotype

Contents of /trunk/dudley/src/Mesh_tri3.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3892 - (show annotations)
Tue Apr 10 08:57:23 2012 UTC (6 years, 11 months ago) by jfenwick
File MIME type: text/plain
File size: 9849 byte(s)
Merged changes across from the attempt2 branch.
This version builds and passes python2 tests.
It also passes most python3 tests.



1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 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 /* Dudley: generates triangular meshes by splitting rectangles */
17
18 /* Generates a numElements[0] x numElements[1] x 2 mesh with first order elements (Tri3) in the rectangle */
19 /* [0,Length[0]] x [0,Length[1]]. order is the desired accuracy of the integration scheme. */
20
21 /**************************************************************/
22
23 #include "TriangularMesh.h"
24
25 Dudley_Mesh *Dudley_TriangularMesh_Tri3(dim_t * numElements,
26 double *Length, index_t order, index_t reduced_order, bool_t optimize)
27 {
28 #define N_PER_E 1
29 #define DIM 2
30 dim_t N0, N1, NE0, NE1, i0, i1, Nstride0 = 0, Nstride1 = 0, local_NE0, local_NE1, local_N0 = 0, local_N1 = 0;
31 index_t offset0 = 0, offset1 = 0, e_offset0 = 0, e_offset1 = 0;
32 dim_t totalNECount, faceNECount, NDOF0 = 0, NDOF1 = 0, NFaceElements;
33 index_t myRank;
34 Dudley_Mesh *out;
35 Esys_MPIInfo *mpi_info = NULL;
36 char name[50];
37 const int LEFTTAG = 1; /* boundary x1=0 */
38 const int RIGHTTAG = 2; /* boundary x1=1 */
39 const int BOTTOMTAG = 10; /* boundary x2=0 */
40 const int TOPTAG = 20; /* boundary x2=1 */
41
42 #ifdef Dudley_TRACE
43 double time0 = Dudley_timer();
44 #endif
45
46 /* get MPI information */
47 mpi_info = Esys_MPIInfo_alloc(MPI_COMM_WORLD);
48 if (!Dudley_noError())
49 {
50 return NULL;
51 }
52 myRank = mpi_info->rank;
53
54 /* set up the global dimensions of the mesh */
55
56 NE0 = MAX(1, numElements[0]);
57 NE1 = MAX(1, numElements[1]);
58 N0 = N_PER_E * NE0 + 1;
59 N1 = N_PER_E * NE1 + 1;
60
61 /* This code was originally copied from Finley's Rec4 constructor.
62 NE? refers to the number of rectangular elements in each direction.
63 The number of nodes produced is the same but the number of non-face elements
64 will double.
65 */
66
67 /* allocate mesh: */
68 sprintf(name, "Triangular %d x %d (x 2) mesh", N0, N1);
69 out = Dudley_Mesh_alloc(name, DIM, mpi_info);
70 if (!Dudley_noError())
71 {
72 Esys_MPIInfo_free(mpi_info);
73 return NULL;
74 }
75 if (Dudley_noError())
76 {
77
78 Dudley_Mesh_setPoints(out, Dudley_ElementFile_alloc(Dudley_Point1, mpi_info));
79 Dudley_Mesh_setFaceElements(out, Dudley_ElementFile_alloc(Dudley_Line2, mpi_info));
80 Dudley_Mesh_setElements(out, Dudley_ElementFile_alloc(Dudley_Tri3, mpi_info));
81 Nstride0 = 1;
82 Nstride1 = N0;
83 if (N1 == MAX(N0, N1))
84 {
85 local_NE0 = NE0;
86 e_offset0 = 0;
87 Esys_MPIInfo_Split(mpi_info, NE1, &local_NE1, &e_offset1);
88 }
89 else
90 {
91 Esys_MPIInfo_Split(mpi_info, NE0, &local_NE0, &e_offset0);
92 local_NE1 = NE1;
93 e_offset1 = 0;
94 }
95 offset0 = e_offset0 * N_PER_E;
96 offset1 = e_offset1 * N_PER_E;
97 local_N0 = local_NE0 > 0 ? local_NE0 * N_PER_E + 1 : 0;
98 local_N1 = local_NE1 > 0 ? local_NE1 * N_PER_E + 1 : 0;
99
100 /* get the number of surface elements */
101
102 NFaceElements = 0;
103 if (local_NE0 > 0)
104 {
105 NDOF0 = N0;
106 if (e_offset0 == 0)
107 NFaceElements += local_NE1;
108 if (local_NE0 + e_offset0 == NE0)
109 NFaceElements += local_NE1;
110 }
111 else
112 {
113 NDOF0 = N0 - 1;
114 }
115 if (local_NE1 > 0)
116 {
117 NDOF1 = N1;
118 if (e_offset1 == 0)
119 NFaceElements += local_NE0;
120 if (local_NE1 + e_offset1 == NE1)
121 NFaceElements += local_NE0;
122 }
123 else
124 {
125 NDOF1 = N1 - 1;
126 }
127
128 /* allocate tables: */
129
130 Dudley_NodeFile_allocTable(out->Nodes, local_N0 * local_N1);
131
132 /* This code was oringinally copied from Finley's rec4 generator
133 We double these numbers because each "rectangle" will be split into
134 two triangles. So the number of nodes is the same but the
135 number of elements will double */
136 Dudley_ElementFile_allocTable(out->Elements, local_NE0 * local_NE1 * 2);
137 Dudley_ElementFile_allocTable(out->FaceElements, NFaceElements);
138
139 }
140 if (Dudley_noError())
141 {
142 dim_t NN;
143 index_t global_adjustment;
144 /* create nodes */
145 #pragma omp parallel for private(i0,i1)
146 for (i1 = 0; i1 < local_N1; i1++)
147 {
148 for (i0 = 0; i0 < local_N0; i0++)
149 {
150 dim_t k = i0 + local_N0 * i1;
151 dim_t global_i0 = i0 + offset0;
152 dim_t 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) + Nstride1 * (global_i1 % NDOF1);
158 }
159 }
160 /* set the elements: */
161 NN = out->Elements->numNodes;
162 global_adjustment = (offset0 + offset1) % 2;
163 #pragma omp parallel for private(i0,i1)
164 for (i1 = 0; i1 < local_NE1; i1++)
165 {
166 for (i0 = 0; i0 < local_NE0; i0++)
167 {
168 index_t a, b, c, d;
169 /* we will split this "rectangle" into two triangles */
170 dim_t k = 2 * (i0 + local_NE0 * i1);
171 index_t node0 = Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (i1 + e_offset1);
172
173 out->Elements->Id[k] = 2 * ((i0 + e_offset0) + NE0 * (i1 + e_offset1));
174 out->Elements->Tag[k] = 0;
175 out->Elements->Owner[k] = myRank;
176 out->Elements->Id[k + 1] = out->Elements->Id[k] + 1;
177 out->Elements->Tag[k + 1] = 0;
178 out->Elements->Owner[k + 1] = myRank;
179
180 /* a,b,c,d gives the nodes in the rectangle in clockwise order */
181 a = node0; b = node0 + Nstride0; c = node0 + Nstride1 + Nstride0; d = node0 + Nstride1;
182 /* For a little bit of variety */
183 if ((global_adjustment + node0) % 2)
184 {
185 out->Elements->Nodes[INDEX2(0, k, NN)] = a;
186 out->Elements->Nodes[INDEX2(1, k, NN)] = b;
187 out->Elements->Nodes[INDEX2(2, k, NN)] = d;
188 out->Elements->Nodes[INDEX2(0, k + 1, NN)] = b;
189 out->Elements->Nodes[INDEX2(1, k + 1, NN)] = c;
190 out->Elements->Nodes[INDEX2(2, k + 1, NN)] = d;
191 }
192 else
193 {
194 out->Elements->Nodes[INDEX2(0, k, NN)] = a;
195 out->Elements->Nodes[INDEX2(1, k, NN)] = b;
196 out->Elements->Nodes[INDEX2(2, k, NN)] = c;
197 out->Elements->Nodes[INDEX2(0, k + 1, NN)] = a;
198 out->Elements->Nodes[INDEX2(1, k + 1, NN)] = c;
199 out->Elements->Nodes[INDEX2(2, k + 1, NN)] = d;
200 }
201 }
202 }
203 /* face elements */
204 NN = out->FaceElements->numNodes;
205 totalNECount = 2 * NE0 * NE1; /* because we have split the rectangles */
206 faceNECount = 0;
207 if (local_NE0 > 0)
208 {
209 /* ** elements on boundary 001 (x1=0): */
210
211 if (e_offset0 == 0)
212 {
213 #pragma omp parallel for private(i1)
214 for (i1 = 0; i1 < local_NE1; i1++)
215 {
216
217 dim_t k = i1 + faceNECount;
218 index_t node0 = Nstride1 * N_PER_E * (i1 + e_offset1);
219
220 out->FaceElements->Id[k] = i1 + e_offset1 + totalNECount;
221 out->FaceElements->Tag[k] = LEFTTAG;
222 out->FaceElements->Owner[k] = myRank;
223 out->FaceElements->Nodes[INDEX2(0, k, NN)] = node0 + Nstride1;
224 out->FaceElements->Nodes[INDEX2(1, k, NN)] = node0;
225 }
226 faceNECount += local_NE1;
227 }
228 totalNECount += NE1;
229 /* ** elements on boundary 002 (x1=1): */
230 if (local_NE0 + e_offset0 == NE0)
231 {
232 #pragma omp parallel for private(i1)
233 for (i1 = 0; i1 < local_NE1; i1++)
234 {
235 dim_t k = i1 + faceNECount;
236 index_t node0 = Nstride0 * N_PER_E * (NE0 - 1) + Nstride1 * N_PER_E * (i1 + e_offset1);
237
238 out->FaceElements->Id[k] = (i1 + e_offset1) + totalNECount;
239 out->FaceElements->Tag[k] = RIGHTTAG;
240 out->FaceElements->Owner[k] = myRank;
241 out->FaceElements->Nodes[INDEX2(0, k, NN)] = node0 + Nstride0;
242 out->FaceElements->Nodes[INDEX2(1, k, NN)] = node0 + Nstride1 + Nstride0;
243 }
244 faceNECount += local_NE1;
245 }
246 totalNECount += NE1;
247 }
248 if (local_NE1 > 0)
249 {
250 /* ** elements on boundary 010 (x2=0): */
251 if (e_offset1 == 0)
252 {
253 #pragma omp parallel for private(i0)
254 for (i0 = 0; i0 < local_NE0; i0++)
255 {
256 dim_t k = i0 + faceNECount;
257 index_t node0 = Nstride0 * N_PER_E * (i0 + e_offset0);
258
259 out->FaceElements->Id[k] = e_offset0 + i0 + totalNECount;
260 out->FaceElements->Tag[k] = BOTTOMTAG;
261 out->FaceElements->Owner[k] = myRank;
262
263 out->FaceElements->Nodes[INDEX2(0, k, NN)] = node0;
264 out->FaceElements->Nodes[INDEX2(1, k, NN)] = node0 + Nstride0;
265 }
266 faceNECount += local_NE0;
267 }
268 totalNECount += NE0;
269 /* ** elements on boundary 020 (x2=1): */
270 if (local_NE1 + e_offset1 == NE1)
271 {
272 #pragma omp parallel for private(i0)
273 for (i0 = 0; i0 < local_NE0; i0++)
274 {
275 dim_t k = i0 + faceNECount;
276 index_t node0 = Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (NE1 - 1);
277
278 out->FaceElements->Id[k] = i0 + e_offset0 + totalNECount;
279 out->FaceElements->Tag[k] = TOPTAG;
280 out->FaceElements->Owner[k] = myRank;
281
282 out->FaceElements->Nodes[INDEX2(0, k, NN)] = node0 + Nstride1 + Nstride0;
283 out->FaceElements->Nodes[INDEX2(1, k, NN)] = node0 + Nstride1;
284 /*printf("E=%d: %d=%d %d=%d\n",k,INDEX2(0,k,NN),out->FaceElements->Nodes[INDEX2(0,k,NN)],
285 INDEX2(1,k,NN),out->FaceElements->Nodes[INDEX2(1,k,NN)]); */
286 }
287 faceNECount += local_NE0;
288 }
289 totalNECount += NE0;
290 }
291 }
292 if (Dudley_noError())
293 {
294 /* add tag names */
295 Dudley_Mesh_addTagMap(out, "top", TOPTAG);
296 Dudley_Mesh_addTagMap(out, "bottom", BOTTOMTAG);
297 Dudley_Mesh_addTagMap(out, "left", LEFTTAG);
298 Dudley_Mesh_addTagMap(out, "right", RIGHTTAG);
299 }
300 /* prepare mesh for further calculatuions: */
301 if (Dudley_noError())
302 {
303 Dudley_Mesh_resolveNodeIds(out);
304 }
305 if (Dudley_noError())
306 {
307 Dudley_Mesh_prepare(out, optimize);
308 }
309
310 /* free up memory */
311 Esys_MPIInfo_free(mpi_info);
312
313 return out;
314 }

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/pasowrap/dudley/src/Mesh_tri3.c:3661-3674 /branches/py3_attempt2/dudley/src/Mesh_tri3.c:3871-3891 /branches/ripleygmg_from_3668/dudley/src/Mesh_tri3.c:3669-3791 /branches/scons_revamp_from_3210/dudley/src/Mesh_tri3.c:3212-3243

  ViewVC Help
Powered by ViewVC 1.1.26