PLaSK library
Loading...
Searching...
No Matches
ellipse.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 "ellipse.hpp"
15#include "../manager.hpp"
16#include "../math.hpp"
17#include "reader.hpp"
18
19#define PLASK_ELLIPSE2D_NAME "ellipse"
20
21namespace plask {
22
23const char* Ellipse::NAME = PLASK_ELLIPSE2D_NAME;
24
25std::string Ellipse::getTypeName() const { return NAME; }
26
27Ellipse::Ellipse(double rx, double ry, const shared_ptr<plask::Material>& material)
28 : GeometryObjectLeaf(material), radius0(std::max(rx, 0.)), radius1(std::max(ry, 0.)) {}
29
32
33typename Ellipse::Box Ellipse::getBoundingBox() const { return Ellipse::Box(vec(-radius0, -radius1), vec(radius0, radius1)); }
34
35bool Ellipse::contains(const typename Ellipse::DVec& p) const {
36 double x = p.c0 / radius0, y = p.c1 / radius1;
37 return x * x + y * y <= 1.;
38}
39
40void Ellipse::writeXMLAttr(XMLWriter::Element& dest_xml_object, const AxisNames& axes) const {
41 GeometryObjectLeaf::writeXMLAttr(dest_xml_object, axes);
42 this->materialProvider->writeXML(dest_xml_object, axes).attr("radius0", this->radius0).attr("radius1", this->radius1);
43}
44
45void Ellipse::addPointsAlongToSet(std::set<double>& points,
46 Primitive<3>::Direction direction,
47 unsigned max_steps,
48 double min_step_size) const {
49 assert(int(direction) >= 1 && int(direction) <= 3);
50 double radius = (direction == Primitive<3>::DIRECTION_VERT) ? radius1 : radius0;
51 if (this->max_steps) max_steps = this->max_steps;
52 if (this->min_step_size) min_step_size = this->min_step_size;
53 unsigned steps = min(unsigned(2. * radius / min_step_size), max_steps);
54 double step = 2. * radius / steps;
55 for (unsigned i = 0; i <= steps; ++i) points.insert(i * step - radius);
56}
57
58void Ellipse::addLineSegmentsToSet(std::set<typename GeometryObjectD<2>::LineSegment>& segments,
59 unsigned max_steps,
60 double min_step_size) const {
61 typedef typename GeometryObjectD<2>::LineSegment Segment;
62 if (this->max_steps) max_steps = this->max_steps;
63 if (this->min_step_size) min_step_size = this->min_step_size;
64 if (materialProvider->isUniform(Primitive<3>::DIRECTION_VERT)) {
65 unsigned steps = min(unsigned(M_PI * min(radius0, radius1) / min_step_size), max_steps);
66 double dphi = M_PI / steps;
67 double x0 = radius0, y0 = 0;
68 for (unsigned i = 1; i <= (steps + 1) / 2; ++i) {
69 double phi = dphi * i;
70 double x1 = radius0 * cos(phi), y1 = radius1 * sin(phi);
71 segments.insert(Segment(DVec(-x0, -y0), DVec(-x1, -y1)));
72 segments.insert(Segment(DVec(x0, -y0), DVec(x1, -y1)));
73 segments.insert(Segment(DVec(-x0, y0), DVec(-x1, y1)));
74 segments.insert(Segment(DVec(x0, y0), DVec(x1, y1)));
75 if (x1 >= 0 && !materialProvider->isUniform(Primitive<3>::DIRECTION_TRAN)) {
76 segments.insert(Segment(DVec(-x1, -y1), DVec(-x1, y1)));
77 segments.insert(Segment(DVec(x1, -y1), DVec(x1, y1)));
78 }
79 x0 = x1;
80 y0 = y1;
81 }
82 } else {
83 // If material is not uniform vertically, we use uniform division in vertical direction
84 unsigned steps = min(unsigned(2. * radius1 / min_step_size), max_steps);
85 double step = 2. * radius1 / steps;
86 double x0 = sqrt(0.5 * step * radius0); // x0 = r sin(φ/2) = r √[(1–cosφ)/2], cosφ = (r-s)/r = 1 – s/r
87 double xr0 = x0 / radius0;
88 double y0 = radius1 * sqrt(1. - xr0 * xr0);
89 segments.insert(Segment(DVec(0., -radius1), DVec(-x0, -y0)));
90 segments.insert(Segment(DVec(0., -radius1), DVec(x0, -y0)));
91 segments.insert(Segment(DVec(0., radius1), DVec(-x0, y0)));
92 segments.insert(Segment(DVec(0., radius1), DVec(x0, y0)));
93 if (!materialProvider->isUniform(Primitive<3>::DIRECTION_TRAN)) {
94 segments.insert(Segment(DVec(-x0, -y0), DVec(-x0, y0)));
95 segments.insert(Segment(DVec(0., -radius1), DVec(0., radius1)));
96 segments.insert(Segment(DVec(x0, -y0), DVec(x0, y0)));
97 }
98 for (unsigned i = 1; i <= (steps + 1) / 2; ++i) {
99 double y1 = radius1 - i * step;
100 double yr1 = y1 / radius1;
101 double x1 = radius0 * sqrt(1. - yr1 * yr1);
102 segments.insert(Segment(DVec(-x1, -y1), DVec(x1, -y1)));
103 segments.insert(Segment(DVec(-x1, y1), DVec(x1, y1)));
104 segments.insert(Segment(DVec(-x0, -y0), DVec(-x1, -y1)));
105 segments.insert(Segment(DVec(x0, -y0), DVec(x1, -y1)));
106 segments.insert(Segment(DVec(-x0, y0), DVec(-x1, y1)));
107 segments.insert(Segment(DVec(x0, y0), DVec(x1, y1)));
108 if (!materialProvider->isUniform(Primitive<3>::DIRECTION_TRAN)) {
109 segments.insert(Segment(DVec(x1, -y1), DVec(x1, y1)));
110 segments.insert(Segment(DVec(-x1, -y1), DVec(-x1, y1)));
111 }
112 x0 = x1;
113 y0 = y1;
114 }
115 }
116}
117
118shared_ptr<GeometryObject> read_ellipse(GeometryReader& reader) {
119 shared_ptr<Ellipse> ellipse =
120 reader.manager.draft
121 ? plask::make_shared<Ellipse>(reader.source.getAttribute("radius0", 0.0), reader.source.getAttribute("radius1", 0.0))
122 : plask::make_shared<Ellipse>(reader.source.requireAttribute<double>("radius0"),
123 reader.source.requireAttribute<double>("radius1"));
124 ellipse->readMaterial(reader);
125 reader.source.requireTagEnd();
126 return ellipse;
127}
128
129static GeometryReader::RegisterObjectReader ellipse_reader(PLASK_ELLIPSE2D_NAME, read_ellipse);
130
131} // namespace plask