Skip to content

Commit

Permalink
Merge branch 'master' into fluid/div-stable-element-p2-p1
Browse files Browse the repository at this point in the history
  • Loading branch information
@rubenzorrilla
rubenzorrilla committed May 14, 2024
2 parents 35e5cf8 + c2a05d7 commit ba4589e
Show file tree
Hide file tree
Showing 214 changed files with 125,280 additions and 118,781 deletions.
1 change: 1 addition & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,6 +80,7 @@ kratos/scripts/wheels/README.md

# Most used editor files
*.vscode
*.idx
remote_*.md

/compile_commands.json
Expand Down
3 changes: 2 additions & 1 deletion applications/GeoMechanicsApplication/.git-blame-ignore-revs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,4 +3,5 @@
ab13f86f8770ccaaa23ef69c63f14f443da8bea0
206d114e8082ce70ebf375e2a4d4e84f366c6696
f051af1a6a180c67b27c0deb29062999629ededc
3ff1502c8993f1ed937f4016045dcb26dd7e271d
3ff1502c8993f1ed937f4016045dcb26dd7e271d
8776deb6c34974b51c2cf2c53d347bf55ec583d6
2 changes: 2 additions & 0 deletions applications/GeoMechanicsApplication/custom_conditions/T_condition.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,10 +93,12 @@ Condition::DofsVectorType GeoTCondition<TDim, TNumNodes>::GetDofs() const
return Geo::DofUtilities::ExtractDofsFromNodes(this->GetGeometry(), TEMPERATURE);
}

template class GeoTCondition<2, 1>;
template class GeoTCondition<2, 2>;
template class GeoTCondition<2, 3>;
template class GeoTCondition<2, 4>;
template class GeoTCondition<2, 5>;
template class GeoTCondition<3, 1>;
template class GeoTCondition<3, 3>;
template class GeoTCondition<3, 4>;
template class GeoTCondition<3, 6>;
Expand Down
3 changes: 1 addition & 2 deletions ...GeoMechanicsApplication/custom_conditions/surface_normal_load_3D_diff_order_condition.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,10 +12,9 @@
// Vahid Galavi
//


// Project includes
#include "custom_conditions/surface_normal_load_3D_diff_order_condition.hpp"
#include "custom_utilities/math_utilities.hpp"
#include "custom_utilities/math_utilities.h"

namespace Kratos
{
Expand Down
50 changes: 50 additions & 0 deletions applications/GeoMechanicsApplication/custom_conditions/thermal_point_flux_condition.cpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,50 @@
// KRATOS___
// // ) )
// // ___ ___
// // ____ //___) ) // ) )
// // / / // // / /
// ((____/ / ((____ ((___/ / MECHANICS
//
// License: geo_mechanics_application/license.txt
//
//
// Main authors: Mohamed Nabi
// John van Esch
//

#include "custom_conditions/thermal_point_flux_condition.h"

namespace Kratos
{

template <unsigned int TDim, unsigned int TNumNodes>
GeoThermalPointFluxCondition<TDim, TNumNodes>::GeoThermalPointFluxCondition()
: GeoTCondition<TDim, TNumNodes>()
{
}

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

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

template <unsigned int TDim, unsigned int TNumNodes>
void GeoThermalPointFluxCondition<TDim, TNumNodes>::CalculateRHS(Vector& rRightHandSideVector, const ProcessInfo&)
{
rRightHandSideVector[0] = this->GetGeometry()[0].FastGetSolutionStepValue(NORMAL_HEAT_FLUX);
}

template class GeoThermalPointFluxCondition<2, 1>;
template class GeoThermalPointFluxCondition<3, 1>;

} // namespace Kratos
64 changes: 64 additions & 0 deletions applications/GeoMechanicsApplication/custom_conditions/thermal_point_flux_condition.h
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,64 @@
// KRATOS___
// // ) )
// // ___ ___
// // ____ //___) ) // ) )
// // / / // // / /
// ((____/ / ((____ ((___/ / MECHANICS
//
// License: geo_mechanics_application/license.txt
//
//
// Main authors: Mohamed Nabi
// John van Esch
//

#pragma once

#include "custom_conditions/T_condition.h"
#include "includes/serializer.h"

namespace Kratos
{

template <unsigned int TDim, unsigned int TNumNodes>
class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoThermalPointFluxCondition
: public GeoTCondition<TDim, TNumNodes>
{
public:
using GeometryType = Geometry<Node>;
using PropertiesType = Properties;
using NodesArrayType = GeometryType::PointsArrayType;
using BaseType = GeoTCondition<TDim, TNumNodes>;

KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION(GeoThermalPointFluxCondition);

GeoThermalPointFluxCondition();

GeoThermalPointFluxCondition(IndexType NewId, GeometryType::Pointer pGeometry);

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

Condition::Pointer Create(IndexType NewId, const NodesArrayType& rThisNodes, PropertiesType::Pointer pProperties) const override
{
return Kratos::make_intrusive<GeoThermalPointFluxCondition>(
NewId, this->GetGeometry().Create(rThisNodes), pProperties);
}

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

private:
friend class Serializer;

void save(Serializer& rSerializer) const override
{
KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, BaseType)
}

void load(Serializer& rSerializer) override
{
KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, BaseType)
}
};

} // namespace Kratos
50 changes: 25 additions & 25 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
// Application includes
#include "custom_elements/U_Pw_base_element.hpp"
#include "custom_utilities/dof_utilities.h"
#include "custom_utilities/equation_of_motion_utilities.h"

namespace Kratos
{
Expand Down Expand Up @@ -322,31 +323,19 @@ void UPwBaseElement<TDim, TNumNodes>::CalculateDampingMatrix(MatrixType&
{
KRATOS_TRY

// Rayleigh Method: Damping Matrix = alpha*M + beta*K

const unsigned int N_DOF = this->GetNumberOfDOF();

// Compute Mass Matrix
MatrixType MassMatrix(N_DOF, N_DOF);

this->CalculateMassMatrix(MassMatrix, rCurrentProcessInfo);

// Compute Stiffness matrix
MatrixType StiffnessMatrix(N_DOF, N_DOF);
MatrixType mass_matrix(N_DOF, N_DOF);
this->CalculateMassMatrix(mass_matrix, rCurrentProcessInfo);

this->CalculateMaterialStiffnessMatrix(StiffnessMatrix, rCurrentProcessInfo);

// Compute Damping Matrix
if (rDampingMatrix.size1() != N_DOF) rDampingMatrix.resize(N_DOF, N_DOF, false);
noalias(rDampingMatrix) = ZeroMatrix(N_DOF, N_DOF);

const PropertiesType& rProp = this->GetProperties();
MatrixType stiffness_matrix(N_DOF, N_DOF);
this->CalculateMaterialStiffnessMatrix(stiffness_matrix, rCurrentProcessInfo);

if (rProp.Has(RAYLEIGH_ALPHA)) noalias(rDampingMatrix) += rProp[RAYLEIGH_ALPHA] * MassMatrix;
else noalias(rDampingMatrix) += rCurrentProcessInfo[RAYLEIGH_ALPHA] * MassMatrix;

if (rProp.Has(RAYLEIGH_BETA)) noalias(rDampingMatrix) += rProp[RAYLEIGH_BETA] * StiffnessMatrix;
else noalias(rDampingMatrix) += rCurrentProcessInfo[RAYLEIGH_BETA] * StiffnessMatrix;
const PropertiesType& r_prop = this->GetProperties();
rDampingMatrix = GeoEquationOfMotionUtilities::CalculateDampingMatrix(
r_prop.Has(RAYLEIGH_ALPHA) ? r_prop[RAYLEIGH_ALPHA] : rCurrentProcessInfo[RAYLEIGH_ALPHA],
r_prop.Has(RAYLEIGH_BETA) ? r_prop[RAYLEIGH_BETA] : rCurrentProcessInfo[RAYLEIGH_BETA],
mass_matrix, stiffness_matrix);

KRATOS_CATCH("")
}
Expand Down Expand Up @@ -527,12 +516,23 @@ void UPwBaseElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLeftHandS

//----------------------------------------------------------------------------------------
template <unsigned int TDim, unsigned int TNumNodes>
double UPwBaseElement<TDim, TNumNodes>::CalculateIntegrationCoefficient(
const GeometryType::IntegrationPointsArrayType& IntegrationPoints, unsigned int PointNumber, double detJ)
double UPwBaseElement<TDim, TNumNodes>::CalculateIntegrationCoefficient(const GeometryType::IntegrationPointType& rIntegrationPoint,
double detJ) const

{
return mpStressStatePolicy->CalculateIntegrationCoefficient(IntegrationPoints[PointNumber],
detJ, GetGeometry());
return mpStressStatePolicy->CalculateIntegrationCoefficient(rIntegrationPoint, detJ, GetGeometry());
}

template <unsigned int TDim, unsigned int TNumNodes>
std::vector<double> UPwBaseElement<TDim, TNumNodes>::CalculateIntegrationCoefficients(
const GeometryType::IntegrationPointsArrayType& rIntegrationPoints, const Vector& rDetJs) const
{
auto result = std::vector<double>{};
std::transform(rIntegrationPoints.begin(), rIntegrationPoints.end(), rDetJs.begin(),
std::back_inserter(result), [this](const auto& rIntegrationPoint, const auto& rDetJ) {
return this->CalculateIntegrationCoefficient(rIntegrationPoint, rDetJ);
});
return result;
}

//----------------------------------------------------------------------------------------
Expand Down
7 changes: 4 additions & 3 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -164,9 +164,10 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwBaseElement : public Element
bool CalculateStiffnessMatrixFlag,
bool CalculateResidualVectorFlag);

virtual double CalculateIntegrationCoefficient(const GeometryType::IntegrationPointsArrayType& IntegrationPoints,
unsigned int PointNumber,
double detJ);
double CalculateIntegrationCoefficient(const GeometryType::IntegrationPointType& rIntegrationPoint,
double detJ) const;
std::vector<double> CalculateIntegrationCoefficients(const GeometryType::IntegrationPointsArrayType& rIntegrationPoints,
const Vector& rDetJs) const;

void CalculateDerivativesOnInitialConfiguration(
double& detJ, Matrix& J0, Matrix& InvJ0, Matrix& DN_DX, unsigned int PointNumber) const;
Expand Down
45 changes: 38 additions & 7 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_FIC_element.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@

// Application includes
#include "custom_elements/U_Pw_small_strain_FIC_element.hpp"
#include "custom_utilities/math_utilities.h"

namespace Kratos
{
Expand Down Expand Up @@ -136,12 +137,22 @@ void UPwSmallStrainFICElement<TDim, TNumNodes>::InitializeNonLinearIteration(con

Vector StressVector(VoigtSize);

const auto b_matrices = this->CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer);
const auto deformation_gradients = this->CalculateDeformationGradients();
const auto determinants_of_deformation_gradients =
GeoMechanicsMathUtilities::CalculateDeterminants(deformation_gradients);
const auto strain_vectors = this->CalculateStrains(
deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain);

// Loop over integration points
for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
this->CalculateKinematics(Variables, GPoint);
Variables.B = b_matrices[GPoint];

// Compute infinitessimal strain
this->CalculateStrain(Variables, GPoint);
Variables.F = deformation_gradients[GPoint];
Variables.detF = determinants_of_deformation_gradients[GPoint];
Variables.StrainVector = strain_vectors[GPoint];

// set gauss points variables to constitutivelaw parameters
this->SetConstitutiveParameters(Variables, ConstitutiveParameters);
Expand Down Expand Up @@ -187,12 +198,22 @@ void UPwSmallStrainFICElement<TDim, TNumNodes>::FinalizeNonLinearIteration(const

Vector StressVector(VoigtSize);

const auto b_matrices = this->CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer);
const auto deformation_gradients = this->CalculateDeformationGradients();
const auto determinants_of_deformation_gradients =
GeoMechanicsMathUtilities::CalculateDeterminants(deformation_gradients);
const auto strain_vectors = this->CalculateStrains(
deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain);

// Loop over integration points
for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
this->CalculateKinematics(Variables, GPoint);
Variables.B = b_matrices[GPoint];

// Compute infinitessimal strain
this->CalculateStrain(Variables, GPoint);
Variables.F = deformation_gradients[GPoint];
Variables.detF = determinants_of_deformation_gradients[GPoint];
Variables.StrainVector = strain_vectors[GPoint];

// set gauss points variables to constitutivelaw parameters
this->SetConstitutiveParameters(Variables, ConstitutiveParameters);
Expand Down Expand Up @@ -438,13 +459,25 @@ void UPwSmallStrainFICElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLeftHa

const bool hasBiotCoefficient = Prop.Has(BIOT_COEFFICIENT);

const auto b_matrices = this->CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer);
const auto integration_coefficients =
this->CalculateIntegrationCoefficients(IntegrationPoints, Variables.detJContainer);
const auto deformation_gradients = this->CalculateDeformationGradients();
const auto determinants_of_deformation_gradients =
GeoMechanicsMathUtilities::CalculateDeterminants(deformation_gradients);
const auto strain_vectors = this->CalculateStrains(
deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain);

// Loop over integration points
for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
// Compute Np, GradNpT, B and StrainVector
this->CalculateKinematics(Variables, GPoint);
Variables.B = b_matrices[GPoint];

// Compute infinitessimal strain
this->CalculateStrain(Variables, GPoint);
Variables.F = deformation_gradients[GPoint];
Variables.detF = determinants_of_deformation_gradients[GPoint];
Variables.StrainVector = strain_vectors[GPoint];

// set gauss points variables to constitutivelaw parameters
this->SetConstitutiveParameters(Variables, ConstitutiveParameters);
Expand All @@ -468,12 +501,10 @@ void UPwSmallStrainFICElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLeftHa
// calculate Bulk modulus from stiffness matrix
this->InitializeBiotCoefficients(Variables, hasBiotCoefficient);

// Compute weighting coefficient for integration
Variables.IntegrationCoefficient =
this->CalculateIntegrationCoefficient(IntegrationPoints, GPoint, Variables.detJ);
Variables.IntegrationCoefficient = integration_coefficients[GPoint];

Variables.IntegrationCoefficientInitialConfiguration = this->CalculateIntegrationCoefficient(
IntegrationPoints, GPoint, Variables.detJInitialConfiguration);
IntegrationPoints[GPoint], Variables.detJInitialConfiguration);

if (CalculateStiffnessMatrixFlag) {
// Contributions to the left hand side
Expand Down
Loading

0 comments on commit ba4589e

Please sign in to comment.