MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /* 00002 * create_dp.cpp 00003 * 00004 * Created on: Dec 12, 2013 00005 * just to add a "DP" tag to an already partitioned mesh, based on some formula 00006 * it is launched in serial 00007 */ 00008 #include "moab/Core.hpp" 00009 #include "moab/Interface.hpp" 00010 #include <iostream> 00011 #include <cmath> 00012 #include <cassert> 00013 00014 #include "moab/IntxMesh/IntxUtils.hpp" 00015 #include "IntxUtilsCSLAM.hpp" 00016 00017 #include <TestUtil.hpp> 00018 00019 using namespace moab; 00020 00021 double radius = 1.; // in m: 6371220. 00022 int field_type = 1; 00023 00024 ErrorCode add_field_value( Interface& mb ) 00025 { 00026 Tag tagTracer = 0; 00027 std::string tag_name( "Tracer" ); 00028 ErrorCode rval = mb.tag_get_handle( tag_name.c_str(), 1, MB_TYPE_DOUBLE, tagTracer, MB_TAG_DENSE | MB_TAG_CREAT );CHECK_ERR( rval ); 00029 00030 // tagElem is the average computed at each element, from nodal values 00031 Tag tagElem = 0; 00032 std::string tag_name2( "TracerAverage" ); 00033 rval = mb.tag_get_handle( tag_name2.c_str(), 1, MB_TYPE_DOUBLE, tagElem, MB_TAG_DENSE | MB_TAG_CREAT );CHECK_ERR( rval ); 00034 00035 Tag tagArea = 0; 00036 std::string tag_name4( "Area" ); 00037 rval = mb.tag_get_handle( tag_name4.c_str(), 1, MB_TYPE_DOUBLE, tagArea, MB_TAG_DENSE | MB_TAG_CREAT );CHECK_ERR( rval ); 00038 00039 /* 00040 * get all plys first, then vertices, then move them on the surface of the sphere 00041 * radius is 1., most of the time 00042 * 00043 */ 00044 Range polygons; 00045 rval = mb.get_entities_by_dimension( 0, 2, polygons ); 00046 if( MB_SUCCESS != rval ) return rval; 00047 00048 Range connecVerts; 00049 rval = mb.get_connectivity( polygons, connecVerts ); 00050 if( MB_SUCCESS != rval ) return rval; 00051 00052 void* data; // pointer to the LOC in memory, for each vertex 00053 int count; 00054 00055 rval = mb.tag_iterate( tagTracer, connecVerts.begin(), connecVerts.end(), count, data );CHECK_ERR( rval ); 00056 // here we are checking contiguity 00057 assert( count == (int)connecVerts.size() ); 00058 double* ptr_DP = (double*)data; 00059 // lambda is for longitude, theta for latitude 00060 // param will be: (la1, te1), (la2, te2), b, c; hmax=1, r=1/2 00061 // nondivergent flow, page 5, case 1, (la1, te1) = (M_PI, M_PI/3) 00062 // (la2, te2) = (M_PI, -M_PI/3) 00063 // la1, te1 la2 te2 b c hmax r 00064 if( field_type == 1 ) // quasi smooth 00065 { 00066 double params[] = { M_PI, M_PI / 3, M_PI, -M_PI / 3, 0.1, 0.9, 1., 0.5 }; 00067 for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit ) 00068 { 00069 EntityHandle oldV = *vit; 00070 CartVect posi; 00071 rval = mb.get_coords( &oldV, 1, &( posi[0] ) );CHECK_ERR( rval ); 00072 00073 moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( posi ); 00074 00075 ptr_DP[0] = IntxUtilsCSLAM::quasi_smooth_field( sphCoord.lon, sphCoord.lat, params ); 00076 ; 00077 00078 ptr_DP++; // increment to the next node 00079 } 00080 } 00081 else if( 2 == field_type ) // smooth 00082 { 00083 CartVect p1, p2; 00084 moab::IntxUtils::SphereCoords spr; 00085 spr.R = 1; 00086 spr.lat = M_PI / 3; 00087 spr.lon = M_PI; 00088 p1 = moab::IntxUtils::spherical_to_cart( spr ); 00089 spr.lat = -M_PI / 3; 00090 p2 = moab::IntxUtils::spherical_to_cart( spr ); 00091 // x1, y1, z1, x2, y2, z2, h_max, b0 00092 double params[] = { p1[0], p1[1], p1[2], p2[0], p2[1], p2[2], 1, 5. }; 00093 for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit ) 00094 { 00095 EntityHandle oldV = *vit; 00096 CartVect posi; 00097 rval = mb.get_coords( &oldV, 1, &( posi[0] ) );CHECK_ERR( rval ); 00098 00099 moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( posi ); 00100 00101 ptr_DP[0] = IntxUtilsCSLAM::smooth_field( sphCoord.lon, sphCoord.lat, params ); 00102 ; 00103 00104 ptr_DP++; // increment to the next node 00105 } 00106 } 00107 else if( 3 == field_type ) // slotted 00108 { 00109 // la1, te1, la2, te2, b, c, r 00110 double params[] = { M_PI, M_PI / 3, M_PI, -M_PI / 3, 0.1, 0.9, 0.5 }; // no h_max 00111 for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit ) 00112 { 00113 EntityHandle oldV = *vit; 00114 CartVect posi; 00115 rval = mb.get_coords( &oldV, 1, &( posi[0] ) );CHECK_ERR( rval ); 00116 00117 moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( posi ); 00118 00119 ptr_DP[0] = IntxUtilsCSLAM::slotted_cylinder_field( sphCoord.lon, sphCoord.lat, params ); 00120 ; 00121 00122 ptr_DP++; // increment to the next node 00123 } 00124 } 00125 00126 // add average value for quad/polygon (average corners) 00127 // do some averages 00128 00129 Range::iterator iter = polygons.begin(); 00130 // double local_mass = 0.; // this is total mass on one proc 00131 while( iter != polygons.end() ) 00132 { 00133 rval = mb.tag_iterate( tagElem, iter, polygons.end(), count, data );CHECK_ERR( rval ); 00134 double* ptr = (double*)data; 00135 00136 rval = mb.tag_iterate( tagArea, iter, polygons.end(), count, data );CHECK_ERR( rval ); 00137 double* ptrArea = (double*)data; 00138 for( int i = 0; i < count; i++, ++iter, ptr++, ptrArea++ ) 00139 { 00140 const moab::EntityHandle* conn = NULL; 00141 int num_nodes = 0; 00142 rval = mb.get_connectivity( *iter, conn, num_nodes );CHECK_ERR( rval ); 00143 if( num_nodes == 0 ) return MB_FAILURE; 00144 std::vector< double > nodeVals( num_nodes ); 00145 double average = 0.; 00146 rval = mb.tag_get_data( tagTracer, conn, num_nodes, &nodeVals[0] );CHECK_ERR( rval ); 00147 for( int j = 0; j < num_nodes; j++ ) 00148 average += nodeVals[j]; 00149 average /= num_nodes; 00150 *ptr = average; 00151 00152 // now get area 00153 std::vector< double > coords; 00154 coords.resize( 3 * num_nodes ); 00155 rval = mb.get_coords( conn, num_nodes, &coords[0] );CHECK_ERR( rval ); 00156 IntxAreaUtils sphAreaUtils; 00157 *ptrArea = sphAreaUtils.area_spherical_polygon( &coords[0], num_nodes, radius ); 00158 00159 // we should have used some 00160 // total mass: 00161 // local_mass += *ptrArea * average; 00162 } 00163 } 00164 00165 // now we can delete the tags? not yet 00166 return MB_SUCCESS; 00167 } 00168 00169 int main( int argc, char** argv ) 00170 { 00171 00172 if( argc < 3 ) 00173 { 00174 std::cout << " usage: create_dp <input> <output> -t <time> -dt <delta_t> [-skipdp] -f " 00175 "<field> -h \n"; 00176 return 1; 00177 } 00178 00179 bool skip = false; 00180 double dt = 0.1; 00181 double t = 0.1; // corresponding to diffusion first step 00182 00183 int index = 2; 00184 char* input_mesh1 = argv[1]; 00185 char* output = argv[2]; 00186 while( index < argc ) 00187 { 00188 if( !strcmp( argv[index], "-t" ) ) // this is for radius to project 00189 { t = atof( argv[++index] ); } 00190 if( !strcmp( argv[index], "-dt" ) ) // delete partition sets 00191 { dt = atof( argv[++index] ); } 00192 00193 if( !strcmp( argv[index], "-h" ) ) 00194 { 00195 std::cout << " usage: create_dp <input> <output> -t <time> -dt <delta_t> [-skipdp] " 00196 "-f <field> -h \n"; 00197 return 1; 00198 } 00199 if( !strcmp( argv[index], "-f" ) ) // delete partition sets 00200 { field_type = atoi( argv[++index] ); } 00201 00202 if( !strcmp( argv[index], "-skipdp" ) ) { skip = true; } 00203 index++; 00204 } 00205 00206 Core moab; 00207 Interface& mb = moab; 00208 00209 ErrorCode rval; 00210 00211 rval = mb.load_mesh( input_mesh1 ); 00212 00213 std::cout << " -t " << t << " -dt " << dt << " input: " << input_mesh1 << " output: " << output << "\n"; 00214 00215 // skip if we need for DP (already existing) 00216 if( skip ) { std::cout << " do not add DP tag \n"; } 00217 else 00218 { 00219 Range verts; 00220 rval = mb.get_entities_by_dimension( 0, 0, verts ); 00221 if( MB_SUCCESS != rval ) return 1; 00222 00223 double *x_ptr, *y_ptr, *z_ptr; 00224 int count; 00225 rval = mb.coords_iterate( verts.begin(), verts.end(), x_ptr, y_ptr, z_ptr, count ); 00226 if( MB_SUCCESS != rval ) return 1; 00227 assert( count == (int)verts.size() ); // should end up with just one contiguous chunk of vertices 00228 00229 Tag tagh = 0; 00230 std::string tag_name( "DP" ); 00231 rval = mb.tag_get_handle( tag_name.c_str(), 3, MB_TYPE_DOUBLE, tagh, MB_TAG_DENSE | MB_TAG_CREAT );CHECK_ERR( rval ); 00232 void* data; // pointer to the LOC in memory, for each vertex 00233 int count_tag; 00234 00235 rval = mb.tag_iterate( tagh, verts.begin(), verts.end(), count_tag, data );CHECK_ERR( rval ); 00236 // here we are checking contiguity 00237 assert( count_tag == (int)verts.size() ); 00238 double* ptr_DP = (double*)data; 00239 00240 for( int v = 0; v < count; v++ ) 00241 { 00242 // EntityHandle v = verts[v]; 00243 CartVect pos( x_ptr[v], y_ptr[v], z_ptr[v] ); 00244 CartVect newPos; 00245 IntxUtilsCSLAM::departure_point_case1( pos, t, dt, newPos ); 00246 ptr_DP[0] = newPos[0]; 00247 ptr_DP[1] = newPos[1]; 00248 ptr_DP[2] = newPos[2]; 00249 ptr_DP += 3; // increment to the next vertex 00250 } 00251 } 00252 00253 rval = add_field_value( mb ); 00254 00255 mb.write_file( output ); 00256 00257 return 0; 00258 }