/[escript]/trunk/finley/src/Mesh_rec8.cpp
ViewVC logotype

Contents of /trunk/finley/src/Mesh_rec8.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6939 - (show annotations)
Mon Jan 20 03:37:18 2020 UTC (4 months, 1 week ago) by uqaeller
File size: 13857 byte(s)
Updated the copyright header.


1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2020 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
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-2017 by Centre for Geoscience Computing (GeoComp)
14 * Development from 2019 by School of Earth and Environmental Sciences
15 **
16 *****************************************************************************/
17
18 #include "FinleyDomain.h"
19
20 #include <escript/index.h>
21
22 using escript::DataTypes::real_t;
23
24 namespace finley {
25
26 escript::Domain_ptr FinleyDomain::createRec8(dim_t NE0, dim_t NE1,
27 real_t l0, real_t l1,
28 bool periodic0, bool periodic1,
29 int order, int reducedOrder,
30 bool useElementsOnFace,
31 bool useFullElementOrder,
32 bool useMacroElements,
33 bool optimize,
34 escript::JMPI mpiInfo)
35 {
36 const int N_PER_E = 2;
37 const int DIM = 2;
38 const int LEFTTAG = 1; // boundary x1=0
39 const int RIGHTTAG = 2; // boundary x1=1
40 const int BOTTOMTAG = 10; // boundary x2=0
41 const int TOPTAG = 20; // boundary x2=1
42 dim_t Nstride0 = 0, Nstride1 = 0, local_NE0, local_NE1;
43 index_t e_offset0 = 0, e_offset1 = 0;
44 const bool generateAllNodes = useFullElementOrder || useMacroElements;
45
46 const int myRank = mpiInfo->rank;
47
48 // set up the global dimensions of the mesh
49 NE0 = std::max((dim_t)1, NE0);
50 NE1 = std::max((dim_t)1, NE1);
51 const dim_t N0 = N_PER_E*NE0 + 1;
52 const dim_t N1 = N_PER_E*NE1 + 1;
53
54 // allocate mesh
55 std::stringstream name;
56 name << "Rectangular " << N0 << " x " << N1 << " mesh";
57 FinleyDomain* out = new FinleyDomain(name.str(), DIM, mpiInfo);
58
59 const_ReferenceElementSet_ptr refPoints, refContactElements, refFaceElements, refElements;
60 if (generateAllNodes) {
61 if (useMacroElements) {
62 refElements.reset(new ReferenceElementSet(Rec9Macro, order, reducedOrder));
63 } else {
64 refElements.reset(new ReferenceElementSet(Rec9, order, reducedOrder));
65 }
66 if (useElementsOnFace) {
67 throw escript::NotImplementedError("rich elements for Rec9 elements are not supported yet.");
68 } else {
69 if (useMacroElements) {
70 refFaceElements.reset(new ReferenceElementSet(Line3Macro, order, reducedOrder));
71 } else {
72 refFaceElements.reset(new ReferenceElementSet(Line3, order, reducedOrder));
73 }
74 refContactElements.reset(new ReferenceElementSet(Line3_Contact, order, reducedOrder));
75 }
76 } else { // !generateAllNodes
77 refElements.reset(new ReferenceElementSet(Rec8, order, reducedOrder));
78 if (useElementsOnFace) {
79 refFaceElements.reset(new ReferenceElementSet(Rec8Face, order, reducedOrder));
80 refContactElements.reset(new ReferenceElementSet(Rec8Face_Contact, order, reducedOrder));
81 } else {
82 refFaceElements.reset(new ReferenceElementSet(Line3, order, reducedOrder));
83 refContactElements.reset(new ReferenceElementSet(Line3_Contact, order, reducedOrder));
84 }
85 }
86 refPoints.reset(new ReferenceElementSet(Point1, order, reducedOrder));
87 if (!refPoints->referenceElement)
88 throw FinleyException("ERRRRORRRRR!!");
89
90 ElementFile* elements = new ElementFile(refElements, mpiInfo);
91 out->setElements(elements);
92 ElementFile* faces = new ElementFile(refFaceElements, mpiInfo);
93 out->setFaceElements(faces);
94 out->setContactElements(new ElementFile(refContactElements, mpiInfo));
95 out->setPoints(new ElementFile(refPoints, mpiInfo));
96
97 // work out the largest dimension
98 if (N1 == std::max(N0, N1)) {
99 Nstride0 = 1;
100 Nstride1 = N0;
101 local_NE0 = NE0;
102 e_offset0 = 0;
103 mpiInfo->split(NE1, &local_NE1, &e_offset1);
104 } else {
105 Nstride0 = N1;
106 Nstride1 = 1;
107 mpiInfo->split(NE0, &local_NE0, &e_offset0);
108 local_NE1 = NE1;
109 e_offset1 = 0;
110 }
111 const index_t offset0 = e_offset0 * N_PER_E;
112 const index_t offset1 = e_offset1 * N_PER_E;
113 const dim_t local_N0 = local_NE0 > 0 ? local_NE0*N_PER_E+1 : 0;
114 const dim_t local_N1 = local_NE1 > 0 ? local_NE1*N_PER_E+1 : 0;
115 dim_t NDOF0 = 0, NDOF1 = 0;
116
117 // get the number of surface elements
118 dim_t NFaceElements = 0;
119 if (!periodic0 && local_NE0 > 0) {
120 NDOF0 = N0;
121 if (e_offset0 == 0)
122 NFaceElements += local_NE1;
123 if (local_NE0 + e_offset0 == NE0)
124 NFaceElements += local_NE1;
125 } else {
126 NDOF0 = N0 - 1;
127 }
128
129 if (!periodic1 && local_NE1 > 0) {
130 NDOF1 = N1;
131 if (e_offset1 == 0)
132 NFaceElements += local_NE0;
133 if (local_NE1 + e_offset1 == NE1)
134 NFaceElements += local_NE0;
135 } else {
136 NDOF1 = N1 - 1;
137 }
138
139 // allocate tables
140 NodeFile* nodes = out->getNodes();
141 nodes->allocTable(local_N0 * local_N1);
142 elements->allocTable(local_NE0 * local_NE1);
143 faces->allocTable(NFaceElements);
144
145 // create nodes
146 #pragma omp parallel for
147 for (index_t i1 = 0; i1 < local_N1; i1++) {
148 for (index_t i0 = 0; i0 < local_N0; i0++) {
149 const dim_t k = i0 + local_N0 * i1;
150 const index_t global_i0 = i0 + offset0;
151 const index_t global_i1 = i1 + offset1;
152 nodes->Coordinates[INDEX2(0, k, DIM)] = (real_t)global_i0 / (real_t)(N0 - 1) * l0;
153 nodes->Coordinates[INDEX2(1, k, DIM)] = (real_t)global_i1 / (real_t)(N1 - 1) * l1;
154 nodes->Id[k] = Nstride0 * global_i0 + Nstride1 * global_i1;
155 nodes->Tag[k] = 0;
156 nodes->globalDegreesOfFreedom[k] = Nstride0 * (global_i0 % NDOF0)
157 + Nstride1 * (global_i1 % NDOF1);
158 }
159 }
160
161 // set the elements
162 dim_t NN = elements->numNodes;
163 #pragma omp parallel for
164 for (index_t i1 = 0; i1 < local_NE1; i1++) {
165 for (index_t i0 = 0; i0 < local_NE0; i0++) {
166 const dim_t k = i0 + local_NE0 * i1;
167 const index_t node0 = Nstride0 * N_PER_E * (i0 + e_offset0)
168 + Nstride1 * N_PER_E * (i1 + e_offset1);
169
170 elements->Id[k] = (i0 + e_offset0) + NE0*(i1 + e_offset1);
171 elements->Tag[k] = 0;
172 elements->Owner[k] = myRank;
173
174 elements->Nodes[INDEX2(0, k, NN)] = node0;
175 elements->Nodes[INDEX2(1, k, NN)] = node0 + 2 * Nstride0;
176 elements->Nodes[INDEX2(2, k, NN)] = node0 + 2 * Nstride1 + 2 * Nstride0;
177 elements->Nodes[INDEX2(3, k, NN)] = node0 + 2 * Nstride1;
178 elements->Nodes[INDEX2(4, k, NN)] = node0 + 1 * Nstride0;
179 elements->Nodes[INDEX2(5, k, NN)] = node0 + Nstride1 + 2 * Nstride0;
180 elements->Nodes[INDEX2(6, k, NN)] = node0 + 2 * Nstride1 + Nstride0;
181 elements->Nodes[INDEX2(7, k, NN)] = node0 + Nstride1;
182 if (generateAllNodes) {
183 elements->Nodes[INDEX2(8, k, NN)] = node0 + Nstride1 + Nstride0;
184 }
185 }
186 }
187
188 // face elements
189 NN = faces->numNodes;
190 dim_t totalNECount = NE0 * NE1;
191 dim_t faceNECount = 0;
192 index_t* eNodes = faces->Nodes;
193
194 if (!periodic0 && local_NE0 > 0) {
195 // ** elements on boundary 001 (x1=0)
196 if (e_offset0 == 0) {
197 #pragma omp parallel for
198 for (index_t i1 = 0; i1 < local_NE1; i1++) {
199 const dim_t k = i1 + faceNECount;
200 const index_t node0 = Nstride1 * N_PER_E * (i1 + e_offset1);
201
202 faces->Id[k] = i1 + e_offset1 + totalNECount;
203 faces->Tag[k] = LEFTTAG;
204 faces->Owner[k] = myRank;
205 if (useElementsOnFace) {
206 eNodes[INDEX2(0, k, NN)] = node0 + 2 * Nstride1;
207 eNodes[INDEX2(1, k, NN)] = node0;
208 eNodes[INDEX2(2, k, NN)] = node0 + 2 * Nstride0;
209 eNodes[INDEX2(3, k, NN)] = node0 + 2 * Nstride1 + 2 * Nstride0;
210 eNodes[INDEX2(4, k, NN)] = node0 + Nstride1;
211 eNodes[INDEX2(5, k, NN)] = node0 + Nstride0;
212 eNodes[INDEX2(6, k, NN)] = node0 + Nstride1 + 2 * Nstride0;
213 eNodes[INDEX2(7, k, NN)] = node0 + 2 * Nstride1 + Nstride0;
214 } else {
215 eNodes[INDEX2(0, k, NN)] = node0 + 2 * Nstride1;
216 eNodes[INDEX2(1, k, NN)] = node0;
217 eNodes[INDEX2(2, k, NN)] = node0 + Nstride1;
218 }
219 }
220 faceNECount += local_NE1;
221 }
222 totalNECount += NE1;
223 // ** elements on boundary 002 (x1=1)
224 if (local_NE0 + e_offset0 == NE0) {
225 #pragma omp parallel for
226 for (index_t i1 = 0; i1 < local_NE1; i1++) {
227 const dim_t k = i1 + faceNECount;
228 const index_t node0 = Nstride0 * N_PER_E * (NE0 - 1)
229 + Nstride1 * N_PER_E * (i1 + e_offset1);
230
231 faces->Id[k] = (i1 + e_offset1) + totalNECount;
232 faces->Tag[k] = RIGHTTAG;
233 faces->Owner[k] = myRank;
234 if (useElementsOnFace) {
235 eNodes[INDEX2(0, k, NN)] = node0 + 2 * Nstride0;
236 eNodes[INDEX2(1, k, NN)] = node0 + 2 * Nstride1 + 2 * Nstride0;
237 eNodes[INDEX2(2, k, NN)] = node0 + 2 * Nstride1;
238 eNodes[INDEX2(3, k, NN)] = node0;
239 eNodes[INDEX2(4, k, NN)] = node0 + Nstride1 + 2 * Nstride0;
240 eNodes[INDEX2(5, k, NN)] = node0 + 2 * Nstride1 + Nstride0;
241 eNodes[INDEX2(6, k, NN)] = node0 + Nstride1;
242 eNodes[INDEX2(7, k, NN)] = node0 + Nstride0;
243 } else {
244 eNodes[INDEX2(0, k, NN)] = node0 + 2 * Nstride0;
245 eNodes[INDEX2(1, k, NN)] = node0 + 2 * Nstride1 + 2 * Nstride0;
246 eNodes[INDEX2(2, k, NN)] = node0 + Nstride1 + 2 * Nstride0;
247 }
248 }
249 faceNECount += local_NE1;
250 }
251 totalNECount += NE1;
252 }
253
254 if (!periodic1 && local_NE1 > 0) {
255 // ** elements on boundary 010 (x2=0)
256 if (e_offset1 == 0) {
257 #pragma omp parallel for
258 for (index_t i0 = 0; i0 < local_NE0; i0++) {
259 const dim_t k = i0 + faceNECount;
260 const index_t node0 = Nstride0 * N_PER_E * (i0 + e_offset0);
261 faces->Id[k] = e_offset0 + i0 + totalNECount;
262 faces->Tag[k] = BOTTOMTAG;
263 faces->Owner[k] = myRank;
264 if (useElementsOnFace) {
265 eNodes[INDEX2(0, k, NN)] = node0;
266 eNodes[INDEX2(1, k, NN)] = node0 + 2 * Nstride0;
267 eNodes[INDEX2(2, k, NN)] = node0 + 2 * Nstride1 + 2 * Nstride0;
268 eNodes[INDEX2(3, k, NN)] = node0 + 2 * Nstride1;
269 eNodes[INDEX2(4, k, NN)] = node0 + Nstride0;
270 eNodes[INDEX2(5, k, NN)] = node0 + Nstride1 + 2 * Nstride0;
271 eNodes[INDEX2(6, k, NN)] = node0 + 2 * Nstride1 + Nstride0;
272 eNodes[INDEX2(7, k, NN)] = node0 + Nstride1;
273 } else {
274 eNodes[INDEX2(0, k, NN)] = node0;
275 eNodes[INDEX2(1, k, NN)] = node0 + 2 * Nstride0;
276 eNodes[INDEX2(2, k, NN)] = node0 + Nstride0;
277 }
278 }
279 faceNECount += local_NE0;
280 }
281 totalNECount += NE0;
282 // ** elements on boundary 020 (x2=1)
283 if (local_NE1 + e_offset1 == NE1) {
284 #pragma omp parallel for
285 for (index_t i0 = 0; i0 < local_NE0; i0++) {
286 const dim_t k = i0 + faceNECount;
287 const index_t node0 = Nstride0 * N_PER_E * (i0 + e_offset0)
288 + Nstride1 * N_PER_E * (NE1 - 1);
289
290 faces->Id[k] = i0 + e_offset0 + totalNECount;
291 faces->Tag[k] = TOPTAG;
292 faces->Owner[k] = myRank;
293 if (useElementsOnFace) {
294 eNodes[INDEX2(0, k, NN)] = node0 + 2 * Nstride1 + 2 * Nstride0;
295 eNodes[INDEX2(1, k, NN)] = node0 + 2 * Nstride1;
296 eNodes[INDEX2(2, k, NN)] = node0;
297 eNodes[INDEX2(3, k, NN)] = node0 + 2 * Nstride0;
298 eNodes[INDEX2(4, k, NN)] = node0 + 2 * Nstride1 + Nstride0;
299 eNodes[INDEX2(5, k, NN)] = node0 + Nstride1;
300 eNodes[INDEX2(6, k, NN)] = node0 + Nstride0;
301 eNodes[INDEX2(7, k, NN)] = node0 + Nstride1 + 2 * Nstride0;
302 } else {
303 eNodes[INDEX2(0, k, NN)] = node0 + 2 * Nstride1 + 2 * Nstride0;
304 eNodes[INDEX2(1, k, NN)] = node0 + 2 * Nstride1;
305 eNodes[INDEX2(2, k, NN)] = node0 + 2 * Nstride1 + Nstride0;
306 }
307 }
308 faceNECount += local_NE0;
309 }
310 totalNECount += NE0;
311 }
312
313 // add tag names
314 out->setTagMap("top", TOPTAG);
315 out->setTagMap("bottom", BOTTOMTAG);
316 out->setTagMap("left", LEFTTAG);
317 out->setTagMap("right", RIGHTTAG);
318
319 // prepare mesh for further calculations
320 out->resolveNodeIds();
321 out->prepare(optimize);
322 return out->getPtr();
323 }
324
325 } // namespace finley
326

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/4.0fordebian/finley/src/Mesh_rec8.cpp:5567-5588 /branches/lapack2681/finley/src/Mesh_rec8.cpp:2682-2741 /branches/pasowrap/finley/src/Mesh_rec8.cpp:3661-3674 /branches/py3_attempt2/finley/src/Mesh_rec8.cpp:3871-3891 /branches/restext/finley/src/Mesh_rec8.cpp:2610-2624 /branches/ripleygmg_from_3668/finley/src/Mesh_rec8.cpp:3669-3791 /branches/stage3.0/finley/src/Mesh_rec8.cpp:2569-2590 /branches/symbolic_from_3470/finley/src/Mesh_rec8.cpp:3471-3974 /branches/trilinos_from_5897/finley/src/Mesh_rec8.cpp:5898-6118 /release/3.0/finley/src/Mesh_rec8.cpp:2591-2601 /release/4.0/finley/src/Mesh_rec8.cpp:5380-5406 /trunk/finley/src/Mesh_rec8.cpp:4257-4344

  ViewVC Help
Powered by ViewVC 1.1.26