Skip to content

Commit

Permalink
[GeoMechanicsApplication] Remove swapping faces in condition process (#…
Browse files Browse the repository at this point in the history 
…12227)

* Removed face swapping for Triangle3D3, and 3D6, Fixed the existing test cases
* Removed face swapping for Quadrilateral3D4, and Quadrilateral3D8, Added 2 test cases to thermal
* Fix sign of the normal stresses for 3D cases
* Added new test to the test_GeoMechanicsApplication.py script
---------
Co-authored-by: mnabideltares <[email protected]>
  • Loading branch information
@rfaasse
rfaasse authored Mar 28, 2024
1 parent f84694f commit ca1560e
Show file tree
Hide file tree
Showing 28 changed files with 5,408 additions and 172 deletions.
6 changes: 5 additions & 1 deletion applications/GeoMechanicsApplication/custom_conditions/U_Pw_normal_face_load_condition.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,7 @@ void UPwNormalFaceLoadCondition<TDim, TNumNodes>::CalculateTractionVector(
const unsigned int& GPoint)
{
Vector NormalVector = ZeroVector(TDim);
const double NormalStress = MathUtils<>::Dot(row(NContainer, GPoint), Variables.NormalStressVector);
double NormalStress = MathUtils<>::Dot(row(NContainer, GPoint), Variables.NormalStressVector);

if constexpr (TDim == 2) {
const double TangentialStress = MathUtils<>::Dot(row(NContainer, GPoint), Variables.TangentialStressVector);
Expand All @@ -107,6 +107,10 @@ void UPwNormalFaceLoadCondition<TDim, TNumNodes>::CalculateTractionVector(
rTractionVector[1] = NormalStress * NormalVector[0] + TangentialStress * NormalVector[1];
}
else if constexpr (TDim == 3) {
// Since the normal vector is pointing outwards for the 3D conditions, the normal stress
// should switch sign, such that positive normal contact stress is defined inwards.
NormalStress *= -1;

MathUtils<double>::CrossProduct(NormalVector, column(Jacobian, 0), column(Jacobian, 1));
rTractionVector = NormalStress * NormalVector;
}
Expand Down
42 changes: 17 additions & 25 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 @@ -15,6 +15,7 @@

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

namespace Kratos
{
Expand All @@ -36,37 +37,28 @@ Condition::Pointer SurfaceNormalLoad3DDiffOrderCondition::Create(IndexType NewId
return Condition::Pointer(new SurfaceNormalLoad3DDiffOrderCondition(NewId, GetGeometry().Create(ThisNodes), pProperties));
}

//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
void SurfaceNormalLoad3DDiffOrderCondition::
CalculateConditionVector(ConditionVariables& rVariables, unsigned int PointNumber)
void SurfaceNormalLoad3DDiffOrderCondition::CalculateConditionVector(ConditionVariables& rVariables,
unsigned int PointNumber)
{
KRATOS_TRY

double NormalVector[3];

NormalVector[0] = rVariables.JContainer[PointNumber](1,0) * rVariables.JContainer[PointNumber](2,1) -
rVariables.JContainer[PointNumber](2,0) * rVariables.JContainer[PointNumber](1,1);

NormalVector[1] = rVariables.JContainer[PointNumber](2,0) * rVariables.JContainer[PointNumber](0,1) -
rVariables.JContainer[PointNumber](0,0) * rVariables.JContainer[PointNumber](2,1);
Vector normal_vector(3);
MathUtils<double>::CrossProduct(normal_vector, column(rVariables.JContainer[PointNumber], 0),
column(rVariables.JContainer[PointNumber], 1));

NormalVector[2] = rVariables.JContainer[PointNumber](0,0) * rVariables.JContainer[PointNumber](1,1) -
rVariables.JContainer[PointNumber](1,0) * rVariables.JContainer[PointNumber](0,1);

const GeometryType& rGeom = GetGeometry();
const SizeType NumUNodes = rGeom.PointsNumber();
double NormalStress = 0.0;
rVariables.ConditionVector.resize(3,false);
const auto& r_geometry = GetGeometry();

for ( SizeType i = 0; i < NumUNodes; ++i ) {
NormalStress += rVariables.Nu[i]*rGeom[i].FastGetSolutionStepValue(NORMAL_CONTACT_STRESS);
}
Vector normal_stresses(r_geometry.PointsNumber());
std::transform(r_geometry.begin(), r_geometry.end(), normal_stresses.begin(), [](const auto& r_node) {
return r_node.FastGetSolutionStepValue(NORMAL_CONTACT_STRESS);
});

rVariables.ConditionVector[0] = NormalStress * NormalVector[0];
rVariables.ConditionVector[1] = NormalStress * NormalVector[1];
rVariables.ConditionVector[2] = NormalStress * NormalVector[2];
// Since the normal vector is pointing outwards for the 3D conditions, the normal stress
// should switch sign, such that positive normal contact stress is defined inwards.
const double normal_stress = -1 * MathUtils<>::Dot(rVariables.Nu, normal_stresses);
rVariables.ConditionVector = normal_stress * normal_vector;

KRATOS_CATCH( "" )
KRATOS_CATCH("")
}

//----------------------------------------------------------------------------------------
Expand All @@ -92,7 +84,7 @@ void SurfaceNormalLoad3DDiffOrderCondition::

for ( SizeType i = 0; i < NumUNodes; ++i ) {
Index = i * 3;

rRightHandSideVector[Index] += rVariables.Nu[i] * rVariables.ConditionVector[0] * rVariables.IntegrationCoefficient;
rRightHandSideVector[Index+1] += rVariables.Nu[i] * rVariables.ConditionVector[1] * rVariables.IntegrationCoefficient;
rRightHandSideVector[Index+2] += rVariables.Nu[i] * rVariables.ConditionVector[2] * rVariables.IntegrationCoefficient;
Expand Down
24 changes: 0 additions & 24 deletions ...eoMechanicsApplication/custom_processes/find_neighbour_elements_of_conditions_process.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,30 +36,6 @@ void FindNeighbourElementsOfConditionsProcess::Execute()
Ids[i] = rGeometry[i].Id();
}

if (rGeometry.GetGeometryType() == GeometryData::KratosGeometryType::Kratos_Triangle3D3) {
// reverse ordering to be consistent with face ordering
std::swap(Ids[1], Ids[2]);
}

if (rGeometry.GetGeometryType() == GeometryData::KratosGeometryType::Kratos_Triangle3D6) {
// reverse ordering to be consistent with face ordering
std::swap(Ids[1], Ids[2]);
std::swap(Ids[3], Ids[5]);
}

if (rGeometry.GetGeometryType() == GeometryData::KratosGeometryType::Kratos_Quadrilateral3D4) {
// reverse ordering to be consistent with face ordering
std::swap(Ids[0], Ids[3]);
std::swap(Ids[1], Ids[2]);
}

if (rGeometry.GetGeometryType() == GeometryData::KratosGeometryType::Kratos_Quadrilateral3D8) {
// reverse ordering to be consistent with face ordering
std::swap(Ids[0], Ids[3]);
std::swap(Ids[1], Ids[2]);
std::swap(Ids[4], Ids[6]);
}

// adds to the map
FacesMap.insert( hashmap::value_type(Ids, std::vector<Condition::Pointer>({*itCond.base()})) );

Expand Down
29 changes: 29 additions & 0 deletions applications/GeoMechanicsApplication/tests/hexa_8n_normal_load/MaterialParameters.json
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,29 @@
{
"properties": [{
"model_part_name": "PorousDomain.Soil",
"properties_id": 1,
"Material": {
"constitutive_law": {
"name" : "LinearElastic3DLaw"
},
"Variables": {
"IGNORE_UNDRAINED" : true,
"YOUNG_MODULUS" : 3.0e7,
"POISSON_RATIO" : 0,
"DENSITY_SOLID" : 0,
"DENSITY_WATER" : 1.0e3,
"POROSITY" : 0.3,
"BULK_MODULUS_SOLID" : 1.0e16,
"BULK_MODULUS_FLUID" : 2.0e-30,
"PERMEABILITY_XX" : 4.5e-30,
"PERMEABILITY_YY" : 4.5e-30,
"PERMEABILITY_ZZ" : 4.5e-30,
"PERMEABILITY_XY" : 0.0,
"PERMEABILITY_YZ" : 0.0,
"PERMEABILITY_ZX" : 0.0,
"DYNAMIC_VISCOSITY" : 1.0e-3
},
"Tables": {}
}
}]
}
Loading

0 comments on commit ca1560e

Please sign in to comment.