MOAB: Mesh Oriented datABase  (version 5.4.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             {
00108                 gtol = atof( argv[++index] );
00109             }
00110             if( !strcmp( argv[index], "-cube" ) )
00111             {
00112                 CubeSide = atof( argv[++index] );
00113             }
00114             if( !strcmp( argv[index], "-dt" ) )
00115             {
00116                 delta_t = atof( argv[++index] );
00117             }
00118             if( !strcmp( argv[index], "-input" ) )
00119             {
00120                 filename_mesh1 = argv[++index];
00121             }
00122             index++;
00123         }
00124     }
00125     std::cout << " case 1: use -gtol " << gtol << " -dt " << delta_t << " -cube " << CubeSide << " -input "
00126               << filename_mesh1 << "\n";
00127 
00128     Core moab;
00129     Interface& mb = moab;
00130     EntityHandle euler_set;
00131     ErrorCode rval;
00132     rval = mb.create_meshset( MESHSET_SET, euler_set );
00133     if( MB_SUCCESS != rval ) return 1;
00134 
00135     rval = mb.load_file( filename_mesh1, &euler_set );
00136 
00137     if( MB_SUCCESS != rval ) return 1;
00138 
00139     // everybody will get a DP tag, including the non owned entities; so exchange tags is not
00140     // required for LOC (here)
00141     EntityHandle lagrange_set;
00142     rval = manufacture_lagrange_mesh_on_sphere( &mb, euler_set, lagrange_set );
00143     if( MB_SUCCESS != rval ) return 1;
00144     rval = mb.write_file( "lagrIni.h5m", 0, 0, &lagrange_set, 1 );
00145     if( MB_SUCCESS != rval ) std::cout << "can't write lagr set\n";
00146 
00147     rval = moab::IntxUtils::enforce_convexity( &mb, lagrange_set );
00148     if( MB_SUCCESS != rval ) return 1;
00149 
00150     rval = mb.write_file( "lagr.h5m", 0, 0, &lagrange_set, 1 );
00151     if( MB_SUCCESS != rval ) std::cout << "can't write lagr set\n";
00152 
00153     Intx2MeshOnSphere worker( &mb );
00154 
00155     double radius = CubeSide / 2 * sqrt( 3. );  // input
00156 
00157     worker.set_radius_source_mesh( radius );
00158     worker.set_radius_destination_mesh( radius );
00159 
00160     // worker.SetEntityType(MBQUAD);
00161 
00162     worker.set_error_tolerance( gtol );
00163     std::cout << "error tolerance epsilon_1=" << gtol << "\n";
00164 
00165     EntityHandle outputSet;
00166     rval = mb.create_meshset( MESHSET_SET, outputSet );
00167     if( MB_SUCCESS != rval ) return 1;
00168     rval = worker.FindMaxEdges( lagrange_set, euler_set );
00169     if( MB_SUCCESS != rval ) return 1;
00170     rval = worker.intersect_meshes( lagrange_set, euler_set, outputSet );
00171     if( MB_SUCCESS != rval ) return 1;
00172 
00173     // std::string opts_write("");
00174     std::stringstream outf;
00175     outf << "intersect1"
00176          << ".h5m";
00177     rval = mb.write_file( outf.str().c_str(), 0, 0, &outputSet, 1 );
00178     if( MB_SUCCESS != rval ) std::cout << "can't write output\n";
00179 
00180     moab::IntxAreaUtils sphAreaUtils;
00181     double intx_area    = sphAreaUtils.area_on_sphere( &mb, outputSet, radius );
00182     double arrival_area = sphAreaUtils.area_on_sphere( &mb, euler_set, radius );
00183     std::cout << " Arrival area: " << arrival_area << "  intersection area:" << intx_area
00184               << " rel error: " << fabs( ( intx_area - arrival_area ) / arrival_area ) << "\n";
00185     if( MB_SUCCESS != rval ) return 1;
00186 
00187     return 0;
00188 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines