MOAB: Mesh Oriented datABase  (version 5.2.1)
addfield.cpp
Go to the documentation of this file.
00001 #include <iostream>
00002 #include "moab/Core.hpp"
00003 #include <math.h>
00004 #include <cstdlib>
00005 #include <assert.h>
00006 
00007 using namespace std;
00008 
00009 namespace moab
00010 {
00011 
00012 // make a nice picture. tweak here
00013 // 217: x,y -> [-8, 8] /// z -> [-65, 65]
00014 double physField( double x, double y, double z, double factor )
00015 {
00016 
00017     // 1/r^2 decay from {0,0,0}
00018     // tjt - changing to 1/r
00019     /*
00020   double scale = 1.0/factor;
00021   out = fabs(x*scale) + fabs(y*scale) + fabs(z*scale);
00022   out += 1e-1; // clamp
00023   out = 1/out;
00024     */
00025     double out = factor * sqrt( x * x + y * y + z * z );
00026 
00027     return out;
00028 }
00029 
00030 void putElementField( Interface* mbi, const char* tagname, double factor )
00031 {
00032     Range elems;
00033 
00034     mbi->get_entities_by_dimension( 0, 3, elems );
00035 
00036     const double defVal = 0.;
00037     Tag fieldTag;
00038     ErrorCode rval = mbi->tag_get_handle( tagname, 1, MB_TYPE_DOUBLE, fieldTag, MB_TAG_DENSE | MB_TAG_CREAT, &defVal );
00039     assert( MB_SUCCESS == rval );
00040 #ifdef NDEBUG
00041     if( MB_SUCCESS == rval ) {};  // Line to avoid compiler warning about unused variable
00042 #endif
00043 
00044     int numElems = elems.size();
00045 
00046     for( int i = 0; i < numElems; i++ )
00047     {
00048         // cout << elems[i] << endl;
00049         EntityHandle elem = elems[i];
00050 
00051         double xyz[3];
00052         mbi->get_coords( &elem, 1, xyz );
00053 
00054         double fieldValue = physField( xyz[0], xyz[1], xyz[2], factor );
00055 
00056         mbi->tag_set_data( fieldTag, &elem, 1, &fieldValue );
00057     }
00058 }
00059 
00060 void putSpectralElementField( Interface* mbi, int dim, int np, const char* tagname, double factor )
00061 {
00062     Range elems;
00063 
00064     mbi->get_entities_by_dimension( 0, dim, elems );
00065 
00066     const int ndofperE = np * np;
00067     double* defDouble  = new double[ndofperE];
00068     Tag fieldTag;
00069     ErrorCode rval =
00070         mbi->tag_get_handle( tagname, ndofperE, MB_TYPE_DOUBLE, fieldTag, MB_TAG_DENSE | MB_TAG_CREAT, &defDouble );
00071     std::cout << "rval = " << rval << std::endl;
00072     assert( MB_SUCCESS == rval || rval == MB_ALREADY_ALLOCATED );
00073 #ifdef NDEBUG
00074     if( MB_SUCCESS == rval || rval == MB_ALREADY_ALLOCATED ) {
00075     };  // Line to avoid compiler warning about unused variable
00076 #endif
00077 
00078     int numElems = elems.size();
00079     std::vector< double > tData( ndofperE, 0.0 );
00080 
00081     for( int i = 0; i < numElems; i++ )
00082     {
00083         EntityHandle elem = elems[i];
00084 
00085         double xyz[3];
00086         mbi->get_coords( &elem, 1, xyz );
00087 
00088         double fieldValue = physField( xyz[0], xyz[1], xyz[2], factor );
00089         for( int j = 0; j < ndofperE; ++j )
00090             tData[j] = fieldValue;  // same value evaluated at the center; could modify it to use
00091                                     // GLL points ?
00092 
00093         mbi->tag_set_data( fieldTag, &elem, 1, &tData[0] );
00094     }
00095 
00096     delete[] defDouble;
00097     tData.clear();
00098 }
00099 
00100 void putVertexField( Interface* mbi, const char* tagname, double factor )
00101 {
00102     Range verts;
00103 
00104     mbi->get_entities_by_type( 0, MBVERTEX, verts );
00105 
00106     const double defVal = 0.;
00107     Tag fieldTag;
00108     mbi->tag_get_handle( tagname, 1, MB_TYPE_DOUBLE, fieldTag, MB_TAG_DENSE | MB_TAG_CREAT, &defVal );
00109 
00110     int numVerts = verts.size();
00111     for( int i = 0; i < numVerts; i++ )
00112     {
00113         EntityHandle vert = verts[i];  //?
00114 
00115         double vertPos[3];
00116         mbi->get_coords( &vert, 1, vertPos );
00117 
00118         double fieldValue = physField( vertPos[0], vertPos[1], vertPos[2], factor );
00119 
00120         mbi->tag_set_data( fieldTag, &vert, 1, &fieldValue );
00121     }
00122 }
00123 
00124 }  // namespace moab
00125 
00126 // Using moab instead of imesh for dev speed
00127 int main( int argc, char** argv )
00128 {
00129 
00130     using namespace moab;
00131 
00132     Interface* mbi = new Core();
00133 
00134     if( argc < 3 )
00135     {
00136         cout << "Usage: " << argv[0] << " <infile> <outfile> [factor]\n"
00137              << "Writes both vertex and element fields.\n";
00138         return 0;
00139     }
00140 
00141     mbi->load_mesh( argv[1] );
00142 
00143     double factor = 1.0;
00144     if( argc == 4 ) factor = atof( argv[3] );
00145 
00146     putVertexField( mbi, "vertex_field", factor );
00147     putElementField( mbi, "element_field", factor );
00148     // putSpectralElementField(mbi, 2, 4, "spectral_element_field", factor);
00149     putSpectralElementField( mbi, 2, 4, "a2oTAG", factor );
00150 
00151     ErrorCode result = mbi->write_mesh( argv[2] );
00152     if( MB_SUCCESS == result )
00153         cout << "wrote " << argv[2] << endl;
00154     else
00155         cout << "Failed to write " << argv[2] << endl;
00156 
00157     //  vector<double> coords;
00158     //  mbi->get_vertex_coordinates(coords);
00159     //  double xavg = 0;
00160     //  for (int i = 0; i < coords.size()/3; i++) xavg += coords[i];
00161     //  cout << xavg << endl;
00162 
00163     return 1;
00164 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines