PLaSK library

Loading...
Searching...
No Matches
prism.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 "prism.hpp"
15#include "reader.hpp"
16
17#include "../manager.hpp"
18
19#define PLASK_TRIANGULAR_PRISM_NAME "triangular-prism"
20#define PLASK_PRISM_NAME "prism"
21
22 namespace plask {
23
24 const char* TriangularPrism::NAME = PLASK_TRIANGULAR_PRISM_NAME;
25
26std::string TriangularPrism::getTypeName() const { return NAME; }
27
28 TriangularPrism::TriangularPrism(const TriangularPrism::Vec2& p0, const TriangularPrism::Vec2& p1, double height, const shared_ptr<Material>& material)
29 : BaseClass(material), p0(p0), p1(p1), height(height) {}
30
31 TriangularPrism::TriangularPrism(const TriangularPrism::Vec2& p0,
32 const TriangularPrism::Vec2& p1,
33 double height,
34 shared_ptr<MaterialsDB::MixedCompositionFactory> materialTopBottom)
35 : BaseClass(materialTopBottom), p0(p0), p1(p1), height(height) {}
36
37 Box3D TriangularPrism::getBoundingBox() const {
38 return Box3D(std::min(std::min(p0.c0, p1.c0), 0.0), std::min(std::min(p0.c1, p1.c1), 0.0), 0.,
39 std::max(std::max(p0.c0, p1.c0), 0.0), std::max(std::max(p0.c1, p1.c1), 0.0), height);
40}
41
42 inline static double sign(const Vec<3, double>& p1, const LateralVec<double>& p2, const LateralVec<double>& p3) {
43 return (p1.c0 - p3.c0) * (p2.c1 - p3.c1) - (p2.c0 - p3.c0) * (p1.c1 - p3.c1);
44}
45
46// Like sign, but with p3 = (0, 0)
47 inline static double sign0(const Vec<3, double>& p1, const LateralVec<double>& p2) {
48 return (p1.c0) * (p2.c1) - (p2.c0) * (p1.c1);
49}
50
51 bool TriangularPrism::contains(const TriangularPrism::DVec& p) const {
52 if (p.c2 < 0 || p.c2 > height) return false;
53 // algorithm comes from:
54 // http://stackoverflow.com/questions/2049582/how-to-determine-a-point-in-a-triangle
55 // with: v1 -> p0, v2 -> p1, v3 -> (0, 0)
56 // maybe barycentric method would be better?
57 bool b1 = sign(p, p0, p1) < 0.0;
58 bool b2 = sign0(p, p1) < 0.0;
59 return (b1 == b2) && (b2 == (sign(p, Primitive<2>::ZERO_VEC, p0) < 0.0));
60}
61
62 void TriangularPrism::writeXMLAttr(XMLWriter::Element& dest_xml_object, const AxisNames& axes) const {
63 BaseClass::writeXMLAttr(dest_xml_object, axes);
64 materialProvider->writeXML(dest_xml_object, axes)
65 .attr("a" + axes.getNameForLong(), p0.tran())
66 .attr("a" + axes.getNameForTran(), p0.vert())
67 .attr("b" + axes.getNameForLong(), p1.tran())
68 .attr("b" + axes.getNameForTran(), p1.vert())
69 .attr("height", height);
70}
71
72 void TriangularPrism::addPointsAlongToSet(std::set<double>& points,
73 Primitive<3>::Direction direction,
74 unsigned max_steps,
75 double min_step_size) const {
76 if (direction == Primitive<3>::DIRECTION_VERT) {
77 if (materialProvider->isUniform(Primitive<3>::DIRECTION_VERT)) {
78 points.insert(0);
79 points.insert(height);
80 } else {
81 if (this->max_steps) max_steps = this->max_steps;
82 if (this->min_step_size) min_step_size = this->min_step_size;
83 unsigned steps = min(unsigned(height / min_step_size), max_steps);
84 double step = height / steps;
85 for (unsigned i = 0; i <= steps; ++i) points.insert(i * step);
86 }
87 } else {
88 if (this->max_steps) max_steps = this->max_steps;
89 if (this->min_step_size) min_step_size = this->min_step_size;
90
91 double x[3] = {0., p0[int(direction)], p1[int(direction)]};
92
93 // Sort x
94 if (x[2] < x[0]) std::swap(x[0], x[2]);
95 if (x[1] > x[2])
96 std::swap(x[1], x[2]);
97 else if (x[1] < x[0])
98 std::swap(x[1], x[0]);
99
100 for (int i = 0; i < 3; ++i) points.insert(x[i]);
101 double dx02 = x[2] - x[0];
102 if (dx02 == 0) return;
103
104 for (int i = 0; i < 2; ++i) {
105 double dx = x[i + 1] - x[i];
106 PLASK_NO_CONVERSION_WARNING_BEGIN
107 unsigned maxsteps = max_steps * (dx / dx02);
108 PLASK_NO_WARNING_END
109 unsigned steps = min(unsigned(dx / min_step_size), maxsteps);
110 double step = dx / steps;
111 for (unsigned j = 1; j < steps; ++j) points.insert(x[i] + j * step);
112 }
113 }
114}
115
116 void TriangularPrism::addLineSegmentsToSet(std::set<typename GeometryObjectD<3>::LineSegment>& segments,
117 unsigned max_steps,
118 double min_step_size) const {
119 if (!materialProvider->isUniform(Primitive<3>::DIRECTION_LONG))
120 throw NotImplemented("prismatic mesh for prisms non-uniform in longitudinal direction");
121 if (!materialProvider->isUniform(Primitive<3>::DIRECTION_TRAN))
122 throw NotImplemented("prismatic mesh for prisms non-uniform in transverse direction");
123 std::set<double> vert;
124 addPointsAlongToSet(vert, Primitive<3>::DIRECTION_VERT, max_steps, min_step_size);
125 typedef typename GeometryObjectD<3>::LineSegment Segment;
126 double pv = 0.;
127 for (double v : vert) {
128 segments.insert(Segment(DVec(0., 0., v), DVec(p0[0], p0[1], v)));
129 segments.insert(Segment(DVec(0., 0., v), DVec(p1[0], p1[1], v)));
130 segments.insert(Segment(DVec(p0[0], p0[1], v), DVec(p1[0], p1[1], v)));
131 if (v != 0.) {
132 segments.insert(Segment(DVec(0., 0., pv), DVec(0., 0., v)));
133 segments.insert(Segment(DVec(p0[0], p0[1], pv), DVec(p0[0], p0[1], v)));
134 segments.insert(Segment(DVec(p1[0], p1[1], pv), DVec(p1[0], p1[1], v)));
135 }
136 pv = v;
137 }
138}
139
140 shared_ptr<GeometryObject> read_triangular_prism(GeometryReader& reader) {
141 shared_ptr<TriangularPrism> prism(new TriangularPrism());
142 if (reader.manager.draft) {
143 prism->p0.c0 = reader.source.getAttribute("a" + reader.getAxisLongName(), 0.0);
144 prism->p0.c1 = reader.source.getAttribute("a" + reader.getAxisTranName(), 0.0);
145 prism->p1.c0 = reader.source.getAttribute("b" + reader.getAxisLongName(), 0.0);
146 prism->p1.c1 = reader.source.getAttribute("b" + reader.getAxisTranName(), 0.0);
147 prism->height = reader.source.getAttribute("height", 0.0);
148 } else {
149 prism->p0.c0 = reader.source.requireAttribute<double>("a" + reader.getAxisLongName());
150 prism->p0.c1 = reader.source.requireAttribute<double>("a" + reader.getAxisTranName());
151 prism->p1.c0 = reader.source.requireAttribute<double>("b" + reader.getAxisLongName());
152 prism->p1.c1 = reader.source.requireAttribute<double>("b" + reader.getAxisTranName());
153 prism->height = reader.source.requireAttribute<double>("height");
154 }
155 prism->readMaterial(reader);
156 reader.source.requireTagEnd();
157 return prism;
158}
159
160 static GeometryReader::RegisterObjectReader prism_reader(PLASK_TRIANGULAR_PRISM_NAME, read_triangular_prism);
161
163
164 const char* Prism::NAME = PLASK_PRISM_NAME;
165
166std::string Prism::getTypeName() const { return NAME; }
167
168 Prism::Prism(double height, const std::vector<LateralVec<double>>& vertices, const shared_ptr<Material>& material)
169 : BaseClass(material), height(height), vertices(vertices) {}
170
171 Prism::Prism(double height, std::vector<LateralVec<double>>&& vertices, const shared_ptr<Material>&& material)
172 : BaseClass(material), height(height), vertices(std::move(vertices)) {}
173
174 Prism::Prism(double height, std::initializer_list<LateralVec<double>> vertices, const shared_ptr<Material>& material)
175 : BaseClass(material), height(height), vertices(vertices) {}
176
177 Prism::Prism(double height,
178 const std::vector<LateralVec<double>>& vertices,
179 shared_ptr<MaterialsDB::MixedCompositionFactory> materialTopBottom)
180 : BaseClass(materialTopBottom), height(height), vertices(vertices) {}
181
182 Prism::Prism(double height,
183 std::vector<LateralVec<double>>&& vertices,
184 shared_ptr<MaterialsDB::MixedCompositionFactory> materialTopBottom)
185 : BaseClass(materialTopBottom), height(height), vertices(std::move(vertices)) {}
186
187 Prism::Prism(double height,
188 std::initializer_list<LateralVec<double>> vertices,
189 shared_ptr<MaterialsDB::MixedCompositionFactory> materialTopBottom)
190 : BaseClass(materialTopBottom), height(height), vertices(vertices) {}
191
192 void Prism::validate() const {
193 if (vertices.size() < 3) {
194 throw GeometryException("polygon has less than 3 vertices");
195 }
196 // if (!checkSegments()) {
197 // throw GeometryException("polygon has intersecting segments");
198 // }
199}
200
201// bool Prism::checkSegments() const {
202// for (size_t i = 0; i < vertices.size(); ++i) {
203// LateralVec<double> a1 = vertices[i];
204// LateralVec<double> b1 = vertices[(i + 1) % vertices.size()];
205// double min1x = std::min(a1.c0, b1.c0), max1x = std::max(a1.c0, b1.c0);
206// double min1y = std::min(a1.c1, b1.c1), max1y = std::max(a1.c1, b1.c1);
207// double d1x = b1.c0 - a1.c0;
208// double d1y = b1.c1 - a1.c1;
209// double det1 = a1.c0 * b1.c1 - a1.c1 * b1.c0;
210// for (size_t j = i + 2; j < vertices.size() - (i ? 0 : 1); ++j) {
211// LateralVec<double> a2 = vertices[j];
212// LateralVec<double> b2 = vertices[(j + 1) % vertices.size()];
213// double min2x = std::min(a2.c0, b2.c0), max2x = std::max(a2.c0, b2.c0);
214// double min2y = std::min(a2.c1, b2.c1), max2y = std::max(a2.c1, b2.c1);
215// if (max2x < min1x || max1x < min2x || max2y < min1y || max1y < min2y) continue;
216// double d2x = b2.c0 - a2.c0;
217// double d2y = b2.c1 - a2.c1;
218// double det = d1x * d2y - d2x * d1y;
219// double det2 = a2.c0 * b2.c1 - a2.c1 * b2.c0;
220// if (det == 0) continue;
221// double x = (d1x * det2 - d2x * det1) / det;
222// double y = (d1y * det2 - d2y * det1) / det;
223// if (x >= min1x && x <= max1x && x >= min2x && x <= max2x && y >= min1y && y <= max1y && y >= min2y && y <= max2y)
224// return false;
225// }
226// }
227// return true;
228// }
229
230 Prism::Box Prism::getBoundingBox() const {
231 if (vertices.empty()) return Box(DVec(0, 0, 0), DVec(0, 0, 0));
232 double min_x = vertices[0].c0;
233 double max_x = vertices[0].c0;
234 double min_y = vertices[0].c1;
235 double max_y = vertices[0].c1;
236 for (const LateralVec<double>& v : vertices) {
237 min_x = std::min(min_x, v.c0);
238 max_x = std::max(max_x, v.c0);
239 min_y = std::min(min_y, v.c1);
240 max_y = std::max(max_y, v.c1);
241 }
242 return Box(DVec(min_x, min_y, 0), DVec(max_x, max_y, height));
243}
244
245 bool Prism::contains(const DVec& p) const {
246 if (vertices.size() < 3) return false;
247 if (p.c2 < 0 || p.c2 > height) return false;
248 int n = vertices.size();
249 int i, j;
250 int c = 0;
251 for (i = 0, j = n - 1; i < n; j = i++) {
252 if (((vertices[i].c1 > p.c1) != (vertices[j].c1 > p.c1)) &&
253 (p.c0 <
254 (vertices[j].c0 - vertices[i].c0) * (p.c1 - vertices[i].c1) / (vertices[j].c1 - vertices[i].c1) + vertices[i].c0))
255 c += (vertices[i].c1 > vertices[j].c1) ? 1 : -1;
256 }
257 return c;
258}
259
260 void Prism::addPointsAlongToSet(std::set<double>& points,
261 Primitive<3>::Direction direction,
262 unsigned max_steps,
263 double min_step_size) const {
264 if (direction == Primitive<3>::Direction::DIRECTION_VERT) {
265 if (materialProvider->isUniform(Primitive<3>::DIRECTION_VERT)) {
266 points.insert(0);
267 points.insert(height);
268 } else {
269 if (this->max_steps) max_steps = this->max_steps;
270 if (this->min_step_size) min_step_size = this->min_step_size;
271 unsigned steps = min(unsigned(height / min_step_size), max_steps);
272 double step = height / steps;
273 for (unsigned i = 0; i <= steps; ++i) points.insert(i * step);
274 }
275 return;
276 }
277
278 if (vertices.size() < 3) return;
279 std::set<double> vert;
280 for (const LateralVec<double>& v : vertices) {
281 vert.insert(v[int(direction)]);
282 }
283 for (std::set<double>::const_iterator b = vert.begin(), a = b++; b != vert.end(); ++a, ++b) {
284 double d = *b - *a;
285 unsigned steps = std::max(1u, static_cast<unsigned>(d / min_step_size));
286 steps = std::min(steps, max_steps);
287 double step = d / steps;
288 for (unsigned i = 0; i <= steps; ++i) points.insert(*a + i * step);
289 }
290}
291
292 void Prism::addLineSegmentsToSet(std::set<typename GeometryObjectD<3>::LineSegment>& segments,
293 unsigned max_steps,
294 double min_step_size) const {
295 if (vertices.size() < 3) return;
296 typedef typename GeometryObjectD<3>::LineSegment Segment;
297 for (size_t i = 0; i < vertices.size(); ++i) {
298 LateralVec<double> a = vertices[i];
299 LateralVec<double> b = vertices[(i + 1) % vertices.size()];
300 LateralVec<double> ab = b - a;
301 double d = std::sqrt(dot(ab, ab));
302 unsigned steps = std::max(1u, static_cast<unsigned>(d / min_step_size));
303 steps = std::min(steps, max_steps);
304 LateralVec<double> p0 = a;
305 for (unsigned j = 1; j <= steps; ++j) {
306 segments.insert(Segment(DVec(p0.c0, p0.c1, 0), DVec(p0.c0, p0.c1, height)));
307 double t = static_cast<double>(j) / steps;
308 LateralVec<double> p = a * (1 - t) + b * t;
309 segments.insert(Segment(DVec(p0.c0, p0.c1, 0), DVec(p.c0, p.c1, 0)));
310 segments.insert(Segment(DVec(p0.c0, p0.c1, height), DVec(p.c0, p.c1, height)));
311 }
312 }
313}
314
315 void Prism::writeXMLAttr(XMLWriter::Element& dest_xml_object, const AxisNames& axes) const {
316 BaseClass::writeXMLAttr(dest_xml_object, axes);
317 materialProvider->writeXML(dest_xml_object, axes).attr("height", height);
318 if (vertices.empty()) return;
319 std::string vertices_str;
320 const char* sep = "";
321 for (const LateralVec<double>& v : vertices) {
322 vertices_str += sep;
323 vertices_str += str(v.c0) + " " + str(v.c1);
324 sep = "; ";
325 }
326 dest_xml_object.writeText(vertices_str);
327}
328
329 shared_ptr<GeometryObject> read_prism(GeometryReader& reader) {
330
331 // DEPRECATED: use TriangularPrism instead
332 if (reader.source.hasAttribute("a" + reader.getAxisLongName()) ||
333 reader.source.hasAttribute("a" + reader.getAxisTranName()) ||
334 reader.source.hasAttribute("b" + reader.getAxisLongName()) ||
335 reader.source.hasAttribute("b" + reader.getAxisTranName())) {
336 writelog(LOG_WARNING, "<prism> with vertices a and b is deprecated, use <triangular-prism> instead");
337 return read_triangular_prism(reader);
338 }
339
340 shared_ptr<Prism> prism = make_shared<Prism>();
341 prism->readMaterial(reader);
342 if (reader.manager.draft)
343 prism->height = reader.source.getAttribute("height", 0.0);
344 else
345 prism->height = reader.source.requireAttribute<double>("height");
346
347 std::string vertex_spec = reader.source.requireTextInCurrentTag();
348 if (reader.source.attributeFilter) vertex_spec = reader.source.attributeFilter(vertex_spec);
349 std::vector<LateralVec<double>> vertices;
350 boost::tokenizer<boost::char_separator<char>> tokens(vertex_spec, boost::char_separator<char>(" \t\n\r", ";"));
351 int vi = 0;
352 for (const std::string& t : tokens) {
353 if (t == ";") { // end of point or segment
354 if (vi != 2) throw Exception("each vertex must have two coordinates");
355 vi = 0;
356 } else { // end of point coordinate
357 if (vi == 2) throw Exception("end of vertex (\";\") was expected, but got \"{0}\"", t);
358 if (vi == 0) vertices.emplace_back();
359 try {
360 vertices.back()[vi++] = boost::lexical_cast<double>(t);
361 } catch (const boost::bad_lexical_cast&) {
362 throw Exception("bad vertex coordinate: {0}", t);
363 }
364 }
365 }
366 prism->vertices = std::move(vertices);
367 if (!reader.manager.draft) prism->validate();
368 return prism;
369}
370
371 static GeometryReader::RegisterObjectReader polygon_reader(PLASK_PRISM_NAME, read_prism);
372
373} // namespace plask