You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
447 lines
17 KiB
447 lines
17 KiB
// // 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 <assert.h>
|
|
#include "triangulate.h"
|
|
#include "clipper.h"
|
|
#include "Constants.h"
|
|
#include <set>
|
|
#include <map>
|
|
#include <algorithm>
|
|
#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<fp_t>* NC_CouplingMatArray;
|
|
|
|
// ///////
|
|
denseMat<fp_t>* MassE;
|
|
denseMat<fp_t>* MassH;
|
|
denseMat<fp_t>* StiffE;
|
|
denseMat<fp_t>* StiffH;
|
|
denseMat<fp_t>* BiiE;
|
|
denseMat<fp_t>* BiiH;
|
|
denseMat<fp_t>* FiiE;
|
|
denseMat<fp_t>* FiiH;
|
|
denseMat<fp_t>* FijE;
|
|
denseMat<fp_t>* FijH;
|
|
denseMat<fp_t>* BijE;
|
|
denseMat<fp_t>* BijH;
|
|
|
|
denseMat<fp_t>* LHSInvE;
|
|
denseMat<fp_t>* LHSInvH;
|
|
denseMat<fp_t>* RHSLoc1E;
|
|
denseMat<fp_t>* RHSLoc1H;
|
|
denseMat<fp_t>* RHSLoc2E;
|
|
denseMat<fp_t>* RHSLoc2H;
|
|
denseMat<fp_t>* RHSCoup1E;
|
|
denseMat<fp_t>* RHSCoup1H;
|
|
denseMat<fp_t>* RHSCoup2E;
|
|
denseMat<fp_t>* RHSCoup2H;
|
|
|
|
vector<map<int, bool>> OverlapMapVec;
|
|
vector<map<int, bool>> IntersectionMapVec;
|
|
vector<map<int, int>> 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<fp_t>* matA;
|
|
denseMat<fp_t>* matB;
|
|
denseMat<fp_t>* matC;
|
|
denseMat<fp_t>* matD;
|
|
denseMat<fp_t>* matF;
|
|
|
|
|
|
|
|
denseMat<fp_t>* MassI_E;
|
|
denseMat<fp_t>* ABC_tetBe;
|
|
|
|
denseMat<fp_t>* Mass_epA_E;
|
|
denseMat<fp_t>* Mass_epB_E;
|
|
denseMat<fp_t>* Mass_epC_E;
|
|
denseMat<fp_t>* Mass_D_E;
|
|
denseMat<fp_t>* Mass_F_E;
|
|
|
|
|
|
denseMat<fp_t>* MassI_H;
|
|
denseMat<fp_t>* ABC_tetBh;
|
|
|
|
denseMat<fp_t>* Mass_muA_H;
|
|
denseMat<fp_t>* Mass_muB_H;
|
|
denseMat<fp_t>* Mass_muC_H;
|
|
denseMat<fp_t>* Mass_D_H;
|
|
denseMat<fp_t>* Mass_F_H;
|
|
|
|
|
|
// In Main Equation
|
|
denseMat<fp_t>* RHSAuxE;
|
|
denseMat<fp_t>* RHSAuxH;
|
|
|
|
// In ADE
|
|
denseMat<fp_t>* LHSInvAuxM;
|
|
denseMat<fp_t>* RHSLocAuxM1;
|
|
denseMat<fp_t>* RHSLocAuxM2;
|
|
|
|
denseMat<fp_t>* LHSInvAuxJ;
|
|
denseMat<fp_t>* RHSLocAuxJ1;
|
|
denseMat<fp_t>* 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<fp_t>* 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<fp_t> *, ArrayFP<fp_t> *);
|
|
void set_hinc(fp_t, ArrayFP<fp_t> *, ArrayFP<fp_t> *);
|
|
|
|
// 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>*, 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>*, 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>*, 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>*, 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<fp_t> *);
|
|
void Calculate_S_Matrix_Numeric(denseMat<fp_t> *);
|
|
|
|
void ApplyPEC_Vector(ArrayFP<fp_t>*, int);
|
|
void ApplyPMC_Vector(ArrayFP<fp_t>*, int);
|
|
void ApplyPEC(denseMat<fp_t>*, int);
|
|
void ApplyPMC(denseMat<fp_t>*, int);
|
|
void ApplyPEC_E(denseMat<fp_t>*, int);
|
|
void ApplyPMC_E(denseMat<fp_t>*, int);
|
|
void ApplyPEC_H(denseMat<fp_t>*, int);
|
|
void ApplyPMC_H(denseMat<fp_t>*, int);
|
|
void ApplyPMC_Neigh_E(denseMat<fp_t>*, tetra*, int);
|
|
void ApplyPEC_Neigh_H(denseMat<fp_t>*, 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<fp_t>& En_1, ArrayFP<fp_t>& En, ArrayFP<fp_t>& Hn_12, fp_t dt, fp_t t);
|
|
void LocalFaceToTetraMapH_NMF1(ArrayFP<fp_t>& Hn_32, ArrayFP<fp_t>& En_1, ArrayFP<fp_t>& Hn_12, fp_t dt, fp_t t);
|
|
|
|
void Excitation_E_BLAS(ArrayFP<fp_t>& En_1, fp_t t, fp_t dt, ArrayFP<fp_t>&);
|
|
void Excitation_H_BLAS(ArrayFP<fp_t>& Hn_32, fp_t t, fp_t dt, ArrayFP<fp_t>&);
|
|
|
|
// 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<fp_t>*& outputMassMat,
|
|
const denseMat<fp_t>* inputTensorMat,
|
|
const tensor& materialTensor,
|
|
bool isElectricField
|
|
);
|
|
void Calculate_Mass_Material_Matrix_Numeric(
|
|
denseMat<fp_t>*& outputMassMat,
|
|
const denseMat<fp_t>* 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<fp_t>*& outputMassMat, // Output mass matrix to fill
|
|
const denseMat<fp_t>* 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<fp_t>* J, denseMat<fp_t>* 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<fp_t>* E, denseMat<fp_t>* S, denseMat<fp_t>* EtSE);
|
|
|
|
|
|
};
|
|
|
|
#endif
|
|
|
|
|