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;
}