MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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 <time.h> 00016 #include <stdlib.h> 00017 #include <stdio.h> 00018 #include <string.h> 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 * eps / 2; 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_ */