v0.14.0
ClosetPointProjection.cpp
Go to the documentation of this file.
1 /**
2  * \file ClosetPointProjection.cpp
3  * \ingroup adoc_plasticity
4  *
5  * \brief Implementation of nested integration algorithm for plasticity
6  *
7  */
8 
9 #include <MoFEM.hpp>
10 using namespace MoFEM;
11 
12 #include <ADOLCPlasticity.hpp>
13 namespace ADOLCPlasticity {
14 
15 ClosestPointProjection::ClosestPointProjection() {
16 
17  if (!LogManager::checkIfChannelExist("ADOLCPlasticityWold")) {
18  auto core_log = logging::core::get();
19 
20  core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(),
21  "ADOLCPlasticityWold"));
22  core_log->add_sink(LogManager::createSink(LogManager::getStrmSync(),
23  "ADOLCPlasticitySync"));
24  core_log->add_sink(LogManager::createSink(LogManager::getStrmSelf(),
25  "ADOLCPlasticitySelf"));
26 
27  LogManager::setLog("ADOLCPlasticityWold");
28  LogManager::setLog("ADOLCPlasticitySync");
29  LogManager::setLog("ADOLCPlasticitySelf");
30 
31  MOFEM_LOG_TAG("ADOLCPlasticityWold", "ADOL-C Plasticity");
32  MOFEM_LOG_TAG("ADOLCPlasticitySync", "ADOL-C Plasticity");
33  MOFEM_LOG_TAG("ADOLCPlasticitySelf", "ADOL-C Plasticity");
34  }
35 
36  tapesTags = {0, 1, 2}; //< set tapes default tags
37 
38 }
39 
40 MoFEMErrorCode ClosestPointProjection::recordW() {
41  MoFEMFunctionBegin;
42 
43  auto plastic_strain = getPlasticStrain();
44  auto total_strain = getTotalStrain();
45  auto internal_variables = getInternalVariables();
46 
47  trace_on(tapesTags[ClosestPointProjection::W]);
48  {
49  // active variables
50  a_plasticStrain.resize(6, false);
51  for (int ii = 0; ii < 6; ii++) {
52  a_plasticStrain[ii] <<= plastic_strain[ii];
53  }
54  a_sTrain.resize(6, false);
55  for (int ii = 0; ii < 6; ii++) {
56  a_sTrain[ii] <<= total_strain[ii];
57  }
58  a_internalVariables.resize(nbInternalVariables, false);
59  for (unsigned int ii = 0; ii < nbInternalVariables; ii++) {
60  a_internalVariables[ii] <<= internal_variables[ii];
61  }
62  // evaluete functions
63  CHKERR freeHemholtzEnergy();
64  a_w >>= w;
65  }
66  trace_off();
67  MoFEMFunctionReturn(0);
68 }
69 
70 MoFEMErrorCode ClosestPointProjection::recordY() {
71  MoFEMFunctionBegin;
72 
73  auto stress = getStress();
74  auto internal_fluxes = getInternalFluxes();
75 
76  trace_on(tapesTags[ClosestPointProjection::Y]);
77  {
78  // active variables
79  a_sTress.resize(6, false);
80  for (int ii = 0; ii < 6; ii++) {
81  a_sTress[ii] <<= stress[ii];
82  }
83  a_internalFluxes.resize(nbInternalVariables, false);
84  for (unsigned int ii = 0; ii < nbInternalVariables; ii++) {
85  a_internalFluxes[ii] <<= internal_fluxes[ii];
86  }
87  // evaluete functions
88  CHKERR yieldFunction();
89  // return variables
90  a_y >>= y;
91  }
92  trace_off();
93  MoFEMFunctionReturn(0);
94 }
95 
96 MoFEMErrorCode ClosestPointProjection::recordH() {
97  MoFEMFunctionBegin;
98 
99  auto stress = getStress();
100  auto internal_fluxes = getInternalFluxes();
101 
102  trace_on(tapesTags[ClosestPointProjection::H]);
103  {
104  // active variables
105  a_sTress.resize(6, false);
106  for (int ii = 0; ii < 6; ii++) {
107  a_sTress[ii] <<= stress[ii];
108  }
109  a_internalFluxes.resize(nbInternalVariables, false);
110  for (unsigned int ii = 0; ii < nbInternalVariables; ii++) {
111  a_internalFluxes[ii] <<= internal_fluxes[ii];
112  }
113  // evaluete functions
114  CHKERR flowPotential();
115  // return variables
116  a_h >>= h;
117  }
118  trace_off();
119  MoFEMFunctionReturn(0);
120 }
121 
122 MoFEMErrorCode ClosestPointProjection::playW() {
123  MoFEMFunctionBegin;
124  int r;
125  int adloc_return_value = 0;
126  r = ::function(tapesTags[ClosestPointProjection::W], 1,
127  activeVariablesW.size(), &activeVariablesW[0], &w);
128  if (r < adloc_return_value) {
129  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
130  "ADOL-C function evaluation with error");
131  }
132  r = gradient(tapesTags[ClosestPointProjection::W], activeVariablesW.size(),
133  &activeVariablesW[0], &gradientW[0]);
134  if (r < adloc_return_value) {
135  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
136  "ADOL-C function evaluation with error");
137  }
138  hessianW.resize(activeVariablesW.size(), activeVariablesW.size(), false);
139  hessianW.clear();
140  vector<double *> hessian_w(activeVariablesW.size());
141  for (unsigned int dd = 0; dd < activeVariablesW.size(); dd++) {
142  hessian_w[dd] = &hessianW(dd, 0);
143  }
144  r = hessian(tapesTags[ClosestPointProjection::W], activeVariablesW.size(),
145  &activeVariablesW[0], &hessian_w[0]);
146  if (r < adloc_return_value) {
147  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
148  "ADOL-C function evaluation with error");
149  }
150  partialWStrainPlasticStrain.resize(6, 6, false);
151  for (int ii = 0; ii < 6; ii++) {
152  for (int jj = 0; jj < 6; jj++) {
153  partialWStrainPlasticStrain(ii, jj) = hessianW(6 + ii, jj);
154  }
155  }
156  C.resize(6, 6, false);
157  for (int ii = 0; ii < 6; ii++) {
158  for (int jj = 0; jj <= ii; jj++) {
159  C(ii, jj) = hessianW(6 + ii, 6 + jj);
160  }
161  }
162  D.resize(nbInternalVariables, nbInternalVariables, false);
163  for (unsigned int ii = 0; ii < nbInternalVariables; ii++) {
164  for (unsigned int jj = 0; jj <= ii; jj++) {
165  D(ii, jj) = hessianW(12 + ii, 12 + jj);
166  }
167  }
168  MoFEMFunctionReturn(0);
169 }
170 
171 MoFEMErrorCode
172 ClosestPointProjection::playW_NoHessian() {
173  MoFEMFunctionBegin;
174  int r;
175  int adloc_return_value = 0;
176  r = ::function(tapesTags[ClosestPointProjection::W], 1,
177  activeVariablesW.size(), &activeVariablesW[0], &w);
178  if (r < adloc_return_value) {
179  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
180  "ADOL-C function evaluation with error");
181  }
182  r = gradient(tapesTags[ClosestPointProjection::W], activeVariablesW.size(),
183  &activeVariablesW[0], &gradientW[0]);
184  if (r < adloc_return_value) {
185  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
186  "ADOL-C function evaluation with error");
187  }
188  MoFEMFunctionReturn(0);
189 }
190 
191 MoFEMErrorCode ClosestPointProjection::playY() {
192  MoFEMFunctionBegin;
193  auto active_variables_yh = getActiveVariablesYH();
194  int r;
195  int adloc_return_value = 0;
196  r = ::function(tapesTags[ClosestPointProjection::Y], 1,
197  active_variables_yh.size(), &active_variables_yh[0], &y);
198  if (r < adloc_return_value) {
199  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
200  "ADOL-C function evaluation with error");
201  }
202  gradientY.resize(active_variables_yh.size());
203  r = gradient(tapesTags[ClosestPointProjection::Y], active_variables_yh.size(),
204  &active_variables_yh[0], &gradientY[0]);
205  if (r < adloc_return_value) {
206  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
207  "ADOL-C function evaluation with error");
208  }
209  partialYSigma =
210  VectorAdaptor(6, ublas::shallow_array_adaptor<double>(6, &gradientY[0]));
211  partialYFlux = VectorAdaptor(
212  nbInternalVariables,
213  ublas::shallow_array_adaptor<double>(nbInternalVariables, &gradientY[6]));
214  MoFEMFunctionReturn(0);
215 }
216 MoFEMErrorCode
217 ClosestPointProjection::playY_NoGradient() {
218  MoFEMFunctionBegin;
219  auto active_variables_yh = getActiveVariablesYH();
220  int r;
221  int adloc_return_value = 0;
222  r = ::function(tapesTags[ClosestPointProjection::Y], 1,
223  active_variables_yh.size(), &active_variables_yh[0], &y);
224  if (r < adloc_return_value) {
225  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
226  "ADOL-C function evaluation with error");
227  }
228  MoFEMFunctionReturn(0);
229 }
230 MoFEMErrorCode ClosestPointProjection::playH() {
231  MoFEMFunctionBegin;
232  auto active_variables_yh = getActiveVariablesYH();
233  int r;
234  int adloc_return_value = 0;
235  r = ::function(tapesTags[ClosestPointProjection::H], 1,
236  active_variables_yh.size(), &active_variables_yh[0], &h);
237  if (r < adloc_return_value) {
238  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
239  "ADOL-C function evaluation with error");
240  }
241  gradientH.resize(active_variables_yh.size());
242  r = gradient(tapesTags[ClosestPointProjection::H], active_variables_yh.size(),
243  &active_variables_yh[0], &gradientH[0]);
244  if (r < adloc_return_value) {
245  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
246  "ADOL-C function evaluation with error");
247  }
248  hessianH.resize(6 + nbInternalVariables, 6 + nbInternalVariables, false);
249  hessianH.clear();
250  vector<double *> hessian_h(active_variables_yh.size());
251  for (int dd = 0; dd < active_variables_yh.size(); dd++) {
252  hessian_h[dd] = &hessianH(dd, 0);
253  }
254  r = hessian(tapesTags[ClosestPointProjection::H], active_variables_yh.size(),
255  &active_variables_yh[0], &hessian_h[0]);
256  if (r < adloc_return_value) {
257  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
258  "ADOL-C function evaluation with error");
259  }
260  partialHSigma =
261  VectorAdaptor(6, ublas::shallow_array_adaptor<double>(6, &gradientH[0]));
262  partialHFlux = VectorAdaptor(
263  nbInternalVariables,
264  ublas::shallow_array_adaptor<double>(nbInternalVariables, &gradientH[6]));
265  partialHFlux *= -1;
266  partial2HSigma.resize(6, 6, false);
267  for (int ii = 0; ii < 6; ii++) {
268  for (int jj = 0; jj <= ii; jj++) {
269  partial2HSigma(ii, jj) = hessianH(ii, jj);
270  }
271  }
272  partial2HFlux.resize(nbInternalVariables, nbInternalVariables, false);
273  for (unsigned int ii = 0; ii < nbInternalVariables; ii++) {
274  for (unsigned int jj = 0; jj <= ii; jj++) {
275  partial2HFlux(ii, jj) = -hessianH(6 + ii, 6 + jj);
276  }
277  }
278  partial2HSigmaFlux.resize(6, nbInternalVariables, false);
279  partial2HSigmaFlux.clear();
280  for (int ii = 0; ii < 6; ii++) {
281  for (unsigned int jj = 0; jj < nbInternalVariables; jj++) {
282  partial2HSigmaFlux(ii, jj) = -hessianH(6 + jj, ii);
283  }
284  }
285  MoFEMFunctionReturn(0);
286 }
287 MoFEMErrorCode ClosestPointProjection::playH_NoHessian() {
288  MoFEMFunctionBegin;
289  auto active_variables_yh = getActiveVariablesYH();
290  int r;
291  int adloc_return_value = 0;
292  r = ::function(tapesTags[ClosestPointProjection::H], 1,
293  active_variables_yh.size(), &active_variables_yh[0], &h);
294  if (r < adloc_return_value) {
295  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
296  "ADOL-C function evaluation with error");
297  }
298  gradientH.resize(active_variables_yh.size());
299  r = gradient(tapesTags[ClosestPointProjection::H], active_variables_yh.size(),
300  &active_variables_yh[0], &gradientH[0]);
301  if (r < adloc_return_value) {
302  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
303  "ADOL-C function evaluation with error");
304  }
305  partialHSigma =
306  VectorAdaptor(6, ublas::shallow_array_adaptor<double>(6, &gradientH[0]));
307  partialHFlux = VectorAdaptor(
308  nbInternalVariables,
309  ublas::shallow_array_adaptor<double>(nbInternalVariables, &gradientH[6]));
310  partialHFlux *= -1;
311  MoFEMFunctionReturn(0);
312 }
313 
314 MoFEMErrorCode ClosestPointProjection::createMatAVecR() {
315  MoFEMFunctionBegin;
316  PetscInt n;
317  n = 1 + 6 + nbInternalVariables;
318  dataA.resize(n, n, false);
319  Mat A_tmp;
320  Vec R_tmp, Chi_tmp, dChi_tmp;
321  CHKERR MatCreateSeqDense(PETSC_COMM_SELF, n, n, &dataA(0, 0), &A_tmp);
322  CHKERR VecCreateSeq(PETSC_COMM_SELF, n, &R_tmp);
323  CHKERR VecCreateSeq(PETSC_COMM_SELF, n, &Chi_tmp);
324  CHKERR VecCreateSeq(PETSC_COMM_SELF, n, &dChi_tmp);
325 
326  A = SmartPetscObj<Mat>(A_tmp);
327  R = SmartPetscObj<Vec>(R_tmp);
328  Chi = SmartPetscObj<Vec>(Chi_tmp);
329  dChi = SmartPetscObj<Vec>(dChi_tmp);
330 
331  MoFEMFunctionReturn(0);
332 }
333 
334 MoFEMErrorCode
335 ClosestPointProjection::evaluatePotentials() {
336  MoFEMFunctionBegin;
337  CHKERR recordW();
338  CHKERR recordY();
339  CHKERR recordH();
340  MoFEMFunctionReturn(0);
341 }
342 
343 MoFEMErrorCode ClosestPointProjection::playPotentials() {
344  MoFEMFunctionBegin;
345  CHKERR playW();
346  CHKERR playY();
347  CHKERR playH();
348  MoFEMFunctionReturn(0);
349 }
350 
351 MoFEMErrorCode ClosestPointProjection::playPotentials_NoHessian() {
352  MoFEMFunctionBegin;
353  CHKERR playW_NoHessian();
354  CHKERR playY_NoGradient();
355  CHKERR playH_NoHessian();
356  MoFEMFunctionReturn(0);
357 }
358 
359 MoFEMErrorCode
360 ClosestPointProjection::calculateR(Vec R) {
361  MoFEMFunctionBegin;
362 
363  auto plastic_strain = getPlasticStrain();
364  auto internal_variables = getInternalVariables();
365 
366  // Residual
367  double *array;
368  CHKERR VecGetArray(R, &array);
369  for (int ii = 0; ii < 6; ii++) {
370  array[ii] = -plastic_strain[ii] + plasticStrain0[ii] +
371  deltaGamma * partialHSigma[ii];
372  }
373  for (unsigned int ii = 0; ii < nbInternalVariables; ii++) {
374  array[6 + ii] = -internal_variables[ii] + internalVariables0[ii] +
375  deltaGamma * partialHFlux[ii];
376  }
377  array[6 + nbInternalVariables] = y;
378  CHKERR VecRestoreArray(R, &array);
379  MoFEMFunctionReturn(0);
380 }
381 
382 MoFEMErrorCode ClosestPointProjection::calculateA() {
383  MoFEMFunctionBegin;
384  // matrix A
385  CHKERR MatZeroEntries(A);
386  // row 0 (Strain)
387  MatrixDouble a00 = prod(partial2HSigma, partialWStrainPlasticStrain);
388  MatrixDouble a01 = prod(partial2HSigmaFlux, D);
389  for (int ii = 0; ii < 6; ii++) {
390  for (int jj = 0; jj < 6; jj++) {
391  dataA(ii, jj) = deltaGamma * a00(ii, jj);
392  }
393  for (unsigned int jj = 0; jj < nbInternalVariables; jj++) {
394  dataA(ii, 6 + jj) = deltaGamma * a01(ii, jj);
395  }
396  dataA(ii, ii) -= 1;
397  }
398  for (int ii = 0; ii < 6; ii++) {
399  dataA(ii, 6 + nbInternalVariables) = partialHSigma[ii];
400  }
401  // row 1 (Fluxes)
402  MatrixDouble a11 = prod(partial2HFlux, D);
403  MatrixDouble a10 =
404  prod(trans(partial2HSigmaFlux), partialWStrainPlasticStrain);
405  for (unsigned int ii = 0; ii < nbInternalVariables; ii++) {
406  for (int jj = 0; jj < 6; jj++) {
407  dataA(6 + ii, jj) += deltaGamma * a10(ii, jj);
408  }
409  for (unsigned int jj = 0; jj < nbInternalVariables; jj++) {
410  dataA(6 + ii, 6 + jj) += deltaGamma * a11(ii, jj);
411  }
412  dataA(6 + ii, 6 + ii) -= 1;
413  }
414  for (unsigned int ii = 0; ii < nbInternalVariables; ii++) {
415  dataA(6 + ii, 6 + nbInternalVariables) = partialHFlux[ii];
416  }
417  // row 3 (Plastic multiplier)
418  VectorDouble partial_y_sigma_c =
419  prod(trans(partialWStrainPlasticStrain), partialYSigma);
420  VectorDouble partial_y_flux_d = prod(trans(D), partialYFlux);
421  for (unsigned int dd = 0; dd < partial_y_sigma_c.size(); dd++) {
422  dataA(6 + nbInternalVariables, dd) = partial_y_sigma_c[dd];
423  }
424  for (unsigned int dd = 0; dd < partial_y_flux_d.size(); dd++) {
425  dataA(6 + nbInternalVariables, 6 + dd) = partial_y_flux_d[dd];
426  }
427  dataA(6 + nbInternalVariables, 6 + nbInternalVariables) = 0;
428  CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
429  CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
430  MoFEMFunctionReturn(0);
431 }
432 
433 MoFEMErrorCode ClosestPointProjection::snesCreate() {
434  MoFEMFunctionBegin;
435  sNes = createSNES(PETSC_COMM_SELF);
436  CHKERR SNESSetFunction(sNes, R, ADOLCPlasticityRes, (void *)this);
437  CHKERR SNESSetJacobian(sNes, A, A, ADOLCPlasticityJac, (void *)this);
438 
439  KSP ksp;
440  PC pc;
441  SNESLineSearch line_search;
442  CHKERR SNESGetKSP(sNes, &ksp);
443  CHKERR KSPGetPC(ksp, &pc);
444  CHKERR KSPSetType(ksp, KSPPREONLY);
445  CHKERR PCSetType(pc, PCLU);
446  CHKERR SNESSetTolerances(sNes, 1e-8, 1e-6, 1e-12, 20, 1000);
447  CHKERR SNESGetLineSearch(sNes, &line_search);
448  // CHKERR SNESLineSearchSetType(line_search, SNESLINESEARCHBASIC);
449  CHKERR SNESLineSearchSetType(line_search, SNESLINESEARCHBT);
450  // CHKERR SNESLineSearchSetType(line_search, SNESLINESEARCHL2);
451 
452  CHKERR SNESAppendOptionsPrefix(sNes, "cp_");
453  CHKERR SNESSetFromOptions(sNes);
454 
455  MoFEMFunctionReturn(0);
456 }
457 
458 MoFEMErrorCode
459 ClosestPointProjection::solveClosetProjection() {
460  MoFEMFunctionBegin;
461 
462  auto plastic_strain = getPlasticStrain();
463  auto internal_variables = getInternalVariables();
464 
465  {
466  double *array;
467  CHKERR VecGetArray(Chi, &array);
468  for (int ii = 0; ii < 6; ii++) {
469  array[ii] = plastic_strain[ii];
470  }
471  int nb_internal_variables = nbInternalVariables;
472  for (unsigned int ii = 0; ii < nb_internal_variables; ii++) {
473  array[6 + ii] = internal_variables[ii];
474  }
475  deltaGamma = 0;
476  CHKERR VecRestoreArray(Chi, &array);
477  }
478 
479  CHKERR SNESSolve(sNes, PETSC_NULL, Chi);
480  SNESConvergedReason reason;
481  CHKERR SNESGetConvergedReason(sNes, &reason);
482  if (reason < 0) {
483  int its;
484  CHKERR SNESGetIterationNumber(sNes, &its);
485  MOFEM_LOG_C(
486  "ADOLCPlasticitySelf", Sev::warning,
487  "Plasticity Closet Point Projection - number of Newton iterations = %d",
488  its);
489  plastic_strain = plasticStrain0;
490  internal_variables = internalVariables0;
491  deltaGamma = 0;
492  } else {
493  double *array;
494  CHKERR VecGetArray(Chi, &array);
495  for (int ii = 0; ii < 6; ii++) {
496  plastic_strain[ii] = array[ii];
497  }
498  int nb_internal_variables = nbInternalVariables;
499  for (unsigned int ii = 0; ii < nb_internal_variables; ii++) {
500  internal_variables[ii] = array[6 + ii];
501  }
502  deltaGamma = array[6 + nb_internal_variables];
503  CHKERR VecRestoreArray(Chi, &array);
504  }
505  MoFEMFunctionReturn(0);
506 }
507 
508 MoFEMErrorCode
509 ClosestPointProjection::consistentTangent() {
510  MoFEMFunctionBegin;
511  CHKERR playPotentials();
512  CHKERR calculateA();
513  Ep.resize(6, 6, false);
514  MatrixDouble ep_row; // inrement internal varaibles as result of incremen of
515  // total strain
516  ep_row = deltaGamma * prod(partial2HSigma, C);
517  MatrixDouble alpha_row;
518  alpha_row = deltaGamma * prod(trans(partial2HSigmaFlux), C);
519  VectorDouble y_row;
520  y_row = prod(C, partialYSigma);
521  const int n = 6 + nbInternalVariables;
522  // Iterate strains
523  for (int dd = 0; dd < 6; dd++) {
524  CHKERR VecZeroEntries(R);
525  double *array;
526  CHKERR VecGetArray(R, &array);
527  {
528  for (auto ii = 0; ii < 6; ii++) {
529  array[ii] = ep_row(ii, dd);
530  }
531  for (auto ii = 0; ii < nbInternalVariables; ii++) {
532  array[6 + ii] = alpha_row(ii, dd);
533  }
534  array[n] = y_row[dd];
535  }
536  CHKERR VecRestoreArray(R, &array);
537  CHKERR VecAssemblyBegin(R);
538  CHKERR VecAssemblyEnd(R);
539 
540  KSP ksp;
541  CHKERR SNESGetKSP(sNes, &ksp);
542  CHKERR KSPSolve(ksp, R, dChi);
543  CHKERR VecGetArray(dChi, &array);
544  {
545  for (auto ii = 0; ii < 6; ii++) {
546  Ep(ii, dd) = array[ii];
547  }
548  }
549  CHKERR VecRestoreArray(dChi, &array);
550  }
551  Cp.resize(6, 6, false);
552  noalias(Cp) = prod(C, Ep);
553  Cep.resize(6, 6, false);
554  noalias(Cep) = C + Cp;
555  MoFEMFunctionReturn(0);
556 }
557 
558 MoFEMErrorCode ADOLCPlasticityRes(SNES snes, Vec chi, Vec r, void *ctx) {
559  MoFEMFunctionBegin;
560  ClosestPointProjection *cp;
561  cp = (ClosestPointProjection *)ctx;
562 
563  auto plastic_strain = cp->getPlasticStrain();
564  auto internal_variables = cp->getInternalVariables();
565 
566  {
567  const double *array;
568  CHKERR VecGetArrayRead(chi, &array);
569  for (auto ii = 0; ii < 6; ii++) {
570  plastic_strain[ii] = array[ii];
571  }
572  for (auto ii = 0; ii < cp->nbInternalVariables; ii++) {
573  internal_variables[ii] = array[6 + ii];
574  }
575  cp->deltaGamma = array[6 + cp->nbInternalVariables];
576  CHKERR VecRestoreArrayRead(chi, &array);
577  }
578  CHKERR cp->playPotentials_NoHessian();
579  CHKERR cp->calculateR(r);
580 
581  MoFEMFunctionReturn(0);
582 }
583 
584 MoFEMErrorCode ADOLCPlasticityJac(SNES snes, Vec chi, Mat A, Mat,
585  void *ctx) {
586  MoFEMFunctionBegin;
587  ClosestPointProjection *cp;
588  cp = (ClosestPointProjection *)ctx;
589 
590  auto plastic_strain = cp->getPlasticStrain();
591  auto internal_variables = cp->getInternalVariables();
592 
593  {
594  const double *array;
595  CHKERR VecGetArrayRead(chi, &array);
596  for (auto ii = 0; ii < 6; ii++) {
597  plastic_strain[ii] = array[ii];
598  }
599  for (auto ii = 0; ii < cp->nbInternalVariables; ii++) {
600  internal_variables[ii] = array[6 + ii];
601  }
602  cp->deltaGamma = array[6 + cp->nbInternalVariables];
603  CHKERR VecRestoreArrayRead(chi, &array);
604  }
605  CHKERR cp->playPotentials();
606  CHKERR cp->calculateA();
607  MoFEMFunctionReturn(0);
608 }
609 
610 MoFEMErrorCode ClosestPointProjection::recordTapes() {
611  MoFEMFunctionBegin;
612  VectorDouble active(12 + nbInternalVariables);
613  VectorDouble gradient(12 + nbInternalVariables);
614  auto tmp_active = getVectorAdaptor(&(active[0]), 12 + nbInternalVariables);
615  auto tmp_gradient =
616  getVectorAdaptor(&(gradient[0]), 12 + nbInternalVariables);
617  activeVariablesW.swap(tmp_active);
618  gradientW.swap(tmp_gradient);
619  CHKERR evaluatePotentials();
620  MoFEMFunctionReturn(0);
621 }
622 
623 } // namespace ADOLCPlasticity
ADOLCPlasticity::ClosestPointProjection::getInternalVariables
VectorAdaptor getInternalVariables()
Definition: ADOLCPlasticity.hpp:155
MoFEM::LogManager::checkIfChannelExist
static bool checkIfChannelExist(const std::string channel)
Check if channel exist.
Definition: LogManager.cpp:404
PlasticOps::w
double w(double eqiv, double dot_tau, double f, double sigma_y, double sigma_Y)
Definition: PlasticOpsGeneric.hpp:276
ADOLCPlasticity::ClosestPointProjection::deltaGamma
double deltaGamma
Definition: ADOLCPlasticity.hpp:180
MoFEM::createSNES
auto createSNES(MPI_Comm comm)
Definition: PetscSmartObj.hpp:251
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM.hpp
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
W
double W
Definition: free_surface.cpp:166
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
sdf.r
int r
Definition: sdf.py:8
ADOLCPlasticity::ClosestPointProjection::calculateR
MoFEMErrorCode calculateR(Vec R)
Definition: ClosetPointProjection.cpp:360
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
ADOLCPlasticity::ClosestPointProjection
Closest point projection algorithm.
Definition: ADOLCPlasticity.hpp:139
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
ADOLCPlasticity::ClosestPointProjection::getPlasticStrain
VectorAdaptor getPlasticStrain()
Definition: ADOLCPlasticity.hpp:149
ADOLCPlasticity::ClosestPointProjection::playPotentials
MoFEMErrorCode playPotentials()
Definition: ClosetPointProjection.cpp:343
R
@ R
Definition: free_surface.cpp:394
MoFEM::getVectorAdaptor
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:31
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
MoFEM::LogManager::getStrmSync
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
Definition: LogManager.cpp:348
h
double h
Definition: photon_diffusion.cpp:60
ADOLCPlasticity::ClosestPointProjection::playPotentials_NoHessian
MoFEMErrorCode playPotentials_NoHessian()
Definition: ClosetPointProjection.cpp:351
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
MoFEM::LogManager::getStrmSelf
static boost::shared_ptr< std::ostream > getStrmSelf()
Get the strm self object.
Definition: LogManager.cpp:340
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
ADOLCPlasticity
Definition: ADOLCPlasticity.hpp:24
MoFEM::Types::VectorAdaptor
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:115
ADOLCPlasticity::ClosestPointProjection::nbInternalVariables
int nbInternalVariables
Definition: ADOLCPlasticity.hpp:141
ADOLCPlasticity::ADOLCPlasticityRes
MoFEMErrorCode ADOLCPlasticityRes(SNES snes, Vec chi, Vec r, void *ctx)
Internal SNES function used at integration points to calulate stress.
Definition: ClosetPointProjection.cpp:558
convert.n
n
Definition: convert.py:82
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
ADOLCPlasticity.hpp
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
ADOLCPlasticity::ClosestPointProjection::calculateA
MoFEMErrorCode calculateA()
Definition: ClosetPointProjection.cpp:382
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
MoFEM::SmartPetscObj< Mat >
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
H
double H
Hardening.
Definition: plastic.cpp:175
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
ADOLCPlasticity::ADOLCPlasticityJac
MoFEMErrorCode ADOLCPlasticityJac(SNES, Vec, Mat, Mat, void *ctx)
Internal SNES function used at integration points to calulate tangent matrix.
Definition: ClosetPointProjection.cpp:584
Generated by Doxygen 1.8.17 and hosted at University of Glasgow