|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
#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) |
Definition at line 13 of file IntxUtilsCSLAM.hpp.
| ErrorCode IntxUtilsCSLAM::create_span_quads | ( | moab::Interface * | mb, |
| moab::EntityHandle | euler_set, | ||
| int | rank | ||
| ) | [static] |
Definition at line 155 of file IntxUtilsCSLAM.cpp.
References moab::Interface::add_entities(), moab::Range::begin(), moab::Interface::create_meshset(), moab::Interface::delete_entities(), moab::Range::end(), ErrorCode, moab::Interface::get_adjacencies(), moab::Interface::get_connectivity(), moab::ReadUtilIface::get_element_connect(), moab::Interface::get_entities_by_dimension(), moab::ReadUtilIface::get_node_coords(), 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::Interface::tag_get_data(), moab::Interface::tag_get_handle(), moab::Interface::tag_set_data(), moab::Interface::UNION, and moab::Interface::write_file().
{
// 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;
}
| ErrorCode IntxUtilsCSLAM::deep_copy_set | ( | moab::Interface * | mb, |
| moab::EntityHandle | source, | ||
| moab::EntityHandle | dest | ||
| ) | [static] |
Definition at line 273 of file IntxUtilsCSLAM.cpp.
References moab::Interface::add_entities(), moab::Range::begin(), CORRTAGNAME, moab::Interface::create_element(), moab::Interface::create_vertex(), moab::dum, moab::Range::end(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), moab::Interface::globalId_tag(), MB_CHK_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_HANDLE, moab::Interface::tag_get_data(), moab::Interface::tag_get_handle(), moab::Interface::tag_set_data(), and moab::Interface::type_from_handle().
Referenced by main().
{
// 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;
}
| 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.
References moab::IntxUtils::cart_to_spherical(), moab::IntxUtils::SphereCoords::lat, moab::IntxUtils::SphereCoords::lon, moab::IntxUtils::SphereCoords::R, moab::IntxUtils::spherical_to_cart(), and T.
Referenced by compute_tracer_case1(), main(), and manufacture_lagrange_mesh_on_sphere().
{
// 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;
}
| double IntxUtilsCSLAM::quasi_smooth_field | ( | double | lam, |
| double | tet, | ||
| double * | params | ||
| ) | [static] |
Definition at line 13 of file IntxUtilsCSLAM.cpp.
References moab::IntxUtils::distance_on_sphere().
Referenced by add_field_value(), and set_density().
{
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;
}
| double IntxUtilsCSLAM::slotted_cylinder_field | ( | double | lam, |
| double | tet, | ||
| double * | params | ||
| ) | [static] |
Definition at line 57 of file IntxUtilsCSLAM.cpp.
References moab::IntxUtils::distance_on_sphere().
Referenced by add_field_value(), and set_density().
{
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;
}
| double IntxUtilsCSLAM::smooth_field | ( | double | lam, |
| double | tet, | ||
| double * | params | ||
| ) | [static] |
Definition at line 40 of file IntxUtilsCSLAM.cpp.
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().
{
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 ) );
}
| void IntxUtilsCSLAM::velocity_case1 | ( | moab::CartVect & | arrival_point, |
| double | t, | ||
| moab::CartVect & | velo | ||
| ) | [static] |
Definition at line 134 of file IntxUtilsCSLAM.cpp.
References moab::IntxUtils::cart_to_spherical(), moab::IntxUtils::SphereCoords::lat, moab::IntxUtils::SphereCoords::lon, and T.
Referenced by compute_velocity_case1().
{
// 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;
}
1.7.6.1