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 "SplineUtils.h"
11 : #include "MooseError.h"
12 : #include "libMeshReducedNamespace.h"
13 : #include "MooseUtils.h"
14 : #include <algorithm>
15 :
16 : namespace SplineUtils
17 : {
18 : std::vector<Point>
19 52 : circularControlPoints(const libMesh::Point & start_point,
20 : const libMesh::Point & end_point,
21 : const libMesh::RealVectorValue & parallel_direction,
22 : const unsigned int num_cps)
23 : {
24 : // establish the direction vector between the two points
25 52 : const auto orthogonal_direction = start_point - end_point;
26 :
27 : // circular control points will only honor the derivatives if the points are colinear along the
28 : // orthogonal direction
29 52 : if (!(MooseUtils::absoluteFuzzyEqual(parallel_direction.unit() * orthogonal_direction.unit(), 0)))
30 8 : mooseError("Parallel lines cannot be interpolated using points along a circle! Direction "
31 : "vectors are incompatible with point locations.\nStart point: ",
32 : start_point,
33 : "\nEnd point: ",
34 : end_point,
35 : "\nDirection: ",
36 : parallel_direction);
37 :
38 : // define the distance between the two points
39 44 : const auto distance_between = orthogonal_direction.norm();
40 :
41 : // connecting circle is assumed to be between each point along the orthogonal_direction vector
42 44 : const auto radius = distance_between / 2.0;
43 :
44 : // normalize vectors
45 44 : const auto unit_normal = orthogonal_direction / distance_between;
46 44 : const auto unit_parallel = parallel_direction / parallel_direction.norm();
47 :
48 : // calculate circle center
49 44 : const auto center_point = end_point + radius * unit_normal;
50 :
51 : // initialize vector of control points
52 44 : std::vector<Point> control_points;
53 :
54 : // loop over the number of control points to generate the control points
55 : mooseAssert(num_cps > 1, "This routine does not support a single control point");
56 1064 : for (const auto i : make_range(num_cps))
57 : {
58 1020 : const auto t = (Real)i / (Real)(num_cps - 1);
59 2040 : control_points.push_back(center_point + radius * (std::cos(t * M_PI) * unit_normal +
60 3060 : std::sin(t * M_PI) * unit_parallel));
61 : }
62 :
63 88 : return control_points;
64 0 : }
65 :
66 : std::vector<Point>
67 562 : bSplineControlPoints(const libMesh::Point & start_point,
68 : const libMesh::Point & end_point,
69 : const libMesh::RealVectorValue & start_direction,
70 : const libMesh::RealVectorValue & end_direction,
71 : const unsigned int cps_per_half,
72 : const libMesh::Real sharpness)
73 : {
74 : // check that start_point is different from end_point
75 562 : if ((start_point - end_point).norm_sq() < libMesh::TOLERANCE * libMesh::TOLERANCE)
76 2 : mooseError("Start and end points must be different.");
77 : // check that neither direction is the zero vector
78 1116 : if (start_direction.norm_sq() < libMesh::TOLERANCE * libMesh::TOLERANCE ||
79 556 : end_direction.norm_sq() < libMesh::TOLERANCE * libMesh::TOLERANCE)
80 6 : mooseError("Direction vectors must have a nonzero magnitude.");
81 : mooseAssert(sharpness <= 1.0 && sharpness >= 0, "Sharpness must be in [0,1]!");
82 :
83 : // check if directions are parallel, use a special case if so
84 554 : const bool parallel = ((start_direction.cross(end_direction)).norm() < libMesh::TOLERANCE);
85 554 : if (parallel)
86 : {
87 : // We should not revert to a circle for this case
88 48 : if (MooseUtils::absoluteFuzzyEqual(start_direction.unit(), -end_direction.unit()))
89 0 : mooseError("Start and end directions are opposite, this is not currently supported");
90 48 : mooseWarning("Directions are parallel! Attempting to use circular control points...");
91 48 : unsigned int num_cps = 2 * cps_per_half + 1;
92 48 : if (num_cps < 30)
93 48 : mooseWarning("Number of control points required for circular control points is much greater "
94 96 : "than for BSplines. Current num_cps: " +
95 96 : std::to_string(num_cps));
96 48 : return SplineUtils::circularControlPoints(start_point, end_point, start_direction, num_cps);
97 : }
98 :
99 : // find closest points --> these will be identical if the extrapolated lines intersect
100 : const auto closest_points_vector =
101 506 : SplineUtils::closestPoints(start_point, end_point, start_direction, end_direction);
102 506 : const auto & closest_point_1 = closest_points_vector.first;
103 506 : const auto & closest_point_2 = closest_points_vector.second;
104 :
105 : std::vector<Point> first_half = SplineUtils::controlPointsAlongLine(
106 506 : start_point, closest_point_1, start_direction, sharpness, cps_per_half);
107 : std::vector<Point> second_half = SplineUtils::controlPointsAlongLine(
108 506 : end_point, closest_point_2, end_direction, sharpness, cps_per_half, true);
109 :
110 506 : std::vector<Point> control_points;
111 :
112 : // put it all together
113 2000 : for (const auto i : index_range(first_half))
114 1494 : control_points.push_back(first_half[i]);
115 2000 : for (const auto i : index_range(second_half))
116 1494 : control_points.push_back(second_half[i]);
117 :
118 506 : return control_points;
119 506 : }
120 :
121 : std::vector<Point>
122 1016 : controlPointsAlongLine(const libMesh::Point & start_point,
123 : const libMesh::Point & end_point,
124 : const libMesh::RealVectorValue & direction_vector,
125 : const libMesh::Real sharpness,
126 : const unsigned int num_cps,
127 : const bool reverse_order)
128 : {
129 : // initialize a vector of control points with the starting point inserted
130 1016 : std::vector<Point> control_points;
131 :
132 : mooseAssert(num_cps != 0, "Number of control points cannot be zero!");
133 1016 : if (num_cps == 1)
134 2 : control_points.push_back(start_point);
135 :
136 : // create the control points by varying the sharpness linearly
137 : else
138 : {
139 4008 : for (const auto i : libMesh::make_range(num_cps))
140 : {
141 2994 : const auto point_sharpness = (Real)i / (Real)(num_cps - 1) * sharpness;
142 2994 : control_points.push_back(
143 5988 : SplineUtils::makeControlPoint(start_point, end_point, direction_vector, point_sharpness));
144 : }
145 :
146 : // reverse_order parameter only set to true when required
147 1014 : if (reverse_order)
148 506 : std::reverse(control_points.begin(), control_points.end());
149 : }
150 1016 : return control_points;
151 0 : }
152 :
153 : libMesh::Point
154 2994 : makeControlPoint(const libMesh::Point & start_point,
155 : const libMesh::Point & end_point,
156 : const libMesh::RealVectorValue & direction_vector,
157 : const libMesh::Real sharpness)
158 : {
159 : // check that sharpness value is ok
160 : mooseAssert(sharpness <= 1.0 && sharpness >= 0, "Sharpness must be in [0,1]!");
161 :
162 : // normalize direction vector
163 2994 : const auto unit_direction = direction_vector.unit();
164 :
165 : // calculate the distance between points
166 2994 : const auto distance = (end_point - start_point).norm();
167 2994 : const auto distance_sq = distance * distance;
168 :
169 : // initialize return quantity and attempt to calculate control vector
170 2994 : libMesh::Point control_point = start_point + sharpness * distance * unit_direction;
171 :
172 : // check if we went the wrong way, further from the end_point
173 2994 : if ((end_point - control_point).norm_sq() > distance_sq)
174 : {
175 0 : control_point = start_point - sharpness * distance * unit_direction;
176 : mooseAssert((end_point - control_point).norm_sq() <= distance_sq,
177 : "Also going further from end point");
178 : }
179 :
180 5988 : return control_point;
181 : }
182 :
183 : std::pair<Point, Point>
184 512 : closestPoints(const libMesh::Point & point_1,
185 : const libMesh::Point & point_2,
186 : const libMesh::RealVectorValue & direction_1,
187 : const libMesh::RealVectorValue & direction_2)
188 : {
189 512 : const auto n_vec = direction_1.cross(direction_2);
190 :
191 : // check if n_vec is the zero vector (indicates parallel lines)
192 512 : if (n_vec.norm_sq() < libMesh::TOLERANCE * libMesh::TOLERANCE)
193 0 : mooseError("Lines are parallel! Infinitely many closest points exist.");
194 :
195 512 : const auto n1_vec = direction_1.cross(n_vec);
196 512 : const auto n2_vec = direction_2.cross(n_vec);
197 :
198 : // calculate closest points
199 : // take absolute value to ensure specified directions are honored
200 : libMesh::Real point_a_dir_scalar =
201 512 : std::abs((point_2 - point_1) * n2_vec / (direction_1 * n2_vec));
202 : libMesh::Real point_b_dir_scalar =
203 512 : std::abs((point_1 - point_2) * n1_vec / (direction_2 * n1_vec));
204 512 : libMesh::Point point_a = point_1 + point_a_dir_scalar * direction_1;
205 512 : libMesh::Point point_b = point_2 + point_b_dir_scalar * direction_2;
206 :
207 : // points will be returned as a vector containing two points
208 512 : const auto closest_points = std::make_pair(point_a, point_b);
209 :
210 1024 : return closest_points;
211 : }
212 : }
|
