sphere.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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
37  _rad(-1.)
38 {
39 }
40 
41 
42 
44  const Real r)
45 {
46  libmesh_assert_greater (r, 0.);
47 
48  this->create_from_center_radius (c, r);
49 }
50 
51 
52 
53 Sphere::Sphere (const Sphere& other_sphere) :
54  Surface()
55 {
56  this->create_from_center_radius (other_sphere.center(),
57  other_sphere.radius());
58 }
59 
60 
61 
63  const Point& pb,
64  const Point& pc,
65  const Point& pd)
66 {
67  Point pad = pa - pd;
68  Point pbd = pb - pd;
69  Point pcd = pc - pd;
70 
71  TensorValue<Real> T(pad,pbd,pcd);
72 
73  Real D = T.det();
74 
75  // The points had better not be coplanar
76  libmesh_assert_greater (std::abs(D), 1e-12);
77 
78  Real e = 0.5*(pa.size_sq() - pd.size_sq());
79  Real f = 0.5*(pb.size_sq() - pd.size_sq());
80  Real g = 0.5*(pc.size_sq() - pd.size_sq());
81 
82  TensorValue<Real> T1(e,pad(1),pad(2),
83  f,pbd(1),pbd(2),
84  g,pcd(1),pcd(2));
85  Real sx = T1.det()/D;
86 
87  TensorValue<Real> T2(pad(0),e,pad(2),
88  pbd(0),f,pbd(2),
89  pcd(0),g,pcd(2));
90  Real sy = T2.det()/D;
91 
92  TensorValue<Real> T3(pad(0),pad(1),e,
93  pbd(0),pbd(1),f,
94  pcd(0),pcd(1),g);
95  Real sz = T3.det()/D;
96 
97  Point c(sx,sy,sz);
98  Real r = (c-pa).size();
99 
100  this->create_from_center_radius(c,r);
101 }
102 
103 
104 
106 {
107 }
108 
109 
110 
112 {
113  this->center() = c;
114  this->radius() = r;
115 
116  libmesh_assert_greater (this->radius(), 0.);
117 }
118 
119 
120 
121 bool Sphere::intersects (const Sphere& other_sphere) const
122 {
123  return distance(other_sphere) < 0 ? true : false;
124 }
125 
126 
127 
128 Real Sphere::distance (const Sphere& other_sphere) const
129 {
130  libmesh_assert_greater ( this->radius(), 0. );
131  libmesh_assert_greater ( other_sphere.radius(), 0. );
132 
133  const Real the_distance = (this->center() - other_sphere.center()).size();
134 
135  return the_distance - (this->radius() + other_sphere.radius());
136 }
137 
138 
139 
140 bool Sphere::above_surface (const Point& p) const
141 {
142  libmesh_assert_greater (this->radius(), 0.);
143 
144  // create a vector from the center to the point.
145  const Point w = p - this->center();
146 
147  if (w.size() > this->radius())
148  return true;
149 
150  return false;
151 }
152 
153 
154 
155 bool Sphere::below_surface (const Point& p) const
156 {
157  libmesh_assert_greater (this->radius(), 0.);
158 
159  return ( !this->above_surface (p) );
160 }
161 
162 
163 
164 bool Sphere::on_surface (const Point& p) const
165 {
166  libmesh_assert_greater (this->radius(), 0.);
167 
168  // Create a vector from the center to the point.
169  const Point w = p - this->center();
170 
171  // if the size of that vector is the same as the radius() then
172  // the point is on the surface.
173  if (std::abs(w.size() - this->radius()) < 1.e-10)
174  return true;
175 
176  return false;
177 }
178 
179 
180 
182 {
183  libmesh_assert_greater (this->radius(), 0.);
184 
185  // get the normal from the surface in the direction
186  // of p
187  Point normal = this->unit_normal (p);
188 
189  // The closest point on the sphere is in the direction
190  // of the normal a distance r from the center.
191  const Point cp = this->center() + normal*this->radius();
192 
193  return cp;
194 }
195 
196 
197 
199 {
200  libmesh_assert_greater (this->radius(), 0.);
201 
202  libmesh_assert_not_equal_to (p, this->center());
203 
204  // Create a vector from the center to the point
205  Point n = p - this->center();
206 
207  return n.unit();
208 }
209 
210 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:07 UTC

Hosted By:
SourceForge.net Logo