create_dp.cpp

Go to the documentation of this file.

/*
 * create_dp.cpp
 *
 * Created on: Dec 12, 2013
 * just to add a "DP" tag to an already partitioned mesh, based on some formula
 * it is launched in serial
 */
#include "moab/Core.hpp"
#include "moab/Interface.hpp"
#include <iostream>
#include <cmath>
#include <cassert>
 
#include "moab/IntxMesh/IntxUtils.hpp"
#include "IntxUtilsCSLAM.hpp"
 
#include <TestUtil.hpp>
 
using namespace moab;
 
double radius = 1.; // in m: 6371220.
int field_type = 1;
 
ErrorCode add_field_value( Interface& mb )
{
 Tag tagTracer = 0;
 std::string tag_name( "Tracer" );
 ErrorCode rval = mb.tag_get_handle( tag_name.c_str(), 1, MB_TYPE_DOUBLE, tagTracer, MB_TAG_DENSE | MB_TAG_CREAT );CHECK_ERR( rval );
 
 // tagElem is the average computed at each element, from nodal values
 Tag tagElem = 0;
 std::string tag_name2( "TracerAverage" );
 rval = mb.tag_get_handle( tag_name2.c_str(), 1, MB_TYPE_DOUBLE, tagElem, MB_TAG_DENSE | MB_TAG_CREAT );CHECK_ERR( rval );
 
 Tag tagArea = 0;
 std::string tag_name4( "Area" );
 rval = mb.tag_get_handle( tag_name4.c_str(), 1, MB_TYPE_DOUBLE, tagArea, MB_TAG_DENSE | MB_TAG_CREAT );CHECK_ERR( rval );
 
 /*
 * get all plys first, then vertices, then move them on the surface of the sphere
 * radius is 1., most of the time
 *
 */
 Range polygons;
 rval = mb.get_entities_by_dimension( 0, 2, polygons );
 if( MB_SUCCESS != rval ) return rval;
 
 Range connecVerts;
 rval = mb.get_connectivity( polygons, connecVerts );
 if( MB_SUCCESS != rval ) return rval;
 
 void* data; // pointer to the LOC in memory, for each vertex
 int count;
 
 rval = mb.tag_iterate( tagTracer, connecVerts.begin(), connecVerts.end(), count, data );CHECK_ERR( rval );
 // here we are checking contiguity
 assert( count == (int)connecVerts.size() );
 double* ptr_DP = (double*)data;
 // lambda is for longitude, theta for latitude
 // param will be: (la1, te1), (la2, te2), b, c; hmax=1, r=1/2
 // nondivergent flow, page 5, case 1, (la1, te1) = (M_PI, M_PI/3)
 // (la2, te2) = (M_PI, -M_PI/3)
 // la1, te1 la2 te2 b c hmax r
 if( field_type == 1 ) // quasi smooth
 {
 double params[] = { M_PI, M_PI / 3, M_PI, -M_PI / 3, 0.1, 0.9, 1., 0.5 };
 for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
 {
 EntityHandle oldV = *vit;
 CartVect posi;
 rval = mb.get_coords( &oldV, 1, &( posi[0] ) );CHECK_ERR( rval );
 
 moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( posi );
 
 ptr_DP[0] = IntxUtilsCSLAM::quasi_smooth_field( sphCoord.lon, sphCoord.lat, params );
 ;
 
 ptr_DP++; // increment to the next node
 }
 }
 else if( 2 == field_type ) // smooth
 {
 CartVect p1, p2;
 moab::IntxUtils::SphereCoords spr;
 spr.R = 1;
 spr.lat = M_PI / 3;
 spr.lon = M_PI;
 p1 = moab::IntxUtils::spherical_to_cart( spr );
 spr.lat = -M_PI / 3;
 p2 = moab::IntxUtils::spherical_to_cart( spr );
 // x1, y1, z1, x2, y2, z2, h_max, b0
 double params[] = { p1[0], p1[1], p1[2], p2[0], p2[1], p2[2], 1, 5. };
 for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
 {
 EntityHandle oldV = *vit;
 CartVect posi;
 rval = mb.get_coords( &oldV, 1, &( posi[0] ) );CHECK_ERR( rval );
 
 moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( posi );
 
 ptr_DP[0] = IntxUtilsCSLAM::smooth_field( sphCoord.lon, sphCoord.lat, params );
 ;
 
 ptr_DP++; // increment to the next node
 }
 }
 else if( 3 == field_type ) // slotted
 {
 // la1, te1, la2, te2, b, c, r
 double params[] = { M_PI, M_PI / 3, M_PI, -M_PI / 3, 0.1, 0.9, 0.5 }; // no h_max
 for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
 {
 EntityHandle oldV = *vit;
 CartVect posi;
 rval = mb.get_coords( &oldV, 1, &( posi[0] ) );CHECK_ERR( rval );
 
 moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( posi );
 
 ptr_DP[0] = IntxUtilsCSLAM::slotted_cylinder_field( sphCoord.lon, sphCoord.lat, params );
 ;
 
 ptr_DP++; // increment to the next node
 }
 }
 
 // add average value for quad/polygon (average corners)
 // do some averages
 
 Range::iterator iter = polygons.begin();
 // double local_mass = 0.; // this is total mass on one proc
 while( iter != polygons.end() )
 {
 rval = mb.tag_iterate( tagElem, iter, polygons.end(), count, data );CHECK_ERR( rval );
 double* ptr = (double*)data;
 
 rval = mb.tag_iterate( tagArea, iter, polygons.end(), count, data );CHECK_ERR( rval );
 double* ptrArea = (double*)data;
 for( int i = 0; i < count; i++, ++iter, ptr++, ptrArea++ )
 {
 const moab::EntityHandle* conn = NULL;
 int num_nodes = 0;
 rval = mb.get_connectivity( *iter, conn, num_nodes );CHECK_ERR( rval );
 if( num_nodes == 0 ) return MB_FAILURE;
 std::vector< double > nodeVals( num_nodes );
 double average = 0.;
 rval = mb.tag_get_data( tagTracer, conn, num_nodes, &nodeVals[0] );CHECK_ERR( rval );
 for( int j = 0; j < num_nodes; j++ )
 average += nodeVals[j];
 average /= num_nodes;
 *ptr = average;
 
 // now get area
 std::vector< double > coords;
 coords.resize( 3 * num_nodes );
 rval = mb.get_coords( conn, num_nodes, &coords[0] );CHECK_ERR( rval );
 IntxAreaUtils sphAreaUtils;
 *ptrArea = sphAreaUtils.area_spherical_polygon( &coords[0], num_nodes, radius );
 
 // we should have used some
 // total mass:
 // local_mass += *ptrArea * average;
 }
 }
 
 // now we can delete the tags? not yet
 return MB_SUCCESS;
}
 
