MOAB: Mesh Oriented datABase  (version 5.3.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, EntityHandle src, double* P, int& nP,
00088                                                            double& area, int markb[MAXEDGES], int markr[MAXEDGES],
00089                                                            int& nsidesSrc, int& nsidesTgt,
00090                                                            bool check_boxes_first = false ) = 0;
00091 
00092     // this is also abstract
00093     virtual ErrorCode findNodes( EntityHandle tgt, int nsTgt, EntityHandle src, int nsSrc, double* iP, int nP ) = 0;
00094 
00095     // this is also computing the area of the tgt cell in plane (gnomonic plane for sphere)
00096     // setting the local variables:
00097     // this will be done once per tgt cell, and once per localQueue for src cells
00098     //  const EntityHandle * tgtConn;
00099     // CartVect redCoords[MAXEDGES];
00100     // double redCoords2D[MAXEDGES2]; // these are in plane
00101 
00102     virtual double setup_tgt_cell( EntityHandle tgt, int& nsTgt ) = 0;
00103 
00104     virtual ErrorCode FindMaxEdgesInSet( EntityHandle eset, int& max_edges );
00105     virtual ErrorCode FindMaxEdges( EntityHandle set1,
00106                                     EntityHandle set2 );  // this needs to be called before any
00107                                                           // covering communication in parallel
00108 
00109     virtual ErrorCode createTags();
00110 
00111     ErrorCode DetermineOrderedNeighbors( EntityHandle inputSet, int max_edges, Tag& neighTag );
00112 
00113     void set_error_tolerance( double eps )
00114     {
00115         epsilon_1    = eps;
00116         epsilon_area = eps * sqrt( eps );
00117     }
00118 
00119 #ifdef MOAB_HAVE_MPI
00120     void set_parallel_comm( moab::ParallelComm* pcomm )
00121     {
00122         parcomm = pcomm;
00123     }
00124 #endif
00125     // void SetEntityType (EntityType tp) { type=tp;}
00126 
00127     // clean some memory allocated
00128     void clean();
00129 
00130     void set_box_error( double berror )
00131     {
00132         box_error = berror;
00133     }
00134 
00135     ErrorCode create_departure_mesh_2nd_alg( EntityHandle& euler_set, EntityHandle& covering_lagr_set );
00136 
00137     // in this method, used in parallel, each departure elements are already created, and at their
00138     // positions the covering_set is output, will contain the departure cells that cover the euler
00139     // set; some of these departure cells might come from different processors so the covering_set
00140     // contains some elements from lagr_set and some elements that come from other procs we need to
00141     // keep track of what processors "sent" the elements so we know were to send back the info about
00142     // the tracers masses
00143 
00144     ErrorCode create_departure_mesh_3rd_alg( EntityHandle& lagr_set, EntityHandle& covering_set );
00145 
00146     /* in this method, used in parallel, each cell from first mesh need to be sent to the
00147      * tasks whose second mesh bounding boxes intersects bounding box of the cell
00148      * this method assumes that the bounding boxes for the second mesh were computed already in a
00149      * previous method (called build_processor_euler_boxes)
00150      *
00151      * the covering_set is output,  will cover the second mesh (set) from the task;
00152      * Some of the cells in the first mesh will be sent to multiple processors; we keep track of
00153      * them using the global id of the  vertices and global ids of the cells (we do not want to
00154      * create multiple vertices with the same ids). The global id of the cells are needed just for
00155      * debugging, the cells cannot come from 2 different tasks, but the vertices can
00156      *
00157      * Right now, we use crystal router, but an improvement might be to use direct communication
00158      * (send_entities) on parallel comm
00159      *
00160      * param initial_distributed_set (IN) : the initial distribution of the first mesh (set)
00161      *
00162      * param (OUT) : the covering set in first mesh , which completely covers the second mesh set
00163      */
00164 #ifdef MOAB_HAVE_MPI
00165 
00166     virtual ErrorCode build_processor_euler_boxes( EntityHandle euler_set, Range& local_verts );
00167 #endif
00168     void correct_polygon( EntityHandle* foundIds, int& nP );
00169 #ifdef MOAB_HAVE_MPI
00170     // share vertices between the intersection target domains
00171     ErrorCode resolve_intersection_sharing();
00172 #endif
00173 #ifdef ENABLE_DEBUG
00174     void enable_debug()
00175     {
00176         dbg_1 = 1;
00177     }
00178     void disable_debug()
00179     {
00180         dbg_1 = 0;
00181     }
00182 #endif
00183 
00184 #ifdef MOAB_HAVE_TEMPESTREMAP
00185     friend class TempestRemapper;
00186 #endif
00187 
00188   protected:  // so it can be accessed in derived classes, InPlane and OnSphere
00189     Interface* mb;
00190 
00191     EntityHandle mbs1;
00192     EntityHandle mbs2;
00193     Range rs1;  // range set 1 (departure set, lagrange set, src set, manufactured set, target mesh)
00194     Range rs2;  // range set 2 (arrival set, euler set, tgt set, initial set, source mesh)
00195 
00196     EntityHandle outSet;  // will contain intersection
00197     Tag gid;              // global id tag will be used to set the parents of the intersection cell
00198 
00199     // tags used in computation, advancing front
00200     Tag TgtFlagTag;  // to mark tgt quads already considered
00201 
00202     Range TgtEdges;  //
00203 
00204     // tgt parent and src parent tags
00205     // these will be on the out mesh
00206     Tag tgtParentTag;
00207     Tag srcParentTag;
00208     Tag countTag;
00209 
00210     Tag srcNeighTag;  // will store neighbors for navigating easily in advancing front, for first
00211                       // mesh (src, target, lagrange)
00212     Tag tgtNeighTag;  // will store neighbors for navigating easily in advancing front, for second
00213                       // mesh (tgt, source, euler)
00214 
00215     Tag neighTgtEdgeTag;  // will store edge borders for each tgt cell
00216 
00217     Tag orgSendProcTag;  /// for coverage mesh, will store the original sender
00218 
00219     // EntityType type; // this will be tri, quad or MBPOLYGON...
00220 
00221     const EntityHandle* tgtConn;
00222     const EntityHandle* srcConn;
00223     CartVect tgtCoords[MAXEDGES];
00224     CartVect srcCoords[MAXEDGES];
00225     double tgtCoords2D[MAXEDGES2];  // these are in plane
00226     double srcCoords2D[MAXEDGES2];  // these are in plane
00227 
00228 #ifdef ENABLE_DEBUG
00229     static int dbg_1;
00230     std::ofstream mout_1[6];  // some debug files
00231 #endif
00232     // for each tgt edge, we keep a vector of extra nodes, coming from intersections
00233     // use the index in TgtEdges range
00234     // so the extra nodes on each tgt edge are kept track of
00235     // only entity handles are in the vector, not the actual coordinates;
00236     // actual coordinates are retrieved every time, which could be expensive
00237     // maybe we should store the coordinates too, along with entity handles (more memory used,
00238     // faster to retrieve)
00239     std::vector< std::vector< EntityHandle >* > extraNodesVec;
00240 
00241     double epsilon_1;
00242     double epsilon_area;
00243 
00244     std::vector< double > allBoxes;
00245     double box_error;
00246     /* \brief Local root of the kdtree */
00247     EntityHandle localRoot;
00248     Range localEnts;  // this range is for local elements of interest, euler cells, or "first mesh"
00249     unsigned int my_rank;
00250 
00251 #ifdef MOAB_HAVE_MPI
00252     ParallelComm* parcomm;
00253     TupleList* remote_cells;                       // not used anymore for communication, just a container
00254     TupleList* remote_cells_with_tracers;          // these will be used now to update tracers on remote procs
00255     std::map< int, EntityHandle > globalID_to_eh;  // needed for parallel, mostly
00256 #endif
00257     int max_edges_1;  // maximum number of edges in the lagrange set (first set, src)
00258     int max_edges_2;  // maximum number of edges in the euler set (second set, tgt)
00259     int counting;
00260 };
00261 
00262 } /* namespace moab */
00263 #endif /* INTX2MESH_HPP_ */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines