MOAB: Mesh Oriented datABase  (version 5.4.0)
Intx2Mesh.hpp
Go to the documentation of this file.
00001 /*
00002  * Intx2Mesh.hpp
00003  *
00004  *  Created on: Oct 2, 2012
00005  *
00006  */
00007 
00008 #ifndef INTX2MESH_HPP_
00009 #define INTX2MESH_HPP_
00010 
00011 #include <iostream>
00012 #include <sstream>
00013 #include <fstream>
00014 #include <map>
00015 #include <ctime>
00016 #include <cstdlib>
00017 #include <cstdio>
00018 #include <cstring>
00019 #include "moab/Core.hpp"
00020 #include "moab/Interface.hpp"
00021 #include "moab/Range.hpp"
00022 #include "moab/CartVect.hpp"
00023 #include "moab/IntxMesh/IntxUtils.hpp"
00024 
00025 // #define ENABLE_DEBUG
00026 
00027 // maximum number of edges on each convex polygon of interest
00028 #define MAXEDGES    10
00029 #define MAXEDGES2   20  // used for coordinates in plane
00030 #define CORRTAGNAME "__correspondent"
00031 
00032 namespace moab
00033 {
00034 
00035 #define ERRORR( rval, str )           \
00036     if( MB_SUCCESS != ( rval ) )      \
00037     {                                 \
00038         std::cout << ( str ) << "\n"; \
00039         return rval;                  \
00040     }
00041 
00042 #define ERRORV( rval, str )           \
00043     if( MB_SUCCESS != ( rval ) )      \
00044     {                                 \
00045         std::cout << ( str ) << "\n"; \
00046         return;                       \
00047     }
00048 
00049 #ifdef MOAB_HAVE_MPI
00050 // forward declarations
00051 class ParallelComm;
00052 class TupleList;
00053 #endif
00054 
00055 class Intx2Mesh
00056 {
00057   public:
00058     Intx2Mesh( Interface* mbimpl );
00059 
00060     virtual ~Intx2Mesh();
00061 
00062     /*
00063      * computes intersection between 2 sets;
00064      * set 2 (mbs2) should be completely covered by set mbs1, as all elements
00065      * from set 2 will be decomposed in polygons ; an area verification is done, and
00066      * will signal elements from set 2 that are not "recovered"
00067      *
00068      * the resulting polygons are inserted in output set; total area from output set should be
00069      * equal to total area of elements in set 2 (check that!)
00070      */
00071     ErrorCode intersect_meshes( EntityHandle mbs1, EntityHandle mbs2, EntityHandle& outputSet );
00072 
00073     /*
00074      *  slower intx, use kd tree only, no adjacency, no adv front
00075      */
00076     ErrorCode intersect_meshes_kdtree( EntityHandle mbset1, EntityHandle mbset2, EntityHandle& outputSet );
00077 
00078     // mark could be (3 or 4, depending on type: ) no, it could go to 10
00079     // no, it will be MAXEDGES = 10
00080     // this is pure abstract, this needs to be implemented by
00081     // all derivations
00082     // the max number of intersection points could be 2*MAXEDGES
00083     // so P must be dimensioned to 4*MAXEDGES (2*2*MAXEDGES)
00084     // so, if you intersect 2 convex polygons with MAXEDGES , you will get a convex polygon
00085     // with 2*MAXEDGES, at most
00086     // will also return the number of nodes of tgt and src elements
00087     virtual ErrorCode computeIntersectionBetweenTgtAndSrc( EntityHandle tgt,
00088                                                            EntityHandle src,
00089                                                            double* P,
00090                                                            int& nP,
00091                                                            double& area,
00092                                                            int markb[MAXEDGES],
00093                                                            int markr[MAXEDGES],
00094                                                            int& nsidesSrc,
00095                                                            int& nsidesTgt,
00096                                                            bool check_boxes_first = false ) = 0;
00097 
00098     // this is also abstract
00099     virtual ErrorCode findNodes( EntityHandle tgt, int nsTgt, EntityHandle src, int nsSrc, double* iP, int nP ) = 0;
00100 
00101     // this is also computing the area of the tgt cell in plane (gnomonic plane for sphere)
00102     // setting the local variables:
00103     // this will be done once per tgt cell, and once per localQueue for src cells
00104     //  const EntityHandle * tgtConn;
00105     // CartVect redCoords[MAXEDGES];
00106     // double redCoords2D[MAXEDGES2]; // these are in plane
00107 
00108     virtual double setup_tgt_cell( EntityHandle tgt, int& nsTgt ) = 0;
00109 
00110     virtual ErrorCode FindMaxEdgesInSet( EntityHandle eset, int& max_edges );
00111     virtual ErrorCode FindMaxEdges( EntityHandle set1,
00112                                     EntityHandle set2 );  // this needs to be called before any
00113                                                           // covering communication in parallel
00114 
00115     virtual ErrorCode createTags();
00116 
00117     ErrorCode DetermineOrderedNeighbors( EntityHandle inputSet, int max_edges, Tag& neighTag );
00118 
00119     void set_error_tolerance( double eps )
00120     {
00121         epsilon_1    = eps;
00122         epsilon_area = eps * sqrt( eps );
00123     }
00124 
00125 #ifdef MOAB_HAVE_MPI
00126     void set_parallel_comm( moab::ParallelComm* pcomm )
00127     {
00128         parcomm = pcomm;
00129     }
00130 #endif
00131     // void SetEntityType (EntityType tp) { type=tp;}
00132 
00133     // clean some memory allocated
00134     void clean();
00135 
00136     void set_box_error( double berror )
00137     {
00138         box_error = berror;
00139     }
00140 
00141     ErrorCode create_departure_mesh_2nd_alg( EntityHandle& euler_set, EntityHandle& covering_lagr_set );
00142 
00143     // in this method, used in parallel, each departure elements are already created, and at their
00144     // positions the covering_set is output, will contain the departure cells that cover the euler
00145     // set; some of these departure cells might come from different processors so the covering_set
00146     // contains some elements from lagr_set and some elements that come from other procs we need to
00147     // keep track of what processors "sent" the elements so we know were to send back the info about
00148     // the tracers masses
00149 
00150     ErrorCode create_departure_mesh_3rd_alg( EntityHandle& lagr_set, EntityHandle& covering_set );
00151 
00152     /* in this method, used in parallel, each cell from first mesh need to be sent to the
00153      * tasks whose second mesh bounding boxes intersects bounding box of the cell
00154      * this method assumes that the bounding boxes for the second mesh were computed already in a
00155      * previous method (called build_processor_euler_boxes)
00156      *
00157      * the covering_set is output,  will cover the second mesh (set) from the task;
00158      * Some of the cells in the first mesh will be sent to multiple processors; we keep track of
00159      * them using the global id of the  vertices and global ids of the cells (we do not want to
00160      * create multiple vertices with the same ids). The global id of the cells are needed just for
00161      * debugging, the cells cannot come from 2 different tasks, but the vertices can
00162      *
00163      * Right now, we use crystal router, but an improvement might be to use direct communication
00164      * (send_entities) on parallel comm
00165      *
00166      * param initial_distributed_set (IN) : the initial distribution of the first mesh (set)
00167      *
00168      * param (OUT) : the covering set in first mesh , which completely covers the second mesh set
00169      */
00170 #ifdef MOAB_HAVE_MPI
00171 
00172     virtual ErrorCode build_processor_euler_boxes( EntityHandle euler_set, Range& local_verts );
00173 #endif
00174     void correct_polygon( EntityHandle* foundIds, int& nP );
00175 #ifdef MOAB_HAVE_MPI
00176     // share vertices between the intersection target domains
00177     ErrorCode resolve_intersection_sharing();
00178 #endif
00179 #ifdef ENABLE_DEBUG
00180     void enable_debug()
00181     {
00182         dbg_1 = 1;
00183     }
00184     void disable_debug()
00185     {
00186         dbg_1 = 0;
00187     }
00188 #endif
00189 
00190 #ifdef MOAB_HAVE_TEMPESTREMAP
00191     friend class TempestRemapper;
00192 #endif
00193 
00194   protected:  // so it can be accessed in derived classes, InPlane and OnSphere
00195     Interface* mb;
00196 
00197     EntityHandle mbs1;
00198     EntityHandle mbs2;
00199     Range rs1;  // range set 1 (departure set, lagrange set, src set, manufactured set, target mesh)
00200     Range rs2;  // range set 2 (arrival set, euler set, tgt set, initial set, source mesh)
00201 
00202     EntityHandle outSet;  // will contain intersection
00203     Tag gid;              // global id tag will be used to set the parents of the intersection cell
00204 
00205     // tags used in computation, advancing front
00206     Tag TgtFlagTag;  // to mark tgt quads already considered
00207 
00208     Range TgtEdges;  //
00209 
00210     // tgt parent and src parent tags
00211     // these will be on the out mesh
00212     Tag tgtParentTag;
00213     Tag srcParentTag;
00214     Tag countTag;
00215 
00216     Tag srcNeighTag;  // will store neighbors for navigating easily in advancing front, for first
00217                       // mesh (src, target, lagrange)
00218     Tag tgtNeighTag;  // will store neighbors for navigating easily in advancing front, for second
00219                       // mesh (tgt, source, euler)
00220 
00221     Tag neighTgtEdgeTag;  // will store edge borders for each tgt cell
00222 
00223     Tag orgSendProcTag;  /// for coverage mesh, will store the original sender
00224 
00225     // EntityType type; // this will be tri, quad or MBPOLYGON...
00226 
00227     const EntityHandle* tgtConn;
00228     const EntityHandle* srcConn;
00229     CartVect tgtCoords[MAXEDGES];
00230     CartVect srcCoords[MAXEDGES];
00231     double tgtCoords2D[MAXEDGES2];  // these are in plane
00232     double srcCoords2D[MAXEDGES2];  // these are in plane
00233 
00234 #ifdef ENABLE_DEBUG
00235     static int dbg_1;
00236     std::ofstream mout_1[6];  // some debug files
00237 #endif
00238     // for each tgt edge, we keep a vector of extra nodes, coming from intersections
00239     // use the index in TgtEdges range
00240     // so the extra nodes on each tgt edge are kept track of
00241     // only entity handles are in the vector, not the actual coordinates;
00242     // actual coordinates are retrieved every time, which could be expensive
00243     // maybe we should store the coordinates too, along with entity handles (more memory used,
00244     // faster to retrieve)
00245     std::vector< std::vector< EntityHandle >* > extraNodesVec;
00246 
00247     double epsilon_1;
00248     double epsilon_area;
00249 
00250     std::vector< double > allBoxes;
00251     double box_error;
00252     /* \brief Local root of the kdtree */
00253     EntityHandle localRoot;
00254     Range localEnts;  // this range is for local elements of interest, euler cells, or "first mesh"
00255     unsigned int my_rank;
00256 
00257 #ifdef MOAB_HAVE_MPI
00258     ParallelComm* parcomm;
00259     TupleList* remote_cells;                       // not used anymore for communication, just a container
00260     TupleList* remote_cells_with_tracers;          // these will be used now to update tracers on remote procs
00261     std::map< int, EntityHandle > globalID_to_eh;  // needed for parallel, mostly
00262 #endif
00263     int max_edges_1;  // maximum number of edges in the lagrange set (first set, src)
00264     int max_edges_2;  // maximum number of edges in the euler set (second set, tgt)
00265     int counting;
00266 };
00267 
00268 } /* namespace moab */
00269 #endif /* INTX2MESH_HPP_ */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines