MOAB: Mesh Oriented datABase  (version 5.2.1)
hireconst_test_parallel.cpp
Go to the documentation of this file.
00001 /*This unit test is for high order reconstruction capability, which could be used for mesh
00002  * refinement*/
00003 #include <iostream>
00004 #include <string>
00005 #include <sstream>
00006 #include <sys/time.h>
00007 #include <vector>
00008 #include <algorithm>
00009 #include "moab/Core.hpp"
00010 #include "moab/Range.hpp"
00011 #include "moab/CartVect.hpp"
00012 #include "moab/MeshTopoUtil.hpp"
00013 #include "moab/NestedRefine.hpp"
00014 #include "moab/DiscreteGeometry/HiReconstruction.hpp"
00015 #include "TestUtil.hpp"
00016 #include <math.h>
00017 
00018 #ifdef MOAB_HAVE_MPI
00019 #include "moab/ParallelComm.hpp"
00020 #include "MBParallelConventions.h"
00021 #include "ReadParallel.hpp"
00022 #include "moab/FileOptions.hpp"
00023 #include "MBTagConventions.hpp"
00024 #include "moab_mpi.h"
00025 #endif
00026 
00027 using namespace moab;
00028 
00029 #ifdef MOAB_HAVE_MPI
00030 std::string read_options;
00031 #endif
00032 
00033 #ifdef MOAB_HAVE_HDF5
00034 #undef MOAB_HAVE_HDF5
00035 #endif
00036 ErrorCode load_meshset_hirec( const char* infile, Interface* mbimpl, EntityHandle& meshset, ParallelComm*& pc,
00037                               const int degree = 0, const int dim = 2 );
00038 ErrorCode test_mesh( const char* infile, const int degree, const bool interp, const int dim );
00039 
00040 void compute_linear_coords( const int nvpe, double* elemcoords, double* naturals, double* linearcoords );
00041 
00042 void usage()
00043 {
00044     std::cout << "usage: mpirun -np <number of processors> ./hireconst_test_parallel <mesh file> "
00045                  "-degree <degree> -interp <0=least square, 1=interpolation> -dim <mesh dimension>"
00046               << std::endl;
00047 }
00048 
00049 int main( int argc, char* argv[] )
00050 {
00051 #ifdef MOAB_HAVE_MPI
00052     MPI_Init( &argc, &argv );
00053     int nprocs, rank;
00054     MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
00055     MPI_Comm_rank( MPI_COMM_WORLD, &rank );
00056 #endif
00057     int degree = 3, dim = 2;
00058     bool interp = false;
00059     ErrorCode error;
00060 
00061 #ifdef MOAB_HAVE_HDF5
00062     std::string infile = TestDir + "/mbcslam/fine4.h5m";
00063 #else
00064     std::string infile = TestDir + "/sphere_quads_20.vtk";
00065 #endif
00066 
00067     if( argc == 1 )
00068     {
00069         usage();
00070         std::cout << "Using default arguments: ./hireconst_test_parallel " << infile << " -degree 3 -interp 0 -dim 2"
00071                   << std::endl;
00072     }
00073     else
00074     {
00075         infile      = argv[1];
00076         bool hasdim = false;
00077 
00078         for( int i = 2; i < argc; ++i )
00079         {
00080             if( i + 1 != argc )
00081             {
00082                 if( std::string( argv[i] ) == "-degree" ) { degree = atoi( argv[++i] ); }
00083                 else if( std::string( argv[i] ) == "-interp" )
00084                 {
00085                     interp = atoi( argv[++i] );
00086                 }
00087                 else if( std::string( argv[i] ) == "-dim" )
00088                 {
00089                     dim    = atoi( argv[++i] );
00090                     hasdim = true;
00091                 }
00092                 else
00093                 {
00094 #ifdef MOAB_HAVE_MPI
00095 
00096                     if( 0 == rank ) { usage(); }
00097 
00098                     MPI_Finalize();
00099 #else
00100                     usage();
00101 #endif
00102                     return 0;
00103                 }
00104             }
00105         }
00106 
00107         if( !hasdim )
00108         {
00109 #ifdef MOAB_HAVE_MPI
00110 
00111             if( 0 == rank )
00112             { std::cout << "Dimension of input mesh should be provided, positive and less than 3" << std::endl; }
00113 
00114 #else
00115             std::cout << "Dimension of input mesh should be provided, positive and less than 3" << std::endl;
00116 #endif
00117             return 0;
00118         }
00119 
00120         if( degree <= 0 || dim > 2 || dim <= 0 )
00121         {
00122 #ifdef MOAB_HAVE_MPI
00123 
00124             if( 0 == rank )
00125             {
00126                 std::cout << "Input degree should be positive number;\n";
00127                 std::cout << "Input dimesion should be positive and less than 3;" << std::endl;
00128             }
00129 
00130 #else
00131             std::cout << "Input degree should be positive number;\n";
00132             std::cout << "Input dimesion should be positive and less than 3;" << std::endl;
00133 #endif
00134             return 0;
00135         }
00136 
00137 #ifdef MOAB_HAVE_MPI
00138 
00139         if( 0 == rank )
00140         {
00141             std::cout << "Testing on " << infile << " with dimension " << dim << "\n";
00142             std::string opts = interp ? "interpolation" : "least square fitting";
00143             std::cout << "High order reconstruction with degree " << degree << " " << opts << std::endl;
00144         }
00145 
00146 #else
00147         std::cout << "Testing on " << infile << " with dimension " << dim << "\n";
00148         std::string opts = interp ? "interpolation" : "least square fitting";
00149         std::cout << "High order reconstruction with degree " << degree << " " << opts << std::endl;
00150 #endif
00151     }
00152 
00153     error = test_mesh( infile.c_str(), degree, interp, dim );MB_CHK_ERR( error );
00154 
00155 #ifdef MOAB_HAVE_MPI
00156     MPI_Finalize();
00157 #endif
00158 }
00159 
00160 ErrorCode load_meshset_hirec( const char* infile, Interface* mbimpl, EntityHandle& meshset, ParallelComm*& pc,
00161                               const int degree, const int dim )
00162 {
00163     ErrorCode error;
00164     error = mbimpl->create_meshset( moab::MESHSET_SET, meshset );MB_CHK_ERR( error );
00165 #ifdef MOAB_HAVE_MPI
00166     int nprocs, rank;
00167     MPI_Comm comm = MPI_COMM_WORLD;
00168     MPI_Comm_size( comm, &nprocs );
00169     MPI_Comm_rank( comm, &rank );
00170     EntityHandle partnset;
00171     error = mbimpl->create_meshset( moab::MESHSET_SET, partnset );MB_CHK_ERR( error );
00172 
00173     if( nprocs > 1 ) { pc = moab::ParallelComm::get_pcomm( mbimpl, partnset, &comm ); }
00174 
00175     if( nprocs > 1 )
00176     {
00177         int nghlayers           = degree > 0 ? HiReconstruction::estimate_num_ghost_layers( degree, true ) : 0;
00178         std::string part_method = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;";
00179 #ifndef MOAB_HAVE_HDF5
00180         part_method = "PARALLEL=BCAST_DELETE;PARTITION=TRIVIAL;";
00181 #endif
00182 
00183         if( nghlayers )
00184         {
00185             // get ghost layers
00186             if( dim == 2 ) { read_options = part_method + ";PARALLEL_RESOLVE_SHARED_ENTS;PARALLEL_GHOSTS=2.0."; }
00187             else if( dim == 1 )
00188             {
00189                 read_options = part_method + ";PARALLEL_RESOLVE_SHARED_ENTS;PARALLEL_GHOSTS=1.0.";
00190             }
00191             else
00192             {
00193                 MB_SET_ERR( MB_FAILURE, "unsupported dimension" );
00194             }
00195 
00196             read_options += (char)( '0' + nghlayers );
00197         }
00198         else
00199         {
00200             read_options = part_method + ";PARALLEL_RESOLVE_SHARED_ENTS;";
00201         }
00202 
00203         error = mbimpl->load_file( infile, &meshset, read_options.c_str() );MB_CHK_ERR( error );
00204     }
00205     else
00206     {
00207         error = mbimpl->load_file( infile, &meshset );MB_CHK_ERR( error );
00208     }
00209 
00210 #else
00211     assert( !pc && degree && dim );
00212     error = mbimpl->load_file( infile, &meshset );MB_CHK_ERR( error );
00213 #endif
00214     return error;
00215 }
00216 
00217 ErrorCode test_mesh( const char* infile, const int degree, const bool interp, const int dim )
00218 {
00219     Core moab;
00220     Interface* mbimpl = &moab;
00221     ParallelComm* pc  = NULL;
00222     EntityHandle meshset;
00223 #ifdef MOAB_HAVE_MPI
00224     int nprocs, rank;
00225     MPI_Comm comm = MPI_COMM_WORLD;
00226     MPI_Comm_size( comm, &nprocs );
00227     MPI_Comm_rank( comm, &rank );
00228 #endif
00229 
00230     ErrorCode error;
00231     // mesh will be loaded and communicator pc will be updated
00232     error = load_meshset_hirec( infile, mbimpl, meshset, pc, degree, dim );MB_CHK_ERR( error );
00233     // initialize
00234     HiReconstruction hirec( dynamic_cast< Core* >( mbimpl ), pc, meshset );
00235     Range elems, elems_owned;
00236     error = mbimpl->get_entities_by_dimension( meshset, dim, elems );MB_CHK_ERR( error );
00237     int nelems = elems.size();
00238 
00239 #ifdef MOAB_HAVE_MPI
00240 
00241     if( pc )
00242     {
00243         error = pc->filter_pstatus( elems, PSTATUS_GHOST, PSTATUS_NOT, -1, &elems_owned );MB_CHK_ERR( error );
00244     }
00245     else
00246     {
00247         elems_owned = elems;
00248     }
00249 
00250 #endif
00251 
00252 #ifdef MOAB_HAVE_MPI
00253     std::cout << "Mesh has " << nelems << " elements on Processor " << rank << " in total;";
00254     std::cout << elems_owned.size() << " of which are locally owned elements" << std::endl;
00255 #else
00256     std::cout << "Mesh has " << nelems << " elements" << std::endl;
00257 #endif
00258 
00259     // reconstruction
00260     if( dim == 2 )
00261     {
00262         error = hirec.reconstruct3D_surf_geom( degree, interp, false );MB_CHK_ERR( error );
00263     }
00264     else if( dim == 1 )
00265     {
00266         error = hirec.reconstruct3D_curve_geom( degree, interp, false );MB_CHK_ERR( error );
00267     }
00268 
00269 #ifdef MOAB_HAVE_MPI
00270     std::cout << "HiRec has been done on Processor " << rank << std::endl;
00271 #else
00272     std::cout << "HiRec has been done " << std::endl;
00273 #endif
00274     // fitting
00275     double mxdist = 0;
00276 
00277     for( Range::iterator ielem = elems_owned.begin(); ielem != elems_owned.end(); ++ielem )
00278     {
00279         int nvpe;
00280         const EntityHandle* conn;
00281         error = mbimpl->get_connectivity( *ielem, conn, nvpe );MB_CHK_ERR( error );
00282         double w = 1.0 / (double)nvpe;
00283         std::vector< double > naturalcoords2fit( nvpe, w );
00284         CartVect newcoords, linearcoords;
00285         error = hirec.hiproj_walf_in_element( *ielem, nvpe, 1, &( naturalcoords2fit[0] ), newcoords.array() );
00286 
00287         if( MB_FAILURE == error ) { continue; }
00288 
00289         std::vector< double > coords( 3 * nvpe );
00290         error = mbimpl->get_coords( conn, nvpe, &( coords[0] ) );MB_CHK_ERR( error );
00291         compute_linear_coords( nvpe, &( coords[0] ), &( naturalcoords2fit[0] ), linearcoords.array() );
00292         CartVect nlcoords = newcoords - linearcoords;
00293         mxdist            = std::max( mxdist, nlcoords.length() );
00294         /*#ifdef MOAB_HAVE_MPI
00295             std::cout << "Error on element " << *ielem << " is " << nlcoords.length() << "on
00296         Processor " << rank << std::endl; #else std::cout << "Error on element " << *ielem << " is "
00297         << nlcoords.length() << std::endl; #endif*/
00298     }
00299 
00300 #ifdef MOAB_HAVE_MPI
00301     std::cout << "Maximum projection lift is " << mxdist << " on Processor " << rank << std::endl;
00302 #else
00303     std::cout << "Maximum projection lift is " << mxdist << std::endl;
00304 #endif
00305     return error;
00306 }
00307 
00308 void compute_linear_coords( const int nvpe, double* elemcoords, double* naturals, double* linearcoords )
00309 {
00310     assert( elemcoords && linearcoords );
00311 
00312     for( int i = 0; i < 3; ++i )
00313     {
00314         linearcoords[i] = 0;
00315 
00316         for( int j = 0; j < nvpe; ++j )
00317         {
00318             linearcoords[i] += naturals[j] * elemcoords[3 * j + i];
00319         }
00320     }
00321 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines