mbex2.cpp

Go to the documentation of this file.

/// \file mbex2.cpp
///
/// \author Milad Fatenejad
///
/// \brief beginner tutorial, example 2: Demonstrates loading a mesh
/// from a file, finding coordinate locations, connectivity
/// information...
///
/// In this example, we read in the VTK file (mbex1.vtk) generated
/// in example 1, pick one of the hexahedrons, find out the vertexes
/// that define it (connectivity), and find the coordinates of those
/// vertexes. Finally, we will rotate the mesh a little and write it
/// out to a new file.
 
// The moab/Core.hpp header file is needed for all MOAB work...
#include "moab/Core.hpp"
 
#include <iostream>
#include <iomanip>
#include <cmath>
 
// ****************
// * *
// * main *
// * *
// ****************
int main()
{
 moab::ErrorCode rval;
 moab::Core mbcore;
 moab::Interface& mbint = mbcore;
 
 // First, lets load the file from the previous example. Note that
 // you must compile and run the previous example to get mbex1.vtk!
 rval = mbint.load_file( "mbex1.vtk" );MB_CHK_SET_ERR( rval, "load_file failed" );
 
 // *********************
 // * Introspection *
 // *********************
 
 // We can now access everything about the mesh through mbint. Lets
 // pick an element and find some information about it. In this case,
 // we know that all of the elements in the mesh are hexahedrons, so
 // we will ask MOAB for a range containing the handle for all hexes:
 moab::Range hex_range;
 rval = mbint.get_entities_by_type( 0, moab::MBHEX, hex_range );MB_CHK_SET_ERR( rval, "get_entities_by_type(HEX) failed" );
 
 std::cout << "Loaded a mesh containing: " << hex_range << std::endl;
 
 // Let's analyze one of the hexes (in this case, I picked hex three):
 moab::EntityHandle handle = hex_range[3];
 
 // Find out the eight vertexes that define this particular hex:
 moab::Range connectivity;
 rval = mbint.get_connectivity( &handle, 1, connectivity );MB_CHK_SET_ERR( rval, "get_connectivity failed" );
 
 // Connectivity should now contain the handles for the eight
 // vertexes of interest. Lets print the handle and coordinate of
 // each vertex. Note how we iterate through the range just like with
 // STL containers...
 double coord[3];
 moab::Range::iterator iter;
 
 std::cout << std::setw( 6 ) << "Handle" << std::setw( 10 ) << "X" << std::setw( 10 ) << "Y" << std::setw( 10 )
 << "Z" << std::endl;
 
 for( iter = connectivity.begin(); iter != connectivity.end(); ++iter )
 {
 rval = mbint.get_coords( &( *iter ), 1, coord );MB_CHK_SET_ERR( rval, "get_coords" );
 
 // Print the entity handle followed by the x, y, z coordinate of
 // the vertex:
 std::cout << std::setw( 6 ) << *iter << std::setw( 10 ) << coord[0] << std::setw( 10 ) << coord[1]
 << std::setw( 10 ) << coord[2] << std::endl;
 }
 
 // ***********************
 // * Rotate the Mesh *
 // ***********************
 
 // Now, let's rotate the mesh about the z-axis by pi/4 radians. Do
 // to this, we will iterate through all of the vertexes
 
 // Start by getting a range containing all of the vertex handles:
 moab::Range vertex_range;
 rval = mbint.get_entities_by_type( 0, moab::MBVERTEX, vertex_range );MB_CHK_SET_ERR( rval, "get_entities_by_type(VERTEX) failed" );
 
 // Get the coordinates of all of the vertexes:
 std::vector< double > vertex_coords( 3 * vertex_range.size() );
 rval = mbint.get_coords( vertex_range, vertex_coords.data() );MB_CHK_SET_ERR( rval, "get_coords" );
 
 unsigned count = 0;
 const double PI = 3.14159265359;
 const double ANGLE = PI / 4;
 for( iter = vertex_range.begin(); iter != vertex_range.end(); ++iter )
 {
 
 // Save the old coordinates:
 double x = vertex_coords[count + 0];
 double y = vertex_coords[count + 1];
 // double z = vertex_coords[count+2];
 
 // Apply the rotation:
 vertex_coords[count + 0] = x * std::cos( ANGLE ) - y * std::sin( ANGLE );
 vertex_coords[count + 1] = x * std::sin( ANGLE ) + y * std::cos( ANGLE );
 
 count += 3;
 }
 
 // Now, the vertex_coords vector contains all of the updated
 // coordinates. Let's push these back into MOAB:
 mbint.set_coords( vertex_range, vertex_coords.data() );MB_CHK_SET_ERR( rval, "set_coords" );
 
 // **************************
 // * Write Mesh to File *
 // **************************
 
 rval = mbint.write_file( "mbex2.vtk" );MB_CHK_SET_ERR( rval, "write_file(mbex2.vtk) failed" );
 
 // Now you can open your favorite visualization tool (VisIt or ParaView)
 // and look at the original (mbex1.vtk) and rotated (mbex2.vtk) meshes.
 
 return 0;
}