Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library 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 GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // C++ includes
21 : #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
22 : #include <cmath> // for std::abs
23 :
24 :
25 : // Local includes
26 : #include "libmesh/tensor_value.h"
27 : #include "libmesh/sphere.h"
28 :
29 : namespace libMesh
30 : {
31 :
32 :
33 :
34 : // ------------------------------------------------------------
35 : // Sphere class member functions
36 0 : Sphere::Sphere () :
37 0 : _rad(-1.)
38 : {
39 0 : }
40 :
41 :
42 :
43 1074 : Sphere::Sphere (const Point & c,
44 1074 : const Real r)
45 : {
46 34 : libmesh_assert_greater (r, 0.);
47 :
48 1074 : this->create_from_center_radius (c, r);
49 1074 : }
50 :
51 :
52 :
53 0 : Sphere::Sphere (const Sphere & other_sphere) :
54 0 : Surface()
55 : {
56 0 : this->create_from_center_radius (other_sphere.center(),
57 : other_sphere.radius());
58 0 : }
59 :
60 :
61 :
62 0 : Sphere::Sphere(const Point & pa,
63 : const Point & pb,
64 : const Point & pc,
65 0 : const Point & pd)
66 : {
67 : #if LIBMESH_DIM > 2
68 0 : Point pad = pa - pd;
69 0 : Point pbd = pb - pd;
70 0 : Point pcd = pc - pd;
71 :
72 0 : TensorValue<Real> T(pad,pbd,pcd);
73 :
74 0 : Real D = T.det();
75 :
76 : // The points had better not be coplanar
77 0 : libmesh_assert_greater (std::abs(D), 1e-12);
78 :
79 0 : Real e = 0.5*(pa.norm_sq() - pd.norm_sq());
80 0 : Real f = 0.5*(pb.norm_sq() - pd.norm_sq());
81 0 : Real g = 0.5*(pc.norm_sq() - pd.norm_sq());
82 :
83 0 : TensorValue<Real> T1(e,pad(1),pad(2),
84 0 : f,pbd(1),pbd(2),
85 0 : g,pcd(1),pcd(2));
86 0 : Real sx = T1.det()/D;
87 :
88 0 : TensorValue<Real> T2(pad(0),e,pad(2),
89 0 : pbd(0),f,pbd(2),
90 0 : pcd(0),g,pcd(2));
91 0 : Real sy = T2.det()/D;
92 :
93 0 : TensorValue<Real> T3(pad(0),pad(1),e,
94 0 : pbd(0),pbd(1),f,
95 0 : pcd(0),pcd(1),g);
96 0 : Real sz = T3.det()/D;
97 :
98 0 : Point c(sx,sy,sz);
99 0 : Real r = (c-pa).norm();
100 :
101 0 : this->create_from_center_radius(c,r);
102 : #else // LIBMESH_DIM > 2
103 : libmesh_ignore(pa, pb, pc, pd);
104 : libmesh_not_implemented();
105 : #endif
106 0 : }
107 :
108 :
109 :
110 1006 : Sphere::~Sphere () = default;
111 :
112 :
113 :
114 1074 : void Sphere::create_from_center_radius (const Point & c, const Real r)
115 : {
116 1074 : this->center() = c;
117 1074 : this->radius() = r;
118 :
119 34 : libmesh_assert_greater (this->radius(), 0.);
120 1074 : }
121 :
122 :
123 :
124 0 : bool Sphere::intersects (const Sphere & other_sphere) const
125 : {
126 0 : return distance(other_sphere) < 0 ? true : false;
127 : }
128 :
129 :
130 :
131 0 : Real Sphere::distance (const Sphere & other_sphere) const
132 : {
133 0 : libmesh_assert_greater ( this->radius(), 0. );
134 0 : libmesh_assert_greater ( other_sphere.radius(), 0. );
135 :
136 0 : const Real the_distance = (this->center() - other_sphere.center()).norm();
137 :
138 0 : return the_distance - (this->radius() + other_sphere.radius());
139 : }
140 :
141 :
142 :
143 0 : bool Sphere::above_surface (const Point & p) const
144 : {
145 0 : libmesh_assert_greater (this->radius(), 0.);
146 :
147 : // create a vector from the center to the point.
148 0 : const Point w = p - this->center();
149 :
150 0 : if (w.norm() > this->radius())
151 0 : return true;
152 :
153 0 : return false;
154 : }
155 :
156 :
157 :
158 0 : bool Sphere::below_surface (const Point & p) const
159 : {
160 0 : libmesh_assert_greater (this->radius(), 0.);
161 :
162 0 : return ( !this->above_surface (p) );
163 : }
164 :
165 :
166 :
167 0 : bool Sphere::on_surface (const Point & p) const
168 : {
169 0 : libmesh_assert_greater (this->radius(), 0.);
170 :
171 : // Create a vector from the center to the point.
172 0 : const Point w = p - this->center();
173 :
174 : // if the size of that vector is the same as the radius() then
175 : // the point is on the surface.
176 0 : if (std::abs(w.norm() - this->radius()) < 1.e-10)
177 0 : return true;
178 :
179 0 : return false;
180 : }
181 :
182 :
183 :
184 259586 : Point Sphere::closest_point (const Point & p) const
185 : {
186 11430 : libmesh_assert_greater (this->radius(), 0.);
187 :
188 : // get the normal from the surface in the direction
189 : // of p
190 259586 : Point normal = this->unit_normal (p);
191 :
192 : // The closest point on the sphere is in the direction
193 : // of the normal a distance r from the center.
194 22860 : const Point cp = this->center() + normal*this->radius();
195 :
196 271016 : return cp;
197 : }
198 :
199 :
200 :
201 259586 : Point Sphere::unit_normal (const Point & p) const
202 : {
203 11430 : libmesh_assert_greater (this->radius(), 0.);
204 :
205 11430 : libmesh_assert_not_equal_to (p, this->center());
206 :
207 : // Create a vector from the center to the point
208 11430 : Point n = p - this->center();
209 :
210 259586 : return n.unit();
211 : }
212 :
213 : } // namespace libMesh
|
