BoundaryContinents.cpp File Reference

Example demonstrating creation of continent boundary edge meshes. More...

#include "moab/Core.hpp"
#include "moab/ReadUtilIface.hpp"
#include "moab/CartVect.hpp"
#include <iostream>
#include <fstream>

Include dependency graph for BoundaryContinents.cpp:

Go to the source code of this file.

Functions
double	getLat (CartVect p)
double	getLon (CartVect p)
int	main (int argc, char **argv)

Variables
string	bd_name = string( MESH_DIR ) + string( "/../examples/earth/boundary_points.dat" )
string	loops = string( MESH_DIR ) + string( "/../examples/earth/SaveLoopCounts" )

Detailed Description

Example demonstrating creation of continent boundary edge meshes.

This example shows how to:

• Load continent boundary point data and loop information
• Create edge meshes from boundary point sequences
• Convert 3D spherical coordinates to 2D lat-lon coordinates
• Handle periodic boundary conditions for longitude wrapping
• Remove long edges that cross the date line
• Create mesh sets for boundary visualization
• Write boundary meshes in VTK format

This tool is useful for creating continent boundary meshes for climate and geophysical applications.

Author

MOAB Development Team

Date

2024

Definition in file BoundaryContinents.cpp.

Function Documentation

getLat()

double getLat

(

CartVect

p

)

Examples

BoundaryContinents.cpp.

Definition at line 44 of file BoundaryContinents.cpp.

{
    p.normalize();
    return asin( p[2] );
}

References moab::CartVect::normalize().

Referenced by main().

getLon()

double getLon

(

CartVect

p

)

Examples

BoundaryContinents.cpp.

Definition at line 49 of file BoundaryContinents.cpp.

{
    p.normalize();
    double lon;
 
    lon = atan2( p[1], p[0] );
    if( lon < -2.95 )
    {  // this is Behring Strait :)
        return 2.0 * M_PI + lon;
    }
    else
    {
        return lon;
    }
}

References moab::CartVect::normalize().

Referenced by main().

main()

int main	(	int	argc,
		char **	argv
	)

get loops for boundaries

Examples

BoundaryContinents.cpp.

Definition at line 65 of file BoundaryContinents.cpp.

{
    // Get MOAB instance
    Interface* mb = new( std::nothrow ) Core;
    if( NULL == mb ) return 1;
 
    if( argc > 1 )
    {
        // User has input file for boundary describing continents
        bd_name = argv[1];
    }
    if( argc > 2 )
    {
        // User has input file for loops for continents/islands
        loops = argv[2];
    }
 
    std::vector< double > coords;
    ifstream bdFile( bd_name.c_str() );
    if( !bdFile )
    {
        cout << endl << "Failed to open file " << bd_name << endl;
        ;
        return 1;
    }
    coords.resize( 0 );
    while( !bdFile.eof() )
    {
        double co;
        bdFile >> co;
        coords.push_back( co );
    }
    cout << " number of boundary points:" << coords.size() / 3 << "\n";
    /// get loops for boundaries
    vector< int > loopsindx;
    ifstream loopFile( loops.c_str() );
    while( !loopFile.eof() )
    {
        int indx;
        loopFile >> indx;
        loopsindx.push_back( indx );
    }
    cout << "number of loops: " << loopsindx.size() / 2 << "\n";
    ReadUtilIface* readMeshIface;
    mb->query_interface( readMeshIface );
    Range verts;
 
    MB_CHK_SET_ERR( mb->create_vertices( &coords[0], coords.size() / 3, verts ), "do not create boundary vertices" );
 
    Tag gid = mb->globalId_tag();
 
    int num_edges = 0;
    for( size_t i = 0; i < loopsindx.size() / 2; i++ )
    {
        num_edges += ( loopsindx[2 * i + 1] - loopsindx[2 * i] + 1 );
    }
    EntityHandle handle;
    EntityHandle* conn_array;
    MB_CHK_SET_ERR( readMeshIface->get_element_connect( num_edges, 2, MBEDGE, 1, handle, conn_array ), "do not elems" );
    int i1 = 0;  // index in edge connectivity
    for( size_t i = 0; i < loopsindx.size() / 2; i++ )
    {
        for( int j = loopsindx[2 * i]; j <= loopsindx[2 * i + 1]; j++ )
        {
            int j2;
            j2 = j;
            if( j == loopsindx[2 * i + 1] ) j2 = loopsindx[2 * i] - 1;  // close the loop
            conn_array[2 * i1]     = verts[j - 1];                      // vertices for sure start at 1
            conn_array[2 * i1 + 1] = verts[j2];                         // vertices for sure start at 1
            i1++;
        }
        // the last vertex is first one in the loop
    }
 
    Range bedges( handle, handle + num_edges - 1 );
    EntityHandle bound_set;
    MB_CHK_SET_ERR( mb->create_meshset( MESHSET_SET, bound_set ), "Can't create boundary edges set" );
 
    MB_CHK_SET_ERR( mb->add_entities( bound_set, bedges ), "Can't add edges to boundary set" );
 
    // set global ids for vertices and edges
    vector< int > gids;
    gids.resize( verts.size() );
    for( int j = 0; j < (int)verts.size(); j++ )
    {
        gids[j] = j + 1;
    }
    MB_CHK_SET_ERR( mb->tag_set_data( gid, verts, &gids[0] ), "Can't set global ids on verts" );
    // we have less edges than vertices, we can reuse the array for edges too
    MB_CHK_SET_ERR( mb->tag_set_data( gid, bedges, &gids[0] ), "Can't set global ids on edges" );
 
    mb->write_file( "bound.vtk", 0, 0, &bound_set, 1 );
 
    // now project the mesh in 2d plane
    // get all the vertices coordinates
    std::vector< CartVect > co3;
    co3.resize( verts.size() );
    MB_CHK_SET_ERR( mb->get_coords( verts, &( co3[0][0] ) ), "Can't get vertex coords" );
    for( size_t i = 0; i < verts.size(); i++ )
    {
        CartVect p  = co3[i];
        double lat1 = getLat( p );
        double lon1 = getLon( p );
        co3[i]      = CartVect( lon1, lat1, 0. );
    }
    MB_CHK_SET_ERR( mb->set_coords( verts, &( co3[0][0] ) ), "Can't set new vertex coords" );
    // remove edges in 2d that are too long (longer than 6; they are on the other side..., periodic)
 
    Range longEdges;
    for( Range::iterator eit = bedges.begin(); eit != bedges.end(); ++eit )
    {
        EntityHandle eh = *eit;
        const EntityHandle* conn;
        int nconn;
        MB_CHK_ERR( mb->get_connectivity( eh, conn, nconn ) );
        // get coordinates of nodes
        CartVect c2[2];
        MB_CHK_ERR( mb->get_coords( conn, 2, &( c2[0][0] ) ) );
        double dist = ( c2[1] - c2[0] ).length_squared();
        if( dist > 36. )  // too long edge in 2d, remove it
            longEdges.insert( eh );
    }
    MB_CHK_ERR( mb->delete_entities( longEdges ) );
 
    mb->write_file( "bound2d.vtk" );
    delete mb;
 
    return 0;
}

References moab::Core::add_entities(), bd_name, moab::Range::begin(), moab::Core::create_meshset(), moab::Core::create_vertices(), moab::Core::delete_entities(), moab::Range::end(), moab::Core::get_connectivity(), moab::Core::get_coords(), moab::ReadUtilIface::get_element_connect(), getLat(), getLon(), moab::Core::globalId_tag(), moab::Range::insert(), length_squared(), loops, mb, MB_CHK_ERR, MB_CHK_SET_ERR, MBEDGE, MESHSET_SET, moab::Interface::query_interface(), moab::Core::set_coords(), moab::Range::size(), moab::Core::tag_set_data(), and moab::Core::write_file().

Variable Documentation

bd_name

string bd_name = string( MESH_DIR ) + string( "/../examples/earth/boundary_points.dat" )

Examples

BoundaryContinents.cpp.

Definition at line 41 of file BoundaryContinents.cpp.

Referenced by main().

loops

string loops = string( MESH_DIR ) + string( "/../examples/earth/SaveLoopCounts" )

Examples

BoundaryContinents.cpp.

Definition at line 42 of file BoundaryContinents.cpp.

Referenced by main().