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) kraftche@cae.wisc.edu
24
25 ***************************************************************** */
26
27
28 /** \file main.cpp
29 * \brief
30 * \author Jason Kraftcheck
31 */
32
33 #include "TestUtil.hpp"
34 #include "MeshImpl.hpp"
35 #include "PlanarDomain.hpp"
36 #include "TriLagrangeShape.hpp"
37 #include "QuadLagrangeShape.hpp"
38 #include "IdealShapeTarget.hpp"
39 #include "TShapeNB1.hpp"
40 #include "TQualityMetric.hpp"
41 #include "SteepestDescent.hpp"
42 #include "SlaveBoundaryVertices.hpp"
43 #include "PMeanPTemplate.hpp"
44 #include "InstructionQueue.hpp"
45 #include "MsqVertex.hpp"
46 #include "MsqError.hpp"
47 #include "IdealShapeTarget.hpp"
48 #include "QualityAssessor.hpp"
49 #include "PatchData.hpp"
50
51 #include <iostream>
52 #include <algorithm>
53 #include <assert.h>
54
55 using namespace MBMesquite;
56
57 std::string DEFAULT_INPUT_FILE = std::string ( STRINGIFY(SRCDIR) ) + "/input.vtk";
58
usage(const char * argv0)59 void usage( const char* argv0 )
60 {
61 std::cout << "Usage: " << argv0 << " [<input_file>] [-o <output_file_base>]" << std::endl;
62 exit(1);
63 }
64
65 int check_slaved_coords( Mesh& mesh, Settings::HigherOrderSlaveMode mode,
66 std::string name, bool have_slaved_flag,
67 MsqError& err );
68
69 // compare vertex coodinates between two topologically equivalent meshes
70 int compare_node_coords( Mesh& mesh1, Mesh& mesh2, MsqError& err );
71
72 // verify that no element corner node is marked as slaved
73 int check_no_slaved_corners( Mesh& mesh, MsqError& err );
74
75 // compare slaved flags in mesh to slaved state in global patch data
76 int check_global_patch_slaved( Mesh& mesh, MsqError& err );
77
78 // tag vertices marked as slaved in the PatchData
79 void tag_patch_slaved( Mesh& mesh,
80 Settings::HigherOrderSlaveMode mode,
81 MsqError& err );
82
main(int argc,char * argv[])83 int main( int argc, char* argv[] )
84 {
85 const char* input_file = DEFAULT_INPUT_FILE.c_str();
86 const char* output_file_base = 0;
87 bool expect_output_base = false;
88 for (int i = 1; i < argc; ++i) {
89 if (expect_output_base) {
90 output_file_base = argv[i];
91 expect_output_base = false;
92 }
93 else if (!strcmp(argv[i],"-o"))
94 expect_output_base = true;
95 else if (DEFAULT_INPUT_FILE.compare(input_file))
96 usage(argv[0]);
97 else
98 input_file = argv[i];
99 }
100 if (expect_output_base)
101 usage(argv[0]);
102
103 MsqPrintError err(std::cerr);
104 SlaveBoundaryVertices slaver(1);
105 TShapeNB1 tmetric;
106 IdealShapeTarget target;
107 TQualityMetric metric( &target, &tmetric );
108 PMeanPTemplate of( 1.0, &metric );
109 SteepestDescent improver( &of );
110 TerminationCriterion inner;
111 inner.add_absolute_vertex_movement( 1e-3 );
112 improver.set_inner_termination_criterion( &inner );
113 QualityAssessor assess( &metric );
114 InstructionQueue q;
115 q.set_master_quality_improver( &improver, err );
116 q.add_quality_assessor( &assess, err );
117
118 TriLagrangeShape trishape;
119 QuadLagrangeShape quadshape;
120 q.set_mapping_function( &trishape );
121 q.set_mapping_function( &quadshape );
122
123 const int NUM_MODES = 4;
124
125 Settings::HigherOrderSlaveMode modes[NUM_MODES] =
126 { Settings::SLAVE_NONE,
127 Settings::SLAVE_ALL,
128 Settings::SLAVE_CALCULATED,
129 Settings::SLAVE_FLAG
130 };
131
132 std::string names[NUM_MODES] =
133 { "NONE",
134 "ALL",
135 "CALCULATED",
136 "FLAG" };
137
138 MeshImpl meshes[NUM_MODES];
139 std::vector<MeshDomainAssoc> meshes_and_domains;
140
141 bool have_slaved_flag = true;
142 std::vector<bool> flag(1);
143 for (int i = 0; i < NUM_MODES; ++i) {
144 std::cout << std::endl
145 << "-----------------------------------------------" << std::endl
146 << " Mode: " << names[i] << std::endl
147 << "-----------------------------------------------" << std::endl;
148
149 meshes[i].read_vtk( input_file, err );
150 if (err) return 1;
151
152 if (modes[i] == Settings::SLAVE_CALCULATED) {
153 q.add_vertex_slaver( &slaver, err );
154 }
155 else if (modes[i] == Settings::SLAVE_FLAG) {
156 std::vector<Mesh::VertexHandle> verts;
157 meshes[i].get_all_vertices( verts, err );
158 if (err) return 1;
159 meshes[i].vertices_get_slaved_flag( arrptr(verts), flag, 1, err );
160 if (err) {
161 have_slaved_flag = false;
162 std::cout << "Skipped because input file does not contain slaved attribute" << std::endl;
163 err.clear();
164 continue;
165 }
166 }
167
168
169 if (have_slaved_flag && modes[i] == Settings::SLAVE_FLAG) {
170 if (!check_no_slaved_corners( meshes[i], err ) || err)
171 return 1;
172 if (!check_global_patch_slaved( meshes[i], err ) || err)
173 return 1;
174 }
175
176
177 PlanarDomain plane;
178 plane.fit_vertices( &meshes[i], err );
179 if (err) return 1;
180
181 q.set_slaved_ho_node_mode( modes[i] );
182 meshes_and_domains.push_back(MeshDomainAssoc(&meshes[i], &plane));
183 q.run_instructions( &meshes_and_domains[i], err );
184 if (err) return 1;
185
186 if (modes[i] == Settings::SLAVE_CALCULATED) {
187 q.remove_vertex_slaver( &slaver, err );
188 }
189
190 if (output_file_base) {
191 tag_patch_slaved( meshes[i], modes[i], err );
192 std::string name(output_file_base);
193 name += "-";
194 name += names[i];
195 name += ".vtk";
196 meshes[i].write_vtk( name.c_str(), err );
197 if (err) return 1;
198 }
199 }
200
201 int exit_code = 0;
202 if (!DEFAULT_INPUT_FILE.compare(input_file)) {
203 for (int i = 0; i < NUM_MODES; ++i) {
204 std::cout << std::endl
205 << "-----------------------------------------------" << std::endl
206 << " Mode: " << names[i] << std::endl
207 << "-----------------------------------------------" << std::endl;
208
209 exit_code += check_slaved_coords( meshes[i], modes[i], names[i], have_slaved_flag, err );
210 if (err) return 1;
211 }
212
213 // flags should correspond to same slaved nodes as calculated,
214 // so resulting meshes should be identical.
215 if (have_slaved_flag) {
216 int flag_idx = std::find( modes, modes+NUM_MODES, Settings::SLAVE_FLAG ) - modes;
217 int calc_idx = std::find( modes, modes+NUM_MODES, Settings::SLAVE_CALCULATED ) - modes;
218 exit_code += compare_node_coords( meshes[flag_idx], meshes[calc_idx], err );
219 if (err) return 1;
220 }
221 }
222
223 return exit_code;
224 }
225
get_slaved_coords(Mesh & mesh,Mesh::VertexHandle vertex,MsqError & err)226 Vector3D get_slaved_coords( Mesh& mesh, Mesh::VertexHandle vertex, MsqError& err )
227 {
228 std::vector<Mesh::ElementHandle> elem;
229 std::vector<size_t> off;
230 mesh.vertices_get_attached_elements( &vertex, 1, elem, off, err );
231 if (MSQ_CHKERR(err)) return Vector3D(0.0);
232
233 std::vector<Mesh::VertexHandle> verts;
234 mesh.elements_get_attached_vertices( arrptr(elem), 1, verts, off, err );
235 if (MSQ_CHKERR(err)) return Vector3D(0.0);
236
237 EntityTopology type;
238 mesh.elements_get_topologies( arrptr(elem), &type, 1, err );
239 if (MSQ_CHKERR(err)) return Vector3D(0.0);
240
241 size_t idx = std::find( verts.begin(), verts.end(), vertex ) - verts.begin();
242 unsigned side_dim, side_num;
243 TopologyInfo::side_from_higher_order( type, verts.size(), idx, side_dim, side_num, err);
244 if (MSQ_CHKERR(err)) return Vector3D(0.0);
245
246 // just return the mean of the corner vertices defining the side.
247 // this should always be correct for mid-edge nodes but is a bit
248 // dubious for mid-face or mid-region nodes. But our default input
249 // file (the only file for which we should be in this function) will
250 // have only mid-edge nodes.
251 unsigned n;
252 const unsigned* side_vtx = TopologyInfo::side_vertices( type, side_dim, side_num, n, err );
253 if (MSQ_CHKERR(err))
254 return Vector3D(0.0);
255
256 Vector3D sum(0.0);
257 for (unsigned i = 0; i < n; ++i) {
258 MsqVertex coords;
259 Mesh::VertexHandle vtx = verts[side_vtx[i]];
260 mesh.vertices_get_coordinates( &vtx, &coords, 1, err );
261 if (MSQ_CHKERR(err)) return Vector3D(0.0);
262 sum += coords;
263 }
264 return sum / n;
265 }
266
check_slaved_coords(Mesh & mesh,Settings::HigherOrderSlaveMode mode,std::string name,bool have_slaved_flag,MsqError & err)267 int check_slaved_coords( Mesh& mesh,
268 Settings::HigherOrderSlaveMode mode,
269 std::string name,
270 bool have_slaved_flag,
271 MsqError& err )
272 {
273 const double EPSILON = 1e-4;
274
275 // We can distinguish between nodes near the boundary and those
276 // further inside based on the slaved attribute flag in the
277 // the default input file. The default input file should always
278 // have it defined
279 if (!have_slaved_flag) {
280 std::cerr << "slaved flag not specified in input file. Cannot validate results." << std::endl;
281 return 1;
282 }
283
284 // Get list of all higher-order nodes
285 std::vector<Mesh::ElementHandle> elems;
286 std::vector<Mesh::VertexHandle> verts;
287 std::vector<size_t> offsets;
288 mesh.get_all_elements( elems, err ); MSQ_ERRZERO(err);
289 std::vector<EntityTopology> types(elems.size());
290 mesh.elements_get_topologies( arrptr(elems), arrptr(types), elems.size(), err );
291 MSQ_ERRZERO(err);
292 mesh.elements_get_attached_vertices( arrptr(elems), elems.size(),
293 verts, offsets, err );
294 MSQ_ERRZERO(err);
295 std::vector<Mesh::VertexHandle>::iterator r, e, w = verts.begin();
296 for (size_t i = 0; i < elems.size(); ++i) {
297 r = verts.begin() + offsets[i] + TopologyInfo::corners( types[i] );
298 e = verts.begin() + offsets[i+1];
299 w = std::copy( r, e, w );
300 }
301 std::sort( verts.begin(), w );
302 verts.erase( std::unique( verts.begin(), w ), verts.end() );
303
304 // Get lists of slaved and non-slaved free vertices, where 'slaved'
305 // are those vertices far from the boundary that would be slaved if
306 // mode == SLAVE_FLAG, not those that were actually slaved during
307 // the optimization
308 std::vector<bool> fixed, slaved;
309 mesh.vertices_get_fixed_flag( arrptr(verts), fixed, verts.size(), err );
310 if (MSQ_CHKERR(err)) {
311 return 1;
312 }
313 mesh.vertices_get_slaved_flag( arrptr(verts), slaved, verts.size(), err );
314 if (MSQ_CHKERR(err)) {
315 return 1;
316 }
317 std::vector<Mesh::VertexHandle> free, slave;
318 for (size_t i = 0; i < verts.size(); ++i) {
319 if (!fixed[i] && slaved[i])
320 slave.push_back( verts[i] );
321 else if (!fixed[i] && !slaved[i])
322 free.push_back( verts[i] );
323 }
324
325 // get all coordinates
326 std::vector<MsqVertex> free_coords(free.size()), slave_coords(slave.size());
327 mesh.vertices_get_coordinates( arrptr(free), arrptr(free_coords), free.size(), err );
328 MSQ_ERRZERO(err);
329 mesh.vertices_get_coordinates( arrptr(slave), arrptr(slave_coords), slave.size(), err );
330 MSQ_ERRZERO(err);
331
332 int error_count = 0;
333
334 // Expect vertices near boundary to be at slave vertex positions
335 // if mode was SLAVE_ALL
336 if (mode == Settings::SLAVE_ALL) {
337 for (size_t i = 0; i < free.size(); ++i) {
338 Vector3D exp = get_slaved_coords( mesh, free[i], err );
339 MSQ_ERRZERO(err);
340 exp -= free_coords[i];
341 if (exp.length() > EPSILON) {
342 std::cerr << "Slaved vertex " << (size_t)free[i] << " not at expected slaved location" << std::endl;
343 ++error_count;
344 }
345 }
346 }
347 // Otherwise, given the default input mesh, at least some of the
348 // vertices should be somewhere other than the slaved position
349 else {
350 int not_at_slaved_count = 0;
351
352 for (size_t i = 0; i < free.size(); ++i) {
353 Vector3D exp = get_slaved_coords( mesh, free[i], err );
354 MSQ_ERRZERO(err);
355 exp -= free_coords[i];
356 if (exp.length() > EPSILON) {
357 ++not_at_slaved_count;
358 }
359 }
360
361 if (0 == not_at_slaved_count) {
362 std::cerr << "All non-slaved vertices at slaved vertex locations" << std::endl;
363 error_count += free.size();
364 }
365 }
366
367 // expect all interior vertices to be at roughly the slaved location
368 if (mode != Settings::SLAVE_NONE) {
369 for (size_t i = 0; i < slave.size(); ++i) {
370 Vector3D exp = get_slaved_coords( mesh, slave[i], err );
371 MSQ_ERRZERO(err);
372 exp -= slave_coords[i];
373 if (exp.length() > EPSILON) {
374 std::cerr << "Interior vertex " << (size_t)slave[i] << " not at expected slaved location" << std::endl;
375 ++error_count;
376 }
377 }
378 }
379
380 return error_count;
381 }
382
compare_node_coords(Mesh & mesh1,Mesh & mesh2,MsqError & err)383 int compare_node_coords( Mesh& mesh1, Mesh& mesh2, MsqError& err )
384 {
385 const double EPSILON = 1e-4;
386
387 std::vector<Mesh::VertexHandle> verts1, verts2;
388 mesh1.get_all_vertices( verts1, err ); MSQ_ERRZERO(err);
389 mesh2.get_all_vertices( verts2, err ); MSQ_ERRZERO(err);
390 std::vector<MsqVertex> coords1(verts1.size()), coords2(verts2.size());
391 mesh1.vertices_get_coordinates( arrptr(verts1), arrptr(coords1), verts1.size(), err );
392 MSQ_ERRZERO(err);
393 mesh2.vertices_get_coordinates( arrptr(verts2), arrptr(coords2), verts2.size(), err );
394 MSQ_ERRZERO(err);
395
396 int error_count = 0;
397 assert(verts1.size() == verts2.size());
398 for (size_t i = 0; i < verts1.size(); ++i) {
399 assert( verts1[i] == verts2[i] );
400 Vector3D diff = coords1[i] - coords2[i];
401 if (diff.length() > EPSILON) {
402 std::cerr << "Vertex coordinates differ between calculated and flagged "
403 "meshes for vertex " << (size_t)verts1[i] << std::endl;
404 ++error_count;
405 }
406 }
407
408 return error_count;
409 }
410
check_no_slaved_corners(Mesh & mesh,MsqError & err)411 int check_no_slaved_corners( Mesh& mesh, MsqError& err )
412 {
413 std::vector<Mesh::ElementHandle> elems;
414 std::vector<Mesh::VertexHandle> verts;
415 std::vector<EntityTopology> types;
416 std::vector<size_t> offsets;
417 mesh.get_all_elements( elems, err ); MSQ_ERRZERO(err);
418 types.resize( elems.size() );
419 mesh.elements_get_topologies( arrptr(elems), arrptr(types), elems.size(), err ); MSQ_ERRZERO(err);
420 mesh.elements_get_attached_vertices( arrptr(elems), elems.size(), verts, offsets, err );
421 MSQ_ERRZERO(err);
422
423 std::vector<bool> slaved;
424 mesh.vertices_get_slaved_flag( arrptr(verts), slaved, verts.size(), err );
425 MSQ_ERRZERO(err);
426
427 int error_count = 0;
428 for (size_t i = 0; i < elems.size(); ++i) {
429 unsigned n = TopologyInfo::corners( types[i] );
430 for (unsigned j = 0; j < n; ++j) {
431 if (slaved[offsets[i]+j]) {
432 std::cerr << "Element " << (size_t)elems[i] << " corner " << j << " is slaved" <<std::endl;
433 ++error_count;
434 }
435 }
436 }
437
438 return error_count == 0;
439 }
440
check_global_patch_slaved(Mesh & mesh,MsqError & err)441 int check_global_patch_slaved( Mesh& mesh, MsqError& err )
442 {
443 Settings s;
444 s.set_slaved_ho_node_mode( Settings::SLAVE_FLAG );
445 MeshDomainAssoc mesh_and_domain = MeshDomainAssoc(&mesh, 0);
446 Instruction::initialize_vertex_byte( &mesh_and_domain, &s, err ); MSQ_ERRZERO(err);
447
448 PatchData pd;
449 pd.attach_settings( &s );
450 pd.set_mesh( &mesh );
451 pd.fill_global_patch( err ); MSQ_ERRZERO(err);
452
453 std::vector<bool> fixed, slaved;
454 mesh.vertices_get_fixed_flag( pd.get_vertex_handles_array(),
455 fixed, pd.num_nodes(), err );
456 MSQ_ERRZERO(err);
457 mesh.vertices_get_slaved_flag( pd.get_vertex_handles_array(),
458 slaved, pd.num_nodes(), err );
459 MSQ_ERRZERO(err);
460
461 const size_t first_free = 0;
462 const size_t first_slaved = pd.num_free_vertices();
463 const size_t first_fixed = pd.num_free_vertices() + pd.num_slave_vertices();
464 int error_count = 0;
465 for (size_t i = first_free; i < first_slaved; ++i) {
466 if (fixed[i]) {
467 std::cerr << "Vertex " << (size_t)pd.get_vertex_handles_array()[i]
468 << " is fixed in mesh but free in PatchData" << std::endl;
469 ++error_count;
470 }
471 if (slaved[i]) {
472 std::cerr << "Vertex " << (size_t)pd.get_vertex_handles_array()[i]
473 << " is slaved in mesh but free in PatchData" << std::endl;
474 ++error_count;
475 }
476 }
477 for (size_t i = first_slaved; i < first_fixed; ++i) {
478 if (fixed[i]) {
479 std::cerr << "Vertex " << (size_t)pd.get_vertex_handles_array()[i]
480 << " is fixed in mesh but slaved in PatchData" << std::endl;
481 ++error_count;
482 }
483 else if (!slaved[i]) {
484 std::cerr << "Vertex " << (size_t)pd.get_vertex_handles_array()[i]
485 << " is free in Mesh but slaved in PatchData" << std::endl;
486 ++error_count;
487 }
488 }
489 for (size_t i = first_fixed; i < pd.num_nodes(); ++i) {
490 if (!fixed[i]) {
491 std::cerr << "Vertex " << (size_t)pd.get_vertex_handles_array()[i]
492 << " is not fixed in mesh but is in PatchData" << std::endl;
493 ++error_count;
494 }
495 }
496 return 0 == error_count;
497 }
498
tag_patch_slaved(Mesh & mesh,Settings::HigherOrderSlaveMode mode,MsqError & err)499 void tag_patch_slaved( Mesh& mesh,
500 Settings::HigherOrderSlaveMode mode,
501 MsqError& err )
502 {
503 int zero = 0;
504 TagHandle tag = mesh.tag_create( "pd_slaved", Mesh::INT, 1, &zero, err );
505 MSQ_ERRRTN(err);
506
507 Settings s;
508 s.set_slaved_ho_node_mode( mode );
509 PatchData pd;
510 pd.attach_settings( &s );
511 pd.set_mesh( &mesh );
512 pd.fill_global_patch( err );
513 MSQ_ERRRTN(err);
514
515 const Mesh::VertexHandle* verts = pd.get_vertex_handles_array() + pd.num_free_vertices();
516 std::vector<int> ones( pd.num_slave_vertices(), 1 );
517 mesh.tag_set_vertex_data( tag, pd.num_slave_vertices(), verts, arrptr(ones), err );
518 MSQ_ERRRTN(err);
519 }
520