Skip to content

Commit

Permalink
Merge branch 'RefactorSpatialPosition' into 'master'
Browse files Browse the repository at this point in the history 
Refactor SpatialPosition for more efficient memory usage

See merge request ogs/ogs!5212
  • Loading branch information
Dmitry Yu. Naumov committed Jan 30, 2025
2 parents f61a95e + 0ad7569 commit ffa4777
Show file tree
Hide file tree
Showing 48 changed files with 85 additions and 196 deletions.
5 changes: 2 additions & 3 deletions ParameterLib/Parameter.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -170,12 +170,11 @@ struct Parameter : public ParameterBase
n_components);

// Column vector of values, copied for each node.
SpatialPosition x_position;
auto const nodes = element.getNodes();
for (int i = 0; i < n_nodes; ++i)
{
x_position.setAll(
nodes[i]->getID(), element.getID(), std::nullopt, *nodes[i]);
SpatialPosition const x_position{nodes[i]->getID(), element.getID(),
*nodes[i]};
auto const& values = this->operator()(t, x_position);
auto const row_values =
Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, 1> const>(
Expand Down
84 changes: 41 additions & 43 deletions ParameterLib/SpatialPosition.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@

#pragma once

#include <bitset>
#include <optional>

#include "MathLib/Point3d.h"
Expand All @@ -21,7 +22,7 @@ namespace ParameterLib
//!
//! The setters of this class make sure that only compatible information can be
//! stored at the same time; e.g., it is not possible to specify an element ID
//! and a node ID at the same time (the setAll() method being an exception to
//! and a node ID at the same time (the constructor being an exception to
//! that rule).
class SpatialPosition
{
Expand All @@ -30,73 +31,70 @@ class SpatialPosition

SpatialPosition(std::optional<std::size_t> const& node_id,
std::optional<std::size_t> const& element_id,
std::optional<unsigned> const& integration_point,
std::optional<MathLib::Point3d> const& coordinates)
: _node_id(node_id),
_element_id(element_id),
_integration_point(integration_point),
_coordinates(coordinates)
{
if (node_id)
{
_node_id = *node_id;
flags.set(node_bit);
}

if (element_id)
{
_element_id = *element_id;
flags.set(element_bit);
}
if (coordinates)
{
_coordinates = *coordinates;
flags.set(coordinates_bit);
}
}

std::optional<std::size_t> getNodeID() const { return _node_id; }
std::optional<std::size_t> getElementID() const { return _element_id; }
std::optional<unsigned> getIntegrationPoint() const
std::optional<std::size_t> getNodeID() const
{
return flags[node_bit] ? std::make_optional(_node_id) : std::nullopt;
}
std::optional<std::size_t> getElementID() const
{
return _integration_point;
return flags[element_bit] ? std::make_optional(_element_id)
: std::nullopt;
}
std::optional<MathLib::Point3d> const& getCoordinates() const
std::optional<MathLib::Point3d> const getCoordinates() const
{
return _coordinates;
return flags[coordinates_bit] ? std::make_optional(_coordinates)
: std::nullopt;
}

void setNodeID(std::size_t node_id)
{
clear();
flags.reset();
flags.set(node_bit);
_node_id = node_id;
}

void setElementID(std::size_t element_id)
{
clear();
flags.reset();
flags.set(element_bit);
_element_id = element_id;
}

void setIntegrationPoint(unsigned integration_point)
{
assert(_element_id);
_integration_point = integration_point;
}

void setCoordinates(MathLib::Point3d const& coordinates)
{
_coordinates = coordinates;
}

void setAll(std::optional<std::size_t> const& node_id,
std::optional<std::size_t> const& element_id,
std::optional<unsigned> const& integration_point,
std::optional<MathLib::Point3d> const& coordinates)
{
_node_id = node_id;
_element_id = element_id;
_integration_point = integration_point;
_coordinates = coordinates;
}

void clear()
{
_node_id = std::nullopt;
_element_id = std::nullopt;
_integration_point = std::nullopt;
_coordinates = std::nullopt;
flags.set(coordinates_bit);
}

private:
std::optional<std::size_t> _node_id;
std::optional<std::size_t> _element_id;
std::optional<unsigned> _integration_point;
std::optional<MathLib::Point3d> _coordinates;
std::size_t _node_id = 0;
std::size_t _element_id = 0;
MathLib::Point3d _coordinates{};

std::bitset<3> flags{};
static constexpr std::size_t node_bit = 0;
static constexpr std::size_t element_bit = 1;
static constexpr std::size_t coordinates_bit = 2;
};

} // namespace ParameterLib
8 changes: 4 additions & 4 deletions ProcessLib/BoundaryConditionAndSourceTerm/ConstraintDirichletBoundaryCondition.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -128,8 +128,6 @@ void ConstraintDirichletBoundaryCondition::getEssentialBCValues(
const double t, const GlobalVector& /*x*/,
NumLib::IndexValueVector<GlobalIndexType>& bc_values) const
{
ParameterLib::SpatialPosition pos;

bc_values.ids.clear();
bc_values.values.clear();

Expand All @@ -153,8 +151,10 @@ void ConstraintDirichletBoundaryCondition::getEssentialBCValues(
unsigned const number_nodes = boundary_element->getNumberOfNodes();
for (unsigned i = 0; i < number_nodes; ++i)
{
auto const id = boundary_element->getNode(i)->getID();
pos.setAll(id, boundary_element->getID(), {}, {});
auto const& node = *boundary_element->getNode(i);
auto const id = node.getID();
ParameterLib::SpatialPosition const pos{
id, boundary_element->getID(), node};

MeshLib::Location l(_bc_mesh.getID(), MeshLib::MeshItemType::Node,
id);
Expand Down
2 changes: 1 addition & 1 deletion ProcessLib/BoundaryConditionAndSourceTerm/NeumannBoundaryConditionLocalAssembler.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,7 @@ class NeumannBoundaryConditionLocalAssembler final
auto const& w = ip_data.weight;

ParameterLib::SpatialPosition const position{
std::nullopt, Base::_element.getID(), ip,
std::nullopt, Base::_element.getID(),
MathLib::Point3d(
NumLib::interpolateCoordinates<ShapeFunction,
ShapeMatricesType>(
Expand Down
2 changes: 1 addition & 1 deletion ProcessLib/BoundaryConditionAndSourceTerm/RobinBoundaryConditionLocalAssembler.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,7 @@ class RobinBoundaryConditionLocalAssembler final
auto const& w = ip_data.weight;

ParameterLib::SpatialPosition const position{
std::nullopt, Base::_element.getID(), ip,
std::nullopt, Base::_element.getID(),
MathLib::Point3d(
NumLib::interpolateCoordinates<ShapeFunction,
ShapeMatricesType>(
Expand Down
2 changes: 1 addition & 1 deletion ProcessLib/BoundaryConditionAndSourceTerm/VolumetricSourceTermFEM.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,7 @@ class VolumetricSourceTermLocalAssembler final
auto const& w = _ip_data[ip].integration_weight;

ParameterLib::SpatialPosition const pos{
std::nullopt, _element.getID(), ip,
std::nullopt, _element.getID(),
MathLib::Point3d(
NumLib::interpolateCoordinates<ShapeFunction,
ShapeMatricesType>(_element,
Expand Down
20 changes: 0 additions & 20 deletions ProcessLib/ComponentTransport/ComponentTransportFEM.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -281,8 +281,6 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface
shape_matrices[ip].integralMeasure *
shape_matrices[ip].detJ * aperture_size);

pos.setIntegrationPoint(ip);

_ip_data[ip].porosity =
medium[MaterialPropertyLib::PropertyType::porosity]
.template initialValue<double>(
Expand Down Expand Up @@ -606,8 +604,6 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface

for (unsigned ip(0); ip < n_integration_points; ++ip)
{
pos.setIntegrationPoint(ip);

auto& ip_data = _ip_data[ip];
auto const& dNdx = ip_data.dNdx;
auto const& N = Ns[ip];
Expand Down Expand Up @@ -867,8 +863,6 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface

for (unsigned ip(0); ip < n_integration_points; ++ip)
{
pos.setIntegrationPoint(ip);

auto& ip_data = _ip_data[ip];
auto const& dNdx = ip_data.dNdx;
auto const& w = ip_data.integration_weight;
Expand Down Expand Up @@ -995,8 +989,6 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface

for (unsigned ip(0); ip < n_integration_points; ip++)
{
pos.setIntegrationPoint(ip);

auto const& ip_data = this->_ip_data[ip];
auto const& dNdx = ip_data.dNdx;
auto const& w = ip_data.integration_weight;
Expand Down Expand Up @@ -1142,8 +1134,6 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface

for (unsigned ip(0); ip < n_integration_points; ++ip)
{
pos.setIntegrationPoint(ip);

auto& ip_data = _ip_data[ip];
auto const& dNdx = ip_data.dNdx;
auto const& w = ip_data.integration_weight;
Expand Down Expand Up @@ -1342,8 +1332,6 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface

for (unsigned ip(0); ip < n_integration_points; ++ip)
{
pos.setIntegrationPoint(ip);

auto& ip_data = _ip_data[ip];
auto const& dNdx = ip_data.dNdx;
auto const& w = ip_data.integration_weight;
Expand Down Expand Up @@ -1468,8 +1456,6 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface

for (unsigned ip(0); ip < n_integration_points; ++ip)
{
pos.setIntegrationPoint(ip);

auto& ip_data = _ip_data[ip];
auto const& dNdx = ip_data.dNdx;
auto const& w = ip_data.integration_weight;
Expand Down Expand Up @@ -1597,8 +1583,6 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface

for (unsigned ip(0); ip < n_integration_points; ++ip)
{
pos.setIntegrationPoint(ip);

auto& ip_data = _ip_data[ip];
auto const w = ip_data.integration_weight;
auto const& N = Ns[ip];
Expand Down Expand Up @@ -1723,8 +1707,6 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface
auto const& N = Ns[ip];
auto const& porosity = ip_data.porosity;

pos.setIntegrationPoint(ip);

double C_int_pt = 0.0;
double p_int_pt = 0.0;

Expand Down Expand Up @@ -1953,8 +1935,6 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface
auto const& N = Ns[ip];
auto const& phi = ip_data.porosity;

pos.setIntegrationPoint(ip);

double const p_ip = N.dot(p);
double const c_ip = N.dot(c);

Expand Down
2 changes: 0 additions & 2 deletions ProcessLib/HT/HTFEM.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -277,8 +277,6 @@ class HTFEM : public HTLocalAssemblerInterface
auto const& dNdx = ip_data.dNdx;
auto const& N = Ns[ip];

pos.setIntegrationPoint(ip);

double T_int_pt = 0.0;
double p_int_pt = 0.0;
NumLib::shapeFunctionInterpolate(local_p, N, p_int_pt);
Expand Down
2 changes: 0 additions & 2 deletions ProcessLib/HT/MonolithicHTFEM.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,8 +127,6 @@ class MonolithicHTFEM : public HTFEM<ShapeFunction, GlobalDim>

for (unsigned ip(0); ip < n_integration_points; ip++)
{
pos.setIntegrationPoint(ip);

auto const& ip_data = this->_ip_data[ip];
auto const& dNdx = ip_data.dNdx;
auto const& N = Ns[ip];
Expand Down
4 changes: 0 additions & 4 deletions ProcessLib/HT/StaggeredHTFEM-impl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,8 +85,6 @@ void StaggeredHTFEM<ShapeFunction, GlobalDim>::assembleHydraulicEquation(

for (unsigned ip(0); ip < n_integration_points; ip++)
{
pos.setIntegrationPoint(ip);

auto const& ip_data = this->_ip_data[ip];
auto const& dNdx = ip_data.dNdx;
auto const& N = Ns[ip];
Expand Down Expand Up @@ -218,8 +216,6 @@ void StaggeredHTFEM<ShapeFunction, GlobalDim>::assembleHeatTransportEquation(

for (unsigned ip(0); ip < n_integration_points; ip++)
{
pos.setIntegrationPoint(ip);

auto const& ip_data = this->_ip_data[ip];
auto const& dNdx = ip_data.dNdx;
auto const& N = Ns[ip];
Expand Down
4 changes: 1 addition & 3 deletions ProcessLib/HeatConduction/HeatConductionFEM.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,7 +109,7 @@ class LocalAssemblerData : public HeatConductionLocalAssemblerInterface
{
auto const& sm = _shape_matrices[ip];
ParameterLib::SpatialPosition const pos{
std::nullopt, _element.getID(), ip,
std::nullopt, _element.getID(),
MathLib::Point3d(
NumLib::interpolateCoordinates<ShapeFunction,
ShapeMatricesType>(_element,
Expand Down Expand Up @@ -188,7 +188,6 @@ class LocalAssemblerData : public HeatConductionLocalAssemblerInterface

for (unsigned ip = 0; ip < n_integration_points; ip++)
{
pos.setIntegrationPoint(ip);
auto const& sm = _shape_matrices[ip];
double const w =
_integration_method.getWeightedPoint(ip).getWeight() * sm.detJ *
Expand Down Expand Up @@ -269,7 +268,6 @@ class LocalAssemblerData : public HeatConductionLocalAssemblerInterface

for (unsigned ip = 0; ip < n_integration_points; ip++)
{
pos.setIntegrationPoint(ip);
auto const& sm = _shape_matrices[ip];
// get the local temperature and put it in the variable array for
// access in MPL
Expand Down
3 changes: 0 additions & 3 deletions ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil-impl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,8 +55,6 @@ HeatTransportBHELocalAssemblerSoil<ShapeFunction>::
// ip data initialization
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
x_position.setIntegrationPoint(ip);

// create the class IntegrationPointDataBHE in place
auto const& sm = _shape_matrices[ip];
double const w = _integration_method.getWeightedPoint(ip).getWeight() *
Expand Down Expand Up @@ -96,7 +94,6 @@ void HeatTransportBHELocalAssemblerSoil<ShapeFunction>::assemble(

for (unsigned ip = 0; ip < n_integration_points; ip++)
{
pos.setIntegrationPoint(ip);
auto& ip_data = _ip_data[ip];
auto const& N = ip_data.N;
auto const& dNdx = ip_data.dNdx;
Expand Down
Loading

0 comments on commit ffa4777

Please sign in to comment.