1 /* *****************************************************************
2 MESQUITE -- The Mesh Quality Improvement Toolkit
3
4 Copyright 2006 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 retian 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 (2006) kraftche@cae.wisc.edu
24
25 ***************************************************************** */
26
27 #include "Mesquite.hpp"
28 #include <iostream>
29 using std::cout;
30 using std::cerr;
31 using std::endl;
32 using std::ostream;
33 #include <memory>
34 using std::auto_ptr;
35 #include <ctype.h>
36
37 #include "MeshImpl.hpp"
38 #include "MsqError.hpp"
39 #include "InstructionQueue.hpp"
40 #include "TerminationCriterion.hpp"
41 #include "QualityAssessor.hpp"
42 #include "PMeanPTemplate.hpp"
43 #include "PatchPowerMeanP.hpp"
44 #include "ConjugateGradient.hpp"
45 #include "PlanarDomain.hpp"
46 #include "IdealShapeTarget.hpp"
47 #include "ConditionNumberQualityMetric.hpp"
48 #include "ReferenceMesh.hpp"
49 #include "RefMeshTargetCalculator.hpp"
50 #include "TestUtil.hpp"
51
52 #include "TQualityMetric.hpp"
53 #include "ElementPMeanP.hpp"
54 #include "VertexPMeanP.hpp"
55
56 #include "TSquared.hpp"
57 #include "TShapeNB1.hpp"
58 #include "TShapeB1.hpp"
59 #include "TShapeOrientNB1.hpp"
60 #include "TShapeOrientNB2.hpp"
61 #include "TShapeOrientB1.hpp"
62 #include "TShapeSize2DNB1.hpp"
63 #include "TShapeSizeB1.hpp"
64 #include "TShapeSize2DB2.hpp"
65 #include "TShapeSizeB3.hpp"
66 #include "TShapeSizeOrientNB1.hpp"
67 #include "TShapeSizeOrientB1.hpp"
68 #include "TShapeSizeOrientB2.hpp"
69
70 using namespace MBMesquite;
71
72 static const struct { TMetric* u; const char* n; }
73 metrics[] = { { new TSquared, "TSquared" },
74 { new TShapeNB1, "Shape" },
75 { new TShapeB1, "ShapeBarrier" },
76 { new TShapeOrientNB1, "ShapeOrient1" },
77 { new TShapeOrientNB2, "ShapeOrient2" },
78 { new TShapeOrientB1, "ShapeOrientBarrier" },
79 { new TShapeSize2DNB1, "ShapeSize" },
80 { new TShapeSizeB1, "ShapeSizeBarrier1" },
81 { new TShapeSize2DB2, "ShapeSizeBarrier2" },
82 { new TShapeSizeB3, "ShapeSizeBarrier3" },
83 { new TShapeSizeOrientB1, "ShapeSizeOrient1" },
84 { new TShapeSizeOrientB1, "ShapeSizeOrientBarrier1"},
85 { new TShapeSizeOrientB2, "ShapeSizeOrientBarrier2"},
86 { 0, 0 } };
87
88 enum AveragingScheme { NONE = 0, ELEMENT, VERTEX, PATCH };
89 const char* const averaging_names[] = { "none", "element", "vertex", "patch" };
90
91 // default values
92 const double DEFAULT_OF_POWER = 1.0;
93 const unsigned DEFAULT_METRIC_IDX = 0;
94 const AveragingScheme DEFAULT_AVG_SCHEME = NONE;
95 std::string DEFAULT_INPUT_FILE = TestDir + "/2D/vtk/quads/untangled/quads_4by2_bad.vtk";
96 std::string DEFAULT_OUTPUT_FILE = "./out.vtk";
97
98 static PlanarDomain make_domain( Mesh* mesh, MsqError& );
99
do_smoother(const char * input_file,const char * output_file,const char * ref_mesh_file,double of_power,unsigned metric_idx,AveragingScheme avg_scheme)100 static int do_smoother( const char* input_file,
101 const char* output_file,
102 const char* ref_mesh_file,
103 double of_power,
104 unsigned metric_idx,
105 AveragingScheme avg_scheme )
106 {
107 MsqPrintError err(cerr);
108
109 TMetric *const target_metric = metrics[metric_idx].u;
110 cout << "Input file: " << input_file << endl;
111 cout << "Metric: ";
112 if (avg_scheme != NONE)
113 cout << averaging_names[avg_scheme] << " average of ";
114 cout << metrics[metric_idx].n << endl;
115 cout << "Of Power: " << of_power << endl;
116
117
118 auto_ptr<TargetCalculator> tc;
119 auto_ptr<MeshImpl> ref_mesh_impl;
120 auto_ptr<ReferenceMesh> ref_mesh;
121 if (ref_mesh_file) {
122 ref_mesh_impl.reset(new MeshImpl);
123 ref_mesh_impl->read_vtk( ref_mesh_file, err );
124 if (MSQ_CHKERR(err)) return 2;
125 ref_mesh.reset( new ReferenceMesh( ref_mesh_impl.get() ));
126 tc.reset( new RefMeshTargetCalculator( ref_mesh.get() ) );
127 }
128 else {
129 tc.reset( new IdealShapeTarget( ) );
130 }
131
132 TQualityMetric jacobian_metric( tc.get(), target_metric );
133 ElementPMeanP elem_avg( of_power, &jacobian_metric );
134 VertexPMeanP vtx_avg( of_power, &jacobian_metric );
135 QualityMetric* mmetrics[] = { &jacobian_metric, &elem_avg, &vtx_avg, &jacobian_metric };
136 QualityMetric* metric = mmetrics[avg_scheme];
137
138 TerminationCriterion outer, inner;
139 outer.add_iteration_limit( 1 );
140 inner.add_absolute_vertex_movement( 1e-4 );
141 inner.add_iteration_limit( 100 );
142
143 PMeanPTemplate obj1( of_power, metric );
144 PatchPowerMeanP obj2( of_power, metric );
145 ObjectiveFunction& objective = *((avg_scheme == PATCH) ? (ObjectiveFunction*)&obj2 : (ObjectiveFunction*)&obj1);
146
147 ConjugateGradient solver( &objective, err );
148 if (MSQ_CHKERR(err)) return 1;
149 solver.set_inner_termination_criterion( &inner );
150 solver.set_outer_termination_criterion( &outer );
151 solver.use_global_patch();
152
153 ConditionNumberQualityMetric qm_metric;
154 QualityAssessor before_assessor;
155 QualityAssessor after_assessor;
156 before_assessor.add_quality_assessment( metric, 10);
157 before_assessor.add_quality_assessment( &qm_metric );
158 after_assessor.add_quality_assessment( metric, 10 );
159
160 InstructionQueue q;
161 q.add_quality_assessor( &before_assessor, err );
162 q.set_master_quality_improver( &solver, err );
163 q.add_quality_assessor( &after_assessor, err );
164
165 MeshImpl mesh;
166 mesh.read_vtk( input_file, err );
167 if (MSQ_CHKERR(err)) return 2;
168 PlanarDomain geom = make_domain( &mesh, err );
169 if (MSQ_CHKERR(err)) return 1;
170 MeshDomainAssoc mesh_and_domain = MeshDomainAssoc(&mesh, &geom);
171 q.run_instructions( &mesh_and_domain, err );
172 if (MSQ_CHKERR(err)) return 3;
173 mesh.write_vtk( output_file, err );
174 if (MSQ_CHKERR(err)) return 2;
175 cout << "Wrote: " << output_file << endl;
176
177 before_assessor.scale_histograms(&after_assessor);
178
179 return 0;
180 }
181
make_domain(Mesh * mesh,MsqError & err)182 static PlanarDomain make_domain( Mesh* mesh, MsqError& err )
183 {
184 // calculate bounding box of mesh vertices
185 Vector3D minimum( HUGE_VAL, HUGE_VAL, HUGE_VAL );
186 Vector3D maximum( -HUGE_VAL, -HUGE_VAL, -HUGE_VAL );
187 std::vector<Mesh::VertexHandle> vertices;
188 mesh->get_all_vertices( vertices, err );
189 if (MSQ_CHKERR(err)) { return PlanarDomain( minimum, maximum ); }
190 if (vertices.empty()) {
191 std::cerr << "Mesh contains no vertices" << std::endl;
192 exit(1);
193 }
194 std::vector<MsqVertex> coords(vertices.size());
195 mesh->vertices_get_coordinates( arrptr(vertices), arrptr(coords), vertices.size(), err );
196 if (MSQ_CHKERR(err)) { return PlanarDomain( minimum, maximum ); }
197 std::vector<MsqVertex>::const_iterator i;
198 for (i = coords.begin(); i != coords.end(); ++i) {
199 const MsqVertex& v = *i;
200 for (unsigned j = 0; j < 3; ++j) {
201 if (v[j] < minimum[j])
202 minimum[j] = v[j];
203 if (v[j] > maximum[j])
204 maximum[j] = v[j];
205 }
206 }
207
208 // Look for a "zero" plane
209 int k;
210 maximum -= minimum;
211 for (k = 2; k >= 0 && maximum[k] > 1e-6; --k);
212 if (k < 0) {
213 MSQ_SETERR(err)("Cannot determine plane of mesh.", MsqError::INVALID_STATE );
214 return PlanarDomain( minimum, maximum );
215 }
216
217 Vector3D point( 0.0, 0.0, 0.0 ), normal( 0.0, 0.0, 0.0 );
218 normal[k] = 1.0;
219 point[k] = minimum[k];
220 return PlanarDomain( normal, point );
221 }
222
223
usage(const char * argv0,bool help=false)224 static void usage( const char* argv0, bool help = false )
225 {
226 ostream& str = help ? cout : cerr;
227 str << argv0 << " [-p <power>] [-m metric_name] [-a averaging] [-e]"
228 << " -r [ref_mesh] [input_file] [output_file]" << endl
229 << argv0 << " <-l|-h>" << std::endl;
230 if (help) {
231 str << " -p : specify exponent value for p-mean^p OF template (default: " << DEFAULT_OF_POWER << ")" << endl
232 << " -m : specify 2D target metric to use (default: " << metrics[DEFAULT_METRIC_IDX].n << ")" << endl
233 << " -a : specify metric averaging scheme (default: " << averaging_names[DEFAULT_AVG_SCHEME] << ")" << endl
234 << " (none,vertex,element,patch)" << endl
235 << " -e : sample at mid-edge points (default is corners only)" << endl
236 << " -r : use reference mesh instead of ideal elements for targets" << endl
237 << " -l : list available metrics" << endl
238 << " -h : this help output" << endl
239 << " default input file: " << DEFAULT_INPUT_FILE << endl
240 << " default output file: " << DEFAULT_OUTPUT_FILE << endl
241 << endl;
242 }
243 exit (!help);
244 }
245
check_next_arg(int argc,char * argv[],int & i)246 static void check_next_arg( int argc, char* argv[], int& i )
247 {
248 if (i == argc) {
249 cerr << "Expected argument following \"" << argv[i] << '"' << endl;
250 usage( argv[0] );
251 }
252 ++i;
253 }
254
parse_double(int argc,char * argv[],int & i)255 static double parse_double( int argc, char* argv[], int& i )
256 {
257 check_next_arg( argc, argv, i );
258 char* endptr;
259 double result = strtod( argv[i], &endptr );
260 if (endptr == argv[i] || *endptr) {
261 cerr << "Expected real number following \"" << argv[i-1] << '"' << endl;
262 usage( argv[0] );
263 }
264 return result;
265 }
266
comp_string_start(const char * p,const char * s)267 static int comp_string_start( const char* p, const char* s ) {
268 int i;
269 for (i = 0; p[i]; ++i)
270 if (tolower(p[i]) != tolower(s[i]))
271 return 0;
272 return s[i] ? -1 : 1;
273 }
274
parse_averaging(int argc,char * argv[],int & i)275 static AveragingScheme parse_averaging( int argc, char* argv[], int& i )
276 {
277 check_next_arg( argc, argv, i );
278 for (unsigned j = 0; j < 4; ++j)
279 if (comp_string_start( argv[i], averaging_names[j]))
280 return (AveragingScheme)j;
281 cerr << "Expected one of { ";
282 for (unsigned j = 0; j < 4; ++j)
283 cerr << '"' << averaging_names[j] << "\" ";
284 cerr << "} following \"" << argv[i-1] << '"' << endl;
285 usage( argv[0] );
286 return NONE;
287 }
288
parse_metric(int argc,char * argv[],int & i)289 static unsigned parse_metric( int argc, char* argv[], int& i )
290 {
291 check_next_arg( argc, argv, i );
292 unsigned part = 0, all = 0, count = 0, have_all = 0;
293 for (unsigned j = 0; metrics[j].u; ++j) {
294 if (unsigned k = comp_string_start( argv[i], metrics[j].n)) {
295 if (k > 0) {
296 all = j;
297 have_all = 1;
298 }
299 else {
300 part = j;
301 ++count;
302 }
303 }
304 }
305
306 if (have_all)
307 return all;
308
309 if (count) {
310 if (count == 1)
311 return part;
312 cerr << "Ambiguous metric name: \"" << argv[i] << '"' << endl;
313 usage( argv[0] );
314 }
315
316 cerr << "Invalid metric name following \"" << argv[i-1] << "\" option" << endl;
317 usage( argv[0] );
318 return (unsigned)-1;
319 }
320
list_metrics()321 static void list_metrics( )
322 {
323 cout << "Available metrics:" << endl;
324 for (unsigned j = 0; metrics[j].u; ++j)
325 cout << "\t" << metrics[j].n << endl;
326 exit(0);
327 }
328
main(int argc,char * argv[])329 int main( int argc, char* argv[] )
330 {
331 MsqPrintError err(cout);
332
333 // CL settings
334 double of_power = DEFAULT_OF_POWER;
335 unsigned metric_idx = DEFAULT_METRIC_IDX;
336 AveragingScheme avg_scheme = DEFAULT_AVG_SCHEME;
337 const char* input_file = 0;
338 const char* output_file = 0;
339 const char* ref_mesh_file = 0;
340
341 bool proc_opts = true;
342 for (int i = 1; i < argc; ++i) {
343 if (!proc_opts || argv[i][0] != '-') {
344 if (output_file) {
345 cerr << "Unexpected file name: \"" << argv[i] << '"' << endl;
346 usage(argv[0]);
347 }
348 else if (input_file)
349 output_file = argv[i];
350 else
351 input_file = argv[i];
352 continue;
353 }
354
355 if (!argv[i][1] || argv[i][2]) {
356 cerr << "Invalid option: \"" << argv[i] << '"' << endl;
357 usage(argv[0]);
358 }
359
360 switch( argv[i][1] ) {
361 case 'p': of_power = parse_double ( argc, argv, i ); break;
362 case 'm': metric_idx = parse_metric ( argc, argv, i ); break;
363 case 'a': avg_scheme = parse_averaging( argc, argv, i ); break;
364 case 'r': check_next_arg(argc,argv,i); ref_mesh_file=argv[i]; break;
365 case '-': proc_opts = false; break;
366 case 'h': usage( argv[0], true ); break;
367 case 'l': list_metrics(); break;
368 default:
369 cerr << "Invalid option: \"" << argv[i] << '"' << endl;
370 usage(argv[0]);
371 }
372 }
373
374 if (!input_file)
375 input_file = DEFAULT_INPUT_FILE.c_str();
376 if (!output_file)
377 output_file = DEFAULT_OUTPUT_FILE.c_str();
378
379 return do_smoother( input_file,
380 output_file,
381 ref_mesh_file,
382 of_power,
383 metric_idx,
384 avg_scheme );
385 }
386