MOAB: Mesh Oriented datABase  (version 5.4.1)
jacobi_test.cpp
Go to the documentation of this file.
00001 /**\file jacobi_1.cpp
00002  *\brief simple convergence test for Jacobi optimization
00003  *\author Jason Kraftcheck <kraftche@cae.wisc.edu>
00004  *
00005  *Create simple quad mesh and test for convergence.
00006  */
00007
00008 #include "ArrayMesh.hpp"
00009 #include "PlanarDomain.hpp"
00010 #include "IdealShapeTarget.hpp"
00011 #include "TShapeB1.hpp"
00012 #include "TQualityMetric.hpp"
00013 #include "PMeanPTemplate.hpp"
00014 #include "SteepestDescent.hpp"
00015 #include "TerminationCriterion.hpp"
00016 #include "InstructionQueue.hpp"
00017 #include "MsqError.hpp"
00018 #include "QualityAssessor.hpp"
00019 #include "MeshWriter.hpp"
00020
00021 #include <cstdio>
00022 #include <cstring>
00023 #include <cstdlib>
00024
00025 using namespace MBMesquite;
00026
00027 /* Mesh:
00028
00029   12---13---14---15
00030   |    |    |    |
00031   |    |    |    |
00032   8----9----10---11
00033   |    |    |    |
00034   |    |    |    |
00035   4----5----6----7
00036   |    |    |    |
00037   |    |    |    |
00038   0----1----2----3
00039
00040 */
00041
00042 double dist( double* p1, double* p2 )
00043 {
00044     double d1 = p1[0] - p2[0];
00045     double d2 = p1[1] - p2[1];
00046     double d3 = p1[2] - p2[2];
00047     return sqrt( d1 * d2 * d3 );
00048 }
00049
00050 void usage( const char* argv0 )
00051 {
00052     fprintf( stderr, "Usage: %s [-i <init_vtk_file>] [-f <final_vtk_file>]\n", argv0 );
00053     exit( 1 );
00054 }
00055
00056 int main( int argc, char* argv[] )
00057 {
00058     const char* initial_mesh_file = 0;
00059     const char* final_mesh_file   = 0;
00060     for( int i = 1; i < argc; ++i )
00061     {
00062         if( !strcmp( argv[i], "-f" ) )
00063         {
00064             ++i;
00065             if( i == argc ) usage( argv[0] );
00066             final_mesh_file = argv[i];
00067         }
00068         else if( !strcmp( argv[i], "-i" ) )
00069         {
00070             ++i;
00071             if( i == argc ) usage( argv[0] );
00072             initial_mesh_file = argv[i];
00073         }
00074         else
00075             usage( argv[i] );
00076     }
00077
00078     // create mesh
00079     const int intervals  = 3;
00080     const double perturb = 0.3;
00081     const int nvtx       = ( intervals + 1 ) * ( intervals + 1 );
00082     double coords[nvtx * 3];
00083     double exp_coords[nvtx * 3];
00084     int fixed[nvtx];
00085     for( int i = 0; i < nvtx; ++i )
00086     {
00087         double* c = coords + 3 * i;
00088         double* e = exp_coords + 3 * i;
00089         int row   = i / ( intervals + 1 );
00090         int col   = i % ( intervals + 1 );
00091         double xdelta, ydelta;
00092         if( row > 0 && row < intervals && col > 0 && col < intervals )
00093         {
00094             fixed[i] = 0;
00095             xdelta   = row % 2 ? -perturb : 0;
00096             ydelta   = col % 2 ? perturb : -perturb;
00097         }
00098         else
00099         {
00100             fixed[i] = 1;
00101             xdelta = ydelta = 0.0;
00102         }
00103         c[0] = col + xdelta;
00104         c[1] = row + ydelta;
00105         c[2] = 0.0;
00106         e[0] = col;
00107         e[1] = row;
00108         e[2] = 0.0;
00109     }
00110     const int nquad = intervals * intervals;
00111     unsigned long conn[nquad * 4];
00112     for( int i = 0; i < nquad; ++i )
00113     {
00114         unsigned long* c = conn + 4 * i;
00115         int row          = i / intervals;
00116         int col          = i % intervals;
00117         int n0           = ( intervals + 1 ) * row + col;
00118         c[0]             = n0;
00119         c[1]             = n0 + 1;
00120         c[2]             = n0 + intervals + 2;
00121         c[3]             = n0 + intervals + 1;
00122     }
00123     MsqPrintError err( std::cerr );
00124     ArrayMesh mesh( 3, nvtx, coords, fixed, nquad, QUADRILATERAL, conn );
00125     PlanarDomain zplane( PlanarDomain::XY );
00126     if( initial_mesh_file )
00127     {
00128         MeshWriter::write_vtk( &mesh, initial_mesh_file, err );
00129         if( err )
00130         {
00131             fprintf( stderr, "%s: failed to write file\n", initial_mesh_file );
00132             return 1;
00133         }
00134     }
00135
00136     // do optimization
00137     const double eps = 0.01;
00138     IdealShapeTarget w;
00139     TShapeB1 mu;
00140     TQualityMetric metric( &w, &mu );
00141     PMeanPTemplate func( 1.0, &metric );
00142     SteepestDescent solver( &func );
00143     solver.use_element_on_vertex_patch();
00144     solver.do_jacobi_optimization();
00145     TerminationCriterion inner, outer;
00146     inner.add_absolute_vertex_movement( 0.5 * eps );
00147     outer.add_absolute_vertex_movement( 0.5 * eps );
00148     QualityAssessor qa( &metric );
00149     InstructionQueue q;
00150     q.add_quality_assessor( &qa, err );
00151     q.set_master_quality_improver( &solver, err );
00152     q.add_quality_assessor( &qa, err );
00153     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, &zplane );
00154     q.run_instructions( &mesh_and_domain, err );
00155     if( err ) return 2;
00156     if( final_mesh_file )
00157     {
00158         MeshWriter::write_vtk( &mesh, final_mesh_file, err );
00159         if( err )
00160         {
00161             fprintf( stderr, "%s: failed to write file\n", final_mesh_file );
00162             return 1;
00163         }
00164     }
00165
00166     // check final vertex positions
00167     int invalid = 0;
00168     for( int i = 0; i < nvtx; ++i )
00169     {
00170         if( dist( coords + 3 * i, exp_coords + 3 * i ) > eps )
00171         {
00172             ++invalid;
00173             printf( "Vertex %d at (%f,%f,%f), expected at (%f,%f,%f)\n", i, coords[3 * i], coords[3 * i + 1],
00174                     coords[3 * i + 2], exp_coords[3 * i], exp_coords[3 * i + 1], exp_coords[3 * i + 2] );
00175         }
00176     }
00177
00178     return invalid ? 2 : 0;
00179 }