PLaSK library

Loading...
Searching...
No Matches
rectangular_masked3d.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 "rectangular_masked3d.hpp"
15
16 namespace plask {
17
18 void RectangularMaskedMesh3D::reset(const RectangularMaskedMesh3D::Predicate &predicate) {
19 RectangularMaskedMeshBase<3>::reset();
20 initNodesAndElements(predicate);
21}
22
23 RectangularMaskedMesh3D::RectangularMaskedMesh3D(const RectangularMesh<3> &rectangularMesh, const RectangularMaskedMesh3D::Predicate &predicate, bool clone_axes)
24 : RectangularMaskedMeshBase(rectangularMesh, clone_axes)
25{
26 initNodesAndElements(predicate);
27}
28
29 void RectangularMaskedMesh3D::reset(const RectangularMesh<3> &rectangularMesh, const RectangularMaskedMesh3D::Predicate &predicate, bool clone_axes)
30{
31 this->fullMesh.reset(rectangularMesh, clone_axes);
32 reset(predicate);
33}
34
35 RectangularMaskedMesh3D::RectangularMaskedMesh3D(const RectangularMesh<DIM> &rectangularMesh, RectangularMaskedMeshBase::Set nodeSet, bool clone_axes)
36 : RectangularMaskedMeshBase(rectangularMesh, std::move(nodeSet), clone_axes)
37{
38}
39
40 void RectangularMaskedMesh3D::initNodesAndElements(const RectangularMaskedMesh3D::Predicate &predicate)
41{
42 for (auto el_it = this->fullMesh.elements().begin(); el_it != this->fullMesh.elements().end(); ++el_it)
43 if (predicate(*el_it)) {
44 elementSet.push_back(el_it.index);
45 nodeSet.insert(el_it->getLoLoLoIndex());
46
47 nodeSet.insert(el_it->getUpLoLoIndex());
48 nodeSet.insert(el_it->getLoUpLoIndex());
49 nodeSet.insert(el_it->getLoLoUpIndex());
50
51 nodeSet.insert(el_it->getLoUpUpIndex());
52 nodeSet.insert(el_it->getUpLoUpIndex());
53 nodeSet.insert(el_it->getUpUpLoIndex());
54
55 nodeSet.push_back(el_it->getUpUpUpIndex());
56 /*boundaryIndex[0].improveLo(el_it->getLowerIndex0());
57 boundaryIndex[0].improveUp(el_it->getUpperIndex0());
58 boundaryIndex[1].improveLo(el_it->getLowerIndex1());
59 boundaryIndex[1].improveUp(el_it->getUpperIndex1());
60 boundaryIndex[2].improveLo(el_it->getLowerIndex2());
61 boundaryIndex[2].improveUp(el_it->getUpperIndex2());*/ // this is initialized lazy
62 }
63 nodeSet.shrink_to_fit();
64 elementSet.shrink_to_fit();
65 elementSetInitialized = true;
66}
67
68 bool RectangularMaskedMesh3D::prepareInterpolation(const Vec<3> &point, Vec<3> &wrapped_point,
69 std::size_t& index0_lo, std::size_t& index0_hi,
70 std::size_t& index1_lo, std::size_t& index1_hi,
71 std::size_t& index2_lo, std::size_t& index2_hi,
72 const InterpolationFlags &flags) const {
73 wrapped_point = flags.wrap(point);
74
75 if (!canBeIncluded(wrapped_point)) return false;
76
77 findIndexes(*fullMesh.axis[0], wrapped_point.c0, index0_lo, index0_hi);
78 findIndexes(*fullMesh.axis[1], wrapped_point.c1, index1_lo, index1_hi);
79 findIndexes(*fullMesh.axis[2], wrapped_point.c2, index2_lo, index2_hi);
80 assert(index0_hi == index0_lo + 1);
81 assert(index1_hi == index1_lo + 1);
82 assert(index2_hi == index2_lo + 1);
83
84 double lo0 = fullMesh.axis[0]->at(index0_lo), hi0 = fullMesh.axis[0]->at(index0_hi),
85 lo1 = fullMesh.axis[1]->at(index1_lo), hi1 = fullMesh.axis[1]->at(index1_hi),
86 lo2 = fullMesh.axis[2]->at(index2_lo), hi2 = fullMesh.axis[2]->at(index2_hi);
87
88 ensureHasElements();
89 for (char i2 = 0; i2 < 2; ++i2) {
90 for (char i1 = 0; i1 < 2; ++i1) {
91 for (char i0 = 0; i0 < 2; ++i0) {
92 if (elementSet.includes(fullMesh.getElementIndexFromLowIndexes(index0_lo, index1_lo, index2_lo))) {
93 index0_hi = index0_lo + 1; index1_hi = index1_lo + 1; index2_hi = index2_lo + 1;
94 return true;
95 }
96 if (index0_lo > 0 && lo0 <= wrapped_point.c0 && wrapped_point.c0 < lo0+MIN_DISTANCE) index0_lo = index0_hi-2;
97 else if (index0_lo < fullMesh.axis[0]->size()-2 && hi0-MIN_DISTANCE < wrapped_point.c0 && wrapped_point.c0 <= hi0) index0_lo = index0_hi;
98 else break;
99 }
100 index0_lo = index0_hi - 1;
101 if (index1_lo > 0 && lo1 <= wrapped_point.c1 && wrapped_point.c1 < lo1+MIN_DISTANCE) index1_lo = index1_hi-2;
102 else if (index1_lo < fullMesh.axis[1]->size()-2 && hi1-MIN_DISTANCE < wrapped_point.c1 && wrapped_point.c1 <= hi1) index1_lo = index1_hi;
103 else break;
104 }
105 index1_lo = index1_hi - 1;
106 if (index2_lo > 0 && lo2 <= wrapped_point.c2 && wrapped_point.c2 < lo2+MIN_DISTANCE) index2_lo = index2_hi-2;
107 else if (index2_lo < fullMesh.axis[2]->size()-2 && hi2-MIN_DISTANCE < wrapped_point.c2 && wrapped_point.c2 <= hi1) index2_lo = index2_hi;
108 else break;
109 }
110
111 return false;
112}
113
114 BoundaryNodeSet RectangularMaskedMesh3D::createIndex0BoundaryAtLine(std::size_t line_nr_axis0, std::size_t index1Begin, std::size_t index1End, std::size_t index2Begin, std::size_t index2End) const
115{
116 if (this->fullMesh.isChangeSlower(1, 2))
117 return new BoundaryNodeSetImpl<1, 2>(*this, line_nr_axis0, index1Begin, index2Begin, index1End, index2End);
118 else
119 return new BoundaryNodeSetImpl<2, 1>(*this, line_nr_axis0, index1Begin, index2Begin, index2End, index1End);
120}
121
122 BoundaryNodeSet RectangularMaskedMesh3D::createIndex0BoundaryAtLine(std::size_t line_nr_axis0) const {
123 ensureHasBoundaryIndex();
124 return createIndex0BoundaryAtLine(line_nr_axis0, boundaryIndex[1].lo, boundaryIndex[1].up+1, boundaryIndex[2].lo, boundaryIndex[2].up+1);
125}
126
127 BoundaryNodeSet RectangularMaskedMesh3D::createIndex1BoundaryAtLine(std::size_t line_nr_axis1, std::size_t index0Begin, std::size_t index0End, std::size_t index2Begin, std::size_t index2End) const
128{
129 if (this->fullMesh.isChangeSlower(0, 2))
130 return new BoundaryNodeSetImpl<0, 2>(*this, index0Begin, line_nr_axis1, index2Begin, index0End, index2End);
131 else
132 return new BoundaryNodeSetImpl<2, 0>(*this, index0Begin, line_nr_axis1, index2Begin, index2End, index0End);
133}
134
135 BoundaryNodeSet RectangularMaskedMesh3D::createIndex1BoundaryAtLine(std::size_t line_nr_axis1) const {
136 ensureHasBoundaryIndex();
137 return createIndex1BoundaryAtLine(line_nr_axis1, boundaryIndex[0].lo, boundaryIndex[0].up+1, boundaryIndex[2].lo, boundaryIndex[2].up+1);
138}
139
140 BoundaryNodeSet RectangularMaskedMesh3D::createIndex2BoundaryAtLine(std::size_t line_nr_axis2, std::size_t index0Begin, std::size_t index0End, std::size_t index1Begin, std::size_t index1End) const
141{
142 if (this->fullMesh.isChangeSlower(0, 1))
143 return new BoundaryNodeSetImpl<0, 1>(*this, index0Begin, index1Begin, line_nr_axis2, index0End, index1End);
144 else
145 return new BoundaryNodeSetImpl<1, 0>(*this, index0Begin, index1Begin, line_nr_axis2, index1End, index0End);
146}
147
148 BoundaryNodeSet RectangularMaskedMesh3D::createIndex2BoundaryAtLine(std::size_t line_nr_axis2) const {
149 ensureHasBoundaryIndex();
150 return createIndex2BoundaryAtLine(line_nr_axis2, boundaryIndex[0].lo, boundaryIndex[0].up+1, boundaryIndex[1].lo, boundaryIndex[1].up+1);
151}
152
153 BoundaryNodeSet RectangularMaskedMesh3D::createBackBoundary() const {
154 return createIndex0BoundaryAtLine(ensureHasBoundaryIndex()[0].lo);
155}
156
157 BoundaryNodeSet RectangularMaskedMesh3D::createFrontBoundary() const {
158 return createIndex0BoundaryAtLine(ensureHasBoundaryIndex()[0].up);
159}
160
161 BoundaryNodeSet RectangularMaskedMesh3D::createLeftBoundary() const {
162 return createIndex1BoundaryAtLine(ensureHasBoundaryIndex()[1].lo);
163}
164
165 BoundaryNodeSet RectangularMaskedMesh3D::createRightBoundary() const {
166 return createIndex1BoundaryAtLine(ensureHasBoundaryIndex()[1].up);
167}
168
169 BoundaryNodeSet RectangularMaskedMesh3D::createBottomBoundary() const {
170 return createIndex2BoundaryAtLine(ensureHasBoundaryIndex()[2].lo);
171}
172
173 BoundaryNodeSet RectangularMaskedMesh3D::createTopBoundary() const {
174 return createIndex2BoundaryAtLine(ensureHasBoundaryIndex()[2].up);
175}
176
177 BoundaryNodeSet RectangularMaskedMesh3D::createBackOfBoundary(const Box3D &box) const {
178 std::size_t line, begInd1, endInd1, begInd2, endInd2;
179 if (details::getLineLo(line, *fullMesh.axis[0], box.lower.c0, box.upper.c0) &&
180 details::getIndexesInBounds(begInd1, endInd1, *fullMesh.axis[1], box.lower.c1, box.upper.c1) &&
181 details::getIndexesInBounds(begInd2, endInd2, *fullMesh.axis[2], box.lower.c2, box.upper.c2))
182 return createIndex0BoundaryAtLine(line, begInd1, endInd1, begInd2, endInd2);
183 else
184 return new EmptyBoundaryImpl();
185}
186
187 BoundaryNodeSet RectangularMaskedMesh3D::createFrontOfBoundary(const Box3D &box) const {
188 std::size_t line, begInd1, endInd1, begInd2, endInd2;
189 if (details::getLineHi(line, *fullMesh.axis[0], box.lower.c0, box.upper.c0) &&
190 details::getIndexesInBounds(begInd1, endInd1, *fullMesh.axis[1], box.lower.c1, box.upper.c1) &&
191 details::getIndexesInBounds(begInd2, endInd2, *fullMesh.axis[2], box.lower.c2, box.upper.c2))
192 return createIndex0BoundaryAtLine(line, begInd1, endInd1, begInd2, endInd2);
193 else
194 return new EmptyBoundaryImpl();
195}
196
197 BoundaryNodeSet RectangularMaskedMesh3D::createLeftOfBoundary(const Box3D &box) const {
198 std::size_t line, begInd0, endInd0, begInd2, endInd2;
199 if (details::getLineLo(line, *fullMesh.axis[1], box.lower.c1, box.upper.c1) &&
200 details::getIndexesInBounds(begInd0, endInd0, *fullMesh.axis[0], box.lower.c0, box.upper.c0) &&
201 details::getIndexesInBounds(begInd2, endInd2, *fullMesh.axis[2], box.lower.c2, box.upper.c2))
202 return createIndex1BoundaryAtLine(line, begInd0, endInd0, begInd2, endInd2);
203 else
204 return new EmptyBoundaryImpl();
205}
206
207 BoundaryNodeSet RectangularMaskedMesh3D::createRightOfBoundary(const Box3D &box) const {
208 std::size_t line, begInd0, endInd0, begInd2, endInd2;
209 if (details::getLineHi(line, *fullMesh.axis[1], box.lower.c1, box.upper.c1) &&
210 details::getIndexesInBounds(begInd0, endInd0, *fullMesh.axis[0], box.lower.c0, box.upper.c0) &&
211 details::getIndexesInBounds(begInd2, endInd2, *fullMesh.axis[2], box.lower.c2, box.upper.c2))
212 return createIndex1BoundaryAtLine(line, begInd0, endInd0, begInd2, endInd2);
213 else
214 return new EmptyBoundaryImpl();
215}
216
217 BoundaryNodeSet RectangularMaskedMesh3D::createBottomOfBoundary(const Box3D &box) const {
218 std::size_t line, begInd0, endInd0, begInd1, endInd1;
219 if (details::getLineLo(line, *fullMesh.axis[2], box.lower.c2, box.upper.c2) &&
220 details::getIndexesInBounds(begInd0, endInd0, *fullMesh.axis[0], box.lower.c0, box.upper.c0) &&
221 details::getIndexesInBounds(begInd1, endInd1, *fullMesh.axis[1], box.lower.c1, box.upper.c1))
222 return createIndex2BoundaryAtLine(line, begInd0, endInd0, begInd1, endInd1);
223 else
224 return new EmptyBoundaryImpl();
225}
226
227 BoundaryNodeSet RectangularMaskedMesh3D::createTopOfBoundary(const Box3D &box) const {
228 std::size_t line, begInd0, endInd0, begInd1, endInd1;
229 if (details::getLineHi(line, *fullMesh.axis[2], box.lower.c2, box.upper.c2) &&
230 details::getIndexesInBounds(begInd0, endInd0, *fullMesh.axis[0], box.lower.c0, box.upper.c0) &&
231 details::getIndexesInBounds(begInd1, endInd1, *fullMesh.axis[1], box.lower.c1, box.upper.c1))
232 return createIndex2BoundaryAtLine(line, begInd0, endInd0, begInd1, endInd1);
233 else
234 return new EmptyBoundaryImpl();
235}
236
237} // namespace plask