www.mooseframework.org
AbaqusUmatMaterial.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
8 #include "AbaqusUmatMaterial.h"
9 #include "Factory.h"
10 #include "MooseMesh.h"
11 
12 #include <dlfcn.h>
13 #define QUOTE(macro) stringifyName(macro)
14 
15 template <>
16 InputParameters
18 {
19  InputParameters params = validParams<SolidModel>();
20  params.addRequiredParam<FileName>(
21  "plugin", "The path to the compiled dynamic library for the plugin you want to use");
22  params.addRequiredParam<std::vector<Real>>("mechanical_constants",
23  "Mechanical Material Properties");
24  params.addParam<std::vector<Real>>("thermal_constants", "Thermal Material Properties");
25  params.addRequiredParam<unsigned int>("num_state_vars",
26  "The number of state variables this UMAT is going to use");
27  return params;
28 }
29 
30 AbaqusUmatMaterial::AbaqusUmatMaterial(const InputParameters & parameters)
31  : SolidModel(parameters),
32  _plugin(getParam<FileName>("plugin")),
33  _mechanical_constants(getParam<std::vector<Real>>("mechanical_constants")),
34  _thermal_constants(getParam<std::vector<Real>>("thermal_constants")),
35  _num_state_vars(getParam<unsigned int>("num_state_vars")),
36  _grad_disp_x(coupledGradient("disp_x")),
37  _grad_disp_y(coupledGradient("disp_y")),
38  _grad_disp_z(coupledGradient("disp_z")),
39  _grad_disp_x_old(coupledGradientOld("disp_x")),
40  _grad_disp_y_old(coupledGradientOld("disp_y")),
41  _grad_disp_z_old(coupledGradientOld("disp_z")),
42  _state_var(declareProperty<std::vector<Real>>("state_var")),
43  _state_var_old(getMaterialPropertyOld<std::vector<Real>>("state_var")),
44  _Fbar(declareProperty<ColumnMajorMatrix>("Fbar")),
45  _Fbar_old(getMaterialPropertyOld<ColumnMajorMatrix>("Fbar")),
46  _elastic_strain_energy(declareProperty<Real>("elastic_strain_energy")),
47  _plastic_dissipation(declareProperty<Real>("plastic_dissipation")),
48  _creep_dissipation(declareProperty<Real>("creep_dissipation"))
49 {
50 #if defined(METHOD)
51  _plugin += std::string("-") + QUOTE(METHOD) + ".plugin";
52 #endif
53 
54  // Size and create full (mechanical+thermal) material property array
56  std::vector<Real> props_array(_num_props);
57  for (unsigned int i = 0; i < _mechanical_constants.size(); ++i)
58  props_array[i] = _mechanical_constants[i];
59  for (unsigned int i = _mechanical_constants.size(); i < _num_props; ++i)
60  props_array[i] = _thermal_constants[i];
61 
62  // Read mesh dimension and size UMAT arrays
63  if (_mesh.dimension() == 3) // 3D case
64  {
65  _NTENS = 6; // Size of the stress or strain component array (NDI+NSHR)
66  _NSHR = 3; // Number of engineering shear stress components
67  _NDI = 3; // Number of direct stress components (always 3)
68  }
69  else if (_mesh.dimension() == 2) // 2D case
70  {
71  _NTENS = 4;
72  _NSHR = 1;
73  _NDI = 3;
74  }
75 
76  _STATEV = new Real[_num_state_vars];
77  _DDSDDT = new Real[_NTENS];
78  _DRPLDE = new Real[_NTENS];
79  _STRAN = new Real[_NTENS];
80  _DFGRD0 = new Real[9];
81  _DFGRD1 = new Real[9];
82  _STRESS = new Real[_NTENS];
83  _DDSDDE = new Real[_NTENS * _NTENS];
84  _DSTRAN = new Real[_NTENS];
85  _PROPS = new Real[_num_props];
86 
87  for (unsigned int i = 0; i < _num_state_vars; ++i)
88  {
89  _STATEV[i] = 0.0;
90  }
91 
92  for (int i = 0; i < _NTENS; ++i)
93  {
94  _DDSDDT[i] = 0.0;
95  _DRPLDE[i] = 0.0;
96  _STRAN[i] = 0.0;
97  _STRESS[i] = 0.0;
98  _DSTRAN[i] = 0.0;
99  }
100 
101  for (unsigned int i = 0; i < 9; ++i)
102  {
103  _DFGRD0[i] = 0.0;
104  _DFGRD1[i] = 0.0;
105  }
106 
107  for (int i = 0; i < _NTENS * _NTENS; ++i)
108  {
109  _DDSDDE[i] = 0.0;
110  }
111 
112  // Assign materials properties from vector form into an array
113  for (unsigned int i = 0; i < _num_props; ++i)
114  {
115  _PROPS[i] = props_array[i];
116  }
117 
118  // Size UMAT state variable (NSTATV) and material constant (NPROPS) arrays
121 
122  // Open the library
123  _handle = dlopen(_plugin.c_str(), RTLD_LAZY);
124 
125  if (!_handle)
126  {
127  std::ostringstream error;
128  error << "Cannot open library: " << dlerror() << '\n';
129  mooseError(error.str());
130  }
131 
132  // Reset errors
133  dlerror();
134 
135  // Snag the function pointer from the library
136  {
137  void * pointer = dlsym(_handle, "umat_");
138  _umat = *reinterpret_cast<umat_t *>(&pointer);
139  }
140 
141  // Catch errors
142  const char * dlsym_error = dlerror();
143  if (dlsym_error)
144  {
145  dlclose(_handle);
146  std::ostringstream error;
147  error << "Cannot load symbol 'umat_': " << dlsym_error << '\n';
148  mooseError(error.str());
149  }
150 }
151 
153 {
154  delete _STATEV;
155  delete _DDSDDT;
156  delete _DRPLDE;
157  delete _STRAN;
158  delete _DFGRD0;
159  delete _DFGRD1;
160  delete _STRESS;
161  delete _DDSDDE;
162  delete _DSTRAN;
163  delete _PROPS;
164 
165  dlclose(_handle);
166 }
167 
168 void
170 {
171  // Initialize state variable vector
172  _state_var[_qp].resize(_num_state_vars);
173  for (unsigned int i = 0; i < _num_state_vars; ++i)
174  _state_var[_qp][i] = 0.0;
175 }
176 
177 void
179 {
180  // Calculate deformation gradient - modeled from "solid_mechanics/src/materials/Nonlinear3D.C"
181  // Fbar = 1 + grad(u(k))
182  ColumnMajorMatrix Fbar;
183 
184  Fbar(0, 0) = _grad_disp_x[_qp](0);
185  Fbar(0, 1) = _grad_disp_x[_qp](1);
186  Fbar(0, 2) = _grad_disp_x[_qp](2);
187  Fbar(1, 0) = _grad_disp_y[_qp](0);
188  Fbar(1, 1) = _grad_disp_y[_qp](1);
189  Fbar(1, 2) = _grad_disp_y[_qp](2);
190  Fbar(2, 0) = _grad_disp_z[_qp](0);
191  Fbar(2, 1) = _grad_disp_z[_qp](1);
192  Fbar(2, 2) = _grad_disp_z[_qp](2);
193 
194  Fbar.addDiag(1);
195  _Fbar[_qp] = Fbar;
196 
197  Real myDFGRD0[9] = {_Fbar_old[_qp](0, 0),
198  _Fbar_old[_qp](1, 0),
199  _Fbar_old[_qp](2, 0),
200  _Fbar_old[_qp](0, 1),
201  _Fbar_old[_qp](1, 1),
202  _Fbar_old[_qp](2, 1),
203  _Fbar_old[_qp](0, 2),
204  _Fbar_old[_qp](1, 2),
205  _Fbar_old[_qp](2, 2)};
206  Real myDFGRD1[9] = {_Fbar[_qp](0, 0),
207  _Fbar[_qp](1, 0),
208  _Fbar[_qp](2, 0),
209  _Fbar[_qp](0, 1),
210  _Fbar[_qp](1, 1),
211  _Fbar[_qp](2, 1),
212  _Fbar[_qp](0, 2),
213  _Fbar[_qp](1, 2),
214  _Fbar[_qp](2, 2)};
215 
216  for (unsigned int i = 0; i < 9; ++i)
217  {
218  _DFGRD0[i] = myDFGRD0[i];
219  _DFGRD1[i] = myDFGRD1[i];
220  }
221 
222  // Recover "old" state variables
223  for (unsigned int i = 0; i < _num_state_vars; ++i)
224  _STATEV[i] = _state_var_old[_qp][i];
225 
226  // Pass through updated stress, total strain, and strain increment arrays
227  for (int i = 0; i < _NTENS; ++i)
228  {
229  _STRESS[i] = _stress_old.component(i);
230  _STRAN[i] = _total_strain[_qp].component(i);
232  }
233 
234  // Pass through step , time, and coordinate system information
235  _KSTEP = _t_step; // Step number
236  _TIME[0] = _t; // Value of step time at the beginning of the current increment - Check
237  _TIME[1] = _t - _dt; // Value of total time at the beginning of the current increment - Check
238  _DTIME = _dt; // Time increment
239  for (unsigned int i = 0; i < 3; ++i) // Loop current coordinates in UMAT COORDS
240  _COORDS[i] = _coord[i];
241 
242  // Connection to extern statement
243  _umat(_STRESS,
244  _STATEV,
245  _DDSDDE,
246  &_SSE,
247  &_SPD,
248  &_SCD,
249  &_RPL,
250  _DDSDDT,
251  _DRPLDE,
252  &_DRPLDT,
253  _STRAN,
254  _DSTRAN,
255  _TIME,
256  &_DTIME,
257  &_TEMP,
258  &_DTEMP,
259  _PREDEF,
260  _DPRED,
261  &_CMNAME,
262  &_NDI,
263  &_NSHR,
264  &_NTENS,
265  &_NSTATV,
266  _PROPS,
267  &_NPROPS,
268  _COORDS,
269  _DROT,
270  &_PNEWDT,
271  &_CELENT,
272  _DFGRD0,
273  _DFGRD1,
274  &_NOEL,
275  &_NPT,
276  &_LAYER,
277  &_KSPT,
278  &_KSTEP,
279  &_KINC);
280 
281  // Energy outputs
283  _plastic_dissipation[_qp] = _SPD;
284  _creep_dissipation[_qp] = _SCD;
285 
286  // Update state variables
287  for (unsigned int i = 0; i < _num_state_vars; ++i)
288  _state_var[_qp][i] = _STATEV[i];
289 
290  // Get new stress tensor - UMAT should update stress
291  SymmTensor stressnew(_STRESS[0], _STRESS[1], _STRESS[2], _STRESS[3], _STRESS[4], _STRESS[5]);
292  _stress[_qp] = stressnew;
293 }
const VariableGradient & _grad_disp_z
const VariableGradient & _grad_disp_y
InputParameters validParams< AbaqusUmatMaterial >()
const VariableGradient & _grad_disp_x
MaterialProperty< SymmTensor > & _total_strain
Definition: SolidModel.h:114
MaterialProperty< std::vector< Real > > & _state_var
SolidModel is the base class for all this module&#39;s solid mechanics material models.
Definition: SolidModel.h:31
const MaterialProperty< std::vector< Real > > & _state_var_old
virtual void computeStress()
Compute the stress (sigma += deltaSigma)
MaterialProperty< Real > & _elastic_strain_energy
std::vector< Real > _mechanical_constants
void(* umat_t)(Real STRESS[], Real STATEV[], Real DDSDDE[], Real *SSE, Real *SPD, Real *SCD, Real *RPL, Real DDSDDT[], Real DRPLDE[], Real *DRPLDT, Real STRAN[], Real DSTRAN[], Real TIME[], Real *DTIME, Real *TEMP, Real *DTEMP, Real PREDEF[], Real DPRED[], Real *CMNAME, int *NDI, int *NSHR, int *NTENS, int *NSTATV, Real PROPS[], int *NPROPS, Real COORDS[], Real DROT[][3], Real *PNEWDT, Real *CELENT, Real DFGRD0[], Real DFGRD1[], int *NOEL, int *NPT, int *LAYER, int *KSPT, int *KSTEP, int *KINC)
MaterialProperty< Real > & _plastic_dissipation
InputParameters validParams< SolidModel >()
Definition: SolidModel.C:27
std::vector< Real > _thermal_constants
SymmTensor _strain_increment
Definition: SolidModel.h:143
SymmTensor _stress_old
Definition: SolidModel.h:112
virtual void initQpStatefulProperties()
MaterialProperty< Real > & _creep_dissipation
AbaqusUmatMaterial(const InputParameters &parameters)
MaterialProperty< SymmTensor > & _stress
Definition: SolidModel.h:106
unsigned int _num_state_vars
MaterialProperty< ColumnMajorMatrix > & _Fbar
const MaterialProperty< ColumnMajorMatrix > & _Fbar_old
Real component(unsigned int i) const
Definition: SymmTensor.h:97