PLaSK library

Loading...
Searching...
No Matches
extruded_triangular3d.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 "extruded_triangular3d.hpp"
15
16#include <boost/range/irange.hpp>
17
18#include "../utils/interpolation.hpp"
19#include "rectangular_common.hpp"
20
21 namespace plask {
22
23 inline Vec<3, double> from_longTran_vert(const Vec<2, double>& longTran, const double& vert) {
24 return vec(longTran.c0, longTran.c1, vert);
25}
26
27 inline Vec<2, double> to_longTran(const Vec<3, double>& longTranVert) {
28 return Vec<2, double>(longTranVert.c0, longTranVert.c1);
29}
30
31 inline Box2D to_longTran(const Box3D& box) {
32 return Box2D(to_longTran(box.lower), to_longTran(box.upper));
33}
34
35 ExtrudedTriangularMesh3D::Element::Element(const ExtrudedTriangularMesh3D &mesh, std::size_t elementIndex)
36 : mesh(mesh)
37{
38 if (mesh.vertFastest) {
39 const std::size_t seg_size = mesh.vertAxis->size()-1;
40 longTranIndex = elementIndex / seg_size;
41 vertIndex = elementIndex % seg_size;
42 } else {
43 const std::size_t seg_size = mesh.longTranMesh.getElementsCount();
44 longTranIndex = elementIndex % seg_size;
45 vertIndex = elementIndex / seg_size;
46 }
47}
48
49 Vec<3, double> ExtrudedTriangularMesh3D::Element::getMidpoint() const {
50 return from_longTran_vert(
51 longTranElement().getMidpoint(),
52 (mesh.vertAxis->at(vertIndex) + mesh.vertAxis->at(vertIndex+1)) * 0.5
53 );
54}
55
56 bool ExtrudedTriangularMesh3D::Element::contains(Vec<3, double> p) const {
57 return mesh.vertAxis->at(vertIndex) <= p.vert() && p.vert() <= mesh.vertAxis->at(vertIndex+1)
58 && longTranElement().contains(to_longTran(p));
59}
60
61 Box3D ExtrudedTriangularMesh3D::Element::getBoundingBox() const {
62 Box2D ltBox = longTranElement().getBoundingBox();
63 return Box3D(
64 from_longTran_vert(ltBox.lower, mesh.vertAxis->at(vertIndex)),
65 from_longTran_vert(ltBox.upper, mesh.vertAxis->at(vertIndex+1))
66 );
67}
68
69 Vec<3, double> ExtrudedTriangularMesh3D::at(std::size_t index) const {
70 if (vertFastest) {
71 const std::size_t seg_size = vertAxis->size();
72 return at(index / seg_size, index % seg_size);
73 } else {
74 const std::size_t seg_size = longTranMesh.size();
75 return at(index % seg_size, index / seg_size);
76 }
77}
78
79std::size_t ExtrudedTriangularMesh3D::size() const {
80 return longTranMesh.size() * vertAxis->size();
81}
82
83 bool ExtrudedTriangularMesh3D::empty() const {
84 return longTranMesh.empty() || vertAxis->empty();
85}
86
87 void ExtrudedTriangularMesh3D::writeXML(XMLElement &object) const {
88 object.attr("type", "extrudedtriangular3d");
89 { auto a = object.addTag("vert"); vertAxis->writeXML(a); }
90 { auto a = object.addTag("long_tran"); longTranMesh.writeXML(a); }
91}
92
93 static shared_ptr<Mesh> readExtrudedTriangularMesh3D(XMLReader& reader) {
94 shared_ptr<ExtrudedTriangularMesh3D> result = plask::make_shared<ExtrudedTriangularMesh3D>();
95
96 XMLReader::CheckTagDuplication dub_check;
97 for (int i = 0; i < 2; ++i) {
98 reader.requireTag();
99 std::string node = reader.getNodeName();
100 dub_check(std::string("<mesh>"), node);
101 if (node == "vert")
102 result->vertAxis = readMeshAxis(reader);
103 else if (node == "long_tran")
104 result->longTranMesh = TriangularMesh2D::read(reader);
105 else
106 throw XMLUnexpectedElementException(reader, "<vert> or <long_tran>");
107 }
108 reader.requireTagEnd();
109 return result;
110}
111
112 static RegisterMeshReader extrudedtriangular3d_reader("extrudedtriangular3d", readExtrudedTriangularMesh3D);
113
114 Vec<3, double> ExtrudedTriangularMesh3D::at(std::size_t longTranIndex, std::size_t vertIndex) const {
115 return from_longTran_vert(longTranMesh[longTranIndex], vertAxis->at(vertIndex));
116}
117
118std::pair<std::size_t, std::size_t> ExtrudedTriangularMesh3D::longTranAndVertIndices(std::size_t index) const {
119 if (vertFastest) {
120 const std::size_t seg_size = vertAxis->size();
121 return std::pair<std::size_t, std::size_t>(index / seg_size, index % seg_size);
122 } else {
123 const std::size_t seg_size = longTranMesh.size();
124 return std::pair<std::size_t, std::size_t>(index % seg_size, index / seg_size);
125 }
126}
127
128std::size_t ExtrudedTriangularMesh3D::vertIndex(std::size_t index) const {
129 return vertFastest ? index % vertAxis->size() : index / longTranMesh.size();
130}
131
132 bool ExtrudedTriangularMesh3D::operator==(const ExtrudedTriangularMesh3D &to_compare) const {
133 if (empty()) return to_compare.empty();
134 // this mesh is not empty here
135 if (*vertAxis != *to_compare.vertAxis || longTranMesh != to_compare.longTranMesh) return false;
136 // here axes are equal
137 return vertFastest == to_compare.vertFastest || vertAxis->size() == 1 || longTranMesh.size() == 1;
138}
139
140 bool ExtrudedTriangularMesh3D::hasSameNodes(const MeshD<3> &to_compare) const {
141 if (const ExtrudedTriangularMesh3D* c = dynamic_cast<const ExtrudedTriangularMesh3D*>(&to_compare))
142 return *this == *c; // this will call == operator from ExtrudedTriangularMesh3D
143 return MeshD<3>::hasSameNodes(to_compare);
144}
145
146 TriangularMesh2D::SegmentsCounts ExtrudedTriangularMesh3D::countSegmentsIn(
147 std::size_t layer,
148 const GeometryD<3> &geometry,
149 const GeometryObject &object,
150 const PathHints *path) const
151{
152 TriangularMesh2D::SegmentsCounts result;
153 for (const auto el: this->longTranMesh.elements())
154 if (geometry.objectIncludes(object, path, this->at(el.getNodeIndex(0), layer)) &&
155 geometry.objectIncludes(object, path, this->at(el.getNodeIndex(1), layer)) &&
156 geometry.objectIncludes(object, path, this->at(el.getNodeIndex(2), layer)))
157 this->longTranMesh.countSegmentsOf(result, el);
158 return result;
159}
160
161 template<ExtrudedTriangularMesh3D::SideBoundaryDir boundaryDir>
162std::set<std::size_t> ExtrudedTriangularMesh3D::boundaryNodes(
163 const ExtrudedTriangularMesh3D::LayersIntervalSet& layers,
164 const GeometryD<3>& geometry,
165 const GeometryObject& object,
166 const PathHints *path) const
167{
168 std::set<std::size_t> result;
169 for (ExtrudedTriangularMesh3D::LayersInterval layer_interval: layers) {
170 for (std::size_t layer = layer_interval.lower(); layer < layer_interval.upper(); ++layer) {
171 for (std::size_t longTranNode: this->longTranMesh.boundaryNodes<ExtrudedTriangularMesh3D::boundaryDir3Dto2D(boundaryDir)>(countSegmentsIn(layer, geometry, object, path)))
172 result.insert(index(longTranNode, layer));
173 }
174 }
175 return result;
176}
177
178ExtrudedTriangularMesh3D::LayersInterval ExtrudedTriangularMesh3D::layersIn(const Box3D &box) const {
179 return LayersInterval(this->vertAxis->findIndex(box.lower.vert()), this->vertAxis->findUpIndex(box.upper.vert()));
180}
181
182ExtrudedTriangularMesh3D::LayersIntervalSet ExtrudedTriangularMesh3D::layersIn(const std::vector<Box3D>& boxes) const {
183 LayersIntervalSet layers;
184 for (const Box3D& box: boxes) {
185 LayersInterval interval = layersIn(box);
186 if (interval.lower() < interval.upper()) // if interval is not empty
187 layers.add(interval);
188 }
189 return layers;
190}
191
192 template <ExtrudedTriangularMesh3D::SideBoundaryDir boundaryDir>
193 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getObjBoundary(shared_ptr<const GeometryObject> object, const PathHints &path) {
194 return Boundary( [=](const ExtrudedTriangularMesh3D& mesh, const shared_ptr<const GeometryD<3>>& geometry) {
195 if (mesh.empty()) return BoundaryNodeSet(new EmptyBoundaryImpl());
196 LayersIntervalSet layers = mesh.layersIn(geometry->getObjectBoundingBoxes(object, path));
197 if (layers.empty()) return BoundaryNodeSet(new EmptyBoundaryImpl());
198 return BoundaryNodeSet(new StdSetBoundaryImpl(mesh.boundaryNodes<boundaryDir>(layers, *geometry, *object, &path)));
199 } );
200}
201
202 template <ExtrudedTriangularMesh3D::SideBoundaryDir boundaryDir>
203 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getObjBoundary(shared_ptr<const GeometryObject> object) {
204 return Boundary( [=](const ExtrudedTriangularMesh3D& mesh, const shared_ptr<const GeometryD<3>>& geometry) {
205 if (mesh.empty()) return BoundaryNodeSet(new EmptyBoundaryImpl());
206 LayersIntervalSet layers = mesh.layersIn(geometry->getObjectBoundingBoxes(object));
207 if (layers.empty()) return BoundaryNodeSet(new EmptyBoundaryImpl());
208 return BoundaryNodeSet(new StdSetBoundaryImpl(mesh.boundaryNodes<boundaryDir>(layers, *geometry, *object)));
209 } );
210}
211
212 template<ExtrudedTriangularMesh3D::SideBoundaryDir boundaryDir>
213 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getMeshBoundary()
214{
215 return Boundary( [=](const ExtrudedTriangularMesh3D& mesh, const shared_ptr<const GeometryD<3>>&) {
216 if (mesh.empty()) return BoundaryNodeSet(new EmptyBoundaryImpl());
217 return BoundaryNodeSet(new ExtrudedTriangularBoundaryImpl(
218 mesh,
219 mesh.longTranMesh.boundaryNodes<ExtrudedTriangularMesh3D::boundaryDir3Dto2D(boundaryDir)>(mesh.longTranMesh.countSegments()),
220 LayersInterval(0, mesh.vertAxis->size()-1)));
221 } );
222}
223
224 template<ExtrudedTriangularMesh3D::SideBoundaryDir boundaryDir>
225 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getBoxBoundary(const Box3D &box)
226{
227 return Boundary( [=](const ExtrudedTriangularMesh3D& mesh, const shared_ptr<const GeometryD<3>>&) {
228 if (mesh.empty()) return BoundaryNodeSet(new EmptyBoundaryImpl());
229 LayersInterval layers = mesh.layersIn(box);
230 if (layers.lower() >= layers.upper()) return BoundaryNodeSet(new EmptyBoundaryImpl());
231 return BoundaryNodeSet(new ExtrudedTriangularBoundaryImpl(
232 mesh,
233 mesh.longTranMesh.boundaryNodes<ExtrudedTriangularMesh3D::boundaryDir3Dto2D(boundaryDir)>(mesh.longTranMesh.countSegmentsIn(to_longTran(box))),
234 layers));
235 } );
236}
237
238BoundaryNodeSet ExtrudedTriangularMesh3D::topOrBottomBoundaryNodeSet(const Box3D &box, bool top) const {
239 LayersInterval layers = layersIn(box);
240 if (layers.lower() >= layers.upper()) return new EmptyBoundaryImpl();
241 const std::size_t layer = top ? layers.upper()-1 : layers.lower();
242 std::set<std::size_t> nodes3d;
243 Box2D box2d = to_longTran(box);
244 for (std::size_t index2d = 0; index2d < longTranMesh.size(); ++index2d)
245 if (box2d.contains(longTranMesh[index2d]))
246 nodes3d.insert(index(index2d, layer));
247 return new StdSetBoundaryImpl(std::move(nodes3d));
248}
249
250BoundaryNodeSet ExtrudedTriangularMesh3D::topOrBottomBoundaryNodeSet(const GeometryD<3>& geometry, const GeometryObject& object, const PathHints *path, bool top) const {
251 if (this->empty()) return BoundaryNodeSet(new EmptyBoundaryImpl());
252 LayersIntervalSet layers = layersIn(geometry.getObjectBoundingBoxes(object, path));
253 if (layers.empty()) return new EmptyBoundaryImpl();
254 std::unique_ptr<std::size_t[]> index2d_to_layer(new std::size_t[this->longTranMesh.size()]);
255 std::fill_n(index2d_to_layer.get(), this->longTranMesh.size(), std::numeric_limits<std::size_t>::max());
256 if (top) {
257 for (ExtrudedTriangularMesh3D::LayersInterval layer_interval: layers)
258 for (std::size_t layer = layer_interval.lower(); layer < layer_interval.upper(); ++layer)
259 for (std::size_t node2d_index = 0; node2d_index < longTranMesh.size(); ++node2d_index)
260 if (geometry.objectIncludes(object, path, at(node2d_index, layer)))
261 index2d_to_layer[node2d_index] = layer;
262 } else {
263 for (auto layers_it = layers.rbegin(); layers_it != layers.rend(); ++layers_it)
264 for (std::size_t layer = layers_it->upper(); layer-- > layers_it->lower(); )
265 for (std::size_t node2d_index = 0; node2d_index < longTranMesh.size(); ++node2d_index)
266 if (geometry.objectIncludes(object, path, at(node2d_index, layer)))
267 index2d_to_layer[node2d_index] = layer;
268 }
269 std::set<std::size_t> nodes3d;
270 for (std::size_t node2d_index = 0; node2d_index < longTranMesh.size(); ++node2d_index)
271 if (index2d_to_layer[node2d_index] != std::numeric_limits<std::size_t>::max())
272 nodes3d.insert(index(node2d_index, index2d_to_layer[node2d_index]));
273 return new StdSetBoundaryImpl(std::move(nodes3d));
274}
275
276 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getBackBoundary() {
277 return getMeshBoundary<SideBoundaryDir::BACK>();
278}
279
280 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getFrontBoundary() {
281 return getMeshBoundary<SideBoundaryDir::FRONT>();
282}
283
284 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getLeftBoundary() {
285 return getMeshBoundary<SideBoundaryDir::LEFT>();
286}
287
288 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getRightBoundary() {
289 return getMeshBoundary<SideBoundaryDir::RIGHT>();
290}
291
292 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getBottomBoundary() {
293 return Boundary( [](const ExtrudedTriangularMesh3D& mesh, const shared_ptr<const GeometryD<3>>&) {
294 if (mesh.empty()) return BoundaryNodeSet(new EmptyBoundaryImpl());
295 return BoundaryNodeSet(new ExtrudedTriangularWholeLayerBoundaryImpl(mesh, 0));
296 } );
297}
298
299 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getTopBoundary() {
300 return Boundary( [](const ExtrudedTriangularMesh3D& mesh, const shared_ptr<const GeometryD<3>>&) {
301 if (mesh.empty()) return BoundaryNodeSet(new EmptyBoundaryImpl());
302 return BoundaryNodeSet(new ExtrudedTriangularWholeLayerBoundaryImpl(mesh, mesh.vertAxis->size()-1));
303 } );
304}
305
306 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getAllSidesBoundary() {
307 return getMeshBoundary<SideBoundaryDir::ALL>();
308}
309
310 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getBackOfBoundary(const Box3D &box) {
311 return getBoxBoundary<SideBoundaryDir::BACK>(box);
312}
313
314 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getFrontOfBoundary(const Box3D &box) {
315 return getBoxBoundary<SideBoundaryDir::FRONT>(box);
316}
317
318 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getLeftOfBoundary(const Box3D &box) {
319 return getBoxBoundary<SideBoundaryDir::LEFT>(box);
320}
321
322 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getRightOfBoundary(const Box3D &box) {
323 return getBoxBoundary<SideBoundaryDir::RIGHT>(box);
324}
325
326 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getBottomOfBoundary(const Box3D &box) {
327 return Boundary( [=](const ExtrudedTriangularMesh3D& mesh, const shared_ptr<const GeometryD<3>>&) {
328 return mesh.topOrBottomBoundaryNodeSet(box, false);
329 } );
330}
331
332 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getTopOfBoundary(const Box3D &box) {
333 return Boundary( [=](const ExtrudedTriangularMesh3D& mesh, const shared_ptr<const GeometryD<3>>&) {
334 return mesh.topOrBottomBoundaryNodeSet(box, true);
335 } );
336}
337
338 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getAllSidesOfBoundary(const Box3D &box) {
339 return getBoxBoundary<SideBoundaryDir::ALL>(box);
340}
341
342
343 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getBackOfBoundary(shared_ptr<const GeometryObject> object, const PathHints &path) {
344 return getObjBoundary<SideBoundaryDir::BACK>(object, path);
345}
346 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getBackOfBoundary(shared_ptr<const GeometryObject> object) {
347 return getObjBoundary<SideBoundaryDir::BACK>(object);
348}
349
350 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getFrontOfBoundary(shared_ptr<const GeometryObject> object, const PathHints &path) {
351 return getObjBoundary<SideBoundaryDir::FRONT>(object, path);
352}
353 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getFrontOfBoundary(shared_ptr<const GeometryObject> object) {
354 return getObjBoundary<SideBoundaryDir::FRONT>(object);
355}
356
357 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getLeftOfBoundary(shared_ptr<const GeometryObject> object, const PathHints &path) {
358 return getObjBoundary<SideBoundaryDir::LEFT>(object, path);
359}
360 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getLeftOfBoundary(shared_ptr<const GeometryObject> object) {
361 return getObjBoundary<SideBoundaryDir::LEFT>(object);
362}
363
364 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getRightOfBoundary(shared_ptr<const GeometryObject> object, const PathHints &path) {
365 return getObjBoundary<SideBoundaryDir::RIGHT>(object, path);
366}
367 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getRightOfBoundary(shared_ptr<const GeometryObject> object) {
368 return getObjBoundary<SideBoundaryDir::RIGHT>(object);
369}
370
371 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getAllSidesBoundaryIn(shared_ptr<const GeometryObject> object, const PathHints &path) {
372 return getObjBoundary<SideBoundaryDir::ALL>(object, path);
373}
374 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getAllSidesBoundaryIn(shared_ptr<const GeometryObject> object) {
375 return getObjBoundary<SideBoundaryDir::ALL>(object);
376}
377
378 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getTopOfBoundary(shared_ptr<const GeometryObject> object, const PathHints &path) {
379 return Boundary( [=](const ExtrudedTriangularMesh3D& mesh, const shared_ptr<const GeometryD<3>>& geometry) {
380 return mesh.topOrBottomBoundaryNodeSet(*geometry, *object, &path, true);
381 } );
382}
383
384 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getTopOfBoundary(shared_ptr<const GeometryObject> object) {
385 return Boundary( [=](const ExtrudedTriangularMesh3D& mesh, const shared_ptr<const GeometryD<3>>& geometry) {
386 return mesh.topOrBottomBoundaryNodeSet(*geometry, *object, nullptr, true);
387 } );
388}
389
390 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getBottomOfBoundary(shared_ptr<const GeometryObject> object, const PathHints &path) {
391 return Boundary( [=](const ExtrudedTriangularMesh3D& mesh, const shared_ptr<const GeometryD<3>>& geometry) {
392 return mesh.topOrBottomBoundaryNodeSet(*geometry, *object, &path, false);
393 } );
394}
395
396 ExtrudedTriangularMesh3D::Boundary ExtrudedTriangularMesh3D::getBottomOfBoundary(shared_ptr<const GeometryObject> object) {
397 return Boundary( [=](const ExtrudedTriangularMesh3D& mesh, const shared_ptr<const GeometryD<3>>& geometry) {
398 return mesh.topOrBottomBoundaryNodeSet(*geometry, *object, nullptr, false);
399 } );
400}
401
402
403// ------------------ Nearest Neighbor interpolation ---------------------
404
405 template<typename DstT, typename SrcT>
406 NearestNeighborExtrudedTriangularMesh3DLazyDataImpl<DstT, SrcT>::NearestNeighborExtrudedTriangularMesh3DLazyDataImpl(
407 const shared_ptr<const ExtrudedTriangularMesh3D> &src_mesh,
408 const DataVector<const SrcT> &src_vec,
409 const shared_ptr<const MeshD<3> > &dst_mesh,
410 const InterpolationFlags &flags)
411 : InterpolatedLazyDataImpl<DstT, ExtrudedTriangularMesh3D, const SrcT>(src_mesh, src_vec, dst_mesh, flags),
412 nodesIndex(boost::irange(std::size_t(0), src_mesh->size()),
413 typename RtreeOfTriangularMesh2DNodes::parameters_type(),
414 TriangularMesh2DGetterForRtree(&src_mesh->longTranMesh))
415{
416}
417
418 template <typename DstT, typename SrcT>
419 DstT NearestNeighborExtrudedTriangularMesh3DLazyDataImpl<DstT, SrcT>::at(std::size_t index) const
420{
421 const auto point = this->dst_mesh->at(index);
422 const auto wrapped_point = this->flags.wrap(point);
423 const auto wrapped_longTran = to_longTran(wrapped_point);
424 for (auto v: nodesIndex | boost::geometry::index::adaptors::queried(boost::geometry::index::nearest(wrapped_longTran, 1)))
425 return this->flags.postprocess(point,
426 this->src_vec[
427 this->src_mesh->index(
428 v, this->src_mesh->vertAxis->findNearestIndex(wrapped_point.vert())
429 )
430 ]
431 );
432 return DstT();
433}
434
435 template struct PLASK_API NearestNeighborExtrudedTriangularMesh3DLazyDataImpl<double, double>;
436 template struct PLASK_API NearestNeighborExtrudedTriangularMesh3DLazyDataImpl<dcomplex, dcomplex>;
437 template struct PLASK_API NearestNeighborExtrudedTriangularMesh3DLazyDataImpl<Vec<2,double>, Vec<2,double>>;
438 template struct PLASK_API NearestNeighborExtrudedTriangularMesh3DLazyDataImpl<Vec<2,dcomplex>, Vec<2,dcomplex>>;
439 template struct PLASK_API NearestNeighborExtrudedTriangularMesh3DLazyDataImpl<Vec<3,double>, Vec<3,double>>;
440 template struct PLASK_API NearestNeighborExtrudedTriangularMesh3DLazyDataImpl<Vec<3,dcomplex>, Vec<3,dcomplex>>;
441 template struct PLASK_API NearestNeighborExtrudedTriangularMesh3DLazyDataImpl<Tensor2<double>, Tensor2<double>>;
442 template struct PLASK_API NearestNeighborExtrudedTriangularMesh3DLazyDataImpl<Tensor2<dcomplex>, Tensor2<dcomplex>>;
443 template struct PLASK_API NearestNeighborExtrudedTriangularMesh3DLazyDataImpl<Tensor3<double>, Tensor3<double>>;
444 template struct PLASK_API NearestNeighborExtrudedTriangularMesh3DLazyDataImpl<Tensor3<dcomplex>, Tensor3<dcomplex>>;
445
446
447// ------------------ Barycentric / Linear interpolation ---------------------
448
449 template<typename DstT, typename SrcT>
450 BarycentricExtrudedTriangularMesh3DLazyDataImpl<DstT, SrcT>::BarycentricExtrudedTriangularMesh3DLazyDataImpl(
451 const shared_ptr<const ExtrudedTriangularMesh3D> &src_mesh,
452 const DataVector<const SrcT> &src_vec,
453 const shared_ptr<const MeshD<3> > &dst_mesh,
454 const InterpolationFlags &flags
455 )
456 : InterpolatedLazyDataImpl<DstT, ExtrudedTriangularMesh3D, const SrcT>(src_mesh, src_vec, dst_mesh, flags),
457 elementIndex(src_mesh->longTranMesh)
458{
459}
460
461 template <typename DstT, typename SrcT>
462 DstT BarycentricExtrudedTriangularMesh3DLazyDataImpl<DstT, SrcT>::at(std::size_t index) const {
463 const auto point = this->dst_mesh->at(index);
464 const auto wrapped_point = this->flags.wrap(point);
465 const auto wrapped_longTran = to_longTran(wrapped_point);
466
467 for (auto v: elementIndex.rtree | boost::geometry::index::adaptors::queried(boost::geometry::index::intersects(wrapped_longTran))) {
468 const auto el = this->src_mesh->longTranMesh.getElement(v.second);
469 const auto b = el.barycentric(wrapped_longTran);
470 if (b.c0 < 0.0 || b.c1 < 0.0 || b.c2 < 0.0) continue; // wrapped_longTran is outside of the triangle
471
472 const std::size_t lonTranIndex0 = el.getNodeIndex(0),
473 lonTranIndex1 = el.getNodeIndex(1),
474 lonTranIndex2 = el.getNodeIndex(2);
475
476 size_t index_vert_lo, index_vert_hi;
477 double vert_lo, vert_hi;
478 bool invert_vert_lo, invert_vert_hi;
479 prepareInterpolationForAxis(
480 *this->src_mesh->vertAxis, this->flags,
481 wrapped_point.vert(), 2 /*index of vert axis*/,
482 index_vert_lo, index_vert_hi,
483 vert_lo, vert_hi,
484 invert_vert_lo, invert_vert_hi);
485
486 typename std::remove_const<typename std::remove_reference<decltype(this->src_vec[0])>::type>::type
487 data_lo = b.c0 * this->src_vec[this->src_mesh->index(lonTranIndex0, index_vert_lo)] +
488 b.c1 * this->src_vec[this->src_mesh->index(lonTranIndex1, index_vert_lo)] +
489 b.c2 * this->src_vec[this->src_mesh->index(lonTranIndex2, index_vert_lo)],
490 data_up = b.c0 * this->src_vec[this->src_mesh->index(lonTranIndex0, index_vert_hi)] +
491 b.c1 * this->src_vec[this->src_mesh->index(lonTranIndex1, index_vert_hi)] +
492 b.c2 * this->src_vec[this->src_mesh->index(lonTranIndex2, index_vert_hi)];
493
494 if (invert_vert_lo) data_lo = this->flags.reflect(2, data_lo);
495 if (invert_vert_hi) data_up = this->flags.reflect(2, data_up);
496
497 return this->flags.postprocess(point, interpolation::linear(vert_lo, data_lo, vert_hi, data_up, wrapped_point.vert()));
498 }
499 return NaN<decltype(this->src_vec[0])>();
500}
501
502 template struct PLASK_API BarycentricExtrudedTriangularMesh3DLazyDataImpl<double, double>;
503 template struct PLASK_API BarycentricExtrudedTriangularMesh3DLazyDataImpl<dcomplex, dcomplex>;
504 template struct PLASK_API BarycentricExtrudedTriangularMesh3DLazyDataImpl<Vec<2,double>, Vec<2,double>>;
505 template struct PLASK_API BarycentricExtrudedTriangularMesh3DLazyDataImpl<Vec<2,dcomplex>, Vec<2,dcomplex>>;
506 template struct PLASK_API BarycentricExtrudedTriangularMesh3DLazyDataImpl<Vec<3,double>, Vec<3,double>>;
507 template struct PLASK_API BarycentricExtrudedTriangularMesh3DLazyDataImpl<Vec<3,dcomplex>, Vec<3,dcomplex>>;
508 template struct PLASK_API BarycentricExtrudedTriangularMesh3DLazyDataImpl<Tensor2<double>, Tensor2<double>>;
509 template struct PLASK_API BarycentricExtrudedTriangularMesh3DLazyDataImpl<Tensor2<dcomplex>, Tensor2<dcomplex>>;
510 template struct PLASK_API BarycentricExtrudedTriangularMesh3DLazyDataImpl<Tensor3<double>, Tensor3<double>>;
511 template struct PLASK_API BarycentricExtrudedTriangularMesh3DLazyDataImpl<Tensor3<dcomplex>, Tensor3<dcomplex>>;
512
513
514// ------------------ Element mesh Nearest Neighbor interpolation ---------------------
515
516 template<typename DstT, typename SrcT>
517 NearestNeighborElementExtrudedTriangularMesh3DLazyDataImpl<DstT, SrcT>::NearestNeighborElementExtrudedTriangularMesh3DLazyDataImpl(
518 const shared_ptr<const ExtrudedTriangularMesh3D::ElementMesh> &src_mesh,
519 const DataVector<const SrcT> &src_vec,
520 const shared_ptr<const MeshD<3> > &dst_mesh,
521 const InterpolationFlags &flags)
522 : InterpolatedLazyDataImpl<DstT, ExtrudedTriangularMesh3D::ElementMesh, const SrcT>(src_mesh, src_vec, dst_mesh, flags),
523 elementIndex(src_mesh->getOriginalMesh().longTranMesh)
524{
525}
526
527 template<typename DstT, typename SrcT>
528 DstT NearestNeighborElementExtrudedTriangularMesh3DLazyDataImpl<DstT, SrcT>::at(std::size_t index) const {
529 const auto point = this->dst_mesh->at(index);
530 const auto wrapped_point = this->flags.wrap(point);
531
532 const ExtrudedTriangularMesh3D& orginal_src_mesh = this->src_mesh->getOriginalMesh();
533 const MeshAxis& vertAxis = *orginal_src_mesh.vertAxis;
534 if (wrapped_point.vert() < vertAxis[0] || vertAxis[vertAxis.size()] < wrapped_point.vert())
535 return NaN<decltype(this->src_vec[0])>();
536
537 const auto wrapped_longTran = to_longTran(wrapped_point);
538
539 std::size_t longTran_element_index = this->elementIndex.getIndex(wrapped_longTran);
540 if (longTran_element_index == TriangularMesh2D::ElementIndex::INDEX_NOT_FOUND)
541 return NaN<decltype(this->src_vec[0])>();
542
543 return this->flags.postprocess(point,
544 this->src_vec[orginal_src_mesh.elementIndex(
545 longTran_element_index,
546 vertAxis.findUpIndex(wrapped_point.vert())-1
547 )]);
548}
549
550 template struct PLASK_API NearestNeighborElementExtrudedTriangularMesh3DLazyDataImpl<double, double>;
551 template struct PLASK_API NearestNeighborElementExtrudedTriangularMesh3DLazyDataImpl<dcomplex, dcomplex>;
552 template struct PLASK_API NearestNeighborElementExtrudedTriangularMesh3DLazyDataImpl<Vec<2,double>, Vec<2,double>>;
553 template struct PLASK_API NearestNeighborElementExtrudedTriangularMesh3DLazyDataImpl<Vec<2,dcomplex>, Vec<2,dcomplex>>;
554 template struct PLASK_API NearestNeighborElementExtrudedTriangularMesh3DLazyDataImpl<Vec<3,double>, Vec<3,double>>;
555 template struct PLASK_API NearestNeighborElementExtrudedTriangularMesh3DLazyDataImpl<Vec<3,dcomplex>, Vec<3,dcomplex>>;
556 template struct PLASK_API NearestNeighborElementExtrudedTriangularMesh3DLazyDataImpl<Tensor2<double>, Tensor2<double>>;
557 template struct PLASK_API NearestNeighborElementExtrudedTriangularMesh3DLazyDataImpl<Tensor2<dcomplex>, Tensor2<dcomplex>>;
558 template struct PLASK_API NearestNeighborElementExtrudedTriangularMesh3DLazyDataImpl<Tensor3<double>, Tensor3<double>>;
559 template struct PLASK_API NearestNeighborElementExtrudedTriangularMesh3DLazyDataImpl<Tensor3<dcomplex>, Tensor3<dcomplex>>;
560
561 bool ExtrudedTriangularMesh3D::ElementMesh::hasSameNodes(const MeshD<3> &to_compare) const {
562 if (const ElementMesh* c = dynamic_cast<const ElementMesh*>(&to_compare))
563 if (this->originalMesh == c->originalMesh) return true;
564 return MeshD<3>::hasSameNodes(to_compare);
565}
566
567
568
569
570
571
572
573
574} // namespace plask