/[escript]/branches/trilinos_from_5897/dudley/src/Mesh_tri3.cpp
ViewVC logotype

Contents of /branches/trilinos_from_5897/dudley/src/Mesh_tri3.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6079 - (show annotations)
Mon Mar 21 12:22:38 2016 UTC (2 years, 11 months ago) by caltinay
File size: 9672 byte(s)
Big commit - making dudley much more like finley to make it more
managable. Fixed quite a few issues that had been fixed in finley.
Disposed of all ReducedNode/ReducedDOF entities that dudley never supported.
Compiles and passes tests.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17 #include "TriangularMesh.h"
18
19 using escript::DataTypes::real_t;
20
21 namespace dudley {
22
23 Mesh* TriangularMesh_Tri3(const dim_t* numElements, const double* length,
24 bool optimize, escript::JMPI mpiInfo)
25 {
26 const int DIM = 2;
27 const int LEFTTAG = 1; // boundary x1=0
28 const int RIGHTTAG = 2; // boundary x1=1
29 const int BOTTOMTAG = 10; // boundary x2=0
30 const int TOPTAG = 20; // boundary x2=1
31
32 #ifdef Dudley_TRACE
33 double time0 = Dudley_timer();
34 #endif
35
36 const int myRank = mpiInfo->rank;
37
38 // set up the global dimensions of the mesh
39 const dim_t NE0 = std::max(1, numElements[0]);
40 const dim_t NE1 = std::max(1, numElements[1]);
41 const dim_t N0 = NE0 + 1;
42 const dim_t N1 = NE1 + 1;
43
44 // This code was originally copied from Finley's Rec4 constructor.
45 // NE? refers to the number of rectangular elements in each direction.
46 // The number of nodes produced is the same but the number of non-face
47 // elements will double since each "rectangle" is split into two triangles.
48
49 // allocate mesh
50 std::stringstream name;
51 name << "Triangular " << N0 << " x " << N1 << " (x 2) mesh";
52 Mesh* out = new Mesh(name.str(), DIM, mpiInfo);
53
54 out->setPoints(new ElementFile(Dudley_Point1, mpiInfo));
55 out->setFaceElements(new ElementFile(Dudley_Line2, mpiInfo));
56 out->setElements(new ElementFile(Dudley_Tri3, mpiInfo));
57
58 const dim_t Nstride0 = 1;
59 const dim_t Nstride1 = N0;
60 dim_t local_NE0, local_NE1;
61 index_t e_offset0 = 0, e_offset1 = 0;
62 if (N1 == std::max(N0, N1)) {
63 local_NE0 = NE0;
64 e_offset0 = 0;
65 mpiInfo->split(NE1, &local_NE1, &e_offset1);
66 } else {
67 mpiInfo->split(NE0, &local_NE0, &e_offset0);
68 local_NE1 = NE1;
69 e_offset1 = 0;
70 }
71 const index_t offset0 = e_offset0;
72 const index_t offset1 = e_offset1;
73 const dim_t local_N0 = local_NE0 > 0 ? local_NE0 + 1 : 0;
74 const dim_t local_N1 = local_NE1 > 0 ? local_NE1 + 1 : 0;
75
76 // get the number of surface elements
77 dim_t NFaceElements = 0;
78 dim_t NDOF0, NDOF1;
79 if (local_NE0 > 0) {
80 NDOF0 = N0;
81 if (e_offset0 == 0)
82 NFaceElements += local_NE1;
83 if (local_NE0 + e_offset0 == NE0)
84 NFaceElements += local_NE1;
85 } else {
86 NDOF0 = N0 - 1;
87 }
88 if (local_NE1 > 0) {
89 NDOF1 = N1;
90 if (e_offset1 == 0)
91 NFaceElements += local_NE0;
92 if (local_NE1 + e_offset1 == NE1)
93 NFaceElements += local_NE0;
94 } else {
95 NDOF1 = N1 - 1;
96 }
97
98 out->Nodes->allocTable(local_N0 * local_N1);
99 out->Elements->allocTable(local_NE0 * local_NE1 * 2);
100 out->FaceElements->allocTable(NFaceElements);
101
102 // create nodes
103 #pragma omp parallel for
104 for (index_t i1 = 0; i1 < local_N1; i1++) {
105 for (index_t i0 = 0; i0 < local_N0; i0++) {
106 const dim_t k = i0 + local_N0 * i1;
107 const dim_t global_i0 = i0 + offset0;
108 const dim_t global_i1 = i1 + offset1;
109 out->Nodes->Coordinates[INDEX2(0, k, DIM)] = (real_t)global_i0 / (real_t)(N0 - 1) * length[0];
110 out->Nodes->Coordinates[INDEX2(1, k, DIM)] = (real_t)global_i1 / (real_t)(N1 - 1) * length[1];
111 out->Nodes->Id[k] = Nstride0 * global_i0 + Nstride1 * global_i1;
112 out->Nodes->Tag[k] = 0;
113 out->Nodes->globalDegreesOfFreedom[k] =
114 Nstride0 * (global_i0 % NDOF0)
115 + Nstride1 * (global_i1 % NDOF1);
116 }
117 }
118
119 // set the elements
120 dim_t NN = out->Elements->numNodes;
121 const index_t global_adjustment = (offset0 + offset1) % 2;
122
123 #pragma omp parallel for
124 for (index_t i1 = 0; i1 < local_NE1; i1++) {
125 for (index_t i0 = 0; i0 < local_NE0; i0++) {
126 // we will split this "rectangle" into two triangles
127 const dim_t k = 2 * (i0 + local_NE0 * i1);
128 const index_t node0 = Nstride0 * (i0 + e_offset0)
129 + Nstride1 * (i1 + e_offset1);
130
131 out->Elements->Id[k] = 2 * (i0 + e_offset0) + NE0*(i1 + e_offset1);
132 out->Elements->Tag[k] = 0;
133 out->Elements->Owner[k] = myRank;
134 out->Elements->Id[k + 1] = out->Elements->Id[k] + 1;
135 out->Elements->Tag[k + 1] = 0;
136 out->Elements->Owner[k + 1] = myRank;
137
138 // a,b,c,d gives the nodes of the rectangle in clockwise order
139 const index_t a = node0;
140 const index_t b = node0 + Nstride0;
141 const index_t c = node0 + Nstride1 + Nstride0;
142 const index_t d = node0 + Nstride1;
143 // For a little bit of variety
144 if ((global_adjustment + node0) % 2) {
145 out->Elements->Nodes[INDEX2(0, k, NN)] = a;
146 out->Elements->Nodes[INDEX2(1, k, NN)] = b;
147 out->Elements->Nodes[INDEX2(2, k, NN)] = d;
148 out->Elements->Nodes[INDEX2(0, k + 1, NN)] = b;
149 out->Elements->Nodes[INDEX2(1, k + 1, NN)] = c;
150 out->Elements->Nodes[INDEX2(2, k + 1, NN)] = d;
151 } else {
152 out->Elements->Nodes[INDEX2(0, k, NN)] = a;
153 out->Elements->Nodes[INDEX2(1, k, NN)] = b;
154 out->Elements->Nodes[INDEX2(2, k, NN)] = c;
155 out->Elements->Nodes[INDEX2(0, k + 1, NN)] = a;
156 out->Elements->Nodes[INDEX2(1, k + 1, NN)] = c;
157 out->Elements->Nodes[INDEX2(2, k + 1, NN)] = d;
158 }
159 }
160 }
161
162 // face elements
163 NN = out->FaceElements->numNodes;
164 dim_t totalNECount = 2 * NE0 * NE1; // because we have split the rectangles
165 dim_t faceNECount = 0;
166 if (local_NE0 > 0) {
167 // ** elements on boundary 001 (x1=0)
168 if (e_offset0 == 0) {
169 #pragma omp parallel for
170 for (index_t i1 = 0; i1 < local_NE1; i1++) {
171 const dim_t k = i1 + faceNECount;
172 const index_t node0 = Nstride1 * (i1 + e_offset1);
173 out->FaceElements->Id[k] = i1 + e_offset1 + totalNECount;
174 out->FaceElements->Tag[k] = LEFTTAG;
175 out->FaceElements->Owner[k] = myRank;
176 out->FaceElements->Nodes[INDEX2(0, k, NN)] = node0 + Nstride1;
177 out->FaceElements->Nodes[INDEX2(1, k, NN)] = node0;
178 }
179 faceNECount += local_NE1;
180 }
181 totalNECount += NE1;
182 // ** elements on boundary 002 (x1=1)
183 if (local_NE0 + e_offset0 == NE0) {
184 #pragma omp parallel for
185 for (index_t i1 = 0; i1 < local_NE1; i1++) {
186 const dim_t k = i1 + faceNECount;
187 const index_t node0 = Nstride0 * (NE0 - 1)
188 + Nstride1 * (i1 + e_offset1);
189
190 out->FaceElements->Id[k] = (i1 + e_offset1) + totalNECount;
191 out->FaceElements->Tag[k] = RIGHTTAG;
192 out->FaceElements->Owner[k] = myRank;
193 out->FaceElements->Nodes[INDEX2(0, k, NN)] = node0 + Nstride0;
194 out->FaceElements->Nodes[INDEX2(1, k, NN)] = node0 + Nstride1 + Nstride0;
195 }
196 faceNECount += local_NE1;
197 }
198 totalNECount += NE1;
199 }
200 if (local_NE1 > 0) {
201 // ** elements on boundary 010 (x2=0)
202 if (e_offset1 == 0) {
203 #pragma omp parallel for
204 for (index_t i0 = 0; i0 < local_NE0; i0++) {
205 const dim_t k = i0 + faceNECount;
206 const index_t node0 = Nstride0 * (i0 + e_offset0);
207 out->FaceElements->Id[k] = e_offset0 + i0 + totalNECount;
208 out->FaceElements->Tag[k] = BOTTOMTAG;
209 out->FaceElements->Owner[k] = myRank;
210 out->FaceElements->Nodes[INDEX2(0, k, NN)] = node0;
211 out->FaceElements->Nodes[INDEX2(1, k, NN)] = node0 + Nstride0;
212 }
213 faceNECount += local_NE0;
214 }
215 totalNECount += NE0;
216 // ** elements on boundary 020 (x2=1)
217 if (local_NE1 + e_offset1 == NE1) {
218 #pragma omp parallel for
219 for (index_t i0 = 0; i0 < local_NE0; i0++) {
220 const dim_t k = i0 + faceNECount;
221 const index_t node0 = Nstride0 * (i0 + e_offset0)
222 + Nstride1 * (NE1 - 1);
223
224 out->FaceElements->Id[k] = i0 + e_offset0 + totalNECount;
225 out->FaceElements->Tag[k] = TOPTAG;
226 out->FaceElements->Owner[k] = myRank;
227 out->FaceElements->Nodes[INDEX2(0, k, NN)] = node0 + Nstride1 + Nstride0;
228 out->FaceElements->Nodes[INDEX2(1, k, NN)] = node0 + Nstride1;
229 }
230 faceNECount += local_NE0;
231 }
232 totalNECount += NE0;
233 }
234
235 // add tag names
236 out->addTagMap("top", TOPTAG);
237 out->addTagMap("bottom", BOTTOMTAG);
238 out->addTagMap("left", LEFTTAG);
239 out->addTagMap("right", RIGHTTAG);
240
241 // prepare mesh for further calculations
242 out->resolveNodeIds();
243 out->prepare(optimize);
244 return out;
245 }
246
247 } // namespace dudley
248

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/4.0fordebian/dudley/src/Mesh_tri3.cpp:5567-5588 /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 /branches/symbolic_from_3470/dudley/src/Mesh_tri3.c:3471-3974 /branches/symbolic_from_3470/ripley/test/python/dudley/src/Mesh_tri3.c:3517-3974 /release/4.0/dudley/src/Mesh_tri3.cpp:5380-5406 /trunk/dudley/src/Mesh_tri3.cpp:4257-4344,5898-6007 /trunk/ripley/test/python/dudley/src/Mesh_tri3.c:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26