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