Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2006 Lawrence Livermore National Laboratory. Under
5 : : the terms of Contract B545069 with the University of Wisconsin --
6 : : Madison, Lawrence Livermore National Laboratory retains certain
7 : : rights in 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 AddQualityMetric.cpp
28 : : * \brief
29 : : * \author Jason Kraftcheck
30 : : */
31 : :
32 : : #include "Mesquite.hpp"
33 : : #include "AddQualityMetric.hpp"
34 : : #include "Vector3D.hpp"
35 : : #include "Matrix3D.hpp"
36 : : #include "MsqError.hpp"
37 : :
38 : : namespace MBMesquite
39 : : {
40 : :
41 : 1 : AddQualityMetric::AddQualityMetric( QualityMetric* qm1, QualityMetric* qm2, MsqError& err )
42 [ + - ][ + - ]: 1 : : metric1( *qm1 ), metric2( *qm2 )
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
43 : : {
44 [ + - ][ + - ]: 1 : if( qm1->get_metric_type() != qm2->get_metric_type() || qm1->get_negate_flag() != qm2->get_negate_flag() )
[ + - ][ + - ]
[ + - ][ - + ]
[ - + ]
45 [ # # ][ # # ]: 0 : { MSQ_SETERR( err )( "Incompatible metrics", MsqError::INVALID_ARG ); }
46 : 1 : }
47 : :
48 [ - + ]: 2 : AddQualityMetric::~AddQualityMetric() {}
49 : :
50 : 3 : QualityMetric::MetricType AddQualityMetric::get_metric_type() const
51 : : {
52 : 3 : return metric1.get_metric_type();
53 : : }
54 : :
55 : 1 : std::string AddQualityMetric::get_name() const
56 : : {
57 [ + - ][ + - ]: 1 : return std::string( "Sum(" ) + metric1.get_name() + ", " + metric2.get_name() + ")";
[ + - ][ + - ]
[ + - ][ + - ]
58 : : }
59 : :
60 : 14 : int AddQualityMetric::get_negate_flag() const
61 : : {
62 : 14 : return metric1.get_negate_flag();
63 : : }
64 : :
65 : 47288 : void AddQualityMetric::get_evaluations( PatchData& pd, std::vector< size_t >& handles, bool free_only, MsqError& err )
66 : : {
67 [ - + ][ # # ]: 47288 : metric1.get_evaluations( pd, handles, free_only, err );MSQ_ERRRTN( err );
[ - + ]
68 [ - + ][ # # ]: 47288 : metric2.get_evaluations( pd, mHandles, free_only, err );MSQ_ERRRTN( err );
[ - + ]
69 [ - + ][ # # ]: 47288 : if( handles != mHandles ) { MSQ_SETERR( err )( "Incompatible metrics", MsqError::INVALID_STATE ); }
70 : : }
71 : :
72 : 171574 : bool AddQualityMetric::evaluate( PatchData& pd, size_t handle, double& value, MsqError& err )
73 : : {
74 : : double val1, val2;
75 : : bool rval1, rval2;
76 [ + - ]: 171574 : rval1 = metric1.evaluate( pd, handle, val1, err );
77 [ + - ][ + + ]: 171574 : MSQ_ERRZERO( err );
[ + - ][ + - ]
[ + + ]
78 [ + - ]: 171567 : rval2 = metric2.evaluate( pd, handle, val2, err );
79 [ + - ][ - + ]: 171567 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
80 : 171567 : value = val1 + val2;
81 [ + - ][ + - ]: 171574 : return rval1 && rval2;
82 : : }
83 : :
84 : 0 : bool AddQualityMetric::evaluate_with_indices( PatchData& pd, size_t handle, double& value,
85 : : std::vector< size_t >& indices, MsqError& err )
86 : : {
87 : : double val1, val2;
88 : : bool rval1, rval2;
89 [ # # ]: 0 : rval1 = metric1.evaluate_with_indices( pd, handle, val1, indices1, err );
90 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
91 [ # # ]: 0 : rval2 = metric2.evaluate_with_indices( pd, handle, val2, indices2, err );
92 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
93 : :
94 : 0 : indices.clear();
95 [ # # ]: 0 : std::sort( indices1.begin(), indices1.end() );
96 [ # # ]: 0 : std::sort( indices2.begin(), indices2.end() );
97 [ # # ][ # # ]: 0 : std::set_union( indices1.begin(), indices1.end(), indices2.begin(), indices2.end(), std::back_inserter( indices ) );
98 : :
99 : 0 : value = val1 + val2;
100 [ # # ][ # # ]: 0 : return rval1 && rval2;
101 : : }
102 : :
103 : 0 : bool AddQualityMetric::evaluate_with_gradient( PatchData& pd, size_t handle, double& value,
104 : : std::vector< size_t >& indices, std::vector< Vector3D >& gradient,
105 : : MsqError& err )
106 : : {
107 : 0 : std::vector< size_t >::iterator i;
108 : : size_t j;
109 : : double val1, val2;
110 : : bool rval1, rval2;
111 [ # # ]: 0 : rval1 = metric1.evaluate_with_gradient( pd, handle, val1, indices1, grad1, err );
112 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
113 [ # # ]: 0 : rval2 = metric2.evaluate_with_gradient( pd, handle, val2, indices2, grad2, err );
114 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
115 : :
116 [ # # ]: 0 : indices.resize( indices1.size() + indices2.size() );
117 [ # # ]: 0 : i = std::copy( indices1.begin(), indices1.end(), indices.begin() );
118 [ # # ]: 0 : std::copy( indices2.begin(), indices2.end(), i );
119 [ # # ]: 0 : std::sort( indices.begin(), indices.end() );
120 [ # # ][ # # ]: 0 : indices.erase( std::unique( indices.begin(), indices.end() ), indices.end() );
121 : :
122 : 0 : gradient.clear();
123 [ # # ][ # # ]: 0 : gradient.resize( indices.size(), Vector3D( 0.0 ) );
124 [ # # ]: 0 : for( j = 0; j < indices1.size(); ++j )
125 : : {
126 [ # # ][ # # ]: 0 : i = std::lower_bound( indices.begin(), indices.end(), indices1[j] );
127 [ # # ]: 0 : size_t k = i - indices.begin();
128 [ # # ][ # # ]: 0 : gradient[k] += grad1[j];
[ # # ]
129 : : }
130 [ # # ]: 0 : for( j = 0; j < indices2.size(); ++j )
131 : : {
132 [ # # ][ # # ]: 0 : i = std::lower_bound( indices.begin(), indices.end(), indices2[j] );
133 [ # # ]: 0 : size_t k = i - indices.begin();
134 [ # # ][ # # ]: 0 : gradient[k] += grad2[j];
[ # # ]
135 : : }
136 : :
137 : 0 : value = val1 + val2;
138 [ # # ][ # # ]: 0 : return rval1 && rval2;
139 : : }
140 : :
141 : 118180 : bool AddQualityMetric::evaluate_with_Hessian( PatchData& pd, size_t handle, double& value,
142 : : std::vector< size_t >& indices, std::vector< Vector3D >& gradient,
143 : : std::vector< Matrix3D >& Hessian, MsqError& err )
144 : : {
145 : 118180 : std::vector< size_t >::iterator i;
146 : : size_t j, r, c, n, h;
147 : : double val1, val2;
148 : : bool rval1, rval2;
149 [ + - ]: 118180 : rval1 = metric1.evaluate_with_Hessian( pd, handle, val1, indices1, grad1, Hess1, err );
150 [ + - ][ - + ]: 118180 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
151 [ + - ]: 118180 : rval2 = metric2.evaluate_with_Hessian( pd, handle, val2, indices2, grad2, Hess2, err );
152 [ + - ][ - + ]: 118180 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
153 : :
154 [ + - ]: 118180 : indices.resize( indices1.size() + indices2.size() );
155 [ + - ]: 118180 : i = std::copy( indices1.begin(), indices1.end(), indices.begin() );
156 [ + - ]: 118180 : std::copy( indices2.begin(), indices2.end(), i );
157 [ + - ]: 118180 : std::sort( indices.begin(), indices.end() );
158 [ + - ][ + - ]: 118180 : indices.erase( std::unique( indices.begin(), indices.end() ), indices.end() );
159 : :
160 : 118180 : gradient.clear();
161 [ + - ][ + - ]: 118180 : gradient.resize( indices.size(), Vector3D( 0.0 ) );
162 [ + + ]: 496370 : for( j = 0; j < indices1.size(); ++j )
163 : : {
164 [ + - ][ + - ]: 378190 : i = std::lower_bound( indices.begin(), indices.end(), indices1[j] );
165 [ + - ][ + - ]: 378190 : indices1[j] = i - indices.begin();
166 [ + - ][ + - ]: 378190 : gradient[indices1[j]] += grad1[j];
[ + - ][ + - ]
167 : : }
168 [ + + ]: 496370 : for( j = 0; j < indices2.size(); ++j )
169 : : {
170 [ + - ][ + - ]: 378190 : i = std::lower_bound( indices.begin(), indices.end(), indices2[j] );
171 [ + - ][ + - ]: 378190 : indices2[j] = i - indices.begin();
172 [ + - ][ + - ]: 378190 : gradient[indices2[j]] += grad2[j];
[ + - ][ + - ]
173 : : }
174 : :
175 : 118180 : const size_t N = indices.size();
176 : 118180 : Hessian.clear();
177 [ + - ][ + - ]: 118180 : Hessian.resize( N * ( N + 1 ) / 2, Matrix3D( 0.0 ) );
178 : :
179 : 118180 : n = indices1.size();
180 : 118180 : h = 0;
181 [ + + ]: 496370 : for( r = 0; r < n; ++r )
182 : : {
183 [ + - ]: 378190 : const size_t nr = indices1[r];
184 [ + + ]: 1244005 : for( c = r; c < n; ++c )
185 : : {
186 [ + - ]: 865815 : const size_t nc = indices1[c];
187 [ + + ]: 865815 : if( nr <= nc )
188 [ + - ][ + - ]: 612415 : Hessian[N * nr - nr * ( nr + 1 ) / 2 + nc] += Hess1[h++];
[ + - ]
189 : : else
190 [ + - ][ + - ]: 253400 : Hessian[N * nc - nc * ( nc + 1 ) / 2 + nr].plus_transpose_equal( Hess1[h++] );
[ + - ]
191 : : }
192 : : }
193 : :
194 : 118180 : n = indices2.size();
195 : 118180 : h = 0;
196 [ + + ]: 496370 : for( r = 0; r < n; ++r )
197 : : {
198 [ + - ]: 378190 : const size_t nr = indices2[r];
199 [ + + ]: 1244005 : for( c = r; c < n; ++c )
200 : : {
201 [ + - ]: 865815 : const size_t nc = indices2[c];
202 [ + + ]: 865815 : if( nr <= nc )
203 [ + - ][ + - ]: 612415 : Hessian[N * nr - nr * ( nr + 1 ) / 2 + nc] += Hess2[h++];
[ + - ]
204 : : else
205 [ + - ][ + - ]: 253400 : Hessian[N * nc - nc * ( nc + 1 ) / 2 + nr].plus_transpose_equal( Hess2[h++] );
[ + - ]
206 : : }
207 : : }
208 : :
209 : 118180 : value = val1 + val2;
210 [ + - ][ + - ]: 118180 : return rval1 && rval2;
211 : : }
212 : :
213 [ + - ][ + - ]: 4 : } // namespace MBMesquite
|
