Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "CircularBoundaryCorrectionGenerator.h"
11 :
12 : #include "MooseUtils.h"
13 : #include "LinearInterpolation.h"
14 : #include "FillBetweenPointVectorsTools.h"
15 : #include "CastUniquePointer.h"
16 :
17 : // C++ includes
18 : #include <cmath> // for atan2
19 :
20 : registerMooseObject("MooseApp", CircularBoundaryCorrectionGenerator);
21 :
22 : InputParameters
23 14381 : CircularBoundaryCorrectionGenerator::validParams()
24 : {
25 14381 : InputParameters params = MeshGenerator::validParams();
26 :
27 14381 : params.addRequiredParam<MeshGeneratorName>("input", "The input mesh to be modified.");
28 :
29 43143 : params.addParam<bool>("move_end_nodes_in_span_direction",
30 28762 : false,
31 : "Whether to move the end nodes in the span direction of the arc or not.");
32 :
33 14381 : params.addRequiredParam<std::vector<BoundaryName>>("input_mesh_circular_boundaries",
34 : "Names of the circular boundaries.");
35 :
36 14381 : params.addRangeCheckedParam<std::vector<Real>>(
37 : "transition_layer_ratios",
38 : "transition_layer_ratios>0&transition_layer_ratios<=1.0",
39 : "Ratios to radii as the transition layers on both sides of the original radii (if this "
40 : "parameter is not provided, 0.01 will be used).");
41 :
42 14381 : params.addRangeCheckedParam<std::vector<Real>>(
43 : "custom_circular_tolerance",
44 : "custom_circular_tolerance>0",
45 : "Custom tolerances (L2 norm) used to verify whether the provided boundaries are circular "
46 : "(default is 1.0e-6).");
47 :
48 14381 : params.addClassDescription(
49 : "This CircularBoundaryCorrectionGenerator object is designed to correct full "
50 : "or partial circular boundaries in a 2D mesh to preserve areas.");
51 :
52 14381 : return params;
53 0 : }
54 :
55 62 : CircularBoundaryCorrectionGenerator::CircularBoundaryCorrectionGenerator(
56 62 : const InputParameters & parameters)
57 : : MeshGenerator(parameters),
58 62 : _input_name(getParam<MeshGeneratorName>("input")),
59 62 : _input_mesh_circular_boundaries(
60 : getParam<std::vector<BoundaryName>>("input_mesh_circular_boundaries")),
61 240 : _transition_layer_ratios(isParamValid("transition_layer_ratios")
62 62 : ? getParam<std::vector<Real>>("transition_layer_ratios")
63 116 : : std::vector<Real>(_input_mesh_circular_boundaries.size(), 0.01)),
64 222 : _custom_circular_tolerance(
65 124 : isParamValid("custom_circular_tolerance")
66 124 : ? (getParam<std::vector<Real>>("custom_circular_tolerance").size() == 1
67 44 : ? std::vector<Real>(
68 : _input_mesh_circular_boundaries.size(),
69 102 : getParam<std::vector<Real>>("custom_circular_tolerance").front())
70 : : getParam<std::vector<Real>>("custom_circular_tolerance"))
71 80 : : std::vector<Real>(_input_mesh_circular_boundaries.size(), 1.0e-6)),
72 62 : _move_end_nodes_in_span_direction(getParam<bool>("move_end_nodes_in_span_direction")),
73 186 : _input(getMeshByName(_input_name))
74 : {
75 62 : if (_transition_layer_ratios.size() != _input_mesh_circular_boundaries.size())
76 4 : paramError("transition_layer_ratios",
77 : "this parameter must have the same length as 'input_mesh_circular_boundaries'.");
78 58 : if (_custom_circular_tolerance.size() != _input_mesh_circular_boundaries.size())
79 4 : paramError("custom_circular_tolerance",
80 : "if provided, this parameter must have either a unity length or the same length as "
81 : "'input_mesh_circular_boundaries'.");
82 54 : }
83 :
84 : std::unique_ptr<MeshBase>
85 54 : CircularBoundaryCorrectionGenerator::generate()
86 : {
87 54 : auto input_mesh = dynamic_cast<ReplicatedMesh *>(_input.get());
88 54 : if (!input_mesh)
89 0 : paramError("input", "Input is not a replicated mesh, which is required");
90 :
91 : // Check if the input mesh has the given boundaries
92 130 : for (const auto & bd : _input_mesh_circular_boundaries)
93 : {
94 80 : if (!MooseMeshUtils::hasBoundaryName(*input_mesh, bd))
95 4 : paramError("input_mesh_circular_boundaries",
96 4 : "the boundary '" + bd + "' does not exist in the input mesh.");
97 : }
98 : _input_mesh_circular_bids =
99 50 : MooseMeshUtils::getBoundaryIDs(*input_mesh, _input_mesh_circular_boundaries, false);
100 :
101 50 : auto side_list = input_mesh->get_boundary_info().build_side_list();
102 50 : input_mesh->get_boundary_info().build_node_list_from_side_list();
103 50 : auto node_list = input_mesh->get_boundary_info().build_node_list();
104 :
105 50 : std::vector<std::vector<Point>> input_circ_bds_pts(_input_mesh_circular_bids.size());
106 : // Loop over all the boundary nodes and store the nodes (points) on each of the circular
107 : // boundaries to be corrected
108 1832 : for (unsigned int i = 0; i < node_list.size(); ++i)
109 : {
110 1786 : auto node_list_it = std::find(_input_mesh_circular_bids.begin(),
111 : _input_mesh_circular_bids.end(),
112 1786 : std::get<1>(node_list[i]));
113 1786 : if (node_list_it != _input_mesh_circular_bids.end())
114 : {
115 842 : input_circ_bds_pts[std::distance(_input_mesh_circular_bids.begin(), node_list_it)].push_back(
116 842 : input_mesh->point(std::get<0>(node_list[i])));
117 842 : if (!MooseUtils::absoluteFuzzyEqual(
118 842 : input_circ_bds_pts[std::distance(_input_mesh_circular_bids.begin(), node_list_it)]
119 842 : .back()(2),
120 1684 : 0.0))
121 4 : paramError("input_mesh_circular_boundaries", "Only boundaries in XY plane are supported.");
122 : }
123 : }
124 : std::vector<std::vector<std::pair<Point, Point>>> input_circ_bds_sds(
125 46 : _input_mesh_circular_bids.size());
126 : std::vector<std::vector<std::pair<dof_id_type, dof_id_type>>> input_circ_bds_sds_ids(
127 46 : _input_mesh_circular_bids.size());
128 : // Loop over all the boundary sides and store the sides (node pairs) on each of the circular
129 : // boundaries to be corrected
130 1800 : for (unsigned int i = 0; i < side_list.size(); ++i)
131 : {
132 1754 : auto side_list_it = std::find(_input_mesh_circular_bids.begin(),
133 : _input_mesh_circular_bids.end(),
134 1754 : std::get<2>(side_list[i]));
135 1754 : if (side_list_it != _input_mesh_circular_bids.end())
136 : {
137 : const auto side =
138 814 : input_mesh->elem_ptr(std::get<0>(side_list[i]))->side_ptr(std::get<1>(side_list[i]));
139 : // In case the sideset contains a pair of duplicated sides with different orientations, we
140 : // only store the first side found.
141 1628 : if (std::find(
142 814 : input_circ_bds_sds_ids[std::distance(_input_mesh_circular_bids.begin(), side_list_it)]
143 : .begin(),
144 814 : input_circ_bds_sds_ids[std::distance(_input_mesh_circular_bids.begin(), side_list_it)]
145 : .end(),
146 1628 : std::make_pair(side->node_id(0), side->node_id(1))) ==
147 814 : input_circ_bds_sds_ids[std::distance(_input_mesh_circular_bids.begin(), side_list_it)]
148 2442 : .end() &&
149 1628 : std::find(
150 814 : input_circ_bds_sds_ids[std::distance(_input_mesh_circular_bids.begin(), side_list_it)]
151 : .begin(),
152 814 : input_circ_bds_sds_ids[std::distance(_input_mesh_circular_bids.begin(), side_list_it)]
153 : .end(),
154 1628 : std::make_pair(side->node_id(1), side->node_id(0))) ==
155 814 : input_circ_bds_sds_ids[std::distance(_input_mesh_circular_bids.begin(), side_list_it)]
156 1628 : .end())
157 : {
158 814 : input_circ_bds_sds[std::distance(_input_mesh_circular_bids.begin(), side_list_it)]
159 814 : .push_back(std::make_pair(side->point(0), side->point(1)));
160 814 : input_circ_bds_sds_ids[std::distance(_input_mesh_circular_bids.begin(), side_list_it)]
161 814 : .push_back(std::make_pair(side->node_id(0), side->node_id(1)));
162 : }
163 814 : }
164 : }
165 :
166 : // Computes the correction coefficient and apply it;
167 : // Also records the nodes that have been modified to avert double-correction
168 46 : std::set<dof_id_type> mod_node_list;
169 : // And check if any of the boundaries are partial circles
170 46 : bool has_partial_circle = false;
171 98 : for (unsigned int i = 0; i < input_circ_bds_pts.size(); i++)
172 : {
173 64 : auto & input_circ_bd_pts = input_circ_bds_pts[i];
174 64 : if (input_circ_bd_pts.size() < 3)
175 4 : paramError("input_mesh_circular_boundaries",
176 : "too few nodes in one of the provided boundaries.");
177 : Real bdry_rad;
178 : const Point boundary_origin =
179 60 : circularCenterCalculator(input_circ_bd_pts, bdry_rad, _custom_circular_tolerance[i]);
180 :
181 : // Check if the boundary is a full circle or partial
182 : // Also make an ordered array of nodes
183 : Real dummy_max_rad;
184 56 : std::vector<dof_id_type> ordered_node_list;
185 : bool is_bdry_closed;
186 112 : FillBetweenPointVectorsTools::isClosedLoop(*input_mesh,
187 : dummy_max_rad,
188 : ordered_node_list,
189 56 : input_circ_bds_sds_ids[i],
190 : boundary_origin,
191 : "boundary",
192 : is_bdry_closed,
193 : true);
194 :
195 56 : has_partial_circle = has_partial_circle || !is_bdry_closed;
196 : // If the user selects to move the end nodes of a partial circular boundary, we need to
197 : // calculate the displacement of the end nodes, which differs from the displacement of the
198 : // other nodes
199 : Real end_node_disp;
200 : // c_coeff is the coefficient used to scale the azimuthal angles
201 : // corr_factor > 1 means the partial boundary is less than a half circle
202 : // corr_factor < 1 means the partial boundary is more than a half circle
203 : // corr_factor = 1 means the boundary is a full circle or a half circle
204 : Real c_coeff;
205 56 : const Real corr_factor = generateRadialCorrectionFactor(input_circ_bds_sds[i],
206 : boundary_origin,
207 : is_bdry_closed,
208 56 : _move_end_nodes_in_span_direction,
209 : c_coeff,
210 : end_node_disp);
211 : // Radial range that within which the nodes will be moved
212 56 : const Real transition_layer_thickness = _transition_layer_ratios[i] * bdry_rad;
213 : // For a partial boundary, we take out the start and end nodes from the ordered node list
214 56 : const Point pt_start = input_mesh->point(ordered_node_list.front());
215 56 : const Point pt_end = input_mesh->point(ordered_node_list.back());
216 : // The direction of the span of the partial boundary (which is an arc)
217 56 : const Point span_direction = is_bdry_closed ? Point(0.0, 0.0, 0.0) : (pt_start - pt_end).unit();
218 : // Although these are also calculated for a full circular boundary, they are not used
219 :
220 : // Find the center of the circular boundary's azimuthal angle and the half of the azimuthal
221 : // range
222 : Real center_bdry_azi;
223 : Real half_azi_range;
224 56 : if (MooseUtils::absoluteFuzzyEqual(c_coeff, 1.0))
225 : {
226 : // For a half circle, the center boundary azimuthal angle is the starting node's azimuthal
227 : // plus or minus pi/2
228 46 : const Point pt_second = input_mesh->point(ordered_node_list[1]);
229 : const Real pt_start_azi =
230 46 : std::atan2(pt_start(1) - boundary_origin(1), pt_start(0) - boundary_origin(0));
231 : const Real pt_second_azi =
232 46 : std::atan2(pt_second(1) - boundary_origin(1), pt_second(0) - boundary_origin(0));
233 : // Figure out whether it is plus or minus pi/2
234 46 : const Real azi_direction =
235 46 : (pt_second_azi - pt_start_azi > 0 or pt_second_azi - pt_start_azi < -M_PI) ? 1.0 : -1.0;
236 46 : center_bdry_azi = pt_start_azi + azi_direction * M_PI / 2.0;
237 46 : half_azi_range = M_PI / 2.0;
238 : }
239 : else
240 : {
241 : const Point center_bdry =
242 10 : (c_coeff > 1.0 ? 1.0 : -1.0) * (pt_start + pt_end) / 2.0 - boundary_origin;
243 : const Real pt_start_azi =
244 10 : std::atan2(pt_start(1) - boundary_origin(1), pt_start(0) - boundary_origin(0));
245 10 : center_bdry_azi = std::atan2(center_bdry(1), center_bdry(0));
246 20 : half_azi_range = std::abs(pt_start_azi - center_bdry_azi > M_PI
247 10 : ? pt_start_azi - center_bdry_azi - 2 * M_PI
248 0 : : (pt_start_azi - center_bdry_azi < -M_PI
249 0 : ? pt_start_azi - center_bdry_azi + 2 * M_PI
250 : : pt_start_azi - center_bdry_azi));
251 : }
252 : // Loop over all the nodes in the mesh to move the nodes in the transition layer so that the
253 : // circular boundary is corrected
254 1936 : for (auto & node : input_mesh->node_ptr_range())
255 : {
256 : // If the end nodes are set to move in the span direction, we should do so
257 1902 : if (node->id() == ordered_node_list.front() && _move_end_nodes_in_span_direction &&
258 18 : end_node_disp >= 0.0)
259 : {
260 10 : if (!mod_node_list.emplace((*node).id()).second)
261 0 : paramError("transition_layer_ratios", "the transition layers are overlapped.");
262 10 : (*node) = (*node) + end_node_disp * bdry_rad * span_direction;
263 1218 : continue;
264 : }
265 1892 : if (node->id() == ordered_node_list.back() && _move_end_nodes_in_span_direction &&
266 18 : end_node_disp >= 0.0)
267 : {
268 10 : if (!mod_node_list.emplace((*node).id()).second)
269 0 : paramError("transition_layer_ratios", "the transition layers are overlapped.");
270 10 : (*node) = (*node) - end_node_disp * bdry_rad * span_direction;
271 10 : continue;
272 : }
273 1864 : const Point tmp_pt = (*node) - boundary_origin;
274 1864 : const Real tmp_pt_azi = std::atan2(tmp_pt(1), tmp_pt(0));
275 1864 : const Real angle_diff =
276 1864 : tmp_pt_azi - center_bdry_azi > M_PI
277 3496 : ? tmp_pt_azi - center_bdry_azi - 2.0 * M_PI
278 1632 : : (tmp_pt_azi - center_bdry_azi < -M_PI ? tmp_pt_azi - center_bdry_azi + 2.0 * M_PI
279 : : tmp_pt_azi - center_bdry_azi);
280 1864 : const Real pt_rad = tmp_pt.norm();
281 : // Skip the node if it is not within the radial and azimuthal range
282 900 : if ((!is_bdry_closed && MooseUtils::absoluteFuzzyGreaterThan(
283 900 : std::abs(angle_diff), half_azi_range, libMesh::TOLERANCE)) ||
284 4390 : pt_rad < bdry_rad - transition_layer_thickness ||
285 1626 : pt_rad > bdry_rad + transition_layer_thickness)
286 1198 : continue;
287 666 : const Real tmp_pt_azi_new =
288 666 : (_move_end_nodes_in_span_direction ? c_coeff : 1.0) * angle_diff + center_bdry_azi;
289 666 : const Real mod_corr_factor =
290 666 : pt_rad <= bdry_rad
291 666 : ? (1.0 + (corr_factor - 1.0) * (pt_rad - bdry_rad + transition_layer_thickness) /
292 : transition_layer_thickness)
293 134 : : (1.0 + (corr_factor - 1.0) * (bdry_rad + transition_layer_thickness - pt_rad) /
294 : transition_layer_thickness);
295 666 : if (!MooseUtils::absoluteFuzzyEqual(mod_corr_factor, 1.0))
296 : {
297 666 : if (!mod_node_list.emplace((*node).id()).second)
298 4 : paramError("transition_layer_ratios", "the transition layers are overlapped.");
299 : }
300 : const Point tmp_pt_alt =
301 662 : Point(pt_rad * std::cos(tmp_pt_azi_new), pt_rad * std::sin(tmp_pt_azi_new), tmp_pt(2)) *
302 662 : mod_corr_factor;
303 662 : (*node) = tmp_pt_alt + boundary_origin;
304 52 : }
305 52 : }
306 34 : if (!has_partial_circle && _move_end_nodes_in_span_direction)
307 4 : paramError("move_end_nodes_in_span_direction",
308 : "all the boundaries are closed, this parameter should be be set as 'true'.");
309 : // Loop over all the elements to check if any of them are inverted due to the correction
310 1400 : for (auto & elem : input_mesh->element_ptr_range())
311 1370 : if (elem->volume() < 0.0)
312 0 : paramError("transition_layer_ratios",
313 : "some elements are inverted during circular boundary correction, consider tuning "
314 30 : "this parameter.");
315 :
316 60 : return dynamic_pointer_cast<MeshBase>(_input);
317 30 : }
318 :
319 : Point
320 60 : CircularBoundaryCorrectionGenerator::circularCenterCalculator(const std::vector<Point> & pts_list,
321 : Real & radius,
322 : const Real tol) const
323 : {
324 : // Usually, just using the first three points would work
325 : // Here a more complex selection is made in case the first three points are too close to each
326 : // other.
327 : // Step 1: find the farthest point from the first point
328 60 : Real tmp_dis(0.0);
329 60 : unsigned int farthest_pt_id(0);
330 830 : for (unsigned int i = 1; i < pts_list.size(); i++)
331 770 : if ((pts_list[i] - pts_list[0]).norm() > tmp_dis)
332 : {
333 458 : tmp_dis = (pts_list[i] - pts_list[0]).norm();
334 458 : farthest_pt_id = i;
335 : }
336 : // Step 2: find the third point that is nearly in the middle of the first two
337 60 : unsigned int mid_pt_id(0);
338 830 : for (unsigned int i = 1; i < pts_list.size(); i++)
339 1480 : if (i != farthest_pt_id && std::abs((pts_list[farthest_pt_id] - pts_list[i]).norm() -
340 1480 : (pts_list[0] - pts_list[i]).norm()) < tmp_dis)
341 : {
342 244 : tmp_dis = std::abs((pts_list[farthest_pt_id] - pts_list[i]).norm() -
343 244 : (pts_list[0] - pts_list[i]).norm());
344 244 : mid_pt_id = i;
345 : }
346 :
347 60 : const Real x12 = pts_list[0](0) - pts_list[farthest_pt_id](0);
348 60 : const Real x13 = pts_list[0](0) - pts_list[mid_pt_id](0);
349 :
350 60 : const Real y12 = pts_list[0](1) - pts_list[farthest_pt_id](1);
351 60 : const Real y13 = pts_list[0](1) - pts_list[mid_pt_id](1);
352 :
353 : const Real sx13 =
354 60 : pts_list[0](0) * pts_list[0](0) - pts_list[mid_pt_id](0) * pts_list[mid_pt_id](0);
355 : const Real sy13 =
356 60 : pts_list[0](1) * pts_list[0](1) - pts_list[mid_pt_id](1) * pts_list[mid_pt_id](1);
357 :
358 : const Real sx21 =
359 60 : pts_list[farthest_pt_id](0) * pts_list[farthest_pt_id](0) - pts_list[0](0) * pts_list[0](0);
360 : const Real sy21 =
361 60 : pts_list[farthest_pt_id](1) * pts_list[farthest_pt_id](1) - pts_list[0](1) * pts_list[0](1);
362 :
363 60 : const Real f =
364 60 : (sx13 * x12 + sy13 * x12 + sx21 * x13 + sy21 * x13) / (2.0 * (-y13 * x12 + y12 * x13));
365 60 : const Real g =
366 60 : (sx13 * y12 + sy13 * y12 + sx21 * y13 + sy21 * y13) / (2.0 * (-x13 * y12 + x12 * y13));
367 :
368 60 : const Real c = -pts_list[0](0) * pts_list[0](0) - pts_list[0](1) * pts_list[0](1) -
369 60 : 2.0 * g * pts_list[0](0) - 2.0 * f * pts_list[0](1);
370 60 : const Real r2 = f * f + g * g - c;
371 :
372 60 : radius = std::sqrt(r2);
373 60 : const Point circ_origin = Point(-g, -f, 0.0);
374 :
375 734 : for (const auto & pt : pts_list)
376 678 : if (!MooseUtils::absoluteFuzzyEqual((pt - circ_origin).norm(), radius, tol))
377 4 : paramError("input_mesh_circular_boundaries",
378 4 : "the boundary provided is not circular. The generator found point " +
379 8 : Moose::stringify(circ_origin) +
380 4 : " to be the likely circle origin to the boundary, with a radius of " +
381 8 : std::to_string(radius) + ". However, boundary node " + Moose::stringify(pt) +
382 8 : " is at a radius of " + std::to_string((pt - circ_origin).norm()));
383 56 : return circ_origin;
384 : }
385 :
386 : Real
387 56 : CircularBoundaryCorrectionGenerator::generateRadialCorrectionFactor(
388 : const std::vector<std::pair<Point, Point>> & bd_side_list,
389 : const Point & circle_center,
390 : const bool is_closed_loop,
391 : const bool move_end_nodes_in_span_direction,
392 : Real & c_coeff,
393 : Real & end_node_disp) const
394 : {
395 56 : if (bd_side_list.empty())
396 0 : mooseError("The 'bd_side_list' argument of "
397 : "CircularBoundaryCorrectionGenerator::generateRadialCorrectionFactor is empty.");
398 56 : Real acc_area(0.0);
399 56 : Real acc_azi(0.0);
400 56 : std::vector<Real> d_azi_list;
401 706 : for (const auto & bd_side : bd_side_list)
402 : {
403 650 : const Point v1 = bd_side.first - circle_center;
404 650 : const Point v2 = bd_side.second - circle_center;
405 : // Use cross to calculate r * r * sin(d_azi_i) and then normalized by r * r to get sin(d_azi_i)
406 650 : acc_area += v1.cross(v2).norm() / v1.norm() / v2.norm();
407 650 : const Real azi1 = std::atan2(v1(1), v1(0));
408 650 : const Real azi2 = std::atan2(v2(1), v2(0));
409 : const Real d_azi =
410 650 : std::abs(azi1 - azi2) > M_PI ? 2 * M_PI - std::abs(azi1 - azi2) : std::abs(azi1 - azi2);
411 650 : acc_azi += d_azi;
412 650 : d_azi_list.push_back(d_azi);
413 : }
414 :
415 56 : if (!MooseUtils::absoluteFuzzyEqual(acc_azi, M_PI) && !is_closed_loop &&
416 : move_end_nodes_in_span_direction)
417 : {
418 10 : const Real k_1 = 1.0 + std::cos(acc_azi);
419 10 : const Real k_2 = acc_azi - std::sin(acc_azi);
420 :
421 10 : unsigned int ct = 0;
422 10 : Real diff = 1.0;
423 10 : Real c_old = 1.0;
424 10 : Real c_new = 1.0;
425 30 : while (diff >= libMesh::TOLERANCE && ct < 100)
426 : {
427 20 : c_old = c_new;
428 20 : const Real gc = k_2 + k_2 * std::cos(acc_azi * c_old) + k_1 * std::sin(acc_azi * c_old) -
429 20 : k_1 * sineSummation(d_azi_list, c_old);
430 20 : const Real gcp = -k_2 * acc_azi * std::sin(acc_azi * c_old) +
431 20 : k_1 * acc_azi * std::cos(acc_azi * c_old) -
432 20 : k_1 * sinePrimeSummation(d_azi_list, c_old);
433 20 : c_new = c_old - gc / gcp;
434 20 : diff = std::abs(c_new - c_old);
435 20 : ct++;
436 : }
437 :
438 10 : if (ct >= 100)
439 0 : mooseError("Newton-Raphson method did not converge for generating the radial correction "
440 : "factor for circular area preservation.");
441 :
442 10 : const Real norm_corr_factor = std::cos(0.5 * acc_azi) / std::cos(0.5 * acc_azi * c_new);
443 10 : end_node_disp = norm_corr_factor * std::sin(0.5 * acc_azi * c_new) - std::sin(0.5 * acc_azi);
444 10 : c_coeff = c_new;
445 :
446 10 : return norm_corr_factor;
447 : }
448 46 : end_node_disp = -1.0;
449 46 : c_coeff = 1.0;
450 46 : return std::sqrt(acc_azi / acc_area);
451 56 : }
452 :
453 : Real
454 20 : CircularBoundaryCorrectionGenerator::sineSummation(const std::vector<Real> & th_list,
455 : const Real coeff) const
456 : {
457 20 : Real acc(0.0);
458 220 : for (const auto & th : th_list)
459 200 : acc += std::sin(th * coeff);
460 20 : return acc;
461 : }
462 :
463 : Real
464 20 : CircularBoundaryCorrectionGenerator::sinePrimeSummation(const std::vector<Real> & th_list,
465 : const Real coeff) const
466 : {
467 : // Note that this calculates the derivative with regards to coeff
468 20 : Real acc(0.0);
469 220 : for (const auto & th : th_list)
470 200 : acc += th * std::cos(th * coeff);
471 20 : return acc;
472 : }
|
