 |
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
|
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
00001 #include
00002 #include
00003 #include
00004 #include
00005 #include
00006 #include
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 }
1.7.6.1.
Dark theme by Tilen Majerle. All rights reserved.