10 #if !defined(KRATOS_ELEMENT_UTILITIES_H_INCLUDED )
11 #define KRATOS_ELEMENT_UTILITIES_H_INCLUDED
18 #include "custom_utilities/solid_mechanics_math_utilities.hpp"
53 if(rDeltaPosition.size1() != number_of_nodes || rDeltaPosition.size2() !=
dimension)
58 if( rGeometry[0].SolutionStepsDataHas(STEP_DISPLACEMENT) )
62 const array_1d<double, 3 > & CurrentStepDisplacement = rGeometry[
i].FastGetSolutionStepValue(STEP_DISPLACEMENT,0);
66 rDeltaPosition(
i,
j) = CurrentStepDisplacement[
j];
75 const array_1d<double, 3 > & CurrentDisplacement = rGeometry[
i].FastGetSolutionStepValue(DISPLACEMENT);
76 const array_1d<double, 3 > & PreviousDisplacement = rGeometry[
i].FastGetSolutionStepValue(DISPLACEMENT,1);
80 rDeltaPosition(
i,
j) = CurrentDisplacement[
j]-PreviousDisplacement[
j];
98 if(rDeltaPosition.size1() != number_of_nodes || rDeltaPosition.size2() !=
dimension)
105 const array_1d<double, 3 > & CurrentDisplacement = rGeometry[
i].FastGetSolutionStepValue(DISPLACEMENT);
109 rDeltaPosition(
i,
j) = CurrentDisplacement[
j];
129 if( rVelocityGradient.size1() !=
dimension || rVelocityGradient.size2() !=
dimension )
140 const array_1d<double,3>& rPreviousVelocity = rGeometry[
i].FastGetSolutionStepValue(VELOCITY,1);
141 const array_1d<double,3>& rCurrentVelocity = rGeometry[
i].FastGetSolutionStepValue(VELOCITY);
142 rVelocityGradient ( 0 , 0 ) += (rCurrentVelocity[0] *
Alpha + rPreviousVelocity[0]* (1.0-
Alpha))*rDN_DX (
i , 0 );
143 rVelocityGradient ( 0 , 1 ) += (rCurrentVelocity[0] *
Alpha + rPreviousVelocity[0]* (1.0-
Alpha))*rDN_DX (
i , 1 );
144 rVelocityGradient ( 1 , 0 ) += (rCurrentVelocity[1] *
Alpha + rPreviousVelocity[1]* (1.0-
Alpha))*rDN_DX (
i , 0 );
145 rVelocityGradient ( 1 , 1 ) += (rCurrentVelocity[1] *
Alpha + rPreviousVelocity[1]* (1.0-
Alpha))*rDN_DX (
i , 1 );
152 const array_1d<double,3>& rPreviousVelocity = rGeometry[
i].FastGetSolutionStepValue(VELOCITY,1);
153 const array_1d<double,3>& rCurrentVelocity = rGeometry[
i].FastGetSolutionStepValue(VELOCITY);
154 rVelocityGradient ( 0 , 0 ) += (rCurrentVelocity[0] *
Alpha + rPreviousVelocity[0]* (1.0-
Alpha))*rDN_DX (
i , 0 );
155 rVelocityGradient ( 0 , 1 ) += (rCurrentVelocity[0] *
Alpha + rPreviousVelocity[0]* (1.0-
Alpha))*rDN_DX (
i , 1 );
156 rVelocityGradient ( 0 , 2 ) += (rCurrentVelocity[0] *
Alpha + rPreviousVelocity[0]* (1.0-
Alpha))*rDN_DX (
i , 2 );
157 rVelocityGradient ( 1 , 0 ) += (rCurrentVelocity[1] *
Alpha + rPreviousVelocity[1]* (1.0-
Alpha))*rDN_DX (
i , 0 );
158 rVelocityGradient ( 1 , 1 ) += (rCurrentVelocity[1] *
Alpha + rPreviousVelocity[1]* (1.0-
Alpha))*rDN_DX (
i , 1 );
159 rVelocityGradient ( 1 , 2 ) += (rCurrentVelocity[1] *
Alpha + rPreviousVelocity[1]* (1.0-
Alpha))*rDN_DX (
i , 2 );
160 rVelocityGradient ( 2 , 0 ) += (rCurrentVelocity[2] *
Alpha + rPreviousVelocity[2]* (1.0-
Alpha))*rDN_DX (
i , 0 );
161 rVelocityGradient ( 2 , 1 ) += (rCurrentVelocity[2] *
Alpha + rPreviousVelocity[2]* (1.0-
Alpha))*rDN_DX (
i , 1 );
162 rVelocityGradient ( 2 , 2 ) += (rCurrentVelocity[2] *
Alpha + rPreviousVelocity[2]* (1.0-
Alpha))*rDN_DX (
i , 2 );
173 const array_1d<double,3>& rCurrentVelocity = rGeometry[
i].FastGetSolutionStepValue(VELOCITY);
174 rVelocityGradient ( 0 , 0 ) += rCurrentVelocity[0]*rDN_DX (
i , 0 );
175 rVelocityGradient ( 0 , 1 ) += rCurrentVelocity[0]*rDN_DX (
i , 1 );
176 rVelocityGradient ( 1 , 0 ) += rCurrentVelocity[1]*rDN_DX (
i , 0 );
177 rVelocityGradient ( 1 , 1 ) += rCurrentVelocity[1]*rDN_DX (
i , 1 );
185 const array_1d<double,3>& rCurrentVelocity = rGeometry[
i].FastGetSolutionStepValue(VELOCITY);
186 rVelocityGradient ( 0 , 0 ) += rCurrentVelocity[0]*rDN_DX (
i , 0 );
187 rVelocityGradient ( 0 , 1 ) += rCurrentVelocity[0]*rDN_DX (
i , 1 );
188 rVelocityGradient ( 0 , 2 ) += rCurrentVelocity[0]*rDN_DX (
i , 2 );
189 rVelocityGradient ( 1 , 0 ) += rCurrentVelocity[1]*rDN_DX (
i , 0 );
190 rVelocityGradient ( 1 , 1 ) += rCurrentVelocity[1]*rDN_DX (
i , 1 );
191 rVelocityGradient ( 1 , 2 ) += rCurrentVelocity[1]*rDN_DX (
i , 2 );
192 rVelocityGradient ( 2 , 0 ) += rCurrentVelocity[2]*rDN_DX (
i , 0 );
193 rVelocityGradient ( 2 , 1 ) += rCurrentVelocity[2]*rDN_DX (
i , 1 );
194 rVelocityGradient ( 2 , 2 ) += rCurrentVelocity[2]*rDN_DX (
i , 2 );
214 if( rDeformationGradient.size1() !=
dimension || rDeformationGradient.size2() !=
dimension )
223 rDeformationGradient ( 0 , 0 ) += rDeltaPosition(
i,0)*rDN_DX (
i , 0 );
224 rDeformationGradient ( 0 , 1 ) += rDeltaPosition(
i,0)*rDN_DX (
i , 1 );
225 rDeformationGradient ( 1 , 0 ) += rDeltaPosition(
i,1)*rDN_DX (
i , 0 );
226 rDeformationGradient ( 1 , 1 ) += rDeltaPosition(
i,1)*rDN_DX (
i , 1 );
234 rDeformationGradient ( 0 , 0 ) += rDeltaPosition(
i,0)*rDN_DX (
i , 0 );
235 rDeformationGradient ( 0 , 1 ) += rDeltaPosition(
i,0)*rDN_DX (
i , 1 );
236 rDeformationGradient ( 0 , 2 ) += rDeltaPosition(
i,0)*rDN_DX (
i , 2 );
237 rDeformationGradient ( 1 , 0 ) += rDeltaPosition(
i,1)*rDN_DX (
i , 0 );
238 rDeformationGradient ( 1 , 1 ) += rDeltaPosition(
i,1)*rDN_DX (
i , 1 );
239 rDeformationGradient ( 1 , 2 ) += rDeltaPosition(
i,1)*rDN_DX (
i , 2 );
240 rDeformationGradient ( 2 , 0 ) += rDeltaPosition(
i,2)*rDN_DX (
i , 0 );
241 rDeformationGradient ( 2 , 1 ) += rDeltaPosition(
i,2)*rDN_DX (
i , 1 );
242 rDeformationGradient ( 2 , 2 ) += rDeltaPosition(
i,2)*rDN_DX (
i , 2 );
247 KRATOS_ERROR <<
" something is wrong with the dimension when computing the Deformation Gradient " << std::endl;
276 const array_1d<double,3>& rPreviousVelocity = rGeometry[
i].FastGetSolutionStepValue(VELOCITY,1);
277 const array_1d<double,3>& rCurrentVelocity = rGeometry[
i].FastGetSolutionStepValue(VELOCITY);
279 Velocity = rCurrentVelocity *
Alpha + rPreviousVelocity * (1.0-
Alpha);
281 rVelocityGradient[0] += Velocity[0]*rDN_DX (
i , 0 );
282 rVelocityGradient[1] += Velocity[1]*rDN_DX (
i , 1 );
283 rVelocityGradient[2] += Velocity[0]*rDN_DX (
i , 1 );
284 rVelocityGradient[3] += Velocity[1]*rDN_DX (
i , 0 );
292 const array_1d<double,3>& rPreviousVelocity = rGeometry[
i].FastGetSolutionStepValue(VELOCITY,1);
293 const array_1d<double,3>& rCurrentVelocity = rGeometry[
i].FastGetSolutionStepValue(VELOCITY);
295 Velocity = rCurrentVelocity *
Alpha + rPreviousVelocity * (1.0-
Alpha);
297 rVelocityGradient[0] += Velocity[0]*rDN_DX (
i , 0 );
298 rVelocityGradient[1] += Velocity[1]*rDN_DX (
i , 1 );
299 rVelocityGradient[2] += Velocity[2]*rDN_DX (
i , 2 );
301 rVelocityGradient[3] += Velocity[0]*rDN_DX (
i , 1 );
302 rVelocityGradient[4] += Velocity[1]*rDN_DX (
i , 2 );
303 rVelocityGradient[5] += Velocity[2]*rDN_DX (
i , 0 );
305 rVelocityGradient[6] += Velocity[1]*rDN_DX (
i , 0 );
306 rVelocityGradient[7] += Velocity[2]*rDN_DX (
i , 1 );
307 rVelocityGradient[8] += Velocity[0]*rDN_DX (
i , 2 );
313 KRATOS_ERROR <<
" something is wrong with the dimension when computing symmetric velocity gradient " << std::endl;
325 Vector& rSymmetricVelocityGradientVector,
329 if( rDimension == 2 )
331 if ( rSymmetricVelocityGradientVector.size() != 3 ) rSymmetricVelocityGradientVector.
resize( 3,
false );
333 rSymmetricVelocityGradientVector[0] = rVelocityGradientMatrix( 0, 0 );
334 rSymmetricVelocityGradientVector[1] = rVelocityGradientMatrix( 1, 1 );
335 rSymmetricVelocityGradientVector[2] = 0.5 * (rVelocityGradientMatrix( 0, 1 ) + rVelocityGradientMatrix( 1, 0 ));
338 else if( rDimension == 3 )
340 if ( rSymmetricVelocityGradientVector.size() != 6 ) rSymmetricVelocityGradientVector.
resize( 6,
false );
342 rSymmetricVelocityGradientVector[0] = rVelocityGradientMatrix( 0, 0 );
343 rSymmetricVelocityGradientVector[1] = rVelocityGradientMatrix( 1, 1 );
344 rSymmetricVelocityGradientVector[2] = rVelocityGradientMatrix( 2, 2 );
345 rSymmetricVelocityGradientVector[3] = 0.5 * ( rVelocityGradientMatrix( 0, 1 ) + rVelocityGradientMatrix( 1, 0 ) );
346 rSymmetricVelocityGradientVector[4] = 0.5 * ( rVelocityGradientMatrix( 1, 2 ) + rVelocityGradientMatrix( 2, 1 ) );
347 rSymmetricVelocityGradientVector[5] = 0.5 * ( rVelocityGradientMatrix( 0, 2 ) + rVelocityGradientMatrix( 2, 0 ) );
352 KRATOS_ERROR <<
" something is wrong with the dimension symmetric velocity gradient " << std::endl;
363 Vector& rSkewSymmetricVelocityGradientVector,
367 if( rDimension == 2 )
369 if ( rSkewSymmetricVelocityGradientVector.size() != 3 ) rSkewSymmetricVelocityGradientVector.
resize( 3,
false );
371 rSkewSymmetricVelocityGradientVector[0] = 0.0;
372 rSkewSymmetricVelocityGradientVector[1] = 0.0;
373 rSkewSymmetricVelocityGradientVector[2] = 0.5 * (rVelocityGradientMatrix( 0, 1 ) - rVelocityGradientMatrix( 1, 0 ));
376 else if( rDimension == 3 )
378 if ( rSkewSymmetricVelocityGradientVector.size() != 6 ) rSkewSymmetricVelocityGradientVector.
resize( 6,
false );
380 rSkewSymmetricVelocityGradientVector[0] = 0.0;
381 rSkewSymmetricVelocityGradientVector[1] = 0.0;
382 rSkewSymmetricVelocityGradientVector[2] = 0.0;
383 rSkewSymmetricVelocityGradientVector[3] = 0.5 * ( rVelocityGradientMatrix( 0, 1 ) - rVelocityGradientMatrix( 1, 0 ) );
384 rSkewSymmetricVelocityGradientVector[4] = 0.5 * ( rVelocityGradientMatrix( 1, 2 ) - rVelocityGradientMatrix( 2, 1 ) );
385 rSkewSymmetricVelocityGradientVector[5] = 0.5 * ( rVelocityGradientMatrix( 0, 2 ) - rVelocityGradientMatrix( 2, 0 ) );
390 KRATOS_ERROR <<
" something is wrong with the dimension skew symmetric velocity gradient " << std::endl;
409 if ( rDeformationMatrix.size1() != voigt_size || rDeformationMatrix.size2() !=
dimension*number_of_nodes )
410 rDeformationMatrix.
resize(voigt_size,
dimension*number_of_nodes,
false );
418 unsigned int index = 2 *
i;
420 rDeformationMatrix( 0, index + 0 ) = rDN_DX(
i, 0 );
421 rDeformationMatrix( 0, index + 1 ) = 0.0;
422 rDeformationMatrix( 1, index + 0 ) = 0.0;
423 rDeformationMatrix( 1, index + 1 ) = rDN_DX(
i, 1 );
424 rDeformationMatrix( 2, index + 0 ) = rDN_DX(
i, 1 );
425 rDeformationMatrix( 2, index + 1 ) = rDN_DX(
i, 0 );
433 unsigned int index = 3 *
i;
435 rDeformationMatrix( 0, index + 0 ) = rDN_DX(
i, 0 );
436 rDeformationMatrix( 0, index + 1 ) = 0.0;
437 rDeformationMatrix( 0, index + 2 ) = 0.0;
439 rDeformationMatrix( 1, index + 0 ) = 0.0;
440 rDeformationMatrix( 1, index + 1 ) = rDN_DX(
i, 1 );
441 rDeformationMatrix( 1, index + 2 ) = 0.0;
443 rDeformationMatrix( 2, index + 0 ) = 0.0;
444 rDeformationMatrix( 2, index + 1 ) = 0.0;
445 rDeformationMatrix( 2, index + 2 ) = rDN_DX(
i, 2 );
447 rDeformationMatrix( 3, index + 0 ) = rDN_DX(
i, 1 );
448 rDeformationMatrix( 3, index + 1 ) = rDN_DX(
i, 0 );
449 rDeformationMatrix( 3, index + 2 ) = 0.0;
451 rDeformationMatrix( 4, index + 0 ) = 0.0;
452 rDeformationMatrix( 4, index + 1 ) = rDN_DX(
i, 2 );
453 rDeformationMatrix( 4, index + 2 ) = rDN_DX(
i, 1 );
455 rDeformationMatrix( 5, index + 0 ) = rDN_DX(
i, 2 );
456 rDeformationMatrix( 5, index + 1 ) = 0.0;
457 rDeformationMatrix( 5, index + 2 ) = rDN_DX(
i, 0 );
463 KRATOS_ERROR <<
" something is wrong with the dimension when computing linear DeformationMatrix " << std::endl;
480 for(
unsigned int i=0;
i<LocalStressTensor.size1();
i++)
482 for(
unsigned int j=0;
j<LocalStressTensor.size2();
j++)
484 StressTensor(
i,
j) = LocalStressTensor(
i,
j);
488 double StressNorm = ((StressTensor(0,0)*StressTensor(0,0))+(StressTensor(1,1)*StressTensor(1,1))+(StressTensor(2,2)*StressTensor(2,2))+
489 (StressTensor(0,1)*StressTensor(0,1))+(StressTensor(0,2)*StressTensor(0,2))+(StressTensor(1,2)*StressTensor(1,2))+
490 (StressTensor(1,0)*StressTensor(1,0))+(StressTensor(2,0)*StressTensor(2,0))+(StressTensor(2,1)*StressTensor(2,1)));
492 StressNorm = sqrt(StressNorm);
510 for(
unsigned int i=0;
i<LocalStressTensor.size1();
i++)
512 for(
unsigned int j=0;
j<LocalStressTensor.size2();
j++)
514 StressTensor(
i,
j) = LocalStressTensor(
i,
j);
519 double SigmaEquivalent = (0.5)*((StressTensor(0,0)-StressTensor(1,1))*((StressTensor(0,0)-StressTensor(1,1)))+
520 (StressTensor(1,1)-StressTensor(2,2))*((StressTensor(1,1)-StressTensor(2,2)))+
521 (StressTensor(2,2)-StressTensor(0,0))*((StressTensor(2,2)-StressTensor(0,0)))+
522 6*(StressTensor(0,1)*StressTensor(1,0)+StressTensor(1,2)*StressTensor(2,1)+StressTensor(2,0)*StressTensor(0,2)));
524 if( SigmaEquivalent < 0 )
527 SigmaEquivalent = sqrt(SigmaEquivalent);
529 return SigmaEquivalent;
546 for(
unsigned int i=0;
i<LocalStressTensor.size1();
i++)
548 for(
unsigned int j=0;
j<LocalStressTensor.size2();
j++)
550 StressTensor(
i,
j) = LocalStressTensor(
i,
j);
555 double tolerance = 1
e-10;
557 double NormStress = 0.00;
561 Vector PrincipalStress(3);
587 for(
unsigned int i=0;
i<StressTensor.size1();
i++)
588 MainStresses[
i]=StressTensor(
i,
i);
592 for(
unsigned int i=0;
i<MainStresses.size();
i++)
593 PrincipalStress[
i]=MainStresses[
i];
596 double SigmaEquivalent = (0.5)*((PrincipalStress[0]-PrincipalStress[1])*(PrincipalStress[0]-PrincipalStress[1]) +
597 (PrincipalStress[1]-PrincipalStress[2])*(PrincipalStress[1]-PrincipalStress[2]) +
598 (PrincipalStress[2]-PrincipalStress[0])*(PrincipalStress[2]-PrincipalStress[0]));
601 SigmaEquivalent = sqrt(SigmaEquivalent);
603 return SigmaEquivalent;
619 for(
unsigned int i=0;
i<StressTensor.size1();
i++)
621 if (fabs(StressTensor(
i,
i))<1
e-10)
623 StressTensor(
i,
i) = 1
e-10;
637 for(
unsigned int i=0;
i<StressTensor.size1();
i++)
639 for(
unsigned int j=0;
j<StressTensor.size2();
j++)
643 if (fabs(StressTensor(
i,
j))>1
e-10)
static void CalculateDeltaPosition(Matrix &rDeltaPosition, const GeometryType &rGeometry)
Calculate Delta Position.
Definition: element_utilities.hpp:47
static double CalculateVonMises(const Vector &rStressVector)
Calculate VonMises stress.
Definition: element_utilities.hpp:504
static void CalculateVelocityGradientVector(Vector &rVelocityGradient, const GeometryType &rGeometry, const Matrix &rDN_DX, const double Alpha=1.0)
Calculate the VelocityGradient vector (no voigt form)
Definition: element_utilities.hpp:259
static void CalculateLinearDeformationMatrix(Matrix &rDeformationMatrix, const GeometryType &rGeometry, const Matrix &rDN_DX)
Calculate Linear deformation matrix BL.
Definition: element_utilities.hpp:402
Node NodeType
definition of node type (default is: Node)
Definition: element_utilities.hpp:35
static void CheckZeroDiagonalComponents(Matrix &StressTensor)
Check and correct diagonal terms in the stress tensor.
Definition: element_utilities.hpp:616
static void CalculateDeformationGradient(Matrix &rDeformationGradient, const GeometryType &rGeometry, const Matrix &rDN_DX, const Matrix &rDeltaPosition)
Calculate the Deformation Gradient Tensor.
Definition: element_utilities.hpp:208
static double CalculateVonMisesUsingPrincipalStresses(const Vector &rStressVector)
Calculate VonMises stress.
Definition: element_utilities.hpp:537
static void CalculateSymmetricVelocityGradientVector(const Matrix &rVelocityGradientMatrix, Vector &rSymmetricVelocityGradientVector, const SizeType &rDimension)
Calculate the symmetric VelocityGradient vector.
Definition: element_utilities.hpp:324
static void CalculateVelocityGradient(Matrix &rVelocityGradient, const GeometryType &rGeometry, const Matrix &rDN_DX, const double Alpha=1.0)
Calculate Norm of stresses.VelocityGradient.
Definition: element_utilities.hpp:122
static bool CheckPrincipalStresses(Matrix &StressTensor)
Check no zero diagonal terms in the diagonalized stress tensor.
Definition: element_utilities.hpp:633
static void CalculateSkewSymmetricVelocityGradientVector(const Matrix &rVelocityGradientMatrix, Vector &rSkewSymmetricVelocityGradientVector, const SizeType &rDimension)
Calculate the skew-symmetric VelocityGradient vector.
Definition: element_utilities.hpp:362
Geometry< NodeType > GeometryType
definition of the geometry type with given NodeType
Definition: element_utilities.hpp:38
std::size_t SizeType
Definition: element_utilities.hpp:32
static double CalculateStressNorm(const Vector &rStressVector)
Calculate Norm of stresses.
Definition: element_utilities.hpp:474
static void CalculateTotalDeltaPosition(Matrix &rDeltaPosition, const GeometryType &rGeometry)
Calculate Total Delta Position.
Definition: element_utilities.hpp:93
Geometry base class.
Definition: geometry.h:71
SizeType PointsNumber() const
Definition: geometry.h:528
SizeType WorkingSpaceDimension() const
Definition: geometry.h:1287
void resize(std::size_t NewSize1, std::size_t NewSize2, bool preserve=0)
Definition: amatrix_interface.h:224
static TMatrixType StressVectorToTensor(const TVector &rStressVector)
Transforms a stess vector into a matrix. Stresses are assumed to be stored in the following way: for...
Definition: math_utils.h:1174
This class defines the node.
Definition: node.h:65
static Vector EigenValues(const Matrix &A, double tolerance, double zero)
Definition: solid_mechanics_math_utilities.hpp:142
static double NormTensor(Matrix &Tensor)
Definition: solid_mechanics_math_utilities.hpp:1179
#define KRATOS_ERROR
Definition: exception.h:161
int main()
Definition: mesh_to_clu_converter.cpp:12
This class includes several utilities necessaries for the computation of the different elements.
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
KratosZeroVector< double > ZeroVector
Definition: amatrix_interface.h:561
AMatrix::IdentityMatrix< double > IdentityMatrix
Definition: amatrix_interface.h:564
KratosZeroMatrix< double > ZeroMatrix
Definition: amatrix_interface.h:559
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
int dimension
Definition: isotropic_damage_automatic_differentiation.py:123
def Alpha(n, j)
Definition: quadrature.py:93
int j
Definition: quadrature.py:648
integer i
Definition: TensorModule.f:17
zero
Definition: test_pureconvectionsolver_benchmarking.py:94
e
Definition: run_cpp_mpi_tests.py:31