MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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 00026 ***************************************************************** */ 00027 // -*- Mode : c++; tab-width: 3; c-tab-always-indent: t; indent-tabs-mode: nil; c-basic-offset: 3 00028 // -*- 00029 00030 /*! \file ConditionNumberFunctions.hpp 00031 00032 Header file for the MBMesquite::ShapeQualityMetric class 00033 00034 \author Thomas Leurent 00035 \date 2002-09-01 00036 */ 00037 00038 #ifndef ConditionNumberFunctions_hpp 00039 #define ConditionNumberFunctions_hpp 00040 00041 #include "Mesquite.hpp" 00042 #include "MsqError.hpp" 00043 #include "QualityMetric.hpp" 00044 #include "PatchData.hpp" 00045 00046 namespace MBMesquite 00047 { 00048 static inline bool condition_number_2d( const Vector3D temp_vec[], size_t e_ind, PatchData& pd, double& fval, 00049 MsqError& err ) 00050 { 00051 // norm squared of J 00052 double term1 = temp_vec[0] % temp_vec[0] + temp_vec[1] % temp_vec[1]; 00053 00054 Vector3D unit_surf_norm; 00055 pd.get_domain_normal_at_element( e_ind, unit_surf_norm, err ); 00056 MSQ_ERRZERO( err ); 00057 unit_surf_norm.normalize(); 00058 00059 // det J 00060 double temp_var = unit_surf_norm % ( temp_vec[0] * temp_vec[1] ); 00061 00062 // revert to old, non-barrier form 00063 if( temp_var <= 0.0 ) return false; 00064 fval = term1 / ( 2.0 * temp_var ); 00065 return true; 00066 00067 /* 00068 double h; 00069 double delta=pd.get_barrier_delta(err); MSQ_ERRZERO(err); 00070 00071 // Note: technically, we want delta=eta*tau-max 00072 // whereas the function above gives delta=eta*alpha-max 00073 // 00074 // Because the only requirement on eta is eta << 1, 00075 // and because tau-max = alpha-max/0.707 we can 00076 // ignore the discrepancy 00077 00078 if (delta==0) { 00079 if (temp_var < MSQ_DBL_MIN ) { 00080 return false; 00081 } 00082 else { 00083 h=temp_var; 00084 } 00085 00086 // Note: when delta=0, the vertex_barrier_function 00087 // formally gives h=temp_var as well. 00088 // We just do it this way to avoid any 00089 // roundoff issues. 00090 // Also: when delta=0, this metric is identical 00091 // to the original condition number with 00092 // the barrier at temp_var=0 00093 } 00094 else { 00095 h = QualityMetric::vertex_barrier_function(temp_var,delta); 00096 00097 if (h<MSQ_DBL_MIN && fabs(temp_var) > MSQ_DBL_MIN ) { 00098 h = delta*delta/fabs(temp_var); } 00099 // Note: Analytically, h is strictly positive, but 00100 // it can be zero numerically if temp_var 00101 // is a large negative number 00102 // In the case h=0, we use a different analytic 00103 // approximation to compute h. 00104 } 00105 00106 if (h<MSQ_DBL_MIN) { 00107 MSQ_SETERR(err)( "Barrier function is zero due to excessively large " 00108 "negative area compared to delta. /n Try to untangle " 00109 "mesh another way. ", MsqError::INVALID_MESH); 00110 return false; 00111 } 00112 00113 fval=term1/(2*h); 00114 00115 if (fval>MSQ_MAX_CAP) { 00116 fval=MSQ_MAX_CAP; 00117 } 00118 return true; 00119 */ 00120 } 00121 00122 //} //namespace 00123 00124 static inline bool condition_number_3d( const Vector3D temp_vec[], PatchData& /*pd*/, double& fval, MsqError& /*err*/ ) 00125 { 00126 // norm squared of J 00127 double term1 = temp_vec[0] % temp_vec[0] + temp_vec[1] % temp_vec[1] + temp_vec[2] % temp_vec[2]; 00128 // norm squared of adjoint of J 00129 double term2 = ( temp_vec[0] * temp_vec[1] ) % ( temp_vec[0] * temp_vec[1] ) + 00130 ( temp_vec[1] * temp_vec[2] ) % ( temp_vec[1] * temp_vec[2] ) + 00131 ( temp_vec[2] * temp_vec[0] ) % ( temp_vec[2] * temp_vec[0] ); 00132 // det of J 00133 double temp_var = temp_vec[0] % ( temp_vec[1] * temp_vec[2] ); 00134 00135 // revert to old, non-barrier formulation 00136 if( temp_var <= 0.0 ) return false; 00137 fval = sqrt( term1 * term2 ) / ( 3.0 * temp_var ); 00138 return true; 00139 00140 /* 00141 double h; 00142 double delta=pd.get_barrier_delta(err); MSQ_ERRZERO(err); 00143 00144 // Note: technically, we want delta=eta*tau-max 00145 // whereas the function above gives delta=eta*alpha-max 00146 // 00147 // Because the only requirement on eta is eta << 1, 00148 // and because tau-max = alpha-max/0.707 we can 00149 // ignore the discrepancy 00150 00151 if (delta==0) { 00152 if (temp_var < MSQ_DBL_MIN ) { 00153 return false; 00154 } 00155 else { 00156 h=temp_var; 00157 } 00158 00159 // Note: when delta=0, the vertex_barrier_function 00160 // formally gives h=temp_var as well. 00161 // We just do it this way to avoid any 00162 // roundoff issues. 00163 // Also: when delta=0, this metric is identical 00164 // to the original condition number with 00165 // the barrier at temp_var=0 00166 00167 } 00168 else { 00169 h = QualityMetric::vertex_barrier_function(temp_var,delta); 00170 00171 if (h<MSQ_DBL_MIN && fabs(temp_var) > MSQ_DBL_MIN ) { 00172 h = delta*delta/fabs(temp_var); } 00173 00174 // Note: Analytically, h is strictly positive, but 00175 // it can be zero numerically if temp_var 00176 // is a large negative number 00177 // In the h=0, we use a different analytic 00178 // approximation to compute h. 00179 } 00180 if (h<MSQ_DBL_MIN) { 00181 MSQ_SETERR(err)("Barrier function is zero due to excessively large " 00182 "negative area compared to delta. /n Try to untangle " 00183 "mesh another way. ", MsqError::INVALID_MESH); 00184 return false; 00185 } 00186 00187 fval=sqrt(term1*term2)/(3*h); 00188 00189 if (fval>MSQ_MAX_CAP) { 00190 fval=MSQ_MAX_CAP; 00191 } 00192 return true; 00193 */ 00194 } 00195 00196 } // namespace MBMesquite 00197 00198 #endif // ConditionNumberFunctions_hpp