10 using namespace MoFEM;
15 ClosestPointProjection::ClosestPointProjection() {
18 auto core_log = logging::core::get();
21 "ADOLCPlasticityWold"));
23 "ADOLCPlasticitySync"));
25 "ADOLCPlasticitySelf"));
36 tapesTags = {0, 1, 2};
43 auto plastic_strain = getPlasticStrain();
44 auto total_strain = getTotalStrain();
45 auto internal_variables = getInternalVariables();
50 a_plasticStrain.resize(6,
false);
51 for (
int ii = 0; ii < 6; ii++) {
52 a_plasticStrain[ii] <<= plastic_strain[ii];
54 a_sTrain.resize(6,
false);
55 for (
int ii = 0; ii < 6; ii++) {
56 a_sTrain[ii] <<= total_strain[ii];
58 a_internalVariables.resize(nbInternalVariables,
false);
59 for (
unsigned int ii = 0; ii < nbInternalVariables; ii++) {
60 a_internalVariables[ii] <<= internal_variables[ii];
63 CHKERR freeHemholtzEnergy();
73 auto stress = getStress();
74 auto internal_fluxes = getInternalFluxes();
76 trace_on(tapesTags[ClosestPointProjection::Y]);
79 a_sTress.resize(6,
false);
80 for (
int ii = 0; ii < 6; ii++) {
81 a_sTress[ii] <<= stress[ii];
83 a_internalFluxes.resize(nbInternalVariables,
false);
84 for (
unsigned int ii = 0; ii < nbInternalVariables; ii++) {
85 a_internalFluxes[ii] <<= internal_fluxes[ii];
99 auto stress = getStress();
100 auto internal_fluxes = getInternalFluxes();
105 a_sTress.resize(6,
false);
106 for (
int ii = 0; ii < 6; ii++) {
107 a_sTress[ii] <<= stress[ii];
109 a_internalFluxes.resize(nbInternalVariables,
false);
110 for (
unsigned int ii = 0; ii < nbInternalVariables; ii++) {
111 a_internalFluxes[ii] <<= internal_fluxes[ii];
125 int adloc_return_value = 0;
127 activeVariablesW.size(), &activeVariablesW[0], &
w);
128 if (
r < adloc_return_value) {
130 "ADOL-C function evaluation with error");
133 &activeVariablesW[0], &gradientW[0]);
134 if (
r < adloc_return_value) {
136 "ADOL-C function evaluation with error");
138 hessianW.resize(activeVariablesW.size(), activeVariablesW.size(),
false);
140 vector<double *> hessian_w(activeVariablesW.size());
141 for (
unsigned int dd = 0;
dd < activeVariablesW.size();
dd++) {
142 hessian_w[
dd] = &hessianW(
dd, 0);
145 &activeVariablesW[0], &hessian_w[0]);
146 if (
r < adloc_return_value) {
148 "ADOL-C function evaluation with error");
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);
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);
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);
172 ClosestPointProjection::playW_NoHessian() {
175 int adloc_return_value = 0;
177 activeVariablesW.size(), &activeVariablesW[0], &
w);
178 if (
r < adloc_return_value) {
180 "ADOL-C function evaluation with error");
183 &activeVariablesW[0], &gradientW[0]);
184 if (
r < adloc_return_value) {
186 "ADOL-C function evaluation with error");
193 auto active_variables_yh = getActiveVariablesYH();
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) {
200 "ADOL-C function evaluation with error");
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) {
207 "ADOL-C function evaluation with error");
210 VectorAdaptor(6, ublas::shallow_array_adaptor<double>(6, &gradientY[0]));
213 ublas::shallow_array_adaptor<double>(nbInternalVariables, &gradientY[6]));
217 ClosestPointProjection::playY_NoGradient() {
219 auto active_variables_yh = getActiveVariablesYH();
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) {
226 "ADOL-C function evaluation with error");
232 auto active_variables_yh = getActiveVariablesYH();
234 int adloc_return_value = 0;
236 active_variables_yh.size(), &active_variables_yh[0], &
h);
237 if (
r < adloc_return_value) {
239 "ADOL-C function evaluation with error");
241 gradientH.resize(active_variables_yh.size());
243 &active_variables_yh[0], &gradientH[0]);
244 if (
r < adloc_return_value) {
246 "ADOL-C function evaluation with error");
248 hessianH.resize(6 + nbInternalVariables, 6 + nbInternalVariables,
false);
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);
255 &active_variables_yh[0], &hessian_h[0]);
256 if (
r < adloc_return_value) {
258 "ADOL-C function evaluation with error");
261 VectorAdaptor(6, ublas::shallow_array_adaptor<double>(6, &gradientH[0]));
264 ublas::shallow_array_adaptor<double>(nbInternalVariables, &gradientH[6]));
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);
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);
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);
289 auto active_variables_yh = getActiveVariablesYH();
291 int adloc_return_value = 0;
293 active_variables_yh.size(), &active_variables_yh[0], &
h);
294 if (
r < adloc_return_value) {
296 "ADOL-C function evaluation with error");
298 gradientH.resize(active_variables_yh.size());
300 &active_variables_yh[0], &gradientH[0]);
301 if (
r < adloc_return_value) {
303 "ADOL-C function evaluation with error");
306 VectorAdaptor(6, ublas::shallow_array_adaptor<double>(6, &gradientH[0]));
309 ublas::shallow_array_adaptor<double>(nbInternalVariables, &gradientH[6]));
317 n = 1 + 6 + nbInternalVariables;
318 dataA.resize(
n,
n,
false);
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);
335 ClosestPointProjection::evaluatePotentials() {
354 CHKERR playY_NoGradient();
360 ClosestPointProjection::calculateR(
Vec R) {
363 auto plastic_strain = getPlasticStrain();
364 auto internal_variables = getInternalVariables();
369 for (
int ii = 0; ii < 6; ii++) {
370 array[ii] = -plastic_strain[ii] + plasticStrain0[ii] +
371 deltaGamma * partialHSigma[ii];
373 for (
unsigned int ii = 0; ii < nbInternalVariables; ii++) {
374 array[6 + ii] = -internal_variables[ii] + internalVariables0[ii] +
375 deltaGamma * partialHFlux[ii];
377 array[6 + nbInternalVariables] = y;
378 CHKERR VecRestoreArray(
R, &array);
387 MatrixDouble a00 = prod(partial2HSigma, partialWStrainPlasticStrain);
389 for (
int ii = 0; ii < 6; ii++) {
390 for (
int jj = 0; jj < 6; jj++) {
391 dataA(ii, jj) = deltaGamma * a00(ii, jj);
393 for (
unsigned int jj = 0; jj < nbInternalVariables; jj++) {
394 dataA(ii, 6 + jj) = deltaGamma * a01(ii, jj);
398 for (
int ii = 0; ii < 6; ii++) {
399 dataA(ii, 6 + nbInternalVariables) = partialHSigma[ii];
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);
409 for (
unsigned int jj = 0; jj < nbInternalVariables; jj++) {
410 dataA(6 + ii, 6 + jj) += deltaGamma * a11(ii, jj);
412 dataA(6 + ii, 6 + ii) -= 1;
414 for (
unsigned int ii = 0; ii < nbInternalVariables; ii++) {
415 dataA(6 + ii, 6 + nbInternalVariables) = partialHFlux[ii];
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];
424 for (
unsigned int dd = 0;
dd < partial_y_flux_d.size();
dd++) {
425 dataA(6 + nbInternalVariables, 6 +
dd) = partial_y_flux_d[
dd];
427 dataA(6 + nbInternalVariables, 6 + nbInternalVariables) = 0;
428 CHKERR MatAssemblyBegin(
A, MAT_FINAL_ASSEMBLY);
429 CHKERR MatAssemblyEnd(
A, MAT_FINAL_ASSEMBLY);
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);
449 CHKERR SNESLineSearchSetType(line_search, SNESLINESEARCHBT);
452 CHKERR SNESAppendOptionsPrefix(sNes,
"cp_");
453 CHKERR SNESSetFromOptions(sNes);
459 ClosestPointProjection::solveClosetProjection() {
462 auto plastic_strain = getPlasticStrain();
463 auto internal_variables = getInternalVariables();
467 CHKERR VecGetArray(Chi, &array);
468 for (
int ii = 0; ii < 6; ii++) {
469 array[ii] = plastic_strain[ii];
471 int nb_internal_variables = nbInternalVariables;
472 for (
unsigned int ii = 0; ii < nb_internal_variables; ii++) {
473 array[6 + ii] = internal_variables[ii];
476 CHKERR VecRestoreArray(Chi, &array);
479 CHKERR SNESSolve(sNes, PETSC_NULL, Chi);
480 SNESConvergedReason reason;
481 CHKERR SNESGetConvergedReason(sNes, &reason);
484 CHKERR SNESGetIterationNumber(sNes, &its);
486 "ADOLCPlasticitySelf", Sev::warning,
487 "Plasticity Closet Point Projection - number of Newton iterations = %d",
489 plastic_strain = plasticStrain0;
490 internal_variables = internalVariables0;
494 CHKERR VecGetArray(Chi, &array);
495 for (
int ii = 0; ii < 6; ii++) {
496 plastic_strain[ii] = array[ii];
498 int nb_internal_variables = nbInternalVariables;
499 for (
unsigned int ii = 0; ii < nb_internal_variables; ii++) {
500 internal_variables[ii] = array[6 + ii];
502 deltaGamma = array[6 + nb_internal_variables];
503 CHKERR VecRestoreArray(Chi, &array);
509 ClosestPointProjection::consistentTangent() {
513 Ep.resize(6, 6,
false);
516 ep_row = deltaGamma * prod(partial2HSigma, C);
518 alpha_row = deltaGamma * prod(trans(partial2HSigmaFlux), C);
520 y_row = prod(C, partialYSigma);
521 const int n = 6 + nbInternalVariables;
523 for (
int dd = 0;
dd < 6;
dd++) {
528 for (
auto ii = 0; ii < 6; ii++) {
529 array[ii] = ep_row(ii,
dd);
531 for (
auto ii = 0; ii < nbInternalVariables; ii++) {
532 array[6 + ii] = alpha_row(ii,
dd);
534 array[
n] = y_row[
dd];
536 CHKERR VecRestoreArray(
R, &array);
541 CHKERR SNESGetKSP(sNes, &ksp);
543 CHKERR VecGetArray(dChi, &array);
545 for (
auto ii = 0; ii < 6; ii++) {
546 Ep(ii,
dd) = array[ii];
549 CHKERR VecRestoreArray(dChi, &array);
551 Cp.resize(6, 6,
false);
552 noalias(Cp) = prod(C, Ep);
553 Cep.resize(6, 6,
false);
554 noalias(Cep) = C + Cp;
568 CHKERR VecGetArrayRead(chi, &array);
569 for (
auto ii = 0; ii < 6; ii++) {
570 plastic_strain[ii] = array[ii];
573 internal_variables[ii] = array[6 + ii];
576 CHKERR VecRestoreArrayRead(chi, &array);
595 CHKERR VecGetArrayRead(chi, &array);
596 for (
auto ii = 0; ii < 6; ii++) {
597 plastic_strain[ii] = array[ii];
600 internal_variables[ii] = array[6 + ii];
603 CHKERR VecRestoreArrayRead(chi, &array);
617 activeVariablesW.swap(tmp_active);
618 gradientW.swap(tmp_gradient);
619 CHKERR evaluatePotentials();