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