MOAB: Mesh Oriented datABase  (version 5.4.1)
hireconst_test.cpp File Reference
#include <iostream>
#include <string>
#include <sstream>
#include <ctime>
#include <vector>
#include <algorithm>
#include "moab/Core.hpp"
#include "moab/Range.hpp"
#include "moab/MeshTopoUtil.hpp"
#include "moab/NestedRefine.hpp"
#include "moab/DiscreteGeometry/DGMSolver.hpp"
#include "moab/DiscreteGeometry/HiReconstruction.hpp"
#include "TestUtil.hpp"
#include <cmath>
+ Include dependency graph for hireconst_test.cpp:

Go to the source code of this file.

Functions

ErrorCode load_meshset_hirec (const char *infile, Interface *mbimpl, EntityHandle &meshset, ParallelComm *&pc, const int degree=0, const int dim=2)
ErrorCode test_mesh (const char *infile, const int degree, const bool interp, const int dim)
void compute_linear_coords (const int nvpe, double *elemcoords, double *naturals, double *linearcoords)
ErrorCode create_unitsq_tris (Interface *mbImpl, size_t n, std::vector< EntityHandle > &tris)
ErrorCode create_unitsq_quads (Interface *mbImpl, size_t n, std::vector< EntityHandle > &quads)
ErrorCode test_unitsq_tris ()
ErrorCode test_unitsq_quads ()
ErrorCode test_unitsphere ()
ErrorCode test_unitcircle ()
ErrorCode exact_error_torus (const double R, const double r, const double center[3], int npnts, double *pnt, double &error_l1, double &error_l2, double &error_li)
ErrorCode project_exact_torus (Interface *mbImpl, EntityHandle meshset, int dim, const double R, const double r, const double center[3])
int main (int argc, char *argv[])
ErrorCode project_exact_torus (Interface *mbImpl, EntityHandle meshset, int dim, const double R, const double r, const double center[])
ErrorCode exact_error_torus (const double R, const double r, const double center[], int npnts, double *pnts, double &error_l1, double &error_l2, double &error_li)

Function Documentation

void compute_linear_coords ( const int  nvpe,
double *  elemcoords,
double *  naturals,
double *  linearcoords 
)

Definition at line 418 of file hireconst_test.cpp.

Referenced by test_mesh(), test_unitcircle(), test_unitsphere(), test_unitsq_quads(), and test_unitsq_tris().

{
    assert( elemcoords && linearcoords );

    for( int i = 0; i < 3; ++i )
    {
        linearcoords[i] = 0;

        for( int j = 0; j < nvpe; ++j )
        {
            linearcoords[i] += naturals[j] * elemcoords[3 * j + i];
        }
    }
}
ErrorCode create_unitsq_quads ( Interface mbImpl,
size_t  n,
std::vector< EntityHandle > &  quads 
)

Definition at line 318 of file hireconst_test.cpp.

References moab::Interface::create_element(), moab::Interface::create_vertex(), moab::error(), ErrorCode, MB_CHK_ERR, MB_SET_ERR, and MBQUAD.

Referenced by test_unitsq_quads().

{
    if( n < 2 )
    {
        MB_SET_ERR( MB_FAILURE, "n must be at least 2" );
    }

    ErrorCode error;
    std::vector< EntityHandle > verts( n * n );
    size_t istr = quads.size();
    quads.resize( istr + ( n - 1 ) * ( n - 1 ) );
    double istep = 1.0 / (double)( n - 1 );

    for( size_t j = 0; j < n; ++j )
    {
        for( size_t i = 0; i < n; ++i )
        {
            double coord[3] = { i * istep, j * istep, 0 };
            error           = mbImpl->create_vertex( coord, verts[j * n + i] );MB_CHK_ERR( error );
        }
    }

    for( size_t jj = 0; jj < n - 1; ++jj )
    {
        for( size_t ii = 0; ii < n - 1; ++ii )
        {
            EntityHandle conn[4] = { verts[jj * n + ii], verts[jj * n + ii + 1], verts[( jj + 1 ) * n + ii + 1],
                                     verts[( jj + 1 ) * n + ii] };
            error                = mbImpl->create_element( MBQUAD, conn, 4, quads[istr + jj * ( n - 1 ) + ii] );MB_CHK_ERR( error );
        }
    }

    return error;
}
ErrorCode create_unitsq_tris ( Interface mbImpl,
size_t  n,
std::vector< EntityHandle > &  tris 
)

