diff --git a/include/DataHandling/TetrahedronMeshManager.hpp b/include/DataHandling/TetrahedronMeshManager.hpp index f059be8..f476441 100644 --- a/include/DataHandling/TetrahedronMeshManager.hpp +++ b/include/DataHandling/TetrahedronMeshManager.hpp @@ -4,6 +4,7 @@ #include #include +#include "FiniteElementMethod/FEMTypes.hpp" #include "Geometry/GeometryTypes.hpp" #include "Utilities/Utilities.hpp" @@ -51,6 +52,10 @@ class TetrahedronMeshManager final std::vector>>> m_elementsPerProc; ///< Elements partitioned by process (used for MPI distribution). std::vector> m_nodeTagsPerProc; ///< Node tags partitioned by process (used for MPI distribution). + DynRankViewHost m_tetraVertices; ///< A view storing coordinates of all tetrahedron vertices. + DynRankViewHost m_tetraGlobalIds; ///< A view mapping local tetrahedron index to its global ID. + std::unordered_map m_globalToLocalTetraId; ///< Map from global tetra ID to local index. + /** * @brief Reads and partitions the mesh file. Only called by the root process (rank 0). * @param mesh_filename Path to the mesh file to be read and partitioned. @@ -72,6 +77,16 @@ class TetrahedronMeshManager final /// @brief Constructs the local mesh representation using the received data. void _constructLocalMesh(); + /** + * @brief Build DynRankViews and maps for easier index access. + * + * This method populates: + * - `m_tetraVertices`: A DynRankViewHost with shape [numTets,4,3], storing coordinates. + * - `m_tetraGlobalIds`: A DynRankViewHost mapping local indices to global tetra IDs. + * - `m_globalToLocalTetraId`: An unordered_map for global-to-local tetra ID lookup. + */ + void _buildViews(); + public: using NodeData = TetrahedronMeshManager::TetrahedronData::NodeData; using TetrahedronData = TetrahedronMeshManager::TetrahedronData; @@ -259,6 +274,62 @@ class TetrahedronMeshManager final */ [[nodiscard("Tetrahedron centers are vital for geometric and physical calculations.")]] std::map getTetrahedronCenters() const; + + /** + * @brief Get the local tetrahedron ID given a global tetrahedron ID. + * + * @param globalId The global ID of the tetrahedron (from Gmsh). + * @return The local index (0-based) of the tetrahedron. + * @throw std::out_of_range if globalId is not found. + */ + size_t getLocalTetraId(size_t globalId) const { return m_globalToLocalTetraId.at(globalId); } + + /** + * @brief Access the DynRankView of tetrahedron vertices. + * + * @return A copy of the DynRankViewHost that stores tetrahedron vertices. + */ + inline DynRankViewHost getTetraVerticesView() const noexcept { return m_tetraVertices; } + + /** + * @brief Access the DynRankView of global tetrahedron IDs. + * + * @return A copy of the DynRankViewHost that stores global IDs for each local tetrahedron index. + */ + inline DynRankViewHost getTetraGlobalIdsView() const noexcept { return m_tetraGlobalIds; } + + /** + * @brief Compute local stiffness matrices and assign nablaPhi for each node of each tetrahedron. + * + * This method internally computes jacobians, inverse jacobians, cell measures, + * transforms the reference gradients of the basis functions, and integrates + * to form local stiffness matrices. It also assigns nablaPhi to each node. + * + * @param basis The chosen basis functions object (already constructed outside). + * @param cubPoints The cubature points. + * @param cubWeights The cubature weights. + * @param numBasisFunctions The number of basis functions per tetrahedron. + * @param numCubaturePoints The number of cubature points. + * @return DynRankView with shape (numTetrahedrons, numBasisFunctions, numBasisFunctions) + * containing the local stiffness matrices. + */ + DynRankView computeLocalStiffnessMatricesAndNablaPhi( + Intrepid2::Basis *basis, + DynRankView const &cubPoints, + DynRankView const &cubWeights, + size_t numBasisFunctions, + size_t numCubaturePoints); + + /** + * @brief Compute the electric field for each tetrahedron based on assigned potentials and nablaPhi. + * + * This method loops over all tetrahedrons, computes: + * \f$ E_{\text{cell}} = \sum_i (\varphi_i * \nabla \varphi_i) \f$ + * for each tetrahedron, and updates the tetrahedron's electricField field. + * + * If either potential or nablaPhi is not set for a node, a warning is issued and that node is skipped in the sum. + */ + void computeElectricFields(); }; #endif // !TETRAHEDRONMESHMANAGER_HPP diff --git a/include/FiniteElementMethod/GSMAssembler.hpp b/include/FiniteElementMethod/GSMAssembler.hpp index 10bb4ec..24aaf0f 100644 --- a/include/FiniteElementMethod/GSMAssembler.hpp +++ b/include/FiniteElementMethod/GSMAssembler.hpp @@ -45,19 +45,6 @@ class GSMAssembler CubatureManager m_cubature_manager; ///< Cubature manager (manages cubature points and weights based on cell type). MatrixManager m_matrix_manager; ///< Matrix manager (manages initialization and filling of matrices). - /** - * @brief Retrieves the vertices of all tetrahedrons in the mesh. - * - * This function extracts the coordinates of the vertices for each tetrahedron - * in the mesh and stores them in a multi-dimensional array (Dynamic Rank View). - * The dimensions of the array are [number of tetrahedrons] x [4 vertices] x [3 coordinates (x, y, z)]. - * - * @return A multi-dimensional array containing the vertices of all tetrahedrons. - * Each tetrahedron is represented by its four vertices, and each vertex has three coordinates (x, y, z). - * @throw std::runtime_error if an error occurs during the extraction of vertices. - */ - DynRankViewHost _getTetrahedronVertices(); - /** * @brief Retrieves matrix entries from calculated local stiffness matrices. * @return A vector of matrix entries, each containing global row, column, and value. diff --git a/src/DataHandling/TetrahedronMeshManager.cpp b/src/DataHandling/TetrahedronMeshManager.cpp index e3577f7..5e55664 100644 --- a/src/DataHandling/TetrahedronMeshManager.cpp +++ b/src/DataHandling/TetrahedronMeshManager.cpp @@ -5,6 +5,9 @@ #endif #include "DataHandling/TetrahedronMeshManager.hpp" +#include "FiniteElementMethod/Cell/CellSelector.hpp" +#include "FiniteElementMethod/FEMLimits.hpp" +#include "Geometry/MathVector.hpp" Point TetrahedronMeshManager::TetrahedronData::getTetrahedronCenter() const { @@ -225,6 +228,28 @@ void TetrahedronMeshManager::_constructLocalMesh() } } +void TetrahedronMeshManager::_buildViews() +{ + m_tetraVertices = DynRankViewHost("tetraVertices", getNumTetrahedrons(), 4, 3); + m_tetraGlobalIds = DynRankViewHost("tetraGlobalIds", getNumTetrahedrons()); + + for (size_t localId{}; localId < getNumTetrahedrons(); ++localId) + { + auto const &data{m_meshComponents[localId]}; + m_tetraGlobalIds(localId) = data.globalTetraId; + for (short node{}; node < 4; ++node) + { + m_tetraVertices(localId, node, 0) = data.nodes[node].nodeCoords.x(); + m_tetraVertices(localId, node, 1) = data.nodes[node].nodeCoords.y(); + m_tetraVertices(localId, node, 2) = data.nodes[node].nodeCoords.z(); + } + } + + m_globalToLocalTetraId.reserve(getNumTetrahedrons()); + for (size_t localId{}; localId < getNumTetrahedrons(); ++localId) + m_globalToLocalTetraId[m_tetraGlobalIds(localId)] = localId; +} + TetrahedronMeshManager::TetrahedronMeshManager(std::string_view mesh_filename) { util::check_gmsh_mesh_file(mesh_filename); @@ -245,6 +270,7 @@ TetrahedronMeshManager::TetrahedronMeshManager(std::string_view mesh_filename) _receiveData(); _constructLocalMesh(); + _buildViews(); } void TetrahedronMeshManager::print() const noexcept @@ -305,7 +331,7 @@ std::optional TetrahedronMeshManager::g { return data.globalTetraId == globalTetrahedronId; })}; #else auto it{std::find_if(m_meshComponents.cbegin(), m_meshComponents.cend(), [globalTetrahedronId](const TetrahedronData &data) - { return data.globalTetraId == globalTetrahedronId; })}; + { return data.globalTetraId == globalTetrahedronId; })}; #endif if (it != m_meshComponents.cend()) return *it; @@ -319,7 +345,7 @@ void TetrahedronMeshManager::assignNablaPhi(size_t tetrahedronId, size_t nodeId, { return data.globalTetraId == tetrahedronId; })}; #else auto it{std::find_if(m_meshComponents.begin(), m_meshComponents.end(), [tetrahedronId](const TetrahedronData &data) - { return data.globalTetraId == tetrahedronId; })}; + { return data.globalTetraId == tetrahedronId; })}; #endif if (it != m_meshComponents.cend()) { @@ -349,7 +375,7 @@ void TetrahedronMeshManager::assignElectricField(size_t tetrahedronId, Point con { return data.globalTetraId == tetrahedronId; })}; #else auto it{std::find_if(m_meshComponents.begin(), m_meshComponents.end(), [tetrahedronId](const TetrahedronData &data) - { return data.globalTetraId == tetrahedronId; })}; + { return data.globalTetraId == tetrahedronId; })}; #endif if (it != m_meshComponents.end()) { @@ -400,3 +426,110 @@ std::map TetrahedronMeshManager::getTetrahedronCenters() const } return tetraCentres; } + +DynRankView TetrahedronMeshManager::computeLocalStiffnessMatricesAndNablaPhi( + Intrepid2::Basis *basis, + DynRankView const &cubPoints, + DynRankView const &cubWeights, + size_t numBasisFunctions, + size_t numCubaturePoints) +{ + auto const numTetrahedrons{getNumTetrahedrons()}; + + // 1. Reference basis grads + DynRankView referenceBasisGrads("referenceBasisGrads", numBasisFunctions, numCubaturePoints, FEM_LIMITS_DEFAULT_SPACE_DIMENSION); + basis->getValues(referenceBasisGrads, cubPoints, Intrepid2::OPERATOR_GRAD); + +#ifdef USE_CUDA + auto vertices{Kokkos::create_mirror_view_and_copy(MemorySpaceDevice(), m_tetraVertices)}; // From host to device. +#else + DynRankView vertices(m_tetraVertices); +#endif + + // 2. Jacobians and measures + DynRankView jacobians("jacobians", numTetrahedrons, numCubaturePoints, FEM_LIMITS_DEFAULT_SPACE_DIMENSION, FEM_LIMITS_DEFAULT_SPACE_DIMENSION); + Intrepid2::CellTools::setJacobian(jacobians, cubPoints, vertices, CellSelector::get(CellType::Tetrahedron)); + + DynRankView invJacobians("invJacobians", numTetrahedrons, numCubaturePoints, FEM_LIMITS_DEFAULT_SPACE_DIMENSION, FEM_LIMITS_DEFAULT_SPACE_DIMENSION); + Intrepid2::CellTools::setJacobianInv(invJacobians, jacobians); + + DynRankView jacobiansDet("jacobiansDet", numTetrahedrons, numCubaturePoints); + Intrepid2::CellTools::setJacobianDet(jacobiansDet, jacobians); + + DynRankView cellMeasures("cellMeasures", numTetrahedrons, numCubaturePoints); + Intrepid2::FunctionSpaceTools::computeCellMeasure(cellMeasures, jacobiansDet, cubWeights); + + // 3. Transform reference grads + DynRankView transformedBasisGradients("transformedBasisGradients", numTetrahedrons, numBasisFunctions, numCubaturePoints, FEM_LIMITS_DEFAULT_SPACE_DIMENSION); + Intrepid2::FunctionSpaceTools::HGRADtransformGRAD(transformedBasisGradients, invJacobians, referenceBasisGrads); + + // 4. Weighted basis grads + DynRankView weightedBasisGrads("weightedBasisGrads", numTetrahedrons, numBasisFunctions, numCubaturePoints, FEM_LIMITS_DEFAULT_SPACE_DIMENSION); + Intrepid2::FunctionSpaceTools::multiplyMeasure(weightedBasisGrads, cellMeasures, transformedBasisGradients); + + // 5. Integrate local stiffness matrices + DynRankView localStiffnessMatrices("localStiffnessMatrices", numTetrahedrons, numBasisFunctions, numBasisFunctions); + Intrepid2::FunctionSpaceTools::integrate(localStiffnessMatrices, weightedBasisGrads, transformedBasisGradients); + +#ifdef USE_CUDA + auto localStiffnessMatrices_host{Kokkos::create_mirror_view_and_copy(MemorySpaceHost(), localStiffnessMatrices)}; + auto weightedBasisGrads_host{Kokkos::create_mirror_view_and_copy(MemorySpaceHost(), weightedBasisGrads)}; +#else + auto &weightedBasisGrads_host{weightedBasisGrads}; +#endif + + // 6. Assign nablaPhi + for (size_t localTetraId{}; localTetraId < numTetrahedrons; ++localTetraId) + { + auto globalTetraId{getGlobalTetraId(localTetraId)}; + auto const &nodes{getTetrahedronNodes(localTetraId)}; + + if (numBasisFunctions > nodes.size()) + { + WARNINGMSG("Basis function count exceeds the number of nodes."); + } + + for (short localNodeId{}; localNodeId < (short)numBasisFunctions; ++localNodeId) + { + auto globalNodeId{nodes[localNodeId].globalNodeId}; + Point grad(weightedBasisGrads_host(localTetraId, localNodeId, 0, 0), + weightedBasisGrads_host(localTetraId, localNodeId, 0, 1), + weightedBasisGrads_host(localTetraId, localNodeId, 0, 2)); + + assignNablaPhi(globalTetraId, globalNodeId, grad); + } + } + + return localStiffnessMatrices; +} + +void TetrahedronMeshManager::computeElectricFields() +{ + // We have a mapping of (Tetrahedron ID -> (Node ID -> Gradient of basis function)). + // To obtain the electric field of the cell, we need to sum over all nodes the product + // of the node's potential (φ_i) and the node's basis function gradient (∇φ_i). + // + // The formula is: + // E_cell = Σ(φ_i * ∇φ_i) + // where i is the global index of the node. + for (auto &tetra : m_meshComponents) + { + ElectricField electricField{}; + bool missingData{}; + for (auto const &node : tetra.nodes) + { + if (node.potential && node.nablaPhi) + electricField += ElectricField(node.nablaPhi->x(), node.nablaPhi->y(), node.nablaPhi->z()) * node.potential.value(); + else + missingData = true; + } + + if (missingData) + { + WARNINGMSG(util::stringify("Warning: Node potential or nablaPhi is not set for one or more nodes of tetrahedron ", + tetra.globalTetraId, ". Those nodes are skipped in electric field computation.\n")); + } + + tetra.electricField = Point(electricField.getX(), electricField.getY(), electricField.getZ()); + } +} diff --git a/src/FiniteElementMethod/GSMAssembler.cpp b/src/FiniteElementMethod/GSMAssembler.cpp index c4e12f1..0214625 100644 --- a/src/FiniteElementMethod/GSMAssembler.cpp +++ b/src/FiniteElementMethod/GSMAssembler.cpp @@ -5,37 +5,6 @@ #include "FiniteElementMethod/FEMCheckers.hpp" #include "FiniteElementMethod/FEMLimits.hpp" -DynRankViewHost GSMAssembler::_getTetrahedronVertices() -{ - try - { - DynRankViewHost vertices("vertices", m_meshManager.getNumTetrahedrons(), FEM_LIMITS_DEFAULT_TETRAHEDRON_VERTICES_COUNT, FEM_LIMITS_DEFAULT_SPACE_DIMENSION); - size_t i{}; - for (auto const &meshParam : m_meshManager.getMeshComponents()) - { - auto tetrahedron{meshParam.tetrahedron}; - for (short node{}; node < FEM_LIMITS_DEFAULT_TETRAHEDRON_VERTICES_COUNT; ++node) - { - vertices(i, node, 0) = tetrahedron.vertex(node).x(); - vertices(i, node, 1) = tetrahedron.vertex(node).y(); - vertices(i, node, 2) = tetrahedron.vertex(node).z(); - } - ++i; - } - return vertices; - } - catch (std::exception const &ex) - { - ERRMSG(ex.what()); - } - catch (...) - { - ERRMSG("Unknown error"); - } - WARNINGMSG("Returning empty multidimensional array with vertices of the all tetrahedrons from the mesh"); - return DynRankViewHost(); -} - std::vector GSMAssembler::_getMatrixEntries() { TetrahedronIndicesVector globalNodeIndicesPerElement; @@ -100,84 +69,19 @@ DynRankView GSMAssembler::_computeLocalStiffnessMatrices() { try { - auto const numTetrahedrons{m_meshManager.getNumTetrahedrons()}; auto const numBasisFunctions{m_cubature_manager.getCountBasisFunctions()}; auto const numCubaturePoints{m_cubature_manager.getCountCubaturePoints()}; auto const &cubPoints{m_cubature_manager.getCubaturePoints()}; auto const &cubWeights{m_cubature_manager.getCubatureWeights()}; - // 1. Calculating basis gradients. - DynRankView referenceBasisGrads("referenceBasisGrads", numBasisFunctions, numCubaturePoints, FEM_LIMITS_DEFAULT_SPACE_DIMENSION); auto basis{BasisSelector::template get(CellType::Tetrahedron, m_polynom_order)}; - basis->getValues(referenceBasisGrads, cubPoints, Intrepid2::OPERATOR_GRAD); - -#ifdef USE_CUDA - auto vertices{Kokkos::create_mirror_view_and_copy(MemorySpaceDevice(), _getTetrahedronVertices())}; // From host to device. -#else - DynRankView vertices(_getTetrahedronVertices()); -#endif - - // 2. Computing cell jacobians, inversed jacobians and jacobian determinants to get cell measure. - DynRankView jacobians("jacobians", numTetrahedrons, numCubaturePoints, FEM_LIMITS_DEFAULT_SPACE_DIMENSION, FEM_LIMITS_DEFAULT_SPACE_DIMENSION); - Intrepid2::CellTools::setJacobian(jacobians, cubPoints, vertices, CellSelector::get(CellType::Tetrahedron)); - - DynRankView invJacobians("invJacobians", numTetrahedrons, numCubaturePoints, FEM_LIMITS_DEFAULT_SPACE_DIMENSION, FEM_LIMITS_DEFAULT_SPACE_DIMENSION); - Intrepid2::CellTools::setJacobianInv(invJacobians, jacobians); - - DynRankView jacobiansDet("jacobiansDet", numTetrahedrons, numCubaturePoints); - Intrepid2::CellTools::setJacobianDet(jacobiansDet, jacobians); - - DynRankView cellMeasures("cellMeasures", numTetrahedrons, numCubaturePoints); - Intrepid2::FunctionSpaceTools::computeCellMeasure(cellMeasures, jacobiansDet, cubWeights); - - // 3. Transforming reference basis gradients to physical frame. - DynRankView transformedBasisGradients("transformedBasisGradients", numTetrahedrons, numBasisFunctions, numCubaturePoints, FEM_LIMITS_DEFAULT_SPACE_DIMENSION); - Intrepid2::FunctionSpaceTools::HGRADtransformGRAD(transformedBasisGradients, invJacobians, referenceBasisGrads); - - // 4. Multiply transformed basis gradients by cell measures to get weighted gradients. - DynRankView weightedBasisGrads("weightedBasisGrads", numTetrahedrons, numBasisFunctions, numCubaturePoints, FEM_LIMITS_DEFAULT_SPACE_DIMENSION); - Intrepid2::FunctionSpaceTools::multiplyMeasure(weightedBasisGrads, cellMeasures, transformedBasisGradients); - - // 5. Integrate to get local stiffness matrices for workset cells. - DynRankView localStiffnessMatrices("localStiffnessMatrices", numTetrahedrons, numBasisFunctions, numBasisFunctions); - Intrepid2::FunctionSpaceTools::integrate(localStiffnessMatrices, weightedBasisGrads, transformedBasisGradients); - -#ifdef USE_CUDA - auto localStiffnessMatrices_host{Kokkos::create_mirror_view_and_copy(MemorySpaceHost(), localStiffnessMatrices)}; // From device to host. - auto weightedBasisGrads_host{Kokkos::create_mirror_view_and_copy(MemorySpaceHost(), weightedBasisGrads)}; // From device to host. -#endif - - // 6. Filling map with basis grads on each node on each tetrahedron. - for (size_t localTetraId{}; localTetraId < numTetrahedrons; ++localTetraId) - { - auto globalTetraId{m_meshManager.getGlobalTetraId(localTetraId)}; - auto const &nodes{m_meshManager.getTetrahedronNodes(localTetraId)}; - if (numBasisFunctions > nodes.size()) - { - ERRMSG("Basis function count exceeds the number of nodes."); - } - - for (short localNodeId{}; localNodeId < numBasisFunctions; ++localNodeId) - { - auto globalNodeId{m_meshManager.getGlobalNodeId(localTetraId, localNodeId)}; - -#ifdef USE_CUDA - Point grad(weightedBasisGrads_host(localTetraId, localNodeId, 0, 0), - weightedBasisGrads_host(localTetraId, localNodeId, 0, 1), - weightedBasisGrads_host(localTetraId, localNodeId, 0, 2)); -#else - Point grad(weightedBasisGrads(localTetraId, localNodeId, 0, 0), - weightedBasisGrads(localTetraId, localNodeId, 0, 1), - weightedBasisGrads(localTetraId, localNodeId, 0, 2)); -#endif - - /// As we have polynomial order \( = 1 \), all the values from \( \nabla \varphi \) - /// in all cubature points are the same, so we can add only one row from each \( \nabla \varphi \). - m_meshManager.assignNablaPhi(globalTetraId, globalNodeId, grad); - } - } - return localStiffnessMatrices; + return m_meshManager.computeLocalStiffnessMatricesAndNablaPhi( + basis.get(), + cubPoints, + cubWeights, + numBasisFunctions, + numCubaturePoints); } catch (std::exception const &ex) { diff --git a/src/FiniteElementMethod/MatrixEquationSolver.cpp b/src/FiniteElementMethod/MatrixEquationSolver.cpp index bb4dc95..d732c2d 100644 --- a/src/FiniteElementMethod/MatrixEquationSolver.cpp +++ b/src/FiniteElementMethod/MatrixEquationSolver.cpp @@ -53,29 +53,8 @@ void MatrixEquationSolver::calculateElectricField() { try { - // Filling node potentials before calculating the electric field in the cell. fillNodesPotential(); - - // We have map: (Tetrahedron ID | map). - // To get electric field of the cell we just need to accumulate all the basis func grads for each node for each tetrahedron: - // \f$ E_{\text{cell}} = \sum (\varphi_i \cdot \nabla \varphi_i) \f$, - // where \( i \) is the global index of the node. - for (auto const &tetrahedronData : m_assembler->getMeshManager().getMeshComponents()) - { - ElectricField electricField{}; - for (auto const &node : tetrahedronData.nodes) - { - if (node.potential && node.nablaPhi) - electricField += ElectricField(node.nablaPhi.value().x(), node.nablaPhi.value().y(), node.nablaPhi.value().z()) * - node.potential.value(); - else - { - WARNINGMSG(util::stringify("Node potential or nablaPhi is not set for the ", - node.globalNodeId, " vertex of the ", tetrahedronData.globalTetraId, " tetrahedron")); - } - } - m_assembler->getMeshManager().assignElectricField(tetrahedronData.globalTetraId, Point(electricField.getX(), electricField.getY(), electricField.getZ())); - } + m_assembler->getMeshManager().computeElectricFields(); } catch (std::exception const &ex) {