MOAB: Mesh Oriented datABase
(version 5.3.1)
|
00001 #include <iostream> 00002 #include "moab/Core.hpp" 00003 #include <cmath> 00004 #include <cstdlib> 00005 #include <cassert> 00006 00007 using namespace std; 00008 00009 namespace moab 00010 { 00011 00012 // make a nice picture. tweak here 00013 // 217: x,y -> [-8, 8] /// z -> [-65, 65] 00014 double physField( double x, double y, double z, double factor ) 00015 { 00016 00017 // 1/r^2 decay from {0,0,0} 00018 // tjt - changing to 1/r 00019 /* 00020 double scale = 1.0/factor; 00021 out = fabs(x*scale) + fabs(y*scale) + fabs(z*scale); 00022 out += 1e-1; // clamp 00023 out = 1/out; 00024 */ 00025 double out = factor * sqrt( x * x + y * y + z * z ); 00026 00027 return out; 00028 } 00029 00030 void putElementField( Interface* mbi, const char* tagname, double factor ) 00031 { 00032 Range elems; 00033 00034 mbi->get_entities_by_dimension( 0, 3, elems ); 00035 00036 const double defVal = 0.; 00037 Tag fieldTag; 00038 ErrorCode rval = mbi->tag_get_handle( tagname, 1, MB_TYPE_DOUBLE, fieldTag, MB_TAG_DENSE | MB_TAG_CREAT, &defVal ); 00039 assert( MB_SUCCESS == rval ); 00040 #ifdef NDEBUG 00041 if( MB_SUCCESS == rval ) {}; // Line to avoid compiler warning about unused variable 00042 #endif 00043 00044 int numElems = elems.size(); 00045 00046 for( int i = 0; i < numElems; i++ ) 00047 { 00048 // cout << elems[i] << endl; 00049 EntityHandle elem = elems[i]; 00050 00051 double xyz[3]; 00052 mbi->get_coords( &elem, 1, xyz ); 00053 00054 double fieldValue = physField( xyz[0], xyz[1], xyz[2], factor ); 00055 00056 mbi->tag_set_data( fieldTag, &elem, 1, &fieldValue ); 00057 } 00058 } 00059 00060 void putSpectralElementField( Interface* mbi, int dim, int np, const char* tagname, double factor ) 00061 { 00062 Range elems; 00063 00064 mbi->get_entities_by_dimension( 0, dim, elems ); 00065 00066 const int ndofperE = np * np; 00067 double* defDouble = new double[ndofperE]; 00068 Tag fieldTag; 00069 ErrorCode rval = 00070 mbi->tag_get_handle( tagname, ndofperE, MB_TYPE_DOUBLE, fieldTag, MB_TAG_DENSE | MB_TAG_CREAT, &defDouble ); 00071 std::cout << "rval = " << rval << std::endl; 00072 assert( MB_SUCCESS == rval || rval == MB_ALREADY_ALLOCATED ); 00073 #ifdef NDEBUG 00074 if( MB_SUCCESS == rval || rval == MB_ALREADY_ALLOCATED ) { 00075 }; // Line to avoid compiler warning about unused variable 00076 #endif 00077 00078 int numElems = elems.size(); 00079 std::vector< double > tData( ndofperE, 0.0 ); 00080 00081 for( int i = 0; i < numElems; i++ ) 00082 { 00083 EntityHandle elem = elems[i]; 00084 00085 double xyz[3]; 00086 mbi->get_coords( &elem, 1, xyz ); 00087 00088 double fieldValue = physField( xyz[0], xyz[1], xyz[2], factor ); 00089 for( int j = 0; j < ndofperE; ++j ) 00090 tData[j] = fieldValue; // same value evaluated at the center; could modify it to use 00091 // GLL points ? 00092 00093 mbi->tag_set_data( fieldTag, &elem, 1, &tData[0] ); 00094 } 00095 00096 delete[] defDouble; 00097 tData.clear(); 00098 } 00099 00100 void putVertexField( Interface* mbi, const char* tagname, double factor ) 00101 { 00102 Range verts; 00103 00104 mbi->get_entities_by_type( 0, MBVERTEX, verts ); 00105 00106 const double defVal = 0.; 00107 Tag fieldTag; 00108 mbi->tag_get_handle( tagname, 1, MB_TYPE_DOUBLE, fieldTag, MB_TAG_DENSE | MB_TAG_CREAT, &defVal ); 00109 00110 int numVerts = verts.size(); 00111 for( int i = 0; i < numVerts; i++ ) 00112 { 00113 EntityHandle vert = verts[i]; //? 00114 00115 double vertPos[3]; 00116 mbi->get_coords( &vert, 1, vertPos ); 00117 00118 double fieldValue = physField( vertPos[0], vertPos[1], vertPos[2], factor ); 00119 00120 mbi->tag_set_data( fieldTag, &vert, 1, &fieldValue ); 00121 } 00122 } 00123 00124 } // namespace moab 00125 00126 // Using moab instead of imesh for dev speed 00127 int main( int argc, char** argv ) 00128 { 00129 00130 using namespace moab; 00131 00132 Interface* mbi = new Core(); 00133 00134 if( argc < 3 ) 00135 { 00136 cout << "Usage: " << argv[0] << " <infile> <outfile> [factor]\n" 00137 << "Writes both vertex and element fields.\n"; 00138 return 0; 00139 } 00140 00141 mbi->load_mesh( argv[1] ); 00142 00143 double factor = 1.0; 00144 if( argc == 4 ) factor = atof( argv[3] ); 00145 00146 putVertexField( mbi, "vertex_field", factor ); 00147 putElementField( mbi, "element_field", factor ); 00148 // putSpectralElementField(mbi, 2, 4, "spectral_element_field", factor); 00149 putSpectralElementField( mbi, 2, 4, "a2oTAG", factor ); 00150 00151 ErrorCode result = mbi->write_mesh( argv[2] ); 00152 if( MB_SUCCESS == result ) 00153 cout << "wrote " << argv[2] << endl; 00154 else 00155 cout << "Failed to write " << argv[2] << endl; 00156 00157 // vector<double> coords; 00158 // mbi->get_vertex_coordinates(coords); 00159 // double xavg = 0; 00160 // for (int i = 0; i < coords.size()/3; i++) xavg += coords[i]; 00161 // cout << xavg << endl; 00162 00163 return 1; 00164 }