PLaSK library

Loading...
Searching...
No Matches
axis1d.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 "axis1d.hpp"
15#include "ordered1d.hpp"
16
17#include "../utils/stl.hpp"
18
19 namespace plask {
20
21
22//enable_shared_from_this for Mesh (for getMidpointAxis impl. and change to shared_ptr) ???
23 class MidpointAxis: public MeshAxis {
24
25 //shared_ptr<MeshAxis> wrapped;
26 const MeshAxis& wrapped;
27
28 public:
29
30 //MidpointAxis(shared_ptr<const MeshAxis> wrapped = nullptr): wrapped(nullptr) { setWrapped(wrapped); }
31 MidpointAxis(const MeshAxis& wrapped): wrapped(wrapped) { }
32
33 //shared_ptr<const MeshAxis > getWrapped() const;
34
35 //void setWrapped(shared_ptr<const MeshAxis > wrapped);
36
37 //void clear() override { setWrapped(nullptr); }
38
39 std::size_t size() const override;
40
41 double at(std::size_t index) const override;
42
43 bool isIncreasing() const override;
44};
45
46std::size_t MidpointAxis::size() const {
47 std::size_t wrapped_size = wrapped.size();
48 return wrapped_size ? wrapped_size - 1 : 0;
49}
50
51 double MidpointAxis::at(std::size_t index) const {
52 return (wrapped.at(index) + wrapped.at(index+1)) * 0.5;
53}
54
55 bool MidpointAxis::isIncreasing() const {
56 return wrapped.isIncreasing();
57}
58
59
60// -------------- MeshAxis ---------------------------------------------
61
62shared_ptr<MeshAxis> MeshAxis::clone() const {
63 //return plask::make_shared<MidpointAxis>(wrapped);
64 return plask::make_shared<OrderedAxis>(*this);
65}
66
67std::size_t MeshAxis::findIndex(double to_find) const {
68 return std::lower_bound(begin(), end(), to_find).index;
69}
70
71std::size_t MeshAxis::findUpIndex(double to_find) const {
72 return std::upper_bound(begin(), end(), to_find).index;
73}
74
75std::size_t MeshAxis::findNearestIndex(double to_find) const {
76 return find_nearest_binary(begin(), end(), to_find).index;
77}
78
79shared_ptr<MeshAxis> MeshAxis::getMidpointAxis() const {
80 beforeCalcMidpointMesh();
81 return plask::make_shared<MidpointAxis>(*this)->clone();
82}
83
84 void MeshAxis::beforeCalcMidpointMesh() const {
85 if (this->size() < 2)
86 throw BadMesh("getMidpointAxis", "at least two points are required");
87}
88
89 PLASK_API void prepareNearestNeighborInterpolationForAxis(const MeshAxis& axis, const InterpolationFlags& flags, double& wrapped_point_coord, int axis_nr) {
90 if (flags.periodic(axis_nr) && !flags.symmetric(axis_nr)) {
91 if (wrapped_point_coord < axis.at(0)) {
92 if (axis.at(0) - wrapped_point_coord > wrapped_point_coord - flags.low(axis_nr) + flags.high(axis_nr) - axis.at(axis.size()-1)) wrapped_point_coord = axis.at(axis.size()-1);
93 } else if (wrapped_point_coord > axis.at(axis.size()-1)) {
94 if (wrapped_point_coord - axis.at(axis.size()-1) > flags.high(axis_nr) - wrapped_point_coord + axis.at(0) - flags.low(axis_nr)) wrapped_point_coord = axis.at(0);
95 }
96 }
97}
98
99 PLASK_API void prepareInterpolationForAxis(const MeshAxis& axis, const InterpolationFlags& flags, double wrapped_point_coord, int axis_nr, std::size_t& index_lo, std::size_t& index_hi, double& lo, double& hi, bool& invert_lo, bool& invert_hi) {
100 index_hi = axis.findUpIndex(wrapped_point_coord);
101 invert_lo = false; invert_hi = false;
102 if (index_hi == 0) {
103 if (flags.symmetric(axis_nr)) {
104 index_lo = 0;
105 lo = axis.at(0);
106 if (lo > 0.) {
107 lo = - lo;
108 invert_lo = true;
109 } else if (flags.periodic(axis_nr)) {
110 lo = 2. * flags.low(axis_nr) - lo;
111 invert_lo = true;
112 } else {
113 lo -= 1.;
114 }
115 } else if (flags.periodic(axis_nr)) {
116 index_lo = axis.size() - 1;
117 lo = axis.at(index_lo) - flags.high(axis_nr) + flags.low(axis_nr);
118 } else {
119 index_lo = 0;
120 lo = axis.at(0) - 1.;
121 }
122 } else {
123 index_lo = index_hi - 1;
124 lo = axis.at(index_lo);
125 }
126 if (index_hi == axis.size()) {
127 if (flags.symmetric(axis_nr)) {
128 --index_hi;
129 hi = axis.at(index_hi);
130 if (hi < 0.) {
131 hi = - hi;
132 invert_hi = true;
133 } else if (flags.periodic(axis_nr)) {
134 if (wrapped_point_coord == flags.high(axis_nr)) {
135 index_lo = index_hi - 1;
136 lo = lo = axis.at(index_lo);
137 } else {
138 lo = 2. * flags.high(axis_nr) - hi;
139 invert_hi = true;
140 }
141 } else {
142 hi += 1.;
143 }
144 } else if (flags.periodic(axis_nr)) {
145 index_hi = 0;
146 hi = axis.at(0) + flags.high(axis_nr) - flags.low(axis_nr);
147 if (hi == lo) hi += 1e-6;
148 } else {
149 --index_hi;
150 hi = axis.at(index_hi) + 1.;
151 }
152 } else {
153 hi = axis.at(index_hi);
154 }
155}
156
157} // namespace plask