PLaSK library

Loading...
Searching...
No Matches
rectangular2d.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__RECTANGULAR2D_H
15#define PLASK__RECTANGULAR2D_H
16
21#include <iterator>
22
23#include "rectangular_common.hpp"
24#include "../utils/interpolation.hpp"
25#include "../geometry/object.hpp"
26#include "../geometry/space.hpp"
27#include "../math.hpp"
28
29 namespace plask {
30
39 class PLASK_API RectangularMesh2D: public RectangularMeshBase2D {
40
41 typedef std::size_t index_ft(const RectangularMesh2D* mesh, std::size_t axis0_index, std::size_t axis1_index);
42 typedef std::size_t index01_ft(const RectangularMesh2D* mesh, std::size_t mesh_index);
43
44 // Our own virtual table, changeable at run-time:
45 index_ft* index_f;
46 index01_ft* index0_f;
47 index01_ft* index1_f;
48
49 const shared_ptr<MeshAxis>* minor_axis;
50 const shared_ptr<MeshAxis>* major_axis;
51
52 void onAxisChanged(Event& e);
53
54 void setChangeSignal(const shared_ptr<MeshAxis>& axis) { if (axis) axis->changedConnectMethod(this, &RectangularMesh2D::onAxisChanged); }
55 void unsetChangeSignal(const shared_ptr<MeshAxis>& axis) { if (axis) axis->changedDisconnectMethod(this, &RectangularMesh2D::onAxisChanged); }
56
57 public:
58
62 class PLASK_API Element {
63 const RectangularMesh2D& mesh;
64 std::size_t index0, index1; // probably this form allows us to do most operation the fastest in average, low indexes of element corner or just element indexes
65
66 public:
67
73 Element(const RectangularMesh2D& mesh, std::size_t index0, std::size_t index1): mesh(mesh), index0(index0), index1(index1) {}
74
80 Element(const RectangularMesh2D& mesh, std::size_t elementIndex): mesh(mesh) {
81 std::size_t v = mesh.getElementMeshLowIndex(elementIndex);
82 index0 = mesh.index0(v);
83 index1 = mesh.index1(v);
84 }
85
87 inline std::size_t getIndex0() const { return index0; }
88
90 inline std::size_t getIndex1() const { return index1; }
91
93 inline std::size_t getLowerIndex0() const { return index0; }
94
96 inline std::size_t getLowerIndex1() const { return index1; }
97
99 inline double getLower0() const { return mesh.axis[0]->at(index0); }
100
102 inline double getLower1() const { return mesh.axis[1]->at(index1); }
103
105 inline std::size_t getUpperIndex0() const { return index0+1; }
106
108 inline std::size_t getUpperIndex1() const { return index1+1; }
109
111 inline double getUpper0() const { return mesh.axis[0]->at(getUpperIndex0()); }
112
114 inline double getUpper1() const { return mesh.axis[1]->at(getUpperIndex1()); }
115
117 inline double getSize0() const { return getUpper0() - getLower0(); }
118
120 inline double getSize1() const { return getUpper1() - getLower1(); }
121
123 inline Vec<2, double> getSize() const { return getUpUp() - getLoLo(); }
124
126 inline Vec<2, double> getMidpoint() const { return mesh.getElementMidpoint(index0, index1); }
127
129 inline std::size_t getIndex() const { return mesh.getElementIndexFromLowIndexes(getLowerIndex0(), getLowerIndex1()); }
130
132 inline Box2D toBox() const { return mesh.getElementBox(index0, index1); }
133
139 bool contains(Vec<2, double> p) const { return toBox().contains(p); }
140
142 inline double getVolume() const { return getSize0() * getSize1(); }
143
145 inline double getArea() const { return getVolume(); }
146
148 inline std::size_t getLoLoIndex() const { return mesh.index(getLowerIndex0(), getLowerIndex1()); }
149
151 inline std::size_t getLoUpIndex() const { return mesh.index(getLowerIndex0(), getUpperIndex1()); }
152
154 inline std::size_t getUpLoIndex() const { return mesh.index(getUpperIndex0(), getLowerIndex1()); }
155
157 inline std::size_t getUpUpIndex() const { return mesh.index(getUpperIndex0(), getUpperIndex1()); }
158
160 inline Vec<2, double> getLoLo() const { return mesh(getLowerIndex0(), getLowerIndex1()); }
161
163 inline Vec<2, double> getLoUp() const { return mesh(getLowerIndex0(), getUpperIndex1()); }
164
166 inline Vec<2, double> getUpLo() const { return mesh(getUpperIndex0(), getLowerIndex1()); }
167
169 inline Vec<2, double> getUpUp() const { return mesh(getUpperIndex0(), getUpperIndex1()); }
170
171 };
172
178 class PLASK_API Elements {
179
180 static inline Element deref(const RectangularMesh2D& mesh, std::size_t index) { return mesh.getElement(index); }
181 public:
182 typedef IndexedIterator<const RectangularMesh2D, Element, deref> const_iterator;
183 typedef const_iterator iterator;
184
185 const RectangularMesh2D* mesh;
186
187 Elements(const RectangularMesh2D* mesh): mesh(mesh) {}
188
194 Element operator[](std::size_t i) const { return Element(*mesh, i); }
195
201 Element operator()(std::size_t i0, std::size_t i1) const { return Element(*mesh, i0, i1); }
202
207 std::size_t size() const { return mesh->getElementsCount(); }
208
210 const_iterator begin() const { return const_iterator(mesh, 0); }
211
213 const_iterator end() const { return const_iterator(mesh, size()); }
214
215 };
216
222 class ElementMesh;
223
225 const shared_ptr<MeshAxis> axis[2];
226
228 Elements elements() const { return Elements(this); }
229 Elements getElements() const { return elements(); }
230
231 Element element(std::size_t i0, std::size_t i1) const { return Element(*this, i0, i1); }
232 Element getElement(std::size_t i0, std::size_t i1) const { return element(i0, i1); }
233
239 Element element(std::size_t i) const { return Element(*this, i); }
240
246 Element getElement(std::size_t i) const { return element(i); }
247
256 enum IterationOrder { ORDER_10, ORDER_01 };
257
263 void setIterationOrder(IterationOrder order);
264
270 IterationOrder getIterationOrder() const;
271
276 const char* getIterationOrderAsArray() const {
277 return &"\1\0\0\1"[getIterationOrder() * 2];
278 // | | |
279 };
280
288 const char* getAxisToIterationOrder() const {
289 return &"\1\0\0\1"[getIterationOrder() * 2];
290 // | | |
291 };
292
296 void setOptimalIterationOrder() {
297 setIterationOrder(axis[0]->size() > axis[1]->size() ? ORDER_01 : ORDER_10);
298 }
299
304 explicit RectangularMesh2D(IterationOrder iterationOrder = ORDER_01);
305
313 RectangularMesh2D(shared_ptr<MeshAxis> axis0, shared_ptr<MeshAxis> axis1, IterationOrder iterationOrder = ORDER_01);
314
321 void reset(shared_ptr<MeshAxis> axis0, shared_ptr<MeshAxis> axis1, IterationOrder iterationOrder = ORDER_01);
322
328 RectangularMesh2D(const RectangularMesh2D& src, bool clone_axes = false);
329
335 void reset(const RectangularMesh2D& src, bool clone_axes = false);
336
337 RectangularMesh2D& operator=(const RectangularMesh2D& src) { reset(src, true); return *this; }
338
339 RectangularMesh2D& operator=(RectangularMesh2D&& src) { reset(src, false); return *this; }
340
341 ~RectangularMesh2D();
342
349 void setAxis(std::size_t axis_nr, shared_ptr<MeshAxis> new_val, bool fireResized = true);
350
351 const shared_ptr<MeshAxis> getAxis0() const { return axis[0]; }
352
353 void setAxis0(shared_ptr<MeshAxis> a0) { setAxis(0, a0); }
354
355 const shared_ptr<MeshAxis> getAxis1() const { return axis[1]; }
356
357 void setAxis1(shared_ptr<MeshAxis> a1) { setAxis(1, a1); }
358
359
360
361 /*
362 * Construct mesh with is based on given 1D meshes
363 *
364 * @param mesh0 mesh for the first coordinate, or constructor argument for the first coordinate mesh
365 * @param mesh1 mesh for the second coordinate, or constructor argument for the second coordinate mesh
366 * @param iterationOrder iteration order
367 */
368 /*template <typename Mesh0CtorArg, typename Mesh1CtorArg>
369 RectangularMesh(Mesh0CtorArg&& mesh0, Mesh1CtorArg&& mesh1, IterationOrder iterationOrder = ORDER_10):
370 axis0(std::forward<Mesh0CtorArg>(mesh0)), axis1(std::forward<Mesh1CtorArg>(mesh1)) elements(this) {
371 axis0->owner = this; axis[1]->owner = this;
372 setIterationOrder(iterationOrder); }*/
373
378 const shared_ptr<MeshAxis>& tran() const { return axis[0]; }
379
380 void setTran(shared_ptr<MeshAxis> a0) { setAxis0(a0); }
381
386 const shared_ptr<MeshAxis>& vert() const { return axis[1]; }
387
388 void setVert(shared_ptr<MeshAxis> a1) { setAxis1(a1); }
389
394 const shared_ptr<MeshAxis>& ee_x() const { return axis[0]; }
395
400 const shared_ptr<MeshAxis>& ee_y() const { return axis[1]; }
401
406 const shared_ptr<MeshAxis>& rad_r() const { return axis[0]; }
407
412 const shared_ptr<MeshAxis>& rad_z() const { return axis[1]; }
413
419 const shared_ptr<MeshAxis>& getAxis(size_t n) const {
420 if (n >= 2) throw Exception("bad axis number");
421 return axis[n];
422 }
423
425 inline const shared_ptr<MeshAxis> majorAxis() const {
426 return *major_axis;
427 }
428
430 inline std::size_t majorAxisIndex() const {
431 return getIterationOrder() == IterationOrder::ORDER_01 ? 0 : 1;
432 }
433
435 inline const shared_ptr<MeshAxis> minorAxis() const {
436 return *minor_axis;
437 }
438
440 inline std::size_t minorAxisIndex() const {
441 return getIterationOrder() == IterationOrder::ORDER_01 ? 1 : 0;
442 }
443
449 bool operator==(const RectangularMesh2D& to_compare) const {
450 if (empty()) return to_compare.empty();
451 // here we know that axes are non-empty
452 if (*axis[0] != *to_compare.axis[0] || *axis[1] != *to_compare.axis[1]) return false;
453 // here we know that axes are equal
454 return this->getIterationOrder() == to_compare.getIterationOrder() || (axis[0]->size() == 1) || (axis[1]->size() == 1);
455 }
456
457 bool operator!=(const RectangularMesh2D& to_compare) const {
458 return !(*this == to_compare);
459 }
460
465 std::size_t size() const override { return axis[0]->size() * axis[1]->size(); }
466
471 std::size_t getMaxSize() const { return std::max(axis[0]->size(), axis[1]->size()); }
472
477 std::size_t getMinSize() const { return std::min(axis[0]->size(), axis[1]->size()); }
478
483 void writeXML(XMLElement& object) const override;
484
486 bool empty() const override { return axis[0]->empty() || axis[1]->empty(); }
487
494 inline std::size_t index(std::size_t axis0_index, std::size_t axis1_index) const {
495 return index_f(this, axis0_index, axis1_index);
496 }
497
503 inline std::size_t index(const Vec<2, std::size_t>& indexes) const {
504 return index(indexes[0], indexes[1]);
505 }
506
512 inline std::size_t index0(std::size_t mesh_index) const {
513 return index0_f(this, mesh_index);
514 }
515
521 inline std::size_t index1(std::size_t mesh_index) const {
522 return index1_f(this, mesh_index);
523 }
524
530 inline Vec<2, std::size_t> indexes(std::size_t mesh_index) const {
531 return Vec<2, std::size_t>(index0(mesh_index), index1(mesh_index));
532 }
533
539 inline std::size_t majorIndex(std::size_t mesh_index) const {
540 return mesh_index / (*minor_axis)->size();
541 }
542
548 inline std::size_t minorIndex(std::size_t mesh_index) const {
549 return mesh_index % (*minor_axis)->size();
550 }
551
558 inline Vec<2, double> at(std::size_t index0, std::size_t index1) const {
559 return Vec<2, double>(axis[0]->at(index0), axis[1]->at(index1));
560 }
561
567 Vec<2, double> at(std::size_t index) const override {
568 return Vec<2, double>(axis[0]->at(index0(index)), axis[1]->at(index1(index)));
569 }
570
577 inline Vec<2,double> operator[](std::size_t index) const {
578 return Vec<2, double>(axis[0]->at(index0(index)), axis[1]->at(index1(index)));
579 }
580
587 inline Vec<2,double> operator()(std::size_t axis0_index, std::size_t axis1_index) const {
588 return Vec<2, double>(axis[0]->at(axis0_index), axis[1]->at(axis1_index));
589 }
590
594 /*void clear() {
595 axis0->clear();
596 axis[1]->clear();
597 }*/
598
603 shared_ptr<RectangularMesh2D::ElementMesh> getElementMesh() const;
604
611 template <typename RandomAccessContainer>
612 auto interpolateLinear(const RandomAccessContainer& data, const Vec<2>& point, const InterpolationFlags& flags) const
613 -> typename std::remove_reference<decltype(data[0])>::type
614 {
615 auto p = flags.wrap(point);
616
617 size_t index0_lo, index0_hi;
618 double left, right;
619 bool invert_left, invert_right;
620 prepareInterpolationForAxis(*axis[0], flags, p.c0, 0, index0_lo, index0_hi, left, right, invert_left, invert_right);
621
622 size_t index1_lo, index1_hi;
623 double bottom, top;
624 bool invert_bottom, invert_top;
625 prepareInterpolationForAxis(*axis[1], flags, p.c1, 1, index1_lo, index1_hi, bottom, top, invert_bottom, invert_top);
626
627 typename std::remove_const<typename std::remove_reference<decltype(data[0])>::type>::type
628 data_lb = data[index(index0_lo, index1_lo)],
629 data_rb = data[index(index0_hi, index1_lo)],
630 data_rt = data[index(index0_hi, index1_hi)],
631 data_lt = data[index(index0_lo, index1_hi)];
632 if (invert_left) { data_lb = flags.reflect(0, data_lb); data_lt = flags.reflect(0, data_lt); }
633 if (invert_right) { data_rb = flags.reflect(0, data_rb); data_rt = flags.reflect(0, data_rt); }
634 if (invert_top) { data_lt = flags.reflect(1, data_lt); data_rt = flags.reflect(1, data_rt); }
635 if (invert_bottom) { data_lb = flags.reflect(1, data_lb); data_rb = flags.reflect(1, data_rb); }
636
637 return flags.postprocess(point, interpolation::bilinear(left, right, bottom, top, data_lb, data_rb, data_rt, data_lt, p.c0, p.c1));
638 }
639
646 template <typename RandomAccessContainer>
647 auto interpolateNearestNeighbor(const RandomAccessContainer& data, const Vec<2>& point, const InterpolationFlags& flags) const
648 -> typename std::remove_reference<decltype(data[0])>::type {
649 auto p = flags.wrap(point);
650 prepareNearestNeighborInterpolationForAxis(*axis[0], flags, p.c0, 0);
651 prepareNearestNeighborInterpolationForAxis(*axis[1], flags, p.c1, 1);
652 return flags.postprocess(point, data[this->index(axis[0]->findNearestIndex(p.c0), axis[1]->findNearestIndex(p.c1))]);
653 }
654
659 std::size_t getElementsCount0() const {
660 const std::size_t s = axis[0]->size();
661 return (s != 0)? s-1 : 0;
662 }
663
668 std::size_t getElementsCount1() const {
669 const std::size_t s = axis[1]->size();
670 return (s != 0)? s-1 : 0;
671 }
672
677 std::size_t getElementsCount() const {
678 return (axis[0]->size() != 0 && axis[1]->size() != 0)?
679 (axis[0]->size()-1) * (axis[1]->size()-1) : 0;
680 }
681
687 bool isLowIndexOfElement(std::size_t meshIndex) const {
688 return index0(meshIndex) + 1 < axis[0]->size() && index1(meshIndex) + 1 < axis[1]->size();
689 }
690
696 std::size_t getElementIndexFromLowIndex(std::size_t mesh_index_of_el_bottom_left) const {
697 return mesh_index_of_el_bottom_left - mesh_index_of_el_bottom_left / (*minor_axis)->size();
698 }
699
706 std::size_t getElementIndexFromLowIndexes(std::size_t axis0_index, std::size_t axis1_index) const {
707 return getElementIndexFromLowIndex(index(axis0_index, axis1_index));
708 }
709
715 std::size_t getElementMeshLowIndex(std::size_t element_index) const {
716 return element_index + (element_index / ((*minor_axis)->size()-1));
717 }
718
725 Vec<2, std::size_t> getElementMeshLowIndexes(std::size_t element_index) const {
726 std::size_t bl_index = getElementMeshLowIndex(element_index);
727 return Vec<2, std::size_t>(index0(bl_index), index1(bl_index));
728 }
729
735 double getElementArea(std::size_t index0, std::size_t index1) const {
736 return (axis[0]->at(index0+1) - axis[0]->at(index0)) * (axis[1]->at(index1+1) - axis[1]->at(index1));
737 }
738
744 double getElementArea(std::size_t element_index) const {
745 std::size_t bl_index = getElementMeshLowIndex(element_index);
746 return getElementArea(index0(bl_index), index1(bl_index));
747 }
748
754 double getElementMidpoint0(std::size_t index0) const { return 0.5 * (axis[0]->at(index0) + axis[0]->at(index0+1)); }
755
761 double getElementMidpoint1(std::size_t index1) const { return 0.5 * (axis[1]->at(index1) + axis[1]->at(index1+1)); }
762
768 Vec<2, double> getElementMidpoint(std::size_t index0, std::size_t index1) const {
769 return vec(getElementMidpoint0(index0), getElementMidpoint1(index1));
770 }
771
777 Vec<2, double> getElementMidpoint(std::size_t element_index) const {
778 std::size_t bl_index = getElementMeshLowIndex(element_index);
779 return getElementMidpoint(index0(bl_index), index1(bl_index));
780 }
781
787 Box2D getElementBox(std::size_t index0, std::size_t index1) const {
788 return Box2D(axis[0]->at(index0), axis[1]->at(index1), axis[0]->at(index0+1), axis[1]->at(index1+1));
789 }
790
796 Box2D getElementBox(std::size_t element_index) const {
797 std::size_t bl_index = getElementMeshLowIndex(element_index);
798 return getElementBox(index0(bl_index), index1(bl_index));
799 }
800
801 protected:
802
803 bool hasSameNodes(const MeshD<2> &to_compare) const override;
804
805
806 private:
807
808 // Common code for: left, right, bottom, top boundaries:
809 struct BoundaryIteratorImpl: public BoundaryNodeSetImpl::IteratorImpl {
810
811 const RectangularMesh2D &mesh;
812
813 std::size_t line;
814
815 std::size_t index;
816
817 BoundaryIteratorImpl(const RectangularMesh2D& mesh, std::size_t line, std::size_t index): mesh(mesh), line(line), index(index) {}
818
819 void increment() override { ++index; }
820
821 bool equal(const typename BoundaryNodeSetImpl::IteratorImpl& other) const override {
822 return index == static_cast<const BoundaryIteratorImpl&>(other).index;
823 }
824
825 };
826
827 // iterator over vertical line (from bottom to top). for left and right boundaries
828 struct VerticalIteratorImpl: public BoundaryIteratorImpl {
829
830 VerticalIteratorImpl(const RectangularMesh2D& mesh, std::size_t line, std::size_t index): BoundaryIteratorImpl(mesh, line, index) {}
831
832 std::size_t dereference() const override { return this->mesh.index(this->line, this->index); }
833
834 std::unique_ptr<typename BoundaryNodeSetImpl::IteratorImpl> clone() const override {
835 return std::unique_ptr<typename BoundaryNodeSetImpl::IteratorImpl>(new VerticalIteratorImpl(*this));
836 }
837 };
838
839 // iterator over horizonstal line (from left to right), for bottom and top boundaries
840 struct HorizontalIteratorImpl: public BoundaryIteratorImpl {
841
842 HorizontalIteratorImpl(const RectangularMesh2D& mesh, std::size_t line, std::size_t index): BoundaryIteratorImpl(mesh, line, index) {}
843
844 std::size_t dereference() const override { return this->mesh.index(this->index, this->line); }
845
846 std::unique_ptr<typename BoundaryNodeSetImpl::IteratorImpl> clone() const override {
847 return std::unique_ptr<typename BoundaryNodeSetImpl::IteratorImpl>(new HorizontalIteratorImpl(*this));
848 }
849 };
850
851 struct VerticalBoundary: public BoundaryNodeSetWithMeshImpl<RectangularMesh2D> {
852
853 typedef typename BoundaryNodeSetImpl::Iterator Iterator;
854
855 std::size_t line;
856
857 VerticalBoundary(const RectangularMesh2D& mesh, std::size_t line_axis0): BoundaryNodeSetWithMeshImpl<RectangularMesh2D>(mesh), line(line_axis0) {}
858
859 //virtual LeftBoundary* clone() const { return new LeftBoundary(); }
860
861 bool contains(std::size_t mesh_index) const override {
862 return this->mesh.index0(mesh_index) == line;
863 }
864
865 Iterator begin() const override {
866 return Iterator(new VerticalIteratorImpl(this->mesh, line, 0));
867 }
868
869 Iterator end() const override {
870 return Iterator(new VerticalIteratorImpl(this->mesh, line, this->mesh.axis[1]->size()));
871 }
872
873 std::size_t size() const override {
874 return this->mesh.axis[1]->size();
875 }
876 };
877
878
879 struct VerticalBoundaryInRange: public BoundaryNodeSetWithMeshImpl<RectangularMesh2D> {
880
881 typedef typename BoundaryNodeSetImpl::Iterator Iterator;
882
883 std::size_t line, beginInLineIndex, endInLineIndex;
884
885 VerticalBoundaryInRange(const RectangularMesh2D& mesh, std::size_t line_axis0, std::size_t beginInLineIndex, std::size_t endInLineIndex)
886 : BoundaryNodeSetWithMeshImpl<RectangularMesh2D>(mesh), line(line_axis0), beginInLineIndex(beginInLineIndex), endInLineIndex(endInLineIndex) {}
887
888 //virtual LeftBoundary* clone() const { return new LeftBoundary(); }
889
890 bool contains(std::size_t mesh_index) const override {
891 return this->mesh.index0(mesh_index) == line && in_range(this->mesh.index1(mesh_index), beginInLineIndex, endInLineIndex);
892 }
893
894 Iterator begin() const override {
895 return Iterator(new VerticalIteratorImpl(this->mesh, line, beginInLineIndex));
896 }
897
898 Iterator end() const override {
899 return Iterator(new VerticalIteratorImpl(this->mesh, line, endInLineIndex));
900 }
901
902 std::size_t size() const override {
903 return endInLineIndex - beginInLineIndex;
904 }
905 };
906
907 struct HorizontalBoundary: public BoundaryNodeSetWithMeshImpl<RectangularMesh2D> {
908
909 typedef typename BoundaryNodeSetImpl::Iterator Iterator;
910
911 std::size_t line;
912
913 HorizontalBoundary(const RectangularMesh2D& mesh, std::size_t line_axis1): BoundaryNodeSetWithMeshImpl<RectangularMesh2D>(mesh), line(line_axis1) {}
914
915 //virtual TopBoundary* clone() const { return new TopBoundary(); }
916
917 bool contains(std::size_t mesh_index) const override {
918 return this->mesh.index1(mesh_index) == line;
919 }
920
921 Iterator begin() const override {
922 return Iterator(new HorizontalIteratorImpl(this->mesh, line, 0));
923 }
924
925 Iterator end() const override {
926 return Iterator(new HorizontalIteratorImpl(this->mesh, line, this->mesh.axis[0]->size()));
927 }
928
929 std::size_t size() const override {
930 return this->mesh.axis[0]->size();
931 }
932 };
933
934 struct HorizontalBoundaryInRange: public BoundaryNodeSetWithMeshImpl<RectangularMesh2D> {
935
936 typedef typename BoundaryNodeSetImpl::Iterator Iterator;
937
938 std::size_t line, beginInLineIndex, endInLineIndex;
939
940 HorizontalBoundaryInRange(const RectangularMesh2D& mesh, std::size_t line_axis1, std::size_t beginInLineIndex, std::size_t endInLineIndex)
941 : BoundaryNodeSetWithMeshImpl<RectangularMesh2D>(mesh), line(line_axis1), beginInLineIndex(beginInLineIndex), endInLineIndex(endInLineIndex) {}
942 //virtual TopBoundary* clone() const { return new TopBoundary(); }
943
944 bool contains(std::size_t mesh_index) const override {
945 return this->mesh.index1(mesh_index) == line && in_range(this->mesh.index0(mesh_index), beginInLineIndex, endInLineIndex);
946 }
947
948 Iterator begin() const override {
949 return Iterator(new HorizontalIteratorImpl(this->mesh, line, beginInLineIndex));
950 }
951
952 Iterator end() const override {
953 return Iterator(new HorizontalIteratorImpl(this->mesh, line, endInLineIndex));
954 }
955
956 std::size_t size() const override {
957 return endInLineIndex - beginInLineIndex;
958 }
959 };
960
961 //TODO
962 /*struct HorizontalLineBoundary: public BoundaryLogicImpl<RectangularMesh2D> {
963
964 double height;
965
966 bool contains(const RectangularMesh &mesh, std::size_t mesh_index) const {
967 return mesh.index1(mesh_index) == mesh.axis[1]->findNearestIndex(height);
968 }
969 };*/
970
971 public: // boundaries:
972
973 BoundaryNodeSet createVerticalBoundaryAtLine(std::size_t line_nr_axis0) const override;
974
975 BoundaryNodeSet createVerticalBoundaryAtLine(std::size_t line_nr_axis0, std::size_t indexBegin, std::size_t indexEnd) const override;
976
977 BoundaryNodeSet createVerticalBoundaryNear(double axis0_coord) const override;
978
979 BoundaryNodeSet createVerticalBoundaryNear(double axis0_coord, double from, double to) const override;
980
981 BoundaryNodeSet createLeftBoundary() const override;
982
983 BoundaryNodeSet createRightBoundary() const override;
984
985 BoundaryNodeSet createLeftOfBoundary(const Box2D& box) const override;
986
987 BoundaryNodeSet createRightOfBoundary(const Box2D& box) const override;
988
989 BoundaryNodeSet createBottomOfBoundary(const Box2D& box) const override;
990
991 BoundaryNodeSet createTopOfBoundary(const Box2D& box) const override;
992
993 BoundaryNodeSet createHorizontalBoundaryAtLine(std::size_t line_nr_axis1) const override;
994
995 BoundaryNodeSet createHorizontalBoundaryAtLine(std::size_t line_nr_axis1, std::size_t indexBegin, std::size_t indexEnd) const override;
996
997 BoundaryNodeSet createHorizontalBoundaryNear(double axis1_coord) const override;
998
999 BoundaryNodeSet createHorizontalBoundaryNear(double axis1_coord, double from, double to) const override;
1000
1001 BoundaryNodeSet createTopBoundary() const override;
1002
1003 BoundaryNodeSet createBottomBoundary() const override;
1004};
1005
1006
1007 class RectangularMesh2D::ElementMesh: public RectangularMesh2D {
1008
1010 const RectangularMesh2D* originalMesh;
1011
1012 public:
1013 template <typename... Args>
1014 ElementMesh(const RectangularMesh2D* originalMesh, Args... args): RectangularMesh2D(args...), originalMesh(originalMesh) {}
1015
1023 template <typename RandomAccessContainer>
1024 auto interpolateNearestNeighbor(const RandomAccessContainer& data, const Vec<2>& point, const InterpolationFlags& flags) const
1025 -> typename std::remove_reference<decltype(data[0])>::type {
1026 auto p = flags.wrap(point);
1027 prepareNearestNeighborInterpolationForAxis(*originalMesh->axis[0], flags, p.c0, 0);
1028 prepareNearestNeighborInterpolationForAxis(*originalMesh->axis[1], flags, p.c1, 1);
1029 size_t i0 = originalMesh->axis[0]->findUpIndex(p.c0), i1 = originalMesh->axis[1]->findUpIndex(p.c1);
1030 if (i0 == originalMesh->axis[0]->size()) --i0;
1031 if (i0 != 0) --i0;
1032 if (i1 == originalMesh->axis[1]->size()) --i1;
1033 if (i1 != 0) --i1;
1034 return flags.postprocess(point, data[this->index(i0, i1)]);
1035 }
1036
1037 protected:
1038
1039 bool hasSameNodes(const MeshD<2> &to_compare) const override;
1040
1041};
1042
1043
1044 template <typename SrcT, typename DstT>
1045 struct InterpolationAlgorithm<RectangularMesh2D, SrcT, DstT, INTERPOLATION_LINEAR> {
1046 static LazyData<DstT> interpolate(const shared_ptr<const RectangularMesh2D>& src_mesh, const DataVector<const SrcT>& src_vec,
1047 const shared_ptr<const MeshD<2>>& dst_mesh, const InterpolationFlags& flags) {
1048 if (src_mesh->axis[0]->size() == 0 || src_mesh->axis[1]->size() == 0)
1049 throw BadMesh("interpolate", "source mesh empty");
1050 return new LinearInterpolatedLazyDataImpl< DstT, RectangularMesh2D, SrcT >(src_mesh, src_vec, dst_mesh, flags);
1051 }
1052};
1053
1054 template <typename SrcT, typename DstT>
1055 struct InterpolationAlgorithm<RectangularMesh2D, SrcT, DstT, INTERPOLATION_NEAREST> {
1056 static LazyData<DstT> interpolate(const shared_ptr<const RectangularMesh2D>& src_mesh, const DataVector<const SrcT>& src_vec,
1057 const shared_ptr<const MeshD<2>>& dst_mesh, const InterpolationFlags& flags) {
1058 if (src_mesh->axis[0]->size() == 0 || src_mesh->axis[1]->size() == 0)
1059 throw BadMesh("interpolate", "source mesh empty");
1060 return new NearestNeighborInterpolatedLazyDataImpl< DstT, RectangularMesh2D, SrcT >(src_mesh, src_vec, dst_mesh, flags);
1061 }
1062};
1063
1064
1065 template <typename SrcT, typename DstT, InterpolationMethod method>
1066 struct InterpolationAlgorithm<typename std::enable_if<method != INTERPOLATION_DEFAULT, RectangularMesh2D::ElementMesh>::type, SrcT, DstT, method> {
1067 static LazyData<DstT> interpolate(const shared_ptr<const RectangularMesh2D::ElementMesh>& src_mesh, const DataVector<const SrcT>& src_vec,
1068 const shared_ptr<const MeshD<2>>& dst_mesh, const InterpolationFlags& flags) {
1069 return InterpolationAlgorithm<RectangularMesh2D, SrcT, DstT, method>::interpolate(src_mesh, src_vec, dst_mesh, flags);
1070 }
1071};
1072
1073 template <typename SrcT, typename DstT>
1074 struct InterpolationAlgorithm<RectangularMesh2D::ElementMesh, SrcT, DstT, INTERPOLATION_NEAREST> {
1075 static LazyData<DstT> interpolate(const shared_ptr<const RectangularMesh2D::ElementMesh>& src_mesh, const DataVector<const SrcT>& src_vec,
1076 const shared_ptr<const MeshD<2>>& dst_mesh, const InterpolationFlags& flags) {
1077 if (src_mesh->axis[0]->size() == 0 || src_mesh->axis[1]->size() == 0)
1078 throw BadMesh("interpolate", "source mesh empty");
1079 return new NearestNeighborInterpolatedLazyDataImpl<DstT, RectangularMesh2D::ElementMesh, SrcT>(src_mesh, src_vec, dst_mesh, flags);
1080 }
1081};
1082
1083
1089 PLASK_API shared_ptr<RectangularMesh2D> make_rectangular_mesh(const RectangularMesh2D& to_copy);
1090 inline shared_ptr<RectangularMesh2D> make_rectangular_mesh(shared_ptr<const RectangularMesh2D> to_copy) { return make_rectangular_mesh(*to_copy); }
1091
1092} // namespace plask
1093
1094#include "rectangular_spline.hpp"
1095
1096#endif // PLASK__RECTANGULAR2D_H