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 : :
26 : : ***************************************************************** */
27 : : // -*- Mode : c++; tab-width: 3; c-tab-always-indent: t; indent-tabs-mode: nil; c-basic-offset: 3
28 : : // -*-
29 : :
30 : : /*! \file ConditionNumberFunctions.hpp
31 : :
32 : : Header file for the MBMesquite::ShapeQualityMetric class
33 : :
34 : : \author Thomas Leurent
35 : : \date 2002-09-01
36 : : */
37 : :
38 : : #ifndef ConditionNumberFunctions_hpp
39 : : #define ConditionNumberFunctions_hpp
40 : :
41 : : #include "Mesquite.hpp"
42 : : #include "MsqError.hpp"
43 : : #include "QualityMetric.hpp"
44 : : #include "PatchData.hpp"
45 : :
46 : : namespace MBMesquite
47 : : {
48 : 2161 : static inline bool condition_number_2d( const Vector3D temp_vec[], size_t e_ind, PatchData& pd, double& fval,
49 : : MsqError& err )
50 : : {
51 : : // norm squared of J
52 [ + - ][ + - ]: 2161 : double term1 = temp_vec[0] % temp_vec[0] + temp_vec[1] % temp_vec[1];
53 : :
54 [ + - ]: 2161 : Vector3D unit_surf_norm;
55 [ + - ]: 2161 : pd.get_domain_normal_at_element( e_ind, unit_surf_norm, err );
56 [ + - ][ - + ]: 2161 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
57 [ + - ]: 2161 : unit_surf_norm.normalize();
58 : :
59 : : // det J
60 [ + - ][ + - ]: 2161 : double temp_var = unit_surf_norm % ( temp_vec[0] * temp_vec[1] );
61 : :
62 : : // revert to old, non-barrier form
63 [ + + ]: 2161 : if( temp_var <= 0.0 ) return false;
64 : 2154 : fval = term1 / ( 2.0 * temp_var );
65 : 2161 : return true;
66 : :
67 : : /*
68 : : double h;
69 : : double delta=pd.get_barrier_delta(err); MSQ_ERRZERO(err);
70 : :
71 : : // Note: technically, we want delta=eta*tau-max
72 : : // whereas the function above gives delta=eta*alpha-max
73 : : //
74 : : // Because the only requirement on eta is eta << 1,
75 : : // and because tau-max = alpha-max/0.707 we can
76 : : // ignore the discrepancy
77 : :
78 : : if (delta==0) {
79 : : if (temp_var < MSQ_DBL_MIN ) {
80 : : return false;
81 : : }
82 : : else {
83 : : h=temp_var;
84 : : }
85 : :
86 : : // Note: when delta=0, the vertex_barrier_function
87 : : // formally gives h=temp_var as well.
88 : : // We just do it this way to avoid any
89 : : // roundoff issues.
90 : : // Also: when delta=0, this metric is identical
91 : : // to the original condition number with
92 : : // the barrier at temp_var=0
93 : : }
94 : : else {
95 : : h = QualityMetric::vertex_barrier_function(temp_var,delta);
96 : :
97 : : if (h<MSQ_DBL_MIN && fabs(temp_var) > MSQ_DBL_MIN ) {
98 : : h = delta*delta/fabs(temp_var); }
99 : : // Note: Analytically, h is strictly positive, but
100 : : // it can be zero numerically if temp_var
101 : : // is a large negative number
102 : : // In the case h=0, we use a different analytic
103 : : // approximation to compute h.
104 : : }
105 : :
106 : : if (h<MSQ_DBL_MIN) {
107 : : MSQ_SETERR(err)( "Barrier function is zero due to excessively large "
108 : : "negative area compared to delta. /n Try to untangle "
109 : : "mesh another way. ", MsqError::INVALID_MESH);
110 : : return false;
111 : : }
112 : :
113 : : fval=term1/(2*h);
114 : :
115 : : if (fval>MSQ_MAX_CAP) {
116 : : fval=MSQ_MAX_CAP;
117 : : }
118 : : return true;
119 : : */
120 : : }
121 : :
122 : : //} //namespace
123 : :
124 : 2029860 : static inline bool condition_number_3d( const Vector3D temp_vec[], PatchData& /*pd*/, double& fval, MsqError& /*err*/ )
125 : : {
126 : : // norm squared of J
127 : 2029860 : double term1 = temp_vec[0] % temp_vec[0] + temp_vec[1] % temp_vec[1] + temp_vec[2] % temp_vec[2];
128 : : // norm squared of adjoint of J
129 [ + - ][ + - ]: 2029860 : double term2 = ( temp_vec[0] * temp_vec[1] ) % ( temp_vec[0] * temp_vec[1] ) +
130 [ + - ][ + - ]: 2029860 : ( temp_vec[1] * temp_vec[2] ) % ( temp_vec[1] * temp_vec[2] ) +
[ + - ]
131 [ + - ][ + - ]: 2029860 : ( temp_vec[2] * temp_vec[0] ) % ( temp_vec[2] * temp_vec[0] );
[ + - ]
132 : : // det of J
133 [ + - ]: 2029860 : double temp_var = temp_vec[0] % ( temp_vec[1] * temp_vec[2] );
134 : :
135 : : // revert to old, non-barrier formulation
136 [ - + ]: 2029860 : if( temp_var <= 0.0 ) return false;
137 : 2029860 : fval = sqrt( term1 * term2 ) / ( 3.0 * temp_var );
138 : 2029860 : return true;
139 : :
140 : : /*
141 : : double h;
142 : : double delta=pd.get_barrier_delta(err); MSQ_ERRZERO(err);
143 : :
144 : : // Note: technically, we want delta=eta*tau-max
145 : : // whereas the function above gives delta=eta*alpha-max
146 : : //
147 : : // Because the only requirement on eta is eta << 1,
148 : : // and because tau-max = alpha-max/0.707 we can
149 : : // ignore the discrepancy
150 : :
151 : : if (delta==0) {
152 : : if (temp_var < MSQ_DBL_MIN ) {
153 : : return false;
154 : : }
155 : : else {
156 : : h=temp_var;
157 : : }
158 : :
159 : : // Note: when delta=0, the vertex_barrier_function
160 : : // formally gives h=temp_var as well.
161 : : // We just do it this way to avoid any
162 : : // roundoff issues.
163 : : // Also: when delta=0, this metric is identical
164 : : // to the original condition number with
165 : : // the barrier at temp_var=0
166 : :
167 : : }
168 : : else {
169 : : h = QualityMetric::vertex_barrier_function(temp_var,delta);
170 : :
171 : : if (h<MSQ_DBL_MIN && fabs(temp_var) > MSQ_DBL_MIN ) {
172 : : h = delta*delta/fabs(temp_var); }
173 : :
174 : : // Note: Analytically, h is strictly positive, but
175 : : // it can be zero numerically if temp_var
176 : : // is a large negative number
177 : : // In the h=0, we use a different analytic
178 : : // approximation to compute h.
179 : : }
180 : : if (h<MSQ_DBL_MIN) {
181 : : MSQ_SETERR(err)("Barrier function is zero due to excessively large "
182 : : "negative area compared to delta. /n Try to untangle "
183 : : "mesh another way. ", MsqError::INVALID_MESH);
184 : : return false;
185 : : }
186 : :
187 : : fval=sqrt(term1*term2)/(3*h);
188 : :
189 : : if (fval>MSQ_MAX_CAP) {
190 : : fval=MSQ_MAX_CAP;
191 : : }
192 : : return true;
193 : : */
194 : : }
195 : :
196 : : } // namespace MBMesquite
197 : :
198 : : #endif // ConditionNumberFunctions_hpp
|
