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
13
14 double physField( double x, double y, double z, double factor )
15 {
16
17
18
19
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 };
44 #endif
45
46 int numElems = elems.size();
47
48 for( int i = 0; i < numElems; i++ )
49 {
50
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 };
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;
94
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 }
128
129
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
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
161
162
163
164
165
166 return 1;
167 }