Loading [MathJax]/extensions/tex2jax.js
Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
addfield.cpp
Go to the documentation of this file.
1 #include <iostream> 2 #include "moab/Core.hpp" 3 #include <cmath> 4 #include <cstdlib> 5 #include <cassert> 6  7 using namespace std; 8  9 namespace moab 10 { 11  12 // make a nice picture. tweak here 13 // 217: x,y -> [-8, 8] /// z -> [-65, 65] 14 double physField( double x, double y, double z, double factor ) 15 { 16  17  // 1/r^2 decay from {0,0,0} 18  // tjt - changing to 1/r 19  /* 20  double scale = 1.0/factor; 21  out = fabs(x*scale) + fabs(y*scale) + fabs(z*scale); 22  out += 1e-1; // clamp 23  out = 1/out; 24  */ 25  double out = factor * sqrt( x * x + y * y + z * z ); 26  27  return out; 28 } 29  30 void putElementField( Interface* mbi, const char* tagname, double factor ) 31 { 32  Range elems; 33  34  mbi->get_entities_by_dimension( 0, 3, elems ); 35  36  const double defVal = 0.; 37  Tag fieldTag; 38  ErrorCode rval = mbi->tag_get_handle( tagname, 1, MB_TYPE_DOUBLE, fieldTag, MB_TAG_DENSE | MB_TAG_CREAT, &defVal ); 39  assert( MB_SUCCESS == rval ); 40 #ifdef NDEBUG 41  if( MB_SUCCESS == rval ) 42  { 43  }; // Line to avoid compiler warning about unused variable 44 #endif 45  46  int numElems = elems.size(); 47  48  for( int i = 0; i < numElems; i++ ) 49  { 50  // cout << elems[i] << endl; 51  EntityHandle elem = elems[i]; 52  53  double xyz[3]; 54  mbi->get_coords( &elem, 1, xyz ); 55  56  double fieldValue = physField( xyz[0], xyz[1], xyz[2], factor ); 57  58  mbi->tag_set_data( fieldTag, &elem, 1, &fieldValue ); 59  } 60 } 61  62 void putSpectralElementField( Interface* mbi, int dim, int np, const char* tagname, double factor ) 63 { 64  Range elems; 65  66  mbi->get_entities_by_dimension( 0, dim, elems ); 67  68  const int ndofperE = np * np; 69  double* defDouble = new double[ndofperE]; 70  Tag fieldTag; 71  ErrorCode rval = 72  mbi->tag_get_handle( tagname, ndofperE, MB_TYPE_DOUBLE, fieldTag, MB_TAG_DENSE | MB_TAG_CREAT, &defDouble ); 73  std::cout << "rval = " << rval << std::endl; 74  assert( MB_SUCCESS == rval || rval == MB_ALREADY_ALLOCATED ); 75 #ifdef NDEBUG 76  if( MB_SUCCESS == rval || rval == MB_ALREADY_ALLOCATED ) 77  { 78  }; // Line to avoid compiler warning about unused variable 79 #endif 80  81  int numElems = elems.size(); 82  std::vector< double > tData( ndofperE, 0.0 ); 83  84  for( int i = 0; i < numElems; i++ ) 85  { 86  EntityHandle elem = elems[i]; 87  88  double xyz[3]; 89  mbi->get_coords( &elem, 1, xyz ); 90  91  double fieldValue = physField( xyz[0], xyz[1], xyz[2], factor ); 92  for( int j = 0; j < ndofperE; ++j ) 93  tData[j] = fieldValue; // same value evaluated at the center; could modify it to use 94  // GLL points ? 95  96  mbi->tag_set_data( fieldTag, &elem, 1, &tData[0] ); 97  } 98  99  delete[] defDouble; 100  tData.clear(); 101 } 102  103 void putVertexField( Interface* mbi, const char* tagname, double factor ) 104 { 105  Range verts; 106  107  mbi->get_entities_by_type( 0, MBVERTEX, verts ); 108  109  const double defVal = 0.; 110  Tag fieldTag; 111  mbi->tag_get_handle( tagname, 1, MB_TYPE_DOUBLE, fieldTag, MB_TAG_DENSE | MB_TAG_CREAT, &defVal ); 112  113  int numVerts = verts.size(); 114  for( int i = 0; i < numVerts; i++ ) 115  { 116  EntityHandle vert = verts[i]; //? 117  118  double vertPos[3]; 119  mbi->get_coords( &vert, 1, vertPos ); 120  121  double fieldValue = physField( vertPos[0], vertPos[1], vertPos[2], factor ); 122  123  mbi->tag_set_data( fieldTag, &vert, 1, &fieldValue ); 124  } 125 } 126  127 } // namespace moab 128  129 // Using moab instead of imesh for dev speed 130 int main( int argc, char** argv ) 131 { 132  133  using namespace moab; 134  135  Interface* mbi = new Core(); 136  137  if( argc < 3 ) 138  { 139  cout << "Usage: " << argv[0] << " <infile> <outfile> [factor]\n" 140  << "Writes both vertex and element fields.\n"; 141  return 0; 142  } 143  144  mbi->load_mesh( argv[1] ); 145  146  double factor = 1.0; 147  if( argc == 4 ) factor = atof( argv[3] ); 148  149  putVertexField( mbi, "vertex_field", factor ); 150  putElementField( mbi, "element_field", factor ); 151  // putSpectralElementField(mbi, 2, 4, "spectral_element_field", factor); 152  putSpectralElementField( mbi, 2, 4, "a2oTAG", factor ); 153  154  ErrorCode result = mbi->write_mesh( argv[2] ); 155  if( MB_SUCCESS == result ) 156  cout << "wrote " << argv[2] << endl; 157  else 158  cout << "Failed to write " << argv[2] << endl; 159  160  // vector<double> coords; 161  // mbi->get_vertex_coordinates(coords); 162  // double xavg = 0; 163  // for (int i = 0; i < coords.size()/3; i++) xavg += coords[i]; 164  // cout << xavg << endl; 165  166  return 1; 167 }