IntxUtilsCSLAM Class Reference

#include <IntxUtilsCSLAM.hpp>

Static Public Member Functions
static double	quasi_smooth_field (double lam, double tet, double *params)
static double	smooth_field (double lam, double tet, double *params)
static double	slotted_cylinder_field (double lam, double tet, double *params)
static void	departure_point_case1 (moab::CartVect &arrival_point, double t, double delta_t, moab::CartVect &departure_point)
static void	velocity_case1 (moab::CartVect &arrival_point, double t, moab::CartVect &velo)
static moab::ErrorCode	create_span_quads (moab::Interface *mb, moab::EntityHandle euler_set, int rank)
static moab::ErrorCode	deep_copy_set (moab::Interface *mb, moab::EntityHandle source, moab::EntityHandle dest)

Detailed Description

Definition at line 13 of file IntxUtilsCSLAM.hpp.

Member Function Documentation

create_span_quads()

ErrorCode IntxUtilsCSLAM::create_span_quads ( moab::Interface *  mb, moab::EntityHandle  euler_set, int  rank  )

static

Definition at line 155 of file IntxUtilsCSLAM.cpp.

{
 // first get all edges adjacent to polygons
 Tag dpTag = 0;
 std::string tag_name( "DP" );
 ErrorCode rval = mb->tag_get_handle( tag_name.c_str(), 3, MB_TYPE_DOUBLE, dpTag, MB_TAG_DENSE );
 // if the tag does not exist, get out early
 if( rval != MB_SUCCESS ) return rval;
 Range polygons;
 rval = mb->get_entities_by_dimension( euler_set, 2, polygons );
 if( MB_SUCCESS != rval ) return rval;
 Range iniEdges;
 rval = mb->get_adjacencies( polygons, 1, false, iniEdges, Interface::UNION );
 if( MB_SUCCESS != rval ) return rval;
 // now create some if missing
 Range allEdges;
 rval = mb->get_adjacencies( polygons, 1, true, allEdges, Interface::UNION );
 if( MB_SUCCESS != rval ) return rval;
 // create the vertices at the DP points, and the quads after that
 Range verts;
 rval = mb->get_connectivity( polygons, verts );
 if( MB_SUCCESS != rval ) return rval;
 int num_verts = (int)verts.size();
 // now see the departure points; to what boxes should we send them?
 std::vector< double > dep_points( 3 * num_verts );
 rval = mb->tag_get_data( dpTag, verts, (void*)&dep_points[0] );
 if( MB_SUCCESS != rval ) return rval;
 
 // create vertices corresponding to dp locations
 ReadUtilIface* read_iface;
 rval = mb->query_interface( read_iface );
 if( MB_SUCCESS != rval ) return rval;
 std::vector< double* > coords;
 EntityHandle start_vert, start_elem, *connect;
 // create verts, num is 2(nquads+1) because they're in a 1d row; will initialize coords in loop
 // over quads later
 rval = read_iface->get_node_coords( 3, num_verts, 0, start_vert, coords );
 if( MB_SUCCESS != rval ) return rval;
 // fill it up
 // Cppcheck warning (false positive): variable coords is assigned a value that is never used
 for( int i = 0; i < num_verts; i++ )
 {
 // block from interleaved
 coords[0][i] = dep_points[3 * i];
 coords[1][i] = dep_points[3 * i + 1];
 coords[2][i] = dep_points[3 * i + 2];
 }
 // create quads; one quad for each edge
 rval = read_iface->get_element_connect( allEdges.size(), 4, MBQUAD, 0, start_elem, connect );
 if( MB_SUCCESS != rval ) return rval;
 
 const EntityHandle* edge_conn = NULL;
 int quad_index = 0;
 EntityHandle firstVertHandle = verts[0]; // assume vertices are contiguous...
 for( Range::iterator eit = allEdges.begin(); eit != allEdges.end(); ++eit, quad_index++ )
 {
 EntityHandle edge = *eit;
 int num_nodes;
 rval = mb->get_connectivity( edge, edge_conn, num_nodes );
 if( MB_SUCCESS != rval ) return rval;
 connect[quad_index * 4] = edge_conn[0];
 connect[quad_index * 4 + 1] = edge_conn[1];
 
 // maybe some indexing in range?
 connect[quad_index * 4 + 2] = start_vert + edge_conn[1] - firstVertHandle;
 connect[quad_index * 4 + 3] = start_vert + edge_conn[0] - firstVertHandle;
 }
 
 Range quads( start_elem, start_elem + allEdges.size() - 1 );
 EntityHandle outSet;
 rval = mb->create_meshset( MESHSET_SET, outSet );
 if( MB_SUCCESS != rval ) return rval;
 mb->add_entities( outSet, quads );
 
 Tag colTag;
 rval = mb->tag_get_handle( "COLOR_ID", 1, MB_TYPE_INTEGER, colTag, MB_TAG_DENSE | MB_TAG_CREAT );
 if( MB_SUCCESS != rval ) return rval;
 int j = 1;
 for( Range::iterator itq = quads.begin(); itq != quads.end(); ++itq, j++ )
 {
 EntityHandle q = *itq;
 rval = mb->tag_set_data( colTag, &q, 1, &j );
 if( MB_SUCCESS != rval ) return rval;
 }
 std::stringstream outf;
 outf << "SpanQuads" << rank << ".h5m";
 rval = mb->write_file( outf.str().c_str(), 0, 0, &outSet, 1 );
 if( MB_SUCCESS != rval ) return rval;
 EntityHandle outSet2;
 rval = mb->create_meshset( MESHSET_SET, outSet2 );
 if( MB_SUCCESS != rval ) return rval;
 
 Range quadEdges;
 rval = mb->get_adjacencies( quads, 1, true, quadEdges, Interface::UNION );
 if( MB_SUCCESS != rval ) return rval;
 mb->add_entities( outSet2, quadEdges );
 
 std::stringstream outf2;
 outf2 << "SpanEdges" << rank << ".h5m";
 rval = mb->write_file( outf2.str().c_str(), 0, 0, &outSet2, 1 );
 if( MB_SUCCESS != rval ) return rval;
 
 // maybe some clean up
 mb->delete_entities( &outSet, 1 );
 mb->delete_entities( &outSet2, 1 );
 mb->delete_entities( quads );
 Range new_edges = subtract( allEdges, iniEdges );
 mb->delete_entities( new_edges );
 new_edges = subtract( quadEdges, iniEdges );
 mb->delete_entities( new_edges );
 Range new_verts( start_vert, start_vert + num_verts );
 mb->delete_entities( new_verts );
 
 return MB_SUCCESS;
}

References moab::Core::add_entities(), moab::Range::begin(), moab::Core::create_meshset(), moab::Core::delete_entities(), moab::Range::end(), ErrorCode, moab::Core::get_adjacencies(), moab::Core::get_connectivity(), moab::ReadUtilIface::get_element_connect(), moab::Core::get_entities_by_dimension(), moab::ReadUtilIface::get_node_coords(), mb, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_DOUBLE, MB_TYPE_INTEGER, MBQUAD, MESHSET_SET, moab::Interface::query_interface(), moab::Range::size(), moab::subtract(), moab::Core::tag_get_data(), moab::Core::tag_get_handle(), moab::Core::tag_set_data(), moab::Interface::UNION, and moab::Core::write_file().

deep_copy_set()

ErrorCode IntxUtilsCSLAM::deep_copy_set ( moab::Interface *  mb, moab::EntityHandle  source, moab::EntityHandle  dest  )

static

Definition at line 273 of file IntxUtilsCSLAM.cpp.

{
 // create the handle tag for the corresponding element / vertex
 
 EntityHandle dum = 0;
 Tag corrTag = 0; // it will be created here
 ErrorCode rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );MB_CHK_ERR( rval );
 
 // give the same global id to new verts and cells created in the lagr(departure) mesh
 Tag gid = mb->globalId_tag();
 
 Range polys;
 rval = mb->get_entities_by_dimension( source_set, 2, polys );MB_CHK_ERR( rval );
 
 Range connecVerts;
 rval = mb->get_connectivity( polys, connecVerts );MB_CHK_ERR( rval );
 
 std::map< EntityHandle, EntityHandle > newNodes;
 for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
 {
 EntityHandle oldV = *vit;
 CartVect posi;
 rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
 int global_id;
 rval = mb->tag_get_data( gid, &oldV, 1, &global_id );MB_CHK_ERR( rval );
 
 EntityHandle new_vert;
 // duplicate the position
 rval = mb->create_vertex( &( posi[0] ), new_vert );MB_CHK_ERR( rval );
 newNodes[oldV] = new_vert;
 
 // set also the correspondent tag :)
 rval = mb->tag_set_data( corrTag, &oldV, 1, &new_vert );MB_CHK_ERR( rval );
 // also the other side
 // need to check if we really need this; the new vertex will never need the old vertex
 // we have the global id which is the same
 rval = mb->tag_set_data( corrTag, &new_vert, 1, &oldV );MB_CHK_ERR( rval );
 // set the global id on the corresponding vertex the same as the initial vertex
 rval = mb->tag_set_data( gid, &new_vert, 1, &global_id );MB_CHK_ERR( rval );
 }
 
 for( Range::iterator it = polys.begin(); it != polys.end(); ++it )
 {
 EntityHandle q = *it;
 int nnodes;
 const EntityHandle* conn;
 rval = mb->get_connectivity( q, conn, nnodes );MB_CHK_ERR( rval );
 
 int global_id;
 rval = mb->tag_get_data( gid, &q, 1, &global_id );MB_CHK_ERR( rval );
 EntityType typeElem = mb->type_from_handle( q );
 std::vector< EntityHandle > new_conn( nnodes );
 for( int i = 0; i < nnodes; i++ )
 {
 EntityHandle v1 = conn[i];
 new_conn[i] = newNodes[v1];
 }
 EntityHandle newElement;
 rval = mb->create_element( typeElem, &new_conn[0], nnodes, newElement );MB_CHK_ERR( rval );
 
 // set the corresponding tag; not sure we need this one, from old to new
 rval = mb->tag_set_data( corrTag, &q, 1, &newElement );MB_CHK_ERR( rval );
 rval = mb->tag_set_data( corrTag, &newElement, 1, &q );MB_CHK_ERR( rval );
 
 // set the global id
 rval = mb->tag_set_data( gid, &newElement, 1, &global_id );MB_CHK_ERR( rval );
 rval = mb->add_entities( dest_set, &newElement, 1 );MB_CHK_ERR( rval );
 }
 
 return MB_SUCCESS;
}

