PLaSK library
Loading...
Searching...
No Matches
circle.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 "circle.hpp"
15#include "../math.hpp"
16#include "../manager.hpp"
17#include "reader.hpp"
18
19#define PLASK_CIRCLE2D_NAME "circle"
20#define PLASK_CIRCLE3D_NAME "sphere"
21
22namespace plask {
23
24template <int dim> const char* Circle<dim>::NAME = dim == 2 ? PLASK_CIRCLE2D_NAME : PLASK_CIRCLE3D_NAME;
25
26template <int dim> std::string Circle<dim>::getTypeName() const { return NAME; }
27
28template <int dim>
29Circle<dim>::Circle(double radius, const shared_ptr<plask::Material>& material)
30 : GeometryObjectLeaf<dim>(material), radius(radius) {
31 if (radius < 0.) radius = 0.;
32}
33
34template <> typename Circle<2>::Box Circle<2>::getBoundingBox() const {
35 return Circle<2>::Box(vec(-radius, -radius), vec(radius, radius));
36}
37
38template <> typename Circle<3>::Box Circle<3>::getBoundingBox() const {
39 return Circle<3>::Box(vec(-radius, -radius, -radius), vec(radius, radius, radius));
40}
41
42template <int dim> bool Circle<dim>::contains(const typename Circle<dim>::DVec& p) const {
43 return abs2(p) <= radius * radius;
44}
45
46template <int dim> void Circle<dim>::writeXMLAttr(XMLWriter::Element& dest_xml_object, const AxisNames& axes) const {
47 GeometryObjectLeaf<dim>::writeXMLAttr(dest_xml_object, axes);
48 this->materialProvider->writeXML(dest_xml_object, axes).attr("radius", this->radius);
49}
50
51template <int dim>
52void Circle<dim>::addPointsAlongToSet(std::set<double>& points,
53 Primitive<3>::Direction direction,
54 unsigned max_steps,
55 double min_step_size) const {
56 assert(int(direction) >= 3 - dim && int(direction) <= 3);
57 if (this->max_steps) max_steps = this->max_steps;
58 if (this->min_step_size) min_step_size = this->min_step_size;
59 unsigned steps = min(unsigned(2. * radius / min_step_size), max_steps);
60 double step = 2. * radius / steps;
61 for (unsigned i = 0; i <= steps; ++i) points.insert(i * step - radius);
62}
63
64template <>
65void Circle<2>::addLineSegmentsToSet(std::set<typename GeometryObjectD<2>::LineSegment>& segments,
66 unsigned max_steps,
67 double min_step_size) const {
68 typedef typename GeometryObjectD<2>::LineSegment Segment;
69 if (this->max_steps) max_steps = this->max_steps;
70 if (this->min_step_size) min_step_size = this->min_step_size;
71 if (materialProvider->isUniform(Primitive<3>::DIRECTION_VERT)) {
72 unsigned steps = min(unsigned(M_PI * radius / min_step_size), max_steps);
73 double dphi = M_PI / steps;
74 double x0 = radius, y0 = 0;
75 for (unsigned i = 1; i <= (steps + 1) / 2; ++i) {
76 double phi = dphi * i;
77 double x1 = radius * cos(phi), y1 = radius * sin(phi);
78 segments.insert(Segment(DVec(-x0, -y0), DVec(-x1, -y1)));
79 segments.insert(Segment(DVec(x0, -y0), DVec(x1, -y1)));
80 segments.insert(Segment(DVec(-x0, y0), DVec(-x1, y1)));
81 segments.insert(Segment(DVec(x0, y0), DVec(x1, y1)));
82 if (x1 >= 0 && !materialProvider->isUniform(Primitive<3>::DIRECTION_TRAN)) {
83 segments.insert(Segment(DVec(-x1, -y1), DVec(-x1, y1)));
84 segments.insert(Segment(DVec(x1, -y1), DVec(x1, y1)));
85 }
86 x0 = x1;
87 y0 = y1;
88 }
89 } else {
90 // If material is not uniform wertically, we use uniform division in vertical direction
91 unsigned steps = min(unsigned(2. * radius / min_step_size), max_steps);
92 double step = 2. * radius / steps;
93 double radius2 = radius * radius;
94 double x0 = sqrt(0.5 * step * radius); // x0 = r sin(φ/2) = r √[(1–cosφ)/2], cosφ = (r-s)/r = 1 – s/r
95 double y0 = sqrt(radius2 - x0 * x0);
96 segments.insert(Segment(DVec(0., -radius), DVec(-x0, -y0)));
97 segments.insert(Segment(DVec(0., -radius), DVec(x0, -y0)));
98 segments.insert(Segment(DVec(0., radius), DVec(-x0, y0)));
99 segments.insert(Segment(DVec(0., radius), DVec(x0, y0)));
100 if (!materialProvider->isUniform(Primitive<3>::DIRECTION_TRAN)) {
101 segments.insert(Segment(DVec(-x0, -y0), DVec(-x0, y0)));
102 segments.insert(Segment(DVec(0., -radius), DVec(0., radius)));
103 segments.insert(Segment(DVec(x0, -y0), DVec(x0, y0)));
104 }
105 for (unsigned i = 1; i <= (steps + 1) / 2; ++i) {
106 double y1 = radius - i * step;
107 double x1 = sqrt(radius2 - y1 * y1);
108 segments.insert(Segment(DVec(-x1, -y1), DVec(x1, -y1)));
109 segments.insert(Segment(DVec(-x1, y1), DVec(x1, y1)));
110 segments.insert(Segment(DVec(-x0, -y0), DVec(-x1, -y1)));
111 segments.insert(Segment(DVec(x0, -y0), DVec(x1, -y1)));
112 segments.insert(Segment(DVec(-x0, y0), DVec(-x1, y1)));
113 segments.insert(Segment(DVec(x0, y0), DVec(x1, y1)));
114 if (!materialProvider->isUniform(Primitive<3>::DIRECTION_TRAN)) {
115 segments.insert(Segment(DVec(x1, -y1), DVec(x1, y1)));
116 segments.insert(Segment(DVec(-x1, -y1), DVec(-x1, y1)));
117 }
118 x0 = x1;
119 y0 = y1;
120 }
121 }
122}
123
124template <>
125void Circle<3>::addLineSegmentsToSet(std::set<typename GeometryObjectD<3>::LineSegment>& segments,
126 unsigned max_steps,
127 double min_step_size) const {
128 if (!materialProvider->isUniform(Primitive<3>::DIRECTION_LONG) ||
129 materialProvider->isUniform(Primitive<3>::DIRECTION_TRAN)) {
130 throw NotImplemented("triangular mesh for sphere non-uniform in any horizontal direction");
131 }
132 typedef typename GeometryObjectD<3>::LineSegment Segment;
133 if (this->max_steps) max_steps = this->max_steps;
134 if (this->min_step_size) min_step_size = this->min_step_size;
135 // if (materialProvider->isUniform(Primitive<3>::DIRECTION_VERT)) {
136 unsigned steps = min(unsigned(M_PI * radius / min_step_size), max_steps);
137 double dphi = M_PI / steps;
138 double r0 = radius, z0 = 0;
139 for (unsigned i = 0; i <= (steps + 1) / 2; ++i) {
140 double theta = dphi * i;
141 double r1 = radius * cos(theta), z1 = radius * sin(theta);
142 double x00 = r0, y00 = 0., x10 = r1, y10 = 0.;
143 for (unsigned j = 1; j <= 2 * steps; ++j) {
144 double phi = j * dphi;
145 double x01 = r0 * cos(phi), y01 = r0 * sin(phi), x11 = r1 * cos(phi), y11 = r1 * sin(phi);
146 if (i != 0) {
147 segments.insert(Segment(DVec(x01, y01, z0), DVec(x11, y11, z1)));
148 segments.insert(Segment(DVec(x01, y01, -z0), DVec(x11, y11, -z1)));
149 }
150 if (abs(r1) > SMALL) {
151 segments.insert(Segment(DVec(x10, y10, z1), DVec(x11, y11, z1)));
152 segments.insert(Segment(DVec(x10, y10, -z1), DVec(x11, y11, -z1)));
153 }
154 x00 = x01, y00 = y01, x10 = x11, y10 = y11;
155 }
156 r0 = r1;
157 z0 = z1;
158 }
159}
160
161template <int dim> shared_ptr<GeometryObject> read_circle(GeometryReader& reader) {
162 shared_ptr<Circle<dim>> circle =
163 plask::make_shared<Circle<dim>>(reader.manager.draft ? reader.source.getAttribute("radius", 0.0)
164 : reader.source.requireAttribute<double>("radius"));
165 circle->readMaterial(reader);
166 reader.source.requireTagEnd();
167 return circle;
168}
169
170template struct PLASK_API Circle<2>;
171template struct PLASK_API Circle<3>;
172
173static GeometryReader::RegisterObjectReader circle_reader(PLASK_CIRCLE2D_NAME, read_circle<2>);
174static GeometryReader::RegisterObjectReader sphere_reader(PLASK_CIRCLE3D_NAME, read_circle<3>);
175
176} // namespace plask