MOAB: Mesh Oriented datABase  (version 5.3.1)
mbex4.cpp
Go to the documentation of this file.
00001 /// \file mbex4.cpp
00002 ///
00004 ///
00005 /// \brief beginner tutorial, example 4: Create a 2D structured mesh
00006 ///        and set some tag data on it
00007 ///
00008 /// In this example, we create a 2D structured mesh (actually a 3D
00009 /// mesh made of quads) and then we will actually set data on the
00010 /// mesh. Tags represent data that is attached to entities. In this
00011 /// example, we will create two tags:
00012 ///
00013 /// -# A "temperature" tag that is a single double precision number
00014 ///    attached to each quad.
00015 /// -# A "velocity" tag that is an array of 2 double precision numbers
00016 ///    attached to each vertex.
00017 ///
00018 /// We will write the mesh out to a file, and you can visualize the
00019 /// data using your favorite tool.
00020 ///
00021 /// In this example, I am demonstrating these operations in the
00022 /// clearest possible way - not using the most efficient method. MOAB
00023 /// has been designed so that users do not need to sacrifice
00024 /// performance - future examples will demonstrate how to
00025 /// access/manipulate the mesh using the fastest methods possible.
00026
00027 // The moab/Core.hpp header file is needed for all MOAB work...
00028 #include "moab/Core.hpp"
00029
00030 // The moab/ScdInterface.hpp contains code which defines the moab
00031 // structured mesh interface
00032 #include "moab/ScdInterface.hpp"
00033
00034 #include <iostream>
00035 #include <cmath>
00036
00037 // ****************
00038 // *              *
00039 // *     main     *
00040 // *              *
00041 // ****************
00042 int main()
00043 {
00044     moab::ErrorCode rval;
00045     moab::Core mbint;
00046
00047     // ***********************
00048     // *   Create the Mesh   *
00049     // ***********************
00050
00051     // First, lets make the mesh. It will be a 100 by 100 uniform grid
00052     // (there will be 100x100 quads, 101x101 vertexes) with dx = dy =
00053     // 0.1. Unlike the previous example, we will first make the mesh,
00054     // then set the coordinates one at a time.
00055
00056     const unsigned NI = 100;
00057     const unsigned NJ = 100;
00058
00059     // moab::ScdInterface is the structured mesh interface class for
00060     // MOAB.
00061     moab::ScdInterface* scdint;
00062
00063     // Tell MOAB that our mesh is structured:
00064     rval = mbint.query_interface( scdint );MB_CHK_SET_ERR( rval, "mbint.query_interface failed" );
00065
00066     // Create the mesh:
00067     moab::ScdBox* scdbox = NULL;
00068     rval = scdint->construct_box( moab::HomCoord( 0, 0, 0 ), moab::HomCoord( NI, NJ, 0 ), NULL, 0, scdbox );MB_CHK_SET_ERR( rval, "scdint->construct_box failed" );
00069
00070     // MOAB knows to make quads instead of hexes because the last start
00071     // and end indexes are the same (0). Note that it is still a "3D"
00072     // mesh because each vertex coordinate is still defined using three
00073     // numbers although every element in the mesh is a quadrilateral.
00074
00075     // ******************************
00076     // *   Set Vertex Coordinates   *
00077     // ******************************
00078
00079     // The "NULL" and "0" arguments in the call to construct_box are
00080     // where we could specify the vertex coordinates. Since we didn't
00081     // give any coordinates, every vertex is given a position of 0,0,0
00082     // by default. Now we will set the vertex coordinates...
00083
00084     const double DX = 0.1;
00085     const double DY = 0.1;
00086
00087     for( unsigned i = 0; i < NI + 1; i++ )
00088         for( unsigned j = 0; j < NJ + 1; j++ )
00089         {
00090             // First, get the entity handle:
00091             moab::EntityHandle handle = scdbox->get_vertex( i, j );
00092
00093             // Compute the coordinate:
00094             double coord[3] = { DX * i, DY * j, 0.0 };
00095
00096             // Change the coordinate of the vertex:
00097             mbint.set_coords( &handle, 1, coord );
00098         }
00099
00100     // *******************
00101     // *   Attach Tags   *
00102     // *******************
00103
00104     // The vertex coordinates have been defined, now let's attach some
00105     // data to the mesh. In MOAB this is done using "tags". Tags are
00106     // little bits of information that can be attached to any mesh
00107     // entity. In our example, we want to create two tags. The
00108     // "temperature" tag will be attached to each quad and will be 1
00109     // double. The "velocity" tag will be attached to each vertex and
00110     // will be an array of two doubles.
00111
00112     // zero and twozeros represent the initial tag
00113
00114     // Create the tags:
00115     moab::Tag temp_tag;
00116     double temp_default_value = 0.0;
00117     rval                      = mbint.tag_get_handle( "temperature", 1, moab::MB_TYPE_DOUBLE, temp_tag,
00118                                  moab::MB_TAG_DENSE | moab::MB_TAG_CREAT, &temp_default_value );MB_CHK_SET_ERR( rval, "mbint.tag_get_handle(temperature) failed" );
00119
00120     moab::Tag vel_tag;
00121     double vel_default_value[2] = { 0.0, 0.0 };
00122     rval = mbint.tag_get_handle( "velocity", 2, moab::MB_TYPE_DOUBLE, vel_tag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT,
00123                                  vel_default_value );MB_CHK_SET_ERR( rval, "mbint.tag_get_handle(velocity) failed" );
00124
00125     // Note that when we created each tag, we specified two flags:
00126     //
00127     // The moab::MB_TAG_DENSE flag tells MOAB that this is a dense
00128     // tag. Dense tags will get automatically assigned to entities which
00129     // have continuous handles. For this example, this means that once
00130     // we set a tag on one vertex, memory will be allocated for
00131     // assigning the tag to all vertexes. The same is true of
00132     // quads. Dense tags are much more efficient when assigning tags to
00133     // lots of entities. If you only want to assign a tag to a few
00134     // entities, it is more efficient to use sparse tags
00135     // (moab::MB_TAG_SPARSE).
00136     //
00137     // The moab::MB_TAG_CREAT flag tells MOAB to create the tag if it
00139
00140     // The tags have now been created, now we have to attach them to
00141     // entities and set their values. NOTE: I am going to do this in a
00142     // manner which emphasizes clarity - this is not the most efficient
00143     // approach - that will be saved for later tutorial examples.
00144
00145     // Loop through each quad and set the temperature:
00146     for( unsigned i = 0; i < NI; i++ )
00147         for( unsigned j = 0; j < NJ; j++ )
00148         {
00149             // Get the handle for this quad:
00150             moab::EntityHandle handle = scdbox->get_element( i, j );
00151
00152             // Compute the temperature...
00153             double xc          = DX * ( i + 0.5 );
00154             double yc          = DY * ( j + 0.5 );
00155             double r           = std::sqrt( xc * xc + yc * yc );
00156             double temperature = std::exp( -0.5 * r );
00157
00158             // Set the temperature on a single quad:
00159             rval = mbint.tag_set_data( temp_tag, &handle, 1, &temperature );MB_CHK_SET_ERR( rval, "mbint.tag_set_data(temp_tag) failed" );
00160         }
00161
00162     // Loop through each vertex and set the velocity:
00163     for( unsigned i = 0; i < NI + 1; i++ )
00164         for( unsigned j = 0; j < NJ + 1; j++ )
00165         {
00166             // Get the handle for this vertex:
00167             moab::EntityHandle handle = scdbox->get_vertex( i, j );
00168             double velocity[2]        = { i, j };
00169
00170             // Set the velocity on a vertex:
00171             rval = mbint.tag_set_data( vel_tag, &handle, 1, velocity );MB_CHK_SET_ERR( rval, "mbint.tag_set_data(vel_tag) failed" );
00172         }
00173
00174     // ***************************
00175     // *   Write Mesh to Files   *
00176     // ***************************
00177
00178     // NOTE: Some visualization software (such as VisIt) may not
00179     // interpret the velocity tag as a vector and you may not be able to
00180     // plot it. But you should be able to plot the temperature on top of
00181     // the mesh.
00182
00183     rval = mbint.write_file( "mbex4.vtk" );MB_CHK_SET_ERR( rval, "write_file(mbex4.vtk) failed" );
00184
00185     return 0;
00186 }