v0.15.0
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;
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()
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
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.