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