Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
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 }