// // TETRA.H #pragma once // Added by qi jian #ifndef TETRA_H #define TETRA_H #include "node.h" #include "edge.h" #include "face.h" #include "coord.h" #include "material.h" #include "densemat.h" #include "array.h" #include "gpc.h" #include #include "triangulate.h" #include "clipper.h" #include "Constants.h" #include #include #include #include "../dgtd-performance.hpp" #define InteriorElement 1 #define BoundaryElementABC 2 #define BoundaryElementR 3 #define BoundaryElementL 4 #define BoundaryElementC 5 #define BoundaryElementD 6 class face; class Octree; // forward declare | added by qi jian class tetra{ friend class FemGrp; friend class face; friend class Octree; // Add by qi jian private: int edgeEDOF; //edge that are NOT PEC nor (PEC & PMC) int edgeHDOF; //edge that are NOT PMC nor (PEC & PMC) int faceEDOF; //faces that are NOT PEC int faceHDOF; //faces that are NOT PMC int cnt; // position of the tetra in tetArray //cnt is the tetra local number.(i.e. the number in the partiotn it belongs to, for MPI_CO use) int objNum; // number of the type of prop of the tetra int nbc[NumOfFaces]; //bc of the faces corresponding with the numbers in the .bc int LocalEDOF; //DOFS_PER_EDGE * edgeEDOF + DOFS_PER_FACE(1+2) * faceEDOF int LocalHDOF; //DOFS_PER_EDGE * edgeHDOF + DOFS_PER_FACE(1+2) * faceHDOF int LocalOffsetE; //Position where the basis functions of each tetra for E field would start whithout all the basis functions corresponding to PEC edges or surfaces int LocalOffsetH; //Position where the basis functions of each tetra for H field would start whithout all the basis functions corresponding to PMC edges or surfaces int TetrahedronFlag; //Is ABC tetra (1 = yes) bool ExcitationFlag; //is Planwave or port int PolyOrderFlag; int MinimumPoly; int LTS_Flag; //Class number of the tetra bool isNonConformal; //Is the tetra on a Non Conformal BC? (true if yes) int regularGroup; //regular group it belongs int nDipoles; // number of dipoles in the tetra; int NeighNum; //number of neighbours int NeighNCNum; //number of non-conforming neighs int Conduct_Flag; //Unles sigma[0][0] == 0, Conduct_Flag = 1 int* mapPEC; int* mapPMC; int cntpec_; int cntpmc_; /*Order in LocMapE/LocMapH: - Order0: Edge0 Edge1 Edge2 Edge3 Edge4 Edge5 - Order1: Edge0 Edge1 Edge2 Edge3 Edge4 Edge5 - Order2: Face0_1 Face0_2 Face1_1 Face1_2 Face2_1 Face2_2 Face3_1 Face3_2 - Edge0 Edge1 Edge2 Edge3 Edge4 Edge5 - Face0_3 Face1_3 Face2_3 Face3_3 The functions are hierarchical so the smaller oreders are included in the higher ones */ int* LocMapE; //Array with numbers in the positions where the edge/face is not pec, if it's PEC the value in this position is -1 int* LocMapH; //Array with numbers in the positions where the edge/face is not pmc, if it's PMC the value in this position is -1 fp_t Class_dt; //Time step in the class where the tetra belongs fp_t Stability_dt; //Time step of the tetra without LTS node* nd[NumOfNodes]; //nodes of the tetra edge* eg[NumOfEdges]; //edges of the tetra ordered by getGlobalCnt face* fc[NumOfFaces]; //faces of the tetra tetra** NeighTetra; //neigh tetras for this element Nullpnt in an element if the face just have 1 tetra (Also include non-conforming ones) int* NeighTetraFaceMap; //Faces that correspond to each neighbour tetra** NeighNCTetra; //neigh tetras that are non-conforming denseMat* NC_CouplingMatArray; // /////// denseMat* MassE; denseMat* MassH; denseMat* StiffE; denseMat* StiffH; denseMat* BiiE; denseMat* BiiH; denseMat* FiiE; denseMat* FiiH; denseMat* FijE; denseMat* FijH; denseMat* BijE; denseMat* BijH; denseMat* LHSInvE; denseMat* LHSInvH; denseMat* RHSLoc1E; denseMat* RHSLoc1H; denseMat* RHSLoc2E; denseMat* RHSLoc2H; denseMat* RHSCoup1E; denseMat* RHSCoup1H; denseMat* RHSCoup2E; denseMat* RHSCoup2H; vector> OverlapMapVec; vector> IntersectionMapVec; vector> IntersectionFaceMapVec; Material* mat; //type of prop of the tetra fp_t GAMMA = 1; // Gamma(penalty) -> 1 upwind , 0 central fp_t vol_; //volume of the tetra vtr gradZeta_[4]; // inward-pointing face normal vectors with magnitude equal to face area // Non-matrix free,ddm variables int *EdofMap; int *HdofMap; // PML Variables bool PML_Flag; denseMat* matA; denseMat* matB; denseMat* matC; denseMat* matD; denseMat* matF; denseMat* MassI_E; denseMat* ABC_tetBe; denseMat* Mass_epA_E; denseMat* Mass_epB_E; denseMat* Mass_epC_E; denseMat* Mass_D_E; denseMat* Mass_F_E; denseMat* MassI_H; denseMat* ABC_tetBh; denseMat* Mass_muA_H; denseMat* Mass_muB_H; denseMat* Mass_muC_H; denseMat* Mass_D_H; denseMat* Mass_F_H; // In Main Equation denseMat* RHSAuxE; denseMat* RHSAuxH; // In ADE denseMat* LHSInvAuxM; denseMat* RHSLocAuxM1; denseMat* RHSLocAuxM2; denseMat* LHSInvAuxJ; denseMat* RHSLocAuxJ1; denseMat* RHSLocAuxJ2; fp_t sigma_1, sigma_2, sigma_3, max_conductivity; bool scattering_region; // Flag to indicate if the tetra is in a scattering region denseMat* Jrot; public: // constructor tetra(int = 0); ~tetra(); // operator tetra& operator = (const tetra &); int operator==(const tetra& right)const; // set functions void set_objNum(const int); void set_nbc(const int, const int, const int, const int); void set_node(node *, node *, node *, node *); void setEdge(edge *, int); void setFace(face *, int); void reArrange(); void setcnt(int nn){cnt = nn;} void setIsNC(bool isNC){isNonConformal = isNC;} void setRegularGroup(int n){regularGroup = n;} void set_nDipoles(int value){nDipoles = value;} void increaseNDipoles(){nDipoles++;} // Added for DG void set_mat(Material* m){mat = m;} void set_TetrahedronFlag(); void set_ConductivityFlag(); void set_PolyOrderFlagDebug(int n); void set_flux_GAMMA(const fp_t GAMMA_){GAMMA = GAMMA_;} void set_MinimumPoly(int n){MinimumPoly = n;} void set_LTS_Flag(int n){LTS_Flag = n;} void set_Class_dt(fp_t dtval){Class_dt = dtval;} void set_Stability_dt(fp_t dtval){Stability_dt = dtval;} //void set_NeighborTetra(tetra* tetARRAY, int* ncARRAY, int nonConformalCNT); // Modified by qi jian void set_NeighborTetra(tetra* tetARRAY, int* ncARRAY, int nonConformalCNT, Octree* octree_ptr, double min_AABB_size); void set_LocalOffsetE(int n){LocalOffsetE = n;} void set_LocalOffsetH(int n){LocalOffsetH = n;} void SetFacePEC(); void SetFacePMC(); void set_cntpec_(int value){cntpec_ = value;} void set_cntpmc_(int value){cntpmc_ = value;} void set_einc(fp_t, ArrayFP *, ArrayFP *); void set_hinc(fp_t, ArrayFP *, ArrayFP *); // get functions int getobjNum(){return objNum;} int getbc(int n){return nbc[n];} bool getIsNC(){return isNonConformal;} int getRegularGroup(){return regularGroup;} int get_nDipoles(){return nDipoles;} int getEdgeID(int); int getFaceID(int); int getcnt(){return cnt;} int getOrientation(); int get_ExcitationFlag(){return ExcitationFlag;} int get_TetrahedronFlag(){return TetrahedronFlag;} int get_PolyOrderFlag(){return PolyOrderFlag;} int get_MinimumPoly(){return MinimumPoly;} int get_ConductivityFlag(){return Conduct_Flag;} edge* getEdgePtr(int); face* getFacePtr(int); node* getNode(int n) {return nd[n];} Material* getMat(){return mat;} tetra* get_NeighborTetra(int n) {return NeighTetra[n];} int get_LocalOffsetE(){return LocalOffsetE;} int get_LocalOffsetH(){return LocalOffsetH;} int get_LTS_Flag(){return LTS_Flag;} fp_t get_Class_dt(){return Class_dt;} fp_t get_Stability_dt(){return Stability_dt;} int get_cntpec_(){return cntpec_;} int get_cntpmc_(){return cntpmc_;} int get_LocMapEntry_E(int l){return LocMapE[l];} int get_LocMapEntry_H(int l){return LocMapH[l];} int get_NeighNum(){return NeighNum;} void TetraBasisW(fp_t *zeta, int p, vtr *W); void TetraBasisCurlW(fp_t *zeta, int p, vtr *W); vtr get_centroid(); //NonConformal Operations bool triangle_overlaps_triangle(face* LocalFace, face* NeighFace); bool triangle_overlaps_triangleCO(face* LocalFace, face* NeighFace); // Elementrary Matrices procedures void MapPEC(); void MapPMC(); void CountDOF_E(); void CountDOF_H(); void Local_DG_mapE(int *, int); void Local_DG_mapH(int *, int); void geometry(vtr *, vtr *, fp_t *); void makeCoeff(vtr *, tensor, fp_t [][3]); void makeCoeff_curl(vtr *, fp_t [][3], fp_t& scale); // Added for DG void Calculate_M_Matrix_E(); void Calculate_R_Matrix_E(denseMat*, fp_t); void Calculate_Bii_Matrix_E(); void Calculate_Bij_Matrix_E(); void Calculate_S_Matrix_E(); void Calculate_Fii_Matrix_E(); void Calculate_Fij_Matrix_E(); void SetUp_LocalFaceToTetraMapE_NMF1(fp_t dt); void Calculate_M_Matrix_E_Numeric(); void Calculate_R_Matrix_E_Numeric(denseMat*, fp_t); void Calculate_Bii_Matrix_E_Numeric(); void Calculate_Bij_Matrix_E_Numeric(); void Calculate_S_Matrix_E_Numeric(); void Calculate_Fii_Matrix_E_Numeric(); void Calculate_Fij_Matrix_E_Numeric(); void SetUp_LocalFaceToTetraMapE_NMF1_Numeric(fp_t dt); void Calculate_M_Matrix_H(); void Calculate_R_Matrix_H(denseMat*, fp_t); void Calculate_Bii_Matrix_H(); void Calculate_Bij_Matrix_H(); void Calculate_S_Matrix_H(); void Calculate_Fii_Matrix_H(); void Calculate_Fij_Matrix_H(); void SetUp_LocalFaceToTetraMapH_NMF1(fp_t dt); void Calculate_M_Matrix_H_Numeric(); void Calculate_R_Matrix_H_Numeric(denseMat*, fp_t); void Calculate_Bii_Matrix_H_Numeric(); void Calculate_Bij_Matrix_H_Numeric(); void Calculate_S_Matrix_H_Numeric(); void Calculate_Fii_Matrix_H_Numeric(); void Calculate_Fij_Matrix_H_Numeric(); void SetUp_LocalFaceToTetraMapH_NMF1_Numeric(fp_t dt); void Calculate_S_Matrix(denseMat *); void Calculate_S_Matrix_Numeric(denseMat *); void ApplyPEC_Vector(ArrayFP*, int); void ApplyPMC_Vector(ArrayFP*, int); void ApplyPEC(denseMat*, int); void ApplyPMC(denseMat*, int); void ApplyPEC_E(denseMat*, int); void ApplyPMC_E(denseMat*, int); void ApplyPEC_H(denseMat*, int); void ApplyPMC_H(denseMat*, int); void ApplyPMC_Neigh_E(denseMat*, tetra*, int); void ApplyPEC_Neigh_H(denseMat*, tetra*, int); void FreeNonExcitatonTetra(); // Time-Marching procedures void TimeStepEstimate (fp_t &dtEstim, fp_t& V_P); // DEVICE int* get_NeighMapsE (); int* get_NeighMapsH (); // Matrix Free Implementation void SetUpMatrixFree(); void LocalFaceToTetraMapE_NMF1(ArrayFP& En_1, ArrayFP& En, ArrayFP& Hn_12, fp_t dt, fp_t t); void LocalFaceToTetraMapH_NMF1(ArrayFP& Hn_32, ArrayFP& En_1, ArrayFP& Hn_12, fp_t dt, fp_t t); void Excitation_E_BLAS(ArrayFP& En_1, fp_t t, fp_t dt, ArrayFP&); void Excitation_H_BLAS(ArrayFP& Hn_32, fp_t t, fp_t dt, ArrayFP&); // print void print(); vtr Center(){ vtr vv = (nd[0]->getCoord() + nd[1]->getCoord() + nd[2]->getCoord() + nd[3]->getCoord()) * (1. / 4.); return vv; } void setEHGlobalMap(int i,int ide,int idh); void allocDofMap(); void freeDofMap(); int* get_LocalEMap(){return LocMapE;} int* get_LocalHMap(){return LocMapH;} int get_LocalEdof(){return LocalEDOF;} int get_LocalHdof(){return LocalHDOF;} void EstimateFLOs(int numLTSClass); bool essentiallyEqual(double a, double b, double epsilon); void SetUpInterfaceMatrices_NC(); void SetUpInterfaceMaps_NC(); bool CheckIntersection(face* LocalFace, face* NeighFace); face* CreateInterSectionMesh(face* LocalFace, face* NeighFace, int& tcount); int getNeighFace(tetra* neigh); void prepare_GPU(fp_t_ts* LHS_InvE_init, fp_t_ts* LHS_InvH_init, fp_t_ts* LE_Matrices, fp_t_ts* LH_Matrices, bool needInverse); void prepareCuBLAS(fp_t_ts* Loc1E_h, fp_t_ts* Loc2E_h, fp_t_ts* Neigh1E_h, fp_t_ts* Neigh2E_h, fp_t_ts* InvE_h, fp_t_ts* Loc1H_h, fp_t_ts* Loc2H_h, fp_t_ts* Neigh1H_h, fp_t_ts* Neigh2H_h, fp_t_ts* InvH_h); void invxCol(fp_t* col, bool Efield); // PML Functions void set_PML_Flag(bool flag) { PML_Flag = flag; } bool get_PML_Flag() const { return PML_Flag; } void prepareCuBLAS_PML(fp_t_ts* Loc1E_h, fp_t_ts* Loc2E_h, fp_t_ts* Neigh1E_h, fp_t_ts* Neigh2E_h, fp_t_ts* Loc1H_h, fp_t_ts* Loc2H_h, fp_t_ts* Neigh1H_h, fp_t_ts* Neigh2H_h, fp_t_ts* AuxE_h, fp_t_ts* AuxH_h, fp_t_ts* AuxM1_h, fp_t_ts* AuxM2_h, fp_t_ts* AuxJ1_h, fp_t_ts* AuxJ2_h); void Calculate_M_Matrix_I_E_Numeric(); void Calculate_ABC_E_Numeric(); void Calculate_M_Matrix_I_H_Numeric(); void Calculate_ABC_H_Numeric(); void set_Conductivity_Profile_Planar(fp_t planewave_xmin, fp_t planewave_ymin, fp_t planewave_zmin, fp_t planewave_xmax, fp_t planewave_ymax, fp_t planewave_zmax); void set_Conductivity_Profile_Conformal(fp_t Rx, fp_t Ry, fp_t Rz); void SetUp_LocalFaceToTetraMapE_NMF1_PML(fp_t dt); void SetUp_LocalFaceToTetraMapH_NMF1_PML(fp_t dt); void Calculate_M_Matrix_I_E(); void Calculate_ABC_E(); void Calculate_M_Matrix_I_H(); void Calculate_ABC_H(); void Calculate_Mass_Material_Matrix( denseMat*& outputMassMat, const denseMat* inputTensorMat, const tensor& materialTensor, bool isElectricField ); void Calculate_Mass_Material_Matrix_Numeric( denseMat*& outputMassMat, const denseMat* inputTensorMat, const tensor& materialTensor, bool isElectricField ); void set_Conductivity_Profile_Spherical(fp_t Rx, fp_t Ry, fp_t Rz); void set_Conductivity_smoothProfile_Planar(fp_t planewave_xmin, fp_t planewave_ymin, fp_t planewave_zmin, fp_t planewave_xmax, fp_t planewave_ymax, fp_t planewave_zmax, vtr rc, vtr gaussPointCoord, string mode); void getCartesianCoord(fp_t *zeta, vtr &pnt); void Calculate_Mass_Material_Matrix_Vary_Numeric( denseMat*& outputMassMat, // Output mass matrix to fill const denseMat* inputTensorMat, // Input tensor like A1, A3 const tensor& materialTensor, // Material tensor (e.g., epsr, mur) bool isElectricField, // true → use E-field space, false → H-field space string mode, fp_t planewave_xmin, fp_t planewave_ymin, fp_t planewave_zmin, fp_t planewave_xmax, fp_t planewave_ymax, fp_t planewave_zmax, fp_t radius_x, fp_t radius_y, fp_t radius_z ); void LocalRadiiCurvature(vtr& r, fp_t radius_x, fp_t radius_y, fp_t radius_z, fp_t& r01, fp_t& r02, vtr& u1, vtr& u2, vtr& u3); void ApplyConformalPMLRotation(denseMat* J, denseMat* tensorU1U2U3); void set_Conductivity_smoothProfile_conformal(fp_t radius_x, fp_t radius_y, fp_t radius_z, vtr rc, vtr gaussPointCoord, string mode); void Et_S_E(denseMat* E, denseMat* S, denseMat* EtSE); }; #endif