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.
isochoric_ogden_model.hpp
Go to the documentation of this file.
1 //
2 // Project Name: KratosConstitutiveModelsApplication $
3 // Created by: $Author: JMCarbonell $
4 // Last modified by: $Co-Author: $
5 // Date: $Date: April 2017 $
6 // Revision: $Revision: 0.0 $
7 //
8 //
9 
10 #if !defined(KRATOS_ISOCHORIC_OGDEN_MODEL_H_INCLUDED )
11 #define KRATOS_ISOCHORIC_OGDEN_MODEL_H_INCLUDED
12 
13 // System includes
14 
15 // External includes
16 
17 // Project includes
18 #include "custom_models/elasticity_models/ogden_model.hpp"
19 
20 namespace Kratos
21 {
24 
27 
31 
35 
39 
43 
45 
47  class IsochoricOgdenModel : public OgdenModel
48  {
49  public:
50 
53 
55  KRATOS_CLASS_POINTER_DEFINITION( IsochoricOgdenModel );
56 
60 
62  IsochoricOgdenModel() : OgdenModel() {}
63 
65  IsochoricOgdenModel(IsochoricOgdenModel const& rOther) : OgdenModel(rOther) {}
66 
68  IsochoricOgdenModel& operator=(IsochoricOgdenModel const& rOther)
69  {
70  OgdenModel::operator=(rOther);
71  return *this;
72  }
73 
75  ConstitutiveModel::Pointer Clone() const override
76  {
77  return Kratos::make_shared<IsochoricOgdenModel>(*this);
78  }
79 
81  ~IsochoricOgdenModel() override {}
82 
83 
87 
88 
92 
93 
94  void CalculateStrainEnergy(ModelDataType& rValues, double& rDensityFunction) override
95  {
96  KRATOS_TRY
97 
98  HyperElasticDataType Variables;
99  this->CalculateStrainData(rValues, Variables);
100 
101  rDensityFunction = 0;
102  this->CalculateAndAddIsochoricStrainEnergy( Variables, rDensityFunction );
103  this->CalculateAndAddVolumetricStrainEnergy( Variables, rDensityFunction );
104 
105 
106  KRATOS_CATCH(" ")
107  }
108 
109 
110  void CalculateStressTensor(ModelDataType& rValues, MatrixType& rStressMatrix) override
111  {
112  KRATOS_TRY
113 
114  HyperElasticDataType Variables;
115  this->CalculateStrainData(rValues,Variables);
116 
117  this->CalculateAndAddIsochoricStressTensor(Variables, rStressMatrix);
118 
119  rValues.StressMatrix = rStressMatrix; //store isochoric stress matrix as StressMatrix
120 
121  this->CalculateAndAddVolumetricStressTensor(Variables, rStressMatrix);
122 
123  Variables.State().Set(ConstitutiveModelData::STRESS_COMPUTED);
124 
125  KRATOS_CATCH(" ")
126  }
127 
128 
129  void CalculateConstitutiveTensor(ModelDataType& rValues, Matrix& rConstitutiveMatrix) override
130  {
131  KRATOS_TRY
132 
133  //Initialize ConstitutiveMatrix
134  HyperElasticDataType Variables;
135  this->CalculateStrainData(rValues,Variables);
136 
137  //Calculate Constitutive Matrix
138  this->CalculateAndAddConstitutiveTensor(Variables,rConstitutiveMatrix);
139 
140  KRATOS_CATCH(" ")
141  }
142 
143 
144  void CalculateStressAndConstitutiveTensors(ModelDataType& rValues, MatrixType& rStressMatrix, Matrix& rConstitutiveMatrix) override
145  {
146  KRATOS_TRY
147 
148  HyperElasticDataType Variables;
149  this->CalculateStrainData(rValues,Variables);
150 
151  //Calculate Stress Matrix
152  this->CalculateAndAddIsochoricStressTensor(Variables, rStressMatrix);
153 
154  rValues.StressMatrix = rStressMatrix; //store isochoric stress matrix as StressMatrix
155 
156  this->CalculateAndAddVolumetricStressTensor(Variables, rStressMatrix);
157 
158  //Calculate Constitutive Matrix
159  this->CalculateAndAddConstitutiveTensor(Variables,rConstitutiveMatrix);
160  //this->CalculateAndAddPerturbedConstitutiveTensor(Variables,rConstitutiveMatrix);
161 
162  KRATOS_CATCH(" ")
163  }
164 
168 
169 
173 
174 
178 
180  std::string Info() const override
181  {
182  std::stringstream buffer;
183  buffer << "IsochoricOgdenModel";
184  return buffer.str();
185  }
186 
188  void PrintInfo(std::ostream& rOStream) const override
189  {
190  rOStream << "IsochoricOgdenModel";
191  }
192 
194  void PrintData(std::ostream& rOStream) const override
195  {
196  rOStream << "IsochoricOgdenModel Data";
197  }
198 
202 
203 
205 
206  protected:
207 
210 
211 
215 
216 
220 
221 
225 
226  void CalculateStrainData(ModelDataType& rValues, HyperElasticDataType& rVariables) override
227  {
228  KRATOS_TRY
229 
230  OgdenModel::CalculateStrainData(rValues,rVariables);
231 
232  //Isochoric eigenvalues
233  for(unsigned int i=0; i<3; i++)
234  {
235  rVariables.Strain.Eigen.Values[i] = rVariables.Strain.Eigen.Values[i] / std::pow(rVariables.Strain.Invariants.J, 1.0/3.0);
236  }
237 
238  //Calculate Invariants
239  this->CalculateInvariants(rVariables);
240 
241  //Algorithmic moduli factors
242  this->CalculateScalingFactors(rVariables);
243 
244  //strain check
245  // double D = 0;
246  // MatrixType maxma;
247  // MatrixType MaxMa;
248  // for(unsigned int i=0; i<3; i++)
249  // {
250  // noalias(maxma) = ZeroMatrix(3,3);
251  // noalias(MaxMa) = ZeroMatrix(3,3);
252  // const double& lambda = rVariables.Strain.Eigen.Values[i];
253 
254  // D = 2.0 * lambda*lambda*lambda*lambda - rVariables.Strain.Invariants.I1 * lambda*lambda + rVariables.Strain.Invariants.I3 / (lambda*lambda);
255 
256  // array_1d<double,3> EigenVector;
257  // noalias(EigenVector) = matrix_row<const MatrixType>(rVariables.Strain.Eigen.Vectors,i);
258 
259  // std::cout<<" naxna "<<outer_prod(EigenVector,EigenVector)<<std::endl;
260 
261  // if( D!= 0 ){
262  // noalias(maxma)=(prod(rVariables.Strain.Matrix,rVariables.Strain.Matrix) - (rVariables.Strain.Invariants.I1-rVariables.Strain.Eigen.Values[i]*rVariables.Strain.Eigen.Values[i]) * rVariables.Strain.Matrix + (rVariables.Strain.Invariants.I3 /(rVariables.Strain.Eigen.Values[i]*rVariables.Strain.Eigen.Values[i])) * this->msIdentityMatrix)/D;
263 
264  // noalias(MaxMa)=(rVariables.Strain.Matrix - (rVariables.Strain.Invariants.I1-rVariables.Strain.Eigen.Values[i]*rVariables.Strain.Eigen.Values[i]) * this->msIdentityMatrix + (rVariables.Strain.Invariants.I3 /(rVariables.Strain.Eigen.Values[i]*rVariables.Strain.Eigen.Values[i])) * rVariables.Strain.InverseMatrix)*(rVariables.Strain.Eigen.Values[i]*rVariables.Strain.Eigen.Values[i])/D;
265 
266  // }
267  // std::cout<<" maxma "<<maxma<<std::endl;
268  // std::cout<<" MaxMa "<<MaxMa<<std::endl;
269 
270  // }
271 
272 
273  KRATOS_CATCH(" ")
274  }
275 
276 
277  void CalculateInvariants(HyperElasticDataType& rVariables) override
278  {
279  KRATOS_TRY
280 
281  //invariants
282  rVariables.Strain.Invariants.I1 = rVariables.Strain.Eigen.Values[0] * rVariables.Strain.Eigen.Values[0] +
283  rVariables.Strain.Eigen.Values[1] * rVariables.Strain.Eigen.Values[1] +
284  rVariables.Strain.Eigen.Values[2] * rVariables.Strain.Eigen.Values[2];
285 
286  rVariables.Strain.Invariants.I2 = rVariables.Strain.Eigen.Values[1] * rVariables.Strain.Eigen.Values[1] *
287  rVariables.Strain.Eigen.Values[2] * rVariables.Strain.Eigen.Values[2] +
288 
289  rVariables.Strain.Eigen.Values[2] * rVariables.Strain.Eigen.Values[2] *
290  rVariables.Strain.Eigen.Values[0] * rVariables.Strain.Eigen.Values[0] +
291 
292  rVariables.Strain.Eigen.Values[0] * rVariables.Strain.Eigen.Values[0] *
293  rVariables.Strain.Eigen.Values[1] * rVariables.Strain.Eigen.Values[1];
294 
295  rVariables.Strain.Invariants.I3 = rVariables.Strain.Eigen.Values[0] * rVariables.Strain.Eigen.Values[0] *
296  rVariables.Strain.Eigen.Values[1] * rVariables.Strain.Eigen.Values[1] *
297  rVariables.Strain.Eigen.Values[2] * rVariables.Strain.Eigen.Values[2];
298 
299 
300  //jacobian
301  rVariables.Strain.Invariants.J = rVariables.GetModelData().GetTotalDeformationDet();
302  rVariables.Strain.Invariants.J_13 = std::pow(rVariables.Strain.Invariants.J,(-1.0/3.0));
303 
304 
305  //rVariables.Strain.Invariants.I3 = rVariables.Strain.Invariants.J * rVariables.Strain.Invariants.J; //for volumetric consistency
306 
307  //std::cout<<" Strain.Invariants [I1:"<<rVariables.Strain.Invariants.I1<<" I2:"<<rVariables.Strain.Invariants.I2<<" I3:"<<rVariables.Strain.Invariants.I3<<"] J:"<<rVariables.Strain.Invariants.J<<std::endl;
308  KRATOS_CATCH(" ")
309  }
310 
311 
312 
313  void CalculateAndAddIsochoricStressTensor(HyperElasticDataType& rVariables, MatrixType& rStressMatrix) override
314  {
315  KRATOS_TRY
316 
317  const ModelDataType& rModelData = rVariables.GetModelData();
318  const StressMeasureType& rStressMeasure = rModelData.GetStressMeasure();
319 
320  array_1d<double,3> MainStresses;
321  this->CalculateMainStresses(rVariables,MainStresses);
322 
323  if( rStressMeasure == ConstitutiveModelData::StressMeasureType::StressMeasure_PK2 ){ //Strain.Matrix = RightCauchyGreen (C)
324 
325  array_1d<double,3> EigenVector;
326  for(unsigned int i=0; i<3; i++)
327  {
328  noalias(EigenVector) = matrix_row<const MatrixType>(rVariables.Strain.Eigen.Vectors,i);
329  EigenVector /= rVariables.Strain.Eigen.Values[i];
330  noalias(rStressMatrix) += MainStresses[i] * outer_prod(EigenVector,EigenVector);
331  }
332 
333  }
334  else if( rStressMeasure == ConstitutiveModelData::StressMeasureType::StressMeasure_Kirchhoff ){ //Strain.Matrix = LeftCauchyGreen (b)
335 
336  array_1d<double,3> EigenVector;
337  for(unsigned int i=0; i<3; i++)
338  {
339  noalias(EigenVector) = matrix_row<const MatrixType>(rVariables.Strain.Eigen.Vectors,i);
340  noalias(rStressMatrix) += MainStresses[i] * outer_prod(EigenVector,EigenVector);
341  }
342 
343  }
344 
345 
346  KRATOS_CATCH(" ")
347  }
348 
349 
350  void CalculateAndAddVolumetricStressTensor(HyperElasticDataType& rVariables, MatrixType& rStressMatrix) override
351  {
352  KRATOS_TRY
353 
354  const ModelDataType& rModelData = rVariables.GetModelData();
355  const StressMeasureType& rStressMeasure = rModelData.GetStressMeasure();
356 
357  MatrixType StressMatrix;
358  if( rStressMeasure == ConstitutiveModelData::StressMeasureType::StressMeasure_PK2 ){ //Variables.Strain.Matrix = RightCauchyGreen (C)
359  StressMatrix = GetJRightCauchyGreenDerivative(rVariables.Strain,StressMatrix);
360  StressMatrix *= rVariables.Factors.Alpha4;
361 
362  StressMatrix *= 2.0;
363 
364  noalias(rStressMatrix) += StressMatrix;
365  }
366  else if( rStressMeasure == ConstitutiveModelData::StressMeasureType::StressMeasure_Kirchhoff ){ //Variables.Strain.Matrix = LeftCauchyGreen (b)
367 
368  StressMatrix = GetJLeftCauchyGreenDerivative(rVariables.Strain,StressMatrix);
369  StressMatrix *= rVariables.Factors.Alpha4;
370 
371  StressMatrix *= 2.0;
372 
373  noalias(rStressMatrix) += StressMatrix;
374  }
375 
376  KRATOS_CATCH(" ")
377  }
378 
379 
380  //************************************************************************************
381  //************************************************************************************
382 
383  void CalculateMainStresses(HyperElasticDataType& rVariables, array_1d<double,3>& rMainStresses) override
384  {
385  KRATOS_TRY
386 
387 
388  const MaterialDataType& rMaterial = rVariables.GetMaterialParameters();
389  const std::vector<double>& rModelParameters = rMaterial.GetModelParameters(); //nu values, lambda values
390 
391  unsigned int size = (rModelParameters.size()/2.0);
392  double athird = 1.0/3.0;
393 
394  for(unsigned int i=0; i<3; i++)
395  {
396  for(unsigned int p=0; p<size; p++)
397  {
398  const double& mu_p = rModelParameters[p];
399  const double& alpha_p = rModelParameters[p+size];
400  rMainStresses[i] += (mu_p) * ( std::pow(rVariables.Strain.Eigen.Values[i],alpha_p) - athird * ( std::pow(rVariables.Strain.Eigen.Values[0],alpha_p) + std::pow(rVariables.Strain.Eigen.Values[1],alpha_p) + std::pow(rVariables.Strain.Eigen.Values[2],alpha_p) ) );
401  }
402 
403  }
404 
405  KRATOS_CATCH(" ")
406  }
407 
408 
409  //************************************************************************************
410  //************************************************************************************
411 
412  void CalculateMainStressDerivatives(HyperElasticDataType& rVariables, MatrixType& rStressDerivatives) override
413  {
414  KRATOS_TRY
415 
416  //Isochoric eigenvalues
417 
418  const MaterialDataType& rMaterial = rVariables.GetMaterialParameters();
419  const std::vector<double>& rModelParameters = rMaterial.GetModelParameters(); //nu values, lambda values
420 
421  unsigned int size = (rModelParameters.size()/2.0);
422  double athird = 1.0/3.0;
423 
424  for(unsigned int i=0; i<3; i++)
425  {
426  for(unsigned int j=0; j<3; j++)
427  {
428 
429  for(unsigned int p=0; p<size; p++)
430  {
431  const double& mu_p = rModelParameters[p];
432  const double& alpha_p = rModelParameters[p+size];
433 
434  rStressDerivatives(i,j) += mu_p * alpha_p * ( athird * (std::pow(rVariables.Strain.Eigen.Values[i],alpha_p) + athird * ( std::pow(rVariables.Strain.Eigen.Values[0],alpha_p) + std::pow(rVariables.Strain.Eigen.Values[1],alpha_p) + std::pow(rVariables.Strain.Eigen.Values[2],alpha_p) ) ) );
435 
436  if( i != j ){
437  rStressDerivatives(i,j) -= mu_p * alpha_p * ( athird * ( 2.0 * std::pow(rVariables.Strain.Eigen.Values[i],alpha_p) + std::pow(rVariables.Strain.Eigen.Values[j],alpha_p) ) );
438  }
439  }
440  }
441 
442  }
443 
444 
445  KRATOS_CATCH(" ")
446  }
447 
448  //************************************************************************************
449  //************************************************************************************
450 
451  virtual void CalculateAndAddPerturbedConstitutiveTensor(HyperElasticDataType& rVariables, Matrix& rConstitutiveMatrix)
452  {
453  KRATOS_TRY
454 
455  ModelDataType Values = rVariables.GetModelData();
456 
457  // double& TotalDeterminant = Values.rConstitutiveLawData().TotalDeformationDet;
458  MatrixType& DeltaDeformationMatrix = Values.rConstitutiveLawData().DeltaDeformationMatrix;
459  MatrixType& TotalDeformationMatrix = Values.rConstitutiveLawData().TotalDeformationMatrix;
460 
461  MatrixType StressMatrix;
462  noalias(StressMatrix) = ZeroMatrix(3,3);
463 
464  const SizeType& rVoigtSize = Values.GetVoigtSize();
465  const VoigtIndexType& rIndexVoigtTensor = Values.GetVoigtIndexTensor();
466 
467  Vector StressVectorI(rVoigtSize);
468  Vector StressVectorII(rVoigtSize);
469 
470  double value = 0;
471  for( unsigned int i=0; i<rVoigtSize; i++)
472  {
473  value = rVariables.GetModelData().GetDeltaDeformationMatrix()(rIndexVoigtTensor[i][0],rIndexVoigtTensor[i][1]);
474  double deltavalue = 1e-10;
475  if( value !=0 )
476  deltavalue = value * 1e-8;
477 
478 
479  //Calculate stress
480  DeltaDeformationMatrix = rVariables.GetModelData().GetDeltaDeformationMatrix();
481  TotalDeformationMatrix = rVariables.GetModelData().GetTotalDeformationMatrix();
482 
483  DeltaDeformationMatrix(rIndexVoigtTensor[i][0],rIndexVoigtTensor[i][1]) += deltavalue;
484  //TotalDeformationMatrix(rIndexVoigtTensor[i][0],rIndexVoigtTensor[i][1]) += deltavalue;
485  //TotalDeterminant = MathUtils<double>::Det(TotalDeformationMatrix);
486 
487  //std::cout<<" Det "<<TotalDeterminant<<" DeltaF "<<DeltaDeformationMatrix<<" TotalDet "<<TotalDeformationMatrix<<std::endl;
488 
489  this->CalculateStressTensor(Values, StressMatrix);
490  StressVectorI = ConstitutiveModelUtilities::StressTensorToVector(StressMatrix, StressVectorI);
491 
492  //Calculate elemental system
493  DeltaDeformationMatrix = rVariables.GetModelData().GetDeltaDeformationMatrix();
494  TotalDeformationMatrix = rVariables.GetModelData().GetTotalDeformationMatrix();
495 
496  DeltaDeformationMatrix(rIndexVoigtTensor[i][0],rIndexVoigtTensor[i][1]) -= deltavalue;
497  //TotalDeformationMatrix(rIndexVoigtTensor[i][0],rIndexVoigtTensor[i][1]) -= deltavalue;
498  //TotalDeterminant = MathUtils<double>::Det(TotalDeformationMatrix);
499 
500  this->CalculateStressTensor(Values, StressMatrix);
501  StressVectorII = ConstitutiveModelUtilities::StressTensorToVector(StressMatrix, StressVectorII);
502 
503  //std::cout<<" StressVector I "<<StressVectorI<<std::endl;
504  //std::cout<<" StressVector II "<<StressVectorII<<std::endl;
505 
506  for( unsigned int j=0; j<rVoigtSize; j++)
507  {
508  rConstitutiveMatrix(j,i) = (-1) * (StressVectorI[j] - StressVectorII[j]) / (2.0*deltavalue);
509  }
510 
511  }
512 
513  //std::cout<<" PerturbedConstitutiveMatrix "<<rConstitutiveMatrix<<std::endl;
514 
515  KRATOS_CATCH(" ")
516  }
517 
518  //************************************************************************************
519  //************************************************************************************
520 
521  void CalculateAndAddConstitutiveTensor(HyperElasticDataType& rVariables, Matrix& rConstitutiveMatrix) override
522  {
523  KRATOS_TRY
524 
525  //Calculate Ogden ConstitutiveMatrix
526  const ModelDataType& rModelData = rVariables.GetModelData();
527  const SizeType& rVoigtSize = rModelData.GetVoigtSize();
528  const VoigtIndexType& rIndexVoigtTensor = rModelData.GetVoigtIndexTensor();
529 
530  //Calculate Stress main streches derivatives
531  MatrixType StressDerivatives;
532  noalias(StressDerivatives)=ZeroMatrix(3,3);
533  this->CalculateMainStressDerivatives(rVariables, StressDerivatives);
534 
535  array_1d<double,3> StressEigenValues;
536  noalias(StressEigenValues)=ZeroVector(3);
537  this->CalculateMainStresses(rVariables, StressEigenValues);
538 
539  //Calculate constitutive components
540  for(SizeType i=0; i<rVoigtSize; i++)
541  {
542  for(SizeType j=0; j<rVoigtSize; j++)
543  {
544  // rConstitutiveMatrix(i,j) = this->AddIsochoricConstitutiveComponent(rVariables,rConstitutiveMatrix(i,j),
545  // StressDerivatives,StressEigenValues,
546  // rIndexVoigtTensor[i][0],rIndexVoigtTensor[i][1],
547  // rIndexVoigtTensor[j][0],rIndexVoigtTensor[j][1]);
548  rConstitutiveMatrix(i,j) = this->AddIsochoricConstitutiveComponent(rVariables,rConstitutiveMatrix(i,j),
549  StressEigenValues,
550  rIndexVoigtTensor[i][0],rIndexVoigtTensor[i][1],
551  rIndexVoigtTensor[j][0],rIndexVoigtTensor[j][1]);
552 
553  //std::cout<<" iso Cij "<<rConstitutiveMatrix(i,j)<<" "<<i<<" "<<j<<std::endl;
554  rConstitutiveMatrix(i,j) = this->AddVolumetricConstitutiveComponent(rVariables,rConstitutiveMatrix(i,j),
555  rIndexVoigtTensor[i][0],rIndexVoigtTensor[i][1],
556  rIndexVoigtTensor[j][0],rIndexVoigtTensor[j][1]);
557  //std::cout<<" vol Cij "<<rConstitutiveMatrix(i,j)<<" "<<i<<" "<<j<<std::endl;
558  }
559 
560  }
561 
562  //std::cout<<" ConstitutiveMatrix "<<rConstitutiveMatrix<<std::endl;
563 
564  rVariables.State().Set(ConstitutiveModelData::CONSTITUTIVE_MATRIX_COMPUTED,true);
565 
566 
567  KRATOS_CATCH(" ")
568  }
569 
570  //************************************************************************************
571  //************************************************************************************
572 
573  double& CalculateStressDerivativesI(HyperElasticDataType& rVariables, double& rValue,
574  const unsigned int& i, const unsigned int& j) override
575  {
576  KRATOS_TRY
577 
578  //Calculate Ogden main stress derivatives
579  const MaterialDataType& rMaterial = rVariables.GetMaterialParameters();
580 
581  const std::vector<double>& rModelParameters = rMaterial.GetModelParameters(); //nu values, lambda values
582 
583  unsigned int size = (rModelParameters.size()/2.0);
584  double athird = 1.0/3.0;
585 
586  rValue = 0;
587  for(unsigned int p=0; p<size; p++)
588  {
589  const double& mu_p = rModelParameters[p];
590  const double& alpha_p = rModelParameters[p+size];
591  double f = athird * ( std::pow(rVariables.Strain.Eigen.Values[0],alpha_p) + std::pow(rVariables.Strain.Eigen.Values[1],alpha_p) + std::pow(rVariables.Strain.Eigen.Values[2],alpha_p) );
592 
593  rValue += athird * (mu_p * alpha_p * ( f - std::pow(rVariables.Strain.Eigen.Values[i],alpha_p) - std::pow(rVariables.Strain.Eigen.Values[j],alpha_p) + 3.0 * std::pow(rVariables.Strain.Eigen.Values[i],alpha_p) * this->msIdentityMatrix(i,j) ) );
594 
595  }
596 
597  return rValue;
598 
599  KRATOS_CATCH(" ")
600  }
601 
602  //************************************************************************************
603  //************************************************************************************
604 
605  double& CalculateStressDerivativesII(HyperElasticDataType& rVariables, double& rValue,
606  const unsigned int& i, const unsigned int& j) override
607  {
608  KRATOS_TRY
609 
610  //Calculate Ogden main stress derivatives
611  const MaterialDataType& rMaterial = rVariables.GetMaterialParameters();
612  const std::vector<double>& rModelParameters = rMaterial.GetModelParameters(); //nu values, lambda values
613 
614  unsigned int size = (rModelParameters.size()/2.0);
615 
616  rValue = 0;
617  for(unsigned int p=0; p<size; p++)
618  {
619  const double& mu_p = rModelParameters[p];
620  const double& alpha_p = rModelParameters[p+size];
621 
622  rValue += 0.5 * mu_p * alpha_p * std::pow(rVariables.Strain.Eigen.Values[i],alpha_p) * ( 1.0 - this->msIdentityMatrix(i,j) );
623 
624  }
625 
626  return rValue;
627 
628  KRATOS_CATCH(" ")
629  }
630 
631  //************************************************************************************
632  //************************************************************************************
633 
634  virtual double& AddIsochoricConstitutiveComponent(HyperElasticDataType& rVariables, double &rCabcd,
635  const array_1d<double,3>& rStressEigenValues,
636  const unsigned int& a, const unsigned int& b,
637  const unsigned int& c, const unsigned int& d) //do not override
638  {
639  KRATOS_TRY
640 
641  //Calculate Ogden ConstitutiveMatrix
642  double Cabcd = 0;
643  if( a == b && c == d ){
644 
645  Cabcd = CalculateStressDerivativesI(rVariables,Cabcd,a,c);
646  rCabcd += Cabcd - 2.0 * rStressEigenValues[a] * this->msIdentityMatrix(a,c);
647  }
648  else if( a == c && b == d ){
649 
650  Cabcd = CalculateStressDerivativesII(rVariables,Cabcd,a,b);
651  rCabcd = Cabcd - rStressEigenValues[a];
652 
653  }
654 
655  return rCabcd;
656 
657  KRATOS_CATCH(" ")
658  }
659 
660 
661  // virtual double& AddIsochoricConstitutiveComponent(HyperElasticDataType& rVariables, double &rCabcd,
662  // const MatrixType& rStressDerivatives, const array_1d<double,3>& rStressEigenValues,
663  // const unsigned int& a, const unsigned int& b,
664  // const unsigned int& c, const unsigned int& d) //do not override
665  // {
666  // KRATOS_TRY
667 
668  // const ModelDataType& rModelData = rVariables.GetModelData();
669  // const StressMeasureType& rStressMeasure = rModelData.GetStressMeasure();
670 
671  // double Dabcd = 0;
672  // double Cabcd = 0;
673 
674  // unsigned int option = 0;
675  // array_1d<unsigned int,3> Order;
676 
677  // this->GetEigenCoincidence(rVariables.Strain.Eigen.Values,Order,option);
678 
679  // if( option == 1 ){ //all eigen values are the different
680 
681  // array_1d<double,3> EigenVectorA;
682  // array_1d<double,3> EigenVectorB;
683 
684  // if( rStressMeasure == ConstitutiveModelData::StressMeasureType::StressMeasure_PK2 ){ //Variables.Strain.Matrix = RightCauchyGreen (C)
685  // for(unsigned int i=0; i<3; i++)
686  // {
687  // noalias(EigenVectorA) = matrix_row<const MatrixType>(rVariables.Strain.Eigen.Vectors,i);
688  // EigenVectorA /= rVariables.Strain.Eigen.Values[i];
689  // for(unsigned int j=0; j<3; j++)
690  // {
691 
692  // noalias(EigenVectorB) = matrix_row<const MatrixType>(rVariables.Strain.Eigen.Vectors,j);
693  // EigenVectorB /= rVariables.Strain.Eigen.Values[j];
694 
695  // Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(EigenVectorA,EigenVectorB,Dabcd,a,b,c,d);
696 
697  // Cabcd += rStressDerivatives(i,j) * Dabcd;
698  // }
699 
700  // Dabcd = GetEigenProductRightCauchyGreenDerivative(rVariables,i,Dabcd,a,b,c,d);
701  // Cabcd += 2.0 * rStressEigenValues[i] * Dabcd;
702  // //std::cout<<" Cabcd "<<Cabcd<<" Dabcd "<<Dabcd<<" "<<a<<" "<<b<<" "<<c<<" "<<d<<std::endl;
703  // }
704 
705  // }
706  // else if( rStressMeasure == ConstitutiveModelData::StressMeasureType::StressMeasure_Kirchhoff ){ //Variables.Strain.M
707 
708  // for(unsigned int i=0; i<3; i++)
709  // {
710  // noalias(EigenVectorA) = matrix_row<const MatrixType>(rVariables.Strain.Eigen.Vectors,i);
711  // for(unsigned int j=0; j<3; j++)
712  // {
713  // noalias(EigenVectorB) = matrix_row<const MatrixType>(rVariables.Strain.Eigen.Vectors,j);
714  // Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(EigenVectorA,EigenVectorB,Dabcd,a,b,c,d);
715  // Cabcd += rStressDerivatives(i,j) * Dabcd;
716  // }
717 
718  // Dabcd = GetEigenProductLeftCauchyGreenDerivative(rVariables,i,Dabcd,a,b,c,d);
719  // Cabcd += 2.0 * rStressEigenValues[i] * Dabcd;
720  // //std::cout<<" Cabcd "<<Cabcd<<" Dabcd "<<Dabcd<<" "<<a<<" "<<b<<" "<<c<<" "<<d<<std::endl;
721  // }
722 
723  // }
724  // }
725  // else if( option == 2 ){ //some eigen values are the same some are different
726 
727  // //std::cout<<" option 2 active "<<std::endl;
728 
729  // array_1d<double,3> EigenVector;
730  // MatrixType EigenOperation;
731 
732  // if( rStressMeasure == ConstitutiveModelData::StressMeasureType::StressMeasure_PK2 ){ //Variables.Strain.Matrix = RightCauchyGreen (C)
733  // noalias(EigenVector) = matrix_row<const MatrixType>(rVariables.Strain.Eigen.Vectors,rStressEigenValues[Order[0]]);
734  // EigenVector /= rVariables.Strain.Eigen.Values[Order[0]];
735 
736  // noalias(EigenOperation) = this->msIdentityMatrix-outer_prod(EigenVector,EigenVector);
737 
738  // Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderUnitTensor(this->msIdentityMatrix,Dabcd,a,b,c,d);
739  // Cabcd -= 2.0 * rStressEigenValues[Order[2]] * Dabcd;
740 
741  // Dabcd = GetEigenProductRightCauchyGreenDerivative(rVariables,Order[0],Dabcd,a,b,c,d);
742 
743  // Cabcd += 2.0 * (rStressEigenValues[Order[0]]-rStressEigenValues[Order[2]])* Dabcd;
744 
745  // Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(EigenOperation,EigenOperation,Dabcd,a,b,c,d);
746  // Cabcd += rStressDerivatives(Order[2],Order[2]) * Dabcd;
747 
748  // Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(EigenVector,EigenVector,Dabcd,a,b,c,d);
749  // Cabcd += rStressDerivatives(Order[0],Order[0]) * Dabcd;
750 
751  // Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(EigenVector,EigenOperation,Dabcd,a,b,c,d);
752  // Cabcd += rStressDerivatives(Order[2],Order[0]) * Dabcd;
753 
754  // Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(EigenOperation,EigenVector,Dabcd,a,b,c,d);
755  // Cabcd += rStressDerivatives(Order[2],Order[0]) * Dabcd;
756 
757  // }
758  // else if( rStressMeasure == ConstitutiveModelData::StressMeasureType::StressMeasure_Kirchhoff ){ //Variables.Strain.M
759  // noalias(EigenVector) = matrix_row<const MatrixType>(rVariables.Strain.Eigen.Vectors,Order[0]);
760  // noalias(EigenOperation) = this->msIdentityMatrix-outer_prod(EigenVector,EigenVector);
761 
762  // Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderUnitTensor(this->msIdentityMatrix,Dabcd,a,b,c,d);
763  // Cabcd -= 2.0 * rStressEigenValues[Order[2]] * Dabcd;
764 
765  // Dabcd = GetEigenProductLeftCauchyGreenDerivative(rVariables,Order[0],Dabcd,a,b,c,d);
766 
767  // Cabcd += 2.0 * (rStressEigenValues[Order[0]]-rStressEigenValues[Order[2]])* Dabcd;
768 
769  // Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(EigenOperation,EigenOperation,Dabcd,a,b,c,d);
770  // Cabcd += rStressDerivatives(Order[2],Order[2]) * Dabcd;
771 
772  // Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(EigenVector,EigenVector,Dabcd,a,b,c,d);
773  // Cabcd += rStressDerivatives(Order[0],Order[0]) * Dabcd;
774 
775  // Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(EigenVector,EigenOperation,Dabcd,a,b,c,d);
776  // Cabcd += rStressDerivatives(Order[2],Order[0]) * Dabcd;
777 
778  // Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(EigenOperation,EigenVector,Dabcd,a,b,c,d);
779  // Cabcd += rStressDerivatives(Order[2],Order[0]) * Dabcd;
780  // }
781 
782  // }
783  // else if( option == 3 ){ //all eigen values are the same
784 
785  // const MaterialDataType& rMaterial = rVariables.GetMaterialParameters();
786  // const std::vector<double>& rModelParameters = rMaterial.GetModelParameters(); //nu values, lambda values
787 
788  // unsigned int size = (rModelParameters.size()/2.0);
789  // double Gamma = 0;
790  // for(unsigned int p=0; p<size; p++)
791  // {
792  // const double& mu_p = rModelParameters[p];
793  // const double& alpha_p = rModelParameters[p+size];
794 
795  // Gamma += mu_p * std::pow(rVariables.Strain.Eigen.Values[0],alpha_p);
796  // }
797 
798  // if( rStressMeasure == ConstitutiveModelData::StressMeasureType::StressMeasure_PK2 ){ //Variables.Strain.Matrix = RightCauchyGreen (C)
799 
800  // // Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensor(rVariables.Strain.InverseMatrix,Dabcd,a,b,c,d);
801  // // rCabcd -= Dabcd;
802 
803  // // Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(rVariables.Strain.InverseMatrix,rVariables.Strain.InverseMatrix,Dabcd,a,b,c,d);
804  // // rCabcd += (1.0/3.0) * Dabcd;
805  // Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderUnitTensor(this->msIdentityMatrix,Dabcd,a,b,c,d);
806  // Cabcd += Dabcd;
807 
808  // Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(this->msIdentityMatrix,this->msIdentityMatrix,Dabcd,a,b,c,d);
809  // Cabcd -= (1.0/3.0) * Dabcd;
810 
811  // Cabcd *= Gamma;
812  // }
813  // else if( rStressMeasure == ConstitutiveModelData::StressMeasureType::StressMeasure_Kirchhoff ){ //Variables.Strain.M
814 
815  // Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderUnitTensor(this->msIdentityMatrix,Dabcd,a,b,c,d);
816  // Cabcd += Dabcd;
817 
818  // Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(this->msIdentityMatrix,this->msIdentityMatrix,Dabcd,a,b,c,d);
819  // Cabcd -= (1.0/3.0) * Dabcd;
820 
821  // Cabcd *= Gamma;
822  // }
823  // }
824 
825  // rCabcd += Cabcd;
826 
827  // return rCabcd;
828 
829  // KRATOS_CATCH(" ")
830  // }
831 
832 
833 
834  virtual double& GetEigenProductRightCauchyGreenDerivative(HyperElasticDataType& rVariables, const unsigned int& i, double &rCabcd,
835  const unsigned int& a, const unsigned int& b,
836  const unsigned int& c, const unsigned int& d)
837  {
838  KRATOS_TRY
839 
840  const double& lambda = rVariables.Strain.Eigen.Values[i];
841 
842  double D = 2.0 * lambda*lambda*lambda*lambda - rVariables.Strain.Invariants.I1 * lambda*lambda + rVariables.Strain.Invariants.I3 / (lambda*lambda);
843 
844  double dD = 8.0 * lambda*lambda*lambda - 2.0 * rVariables.Strain.Invariants.I1 * lambda - 2.0 * rVariables.Strain.Invariants.I3 / (lambda*lambda*lambda);
845 
846  double Dabcd = 0;
847 
848  array_1d<double,3> EigenVector;
849  noalias(EigenVector) = matrix_row<const MatrixType>(rVariables.Strain.Eigen.Vectors,i);
850  EigenVector /= rVariables.Strain.Eigen.Values[i];
851 
852  Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderUnitTensor(this->msIdentityMatrix,Dabcd,a,b,c,d);
853  rCabcd += Dabcd;
854  Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(this->msIdentityMatrix,this->msIdentityMatrix,Dabcd,a,b,c,d);
855  rCabcd -= Dabcd;
856  Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(rVariables.Strain.InverseMatrix,rVariables.Strain.InverseMatrix,Dabcd,a,b,c,d);
857  rCabcd += Dabcd * rVariables.Strain.Invariants.I3 / (lambda*lambda);
858  Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensor(rVariables.Strain.InverseMatrix,Dabcd,a,b,c,d);
859  rCabcd -= Dabcd * rVariables.Strain.Invariants.I3 / (lambda*lambda);
860  Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(this->msIdentityMatrix,EigenVector,Dabcd,a,b,c,d);
861  rCabcd += (lambda*lambda) * Dabcd;
862  Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(EigenVector,this->msIdentityMatrix,Dabcd,a,b,c,d);
863  rCabcd += (lambda*lambda) * Dabcd;
864  Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(EigenVector,EigenVector,Dabcd,a,b,c,d);
865  rCabcd -= 0.5 * dD * lambda * Dabcd;
866  Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(rVariables.Strain.InverseMatrix,EigenVector,Dabcd,a,b,c,d);
867  rCabcd -= rVariables.Strain.Invariants.I3 * Dabcd / (lambda*lambda);
868  Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(EigenVector,rVariables.Strain.InverseMatrix,Dabcd,a,b,c,d);
869  rCabcd -= rVariables.Strain.Invariants.I3 * Dabcd / (lambda*lambda);
870 
871  if( D != 0)
872  rCabcd /= D;
873 
874  return rCabcd;
875 
876  KRATOS_CATCH(" ")
877 
878  }
879 
880  virtual double& GetEigenProductLeftCauchyGreenDerivative(HyperElasticDataType& rVariables, const unsigned int& i, double &rCabcd,
881  const unsigned int& a, const unsigned int& b,
882  const unsigned int& c, const unsigned int& d)
883  {
884  KRATOS_TRY
885 
886  const double& lambda = rVariables.Strain.Eigen.Values[i];
887 
888  double D = 2.0 * lambda*lambda*lambda*lambda - rVariables.Strain.Invariants.I1 * lambda*lambda + rVariables.Strain.Invariants.I3 / (lambda*lambda);
889 
890  double dD = 8.0 * lambda*lambda*lambda - 2.0 * rVariables.Strain.Invariants.I1 * lambda - 2.0 * rVariables.Strain.Invariants.I3 / (lambda*lambda*lambda);
891 
892  double Dabcd = 0;
893 
894  array_1d<double,3> EigenVector;
895  noalias(EigenVector) = matrix_row<const MatrixType>(rVariables.Strain.Eigen.Vectors,i);
896 
897  Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensor(rVariables.Strain.Matrix,Dabcd,a,b,c,d);
898  rCabcd += Dabcd;
899  Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(rVariables.Strain.Matrix,rVariables.Strain.Matrix,Dabcd,a,b,c,d);
900  rCabcd -= Dabcd;
901  Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(this->msIdentityMatrix,this->msIdentityMatrix,Dabcd,a,b,c,d);
902  rCabcd += Dabcd * rVariables.Strain.Invariants.I3 / (lambda*lambda);
903  Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderUnitTensor(this->msIdentityMatrix,Dabcd,a,b,c,d);
904  rCabcd -= Dabcd * rVariables.Strain.Invariants.I3 / (lambda*lambda);
905  Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(rVariables.Strain.Matrix,EigenVector,Dabcd,a,b,c,d);
906  rCabcd += (lambda*lambda) * Dabcd;
907  Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(EigenVector,rVariables.Strain.Matrix,Dabcd,a,b,c,d);
908  rCabcd += (lambda*lambda) * Dabcd;
909  Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(EigenVector,EigenVector,Dabcd,a,b,c,d);
910  rCabcd -= 0.5 * dD * lambda * Dabcd;
911  Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(this->msIdentityMatrix,EigenVector,Dabcd,a,b,c,d);
912  rCabcd -= rVariables.Strain.Invariants.I3 * Dabcd / (lambda*lambda);
913  Dabcd = ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct(EigenVector,this->msIdentityMatrix,Dabcd,a,b,c,d);
914  rCabcd -= rVariables.Strain.Invariants.I3 * Dabcd / (lambda*lambda);
915 
916  if( D != 0)
917  rCabcd /= D;
918 
919  return rCabcd;
920 
921  KRATOS_CATCH(" ")
922 
923  }
924 
925  double& AddVolumetricConstitutiveComponent(HyperElasticDataType& rVariables, double &rCabcd,
926  const unsigned int& a, const unsigned int& b,
927  const unsigned int& c, const unsigned int& d) override
928  {
929  KRATOS_TRY
930 
931  const ModelDataType& rModelData = rVariables.GetModelData();
932  const StressMeasureType& rStressMeasure = rModelData.GetStressMeasure();
933 
934  double Dabcd = 0;
935  double Cabcd = 0;
936 
937  if( rStressMeasure == ConstitutiveModelData::StressMeasureType::StressMeasure_PK2 ){ //Variables.Strain.Matrix = RightCauchyGreen (C)
938 
939  //2nd derivatives
940  Dabcd = GetJRightCauchyGreen2ndDerivative(rVariables.Strain,Dabcd,a,b,c,d);
941  Cabcd += rVariables.Factors.Alpha4 * Dabcd;
942 
943  //1st derivatives
944  Dabcd = GetJRightCauchyGreenSquare1stDerivative(rVariables.Strain,Dabcd,a,b,c,d);
945  Cabcd += rVariables.Factors.Beta4 * Dabcd;
946 
947  Cabcd *= 4.0;
948  }
949  else if( rStressMeasure == ConstitutiveModelData::StressMeasureType::StressMeasure_Kirchhoff ){ //Variables.Strain.Matrix = LeftCauchyGreen (b)
950  //2nd derivatives
951  Dabcd = GetJLeftCauchyGreen2ndDerivative(rVariables.Strain,Dabcd,a,b,c,d);
952  Cabcd += rVariables.Factors.Alpha4 * Dabcd;
953 
954  //1st derivatives
955  Dabcd = GetJLeftCauchyGreenSquare1stDerivative(rVariables.Strain,Dabcd,a,b,c,d);
956  Cabcd += rVariables.Factors.Beta4 * Dabcd;
957 
958  Cabcd *= 4.0;
959 
960  }
961 
962  rCabcd += Cabcd;
963 
964  return rCabcd;
965 
966  KRATOS_CATCH(" ")
967  }
968 
969 
970 
971  void CalculateScalingFactors(HyperElasticDataType& rVariables) override
972  {
973  KRATOS_TRY
974 
975  rVariables.Factors.Alpha4 = this->GetVolumetricFunction1stJDerivative(rVariables,rVariables.Factors.Alpha4);
976  rVariables.Factors.Beta4 = this->GetVolumetricFunction2ndJDerivative(rVariables,rVariables.Factors.Beta4);
977 
978  KRATOS_CATCH(" ")
979  }
980 
981 
982  //************// W
983 
984  void CalculateAndAddIsochoricStrainEnergy(HyperElasticDataType& rVariables, double& rIsochoricDensityFunction) override
985  {
986  KRATOS_TRY
987 
988  const MaterialDataType& rMaterial = rVariables.GetMaterialParameters();
989 
990  const std::vector<double>& rModelParameters = rMaterial.GetModelParameters(); //nu values, lambda values
991 
992  unsigned int size = (rModelParameters.size()/2.0);
993 
994  for(unsigned int p=0; p<size; p++)
995  {
996  const double& mu_p = rModelParameters[p];
997  const double& alpha_p = rModelParameters[p+size];
998  rIsochoricDensityFunction += (mu_p/alpha_p) * ( std::pow(rVariables.Strain.Eigen.Values[0],alpha_p) + std::pow(rVariables.Strain.Eigen.Values[1],alpha_p) + std::pow(rVariables.Strain.Eigen.Values[2],alpha_p) - 3.0 );
999  }
1000 
1001  KRATOS_CATCH(" ")
1002  }
1003 
1004 
1005  void CalculateAndAddVolumetricStrainEnergy(HyperElasticDataType& rVariables, double& rVolumetricDensityFunction) override
1006  {
1007  KRATOS_TRY
1008 
1009  const MaterialDataType& rMaterial = rVariables.GetMaterialParameters();
1010 
1011  //energy function "U(J) = (K/2)*(lnJ)²"
1012  rVolumetricDensityFunction += rMaterial.GetBulkModulus() * 0.5 * pow(std::log(rVariables.Strain.Invariants.J),2);
1013 
1014  KRATOS_CATCH(" ")
1015  }
1016 
1017  //************// dW
1018 
1019  double& GetVolumetricFunction1stJDerivative(HyperElasticDataType& rVariables, double& rDerivative) override //dU/dJ
1020  {
1021  KRATOS_TRY
1022 
1023  // const ModelDataType& rValues = rVariables.GetModelData();
1024 
1025  // rDerivative = rValues.GetPressure();
1026 
1027  // return rDerivative;
1028 
1029  const MaterialDataType& rMaterial = rVariables.GetMaterialParameters();
1030 
1031  //derivative of "U(J) = (K/2)*ln(J)²"
1032  //dU(J)/dJ = (K)*(lnJ/J)
1033  rDerivative = rMaterial.GetBulkModulus() * std::log( rVariables.Strain.Invariants.J );
1034 
1035  rDerivative /= rVariables.Strain.Invariants.J;
1036 
1037  return rDerivative;
1038 
1039  KRATOS_CATCH(" ")
1040  };
1041 
1042 
1043  double& GetVolumetricFunction2ndJDerivative(HyperElasticDataType& rVariables, double& rDerivative) override //ddU/dJdJ
1044  {
1045  KRATOS_TRY
1046 
1047  // rDerivative = 0.0;
1048 
1049  // return rDerivative;
1050 
1051  const MaterialDataType& rMaterial = rVariables.GetMaterialParameters();
1052 
1053  //derivative of "dU(J)/dJ = (K)*(lnJ/J)"
1054  //ddU(J)/dJdJ = (K)*(1-lnJ)/J²
1055  rDerivative = rMaterial.GetBulkModulus() * (1.0 -std::log(rVariables.Strain.Invariants.J)) / (rVariables.Strain.Invariants.J * rVariables.Strain.Invariants.J);
1056 
1057  return rDerivative;
1058 
1059  KRATOS_CATCH(" ")
1060  };
1061 
1062 
1066 
1067 
1071 
1072 
1076 
1077 
1079 
1080  private:
1081 
1084 
1085 
1089 
1090 
1094 
1095 
1099 
1100  using HyperElasticModel::AddIsochoricConstitutiveComponent;
1101 
1105 
1106 
1110  friend class Serializer;
1111 
1112 
1113  void save(Serializer& rSerializer) const override
1114  {
1115  KRATOS_SERIALIZE_SAVE_BASE_CLASS( rSerializer, OgdenModel )
1116  }
1117 
1118  void load(Serializer& rSerializer) override
1119  {
1120  KRATOS_SERIALIZE_LOAD_BASE_CLASS( rSerializer, OgdenModel )
1121  }
1122 
1126 
1127 
1131 
1133 
1134  }; // Class IsochoricOgdenModel
1135 
1137 
1140 
1141 
1145 
1146 
1148 
1150 
1151 } // namespace Kratos.
1152 
1153 #endif // KRATOS_ISOCHORIC_OGDEN_MODEL_H_INCLUDED defined
Kratos::ConstitutiveModelData::StressMeasureType
StressMeasureType
Definition: constitutive_model_data.hpp:83
Kratos::ConstitutiveModel::VoigtIndexType
const unsigned int(*)[2] VoigtIndexType
Definition: constitutive_model.hpp:55
Kratos::ConstitutiveModel::SizeType
ConstitutiveModelData::SizeType SizeType
Definition: constitutive_model.hpp:66
Kratos::ConstitutiveModelUtilities::CalculateFourthOrderUnitTensor
static double & CalculateFourthOrderUnitTensor(double &rValue, const unsigned int &a, const unsigned int &b, const unsigned int &c, const unsigned int &d)
Definition: constitutive_model_utilities.hpp:1140
Kratos::ConstitutiveModelUtilities::CalculateFourthOrderTensorProduct
static double & CalculateFourthOrderTensorProduct(const MatrixType &rMatrixA, const MatrixType &rMatrixB, double &rValue, const unsigned int &a, const unsigned int &b, const unsigned int &c, const unsigned int &d)
Definition: constitutive_model_utilities.hpp:1250
Kratos::ConstitutiveModelUtilities::StressTensorToVector
static Vector & StressTensorToVector(const MatrixType &rStressTensor, Vector &rStressVector)
Definition: constitutive_model_utilities.hpp:843
Kratos::ConstitutiveModelUtilities::CalculateFourthOrderTensor
static double & CalculateFourthOrderTensor(const MatrixType &rMatrix, double &rValue, const double &a, const double &b, const double &c, const double &d)
Definition: constitutive_model_utilities.hpp:1187
Kratos::Flags::Set
void Set(const Flags ThisFlag)
Definition: flags.cpp:33
Kratos::HyperElasticModel::GetJLeftCauchyGreenSquare1stDerivative
double & GetJLeftCauchyGreenSquare1stDerivative(const StrainData &rStrain, double &rDerivative, const double &a, const double &b, const double &c, const double &d)
Definition: hyper_elastic_model.cpp:662
Kratos::HyperElasticModel::GetJLeftCauchyGreenDerivative
MatrixType & GetJLeftCauchyGreenDerivative(const StrainData &rStrain, MatrixType &rDerivative)
Definition: hyper_elastic_model.cpp:635
Kratos::HyperElasticModel::GetJRightCauchyGreenDerivative
MatrixType & GetJRightCauchyGreenDerivative(const StrainData &rStrain, MatrixType &rDerivative)
Definition: hyper_elastic_model.cpp:565
Kratos::HyperElasticModel::GetJRightCauchyGreenSquare1stDerivative
double & GetJRightCauchyGreenSquare1stDerivative(const StrainData &rStrain, double &rDerivative, const double &a, const double &b, const double &c, const double &d)
dJ/dC
Definition: hyper_elastic_model.cpp:592
Kratos::HyperElasticModel::msIdentityMatrix
static const MatrixType msIdentityMatrix
Definition: hyper_elastic_model.hpp:322
Kratos::HyperElasticModel::GetJLeftCauchyGreen2ndDerivative
double & GetJLeftCauchyGreen2ndDerivative(const StrainData &rStrain, double &rDerivative, const double &a, const double &b, const double &c, const double &d)
Definition: hyper_elastic_model.cpp:679
Kratos::HyperElasticModel::AddIsochoricConstitutiveComponent
virtual double & AddIsochoricConstitutiveComponent(HyperElasticDataType &rVariables, double &rCabcd, const unsigned int &a, const unsigned int &b, const unsigned int &c, const unsigned int &d)
Definition: hyper_elastic_model.cpp:468
Kratos::HyperElasticModel::GetJRightCauchyGreen2ndDerivative
double & GetJRightCauchyGreen2ndDerivative(const StrainData &rStrain, double &rDerivative, const double &a, const double &b, const double &c, const double &d)
Definition: hyper_elastic_model.cpp:610
Kratos::Internals::Matrix
Definition: amatrix_interface.h:41
Kratos::IsochoricOgdenModel
Short class definition.
Definition: isochoric_ogden_model.hpp:48
Kratos::IsochoricOgdenModel::CalculateAndAddVolumetricStressTensor
void CalculateAndAddVolumetricStressTensor(HyperElasticDataType &rVariables, MatrixType &rStressMatrix) override
Definition: isochoric_ogden_model.hpp:350
Kratos::IsochoricOgdenModel::Info
std::string Info() const override
Turn back information as a string.
Definition: isochoric_ogden_model.hpp:180
Kratos::IsochoricOgdenModel::CalculateStrainEnergy
void CalculateStrainEnergy(ModelDataType &rValues, double &rDensityFunction) override
Definition: isochoric_ogden_model.hpp:94
Kratos::IsochoricOgdenModel::~IsochoricOgdenModel
~IsochoricOgdenModel() override
Destructor.
Definition: isochoric_ogden_model.hpp:81
Kratos::IsochoricOgdenModel::IsochoricOgdenModel
IsochoricOgdenModel(IsochoricOgdenModel const &rOther)
Copy constructor.
Definition: isochoric_ogden_model.hpp:65
Kratos::IsochoricOgdenModel::CalculateConstitutiveTensor
void CalculateConstitutiveTensor(ModelDataType &rValues, Matrix &rConstitutiveMatrix) override
Definition: isochoric_ogden_model.hpp:129
Kratos::IsochoricOgdenModel::Clone
ConstitutiveModel::Pointer Clone() const override
Clone.
Definition: isochoric_ogden_model.hpp:75
Kratos::IsochoricOgdenModel::CalculateAndAddIsochoricStrainEnergy
void CalculateAndAddIsochoricStrainEnergy(HyperElasticDataType &rVariables, double &rIsochoricDensityFunction) override
Definition: isochoric_ogden_model.hpp:984
Kratos::IsochoricOgdenModel::CalculateScalingFactors
void CalculateScalingFactors(HyperElasticDataType &rVariables) override
Definition: isochoric_ogden_model.hpp:971
Kratos::IsochoricOgdenModel::CalculateAndAddPerturbedConstitutiveTensor
virtual void CalculateAndAddPerturbedConstitutiveTensor(HyperElasticDataType &rVariables, Matrix &rConstitutiveMatrix)
Definition: isochoric_ogden_model.hpp:451
Kratos::IsochoricOgdenModel::CalculateAndAddVolumetricStrainEnergy
void CalculateAndAddVolumetricStrainEnergy(HyperElasticDataType &rVariables, double &rVolumetricDensityFunction) override
Definition: isochoric_ogden_model.hpp:1005
Kratos::IsochoricOgdenModel::CalculateStressAndConstitutiveTensors
void CalculateStressAndConstitutiveTensors(ModelDataType &rValues, MatrixType &rStressMatrix, Matrix &rConstitutiveMatrix) override
Definition: isochoric_ogden_model.hpp:144
Kratos::IsochoricOgdenModel::CalculateAndAddConstitutiveTensor
void CalculateAndAddConstitutiveTensor(HyperElasticDataType &rVariables, Matrix &rConstitutiveMatrix) override
Definition: isochoric_ogden_model.hpp:521
Kratos::IsochoricOgdenModel::CalculateStressDerivativesII
double & CalculateStressDerivativesII(HyperElasticDataType &rVariables, double &rValue, const unsigned int &i, const unsigned int &j) override
Definition: isochoric_ogden_model.hpp:605
Kratos::IsochoricOgdenModel::AddIsochoricConstitutiveComponent
virtual double & AddIsochoricConstitutiveComponent(HyperElasticDataType &rVariables, double &rCabcd, const array_1d< double, 3 > &rStressEigenValues, const unsigned int &a, const unsigned int &b, const unsigned int &c, const unsigned int &d)
Definition: isochoric_ogden_model.hpp:634
Kratos::IsochoricOgdenModel::PrintData
void PrintData(std::ostream &rOStream) const override
Print object's data.
Definition: isochoric_ogden_model.hpp:194
Kratos::IsochoricOgdenModel::operator=
IsochoricOgdenModel & operator=(IsochoricOgdenModel const &rOther)
Assignment operator.
Definition: isochoric_ogden_model.hpp:68
Kratos::IsochoricOgdenModel::GetVolumetricFunction1stJDerivative
double & GetVolumetricFunction1stJDerivative(HyperElasticDataType &rVariables, double &rDerivative) override
Definition: isochoric_ogden_model.hpp:1019
Kratos::IsochoricOgdenModel::CalculateStressTensor
void CalculateStressTensor(ModelDataType &rValues, MatrixType &rStressMatrix) override
Definition: isochoric_ogden_model.hpp:110
Kratos::IsochoricOgdenModel::CalculateStrainData
void CalculateStrainData(ModelDataType &rValues, HyperElasticDataType &rVariables) override
Definition: isochoric_ogden_model.hpp:226
Kratos::IsochoricOgdenModel::PrintInfo
void PrintInfo(std::ostream &rOStream) const override
Print information about this object.
Definition: isochoric_ogden_model.hpp:188
Kratos::IsochoricOgdenModel::CalculateAndAddIsochoricStressTensor
void CalculateAndAddIsochoricStressTensor(HyperElasticDataType &rVariables, MatrixType &rStressMatrix) override
Definition: isochoric_ogden_model.hpp:313
Kratos::IsochoricOgdenModel::CalculateInvariants
void CalculateInvariants(HyperElasticDataType &rVariables) override
Definition: isochoric_ogden_model.hpp:277
Kratos::IsochoricOgdenModel::KRATOS_CLASS_POINTER_DEFINITION
KRATOS_CLASS_POINTER_DEFINITION(IsochoricOgdenModel)
Pointer definition of IsochoricOgdenModel.
Kratos::IsochoricOgdenModel::CalculateMainStresses
void CalculateMainStresses(HyperElasticDataType &rVariables, array_1d< double, 3 > &rMainStresses) override
Definition: isochoric_ogden_model.hpp:383
Kratos::IsochoricOgdenModel::GetVolumetricFunction2ndJDerivative
double & GetVolumetricFunction2ndJDerivative(HyperElasticDataType &rVariables, double &rDerivative) override
Definition: isochoric_ogden_model.hpp:1043
Kratos::IsochoricOgdenModel::AddVolumetricConstitutiveComponent
double & AddVolumetricConstitutiveComponent(HyperElasticDataType &rVariables, double &rCabcd, const unsigned int &a, const unsigned int &b, const unsigned int &c, const unsigned int &d) override
Definition: isochoric_ogden_model.hpp:925
Kratos::IsochoricOgdenModel::CalculateMainStressDerivatives
void CalculateMainStressDerivatives(HyperElasticDataType &rVariables, MatrixType &rStressDerivatives) override
Definition: isochoric_ogden_model.hpp:412
Kratos::IsochoricOgdenModel::GetEigenProductLeftCauchyGreenDerivative
virtual double & GetEigenProductLeftCauchyGreenDerivative(HyperElasticDataType &rVariables, const unsigned int &i, double &rCabcd, const unsigned int &a, const unsigned int &b, const unsigned int &c, const unsigned int &d)
Definition: isochoric_ogden_model.hpp:880
Kratos::IsochoricOgdenModel::CalculateStressDerivativesI
double & CalculateStressDerivativesI(HyperElasticDataType &rVariables, double &rValue, const unsigned int &i, const unsigned int &j) override
Definition: isochoric_ogden_model.hpp:573
Kratos::IsochoricOgdenModel::GetEigenProductRightCauchyGreenDerivative
virtual double & GetEigenProductRightCauchyGreenDerivative(HyperElasticDataType &rVariables, const unsigned int &i, double &rCabcd, const unsigned int &a, const unsigned int &b, const unsigned int &c, const unsigned int &d)
Definition: isochoric_ogden_model.hpp:834
Kratos::IsochoricOgdenModel::IsochoricOgdenModel
IsochoricOgdenModel()
Default constructor.
Definition: isochoric_ogden_model.hpp:62
Kratos::OgdenModel
Short class definition.
Definition: ogden_model.hpp:50
Kratos::OgdenModel::operator=
OgdenModel & operator=(OgdenModel const &rOther)
Assignment operator.
Definition: ogden_model.cpp:47
Kratos::OgdenModel::CalculateStrainData
void CalculateStrainData(ModelDataType &rValues, HyperElasticDataType &rVariables) override
Definition: ogden_model.cpp:209
Kratos::Serializer
The serialization consists in storing the state of an object into a storage format like data file or ...
Definition: serializer.h:123
Kratos::array_1d< double, 3 >
KRATOS_SERIALIZE_SAVE_BASE_CLASS
#define KRATOS_SERIALIZE_SAVE_BASE_CLASS(Serializer, BaseType)
Definition: define.h:812
KRATOS_CATCH
#define KRATOS_CATCH(MoreInfo)
Definition: define.h:110
KRATOS_TRY
#define KRATOS_TRY
Definition: define.h:109
KRATOS_SERIALIZE_LOAD_BASE_CLASS
#define KRATOS_SERIALIZE_LOAD_BASE_CLASS(Serializer, BaseType)
Definition: define.h:815
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::ZeroMatrix
KratosZeroMatrix< double > ZeroMatrix
Definition: amatrix_interface.h:559
Kratos::outer_prod
AMatrix::VectorOuterProductExpression< TExpression1Type, TExpression2Type > outer_prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:582
Kratos::noalias
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
generate_convection_diffusion_explicit_element.f
f
Definition: generate_convection_diffusion_explicit_element.py:112
generate_stokes_twofluid_element.a
a
Definition: generate_stokes_twofluid_element.py:77
generate_total_lagrangian_mixed_volumetric_strain_element.b
b
Definition: generate_total_lagrangian_mixed_volumetric_strain_element.py:31
generate_weakly_compressible_navier_stokes_element.c
c
Definition: generate_weakly_compressible_navier_stokes_element.py:108
ode_solve.load
def load(f)
Definition: ode_solve.py:307
ode_solve.d
int d
Definition: ode_solve.py:397
quadrature.j
int j
Definition: quadrature.py:648
sensitivityMatrix.p
p
Definition: sensitivityMatrix.py:52
tensors::i
integer i
Definition: TensorModule.f:17
testing.run_cpp_mpi_tests.e
e
Definition: run_cpp_mpi_tests.py:31
ogden_model.hpp
Kratos::ConstitutiveModelData::MaterialData
Definition: constitutive_model_data.hpp:92
Kratos::ConstitutiveModelData::MaterialData::GetBulkModulus
const double & GetBulkModulus() const
Definition: constitutive_model_data.hpp:114
Kratos::ConstitutiveModelData::MaterialData::GetModelParameters
const std::vector< double > & GetModelParameters() const
Definition: constitutive_model_data.hpp:116
Kratos::ConstitutiveModelData::ModelData
Definition: constitutive_model_data.hpp:383
Kratos::ConstitutiveModelData::ModelData::GetDeltaDeformationMatrix
const MatrixType & GetDeltaDeformationMatrix() const
Definition: constitutive_model_data.hpp:455
Kratos::ConstitutiveModelData::ModelData::GetStressMeasure
const StressMeasureType & GetStressMeasure() const
Definition: constitutive_model_data.hpp:452
Kratos::ConstitutiveModelData::ModelData::GetVoigtSize
const SizeType & GetVoigtSize() const
Definition: constitutive_model_data.hpp:436
Kratos::ConstitutiveModelData::ModelData::GetTotalDeformationDet
const double & GetTotalDeformationDet() const
Definition: constitutive_model_data.hpp:449
Kratos::ConstitutiveModelData::ModelData::GetTotalDeformationMatrix
const MatrixType & GetTotalDeformationMatrix() const
Definition: constitutive_model_data.hpp:456
Kratos::ConstitutiveModelData::ModelData::StressMatrix
MatrixType StressMatrix
Definition: constitutive_model_data.hpp:401
Kratos::HyperElasticModel::HyperElasticFactors::Alpha4
double Alpha4
Definition: hyper_elastic_model.hpp:71
Kratos::HyperElasticModel::HyperElasticFactors::Beta4
double Beta4
Definition: hyper_elastic_model.hpp:76
Kratos::HyperElasticModel::HyperElasticModelData
Definition: hyper_elastic_model.hpp:108
Kratos::HyperElasticModel::HyperElasticModelData::Strain
StrainData Strain
Definition: hyper_elastic_model.hpp:117
Kratos::HyperElasticModel::HyperElasticModelData::State
Flags & State()
Definition: hyper_elastic_model.hpp:128
Kratos::HyperElasticModel::HyperElasticModelData::Factors
HyperElasticFactors Factors
Definition: hyper_elastic_model.hpp:116
Kratos::HyperElasticModel::HyperElasticModelData::GetMaterialParameters
const MaterialDataType & GetMaterialParameters() const
Definition: hyper_elastic_model.hpp:125
Kratos::HyperElasticModel::HyperElasticModelData::GetModelData
const ModelDataType & GetModelData() const
Definition: hyper_elastic_model.hpp:124
Kratos::HyperElasticModel::StrainData::Eigen
StrainEigenData Eigen
Definition: hyper_elastic_model.hpp:99
Kratos::HyperElasticModel::StrainData::InverseMatrix
MatrixType InverseMatrix
Definition: hyper_elastic_model.hpp:102
Kratos::HyperElasticModel::StrainData::Invariants
StrainInvariants Invariants
Definition: hyper_elastic_model.hpp:98
Kratos::HyperElasticModel::StrainData::Matrix
MatrixType Matrix
Definition: hyper_elastic_model.hpp:101
Kratos::HyperElasticModel::StrainEigenData::Vectors
MatrixType Vectors
Definition: hyper_elastic_model.hpp:91
Kratos::HyperElasticModel::StrainEigenData::Values
array_1d< double, 3 > Values
Definition: hyper_elastic_model.hpp:90
Kratos::HyperElasticModel::StrainInvariants::I2
double I2
Definition: hyper_elastic_model.hpp:57
Kratos::HyperElasticModel::StrainInvariants::J
double J
Definition: hyper_elastic_model.hpp:60
Kratos::HyperElasticModel::StrainInvariants::I1
double I1
Definition: hyper_elastic_model.hpp:56
Kratos::HyperElasticModel::StrainInvariants::J_13
double J_13
Definition: hyper_elastic_model.hpp:61
Kratos::HyperElasticModel::StrainInvariants::I3
double I3
Definition: hyper_elastic_model.hpp:58
Generated on Sat Feb 3 2024 01:54:35 for KratosMultiphysics by doxygen 1.9.1