Skip to content

Commit

Permalink
[GeoMechanicsApplication] Small refactoring of T_normal_flux_condition (
Browse files Browse the repository at this point in the history 
#11871)

* Renamed fluxes to reflect their purpose. Factored out function and struct.

* Naming of local variables.

* Review comment of Gennady taken up.
  • Loading branch information
@WPK4FEM
WPK4FEM authored Dec 4, 2023
1 parent 16cf50b commit e434e2c
Show file tree
Hide file tree
Showing 2 changed files with 23 additions and 51 deletions.
48 changes: 17 additions & 31 deletions applications/GeoMechanicsApplication/custom_conditions/T_normal_flux_condition.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,15 +26,15 @@ GeoTNormalFluxCondition<TDim, TNumNodes>::GeoTNormalFluxCondition()

template <unsigned int TDim, unsigned int TNumNodes>
GeoTNormalFluxCondition<TDim, TNumNodes>::GeoTNormalFluxCondition(IndexType NewId,
GeometryType::Pointer pGeometry)
GeometryType::Pointer pGeometry)
: GeoTCondition<TDim, TNumNodes>(NewId, pGeometry)
{
}

template <unsigned int TDim, unsigned int TNumNodes>
GeoTNormalFluxCondition<TDim, TNumNodes>::GeoTNormalFluxCondition(IndexType NewId,
GeometryType::Pointer pGeometry,
PropertiesType::Pointer pProperties)
GeometryType::Pointer pGeometry,
PropertiesType::Pointer pProperties)
: GeoTCondition<TDim, TNumNodes>(NewId, pGeometry, pProperties)
{
}
Expand All @@ -43,17 +43,16 @@ template <unsigned int TDim, unsigned int TNumNodes>
GeoTNormalFluxCondition<TDim, TNumNodes>::~GeoTNormalFluxCondition() = default;

template <unsigned int TDim, unsigned int TNumNodes>
void GeoTNormalFluxCondition<TDim, TNumNodes>::CalculateRHS(VectorType& rRightHandSideVector,
const ProcessInfo& rCurrentProcessInfo)
void GeoTNormalFluxCondition<TDim, TNumNodes>::CalculateRHS(Vector& rRightHandSideVector,
const ProcessInfo& rCurrentProcessInfo)
{
const GeometryType& r_geom = this->GetGeometry();
const GeometryType::IntegrationPointsArrayType& r_integration_points =
r_geom.IntegrationPoints(this->GetIntegrationMethod());
const unsigned int num_integration_points = r_integration_points.size();
const unsigned int local_dim = r_geom.LocalSpaceDimension();

const Matrix& r_N_container =
r_geom.ShapeFunctionsValues(this->GetIntegrationMethod());
const Matrix& r_N_container = r_geom.ShapeFunctionsValues(this->GetIntegrationMethod());
GeometryType::JacobiansType j_container(num_integration_points);

for (auto& x : j_container) {
Expand All @@ -66,40 +65,27 @@ void GeoTNormalFluxCondition<TDim, TNumNodes>::CalculateRHS(VectorType& rRightHa
for (unsigned int i = 0; i < TNumNodes; ++i) {
normal_flux_vector[i] = r_geom[i].FastGetSolutionStepValue(NORMAL_HEAT_FLUX);
}
NormalFluxVariables variables;

for (unsigned int integration_point = 0;
integration_point < num_integration_points; ++integration_point) {
for (unsigned int integration_point = 0; integration_point < num_integration_points; ++integration_point) {
// Obtain N
noalias(variables.N) = row(r_N_container, integration_point);
auto N = row(r_N_container, integration_point);

// Interpolation of nodal normal flux to integration point normal flux.
variables.NormalFluxOnNode = MathUtils<>::Dot(variables.N, normal_flux_vector);
auto normal_flux_on_integration_point = MathUtils<>::Dot(N, normal_flux_vector);

// Compute weighting coefficient for integration
this->CalculateIntegrationCoefficient(
variables.IntegrationCoefficient, j_container[integration_point],
r_integration_points[integration_point].Weight());
auto weighting_integration_coefficient = CalculateIntegrationCoefficient(j_container[integration_point],
r_integration_points[integration_point].Weight());

// Contributions to the right hand side
this->CalculateAndAddRHS(rRightHandSideVector, variables);
auto normal_flux_on_DOF = normal_flux_on_integration_point * N * weighting_integration_coefficient;
GeoElementUtilities::AssemblePBlockVector<0, TNumNodes>(rRightHandSideVector,
normal_flux_on_DOF);
}
}

template <unsigned int TDim, unsigned int TNumNodes>
void GeoTNormalFluxCondition<TDim, TNumNodes>::CalculateAndAddRHS(VectorType& rRightHandSideVector,
NormalFluxVariables& rVariables)
{
noalias(rVariables.NormalFluxOnIntegPoints) =
rVariables.NormalFluxOnNode * rVariables.N * rVariables.IntegrationCoefficient;

GeoElementUtilities::AssemblePBlockVector<0, TNumNodes>(
rRightHandSideVector, rVariables.NormalFluxOnIntegPoints);
}

template <unsigned int TDim, unsigned int TNumNodes>
void GeoTNormalFluxCondition<TDim, TNumNodes>::CalculateIntegrationCoefficient(
double& rIntegrationCoefficient, const Matrix& rJacobian, double Weight)
double GeoTNormalFluxCondition<TDim, TNumNodes>::CalculateIntegrationCoefficient(const Matrix& rJacobian,
double Weight)
{
Vector normal_vector = ZeroVector(TDim);

Expand All @@ -110,7 +96,7 @@ void GeoTNormalFluxCondition<TDim, TNumNodes>::CalculateIntegrationCoefficient(
MathUtils<double>::CrossProduct(normal_vector, column(rJacobian, 0),
column(rJacobian, 1));
}
rIntegrationCoefficient = Weight * MathUtils<double>::Norm(normal_vector);
return Weight * MathUtils<double>::Norm(normal_vector);
}

template class GeoTNormalFluxCondition<2, 2>;
Expand Down
26 changes: 6 additions & 20 deletions applications/GeoMechanicsApplication/custom_conditions/T_normal_flux_condition.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,8 +17,6 @@
#include "custom_conditions/T_condition.h"
#include "custom_utilities/condition_utilities.hpp"
#include "custom_utilities/element_utilities.hpp"
#include "geo_mechanics_application_variables.h"
#include "includes/condition.h"
#include "includes/serializer.h"

namespace Kratos {
Expand All @@ -27,11 +25,9 @@ template <unsigned int TDim, unsigned int TNumNodes>
class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoTNormalFluxCondition
: public GeoTCondition<TDim, TNumNodes> {
public:
using NodeType = Node;
using GeometryType = Geometry<NodeType>;
using GeometryType = Geometry<Node>;
using PropertiesType = Properties;
using NodesArrayType = GeometryType::PointsArrayType;
using VectorType = Vector;

KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION(GeoTNormalFluxCondition);

Expand All @@ -40,8 +36,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoTNormalFluxCondition
GeoTNormalFluxCondition(IndexType NewId, GeometryType::Pointer pGeometry);

GeoTNormalFluxCondition(IndexType NewId,
GeometryType::Pointer pGeometry,
PropertiesType::Pointer pProperties);
GeometryType::Pointer pGeometry,
PropertiesType::Pointer pProperties);

~GeoTNormalFluxCondition() override;

Expand All @@ -54,21 +50,11 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoTNormalFluxCondition
}

protected:
struct NormalFluxVariables {
double NormalFluxOnNode;
double IntegrationCoefficient;
array_1d<double, TNumNodes> N;
array_1d<double, TNumNodes> NormalFluxOnIntegPoints;
};

void CalculateRHS(VectorType& rRightHandSideVector,
void CalculateRHS(Vector& rRightHandSideVector,
const ProcessInfo& CurrentProcessInfo) override;

void CalculateAndAddRHS(VectorType& rRightHandSideVector, NormalFluxVariables& rVariables);

virtual void CalculateIntegrationCoefficient(double& rIntegrationCoefficient,
const Matrix& rJacobian,
double Weight);
double CalculateIntegrationCoefficient(const Matrix& rJacobian,
double Weight);

private:
friend class Serializer;
Expand Down

0 comments on commit e434e2c

Please sign in to comment.