MOAB: Mesh Oriented datABase  (version 5.2.1)
high_aspect_ratio_test.cpp File Reference
#include "Mesquite.hpp"
#include "MsqError.hpp"
#include "MeshImpl.hpp"
#include "XYRectangle.hpp"
#include "ReferenceMesh.hpp"
#include "RefMeshTargetCalculator.hpp"
#include "TShapeB1.hpp"
#include "TQualityMetric.hpp"
#include "TerminationCriterion.hpp"
#include "QualityAssessor.hpp"
#include "PMeanPTemplate.hpp"
#include "FeasibleNewton.hpp"
#include "ConjugateGradient.hpp"
#include "InstructionQueue.hpp"
#include <iostream>
#include <fstream>
#include <cstdlib>
#include <cstdio>
+ Include dependency graph for high_aspect_ratio_test.cpp:

Go to the source code of this file.

Classes

struct  MeshParams
 Struct in which to store mesh description. More...

Defines

#define CHECKERR

Enumerations

enum  ExitCodes {
  NO_ERROR = 0, USAGE_ERROR = 1, NOT_AT_TARGET = 2, FAR_FROM_TARGET = 3,
  WRONG_DIRECTION = 4, DEGENERATE_ELEMENT = 5, INVERTED_ELEMENT = 6, LAST_EXIT_CODE = INVERTED_ELEMENT
}
enum  ParseState {
  OPEN, EXPECTING_M, EXPECTING_R, EXPECTING_O,
  EXPECTING_I
}

Functions

std::ostream & operator<< (std::ostream &str, const MeshParams &p)
void usage (const char *argv0, bool brief=true)
void create_input_mesh (const MeshParams &params, bool all_fixed, MeshImpl &mesh, MsqError &err)
void parse_options (char *argv[], int argc, MeshParams &mesh, MeshParams &ref, std::string &output_file, bool &fixed_boundary, TerminationCriterion::TimeStepFileType &write_timesteps, bool &use_feas_newt, int &num_iterations)
std::string base_name (std::string filename)
int main (int argc, char *argv[])
void parse_mesh_params (const char *argv, const char *arg, MeshParams &result)

Variables

const char default_out_file [] = "high_aspect.vtk"
 Default values for parameters.
const MeshParams default_mesh = { 0.6, 0.3, 2.0, 1.0 }
const MeshParams default_ref = { 0.5, 0.5, 1.0, 1.0 }
const int INNER_ITERATES = 1
const int OUTER_ITERATES = 10
const char * temp_file = "high_aspect_input.vtk"

Define Documentation

#define CHECKERR
Value:
if( err )                                                                                           \
    {                                                                                                   \
        std::cerr << "Internal error at line " << __LINE__ << ":" << std::endl << ( err ) << std::endl; \
        std::exit( LAST_EXIT_CODE + ( err ).error_code() );                                             \
    }

Definition at line 115 of file high_aspect_ratio_test.cpp.

Referenced by main().


Enumeration Type Documentation

enum ExitCodes
Enumerator:
NO_ERROR 
USAGE_ERROR 
NOT_AT_TARGET 
FAR_FROM_TARGET 
WRONG_DIRECTION 
DEGENERATE_ELEMENT 
INVERTED_ELEMENT 
LAST_EXIT_CODE 

Definition at line 74 of file high_aspect_ratio_test.cpp.

enum ParseState
Enumerator:
OPEN 
EXPECTING_M 
EXPECTING_R 
EXPECTING_O 
EXPECTING_I 

Definition at line 284 of file high_aspect_ratio_test.cpp.


Function Documentation

std::string base_name ( std::string  filename)

Definition at line 418 of file high_aspect_ratio_test.cpp.

References filename.

Referenced by main().

{
    std::string::size_type i = filename.rfind( "." );
    if( !i || i == std::string::npos )
        return filename;
    else
        return filename.substr( 0, i );
}
void create_input_mesh ( const MeshParams params,
bool  all_fixed,
MeshImpl mesh,
MsqError err 
)