References moab::Core::add_entities(), moab::Range::begin(), CORRTAGNAME, moab::Core::create_element(), moab::Core::create_vertex(), moab::dum, moab::Range::end(), ErrorCode, moab::Core::get_connectivity(), moab::Core::get_coords(), moab::Core::get_entities_by_dimension(), moab::Core::globalId_tag(), mb, MB_CHK_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_HANDLE, moab::Core::tag_get_data(), moab::Core::tag_get_handle(), moab::Core::tag_set_data(), and moab::Core::type_from_handle().

Referenced by main().

departure_point_case1()

void IntxUtilsCSLAM::departure_point_case1 ( moab::CartVect &  arrival_point, double  t, double  delta_t, moab::CartVect &  departure_point  )

static

Definition at line 91 of file IntxUtilsCSLAM.cpp.

{
 // always assume radius is 1 here?
 moab::IntxUtils::SphereCoords sph = moab::IntxUtils::cart_to_spherical( arrival_point );
 double T = 5; // duration of integration (5 units)
 double k = 2.4; // flow parameter
 /* radius needs to be within some range */
 double sl2 = sin( sph.lon / 2 );
 double pit = M_PI * t / T;
 double omega = M_PI / T;
 double costetha = cos( sph.lat );
 // double u = k * sl2*sl2 * sin(2*sph.lat) * cos(pit);
 double v = k * sin( sph.lon ) * costetha * cos( pit );
 // double psi = k * sl2 * sl2 *costetha * costetha * cos(pit);
 double u_tilda = 2 * k * sl2 * sl2 * sin( sph.lat ) * cos( pit );
 
 // formula 35, page 8
 // this will approximate dep point using a Taylor series with up to second derivative
 // this will be O(delta_t^3) exact.
 double lon_dep =
 sph.lon - delta_t * u_tilda -
 delta_t * delta_t * k * sl2 *
 ( sl2 * sin( sph.lat ) * sin( pit ) * omega - u_tilda * sin( sph.lat ) * cos( pit ) * cos( sph.lon / 2 ) -
 v * sl2 * costetha * cos( pit ) );
 // formula 36, page 8 again
 double lat_dep = sph.lat - delta_t * v -
 delta_t * delta_t / 4 * k *
 ( sin( sph.lon ) * cos( sph.lat ) * sin( pit ) * omega -
 u_tilda * cos( sph.lon ) * cos( sph.lat ) * cos( pit ) +
 v * sin( sph.lon ) * sin( sph.lat ) * cos( pit ) );
 
 moab::IntxUtils::SphereCoords sph_dep;
 sph_dep.R = 1.; // radius
 sph_dep.lat = lat_dep;
 sph_dep.lon = lon_dep;
 
 departure_point = moab::IntxUtils::spherical_to_cart( sph_dep );
 return;
}

References moab::IntxUtils::cart_to_spherical(), delta_t, moab::IntxUtils::SphereCoords::lat, moab::IntxUtils::SphereCoords::lon, moab::IntxUtils::SphereCoords::R, moab::IntxUtils::spherical_to_cart(), t, and T.

Referenced by compute_tracer_case1(), main(), and manufacture_lagrange_mesh_on_sphere().

quasi_smooth_field()

double IntxUtilsCSLAM::quasi_smooth_field ( double  lam, double  tet, double *  params  )

static

