Skip to content

Commit 06d5fdb

Browse files
avdg81avdg81
authored
[GeoMechanicsApplication] Nodal extrapolation for interface elements (#14061)
Extended the nodal extrapolator such that it also supports interface elements (line and surface ones). Changes include: - The interface geometry can now return its mid-geometry through the member functions `pGetGeometryPart` and `GetGeometryPart` (the provided part index is not being used). - Since `pGetGeometryPart` returns a shared pointer to the underlying mid-geometry, the data member that keeps the mid-geometry also had to be changed to a shared pointer. - Added overrides for members `GetGeometryFamily` and `GetGeometryOrderType` for class template `InterfaceGeometry`, which forward the request to the underlying mid-geometry. - Any interface element can now return a list of its integration points (regardless of the adopted integration scheme). - Any interface element can now return its mid-geometry. - When calculating the linear extrapolation matrix, we now only need to provide it with the element itself rather than the element's geometry and integration method. This change enables us to define a unified API that works for all elements (recall that the integration method could not be easily extended towards ones that are commonly used by interface elements). - Added a utility function that returns the integration points of any element, including interface elements. - Added a utility function that calculates the determinants of the Jacobian matrices of any element, including interface elements. - Extended the nodal linear extrapolator such that it also supports any interface element. - Made a utility function that generates nodes public. - Added element setup functions that create linear and quadratic quadrilateral elements. - Extended several existing unit tests to cover the newly added overrides for interface geometry (i.e. getting the geometry's family as well as its order). - Added new unit tests to cover the proper construction of the mid-geometry of any interface geometry. - The unit tests of the linear nodal extrapolator have been updated such that they work with elements rather than geometries. - Added a helper function that checks whether two given lists of points are equal to within a provided tolerance on a pair-by-pair basis. - Reuse the UBlas utility functions to create expected matrices and vectors.
1 parent 027c203 commit 06d5fdb

25 files changed

Lines changed: 646 additions & 358 deletions

‎applications/GeoMechanicsApplication/custom_conditions/U_Pw_normal_lysmer_absorbing_condition.cpp‎

Lines changed: 2 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,7 @@
1515
#include "custom_conditions/U_Pw_normal_lysmer_absorbing_condition.h"
1616
#include "custom_utilities/condition_utilities.hpp"
1717
#include "custom_utilities/dof_utilities.h"
18-
#include "custom_utilities/linear_nodal_extrapolator.h"
18+
#include "custom_utilities/extrapolation_utilities.h"
1919

2020
namespace Kratos
2121
{
@@ -276,14 +276,6 @@ void UPwLysmerAbsorbingCondition<TDim, TNumNodes>::CalculateNodalStiffnessMatrix
276276
}
277277
}
278278

279-
template <unsigned int TDim, unsigned int TNumNodes>
280-
Matrix UPwLysmerAbsorbingCondition<TDim, TNumNodes>::CalculateExtrapolationMatrixNeighbour(const Element& rNeighbourElement)
281-
{
282-
LinearNodalExtrapolator extrapolator;
283-
return extrapolator.CalculateElementExtrapolationMatrix(
284-
rNeighbourElement.GetGeometry(), rNeighbourElement.GetIntegrationMethod());
285-
}
286-
287279
template <unsigned int TDim, unsigned int TNumNodes>
288280
void UPwLysmerAbsorbingCondition<TDim, TNumNodes>::GetNeighbourElementVariables(
289281
NormalLysmerAbsorbingVariables& rVariables, const ProcessInfo& rCurrentProcessInfo)
@@ -337,7 +329,7 @@ void UPwLysmerAbsorbingCondition<TDim, TNumNodes>::GetNeighbourElementVariables(
337329
(1.0 - prop_neighbour[POROSITY]) * prop_neighbour[DENSITY_SOLID];
338330
}
339331

340-
Matrix extrapolation_matrix = CalculateExtrapolationMatrixNeighbour(r_neighbour_element);
332+
Matrix extrapolation_matrix = ExtrapolationUtilities::CalculateExtrapolationMatrix(r_neighbour_element);
341333

342334
// project parameters on neighbour nodes
343335
Vector Ec_nodes_neighbour = prod(extrapolation_matrix, confined_stiffness_vector);

‎applications/GeoMechanicsApplication/custom_conditions/U_Pw_normal_lysmer_absorbing_condition.h‎

Lines changed: 0 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -169,12 +169,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwLysmerAbsorbingCondition
169169
void CalculateNodalStiffnessMatrix(NormalLysmerAbsorbingVariables& rVariables,
170170
const Element::GeometryType& rGeom);
171171

172-
/**
173-
* @brief Calculates the extrapolation matrix for neighbour elements.Values from integration points are extrapolated to the nodes
174-
* @param rNeighbourElement The neighbouring element of the condition
175-
*/
176-
Matrix CalculateExtrapolationMatrixNeighbour(const Element& rNeighbourElement);
177-
178172
private:
179173
using hashmap =
180174
std::unordered_multimap<DenseVector<int>, std::vector<Condition::Pointer>, KeyHasherRange<DenseVector<int>>, KeyComparorRange<DenseVector<int>>>;

‎applications/GeoMechanicsApplication/custom_conditions/surface_load_3D_diff_order_condition.cpp‎

Lines changed: 2 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -94,10 +94,8 @@ void SurfaceLoad3DDiffOrderCondition::save(Serializer& rSerializer) const
9494
KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, GeneralUPwDiffOrderCondition)
9595
}
9696

97-
void SurfaceLoad3DDiffOrderCondition::load(Serializer& rSerializer)
98-
{
99-
KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, GeneralUPwDiffOrderCondition)
100-
}
97+
void SurfaceLoad3DDiffOrderCondition::load(Serializer& rSerializer){
98+
KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, GeneralUPwDiffOrderCondition)}
10199

102100
std::string SurfaceLoad3DDiffOrderCondition::Info() const
103101
{

‎applications/GeoMechanicsApplication/custom_conditions/surface_normal_fluid_flux_3D_diff_order_condition.cpp‎

Lines changed: 2 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -80,10 +80,8 @@ void SurfaceNormalFluidFlux3DDiffOrderCondition::save(Serializer& rSerializer) c
8080
KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, SurfaceLoad3DDiffOrderCondition)
8181
}
8282

83-
void SurfaceNormalFluidFlux3DDiffOrderCondition::load(Serializer& rSerializer)
84-
{
85-
KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, SurfaceLoad3DDiffOrderCondition)
86-
}
83+
void SurfaceNormalFluidFlux3DDiffOrderCondition::load(Serializer& rSerializer){
84+
KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, SurfaceLoad3DDiffOrderCondition)}
8785

8886
std::string SurfaceNormalFluidFlux3DDiffOrderCondition::Info() const
8987
{

‎applications/GeoMechanicsApplication/custom_conditions/surface_normal_load_3D_diff_order_condition.cpp‎

Lines changed: 2 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -107,10 +107,8 @@ void SurfaceNormalLoad3DDiffOrderCondition::save(Serializer& rSerializer) const
107107
KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, SurfaceLoad3DDiffOrderCondition)
108108
}
109109

110-
void SurfaceNormalLoad3DDiffOrderCondition::load(Serializer& rSerializer)
111-
{
112-
KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, SurfaceLoad3DDiffOrderCondition)
113-
}
110+
void SurfaceNormalLoad3DDiffOrderCondition::load(Serializer& rSerializer){
111+
KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, SurfaceLoad3DDiffOrderCondition)}
114112

