www.mooseframework.org
MultiPlasticityLinearSystem.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 /****************************************************************/
8 
9 // Following is for perturbing distances in eliminating linearly-dependent directions
10 #include "MooseRandom.h"
11 
12 // Following is used to access PETSc's LAPACK routines
13 #include <petscblaslapack.h>
14 
15 template <>
16 InputParameters
18 {
19  InputParameters params = validParams<MultiPlasticityRawComponentAssembler>();
20  params.addRangeCheckedParam<Real>("linear_dependent",
21  1E-4,
22  "linear_dependent>=0 & linear_dependent<1",
23  "Flow directions are considered linearly dependent if the "
24  "smallest singular value is less than linear_dependent times "
25  "the largest singular value");
26  return params;
27 }
28 
31  _svd_tol(_params.get<Real>("linear_dependent")),
32  _min_f_tol(-1.0)
33 {
34  for (unsigned model = 0; model < _num_models; ++model)
35  if (_min_f_tol == -1.0 || _min_f_tol > _f[model]->_f_tol)
36  _min_f_tol = _f[model]->_f_tol;
37 
38  MooseRandom::seed(0);
39 }
40 
41 int
42 MultiPlasticityLinearSystem::singularValuesOfR(const std::vector<RankTwoTensor> & r,
43  std::vector<Real> & s)
44 {
45  int bm = r.size();
46  int bn = 6;
47 
48  s.resize(std::min(bm, bn));
49 
50  // prepare for gesvd or gesdd routine provided by PETSc
51  // Want to find the singular values of matrix
52  // ( r[0](0,0) r[0](0,1) r[0](0,2) r[0](1,1) r[0](1,2) r[0](2,2) )
53  // ( r[1](0,0) r[1](0,1) r[1](0,2) r[1](1,1) r[1](1,2) r[1](2,2) )
54  // a = ( r[2](0,0) r[2](0,1) r[2](0,2) r[2](1,1) r[2](1,2) r[2](2,2) )
55  // ( r[3](0,0) r[3](0,1) r[3](0,2) r[3](1,1) r[3](1,2) r[3](2,2) )
56  // ( r[4](0,0) r[4](0,1) r[4](0,2) r[4](1,1) r[4](1,2) r[4](2,2) )
57  // bm = 5
58 
59  std::vector<double> a(bm * 6);
60  // Fill in the a "matrix" by going down columns
61  unsigned ind = 0;
62  for (int col = 0; col < 3; ++col)
63  for (int row = 0; row < bm; ++row)
64  a[ind++] = r[row](0, col);
65  for (int col = 3; col < 5; ++col)
66  for (int row = 0; row < bm; ++row)
67  a[ind++] = r[row](1, col - 2);
68  for (int row = 0; row < bm; ++row)
69  a[ind++] = r[row](2, 2);
70 
71  // u and vt are dummy variables because they won't
72  // get referenced due to the "N" and "N" choices
73  int sizeu = 1;
74  std::vector<double> u(sizeu);
75  int sizevt = 1;
76  std::vector<double> vt(sizevt);
77 
78  int sizework = 16 * (bm + 6); // this is above the lowerbound specified in the LAPACK doco
79  std::vector<double> work(sizework);
80 
81  int info;
82 
83  LAPACKgesvd_("N",
84  "N",
85  &bm,
86  &bn,
87  &a[0],
88  &bm,
89  &s[0],
90  &u[0],
91  &sizeu,
92  &vt[0],
93  &sizevt,
94  &work[0],
95  &sizework,
96  &info);
97 
98  return info;
99 }
100 
101 void
103  const std::vector<Real> & intnl,
104  const std::vector<Real> & f,
105  const std::vector<RankTwoTensor> & r,
106  const std::vector<bool> & active,
107  std::vector<bool> & deactivated_due_to_ld)
108 {
109  deactivated_due_to_ld.resize(_num_surfaces, false);
110 
111  unsigned num_active = r.size();
112 
113  if (num_active <= 1)
114  return;
115 
116  std::vector<double> s;
117  int info = singularValuesOfR(r, s);
118  if (info != 0)
119  mooseError("In finding the SVD in the return-map algorithm, the PETSC LAPACK gesvd routine "
120  "returned with error code ",
121  info);
122 
123  // num_lin_dep are the number of linearly dependent
124  // "r vectors", if num_active <= 6
125  unsigned int num_lin_dep = 0;
126 
127  unsigned i = s.size();
128  while (i-- > 0)
129  if (s[i] < _svd_tol * s[0])
130  num_lin_dep++;
131  else
132  break;
133 
134  if (num_lin_dep == 0 && num_active <= 6)
135  return;
136 
137  // From here on, some flow directions are linearly dependent
138 
139  // Find the signed "distance" of the current (stress, internal) configuration
140  // from the yield surfaces. This distance will not be precise, but
141  // i want to preferentially deactivate yield surfaces that are close
142  // to the current stress point.
143  std::vector<RankTwoTensor> df_dstress;
144  dyieldFunction_dstress(stress, intnl, active, df_dstress);
145 
146  typedef std::pair<Real, unsigned> pair_for_sorting;
147  std::vector<pair_for_sorting> dist(num_active);
148  for (unsigned i = 0; i < num_active; ++i)
149  {
150  dist[i].first = f[i] / df_dstress[i].L2norm();
151  dist[i].second = i;
152  }
153  std::sort(dist.begin(), dist.end()); // sorted in ascending order
154 
155  // There is a potential problem when we have equal f[i], for it can give oscillations
156  bool equals_detected = false;
157  for (unsigned i = 0; i < num_active - 1; ++i)
158  if (std::abs(dist[i].first - dist[i + 1].first) < _min_f_tol)
159  {
160  equals_detected = true;
161  dist[i].first += _min_f_tol * (MooseRandom::rand() - 0.5);
162  }
163  if (equals_detected)
164  std::sort(dist.begin(), dist.end()); // sorted in ascending order
165 
166  std::vector<bool> scheduled_for_deactivation;
167  scheduled_for_deactivation.assign(num_active, false);
168 
169  // In the following loop we go through all the flow directions, from
170  // the one with the largest dist, to the one with the smallest dist,
171  // adding them one-by-one into r_tmp. Upon each addition we check
172  // for linear-dependence. if LD is found, we schedule the most
173  // recently added flow direction for deactivation, and pop it
174  // back off r_tmp
175  unsigned current_yf;
176  current_yf = dist[num_active - 1].second;
177  // the one with largest dist
178  std::vector<RankTwoTensor> r_tmp = {r[current_yf]};
179 
180  unsigned num_kept_active = 1;
181  for (unsigned yf_to_try = 2; yf_to_try <= num_active; ++yf_to_try)
182  {
183  current_yf = dist[num_active - yf_to_try].second;
184  if (num_active == 2) // shortcut to we don't have to singularValuesOfR
185  scheduled_for_deactivation[current_yf] = true;
186  else if (num_kept_active >= 6) // shortcut to we don't have to singularValuesOfR: there can
187  // never be > 6 linearly-independent r vectors
188  scheduled_for_deactivation[current_yf] = true;
189  else
190  {
191  r_tmp.push_back(r[current_yf]);
192  info = singularValuesOfR(r_tmp, s);
193  if (info != 0)
194  mooseError("In finding the SVD in the return-map algorithm, the PETSC LAPACK gesvd routine "
195  "returned with error code ",
196  info);
197  if (s[s.size() - 1] < _svd_tol * s[0])
198  {
199  scheduled_for_deactivation[current_yf] = true;
200  r_tmp.pop_back();
201  num_lin_dep--;
202  }
203  else
204  num_kept_active++;
205  if (num_lin_dep == 0 && num_active <= 6)
206  // have taken out all the vectors that were linearly dependent
207  // so no point continuing
208  break;
209  }
210  }
211 
212  unsigned int old_active_number = 0;
213  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
214  if (active[surface])
215  {
216  if (scheduled_for_deactivation[old_active_number])
217  deactivated_due_to_ld[surface] = true;
218  old_active_number++;
219  }
220 }
221 
222 void
224  const std::vector<Real> & intnl_old,
225  const std::vector<Real> & intnl,
226  const std::vector<Real> & pm,
227  const RankTwoTensor & delta_dp,
228  std::vector<Real> & f,
229  std::vector<RankTwoTensor> & r,
230  RankTwoTensor & epp,
231  std::vector<Real> & ic,
232  const std::vector<bool> & active)
233 {
234  // see comments at the start of .h file
235 
236  mooseAssert(intnl_old.size() == _num_models,
237  "Size of intnl_old is " << intnl_old.size()
238  << " which is incorrect in calculateConstraints");
239  mooseAssert(intnl.size() == _num_models,
240  "Size of intnl is " << intnl.size() << " which is incorrect in calculateConstraints");
241  mooseAssert(pm.size() == _num_surfaces,
242  "Size of pm is " << pm.size() << " which is incorrect in calculateConstraints");
243  mooseAssert(active.size() == _num_surfaces,
244  "Size of active is " << active.size()
245  << " which is incorrect in calculateConstraints");
246 
247  // yield functions
248  yieldFunction(stress, intnl, active, f);
249 
250  // flow directions and "epp"
251  flowPotential(stress, intnl, active, r);
252  epp = RankTwoTensor();
253  unsigned ind = 0;
254  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
255  if (active[surface])
256  epp += pm[surface] * r[ind++]; // note, even the deactivated_due_to_ld must get added in
257  epp -= delta_dp;
258 
259  // internal constraints
260  std::vector<Real> h;
261  hardPotential(stress, intnl, active, h);
262  ic.resize(0);
263  ind = 0;
264  std::vector<unsigned int> active_surfaces;
265  std::vector<unsigned int>::iterator active_surface;
266  for (unsigned model = 0; model < _num_models; ++model)
267  {
268  activeSurfaces(model, active, active_surfaces);
269  if (active_surfaces.size() > 0)
270  {
271  // some surfaces are active in this model, so must form an internal constraint
272  ic.push_back(intnl[model] - intnl_old[model]);
273  for (active_surface = active_surfaces.begin(); active_surface != active_surfaces.end();
274  ++active_surface)
275  ic[ic.size() - 1] += pm[*active_surface] * h[ind++]; // we know the correct one is h[ind]
276  // since it was constructed in the same
277  // manner
278  }
279  }
280 }
281 
282 void
283 MultiPlasticityLinearSystem::calculateRHS(const RankTwoTensor & stress,
284  const std::vector<Real> & intnl_old,
285  const std::vector<Real> & intnl,
286  const std::vector<Real> & pm,
287  const RankTwoTensor & delta_dp,
288  std::vector<Real> & rhs,
289  const std::vector<bool> & active,
290  bool eliminate_ld,
291  std::vector<bool> & deactivated_due_to_ld)
292 {
293  // see comments at the start of .h file
294 
295  mooseAssert(intnl_old.size() == _num_models,
296  "Size of intnl_old is " << intnl_old.size() << " which is incorrect in calculateRHS");
297  mooseAssert(intnl.size() == _num_models,
298  "Size of intnl is " << intnl.size() << " which is incorrect in calculateRHS");
299  mooseAssert(pm.size() == _num_surfaces,
300  "Size of pm is " << pm.size() << " which is incorrect in calculateRHS");
301  mooseAssert(active.size() == _num_surfaces,
302  "Size of active is " << active.size() << " which is incorrect in calculateRHS");
303 
304  std::vector<Real> f; // the yield functions
305  RankTwoTensor epp; // the plastic-strain constraint ("direction constraint")
306  std::vector<Real> ic; // the "internal constraints"
307 
308  std::vector<RankTwoTensor> r;
309  calculateConstraints(stress, intnl_old, intnl, pm, delta_dp, f, r, epp, ic, active);
310 
311  if (eliminate_ld)
312  eliminateLinearDependence(stress, intnl, f, r, active, deactivated_due_to_ld);
313  else
314  deactivated_due_to_ld.assign(_num_surfaces, false);
315 
316  std::vector<bool> active_not_deact(_num_surfaces);
317  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
318  active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
319 
320  unsigned num_active_f = 0;
321  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
322  if (active_not_deact[surface])
323  num_active_f++;
324 
325  unsigned num_active_ic = 0;
326  for (unsigned model = 0; model < _num_models; ++model)
327  if (anyActiveSurfaces(model, active_not_deact))
328  num_active_ic++;
329 
330  unsigned int dim = 3;
331  unsigned int system_size = 6 + num_active_f + num_active_ic; // "6" comes from symmeterizing epp,
332  // num_active_f comes from "f",
333  // num_active_f comes from "ic"
334 
335  rhs.resize(system_size);
336 
337  unsigned ind = 0;
338  for (unsigned i = 0; i < dim; ++i)
339  for (unsigned j = 0; j <= i; ++j)
340  rhs[ind++] = -epp(i, j);
341  unsigned active_surface = 0;
342  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
343  if (active[surface])
344  {
345  if (!deactivated_due_to_ld[surface])
346  rhs[ind++] = -f[active_surface];
347  active_surface++;
348  }
349  unsigned active_model = 0;
350  for (unsigned model = 0; model < _num_models; ++model)
351  if (anyActiveSurfaces(model, active))
352  {
353  if (anyActiveSurfaces(model, active_not_deact))
354  rhs[ind++] = -ic[active_model];
355  active_model++;
356  }
357 
358  mooseAssert(ind == system_size, "Incorrect filling of the rhs in calculateRHS");
359 }
360 
361 void
363  const std::vector<Real> & intnl,
364  const std::vector<Real> & pm,
365  const RankFourTensor & E_inv,
366  const std::vector<bool> & active,
367  const std::vector<bool> & deactivated_due_to_ld,
368  std::vector<std::vector<Real>> & jac)
369 {
370  // see comments at the start of .h file
371 
372  mooseAssert(intnl.size() == _num_models,
373  "Size of intnl is " << intnl.size() << " which is incorrect in calculateJacobian");
374  mooseAssert(pm.size() == _num_surfaces,
375  "Size of pm is " << pm.size() << " which is incorrect in calculateJacobian");
376  mooseAssert(active.size() == _num_surfaces,
377  "Size of active is " << active.size() << " which is incorrect in calculateJacobian");
378  mooseAssert(deactivated_due_to_ld.size() == _num_surfaces,
379  "Size of deactivated_due_to_ld is " << deactivated_due_to_ld.size()
380  << " which is incorrect in calculateJacobian");
381 
382  unsigned ind = 0;
383  unsigned active_surface_ind = 0;
384 
385  std::vector<bool> active_surface(_num_surfaces); // active and not deactivated_due_to_ld
386  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
387  active_surface[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
388  unsigned num_active_surface = 0;
389  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
390  if (active_surface[surface])
391  num_active_surface++;
392 
393  std::vector<bool> active_model(
394  _num_models); // whether a model has surfaces that are active and not deactivated_due_to_ld
395  for (unsigned model = 0; model < _num_models; ++model)
396  active_model[model] = anyActiveSurfaces(model, active_surface);
397 
398  unsigned num_active_model = 0;
399  for (unsigned model = 0; model < _num_models; ++model)
400  if (active_model[model])
401  num_active_model++;
402 
403  ind = 0;
404  std::vector<unsigned int> active_model_index(_num_models);
405  for (unsigned model = 0; model < _num_models; ++model)
406  if (active_model[model])
407  active_model_index[model] = ind++;
408  else
409  active_model_index[model] =
410  _num_models + 1; // just a dummy, that will probably cause a crash if something goes wrong
411 
412  std::vector<RankTwoTensor> df_dstress;
413  dyieldFunction_dstress(stress, intnl, active_surface, df_dstress);
414 
415  std::vector<Real> df_dintnl;
416  dyieldFunction_dintnl(stress, intnl, active_surface, df_dintnl);
417 
418  std::vector<RankTwoTensor> r;
419  flowPotential(stress, intnl, active, r);
420 
421  std::vector<RankFourTensor> dr_dstress;
422  dflowPotential_dstress(stress, intnl, active, dr_dstress);
423 
424  std::vector<RankTwoTensor> dr_dintnl;
425  dflowPotential_dintnl(stress, intnl, active, dr_dintnl);
426 
427  std::vector<Real> h;
428  hardPotential(stress, intnl, active, h);
429 
430  std::vector<RankTwoTensor> dh_dstress;
431  dhardPotential_dstress(stress, intnl, active, dh_dstress);
432 
433  std::vector<Real> dh_dintnl;
434  dhardPotential_dintnl(stress, intnl, active, dh_dintnl);
435 
436  // d(epp)/dstress = sum_{active alpha} pm[alpha]*dr_dstress
437  RankFourTensor depp_dstress;
438  ind = 0;
439  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
440  if (active[surface]) // includes deactivated_due_to_ld
441  depp_dstress += pm[surface] * dr_dstress[ind++];
442  depp_dstress += E_inv;
443 
444  // d(epp)/dpm_{active_surface_index} = r_{active_surface_index}
445  std::vector<RankTwoTensor> depp_dpm;
446  depp_dpm.resize(num_active_surface);
447  ind = 0;
448  active_surface_ind = 0;
449  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
450  {
451  if (active[surface])
452  {
453  if (active_surface[surface]) // do not include the deactived_due_to_ld, since their pm are not
454  // dofs in the NR
455  depp_dpm[active_surface_ind++] = r[ind];
456  ind++;
457  }
458  }
459 
460  // d(epp)/dintnl_{active_model_index} = sum(pm[asdf]*dr_dintnl[fdsa])
461  std::vector<RankTwoTensor> depp_dintnl;
462  depp_dintnl.assign(num_active_model, RankTwoTensor());
463  ind = 0;
464  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
465  {
466  if (active[surface])
467  {
468  unsigned int model_num = modelNumber(surface);
469  if (active_model[model_num]) // only include models with surfaces which are still active after
470  // deactivated_due_to_ld
471  depp_dintnl[active_model_index[model_num]] += pm[surface] * dr_dintnl[ind];
472  ind++;
473  }
474  }
475 
476  // df_dstress has been calculated above
477  // df_dpm is always zero
478  // df_dintnl has been calculated above, but only the active_surface+active_model stuff needs to be
479  // included in Jacobian: see below
480 
481  std::vector<RankTwoTensor> dic_dstress;
482  dic_dstress.assign(num_active_model, RankTwoTensor());
483  ind = 0;
484  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
485  {
486  if (active[surface])
487  {
488  unsigned int model_num = modelNumber(surface);
489  if (active_model[model_num]) // only include ic for models with active_surface (ie, if model
490  // only contains deactivated_due_to_ld don't include it)
491  dic_dstress[active_model_index[model_num]] += pm[surface] * dh_dstress[ind];
492  ind++;
493  }
494  }
495 
496  std::vector<std::vector<Real>> dic_dpm;
497  dic_dpm.resize(num_active_model);
498  ind = 0;
499  active_surface_ind = 0;
500  for (unsigned model = 0; model < num_active_model; ++model)
501  dic_dpm[model].assign(num_active_surface, 0);
502  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
503  {
504  if (active[surface])
505  {
506  if (active_surface[surface]) // only take derivs wrt active-but-not-deactivated_due_to_ld pm
507  {
508  unsigned int model_num = modelNumber(surface);
509  // if (active_model[model_num]) // do not need this check as if the surface has
510  // active_surface, the model must be deemed active!
511  dic_dpm[active_model_index[model_num]][active_surface_ind] = h[ind];
512  active_surface_ind++;
513  }
514  ind++;
515  }
516  }
517 
518  std::vector<std::vector<Real>> dic_dintnl;
519  dic_dintnl.resize(num_active_model);
520  for (unsigned model = 0; model < num_active_model; ++model)
521  {
522  dic_dintnl[model].assign(num_active_model, 0);
523  dic_dintnl[model][model] = 1; // deriv wrt internal parameter
524  }
525  ind = 0;
526  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
527  {
528  if (active[surface])
529  {
530  unsigned int model_num = modelNumber(surface);
531  if (active_model[model_num]) // only the models that contain surfaces that are still active
532  // after deactivation_due_to_ld
533  dic_dintnl[active_model_index[model_num]][active_model_index[model_num]] +=
534  pm[surface] * dh_dintnl[ind];
535  ind++;
536  }
537  }
538 
539  unsigned int dim = 3;
540  unsigned int system_size =
541  6 + num_active_surface + num_active_model; // "6" comes from symmeterizing epp
542  jac.resize(system_size);
543  for (unsigned i = 0; i < system_size; ++i)
544  jac[i].assign(system_size, 0);
545 
546  unsigned int row_num = 0;
547  unsigned int col_num = 0;
548  for (unsigned i = 0; i < dim; ++i)
549  for (unsigned j = 0; j <= i; ++j)
550  {
551  for (unsigned k = 0; k < dim; ++k)
552  for (unsigned l = 0; l <= k; ++l)
553  jac[col_num][row_num++] =
554  depp_dstress(i, j, k, l) +
555  (k != l ? depp_dstress(i, j, l, k)
556  : 0); // extra part is needed because i assume dstress(i, j) = dstress(j, i)
557  for (unsigned surface = 0; surface < num_active_surface; ++surface)
558  jac[col_num][row_num++] = depp_dpm[surface](i, j);
559  for (unsigned a = 0; a < num_active_model; ++a)
560  jac[col_num][row_num++] = depp_dintnl[a](i, j);
561  row_num = 0;
562  col_num++;
563  }
564 
565  ind = 0;
566  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
567  if (active_surface[surface])
568  {
569  for (unsigned k = 0; k < dim; ++k)
570  for (unsigned l = 0; l <= k; ++l)
571  jac[col_num][row_num++] =
572  df_dstress[ind](k, l) +
573  (k != l ? df_dstress[ind](l, k)
574  : 0); // extra part is needed because i assume dstress(i, j) = dstress(j, i)
575  for (unsigned beta = 0; beta < num_active_surface; ++beta)
576  jac[col_num][row_num++] = 0; // df_dpm
577  for (unsigned model = 0; model < _num_models; ++model)
578  if (active_model[model]) // only use df_dintnl for models in active_model
579  {
580  if (modelNumber(surface) == model)
581  jac[col_num][row_num++] = df_dintnl[ind];
582  else
583  jac[col_num][row_num++] = 0;
584  }
585  ind++;
586  row_num = 0;
587  col_num++;
588  }
589 
590  for (unsigned a = 0; a < num_active_model; ++a)
591  {
592  for (unsigned k = 0; k < dim; ++k)
593  for (unsigned l = 0; l <= k; ++l)
594  jac[col_num][row_num++] =
595  dic_dstress[a](k, l) +
596  (k != l ? dic_dstress[a](l, k)
597  : 0); // extra part is needed because i assume dstress(i, j) = dstress(j, i)
598  for (unsigned alpha = 0; alpha < num_active_surface; ++alpha)
599  jac[col_num][row_num++] = dic_dpm[a][alpha];
600  for (unsigned b = 0; b < num_active_model; ++b)
601  jac[col_num][row_num++] = dic_dintnl[a][b];
602  row_num = 0;
603  col_num++;
604  }
605 
606  mooseAssert(col_num == system_size, "Incorrect filling of cols in Jacobian");
607 }
608 
609 void
610 MultiPlasticityLinearSystem::nrStep(const RankTwoTensor & stress,
611  const std::vector<Real> & intnl_old,
612  const std::vector<Real> & intnl,
613  const std::vector<Real> & pm,
614  const RankFourTensor & E_inv,
615  const RankTwoTensor & delta_dp,
616  RankTwoTensor & dstress,
617  std::vector<Real> & dpm,
618  std::vector<Real> & dintnl,
619  const std::vector<bool> & active,
620  std::vector<bool> & deactivated_due_to_ld)
621 {
622  // Calculate RHS and Jacobian
623  std::vector<Real> rhs;
624  calculateRHS(stress, intnl_old, intnl, pm, delta_dp, rhs, active, true, deactivated_due_to_ld);
625 
626  std::vector<std::vector<Real>> jac;
627  calculateJacobian(stress, intnl, pm, E_inv, active, deactivated_due_to_ld, jac);
628 
629  // prepare for LAPACKgesv_ routine provided by PETSc
630  int system_size = rhs.size();
631 
632  std::vector<double> a(system_size * system_size);
633  // Fill in the a "matrix" by going down columns
634  unsigned ind = 0;
635  for (int col = 0; col < system_size; ++col)
636  for (int row = 0; row < system_size; ++row)
637  a[ind++] = jac[row][col];
638 
639  int nrhs = 1;
640  std::vector<int> ipiv(system_size);
641  int info;
642  LAPACKgesv_(&system_size, &nrhs, &a[0], &system_size, &ipiv[0], &rhs[0], &system_size, &info);
643 
644  if (info != 0)
645  mooseError("In solving the linear system in a Newton-Raphson process, the PETSC LAPACK gsev "
646  "routine returned with error code ",
647  info);
648 
649  // Extract the results back to dstress, dpm and dintnl
650  std::vector<bool> active_not_deact(_num_surfaces);
651  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
652  active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
653 
654  unsigned int dim = 3;
655  ind = 0;
656 
657  for (unsigned i = 0; i < dim; ++i)
658  for (unsigned j = 0; j <= i; ++j)
659  dstress(i, j) = dstress(j, i) = rhs[ind++];
660  dpm.assign(_num_surfaces, 0);
661  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
662  if (active_not_deact[surface])
663  dpm[surface] = rhs[ind++];
664  dintnl.assign(_num_models, 0);
665  for (unsigned model = 0; model < _num_models; ++model)
666  if (anyActiveSurfaces(model, active_not_deact))
667  dintnl[model] = rhs[ind++];
668 
669  mooseAssert(static_cast<int>(ind) == system_size,
670  "Incorrect extracting of changes from NR solution in nrStep");
671 }
virtual void calculateRHS(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, std::vector< Real > &rhs, const std::vector< bool > &active, bool eliminate_ld, std::vector< bool > &deactivated_due_to_ld)
Calculate the RHS which is rhs = -(epp(0,0), epp(1,0), epp(1,1), epp(2,0), epp(2,1), epp(2,2), f[0], f[1], ..., f[num_f], ic[0], ic[1], ..., ic[num_ic])
virtual int singularValuesOfR(const std::vector< RankTwoTensor > &r, std::vector< Real > &s)
Performs a singular-value decomposition of r and returns the singular values.
virtual void dyieldFunction_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &df_dstress)
The derivative of the active yield function(s) with respect to stress.
Real _svd_tol
Tolerance on the minimum ratio of singular values before flow-directions are deemed linearly dependen...
unsigned int modelNumber(unsigned int surface)
returns the model number, given the surface number
virtual void dhardPotential_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &dh_dintnl)
The derivative of the active hardening potentials with respect to the active internal parameters...
virtual void calculateJacobian(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankFourTensor &E_inv, const std::vector< bool > &active, const std::vector< bool > &deactivated_due_to_ld, std::vector< std::vector< Real >> &jac)
d(rhs)/d(dof)
virtual void eliminateLinearDependence(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< Real > &f, const std::vector< RankTwoTensor > &r, const std::vector< bool > &active, std::vector< bool > &deactivated_due_to_ld)
Performs a number of singular-value decompositions to check for linear-dependence of the active direc...
InputParameters validParams< MultiPlasticityRawComponentAssembler >()
virtual void yieldFunction(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &f)
The active yield function(s)
virtual void dhardPotential_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &dh_dstress)
The derivative of the active hardening potentials with respect to stress By assumption in the Userobj...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
MultiPlasticityRawComponentAssembler holds and computes yield functions, flow directions, etc, for use in FiniteStrainMultiPlasticity.
MultiPlasticityLinearSystem(const MooseObject *moose_object)
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
virtual void dflowPotential_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &dr_dintnl)
The derivative of the active flow potentials with respect to the active internal parameters The UserO...
bool anyActiveSurfaces(int model, const std::vector< bool > &active)
returns true if any internal surfaces of the given model are active according to &#39;active&#39; ...
virtual void flowPotential(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &r)
The active flow potential(s) - one for each yield function.
unsigned int _num_models
Number of plastic models for this material.
virtual void dflowPotential_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankFourTensor > &dr_dstress)
The derivative of the active flow potential(s) with respect to stress.
InputParameters validParams< MultiPlasticityLinearSystem >()
virtual void hardPotential(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &h)
The active hardening potentials (one for each internal parameter and for each yield function) by assu...
virtual void dyieldFunction_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &df_dintnl)
The derivative of active yield function(s) with respect to their internal parameters (the user object...
virtual void nrStep(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankFourTensor &E_inv, const RankTwoTensor &delta_dp, RankTwoTensor &dstress, std::vector< Real > &dpm, std::vector< Real > &dintnl, const std::vector< bool > &active, std::vector< bool > &deactivated_due_to_ld)
Performs one Newton-Raphson step.
virtual void calculateConstraints(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, std::vector< Real > &f, std::vector< RankTwoTensor > &r, RankTwoTensor &epp, std::vector< Real > &ic, const std::vector< bool > &active)
The constraints.
void activeSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces)
Returns the external surface number(s) of the active surfaces of the given model This may be of size=...
Real _min_f_tol
Minimum value of the _f_tol parameters for the Yield Function User Objects.