PLaSK library
Loading...
Searching...
No Matches
equilateral3d.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__MESH_EQUILATERAL3D_H
15#define PLASK__MESH_EQUILATERAL3D_H
16
21#include "rectilinear3d.hpp"
22
23namespace plask {
24
25class PLASK_API EquilateralMesh3D: public RectilinearMesh3D {
26
27 protected:
28 double trans[9];
29
30 private:
31 double inv[9];
32
34 void findInverse() {
35 double idet = 1. / (trans[0] * (trans[4]*trans[8] - trans[5]*trans[7]) +
36 trans[1] * (trans[5]*trans[6] - trans[3]*trans[8]) +
37 trans[2] * (trans[3]*trans[7] - trans[4]*trans[6]));
38 inv[0] = idet * ( trans[4]*trans[8] - trans[5]*trans[7]);
39 inv[1] = idet * (-trans[1]*trans[8] + trans[2]*trans[7]);
40 inv[2] = idet * ( trans[1]*trans[5] - trans[2]*trans[4]);
41 inv[3] = idet * (-trans[3]*trans[8] + trans[5]*trans[6]);
42 inv[4] = idet * ( trans[0]*trans[8] - trans[2]*trans[6]);
43 inv[5] = idet * (-trans[0]*trans[5] + trans[2]*trans[3]);
44 inv[6] = idet * ( trans[3]*trans[7] - trans[4]*trans[6]);
45 inv[7] = idet * (-trans[0]*trans[7] + trans[1]*trans[6]);
46 inv[8] = idet * ( trans[0]*trans[4] - trans[1]*trans[3]);
47 }
48
49 public:
50
52 struct Transformed: public MeshD<3> {
53 const EquilateralMesh3D* src;
54 shared_ptr<const MeshD<3>> dst;
55 Transformed(const EquilateralMesh3D* src, const shared_ptr<const MeshD<3>>& dst):
56 src(src), dst(dst) {}
57 Transformed(const shared_ptr<const EquilateralMesh3D>& src, const shared_ptr<const MeshD<3>>& dst):
58 src(src.get()), dst(dst) {}
59 size_t size() const override { return dst->size(); }
60 Vec<3,double> at(size_t index) const override {
61 return src->toMeshCoords(dst->at(index));
62 }
63 };
64
65 typedef RectilinearMesh3D::ElementMesh<EquilateralMesh3D> ElementMesh;
66
74 explicit EquilateralMesh3D(IterationOrder iterationOrder = ORDER_012,
75 Vec<3> vec0 = vec(1., 0., 0.), Vec<3> vec1 = vec(0., 1., 0.), Vec<3> vec2 = vec(0., 0., 1.));
76
88 EquilateralMesh3D(shared_ptr<MeshAxis> mesh0, shared_ptr<MeshAxis> mesh1, shared_ptr<MeshAxis> mesh2, IterationOrder iterationOrder = ORDER_012,
89 Vec<3> vec0 = vec(1., 0., 0.), Vec<3> vec1 = vec(0., 1., 0.), Vec<3> vec2 = vec(0., 0., 1.));
90
92 Vec<3,double> getVec0() const {
93 return vec<double>(trans[0], trans[3], trans[6]);
94 }
95
97 Vec<3,double> getVec1() const {
98 return vec<double>(trans[1], trans[4], trans[7]);
99 }
100
102 Vec<3,double> getVec2() const {
103 return vec<double>(trans[2], trans[5], trans[8]);
104 }
105
107 Vec<3,double> getInvVec0() const {
108 return vec<double>(inv[0], inv[1], inv[2]);
109 }
110
112 Vec<3,double> getInvVec1() const {
113 return vec<double>(inv[3], inv[4], inv[5]);
114 }
115
117 Vec<3,double> getInvVec2() const {
118 return vec<double>(inv[6], inv[7], inv[8]);
119 }
120
122 void setVec0(Vec<3,double> vec0) {
123 trans[0] = vec0.c0;
124 trans[3] = vec0.c1;
125 trans[6] = vec0.c2;
126 findInverse();
127 fireChanged(Event::EVENT_USER_DEFINED);
128 }
129
131 void setVec1(Vec<3,double> vec1) {
132 trans[1] = vec1.c0;
133 trans[4] = vec1.c1;
134 trans[7] = vec1.c2;
135 findInverse();
136 fireChanged(Event::EVENT_USER_DEFINED);
137 }
138
140 void setVec2(Vec<3,double> vec2) {
141 trans[2] = vec2.c0;
142 trans[5] = vec2.c1;
143 trans[8] = vec2.c2;
144 findInverse();
145 fireChanged();
146 }
147
153 Vec<3,double> toMeshCoords(Vec<3,double> point) const {
154 return Vec<3,double>(inv[0] * point.c0 + inv[1] * point.c1 + inv[2] * point.c2,
155 inv[3] * point.c0 + inv[4] * point.c1 + inv[5] * point.c2,
156 inv[6] * point.c0 + inv[7] * point.c1 + inv[8] * point.c2);
157 }
158
164 Vec<3,double> toMeshCoords(double c0, double c1, double c2) const {
165 return Vec<3,double>(inv[0] * c0 + inv[1] * c1 + inv[2] * c2,
166 inv[3] * c0 + inv[4] * c1 + inv[5] * c2,
167 inv[6] * c0 + inv[7] * c1 + inv[8] * c2);
168 }
169
175 inline Vec<3,double> fromMeshCoords(Vec<3,double> coords) const {
176 return Vec<3,double>(trans[0] * coords.c0 + trans[1] * coords.c1 + trans[2] * coords.c2,
177 trans[3] * coords.c0 + trans[4] * coords.c1 + trans[5] * coords.c2,
178 trans[6] * coords.c0 + trans[7] * coords.c1 + trans[8] * coords.c2);
179 }
180
186 inline Vec<3,double> fromMeshCoords(double c0, double c1, double c2) const {
187 return Vec<3,double>(trans[0] * c0 + trans[1] * c1 + trans[2] * c2,
188 trans[3] * c0 + trans[4] * c1 + trans[5] * c2,
189 trans[6] * c0 + trans[7] * c1 + trans[8] * c2);
190 }
191
199 Vec<3, double> at(std::size_t index0, std::size_t index1, std::size_t index2) const override {
200 return fromMeshCoords(axis[0]->at(index0), axis[1]->at(index1), axis[2]->at(index2));
201 }
202
207 shared_ptr<EquilateralMesh3D::ElementMesh> getElementMesh() const;
208
214 Vec<3,double> getElementMidpoint(std::size_t index0, std::size_t index1, std::size_t index2) const override {
215 return fromMeshCoords(getElementMidpoint0(index0), getElementMidpoint1(index1), getElementMidpoint2(index2));
216 }
217};
218
219
220template <typename SrcT, typename DstT, InterpolationMethod method>
221struct InterpolationAlgorithm<typename std::enable_if<method != INTERPOLATION_DEFAULT, EquilateralMesh3D>::type, SrcT, DstT, method> {
222 static LazyData<DstT> interpolate(const shared_ptr<const EquilateralMesh3D>& src_mesh, const DataVector<const SrcT>& src_vec,
223 const shared_ptr<const MeshD<3>>& dst_mesh, const InterpolationFlags& flags) {
224 //TODO Proprably we need to prohibit symmetry (or warn that its meaning is different)
225 return InterpolationAlgorithm<RectilinearMesh3D, SrcT, DstT, method>::interpolate(src_mesh, src_vec, EquilateralMesh3D::Transformed(src_mesh, dst_mesh), flags);
226 }
227};
228
229template <typename SrcT, typename DstT, InterpolationMethod method>
230struct InterpolationAlgorithm<typename std::enable_if<method != INTERPOLATION_DEFAULT, EquilateralMesh3D::ElementMesh>::type, SrcT, DstT, method> {
231 static LazyData<DstT> interpolate(const shared_ptr<const EquilateralMesh3D::ElementMesh>& src_mesh, const DataVector<const SrcT>& src_vec,
232 const shared_ptr<const MeshD<3>>& dst_mesh, const InterpolationFlags& flags) {
233 //TODO Proprably we need to prohibit symmetry (or warn that its meaning is different)
234 return InterpolationAlgorithm<RectilinearMesh3D, SrcT, DstT, method>::interpolate(src_mesh, src_vec, EquilateralMesh3D::Transformed(src_mesh, dst_mesh), flags);
235 }
236};
237
238template <typename SrcT, typename DstT>
239struct InterpolationAlgorithm<EquilateralMesh3D::ElementMesh, SrcT, DstT, INTERPOLATION_NEAREST> {
240 static LazyData<DstT> interpolate(const shared_ptr<const EquilateralMesh3D::ElementMesh>& src_mesh, const DataVector<const SrcT>& src_vec,
241 const shared_ptr<const MeshD<3>>& dst_mesh, const InterpolationFlags& flags) {
242 if (src_mesh->axis[0]->size() == 0 || src_mesh->axis[1]->size() == 0)
243 throw BadMesh("interpolate", "source mesh empty");
244 //TODO Proprably we need to prohibit symmetry (or warn that its meaning is different)
245 return new NearestNeighborInterpolatedLazyDataImpl<DstT, EquilateralMesh3D::ElementMesh, SrcT>(src_mesh, src_vec, EquilateralMesh3D::Transformed(src_mesh, dst_mesh), flags);
246 }
247};
248
249} // namespace plask
250
251#endif // PLASK__SMESH_EQUILATERAL3D_H