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