From 00ae3649eb1debd8a34723a19672a748fce71a5c Mon Sep 17 00:00:00 2001 From: Vicente Mataix Ferrandiz Date: Fri, 3 Nov 2023 23:46:03 +0100 Subject: [PATCH 1/4] Define the `PointLocalCoordinatesTetrahedra3D4N` --- kratos/utilities/geometry_utilities.h | 108 ++++++++++++++++++++++++++ 1 file changed, 108 insertions(+) diff --git a/kratos/utilities/geometry_utilities.h b/kratos/utilities/geometry_utilities.h index a5980e289950..7b8a019cf941 100644 --- a/kratos/utilities/geometry_utilities.h +++ b/kratos/utilities/geometry_utilities.h @@ -55,6 +55,114 @@ class KRATOS_API(KRATOS_CORE) GeometryUtils */ static std::string GetGeometryName(const GeometryData::KratosGeometryType TypeOfGeometry); + /** + * @brief Returns the local coordinates of a given arbitrary point for a given linear tetrahedra + * @details Based on https://www.colorado.edu/engineering/CAS/courses.d/AFEM.d/AFEM.Ch09.d/AFEM.Ch09.pdf. Section 9.1.6 + * @param rGeometry The geometry to be considered + * @param rResult The vector containing the local coordinates of the point + * @param rPoint The point in global coordinates + * @return The vector containing the local coordinates of the point + */ + template + static inline typename TGeometryType::CoordinatesArrayType& PointLocalCoordinatesTetrahedra3D4N( + const TGeometryType& rGeometry, + typename TGeometryType::CoordinatesArrayType& rResult, + const typename TGeometryType::CoordinatesArrayType& rPoint + ) + { + // Compute RHS + array_1d X; + X[0] = 1.0; + X[1] = rPoint[0]; + X[2] = rPoint[1]; + X[3] = rPoint[2]; + + // Auxiliary coordinates + const auto& r_coordinates_0 = rGeometry[0].Coordinates(); + const auto& r_coordinates_1 = rGeometry[1].Coordinates(); + const auto& r_coordinates_2 = rGeometry[2].Coordinates(); + const auto& r_coordinates_3 = rGeometry[3].Coordinates(); + const double x1 = r_coordinates_0[0]; + const double y1 = r_coordinates_0[1]; + const double z1 = r_coordinates_0[2]; + const double x2 = r_coordinates_1[0]; + const double y2 = r_coordinates_1[1]; + const double z2 = r_coordinates_1[2]; + const double x3 = r_coordinates_2[0]; + const double y3 = r_coordinates_2[1]; + const double z3 = r_coordinates_2[2]; + const double x4 = r_coordinates_3[0]; + const double y4 = r_coordinates_3[1]; + const double z4 = r_coordinates_3[2]; + + // Auxiliary diff + const double x12 = x1 - x2; + const double x13 = x1 - x3; + const double x14 = x1 - x4; + const double x21 = x2 - x1; + const double x24 = x2 - x4; + const double x31 = x3 - x1; + const double x32 = x3 - x2; + const double x34 = x3 - x4; + const double x42 = x4 - x2; + const double x43 = x4 - x3; + const double y12 = y1 - y2; + const double y13 = y1 - y3; + const double y14 = y1 - y4; + const double y21 = y2 - y1; + const double y24 = y2 - y4; + const double y31 = y3 - y1; + const double y32 = y3 - y2; + const double y34 = y3 - y4; + const double y42 = y4 - y2; + const double y43 = y4 - y3; + const double z12 = z1 - z2; + const double z13 = z1 - z3; + const double z14 = z1 - z4; + const double z21 = z2 - z1; + const double z24 = z2 - z4; + const double z31 = z3 - z1; + const double z32 = z3 - z2; + const double z34 = z3 - z4; + const double z42 = z4 - z2; + const double z43 = z4 - z3; + + // Compute LHS + BoundedMatrix invJ; + const double aux_volume = 1.0/(6.0*rGeometry.Volume()); + invJ(0,0) = aux_volume * (x2*(y3*z4-y4*z3)+x3*(y4*z2-y2*z4)+x4*(y2*z3-y3*z2)); + invJ(1,0) = aux_volume * (x1*(y4*z3-y3*z4)+x3*(y1*z4-y4*z1)+x4*(y3*z1-y1*z3)); + invJ(2,0) = aux_volume * (x1*(y2*z4-y4*z2)+x2*(y4*z1-y1*z4)+x4*(y1*z2-y2*z1)); + invJ(3,0) = aux_volume * (x1*(y3*z2-y2*z3)+x2*(y1*z3-y3*z1)+x3*(y2*z1-y1*z2)); + invJ(0,1) = aux_volume * (y42*z32 - y32*z42); + invJ(1,1) = aux_volume * (y31*z43 - y34*z13); + invJ(2,1) = aux_volume * (y24*z14 - y14*z24); + invJ(3,1) = aux_volume * (y13*z21 - y12*z31); + invJ(0,2) = aux_volume * (x32*z42 - x42*z32); + invJ(1,2) = aux_volume * (x43*z31 - x13*z34); + invJ(2,2) = aux_volume * (x14*z24 - x24*z14); + invJ(3,2) = aux_volume * (x21*z13 - x31*z12); + invJ(0,3) = aux_volume * (x42*y32 - x32*y42); + invJ(1,3) = aux_volume * (x31*y43 - x34*y13); + invJ(2,3) = aux_volume * (x24*y14 - x14*y24); + invJ(3,3) = aux_volume * (x13*y21 - x12*y31); + + // Compute result + const array_1d result = prod(invJ, X); + + // Resize if needed + if (rResult.size() != 3) { + rResult.resize(3,false); + } + + // Copy result + rResult[0] = result[1]; + rResult[1] = result[2]; + rResult[2] = result[3]; + + return rResult; + } + /** * @brief This function is designed to compute the shape function derivatives, shape functions and volume in 3D * @param rGeometry it is the array of nodes. It is expected to be a tetrahedra From fcc66ef8fdb214962eca58776c5cf0e2684d586e Mon Sep 17 00:00:00 2001 From: Vicente Mataix Ferrandiz Date: Fri, 3 Nov 2023 23:46:20 +0100 Subject: [PATCH 2/4] Replace `PointLocalCoordinatesTetrahedra3D4N` in `Tetrahedra3D4N` --- kratos/geometries/tetrahedra_3d_4.h | 85 +---------------------------- 1 file changed, 1 insertion(+), 84 deletions(-) diff --git a/kratos/geometries/tetrahedra_3d_4.h b/kratos/geometries/tetrahedra_3d_4.h index 1f16842282e9..0eb1be06ab4b 100644 --- a/kratos/geometries/tetrahedra_3d_4.h +++ b/kratos/geometries/tetrahedra_3d_4.h @@ -961,7 +961,6 @@ template class Tetrahedra3D4 : public Geometry /** * @brief Returns the local coordinates of a given arbitrary point - * @details Based on https://www.colorado.edu/engineering/CAS/courses.d/AFEM.d/AFEM.Ch09.d/AFEM.Ch09.pdf. Section 9.1.6 * @param rResult The vector containing the local coordinates of the point * @param rPoint The point in global coordinates * @return The vector containing the local coordinates of the point @@ -971,89 +970,7 @@ template class Tetrahedra3D4 : public Geometry const CoordinatesArrayType& rPoint ) const override { - // Compute RHS - array_1d X; - X[0] = 1.0; - X[1] = rPoint[0]; - X[2] = rPoint[1]; - X[3] = rPoint[2]; - - // Auxiliar coordinates - const double x1 = this->GetPoint( 0 ).X(); - const double x2 = this->GetPoint( 1 ).X(); - const double x3 = this->GetPoint( 2 ).X(); - const double x4 = this->GetPoint( 3 ).X(); - const double y1 = this->GetPoint( 0 ).Y(); - const double y2 = this->GetPoint( 1 ).Y(); - const double y3 = this->GetPoint( 2 ).Y(); - const double y4 = this->GetPoint( 3 ).Y(); - const double z1 = this->GetPoint( 0 ).Z(); - const double z2 = this->GetPoint( 1 ).Z(); - const double z3 = this->GetPoint( 2 ).Z(); - const double z4 = this->GetPoint( 3 ).Z(); - - // Auxiliar diff - const double x12 = x1 - x2; - const double x13 = x1 - x3; - const double x14 = x1 - x4; - const double x21 = x2 - x1; - const double x24 = x2 - x4; - const double x31 = x3 - x1; - const double x32 = x3 - x2; - const double x34 = x3 - x4; - const double x42 = x4 - x2; - const double x43 = x4 - x3; - const double y12 = y1 - y2; - const double y13 = y1 - y3; - const double y14 = y1 - y4; - const double y21 = y2 - y1; - const double y24 = y2 - y4; - const double y31 = y3 - y1; - const double y32 = y3 - y2; - const double y34 = y3 - y4; - const double y42 = y4 - y2; - const double y43 = y4 - y3; - const double z12 = z1 - z2; - const double z13 = z1 - z3; - const double z14 = z1 - z4; - const double z21 = z2 - z1; - const double z24 = z2 - z4; - const double z31 = z3 - z1; - const double z32 = z3 - z2; - const double z34 = z3 - z4; - const double z42 = z4 - z2; - const double z43 = z4 - z3; - - // Compute LHS - BoundedMatrix invJ; - const double aux_volume = 1.0/(6.0*this->Volume()); - invJ(0,0) = aux_volume * (x2*(y3*z4-y4*z3)+x3*(y4*z2-y2*z4)+x4*(y2*z3-y3*z2)); - invJ(1,0) = aux_volume * (x1*(y4*z3-y3*z4)+x3*(y1*z4-y4*z1)+x4*(y3*z1-y1*z3)); - invJ(2,0) = aux_volume * (x1*(y2*z4-y4*z2)+x2*(y4*z1-y1*z4)+x4*(y1*z2-y2*z1)); - invJ(3,0) = aux_volume * (x1*(y3*z2-y2*z3)+x2*(y1*z3-y3*z1)+x3*(y2*z1-y1*z2)); - invJ(0,1) = aux_volume * (y42*z32 - y32*z42); - invJ(1,1) = aux_volume * (y31*z43 - y34*z13); - invJ(2,1) = aux_volume * (y24*z14 - y14*z24); - invJ(3,1) = aux_volume * (y13*z21 - y12*z31); - invJ(0,2) = aux_volume * (x32*z42 - x42*z32); - invJ(1,2) = aux_volume * (x43*z31 - x13*z34); - invJ(2,2) = aux_volume * (x14*z24 - x24*z14); - invJ(3,2) = aux_volume * (x21*z13 - x31*z12); - invJ(0,3) = aux_volume * (x42*y32 - x32*y42); - invJ(1,3) = aux_volume * (x31*y43 - x34*y13); - invJ(2,3) = aux_volume * (x24*y14 - x14*y24); - invJ(3,3) = aux_volume * (x13*y21 - x12*y31); - - const array_1d result = prod(invJ, X); - - if (rResult.size() != 3) - rResult.resize(3, false); - - rResult[0] = result[1]; - rResult[1] = result[2]; - rResult[2] = result[3]; - - return rResult; + return GeometryUtils::PointLocalCoordinatesTetrahedra3D4N(*this, rResult, rPoint); } /** From 35a65207717d58be1ccaaa3a1e8958a75faa035b Mon Sep 17 00:00:00 2001 From: Vicente Mataix Ferrandiz Date: Fri, 3 Nov 2023 23:46:52 +0100 Subject: [PATCH 3/4] Use `PointLocalCoordinatesTetrahedra3D4N` in `Tetrahedra3D10` when planar face --- kratos/geometries/tetrahedra_3d_10.h | 337 ++++++++++++--------------- 1 file changed, 150 insertions(+), 187 deletions(-) diff --git a/kratos/geometries/tetrahedra_3d_10.h b/kratos/geometries/tetrahedra_3d_10.h index 13f61baa0742..b133b7f225b3 100644 --- a/kratos/geometries/tetrahedra_3d_10.h +++ b/kratos/geometries/tetrahedra_3d_10.h @@ -24,11 +24,30 @@ // Project includes #include "geometries/triangle_3d_6.h" #include "geometries/tetrahedra_3d_4.h" +#include "utilities/geometry_utilities.h" #include "utilities/integration_utilities.h" #include "integration/tetrahedron_gauss_legendre_integration_points.h" namespace Kratos { +///@name Kratos Globals +///@{ + +///@} +///@name Type Definitions +///@{ + +///@} +///@name Enum's +///@{ + +///@} +///@name Functions +///@{ + +///@} +///@name Kratos Classes +///@{ /** * @class Tetrahedra3D10 * @ingroup KratosCore @@ -50,157 +69,96 @@ namespace Kratos * @author Janosch Stascheit * @author Felix Nagel */ -template class Tetrahedra3D10 : public Geometry +template +class Tetrahedra3D10 + : public Geometry { public: - /** - * Type Definitions - */ + ///@name Type Definitions + ///@{ - /** - * Geometry as base class. - */ - typedef Geometry BaseType; + /// Geometry as base class. + using BaseType = Geometry; - /** - * Type of edge and face geometries - */ - typedef Line3D3 EdgeType; - typedef Triangle3D6 FaceType; + /// Type of edge and face geometries + using EdgeType = Line3D3; + using FaceType = Triangle3D6; - /** - * Pointer definition of Tetrahedra3D10 - */ + /// Pointer definition of Tetrahedra3D10 KRATOS_CLASS_POINTER_DEFINITION( Tetrahedra3D10 ); - /** - * Integration methods implemented in geometry. - */ - typedef GeometryData::IntegrationMethod IntegrationMethod; - - /** - * A Vector of counted pointers to Geometries. Used for - * returning edges of the geometry. - */ - typedef typename BaseType::GeometriesArrayType GeometriesArrayType; + /// Integration methods implemented in geometry. + using IntegrationMethod = GeometryData::IntegrationMethod; - /** - * Redefinition of template parameter TPointType. - */ - typedef TPointType PointType; + /// A Vector of counted pointers to Geometries. Used for + /// returning edges of the geometry. + using GeometriesArrayType = typename BaseType::GeometriesArrayType; - /** - * Type used for indexing in geometry class.std::size_t used for indexing - * point or integration point access methods and also all other - * methods which need point or integration point index. - */ - typedef typename BaseType::IndexType IndexType; + /// Redefinition of template parameter TPointType. + using PointType = TPointType; + /// Type used for indexing in geometry class. std::size_t used for indexing + /// point or integration point access methods and also all other + /// methods which need point or integration point index. + using IndexType = typename BaseType::IndexType; - /** - * This typed used to return size or dimension in - * geometry. Dimension, WorkingDimension, PointsNumber and - * ... return this type as their results. - */ - typedef typename BaseType::SizeType SizeType; - - /** - * Array of counted pointers to point. This type used to hold - * geometry's points. - */ - typedef typename BaseType::PointsArrayType PointsArrayType; - - /** - * This type used for representing an integration point in - * geometry. This integration point is a point with an - * additional weight component. - */ - typedef typename BaseType::IntegrationPointType IntegrationPointType; + /// This type used to return size or dimension in + /// geometry. Dimension, WorkingDimension, PointsNumber, and + /// ... return this type as their results. + using SizeType = typename BaseType::SizeType; - /** - * A Vector of IntegrationPointType which used to hold - * integration points related to an integration - * method. IntegrationPoints functions used this type to return - * their results. - */ - typedef typename BaseType::IntegrationPointsArrayType IntegrationPointsArrayType; + /// Array of counted pointers to point. This type used to hold + /// geometry's points. + using PointsArrayType = typename BaseType::PointsArrayType; - /** - * A Vector of IntegrationPointsArrayType which used to hold - * integration points related to different integration method - * implemented in geometry. - */ - typedef typename BaseType::IntegrationPointsContainerType IntegrationPointsContainerType; + /// This type used for representing an integration point in + /// geometry. This integration point is a point with an + /// additional weight component. + using IntegrationPointType = typename BaseType::IntegrationPointType; - /** - * A third order tensor used as shape functions' values - * container. - */ - typedef typename BaseType::ShapeFunctionsValuesContainerType - ShapeFunctionsValuesContainerType; + /// A Vector of IntegrationPointType which used to hold + /// integration points related to an integration + /// method. IntegrationPoints functions used this type to return + /// their results. + using IntegrationPointsArrayType = typename BaseType::IntegrationPointsArrayType; - /** - * A fourth order tensor used as shape functions' local - * gradients container in geometry. - */ - typedef typename BaseType::ShapeFunctionsLocalGradientsContainerType - ShapeFunctionsLocalGradientsContainerType; + /// A Vector of IntegrationPointsArrayType which used to hold + /// integration points related to different integration method + /// implemented in geometry. + using IntegrationPointsContainerType = typename BaseType::IntegrationPointsContainerType; - /** - * A third order tensor to hold jacobian matrices evaluated at - * integration points. Jacobian and InverseOfJacobian functions - * return this type as their result. - */ - typedef typename BaseType::JacobiansType JacobiansType; + /// A third order tensor used as shape functions' values + /// container. + using ShapeFunctionsValuesContainerType = typename BaseType::ShapeFunctionsValuesContainerType; - /** - * A third order tensor to hold shape functions' local - * gradients. ShapefunctionsLocalGradients function return this - * type as its result. - */ - typedef typename BaseType::ShapeFunctionsGradientsType ShapeFunctionsGradientsType; + /// A fourth order tensor used as shape functions' local + /// gradients container in geometry. + using ShapeFunctionsLocalGradientsContainerType = typename BaseType::ShapeFunctionsLocalGradientsContainerType; - /** - * Type of the normal vector used for normal to edges in geomety. - */ - typedef typename BaseType::NormalType NormalType; + /// A third order tensor to hold Jacobian matrices evaluated at + /// integration points. Jacobian and InverseOfJacobian functions + /// return this type as their result. + using JacobiansType = typename BaseType::JacobiansType; - /** - * Type of coordinates array - */ - typedef typename BaseType::CoordinatesArrayType CoordinatesArrayType; + /// A third order tensor to hold shape functions' local + /// gradients. ShapefunctionsLocalGradients function return this + /// type as its result. + using ShapeFunctionsGradientsType = typename BaseType::ShapeFunctionsGradientsType; - /** - * Type of Matrix - */ - typedef Matrix MatrixType; + /// Type of the normal vector used for normal to edges in geometry. + using NormalType = typename BaseType::NormalType; + /// Type of coordinates array + using CoordinatesArrayType = typename BaseType::CoordinatesArrayType; - /** - * Life Cycle - */ + /// Type of Matrix + using MatrixType = Matrix; -// Tetrahedra3D10( const PointType& Point1, const PointType& Point2, -// const PointType& Point3, const PointType& Point4 , -// const PointType& Point5, const PointType& Point6, -// const PointType& Point7, const PointType& Point8 , -// const PointType& Point9, const PointType& Point10 -// ) -// : BaseType( PointsArrayType(), &msGeometryData ) -// { -// this->Points().reserve( 10 ); -// this->Points().push_back( typename PointType::Pointer( new PointType( Point1 ) ) ); -// this->Points().push_back( typename PointType::Pointer( new PointType( Point2 ) ) ); -// this->Points().push_back( typename PointType::Pointer( new PointType( Point3 ) ) ); -// this->Points().push_back( typename PointType::Pointer( new PointType( Point4 ) ) ); -// this->Points().push_back( typename PointType::Pointer( new PointType( Point5 ) ) ); -// this->Points().push_back( typename PointType::Pointer( new PointType( Point6 ) ) ); -// this->Points().push_back( typename PointType::Pointer( new PointType( Point7 ) ) ); -// this->Points().push_back( typename PointType::Pointer( new PointType( Point8 ) ) ); -// this->Points().push_back( typename PointType::Pointer( new PointType( Point9 ) ) ); -// this->Points().push_back( typename PointType::Pointer( new PointType( Point10 ) ) ); -// } + ///@} + ///@name Life Cycle + ///@{ + /// Default constructor. Tetrahedra3D10( typename PointType::Pointer pPoint1, typename PointType::Pointer pPoint2, typename PointType::Pointer pPoint3, @@ -227,11 +185,14 @@ template class Tetrahedra3D10 : public Geometry this->Points().push_back( pPoint10 ); } + /** + * @brief Constructor with points array + * @param ThisPoints The points of the current geometry + */ Tetrahedra3D10( const PointsArrayType& ThisPoints ) : BaseType( ThisPoints, &msGeometryData ) { - if ( this->PointsNumber() != 10 ) - KRATOS_ERROR << "Invalid points number. Expected 10, given " << this->PointsNumber() << std::endl; + KRATOS_ERROR_IF( this->PointsNumber() != 10 ) << "Invalid points number. Expected 10, given " << this->PointsNumber() << std::endl; } /// Constructor with Geometry Id @@ -296,9 +257,9 @@ template class Tetrahedra3D10 : public Geometry return GeometryData::KratosGeometryType::Kratos_Tetrahedra3D10; } - /** - * Operators - */ + ///@} + ///@name Operators + ///@{ /** * Assignment operator. @@ -336,7 +297,6 @@ template class Tetrahedra3D10 : public Geometry return *this; } - ///@} ///@name Operations ///@{ @@ -418,7 +378,7 @@ template class Tetrahedra3D10 : public Geometry */ double Length() const override { - return sqrt( fabs( this->DeterminantOfJacobian( PointType() ) ) ); + return std::sqrt( std::abs( this->DeterminantOfJacobian( PointType() ) ) ); } /** @@ -463,12 +423,29 @@ template class Tetrahedra3D10 : public Geometry * @see Length() * @see Area() * @see Volume() - * - * :TODO: might be necessary to reimplement */ double DomainSize() const override { - return Volume(); + return Volume(); + } + + /** + * @brief Returns the local coordinates of a given arbitrary point + * @param rResult The vector containing the local coordinates of the point + * @param rPoint The point in global coordinates + * @return The vector containing the local coordinates of the point + */ + CoordinatesArrayType& PointLocalCoordinates( + CoordinatesArrayType& rResult, + const CoordinatesArrayType& rPoint + ) const override + { + // Using linear approximation for planar faces + if (this->FacesArePlanar()) { + return GeometryUtils::PointLocalCoordinatesTetrahedra3D4N(*this, rResult, rPoint); + } else { + return BaseType::PointLocalCoordinates( rResult, rPoint ); + } } /** @@ -487,14 +464,10 @@ template class Tetrahedra3D10 : public Geometry { this->PointLocalCoordinates( rResult, rPoint ); - if ( (rResult[0] >= (0.0 - Tolerance)) && (rResult[0] <= (1.0 + Tolerance)) ) - { - if ( (rResult[1] >= (0.0 - Tolerance)) && (rResult[1] <= (1.0 + Tolerance)) ) - { - if ( (rResult[2] >= (0.0 - Tolerance)) && (rResult[2] <= (1.0 + Tolerance)) ) - { - if ((( 1.0 - ( rResult[0] + rResult[1] + rResult[2] ) ) >= (0.0 - Tolerance) ) && (( 1.0 - ( rResult[0] + rResult[1] + rResult[2] ) ) <= (1.0 + Tolerance) ) ) - { + if ( (rResult[0] >= (0.0 - Tolerance)) && (rResult[0] <= (1.0 + Tolerance)) ) { + if ( (rResult[1] >= (0.0 - Tolerance)) && (rResult[1] <= (1.0 + Tolerance)) ) { + if ( (rResult[2] >= (0.0 - Tolerance)) && (rResult[2] <= (1.0 + Tolerance)) ) { + if ((( 1.0 - ( rResult[0] + rResult[1] + rResult[2] ) ) >= (0.0 - Tolerance) ) && (( 1.0 - ( rResult[0] + rResult[1] + rResult[2] ) ) <= (1.0 + Tolerance) ) ) { return true; } } @@ -675,10 +648,9 @@ template class Tetrahedra3D10 : public Geometry return rResult; } - - /** - * Shape Function - */ + ///@} + ///@name Shape Function + ///@{ Vector& ShapeFunctionsValues(Vector &rResult, const CoordinatesArrayType& rCoordinates) const override { @@ -781,7 +753,9 @@ template class Tetrahedra3D10 : public Geometry */ bool HasIntersection(const Point& rLowPoint, const Point& rHighPoint) const override { + // Check if the faces are planar if (this->FacesArePlanar()) { + // TODO: Move implementation to GeometryUtils to avoid creating a new tetrahedra return Tetrahedra3D4( this->pGetPoint(0), this->pGetPoint(1), @@ -793,7 +767,6 @@ template class Tetrahedra3D10 : public Geometry return false; } - ///@} ///@name Spatial Operations ///@{ @@ -884,23 +857,14 @@ template class Tetrahedra3D10 : public Geometry rOStream << " Jacobian in the origin\t : " << jacobian; } } - -protected: - - /** - * there are no protected class members - */ - private: + ///@name Static Member Variables + ///@{ - /** - * Static Member Variables - */ static const GeometryData msGeometryData; static const GeometryDimension msGeometryDimension; - ///@} ///@name Serialization ///@{ @@ -917,12 +881,15 @@ template class Tetrahedra3D10 : public Geometry KRATOS_SERIALIZE_LOAD_BASE_CLASS( rSerializer, BaseType ); } - Tetrahedra3D10(): BaseType( PointsArrayType(), &msGeometryData ) {} + ///@} + ///@name Private Life Cycle + ///@{ + Tetrahedra3D10(): BaseType( PointsArrayType(), &msGeometryData ) {} - /** - * Private Operations - */ + ///@} + ///@name Private Operations + ///@{ /** * @brief Returns vector of shape function values at local coordinate. @@ -1101,11 +1068,9 @@ template class Tetrahedra3D10 : public Geometry } /** - * Checks if faces are planar. We iterate for all edges and check + * @brief Checks if faces are planar. We iterate for all edges and check * that the sum of 0-2 and 2-1 segments is no bigger than 0-1. - * * @return bool faces are planar or not - * */ bool FacesArePlanar() const { @@ -1121,35 +1086,32 @@ template class Tetrahedra3D10 : public Geometry return true; } - - /** - * Private Friends - */ + ///@} + ///@name Private Friends + ///@{ template friend class Tetrahedra3D10; + ///@} + ///@name Un accessible methods + ///@{ - /** - * Un accessible methods - */ - - + ///@} };// Class Tetrahedra3D10 +///@} +///@name Type Definitions +///@{ -/** - * Input and output - */ +///@} +///@name Input and output +///@{ -/** - * input stream function - */ +/// input stream function template inline std::istream& operator >> ( std::istream& rIStream, Tetrahedra3D10& rThis ); -/** - * output stream function - */ +/// output stream function template inline std::ostream& operator << ( std::ostream& rOStream, const Tetrahedra3D10& rThis ) { @@ -1160,7 +1122,6 @@ template inline std::ostream& operator << ( return rOStream; } - template const GeometryData Tetrahedra3D10::msGeometryData( &msGeometryDimension, @@ -1173,4 +1134,6 @@ GeometryData Tetrahedra3D10::msGeometryData( template const GeometryDimension Tetrahedra3D10::msGeometryDimension(3, 3); +///@} + }// namespace Kratos. From 34543ba9b21ea9cb3a6437b1d01058c296ce574e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vicente=20Mataix=20Ferr=C3=A1ndiz?= Date: Mon, 6 Nov 2023 14:17:59 +0100 Subject: [PATCH 4/4] Suggestions --- kratos/geometries/tetrahedra_3d_10.h | 2 +- kratos/geometries/tetrahedra_3d_4.h | 2 +- kratos/utilities/geometry_utilities.h | 5 ++++- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/kratos/geometries/tetrahedra_3d_10.h b/kratos/geometries/tetrahedra_3d_10.h index b133b7f225b3..0d027c95fbe2 100644 --- a/kratos/geometries/tetrahedra_3d_10.h +++ b/kratos/geometries/tetrahedra_3d_10.h @@ -442,7 +442,7 @@ class Tetrahedra3D10 { // Using linear approximation for planar faces if (this->FacesArePlanar()) { - return GeometryUtils::PointLocalCoordinatesTetrahedra3D4N(*this, rResult, rPoint); + return GeometryUtils::PointLocalCoordinatesPlanarFaceTetrahedra(*this, rResult, rPoint); } else { return BaseType::PointLocalCoordinates( rResult, rPoint ); } diff --git a/kratos/geometries/tetrahedra_3d_4.h b/kratos/geometries/tetrahedra_3d_4.h index 0eb1be06ab4b..598fead0f89f 100644 --- a/kratos/geometries/tetrahedra_3d_4.h +++ b/kratos/geometries/tetrahedra_3d_4.h @@ -970,7 +970,7 @@ template class Tetrahedra3D4 : public Geometry const CoordinatesArrayType& rPoint ) const override { - return GeometryUtils::PointLocalCoordinatesTetrahedra3D4N(*this, rResult, rPoint); + return GeometryUtils::PointLocalCoordinatesPlanarFaceTetrahedra(*this, rResult, rPoint); } /** diff --git a/kratos/utilities/geometry_utilities.h b/kratos/utilities/geometry_utilities.h index 7b8a019cf941..bb3389333d8f 100644 --- a/kratos/utilities/geometry_utilities.h +++ b/kratos/utilities/geometry_utilities.h @@ -64,12 +64,15 @@ class KRATOS_API(KRATOS_CORE) GeometryUtils * @return The vector containing the local coordinates of the point */ template - static inline typename TGeometryType::CoordinatesArrayType& PointLocalCoordinatesTetrahedra3D4N( + static inline typename TGeometryType::CoordinatesArrayType& PointLocalCoordinatesPlanarFaceTetrahedra( const TGeometryType& rGeometry, typename TGeometryType::CoordinatesArrayType& rResult, const typename TGeometryType::CoordinatesArrayType& rPoint ) { + // Debug check that it is at least a tetrahedra + KRATOS_DEBUG_ERROR_IF_NOT(rGeometry.GetGeometryFamily() == GeometryData::KratosGeometryFamily::Kratos_Tetrahedra) << "Geometry should be a tetrahedra in order to use PointLocalCoordinatesPlanarFaceTetrahedra" << std::endl; + // Compute RHS array_1d X; X[0] = 1.0;