v0.13.1
SmallStrainPlasticity.cpp
Go to the documentation of this file.
1/** \fiele SmallStrainPlasticity.cpp
2 * \brief Operators and data structures for small strain plasticity
3 * \ingroup small_strain_plasticity
4 *
5 * Implementation of small strain plasticity
6 *
7 * This file is part of MoFEM.
8 * MoFEM is free software: you can redistribute it and/or modify it under
9 * the terms of the GNU Lesser General Public License as published by the
10 * Free Software Foundation, either version 3 of the License, or (at your
11 * option) any later version.
12 *
13 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 * License for more details.
17 *
18 * You should have received a copy of the GNU Lesser General Public
19 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>.
20*/
21
22#include <MoFEM.hpp>
23using namespace MoFEM;
24
25#include <boost/numeric/ublas/vector_proxy.hpp>
26#include <boost/numeric/ublas/matrix.hpp>
27#include <boost/numeric/ublas/matrix_proxy.hpp>
28#include <boost/numeric/ublas/vector.hpp>
29#include <boost/numeric/ublas/symmetric.hpp>
30
31#include <adolc/adolc.h>
33
35 PetscFunctionBegin;
36 int def_length = 0;
37 rval = mField.get_moab().tag_get_handle(
38 "PLASTIC_STRAIN",def_length,MB_TYPE_DOUBLE,
39 thPlasticStrain,MB_TAG_CREAT|MB_TAG_SPARSE|MB_TAG_VARLEN,NULL
41 rval = mField.get_moab().tag_get_handle(
42 "INTERNAL_VARIABLES",def_length,MB_TYPE_DOUBLE,
43 thInternalVariables,MB_TAG_CREAT|MB_TAG_SPARSE|MB_TAG_VARLEN,NULL
45 PetscFunctionReturn(0);
46}
47
49 vector<VectorDouble > &values_at_gauss_pts,
50 vector<MatrixDouble > &gardient_at_gauss_pts
51):
53valuesAtGaussPts(values_at_gauss_pts),
54gradientAtGaussPts(gardient_at_gauss_pts),
55zeroAtType(MBVERTEX) {
56}
57
60) {
61 PetscFunctionBegin;
62 try {
63 int nb_dofs = data.getFieldData().size();
64 if(nb_dofs == 0) {
65 PetscFunctionReturn(0);
66 }
67 int nb_gauss_pts = data.getN().size1();
68 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
69 //initialize
70 if(type == zeroAtType) {
71 valuesAtGaussPts.resize(nb_gauss_pts);
72 gradientAtGaussPts.resize(nb_gauss_pts);
73 for(int gg = 0;gg<nb_gauss_pts;gg++) {
74 valuesAtGaussPts[gg].resize(rank,false);
75 gradientAtGaussPts[gg].resize(rank,3,false);
76 valuesAtGaussPts[gg].clear();
77 gradientAtGaussPts[gg].clear();
78 }
79 }
80 VectorDouble& values = data.getFieldData();
81 for(int gg = 0;gg<nb_gauss_pts;gg++) {
82 VectorDouble &values_at_gauss_pts = valuesAtGaussPts[gg];
83 MatrixDouble &gradient_at_gauss_pts = gradientAtGaussPts[gg];
84 VectorAdaptor N = data.getN(gg,nb_dofs/rank);
85 MatrixAdaptor diffN = data.getDiffN(gg,nb_dofs/rank);
86 for(int dd = 0;dd<nb_dofs/rank;dd++) {
87 for(int rr1 = 0;rr1<rank;rr1++) {
88 const double v = values[rank*dd+rr1];
89 values_at_gauss_pts[rr1] += N[dd]*v;
90 for(int rr2 = 0;rr2<3;rr2++) {
91 gradient_at_gauss_pts(rr1,rr2) += diffN(dd,rr2)*v;
92 }
93 }
94 }
95 }
96 // cerr << rowFieldName << endl;
97 // cerr << side << " " << type << endl;
98 // cerr << values << endl;
99 // cerr << valuesAtGaussPts[0] << endl;
100 // cerr << gradientAtGaussPts[0] << endl;
101 // cerr << endl;
102 } catch (const std::exception& ex) {
103 ostringstream ss;
104 ss << "throw in method: " << ex.what() << endl;
105 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
106 }
107 PetscFunctionReturn(0);
108}
109
111 const string field_name,CommonData &common_data
112):
114 field_name,common_data.dataAtGaussPts[field_name],common_data.gradAtGaussPts[field_name]
115) {
116}
117
119 string field_name,
120 CommonData &common_data
121):
123commonData(common_data) {
124}
125
128) {
129 PetscFunctionBegin;
130
131 try {
132 //do it only once, no need to repeat this for edges,faces or tets
133 if(type != MBVERTEX) PetscFunctionReturn(0);
134
135 int nb_dofs = data.getFieldData().size();
136 if(nb_dofs==0) PetscFunctionReturn(0);
137 int nb_gauss_pts = data.getN().size1();
138
139 for(int gg = 0;gg<nb_gauss_pts;gg++) {
140
141 if(commonData.deltaGamma[gg]>0) {
142
143 VectorAdaptor plastic_strain = VectorAdaptor(
144 6,ublas::shallow_array_adaptor<double>(6,&commonData.plasticStrainPtr[gg*6])
145 );
146 int nb_internal_variables = commonData.internalVariables[gg].size();
147 VectorAdaptor internal_variables = VectorAdaptor(
148 nb_internal_variables,
149 ublas::shallow_array_adaptor<double>(
150 nb_internal_variables,&commonData.internalVariablesPtr[gg*nb_internal_variables]
151 )
152 );
153
154 // cerr << plastic_strain << " " << commonData.plasticStrain[gg] << endl;
155 // cerr << internal_variables << " " << commonData.internalVariables[gg] << endl;
156
157 plastic_strain = commonData.plasticStrain[gg];
158 internal_variables = commonData.internalVariables[gg];
159 }
160
161 }
162
163 } catch (const std::exception& ex) {
164 ostringstream ss;
165 ss << "throw in method: " << ex.what() << endl;
166 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
167 }
168
169 PetscFunctionReturn(0);
170}
171
173 PetscFunctionBegin;
174 unsigned int nb_dofs = diffN.size1();
175 if(B.size2()!=3*nb_dofs) {
176 B.resize(6,3*nb_dofs,false);
177 }
178 B.clear();
179 for(int dd = 0;dd<nb_dofs;dd++) {
180 const double diff[] = {
181 diffN(dd,0), diffN(dd,1), diffN(dd,2)
182 };
183 const int dd3 = 3*dd;
184 for(int rr = 0;rr<3;rr++) {
185 B(rr,dd3+rr) = diff[rr];
186 }
187 // gamma_xy
188 B(3,dd3+0) = diff[1];
189 B(3,dd3+1) = diff[0];
190 // gamma_yz
191 B(4,dd3+1) = diff[2];
192 B(4,dd3+2) = diff[1];
193 // gamma_xz
194 B(5,dd3+0) = diff[2];
195 B(5,dd3+2) = diff[0];
196 }
197 PetscFunctionReturn(0);
198}
199
201 const MatrixAdaptor &diffN,MatrixDouble &volB,double alpha
202) {
203 PetscFunctionBegin;
204 int nb_dofs = diffN.size1();
205 for(int dd = 0;dd<nb_dofs;dd++) {
206 const double diff[] = {
207 diffN(dd,0), diffN(dd,1), diffN(dd,2)
208 };
209 const int dd3 = 3*dd;
210 for(int rr1 = 0;rr1<3;rr1++) {
211 for(int rr2 = 0;rr2<3;rr2++) {
212 volB(rr1,dd3+rr2) += alpha*diff[rr2]/3.;
213 }
214 }
215 }
216 PetscFunctionReturn(0);
217}
218
220 string field_name,
221 CommonData &common_data
222):
224commonData(common_data) {
225}
226
229) {
230 PetscFunctionBegin;
231
232 try {
233
234 int nb_dofs = data.getFieldData().size();
235 if(nb_dofs==0) PetscFunctionReturn(0);
236 int nb_gauss_pts = data.getN().size1();
237 nF.resize(nb_dofs,false);
238 nF.clear();
239
240 if(commonData.bBar) {
241 double v = 0;
242 volB.resize(6,nb_dofs,false);
243 volB.clear();
244 for(int gg = 0;gg<nb_gauss_pts;gg++) {
245 double val = getVolume()*getGaussPts()(3,gg);
246 v += val;
247 const MatrixAdaptor &diffN = data.getDiffN(gg,nb_dofs/3);
248 ierr = addVolumetricB(diffN,volB,val); CHKERRQ(ierr);
249 }
250 volB /= v;
251 }
252
253 for(int gg = 0;gg<nb_gauss_pts;gg++) {
254
255 double val = getVolume()*getGaussPts()(3,gg);
256 const MatrixAdaptor &diffN = data.getDiffN(gg,nb_dofs/3);
257 ierr = makeB(diffN,B); CHKERRQ(ierr);
258 if(commonData.bBar) {
259 ierr = addVolumetricB(diffN,B,-1); CHKERRQ(ierr);
260 B += volB;
261 }
262 // cerr << gg << endl;
263 // cerr << diffN << endl;
264 // cerr << B << endl;
265 // cerr << commonData.sTress[gg] << endl;
266 // cerr << prod(trans(B),commonData.sTress[gg]) << endl;
267
268 noalias(nF) += val*prod(trans(B),commonData.sTress[gg]);
269 // cerr << "nF " << nF << endl;
270
271 }
272
273 // cerr << side << " " << type << " " << data.getN().size1() << " "<< data.getN().size2() << endl;
274
275 int *indices_ptr;
276 indices_ptr = &data.getIndices()[0];
277
279 getFEMethod()->snes_f,
280 nb_dofs,
281 indices_ptr,
282 &nF[0],
283 ADD_VALUES
284 ); CHKERRQ(ierr);
285
286
287 } catch (const std::exception& ex) {
288 ostringstream ss;
289 ss << "throw in method: " << ex.what() << endl;
290 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
291 }
292
293 PetscFunctionReturn(0);
294}
295
297 string field_name,
298 CommonData &common_data
299):
301commonData(common_data) {
302 // sYmm = false;
303}
304
306 int row_side,int col_side,
307 EntityType row_type,EntityType col_type,
310) {
311 PetscFunctionBegin;
312 try {
313
314 int nb_dofs_row = row_data.getFieldData().size();
315 if(nb_dofs_row==0) PetscFunctionReturn(0);
316 int nb_dofs_col = col_data.getFieldData().size();
317 if(nb_dofs_col==0) PetscFunctionReturn(0);
318
319 K.resize(nb_dofs_row,nb_dofs_col,false);
320 K.clear();
321
322 int nb_gauss_pts = row_data.getN().size1();
323
324 if(commonData.bBar) {
325 double v = 0;
326 rowVolB.resize(6,nb_dofs_row,false);
327 rowVolB.clear();
328 colVolB.resize(6,nb_dofs_col,false);
329 colVolB.clear();
330 for(int gg = 0;gg<nb_gauss_pts;gg++) {
331 double val = getVolume()*getGaussPts()(3,gg);
332 v += val;
333 const MatrixAdaptor &diffN_row = row_data.getDiffN(gg,nb_dofs_row/3);
334 const MatrixAdaptor &diffN_col = col_data.getDiffN(gg,nb_dofs_col/3);
335 ierr = addVolumetricB(diffN_row,rowVolB,val); CHKERRQ(ierr);
336 ierr = addVolumetricB(diffN_col,colVolB,val); CHKERRQ(ierr);
337 }
338 rowVolB /= v;
339 colVolB /= v;
340 }
341
342 for(int gg = 0;gg<nb_gauss_pts;gg++) {
343
344 double val = getVolume()*getGaussPts()(3,gg);
345
346 const MatrixAdaptor &diffN_row = row_data.getDiffN(gg,nb_dofs_row/3);
347 const MatrixAdaptor &diffN_col = col_data.getDiffN(gg,nb_dofs_col/3);
348 ierr = makeB(diffN_row,rowB); CHKERRQ(ierr);
349 ierr = makeB(diffN_col,colB); CHKERRQ(ierr);
350 if(commonData.bBar) {
351 ierr = addVolumetricB(diffN_row,rowB,-1); CHKERRQ(ierr); // B-bar
352 ierr = addVolumetricB(diffN_col,colB,-1); CHKERRQ(ierr);
353 rowB += rowVolB;
354 colB += colVolB;
355 }
356
357 CB.resize(6,nb_dofs_col,false);
358 cblas_dgemm(
359 CblasRowMajor,CblasNoTrans,CblasNoTrans,
360 6,nb_dofs_col,6,
361 1.,
363 &colB(0,0),nb_dofs_col,
364 0,
365 &CB(0,0),nb_dofs_col
366 );
367 cblas_dgemm(
368 CblasRowMajor,CblasTrans,CblasNoTrans,
369 nb_dofs_row,nb_dofs_col,6,
370 val,
371 &rowB(0,0),nb_dofs_row,
372 &CB(0,0),nb_dofs_col,
373 1,
374 &K(0,0),nb_dofs_col
375 );
376 //noalias(CB) = prod(commonData.materialTangentOperator[gg],colB);
377 //noalias(K) += val*prod(trans(rowB),CB);
378
379 }
380
381 // cerr << K << endl;
382 // cerr << row_data.getIndices() << endl;
383 // cerr << col_data.getIndices() << endl;
384 // cerr << getFEMethod()->snes_A << endl;
385 // cerr << getFEMethod()->snes_B << endl;
386
387 int *row_indices_ptr,*col_indices_ptr;
388 row_indices_ptr = &row_data.getIndices()[0];
389 col_indices_ptr = &col_data.getIndices()[0];
390
392 getFEMethod()->snes_B,
393 nb_dofs_row,row_indices_ptr,
394 nb_dofs_col,col_indices_ptr,
395 &K(0,0),ADD_VALUES
396 ); CHKERRQ(ierr);
397
398 //is symmetric
399 if(row_side != col_side || row_type != col_type) {
400
401 transK.resize(nb_dofs_col,nb_dofs_row,false);
402 noalias(transK) = trans(K);
404 getFEMethod()->snes_B,
405 nb_dofs_col,col_indices_ptr,
406 nb_dofs_row,row_indices_ptr,
407 &transK(0,0),ADD_VALUES
408 ); CHKERRQ(ierr);
409
410 }
411
412 } catch (const std::exception& ex) {
413 ostringstream ss;
414 ss << "throw in method: " << ex.what() << endl;
415 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
416 }
417 PetscFunctionReturn(0);
418}
419
421 PetscFunctionBegin;
422
423 rval = mField.get_moab().tag_get_handle(
424 "PLASTIC_STRAIN",thPlasticStrain
425 ); CHKERRQ_MOAB(rval);
426 rval = mField.get_moab().tag_get_handle(
427 "INTERNAL_VARIABLES",thInternalVariables
428 ); CHKERRQ_MOAB(rval);
429 PetscFunctionReturn(0);
430}
431
433 const EntityHandle tet,const int nb_gauss_pts,const int nb_internal_variables
434) {
435 PetscFunctionBegin;
436
438 {
439 rval = mField.get_moab().tag_get_by_ptr(
441 );
442 if(rval != MB_SUCCESS || commonData.plasticStrainSize != 6*nb_gauss_pts) {
443 v.resize(6*nb_gauss_pts,false);
444 v.clear();
445 int tag_size[1];
446 tag_size[0] = v.size();
447 void const* tag_data[] = { &v[0] };
448 rval = mField.get_moab().tag_set_by_ptr(
449 thPlasticStrain,&tet,1,tag_data,tag_size
450 ); CHKERRQ_MOAB(rval);
451 rval = mField.get_moab().tag_get_by_ptr(
453 ); CHKERRQ_MOAB(rval);
454 }
455 }
456 if(nb_internal_variables>0) {
457 rval = mField.get_moab().tag_get_by_ptr(
459 );
460 if(rval != MB_SUCCESS || commonData.internalVariablesSize != nb_internal_variables*nb_gauss_pts) {
461 v.resize(nb_internal_variables*nb_gauss_pts,false);
462 v.clear();
463 int tag_size[1];
464 tag_size[0] = v.size();
465 void const* tag_data[] = { &v[0] };
466 rval = mField.get_moab().tag_set_by_ptr(
467 thInternalVariables,&tet,1,tag_data,tag_size
468 ); CHKERRQ_MOAB(rval);
469 rval = mField.get_moab().tag_get_by_ptr(
471 );
472 }
473 }
474 PetscFunctionReturn(0);
475}
476
479) {
480 PetscFunctionBegin;
481
482 try {
483
484 //do it only once, no need to repeat this for edges,faces or tets
485 if(type != MBVERTEX) PetscFunctionReturn(0);
486
487 int nb_dofs = data.getFieldData().size();
488 if(nb_dofs==0) PetscFunctionReturn(0);
489
490 int nb_gauss_pts = data.getN().size1();
491 int nb_internal_variables = cP.internalVariables.size();
492 if(!initCp) {
493 ierr = getTags(); CHKERRQ(ierr);
494 }
495 EntityHandle tet = getNumeredEntFiniteElementPtr()->getEnt();
496 ierr = setTagsData(tet,nb_gauss_pts,nb_internal_variables); CHKERRQ(ierr);
497
498 commonData.sTress.resize(nb_gauss_pts);
499 commonData.materialTangentOperator.resize(nb_gauss_pts);
500 commonData.plasticStrain.resize(nb_gauss_pts);
501 commonData.internalVariables.resize(nb_gauss_pts);
502 commonData.internalFluxes.resize(nb_gauss_pts);
503 commonData.deltaGamma.resize(nb_gauss_pts);
504
505 // B-bar
506 double ave_tr_strain = 0;
507 if(commonData.bBar) {
508 double v = 0;
509 for(int gg = 0;gg<nb_gauss_pts;gg++) {
510 double val = getVolume()*getGaussPts()(3,gg);
511 v += val;
512 for(int ii = 0;ii<3;ii++) {
513 ave_tr_strain += val*commonData.gradAtGaussPts[rowFieldName][gg](ii,ii)/3.;
514 }
515 // cerr << commonData.gradAtGaussPts[rowFieldName][gg] << endl;
516 }
517 ave_tr_strain /= v;
518 // cerr << "v " << v << " ave_tr_strain " << ave_tr_strain << endl;
519 }
520
521 for(int gg = 0;gg<nb_gauss_pts;gg++) {
522
523 VectorAdaptor plastic_strain = VectorAdaptor(
524 6,ublas::shallow_array_adaptor<double>(6,&commonData.plasticStrainPtr[gg*6])
525 );
526 VectorAdaptor internal_variables = VectorAdaptor(
527 nb_internal_variables,
528 ublas::shallow_array_adaptor<double>(
529 nb_internal_variables,&commonData.internalVariablesPtr[gg*nb_internal_variables]
530 )
531 );
532
533 cP.sTrain.resize(6,false);
534 {
535 double tr_strain = 0;
536 for(int ii = 0;ii<3;ii++) {
537 cP.sTrain[ii] = commonData.gradAtGaussPts[rowFieldName][gg](ii,ii);
538 tr_strain += cP.sTrain[ii]/3;
539 }
540 if(commonData.bBar) {
541 for(int ii = 0;ii<3;ii++) {
542 cP.sTrain[ii] += ave_tr_strain-tr_strain;
543 }
544 }
545 cP.sTrain[3] =
546 commonData.gradAtGaussPts[rowFieldName][gg](0,1)+
547 commonData.gradAtGaussPts[rowFieldName][gg](1,0);
548 cP.sTrain[4] =
549 commonData.gradAtGaussPts[rowFieldName][gg](1,2)+
550 commonData.gradAtGaussPts[rowFieldName][gg](2,1);
551 cP.sTrain[5] =
552 commonData.gradAtGaussPts[rowFieldName][gg](0,2)+
553 commonData.gradAtGaussPts[rowFieldName][gg](2,0);
554 // cerr << "Grad " << commonData.gradAtGaussPts[rowFieldName][gg] << endl;
555 // cerr << cP.sTrain << endl;
556 }
557
558
559 {
560 cP.plasticStrain0.resize(6,false);
561 noalias(cP.plasticStrain0) = plastic_strain;
562 cP.internalVariables0.resize(nb_internal_variables,false);
563 noalias(cP.internalVariables0) = internal_variables;
564 cP.plasticStrain.resize(6,false);
565 noalias(cP.plasticStrain) = plastic_strain;
566 cP.internalVariables.resize(nb_internal_variables,false);
567 noalias(cP.internalVariables) = internal_variables;
568 cP.deltaGamma = 0;
569 commonData.sTress[gg].clear();
570 commonData.internalFluxes[gg].clear();
571 }
572
573 if(!initCp) {
574 cP.gG = 0;
575 ierr = cP.evaluatePotentials(); CHKERRQ(ierr);
576 } else {
577 cP.gG = 1;
578 ierr = cP.setActiveVariablesW(); CHKERRQ(ierr);
579 if(
580 getFEMethod()->snes_ctx ==
581 SnesMethod::CTX_SNESSETJACOBIAN &&
582 (!isLinear || !hessianWCalculated)
583 ) {
584 ierr = cP.pLayW(); CHKERRQ(ierr);
585 hessianWCalculated = true;
586 } else {
587 ierr = cP.pLayW_NoHessian(); CHKERRQ(ierr);
588 }
589 ierr = cP.setActiveVariablesYH(); CHKERRQ(ierr);
590 ierr = cP.pLayY_NoGradient(); CHKERRQ(ierr);
591 // cerr << "C " << cP.C << endl;
592 }
593
594 double y = cP.y;
595 if(y>0) {
596 // ierr = PCFactorSetMatSolverPackage(cP.pC,MATSOLVERPETSC); CHKERRQ(ierr);
597 // ierr = KSPMonitorCancel(cP.kSp); CHKERRQ(ierr);
598 // ierr = SNESMonitorCancel(cP.sNes); CHKERRQ(ierr);
599 // ierr = SNESSetTolerances(cP.sNes,cP.tOl,1e-12,0,20,1000); CHKERRQ(ierr);
600 // ierr = SNESLineSearchSetType(cP.lineSearch,SNESLINESEARCHBT); CHKERRQ(ierr);
601 ierr = cP.solveColasetProjection(); CHKERRQ(ierr);
602 }
603
604 {
605 commonData.sTress[gg].resize(6,false);
606 noalias(commonData.sTress[gg]) = cP.sTress;
607 commonData.plasticStrain[gg].resize(6,false);
608 noalias(commonData.plasticStrain[gg]) = cP.plasticStrain;
609 commonData.internalVariables[gg].resize(nb_internal_variables,false);
610 noalias(commonData.internalVariables[gg]) = cP.internalVariables;
611 commonData.internalFluxes[gg].resize(nb_internal_variables,false);
612 noalias(commonData.internalFluxes[gg]) = cP.internalFluxes;
613 commonData.deltaGamma[gg] = cP.deltaGamma;
614 }
615
616 if(getFEMethod()->snes_ctx == SnesMethod::CTX_SNESSETJACOBIAN) {
617 int iter;
618 ierr = SNESGetIterationNumber(getFEMethod()->snes,&iter); CHKERRQ(ierr);
619 commonData.materialTangentOperator[gg].resize(6,6,false);
620 if(iter>0 && cP.deltaGamma>0) {
621 ierr = cP.consistentTangent(); CHKERRQ(ierr);
622 noalias(commonData.materialTangentOperator[gg]) = cP.Cep;
623 } else {
624 noalias(commonData.materialTangentOperator[gg]) = cP.C;
625 }
626 }
627
628 if(getFEMethod()->snes_ctx == SnesMethod::CTX_SNESSETFUNCTION) {
629 int iter;
630 ierr = SNESGetIterationNumber(getFEMethod()->snes,&iter); CHKERRQ(ierr);
631 if(iter>1) {
632 ierr = SNESSetLagJacobian(getFEMethod()->snes,1); CHKERRQ(ierr);
633 }
634 }
635
636
637 initCp = true;
638
639 }
640
641 } catch (const std::exception& ex) {
642 ostringstream ss;
643 ss << "throw in method: " << ex.what() << endl;
644 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
645 }
646
647 PetscFunctionReturn(0);
648}
649
651 PetscFunctionBegin;
652 activeVariablesW.resize(12+internalVariables.size(),false);
653 for(int ii = 0;ii<6;ii++) {
654 activeVariablesW[ii] = sTrain[ii];
655 }
656 for(int ii = 0;ii<6;ii++) {
657 activeVariablesW[6+ii] = plasticStrain[ii];
658 }
659 for(unsigned int ii = 0;ii<internalVariables.size();ii++) {
660 activeVariablesW[12+ii] = internalVariables[ii];
661 }
662 PetscFunctionReturn(0);
663}
664
666 PetscFunctionBegin;
667 activeVariablesYH.resize(6+internalVariables.size(),false);
668 for(int ii = 0;ii<6;ii++) {
669 activeVariablesYH[ii] = sTress[ii];
670 }
671 for(unsigned int ii = 0;ii<internalVariables.size();ii++) {
672 activeVariablesYH[6+ii] = internalFluxes[ii];
673 }
674 PetscFunctionReturn(0);
675}
676
678 PetscFunctionBegin;
679 trace_on(tAgs[0]);
680 {
681 //active variables
682 a_sTrain.resize(6,false);
683 for(int ii = 0;ii<6;ii++) {
684 a_sTrain[ii] <<= sTrain[ii];
685 }
686 a_plasticStrain.resize(6,false);
687 for(int ii = 0;ii<6;ii++) {
688 a_plasticStrain[ii] <<= plasticStrain[ii];
689 }
690 a_internalVariables.resize(internalVariables.size(),false);
691 for(unsigned int ii = 0;ii<internalVariables.size();ii++) {
692 a_internalVariables[ii] <<= internalVariables[ii];
693 }
694 //evaluete functions
695 ierr = freeHemholtzEnergy(); CHKERRQ(ierr);
696 //return variables
697 a_w >>= w;
698 }
699 trace_off();
700 PetscFunctionReturn(0);
701}
702
704 PetscFunctionBegin;
705 trace_on(tAgs[1]);
706 {
707 //active variables
708 a_sTress.resize(6,false);
709 for(int ii = 0;ii<6;ii++) {
710 a_sTress[ii] <<= sTress[ii];
711 }
712 a_internalFluxes.resize(internalVariables.size(),false);
713 for(unsigned int ii = 0;ii<internalVariables.size();ii++) {
714 a_internalFluxes[ii] <<= internalFluxes[ii];
715 }
716 //evaluete functions
717 ierr = yieldFunction(); CHKERRQ(ierr);
718 //return variables
719 a_y >>= y;
720 }
721 trace_off();
722 PetscFunctionReturn(0);
723}
724
726 PetscFunctionBegin;
727 trace_on(tAgs[2]);
728 {
729 //active variables
730 a_sTress.resize(6,false);
731 for(int ii = 0;ii<6;ii++) {
732 a_sTress[ii] <<= sTress[ii];
733 }
734 a_internalFluxes.resize(internalVariables.size(),false);
735 for(unsigned int ii = 0;ii<internalVariables.size();ii++) {
736 a_internalFluxes[ii] <<= internalFluxes[ii];
737 }
738 //evaluete functions
739 ierr = flowPotential(); CHKERRQ(ierr);
740 //return variables
741 a_h >>= h;
742 }
743 trace_off();
744 PetscFunctionReturn(0);
745}
746
748 PetscFunctionBegin;
749 int r;
750 int adloc_return_value = 0;
751 r = ::function(tAgs[0],1,activeVariablesW.size(),&activeVariablesW[0],&w);
752 if(r<adloc_return_value) {
753 SETERRQ(PETSC_COMM_SELF,MOFEM_OPERATION_UNSUCCESSFUL,"ADOL-C function evaluation with error");
754 }
755 gradientW.resize(activeVariablesW.size(),false);
756 r = gradient(tAgs[0],activeVariablesW.size(),&activeVariablesW[0],&gradientW[0]);
757 if(r<adloc_return_value) {
758 SETERRQ(PETSC_COMM_SELF,MOFEM_OPERATION_UNSUCCESSFUL,"ADOL-C function evaluation with error");
759 }
760 hessianW.resize(activeVariablesW.size(),activeVariablesW.size(),false);
761 hessianW.clear();
762 vector<double*> hessian_w(activeVariablesW.size());
763 for(unsigned int dd = 0;dd<activeVariablesW.size();dd++) {
764 hessian_w[dd] = &hessianW(dd,0);
765 }
766 r = hessian(tAgs[0],activeVariablesW.size(),&activeVariablesW[0],&hessian_w[0]);
767 if(r<adloc_return_value) {
768 SETERRQ(PETSC_COMM_SELF,MOFEM_OPERATION_UNSUCCESSFUL,"ADOL-C function evaluation with error");
769 }
770 sTress = VectorAdaptor(6,ublas::shallow_array_adaptor<double>(6,&gradientW[0]));
771 internalFluxes = VectorAdaptor(
772 internalVariables.size(),
773 ublas::shallow_array_adaptor<double>(
774 internalVariables.size(),
775 &gradientW[12]
776 )
777 );
778 C.resize(6,6,false);
779 for(int ii = 0;ii<6;ii++) {
780 for(int jj = 0;jj<=ii;jj++) {
781 C(ii,jj) = hessianW(ii,jj);
782 }
783 }
784 partialWStrainPlasticStrain.resize(6,6,false);
785 for(int ii = 0;ii<6;ii++) {
786 for(int jj = 0;jj<6;jj++) {
787 partialWStrainPlasticStrain(ii,jj) = hessianW(6+jj,ii);
788 }
789 }
790 D.resize(internalVariables.size(),internalVariables.size(),false);
791 for(unsigned int ii = 0;ii<internalVariables.size();ii++) {
792 for(unsigned int jj = 0;jj<=ii;jj++) {
793 D(ii,jj) = hessianW(12+ii,12+jj);
794 }
795 }
796 PetscFunctionReturn(0);
797}
799 PetscFunctionBegin;
800 int r;
801 int adloc_return_value = 0;
802 r = ::function(tAgs[0],1,activeVariablesW.size(),&activeVariablesW[0],&w);
803 if(r<adloc_return_value) {
804 SETERRQ(PETSC_COMM_SELF,MOFEM_OPERATION_UNSUCCESSFUL,"ADOL-C function evaluation with error");
805 }
806 gradientW.resize(activeVariablesW.size(),false);
807 r = gradient(tAgs[0],activeVariablesW.size(),&activeVariablesW[0],&gradientW[0]);
808 if(r<adloc_return_value) {
809 SETERRQ(PETSC_COMM_SELF,MOFEM_OPERATION_UNSUCCESSFUL,"ADOL-C function evaluation with error");
810 }
811 sTress = VectorAdaptor(6,ublas::shallow_array_adaptor<double>(6,&gradientW[0]));
812 internalFluxes = VectorAdaptor(
813 internalVariables.size(),
814 ublas::shallow_array_adaptor<double>(
815 internalVariables.size(),
816 &gradientW[12]
817 )
818 );
819 PetscFunctionReturn(0);
820}
821
823 PetscFunctionBegin;
824 int r;
825 int adloc_return_value = 0;
826 r = ::function(tAgs[1],1,activeVariablesYH.size(),&activeVariablesYH[0],&y);
827 if(r<adloc_return_value) {
828 SETERRQ(PETSC_COMM_SELF,MOFEM_OPERATION_UNSUCCESSFUL,"ADOL-C function evaluation with error");
829 }
830 gradientY.resize(activeVariablesYH.size());
831 r = gradient(tAgs[1],activeVariablesYH.size(),&activeVariablesYH[0],&gradientY[0]);
832 if(r<adloc_return_value) {
833 SETERRQ(PETSC_COMM_SELF,MOFEM_OPERATION_UNSUCCESSFUL,"ADOL-C function evaluation with error");
834 }
835 partialYSigma = VectorAdaptor(6,ublas::shallow_array_adaptor<double>(6,&gradientY[0]));
836 partialYFlux = VectorAdaptor(
837 internalVariables.size(),
838 ublas::shallow_array_adaptor<double>(
839 internalVariables.size(),
840 &gradientY[6]
841 )
842 );
843 PetscFunctionReturn(0);
844}
846 PetscFunctionBegin;
847 int r;
848 int adloc_return_value = 0;
849 r = ::function(tAgs[1],1,activeVariablesYH.size(),&activeVariablesYH[0],&y);
850 if(r<adloc_return_value) {
851 SETERRQ(PETSC_COMM_SELF,MOFEM_OPERATION_UNSUCCESSFUL,"ADOL-C function evaluation with error");
852 }
853 PetscFunctionReturn(0);
854}
856 PetscFunctionBegin;
857 int r;
858 int adloc_return_value = 0;
859 r = ::function(tAgs[2],1,activeVariablesYH.size(),&activeVariablesYH[0],&h);
860 if(r<adloc_return_value) {
861 SETERRQ(PETSC_COMM_SELF,MOFEM_OPERATION_UNSUCCESSFUL,"ADOL-C function evaluation with error");
862 }
863 gradientH.resize(activeVariablesYH.size());
864 r = gradient(tAgs[2],activeVariablesYH.size(),&activeVariablesYH[0],&gradientH[0]);
865 if(r<adloc_return_value) {
866 SETERRQ(PETSC_COMM_SELF,MOFEM_OPERATION_UNSUCCESSFUL,"ADOL-C function evaluation with error");
867 }
868 hessianH.resize(6+internalVariables.size(),6+internalVariables.size(),false);
869 hessianH.clear();
870 vector<double*> hessian_h(activeVariablesYH.size());
871 for(int dd = 0;dd<activeVariablesYH.size();dd++) {
872 hessian_h[dd] = &hessianH(dd,0);
873 }
874 r = hessian(tAgs[2],activeVariablesYH.size(),&activeVariablesYH[0],&hessian_h[0]);
875 if(r<adloc_return_value) {
876 SETERRQ(PETSC_COMM_SELF,MOFEM_OPERATION_UNSUCCESSFUL,"ADOL-C function evaluation with error");
877 }
878 partialHSigma = VectorAdaptor(6,ublas::shallow_array_adaptor<double>(6,&gradientH[0]));
879 partialHFlux = VectorAdaptor(
880 internalVariables.size(),
881 ublas::shallow_array_adaptor<double>(
882 internalVariables.size(),
883 &gradientH[6]
884 )
885 );
886 partialHFlux *= -1;
887 partial2HSigma.resize(6,6,false);
888 for(int ii = 0;ii<6;ii++) {
889 for(int jj = 0;jj<=ii;jj++) {
890 partial2HSigma(ii,jj) = hessianH(ii,jj);
891 }
892 }
893 partial2HFlux.resize(internalVariables.size(),internalVariables.size(),false);
894 for(unsigned int ii = 0;ii<internalVariables.size();ii++) {
895 for(unsigned int jj = 0;jj<=ii;jj++) {
896 partial2HFlux(ii,jj) = -hessianH(6+ii,6+jj);
897 }
898 }
899 partial2HSigmaFlux.resize(6,internalVariables.size(),false);
900 partial2HSigmaFlux.clear();
901 for(int ii = 0;ii<6;ii++) {
902 for(unsigned int jj = 0;jj<internalVariables.size();jj++) {
903 partial2HSigmaFlux(ii,jj) = -hessianH(6+jj,ii);
904 }
905 }
906 PetscFunctionReturn(0);
907}
908
910 PetscFunctionBegin;
911 PetscInt n;
912 n = 1+6+internalVariables.size();
913 dataA.resize(n,n,false);
914 ierr = MatCreateSeqDense(PETSC_COMM_SELF,n,n,&dataA(0,0),&A); CHKERRQ(ierr);
915 ierr = VecCreateSeq(PETSC_COMM_SELF,n,&R); CHKERRQ(ierr);
916 ierr = VecCreateSeq(PETSC_COMM_SELF,n,&Chi); CHKERRQ(ierr);
917 ierr = VecCreateSeq(PETSC_COMM_SELF,n,&dChi); CHKERRQ(ierr);
918 PetscFunctionReturn(0);
919}
920
922 PetscFunctionBegin;
923 ierr = MatDestroy(&A); CHKERRQ(ierr);
924 ierr = VecDestroy(&R); CHKERRQ(ierr);
925 ierr = VecDestroy(&Chi); CHKERRQ(ierr);
926 ierr = VecDestroy(&dChi); CHKERRQ(ierr);
927 PetscFunctionReturn(0);
928}
929
931 PetscFunctionBegin;
932 if(gG == 0) {
933 ierr = rEcordW(); CHKERRQ(ierr);
934 }
935 ierr = setActiveVariablesW(); CHKERRQ(ierr);
936 ierr = pLayW(); CHKERRQ(ierr);
937 if(gG == 0) {
938 ierr = rEcordY(); CHKERRQ(ierr);
939 ierr = rEcordH(); CHKERRQ(ierr);
940 }
941 ierr = setActiveVariablesYH(); CHKERRQ(ierr);
942 ierr = pLayY(); CHKERRQ(ierr);
943 ierr = pLayH(); CHKERRQ(ierr);
944 PetscFunctionReturn(0);
945}
946
948 PetscFunctionBegin;
949 // Residual
950 double *array;
951 ierr = VecGetArray(R,&array); CHKERRQ(ierr);
952 for(int ii = 0;ii<6;ii++) {
953 array[ii] = -plasticStrain[ii] + plasticStrain0[ii] + deltaGamma*partialHSigma[ii];
954 }
955 for(unsigned int ii = 0;ii<internalVariables.size();ii++) {
956 array[6+ii] = -internalVariables[ii] + internalVariables0[ii] + deltaGamma*partialHFlux[ii];
957 }
958 array[6+internalVariables.size()] = y;
959 ierr = VecRestoreArray(R,&array); CHKERRQ(ierr);
960 PetscFunctionReturn(0);
961}
962
964 PetscFunctionBegin;
965 // matrix A
966 ierr = MatZeroEntries(A); CHKERRQ(ierr);
967 // row 0
968 MatrixDouble a00 = prod(partial2HSigma,partialWStrainPlasticStrain);
969 MatrixDouble a01 = prod(partial2HSigmaFlux,D);
970 for(int ii = 0;ii<6;ii++) {
971 for(int jj = 0;jj<6;jj++) {
972 dataA(ii,jj) = deltaGamma*a00(ii,jj);
973 }
974 for(unsigned int jj = 0;jj<internalVariables.size();jj++) {
975 dataA(ii,6+jj) = deltaGamma*a01(ii,jj);
976 }
977 dataA(ii,ii) -= 1;
978 }
979 for(int ii = 0;ii<6;ii++) {
980 dataA(ii,6+internalVariables.size()) = partialHSigma[ii];
981 }
982 // row 1
983 MatrixDouble a11 = prod(partial2HFlux,D);
984 MatrixDouble a10 = prod(trans(partial2HSigmaFlux),partialWStrainPlasticStrain);
985 for(unsigned int ii = 0;ii<internalVariables.size();ii++) {
986 for(int jj = 0;jj<6;jj++) {
987 dataA(6+ii,jj) += deltaGamma*a10(ii,jj);
988 }
989 for(unsigned int jj = 0;jj<internalVariables.size();jj++) {
990 dataA(6+ii,6+jj) += deltaGamma*a11(ii,jj);
991 }
992 dataA(6+ii,6+ii) -= 1;
993 }
994 for(unsigned int ii = 0;ii<internalVariables.size();ii++) {
995 dataA(6+ii,6+internalVariables.size()) = partialHFlux[ii];
996 }
997 // row 3
998 VectorDouble partial_y_sigma_c = prod(trans(partialWStrainPlasticStrain),partialYSigma);
999 VectorDouble partial_y_flux_d = prod(trans(D),partialYFlux);
1000 for(unsigned int dd = 0;dd<partial_y_sigma_c.size();dd++) {
1001 dataA(6+internalVariables.size(),dd) = partial_y_sigma_c[dd];
1002 }
1003 for(unsigned int dd = 0;dd<partial_y_flux_d.size();dd++) {
1004 dataA(6+internalVariables.size(),6+dd) = partial_y_flux_d[dd];
1005 }
1006 dataA(6+internalVariables.size(),6+internalVariables.size()) = 0;
1007 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
1008 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
1009 PetscFunctionReturn(0);
1010}
1011
1013 PetscFunctionBegin;
1014 ierr = SNESCreate(PETSC_COMM_SELF,&sNes); CHKERRQ(ierr);
1015 ierr = SNESSetFunction(sNes,R,SmallStrainPlasticityfRes,(void*)this); CHKERRQ(ierr);
1016 ierr = SNESSetJacobian(sNes,A,A,SmallStrainPlasticityfJac,(void*)this); CHKERRQ(ierr);
1017 ierr = SNESGetKSP(sNes,&kSp); CHKERRQ(ierr);
1018 ierr = KSPGetPC(kSp,&pC); CHKERRQ(ierr);
1019 ierr = KSPSetType(kSp,KSPPREONLY); CHKERRQ(ierr);
1020 ierr = PCSetType(pC,PCLU); CHKERRQ(ierr);
1021 ierr = SNESSetTolerances(sNes,tOl,1e-12,0,20,1000); CHKERRQ(ierr);
1022 ierr = SNESGetLineSearch(sNes,&lineSearch); CHKERRQ(ierr); CHKERRQ(ierr);
1023 //ierr = SNESLineSearchSetType(lineSearch,SNESLINESEARCHBASIC); CHKERRQ(ierr);
1024 ierr = SNESLineSearchSetType(lineSearch,SNESLINESEARCHBT); CHKERRQ(ierr);
1025
1026 PetscFunctionReturn(0);
1027}
1028
1030 PetscFunctionBegin;
1031 ierr = SNESDestroy(&sNes); CHKERRQ(ierr);
1032 PetscFunctionReturn(0);
1033}
1034
1036 PetscFunctionBegin;
1037 {
1038 double *array;
1039 ierr = VecGetArray(Chi,&array); CHKERRQ(ierr);
1040 for(int ii = 0;ii<6;ii++) {
1041 array[ii] = plasticStrain[ii];
1042 }
1043 int nb_internal_variables = internalVariables.size();
1044 for(unsigned int ii = 0;ii<nb_internal_variables;ii++) {
1045 array[6+ii] = internalVariables[ii];
1046 }
1047 deltaGamma = 0;
1048 ierr = VecRestoreArray(Chi,&array); CHKERRQ(ierr);
1049 }
1050 ierr = SNESSolve(sNes,PETSC_NULL,Chi); CHKERRQ(ierr);
1051 SNESConvergedReason reason;
1052 ierr = SNESGetConvergedReason(sNes,&reason); CHKERRQ(ierr);
1053 if(reason < 0) {
1054 int its;
1055 ierr = SNESGetIterationNumber(sNes,&its); CHKERRQ(ierr);
1056// ierr = PetscPrintf(
1057// PETSC_COMM_SELF,
1058// "** WARNING Plasticity Closet Point Projection - number of Newton iterations = %D\n",
1059// its
1060// ); CHKERRQ(ierr);
1061 noalias(plasticStrain) = plasticStrain0;
1062 noalias(internalVariables) = internalVariables0;
1063 deltaGamma = 0;
1064 } else {
1065 double *array;
1066 ierr = VecGetArray(Chi,&array); CHKERRQ(ierr);
1067 for(int ii = 0;ii<6;ii++) {
1068 plasticStrain[ii] = array[ii];
1069 }
1070 int nb_internal_variables = internalVariables.size();
1071 for(unsigned int ii = 0;ii<nb_internal_variables;ii++) {
1072 internalVariables[ii] = array[6+ii];
1073 }
1074 deltaGamma = array[6+nb_internal_variables];
1075 ierr = VecRestoreArray(Chi,&array); CHKERRQ(ierr);
1076 }
1077 PetscFunctionReturn(0);
1078}
1079
1081 PetscFunctionBegin;
1082 ierr = evaluatePotentials(); CHKERRQ(ierr);
1083 ierr = cAlculateA(); CHKERRQ(ierr);
1084 Ep.resize(6,6,false);
1085 Lp.resize(internalVariables.size(),6,false);
1086 Yp.resize(6,false);
1087 MatrixDouble ep_row;
1088 ep_row = deltaGamma*prod(partial2HSigma,C);
1089 MatrixDouble alpha_row;
1090 alpha_row = deltaGamma*prod(trans(partial2HSigmaFlux),C);
1091 VectorDouble y_row;
1092 y_row = prod(trans(C),partialYSigma);
1093 const int n = 6+internalVariables.size();
1094 for(int dd = 0;dd<6;dd++) {
1095 ierr = VecZeroEntries(R); CHKERRQ(ierr);
1096 double *array;
1097 ierr = VecGetArray(R,&array); CHKERRQ(ierr);
1098 {
1099 for(int ii = 0;ii<6;ii++) {
1100 array[ii] = ep_row(ii,dd);
1101 }
1102 for(unsigned int ii = 0;ii<internalVariables.size();ii++) {
1103 array[6+ii] = alpha_row(ii,dd);
1104 }
1105 array[n] = y_row[dd];
1106 }
1107 ierr = VecRestoreArray(R,&array); CHKERRQ(ierr);
1108 ierr = VecAssemblyBegin(R); CHKERRQ(ierr);
1109 ierr = VecAssemblyEnd(R); CHKERRQ(ierr);
1110 ierr = KSPSolve(kSp,R,dChi); CHKERRQ(ierr);
1111 ierr = VecGetArray(dChi,&array); CHKERRQ(ierr);
1112 {
1113 for(int ii = 0;ii<6;ii++) {
1114 Ep(ii,dd) = array[ii];
1115 }
1116 for(unsigned int ii = 0;ii<internalVariables.size();ii++) {
1117 Lp(ii,dd) = array[6+ii];
1118 }
1119 Yp[dd] = array[n];
1120 }
1121 ierr = VecRestoreArray(dChi,&array); CHKERRQ(ierr);
1122 }
1123 Cp.resize(6,6,false);
1124 noalias(Cp) = prod(C,Ep);
1125 Cep.resize(6,6,false);
1126 noalias(Cep) = C+Cp;
1127 // cout << endl;
1128 // cout << "A " << dataA << endl;
1129 // cout << "C " << C << endl;
1130 // cout << "Ep " << Ep << endl;
1131 // // cout << "Lp " << Lp << endl;
1132 // // cout << "Yp " << Yp << endl;
1133 // cout << "Cp " << Cp << endl;
1134 // cout << "Cep " << Cep << endl;
1135 PetscFunctionReturn(0);
1136}
1137
1138
1139PetscErrorCode SmallStrainPlasticityfRes(SNES snes,Vec chi,Vec r,void *ctx) {
1140 PetscFunctionBegin;
1143
1144 //ierr = VecView(chi,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
1145 {
1146 #if PETSC_VERSION_GE(3,6,0)
1147 const double *array;
1148 ierr = VecGetArrayRead(chi,&array); CHKERRQ(ierr);
1149 #else
1150 const double *array;
1151 ierr = VecGetArray(chi,&array); CHKERRQ(ierr);
1152 #endif
1153 for(int ii = 0;ii<6;ii++) {
1154 cp->plasticStrain[ii] = array[ii];
1155 }
1156 for(unsigned int ii = 0;ii<cp->internalVariables.size();ii++) {
1157 cp->internalVariables[ii] = array[6+ii];
1158 }
1159 cp->deltaGamma = array[6+cp->internalVariables.size()];
1160 #if PETSC_VERSION_GE(3,6,0)
1161 ierr = VecRestoreArrayRead(chi,&array); CHKERRQ(ierr);
1162 #else
1163 ierr = VecRestoreArray(chi,&array); CHKERRQ(ierr);
1164 #endif
1165 }
1166 ierr = cp->evaluatePotentials(); CHKERRQ(ierr);
1167 ierr = cp->cAlculateR(r); CHKERRQ(ierr);
1168 // ierr = VecView(r,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
1169 // cout << "sTress " << cp->sTress << endl;
1170 // cout << "internalVariables " << cp->internalVariables << endl;
1171 // cout << "C " << cp->C << endl;
1172 // cout << "D " << cp->D << endl;
1173 // cout << "y " << cp->y << endl;
1174 // cout << "h " << cp->h << endl;
1175 // cout << "partialYSigma " << cp->partialYSigma << endl;
1176 // cout << "partialYFlux " << cp->partialYFlux << endl;
1177 // cout << "partialHSigma " << cp->partialHSigma << endl;
1178 // cout << "partialHFlux " << cp->partialHFlux << endl;
1179 // cout << "partial2HSigma " << cp->partial2HSigma << endl;
1180 // cout << "partial2HFlux " << cp->partial2HFlux << endl;
1181 // cout << "partial2HSigmaFlux " << cp->partial2HSigmaFlux << endl;
1182 // cout << "C " << cp->C << endl;
1183 // cout << "D " << cp->D << endl;
1184 // cout << "partialWStrainPlasticStrain " << cp->partialWStrainPlasticStrain << endl;
1185 // cout << "plasticStrain " << cp->plasticStrain << endl;
1186 // cout << "internalVariables " << cp->internalVariables << endl;
1187 // cout << "delta plasticStrain " << cp->plasticStrain-cp->plasticStrain0 << endl;
1188 // cout << "delta internalVariables " << cp->internalVariables-cp->internalVariables << endl;
1189 PetscFunctionReturn(0);
1190}
1191
1192PetscErrorCode SmallStrainPlasticityfJac(SNES snes,Vec chi,Mat A,Mat,void *ctx) {
1193 PetscFunctionBegin;
1196
1197 {
1198 #if PETSC_VERSION_GE(3,6,0)
1199 const double *array;
1200 ierr = VecGetArrayRead(chi,&array); CHKERRQ(ierr);
1201 #else
1202 double *array;
1203 ierr = VecGetArray(chi,&array); CHKERRQ(ierr);
1204 #endif
1205 for(int ii = 0;ii<6;ii++) {
1206 cp->plasticStrain[ii] = array[ii];
1207 }
1208 for(unsigned int ii = 0;ii<cp->internalVariables.size();ii++) {
1209 cp->internalVariables[ii] = array[6+ii];
1210 }
1211 cp->deltaGamma = array[6+cp->internalVariables.size()];
1212 #if PETSC_VERSION_GE(3,6,0)
1213 ierr = VecRestoreArrayRead(chi,&array); CHKERRQ(ierr);
1214 #else
1215 ierr = VecRestoreArray(chi,&array); CHKERRQ(ierr);
1216 #endif
1217 }
1218 ierr = cp->evaluatePotentials(); CHKERRQ(ierr);
1219 ierr = cp->cAlculateA(); CHKERRQ(ierr);
1220 PetscFunctionReturn(0);
1221}
static Index< 'K', 3 > K
ForcesAndSourcesCore::UserDataOperator UserDataOperator
PetscErrorCode SmallStrainPlasticityfRes(SNES snes, Vec chi, Vec r, void *ctx)
PetscErrorCode SmallStrainPlasticityfJac(SNES snes, Vec chi, Mat A, Mat, void *ctx)
Operators and data structures for small strain plasticity.
EntitiesFieldData::EntData EntData
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:454
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
FTensor::Index< 'n', SPACE_DIM > n
const double v
phase velocity of light in medium (cm/ns)
const FTensor::Tensor2< T, Dim, Dim > Vec
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
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:115
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
Definition: Types.hpp:132
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
double w(const double g, const double t)
const double D
diffusivity
const double r
rate factor
double A
double h
constexpr auto field_name
const int N
Definition: speed_test.cpp:3
virtual moab::Interface & get_moab()=0
common data used by volume elements
map< string, vector< MatrixDouble > > gradAtGaussPts
vector< MatrixDouble > materialTangentOperator
PetscErrorCode makeB(const MatrixAdaptor &diffN, MatrixDouble &B)
PetscErrorCode addVolumetricB(const MatrixAdaptor &diffN, MatrixDouble &volB, double alpha)
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
OpAssembleLhs(string field_name, CommonData &common_data)
OpAssembleRhs(string field_name, CommonData &common_data)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
PetscErrorCode setTagsData(const EntityHandle tet, const int nb_gauss_pts, const int nb_internal_variables)
OpGetCommonDataAtGaussPts(const string field_name, CommonData &common_data)
OpGetDataAtGaussPts(const string field_name, vector< VectorDouble > &values_at_gauss_pts, vector< MatrixDouble > &gardient_at_gauss_pts)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
operator calculating field data and gradients
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpUpdate(string field_name, CommonData &common_data)
PetscErrorCode createInternalVariablesTag()