PLaSK library

Loading...
Searching...
No Matches
rectangular_masked_common.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_MASKED_COMMON_H
15#define PLASK__RECTANGULAR_MASKED_COMMON_H
16
17#include <functional>
18
19#include "rectangular.hpp"
20#include "../utils/numbers_set.hpp"
21
22#include <boost/thread.hpp>
23
24
25 namespace plask {
26
32 template <int DIM>
33 struct RectangularMaskedMeshBase: public RectangularMeshBase<DIM> {
34
36 constexpr static double MIN_DISTANCE = 1e-9; // 1 femtometer
37
39 RectangularMesh<DIM> fullMesh;
40
41 using typename MeshD<DIM>::LocalCoords;
42
43 protected:
44
45 //typedef CompressedSetOfNumbers<std::uint32_t> Set;
46 typedef CompressedSetOfNumbers<std::size_t> Set;
47
49 Set nodeSet;
50
52 Set elementSet;
53
55 struct BoundaryIndexForAxis {
56 std::size_t lo, up;
57 void improveLo(std::size_t i) { if (i < lo) lo = i; }
58 void improveUp(std::size_t i) { if (i > up) up = i; }
59 };
60 typedef BoundaryIndexForAxis BoundaryIndex[DIM];
61 BoundaryIndex boundaryIndex;
62
70 static void findIndexes(const MeshAxis& axis, double wrapped_point_coord, std::size_t& index_lo, std::size_t& index_hi) {
71 index_hi = axis.findUpIndex(wrapped_point_coord);
72 if (index_hi == axis.size()) --index_hi; // p.c0 == axis0->at(axis0->size()-1)
73 else if (index_hi == 0) index_hi = 1;
74 index_lo = index_hi - 1;
75 }
76
84 static std::size_t nearest(double p, const MeshAxis& axis, std::size_t index_lo, std::size_t index_hi) {
85 return p - axis.at(index_lo) <= axis.at(index_hi) - p ? index_lo : index_hi;
86 }
87
91 template <typename MaskedMeshType>
92 struct ElementsBase {
93
94 using Element = typename MaskedMeshType::Element;
95
103 class const_iterator: public Set::ConstIteratorFacade<const_iterator, Element> {
104
105 const MaskedMeshType* maskedMesh;
106
107 Element dereference() const {
108 return Element(*maskedMesh, this->getIndex(), this->getNumber());
109 }
110
111 friend class boost::iterator_core_access;
112
113 public:
114
115 template <typename... CtorArgs>
116 explicit const_iterator(const MaskedMeshType& maskedMesh, CtorArgs&&... ctorArgs)
117 : Set::ConstIteratorFacade<const_iterator, Element>(std::forward<CtorArgs>(ctorArgs)...), maskedMesh(&maskedMesh) {}
118
119 const Set& set() const { return maskedMesh->elementSet; }
120 };
121
123 typedef const_iterator iterator;
124
125 const MaskedMeshType* maskedMesh;
126
127 explicit ElementsBase(const MaskedMeshType& maskedMesh): maskedMesh(&maskedMesh) {}
128
133 std::size_t size() const { return maskedMesh->getElementsCount(); }
134
136 const_iterator begin() const { return const_iterator(*maskedMesh, 0, maskedMesh->elementSet.segments.begin()); }
137
139 const_iterator end() const { return const_iterator(*maskedMesh, size(), maskedMesh->elementSet.segments.end()); }
140
146 Element operator[](std::size_t i) const { return Element(*maskedMesh, i); }
147
148 }; // struct Elements
149
153 template <typename MaskedMeshType>
154 struct ElementMeshBase: MeshD<DIM> {
155
156 using Element = typename MaskedMeshType::Element;
157
159 class const_iterator: public Set::ConstIteratorFacade<const_iterator, LocalCoords> {
160
161 const MaskedMeshType* originalMesh;
162
163 LocalCoords dereference() const {
164 return Element(*originalMesh, this->getIndex(), this->getNumber()).getMidpoint();
165 }
166
167 friend class boost::iterator_core_access;
168
169 public:
170
171 template <typename... CtorArgs>
172 explicit const_iterator(const MaskedMeshType& originalMesh, CtorArgs&&... ctorArgs)
173 : Set::ConstIteratorFacade<const_iterator, LocalCoords>(std::forward<CtorArgs>(ctorArgs)...), originalMesh(&originalMesh) {}
174
175 const Set& set() const { return originalMesh->elementSet; }
176 };
177
179 typedef const_iterator iterator;
180
182 const_iterator begin() const { return const_iterator(*originalMesh, 0, originalMesh->elementSet.segments.begin()); }
183
185 const_iterator end() const { return const_iterator(*originalMesh, size(), originalMesh->elementSet.segments.end()); }
186
187 const MaskedMeshType* originalMesh;
188 RectangularMesh<DIM> fullMesh;
189
190 explicit ElementMeshBase(const MaskedMeshType* originalMesh):
191 originalMesh(originalMesh), fullMesh(*originalMesh->fullMesh.getElementMesh()) {}
192
193 explicit ElementMeshBase(const MaskedMeshType& originalMesh): ElementMeshBase(&originalMesh) {}
194
199 std::size_t size() const override { return originalMesh->getElementsCount(); }
200
201 LocalCoords at(std::size_t index) const override {
202 return fullMesh.at(originalMesh->elementSet.at(index));
203 }
204
205 bool empty() const override { return originalMesh->elementSet.empty(); }
206
207 }; // ElementMeshBase
208
209 void resetBoundyIndex() {
210 for (int d = 0; d < DIM; ++d) { // prepare for finding indexes by subclass constructor:
211 boundaryIndex[d].lo = this->fullMesh.axis[d]->size()-1;
212 boundaryIndex[d].up = 0;
213 }
214 boundaryIndexInitialized = false;
215 }
216
218 void reset() {
219 nodeSet.clear();
220 elementSet.clear();
221 resetBoundyIndex();
222 }
223
224 public:
225
227 enum:std::size_t { NOT_INCLUDED = Set::NOT_INCLUDED };
228
230 RectangularMaskedMeshBase() = default;
231
233 RectangularMaskedMeshBase(const RectangularMesh<DIM>& rectangularMesh, Set nodeSet, bool clone_axes = false)
234 : fullMesh(rectangularMesh, clone_axes), nodeSet(std::move(nodeSet)), elementSetInitialized(false)
235 {
236 resetBoundyIndex();
237 }
238
244 RectangularMaskedMeshBase(const RectangularMesh<DIM>& rectangularMesh, bool clone_axes = false)
245 : fullMesh(rectangularMesh, clone_axes)
246 {
247 resetBoundyIndex();
248 }
249
260 class const_iterator: public CompressedSetOfNumbers<std::size_t>::ConstIteratorFacade<const_iterator, LocalCoords> {
261
262 friend class boost::iterator_core_access;
263
264 const RectangularMaskedMeshBase* mesh;
265
266 LocalCoords dereference() const {
267 return mesh->fullMesh.at(this->getNumber());
268 }
269
270 public:
271
272 template <typename... CtorArgs>
273 explicit const_iterator(const RectangularMaskedMeshBase& mesh, CtorArgs&&... ctorArgs)
274 : CompressedSetOfNumbers<std::size_t>::ConstIteratorFacade<const_iterator, LocalCoords>(std::forward<CtorArgs>(ctorArgs)...), mesh(&mesh) {}
275
276 const CompressedSetOfNumbers<std::size_t>& set() const { return mesh->nodeSet; }
277 };
278
280 typedef const_iterator iterator;
281
282 const_iterator begin() const { return const_iterator(*this, 0, nodeSet.segments.begin()); }
283 const_iterator end() const { return const_iterator(*this, size(), nodeSet.segments.end()); }
284
285 LocalCoords at(std::size_t index) const override {
286 return fullMesh.at(nodeSet.at(index));
287 }
288
289 std::size_t size() const override { return nodeSet.size(); }
290
291 bool empty() const override { return nodeSet.empty(); }
292
296 bool full() const { return nodeSet.size() == fullMesh.size(); }
297
303 void selectAll(const RectangularMesh<DIM>& rectangularMesh, bool clone_axes = false) {
304 fullMesh.reset(rectangularMesh, clone_axes);
305 selectAll();
306 }
307
311 void selectAll() {
312 this->nodeSet.assignRange(fullMesh.size());
313 this->elementSet.assignRange(fullMesh.getElementsCount());
314 elementSetInitialized = true;
315 for (int d = 0; d < DIM; ++d) {
316 boundaryIndex[d].lo = 0;
317 boundaryIndex[d].up = fullMesh.axis[d]->size()-1;
318 }
319 boundaryIndexInitialized = true;
320 }
321
327 inline std::size_t index(const Vec<DIM, std::size_t>& indexes) const {
328 return nodeSet.indexOf(fullMesh.index(indexes));
329 }
330
336 inline std::size_t index0(std::size_t mesh_index) const {
337 return fullMesh.index0(nodeSet.at(mesh_index));
338 }
339
345 inline std::size_t index1(std::size_t mesh_index) const {
346 return fullMesh.index1(nodeSet.at(mesh_index));
347 }
348
354 inline Vec<DIM, std::size_t> indexes(std::size_t mesh_index) const {
355 return fullMesh.indexes(nodeSet.at(mesh_index));
356 }
357
363 inline std::size_t majorIndex(std::size_t mesh_index) const {
364 return fullMesh.majorIndex(nodeSet.at(mesh_index));
365 }
366
372 inline std::size_t minorIndex(std::size_t mesh_index) const {
373 return fullMesh.minorIndex(nodeSet.at(mesh_index));
374 }
375
380 std::size_t getElementsCount0() const {
381 return fullMesh.getElementsCount0();
382 }
383
388 std::size_t getElementsCount1() const {
389 return fullMesh.getElementsCount1();
390 }
391
396 std::size_t getElementsCount() const {
397 return ensureHasElements().size();
398 }
399
405 std::size_t getElementIndexFromLowIndex(std::size_t mesh_index_of_el_bottom_left) const {
406 return ensureHasElements().indexOf(fullMesh.getElementIndexFromLowIndex(nodeSet.at(mesh_index_of_el_bottom_left)));
407 }
408
414 std::size_t getElementMeshLowIndex(std::size_t element_index) const {
415 return nodeSet.indexOf(fullMesh.getElementMeshLowIndex(ensureHasElements().at(element_index)));
416 }
417
424 Vec<DIM, std::size_t> getElementMeshLowIndexes(std::size_t element_index) const {
425 return fullMesh.getElementMeshLowIndexes(ensureHasElements().at(element_index));
426 }
427
433 double getElementArea(std::size_t element_index) const {
434 return fullMesh.getElementArea(ensureHasElements().at(element_index));
435 }
436
442 double getElementMidpoint0(std::size_t index0) const {
443 return fullMesh.getElementMidpoint0(index0);
444 }
445
451 double getElementMidpoint1(std::size_t index1) const {
452 return fullMesh.getElementMidpoint1(index1);
453 }
454
460 Vec<DIM, double> getElementMidpoint(std::size_t element_index) const {
461 return fullMesh.getElementMidpoint(ensureHasElements().at(element_index));
462 }
463
469 typename Primitive<DIM>::Box getElementBox(std::size_t element_index) const {
470 return fullMesh.getElementBox(ensureHasElements().at(element_index));
471 }
472
473 protected: // constructing elementSet from nodes set (element is chosen when all its vertices are chosen) on-deamand
474
476 DontCopyThisField<boost::mutex> writeMutex;
477
479 bool elementSetInitialized = true;
480
482 bool boundaryIndexInitialized;
483
484 private:
485
486 /*bool restVerticesIncluded(const RectangularMesh2D::Element& el) const {
487 return //nodeSet.includes(el.getLoLoIndex()) &&
488 nodeSet.includes(el.getUpLoIndex()) &&
489 nodeSet.includes(el.getLoUpIndex()) &&
490 nodeSet.includes(el.getUpUpIndex());
491 }
492 bool restVerticesIncluded(const RectangularMesh3D::Element& el) const {
493 return //nodeSet.includes(el.getLoLoLoIndex()) &&
494 nodeSet.includes(el.getUpLoLoIndex()) &&
495 nodeSet.includes(el.getLoUpLoIndex()) &&
496 nodeSet.includes(el.getLoLoUpIndex()) &&
497 nodeSet.includes(el.getLoUpUpIndex()) &&
498 nodeSet.includes(el.getUpLoUpIndex()) &&
499 nodeSet.includes(el.getUpUpLoIndex()) &&
500 nodeSet.includes(el.getUpUpUpIndex());
501 }
502 void calculateElements() {
503 boost::lock_guard<boost::mutex> lock((boost::mutex&)writeElementSet);
504 if (elementSetInitialized) return; // another thread has initialized elementSet just when we waited for mutex
505 // TODO faster implementation
506 for (auto lowNode: nodeSet) {
507 if (!fullMesh.isLowIndexOfElement(lowNode)) continue;
508 auto elementIndex = fullMesh.getElementIndexFromLowIndex(lowNode);
509 if (restVerticesIncluded(fullMesh.getElement(elementIndex))) elementSet.push_back(elementIndex);
510 }
511 elementSetInitialized = true;
512 }*/ // ^ generic algorithm, slow but probably correct
513
514 template <int d = DIM>
515 typename std::enable_if<d == 2>::type calculateElements() {
516 boost::lock_guard<boost::mutex> lock((boost::mutex&)writeMutex);
517 if (elementSetInitialized) return; // another thread has initialized elementSet just when we waited for mutex
518
519 if (fullMesh.axis[0]->size() <= 1 || fullMesh.axis[1]->size() <= 1) {
520 elementSetInitialized = true;
521 return;
522 }
523
524 elementSet = nodeSet.transformed([] (std::size_t&, std::size_t& e) { --e; }); // same as nodeSet.intersected(nodeSet.shiftedLeft(1))
525 auto minor_axis_size = fullMesh.minorAxis()->size();
526 elementSet = elementSet.intersection(elementSet.shiftedLeft(minor_axis_size));
527 // now elementSet includes all low indexes which have other corners of elements (plus some indexes in the last column)
528 // we have to transform low indexes to indexes of elements:
529 elementSet = elementSet.transformed([&, minor_axis_size] (std::size_t& b, std::size_t& e) { // here: 0 <= b < e
530 if (e % minor_axis_size == 0) --e; // end of segment cannot lie at the last column, as getElementIndexFromLowIndex confuses last column with the first element in the next row
531 // no need to fix b (if b % minor_axis_size == 0) as rounding in getElementIndexFromLowIndex will give good value (as getElementIndexFromLowIndex(b)==getElementIndexFromLowIndex(b+1))
532 b = fullMesh.getElementIndexFromLowIndex(b);
533 e = fullMesh.getElementIndexFromLowIndex(e);
534 });
535
536 elementSetInitialized = true;
537 }
538
539 template <int d = DIM>
540 typename std::enable_if<d == 3>::type calculateElements() {
541 boost::lock_guard<boost::mutex> lock((boost::mutex&)writeMutex);
542 if (elementSetInitialized) return; // another thread has initialized elementSet just when we waited for mutex
543
544 if (fullMesh.axis[0]->size() <= 1 || fullMesh.axis[1]->size() <= 1 || fullMesh.axis[2]->size() <= 1) {
545 elementSetInitialized = true;
546 return;
547 }
548
549 elementSet = nodeSet.transformed([] (std::size_t&, std::size_t& e) { --e; }); // same as nodeSet.intersected(nodeSet.shiftedLeft(1))
550 auto minor_axis_size = fullMesh.minorAxis()->size();
551 elementSet = elementSet.intersection(elementSet.shiftedLeft(minor_axis_size));
552 auto medium_axis_size = fullMesh.mediumAxis()->size();
553 elementSet = elementSet.intersection(elementSet.shiftedLeft(minor_axis_size*medium_axis_size));
554 // now elementSet includes all low indexes which have other corners of elements (plus some indexes in the last column)
555 // we have to transform low indexes to indexes of elements:
556 elementSet = elementSet.transformed([&, minor_axis_size, medium_axis_size] (std::size_t& b, std::size_t& e) { // here: 0 <= b < e
557 const std::size_t b_div = b / minor_axis_size;
558 if (b_div % medium_axis_size == (medium_axis_size-1))
559 b = (b_div+1) * minor_axis_size; // first index in the next plane
560 b = fullMesh.getElementIndexFromLowIndex(b);
561
562 const std::size_t e_div = (e-1) / minor_axis_size; // e-1 is the last index of the range
563 if (e_div % medium_axis_size == (medium_axis_size-1))
564 e = e_div * minor_axis_size - 1; // -1 to not be divisible by minor_axis_size
565 else if (e % minor_axis_size == 0) --e;
566 e = fullMesh.getElementIndexFromLowIndex(e);
567 });
568
569 elementSetInitialized = true;
570 }
571
572 template <int d = DIM>
573 typename std::enable_if<d == 2>::type calculateBoundaryIndex() {
574 boost::lock_guard<boost::mutex> lock((boost::mutex&)writeMutex);
575 if (boundaryIndexInitialized) return; // another thread has initialized boundaryIndex just when we waited for mutex
576
577 const auto minor = fullMesh.minorAxisIndex();
578 const auto major = fullMesh.majorAxisIndex();
579 nodeSet.forEachSegment([this, minor, major] (std::size_t b, std::size_t e) {
580 const auto indexes_f = fullMesh.indexes(b);
581 const auto indexes_l = fullMesh.indexes(e-1);
582 if (indexes_f[major] != indexes_l[major]) {
583 boundaryIndex[minor].lo = 0;
584 boundaryIndex[minor].up = fullMesh.minorAxis()->size()-1;
585 } else {
586 boundaryIndex[minor].improveLo(indexes_f[minor]);
587 boundaryIndex[minor].improveUp(indexes_l[minor]);
588 }
589 boundaryIndex[major].improveLo(indexes_f[major]);
590 boundaryIndex[major].improveUp(indexes_l[major]);
591 });
592
593 boundaryIndexInitialized = true;
594 }
595
596 template <int d = DIM>
597 typename std::enable_if<d == 3>::type calculateBoundaryIndex() {
598 boost::lock_guard<boost::mutex> lock((boost::mutex&)writeMutex);
599 if (boundaryIndexInitialized) return; // another thread has initialized boundaryIndex just when we waited for mutex
600
601 const auto minor = fullMesh.minorAxisIndex();
602 const auto medium = fullMesh.mediumAxisIndex();
603 const auto major = fullMesh.majorAxisIndex();
604 nodeSet.forEachSegment([this, minor, medium, major] (std::size_t b, std::size_t e) {
605 const auto indexes_f = fullMesh.indexes(b);
606 const auto indexes_l = fullMesh.indexes(e-1);
607 if (indexes_f[major] != indexes_l[major]) {
608 boundaryIndex[minor].lo = 0;
609 boundaryIndex[minor].up = fullMesh.minorAxis()->size()-1;
610 boundaryIndex[medium].lo = 0;
611 boundaryIndex[medium].up = fullMesh.mediumAxis()->size()-1;
612 } else {
613 if (indexes_f[medium] != indexes_l[medium]) {
614 boundaryIndex[minor].lo = 0;
615 boundaryIndex[minor].up = fullMesh.minorAxis()->size()-1;
616 } else { // here: indexes_f[minor] <= indexes_l[minor]
617 boundaryIndex[minor].improveLo(indexes_f[minor]);
618 boundaryIndex[minor].improveUp(indexes_l[minor]);
619 }
620 // indexes_f[major] == indexes_l[major] => indexes_f[medium] <= indexes_l[medium]
621 boundaryIndex[medium].improveLo(indexes_f[medium]);
622 boundaryIndex[medium].improveUp(indexes_l[medium]);
623 }
624 boundaryIndex[major].improveLo(indexes_f[major]);
625 boundaryIndex[major].improveUp(indexes_l[major]);
626 });
627
628 boundaryIndexInitialized = true;
629 }
630
631 protected:
636 const Set& ensureHasElements() const {
637 if (!elementSetInitialized) const_cast<RectangularMaskedMeshBase<DIM>*>(this)->calculateElements();
638 return elementSet;
639 }
640
645 const BoundaryIndex& ensureHasBoundaryIndex() const {
646 if (!boundaryIndexInitialized) const_cast<RectangularMaskedMeshBase<DIM>*>(this)->calculateBoundaryIndex();
647 return boundaryIndex;
648 }
649};
650
651
652} // namespace plask
653
654#endif // PLASK__RECTANGULAR_MASKED_COMMON_H