PLaSK library

Loading...
Searching...
No Matches
triangular2d.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 "triangular2d.hpp"
15
16#include <boost/range/irange.hpp>
17#include <boost/icl/interval_map.hpp>
18#include <unordered_map>
19
20#include "../utils/interpolation.hpp"
21#include "rectangular_common.hpp"
22
23 namespace plask {
24
25 Vec<3, double> TriangularMesh2D::Element::barycentric(Vec<2, double> p) const {
26 // formula comes from https://codeplea.com/triangular-interpolation
27 // but it is modified a bit, to reuse more diffs and call cross (which can be optimized)
28 Vec<2, double>
29 diff_2_3 = getNode(1) - getNode(2),
30 diff_p_3 = p - getNode(2),
31 diff_1_3 = getNode(0) - getNode(2);
32 const double den = cross(diff_1_3, diff_2_3); // diff_2_3.c1 * diff_1_3.c0 - diff_2_3.c0 * diff_1_3.c1
33 Vec<3, double> res;
34 res.c0 = cross(diff_p_3, diff_2_3) / den; // (diff_2_3.c1 * diff_p_3.c0 - diff_2_3.c0 * diff_p_3.c1) / den
35 res.c1 = cross(diff_1_3, diff_p_3) / den; // (- diff_1_3.c1 * diff_p_3.c0 + diff_1_3.c0 * diff_p_3.c1) / den
36 res.c2 = 1.0 - res.c0 - res.c1;
37 return res;
38}
39
40 Box2D TriangularMesh2D::Element::getBoundingBox() const {
41 LocalCoords lo = getNode(0);
42 LocalCoords up = lo;
43 const LocalCoords B = getNode(1);
44 if (B.c0 < lo.c0) lo.c0 = B.c0; else up.c0 = B.c0;
45 if (B.c1 < lo.c1) lo.c1 = B.c1; else up.c1 = B.c1;
46 const LocalCoords C = getNode(2);
47 if (C.c0 < lo.c0) lo.c0 = C.c0; else if (C.c0 > up.c0) up.c0 = C.c0;
48 if (C.c1 < lo.c1) lo.c1 = C.c1; else if (C.c1 > up.c1) up.c1 = C.c1;
49 return Box2D(lo, up);
50}
51
52 TriangularMesh2D::Builder::Builder(TriangularMesh2D &mesh): mesh(mesh) {
53 for (std::size_t i = 0; i < mesh.nodes.size(); ++i)
54 this->indexOfNode[mesh.nodes[i]] = i;
55}
56
57 TriangularMesh2D::Builder::Builder(TriangularMesh2D &mesh, std::size_t predicted_number_of_elements, std::size_t predicted_number_of_nodes)
58 : Builder(mesh)
59{
60 mesh.elementNodes.reserve(mesh.elementNodes.size() + predicted_number_of_elements);
61 mesh.nodes.reserve(mesh.nodes.size() + predicted_number_of_nodes);
62}
63
64 TriangularMesh2D::Builder::~Builder() {
65 mesh.elementNodes.shrink_to_fit();
66 mesh.nodes.shrink_to_fit();
67}
68
69 TriangularMesh2D::Builder &TriangularMesh2D::Builder::add(TriangularMesh2D::LocalCoords p1, TriangularMesh2D::LocalCoords p2, TriangularMesh2D::LocalCoords p3) {
70 mesh.elementNodes.push_back({addNode(p1), addNode(p2), addNode(p3)});
71 return *this;
72}
73
74std::size_t TriangularMesh2D::Builder::addNode(TriangularMesh2D::LocalCoords node) {
75 auto it = this->indexOfNode.emplace(node, mesh.nodes.size());
76 if (it.second) // new element has been appended to the map
77 this->mesh.nodes.push_back(node);
78 return it.first->second; // an index of the node (inserted or found)
79}
80
81 struct ElementIndexValueGetter {
82 TriangularMesh2D::ElementIndex::Rtree::value_type operator()(std::size_t index) const {
83 TriangularMesh2D::Element el = this->src_mesh->getElement(index);
84 const auto n0 = el.getNode(0);
85 const auto n1 = el.getNode(1);
86 const auto n2 = el.getNode(2);
87 return std::make_pair(
88 TriangularMesh2D::ElementIndex::Box(
89 vec(std::min(std::min(n0.c0, n1.c0), n2.c0), std::min(std::min(n0.c1, n1.c1), n2.c1)),
90 vec(std::max(std::max(n0.c0, n1.c0), n2.c0), std::max(std::max(n0.c1, n1.c1), n2.c1))
91 ),
92 index);
93 }
94
95 const TriangularMesh2D* src_mesh;
96
97 ElementIndexValueGetter(const TriangularMesh2D& src_mesh): src_mesh(&src_mesh) {}
98};
99
100 TriangularMesh2D::ElementIndex::ElementIndex(const TriangularMesh2D &mesh)
101 : mesh(mesh),
102 rtree(makeFunctorIndexedIterator(ElementIndexValueGetter(mesh), 0),
103 makeFunctorIndexedIterator(ElementIndexValueGetter(mesh), mesh.getElementsCount()))
104{}
105
106std::size_t TriangularMesh2D::ElementIndex::getIndex(Vec<2, double> p) const
107{
108 for (auto v: rtree | boost::geometry::index::adaptors::queried(boost::geometry::index::intersects(p))) {
109 const Element el = mesh.getElement(v.second);
110 if (el.contains(p)) return v.second;
111 }
112 return INDEX_NOT_FOUND;
113}
114
115 optional<TriangularMesh2D::Element> TriangularMesh2D::ElementIndex::getElement(Vec<2, double> p) const {
116 for (auto v: rtree | boost::geometry::index::adaptors::queried(boost::geometry::index::intersects(p))) {
117 const Element el = mesh.getElement(v.second);
118 if (el.contains(p)) return el;
119 }
120 return optional<TriangularMesh2D::Element>();
121}
122
123 TriangularMesh2D TriangularMesh2D::masked(const TriangularMesh2D::Predicate &predicate) const {
124 TriangularMesh2D result;
125 Builder builder(result, elementNodes.size(), nodes.size());
126 for (auto el: elements())
127 if (predicate(el)) builder.add(el);
128 return result;
129}
130
131 void TriangularMesh2D::writeXML(XMLElement &object) const {
132 object.attr("type", "triangular2d");
133 for (auto & node: nodes)
134 object.addTag("node").attr("tran", node.tran()).attr("vert", node.vert());
135 for (auto & el: elementNodes)
136 object.addTag("element").attr("a", el[0]).attr("b", el[1]).attr("c", el[2]);
137}
138
139 inline TriangularMesh2D::Segment makeSegment(std::size_t a, std::size_t b) {
140 return a < b ? std::make_pair(a, b) : std::make_pair(b, a);
141}
142
143 void TriangularMesh2D::countSegmentsOf(TriangularMesh2D::SegmentsCounts& counter, const TriangularMesh2D::Element& el) {
144 ++counter[makeSegment(el.getNodeIndex(0), el.getNodeIndex(1))];
145 ++counter[makeSegment(el.getNodeIndex(1), el.getNodeIndex(2))];
146 ++counter[makeSegment(el.getNodeIndex(2), el.getNodeIndex(0))];
147}
148
149 TriangularMesh2D::SegmentsCounts TriangularMesh2D::countSegments() const {
150 TriangularMesh2D::SegmentsCounts result;
151 for (const auto el: elements()) countSegmentsOf(result, el);
152 return result;
153}
154
155 TriangularMesh2D::SegmentsCounts TriangularMesh2D::countSegmentsIn(const Box2D &box) const {
156 TriangularMesh2D::SegmentsCounts result;
157 for (const auto el: elements())
158 if (box.contains(el.getNode(0)) && box.contains(el.getNode(1)) && box.contains(el.getNode(2)))
159 countSegmentsOf(result, el);
160 return result;
161}
162
163 TriangularMesh2D::SegmentsCounts TriangularMesh2D::countSegmentsIn(const std::vector<Box2D>& boxes) const {
164 TriangularMesh2D::SegmentsCounts result;
165 for (const auto el: elements()) {
166 bool vertex0_included = false, vertex1_included = false, vertex2_included = false;
167 for (auto& box: boxes) {
168 if (!vertex0_included) vertex0_included = box.contains(el.getNode(0));
169 if (!vertex1_included) vertex1_included = box.contains(el.getNode(1));
170 if (!vertex2_included) vertex2_included = box.contains(el.getNode(2));
171 if (vertex0_included && vertex1_included && vertex2_included) {
172 countSegmentsOf(result, el);
173 break;
174 }
175 }
176 }
177 return result;
178}
179
180 TriangularMesh2D::SegmentsCounts TriangularMesh2D::countSegmentsIn(const GeometryD<2> &geometry, const GeometryObject &object, const PathHints *path) const {
181 TriangularMesh2D::SegmentsCounts result;
182 for (const auto el: elements())
183 if (geometry.objectIncludes(object, path, el.getNode(0)) &&
184 geometry.objectIncludes(object, path, el.getNode(1)) &&
185 geometry.objectIncludes(object, path, el.getNode(2)))
186 countSegmentsOf(result, el);
187 return result;
188}
189
190std::set<std::size_t> TriangularMesh2D::allBoundaryNodes(const TriangularMesh2D::SegmentsCounts& segmentsCount) {
191 std::set<std::size_t> result;
192 for (const std::pair<TriangularMesh2D::Segment, std::size_t>& s: segmentsCount)
193 if (s.second == 1) {
194 result.insert(s.first.first);
195 result.insert(s.first.second);
196 }
197 return result;
198}
199
200 template <int DIR, template<class> class Compare = std::less> // DIR - oś prostopadła do przedziałów trzymanych w boost::icl::interval_map
201 struct SegmentSetMember: private Compare<TriangularMesh2D::Segment> {
202
203 constexpr static int SEG_DIR = 1 - DIR; // kierunek rozciągania się przedziałów w interval_map
204
205 TriangularMesh2D::Segment segment;
206 double min, max; // mniejsza i większa współrzędna (w osi DIR) końców segmentu
207
208 SegmentSetMember() = default;
209
210 SegmentSetMember(const TriangularMesh2D& mesh, const TriangularMesh2D::Segment& segment)
211 : segment(segment),
212 min(mesh.nodes[segment.first][DIR]),
213 max(mesh.nodes[segment.second][DIR])
214 {
215 //if (max < min) std::swap(min, max);
216 if (Compare<double>()(max, min)) // max < min ?
217 std::swap(min, max);
218 }
219
220 bool operator<(const SegmentSetMember& other) const {
221 return Compare<TriangularMesh2D::Segment>::operator()(this->segment, other.segment);
222 //this->segment < other.segment;
223 }
224
225 bool operator==(const SegmentSetMember& other) const {
226 return this->segment == other.segment;
227 }
228
235 bool dominates(const TriangularMesh2D& mesh, const Vec<2, double>& point) const {
236 // this two cases are supported also by interpolation, but they are common and we want to ensure to check them exact:
237 const Vec<2, double>& f = mesh.nodes[segment.first];
238 if (point[SEG_DIR] == f[SEG_DIR])
239 return Compare<double>()(point[DIR], f[DIR]); //point[DIR] < f[DIR];
240 const Vec<2, double>& s = mesh.nodes[segment.second];
241 if (point[SEG_DIR] == s[SEG_DIR])
242 return Compare<double>()(point[DIR], s[DIR]); //point[DIR] < f[DIR];
243
244 if (f[SEG_DIR] < s[SEG_DIR])
245 return Compare<double>()(point[DIR], interpolation::linear(f[SEG_DIR], f[DIR], s[SEG_DIR], s[DIR], point[SEG_DIR]));
246 //point[DIR] < interpolation::linear(f[SEG_DIR], f[DIR], s[SEG_DIR], s[DIR], point[SEG_DIR]);
247 else // same as above, but f<->s are swapped
248 return Compare<double>()(point[DIR], interpolation::linear(s[SEG_DIR], s[DIR], f[SEG_DIR], f[DIR], point[SEG_DIR]));
249 }
250};
251
252// DIR - oś prostopadła do przedziałów trzymanych w boost::icl::interval_map
253 template<int DIR, template<class> class Compare = std::less>
254 struct SegmentSet: Compare<double> {
255 SegmentSet() = default;
256
257 SegmentSet(const SegmentSetMember<DIR, Compare>& to_append)
258 : maxOfMins(to_append.min) { set.insert(to_append); }
259
260 SegmentSet(const TriangularMesh2D& mesh, const TriangularMesh2D::Segment& segment)
261 : SegmentSet(SegmentSetMember<DIR, Compare>(mesh, segment)) {}
262
263 SegmentSet& operator+=(const SegmentSet& right) {
264 //if (maxOfMins < right.maxOfMins) {
265 if (Compare<double>::operator()(maxOfMins, right.maxOfMins)) {
266 SetT this_set_backup = std::move(set);
267 set = right.set;
268 maxOfMins = right.maxOfMins;
269 insert(this_set_backup);
270 } else
271 insert(right.set);
272 return *this;
273 }
274
275 bool operator==(const SegmentSet& other) const {
276 return this->set == other.set;
277 }
278
279 // check if point is dominated by any segment in set represented by this
280 bool dominates(const TriangularMesh2D& mesh, const Vec<2, double>& point) const {
281 for (const SegmentSetMember<DIR, Compare>& s: set)
282 if (s.dominates(mesh, point)) return true;
283 return false;
284 }
285
286 private:
287 typedef std::set<SegmentSetMember<DIR, Compare>, Compare<SegmentSetMember<DIR, Compare>>> SetT;
288
289 void insert(const SetT& right) {
290 for (auto& s: right)
291 if (!Compare<double>::operator()(s.max, maxOfMins))
292 set.insert(s);
293 //if (s.max >= maxOfMins) set.insert(s);
294 }
295
296 SetT set;
297
303 double maxOfMins;
304};
305
306 template<int SEG_DIR, template<class> class Compare>
307std::set<std::size_t> TriangularMesh2D::dirBoundaryNodes(const TriangularMesh2D::SegmentsCounts &segmentsCount) const {
308 typedef SegmentSet<1-SEG_DIR, Compare> SegmentSetT;
309 std::set<std::size_t> result;
310 boost::icl::interval_map<double, SegmentSetT> non_dominated;
311 for (const std::pair<TriangularMesh2D::Segment, std::size_t>& s: segmentsCount)
312 if (s.second == 1) {
313 const TriangularMesh2D::Segment& seg = s.first;
314 non_dominated += std::make_pair(
315 boost::icl::interval<double>::closed(this->nodes[seg.first][SEG_DIR], this->nodes[seg.second][SEG_DIR]),
316 SegmentSetT(*this, seg)
317 );
318 result.insert(s.first.first);
319 result.insert(s.first.second);
320 }
321 // remove nodes which are dominated:
322 for (auto it = result.begin(); it != result.end(); ) {
323 const auto& node_coords = this->nodes[*it];
324 if (non_dominated.find(node_coords[SEG_DIR])->second.dominates(*this, node_coords))
325 it = result.erase(it);
326 else
327 ++it;
328 }
329 return result;
330}
331
332 TriangularMesh2D::Boundary TriangularMesh2D::getRightBoundary() {
333 return Boundary( [](const TriangularMesh2D& mesh, const shared_ptr<const GeometryD<2>>&) {
334 return BoundaryNodeSet(new StdSetBoundaryImpl(mesh.boundaryNodes<BoundaryDir::RIGHT>(mesh.countSegments())));
335 } );
336}
337
338 TriangularMesh2D::Boundary TriangularMesh2D::getTopBoundary() {
339 return Boundary( [](const TriangularMesh2D& mesh, const shared_ptr<const GeometryD<2>>&) {
340 return BoundaryNodeSet(new StdSetBoundaryImpl(mesh.boundaryNodes<BoundaryDir::TOP>(mesh.countSegments())));
341 } );
342}
343
344 TriangularMesh2D::Boundary TriangularMesh2D::getLeftBoundary() {
345 return Boundary( [](const TriangularMesh2D& mesh, const shared_ptr<const GeometryD<2>>&) {
346 return BoundaryNodeSet(new StdSetBoundaryImpl(mesh.boundaryNodes<BoundaryDir::LEFT>(mesh.countSegments())));
347 } );
348}
349
350 TriangularMesh2D::Boundary TriangularMesh2D::getBottomBoundary() {
351 return Boundary( [](const TriangularMesh2D& mesh, const shared_ptr<const GeometryD<2>>&) {
352 return BoundaryNodeSet(new StdSetBoundaryImpl(mesh.boundaryNodes<BoundaryDir::BOTTOM>(mesh.countSegments())));
353 } );
354}
355
356 TriangularMesh2D::Boundary TriangularMesh2D::getRightOfBoundary(const Box2D &box) {
357 return Boundary( [box](const TriangularMesh2D& mesh, const shared_ptr<const GeometryD<2>>&) {
358 return BoundaryNodeSet(new StdSetBoundaryImpl(mesh.boundaryNodes<BoundaryDir::RIGHT>(mesh.countSegmentsIn(box))));
359 } );
360}
361
362 TriangularMesh2D::Boundary TriangularMesh2D::getTopOfBoundary(const Box2D &box) {
363 return Boundary( [box](const TriangularMesh2D& mesh, const shared_ptr<const GeometryD<2>>&) {
364 return BoundaryNodeSet(new StdSetBoundaryImpl(mesh.boundaryNodes<BoundaryDir::TOP>(mesh.countSegmentsIn(box))));
365 } );
366}
367
368 TriangularMesh2D::Boundary TriangularMesh2D::getLeftOfBoundary(const Box2D &box) {
369 return Boundary( [box](const TriangularMesh2D& mesh, const shared_ptr<const GeometryD<2>>&) {
370 return BoundaryNodeSet(new StdSetBoundaryImpl(mesh.boundaryNodes<BoundaryDir::LEFT>(mesh.countSegmentsIn(box))));
371 } );
372}
373
374 TriangularMesh2D::Boundary TriangularMesh2D::getBottomOfBoundary(const Box2D &box) {
375 return Boundary( [box](const TriangularMesh2D& mesh, const shared_ptr<const GeometryD<2>>&) {
376 return BoundaryNodeSet(new StdSetBoundaryImpl(mesh.boundaryNodes<BoundaryDir::BOTTOM>(mesh.countSegmentsIn(box))));
377 } );
378}
379
380 TriangularMesh2D::Boundary TriangularMesh2D::getRightOfBoundary(const std::vector<Box2D> &boxes) {
381 return Boundary( [boxes](const TriangularMesh2D& mesh, const shared_ptr<const GeometryD<2>>&) {
382 return BoundaryNodeSet(new StdSetBoundaryImpl(mesh.boundaryNodes<BoundaryDir::RIGHT>(mesh.countSegmentsIn(boxes))));
383 } );
384}
385
386 TriangularMesh2D::Boundary TriangularMesh2D::getTopOfBoundary(const std::vector<Box2D> &boxes) {
387 return Boundary( [boxes](const TriangularMesh2D& mesh, const shared_ptr<const GeometryD<2>>&) {
388 return BoundaryNodeSet(new StdSetBoundaryImpl(mesh.boundaryNodes<BoundaryDir::TOP>(mesh.countSegmentsIn(boxes))));
389 } );
390}
391
392 TriangularMesh2D::Boundary TriangularMesh2D::getLeftOfBoundary(const std::vector<Box2D> &boxes) {
393 return Boundary( [boxes](const TriangularMesh2D& mesh, const shared_ptr<const GeometryD<2>>&) {
394 return BoundaryNodeSet(new StdSetBoundaryImpl(mesh.boundaryNodes<BoundaryDir::LEFT>(mesh.countSegmentsIn(boxes))));
395 } );
396}
397
398 TriangularMesh2D::Boundary TriangularMesh2D::getBottomOfBoundary(const std::vector<Box2D> &boxes) {
399 return Boundary( [boxes](const TriangularMesh2D& mesh, const shared_ptr<const GeometryD<2>>&) {
400 return BoundaryNodeSet(new StdSetBoundaryImpl(mesh.boundaryNodes<BoundaryDir::BOTTOM>(mesh.countSegmentsIn(boxes))));
401 } );
402}
403
404 TriangularMesh2D::Boundary TriangularMesh2D::getRightOfBoundary(shared_ptr<const GeometryObject> object) {
405 return Boundary( [=](const TriangularMesh2D& mesh, const shared_ptr<const GeometryD<2>>& geom) {
406 return BoundaryNodeSet(new StdSetBoundaryImpl(mesh.boundaryNodes<BoundaryDir::RIGHT>(mesh.countSegmentsIn(*geom, *object))));
407 } );
408}
409
410 TriangularMesh2D::Boundary TriangularMesh2D::getRightOfBoundary(shared_ptr<const GeometryObject> object, const PathHints &path) {
411 return Boundary( [=](const TriangularMesh2D& mesh, const shared_ptr<const GeometryD<2>>& geom) {
412 return BoundaryNodeSet(new StdSetBoundaryImpl(mesh.boundaryNodes<BoundaryDir::RIGHT>(mesh.countSegmentsIn(*geom, *object, &path))));
413 } );
414}
415
416 TriangularMesh2D::Boundary TriangularMesh2D::getTopOfBoundary(shared_ptr<const GeometryObject> object) {
417 return Boundary( [=](const TriangularMesh2D& mesh, const shared_ptr<const GeometryD<2>>& geom) {
418 return BoundaryNodeSet(new StdSetBoundaryImpl(mesh.boundaryNodes<BoundaryDir::TOP>(mesh.countSegmentsIn(*geom, *object))));
419 } );
420}
421
422 TriangularMesh2D::Boundary TriangularMesh2D::getTopOfBoundary(shared_ptr<const GeometryObject> object, const PathHints &path) {
423 return Boundary( [=](const TriangularMesh2D& mesh, const shared_ptr<const GeometryD<2>>& geom) {
424 return BoundaryNodeSet(new StdSetBoundaryImpl(mesh.boundaryNodes<BoundaryDir::TOP>(mesh.countSegmentsIn(*geom, *object, &path))));
425 } );
426}
427
428 TriangularMesh2D::Boundary TriangularMesh2D::getLeftOfBoundary(shared_ptr<const GeometryObject> object) {
429 return Boundary( [=](const TriangularMesh2D& mesh, const shared_ptr<const GeometryD<2>>& geom) {
430 return BoundaryNodeSet(new StdSetBoundaryImpl(mesh.boundaryNodes<BoundaryDir::LEFT>(mesh.countSegmentsIn(*geom, *object))));
431 } );
432}
433
434 TriangularMesh2D::Boundary TriangularMesh2D::getLeftOfBoundary(shared_ptr<const GeometryObject> object, const PathHints &path) {
435 return Boundary( [=](const TriangularMesh2D& mesh, const shared_ptr<const GeometryD<2>>& geom) {
436 return BoundaryNodeSet(new StdSetBoundaryImpl(mesh.boundaryNodes<BoundaryDir::LEFT>(mesh.countSegmentsIn(*geom, *object, &path))));
437 } );
438}
439
440 TriangularMesh2D::Boundary TriangularMesh2D::getBottomOfBoundary(shared_ptr<const GeometryObject> object) {
441 return Boundary( [=](const TriangularMesh2D& mesh, const shared_ptr<const GeometryD<2>>& geom) {
442 return BoundaryNodeSet(new StdSetBoundaryImpl(mesh.boundaryNodes<BoundaryDir::BOTTOM>(mesh.countSegmentsIn(*geom, *object))));
443 } );
444}
445
446 TriangularMesh2D::Boundary TriangularMesh2D::getBottomOfBoundary(shared_ptr<const GeometryObject> object, const PathHints &path) {
447 return Boundary( [=](const TriangularMesh2D& mesh, const shared_ptr<const GeometryD<2>>& geom) {
448 return BoundaryNodeSet(new StdSetBoundaryImpl(mesh.boundaryNodes<BoundaryDir::BOTTOM>(mesh.countSegmentsIn(*geom, *object, &path))));
449 } );
450}
451
452 TriangularMesh2D::Boundary TriangularMesh2D::getBoundary(const std::string &boundary_desc) {
453 if (boundary_desc == "bottom") return getBottomBoundary();
454 if (boundary_desc == "left") return getLeftBoundary();
455 if (boundary_desc == "right") return getRightBoundary();
456 if (boundary_desc == "top") return getTopBoundary();
457 if (boundary_desc == "all") return getAllBoundary();
458 return Boundary();
459}
460
461 TriangularMesh2D::Boundary TriangularMesh2D::getBoundary(XMLReader &boundary_desc, Manager &manager) {
462 auto side = boundary_desc.getAttribute("side");
463 if (side) {
464 if (*side == "bottom")
465 return details::parseBoundaryFromXML<Boundary, 2>(boundary_desc, manager, &getBottomBoundary, &getBottomOfBoundary);
466 if (*side == "left")
467 return details::parseBoundaryFromXML<Boundary, 2>(boundary_desc, manager, &getLeftBoundary, &getLeftOfBoundary);
468 if (*side == "right")
469 return details::parseBoundaryFromXML<Boundary, 2>(boundary_desc, manager, &getRightBoundary, &getRightOfBoundary);
470 if (*side == "top")
471 return details::parseBoundaryFromXML<Boundary, 2>(boundary_desc, manager, &getTopBoundary, &getTopOfBoundary);
472 if (*side == "all")
473 return details::parseBoundaryFromXML<Boundary, 2>(boundary_desc, manager, &getAllBoundary, &getAllBoundaryIn);
474 throw XMLBadAttrException(boundary_desc, "side", *side);
475 }
476 return Boundary();
477}
478
479std::size_t readTriangularMesh2D_readNodeIndex(XMLReader& reader, const char* attrName, std::size_t nodes_size) {
480 std::size_t result = reader.requireAttribute<std::size_t>(attrName);
481 if (result >= nodes_size)
482 reader.throwException(format("{} in <element> equals {} and is out of range [0, {})", attrName, result, nodes_size));
483 return result;
484}
485
486 TriangularMesh2D TriangularMesh2D::read(XMLReader& reader) {
487 TriangularMesh2D result;
488
489 if (reader.requireTagOrEnd()) { // has tag?
490 std::string tag_name = reader.getNodeName();
491 if (tag_name == "triangle") { // sequence of triangles
492 TriangularMesh2D::Builder builder(result);
493 do {
494 builder.add(
495 vec(reader.requireAttribute<double>("a0"), reader.requireAttribute<double>("a1")),
496 vec(reader.requireAttribute<double>("b0"), reader.requireAttribute<double>("b1")),
497 vec(reader.requireAttribute<double>("c0"), reader.requireAttribute<double>("c1"))
498 );
499 reader.requireTagEnd();
500 } while (reader.requireTagOrEnd("triangle"));
501 } else if (tag_name == "node") {
502 result.nodes.emplace_back(reader.requireAttribute<double>("tran"), reader.requireAttribute<double>("vert"));
503 reader.requireTagEnd();
504 bool accept_nodes = true; // accept <node> and <element> tags if true, else accept only <element>
505 while (reader.requireTagOrEnd()) {
506 std::string tag_name = reader.getNodeName();
507 if (accept_nodes && tag_name == "node") {
508 result.nodes.emplace_back(reader.requireAttribute<double>("tran"), reader.requireAttribute<double>("vert"));
509 reader.requireTagEnd();
510 } else if (tag_name == "element") {
511 result.elementNodes.push_back({
512 readTriangularMesh2D_readNodeIndex(reader, "a", result.nodes.size()),
513 readTriangularMesh2D_readNodeIndex(reader, "b", result.nodes.size()),
514 readTriangularMesh2D_readNodeIndex(reader, "c", result.nodes.size())
515 });
516 reader.requireTagEnd();
517 accept_nodes = false;
518 }
519 }
520 } else
521 reader.throwUnexpectedElementException("expected <triangle> or <node> tag, got <" + tag_name + ">");
522 }
523
524 return result;
525}
526
527 static shared_ptr<Mesh> readTriangularMesh2D(XMLReader& reader) {
528 return plask::make_shared<TriangularMesh2D>(TriangularMesh2D::read(reader));
529}
530
531 static RegisterMeshReader rectangular2d_reader("triangular2d", readTriangularMesh2D);
532
533 template PLASK_API std::set<std::size_t> TriangularMesh2D::boundaryNodes<TriangularMesh2D::BoundaryDir::TOP>(const TriangularMesh2D::SegmentsCounts& segmentsCount) const;
534 template PLASK_API std::set<std::size_t> TriangularMesh2D::boundaryNodes<TriangularMesh2D::BoundaryDir::LEFT>(const TriangularMesh2D::SegmentsCounts& segmentsCount) const;
535 template PLASK_API std::set<std::size_t> TriangularMesh2D::boundaryNodes<TriangularMesh2D::BoundaryDir::RIGHT>(const TriangularMesh2D::SegmentsCounts& segmentsCount) const;
536 template PLASK_API std::set<std::size_t> TriangularMesh2D::boundaryNodes<TriangularMesh2D::BoundaryDir::BOTTOM>(const TriangularMesh2D::SegmentsCounts& segmentsCount) const;
537 template PLASK_API std::set<std::size_t> TriangularMesh2D::boundaryNodes<TriangularMesh2D::BoundaryDir::ALL>(const TriangularMesh2D::SegmentsCounts& segmentsCount) const;
538
539// ------------------ Nearest Neighbor interpolation ---------------------
540
541 template<typename DstT, typename SrcT>
542 NearestNeighborTriangularMesh2DLazyDataImpl<DstT, SrcT>::NearestNeighborTriangularMesh2DLazyDataImpl(
543 const shared_ptr<const TriangularMesh2D>& src_mesh,
544 const DataVector<const SrcT>& src_vec,
545 const shared_ptr<const MeshD<2>>& dst_mesh,
546 const InterpolationFlags& flags)
547 : InterpolatedLazyDataImpl<DstT, TriangularMesh2D, const SrcT>(src_mesh, src_vec, dst_mesh, flags),
548 nodesIndex(boost::irange(std::size_t(0), src_mesh->size()),
549 typename RtreeOfTriangularMesh2DNodes::parameters_type(),
550 TriangularMesh2DGetterForRtree(src_mesh.get()))
551{
552}
553
554 template <typename DstT, typename SrcT>
555 DstT NearestNeighborTriangularMesh2DLazyDataImpl<DstT, SrcT>::at(std::size_t index) const
556{
557 auto point = this->dst_mesh->at(index);
558 auto wrapped_point = this->flags.wrap(point);
559 for (auto v: nodesIndex | boost::geometry::index::adaptors::queried(boost::geometry::index::nearest(wrapped_point, 1)))
560 return this->flags.postprocess(point, this->src_vec[v]);
561 return DstT();
562}
563
564 template struct PLASK_API NearestNeighborTriangularMesh2DLazyDataImpl<double, double>;
565 template struct PLASK_API NearestNeighborTriangularMesh2DLazyDataImpl<dcomplex, dcomplex>;
566 template struct PLASK_API NearestNeighborTriangularMesh2DLazyDataImpl<Vec<2,double>, Vec<2,double>>;
567 template struct PLASK_API NearestNeighborTriangularMesh2DLazyDataImpl<Vec<2,dcomplex>, Vec<2,dcomplex>>;
568 template struct PLASK_API NearestNeighborTriangularMesh2DLazyDataImpl<Vec<3,double>, Vec<3,double>>;
569 template struct PLASK_API NearestNeighborTriangularMesh2DLazyDataImpl<Vec<3,dcomplex>, Vec<3,dcomplex>>;
570 template struct PLASK_API NearestNeighborTriangularMesh2DLazyDataImpl<Tensor2<double>, Tensor2<double>>;
571 template struct PLASK_API NearestNeighborTriangularMesh2DLazyDataImpl<Tensor2<dcomplex>, Tensor2<dcomplex>>;
572 template struct PLASK_API NearestNeighborTriangularMesh2DLazyDataImpl<Tensor3<double>, Tensor3<double>>;
573 template struct PLASK_API NearestNeighborTriangularMesh2DLazyDataImpl<Tensor3<dcomplex>, Tensor3<dcomplex>>;
574
575
576// ------------------ Barycentric / Linear interpolation ---------------------
577
578 template<typename DstT, typename SrcT>
579 BarycentricTriangularMesh2DLazyDataImpl<DstT, SrcT>::BarycentricTriangularMesh2DLazyDataImpl(const shared_ptr<const TriangularMesh2D> &src_mesh, const DataVector<const SrcT> &src_vec, const shared_ptr<const MeshD<2> > &dst_mesh, const InterpolationFlags &flags)
580 : InterpolatedLazyDataImpl<DstT, TriangularMesh2D, const SrcT>(src_mesh, src_vec, dst_mesh, flags),
581 elementIndex(*src_mesh)
582{
583}
584
585 template <typename DstT, typename SrcT>
586 DstT BarycentricTriangularMesh2DLazyDataImpl<DstT, SrcT>::at(std::size_t index) const
587{
588 auto point = this->dst_mesh->at(index);
589 auto wrapped_point = this->flags.wrap(point);
590 for (auto v: elementIndex.rtree | boost::geometry::index::adaptors::queried(boost::geometry::index::intersects(wrapped_point))) {
591 const auto el = this->src_mesh->getElement(v.second);
592 const auto b = el.barycentric(wrapped_point);
593 if (b.c0 < 0.0 || b.c1 < 0.0 || b.c2 < 0.0) continue; // wrapped_point is outside of the triangle
594 return this->flags.postprocess(point,
595 b.c0*this->src_vec[el.getNodeIndex(0)] +
596 b.c1*this->src_vec[el.getNodeIndex(1)] +
597 b.c2*this->src_vec[el.getNodeIndex(2)]);
598 }
599 return NaN<decltype(this->src_vec[0])>();
600}
601
602 template struct PLASK_API BarycentricTriangularMesh2DLazyDataImpl<double, double>;
603 template struct PLASK_API BarycentricTriangularMesh2DLazyDataImpl<dcomplex, dcomplex>;
604 template struct PLASK_API BarycentricTriangularMesh2DLazyDataImpl<Vec<2,double>, Vec<2,double>>;
605 template struct PLASK_API BarycentricTriangularMesh2DLazyDataImpl<Vec<2,dcomplex>, Vec<2,dcomplex>>;
606 template struct PLASK_API BarycentricTriangularMesh2DLazyDataImpl<Vec<3,double>, Vec<3,double>>;
607 template struct PLASK_API BarycentricTriangularMesh2DLazyDataImpl<Vec<3,dcomplex>, Vec<3,dcomplex>>;
608 template struct PLASK_API BarycentricTriangularMesh2DLazyDataImpl<Tensor2<double>, Tensor2<double>>;
609 template struct PLASK_API BarycentricTriangularMesh2DLazyDataImpl<Tensor2<dcomplex>, Tensor2<dcomplex>>;
610 template struct PLASK_API BarycentricTriangularMesh2DLazyDataImpl<Tensor3<double>, Tensor3<double>>;
611 template struct PLASK_API BarycentricTriangularMesh2DLazyDataImpl<Tensor3<dcomplex>, Tensor3<dcomplex>>;
612
613
614// ------------------ Element mesh Nearest Neighbor interpolation ---------------------
615
616 template<typename DstT, typename SrcT>
617 NearestNeighborElementTriangularMesh2DLazyDataImpl<DstT, SrcT>::NearestNeighborElementTriangularMesh2DLazyDataImpl(
618 const shared_ptr<const TriangularMesh2D::ElementMesh>& src_mesh,
619 const DataVector<const SrcT>& src_vec,
620 const shared_ptr<const MeshD<2>>& dst_mesh,
621 const InterpolationFlags& flags)
622 : InterpolatedLazyDataImpl<DstT, TriangularMesh2D::ElementMesh, const SrcT>(src_mesh, src_vec, dst_mesh, flags),
623 elementIndex(src_mesh->getOriginalMesh())
624{
625}
626
627 template<typename DstT, typename SrcT>
628 DstT NearestNeighborElementTriangularMesh2DLazyDataImpl<DstT, SrcT>::at(std::size_t index) const {
629 auto point = this->dst_mesh->at(index);
630 auto wrapped_point = this->flags.wrap(point);
631 std::size_t element_index = elementIndex.getIndex(wrapped_point);
632 if (element_index == TriangularMesh2D::ElementIndex::INDEX_NOT_FOUND)
633 return NaN<decltype(this->src_vec[0])>();
634 return this->flags.postprocess(point, this->src_vec[element_index]);
635}
636
637 template struct PLASK_API NearestNeighborElementTriangularMesh2DLazyDataImpl<double, double>;
638 template struct PLASK_API NearestNeighborElementTriangularMesh2DLazyDataImpl<dcomplex, dcomplex>;
639 template struct PLASK_API NearestNeighborElementTriangularMesh2DLazyDataImpl<Vec<2,double>, Vec<2,double>>;
640 template struct PLASK_API NearestNeighborElementTriangularMesh2DLazyDataImpl<Vec<2,dcomplex>, Vec<2,dcomplex>>;
641 template struct PLASK_API NearestNeighborElementTriangularMesh2DLazyDataImpl<Vec<3,double>, Vec<3,double>>;
642 template struct PLASK_API NearestNeighborElementTriangularMesh2DLazyDataImpl<Vec<3,dcomplex>, Vec<3,dcomplex>>;
643 template struct PLASK_API NearestNeighborElementTriangularMesh2DLazyDataImpl<Tensor2<double>, Tensor2<double>>;
644 template struct PLASK_API NearestNeighborElementTriangularMesh2DLazyDataImpl<Tensor2<dcomplex>, Tensor2<dcomplex>>;
645 template struct PLASK_API NearestNeighborElementTriangularMesh2DLazyDataImpl<Tensor3<double>, Tensor3<double>>;
646 template struct PLASK_API NearestNeighborElementTriangularMesh2DLazyDataImpl<Tensor3<dcomplex>, Tensor3<dcomplex>>;
647
648 bool TriangularMesh2D::ElementMesh::hasSameNodes(const MeshD<2> &to_compare) const {
649 if (const ElementMesh* c = dynamic_cast<const ElementMesh*>(&to_compare))
650 if (this->originalMesh == c->originalMesh) return true;
651 return MeshD<2>::hasSameNodes(to_compare);
652}
653
654} // namespace plask