PLaSK library
Loading...
Searching...
No Matches
rectilinear3d.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__RECTILINEAR3D_H
15#define PLASK__RECTILINEAR3D_H
16
21#include <type_traits>
22
23#include "rectangular_common.hpp"
24#include "axis1d.hpp"
25#include "../utils/interpolation.hpp"
26
27
28namespace plask {
29
39class PLASK_API RectilinearMesh3D: public RectangularMeshBase3D /*MeshD<3>*/ {
40
41 typedef std::size_t index_ft(const RectilinearMesh3D* mesh, std::size_t c0_index, std::size_t c1_index, std::size_t c2_index);
42 typedef std::size_t index012_ft(const RectilinearMesh3D* mesh, std::size_t mesh_index);
43
44 // our own virtual table, changeable at run-time:
45 index_ft* index_f;
46 index012_ft* index0_f;
47 index012_ft* index1_f;
48 index012_ft* index2_f;
49 const shared_ptr<MeshAxis>* minor_axis;
50 const shared_ptr<MeshAxis>* medium_axis;
51 const shared_ptr<MeshAxis>* major_axis;
52
53 void onAxisChanged(Event& e);
54
55 void setChangeSignal(const shared_ptr<MeshAxis>& axis) { if (axis) axis->changedConnectMethod(this, &RectilinearMesh3D::onAxisChanged); }
56 void unsetChangeSignal(const shared_ptr<MeshAxis>& axis) { if (axis) axis->changedDisconnectMethod(this, &RectilinearMesh3D::onAxisChanged); }
57
58 public:
59
63 class PLASK_API Element {
64 const RectilinearMesh3D& mesh;
65 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
66
67 public:
68
74 Element(const RectilinearMesh3D& mesh, std::size_t index0, std::size_t index1, std::size_t index2): mesh(mesh), index0(index0), index1(index1), index2(index2) {}
75
81 Element(const RectilinearMesh3D& mesh, std::size_t elementIndex): mesh(mesh) {
82 std::size_t v = mesh.getElementMeshLowIndex(elementIndex);
83 index0 = mesh.index0(v);
84 index1 = mesh.index1(v);
85 index2 = mesh.index2(v);
86 }
87
89 inline std::size_t getIndex0() const { return index0; }
90
92 inline std::size_t getIndex1() const { return index1; }
93
95 inline std::size_t getIndex2() const { return index2; }
96
98 inline std::size_t getLowerIndex0() const { return index0; }
99
101 inline std::size_t getLowerIndex1() const { return index1; }
102
104 inline std::size_t getLowerIndex2() const { return index2; }
105
107 inline double getLower0() const { return mesh.axis[0]->at(index0); }
108
110 inline double getLower1() const { return mesh.axis[1]->at(index1); }
111
113 inline double getLower2() const { return mesh.axis[2]->at(index2); }
114
116 inline std::size_t getUpperIndex0() const { return index0+1; }
117
119 inline std::size_t getUpperIndex1() const { return index1+1; }
120
122 inline std::size_t getUpperIndex2() const { return index2+1; }
123
125 inline double getUpper0() const { return mesh.axis[0]->at(getUpperIndex0()); }
126
128 inline double getUpper1() const { return mesh.axis[1]->at(getUpperIndex1()); }
129
131 inline double getUpper2() const { return mesh.axis[2]->at(getUpperIndex2()); }
132
134 inline double getSize0() const { return getUpper0() - getLower0(); }
135
137 inline double getSize1() const { return getUpper1() - getLower1(); }
138
140 inline double getSize2() const { return getUpper2() - getLower2(); }
141
143 inline Vec<3, double> getSize() const { return getUpUpUp() - getLoLoLo(); }
144
146 inline Vec<3, double> getMidpoint() const { return mesh.getElementMidpoint(index0, index1, index2); }
147
149 inline std::size_t getIndex() const { return mesh.getElementIndexFromLowIndex(getLoLoLoIndex()); }
150
152 inline std::size_t getLoLoLoIndex() const { return mesh.index(getLowerIndex0(), getLowerIndex1(), getLowerIndex2()); }
153
155 inline std::size_t getUpLoLoIndex() const { return mesh.index(getUpperIndex0(), getLowerIndex1(), getLowerIndex2()); }
156
158 inline std::size_t getLoUpLoIndex() const { return mesh.index(getLowerIndex0(), getUpperIndex1(), getLowerIndex2()); }
159
161 inline std::size_t getUpUpLoIndex() const { return mesh.index(getUpperIndex0(), getUpperIndex1(), getLowerIndex2()); }
162
164 inline std::size_t getLoLoUpIndex() const { return mesh.index(getLowerIndex0(), getLowerIndex1(), getUpperIndex2()); }
165
167 inline std::size_t getUpLoUpIndex() const { return mesh.index(getUpperIndex0(), getLowerIndex1(), getUpperIndex2()); }
168
170 inline std::size_t getLoUpUpIndex() const { return mesh.index(getLowerIndex0(), getUpperIndex1(), getUpperIndex2()); }
171
173 inline std::size_t getUpUpUpIndex() const { return mesh.index(getUpperIndex0(), getUpperIndex1(), getUpperIndex2()); }
174
176 inline Vec<3, double> getLoLoLo() const { return mesh(getLowerIndex0(), getLowerIndex1(), getLowerIndex2()); }
177
179 inline Vec<3, double> getUpLoLo() const { return mesh(getUpperIndex0(), getLowerIndex1(), getLowerIndex2()); }
180
182 inline Vec<3, double> getLoUpLo() const { return mesh(getLowerIndex0(), getUpperIndex1(), getLowerIndex2()); }
183
185 inline Vec<3, double> getUpUpLo() const { return mesh(getUpperIndex0(), getUpperIndex1(), getLowerIndex2()); }
186
188 inline Vec<3, double> getLoLoUp() const { return mesh(getLowerIndex0(), getLowerIndex1(), getUpperIndex2()); }
189
191 inline Vec<3, double> getUpLoUp() const { return mesh(getUpperIndex0(), getLowerIndex1(), getUpperIndex2()); }
192
194 inline Vec<3, double> getLoUpUp() const { return mesh(getLowerIndex0(), getUpperIndex1(), getUpperIndex2()); }
195
197 inline Vec<3, double> getUpUpUp() const { return mesh(getUpperIndex0(), getUpperIndex1(), getUpperIndex2()); }
198
200 inline Box3D toBox() const { return mesh.getElementBox(index0, index1, index2); }
201
207 bool contains(Vec<3, double> p) const { return toBox().contains(p); }
208
209 };
210
216 struct PLASK_API Elements {
217
218 static inline Element deref(const RectilinearMesh3D& mesh, std::size_t index) { return mesh.getElement(index); }
219 public:
220 typedef IndexedIterator<const RectilinearMesh3D, Element, deref> const_iterator;
221 typedef const_iterator iterator;
222
224 const RectilinearMesh3D* mesh;
225
230 Elements(const RectilinearMesh3D* mesh): mesh(mesh) {}
231
237 Element operator[](std::size_t i) const { return Element(*mesh, i); }
238
244 Element operator()(std::size_t i0, std::size_t i1, std::size_t i2) const { return Element(*mesh, i0, i1, i2); }
245
250 std::size_t size() const { return mesh->getElementsCount(); }
251
253 const_iterator begin() const { return const_iterator(mesh, 0); }
254
256 const_iterator end() const { return const_iterator(mesh, size()); }
257
258 };
259
265 template <typename BaseMeshT>
266 struct PLASK_API ElementMesh: BaseMeshT {
267
269 const BaseMeshT* originalMesh;
270
271 template <typename... Args>
272 ElementMesh(const BaseMeshT* originalMesh, Args... args): BaseMeshT(args...), originalMesh(originalMesh) {}
273
281 template <typename RandomAccessContainer>
282 auto interpolateNearestNeighbor(const RandomAccessContainer& data, const Vec<3>& point, const InterpolationFlags& flags) const
283 -> typename std::remove_reference<decltype(data[0])>::type {
284 auto p = flags.wrap(point);
285 prepareNearestNeighborInterpolationForAxis(*originalMesh->axis[0], flags, p.c0, 0);
286 prepareNearestNeighborInterpolationForAxis(*originalMesh->axis[1], flags, p.c1, 1);
287 prepareNearestNeighborInterpolationForAxis(*originalMesh->axis[2], flags, p.c2, 2);
288 size_t i0 = originalMesh->axis[0]->findUpIndex(p.c0),
289 i1 = originalMesh->axis[0]->findUpIndex(p.c1),
290 i2 = originalMesh->axis[2]->findUpIndex(p.c2);
291 if (i0 == originalMesh->axis[0]->size()) --i0;
292 if (i0 != 0) --i0;
293 if (i1 == originalMesh->axis[1]->size()) --i1;
294 if (i1 != 0) --i1;
295 if (i2 == originalMesh->axis[2]->size()) --i2;
296 if (i2 != 0) --i2;
297 return flags.postprocess(point, data[this->index(i0, i1, i2)]);
298 }
299
300 protected:
301
302 bool hasSameNodes(const MeshD<3> &to_compare) const override {
303 if (const ElementMesh* c = dynamic_cast<const ElementMesh*>(&to_compare))
304 if (this->originalMesh == c->originalMesh) return true;
305 return BaseMeshT::hasSameNodes(to_compare);
306 }
307 };
308
310 const shared_ptr<MeshAxis> axis[3];
311
313 Elements elements() const { return Elements(this); }
314 Elements getElements() const { return elements(); }
315
316 Element element(std::size_t i0, std::size_t i1, std::size_t i2) const { return Element(*this, i0, i1, i2); }
317 Element getElement(std::size_t i0, std::size_t i1, std::size_t i2) const { return element(i0, i1, i2); }
318
324 Element element(std::size_t i) const { return Element(*this, i); }
325
331 Element getElement(std::size_t i) const { return element(i); }
332
338 Box3D getElementBox(std::size_t index0, std::size_t index1, std::size_t index2) const {
339 return Box3D(axis[0]->at(index0), axis[1]->at(index1), axis[2]->at(index2),
340 axis[0]->at(index0+1), axis[1]->at(index1+1), axis[2]->at(index2+1));
341 }
342
349 enum IterationOrder { ORDER_012, ORDER_021, ORDER_102, ORDER_120, ORDER_201, ORDER_210 };
350
356 void setIterationOrder(IterationOrder order);
357
363 IterationOrder getIterationOrder() const;
364
369 const char* getIterationOrderAsArray() const {
370 return &"\0\1\2\0\2\1\1\0\2\1\2\0\2\0\1\2\1\0"[getIterationOrder() * 3];
371 // |0 1 2|0 1 2|0 1 2|0 1 2|0 1 2|0 1 2|
372 };
373
381 const char* getAxisToIterationOrder() const {
382 return &"\0\1\2\0\2\1\1\0\2\2\0\1\1\2\0\2\1\0"[getIterationOrder() * 3];
383 // |0 1 2|0 1 2|0 1 2|0 1 2|0 1 2|0 1 2|
384 };
385
391 bool isChangeSlower(std::size_t axis_index1, std::size_t axis_index2) const {
392 return getAxisToIterationOrder()[axis_index1] < getAxisToIterationOrder()[axis_index2];
393 }
394
398 void setOptimalIterationOrder();
399
404 explicit RectilinearMesh3D(IterationOrder iterationOrder = ORDER_012);
405
414 RectilinearMesh3D(shared_ptr<MeshAxis> mesh0, shared_ptr<MeshAxis> mesh1, shared_ptr<MeshAxis> mesh2, IterationOrder iterationOrder = ORDER_012);
415
423 void reset(shared_ptr<MeshAxis> mesh0, shared_ptr<MeshAxis> mesh1, shared_ptr<MeshAxis> mesh2, IterationOrder iterationOrder = ORDER_012);
424
430 RectilinearMesh3D(const RectilinearMesh3D& src, bool clone_axes = false);
431
436 void reset(const RectilinearMesh3D& src, bool clone_axes = false);
437
438 RectilinearMesh3D& operator=(const RectilinearMesh3D& src) { reset(src, true); return *this; }
439
440 RectilinearMesh3D& operator=(RectilinearMesh3D&& src) { reset(src, false); return *this; }
441
442 ~RectilinearMesh3D();
443
450 void setAxis(std::size_t axis_nr, shared_ptr<MeshAxis> new_val, bool fireResized = true);
451
452 const shared_ptr<MeshAxis> getAxis0() const { return axis[0]; }
453
454 void setAxis0(shared_ptr<MeshAxis> a0) { setAxis(0, a0); }
455
456 const shared_ptr<MeshAxis> getAxis1() const { return axis[1]; }
457
458 void setAxis1(shared_ptr<MeshAxis> a1) { setAxis(1, a1); }
459
460 const shared_ptr<MeshAxis> getAxis2() const { return axis[2]; }
461
462 void setAxis2(shared_ptr<MeshAxis> a2) { setAxis(2, a2); }
463
468 const shared_ptr<MeshAxis>& getAxis(size_t n) const {
469 if (n >= 3) throw Exception("bad axis number");
470 return axis[n];
471 }
472
474 inline const shared_ptr<MeshAxis> majorAxis() const {
475 return *major_axis;
476 }
477
479 inline std::size_t majorAxisIndex() const {
480 if (major_axis == &axis[0]) return 0;
481 if (major_axis == &axis[1]) return 1;
482 return 2;
483 }
484
486 inline const shared_ptr<MeshAxis> mediumAxis() const {
487 return *medium_axis;
488 }
489
491 inline std::size_t mediumAxisIndex() const {
492 if (medium_axis == &axis[0]) return 0;
493 if (medium_axis == &axis[1]) return 1;
494 return 2;
495 }
496
498 inline const shared_ptr<MeshAxis> minorAxis() const {
499 return *minor_axis;
500 }
501
503 inline std::size_t minorAxisIndex() const {
504 if (minor_axis == &axis[0]) return 0;
505 if (minor_axis == &axis[1]) return 1;
506 return 2;
507 }
508
514 bool operator==(const RectilinearMesh3D& to_compare) const {
515 if (empty()) return to_compare.empty();
516 // here we know that axes are non-empty
517 if (*axis[0] != *to_compare.axis[0] || *axis[1] != *to_compare.axis[1] || *axis[2] != *to_compare.axis[2]) return false;
518 // here we know that axes are equal
519 return this->getIterationOrder() == to_compare.getIterationOrder() ||
520 (unsigned(axis[0]->size() == 1) + unsigned(axis[1]->size() == 1) + unsigned(axis[2]->size() == 1)) >= 2;
521 }
522
523 bool operator!=(const RectilinearMesh3D& to_compare) const {
524 return !(*this == to_compare);
525 }
526
531 std::size_t size() const override { return axis[0]->size() * axis[1]->size() * axis[2]->size(); }
532
534 bool empty() const override { return axis[0]->empty() || axis[1]->empty() || axis[2]->empty(); }
535
543 virtual Vec<3, double> at(std::size_t index0, std::size_t index1, std::size_t index2) const = 0;
544
550 Vec<3, double> at(std::size_t index) const override {
551 return at(index0(index), index1(index), index2(index));
552 }
553
560 inline Vec<3, double> operator[](std::size_t index) const {
561 return at(index0(index), index1(index), index2(index));
562 }
563
571 inline Vec<3, double> operator()(std::size_t index0, std::size_t index1, std::size_t index2) const {
572 return at(index0, index1, index2);
573 }
574
580 double getElementMidpoint0(std::size_t index0) const { return 0.5 * (axis[0]->at(index0) + axis[0]->at(index0+1)); }
581
587 double getElementMidpoint1(std::size_t index1) const { return 0.5 * (axis[1]->at(index1) + axis[1]->at(index1+1)); }
588
594 double getElementMidpoint2(std::size_t index2) const { return 0.5 * (axis[2]->at(index2) + axis[2]->at(index2+1)); }
595
601 virtual Vec<3, double> getElementMidpoint(std::size_t index0, std::size_t index1, std::size_t index2) const = 0;
602
610 std::size_t index(std::size_t c0_index, std::size_t c1_index, std::size_t c2_index) const {
611 return index_f(this, c0_index, c1_index, c2_index);
612 }
613
619 inline std::size_t index(const Vec<3, std::size_t>& indexes) const {
620 return index(indexes[0], indexes[1], indexes[2]);
621 }
622
628 inline std::size_t index0(std::size_t mesh_index) const {
629 return index0_f(this, mesh_index);
630 }
631
637 inline std::size_t index1(std::size_t mesh_index) const {
638 return index1_f(this, mesh_index);
639 }
640
646 inline std::size_t index2(std::size_t mesh_index) const {
647 return index2_f(this, mesh_index);
648 }
649
655 template <std::size_t axis_nr>
656 inline std::size_t index_axis(std::size_t mesh_index) const {
657 if (axis_nr == 0) return index0(mesh_index);
658 if (axis_nr == 1) return index1(mesh_index);
659 if (axis_nr == 2) return index2(mesh_index);
660 assert(false);
661 }
662
668 inline Vec<3, std::size_t> indexes(std::size_t mesh_index) const {
669 return Vec<3, std::size_t>(index0(mesh_index), index1(mesh_index), index2(mesh_index));
670 }
671
677 inline std::size_t majorIndex(std::size_t mesh_index) const {
678 return mesh_index / (*minor_axis)->size() / (*medium_axis)->size();
679 }
680
686 inline std::size_t middleIndex(std::size_t mesh_index) const {
687 return (mesh_index / (*minor_axis)->size()) % (*medium_axis)->size();
688 }
689
695 inline std::size_t minorIndex(std::size_t mesh_index) const {
696 return mesh_index % (*minor_axis)->size();
697 }
698
703 std::size_t getElementsCount0() const {
704 const std::size_t s = axis[0]->size();
705 return (s != 0)? s-1 : 0;
706 }
707
712 std::size_t getElementsCount1() const {
713 const std::size_t s = axis[1]->size();
714 return (s != 0)? s-1 : 0;
715 }
716
721 size_t getElementsCount2() const {
722 const std::size_t s = axis[2]->size();
723 return (s != 0)? s-1 : 0;
724 }
725
730 size_t getElementsCount() const {
731 return size_t(std::max(int(axis[0]->size())-1, 0) * std::max(int(axis[1]->size())-1, 0) * std::max(int(axis[2]->size())-1, 0));
732 }
733
739 std::size_t getElementMeshLowIndex(std::size_t element_index) const {
740 const std::size_t minor_size_minus_1 = (*minor_axis)->size()-1;
741 const std::size_t elements_per_level = minor_size_minus_1 * ((*medium_axis)->size()-1);
742 return element_index + (element_index / elements_per_level) * ((*medium_axis)->size() + minor_size_minus_1)
743 + (element_index % elements_per_level) / minor_size_minus_1;
744 }
745
751 std::size_t getElementIndexFromLowIndexes(std::size_t axis0_index, std::size_t axis1_index, std::size_t axis2_index) const {
752 return getElementIndexFromLowIndex(index(axis0_index, axis1_index, axis2_index));
753 }
754
760 bool isLowIndexOfElement(std::size_t meshIndex) const {
761 return index0(meshIndex) + 1 < axis[0]->size() && index1(meshIndex) + 1 < axis[1]->size() && index2(meshIndex) + 1 < axis[2]->size();
762 }
763
769 std::size_t getElementIndexFromLowIndex(std::size_t mesh_index_of_el_bottom_left) const {
770 const std::size_t verticles_per_level = (*minor_axis)->size() * (*medium_axis)->size();
771 return mesh_index_of_el_bottom_left - (mesh_index_of_el_bottom_left / verticles_per_level) * ((*medium_axis)->size() + (*minor_axis)->size() - 1)
772 - (mesh_index_of_el_bottom_left % verticles_per_level) / (*minor_axis)->size();
773 }
774
780 std::size_t getElementIndexFromLowIndex(std::size_t axis0_index, std::size_t axis1_index, std::size_t axis2_index) const {
781 return getElementIndexFromLowIndex(index(axis0_index, axis1_index, axis2_index));
782 }
783
790 Vec<3, std::size_t> getElementMeshLowIndexes(std::size_t element_index) const {
791 std::size_t bl_index = getElementMeshLowIndex(element_index);
792 return Vec<3, std::size_t>(index0(bl_index), index1(bl_index), index2(bl_index));
793 }
794
801 template <typename RandomAccessContainer>
802 auto interpolateLinear(const RandomAccessContainer& data, const Vec<3>& point, const InterpolationFlags& flags) const
803 -> typename std::remove_reference<decltype(data[0])>::type
804 {
805 auto p = flags.wrap(point);
806
807 size_t index0_lo, index0_hi;
808 double back, front;
809 bool invert_back, invert_front;
810 prepareInterpolationForAxis(*axis[0], flags, p.c0, 0, index0_lo, index0_hi, back, front, invert_back, invert_front);
811
812 size_t index1_lo, index1_hi;
813 double left, right;
814 bool invert_left, invert_right;
815 prepareInterpolationForAxis(*axis[1], flags, p.c1, 1, index1_lo, index1_hi, left, right, invert_left, invert_right);
816
817 size_t index2_lo, index2_hi;
818 double bottom, top;
819 bool invert_bottom, invert_top;
820 prepareInterpolationForAxis(*axis[2], flags, p.c2, 2, index2_lo, index2_hi, bottom, top, invert_bottom, invert_top);
821
822 // all indexes are in bounds
823 typename std::remove_const<typename std::remove_reference<decltype(data[0])>::type>::type
824 data_lll = data[index(index0_lo, index1_lo, index2_lo)],
825 data_hll = data[index(index0_hi, index1_lo, index2_lo)],
826 data_hhl = data[index(index0_hi, index1_hi, index2_lo)],
827 data_lhl = data[index(index0_lo, index1_hi, index2_lo)],
828 data_llh = data[index(index0_lo, index1_lo, index2_hi)],
829 data_hlh = data[index(index0_hi, index1_lo, index2_hi)],
830 data_hhh = data[index(index0_hi, index1_hi, index2_hi)],
831 data_lhh = data[index(index0_lo, index1_hi, index2_hi)];
832
833 if (invert_back) { data_lll = flags.reflect(0, data_lll); data_llh = flags.reflect(0, data_llh); data_lhl = flags.reflect(0, data_lhl); data_lhh = flags.reflect(0, data_lhh); }
834 if (invert_front) { data_hll = flags.reflect(0, data_hll); data_llh = flags.reflect(0, data_hlh); data_lhl = flags.reflect(0, data_hhl); data_lhh = flags.reflect(0, data_hhh); }
835 if (invert_left) { data_lll = flags.reflect(1, data_lll); data_llh = flags.reflect(1, data_llh); data_hll = flags.reflect(1, data_hll); data_hlh = flags.reflect(1, data_hlh); }
836 if (invert_right) { data_lhl = flags.reflect(1, data_lhl); data_llh = flags.reflect(1, data_lhh); data_hll = flags.reflect(1, data_hhl); data_hlh = flags.reflect(1, data_hhh); }
837 if (invert_bottom) { data_lll = flags.reflect(2, data_lll); data_lhl = flags.reflect(2, data_lhl); data_hll = flags.reflect(2, data_hll); data_hhl = flags.reflect(2, data_hhl); }
838 if (invert_top) { data_llh = flags.reflect(2, data_llh); data_lhl = flags.reflect(2, data_lhh); data_hll = flags.reflect(2, data_hlh); data_hhl = flags.reflect(2, data_hhh); }
839
840 return flags.postprocess(point,
841 interpolation::trilinear(back, front, left, right, bottom, top,
842 data_lll, data_hll, data_hhl, data_lhl, data_llh, data_hlh, data_hhh, data_lhh,
843 p.c0, p.c1, p.c2));
844 }
845
852 template <typename RandomAccessContainer>
853 auto interpolateNearestNeighbor(const RandomAccessContainer& data, Vec<3> point, const InterpolationFlags& flags) const
854 -> typename std::remove_reference<decltype(data[0])>::type {
855 auto p = flags.wrap(point);
856 prepareNearestNeighborInterpolationForAxis(*axis[0], flags, p.c0, 0);
857 prepareNearestNeighborInterpolationForAxis(*axis[1], flags, p.c1, 1);
858 prepareNearestNeighborInterpolationForAxis(*axis[2], flags, p.c2, 2);
859 return flags.postprocess(point, data[this->index(axis[0]->findNearestIndex(p.c0), axis[1]->findNearestIndex(p.c1), axis[2]->findNearestIndex(p.c2))]);
860 }
861
862protected:
863 bool hasSameNodes(const MeshD<3> &to_compare) const override {
864 if (const RectilinearMesh3D* c = dynamic_cast<const RectilinearMesh3D*>(&to_compare))
865 return *this == *c; // this will call == operator from RectilinearMesh3D
866 return RectangularMeshBase3D::hasSameNodes(to_compare);
867 }
868
869private:
870
878 template <int CHANGE_DIR_SLOWER, int CHANGE_DIR_FASTER>
879 struct BoundaryIteratorImpl: public plask::BoundaryNodeSetImpl::IteratorImpl {
880
881 const RectilinearMesh3D &mesh;
882
884 Vec<3, std::size_t> index;
885
887 const std::size_t indexFasterBegin, indexFasterEnd, indexSlowerEnd;
888
889 public:
890 BoundaryIteratorImpl(const RectilinearMesh3D& mesh, Vec<3, std::size_t> index, std::size_t indexSlowerEnd, std::size_t indexFasterEnd)
891 : mesh(mesh), index(index), indexFasterBegin(index[CHANGE_DIR_FASTER]), indexFasterEnd(indexFasterEnd), indexSlowerEnd(indexSlowerEnd)
892 {
893 }
894
895 void increment() override {
896 ++index[CHANGE_DIR_FASTER];
897 if (index[CHANGE_DIR_FASTER] == indexFasterEnd) {
898 index[CHANGE_DIR_FASTER] = indexFasterBegin;
899 ++index[CHANGE_DIR_SLOWER];
900 }
901 }
902
903 bool equal(const plask::BoundaryNodeSetImpl::IteratorImpl& other) const override {
904 return index == static_cast<const BoundaryIteratorImpl&>(other).index;
905 }
906
907 std::size_t dereference() const override {
908 return mesh.index(index);
909 }
910
911 std::unique_ptr<plask::BoundaryNodeSetImpl::IteratorImpl> clone() const override {
912 return std::unique_ptr<plask::BoundaryNodeSetImpl::IteratorImpl>(new BoundaryIteratorImpl<CHANGE_DIR_SLOWER, CHANGE_DIR_FASTER>(*this));
913 }
914
915 };
916
917 template <int CHANGE_DIR_SLOWER, int CHANGE_DIR_FASTER>
918 struct BoundaryNodeSetImpl: public BoundaryNodeSetWithMeshImpl<RectilinearMesh3D> {
919
920 static constexpr int FIXED_DIR = 3 - CHANGE_DIR_SLOWER - CHANGE_DIR_FASTER;
921
922 using typename BoundaryNodeSetWithMeshImpl<RectilinearMesh3D>::const_iterator;
923
925 std::size_t level;
926
928 std::size_t indexFasterEnd, indexSlowerEnd;
929
930 BoundaryNodeSetImpl(const RectilinearMesh3D& mesh, std::size_t level)
931 : BoundaryNodeSetWithMeshImpl<RectilinearMesh3D>(mesh), level(level) {}
932
933 bool contains(std::size_t mesh_index) const override {
934 return mesh_index < this->mesh.size() && this->mesh.template index_axis<FIXED_DIR>(mesh_index) == level;
935 }
936
937 const_iterator begin() const override {
938 Vec<3, std::size_t> index_begin(0, 0, 0);
939 index_begin[FIXED_DIR] = level;
940 return Iterator(new BoundaryIteratorImpl<CHANGE_DIR_SLOWER, CHANGE_DIR_FASTER>(this->mesh, index_begin,
941 this->mesh.axis[CHANGE_DIR_SLOWER]->size(),
942 this->mesh.axis[CHANGE_DIR_FASTER]->size()));
943 }
944
945 const_iterator end() const override {
946 Vec<3, std::size_t> index_end;
947 index_end[CHANGE_DIR_SLOWER] = this->mesh.axis[CHANGE_DIR_SLOWER]->size();
948 index_end[FIXED_DIR] = level;
949 index_end[CHANGE_DIR_FASTER] = 0;
950 return Iterator(new BoundaryIteratorImpl<CHANGE_DIR_SLOWER, CHANGE_DIR_FASTER>(this->mesh, index_end,
951 this->mesh.axis[CHANGE_DIR_SLOWER]->size(),
952 this->mesh.axis[CHANGE_DIR_FASTER]->size()));
953 }
954
955 std::size_t size() const override {
956 return this->mesh.axis[CHANGE_DIR_FASTER]->size() * this->mesh.axis[CHANGE_DIR_SLOWER]->size();
957 }
958
959 bool empty() const override {
960 return this->mesh.axis[CHANGE_DIR_FASTER]->empty() || this->mesh.axis[CHANGE_DIR_SLOWER]->empty();
961 }
962 };
963
964 template <int CHANGE_DIR_SLOWER, int CHANGE_DIR_FASTER>
965 struct BoundaryNodeSetRangeImpl: public BoundaryNodeSetWithMeshImpl<RectilinearMesh3D> {
966
967 static constexpr int FIXED_DIR = 3 - CHANGE_DIR_SLOWER - CHANGE_DIR_FASTER;
968
969 typedef typename BoundaryNodeSetImpl::Iterator Iterator;
970
972 Vec<3, std::size_t> index;
973
975 std::size_t indexFasterEnd, indexSlowerEnd;
976
977 BoundaryNodeSetRangeImpl(const RectilinearMesh3D& mesh, Vec<3, std::size_t> index, std::size_t indexSlowerEnd, std::size_t indexFasterEnd)
978 : BoundaryNodeSetWithMeshImpl<RectilinearMesh3D>(mesh), index(index), indexFasterEnd(indexFasterEnd), indexSlowerEnd(indexSlowerEnd) {}
979
980 BoundaryNodeSetRangeImpl(const RectilinearMesh3D& mesh, std::size_t index0, std::size_t index1, std::size_t index2, std::size_t indexSlowerEnd, std::size_t indexFasterEnd)
981 : BoundaryNodeSetWithMeshImpl<RectilinearMesh3D>(mesh), index(index0, index1, index2), indexFasterEnd(indexFasterEnd), indexSlowerEnd(indexSlowerEnd) {}
982
983 bool contains(std::size_t mesh_index) const override {
984 if (mesh_index >= this->mesh.size()) return false;
985 Vec<3, std::size_t> mesh_indexes = this->mesh.indexes(mesh_index);
986 return mesh_indexes[FIXED_DIR] == index[FIXED_DIR] &&
987 (index[CHANGE_DIR_FASTER] <= mesh_indexes[CHANGE_DIR_FASTER] && mesh_indexes[CHANGE_DIR_FASTER] < indexFasterEnd) &&
988 (index[CHANGE_DIR_SLOWER] <= mesh_indexes[CHANGE_DIR_SLOWER] && mesh_indexes[CHANGE_DIR_SLOWER] < indexSlowerEnd);
989 }
990
991 const_iterator begin() const override {
992 return Iterator(new BoundaryIteratorImpl<CHANGE_DIR_SLOWER, CHANGE_DIR_FASTER>(this->mesh, index, indexSlowerEnd, indexFasterEnd));
993 }
994
995 const_iterator end() const override {
996 Vec<3, std::size_t> index_end = index;
997 index_end[CHANGE_DIR_SLOWER] = indexSlowerEnd;
998 return Iterator(new BoundaryIteratorImpl<CHANGE_DIR_SLOWER, CHANGE_DIR_FASTER>(this->mesh, index_end, indexSlowerEnd, indexFasterEnd));
999 }
1000
1001 std::size_t size() const override {
1002 return (indexFasterEnd - index[CHANGE_DIR_FASTER]) * (indexSlowerEnd - index[CHANGE_DIR_SLOWER]);
1003 }
1004
1005 bool empty() const override {
1006 return indexFasterEnd == index[CHANGE_DIR_FASTER] || indexSlowerEnd == index[CHANGE_DIR_SLOWER];
1007 }
1008 };
1009
1010 public:
1011
1012 BoundaryNodeSet createIndex0BoundaryAtLine(std::size_t line_nr_axis0) const override;
1013
1014 BoundaryNodeSet createBackBoundary() const override;
1015
1016 BoundaryNodeSet createFrontBoundary() const override;
1017
1018 BoundaryNodeSet createIndex1BoundaryAtLine(std::size_t line_nr_axis1) const override;
1019
1020 BoundaryNodeSet createLeftBoundary() const override;
1021
1022 BoundaryNodeSet createRightBoundary() const override;
1023
1024 BoundaryNodeSet createIndex2BoundaryAtLine(std::size_t line_nr_axis2) const override;
1025
1026 BoundaryNodeSet createBottomBoundary() const override;
1027
1028 BoundaryNodeSet createTopBoundary() const override;
1029
1030 BoundaryNodeSet createIndex0BoundaryAtLine(std::size_t line_nr_axis0,
1031 std::size_t index1Begin, std::size_t index1End,
1032 std::size_t index2Begin, std::size_t index2End) const override;
1033
1034 BoundaryNodeSet createBackOfBoundary(const Box3D& box) const override;
1035
1036 BoundaryNodeSet createFrontOfBoundary(const Box3D& box) const override;
1037
1038 BoundaryNodeSet createIndex1BoundaryAtLine(std::size_t line_nr_axis1,
1039 std::size_t index0Begin, std::size_t index0End,
1040 std::size_t index2Begin, std::size_t index2End) const override;
1041
1042 BoundaryNodeSet createLeftOfBoundary(const Box3D& box) const override;
1043
1044 BoundaryNodeSet createRightOfBoundary(const Box3D& box) const override;
1045
1046 BoundaryNodeSet createIndex2BoundaryAtLine(std::size_t line_nr_axis2,
1047 std::size_t index0Begin, std::size_t index0End,
1048 std::size_t index1Begin, std::size_t index1End) const override;
1049
1050 BoundaryNodeSet createBottomOfBoundary(const Box3D& box) const override;
1051
1052 BoundaryNodeSet createTopOfBoundary(const Box3D& box) const override;
1053};
1054
1055
1056} // namespace plask
1057
1058#endif // PLASK__RECTILINEAR3D_H