PLaSK library
Loading...
Searching...
No Matches
triangle.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 "triangle.hpp"
15#include "reader.hpp"
16
17#include "../manager.hpp"
18
19#define PLASK_TRIANGLE_NAME "triangle"
20
21namespace plask {
22
23const char* Triangle::NAME = PLASK_TRIANGLE_NAME;
24
25std::string Triangle::getTypeName() const { return NAME; }
26
27Triangle::Triangle(const Triangle::DVec& p0, const Triangle::DVec& p1, const shared_ptr<Material>& material)
28 : BaseClass(material), p0(p0), p1(p1) {}
29
34
35Box2D Triangle::getBoundingBox() const {
36 return Box2D(std::min(std::min(p0.c0, p1.c0), 0.0), std::min(std::min(p0.c1, p1.c1), 0.0),
37 std::max(std::max(p0.c0, p1.c0), 0.0), std::max(std::max(p0.c1, p1.c1), 0.0));
38}
39
40inline static double sign(const Vec<2, double>& p1, const Vec<2, double>& p2, const Vec<2, double>& p3) {
41 return (p1.c0 - p3.c0) * (p2.c1 - p3.c1) - (p2.c0 - p3.c0) * (p1.c1 - p3.c1);
42}
43
44// Like sign, but with p3 = (0, 0)
45inline static double sign0(const Vec<2, double>& p1, const Vec<2, double>& p2) {
46 return (p1.c0) * (p2.c1) - (p2.c0) * (p1.c1);
47}
48
49bool Triangle::contains(const Triangle::DVec& p) const {
50 // algorithm comes from:
51 // http://stackoverflow.com/questions/2049582/how-to-determine-a-point-in-a-triangle
52 // with: v1 -> p0, v2 -> p1, v3 -> (0, 0)
53 // maybe barycentric method would be better?
54 bool b1 = sign(p, p0, p1) < 0.0;
55 bool b2 = sign0(p, p1) < 0.0;
56 return (b1 == b2) && (b2 == (sign(p, Primitive<2>::ZERO_VEC, p0) < 0.0));
57}
58
59void Triangle::addPointsAlongToSet(std::set<double>& points,
60 Primitive<3>::Direction direction,
61 unsigned max_steps,
62 double min_step_size) const {
63 assert(0 < int(direction) && int(direction) < 3);
64 if (this->max_steps) max_steps = this->max_steps;
65 if (this->min_step_size) min_step_size = this->min_step_size;
66
67 double x[3] = {0., p0[int(direction) - 1], p1[int(direction) - 1]};
68
69 // Sort x
70 if (x[2] < x[0]) std::swap(x[0], x[2]);
71 if (x[1] > x[2])
72 std::swap(x[1], x[2]);
73 else if (x[1] < x[0])
74 std::swap(x[1], x[0]);
75
76 for (int i = 0; i < 3; ++i) points.insert(x[i]);
77 double dx02 = x[2] - x[0];
78 if (dx02 == 0) return;
79
80 for (int i = 0; i < 2; ++i) {
81 double dx = x[i + 1] - x[i];
82PLASK_NO_CONVERSION_WARNING_BEGIN
83 unsigned maxsteps = max_steps * (dx / dx02);
84PLASK_NO_WARNING_END
85 unsigned steps = min(unsigned(dx / min_step_size), maxsteps);
86 double step = dx / steps;
87 for (unsigned j = 1; j < steps; ++j) points.insert(x[i] + j * step);
88 }
89}
90
91void Triangle::addLineSegmentsToSet(std::set<typename GeometryObjectD<2>::LineSegment>& segments,
92 unsigned max_steps,
93 double min_step_size) const {
94 if (!materialProvider->isUniform(Primitive<3>::DIRECTION_TRAN))
95 throw NotImplemented("triangular mesh for triangles non-uniform in transverse direction");
96 typedef typename GeometryObjectD<2>::LineSegment Segment;
97 if (materialProvider->isUniform(Primitive<3>::DIRECTION_VERT)) {
98 segments.insert(Segment(Primitive<2>::ZERO_VEC, p0));
99 segments.insert(Segment(Primitive<2>::ZERO_VEC, p1));
100 segments.insert(Segment(p0, p1));
101 } else {
102 if (this->max_steps) max_steps = this->max_steps;
103 if (this->min_step_size) min_step_size = this->min_step_size;
104
105 // Here we replace x and y for simplicity of analysis
106 double x[3] = {0., p0[1], p1[1]};
107 double y[3] = {0., p0[0], p1[0]};
108
109 // Sort x
110 if (x[2] < x[0]) {
111 std::swap(x[0], x[2]);
112 std::swap(y[0], y[2]);
113 }
114 if (x[1] > x[2]) {
115 std::swap(x[1], x[2]);
116 std::swap(y[1], y[2]);
117 } else if (x[1] < x[0]) {
118 std::swap(x[1], x[0]);
119 std::swap(y[1], y[0]);
120 }
121
122 double dx02 = x[2] - x[0];
123 if (dx02 == 0) return;
124
125 double a2 = (y[2] - y[0]) / dx02, b2 = (x[2] * y[0] - x[0] * y[2]) / dx02;
126
127 DVec d1(y[0], x[0]), d2(y[0], x[0]);
128 for (int i = 0; i < 2; ++i) {
129 double dx = x[i + 1] - x[i];
130 double a1 = (y[i + 1] - y[i]) / dx, b1 = (x[i + 1] * y[i] - x[i] * y[i + 1]) / dx;
131PLASK_NO_CONVERSION_WARNING_BEGIN
132 unsigned maxsteps = max_steps * (dx / dx02);
133PLASK_NO_WARNING_END
134 unsigned steps = min(unsigned(dx / min_step_size), maxsteps);
135 if (steps < 2) continue;
136 double step = dx / steps;
137 for (unsigned j = 0; j < steps; ++j) {
138 double t = x[i] + j * step;
139 DVec e1(a1 * t + b1, t), e2(a2 * t + b2, t);
140 if (i != 0 || j != 0) {
141 segments.insert(Segment(d1, e1));
142 segments.insert(Segment(d2, e2));
143 segments.insert(Segment(e1, e2));
144 }
145 d1 = e1;
146 d2 = e2;
147 }
148 }
149 DVec e(y[2], x[2]);
150 segments.insert(Segment(d1, e));
151 segments.insert(Segment(d2, e));
152 }
153}
154
155void Triangle::writeXMLAttr(XMLWriter::Element& dest_xml_object, const AxisNames& axes) const {
156 BaseClass::writeXMLAttr(dest_xml_object, axes);
157 materialProvider->writeXML(dest_xml_object, axes)
158 .attr("a" + axes.getNameForTran(), p0.tran())
159 .attr("a" + axes.getNameForVert(), p0.vert())
160 .attr("b" + axes.getNameForTran(), p1.tran())
161 .attr("b" + axes.getNameForVert(), p1.vert());
162}
163
164shared_ptr<GeometryObject> read_triangle(GeometryReader& reader) {
165 shared_ptr<Triangle> triangle(new Triangle());
166 if (reader.manager.draft) {
167 triangle->p0.tran() = reader.source.getAttribute("a" + reader.getAxisTranName(), 0.0);
168 triangle->p0.vert() = reader.source.getAttribute("a" + reader.getAxisVertName(), 0.0);
169 triangle->p1.tran() = reader.source.getAttribute("b" + reader.getAxisTranName(), 0.0);
170 triangle->p1.vert() = reader.source.getAttribute("b" + reader.getAxisVertName(), 0.0);
171 } else {
172 triangle->p0.tran() = reader.source.requireAttribute<double>("a" + reader.getAxisTranName());
173 triangle->p0.vert() = reader.source.requireAttribute<double>("a" + reader.getAxisVertName());
174 triangle->p1.tran() = reader.source.requireAttribute<double>("b" + reader.getAxisTranName());
175 triangle->p1.vert() = reader.source.requireAttribute<double>("b" + reader.getAxisVertName());
176 }
177 triangle->readMaterial(reader);
178 reader.source.requireTagEnd();
179 return triangle;
180}
181
182static GeometryReader::RegisterObjectReader triangle_reader(PLASK_TRIANGLE_NAME, read_triangle);
183
184} // namespace plask