PLaSK library

Loading...
Searching...
No Matches
ordered1d.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__RECTILINEAR1D_H
15#define PLASK__RECTILINEAR1D_H
16
21#include <vector>
22#include <algorithm>
23#include <initializer_list>
24
25#include "mesh.hpp"
26#include "../vec.hpp"
27#include "../log/log.hpp"
28#include "../utils/iterators.hpp"
29#include "../utils/interpolation.hpp"
30
31#include "axis1d.hpp"
32
33 namespace plask {
34
38 class PLASK_API OrderedAxis: public MeshAxis {
39
41 std::vector<double> points;
42
43 void sortPointsAndRemoveNonUnique(double min_dist);
44
45 public:
46
47 struct WarningOff {
48 OrderedAxis* axis;
49 bool prev_state;
50 WarningOff(OrderedAxis& axis): axis(&axis), prev_state(axis.warn_too_close) { axis.warn_too_close = false; }
51 WarningOff(OrderedAxis* axis): axis(axis), prev_state(axis->warn_too_close) { axis->warn_too_close = false; }
52 WarningOff(const shared_ptr<OrderedAxis>& axis): axis(axis.get()), prev_state(axis->warn_too_close) { axis->warn_too_close = false; }
53 ~WarningOff() { axis->warn_too_close = prev_state; }
54 };
55
57 bool warn_too_close;
58
60 constexpr static double MIN_DISTANCE = 1e-6; // 1 picometer
61
63 typedef double PointType;
64
66 typedef std::vector<double>::const_iterator native_const_iterator;
67
69 native_const_iterator begin() const { return points.begin(); }
70
72 native_const_iterator end() const { return points.end(); }
73
75 const std::vector<double>& getPointsVector() const { return points; }
76
85 native_const_iterator find(double to_find) const;
86
87 std::size_t findIndex(double to_find) const override { return std::size_t(find(to_find) - begin()); }
88
95 native_const_iterator findUp(double to_find) const;
96
97 std::size_t findUpIndex(double to_find) const override { return std::size_t(findUp(to_find) - begin()); }
98
104 native_const_iterator findNearest(double to_find) const;
105
106 std::size_t findNearestIndex(double to_find) const override { return findNearest(to_find) - begin(); }
107
109 OrderedAxis(): warn_too_close(true) {}
110
112 OrderedAxis(const OrderedAxis& src): points(src.points), warn_too_close(true) {}
113
115 OrderedAxis(OrderedAxis&& src): points(std::move(src.points)), warn_too_close(true) {}
116
118 OrderedAxis(const MeshAxis& src): points(src.size()), warn_too_close(true) {
119 if (src.isIncreasing())
120 std::copy(src.begin(), src.end(), points.begin());
121 else
122 std::reverse_copy(src.begin(), src.end(), points.begin());
123 }
124
131 OrderedAxis(std::initializer_list<PointType> points, double min_dist=MIN_DISTANCE);
132
139 OrderedAxis(const std::vector<PointType>& points, double min_dist=MIN_DISTANCE);
140
147 OrderedAxis(std::vector<PointType>&& points, double min_dist=MIN_DISTANCE);
148
150 OrderedAxis& operator=(const OrderedAxis& src);
151
153 OrderedAxis& operator=(OrderedAxis&& src);
154
156 OrderedAxis& operator=(const MeshAxis& src);
157
164 bool operator==(const OrderedAxis& to_compare) const;
165
166 void writeXML(XMLElement& object) const override;
167
169 std::size_t size() const override { return points.size(); }
170
171 double at(std::size_t index) const override { assert(index < points.size()); return points[index]; }
172
173 // @return true only if there are no points in mesh
174 //bool empty() const { return points.empty(); }
175
184 bool addPoint(double new_node_cord, double min_dist);
185
193 bool addPoint(double new_node_cord) {
194 return addPoint(new_node_cord, MIN_DISTANCE);
195 }
196
201 void removePoint(std::size_t index);
202
208 void removePoints(std::size_t start, std::size_t stop);
209
216 void removePoints(std::size_t start, std::size_t stop, std::ptrdiff_t step);
217
226 template <typename IteratorT>
227 void addOrderedPoints(IteratorT begin, IteratorT end, std::size_t points_count_hint, double min_dist=MIN_DISTANCE);
228
235 //TODO use iterator traits and write version for input iterator
236 template <typename RandomAccessIteratorT>
237 void addOrderedPoints(RandomAccessIteratorT begin, RandomAccessIteratorT end) { addOrderedPoints(begin, end, end - begin); }
238
246 void addPointsLinear(double first, double last, std::size_t points_count);
247
253 //const double& operator[](std::size_t index) const { return points[index]; }
254
258 void clear();
259
260 shared_ptr<MeshAxis> clone() const override;
261
262 bool isIncreasing() const override { return true; }
263
270 template <typename RandomAccessContainer>
271 auto interpolateLinear(const RandomAccessContainer& data, double point) const -> typename std::remove_reference<decltype(data[0])>::type;
272
273};
274
275// OrderedAxis method templates implementation
276 template <typename RandomAccessContainer>
277 inline auto OrderedAxis::interpolateLinear(const RandomAccessContainer& data, double point) const -> typename std::remove_reference<decltype(data[0])>::type {
278 std::size_t index = findIndex(point);
279 if (index == size()) return data[index-1]; //TODO what should it do if mesh is empty?
280 if (index == 0 || points[index] == point) return data[index]; // hit exactly
281 // here: points[index-1] < point < points[index]
282 return interpolation::linear(points[index-1], data[index-1], points[index], data[index], point);
283}
284
285 template <typename IteratorT>
286 inline void OrderedAxis::addOrderedPoints(IteratorT begin, IteratorT end, std::size_t points_count_hint, double min_dist) {
287 std::vector<double> result;
288 result.reserve(this->size() + points_count_hint);
289 std::set_union(this->points.begin(), this->points.end(), begin, end, std::back_inserter(result));
290 this->points = std::move(result);
291 // Remove points too close to each other
292 auto almost_equal = [min_dist, this](const double& x, const double& y) -> bool {
293 bool remove = std::abs(x-y) < min_dist;
294 if (warn_too_close && remove) writelog(LOG_WARNING, "Points in ordered mesh too close, skipping point at {0}", y);
295 return remove;
296 };
297 this->points.erase(std::unique(this->points.begin(), this->points.end(), almost_equal), this->points.end());
298 fireResized();
299}
300
301 typedef OrderedAxis OrderedMesh1D;
302
303 shared_ptr<OrderedMesh1D> readRectilinearMeshAxis(XMLReader& reader);
304
305} // namespace plask
306
307#endif // PLASK__RECTILINEAR1D_H