PLaSK library

Loading...
Searching...
No Matches
rectangular_masked3d.hpp
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#ifndef PLASK__RECTANGULAR_MASKED3D_H
15#define PLASK__RECTANGULAR_MASKED3D_H
16
17#include "rectangular_masked_common.hpp"
18
19 namespace plask {
20
28 struct PLASK_API RectangularMaskedMesh3D: public RectangularMaskedMeshBase<3> {
29
35 inline std::size_t index2(std::size_t mesh_index) const { // method missing in the base as it is specific for 3D
36 return fullMesh.index2(nodeSet.at(mesh_index));
37 }
38
44 inline std::size_t middleIndex(std::size_t mesh_index) const { // method missing in the base as it is specific for 3D
45 return fullMesh.middleIndex(nodeSet.at(mesh_index));
46 }
47
52 std::size_t getElementsCount2() const { // method missing in the base as it is specific for 3D
53 return fullMesh.getElementsCount2();
54 }
55
61 double getElementMidpoint2(std::size_t index2) const { // method missing in the base as it is specific for 3D
62 return fullMesh.getElementMidpoint2(index2);
63 }
64
66 typedef std::function<bool(const RectangularMesh3D::Element&)> Predicate;
67
68 class PLASK_API Element {
69
70 const RectangularMaskedMesh3D& maskedMesh;
71
72 //std::uint32_t elementNumber; ///< index of element in original mesh
73 std::size_t index0, index1, index2; // probably this form allows to do most operation fastest in average, low indexes of element corner or just element indexes
74
76 mutable std::size_t elementIndex;
77
78 const RectangularMesh<3>& fullMesh() const { return maskedMesh.fullMesh; }
79
80 public:
81
82 enum: std::size_t { UNKNOWN_ELEMENT_INDEX = std::numeric_limits<std::size_t>::max() };
83
84 Element(const RectangularMaskedMesh3D& maskedMesh, std::size_t elementIndex, std::size_t index0, std::size_t index1, std::size_t index2)
85 : maskedMesh(maskedMesh), index0(index0), index1(index1), index2(index2), elementIndex(elementIndex)
86 {
87 }
88
89 Element(const RectangularMaskedMesh3D& maskedMesh, std::size_t elementIndex, std::size_t elementIndexOfFullMesh)
90 : maskedMesh(maskedMesh), elementIndex(elementIndex)
91 {
92 const std::size_t v = fullMesh().getElementMeshLowIndex(elementIndexOfFullMesh);
93 index0 = fullMesh().index0(v);
94 index1 = fullMesh().index1(v);
95 index2 = fullMesh().index2(v);
96 }
97
98 Element(const RectangularMaskedMesh3D& maskedMesh, std::size_t elementIndex)
99 : Element(maskedMesh, elementIndex, maskedMesh.elementSet.at(elementIndex))
100 {}
101
102
104 inline std::size_t getIndex0() const { return index0; }
105
107 inline std::size_t getIndex1() const { return index1; }
108
110 inline std::size_t getIndex2() const { return index2; }
111
113 inline std::size_t getLowerIndex0() const { return index0; }
114
116 inline std::size_t getLowerIndex1() const { return index1; }
117
119 inline std::size_t getLowerIndex2() const { return index2; }
120
122 inline double getLower0() const { return fullMesh().axis[0]->at(index0); }
123
125 inline double getLower1() const { return fullMesh().axis[1]->at(index1); }
126
128 inline double getLower2() const { return fullMesh().axis[2]->at(index2); }
129
131 inline std::size_t getUpperIndex0() const { return index0+1; }
132
134 inline std::size_t getUpperIndex1() const { return index1+1; }
135
137 inline std::size_t getUpperIndex2() const { return index2+1; }
138
140 inline double getUpper0() const { return fullMesh().axis[0]->at(getUpperIndex0()); }
141
143 inline double getUpper1() const { return fullMesh().axis[1]->at(getUpperIndex1()); }
144
146 inline double getUpper2() const { return fullMesh().axis[2]->at(getUpperIndex2()); }
147
149 inline double getSize0() const { return getUpper0() - getLower0(); }
150
152 inline double getSize1() const { return getUpper1() - getLower1(); }
153
155 inline double getSize2() const { return getUpper2() - getLower2(); }
156
158 inline Vec<3, double> getSize() const { return getUpUpUp() - getLoLoLo(); }
159
161 inline Vec<3, double> getMidpoint() const { return maskedMesh.getElementMidpoint(index0, index1, index2); }
162
164 inline std::size_t getIndex() const {
165 if (elementIndex == UNKNOWN_ELEMENT_INDEX)
166 elementIndex = maskedMesh.getElementIndexFromLowIndexes(getLowerIndex0(), getLowerIndex1(), getLowerIndex2());
167 return elementIndex;
168 }
169
171 inline Box3D toBox() const { return maskedMesh.getElementBox(index0, index1, index2); }
172
174 inline double getVolume() const { return getSize0() * getSize1() * getSize2(); }
175
177 inline double getArea() const { return getVolume(); }
178
180 inline std::size_t getLoLoLoIndex() const { return maskedMesh.index(getLowerIndex0(), getLowerIndex1(), getLowerIndex2()); }
181
183 inline std::size_t getUpLoLoIndex() const { return maskedMesh.index(getUpperIndex0(), getLowerIndex1(), getLowerIndex2()); }
184
186 inline std::size_t getLoUpLoIndex() const { return maskedMesh.index(getLowerIndex0(), getUpperIndex1(), getLowerIndex2()); }
187
189 inline std::size_t getUpUpLoIndex() const { return maskedMesh.index(getUpperIndex0(), getUpperIndex1(), getLowerIndex2()); }
190
192 inline std::size_t getLoLoUpIndex() const { return maskedMesh.index(getLowerIndex0(), getLowerIndex1(), getUpperIndex2()); }
193
195 inline std::size_t getUpLoUpIndex() const { return maskedMesh.index(getUpperIndex0(), getLowerIndex1(), getUpperIndex2()); }
196
198 inline std::size_t getLoUpUpIndex() const { return maskedMesh.index(getLowerIndex0(), getUpperIndex1(), getUpperIndex2()); }
199
201 inline std::size_t getUpUpUpIndex() const { return maskedMesh.index(getUpperIndex0(), getUpperIndex1(), getUpperIndex2()); }
202
204 inline Vec<3, double> getLoLoLo() const { return maskedMesh(getLowerIndex0(), getLowerIndex1(), getLowerIndex2()); }
205
207 inline Vec<3, double> getUpLoLo() const { return maskedMesh(getUpperIndex0(), getLowerIndex1(), getLowerIndex2()); }
208
210 inline Vec<3, double> getLoUpLo() const { return maskedMesh(getLowerIndex0(), getUpperIndex1(), getLowerIndex2()); }
211
213 inline Vec<3, double> getUpUpLo() const { return maskedMesh(getUpperIndex0(), getUpperIndex1(), getLowerIndex2()); }
214
216 inline Vec<3, double> getLoLoUp() const { return maskedMesh(getLowerIndex0(), getLowerIndex1(), getUpperIndex2()); }
217
219 inline Vec<3, double> getUpLoUp() const { return maskedMesh(getUpperIndex0(), getLowerIndex1(), getUpperIndex2()); }
220
222 inline Vec<3, double> getLoUpUp() const { return maskedMesh(getLowerIndex0(), getUpperIndex1(), getUpperIndex2()); }
223
225 inline Vec<3, double> getUpUpUp() const { return maskedMesh(getUpperIndex0(), getUpperIndex1(), getUpperIndex2()); }
226
227 }; // class Element
228
229 struct PLASK_API Elements: ElementsBase<RectangularMaskedMesh3D> {
230
231 explicit Elements(const RectangularMaskedMesh3D& mesh): ElementsBase(mesh) { mesh.ensureHasElements(); }
232
233 Element operator()(std::size_t i0, std::size_t i1, std::size_t i2) const { return Element(*maskedMesh, Element::UNKNOWN_ELEMENT_INDEX, i0, i1, i2); }
234
235 }; // struct Elements
236
238 struct PLASK_API ElementMesh: ElementMeshBase<RectangularMaskedMesh3D> {
239
240 explicit ElementMesh(const RectangularMaskedMesh3D* originalMesh): ElementMeshBase<RectangularMaskedMesh3D>(originalMesh) {}
241
242 // Convert to recctangular masked mesh
243 RectangularMaskedMesh3D toMasked() const {
244 return RectangularMaskedMesh3D(fullMesh, originalMesh->elementSet);
245 }
246
247 // Convert to recctangular masked mesh
248 operator RectangularMaskedMesh3D() const {
249 return RectangularMaskedMesh3D(fullMesh, originalMesh->elementSet);
250 }
251
259 inline std::size_t index(std::size_t axis0_index, std::size_t axis1_index, std::size_t axis2_index) const {
260 return originalMesh->elementSet.indexOf(fullMesh.index(axis0_index, axis1_index, axis2_index));
261 }
262
263 bool prepareInterpolation(const Vec<3>& point, Vec<3>& wrapped_point,
264 std::size_t& index0_lo, std::size_t& index0_hi,
265 std::size_t& index1_lo, std::size_t& index1_hi,
266 std::size_t& index2_lo, std::size_t& index2_hi,
267 const InterpolationFlags& flags) const {
268 return originalMesh->prepareInterpolation(point, wrapped_point, index0_lo, index0_hi, index1_lo, index1_hi, index2_lo, index2_hi, flags);
269 }
270
277 template <typename RandomAccessContainer>
278 auto interpolateLinear(const RandomAccessContainer& data, const Vec<3>& point, const InterpolationFlags& flags) const
279 -> typename std::remove_reference<decltype(data[0])>::type {
280 typedef typename std::remove_reference<decltype(data[0])>::type DataT;
281 Vec<3> p;
282 size_t index0, index0_hi, index1, index1_hi, index2, index2_hi;
283
284 if (!prepareInterpolation(point, p, index0, index0_hi, index1, index1_hi, index2, index2_hi, flags))
285 return NaN<decltype(data[0])>();
286
287 Vec<3> pa = fullMesh.at(index0, index1, index2);
288
289 size_t step0 = (p.c0 < pa.c0)?
290 (index0 == 0)? 0 : -1 :
291 (index0_hi == fullMesh.axis[0]->size())? 0 : 1;
292 size_t step1 = (p.c1 < pa.c1)?
293 (index1 == 0)? 0 : -1 :
294 (index1_hi == fullMesh.axis[1]->size())? 0 : 1;
295 size_t step2 = (p.c2 < pa.c2)?
296 (index2 == 0)? 0 : -1 :
297 (index2_hi == fullMesh.axis[2]->size())? 0 : 1;
298
299 size_t index_aaa = index(index0, index1, index2), index_aab, index_aba, index_abb,
300 index_baa, index_bab, index_bba, index_bbb;
301
302 typename std::remove_const<DataT>::type data_aaa = data[index_aaa], data_aab, data_aba, data_abb,
303 data_baa, data_bab, data_bba, data_bbb;
304
305 if (step0 == 0 && step1 == 0 && step2 == 0) {
306 index_aab = index_aba = index_abb = index_baa = index_bab = index_bba = index_bbb = index_aaa;
307 data_aab = data_aba = data_abb = data_baa = data_bab = data_bba = data_bbb = data_aaa;
308 } else {
309 index_aab = index(index0, index1, index2+step2);
310 index_aba = index(index0, index1+step1, index2);
311 index_abb = index(index0, index1+step1, index2+step2);
312 index_baa = index(index0+step0, index1, index2);
313 index_bab = index(index0+step0, index1, index2+step2);
314 index_bba = index(index0+step0, index1+step1, index2);
315 index_bbb = index(index0+step0, index1+step1, index2+step2);
316 data_aab = (index_aab != Element::UNKNOWN_ELEMENT_INDEX)? data[index_aab] : data_aaa;
317 data_aba = (index_aba != Element::UNKNOWN_ELEMENT_INDEX)? data[index_aba] : data_aaa;
318 data_baa = (index_baa != Element::UNKNOWN_ELEMENT_INDEX)? data[index_baa] : data_aaa;
319 data_abb = (index_abb != Element::UNKNOWN_ELEMENT_INDEX)? data[index_abb] : data_aab + data_aba - data_aaa;
320 data_bab = (index_bab != Element::UNKNOWN_ELEMENT_INDEX)? data[index_bab] : data_baa + data_aab - data_aaa;
321 data_bba = (index_bba != Element::UNKNOWN_ELEMENT_INDEX)? data[index_bba] : data_aba + data_baa - data_aaa;
322 data_bbb = (index_bbb != Element::UNKNOWN_ELEMENT_INDEX)? data[index_bbb] : data_aab + data_aba + data_baa - 2. * data_aaa;
323 }
324
325 Vec<3> pb = fullMesh.at(index0+step0, index1+step1, index2+step2);
326 if (step0 == 0) pb.c0 += 1.;
327 if (step1 == 0) pb.c1 += 1.;
328 if (step2 == 0) pb.c2 += 1.;
329
330 return flags.postprocess(point,
331 interpolation::trilinear(pa.c0, pb.c0, pa.c1, pb.c1, pa.c2, pb.c2,
332 data_aaa, data_baa, data_bba, data_aba,
333 data_aab, data_bab, data_bbb, data_abb,
334 p.c0, p.c1, p.c2));
335 }
336
343 template <typename RandomAccessContainer>
344 auto interpolateNearestNeighbor(const RandomAccessContainer& data, const Vec<3>& point, const InterpolationFlags& flags) const
345 -> typename std::remove_reference<decltype(data[0])>::type {
346 Vec<3> wrapped_point;
347 std::size_t index0_lo, index0_hi, index1_lo, index1_hi, index2_lo, index2_hi;
348
349 if (!originalMesh->prepareInterpolation(point, wrapped_point, index0_lo, index0_hi, index1_lo, index1_hi, index2_lo, index2_hi, flags))
350 return NaN<decltype(data[0])>();
351
352 return flags.postprocess(point, data[this->index(index0_lo, index1_lo, index2_lo)]);
353 }
354
355 }; // struct ElementMesh
356
360 RectangularMaskedMesh3D() = default;
361
366 void reset(const Predicate& predicate);
367
375 RectangularMaskedMesh3D(const RectangularMesh<3>& fullMesh, const Predicate& predicate, bool clone_axes = false);
376
384 void reset(const RectangularMesh<3>& fullMesh, const Predicate& predicate, bool clone_axes = false);
385
394 RectangularMaskedMesh3D(const RectangularMesh<3>& fullMesh, const GeometryD<3>& geom, const std::function<bool(shared_ptr<const Material>)> materialPredicate, bool clone_axes = false)
395 : RectangularMaskedMesh3D(fullMesh,
396 [&](const RectangularMesh3D::Element& el) { return materialPredicate(geom.getMaterial(el.getMidpoint())); },
397 clone_axes)
398 {
399 }
400
409 void reset(const RectangularMesh<3>& rectangularMesh, const GeometryD<3>& geom, const std::function<bool(shared_ptr<const Material>)> materialPredicate, bool clone_axes = false) {
410 reset(rectangularMesh,
411 [&](const RectangularMesh3D::Element& el) { return materialPredicate(geom.getMaterial(el.getMidpoint())); },
412 clone_axes);
413 }
414
425 RectangularMaskedMesh3D(const RectangularMesh<3>& rectangularMesh, const GeometryD<3>& geom, unsigned materialKinds, bool clone_axes = false)
426 : RectangularMaskedMesh3D(rectangularMesh,
427 [&](const RectangularMesh3D::Element& el) { return (geom.getMaterial(el.getMidpoint())->kind() & materialKinds) != 0; },
428 clone_axes)
429 {
430 }
431
442 void reset(const RectangularMesh<3>& rectangularMesh, const GeometryD<3>& geom, unsigned materialKinds, bool clone_axes = false) {
443 reset(rectangularMesh,
444 [&](const RectangularMesh3D::Element& el) { return (geom.getMaterial(el.getMidpoint())->kind() & materialKinds) != 0; },
445 clone_axes);
446 }
447
456 RectangularMaskedMesh3D(const RectangularMesh<DIM>& rectangularMesh, Set nodeSet, bool clone_axes = false);
457
458 Elements elements() const { return Elements(*this); }
459 Elements getElements() const { return elements(); }
460
461 Element element(std::size_t i0, std::size_t i1, std::size_t i2) const { ensureHasElements(); return Element(*this, Element::UNKNOWN_ELEMENT_INDEX, i0, i1, i2); }
462 Element getElement(std::size_t i0, std::size_t i1, std::size_t i2) const { return element(i0, i1, i2); }
463
469 Element element(std::size_t i) const { ensureHasElements(); return Element(*this, i); }
470
476 Element getElement(std::size_t i) const { return element(i); }
477
485 inline std::size_t index(std::size_t axis0_index, std::size_t axis1_index, std::size_t axis2_index) const {
486 return nodeSet.indexOf(fullMesh.index(axis0_index, axis1_index, axis2_index));
487 }
488
489 using RectangularMaskedMeshBase<3>::index;
490 using RectangularMaskedMeshBase<3>::at;
491
499 inline Vec<3, double> at(std::size_t index0, std::size_t index1, std::size_t index2) const {
500 return fullMesh.at(index0, index1, index2);
501 }
502
510 inline Vec<3, double> operator()(std::size_t axis0_index, std::size_t axis1_index, std::size_t axis2_index) const {
511 return fullMesh.operator()(axis0_index, axis1_index, axis2_index);
512 }
513
518 shared_ptr<RectangularMaskedMesh3D::ElementMesh> getElementMesh() const {
519 return make_shared<RectangularMaskedMesh3D::ElementMesh>(this);
520 }
521
522 private:
523
524 void initNodesAndElements(const RectangularMaskedMesh3D::Predicate &predicate);
525
526 bool canBeIncluded(const Vec<3>& point) const {
527 return
528 fullMesh.axis[0]->at(0) - point[0] < MIN_DISTANCE && point[0] - fullMesh.axis[0]->at(fullMesh.axis[0]->size()-1) < MIN_DISTANCE &&
529 fullMesh.axis[1]->at(0) - point[1] < MIN_DISTANCE && point[1] - fullMesh.axis[1]->at(fullMesh.axis[1]->size()-1) < MIN_DISTANCE &&
530 fullMesh.axis[2]->at(0) - point[2] < MIN_DISTANCE && point[2] - fullMesh.axis[2]->at(fullMesh.axis[2]->size()-1) < MIN_DISTANCE;
531 }
532
533 public:
534
544 bool prepareInterpolation(const Vec<3>& point, Vec<3>& wrapped_point,
545 std::size_t& index0_lo, std::size_t& index0_hi,
546 std::size_t& index1_lo, std::size_t& index1_hi,
547 std::size_t& index2_lo, std::size_t& index2_hi,
548 const InterpolationFlags& flags) const;
555 template <typename RandomAccessContainer>
556 auto interpolateLinear(const RandomAccessContainer& data, const Vec<3>& point, const InterpolationFlags& flags) const
557 -> typename std::remove_reference<decltype(data[0])>::type
558 {
559 Vec<3> wrapped_point;
560 std::size_t index0_lo, index0_hi, index1_lo, index1_hi, index2_lo, index2_hi;
561
562 if (!prepareInterpolation(point, wrapped_point, index0_lo, index0_hi, index1_lo, index1_hi, index2_lo, index2_hi, flags))
563 return NaN<decltype(data[0])>();
564
565 return flags.postprocess(point,
566 interpolation::trilinear(
567 fullMesh.axis[0]->at(index0_lo), fullMesh.axis[0]->at(index0_hi),
568 fullMesh.axis[1]->at(index1_lo), fullMesh.axis[1]->at(index1_hi),
569 fullMesh.axis[2]->at(index2_lo), fullMesh.axis[2]->at(index2_hi),
570 data[index(index0_lo, index1_lo, index2_lo)],
571 data[index(index0_hi, index1_lo, index2_lo)],
572 data[index(index0_hi, index1_hi, index2_lo)],
573 data[index(index0_lo, index1_hi, index2_lo)],
574 data[index(index0_lo, index1_lo, index2_hi)],
575 data[index(index0_hi, index1_lo, index2_hi)],
576 data[index(index0_hi, index1_hi, index2_hi)],
577 data[index(index0_lo, index1_hi, index2_hi)],
578 wrapped_point.c0, wrapped_point.c1, wrapped_point.c2));
579 }
580
587 template <typename RandomAccessContainer>
588 auto interpolateNearestNeighbor(const RandomAccessContainer& data, const Vec<3>& point, const InterpolationFlags& flags) const
589 -> typename std::remove_reference<decltype(data[0])>::type
590 {
591 Vec<3> wrapped_point;
592 std::size_t index0_lo, index0_hi, index1_lo, index1_hi, index2_lo, index2_hi;
593
594 if (!prepareInterpolation(point, wrapped_point, index0_lo, index0_hi, index1_lo, index1_hi, index2_lo, index2_hi, flags))
595 return NaN<decltype(data[0])>();
596
597 return flags.postprocess(point,
598 data[this->index(
599 nearest(wrapped_point.c0, *fullMesh.axis[0], index0_lo, index0_hi),
600 nearest(wrapped_point.c1, *fullMesh.axis[1], index1_lo, index1_hi),
601 nearest(wrapped_point.c2, *fullMesh.axis[2], index2_lo, index2_hi)
602 )]);
603 }
604
612 std::size_t getElementIndexFromLowIndexes(std::size_t axis0_index, std::size_t axis1_index, std::size_t axis2_index) const {
613 return ensureHasElements().indexOf(fullMesh.getElementIndexFromLowIndexes(axis0_index, axis1_index, axis2_index));
614 }
615
621 double getElementArea(std::size_t index0, std::size_t index1, std::size_t index2) const {
622 return fullMesh.getElementArea(index0, index1, index2);
623 }
624
630 Vec<3, double> getElementMidpoint(std::size_t index0, std::size_t index1, std::size_t index2) const {
631 return fullMesh.getElementMidpoint(index0, index1, index2);
632 }
633
639 Box3D getElementBox(std::size_t index0, std::size_t index1, std::size_t index2) const {
640 return fullMesh.getElementBox(index0, index1, index2);
641 }
642
643
644 protected:
645
653 template <int CHANGE_DIR_SLOWER, int CHANGE_DIR_FASTER>
654 struct BoundaryIteratorImpl: public plask::BoundaryNodeSetImpl::IteratorImpl {
655
656 const RectangularMaskedMeshBase<3> &mesh;
657
659 Vec<3, std::size_t> index;
660
662 const std::size_t indexFasterBegin, indexFasterEnd, indexSlowerEnd;
663
664 private:
666 void naiveIncrement() {
667 ++index[CHANGE_DIR_FASTER];
668 if (index[CHANGE_DIR_FASTER] == indexFasterEnd) {
669 index[CHANGE_DIR_FASTER] = indexFasterBegin;
670 ++index[CHANGE_DIR_SLOWER];
671 }
672 }
673 public:
674 BoundaryIteratorImpl(const RectangularMaskedMeshBase<3>& mesh, Vec<3, std::size_t> index, std::size_t indexSlowerEnd, std::size_t indexFasterEnd)
675 : mesh(mesh), index(index), indexFasterBegin(index[CHANGE_DIR_FASTER]), indexFasterEnd(indexFasterEnd), indexSlowerEnd(indexSlowerEnd)
676 {
677 // go to the first index existed in order to make dereference possible:
678 while (this->index[CHANGE_DIR_SLOWER] < indexSlowerEnd && mesh.index(this->index) == NOT_INCLUDED)
679 naiveIncrement();
680 }
681
682 void increment() override {
683 do {
684 naiveIncrement();
685 } while (index[CHANGE_DIR_SLOWER] < indexSlowerEnd && mesh.index(index) == NOT_INCLUDED);
686 }
687
688 bool equal(const plask::BoundaryNodeSetImpl::IteratorImpl& other) const override {
689 return index == static_cast<const BoundaryIteratorImpl&>(other).index;
690 }
691
692 std::size_t dereference() const override {
693 return mesh.index(index);
694 }
695
696 std::unique_ptr<plask::BoundaryNodeSetImpl::IteratorImpl> clone() const override {
697 return std::unique_ptr<plask::BoundaryNodeSetImpl::IteratorImpl>(new BoundaryIteratorImpl<CHANGE_DIR_SLOWER, CHANGE_DIR_FASTER>(*this));
698 }
699
700 };
701
702 template <int CHANGE_DIR_SLOWER, int CHANGE_DIR_FASTER>
703 struct BoundaryNodeSetImpl: public BoundaryNodeSetWithMeshImpl<RectangularMaskedMeshBase<3>> {
704
705 static constexpr int FIXED_DIR = 3 - CHANGE_DIR_SLOWER - CHANGE_DIR_FASTER;
706
707 using typename BoundaryNodeSetWithMeshImpl<RectangularMaskedMeshBase<3>>::const_iterator;
708
710 Vec<3, std::size_t> index;
711
713 std::size_t indexFasterEnd, indexSlowerEnd;
714
715 BoundaryNodeSetImpl(const RectangularMaskedMeshBase<DIM>& mesh, Vec<3, std::size_t> index, std::size_t indexSlowerEnd, std::size_t indexFasterEnd)
716 : BoundaryNodeSetWithMeshImpl<RectangularMaskedMeshBase<3>>(mesh), index(index), indexFasterEnd(indexFasterEnd), indexSlowerEnd(indexSlowerEnd) {}
717
718 BoundaryNodeSetImpl(const RectangularMaskedMeshBase<DIM>& mesh, std::size_t index0, std::size_t index1, std::size_t index2, std::size_t indexSlowerEnd, std::size_t indexFasterEnd)
719 : BoundaryNodeSetWithMeshImpl<RectangularMaskedMeshBase<3>>(mesh), index(index0, index1, index2), indexFasterEnd(indexFasterEnd), indexSlowerEnd(indexSlowerEnd) {}
720
721 bool contains(std::size_t mesh_index) const override {
722 if (mesh_index >= this->mesh.size()) return false;
723 Vec<3, std::size_t> mesh_indexes = this->mesh.indexes(mesh_index);
724 return mesh_indexes[FIXED_DIR] == index[FIXED_DIR] &&
725 (index[CHANGE_DIR_FASTER] <= mesh_indexes[CHANGE_DIR_FASTER] && mesh_indexes[CHANGE_DIR_FASTER] < indexFasterEnd) &&
726 (index[CHANGE_DIR_SLOWER] <= mesh_indexes[CHANGE_DIR_SLOWER] && mesh_indexes[CHANGE_DIR_SLOWER] < indexSlowerEnd);
727 }
728
729 const_iterator begin() const override {
730 return Iterator(new BoundaryIteratorImpl<CHANGE_DIR_SLOWER, CHANGE_DIR_FASTER>(this->mesh, index, indexSlowerEnd, indexFasterEnd));
731 }
732
733 const_iterator end() const override {
734 Vec<3, std::size_t> index_end = index;
735 index_end[CHANGE_DIR_SLOWER] = indexSlowerEnd;
736 return Iterator(new BoundaryIteratorImpl<CHANGE_DIR_SLOWER, CHANGE_DIR_FASTER>(this->mesh, index_end, indexSlowerEnd, indexFasterEnd));
737 }
738 };
739
740 public: // boundaries:
741
742 BoundaryNodeSet createIndex0BoundaryAtLine(std::size_t line_nr_axis0,
743 std::size_t index1Begin, std::size_t index1End,
744 std::size_t index2Begin, std::size_t index2End
745 ) const override;
746
747 BoundaryNodeSet createIndex0BoundaryAtLine(std::size_t line_nr_axis0) const override;
748
749 BoundaryNodeSet createIndex1BoundaryAtLine(std::size_t line_nr_axis1,
750 std::size_t index0Begin, std::size_t index0End,
751 std::size_t index2Begin, std::size_t index2End
752 ) const override;
753
754 BoundaryNodeSet createIndex1BoundaryAtLine(std::size_t line_nr_axis1) const override;
755
756 BoundaryNodeSet createIndex2BoundaryAtLine(std::size_t line_nr_axis2,
757 std::size_t index0Begin, std::size_t index0End,
758 std::size_t index1Begin, std::size_t index1End
759 ) const override;
760
761 BoundaryNodeSet createIndex2BoundaryAtLine(std::size_t line_nr_axis2) const override;
762
763 BoundaryNodeSet createBackBoundary() const override;
764
765 BoundaryNodeSet createFrontBoundary() const override;
766
767 BoundaryNodeSet createLeftBoundary() const override;
768
769 BoundaryNodeSet createRightBoundary() const override;
770
771 BoundaryNodeSet createBottomBoundary() const override;
772
773 BoundaryNodeSet createTopBoundary() const override;
774
775 BoundaryNodeSet createBackOfBoundary(const Box3D& box) const override;
776
777 BoundaryNodeSet createFrontOfBoundary(const Box3D& box) const override;
778
779 BoundaryNodeSet createLeftOfBoundary(const Box3D& box) const override;
780
781 BoundaryNodeSet createRightOfBoundary(const Box3D& box) const override;
782
783 BoundaryNodeSet createBottomOfBoundary(const Box3D& box) const override;
784
785 BoundaryNodeSet createTopOfBoundary(const Box3D& box) const override;
786};
787
788 template <typename SrcT, typename DstT>
789 struct InterpolationAlgorithm<RectangularMaskedMesh3D, SrcT, DstT, INTERPOLATION_LINEAR> {
790 static LazyData<DstT> interpolate(const shared_ptr<const RectangularMaskedMesh3D>& src_mesh, const DataVector<const SrcT>& src_vec,
791 const shared_ptr<const MeshD<3>>& dst_mesh, const InterpolationFlags& flags) {
792 if (src_mesh->empty()) throw BadMesh("interpolate", "source mesh empty");
793 return new LinearInterpolatedLazyDataImpl<DstT, RectangularMaskedMesh3D, SrcT>(src_mesh, src_vec, dst_mesh, flags);
794 }
795};
796
797 template <typename SrcT, typename DstT>
798 struct InterpolationAlgorithm<RectangularMaskedMesh3D, SrcT, DstT, INTERPOLATION_NEAREST> {
799 static LazyData<DstT> interpolate(const shared_ptr<const RectangularMaskedMesh3D>& src_mesh, const DataVector<const SrcT>& src_vec,
800 const shared_ptr<const MeshD<3>>& dst_mesh, const InterpolationFlags& flags) {
801 if (src_mesh->empty()) throw BadMesh("interpolate", "source mesh empty");
802 return new NearestNeighborInterpolatedLazyDataImpl<DstT, RectangularMaskedMesh3D, SrcT>(src_mesh, src_vec, dst_mesh, flags);
803 }
804};
805
806
807 template <typename SrcT, typename DstT>
808 struct InterpolationAlgorithm<RectangularMaskedMesh3D::ElementMesh, SrcT, DstT, INTERPOLATION_LINEAR> {
809 static LazyData<DstT> interpolate(const shared_ptr<const RectangularMaskedMesh3D::ElementMesh>& src_mesh, const DataVector<const SrcT>& src_vec,
810 const shared_ptr<const MeshD<3>>& dst_mesh, const InterpolationFlags& flags) {
811 if (src_mesh->empty()) throw BadMesh("interpolate", "source mesh empty");
812 return new LinearInterpolatedLazyDataImpl<DstT, RectangularMaskedMesh3D::ElementMesh, SrcT>(src_mesh, src_vec, dst_mesh, flags);
813 }
814};
815
816 template <typename SrcT, typename DstT>
817 struct InterpolationAlgorithm<RectangularMaskedMesh3D::ElementMesh, SrcT, DstT, INTERPOLATION_NEAREST> {
818 static LazyData<DstT> interpolate(const shared_ptr<const RectangularMaskedMesh3D::ElementMesh>& src_mesh, const DataVector<const SrcT>& src_vec,
819 const shared_ptr<const MeshD<3>>& dst_mesh, const InterpolationFlags& flags) {
820 if (src_mesh->empty()) throw BadMesh("interpolate", "source mesh empty");
821 return new NearestNeighborInterpolatedLazyDataImpl<DstT, RectangularMaskedMesh3D::ElementMesh, SrcT>(src_mesh, src_vec, dst_mesh, flags);
822 }
823};
824
825} // namespace plask
826
827#endif // PLASK__RECTANGULAR_MASKED3D_H