1 /* 2 * Intx2Mesh.hpp 3 * 4 * Created on: Oct 2, 2012 5 * 6 */ 7 8 #ifndef INTX2MESH_HPP_ 9 #define INTX2MESH_HPP_ 10 11 #include <iostream> 12 #include <sstream> 13 #include <fstream> 14 #include <map> 15 #include <time.h> 16 #include <stdlib.h> 17 #include <stdio.h> 18 #include <string.h> 19 #include "moab/Core.hpp" 20 #include "moab/Interface.hpp" 21 #include "moab/Range.hpp" 22 #include "moab/CartVect.hpp" 23 24 // #define ENABLE_DEBUG 25 26 // maximum number of edges on each convex polygon of interest 27 #define MAXEDGES 10 28 #define MAXEDGES2 20 // used for coordinates in plane 29 #define CORRTAGNAME "__correspondent" 30 31 namespace moab { 32 33 #define ERRORR(rval, str) \ 34 if (MB_SUCCESS != rval) \ 35 {std::cout << str << "\n"; return rval;} 36 37 #define ERRORV(rval, str) \ 38 if (MB_SUCCESS != rval) \ 39 {std::cout << str << "\n"; return ;} 40 41 #ifdef MOAB_HAVE_MPI 42 // forward declarations 43 class ParallelComm; 44 class TupleList; 45 #endif 46 47 class Intx2Mesh 48 { 49 public: 50 Intx2Mesh(Interface * mbimpl); 51 virtual ~Intx2Mesh(); 52 53 /* 54 * computes intersection between 2 sets; 55 * set 2 (mbs2) should be completely covered by set mbs1, as all elements 56 * from set 2 will be decomposed in polygons ; an area verification is done, and 57 * will signal elements from set 2 that are not "recovered" 58 * 59 * the resulting polygons are inserted in output set; total area from output set should be 60 * equal to total area of elements in set 2 (check that!) 61 */ 62 ErrorCode intersect_meshes(EntityHandle mbs1, EntityHandle mbs2, 63 EntityHandle & outputSet); 64 65 // mark could be (3 or 4, depending on type: ) no, it could go to 10 66 // no, it will be MAXEDGES = 10 67 // this is pure abstract, this needs to be implemented by 68 // all derivations 69 // the max number of intersection points could be 2*MAXEDGES 70 // so P must be dimensioned to 4*MAXEDGES (2*2*MAXEDGES) 71 // so, if you intersect 2 convex polygons with MAXEDGES , you will get a convex polygon 72 // with 2*MAXEDGES, at most 73 // will also return the number of nodes of red and blue elements 74 virtual ErrorCode computeIntersectionBetweenRedAndBlue(EntityHandle red, 75 EntityHandle blue, double * P, int & nP, double & area, 76 int markb[MAXEDGES], int markr[MAXEDGES], int & nsidesBlue, 77 int & nsidesRed, bool check_boxes_first=false)=0; 78 79 // this is also abstract 80 virtual ErrorCode findNodes(EntityHandle red, int nsRed, EntityHandle blue, int nsBlue, 81 double * iP, int nP)=0; 82 83 // this is also computing the area of the red cell in plane (gnomonic plane for sphere) 84 // setting the local variables: 85 // this will be done once per red cell, and once per localQueue for blue cells 86 // const EntityHandle * redConn; 87 // CartVect redCoords[MAXEDGES]; 88 // double redCoords2D[MAXEDGES2]; // these are in plane 89 90 virtual double setup_red_cell(EntityHandle red, int & nsRed)= 0; 91 92 virtual ErrorCode FindMaxEdgesInSet(EntityHandle eset, int & max_edges); 93 virtual ErrorCode FindMaxEdges(EntityHandle set1, EntityHandle set2); // this needs to be called before any covering communication in parallel 94 95 virtual ErrorCode createTags(); 96 97 ErrorCode DetermineOrderedNeighbors(EntityHandle inputSet, int max_edges, Tag & neighTag); 98 set_error_tolerance(double eps)99 void set_error_tolerance(double eps) { epsilon_1=eps; epsilon_area = eps*eps/2;} 100 101 #ifdef MOAB_HAVE_MPI set_parallel_comm(moab::ParallelComm * pcomm)102 void set_parallel_comm(moab::ParallelComm* pcomm) { parcomm = pcomm; } 103 #endif 104 //void SetEntityType (EntityType tp) { type=tp;} 105 106 // clean some memory allocated 107 void clean(); 108 set_box_error(double berror)109 void set_box_error(double berror) 110 {box_error = berror;} 111 112 ErrorCode create_departure_mesh_2nd_alg(EntityHandle & euler_set, EntityHandle & covering_lagr_set); 113 114 // in this method, used in parallel, each departure elements are already created, and at their positions 115 // the covering_set is output, will contain the departure cells that cover the euler set; some of these 116 // departure cells might come from different processors 117 // so the covering_set contains some elements from lagr_set and some elements that come from other procs 118 // we need to keep track of what processors "sent" the elements so we know were to 119 // send back the info about the tracers masses 120 121 ErrorCode create_departure_mesh_3rd_alg(EntityHandle & lagr_set, EntityHandle & covering_set); 122 123 /* in this method, used in parallel, each cell from first mesh need to be sent to the 124 * tasks whose second mesh bounding boxes intersects bounding box of the cell 125 * this method assumes that the bounding boxes for the second mesh were computed already in a previous method 126 * (called build_processor_euler_boxes) 127 * 128 * the covering_set is output, will cover the second mesh (set) from the task; 129 * Some of the cells in the first mesh will be sent to multiple processors; we keep track of them using the global id 130 * of the vertices and global ids of the cells (we do not want to create multiple vertices with the 131 * same ids). The global id of the cells are needed just for debugging, the cells cannot come from 2 different 132 * tasks, but the vertices can 133 * 134 * Right now, we use crystal router, but an improvement might be to use direct communication (send_entities) 135 * on parallel comm 136 * 137 * param initial_distributed_set (IN) : the initial distribution of the first mesh (set) 138 * 139 * param (OUT) : the covering set in first mesh , which completely covers the second mesh set 140 */ 141 #ifdef MOAB_HAVE_MPI 142 143 virtual ErrorCode build_processor_euler_boxes(EntityHandle euler_set, Range & local_verts); 144 #endif 145 void correct_polygon(EntityHandle * foundIds, int & nP); 146 #ifdef MOAB_HAVE_MPI 147 ErrorCode correct_intersection_points_positions(); 148 #endif 149 #ifdef ENABLE_DEBUG enable_debug()150 void enable_debug() {dbg_1 = 1;} disable_debug()151 void disable_debug() {dbg_1 = 0;} 152 #endif 153 protected: // so it can be accessed in derived classes, InPlane and OnSphere 154 Interface * mb; 155 156 EntityHandle mbs1; 157 EntityHandle mbs2; 158 Range rs1;// range set 1 (departure set, lagrange set, blue set, manufactured set, target mesh) 159 Range rs2;// range set 2 (arrival set, euler set, red set, initial set, source mesh) 160 161 EntityHandle outSet; // will contain intersection 162 Tag gid; // global id tag will be used to set the parents of the intersection cell 163 164 // tags used in computation, advancing front 165 Tag RedFlagTag; // to mark red quads already considered 166 167 Range RedEdges; // 168 169 // red parent and blue parent tags 170 // these will be on the out mesh 171 Tag redParentTag; 172 Tag blueParentTag; 173 Tag countTag; 174 175 Tag blueNeighTag; // will store neighbors for navigating easily in advancing front, for first mesh (blue, target, lagrange) 176 Tag redNeighTag; // will store neighbors for navigating easily in advancing front, for second mesh (red, source, euler) 177 178 Tag neighRedEdgeTag; // will store edge borders for each red cell 179 180 Tag orgSendProcTag; /// for coverage mesh, will store the original sender 181 182 //EntityType type; // this will be tri, quad or MBPOLYGON... 183 184 const EntityHandle * redConn; 185 const EntityHandle * blueConn; 186 CartVect redCoords[MAXEDGES]; 187 CartVect blueCoords[MAXEDGES]; 188 double redCoords2D[MAXEDGES2]; // these are in plane 189 double blueCoords2D[MAXEDGES2]; // these are in plane 190 191 #ifdef ENABLE_DEBUG 192 static int dbg_1; 193 std::ofstream mout_1[6]; // some debug files 194 #endif 195 // for each red edge, we keep a vector of extra nodes, coming from intersections 196 // use the index in RedEdges range 197 // so the extra nodes on each red edge are kept track of 198 // only entity handles are in the vector, not the actual coordinates; 199 // actual coordinates are retrieved every time, which could be expensive 200 // maybe we should store the coordinates too, along with entity handles (more memory used, faster to retrieve) 201 std::vector<std::vector<EntityHandle> *> extraNodesVec; 202 203 double epsilon_1; 204 double epsilon_area; 205 206 std::vector<double> allBoxes; 207 double box_error; 208 /* \brief Local root of the kdtree */ 209 EntityHandle localRoot; 210 Range localEnts;// this range is for local elements of interest, euler cells, or "first mesh" 211 unsigned int my_rank; 212 213 #ifdef MOAB_HAVE_MPI 214 ParallelComm * parcomm; 215 TupleList * remote_cells; // not used anymore for communication, just a container 216 TupleList * remote_cells_with_tracers; // these will be used now to update tracers on remote procs 217 std::map<int, EntityHandle> globalID_to_eh;// needed for parallel, mostly 218 #endif 219 int max_edges_1; // maximum number of edges in the lagrange set (first set, blue) 220 int max_edges_2; // maximum number of edges in the euler set (second set, red) 221 int counting; 222 223 }; 224 225 } /* namespace moab */ 226 #endif /* INTX2MESH_HPP_ */ 227