115113
std::string SurfaceNormalLoad3DDiffOrderCondition::Info() const
116114
{

‎applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_FIC_element.cpp‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -22,8 +22,7 @@ namespace
2222
{
2323
auto CalculateSquareExtrapolationMatrix(const Kratos::Element* pElement)
2424
{
25-
const auto extrapolation_matrix = Kratos::ExtrapolationUtilities::CalculateExtrapolationMatrix(
26-
pElement->GetGeometry(), pElement->GetIntegrationMethod(), pElement->Id());
25+
const auto extrapolation_matrix = Kratos::ExtrapolationUtilities::CalculateExtrapolationMatrix(*pElement);
2726
if (extrapolation_matrix.size1() != extrapolation_matrix.size2()) {
2827
KRATOS_ERROR << "Extrapolation matrix is not square for element id " << pElement->Id() << std::endl;
2928
}

‎applications/GeoMechanicsApplication/custom_elements/interface_element.cpp‎

Lines changed: 22 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -92,10 +92,10 @@ InterfaceElement::InterfaceElement(IndexType NewId,
9292
void InterfaceElement::MakeIntegrationSchemeAndAssignFunction()
9393
{
9494
if (GetGeometry().LocalSpaceDimension() == 1) {
95-
mIntegrationScheme = std::make_unique<LobattoIntegrationScheme>(GetGeometry().PointsNumber() / 2);
95+
mpIntegrationScheme = std::make_unique<LobattoIntegrationScheme>(GetGeometry().PointsNumber() / 2);
9696
mfpCalculateRotationMatrix = GeometryUtilities::Calculate2DRotationMatrixForLineGeometry;
9797
} else {
98-
mIntegrationScheme = std::make_unique<LumpedIntegrationScheme>(GetGeometry().PointsNumber() / 2);
98+
mpIntegrationScheme = std::make_unique<LumpedIntegrationScheme>(GetGeometry().PointsNumber() / 2);
9999
mfpCalculateRotationMatrix = GeometryUtilities::Calculate3DRotationMatrixForPlaneGeometry;
100100
}
101101
}
@@ -187,7 +187,7 @@ void InterfaceElement::Initialize(const ProcessInfo& rCurrentProcessInfo)
187187
Element::Initialize(rCurrentProcessInfo);
188188

189189
mConstitutiveLaws.clear();
190-
for (auto i = std::size_t{0}; i < mIntegrationScheme->GetNumberOfIntegrationPoints(); ++i) {
190+
for (auto i = std::size_t{0}; i < mpIntegrationScheme->GetNumberOfIntegrationPoints(); ++i) {
191191
mConstitutiveLaws.push_back(GetProperties()[CONSTITUTIVE_LAW]->Clone());
192192
}
193193
// Only interpolate when neighbouring elements that provide nodal stresses were found
@@ -205,13 +205,12 @@ void InterfaceElement::Initialize(const ProcessInfo& rCurrentProcessInfo)
205205
r_neighbour_element.CalculateOnIntegrationPoints(
206206
TOTAL_STRESS_VECTOR, neighbour_cauchy_stresses, rCurrentProcessInfo);
207207
interface_nodal_cauchy_stresses = ExtrapolationUtilities::CalculateNodalVectors(
208-
interface_node_ids, r_neighbour_element.GetGeometry(), r_neighbour_element.GetIntegrationMethod(),
209-
neighbour_cauchy_stresses, r_neighbour_element.Id());
208+
interface_node_ids, r_neighbour_element, neighbour_cauchy_stresses);
210209
InterpolateNodalStressesToInitialTractions(interface_nodal_cauchy_stresses);
211210
}
212211
const auto shape_function_values_at_integration_points =
213212
GeoElementUtilities::EvaluateShapeFunctionsAtIntegrationPoints(
214-
mIntegrationScheme->GetIntegrationPoints(), GetGeometry());
213+
mpIntegrationScheme->GetIntegrationPoints(), GetGeometry());
215214
for (auto i = std::size_t{0}; i < mConstitutiveLaws.size(); ++i) {
216215
mConstitutiveLaws[i]->InitializeMaterial(GetProperties(), GetGeometry(),
217216
shape_function_values_at_integration_points[i]);
@@ -224,8 +223,8 @@ int InterfaceElement::Check(const ProcessInfo& rCurrentProcessInfo) const
224223
if (error != 0) return error;
225224

226225
if (this->IsActive()) {
227-
KRATOS_ERROR_IF(mIntegrationScheme->GetNumberOfIntegrationPoints() != mConstitutiveLaws.size())
228-
<< "Number of integration points (" << mIntegrationScheme->GetNumberOfIntegrationPoints()
226+
KRATOS_ERROR_IF(mpIntegrationScheme->GetNumberOfIntegrationPoints() != mConstitutiveLaws.size())
227+
<< "Number of integration points (" << mpIntegrationScheme->GetNumberOfIntegrationPoints()
229228
<< ") and constitutive laws (" << mConstitutiveLaws.size() << ") do not match.\n";
230229

231230
const auto r_properties = GetProperties();
@@ -239,6 +238,17 @@ int InterfaceElement::Check(const ProcessInfo& rCurrentProcessInfo) const
239238
return 0;
240239
}
241240

241+
const IntegrationScheme& InterfaceElement::GetIntegrationScheme() const
242+
{
243+
return *mpIntegrationScheme;
244+
}
245+
246+
const Geometry<Node>& InterfaceElement::GetMidGeometry() const
247+
{
248+
constexpr auto unused_part_index = std::size_t{0};
249+
return GetGeometry().GetGeometryPart(unused_part_index);
250+
}
251+
242252
Element::DofsVectorType InterfaceElement::GetDofs() const
243253
{
244254
const auto no_Pw_geometry_yet = Geometry<Node>{};
@@ -248,7 +258,7 @@ Element::DofsVectorType InterfaceElement::GetDofs() const
248258

249259
std::vector<Matrix> InterfaceElement::CalculateLocalBMatricesAtIntegrationPoints() const
250260
{
251-
const auto& r_integration_points = mIntegrationScheme->GetIntegrationPoints();
261+
const auto& r_integration_points = mpIntegrationScheme->GetIntegrationPoints();
252262
const auto shape_function_values_at_integration_points =
253263
GeoElementUtilities::EvaluateShapeFunctionsAtIntegrationPoints(r_integration_points, GetGeometry());
254264

@@ -274,8 +284,8 @@ std::vector<Matrix> InterfaceElement::CalculateLocalBMatricesAtIntegrationPoints
274284
std::vector<double> InterfaceElement::CalculateIntegrationCoefficients() const
275285
{
276286
const auto determinants_of_jacobian = CalculateDeterminantsOfJacobiansAtIntegrationPoints(
277-
mIntegrationScheme->GetIntegrationPoints(), GetGeometry());
278-
return mIntegrationCoefficientsCalculator.Run<>(mIntegrationScheme->GetIntegrationPoints(),
287+
mpIntegrationScheme->GetIntegrationPoints(), GetGeometry());
288+
return mIntegrationCoefficientsCalculator.Run<>(mpIntegrationScheme->GetIntegrationPoints(),
279289
determinants_of_jacobian, this);
280290
}
281291

@@ -355,7 +365,7 @@ void InterfaceElement::InterpolateNodalStressesToInitialTractions(const std::vec
355365
}
356366

357367
std::size_t integration_point_index = 0;
358-
for (const auto& r_integration_point : mIntegrationScheme->GetIntegrationPoints()) {
368+
for (const auto& r_integration_point : mpIntegrationScheme->GetIntegrationPoints()) {
359369
const auto integration_point_stress =
360370
InterpolateNodalStressToIntegrationPoints(r_integration_point, nodal_stresses);
361371

‎applications/GeoMechanicsApplication/custom_elements/interface_element.h‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -66,6 +66,10 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) InterfaceElement : public Element
6666
void Initialize(const ProcessInfo& rCurrentProcessInfo) override;
6767
int Check(const ProcessInfo& rCurrentProcessInfo) const override;
6868

69+
const IntegrationScheme& GetIntegrationScheme() const;
70+
71+
const Geometry<Node>& GetMidGeometry() const;
72+
6973
private:
7074
InterfaceElement() = default;
7175

@@ -88,7 +92,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) InterfaceElement : public Element
8892
Vector ConvertLocalStressToTraction(const Matrix& rLocalStress) const;
8993
std::function<Matrix(const Geometry<Node>&, const array_1d<double, 3>&)> mfpCalculateRotationMatrix;
9094

91-
std::unique_ptr<IntegrationScheme> mIntegrationScheme;
95+
std::unique_ptr<IntegrationScheme> mpIntegrationScheme;
9296
std::unique_ptr<StressStatePolicy> mpStressStatePolicy;
9397
std::vector<ConstitutiveLaw::Pointer> mConstitutiveLaws;
9498
IntegrationCoefficientsCalculator mIntegrationCoefficientsCalculator;

‎applications/GeoMechanicsApplication/custom_geometries/interface_geometry.h‎

Lines changed: 36 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -44,8 +44,8 @@ class InterfaceGeometry : public Geometry<Node>
4444
(rThisPoints.size() == 8) || (rThisPoints.size() == 16))
4545
<< "Number of nodes must be 2+2, 3+3, 6+6, 4+4 or 8+8\n";
4646

47-
mMidGeometry = std::make_unique<MidGeometryType>(CreatePointsOfMidGeometry());
48-
this->SetGeometryData(&mMidGeometry->GetGeometryData());
47+
mpMidGeometry = std::make_shared<MidGeometryType>(CreatePointsOfMidGeometry());
48+
this->SetGeometryData(&mpMidGeometry->GetGeometryData());
4949
}
5050

5151
[[nodiscard]] BaseType::Pointer Create(const PointsArrayType& rThisPoints) const override
@@ -59,62 +59,78 @@ class InterfaceGeometry : public Geometry<Node>
5959
return std::make_shared<InterfaceGeometry>(NewGeometryId, rThisPoints);
6060
}
6161

62-
[[nodiscard]] double Area() const override { return mMidGeometry->Area(); }
62+
[[nodiscard]] GeometryData::KratosGeometryFamily GetGeometryFamily() const override
63+
{
64+
return mpMidGeometry->GetGeometryFamily();
65+
}
66+
67+
[[nodiscard]] GeometryData::KratosGeometryOrderType GetGeometryOrderType() const override
68+
{
69+
return mpMidGeometry->GetGeometryOrderType();
70+
}
71+
72+
[[nodiscard]] double Area() const override { return mpMidGeometry->Area(); }
6373

6474
[[nodiscard]] double ShapeFunctionValue(IndexType ShapeFunctionIndex,
6575
const CoordinatesArrayType& rLocalCoordinate) const override
6676
{
67-
return mMidGeometry->ShapeFunctionValue(ShapeFunctionIndex, rLocalCoordinate);
77+
return mpMidGeometry->ShapeFunctionValue(ShapeFunctionIndex, rLocalCoordinate);
6878
}
6979

7080
Vector& ShapeFunctionsValues(Vector& rResult, const CoordinatesArrayType& rLocalCoordinate) const override
7181
{
72-
return mMidGeometry->ShapeFunctionsValues(rResult, rLocalCoordinate);
82+
return mpMidGeometry->ShapeFunctionsValues(rResult, rLocalCoordinate);
7383
}
7484

7585
Matrix& ShapeFunctionsLocalGradients(Matrix& rResult, const CoordinatesArrayType& rLocalCoordinate) const override
7686
{
77-
return mMidGeometry->ShapeFunctionsLocalGradients(rResult, rLocalCoordinate);
87+
return mpMidGeometry->ShapeFunctionsLocalGradients(rResult, rLocalCoordinate);
7888
}
7989

8090
Matrix& Jacobian(Matrix& rResult, const CoordinatesArrayType& rLocalCoordinate) const override
8191
{
82-
return mMidGeometry->Jacobian(rResult, rLocalCoordinate);
92+
return mpMidGeometry->Jacobian(rResult, rLocalCoordinate);
8393
}
8494

8595
[[nodiscard]] double DeterminantOfJacobian(const CoordinatesArrayType& rLocalCoordinate) const override
8696
{
87-
return mMidGeometry->DeterminantOfJacobian(rLocalCoordinate);
97+
return mpMidGeometry->DeterminantOfJacobian(rLocalCoordinate);
8898
}
8999

90100
Matrix& InverseOfJacobian(Matrix& rResult, const CoordinatesArrayType& rLocalCoordinate) const override
91101
{
92102
KRATOS_ERROR << "Inverse of Jacobian is not implemented for the interface geometry\n";
93103
}
94104

95-
[[nodiscard]] double Length() const override { return mMidGeometry->Length(); }
105+
[[nodiscard]] double Length() const override { return mpMidGeometry->Length(); }
96106

97-
[[nodiscard]] double DomainSize() const override { return mMidGeometry->DomainSize(); }
107+
[[nodiscard]] double DomainSize() const override { return mpMidGeometry->DomainSize(); }
98108

99109
[[nodiscard]] std::string Info() const override
100110
{
101-
return "An interface geometry consisting of two sub-geometries with Info: " + mMidGeometry->Info();
111+
return "An interface geometry consisting of two sub-geometries with Info: " + mpMidGeometry->Info();
102112
}
103113

104114
CoordinatesArrayType& PointLocalCoordinates(CoordinatesArrayType& rResult,
105115
const CoordinatesArrayType& rGlobalCoordinate) const override
106116
{
107-
return mMidGeometry->PointLocalCoordinates(rResult, rGlobalCoordinate);
117+
return mpMidGeometry->PointLocalCoordinates(rResult, rGlobalCoordinate);
108118
}
109119

110120
Matrix& PointsLocalCoordinates(Matrix& rResult) const override
111121
{
112-
return mMidGeometry->PointsLocalCoordinates(rResult);
122+
return mpMidGeometry->PointsLocalCoordinates(rResult);
123+
}
124+
125+
// The way for client code to access the mid-geometry
126+
const GeometryType::Pointer pGetGeometryPart(const IndexType) const override
127+
{
128+
return mpMidGeometry;
113129
}
114130

115131
void PrintInfo(std::ostream& rOStream) const override { rOStream << Info(); }
116132

117-
void PrintData(std::ostream& rOStream) const override { mMidGeometry->PrintData(rOStream); }
133+
void PrintData(std::ostream& rOStream) const override { mpMidGeometry->PrintData(rOStream); }
118134

119135
array_1d<double, 3> Normal(IndexType IntegrationPointIndex) const override
120136
{
@@ -207,7 +223,7 @@ class InterfaceGeometry : public Geometry<Node>
207223

208224
GeometriesArrayType GenerateEdges() const override
209225
{
210-
KRATOS_ERROR_IF_NOT(mMidGeometry->GetGeometryFamily() == GeometryData::KratosGeometryFamily::Kratos_Linear)
226+
KRATOS_ERROR_IF_NOT(mpMidGeometry->GetGeometryFamily() == GeometryData::KratosGeometryFamily::Kratos_Linear)
211227
<< "Edges can only be generated for line geometries. This is a surface interface "
212228
"geometry, which does not support edges.\n";
213229

@@ -216,7 +232,7 @@ class InterfaceGeometry : public Geometry<Node>
216232

217233
GeometriesArrayType GenerateFaces() const override
218234
{
219-
KRATOS_ERROR_IF(mMidGeometry->GetGeometryFamily() == GeometryData::KratosGeometryFamily::Kratos_Linear)
235+
KRATOS_ERROR_IF(mpMidGeometry->GetGeometryFamily() == GeometryData::KratosGeometryFamily::Kratos_Linear)
220236
<< "Faces can only be generated for surface geometries. This is a line "
221237
"interface geometry, which does not support faces.\n";
222238

@@ -225,7 +241,7 @@ class InterfaceGeometry : public Geometry<Node>
225241

226242
GeometriesArrayType GenerateBoundariesEntities() const override
227243
{
228-
switch (mMidGeometry->GetGeometryFamily()) {
244+
switch (mpMidGeometry->GetGeometryFamily()) {
229245
using enum GeometryData::KratosGeometryFamily;
230246
case Kratos_Linear:
231247
return this->GenerateEdges();
@@ -302,8 +318,8 @@ class InterfaceGeometry : public Geometry<Node>
302318
// The second side is defined by the second half of the element nodes. However, the
303319
// nodes must be traversed in opposite direction.
304320
auto nodes_of_second_side = PointerVector<Node>{begin_of_second_side, points.ptr_end()};
305-
GeometryUtilities::ReverseNodes(nodes_of_second_side, mMidGeometry->GetGeometryFamily(),
306-
mMidGeometry->GetGeometryOrderType());
321+
GeometryUtilities::ReverseNodes(nodes_of_second_side, this->GetGeometryFamily(),
322+
this->GetGeometryOrderType());
307323

308324
auto result = GeometriesArrayType{};
309325
result.push_back(std::make_shared<MidGeometryType>(nodes_of_first_side));
@@ -317,7 +333,7 @@ class InterfaceGeometry : public Geometry<Node>
317333
"schemes.\n";
318334
}
319335

320-
std::unique_ptr<BaseType> mMidGeometry;
336+
std::shared_ptr<BaseType> mpMidGeometry;
321337
};
322338

323339
} // namespace Kratos

0 commit comments

Comments
 (0)