Definition at line 278 of file hireconst_test.cpp.

References moab::Interface::create_element(), moab::Interface::create_vertex(), moab::error(), ErrorCode, MB_CHK_ERR, MB_SET_ERR, and MBTRI.

Referenced by test_unitsq_tris().

{
    if( n < 2 )
    {
        MB_SET_ERR( MB_FAILURE, "n must be at least 2" );
    }

    ErrorCode error;
    std::vector< EntityHandle > verts( n * n );
    size_t istr = tris.size();
    tris.resize( istr + 2 * ( n - 1 ) * ( n - 1 ) );
    double istep = 1.0 / (double)( n - 1 );

    for( size_t j = 0; j < n; ++j )
    {
        for( size_t i = 0; i < n; ++i )
        {
            double coord[3] = { i * istep, j * istep, 0 };
            EntityHandle temp;
            error = mbImpl->create_vertex( coord, temp );MB_CHK_ERR( error );
            verts[j * n + i] = temp;
        }
    }

    for( size_t jj = 0; jj < n - 1; ++jj )
    {
        for( size_t ii = 0; ii < n - 1; ++ii )
        {
            EntityHandle conn[3] = { verts[jj * n + ii], verts[( jj + 1 ) * n + ii + 1], verts[( jj + 1 ) * n + ii] };
            error                = mbImpl->create_element( MBTRI, conn, 3, tris[istr + 2 * jj * ( n - 1 ) + 2 * ii] );MB_CHK_ERR( error );
            conn[0] = verts[jj * n + ii];
            conn[1] = verts[jj * n + ii + 1];
            conn[2] = verts[( jj + 1 ) * n + ii + 1];
            error   = mbImpl->create_element( MBTRI, conn, 3, tris[istr + 2 * jj * ( n - 1 ) + 2 * ii + 1] );MB_CHK_ERR( error );
        }
    }

    return error;
}
ErrorCode exact_error_torus ( const double  R,
const double  r,
const double  center[3],
int  npnts,
double *  pnt,
double &  error_l1,
double &  error_l2,
double &  error_li 
)

Referenced by test_mesh().

ErrorCode exact_error_torus ( const double  R,
const double  r,
const double  center[],
int  npnts,
double *  pnts,
double &  error_l1,
double &  error_l2,
double &  error_li 
)

Definition at line 691 of file hireconst_test.cpp.

References MB_SUCCESS.

{
    error_l1 = 0;
    error_l2 = 0;
    error_li = 0;
    double x = 0, y = 0, z = 0;
    double error_single = 0;

    for( int i = 0; i < npnts; i++ )
    {
        x            = pnts[3 * i] - center[0];
        y            = pnts[3 * i + 1] - center[1];
        z            = pnts[3 * i + 2] - center[2];
        error_single = fabs( sqrt( pow( sqrt( pow( x, 2 ) + pow( y, 2 ) ) - R, 2 ) + pow( z, 2 ) ) - r );
        error_l1     = error_l1 + error_single;
        error_l2     = error_l2 + error_single * error_single;

        if( error_li < error_single )
        {
            error_li = error_single;
        }
    }

    error_l1 = error_l1 / (double)npnts;
    error_l2 = sqrt( error_l2 / (double)npnts );
    return MB_SUCCESS;
}
ErrorCode load_meshset_hirec ( const char *  infile,
Interface mbimpl,
EntityHandle meshset,
ParallelComm *&  pc,
const int  degree = 0,
const int  dim = 2 
)

Definition at line 153 of file hireconst_test.cpp.

References moab::Interface::create_meshset(), moab::error(), ErrorCode, moab::HiReconstruction::estimate_num_ghost_layers(), moab::ParallelComm::get_pcomm(), moab::Interface::load_file(), MB_CHK_ERR, MB_SET_ERR, MESHSET_SET, MPI_COMM_WORLD, rank, and read_options.

Referenced by closedsurface_uref_hirec_convergence_study(), test_mesh(), test_unitcircle(), and test_unitsphere().

{
    ErrorCode error;
    error = mbimpl->create_meshset( moab::MESHSET_SET, meshset );MB_CHK_ERR( error );
#ifdef MOAB_HAVE_MPI
    int nprocs, rank;
    MPI_Comm comm = MPI_COMM_WORLD;
    MPI_Comm_size( comm, &nprocs );
    MPI_Comm_rank( comm, &rank );
    EntityHandle partnset;
    error = mbimpl->create_meshset( moab::MESHSET_SET, partnset );MB_CHK_ERR( error );

    if( nprocs > 1 )
    {
        pc = moab::ParallelComm::get_pcomm( mbimpl, partnset, &comm );
    }

    if( nprocs > 1 )
    {
        int nghlayers = degree > 0 ? HiReconstruction::estimate_num_ghost_layers( degree, true ) : 0;

        if( nghlayers )
        {
            // get ghost layers
            if( dim == 2 )
            {
                read_options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_"
                               "SHARED_ENTS;PARALLEL_GHOST=2.0.";
            }
            else if( dim == 1 )
            {
                read_options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_"
                               "SHARED_ENTS;PARALLEL_GHOST=1.0.";
            }
            else
            {
                MB_SET_ERR( MB_FAILURE, "unsupported dimension" );
            }

            read_options += ( '0' + nghlayers );
        }
        else
        {
            read_options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;";
        }

        error = mbimpl->load_file( infile, &meshset, read_options.c_str() );MB_CHK_ERR( error );
    }
    else
    {
        error = mbimpl->load_file( infile, &meshset );MB_CHK_ERR( error );
    }

#else
    assert( !pc && degree && dim );
    error = mbimpl->load_file( infile, &meshset );MB_CHK_ERR( error );
#endif
    return error;
}
int main ( int  argc,
char *  argv[] 
)

Definition at line 70 of file hireconst_test.cpp.

References dim, moab::error(), ErrorCode, MB_CHK_ERR, MPI_COMM_WORLD, rank, test_mesh(), test_unitcircle(), test_unitsphere(), test_unitsq_quads(), and test_unitsq_tris().

{
#ifdef MOAB_HAVE_MPI
    MPI_Init( &argc, &argv );
    int nprocs, rank;
    MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );
#endif
    std::string infile;
    int degree = 2, dim = 0;
    bool interp = false;
    ErrorCode error;

    if( argc == 1 )
    {
        error = test_unitsq_tris();MB_CHK_ERR( error );
        error = test_unitsq_quads();MB_CHK_ERR( error );
        error = test_unitsphere();MB_CHK_ERR( error );
        error = test_unitcircle();MB_CHK_ERR( error );
#ifdef MOAB_HAVE_MPI
        MPI_Finalize();
#endif
        return 0;
    }
    else
    {
        infile      = std::string( argv[1] );
        bool hasdim = false;

        for( int i = 2; i < argc; ++i )
        {
            if( i + 1 != argc )
            {
                if( std::string( argv[i] ) == "-degree" )
                {
                    degree = atoi( argv[++i] );
                }
                else if( std::string( argv[i] ) == "-interp" )
                {
                    interp = atoi( argv[++i] );
                }
                else if( std::string( argv[i] ) == "-dim" )
                {
                    dim    = atoi( argv[++i] );
                    hasdim = true;
                }
                else
                {
                    std::cout << argv[i] << std::endl;
                    std::cout << "usage: " << argv[0]
                              << " <mesh file> -degree <degree> -interp <0=least square, "
                                 "1=interpolation> -dim <mesh dimension>"
                              << std::endl;
                    return 0;
                }
            }
        }

        if( !hasdim )
        {
            std::cout << "Dimension of input mesh should be provided, positive and less than 3" << std::endl;
            return 0;
        }

        if( degree <= 0 || dim > 2 || dim <= 0 )
        {
            std::cout << "Input degree should be positive number;\n";
            std::cout << "Input dimesion should be positive and less than 3;" << std::endl;
            return -1;
        }

        std::cout << "Testing on " << infile << " with dimension " << dim << "\n";
        std::string opts = interp ? "interpolation" : "least square fitting";
        std::cout << "High order reconstruction with degree " << degree << " " << opts << std::endl;
    }

    error = test_mesh( infile.c_str(), degree, interp, dim );MB_CHK_ERR( error );
#ifdef MOAB_HAVE_MPI
    MPI_Finalize();
#endif
    return 0;
}
ErrorCode project_exact_torus ( Interface mbImpl,
EntityHandle  meshset,
int  dim,
const double  R,
const double  r,
const double  center[3] 
)

Referenced by test_mesh().

ErrorCode project_exact_torus ( Interface mbImpl,
EntityHandle  meshset,
int  dim,
const double  R,
const double  r,
const double  center[] 
)

Definition at line 653 of file hireconst_test.cpp.

References CHECK_ERR, moab::error(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), MB_SUCCESS, moab::Interface::set_coords(), and moab::Range::size().

{
    ErrorCode error;
    Range elems, verts;
    error = mbImpl->get_entities_by_dimension( meshset, dim, elems );CHECK_ERR( error );
    error = mbImpl->get_connectivity( elems, verts );CHECK_ERR( error );
    double pnts[3] = { 0, 0, 0 }, cnt[3] = { 0, 0, 0 }, nrms[3] = { 0, 0, 0 };
    double x = 0, y = 0, z = 0, d1 = 0, d2 = 0;

    for( int i = 0; i < (int)verts.size(); i++ )
    {
        EntityHandle v = verts[i];
        error          = mbImpl->get_coords( &v, 1, &pnts[0] );CHECK_ERR( error );
        x       = pnts[0] - center[0];
        y       = pnts[1] - center[0];
        z       = pnts[2] - center[2];
        d1      = sqrt( x * x + y * y );
        cnt[0]  = R * x / d1;
        cnt[1]  = R * y / d1;
        cnt[2]  = 0;
        d2      = sqrt( ( x - cnt[0] ) * ( x - cnt[0] ) + ( y - cnt[1] ) * ( y - cnt[1] ) + z * z );
        nrms[0] = ( x - cnt[0] ) / d2;
        nrms[1] = ( y - cnt[1] ) / d2;
        nrms[2] = z / d2;
        pnts[0] = cnt[0] + r * nrms[0];
        pnts[1] = cnt[1] + r * nrms[1];
        pnts[2] = cnt[2] + r * nrms[2];
        error   = mbImpl->set_coords( &v, 1, &pnts[0] );CHECK_ERR( error );
    }

    return MB_SUCCESS;
}
ErrorCode test_mesh ( const char *  infile,
const int  degree,
const bool  interp,
const int  dim 
)

Definition at line 218 of file hireconst_test.cpp.

References moab::Range::begin(), center(), CHECK_ERR, compute_linear_coords(), moab::Range::end(), moab::error(), ErrorCode, exact_error_torus(), moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), moab::HiReconstruction::hiproj_walf_in_element(), load_meshset_hirec(), MB_CHK_ERR, project_exact_torus(), moab::R, moab::HiReconstruction::reconstruct3D_curve_geom(), moab::HiReconstruction::reconstruct3D_surf_geom(), moab::Range::size(), and moab::DGMSolver::vec_distance().

Referenced by main(), test_1D(), test_2D(), and test_3D().

