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 TShapeSize2DNB1.cpp
28 : : * \brief
29 : : * \author Jason Kraftcheck
30 : : */
31 : :
32 : : #include "Mesquite.hpp"
33 : : #include "TShapeSize2DNB1.hpp"
34 : : #include "MsqMatrix.hpp"
35 : : #include "TMPDerivs.hpp"
36 : :
37 : : #include <iostream>
38 : :
39 : : namespace MBMesquite
40 : : {
41 : :
42 : 6 : std::string TShapeSize2DNB1::get_name() const
43 : : {
44 [ + - ]: 6 : return "TShapeSize2DNB1";
45 : : }
46 : :
47 [ - + ]: 46 : TShapeSize2DNB1::~TShapeSize2DNB1() {}
48 : :
49 : : /** \f$ \mu(T) = |T|^2 - 2 \psi(T) + 2 \f$
50 : : * \f$ \psi(T) = \sqrt{|T|^2 + 2 \tau} \f$
51 : : * \f$ \tau = det(T) \f$
52 : : */
53 : 927392 : bool TShapeSize2DNB1::evaluate( const MsqMatrix< 2, 2 >& T, double& result, MsqError& /*err*/ )
54 : : {
55 [ + - ]: 927392 : double frob_sqr = sqr_Frobenius( T );
56 [ + - ]: 927392 : double psi = sqrt( frob_sqr + 2.0 * det( T ) );
57 : :
58 : 927392 : MsqMatrix< 2, 2 > Tdelta( T );
59 [ - + ]: 927392 : while( fabs( psi ) < DBL_EPSILON )
60 : : {
61 [ # # ]: 0 : Tdelta( 0, 0 ) += 1e-12;
62 [ # # ]: 0 : Tdelta( 1, 1 ) += 1e-12;
63 [ # # ]: 0 : frob_sqr = sqr_Frobenius( Tdelta );
64 [ # # ]: 0 : psi = sqrt( frob_sqr + 2.0 * det( Tdelta ) );
65 : : }
66 : :
67 : 927392 : result = frob_sqr - 2.0 * psi + 2.0;
68 : 927392 : return true;
69 : : }
70 : :
71 : 462976 : bool TShapeSize2DNB1::evaluate_with_grad( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& deriv_wrt_T,
72 : : MsqError& /*err*/ )
73 : : {
74 [ + - ]: 462976 : double frob_sqr = sqr_Frobenius( T );
75 [ + - ]: 462976 : double psi = sqrt( frob_sqr + 2.0 * det( T ) );
76 : :
77 : 462976 : MsqMatrix< 2, 2 > Tdelta( T );
78 [ - + ]: 462976 : while( fabs( psi ) < DBL_EPSILON )
79 : : {
80 [ # # ]: 0 : Tdelta( 0, 0 ) += 1e-12;
81 [ # # ]: 0 : Tdelta( 1, 1 ) += 1e-12;
82 [ # # ]: 0 : frob_sqr = sqr_Frobenius( Tdelta );
83 [ # # ]: 0 : psi = sqrt( frob_sqr + 2.0 * det( Tdelta ) );
84 : : }
85 : :
86 : 462976 : result = frob_sqr - 2.0 * psi + 2.0;
87 : :
88 : 462976 : deriv_wrt_T = T;
89 [ + - ]: 462976 : if( psi > 1e-50 )
90 : : {
91 [ + - ]: 462976 : deriv_wrt_T *= ( 1.0 - 1.0 / psi );
92 [ + - ][ + - ]: 462976 : deriv_wrt_T -= 1.0 / psi * transpose_adj( T );
[ + - ]
93 [ + - ]: 462976 : deriv_wrt_T *= 2;
94 : : }
95 : : else
96 : : {
97 [ # # ][ # # ]: 0 : std::cout << "Warning: Division by zero avoided in TShapeSize2DNB2::evaluate_with_grad()" << std::endl;
98 : : }
99 : :
100 : 462976 : return true;
101 : : }
102 : :
103 : 0 : bool TShapeSize2DNB1::evaluate_with_hess( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& deriv_wrt_T,
104 : : MsqMatrix< 2, 2 > second[3], MsqError& /*err*/ )
105 : : {
106 [ # # ]: 0 : double frob_sqr = sqr_Frobenius( T );
107 [ # # ]: 0 : double psi = sqrt( frob_sqr + 2.0 * det( T ) );
108 : :
109 : 0 : MsqMatrix< 2, 2 > Tdelta( T );
110 [ # # ]: 0 : while( fabs( psi ) < DBL_EPSILON )
111 : : {
112 [ # # ]: 0 : Tdelta( 0, 0 ) += 1e-12;
113 [ # # ]: 0 : Tdelta( 1, 1 ) += 1e-12;
114 [ # # ]: 0 : frob_sqr = sqr_Frobenius( Tdelta );
115 [ # # ]: 0 : psi = sqrt( frob_sqr + 2.0 * det( Tdelta ) );
116 : : }
117 : :
118 : 0 : result = frob_sqr - 2.0 * psi + 2.0;
119 : :
120 : 0 : deriv_wrt_T = T;
121 [ # # ]: 0 : if( psi > 1e-50 )
122 : : {
123 [ # # ]: 0 : deriv_wrt_T *= ( 1.0 - 1.0 / psi );
124 [ # # ][ # # ]: 0 : deriv_wrt_T -= 1.0 / psi * transpose_adj( T );
[ # # ]
125 [ # # ]: 0 : deriv_wrt_T *= 2;
126 : : }
127 : : else
128 : : {
129 [ # # ][ # # ]: 0 : std::cout << "Warning: Division by zero avoided in TShapeSize2DNB2::evaluate_with_hess()" << std::endl;
130 : : }
131 : :
132 [ # # ]: 0 : set_scaled_2nd_deriv_wrt_psi( second, -2.0, psi, T );
133 [ # # ]: 0 : pluseq_scaled_I( second, 2 );
134 : :
135 : 0 : return true;
136 : : }
137 : :
138 [ + - ][ + - ]: 20 : } // namespace MBMesquite
|
