Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2004 Sandia Corporation and Argonne National
5 : : Laboratory. Under the terms of Contract DE-AC04-94AL85000
6 : : with Sandia Corporation, the U.S. Government retains certain
7 : : rights in this software.
8 : :
9 : : This library is free software; you can redistribute it and/or
10 : : modify it under the terms of the GNU Lesser General Public
11 : : License as published by the Free Software Foundation; either
12 : : version 2.1 of the License, or (at your option) any later version.
13 : :
14 : : This library is distributed in the hope that it will be useful,
15 : : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 : : Lesser General Public License for more details.
18 : :
19 : : You should have received a copy of the GNU Lesser General Public License
20 : : (lgpl.txt) along with this library; if not, write to the Free Software
21 : : Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 : :
23 : : [email protected], [email protected], [email protected],
24 : : [email protected], [email protected], [email protected],
25 : : [email protected]
26 : :
27 : : ***************************************************************** */
28 : : #include "Mesquite.hpp"
29 : : #include "SphericalDomain.hpp"
30 : : #include "Vector3D.hpp"
31 : : #include "MsqError.hpp"
32 : : #include "MsqVertex.hpp"
33 : : #include "DomainUtil.hpp"
34 : : #include "MsqMatrix.hpp"
35 : : #include "moab/Util.hpp"
36 : :
37 : : #ifdef MSQ_HAVE_IEEEFP_H
38 : : #include <ieeefp.h>
39 : : #endif
40 : :
41 : : #include <algorithm>
42 : :
43 [ - + ]: 4 : MBMesquite::SphericalDomain::~SphericalDomain() {}
44 : :
45 : 0 : void MBMesquite::SphericalDomain::snap_to( Mesh::VertexHandle /*entity_handle*/, Vector3D& coordinate ) const
46 : : {
47 : : // Get vector center to coordinate, store in coordinate.
48 : 0 : coordinate -= mCenter;
49 : : // Get distance from center of sphere
50 : 0 : double len = coordinate.length();
51 : : // Scale vector to have length of radius
52 : 0 : coordinate *= mRadius / len;
53 : : // If was at center, return arbitrary position on sphere
54 : : // (all possitions are equally close)
55 [ # # ]: 0 : if( !moab::Util::is_finite( coordinate.x() ) ) coordinate.set( mRadius, 0.0, 0.0 );
56 : : // Get position from vector
57 : 0 : coordinate += mCenter;
58 : 0 : }
59 : :
60 : 103692 : void MBMesquite::SphericalDomain::vertex_normal_at( Mesh::VertexHandle /*entity_handle*/, Vector3D& coordinate ) const
61 : : {
62 : : // normal is vector from center to input position
63 : 103692 : coordinate -= mCenter;
64 : : // make it a unit vector
65 : 103692 : double length = coordinate.length();
66 : 103692 : coordinate /= length;
67 : : // if input position was at center, choose same position
68 : : // on sphere as snap_to.
69 [ - + ]: 103692 : if( !moab::Util::is_finite( coordinate.x() ) ) coordinate.set( 1.0, 0.0, 0.0 );
70 : 103692 : }
71 : 0 : void MBMesquite::SphericalDomain::element_normal_at( Mesh::ElementHandle h, Vector3D& coordinate ) const
72 : : {
73 : 0 : SphericalDomain::vertex_normal_at( h, coordinate );
74 : 0 : }
75 : :
76 : 20742 : void MBMesquite::SphericalDomain::vertex_normal_at( const MBMesquite::Mesh::VertexHandle* handle,
77 : : MBMesquite::Vector3D coords[], unsigned count,
78 : : MBMesquite::MsqError& ) const
79 : : {
80 [ + + ]: 124434 : for( unsigned i = 0; i < count; ++i )
81 : 103692 : vertex_normal_at( handle[i], coords[i] );
82 : 20742 : }
83 : :
84 : 172902 : void MBMesquite::SphericalDomain::closest_point( MBMesquite::Mesh::VertexHandle, const MBMesquite::Vector3D& position,
85 : : MBMesquite::Vector3D& closest, MBMesquite::Vector3D& normal,
86 : : MBMesquite::MsqError& ) const
87 : : {
88 [ + - ]: 172902 : normal = position - mCenter;
89 : 172902 : normal.normalize();
90 [ - + ]: 172902 : if( !moab::Util::is_finite( normal.x() ) ) normal.set( 1.0, 0.0, 0.0 );
91 [ + - ][ + - ]: 172902 : closest = mCenter + mRadius * normal;
92 : 172902 : }
93 : :
94 : 20742 : void MBMesquite::SphericalDomain::domain_DoF( const Mesh::VertexHandle*, unsigned short* dof_array, size_t num_vertices,
95 : : MsqError& ) const
96 : : {
97 [ + - ]: 20742 : std::fill( dof_array, dof_array + num_vertices, 2 );
98 : 20742 : }
99 : :
100 : 0 : void MBMesquite::SphericalDomain::fit_vertices( Mesh* mesh, MsqError& err, double epsilon )
101 : : {
102 [ # # ]: 0 : std::vector< Mesh::VertexHandle > verts;
103 [ # # ]: 0 : mesh->get_all_vertices( verts, err );
104 [ # # ][ # # ]: 0 : if( !MSQ_CHKERR( err ) ) fit_vertices( mesh, arrptr( verts ), verts.size(), err, epsilon );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
105 : 0 : }
106 : :
107 : 0 : void MBMesquite::SphericalDomain::fit_vertices( Mesh* mesh, const Mesh::VertexHandle* verts, size_t num_verts,
108 : : MsqError& err, double epsilon )
109 : : {
110 [ # # ]: 0 : std::vector< MsqVertex > coords( num_verts );
111 [ # # ][ # # ]: 0 : mesh->vertices_get_coordinates( verts, arrptr( coords ), num_verts, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
112 : :
113 [ # # ][ # # ]: 0 : if( epsilon <= 0.0 ) epsilon = DomainUtil::default_tolerance( arrptr( coords ), num_verts );
[ # # ]
114 : :
115 [ # # ][ # # ]: 0 : Vector3D pts[4];
116 [ # # ][ # # ]: 0 : if( !DomainUtil::non_coplanar_vertices( arrptr( coords ), num_verts, pts, epsilon ) )
[ # # ]
117 : : {
118 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "All vertices are co-planar", MsqError::INVALID_MESH );
119 : 0 : return;
120 : : }
121 : :
122 : : // solve deterinant form of four-point sphere
123 : :
124 : : // Define the bottom 4 rows of a 5x5 matrix. The top
125 : : // row contains the variables we are solving for, so just
126 : : // fill it with ones.
127 : : const double M_vals[25] = { 1,
128 : : 1,
129 : : 1,
130 : : 1,
131 : : 1,
132 [ # # ]: 0 : pts[0] % pts[0],
133 [ # # ]: 0 : pts[0][0],
134 [ # # ]: 0 : pts[0][1],
135 [ # # ]: 0 : pts[0][2],
136 : : 1,
137 [ # # ]: 0 : pts[1] % pts[1],
138 [ # # ]: 0 : pts[1][0],
139 [ # # ]: 0 : pts[1][1],
140 [ # # ]: 0 : pts[1][2],
141 : : 1,
142 [ # # ]: 0 : pts[2] % pts[2],
143 [ # # ]: 0 : pts[2][0],
144 [ # # ]: 0 : pts[2][1],
145 [ # # ]: 0 : pts[2][2],
146 : : 1,
147 [ # # ]: 0 : pts[3] % pts[3],
148 [ # # ]: 0 : pts[3][0],
149 [ # # ]: 0 : pts[3][1],
150 [ # # ]: 0 : pts[3][2],
151 : 0 : 1 };
152 [ # # ]: 0 : MsqMatrix< 5, 5 > M( M_vals );
153 [ # # ][ # # ]: 0 : double M11 = det( MsqMatrix< 4, 4 >( M, 0, 0 ) );
154 [ # # ][ # # ]: 0 : double M12 = det( MsqMatrix< 4, 4 >( M, 0, 1 ) );
155 [ # # ][ # # ]: 0 : double M13 = det( MsqMatrix< 4, 4 >( M, 0, 2 ) );
156 [ # # ][ # # ]: 0 : double M14 = det( MsqMatrix< 4, 4 >( M, 0, 3 ) );
157 [ # # ][ # # ]: 0 : double M15 = det( MsqMatrix< 4, 4 >( M, 0, 4 ) );
158 : :
159 : : // define the sphere
160 [ # # ]: 0 : Vector3D cent( 0.5 * M12 / M11, -0.5 * M13 / M11, 0.5 * M14 / M11 );
161 [ # # ][ # # ]: 0 : this->set_sphere( cent, sqrt( cent % cent - M15 / M11 ) );
[ # # ]
162 [ + - ][ + - ]: 12 : }
|