{
    Core moab;
    Interface* mbimpl = &moab;
    ParallelComm* pc  = NULL;
    EntityHandle meshset;
    // load mesh file
    ErrorCode error;
    error = load_meshset_hirec( infile, mbimpl, meshset, pc, degree, dim );MB_CHK_ERR( error );
    // project to exact surface: torus
    double center[3] = { 0, 0, 0 };
    double R = 1, r = 0.3;
    error = project_exact_torus( mbimpl, meshset, dim, R, r, center );CHECK_ERR( error );
    // initialize
    HiReconstruction hirec( dynamic_cast< Core* >( mbimpl ), pc, meshset );
    Range elems;
    error = mbimpl->get_entities_by_dimension( meshset, dim, elems );MB_CHK_ERR( error );

    // reconstruction
    if( dim == 2 )
    {
        // error = hirec.reconstruct3D_surf_geom(degree, interp, false); MB_CHK_ERR(error);
        error = hirec.reconstruct3D_surf_geom( degree, interp, true );MB_CHK_ERR( error );
    }
    else if( dim == 1 )
    {
        error = hirec.reconstruct3D_curve_geom( degree, interp, true );MB_CHK_ERR( error );
    }

    // fitting
    double mxdist = 0, errl1 = 0, errl2 = 0, errli = 0;
    double* pnts_proj = new double[elems.size() * 3];

    for( Range::iterator ielem = elems.begin(); ielem != elems.end(); ++ielem )
    {
        int nvpe;
        const EntityHandle* conn;
        error = mbimpl->get_connectivity( *ielem, conn, nvpe );MB_CHK_ERR( error );
        double w = 1.0 / (double)nvpe;
        std::vector< double > naturalcoords2fit( nvpe, w );
        double newcoords[3], linearcoords[3];
        error = hirec.hiproj_walf_in_element( *ielem, nvpe, 1, &( naturalcoords2fit[0] ), newcoords );MB_CHK_ERR( error );
        pnts_proj[3 * ( *ielem - *elems.begin() )]     = newcoords[0];
        pnts_proj[3 * ( *ielem - *elems.begin() ) + 1] = newcoords[1];
        pnts_proj[3 * ( *ielem - *elems.begin() ) + 2] = newcoords[2];
        std::vector< double > coords( 3 * nvpe );
        error = mbimpl->get_coords( conn, nvpe, &( coords[0] ) );MB_CHK_ERR( error );
        compute_linear_coords( nvpe, &( coords[0] ), &( naturalcoords2fit[0] ), linearcoords );
        mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
    }

    delete[] pnts_proj;
    // compute error for torus
    error = exact_error_torus( R, r, center, (int)elems.size(), pnts_proj, errl1, errl2, errli );MB_CHK_ERR( error );
    std::cout << "Errors using exact torus for degree " << degree << " fit : L1 = " << errl1 << ", L2 = " << errl2
              << ", Linf = " << errli << std::endl;
    std::cout << "Maximum projection lift is " << mxdist << std::endl;
    return error;
}

Definition at line 574 of file hireconst_test.cpp.

References moab::Range::begin(), compute_linear_coords(), dim, moab::Range::end(), moab::error(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), moab::HiReconstruction::hiproj_walf_in_element(), load_meshset_hirec(), MB_CHK_ERR, moab::HiReconstruction::reconstruct3D_curve_geom(), moab::DGMSolver::vec_2norm(), and moab::DGMSolver::vec_distance().

Referenced by main().

{
    // path to test files
    int nfiles              = 4;
    std::string filenames[] = { TestDir + "unittest/circle_3.vtk", TestDir + "unittest/circle_4.vtk",
                                TestDir + "unittest/circle_10.vtk", TestDir + "unittest/circle_20.vtk" };
    ErrorCode error;
    int maxdeg = 6;

    for( int ifile = 0; ifile < nfiles; ++ifile )
    {
        Core moab;
        Interface* mbimpl = &moab;
        ParallelComm* pc  = 0;
        EntityHandle meshset;
        int dim = 1;
        // load file
        error = load_meshset_hirec( filenames[ifile].c_str(), mbimpl, meshset, pc, maxdeg, dim );MB_CHK_ERR( error );
        // initialize
        HiReconstruction hirec( &moab, pc, meshset );
        Range edges;
        error = mbimpl->get_entities_by_dimension( meshset, dim, edges );MB_CHK_ERR( error );

        // reconstruction
        for( int degree = 1; degree <= maxdeg; ++degree )
        {
            hirec.reconstruct3D_curve_geom( degree, true, false, true );
            // fitting
            double mxdist = 0, mxerr = 0;

            for( Range::iterator iedge = edges.begin(); iedge != edges.end(); ++iedge )
            {
                int nvpe;
                const EntityHandle* conn;
                error = mbimpl->get_connectivity( *iedge, conn, nvpe );MB_CHK_ERR( error );
                double w = 1.0 / (double)nvpe;
                std::vector< double > naturalcoords2fit( nvpe, w );
                double newcoords[3], linearcoords[3];
                error = hirec.hiproj_walf_in_element( *iedge, nvpe, 1, &( naturalcoords2fit[0] ), newcoords );MB_CHK_ERR( error );
                std::vector< double > coords( 3 * nvpe );
                error = mbimpl->get_coords( conn, nvpe, &( coords[0] ) );MB_CHK_ERR( error );
                compute_linear_coords( nvpe, &( coords[0] ), &( naturalcoords2fit[0] ), linearcoords );
                mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
                mxerr  = std::max( mxerr, fabs( DGMSolver::vec_2norm( 3, newcoords ) - 1 ) );
            }

            std::cout << filenames[ifile] << ": unit circle"
                      << " degree= " << degree << " interpolation:\n";
            std::cout << "maximum projection lift is " << mxdist << ", maximum error is " << mxerr << std::endl;
            mxdist = 0;
            mxerr  = 0;
            hirec.reconstruct3D_curve_geom( degree, false, false, true );

            // fitting
            for( Range::iterator iedge = edges.begin(); iedge != edges.end(); ++iedge )
            {
                int nvpe;
                const EntityHandle* conn;
                error = mbimpl->get_connectivity( *iedge, conn, nvpe );MB_CHK_ERR( error );
                double w = 1.0 / (double)nvpe;
                std::vector< double > naturalcoords2fit( nvpe, w );
                double newcoords[3], linearcoords[3];
                error = hirec.hiproj_walf_in_element( *iedge, nvpe, 1, &( naturalcoords2fit[0] ), newcoords );MB_CHK_ERR( error );
                std::vector< double > coords( 3 * nvpe );
                error = mbimpl->get_coords( conn, nvpe, &( coords[0] ) );MB_CHK_ERR( error );
                compute_linear_coords( nvpe, &( coords[0] ), &( naturalcoords2fit[0] ), linearcoords );
                mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
                mxerr  = std::max( mxerr, fabs( DGMSolver::vec_2norm( 3, newcoords ) - 1 ) );
            }

            std::cout << filenames[ifile] << ": unit circle"
                      << " degree= " << degree << " least square:\n";
            std::cout << "maximum projection lift is " << mxdist << ", maximum error is " << mxerr << std::endl;
        }
    }

    return error;
}

Definition at line 496 of file hireconst_test.cpp.

References moab::Range::begin(), compute_linear_coords(), moab::Range::end(), moab::error(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), moab::HiReconstruction::hiproj_walf_in_element(), load_meshset_hirec(), MB_CHK_ERR, moab::HiReconstruction::reconstruct3D_surf_geom(), moab::DGMSolver::vec_2norm(), and moab::DGMSolver::vec_distance().

Referenced by main().

{
    // path to test files
    int nfiles              = 4;
    std::string filenames[] = { TestDir + "unittest/sphere_tris_5.vtk", TestDir + "unittest/sphere_tris_20.vtk",
                                TestDir + "unittest/sphere_quads_5.vtk", TestDir + "unittest/sphere_quads_20.vtk" };
    ErrorCode error;
    int maxdeg = 6;

    for( int ifile = 0; ifile < nfiles; ++ifile )
    {
        Core moab;
        Interface* mbimpl = &moab;
        ParallelComm* pc  = NULL;
        EntityHandle meshset;
        // load file
        error = load_meshset_hirec( filenames[ifile].c_str(), mbimpl, meshset, pc, maxdeg );MB_CHK_ERR( error );
        // initialize
        HiReconstruction hirec( &moab, pc, meshset );
        Range elems;
        error = mbimpl->get_entities_by_dimension( meshset, 2, elems );MB_CHK_ERR( error );

        // reconstruction
        for( int degree = 1; degree <= maxdeg; ++degree )
        {
            hirec.reconstruct3D_surf_geom( degree, true, false, true );
            // fitting
            double mxdist = 0, mxerr = 0;

            for( Range::iterator ielem = elems.begin(); ielem != elems.end(); ++ielem )
            {
                int nvpe;
                const EntityHandle* conn;
                error = mbimpl->get_connectivity( *ielem, conn, nvpe );MB_CHK_ERR( error );
                double w = 1.0 / (double)nvpe;
                std::vector< double > naturalcoords2fit( nvpe, w );
                double newcoords[3], linearcoords[3];
                error = hirec.hiproj_walf_in_element( *ielem, nvpe, 1, &( naturalcoords2fit[0] ), newcoords );MB_CHK_ERR( error );
                std::vector< double > coords( 3 * nvpe );
                error = mbimpl->get_coords( conn, nvpe, &( coords[0] ) );MB_CHK_ERR( error );
                compute_linear_coords( nvpe, &( coords[0] ), &( naturalcoords2fit[0] ), linearcoords );
                mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
                mxerr  = std::max( mxerr, fabs( DGMSolver::vec_2norm( 3, newcoords ) - 1 ) );
            }

            std::cout << filenames[ifile] << ": unit sphere"
                      << " degree= " << degree << " interpolation:\n";
            std::cout << "maximum projection lift is " << mxdist << ", maximum error is " << mxerr << std::endl;
            mxdist = 0;
            mxerr  = 0;
            hirec.reconstruct3D_surf_geom( degree, false, false, true );

            // fitting
            for( Range::iterator ielem = elems.begin(); ielem != elems.end(); ++ielem )
            {
                int nvpe;
                const EntityHandle* conn;
                error = mbimpl->get_connectivity( *ielem, conn, nvpe );MB_CHK_ERR( error );
                double w = 1.0 / (double)nvpe;
                std::vector< double > naturalcoords2fit( nvpe, w );
                double newcoords[3], linearcoords[3];
                error = hirec.hiproj_walf_in_element( *ielem, nvpe, 1, &( naturalcoords2fit[0] ), newcoords );MB_CHK_ERR( error );
                std::vector< double > coords( 3 * nvpe );
                error = mbimpl->get_coords( conn, nvpe, &( coords[0] ) );MB_CHK_ERR( error );
                compute_linear_coords( nvpe, &( coords[0] ), &( naturalcoords2fit[0] ), linearcoords );
                mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
                mxerr  = std::max( mxerr, fabs( DGMSolver::vec_2norm( 3, newcoords ) - 1 ) );
            }

            std::cout << filenames[ifile] << ": unit sphere"
                      << " degree= " << degree << " least square:\n";
            std::cout << "maximum projection lift is " << mxdist << ", maximum error is " << mxerr << std::endl;
        }
    }

    return error;
}

Definition at line 433 of file hireconst_test.cpp.

References compute_linear_coords(), create_unitsq_quads(), moab::error(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::HiReconstruction::hiproj_walf_in_element(), MB_CHK_ERR, moab::HiReconstruction::reconstruct3D_surf_geom(), and moab::DGMSolver::vec_distance().

Referenced by main().

{
    ErrorCode error;

    for( size_t n = 2; n <= 8; ++n )
    {
        Core moab;
        Interface* mbImpl = &moab;
        std::vector< EntityHandle > quads;
        error = create_unitsq_quads( mbImpl, n, quads );MB_CHK_ERR( error );
        EntityHandle meshIn = 0;
        HiReconstruction hirec( dynamic_cast< Core* >( mbImpl ), 0, meshIn );

        for( int degree = 1; degree <= 6; ++degree )
        {
            // reconstruct geometry, interpolation
            hirec.reconstruct3D_surf_geom( degree, true, false, true );
            // test fitting result
            double mxdist = 0;

            for( size_t iquad = 0; iquad < quads.size(); ++iquad )
            {
                const int nvpe                 = 4;
                double w                       = 1.0 / (double)nvpe;
                double naturalcoords2fit[nvpe] = { w, w, w, w }, newcoords[3];
                error = hirec.hiproj_walf_in_element( quads[iquad], nvpe, 1, naturalcoords2fit, newcoords );MB_CHK_ERR( error );
                std::vector< EntityHandle > conn;
                error = mbImpl->get_connectivity( &( quads[iquad] ), 1, conn );MB_CHK_ERR( error );
                double coords[3 * nvpe], linearcoords[3];
                error = mbImpl->get_coords( &( conn[0] ), nvpe, coords );MB_CHK_ERR( error );
                compute_linear_coords( nvpe, coords, naturalcoords2fit, linearcoords );
                mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
            }

            std::cout << "quadrilateral unit square n= " << n << " degree= " << degree << " interpolation:\n";
            std::cout << "maximum projection lift is " << mxdist << std::endl;
            mxdist = 0;
            // reconstruct geometry, least square fitting
            hirec.reconstruct3D_surf_geom( degree, false, false, true );

            // test fitting result
            for( size_t iquad = 0; iquad < quads.size(); ++iquad )
            {
                const int nvpe                 = 4;
                double w                       = 1.0 / (double)nvpe;
                double naturalcoords2fit[nvpe] = { w, w, w, w }, newcoords[3];
                error = hirec.hiproj_walf_in_element( quads[iquad], nvpe, 1, naturalcoords2fit, newcoords );MB_CHK_ERR( error );
                std::vector< EntityHandle > conn;
                error = mbImpl->get_connectivity( &( quads[iquad] ), 1, conn );MB_CHK_ERR( error );
                double coords[3 * nvpe], linearcoords[3];
                error = mbImpl->get_coords( &( conn[0] ), nvpe, coords );MB_CHK_ERR( error );
                compute_linear_coords( nvpe, coords, naturalcoords2fit, linearcoords );
                mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
            }

            std::cout << "quadrilateral unit square n= " << n << " degree= " << degree << " least square:\n";
            std::cout << "maximum projection lift is " << mxdist << std::endl;
        }
    }

    return error;
}

Definition at line 353 of file hireconst_test.cpp.

References compute_linear_coords(), create_unitsq_tris(), moab::error(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::HiReconstruction::hiproj_walf_in_element(), MB_CHK_ERR, moab::HiReconstruction::reconstruct3D_surf_geom(), and moab::DGMSolver::vec_distance().

Referenced by main().

{
    ErrorCode error;

    for( size_t n = 2; n <= 2; ++n )
    {
        Core moab;
        Interface* mbImpl = &moab;
        std::vector< EntityHandle > tris;
        error = create_unitsq_tris( mbImpl, n, tris );MB_CHK_ERR( error );
        EntityHandle meshIn = 0;
        HiReconstruction hirec( dynamic_cast< Core* >( mbImpl ), 0, meshIn );

        for( int degree = 2; degree <= 2; ++degree )
        {
            // reconstruct geometry, interpolation
            hirec.reconstruct3D_surf_geom( degree, true, true, true );
            // test fitting result
            double mxdist = 0;

            for( size_t itri = 0; itri < tris.size(); ++itri )
            {
                const int nvpe                 = 3;
                double naturalcoords2fit[nvpe] = { 1.0 / (double)nvpe, 1.0 / (double)nvpe, 1.0 / (double)nvpe },
                       newcoords[3];
                error = hirec.hiproj_walf_in_element( tris[itri], nvpe, 1, naturalcoords2fit, newcoords );MB_CHK_ERR( error );
                std::vector< EntityHandle > conn;
                error = mbImpl->get_connectivity( &( tris[itri] ), 1, conn );MB_CHK_ERR( error );
                double coords[3 * nvpe], linearcoords[3];
                error = mbImpl->get_coords( &( conn[0] ), nvpe, coords );MB_CHK_ERR( error );
                compute_linear_coords( nvpe, coords, naturalcoords2fit, linearcoords );
                mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
            }

            std::cout << "triangulated unit square n= " << n << " degree= " << degree << " interpolation:\n";
            std::cout << "maximum projection lift is " << mxdist << std::endl;
            // for debug
            // return error;
            mxdist = 0;
            // reconstruct geometry, least square fitting
            hirec.reconstruct3D_surf_geom( degree, false, false, true );

            // test fitting result
            for( size_t itri = 0; itri < tris.size(); ++itri )
            {
                const int nvpe                 = 3;
                double naturalcoords2fit[nvpe] = { 1.0 / (double)nvpe, 1.0 / (double)nvpe, 1.0 / (double)nvpe },
                       newcoords[3];
                error = hirec.hiproj_walf_in_element( tris[itri], nvpe, 1, naturalcoords2fit, newcoords );MB_CHK_ERR( error );
                std::vector< EntityHandle > conn;
                error = mbImpl->get_connectivity( &( tris[itri] ), 1, conn );MB_CHK_ERR( error );
                double coords[3 * nvpe], linearcoords[3];
                error = mbImpl->get_coords( &( conn[0] ), nvpe, coords );MB_CHK_ERR( error );
                compute_linear_coords( nvpe, coords, naturalcoords2fit, linearcoords );
                mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
            }

            std::cout << "unit square n= " << n << " degree= " << degree << " least square:\n";
            std::cout << "maximum projection lift is " << mxdist << std::endl;
        }
    }

    return error;
}
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines