Skip to content
Merged
Show file tree
Hide file tree
Changes from 19 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Original file line number Diff line number Diff line change
@@ -0,0 +1,27 @@
// KRATOS___
// // ) )
// // ___ ___
// // ____ //___) ) // ) )
// // / / // // / /
// ((____/ / ((____ ((___/ / MECHANICS
//
// License: geo_mechanics_application/license.txt
//
// Main authors: Mohamed Nabi
// John van Esch
// Gennady Markelov
//

#include "custom_elements/transient_Pw_line_element.h"

namespace Kratos
{

template class TransientPwLineElement<2, 2>;
template class TransientPwLineElement<2, 3>;
template class TransientPwLineElement<2, 4>;
template class TransientPwLineElement<2, 5>;
template class TransientPwLineElement<3, 2>;
template class TransientPwLineElement<3, 3>;

} // Namespace Kratos

Large diffs are not rendered by default.

Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,6 @@ class GeoElementUtilities
}
}

//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
template< unsigned int TDim, unsigned int TNumNodes >
static inline void InterpolateVariableWithComponents(array_1d<double, TDim>& rVector,
const Matrix& NContainer,
Expand All @@ -61,7 +60,6 @@ class GeoElementUtilities
}
}

//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
template< unsigned int TDof, unsigned int TNumNodes >
static inline void InterpolateVariableWithComponents(Vector& rVector,
const Matrix& NContainer,
Expand All @@ -81,7 +79,6 @@ class GeoElementUtilities
}
}

//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
static inline void FillArray1dOutput(array_1d<double,3>& rOutputValue,
const array_1d<double,2>& ComputedValue)
{
Expand All @@ -90,7 +87,7 @@ class GeoElementUtilities
rOutputValue[2] = 0.0;
}

//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

static inline void FillArray1dOutput(array_1d<double,3>& rOutputValue,
const array_1d<double,3>& ComputedValue)
{
Expand All @@ -99,7 +96,6 @@ class GeoElementUtilities
rOutputValue[2] = ComputedValue[2];
}

//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
template< unsigned int TDim, unsigned int TNumNodes >
static inline void GetNodalVariableVector(array_1d<double,TDim*TNumNodes>& rNodalVariableVector,
const Element::GeometryType& Geom,
Expand All @@ -116,7 +112,6 @@ class GeoElementUtilities
}
}

//----------------------------------------------------------------------------------------
template< unsigned int TDim, unsigned int TNumNodes >
static inline void GetNodalVariableMatrix(Matrix& rNodalVariableMatrix,
const Element::GeometryType& rGeom,
Expand All @@ -133,7 +128,14 @@ class GeoElementUtilities
}
}

//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
static inline void FillPermeabilityMatrix(BoundedMatrix<double, 1, 1>& rPermeabilityMatrix,
const Element::PropertiesType& Prop)
{
// 1D
rPermeabilityMatrix(0, 0) = Prop[PERMEABILITY_XX];
}


static inline void FillPermeabilityMatrix(BoundedMatrix<double,2,2>& rPermeabilityMatrix,
const Element::PropertiesType& Prop)
{
Expand All @@ -145,7 +147,6 @@ class GeoElementUtilities
rPermeabilityMatrix(1,0) = rPermeabilityMatrix(0,1);
}

//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
static inline void FillPermeabilityMatrix(BoundedMatrix<double,3,3>& rPermeabilityMatrix,
const Element::PropertiesType& Prop)
{
Expand All @@ -164,8 +165,6 @@ class GeoElementUtilities
rPermeabilityMatrix(0,2) = rPermeabilityMatrix(2,0);
}


//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
static inline void InvertMatrix2(BoundedMatrix<double,2,2>& rInvertedMatrix,
const BoundedMatrix<double,2,2>& InputMatrix,
double &InputMatrixDet)
Expand All @@ -188,7 +187,6 @@ class GeoElementUtilities
KRATOS_CATCH("")
}

//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
static inline void InvertMatrix2(BoundedMatrix<double,2,2>& rInvertedMatrix,
const BoundedMatrix<double,2,2>& InputMatrix)
{
Expand All @@ -201,7 +199,6 @@ class GeoElementUtilities
KRATOS_CATCH("")
}

//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
template< unsigned int TDim>
static inline void AssembleDensityMatrix(BoundedMatrix<double,TDim+1, TDim+1> &DensityMatrix,
double Density)
Expand Down Expand Up @@ -269,7 +266,6 @@ class GeoElementUtilities
AddVectorAtPosition(rPBlockVector, rRightHandSideVector, offset);
}

//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
static inline void CalculateNewtonCotesLocalShapeFunctionsGradients(BoundedMatrix<double,2,2>& DN_DeContainer)
{
//Line 2-noded
Expand All @@ -286,7 +282,6 @@ class GeoElementUtilities
}
}

//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
static inline void CalculateNewtonCotesLocalShapeFunctionsGradients(BoundedMatrix<double,3,3>& DN_DeContainer)
{
//Line 3-noded
Expand All @@ -304,7 +299,6 @@ class GeoElementUtilities
}
}

//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
static inline void CalculateNewtonCotesShapeFunctions(BoundedMatrix<double,3,3>& NContainer)
{
//Line 3-noded
Expand All @@ -319,7 +313,6 @@ class GeoElementUtilities
}
}

//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
static inline void CalculateEquallyDistributedPointsLineShapeFunctions3N(Matrix& NContainer)
{
//Line 3-noded
Expand All @@ -339,7 +332,6 @@ class GeoElementUtilities
}
}

//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
static inline void CalculateEquallyDistributedPointsLineGradientShapeFunctions3N(GeometryData::ShapeFunctionsGradientsType& DN_DeContainer)
{
//Line 3-noded
Expand Down Expand Up @@ -388,8 +380,6 @@ class GeoElementUtilities
return Circumference;
}


//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
static inline void CalculateExtrapolationMatrixTriangle(Matrix& rExtrapolationMatrix, const GeometryData::IntegrationMethod& rIntegrationMethod)
{
/// The matrix contains the shape functions at each GP evaluated at each node.
Expand Down Expand Up @@ -427,7 +417,6 @@ class GeoElementUtilities

}

//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
static inline void CalculateExtrapolationMatrixQuad(Matrix& rExtrapolationMatrix, const GeometryData::IntegrationMethod& rIntegrationMethod)
{
if (rIntegrationMethod == GeometryData::IntegrationMethod::GI_GAUSS_1)
Expand Down Expand Up @@ -462,7 +451,6 @@ class GeoElementUtilities
}
}


static inline void CalculateExtrapolationMatrixTetra(Matrix& rExtrapolationMatrix, const GeometryData::IntegrationMethod& rIntegrationMethod)
{
if (rIntegrationMethod == GeometryData::IntegrationMethod::GI_GAUSS_1)
Expand Down Expand Up @@ -498,8 +486,6 @@ class GeoElementUtilities

}


//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
static inline void CalculateExtrapolationMatrixHexa(Matrix& rExtrapolationMatrix, const GeometryData::IntegrationMethod& rIntegrationMethod)
{
if (rIntegrationMethod == GeometryData::IntegrationMethod::GI_GAUSS_1)
Expand Down Expand Up @@ -557,9 +543,6 @@ class GeoElementUtilities
}
}


//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

static Vector CalculateNodalHydraulicHeadFromWaterPressures(const GeometryType& rGeom, const Properties& rProp)
{
const auto NumericalLimit = std::numeric_limits<double>::epsilon();
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,13 @@ void KratosGeoMechanicsApplication::Register() {
KRATOS_REGISTER_ELEMENT("TransientPwElement3D20N", mTransientPwElement3D20N)
KRATOS_REGISTER_ELEMENT("TransientPwElement3D27N", mTransientPwElement3D27N)

KRATOS_REGISTER_ELEMENT("TransientPwLineElement2D2N", mTransientPwLineElement2D2N)
KRATOS_REGISTER_ELEMENT("TransientPwLineElement2D3N", mTransientPwLineElement2D3N)
KRATOS_REGISTER_ELEMENT("TransientPwLineElement2D4N", mTransientPwLineElement2D4N)
KRATOS_REGISTER_ELEMENT("TransientPwLineElement2D5N", mTransientPwLineElement2D5N)
KRATOS_REGISTER_ELEMENT("TransientPwLineElement3D2N", mTransientPwLineElement3D2N)
KRATOS_REGISTER_ELEMENT("TransientPwLineElement3D3N", mTransientPwLineElement3D3N)

KRATOS_REGISTER_ELEMENT("TransientPwInterfaceElement2D4N", mTransientPwInterfaceElement2D4N)
KRATOS_REGISTER_ELEMENT("TransientPwInterfaceElement3D6N", mTransientPwInterfaceElement3D6N)
KRATOS_REGISTER_ELEMENT("TransientPwInterfaceElement3D8N", mTransientPwInterfaceElement3D8N)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -94,6 +94,7 @@
#include "custom_elements/transient_Pw_interface_element.hpp"
#include "custom_elements/undrained_U_Pw_small_strain_element.hpp"
#include "custom_elements/updated_lagrangian_U_Pw_diff_order_element.hpp"
#include "custom_elements/transient_Pw_line_element.h"

// Element policies
#include "custom_elements/axisymmetric_stress_state.h"
Expand Down Expand Up @@ -288,6 +289,13 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) KratosGeoMechanicsApplication : publ
const TransientPwElement<3,20> mTransientPwElement3D20N{ 0, Kratos::make_shared< Hexahedra3D20 <NodeType> >(Element::GeometryType::PointsArrayType(20)), std::make_unique<ThreeDimensionalStressState>() };
const TransientPwElement<3,27> mTransientPwElement3D27N{ 0, Kratos::make_shared< Hexahedra3D27 <NodeType> >(Element::GeometryType::PointsArrayType(27)), std::make_unique<ThreeDimensionalStressState>() };

const TransientPwLineElement<2, 2> mTransientPwLineElement2D2N{ 0, Kratos::make_shared<Line2D2<NodeType>>(Element::GeometryType::PointsArrayType(2)) };
const TransientPwLineElement<2, 3> mTransientPwLineElement2D3N{ 0, Kratos::make_shared<Line2D3<NodeType>>(Element::GeometryType::PointsArrayType(3)) };
const TransientPwLineElement<2, 4> mTransientPwLineElement2D4N{ 0, Kratos::make_shared<Line2D4<NodeType>>(Element::GeometryType::PointsArrayType(4)) };
const TransientPwLineElement<2, 5> mTransientPwLineElement2D5N{ 0, Kratos::make_shared<Line2D5<NodeType>>(Element::GeometryType::PointsArrayType(5)) };
const TransientPwLineElement<3, 2> mTransientPwLineElement3D2N{ 0, Kratos::make_shared<Line3D2<NodeType>>(Element::GeometryType::PointsArrayType(2))};
const TransientPwLineElement<3, 3> mTransientPwLineElement3D3N{ 0, Kratos::make_shared<Line3D3<NodeType>>(Element::GeometryType::PointsArrayType(3))};

const TransientPwInterfaceElement<2,4> mTransientPwInterfaceElement2D4N{ 0, Kratos::make_shared< QuadrilateralInterface2D4 <NodeType> >(Element::GeometryType::PointsArrayType(4)), std::make_unique<PlaneStrainStressState>() };
const TransientPwInterfaceElement<3,6> mTransientPwInterfaceElement3D6N{ 0, Kratos::make_shared< PrismInterface3D6 <NodeType> >(Element::GeometryType::PointsArrayType(6)), std::make_unique<ThreeDimensionalStressState>() };
const TransientPwInterfaceElement<3,8> mTransientPwInterfaceElement3D8N{ 0, Kratos::make_shared< HexahedraInterface3D8 <NodeType> >(Element::GeometryType::PointsArrayType(8)), std::make_unique<ThreeDimensionalStressState>() };
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@
from test_prescribed_derivatives import KratosGeoMechanicsPrescribedDerivatives
from test_dirichlet_u import KratosGeoMechanicsDirichletUTests
from test_normal_load_on_hexa_element import KratosGeoMechanicsNormalLoadHexaTests
from test_pressure_line_element import KratosGeoMechanicsTransientPressureLineElementTests

def AssembleTestSuites():
''' Populates the test suites to run.
Expand Down Expand Up @@ -104,7 +105,8 @@ def AssembleTestSuites():
TestSellmeijersRule,
TestElementaryGroundWaterFlow,
KratosGeoMechanicsTransientThermalTests,
KratosGeoMechanicsTimeIntegrationTests
KratosGeoMechanicsTimeIntegrationTests,
KratosGeoMechanicsTransientPressureLineElementTests
]
night_test_cases.extend(small_test_cases)

Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,63 @@
import os

import KratosMultiphysics.KratosUnittest as KratosUnittest
import test_helper

class KratosGeoMechanicsTransientPressureLineElementTests(KratosUnittest.TestCase):
"""
This class contains benchmark tests which are checked with the regression on a previously obtained value.
"""
etalon_value1 = -20000.0

def setUp(self):
# Code here will be placed BEFORE every test in this TestCase.
pass

def tearDown(self):
# Code here will be placed AFTER every test in this TestCase.
pass

def check_water_pressure(self, test_name, etalon_value):
file_path = test_helper.get_file_path(os.path.join('test_pressure_line_element', test_name))
simulation = test_helper.run_kratos(file_path)
pressure = test_helper.get_water_pressure(simulation)
self.assertAlmostEqual(etalon_value, pressure[2])

def test_oblique_line_element2D2N(self):
self.check_water_pressure("test_oblique_line_element2D2N", self.etalon_value1)

def test_oblique_line_element2D3N(self):
self.check_water_pressure("test_oblique_line_element2D3N", self.etalon_value1)

def test_oblique_line_element2D4N(self):
self.check_water_pressure("test_oblique_line_element2D4N", self.etalon_value1)

def test_oblique_line_element2D5N(self):
self.check_water_pressure("test_oblique_line_element2D5N", self.etalon_value1)

def test_vertical_line_element2D2N(self):
self.check_water_pressure("test_vertical_line_element2D2N", self.etalon_value1)

def test_vertical_line_element2D3N(self):
self.check_water_pressure("test_vertical_line_element2D3N", self.etalon_value1)

def test_vertical_line_element2D4N(self):
self.check_water_pressure("test_vertical_line_element2D4N", self.etalon_value1)

def test_vertical_line_element2D5N(self):
self.check_water_pressure("test_vertical_line_element2D5N", self.etalon_value1)

def test_oblique_line_element3D2N(self):
self.check_water_pressure("test_oblique_line_element3D2N", self.etalon_value1)

def test_oblique_line_element3D3N(self):
self.check_water_pressure("test_oblique_line_element3D3N", self.etalon_value1)

def test_vertical_line_element3D2N(self):
self.check_water_pressure("test_vertical_line_element3D2N", self.etalon_value1)

def test_vertical_line_element3D3N(self):
self.check_water_pressure("test_vertical_line_element3D3N", self.etalon_value1)

if __name__ == '__main__':
KratosUnittest.main()
Original file line number Diff line number Diff line change
@@ -0,0 +1,29 @@
{
"properties": [{
"model_part_name" : "PorousDomain.filter",
"properties_id" : 1,
"Material" : {
"constitutive_law" : {
"name" : "GeoLinearElasticPlaneStrain2DLaw"
},
"Variables" : {
"IGNORE_UNDRAINED" : false,
"YOUNG_MODULUS" : 1.000000e+07,
"POISSON_RATIO" : 0.000000e+00,
"DENSITY_SOLID" : 2.650000e+03,
"DENSITY_WATER" : 1.000000e+03,
"POROSITY" : 1.000000e-01,
"BULK_MODULUS_SOLID" : 9.000000e+19,
"BULK_MODULUS_FLUID" : 1.000000e+20,
"PERMEABILITY_XX" : 9.084000e-06,
"DYNAMIC_VISCOSITY" : 1.000000e-03,
"BIOT_COEFFICIENT" : 1.000000e+00,
"RETENTION_LAW" : "SaturatedLaw",
"SATURATED_SATURATION" : 1.000000e+00,
"RESIDUAL_SATURATION" : 0.000000e+00
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is RESIDUAL_SATURATION needed for the SaturatedLaw? The 3D materials parameters files shows that it is possible without.

},
"Tables": {
}
}
}]
}
Original file line number Diff line number Diff line change
@@ -0,0 +1,28 @@
{
"properties": [{
"model_part_name" : "PorousDomain.filter",
"properties_id" : 1,
"Material" : {
"constitutive_law" : {
"name" : "LinearElasticK03DLaw"
},
"Variables" : {
"IGNORE_UNDRAINED" : false,
"YOUNG_MODULUS" : 1.000000e+07,
"POISSON_RATIO" : 0.000000e+00,
"DENSITY_SOLID" : 2.650000e+03,
"DENSITY_WATER" : 1.000000e+03,
"POROSITY" : 1.000000e-01,
"BULK_MODULUS_SOLID" : 9.000000e+19,
"BULK_MODULUS_FLUID" : 1.000000e+20,
"PERMEABILITY_XX" : 9.084000e-06,
"DYNAMIC_VISCOSITY" : 1.000000e-03,
"BIOT_COEFFICIENT" : 1.000000e+00,
"RETENTION_LAW" : "SaturatedLaw",
"SATURATED_SATURATION" : 1.000000e+00
},
"Tables": {
}
}
}]
}
Loading