Intx2Mesh.hpp

Go to the documentation of this file.

/*
 * Intx2Mesh.hpp
 *
 * Created on: Oct 2, 2012
 *
 */
 
#ifndef INTX2MESH_HPP_
#define INTX2MESH_HPP_
 
#include <iostream>
#include <sstream>
#include <fstream>
#include <map>
#include <ctime>
#include <cstdlib>
#include <cstdio>
#include <cstring>
#include "moab/Core.hpp"
#include "moab/Interface.hpp"
#include "moab/Range.hpp"
#include "moab/CartVect.hpp"
#include "moab/IntxMesh/IntxUtils.hpp"
 
// #define ENABLE_DEBUG
 
// maximum number of edges on each convex polygon of interest
#define MAXEDGES 10
#define MAXEDGES2 20 // used for coordinates in plane
#define CORRTAGNAME "__correspondent"
 
namespace moab
{
 
#define ERRORR( rval, str ) \
 if( MB_SUCCESS != ( rval ) ) \
 { \
 std::cout << ( str ) << "\n"; \
 return rval; \
 }
 
#define ERRORV( rval, str ) \
 if( MB_SUCCESS != ( rval ) ) \
 { \
 std::cout << ( str ) << "\n"; \
 return; \
 }
 
#ifdef MOAB_HAVE_MPI
// forward declarations
class ParallelComm;
class TupleList;
#endif
 
class Intx2Mesh
{
 public:
 Intx2Mesh( Interface* mbimpl );
 
 virtual ~Intx2Mesh();
 
 /*
 * computes intersection between 2 sets;
 * set 2 (mbs2) should be completely covered by set mbs1, as all elements
 * from set 2 will be decomposed in polygons ; an area verification is done, and
 * will signal elements from set 2 that are not "recovered"
 *
 * the resulting polygons are inserted in output set; total area from output set should be
 * equal to total area of elements in set 2 (check that!)
 */
 ErrorCode intersect_meshes( EntityHandle mbs1, EntityHandle mbs2, EntityHandle& outputSet );
 
 /*
 * slower intx, use kd tree only, no adjacency, no adv front
 */
 ErrorCode intersect_meshes_kdtree( EntityHandle mbset1, EntityHandle mbset2, EntityHandle& outputSet );
 
 // mark could be (3 or 4, depending on type: ) no, it could go to 10
 // no, it will be MAXEDGES = 10
 // this is pure abstract, this needs to be implemented by
 // all derivations
 // the max number of intersection points could be 2*MAXEDGES
 // so P must be dimensioned to 4*MAXEDGES (2*2*MAXEDGES)
 // so, if you intersect 2 convex polygons with MAXEDGES , you will get a convex polygon
 // with 2*MAXEDGES, at most
 // will also return the number of nodes of tgt and src elements
 virtual ErrorCode computeIntersectionBetweenTgtAndSrc( EntityHandle tgt,
 EntityHandle src,
 double* P,
 int& nP,
 double& area,
 int markb[MAXEDGES],
 int markr[MAXEDGES],
 int& nsidesSrc,
 int& nsidesTgt,
 bool check_boxes_first = false ) = 0;
 
 // this is also abstract
 virtual ErrorCode findNodes( EntityHandle tgt, int nsTgt, EntityHandle src, int nsSrc, double* iP, int nP ) = 0;
 
 // this is also computing the area of the tgt cell in plane (gnomonic plane for sphere)
 // setting the local variables:
 // this will be done once per tgt cell, and once per localQueue for src cells
 // const EntityHandle * tgtConn;
 // CartVect redCoords[MAXEDGES];
 // double redCoords2D[MAXEDGES2]; // these are in plane
 
 virtual double setup_tgt_cell( EntityHandle tgt, int& nsTgt ) = 0;
 
 virtual ErrorCode FindMaxEdgesInSet( EntityHandle eset, int& max_edges );
 virtual ErrorCode FindMaxEdges( EntityHandle set1,
 EntityHandle set2 ); // this needs to be called before any
 // covering communication in parallel
 
 virtual ErrorCode createTags();
 
 virtual ErrorCode filterByMask(Range & cells);
 
 ErrorCode DetermineOrderedNeighbors( EntityHandle inputSet, int max_edges, Tag& neighTag );
 
 void set_error_tolerance( double eps )
 {
 epsilon_1 = eps;
 epsilon_area = eps * sqrt( eps );
 }
 
#ifdef MOAB_HAVE_MPI
 void set_parallel_comm( moab::ParallelComm* pcomm )
 {
 parcomm = pcomm;
 }
#endif
 // void SetEntityType (EntityType tp) { type=tp;}
 
 // clean some memory allocated
 void clean();
 
 void set_box_error( double berror )
 {
 box_error = berror;
 }
 
 ErrorCode create_departure_mesh_2nd_alg( EntityHandle& euler_set, EntityHandle& covering_lagr_set );
 
 // in this method, used in parallel, each departure elements are already created, and at their
 // positions the covering_set is output, will contain the departure cells that cover the euler
 // set; some of these departure cells might come from different processors so the covering_set
 // contains some elements from lagr_set and some elements that come from other procs we need to
 // keep track of what processors "sent" the elements so we know were to send back the info about
 // the tracers masses
 
 ErrorCode create_departure_mesh_3rd_alg( EntityHandle& lagr_set, EntityHandle& covering_set );
 
 /* in this method, used in parallel, each cell from first mesh need to be sent to the
 * tasks whose second mesh bounding boxes intersects bounding box of the cell
 * this method assumes that the bounding boxes for the second mesh were computed already in a
 * previous method (called build_processor_euler_boxes)
 *
 * the covering_set is output, will cover the second mesh (set) from the task;
 * Some of the cells in the first mesh will be sent to multiple processors; we keep track of
 * them using the global id of the vertices and global ids of the cells (we do not want to
 * create multiple vertices with the same ids). The global id of the cells are needed just for
 * debugging, the cells cannot come from 2 different tasks, but the vertices can
 *
 * Right now, we use crystal router, but an improvement might be to use direct communication
 * (send_entities) on parallel comm
 *
 * param initial_distributed_set (IN) : the initial distribution of the first mesh (set)
 *
 * param (OUT) : the covering set in first mesh , which completely covers the second mesh set
 */
#ifdef MOAB_HAVE_MPI
 
 virtual ErrorCode build_processor_euler_boxes( EntityHandle euler_set, Range& local_verts, bool gnomonic = true );
#endif
 void correct_polygon( EntityHandle* foundIds, int& nP );
#ifdef MOAB_HAVE_MPI
 // share vertices between the intersection target domains
 ErrorCode resolve_intersection_sharing();
#endif
#ifdef ENABLE_DEBUG
 void enable_debug()
 {
 dbg_1 = 1;
 }
 void disable_debug()
 {
 dbg_1 = 0;
 }
#endif
 
#ifdef MOAB_HAVE_TEMPESTREMAP
 friend class TempestRemapper;
#endif
 
 protected: // so it can be accessed in derived classes, InPlane and OnSphere
 Interface* mb;
 
 EntityHandle mbs1;
 EntityHandle mbs2;
 Range rs1; // range set 1 (departure set, lagrange set, src set, manufactured set, target mesh)
 Range rs2; // range set 2 (arrival set, euler set, tgt set, initial set, source mesh)
 
 EntityHandle outSet; // will contain intersection
 Tag gid; // global id tag will be used to set the parents of the intersection cell
 
 // tags used in computation, advancing front
 Tag TgtFlagTag; // to mark tgt quads already considered
 
 Range TgtEdges; //
 
 // tgt parent and src parent tags
 // these will be on the out mesh
 Tag tgtParentTag;
 Tag srcParentTag;
 Tag countTag;
 
 Tag srcNeighTag; // will store neighbors for navigating easily in advancing front, for first
 // mesh (src, target, lagrange)
 Tag tgtNeighTag; // will store neighbors for navigating easily in advancing front, for second
 // mesh (tgt, source, euler)
 
 Tag neighTgtEdgeTag; // will store edge borders for each tgt cell
 
 Tag orgSendProcTag; /// for coverage mesh, will store the original sender
 Tag imaskTag; // if it exists, use it for filtering the source or target cells
 
 // EntityType type; // this will be tri, quad or MBPOLYGON...
 
 const EntityHandle* tgtConn;
 const EntityHandle* srcConn;
 CartVect tgtCoords[MAXEDGES];
 CartVect srcCoords[MAXEDGES];
 double tgtCoords2D[MAXEDGES2]; // these are in plane
 double srcCoords2D[MAXEDGES2]; // these are in plane
 
#ifdef ENABLE_DEBUG
 static int dbg_1;
 std::ofstream mout_1[6]; // some debug files
#endif
 // for each tgt edge, we keep a vector of extra nodes, coming from intersections
 // use the index in TgtEdges range
 // so the extra nodes on each tgt edge are kept track of
 // only entity handles are in the vector, not the actual coordinates;
 // actual coordinates are retrieved every time, which could be expensive
 // maybe we should store the coordinates too, along with entity handles (more memory used,
 // faster to retrieve)
 std::vector< std::vector< EntityHandle >* > extraNodesVec;
 
 double epsilon_1;
 double epsilon_area;
 
 std::vector< double > allBoxes;
 double box_error;
 /* \brief Local root of the kdtree */
 EntityHandle localRoot;
 Range localEnts; // this range is for local elements of interest, euler cells, or "first mesh"
 unsigned int my_rank;
 
#ifdef MOAB_HAVE_MPI
 ParallelComm* parcomm;
 TupleList* remote_cells; // not used anymore for communication, just a container
 TupleList* remote_cells_with_tracers; // these will be used now to update tracers on remote procs
 std::map< int, EntityHandle > globalID_to_eh; // needed for parallel, mostly
#endif
 int max_edges_1; // maximum number of edges in the lagrange set (first set, src)
 int max_edges_2; // maximum number of edges in the euler set (second set, tgt)
 int counting;
};
 
} /* namespace moab */
#endif /* INTX2MESH_HPP_ */