MOAB: Mesh Oriented datABase  (version 5.2.1)
TetDihedralWeight.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2009 Sandia National Laboratories.  Developed at the
00005     University of Wisconsin--Madison under SNL contract number
00006     624796.  The U.S. Government and the University of Wisconsin
00007     retain certain rights to 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     (2009) kraftche@cae.wisc.edu
00024 
00025   ***************************************************************** */
00026 
00027 /** \file TetDihedralWeight.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 #ifdef WIN32               /* windows */
00032 #define _USE_MATH_DEFINES  // For M_PI
00033 #endif
00034 #include "Mesquite.hpp"
00035 #include "TetDihedralWeight.hpp"
00036 #include "PatchData.hpp"
00037 #include "ReferenceMesh.hpp"
00038 #include <algorithm>
00039 
00040 namespace MBMesquite
00041 {
00042 
00043 static inline double da( double dot )
00044 {
00045     return 180 - ( 180 / M_PI ) * acos( dot );
00046 }
00047 
00048 double TetDihedralWeight::get_weight( PatchData& pd, size_t element, Sample, MsqError& err )
00049 {
00050     const double eps = 1e-10;
00051 
00052     MsqMeshEntity& elem = pd.element_by_index( element );
00053     if( elem.get_element_type() != TETRAHEDRON )
00054     {
00055         MSQ_SETERR( err )( MsqError::UNSUPPORTED_ELEMENT );
00056         return 0.0;
00057     }
00058 
00059     const size_t* indices = elem.get_vertex_index_array();
00060     Vector3D v01, v02, v31, v32;
00061     if( refMesh )
00062     {
00063         const Mesh::VertexHandle* vtx_hdl = pd.get_vertex_handles_array();
00064         Mesh::VertexHandle handles[]      = { vtx_hdl[indices[0]], vtx_hdl[indices[1]], vtx_hdl[indices[2]],
00065                                          vtx_hdl[indices[3]] };
00066         Vector3D coords[4];
00067         refMesh->get_reference_vertex_coordinates( handles, 4, coords, err );
00068         MSQ_ERRZERO( err );
00069 
00070         v01 = coords[1] - coords[0];
00071         v02 = coords[2] - coords[0];
00072         v31 = coords[1] - coords[3];
00073         v32 = coords[2] - coords[3];
00074     }
00075     else
00076     {
00077         const MsqVertex* coords = pd.get_vertex_array();
00078         v01                     = coords[indices[1]] - coords[indices[0]];
00079         v02                     = coords[indices[2]] - coords[indices[0]];
00080         v31                     = coords[indices[1]] - coords[indices[3]];
00081         v32                     = coords[indices[2]] - coords[indices[3]];
00082     }
00083 
00084     Vector3D n012 = v02 * v01;
00085     Vector3D n013 = v31 * v01;
00086     Vector3D n023 = v02 * v32;
00087     Vector3D n123 = v31 * v32;
00088 
00089     // normalize face vectors.
00090     double l012 = n012.length();
00091     double l013 = n013.length();
00092     double l023 = n023.length();
00093     double l123 = n123.length();
00094     n012 *= ( l012 < eps ) ? 0.0 : 1.0 / l012;
00095     n013 *= ( l013 < eps ) ? 0.0 : 1.0 / l013;
00096     n023 *= ( l023 < eps ) ? 0.0 : 1.0 / l023;
00097     n123 *= ( l123 < eps ) ? 0.0 : 1.0 / l123;
00098 
00099     // calculate dihedral handles for each edge
00100     double ds[] = { da( n012 % n013 ), da( n012 % n123 ), da( n012 % n023 ),
00101                     da( n013 % n023 ), da( n013 % n123 ), da( n023 % n123 ) };
00102 
00103     // calculate weight from max dihedral handle
00104     double d = *std::max_element( ds, ds + 6 );
00105     return 1 / ( 1 + exp( -mA * ( d - mCutoff ) ) );
00106 }
00107 
00108 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines