Skip to content
Merged
Show file tree
Hide file tree
Changes from all 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
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.

I would propose to remove all unused functions from this file after you have merged this PR.

Original file line number Diff line number Diff line change
Expand Up @@ -39,23 +39,11 @@ class GeoMechanicsMathUtilities
* @name type definitions
* @{
*/
typedef Matrix MatrixType;

typedef Vector VectorType;

typedef unsigned int IndexType;

typedef unsigned int SizeType;

typedef MathUtils<TDataType> MathUtilsType;

typedef DenseVector<Vector> Second_Order_Tensor;

typedef DenseVector<Second_Order_Tensor> Third_Order_Tensor;

typedef DenseVector<DenseVector<Matrix>> Fourth_Order_Tensor;

typedef matrix<Second_Order_Tensor> Matrix_Second_Tensor;
using MatrixType = Matrix;
using VectorType = Vector;
using IndexType = unsigned int;
using SizeType = unsigned int;
using Fourth_Order_Tensor = DenseVector<DenseVector<Matrix>>;
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.

Perhaps for consistency, we should rename this alias to FourthOrderTensorType?


/**
* @}
Expand Down Expand Up @@ -105,7 +93,7 @@ class GeoMechanicsMathUtilities
}

if (p < 0)
return false; // in this case the square roots below will be negative. This substitutes with better efficiency lines 107-110
Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

This one I changed since due to the clang-formatting the lines shifted

return false; // in this case the square roots below will be negative. This substitutes with better efficiency lines 100-102

solution(0) = -sqrt(-4.0 / 3.0 * p) *
cos(1.0 / 3.0 * acos(-q / 2.0 * sqrt(-27.0 / (p * p * p))) + Globals::Pi / 3.0) -
Expand All @@ -116,11 +104,6 @@ class GeoMechanicsMathUtilities
cos(1.0 / 3.0 * acos(-q / 2.0 * sqrt(-27.0 / (p * p * p))) - Globals::Pi / 3.0) -
b / (3 * a);

// if(std::isnan<double>(solution(0)) || std::isnan<double>(solution(1))|| std::isnan<double>(solution(2)))
// {
// return false;
// }

return true;
}

Expand Down Expand Up @@ -260,7 +243,6 @@ class GeoMechanicsMathUtilities
B(1, 2) = inv_p * A(1, 2);
B(2, 1) = inv_p * A(2, 1);

// r = det(B) / 2
double r = 0.5 * (B(0, 0) * B(1, 1) * B(2, 2) + B(0, 1) * B(1, 2) * B(2, 0) +
B(1, 0) * B(2, 1) * B(0, 2) - B(2, 0) * B(1, 1) * B(0, 2) -
B(1, 0) * B(0, 1) * B(2, 2) - B(0, 0) * B(2, 1) * B(1, 2));
Expand Down Expand Up @@ -506,9 +488,6 @@ class GeoMechanicsMathUtilities
Rotation(index1, index1) = c;
Rotation(index2, index2) = c;

// Help.resize(A.size1(),A.size1(),false);
// noalias(Help)=ZeroMatrix(A.size1(),A.size1());

noalias(VDummy) = ZeroMatrix(Help.size1(), Help.size2());

for (unsigned int i = 0; i < Help.size1(); i++) {
Expand All @@ -519,30 +498,6 @@ class GeoMechanicsMathUtilities
}
}
V = VDummy;

// Matrix VTA(3,3);
// noalias(VTA) = ZeroMatrix(3,3);
// for(int i=0; i< Help.size1(); i++)
// {
// for(int j=0; j< Help.size1(); j++)
// {
// for(int k=0; k< Help.size1(); k++)
// {
// VTA(i,j) += V(k,i)*A(k,j);
// }
// }
// }
//
// for(int i=0; i< Help.size1(); i++)
// {
// for(int j=0; j< Help.size1(); j++)
// {
// for(int k=0; k< Help.size1(); k++)
// {
// Help(i,j) += VTA(i,k)*V(k,j);
// }
// }
// }
}

if (!(is_converged)) {
Expand Down Expand Up @@ -711,11 +666,11 @@ class GeoMechanicsMathUtilities

static inline void TensorToMatrix(Fourth_Order_Tensor& Tensor, Matrix& Matrix)
{
// Simetrias seguras
// Cijkl = Cjilk;
// Cijkl = Cklji;
// Symmetric fourth order tensor:
// Cijkl = Cjilk
// Cijkl = Cklji
if (Tensor[0].size() == 3) {
// Tensor de cuarto orden cuyos componentes correspondes a una matriz de 3x3
// Fourth-order tensor whose components correspond to a 3x3 matrix
if (Matrix.size1() != 6 || Matrix.size2() != 6) Matrix.resize(6, 6, false);
Matrix(0, 0) = Tensor[0][0](0, 0);
Matrix(0, 1) = Tensor[0][0](1, 1);
Expand Down Expand Up @@ -759,7 +714,7 @@ class GeoMechanicsMathUtilities
Matrix(5, 4) = Tensor[1][2](0, 2);
Matrix(5, 5) = Tensor[1][2](1, 2);
} else {
// Tensor de cuarto orden cuyos componentes correspondes a una matriz de 2x2
// Fourth-order tensor whose components correspond to a 2x2 matrix
if (Matrix.size1() != 3 || Matrix.size2() != 3) Matrix.resize(3, 3, false);
Matrix(0, 0) = Tensor[0][0](0, 0);
Matrix(0, 1) = Tensor[0][0](1, 1);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -28,16 +28,16 @@ namespace Kratos
/**
* @class GeoStaticCondensationUtility
*
* @brief This utilitiy condenses given degrees of freedom from any element stiffness matrix to model e.g. hinges
* @brief This utility condenses given degrees of freedom from any element stiffness matrix to model e.g. hinges
*
*/

class GeoStaticCondensationUtility
Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

It seems this was copied from structural mechanics, we might just be able to re-use that one.

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.

Yes, please :-)

{
public:
typedef Element ElementType;
typedef std::size_t SizeType;
typedef Matrix MatrixType;
using ElementType = Element;
using SizeType = std::size_t;
using MatrixType = Matrix;

/**
* @brief This function is the main operation of this utility. It sorts the reference matrix w.r.t. the given dofs and condenses the reference matrix by using the following inputs:
Expand All @@ -64,7 +64,7 @@ class GeoStaticCondensationUtility
KRATOS_ERROR_IF(std::abs(detK22) < numerical_limit)
<< "Element " << rTheElement.Id() << " is singular !" << std::endl;

// 2.) K_cond = K11 - K12*inv(K22)*K21
// 2.) K_cond -> K11 - K12*inv(K22)*K21
K_temp = prod(K_temp, sub_matrices[2]);
K_temp = prod(sub_matrices[1], K_temp);
K_temp = sub_matrices[0] - K_temp;
Expand Down Expand Up @@ -259,7 +259,6 @@ class GeoStaticCondensationUtility
}
}
}
// rTheElement.GlobalizeVector(rValues); // globalize local displacements -> global lvl
KRATOS_CATCH("")
}

Expand Down