moab
structuredmesh.cpp
Go to the documentation of this file.
00001 /* \example structuredmesh structuredmesh.cpp
00002  * \brief Show creation and query of structured mesh through MOAB's structured mesh interface.
00003  * This is a serial example showing creation and query of a 3D structured mesh.  A single
00004  * rectangular brick of hex elements is created, then referenced by its ijk parameterization.
00005  * 1D and 2D examples could be made simply by changing the dimension parameter and the number
00006  * of ijk parameters passed into the MOAB functions.
00007  *
00008  * This example:
00009  * 0. Instantiate MOAB and get the structured mesh interface
00010  * 1. Creates a IxJxK structured mesh, which includes I*J*K vertices and (I-1)*(J-1)*(K-1) hexes.
00011  * 2. Get the vertices and hexes from moab and check their numbers against I*J*K and (I-1)*(J-1)*(K-1), resp.
00012  * 3. Loop over elements in 3 nested loops over i, j, k; for each (i,j,k):
00013  * 3a. Get the element corresponding to (i,j,k)
00014  * 3b. Get the connectivity of the element
00015  * 3c. Get the coordinates of the vertices comprising that element
00016  * 4. Release the structured mesh interface and destroy the MOAB instance
00017  *
00018  * To run: ./structuredmesh <dim> <N>
00019  * (default values so can run w/ no user interaction)
00020  */
00021 
00022 #include "moab/Core.hpp"
00023 #include "moab/ScdInterface.hpp"
00024 #include <iostream>
00025 #include <vector>
00026 
00027 using namespace moab;
00028 
00029 int main(int argc, char **argv) 
00030 {
00031   argv[0] = argv[argc - argc]; // To remove the warnings about unused parameters
00032   int I, J, K;
00033     // progoptions?
00034   std::cout << "Enter I, J, K... " << std::endl;
00035   std::cin >> I >> J >> K;
00036   
00037     // 0. Instantiate MOAB and get the structured mesh interface
00038   Interface *mb = new Core();
00039   ScdInterface *scdiface;
00040   ErrorCode rval = mb->query_interface(scdiface); // get a ScdInterface object through moab instance
00041   if (MB_SUCCESS != rval) return rval;
00042   
00043     // 1. Creates a IxJxK structured mesh, which includes I*J*K vertices and (I-1)*(J-1)*(K-1) hexes.
00044   ScdBox *box;
00045   rval = scdiface->construct_box(HomCoord(0, 0, 0), HomCoord(I-1, J-1, K-1), // low, high box corners in parametric space
00046                                  NULL, 0, // NULL coords vector and 0 coords (don't specify coords for now)
00047                                  box);    // box is the structured box object providing the parametric
00048                                           // structured mesh interface for this rectangle of elements
00049   if (MB_SUCCESS != rval) return rval;
00050 
00051     // 2. Get the vertices and hexes from moab and check their numbers against I*J*K and (I-1)*(J-1)*(K-1), resp.
00052   Range verts, hexes;
00053   rval = mb->get_entities_by_dimension(0, 0, verts); // first '0' specifies "root set", or entire MOAB instance, second the entity dimension being requested
00054   if (MB_SUCCESS != rval) return rval;
00055   rval = mb->get_entities_by_dimension(0, 3, hexes);
00056   if (MB_SUCCESS != rval) return rval;
00057 
00058   if ((I-1)*(J-1)*(K-1) == (int) hexes.size() && I*J*K == (int) verts.size())
00059     std::cout << "Created " << hexes.size() << " hexes and " << verts.size() << " vertices." << std::endl;
00060   else
00061     std::cout << "Created the wrong number of vertices or hexes!" << std::endl;
00062   
00063     // 3. Loop over elements in 3 nested loops over i, j, k; for each (i,j,k):
00064   std::vector<double> coords(8*3);
00065   std::vector<EntityHandle> connect;
00066   for (int k = 0; k < K-1; k++) {
00067     for (int j = 0; j < J-1; j++) {
00068       for (int i = 0; i < I-1; i++) {
00069           // 3a. Get the element corresponding to (i,j,k)
00070         EntityHandle ehandle = box->get_element(i, j, k);
00071         if (0 == ehandle) return MB_FAILURE;
00072           // 3b. Get the connectivity of the element
00073         rval = mb->get_connectivity(&ehandle, 1, connect); // get the connectivity, in canonical order
00074         if (MB_SUCCESS != rval) return rval;
00075           // 3c. Get the coordinates of the vertices comprising that element
00076         rval = mb->get_coords(connect.data(), connect.size(), coords.data()); // get the coordinates of those vertices
00077         if (MB_SUCCESS != rval) return rval;
00078       }
00079     }
00080   }
00081 
00082     // 4. Release the structured mesh interface and destroy the MOAB instance
00083   mb->release_interface(scdiface); // tell MOAB we're done with the ScdInterface
00084   delete mb;
00085   
00086   return 0;
00087 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines