OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
abaqususermaterial.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 "abaqususermaterial.h"
36 #include "gausspoint.h"
37 #include "classfactory.h"
38 #include "dynamicinputrecord.h"
39 
40 #ifdef _WIN32 //_MSC_VER and __MINGW32__ included
41  #include <Windows.h>
42 #else
43  #include <dlfcn.h>
44 #endif
45 
46 #include <cstring>
47 
48 namespace oofem {
49 REGISTER_Material(AbaqusUserMaterial);
50 
52 
54  StructuralMaterial(n, d),
55  umatobj(NULL), umat(NULL),
56  mStressInterpretation(0),
57  mUseNumericalTangent(false),
58  mPerturbation(1.0e-7)
59 { }
60 
62 {
63 #ifdef _WIN32
64  if ( this->umatobj ) {
65  FreeLibrary( ( HMODULE ) this->umatobj );
66  }
67 #else
68  if ( this->umatobj ) {
69  dlclose(this->umatobj);
70  }
71 
72 #endif
73 }
74 
76 {
77  IRResultType result;
78  std :: string umatname;
79 
81  if ( result != IRRT_OK ) return result;
82 
86  umatname = "umat";
88  strncpy(this->cmname, umatname.c_str(), 80);
90 
91 #ifdef _WIN32
92  this->umatobj = ( void * ) LoadLibrary( filename.c_str() );
94  if ( !this->umatobj ) {
95  DWORD dlresult = GetLastError(); //works for MinGW 32bit and MSVC
96  OOFEM_ERROR("Couldn't load \"%s\",\nerror code = %d", filename.c_str(), dlresult);
97  }
98 
99  // * ( void ** )( & this->umat ) = GetProcAddress( ( HMODULE ) this->umatobj, "umat_" );
100  * ( FARPROC * ) ( & this->umat ) = GetProcAddress( ( HMODULE ) this->umatobj, "umat_" ); //works for MinGW 32bit
101  if ( !this->umat ) {
102  // char *dlresult = GetLastError();
103  DWORD dlresult = GetLastError(); //works for MinGW 32bit
104  OOFEM_ERROR("Couldn't load symbol umat,\nerror code: %d\n", dlresult);
105  }
106 
107 #else
108  this->umatobj = dlopen(filename.c_str(), RTLD_NOW);
109  if ( !this->umatobj ) {
110  OOFEM_ERROR("couldn't load \"%s\",\ndlerror: %s", filename.c_str(), dlerror() );
111  }
112 
113  * ( void ** ) ( & this->umat ) = dlsym(this->umatobj, "umat_");
114 
115  char *dlresult = dlerror();
116  if ( dlresult ) {
117  OOFEM_ERROR("couldn't load symbol umat,\ndlerror: %s\n", dlresult);
118  }
119 
120 #endif
121 
123  mUseNumericalTangent = true;
124  }
125 
128  printf("mPerturbation: %e\n", mPerturbation);
129  }
130 
131  return IRRT_OK;
132 }
133 
135 {
137 
141  input.setField(std::string(this->cmname), _IFT_AbaqusUserMaterial_name);
142 }
143 
145 {
146  return new AbaqusUserMaterialStatus(n++, this->giveDomain(), gp, this->numState);
147 }
148 
149 
151  MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
152 {
153  AbaqusUserMaterialStatus *ms = dynamic_cast< AbaqusUserMaterialStatus * >( this->giveStatus(gp) );
154  if ( !ms->hasTangent() ) {
155  // Evaluating the function once, so that the tangent can be obtained.
156  FloatArray stress(6), strain(6);
157  strain.zero();
158  this->giveRealStressVector_3d(stress, gp, strain, tStep);
159  }
160 
161  answer = ms->giveTempTangent();
162 
163 #if 0
164  double h = 1e-7;
165  FloatArray strain, strainh, stress, stressh;
166  strain = static_cast< StructuralMaterialStatus * >( gp->giveMaterialStatus() )->giveTempStrainVector();
167  stress = static_cast< StructuralMaterialStatus * >( gp->giveMaterialStatus() )->giveTempStressVector();
168  FloatMatrix En( strain.giveSize(), strain.giveSize() );
169  for ( int i = 1; i <= strain.giveSize(); ++i ) {
170  strainh = strain;
171  strainh.at(i) += h;
172  this->giveRealStressVector_3d(stressh, form, gp, strainh, tStep);
173  stressh.subtract(stress);
174  stressh.times(1.0 / h);
175  En.setColumn(stressh, i);
176  }
177  this->giveRealStressVector_3d(stress, form, gp, strain, tStep);
178 
179  printf("En = ");
180  En.printYourself();
181  printf("Tangent = ");
182  answer.printYourself();
183 #endif
184 }
185 
186 
188  MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
189 {
190  AbaqusUserMaterialStatus *ms = dynamic_cast< AbaqusUserMaterialStatus * >( this->giveStatus(gp) );
191  if ( !ms->hasTangent() ) {
192  // Evaluating the function once, so that the tangent can be obtained.
193  FloatArray stress, vF;
194  vF = ms->giveTempFVector();
195  this->giveFirstPKStressVector_3d(stress, gp, vF, tStep);
196  }
197 
198  if ( !mUseNumericalTangent ) {
199  // if(mStressInterpretation == 0) {
200  answer = ms->giveTempTangent();
201  /*
202  * }
203  * else {
204  * // The Abaqus Documentation of User Subroutines for UMAT Section 1.1.31 says that DDSDDE is defined as
205  * // partial(Delta(sigma))/partial(Delta(epsilon)).
206  * FloatMatrix dSdE;
207  * dSdE = ms->giveTempTangent();
208  * this->give_dPdF_from(dSdE, answer, gp);
209  * }
210  */
211  } else {
212  double h = mPerturbation;
213  FloatArray vF, vF_h, stress, stressh;
214  vF = ms->giveTempFVector();
215  stress = ( ( StructuralMaterialStatus * ) gp->giveMaterialStatus() )->giveTempPVector();
216  FloatMatrix En(9, 9);
217  for ( int i = 1; i <= 9; ++i ) {
218  vF_h = vF;
219  vF_h.at(i) += h;
220  this->giveFirstPKStressVector_3d(stressh, gp, vF_h, tStep);
221  stressh.subtract(stress);
222  stressh.times(1.0 / h);
223  En.setColumn(stressh, i);
224  }
225 
226  // Reset
227  this->giveFirstPKStressVector_3d(stressh, gp, vF, tStep);
228 
229  /*
230  * printf("En = ");
231  * En.printYourself();
232  *
233  * answer = ms->giveTempTangent();
234  * printf("Tangent = ");
235  * answer.printYourself();
236  *
237  * FloatMatrix diff = En;
238  * diff.subtract(answer);
239  * printf("diff: "); diff.printYourself();
240  */
241 
242  answer = En;
243  }
244 }
245 
247  MatResponseMode mmode, GaussPoint *gp,
248  TimeStep *tStep)
249 {
250  FloatMatrix dPdF3d;
251  this->give3dMaterialStiffnessMatrix_dPdF(dPdF3d, mmode, gp, tStep);
252  StructuralMaterial :: giveReducedMatrixForm(answer, dPdF3d, _PlaneStrain);
253 }
254 
256  const FloatArray &reducedStrain, TimeStep *tStep)
257 {
258  AbaqusUserMaterialStatus *ms = static_cast< AbaqusUserMaterialStatus * >( this->giveStatus(gp) );
259  /* User-defined material name, left justified. Some internal material models are given names starting with
260  * the “ABQ_” character string. To avoid conflict, you should not use “ABQ_” as the leading string for
261  * CMNAME. */
262  //char cmname[80];
263 
264  // Sizes of the tensors.
265  int ndi = 3;
266  int nshr = 3;
267 
268  int ntens = ndi + nshr;
269  FloatArray strain = ms->giveStrainVector();
270  FloatArray stress = ms->giveStressVector();
271  // adding the initial stress
272  stress.add(initialStress);
273  FloatArray strainIncrement;
274  strainIncrement.beDifferenceOf(reducedStrain, strain);
275  FloatArray state = ms->giveStateVector();
276  FloatMatrix jacobian(ntens, ntens);
277  int numProperties = this->properties.giveSize();
278 
279  // Temperature and increment
280  double temp = 0.0, dtemp = 0.0;
281 
282  // Times and increment
283  double dtime = tStep->giveTimeIncrement();
285  double time [ 2 ] = {
286  tStep->giveTargetTime() - dtime, tStep->giveTargetTime()
287  };
288  double pnewdt = 1.0;
289 
290  /* Specific elastic strain energy, plastic dissipation, and “creep” dissipation, respectively. These are passed
291  * in as the values at the start of the increment and should be updated to the corresponding specific energy
292  * values at the end of the increment. They have no effect on the solution, except that they are used for
293  * energy output. */
294  double sse, spd, scd;
295 
296  // Outputs only in a fully coupled thermal-stress analysis:
297  double rpl = 0.0; // Volumetric heat generation per unit time at the end of the increment caused by mechanical working of the material.
298  FloatArray ddsddt(ntens); // Variation of the stress increments with respect to the temperature.
299  FloatArray drplde(ntens); // Variation of RPL with respect to the strain increments.
300  double drpldt = 0.0; // Variation of RPL with respect to the temperature.
301 
302  /* An array containing the coordinates of this point. These are the current coordinates if geometric
303  * nonlinearity is accounted for during the step (see “Procedures: overview,” Section 6.1.1); otherwise,
304  * the array contains the original coordinates of the point */
305  FloatArray coords;
307 
308  /* Rotation increment matrix. This matrix represents the increment of rigid body rotation of the basis
309  * system in which the components of stress (STRESS) and strain (STRAN) are stored. It is provided so
310  * that vector- or tensor-valued state variables can be rotated appropriately in this subroutine: stress and
311  * strain components are already rotated by this amount before UMAT is called. This matrix is passed in
312  * as a unit matrix for small-displacement analysis and for large-displacement analysis if the basis system
313  * for the material point rotates with the material (as in a shell element or when a local orientation is used).*/
314  FloatMatrix drot(3, 3);
315  drot.beUnitMatrix();
316 
317  /* Characteristic element length, which is a typical length of a line across an element for a first-order
318  * element; it is half of the same typical length for a second-order element. For beams and trusses it is a
319  * characteristic length along the element axis. For membranes and shells it is a characteristic length in
320  * the reference surface. For axisymmetric elements it is a characteristic length in the
321  * plane only.
322  * For cohesive elements it is equal to the constitutive thickness.*/
323  double celent = 0.0;
324 
325  /* Array containing the deformation gradient at the beginning of the increment. See the discussion
326  * regarding the availability of the deformation gradient for various element types. */
327  FloatMatrix dfgrd0(3, 3);
328  dfgrd0.beUnitMatrix();
329  /* Array containing the deformation gradient at the end of the increment. The components of this array
330  * are set to zero if nonlinear geometric effects are not included in the step definition associated with
331  * this increment. See the discussion regarding the availability of the deformation gradient for various
332  * element types. */
333  FloatMatrix dfgrd1(3, 3);
334  dfgrd1.beUnitMatrix();
335 
336  int noel = gp->giveElement()->giveNumber(); // Element number.
337  int npt = 0; // Integration point number.
338 
339  int layer = 0; // Layer number (for composite shells and layered solids)..
340  int kspt = 0; // Section point number within the current layer.
341  int kstep = 0; // Step number.
342  int kinc = 0; // Increment number.
343 
345  double predef;
346  double dpred;
347 
348  // Change to Abaqus's component order
349  stress.changeComponentOrder();
350  strain.changeComponentOrder();
351  strainIncrement.changeComponentOrder();
352 
353  // printf("stress oofem: "); stress.printYourself();
354 
355  OOFEM_LOG_DEBUG("AbaqusUserMaterial :: giveRealStressVector_3d - Calling subroutine");
356  this->umat(stress.givePointer(), // STRESS
357  state.givePointer(), // STATEV
358  jacobian.givePointer(), // DDSDDE
359  & sse, // SSE
360  & spd, // SPD
361  & scd, // SCD
362  & rpl, // RPL
363  ddsddt.givePointer(), // DDSDDT
364  drplde.givePointer(), // DRPLDE
365  & drpldt, // DRPLDT
366  strain.givePointer(), // STRAN
367  strainIncrement.givePointer(), // DSTRAN
368  time, // TIME
369  & dtime, // DTIME
370  & temp, // TEMP
371  & dtemp, // DTEMP
372  & predef, // PREDEF
373  & dpred, // DPRED
374  this->cmname, // CMNAME
375  & ndi, // NDI
376  & nshr, // NSHR
377  & ntens, // NTENS
378  & numState, // NSTATV
379  properties.givePointer(), // PROPS
380  & numProperties, // NPROPS
381  coords.givePointer(), // COORDS
382  drot.givePointer(), // DROT
383  & pnewdt, // PNEWDT
384  & celent, // CELENT
385  dfgrd0.givePointer(), // DFGRD0
386  dfgrd1.givePointer(), // DFGRD1
387  & noel, // NOEL
388  & npt, // NPT
389  & layer, // LAYER
390  & kspt, // KSPT
391  & kstep, // KSTEP
392  & kinc); // KINC
393 
394  // Change to OOFEM's component order
395  stress.changeComponentOrder();
396  // subtracking the initial stress
397  stress.subtract(initialStress);
398  strain.changeComponentOrder();
399  strainIncrement.changeComponentOrder();
400  jacobian.changeComponentOrder();
401 
402  ms->letTempStrainVectorBe(reducedStrain);
403  ms->letTempStressVectorBe(stress);
404  ms->letTempStateVectorBe(state);
405  ms->letTempTangentBe(jacobian);
406  answer = stress;
407 
408  OOFEM_LOG_DEBUG("AbaqusUserMaterial :: giveRealStressVector - Calling subroutine was successful");
409 }
410 
411 
413  const FloatArray &vF, TimeStep *tStep)
414 {
415  AbaqusUserMaterialStatus *ms = static_cast< AbaqusUserMaterialStatus * >( this->giveStatus(gp) );
416  /* User-defined material name, left justified. Some internal material models are given names starting with
417  * the “ABQ_” character string. To avoid conflict, you should not use “ABQ_” as the leading string for
418  * CMNAME. */
419  //char cmname[80];
420 
421  // Sizes of the tensors.
422  int ndi = 3;
423  int nshr = 3;
424 
425  int ntens = ndi + nshr;
426  FloatArray stress(9); // PK1
427  FloatArray strainIncrement;
428 
429  // compute Green-Lagrange strain
430  FloatArray strain;
431  FloatArray vS;
432  FloatMatrix F, E;
433  F.beMatrixForm(vF);
434  E.beTProductOf(F, F);
435  E.at(1, 1) -= 1.0;
436  E.at(2, 2) -= 1.0;
437  E.at(3, 3) -= 1.0;
438  E.times(0.5);
439  strain.beSymVectorFormOfStrain(E);
440 
441  strainIncrement.beDifferenceOf(strain, ms->giveStrainVector());
442  FloatArray state = ms->giveStateVector();
443  FloatMatrix jacobian(9, 9); // dPdF
444  int numProperties = this->properties.giveSize();
445 
446  // Temperature and increment
447  double temp = 0.0, dtemp = 0.0;
448 
449  // Times and increment
450  double dtime = tStep->giveTimeIncrement();
452  double time [ 2 ] = {
453  tStep->giveTargetTime() - dtime, tStep->giveTargetTime()
454  };
455  double pnewdt = 1.0;
456 
457  /* Specific elastic strain energy, plastic dissipation, and “creep” dissipation, respectively. These are passed
458  * in as the values at the start of the increment and should be updated to the corresponding specific energy
459  * values at the end of the increment. They have no effect on the solution, except that they are used for
460  * energy output. */
461  double sse, spd, scd;
462 
463  // Outputs only in a fully coupled thermal-stress analysis:
464  double rpl = 0.0; // Volumetric heat generation per unit time at the end of the increment caused by mechanical working of the material.
465  FloatArray ddsddt(ntens); // Variation of the stress increments with respect to the temperature.
466  FloatArray drplde(ntens); // Variation of RPL with respect to the strain increments.
467  double drpldt = 0.0; // Variation of RPL with respect to the temperature.
468 
469  /* An array containing the coordinates of this point. These are the current coordinates if geometric
470  * nonlinearity is accounted for during the step (see “Procedures: overview,” Section 6.1.1); otherwise,
471  * the array contains the original coordinates of the point */
472  FloatArray coords;
474 
475  /* Rotation increment matrix. This matrix represents the increment of rigid body rotation of the basis
476  * system in which the components of stress (STRESS) and strain (STRAN) are stored. It is provided so
477  * that vector- or tensor-valued state variables can be rotated appropriately in this subroutine: stress and
478  * strain components are already rotated by this amount before UMAT is called. This matrix is passed in
479  * as a unit matrix for small-displacement analysis and for large-displacement analysis if the basis system
480  * for the material point rotates with the material (as in a shell element or when a local orientation is used).*/
481  FloatMatrix drot(3, 3);
482  drot.beUnitMatrix();
483 
484  /* Characteristic element length, which is a typical length of a line across an element for a first-order
485  * element; it is half of the same typical length for a second-order element. For beams and trusses it is a
486  * characteristic length along the element axis. For membranes and shells it is a characteristic length in
487  * the reference surface. For axisymmetric elements it is a characteristic length in the
488  * plane only.
489  * For cohesive elements it is equal to the constitutive thickness.*/
490  double celent = 0.0;
491 
492  /* Array containing the deformation gradient at the beginning of the increment. See the discussion
493  * regarding the availability of the deformation gradient for various element types. */
494  FloatMatrix dfgrd0(3, 3);
495  /* Array containing the deformation gradient at the end of the increment. The components of this array
496  * are set to zero if nonlinear geometric effects are not included in the step definition associated with
497  * this increment. See the discussion regarding the availability of the deformation gradient for various
498  * element types. */
499  FloatMatrix dfgrd1(3, 3);
500 
501  dfgrd0.beMatrixForm( ms->giveFVector() );
502  dfgrd1.beMatrixForm(vF);
503 
504  int noel = gp->giveElement()->giveNumber(); // Element number.
505  int npt = 0; // Integration point number.
506 
507  // We intentionally ignore the layer number since that is handled by the layered cross-section in OOFEM.
508  int layer = 0; // Layer number (for composite shells and layered solids)..
509  int kspt = 0; // Section point number within the current layer.
510  int kstep = tStep->giveMetaStepNumber(); // Step number.
511  int kinc = 0; // Increment number.
512 
514  double predef;
515  double dpred;
516 
517  // Change to Abaqus's component order
518  stress.changeComponentOrder();
519  strain.changeComponentOrder();
520  strainIncrement.changeComponentOrder();
521 
522  OOFEM_LOG_DEBUG("AbaqusUserMaterial :: giveRealStressVector - Calling subroutine");
523  this->umat(stress.givePointer(), // STRESS
524  state.givePointer(), // STATEV
525  jacobian.givePointer(), // DDSDDE
526  & sse, // SSE
527  & spd, // SPD
528  & scd, // SCD
529  & rpl, // RPL
530  ddsddt.givePointer(), // DDSDDT
531  drplde.givePointer(), // DRPLDE
532  & drpldt, // DRPLDT
533  strain.givePointer(), // STRAN
534  strainIncrement.givePointer(), // DSTRAN
535  time, // TIME
536  & dtime, // DTIME
537  & temp, // TEMP
538  & dtemp, // DTEMP
539  & predef, // PREDEF
540  & dpred, // DPRED
541  this->cmname, // CMNAME
542  & ndi, // NDI
543  & nshr, // NSHR
544  & ntens, // NTENS
545  & numState, // NSTATV
546  properties.givePointer(), // PROPS
547  & numProperties, // NPROPS
548  coords.givePointer(), // COORDS
549  drot.givePointer(), // DROT
550  & pnewdt, // PNEWDT
551  & celent, // CELENT
552  dfgrd0.givePointer(), // DFGRD0
553  dfgrd1.givePointer(), // DFGRD1
554  & noel, // NOEL
555  & npt, // NPT
556  & layer, // LAYER
557  & kspt, // KSPT
558  & kstep, // KSTEP
559  & kinc); // KINC
560 
561 
562  // Change to OOFEM's component order
563  stress.changeComponentOrder();
564  strain.changeComponentOrder();
565  strainIncrement.changeComponentOrder();
566  jacobian.changeComponentOrder();
567 
568 
569  if ( mStressInterpretation == 0 ) {
570  FloatMatrix P, cauchyStress;
571  P.beMatrixForm(stress);
572 
573  FloatArray vP;
574  vP.beVectorForm(P);
575 
576  cauchyStress.beProductTOf(P, F);
577  cauchyStress.times( 1.0 / F.giveDeterminant() );
578 
579  FloatArray vCauchyStress;
580  vCauchyStress.beSymVectorForm(cauchyStress);
581 
582  ms->letTempStrainVectorBe(strain);
583  ms->letTempStressVectorBe(vCauchyStress);
584  ms->letTempStateVectorBe(state);
585  ms->letTempTangentBe(jacobian);
586  ms->letTempPVectorBe(vP);
587  ms->letTempFVectorBe(vF);
588 
589  answer = vP;
590  } else {
591  FloatArray vP;
592  FloatMatrix P, sigma, F_inv;
593  F_inv.beInverseOf(F);
594 
595  sigma.beMatrixForm(stress);
596  P.beProductOf(F, sigma);
597  vP.beVectorForm(P);
598  answer = vP;
599 
600  // Convert from sigma to S
601  FloatMatrix S;
602  S.beProductOf(F_inv, P);
603  vS.beSymVectorForm(S);
604 
605  // @todo Should convert from dsigmadE to dSdE here
606  // L2=detF*matmul( matmul ( inv(op_a_V9(F,F), cm-op_a_V9(ident,Tau)-op_b_V9(Tau,ident) ), inv(op_a_V9(Ft,Ft)))
607 
608  ms->letTempStrainVectorBe(strain);
609  ms->letTempStressVectorBe(vS);
610  ms->letTempStateVectorBe(state);
611  ms->letTempTangentBe(jacobian);
612  ms->letTempPVectorBe(vP);
613  ms->letTempFVectorBe(vF);
614  }
615 
616  OOFEM_LOG_DEBUG("AbaqusUserMaterial :: giveRealStressVector_3d - Calling subroutine was successful");
617 }
618 
619 
621 {
622  AbaqusUserMaterialStatus *ms = static_cast< AbaqusUserMaterialStatus * >( this->giveStatus(gp) );
623  if ( type == IST_Undefined || type == IST_AbaqusStateVector ) {
624  // The undefined value is used to just dump the entire state vector.
625  answer = ms->giveStateVector();
626  return 1;
627  } else if ( type == IST_StressTensor ) {
628  answer = ms->giveStressVector();
629  answer.add(initialStress);
630  return 1;
631  } else {
632  return StructuralMaterial :: giveIPValue(answer, gp, type, tStep);
633  }
634 }
635 
636 
637 
639 {
641  tempStateVector = stateVector;
642 }
643 
645  StructuralMaterialStatus(n, d, gp),
646  numState(numState), stateVector(numState), tempStateVector(numState), hasTangentFlag(false)
647 {
648  strainVector.resize(6);
649  strainVector.zero();
650 }
651 
653 {
656 }
657 
659 // Prints the strains and stresses on the data file.
660 {
662 
663  fprintf(File, " stateVector ");
664  FloatArray state = this->giveStateVector();
665  for ( auto &var : state ) {
666  fprintf( File, " % .4e", var );
667  }
668 
669  fprintf(File, "\n");
670 }
671 
672 
673 
674 } // end namespace oofem
double giveDeterminant() const
Returns the trace (sum of diagonal components) of the receiver.
Definition: floatmatrix.C:1408
virtual void updateYourself(TimeStep *tStep)
Update equilibrium history variables according to temp-variables.
Definition: structuralms.C:96
InternalStateType
Type representing the physical meaning of element or constitutive model internal variable.
void setField(int item, InputFieldType id)
void subtract(const FloatArray &src)
Subtracts array src to receiver.
Definition: floatarray.C:258
#define _IFT_AbaqusUserMaterial_initialStress
void letTempStrainVectorBe(const FloatArray &v)
Assigns tempStrainVector to given vector v.
Definition: structuralms.h:137
virtual MaterialStatus * giveStatus(GaussPoint *gp) const
Returns material status of receiver in given integration point.
Definition: material.C:244
virtual void giveRealStressVector_3d(FloatArray &answer, GaussPoint *gp, const FloatArray &reducedStrain, TimeStep *tStep)
Default implementation relies on giveRealStressVector for second Piola-Kirchoff stress.
Class and object Domain.
Definition: domain.h:115
void letTempFVectorBe(const FloatArray &v)
Assigns tempFVector to given vector v.
Definition: structuralms.h:143
double mPerturbation
Size of perturbation if numerical tangent is used.
FloatArray & letTempStateVectorBe(FloatArray &s)
AbaqusUserMaterial(int n, Domain *d)
Constructor.
FloatArray tempStateVector
Temporary state vector.
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in Reduced form.
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
int numState
Size of the state vector.
This class implements a structural material status information.
Definition: structuralms.h:65
char cmname[80]
Name for material routine.
virtual void printOutputAt(FILE *file, TimeStep *tStep)
Print receiver&#39;s output to given stream.
virtual void give3dMaterialStiffnessMatrix_dPdF(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
virtual bool hasField(InputFieldType id)=0
Returns true if record contains field identified by idString keyword.
void beVectorForm(const FloatMatrix &aMatrix)
Reciever will be a vector with 9 components formed from a 3x3 matrix.
Definition: floatarray.C:992
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in Reduced form.
double giveTargetTime()
Returns target time.
Definition: timestep.h:146
virtual void giveInputRecord(DynamicInputRecord &input)
Setups the input record string of receiver.
Element * giveElement()
Returns corresponding element to receiver.
Definition: gausspoint.h:188
#define S(p)
Definition: mdm.C:481
AbaqusUserMaterialStatus(int n, Domain *d, GaussPoint *gp, int numState)
Constructor.
const double * givePointer() const
Exposes the internal values of the matrix.
Definition: floatmatrix.h:558
#define _IFT_AbaqusUserMaterial_numericalTangentPerturbation
virtual void printOutputAt(FILE *file, TimeStep *tStep)
Print receiver&#39;s output to given stream.
Definition: structuralms.C:73
#define OOFEM_LOG_DEBUG(...)
Definition: logger.h:128
MatResponseMode
Describes the character of characteristic material matrix.
virtual void give3dMaterialStiffnessMatrix(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
Computes full 3d material stiffness matrix at given integration point, time, respecting load history ...
#define _IFT_AbaqusUserMaterial_numState
void changeComponentOrder()
Swaps the fourth and sixth index in the array.
Definition: floatarray.C:1055
void beDifferenceOf(const FloatArray &a, const FloatArray &b)
Sets receiver to be a - b.
Definition: floatarray.C:341
double giveTimeIncrement()
Returns solution step associated time increment.
Definition: timestep.h:150
static int n
Gausspoint counter.
void beMatrixForm(const FloatArray &aArray)
Definition: floatmatrix.C:1657
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
FloatArray stateVector
General state vector.
bool mUseNumericalTangent
Flag to determine if numerical tangent should be used.
#define E(p)
Definition: mdm.C:368
#define _IFT_AbaqusUserMaterial_numericalTangent
#define OOFEM_ERROR(...)
Definition: error.h:61
const FloatArray & giveStateVector() const
void times(double f)
Multiplies receiver by factor f.
Definition: floatmatrix.C:1594
virtual void giveFirstPKStressVector_3d(FloatArray &answer, GaussPoint *gp, const FloatArray &reducedF, TimeStep *tStep)
Default implementation relies on giveRealStressVector for second Piola-Kirchoff stress.
std::string filename
Name of the file that contains the umat function.
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
#define _IFT_AbaqusUserMaterial_name
virtual ~AbaqusUserMaterial()
Destructor.
virtual IRResultType initializeFrom(InputRecord *ir)
Reads the following values;.
Abstract base class representing a material status information.
Definition: matstatus.h:84
virtual void initTempStatus()
Initializes the temporary internal variables, describing the current state according to previously re...
Class representing vector of real numbers.
Definition: floatarray.h:82
int mStressInterpretation
Flag to determine how the stress and Jacobian are interpreted.
FloatArray initialStress
Initial stress.
FloatArray strainVector
Equilibrated strain vector in reduced form.
Definition: structuralms.h:69
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
FloatArray properties
Material properties.
void changeComponentOrder()
Swaps the indices in the 6x6 matrix such that it converts between OOFEM&#39;s and Abaqus&#39; way of writing ...
Definition: floatmatrix.C:1691
void letTempStressVectorBe(const FloatArray &v)
Assigns tempStressVector to given vector v.
Definition: structuralms.h:135
const FloatArray & giveStressVector() const
Returns the const pointer to receiver&#39;s stress vector.
Definition: structuralms.h:107
IntegrationPointStatus * giveMaterialStatus()
Returns reference to associated material status (NULL if not defined).
Definition: gausspoint.h:205
void * umatobj
Dynamically loaded umat.
const FloatArray & giveFVector() const
Returns the const pointer to receiver&#39;s deformation gradient vector.
Definition: structuralms.h:113
int giveMetaStepNumber()
Returns receiver&#39;s meta step number.
Definition: timestep.h:135
Class representing the general Input Record.
Definition: inputrecord.h:101
static void giveReducedMatrixForm(FloatMatrix &answer, const FloatMatrix &full, MaterialMode matMode)
Converts the full symmetric Voigt matrix (4th order tensor) to reduced form.
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
void beTProductOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:367
const FloatArray & giveTempFVector() const
Returns the const pointer to receiver&#39;s temporary deformation gradient vector.
Definition: structuralms.h:123
void setColumn(const FloatArray &src, int c)
Sets the values of the matrix in specified column.
Definition: floatmatrix.C:648
Class representing the a dynamic Input Record.
void times(double s)
Multiplies receiver with scalar.
Definition: floatarray.C:818
const double * givePointer() const
Gives the pointer to the raw data, breaking encapsulation.
Definition: floatarray.h:469
virtual void updateYourself(TimeStep *tStep)
Update equilibrium history variables according to temp-variables.
void printYourself() const
Prints matrix to stdout. Useful for debugging.
Definition: floatmatrix.C:1458
void beSymVectorForm(const FloatMatrix &aMatrix)
Reciever will be a vector with 6 components formed from a 3x3 matrix.
Definition: floatarray.C:1035
Abstract base class for all "structural" constitutive models.
void beProductTOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:397
virtual void initTempStatus()
Initializes the temporary internal variables, describing the current state according to previously re...
Definition: structuralms.C:108
Domain * giveDomain() const
Definition: femcmpnn.h:100
void beUnitMatrix()
Sets receiver to unity matrix.
Definition: floatmatrix.C:1332
virtual void giveInputRecord(DynamicInputRecord &input)
Setups the input record string of receiver.
REGISTER_Material(DummyMaterial)
void beSymVectorFormOfStrain(const FloatMatrix &aMatrix)
Definition: floatarray.C:1014
#define IR_GIVE_OPTIONAL_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:78
void beProductOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:337
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.
#define _IFT_AbaqusUserMaterial_properties
void(* umat)(double *stress, double *statev, double *ddsdde, double *sse, double *spd, double *scd, double *rpl, double *ddsddt, double *drplde, double *drpldt, double *stran, double *dstran, double time[2], double *dtime, double *temp, double *dtemp, double predef[1], double dpred[1], char cmname[80], int *ndi, int *nshr, int *ntens, int *nstatv, double *props, int *nprops, double coords[3], double *drot, double *pnewdt, double *celent, double *dfgrd0, double *dfgrd1, int *noel, int *npt, int *layer, int *kspt, int *kstep, int *kinc)
Pointer to the dynamically loaded umat-function (translated to C)
#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
void beInverseOf(const FloatMatrix &src)
Modifies receiver to become inverse of given parameter.
Definition: floatmatrix.C:835
const FloatArray & giveStrainVector() const
Returns the const pointer to receiver&#39;s strain vector.
Definition: structuralms.h:105
virtual MaterialStatus * CreateStatus(GaussPoint *gp) const
Creates new copy of associated status and inserts it into given integration point.
Class representing integration point in finite element program.
Definition: gausspoint.h:93
#define _IFT_AbaqusUserMaterial_userMaterial
Class representing solution step.
Definition: timestep.h:80
virtual void givePlaneStrainStiffMtrx_dPdF(FloatMatrix &answer, MatResponseMode mmode, GaussPoint *gp, TimeStep *tStep)
virtual int computeGlobalCoordinates(FloatArray &answer, const FloatArray &lcoords)
Computes the global coordinates from given element&#39;s local coordinates.
Definition: element.C:1207
void add(const FloatArray &src)
Adds array src to receiver.
Definition: floatarray.C:156
const FloatArray & giveNaturalCoordinates()
Returns coordinate array of receiver.
Definition: gausspoint.h:138
void letTempPVectorBe(const FloatArray &v)
Assigns tempPVector to given vector v.
Definition: structuralms.h:139
void resize(int s)
Resizes receiver towards requested size.
Definition: floatarray.C:631
const FloatMatrix & giveTempTangent()

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