Definition at line 379 of file high_aspect_ratio_test.cpp.

References moab::F, MeshParams::h, MSQ_CHKERR, MBMesquite::MeshImpl::read_vtk(), temp_file, MeshParams::w, MeshParams::x, MeshParams::y, and z.

Referenced by main().

{
    const double z = 0;
    const int F    = all_fixed;
    std::ofstream vtkfile( temp_file );
    double bx = all_fixed ? 0.5 * p.w : p.x;
    double by = all_fixed ? 0.5 * p.h : p.y;
    vtkfile << "# vtk DataFile Version 3.0" << std::endl
            << "Mesquite High Aspect Ratio test" << std::endl
            << "ASCII" << std::endl
            << "DATASET UNSTRUCTURED_GRID" << std::endl
            << "POINTS 9 float" << std::endl
            << 0.0 << ' ' << 0.0 << ' ' << z << std::endl
            << bx << ' ' << 0.0 << ' ' << z << std::endl
            << p.w << ' ' << 0.0 << ' ' << z << std::endl
            << 0.0 << ' ' << by << ' ' << z << std::endl
            << p.x << ' ' << p.y << ' ' << z << std::endl
            << p.w << ' ' << by << ' ' << z << std::endl
            << 0.0 << ' ' << p.h << ' ' << z << std::endl
            << bx << ' ' << p.h << ' ' << z << std::endl
            << p.w << ' ' << p.h << ' ' << z << std::endl
            << "CELLS 4 20" << std::endl
            << "4 0 1 4 3" << std::endl
            << "4 1 2 5 4" << std::endl
            << "4 4 5 8 7" << std::endl
            << "4 3 4 7 6" << std::endl
            << "CELL_TYPES 4" << std::endl
            << "9 9 9 9" << std::endl
            << "POINT_DATA 9" << std::endl
            << "SCALARS fixed int" << std::endl
            << "LOOKUP_TABLE default" << std::endl
            << "1 " << F << " 1" << std::endl
            << F << " 0 " << F << std::endl
            << "1 " << F << " 1" << std::endl;
    vtkfile.close();
    mesh.read_vtk( temp_file, err );
    std::remove( temp_file );MSQ_CHKERR( err );
}
int main ( int  argc,
char *  argv[] 
)

Definition at line 148 of file high_aspect_ratio_test.cpp.

References MBMesquite::TerminationCriterion::add_iteration_limit(), MBMesquite::InstructionQueue::add_quality_assessor(), MBMesquite::arrptr(), base_name(), CHECKERR, create_input_mesh(), DEGENERATE_ELEMENT, dof(), MBMesquite::XYRectangle::domain_DoF(), EPS, FAR_FROM_TARGET, MBMesquite::MeshImpl::get_all_vertices(), MBMesquite::QualityAssessor::get_inverted_element_count(), MeshParams::h, init(), MBMesquite::inner(), INNER_ITERATES, MBMesquite::inv(), INVERTED_ELEMENT, MBMesquite::length(), mesh, NOT_AT_TARGET, MBMesquite::TerminationCriterion::NOTYPE, MBMesquite::outer(), parse_options(), MBMesquite::IQInterface::run_instructions(), MBMesquite::QualityImprover::set_inner_termination_criterion(), MBMesquite::InstructionQueue::set_master_quality_improver(), MBMesquite::QualityImprover::set_outer_termination_criterion(), MBMesquite::XYRectangle::setup(), USAGE_ERROR, MBMesquite::PatchSetUser::use_element_on_vertex_patch(), MBMesquite::MeshImpl::vertices_get_coordinates(), MeshParams::w, MBMesquite::TerminationCriterion::write_mesh_steps(), MBMesquite::MeshImpl::write_vtk(), WRONG_DIRECTION, MeshParams::x, and MeshParams::y.

{
    MeshParams input_params, reference_params;
    bool fixed_boundary_vertices, feas_newt_solver;
    TerminationCriterion::TimeStepFileType write_timestep_files;
    std::string output_file_name;
    int num_iterations;

    parse_options( argv, argc, input_params, reference_params, output_file_name, fixed_boundary_vertices,
                   write_timestep_files, feas_newt_solver, num_iterations );

    MsqError err;
    MeshImpl mesh, refmesh;
    XYRectangle domain( input_params.w, input_params.h );
    create_input_mesh( input_params, fixed_boundary_vertices, mesh, err );
    CHECKERR
    create_input_mesh( reference_params, fixed_boundary_vertices, refmesh, err );
    CHECKERR
    domain.setup( &mesh, err );
    CHECKERR

    ReferenceMesh rmesh( &refmesh );
    RefMeshTargetCalculator tc( &rmesh );
    TShapeB1 tm;
    TQualityMetric qm( &tc, &tm );

    PMeanPTemplate of( 1.0, &qm );
    ConjugateGradient cg( &of );
    cg.use_element_on_vertex_patch();
    FeasibleNewton fn( &of );
    fn.use_element_on_vertex_patch();
    VertexMover* solver = feas_newt_solver ? (VertexMover*)&fn : (VertexMover*)&cg;

    TerminationCriterion inner, outer;
    inner.add_iteration_limit( INNER_ITERATES );
    outer.add_iteration_limit( num_iterations );
    if( write_timestep_files != TerminationCriterion::NOTYPE )
        outer.write_mesh_steps( base_name( output_file_name ).c_str(), write_timestep_files );
    solver->set_inner_termination_criterion( &inner );
    solver->set_outer_termination_criterion( &outer );

    QualityAssessor qa( &qm );
    InstructionQueue q;
    q.add_quality_assessor( &qa, err );
    q.set_master_quality_improver( solver, err );
    q.add_quality_assessor( &qa, err );

    MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, &domain );
    q.run_instructions( &mesh_and_domain, err );
    CHECKERR

    mesh.write_vtk( output_file_name.c_str(), err );
    CHECKERR

    // check for inverted elements
    int inv, unk;
    qa.get_inverted_element_count( inv, unk, err );
    if( inv )
    {
        std::cerr << inv << " inverted elements in final mesh" << std::endl;
        return INVERTED_ELEMENT;
    }
    else if( unk )
    {
        std::cerr << unk << " degenerate elements in final mesh" << std::endl;
        return DEGENERATE_ELEMENT;
    }

    // find the free vertex
    std::vector< Mesh::VertexHandle > vertices;
    mesh.get_all_vertices( vertices, err );
    if( vertices.empty() )
    {
        std::cerr << "Mesh contains no vertices" << std::endl;
        return USAGE_ERROR;
    }
    std::vector< unsigned short > dof( vertices.size(), 0 );
    domain.domain_DoF( arrptr( vertices ), arrptr( dof ), vertices.size(), err );
    CHECKERR
    int idx                              = std::find( dof.begin(), dof.end(), 2 ) - dof.begin();
    const Mesh::VertexHandle free_vertex = vertices[idx];
    MsqVertex coords;
    mesh.vertices_get_coordinates( &free_vertex, &coords, 1, err );
    CHECKERR

    // calculate optimal position for vertex
    const double xf = reference_params.x / reference_params.w;
    const double yf = reference_params.y / reference_params.h;
    Vector3D expect( xf * input_params.w, yf * input_params.h, 0 );

    // Test that we aren't further from the expected location
    // than when we started.
    const Vector3D init( input_params.x, input_params.y, 0 );
    if( ( coords - expect ).length() > ( init - expect ).length() )
    {
        std::cerr << "Vertex moved away from expected optimal position: "
                  << "(" << coords[0] << ", " << coords[1] << std::endl;
        return WRONG_DIRECTION;
    }

    // check if vertex moved MIN_FRACT of the way from the original position
    // to the desired one in the allowed iterations
    const double MIN_FRACT = 0.2;  // 20% of the way in 10 iterations
    const double fract     = ( coords - init ).length() / ( expect - init ).length();
    if( fract < MIN_FRACT )
    {
        std::cerr << "Vertex far from optimimal location" << std::endl
                  << "  Expected: (" << expect[0] << ", " << expect[1] << ", " << expect[2] << ")" << std::endl
                  << "  Actual:   (" << coords[0] << ", " << coords[1] << ", " << coords[2] << ")" << std::endl;
        return FAR_FROM_TARGET;
    }

    // check if vertex is at destired location
    const double EPS = 5e-2;  // allow a little leway
    if( fabs( coords[0] - expect[0] ) > EPS * input_params.w || fabs( coords[1] - expect[1] ) > EPS * input_params.h ||
        fabs( expect[2] ) > EPS )
    {
        std::cerr << "Vertex not at optimimal location" << std::endl
                  << "  Expected: (" << expect[0] << ", " << expect[1] << ", " << expect[2] << ")" << std::endl
                  << "  Actual:   (" << coords[0] << ", " << coords[1] << ", " << coords[2] << ")" << std::endl;
        return NOT_AT_TARGET;
    }

    return 0;
}
std::ostream& operator<< ( std::ostream &  str,
const MeshParams p 
)

Definition at line 61 of file high_aspect_ratio_test.cpp.

References MeshParams::h, MeshParams::w, MeshParams::x, and MeshParams::y.

{
    return str << p.x << ',' << p.y << ',' << p.w << ',' << p.h;
}
void parse_mesh_params ( const char *  argv,
const char *  arg,
MeshParams result 
)

Definition at line 274 of file high_aspect_ratio_test.cpp.

References MeshParams::h, usage, MeshParams::w, MeshParams::x, and MeshParams::y.

Referenced by parse_options().

{
    int c = std::sscanf( arg, "%lf,%lf,%lf,%lf", &result.x, &result.y, &result.w, &result.h );
    if( c != 2 && c != 4 )
    {
        std::cerr << "Error parsing mesh dimensions: \"" << arg << '"' << std::endl;
        usage( argv );
    }
}
void parse_options ( char *  argv[],
int  argc,
MeshParams mesh,
MeshParams ref,
std::string &  output_file,
bool &  fixed_boundary,
TerminationCriterion::TimeStepFileType write_timesteps,
bool &  use_feas_newt,
int &  num_iterations 
)

Definition at line 292 of file high_aspect_ratio_test.cpp.

References default_mesh, default_out_file, default_ref, EXPECTING_I, EXPECTING_M, EXPECTING_O, EXPECTING_R, MBMesquite::TerminationCriterion::GNUPLOT, MBMesquite::TerminationCriterion::NOTYPE, OPEN, OUTER_ITERATES, parse_mesh_params(), usage, and MBMesquite::TerminationCriterion::VTK.

Referenced by main().