int main( int argc, char** argv )
{
 
 if( argc < 3 )
 {
 std::cout << " usage: create_dp <input> <output> -t <time> -dt <delta_t> [-skipdp] -f "
 "<field> -h \n";
 return 1;
 }
 
 bool skip = false;
 double dt = 0.1;
 double t = 0.1; // corresponding to diffusion first step
 
 int index = 2;
 char* input_mesh1 = argv[1];
 char* output = argv[2];
 while( index < argc )
 {
 if( !strcmp( argv[index], "-t" ) ) // this is for radius to project
 {
 t = atof( argv[++index] );
 }
 if( !strcmp( argv[index], "-dt" ) ) // delete partition sets
 {
 dt = atof( argv[++index] );
 }
 
 if( !strcmp( argv[index], "-h" ) )
 {
 std::cout << " usage: create_dp <input> <output> -t <time> -dt <delta_t> [-skipdp] "
 "-f <field> -h \n";
 return 1;
 }
 if( !strcmp( argv[index], "-f" ) ) // delete partition sets
 {
 field_type = atoi( argv[++index] );
 }
 
 if( !strcmp( argv[index], "-skipdp" ) )
 {
 skip = true;
 }
 index++;
 }
 
 Core moab;
 Interface& mb = moab;
 
 ErrorCode rval;
 
 rval = mb.load_mesh( input_mesh1 );
 
 std::cout << " -t " << t << " -dt " << dt << " input: " << input_mesh1 << " output: " << output << "\n";
 
 // skip if we need for DP (already existing)
 if( skip )
 {
 std::cout << " do not add DP tag \n";
 }
 else
 {
 Range verts;
 rval = mb.get_entities_by_dimension( 0, 0, verts );
 if( MB_SUCCESS != rval ) return 1;
 
 double *x_ptr, *y_ptr, *z_ptr;
 int count;
 rval = mb.coords_iterate( verts.begin(), verts.end(), x_ptr, y_ptr, z_ptr, count );
 if( MB_SUCCESS != rval ) return 1;
 assert( count == (int)verts.size() ); // should end up with just one contiguous chunk of vertices
 
 Tag tagh = 0;
 std::string tag_name( "DP" );
 rval = mb.tag_get_handle( tag_name.c_str(), 3, MB_TYPE_DOUBLE, tagh, MB_TAG_DENSE | MB_TAG_CREAT );CHECK_ERR( rval );
 void* data; // pointer to the LOC in memory, for each vertex
 int count_tag;
 
 rval = mb.tag_iterate( tagh, verts.begin(), verts.end(), count_tag, data );CHECK_ERR( rval );
 // here we are checking contiguity
 assert( count_tag == (int)verts.size() );
 double* ptr_DP = (double*)data;
 
 for( int v = 0; v < count; v++ )
 {
 // EntityHandle v = verts[v];
 CartVect pos( x_ptr[v], y_ptr[v], z_ptr[v] );
 CartVect newPos;
 IntxUtilsCSLAM::departure_point_case1( pos, t, dt, newPos );
 ptr_DP[0] = newPos[0];
 ptr_DP[1] = newPos[1];
 ptr_DP[2] = newPos[2];
 ptr_DP += 3; // increment to the next vertex
 }
 }
 
 rval = add_field_value( mb );
 
 mb.write_file( output );
 
 return 0;
}