v0.13.0
SmallTransverselyIsotropic.hpp
Go to the documentation of this file.
1 /** \file SmallStrainTranverslyIsotropic.hpp
2  * \ingroup nonlinear_elastic_elem
3  * \brief Implementation of linear transverse isotropic material.
4  */
5 
6 /* This file is part of MoFEM.
7  * MoFEM is free software: you can redistribute it and/or modify it under
8  * the terms of the GNU Lesser General Public License as published by the
9  * Free Software Foundation, either version 3 of the License, or (at your
10  * option) any later version.
11  *
12  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
13  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  * License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
19 
20 #ifndef __SMALLSTRAINTRANVERSLYISOTROPIC_HPP__
21 #define __SMALLSTRAINTRANVERSLYISOTROPIC_HPP__
22 
23 #ifndef WITH_ADOL_C
24  #error "MoFEM need to be compiled with ADOL-C"
25 #endif
26 
27 /** \brief Hook equation
28  * \ingroup nonlinear_elastic_elem
29  */
30 template<typename TYPE>
32 
34 
35 
36 
37  ublas::matrix<TYPE> sTrain;
38  ublas::vector<TYPE> voightStrain;
40 
43  sTrain.resize(3,3,false);
44  noalias(sTrain) = this->F;
45  for(int dd = 0;dd<3;dd++) {
46  sTrain(dd,dd) -= 1;
47  }
48  sTrain += trans(sTrain);
49  sTrain *= 0.5;
50  voightStrain.resize(6,false);
51  voightStrain[0] = sTrain(0,0);
52  voightStrain[1] = sTrain(1,1);
53  voightStrain[2] = sTrain(2,2);
54  voightStrain[3] = 2*sTrain(0,1);
55  voightStrain[4] = 2*sTrain(1,2);
56  voightStrain[5] = 2*sTrain(2,0);
58  }
59 
60  double nu_p, nu_pz, E_p, E_z, G_zp;
61  ublas::symmetric_matrix<TYPE,ublas::upper> localStiffnessMatrix;
64  double nu_zp=(nu_pz*E_z)/E_p;
65  double delta=((1+nu_p)*(1-nu_p-(2*nu_pz*nu_zp)))/(E_p*E_p*E_z);
66 
67  localStiffnessMatrix.resize(6);
68  localStiffnessMatrix.clear();
71 
74 
75  localStiffnessMatrix(3,3)=E_p/(2*(1+nu_p));
77 
79  }
80 
81  ublas::matrix<TYPE> aARotMat;
82  ublas::vector<TYPE> axVector;
84 
85  /**
86  * \brief Function to Calculate the Rotation Matrix at a given axis and angle of rotation
87 
88  * This function computes the rotational matrix for a given axis of rotation
89  * and angle of rotation about that angle <br>
90 
91  *\param axVector A vector representing the axis of rotation
92  *\param axAngle Angle of rotation along the axis (in radians)
93  */
96 
97  aARotMat.resize(3,3,false);
98  aARotMat.clear();
99 
100  TYPE norm = sqrt(pow(axVector[0],2) + pow(axVector[1],2) + pow(axVector[2],2));
101 
102  aARotMat(0,0) = 1-((1-cos(axAngle))*(pow(axVector[1],2)+pow(axVector[2],2))/pow(norm,2));
103  aARotMat(1,1) = 1-((1-cos(axAngle))*(pow(axVector[0],2)+pow(axVector[2],2))/pow(norm,2));
104  aARotMat(2,2) = 1-((1-cos(axAngle))*(pow(axVector[0],2)+pow(axVector[1],2))/pow(norm,2));
105 
106  aARotMat(0,1) = ((1-cos(axAngle))*axVector[0]*axVector[1]-norm*axVector[2]*sin(axAngle))/pow(norm,2);
107  aARotMat(1,0) = ((1-cos(axAngle))*axVector[0]*axVector[1]+norm*axVector[2]*sin(axAngle))/pow(norm,2);
108 
109  aARotMat(0,2) = ((1-cos(axAngle))*axVector[0]*axVector[2]+norm*axVector[1]*sin(axAngle))/pow(norm,2);
110  aARotMat(2,0) = ((1-cos(axAngle))*axVector[0]*axVector[2]-norm*axVector[1]*sin(axAngle))/pow(norm,2);
111 
112  aARotMat(1,2) = ((1-cos(axAngle))*axVector[1]*axVector[2]-norm*axVector[0]*sin(axAngle))/pow(norm,2);
113  aARotMat(2,1) = ((1-cos(axAngle))*axVector[1]*axVector[2]+norm*axVector[0]*sin(axAngle))/pow(norm,2);
114 
116  }
117 
118  ublas::matrix<TYPE> stressRotMat;
119 
120  /**
121  * \brief Function to Calculate Stress Transformation Matrix
122  * This function computes the stress transformation Matrix at a give axis and angle of rotation <br>
123  * One can also output the axis/angle rotational Matrix
124  */
127 
128  stressRotMat.resize(6,6,false);
129  stressRotMat.clear();
130 
131  stressRotMat(0, 0) = aARotMat(0,0) * aARotMat(0,0);
132  stressRotMat(0, 1) = aARotMat(1,0) * aARotMat(1,0);
133  stressRotMat(0, 2) = aARotMat(2,0) * aARotMat(2,0);
134  stressRotMat(0, 3) = 2.0 * aARotMat(1,0) * aARotMat(0,0);
135  stressRotMat(0, 4) = 2.0 * aARotMat(2,0) * aARotMat(1,0);
136  stressRotMat(0, 5) = 2.0 * aARotMat(0,0) * aARotMat(2,0);
137 
138  stressRotMat(1, 0) = aARotMat(0,1) * aARotMat(0,1);
139  stressRotMat(1, 1) = aARotMat(1,1) * aARotMat(1,1);
140  stressRotMat(1, 2) = aARotMat(2,1) * aARotMat(2,1);
141  stressRotMat(1, 3) = 2.0 * aARotMat(1,1) * aARotMat(0,1);
142  stressRotMat(1, 4) = 2.0 * aARotMat(2,1) * aARotMat(1,1);
143  stressRotMat(1, 5) = 2.0 * aARotMat(0,1) * aARotMat(2,1);
144 
145  stressRotMat(2, 0) = aARotMat(0,2) * aARotMat(0,2);
146  stressRotMat(2, 1) = aARotMat(1,2) * aARotMat(1,2);
147  stressRotMat(2, 2) = aARotMat(2,2) * aARotMat(2,2);
148  stressRotMat(2, 3) = 2.0 * aARotMat(1,2) * aARotMat(0,2);
149  stressRotMat(2, 4) = 2.0 * aARotMat(2,2) * aARotMat(1,2);
150  stressRotMat(2, 5) = 2.0 * aARotMat(0,2) * aARotMat(2,2);
151 
152  stressRotMat(3, 0) = aARotMat(0,1) * aARotMat(0,0);
153  stressRotMat(3, 1) = aARotMat(1,1) * aARotMat(1,0);
154  stressRotMat(3, 2) = aARotMat(2,1) * aARotMat(2,0);
155  stressRotMat(3, 3) = ( aARotMat(1,1) * aARotMat(0,0) + aARotMat(0,1) * aARotMat(1,0) );
156  stressRotMat(3, 4) = ( aARotMat(2,1) * aARotMat(1,0) + aARotMat(1,1) * aARotMat(2,0) );
157  stressRotMat(3, 5) = ( aARotMat(0,1) * aARotMat(2,0) + aARotMat(2,1) * aARotMat(0,0) );
158 
159  stressRotMat(4, 0) = aARotMat(0,2) * aARotMat(0,1);
160  stressRotMat(4, 1) = aARotMat(1,2) * aARotMat(1,1);
161  stressRotMat(4, 2) = aARotMat(2,2) * aARotMat(2,1);
162  stressRotMat(4, 3) = ( aARotMat(1,2) * aARotMat(0,1) + aARotMat(0,2) * aARotMat(1,1) );
163  stressRotMat(4, 4) = ( aARotMat(2,2) * aARotMat(1,1) + aARotMat(1,2) * aARotMat(2,1) );
164  stressRotMat(4, 5) = ( aARotMat(0,2) * aARotMat(2,1) + aARotMat(2,2) * aARotMat(0,1) );
165 
166  stressRotMat(5, 0) = aARotMat(0,0) * aARotMat(0,2);
167  stressRotMat(5, 1) = aARotMat(1,0) * aARotMat(1,2);
168  stressRotMat(5, 2) = aARotMat(2,0) * aARotMat(2,2);
169  stressRotMat(5, 3) = ( aARotMat(1,0) * aARotMat(0,2) + aARotMat(0,0) * aARotMat(1,2) );
170  stressRotMat(5, 4) = ( aARotMat(2,0) * aARotMat(1,2) + aARotMat(1,0) * aARotMat(2,2) );
171  stressRotMat(5, 5) = ( aARotMat(0,0) * aARotMat(2,2) + aARotMat(2,0) * aARotMat(0,2) );
172 
174  }
175 
176  ublas::matrix<TYPE> strainRotMat;
177 
178  /**
179  * \brief Function to Calculate Strain Transformation Matrix<br>
180  * This function computes the strain transformation Matrix at a give axis and angle of rotation <br>
181  * One can also output the axis/angle rotational Matrix
182  */
185 
186  strainRotMat.resize(6,6,false);
187  strainRotMat.clear();
188 
189  strainRotMat(0, 0) = aARotMat(0,0) * aARotMat(0,0);
190  strainRotMat(0, 1) = aARotMat(1,0) * aARotMat(1,0);
191  strainRotMat(0, 2) = aARotMat(2,0) * aARotMat(2,0);
192  strainRotMat(0, 3) = aARotMat(1,0) * aARotMat(0,0);
193  strainRotMat(0, 4) = aARotMat(2,0) * aARotMat(1,0);
194  strainRotMat(0, 5) = aARotMat(0,0) * aARotMat(2,0);
195 
196  strainRotMat(1, 0) = aARotMat(0,1) * aARotMat(0,1);
197  strainRotMat(1, 1) = aARotMat(1,1) * aARotMat(1,1);
198  strainRotMat(1, 2) = aARotMat(2,1) * aARotMat(2,1);
199  strainRotMat(1, 3) = aARotMat(1,1) * aARotMat(0,1);
200  strainRotMat(1, 4) = aARotMat(2,1) * aARotMat(1,1);
201  strainRotMat(1, 5) = aARotMat(0,1) * aARotMat(2,1);
202 
203  strainRotMat(2, 0) = aARotMat(0,2) * aARotMat(0,2);
204  strainRotMat(2, 1) = aARotMat(1,2) * aARotMat(1,2);
205  strainRotMat(2, 2) = aARotMat(2,2) * aARotMat(2,2);
206  strainRotMat(2, 3) = aARotMat(1,2) * aARotMat(0,2);
207  strainRotMat(2, 4) = aARotMat(2,2) * aARotMat(1,2);
208  strainRotMat(2, 5) = aARotMat(0,2) * aARotMat(2,2);
209 
210  strainRotMat(3, 0) = 2.0 * aARotMat(0,1) * aARotMat(0,0);
211  strainRotMat(3, 1) = 2.0 * aARotMat(1,1) * aARotMat(1,0);
212  strainRotMat(3, 2) = 2.0 * aARotMat(2,1) * aARotMat(2,0);
213  strainRotMat(3, 3) = ( aARotMat(1,1) * aARotMat(0,0) + aARotMat(0,1) * aARotMat(1,0) );
214  strainRotMat(3, 4) = ( aARotMat(2,1) * aARotMat(1,0) + aARotMat(1,1) * aARotMat(2,0) );
215  strainRotMat(3, 5) = ( aARotMat(0,1) * aARotMat(2,0) + aARotMat(2,1) * aARotMat(0,0) );
216 
217  strainRotMat(4, 0) = 2.0 * aARotMat(0,2) * aARotMat(0,1);
218  strainRotMat(4, 1) = 2.0 * aARotMat(1,2) * aARotMat(1,1);
219  strainRotMat(4, 2) = 2.0 * aARotMat(2,2) * aARotMat(2,1);
220  strainRotMat(4, 3) = ( aARotMat(1,2) * aARotMat(0,1) + aARotMat(0,2) * aARotMat(1,1) );
221  strainRotMat(4, 4) = ( aARotMat(2,2) * aARotMat(1,1) + aARotMat(1,2) * aARotMat(2,1) );
222  strainRotMat(4, 5) = ( aARotMat(0,2) * aARotMat(2,1) + aARotMat(2,2) * aARotMat(0,1) );
223 
224  strainRotMat(5, 0) = 2.0 * aARotMat(0,0) * aARotMat(0,2);
225  strainRotMat(5, 1) = 2.0 * aARotMat(1,0) * aARotMat(1,2);
226  strainRotMat(5, 2) = 2.0 * aARotMat(2,0) * aARotMat(2,2);
227  strainRotMat(5, 3) = ( aARotMat(1,0) * aARotMat(0,2) + aARotMat(0,0) * aARotMat(1,2) );
228  strainRotMat(5, 4) = ( aARotMat(2,0) * aARotMat(1,2) + aARotMat(1,0) * aARotMat(2,2) );
229  strainRotMat(5, 5) = ( aARotMat(0,0) * aARotMat(2,2) + aARotMat(2,0) * aARotMat(0,2) );
230 
232  }
233 
234  ublas::matrix<TYPE> dR;
235  ublas::matrix<TYPE> globalStiffnessMatrix;
236 
239 
240  dR.resize(6,6,false);
241  noalias(dR) = prod(localStiffnessMatrix,strainRotMat);
242  globalStiffnessMatrix.resize(6,6,false);
243  noalias(globalStiffnessMatrix) = prod(stressRotMat,dR);
244 
246  }
247 
251  }
252 
253  ublas::vector<TYPE> voigtStress;
254 
255  /** \brief Calculate global stress
256 
257  This is small strain approach, i.e. Piola stress
258  is like a Cauchy stress, since configurations are notation
259  distinguished.
260 
261  */
263  const NonlinearElasticElement::BlockData block_data,
264  boost::shared_ptr<const NumeredEntFiniteElement> fe_ptr
265  ) {
272  axAngle = -axAngle;
276 
277  voigtStress.resize(6,false);
279  this->P.resize(3,3,false);
280  this->P(0,0) = voigtStress[0];
281  this->P(1,1) = voigtStress[1];
282  this->P(2,2) = voigtStress[2];
283  this->P(0,1) = voigtStress[3];
284  this->P(1,2) = voigtStress[4];
285  this->P(0,2) = voigtStress[5];
286  this->P(1,0) = this->P(0,1);
287  this->P(2,1) = this->P(1,2);
288  this->P(2,0) = this->P(0,2);
289 
291  }
292 
293  /** \brief calculate density of strain energy
294  *
295  */
297  const NonlinearElasticElement::BlockData block_data,
298  boost::shared_ptr<const NumeredEntFiniteElement> fe_ptr
299 ) {
301 
307  axAngle = -axAngle;
311 
312  voigtStress.resize(6,false);
314  this->eNergy = 0.5*inner_prod(voigtStress,voightStrain);
316  }
317 
321 
324 
325  try {
326 
327  int gg = this->gG; // number of integration point
328  MatrixDouble &phi = (this->commonDataPtr->gradAtGaussPts["POTENTIAL_FIELD"][gg]);
329  normalizedPhi.resize(3,false);
330  double nrm2_phi = sqrt(pow(phi(0,0),2)+pow(phi(0,1),2)+pow(phi(0,2),2));
331  for(int ii = 0;ii<3;ii++) {
332  normalizedPhi[ii] = -phi(0,ii)/nrm2_phi;
333  }
334 
335  axVectorDouble.resize(3,false);
336  const double zVec[3]={ 0.0,0.0,1.0 };
337  axVectorDouble[0] = normalizedPhi[1]*zVec[2]-normalizedPhi[2]*zVec[1];
338  axVectorDouble[1] = normalizedPhi[2]*zVec[0]-normalizedPhi[0]*zVec[2];
339  axVectorDouble[2] = normalizedPhi[0]*zVec[1]-normalizedPhi[1]*zVec[0];
340  double nrm2_ax_vector = norm_2(axVectorDouble);
341  const double eps = 1e-12;
342  if(nrm2_ax_vector<eps) {
343  axVectorDouble[0] = 0;
344  axVectorDouble[1] = 0;
345  axVectorDouble[2] = 1;
346  nrm2_ax_vector = 1;
347  }
348  axAngleDouble = asin(nrm2_ax_vector);
349 
350  } catch (const std::exception& ex) {
351  std::ostringstream ss;
352  ss << "throw in method: " << ex.what() << std::endl;
353  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
354  }
355 
357  }
358 
359 };
360 
361 
363 
365 
368 
369  try {
370 
372  axVector.resize(3,false);
373  axVector[0] = axVectorDouble[0];
374  axVector[1] = axVectorDouble[1];
375  axVector[2] = axVectorDouble[2];
377 
378  } catch (const std::exception& ex) {
379  std::ostringstream ss;
380  ss << "throw in method: " << ex.what() << std::endl;
381  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
382  }
384  }
385 
387  std::map<std::string,std::vector<VectorDouble > > &field_map,
388  std::map<std::string,std::vector<MatrixDouble > > &grad_map
389  ) {
391  int nb_gauss_pts = grad_map["POTENTIAL_FIELD"].size();
392  this->commonDataPtr->gradAtGaussPts["POTENTIAL_FIELD"].resize(nb_gauss_pts);
393  for(int gg = 0;gg<nb_gauss_pts;gg++) {
394  this->commonDataPtr->gradAtGaussPts["POTENTIAL_FIELD"][gg].resize(1,3,false);
395  for(int ii = 0;ii<3;ii++) {
396  this->commonDataPtr->gradAtGaussPts["POTENTIAL_FIELD"][gg](0,ii) =
397  ((grad_map["POTENTIAL_FIELD"])[gg])(0,ii);
398  }
399  }
401  }
402 
403 };
404 
406 
408 
410 
412  int &nb_active_variables
413  ) {
415 
416  try {
417 
419  axVector.resize(3,false);
420  axVector[0] <<= axVectorDouble[0];
421  axVector[1] <<= axVectorDouble[1];
422  axVector[2] <<= axVectorDouble[2];
424  nbActiveVariables0 = nb_active_variables;
425  nb_active_variables += 4;
426 
427  } catch (const std::exception& ex) {
428  std::ostringstream ss;
429  ss << "throw in method: " << ex.what() << std::endl;
430  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
431  }
432 
434  }
435 
437  VectorDouble &active_variables) {
439 
440  try {
441 
442  int shift = nbActiveVariables0; // is a number of elements in F
444  active_variables[shift+0] = axVectorDouble[0];
445  active_variables[shift+1] = axVectorDouble[1];
446  active_variables[shift+2] = axVectorDouble[2];
447  active_variables[shift+3] = axAngleDouble;
448 
449  } catch (const std::exception& ex) {
450  std::ostringstream ss;
451  ss << "throw in method: " << ex.what() << std::endl;
452  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
453  }
454 
456  }
457 
458 };
459 
460 #endif //__SMALLSTRAINTRANVERSLYISOTROPIC_HPP__
static const double eps
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
static double phi
static constexpr double delta
data for calculation heat conductivity and heat capacity elements
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
Implementation of elastic (non-linear) St. Kirchhoff equation.
structure grouping operators and data used for calculation of nonlinear elastic element
virtual MoFEMErrorCode setUserActiveVariables(VectorDouble &active_variables)
Add additional independent variables More complex physical models depend on gradient of defamation an...
virtual MoFEMErrorCode setUserActiveVariables(int &nb_active_variables)
add additional active variables
virtual MoFEMErrorCode getDataOnPostProcessor(std::map< std::string, std::vector< VectorDouble > > &field_map, std::map< std::string, std::vector< MatrixDouble > > &grad_map)
virtual MoFEMErrorCode calculateElasticEnergy(const NonlinearElasticElement::BlockData block_data, boost::shared_ptr< const NumeredEntFiniteElement > fe_ptr)
calculate density of strain energy
virtual MoFEMErrorCode calculateP_PiolaKirchhoffI(const NonlinearElasticElement::BlockData block_data, boost::shared_ptr< const NumeredEntFiniteElement > fe_ptr)
Calculate global stress.
MoFEMErrorCode stressTransformation()
Function to Calculate Stress Transformation Matrix This function computes the stress transformation M...
ublas::symmetric_matrix< TYPE, ublas::upper > localStiffnessMatrix
MoFEMErrorCode strainTransformation()
Function to Calculate Strain Transformation Matrix This function computes the strain transformation ...
MoFEMErrorCode calculateAxisAngleRotationalMatrix()
Function to Calculate the Rotation Matrix at a given axis and angle of rotation.
@ TYPE
Definition: inflate.h:32