PLaSK library
Loading...
Searching...
No Matches
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 "cylinder.hpp"
15#include "../manager.hpp"
16#include "reader.hpp"
17
18#define PLASK_CYLINDER_NAME "cylinder"
19#define PLASK_HOLLOW_CYLINDER_NAME "tube"
20
21namespace plask {
22
24
25const char* Cylinder::NAME = PLASK_CYLINDER_NAME;
26
27Cylinder::Cylinder(double radius, double height, const shared_ptr<Material>& material)
28 : GeometryObjectLeaf<3>(material), radius(std::max(radius, 0.)), height(std::max(height, 0.)) {}
29
30Cylinder::Cylinder(double radius, double height, shared_ptr<MaterialsDB::MixedCompositionFactory> materialTopBottom)
31 : GeometryObjectLeaf<3>(materialTopBottom), radius(std::max(radius, 0.)), height(std::max(height, 0.)) {}
32
33Cylinder::Cylinder(const Cylinder& src) : GeometryObjectLeaf<3>(src), radius(src.radius), height(src.height) {}
34
35Cylinder::Box Cylinder::getBoundingBox() const { return Box(vec(-radius, -radius, 0.0), vec(radius, radius, height)); }
36
37bool Cylinder::contains(const Cylinder::DVec& p) const {
38 return 0.0 <= p.vert() && p.vert() <= height && std::fma(p.lon(), p.lon(), p.tran() * p.tran()) <= radius * radius;
39}
40
45
46void Cylinder::addPointsAlongToSet(std::set<double>& points,
47 Primitive<3>::Direction direction,
48 unsigned max_steps,
49 double min_step_size) const {
50 if (direction == Primitive<3>::DIRECTION_VERT) {
51 if (materialProvider->isUniform(Primitive<3>::DIRECTION_VERT)) {
52 points.insert(0);
53 points.insert(height);
54 } else {
55 if (this->max_steps) max_steps = this->max_steps;
56 if (this->min_step_size) min_step_size = this->min_step_size;
57 unsigned steps = min(unsigned(height / min_step_size), max_steps);
58 double step = height / steps;
59 for (unsigned i = 0; i <= steps; ++i) points.insert(i * step);
60 }
61 } else {
62 if (this->max_steps) max_steps = this->max_steps;
63 if (this->min_step_size) min_step_size = this->min_step_size;
64 unsigned steps = min(unsigned(2. * radius / min_step_size), max_steps);
65 double step = 2. * radius / steps;
66 for (unsigned i = 0; i <= steps; ++i) points.insert(i * step - radius);
67 }
68}
69
70void Cylinder::addLineSegmentsToSet(std::set<typename GeometryObjectD<3>::LineSegment>& segments,
71 unsigned max_steps,
72 double min_step_size) const {
73 typedef typename GeometryObjectD<3>::LineSegment Segment;
74 if (this->max_steps) max_steps = this->max_steps;
75 if (this->min_step_size) min_step_size = this->min_step_size;
76 unsigned steps = min(unsigned(M_PI * radius / min_step_size), max_steps);
77 double dphi = M_PI / steps;
78
79 // Make vertical axis
80 std::vector<double> zz;
81 if (materialProvider->isUniform(Primitive<3>::DIRECTION_VERT)) {
82 zz.reserve(2);
83 zz.push_back(0);
84 zz.push_back(height);
85 } else {
86 unsigned zsteps = min(unsigned(height / min_step_size), max_steps);
87 double zstep = height / zsteps;
88 zz.reserve(zsteps + 1);
89 for (unsigned i = 0; i <= zsteps; ++i) zz.push_back(i * zstep);
90 }
91
92 double z0 = 0.;
93 for (double z1 : zz) {
94 double x0 = radius, y0 = 0;
95 if (z1 != 0.) {
96 segments.insert(Segment(DVec(-x0, -y0, z0), DVec(-x0, -y0, z1)));
97 segments.insert(Segment(DVec(x0, -y0, z0), DVec(x0, -y0, z1)));
98 segments.insert(Segment(DVec(-x0, y0, z0), DVec(-x0, y0, z1)));
99 segments.insert(Segment(DVec(x0, y0, z0), DVec(x0, y0, z1)));
100 }
101 for (unsigned i = 1; i <= (steps + 1) / 2; ++i) {
102 double phi = dphi * i;
103 double x1 = radius * cos(phi), y1 = radius * sin(phi);
104 segments.insert(Segment(DVec(-x0, -y0, z1), DVec(-x1, -y1, z1)));
105 segments.insert(Segment(DVec(x0, -y0, z1), DVec(x1, -y1, z1)));
106 segments.insert(Segment(DVec(-x0, y0, z1), DVec(-x1, y1, z1)));
107 segments.insert(Segment(DVec(x0, y0, z1), DVec(x1, y1, z1)));
108 if (z1 != 0.) {
109 segments.insert(Segment(DVec(-x1, -y1, z0), DVec(-x1, -y1, z1)));
110 segments.insert(Segment(DVec(x1, -y1, z0), DVec(x1, -y1, z1)));
111 segments.insert(Segment(DVec(-x1, y1, z0), DVec(-x1, y1, z1)));
112 segments.insert(Segment(DVec(x1, y1, z0), DVec(x1, y1, z1)));
113 }
114 if (x1 >= 0 && !materialProvider->isUniform(Primitive<3>::DIRECTION_LONG)) {
115 segments.insert(Segment(DVec(-x1, -y1, z1), DVec(-x1, y1, z1)));
116 segments.insert(Segment(DVec(x1, -y1, z1), DVec(x1, y1, z1)));
117 }
118 if (!materialProvider->isUniform(Primitive<3>::DIRECTION_TRAN)) {
119 segments.insert(Segment(DVec(-x1, -y1, z1), DVec(x1, -y1, z1)));
120 segments.insert(Segment(DVec(-x1, y1, z1), DVec(x1, y1, z1)));
121 }
122 x0 = x1;
123 y0 = y1;
124 }
125 z0 = z1;
126 }
127}
128
129shared_ptr<GeometryObject> read_cylinder(GeometryReader& reader) {
130 shared_ptr<Cylinder> result(new Cylinder(
131 reader.manager.draft ? reader.source.getAttribute("radius", 0.0) : reader.source.requireAttribute<double>("radius"),
132 reader.manager.draft ? reader.source.getAttribute("height", 0.0) : reader.source.requireAttribute<double>("height")));
133 result->readMaterial(reader);
134 reader.source.requireTagEnd();
135 return result;
136}
137
138static GeometryReader::RegisterObjectReader cylinder_reader(PLASK_CYLINDER_NAME, read_cylinder);
139
141
142const char* HollowCylinder::NAME = PLASK_CYLINDER_NAME;
143
144HollowCylinder::HollowCylinder(double inner_radius, double outer_radius, double height, const shared_ptr<Material>& material)
145 : GeometryObjectLeaf<3>(material),
146 inner_radius(std::max(inner_radius, 0.)),
147 outer_radius(std::max(outer_radius, 0.)),
148 height(std::max(height, 0.)) {
149 if (inner_radius > outer_radius) throw BadInput("Tube", "Inner radius must be less than outer radius");
150}
151
152HollowCylinder::HollowCylinder(double inner_radius,
153 double outer_radius,
154 double height,
155 shared_ptr<MaterialsDB::MixedCompositionFactory> materialTopBottom)
156 : GeometryObjectLeaf<3>(materialTopBottom),
157 inner_radius(std::max(inner_radius, 0.)),
158 outer_radius(std::max(outer_radius, 0.)),
159 height(std::max(height, 0.)) {
160 if (inner_radius > outer_radius) throw BadInput("Tube", "Inner radius must be less than outer radius");
161}
162
163HollowCylinder::HollowCylinder(const HollowCylinder& src)
164 : GeometryObjectLeaf<3>(src), inner_radius(src.inner_radius), outer_radius(src.outer_radius), height(src.height) {}
165
169
170bool HollowCylinder::contains(const HollowCylinder::DVec& p) const {
171 return 0.0 <= p.vert() && p.vert() <= height &&
172 std::fma(p.lon(), p.lon(), p.tran() * p.tran()) <= outer_radius * outer_radius &&
173 std::fma(p.lon(), p.lon(), p.tran() * p.tran()) >= inner_radius * inner_radius;
174}
175
176void HollowCylinder::writeXMLAttr(XMLWriter::Element& dest_xml_object, const AxisNames& axes) const {
177 GeometryObjectLeaf<3>::writeXMLAttr(dest_xml_object, axes);
178 materialProvider->writeXML(dest_xml_object, axes)
179 .attr("inner-radius", inner_radius)
180 .attr("outer-radius", outer_radius)
181 .attr("height", height);
182}
183
184void HollowCylinder::addPointsAlongToSet(std::set<double>& points,
185 Primitive<3>::Direction direction,
186 unsigned max_steps,
187 double min_step_size) const {
188 if (direction == Primitive<3>::DIRECTION_VERT) {
189 if (materialProvider->isUniform(Primitive<3>::DIRECTION_VERT)) {
190 points.insert(0);
191 points.insert(height);
192 } else {
193 if (this->max_steps) max_steps = this->max_steps;
194 if (this->min_step_size) min_step_size = this->min_step_size;
195 unsigned steps = min(unsigned(height / min_step_size), max_steps);
196 double step = height / steps;
197 for (unsigned i = 0; i <= steps; ++i) points.insert(i * step);
198 }
199 } else {
200 if (this->max_steps) max_steps = this->max_steps;
201 if (this->min_step_size) min_step_size = this->min_step_size;
202 int max_steps_2 = int(round(max_steps * inner_radius / outer_radius));
203 int max_steps_1 = max_steps - max_steps_2 / 2;
204 int steps1 = min(int((outer_radius - inner_radius) / min_step_size), max_steps_1);
205 int steps2 = min(int(2 * inner_radius / min_step_size), max_steps_2);
206 double step1 = (outer_radius - inner_radius) / steps1;
207 double step2 = 2. * inner_radius / steps2;
208 for (int i = 0; i < steps1; ++i) points.insert(i * step1 - outer_radius);
209 for (int i = 0; i < steps2; ++i) points.insert(i * step2 - inner_radius);
210 for (int i = steps1; i >= 0; ++i) points.insert(outer_radius - i * step1);
211 }
212}
213
214void HollowCylinder::addLineSegmentsToSet(std::set<typename GeometryObjectD<3>::LineSegment>& segments,
215 unsigned max_steps,
216 double min_step_size) const {
217 typedef typename GeometryObjectD<3>::LineSegment Segment;
218 if (this->max_steps) max_steps = this->max_steps;
219 if (this->min_step_size) min_step_size = this->min_step_size;
220 unsigned steps = min(unsigned(M_PI * inner_radius / min_step_size), max_steps);
221 double dphi = M_PI / steps;
222
223 // Make vertical axis
224 std::vector<double> zz;
225 if (materialProvider->isUniform(Primitive<3>::DIRECTION_VERT)) {
226 zz.reserve(2);
227 zz.push_back(0);
228 zz.push_back(height);
229 } else {
230 unsigned zsteps = min(unsigned(height / min_step_size), max_steps);
231 double zstep = height / zsteps;
232 zz.reserve(zsteps + 1);
233 for (unsigned i = 0; i <= zsteps; ++i) zz.push_back(i * zstep);
234 }
235
236 double z0 = 0.;
237 for (double z1 : zz) {
238 double x0 = outer_radius, y0 = 0;
239 if (z1 != 0.) {
240 segments.insert(Segment(DVec(-x0, -y0, z0), DVec(-x0, -y0, z1)));
241 segments.insert(Segment(DVec(x0, -y0, z0), DVec(x0, -y0, z1)));
242 segments.insert(Segment(DVec(-x0, y0, z0), DVec(-x0, y0, z1)));
243 segments.insert(Segment(DVec(x0, y0, z0), DVec(x0, y0, z1)));
244 }
245 for (unsigned i = 1; i <= (steps + 1) / 2; ++i) {
246 double phi = dphi * i;
247 double x1 = outer_radius * cos(phi), y1 = outer_radius * sin(phi);
248 segments.insert(Segment(DVec(-x0, -y0, z1), DVec(-x1, -y1, z1)));
249 segments.insert(Segment(DVec(x0, -y0, z1), DVec(x1, -y1, z1)));
250 segments.insert(Segment(DVec(-x0, y0, z1), DVec(-x1, y1, z1)));
251 segments.insert(Segment(DVec(x0, y0, z1), DVec(x1, y1, z1)));
252 if (z1 != 0.) {
253 segments.insert(Segment(DVec(-x1, -y1, z0), DVec(-x1, -y1, z1)));
254 segments.insert(Segment(DVec(x1, -y1, z0), DVec(x1, -y1, z1)));
255 segments.insert(Segment(DVec(-x1, y1, z0), DVec(-x1, y1, z1)));
256 segments.insert(Segment(DVec(x1, y1, z0), DVec(x1, y1, z1)));
257 }
258 if (x1 >= 0 && !materialProvider->isUniform(Primitive<3>::DIRECTION_LONG)) {
259 segments.insert(Segment(DVec(-x1, -y1, z1), DVec(-x1, y1, z1)));
260 segments.insert(Segment(DVec(x1, -y1, z1), DVec(x1, y1, z1)));
261 }
262 if (!materialProvider->isUniform(Primitive<3>::DIRECTION_TRAN)) {
263 segments.insert(Segment(DVec(-x1, -y1, z1), DVec(x1, -y1, z1)));
264 segments.insert(Segment(DVec(-x1, y1, z1), DVec(x1, y1, z1)));
265 }
266 x0 = x1;
267 y0 = y1;
268 }
269 z0 = z1;
270 }
271
272 z0 = 0.;
273 for (double z1 : zz) {
274 double x0 = inner_radius, y0 = 0;
275 if (z1 != 0.) {
276 segments.insert(Segment(DVec(-x0, -y0, z0), DVec(-x0, -y0, z1)));
277 segments.insert(Segment(DVec(x0, -y0, z0), DVec(x0, -y0, z1)));
278 segments.insert(Segment(DVec(-x0, y0, z0), DVec(-x0, y0, z1)));
279 segments.insert(Segment(DVec(x0, y0, z0), DVec(x0, y0, z1)));
280 }
281 for (unsigned i = 1; i <= (steps + 1) / 2; ++i) {
282 double phi = dphi * i;
283 double x1 = inner_radius * cos(phi), y1 = inner_radius * sin(phi);
284 segments.insert(Segment(DVec(-x0, -y0, z1), DVec(-x1, -y1, z1)));
285 segments.insert(Segment(DVec(x0, -y0, z1), DVec(x1, -y1, z1)));
286 segments.insert(Segment(DVec(-x0, y0, z1), DVec(-x1, y1, z1)));
287 segments.insert(Segment(DVec(x0, y0, z1), DVec(x1, y1, z1)));
288 if (z1 != 0.) {
289 segments.insert(Segment(DVec(-x1, -y1, z0), DVec(-x1, -y1, z1)));
290 segments.insert(Segment(DVec(x1, -y1, z0), DVec(x1, -y1, z1)));
291 segments.insert(Segment(DVec(-x1, y1, z0), DVec(-x1, y1, z1)));
292 segments.insert(Segment(DVec(x1, y1, z0), DVec(x1, y1, z1)));
293 }
294 if (x1 >= 0 && !materialProvider->isUniform(Primitive<3>::DIRECTION_LONG)) {
295 segments.insert(Segment(DVec(-x1, -y1, z1), DVec(-x1, y1, z1)));
296 segments.insert(Segment(DVec(x1, -y1, z1), DVec(x1, y1, z1)));
297 }
298 if (!materialProvider->isUniform(Primitive<3>::DIRECTION_TRAN)) {
299 segments.insert(Segment(DVec(-x1, -y1, z1), DVec(x1, -y1, z1)));
300 segments.insert(Segment(DVec(-x1, y1, z1), DVec(x1, y1, z1)));
301 }
302 x0 = x1;
303 y0 = y1;
304 }
305 z0 = z1;
306 }
307}
308
309shared_ptr<GeometryObject> read_hollow_cylinder(GeometryReader& reader) {
310 double inner_radius = reader.manager.draft ? reader.source.getAttribute("inner-radius", 0.0)
311 : reader.source.requireAttribute<double>("inner-radius");
312 double outer_radius = reader.manager.draft ? reader.source.getAttribute("outer-radius", 0.0)
313 : reader.source.requireAttribute<double>("outer-radius");
314 if (reader.manager.draft && inner_radius > outer_radius) inner_radius = outer_radius;
315 shared_ptr<HollowCylinder> result(new HollowCylinder(
316 inner_radius, outer_radius,
317 reader.manager.draft ? reader.source.getAttribute("height", 0.0) : reader.source.requireAttribute<double>("height")));
318 result->readMaterial(reader);
319 reader.source.requireTagEnd();
320 return result;
321}
322
323static GeometryReader::RegisterObjectReader hollow_cylinder_reader(PLASK_HOLLOW_CYLINDER_NAME, read_hollow_cylinder);
324
326
327} // namespace plask