Definition at line 13 of file IntxUtilsCSLAM.cpp.

{
 double la1 = params[0];
 double te1 = params[1];
 double la2 = params[2];
 double te2 = params[3];
 double b = params[4];
 double c = params[5];
 double hmax = params[6]; // 1;
 double r = params[7]; // 0.5;
 double r1 = moab::IntxUtils::distance_on_sphere( lam, tet, la1, te1 );
 double r2 = moab::IntxUtils::distance_on_sphere( lam, tet, la2, te2 );
 double value = b;
 if( r1 < r )
 {
 value += c * hmax / 2 * ( 1 + cos( M_PI * r1 / r ) );
 }
 if( r2 < r )
 {
 value += c * hmax / 2 * ( 1 + cos( M_PI * r2 / r ) );
 }
 return value;
}

References moab::IntxUtils::distance_on_sphere().

Referenced by add_field_value(), and set_density().

slotted_cylinder_field()

double IntxUtilsCSLAM::slotted_cylinder_field ( double  lam, double  tet, double *  params  )

static

Definition at line 57 of file IntxUtilsCSLAM.cpp.

{
 double la1 = params[0];
 double te1 = params[1];
 double la2 = params[2];
 double te2 = params[3];
 double b = params[4];
 double c = params[5];
 // double hmax = params[6]; // 1;
 double r = params[6]; // 0.5;
 double r1 = moab::IntxUtils::distance_on_sphere( lam, tet, la1, te1 );
 double r2 = moab::IntxUtils::distance_on_sphere( lam, tet, la2, te2 );
 double value = b;
 double d1 = fabs( lam - la1 );
 double d2 = fabs( lam - la2 );
 double rp6 = r / 6;
 double rt5p12 = r * 5 / 12;
 
 if( r1 <= r && d1 >= rp6 ) value = c;
 if( r2 <= r && d2 >= rp6 ) value = c;
 if( r1 <= r && d1 < rp6 && tet - te1 < -rt5p12 ) value = c;
 if( r2 <= r && d2 < rp6 && tet - te2 > rt5p12 ) value = c;
 
 return value;
}

References moab::IntxUtils::distance_on_sphere().

Referenced by add_field_value(), and set_density().

smooth_field()

double IntxUtilsCSLAM::smooth_field ( double  lam, double  tet, double *  params  )

static

Definition at line 40 of file IntxUtilsCSLAM.cpp.

{
 moab::IntxUtils::SphereCoords sc;
 sc.R = 1.;
 sc.lat = tet;
 sc.lon = lam;
 double hmax = params[6];
 double b0 = params[7];
 CartVect xyz = moab::IntxUtils::spherical_to_cart( sc );
 CartVect c1( params );
 CartVect c2( params + 3 );
 double expo1 = -b0 * ( xyz - c1 ).length_squared();
 double expo2 = -b0 * ( xyz - c2 ).length_squared();
 return hmax * ( exp( expo1 ) + exp( expo2 ) );
}

References moab::IntxUtils::SphereCoords::lat, length_squared(), moab::IntxUtils::SphereCoords::lon, moab::IntxUtils::SphereCoords::R, and moab::IntxUtils::spherical_to_cart().

Referenced by add_field_value(), and set_density().

velocity_case1()

void IntxUtilsCSLAM::velocity_case1 ( moab::CartVect &  arrival_point, double  t, moab::CartVect &  velo  )

static

Definition at line 134 of file IntxUtilsCSLAM.cpp.

{
 // always assume radius is 1 here?
 moab::IntxUtils::SphereCoords sph = moab::IntxUtils::cart_to_spherical( arrival_point );
 double T = 5; // duration of integration (5 units)
 double k = 2.4; // flow parameter
 /* radius needs to be within some range */
 double sl2 = sin( sph.lon / 2 );
 double pit = M_PI * t / T;
 // double omega = M_PI/T;
 double coslat = cos( sph.lat );
 double sinlat = sin( sph.lat );
 double sinlon = sin( sph.lon );
 double coslon = cos( sph.lon );
 double u = k * sl2 * sl2 * sin( 2 * sph.lat ) * cos( pit );
 double v = k / 2 * sinlon * coslat * cos( pit );
 velo[0] = -u * sinlon - v * sinlat * coslon;
 velo[1] = u * coslon - v * sinlat * sinlon;
 velo[2] = v * coslat;
}

References moab::IntxUtils::cart_to_spherical(), moab::IntxUtils::SphereCoords::lat, moab::IntxUtils::SphereCoords::lon, t, and T.

Referenced by compute_velocity_case1().

---

The documentation for this class was generated from the following files:

• IntxUtilsCSLAM.hpp
• IntxUtilsCSLAM.cpp