OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
concretefcm.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 - 2015 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 "concretefcm.h"
36 #include "gausspoint.h"
37 #include "mathfem.h"
38 #include "classfactory.h"
39 #include "contextioerr.h"
40 
41 namespace oofem {
42 REGISTER_Material(ConcreteFCM);
43 
45 {
46  Gf = Ft = 0.;
50 }
51 
54 {
55  IRResultType result; // Required by IR_GIVE_FIELD macro
56 
57  result = FCMMaterial :: initializeFrom(ir);
58  if ( result != IRRT_OK ) {
59  return result;
60  }
61 
63  if ( result != IRRT_OK ) {
64  return result;
65  }
66 
68  if ( result != IRRT_OK ) {
69  return result;
70  }
71 
72  // type of TENSION SOFTENING
73  int softening = 0;
75  switch ( softening ) {
76  case 0:
77  softType = ST_NONE;
78  Ft = 1.e50;
79  Gf = 1.e50;
80  break;
81 
82  case 1:
86  break;
87 
88  case 2:
92  break;
93 
94  case 3:
98  break;
99 
100  case 4:
105 
106  if ( soft_w.at(1) > 0. ) {
107  OOFEM_WARNING("The first entry in soft_w must be 0.");
108  return IRRT_BAD_FORMAT;
109  }
110 
111  if ( soft_w.giveSize() != soft_function_w.giveSize() ) {
112  OOFEM_WARNING("the size of 'soft_w' and 'soft(w)' must be the same");
113  return IRRT_BAD_FORMAT;
114  }
115  break;
116 
117  case 5:
120  this->eps_f = 1.;
123 
124  if ( this->giveCrackSpacing() > 0. ) {
125  OOFEM_ERROR("crack spacing must not be defined in strain-defined materials (does not make sense)");
126  }
127 
128  break;
129 
130  case 6:
135 
136  if ( soft_eps.at(1) > 0. ) {
137  OOFEM_WARNING("The first entry in soft_eps must be 0.");
138  return IRRT_BAD_FORMAT;
139  }
140 
142  OOFEM_WARNING("the size of 'soft_eps' and 'soft(eps)' must be the same");
143  return IRRT_BAD_FORMAT;
144  }
145 
146  if ( this->giveCrackSpacing() > 0. ) {
147  OOFEM_ERROR("crack spacing must not be defined in strain-defined materials (does not make sense)");
148  }
149 
150  break;
151 
152  default:
153  OOFEM_WARNING("Softening type number %d is unknown", softType);
154  return IRRT_BAD_FORMAT;
155  }
156 
157 
158  // type of SHEAR STIFFNESS REDUCTION
159  int shear = 0;
161  switch ( shear ) {
162  case 0:
164 
165  break;
166 
167  case 1:
169  beta = 0.01;
171 
172  if ( beta >= 1. ) {
173  OOFEM_ERROR("Parameter beta %f must be smaller than one", beta);
174  }
175 
176  break;
177 
178  case 2:
179  // todo: resolve problems with convergence when crack is wide open
180  // this does not happen when beta is used instead sf
182  sf = 20.;
184  break;
185 
186  case 3: // w is the maximum reached crack opening, not the current value!
190 
191  if ( beta_w.giveSize() != beta_function.giveSize() ) {
192  OOFEM_WARNING("the size of 'beta_w' and 'beta(w)' must be the same");
193  return IRRT_BAD_FORMAT;
194  }
195  break;
196 
197 
198  default:
199  OOFEM_WARNING("Shear stiffness reduction type number %d is unknown", softType);
200  return IRRT_BAD_FORMAT;
201  }
202 
203 
204  // type of SHEAR STRENGTH CUTOFF
205  shear = 0;
206 
208  switch ( shear ) {
209  case 0:
211  break;
212 
213  case 1:
215  break;
216 
217  case 2:
219 
223 
224  break;
225 
226  default:
227  OOFEM_WARNING("Shear stiffness reduction type number %d is unknown", softType);
228  return IRRT_BAD_FORMAT;
229  }
230 
231  return IRRT_OK;
232 }
233 
234 
235 double
237 //
238 // returns current cracking modulus according to crackStrain for i-th
239 // crackplane
240 //
241 {
242  ConcreteFCMStatus *status = static_cast< ConcreteFCMStatus * >( this->giveStatus(gp) );
243 
244  double Cf = 0.;
245  double ef, emax, c1, c2, Le, Ncr, w;
246 
247  double Ft = this->giveTensileStrength(gp);
248  double Gf = this->giveFractureEnergy(gp);
249  double ec = status->giveTempCrackStrain(i);
250  double E = this->computeOverallElasticStiffness();
251 
252  if ( status->giveTempCrackStatus(i) == pscm_NONE ) {
253  return E * fcm_BIGNUMBER;
254  }
255 
256  // very high stiffness also for compression
257  if ( status->giveTempCrackStatus(i) == pscm_CLOSED ) {
258  return E * fcm_BIGNUMBER;
259  }
260 
261  emax = max(status->giveMaxCrackStrain(i), ec);
262 
263  Le = status->giveCharLength(i);
264 
265  // modification to take into account crack spacing - makes sense only in crack-opening-based materials
266  Ncr = this->giveNumberOfCracksInDirection(gp, i);
267 
268  ec /= Ncr;
269  emax /= Ncr;
270 
271 
272  if ( ( rMode == TangentStiffness ) && ( ec >= emax ) ) { // tangent and NOT unloading
273  if ( softType == ST_NONE ) {
274  OOFEM_ERROR("For unknown reason the tensile strength has been exceeded and cracking has been initiated!");
275  } else if ( softType == ST_Exponential ) {
276  ef = Gf / ( Ft * Le );
277 
278  if ( emax == 0. ) { // just initiated and closing crack
279  Cf = -Ft / ef;
280  } else { // further softening - negative stiffness
281  Cf = -Ft / ef *exp(-ec / ef);
282  }
283  } else if ( softType == ST_Linear ) {
284  // fracturing strain at zero traction
285  ef = 2. * Gf / ( Ft * Le );
286 
287  if ( ( ec >= ef ) || ( emax >= ef ) ) {
288  // fully open crack - no stiffness
289  Cf = 0.;
290  } else { // further softening - negative stiffness
291  Cf = -Ft / ef;
292  }
293  } else if ( softType == ST_Hordijk ) {
294  c1 = 3.;
295  c2 = 6.93;
296 
297  ef = 5.14 * Gf / Ft / Le;
298 
299  if ( ( ec >= ef ) || ( emax >= ef ) ) {
300  Cf = 0.;
301  } else if ( emax == 0. ) { // just initiated and closing crack
302  Cf = Ft * ( -c2 / ef - ( exp(-c2) * ( c1 * c1 * c1 + 1. ) ) / ef );
303  } else { // further softening - negative stiffness
304  Cf = Ft * ( ( 3. * c1 * c1 * c1 * ec * ec * exp(-( c2 * ec ) / ef) ) / ( ef * ef * ef ) - ( c2 * exp(-( c2 * ec ) / ef) * ( c1 * c1 * c1 * ec * ec * ec + ef * ef * ef ) ) / ( ef * ef * ef * ef ) - ( exp(-c2) * ( c1 * c1 * c1 + 1. ) ) / ef );
305  }
306  } else if ( softType == ST_UserDefinedCrack ) {
307  w = ec * Le;
308 
309  if ( w > soft_w.at(soft_w.giveSize() ) ) {
310  Cf = 0.;
311 #ifdef DEBUG
312  OOFEM_WARNING("Crack width is larger than the last user-defined value in the traction-opening law.");
313 #endif
314  } else if ( emax == 0. ) {
315  Cf = this->Ft * ( soft_function_w.at(2) - soft_function_w.at(1) ) / ( ( soft_w.at(2) - soft_w.at(1) ) / Le );
316  } else { // softening
317  for ( int i = 1; i <= soft_w.giveSize(); i++ ) {
318  if ( ( w - soft_w.at(i) ) < fcm_SMALL_STRAIN ) {
319  if ( i == 1 ) {
320  Cf = this->Ft * ( soft_function_w.at(2) - soft_function_w.at(1) ) / ( ( soft_w.at(2) - soft_w.at(1) ) / Le );
321  break;
322  } else {
323  Cf = this->Ft * ( soft_function_w.at(i) - soft_function_w.at(i - 1) ) / ( ( soft_w.at(i) - soft_w.at(i - 1) ) / Le );
324  break;
325  }
326  }
327  }
328  }
329  } else if ( softType == ST_LinearHardeningStrain ) {
330  if ( ( ec >= this->eps_f ) || ( emax >= this->eps_f ) ) {
331  Cf = 0.;
332  } else {
333  // hardening
334  Cf = this->H;
335  }
336  } else if ( softType == ST_UserDefinedStrain ) {
337  if ( ec > soft_eps.at(soft_eps.giveSize() ) ) {
338  Cf = 0.;
339 #ifdef DEBUG
340  OOFEM_WARNING("Crack strain is larger than the last user-defined value in the stress-crack-strain law.");
341 #endif
342  } else if ( emax == 0. ) {
343  Cf = this->Ft * ( soft_function_eps.at(2) - soft_function_eps.at(1) ) / ( soft_eps.at(2) - soft_eps.at(1) );
344  } else { // softening
345  for ( int i = 1; i <= soft_eps.giveSize(); i++ ) {
346  if ( ( ec - soft_eps.at(i) ) < fcm_SMALL_STRAIN ) {
347  if ( i == 1 ) {
348  Cf = this->Ft * ( soft_function_eps.at(2) - soft_function_eps.at(1) ) / ( soft_eps.at(2) - soft_eps.at(1) );
349  break;
350  } else {
351  Cf = this->Ft * ( soft_function_eps.at(i) - soft_function_eps.at(i - 1) ) / ( soft_eps.at(i) - soft_eps.at(i - 1) );
352  break;
353  }
354  }
355  }
356  }
357  } else {
358  OOFEM_ERROR("Unknown Softening Mode");
359  }
360  } else { // secant stiffness OR tangent stiffness and unloading
361  if ( emax <= 0. ) { // just initiated crack that wants to close
362  return E * fcm_BIGNUMBER;
363  }
364 
365  Cf = ConcreteFCM :: giveNormalCrackingStress(gp, emax * Ncr, i); // gets stress
366  Cf /= ( emax ); // converted into stiffness
367  }
368 
369  // to take into account crack-spacing
370  return Cf / Ncr;
371 }
372 
373 double
375 //
376 // returns receivers Normal Stress in crack i for given cracking strain
377 //
378 {
379  ConcreteFCMStatus *status = static_cast< ConcreteFCMStatus * >( this->giveStatus(gp) );
380 
381  double traction;
382  double Le, ef, emax, w, E, c1, c2;
383  double Gf = this->giveFractureEnergy(gp);
384  double Ft = this->giveTensileStrength(gp);
385 
386  if ( status->giveTempCrackStatus(i) == pscm_CLOSED ) {
387  E = this->computeOverallElasticStiffness();
388  return E * ec;
389  }
390 
391  Le = status->giveCharLength(i);
392  emax = max(status->giveMaxCrackStrain(i), ec);
393 
394  // to take into account crack spacing
395  ec /= this->giveNumberOfCracksInDirection(gp, i);
396  emax /= this->giveNumberOfCracksInDirection(gp, i);
397 
398 
399  if ( softType == ST_NONE ) {
400  OOFEM_ERROR("For unknown reason the tensile strength has been exceeded and cracking has been initiated!");
401  } else if ( softType == ST_Exponential ) {
402  ef = Gf / ( Ft * Le );
403  if ( ec >= emax ) {
404  // further softening
405  traction = Ft * exp(-ec / ef);
406  } else {
407  // crack closing
408  // or unloading or reloading regime
409  traction = Ft * ec / emax *exp(-emax / ef);
410  }
411  } else if ( softType == ST_Linear ) {
412  // fracturing strain at zero traction
413  ef = 2. * Gf / ( Ft * Le );
414 
415  if ( ( ec >= ef ) || ( emax >= ef ) ) {
416  // fully open crack - no stiffness
417  traction = 0.;
418  } else if ( ec >= emax ) {
419  // further softening
420  traction = Ft - Ft * ec / ef;
421  } else if ( ec <= 0. ) {
422  traction = 0.;
423  } else {
424  traction = Ft * ec * ( ef - emax ) / ( emax * ef );
425  }
426  } else if ( softType == ST_Hordijk ) {
427  c1 = 3.;
428  c2 = 6.93;
429 
430  ef = 5.14 * Gf / Ft / Le;
431 
432  if ( ( ec >= ef ) || ( emax >= ef ) ) {
433  // fully open crack - no stiffness
434  traction = 0.;
435  } else if ( ec >= emax ) {
436  // further softening
437  traction = Ft * ( ( 1. + pow( ( c1 * ec / ef ), 3. ) ) * exp(-c2 * ec / ef) - ec / ef * ( 1. + c1 * c1 * c1 ) * exp(-c2) );
438  } else if ( ec <= 0. ) {
439  traction = 0.;
440  } else {
441  traction = Ft * ec / emax * ( ( 1. + pow( ( c1 * emax / ef ), 3. ) ) * exp(-c2 * emax / ef) - emax / ef * ( 1. + c1 * c1 * c1 ) * exp(-c2) );
442  }
443  } else if ( softType == ST_UserDefinedCrack ) {
444  w = ec * Le;
445 
446  traction = 0.;
447 
448  if ( w > soft_w.at(soft_w.giveSize() ) ) {
449  traction = 0.;
450 
451 #ifdef DEBUG
452  OOFEM_WARNING("Crack width is larger than the last user-defined value in the traction-opening law.");
453 #endif
454 
455  } else if ( ec >= emax ) { // softening
456  for ( int i = 1; i <= soft_w.giveSize(); i++ ) {
457  if ( ( w - soft_w.at(i) ) < fcm_SMALL_STRAIN ) {
458  if ( i == 1 ) { // eps \approx 0 bound
459  traction = soft_function_w.at(i) * this->Ft;
460  break;
461  } else {
462  traction = soft_function_w.at(i - 1) + ( soft_function_w.at(i) - soft_function_w.at(i - 1) ) / ( soft_w.at(i) - soft_w.at(i - 1) ) * ( w - soft_w.at(i - 1) );
463  traction *= this->Ft;
464  break;
465  }
466  }
467  }
468  } else if ( ec <= 0. ) { // closing
469  traction = 0.;
470  } else { //unlo-relo
471  w = emax * Le; // lower stiffness
472 
473  for ( int i = 1; i <= soft_w.giveSize(); i++ ) {
474  if ( ( w - soft_w.at(i) ) < fcm_SMALL_STRAIN ) {
475  if ( i == 1 ) { // eps \approx 0
476  traction = soft_function_w.at(i) * this->Ft * ec / emax;
477  break;
478  } else {
479  traction = soft_function_w.at(i - 1) + ( soft_function_w.at(i) - soft_function_w.at(i - 1) ) / ( soft_w.at(i) - soft_w.at(i - 1) ) * ( w - soft_w.at(i - 1) );
480  traction *= this->Ft * ec / emax;
481  break;
482  }
483  }
484  }
485  }
486  } else if ( softType == ST_LinearHardeningStrain ) {
487  if ( ( ec >= this->eps_f ) || ( emax >= this->eps_f ) ) {
488  traction = 0.;
489  } else if ( emax == 0. ) {
490  traction = this->Ft;
491  } else {
492  // unloading or reloading regime
493  // emax * Le = crack width
494  // traction = ( this->Ft + emax * Le * this->H ) * ec / emax;
495  // written in terms of strain and not crack width
496  traction = ( this->Ft + emax * this->H ) * ec / emax;
497  }
498  } else if ( softType == ST_UserDefinedStrain ) {
499  if ( emax > soft_eps.at(soft_eps.giveSize() ) ) {
500  traction = 0.;
501 #ifdef DEBUG
502  OOFEM_WARNING("Cracking strain is larger than the last user-defined value in the traction-cracking-strain law.");
503 #endif
504  } else if ( ec >= emax ) { // softening
505  for ( int i = 1; i <= soft_eps.giveSize(); i++ ) {
506  if ( ( ec - soft_eps.at(i) ) < fcm_SMALL_STRAIN ) {
507  if ( i == 1 ) { // eps \approx 0 bound
508  traction = soft_function_eps.at(i) * this->Ft;
509  break;
510  } else {
511  traction = soft_function_eps.at(i - 1) + ( soft_function_eps.at(i) - soft_function_eps.at(i - 1) ) / ( soft_eps.at(i) - soft_eps.at(i - 1) ) * ( ec - soft_eps.at(i - 1) );
512  traction *= this->Ft;
513  break;
514  }
515  }
516  }
517  } else if ( ec <= 0. ) { // closing
518  traction = 0.;
519  } else { //unlo-relo
520  for ( int i = 1; i <= soft_eps.giveSize(); i++ ) {
521  if ( ( emax - soft_eps.at(i) ) < fcm_SMALL_STRAIN ) {
522  if ( i == 1 ) { // eps \approx 0
523  traction = soft_function_eps.at(i) * this->Ft * ec / emax;
524  break;
525  } else {
526  traction = soft_function_eps.at(i - 1) + ( soft_function_eps.at(i) - soft_function_eps.at(i - 1) ) / ( soft_eps.at(i) - soft_eps.at(i - 1) ) * ( emax - soft_eps.at(i - 1) );
527  traction *= this->Ft * ec / emax;
528  break;
529  }
530  }
531  }
532  }
533  } else {
534  OOFEM_ERROR("Unknown Softening Mode");
535  }
536 
537  return traction;
538 }
539 
540 double
542 
543 {
544  double G, Geff, D2tot, N;
545  int crackA, crackB;
546 
548 
549  if ( this->isIntactForShear(gp, shearDirection) ) {
550  Geff = G;
551  } else {
552  if ( this->shearType == SHR_NONE ) {
553  Geff = G;
554  } else if ( this->shearType == SHR_Const_ShearRetFactor ) {
555  if ( multipleCrackShear ) {
556  if ( shearDirection == 4 ) {
557  crackA = 2;
558  crackB = 3;
559  } else if ( shearDirection == 5 ) {
560  crackA = 1;
561  crackB = 3;
562  } else if ( shearDirection == 6 ) {
563  crackA = 1;
564  crackB = 2;
565  } else {
566  OOFEM_ERROR("Unexpected value of index i (4, 5, 6 permitted only)");
567  crackA = crackB = 0; // happy compiler
568  }
569 
570  // total number of parallel cracks in both directions
571  N = this->giveNumberOfCracksInDirection(gp, crackA) + this->giveNumberOfCracksInDirection(gp, crackB);
572  } else {
573  N = this->giveNumberOfCracksForShearDirection(gp, shearDirection);
574  }
575 
576  Geff = G * this->beta / ( N - this->beta * ( N - 1 ) ); // for N = 1... Geff = G * beta
577  } else if ( this->shearType == SHR_Const_ShearFactorCoeff ) {
578  // contributions from shear terms in: De - De (De + Dcr)^-1 De
579  D2tot = this->computeTotalD2Modulus(gp, shearDirection);
580  // number of parallel cracks already taken into account
581  Geff = G * D2tot / ( G + D2tot );
582  } else if ( this->shearType == SHR_UserDefined_ShearRetFactor ) {
583  // contributions from shear terms in: De - De (De + Dcr)^-1 De
584  D2tot = this->computeTotalD2Modulus(gp, shearDirection);
585  // number of parallel cracks already taken into account
586  Geff = G * D2tot / ( G + D2tot );
587  } else {
588  OOFEM_ERROR("Unknown Shear Mode");
589  Geff = 0.; // happy compiler
590  }
591  }
592 
593  return Geff;
594 }
595 
596 
597 double
599 {
600  double D2, N, G, E, w, beta_m;
601 
602  if ( icrack >= 4 ) {
603  OOFEM_ERROR("Unexpected crack number");
604  }
605 
606  N = this->giveNumberOfCracksInDirection(gp, icrack);
607 
608  if ( this->isIntact(gp, icrack) || ( this->shearType == SHR_NONE ) ) {
609  E = this->computeOverallElasticStiffness();
610  D2 = E * fcm_BIGNUMBER;
611  } else {
614  D2 = G * this->beta / ( 1. - this->beta );
615  D2 /= N;
616  } else if ( shearType == SHR_Const_ShearFactorCoeff ) {
617  // for the number of parallel cracks is taken care of in "giveCrackingModulus"
618  // however, it is essential to reduce the final modulus due to the presence of multiple cracks
619  D2 = this->sf * this->giveCrackingModulus(SecantStiffness, gp, icrack);
620  D2 /= N;
621  } else if ( this->shearType == SHR_UserDefined_ShearRetFactor ) {
622  // for the number of parallel cracks is taken care of in the evaluation of normalcrackopening
623  // however, it is essential to reduce the final modulus due to the presence of multiple cracks
624 
626  w = max( this->computeNormalCrackOpening(gp, icrack), this->computeMaxNormalCrackOpening(gp, icrack) );
627 
629 
630  beta_m = 0.;
631  for ( int i = 1; i <= this->beta_w.giveSize(); i++ ) {
632  if ( ( w - this->beta_w.at(i) ) < fcm_SMALL_STRAIN ) {
633  if ( i == 1 ) {
634  beta_m = this->beta_function.at(i);
635  } else {
636  beta_m = this->beta_function.at(i - 1) + ( this->beta_function.at(i) - this->beta_function.at(i - 1) ) / ( this->beta_w.at(i) - this->beta_w.at(i - 1) ) * ( w - this->beta_w.at(i - 1) );
637  }
638 
639  beta_m = min(1., beta_m);
640  break;
641  }
642  }
643 
644  if ( beta_m >= 1. ) {
645  D2 = G * fcm_BIGNUMBER;
646  } else {
647  D2 = G * beta_m / ( 1. - beta_m );
648  D2 /= N;
649  }
650  } else {
651  D2 = 0;
652  OOFEM_ERROR("Unknown Softening Mode");
653  }
654  }
655 
656  return D2;
657 }
658 
659 
660 void
662 {
663  ConcreteFCMStatus *status = static_cast< ConcreteFCMStatus * >( this->giveStatus(gp) );
664 
665  double E, Gf, Ft, Le;
666  double ef = 0.;
667  double slope;
668  double efail, c1, c2;
669 
670  Le = status->giveCharLength(i);
671 
672  Gf = this->giveFractureEnergy(gp);
673  Ft = this->giveTensileStrength(gp);
674 
675  E = this->computeOverallElasticStiffness();
676 
677  if ( softType == ST_NONE ) {
678  OOFEM_ERROR("For unknown reason the tensile strength has been exceeded and cracking has been initiated!");
679  } else if ( softType == ST_Exponential ) {
680  ef = Gf / ( Ft * Le );
681  } else if ( softType == ST_Linear ) {
682  ef = 2. * Gf / ( Ft * Le );
683  } else if ( softType == ST_Hordijk ) {
684  efail = 5.14 * Gf / Ft / Le;
685  c1 = 3.;
686  c2 = 6.93;
687  // ef = ft / slope of the softening branch at the peak stress
688  ef = 1. / ( c2 / efail + ( exp(-c2) * ( c1 * c1 * c1 + 1. ) ) / efail );
689  } else if ( softType == ST_UserDefinedCrack ) {
690  ef = 1.e20;
691 
692  for ( int i = 1; i < soft_w.giveSize(); i++ ) {
693  slope = ( soft_function_w.at(i + 1) - soft_function_w.at(i) ) / ( ( soft_w.at(i + 1) - soft_w.at(i) ) / Le );
694  if ( slope < 0. ) {
695  ef = min(fabs(1. / slope), ef);
696  }
697  }
698  } else if ( softType == ST_LinearHardeningStrain ) {
699  return;
700  } else if ( softType == ST_UserDefinedStrain ) {
701  ef = 1.e20;
702 
703  for ( int i = 1; i < soft_eps.giveSize(); i++ ) {
704  slope = ( soft_function_eps.at(i + 1) - soft_function_eps.at(i) ) / ( soft_eps.at(i + 1) - soft_eps.at(i) );
705  if ( slope < 0. ) {
706  ef = min(fabs(1. / slope), ef);
707  }
708  }
709  } else {
710  OOFEM_ERROR("Unknown Softening Mode");
711  }
712 
713  if ( ef <= Ft / E ) {
714  OOFEM_ERROR("ef %e < e0 %e, this leads to material snapback in element %d, characteristic length %f", ef, Ft / E, gp->giveElement()->giveNumber(), Le);
715  }
716 }
717 
718 
719 double
721  int dir_1, dir_2;
722  double maxTau, crackOpening, wmax, scale;
723 
724  if ( this->isIntactForShear(gp, i) ) {
725  return Ft * fcm_BIGNUMBER;
726  }
727 
728  if ( this->shearStrengthType == SHS_NONE ) {
729  maxTau = Ft * fcm_BIGNUMBER;
730  } else if ( this->shearStrengthType == SHS_Const_Ft ) {
731  maxTau = Ft;
732  } else if ( this->shearStrengthType == SHS_Collins_Interlock ) {
733  if ( i == 4 ) { // y-z
734  dir_1 = 2;
735  dir_2 = 3;
736  } else if ( i == 5 ) { // x-z
737  dir_1 = 1;
738  dir_2 = 3;
739  } else if ( i == 6 ) { // x-y
740  dir_1 = 1;
741  dir_2 = 2;
742  } else {
743  OOFEM_ERROR("Unexpected number for shear stress (must be either 4, 5 or 6).");
744  dir_1 = dir_2 = 0; // happy compiler
745  }
746 
747  // from temporary strain
748  crackOpening = max( this->computeNormalCrackOpening(gp, dir_1), this->computeNormalCrackOpening(gp, dir_2) );
749 
750  // from the entire history
751  wmax = max( this->computeMaxNormalCrackOpening(gp, dir_1), this->computeMaxNormalCrackOpening(gp, dir_2) );
752 
753  crackOpening = max(wmax, crackOpening);
754 
755  scale = 1000. / lengthScale;
756  maxTau = 0.18 * sqrt(this->fc) / ( 0.31 + 24. * crackOpening * scale / ( this->ag * scale + 16. ) );
757  } else {
758  OOFEM_ERROR("Unexpected shearStrengthType");
759  maxTau = 0.; // happy compiler
760  }
761 
762  return maxTau;
763 }
764 
765 
766 double
767 ConcreteFCM :: give(int aProperty, GaussPoint *gp)
768 {
769  double answer;
770  if ( RandomMaterialExtensionInterface :: give(aProperty, gp, answer) ) {
771  return answer;
772  } else if ( aProperty == gf_ID ) {
773  return this->Gf;
774  } else if ( aProperty == ft_strength ) {
775  return this->Ft;
776  } else {
777  return FCMMaterial :: give(aProperty, gp);
778  }
779 }
780 
781 int
783 {
784  return FCMMaterial :: giveIPValue(answer, gp, type, tStep);
785 }
786 
787 
788 
791 {
792  MaterialStatus *status = static_cast< MaterialStatus * >( gp->giveMaterialStatus() );
793  if ( status == NULL ) {
794  // create a new one
795  status = this->CreateStatus(gp);
796 
797  if ( status != NULL ) {
798  gp->setMaterialStatus( status, this->giveNumber() );
799  this->_generateStatusVariables(gp);
800  }
801  }
802 
803  return status;
804 }
805 
806 
807 
808 
809 
811 // CONCRETE FCM STATUS ///
813 
816 {}
817 
818 
820 { }
821 
822 
823 
824 void
826 {
828 }
829 
830 
831 
832 
833 
834 void
836 //
837 // initializes temp variables according to variables form previous equlibrium state.
838 //
839 {
841 }
842 
843 
844 
845 void
847 //
848 // updates variables (nonTemp variables describing situation at previous equilibrium state)
849 // after a new equilibrium state has been reached
850 // temporary variables correspond to newly reched equilibrium.
851 //
852 {
854 }
855 
856 Interface *
858 {
860  return static_cast< RandomMaterialStatusExtensionInterface * >(this);
861  } else {
862  return NULL;
863  }
864 }
865 
866 
867 
870 //
871 // saves full information stored in this Status
872 // no temp variables stored
873 //
874 {
875  contextIOResultType iores;
876 
877  // save parent class status
878  if ( ( iores = FCMMaterialStatus :: saveContext(stream, mode, obj) ) != CIO_OK ) {
879  THROW_CIOERR(iores);
880  }
881 
882  return CIO_OK;
883 }
884 
887 //
888 // restores full information stored in stream to this Status
889 //
890 {
891  contextIOResultType iores;
892 
893  // read parent class status
894  if ( ( iores = FCMMaterialStatus :: restoreContext(stream, mode, obj) ) != CIO_OK ) {
895  THROW_CIOERR(iores);
896  }
897 
898  return CIO_OK; // return succes
899 }
900 } // end namespace oofem
FloatArray soft_eps
user-defined softening (traction-strain)
Definition: concretefcm.h:143
InternalStateType
Type representing the physical meaning of element or constitutive model internal variable.
#define fcm_BIGNUMBER
Definition: fcm.h:61
bool multipleCrackShear
if true = takes shear compliance of all cracks, false = only dominant crack contribution, default value is false
Definition: fcm.h:274
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: fcm.C:1424
virtual double computeOverallElasticStiffness(void)
returns overall Young&#39;s modulus
Definition: fcm.h:344
double H
hardening modulus
Definition: concretefcm.h:148
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in Reduced form.
Definition: fcm.C:1527
Class and object Domain.
Definition: domain.h:115
#define _IFT_ConcreteFCM_softType
Definition: concretefcm.h:44
IntegrationPointStatus * setMaterialStatus(IntegrationPointStatus *ptr, int n)
Sets Material status managed by receiver.
Definition: gausspoint.h:213
virtual bool isIntact(GaussPoint *gp, int icrack)
returns true for closed or no crack (i = 1, 2, 3)
Definition: fcm.C:720
This class implements associated Material Status to FCMMaterial (fixed crack material).
Definition: fcm.h:68
#define ft_strength
Definition: matconst.h:91
#define _IFT_ConcreteFCM_shearStrengthType
Definition: concretefcm.h:46
The purpose of DataStream abstract class is to allow to store/restore context to different streams...
Definition: datastream.h:54
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj=NULL)
saves current context(state) into stream
Definition: fcm.C:1951
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
#define pscm_CLOSED
Definition: fcm.h:58
virtual double giveCrackSpacing(void)
returns either user-provided value of crack spacing or a value computed from composition ...
Definition: fcm.C:1475
This class implements a Fixed Crack Model for fracture (after initiation the crack directions cannot ...
Definition: fcm.h:197
#define _IFT_ConcreteFCM_gf
Definition: concretefcm.h:47
This class manages the status of ConcreteFCM.
Definition: concretefcm.h:68
#define pscm_NONE
Definition: fcm.h:54
virtual double giveNormalCrackingStress(GaussPoint *gp, double eps_cr, int i)
computes normal stress associated with i-th crack direction
Definition: concretefcm.C:374
virtual Interface * giveInterface(InterfaceType it)
Interface requesting service.
Definition: concretefcm.C:857
Element * giveElement()
Returns corresponding element to receiver.
Definition: gausspoint.h:188
double giveCharLength(int icrack) const
returns characteristic length associated with i-th crack direction
Definition: fcm.h:158
#define _IFT_ConcreteFCM_soft_w
Definition: concretefcm.h:54
#define _IFT_ConcreteFCM_soft_function_w
Definition: concretefcm.h:55
#define _IFT_ConcreteFCM_H
Definition: concretefcm.h:60
MatResponseMode
Describes the character of characteristic material matrix.
virtual double computeMaxNormalCrackOpening(GaussPoint *gp, int i)
uses maximum equilibrated cracking strain and characteristic length to obtain the maximum reached cra...
Definition: fcm.C:788
virtual double computeD2ModulusForCrack(GaussPoint *gp, int icrack)
shear modulus for a given crack plane (1, 2, 3)
Definition: concretefcm.C:598
#define _IFT_ConcreteFCM_soft_function_eps
Definition: concretefcm.h:57
#define THROW_CIOERR(e)
Definition: contextioerr.h:61
ConcreteFCM(int n, Domain *d)
Definition: concretefcm.C:44
ShearStrengthType shearStrengthType
Definition: concretefcm.h:186
virtual double computeTotalD2Modulus(GaussPoint *gp, int i)
shear modulus for a given shear direction (4, 5, 6)
Definition: fcm.C:1378
SofteningType softType
Definition: concretefcm.h:178
void _generateStatusVariables(GaussPoint *) const
Sets up (generates) the variables identified in randVariables array using generators given in randomV...
virtual double give(int aProperty, GaussPoint *gp)
Returns the value of material property &#39;aProperty&#39;.
Definition: fcm.C:1467
virtual void initTempStatus()
initializes temporary status
Definition: fcm.C:1891
virtual void checkSnapBack(GaussPoint *gp, int crack)
checks possible snap-back
Definition: concretefcm.C:661
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj=NULL)
saves current context(state) into stream
Definition: concretefcm.C:869
double giveTempCrackStrain(int icrack) const
returns i-th component of the crack strain vector (temporary)
Definition: fcm.h:132
ConcreteFCMStatus(int n, Domain *d, GaussPoint *g)
Definition: concretefcm.C:814
double beta
shear retention factor
Definition: concretefcm.h:128
virtual double giveFractureEnergy(GaussPoint *gp)
returns fracture energy (can be random)
Definition: concretefcm.h:156
#define E(p)
Definition: mdm.C:368
#define _IFT_ConcreteFCM_ft
Definition: concretefcm.h:48
This class implements an isotropic linear elastic material in a finite element problem.
#define _IFT_ConcreteFCM_beta
Definition: concretefcm.h:49
#define OOFEM_ERROR(...)
Definition: error.h:61
#define gf_ID
Definition: matconst.h:85
#define _IFT_ConcreteFCM_fc
Definition: concretefcm.h:51
#define _IFT_ConcreteFCM_lengthScale
Definition: concretefcm.h:53
#define _IFT_ConcreteFCM_shearType
Definition: concretefcm.h:45
Abstract base class for all random materials.
virtual MaterialStatus * CreateStatus(GaussPoint *gp) const
Creates new copy of associated status and inserts it into given integration point.
Definition: concretefcm.h:112
FloatArray beta_w
user-defined shear retention factor (with respect to crack opening)
Definition: concretefcm.h:145
#define _IFT_ConcreteFCM_soft_eps
Definition: concretefcm.h:56
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj=NULL)
restores context(state) from stream
Definition: fcm.C:2010
double ag
Collins&#39; aggregate interlock: aggregate diameter in appropriate units (same as FE mesh) ...
Definition: concretefcm.h:136
#define _IFT_ConcreteFCM_sf
Definition: concretefcm.h:50
#define N(p, q)
Definition: mdm.C:367
bool give(int key, GaussPoint *gp, double &value)
Returns the property in associated status of given integration point if defined.
ShearRetentionType shearType
Definition: concretefcm.h:182
Abstract base class for all random constitutive model statuses.
virtual double computeEffectiveShearModulus(GaussPoint *gp, int i)
returns Geff which is necessary in the global stiffness matrix
Definition: concretefcm.C:541
double eps_f
strain at failure
Definition: concretefcm.h:150
double Ft
Tensile strenght.
Definition: concretefcm.h:126
virtual MaterialStatus * giveStatus(GaussPoint *gp) const
Returns material status of receiver in given integration point.
Definition: concretefcm.C:790
Abstract base class representing a material status information.
Definition: matstatus.h:84
Class representing vector of real numbers.
Definition: floatarray.h:82
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
virtual double maxShearStress(GaussPoint *gp, int i)
computes the maximum value of the shear stress; if the shear stress exceeds this value, it is cropped
Definition: concretefcm.C:720
virtual double giveNumberOfCracksForShearDirection(GaussPoint *gp, int i)
returns number of cracks for given shear direction (i = 4, 5, 6) which is treated as the maximum of t...
Definition: fcm.C:1501
const IntArray & giveTempCrackStatus()
returns vector of temporary crack statuses
Definition: fcm.h:117
IRResultType
Type defining the return values of InputRecord reading operations.
Definition: irresulttype.h:47
virtual void printOutputAt(FILE *file, TimeStep *tStep)
Writes information into the output file.
Definition: concretefcm.C:825
virtual void initTempStatus()
initializes temporary status
Definition: concretefcm.C:835
IntegrationPointStatus * giveMaterialStatus()
Returns reference to associated material status (NULL if not defined).
Definition: gausspoint.h:205
double lengthScale
Collins&#39; aggregate interlock: 1 for meter, 1000 for analysis in mm.
Definition: concretefcm.h:138
#define _IFT_ConcreteFCM_beta_function
Definition: concretefcm.h:59
Class representing the general Input Record.
Definition: inputrecord.h:101
#define _IFT_ConcreteFCM_ag
Definition: concretefcm.h:52
IsotropicLinearElasticMaterial * linearElasticMaterial
Definition: fcm.h:200
double giveMaxCrackStrain(int icrack)
returns maximum crack strain for the i-th crack (equilibrated value)
Definition: fcm.h:107
Class Interface.
Definition: interface.h:82
#define _IFT_ConcreteFCM_eps_f
Definition: concretefcm.h:61
virtual bool isIntactForShear(GaussPoint *gp, int i)
returns true for closed or no cracks associated to given shear direction (i = 4, 5, 6)
Definition: fcm.C:738
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj=NULL)
restores context(state) from stream
Definition: concretefcm.C:886
virtual double giveNumberOfCracksInDirection(GaussPoint *gp, int iCrack)
returns number of fictiotious parallel cracks in the direction of i-th crack
Definition: fcm.C:1482
virtual void printOutputAt(FILE *file, TimeStep *tStep)
prints the output into the output file
Definition: fcm.C:1784
virtual void updateYourself(TimeStep *tStep)
replaces equilibrated values with temporary values
Definition: fcm.C:1906
long ContextMode
Context mode (mask), defining the type of information written/read to/from context.
Definition: contextmode.h:43
InterfaceType
Enumerative type, used to identify interface type.
Definition: interfacetype.h:43
int min(int i, int j)
Returns smaller value from two given decimals.
Definition: mathfem.h:59
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in Reduced form.
Definition: concretefcm.C:782
REGISTER_Material(DummyMaterial)
#define IR_GIVE_OPTIONAL_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:78
FloatArray beta_function
Definition: concretefcm.h:145
int giveSize() const
Returns the size of receiver.
Definition: floatarray.h:218
virtual double computeNormalCrackOpening(GaussPoint *gp, int i)
uses temporary cracking strain and characteristic length to obtain the crack opening ...
Definition: fcm.C:769
virtual void updateYourself(TimeStep *tStep)
replaces equilibrated values with temporary values
Definition: concretefcm.C:846
the oofem namespace is to define a context or scope in which all oofem names are defined.
#define IR_GIVE_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:69
int giveNumber() const
Definition: femcmpnn.h:107
virtual double give(int aProperty, GaussPoint *gp)
Returns the value of material property &#39;aProperty&#39;.
Definition: concretefcm.C:767
#define fcm_SMALL_STRAIN
Definition: fcm.h:60
Class representing integration point in finite element program.
Definition: gausspoint.h:93
#define OOFEM_WARNING(...)
Definition: error.h:62
IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Class representing solution step.
Definition: timestep.h:80
virtual double giveCrackingModulus(MatResponseMode rMode, GaussPoint *gp, int i)
returns stiffness in the normal direction of the i-th crack
Definition: concretefcm.C:236
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: concretefcm.C:53
double sf
shear factor
Definition: concretefcm.h:130
FloatArray soft_w
user-defined softening (traction-COD)
Definition: concretefcm.h:141
double fc
Collins&#39; aggregate interlock: compressive strength in MPa.
Definition: concretefcm.h:134
FloatArray soft_function_eps
Definition: concretefcm.h:143
FloatArray soft_function_w
Definition: concretefcm.h:141
#define _IFT_ConcreteFCM_beta_w
Definition: concretefcm.h:58
virtual double computeOverallElasticShearModulus(void)
returns overall shear modulus
Definition: fcm.h:347
double Gf
Fracture energy.
Definition: concretefcm.h:124
virtual double giveTensileStrength(GaussPoint *gp)
returns tensile strength (can be random)
Definition: concretefcm.h:153

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:28 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011