Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2010 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 : : (2010) [email protected]
24 : :
25 : : ***************************************************************** */
26 : :
27 : : /** \file DeformingDomainWrapper.cpp
28 : : * \brief Implement DeformingDomainWrapper class
29 : : * \author Jason Kraftcheck
30 : : */
31 : :
32 : : #include "DeformingDomainWrapper.hpp"
33 : : #include "TagVertexMesh.hpp"
34 : : #include "ReferenceMesh.hpp"
35 : : #include "RefMeshTargetCalculator.hpp"
36 : : #include "InstructionQueue.hpp"
37 : : #include "QualityImprover.hpp"
38 : : #include "SteepestDescent.hpp"
39 : : #include "TerminationCriterion.hpp"
40 : : #include "PMeanPTemplate.hpp"
41 : : #include "ElementPMeanP.hpp"
42 : : #include "TQualityMetric.hpp"
43 : : #include "TMixed.hpp"
44 : : #include "TShapeNB1.hpp"
45 : : #include "TShapeSize2DNB1.hpp"
46 : : #include "TShapeSize3DNB1.hpp"
47 : : #include "TShapeSizeOrientNB1.hpp"
48 : : #include "MeshInterface.hpp"
49 : : #include "MsqError.hpp"
50 : : #include "MsqVertex.hpp"
51 : : #include "QualityAssessor.hpp"
52 : : #include "CurveDomain.hpp"
53 : :
54 : : #include <numeric> // for std::accumulate
55 : :
56 : : namespace MBMesquite
57 : : {
58 : :
59 : : const DeformingDomainWrapper::MeshCharacteristic DEFAULT_METRIC_TYPE = DeformingDomainWrapper::SHAPE;
60 : :
61 : : const bool DEFAULT_CULLING = true;
62 : : const double DEFAULT_CPU_TIME = 0.0;
63 : : const double DEFAULT_MOVEMENT_FACTOR = 0.01;
64 : : const int DEFAULT_INNER_ITERATIONS = 2;
65 : : const char DEFAULT_CURVE_TAG[] = "MesquiteCurveFraction";
66 : : const DeformingCurveSmoother::Scheme DEFAULT_CURVE_TYPE = DeformingCurveSmoother::PROPORTIONAL;
67 : :
68 : 3 : DeformingDomainWrapper::DeformingDomainWrapper()
69 : : : metricType( DEFAULT_METRIC_TYPE ), initVertexCoords( "" ), doCulling( DEFAULT_CULLING ),
70 [ + - ]: 3 : cpuTime( DEFAULT_CPU_TIME ), movementFactor( DEFAULT_MOVEMENT_FACTOR )
71 : : {
72 : 3 : }
73 : :
74 [ - + ]: 6 : DeformingDomainWrapper::~DeformingDomainWrapper() {}
75 : :
76 : 3 : void DeformingDomainWrapper::store_initial_mesh( Mesh* mesh, MsqError& err )
77 : : {
78 [ + - ][ + - ]: 3 : TagVertexMesh tool( err, mesh, false, initVertexCoords );MSQ_ERRRTN( err );
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
79 [ + - ][ + - ]: 6 : InstructionQueue q;
80 [ + - ][ + - ]: 3 : q.add_tag_vertex_mesh( &tool, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
81 [ + - ][ + - ]: 3 : q.run_instructions( mesh, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
82 [ + - ][ + - ]: 6 : if( initVertexCoords.empty() ) initVertexCoords = tool.get_tag_name();
[ + - ][ + - ]
83 : : }
84 : :
85 : 1 : void DeformingDomainWrapper::set_vertex_movement_limit_factor( double f )
86 : : {
87 [ - + ]: 1 : assert( f > 0.0 );
88 : 1 : movementFactor = f;
89 : 1 : }
90 : :
91 : 3 : void DeformingDomainWrapper::run_wrapper( MeshDomainAssoc* mesh_and_domain, ParallelMesh* pmesh, Settings* settings,
92 : : QualityAssessor* qa, MsqError& err )
93 : : {
94 [ - + ]: 3 : if( movementFactor <= 0 )
95 : : {
96 : : MSQ_SETERR( err )
97 : : ( MsqError::INVALID_STATE, "Optimization will not terminate with non-positive movement factor %f",
98 [ # # ][ # # ]: 0 : movementFactor );
99 : 3 : return;
100 : : }
101 : :
102 [ + - ]: 3 : Mesh* mesh = mesh_and_domain->get_mesh();
103 [ + - ]: 3 : MeshDomain* geom = mesh_and_domain->get_domain();
104 : :
105 : : // Move initial mesh to domain in case caller did not
106 [ + - ][ + - ]: 3 : move_to_domain( mesh, geom, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
107 : :
108 : 3 : const double P = 1.0;
109 [ + - ][ + - ]: 3 : TagVertexMesh init_mesh( err, mesh );MSQ_ERRRTN( err );
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
110 [ + - ][ + - ]: 6 : ReferenceMesh ref_mesh( &init_mesh );
111 [ + - ][ + - ]: 6 : RefMeshTargetCalculator W( &ref_mesh );
112 : :
113 [ + - ]: 3 : TShapeNB1 mu_s;
114 [ + - ]: 3 : TShapeSize2DNB1 mu_2d_ss;
115 [ + - ][ + - ]: 6 : TShapeSize3DNB1 mu_3d_ss;
116 [ + - ][ + - ]: 6 : TMixed mu_ss( &mu_2d_ss, &mu_3d_ss );
117 [ + - ]: 3 : TShapeSizeOrientNB1 mu_sso;
118 : 3 : TMetric* mu = 0;
119 [ + - - - ]: 3 : switch( metricType )
120 : : {
121 : : case SHAPE:
122 : 3 : mu = &mu_s;
123 : 3 : break;
124 : : case SHAPE_SIZE:
125 : 0 : mu = &mu_ss;
126 : 0 : break;
127 : : case SHAPE_SIZE_ORIENT:
128 : 0 : mu = &mu_sso;
129 : 0 : break;
130 : : }
131 : :
132 [ + - ]: 3 : TQualityMetric sample_metric( &W, mu );
133 [ + - ][ + - ]: 6 : ElementPMeanP elem_metric( P, &sample_metric );
134 [ + - ][ + - ]: 6 : PMeanPTemplate obj_func( P, &elem_metric );
135 : :
136 [ + - ][ + - ]: 6 : TerminationCriterion inner, outer;
[ + - ][ + - ]
[ + - ][ + - ]
137 [ + - ][ + - ]: 6 : SteepestDescent improver( &obj_func );
138 [ + - ]: 3 : improver.use_element_on_vertex_patch(); // Nash optimization
139 [ + - ]: 3 : inner.add_iteration_limit( DEFAULT_INNER_ITERATIONS );
140 [ + - ]: 3 : if( doCulling )
141 [ + - ]: 3 : inner.cull_on_absolute_vertex_movement_edge_length( movementFactor );
142 : : else
143 [ # # ]: 0 : outer.add_absolute_vertex_movement_edge_length( movementFactor );
144 [ - + ][ # # ]: 3 : if( cpuTime > 0.0 ) outer.add_cpu_time( cpuTime );
145 [ + - ]: 3 : improver.set_inner_termination_criterion( &inner );
146 [ + - ]: 3 : improver.set_outer_termination_criterion( &outer );
147 : :
148 [ + - ]: 3 : qa->add_quality_assessment( &elem_metric );
149 [ + - ][ + - ]: 6 : InstructionQueue q;
150 [ + - ][ + - ]: 3 : q.add_quality_assessor( qa, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
151 [ + - ][ + - ]: 3 : q.set_master_quality_improver( &improver, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
152 [ + - ][ + - ]: 3 : q.add_quality_assessor( qa, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
153 [ + - ][ + - ]: 9 : q.run_common( mesh_and_domain, pmesh, settings, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
[ + - ][ + - ]
154 : : }
155 : :
156 : 3 : void DeformingDomainWrapper::move_to_domain( Mesh* mesh, MeshDomain* geom, MsqError& err )
157 : : {
158 [ + - ]: 3 : std::vector< Mesh::VertexHandle > verts;
159 [ + - ][ + - ]: 3 : mesh->get_all_vertices( verts, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
160 : :
161 [ + - ]: 3 : MsqVertex coords;
162 : 3 : std::vector< Mesh::VertexHandle >::const_iterator i;
163 [ + - ][ + - ]: 744 : for( i = verts.begin(); i != verts.end(); ++i )
[ + - ][ + + ]
[ + - ]
164 : : {
165 [ + - ][ + - ]: 741 : mesh->vertices_get_coordinates( &*i, &coords, 1, err );MSQ_ERRRTN( err );
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
166 [ + - ][ + - ]: 741 : geom->snap_to( *i, coords );
167 [ + - ][ + - ]: 741 : mesh->vertex_set_coordinates( *i, coords, err );MSQ_ERRRTN( err );
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
168 : 3 : }
169 : : }
170 : :
171 [ + - ]: 6 : DeformingCurveSmoother::DeformingCurveSmoother() : metricType( DEFAULT_CURVE_TYPE ), initFractTag( DEFAULT_CURVE_TAG )
172 : : {
173 : 3 : }
174 : :
175 : 6 : DeformingCurveSmoother::~DeformingCurveSmoother() {}
176 : :
177 : 12 : void DeformingCurveSmoother::store_initial_mesh( Mesh* mesh, const Mesh::VertexHandle* verts, int nverts, CurveDomain*,
178 : : MsqError& err )
179 : : {
180 [ - + ]: 12 : if( nverts < 2 )
181 : : {
182 : : MSQ_SETERR( err )
183 [ # # ][ # # ]: 0 : ( "Invalid curve mesh. Need at least two end vertices.", MsqError::INVALID_MESH );
184 : 12 : return;
185 : : }
186 : :
187 : : // get edge lengths
188 [ + - ]: 12 : std::vector< double > vals( nverts - 1 );
189 [ + - ][ + + ]: 36 : MsqVertex coords[2];
190 : 12 : int prev = 0;
191 [ + - ][ + - ]: 12 : mesh->vertices_get_coordinates( verts, coords + prev, 1, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
192 [ + + ]: 192 : for( int i = 1; i < nverts; ++i )
193 : : {
194 : 180 : int next = 1 - prev;
195 [ + - ][ + - ]: 180 : mesh->vertices_get_coordinates( verts + i, coords + next, 1, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
196 [ + - ][ + - ]: 180 : vals[i - 1] = ( coords[0] - coords[1] ).length();
[ + - ]
197 : 180 : prev = next;
198 : : }
199 : :
200 : : // convert to length fraction before each iterior vertex
201 : : // (sum of lengths of adjacent edges over total curve length)
202 [ + - ]: 12 : const double total = std::accumulate( vals.begin(), vals.end(), 0.0 );
203 [ + + ]: 180 : for( int i = 1; i < nverts - 1; ++i )
204 [ + - ][ + - ]: 168 : vals[i - 1] = vals[i - 1] / total;
205 [ + - ]: 12 : vals.resize( nverts - 2 );
206 : :
207 : : // create tag
208 [ + - ]: 12 : TagHandle tag = mesh->tag_create( initFractTag, Mesh::DOUBLE, 1, 0, err );
209 [ + - ][ + + ]: 12 : if( err.error_code() == MsqError::TAG_ALREADY_EXISTS )
210 : : {
211 [ + - ]: 9 : err.clear();
212 [ + - ]: 9 : tag = get_tag( mesh, err );
213 : : }
214 [ + - ][ - + ]: 12 : MSQ_ERRRTN( err );
[ # # ][ # # ]
[ - + ]
215 : :
216 : : // store tag data on interior vertices
217 [ + - ][ + - ]: 12 : mesh->tag_set_vertex_data( tag, nverts - 2, verts + 1, &vals[0], err );MSQ_ERRRTN( err );
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ][ + - ]
218 : : }
219 : :
220 : 12 : void DeformingCurveSmoother::smooth_curve( Mesh* mesh, const Mesh::VertexHandle* verts, int nverts, CurveDomain* geom,
221 : : Scheme, MsqError& err )
222 : : {
223 [ - + ]: 12 : if( nverts < 2 )
224 : : {
225 : : MSQ_SETERR( err )
226 [ # # ][ # # ]: 0 : ( "Invalid curve mesh. Need at least two end vertices.", MsqError::INVALID_MESH );
227 : 12 : return;
228 : : }
229 : :
230 : : // Verify that end vertices are on curve. We cannot snap end vertices
231 : : // to ends of curve for application because we don't know where the
232 : : // ends of the curve are.
233 [ + - ][ + + ]: 36 : MsqVertex coords[2], coords2;
[ + - ]
234 : 12 : Mesh::VertexHandle ends[2] = { verts[0], verts[nverts - 1] };
235 [ + + ]: 36 : for( int i = 0; i < 2; ++i )
236 : : {
237 [ + - ][ + - ]: 24 : mesh->vertices_get_coordinates( ends + i, coords + i, 1, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
238 [ + - ][ + - ]: 24 : geom->position_from_length( coords[i].to_array(), 0, coords2.to_array(), err );MSQ_ERRRTN( err );
[ + - ][ + - ]
[ - + ][ # # ]
[ # # ][ - + ]
239 [ + - ][ + - ]: 24 : if( ( coords[i] - coords2 ).length_squared() > DBL_EPSILON )
[ - + ]
240 : : {
241 : : MSQ_SETERR( err )
242 [ # # ][ # # ]: 0 : ( "Curve end vertices do not line on curve. Move ends to curve end points first", MsqError::INVALID_MESH );
243 : 0 : return;
244 : : }
245 : : }
246 [ + - ][ + - ]: 12 : const double total = geom->arc_length( coords[0].to_array(), coords[1].to_array(), err );MSQ_ERRRTN( err );
[ + - ][ + - ]
[ - + ][ # # ]
[ # # ][ - + ]
247 : :
248 [ + - ]: 12 : std::vector< double > vals( nverts - 1 );
249 [ - + ]: 12 : if( metricType == EQUAL )
250 : : {
251 [ # # ]: 0 : std::fill( vals.begin(), vals.end(), 1.0 / ( nverts - 1 ) );
252 : : // fracsum = 1.0;
253 : : }
254 : : else
255 : : { // metricType == PROPORTIONAL
256 [ + - ][ + - ]: 12 : TagHandle tag = get_tag( mesh, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
257 [ + - ][ + - ]: 12 : mesh->tag_get_vertex_data( tag, nverts - 2, verts + 1, &vals[0], err );MSQ_ERRRTN( err );
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
258 [ + - ][ + - ]: 12 : double sum = std::accumulate( vals.begin(), vals.end() - 1, 0.0 );
259 [ + - ]: 12 : if( 1.0 - sum > 1e-8 )
260 [ + - ]: 12 : vals.back() = 1.0 - sum;
261 : : else
262 : : {
263 [ # # ][ # # ]: 0 : vals.back() = *std::min_element( vals.begin(), vals.end() - 1 );
[ # # ][ # # ]
264 [ # # ]: 0 : sum += vals.back();
265 [ # # ]: 0 : for( size_t i = 0; i < vals.size(); ++i )
266 [ # # ]: 0 : vals[i] /= sum;
267 : : }
268 : : }
269 : :
270 : 12 : double frac_sum = 0.0;
271 [ + + ][ + - ]: 180 : for( int i = 1; i < nverts - 1; ++i )
272 : : {
273 [ + - ]: 168 : frac_sum += vals[i - 1];
274 [ + - ][ + - ]: 168 : geom->position_from_length( coords[0].to_array(), total * frac_sum, coords[1].to_array(), err );MSQ_ERRRTN( err );
[ + - ][ + - ]
[ - + ][ # # ]
[ # # ][ - + ]
275 [ + - ][ + - ]: 168 : mesh->vertex_set_coordinates( verts[i], coords[1], err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
276 : 12 : }
277 : : }
278 : :
279 : 21 : TagHandle DeformingCurveSmoother::get_tag( Mesh* mesh, MsqError& err )
280 : : {
281 [ + - ]: 21 : TagHandle h = mesh->tag_get( initFractTag, err );
282 [ + - ][ - + ]: 21 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
283 [ + - ]: 21 : std::string name;
284 : : Mesh::TagType type;
285 : : unsigned len;
286 [ + - ]: 21 : mesh->tag_properties( h, name, type, len, err );
287 [ + - ][ - + ]: 21 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
288 [ + - ][ - + ]: 21 : if( type != Mesh::DOUBLE || len != 1 )
289 : : {
290 : : MSQ_SETERR( err )
291 [ # # ][ # # ]: 0 : ( MsqError::INVALID_MESH, "Tag \"%s\" exists but is not 1 double value per vertex.", initFractTag.c_str() );
292 : : }
293 : 21 : return h;
294 : : }
295 : :
296 [ + - ][ + - ]: 8 : } // namespace MBMesquite
|
