Line data Source code
1 : /*
2 : * Copyright (c) 2016: G-CSC, Goethe University Frankfurt
3 : * Author: Markus Breit
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #include "elliptic_cylinder_projector.h"
34 :
35 : #include "common/math/misc/math_util.h"
36 : #include "refinement_projector.h"
37 :
38 :
39 : namespace ug {
40 :
41 :
42 0 : EllipticCylinderProjector::EllipticCylinderProjector()
43 : : m_center(0, 0, 0),
44 : m_cylinder_axis(0, 0, 1),
45 : m_ellipse_axis1(1, 0, 0),
46 : m_ellipse_axis2(0, 1, 0),
47 0 : m_radius(-1),
48 0 : m_influenceRadius(-1),
49 : m_ellipseNormal(0, 0, 1),
50 0 : m_a(1.0),
51 0 : m_b(1.0)
52 0 : {}
53 :
54 :
55 0 : EllipticCylinderProjector::EllipticCylinderProjector
56 : (
57 : const vector3& center,
58 : const vector3& cylAxis,
59 : const vector3& ellipseAxis1,
60 : const vector3& ellipseAxis2
61 0 : )
62 : : m_center(center),
63 : m_cylinder_axis(cylAxis),
64 : m_ellipse_axis1(ellipseAxis1),
65 : m_ellipse_axis2(ellipseAxis2),
66 0 : m_radius(-1),
67 0 : m_influenceRadius(-1),
68 0 : m_a(VecLength(ellipseAxis1)),
69 0 : m_b(VecLength(ellipseAxis2))
70 : {
71 0 : VecCross(m_ellipseNormal, m_ellipse_axis1, m_ellipse_axis2);
72 0 : }
73 :
74 :
75 0 : EllipticCylinderProjector::EllipticCylinderProjector
76 : (
77 : const vector3& center,
78 : const vector3& cylAxis,
79 : const vector3& ellipseAxis1,
80 : const vector3& ellipseAxis2,
81 : number radius
82 0 : )
83 : : m_center(center),
84 : m_cylinder_axis(cylAxis),
85 : m_ellipse_axis1(ellipseAxis1),
86 : m_ellipse_axis2(ellipseAxis2),
87 0 : m_radius(radius),
88 0 : m_influenceRadius(-1),
89 0 : m_a(VecLength(ellipseAxis1)),
90 0 : m_b(VecLength(ellipseAxis2))
91 : {
92 : UG_LOGN("Constructor5");
93 0 : VecCross(m_ellipseNormal, m_ellipse_axis1, m_ellipse_axis2);
94 0 : }
95 :
96 :
97 0 : EllipticCylinderProjector::EllipticCylinderProjector
98 : (
99 : const vector3& center,
100 : const vector3& cylAxis,
101 : const vector3& ellipseAxis1,
102 : const vector3& ellipseAxis2,
103 : number radius,
104 : number influenceRadius
105 0 : )
106 : : m_center(center),
107 : m_cylinder_axis(cylAxis),
108 : m_ellipse_axis1(ellipseAxis1),
109 : m_ellipse_axis2(ellipseAxis2),
110 0 : m_radius(radius),
111 0 : m_influenceRadius(influenceRadius),
112 0 : m_a(VecLength(ellipseAxis1)),
113 0 : m_b(VecLength(ellipseAxis2))
114 : {
115 0 : VecCross(m_ellipseNormal, m_ellipse_axis1, m_ellipse_axis2);
116 0 : }
117 :
118 :
119 0 : EllipticCylinderProjector::EllipticCylinderProjector
120 : (
121 : SPIGeometry3d geometry,
122 : const vector3& center,
123 : const vector3& cylAxis,
124 : const vector3& ellipseAxis1,
125 : const vector3& ellipseAxis2,
126 : number radius,
127 : number influenceRadius
128 0 : )
129 : : RefinementProjector(geometry),
130 : m_center(center),
131 : m_cylinder_axis(cylAxis),
132 : m_ellipse_axis1(ellipseAxis1),
133 : m_ellipse_axis2(ellipseAxis2),
134 0 : m_radius(radius),
135 0 : m_influenceRadius(influenceRadius),
136 0 : m_a(VecLength(ellipseAxis1)),
137 0 : m_b(VecLength(ellipseAxis2))
138 : {
139 0 : VecCross(m_ellipseNormal, m_ellipse_axis1, m_ellipse_axis2);
140 0 : }
141 :
142 :
143 0 : EllipticCylinderProjector::~EllipticCylinderProjector()
144 0 : {}
145 :
146 :
147 :
148 0 : void EllipticCylinderProjector::set_center(const vector3& center)
149 : {
150 : m_center = center;
151 0 : }
152 :
153 0 : const vector3& EllipticCylinderProjector::center() const
154 : {
155 0 : return m_center;
156 : }
157 :
158 :
159 0 : void EllipticCylinderProjector::set_cylinder_axis(const vector3& axis)
160 : {
161 : m_cylinder_axis = axis;
162 0 : }
163 :
164 0 : const vector3& EllipticCylinderProjector::cylinder_axis() const
165 : {
166 0 : return m_cylinder_axis;
167 : }
168 :
169 :
170 0 : void EllipticCylinderProjector::set_ellipse_axis1(const vector3& axis)
171 : {
172 : m_ellipse_axis1 = axis;
173 0 : VecCross(m_ellipseNormal, m_ellipse_axis1, m_ellipse_axis2);
174 0 : m_a = VecLength(axis);
175 0 : }
176 :
177 0 : const vector3& EllipticCylinderProjector::ellipse_axis1() const
178 : {
179 0 : return m_ellipse_axis1;
180 : }
181 :
182 :
183 0 : void EllipticCylinderProjector::set_ellipse_axis2(const vector3& axis)
184 : {
185 : m_ellipse_axis2 = axis;
186 0 : VecCross(m_ellipseNormal, m_ellipse_axis1, m_ellipse_axis2);
187 0 : m_b = VecLength(axis);
188 0 : }
189 :
190 0 : const vector3& EllipticCylinderProjector::ellipse_axis2() const
191 : {
192 0 : return m_ellipse_axis2;
193 : }
194 :
195 :
196 0 : void EllipticCylinderProjector::set_radius(number radius)
197 : {
198 0 : m_radius = radius;
199 0 : }
200 :
201 0 : number EllipticCylinderProjector::radius() const
202 : {
203 0 : return m_radius;
204 : }
205 :
206 :
207 0 : void EllipticCylinderProjector::set_influence_radius(number influenceRadius)
208 : {
209 0 : m_influenceRadius = influenceRadius;
210 0 : }
211 :
212 0 : number EllipticCylinderProjector::influence_radius() const
213 : {
214 0 : return m_influenceRadius;
215 : }
216 :
217 :
218 :
219 0 : number EllipticCylinderProjector::new_vertex(Vertex* vrt, Edge* parent)
220 : {
221 0 : return perform_projection(vrt, parent);
222 : }
223 :
224 0 : number EllipticCylinderProjector::new_vertex(Vertex* vrt, Face* parent)
225 : {
226 0 : return perform_projection(vrt, parent);
227 : }
228 :
229 0 : number EllipticCylinderProjector::new_vertex(Vertex* vrt, Volume* parent)
230 : {
231 0 : return perform_projection(vrt, parent);
232 : }
233 :
234 :
235 :
236 0 : number EllipticCylinderProjector::radial_ellipse_coord(const vector3& v)
237 : {
238 : // project coordinates along cylinder axis to ellipse plane
239 : vector3 proj2Ellipse;
240 : number dummy;
241 0 : RayPlaneIntersection(proj2Ellipse, dummy, v, m_cylinder_axis, m_center, m_ellipseNormal);
242 :
243 : // get x and y coordinates
244 0 : const number x = VecDot(m_ellipse_axis1, proj2Ellipse) / m_a;
245 0 : const number y = VecDot(m_ellipse_axis2, proj2Ellipse) / m_b;
246 :
247 0 : return sqrt(x*x/(m_a*m_a) + y*y/(m_b*m_b));
248 : }
249 :
250 :
251 0 : number EllipticCylinderProjector::scale_point_to_radius(vector3& vIO, number r)
252 : {
253 : // project coordinates along cylinder axis to ellipse plane
254 : vector3 proj2Ellipse;
255 : number dummy;
256 0 : RayPlaneIntersection(proj2Ellipse, dummy, vIO, m_cylinder_axis, m_center, m_ellipseNormal);
257 :
258 : VecSubtract(vIO, vIO, proj2Ellipse);
259 :
260 : // current radius
261 0 : const number x = VecDot(m_ellipse_axis1, proj2Ellipse) / m_a;
262 0 : const number y = VecDot(m_ellipse_axis2, proj2Ellipse) / m_b;
263 0 : const number rcur = sqrt(x*x/(m_a*m_a) + y*y/(m_b*m_b));
264 :
265 : // scale and undo projection
266 0 : if (rcur > SMALL * r)
267 0 : VecScaleAdd(vIO, 1.0, vIO, r / rcur, proj2Ellipse);
268 : else
269 : // if current position is in the center of the cylinder, leave it there
270 : VecAdd(vIO, vIO, m_center);
271 :
272 0 : return rcur;
273 : }
274 :
275 :
276 : template <class TElem>
277 0 : number EllipticCylinderProjector::perform_projection(Vertex* vrt, TElem* parent)
278 : {
279 : // General method:
280 : // The "radius" of all parent vertices is averaged
281 : // and used for the radius of the projected child.
282 : // The (Cartesian) coordinates of all parents are averaged.
283 : // The resulting coordinates are scaled around the axis
284 : // within the ellipse plane to reach the previously determined radius.
285 : // TODO: It might be better to also average arc lengths and axial coords of parents
286 : // and project the child to averaged radius, arc length and axial coords.
287 :
288 : // average parent vertex positions and radii
289 0 : typename TElem::ConstVertexArray vrts = parent->vertices();
290 0 : const size_t numVrts = parent->num_vertices();
291 :
292 0 : if (numVrts == 0)
293 : {
294 0 : set_pos(vrt, vector3(0, 0, 0));
295 0 : return 1;
296 : }
297 :
298 : number avgR = 0;
299 : vector3 proj(0, 0, 0);
300 0 : for (size_t i = 0; i < numVrts; ++i)
301 : {
302 0 : const vector3& p = pos(vrts[i]);
303 0 : avgR += radial_ellipse_coord(p);
304 : proj += p;
305 : }
306 0 : avgR /= numVrts;
307 0 : VecScale(proj, proj, 1.0 / numVrts);
308 :
309 : // move averaged position to new radius
310 0 : const number curR = scale_point_to_radius(proj, avgR);
311 : set_pos(vrt, proj);
312 :
313 0 : if (m_influenceRadius > 0)
314 : {
315 0 : if (m_radius > m_influenceRadius)
316 : {
317 0 : const number dist = m_radius - m_influenceRadius;
318 0 : return clip<number>((curR - m_influenceRadius) / dist, 0, 1);
319 : }
320 0 : else if (m_radius >= 0)
321 : {
322 0 : const number dist = m_influenceRadius - m_radius;
323 0 : if (dist > 0)
324 0 : return clip<number>(1 - (curR - m_radius) / dist, 0, 1);
325 0 : return curR < m_radius ? 1 : 0;
326 : }
327 : else
328 0 : return clip<number>(1 - curR / m_influenceRadius, 0, 1);
329 : }
330 :
331 : return 1;
332 : }
333 :
334 :
335 : } // namespace ug
|
