MOAB: Mesh Oriented datABase  (version 5.3.0)
create_dp.cpp
Go to the documentation of this file.
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         {
00190             t = atof( argv[++index] );
00191         }
00192         if( !strcmp( argv[index], "-dt" ) )  // delete partition sets
00193         {
00194             dt = atof( argv[++index] );
00195         }
00196 
00197         if( !strcmp( argv[index], "-h" ) )
00198         {
00199             std::cout << " usage: create_dp <input> <output> -t <time>  -dt <delta_t>  [-skipdp] "
00200                          "-f <field>  -h \n";
00201             return 1;
00202         }
00203         if( !strcmp( argv[index], "-f" ) )  // delete partition sets
00204         {
00205             field_type = atoi( argv[++index] );
00206         }
00207 
00208         if( !strcmp( argv[index], "-skipdp" ) ) { skip = true; }
00209         index++;
00210     }
00211 
00212     Core moab;
00213     Interface& mb = moab;
00214 
00215     ErrorCode rval;
00216 
00217     rval = mb.load_mesh( input_mesh1 );
00218 
00219     std::cout << " -t " << t << " -dt " << dt << " input: " << input_mesh1 << "  output: " << output << "\n";
00220 
00221     // skip if we need for DP (already existing)
00222     if( skip ) { std::cout << " do not add DP tag \n"; }
00223     else
00224     {
00225         Range verts;
00226         rval = mb.get_entities_by_dimension( 0, 0, verts );
00227         if( MB_SUCCESS != rval ) return 1;
00228 
00229         double *x_ptr, *y_ptr, *z_ptr;
00230         int count;
00231         rval = mb.coords_iterate( verts.begin(), verts.end(), x_ptr, y_ptr, z_ptr, count );
00232         if( MB_SUCCESS != rval ) return 1;
00233         assert( count == (int)verts.size() );  // should end up with just one contiguous chunk of vertices
00234 
00235         Tag tagh = 0;
00236         std::string tag_name( "DP" );
00237         rval = mb.tag_get_handle( tag_name.c_str(), 3, MB_TYPE_DOUBLE, tagh, MB_TAG_DENSE | MB_TAG_CREAT );CHECK_ERR( rval );
00238         void* data;  // pointer to the LOC in memory, for each vertex
00239         int count_tag;
00240 
00241         rval = mb.tag_iterate( tagh, verts.begin(), verts.end(), count_tag, data );CHECK_ERR( rval );
00242         // here we are checking contiguity
00243         assert( count_tag == (int)verts.size() );
00244         double* ptr_DP = (double*)data;
00245 
00246         for( int v = 0; v < count; v++ )
00247         {
00248             // EntityHandle v = verts[v];
00249             CartVect pos( x_ptr[v], y_ptr[v], z_ptr[v] );
00250             CartVect newPos;
00251             IntxUtilsCSLAM::departure_point_case1( pos, t, dt, newPos );
00252             ptr_DP[0] = newPos[0];
00253             ptr_DP[1] = newPos[1];
00254             ptr_DP[2] = newPos[2];
00255             ptr_DP += 3;  // increment to the next vertex
00256         }
00257     }
00258 
00259     rval = add_field_value( mb );
00260 
00261     mb.write_file( output );
00262 
00263     return 0;
00264 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines