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