Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
linear_remap.cpp
Go to the documentation of this file.
00001 #include <iostream>
00002 #include <cmath>
00003 
00004 #include "moab/Core.hpp"
00005 #include "moab/Interface.hpp"
00006 #include "moab/IntxMesh/Intx2MeshOnSphere.hpp"
00007 #include "moab/ParallelComm.hpp"
00008 #include "moab/ProgOptions.hpp"
00009 #include "MBParallelConventions.h"
00010 #include "moab/ReadUtilIface.hpp"
00011 #include "MBTagConventions.hpp"
00012 
00013 #include "moab/IntxMesh/IntxUtils.hpp"
00014 #include "IntxUtilsCSLAM.hpp"
00015 
00016 #include "TestUtil.hpp"
00017 
00018 // std::string file_name("./uniform_15.g");
00019 // std::string file_name("./eulerHomme.vtk");
00020 
00021 using namespace moab;
00022 
00023 // computes cell barycenters in Cartesian X,Y,Z coordinates
00024 void get_barycenters( moab::Interface* mb, moab::EntityHandle set, moab::Tag& planeTag, moab::Tag& barycenterTag );
00025 
00026 // computes gnomonic plane for each cell of the Eulerian mesh
00027 void get_gnomonic_plane( moab::Interface* mb, moab::EntityHandle set, moab::Tag& planeTag );
00028 
00029 // computes coefficients (A,B,C) for linear reconstruction in gnomonic coordinates: rho(x,y) = Ax +
00030 // By + C
00031 void get_linear_reconstruction( moab::Interface* mb,
00032                                 moab::EntityHandle set,
00033                                 moab::Tag& rhoTag,
00034                                 moab::Tag& planeTag,
00035                                 moab::Tag& barycenterTag,
00036                                 moab::Tag& linearCoefTag );
00037 
00038 // evaluates the integral of rho(x,y) over cell, should be equal to cell average rho
00039 void test_linear_reconstruction( moab::Interface* mb,
00040                                  moab::EntityHandle set,
00041                                  moab::Tag& rhoTag,
00042                                  moab::Tag& planeTag,
00043                                  moab::Tag& barycenterTag,
00044                                  moab::Tag& linearCoefTag );
00045 
00046 // set density function
00047 moab::ErrorCode add_field_value( moab::Interface* mb,
00048                                  moab::EntityHandle euler_set,
00049                                  int rank,
00050                                  moab::Tag& tagTracer,
00051                                  moab::Tag& tagElem,
00052                                  moab::Tag& tagArea,
00053                                  int field_type );
00054 
00055 // functions that implement gnomonic projection as in Homme (different from Moab implementation)
00056 int gnomonic_projection_test( const moab::CartVect& pos, double R, int plane, double& c1, double& c2 );
00057 int reverse_gnomonic_projection_test( const double& c1, const double& c2, double R, int plane, moab::CartVect& pos );
00058 void decide_gnomonic_plane_test( const CartVect& pos, int& plane );
00059 
00060 double radius = 1.;  // in m:  6371220.
00061 
00062 int main( int argc, char* argv[] )
00063 {
00064 
00065     // set up MOAB interface and parallel communication
00066     MPI_Init( &argc, &argv );
00067     moab::Core moab;
00068     moab::Interface& mb = moab;
00069     moab::ParallelComm mb_pcomm( &mb, MPI_COMM_WORLD );
00070 
00071     // int rank = mb_pcomm->proc_config().proc_rank();
00072     int rank = mb_pcomm.proc_config().proc_rank();
00073 
00074     // create meshset
00075     moab::EntityHandle euler_set;
00076     moab::ErrorCode rval = mb.create_meshset( MESHSET_SET, euler_set );CHECK_ERR( rval );
00077 
00078     // std::stringstream opts;
00079     // opts <<
00080     // "PARALLEL=READ_PART;PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;GATHER_SET=0;PARTITION_METHOD=TRIVIAL_PARTITION;VARIABLE=";
00081     // opts << "PARALLEL=READ_PART;PARTITION;PARALLEL_RESOLVE_SHARED_ENTS";
00082     std::string fileN = TestDir + "unittest/mbcslam/fine4.h5m";
00083 
00084     rval = mb.load_file( fileN.c_str(), &euler_set );CHECK_ERR( rval );
00085 
00086     // Create tag for cell density
00087     moab::Tag rhoTag = 0;
00088     rval = mb.tag_get_handle( "Density", 1, moab::MB_TYPE_DOUBLE, rhoTag, moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );CHECK_ERR( rval );
00089     moab::Tag rhoNodeTag = 0;
00090     rval                 = mb.tag_get_handle( "DensityNode", 1, moab::MB_TYPE_DOUBLE, rhoNodeTag,
00091                                               moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );CHECK_ERR( rval );
00092     // Create tag for cell area
00093     moab::Tag areaTag = 0;
00094     rval = mb.tag_get_handle( "Area", 1, moab::MB_TYPE_DOUBLE, areaTag, moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );CHECK_ERR( rval );
00095     // Create tag for cell barycenters in 3D Cartesian space
00096     moab::Tag barycenterTag = 0;
00097     rval                    = mb.tag_get_handle( "CellBarycenter", 3, moab::MB_TYPE_DOUBLE, barycenterTag,
00098                                                  moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );CHECK_ERR( rval );
00099     // Create tag for cell density reconstruction coefficients
00100     moab::Tag coefRhoTag = 0;
00101     rval                 = mb.tag_get_handle( "LinearCoefRho", 3, moab::MB_TYPE_DOUBLE, coefRhoTag,
00102                                               moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );CHECK_ERR( rval );
00103     // Create tag for index of gnomonic plane for each cell
00104     moab::Tag planeTag = 0;
00105     rval               = mb.tag_get_handle( "gnomonicPlane", 1, moab::MB_TYPE_INTEGER, planeTag,
00106                                             moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );CHECK_ERR( rval );
00107     // Set density distributions
00108     rval = add_field_value( &mb, euler_set, rank, rhoNodeTag, rhoTag, areaTag, 2 );CHECK_ERR( rval );
00109     // get cell plane
00110     get_gnomonic_plane( &mb, euler_set, planeTag );
00111 
00112     // get cell barycenters
00113     get_barycenters( &mb, euler_set, planeTag, barycenterTag );
00114 
00115     // get linear reconstruction
00116     get_linear_reconstruction( &mb, euler_set, rhoTag, planeTag, barycenterTag, coefRhoTag );
00117 
00118     // test linear reconstruction
00119     test_linear_reconstruction( &mb, euler_set, rhoTag, planeTag, barycenterTag, coefRhoTag );
00120 
00121     MPI_Finalize();
00122     return 0;
00123 }
00124 
00125 void get_gnomonic_plane( moab::Interface* mb, moab::EntityHandle set, moab::Tag& planeTag )
00126 {
00127     // get all entities of dimension 2
00128     // then get the connectivity, etc
00129     moab::Range cells;
00130     ErrorCode rval = mb->get_entities_by_dimension( set, 2, cells );
00131     if( MB_SUCCESS != rval ) return;
00132 
00133     for( Range::iterator it = cells.begin(); it != cells.end(); ++it )
00134     {
00135         moab::EntityHandle icell = *it;
00136         // get the nodes, then the coordinates
00137         const moab::EntityHandle* verts;
00138         int num_nodes;
00139         rval = mb->get_connectivity( icell, verts, num_nodes );
00140         if( MB_SUCCESS != rval ) return;
00141 
00142         std::vector< double > coords( 3 * num_nodes );
00143         // get coordinates
00144         rval = mb->get_coords( verts, num_nodes, &coords[0] );
00145         if( MB_SUCCESS != rval ) return;
00146 
00147         double centerx = 0;
00148         double centery = 0;
00149         double centerz = 0;
00150         for( int inode = 0; inode < num_nodes; inode++ )
00151         {
00152             centerx += coords[inode * 3] / num_nodes;
00153             centery += coords[inode * 3 + 1] / num_nodes;
00154             centerz += coords[inode * 3 + 2] / num_nodes;
00155         }
00156         double R = std::sqrt( centerx * centerx + centery * centery + centerz * centerz );
00157         centerx  = centerx / R;
00158         centery  = centery / R;
00159         centerz  = centerz / R;
00160         moab::CartVect center( centerx, centery, centerz );
00161         R = 1.0;
00162 
00163         int plane = 0;
00164         moab::IntxUtils::decide_gnomonic_plane( center, plane );
00165         // decide_gnomonic_plane_test(center,plane);
00166 
00167         rval = mb->tag_set_data( planeTag, &icell, 1, &plane );CHECK_ERR( rval );
00168     }
00169     return;
00170 }
00171 
00172 void get_barycenters( moab::Interface* mb, moab::EntityHandle set, moab::Tag& planeTag, moab::Tag& barycenterTag )
00173 {
00174     // get all entities of dimension 2
00175     moab::Range cells;
00176     ErrorCode rval = mb->get_entities_by_dimension( set, 2, cells );
00177     if( MB_SUCCESS != rval ) return;
00178 
00179     // set sphere radius to 1
00180     double R = 1.0;
00181 
00182     for( Range::iterator it = cells.begin(); it != cells.end(); ++it )
00183     {
00184         moab::EntityHandle icell = *it;
00185 
00186         // get the nodes
00187         const moab::EntityHandle* verts;
00188         int num_nodes;
00189         rval = mb->get_connectivity( icell, verts, num_nodes );
00190         if( MB_SUCCESS != rval ) return;
00191 
00192         // get coordinates
00193         std::vector< double > coords( 3 * num_nodes );
00194         rval = mb->get_coords( verts, num_nodes, &coords[0] );
00195         if( MB_SUCCESS != rval ) return;
00196 
00197         // get plane for cell
00198         int plane = 0;
00199         rval      = mb->tag_get_data( planeTag, &icell, 1, &plane );
00200         if( MB_SUCCESS != rval ) return;
00201         std::vector< double > x( num_nodes );
00202         std::vector< double > y( num_nodes );
00203         double area   = 0;
00204         double bary_x = 0;
00205         double bary_y = 0;
00206         for( int inode = 0; inode < num_nodes; inode++ )
00207         {
00208             // radius should be 1.0, but divide by it just in case for now
00209             double rad = sqrt( coords[inode * 3] * coords[inode * 3] + coords[inode * 3 + 1] * coords[inode * 3 + 1] +
00210                                coords[inode * 3 + 2] * coords[inode * 3 + 2] );
00211             CartVect xyzcoord( coords[inode * 3] / rad, coords[inode * 3 + 1] / rad, coords[inode * 3 + 2] / rad );
00212             moab::IntxUtils::gnomonic_projection( xyzcoord, R, plane, x[inode], y[inode] );
00213             // int dum = gnomonic_projection_test(xyzcoord, R, plane, x[inode],y[inode]);
00214         }
00215 
00216         for( int inode = 0; inode < num_nodes; inode++ )
00217         {
00218             int inode2 = inode + 1;
00219             if( inode2 >= num_nodes ) inode2 = 0;
00220             double xmid = 0.5 * ( x[inode] + x[inode2] );
00221             double ymid = 0.5 * ( y[inode] + y[inode2] );
00222             double r1   = sqrt( 1 + x[inode] * x[inode] + y[inode] * y[inode] );
00223             double rm   = sqrt( 1 + xmid * xmid + ymid * ymid );
00224             double r2   = sqrt( 1 + x[inode2] * x[inode2] + y[inode2] * y[inode2] );
00225             double hx   = x[inode2] - x[inode];
00226 
00227             area += hx *
00228                     ( y[inode] / ( r1 * ( 1 + x[inode] * x[inode] ) ) + 4.0 * ymid / ( rm * ( 1 + xmid * xmid ) ) +
00229                       y[inode2] / ( r2 * ( 1 + x[inode2] * x[inode2] ) ) ) /
00230                     6.0;
00231 
00232             bary_x += hx *
00233                       ( x[inode] * y[inode] / ( r1 * ( 1 + x[inode] * x[inode] ) ) +
00234                         4.0 * xmid * ymid / ( rm * ( 1 + xmid * xmid ) ) +
00235                         x[inode2] * y[inode2] / ( r2 * ( 1 + x[inode2] * x[inode2] ) ) ) /
00236                       6.0;
00237 
00238             bary_y += -hx * ( 1.0 / r1 + 4.0 / rm + 1.0 / r2 ) / 6.0;
00239         }
00240 
00241         bary_x = bary_x / area;
00242         bary_y = bary_y / area;
00243 
00244         moab::CartVect barycent;
00245         moab::IntxUtils::reverse_gnomonic_projection( bary_x, bary_y, R, plane, barycent );
00246         // reverse_gnomonic_projection_test(bary_x, bary_y, R, plane, barycent);
00247         std::vector< double > barycenter( 3 );
00248         barycenter[0] = barycent[0];
00249         barycenter[1] = barycent[1];
00250         barycenter[2] = barycent[2];
00251 
00252         rval = mb->tag_set_data( barycenterTag, &icell, 1, &barycenter[0] );CHECK_ERR( rval );
00253     }
00254     return;
00255 }
00256 
00257 void get_linear_reconstruction( moab::Interface* mb,
00258                                 moab::EntityHandle set,
00259                                 moab::Tag& rhoTag,
00260                                 moab::Tag& planeTag,
00261                                 moab::Tag& barycenterTag,
00262                                 moab::Tag& linearCoefTag )
00263 {
00264     // get all entities of dimension 2
00265     Range cells;
00266     ErrorCode rval = mb->get_entities_by_dimension( set, 2, cells );
00267     if( MB_SUCCESS != rval ) return;
00268 
00269     // set sphere radius to 1
00270     double R = 1;
00271 
00272     // Get coefficients of reconstruction (in cubed-sphere coordinates)
00273     // rho(x,y) = Ax + By + C
00274     for( Range::iterator it = cells.begin(); it != cells.end(); ++it )
00275     {
00276         moab::EntityHandle icell = *it;
00277 
00278         // get the nodes, then the coordinates
00279         const moab::EntityHandle* verts;
00280         int num_nodes;
00281         rval = mb->get_connectivity( icell, verts, num_nodes );
00282         if( MB_SUCCESS != rval ) return;
00283 
00284         moab::Range adjacentEdges;
00285         rval = mb->get_adjacencies( &icell, 1, 1, true, adjacentEdges );CHECK_ERR( rval );
00286         // int num_edges = adjacentEdges.size();
00287 
00288         // get adjacent cells from edges
00289         moab::Range adjacentCells;
00290         rval = mb->get_adjacencies( adjacentEdges, 2, true, adjacentCells, Interface::UNION );
00291         if( MB_SUCCESS != rval ) return;
00292         // get plane for cell
00293         int plane = 0;
00294         rval      = mb->tag_get_data( planeTag, &icell, 1, &plane );CHECK_ERR( rval );
00295 
00296         std::vector< double > dx( adjacentCells.size() - 1 );
00297         std::vector< double > dy( adjacentCells.size() - 1 );
00298         std::vector< double > dr( adjacentCells.size() - 1 );
00299         double bary_x;
00300         double bary_y;
00301 
00302         // get barycenter of cell where reconstruction occurs
00303         std::vector< double > barycent( 3 );
00304         rval = mb->tag_get_data( barycenterTag, &icell, 1, &barycent[0] );CHECK_ERR( rval );
00305         CartVect barycenter( barycent[0], barycent[1], barycent[2] );
00306         double cellbaryx = 0;
00307         double cellbaryy = 0;
00308         moab::IntxUtils::gnomonic_projection( barycenter, R, plane, cellbaryx, cellbaryy );
00309         // int rc = gnomonic_projection_test(barycenter, R, plane, cellbaryx, cellbaryy);
00310 
00311         // get density value
00312         double rhocell = 0;
00313         rval           = mb->tag_get_data( rhoTag, &icell, 1, &rhocell );CHECK_ERR( rval );
00314 
00315         // get barycenters of surrounding cells
00316         std::vector< double > cell_barys( 3 * adjacentCells.size() );
00317         rval = mb->tag_get_data( barycenterTag, adjacentCells, &cell_barys[0] );CHECK_ERR( rval );
00318 
00319         // get density of surrounding cells
00320         std::vector< double > rho_vals( adjacentCells.size() );
00321         rval = mb->tag_get_data( rhoTag, adjacentCells, &rho_vals[0] );CHECK_ERR( rval );
00322 
00323         std::size_t jind = 0;
00324         for( std::size_t i = 0; i < adjacentCells.size(); i++ )
00325         {
00326 
00327             if( adjacentCells[i] != icell )
00328             {
00329                 CartVect bary_xyz( cell_barys[i * 3], cell_barys[i * 3 + 1], cell_barys[i * 3 + 2] );
00330                 moab::IntxUtils::gnomonic_projection( bary_xyz, R, plane, bary_x, bary_y );
00331                 // rc = gnomonic_projection_test(bary_xyz, R, plane, bary_x, bary_y);
00332 
00333                 dx[jind] = bary_x - cellbaryx;
00334                 dy[jind] = bary_y - cellbaryy;
00335                 dr[jind] = rho_vals[i] - rhocell;
00336 
00337                 jind++;
00338             }
00339         }
00340 
00341         std::vector< double > linearCoef( 3 );
00342         if( adjacentCells.size() == 5 )
00343         {
00344 
00345             // compute normal equations matrix
00346             double N11 = dx[1] * dx[1] + dx[2] * dx[2] + dx[3] * dx[3] + dx[0] * dx[0];
00347             double N22 = dy[1] * dy[1] + dy[2] * dy[2] + dy[3] * dy[3] + dy[0] * dy[0];
00348             double N12 = dx[1] * dy[1] + dx[2] * dy[2] + dx[3] * dy[3] + dx[0] * dy[0];
00349 
00350             // rhs
00351             double Rx = dx[1] * dr[1] + dx[2] * dr[2] + dx[3] * dr[3] + dx[0] * dr[0];
00352             double Ry = dy[1] * dr[1] + dy[2] * dr[2] + dy[3] * dr[3] + dy[0] * dr[0];
00353 
00354             // determinant
00355             double Det = N11 * N22 - N12 * N12;
00356 
00357             // solution
00358             linearCoef[0] = ( Rx * N22 - Ry * N12 ) / Det;
00359             linearCoef[1] = ( Ry * N11 - Rx * N12 ) / Det;
00360             linearCoef[2] = rhocell - linearCoef[0] * cellbaryx - linearCoef[1] * cellbaryy;
00361         }
00362         else
00363         {
00364             // default to first order
00365             linearCoef[0] = 0.0;
00366             linearCoef[1] = 0.0;
00367             linearCoef[2] = rhocell;
00368             std::cout << "Need 4 adjacent cells for linear reconstruction! \n";
00369         }
00370 
00371         rval = mb->tag_set_data( linearCoefTag, &icell, 1, &linearCoef[0] );CHECK_ERR( rval );
00372     }
00373     return;
00374 }
00375 
00376 void test_linear_reconstruction( moab::Interface* mb,
00377                                  moab::EntityHandle set,
00378                                  moab::Tag& rhoTag,
00379                                  moab::Tag& planeTag,
00380                                  moab::Tag& barycenterTag,
00381                                  moab::Tag& linearCoefTag )
00382 {
00383     // get all entities of dimension 2
00384     Range cells;
00385     ErrorCode rval = mb->get_entities_by_dimension( set, 2, cells );
00386     if( MB_SUCCESS != rval ) return;
00387 
00388     // set sphere radius to 1
00389     double R = 1;
00390 
00391     // For get coefficients for reconstruction (in cubed-sphere coordinates)
00392     for( Range::iterator it = cells.begin(); it != cells.end(); ++it )
00393     {
00394         moab::EntityHandle icell = *it;
00395 
00396         // get the nodes, then the coordinates
00397         const moab::EntityHandle* verts;
00398         int num_nodes;
00399         rval = mb->get_connectivity( icell, verts, num_nodes );
00400         if( MB_SUCCESS != rval ) return;
00401 
00402         // get coordinates
00403         std::vector< double > coords( 3 * num_nodes );
00404         rval = mb->get_coords( verts, num_nodes, &coords[0] );
00405         if( MB_SUCCESS != rval ) return;
00406 
00407         // get plane for cell
00408         int plane = 0;
00409         rval      = mb->tag_get_data( planeTag, &icell, 1, &plane );CHECK_ERR( rval );
00410 
00411         // get vertex coordinates projections
00412         std::vector< double > x( num_nodes );
00413         std::vector< double > y( num_nodes );
00414         for( int inode = 0; inode < num_nodes; inode++ )
00415         {
00416             double rad = sqrt( coords[inode * 3] * coords[inode * 3] + coords[inode * 3 + 1] * coords[inode * 3 + 1] +
00417                                coords[inode * 3 + 2] * coords[inode * 3 + 2] );
00418             CartVect xyzcoord( coords[inode * 3] / rad, coords[inode * 3 + 1] / rad, coords[inode * 3 + 2] / rad );
00419             moab::IntxUtils::gnomonic_projection( xyzcoord, R, plane, x[inode], y[inode] );
00420             // int dum = gnomonic_projection_test(xyzcoord, R, plane, x[inode],y[inode]);
00421         }
00422 
00423         double area  = 0;
00424         double int_x = 0;
00425         double int_y = 0;
00426         for( int inode = 0; inode < num_nodes; inode++ )
00427         {
00428             int inode2 = inode + 1;
00429             if( inode2 >= num_nodes ) inode2 = 0;
00430             double xmid = 0.5 * ( x[inode] + x[inode2] );
00431             double ymid = 0.5 * ( y[inode] + y[inode2] );
00432             double r1   = sqrt( 1 + x[inode] * x[inode] + y[inode] * y[inode] );
00433             double rm   = sqrt( 1 + xmid * xmid + ymid * ymid );
00434             double r2   = sqrt( 1 + x[inode2] * x[inode2] + y[inode2] * y[inode2] );
00435             double hx   = x[inode2] - x[inode];
00436 
00437             area += hx *
00438                     ( y[inode] / ( r1 * ( 1 + x[inode] * x[inode] ) ) + 4.0 * ymid / ( rm * ( 1 + xmid * xmid ) ) +
00439                       y[inode2] / ( r2 * ( 1 + x[inode2] * x[inode2] ) ) ) /
00440                     6.0;
00441 
00442             int_x += hx *
00443                      ( x[inode] * y[inode] / ( r1 * ( 1 + x[inode] * x[inode] ) ) +
00444                        4.0 * xmid * ymid / ( rm * ( 1 + xmid * xmid ) ) +
00445                        x[inode2] * y[inode2] / ( r2 * ( 1 + x[inode2] * x[inode2] ) ) ) /
00446                      6.0;
00447 
00448             int_y += -hx * ( 1.0 / r1 + 4.0 / rm + 1.0 / r2 ) / 6.0;
00449         }
00450 
00451         // get linear coeficients
00452         std::vector< double > rho_coefs( 3 );
00453         rval = mb->tag_get_data( linearCoefTag, &icell, 1, &rho_coefs[0] );CHECK_ERR( rval );
00454 
00455         // get barycenters
00456         std::vector< double > bary( 3 );
00457         rval = mb->tag_get_data( barycenterTag, &icell, 1, &bary[0] );CHECK_ERR( rval );
00458         double bary_x;
00459         double bary_y;
00460         CartVect bary_xyz( bary[0], bary[1], bary[2] );
00461         moab::IntxUtils::gnomonic_projection( bary_xyz, R, plane, bary_x, bary_y );
00462         // int rc = gnomonic_projection_test(bary_xyz, R, plane, bary_x, bary_y);
00463 
00464         // get cell average density
00465         double cell_rho;
00466         rval = mb->tag_get_data( rhoTag, &icell, 1, &cell_rho );CHECK_ERR( rval );
00467 
00468         // ave rho = \int rho^h(x,y) dV / area = (\int (Ax + By + C) dV) / area
00469         double rho_test1 = ( rho_coefs[0] * int_x + rho_coefs[1] * int_y + rho_coefs[2] * area ) / area;
00470 
00471         // ave rho = A*bary_x + B*bary_y + C
00472         double rho_test2 = rho_coefs[0] * bary_x + rho_coefs[1] * bary_y + rho_coefs[2];
00473 
00474         std::cout << cell_rho << "  " << rho_test1 << "  " << rho_test2 << "  " << cell_rho - rho_test1 << "\n";
00475     }
00476     return;
00477 }
00478 
00479 int reverse_gnomonic_projection_test( const double& c1, const double& c2, double R, int plane, CartVect& pos )
00480 {
00481 
00482     double x = c1;
00483     double y = c2;
00484     double r = sqrt( 1.0 + x * x + y * y );
00485 
00486     switch( plane )
00487     {
00488 
00489         case 1: {
00490             pos[0] = R / r * R;
00491             pos[1] = R / r * x;
00492             pos[2] = R / r * y;
00493             break;
00494         }
00495         case 2: {
00496             pos[0] = -R / r * x;
00497             pos[1] = R / r * R;
00498             pos[2] = R / r * y;
00499             break;
00500         }
00501         case 3: {
00502             pos[0] = -R / r * R;
00503             pos[1] = -R / r * x;
00504             pos[2] = R / r * y;
00505             break;
00506         }
00507         case 4: {
00508             pos[0] = R / r * x;
00509             pos[1] = -R / r * R;
00510             pos[2] = R / r * y;
00511             break;
00512         }
00513         case 5: {
00514             pos[0] = R / r * y;
00515             pos[1] = R / r * x;
00516             pos[2] = -R / r * R;
00517             break;
00518         }
00519         case 6: {
00520             pos[0] = -R / r * y;
00521             pos[1] = R / r * x;
00522             pos[2] = R / r * R;
00523             break;
00524         }
00525     }
00526 
00527     return 0;  // no error
00528 }
00529 
00530 void decide_gnomonic_plane_test( const CartVect& pos, int& plane )
00531 {
00532 
00533     // This is from early version of Homme vorticity calculation in parvis
00534     // Poles are reversed from Homme and Iulian version.
00535 
00536     //  Now has been changed for consistency
00537 
00538     double X = pos[0];
00539     double Y = pos[1];
00540     double Z = pos[2];
00541     double R = sqrt( X * X + Y * Y + Z * Z );
00542     X        = X / R;
00543     Y        = Y / R;
00544     Z        = Z / R;
00545 
00546     if( ( Y < X ) & ( Y > -X ) )
00547     {
00548         if( Z > X )
00549         {
00550             plane = 6;
00551         }
00552         else if( Z < -X )
00553         {
00554             plane = 5;
00555         }
00556         else
00557         {
00558             plane = 1;
00559         }
00560     }
00561     else if( ( Y > X ) & ( Y < -X ) )
00562     {
00563         if( Z > -X )
00564         {
00565             plane = 6;
00566         }
00567         else if( Z < X )
00568         {
00569             plane = 5;
00570         }
00571         else
00572         {
00573             plane = 3;
00574         }
00575     }
00576     else if( ( Y > X ) & ( Y > -X ) )
00577     {
00578         if( Z > Y )
00579         {
00580             plane = 6;
00581         }
00582         else if( Z < -Y )
00583         {
00584             plane = 5;
00585         }
00586         else
00587         {
00588             plane = 2;
00589         }
00590     }
00591     else if( ( Y < X ) & ( Y < -X ) )
00592     {
00593         if( Z > -Y )
00594         {
00595             plane = 6;
00596         }
00597         else if( Z < Y )
00598         {
00599             plane = 5;
00600         }
00601         else
00602         {
00603             plane = 4;
00604         }
00605     }
00606     else
00607     {
00608         if( fabs( X ) < Z )
00609         {
00610             plane = 6;
00611         }
00612         else if( Z < -fabs( X ) )
00613         {
00614             plane = 5;
00615         }
00616         else if( ( X > 0 ) & ( Y > 0 ) )
00617         {
00618             plane = 1;
00619         }
00620         else if( ( X < 0 ) & ( Y > 0 ) )
00621         {
00622             plane = 2;
00623         }
00624         else if( ( X < 0 ) & ( Y < 0 ) )
00625         {
00626             plane = 3;
00627         }
00628         else
00629         {
00630             plane = 4;
00631         }
00632     }
00633 
00634     return;
00635 }
00636 
00637 moab::ErrorCode add_field_value( moab::Interface* mb,
00638                                  moab::EntityHandle euler_set,
00639                                  int /*rank*/,
00640                                  moab::Tag& tagTracer,
00641                                  moab::Tag& tagElem,
00642                                  moab::Tag& tagArea,
00643                                  int field_type )
00644 {
00645 
00646     /*
00647      * get all plys first, then vertices, then move them on the surface of the sphere
00648      *  radius is 1., most of the time
00649      *
00650      */
00651     moab::Range polygons;
00652     moab::ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, polygons );
00653     if( MB_SUCCESS != rval ) return rval;
00654 
00655     moab::Range connecVerts;
00656     rval = mb->get_connectivity( polygons, connecVerts );
00657     if( MB_SUCCESS != rval ) return rval;
00658 
00659     void* data;  // pointer to the LOC in memory, for each vertex
00660     int count;
00661 
00662     rval = mb->tag_iterate( tagTracer, connecVerts.begin(), connecVerts.end(), count, data );CHECK_ERR( rval );
00663     // here we are checking contiguity
00664     assert( count == (int)connecVerts.size() );
00665     double* ptr_DP = (double*)data;
00666     // lambda is for longitude, theta for latitude
00667     // param will be: (la1, te1), (la2, te2), b, c; hmax=1, r=1/2
00668     // nondivergent flow, page 5, case 1, (la1, te1) = (M_PI, M_PI/3)
00669     //                                    (la2, te2) = (M_PI, -M_PI/3)
00670     //                 la1,    te1    la2    te2     b     c  hmax  r
00671     if( field_type == 1 )  // quasi smooth
00672     {
00673         double params[] = { 5 * M_PI / 6.0, 0.0, 7 * M_PI / 6, 0.0, 0.1, 0.9, 1., 0.5 };
00674         for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
00675         {
00676             moab::EntityHandle oldV = *vit;
00677             moab::CartVect posi;
00678             rval = mb->get_coords( &oldV, 1, &( posi[0] ) );CHECK_ERR( rval );
00679 
00680             moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( posi );
00681 
00682             ptr_DP[0] = IntxUtilsCSLAM::quasi_smooth_field( sphCoord.lon, sphCoord.lat, params );
00683 
00684             ptr_DP++;  // increment to the next node
00685         }
00686     }
00687     else if( 2 == field_type )  // smooth
00688     {
00689         moab::CartVect p1, p2;
00690         moab::IntxUtils::SphereCoords spr;
00691         spr.R   = 1;
00692         spr.lat = M_PI / 3;
00693         spr.lon = M_PI;
00694         p1      = moab::IntxUtils::spherical_to_cart( spr );
00695         spr.lat = -M_PI / 3;
00696         p2      = moab::IntxUtils::spherical_to_cart( spr );
00697         //                  x1,    y1,     z1,    x2,   y2,    z2,   h_max, b0
00698         double params[] = { p1[0], p1[1], p1[2], p2[0], p2[1], p2[2], 1, 5. };
00699         for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
00700         {
00701             moab::EntityHandle oldV = *vit;
00702             moab::CartVect posi;
00703             rval = mb->get_coords( &oldV, 1, &( posi[0] ) );CHECK_ERR( rval );
00704 
00705             moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( posi );
00706 
00707             ptr_DP[0] = IntxUtilsCSLAM::smooth_field( sphCoord.lon, sphCoord.lat, params );
00708 
00709             ptr_DP++;  // increment to the next node
00710         }
00711     }
00712     else if( 3 == field_type )  // slotted
00713     {
00714         //                   la1, te1,   la2, te2,       b,   c,   r
00715         double params[] = { M_PI, M_PI / 3, M_PI, -M_PI / 3, 0.1, 0.9, 0.5 };  // no h_max
00716         for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
00717         {
00718             moab::EntityHandle oldV = *vit;
00719             moab::CartVect posi;
00720             rval = mb->get_coords( &oldV, 1, &( posi[0] ) );CHECK_ERR( rval );
00721 
00722             moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( posi );
00723 
00724             ptr_DP[0] = IntxUtilsCSLAM::slotted_cylinder_field( sphCoord.lon, sphCoord.lat, params );
00725             ;
00726 
00727             ptr_DP++;  // increment to the next node
00728         }
00729     }
00730     else if( 4 == field_type )  // constant = 1
00731     {
00732         for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
00733         {
00734             /* moab::EntityHandle oldV = *vit;
00735              moab::CartVect posi;
00736              rval = mb->get_coords(&oldV, 1, &(posi[0]));
00737 
00738              moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical(posi);*/
00739 
00740             ptr_DP[0] = 1.0;
00741 
00742             ptr_DP++;  // increment to the next node
00743         }
00744     }
00745     // add average value for quad/polygon (average corners)
00746     // do some averages
00747 
00748     Range::iterator iter = polygons.begin();
00749     while( iter != polygons.end() )
00750     {
00751         rval = mb->tag_iterate( tagElem, iter, polygons.end(), count, data );CHECK_ERR( rval );
00752         double* ptr = (double*)data;
00753 
00754         rval = mb->tag_iterate( tagArea, iter, polygons.end(), count, data );CHECK_ERR( rval );
00755         double* ptrArea = (double*)data;
00756         for( int i = 0; i < count; i++, ++iter, ptr++, ptrArea++ )
00757         {
00758             const moab::EntityHandle* conn = NULL;
00759             int num_nodes                  = 0;
00760             rval                           = mb->get_connectivity( *iter, conn, num_nodes );CHECK_ERR( rval );
00761             if( num_nodes == 0 ) return MB_FAILURE;
00762             std::vector< double > nodeVals( num_nodes );
00763             double average = 0.;
00764             rval           = mb->tag_get_data( tagTracer, conn, num_nodes, &nodeVals[0] );CHECK_ERR( rval );
00765             for( int j = 0; j < num_nodes; j++ )
00766                 average += nodeVals[j];
00767             average /= num_nodes;
00768             *ptr = average;
00769         }
00770     }
00771 
00772     return MB_SUCCESS;
00773 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines