Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2010 Sandia National Laboratories. Developed at the
5 : : University of Wisconsin--Madison under SNL contract number
6 : : 624796. The U.S. Government and the University of Wisconsin
7 : : retain certain rights to 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 : : (2010) [email protected]
24 : :
25 : : ***************************************************************** */
26 : :
27 : : /** \file TMetric.cpp
28 : : * \brief
29 : : * \author Jason Kraftcheck
30 : : */
31 : :
32 : : #include "Mesquite.hpp"
33 : : #include "TMetric.hpp"
34 : : #include "TMetricBarrier.hpp"
35 : : #include "MsqMatrix.hpp"
36 : : #include "MsqError.hpp"
37 : : #include <limits>
38 : :
39 : : namespace MBMesquite
40 : : {
41 : :
42 : : template < unsigned Dim >
43 : 437204 : static inline double do_finite_difference( int r, int c, TMetric* metric, MsqMatrix< Dim, Dim > A, double value,
44 : : MsqError& err )
45 : : {
46 [ # # ][ + - ]: 437204 : const double INITIAL_STEP = std::max( 1e-6, fabs( 1e-14 * value ) );
47 [ # # ][ + - ]: 437204 : const double init = A( r, c );
48 : : bool valid;
49 : : double diff_value;
50 [ # # ][ + - ]: 437204 : for( double step = INITIAL_STEP; step > std::numeric_limits< double >::epsilon(); step *= 0.1 )
51 : : {
52 [ # # ][ + - ]: 437204 : A( r, c ) = init + step;
53 [ # # ][ + - ]: 437204 : valid = metric->evaluate( A, diff_value, err );
54 [ # # ][ # # ]: 437204 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ][ + - ]
[ - + ][ # # ]
[ # # ][ - + ]
55 [ # # ][ + - ]: 437204 : if( valid ) return ( diff_value - value ) / step;
56 : : }
57 : :
58 : : // If we couldn't find a valid step, try stepping in the other
59 : : // direciton
60 [ # # ][ # # ]: 0 : for( double step = INITIAL_STEP; step > std::numeric_limits< double >::epsilon(); step *= 0.1 )
61 : : {
62 [ # # ][ # # ]: 0 : A( r, c ) = init - step;
63 [ # # ][ # # ]: 0 : valid = metric->evaluate( A, diff_value, err );
64 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
65 [ # # ][ # # ]: 0 : if( valid ) return ( value - diff_value ) / step;
66 : : }
67 : : // If that didn't work either, then give up.
68 [ # # ][ # # ]: 0 : MSQ_SETERR( err )
69 [ # # ][ # # ]: 0 : ( "No valid step size for finite difference of 2D target metric.", MsqError::INTERNAL_ERROR );
70 : 437204 : return 0.0;
71 : : }
72 : :
73 : : template < unsigned Dim >
74 : 109318 : static inline bool do_numerical_gradient( TMetric* mu, MsqMatrix< Dim, Dim > A, double& result,
75 : : MsqMatrix< Dim, Dim >& wrt_A, MsqError& err )
76 : : {
77 : : bool valid;
78 : 109318 : valid = mu->evaluate( A, result, err );
79 [ # # ][ # # ]: 109318 : MSQ_ERRZERO( err );
[ # # ][ + + ]
[ + - ][ + + ]
80 [ # # ][ # # ]: 109301 : if( MSQ_CHKERR( err ) || !valid ) return valid;
[ # # ][ # # ]
[ - + ][ # # ]
[ - + ][ - + ]
81 : :
82 : : switch( Dim )
83 : : {
84 : : case 3:
85 : 0 : wrt_A( 0, 2 ) = do_finite_difference( 0, 2, mu, A, result, err );
86 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ]
87 : 0 : wrt_A( 1, 2 ) = do_finite_difference( 1, 2, mu, A, result, err );
88 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ]
89 : 0 : wrt_A( 2, 0 ) = do_finite_difference( 2, 0, mu, A, result, err );
90 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ]
91 : 0 : wrt_A( 2, 1 ) = do_finite_difference( 2, 1, mu, A, result, err );
92 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ]
93 : 0 : wrt_A( 2, 2 ) = do_finite_difference( 2, 2, mu, A, result, err );
94 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ]
95 : : case 2:
96 : 109301 : wrt_A( 0, 1 ) = do_finite_difference( 0, 1, mu, A, result, err );
97 [ # # ][ # # ]: 109301 : MSQ_ERRZERO( err );
[ # # ][ - + ]
[ # # ][ - + ]
98 : 109301 : wrt_A( 1, 0 ) = do_finite_difference( 1, 0, mu, A, result, err );
99 [ # # ][ # # ]: 109301 : MSQ_ERRZERO( err );
[ # # ][ - + ]
[ # # ][ - + ]
100 : 109301 : wrt_A( 1, 1 ) = do_finite_difference( 1, 1, mu, A, result, err );
101 [ # # ][ # # ]: 109301 : MSQ_ERRZERO( err );
[ # # ][ - + ]
[ # # ][ - + ]
102 : : case 1:
103 : 109301 : wrt_A( 0, 0 ) = do_finite_difference( 0, 0, mu, A, result, err );
104 [ # # ][ # # ]: 109301 : MSQ_ERRZERO( err );
[ # # ][ - + ]
[ # # ][ - + ]
105 : 109301 : break;
106 : : default:
107 : : assert( false );
108 : : }
109 : 109301 : return true;
110 : : }
111 : :
112 : : template < unsigned Dim >
113 : 0 : static inline bool do_numerical_hessian( TMetric* metric, MsqMatrix< Dim, Dim > A, double& value,
114 : : MsqMatrix< Dim, Dim >& grad, MsqMatrix< Dim, Dim >* Hess, MsqError& err )
115 : : {
116 : : // zero hessian data
117 : 0 : const int num_block = Dim * ( Dim + 1 ) / 2;
118 [ # # ][ # # ]: 0 : for( int i = 0; i < num_block; ++i )
119 [ # # ][ # # ]: 0 : Hess[i].zero();
120 : :
121 : : // evaluate gradient for input values
122 : : bool valid;
123 [ # # ][ # # ]: 0 : valid = metric->evaluate_with_grad( A, value, grad, err );
124 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) || !valid ) return false;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
125 : :
126 : : // do finite difference for each term of A
127 [ # # ][ # # ]: 0 : const double INITAL_STEP = std::max( 1e-6, fabs( 1e-14 * value ) );
128 : : double value2;
129 [ # # ][ # # ]: 0 : MsqMatrix< Dim, Dim > grad2;
130 [ # # ][ # # ]: 0 : for( unsigned r = 0; r < Dim; ++r )
131 : : { // for each row of A
132 [ # # ][ # # ]: 0 : for( unsigned c = 0; c < Dim; ++c )
133 : : { // for each column of A
134 [ # # ][ # # ]: 0 : const double in_val = A( r, c );
135 : : double step;
136 [ # # ][ # # ]: 0 : for( step = INITAL_STEP; step > std::numeric_limits< double >::epsilon(); step *= 0.1 )
137 : : {
138 [ # # ][ # # ]: 0 : A( r, c ) = in_val + step;
139 [ # # ][ # # ]: 0 : valid = metric->evaluate_with_grad( A, value2, grad2, err );
140 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
141 [ # # ][ # # ]: 0 : if( valid ) break;
142 : : }
143 : :
144 : : // if no valid step size, try step in other direction
145 [ # # ][ # # ]: 0 : if( !valid )
146 : : {
147 [ # # ][ # # ]: 0 : for( step = -INITAL_STEP; step < -std::numeric_limits< double >::epsilon(); step *= 0.1 )
148 : : {
149 [ # # ][ # # ]: 0 : A( r, c ) = in_val + step;
150 [ # # ][ # # ]: 0 : valid = metric->evaluate_with_grad( A, value2, grad2, err );
151 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
152 [ # # ][ # # ]: 0 : if( valid ) break;
153 : : }
154 : :
155 : : // if still no valid step size, give up.
156 [ # # ][ # # ]: 0 : if( !valid )
157 : : {
158 [ # # ][ # # ]: 0 : MSQ_SETERR( err )
159 [ # # ][ # # ]: 0 : ( "No valid step size for finite difference of 2D target metric.", MsqError::INTERNAL_ERROR );
160 : 0 : return false;
161 : : }
162 : : }
163 : :
164 : : // restore A.
165 [ # # ][ # # ]: 0 : A( r, c ) = in_val;
166 : :
167 : : // add values into result matrix
168 : : // values of grad2, in row-major order, are a single 9-value row of the Hessian
169 [ # # ][ # # ]: 0 : grad2 -= grad;
170 [ # # ][ # # ]: 0 : grad2 /= step;
171 [ # # ][ # # ]: 0 : for( unsigned b = 0; b < r; ++b )
172 : : {
173 : 0 : const int idx = Dim * b - b * ( b + 1 ) / 2 + r;
174 [ # # ][ # # ]: 0 : Hess[idx].add_column( c, transpose( grad2.row( b ) ) );
[ # # ][ # # ]
[ # # ][ # # ]
175 : : }
176 [ # # ][ # # ]: 0 : for( unsigned b = r; b < Dim; ++b )
177 : : {
178 : 0 : const int idx = Dim * r - r * ( r + 1 ) / 2 + b;
179 [ # # ][ # # ]: 0 : Hess[idx].add_row( c, grad2.row( b ) );
[ # # ][ # # ]
180 : : }
181 : : } // for (c)
182 : : } // for (r)
183 : :
184 : : // Values in non-diagonal blocks were added twice.
185 [ # # ][ # # ]: 0 : for( unsigned r = 0, h = 1; r < Dim - 1; ++r, ++h )
186 [ # # ][ # # ]: 0 : for( unsigned c = r + 1; c < Dim; ++c, ++h )
187 [ # # ][ # # ]: 0 : Hess[h] *= 0.5;
188 : :
189 : 0 : return true;
190 : : }
191 : :
192 [ - + ]: 300 : TMetric::~TMetric() {}
193 : :
194 : 0 : bool TMetric::evaluate( const MsqMatrix< 2, 2 >& /*T*/, double& /*result*/, MsqError& /*err*/ )
195 : : {
196 : 0 : return false;
197 : : }
198 : :
199 : 0 : bool TMetric::evaluate( const MsqMatrix< 3, 3 >& /*T*/, double& /*result*/, MsqError& /*err*/ )
200 : : {
201 : 0 : return false;
202 : : }
203 : :
204 : 109318 : bool TMetric::evaluate_with_grad( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& wrt_T, MsqError& err )
205 : : {
206 : 109318 : return do_numerical_gradient( this, T, result, wrt_T, err );
207 : : }
208 : :
209 : 0 : bool TMetric::evaluate_with_grad( const MsqMatrix< 3, 3 >& T, double& result, MsqMatrix< 3, 3 >& wrt_T, MsqError& err )
210 : : {
211 : 0 : return do_numerical_gradient( this, T, result, wrt_T, err );
212 : : }
213 : :
214 : 0 : bool TMetric::evaluate_with_hess( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& deriv_wrt_T,
215 : : MsqMatrix< 2, 2 > hess_wrt_T[3], MsqError& err )
216 : : {
217 : 0 : return do_numerical_hessian( this, T, result, deriv_wrt_T, hess_wrt_T, err );
218 : : }
219 : :
220 : 0 : bool TMetric::evaluate_with_hess( const MsqMatrix< 3, 3 >& T, double& result, MsqMatrix< 3, 3 >& deriv_wrt_T,
221 : : MsqMatrix< 3, 3 > hess_wrt_T[6], MsqError& err )
222 : : {
223 : 0 : return do_numerical_hessian( this, T, result, deriv_wrt_T, hess_wrt_T, err );
224 : : }
225 : :
226 [ # # ]: 0 : TMetric2D::~TMetric2D() {}
227 [ # # ]: 0 : TMetric3D::~TMetric3D() {}
228 : :
229 : 0 : bool TMetric2D::evaluate( const MsqMatrix< 3, 3 >&, double&, MsqError& err )
230 : : {
231 : : MSQ_SETERR( err )
232 [ # # ]: 0 : ( "2D target metric cannot be evaluated for volume elements", MsqError::UNSUPPORTED_ELEMENT );
233 : 0 : return false;
234 : : }
235 : :
236 : 0 : bool TMetric3D::evaluate( const MsqMatrix< 2, 2 >&, double&, MsqError& err )
237 : : {
238 : : MSQ_SETERR( err )
239 [ # # ]: 0 : ( "2D target metric cannot be evaluated for volume elements", MsqError::UNSUPPORTED_ELEMENT );
240 : 0 : return false;
241 : : }
242 : :
243 [ + - ][ + - ]: 120 : } // namespace MBMesquite
|
