OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
binghamfluid2.C
Go to the documentation of this file.
1 /*
2  *
3  * ##### ##### ###### ###### ### ###
4  * ## ## ## ## ## ## ## ### ##
5  * ## ## ## ## #### #### ## # ##
6  * ## ## ## ## ## ## ## ##
7  * ## ## ## ## ## ## ## ##
8  * ##### ##### ## ###### ## ##
9  *
10  *
11  * OOFEM : Object Oriented Finite Element Code
12  *
13  * Copyright (C) 1993 - 2013 Borek Patzak
14  *
15  *
16  *
17  * Czech Technical University, Faculty of Civil Engineering,
18  * Department of Structural Mechanics, 166 29 Prague, Czech Republic
19  *
20  * This library is free software; you can redistribute it and/or
21  * modify it under the terms of the GNU Lesser General Public
22  * License as published by the Free Software Foundation; either
23  * version 2.1 of the License, or (at your option) any later version.
24  *
25  * This program is distributed in the hope that it will be useful,
26  * but WITHOUT ANY WARRANTY; without even the implied warranty of
27  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28  * Lesser General Public License for more details.
29  *
30  * You should have received a copy of the GNU Lesser General Public
31  * License along with this library; if not, write to the Free Software
32  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
33  */
34 
35 #include "binghamfluid2.h"
36 #include "fluiddynamicmaterial.h"
37 #include "domain.h"
38 #include "floatmatrix.h"
39 #include "gausspoint.h"
40 #include "engngm.h"
41 #include "mathfem.h"
42 #include "datastream.h"
43 #include "contextioerr.h"
44 #include "dynamicinputrecord.h"
45 #include "classfactory.h"
46 
47 #include <cstdlib>
48 
49 namespace oofem {
50 #define BINGHAM_ALT 1
51 #define BINGHAM_MIN_SHEAR_RATE 1.e-10
52 
53 REGISTER_Material(BinghamFluidMaterial2);
55 // static bool __dummy_BinghamFluidMaterial2_alt __attribute__((unused)) = GiveClassFactory().registerMaterial("binghamfluid2", matCreator< BinghamFluidMaterial2 > );
56 REGISTER_Material_Alt(BinghamFluidMaterial2, binghamfluid2);
57 
59  mu_0(0.),
60  tau_0(0.),
61  tau_c(0.),
62  mu_inf(1.e6),
63  stressGrowthRate(BINGHAM_DEFAULT_STRESS_GROWTH_RATE)
64 { }
65 
66 
69 {
70  IRResultType result; // Required by IR_GIVE_FIELD macro
71 
72  // we use rather object's member data than to store data into slow
73  // key-val dictionary with lot of memory allocations
76  mu_inf = 1.e6;
80  tau_c = tau_0 * mu_inf / ( mu_inf - mu_0 );
81  //tau_c = tau_0;
83 }
84 
85 
86 void
88 {
94 }
95 
96 double
98 {
99  BinghamFluidMaterial2Status *status = static_cast< BinghamFluidMaterial2Status * >( this->giveStatus(gp) );
100  //double temp_tau=status->giveTempDevStressMagnitude();
101  //double temp_gamma=status->giveTempDevStrainMagnitude();
102  //determine actual viscosity
103  //return this->computeActualViscosity(temp_tau, temp_gamma);
104 #ifdef BINGHAM_ALT
105  double gamma = status->giveTempDevStrainMagnitude(); //status->giveTempDevStrainMagnitude();
106  return computeActualViscosity(tau_0, gamma);
107 
108  #if 0
109  const FloatArray &epsd = status->giveTempDeviatoricStrainVector(); //status->giveTempDeviatoricStrainVector();
110  double gamma2 = gamma * gamma;
111  double dmudg, dgde1, dgde2, dgde3, mu;
112  if ( gamma < BINGHAM_MIN_SHEAR_RATE ) {
113  dmudg = dgde1 = dgde2 = dgde3 = 0.0;
114  mu = computeActualViscosity(tau_0, gamma);
115  } else {
116  dmudg = ( -1.0 ) * tau_0 * ( 1.0 - exp(-this->stressGrowthRate * gamma) ) / gamma2 +
117  tau_0 *this->stressGrowthRate *exp(-this->stressGrowthRate *gamma) / gamma;
118  mu = mu_0 + tau_0 * ( 1. - exp(-this->stressGrowthRate * gamma) ) / gamma;
119 
120  dgde1 = 2.0 * fabs( epsd.at(1) ) / gamma;
121  dgde2 = 2.0 * fabs( epsd.at(2) ) / gamma;
122  dgde3 = 1.0 * fabs( epsd.at(3) ) / gamma;
123  }
124 
125  return min( min( ( epsd.at(1) * dmudg * dgde1 + mu ), ( epsd.at(2) * dmudg * dgde2 + mu ) ),
126  ( epsd.at(3) * dmudg * dgde3 + mu ) );
127 
128  #endif
129 
130 #else
131  if ( temp_tau < tau_c ) {
132  return mu_inf;
133  } else {
134  return mu_0;
135  }
136 
137 #endif
138 }
139 
140 
141 double
143 {
144  if ( aProperty == Viscosity ) {
145  return mu_0;
146  } else if ( aProperty == YieldStress ) {
147  return tau_0;
148  } else {
149  return FluidDynamicMaterial :: give(aProperty, gp);
150  }
151 }
152 
153 
156 {
157  return new BinghamFluidMaterial2Status(1, this->giveDomain(), gp);
158 }
159 
160 
161 void
163 {
164  int size = eps.giveSize();
165  answer.resize(size);
166  BinghamFluidMaterial2Status *status = static_cast< BinghamFluidMaterial2Status * >( this->giveStatus(gp) );
167  MaterialMode mmode = gp->giveMaterialMode();
168  FloatArray epsd;
169  double gamma, tau, _nu;
170 
171  // determine actual viscosity
172  this->computeDeviatoricStrain(epsd, eps, mmode);
173  // determine shear strain magnitude
174  gamma = this->computeDevStrainMagnitude(mmode, epsd);
175 
176 #ifdef BINGHAM_ALT
177  _nu = computeActualViscosity(tau_0, gamma);
178  this->computeDeviatoricStress(answer, epsd, _nu, mmode);
179  tau = this->computeDevStressMagnitude(mmode, answer);
180 
181  //printf ("_nu %e gamma %e\n", _nu, gamma);
182 #else
183  // compute trial state
184  this->computeDeviatoricStress(answer, epsd, this->mu_inf, mmode);
185  // check if state allowed
186  tau = this->computeDevStressMagnitude(mmode, answer);
187  if ( tau > this->tau_c ) {
188  _nu = this->computeActualViscosity(tau, gamma);
189  this->computeDeviatoricStress(answer, epsd, _nu, mmode);
190  tau = this->computeDevStressMagnitude(mmode, answer);
191  }
192 
193 #endif
194  // update status
195  status->letTempDeviatoricStrainVectorBe(epsd);
196  status->letDeviatoricStressVectorBe(answer);
197  status->letTempDevStrainMagnitudeBe(gamma);
198  status->letTempDevStressMagnitudeBe(tau);
199 }
200 
201 void
203  TimeStep *tStep)
204 {
205  BinghamFluidMaterial2Status *status = static_cast< BinghamFluidMaterial2Status * >( this->giveStatus(gp) );
206  MaterialMode mmode = gp->giveMaterialMode();
207  const FloatArray &epsd = status->giveTempDeviatoricStrainVector(); //status->giveTempDeviatoricStrainVector();
209  //double tau = status->giveTempDevStressMagnitude();
210  double gamma = status->giveTempDevStrainMagnitude(); //status->giveTempDevStrainMagnitude();
211  // determine actual viscosity
212  double gamma2 = gamma * gamma;
213 
214  if ( mmode == _2dFlow ) {
215  answer.resize(3, 3);
216  answer.zero();
217 
218  double dgde1, dgde2, dgde3;
219  double dmudg, mu;
220 
221 #if 0
222  double _nu = computeActualViscosity(tau_0, gamma);
223 
224  answer.at(1, 1) = answer.at(2, 2) = 2.0 * _nu;
225  answer.at(3, 3) = _nu;
226 
227  //answer.at(1,1) = answer.at(2,2) = answer.at(3,3) = 2.0*_nu;
228  //answer.at(4,4) = _nu;
229  return;
230 #endif
231 
232  if ( ( mode == ElasticStiffness ) || ( mode == SecantStiffness ) ) {
233  double _nu = computeActualViscosity(tau_0, gamma);
234  answer.at(1, 1) = answer.at(2, 2) = 2.0 * _nu;
235  answer.at(3, 3) = _nu;
236  return;
237  } else { // tangent stiffness
238  if ( gamma < BINGHAM_MIN_SHEAR_RATE ) {
239  dmudg = dgde1 = dgde2 = dgde3 = 0.0;
240  mu = computeActualViscosity(tau_0, gamma);
241  } else {
242  dmudg = ( -1.0 ) * tau_0 * ( 1.0 - exp(-this->stressGrowthRate * gamma) ) / gamma2 +
243  tau_0 *this->stressGrowthRate *exp(-this->stressGrowthRate *gamma) / gamma;
244  mu = mu_0 + tau_0 * ( 1. - exp(-this->stressGrowthRate * gamma) ) / gamma;
245 
246 #if 1
247  //dgde1 = 2.0 * fabs( epsd.at(1) ) / gamma;
248  //dgde2 = 2.0 * fabs( epsd.at(2) ) / gamma;
249  //dgde3 = 1.0 * fabs( epsd.at(3) ) / gamma;
250  dgde1 = 0.5 * fabs( epsd.at(1) ) / gamma;
251  dgde2 = 0.5 * fabs( epsd.at(2) ) / gamma;
252  dgde3 = 1.0 * fabs( epsd.at(3) ) / gamma;
253 #else
254  //dgde1 = 2.0 * epsd.at(1) / gamma;
255  //dgde2 = 2.0 * epsd.at(2) / gamma;
256  //dgde3 = 1.0 * epsd.at(3) / gamma;
257 #endif
258  }
259 
260  answer.at(1, 1) = 2.0 * epsd.at(1) * dmudg * dgde1 + 2.0 * mu;
261  answer.at(1, 2) = 2.0 * epsd.at(1) * dmudg * dgde2;
262  answer.at(1, 3) = 2.0 * epsd.at(1) * dmudg * dgde3;
263 
264  answer.at(2, 1) = 2.0 * epsd.at(2) * dmudg * dgde1;
265  answer.at(2, 2) = 2.0 * epsd.at(2) * dmudg * dgde2 + 2.0 * mu;
266  answer.at(2, 3) = 2.0 * epsd.at(2) * dmudg * dgde3;
267 
268  answer.at(3, 1) = epsd.at(3) * dmudg * dgde1;
269  answer.at(3, 2) = epsd.at(3) * dmudg * dgde2;
270  answer.at(3, 3) = epsd.at(3) * dmudg * dgde3 + mu;
271 
272  return;
273  }
274  } else if ( mmode == _2dAxiFlow ) {
275  answer.resize(4, 4);
276  answer.zero();
277 
278 
279  double dgde1, dgde2, dgde3, dgde4;
280  double dmudg, mu;
281 
282  if ( 0 ) {
283  double _nu = computeActualViscosity(tau_0, gamma);
284 
285  answer.at(1, 1) = answer.at(2, 2) = answer.at(3, 3) = 2.0 * _nu;
286  answer.at(4, 4) = _nu;
287 
288  //answer.at(1,1) = answer.at(2,2) = answer.at(3,3) = 2.0*_nu;
289  //answer.at(4,4) = _nu;
290  return;
291  }
292 
293  if ( ( mode == ElasticStiffness ) || ( mode == SecantStiffness ) ) {
294  double _nu = computeActualViscosity(tau_0, gamma);
295  answer.at(1, 1) = answer.at(2, 2) = answer.at(3, 3) = 2.0 * _nu;
296  answer.at(4, 4) = _nu;
297  return;
298  } else { // tangent stiffness
299  if ( gamma < BINGHAM_MIN_SHEAR_RATE ) {
300  dmudg = dgde1 = dgde2 = dgde3 = dgde4 = 0.0;
301  mu = computeActualViscosity(tau_0, gamma);
302  } else {
303  dmudg = ( -1.0 ) * tau_0 * ( 1.0 - exp(-this->stressGrowthRate * gamma) ) / gamma2 +
304  tau_0 *this->stressGrowthRate *exp(-this->stressGrowthRate *gamma) / gamma;
305  mu = mu_0 + tau_0 * ( 1. - exp(-this->stressGrowthRate * gamma) ) / gamma;
306 
307 #if 1
308  dgde1 = 2.0 * fabs( epsd.at(1) ) / gamma;
309  dgde2 = 2.0 * fabs( epsd.at(2) ) / gamma;
310  dgde3 = 2.0 * fabs( epsd.at(3) ) / gamma;
311  dgde4 = 1.0 * fabs( epsd.at(4) ) / gamma;
312 #else
313  dgde1 = 2.0 * epsd.at(1) / gamma;
314  dgde2 = 2.0 * epsd.at(2) / gamma;
315  dgde3 = 2.0 * epsd.at(3) / gamma;
316  dgde4 = 1.0 * epsd.at(4) / gamma;
317 #endif
318  }
319 
320  answer.at(1, 1) = 2.0 * epsd.at(1) * dmudg * dgde1 + 2.0 * mu;
321  answer.at(1, 2) = 2.0 * epsd.at(1) * dmudg * dgde2;
322  answer.at(1, 3) = 2.0 * epsd.at(1) * dmudg * dgde3;
323  answer.at(1, 4) = 2.0 * epsd.at(1) * dmudg * dgde4;
324 
325  answer.at(2, 1) = 2.0 * epsd.at(2) * dmudg * dgde1;
326  answer.at(2, 2) = 2.0 * epsd.at(2) * dmudg * dgde2 + 2.0 * mu;
327  answer.at(2, 3) = 2.0 * epsd.at(2) * dmudg * dgde3;
328  answer.at(2, 4) = 2.0 * epsd.at(2) * dmudg * dgde4;
329 
330  answer.at(3, 1) = 2.0 * epsd.at(3) * dmudg * dgde1;
331  answer.at(3, 2) = 2.0 * epsd.at(3) * dmudg * dgde2;
332  answer.at(3, 3) = 2.0 * epsd.at(3) * dmudg * dgde3 + 2.0 * mu;
333  answer.at(3, 4) = 2.0 * epsd.at(3) * dmudg * dgde4;
334 
335 
336  answer.at(4, 1) = epsd.at(4) * dmudg * dgde1;
337  answer.at(4, 2) = epsd.at(4) * dmudg * dgde2;
338  answer.at(4, 3) = epsd.at(4) * dmudg * dgde3;
339  answer.at(4, 4) = epsd.at(4) * dmudg * dgde4 + mu;
340 
341  return;
342  }
343  } else if ( gp->giveMaterialMode() == _3dFlow ) {
344  answer.resize(6, 6);
345  answer.zero();
346 
347  FloatArray dgde(6);
348  double dmudg, mu;
349 
350  if ( gamma < BINGHAM_MIN_SHEAR_RATE ) {
351  dmudg = 0.0;
352  dgde.zero();
353  mu = computeActualViscosity(tau_0, gamma);
354  } else {
355  dmudg = ( -1.0 ) * tau_0 * ( 1.0 - exp(-this->stressGrowthRate * gamma) ) / gamma2 +
356  tau_0 *this->stressGrowthRate *exp(-this->stressGrowthRate *gamma) / gamma;
357  mu = mu_0 + tau_0 * ( 1. - exp(-this->stressGrowthRate * gamma) ) / gamma;
358 
359  dgde.at(1) = 2.0 * epsd.at(1) / gamma;
360  dgde.at(2) = 2.0 * epsd.at(2) / gamma;
361  dgde.at(3) = 2.0 * epsd.at(3) / gamma;
362  dgde.at(4) = 1.0 * epsd.at(4) / gamma;
363  dgde.at(5) = 1.0 * epsd.at(5) / gamma;
364  dgde.at(6) = 1.0 * epsd.at(6) / gamma;
365  }
366 
367  for ( int i = 1; i <= 6; i++ ) {
368  answer.at(1, i) = fabs( 2.0 * epsd.at(1) * dmudg * dgde.at(i) );
369  answer.at(2, i) = fabs( 2.0 * epsd.at(2) * dmudg * dgde.at(i) );
370  answer.at(3, i) = fabs( 2.0 * epsd.at(3) * dmudg * dgde.at(i) );
371  answer.at(4, i) = fabs( epsd.at(4) * dmudg * dgde.at(i) );
372  answer.at(5, i) = fabs( epsd.at(5) * dmudg * dgde.at(i) );
373  answer.at(6, i) = fabs( epsd.at(6) * dmudg * dgde.at(i) );
374  }
375 
376  answer.at(1, 1) += 2.0 * mu;
377  answer.at(2, 2) += 2.0 * mu;
378  answer.at(3, 3) += 2.0 * mu;
379  answer.at(4, 4) += mu;
380  answer.at(5, 5) += mu;
381  answer.at(6, 6) += mu;
382 
383  return;
384  } else {
385  OOFEM_ERROR("unsupportted material mode");
386  }
387 }
388 
389 
390 int
392 {
394  double scale;
396  propertyDictionary.at('d') /= scale;
397 
399  this->mu_0 /= scale;
400  this->tau_0 /= scale;
401  }
402 
403  return 1;
404 }
405 
406 double
408 {
409 #ifdef BINGHAM_ALT
410  if ( tau_0 > 0.0 ) {
411  shearRate = max(shearRate, BINGHAM_MIN_SHEAR_RATE);
412  return ( mu_0 + tau_0 * ( 1. - exp(-this->stressGrowthRate * shearRate) ) / shearRate );
413  } else {
414  // newtonian flow
415  return mu_0;
416  }
417 
418 #else
419  if ( Tau <= tau_c ) {
420  return this->mu_inf;
421  } else {
422  return mu_0 + ( tau_0 / shearRate );
423  }
424 
425 #endif
426 }
427 
428 
429 double
431 {
432  double _val = 0.0;
433  if ( mmode == _2dFlow ) {
434  // _val = 2.0 * ( epsd.at(1) * epsd.at(1) + epsd.at(2) * epsd.at(2) ) + epsd.at(3) * epsd.at(3);
435  _val = 0.5 * ( epsd.at(1) * epsd.at(1) + epsd.at(2) * epsd.at(2) ) + epsd.at(3) * epsd.at(3);
436  } else if ( mmode == _2dAxiFlow ) {
437  _val = 2.0 * ( epsd.at(1) * epsd.at(1) + epsd.at(2) * epsd.at(2) +
438  epsd.at(3) * epsd.at(3) ) + epsd.at(4) * epsd.at(4);
439  } else if ( mmode == _3dFlow ) {
440  _val = 2.0 * ( epsd.at(1) * epsd.at(1) + epsd.at(2) * epsd.at(2) + epsd.at(3) * epsd.at(3) )
441  + epsd.at(4) * epsd.at(4) + epsd.at(5) * epsd.at(5) + epsd.at(6) * epsd.at(6);
442  } else {
443  OOFEM_ERROR("unsupported material mode");
444  }
445 
446  return sqrt(_val);
447 }
448 
449 double
451 {
452  double _val = 0.0;
453  if ( mmode == _2dFlow ) {
454  _val = 0.5 * ( sigd.at(1) * sigd.at(1) + sigd.at(2) * sigd.at(2) + 2.0 * sigd.at(3) * sigd.at(3) );
455  } else if ( mmode == _2dAxiFlow ) {
456  _val = 0.5 * ( sigd.at(1) * sigd.at(1) +
457  sigd.at(2) * sigd.at(2) +
458  sigd.at(3) * sigd.at(3) +
459  2.0 * sigd.at(4) * sigd.at(4) );
460  } else if ( mmode == _3dFlow ) {
461  _val = 0.5 * ( sigd.at(1) * sigd.at(1) + sigd.at(2) * sigd.at(2) + sigd.at(3) * sigd.at(3) +
462  2.0 * sigd.at(4) * sigd.at(4) + 2.0 * sigd.at(5) * sigd.at(5) + 2.0 * sigd.at(6) * sigd.at(6) );
463  } else {
464  OOFEM_ERROR("unsupported material mode");
465  }
466 
467  return sqrt(_val);
468 }
469 
470 void
472 {
473  if ( mmode == _2dFlow ) {
474  //double ekk=(eps.at(1)+eps.at(2))/3.0;
475  double ekk = 0.0;
476 
477  answer = {
478  eps.at(1) - ekk,
479  eps.at(2) - ekk,
480  eps.at(3)
481  };
482  } else if ( mmode == _2dAxiFlow ) {
483  //double ekk=(eps.at(1)+eps.at(2)+eps.at(3))/3.0;
484  double ekk = 0.0;
485 
486  answer = {
487  eps.at(1) - ekk,
488  eps.at(2) - ekk,
489  eps.at(3) - ekk,
490  eps.at(4)
491  };
492  } else if ( mmode == _3dFlow ) {
493  //double ekk=(eps.at(1)+eps.at(2)+eps.at(3))/3.0;
494  double ekk = 0.0;
495 
496  answer = {
497  eps.at(1) - ekk,
498  eps.at(2) - ekk,
499  eps.at(3) - ekk,
500  eps.at(4),
501  eps.at(5),
502  eps.at(6)
503  };
504  } else {
505  OOFEM_ERROR("unsupported material mode");
506  }
507 }
508 
509 
510 void
512  double _nu, MaterialMode mmode)
513 {
514  if ( mmode == _2dFlow ) {
515  answer = {
516  2.0 * _nu * ( deps.at(1) ),
517  2.0 * _nu * ( deps.at(2) ),
518  deps.at(3) * _nu
519  };
520  } else if ( mmode == _2dAxiFlow ) {
521  answer = {
522  answer.at(1) = 2.0 * _nu * ( deps.at(1) ),
523  answer.at(2) = 2.0 * _nu * ( deps.at(2) ),
524  answer.at(3) = 2.0 * _nu * ( deps.at(3) ),
525  answer.at(4) = deps.at(4) * _nu,
526  };
527  } else if ( mmode == _3dFlow ) {
528  answer = {
529  2.0 * _nu * ( deps.at(1) ),
530  2.0 * _nu * ( deps.at(2) ),
531  2.0 * _nu * ( deps.at(3) ),
532  deps.at(4) * _nu,
533  deps.at(5) * _nu,
534  deps.at(6) * _nu
535  };
536  } else {
537  OOFEM_ERROR("unsupported material mode");
538  }
539 }
540 
541 
544 {
545  MaterialMode mmode = gp->giveMaterialMode();
546  int _size = 0;
547 
550 
551  if ( mmode == _2dFlow ) {
552  _size = 3;
553  } else if ( mmode == _2dAxiFlow ) {
554  _size = 4;
555  } else if ( mmode == _3dFlow ) {
556  _size = 6;
557  } else {
558  OOFEM_ERROR("unsupported material mode");
559  }
560 
566 }
567 
568 void
570 // Prints the strains and stresses on the data file.
571 {
572  fprintf(File, " strains ");
573  for ( double e: deviatoricStrainRateVector ) {
574  fprintf( File, " %.4e", e );
575  }
576 
577  fprintf(File, "\n deviatoric stresses");
578  for ( double e: deviatoricStressVector ) {
579  fprintf( File, " %.4e", e );
580  }
581 
582  fprintf(File, "\n status { gamma %e, tau %e }", devStrainMagnitude, devStressMagnitude);
583 
584  fprintf(File, "\n");
585 }
586 
587 void
589 // Performs end-of-step updates.
590 {
592 
596 }
597 
598 
599 void
601 //
602 // initialize record at the begining of new load step
603 //
604 {
606 
610 }
611 
612 
615 //
616 // saves full ms context (saves state variables, that completely describe
617 // current state)
618 // saving the data in dictionary is left to material (yield crit. level).
619 {
620  contextIOResultType iores;
621 
622  if ( ( iores = FluidDynamicMaterialStatus :: saveContext(stream, mode, obj) ) != CIO_OK ) {
623  THROW_CIOERR(iores);
624  }
625 
626  if ( !stream.write(devStrainMagnitude) ) {
628  }
629 
630  if ( !stream.write(devStressMagnitude) ) {
632  }
633 
634  return CIO_OK;
635 }
636 
637 
640 //
641 // restores full material context (saves state variables, that completely describe
642 // current state)
643 //
644 {
645  contextIOResultType iores;
646 
647  if ( ( iores = FluidDynamicMaterialStatus :: restoreContext(stream, mode, obj) ) != CIO_OK ) {
648  THROW_CIOERR(iores);
649  }
650 
651  if ( !stream.read(devStrainMagnitude) ) {
653  }
654 
655  if ( !stream.read(devStressMagnitude) ) {
657  }
658 
659  return CIO_OK;
660 }
661 
662 
663 #if 0
664 void
666 {
667  BinghamFluidMaterial2Status *status = static_cast< BinghamFluidMaterial2Status * >( this->giveStatus(gp) );
668  const FloatArray &epsd = status->giveTempDeviatoricStrainVector();
669  const FloatArray &sigd = status->giveTempDeviatoricStrainVector()
670  for ( int i = 1; i <= nincr; i++ ) {
671  eps.add(eps_i);
672  computeDeviatoricStressVector(tau, gp, eps, tStep);
673  giveDeviatoricStiffnessMatrix(d, TangentStiffness, gp, tStep);
674  tau_t.beProductOf(d, eps_i);
675  tau_t.add(tau_p);
676  //tau.printYourself();
677  //tau_t.printYourself();
678  //d.printYourself();
679  printf( "%e %e %e %e %e %e %e %e %e\n", eps.at(1), eps.at(2), eps.at(3), tau.at(1), tau.at(2), tau.at(3), tau_t.at(1), tau_t.at(2), tau_t.at(3) );
680  tau_p = tau_t;
681  }
682 }
683 #endif
684 } // end namespace oofem
double computeActualViscosity(double Tau, double shearRate)
void setField(int item, InputFieldType id)
MaterialMode giveMaterialMode()
Returns corresponding material mode of receiver.
Definition: gausspoint.h:191
virtual void updateYourself(TimeStep *)
Update equilibrium history variables according to temp-variables.
Definition: matstatus.h:108
#define Viscosity
Definition: matconst.h:79
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Stores receiver state to output stream.
const FloatArray & giveTempDeviatoricStrainVector()
Definition: binghamfluid2.h:93
virtual MaterialStatus * giveStatus(GaussPoint *gp) const
Returns material status of receiver in given integration point.
Definition: material.C:244
FloatArray deviatoricStressVector
Stress vector in reduced form.
GaussPoint * gp
Associated integration point.
Class and object Domain.
Definition: domain.h:115
Abstract base class for all fluid materials.
Domain * domain
Link to domain object, useful for communicating with other FEM components.
Definition: femcmpnn.h:82
virtual void giveDeviatoricStiffnessMatrix(FloatMatrix &answer, MatResponseMode, GaussPoint *gp, TimeStep *tStep)
Computes the deviatoric stiffness; .
void letTempDevStrainMagnitudeBe(double _val)
Definition: binghamfluid2.h:90
virtual void printOutputAt(FILE *file, TimeStep *tStep)
Print receiver&#39;s output to given stream.
double & at(int aKey)
Returns the value of the pair which key is aKey.
Definition: dictionary.C:106
The purpose of DataStream abstract class is to allow to store/restore context to different streams...
Definition: datastream.h:54
#define _IFT_BinghamFluidMaterial2_mu0
Definition: binghamfluid2.h:47
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
int max(int i, int j)
Returns bigger value form two given decimals.
Definition: mathfem.h:71
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Stores receiver state to output stream.
EngngModel * giveEngngModel()
Returns engineering model to which receiver is associated.
Definition: domain.C:433
#define BINGHAM_DEFAULT_STRESS_GROWTH_RATE
Definition: binghamfluid2.h:56
Dictionary propertyDictionary
Property dictionary.
Definition: material.h:105
General IO error.
double giveTempDevStrainMagnitude() const
Definition: binghamfluid2.h:86
virtual void giveInputRecord(DynamicInputRecord &input)
Setups the input record string of receiver.
Definition: material.C:110
BinghamFluidMaterial2Status(int n, Domain *d, GaussPoint *g)
Constructor - creates new BinghamFluidMaterial2Status with number n, belonging to domain d and Integr...
MaterialMode
Type representing material mode of integration point.
Definition: materialmode.h:89
virtual int read(int *data, int count)=0
Reads count integer values into array pointed by data.
MatResponseMode
Describes the character of characteristic material matrix.
virtual void initTempStatus()
Initializes the temporary internal variables, describing the current state according to previously re...
#define THROW_CIOERR(e)
Definition: contextioerr.h:61
#define _IFT_BinghamFluidMaterial2_tau0
Definition: binghamfluid2.h:48
BinghamFluidMaterial2(int n, Domain *d)
Constructor.
Definition: binghamfluid2.C:58
REGISTER_Material_Alt(ConcreteDPM, concreteidp)
void letDeviatoricStressVectorBe(FloatArray v)
Sets the deviatoric stress.
void letTempDevStressMagnitudeBe(double _val)
Definition: binghamfluid2.h:91
void letTempDeviatoricStrainVectorBe(FloatArray v)
Definition: binghamfluid2.h:94
virtual int write(const int *data, int count)=0
Writes count integer values from array pointed by data.
double devStressMagnitude
Magnitude of deviatoric stresses.
Definition: binghamfluid2.h:67
virtual double give(int aProperty, GaussPoint *gp)
Returns the value of material property &#39;aProperty&#39;.
Definition: material.C:52
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Restores the receiver state previously written in stream.
FloatArray deviatoricStrainRateVector
Strain vector in reduced form.
This class implements a transport material status information.
Class representing material status for Bingham material.
Definition: binghamfluid2.h:61
#define OOFEM_ERROR(...)
Definition: error.h:61
virtual double giveVariableScale(VarScaleType varId)
Returns the scale factor for given variable type.
Definition: engngm.h:1087
void __debug(GaussPoint *gp, TimeStep *tStep)
FloatArray temp_deviatoricStrainVector
Deviatoric stresses and strains (reduced form).
Definition: binghamfluid2.h:69
virtual double give(int aProperty, GaussPoint *gp)
Returns the value of material property &#39;aProperty&#39;.
double computeDevStressMagnitude(MaterialMode mmode, const FloatArray &sigd)
void computeDeviatoricStrain(FloatArray &answer, const FloatArray &eps, MaterialMode mmode)
#define _IFT_BinghamFluidMaterial2_stressGrowthRate
Definition: binghamfluid2.h:50
#define _IFT_BinghamFluidMaterial2_muinf
Definition: binghamfluid2.h:49
virtual void updateYourself(TimeStep *tStep)
Update equilibrium history variables according to temp-variables.
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: binghamfluid2.C:68
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
double computeDevStrainMagnitude(MaterialMode mmode, const FloatArray &epsd)
Abstract base class representing a material status information.
Definition: matstatus.h:84
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Restores the receiver state previously written in stream.
Class representing vector of real numbers.
Definition: floatarray.h:82
virtual MaterialStatus * CreateStatus(GaussPoint *gp) const
Creates new copy of associated status and inserts it into given integration point.
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
IRResultType
Type defining the return values of InputRecord reading operations.
Definition: irresulttype.h:47
void computeDeviatoricStress(FloatArray &answer, const FloatArray &deps, double _nu, MaterialMode mmode)
void resize(int rows, int cols)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1358
Class representing the general Input Record.
Definition: inputrecord.h:101
virtual bool giveEquationScalingFlag()
Returns the Equation scaling flag, which is used to indicate that governing equation(s) are scaled...
Definition: engngm.h:1085
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
double devStrainMagnitude
Magnitude of deviatoric strains.
Definition: binghamfluid2.h:65
Class representing the a dynamic Input Record.
long ContextMode
Context mode (mask), defining the type of information written/read to/from context.
Definition: contextmode.h:43
virtual void initTempStatus()
Initializes the temporary internal variables, describing the current state according to previously re...
#define YieldStress
Definition: matconst.h:80
void zero()
Zeroes all coefficient of receiver.
Definition: floatmatrix.C:1326
double tau_0
Yield stress.
Domain * giveDomain() const
Definition: femcmpnn.h:100
int min(int i, int j)
Returns smaller value from two given decimals.
Definition: mathfem.h:59
#define BINGHAM_MIN_SHEAR_RATE
Definition: binghamfluid2.C:51
REGISTER_Material(DummyMaterial)
virtual double giveEffectiveViscosity(GaussPoint *gp, TimeStep *tStep)
Gives the effective viscosity for the given integration point.
Definition: binghamfluid2.C:97
#define IR_GIVE_OPTIONAL_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:78
int giveSize() const
Returns the size of receiver.
Definition: floatarray.h:218
the oofem namespace is to define a context or scope in which all oofem names are defined.
virtual int checkConsistency()
Allows programmer to test some internal data, before computation begins.
virtual void computeDeviatoricStressVector(FloatArray &answer, GaussPoint *gp, const FloatArray &eps, TimeStep *tStep)
Computes the deviatoric stress vector from given strain.
#define IR_GIVE_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:69
virtual void giveInputRecord(DynamicInputRecord &input)
Setups the input record string of receiver.
Definition: binghamfluid2.C:87
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: material.C:89
Class representing integration point in finite element program.
Definition: gausspoint.h:93
Class representing solution step.
Definition: timestep.h:80
void add(const FloatArray &src)
Adds array src to receiver.
Definition: floatarray.C:156
void resize(int s)
Resizes receiver towards requested size.
Definition: floatarray.C:631

This page is part of the OOFEM documentation. Copyright (c) 2011 Borek Patzak
Project e-mail: info@oofem.org
Generated at Tue Jan 2 2018 20:07:27 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011