addfield.cpp

Go to the documentation of this file.

#include <iostream>
#include "moab/Core.hpp"
#include <cmath>
#include <cstdlib>
#include <cassert>
 
using namespace std;
 
namespace moab
{
 
// make a nice picture. tweak here
// 217: x,y -> [-8, 8] /// z -> [-65, 65]
double physField( double x, double y, double z, double factor )
{
 
 // 1/r^2 decay from {0,0,0}
 // tjt - changing to 1/r
 /*
 double scale = 1.0/factor;
 out = fabs(x*scale) + fabs(y*scale) + fabs(z*scale);
 out += 1e-1; // clamp
 out = 1/out;
 */
 double out = factor * sqrt( x * x + y * y + z * z );
 
 return out;
}
 
void putElementField( Interface* mbi, const char* tagname, double factor )
{
 Range elems;
 
 mbi->get_entities_by_dimension( 0, 3, elems );
 
 const double defVal = 0.;
 Tag fieldTag;
 ErrorCode rval = mbi->tag_get_handle( tagname, 1, MB_TYPE_DOUBLE, fieldTag, MB_TAG_DENSE | MB_TAG_CREAT, &defVal );
 assert( MB_SUCCESS == rval );
#ifdef NDEBUG
 if( MB_SUCCESS == rval )
 {
 }; // Line to avoid compiler warning about unused variable
#endif
 
 int numElems = elems.size();
 
 for( int i = 0; i < numElems; i++ )
 {
 // cout << elems[i] << endl;
 EntityHandle elem = elems[i];
 
 double xyz[3];
 mbi->get_coords( &elem, 1, xyz );
 
 double fieldValue = physField( xyz[0], xyz[1], xyz[2], factor );
 
 mbi->tag_set_data( fieldTag, &elem, 1, &fieldValue );
 }
}
 
void putSpectralElementField( Interface* mbi, int dim, int np, const char* tagname, double factor )
{
 Range elems;
 
 mbi->get_entities_by_dimension( 0, dim, elems );
 
 const int ndofperE = np * np;
 double* defDouble = new double[ndofperE];
 Tag fieldTag;
 ErrorCode rval =
 mbi->tag_get_handle( tagname, ndofperE, MB_TYPE_DOUBLE, fieldTag, MB_TAG_DENSE | MB_TAG_CREAT, &defDouble );
 std::cout << "rval = " << rval << std::endl;
 assert( MB_SUCCESS == rval || rval == MB_ALREADY_ALLOCATED );
#ifdef NDEBUG
 if( MB_SUCCESS == rval || rval == MB_ALREADY_ALLOCATED )
 {
 }; // Line to avoid compiler warning about unused variable
#endif
 
 int numElems = elems.size();
 std::vector< double > tData( ndofperE, 0.0 );
 
 for( int i = 0; i < numElems; i++ )
 {
 EntityHandle elem = elems[i];
 
 double xyz[3];
 mbi->get_coords( &elem, 1, xyz );
 
 double fieldValue = physField( xyz[0], xyz[1], xyz[2], factor );
 for( int j = 0; j < ndofperE; ++j )
 tData[j] = fieldValue; // same value evaluated at the center; could modify it to use
 // GLL points ?
 
 mbi->tag_set_data( fieldTag, &elem, 1, &tData[0] );
 }
 
 delete[] defDouble;
 tData.clear();
}
 
void putVertexField( Interface* mbi, const char* tagname, double factor )
{
 Range verts;
 
 mbi->get_entities_by_type( 0, MBVERTEX, verts );
 
 const double defVal = 0.;
 Tag fieldTag;
 mbi->tag_get_handle( tagname, 1, MB_TYPE_DOUBLE, fieldTag, MB_TAG_DENSE | MB_TAG_CREAT, &defVal );
 
 int numVerts = verts.size();
 for( int i = 0; i < numVerts; i++ )
 {
 EntityHandle vert = verts[i]; //?
 
 double vertPos[3];
 mbi->get_coords( &vert, 1, vertPos );
 
 double fieldValue = physField( vertPos[0], vertPos[1], vertPos[2], factor );
 
 mbi->tag_set_data( fieldTag, &vert, 1, &fieldValue );
 }
}
 
} // namespace moab
 
// Using moab instead of imesh for dev speed
int main( int argc, char** argv )
{
 
 using namespace moab;
 
 Interface* mbi = new Core();
 
 if( argc < 3 )
 {
 cout << "Usage: " << argv[0] << " <infile> <outfile> [factor]\n"
 << "Writes both vertex and element fields.\n";
 return 0;
 }
 
 mbi->load_mesh( argv[1] );
 
 double factor = 1.0;
 if( argc == 4 ) factor = atof( argv[3] );
 
 putVertexField( mbi, "vertex_field", factor );
 putElementField( mbi, "element_field", factor );
 // putSpectralElementField(mbi, 2, 4, "spectral_element_field", factor);
 putSpectralElementField( mbi, 2, 4, "a2oTAG", factor );
 
 ErrorCode result = mbi->write_mesh( argv[2] );
 if( MB_SUCCESS == result )
 cout << "wrote " << argv[2] << endl;
 else
 cout << "Failed to write " << argv[2] << endl;
 
 // vector<double> coords;
 // mbi->get_vertex_coordinates(coords);
 // double xavg = 0;
 // for (int i = 0; i < coords.size()/3; i++) xavg += coords[i];
 // cout << xavg << endl;
 
 return 1;
}