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
26 ***************************************************************** */
27 // -*- Mode : c++; tab-width: 2; c-tab-always-indent: t; indent-tabs-mode: nil; c-basic-offset: 2 -*-
28 //
29 // SUMMARY:
30 // USAGE:
31 //
32 // ORIG-DATE: 19-Feb-02 at 10:57:52
33 // LAST-MOD: 23-Jul-03 at 18:04:37 by Thomas Leurent
34 //
35 //
36 // DESCRIPTION:
37 // ============
38 /*! \file main.cpp
39
40 describe main.cpp here
41
42 */
43 // DESCRIP-END.
44 //
45
46 #include <iostream>
47 using std::cout;
48 using std::endl;
49 #include <cstdlib>
50
51 #include "Mesquite.hpp"
52 #include "MsqError.hpp"
53 #include "MeshImpl.hpp"
54 #include "Vector3D.hpp"
55 #include "InstructionQueue.hpp"
56 #include "PatchData.hpp"
57 #include "TerminationCriterion.hpp"
58 #include "QualityAssessor.hpp"
59 #include "PlanarDomain.hpp"
60 #include "MsqTimer.hpp"
61 #include "TestUtil.hpp"
62
63 // algorythms
64 #include "ConditionNumberQualityMetric.hpp"
65 #include "LInfTemplate.hpp"
66 #include "SteepestDescent.hpp"
67 #include "LaplacianSmoother.hpp"
68 #include "EdgeLengthQualityMetric.hpp"
69 using namespace MBMesquite;
70
71 std::string DEFAULT_INPUT = TestDir + "/2D/vtk/quads/untangled/square_quad_2.vtk";
72
help(const char * argv0)73 void help(const char* argv0)
74 {
75 std::cerr << "Usage: " << argv0 << " [<input_file>] [<output_file>]" << std::endl
76 << " default input file is: " << DEFAULT_INPUT << std::endl
77 << " defualt is no output file" << std::endl
78 << " Warning: input mesh is assumed to lie in Z=5 plane" << std::endl;
79 exit(1);
80 }
81
main(int argc,char * argv[])82 int main(int argc, char* argv[])
83 {
84 const char* input_file = DEFAULT_INPUT.c_str();
85 const char* output_file = NULL;
86 switch (argc) {
87 default:
88 help(argv[0]);
89 case 3:
90 if (!strcmp(argv[2],"-h"))
91 help(argv[0]);
92 output_file = argv[2];
93 case 2:
94 if (!strcmp(argv[1],"-h"))
95 help(argv[0]);
96 input_file = argv[1];
97 case 1:
98 ;
99 }
100
101 /* Read a VTK Mesh file */
102 MsqPrintError err(cout);
103 MBMesquite::MeshImpl mesh;
104 mesh.read_vtk( input_file, err);
105 if (err) return 1;
106
107 // creates an intruction queue
108 InstructionQueue queue1;
109
110 // creates a mean ratio quality metric ...
111 ConditionNumberQualityMetric shape_metric;
112 EdgeLengthQualityMetric lapl_met;
113 lapl_met.set_averaging_method(QualityMetric::RMS);
114
115 // creates the laplacian smoother procedures
116 LaplacianSmoother lapl1;
117 QualityAssessor stop_qa=QualityAssessor(&shape_metric);
118 stop_qa.add_quality_assessment(&lapl_met);
119
120 //**************Set stopping criterion****************
121 TerminationCriterion sc2;
122 sc2.add_iteration_limit( 10 );
123 if (err) return 1;
124 lapl1.set_outer_termination_criterion(&sc2);
125
126 // adds 1 pass of pass1 to mesh_set1
127 queue1.add_quality_assessor(&stop_qa,err);
128 if (err) return 1;
129 queue1.set_master_quality_improver(&lapl1, err);
130 if (err) return 1;
131 queue1.add_quality_assessor(&stop_qa,err);
132 if (err) return 1;
133 // adds 1 passes of pass2 to mesh_set1
134 // mesh_set1.add_quality_pass(pass2);
135
136 //writeVtkMesh("original_mesh", mesh, err); MSQ_CHKERR(err);
137
138 PlanarDomain plane(Vector3D(0,0,1), Vector3D(0,0,5));
139
140 // launches optimization on mesh_set1
141 MeshDomainAssoc mesh_and_domain = MeshDomainAssoc(&mesh, &plane);
142 Timer t;
143 queue1.run_instructions(&mesh_and_domain, err);
144 if (err) return 1;
145 double secs = t.since_birth();
146 std::cout << "Optimization completed in " << secs << " seconds" << std::endl;
147
148 if (output_file) {
149 mesh.write_vtk(output_file, err);
150 if (err) return 1;
151 std::cout << "Wrote file: " << output_file << std::endl;
152 }
153
154 // check that smoother is working:
155 // the one free vertex must be at the origin
156 if (!DEFAULT_INPUT.compare(input_file)) {
157 std::vector<Mesh::VertexHandle> vertices;
158 mesh.get_all_vertices( vertices, err );
159 if (err) return 1;
160
161 std::vector<bool> fixed_flags;
162 mesh.vertices_get_fixed_flag( arrptr(vertices), fixed_flags, vertices.size(), err );
163 if (err) return 1;
164
165 // find one free vertex
166 int idx = -1;
167 for (unsigned i = 0; i < vertices.size(); ++i) {
168 if (fixed_flags[i] == true)
169 continue;
170 if (idx != -1) {
171 std::cerr << "Multiple free vertices in mesh." << std::endl;
172 return 1;
173 }
174 idx = i;
175 }
176
177 if (idx == -1) {
178 std::cerr << "No free vertex in mesh!!!!!" << std::endl;
179 return 1;
180 }
181
182 Mesh::VertexHandle vertex = vertices[idx];
183 MsqVertex coords;
184 mesh.vertices_get_coordinates( &vertex, &coords, 1, err );
185 if (err) return 1;
186
187 // calculate distance from origin
188 double dist = sqrt( coords[0]*coords[0] + coords[1]*coords[1] );
189 if (dist > 1e-8) {
190 std::cerr << "Free vertex not at origin after Laplace smooth." << std::endl
191 << "Expected location: (0,0)" << std::endl
192 << "Actual location: (" << coords[0] << "," << coords[1] << ")" << std::endl;
193 return 2;
194 }
195 }
196
197 return 0;
198 }
199