1 /* *****************************************************************
2 MESQUITE -- The Mesh Quality Improvement Toolkit
3
4 Copyright 2004 Sandia Corporation and Argonne National
5 Laboratory. Under the terms of Contract DE-AC04-94AL85000
6 with Sandia Corporation, the U.S. Government 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 diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov,
24 pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov
25 kraftche@cae.wisc.edu
26
27 ***************************************************************** */
28
29 #define TOL 1e-5
30
31 #include <iostream>
32 using std::cout;
33 using std::endl;
34 #include <stdlib.h>
35
36 #include "Mesquite.hpp"
37 #include "MeshImpl.hpp"
38 #include "MsqError.hpp"
39 #include "InstructionQueue.hpp"
40 #include "TerminationCriterion.hpp"
41 #include "QualityAssessor.hpp"
42 #include "LPtoPTemplate.hpp"
43 #include "LInfTemplate.hpp"
44
45 // algorithms
46 #include "IdealWeightMeanRatio.hpp"
47 #include "IdealWeightInverseMeanRatio.hpp"
48 #include "UntangleBetaQualityMetric.hpp"
49 #include "FeasibleNewton.hpp"
50 #include "ConjugateGradient.hpp"
51 #include "TrustRegion.hpp"
52 #include "ConditionNumberQualityMetric.hpp"
53 #include "TestUtil.hpp"
54
55 using namespace MBMesquite;
56
57 // Use CPPUNIT_ASSERT in code so it's easy to convert to a unit test later.
58 #define CPPUNIT_ASSERT(A) \
59 do { if (!(A)) { \
60 std::cout << "Assertion Failed: " << #A << std::endl; \
61 std::cout << " File: " << __FILE__ << std::endl; \
62 std::cout << " Line: " << __LINE__ << std::endl; \
63 return true; \
64 } } while (false)
65
66
67 // Given a mesh with a single free vertex located at the origin,
68 // move the vertex to the specified position, smooth the mesh,
69 // and verify that the vertex was moved back to the origin by
70 // the smoother.
71 bool smooth_mesh( Mesh* mesh, Mesh* ref_mesh,
72 Mesh::VertexHandle free_vertex_at_origin,
73 Vector3D initial_free_vertex_position,
74 QualityMetric* metric );
75
76 bool smooth_mixed_mesh( const char* filename );
77
main(int argc,char * argv[])78 int main( int argc, char* argv[] )
79 {
80 unsigned i;
81 std::string input_file = TestDir + "/3D/vtk/mixed/tangled/mixed-hex-pyr-tet.vtk";
82 if (argc == 2)
83 input_file = argv[1];
84 else if (argc != 1)
85 {
86 std::cerr << "Invalid arguments.\n";
87 return 2;
88 }
89
90
91 MBMesquite::MsqPrintError err(cout);
92 IdealWeightMeanRatio m1;
93 IdealWeightInverseMeanRatio m2(err);
94 ConditionNumberQualityMetric m3;
95 QualityMetric* metrics[] = { &m1, &m2, &m3, 0 };
96
97 // Read Mesh
98 std::string mesh_file = TestDir + "/3D/vtk/pyramids/untangled/12-pyramid-unit-sphere.vtk";
99 std::string imesh_file = TestDir + "/3D/vtk/pyramids/untangled/12-pyramid-unit-sphere.vtk";
100 MBMesquite::MeshImpl mesh;
101 mesh.read_vtk(mesh_file.c_str(), err);
102 CPPUNIT_ASSERT(!err);
103 MBMesquite::MeshImpl ideal_mesh;
104 ideal_mesh.read_vtk(imesh_file.c_str(), err);
105 CPPUNIT_ASSERT(!err);
106
107 // Check that the mesh read correctly, and contains what is
108 // expected later.
109
110 // Get mesh data
111 // Expecting file to contain 12 pyramid elements constructed
112 // from 15 vertices.
113 std::vector<Mesh::VertexHandle> vert_array;
114 std::vector<Mesh::ElementHandle> elem_array;
115 std::vector<size_t> conn_offsets;
116 mesh.get_all_elements( elem_array, err );
117 CPPUNIT_ASSERT(!err);
118 CPPUNIT_ASSERT( elem_array.size() == 12 );
119 mesh.elements_get_attached_vertices( arrptr(elem_array),
120 elem_array.size(),
121 vert_array,
122 conn_offsets,
123 err );
124 CPPUNIT_ASSERT(!err);
125 CPPUNIT_ASSERT(vert_array.size() == 60);
126 CPPUNIT_ASSERT(conn_offsets.size() == 13);
127 EntityTopology type_array[12];
128 mesh.elements_get_topologies( arrptr(elem_array), type_array, 12, err );
129 CPPUNIT_ASSERT(!err);
130
131 // Verify element types and number of vertices
132 for (i = 0; i < 12; ++i)
133 {
134 CPPUNIT_ASSERT( type_array[i] == PYRAMID );
135 CPPUNIT_ASSERT( conn_offsets[i] == 5*i );
136 }
137
138 // All pyramids should share a common apex, at the
139 // center of the sphere
140 Mesh::VertexHandle apex_handle = vert_array[4];
141 for (i = 1; i < 12; ++i)
142 {
143 CPPUNIT_ASSERT( vert_array[5*i+4] == apex_handle );
144 }
145
146 // Verify that apex is at origin and all other vertices are
147 // on unit sphere
148 MsqVertex vertices[60];
149 mesh.vertices_get_coordinates( arrptr(vert_array), vertices, 60, err );
150 CPPUNIT_ASSERT(!err);
151 for (i = 0; i < 60; ++i)
152 {
153 if (vert_array[i] == apex_handle)
154 CPPUNIT_ASSERT( vertices[i].within_tolerance_box( Vector3D(0,0,0), 1e-6 ) );
155 else
156 CPPUNIT_ASSERT( fabs(1.0 - vertices[i].length()) < 1e-6 );
157 }
158
159 // Try smoothing w/out moving the free vertex and verify that
160 // the smoother didn't move the vertex
161 Vector3D position(0,0,0);
162 for (i = 0; metrics[i] != NULL; ++i)
163 CPPUNIT_ASSERT( !smooth_mesh( &mesh, &ideal_mesh, apex_handle, position, metrics[i] ) );
164
165 // Now try moving the vertex and see if the smoother moves it back
166 // to the origin
167 position.set( 0.1, 0.1, 0.1 );
168 for (i = 0; metrics[i] != NULL; ++i)
169 CPPUNIT_ASSERT( !smooth_mesh( &mesh, &ideal_mesh, apex_handle, position, metrics[i] ) );
170
171 // Now try moving the vertex further and see if the smoother moves it back
172 // to the origin
173 position.set( 0.3, 0.3, 0.3 );
174 for (i = 0; metrics[i] != NULL; ++i)
175 CPPUNIT_ASSERT( !smooth_mesh( &mesh, &ideal_mesh, apex_handle, position, metrics[i] ) );
176
177 // Now try smoothing a real mixed mesh
178 CPPUNIT_ASSERT( !smooth_mixed_mesh( input_file.c_str() ) );
179
180 return 0;
181 }
182
183
smooth_mesh(Mesh * mesh,Mesh *,Mesh::VertexHandle free_vertex_at_origin,Vector3D initial_free_vertex_position,QualityMetric * metric)184 bool smooth_mesh( Mesh* mesh, Mesh* ,
185 Mesh::VertexHandle free_vertex_at_origin,
186 Vector3D initial_free_vertex_position,
187 QualityMetric* metric )
188 {
189 MBMesquite::MsqPrintError err(cout);
190 const Vector3D origin( 0, 0, 0 );
191
192 // print a little output so we know when we died
193 std::cout <<
194 "**************************************************************************"
195 << std::endl <<
196 "* Smoothing..."
197 << std::endl <<
198 "* Metric: " << metric->get_name()
199 << std::endl <<
200 "* Apex position: " << initial_free_vertex_position
201 << std::endl //<<
202 //"**************************************************************************"
203 << std::endl;
204
205
206 // Set free vertex to specified position
207 mesh->vertex_set_coordinates( free_vertex_at_origin,
208 initial_free_vertex_position,
209 err );
210 CPPUNIT_ASSERT(!err);
211
212 // Create an InstructionQueue
213 InstructionQueue Q;
214
215 // Set up objective function
216 LPtoPTemplate obj_func(metric, 1, err);
217 CPPUNIT_ASSERT(!err);
218
219 // Create solver
220 FeasibleNewton solver( &obj_func );
221 CPPUNIT_ASSERT(!err);
222 solver.use_global_patch();
223 CPPUNIT_ASSERT(!err);
224
225 // Set stoping criteria for solver
226 TerminationCriterion tc_inner;
227 tc_inner.add_absolute_vertex_movement( 1e-6 );
228 solver.set_inner_termination_criterion(&tc_inner);
229
230 TerminationCriterion tc_outer;
231 tc_outer.add_iteration_limit( 1 );
232 solver.set_outer_termination_criterion(&tc_outer);
233
234 // Add solver to queue
235 Q.set_master_quality_improver(&solver, err);
236 CPPUNIT_ASSERT(!err);
237
238 // And smooth...
239 Q.run_instructions(mesh, err);
240 CPPUNIT_ASSERT(!err);
241
242 // Verify that vertex was moved back to origin
243 MsqVertex vtx;
244 mesh->vertices_get_coordinates( &free_vertex_at_origin, &vtx, 1, err );
245 CPPUNIT_ASSERT( !err );
246 Vector3D position = vtx;
247
248 // print a little output so we know when we died
249 std::cout //<<
250 //"**************************************************************************"
251 << std::endl <<
252 "* Done Smoothing:"
253 << std::endl <<
254 "* Metric: " << metric->get_name()
255 << std::endl <<
256 "* Apex position: " << position
257 << std::endl <<
258 "**************************************************************************"
259 << std::endl;
260
261 CPPUNIT_ASSERT( position.within_tolerance_box( Vector3D(0,0,0), TOL ) );
262 return false;
263 }
264
265
266
267
smooth_mixed_mesh(const char * filename)268 bool smooth_mixed_mesh( const char* filename )
269 {
270 MBMesquite::MsqPrintError err(cout);
271
272 // print a little output so we know when we died
273 std::cout <<
274 "**************************************************************************"
275 << std::endl <<
276 "* Smoothing: " << filename
277 << std::endl <<
278 "**************************************************************************"
279 << std::endl;
280
281 // The instruction queue to set up
282 InstructionQueue Q;
283
284 // Use numeric approx of derivitives until analytic solutions
285 // are working for pyramids
286 IdealWeightInverseMeanRatio mr_metric(err);
287 //sRI_DFT dft_metric;
288 UntangleBetaQualityMetric un_metric(0);
289 CPPUNIT_ASSERT(!err);
290
291 // Create Mesh object
292 MBMesquite::MeshImpl mesh;
293 mesh.read_vtk(filename, err);
294 CPPUNIT_ASSERT(!err);
295
296 // Set up a preconditioner
297 LInfTemplate pre_obj_func( &un_metric );
298 ConjugateGradient precond( &pre_obj_func, err ); CPPUNIT_ASSERT(!err);
299 precond.use_element_on_vertex_patch();
300 TerminationCriterion pre_term, pre_outer;
301 //pre_term.add_relative_quality_improvement( 0.1 );
302 pre_term .add_iteration_limit( 3 );
303 pre_outer.add_iteration_limit( 1 );
304 CPPUNIT_ASSERT(!err);
305 precond.set_inner_termination_criterion( &pre_term );
306 precond.set_outer_termination_criterion( &pre_outer );
307 //precond.use_element_on_vertex_patch();
308
309 // Set up objective function
310 LPtoPTemplate obj_func(&mr_metric, 1, err);
311 CPPUNIT_ASSERT(!err);
312
313 // Create solver
314 FeasibleNewton solver( &obj_func );
315 CPPUNIT_ASSERT(!err);
316 solver.use_global_patch();
317 CPPUNIT_ASSERT(!err);
318
319 // Set stoping criteria for solver
320 TerminationCriterion tc_inner;
321 tc_inner.add_relative_quality_improvement( 0.25 );
322 solver.set_inner_termination_criterion(&tc_inner);
323
324 TerminationCriterion tc_outer;
325 tc_outer.add_iteration_limit( 1 );
326 CPPUNIT_ASSERT(!err);
327 solver.set_outer_termination_criterion(&tc_outer);
328
329 // Create a QualityAssessor
330 MBMesquite::QualityAssessor qa;
331 qa.add_quality_assessment( &mr_metric );
332 qa.add_quality_assessment( &un_metric );
333 Q.add_quality_assessor( &qa, err );
334 CPPUNIT_ASSERT(!err);
335
336 // Add untangler to queue
337 Q.add_preconditioner( &precond, err ); CPPUNIT_ASSERT(!err);
338 Q.add_quality_assessor( &qa, err );
339 CPPUNIT_ASSERT(!err);
340
341 // Add solver to queue
342 Q.set_master_quality_improver(&solver, err);
343 CPPUNIT_ASSERT(!err);
344 Q.add_quality_assessor( &qa, err );
345 CPPUNIT_ASSERT(!err);
346
347 // And smooth...
348 Q.run_instructions(&mesh, err);
349 CPPUNIT_ASSERT(!err);
350
351 return false;
352 }
353
354