moab
|
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 00024 // these are intersection utils 00025 #include "CslamUtils.hpp" 00026 00027 namespace moab { 00028 00029 #define ERRORR(rval, str) \ 00030 if (MB_SUCCESS != rval) {std::cout << str << "\n"; return rval;} 00031 00032 #define ERRORV(rval, str) \ 00033 if (MB_SUCCESS != rval) {std::cout << str << "\n"; return ;} 00034 00035 #ifdef USE_MPI 00036 // forward declarations 00037 class ParallelComm; 00038 class TupleList; 00039 #endif 00040 00041 class Intx2Mesh 00042 { 00043 public: 00044 Intx2Mesh(Interface * mbimpl); 00045 virtual ~Intx2Mesh(); 00046 00047 ErrorCode intersect_meshes(EntityHandle mbs1, EntityHandle mbs2, 00048 EntityHandle & outputSet); 00049 00050 // mark could be (3 or 4, depending on type: ) no, it could go to 10 00051 // no, it will be MAXEDGES = 10 00052 // this is pure abstract, this needs to be implemented by 00053 // all derivations 00054 // the max number of intersection points could be 2*MAXEDGES 00055 // so P must be dimensioned to 4*MAXEDGES (2*2*MAXEDGES) 00056 // so, if you intersect 2 convex polygons with MAXEDGES , you will get a convex polygon 00057 // with 2*MAXEDGES, at most 00058 // will also return the number of nodes of red and blue elements 00059 virtual int computeIntersectionBetweenRedAndBlue(EntityHandle red, 00060 EntityHandle blue, double * P, int & nP, double & area, 00061 int markb[MAXEDGES], int markr[MAXEDGES], int & nsidesBlue, 00062 int & nsidesRed, bool check_boxes_first=false)=0; 00063 00064 // this is also abstract 00065 virtual int findNodes(EntityHandle red, int nsRed, EntityHandle blue, int nsBlue, 00066 double * iP, int nP)=0; 00067 00068 virtual void createTags(); 00069 ErrorCode GetOrderedNeighbors(EntityHandle set, EntityHandle quad, 00070 EntityHandle neighbors[MAXEDGES]); 00071 00072 void SetErrorTolerance(double eps) { epsilon_1=eps; epsilon_area = eps*eps/2;} 00073 00074 //void SetEntityType (EntityType tp) { type=tp;} 00075 00076 // clean some memory allocated 00077 void clean(); 00078 00079 // this will depend on the problem and element type; return true if on the border edge too 00080 virtual bool is_inside_element(double xyz[3], EntityHandle eh) = 0; 00081 void set_box_error(double berror) 00082 {box_error = berror;} 00083 00084 ErrorCode create_departure_mesh_2nd_alg(EntityHandle & euler_set, EntityHandle & covering_lagr_set); 00085 00086 // in this method, used in parallel, each departure elements are already created, and at their positions 00087 // the covering_set is output, will contain the departure cells that cover the euler set; some of these 00088 // departure cells might come from different processors 00089 // so the covering_set contains some elements from lagr_set and some elements that come from other procs 00090 // we need to keep track of what processors "sent" the elements so we know were to 00091 // send back the info about the tracers masses 00092 00093 ErrorCode create_departure_mesh_3rd_alg(EntityHandle & lagr_set, EntityHandle & covering_set); 00094 00095 ErrorCode build_processor_euler_boxes(EntityHandle euler_set, Range & local_verts); 00096 00097 void correct_polygon(EntityHandle * foundIds, int & nP); 00098 #ifdef USE_MPI 00099 ErrorCode correct_intersection_points_positions(); 00100 #endif 00101 void enable_debug() {dbg_1 = 1;} 00102 void disable_debug() {dbg_1 = 0;} 00103 protected: // so it can be accessed in derived classes, InPlane and OnSphere 00104 Interface * mb; 00105 00106 EntityHandle mbs1; 00107 EntityHandle mbs2; 00108 Range rs1;// range set 1 (departure set, lagrange set, blue set, manufactured set) 00109 Range rs2;// range set 2 (arrival set, euler set, red set, initial set) 00110 00111 EntityHandle outSet; // will contain intersection 00112 00113 // tags used in computation, advanced front 00114 Tag BlueFlagTag; // to mark blue quads already considered 00115 00116 Tag RedFlagTag; // to mark red quads already considered 00117 00118 Range RedEdges; // 00119 00120 // red parent and blue parent tags 00121 // these will be on the out mesh 00122 Tag redParentTag; 00123 Tag blueParentTag; 00124 Tag countTag; 00125 00126 //EntityType type; // this will be tri, quad or MBPOLYGON... 00127 00128 const EntityHandle * redConn; 00129 const EntityHandle * blueConn; 00130 CartVect redCoords[MAXEDGES]; 00131 CartVect blueCoords[MAXEDGES]; 00132 double redCoords2D[MAXEDGES2]; // these are in plane 00133 double blueCoords2D[MAXEDGES2]; // these are in plane 00134 00135 std::ofstream mout_1[6]; // some debug files 00136 int dbg_1; 00137 // for each red edge, we keep a vector of extra nodes, coming from intersections 00138 // use the index in RedEdges range, instead of a map, as before 00139 // std::map<EntityHandle, std::vector<EntityHandle> *> extraNodesMap; 00140 // so the extra nodes on each red edge are kept track of 00141 std::vector<std::vector<EntityHandle> *> extraNodesVec; 00142 00143 double epsilon_1; 00144 double epsilon_area; 00145 00146 std::vector<double> allBoxes; 00147 double box_error; 00148 /* \brief Local root of the kdtree */ 00149 EntityHandle localRoot; 00150 Range localEnts;// this range is for local elements of interest, euler cells 00151 unsigned int my_rank; 00152 00153 #ifdef USE_MPI 00154 ParallelComm * parcomm; 00155 TupleList * remote_cells; 00156 std::map<int, EntityHandle> globalID_to_eh;// needed for parallel, mostly 00157 #endif 00158 int max_edges; // maximum number of edges in the euler set 00159 00160 00161 00162 }; 00163 00164 } /* namespace moab */ 00165 #endif /* INTX2MESH_HPP_ */