|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
#include <iostream>#include <cmath>#include "moab/Core.hpp"#include "moab/Interface.hpp"#include "moab/IntxMesh/Intx2MeshOnSphere.hpp"#include "moab/ParallelComm.hpp"#include "moab/ProgOptions.hpp"#include "MBParallelConventions.h"#include "moab/ReadUtilIface.hpp"#include "MBTagConventions.hpp"#include "moab/IntxMesh/IntxUtils.hpp"#include "IntxUtilsCSLAM.hpp"#include "TestUtil.hpp"
Include dependency graph for linear_remap.cpp:Go to the source code of this file.
| moab::ErrorCode add_field_value | ( | moab::Interface * | mb, |
| moab::EntityHandle | euler_set, | ||
| int | rank, | ||
| moab::Tag & | tagTracer, | ||
| moab::Tag & | tagElem, | ||
| moab::Tag & | tagArea, | ||
| int | field_type | ||
| ) |
Definition at line 637 of file linear_remap.cpp.
References moab::Range::begin(), moab::IntxUtils::cart_to_spherical(), CHECK_ERR, moab::Range::end(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), moab::IntxUtils::SphereCoords::lat, moab::IntxUtils::SphereCoords::lon, MB_SUCCESS, IntxUtilsCSLAM::quasi_smooth_field(), moab::IntxUtils::SphereCoords::R, moab::Range::size(), IntxUtilsCSLAM::slotted_cylinder_field(), IntxUtilsCSLAM::smooth_field(), moab::IntxUtils::spherical_to_cart(), moab::Interface::tag_get_data(), and moab::Interface::tag_iterate().
{
/*
* get all plys first, then vertices, then move them on the surface of the sphere
* radius is 1., most of the time
*
*/
moab::Range polygons;
moab::ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, polygons );
if( MB_SUCCESS != rval ) return rval;
moab::Range connecVerts;
rval = mb->get_connectivity( polygons, connecVerts );
if( MB_SUCCESS != rval ) return rval;
void* data; // pointer to the LOC in memory, for each vertex
int count;
rval = mb->tag_iterate( tagTracer, connecVerts.begin(), connecVerts.end(), count, data );CHECK_ERR( rval );
// here we are checking contiguity
assert( count == (int)connecVerts.size() );
double* ptr_DP = (double*)data;
// lambda is for longitude, theta for latitude
// param will be: (la1, te1), (la2, te2), b, c; hmax=1, r=1/2
// nondivergent flow, page 5, case 1, (la1, te1) = (M_PI, M_PI/3)
// (la2, te2) = (M_PI, -M_PI/3)
// la1, te1 la2 te2 b c hmax r
if( field_type == 1 ) // quasi smooth
{
double params[] = { 5 * M_PI / 6.0, 0.0, 7 * M_PI / 6, 0.0, 0.1, 0.9, 1., 0.5 };
for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
{
moab::EntityHandle oldV = *vit;
moab::CartVect posi;
rval = mb->get_coords( &oldV, 1, &( posi[0] ) );CHECK_ERR( rval );
moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( posi );
ptr_DP[0] = IntxUtilsCSLAM::quasi_smooth_field( sphCoord.lon, sphCoord.lat, params );
ptr_DP++; // increment to the next node
}
}
else if( 2 == field_type ) // smooth
{
moab::CartVect p1, p2;
moab::IntxUtils::SphereCoords spr;
spr.R = 1;
spr.lat = M_PI / 3;
spr.lon = M_PI;
p1 = moab::IntxUtils::spherical_to_cart( spr );
spr.lat = -M_PI / 3;
p2 = moab::IntxUtils::spherical_to_cart( spr );
// x1, y1, z1, x2, y2, z2, h_max, b0
double params[] = { p1[0], p1[1], p1[2], p2[0], p2[1], p2[2], 1, 5. };
for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
{
moab::EntityHandle oldV = *vit;
moab::CartVect posi;
rval = mb->get_coords( &oldV, 1, &( posi[0] ) );CHECK_ERR( rval );
moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( posi );
ptr_DP[0] = IntxUtilsCSLAM::smooth_field( sphCoord.lon, sphCoord.lat, params );
ptr_DP++; // increment to the next node
}
}
else if( 3 == field_type ) // slotted
{
// la1, te1, la2, te2, b, c, r
double params[] = { M_PI, M_PI / 3, M_PI, -M_PI / 3, 0.1, 0.9, 0.5 }; // no h_max
for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
{
moab::EntityHandle oldV = *vit;
moab::CartVect posi;
rval = mb->get_coords( &oldV, 1, &( posi[0] ) );CHECK_ERR( rval );
moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( posi );
ptr_DP[0] = IntxUtilsCSLAM::slotted_cylinder_field( sphCoord.lon, sphCoord.lat, params );
;
ptr_DP++; // increment to the next node
}
}
else if( 4 == field_type ) // constant = 1
{
for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
{
/* moab::EntityHandle oldV = *vit;
moab::CartVect posi;
rval = mb->get_coords(&oldV, 1, &(posi[0]));
moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical(posi);*/
ptr_DP[0] = 1.0;
ptr_DP++; // increment to the next node
}
}
// add average value for quad/polygon (average corners)
// do some averages
Range::iterator iter = polygons.begin();
while( iter != polygons.end() )
{
rval = mb->tag_iterate( tagElem, iter, polygons.end(), count, data );CHECK_ERR( rval );
double* ptr = (double*)data;
rval = mb->tag_iterate( tagArea, iter, polygons.end(), count, data );CHECK_ERR( rval );
double* ptrArea = (double*)data;
for( int i = 0; i < count; i++, ++iter, ptr++, ptrArea++ )
{
const moab::EntityHandle* conn = NULL;
int num_nodes = 0;
rval = mb->get_connectivity( *iter, conn, num_nodes );CHECK_ERR( rval );
if( num_nodes == 0 ) return MB_FAILURE;
std::vector< double > nodeVals( num_nodes );
double average = 0.;
rval = mb->tag_get_data( tagTracer, conn, num_nodes, &nodeVals[0] );CHECK_ERR( rval );
for( int j = 0; j < num_nodes; j++ )
average += nodeVals[j];
average /= num_nodes;
*ptr = average;
}
}
return MB_SUCCESS;
}
| void decide_gnomonic_plane_test | ( | const CartVect & | pos, |
| int & | plane | ||
| ) |
Definition at line 530 of file linear_remap.cpp.
{
// This is from early version of Homme vorticity calculation in parvis
// Poles are reversed from Homme and Iulian version.
// Now has been changed for consistency
double X = pos[0];
double Y = pos[1];
double Z = pos[2];
double R = sqrt( X * X + Y * Y + Z * Z );
X = X / R;
Y = Y / R;
Z = Z / R;
if( ( Y < X ) & ( Y > -X ) )
{
if( Z > X )
{
plane = 6;
}
else if( Z < -X )
{
plane = 5;
}
else
{
plane = 1;
}
}
else if( ( Y > X ) & ( Y < -X ) )
{
if( Z > -X )
{
plane = 6;
}
else if( Z < X )
{
plane = 5;
}
else
{
plane = 3;
}
}
else if( ( Y > X ) & ( Y > -X ) )
{
if( Z > Y )
{
plane = 6;
}
else if( Z < -Y )
{
plane = 5;
}
else
{
plane = 2;
}
}
else if( ( Y < X ) & ( Y < -X ) )
{
if( Z > -Y )
{
plane = 6;
}
else if( Z < Y )
{
plane = 5;
}
else
{
plane = 4;
}
}
else
{
if( fabs( X ) < Z )
{
plane = 6;
}
else if( Z < -fabs( X ) )
{
plane = 5;
}
else if( ( X > 0 ) & ( Y > 0 ) )
{
plane = 1;
}
else if( ( X < 0 ) & ( Y > 0 ) )
{
plane = 2;
}
else if( ( X < 0 ) & ( Y < 0 ) )
{
plane = 3;
}
else
{
plane = 4;
}
}
return;
}
| void get_barycenters | ( | moab::Interface * | mb, |
| moab::EntityHandle | set, | ||
| moab::Tag & | planeTag, | ||
| moab::Tag & | barycenterTag | ||
| ) |
Definition at line 172 of file linear_remap.cpp.
References moab::Range::begin(), CHECK_ERR, moab::Range::end(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), moab::IntxUtils::gnomonic_projection(), MB_SUCCESS, moab::R, moab::IntxUtils::reverse_gnomonic_projection(), moab::Interface::tag_get_data(), and moab::Interface::tag_set_data().
{
// get all entities of dimension 2
moab::Range cells;
ErrorCode rval = mb->get_entities_by_dimension( set, 2, cells );
if( MB_SUCCESS != rval ) return;
// set sphere radius to 1
double R = 1.0;
for( Range::iterator it = cells.begin(); it != cells.end(); ++it )
{
moab::EntityHandle icell = *it;
// get the nodes
const moab::EntityHandle* verts;
int num_nodes;
rval = mb->get_connectivity( icell, verts, num_nodes );
if( MB_SUCCESS != rval ) return;
// get coordinates
std::vector< double > coords( 3 * num_nodes );
rval = mb->get_coords( verts, num_nodes, &coords[0] );
if( MB_SUCCESS != rval ) return;
// get plane for cell
int plane = 0;
rval = mb->tag_get_data( planeTag, &icell, 1, &plane );
if( MB_SUCCESS != rval ) return;
std::vector< double > x( num_nodes );
std::vector< double > y( num_nodes );
double area = 0;
double bary_x = 0;
double bary_y = 0;
for( int inode = 0; inode < num_nodes; inode++ )
{
// radius should be 1.0, but divide by it just in case for now
double rad = sqrt( coords[inode * 3] * coords[inode * 3] + coords[inode * 3 + 1] * coords[inode * 3 + 1] +
coords[inode * 3 + 2] * coords[inode * 3 + 2] );
CartVect xyzcoord( coords[inode * 3] / rad, coords[inode * 3 + 1] / rad, coords[inode * 3 + 2] / rad );
moab::IntxUtils::gnomonic_projection( xyzcoord, R, plane, x[inode], y[inode] );
// int dum = gnomonic_projection_test(xyzcoord, R, plane, x[inode],y[inode]);
}
for( int inode = 0; inode < num_nodes; inode++ )
{
int inode2 = inode + 1;
if( inode2 >= num_nodes ) inode2 = 0;
double xmid = 0.5 * ( x[inode] + x[inode2] );
double ymid = 0.5 * ( y[inode] + y[inode2] );
double r1 = sqrt( 1 + x[inode] * x[inode] + y[inode] * y[inode] );
double rm = sqrt( 1 + xmid * xmid + ymid * ymid );
double r2 = sqrt( 1 + x[inode2] * x[inode2] + y[inode2] * y[inode2] );
double hx = x[inode2] - x[inode];
area += hx *
( y[inode] / ( r1 * ( 1 + x[inode] * x[inode] ) ) + 4.0 * ymid / ( rm * ( 1 + xmid * xmid ) ) +
y[inode2] / ( r2 * ( 1 + x[inode2] * x[inode2] ) ) ) /
6.0;
bary_x += hx *
( x[inode] * y[inode] / ( r1 * ( 1 + x[inode] * x[inode] ) ) +
4.0 * xmid * ymid / ( rm * ( 1 + xmid * xmid ) ) +
x[inode2] * y[inode2] / ( r2 * ( 1 + x[inode2] * x[inode2] ) ) ) /
6.0;
bary_y += -hx * ( 1.0 / r1 + 4.0 / rm + 1.0 / r2 ) / 6.0;
}
bary_x = bary_x / area;
bary_y = bary_y / area;
moab::CartVect barycent;
moab::IntxUtils::reverse_gnomonic_projection( bary_x, bary_y, R, plane, barycent );
// reverse_gnomonic_projection_test(bary_x, bary_y, R, plane, barycent);
std::vector< double > barycenter( 3 );
barycenter[0] = barycent[0];
barycenter[1] = barycent[1];
barycenter[2] = barycent[2];
rval = mb->tag_set_data( barycenterTag, &icell, 1, &barycenter[0] );CHECK_ERR( rval );
}
return;
}
| void get_gnomonic_plane | ( | moab::Interface * | mb, |
| moab::EntityHandle | set, | ||
| moab::Tag & | planeTag | ||
| ) |
| void get_linear_reconstruction | ( | moab::Interface * | mb, |
| moab::EntityHandle | set, | ||
| moab::Tag & | rhoTag, | ||
| moab::Tag & | planeTag, | ||
| moab::Tag & | barycenterTag, | ||
| moab::Tag & | linearCoefTag | ||
| ) |
| int gnomonic_projection_test | ( | const moab::CartVect & | pos, |
| double | R, | ||
| int | plane, | ||
| double & | c1, | ||
| double & | c2 | ||
| ) |
| int main | ( | int | argc, |
| char * | argv[] | ||
| ) |
Definition at line 62 of file linear_remap.cpp.
References add_field_value(), CHECK_ERR, moab::Interface::create_meshset(), ErrorCode, get_barycenters(), get_gnomonic_plane(), get_linear_reconstruction(), moab::Interface::load_file(), mb, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_DOUBLE, MB_TYPE_INTEGER, MESHSET_SET, MPI_COMM_WORLD, moab::ParallelComm::proc_config(), moab::ProcConfig::proc_rank(), rank, moab::Interface::tag_get_handle(), and test_linear_reconstruction().
{
// set up MOAB interface and parallel communication
MPI_Init( &argc, &argv );
moab::Core moab;
moab::Interface& mb = moab;
moab::ParallelComm mb_pcomm( &mb, MPI_COMM_WORLD );
// int rank = mb_pcomm->proc_config().proc_rank();
int rank = mb_pcomm.proc_config().proc_rank();
// create meshset
moab::EntityHandle euler_set;
moab::ErrorCode rval = mb.create_meshset( MESHSET_SET, euler_set );CHECK_ERR( rval );
// std::stringstream opts;
// opts <<
// "PARALLEL=READ_PART;PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;GATHER_SET=0;PARTITION_METHOD=TRIVIAL_PARTITION;VARIABLE=";
// opts << "PARALLEL=READ_PART;PARTITION;PARALLEL_RESOLVE_SHARED_ENTS";
std::string fileN = TestDir + "unittest/mbcslam/fine4.h5m";
rval = mb.load_file( fileN.c_str(), &euler_set );CHECK_ERR( rval );
// Create tag for cell density
moab::Tag rhoTag = 0;
rval = mb.tag_get_handle( "Density", 1, moab::MB_TYPE_DOUBLE, rhoTag, moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );CHECK_ERR( rval );
moab::Tag rhoNodeTag = 0;
rval = mb.tag_get_handle( "DensityNode", 1, moab::MB_TYPE_DOUBLE, rhoNodeTag,
moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );CHECK_ERR( rval );
// Create tag for cell area
moab::Tag areaTag = 0;
rval = mb.tag_get_handle( "Area", 1, moab::MB_TYPE_DOUBLE, areaTag, moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );CHECK_ERR( rval );
// Create tag for cell barycenters in 3D Cartesian space
moab::Tag barycenterTag = 0;
rval = mb.tag_get_handle( "CellBarycenter", 3, moab::MB_TYPE_DOUBLE, barycenterTag,
moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );CHECK_ERR( rval );
// Create tag for cell density reconstruction coefficients
moab::Tag coefRhoTag = 0;
rval = mb.tag_get_handle( "LinearCoefRho", 3, moab::MB_TYPE_DOUBLE, coefRhoTag,
moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );CHECK_ERR( rval );
// Create tag for index of gnomonic plane for each cell
moab::Tag planeTag = 0;
rval = mb.tag_get_handle( "gnomonicPlane", 1, moab::MB_TYPE_INTEGER, planeTag,
moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );CHECK_ERR( rval );
// Set density distributions
rval = add_field_value( &mb, euler_set, rank, rhoNodeTag, rhoTag, areaTag, 2 );CHECK_ERR( rval );
// get cell plane
get_gnomonic_plane( &mb, euler_set, planeTag );
// get cell barycenters
get_barycenters( &mb, euler_set, planeTag, barycenterTag );
// get linear reconstruction
get_linear_reconstruction( &mb, euler_set, rhoTag, planeTag, barycenterTag, coefRhoTag );
// test linear reconstruction
test_linear_reconstruction( &mb, euler_set, rhoTag, planeTag, barycenterTag, coefRhoTag );
MPI_Finalize();
return 0;
}
| int reverse_gnomonic_projection_test | ( | const double & | c1, |
| const double & | c2, | ||
| double | R, | ||
| int | plane, | ||
| moab::CartVect & | pos | ||
| ) |
Definition at line 479 of file linear_remap.cpp.
References moab::R.
{
double x = c1;
double y = c2;
double r = sqrt( 1.0 + x * x + y * y );
switch( plane )
{
case 1: {
pos[0] = R / r * R;
pos[1] = R / r * x;
pos[2] = R / r * y;
break;
}
case 2: {
pos[0] = -R / r * x;
pos[1] = R / r * R;
pos[2] = R / r * y;
break;
}
case 3: {
pos[0] = -R / r * R;
pos[1] = -R / r * x;
pos[2] = R / r * y;
break;
}
case 4: {
pos[0] = R / r * x;
pos[1] = -R / r * R;
pos[2] = R / r * y;
break;
}
case 5: {
pos[0] = R / r * y;
pos[1] = R / r * x;
pos[2] = -R / r * R;
break;
}
case 6: {
pos[0] = -R / r * y;
pos[1] = R / r * x;
pos[2] = R / r * R;
break;
}
}
return 0; // no error
}
| void test_linear_reconstruction | ( | moab::Interface * | mb, |
| moab::EntityHandle | set, | ||
| moab::Tag & | rhoTag, | ||
| moab::Tag & | planeTag, | ||
| moab::Tag & | barycenterTag, | ||
| moab::Tag & | linearCoefTag | ||
| ) |
Definition at line 376 of file linear_remap.cpp.
References moab::Range::begin(), CHECK_ERR, moab::Range::end(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), moab::IntxUtils::gnomonic_projection(), MB_SUCCESS, and moab::Interface::tag_get_data().
Referenced by main().
{
// get all entities of dimension 2
Range cells;
ErrorCode rval = mb->get_entities_by_dimension( set, 2, cells );
if( MB_SUCCESS != rval ) return;
// set sphere radius to 1
double R = 1;
// For get coefficients for reconstruction (in cubed-sphere coordinates)
for( Range::iterator it = cells.begin(); it != cells.end(); ++it )
{
moab::EntityHandle icell = *it;
// get the nodes, then the coordinates
const moab::EntityHandle* verts;
int num_nodes;
rval = mb->get_connectivity( icell, verts, num_nodes );
if( MB_SUCCESS != rval ) return;
// get coordinates
std::vector< double > coords( 3 * num_nodes );
rval = mb->get_coords( verts, num_nodes, &coords[0] );
if( MB_SUCCESS != rval ) return;
// get plane for cell
int plane = 0;
rval = mb->tag_get_data( planeTag, &icell, 1, &plane );CHECK_ERR( rval );
// get vertex coordinates projections
std::vector< double > x( num_nodes );
std::vector< double > y( num_nodes );
for( int inode = 0; inode < num_nodes; inode++ )
{
double rad = sqrt( coords[inode * 3] * coords[inode * 3] + coords[inode * 3 + 1] * coords[inode * 3 + 1] +
coords[inode * 3 + 2] * coords[inode * 3 + 2] );
CartVect xyzcoord( coords[inode * 3] / rad, coords[inode * 3 + 1] / rad, coords[inode * 3 + 2] / rad );
moab::IntxUtils::gnomonic_projection( xyzcoord, R, plane, x[inode], y[inode] );
// int dum = gnomonic_projection_test(xyzcoord, R, plane, x[inode],y[inode]);
}
double area = 0;
double int_x = 0;
double int_y = 0;
for( int inode = 0; inode < num_nodes; inode++ )
{
int inode2 = inode + 1;
if( inode2 >= num_nodes ) inode2 = 0;
double xmid = 0.5 * ( x[inode] + x[inode2] );
double ymid = 0.5 * ( y[inode] + y[inode2] );
double r1 = sqrt( 1 + x[inode] * x[inode] + y[inode] * y[inode] );
double rm = sqrt( 1 + xmid * xmid + ymid * ymid );
double r2 = sqrt( 1 + x[inode2] * x[inode2] + y[inode2] * y[inode2] );
double hx = x[inode2] - x[inode];
area += hx *
( y[inode] / ( r1 * ( 1 + x[inode] * x[inode] ) ) + 4.0 * ymid / ( rm * ( 1 + xmid * xmid ) ) +
y[inode2] / ( r2 * ( 1 + x[inode2] * x[inode2] ) ) ) /
6.0;
int_x += hx *
( x[inode] * y[inode] / ( r1 * ( 1 + x[inode] * x[inode] ) ) +
4.0 * xmid * ymid / ( rm * ( 1 + xmid * xmid ) ) +
x[inode2] * y[inode2] / ( r2 * ( 1 + x[inode2] * x[inode2] ) ) ) /
6.0;
int_y += -hx * ( 1.0 / r1 + 4.0 / rm + 1.0 / r2 ) / 6.0;
}
// get linear coeficients
std::vector< double > rho_coefs( 3 );
rval = mb->tag_get_data( linearCoefTag, &icell, 1, &rho_coefs[0] );CHECK_ERR( rval );
// get barycenters
std::vector< double > bary( 3 );
rval = mb->tag_get_data( barycenterTag, &icell, 1, &bary[0] );CHECK_ERR( rval );
double bary_x;
double bary_y;
CartVect bary_xyz( bary[0], bary[1], bary[2] );
moab::IntxUtils::gnomonic_projection( bary_xyz, R, plane, bary_x, bary_y );
// int rc = gnomonic_projection_test(bary_xyz, R, plane, bary_x, bary_y);
// get cell average density
double cell_rho;
rval = mb->tag_get_data( rhoTag, &icell, 1, &cell_rho );CHECK_ERR( rval );
// ave rho = \int rho^h(x,y) dV / area = (\int (Ax + By + C) dV) / area
double rho_test1 = ( rho_coefs[0] * int_x + rho_coefs[1] * int_y + rho_coefs[2] * area ) / area;
// ave rho = A*bary_x + B*bary_y + C
double rho_test2 = rho_coefs[0] * bary_x + rho_coefs[1] * bary_y + rho_coefs[2];
std::cout << cell_rho << " " << rho_test1 << " " << rho_test2 << " " << cell_rho - rho_test1 << "\n";
}
return;
}
| double radius = 1. |
Definition at line 60 of file linear_remap.cpp.
1.7.6.1