MOAB: Mesh Oriented datABase  (version 5.2.1)
case1_test.cpp
Go to the documentation of this file.
00001 /*
00002  * case1_test.cpp
00003  *
00004  *  Created on: Feb 12, 2013
00005  */
00006 
00007 // copy from par_sph_intx
00008 // will save the LOC tag on the euler nodes? maybe
00009 #include <iostream>
00010 #include <cmath>  // for M_PI
00011 
00012 #include "moab/Core.hpp"
00013 #include "moab/Interface.hpp"
00014 #include "moab/IntxMesh/Intx2MeshOnSphere.hpp"
00015 #include "moab/ProgOptions.hpp"
00016 #include "MBTagConventions.hpp"
00017 
00018 #include "moab/IntxMesh/IntxUtils.hpp"
00019 #include "IntxUtilsCSLAM.hpp"
00020 
00021 #include "TestUtil.hpp"
00022 
00023 using namespace moab;
00024 // some input data
00025 double gtol     = 0.0001;        // this is for geometry tolerance
00026 double CubeSide = 6.;            // the above file starts with cube side 6; radius depends on cube side
00027 double t = 0.1, delta_t = 0.43;  // check the script
00028 
00029 ErrorCode manufacture_lagrange_mesh_on_sphere( Interface* mb, EntityHandle euler_set, EntityHandle& lagr_set )
00030 {
00031     /*
00032      * get all quads first, then vertices, then move them on the surface of the sphere
00033      *  radius is in, it comes from MeshKit/python/examples/manufHomme.py :
00034      *  length = 6.
00035      *  each edge of the cube will be divided using this meshcount
00036      *  meshcount = 11
00037      *   circumscribed sphere radius
00038      *   radius = length * math.sqrt(3) /2
00039      */
00040     double radius = CubeSide / 2 * sqrt( 3. );  // our value depends on cube side
00041     Range quads;
00042     ErrorCode rval = mb->get_entities_by_type( euler_set, MBQUAD, quads );
00043     if( MB_SUCCESS != rval ) return rval;
00044 
00045     Range connecVerts;
00046     rval = mb->get_connectivity( quads, connecVerts );
00047     if( MB_SUCCESS != rval ) return rval;
00048 
00049     // create new set
00050     rval = mb->create_meshset( MESHSET_SET, lagr_set );
00051     if( MB_SUCCESS != rval ) return rval;
00052 
00053     // get the coordinates of the old mesh, and move it around the sphere according to case 1
00054     // now put the vertices in the right place....
00055     // int vix=0; // vertex index in new array
00056 
00057     // first create departure points (vertices in the lagrange mesh)
00058     // then connect them in quads
00059     std::map< EntityHandle, EntityHandle > newNodes;
00060     for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
00061     {
00062         EntityHandle oldV = *vit;
00063         CartVect posi;
00064         rval = mb->get_coords( &oldV, 1, &( posi[0] ) );
00065         if( MB_SUCCESS != rval ) return rval;
00066         // Intx utils, case 1
00067         CartVect newPos;
00068         IntxUtilsCSLAM::departure_point_case1( posi, t, delta_t, newPos );
00069         newPos = radius * newPos;
00070         EntityHandle new_vert;
00071         rval = mb->create_vertex( &( newPos[0] ), new_vert );
00072         if( MB_SUCCESS != rval ) return rval;
00073         newNodes[oldV] = new_vert;
00074     }
00075     for( Range::iterator it = quads.begin(); it != quads.end(); ++it )
00076     {
00077         EntityHandle q = *it;
00078         int nnodes;
00079         const EntityHandle* conn4;
00080         rval = mb->get_connectivity( q, conn4, nnodes );
00081         if( MB_SUCCESS != rval ) return rval;
00082         EntityHandle new_conn[4];
00083         for( int i = 0; i < nnodes; i++ )
00084         {
00085             EntityHandle v1 = conn4[i];
00086             new_conn[i]     = newNodes[v1];
00087         }
00088         EntityHandle new_quad;
00089         rval = mb->create_element( MBQUAD, new_conn, 4, new_quad );
00090         if( MB_SUCCESS != rval ) return rval;
00091         rval = mb->add_entities( lagr_set, &new_quad, 1 );
00092         if( MB_SUCCESS != rval ) return rval;
00093     }
00094 
00095     return rval;
00096 }
00097 int main( int argc, char** argv )
00098 {
00099 
00100     const char* filename_mesh1 = STRINGIFY( MESHDIR ) "/mbcslam/eulerHomme.vtk";
00101     if( argc > 1 )
00102     {
00103         int index = 1;
00104         while( index < argc )
00105         {
00106             if( !strcmp( argv[index], "-gtol" ) )  // this is for geometry tolerance
00107             { gtol = atof( argv[++index] ); }
00108             if( !strcmp( argv[index], "-cube" ) ) { CubeSide = atof( argv[++index] ); }
00109             if( !strcmp( argv[index], "-dt" ) ) { delta_t = atof( argv[++index] ); }
00110             if( !strcmp( argv[index], "-input" ) ) { filename_mesh1 = argv[++index]; }
00111             index++;
00112         }
00113     }
00114     std::cout << " case 1: use -gtol " << gtol << " -dt " << delta_t << " -cube " << CubeSide << " -input "
00115               << filename_mesh1 << "\n";
00116 
00117     Core moab;
00118     Interface& mb = moab;
00119     EntityHandle euler_set;
00120     ErrorCode rval;
00121     rval = mb.create_meshset( MESHSET_SET, euler_set );
00122     if( MB_SUCCESS != rval ) return 1;
00123 
00124     rval = mb.load_file( filename_mesh1, &euler_set );
00125 
00126     if( MB_SUCCESS != rval ) return 1;
00127 
00128     // everybody will get a DP tag, including the non owned entities; so exchange tags is not
00129     // required for LOC (here)
00130     EntityHandle lagrange_set;
00131     rval = manufacture_lagrange_mesh_on_sphere( &mb, euler_set, lagrange_set );
00132     if( MB_SUCCESS != rval ) return 1;
00133     rval = mb.write_file( "lagrIni.h5m", 0, 0, &lagrange_set, 1 );
00134     if( MB_SUCCESS != rval ) std::cout << "can't write lagr set\n";
00135 
00136     rval = moab::IntxUtils::enforce_convexity( &mb, lagrange_set );
00137     if( MB_SUCCESS != rval ) return 1;
00138 
00139     rval = mb.write_file( "lagr.h5m", 0, 0, &lagrange_set, 1 );
00140     if( MB_SUCCESS != rval ) std::cout << "can't write lagr set\n";
00141 
00142     Intx2MeshOnSphere worker( &mb );
00143 
00144     double radius = CubeSide / 2 * sqrt( 3. );  // input
00145 
00146     worker.set_radius_source_mesh( radius );
00147     worker.set_radius_destination_mesh( radius );
00148 
00149     // worker.SetEntityType(MBQUAD);
00150 
00151     worker.set_error_tolerance( gtol );
00152     std::cout << "error tolerance epsilon_1=" << gtol << "\n";
00153 
00154     EntityHandle outputSet;
00155     rval = mb.create_meshset( MESHSET_SET, outputSet );
00156     if( MB_SUCCESS != rval ) return 1;
00157     rval = worker.FindMaxEdges( lagrange_set, euler_set );
00158     if( MB_SUCCESS != rval ) return 1;
00159     rval = worker.intersect_meshes( lagrange_set, euler_set, outputSet );
00160     if( MB_SUCCESS != rval ) return 1;
00161 
00162     // std::string opts_write("");
00163     std::stringstream outf;
00164     outf << "intersect1"
00165          << ".h5m";
00166     rval = mb.write_file( outf.str().c_str(), 0, 0, &outputSet, 1 );
00167     if( MB_SUCCESS != rval ) std::cout << "can't write output\n";
00168 
00169     moab::IntxAreaUtils sphAreaUtils;
00170     double intx_area    = sphAreaUtils.area_on_sphere( &mb, outputSet, radius );
00171     double arrival_area = sphAreaUtils.area_on_sphere( &mb, euler_set, radius );
00172     std::cout << " Arrival area: " << arrival_area << "  intersection area:" << intx_area
00173               << " rel error: " << fabs( ( intx_area - arrival_area ) / arrival_area ) << "\n";
00174     if( MB_SUCCESS != rval ) return 1;
00175 
00176     return 0;
00177 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines