MOAB: Mesh Oriented datABase  (version 5.4.0)
SphericalDomain.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2004 Sandia Corporation and Argonne National
00005     Laboratory.  Under the terms of Contract DE-AC04-94AL85000
00006     with Sandia Corporation, the U.S. Government retains certain
00007     rights in this software.
00008 
00009     This library is free software; you can redistribute it and/or
00010     modify it under the terms of the GNU Lesser General Public
00011     License as published by the Free Software Foundation; either
00012     version 2.1 of the License, or (at your option) any later version.
00013 
00014     This library is distributed in the hope that it will be useful,
00015     but WITHOUT ANY WARRANTY; without even the implied warranty of
00016     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017     Lesser General Public License for more details.
00018 
00019     You should have received a copy of the GNU Lesser General Public License
00020     (lgpl.txt) along with this library; if not, write to the Free Software
00021     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 
00023     diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov,
00024     pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov,
00025     kraftche@cae.wisc.edu
00026 
00027   ***************************************************************** */
00028 #include "Mesquite.hpp"
00029 #include "SphericalDomain.hpp"
00030 #include "Vector3D.hpp"
00031 #include "MsqError.hpp"
00032 #include "MsqVertex.hpp"
00033 #include "DomainUtil.hpp"
00034 #include "MsqMatrix.hpp"
00035 #include "moab/Util.hpp"
00036 
00037 #ifdef MSQ_HAVE_IEEEFP_H
00038 #include <ieeefp.h>
00039 #endif
00040 
00041 #include <algorithm>
00042 
00043 MBMesquite::SphericalDomain::~SphericalDomain() {}
00044 
00045 void MBMesquite::SphericalDomain::snap_to( Mesh::VertexHandle /*entity_handle*/, Vector3D& coordinate ) const
00046 {
00047     // Get vector center to coordinate, store in coordinate.
00048     coordinate -= mCenter;
00049     // Get distance from center of sphere
00050     double len = coordinate.length();
00051     // Scale vector to have length of radius
00052     coordinate *= mRadius / len;
00053     // If was at center, return arbitrary position on sphere
00054     // (all possitions are equally close)
00055     if( !moab::Util::is_finite( coordinate.x() ) ) coordinate.set( mRadius, 0.0, 0.0 );
00056     // Get position from vector
00057     coordinate += mCenter;
00058 }
00059 
00060 void MBMesquite::SphericalDomain::vertex_normal_at( Mesh::VertexHandle /*entity_handle*/, Vector3D& coordinate ) const
00061 {
00062     // normal is vector from center to input position
00063     coordinate -= mCenter;
00064     // make it a unit vector
00065     double length = coordinate.length();
00066     coordinate /= length;
00067     // if input position was at center, choose same position
00068     // on sphere as snap_to.
00069     if( !moab::Util::is_finite( coordinate.x() ) ) coordinate.set( 1.0, 0.0, 0.0 );
00070 }
00071 void MBMesquite::SphericalDomain::element_normal_at( Mesh::ElementHandle h, Vector3D& coordinate ) const
00072 {
00073     SphericalDomain::vertex_normal_at( h, coordinate );
00074 }
00075 
00076 void MBMesquite::SphericalDomain::vertex_normal_at( const MBMesquite::Mesh::VertexHandle* handle,
00077                                                     MBMesquite::Vector3D coords[],
00078                                                     unsigned count,
00079                                                     MBMesquite::MsqError& ) const
00080 {
00081     for( unsigned i = 0; i < count; ++i )
00082         vertex_normal_at( handle[i], coords[i] );
00083 }
00084 
00085 void MBMesquite::SphericalDomain::closest_point( MBMesquite::Mesh::VertexHandle,
00086                                                  const MBMesquite::Vector3D& position,
00087                                                  MBMesquite::Vector3D& closest,
00088                                                  MBMesquite::Vector3D& normal,
00089                                                  MBMesquite::MsqError& ) const
00090 {
00091     normal = position - mCenter;
00092     normal.normalize();
00093     if( !moab::Util::is_finite( normal.x() ) ) normal.set( 1.0, 0.0, 0.0 );
00094     closest = mCenter + mRadius * normal;
00095 }
00096 
00097 void MBMesquite::SphericalDomain::domain_DoF( const Mesh::VertexHandle*,
00098                                               unsigned short* dof_array,
00099                                               size_t num_vertices,
00100                                               MsqError& ) const
00101 {
00102     std::fill( dof_array, dof_array + num_vertices, 2 );
00103 }
00104 
00105 void MBMesquite::SphericalDomain::fit_vertices( Mesh* mesh, MsqError& err, double epsilon )
00106 {
00107     std::vector< Mesh::VertexHandle > verts;
00108     mesh->get_all_vertices( verts, err );
00109     if( !MSQ_CHKERR( err ) ) fit_vertices( mesh, arrptr( verts ), verts.size(), err, epsilon );
00110 }
00111 
00112 void MBMesquite::SphericalDomain::fit_vertices( Mesh* mesh,
00113                                                 const Mesh::VertexHandle* verts,
00114                                                 size_t num_verts,
00115                                                 MsqError& err,
00116                                                 double epsilon )
00117 {
00118     std::vector< MsqVertex > coords( num_verts );
00119     mesh->vertices_get_coordinates( verts, arrptr( coords ), num_verts, err );MSQ_ERRRTN( err );
00120 
00121     if( epsilon <= 0.0 ) epsilon = DomainUtil::default_tolerance( arrptr( coords ), num_verts );
00122 
00123     Vector3D pts[4];
00124     if( !DomainUtil::non_coplanar_vertices( arrptr( coords ), num_verts, pts, epsilon ) )
00125     {
00126         MSQ_SETERR( err )( "All vertices are co-planar", MsqError::INVALID_MESH );
00127         return;
00128     }
00129 
00130     // solve deterinant form of four-point sphere
00131 
00132     // Define the bottom 4 rows of a 5x5 matrix.  The top
00133     // row contains the variables we are solving for, so just
00134     // fill it with ones.
00135     const double M_vals[25] = { 1,
00136                                 1,
00137                                 1,
00138                                 1,
00139                                 1,
00140                                 pts[0] % pts[0],
00141                                 pts[0][0],
00142                                 pts[0][1],
00143                                 pts[0][2],
00144                                 1,
00145                                 pts[1] % pts[1],
00146                                 pts[1][0],
00147                                 pts[1][1],
00148                                 pts[1][2],
00149                                 1,
00150                                 pts[2] % pts[2],
00151                                 pts[2][0],
00152                                 pts[2][1],
00153                                 pts[2][2],
00154                                 1,
00155                                 pts[3] % pts[3],
00156                                 pts[3][0],
00157                                 pts[3][1],
00158                                 pts[3][2],
00159                                 1 };
00160     MsqMatrix< 5, 5 > M( M_vals );
00161     double M11 = det( MsqMatrix< 4, 4 >( M, 0, 0 ) );
00162     double M12 = det( MsqMatrix< 4, 4 >( M, 0, 1 ) );
00163     double M13 = det( MsqMatrix< 4, 4 >( M, 0, 2 ) );
00164     double M14 = det( MsqMatrix< 4, 4 >( M, 0, 3 ) );
00165     double M15 = det( MsqMatrix< 4, 4 >( M, 0, 4 ) );
00166 
00167     // define the sphere
00168     Vector3D cent( 0.5 * M12 / M11, -0.5 * M13 / M11, 0.5 * M14 / M11 );
00169     this->set_sphere( cent, sqrt( cent % cent - M15 / M11 ) );
00170 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines