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 AffineMapMetric.cpp
28 : : * \brief
29 : : * \author Jason Kraftcheck
30 : : */
31 : :
32 : : #include "Mesquite.hpp"
33 : : #include "AffineMapMetric.hpp"
34 : : #include "MsqMatrix.hpp"
35 : : #include "ElementQM.hpp"
36 : : #include "MsqError.hpp"
37 : : #include "Vector3D.hpp"
38 : : #include "PatchData.hpp"
39 : : #include "MappingFunction.hpp"
40 : : #include "WeightCalculator.hpp"
41 : : #include "TargetCalculator.hpp"
42 : : #include "TMetric.hpp"
43 : : #include "TargetMetricUtil.hpp"
44 : :
45 : : #include <functional>
46 : : #include <algorithm>
47 : :
48 : : namespace MBMesquite
49 : : {
50 : :
51 : : const double TRI_XFORM_VALS[] = { 1.0, -1.0 / sqrt( 3.0 ), 0.0, 2.0 / sqrt( 3.0 ) };
52 : 1 : MsqMatrix< 2, 2 > TRI_XFORM( TRI_XFORM_VALS );
53 : :
54 : : const double TET_XFORM_VALS[] = {
55 : : 1.0, -1.0 / sqrt( 3.0 ), -1.0 / sqrt( 6.0 ), 0.0, 2.0 / sqrt( 3.0 ), -1.0 / sqrt( 6.0 ), 0.0, 0.0, sqrt( 3.0 / 2.0 )
56 : : };
57 : 1 : MsqMatrix< 3, 3 > TET_XFORM( TET_XFORM_VALS );
58 : :
59 : 0 : AffineMapMetric::AffineMapMetric( TargetCalculator* tc, WeightCalculator* wc, TMetric* target_metric )
60 : 0 : : targetCalc( tc ), weightCalc( wc ), targetMetric( target_metric )
61 : : {
62 : 0 : }
63 : :
64 : 9 : AffineMapMetric::AffineMapMetric( TargetCalculator* tc, TMetric* target_metric )
65 : 9 : : targetCalc( tc ), weightCalc( 0 ), targetMetric( target_metric )
66 : : {
67 : 9 : }
68 : :
69 : 6281 : int AffineMapMetric::get_negate_flag() const
70 : : {
71 : 6281 : return 1;
72 : : }
73 : :
74 : 0 : std::string AffineMapMetric::get_name() const
75 : : {
76 [ # # ][ # # ]: 0 : return std::string( "AffineMap(" ) + targetMetric->get_name() + ')';
[ # # ]
77 : : }
78 : :
79 : 5720 : void AffineMapMetric::get_evaluations( PatchData& pd, std::vector< size_t >& handles, bool free, MsqError& err )
80 : : {
81 : 5720 : get_sample_pt_evaluations( pd, handles, free, err );
82 : 5720 : }
83 : :
84 : 43008 : void AffineMapMetric::get_element_evaluations( PatchData& pd, size_t p_elem, std::vector< size_t >& handles,
85 : : MsqError& err )
86 : : {
87 : 43008 : get_elem_sample_points( pd, p_elem, handles, err );
88 : 43008 : }
89 : :
90 : 2596324 : bool AffineMapMetric::evaluate( PatchData& pd, size_t p_handle, double& value, MsqError& err )
91 : : {
92 [ + - ]: 2596324 : Sample s = ElemSampleQM::sample( p_handle );
93 [ + - ]: 2596324 : size_t e = ElemSampleQM::elem( p_handle );
94 [ + - ]: 2596324 : MsqMeshEntity& p_elem = pd.element_by_index( e );
95 [ + - ]: 2596324 : EntityTopology type = p_elem.get_element_type();
96 [ + - ]: 2596324 : unsigned edim = TopologyInfo::dimension( type );
97 [ + - ]: 2596324 : const size_t* conn = p_elem.get_vertex_index_array();
98 : :
99 : : // This metric only supports sampling at corners, except for simplices.
100 : : // If element is a simpex, then the Jacobian is constant over a linear
101 : : // element. In this case, always evaluate at any vertex.
102 : : // unsigned corner = s.number;
103 [ + + ]: 2596324 : if( s.dimension != 0 )
104 : : {
105 [ - + ][ # # ]: 1735923 : if( type == TRIANGLE || type == TETRAHEDRON )
106 : : /*corner = 0*/;
107 : : else
108 : : {
109 : : MSQ_SETERR( err )
110 [ # # ][ # # ]: 0 : ( "Invalid sample point for AffineMapMetric", MsqError::UNSUPPORTED_ELEMENT );
111 : 0 : return false;
112 : : }
113 : : }
114 : :
115 : : bool rval;
116 [ - + ]: 2596324 : if( edim == 3 )
117 : : { // 3x3 or 3x2 targets ?
118 [ # # ][ # # ]: 0 : Vector3D c[3] = { Vector3D( 0, 0, 0 ), Vector3D( 0, 0, 0 ), Vector3D( 0, 0, 0 ) };
[ # # ]
119 : : unsigned n;
120 [ # # ]: 0 : const unsigned* adj = TopologyInfo::adjacent_vertices( type, s.number, n );
121 [ # # ][ # # ]: 0 : c[0] = pd.vertex_by_index( conn[adj[0]] ) - pd.vertex_by_index( conn[s.number] );
[ # # ][ # # ]
122 [ # # ][ # # ]: 0 : c[1] = pd.vertex_by_index( conn[adj[1]] ) - pd.vertex_by_index( conn[s.number] );
[ # # ][ # # ]
123 [ # # ][ # # ]: 0 : c[2] = pd.vertex_by_index( conn[adj[2]] ) - pd.vertex_by_index( conn[s.number] );
[ # # ][ # # ]
124 [ # # ]: 0 : MsqMatrix< 3, 3 > A;
125 [ # # ][ # # ]: 0 : A.set_column( 0, MsqMatrix< 3, 1 >( c[0].to_array() ) );
[ # # ]
126 [ # # ][ # # ]: 0 : A.set_column( 1, MsqMatrix< 3, 1 >( c[1].to_array() ) );
[ # # ]
127 [ # # ][ # # ]: 0 : A.set_column( 2, MsqMatrix< 3, 1 >( c[2].to_array() ) );
[ # # ]
128 [ # # ][ # # ]: 0 : if( type == TETRAHEDRON ) A = A * TET_XFORM;
129 : :
130 [ # # ]: 0 : MsqMatrix< 3, 3 > W;
131 [ # # ]: 0 : targetCalc->get_3D_target( pd, e, s, W, err );
132 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
133 [ # # ][ # # ]: 0 : rval = targetMetric->evaluate( A * inverse( W ), value, err );
[ # # ]
134 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
135 : : }
136 : : else
137 : : {
138 [ + - ][ + - ]: 2596324 : Vector3D c[2] = { Vector3D( 0, 0, 0 ), Vector3D( 0, 0, 0 ) };
139 : : unsigned n;
140 [ + - ]: 2596324 : const unsigned* adj = TopologyInfo::adjacent_vertices( type, s.number, n );
141 [ + - ][ + - ]: 2596324 : c[0] = pd.vertex_by_index( conn[adj[0]] ) - pd.vertex_by_index( conn[s.number] );
[ + - ][ + - ]
142 [ + - ][ + - ]: 2596324 : c[1] = pd.vertex_by_index( conn[adj[1]] ) - pd.vertex_by_index( conn[s.number] );
[ + - ][ + - ]
143 [ + - ]: 2596324 : MsqMatrix< 3, 2 > App;
144 [ + - ][ + - ]: 2596324 : App.set_column( 0, MsqMatrix< 3, 1 >( c[0].to_array() ) );
[ + - ]
145 [ + - ][ + - ]: 2596324 : App.set_column( 1, MsqMatrix< 3, 1 >( c[1].to_array() ) );
[ + - ]
146 : :
147 [ + - ]: 2596324 : MsqMatrix< 3, 2 > Wp;
148 [ + - ]: 2596324 : targetCalc->get_surface_target( pd, e, s, Wp, err );
149 [ + - ][ - + ]: 2596324 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
150 : :
151 [ + - ][ + - ]: 2596324 : MsqMatrix< 2, 2 > A, W;
152 [ + - ]: 2596324 : MsqMatrix< 3, 2 > RZ;
153 [ + - ]: 2596324 : surface_to_2d( App, Wp, W, RZ );
154 [ + - ][ + - ]: 2596324 : A = transpose( RZ ) * App;
155 [ + + ][ + - ]: 2596324 : if( type == TRIANGLE ) A = A * TRI_XFORM;
156 : :
157 [ + - ][ + - ]: 2596324 : rval = targetMetric->evaluate( A * inverse( W ), value, err );
[ + - ]
158 [ + - ][ - + ]: 2596324 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
159 : : }
160 : :
161 : : // apply target weight to value
162 [ + + ]: 2596324 : if( weightCalc )
163 : : {
164 [ + - ]: 370575 : double ck = weightCalc->get_weight( pd, e, s, err );
165 [ + - ][ - + ]: 370575 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
166 : 370575 : value *= ck;
167 : : }
168 : 2596324 : return rval;
169 : : }
170 : :
171 : 110331 : bool AffineMapMetric::evaluate_with_indices( PatchData& pd, size_t p_handle, double& value,
172 : : std::vector< size_t >& indices, MsqError& err )
173 : : {
174 [ + - ]: 110331 : Sample s = ElemSampleQM::sample( p_handle );
175 [ + - ]: 110331 : size_t e = ElemSampleQM::elem( p_handle );
176 [ + - ]: 110331 : MsqMeshEntity& p_elem = pd.element_by_index( e );
177 [ + - ]: 110331 : EntityTopology type = p_elem.get_element_type();
178 [ + - ]: 110331 : const size_t* conn = p_elem.get_vertex_index_array();
179 : :
180 : : // this metric only supports sampling at corners
181 [ + + ]: 110331 : if( s.dimension != 0 )
182 : : {
183 [ - + ][ # # ]: 73431 : if( type != TRIANGLE && type != TETRAHEDRON )
184 : : {
185 : : MSQ_SETERR( err )
186 [ # # ][ # # ]: 0 : ( "Invalid sample point for AffineMapMetric", MsqError::UNSUPPORTED_ELEMENT );
187 : 0 : return false;
188 : : }
189 : 73431 : s.dimension = 0;
190 : 73431 : s.number = 0;
191 : : }
192 : :
193 : : unsigned n;
194 [ + - ]: 110331 : const unsigned* adj = TopologyInfo::adjacent_vertices( type, s.number, n );
195 : 110331 : indices.clear();
196 [ + - ][ + + ]: 110331 : if( conn[s.number] < pd.num_free_vertices() ) indices.push_back( conn[s.number] );
[ + - ]
197 [ + + ]: 330993 : for( unsigned i = 0; i < n; ++i )
198 [ + - ][ + + ]: 220662 : if( conn[adj[i]] < pd.num_free_vertices() ) indices.push_back( conn[adj[i]] );
[ + - ]
199 : :
200 [ + - ]: 110331 : return evaluate( pd, p_handle, value, err );
201 : : }
202 : :
203 [ + - ][ + - ]: 4 : } // namespace MBMesquite
|