{
    // begin with defaults
    mesh             = default_mesh;
    ref              = default_ref;
    output_file      = default_out_file;
    fixed_boundary   = false;
    write_timesteps  = TerminationCriterion::NOTYPE;
    feas_newt_solver = false;
    num_iterations   = OUTER_ITERATES;

    // parse CLI args
    ParseState state = OPEN;
    for( int i = 1; i < argc; ++i )
    {
        switch( state )
        {
            case EXPECTING_M:
                parse_mesh_params( argv[0], argv[i], mesh );
                state = OPEN;
                break;
            case EXPECTING_R:
                parse_mesh_params( argv[0], argv[i], ref );
                state = OPEN;
                break;
            case EXPECTING_O:
                output_file = argv[i];
                state       = OPEN;
                break;
            case EXPECTING_I:
                num_iterations = atoi( argv[i] );
                state          = OPEN;
                break;
            case OPEN:
                if( argv[i][0] != '-' || argv[i][1] == '\0' || argv[i][2] != '\0' )
                {
                    std::cerr << "Unexpected argument: \"" << argv[i] << '"' << std::endl;
                    usage( argv[0] );
                }

                switch( argv[i][1] )
                {
                    default:
                        usage( argv[0], true );
                        break;
                    case 'h':
                        usage( argv[0], false );
                        break;
                    case 'o':
                        state = EXPECTING_O;
                        break;
                    case 'f':
                        fixed_boundary = true;
                        break;
                    case 'F':
                        fixed_boundary = false;
                        break;
                    case 't':
                        write_timesteps = TerminationCriterion::VTK;
                        break;
                    case 'T':
                        write_timesteps = TerminationCriterion::GNUPLOT;
                        break;
                    case 'm':
                        state = EXPECTING_M;
                        break;
                    case 'r':
                        state = EXPECTING_R;
                        break;
                    case 'n':
                        feas_newt_solver = true;
                        break;
                    case 'c':
                        feas_newt_solver = false;
                        break;
                    case 'i':
                        state = EXPECTING_I;
                        break;
                }
                break;
        }
    }
}
void usage ( const char *  argv0,
bool  brief = true 
)

Definition at line 86 of file high_aspect_ratio_test.cpp.

References default_mesh, default_out_file, default_ref, NO_ERROR, OUTER_ITERATES, and USAGE_ERROR.

{
    std::ostream& str = brief ? std::cerr : std::cout;

    str << "Usage: " << argv0 << " [-o <output_file>]"
        << " [-f|-F] [-t|-T] [-n|-c] [-i <n>]"
        << " [-m <x>,<y>[,<w>,<h>]]"
        << " [-r <x>,<y>[,<w>,<h>]]" << std::endl;
    if( brief )
    {
        str << "       " << argv0 << " -h" << std::endl;
        std::exit( USAGE_ERROR );
    }

    str << "  -o  Specify output file (default is \"" << default_out_file << "\")" << std::endl
        << "  -f  Fixed boundary vertices" << std::endl
        << "  -F  Free boundary vertices (default)" << std::endl
        << "  -t  Write VTK timesteps" << std::endl
        << "  -T  Write GNUPlot timesteps" << std::endl
        << "  -m  Specify input mesh parameters (default " << default_mesh << ")" << std::endl
        << "  -r  Specify reference mesh parameters (default " << default_ref << ")" << std::endl
        << "  -n  Use FeasibleNewton solver" << std::endl
        << "  -c  Use ConjugateGradient solver (default)" << std::endl
        << "  -i  Specify number of iterations (default:" << OUTER_ITERATES << ")" << std::endl
        << std::endl;

    std::exit( NO_ERROR );
}

Variable Documentation

const MeshParams default_mesh = { 0.6, 0.3, 2.0, 1.0 }

Definition at line 68 of file high_aspect_ratio_test.cpp.

Referenced by parse_options(), and usage().

const char default_out_file[] = "high_aspect.vtk"

Default values for parameters.

Definition at line 67 of file high_aspect_ratio_test.cpp.

Referenced by parse_options(), and usage().

const MeshParams default_ref = { 0.5, 0.5, 1.0, 1.0 }

Definition at line 69 of file high_aspect_ratio_test.cpp.

Referenced by parse_options(), and usage().

const int INNER_ITERATES = 1

Definition at line 71 of file high_aspect_ratio_test.cpp.

Referenced by main().

const int OUTER_ITERATES = 10

Definition at line 72 of file high_aspect_ratio_test.cpp.

Referenced by parse_options(), and usage().

const char* temp_file = "high_aspect_input.vtk"

Definition at line 378 of file high_aspect_ratio_test.cpp.

Referenced by create_input_mesh().

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines