KratosMultiphysics
KRATOS Multiphysics (Kratos) is a framework for building parallel, multi-disciplinary simulation software, aiming at modularity, extensibility, and high performance. Kratos is written in C++, and counts with an extensive Python interface.
element_utilities.hpp
Go to the documentation of this file.
1 //
2 // Project Name: KratosPoromechanicsApplication $
3 // Created by: $Author: JMCarbonell $
4 // Last modified by: $Co-Author: $
5 // Date: $Date: July 2018 $
6 // Revision: $Revision: 0.0 $
7 //
8 //
9 
10 #if !defined(KRATOS_ELEMENT_UTILITIES_H_INCLUDED )
11 #define KRATOS_ELEMENT_UTILITIES_H_INCLUDED
12 
13 // System includes
14 
15 // External includes
16 
17 // Project includes
18 #include "custom_utilities/solid_mechanics_math_utilities.hpp"
19 #include "poromechanics_application_variables.h"
20 
21 namespace Kratos
22 {
23 
24 class ElementUtilities
25 {
26  public:
27 
30 
31  //definition of the size type
32  typedef std::size_t SizeType;
33 
35  typedef Node NodeType;
36 
38  typedef Geometry<NodeType> GeometryType;
40 
41 
47  static inline void CalculateDeltaPosition(Matrix& rDeltaPosition, const GeometryType& rGeometry)
48  {
49 
50  const SizeType number_of_nodes = rGeometry.PointsNumber();
51  const SizeType dimension = rGeometry.WorkingSpaceDimension();
52 
53  if(rDeltaPosition.size1() != number_of_nodes || rDeltaPosition.size2() != dimension)
54  rDeltaPosition.resize(number_of_nodes,dimension,false);
55 
56  //noalias(rDeltaPosition) = ZeroMatrix(number_of_nodes,dimension);
57 
58  if( rGeometry[0].SolutionStepsDataHas(STEP_DISPLACEMENT) )
59  {
60  for ( SizeType i = 0; i < number_of_nodes; i++ )
61  {
62  const array_1d<double, 3 > & CurrentStepDisplacement = rGeometry[i].FastGetSolutionStepValue(STEP_DISPLACEMENT,0);
63 
64  for ( SizeType j = 0; j < dimension; j++ )
65  {
66  rDeltaPosition(i,j) = CurrentStepDisplacement[j];
67  }
68 
69  }
70  }
71  else{
72 
73  for ( SizeType i = 0; i < number_of_nodes; i++ )
74  {
75  const array_1d<double, 3 > & CurrentDisplacement = rGeometry[i].FastGetSolutionStepValue(DISPLACEMENT);
76  const array_1d<double, 3 > & PreviousDisplacement = rGeometry[i].FastGetSolutionStepValue(DISPLACEMENT,1);
77 
78  for ( SizeType j = 0; j < dimension; j++ )
79  {
80  rDeltaPosition(i,j) = CurrentDisplacement[j]-PreviousDisplacement[j];
81  }
82  }
83 
84  }
85 
86  }
87 
93  static inline void CalculateTotalDeltaPosition(Matrix & rDeltaPosition, const GeometryType& rGeometry)
94  {
95  const SizeType number_of_nodes = rGeometry.PointsNumber();
96  const SizeType dimension = rGeometry.WorkingSpaceDimension();
97 
98  if(rDeltaPosition.size1() != number_of_nodes || rDeltaPosition.size2() != dimension)
99  rDeltaPosition.resize(number_of_nodes,dimension,false);
100 
101  //noalias(rDeltaPosition) = ZeroMatrix(number_of_nodes,dimension);
102 
103  for ( SizeType i = 0; i < number_of_nodes; i++ )
104  {
105  const array_1d<double, 3 > & CurrentDisplacement = rGeometry[i].FastGetSolutionStepValue(DISPLACEMENT);
106 
107  for ( SizeType j = 0; j < dimension; j++ )
108  {
109  rDeltaPosition(i,j) = CurrentDisplacement[j];
110  }
111  }
112 
113  }
114 
122  static inline void CalculateVelocityGradient(Matrix& rVelocityGradient, const GeometryType& rGeometry,
123  const Matrix& rDN_DX, const double Alpha = 1.0)
124  {
125 
126  const SizeType number_of_nodes = rGeometry.PointsNumber();
127  const SizeType dimension = rGeometry.WorkingSpaceDimension();
128 
129  if( rVelocityGradient.size1() != dimension || rVelocityGradient.size2() != dimension )
130  rVelocityGradient.resize(dimension,dimension);
131 
132  noalias(rVelocityGradient) = ZeroMatrix(dimension,dimension);
133 
134  if( Alpha != 1.0 ){
135 
136  if( dimension == 2 )
137  {
138  for ( SizeType i = 0; i < number_of_nodes; i++ )
139  {
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 );
146  }
147  }
148  else if( dimension == 3)
149  {
150  for ( SizeType i = 0; i < number_of_nodes; i++ )
151  {
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 );
163  }
164  }
165 
166  }
167  else{
168 
169  if( dimension == 2 )
170  {
171  for ( SizeType i = 0; i < number_of_nodes; i++ )
172  {
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 );
178  }
179 
180  }
181  else if( dimension == 3)
182  {
183  for ( SizeType i = 0; i < number_of_nodes; i++ )
184  {
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 );
195  }
196  }
197  }
198 
199  }
200 
208  static inline void CalculateDeformationGradient(Matrix& rDeformationGradient, const GeometryType& rGeometry,
209  const Matrix& rDN_DX, const Matrix& rDeltaPosition)
210  {
211  const SizeType number_of_nodes = rGeometry.PointsNumber();
212  const SizeType dimension = rGeometry.WorkingSpaceDimension();
213 
214  if( rDeformationGradient.size1() != dimension || rDeformationGradient.size2() != dimension )
215  rDeformationGradient.resize(dimension, dimension, false);
216 
217  noalias(rDeformationGradient) = IdentityMatrix(dimension);
218 
219  if( dimension == 2 )
220  {
221  for ( SizeType i = 0; i < number_of_nodes; i++ )
222  {
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 );
227  }
228 
229  }
230  else if( dimension == 3)
231  {
232  for ( SizeType i = 0; i < number_of_nodes; i++ )
233  {
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 );
243  }
244  }
245  else
246  {
247  KRATOS_ERROR << " something is wrong with the dimension when computing the Deformation Gradient " << std::endl;
248  }
249 
250  }
251 
259  static inline void CalculateVelocityGradientVector(Vector& rVelocityGradient, const GeometryType& rGeometry,
260  const Matrix& rDN_DX, const double Alpha = 1.0)
261  {
262 
263  const SizeType number_of_nodes = rGeometry.PointsNumber();
264  const SizeType dimension = rGeometry.WorkingSpaceDimension();
265 
266  if( rVelocityGradient.size() != dimension*dimension )
267  rVelocityGradient.resize(dimension*dimension);
268 
269  noalias(rVelocityGradient) = ZeroVector(dimension * dimension);
270 
271  array_1d<double,3> Velocity;
272  if( dimension == 2 )
273  {
274  for ( SizeType i = 0; i < number_of_nodes; i++ )
275  {
276  const array_1d<double,3>& rPreviousVelocity = rGeometry[i].FastGetSolutionStepValue(VELOCITY,1);
277  const array_1d<double,3>& rCurrentVelocity = rGeometry[i].FastGetSolutionStepValue(VELOCITY);
278 
279  Velocity = rCurrentVelocity * Alpha + rPreviousVelocity * (1.0-Alpha);
280 
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 );
285  }
286  }
287  else if( dimension == 3)
288  {
289 
290  for ( SizeType i = 0; i < number_of_nodes; i++ )
291  {
292  const array_1d<double,3>& rPreviousVelocity = rGeometry[i].FastGetSolutionStepValue(VELOCITY,1);
293  const array_1d<double,3>& rCurrentVelocity = rGeometry[i].FastGetSolutionStepValue(VELOCITY);
294 
295  Velocity = rCurrentVelocity * Alpha + rPreviousVelocity * (1.0-Alpha);
296 
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 );
300 
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 );
304 
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 );
308  }
309 
310  }
311  else
312  {
313  KRATOS_ERROR << " something is wrong with the dimension when computing symmetric velocity gradient " << std::endl;
314  }
315 
316 
317  }
318 
324  static inline void CalculateSymmetricVelocityGradientVector(const Matrix& rVelocityGradientMatrix,
325  Vector& rSymmetricVelocityGradientVector,
326  const SizeType& rDimension)
327  {
328 
329  if( rDimension == 2 )
330  {
331  if ( rSymmetricVelocityGradientVector.size() != 3 ) rSymmetricVelocityGradientVector.resize( 3, false );
332 
333  rSymmetricVelocityGradientVector[0] = rVelocityGradientMatrix( 0, 0 );
334  rSymmetricVelocityGradientVector[1] = rVelocityGradientMatrix( 1, 1 );
335  rSymmetricVelocityGradientVector[2] = 0.5 * (rVelocityGradientMatrix( 0, 1 ) + rVelocityGradientMatrix( 1, 0 )); // xy
336 
337  }
338  else if( rDimension == 3 )
339  {
340  if ( rSymmetricVelocityGradientVector.size() != 6 ) rSymmetricVelocityGradientVector.resize( 6, false );
341 
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 ) ); // xy
346  rSymmetricVelocityGradientVector[4] = 0.5 * ( rVelocityGradientMatrix( 1, 2 ) + rVelocityGradientMatrix( 2, 1 ) ); // yz
347  rSymmetricVelocityGradientVector[5] = 0.5 * ( rVelocityGradientMatrix( 0, 2 ) + rVelocityGradientMatrix( 2, 0 ) ); // xz
348 
349  }
350  else
351  {
352  KRATOS_ERROR << " something is wrong with the dimension symmetric velocity gradient " << std::endl;
353  }
354 
355  }
356 
362  static inline void CalculateSkewSymmetricVelocityGradientVector(const Matrix& rVelocityGradientMatrix,
363  Vector& rSkewSymmetricVelocityGradientVector,
364  const SizeType& rDimension)
365  {
366 
367  if( rDimension == 2 )
368  {
369  if ( rSkewSymmetricVelocityGradientVector.size() != 3 ) rSkewSymmetricVelocityGradientVector.resize( 3, false );
370 
371  rSkewSymmetricVelocityGradientVector[0] = 0.0;
372  rSkewSymmetricVelocityGradientVector[1] = 0.0;
373  rSkewSymmetricVelocityGradientVector[2] = 0.5 * (rVelocityGradientMatrix( 0, 1 ) - rVelocityGradientMatrix( 1, 0 )); // xy
374 
375  }
376  else if( rDimension == 3 )
377  {
378  if ( rSkewSymmetricVelocityGradientVector.size() != 6 ) rSkewSymmetricVelocityGradientVector.resize( 6, false );
379 
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 ) ); // xy
384  rSkewSymmetricVelocityGradientVector[4] = 0.5 * ( rVelocityGradientMatrix( 1, 2 ) - rVelocityGradientMatrix( 2, 1 ) ); // yz
385  rSkewSymmetricVelocityGradientVector[5] = 0.5 * ( rVelocityGradientMatrix( 0, 2 ) - rVelocityGradientMatrix( 2, 0 ) ); // xz
386 
387  }
388  else
389  {
390  KRATOS_ERROR << " something is wrong with the dimension skew symmetric velocity gradient " << std::endl;
391  }
392 
393  }
394 
395 
402  static inline void CalculateLinearDeformationMatrix(Matrix& rDeformationMatrix, const GeometryType& rGeometry, const Matrix& rDN_DX)
403  {
404 
405  const SizeType number_of_nodes = rGeometry.PointsNumber();
406  const SizeType dimension = rGeometry.WorkingSpaceDimension();
407  unsigned int voigt_size = dimension * (dimension +1) * 0.5;
408 
409  if ( rDeformationMatrix.size1() != voigt_size || rDeformationMatrix.size2() != dimension*number_of_nodes )
410  rDeformationMatrix.resize(voigt_size, dimension*number_of_nodes, false );
411 
412 
413  if( dimension == 2 )
414  {
415 
416  for ( SizeType i = 0; i < number_of_nodes; i++ )
417  {
418  unsigned int index = 2 * i;
419 
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 );
426  }
427 
428  }
429  else if( dimension == 3 )
430  {
431  for ( SizeType i = 0; i < number_of_nodes; i++ )
432  {
433  unsigned int index = 3 * i;
434 
435  rDeformationMatrix( 0, index + 0 ) = rDN_DX( i, 0 );
436  rDeformationMatrix( 0, index + 1 ) = 0.0;
437  rDeformationMatrix( 0, index + 2 ) = 0.0;
438 
439  rDeformationMatrix( 1, index + 0 ) = 0.0;
440  rDeformationMatrix( 1, index + 1 ) = rDN_DX( i, 1 );
441  rDeformationMatrix( 1, index + 2 ) = 0.0;
442 
443  rDeformationMatrix( 2, index + 0 ) = 0.0;
444  rDeformationMatrix( 2, index + 1 ) = 0.0;
445  rDeformationMatrix( 2, index + 2 ) = rDN_DX( i, 2 );
446 
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;
450 
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 );
454 
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 );
458 
459  }
460  }
461  else
462  {
463  KRATOS_ERROR << " something is wrong with the dimension when computing linear DeformationMatrix " << std::endl;
464  }
465 
466  }
467 
468 
474  static inline double CalculateStressNorm(const Vector& rStressVector)
475  {
476  Matrix LocalStressTensor = MathUtils<double>::StressVectorToTensor(rStressVector); //reduced dimension stress tensor
477 
478  Matrix StressTensor(3,3); //3D stress tensor
479  noalias(StressTensor) = ZeroMatrix(3,3);
480  for(unsigned int i=0; i<LocalStressTensor.size1(); i++)
481  {
482  for(unsigned int j=0; j<LocalStressTensor.size2(); j++)
483  {
484  StressTensor(i,j) = LocalStressTensor(i,j);
485  }
486  }
487 
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)));
491 
492  StressNorm = sqrt(StressNorm);
493 
494  return StressNorm;
495 
496  };
497 
498 
504  static inline double CalculateVonMises(const Vector& rStressVector)
505  {
506  Matrix LocalStressTensor = MathUtils<double>::StressVectorToTensor(rStressVector); //reduced dimension stress tensor
507 
508  Matrix StressTensor(3,3); //3D stress tensor
509  noalias(StressTensor) = ZeroMatrix(3,3);
510  for(unsigned int i=0; i<LocalStressTensor.size1(); i++)
511  {
512  for(unsigned int j=0; j<LocalStressTensor.size2(); j++)
513  {
514  StressTensor(i,j) = LocalStressTensor(i,j);
515  }
516  }
517 
518  //in general coordinates:
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)));
523 
524  if( SigmaEquivalent < 0 )
525  SigmaEquivalent = 0;
526 
527  SigmaEquivalent = sqrt(SigmaEquivalent);
528 
529  return SigmaEquivalent;
530  }
531 
537  static inline double CalculateVonMisesUsingPrincipalStresses(const Vector& rStressVector)
538  {
539 
540  //in principal stresses:
541 
542  Matrix LocalStressTensor = MathUtils<double>::StressVectorToTensor(rStressVector); //reduced dimension stress tensor
543 
544  Matrix StressTensor(3,3); //3D stress tensor
545  noalias(StressTensor) = ZeroMatrix(3,3);
546  for(unsigned int i=0; i<LocalStressTensor.size1(); i++)
547  {
548  for(unsigned int j=0; j<LocalStressTensor.size2(); j++)
549  {
550  StressTensor(i,j) = LocalStressTensor(i,j);
551  }
552  }
553 
554 
555  double tolerance = 1e-10;
556  double zero = 1e-10;
557  double NormStress = 0.00;
558 
559  CheckZeroDiagonalComponents (StressTensor);
560 
561  Vector PrincipalStress(3);
562  noalias(PrincipalStress) = ZeroVector(3);
563 
564  NormStress =SolidMechanicsMathUtilities<double>::NormTensor(StressTensor);
565 
566  Vector MainStresses(3);
567  noalias(MainStresses) = ZeroVector(3);
568 
569  bool main_tensor = CheckPrincipalStresses( StressTensor );
570 
571  if(!main_tensor)
572  {
573 
574  if(NormStress>1e-6)
575  {
576  MainStresses = SolidMechanicsMathUtilities<double>::EigenValues(StressTensor,tolerance,zero);
577  }
578  else
579  {
580  noalias(MainStresses) = ZeroVector(3);
581  }
582 
583  }
584  else
585  {
586  noalias(MainStresses) = ZeroVector(3);
587  for(unsigned int i=0; i<StressTensor.size1(); i++)
588  MainStresses[i]=StressTensor(i,i);
589  }
590 
591 
592  for(unsigned int i=0; i<MainStresses.size(); i++)
593  PrincipalStress[i]=MainStresses[i];
594 
595 
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]));
599 
600 
601  SigmaEquivalent = sqrt(SigmaEquivalent);
602 
603  return SigmaEquivalent;
604  }
605 
606 
607 
608 
609 
610  protected:
611 
616  static inline void CheckZeroDiagonalComponents (Matrix& StressTensor)
617  {
618  // No null diagonal terms are accepted in the eigenvalue calculation
619  for(unsigned int i=0; i<StressTensor.size1(); i++)
620  {
621  if (fabs(StressTensor(i,i))<1e-10)
622  {
623  StressTensor(i,i) = 1e-10;
624  }
625  }
626  }
627 
633  static inline bool CheckPrincipalStresses (Matrix& StressTensor)
634  {
635  // No null diagonal terms are accepted in the eigenvalue calculation
636  bool main = true;
637  for(unsigned int i=0; i<StressTensor.size1(); i++)
638  {
639  for(unsigned int j=0; j<StressTensor.size2(); j++)
640  {
641  if(i!=j)
642  {
643  if (fabs(StressTensor(i,j))>1e-10)
644  {
645  main = false;
646  }
647  }
648  }
649  }
650 
651  return main;
652  }
653 
654 };
655 } // namespace Kratos.
656 
657 #endif // KRATOS_ELELMENT_UTILITIES_H_INCLUDED
Kratos::ElementUtilities::CalculateDeltaPosition
static void CalculateDeltaPosition(Matrix &rDeltaPosition, const GeometryType &rGeometry)
Calculate Delta Position.
Definition: element_utilities.hpp:47
Kratos::ElementUtilities::CalculateVonMises
static double CalculateVonMises(const Vector &rStressVector)
Calculate VonMises stress.
Definition: element_utilities.hpp:504
Kratos::ElementUtilities::CalculateVelocityGradientVector
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
Kratos::ElementUtilities::CalculateLinearDeformationMatrix
static void CalculateLinearDeformationMatrix(Matrix &rDeformationMatrix, const GeometryType &rGeometry, const Matrix &rDN_DX)
Calculate Linear deformation matrix BL.
Definition: element_utilities.hpp:402
Kratos::ElementUtilities::NodeType
Node NodeType
definition of node type (default is: Node)
Definition: element_utilities.hpp:35
Kratos::ElementUtilities::CheckZeroDiagonalComponents
static void CheckZeroDiagonalComponents(Matrix &StressTensor)
Check and correct diagonal terms in the stress tensor.
Definition: element_utilities.hpp:616
Kratos::ElementUtilities::CalculateDeformationGradient
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
Kratos::ElementUtilities::CalculateVonMisesUsingPrincipalStresses
static double CalculateVonMisesUsingPrincipalStresses(const Vector &rStressVector)
Calculate VonMises stress.
Definition: element_utilities.hpp:537
Kratos::ElementUtilities::CalculateSymmetricVelocityGradientVector
static void CalculateSymmetricVelocityGradientVector(const Matrix &rVelocityGradientMatrix, Vector &rSymmetricVelocityGradientVector, const SizeType &rDimension)
Calculate the symmetric VelocityGradient vector.
Definition: element_utilities.hpp:324
Kratos::ElementUtilities::CalculateVelocityGradient
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
Kratos::ElementUtilities::CheckPrincipalStresses
static bool CheckPrincipalStresses(Matrix &StressTensor)
Check no zero diagonal terms in the diagonalized stress tensor.
Definition: element_utilities.hpp:633
Kratos::ElementUtilities::CalculateSkewSymmetricVelocityGradientVector
static void CalculateSkewSymmetricVelocityGradientVector(const Matrix &rVelocityGradientMatrix, Vector &rSkewSymmetricVelocityGradientVector, const SizeType &rDimension)
Calculate the skew-symmetric VelocityGradient vector.
Definition: element_utilities.hpp:362
Kratos::ElementUtilities::GeometryType
Geometry< NodeType > GeometryType
definition of the geometry type with given NodeType
Definition: element_utilities.hpp:38
Kratos::ElementUtilities::SizeType
std::size_t SizeType
Definition: element_utilities.hpp:32
Kratos::ElementUtilities::CalculateStressNorm
static double CalculateStressNorm(const Vector &rStressVector)
Calculate Norm of stresses.
Definition: element_utilities.hpp:474
Kratos::ElementUtilities::CalculateTotalDeltaPosition
static void CalculateTotalDeltaPosition(Matrix &rDeltaPosition, const GeometryType &rGeometry)
Calculate Total Delta Position.
Definition: element_utilities.hpp:93
Kratos::Geometry
Geometry base class.
Definition: geometry.h:71
Kratos::Geometry::PointsNumber
SizeType PointsNumber() const
Definition: geometry.h:528
Kratos::Geometry::WorkingSpaceDimension
SizeType WorkingSpaceDimension() const
Definition: geometry.h:1287
Kratos::Internals::Matrix::resize
void resize(std::size_t NewSize1, std::size_t NewSize2, bool preserve=0)
Definition: amatrix_interface.h:224
Kratos::MathUtils::StressVectorToTensor
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
Kratos::Node
This class defines the node.
Definition: node.h:65
Kratos::SolidMechanicsMathUtilities::EigenValues
static Vector EigenValues(const Matrix &A, double tolerance, double zero)
Definition: solid_mechanics_math_utilities.hpp:142
Kratos::SolidMechanicsMathUtilities::NormTensor
static double NormTensor(Matrix &Tensor)
Definition: solid_mechanics_math_utilities.hpp:1179
Kratos::array_1d< double, 3 >
KRATOS_ERROR
#define KRATOS_ERROR
Definition: exception.h:161
main
int main()
Definition: mesh_to_clu_converter.cpp:12
ElementUtilities
This class includes several utilities necessaries for the computation of the different elements.
Kratos
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
Kratos::ZeroVector
KratosZeroVector< double > ZeroVector
Definition: amatrix_interface.h:561
Kratos::IdentityMatrix
AMatrix::IdentityMatrix< double > IdentityMatrix
Definition: amatrix_interface.h:564
Kratos::ZeroMatrix
KratosZeroMatrix< double > ZeroMatrix
Definition: amatrix_interface.h:559
Kratos::noalias
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
isotropic_damage_automatic_differentiation.dimension
int dimension
Definition: isotropic_damage_automatic_differentiation.py:123
quadrature.Alpha
def Alpha(n, j)
Definition: quadrature.py:93
quadrature.j
int j
Definition: quadrature.py:648
tensors::i
integer i
Definition: TensorModule.f:17
test_pureconvectionsolver_benchmarking.zero
zero
Definition: test_pureconvectionsolver_benchmarking.py:94
testing.run_cpp_mpi_tests.e
e
Definition: run_cpp_mpi_tests.py:31
poromechanics_application_variables.h
Generated on Sat Feb 3 2024 02:01:53 for KratosMultiphysics by doxygen 1.9.1