v0.13.2
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_NULL), F(PETSC_NULL),
19 addToRule(1) {
20
21 auto create_vec = [&]() {
22 constexpr int ghosts[] = {0};
23 if (mField.get_comm_rank() == 0) {
24 return createSmartVectorMPI(mField.get_comm(), 1, 1);
25 } else {
26 return createSmartVectorMPI(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_NULL) {
43 snes_B = A;
44 }
45
46 if (F != PETSC_NULL) {
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 const double *array;
65
66 switch (snes_ctx) {
67 case CTX_SNESNONE:
68 CHKERR VecAssemblyBegin(V);
69 CHKERR VecAssemblyEnd(V);
70 CHKERR VecSum(V, &eNergy);
71 break;
72 default:
73 break;
74 }
75
76 CHKERR VolumeElementForcesAndSourcesCore::postProcess();
77
79}
80
82 short int tag)
83 : feRhs(m_field), feLhs(m_field), feEnergy(m_field), mField(m_field),
84 tAg(tag) {}
85
87 const std::string field_name,
88 std::vector<VectorDouble> &values_at_gauss_pts,
89 std::vector<MatrixDouble> &gardient_at_gauss_pts)
92 valuesAtGaussPts(values_at_gauss_pts),
93 gradientAtGaussPts(gardient_at_gauss_pts), zeroAtType(MBVERTEX) {}
94
96 int side, EntityType type, EntitiesFieldData::EntData &data) {
98
99 const int nb_dofs = data.getFieldData().size();
100 const int nb_base_functions = data.getN().size2();
101 if (nb_dofs == 0) {
103 }
104 const int nb_gauss_pts = data.getN().size1();
105 const int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
106
107 // initialize
108 if (type == zeroAtType) {
109 valuesAtGaussPts.resize(nb_gauss_pts);
110 gradientAtGaussPts.resize(nb_gauss_pts);
111 for (int gg = 0; gg != nb_gauss_pts; gg++) {
112 valuesAtGaussPts[gg].resize(rank, false);
113 gradientAtGaussPts[gg].resize(rank, 3, false);
114 }
115 for (int gg = 0; gg != nb_gauss_pts; gg++) {
116 valuesAtGaussPts[gg].clear();
117 gradientAtGaussPts[gg].clear();
118 }
119 }
120
121 auto base_function = data.getFTensor0N();
122 auto diff_base_functions = data.getFTensor1DiffN<3>();
125
126 if (rank == 1) {
127
128 for (int gg = 0; gg != nb_gauss_pts; gg++) {
129 auto field_data = data.getFTensor0FieldData();
130 double &val = valuesAtGaussPts[gg][0];
131 FTensor::Tensor1<double *, 3> grad(&gradientAtGaussPts[gg](0, 0),
132 &gradientAtGaussPts[gg](0, 1),
133 &gradientAtGaussPts[gg](0, 2));
134 int bb = 0;
135 for (; bb != nb_dofs; bb++) {
136 val += base_function * field_data;
137 grad(i) += diff_base_functions(i) * field_data;
138 ++diff_base_functions;
139 ++base_function;
140 ++field_data;
141 }
142 for (; bb != nb_base_functions; bb++) {
143 ++diff_base_functions;
144 ++base_function;
145 }
146 }
147
148 } else if (rank == 3) {
149
150 for (int gg = 0; gg != nb_gauss_pts; gg++) {
151 auto field_data = data.getFTensor1FieldData<3>();
152 FTensor::Tensor1<double *, 3> values(&valuesAtGaussPts[gg][0],
153 &valuesAtGaussPts[gg][1],
154 &valuesAtGaussPts[gg][2]);
156 &gradientAtGaussPts[gg](0, 0), &gradientAtGaussPts[gg](0, 1),
157 &gradientAtGaussPts[gg](0, 2), &gradientAtGaussPts[gg](1, 0),
158 &gradientAtGaussPts[gg](1, 1), &gradientAtGaussPts[gg](1, 2),
159 &gradientAtGaussPts[gg](2, 0), &gradientAtGaussPts[gg](2, 1),
160 &gradientAtGaussPts[gg](2, 2));
161 int bb = 0;
162 for (; bb != nb_dofs / 3; bb++) {
163 values(i) += base_function * field_data(i);
164 gradient(i, j) += field_data(i) * diff_base_functions(j);
165 ++diff_base_functions;
166 ++base_function;
167 ++field_data;
168 }
169 for (; bb != nb_base_functions; bb++) {
170 ++diff_base_functions;
171 ++base_function;
172 }
173 }
174
175 } else {
176 // FIXME: THat part is inefficient
177 VectorDouble &values = data.getFieldData();
178 for (int gg = 0; gg < nb_gauss_pts; gg++) {
179 VectorAdaptor N = data.getN(gg, nb_dofs / rank);
180 MatrixAdaptor diffN = data.getDiffN(gg, nb_dofs / rank);
181 for (int dd = 0; dd < nb_dofs / rank; dd++) {
182 for (int rr1 = 0; rr1 < rank; rr1++) {
183 valuesAtGaussPts[gg][rr1] += N[dd] * values[rank * dd + rr1];
184 for (int rr2 = 0; rr2 < 3; rr2++) {
185 gradientAtGaussPts[gg](rr1, rr2) +=
186 diffN(dd, rr2) * values[rank * dd + rr1];
187 }
188 }
189 }
190 }
191 }
192
194}
195
197 const std::string field_name, CommonData &common_data)
198 : OpGetDataAtGaussPts(field_name, common_data.dataAtGaussPts[field_name],
199 common_data.gradAtGaussPts[field_name]) {}
200
203 BlockData &data, CommonData &common_data,
204 int tag, bool jacobian, bool ale,
205 bool field_disp)
208 dAta(data), commonData(common_data), tAg(tag), adlocReturnValue(0),
209 jAcobian(jacobian), fUnction(!jacobian), aLe(ale), fieldDisp(field_disp) {
210
211}
212
215 const int gg) {
217
218 CHKERR dAta.materialAdoublePtr->calculateP_PiolaKirchhoffI(
219 dAta, getNumeredEntFiniteElementPtr());
220
221 if (aLe) {
222 auto &t_P = dAta.materialAdoublePtr->t_P;
223 auto &t_invH = dAta.materialAdoublePtr->t_invH;
224 t_P(i, j) = t_P(i, k) * t_invH(j, k);
225 t_P(i, j) *= dAta.materialAdoublePtr->detH;
226 }
227
228 commonData.sTress[gg].resize(3, 3, false);
229 for (int dd1 = 0; dd1 < 3; dd1++) {
230 for (int dd2 = 0; dd2 < 3; dd2++) {
231 dAta.materialAdoublePtr->P(dd1, dd2) >>=
232 (commonData.sTress[gg])(dd1, dd2);
233 }
234 }
235
237}
238
241 const int gg) {
243
244 trace_on(tAg, 0);
245
246 dAta.materialAdoublePtr->F.resize(3, 3, false);
247
248 if (!aLe) {
249
250 nbActiveVariables = 0;
251 for (int dd1 = 0; dd1 < 3; dd1++) {
252 for (int dd2 = 0; dd2 < 3; dd2++) {
253 dAta.materialAdoublePtr->F(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
254 if (fieldDisp) {
255 if (dd1 == dd2) {
256 dAta.materialAdoublePtr->F(dd1, dd2) += 1;
257 }
258 }
259 nbActiveVariables++;
260 }
261 }
262
263 } else {
264
265 nbActiveVariables = 0;
266
267 dAta.materialAdoublePtr->h.resize(3, 3, false);
268 for (int dd1 = 0; dd1 < 3; dd1++) {
269 for (int dd2 = 0; dd2 < 3; dd2++) {
270 dAta.materialAdoublePtr->h(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
271 nbActiveVariables++;
272 }
273 }
274
275 dAta.materialAdoublePtr->H.resize(3, 3, false);
276 for (int dd1 = 0; dd1 < 3; dd1++) {
277 for (int dd2 = 0; dd2 < 3; dd2++) {
278 dAta.materialAdoublePtr->H(dd1, dd2) <<= (*ptrH)[gg](dd1, dd2);
279 nbActiveVariables++;
280 }
281 }
282
283 dAta.materialAdoublePtr->detH = determinantTensor3by3(dAta.materialAdoublePtr->H);
284 dAta.materialAdoublePtr->invH.resize(3, 3, false);
285 CHKERR invertTensor3by3(dAta.materialAdoublePtr->H,
286 dAta.materialAdoublePtr->detH,
287 dAta.materialAdoublePtr->invH);
288
289 auto &t_F = dAta.materialAdoublePtr->t_F;
290 auto &t_h = dAta.materialAdoublePtr->t_h;
291 auto &t_invH = dAta.materialAdoublePtr->t_invH;
292
293 t_F(i, j) = t_h(i, k) * t_invH(k, j);
294
295 }
296
297 CHKERR dAta.materialAdoublePtr->setUserActiveVariables(nbActiveVariables);
298 CHKERR calculateStress(gg);
299
300 trace_off();
301
303}
304
308
309 int r;
310
311 if (fUnction) {
312 commonData.sTress[gg].resize(3, 3, false);
313 // play recorder for values
314 r = ::function(tAg, 9, nbActiveVariables, &activeVariables[0],
315 &commonData.sTress[gg](0, 0));
316 if (r < adlocReturnValue) { // function is locally analytic
317 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
318 "ADOL-C function evaluation with error r = %d", r);
319 }
320 }
321
322 if (jAcobian) {
323 commonData.jacStress[gg].resize(9, nbActiveVariables, false);
324 double *jac_ptr[] = {
325 &(commonData.jacStress[gg](0, 0)), &(commonData.jacStress[gg](1, 0)),
326 &(commonData.jacStress[gg](2, 0)), &(commonData.jacStress[gg](3, 0)),
327 &(commonData.jacStress[gg](4, 0)), &(commonData.jacStress[gg](5, 0)),
328 &(commonData.jacStress[gg](6, 0)), &(commonData.jacStress[gg](7, 0)),
329 &(commonData.jacStress[gg](8, 0))};
330 // play recorder for jacobians
331 r = jacobian(tAg, 9, nbActiveVariables, &activeVariables[0], jac_ptr);
332 if (r < adlocReturnValue) {
333 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
334 "ADOL-C function evaluation with error");
335 }
336 }
337
339}
340
342 int row_side, EntityType row_type,
343 EntitiesFieldData::EntData &row_data) {
345
346 // do it only once, no need to repeat this for edges,faces or tets
347 if (row_type != MBVERTEX)
349
350 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
351 dAta.tEts.end()) {
353 }
354
355 int nb_dofs = row_data.getFieldData().size();
356 if (nb_dofs == 0)
358 dAta.materialAdoublePtr->commonDataPtr = &commonData;
359 dAta.materialAdoublePtr->opPtr = this;
360
361 int nb_gauss_pts = row_data.getN().size1();
362 commonData.sTress.resize(nb_gauss_pts);
363 commonData.jacStress.resize(nb_gauss_pts);
364
366 if (aLe) {
368 }
369
370 for (int gg = 0; gg != nb_gauss_pts; gg++) {
371
372 dAta.materialAdoublePtr->gG = gg;
373
374 // Record tag and calculate stress
375 if (recordTagForIntegrationPoint(gg)) {
376 CHKERR recordTag(gg);
377 }
378
379 // Set active variables vector
380 if (jAcobian || (!recordTagForIntegrationPoint(gg))) {
381 activeVariables.resize(nbActiveVariables, false);
382 if (!aLe) {
383 for (int dd1 = 0; dd1 < 3; dd1++) {
384 for (int dd2 = 0; dd2 < 3; dd2++) {
385 activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
386 }
387 }
388 } else {
389 for (int dd1 = 0; dd1 < 3; dd1++) {
390 for (int dd2 = 0; dd2 < 3; dd2++) {
391 activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
392 }
393 }
394 for (int dd1 = 0; dd1 < 3; dd1++) {
395 for (int dd2 = 0; dd2 < 3; dd2++) {
396 activeVariables(9 + dd1 * 3 + dd2) = (*ptrH)[gg](dd1, dd2);
397 }
398 }
399 }
400 CHKERR dAta.materialAdoublePtr->setUserActiveVariables(activeVariables);
401
402 // Play tag and calculate stress or tangent
403 if (jAcobian || (!recordTagForIntegrationPoint(gg))) {
404 CHKERR playTag(gg);
405 }
406 }
407 }
408
410}
411
413 const std::string
414 field_name, ///< field name for spatial positions or displacements
415 BlockData &data, CommonData &common_data, int tag, bool gradient,
416 bool hessian, bool ale, bool field_disp)
419 dAta(data), commonData(common_data), tAg(tag), gRadient(gradient),
420 hEssian(hessian), aLe(ale), fieldDisp(field_disp) {}
421
425 CHKERR dAta.materialAdoublePtr->calculateElasticEnergy(
426 dAta, getNumeredEntFiniteElementPtr());
427 dAta.materialAdoublePtr->eNergy >>= commonData.eNergy[gg];
429}
430
434
435 trace_on(tAg, 0);
436
437 if (!aLe) {
438
439 nbActiveVariables = 0;
440 for (int dd1 = 0; dd1 < 3; dd1++) {
441 for (int dd2 = 0; dd2 < 3; dd2++) {
442 dAta.materialAdoublePtr->F(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
443 if (fieldDisp) {
444 if (dd1 == dd2) {
445 dAta.materialAdoublePtr->F(dd1, dd2) += 1;
446 }
447 }
448 nbActiveVariables++;
449 }
450 }
451
452 } else {
453
454 nbActiveVariables = 0;
455
456 dAta.materialAdoublePtr->h.resize(3, 3, false);
457 for (int dd1 = 0; dd1 < 3; dd1++) {
458 for (int dd2 = 0; dd2 < 3; dd2++) {
459 dAta.materialAdoublePtr->h(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
460 nbActiveVariables++;
461 }
462 }
463
464 dAta.materialAdoublePtr->H.resize(3, 3, false);
465 for (int dd1 = 0; dd1 < 3; dd1++) {
466 for (int dd2 = 0; dd2 < 3; dd2++) {
467 dAta.materialAdoublePtr->H(dd1, dd2) <<= (*ptrH)[gg](dd1, dd2);
468 nbActiveVariables++;
469 }
470 }
471
472 dAta.materialAdoublePtr->detH = determinantTensor3by3(dAta.materialAdoublePtr->H);
473 dAta.materialAdoublePtr->invH.resize(3, 3, false);
474 CHKERR invertTensor3by3(dAta.materialAdoublePtr->H,
475 dAta.materialAdoublePtr->detH,
476 dAta.materialAdoublePtr->invH);
477
478 auto &t_F = dAta.materialAdoublePtr->t_F;
479 auto &t_h = dAta.materialAdoublePtr->t_h;
480 auto &t_invH = dAta.materialAdoublePtr->t_invH;
481
482 t_F(i, j) = t_h(i, k) * t_invH(k, j);
483
484 }
485
486 CHKERR dAta.materialAdoublePtr->setUserActiveVariables(nbActiveVariables);
488
489 trace_off();
490
492}
493
497
498 if (gRadient) {
499 commonData.jacEnergy[gg].resize(nbActiveVariables, false);
500 int r = ::gradient(tAg, nbActiveVariables, &activeVariables[0],
501 &commonData.jacEnergy[gg][0]);
502 if (r < 0) {
503 // That means that energy function is not smooth and derivative
504 // can not be calculated,
505 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
506 "ADOL-C function evaluation with error");
507 }
508 }
509
510 if (hEssian) {
511 commonData.hessianEnergy[gg].resize(nbActiveVariables * nbActiveVariables,
512 false);
513 double *H[nbActiveVariables];
514 for (int n = 0; n != nbActiveVariables; n++) {
515 H[n] = &(commonData.hessianEnergy[gg][n * nbActiveVariables]);
516 }
517 int r = ::hessian(tAg, nbActiveVariables, &*activeVariables.begin(), H);
518 if (r < 0) {
519 // That means that energy function is not smooth and derivative
520 // can not be calculated,
521 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
522 "ADOL-C function evaluation with error");
523 }
524 }
525
527}
528
530 int row_side, EntityType row_type,
531 EntitiesFieldData::EntData &row_data) {
533
534 // do it only once, no need to repeat this for edges,faces or tets
535 if (row_type != MBVERTEX)
537
538 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
539 dAta.tEts.end()) {
541 }
542
543 int nb_dofs = row_data.getFieldData().size();
544 if (nb_dofs == 0)
546 dAta.materialAdoublePtr->commonDataPtr = &commonData;
547 dAta.materialAdoublePtr->opPtr = this;
548
549 int nb_gauss_pts = row_data.getN().size1();
550 commonData.eNergy.resize(nb_gauss_pts);
551 commonData.jacEnergy.resize(nb_gauss_pts);
552
554 if (aLe) {
556 }
557
558 for (int gg = 0; gg != nb_gauss_pts; gg++) {
559
560 dAta.materialAdoublePtr->gG = gg;
561
562 // Record tag and calualte stress
563 if (recordTagForIntegrationPoint(gg)) {
564 CHKERR recordTag(gg);
565 }
566
567 activeVariables.resize(nbActiveVariables, false);
568 if (!aLe) {
569 for (int dd1 = 0; dd1 < 3; dd1++) {
570 for (int dd2 = 0; dd2 < 3; dd2++) {
571 activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
572 }
573 }
574 } else {
575 for (int dd1 = 0; dd1 < 3; dd1++) {
576 for (int dd2 = 0; dd2 < 3; dd2++) {
577 activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
578 }
579 }
580 for (int dd1 = 0; dd1 < 3; dd1++) {
581 for (int dd2 = 0; dd2 < 3; dd2++) {
582 activeVariables(9 + dd1 * 3 + dd2) = (*ptrH)[gg](dd1, dd2);
583 }
584 }
585 }
586 CHKERR dAta.materialAdoublePtr->setUserActiveVariables(activeVariables);
587
588 // Play tag and calculate stress or tangent
589 CHKERR playTag(gg);
590 }
591
593}
594
596 const std::string field_name, BlockData &data, CommonData &common_data)
599 dAta(data), commonData(common_data), aLe(false) {}
600
602 int row_side, EntityType row_type,
603 EntitiesFieldData::EntData &row_data) {
605
606 int nb_dofs = row_data.getIndices().size();
607 int *indices_ptr = &row_data.getIndices()[0];
608 if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
609 iNdices.resize(nb_dofs, false);
610 noalias(iNdices) = row_data.getIndices();
611 indices_ptr = &iNdices[0];
612 VectorDofs &dofs = row_data.getFieldDofs();
613 VectorDofs::iterator dit = dofs.begin();
614 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
615 if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
616 dAta.forcesOnlyOnEntitiesRow.end()) {
617 iNdices[ii] = -1;
618 }
619 }
620 }
621 CHKERR VecSetValues(getFEMethod()->snes_f, nb_dofs, indices_ptr, &nf[0],
622 ADD_VALUES);
624}
625
627 int row_side, EntityType row_type,
628 EntitiesFieldData::EntData &row_data) {
630
631 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
632 dAta.tEts.end()) {
634 }
635
636 const int nb_dofs = row_data.getIndices().size();
637 if (nb_dofs == 0)
639 if ((unsigned int)nb_dofs > 3 * row_data.getN().size2()) {
640 SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
641 }
642 const int nb_base_functions = row_data.getN().size2();
643 const int nb_gauss_pts = row_data.getN().size1();
644
645 nf.resize(nb_dofs, false);
646 nf.clear();
647
648 auto diff_base_functions = row_data.getFTensor1DiffN<3>();
651
652 for (int gg = 0; gg != nb_gauss_pts; gg++) {
653 double val = getVolume() * getGaussPts()(3, gg);
654 MatrixDouble3by3 &stress = commonData.sTress[gg];
656 &stress(0, 0), &stress(0, 1), &stress(0, 2), &stress(1, 0),
657 &stress(1, 1), &stress(1, 2), &stress(2, 0), &stress(2, 1),
658 &stress(2, 2));
659 FTensor::Tensor1<double *, 3> rhs(&nf[0], &nf[1], &nf[2], 3);
660 int bb = 0;
661 for (; bb != nb_dofs / 3; bb++) {
662 rhs(i) += val * t3(i, j) * diff_base_functions(j);
663 ++rhs;
664 ++diff_base_functions;
665 }
666 for (; bb != nb_base_functions; bb++) {
667 ++diff_base_functions;
668 }
669 }
670
671 CHKERR aSemble(row_side, row_type, row_data);
672
674}
675
677 BlockData &data,
678 CommonData &common_data,
679 SmartPetscObj<Vec> ghost_vec,
680 bool field_disp)
683 dAta(data), commonData(common_data), ghostVec(ghost_vec, true),
684 fieldDisp(field_disp) {}
685
687 int row_side, EntityType row_type,
688 EntitiesFieldData::EntData &row_data) {
690
691 if (row_type != MBVERTEX)
693 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
694 dAta.tEts.end()) {
696 }
697
698 std::vector<MatrixDouble> &F =
700 dAta.materialDoublePtr->F.resize(3, 3, false);
701
702 double energy = 0;
703
704 for (unsigned int gg = 0; gg != row_data.getN().size1(); ++gg) {
705 double val = getVolume() * getGaussPts()(3, gg);
706 noalias(dAta.materialDoublePtr->F) = F[gg];
707 if (fieldDisp) {
708 for (int dd = 0; dd < 3; dd++) {
709 dAta.materialDoublePtr->F(dd, dd) += 1;
710 }
711 }
712 int nb_active_variables = 0;
713 CHKERR dAta.materialDoublePtr->setUserActiveVariables(nb_active_variables);
714 CHKERR dAta.materialDoublePtr->calculateElasticEnergy(
715 dAta, getNumeredEntFiniteElementPtr());
716 energy += val * dAta.materialDoublePtr->eNergy;
717 }
718
719 CHKERR VecSetValue(ghostVec, 0, energy, ADD_VALUES);
721}
722
724 const std::string vel_field, const std::string field_name, BlockData &data,
725 CommonData &common_data)
727 vel_field, field_name, UserDataOperator::OPROWCOL),
728 dAta(data), commonData(common_data), aLe(false) {}
729
730template <int S>
732 int gg, MatrixDouble &jac_stress,
733 MatrixDouble &jac) {
734 jac.clear();
738 int nb_col = col_data.getFieldData().size();
739 double *diff_ptr =
740 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
741 // First two indices 'i','j' derivatives of 1st Piola-stress, third index 'k'
742 // is displacement component
744 &jac_stress(3 * 0 + 0, S + 0), &jac_stress(3 * 0 + 0, S + 1),
745 &jac_stress(3 * 0 + 0, S + 2), &jac_stress(3 * 0 + 1, S + 0),
746 &jac_stress(3 * 0 + 1, S + 1), &jac_stress(3 * 0 + 1, S + 2),
747 &jac_stress(3 * 0 + 2, S + 0), &jac_stress(3 * 0 + 2, S + 1),
748 &jac_stress(3 * 0 + 2, S + 2), &jac_stress(3 * 1 + 0, S + 0),
749 &jac_stress(3 * 1 + 0, S + 1), &jac_stress(3 * 1 + 0, S + 2),
750 &jac_stress(3 * 1 + 1, S + 0), &jac_stress(3 * 1 + 1, S + 1),
751 &jac_stress(3 * 1 + 1, S + 2), &jac_stress(3 * 1 + 2, S + 0),
752 &jac_stress(3 * 1 + 2, S + 1), &jac_stress(3 * 1 + 2, S + 2),
753 &jac_stress(3 * 2 + 0, S + 0), &jac_stress(3 * 2 + 0, S + 1),
754 &jac_stress(3 * 2 + 0, S + 2), &jac_stress(3 * 2 + 1, S + 0),
755 &jac_stress(3 * 2 + 1, S + 1), &jac_stress(3 * 2 + 1, S + 2),
756 &jac_stress(3 * 2 + 2, S + 0), &jac_stress(3 * 2 + 2, S + 1),
757 &jac_stress(3 * 2 + 2, S + 2));
759 &jac_stress(3 * 0 + 0, S + 3), &jac_stress(3 * 0 + 0, S + 4),
760 &jac_stress(3 * 0 + 0, S + 5), &jac_stress(3 * 0 + 1, S + 3),
761 &jac_stress(3 * 0 + 1, S + 4), &jac_stress(3 * 0 + 1, S + 5),
762 &jac_stress(3 * 0 + 2, S + 3), &jac_stress(3 * 0 + 2, S + 4),
763 &jac_stress(3 * 0 + 2, S + 5), &jac_stress(3 * 1 + 0, S + 3),
764 &jac_stress(3 * 1 + 0, S + 4), &jac_stress(3 * 1 + 0, S + 5),
765 &jac_stress(3 * 1 + 1, S + 3), &jac_stress(3 * 1 + 1, S + 4),
766 &jac_stress(3 * 1 + 1, S + 5), &jac_stress(3 * 1 + 2, S + 3),
767 &jac_stress(3 * 1 + 2, S + 4), &jac_stress(3 * 1 + 2, S + 5),
768 &jac_stress(3 * 2 + 0, S + 3), &jac_stress(3 * 2 + 0, S + 4),
769 &jac_stress(3 * 2 + 0, S + 5), &jac_stress(3 * 2 + 1, S + 3),
770 &jac_stress(3 * 2 + 1, S + 4), &jac_stress(3 * 2 + 1, S + 5),
771 &jac_stress(3 * 2 + 2, S + 3), &jac_stress(3 * 2 + 2, S + 4),
772 &jac_stress(3 * 2 + 2, S + 5));
774 &jac_stress(3 * 0 + 0, S + 6), &jac_stress(3 * 0 + 0, S + 7),
775 &jac_stress(3 * 0 + 0, S + 8), &jac_stress(3 * 0 + 1, S + 6),
776 &jac_stress(3 * 0 + 1, S + 7), &jac_stress(3 * 0 + 1, S + 8),
777 &jac_stress(3 * 0 + 2, S + 6), &jac_stress(3 * 0 + 2, S + 7),
778 &jac_stress(3 * 0 + 2, S + 8), &jac_stress(3 * 1 + 0, S + 6),
779 &jac_stress(3 * 1 + 0, S + 7), &jac_stress(3 * 1 + 0, S + 8),
780 &jac_stress(3 * 1 + 1, S + 6), &jac_stress(3 * 1 + 1, S + 7),
781 &jac_stress(3 * 1 + 1, S + 8), &jac_stress(3 * 1 + 2, S + 6),
782 &jac_stress(3 * 1 + 2, S + 7), &jac_stress(3 * 1 + 2, S + 8),
783 &jac_stress(3 * 2 + 0, S + 6), &jac_stress(3 * 2 + 0, S + 7),
784 &jac_stress(3 * 2 + 0, S + 8), &jac_stress(3 * 2 + 1, S + 6),
785 &jac_stress(3 * 2 + 1, S + 7), &jac_stress(3 * 2 + 1, S + 8),
786 &jac_stress(3 * 2 + 2, S + 6), &jac_stress(3 * 2 + 2, S + 7),
787 &jac_stress(3 * 2 + 2, S + 8));
788 // Derivate of 1st Piola-stress multiplied by gradient of defamation for
789 // base function (dd) and displacement component (rr)
791 &jac(0, 0), &jac(1, 0), &jac(2, 0), &jac(3, 0), &jac(4, 0), &jac(5, 0),
792 &jac(6, 0), &jac(7, 0), &jac(8, 0));
794 &jac(0, 1), &jac(1, 1), &jac(2, 1), &jac(3, 1), &jac(4, 1), &jac(5, 1),
795 &jac(6, 1), &jac(7, 1), &jac(8, 1));
797 &jac(0, 2), &jac(1, 2), &jac(2, 2), &jac(3, 2), &jac(4, 2), &jac(5, 2),
798 &jac(6, 2), &jac(7, 2), &jac(8, 2));
800 diff_ptr, &diff_ptr[1], &diff_ptr[2]);
801 for (int dd = 0; dd != nb_col / 3; ++dd) {
802 t2_1_0(i, j) += t3_1_0(i, j, k) * diff(k);
803 t2_1_1(i, j) += t3_1_1(i, j, k) * diff(k);
804 t2_1_2(i, j) += t3_1_2(i, j, k) * diff(k);
805 ++t2_1_0;
806 ++t2_1_1;
807 ++t2_1_2;
808 ++diff;
809 }
811}
812
814 EntitiesFieldData::EntData &col_data, int gg) {
815 return get_jac<0>(col_data, gg, commonData.jacStress[gg], jac);
816}
817
819 int row_side, int col_side, EntityType row_type, EntityType col_type,
821 EntitiesFieldData::EntData &col_data) {
823
824 int nb_row = row_data.getIndices().size();
825 int nb_col = col_data.getIndices().size();
826
827 int *row_indices_ptr = &row_data.getIndices()[0];
828 int *col_indices_ptr = &col_data.getIndices()[0];
829
830 if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
831 rowIndices.resize(nb_row, false);
832 noalias(rowIndices) = row_data.getIndices();
833 row_indices_ptr = &rowIndices[0];
834 VectorDofs &dofs = row_data.getFieldDofs();
835 VectorDofs::iterator dit = dofs.begin();
836 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
837 if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
838 dAta.forcesOnlyOnEntitiesRow.end()) {
839 rowIndices[ii] = -1;
840 }
841 }
842 }
843
844 if (!dAta.forcesOnlyOnEntitiesCol.empty()) {
845 colIndices.resize(nb_col, false);
846 noalias(colIndices) = col_data.getIndices();
847 col_indices_ptr = &colIndices[0];
848 VectorDofs &dofs = col_data.getFieldDofs();
849 VectorDofs::iterator dit = dofs.begin();
850 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
851 if (dAta.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
852 dAta.forcesOnlyOnEntitiesCol.end()) {
853 colIndices[ii] = -1;
854 }
855 }
856 }
857
858 CHKERR MatSetValues(getFEMethod()->snes_B, nb_row, row_indices_ptr, nb_col,
859 col_indices_ptr, &k(0, 0), ADD_VALUES);
860
861 // is symmetric
862 if (row_side != col_side || row_type != col_type) {
863
864 row_indices_ptr = &row_data.getIndices()[0];
865 col_indices_ptr = &col_data.getIndices()[0];
866
867 if (!dAta.forcesOnlyOnEntitiesCol.empty()) {
868 rowIndices.resize(nb_row, false);
869 noalias(rowIndices) = row_data.getIndices();
870 row_indices_ptr = &rowIndices[0];
871 VectorDofs &dofs = row_data.getFieldDofs();
872 VectorDofs::iterator dit = dofs.begin();
873 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
874 if (dAta.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
875 dAta.forcesOnlyOnEntitiesCol.end()) {
876 rowIndices[ii] = -1;
877 }
878 }
879 }
880
881 if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
882 colIndices.resize(nb_col, false);
883 noalias(colIndices) = col_data.getIndices();
884 col_indices_ptr = &colIndices[0];
885 VectorDofs &dofs = col_data.getFieldDofs();
886 VectorDofs::iterator dit = dofs.begin();
887 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
888 if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
889 dAta.forcesOnlyOnEntitiesRow.end()) {
890 colIndices[ii] = -1;
891 }
892 }
893 }
894
895 trans_k.resize(nb_col, nb_row, false);
896 noalias(trans_k) = trans(k);
897 CHKERR MatSetValues(getFEMethod()->snes_B, nb_col, col_indices_ptr, nb_row,
898 row_indices_ptr, &trans_k(0, 0), ADD_VALUES);
899 }
900
902}
903
905 int row_side, int col_side, EntityType row_type, EntityType col_type,
907 EntitiesFieldData::EntData &col_data) {
909
910 int nb_row = row_data.getIndices().size();
911 int nb_col = col_data.getIndices().size();
912 if (nb_row == 0)
914 if (nb_col == 0)
916
917 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
918 dAta.tEts.end()) {
920 }
921
922 // const int nb_base_functions = row_data.getN().size2();
923 const int nb_gauss_pts = row_data.getN().size1();
924
928
929 k.resize(nb_row, nb_col, false);
930 k.clear();
931 jac.resize(9, nb_col, false);
932
933 for (int gg = 0; gg != nb_gauss_pts; gg++) {
934 CHKERR getJac(col_data, gg);
935 double val = getVolume() * getGaussPts()(3, gg);
937 &jac(3 * 0 + 0, 0), &jac(3 * 0 + 0, 1), &jac(3 * 0 + 0, 2),
938 &jac(3 * 0 + 1, 0), &jac(3 * 0 + 1, 1), &jac(3 * 0 + 1, 2),
939 &jac(3 * 0 + 2, 0), &jac(3 * 0 + 2, 1), &jac(3 * 0 + 2, 2),
940 &jac(3 * 1 + 0, 0), &jac(3 * 1 + 0, 1), &jac(3 * 1 + 0, 2),
941 &jac(3 * 1 + 1, 0), &jac(3 * 1 + 1, 1), &jac(3 * 1 + 1, 2),
942 &jac(3 * 1 + 2, 0), &jac(3 * 1 + 2, 1), &jac(3 * 1 + 2, 2),
943 &jac(3 * 2 + 0, 0), &jac(3 * 2 + 0, 1), &jac(3 * 2 + 0, 2),
944 &jac(3 * 2 + 1, 0), &jac(3 * 2 + 1, 1), &jac(3 * 2 + 1, 2),
945 &jac(3 * 2 + 2, 0), &jac(3 * 2 + 2, 1), &jac(3 * 2 + 2, 2));
946 for (int cc = 0; cc != nb_col / 3; cc++) {
947 auto diff_base_functions = row_data.getFTensor1DiffN<3>(gg, 0);
949 &k(0, 3 * cc + 0), &k(0, 3 * cc + 1), &k(0, 3 * cc + 2),
950 &k(1, 3 * cc + 0), &k(1, 3 * cc + 1), &k(1, 3 * cc + 2),
951 &k(2, 3 * cc + 0), &k(2, 3 * cc + 1), &k(2, 3 * cc + 2), 3 * nb_col);
952 for (int rr = 0; rr != nb_row / 3; rr++) {
953 lhs(i, j) += val * t3_1(i, m, j) * diff_base_functions(m);
954 ++diff_base_functions;
955 ++lhs;
956 }
957 ++t3_1;
958 }
959 }
960
961 CHKERR aSemble(row_side, col_side, row_type, col_type, row_data, col_data);
962
964}
965
967 const std::string vel_field, const std::string field_name, BlockData &data,
968 CommonData &common_data)
969 : OpLhsPiolaKirchhoff_dx(vel_field, field_name, data, common_data) {
970 sYmm = false;
971}
972
974 EntitiesFieldData::EntData &col_data, int gg) {
975 return get_jac<9>(col_data, gg, commonData.jacStress[gg], jac);
976}
977
979 int row_side, int col_side, EntityType row_type, EntityType col_type,
981 EntitiesFieldData::EntData &col_data) {
983
984 int nb_row = row_data.getIndices().size();
985 int nb_col = col_data.getIndices().size();
986
987 int *row_indices_ptr = &row_data.getIndices()[0];
988 if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
989 rowIndices.resize(nb_row, false);
990 noalias(rowIndices) = row_data.getIndices();
991 row_indices_ptr = &rowIndices[0];
992 VectorDofs &dofs = row_data.getFieldDofs();
993 VectorDofs::iterator dit = dofs.begin();
994 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
995 if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
996 dAta.forcesOnlyOnEntitiesRow.end()) {
997 rowIndices[ii] = -1;
998 }
999 }
1000 }
1001
1002 int *col_indices_ptr = &col_data.getIndices()[0];
1003 if (!dAta.forcesOnlyOnEntitiesCol.empty()) {
1004 colIndices.resize(nb_col, false);
1005 noalias(colIndices) = col_data.getIndices();
1006 col_indices_ptr = &colIndices[0];
1007 VectorDofs &dofs = col_data.getFieldDofs();
1008 VectorDofs::iterator dit = dofs.begin();
1009 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
1010 if (dAta.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
1011 dAta.forcesOnlyOnEntitiesCol.end()) {
1012 colIndices[ii] = -1;
1013 }
1014 }
1015 }
1016
1017 /*for(int dd1 = 0;dd1<k.size1();dd1++) {
1018 for(int dd2 = 0;dd2<k.size2();dd2++) {
1019 if(k(dd1,dd2)!=k(dd1,dd2)) {
1020 SETERRQ(PETSC_COMM_SELF,1,"Wrong result");
1021 }
1022 }
1023 }*/
1024
1025 CHKERR MatSetValues(getFEMethod()->snes_B, nb_row, row_indices_ptr, nb_col,
1026 col_indices_ptr, &k(0, 0), ADD_VALUES);
1027
1029}
1030
1032 const std::string field_name, BlockData &data, CommonData &common_data,
1033 int tag, bool jacobian, bool ale)
1034 : OpJacobianPiolaKirchhoffStress(field_name, data, common_data, tag,
1035 jacobian, ale, false) {}
1036
1039 const int gg) {
1041
1042 CHKERR dAta.materialAdoublePtr->calculatesIGma_EshelbyStress(
1043 dAta, getNumeredEntFiniteElementPtr());
1044 if (aLe) {
1045 auto &t_sIGma = dAta.materialAdoublePtr->t_sIGma;
1046 auto &t_invH = dAta.materialAdoublePtr->t_invH;
1047 t_sIGma(i, j) = t_sIGma(i, k) * t_invH(j, k);
1048 t_sIGma(i, j) *= dAta.materialAdoublePtr->detH;
1049
1050 }
1051 commonData.sTress[gg].resize(3, 3, false);
1052 for (int dd1 = 0; dd1 < 3; dd1++) {
1053 for (int dd2 = 0; dd2 < 3; dd2++) {
1054 dAta.materialAdoublePtr->sIGma(dd1, dd2) >>=
1055 (commonData.sTress[gg])(dd1, dd2);
1056 }
1057 }
1058
1060}
1061
1063 const std::string field_name, BlockData &data, CommonData &common_data)
1064 : OpRhsPiolaKirchhoff(field_name, data, common_data) {}
1065
1067 const std::string vel_field, const std::string field_name, BlockData &data,
1068 CommonData &common_data)
1069 : OpLhsPiolaKirchhoff_dX(vel_field, field_name, data, common_data) {}
1070
1072 EntitiesFieldData::EntData &col_data, int gg) {
1073 return get_jac<0>(col_data, gg, commonData.jacStress[gg], jac);
1074}
1075
1077 const std::string vel_field, const std::string field_name, BlockData &data,
1078 CommonData &common_data)
1079 : OpLhsPiolaKirchhoff_dx(vel_field, field_name, data, common_data) {}
1080
1082 EntitiesFieldData::EntData &col_data, int gg) {
1083 return get_jac<9>(col_data, gg, commonData.jacStress[gg], jac);
1084}
1085
1088 materialDoublePtr,
1090 materialAdoublePtr) {
1092
1093 if (!materialDoublePtr) {
1095 "Pointer for materialDoublePtr not allocated");
1096 }
1097 if (!materialAdoublePtr) {
1099 "Pointer for materialAdoublePtr not allocated");
1100 }
1101
1103 mField, BLOCKSET | MAT_ELASTICSET, it)) {
1104 Mat_Elastic mydata;
1105 CHKERR it->getAttributeDataStructure(mydata);
1106 int id = it->getMeshsetId();
1107 EntityHandle meshset = it->getMeshset();
1108 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET,
1109 setOfBlocks[id].tEts, true);
1110 setOfBlocks[id].iD = id;
1111 setOfBlocks[id].E = mydata.data.Young;
1112 setOfBlocks[id].PoissonRatio = mydata.data.Poisson;
1113 setOfBlocks[id].materialDoublePtr = materialDoublePtr;
1114 setOfBlocks[id].materialAdoublePtr = materialAdoublePtr;
1115 }
1116
1118}
1119
1121 const std::string element_name,
1122 const std::string spatial_position_field_name,
1123 const std::string material_position_field_name, const bool ale) {
1125
1126 CHKERR mField.add_finite_element(element_name, MF_ZERO);
1128 element_name, spatial_position_field_name);
1130 element_name, spatial_position_field_name);
1132 element_name, spatial_position_field_name);
1133 if (mField.check_field(material_position_field_name)) {
1134 if (ale) {
1136 element_name, material_position_field_name);
1138 element_name, material_position_field_name);
1139 }
1141 element_name, material_position_field_name);
1142 }
1143
1144 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
1145 for (; sit != setOfBlocks.end(); sit++) {
1146 CHKERR mField.add_ents_to_finite_element_by_type(sit->second.tEts, MBTET,
1147 element_name);
1148 }
1149
1151}
1152
1154 const std::string spatial_position_field_name,
1155 const std::string material_position_field_name, const bool ale,
1156 const bool field_disp) {
1158
1159 commonData.spatialPositions = spatial_position_field_name;
1160 commonData.meshPositions = material_position_field_name;
1161
1162 // Rhs
1163 feRhs.getOpPtrVector().push_back(
1164 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
1165 if (mField.check_field(material_position_field_name)) {
1167 material_position_field_name, commonData));
1168 }
1169 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
1170 for (; sit != setOfBlocks.end(); sit++) {
1172 spatial_position_field_name, sit->second, commonData, tAg, false, ale,
1173 field_disp));
1175 spatial_position_field_name, sit->second, commonData));
1176 }
1177
1178 // Energy
1179 feEnergy.getOpPtrVector().push_back(
1180 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
1181 if (mField.check_field(material_position_field_name)) {
1183 material_position_field_name, commonData));
1184 }
1185 sit = setOfBlocks.begin();
1186 for (; sit != setOfBlocks.end(); sit++) {
1187 feEnergy.getOpPtrVector().push_back(
1188 new OpEnergy(spatial_position_field_name, sit->second, commonData,
1189 feEnergy.V, field_disp));
1190 }
1191
1192 // Lhs
1193 feLhs.getOpPtrVector().push_back(
1194 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
1195 if (mField.check_field(material_position_field_name)) {
1197 material_position_field_name, commonData));
1198 }
1199 sit = setOfBlocks.begin();
1200 for (; sit != setOfBlocks.end(); sit++) {
1202 spatial_position_field_name, sit->second, commonData, tAg, true, ale,
1203 field_disp));
1205 spatial_position_field_name, spatial_position_field_name, sit->second,
1206 commonData));
1207 }
1208
1210}
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
Definition: definitions.h:98
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:159
@ BLOCKSET
Definition: definitions.h:148
@ 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()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
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_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
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 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.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
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
auto createSmartVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
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.
Definition: Templates.hpp:1323
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1312
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
constexpr AssemblyType A
double H
Definition: plastic.cpp:100
constexpr auto field_name
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 filed data coeffinects.
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_vector< 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::vector< MatrixDouble3by3 > sTress
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
FTensor::Index< 'k', 3 > k
FTensor::Index< 'j', 3 > j
MoFEMErrorCode setBlocks(boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< double > > materialDoublePtr, boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< adouble > > materialAdoublePtr)
MyVolumeFE feEnergy
calculate elastic energy
FTensor::Index< 'i', 3 > i
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.