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