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