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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6939 - (show annotations)
Mon Jan 20 03:37:18 2020 UTC (4 months, 2 weeks ago) by uqaeller
File size: 10711 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::createRec4(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 optimize,
32 escript::JMPI mpiInfo)
33 {
34 const int N_PER_E = 1;
35 const int DIM = 2;
36 const int LEFTTAG = 1; // boundary x1=0
37 const int RIGHTTAG = 2; // boundary x1=1
38 const int BOTTOMTAG = 10; // boundary x2=0
39 const int TOPTAG = 20; // boundary x2=1
40 dim_t Nstride0 = 0, Nstride1 = 0, local_NE0, local_NE1;
41 index_t e_offset0 = 0, e_offset1 = 0;
42
43 const int myRank = mpiInfo->rank;
44
45 // set up the global dimensions of the mesh
46 NE0 = std::max((dim_t)1, NE0);
47 NE1 = std::max((dim_t)1, NE1);
48 const dim_t N0 = N_PER_E*NE0 + 1;
49 const dim_t N1 = N_PER_E*NE1 + 1;
50
51 // allocate mesh
52 std::stringstream name;
53 name << "Rectangular " << N0 << " x " << N1 << " mesh";
54 FinleyDomain* out = new FinleyDomain(name.str(), DIM, mpiInfo);
55
56 const_ReferenceElementSet_ptr refPoints, refContactElements, refFaceElements, refElements;
57 if (useElementsOnFace) {
58 refFaceElements.reset(new ReferenceElementSet(Rec4Face, order, reducedOrder));
59 refContactElements.reset(new ReferenceElementSet(Rec4Face_Contact, order, reducedOrder));
60 } else {
61 refFaceElements.reset(new ReferenceElementSet(Line2, order, reducedOrder));
62 refContactElements.reset(new ReferenceElementSet(Line2_Contact, order, reducedOrder));
63 }
64 refElements.reset(new ReferenceElementSet(Rec4, order, reducedOrder));
65 refPoints.reset(new ReferenceElementSet(Point1, order, reducedOrder));
66
67 ElementFile* elements = new ElementFile(refElements, mpiInfo);
68 out->setElements(elements);
69 ElementFile* faces = new ElementFile(refFaceElements, mpiInfo);
70 out->setFaceElements(faces);
71 out->setContactElements(new ElementFile(refContactElements, mpiInfo));
72 out->setPoints(new ElementFile(refPoints, mpiInfo));
73
74 // work out the largest dimension
75 if (N1 == std::max(N0, N1)) {
76 Nstride0 = 1;
77 Nstride1 = N0;
78 local_NE0 = NE0;
79 e_offset0 = 0;
80 mpiInfo->split(NE1, &local_NE1, &e_offset1);
81 } else {
82 Nstride0 = N1;
83 Nstride1 = 1;
84 mpiInfo->split(NE0, &local_NE0, &e_offset0);
85 local_NE1 = NE1;
86 e_offset1 = 0;
87 }
88 const index_t offset0 = e_offset0 * N_PER_E;
89 const index_t offset1 = e_offset1 * N_PER_E;
90 const dim_t local_N0 = local_NE0 > 0 ? local_NE0*N_PER_E+1 : 0;
91 const dim_t local_N1 = local_NE1 > 0 ? local_NE1*N_PER_E+1 : 0;
92 dim_t NDOF0 = 0, NDOF1 = 0;
93
94 // get the number of surface elements
95 dim_t NFaceElements = 0;
96 if (!periodic0 && local_NE0 > 0) {
97 NDOF0 = N0;
98 if (e_offset0 == 0)
99 NFaceElements += local_NE1;
100 if (local_NE0 + e_offset0 == NE0)
101 NFaceElements += local_NE1;
102 } else {
103 NDOF0 = N0 - 1;
104 }
105
106 if (!periodic1 && local_NE1 > 0) {
107 NDOF1 = N1;
108 if (e_offset1 == 0)
109 NFaceElements += local_NE0;
110 if (local_NE1 + e_offset1 == NE1)
111 NFaceElements += local_NE0;
112 } else {
113 NDOF1 = N1 - 1;
114 }
115
116 // allocate tables
117 NodeFile* nodes = out->getNodes();
118 nodes->allocTable(local_N0 * local_N1);
119 elements->allocTable(local_NE0 * local_NE1);
120 faces->allocTable(NFaceElements);
121
122 // create nodes
123 #pragma omp parallel for
124 for (index_t i1 = 0; i1 < local_N1; i1++) {
125 for (index_t i0 = 0; i0 < local_N0; i0++) {
126 const dim_t k = i0 + local_N0 * i1;
127 const index_t global_i0 = i0 + offset0;
128 const index_t global_i1 = i1 + offset1;
129 nodes->Coordinates[INDEX2(0, k, DIM)] = (real_t)global_i0 / (real_t)(N0 - 1) * l0;
130 nodes->Coordinates[INDEX2(1, k, DIM)] = (real_t)global_i1 / (real_t)(N1 - 1) * l1;
131 nodes->Id[k] = Nstride0 * global_i0 + Nstride1 * global_i1;
132 nodes->Tag[k] = 0;
133 nodes->globalDegreesOfFreedom[k] = Nstride0 * (global_i0 % NDOF0)
134 + Nstride1 * (global_i1 % NDOF1);
135 }
136 }
137
138 // set the elements
139 dim_t NN = elements->numNodes;
140 #pragma omp parallel for
141 for (index_t i1 = 0; i1 < local_NE1; i1++) {
142 for (index_t i0 = 0; i0 < local_NE0; i0++) {
143 const dim_t k = i0 + local_NE0 * i1;
144 const index_t node0 = Nstride0 * N_PER_E * (i0 + e_offset0)
145 + Nstride1 * N_PER_E * (i1 + e_offset1);
146
147 elements->Id[k] = (i0 + e_offset0) + NE0*(i1 + e_offset1);
148 elements->Tag[k] = 0;
149 elements->Owner[k] = myRank;
150
151 elements->Nodes[INDEX2(0, k, NN)] = node0;
152 elements->Nodes[INDEX2(1, k, NN)] = node0 + Nstride0;
153 elements->Nodes[INDEX2(2, k, NN)] = node0 + Nstride1 + Nstride0;
154 elements->Nodes[INDEX2(3, k, NN)] = node0 + Nstride1;
155 }
156 }
157
158 // face elements
159 NN = faces->numNodes;
160 dim_t totalNECount = NE0 * NE1;
161 dim_t faceNECount = 0;
162 index_t* eNodes = faces->Nodes;
163
164 if (!periodic0 && local_NE0 > 0) {
165 // ** elements on boundary 001 (x1=0)
166 if (e_offset0 == 0) {
167 #pragma omp parallel for
168 for (index_t i1 = 0; i1 < local_NE1; i1++) {
169 const dim_t k = i1 + faceNECount;
170 const index_t node0 = Nstride1 * N_PER_E * (i1 + e_offset1);
171
172 faces->Id[k] = i1 + e_offset1 + totalNECount;
173 faces->Tag[k] = LEFTTAG;
174 faces->Owner[k] = myRank;
175 if (useElementsOnFace) {
176 eNodes[INDEX2(0, k, NN)] = node0 + Nstride1;
177 eNodes[INDEX2(1, k, NN)] = node0;
178 eNodes[INDEX2(2, k, NN)] = node0 + Nstride0;
179 eNodes[INDEX2(3, k, NN)] = node0 + Nstride1 + Nstride0;
180 } else {
181 eNodes[INDEX2(0, k, NN)] = node0 + Nstride1;
182 eNodes[INDEX2(1, k, NN)] = node0;
183 }
184 }
185 faceNECount += local_NE1;
186 }
187 totalNECount += NE1;
188 // ** elements on boundary 002 (x1=1)
189 if (local_NE0 + e_offset0 == NE0) {
190 #pragma omp parallel for
191 for (index_t i1 = 0; i1 < local_NE1; i1++) {
192 const dim_t k = i1 + faceNECount;
193 const index_t node0 = Nstride0 * N_PER_E * (NE0 - 1)
194 + Nstride1 * N_PER_E * (i1 + e_offset1);
195
196 faces->Id[k] = (i1 + e_offset1) + totalNECount;
197 faces->Tag[k] = RIGHTTAG;
198 faces->Owner[k] = myRank;
199 if (useElementsOnFace) {
200 eNodes[INDEX2(0, k, NN)] = node0 + Nstride0;
201 eNodes[INDEX2(1, k, NN)] = node0 + Nstride1 + Nstride0;
202 eNodes[INDEX2(2, k, NN)] = node0 + Nstride1;
203 eNodes[INDEX2(3, k, NN)] = node0;
204 } else {
205 eNodes[INDEX2(0, k, NN)] = node0 + Nstride0;
206 eNodes[INDEX2(1, k, NN)] = node0 + Nstride1 + Nstride0;
207 }
208 }
209 faceNECount += local_NE1;
210 }
211 totalNECount += NE1;
212 }
213
214 if (!periodic1 && local_NE1 > 0) {
215 // ** elements on boundary 010 (x2=0)
216 if (e_offset1 == 0) {
217 #pragma omp parallel for
218 for (index_t i0 = 0; i0 < local_NE0; i0++) {
219 const dim_t k = i0 + faceNECount;
220 const index_t node0 = Nstride0 * N_PER_E * (i0 + e_offset0);
221 faces->Id[k] = e_offset0 + i0 + totalNECount;
222 faces->Tag[k] = BOTTOMTAG;
223 faces->Owner[k] = myRank;
224 if (useElementsOnFace) {
225 eNodes[INDEX2(0, k, NN)] = node0;
226 eNodes[INDEX2(1, k, NN)] = node0 + Nstride0;
227 eNodes[INDEX2(2, k, NN)] = node0 + Nstride1 + Nstride0;
228 eNodes[INDEX2(3, k, NN)] = node0 + Nstride1;
229 } else {
230 eNodes[INDEX2(0, k, NN)] = node0;
231 eNodes[INDEX2(1, k, NN)] = node0 + Nstride0;
232 }
233 }
234 faceNECount += local_NE0;
235 }
236 totalNECount += NE0;
237 // ** elements on boundary 020 (x2=1)
238 if (local_NE1 + e_offset1 == NE1) {
239 #pragma omp parallel for
240 for (index_t i0 = 0; i0 < local_NE0; i0++) {
241 const dim_t k = i0 + faceNECount;
242 const index_t node0 = Nstride0 * N_PER_E * (i0 + e_offset0)
243 + Nstride1 * N_PER_E * (NE1 - 1);
244
245 faces->Id[k] = i0 + e_offset0 + totalNECount;
246 faces->Tag[k] = TOPTAG;
247 faces->Owner[k] = myRank;
248 if (useElementsOnFace) {
249 eNodes[INDEX2(0, k, NN)] = node0 + Nstride1 + Nstride0;
250 eNodes[INDEX2(1, k, NN)] = node0 + Nstride1;
251 eNodes[INDEX2(2, k, NN)] = node0;
252 eNodes[INDEX2(3, k, NN)] = node0 + Nstride0;
253 } else {
254 eNodes[INDEX2(0, k, NN)] = node0 + Nstride1 + Nstride0;
255 eNodes[INDEX2(1, k, NN)] = node0 + Nstride1;
256 }
257 }
258 faceNECount += local_NE0;
259 }
260 totalNECount += NE0;
261 }
262
263 // add tag names
264 out->setTagMap("top", TOPTAG);
265 out->setTagMap("bottom", BOTTOMTAG);
266 out->setTagMap("left", LEFTTAG);
267 out->setTagMap("right", RIGHTTAG);
268
269 // prepare mesh for further calculations
270 out->resolveNodeIds();
271 out->prepare(optimize);
272 return out->getPtr();
273 }
274
275 } // namespace finley
276

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26