MOAB: Mesh Oriented datABase  (version 5.4.1)
2d_formulation_test.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2007 Sandia National Laboratories.  Developed at the
00005     University of Wisconsin--Madison under SNL contract number
00006     624796.  The U.S. Government and the University of Wisconsin
00007     retain certain rights to this software.
00008 
00009     This library is free software; you can redistribute it and/or
00010     modify it under the terms of the GNU Lesser General Public
00011     License as published by the Free Software Foundation; either
00012     version 2.1 of the License, or (at your option) any later version.
00013 
00014     This library is distributed in the hope that it will be useful,
00015     but WITHOUT ANY WARRANTY; without even the implied warranty of
00016     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017     Lesser General Public License for more details.
00018 
00019     You should have received a copy of the GNU Lesser General Public License
00020     (lgpl.txt) along with this library; if not, write to the Free Software
00021     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 
00023     (2007) [email protected]
00024 
00025   ***************************************************************** */
00026 
00027 /** \file main.cpp
00028  *  \brief Try examples from "Formulation of a Target-Matrix Paradigm
00029  *         for Mesh Optimization", Patrick Knupp.
00030  *  \author Jason Kraftcheck
00031  */
00032 
00033 #define USE_GLOBAL_PATCH
00034 
00035 #include "TestUtil.hpp"
00036 #include "Mesquite.hpp"
00037 
00038 #include "PMeanPTemplate.hpp"
00039 #include "AffineMapMetric.hpp"
00040 #include "ConjugateGradient.hpp"
00041 #include "TerminationCriterion.hpp"
00042 #include "ElementPMeanP.hpp"
00043 #include "MsqError.hpp"
00044 #include "TSquared.hpp"
00045 #include "MeshImpl.hpp"
00046 #include "PlanarDomain.hpp"
00047 #include "InstructionQueue.hpp"
00048 #include "TargetCalculator.hpp"
00049 #include "MetricWeight.hpp"
00050 #include "InverseMetricWeight.hpp"
00051 #include "TargetWriter.hpp"
00052 #include "WeightReader.hpp"
00053 
00054 #include <iostream>
00055 #include <cstdlib>
00056 
00057 using namespace MBMesquite;
00058 using namespace std;
00059 
00060 const double epsilon     = 2e-2;
00061 const bool write_results = true;
00062 
00063 #define CHKERR( A )            \
00064     if( A )                    \
00065     {                          \
00066         cerr << ( A ) << endl; \
00067         exit( 1 );             \
00068     }
00069 
00070 enum Grouping
00071 {
00072     SAMPLE,
00073     ELEMENT,
00074     QUADRANT,
00075     HALF
00076 };
00077 enum Weight
00078 {
00079     UNIT,
00080     METRIC,
00081     INV_METRIC
00082 };
00083 
00084 class IdentityTarget : public TargetCalculator
00085 {
00086   public:
00087     bool get_3D_target( PatchData&, size_t, Sample, MsqMatrix< 3, 3 >& W_out, MsqError& )
00088     {
00089         W_out = MsqMatrix< 3, 3 >( 1.0 );
00090         return true;
00091     }
00092 
00093     bool get_surface_target( PatchData&, size_t, Sample, MsqMatrix< 3, 2 >& W_out, MsqError& )
00094     {
00095         W_out = MsqMatrix< 3, 2 >( 1.0 );
00096         return true;
00097     }
00098 
00099     bool get_2D_target( PatchData&, size_t, Sample, MsqMatrix< 2, 2 >& W_out, MsqError& )
00100     {
00101         W_out = MsqMatrix< 2, 2 >( 1.0 );
00102         return true;
00103     }
00104 
00105     bool have_surface_orient() const
00106     {
00107         return false;
00108     }
00109 };
00110 
00111 void run_test( Grouping grouping, int of_power, Weight w, const string filename )
00112 {
00113     MsqError err;
00114 
00115     IdentityTarget target;
00116     TSquared target_metric;
00117     AffineMapMetric qual_metric( &target, &target_metric );
00118     ElementPMeanP elem_metric( of_power, &qual_metric );
00119     QualityMetric* qm_ptr = ( grouping == ELEMENT ) ? (QualityMetric*)&elem_metric : (QualityMetric*)&qual_metric;
00120 
00121     PMeanPTemplate OF( of_power, qm_ptr );
00122     ConjugateGradient solver( &OF );
00123     TerminationCriterion tc;
00124     TerminationCriterion itc;
00125     tc.add_absolute_vertex_movement( 1e-4 );
00126     itc.add_iteration_limit( 2 );
00127 #ifdef USE_GLOBAL_PATCH
00128     solver.use_global_patch();
00129     solver.set_inner_termination_criterion( &tc );
00130 #else
00131     solver.use_element_on_vertex_patch();
00132     solver.set_inner_termination_criterion( &itc );
00133     solver.set_outer_termination_criterion( &tc );
00134 #endif
00135 
00136     MeshImpl mesh, expected_mesh;
00137     std::string initfname = std::string( STRINGIFY( SRCDIR ) ) + "/2d_formulation_initial.vtk";
00138     mesh.read_vtk( initfname.c_str(), err );CHKERR( err )
00139     //  expected_mesh.read_vtk( (filename + ".vtk").c_str(), err ); CHKERR(err)
00140 
00141     PlanarDomain plane( PlanarDomain::XY );
00142 
00143     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, &plane );
00144 
00145     MetricWeight mw( &qual_metric );
00146     InverseMetricWeight imw( &qual_metric );
00147     WeightReader reader;
00148     if( w == METRIC )
00149     {
00150         TargetWriter writer( 0, &mw );
00151         InstructionQueue tq;
00152         tq.add_target_calculator( &writer, err );
00153         tq.run_instructions( &mesh_and_domain, err );CHKERR( err );
00154         qual_metric.set_weight_calculator( &reader );
00155     }
00156     else if( w == INV_METRIC )
00157     {
00158         TargetWriter writer( 0, &imw );
00159         InstructionQueue tq;
00160         tq.add_target_calculator( &writer, err );
00161         tq.run_instructions( &mesh_and_domain, err );CHKERR( err );
00162         qual_metric.set_weight_calculator( &reader );
00163     }
00164 
00165     InstructionQueue q;
00166     q.set_master_quality_improver( &solver, err );
00167     q.run_instructions( &mesh_and_domain, err );CHKERR( err )
00168     /*
00169       vector<Mesh::VertexHandle> vemain.cpprts;
00170       vector<MsqVertex> mesh_coords, expected_coords;
00171       mesh.get_all_vertices( verts, err ); CHKERR(err)
00172       mesh_coords.resize(verts.size());
00173       mesh.vertices_get_coordinates( arrptr(verts), arrptr(mesh_coords), verts.size(), err
00174       );CHKERR(err) expected_mesh.get_all_vertices( verts, err ); CHKERR(err)
00175       expected_coords.resize(verts.size());
00176       expected_mesh.vertices_get_coordinates( arrptr(verts), arrptr(expected_coords), verts.size(),
00177       err ); CHKERR(err) if (expected_coords.size() != mesh_coords.size()) { cerr << "Invlid
00178       expected mesh.  Vertex count doesn't match" << endl; exit(1);
00179       }
00180 
00181       unsigned error_count = 0;
00182       for (size_t i = 0; i < mesh_coords.size(); ++i)
00183         if ((expected_coords[i] - mesh_coords[i]).length_squared() > epsilon*epsilon)
00184           ++error_count;
00185 
00186       if (!error_count)
00187         cout << filename << " : SUCCESS" << endl;
00188       else
00189         cout << filename << " : FAILURE (" << error_count
00190              << " vertices differ by more than " << epsilon << ")" << endl;
00191     */
00192     if( write_results ) mesh.write_vtk( ( filename + ".results.vtk" ).c_str(), err );CHKERR( err )
00193 }
00194 
00195 int main()
00196 {
00197     run_test( SAMPLE, 1, UNIT, "1-1" );
00198     run_test( SAMPLE, 2, UNIT, "1-2" );
00199     run_test( SAMPLE, 4, UNIT, "1-4" );
00200     run_test( SAMPLE, 8, UNIT, "1-8" );
00201 
00202     run_test( SAMPLE, 1, UNIT, "2-NW" );
00203     run_test( ELEMENT, 1, UNIT, "2-NE" );
00204 
00205     run_test( SAMPLE, 1, UNIT, "3-Left" );
00206     run_test( SAMPLE, 1, METRIC, "3-Mid" );
00207     run_test( SAMPLE, 1, INV_METRIC, "3-Right" );
00208 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines