#ifndef _SCTL_BOUNDARY_INTEGRAL_HPP_ #define _SCTL_BOUNDARY_INTEGRAL_HPP_ #include #include SCTL_INCLUDE(comm.hpp) #include SCTL_INCLUDE(mem_mgr.hpp) #include SCTL_INCLUDE(vector.hpp) #include #include namespace SCTL_NAMESPACE { template class Matrix; template class ParticleFMM; /** * Abstract base class for an element-list. In addition to the functions * declared in this base class, the derived class must be copy-constructible. */ template class ElementListBase { public: virtual ~ElementListBase() {} /** * Returns the position and normals of the surface nodal points for each * element. * * @param[out] X the position of the node points in array-of-struct * order: {x_1, y_1, z_1, x_2, ..., x_n, y_n, z_n} * * @param[out] Xn the normal vectors of the node points in * array-of-struct order: {nx_1, ny_1, nz_1, nx_2, ..., nx_n, ny_n, nz_n} * * @param[out] element_wise_node_cnt the number of node points * belonging to each element. */ virtual void GetNodeCoord(Vector* X, Vector* Xn, Vector* element_wise_node_cnt) const = 0; /** * Given an accuracy tolerance, returns the quadrature node positions, * the normals at the nodes, the weights and the cut-off distance from * the nodes for computing the far-field potential from the surface (at * target points beyond the cut-off distance). * * @param[out] X the position of the quadrature node points in * array-of-struct order: {x_1, y_1, z_1, x_2, ..., x_n, y_n, z_n} * * @param[out] Xn the normal vectors at the quadrature node points in * array-of-struct order: {nx_1, ny_1, nz_1, nx_2, ..., nx_n, ny_n, nz_n} * * @param[out] wts the weights corresponding to each quadrature node. * * @param[out] dist_far the cut-off distance from each quadrature node * such that quadrature rule will be accurate to the specified tolerance * for target points further away than this distance. * * @param[out] element_wise_node_cnt the number of quadrature node * points belonging to each element. * * @param[in] tol the accuracy tolerance. */ virtual void GetFarFieldNodes(Vector& X, Vector& Xn, Vector& wts, Vector& dist_far, Vector& element_wise_node_cnt, const Real tol) const = 0; /** * Interpolates the density from surface node points to far-field * quadrature node points. * * @param[out] Fout the interpolated density at far-field quadrature * nodes in array-of-struct order. * * @param[in] Fin the input density at surface node points in * array-of-struct order. */ virtual void GetFarFieldDensity(Vector& Fout, const Vector& Fin) const = 0; /** * Apply the transpose of the GetFarFieldDensity() operator for a single * element applied to the column-vectors of Min and the result is * returned in Mout. * * @param[out] Mout the output matrix where the column-vectors are the * result of the application of the transpose operator. * * @param[in] Min the input matrix whose column-vectors are * multiplied by the transpose operator. * * @param[in] elem_idx the index of the element. */ virtual void FarFieldDensityOperatorTranspose(Matrix& Mout, const Matrix& Min, const Long elem_idx) const = 0; /** * Compute self-interaction operator for each element. * * @param[out] M_lst the vector of all self-interaction matrices * (in row-major format). * * @param[in] ker the kernel object. * * @param[in] tol the accuracy tolerance. * * @param[in] trg_dot_prod whether to compute dot product of the potential with the target-normal vector. * * @param[in] self pointer to element-list object. */ template static void SelfInterac(Vector>& M_lst, const Kernel& ker, Real tol, bool trg_dot_prod, const ElementListBase* self); /** * Compute near-interaction operator for a given element-idx and each each target. * * @param[out] M the near-interaction matrix (in row-major format). * * @param[in] Xt the position of the target points in array-of-structure * order: {x_1, y_1, z_1, x_2, ..., x_n, y_n, z_n} * * @param[in] normal_trg the normal at the target points in array-of-structure * order: {nx_1, ny_1, nz_1, nx_2, ..., nx_n, ny_n, nz_n} * * @param[in] ker the kernel object. * * @param[in] tol the accuracy tolerance. * * @param[in] elem_idx the index of the source element. * * @param[in] self pointer to element-list object. */ template static void NearInterac(Matrix& M, const Vector& Xt, const Vector& normal_trg, const Kernel& ker, Real tol, const Long elem_idx, const ElementListBase* self); }; /** * Implements parallel computation of boundary integrals. */ template class BoundaryIntegralOp { static constexpr Integer KDIM0 = Kernel::SrcDim(); static constexpr Integer KDIM1 = Kernel::TrgDim(); static constexpr Integer COORD_DIM = 3; public: BoundaryIntegralOp() = delete; BoundaryIntegralOp(const BoundaryIntegralOp&) = delete; BoundaryIntegralOp& operator= (const BoundaryIntegralOp&) = delete; /** * Constructor * * @param[in] ker the kernel object. * * @param[in] comm the MPI communicator. */ explicit BoundaryIntegralOp(const Kernel& ker, bool trg_normal_dot_prod = false, const Comm& comm = Comm::Self()); /** * Destructor */ ~BoundaryIntegralOp(); /** * Specify accuracy tolerance. */ void SetAccuracy(Real tol); /** * Set kernel functions for FMM translation operators */ template void SetFMMKer(const KerS2M& k_s2m, const KerS2L& k_s2l, const KerS2T& k_s2t, const KerM2M& k_m2m, const KerM2L& k_m2l, const KerM2T& k_m2t, const KerL2L& k_l2l, const KerL2T& k_l2t); /** * Add an element-list. */ template void AddElemList(const ElemLstType& elem_lst, const std::string& name = std::to_string(typeid(ElemLstType).hash_code())); /** * Get const reference to an element-list. */ template const ElemLstType& GetElemList(const std::string& name = std::to_string(typeid(ElemLstType).hash_code())) const; /** * Delete an element-list. */ void DeleteElemList(const std::string& name); /** * Delete an element-list. */ template void DeleteElemList(); /** * Set target point coordinates. * * @param[in] Xtrg the coordinates of target points in array-of-struct * order: {x_1, y_1, z_1, x_2, ..., x_n, y_n, z_n} */ void SetTargetCoord(const Vector& Xtrg); /** * Set target point normals. * * @param[in] Xn_trg the coordinates of target points in array-of-struct * order: {nx_1, ny_1, nz_1, nx_2, ..., nx_n, ny_n, nz_n} */ void SetTargetNormal(const Vector& Xn_trg); /** * Get dimension of the boundary integral operator. Dim(0) is the input * dimension and Dim(1) is the output dimension. */ Long Dim(Integer k) const; /** * Setup the boundary integral operator. */ void Setup() const; /** * Clear setup data. */ void ClearSetup() const; /** * Evaluate the boundary integral operator. */ void ComputePotential(Vector& U, const Vector& F) const; private: void SetupBasic() const; void SetupFar() const; void SetupSelf() const; void SetupNear() const; void ComputeFarField(Vector& U, const Vector& F) const; void ComputeNearInterac(Vector& U, const Vector& F) const; struct ElemLstData { void (*SelfInterac)(Vector>&, const Kernel&, Real, bool, const ElementListBase*); void (*NearInterac)(Matrix&, const Vector&, const Vector&, const Kernel&, Real, const Long, const ElementListBase*); }; std::map*> elem_lst_map; std::map elem_data_map; Vector Xt; // User specified position of target points Vector Xnt; // User specified normal at target points Real tol_; Kernel ker_; bool trg_normal_dot_prod_; Comm comm_; mutable bool setup_flag; mutable Vector elem_lst_name; // name of each element-list (size=Nlst) mutable Vector elem_lst_cnt, elem_lst_dsp; // cnt and dsp of elements for each elem_lst (size=Nlst) mutable Vector elem_nds_cnt, elem_nds_dsp; // cnt and dsp of nodes for each element (size=Nelem) mutable Vector Xsurf; // Position of surface node points (target points for on-surface evaluation) mutable Vector Xn_surf; // Normal at surface node points (normal vector for on-surface evaluation) mutable Vector Xtrg; // Position of target points mutable Vector Xn_trg; // Normal vector at target points mutable bool setup_far_flag; mutable ParticleFMM fmm; mutable Vector elem_nds_cnt_far, elem_nds_dsp_far; // cnt and dsp of far-nodes for each element (size=Nelem) mutable Vector X_far, Xn_far, wts_far; // position, normal and weights for far-field quadrature mutable Vector dist_far; // minimum distance of target points for far-field evaluation mutable Vector F_far; // pre-allocated memory for density in far-field evaluation mutable bool setup_near_flag; mutable Vector Xtrg_near; // position of near-interaction target points sorted by element (size=Nnear*COORD_DIM) mutable Vector Xn_trg_near; // normal at near-interaction target points sorted by element (size=Nnear*COORD_DIM) mutable Vector near_scatter_index; // prmutation vector that takes near-interactions sorted by elem-idx to sorted by trg-idx (size=Nnear) mutable Vector near_trg_cnt, near_trg_dsp; // cnt and dsp of near-interactions for each target (size=Ntrg) mutable Vector near_elem_cnt, near_elem_dsp; // cnt and dsp of near-interaction for each element (size=Nelem) mutable Vector K_near_cnt, K_near_dsp; // cnt and dsp of element wise near-interaction matrix (size=Nelem) mutable Vector K_near; mutable bool setup_self_flag; mutable Vector> K_self; }; } #include SCTL_INCLUDE(boundary_integral.txx) #endif //_SCTL_BOUNDARY_INTEGRAL_HPP_