PLaSK library
Loading...
Searching...
No Matches
generator_rectangular.cpp
Go to the documentation of this file.
1/*
2 * This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
3 * Copyright (c) 2022 Lodz University of Technology
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, version 3.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 */
14#include "plask/log/log.hpp"
15#include "plask/manager.hpp"
16
17#include "generator_rectangular.hpp"
18
19namespace plask {
20
21template <int dim>
22inline static void addPoints(shared_ptr<OrderedAxis> mesh,
23 const shared_ptr<GeometryObjectD<dim>>& geometry,
24 Primitive<3>::Direction dir,
25 double split) {
26 OrderedAxis::WarningOff warning_off(mesh);
27 (void)warning_off; // don't warn about unused variable warning_off
28 std::set<double> pts = geometry->getPointsAlong(dir);
29 if (split == 0) {
30 mesh->addOrderedPoints(pts.begin(), pts.end(), pts.size());
31 } else {
32 std::vector<double> ptsv;
33 ptsv.reserve(2 * pts.size());
34 for (double p : pts) {
35 std::vector<double>::reverse_iterator it;
36 for (it = ptsv.rbegin(); it != ptsv.rend() && *it > p - split; ++it)
37 ;
38 ptsv.insert(it.base(), p - split);
39 ptsv.push_back(p + split);
40 }
41 mesh->addOrderedPoints(ptsv.begin(), ptsv.end(), ptsv.size());
42 }
43}
44
45shared_ptr<OrderedAxis> makeGeometryGrid1D(const shared_ptr<GeometryObjectD<2>>& geometry, double split) {
46 auto mesh = plask::make_shared<OrderedAxis>();
47 addPoints(mesh, geometry, Primitive<3>::DIRECTION_TRAN, split);
48 return mesh;
49}
50
51shared_ptr<MeshD<1>> OrderedMesh1DSimpleGenerator::generate(const shared_ptr<GeometryObjectD<2>>& geometry) {
52 auto mesh = makeGeometryGrid1D(geometry, split ? OrderedAxis::MIN_DISTANCE : 0.);
53 writelog(LOG_DETAIL, "mesh.Rectangular1D.SimpleGenerator: Generating new mesh ({0})", mesh->size());
54 return mesh;
55}
56
57shared_ptr<RectangularMesh<2>> makeGeometryGrid(const shared_ptr<GeometryObjectD<2>>& geometry, double split) {
58 shared_ptr<OrderedAxis> axis0(new OrderedAxis), axis1(new OrderedAxis);
59 addPoints(axis0, geometry, Primitive<3>::DIRECTION_TRAN, split);
60 addPoints(axis1, geometry, Primitive<3>::DIRECTION_VERT, split);
61 shared_ptr<RectangularMesh<2>> mesh = plask::make_shared<RectangularMesh<2>>(std::move(axis0), std::move(axis1));
62 mesh->setOptimalIterationOrder();
63 return mesh;
64}
65
66shared_ptr<MeshD<2>> RectangularMesh2DSimpleGenerator::generate(const shared_ptr<GeometryObjectD<2>>& geometry) {
67 shared_ptr<RectangularMesh<2>> mesh = makeGeometryGrid(geometry, split ? OrderedAxis::MIN_DISTANCE : 0.);
68 writelog(LOG_DETAIL, "mesh.Rectangular2D.SimpleGenerator: Generating new mesh ({0}x{1})", mesh->axis[0]->size(),
69 mesh->axis[1]->size());
70 return mesh;
71}
72
73shared_ptr<MeshD<2>> RectangularMesh2DFrom1DGenerator::generate(const shared_ptr<GeometryObjectD<2>>& geometry) {
74 return plask::make_shared<RectangularMesh<2>>(horizontal_generator->get<MeshAxis>(geometry),
75 makeGeometryGrid(geometry)->axis[1]);
76}
77
78shared_ptr<RectangularMesh<3>> makeGeometryGrid(const shared_ptr<GeometryObjectD<3>>& geometry, double split) {
79 shared_ptr<OrderedAxis> axis0(new OrderedAxis), axis1(new OrderedAxis), axis2(new OrderedAxis);
80 addPoints(axis0, geometry, Primitive<3>::DIRECTION_LONG, split);
81 addPoints(axis1, geometry, Primitive<3>::DIRECTION_TRAN, split);
82 addPoints(axis2, geometry, Primitive<3>::DIRECTION_VERT, split);
83 shared_ptr<RectangularMesh<3>> mesh =
84 plask::make_shared<RectangularMesh<3>>(std::move(axis0), std::move(axis1), std::move(axis2));
85 mesh->setOptimalIterationOrder();
86 return mesh;
87}
88
89shared_ptr<MeshD<3>> RectangularMesh3DSimpleGenerator::generate(const shared_ptr<GeometryObjectD<3>>& geometry) {
90 auto mesh = makeGeometryGrid(geometry, split ? OrderedAxis::MIN_DISTANCE : 0.);
91 writelog(LOG_DETAIL, "mesh.Rectangular3D.SimpleGenerator: Generating new mesh ({0}x{1}x{2})", mesh->axis[0]->size(),
92 mesh->axis[1]->size(), mesh->axis[2]->size());
93 return mesh;
94}
95
96shared_ptr<OrderedAxis> refineAxis(const shared_ptr<MeshAxis>& axis, double spacing) {
97 if (spacing == 0. || isinf(spacing) || isnan(spacing)) return make_shared<OrderedAxis>(*axis);
98 size_t total = 1;
99 for (size_t i = 1; i < axis->size(); ++i) {
100 total += size_t(max(round((axis->at(i) - axis->at(i - 1)) / spacing), 1.));
101 }
102 std::vector<double> points;
103 points.reserve(total);
104 for (size_t i = 1; i < axis->size(); ++i) {
105 const double offset = axis->at(i - 1);
106 const double range = axis->at(i) - offset;
107 const double steps = std::max(round(range / spacing), 1.);
108 const double step = range / steps;
109 for (size_t j = 0, n = std::size_t(steps); j < n; ++j) {
110 points.push_back(offset + double(j) * step);
111 }
112 }
113 points.push_back(axis->at(axis->size() - 1));
114 assert(points.size() == total);
115 return shared_ptr<OrderedAxis>(new OrderedAxis(std::move(points)));
116}
117
118shared_ptr<MeshD<1>> OrderedMesh1DRegularGenerator::generate(const shared_ptr<GeometryObjectD<2>>& geometry) {
119 auto mesh = refineAxis(makeGeometryGrid1D(geometry, split ? OrderedAxis::MIN_DISTANCE : 0.), spacing);
120 writelog(LOG_DETAIL, "mesh.Rectangular1D.RegularGenerator: Generating new mesh ({0})", mesh->size());
121 return mesh;
122}
123
124shared_ptr<MeshD<2>> RectangularMesh2DRegularGenerator::generate(const shared_ptr<GeometryObjectD<2>>& geometry) {
125 auto mesh1 = makeGeometryGrid(geometry, split ? OrderedAxis::MIN_DISTANCE : 0.);
126 auto mesh = make_shared<RectangularMesh<2>>(refineAxis(mesh1->axis[0], spacing0), refineAxis(mesh1->axis[1], spacing1));
127 writelog(LOG_DETAIL, "mesh.Rectangular2D.RegularGenerator: Generating new mesh ({0}x{1})", mesh->axis[0]->size(),
128 mesh->axis[1]->size());
129 return mesh;
130}
131
132shared_ptr<MeshD<3>> RectangularMesh3DRegularGenerator::generate(const shared_ptr<GeometryObjectD<3>>& geometry) {
133 auto mesh1 = makeGeometryGrid(geometry, split ? OrderedAxis::MIN_DISTANCE : 0.);
134 auto mesh = make_shared<RectangularMesh<3>>(refineAxis(mesh1->axis[0], spacing0), refineAxis(mesh1->axis[1], spacing1),
135 refineAxis(mesh1->axis[2], spacing2));
136 writelog(LOG_DETAIL, "mesh.Rectangular3D.RegularGenerator: Generating new mesh ({0}x{1}x{2})", mesh->axis[0]->size(),
137 mesh->axis[1]->size(), mesh->axis[2]->size());
138 return mesh;
139}
140
141template <int dim> std::pair<double, double> RectangularMeshRefinedGenerator<dim>::getMinMax(const shared_ptr<OrderedAxis>& axis) {
142 double min = INFINITY, max = 0;
143 for (size_t i = 1; i != axis->size(); ++i) {
144 double L = axis->at(i) - axis->at(i - 1);
145 if (L < min) min = L;
146 if (L > max) max = L;
147 }
148 return std::pair<double, double>(min, max);
149}
150
151template <int dim> void RectangularMeshRefinedGenerator<dim>::divideLargestSegment(shared_ptr<OrderedAxis> axis) {
152 double max = 0;
153 double newpoint;
154 for (size_t i = 1; i != axis->size(); ++i) {
155 double L = axis->at(i) - axis->at(i - 1);
156 if (L > max) {
157 max = L;
158 newpoint = 0.5 * (axis->at(i - 1) + axis->at(i));
159 }
160 }
161 OrderedAxis::WarningOff warning_off(axis);
162 axis->addPoint(newpoint);
163}
164
165template <int dim>
166shared_ptr<OrderedAxis> RectangularMeshRefinedGenerator<dim>::getAxis(shared_ptr<OrderedAxis> axis,
167 const shared_ptr<GeometryObjectD<DIM>>& geometry,
168 size_t dir) {
169 assert(bool(axis));
170 OrderedAxis::WarningOff warning_off(axis);
171
172 double geometry_lower = geometry->getBoundingBox().lower[dir], geometry_upper = geometry->getBoundingBox().upper[dir];
173
174 // Add refinement points
175 for (auto ref : this->refinements[dir]) {
176 auto object = ref.first.first.lock();
177 if (!object) {
178 writelog(LOG_WARNING, "{}: Refinement defined for object not existing any more", name());
179 } else {
180 auto path = ref.first.second;
181 auto boxes = geometry->getObjectBoundingBoxes(*object, path);
182 auto origins = geometry->getObjectPositions(*object, path);
183 if (boxes.size() == 0)
184 writelog(LOG_WARNING, "DivideGenerator: Refinement defined for object absent from the geometry");
185 auto box = boxes.begin();
186 auto origin = origins.begin();
187 for (; box != boxes.end(); ++box, ++origin) {
188 for (auto x : ref.second) {
189 double zero = (*origin)[dir];
190 double lower = box->lower[dir] - zero;
191 double upper = box->upper[dir] - zero;
192 double x0 = x + zero;
193 if (geometry_lower <= x0 && x0 <= geometry_upper)
194 axis->addPoint(x0);
195 }
196 }
197 }
198 }
199
200 // Have specialization make further axis processing
201 return processAxis(axis, geometry, dir);
202}
203
204template <>
205shared_ptr<MeshD<1>> RectangularMeshRefinedGenerator<1>::generate(const boost::shared_ptr<plask::GeometryObjectD<2>>& geometry) {
206 shared_ptr<OrderedAxis> mesh = makeGeometryGrid1D(geometry);
207 getAxis(mesh, geometry, 0);
208 writelog(LOG_DETAIL, "mesh.Rectilinear1D::{}: Generating new mesh ({:d})", name(), mesh->size());
209 return mesh;
210}
211
212template <>
213shared_ptr<MeshD<2>> RectangularMeshRefinedGenerator<2>::generate(const boost::shared_ptr<plask::GeometryObjectD<2>>& geometry) {
214 auto mesh = makeGeometryGrid(geometry);
215 auto axis0 = dynamic_pointer_cast<OrderedAxis>(mesh->axis[0]), axis1 = dynamic_pointer_cast<OrderedAxis>(mesh->axis[1]);
216 getAxis(axis0, geometry, 0);
217 getAxis(axis1, geometry, 1);
218
219 auto minmax0 = getMinMax(axis0), minmax1 = getMinMax(axis1);
220 double asp0 = minmax0.second / minmax1.first, asp1 = minmax1.second / minmax0.first;
221 double limit = (1 + SMALL) * aspect;
222 while (aspect != 0. && (asp0 > limit || asp1 > limit)) {
223 if (asp0 > aspect) divideLargestSegment(axis0);
224 if (asp1 > aspect) divideLargestSegment(axis1);
225 minmax0 = getMinMax(axis0);
226 minmax1 = getMinMax(axis1);
227 asp0 = minmax0.second / minmax1.first;
228 asp1 = minmax1.second / minmax0.first;
229 }
230
231 mesh->setOptimalIterationOrder();
232 writelog(LOG_DETAIL, "mesh.Rectangular2D::{}: Generating new mesh ({:d}x{:d}, max. aspect {:.0f}:1)", name(),
233 mesh->axis[0]->size(), mesh->axis[1]->size(), max(asp0, asp1));
234 return mesh;
235}
236
237template <>
238shared_ptr<MeshD<3>> RectangularMeshRefinedGenerator<3>::generate(const boost::shared_ptr<plask::GeometryObjectD<3>>& geometry) {
239 auto mesh = makeGeometryGrid(geometry);
240 auto axis0 = dynamic_pointer_cast<OrderedAxis>(mesh->axis[0]), axis1 = dynamic_pointer_cast<OrderedAxis>(mesh->axis[1]),
241 axis2 = dynamic_pointer_cast<OrderedAxis>(mesh->axis[2]);
242 getAxis(axis0, geometry, 0);
243 getAxis(axis1, geometry, 1);
244 getAxis(axis2, geometry, 2);
245
246 auto minmax0 = getMinMax(axis0), minmax1 = getMinMax(axis1), minmax2 = getMinMax(axis2);
247 double asp0 = minmax0.second / min(minmax1.first, minmax2.first), asp1 = minmax1.second / min(minmax0.first, minmax2.first),
248 asp2 = minmax2.second / min(minmax0.first, minmax1.first);
249 double limit = (1 + SMALL) * aspect;
250 while (aspect != 0. && (asp0 > limit || asp1 > limit || asp2 > limit)) {
251 if (asp0 > aspect) divideLargestSegment(axis0);
252 if (asp1 > aspect) divideLargestSegment(axis1);
253 if (asp2 > aspect) divideLargestSegment(axis2);
254 minmax0 = getMinMax(axis0);
255 minmax1 = getMinMax(axis1);
256 minmax2 = getMinMax(axis2);
257 asp0 = minmax0.second / min(minmax1.first, minmax2.first);
258 asp1 = minmax1.second / min(minmax0.first, minmax2.first);
259 asp2 = minmax2.second / min(minmax0.first, minmax1.first);
260 }
261
262 mesh->setOptimalIterationOrder();
263 writelog(LOG_DETAIL, "mesh.Rectangular3D::{}: Generating new mesh ({:d}x{:d}x{:d}, max. aspect {:.0f}:1)", name(),
264 mesh->axis[0]->size(), mesh->axis[1]->size(), mesh->axis[2]->size(), max(asp0, max(asp1, asp2)));
265 return mesh;
266}
267
268template <int dim>
269shared_ptr<OrderedAxis> RectangularMeshDivideGenerator<dim>::processAxis(shared_ptr<OrderedAxis> axis,
270 const shared_ptr<GeometryObjectD<DIM>>& /*geometry*/,
271 size_t dir) {
272 assert(bool(axis));
273 OrderedAxis::WarningOff warning_off(axis);
274 (void)warning_off; // don't warn about unused variable warning_off
275
276 if (pre_divisions[dir] == 0) pre_divisions[dir] = 1;
277 if (post_divisions[dir] == 0) post_divisions[dir] = 1;
278
279 OrderedAxis& result = *axis.get();
280
281 // Pre-divide each object
282 double x = *result.begin();
283 std::vector<double> points;
284 points.reserve((pre_divisions[dir] - 1) * (result.size() - 1));
285 for (auto i = result.begin() + 1; i != result.end(); ++i) {
286 double w = *i - x;
287 for (size_t j = 1; j != pre_divisions[dir]; ++j) points.push_back(x + w * double(j) / double(pre_divisions[dir]));
288 x = *i;
289 }
290 result.addOrderedPoints(points.begin(), points.end());
291
292 // Now ensure, that the grids do not change to quickly
293 if (result.size() > 2 && getGradual(dir)) {
294 size_t end = result.size() - 2;
295 double w_prev = INFINITY, w = result[1] - result[0], w_next = result[2] - result[1];
296 for (size_t i = 0; i <= end;) {
297 bool goon = true;
298 if (w > 2.001 * w_prev) { // .0001 is for border case w == 2*w_prev, to avoid division even in presence of
299 // numerical error
300 if (result.addPoint(0.5 * (result[i] + result[i + 1]))) {
301 ++end;
302 w = w_next = result[i + 1] - result[i];
303 goon = false;
304 }
305 } else if (w > 2.001 * w_next) {
306 if (result.addPoint(0.5 * (result[i] + result[i + 1]))) {
307 ++end;
308 w_next = result[i + 1] - result[i];
309 if (i) {
310 --i;
311 w = w_prev;
312 w_prev = (i == 0) ? INFINITY : result[i] - result[i - 1];
313 } else
314 w = w_next;
315 goon = false;
316 }
317 }
318 if (goon) {
319 ++i;
320 w_prev = w;
321 w = w_next;
322 w_next = (i < end) ? result[i + 2] - result[i + 1] : INFINITY;
323 }
324 }
325 }
326
327 // Finally divide each object in post- division
328 x = *result.begin();
329 points.clear();
330 points.reserve((post_divisions[dir] - 1) * (result.size() - 1));
331 for (auto i = result.begin() + 1; i != result.end(); ++i) {
332 double w = *i - x;
333 for (size_t j = 1; j != post_divisions[dir]; ++j) points.push_back(x + w * double(j) / double(post_divisions[dir]));
334 x = *i;
335 }
336
337 result.addOrderedPoints(points.begin(), points.end());
338
339 return axis;
340}
341
342template <> RectangularMeshSmoothGenerator<1>::RectangularMeshSmoothGenerator() : finestep{0.005}, maxstep{INFINITY}, factor{1.2} {}
343
344template <>
345RectangularMeshSmoothGenerator<2>::RectangularMeshSmoothGenerator()
346 : finestep{0.005, 0.005}, maxstep{INFINITY, INFINITY}, factor{1.2, 1.2} {}
347
348template <>
349RectangularMeshSmoothGenerator<3>::RectangularMeshSmoothGenerator()
350 : finestep{0.005, 0.005, 0.005}, maxstep{INFINITY, INFINITY, INFINITY}, factor{1.2, 1.2, 1.2} {}
351
352template <int dim>
353shared_ptr<OrderedAxis> RectangularMeshSmoothGenerator<dim>::processAxis(shared_ptr<OrderedAxis> axis,
354 const shared_ptr<GeometryObjectD<DIM>>& /*geometry*/,
355 size_t dir) {
356 OrderedAxis::WarningOff warning_off(axis);
357
358 // Next divide each object
359 double x = *axis->begin();
360 std::vector<double> points; // points.reserve(...);
361 for (auto i = axis->begin() + 1; i != axis->end(); ++i) {
362 double width = *i - x;
363 if (width + OrderedAxis::MIN_DISTANCE <= finestep[dir]) continue;
364 if (factor[dir] == 1.) {
365 double m = ceil(width / finestep[dir]);
366 double d = width / m;
367 for (size_t i = 1, n = size_t(m); i < n; ++i) points.push_back(x + double(i) * d);
368 continue;
369 }
370 double logf = log(factor[dir]);
371 double maxm = floor(log(maxstep[dir] / finestep[dir]) / logf + OrderedAxis::MIN_DISTANCE);
372 double m = ceil(log(0.5 * (width - OrderedAxis::MIN_DISTANCE) / finestep[dir] * (factor[dir] - 1.) + 1.) / logf) -
373 1.; // number of points in one half
374 size_t lin = 0;
375 if (m > maxm) {
376 m = maxm;
377 lin = 1;
378 }
379 size_t n = size_t(m);
380 double end = finestep[dir] * (pow(factor[dir], m) - 1.) / (factor[dir] - 1.);
381 double last = finestep[dir] * pow(factor[dir], m);
382 if (lin) {
383 lin = size_t(ceil((width - 2. * end) / last));
384 } else if (width - 2. * end <= last) {
385 lin = 1;
386 } else {
387 lin = 2;
388 }
389 double s = finestep[dir] * 0.5 * width / (end + 0.5 * double(lin) * last);
390 double dx = 0.;
391 for (size_t i = 0; i < n; ++i) {
392 dx += s;
393 s *= factor[dir];
394 points.push_back(x + dx);
395 }
396 for (size_t i = 0; i < lin; ++i) {
397 dx += s;
398 points.push_back(x + dx);
399 }
400 for (size_t i = 1; i < n; ++i) {
401 s /= factor[dir];
402 dx += s;
403 points.push_back(x + dx);
404 }
405 x = *i;
406 }
407 axis->addOrderedPoints(points.begin(), points.end());
408
409 return axis;
410}
411
412template <int dim> void RectangularMeshRefinedGenerator<dim>::fromXML(XMLReader& reader, Manager& manager) {
413 if (reader.getNodeName() == "refinements") {
414 while (reader.requireTagOrEnd()) {
415 if (reader.getNodeName() != "axis0" &&
416 (dim == 1 || (reader.getNodeName() != "axis1" && (dim == 2 || reader.getNodeName() != "axis2")))) {
417 if (dim == 1) throw XMLUnexpectedElementException(reader, "<axis0>");
418 if (dim == 2) throw XMLUnexpectedElementException(reader, "<axis0> or <axis1>");
419 if (dim == 3) throw XMLUnexpectedElementException(reader, "<axis0>, <axis1>, or <axis2>");
420 }
421 PLASK_NO_CONVERSION_WARNING_BEGIN
422 auto direction =
423 (reader.getNodeName() == "axis0") ? typename Primitive<RectangularMeshRefinedGenerator<dim>::DIM>::Direction(0)
424 : (reader.getNodeName() == "axis1") ? typename Primitive<RectangularMeshRefinedGenerator<dim>::DIM>::Direction(1)
425 : typename Primitive<RectangularMeshRefinedGenerator<dim>::DIM>::Direction(2);
426 PLASK_NO_WARNING_END
427 weak_ptr<GeometryObjectD<RectangularMeshRefinedGenerator<dim>::DIM>> object =
428 manager.requireGeometryObject<GeometryObjectD<RectangularMeshRefinedGenerator<dim>::DIM>>(
429 reader.requireAttribute("object"));
430 if (!object.lock()) {
431 if (manager.draft) {
432 writelog(LOG_WARNING, "XML line {:d}: geometry object '{}' does not exist", reader.getLineNr(),
433 reader.requireAttribute("object"));
434 reader.ignoreAllAttributes();
435 reader.requireTagEnd();
436 continue;
437 } else
438 throw XMLException(reader, format("geometry object '{}' does not exist", reader.requireAttribute("object")));
439 }
440 PathHints path;
441 if (auto pathattr = reader.getAttribute("path")) path = manager.requirePathHints(*pathattr);
442 if (auto by = reader.getAttribute<unsigned>("by")) {
443 double objsize = object.lock()->getBoundingBox().size()[unsigned(direction)];
444 for (unsigned i = 1; i < *by; ++i) {
445 double pos = objsize * i / *by;
446 addRefinement(direction, object, path, pos);
447 }
448 } else if (auto every = reader.getAttribute<double>("every")) {
449 double objsize = object.lock()->getBoundingBox().size()[unsigned(direction)];
450 size_t n = int(max(round(objsize / *every), 1.));
451 double step = objsize / n;
452 for (
453 struct {
454 size_t i;
455 double pos;
456 } loop = {1, step};
457 loop.i < n; ++loop.i, loop.pos += step)
458 addRefinement(direction, object, path, loop.pos);
459 } else if (auto pos = reader.getAttribute<double>("at")) {
460 addRefinement(direction, object, path, *pos);
461 } else
462 throw XMLNoAttrException(reader, "at', 'every', or 'by");
463 reader.requireTagEnd();
464 }
465 } else if (reader.getNodeName() == "warnings") {
466 writelog(LOG_WARNING, "XML {:d}: <warnings> tag is deprecated", reader.getLineNr());
467 reader.ignoreAllAttributes();
468 reader.requireTagEnd();
469 } else
470 throw XMLUnexpectedElementException(reader, "proper generator configuration tag");
471}
472
473template struct PLASK_API RectangularMeshRefinedGenerator<1>;
474template struct PLASK_API RectangularMeshRefinedGenerator<2>;
475template struct PLASK_API RectangularMeshRefinedGenerator<3>;
476
477template <typename GeneratorT> static shared_ptr<MeshGenerator> readTrivialGenerator(XMLReader& reader, Manager&) {
478 bool split = false;
479 while (reader.requireTagOrEnd()) {
480 if (reader.getNodeName() == "boundaries") {
481 split = reader.getAttribute<bool>("split", split);
482 reader.requireTagEnd();
483 } else
484 throw XMLUnexpectedElementException(reader, "<boundaries>");
485 }
486 return plask::make_shared<GeneratorT>(split);
487}
488
489static RegisterMeshGeneratorReader ordered_simplegenerator_reader("ordered.simple",
490 readTrivialGenerator<OrderedMesh1DSimpleGenerator>);
491static RegisterMeshGeneratorReader rectangular2d_simplegenerator_reader("rectangular2d.simple",
492 readTrivialGenerator<RectangularMesh2DSimpleGenerator>);
493static RegisterMeshGeneratorReader rectangular3d_simplegenerator_reader("rectangular3d.simple",
494 readTrivialGenerator<RectangularMesh3DSimpleGenerator>);
495
496static shared_ptr<MeshGenerator> readRegularGenerator1(XMLReader& reader, Manager&) {
497 double spacing = INFINITY;
498 bool split = false;
499 while (reader.requireTagOrEnd()) {
500 if (reader.getNodeName() == "spacing") {
501 spacing = reader.getAttribute<double>("every", spacing);
502 reader.requireTagEnd();
503 } else if (reader.getNodeName() == "boundaries") {
504 split = reader.getAttribute<bool>("split", split);
505 reader.requireTagEnd();
506 } else
507 throw XMLUnexpectedElementException(reader, "<spacing>, <boundaries>");
508 }
509 return plask::make_shared<OrderedMesh1DRegularGenerator>(spacing, split);
510}
511
512static shared_ptr<MeshGenerator> readRegularGenerator2(XMLReader& reader, Manager&) {
513 double spacing0 = INFINITY, spacing1 = INFINITY;
514 bool split = false;
515 while (reader.requireTagOrEnd()) {
516 if (reader.getNodeName() == "spacing") {
517 if (reader.hasAttribute("every")) {
518 if (reader.hasAttribute("every0")) throw XMLConflictingAttributesException(reader, "every", "every0");
519 if (reader.hasAttribute("every1")) throw XMLConflictingAttributesException(reader, "every", "every1");
520 spacing0 = spacing1 = reader.requireAttribute<double>("every");
521 } else {
522 spacing0 = reader.getAttribute<double>("every0", spacing0);
523 spacing1 = reader.getAttribute<double>("every1", spacing1);
524 }
525 reader.requireTagEnd();
526 } else if (reader.getNodeName() == "boundaries") {
527 split = reader.getAttribute<bool>("split", split);
528 reader.requireTagEnd();
529 } else
530 throw XMLUnexpectedElementException(reader, "<spacing>, <boundaries>");
531 }
532 return plask::make_shared<RectangularMesh2DRegularGenerator>(spacing0, spacing1, split);
533}
534
535static shared_ptr<MeshGenerator> readRegularGenerator3(XMLReader& reader, Manager&) {
536 double spacing0 = INFINITY, spacing1 = INFINITY, spacing2 = INFINITY;
537 bool split = false;
538 while (reader.requireTagOrEnd()) {
539 if (reader.getNodeName() == "spacing") {
540 if (reader.hasAttribute("every")) {
541 if (reader.hasAttribute("every0")) throw XMLConflictingAttributesException(reader, "every", "every0");
542 if (reader.hasAttribute("every1")) throw XMLConflictingAttributesException(reader, "every", "every1");
543 if (reader.hasAttribute("every2")) throw XMLConflictingAttributesException(reader, "every", "every2");
544 spacing0 = spacing1 = reader.requireAttribute<double>("every");
545 } else {
546 spacing0 = reader.getAttribute<double>("every0", spacing0);
547 spacing1 = reader.getAttribute<double>("every1", spacing1);
548 spacing2 = reader.getAttribute<double>("every2", spacing2);
549 }
550 reader.requireTagEnd();
551 } else if (reader.getNodeName() == "boundaries") {
552 split = reader.getAttribute<bool>("split", split);
553 reader.requireTagEnd();
554 } else
555 throw XMLUnexpectedElementException(reader, "<spacing>, <boundaries>");
556 }
557 return plask::make_shared<RectangularMesh3DRegularGenerator>(spacing0, spacing1, spacing2, split);
558}
559
560static RegisterMeshGeneratorReader ordered_regulargenerator_reader("ordered.regular", readRegularGenerator1);
561static RegisterMeshGeneratorReader rectangular2d_regulargenerator_reader("rectangular2d.regular", readRegularGenerator2);
562static RegisterMeshGeneratorReader rectangular3d_regulargenerator_reader("rectangular3d.regular", readRegularGenerator3);
563
564template <int dim> shared_ptr<MeshGenerator> readRectangularDivideGenerator(XMLReader& reader, Manager& manager) {
565 auto result = plask::make_shared<RectangularMeshDivideGenerator<dim>>();
566
567 std::set<std::string> read;
568 while (reader.requireTagOrEnd()) {
569 if (read.find(reader.getNodeName()) != read.end())
570 throw XMLDuplicatedElementException(std::string("<generator>"), reader.getNodeName());
571 read.insert(reader.getNodeName());
572 if (reader.getNodeName() == "prediv") {
573 plask::optional<size_t> into = reader.getAttribute<size_t>("by");
574 if (into) {
575 if (reader.hasAttribute("by0")) throw XMLConflictingAttributesException(reader, "by", "by0");
576 if (reader.hasAttribute("by1")) throw XMLConflictingAttributesException(reader, "by", "by1");
577 if (reader.hasAttribute("by2")) throw XMLConflictingAttributesException(reader, "by", "by2");
578 for (int i = 0; i < dim; ++i) result->pre_divisions[i] = *into;
579 } else {
580 for (int i = 0; i < dim; ++i) result->pre_divisions[i] = reader.getAttribute<size_t>(format("by{0}", i), 1);
581 }
582 reader.requireTagEnd();
583 } else if (reader.getNodeName() == "postdiv") {
584 plask::optional<size_t> into = reader.getAttribute<size_t>("by");
585 if (into) {
586 if (reader.hasAttribute("by0")) throw XMLConflictingAttributesException(reader, "by", "by0");
587 if (reader.hasAttribute("by1")) throw XMLConflictingAttributesException(reader, "by", "by1");
588 if (reader.hasAttribute("by2")) throw XMLConflictingAttributesException(reader, "by", "by2");
589 for (int i = 0; i < dim; ++i) result->post_divisions[i] = *into;
590 } else {
591 for (int i = 0; i < dim; ++i) result->post_divisions[i] = reader.getAttribute<size_t>(format("by{0}", i), 1);
592 }
593 reader.requireTagEnd();
594 } else if (reader.getNodeName() == "options") {
595 plask::optional<bool> gradual = reader.getAttribute<bool>("gradual");
596 if (gradual) {
597 if (reader.hasAttribute("gradual0")) throw XMLConflictingAttributesException(reader, "gradual", "gradual0");
598 if (reader.hasAttribute("gradual1")) throw XMLConflictingAttributesException(reader, "gradual", "gradual1");
599 if (reader.hasAttribute("gradual2")) throw XMLConflictingAttributesException(reader, "gradual", "gradual2");
600 result->gradual = (*gradual) ? 7 : 0;
601 } else {
602 for (int i = 0; i < dim; ++i) {
603 result->setGradual(i, reader.getAttribute<bool>(format("gradual{0}", i), true));
604 }
605 }
606 result->setAspect(reader.getAttribute<double>("aspect", result->getAspect()));
607 reader.requireTagEnd();
608 } else
609 result->fromXML(reader, manager);
610 }
611 return result;
612}
613
614template struct PLASK_API RectangularMeshDivideGenerator<1>;
615template struct PLASK_API RectangularMeshDivideGenerator<2>;
616template struct PLASK_API RectangularMeshDivideGenerator<3>;
617
618static RegisterMeshGeneratorReader ordered_dividinggenerator_reader("ordered.divide", readRectangularDivideGenerator<1>);
619static RegisterMeshGeneratorReader rectangular2d_dividinggenerator_reader("rectangular2d.divide",
620 readRectangularDivideGenerator<2>);
621static RegisterMeshGeneratorReader rectangular3d_dividinggenerator_reader("rectangular3d.divide",
622 readRectangularDivideGenerator<3>);
623
624template <int dim> shared_ptr<MeshGenerator> readRectangularSmoothGenerator(XMLReader& reader, Manager& manager) {
625 auto result = plask::make_shared<RectangularMeshSmoothGenerator<dim>>();
626
627 std::set<std::string> read;
628 while (reader.requireTagOrEnd()) {
629 if (read.find(reader.getNodeName()) != read.end())
630 throw XMLDuplicatedElementException(std::string("<generator>"), reader.getNodeName());
631 read.insert(reader.getNodeName());
632 if (reader.getNodeName() == "steps") {
633 plask::optional<double> small_op =
634 reader.getAttribute<double>("small"); // don't use small since some windows headers: #define small char
635 if (small_op) {
636 if (reader.hasAttribute("small0")) throw XMLConflictingAttributesException(reader, "small", "small0");
637 if (reader.hasAttribute("small1")) throw XMLConflictingAttributesException(reader, "small", "small1");
638 if (reader.hasAttribute("small2")) throw XMLConflictingAttributesException(reader, "small", "small2");
639 for (int i = 0; i < dim; ++i) result->finestep[i] = *small_op;
640 } else
641 for (int i = 0; i < dim; ++i)
642 result->finestep[i] = reader.getAttribute<double>(format("small{:d}", i), result->finestep[i]);
643 plask::optional<double> large = reader.getAttribute<double>("large");
644 if (large) {
645 if (reader.hasAttribute("large0")) throw XMLConflictingAttributesException(reader, "large", "large0");
646 if (reader.hasAttribute("large1")) throw XMLConflictingAttributesException(reader, "large", "large1");
647 if (reader.hasAttribute("large2")) throw XMLConflictingAttributesException(reader, "large", "large2");
648 for (int i = 0; i < dim; ++i) result->maxstep[i] = *large;
649 } else
650 for (int i = 0; i < dim; ++i)
651 result->maxstep[i] = reader.getAttribute<double>(format("large{:d}", i), result->maxstep[i]);
652 plask::optional<double> factor = reader.getAttribute<double>("factor");
653 if (factor) {
654 if (reader.hasAttribute("factor0")) throw XMLConflictingAttributesException(reader, "factor", "factor0");
655 if (reader.hasAttribute("factor1")) throw XMLConflictingAttributesException(reader, "factor", "factor1");
656 if (reader.hasAttribute("factor2")) throw XMLConflictingAttributesException(reader, "factor", "factor2");
657 for (int i = 0; i < dim; ++i) result->factor[i] = *factor;
658 } else
659 for (int i = 0; i < dim; ++i)
660 result->factor[i] = reader.getAttribute<double>(format("factor{:d}", i), result->factor[i]);
661 reader.requireTagEnd();
662 } else if (reader.getNodeName() == "options") {
663 result->setAspect(reader.getAttribute<double>("aspect", result->getAspect()));
664 reader.requireTagEnd();
665 } else
666 result->fromXML(reader, manager);
667 }
668 return result;
669}
670
671static RegisterMeshGeneratorReader ordered_smoothgenerator_reader("ordered.smooth", readRectangularSmoothGenerator<1>);
672static RegisterMeshGeneratorReader rectangular2d_smoothgenerator_reader("rectangular2d.smooth", readRectangularSmoothGenerator<2>);
673static RegisterMeshGeneratorReader rectangular3d_smoothgenerator_reader("rectangular3d.smooth", readRectangularSmoothGenerator<3>);
674
675template struct PLASK_API RectangularMeshSmoothGenerator<1>;
676template struct PLASK_API RectangularMeshSmoothGenerator<2>;
677template struct PLASK_API RectangularMeshSmoothGenerator<3>;
678
679} // namespace plask