v0.15.0
Loading...
Searching...
No Matches
NonLinearElasticElement.cpp
Go to the documentation of this file.
1/**
2 * \brief Operators and data structures for nonlinear elastic material
3 *
4 * Implementation of nonlinear elastic element.
5 */
6
7
8
9#include <MoFEM.hpp>
10
11using namespace MoFEM;
13
14#include <adolc/adolc.h>
16
18 : VolumeElementForcesAndSourcesCore(m_field), A(PETSC_NULLPTR), F(PETSC_NULLPTR),
19 addToRule(1) {
20
21 auto create_vec = [&]() {
22 constexpr int ghosts[] = {0};
23 if (mField.get_comm_rank() == 0) {
24 return createVectorMPI(mField.get_comm(), 1, 1);
25 } else {
26 return createVectorMPI(mField.get_comm(), 0, 1);
27 }
28 };
29
30 V = create_vec();
31}
32
34 return 2 * (order - 1) + addToRule;
35};
36
39
40 CHKERR VolumeElementForcesAndSourcesCore::preProcess();
41
42 if (A != PETSC_NULLPTR) {
43 snes_B = A;
44 }
45
46 if (F != PETSC_NULLPTR) {
47 snes_f = F;
48 }
49
50 switch (snes_ctx) {
51 case CTX_SNESNONE:
52 CHKERR VecZeroEntries(V);
53 break;
54 default:
55 break;
56 }
57
59}
60
63
64 switch (snes_ctx) {
65 case CTX_SNESNONE:
66 CHKERR VecAssemblyBegin(V);
67 CHKERR VecAssemblyEnd(V);
68 CHKERR VecSum(V, &eNergy);
69 break;
70 default:
71 break;
72 }
73
74 CHKERR VolumeElementForcesAndSourcesCore::postProcess();
75
77}
78
80 short int tag)
81 : feRhs(m_field), feLhs(m_field), feEnergy(m_field), mField(m_field),
82 tAg(tag) {}
83
85 const std::string field_name,
86 std::vector<VectorDouble> &values_at_gauss_pts,
87 std::vector<MatrixDouble> &gardient_at_gauss_pts)
90 valuesAtGaussPts(values_at_gauss_pts),
91 gradientAtGaussPts(gardient_at_gauss_pts), zeroAtType(MBVERTEX) {}
92
94 int side, EntityType type, EntitiesFieldData::EntData &data) {
96
97 const int nb_dofs = data.getFieldData().size();
98 const int nb_base_functions = data.getN().size2();
99 if (nb_dofs == 0) {
101 }
102 const int nb_gauss_pts = data.getN().size1();
103 const int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
104
105 // initialize
106 if (type == zeroAtType) {
107 valuesAtGaussPts.resize(nb_gauss_pts);
108 gradientAtGaussPts.resize(nb_gauss_pts);
109 for (int gg = 0; gg != nb_gauss_pts; gg++) {
110 valuesAtGaussPts[gg].resize(rank, false);
111 gradientAtGaussPts[gg].resize(rank, 3, false);
112 }
113 for (int gg = 0; gg != nb_gauss_pts; gg++) {
114 valuesAtGaussPts[gg].clear();
115 gradientAtGaussPts[gg].clear();
116 }
117 }
118
119 auto base_function = data.getFTensor0N();
120 auto diff_base_functions = data.getFTensor1DiffN<3>();
121 FTensor::Index<'i', 3> i;
122 FTensor::Index<'j', 3> j;
123
124 if (rank == 1) {
125
126 for (int gg = 0; gg != nb_gauss_pts; gg++) {
127 auto field_data = data.getFTensor0FieldData();
128 double &val = valuesAtGaussPts[gg][0];
129 FTensor::Tensor1<double *, 3> grad(&gradientAtGaussPts[gg](0, 0),
130 &gradientAtGaussPts[gg](0, 1),
131 &gradientAtGaussPts[gg](0, 2));
132 int bb = 0;
133 for (; bb != nb_dofs; bb++) {
134 val += base_function * field_data;
135 grad(i) += diff_base_functions(i) * field_data;
136 ++diff_base_functions;
137 ++base_function;
138 ++field_data;
139 }
140 for (; bb != nb_base_functions; bb++) {
141 ++diff_base_functions;
142 ++base_function;
143 }
144 }
145
146 } else if (rank == 3) {
147
148 for (int gg = 0; gg != nb_gauss_pts; gg++) {
149 auto field_data = data.getFTensor1FieldData<3>();
150 FTensor::Tensor1<double *, 3> values(&valuesAtGaussPts[gg][0],
151 &valuesAtGaussPts[gg][1],
152 &valuesAtGaussPts[gg][2]);
154 &gradientAtGaussPts[gg](0, 0), &gradientAtGaussPts[gg](0, 1),
155 &gradientAtGaussPts[gg](0, 2), &gradientAtGaussPts[gg](1, 0),
156 &gradientAtGaussPts[gg](1, 1), &gradientAtGaussPts[gg](1, 2),
157 &gradientAtGaussPts[gg](2, 0), &gradientAtGaussPts[gg](2, 1),
158 &gradientAtGaussPts[gg](2, 2));
159 int bb = 0;
160 for (; bb != nb_dofs / 3; bb++) {
161 values(i) += base_function * field_data(i);
162 gradient(i, j) += field_data(i) * diff_base_functions(j);
163 ++diff_base_functions;
164 ++base_function;
165 ++field_data;
166 }
167 for (; bb != nb_base_functions; bb++) {
168 ++diff_base_functions;
169 ++base_function;
170 }
171 }
172
173 } else {
174 // FIXME: THat part is inefficient
175 VectorDouble &values = data.getFieldData();
176 for (int gg = 0; gg < nb_gauss_pts; gg++) {
177 VectorAdaptor N = data.getN(gg, nb_dofs / rank);
178 MatrixAdaptor diffN = data.getDiffN(gg, nb_dofs / rank);
179 for (int dd = 0; dd < nb_dofs / rank; dd++) {
180 for (int rr1 = 0; rr1 < rank; rr1++) {
181 valuesAtGaussPts[gg][rr1] += N[dd] * values[rank * dd + rr1];
182 for (int rr2 = 0; rr2 < 3; rr2++) {
183 gradientAtGaussPts[gg](rr1, rr2) +=
184 diffN(dd, rr2) * values[rank * dd + rr1];
185 }
186 }
187 }
188 }
189 }
190
192}
193
195 const std::string field_name, CommonData &common_data)
196 : OpGetDataAtGaussPts(field_name, common_data.dataAtGaussPts[field_name],
197 common_data.gradAtGaussPts[field_name]) {}
198
201 BlockData &data, CommonData &common_data,
202 int tag, bool jacobian, bool ale,
203 bool field_disp)
206 dAta(data), commonData(common_data), tAg(tag), adlocReturnValue(0),
207 jAcobian(jacobian), fUnction(!jacobian), aLe(ale), fieldDisp(field_disp) {
208
209}
210
213 const int gg) {
215
216 CHKERR dAta.materialAdoublePtr->calculateP_PiolaKirchhoffI(
217 dAta, getNumeredEntFiniteElementPtr());
218
219 if (aLe) {
220 auto &t_P = dAta.materialAdoublePtr->t_P;
221 auto &t_invH = dAta.materialAdoublePtr->t_invH;
222 t_P(i, j) = t_P(i, k) * t_invH(j, k);
223 t_P(i, j) *= dAta.materialAdoublePtr->detH;
224 }
225
226 commonData.sTress[gg].resize(3, 3, false);
227 for (int dd1 = 0; dd1 < 3; dd1++) {
228 for (int dd2 = 0; dd2 < 3; dd2++) {
229 dAta.materialAdoublePtr->P(dd1, dd2) >>=
230 (commonData.sTress[gg])(dd1, dd2);
231 }
232 }
233
235}
236
239 const int gg) {
241
242 trace_on(tAg, 0);
243
244 dAta.materialAdoublePtr->F.resize(3, 3, false);
245
246 if (!aLe) {
247
248 nbActiveVariables = 0;
249 for (int dd1 = 0; dd1 < 3; dd1++) {
250 for (int dd2 = 0; dd2 < 3; dd2++) {
251 dAta.materialAdoublePtr->F(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
252 if (fieldDisp) {
253 if (dd1 == dd2) {
254 dAta.materialAdoublePtr->F(dd1, dd2) += 1;
255 }
256 }
257 nbActiveVariables++;
258 }
259 }
260
261 } else {
262
263 nbActiveVariables = 0;
264
265 dAta.materialAdoublePtr->h.resize(3, 3, false);
266 for (int dd1 = 0; dd1 < 3; dd1++) {
267 for (int dd2 = 0; dd2 < 3; dd2++) {
268 dAta.materialAdoublePtr->h(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
269 nbActiveVariables++;
270 }
271 }
272
273 dAta.materialAdoublePtr->H.resize(3, 3, false);
274 for (int dd1 = 0; dd1 < 3; dd1++) {
275 for (int dd2 = 0; dd2 < 3; dd2++) {
276 dAta.materialAdoublePtr->H(dd1, dd2) <<= (*ptrH)[gg](dd1, dd2);
277 nbActiveVariables++;
278 }
279 }
280
281 dAta.materialAdoublePtr->detH = determinantTensor3by3(dAta.materialAdoublePtr->H);
282 dAta.materialAdoublePtr->invH.resize(3, 3, false);
283 CHKERR invertTensor3by3(dAta.materialAdoublePtr->H,
284 dAta.materialAdoublePtr->detH,
285 dAta.materialAdoublePtr->invH);
286
287 auto &t_F = dAta.materialAdoublePtr->t_F;
288 auto &t_h = dAta.materialAdoublePtr->t_h;
289 auto &t_invH = dAta.materialAdoublePtr->t_invH;
290
291 t_F(i, j) = t_h(i, k) * t_invH(k, j);
292
293 }
294
295 CHKERR dAta.materialAdoublePtr->setUserActiveVariables(nbActiveVariables);
296 CHKERR calculateStress(gg);
297
298 trace_off();
299
301}
302
306
307 int r;
308
309 if (fUnction) {
310 commonData.sTress[gg].resize(3, 3, false);
311 // play recorder for values
312 r = ::function(tAg, 9, nbActiveVariables, &activeVariables[0],
313 &commonData.sTress[gg](0, 0));
314 if (r < adlocReturnValue) { // function is locally analytic
315 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
316 "ADOL-C function evaluation with error r = %d", r);
317 }
318 }
319
320 if (jAcobian) {
321 commonData.jacStress[gg].resize(9, nbActiveVariables, false);
322 double *jac_ptr[] = {
323 &(commonData.jacStress[gg](0, 0)), &(commonData.jacStress[gg](1, 0)),
324 &(commonData.jacStress[gg](2, 0)), &(commonData.jacStress[gg](3, 0)),
325 &(commonData.jacStress[gg](4, 0)), &(commonData.jacStress[gg](5, 0)),
326 &(commonData.jacStress[gg](6, 0)), &(commonData.jacStress[gg](7, 0)),
327 &(commonData.jacStress[gg](8, 0))};
328 // play recorder for jacobians
329 r = jacobian(tAg, 9, nbActiveVariables, &activeVariables[0], jac_ptr);
330 if (r < adlocReturnValue) {
331 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
332 "ADOL-C function evaluation with error");
333 }
334 }
335
337}
338
340 int row_side, EntityType row_type,
341 EntitiesFieldData::EntData &row_data) {
343
344 // do it only once, no need to repeat this for edges,faces or tets
345 if (row_type != MBVERTEX)
347
348 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
349 dAta.tEts.end()) {
351 }
352
353 int nb_dofs = row_data.getFieldData().size();
354 if (nb_dofs == 0)
356 dAta.materialAdoublePtr->commonDataPtr = &commonData;
357 dAta.materialAdoublePtr->opPtr = this;
358
359 int nb_gauss_pts = row_data.getN().size1();
360 commonData.sTress.resize(nb_gauss_pts);
361 commonData.jacStress.resize(nb_gauss_pts);
362
364 if (aLe) {
366 }
367
368 for (int gg = 0; gg != nb_gauss_pts; gg++) {
369
370 dAta.materialAdoublePtr->gG = gg;
371
372 // Record tag and calculate stress
373 if (recordTagForIntegrationPoint(gg)) {
374 CHKERR recordTag(gg);
375 }
376
377 // Set active variables vector
378 if (jAcobian || (!recordTagForIntegrationPoint(gg))) {
379 activeVariables.resize(nbActiveVariables, false);
380 if (!aLe) {
381 for (int dd1 = 0; dd1 < 3; dd1++) {
382 for (int dd2 = 0; dd2 < 3; dd2++) {
383 activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
384 }
385 }
386 } else {
387 for (int dd1 = 0; dd1 < 3; dd1++) {
388 for (int dd2 = 0; dd2 < 3; dd2++) {
389 activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
390 }
391 }
392 for (int dd1 = 0; dd1 < 3; dd1++) {
393 for (int dd2 = 0; dd2 < 3; dd2++) {
394 activeVariables(9 + dd1 * 3 + dd2) = (*ptrH)[gg](dd1, dd2);
395 }
396 }
397 }
398 CHKERR dAta.materialAdoublePtr->setUserActiveVariables(activeVariables);
399
400 // Play tag and calculate stress or tangent
401 if (jAcobian || (!recordTagForIntegrationPoint(gg))) {
402 CHKERR playTag(gg);
403 }
404 }
405 }
406
408}
409
411 const std::string
412 field_name, ///< field name for spatial positions or displacements
413 BlockData &data, CommonData &common_data, int tag, bool gradient,
414 bool hessian, bool ale, bool field_disp)
417 dAta(data), commonData(common_data), tAg(tag), gRadient(gradient),
418 hEssian(hessian), aLe(ale), fieldDisp(field_disp) {}
419
423 CHKERR dAta.materialAdoublePtr->calculateElasticEnergy(
424 dAta, getNumeredEntFiniteElementPtr());
425 dAta.materialAdoublePtr->eNergy >>= commonData.eNergy[gg];
427}
428
432
433 trace_on(tAg, 0);
434
435 if (!aLe) {
436
437 nbActiveVariables = 0;
438 for (int dd1 = 0; dd1 < 3; dd1++) {
439 for (int dd2 = 0; dd2 < 3; dd2++) {
440 dAta.materialAdoublePtr->F(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
441 if (fieldDisp) {
442 if (dd1 == dd2) {
443 dAta.materialAdoublePtr->F(dd1, dd2) += 1;
444 }
445 }
446 nbActiveVariables++;
447 }
448 }
449
450 } else {
451
452 nbActiveVariables = 0;
453
454 dAta.materialAdoublePtr->h.resize(3, 3, false);
455 for (int dd1 = 0; dd1 < 3; dd1++) {
456 for (int dd2 = 0; dd2 < 3; dd2++) {
457 dAta.materialAdoublePtr->h(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
458 nbActiveVariables++;
459 }
460 }
461
462 dAta.materialAdoublePtr->H.resize(3, 3, false);
463 for (int dd1 = 0; dd1 < 3; dd1++) {
464 for (int dd2 = 0; dd2 < 3; dd2++) {
465 dAta.materialAdoublePtr->H(dd1, dd2) <<= (*ptrH)[gg](dd1, dd2);
466 nbActiveVariables++;
467 }
468 }
469
470 dAta.materialAdoublePtr->detH = determinantTensor3by3(dAta.materialAdoublePtr->H);
471 dAta.materialAdoublePtr->invH.resize(3, 3, false);
472 CHKERR invertTensor3by3(dAta.materialAdoublePtr->H,
473 dAta.materialAdoublePtr->detH,
474 dAta.materialAdoublePtr->invH);
475
476 auto &t_F = dAta.materialAdoublePtr->t_F;
477 auto &t_h = dAta.materialAdoublePtr->t_h;
478 auto &t_invH = dAta.materialAdoublePtr->t_invH;
479
480 t_F(i, j) = t_h(i, k) * t_invH(k, j);
481
482 }
483
484 CHKERR dAta.materialAdoublePtr->setUserActiveVariables(nbActiveVariables);
486
487 trace_off();
488
490}
491
495
496 if (gRadient) {
497 commonData.jacEnergy[gg].resize(nbActiveVariables, false);
498 int r = ::gradient(tAg, nbActiveVariables, &activeVariables[0],
499 &commonData.jacEnergy[gg][0]);
500 if (r < 0) {
501 // That means that energy function is not smooth and derivative
502 // can not be calculated,
503 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
504 "ADOL-C function evaluation with error");
505 }
506 }
507
508 if (hEssian) {
509 commonData.hessianEnergy[gg].resize(nbActiveVariables * nbActiveVariables,
510 false);
511 double *H[nbActiveVariables];
512 for (int n = 0; n != nbActiveVariables; n++) {
513 H[n] = &(commonData.hessianEnergy[gg][n * nbActiveVariables]);
514 }
515 int r = ::hessian(tAg, nbActiveVariables, &*activeVariables.begin(), H);
516 if (r < 0) {
517 // That means that energy function is not smooth and derivative
518 // can not be calculated,
519 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
520 "ADOL-C function evaluation with error");
521 }
522 }
523
525}
526
528 int row_side, EntityType row_type,
529 EntitiesFieldData::EntData &row_data) {
531
532 // do it only once, no need to repeat this for edges,faces or tets
533 if (row_type != MBVERTEX)
535
536 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
537 dAta.tEts.end()) {
539 }
540
541 int nb_dofs = row_data.getFieldData().size();
542 if (nb_dofs == 0)
544 dAta.materialAdoublePtr->commonDataPtr = &commonData;
545 dAta.materialAdoublePtr->opPtr = this;
546
547 int nb_gauss_pts = row_data.getN().size1();
548 commonData.eNergy.resize(nb_gauss_pts);
549 commonData.jacEnergy.resize(nb_gauss_pts);
550
552 if (aLe) {
554 }
555
556 for (int gg = 0; gg != nb_gauss_pts; gg++) {
557
558 dAta.materialAdoublePtr->gG = gg;
559
560 // Record tag and calualte stress
561 if (recordTagForIntegrationPoint(gg)) {
562 CHKERR recordTag(gg);
563 }
564
565 activeVariables.resize(nbActiveVariables, false);
566 if (!aLe) {
567 for (int dd1 = 0; dd1 < 3; dd1++) {
568 for (int dd2 = 0; dd2 < 3; dd2++) {
569 activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
570 }
571 }
572 } else {
573 for (int dd1 = 0; dd1 < 3; dd1++) {
574 for (int dd2 = 0; dd2 < 3; dd2++) {
575 activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
576 }
577 }
578 for (int dd1 = 0; dd1 < 3; dd1++) {
579 for (int dd2 = 0; dd2 < 3; dd2++) {
580 activeVariables(9 + dd1 * 3 + dd2) = (*ptrH)[gg](dd1, dd2);
581 }
582 }
583 }
584 CHKERR dAta.materialAdoublePtr->setUserActiveVariables(activeVariables);
585
586 // Play tag and calculate stress or tangent
587 CHKERR playTag(gg);
588 }
589
591}
592
598
600 int row_side, EntityType row_type,
601 EntitiesFieldData::EntData &row_data) {
603
604 int nb_dofs = row_data.getIndices().size();
605 int *indices_ptr = &row_data.getIndices()[0];
606 if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
607 iNdices.resize(nb_dofs, false);
608 noalias(iNdices) = row_data.getIndices();
609 indices_ptr = &iNdices[0];
610 VectorDofs &dofs = row_data.getFieldDofs();
611 VectorDofs::iterator dit = dofs.begin();
612 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
613 if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
614 dAta.forcesOnlyOnEntitiesRow.end()) {
615 iNdices[ii] = -1;
616 }
617 }
618 }
619 CHKERR VecSetValues(getFEMethod()->snes_f, nb_dofs, indices_ptr, &nf[0],
620 ADD_VALUES);
622}
623
625 int row_side, EntityType row_type,
626 EntitiesFieldData::EntData &row_data) {
628
629 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
630 dAta.tEts.end()) {
632 }
633
634 const int nb_dofs = row_data.getIndices().size();
635 if (nb_dofs == 0)
637 if ((unsigned int)nb_dofs > 3 * row_data.getN().size2()) {
638 SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
639 }
640 const int nb_base_functions = row_data.getN().size2();
641 const int nb_gauss_pts = row_data.getN().size1();
642
643 nf.resize(nb_dofs, false);
644 nf.clear();
645
646 auto diff_base_functions = row_data.getFTensor1DiffN<3>();
647 FTensor::Index<'i', 3> i;
648 FTensor::Index<'j', 3> j;
649
650 for (int gg = 0; gg != nb_gauss_pts; gg++) {
651 double val = getVolume() * getGaussPts()(3, gg);
652 MatrixDouble3by3 &stress = commonData.sTress[gg];
654 &stress(0, 0), &stress(0, 1), &stress(0, 2), &stress(1, 0),
655 &stress(1, 1), &stress(1, 2), &stress(2, 0), &stress(2, 1),
656 &stress(2, 2));
657 FTensor::Tensor1<double *, 3> rhs(&nf[0], &nf[1], &nf[2], 3);
658 int bb = 0;
659 for (; bb != nb_dofs / 3; bb++) {
660 rhs(i) += val * t3(i, j) * diff_base_functions(j);
661 ++rhs;
662 ++diff_base_functions;
663 }
664 for (; bb != nb_base_functions; bb++) {
665 ++diff_base_functions;
666 }
667 }
668
669 CHKERR aSemble(row_side, row_type, row_data);
670
672}
673
675 BlockData &data,
676 CommonData &common_data,
677 SmartPetscObj<Vec> ghost_vec,
678 bool field_disp)
681 dAta(data), commonData(common_data), ghostVec(ghost_vec, true),
682 fieldDisp(field_disp) {}
683
685 int row_side, EntityType row_type,
686 EntitiesFieldData::EntData &row_data) {
688
689 if (row_type != MBVERTEX)
691 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
692 dAta.tEts.end()) {
694 }
695
696 std::vector<MatrixDouble> &F =
698 dAta.materialDoublePtr->F.resize(3, 3, false);
699
700 double energy = 0;
701
702 for (unsigned int gg = 0; gg != row_data.getN().size1(); ++gg) {
703 double val = getVolume() * getGaussPts()(3, gg);
704 noalias(dAta.materialDoublePtr->F) = F[gg];
705 if (fieldDisp) {
706 for (int dd = 0; dd < 3; dd++) {
707 dAta.materialDoublePtr->F(dd, dd) += 1;
708 }
709 }
710 int nb_active_variables = 0;
711 CHKERR dAta.materialDoublePtr->setUserActiveVariables(nb_active_variables);
712 CHKERR dAta.materialDoublePtr->calculateElasticEnergy(
713 dAta, getNumeredEntFiniteElementPtr());
714 energy += val * dAta.materialDoublePtr->eNergy;
715 }
716
717 CHKERR VecSetValue(ghostVec, 0, energy, ADD_VALUES);
719}
720
722 const std::string vel_field, const std::string field_name, BlockData &data,
723 CommonData &common_data)
725 vel_field, field_name, UserDataOperator::OPROWCOL),
726 dAta(data), commonData(common_data), aLe(false) {}
727
728template <int S>
730 int gg, MatrixDouble &jac_stress,
731 MatrixDouble &jac) {
733 jac.clear();
734 FTensor::Index<'i', 3> i;
735 FTensor::Index<'j', 3> j;
736 FTensor::Index<'k', 3> k;
737 int nb_col = col_data.getFieldData().size();
738 double *diff_ptr =
739 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
740 // First two indices 'i','j' derivatives of 1st Piola-stress, third index 'k'
741 // is displacement component
743 &jac_stress(3 * 0 + 0, S + 0), &jac_stress(3 * 0 + 0, S + 1),
744 &jac_stress(3 * 0 + 0, S + 2), &jac_stress(3 * 0 + 1, S + 0),
745 &jac_stress(3 * 0 + 1, S + 1), &jac_stress(3 * 0 + 1, S + 2),
746 &jac_stress(3 * 0 + 2, S + 0), &jac_stress(3 * 0 + 2, S + 1),
747 &jac_stress(3 * 0 + 2, S + 2), &jac_stress(3 * 1 + 0, S + 0),
748 &jac_stress(3 * 1 + 0, S + 1), &jac_stress(3 * 1 + 0, S + 2),
749 &jac_stress(3 * 1 + 1, S + 0), &jac_stress(3 * 1 + 1, S + 1),
750 &jac_stress(3 * 1 + 1, S + 2), &jac_stress(3 * 1 + 2, S + 0),
751 &jac_stress(3 * 1 + 2, S + 1), &jac_stress(3 * 1 + 2, S + 2),
752 &jac_stress(3 * 2 + 0, S + 0), &jac_stress(3 * 2 + 0, S + 1),
753 &jac_stress(3 * 2 + 0, S + 2), &jac_stress(3 * 2 + 1, S + 0),
754 &jac_stress(3 * 2 + 1, S + 1), &jac_stress(3 * 2 + 1, S + 2),
755 &jac_stress(3 * 2 + 2, S + 0), &jac_stress(3 * 2 + 2, S + 1),
756 &jac_stress(3 * 2 + 2, S + 2));
758 &jac_stress(3 * 0 + 0, S + 3), &jac_stress(3 * 0 + 0, S + 4),
759 &jac_stress(3 * 0 + 0, S + 5), &jac_stress(3 * 0 + 1, S + 3),
760 &jac_stress(3 * 0 + 1, S + 4), &jac_stress(3 * 0 + 1, S + 5),
761 &jac_stress(3 * 0 + 2, S + 3), &jac_stress(3 * 0 + 2, S + 4),
762 &jac_stress(3 * 0 + 2, S + 5), &jac_stress(3 * 1 + 0, S + 3),
763 &jac_stress(3 * 1 + 0, S + 4), &jac_stress(3 * 1 + 0, S + 5),
764 &jac_stress(3 * 1 + 1, S + 3), &jac_stress(3 * 1 + 1, S + 4),
765 &jac_stress(3 * 1 + 1, S + 5), &jac_stress(3 * 1 + 2, S + 3),
766 &jac_stress(3 * 1 + 2, S + 4), &jac_stress(3 * 1 + 2, S + 5),
767 &jac_stress(3 * 2 + 0, S + 3), &jac_stress(3 * 2 + 0, S + 4),
768 &jac_stress(3 * 2 + 0, S + 5), &jac_stress(3 * 2 + 1, S + 3),
769 &jac_stress(3 * 2 + 1, S + 4), &jac_stress(3 * 2 + 1, S + 5),
770 &jac_stress(3 * 2 + 2, S + 3), &jac_stress(3 * 2 + 2, S + 4),
771 &jac_stress(3 * 2 + 2, S + 5));
773 &jac_stress(3 * 0 + 0, S + 6), &jac_stress(3 * 0 + 0, S + 7),
774 &jac_stress(3 * 0 + 0, S + 8), &jac_stress(3 * 0 + 1, S + 6),
775 &jac_stress(3 * 0 + 1, S + 7), &jac_stress(3 * 0 + 1, S + 8),
776 &jac_stress(3 * 0 + 2, S + 6), &jac_stress(3 * 0 + 2, S + 7),
777 &jac_stress(3 * 0 + 2, S + 8), &jac_stress(3 * 1 + 0, S + 6),
778 &jac_stress(3 * 1 + 0, S + 7), &jac_stress(3 * 1 + 0, S + 8),
779 &jac_stress(3 * 1 + 1, S + 6), &jac_stress(3 * 1 + 1, S + 7),
780 &jac_stress(3 * 1 + 1, S + 8), &jac_stress(3 * 1 + 2, S + 6),
781 &jac_stress(3 * 1 + 2, S + 7), &jac_stress(3 * 1 + 2, S + 8),
782 &jac_stress(3 * 2 + 0, S + 6), &jac_stress(3 * 2 + 0, S + 7),
783 &jac_stress(3 * 2 + 0, S + 8), &jac_stress(3 * 2 + 1, S + 6),
784 &jac_stress(3 * 2 + 1, S + 7), &jac_stress(3 * 2 + 1, S + 8),
785 &jac_stress(3 * 2 + 2, S + 6), &jac_stress(3 * 2 + 2, S + 7),
786 &jac_stress(3 * 2 + 2, S + 8));
787 // Derivate of 1st Piola-stress multiplied by gradient of defamation for
788 // base function (dd) and displacement component (rr)
790 &jac(0, 0), &jac(1, 0), &jac(2, 0), &jac(3, 0), &jac(4, 0), &jac(5, 0),
791 &jac(6, 0), &jac(7, 0), &jac(8, 0));
793 &jac(0, 1), &jac(1, 1), &jac(2, 1), &jac(3, 1), &jac(4, 1), &jac(5, 1),
794 &jac(6, 1), &jac(7, 1), &jac(8, 1));
796 &jac(0, 2), &jac(1, 2), &jac(2, 2), &jac(3, 2), &jac(4, 2), &jac(5, 2),
797 &jac(6, 2), &jac(7, 2), &jac(8, 2));
799 diff_ptr, &diff_ptr[1], &diff_ptr[2]);
800 for (int dd = 0; dd != nb_col / 3; ++dd) {
801 t2_1_0(i, j) += t3_1_0(i, j, k) * diff(k);
802 t2_1_1(i, j) += t3_1_1(i, j, k) * diff(k);
803 t2_1_2(i, j) += t3_1_2(i, j, k) * diff(k);
804 ++t2_1_0;
805 ++t2_1_1;
806 ++t2_1_2;
807 ++diff;
808 }
810}
811
816
818 int row_side, int col_side, EntityType row_type, EntityType col_type,
820 EntitiesFieldData::EntData &col_data) {
822
823 int nb_row = row_data.getIndices().size();
824 int nb_col = col_data.getIndices().size();
825
826 int *row_indices_ptr = &row_data.getIndices()[0];
827 int *col_indices_ptr = &col_data.getIndices()[0];
828
829 if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
830 rowIndices.resize(nb_row, false);
831 noalias(rowIndices) = row_data.getIndices();
832 row_indices_ptr = &rowIndices[0];
833 VectorDofs &dofs = row_data.getFieldDofs();
834 VectorDofs::iterator dit = dofs.begin();
835 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
836 if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
837 dAta.forcesOnlyOnEntitiesRow.end()) {
838 rowIndices[ii] = -1;
839 }
840 }
841 }
842
843 if (!dAta.forcesOnlyOnEntitiesCol.empty()) {
844 colIndices.resize(nb_col, false);
845 noalias(colIndices) = col_data.getIndices();
846 col_indices_ptr = &colIndices[0];
847 VectorDofs &dofs = col_data.getFieldDofs();
848 VectorDofs::iterator dit = dofs.begin();
849 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
850 if (dAta.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
851 dAta.forcesOnlyOnEntitiesCol.end()) {
852 colIndices[ii] = -1;
853 }
854 }
855 }
856
857 CHKERR MatSetValues(getFEMethod()->snes_B, nb_row, row_indices_ptr, nb_col,
858 col_indices_ptr, &k(0, 0), ADD_VALUES);
859
860 // is symmetric
861 if (row_side != col_side || row_type != col_type) {
862
863 row_indices_ptr = &row_data.getIndices()[0];
864 col_indices_ptr = &col_data.getIndices()[0];
865
866 if (!dAta.forcesOnlyOnEntitiesCol.empty()) {
867 rowIndices.resize(nb_row, false);
868 noalias(rowIndices) = row_data.getIndices();
869 row_indices_ptr = &rowIndices[0];
870 VectorDofs &dofs = row_data.getFieldDofs();
871 VectorDofs::iterator dit = dofs.begin();
872 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
873 if (dAta.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
874 dAta.forcesOnlyOnEntitiesCol.end()) {
875 rowIndices[ii] = -1;
876 }
877 }
878 }
879
880 if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
881 colIndices.resize(nb_col, false);
882 noalias(colIndices) = col_data.getIndices();
883 col_indices_ptr = &colIndices[0];
884 VectorDofs &dofs = col_data.getFieldDofs();
885 VectorDofs::iterator dit = dofs.begin();
886 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
887 if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
888 dAta.forcesOnlyOnEntitiesRow.end()) {
889 colIndices[ii] = -1;
890 }
891 }
892 }
893
894 trans_k.resize(nb_col, nb_row, false);
895 noalias(trans_k) = trans(k);
896 CHKERR MatSetValues(getFEMethod()->snes_B, nb_col, col_indices_ptr, nb_row,
897 row_indices_ptr, &trans_k(0, 0), ADD_VALUES);
898 }
899
901}
902
904 int row_side, int col_side, EntityType row_type, EntityType col_type,
906 EntitiesFieldData::EntData &col_data) {
908
909 int nb_row = row_data.getIndices().size();
910 int nb_col = col_data.getIndices().size();
911 if (nb_row == 0)
913 if (nb_col == 0)
915
916 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
917 dAta.tEts.end()) {
919 }
920
921 // const int nb_base_functions = row_data.getN().size2();
922 const int nb_gauss_pts = row_data.getN().size1();
923
924 FTensor::Index<'i', 3> i;
925 FTensor::Index<'j', 3> j;
926 FTensor::Index<'m', 3> m;
927
928 k.resize(nb_row, nb_col, false);
929 k.clear();
930 jac.resize(9, nb_col, false);
931
932 for (int gg = 0; gg != nb_gauss_pts; gg++) {
933 CHKERR getJac(col_data, gg);
934 double val = getVolume() * getGaussPts()(3, gg);
936 &jac(3 * 0 + 0, 0), &jac(3 * 0 + 0, 1), &jac(3 * 0 + 0, 2),
937 &jac(3 * 0 + 1, 0), &jac(3 * 0 + 1, 1), &jac(3 * 0 + 1, 2),
938 &jac(3 * 0 + 2, 0), &jac(3 * 0 + 2, 1), &jac(3 * 0 + 2, 2),
939 &jac(3 * 1 + 0, 0), &jac(3 * 1 + 0, 1), &jac(3 * 1 + 0, 2),
940 &jac(3 * 1 + 1, 0), &jac(3 * 1 + 1, 1), &jac(3 * 1 + 1, 2),
941 &jac(3 * 1 + 2, 0), &jac(3 * 1 + 2, 1), &jac(3 * 1 + 2, 2),
942 &jac(3 * 2 + 0, 0), &jac(3 * 2 + 0, 1), &jac(3 * 2 + 0, 2),
943 &jac(3 * 2 + 1, 0), &jac(3 * 2 + 1, 1), &jac(3 * 2 + 1, 2),
944 &jac(3 * 2 + 2, 0), &jac(3 * 2 + 2, 1), &jac(3 * 2 + 2, 2));
945 for (int cc = 0; cc != nb_col / 3; cc++) {
946 auto diff_base_functions = row_data.getFTensor1DiffN<3>(gg, 0);
948 &k(0, 3 * cc + 0), &k(0, 3 * cc + 1), &k(0, 3 * cc + 2),
949 &k(1, 3 * cc + 0), &k(1, 3 * cc + 1), &k(1, 3 * cc + 2),
950 &k(2, 3 * cc + 0), &k(2, 3 * cc + 1), &k(2, 3 * cc + 2), 3 * nb_col);
951 for (int rr = 0; rr != nb_row / 3; rr++) {
952 lhs(i, j) += val * t3_1(i, m, j) * diff_base_functions(m);
953 ++diff_base_functions;
954 ++lhs;
955 }
956 ++t3_1;
957 }
958 }
959
960 CHKERR aSemble(row_side, col_side, row_type, col_type, row_data, col_data);
961
963}
964
966 const std::string vel_field, const std::string field_name, BlockData &data,
967 CommonData &common_data)
968 : OpLhsPiolaKirchhoff_dx(vel_field, field_name, data, common_data) {
969 sYmm = false;
970}
971
976
978 int row_side, int col_side, EntityType row_type, EntityType col_type,
980 EntitiesFieldData::EntData &col_data) {
982
983 int nb_row = row_data.getIndices().size();
984 int nb_col = col_data.getIndices().size();
985
986 int *row_indices_ptr = &row_data.getIndices()[0];
987 if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
988 rowIndices.resize(nb_row, false);
989 noalias(rowIndices) = row_data.getIndices();
990 row_indices_ptr = &rowIndices[0];
991 VectorDofs &dofs = row_data.getFieldDofs();
992 VectorDofs::iterator dit = dofs.begin();
993 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
994 if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
995 dAta.forcesOnlyOnEntitiesRow.end()) {
996 rowIndices[ii] = -1;
997 }
998 }
999 }
1000
1001 int *col_indices_ptr = &col_data.getIndices()[0];
1002 if (!dAta.forcesOnlyOnEntitiesCol.empty()) {
1003 colIndices.resize(nb_col, false);
1004 noalias(colIndices) = col_data.getIndices();
1005 col_indices_ptr = &colIndices[0];
1006 VectorDofs &dofs = col_data.getFieldDofs();
1007 VectorDofs::iterator dit = dofs.begin();
1008 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
1009 if (dAta.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
1010 dAta.forcesOnlyOnEntitiesCol.end()) {
1011 colIndices[ii] = -1;
1012 }
1013 }
1014 }
1015
1016 /*for(int dd1 = 0;dd1<k.size1();dd1++) {
1017 for(int dd2 = 0;dd2<k.size2();dd2++) {
1018 if(k(dd1,dd2)!=k(dd1,dd2)) {
1019 SETERRQ(PETSC_COMM_SELF,1,"Wrong result");
1020 }
1021 }
1022 }*/
1023
1024 CHKERR MatSetValues(getFEMethod()->snes_B, nb_row, row_indices_ptr, nb_col,
1025 col_indices_ptr, &k(0, 0), ADD_VALUES);
1026
1028}
1029
1031 const std::string field_name, BlockData &data, CommonData &common_data,
1032 int tag, bool jacobian, bool ale)
1033 : OpJacobianPiolaKirchhoffStress(field_name, data, common_data, tag,
1034 jacobian, ale, false) {}
1035
1038 const int gg) {
1040
1041 CHKERR dAta.materialAdoublePtr->calculatesIGma_EshelbyStress(
1042 dAta, getNumeredEntFiniteElementPtr());
1043 if (aLe) {
1044 auto &t_sIGma = dAta.materialAdoublePtr->t_sIGma;
1045 auto &t_invH = dAta.materialAdoublePtr->t_invH;
1046 t_sIGma(i, j) = t_sIGma(i, k) * t_invH(j, k);
1047 t_sIGma(i, j) *= dAta.materialAdoublePtr->detH;
1048
1049 }
1050 commonData.sTress[gg].resize(3, 3, false);
1051 for (int dd1 = 0; dd1 < 3; dd1++) {
1052 for (int dd2 = 0; dd2 < 3; dd2++) {
1053 dAta.materialAdoublePtr->sIGma(dd1, dd2) >>=
1054 (commonData.sTress[gg])(dd1, dd2);
1055 }
1056 }
1057
1059}
1060
1064
1066 const std::string vel_field, const std::string field_name, BlockData &data,
1067 CommonData &common_data)
1068 : OpLhsPiolaKirchhoff_dX(vel_field, field_name, data, common_data) {}
1069
1074
1076 const std::string vel_field, const std::string field_name, BlockData &data,
1077 CommonData &common_data)
1078 : OpLhsPiolaKirchhoff_dx(vel_field, field_name, data, common_data) {}
1079
1084
1087 materialDoublePtr,
1089 materialAdoublePtr) {
1091
1092 if (!materialDoublePtr) {
1094 "Pointer for materialDoublePtr not allocated");
1095 }
1096 if (!materialAdoublePtr) {
1098 "Pointer for materialAdoublePtr not allocated");
1099 }
1100
1102 mField, BLOCKSET | MAT_ELASTICSET, it)) {
1103 Mat_Elastic mydata;
1104 CHKERR it->getAttributeDataStructure(mydata);
1105 int id = it->getMeshsetId();
1106 EntityHandle meshset = it->getMeshset();
1107 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET,
1108 setOfBlocks[id].tEts, true);
1109 setOfBlocks[id].iD = id;
1110 setOfBlocks[id].E = mydata.data.Young;
1111 setOfBlocks[id].PoissonRatio = mydata.data.Poisson;
1112 setOfBlocks[id].materialDoublePtr = materialDoublePtr;
1113 setOfBlocks[id].materialAdoublePtr = materialAdoublePtr;
1114 }
1115
1117}
1118
1120 const std::string element_name,
1121 const std::string spatial_position_field_name,
1122 const std::string material_position_field_name, const bool ale) {
1124
1125 CHKERR mField.add_finite_element(element_name, MF_ZERO);
1127 element_name, spatial_position_field_name);
1129 element_name, spatial_position_field_name);
1131 element_name, spatial_position_field_name);
1132 if (mField.check_field(material_position_field_name)) {
1133 if (ale) {
1135 element_name, material_position_field_name);
1137 element_name, material_position_field_name);
1138 }
1140 element_name, material_position_field_name);
1141 }
1142
1143 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
1144 for (; sit != setOfBlocks.end(); sit++) {
1145 CHKERR mField.add_ents_to_finite_element_by_type(sit->second.tEts, MBTET,
1146 element_name);
1147 }
1148
1150}
1151
1153 const std::string spatial_position_field_name,
1154 const std::string material_position_field_name, const bool ale,
1155 const bool field_disp) {
1157
1158 commonData.spatialPositions = spatial_position_field_name;
1159 commonData.meshPositions = material_position_field_name;
1160
1161 // Rhs
1162 feRhs.getOpPtrVector().push_back(
1163 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
1164 if (mField.check_field(material_position_field_name)) {
1166 material_position_field_name, commonData));
1167 }
1168 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
1169 for (; sit != setOfBlocks.end(); sit++) {
1171 spatial_position_field_name, sit->second, commonData, tAg, false, ale,
1172 field_disp));
1174 spatial_position_field_name, sit->second, commonData));
1175 }
1176
1177 // Energy
1178 feEnergy.getOpPtrVector().push_back(
1179 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
1180 if (mField.check_field(material_position_field_name)) {
1182 material_position_field_name, commonData));
1183 }
1184 sit = setOfBlocks.begin();
1185 for (; sit != setOfBlocks.end(); sit++) {
1186 feEnergy.getOpPtrVector().push_back(
1187 new OpEnergy(spatial_position_field_name, sit->second, commonData,
1188 feEnergy.V, field_disp));
1189 }
1190
1191 // Lhs
1192 feLhs.getOpPtrVector().push_back(
1193 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
1194 if (mField.check_field(material_position_field_name)) {
1196 material_position_field_name, commonData));
1197 }
1198 sit = setOfBlocks.begin();
1199 for (; sit != setOfBlocks.end(); sit++) {
1201 spatial_position_field_name, sit->second, commonData, tAg, true, ale,
1202 field_disp));
1204 spatial_position_field_name, spatial_position_field_name, sit->second,
1205 commonData));
1206 }
1207
1209}
static MoFEMErrorCode calculateEnergy(DM dm, boost::shared_ptr< map< int, BlockData > > block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, SmartPetscObj< Vec > &v_energy_ptr)
ForcesAndSourcesCore::UserDataOperator UserDataOperator
static MoFEMErrorCode get_jac(EntitiesFieldData::EntData &col_data, int gg, MatrixDouble &jac_stress, MatrixDouble &jac)
Operators and data structures for non-linear elastic analysis.
@ MF_ZERO
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
@ BLOCKSET
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
constexpr int order
@ F
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
virtual bool check_field(const std::string &name) const =0
check if field is in database
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
const double n
refractive index of diffusive medium
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
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 Common.hpp:10
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
constexpr AssemblyType A
double H
Hardening.
Definition plastic.cpp:128
constexpr auto field_name
FTensor::Index< 'm', 3 > m
const int N
Definition speed_test.cpp:3
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from field data coefficients.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
const VectorInt & getIndices() const
Get global indices of dofs on entity.
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Elastic material data structure.
intrusive_ptr for managing petsc objects
FTensor::Tensor2< double *, 3, 3 > & getJac()
get element Jacobian
data for calculation heat conductivity and heat capacity elements
common data used by volume elements
std::vector< MatrixDouble > jacStress
this is simply material tangent operator
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
Implementation of elastic (non-linear) St. Kirchhoff equation.
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode postProcess()
function is run at the end of loop
int getRule(int order)
it is used to calculate nb. of Gauss integration points
OpEnergy(const std::string field_name, BlockData &data, CommonData &common_data, SmartPetscObj< Vec > ghost_vec, bool field_disp)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpGetCommonDataAtGaussPts(const std::string field_name, CommonData &common_data)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
operator calculating deformation gradient
OpGetDataAtGaussPts(const std::string field_name, std::vector< VectorDouble > &values_at_gauss_pts, std::vector< MatrixDouble > &gradient_at_gauss_pts)
virtual MoFEMErrorCode playTag(const int gg)
Play ADOL-C tape.
OpJacobianEnergy(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool gradient, bool hessian, bool ale, bool field_disp)
virtual MoFEMErrorCode recordTag(const int gg)
Record ADOL-C tape.
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
virtual MoFEMErrorCode calculateEnergy(const int gg)
Calculate Paola-Kirchhoff I stress.
MoFEMErrorCode calculateStress(const int gg)
Calculate Paola-Kirchhoff I stress.
OpJacobianEshelbyStress(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian, bool ale)
virtual MoFEMErrorCode recordTag(const int gg)
Record ADOL-C tape.
virtual MoFEMErrorCode playTag(const int gg)
Play ADOL-C tape.
OpJacobianPiolaKirchhoffStress(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian, bool ale, bool field_disp)
Construct operator to calculate Piola-Kirchhoff stress or its derivatives over gradient deformation.
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Calculate stress or jacobian at gauss points.
virtual MoFEMErrorCode calculateStress(const int gg)
Calculate Paola-Kirchhoff I stress.
OpLhsEshelby_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
OpLhsEshelby_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
OpLhsPiolaKirchhoff_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
virtual MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpLhsPiolaKirchhoff_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
OpRhsEshelbyStress(const std::string field_name, BlockData &data, CommonData &common_data)
virtual MoFEMErrorCode aSemble(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpRhsPiolaKirchhoff(const std::string field_name, BlockData &data, CommonData &common_data)
MyVolumeFE feRhs
calculate right hand side for tetrahedral elements
NonlinearElasticElement(MoFEM::Interface &m_field, short int tag)
MoFEMErrorCode addElement(const std::string element_name, const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false)
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
MoFEMErrorCode setBlocks(boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< double > > materialDoublePtr, boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< adouble > > materialAdoublePtr)
MyVolumeFE feEnergy
calculate elastic energy
MoFEMErrorCode setOperators(const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false, const bool field_disp=false)
Set operators to calculate left hand tangent matrix and right hand residual.