Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2007 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 : : (2007) [email protected]
24 : :
25 : : ***************************************************************** */
26 : :
27 : : /** \file TargetMetricUtil.hpp
28 : : * \brief A collection of utility code used by QualtiyMetrics
29 : : * composed of TMP Target Metrics
30 : : * \author Jason Kraftcheck
31 : : */
32 : :
33 : : #ifndef MSQ_TARGET_METRIC_UTIL_HPP
34 : : #define MSQ_TARGET_METRIC_UTIL_HPP
35 : :
36 : : #include "Mesquite.hpp"
37 : : #include "SymMatrix3D.hpp"
38 : : #include <vector>
39 : : #include <assert.h>
40 : :
41 : : namespace MBMesquite
42 : : {
43 : :
44 : : template < unsigned R, unsigned C >
45 : : class MsqMatrix;
46 : : template < unsigned C >
47 : : class MsqVector;
48 : : class PatchData;
49 : : class MsqError;
50 : : class Vector3D;
51 : : class Matrix3D;
52 : :
53 : : /**\brief Calculate R and Z such that \f$W\prime = Z^{-1} W\f$ and
54 : : * \f$A\prime = (RZ)^{-1} A\f$
55 : : *
56 : : * Calculate the matrices required to transform the active and target
57 : : * matrices from the 3x2 surface domain to a 2x2 2D domain.
58 : : *\param A Input: Element Jacobian matrix.
59 : : *\param W_32 Input: Target Jacobian matrix.
60 : : *\param W_22 Output: 2D Target matrix.
61 : : *\param RZ Output: Product of R and Z needed to calculate the 2D
62 : : * element matrix.
63 : : */
64 : : void surface_to_2d( const MsqMatrix< 3, 2 >& A, const MsqMatrix< 3, 2 >& W_32, MsqMatrix< 2, 2 >& W_22,
65 : : MsqMatrix< 3, 2 >& RZ );
66 : : /*
67 : : void surface_to_2d( const MsqMatrix<3,2>& A_in,
68 : : const MsqMatrix<3,2>& W_in,
69 : : MsqMatrix<2,2>& A_out,
70 : : MsqMatrix<2,2>& W_out );
71 : : */
72 : : void get_sample_pt_evaluations( PatchData& pd, std::vector< size_t >& handles, bool free, MsqError& err );
73 : :
74 : : void get_elem_sample_points( PatchData& pd, size_t elem, std::vector< size_t >& handles, MsqError& err );
75 : :
76 : : /**\brief Calculate gradient from derivatives of mapping function terms
77 : : * and derivatives of target metric. */
78 : : template < int DIM >
79 : 7715396 : inline void gradient( size_t num_free_verts, const MsqVector< DIM >* dNdxi, const MsqMatrix< 3, DIM >& dmdA,
80 : : std::vector< Vector3D >& grad )
81 : : {
82 : 7715396 : grad.clear();
83 [ + - + - ]: 7715396 : grad.resize( num_free_verts, Vector3D( 0, 0, 0 ) );
84 [ + + ][ + + ]: 27088663 : for( size_t i = 0; i < num_free_verts; ++i )
85 [ + - ][ + - ]: 19373267 : grad[i] = Vector3D( ( dmdA * dNdxi[i] ).data() );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
86 : 7715396 : }
87 : :
88 : : /**\brief Calculate Hessian from derivatives of mapping function terms
89 : : * and derivatives of target metric. */
90 : : template < int DIM, typename MAT >
91 : 236360 : inline void hessian( size_t num_free_verts, const MsqVector< DIM >* dNdxi, const MsqMatrix< DIM, DIM >* d2mdA2,
92 : : MAT* hess )
93 : : {
94 [ + - ][ + + ]: 3072680 : MsqMatrix< 1, DIM > tmp[DIM][DIM];
[ + + ]
95 : 236360 : size_t h = 0; // index of current Hessian block
96 : :
97 [ + + ]: 992740 : for( size_t i = 0; i < num_free_verts; ++i )
98 : : {
99 : :
100 : : // Populate TMP with vector-matrix procucts common
101 : : // to terms of this Hessian row.
102 [ + - ]: 756380 : const MsqMatrix< 1, DIM >& gi = transpose( dNdxi[i] );
103 : : switch( DIM )
104 : : {
105 : : case 3:
106 [ + - ]: 756380 : tmp[0][2] = gi * d2mdA2[2];
107 [ + - ]: 756380 : tmp[1][2] = gi * d2mdA2[4];
108 [ + - ][ + - ]: 756380 : tmp[2][0] = gi * transpose( d2mdA2[2] );
109 [ + - ][ + - ]: 756380 : tmp[2][1] = gi * transpose( d2mdA2[4] );
110 [ + - ]: 756380 : tmp[2][2] = gi * d2mdA2[5];
111 : : case 2:
112 [ + - ]: 756380 : tmp[0][1] = gi * d2mdA2[1];
113 [ + - ][ + - ]: 756380 : tmp[1][0] = gi * transpose( d2mdA2[1] );
114 [ + - ]: 756380 : tmp[1][1] = gi * d2mdA2[DIM];
115 : : case 1:
116 [ + - ]: 756380 : tmp[0][0] = gi * d2mdA2[0];
117 : : case 0:
118 : 756380 : break;
119 : : default:
120 : : assert( false );
121 : : }
122 : :
123 : : // Calculate Hessian diagonal block
124 : 756380 : MAT& H = hess[h++];
125 : : switch( DIM )
126 : : {
127 : : case 3:
128 [ + - ][ + - ]: 756380 : H( 0, 2 ) = H( 2, 0 ) = tmp[0][2] * transpose( gi );
[ + - ][ + - ]
[ + - ]
129 [ + - ][ + - ]: 756380 : H( 1, 2 ) = H( 2, 1 ) = tmp[1][2] * transpose( gi );
[ + - ][ + - ]
[ + - ]
130 [ + - ][ + - ]: 756380 : H( 2, 2 ) = tmp[2][2] * transpose( gi );
[ + - ][ + - ]
131 : : case 2:
132 [ + - ][ + - ]: 756380 : H( 0, 1 ) = H( 1, 0 ) = tmp[0][1] * transpose( gi );
[ + - ][ + - ]
[ + - ]
133 [ + - ][ + - ]: 756380 : H( 1, 1 ) = tmp[1][1] * transpose( gi );
[ + - ][ + - ]
134 : : case 1:
135 [ + - ][ + - ]: 756380 : H( 0, 0 ) = tmp[0][0] * transpose( gi );
[ + - ][ + - ]
136 : : case 0:
137 : 756380 : break;
138 : : default:
139 : : assert( false );
140 : : }
141 : :
142 : : // Calculate remainder of Hessian row
143 [ + + ]: 1731630 : for( size_t j = i + 1; j < num_free_verts; ++j )
144 : : {
145 : 975250 : MAT& HH = hess[h++];
146 : 975250 : const MsqMatrix< DIM, 1 >& gj = dNdxi[j];
147 : : switch( DIM )
148 : : {
149 : : case 3:
150 [ + - ][ + - ]: 975250 : HH( 0, 2 ) = tmp[0][2] * gj;
[ + - ]
151 [ + - ][ + - ]: 975250 : HH( 1, 2 ) = tmp[1][2] * gj;
[ + - ]
152 [ + - ][ + - ]: 975250 : HH( 2, 0 ) = tmp[2][0] * gj;
[ + - ]
153 [ + - ][ + - ]: 975250 : HH( 2, 1 ) = tmp[2][1] * gj;
[ + - ]
154 [ + - ][ + - ]: 975250 : HH( 2, 2 ) = tmp[2][2] * gj;
[ + - ]
155 : : case 2:
156 [ + - ][ + - ]: 975250 : HH( 0, 1 ) = tmp[0][1] * gj;
[ + - ]
157 [ + - ][ + - ]: 975250 : HH( 1, 0 ) = tmp[1][0] * gj;
[ + - ]
158 [ + - ][ + - ]: 975250 : HH( 1, 1 ) = tmp[1][1] * gj;
[ + - ]
159 : : case 1:
160 [ + - ][ + - ]: 975250 : HH( 0, 0 ) = tmp[0][0] * gj;
[ + - ]
161 : : case 0:
162 : 975250 : break;
163 : : default:
164 : : assert( false );
165 : : }
166 : : }
167 : : }
168 : 236360 : }
169 : :
170 : : /**\brief Calculate Hessian from derivatives of mapping function terms
171 : : * and derivatives of target metric. */
172 : : template < int DIM >
173 : 0 : inline void hessian_diagonal( size_t num_free_verts, const MsqVector< DIM >* dNdxi, const MsqMatrix< DIM, DIM >* d2mdA2,
174 : : SymMatrix3D* diagonal )
175 : : {
176 [ # # ]: 0 : for( size_t i = 0; i < num_free_verts; ++i )
177 : : {
178 : 0 : SymMatrix3D& H = diagonal[i];
179 [ # # ]: 0 : for( unsigned j = 0; j < ( ( DIM ) * ( DIM + 1 ) / 2 ); ++j )
180 [ # # ][ # # ]: 0 : H[j] = transpose( dNdxi[i] ) * d2mdA2[j] * dNdxi[i];
[ # # ]
181 : : }
182 : 0 : }
183 : :
184 : : #ifdef PRINT_INFO
185 : : template < int R, int C >
186 : : inline void write_vect( char n, const MsqMatrix< R, C >& M )
187 : : {
188 : : std::cout << " " << n << ':';
189 : : for( int c = 0; c < C; ++c )
190 : : {
191 : : std::cout << '[';
192 : : for( int r = 0; r < R; ++r )
193 : : std::cout << M( r, c ) << ' ';
194 : : std::cout << ']';
195 : : }
196 : : std::cout << std::endl;
197 : : }
198 : :
199 : : template < int D >
200 : : inline void print_info( size_t elem, Sample sample, const MsqMatrix< 3, D >& A, const MsqMatrix< 3, D >& W,
201 : : const MsqMatrix< D, D >& T )
202 : : {
203 : : std::cout << "Elem " << elem << " Dim " << sample.dimension << " Num " << sample.number << " :" << std::endl;
204 : : write_vect< 3, D >( 'A', A );
205 : : write_vect< 3, D >( 'W', W );
206 : : write_vect< D, D >( 'T', T );
207 : : }
208 : : #endif
209 : :
210 : : } // namespace MBMesquite
211 : :
212 : : #endif
|
