IntxUtilsCSLAM.cpp

Go to the documentation of this file.

#include <cmath>
#include "moab/ReadUtilIface.hpp"
#include "IntxUtilsCSLAM.hpp"

using namespace moab;

#ifndef CORRTAGNAME
#define CORRTAGNAME "__correspondent"
#endif

// page 4 Nair Lauritzen paper
// param will be: (la1, te1), (la2, te2), b, c; hmax=1, r=1/2
double IntxUtilsCSLAM::quasi_smooth_field( double lam, double tet, double* params )
{
    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;
}

// page 4
// params are now x1, y1, ..., y2, z2 (6 params)
// plus h max and b0 (total of 8 params); radius is 1
double IntxUtilsCSLAM::smooth_field( double lam, double tet, double* params )
{
    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 ) );
}

// page 5
double IntxUtilsCSLAM::slotted_cylinder_field( double lam, double tet, double* params )
{
    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;
}

/*
 * based on paper A class of deformational flow test cases for linear transport problems on the
 * sphere longitude lambda [0.. 2*pi) and latitude theta (-pi/2, pi/2) lambda: -> lon (0, 2*pi)
 *  theta : -> lat (-pi/2, po/2)
 *  Radius of the sphere is 1 (if not, everything gets multiplied by 1)
 *
 *  cosine bell: center lambda0, theta0:
 */
void IntxUtilsCSLAM::departure_point_case1( CartVect& arrival_point,
                                            double t,
                                            double delta_t,
                                            CartVect& departure_point )
{
    // 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;
}

void IntxUtilsCSLAM::velocity_case1( CartVect& arrival_point, double t, CartVect& velo )
{
    // 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;
}

ErrorCode IntxUtilsCSLAM::create_span_quads( Interface* mb, EntityHandle euler_set, int rank )
{
    // 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;
}

// this simply copies the one mesh set into another, and sets some correlation tags
// for easy mapping back and forth
ErrorCode IntxUtilsCSLAM::deep_copy_set( Interface* mb, EntityHandle source_set, EntityHandle dest_set )
{
    // 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;
}