v0.13.1
Loading...
Searching...
No Matches
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
7
8#ifndef __SMALLSTRAINTRANVERSLYISOTROPIC_HPP__
9#define __SMALLSTRAINTRANVERSLYISOTROPIC_HPP__
10
11#ifndef WITH_ADOL_C
12 #error "MoFEM need to be compiled with ADOL-C"
13#endif
14
15/** \brief Hook equation
16 * \ingroup nonlinear_elastic_elem
17 */
18template<typename TYPE>
20
22
23
24
25 ublas::matrix<TYPE> sTrain;
26 ublas::vector<TYPE> voightStrain;
27 TYPE tR;
28
29 MoFEMErrorCode calculateStrain() {
31 sTrain.resize(3,3,false);
32 noalias(sTrain) = this->F;
33 for(int dd = 0;dd<3;dd++) {
34 sTrain(dd,dd) -= 1;
35 }
36 sTrain += trans(sTrain);
37 sTrain *= 0.5;
38 voightStrain.resize(6,false);
39 voightStrain[0] = sTrain(0,0);
40 voightStrain[1] = sTrain(1,1);
41 voightStrain[2] = sTrain(2,2);
42 voightStrain[3] = 2*sTrain(0,1);
43 voightStrain[4] = 2*sTrain(1,2);
44 voightStrain[5] = 2*sTrain(2,0);
46 }
47
48 double nu_p, nu_pz, E_p, E_z, G_zp;
49 ublas::symmetric_matrix<TYPE,ublas::upper> localStiffnessMatrix;
52 double nu_zp=(nu_pz*E_z)/E_p;
53 double delta=((1+nu_p)*(1-nu_p-(2*nu_pz*nu_zp)))/(E_p*E_p*E_z);
54
55 localStiffnessMatrix.resize(6);
59
62
63 localStiffnessMatrix(3,3)=E_p/(2*(1+nu_p));
65
67 }
68
69 ublas::matrix<TYPE> aARotMat;
70 ublas::vector<TYPE> axVector;
71 TYPE axAngle;
72
73 /**
74 * \brief Function to Calculate the Rotation Matrix at a given axis and angle of rotation
75
76 * This function computes the rotational matrix for a given axis of rotation
77 * and angle of rotation about that angle <br>
78
79 *\param axVector A vector representing the axis of rotation
80 *\param axAngle Angle of rotation along the axis (in radians)
81 */
84
85 aARotMat.resize(3,3,false);
86 aARotMat.clear();
87
88 TYPE norm = sqrt(pow(axVector[0],2) + pow(axVector[1],2) + pow(axVector[2],2));
89
90 aARotMat(0,0) = 1-((1-cos(axAngle))*(pow(axVector[1],2)+pow(axVector[2],2))/pow(norm,2));
91 aARotMat(1,1) = 1-((1-cos(axAngle))*(pow(axVector[0],2)+pow(axVector[2],2))/pow(norm,2));
92 aARotMat(2,2) = 1-((1-cos(axAngle))*(pow(axVector[0],2)+pow(axVector[1],2))/pow(norm,2));
93
94 aARotMat(0,1) = ((1-cos(axAngle))*axVector[0]*axVector[1]-norm*axVector[2]*sin(axAngle))/pow(norm,2);
95 aARotMat(1,0) = ((1-cos(axAngle))*axVector[0]*axVector[1]+norm*axVector[2]*sin(axAngle))/pow(norm,2);
96
97 aARotMat(0,2) = ((1-cos(axAngle))*axVector[0]*axVector[2]+norm*axVector[1]*sin(axAngle))/pow(norm,2);
98 aARotMat(2,0) = ((1-cos(axAngle))*axVector[0]*axVector[2]-norm*axVector[1]*sin(axAngle))/pow(norm,2);
99
100 aARotMat(1,2) = ((1-cos(axAngle))*axVector[1]*axVector[2]-norm*axVector[0]*sin(axAngle))/pow(norm,2);
101 aARotMat(2,1) = ((1-cos(axAngle))*axVector[1]*axVector[2]+norm*axVector[0]*sin(axAngle))/pow(norm,2);
102
104 }
105
106 ublas::matrix<TYPE> stressRotMat;
107
108 /**
109 * \brief Function to Calculate Stress Transformation Matrix
110 * This function computes the stress transformation Matrix at a give axis and angle of rotation <br>
111 * One can also output the axis/angle rotational Matrix
112 */
113 MoFEMErrorCode stressTransformation() {
115
116 stressRotMat.resize(6,6,false);
117 stressRotMat.clear();
118
119 stressRotMat(0, 0) = aARotMat(0,0) * aARotMat(0,0);
120 stressRotMat(0, 1) = aARotMat(1,0) * aARotMat(1,0);
121 stressRotMat(0, 2) = aARotMat(2,0) * aARotMat(2,0);
122 stressRotMat(0, 3) = 2.0 * aARotMat(1,0) * aARotMat(0,0);
123 stressRotMat(0, 4) = 2.0 * aARotMat(2,0) * aARotMat(1,0);
124 stressRotMat(0, 5) = 2.0 * aARotMat(0,0) * aARotMat(2,0);
125
126 stressRotMat(1, 0) = aARotMat(0,1) * aARotMat(0,1);
127 stressRotMat(1, 1) = aARotMat(1,1) * aARotMat(1,1);
128 stressRotMat(1, 2) = aARotMat(2,1) * aARotMat(2,1);
129 stressRotMat(1, 3) = 2.0 * aARotMat(1,1) * aARotMat(0,1);
130 stressRotMat(1, 4) = 2.0 * aARotMat(2,1) * aARotMat(1,1);
131 stressRotMat(1, 5) = 2.0 * aARotMat(0,1) * aARotMat(2,1);
132
133 stressRotMat(2, 0) = aARotMat(0,2) * aARotMat(0,2);
134 stressRotMat(2, 1) = aARotMat(1,2) * aARotMat(1,2);
135 stressRotMat(2, 2) = aARotMat(2,2) * aARotMat(2,2);
136 stressRotMat(2, 3) = 2.0 * aARotMat(1,2) * aARotMat(0,2);
137 stressRotMat(2, 4) = 2.0 * aARotMat(2,2) * aARotMat(1,2);
138 stressRotMat(2, 5) = 2.0 * aARotMat(0,2) * aARotMat(2,2);
139
140 stressRotMat(3, 0) = aARotMat(0,1) * aARotMat(0,0);
141 stressRotMat(3, 1) = aARotMat(1,1) * aARotMat(1,0);
142 stressRotMat(3, 2) = aARotMat(2,1) * aARotMat(2,0);
143 stressRotMat(3, 3) = ( aARotMat(1,1) * aARotMat(0,0) + aARotMat(0,1) * aARotMat(1,0) );
144 stressRotMat(3, 4) = ( aARotMat(2,1) * aARotMat(1,0) + aARotMat(1,1) * aARotMat(2,0) );
145 stressRotMat(3, 5) = ( aARotMat(0,1) * aARotMat(2,0) + aARotMat(2,1) * aARotMat(0,0) );
146
147 stressRotMat(4, 0) = aARotMat(0,2) * aARotMat(0,1);
148 stressRotMat(4, 1) = aARotMat(1,2) * aARotMat(1,1);
149 stressRotMat(4, 2) = aARotMat(2,2) * aARotMat(2,1);
150 stressRotMat(4, 3) = ( aARotMat(1,2) * aARotMat(0,1) + aARotMat(0,2) * aARotMat(1,1) );
151 stressRotMat(4, 4) = ( aARotMat(2,2) * aARotMat(1,1) + aARotMat(1,2) * aARotMat(2,1) );
152 stressRotMat(4, 5) = ( aARotMat(0,2) * aARotMat(2,1) + aARotMat(2,2) * aARotMat(0,1) );
153
154 stressRotMat(5, 0) = aARotMat(0,0) * aARotMat(0,2);
155 stressRotMat(5, 1) = aARotMat(1,0) * aARotMat(1,2);
156 stressRotMat(5, 2) = aARotMat(2,0) * aARotMat(2,2);
157 stressRotMat(5, 3) = ( aARotMat(1,0) * aARotMat(0,2) + aARotMat(0,0) * aARotMat(1,2) );
158 stressRotMat(5, 4) = ( aARotMat(2,0) * aARotMat(1,2) + aARotMat(1,0) * aARotMat(2,2) );
159 stressRotMat(5, 5) = ( aARotMat(0,0) * aARotMat(2,2) + aARotMat(2,0) * aARotMat(0,2) );
160
162 }
163
164 ublas::matrix<TYPE> strainRotMat;
165
166 /**
167 * \brief Function to Calculate Strain Transformation Matrix<br>
168 * This function computes the strain transformation Matrix at a give axis and angle of rotation <br>
169 * One can also output the axis/angle rotational Matrix
170 */
171 MoFEMErrorCode strainTransformation() {
173
174 strainRotMat.resize(6,6,false);
175 strainRotMat.clear();
176
177 strainRotMat(0, 0) = aARotMat(0,0) * aARotMat(0,0);
178 strainRotMat(0, 1) = aARotMat(1,0) * aARotMat(1,0);
179 strainRotMat(0, 2) = aARotMat(2,0) * aARotMat(2,0);
180 strainRotMat(0, 3) = aARotMat(1,0) * aARotMat(0,0);
181 strainRotMat(0, 4) = aARotMat(2,0) * aARotMat(1,0);
182 strainRotMat(0, 5) = aARotMat(0,0) * aARotMat(2,0);
183
184 strainRotMat(1, 0) = aARotMat(0,1) * aARotMat(0,1);
185 strainRotMat(1, 1) = aARotMat(1,1) * aARotMat(1,1);
186 strainRotMat(1, 2) = aARotMat(2,1) * aARotMat(2,1);
187 strainRotMat(1, 3) = aARotMat(1,1) * aARotMat(0,1);
188 strainRotMat(1, 4) = aARotMat(2,1) * aARotMat(1,1);
189 strainRotMat(1, 5) = aARotMat(0,1) * aARotMat(2,1);
190
191 strainRotMat(2, 0) = aARotMat(0,2) * aARotMat(0,2);
192 strainRotMat(2, 1) = aARotMat(1,2) * aARotMat(1,2);
193 strainRotMat(2, 2) = aARotMat(2,2) * aARotMat(2,2);
194 strainRotMat(2, 3) = aARotMat(1,2) * aARotMat(0,2);
195 strainRotMat(2, 4) = aARotMat(2,2) * aARotMat(1,2);
196 strainRotMat(2, 5) = aARotMat(0,2) * aARotMat(2,2);
197
198 strainRotMat(3, 0) = 2.0 * aARotMat(0,1) * aARotMat(0,0);
199 strainRotMat(3, 1) = 2.0 * aARotMat(1,1) * aARotMat(1,0);
200 strainRotMat(3, 2) = 2.0 * aARotMat(2,1) * aARotMat(2,0);
201 strainRotMat(3, 3) = ( aARotMat(1,1) * aARotMat(0,0) + aARotMat(0,1) * aARotMat(1,0) );
202 strainRotMat(3, 4) = ( aARotMat(2,1) * aARotMat(1,0) + aARotMat(1,1) * aARotMat(2,0) );
203 strainRotMat(3, 5) = ( aARotMat(0,1) * aARotMat(2,0) + aARotMat(2,1) * aARotMat(0,0) );
204
205 strainRotMat(4, 0) = 2.0 * aARotMat(0,2) * aARotMat(0,1);
206 strainRotMat(4, 1) = 2.0 * aARotMat(1,2) * aARotMat(1,1);
207 strainRotMat(4, 2) = 2.0 * aARotMat(2,2) * aARotMat(2,1);
208 strainRotMat(4, 3) = ( aARotMat(1,2) * aARotMat(0,1) + aARotMat(0,2) * aARotMat(1,1) );
209 strainRotMat(4, 4) = ( aARotMat(2,2) * aARotMat(1,1) + aARotMat(1,2) * aARotMat(2,1) );
210 strainRotMat(4, 5) = ( aARotMat(0,2) * aARotMat(2,1) + aARotMat(2,2) * aARotMat(0,1) );
211
212 strainRotMat(5, 0) = 2.0 * aARotMat(0,0) * aARotMat(0,2);
213 strainRotMat(5, 1) = 2.0 * aARotMat(1,0) * aARotMat(1,2);
214 strainRotMat(5, 2) = 2.0 * aARotMat(2,0) * aARotMat(2,2);
215 strainRotMat(5, 3) = ( aARotMat(1,0) * aARotMat(0,2) + aARotMat(0,0) * aARotMat(1,2) );
216 strainRotMat(5, 4) = ( aARotMat(2,0) * aARotMat(1,2) + aARotMat(1,0) * aARotMat(2,2) );
217 strainRotMat(5, 5) = ( aARotMat(0,0) * aARotMat(2,2) + aARotMat(2,0) * aARotMat(0,2) );
218
220 }
221
222 ublas::matrix<TYPE> dR;
223 ublas::matrix<TYPE> globalStiffnessMatrix;
224
227
228 dR.resize(6,6,false);
229 noalias(dR) = prod(localStiffnessMatrix,strainRotMat);
230 globalStiffnessMatrix.resize(6,6,false);
231 noalias(globalStiffnessMatrix) = prod(stressRotMat,dR);
232
234 }
235
236 virtual MoFEMErrorCode calculateAngles() {
239 }
240
241 ublas::vector<TYPE> voigtStress;
242
243 /** \brief Calculate global stress
244
245 This is small strain approach, i.e. Piola stress
246 is like a Cauchy stress, since configurations are notation
247 distinguished.
248
249 */
250 virtual MoFEMErrorCode calculateP_PiolaKirchhoffI(
251 const NonlinearElasticElement::BlockData block_data,
252 boost::shared_ptr<const NumeredEntFiniteElement> fe_ptr
253 ) {
260 axAngle = -axAngle;
264
265 voigtStress.resize(6,false);
267 this->P.resize(3,3,false);
268 this->P(0,0) = voigtStress[0];
269 this->P(1,1) = voigtStress[1];
270 this->P(2,2) = voigtStress[2];
271 this->P(0,1) = voigtStress[3];
272 this->P(1,2) = voigtStress[4];
273 this->P(0,2) = voigtStress[5];
274 this->P(1,0) = this->P(0,1);
275 this->P(2,1) = this->P(1,2);
276 this->P(2,0) = this->P(0,2);
277
279 }
280
281 /** \brief calculate density of strain energy
282 *
283 */
284 virtual MoFEMErrorCode calculateElasticEnergy(
285 const NonlinearElasticElement::BlockData block_data,
286 boost::shared_ptr<const NumeredEntFiniteElement> fe_ptr
287) {
289
295 axAngle = -axAngle;
299
300 voigtStress.resize(6,false);
302 this->eNergy = 0.5*inner_prod(voigtStress,voightStrain);
304 }
305
306 VectorDouble normalizedPhi;
307 VectorDouble axVectorDouble;
309
310 MoFEMErrorCode calculateFibreAngles() {
312
313 try {
314
315 int gg = this->gG; // number of integration point
316 MatrixDouble &phi = (this->commonDataPtr->gradAtGaussPts["POTENTIAL_FIELD"][gg]);
317 normalizedPhi.resize(3,false);
318 double nrm2_phi = sqrt(pow(phi(0,0),2)+pow(phi(0,1),2)+pow(phi(0,2),2));
319 for(int ii = 0;ii<3;ii++) {
320 normalizedPhi[ii] = -phi(0,ii)/nrm2_phi;
321 }
322
323 axVectorDouble.resize(3,false);
324 const double zVec[3]={ 0.0,0.0,1.0 };
325 axVectorDouble[0] = normalizedPhi[1]*zVec[2]-normalizedPhi[2]*zVec[1];
326 axVectorDouble[1] = normalizedPhi[2]*zVec[0]-normalizedPhi[0]*zVec[2];
327 axVectorDouble[2] = normalizedPhi[0]*zVec[1]-normalizedPhi[1]*zVec[0];
328 double nrm2_ax_vector = norm_2(axVectorDouble);
329 const double eps = 1e-12;
330 if(nrm2_ax_vector<eps) {
331 axVectorDouble[0] = 0;
332 axVectorDouble[1] = 0;
333 axVectorDouble[2] = 1;
334 nrm2_ax_vector = 1;
335 }
336 axAngleDouble = asin(nrm2_ax_vector);
337
338 } catch (const std::exception& ex) {
339 std::ostringstream ss;
340 ss << "throw in method: " << ex.what() << std::endl;
341 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
342 }
343
345 }
346
347};
348
349
351
353
354 virtual MoFEMErrorCode calculateAngles() {
356
357 try {
358
360 axVector.resize(3,false);
361 axVector[0] = axVectorDouble[0];
362 axVector[1] = axVectorDouble[1];
363 axVector[2] = axVectorDouble[2];
365
366 } catch (const std::exception& ex) {
367 std::ostringstream ss;
368 ss << "throw in method: " << ex.what() << std::endl;
369 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
370 }
372 }
373
374 virtual MoFEMErrorCode getDataOnPostProcessor(
375 std::map<std::string,std::vector<VectorDouble > > &field_map,
376 std::map<std::string,std::vector<MatrixDouble > > &grad_map
377 ) {
379 int nb_gauss_pts = grad_map["POTENTIAL_FIELD"].size();
380 this->commonDataPtr->gradAtGaussPts["POTENTIAL_FIELD"].resize(nb_gauss_pts);
381 for(int gg = 0;gg<nb_gauss_pts;gg++) {
382 this->commonDataPtr->gradAtGaussPts["POTENTIAL_FIELD"][gg].resize(1,3,false);
383 for(int ii = 0;ii<3;ii++) {
384 this->commonDataPtr->gradAtGaussPts["POTENTIAL_FIELD"][gg](0,ii) =
385 ((grad_map["POTENTIAL_FIELD"])[gg])(0,ii);
386 }
387 }
389 }
390
391};
392
394
396
398
399 virtual MoFEMErrorCode setUserActiveVariables(
400 int &nb_active_variables
401 ) {
403
404 try {
405
407 axVector.resize(3,false);
408 axVector[0] <<= axVectorDouble[0];
409 axVector[1] <<= axVectorDouble[1];
410 axVector[2] <<= axVectorDouble[2];
412 nbActiveVariables0 = nb_active_variables;
413 nb_active_variables += 4;
414
415 } catch (const std::exception& ex) {
416 std::ostringstream ss;
417 ss << "throw in method: " << ex.what() << std::endl;
418 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
419 }
420
422 }
423
424 virtual MoFEMErrorCode setUserActiveVariables(
425 VectorDouble &active_variables) {
427
428 try {
429
430 int shift = nbActiveVariables0; // is a number of elements in F
432 active_variables[shift+0] = axVectorDouble[0];
433 active_variables[shift+1] = axVectorDouble[1];
434 active_variables[shift+2] = axVectorDouble[2];
435 active_variables[shift+3] = axAngleDouble;
436
437 } catch (const std::exception& ex) {
438 std::ostringstream ss;
439 ss << "throw in method: " << ex.what() << std::endl;
440 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
441 }
442
444 }
445
446};
447
448#endif //__SMALLSTRAINTRANVERSLYISOTROPIC_HPP__
static PetscErrorCode ierr
static const double eps
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
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)
Do operations when pre-process.
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.