PLaSK library
Loading...
Searching...
No Matches
elliptic_cylinder.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 "elliptic_cylinder.hpp"
15#include "../manager.hpp"
16#include "reader.hpp"
17
18#define PLASK_ELLIPTIC_CYLINDER_NAME "elliptic-cylinder"
19
20namespace plask {
21
23
24const char* EllipticCylinder::NAME = PLASK_ELLIPTIC_CYLINDER_NAME;
25
26EllipticCylinder::EllipticCylinder(double radius0,
27 double radius1,
28 double angle,
29 double height,
30 const shared_ptr<Material>& material)
31 : GeometryObjectLeaf<3>(material),
32 radius0(std::max(radius0, 0.)),
33 radius1(std::max(radius1, 0.)),
34 sina(sin(angle)),
35 cosa(cos(angle)),
36 height(std::max(height, 0.)) {}
37
38EllipticCylinder::EllipticCylinder(double radius0,
39 double radius1,
40 double angle,
41 double height,
42 const shared_ptr<MaterialsDB::MixedCompositionFactory>& materialTopBottom)
43 : GeometryObjectLeaf<3>(materialTopBottom),
44 radius0(std::max(radius0, 0.)),
45 radius1(std::max(radius1, 0.)),
46 sina(sin(angle)),
47 cosa(cos(angle)),
48 height(std::max(height, 0.)) {}
49
50EllipticCylinder::EllipticCylinder(double radius0, double radius1, double height, const shared_ptr<Material>& material)
51 : GeometryObjectLeaf<3>(material),
52 radius0(std::max(radius0, 0.)),
53 radius1(std::max(radius1, 0.)),
54 sina(0.),
55 cosa(1.),
56 height(std::max(height, 0.)) {}
57
58EllipticCylinder::EllipticCylinder(double radius0,
59 double radius1,
60 double height,
61 const shared_ptr<MaterialsDB::MixedCompositionFactory>& materialTopBottom)
62 : GeometryObjectLeaf<3>(materialTopBottom),
63 radius0(std::max(radius0, 0.)),
64 radius1(std::max(radius1, 0.)),
65 sina(0.),
66 cosa(1.),
67 height(std::max(height, 0.)) {}
68
69EllipticCylinder::EllipticCylinder(const EllipticCylinder& src)
70 : GeometryObjectLeaf<3>(src), radius0(src.radius0), radius1(src.radius1), sina(src.sina), cosa(src.cosa), height(src.height) {}
71
72EllipticCylinder::Box EllipticCylinder::getBoundingBox() const {
73 double c2 = cosa * cosa, s2 = sina * sina;
74 double a2 = radius0 * radius0, b2 = radius1 * radius1;
75 double x = sqrt(c2 * a2 + s2 * b2), y = sqrt(s2 * a2 + c2 * b2);
76 return Box(DVec(-x, -y, 0), DVec(x, y, height));
77}
78
79bool EllipticCylinder::contains(const EllipticCylinder::DVec& p) const {
80 if (0.0 > p.vert() || p.vert() > height) return false;
81 DVec p1 = invT(p);
82 double x = p1.c0 / radius0, y = p1.c1 / radius1;
83 return x * x + y * y <= 1.;
84}
85
86void EllipticCylinder::writeXMLAttr(XMLWriter::Element& dest_xml_object, const AxisNames& axes) const {
87 GeometryObjectLeaf<3>::writeXMLAttr(dest_xml_object, axes);
88 XMLWriter::Element& w = materialProvider->writeXML(dest_xml_object, axes);
89 w.attr("radius0", radius0).attr("radius1", radius1).attr("height", height);
90 if (sina != 0.) w.attr("angle", getAngle());
91}
92
93void EllipticCylinder::addPointsAlongToSet(std::set<double>& points,
94 Primitive<3>::Direction direction,
95 unsigned max_steps,
96 double min_step_size) const {
97 if (direction == Primitive<3>::DIRECTION_VERT) {
98 if (materialProvider->isUniform(Primitive<3>::DIRECTION_VERT)) {
99 points.insert(0);
100 points.insert(height);
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 unsigned steps = min(unsigned(height / min_step_size), max_steps);
105 double step = height / steps;
106 for (unsigned i = 0; i <= steps; ++i) points.insert(i * step);
107 }
108 } else {
109 auto bbox = getBoundingBox();
110 double b0, b1;
111 if (direction == Primitive<3>::DIRECTION_LONG) {
112 b0 = bbox.lower.lon();
113 b1 = bbox.upper.lon();
114 } else {
115 b0 = bbox.lower.tran();
116 b1 = bbox.upper.tran();
117 }
118 double d = b1 - b0;
119 if (this->max_steps) max_steps = this->max_steps;
120 if (this->min_step_size) min_step_size = this->min_step_size;
121 unsigned steps = min(unsigned(d / min_step_size), max_steps);
122 double step = d / steps;
123 for (unsigned i = 0; i <= steps; ++i) points.insert(b0 + i * step);
124 }
125}
126
127void EllipticCylinder::addLineSegmentsToSet(std::set<typename GeometryObjectD<3>::LineSegment>& segments,
128 unsigned max_steps,
129 double min_step_size) const {
130 typedef typename GeometryObjectD<3>::LineSegment Segment;
131 if (this->max_steps) max_steps = this->max_steps;
132 if (this->min_step_size) min_step_size = this->min_step_size;
133 unsigned steps = min(unsigned(M_PI * min(radius0, radius1) / min_step_size), max_steps);
134 double dphi = M_PI / steps;
135
136 // Make vertical axis
137 std::vector<double> zz;
138 if (materialProvider->isUniform(Primitive<3>::DIRECTION_VERT)) {
139 zz.reserve(2);
140 zz.push_back(0);
141 zz.push_back(height);
142 } else {
143 unsigned zsteps = min(unsigned(height / min_step_size), max_steps);
144 double zstep = height / zsteps;
145 zz.reserve(zsteps + 1);
146 for (unsigned i = 0; i <= zsteps; ++i) zz.push_back(i * zstep);
147 }
148
149 double z0 = 0.;
150 for (double z1 : zz) {
151 double x0 = radius0, y0 = 0;
152 if (z1 != 0.) {
153 segments.insert(Segment(TVec(-x0, -y0, z0), TVec(-x0, -y0, z1)));
154 segments.insert(Segment(TVec(x0, -y0, z0), TVec(x0, -y0, z1)));
155 segments.insert(Segment(TVec(-x0, y0, z0), TVec(-x0, y0, z1)));
156 segments.insert(Segment(TVec(x0, y0, z0), TVec(x0, y0, z1)));
157 }
158 for (unsigned i = 1; i <= (steps + 1) / 2; ++i) {
159 double phi = dphi * i;
160 double x1 = radius0 * cos(phi), y1 = radius0 * sin(phi);
161 segments.insert(Segment(TVec(-x0, -y0, z1), TVec(-x1, -y1, z1)));
162 segments.insert(Segment(TVec(x0, -y0, z1), TVec(x1, -y1, z1)));
163 segments.insert(Segment(TVec(-x0, y0, z1), TVec(-x1, y1, z1)));
164 segments.insert(Segment(TVec(x0, y0, z1), TVec(x1, y1, z1)));
165 if (z1 != 0.) {
166 segments.insert(Segment(TVec(-x1, -y1, z0), TVec(-x1, -y1, z1)));
167 segments.insert(Segment(TVec(x1, -y1, z0), TVec(x1, -y1, z1)));
168 segments.insert(Segment(TVec(-x1, y1, z0), TVec(-x1, y1, z1)));
169 segments.insert(Segment(TVec(x1, y1, z0), TVec(x1, y1, z1)));
170 }
171 if (x1 >= 0 && !materialProvider->isUniform(Primitive<3>::DIRECTION_LONG)) {
172 segments.insert(Segment(TVec(-x1, -y1, z1), TVec(-x1, y1, z1)));
173 segments.insert(Segment(TVec(x1, -y1, z1), TVec(x1, y1, z1)));
174 }
175 if (!materialProvider->isUniform(Primitive<3>::DIRECTION_TRAN)) {
176 segments.insert(Segment(TVec(-x1, -y1, z1), TVec(x1, -y1, z1)));
177 segments.insert(Segment(TVec(-x1, y1, z1), TVec(x1, y1, z1)));
178 }
179 x0 = x1;
180 y0 = y1;
181 }
182 z0 = z1;
183 }
184}
185
186shared_ptr<GeometryObject> read_elliptic_cylinder(GeometryReader& reader) {
187 shared_ptr<EllipticCylinder> result =
188 reader.manager.draft ? plask::make_shared<EllipticCylinder>(reader.source.getAttribute("radius0", 0.0),
189 reader.source.getAttribute("radius1", 0.0),
190 M_PI / 180. * reader.source.getAttribute("angle", 0.0),
191 reader.source.getAttribute("height", 0.0))
192 : plask::make_shared<EllipticCylinder>(reader.source.requireAttribute<double>("radius0"),
193 reader.source.requireAttribute<double>("radius1"),
194 M_PI / 180. * reader.source.getAttribute("angle", 0.0),
195 reader.source.requireAttribute<double>("height"));
196 result->readMaterial(reader);
197 reader.source.requireTagEnd();
198 return result;
199}
200
201static GeometryReader::RegisterObjectReader cylinder_reader(PLASK_ELLIPTIC_CYLINDER_NAME, read_elliptic_cylinder);
202
204
205} // namespace plask