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 : : (2010) [email protected]
25 : :
26 : : ***************************************************************** */
27 : :
28 : : /** \file TShapeOrientB1.cpp
29 : : * \brief
30 : : * \author Jason Kraftcheck
31 : : */
32 : :
33 : : #include "Mesquite.hpp"
34 : : #include "TShapeOrientB1.hpp"
35 : : #include "MsqMatrix.hpp"
36 : : #include "MsqError.hpp"
37 : : #include "TMPDerivs.hpp"
38 : : #include "TMPCommon.hpp"
39 : :
40 : : namespace MBMesquite
41 : : {
42 : :
43 : 0 : std::string TShapeOrientB1::get_name() const
44 : : {
45 [ # # ]: 0 : return "TShapeOrientB1";
46 : : }
47 : :
48 [ # # ]: 0 : TShapeOrientB1::~TShapeOrientB1() {}
49 : :
50 : 15446 : bool TShapeOrientB1::evaluate( const MsqMatrix< 2, 2 >& T, double& result, MsqError& err )
51 : : {
52 : 15446 : const double tau = det( T );
53 [ + + ]: 15446 : if( TMetric::invalid_determinant( tau ) )
54 : : { // barrier
55 [ + - ]: 1 : MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
56 : 1 : return false;
57 : : }
58 : 15445 : result = 0.5 / tau * ( Frobenius( T ) - trace( T ) / MSQ_SQRT_TWO );
59 : 15446 : return true;
60 : : }
61 : :
62 : 1869 : bool TShapeOrientB1::evaluate_with_grad( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& deriv_wrt_T,
63 : : MsqError& err )
64 : : {
65 : 1869 : const double norm = Frobenius( T );
66 : 1869 : const double invroot = 1.0 / MSQ_SQRT_TWO;
67 : 1869 : const double tau = det( T );
68 [ + + ]: 1869 : if( TMetric::invalid_determinant( tau ) )
69 : : { // barrier
70 [ + - ]: 1 : MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
71 : 1 : return false;
72 : : }
73 : 1868 : const double inv_tau = 1.0 / tau;
74 : 1868 : const double invnorm = 1.0 / norm;
75 : :
76 : 1868 : result = 0.5 * inv_tau * ( norm - invroot * trace( T ) );
77 : :
78 : 1868 : deriv_wrt_T = invnorm * T;
79 : 1868 : pluseq_scaled_I( deriv_wrt_T, -invroot );
80 : 1868 : deriv_wrt_T *= 0.5;
81 [ + - ][ + - ]: 1868 : deriv_wrt_T -= result * transpose_adj( T );
82 : 1868 : deriv_wrt_T *= inv_tau;
83 : 1869 : return true;
84 : : }
85 : :
86 : 0 : bool TShapeOrientB1::evaluate_with_hess( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& deriv_wrt_T,
87 : : MsqMatrix< 2, 2 > second_wrt_T[3], MsqError& err )
88 : : {
89 [ # # ]: 0 : const double norm = Frobenius( T );
90 : 0 : const double invroot = 1.0 / MSQ_SQRT_TWO;
91 [ # # ]: 0 : const double tau = det( T );
92 [ # # ][ # # ]: 0 : if( TMetric::invalid_determinant( tau ) )
93 : : { // barrier
94 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
95 : 0 : return false;
96 : : }
97 : 0 : const double inv_tau = 1.0 / tau;
98 : 0 : const double invnorm = 1.0 / norm;
99 : :
100 [ # # ]: 0 : const double f = norm - invroot * trace( T );
101 : 0 : result = 0.5 * inv_tau * f;
102 : :
103 [ # # ]: 0 : const MsqMatrix< 2, 2 > adjt = transpose_adj( T );
104 [ # # ]: 0 : deriv_wrt_T = invnorm * T;
105 [ # # ]: 0 : pluseq_scaled_I( deriv_wrt_T, -invroot );
106 [ # # ]: 0 : deriv_wrt_T *= 0.5;
107 [ # # ][ # # ]: 0 : deriv_wrt_T -= result * adjt;
108 [ # # ]: 0 : deriv_wrt_T *= inv_tau;
109 : :
110 : 0 : const double a = 0.5 * inv_tau * invnorm;
111 [ # # ]: 0 : set_scaled_outer_product( second_wrt_T, -a * invnorm * invnorm, T );
112 [ # # ]: 0 : pluseq_scaled_I( second_wrt_T, a );
113 [ # # ]: 0 : pluseq_scaled_outer_product( second_wrt_T, f * inv_tau * inv_tau * inv_tau, adjt );
114 [ # # ]: 0 : pluseq_scaled_2nd_deriv_of_det( second_wrt_T, -0.5 * f * inv_tau * inv_tau, T );
115 [ # # ]: 0 : pluseq_scaled_sum_outer_product( second_wrt_T, -0.5 * inv_tau * inv_tau * invnorm, T, adjt );
116 [ # # ]: 0 : pluseq_scaled_sum_outer_product_I( second_wrt_T, 0.5 * inv_tau * inv_tau * invroot, adjt );
117 : 0 : return true;
118 : : }
119 : :
120 : 0 : bool TShapeOrientB1::evaluate( const MsqMatrix< 3, 3 >& T, double& result, MsqError& err )
121 : : {
122 : 0 : const double tau = det( T );
123 [ # # ]: 0 : if( TMetric::invalid_determinant( tau ) )
124 : : { // barrier
125 [ # # ]: 0 : MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
126 : 0 : return false;
127 : : }
128 : 0 : result = 0.5 / tau * ( Frobenius( T ) - trace( T ) / MSQ_SQRT_THREE );
129 : 0 : return true;
130 : : }
131 : :
132 : 0 : bool TShapeOrientB1::evaluate_with_grad( const MsqMatrix< 3, 3 >& T, double& result, MsqMatrix< 3, 3 >& deriv_wrt_T,
133 : : MsqError& err )
134 : : {
135 : 0 : const double norm = Frobenius( T );
136 : 0 : const double invroot = 1.0 / MSQ_SQRT_THREE;
137 : 0 : const double tau = det( T );
138 [ # # ]: 0 : if( TMetric::invalid_determinant( tau ) )
139 : : { // barrier
140 [ # # ]: 0 : MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
141 : 0 : return false;
142 : : }
143 : 0 : const double inv_tau = 1.0 / tau;
144 : 0 : const double invnorm = 1.0 / norm;
145 : :
146 : 0 : result = 0.5 * inv_tau * ( norm - invroot * trace( T ) );
147 : :
148 : 0 : deriv_wrt_T = invnorm * T;
149 : 0 : pluseq_scaled_I( deriv_wrt_T, -invroot );
150 : 0 : deriv_wrt_T *= 0.5;
151 [ # # ][ # # ]: 0 : deriv_wrt_T -= result * transpose_adj( T );
152 : 0 : deriv_wrt_T *= inv_tau;
153 : 0 : return true;
154 : : }
155 : :
156 : 0 : bool TShapeOrientB1::evaluate_with_hess( const MsqMatrix< 3, 3 >& T, double& result, MsqMatrix< 3, 3 >& deriv_wrt_T,
157 : : MsqMatrix< 3, 3 > second_wrt_T[6], MsqError& err )
158 : : {
159 [ # # ]: 0 : const double norm = Frobenius( T );
160 : 0 : const double invroot = 1.0 / MSQ_SQRT_THREE;
161 [ # # ]: 0 : const double tau = det( T );
162 [ # # ][ # # ]: 0 : if( TMetric::invalid_determinant( tau ) )
163 : : { // barrier
164 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( barrier_violated_msg, MsqError::BARRIER_VIOLATED );
165 : 0 : return false;
166 : : }
167 : 0 : const double inv_tau = 1.0 / tau;
168 : 0 : const double invnorm = 1.0 / norm;
169 : :
170 [ # # ]: 0 : const double f = norm - invroot * trace( T );
171 : 0 : result = 0.5 * inv_tau * f;
172 : :
173 [ # # ]: 0 : const MsqMatrix< 3, 3 > adjt = transpose_adj( T );
174 [ # # ]: 0 : deriv_wrt_T = invnorm * T;
175 [ # # ]: 0 : pluseq_scaled_I( deriv_wrt_T, -invroot );
176 [ # # ]: 0 : deriv_wrt_T *= 0.5;
177 [ # # ][ # # ]: 0 : deriv_wrt_T -= result * adjt;
178 [ # # ]: 0 : deriv_wrt_T *= inv_tau;
179 : :
180 : 0 : const double a = 0.5 * inv_tau * invnorm;
181 [ # # ]: 0 : set_scaled_outer_product( second_wrt_T, -a * invnorm * invnorm, T );
182 [ # # ]: 0 : pluseq_scaled_I( second_wrt_T, a );
183 [ # # ]: 0 : pluseq_scaled_outer_product( second_wrt_T, f * inv_tau * inv_tau * inv_tau, adjt );
184 [ # # ]: 0 : pluseq_scaled_2nd_deriv_of_det( second_wrt_T, -0.5 * f * inv_tau * inv_tau, T );
185 [ # # ]: 0 : pluseq_scaled_sum_outer_product( second_wrt_T, -0.5 * inv_tau * inv_tau * invnorm, T, adjt );
186 [ # # ]: 0 : pluseq_scaled_sum_outer_product_I( second_wrt_T, 0.5 * inv_tau * inv_tau * invroot, adjt );
187 : 0 : return true;
188 : : }
189 : :
190 [ + - ][ + - ]: 8 : } // namespace MBMesquite
|
