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 TShapeSizeB1.cpp
28 : : * \brief
29 : : * \author Jason Kraftcheck
30 : : */
31 : :
32 : : #include "Mesquite.hpp"
33 : : #include "TShapeSizeB1.hpp"
34 : : #include "TMPDerivs.hpp"
35 : : #include "MsqError.hpp"
36 : :
37 : : namespace MBMesquite
38 : : {
39 : :
40 : 6 : std::string TShapeSizeB1::get_name() const
41 : : {
42 [ + - ]: 6 : return "TShapeSizeB1";
43 : : }
44 : :
45 [ - + ]: 12 : TShapeSizeB1::~TShapeSizeB1() {}
46 : :
47 : 955844 : bool TShapeSizeB1::evaluate( const MsqMatrix< 2, 2 >& T, double& result, MsqError& /*err*/ )
48 : : {
49 : 955844 : const double tau = det( T );
50 [ + + ]: 955844 : if( invalid_determinant( tau ) )
51 : : { // barrier
52 : 12 : return false;
53 : : }
54 : :
55 : 955832 : const double nT = sqr_Frobenius( T );
56 : 955832 : const double f = 1 / ( tau * tau );
57 : 955832 : result = ( 1 + f ) * nT - 4;
58 : 955832 : return true;
59 : : }
60 : :
61 : 0 : bool TShapeSizeB1::evaluate( const MsqMatrix< 3, 3 >& T, double& result, MsqError& /*err*/ )
62 : : {
63 : 0 : const double tau = det( T );
64 [ # # ]: 0 : if( invalid_determinant( tau ) )
65 : : { // barrier
66 : 0 : return false;
67 : : }
68 : :
69 : 0 : const double nT = sqr_Frobenius( T );
70 [ # # ]: 0 : const double nadj = sqr_Frobenius( transpose_adj( T ) );
71 : 0 : const double f = 1 / ( tau * tau );
72 : 0 : result = nT + f * nadj - 6;
73 : 0 : return true;
74 : : }
75 : :
76 : 4965548 : bool TShapeSizeB1::evaluate_with_grad( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& deriv_wrt_T,
77 : : MsqError& err )
78 : : {
79 [ + - ]: 4965548 : const double tau = det( T );
80 [ + - ][ - + ]: 4965548 : if( invalid_determinant( tau ) )
81 : : { // barrier
82 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
83 : 0 : return false;
84 : : }
85 : :
86 [ + - ]: 4965548 : const MsqMatrix< 2, 2 > adjt = transpose_adj( T );
87 [ + - ]: 4965548 : const double nT = sqr_Frobenius( T );
88 : 4965548 : const double f = 1 / ( tau * tau );
89 : 4965548 : result = ( 1 + f ) * nT - 4;
90 : :
91 : 4965548 : deriv_wrt_T = T;
92 [ + - ]: 4965548 : deriv_wrt_T *= 2 + 2 * f;
93 [ + - ][ + - ]: 4965548 : deriv_wrt_T -= 2 * f / tau * nT * adjt;
94 : :
95 : 4965548 : return true;
96 : : }
97 : :
98 : 0 : bool TShapeSizeB1::evaluate_with_grad( const MsqMatrix< 3, 3 >& T, double& result, MsqMatrix< 3, 3 >& deriv_wrt_T,
99 : : MsqError& err )
100 : : {
101 [ # # ]: 0 : const double tau = det( T );
102 [ # # ][ # # ]: 0 : if( invalid_determinant( tau ) )
103 : : { // barrier
104 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
105 : 0 : return false;
106 : : }
107 : :
108 [ # # ]: 0 : const MsqMatrix< 3, 3 > adjt = transpose_adj( T );
109 [ # # ]: 0 : const double nT = sqr_Frobenius( T );
110 [ # # ]: 0 : const double nadj = sqr_Frobenius( adjt );
111 : 0 : const double f = 1 / ( tau * tau );
112 : 0 : result = nT + f * nadj - 6;
113 : :
114 : 0 : deriv_wrt_T = T;
115 [ # # ]: 0 : deriv_wrt_T *= ( 1 + f * nT );
116 [ # # ][ # # ]: 0 : deriv_wrt_T -= f * T * transpose( T ) * T;
[ # # ][ # # ]
[ # # ]
117 [ # # ][ # # ]: 0 : deriv_wrt_T -= f / tau * nadj * adjt;
118 [ # # ]: 0 : deriv_wrt_T *= 2;
119 : :
120 : 0 : return true;
121 : : }
122 : :
123 : 0 : bool TShapeSizeB1::evaluate_with_hess( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& deriv_wrt_T,
124 : : MsqMatrix< 2, 2 > second_wrt_T[3], MsqError& err )
125 : : {
126 [ # # ]: 0 : const double tau = det( T );
127 [ # # ][ # # ]: 0 : if( invalid_determinant( tau ) )
128 : : { // barrier
129 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
130 : 0 : return false;
131 : : }
132 : :
133 [ # # ]: 0 : const MsqMatrix< 2, 2 > adjt = transpose_adj( T );
134 [ # # ]: 0 : const double nT = sqr_Frobenius( T );
135 : 0 : const double f = 1 / ( tau * tau );
136 : 0 : result = ( 1 + f ) * nT - 4;
137 : :
138 : 0 : deriv_wrt_T = T;
139 [ # # ]: 0 : deriv_wrt_T *= 2 + 2 * f;
140 [ # # ][ # # ]: 0 : deriv_wrt_T -= 2 * f / tau * nT * adjt;
141 : :
142 [ # # ]: 0 : set_scaled_sum_outer_product( second_wrt_T, -4 * f / tau, T, adjt );
143 [ # # ]: 0 : pluseq_scaled_I( second_wrt_T, 2 + 2 * f );
144 [ # # ]: 0 : pluseq_scaled_outer_product( second_wrt_T, 6 * nT * f * f, adjt );
145 [ # # ]: 0 : pluseq_scaled_2nd_deriv_of_det( second_wrt_T, -2 * nT * f / tau );
146 : :
147 : 0 : return true;
148 : : }
149 : 0 : bool TShapeSizeB1::evaluate_with_hess( const MsqMatrix< 3, 3 >& T, double& result, MsqMatrix< 3, 3 >& deriv_wrt_T,
150 : : MsqMatrix< 3, 3 > second_wrt_T[6], MsqError& err )
151 : : {
152 [ # # ]: 0 : const double tau = det( T );
153 [ # # ][ # # ]: 0 : if( invalid_determinant( tau ) )
154 : : { // barrier
155 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
156 : 0 : return false;
157 : : }
158 : :
159 [ # # ]: 0 : const MsqMatrix< 3, 3 > adjt = transpose_adj( T );
160 [ # # ]: 0 : const double nT = sqr_Frobenius( T );
161 [ # # ]: 0 : const double nadj = sqr_Frobenius( adjt );
162 : 0 : const double f = 1 / ( tau * tau );
163 : 0 : result = nT + f * nadj - 6;
164 : :
165 : : //! \f$ \frac{\partial}{\partial T} |adj T|^2 \f$
166 [ # # ][ # # ]: 0 : const MsqMatrix< 3, 3 > dNadj_dT = 2 * ( nT * T - T * transpose( T ) * T );
[ # # ][ # # ]
[ # # ][ # # ]
167 : 0 : deriv_wrt_T = T;
168 [ # # ][ # # ]: 0 : deriv_wrt_T -= f / tau * nadj * adjt;
169 [ # # ]: 0 : deriv_wrt_T *= 2;
170 [ # # ][ # # ]: 0 : deriv_wrt_T += f * dNadj_dT;
171 : :
172 : : // calculate negative of 2nd wrt T of (|adj T|^2 / tau^2) (sec 3.2.2)
173 [ # # ]: 0 : set_scaled_2nd_deriv_norm_sqr_adj( second_wrt_T, f, T );
174 [ # # ]: 0 : pluseq_scaled_2nd_deriv_of_det( second_wrt_T, -2 * f * f * nadj * tau, T );
175 [ # # ]: 0 : pluseq_scaled_outer_product( second_wrt_T, 6 * f * f * nadj, adjt );
176 [ # # ]: 0 : pluseq_scaled_sum_outer_product( second_wrt_T, -2 * f * f * tau, adjt, dNadj_dT );
177 : : // calculate 2nd wrt T of this metric
178 [ # # ]: 0 : pluseq_scaled_I( second_wrt_T, 2.0 );
179 : :
180 : 0 : return true;
181 : : }
182 : :
183 [ + - ][ + - ]: 16 : } // namespace MBMesquite
|
