spherical_area_test.cpp

Go to the documentation of this file.

/*
 * spherical_area_test.cpp
 *
 *  Created on: Feb 1, 2013
 */
#include <iostream>
#include "moab/Core.hpp"
#include "moab/Interface.hpp"
#include "moab/IntxMesh/IntxUtils.hpp"
#include "TestUtil.hpp"

using namespace moab;

int main( int /* argc*/, char** /* argv[]*/ )
{
    // check command line arg
    const char* filename_mesh = STRINGIFY( MESHDIR ) "/mbcslam/eulerHomme.vtk";

    // read input mesh in a set
    Core moab;
    Interface* mb = &moab;  // global
    EntityHandle sf;
    ErrorCode rval = mb->create_meshset( MESHSET_SET, sf );
    if( MB_SUCCESS != rval ) return 1;

    rval = mb->load_file( filename_mesh, &sf );
    if( MB_SUCCESS != rval ) return 1;

    double R = 6.;  // should be input
    // compare total area with 4*M_PI * R^2

    const double area_sphere = R * R * M_PI * 4.;
    std::cout << "total area of the sphere        :  " << area_sphere << "\n";

    {
        moab::IntxAreaUtils areaAdaptor( moab::IntxAreaUtils::Girard );  // use_lHuiller = true
        double area1 = areaAdaptor.area_on_sphere( mb, sf, R );
        std::cout << "total area with Girard          :  " << area1
                  << " rel error:" << fabs( ( area1 - area_sphere ) / area_sphere ) << "\n";
    }

    {
        moab::IntxAreaUtils areaAdaptor( moab::IntxAreaUtils::lHuiller );
        double area2 = areaAdaptor.area_on_sphere( mb, sf, R );
        std::cout << "total area with l'Huiller       : " << area2
                  << " rel error:" << fabs( ( area2 - area_sphere ) / area_sphere ) << "\n";
    }

#ifdef MOAB_HAVE_TEMPESTREMAP
    {
        moab::IntxAreaUtils areaAdaptor( moab::IntxAreaUtils::GaussQuadrature );
        double area3 = areaAdaptor.area_on_sphere( mb, sf, R );
        std::cout << "total area with GaussQuadrature : " << area3
                  << " rel error:" << fabs( ( area3 - area_sphere ) / area_sphere ) << "\n";
    }
#endif

    return 0;
}