ExtractClose.cpp File Reference

Example demonstrating extraction of mesh elements based on distance from a point. More...

#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"

Include dependency graph for ExtractClose.cpp:

Go to the source code of this file.

Functions
int	main (int argc, char **argv)

Detailed Description

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

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.

Parameters

argc	Number of command line arguments
argv	Command line arguments array

Returns

0 on success, 1 on failure

Definition in file ExtractClose.cpp.

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 54 of file ExtractClose.cpp.

{
    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;
}

References moab::Core::add_entities(), ProgOptions::addOpt(), moab::Range::begin(), center(), moab::Core::create_meshset(), moab::Range::end(), moab::Core::get_coords(), moab::Core::get_entities_by_dimension(), moab::Range::insert(), moab::IntxUtils::SphereCoords::lat, length(), moab::Core::load_file(), moab::IntxUtils::SphereCoords::lon, mb, MB_CHK_ERR, MB_CHK_SET_ERR, MESHSET_SET, ProgOptions::parseCommandLine(), moab::IntxUtils::SphereCoords::R, moab::Range::size(), moab::IntxUtils::spherical_to_cart(), test_file_name, and moab::Core::write_file().