MOAB: Mesh Oriented datABase  (version 5.4.0)
IntxUtilsCSLAM.cpp
Go to the documentation of this file.
00001 #include <cmath>
00002 #include "moab/ReadUtilIface.hpp"
00003 #include "IntxUtilsCSLAM.hpp"
00004 
00005 using namespace moab;
00006 
00007 #ifndef CORRTAGNAME
00008 #define CORRTAGNAME "__correspondent"
00009 #endif
00010 
00011 // page 4 Nair Lauritzen paper
00012 // param will be: (la1, te1), (la2, te2), b, c; hmax=1, r=1/2
00013 double IntxUtilsCSLAM::quasi_smooth_field( double lam, double tet, double* params )
00014 {
00015     double la1   = params[0];
00016     double te1   = params[1];
00017     double la2   = params[2];
00018     double te2   = params[3];
00019     double b     = params[4];
00020     double c     = params[5];
00021     double hmax  = params[6];  // 1;
00022     double r     = params[7];  // 0.5;
00023     double r1    = moab::IntxUtils::distance_on_sphere( lam, tet, la1, te1 );
00024     double r2    = moab::IntxUtils::distance_on_sphere( lam, tet, la2, te2 );
00025     double value = b;
00026     if( r1 < r )
00027     {
00028         value += c * hmax / 2 * ( 1 + cos( M_PI * r1 / r ) );
00029     }
00030     if( r2 < r )
00031     {
00032         value += c * hmax / 2 * ( 1 + cos( M_PI * r2 / r ) );
00033     }
00034     return value;
00035 }
00036 
00037 // page 4
00038 // params are now x1, y1, ..., y2, z2 (6 params)
00039 // plus h max and b0 (total of 8 params); radius is 1
00040 double IntxUtilsCSLAM::smooth_field( double lam, double tet, double* params )
00041 {
00042     moab::IntxUtils::SphereCoords sc;
00043     sc.R         = 1.;
00044     sc.lat       = tet;
00045     sc.lon       = lam;
00046     double hmax  = params[6];
00047     double b0    = params[7];
00048     CartVect xyz = moab::IntxUtils::spherical_to_cart( sc );
00049     CartVect c1( params );
00050     CartVect c2( params + 3 );
00051     double expo1 = -b0 * ( xyz - c1 ).length_squared();
00052     double expo2 = -b0 * ( xyz - c2 ).length_squared();
00053     return hmax * ( exp( expo1 ) + exp( expo2 ) );
00054 }
00055 
00056 // page 5
00057 double IntxUtilsCSLAM::slotted_cylinder_field( double lam, double tet, double* params )
00058 {
00059     double la1 = params[0];
00060     double te1 = params[1];
00061     double la2 = params[2];
00062     double te2 = params[3];
00063     double b   = params[4];
00064     double c   = params[5];
00065     // double hmax = params[6]; // 1;
00066     double r      = params[6];  // 0.5;
00067     double r1     = moab::IntxUtils::distance_on_sphere( lam, tet, la1, te1 );
00068     double r2     = moab::IntxUtils::distance_on_sphere( lam, tet, la2, te2 );
00069     double value  = b;
00070     double d1     = fabs( lam - la1 );
00071     double d2     = fabs( lam - la2 );
00072     double rp6    = r / 6;
00073     double rt5p12 = r * 5 / 12;
00074 
00075     if( r1 <= r && d1 >= rp6 ) value = c;
00076     if( r2 <= r && d2 >= rp6 ) value = c;
00077     if( r1 <= r && d1 < rp6 && tet - te1 < -rt5p12 ) value = c;
00078     if( r2 <= r && d2 < rp6 && tet - te2 > rt5p12 ) value = c;
00079 
00080     return value;
00081 }
00082 
00083 /*
00084  * based on paper A class of deformational flow test cases for linear transport problems on the
00085  * sphere longitude lambda [0.. 2*pi) and latitude theta (-pi/2, pi/2) lambda: -> lon (0, 2*pi)
00086  *  theta : -> lat (-pi/2, po/2)
00087  *  Radius of the sphere is 1 (if not, everything gets multiplied by 1)
00088  *
00089  *  cosine bell: center lambda0, theta0:
00090  */
00091 void IntxUtilsCSLAM::departure_point_case1( CartVect& arrival_point,
00092                                             double t,
00093                                             double delta_t,
00094                                             CartVect& departure_point )
00095 {
00096     // always assume radius is 1 here?
00097     moab::IntxUtils::SphereCoords sph = moab::IntxUtils::cart_to_spherical( arrival_point );
00098     double T                          = 5;    // duration of integration (5 units)
00099     double k                          = 2.4;  // flow parameter
00100     /*     radius needs to be within some range   */
00101     double sl2      = sin( sph.lon / 2 );
00102     double pit      = M_PI * t / T;
00103     double omega    = M_PI / T;
00104     double costetha = cos( sph.lat );
00105     // double u = k * sl2*sl2 * sin(2*sph.lat) * cos(pit);
00106     double v = k * sin( sph.lon ) * costetha * cos( pit );
00107     // double psi = k * sl2 * sl2 *costetha * costetha * cos(pit);
00108     double u_tilda = 2 * k * sl2 * sl2 * sin( sph.lat ) * cos( pit );
00109 
00110     // formula 35, page 8
00111     // this will approximate dep point using a Taylor series with up to second derivative
00112     // this will be O(delta_t^3) exact.
00113     double lon_dep =
00114         sph.lon - delta_t * u_tilda -
00115         delta_t * delta_t * k * sl2 *
00116             ( sl2 * sin( sph.lat ) * sin( pit ) * omega - u_tilda * sin( sph.lat ) * cos( pit ) * cos( sph.lon / 2 ) -
00117               v * sl2 * costetha * cos( pit ) );
00118     // formula 36, page 8 again
00119     double lat_dep = sph.lat - delta_t * v -
00120                      delta_t * delta_t / 4 * k *
00121                          ( sin( sph.lon ) * cos( sph.lat ) * sin( pit ) * omega -
00122                            u_tilda * cos( sph.lon ) * cos( sph.lat ) * cos( pit ) +
00123                            v * sin( sph.lon ) * sin( sph.lat ) * cos( pit ) );
00124 
00125     moab::IntxUtils::SphereCoords sph_dep;
00126     sph_dep.R   = 1.;  // radius
00127     sph_dep.lat = lat_dep;
00128     sph_dep.lon = lon_dep;
00129 
00130     departure_point = moab::IntxUtils::spherical_to_cart( sph_dep );
00131     return;
00132 }
00133 
00134 void IntxUtilsCSLAM::velocity_case1( CartVect& arrival_point, double t, CartVect& velo )
00135 {
00136     // always assume radius is 1 here?
00137     moab::IntxUtils::SphereCoords sph = moab::IntxUtils::cart_to_spherical( arrival_point );
00138     double T                          = 5;    // duration of integration (5 units)
00139     double k                          = 2.4;  // flow parameter
00140     /*     radius needs to be within some range   */
00141     double sl2 = sin( sph.lon / 2 );
00142     double pit = M_PI * t / T;
00143     // double omega = M_PI/T;
00144     double coslat = cos( sph.lat );
00145     double sinlat = sin( sph.lat );
00146     double sinlon = sin( sph.lon );
00147     double coslon = cos( sph.lon );
00148     double u      = k * sl2 * sl2 * sin( 2 * sph.lat ) * cos( pit );
00149     double v      = k / 2 * sinlon * coslat * cos( pit );
00150     velo[0]       = -u * sinlon - v * sinlat * coslon;
00151     velo[1]       = u * coslon - v * sinlat * sinlon;
00152     velo[2]       = v * coslat;
00153 }
00154 
00155 ErrorCode IntxUtilsCSLAM::create_span_quads( Interface* mb, EntityHandle euler_set, int rank )
00156 {
00157     // first get all edges adjacent to polygons
00158     Tag dpTag = 0;
00159     std::string tag_name( "DP" );
00160     ErrorCode rval = mb->tag_get_handle( tag_name.c_str(), 3, MB_TYPE_DOUBLE, dpTag, MB_TAG_DENSE );
00161     // if the tag does not exist, get out early
00162     if( rval != MB_SUCCESS ) return rval;
00163     Range polygons;
00164     rval = mb->get_entities_by_dimension( euler_set, 2, polygons );
00165     if( MB_SUCCESS != rval ) return rval;
00166     Range iniEdges;
00167     rval = mb->get_adjacencies( polygons, 1, false, iniEdges, Interface::UNION );
00168     if( MB_SUCCESS != rval ) return rval;
00169     // now create some if missing
00170     Range allEdges;
00171     rval = mb->get_adjacencies( polygons, 1, true, allEdges, Interface::UNION );
00172     if( MB_SUCCESS != rval ) return rval;
00173     // create the vertices at the DP points, and the quads after that
00174     Range verts;
00175     rval = mb->get_connectivity( polygons, verts );
00176     if( MB_SUCCESS != rval ) return rval;
00177     int num_verts = (int)verts.size();
00178     // now see the departure points; to what boxes should we send them?
00179     std::vector< double > dep_points( 3 * num_verts );
00180     rval = mb->tag_get_data( dpTag, verts, (void*)&dep_points[0] );
00181     if( MB_SUCCESS != rval ) return rval;
00182 
00183     // create vertices corresponding to dp locations
00184     ReadUtilIface* read_iface;
00185     rval = mb->query_interface( read_iface );
00186     if( MB_SUCCESS != rval ) return rval;
00187     std::vector< double* > coords;
00188     EntityHandle start_vert, start_elem, *connect;
00189     // create verts, num is 2(nquads+1) because they're in a 1d row; will initialize coords in loop
00190     // over quads later
00191     rval = read_iface->get_node_coords( 3, num_verts, 0, start_vert, coords );
00192     if( MB_SUCCESS != rval ) return rval;
00193     // fill it up
00194     // Cppcheck warning (false positive): variable coords is assigned a value that is never used
00195     for( int i = 0; i < num_verts; i++ )
00196     {
00197         // block from interleaved
00198         coords[0][i] = dep_points[3 * i];
00199         coords[1][i] = dep_points[3 * i + 1];
00200         coords[2][i] = dep_points[3 * i + 2];
00201     }
00202     // create quads; one quad for each edge
00203     rval = read_iface->get_element_connect( allEdges.size(), 4, MBQUAD, 0, start_elem, connect );
00204     if( MB_SUCCESS != rval ) return rval;
00205 
00206     const EntityHandle* edge_conn = NULL;
00207     int quad_index                = 0;
00208     EntityHandle firstVertHandle  = verts[0];  // assume vertices are contiguous...
00209     for( Range::iterator eit = allEdges.begin(); eit != allEdges.end(); ++eit, quad_index++ )
00210     {
00211         EntityHandle edge = *eit;
00212         int num_nodes;
00213         rval = mb->get_connectivity( edge, edge_conn, num_nodes );
00214         if( MB_SUCCESS != rval ) return rval;
00215         connect[quad_index * 4]     = edge_conn[0];
00216         connect[quad_index * 4 + 1] = edge_conn[1];
00217 
00218         // maybe some indexing in range?
00219         connect[quad_index * 4 + 2] = start_vert + edge_conn[1] - firstVertHandle;
00220         connect[quad_index * 4 + 3] = start_vert + edge_conn[0] - firstVertHandle;
00221     }
00222 
00223     Range quads( start_elem, start_elem + allEdges.size() - 1 );
00224     EntityHandle outSet;
00225     rval = mb->create_meshset( MESHSET_SET, outSet );
00226     if( MB_SUCCESS != rval ) return rval;
00227     mb->add_entities( outSet, quads );
00228 
00229     Tag colTag;
00230     rval = mb->tag_get_handle( "COLOR_ID", 1, MB_TYPE_INTEGER, colTag, MB_TAG_DENSE | MB_TAG_CREAT );
00231     if( MB_SUCCESS != rval ) return rval;
00232     int j = 1;
00233     for( Range::iterator itq = quads.begin(); itq != quads.end(); ++itq, j++ )
00234     {
00235         EntityHandle q = *itq;
00236         rval           = mb->tag_set_data( colTag, &q, 1, &j );
00237         if( MB_SUCCESS != rval ) return rval;
00238     }
00239     std::stringstream outf;
00240     outf << "SpanQuads" << rank << ".h5m";
00241     rval = mb->write_file( outf.str().c_str(), 0, 0, &outSet, 1 );
00242     if( MB_SUCCESS != rval ) return rval;
00243     EntityHandle outSet2;
00244     rval = mb->create_meshset( MESHSET_SET, outSet2 );
00245     if( MB_SUCCESS != rval ) return rval;
00246 
00247     Range quadEdges;
00248     rval = mb->get_adjacencies( quads, 1, true, quadEdges, Interface::UNION );
00249     if( MB_SUCCESS != rval ) return rval;
00250     mb->add_entities( outSet2, quadEdges );
00251 
00252     std::stringstream outf2;
00253     outf2 << "SpanEdges" << rank << ".h5m";
00254     rval = mb->write_file( outf2.str().c_str(), 0, 0, &outSet2, 1 );
00255     if( MB_SUCCESS != rval ) return rval;
00256 
00257     // maybe some clean up
00258     mb->delete_entities( &outSet, 1 );
00259     mb->delete_entities( &outSet2, 1 );
00260     mb->delete_entities( quads );
00261     Range new_edges = subtract( allEdges, iniEdges );
00262     mb->delete_entities( new_edges );
00263     new_edges = subtract( quadEdges, iniEdges );
00264     mb->delete_entities( new_edges );
00265     Range new_verts( start_vert, start_vert + num_verts );
00266     mb->delete_entities( new_verts );
00267 
00268     return MB_SUCCESS;
00269 }
00270 
00271 // this simply copies the one mesh set into another, and sets some correlation tags
00272 // for easy mapping back and forth
00273 ErrorCode IntxUtilsCSLAM::deep_copy_set( Interface* mb, EntityHandle source_set, EntityHandle dest_set )
00274 {
00275     // create the handle tag for the corresponding element / vertex
00276 
00277     EntityHandle dum = 0;
00278     Tag corrTag      = 0;  // it will be created here
00279     ErrorCode rval   = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );MB_CHK_ERR( rval );
00280 
00281     // give the same global id to new verts and cells created in the lagr(departure) mesh
00282     Tag gid = mb->globalId_tag();
00283 
00284     Range polys;
00285     rval = mb->get_entities_by_dimension( source_set, 2, polys );MB_CHK_ERR( rval );
00286 
00287     Range connecVerts;
00288     rval = mb->get_connectivity( polys, connecVerts );MB_CHK_ERR( rval );
00289 
00290     std::map< EntityHandle, EntityHandle > newNodes;
00291     for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
00292     {
00293         EntityHandle oldV = *vit;
00294         CartVect posi;
00295         rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
00296         int global_id;
00297         rval = mb->tag_get_data( gid, &oldV, 1, &global_id );MB_CHK_ERR( rval );
00298 
00299         EntityHandle new_vert;
00300         // duplicate the position
00301         rval = mb->create_vertex( &( posi[0] ), new_vert );MB_CHK_ERR( rval );
00302         newNodes[oldV] = new_vert;
00303 
00304         // set also the correspondent tag :)
00305         rval = mb->tag_set_data( corrTag, &oldV, 1, &new_vert );MB_CHK_ERR( rval );
00306         // also the other side
00307         // need to check if we really need this; the new vertex will never need the old vertex
00308         // we have the global id which is the same
00309         rval = mb->tag_set_data( corrTag, &new_vert, 1, &oldV );MB_CHK_ERR( rval );
00310         // set the global id on the corresponding vertex the same as the initial vertex
00311         rval = mb->tag_set_data( gid, &new_vert, 1, &global_id );MB_CHK_ERR( rval );
00312     }
00313 
00314     for( Range::iterator it = polys.begin(); it != polys.end(); ++it )
00315     {
00316         EntityHandle q = *it;
00317         int nnodes;
00318         const EntityHandle* conn;
00319         rval = mb->get_connectivity( q, conn, nnodes );MB_CHK_ERR( rval );
00320 
00321         int global_id;
00322         rval = mb->tag_get_data( gid, &q, 1, &global_id );MB_CHK_ERR( rval );
00323         EntityType typeElem = mb->type_from_handle( q );
00324         std::vector< EntityHandle > new_conn( nnodes );
00325         for( int i = 0; i < nnodes; i++ )
00326         {
00327             EntityHandle v1 = conn[i];
00328             new_conn[i]     = newNodes[v1];
00329         }
00330         EntityHandle newElement;
00331         rval = mb->create_element( typeElem, &new_conn[0], nnodes, newElement );MB_CHK_ERR( rval );
00332 
00333         // set the corresponding tag; not sure we need this one, from old to new
00334         rval = mb->tag_set_data( corrTag, &q, 1, &newElement );MB_CHK_ERR( rval );
00335         rval = mb->tag_set_data( corrTag, &newElement, 1, &q );MB_CHK_ERR( rval );
00336 
00337         // set the global id
00338         rval = mb->tag_set_data( gid, &newElement, 1, &global_id );MB_CHK_ERR( rval );
00339         rval = mb->add_entities( dest_set, &newElement, 1 );MB_CHK_ERR( rval );
00340     }
00341 
00342     return MB_SUCCESS;
00343 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines