Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2006 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 : : (2006) [email protected]
24 : :
25 : : ***************************************************************** */
26 : :
27 : : /** \file TShapeSize2DB2.cpp
28 : : * \brief
29 : : * \author Jason Kraftcheck
30 : : */
31 : :
32 : : #include "Mesquite.hpp"
33 : : #include "TShapeSize2DB2.hpp"
34 : : #include "MsqMatrix.hpp"
35 : : #include "MsqError.hpp"
36 : : #include "TMPDerivs.hpp"
37 : :
38 : : #include <iostream>
39 : :
40 : : namespace MBMesquite
41 : : {
42 : :
43 : 0 : std::string TShapeSize2DB2::get_name() const
44 : : {
45 [ # # ]: 0 : return "TShapeSize2DB2";
46 : : }
47 : :
48 [ # # ]: 0 : TShapeSize2DB2::~TShapeSize2DB2() {}
49 : :
50 : 212 : bool TShapeSize2DB2::evaluate( const MsqMatrix< 2, 2 >& T, double& result, MsqError& err )
51 : : {
52 : 212 : const double two_det = 2.0 * det( T );
53 [ - + ]: 212 : if( invalid_determinant( two_det ) )
54 : : { // barrier
55 [ # # ]: 0 : MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
56 : 0 : return false;
57 : : }
58 : :
59 : 212 : const double frob_sqr = sqr_Frobenius( T );
60 : 212 : result = ( frob_sqr - 2.0 * sqrt( frob_sqr + two_det ) + 2.0 ) / two_det;
61 : 212 : return true;
62 : : }
63 : :
64 : 11112 : bool TShapeSize2DB2::evaluate_with_grad( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& deriv_wrt_T,
65 : : MsqError& err )
66 : : {
67 [ + - ]: 11112 : const double d = det( T );
68 [ + - ][ - + ]: 11112 : if( invalid_determinant( d ) )
69 : : { // barrier
70 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
71 : 0 : return false;
72 : : }
73 [ + - ]: 11112 : const double frob_sqr = sqr_Frobenius( T );
74 [ + - ]: 11112 : const double psi = sqrt( frob_sqr + 2.0 * det( T ) );
75 : 11112 : const double v = frob_sqr - 2.0 * psi + 2.0;
76 : 11112 : result = v / ( 2 * d );
77 : :
78 : : // deriv of V wrt T
79 [ + - ]: 11112 : MsqMatrix< 2, 2 > adjt = transpose_adj( T );
80 : 11112 : MsqMatrix< 2, 2 > v_wrt_T( T );
81 [ + - ]: 11112 : if( psi > 1e-50 )
82 : : {
83 [ + - ]: 11112 : v_wrt_T *= ( 1.0 - 1.0 / psi );
84 [ + - ][ + - ]: 11112 : v_wrt_T -= 1.0 / psi * adjt;
85 [ + - ]: 11112 : v_wrt_T *= 2;
86 : : }
87 : : else
88 : : {
89 [ # # ][ # # ]: 0 : std::cout << "Warning: Division by zero avoided in TShapeSize2DB2::evaluate_with_grad()" << std::endl;
90 : : }
91 : :
92 : : // deriv of mu wrt T
93 : 11112 : deriv_wrt_T = v_wrt_T;
94 [ + - ]: 11112 : deriv_wrt_T *= 0.5 / d;
95 [ + - ][ + - ]: 11112 : deriv_wrt_T -= v / ( 2 * d * d ) * adjt;
96 : :
97 : 11112 : return true;
98 : : }
99 : :
100 : 0 : bool TShapeSize2DB2::evaluate_with_hess( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& deriv_wrt_T,
101 : : MsqMatrix< 2, 2 > second[3], MsqError& err )
102 : : {
103 [ # # ]: 0 : const double d = det( T );
104 [ # # ][ # # ]: 0 : if( invalid_determinant( d ) )
105 : : { // barrier
106 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
107 : 0 : return false;
108 : : }
109 [ # # ]: 0 : const double frob_sqr = sqr_Frobenius( T );
110 [ # # ]: 0 : const double psi = sqrt( frob_sqr + 2.0 * det( T ) );
111 : 0 : const double v = frob_sqr - 2.0 * psi + 2.0;
112 : 0 : result = v / ( 2 * d );
113 : :
114 : : // deriv of V wrt T
115 [ # # ]: 0 : MsqMatrix< 2, 2 > adjt = transpose_adj( T );
116 : 0 : MsqMatrix< 2, 2 > v_wrt_T( T );
117 [ # # ]: 0 : if( psi > 1e-50 )
118 : : {
119 [ # # ]: 0 : v_wrt_T *= ( 1.0 - 1.0 / psi );
120 [ # # ][ # # ]: 0 : v_wrt_T -= 1.0 / psi * adjt;
121 [ # # ]: 0 : v_wrt_T *= 2;
122 : : }
123 : : else
124 : : {
125 [ # # ][ # # ]: 0 : std::cout << "Warning: Division by zero avoided in TShapeSize2DB2::evaluate_with_hess()" << std::endl;
126 : : }
127 : :
128 : : // deriv of mu wrt T
129 : 0 : deriv_wrt_T = v_wrt_T;
130 [ # # ]: 0 : deriv_wrt_T *= 0.5 / d;
131 [ # # ][ # # ]: 0 : deriv_wrt_T -= v / ( 2 * d * d ) * adjt;
132 : :
133 : : // second of V wrt T
134 [ # # ][ # # ]: 0 : const double s = T( 0, 1 ) - T( 1, 0 );
135 [ # # ]: 0 : const double tr = trace( T );
136 : 0 : const double f = -2.0 / ( psi * psi * psi );
137 [ # # ][ # # ]: 0 : second[0]( 0, 0 ) = second[1]( 0, 1 ) = second[2]( 1, 1 ) = f * s * s;
[ # # ]
138 [ # # ][ # # ]: 0 : second[0]( 0, 1 ) = second[0]( 1, 0 ) = second[1]( 1, 1 ) = -f * s * tr;
[ # # ]
139 [ # # ][ # # ]: 0 : second[1]( 0, 0 ) = second[2]( 0, 1 ) = second[2]( 1, 0 ) = f * s * tr;
[ # # ]
140 [ # # ][ # # ]: 0 : second[0]( 1, 1 ) = second[2]( 0, 0 ) = -( second[1]( 1, 0 ) = -f * tr * tr );
[ # # ]
141 [ # # ]: 0 : pluseq_scaled_I( second, 2 );
142 : :
143 : : // second of mu wrt T
144 : 0 : const double x = 1.0 / ( 2 * d );
145 [ # # ]: 0 : second[0] *= x;
146 [ # # ]: 0 : second[1] *= x;
147 [ # # ]: 0 : second[2] *= x;
148 [ # # ]: 0 : pluseq_scaled_2nd_deriv_of_det( second, v / ( -2 * d * d ) );
149 [ # # ]: 0 : pluseq_scaled_outer_product( second, v / ( d * d * d ), adjt );
150 [ # # ]: 0 : pluseq_scaled_sum_outer_product( second, -1 / ( 2 * d * d ), v_wrt_T, adjt );
151 : :
152 : 0 : return true;
153 : : }
154 : :
155 [ + - ][ + - ]: 8 : } // namespace MBMesquite
|
