ExtractClose.cpp

Go to the documentation of this file.

/**
 * @file ExtractClose.cpp
 * @brief Example demonstrating extraction of mesh elements based on distance from a point
 *
 * This example shows how to:
 * - Load a large mesh file
 * - Extract elements within a specified distance from a point
 * - Handle both Cartesian and spherical coordinate systems
 * - Compute element centroids for distance calculations
 * - Create mesh sets for extracted elements
 * - Write extracted mesh subsets to new files
 *
 * This is particularly useful for debugging large mesh files where
 * you need to examine elements near a specific location without
 * loading the entire mesh into memory.
 *
 * @author MOAB Development Team
 * @date 2024
 *
 
 * \brief This test shows how to extract mesh from a model, based on distance.
 *
 * MOAB's It is needed to extract from large mesh files cells close to some point, where we suspect
 * errors It would be useful for 9Gb input file that we cannot visualize, but we have overlapped
 * elements.
 *
 * @param argc Number of command line arguments
 * @param argv Command line arguments array
 * @return 0 on success, 1 on failure
 */
 
#include <iostream>
#include <cstdlib>
#include <cstdio>
 
#include "moab/Core.hpp"
#include "moab/Interface.hpp"
#include "moab/Range.hpp"
#include "moab/ProgOptions.hpp"
#include "moab/CartVect.hpp"
#include "moab/IntxMesh/IntxUtils.hpp"
 
#ifdef MOAB_HAVE_MPI
#include "moab_mpi.h"
#include "moab/ParallelComm.hpp"
#endif
 
using namespace moab;
using namespace std;
 
#ifdef MOAB_HAVE_HDF5
string test_file_name = string( MESH_DIR ) + string( "/64bricks_512hex_256part.h5m" );
#endif
int main( int argc, char** argv )
{
    ProgOptions opts;
 
    // vertex 25 in filt out  (-0.784768, 0.594952, -0.173698)
    double x = -0.784768, y = 0.594952, z = -0.173698;
    double lat = 28.6551, lon = 61.0204;  // from Charlie Zender's page
    // https://acme-climate.atlassian.net/wiki/spaces/ED/pages/1347092547/Assessing+and+Addressing+Regridding+Weights+in+Map-Files
 
    string inputFile = test_file_name;
    opts.addOpt< string >( "inFile,i", "Specify the input file name string ", &inputFile );
 
    string outFile = "out.h5m";
    opts.addOpt< string >( "outFile,o", "Specify the output file name string ", &outFile );
 
    opts.addOpt< double >( string( "xpos,x" ), string( "x position" ), &x );
    opts.addOpt< double >( string( "ypos,y" ), string( "y position" ), &y );
    opts.addOpt< double >( string( "zpos,z" ), string( "z position" ), &z );
 
    opts.addOpt< double >( string( "latitude,t" ), string( "north latitude in degrees" ), &lat );
    opts.addOpt< double >( string( "longitude,w" ), string( "east longitude in degrees" ), &lon );
 
    bool spherical = false;
    opts.addOpt< void >( "spherical,s", "use spherical coords input (default false) ", &spherical );
 
    double distance = 0.01;
    opts.addOpt< double >( string( "distance,d" ), string( "distance " ), &distance );
 
    opts.parseCommandLine( argc, argv );
 
    int rank         = 0;
    int numProcesses = 1;
    std::string readopts;
 
    // Instantiate
    Core mb;
 
    // Load the file into a new file set
    EntityHandle fileSet;
    MB_CHK_SET_ERR( mb.create_meshset( MESHSET_SET, fileSet ), "Error creating file set" );
    MB_CHK_SET_ERR( mb.load_file( inputFile.c_str(), &fileSet, readopts.c_str() ), "Error loading file" );
 
    if( !rank )
    {
        cout << " reading file " << inputFile << " on " << numProcesses << " task";
        if( numProcesses > 1 ) cout << "s";
        cout << "\n";
    }
    // Get all 2d elements in the file set
    Range elems;
    MB_CHK_SET_ERR( mb.get_entities_by_dimension( fileSet, 2, elems ), "Error getting 2d elements" );
 
    // create a meshset with close elements
    EntityHandle outSet;
 
    // create meshset
    MB_CHK_ERR( mb.create_meshset( MESHSET_SET, outSet ) );
 
    // double sphere radius is 1
    CartVect point( x, y, z );
    if( spherical )
    {
        std::cout << " use spherical coordinates on the sphere of radius 1.\n";
        IntxUtils::SphereCoords pos;
        pos.R   = 1;                 // always on the
        pos.lat = lat * M_PI / 180;  // defined in math.h
        pos.lon = lon * M_PI / 180;
        point   = IntxUtils::spherical_to_cart( pos );
    }
 
    Range closeByCells;
    for( Range::iterator it = elems.begin(); it != elems.end(); it++ )
    {
        EntityHandle cell = *it;
        CartVect center;
        MB_CHK_SET_ERR( mb.get_coords( &cell, 1, &( center[0] ) ), "Can't get cell center coords" );
        double dist = ( center - point ).length();
        if( dist <= distance )
        {
            closeByCells.insert( cell );
        }
    }
 
    MB_CHK_SET_ERR( mb.add_entities( outSet, closeByCells ), "Can't add to entity set" );
 
    int numCells = (int)closeByCells.size();
 
    if( numCells == 0 )
    {
        if( !rank )
            std::cout << " no close cells to the point " << point << " at distance less than " << distance << "\n";
    }
    else
    {
        if( !rank )
            std::cout << " write file " << outFile << " with cells closer than " << distance << " from " << point
                      << "\n";
        string writeOpts;
        if( numProcesses > 1 ) writeOpts = string( "PARALLEL=WRITE_PART;" );
        MB_CHK_SET_ERR( mb.write_file( outFile.c_str(), 0, writeOpts.c_str(), &outSet, 1 ), "Can't write file" );
    }
 
    return 0;
}