1 /*!
2 * \file   EllipticCreep.hxx
3 * \brief  this file implements the EllipticCreep Behaviour.
4 *         File generated by tfel version 2.0
5 * \author Maxime Salvo / Thomas Helfer
6 * \date   9 Octobre 2013
7  */
8 
9 #ifndef LIB_TFELMATERIAL_ELLIPTICCREEP_HXX_
10 #define LIB_TFELMATERIAL_ELLIPTICCREEP_HXX_
11 
12 #include<iostream>
13 #include<algorithm>
14 #include<string>
15 
16 #include"TFEL/Config/TFELConfig.hxx"
17 #include"TFEL/Config/TFELTypes.hxx"
18 #include"TFEL/Metaprogramming/StaticAssert.hxx"
19 #include"TFEL/TypeTraits/IsFundamentalNumericType.hxx"
20 #include"TFEL/TypeTraits/IsReal.hxx"
21 #include"TFEL/Material/MechanicalBehaviour.hxx"
22 #include"TFEL/Material/MechanicalBehaviourTraits.hxx"
23 #include"TFEL/Material/OutOfBoundsPolicy.hxx"
24 #include"TFEL/Material/BoundsCheck.hxx"
25 #include"TFEL/Material/EllipticCreepBehaviourData.hxx"
26 #include"TFEL/Material/EllipticCreepIntegrationData.hxx"
27 
28 #include<limits>
29 #include<algorithm>
30 
31 #include"TFEL/Math/st2tost2.hxx"
32 #include"TFEL/Math/tmatrix.hxx"
33 #include"TFEL/Math/tvector.hxx"
34 #include"TFEL/Math/TinyMatrixSolve.hxx"
35 #include"TFEL/Math/Vector/TinyVectorFromTinyVectorView.hxx"
36 #include"TFEL/Math/Stensor/StensorFromTinyVectorView.hxx"
37 #include"TFEL/Math/ST2toST2/ST2toST2FromTinyMatrixView.hxx"
38 
39 #include"TFEL/Math/ST2toST2/ST2toST2FromTinyMatrixView2.hxx"
40 
41 #include"TFEL/Math/Stensor/StensorFromTinyMatrixColumnView.hxx"
42 #include"TFEL/Math/Stensor/StensorFromTinyMatrixRowView.hxx"
43 #include"TFEL/Math/Stensor/StensorFromTinyMatrixColumnView2.hxx"
44 #include"TFEL/Math/Stensor/StensorFromTinyMatrixRowView2.hxx"
45 #include"TFEL/Math/Vector/TinyVectorOfStensorFromTinyVectorView.hxx"
46 #line 7 "EllipticCreep-1.mfront"
47 #include "TFEL/Material/Lame.hxx"
48 #include"UO2_YoungModulus_Fink1981-mfront.hxx"
49 
50 #include"UO2_PoissonRatio_Fink1981-mfront.hxx"
51 
52 #include"UO2_ThermalExpansion-mfront.hxx"
53 
54 
55 namespace tfel{
56 
57 namespace material{
58 
59 struct EllipticCreepParametersInitializer
60 {
61 static EllipticCreepParametersInitializer&
62 get();
63 
64 double A;
65 double E;
66 double CAf;
67 double referenceTemperatureForThermalExpansion;
68 double epsilon;
69 double theta;
70 double numerical_jacobian_epsilon;
71 unsigned short iterMax;
72 
73 void set(const char* const,const double);
74 
75 void set(const char* const,const unsigned short);
76 
77 private :
78 
79 EllipticCreepParametersInitializer();
80 
81 EllipticCreepParametersInitializer(const EllipticCreepParametersInitializer&);
82 
83 EllipticCreepParametersInitializer&
84 operator=(const EllipticCreepParametersInitializer&);
85 
86 };
87 
88 // Forward Declaration
89 template<ModellingHypothesis::Hypothesis,typename Type,bool use_qt>
90 class EllipticCreep;
91 
92 // Forward Declaration
93 template<ModellingHypothesis::Hypothesis hypothesis,typename Type>
94 std::ostream&
95  operator <<(std::ostream&,const EllipticCreep<hypothesis,Type,false>&);
96 
97 /*!
98 * \class EllipticCreep
99 * \brief This class implements the EllipticCreep behaviour.
100 * \param hypothesis, modelling hypothesis.
101 * \param Type, numerical type.
102 * \author Maxime Salvo / Thomas Helfer
103 * \date   9 Octobre 2013
104 */
105 template<ModellingHypothesis::Hypothesis hypothesis,typename Type>
106 class EllipticCreep<hypothesis,Type,false>
107 : public MechanicalBehaviour<hypothesis,Type,false>,
108 public EllipticCreepBehaviourData<hypothesis,Type,false>,
109 public EllipticCreepIntegrationData<hypothesis,Type,false>
110 {
111 
112 static constexpr unsigned short N = ModellingHypothesisToSpaceDimension<hypothesis>::value;
113 
114 TFEL_STATIC_ASSERT(N==1||N==2||N==3);
115 TFEL_STATIC_ASSERT(tfel::typetraits::IsFundamentalNumericType<Type>::cond);
116 TFEL_STATIC_ASSERT(tfel::typetraits::IsReal<Type>::cond);
117 
118 friend std::ostream& operator<< <>(std::ostream&,const EllipticCreep&);
119 
120 static constexpr unsigned short TVectorSize = N;
121 typedef tfel::math::StensorDimeToSize<N> StensorDimeToSize;
122 static constexpr unsigned short StensorSize = StensorDimeToSize::value;
123 typedef tfel::math::TensorDimeToSize<N> TensorDimeToSize;
124 static constexpr unsigned short TensorSize = TensorDimeToSize::value;
125 
126 typedef unsigned short ushort;
127 typedef tfel::config::Types<N,Type,false> Types;
128 typedef typename Types::real                              real;
129 typedef typename Types::time                              time;
130 typedef typename Types::length                            length;
131 typedef typename Types::frequency                         frequency;
132 typedef typename Types::stress                            stress;
133 typedef typename Types::strain                            strain;
134 typedef typename Types::strainrate                        strainrate;
135 typedef typename Types::stressrate                        stressrate;
136 typedef typename Types::temperature                       temperature;
137 typedef typename Types::thermalexpansion                  thermalexpansion;
138 typedef typename Types::massdensity                       massdensity;
139 typedef typename Types::TVector                           TVector;
140 typedef typename Types::Stensor                           Stensor;
141 typedef typename Types::Stensor4                          Stensor4;
142 typedef typename Types::FrequencyStensor                  FrequencyStensor;
143 typedef typename Types::ForceTVector                      ForceTVector;
144 typedef typename Types::StressStensor                     StressStensor;
145 typedef typename Types::StrainRateStensor                 StressRateStensor;
146 typedef typename Types::DisplacementTVector               DisplacementTVector;
147 typedef typename Types::StrainStensor                     StrainStensor;
148 typedef typename Types::StrainRateStensor                 StrainRateStensor;
149 typedef typename Types::StiffnessTensor                   StiffnessTensor;
150 typedef typename Types::Tensor                            Tensor;
151 typedef typename Types::ThermalExpansionCoefficientTensor ThermalExpansionCoefficientTensor;
152 typedef typename Types::DeformationGradientTensor         DeformationGradientTensor;
153 typedef typename Types::FiniteStrainStiffnessTensor       FiniteStrainStiffnessTensor;
154 typedef StiffnessTensor StiffnessOperator;
155 
156 public :
157 
158 typedef EllipticCreepBehaviourData<hypothesis,Type,false> BehaviourData;
159 typedef EllipticCreepIntegrationData<hypothesis,Type,false> IntegrationData;
160 typedef typename MechanicalBehaviour<hypothesis,Type,false>::SMType            SMType;
161 
162 using MechanicalBehaviour<hypothesis,Type,false>::ELASTIC;
163 using MechanicalBehaviour<hypothesis,Type,false>::SECANTOPERATOR;
164 using MechanicalBehaviour<hypothesis,Type,false>::TANGENTOPERATOR;
165 using MechanicalBehaviour<hypothesis,Type,false>::CONSISTENTTANGENTOPERATOR;
166 using MechanicalBehaviour<hypothesis,Type,false>::NOSTIFFNESSREQUESTED;
167 typedef typename MechanicalBehaviour<hypothesis,Type,false>::IntegrationResult IntegrationResult;
168 using MechanicalBehaviour<hypothesis,Type,false>::SUCCESS;
169 using MechanicalBehaviour<hypothesis,Type,false>::FAILURE;
170 using MechanicalBehaviour<hypothesis,Type,false>::UNRELIABLE_RESULTS;
171 
172 typedef StrainStensor StressFreeExpansionType;
173 
174 private :
175 
176 
177 typename tfel::math::StensorFromTinyVectorView<N,2+StensorSize,0,real>::type deel;
178 #line 13 "EllipticCreep-1.mfront"
179 real& dpv;
180 #line 16 "EllipticCreep-1.mfront"
181 real& dpd;
182 
183 #line 30 "EllipticCreep-1.mfront"
184 real nu;
185 #line 32 "EllipticCreep-1.mfront"
186 stress lambda;
187 #line 34 "EllipticCreep-1.mfront"
188 stress mu;
189 #line 36 "EllipticCreep-1.mfront"
190 real f_t;
191 
192 #line 22 "EllipticCreep-1.mfront"
193 real A;
194 #line 24 "EllipticCreep-1.mfront"
195 real E;
196 #line 26 "EllipticCreep-1.mfront"
197 real CAf;
198 real referenceTemperatureForThermalExpansion;
199 real epsilon;
200 real theta;
201 real numerical_jacobian_epsilon;
202 ushort iterMax;
203 
204 // Jacobian
205 tfel::math::tmatrix<2+StensorSize,2+StensorSize,Type> jacobian;
206 // zeros
207 tfel::math::tvector<2+StensorSize,Type> zeros;
208 
209 // previous zeros
210 tfel::math::tvector<2+StensorSize,Type> zeros_1;
211 
212 // function
213 tfel::math::tvector<2+StensorSize,Type> fzeros;
214 
215 // number of iterations
216 unsigned int iter;
217 
218 void
computeNumericalJacobian(tfel::math::tmatrix<2+StensorSize,2+StensorSize,real> & njacobian)219 computeNumericalJacobian(tfel::math::tmatrix<2+StensorSize,2+StensorSize,real>& njacobian)
220 {
221 using namespace std;
222 using namespace tfel::math;
223 tvector<2+StensorSize,real> tzeros(this->zeros);
224 tvector<2+StensorSize,real> tfzeros(this->fzeros);
225 tmatrix<2+StensorSize,2+StensorSize,real> tjacobian(this->jacobian);
226 for(unsigned short idx = 0; idx!= 2+StensorSize;++idx){
227 this->zeros(idx) -= this->numerical_jacobian_epsilon;
228 this->computeStress();
229 this->computeFdF();
230 this->zeros = tzeros;
231 tvector<2+StensorSize,real> tfzeros2(this->fzeros);
232 this->zeros(idx) += this->numerical_jacobian_epsilon;
233 this->computeStress();
234 this->computeFdF();
235 this->fzeros = (this->fzeros-tfzeros2)/(real(2)*(this->numerical_jacobian_epsilon));
236 for(unsigned short idx2 = 0; idx2!= 2+StensorSize;++idx2){
237 njacobian(idx2,idx) = this->fzeros(idx2);
238 }
239 this->zeros    = tzeros;
240 this->fzeros   = tfzeros;
241 }
242 if(&njacobian!=&(this->jacobian)){
243 this->jacobian = tjacobian;
244 }
245 }
246 
247 void
getPartialJacobianInvert(Stensor4 & partial_jacobian_eel)248 getPartialJacobianInvert(Stensor4& partial_jacobian_eel)
249 {
250 using namespace tfel::math;
251 TinyPermutation<2+StensorSize> permuation;
252 this->computeNumericalJacobian(this->jacobian);
253 TinyMatrixSolve<2+StensorSize,real>::decomp(this->jacobian,permuation);
254 for(unsigned short idx=0;idx!=StensorSize;++idx){
255 tvector<2+StensorSize,real> vect_e(real(0));
256 vect_e(idx) = real(1);
257 TinyMatrixSolve<2+StensorSize,real>::back_substitute(this->jacobian,permuation,vect_e);
258 for(unsigned short idx2=0;idx2!=StensorSize;++idx2){
259 partial_jacobian_eel(idx2,idx)=vect_e(idx2);
260 }
261 }
262 }
263 
264 void
getPartialJacobianInvert(Stensor4 & partial_jacobian_eel,Stensor & partial_jacobian_pv)265 getPartialJacobianInvert(Stensor4& partial_jacobian_eel,
266 Stensor& partial_jacobian_pv)
267 {
268 using namespace tfel::math;
269 TinyPermutation<2+StensorSize> permuation;
270 this->computeNumericalJacobian(this->jacobian);
271 TinyMatrixSolve<2+StensorSize,real>::decomp(this->jacobian,permuation);
272 for(unsigned short idx=0;idx!=StensorSize;++idx){
273 tvector<2+StensorSize,real> vect_e(real(0));
274 vect_e(idx) = real(1);
275 TinyMatrixSolve<2+StensorSize,real>::back_substitute(this->jacobian,permuation,vect_e);
276 for(unsigned short idx2=0;idx2!=StensorSize;++idx2){
277 partial_jacobian_eel(idx2,idx)=vect_e(idx2);
278 }
279 partial_jacobian_pv(idx)=vect_e(StensorSize);
280 }
281 }
282 
283 void
getPartialJacobianInvert(Stensor4 & partial_jacobian_eel,Stensor & partial_jacobian_pv,Stensor & partial_jacobian_pd)284 getPartialJacobianInvert(Stensor4& partial_jacobian_eel,
285 Stensor& partial_jacobian_pv,
286 Stensor& partial_jacobian_pd)
287 {
288 using namespace tfel::math;
289 TinyPermutation<2+StensorSize> permuation;
290 this->computeNumericalJacobian(this->jacobian);
291 TinyMatrixSolve<2+StensorSize,real>::decomp(this->jacobian,permuation);
292 for(unsigned short idx=0;idx!=StensorSize;++idx){
293 tvector<2+StensorSize,real> vect_e(real(0));
294 vect_e(idx) = real(1);
295 TinyMatrixSolve<2+StensorSize,real>::back_substitute(this->jacobian,permuation,vect_e);
296 for(unsigned short idx2=0;idx2!=StensorSize;++idx2){
297 partial_jacobian_eel(idx2,idx)=vect_e(idx2);
298 }
299 partial_jacobian_pv(idx)=vect_e(StensorSize);
300 partial_jacobian_pd(idx)=vect_e(1+StensorSize);
301 }
302 }
303 
304 void
computeStress(void)305 computeStress(void){
306 using namespace std;
307 using namespace tfel::math;
308 using std::vector;
309 using mfront::UO2_YoungModulus_Fink1981;
310 using mfront::UO2_YoungModulus_Fink1981_checkBounds;
311 using mfront::UO2_PoissonRatio_Fink1981;
312 using mfront::UO2_PoissonRatio_Fink1981_checkBounds;
313 using mfront::UO2_ThermalExpansion;
314 using mfront::UO2_ThermalExpansion_checkBounds;
315 #line 51 "EllipticCreep-1.mfront"
316 using namespace tfel :: material :: lame ;
317 #line 53 "EllipticCreep-1.mfront"
318 this->f = this->f_t + this->theta * ( 1 - this->f_t ) / ( 1 + this->theta * this->dpv ) * this->dpv ;
319 #line 55 "EllipticCreep-1.mfront"
320 const stress young = UO2_YoungModulus_Fink1981 ( (this->T+(this->theta)*(this->dT)) , min ( max ( this->f , real ( 0 ) ) , real ( 1 ) ) ) ;
321 #line 56 "EllipticCreep-1.mfront"
322 this->lambda = computeLambda ( young , this->nu ) ;
323 #line 57 "EllipticCreep-1.mfront"
324 this->mu = computeMu ( young , this->nu ) ;
325 #line 59 "EllipticCreep-1.mfront"
326 this->sig = this->lambda * trace ( (this->eel+(this->theta)*(this->deel)) ) * Stensor :: Id ( ) +2 * this->mu * (this->eel+(this->theta)*(this->deel)) ;
327 } // end of EllipticCreep::computeStress
328 
329 void
computeFinalStress(void)330 computeFinalStress(void){
331 using namespace std;
332 using namespace tfel::math;
333 using std::vector;
334 using mfront::UO2_YoungModulus_Fink1981;
335 using mfront::UO2_YoungModulus_Fink1981_checkBounds;
336 using mfront::UO2_PoissonRatio_Fink1981;
337 using mfront::UO2_PoissonRatio_Fink1981_checkBounds;
338 using mfront::UO2_ThermalExpansion;
339 using mfront::UO2_ThermalExpansion_checkBounds;
340 #line 63 "EllipticCreep-1.mfront"
341 using namespace tfel :: material :: lame ;
342 #line 65 "EllipticCreep-1.mfront"
343 this->f = this->f_t + ( 1 - this->f_t ) / ( 1 + this->theta * this->dpv ) * this->dpv ;
344 #line 67 "EllipticCreep-1.mfront"
345 const stress young = UO2_YoungModulus_Fink1981 ( (this->T+this->dT) , min ( max ( this->f , real ( 0 ) ) , real ( 1 ) ) ) ;
346 #line 68 "EllipticCreep-1.mfront"
347 this->lambda = computeLambda ( young , this->nu ) ;
348 #line 69 "EllipticCreep-1.mfront"
349 this->mu = computeMu ( young , this->nu ) ;
350 #line 71 "EllipticCreep-1.mfront"
351 this->sig = this->lambda * trace ( this->eel ) * Stensor :: Id ( ) +2 * this->mu * this->eel ;
352 } // end of EllipticCreep::computeStress
353 
354 
355 /*!
356 * \brief Update internal variables at end of integration
357 */
358 void
updateStateVariables(void)359 updateStateVariables(void){
360 this->eel += this->deel;
361 this->pv += this->dpv;
362 this->pd += this->dpd;
363 }
364 
365 /*!
366 * \brief Update auxiliary state variables at end of integration
367 */
368 void
updateAuxiliaryStateVariables(void)369 updateAuxiliaryStateVariables(void)
370 {}
371 
372 /*!
373 * \brief Default constructor (disabled)
374 */
375 EllipticCreep();
376 
377 /*!
378 * \brief Copy constructor (disabled)
379 */
380 EllipticCreep(const EllipticCreep&);
381 
382 /*!
383 * \brief Assignement operator (disabled)
384 */
385 EllipticCreep& operator = (const EllipticCreep&);
386 
387 public:
388 
389 /*!
390 * \brief Constructor
391 */
EllipticCreep(const EllipticCreepBehaviourData<hypothesis,Type,false> & src1,const EllipticCreepIntegrationData<hypothesis,Type,false> & src2)392 EllipticCreep(const EllipticCreepBehaviourData<hypothesis,Type,false>& src1,
393 const EllipticCreepIntegrationData<hypothesis,Type,false>& src2)
394 : EllipticCreepBehaviourData<hypothesis,Type,false>(src1),
395 EllipticCreepIntegrationData<hypothesis,Type,false>(src2),
396 deel(this->zeros),
397 dpv(this->zeros(StensorSize)),
398 dpd(this->zeros(1+StensorSize)),
399 zeros(real(0)),
400 fzeros(real(0))
401 {
402 using namespace std;
403 using namespace tfel::math;
404 using std::vector;
405 using mfront::UO2_YoungModulus_Fink1981;
406 using mfront::UO2_YoungModulus_Fink1981_checkBounds;
407 using mfront::UO2_PoissonRatio_Fink1981;
408 using mfront::UO2_PoissonRatio_Fink1981_checkBounds;
409 using mfront::UO2_ThermalExpansion;
410 using mfront::UO2_ThermalExpansion_checkBounds;
411 this->A = EllipticCreepParametersInitializer::get().A;
412 this->E = EllipticCreepParametersInitializer::get().E;
413 this->CAf = EllipticCreepParametersInitializer::get().CAf;
414 this->referenceTemperatureForThermalExpansion = EllipticCreepParametersInitializer::get().referenceTemperatureForThermalExpansion;
415 this->epsilon = EllipticCreepParametersInitializer::get().epsilon;
416 this->theta = EllipticCreepParametersInitializer::get().theta;
417 this->numerical_jacobian_epsilon = EllipticCreepParametersInitializer::get().numerical_jacobian_epsilon;
418 this->iterMax = EllipticCreepParametersInitializer::get().iterMax;
419 }
420 
421 /*
422  * \brief constructor for the umat interface
423  *
424  * \param const Type *const UMATdt_, time increment
425  * \param const Type *const UMATT_, temperature
426  * \param const Type *const UMATdT_, temperature increment
427  * \param const Type *const UMATmat, material properties
428  * \param const Type *const UMATint_vars, state variables
429  * \param const Type *const UMAText_vars, external state variables
430  * \param const Type *const UMATdext_vars, external state variables increments
431  */
EllipticCreep(const Type * const UMATdt_,const Type * const UMATT_,const Type * const UMATdT_,const Type * const UMATmat,const Type * const UMATint_vars,const Type * const UMAText_vars,const Type * const UMATdext_vars)432 EllipticCreep(const Type* const UMATdt_,
433 const Type* const UMATT_,const Type* const UMATdT_,
434 const Type* const UMATmat,const Type* const UMATint_vars,
435 const Type* const UMAText_vars,const Type* const UMATdext_vars)
436 : EllipticCreepBehaviourData<hypothesis,Type,false>(UMATT_,UMATmat,
437 UMATint_vars,UMAText_vars),
438 EllipticCreepIntegrationData<hypothesis,Type,false>(UMATdt_,UMATdT_,UMATdext_vars),
439 deel(this->zeros),
440 dpv(this->zeros(StensorSize)),
441 dpd(this->zeros(1+StensorSize)),
442 zeros(real(0)),
443 fzeros(real(0))
444 {
445 using namespace std;
446 using namespace tfel::math;
447 using std::vector;
448 using mfront::UO2_YoungModulus_Fink1981;
449 using mfront::UO2_YoungModulus_Fink1981_checkBounds;
450 using mfront::UO2_PoissonRatio_Fink1981;
451 using mfront::UO2_PoissonRatio_Fink1981_checkBounds;
452 using mfront::UO2_ThermalExpansion;
453 using mfront::UO2_ThermalExpansion_checkBounds;
454 this->A = EllipticCreepParametersInitializer::get().A;
455 this->E = EllipticCreepParametersInitializer::get().E;
456 this->CAf = EllipticCreepParametersInitializer::get().CAf;
457 this->referenceTemperatureForThermalExpansion = EllipticCreepParametersInitializer::get().referenceTemperatureForThermalExpansion;
458 this->epsilon = EllipticCreepParametersInitializer::get().epsilon;
459 this->theta = EllipticCreepParametersInitializer::get().theta;
460 this->numerical_jacobian_epsilon = EllipticCreepParametersInitializer::get().numerical_jacobian_epsilon;
461 this->iterMax = EllipticCreepParametersInitializer::get().iterMax;
462 }
463 
464 void
computeStressFreeExpansion(std::pair<StressFreeExpansionType,StressFreeExpansionType> & dl_l01)465 computeStressFreeExpansion(std::pair<StressFreeExpansionType,StressFreeExpansionType>& dl_l01)
466 {
467 using namespace std;
468 using namespace tfel::math;
469 using std::vector;
470 using mfront::UO2_YoungModulus_Fink1981;
471 using mfront::UO2_YoungModulus_Fink1981_checkBounds;
472 using mfront::UO2_PoissonRatio_Fink1981;
473 using mfront::UO2_PoissonRatio_Fink1981_checkBounds;
474 using mfront::UO2_ThermalExpansion;
475 using mfront::UO2_ThermalExpansion_checkBounds;
476 auto& dl_l0 = dl_l01.first;
477 auto& dl_l1 = dl_l01.second;
478 dl_l0 = StressFreeExpansionType(typename StressFreeExpansionType::value_type(0));
479 dl_l1 = StressFreeExpansionType(typename StressFreeExpansionType::value_type(0));
480 const thermalexpansion alpha_Ti        = UO2_ThermalExpansion();
481 const thermalexpansion alpha_T_t       = UO2_ThermalExpansion();
482 const thermalexpansion alpha_T_t_dt    = UO2_ThermalExpansion();
483 dl_l0[0] = 1/(1+alpha_Ti * (this->referenceTemperatureForThermalExpansion-300))*(alpha_T_t * (this->T-300)-alpha_Ti * (this->referenceTemperatureForThermalExpansion-300));
484 dl_l0[1] = dl_l0[0];
485 dl_l0[2] = dl_l0[0];
486 dl_l1[0] = 1/(1+alpha_Ti * (this->referenceTemperatureForThermalExpansion-300))*(alpha_T_t_dt * (this->T+this->dT-300)-alpha_Ti * (this->referenceTemperatureForThermalExpansion-300));
487 dl_l1[1] = dl_l1[0];
488 dl_l1[2] = dl_l1[0];
489 }
490 
491 /*!
492  * \ brief initialize the behaviour with user code
493  */
initialize(void)494 void initialize(void){
495 using namespace std;
496 using namespace tfel::math;
497 using std::vector;
498 using mfront::UO2_YoungModulus_Fink1981;
499 using mfront::UO2_YoungModulus_Fink1981_checkBounds;
500 using mfront::UO2_PoissonRatio_Fink1981;
501 using mfront::UO2_PoissonRatio_Fink1981_checkBounds;
502 using mfront::UO2_ThermalExpansion;
503 using mfront::UO2_ThermalExpansion_checkBounds;
504 #line 45 "EllipticCreep-1.mfront"
505 this->nu = UO2_PoissonRatio_Fink1981 ( ) ;
506 #line 47 "EllipticCreep-1.mfront"
507 this->f_t = this->f ;
508 }
509 
510 /*!
511 * \brief set the policy for "out of bounds" conditions
512 */
513 void
setOutOfBoundsPolicy(const OutOfBoundsPolicy policy_value)514 setOutOfBoundsPolicy(const OutOfBoundsPolicy policy_value){
515 this->policy = policy_value;
516 } // end of setOutOfBoundsPolicy
517 
518 /*!
519 * \brief set the policy for "out of bounds" conditions
520 */
521 ModellingHypothesis::Hypothesis
getModellingHypothesis(void) const522 getModellingHypothesis(void) const{
523 return hypothesis;
524 } // end of getModellingHypothesis
525 
526 /*!
527 * \brief check bounds
528 */
529 void
checkBounds(void) const530 checkBounds(void) const{
531 } // end of checkBounds
532 
computePredictionOperator(const SMType)533 IntegrationResult computePredictionOperator(const SMType){
534 using namespace std;
535 string msg("EllipticCreep::computePredictionOperator : ");
536 msg +="unimplemented feature";
537 throw(runtime_error(msg));
538 return FAILURE;
539 }
540 
541 /*!
542 * \brief Integrate behaviour law over the time step
543 */
544 IntegrationResult
integrate(const SMType smt)545 integrate(const SMType smt){
546 using namespace std;
547 using namespace tfel::math;
548 real error;
549 bool converge=false;
550 this->iter=0;
551 while((converge==false)&&
552 (this->iter<EllipticCreep::iterMax)){
553 ++(this->iter);
554 this->computeStress();
555 const bool computeFdF_ok = this->computeFdF();
556 if(!computeFdF_ok){
557 if(this->iter==1){
558 return MechanicalBehaviour<hypothesis,Type,false>::FAILURE;
559 } else {
560 const real integrate_one_half = real(1)/real(2);
561 this->zeros -= (this->zeros-this->zeros_1)*integrate_one_half;
562 }
563 } else {
564 this->zeros_1  = this->zeros;
565 error=norm(this->fzeros);
566 converge = ((error)/(real(2+StensorSize))<(this->epsilon));
567 if(!converge){
568 try{
569 this->computeNumericalJacobian(this->jacobian);
570 TinyMatrixSolve<2+StensorSize,real>::exe(this->jacobian,this->fzeros);
571 }
572 catch(LUException&){
573 return MechanicalBehaviour<hypothesis,Type,false>::FAILURE;
574 }
575 this->zeros -= this->fzeros;
576 }
577 }
578 }
579 if(this->iter==this->iterMax){
580 return MechanicalBehaviour<hypothesis,Type,false>::FAILURE;
581 }
582 if(smt!=NOSTIFFNESSREQUESTED){
583 string msg("EllipticCreep::integrate : ");
584 msg +="unimplemented feature";
585 throw(runtime_error(msg));
586 }
587 this->updateStateVariables();
588 this->computeFinalStress();
589 this->updateAuxiliaryStateVariables();
590 return MechanicalBehaviour<hypothesis,Type,false>::SUCCESS;
591 } // end of EllipticCreep::integrate
592 
593 /*!
594 * \brief compute fzeros and jacobian
595 */
596 bool
computeFdF(void)597 computeFdF(void){
598 using namespace std;
599 using namespace tfel::math;
600 using std::vector;
601 using mfront::UO2_YoungModulus_Fink1981;
602 using mfront::UO2_YoungModulus_Fink1981_checkBounds;
603 using mfront::UO2_PoissonRatio_Fink1981;
604 using mfront::UO2_PoissonRatio_Fink1981_checkBounds;
605 using mfront::UO2_ThermalExpansion;
606 using mfront::UO2_ThermalExpansion_checkBounds;
607 typename tfel::math::StensorFromTinyVectorView<N,2+StensorSize,0,real>::type feel(this->fzeros);
608 real& fpv(this->fzeros(StensorSize));
609 real& fpd(this->fzeros(1+StensorSize));
610 // setting f values to zeros
611 this->fzeros = this->zeros;
612 #line 76 "EllipticCreep-1.mfront"
613 const stress pr = trace ( (this->sig) ) / real ( 3 ) ;
614 #line 78 "EllipticCreep-1.mfront"
615 const stress seq = sigmaeq ( (this->sig) ) ;
616 #line 80 "EllipticCreep-1.mfront"
617 const real Af = 2.25 * (this->CAf) * pow ( (this->E) * ( pow ( (this->f) , - (this->E) ) -1 ) , -2 * (this->E) / ( (this->E) + 1 ) ) ;
618 #line 81 "EllipticCreep-1.mfront"
619 const real Bf = ( 1 + 2 * (this->f) / real ( 3 ) ) * pow ( 1 - (this->f) , -2 * (this->E) / ( (this->E) + 1 ) ) ;
620 #line 82 "EllipticCreep-1.mfront"
621 const stress s = sqrt ( Af * pr * pr + Bf * seq * seq ) ;
622 #line 83 "EllipticCreep-1.mfront"
623 if ( s > 1.e-8 * (this->mu) ) {
624 #line 85 "EllipticCreep-1.mfront"
625 real inv_seq ( 0 ) ;
626 #line 86 "EllipticCreep-1.mfront"
627 Stensor n ( real ( 0. ) ) ;
628 #line 87 "EllipticCreep-1.mfront"
629 if ( seq > 1.e-8 * (this->mu) ) {
630 #line 88 "EllipticCreep-1.mfront"
631 inv_seq = 1 / seq ;
632 #line 89 "EllipticCreep-1.mfront"
633 n = 1.5 * deviator ( (this->sig) ) * inv_seq ;
634 #line 90 "EllipticCreep-1.mfront"
635 }
636 #line 91 "EllipticCreep-1.mfront"
637 const real ds_dpr = Af * pr / s ;
638 #line 92 "EllipticCreep-1.mfront"
639 const real ds_dseq = Bf * seq / s ;
640 #line 93 "EllipticCreep-1.mfront"
641 const real dphi_ds = (this->A) * pow ( s , (this->E) ) ;
642 #line 94 "EllipticCreep-1.mfront"
643 const real dphi_dp = dphi_ds * ds_dpr ;
644 #line 95 "EllipticCreep-1.mfront"
645 const real dphi_dseq = dphi_ds * ds_dseq ;
646 #line 97 "EllipticCreep-1.mfront"
647 const real K = (this->lambda) + ( 2 * (this->mu) ) / 3 ;
648 #line 98 "EllipticCreep-1.mfront"
649 fpv -= dphi_dp * (this->dt) ;
650 #line 100 "EllipticCreep-1.mfront"
651 fpd -= dphi_dseq * (this->dt) ;
652 #line 102 "EllipticCreep-1.mfront"
653 feel += ( (this->dpv) / 3 ) * Stensor :: Id ( ) + (this->dpd) * n ;
654 #line 103 "EllipticCreep-1.mfront"
655 }
656 #line 104 "EllipticCreep-1.mfront"
657 feel -= (this->deto) ;
658 return true;
659 }
660 
661 const StiffnessOperator&
getTangentOperator(void) const662 getTangentOperator(void) const{
663 return this->Dt;
664 }
665 
666 real
getTimeStepScalingFactor(void) const667 getTimeStepScalingFactor(void) const{
668 return real(1);
669 }
670 
updateExternalStateVariables(void)671 void updateExternalStateVariables(void){
672 this->eto  += this->deto;
673 
674 this->T   += this->dT;
675 }
676 
677 /*!
678 * \brief Destructor
679 */
~EllipticCreep()680 ~EllipticCreep()
681 {}
682 
683 private:
684 
685 //! Tangent operator;
686 StiffnessOperator Dt;
687 //! policy for treating out of bounds conditions
688 OutOfBoundsPolicy policy;
689 }; // end of EllipticCreep class
690 
691 template<ModellingHypothesis::Hypothesis hypothesis,typename Type>
692 std::ostream&
operator <<(std::ostream & os,const EllipticCreep<hypothesis,Type,false> & b)693 operator <<(std::ostream& os,const EllipticCreep<hypothesis,Type,false>& b)
694 {
695 using namespace std;
696 using namespace tfel::utilities;
697 os << Name<EllipticCreep<hypothesis,Type,false> >::getName() << endl;
698 os << "eto : " << b.eto << endl;
699 os << "deto : " << b.deto << endl;
700 os << "sig : " << b.sig << endl;
701 os << "dt : " << b.dt << endl;
702 os << "T : " << b.T << endl;
703 os << "dT : " << b.dT << endl;
704 os << "eel : " << b.eel << endl;
705 os << "deel : " << b.deel << endl;
706 os << "pv : " << b.pv << endl;
707 os << "dpv : " << b.dpv << endl;
708 os << "pd : " << b.pd << endl;
709 os << "dpd : " << b.dpd << endl;
710 os << "f : " << b.f << endl;
711 os << "nu : " << b.nu << endl;
712 os << "lambda : " << b.lambda << endl;
713 os << "mu : " << b.mu << endl;
714 os << "f_t : " << b.f_t << endl;
715 os << "A : " << b.A << endl;
716 os << "E : " << b.E << endl;
717 os << "CAf : " << b.CAf << endl;
718 os << "referenceTemperatureForThermalExpansion : " << b.referenceTemperatureForThermalExpansion << endl;
719 os << "epsilon : " << b.epsilon << endl;
720 os << "theta : " << b.theta << endl;
721 os << "numerical_jacobian_epsilon : " << b.numerical_jacobian_epsilon << endl;
722 os << "iterMax : " << b.iterMax << endl;
723 return os;
724 }
725 
726 /*!
727 * Partial specialisation for EllipticCreep.
728 */
729 template<ModellingHypothesis::Hypothesis hypothesis,typename Type>
730 class MechanicalBehaviourTraits<EllipticCreep<hypothesis,Type,false> >
731 {
732 static constexpr unsigned short N = ModellingHypothesisToSpaceDimension<hypothesis>::value;
733 static constexpr unsigned short TVectorSize = N;
734 typedef tfel::math::StensorDimeToSize<N> StensorDimeToSize;
735 static constexpr unsigned short StensorSize = StensorDimeToSize::value;
736 typedef tfel::math::TensorDimeToSize<N> TensorDimeToSize;
737 static constexpr unsigned short TensorSize = TensorDimeToSize::value;
738 public:
739 static constexpr bool is_defined = true;
740 static constexpr bool use_quantities = false;
741 static constexpr bool hasStressFreeExpansion = true;
742 static constexpr bool handlesThermalExpansion = true;
743 static constexpr unsigned short dimension = N;
744 typedef Type NumType;
745 static constexpr unsigned short material_properties_nb = 0;
746 static constexpr unsigned short internal_variables_nb  = 3+StensorSize;
747 static constexpr unsigned short external_variables_nb  = 0;
748 static constexpr bool hasConsistantTangentOperator = false;
749 static constexpr bool isConsistantTangentOperatorSymmetric = false;
750 static constexpr bool hasPredictionOperator = false;
751 static constexpr bool hasTimeStepScalingFactor = false;
752 };
753 
754 /*!
755 * Partial specialisation for EllipticCreep.
756 */
757 template<typename Type>
758 class MechanicalBehaviourTraits<EllipticCreep<ModellingHypothesis::AXISYMMETRICALGENERALISEDPLANESTRESS,Type,false> >
759 {
760 public:
761 static constexpr bool is_defined = false;
762 static constexpr bool use_quantities = false;
763 static constexpr bool hasStressFreeExpansion = true;
764 static constexpr bool handlesThermalExpansion = true;
765 static constexpr unsigned short dimension = 0u;
766 typedef Type NumType;
767 static constexpr unsigned short material_properties_nb = 0;
768 static constexpr unsigned short internal_variables_nb  = 0;
769 static constexpr unsigned short external_variables_nb  = 0;
770 static constexpr bool hasConsistantTangentOperator = false;
771 static constexpr bool isConsistantTangentOperatorSymmetric = false;
772 static constexpr bool hasPredictionOperator = false;
773 static constexpr bool hasTimeStepScalingFactor = false;
774 };
775 
776 /*!
777 * Partial specialisation for EllipticCreep.
778 */
779 template<typename Type>
780 class MechanicalBehaviourTraits<EllipticCreep<ModellingHypothesis::PLANESTRESS,Type,false> >
781 {
782 public:
783 static constexpr bool is_defined = false;
784 static constexpr bool use_quantities = false;
785 static constexpr bool hasStressFreeExpansion = true;
786 static constexpr bool handlesThermalExpansion = true;
787 static constexpr unsigned short dimension = 0u;
788 typedef Type NumType;
789 static constexpr unsigned short material_properties_nb = 0;
790 static constexpr unsigned short internal_variables_nb  = 0;
791 static constexpr unsigned short external_variables_nb  = 0;
792 static constexpr bool hasConsistantTangentOperator = false;
793 static constexpr bool isConsistantTangentOperatorSymmetric = false;
794 static constexpr bool hasPredictionOperator = false;
795 static constexpr bool hasTimeStepScalingFactor = false;
796 };
797 
798 } // end of namespace material
799 
800 } // end of namespace tfel
801 
802 namespace tfel{
803 namespace utilities{
804 //! Partial specialisation of the Name class
805 template<tfel::material::ModellingHypothesis::Hypothesis hypothesis,typename Type>
806 struct Name<tfel::material::EllipticCreep<hypothesis,Type,false> >
807 {
808 /*!
809 * \brief  Return the name of the class.
810 * \return the name of the class.
811 * \see    Name.
812 */
813 static std::string
getNametfel::utilities::Name814 getName(void){
815 return std::string("EllipticCreep");
816 }
817 }; // end of struct Name
818 } // end of namespace utilities
819 } // end of namespace tfel
820 
821 #endif /* LIB_TFELMATERIAL_ELLIPTICCREEP_HXX_ */
822