Skip to content

Commit

Permalink
Moved functionality to compute all the geometry aspects to move from …
Browse files Browse the repository at this point in the history 
…reference tetrahedron to physical space tetraherdon + compute their LSMs (Local Stiffness matrices) to the TetrahedronMeshManager class. Also, moved functionality to calculate electric field in cells to this class

Signed-off-by: ViNN280801 <vladislav_semykin01@mail.ru>
  • Loading branch information
@ViNN280801
ViNN280801 committed Dec 12, 2024
1 parent 2f6d5c7 commit 0aca952
Show file tree
Hide file tree
Showing 5 changed files with 214 additions and 140 deletions.
71 changes: 71 additions & 0 deletions include/DataHandling/TetrahedronMeshManager.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
#include <array>
#include <optional>

#include "FiniteElementMethod/FEMTypes.hpp"
#include "Geometry/GeometryTypes.hpp"
#include "Utilities/Utilities.hpp"

Expand Down Expand Up @@ -51,6 +52,10 @@ class TetrahedronMeshManager final
std::vector<std::vector<std::pair<size_t, std::array<size_t, 4>>>> m_elementsPerProc; ///< Elements partitioned by process (used for MPI distribution).
std::vector<std::vector<size_t>> 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<size_t, size_t> 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.
Expand All @@ -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;
Expand Down Expand Up @@ -259,6 +274,62 @@ class TetrahedronMeshManager final
*/
[[nodiscard("Tetrahedron centers are vital for geometric and physical calculations.")]]
std::map<size_t, Point> 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<DeviceType, Scalar, Scalar> *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
13 changes: 0 additions & 13 deletions include/FiniteElementMethod/GSMAssembler.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
139 changes: 136 additions & 3 deletions src/DataHandling/TetrahedronMeshManager.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
{
Expand Down Expand Up @@ -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);
Expand All @@ -245,6 +270,7 @@ TetrahedronMeshManager::TetrahedronMeshManager(std::string_view mesh_filename)
_receiveData();

_constructLocalMesh();
_buildViews();
}

void TetrahedronMeshManager::print() const noexcept
Expand Down Expand Up @@ -305,7 +331,7 @@ std::optional<TetrahedronMeshManager::TetrahedronData> 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;
Expand All @@ -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())
{
Expand Down Expand Up @@ -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())
{
Expand Down Expand Up @@ -400,3 +426,110 @@ std::map<size_t, Point> TetrahedronMeshManager::getTetrahedronCenters() const
}
return tetraCentres;
}

DynRankView TetrahedronMeshManager::computeLocalStiffnessMatricesAndNablaPhi(
Intrepid2::Basis<DeviceType, Scalar, Scalar> *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<DeviceType>::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<DeviceType>::setJacobianInv(invJacobians, jacobians);

DynRankView jacobiansDet("jacobiansDet", numTetrahedrons, numCubaturePoints);
Intrepid2::CellTools<DeviceType>::setJacobianDet(jacobiansDet, jacobians);

DynRankView cellMeasures("cellMeasures", numTetrahedrons, numCubaturePoints);
Intrepid2::FunctionSpaceTools<DeviceType>::computeCellMeasure(cellMeasures, jacobiansDet, cubWeights);

// 3. Transform reference grads
DynRankView transformedBasisGradients("transformedBasisGradients", numTetrahedrons, numBasisFunctions, numCubaturePoints, FEM_LIMITS_DEFAULT_SPACE_DIMENSION);
Intrepid2::FunctionSpaceTools<DeviceType>::HGRADtransformGRAD(transformedBasisGradients, invJacobians, referenceBasisGrads);

// 4. Weighted basis grads
DynRankView weightedBasisGrads("weightedBasisGrads", numTetrahedrons, numBasisFunctions, numCubaturePoints, FEM_LIMITS_DEFAULT_SPACE_DIMENSION);
Intrepid2::FunctionSpaceTools<DeviceType>::multiplyMeasure(weightedBasisGrads, cellMeasures, transformedBasisGradients);

// 5. Integrate local stiffness matrices
DynRankView localStiffnessMatrices("localStiffnessMatrices", numTetrahedrons, numBasisFunctions, numBasisFunctions);
Intrepid2::FunctionSpaceTools<DeviceType>::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());
}
}
Loading

0 comments on commit 0aca952

Please sign in to comment.