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