Loading [MathJax]/extensions/tex2jax.js
Plugins
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
hooke_impl.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2014-2015: G-CSC, Goethe University Frankfurt
3 * Author: Raphael Prohl
4 *
5 * This file is part of UG4.
6 *
7 * UG4 is free software: you can redistribute it and/or modify it under the
8 * terms of the GNU Lesser General Public License version 3 (as published by the
9 * Free Software Foundation) with the following additional attribution
10 * requirements (according to LGPL/GPL v3 §7):
11 *
12 * (1) The following notice must be displayed in the Appropriate Legal Notices
13 * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 *
15 * (2) The following notice must be displayed at a prominent place in the
16 * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 *
18 * (3) The following bibliography is recommended for citation and must be
19 * preserved in all covered files:
20 * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 * parallel geometric multigrid solver on hierarchically distributed grids.
22 * Computing and visualization in science 16, 4 (2013), 151-164"
23 * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 * flexible software system for simulating pde based models on high performance
25 * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 *
27 * This program is distributed in the hope that it will be useful,
28 * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 * GNU Lesser General Public License for more details.
31 */
32
33#ifndef HOOKE_IMPL_H_
34#define HOOKE_IMPL_H_
35
36#include "common/util/string_table_stream.h"
37
38#include "hooke.h"
39
40#include "voigtian_notation.h"
41
42#define PROFILE_HOOKE
43#ifdef PROFILE_HOOKE
44 #define HOOKE_PROFILE_FUNC() PROFILE_FUNC_GROUP("Hooke")
45 #define HOOKE_PROFILE_BEGIN(name) PROFILE_BEGIN_GROUP(name, "Hooke")
46 #define HOOKE_PROFILE_END() PROFILE_END()
47#else
48 #define HOOKE_PROFILE_FUNC()
49 #define HOOKE_PROFILE_BEGIN(name)
50 #define HOOKE_PROFILE_END()
51#endif
52
53namespace ug{
54namespace SmallStrainMechanics{
55
56template <typename TDomain>
57void
58HookeLaw<TDomain>::
59init()
60{
61 // check, if ElasticityTensor is set
62 if (m_spElastTensorFunct.invalid())
63 UG_THROW("No elasticity tensor set in HookeLaw::init()!");
64
65 base_type::m_bInit = true;
66}
67
68template <typename TDomain>
69void
70HookeLaw<TDomain>::
71stressTensor(MathMatrix<dim,dim>& stressTens, const size_t ip,
72 const MathMatrix<dim, dim>& GradU)
73{
74 HOOKE_PROFILE_BEGIN(HookeLaw_stressTensor);
75
76 // get linearized strain tensor (eps) at ip
77 MathMatrix<dim, dim> strainTens;
78 strainTensor(strainTens, GradU);
79
80 // TODO: replace this with general implementation of TensContractMat
81 for(size_t i = 0; i < (size_t) dim; ++i)
82 for(size_t j = 0; j < (size_t) dim; ++j)
83 {
84 stressTens[i][j] = 0.0;
85
86 for(size_t k = 0; k < (size_t) dim; ++k)
87 for(size_t l = 0; l < (size_t) dim; ++l)
88 stressTens[i][j] += (*m_spElastTensorFunct)[i][j][k][l]
89 * strainTens[k][l];
90
91 }
92
93}
94
95template <typename TDomain>
96inline
97SmartPtr<MathTensor4<TDomain::dim,TDomain::dim,TDomain::dim,TDomain::dim> >
98HookeLaw<TDomain>::
99elasticityTensor()
100{
101 HOOKE_PROFILE_BEGIN(HookeLaw_elasticityTensor);
102 return m_spElastTensorFunct;
103}
104
105template <typename TDomain>
106inline
107void
108HookeLaw<TDomain>::
109strainTensor(MathMatrix<dim,dim>& strainTens, const MathMatrix<dim, dim>& GradU)
110{
111 for(size_t i = 0; i < (size_t) dim; ++i)
112 for(size_t j = 0; j < (size_t) dim; ++j)
113 strainTens[i][j] = 0.5 * (GradU[i][j] + GradU[j][i]);
114}
115
116
117
118template <typename TDomain>
119void
120HookeLaw<TDomain>::
121set_elasticity_tensor_orthotropic_plain_stress_E_G_nu(
122const number E1, const number E2, //const number E3,
123const number G12, //const number G13, const number G23,
124const number v12 //const number v13, const number v23
125)
126{
127 if(dim != 2) UG_THROW("Orthotrope Tensor PLAIN STRESS only for 2 dimensions" );
128
129 MathTensor4<dim,dim,dim,dim> elastTensorFunct;
130
131 // setze alle wWerte auf 0
132 for (size_t i = 0; i < (size_t) dim; ++i)
133 for (size_t j = 0; j < (size_t) dim; ++j)
134 for (size_t k = 0; k < (size_t) dim; ++k)
135 for (size_t l = 0; l < (size_t) dim; ++l)
136 elastTensorFunct[i][j][k][l] = 0.0;
137
138
139 const number v21 = (E2/E1) * v12;
140 //const number v31 = (E3/E1) * v13;
141 //const number v32 = (E3/E2) * v23;
142
143 const number D = 1 - v12*v21;
144
145
146 // Tensor mit Werte fuellen
147 // i j k l
148 elastTensorFunct[0][0][0][0] = E1/D; // C11
149 elastTensorFunct[1][1][1][1] = E2/D; // C22
150
151 elastTensorFunct[0][0][1][1] =
152 elastTensorFunct[1][1][0][0] = E2*v12/D; // = C12
153
154 elastTensorFunct[0][1][0][1] =
155 elastTensorFunct[1][0][0][1] =
156 elastTensorFunct[0][1][1][0] =
157 elastTensorFunct[1][0][1][0] = 2*G12; // C66
158
159 // remembering the elasticity tensor
160 SmartPtr<MathTensor4<dim,dim,dim,dim> > spElastTens(new MathTensor4<dim,dim,dim,dim>(elastTensorFunct));
161 m_spElastTensorFunct = spElastTens;
162
163 StringTableStream sts;
164 std::stringstream ss;
165
166 ss << "Hooke constants Orthotrope (PLAIN STRESS):\n";
167 sts.clear();
168 sts << " E1 = " << E1 << "E2 = " << E2 << "\n";
169 sts << " v12 = " << v12 << "v21 = " << v21 << "\n";
170 sts << " G12 = " << G12 << "\n";
171 ss << sts;
172 m_materialConfiguration = ss.str();
173 UG_LOG("\nset_elasticity_tensor_orthotropic " << ConfigShift(m_materialConfiguration) << "\n");
174}
175
176
177template <typename TDomain>
178void
179HookeLaw<TDomain>::
180set_elasticity_tensor_orthotropic_plain_strain_E_G_nu(
181const number E1, const number E2, const number E3,
182const number G12, const number G13, const number G23,
183const number v12, const number v13, const number v23
184)
185{
186 if(dim != 2) UG_THROW("Orthotrope Tensor PLAIN STRAIN only for 2 dimensions" );
187
188 MathTensor4<dim,dim,dim,dim> elastTensorFunct;
189
190 // setze alle wWerte auf 0
191 for (size_t i = 0; i < (size_t) dim; ++i)
192 for (size_t j = 0; j < (size_t) dim; ++j)
193 for (size_t k = 0; k < (size_t) dim; ++k)
194 for (size_t l = 0; l < (size_t) dim; ++l)
195 elastTensorFunct[i][j][k][l] = 0.0;
196
197
198 const number v21 = (E2/E1) * v12;
199 const number v31 = (E3/E1) * v13;
200 const number v32 = (E3/E2) * v23;
201
202 const number D = 1 - v12*v21 - v13*v31 - v23*v32 - 2*v12*v23*v31;
203
204
205 // Tensor mit Werte fuellen
206 // i j k l
207 elastTensorFunct[0][0][0][0] = E1*(1-v23*v32)/D; // C11
208 elastTensorFunct[1][1][1][1] = E2*(1-v13*v31)/D; // C22
209
210 elastTensorFunct[0][0][1][1] =
211 elastTensorFunct[1][1][0][0] = E2*(v13*v32+v12)/D; // = C12
212
213 elastTensorFunct[0][1][0][1] =
214 elastTensorFunct[1][0][0][1] =
215 elastTensorFunct[0][1][1][0] =
216 elastTensorFunct[1][0][1][0] = 2*G12; // C66
217
218 // remembering the elasticity tensor
219 SmartPtr<MathTensor4<dim,dim,dim,dim> > spElastTens(new MathTensor4<dim,dim,dim,dim>(elastTensorFunct));
220 m_spElastTensorFunct = spElastTens;
221
222 StringTableStream sts;
223 std::stringstream ss;
224
225 ss << "Hooke constants Orthotrope (PLAIN STRAIN):\n";
226 sts.clear();
227 sts << " E1 = " << E1 << "E2 = " << E2 << "E3 = " << E3 << "\n";
228 sts << " v12 = " << v12 << "v23 = " << v23 << "v13 = " << v13 << "\n";
229 sts << " G12 = " << G12 << "G23 = " << G23 << "G13 = " << G13 << "\n";
230 ss << sts;
231 m_materialConfiguration = ss.str();
232 UG_LOG("\nset_elasticity_tensor_orthotropic " << ConfigShift(m_materialConfiguration) << "\n");
233}
234
235
236
237template <typename TDomain>
238void
239HookeLaw<TDomain>::
240set_elasticity_tensor_orthotropic_E_G_nu(
241const number E1, const number E2, const number E3,
242const number G12, const number G13, const number G23,
243const number v12, const number v13, const number v23
244)
245{
246 if(dim != 3) UG_THROW("Orthotrope Tensor only for 3 dimensions" );
247
248 MathTensor4<dim,dim,dim,dim> elastTensorFunct;
249
250 VoigtianMatrix<TDomain> voigt;
251 const number v21 = (E2/E1) * v12;
252 const number v31 = (E3/E1) * v13;
253 const number v32 = (E3/E2) * v23;
254
255 const number D = 1 - v12*v21 - v13*v31 - v23*v32 - 2*v12*v23*v31;
256
257 const number C11 = E1*(1-v23*v32)/D;
258 const number C22 = E2*(1-v13*v31)/D;
259 const number C33 = E3*(1-v12*v21)/D;
260
261 const number C12 = E2*(v13*v32+v12)/D;
262 const number C13 = E3*(v12*v23+v13)/D;
263 const number C23 = E3*(v21*v13+v23)/D;
264
265 const number C44 = 2*G23;
266 const number C55 = 2*G13;
267 const number C66 = 2*G12;
268
269
270 voigt.set_orthotropic(C11, C12, C13, C22, C23, C33, C44, C55, C66);
271 voigt.copy_to_tensor(elastTensorFunct);
272
273 // remembering the elasticity tensor
274 SmartPtr<MathTensor4<dim,dim,dim,dim> > spElastTens(new MathTensor4<dim,dim,dim,dim>(elastTensorFunct));
275 m_spElastTensorFunct = spElastTens;
276
277 StringTableStream sts;
278 std::stringstream ss;
279
280 ss << "Hooke constants Orthotrope:\n";
281 sts.clear();
282 sts << " E1 = " << E1 << "E2 = " << E2 << "E3 = " << E3 << "\n";
283 sts << " v12 = " << v12 << "v23 = " << v23 << "v13 = " << v13 << "\n";
284 sts << " G12 = " << G12 << "G23 = " << G23 << "G13 = " << G13 << "\n";
285 ss << sts;
286 m_materialConfiguration = ss.str();
287 UG_LOG("\nset_elasticity_tensor_orthotropic " << ConfigShift(m_materialConfiguration) << "\n");
288}
289
290
291template <typename TDomain>
292void
293HookeLaw<TDomain>::
294set_elasticity_tensor_orthotropic(
295const number C11, const number C12, const number C13,
296const number C22, const number C23, const number C33,
297const number C44, const number C55, const number C66)
298{
299 if(dim != 3) UG_THROW("Orthotrope Tensor only for 3 dimensions" );
300
301 MathTensor4<dim,dim,dim,dim> elastTensorFunct;
302
303 VoigtianMatrix<TDomain> voigt;
304 voigt.set_orthotropic(C11, C12, C13, C22, C23, C33, C44, C55, C66);
305 voigt.copy_to_tensor(elastTensorFunct);
306
307 // remembering the elasticity tensor
308 SmartPtr<MathTensor4<dim,dim,dim,dim> > spElastTens(new MathTensor4<dim,dim,dim,dim>(elastTensorFunct));
309 m_spElastTensorFunct = spElastTens;
310
311 DenseMatrix<FixedArray2<double, 3, 3> > mat;
312 mat(0,0) = C11; mat(0,1) = C12; mat(0,2) = C13;
313 mat(1,0) = C12; mat(1,1) = C22; mat(1,2) = C23;
314 mat(2,0) = C13; mat(2,1) = C23; mat(2,2) = C33;
315
316 Invert(mat);
317 //UG_LOG("S = " mat << "\n");
318 number E1 = 1/mat(0,0), E2 = 1/mat(1,1), E3 = 1/mat(2,2);
319 number v12 = -mat(1, 0)/E1;
320 //number v21 = -mat(0, 1)/E2;
321 number v23 = -mat(2, 1)/E2;
322 //number v32 = -mat(1, 2)/E3;
323 number v13 = -mat(2, 0)/E1;
324 //number v31 = -mat(0, 2)/E3;
325 number G23 = C44;
326 number G13 = C55;
327 number G12 = C66;
328
329 StringTableStream sts;
330 std::stringstream ss;
331 ss << "Orthrope Elasticity Tensor: \n";
332 sts.clear();
333 sts << C11 << C12 << C13 << "\n";
334 sts << C12 << C22 << C23 << "\n";
335 sts << C13 << C23 << C33 << "\n";
336 sts << sts.empty_col(3) << C44 << "\n";
337 sts << sts.empty_col(4) << C55 << "\n";
338 sts << sts.empty_col(5) << C66 << "\n";
339 ss << sts;
340
341 ss << "Hooke constants Orthotrope:\n";
342 sts.clear();
343 sts << " E1 = " << E1 << "E2 = " << E2 << "E3 = " << E3 << "\n";
344 sts << " v12 = " << v12 << "v23 = " << v23 << "v13 = " << v13 << "\n";
345 sts << " G12 = " << G12 << "G23 = " << G23 << "G13 = " << G13 << "\n";
346 ss << sts;
347 m_materialConfiguration = ss.str();
348 UG_LOG("\nset_elasticity_tensor_orthotropic " << ConfigShift(m_materialConfiguration) << "\n");
349}
350
351
352template <typename TDomain>
353void
354HookeLaw<TDomain>::
355set_hooke_elasticity_tensor_plain_stress_E_nu(const number E, const number nu)
356{
357 if(dim != 2) UG_THROW("Plain Stress Tensor only for 2 dimensions" );
358
359 MathTensor4<dim,dim,dim,dim> elastTensorFunct;
360
361 // setze alle wWerte auf 0
362 for (size_t i = 0; i < (size_t) dim; ++i)
363 for (size_t j = 0; j < (size_t) dim; ++j)
364 for (size_t k = 0; k < (size_t) dim; ++k)
365 for (size_t l = 0; l < (size_t) dim; ++l)
366 elastTensorFunct[i][j][k][l] = 0.0;
367
368 // Tensor mit Werte fuellen
369 // i j k l
370 elastTensorFunct[0][0][0][0] = E/(1-nu*nu); // C11
371 elastTensorFunct[1][1][1][1] = E/(1-nu*nu); // C22
372
373 elastTensorFunct[0][0][1][1] =
374 elastTensorFunct[1][1][0][0] = (E*nu)/(1-nu*nu); // = C12 = C21
375
376 elastTensorFunct[0][1][0][1] =
377 elastTensorFunct[1][0][0][1] =
378 elastTensorFunct[0][1][1][0] =
379 elastTensorFunct[1][0][1][0] = E*(1-nu)/(1-nu*nu); // C33
380
381 // remembering the elasticity tensor
382 SmartPtr<MathTensor4<dim,dim,dim,dim> > spElastTens(new MathTensor4<dim,dim,dim,dim>(elastTensorFunct));
383 m_spElastTensorFunct = spElastTens;
384
385 std::stringstream ss;
386 ss << "Hooke Elasticity Isotropic (PLAIN STRESS): \n";
387 ss << " young modulus (Elastizitaetsmodul): " << E << "\n";
388 ss << " poisson ratio (Querkontraktionszahl) v: " << nu << "\n";
389 ss << " Elasticity Tensor = " << elastTensorFunct << "\n";
390 m_materialConfiguration = ss.str();
391}
392
393template <typename TDomain>
394void
395HookeLaw<TDomain>::
396set_hooke_elasticity_tensor_plain_strain_E_nu(const number E, const number nu)
397{
398 if(dim != 2) UG_THROW("Plain Strain Tensor only for 2 dimensions" );
399
400 set_hooke_elasticity_tensor_E_nu(E, nu);
401}
402
403
404template <typename TDomain>
405void
406HookeLaw<TDomain>::
407set_hooke_elasticity_tensor_E_nu(const number E, const number nu)
408{
409 number lambda = (E*nu) / ((1+nu)*(1-2*nu));
410 number mu = E/(2*(1+nu));
411 set_hooke_elasticity_tensor(lambda, mu);
412}
413
414template <typename TDomain>
415void
416HookeLaw<TDomain>::
417set_hooke_elasticity_tensor(const number lambda, const number mu)
418{
419 // sets the 'Hooke'-elasticity tensor for isotropic materials
420 // in 2D this tensor formulation corresponds to the
421 // plane strain assumption for Hookes`s law
422
423 MathTensor4<dim,dim,dim,dim> elastTensorFunct;
424
425 // filling the constant elasticity tensor
426 for (size_t i = 0; i < (size_t) dim; ++i)
427 for (size_t j = 0; j < (size_t) dim; ++j)
428 for (size_t k = 0; k < (size_t) dim; ++k)
429 for (size_t l = 0; l < (size_t) dim; ++l)
430 {
431 elastTensorFunct[i][j][k][l] = 0.0;
432
433 if ((i == j) && (k == l))
434 elastTensorFunct[i][j][k][l] += lambda;
435
436 if ((i == k) && (j == l))
437 elastTensorFunct[i][j][k][l] += mu;
438
439 if ((i == l) && (j == k))
440 elastTensorFunct[i][j][k][l] += mu;
441
442 } //end (l)
443
444 // remembering the elasticity tensor
445 SmartPtr<MathTensor4<dim,dim,dim,dim> > spElastTens(new MathTensor4<dim,dim,dim,dim>(elastTensorFunct));
446 m_spElastTensorFunct = spElastTens;
447
448 number E = mu * (3.0 * lambda + 2.0 * mu)/(lambda + mu);
449 number v = 0.5 * lambda/(lambda + mu);
450 number kappa = lambda + 2.0/3.0 * mu;
451
452 std::stringstream ss;
453 ss << "Hooke Elasticity Isotropic: \n";
454 ss << " Lame`s first constant lambda: " << lambda << "\n";
455 ss << " Lame`s second constant mue (sometimes 'G', shear modulus) (Schubmodul): " << mu << "\n";
456 ss << " This setting equals: \n";
457 ss << " young modulus (Elastizitaetsmodul): " << E << "\n";
458 ss << " poisson ratio (Querkontraktionszahl) v: " << v << "\n";
459 ss << " bulk modulus (Kompressionsmodul): " << kappa << "\n";
460 ss << " Elasticity Tensor = " << elastTensorFunct << "\n";
461 m_materialConfiguration = ss.str();
462}
463
464}// end of namespace SmallStrainMechanics
465}// end of namespace ug
466
467#endif /* HOOKE_IMPL_H_ */
ug::SmallStrainMechanics::VoigtianMatrix
Definition voigtian_notation.h:40
ug::SmallStrainMechanics::VoigtianMatrix::set_orthotropic
void set_orthotropic(const number C11, const number C12, const number C13, const number C22, const number C23, const number C33, const number C44, const number C55, const number C66)
Definition voigtian_notation.h:46
ug::SmallStrainMechanics::VoigtianMatrix::copy_to_tensor
void copy_to_tensor(MathTensor4< dim, dim, dim, dim > &elastTensorFunct)
Definition voigtian_notation.h:77
ug::Invert
bool Invert(DenseMatrix< FixedArray2< double, 1, 1 > > &mat)
ug::SmallStrainMechanics::HookeLaw::set_elasticity_tensor_orthotropic_E_G_nu
void set_elasticity_tensor_orthotropic_E_G_nu(const number E1, const number E2, const number E3, const number G12, const number G13, const number G23, const number v12, const number v13, const number v23)
Definition hooke_impl.h:240
ug::SmallStrainMechanics::HookeLaw::elasticityTensor
virtual SmartPtr< MathTensor4< TDomain::dim, TDomain::dim, TDomain::dim, TDomain::dim > > elasticityTensor()
computes the elasticity tensor; commonly denoted by C
Definition hooke_impl.h:99
ug::SmallStrainMechanics::HookeLaw::set_hooke_elasticity_tensor_plain_strain_E_nu
void set_hooke_elasticity_tensor_plain_strain_E_nu(const number E, const number nu)
Definition hooke_impl.h:396
ug::SmallStrainMechanics::HookeLaw::strainTensor
void strainTensor(MathMatrix< dim, dim > &strainTens, const MathMatrix< dim, dim > &GradU)
Definition hooke_impl.h:109
ug::SmallStrainMechanics::HookeLaw::set_hooke_elasticity_tensor_E_nu
void set_hooke_elasticity_tensor_E_nu(const number E, const number nu)
Definition hooke_impl.h:407
ug::SmallStrainMechanics::HookeLaw::set_hooke_elasticity_tensor
void set_hooke_elasticity_tensor(const number lambda, const number mu)
Definition hooke_impl.h:417
ug::SmallStrainMechanics::HookeLaw::set_elasticity_tensor_orthotropic
void set_elasticity_tensor_orthotropic(const number C11, const number C12, const number C13, const number C22, const number C23, const number C33, const number C44, const number C55, const number C66)
Definition hooke_impl.h:294
ug::SmallStrainMechanics::HookeLaw::set_elasticity_tensor_orthotropic_plain_strain_E_G_nu
void set_elasticity_tensor_orthotropic_plain_strain_E_G_nu(const number E1, const number E2, const number E3, const number G12, const number G13, const number G23, const number v12, const number v13, const number v23)
Definition hooke_impl.h:180
ug::SmallStrainMechanics::HookeLaw::stressTensor
virtual void stressTensor(MathMatrix< dim, dim > &stressTens, const size_t ip, const MathMatrix< dim, dim > &GradU)
computes the cauchy stress tensor sigma at an integration point 'ip'
Definition hooke_impl.h:71
ug::SmallStrainMechanics::HookeLaw::set_hooke_elasticity_tensor_plain_stress_E_nu
void set_hooke_elasticity_tensor_plain_stress_E_nu(const number E, const number nu)
Definition hooke_impl.h:355
ug::SmallStrainMechanics::HookeLaw::set_elasticity_tensor_orthotropic_plain_stress_E_G_nu
void set_elasticity_tensor_orthotropic_plain_stress_E_G_nu(const number E1, const number E2, const number G12, const number v12)
Definition hooke_impl.h:121
ug::SmallStrainMechanics::HookeLaw::init
virtual void init()
Definition hooke_impl.h:59
ug::StringTableStream::clear
void clear()
ug::StringTableStream::empty_col
RepeatedCol empty_col(size_t i)
UG_THROW
#define UG_THROW(msg)
UG_LOG
#define UG_LOG(msg)
number
double number
HOOKE_PROFILE_BEGIN
#define HOOKE_PROFILE_BEGIN(name)
Definition hooke_impl.h:45
ug::ConfigShift
UG_API std::string ConfigShift(std::string s)
string_table_stream.h
voigtian_notation.h