MOAB: Mesh Oriented datABase  (version 5.2.1)
linear_advection.cpp
Go to the documentation of this file.
00001 #include <iostream>
00002 #include <sstream>
00003 #include <ctime>
00004 #include <cstdlib>
00005 #include <string>
00006 #include <cmath>
00007 #include "moab/Core.hpp"
00008 #include "moab/Interface.hpp"
00009 #include "moab/IntxMesh/Intx2MeshOnSphere.hpp"
00010 #include "moab/ParallelComm.hpp"
00011 #include "moab/ProgOptions.hpp"
00012 #include "MBParallelConventions.h"
00013 #include "moab/ReadUtilIface.hpp"
00014 #include "MBTagConventions.hpp"
00015 #include "IntxUtilsCSLAM.hpp"
00016 
00017 #include "TestUtil.hpp"
00018 
00019 // std::string file_name("./uniform_30.g");
00020 // std::string file_name("./uniform_120.g");
00021 // std::string file_name("./eulerHomme.vtk");
00022 
00023 using namespace moab;
00024 
00025 moab::ErrorCode update_density( moab::Interface* mb, moab::EntityHandle euler_set, moab::EntityHandle lagr_set,
00026                                 moab::EntityHandle out_set, moab::Tag& rhoTag, moab::Tag& areaTag,
00027                                 moab::Tag& rhoCoefsTag, moab::Tag& weightsTag, moab::Tag& planeTag );
00028 
00029 moab::ErrorCode get_departure_grid( moab::Interface* mb, moab::EntityHandle euler_set, moab::EntityHandle lagr_set,
00030                                     moab::EntityHandle covering_set, int tStep, moab::Range& connecVerts );
00031 
00032 moab::ErrorCode get_barycenters( moab::Interface* mb, moab::EntityHandle set, moab::Tag& planeTag,
00033                                  moab::Tag& barycenterTag, moab::Tag& areaTag );
00034 moab::ErrorCode get_gnomonic_plane( moab::Interface* mb, moab::EntityHandle set, moab::Tag& planeTag );
00035 moab::ErrorCode get_linear_reconstruction( moab::Interface* mb, moab::EntityHandle set, moab::Tag& rhoTag,
00036                                            moab::Tag& planeTag, moab::Tag& barycenterTag, moab::Tag& linearCoefTag );
00037 moab::ErrorCode get_intersection_weights( moab::Interface* mb, moab::EntityHandle euler_set,
00038                                           moab::EntityHandle lagr_set, moab::EntityHandle intx_set, moab::Tag& planeTag,
00039                                           moab::Tag& weightsTag );
00040 
00041 moab::ErrorCode set_density( moab::Interface* mb, moab::EntityHandle euler_set, moab::Tag& barycenterTag,
00042                              moab::Tag& rhoTag, int field_type );
00043 
00044 moab::ErrorCode create_lagr_mesh( moab::Interface* mb, moab::EntityHandle euler_set, moab::EntityHandle lagr_set );
00045 
00046 // functions to compute departure point locations
00047 void departure_point_swirl( moab::CartVect& arrival_point, double t, double delta_t, moab::CartVect& departure_point );
00048 void departure_point_swirl_rot( moab::CartVect& arrival_point, double t, double delta_t,
00049                                 moab::CartVect& departure_point );
00050 
00051 double gtol = 1.e-9;  // this is for geometry tolerance
00052 
00053 double radius = 1.;
00054 
00055 bool writeFiles    = true;
00056 bool parallelWrite = false;
00057 bool velocity      = false;
00058 
00059 int numSteps = 200;  // number of times with velocity displayed at points
00060 double T     = 5;
00061 
00062 Intx2MeshOnSphere* pworker = NULL;
00063 
00064 int main( int argc, char* argv[] )
00065 {
00066 
00067     MPI_Init( &argc, &argv );
00068 
00069     // set up MOAB interface and parallel communication
00070     moab::Core moab;
00071     moab::Interface& mb = moab;
00072     moab::ParallelComm mb_pcomm( &mb, MPI_COMM_WORLD );
00073 
00074     // int rank = mb_pcomm->proc_config().proc_rank();
00075     int rank = mb_pcomm.proc_config().proc_rank();
00076 
00077     // create meshset
00078     moab::EntityHandle euler_set;
00079     moab::ErrorCode rval = mb.create_meshset( MESHSET_SET, euler_set );MB_CHK_ERR( rval );
00080 
00081     // std::stringstream opts;
00082     // opts <<
00083     // "PARALLEL=READ_PART;PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;GATHER_SET=0;PARTITION_METHOD=TRIVIAL_PARTITION;VARIABLE=";
00084     // opts << "PARALLEL=READ_PART;PARTITION;PARALLEL_RESOLVE_SHARED_ENTS";
00085     // rval = mb.load_file(file_name.c_str(), &euler_set, opts.str().c_str());
00086     std::string fileN = TestDir + "/mbcslam/fine4.h5m";
00087 
00088     rval = mb.load_file( fileN.c_str(), &euler_set );MB_CHK_ERR( rval );
00089 
00090     // Create tag for cell density
00091     moab::Tag rhoTag = 0;
00092     rval = mb.tag_get_handle( "Density", 1, moab::MB_TYPE_DOUBLE, rhoTag, moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );MB_CHK_ERR( rval );
00093 
00094     // Create tag for cell area
00095     moab::Tag areaTag = 0;
00096     rval = mb.tag_get_handle( "Area", 1, moab::MB_TYPE_DOUBLE, areaTag, moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );MB_CHK_ERR( rval );
00097     // Create tag for cell barycenters in 3D Cartesian space
00098     moab::Tag barycenterTag = 0;
00099     rval                    = mb.tag_get_handle( "CellBarycenter", 3, moab::MB_TYPE_DOUBLE, barycenterTag,
00100                               moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );MB_CHK_ERR( rval );
00101     // Create tag for cell density reconstruction coefficients
00102     moab::Tag rhoCoefTag = 0;
00103     rval                 = mb.tag_get_handle( "LinearCoefRho", 3, moab::MB_TYPE_DOUBLE, rhoCoefTag,
00104                               moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );MB_CHK_ERR( rval );
00105     // Create tag for index of gnomonic plane for each cell
00106     moab::Tag planeTag = 0;
00107     rval               = mb.tag_get_handle( "gnomonicPlane", 1, moab::MB_TYPE_INTEGER, planeTag,
00108                               moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );MB_CHK_ERR( rval );
00109     // Create tag for intersection weights
00110     moab::Tag weightsTag = 0;
00111     rval = mb.tag_get_handle( "Weights", 3, moab::MB_TYPE_DOUBLE, weightsTag, moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );MB_CHK_ERR( rval );
00112 
00113     // get cell plane
00114     rval = get_gnomonic_plane( &mb, euler_set, planeTag );MB_CHK_ERR( rval );
00115 
00116     // get cell barycenters (and cell area)
00117     rval = get_barycenters( &mb, euler_set, planeTag, areaTag, barycenterTag );MB_CHK_ERR( rval );
00118 
00119     // Set density distributions
00120     rval = set_density( &mb, euler_set, barycenterTag, rhoTag, 1 );MB_CHK_ERR( rval );
00121 
00122     // Get initial values for use in error computation
00123     moab::Range redEls;
00124     rval = mb.get_entities_by_dimension( euler_set, 2, redEls );MB_CHK_ERR( rval );
00125     std::vector< double > iniValsRho( redEls.size() );
00126     rval = mb.tag_get_data( rhoTag, redEls, &iniValsRho[0] );MB_CHK_ERR( rval );
00127 
00128     // Get Lagrangian set
00129     moab::EntityHandle out_set, lagrange_set, covering_set;
00130     rval = mb.create_meshset( MESHSET_SET, out_set );MB_CHK_ERR( rval );
00131     rval = mb.create_meshset( MESHSET_SET, lagrange_set );MB_CHK_ERR( rval );
00132     rval = mb.create_meshset( MESHSET_SET, covering_set );MB_CHK_ERR( rval );
00133     // rval = create_lagr_mesh(&mb, euler_set, lagrange_set);
00134     rval = IntxUtilsCSLAM::deep_copy_set( &mb, euler_set, lagrange_set );MB_CHK_ERR( rval );
00135     moab::EntityHandle dum = 0;
00136     moab::Tag corrTag;
00137     rval = mb.tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );MB_CHK_ERR( rval );
00138     // Set up intersection of two meshes
00139 
00140     /*
00141         moab::Intx2MeshOnSphere worker(&mb);
00142         worker.SetRadius(radius);
00143         worker.SetErrorTolerance(gtol);
00144     */
00145 
00146     pworker = new Intx2MeshOnSphere( &mb );
00147     pworker->set_error_tolerance( gtol );
00148     pworker->set_radius_source_mesh( radius );
00149     pworker->set_radius_destination_mesh( radius );
00150     pworker->set_box_error( 100 * gtol );
00151 
00152     rval = pworker->FindMaxEdges( euler_set, euler_set );MB_CHK_ERR( rval );
00153 
00154     // these stay fixed for one run
00155     moab::Range local_verts;
00156     rval = pworker->build_processor_euler_boxes( euler_set,
00157                                                  local_verts );  // output also the local_verts
00158     MB_CHK_ERR( rval );
00159     // rval = worker.build_processor_euler_boxes(euler_set, local_verts);// output also the
00160     // local_verts
00161 
00162     // loop over time to update density
00163     for( int ts = 1; ts < numSteps + 1; ts++ )
00164     {
00165 
00166         if( ts == 1 )  // output initial condition
00167         {
00168             std::stringstream newDensity;
00169             newDensity << "Density" << rank << "_" << ts - 1 << ".vtk";
00170             rval = mb.write_file( newDensity.str().c_str(), 0, 0, &euler_set, 1 );
00171         }
00172 
00173         // get linear reconstruction coefficients
00174         get_linear_reconstruction( &mb, euler_set, rhoTag, planeTag, barycenterTag, rhoCoefTag );
00175 
00176         // get depature grid
00177         rval = get_departure_grid( &mb, euler_set, lagrange_set, covering_set, ts, local_verts );MB_CHK_ERR( rval );
00178 
00179         // intersect the meshes
00180         rval = pworker->intersect_meshes( lagrange_set, euler_set, out_set );MB_CHK_ERR( rval );
00181 
00182         // intersection weights (i.e. area, x integral, and y integral over cell intersections)
00183         get_intersection_weights( &mb, euler_set, lagrange_set, out_set, planeTag, weightsTag );
00184 
00185         // update the density
00186         rval =
00187             update_density( &mb, euler_set, lagrange_set, out_set, rhoTag, areaTag, rhoCoefTag, weightsTag, planeTag );MB_CHK_ERR( rval );
00188 
00189         if( writeFiles && ( ts % 5 == 0 ) )  // so if write
00190         {
00191             std::stringstream newDensity;
00192             newDensity << "Density" << rank << "_" << ts << ".vtk";
00193             rval = mb.write_file( newDensity.str().c_str(), 0, 0, &euler_set, 1 );MB_CHK_ERR( rval );
00194         }
00195 
00196         // delete the polygons and elements of out_set
00197         moab::Range allVerts;
00198         rval = mb.get_entities_by_dimension( 0, 0, allVerts );MB_CHK_ERR( rval );
00199 
00200         moab::Range allElems;
00201         rval = mb.get_entities_by_dimension( 0, 2, allElems );MB_CHK_ERR( rval );
00202 
00203         // get Eulerian and lagrangian cells
00204         moab::Range polys;
00205         rval = mb.get_entities_by_dimension( euler_set, 2, polys );MB_CHK_ERR( rval );
00206         rval =
00207             mb.get_entities_by_dimension( lagrange_set, 2, polys );  // do not delete lagr set either, with its vertices
00208         MB_CHK_ERR( rval );
00209 
00210         // add to the connecVerts range all verts, from all initial polys
00211         moab::Range vertsToStay;
00212         rval = mb.get_connectivity( polys, vertsToStay );MB_CHK_ERR( rval );
00213 
00214         moab::Range todeleteVerts = subtract( allVerts, vertsToStay );
00215 
00216         moab::Range todeleteElem = subtract( allElems, polys );
00217         // empty the out mesh set
00218         rval = mb.clear_meshset( &out_set, 1 );MB_CHK_ERR( rval );
00219 
00220         rval = mb.delete_entities( todeleteElem );MB_CHK_ERR( rval );
00221         rval = mb.delete_entities( todeleteVerts );MB_CHK_ERR( rval );
00222         if( rank == 0 ) std::cout << " step: " << ts << "\n";
00223     }
00224 
00225     // final vals and errors
00226     moab::Range::iterator iter = redEls.begin();
00227     double norm1               = 0.;
00228     double norm2               = 0.;
00229     double exact2              = 0.;
00230     double exact1              = 0.;
00231     int count                  = 0;
00232     void* data;
00233     int j = 0;  // index in iniVals
00234     while( iter != redEls.end() )
00235     {
00236         rval = mb.tag_iterate( rhoTag, iter, redEls.end(), count, data );MB_CHK_ERR( rval );
00237         double* ptrTracer = (double*)data;
00238 
00239         rval = mb.tag_iterate( areaTag, iter, redEls.end(), count, data );MB_CHK_ERR( rval );
00240         double* ptrArea = (double*)data;
00241         for( int i = 0; i < count; i++, ++iter, ptrTracer++, ptrArea++, j++ )
00242         {
00243             // double area = *ptrArea;
00244             norm1 += fabs( *ptrTracer - iniValsRho[j] ) * ( *ptrArea );
00245             norm2 += ( *ptrTracer - iniValsRho[j] ) * ( *ptrTracer - iniValsRho[j] ) * ( *ptrArea );
00246             exact1 += ( iniValsRho[j] ) * ( *ptrArea );
00247             exact2 += ( iniValsRho[j] ) * ( iniValsRho[j] ) * ( *ptrArea );
00248         }
00249     }
00250 
00251     double total_norm1  = 0;
00252     double total_norm2  = 0;
00253     double total_exact1 = 0;
00254     double total_exact2 = 0;
00255     int mpi_err         = MPI_Reduce( &norm1, &total_norm1, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
00256     if( mpi_err ) std::cout << " error in MPI_reduce:" << mpi_err << "\n";
00257     mpi_err = MPI_Reduce( &norm2, &total_norm2, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
00258     if( mpi_err ) std::cout << " error in MPI_reduce:" << mpi_err << "\n";
00259     mpi_err = MPI_Reduce( &exact1, &total_exact1, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
00260     if( mpi_err ) std::cout << " error in MPI_reduce:" << mpi_err << "\n";
00261     mpi_err = MPI_Reduce( &exact2, &total_exact2, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
00262     if( mpi_err ) std::cout << " error in MPI_reduce:" << mpi_err << "\n";
00263     if( 0 == rank )
00264         std::cout << " numSteps:" << numSteps << " 1-norm:" << total_norm1 / total_exact1
00265                   << " 2-norm:" << total_norm2 / total_exact2 << "\n";
00266 
00267     MPI_Finalize();
00268     return 0;
00269 }
00270 
00271 moab::ErrorCode get_gnomonic_plane( moab::Interface* mb, moab::EntityHandle set, moab::Tag& planeTag )
00272 {
00273     // get all entities of dimension 2
00274     moab::Range cells;
00275     ErrorCode rval = mb->get_entities_by_dimension( set, 2, cells );MB_CHK_ERR( rval );
00276 
00277     for( Range::iterator it = cells.begin(); it != cells.end(); ++it )
00278     {
00279         moab::EntityHandle icell = *it;
00280 
00281         // get the nodes
00282         const moab::EntityHandle* verts;
00283         int num_nodes;
00284         rval = mb->get_connectivity( icell, verts, num_nodes );MB_CHK_ERR( rval );
00285 
00286         // get coordinates
00287         std::vector< double > coords( 3 * num_nodes );
00288         rval = mb->get_coords( verts, num_nodes, &coords[0] );MB_CHK_ERR( rval );
00289 
00290         // get cell center
00291         double centerx = 0;
00292         double centery = 0;
00293         double centerz = 0;
00294         for( int inode = 0; inode < num_nodes; inode++ )
00295         {
00296             centerx += coords[inode * 3] / num_nodes;
00297             centery += coords[inode * 3 + 1] / num_nodes;
00298             centerz += coords[inode * 3 + 2] / num_nodes;
00299         }
00300         double rad = std::sqrt( centerx * centerx + centery * centery + centerz * centerz );
00301         centerx    = centerx / rad;
00302         centery    = centery / rad;
00303         centerz    = centerz / rad;
00304         moab::CartVect center( centerx, centery, centerz );
00305 
00306         // define gnomonic plane based on cell center coordinates
00307         int plane = 0;
00308         moab::IntxUtils::decide_gnomonic_plane( center, plane );
00309 
00310         rval = mb->tag_set_data( planeTag, &icell, 1, &plane );MB_CHK_ERR( rval );
00311     }
00312 
00313     return moab::MB_SUCCESS;
00314 }
00315 
00316 moab::ErrorCode get_barycenters( moab::Interface* mb, moab::EntityHandle set, moab::Tag& planeTag, moab::Tag& areaTag,
00317                                  moab::Tag& barycenterTag )
00318 {
00319 
00320     // get all entities of dimension 2
00321     moab::Range cells;
00322     ErrorCode rval = mb->get_entities_by_dimension( set, 2, cells );MB_CHK_ERR( rval );
00323 
00324     // set sphere radius to 1
00325     double R = 1.0;
00326 
00327     for( Range::iterator it = cells.begin(); it != cells.end(); ++it )
00328     {
00329         moab::EntityHandle icell = *it;
00330 
00331         // get the nodes
00332         const moab::EntityHandle* verts;
00333         int num_nodes;
00334         rval = mb->get_connectivity( icell, verts, num_nodes );MB_CHK_ERR( rval );
00335 
00336         // get coordinates
00337         std::vector< double > coords( 3 * num_nodes );
00338         rval = mb->get_coords( verts, num_nodes, &coords[0] );MB_CHK_ERR( rval );
00339 
00340         // get gnomonic plane
00341         int plane = 0;
00342         rval      = mb->tag_get_data( planeTag, &icell, 1, &plane );MB_CHK_ERR( rval );
00343 
00344         // get vertex coordinates and project onto gnomonic plane
00345         std::vector< double > x( num_nodes );
00346         std::vector< double > y( num_nodes );
00347         double area   = 0;
00348         double bary_x = 0;
00349         double bary_y = 0;
00350         for( int inode = 0; inode < num_nodes; inode++ )
00351         {
00352             double rad = sqrt( coords[inode * 3] * coords[inode * 3] + coords[inode * 3 + 1] * coords[inode * 3 + 1] +
00353                                coords[inode * 3 + 2] * coords[inode * 3 + 2] );
00354             CartVect xyzcoord( coords[inode * 3] / rad, coords[inode * 3 + 1] / rad, coords[inode * 3 + 2] / rad );
00355             moab::IntxUtils::gnomonic_projection( xyzcoord, R, plane, x[inode], y[inode] );
00356         }
00357 
00358         // integrate over cell to get barycenter in gnomonic coordinates
00359         for( int inode = 0; inode < num_nodes; inode++ )
00360         {
00361             int inode2 = inode + 1;
00362             if( inode2 >= num_nodes ) inode2 = 0;
00363             double xmid = 0.5 * ( x[inode] + x[inode2] );
00364             double ymid = 0.5 * ( y[inode] + y[inode2] );
00365             double r1   = sqrt( 1 + x[inode] * x[inode] + y[inode] * y[inode] );
00366             double rm   = sqrt( 1 + xmid * xmid + ymid * ymid );
00367             double r2   = sqrt( 1 + x[inode2] * x[inode2] + y[inode2] * y[inode2] );
00368             double hx   = x[inode2] - x[inode];
00369 
00370             area += -hx *
00371                     ( y[inode] / ( r1 * ( 1 + x[inode] * x[inode] ) ) + 4.0 * ymid / ( rm * ( 1 + xmid * xmid ) ) +
00372                       y[inode2] / ( r2 * ( 1 + x[inode2] * x[inode2] ) ) ) /
00373                     6.0;
00374 
00375             bary_x += -hx *
00376                       ( x[inode] * y[inode] / ( r1 * ( 1 + x[inode] * x[inode] ) ) +
00377                         4.0 * xmid * ymid / ( rm * ( 1 + xmid * xmid ) ) +
00378                         x[inode2] * y[inode2] / ( r2 * ( 1 + x[inode2] * x[inode2] ) ) ) /
00379                       6.0;
00380 
00381             bary_y += hx * ( 1.0 / r1 + 4.0 / rm + 1.0 / r2 ) / 6.0;
00382         }
00383         bary_x = bary_x / area;
00384         bary_y = bary_y / area;
00385 
00386         // barycenter in Cartesian X,Y,Z coordinates
00387         moab::CartVect barycent;
00388         moab::IntxUtils::reverse_gnomonic_projection( bary_x, bary_y, R, plane, barycent );
00389 
00390         // set barycenter
00391         std::vector< double > barycenter( 3 );
00392         barycenter[0] = barycent[0];
00393         barycenter[1] = barycent[1];
00394         barycenter[2] = barycent[2];
00395         rval          = mb->tag_set_data( barycenterTag, &icell, 1, &barycenter[0] );MB_CHK_ERR( rval );
00396 
00397         // set area
00398         rval = mb->tag_set_data( areaTag, &icell, 1, &area );MB_CHK_ERR( rval );
00399     }
00400 
00401     return moab::MB_SUCCESS;
00402 }
00403 
00404 moab::ErrorCode set_density( moab::Interface* mb, moab::EntityHandle euler_set, moab::Tag& barycenterTag,
00405                              moab::Tag& rhoTag, int field_type )
00406 {
00407     // get cells
00408     moab::Range cells;
00409     moab::ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, cells );MB_CHK_ERR( rval );
00410 
00411     // get barycenters
00412     std::vector< double > cell_barys( 3 * cells.size() );
00413     rval = mb->tag_get_data( barycenterTag, cells, &cell_barys[0] );MB_CHK_ERR( rval );
00414 
00415     // loop over cells
00416     int cell_ind = 0;
00417     for( Range::iterator it = cells.begin(); it != cells.end(); ++it )
00418     {
00419         moab::EntityHandle icell = *it;
00420 
00421         // convert barycenter from 3-D Cartesian to lat/lon
00422         moab::CartVect bary_xyz( cell_barys[cell_ind * 3], cell_barys[cell_ind * 3 + 1], cell_barys[cell_ind * 3 + 2] );
00423         moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( bary_xyz );
00424 
00425         if( field_type == 1 )  // cosine bells
00426         {
00427             //                 lon1,        lat1  lon2    lat2   b    c  hmax  r
00428             double params[] = { 5 * M_PI / 6.0, 0.0, 7 * M_PI / 6, 0.0, 0.1, 0.9, 1., 0.5 };
00429 
00430             double rho_barycent = IntxUtilsCSLAM::quasi_smooth_field( sphCoord.lon, sphCoord.lat, params );
00431             rval                = mb->tag_set_data( rhoTag, &icell, 1, &rho_barycent );MB_CHK_ERR( rval );
00432         }
00433 
00434         if( field_type == 2 )  // Gaussian hills
00435         {
00436             moab::CartVect p1, p2;
00437             moab::IntxUtils::SphereCoords spr;
00438             spr.R   = 1;
00439             spr.lat = M_PI / 3;
00440             spr.lon = M_PI;
00441             p1      = moab::IntxUtils::spherical_to_cart( spr );
00442             spr.lat = -M_PI / 3;
00443             p2      = moab::IntxUtils::spherical_to_cart( spr );
00444             // X1, Y1, Z1, X2, Y2, Z2, ?, ?
00445             double params[] = { p1[0], p1[1], p1[2], p2[0], p2[1], p2[2], 1, 5. };
00446 
00447             double rho_barycent = IntxUtilsCSLAM::smooth_field( sphCoord.lon, sphCoord.lat, params );
00448             rval                = mb->tag_set_data( rhoTag, &icell, 1, &rho_barycent );MB_CHK_ERR( rval );
00449         }
00450 
00451         if( field_type == 3 )  // Zalesak cylinders
00452         {
00453             //                   lon1,      lat1,    lon2,   lat2, b,   c,   r
00454             double params[] = { 5 * M_PI / 6.0, 0.0, 7 * M_PI / 6.0, 0.0, 0.1, 0.9, 0.5 };
00455 
00456             double rho_barycent = IntxUtilsCSLAM::slotted_cylinder_field( sphCoord.lon, sphCoord.lat, params );
00457             rval                = mb->tag_set_data( rhoTag, &icell, 1, &rho_barycent );MB_CHK_ERR( rval );
00458         }
00459 
00460         if( field_type == 4 )  // constant
00461         {
00462             double rho_barycent = 1.0;
00463             rval                = mb->tag_set_data( rhoTag, &icell, 1, &rho_barycent );MB_CHK_ERR( rval );
00464         }
00465 
00466         cell_ind++;
00467     }
00468 
00469     return moab::MB_SUCCESS;
00470 }
00471 
00472 moab::ErrorCode get_linear_reconstruction( moab::Interface* mb, moab::EntityHandle set, moab::Tag& rhoTag,
00473                                            moab::Tag& planeTag, moab::Tag& barycenterTag, moab::Tag& linearCoefTag )
00474 {
00475     // get all entities of dimension 2
00476     Range cells;
00477     ErrorCode rval = mb->get_entities_by_dimension( set, 2, cells );MB_CHK_ERR( rval );
00478 
00479     // Get coefficients for reconstruction (in cubed-sphere coordinates)
00480     for( Range::iterator it = cells.begin(); it != cells.end(); ++it )
00481     {
00482         moab::EntityHandle icell = *it;
00483 
00484         // get the nodes, then the coordinates
00485         const moab::EntityHandle* verts;
00486         int num_nodes;
00487         rval = mb->get_connectivity( icell, verts, num_nodes );MB_CHK_ERR( rval );
00488 
00489         moab::Range adjacentEdges;
00490         rval = mb->get_adjacencies( &icell, 1, 1, true, adjacentEdges );MB_CHK_ERR( rval );
00491 
00492         // get adjacent cells from edges
00493         moab::Range adjacentCells;
00494         rval = mb->get_adjacencies( adjacentEdges, 2, true, adjacentCells, Interface::UNION );MB_CHK_ERR( rval );
00495 
00496         // get gnomonic plane
00497         int plane = 0;
00498         rval      = mb->tag_get_data( planeTag, &icell, 1, &plane );MB_CHK_ERR( rval );
00499 
00500         std::vector< double > dx( adjacentCells.size() - 1 );
00501         std::vector< double > dy( adjacentCells.size() - 1 );
00502         std::vector< double > dr( adjacentCells.size() - 1 );
00503         double bary_x;
00504         double bary_y;
00505 
00506         // get barycenter of cell where reconstruction occurs
00507         double rad = 1;
00508         std::vector< double > barycent( 3 );
00509         rval = mb->tag_get_data( barycenterTag, &icell, 1, &barycent[0] );MB_CHK_ERR( rval );
00510         CartVect barycenter( barycent[0], barycent[1], barycent[2] );
00511         double cellbaryx = 0;
00512         double cellbaryy = 0;
00513         moab::IntxUtils::gnomonic_projection( barycenter, rad, plane, cellbaryx, cellbaryy );
00514 
00515         // get density value
00516         double rhocell = 0;
00517         rval           = mb->tag_get_data( rhoTag, &icell, 1, &rhocell );MB_CHK_ERR( rval );
00518 
00519         // get barycenters of surrounding cells
00520         std::vector< double > cell_barys( 3 * adjacentCells.size() );
00521         rval = mb->tag_get_data( barycenterTag, adjacentCells, &cell_barys[0] );MB_CHK_ERR( rval );
00522 
00523         // get density of surrounding cells
00524         std::vector< double > rho_vals( adjacentCells.size() );
00525         rval = mb->tag_get_data( rhoTag, adjacentCells, &rho_vals[0] );MB_CHK_ERR( rval );
00526 
00527         std::size_t jind = 0;
00528         for( std::size_t i = 0; i < adjacentCells.size(); i++ )
00529         {
00530 
00531             if( adjacentCells[i] != icell )
00532             {
00533 
00534                 CartVect bary_xyz( cell_barys[i * 3], cell_barys[i * 3 + 1], cell_barys[i * 3 + 2] );
00535                 moab::IntxUtils::gnomonic_projection( bary_xyz, rad, plane, bary_x, bary_y );
00536 
00537                 dx[jind] = bary_x - cellbaryx;
00538                 dy[jind] = bary_y - cellbaryy;
00539                 dr[jind] = rho_vals[i] - rhocell;
00540 
00541                 jind++;
00542             }
00543         }
00544 
00545         std::vector< double > linearCoef( 3 );
00546         if( adjacentCells.size() == 5 )
00547         {
00548 
00549             // compute normal equations matrix
00550             double N11 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2] + dx[3] * dx[3];
00551             double N22 = dy[0] * dy[0] + dy[1] * dy[1] + dy[2] * dy[2] + dy[3] * dy[3];
00552             double N12 = dx[0] * dy[0] + dx[1] * dy[1] + dx[2] * dy[2] + dx[3] * dy[3];
00553 
00554             // rhs
00555             double Rx = dx[0] * dr[0] + dx[1] * dr[1] + dx[2] * dr[2] + dx[3] * dr[3];
00556             double Ry = dy[0] * dr[0] + dy[1] * dr[1] + dy[2] * dr[2] + dy[3] * dr[3];
00557 
00558             // determinant
00559             double Det = N11 * N22 - N12 * N12;
00560 
00561             // solution
00562             linearCoef[0] = ( Rx * N22 - Ry * N12 ) / Det;
00563             linearCoef[1] = ( Ry * N11 - Rx * N12 ) / Det;
00564             linearCoef[2] = rhocell - linearCoef[0] * cellbaryx - linearCoef[1] * cellbaryy;
00565         }
00566         else
00567         {
00568             // default to first order
00569             linearCoef[0] = 0.0;
00570             linearCoef[1] = 0.0;
00571             linearCoef[2] = rhocell;
00572             std::cout << "Need 4 adjacent cells for linear reconstruction! \n";
00573         }
00574 
00575         rval = mb->tag_set_data( linearCoefTag, &icell, 1, &linearCoef[0] );MB_CHK_ERR( rval );
00576     }
00577 
00578     return moab::MB_SUCCESS;
00579 }
00580 
00581 moab::ErrorCode get_departure_grid( moab::Interface* mb, moab::EntityHandle euler_set, moab::EntityHandle lagr_set,
00582                                     moab::EntityHandle covering_set, int tStep, Range& connecVerts )
00583 {
00584     EntityHandle dum = 0;
00585     Tag corrTag;
00586     mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE, &dum );
00587 
00588     double t       = tStep * T / numSteps;  // numSteps is global; so is T
00589     double delta_t = T / numSteps;          // this is global too, actually
00590     // double delta_t = 0.0001;
00591     // double t = delta_t;
00592 
00593     Range polys;
00594     ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, polys );MB_CHK_ERR( rval );
00595 
00596     // change coordinates of lagr mesh vertices
00597     for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
00598     {
00599         moab::EntityHandle oldV = *vit;
00600         CartVect posi;
00601         rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
00602 
00603         // Intx utils, case 1
00604         CartVect newPos;
00605         departure_point_swirl_rot( posi, t, delta_t, newPos );
00606         newPos = radius * newPos;  // do we need this? the radius should be 1
00607         moab::EntityHandle new_vert;
00608         rval = mb->tag_get_data( corrTag, &oldV, 1, &new_vert );MB_CHK_ERR( rval );
00609 
00610         // set the new position for the new vertex
00611         rval = mb->set_coords( &new_vert, 1, &( newPos[0] ) );MB_CHK_ERR( rval );
00612     }
00613 
00614     // if in parallel, we have to move some elements to another proc, and receive other cells
00615     // from other procs
00616     rval = pworker->create_departure_mesh_3rd_alg( lagr_set, covering_set );MB_CHK_ERR( rval );
00617 
00618     return rval;
00619 }
00620 
00621 // !!! For now serial !!!
00622 moab::ErrorCode update_density( moab::Interface* mb, moab::EntityHandle euler_set, moab::EntityHandle lagr_set,
00623                                 moab::EntityHandle out_set, moab::Tag& rhoTag, moab::Tag& areaTag,
00624                                 moab::Tag& rhoCoefsTag, moab::Tag& weightsTag, moab::Tag& planeTag )
00625 {
00626     //  moab::ParallelComm * parcomm = ParallelComm::get_pcomm(mb, 0);
00627     ErrorCode rval;
00628     double R               = 1.0;
00629     moab::EntityHandle dum = 0;
00630     Tag corrTag;
00631 
00632     rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE, &dum );MB_CHK_ERR( rval );
00633 
00634     // get all polygons out of out_set; then see where are they coming from
00635     moab::Range polys;
00636     rval = mb->get_entities_by_dimension( out_set, 2, polys );MB_CHK_ERR( rval );
00637 
00638     // get all Lagrangian cells
00639     moab::Range rs1;
00640     rval = mb->get_entities_by_dimension( lagr_set, 2, rs1 );MB_CHK_ERR( rval );
00641 
00642     // get all Eulerian cells
00643     moab::Range rs2;
00644     rval = mb->get_entities_by_dimension( euler_set, 2, rs2 );MB_CHK_ERR( rval );
00645 
00646     // get gnomonic plane for Eulerian cells
00647     std::vector< int > plane( rs2.size() );
00648     rval = mb->tag_get_data( planeTag, rs2, &plane[0] );MB_CHK_ERR( rval );
00649 
00650     // get Eulerian cell reconstruction coefficients
00651     std::vector< double > rhoCoefs( 3 * rs2.size() );
00652     rval = mb->tag_get_data( rhoCoefsTag, rs2, &rhoCoefs[0] );MB_CHK_ERR( rval );
00653 
00654     // get intersection weights
00655     std::vector< double > weights( 3 * polys.size() );
00656     rval = mb->tag_get_data( weightsTag, polys, &weights[0] );MB_CHK_ERR( rval );
00657 
00658     // Initialize the new values
00659     std::vector< double > newValues( rs2.size(), 0. );  // initialize with 0 all of them
00660 
00661     // For each polygon get red/blue parent
00662     moab::Tag tgtParentTag;
00663     moab::Tag srcParentTag;
00664     rval = mb->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag, MB_TAG_DENSE );MB_CHK_ERR( rval );
00665     rval = mb->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag, MB_TAG_DENSE );MB_CHK_ERR( rval );
00666 
00667     // mass_lagr = (\sum_intx \int rho^h(x,y) dV)
00668     // rho_eul^n+1 = mass_lagr/area_eul
00669     double check_intx_area = 0.;
00670     int polyIndex          = 0;
00671     for( Range::iterator it = polys.begin(); it != polys.end(); ++it )
00672     {
00673         moab::EntityHandle poly = *it;
00674         int blueIndex, redIndex;
00675         rval = mb->tag_get_data( srcParentTag, &poly, 1, &blueIndex );MB_CHK_ERR( rval );
00676 
00677         moab::EntityHandle blue = rs1[blueIndex];
00678         rval                    = mb->tag_get_data( tgtParentTag, &poly, 1, &redIndex );MB_CHK_ERR( rval );
00679 
00680         moab::EntityHandle redArr;
00681         rval = mb->tag_get_data( corrTag, &blue, 1, &redArr );MB_CHK_ERR( rval );
00682         int arrRedIndex = rs2.index( redArr );
00683 
00684         // sum into new density values
00685         newValues[arrRedIndex] += rhoCoefs[redIndex * 3] * weights[polyIndex * 3] +
00686                                   rhoCoefs[redIndex * 3 + 1] * weights[polyIndex * 3 + 1] +
00687                                   rhoCoefs[redIndex * 3 + 2] * weights[polyIndex * 3 + 2];
00688 
00689         check_intx_area += weights[polyIndex * 3 + 2];
00690 
00691         polyIndex++;
00692     }
00693 
00694     // now divide by red area (current)
00695     int j                   = 0;
00696     Range::iterator iter    = rs2.begin();
00697     void* data              = NULL;  // used for stored area
00698     int count               = 0;
00699     double total_mass_local = 0.;
00700     while( iter != rs2.end() )
00701     {
00702         rval = mb->tag_iterate( areaTag, iter, rs2.end(), count, data );MB_CHK_ERR( rval );
00703 
00704         double* ptrArea = (double*)data;
00705         for( int i = 0; i < count; i++, ++iter, j++, ptrArea++ )
00706         {
00707             total_mass_local += newValues[j];
00708             newValues[j] /= ( *ptrArea );
00709         }
00710     }
00711 
00712     rval = mb->tag_set_data( rhoTag, rs2, &newValues[0] );MB_CHK_ERR( rval );
00713 
00714     std::cout << "total mass now:" << total_mass_local << "\n";
00715     std::cout << "check: total intersection area: (4 * M_PI * R^2): " << 4 * M_PI * R * R << " " << check_intx_area
00716               << "\n";
00717 
00718     return MB_SUCCESS;
00719 }
00720 
00721 /*
00722  *  Deformational flow
00723  */
00724 void departure_point_swirl( moab::CartVect& arrival_point, double t, double delta_t, moab::CartVect& departure_point )
00725 {
00726 
00727     // always assume radius is 1 here?
00728     moab::IntxUtils::SphereCoords sph = moab::IntxUtils::cart_to_spherical( arrival_point );
00729     double k                          = 2.4;  // flow parameter
00730     /*     radius needs to be within some range   */
00731     double sl2      = sin( sph.lon / 2 );
00732     double pit      = M_PI * t / T;
00733     double omega    = M_PI / T;
00734     double costheta = cos( sph.lat );
00735     // double u = k * sl2*sl2 * sin(2*sph.lat) * cos(pit);
00736     double v = k * sin( sph.lon ) * costheta * cos( pit );
00737     // double psi = k * sl2 * sl2 *costheta * costheta * cos(pit);
00738     double u_tilda = 2 * k * sl2 * sl2 * sin( sph.lat ) * cos( pit );
00739 
00740     // formula 35, page 8
00741     // this will approximate dep point using a Taylor series with up to second derivative
00742     // this will be O(delta_t^3) exact.
00743     double lon_dep =
00744         sph.lon - delta_t * u_tilda -
00745         delta_t * delta_t * k * sl2 *
00746             ( sl2 * sin( sph.lat ) * sin( pit ) * omega - u_tilda * sin( sph.lat ) * cos( pit ) * cos( sph.lon / 2 ) -
00747               v * sl2 * costheta * cos( pit ) );
00748     // formula 36, page 8 again
00749     double lat_dep = sph.lat - delta_t * v -
00750                      delta_t * delta_t / 4 * k *
00751                          ( sin( sph.lon ) * cos( sph.lat ) * sin( pit ) * omega -
00752                            u_tilda * cos( sph.lon ) * cos( sph.lat ) * cos( pit ) +
00753                            v * sin( sph.lon ) * sin( sph.lat ) * cos( pit ) );
00754     moab::IntxUtils::SphereCoords sph_dep;
00755     sph_dep.R   = 1.;  // radius
00756     sph_dep.lat = lat_dep;
00757     sph_dep.lon = lon_dep;
00758 
00759     departure_point = moab::IntxUtils::spherical_to_cart( sph_dep );
00760     return;
00761 }
00762 
00763 /*
00764  *  Deformational flow with rotation
00765  */
00766 void departure_point_swirl_rot( moab::CartVect& arrival_point, double t, double delta_t,
00767                                 moab::CartVect& departure_point )
00768 {
00769 
00770     moab::IntxUtils::SphereCoords sph = moab::IntxUtils::cart_to_spherical( arrival_point );
00771     double omega                      = M_PI / T;
00772     double gt                         = cos( M_PI * t / T );
00773 
00774     double lambda  = sph.lon - 2.0 * omega * t;
00775     double u_tilda = 4.0 * sin( lambda ) * sin( lambda ) * sin( sph.lat ) * gt + 2.0 * omega;
00776     double v       = 2.0 * sin( 2.0 * lambda ) * cos( sph.lat ) * gt;
00777 
00778     double lon_dep = sph.lon - delta_t * u_tilda -
00779                      delta_t * delta_t * 2.0 * sin( lambda ) *
00780                          ( sin( lambda ) * sin( sph.lat ) * sin( omega * t ) * omega -
00781                            sin( lambda ) * cos( sph.lat ) * cos( omega * t ) * v -
00782                            2.0 * cos( lambda ) * sin( sph.lat ) * cos( omega * t ) * u_tilda );
00783 
00784     double lat_dep = sph.lat - delta_t * v -
00785                      delta_t * delta_t * 2.0 *
00786                          ( cos( sph.lat ) * sin( omega * t ) * omega * sin( lambda ) * cos( lambda ) -
00787                            2.0 * u_tilda * cos( sph.lat ) * cos( omega * t ) * cos( lambda ) * cos( lambda ) +
00788                            u_tilda * cos( sph.lat ) * cos( omega * t ) +
00789                            v * sin( sph.lat ) * cos( omega * t ) * sin( lambda ) * cos( lambda ) );
00790 
00791     moab::IntxUtils::SphereCoords sph_dep;
00792     sph_dep.R   = 1.;  // radius
00793     sph_dep.lat = lat_dep;
00794     sph_dep.lon = lon_dep;
00795 
00796     departure_point = moab::IntxUtils::spherical_to_cart( sph_dep );
00797     return;
00798 }
00799 
00800 /*
00801  *  Zonal flow
00802  */
00803 moab::ErrorCode get_intersection_weights( moab::Interface* mb, moab::EntityHandle euler_set,
00804                                           moab::EntityHandle lagr_set, moab::EntityHandle intx_set, moab::Tag& planeTag,
00805                                           moab::Tag& weightsTag )
00806 {
00807     // get all intersection polygons
00808     moab::Range polys;
00809     ErrorCode rval = mb->get_entities_by_dimension( intx_set, 2, polys );MB_CHK_ERR( rval );
00810 
00811     // get all Eulerian cells
00812     moab::Range eul_cells;
00813     rval = mb->get_entities_by_dimension( euler_set, 2, eul_cells );MB_CHK_ERR( rval );
00814 
00815     // get all Lagrangian cells
00816     moab::Range lagr_cells;
00817     rval = mb->get_entities_by_dimension( lagr_set, 2, lagr_cells );MB_CHK_ERR( rval );
00818 
00819     // get tag for Eulerian parent cell of intersection polygon
00820     moab::Tag tgtParentTag;
00821     rval = mb->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag, MB_TAG_DENSE );
00822 
00823     // get tag for Lagrangian parent cell of intersection polygon
00824     moab::Tag srcParentTag;
00825     rval = mb->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag, MB_TAG_DENSE );MB_CHK_ERR( rval );
00826 
00827     // get gnomonic plane for Eulerian cells
00828     std::vector< int > plane( eul_cells.size() );
00829     rval = mb->tag_get_data( planeTag, eul_cells, &plane[0] );MB_CHK_ERR( rval );
00830 
00831     double total_area = 0.;
00832     for( moab::Range::iterator it = polys.begin(); it != polys.end(); ++it )
00833     {
00834         moab::EntityHandle poly = *it;
00835 
00836         // get the nodes
00837         const moab::EntityHandle* verts;
00838         int num_nodes;
00839         rval = mb->get_connectivity( poly, verts, num_nodes );MB_CHK_ERR( rval );
00840 
00841         // get coordinates
00842         std::vector< double > coords( 3 * num_nodes );
00843         rval = mb->get_coords( verts, num_nodes, &coords[0] );MB_CHK_ERR( rval );
00844 
00845         // get index of Eulerian parent cell for polygon
00846         int redIndex;
00847         rval = mb->tag_get_data( tgtParentTag, &poly, 1, &redIndex );MB_CHK_ERR( rval );
00848 
00849         // get index of Lagrangian parent cell for polygon
00850         int blueIndex;
00851         rval = mb->tag_get_data( srcParentTag, &poly, 1, &blueIndex );MB_CHK_ERR( rval );
00852 
00853         std::vector< double > x( num_nodes );
00854         std::vector< double > y( num_nodes );
00855         double poly_area = 0;
00856         double poly_intx = 0;
00857         double poly_inty = 0;
00858         double R         = 1.0;
00859         for( int inode = 0; inode < num_nodes; inode++ )
00860         {
00861             double rad = sqrt( coords[inode * 3] * coords[inode * 3] + coords[inode * 3 + 1] * coords[inode * 3 + 1] +
00862                                coords[inode * 3 + 2] * coords[inode * 3 + 2] );
00863             moab::CartVect xyzcoord( coords[inode * 3] / rad, coords[inode * 3 + 1] / rad,
00864                                      coords[inode * 3 + 2] / rad );
00865             moab::IntxUtils::gnomonic_projection( xyzcoord, R, plane[redIndex], x[inode], y[inode] );
00866         }
00867 
00868         std::vector< double > weights( 3 );
00869         for( int inode = 0; inode < num_nodes; inode++ )
00870         {
00871             int inode2 = inode + 1;
00872             if( inode2 >= num_nodes ) inode2 = 0;
00873             double xmid = 0.5 * ( x[inode] + x[inode2] );
00874             double ymid = 0.5 * ( y[inode] + y[inode2] );
00875             double r1   = sqrt( 1 + x[inode] * x[inode] + y[inode] * y[inode] );
00876             double rm   = sqrt( 1 + xmid * xmid + ymid * ymid );
00877             double r2   = sqrt( 1 + x[inode2] * x[inode2] + y[inode2] * y[inode2] );
00878             double hx   = x[inode2] - x[inode];
00879 
00880             poly_area += -hx *
00881                          ( y[inode] / ( r1 * ( 1 + x[inode] * x[inode] ) ) + 4.0 * ymid / ( rm * ( 1 + xmid * xmid ) ) +
00882                            y[inode2] / ( r2 * ( 1 + x[inode2] * x[inode2] ) ) ) /
00883                          6.0;
00884 
00885             poly_intx += -hx *
00886                          ( x[inode] * y[inode] / ( r1 * ( 1 + x[inode] * x[inode] ) ) +
00887                            4.0 * xmid * ymid / ( rm * ( 1 + xmid * xmid ) ) +
00888                            x[inode2] * y[inode2] / ( r2 * ( 1 + x[inode2] * x[inode2] ) ) ) /
00889                          6.0;
00890 
00891             poly_inty += hx * ( 1.0 / r1 + 4.0 / rm + 1.0 / r2 ) / 6.0;
00892         }
00893         weights[0] = poly_intx;
00894         weights[1] = poly_inty;
00895         weights[2] = poly_area;
00896 
00897         total_area += poly_area;
00898 
00899         rval = mb->tag_set_data( weightsTag, &poly, 1, &weights[0] );MB_CHK_ERR( rval );
00900     }
00901 
00902     std::cout << "polygon area = " << total_area << "\n";
00903 
00904     return moab::MB_SUCCESS;
00905 }
00906 
00907 ErrorCode create_lagr_mesh( moab::Interface* mb, moab::EntityHandle euler_set, moab::EntityHandle lagr_set )
00908 {
00909     // create the handle tag for the corresponding element / vertex
00910 
00911     moab::EntityHandle dum = 0;
00912 
00913     moab::Tag corrTag;
00914     ErrorCode rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );
00915 
00916     MB_CHK_ERR( rval );
00917     // give the same global id to new verts and cells created in the lagr(departure) mesh
00918     moab::Tag gid = mb->globalId_tag();
00919 
00920     moab::Range polys;
00921     rval = mb->get_entities_by_dimension( euler_set, 2, polys );MB_CHK_ERR( rval );
00922 
00923     moab::Range connecVerts;
00924     rval = mb->get_connectivity( polys, connecVerts );MB_CHK_ERR( rval );
00925 
00926     std::map< moab::EntityHandle, moab::EntityHandle > newNodes;
00927     for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
00928     {
00929         moab::EntityHandle oldV = *vit;
00930         moab::CartVect posi;
00931         rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
00932         int global_id;
00933         rval = mb->tag_get_data( gid, &oldV, 1, &global_id );MB_CHK_ERR( rval );
00934         moab::EntityHandle new_vert;
00935         // duplicate the position
00936         rval = mb->create_vertex( &( posi[0] ), new_vert );MB_CHK_ERR( rval );
00937         newNodes[oldV] = new_vert;
00938         // set also the correspondent tag :)
00939         rval = mb->tag_set_data( corrTag, &oldV, 1, &new_vert );MB_CHK_ERR( rval );
00940 
00941         // also the other side
00942         // need to check if we really need this; the new vertex will never need the old vertex
00943         // we have the global id which is the same
00944         /*rval = mb->tag_set_data(corrTag, &new_vert, 1, &oldV);MB_CHK_ERR(rval);*/
00945         // set the global id on the corresponding vertex the same as the initial vertex
00946         rval = mb->tag_set_data( gid, &new_vert, 1, &global_id );MB_CHK_ERR( rval );
00947     }
00948 
00949     for( Range::iterator it = polys.begin(); it != polys.end(); ++it )
00950     {
00951         moab::EntityHandle q = *it;
00952         int global_id;
00953         int nnodes;
00954         const moab::EntityHandle* conn;
00955 
00956         rval = mb->get_connectivity( q, conn, nnodes );MB_CHK_ERR( rval );
00957         rval = mb->tag_get_data( gid, &q, 1, &global_id );MB_CHK_ERR( rval );
00958 
00959         EntityType typeElem = mb->type_from_handle( q );
00960         std::vector< moab::EntityHandle > new_conn( nnodes );
00961         for( int i = 0; i < nnodes; i++ )
00962         {
00963             moab::EntityHandle v1 = conn[i];
00964             new_conn[i]           = newNodes[v1];
00965         }
00966 
00967         moab::EntityHandle newElement;
00968         rval = mb->create_element( typeElem, &new_conn[0], nnodes, newElement );MB_CHK_ERR( rval );
00969         // set the corresponding tag; not sure we need this one, from old to new
00970         /*rval = mb->tag_set_data(corrTag, &q, 1, &newElement);MB_CHK_ERR(rval);*/
00971         rval = mb->tag_set_data( corrTag, &newElement, 1, &q );MB_CHK_ERR( rval );
00972         // set the global id
00973         rval = mb->tag_set_data( gid, &newElement, 1, &global_id );MB_CHK_ERR( rval );
00974 
00975         rval = mb->add_entities( lagr_set, &newElement, 1 );MB_CHK_ERR( rval );
00976     }
00977 
00978     return MB_SUCCESS;
00979 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines