v0.13.1
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 * This file is part of MoFEM.
7 * MoFEM is free software: you can redistribute it and/or modify it under
8 * the terms of the GNU Lesser General Public License as published by the
9 * Free Software Foundation, either version 3 of the License, or (at your
10 * option) any later version.
11 *
12 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
13 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15 * License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
19
20#include <MoFEM.hpp>
21
22using namespace MoFEM;
24
25#include <adolc/adolc.h>
27
29 : VolumeElementForcesAndSourcesCore(m_field), A(PETSC_NULL), F(PETSC_NULL),
30 addToRule(1) {
31
32 auto create_vec = [&]() {
33 constexpr int ghosts[] = {0};
34 if (mField.get_comm_rank() == 0) {
35 return createSmartVectorMPI(mField.get_comm(), 1, 1);
36 } else {
37 return createSmartVectorMPI(mField.get_comm(), 0, 1);
38 }
39 };
40
41 V = create_vec();
42}
43
45 return 2 * (order - 1) + addToRule;
46};
47
50
51 CHKERR VolumeElementForcesAndSourcesCore::preProcess();
52
53 if (A != PETSC_NULL) {
54 snes_B = A;
55 }
56
57 if (F != PETSC_NULL) {
58 snes_f = F;
59 }
60
61 switch (snes_ctx) {
62 case CTX_SNESNONE:
63 CHKERR VecZeroEntries(V);
64 break;
65 default:
66 break;
67 }
68
70}
71
74
75 const double *array;
76
77 switch (snes_ctx) {
78 case CTX_SNESNONE:
79 CHKERR VecAssemblyBegin(V);
80 CHKERR VecAssemblyEnd(V);
81 CHKERR VecSum(V, &eNergy);
82 break;
83 default:
84 break;
85 }
86
87 CHKERR VolumeElementForcesAndSourcesCore::postProcess();
88
90}
91
93 short int tag)
94 : feRhs(m_field), feLhs(m_field), feEnergy(m_field), mField(m_field),
95 tAg(tag) {}
96
98 const std::string field_name,
99 std::vector<VectorDouble> &values_at_gauss_pts,
100 std::vector<MatrixDouble> &gardient_at_gauss_pts)
103 valuesAtGaussPts(values_at_gauss_pts),
104 gradientAtGaussPts(gardient_at_gauss_pts), zeroAtType(MBVERTEX) {}
105
107 int side, EntityType type, EntitiesFieldData::EntData &data) {
109
110 const int nb_dofs = data.getFieldData().size();
111 const int nb_base_functions = data.getN().size2();
112 if (nb_dofs == 0) {
114 }
115 const int nb_gauss_pts = data.getN().size1();
116 const int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
117
118 // initialize
119 if (type == zeroAtType) {
120 valuesAtGaussPts.resize(nb_gauss_pts);
121 gradientAtGaussPts.resize(nb_gauss_pts);
122 for (int gg = 0; gg != nb_gauss_pts; gg++) {
123 valuesAtGaussPts[gg].resize(rank, false);
124 gradientAtGaussPts[gg].resize(rank, 3, false);
125 }
126 for (int gg = 0; gg != nb_gauss_pts; gg++) {
127 valuesAtGaussPts[gg].clear();
128 gradientAtGaussPts[gg].clear();
129 }
130 }
131
132 auto base_function = data.getFTensor0N();
133 auto diff_base_functions = data.getFTensor1DiffN<3>();
136
137 if (rank == 1) {
138
139 for (int gg = 0; gg != nb_gauss_pts; gg++) {
140 auto field_data = data.getFTensor0FieldData();
141 double &val = valuesAtGaussPts[gg][0];
142 FTensor::Tensor1<double *, 3> grad(&gradientAtGaussPts[gg](0, 0),
143 &gradientAtGaussPts[gg](0, 1),
144 &gradientAtGaussPts[gg](0, 2));
145 int bb = 0;
146 for (; bb != nb_dofs; bb++) {
147 val += base_function * field_data;
148 grad(i) += diff_base_functions(i) * field_data;
149 ++diff_base_functions;
150 ++base_function;
151 ++field_data;
152 }
153 for (; bb != nb_base_functions; bb++) {
154 ++diff_base_functions;
155 ++base_function;
156 }
157 }
158
159 } else if (rank == 3) {
160
161 for (int gg = 0; gg != nb_gauss_pts; gg++) {
162 auto field_data = data.getFTensor1FieldData<3>();
163 FTensor::Tensor1<double *, 3> values(&valuesAtGaussPts[gg][0],
164 &valuesAtGaussPts[gg][1],
165 &valuesAtGaussPts[gg][2]);
167 &gradientAtGaussPts[gg](0, 0), &gradientAtGaussPts[gg](0, 1),
168 &gradientAtGaussPts[gg](0, 2), &gradientAtGaussPts[gg](1, 0),
169 &gradientAtGaussPts[gg](1, 1), &gradientAtGaussPts[gg](1, 2),
170 &gradientAtGaussPts[gg](2, 0), &gradientAtGaussPts[gg](2, 1),
171 &gradientAtGaussPts[gg](2, 2));
172 int bb = 0;
173 for (; bb != nb_dofs / 3; bb++) {
174 values(i) += base_function * field_data(i);
175 gradient(i, j) += field_data(i) * diff_base_functions(j);
176 ++diff_base_functions;
177 ++base_function;
178 ++field_data;
179 }
180 for (; bb != nb_base_functions; bb++) {
181 ++diff_base_functions;
182 ++base_function;
183 }
184 }
185
186 } else {
187 // FIXME: THat part is inefficient
188 VectorDouble &values = data.getFieldData();
189 for (int gg = 0; gg < nb_gauss_pts; gg++) {
190 VectorAdaptor N = data.getN(gg, nb_dofs / rank);
191 MatrixAdaptor diffN = data.getDiffN(gg, nb_dofs / rank);
192 for (int dd = 0; dd < nb_dofs / rank; dd++) {
193 for (int rr1 = 0; rr1 < rank; rr1++) {
194 valuesAtGaussPts[gg][rr1] += N[dd] * values[rank * dd + rr1];
195 for (int rr2 = 0; rr2 < 3; rr2++) {
196 gradientAtGaussPts[gg](rr1, rr2) +=
197 diffN(dd, rr2) * values[rank * dd + rr1];
198 }
199 }
200 }
201 }
202 }
203
205}
206
208 const std::string field_name, CommonData &common_data)
209 : OpGetDataAtGaussPts(field_name, common_data.dataAtGaussPts[field_name],
210 common_data.gradAtGaussPts[field_name]) {}
211
214 BlockData &data, CommonData &common_data,
215 int tag, bool jacobian, bool ale,
216 bool field_disp)
219 dAta(data), commonData(common_data), tAg(tag), adlocReturnValue(0),
220 jAcobian(jacobian), fUnction(!jacobian), aLe(ale), fieldDisp(field_disp) {
221
222}
223
226 const int gg) {
228
229 CHKERR dAta.materialAdoublePtr->calculateP_PiolaKirchhoffI(
230 dAta, getNumeredEntFiniteElementPtr());
231
232 if (aLe) {
233 auto &t_P = dAta.materialAdoublePtr->t_P;
234 auto &t_invH = dAta.materialAdoublePtr->t_invH;
235 t_P(i, j) = t_P(i, k) * t_invH(j, k);
236 t_P(i, j) *= dAta.materialAdoublePtr->detH;
237 }
238
239 commonData.sTress[gg].resize(3, 3, false);
240 for (int dd1 = 0; dd1 < 3; dd1++) {
241 for (int dd2 = 0; dd2 < 3; dd2++) {
242 dAta.materialAdoublePtr->P(dd1, dd2) >>=
243 (commonData.sTress[gg])(dd1, dd2);
244 }
245 }
246
248}
249
252 const int gg) {
254
255 trace_on(tAg, 0);
256
257 dAta.materialAdoublePtr->F.resize(3, 3, false);
258
259 if (!aLe) {
260
261 nbActiveVariables = 0;
262 for (int dd1 = 0; dd1 < 3; dd1++) {
263 for (int dd2 = 0; dd2 < 3; dd2++) {
264 dAta.materialAdoublePtr->F(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
265 if (fieldDisp) {
266 if (dd1 == dd2) {
267 dAta.materialAdoublePtr->F(dd1, dd2) += 1;
268 }
269 }
270 nbActiveVariables++;
271 }
272 }
273
274 } else {
275
276 nbActiveVariables = 0;
277
278 dAta.materialAdoublePtr->h.resize(3, 3, false);
279 for (int dd1 = 0; dd1 < 3; dd1++) {
280 for (int dd2 = 0; dd2 < 3; dd2++) {
281 dAta.materialAdoublePtr->h(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
282 nbActiveVariables++;
283 }
284 }
285
286 dAta.materialAdoublePtr->H.resize(3, 3, false);
287 for (int dd1 = 0; dd1 < 3; dd1++) {
288 for (int dd2 = 0; dd2 < 3; dd2++) {
289 dAta.materialAdoublePtr->H(dd1, dd2) <<= (*ptrH)[gg](dd1, dd2);
290 nbActiveVariables++;
291 }
292 }
293
294 dAta.materialAdoublePtr->detH = determinantTensor3by3(dAta.materialAdoublePtr->H);
295 dAta.materialAdoublePtr->invH.resize(3, 3, false);
296 CHKERR invertTensor3by3(dAta.materialAdoublePtr->H,
297 dAta.materialAdoublePtr->detH,
298 dAta.materialAdoublePtr->invH);
299
300 auto &t_F = dAta.materialAdoublePtr->t_F;
301 auto &t_h = dAta.materialAdoublePtr->t_h;
302 auto &t_invH = dAta.materialAdoublePtr->t_invH;
303
304 t_F(i, j) = t_h(i, k) * t_invH(k, j);
305
306 }
307
308 CHKERR dAta.materialAdoublePtr->setUserActiveVariables(nbActiveVariables);
309 CHKERR calculateStress(gg);
310
311 trace_off();
312
314}
315
319
320 int r;
321
322 if (fUnction) {
323 commonData.sTress[gg].resize(3, 3, false);
324 // play recorder for values
325 r = ::function(tAg, 9, nbActiveVariables, &activeVariables[0],
326 &commonData.sTress[gg](0, 0));
327 if (r < adlocReturnValue) { // function is locally analytic
328 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
329 "ADOL-C function evaluation with error r = %d", r);
330 }
331 }
332
333 if (jAcobian) {
334 commonData.jacStress[gg].resize(9, nbActiveVariables, false);
335 double *jac_ptr[] = {
336 &(commonData.jacStress[gg](0, 0)), &(commonData.jacStress[gg](1, 0)),
337 &(commonData.jacStress[gg](2, 0)), &(commonData.jacStress[gg](3, 0)),
338 &(commonData.jacStress[gg](4, 0)), &(commonData.jacStress[gg](5, 0)),
339 &(commonData.jacStress[gg](6, 0)), &(commonData.jacStress[gg](7, 0)),
340 &(commonData.jacStress[gg](8, 0))};
341 // play recorder for jacobians
342 r = jacobian(tAg, 9, nbActiveVariables, &activeVariables[0], jac_ptr);
343 if (r < adlocReturnValue) {
344 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
345 "ADOL-C function evaluation with error");
346 }
347 }
348
350}
351
353 int row_side, EntityType row_type,
354 EntitiesFieldData::EntData &row_data) {
356
357 // do it only once, no need to repeat this for edges,faces or tets
358 if (row_type != MBVERTEX)
360
361 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
362 dAta.tEts.end()) {
364 }
365
366 int nb_dofs = row_data.getFieldData().size();
367 if (nb_dofs == 0)
369 dAta.materialAdoublePtr->commonDataPtr = &commonData;
370 dAta.materialAdoublePtr->opPtr = this;
371
372 int nb_gauss_pts = row_data.getN().size1();
373 commonData.sTress.resize(nb_gauss_pts);
374 commonData.jacStress.resize(nb_gauss_pts);
375
377 if (aLe) {
379 }
380
381 for (int gg = 0; gg != nb_gauss_pts; gg++) {
382
383 dAta.materialAdoublePtr->gG = gg;
384
385 // Record tag and calculate stress
386 if (recordTagForIntegrationPoint(gg)) {
387 CHKERR recordTag(gg);
388 }
389
390 // Set active variables vector
391 if (jAcobian || (!recordTagForIntegrationPoint(gg))) {
392 activeVariables.resize(nbActiveVariables, false);
393 if (!aLe) {
394 for (int dd1 = 0; dd1 < 3; dd1++) {
395 for (int dd2 = 0; dd2 < 3; dd2++) {
396 activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
397 }
398 }
399 } else {
400 for (int dd1 = 0; dd1 < 3; dd1++) {
401 for (int dd2 = 0; dd2 < 3; dd2++) {
402 activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
403 }
404 }
405 for (int dd1 = 0; dd1 < 3; dd1++) {
406 for (int dd2 = 0; dd2 < 3; dd2++) {
407 activeVariables(9 + dd1 * 3 + dd2) = (*ptrH)[gg](dd1, dd2);
408 }
409 }
410 }
411 CHKERR dAta.materialAdoublePtr->setUserActiveVariables(activeVariables);
412
413 // Play tag and calculate stress or tangent
414 if (jAcobian || (!recordTagForIntegrationPoint(gg))) {
415 CHKERR playTag(gg);
416 }
417 }
418 }
419
421}
422
424 const std::string
425 field_name, ///< field name for spatial positions or displacements
426 BlockData &data, CommonData &common_data, int tag, bool gradient,
427 bool hessian, bool ale, bool field_disp)
430 dAta(data), commonData(common_data), tAg(tag), gRadient(gradient),
431 hEssian(hessian), aLe(ale), fieldDisp(field_disp) {}
432
436 CHKERR dAta.materialAdoublePtr->calculateElasticEnergy(
437 dAta, getNumeredEntFiniteElementPtr());
438 dAta.materialAdoublePtr->eNergy >>= commonData.eNergy[gg];
440}
441
445
446 trace_on(tAg, 0);
447
448 if (!aLe) {
449
450 nbActiveVariables = 0;
451 for (int dd1 = 0; dd1 < 3; dd1++) {
452 for (int dd2 = 0; dd2 < 3; dd2++) {
453 dAta.materialAdoublePtr->F(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
454 if (fieldDisp) {
455 if (dd1 == dd2) {
456 dAta.materialAdoublePtr->F(dd1, dd2) += 1;
457 }
458 }
459 nbActiveVariables++;
460 }
461 }
462
463 } else {
464
465 nbActiveVariables = 0;
466
467 dAta.materialAdoublePtr->h.resize(3, 3, false);
468 for (int dd1 = 0; dd1 < 3; dd1++) {
469 for (int dd2 = 0; dd2 < 3; dd2++) {
470 dAta.materialAdoublePtr->h(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
471 nbActiveVariables++;
472 }
473 }
474
475 dAta.materialAdoublePtr->H.resize(3, 3, false);
476 for (int dd1 = 0; dd1 < 3; dd1++) {
477 for (int dd2 = 0; dd2 < 3; dd2++) {
478 dAta.materialAdoublePtr->H(dd1, dd2) <<= (*ptrH)[gg](dd1, dd2);
479 nbActiveVariables++;
480 }
481 }
482
483 dAta.materialAdoublePtr->detH = determinantTensor3by3(dAta.materialAdoublePtr->H);
484 dAta.materialAdoublePtr->invH.resize(3, 3, false);
485 CHKERR invertTensor3by3(dAta.materialAdoublePtr->H,
486 dAta.materialAdoublePtr->detH,
487 dAta.materialAdoublePtr->invH);
488
489 auto &t_F = dAta.materialAdoublePtr->t_F;
490 auto &t_h = dAta.materialAdoublePtr->t_h;
491 auto &t_invH = dAta.materialAdoublePtr->t_invH;
492
493 t_F(i, j) = t_h(i, k) * t_invH(k, j);
494
495 }
496
497 CHKERR dAta.materialAdoublePtr->setUserActiveVariables(nbActiveVariables);
499
500 trace_off();
501
503}
504
508
509 if (gRadient) {
510 commonData.jacEnergy[gg].resize(nbActiveVariables, false);
511 int r = ::gradient(tAg, nbActiveVariables, &activeVariables[0],
512 &commonData.jacEnergy[gg][0]);
513 if (r < 0) {
514 // That means that energy function is not smooth and derivative
515 // can not be calculated,
516 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
517 "ADOL-C function evaluation with error");
518 }
519 }
520
521 if (hEssian) {
522 commonData.hessianEnergy[gg].resize(nbActiveVariables * nbActiveVariables,
523 false);
524 double *H[nbActiveVariables];
525 for (int n = 0; n != nbActiveVariables; n++) {
526 H[n] = &(commonData.hessianEnergy[gg][n * nbActiveVariables]);
527 }
528 int r = ::hessian(tAg, nbActiveVariables, &*activeVariables.begin(), H);
529 if (r < 0) {
530 // That means that energy function is not smooth and derivative
531 // can not be calculated,
532 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
533 "ADOL-C function evaluation with error");
534 }
535 }
536
538}
539
541 int row_side, EntityType row_type,
542 EntitiesFieldData::EntData &row_data) {
544
545 // do it only once, no need to repeat this for edges,faces or tets
546 if (row_type != MBVERTEX)
548
549 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
550 dAta.tEts.end()) {
552 }
553
554 int nb_dofs = row_data.getFieldData().size();
555 if (nb_dofs == 0)
557 dAta.materialAdoublePtr->commonDataPtr = &commonData;
558 dAta.materialAdoublePtr->opPtr = this;
559
560 int nb_gauss_pts = row_data.getN().size1();
561 commonData.eNergy.resize(nb_gauss_pts);
562 commonData.jacEnergy.resize(nb_gauss_pts);
563
565 if (aLe) {
567 }
568
569 for (int gg = 0; gg != nb_gauss_pts; gg++) {
570
571 dAta.materialAdoublePtr->gG = gg;
572
573 // Record tag and calualte stress
574 if (recordTagForIntegrationPoint(gg)) {
575 CHKERR recordTag(gg);
576 }
577
578 activeVariables.resize(nbActiveVariables, false);
579 if (!aLe) {
580 for (int dd1 = 0; dd1 < 3; dd1++) {
581 for (int dd2 = 0; dd2 < 3; dd2++) {
582 activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
583 }
584 }
585 } else {
586 for (int dd1 = 0; dd1 < 3; dd1++) {
587 for (int dd2 = 0; dd2 < 3; dd2++) {
588 activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
589 }
590 }
591 for (int dd1 = 0; dd1 < 3; dd1++) {
592 for (int dd2 = 0; dd2 < 3; dd2++) {
593 activeVariables(9 + dd1 * 3 + dd2) = (*ptrH)[gg](dd1, dd2);
594 }
595 }
596 }
597 CHKERR dAta.materialAdoublePtr->setUserActiveVariables(activeVariables);
598
599 // Play tag and calculate stress or tangent
600 CHKERR playTag(gg);
601 }
602
604}
605
607 const std::string field_name, BlockData &data, CommonData &common_data)
610 dAta(data), commonData(common_data), aLe(false) {}
611
613 int row_side, EntityType row_type,
614 EntitiesFieldData::EntData &row_data) {
616
617 int nb_dofs = row_data.getIndices().size();
618 int *indices_ptr = &row_data.getIndices()[0];
619 if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
620 iNdices.resize(nb_dofs, false);
621 noalias(iNdices) = row_data.getIndices();
622 indices_ptr = &iNdices[0];
623 VectorDofs &dofs = row_data.getFieldDofs();
624 VectorDofs::iterator dit = dofs.begin();
625 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
626 if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
627 dAta.forcesOnlyOnEntitiesRow.end()) {
628 iNdices[ii] = -1;
629 }
630 }
631 }
632 CHKERR VecSetValues(getFEMethod()->snes_f, nb_dofs, indices_ptr, &nf[0],
633 ADD_VALUES);
635}
636
638 int row_side, EntityType row_type,
639 EntitiesFieldData::EntData &row_data) {
641
642 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
643 dAta.tEts.end()) {
645 }
646
647 const int nb_dofs = row_data.getIndices().size();
648 if (nb_dofs == 0)
650 if ((unsigned int)nb_dofs > 3 * row_data.getN().size2()) {
651 SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
652 }
653 const int nb_base_functions = row_data.getN().size2();
654 const int nb_gauss_pts = row_data.getN().size1();
655
656 nf.resize(nb_dofs, false);
657 nf.clear();
658
659 auto diff_base_functions = row_data.getFTensor1DiffN<3>();
662
663 for (int gg = 0; gg != nb_gauss_pts; gg++) {
664 double val = getVolume() * getGaussPts()(3, gg);
665 MatrixDouble3by3 &stress = commonData.sTress[gg];
667 &stress(0, 0), &stress(0, 1), &stress(0, 2), &stress(1, 0),
668 &stress(1, 1), &stress(1, 2), &stress(2, 0), &stress(2, 1),
669 &stress(2, 2));
670 FTensor::Tensor1<double *, 3> rhs(&nf[0], &nf[1], &nf[2], 3);
671 int bb = 0;
672 for (; bb != nb_dofs / 3; bb++) {
673 rhs(i) += val * t3(i, j) * diff_base_functions(j);
674 ++rhs;
675 ++diff_base_functions;
676 }
677 for (; bb != nb_base_functions; bb++) {
678 ++diff_base_functions;
679 }
680 }
681
682 CHKERR aSemble(row_side, row_type, row_data);
683
685}
686
688 BlockData &data,
689 CommonData &common_data,
690 SmartPetscObj<Vec> ghost_vec,
691 bool field_disp)
694 dAta(data), commonData(common_data), ghostVec(ghost_vec, true),
695 fieldDisp(field_disp) {}
696
698 int row_side, EntityType row_type,
699 EntitiesFieldData::EntData &row_data) {
701
702 if (row_type != MBVERTEX)
704 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
705 dAta.tEts.end()) {
707 }
708
709 std::vector<MatrixDouble> &F =
711 dAta.materialDoublePtr->F.resize(3, 3, false);
712
713 double energy = 0;
714
715 for (unsigned int gg = 0; gg != row_data.getN().size1(); ++gg) {
716 double val = getVolume() * getGaussPts()(3, gg);
717 noalias(dAta.materialDoublePtr->F) = F[gg];
718 if (fieldDisp) {
719 for (int dd = 0; dd < 3; dd++) {
720 dAta.materialDoublePtr->F(dd, dd) += 1;
721 }
722 }
723 int nb_active_variables = 0;
724 CHKERR dAta.materialDoublePtr->setUserActiveVariables(nb_active_variables);
725 CHKERR dAta.materialDoublePtr->calculateElasticEnergy(
726 dAta, getNumeredEntFiniteElementPtr());
727 energy += val * dAta.materialDoublePtr->eNergy;
728 }
729
730 CHKERR VecSetValue(ghostVec, 0, energy, ADD_VALUES);
732}
733
735 const std::string vel_field, const std::string field_name, BlockData &data,
736 CommonData &common_data)
738 vel_field, field_name, UserDataOperator::OPROWCOL),
739 dAta(data), commonData(common_data), aLe(false) {}
740
741template <int S>
743 int gg, MatrixDouble &jac_stress,
744 MatrixDouble &jac) {
745 jac.clear();
749 int nb_col = col_data.getFieldData().size();
750 double *diff_ptr =
751 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
752 // First two indices 'i','j' derivatives of 1st Piola-stress, third index 'k'
753 // is displacement component
755 &jac_stress(3 * 0 + 0, S + 0), &jac_stress(3 * 0 + 0, S + 1),
756 &jac_stress(3 * 0 + 0, S + 2), &jac_stress(3 * 0 + 1, S + 0),
757 &jac_stress(3 * 0 + 1, S + 1), &jac_stress(3 * 0 + 1, S + 2),
758 &jac_stress(3 * 0 + 2, S + 0), &jac_stress(3 * 0 + 2, S + 1),
759 &jac_stress(3 * 0 + 2, S + 2), &jac_stress(3 * 1 + 0, S + 0),
760 &jac_stress(3 * 1 + 0, S + 1), &jac_stress(3 * 1 + 0, S + 2),
761 &jac_stress(3 * 1 + 1, S + 0), &jac_stress(3 * 1 + 1, S + 1),
762 &jac_stress(3 * 1 + 1, S + 2), &jac_stress(3 * 1 + 2, S + 0),
763 &jac_stress(3 * 1 + 2, S + 1), &jac_stress(3 * 1 + 2, S + 2),
764 &jac_stress(3 * 2 + 0, S + 0), &jac_stress(3 * 2 + 0, S + 1),
765 &jac_stress(3 * 2 + 0, S + 2), &jac_stress(3 * 2 + 1, S + 0),
766 &jac_stress(3 * 2 + 1, S + 1), &jac_stress(3 * 2 + 1, S + 2),
767 &jac_stress(3 * 2 + 2, S + 0), &jac_stress(3 * 2 + 2, S + 1),
768 &jac_stress(3 * 2 + 2, S + 2));
770 &jac_stress(3 * 0 + 0, S + 3), &jac_stress(3 * 0 + 0, S + 4),
771 &jac_stress(3 * 0 + 0, S + 5), &jac_stress(3 * 0 + 1, S + 3),
772 &jac_stress(3 * 0 + 1, S + 4), &jac_stress(3 * 0 + 1, S + 5),
773 &jac_stress(3 * 0 + 2, S + 3), &jac_stress(3 * 0 + 2, S + 4),
774 &jac_stress(3 * 0 + 2, S + 5), &jac_stress(3 * 1 + 0, S + 3),
775 &jac_stress(3 * 1 + 0, S + 4), &jac_stress(3 * 1 + 0, S + 5),
776 &jac_stress(3 * 1 + 1, S + 3), &jac_stress(3 * 1 + 1, S + 4),
777 &jac_stress(3 * 1 + 1, S + 5), &jac_stress(3 * 1 + 2, S + 3),
778 &jac_stress(3 * 1 + 2, S + 4), &jac_stress(3 * 1 + 2, S + 5),
779 &jac_stress(3 * 2 + 0, S + 3), &jac_stress(3 * 2 + 0, S + 4),
780 &jac_stress(3 * 2 + 0, S + 5), &jac_stress(3 * 2 + 1, S + 3),
781 &jac_stress(3 * 2 + 1, S + 4), &jac_stress(3 * 2 + 1, S + 5),
782 &jac_stress(3 * 2 + 2, S + 3), &jac_stress(3 * 2 + 2, S + 4),
783 &jac_stress(3 * 2 + 2, S + 5));
785 &jac_stress(3 * 0 + 0, S + 6), &jac_stress(3 * 0 + 0, S + 7),
786 &jac_stress(3 * 0 + 0, S + 8), &jac_stress(3 * 0 + 1, S + 6),
787 &jac_stress(3 * 0 + 1, S + 7), &jac_stress(3 * 0 + 1, S + 8),
788 &jac_stress(3 * 0 + 2, S + 6), &jac_stress(3 * 0 + 2, S + 7),
789 &jac_stress(3 * 0 + 2, S + 8), &jac_stress(3 * 1 + 0, S + 6),
790 &jac_stress(3 * 1 + 0, S + 7), &jac_stress(3 * 1 + 0, S + 8),
791 &jac_stress(3 * 1 + 1, S + 6), &jac_stress(3 * 1 + 1, S + 7),
792 &jac_stress(3 * 1 + 1, S + 8), &jac_stress(3 * 1 + 2, S + 6),
793 &jac_stress(3 * 1 + 2, S + 7), &jac_stress(3 * 1 + 2, S + 8),
794 &jac_stress(3 * 2 + 0, S + 6), &jac_stress(3 * 2 + 0, S + 7),
795 &jac_stress(3 * 2 + 0, S + 8), &jac_stress(3 * 2 + 1, S + 6),
796 &jac_stress(3 * 2 + 1, S + 7), &jac_stress(3 * 2 + 1, S + 8),
797 &jac_stress(3 * 2 + 2, S + 6), &jac_stress(3 * 2 + 2, S + 7),
798 &jac_stress(3 * 2 + 2, S + 8));
799 // Derivate of 1st Piola-stress multiplied by gradient of defamation for
800 // base function (dd) and displacement component (rr)
802 &jac(0, 0), &jac(1, 0), &jac(2, 0), &jac(3, 0), &jac(4, 0), &jac(5, 0),
803 &jac(6, 0), &jac(7, 0), &jac(8, 0));
805 &jac(0, 1), &jac(1, 1), &jac(2, 1), &jac(3, 1), &jac(4, 1), &jac(5, 1),
806 &jac(6, 1), &jac(7, 1), &jac(8, 1));
808 &jac(0, 2), &jac(1, 2), &jac(2, 2), &jac(3, 2), &jac(4, 2), &jac(5, 2),
809 &jac(6, 2), &jac(7, 2), &jac(8, 2));
811 diff_ptr, &diff_ptr[1], &diff_ptr[2]);
812 for (int dd = 0; dd != nb_col / 3; ++dd) {
813 t2_1_0(i, j) += t3_1_0(i, j, k) * diff(k);
814 t2_1_1(i, j) += t3_1_1(i, j, k) * diff(k);
815 t2_1_2(i, j) += t3_1_2(i, j, k) * diff(k);
816 ++t2_1_0;
817 ++t2_1_1;
818 ++t2_1_2;
819 ++diff;
820 }
822}
823
825 EntitiesFieldData::EntData &col_data, int gg) {
826 return get_jac<0>(col_data, gg, commonData.jacStress[gg], jac);
827}
828
830 int row_side, int col_side, EntityType row_type, EntityType col_type,
832 EntitiesFieldData::EntData &col_data) {
834
835 int nb_row = row_data.getIndices().size();
836 int nb_col = col_data.getIndices().size();
837
838 int *row_indices_ptr = &row_data.getIndices()[0];
839 int *col_indices_ptr = &col_data.getIndices()[0];
840
841 if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
842 rowIndices.resize(nb_row, false);
843 noalias(rowIndices) = row_data.getIndices();
844 row_indices_ptr = &rowIndices[0];
845 VectorDofs &dofs = row_data.getFieldDofs();
846 VectorDofs::iterator dit = dofs.begin();
847 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
848 if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
849 dAta.forcesOnlyOnEntitiesRow.end()) {
850 rowIndices[ii] = -1;
851 }
852 }
853 }
854
855 if (!dAta.forcesOnlyOnEntitiesCol.empty()) {
856 colIndices.resize(nb_col, false);
857 noalias(colIndices) = col_data.getIndices();
858 col_indices_ptr = &colIndices[0];
859 VectorDofs &dofs = col_data.getFieldDofs();
860 VectorDofs::iterator dit = dofs.begin();
861 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
862 if (dAta.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
863 dAta.forcesOnlyOnEntitiesCol.end()) {
864 colIndices[ii] = -1;
865 }
866 }
867 }
868
869 CHKERR MatSetValues(getFEMethod()->snes_B, nb_row, row_indices_ptr, nb_col,
870 col_indices_ptr, &k(0, 0), ADD_VALUES);
871
872 // is symmetric
873 if (row_side != col_side || row_type != col_type) {
874
875 row_indices_ptr = &row_data.getIndices()[0];
876 col_indices_ptr = &col_data.getIndices()[0];
877
878 if (!dAta.forcesOnlyOnEntitiesCol.empty()) {
879 rowIndices.resize(nb_row, false);
880 noalias(rowIndices) = row_data.getIndices();
881 row_indices_ptr = &rowIndices[0];
882 VectorDofs &dofs = row_data.getFieldDofs();
883 VectorDofs::iterator dit = dofs.begin();
884 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
885 if (dAta.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
886 dAta.forcesOnlyOnEntitiesCol.end()) {
887 rowIndices[ii] = -1;
888 }
889 }
890 }
891
892 if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
893 colIndices.resize(nb_col, false);
894 noalias(colIndices) = col_data.getIndices();
895 col_indices_ptr = &colIndices[0];
896 VectorDofs &dofs = col_data.getFieldDofs();
897 VectorDofs::iterator dit = dofs.begin();
898 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
899 if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
900 dAta.forcesOnlyOnEntitiesRow.end()) {
901 colIndices[ii] = -1;
902 }
903 }
904 }
905
906 trans_k.resize(nb_col, nb_row, false);
907 noalias(trans_k) = trans(k);
908 CHKERR MatSetValues(getFEMethod()->snes_B, nb_col, col_indices_ptr, nb_row,
909 row_indices_ptr, &trans_k(0, 0), ADD_VALUES);
910 }
911
913}
914
916 int row_side, int col_side, EntityType row_type, EntityType col_type,
918 EntitiesFieldData::EntData &col_data) {
920
921 int nb_row = row_data.getIndices().size();
922 int nb_col = col_data.getIndices().size();
923 if (nb_row == 0)
925 if (nb_col == 0)
927
928 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
929 dAta.tEts.end()) {
931 }
932
933 // const int nb_base_functions = row_data.getN().size2();
934 const int nb_gauss_pts = row_data.getN().size1();
935
939
940 k.resize(nb_row, nb_col, false);
941 k.clear();
942 jac.resize(9, nb_col, false);
943
944 for (int gg = 0; gg != nb_gauss_pts; gg++) {
945 CHKERR getJac(col_data, gg);
946 double val = getVolume() * getGaussPts()(3, gg);
948 &jac(3 * 0 + 0, 0), &jac(3 * 0 + 0, 1), &jac(3 * 0 + 0, 2),
949 &jac(3 * 0 + 1, 0), &jac(3 * 0 + 1, 1), &jac(3 * 0 + 1, 2),
950 &jac(3 * 0 + 2, 0), &jac(3 * 0 + 2, 1), &jac(3 * 0 + 2, 2),
951 &jac(3 * 1 + 0, 0), &jac(3 * 1 + 0, 1), &jac(3 * 1 + 0, 2),
952 &jac(3 * 1 + 1, 0), &jac(3 * 1 + 1, 1), &jac(3 * 1 + 1, 2),
953 &jac(3 * 1 + 2, 0), &jac(3 * 1 + 2, 1), &jac(3 * 1 + 2, 2),
954 &jac(3 * 2 + 0, 0), &jac(3 * 2 + 0, 1), &jac(3 * 2 + 0, 2),
955 &jac(3 * 2 + 1, 0), &jac(3 * 2 + 1, 1), &jac(3 * 2 + 1, 2),
956 &jac(3 * 2 + 2, 0), &jac(3 * 2 + 2, 1), &jac(3 * 2 + 2, 2));
957 for (int cc = 0; cc != nb_col / 3; cc++) {
958 auto diff_base_functions = row_data.getFTensor1DiffN<3>(gg, 0);
960 &k(0, 3 * cc + 0), &k(0, 3 * cc + 1), &k(0, 3 * cc + 2),
961 &k(1, 3 * cc + 0), &k(1, 3 * cc + 1), &k(1, 3 * cc + 2),
962 &k(2, 3 * cc + 0), &k(2, 3 * cc + 1), &k(2, 3 * cc + 2), 3 * nb_col);
963 for (int rr = 0; rr != nb_row / 3; rr++) {
964 lhs(i, j) += val * t3_1(i, m, j) * diff_base_functions(m);
965 ++diff_base_functions;
966 ++lhs;
967 }
968 ++t3_1;
969 }
970 }
971
972 CHKERR aSemble(row_side, col_side, row_type, col_type, row_data, col_data);
973
975}
976
978 const std::string vel_field, const std::string field_name, BlockData &data,
979 CommonData &common_data)
980 : OpLhsPiolaKirchhoff_dx(vel_field, field_name, data, common_data) {
981 sYmm = false;
982}
983
985 EntitiesFieldData::EntData &col_data, int gg) {
986 return get_jac<9>(col_data, gg, commonData.jacStress[gg], jac);
987}
988
990 int row_side, int col_side, EntityType row_type, EntityType col_type,
992 EntitiesFieldData::EntData &col_data) {
994
995 int nb_row = row_data.getIndices().size();
996 int nb_col = col_data.getIndices().size();
997
998 int *row_indices_ptr = &row_data.getIndices()[0];
999 if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
1000 rowIndices.resize(nb_row, false);
1001 noalias(rowIndices) = row_data.getIndices();
1002 row_indices_ptr = &rowIndices[0];
1003 VectorDofs &dofs = row_data.getFieldDofs();
1004 VectorDofs::iterator dit = dofs.begin();
1005 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
1006 if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
1007 dAta.forcesOnlyOnEntitiesRow.end()) {
1008 rowIndices[ii] = -1;
1009 }
1010 }
1011 }
1012
1013 int *col_indices_ptr = &col_data.getIndices()[0];
1014 if (!dAta.forcesOnlyOnEntitiesCol.empty()) {
1015 colIndices.resize(nb_col, false);
1016 noalias(colIndices) = col_data.getIndices();
1017 col_indices_ptr = &colIndices[0];
1018 VectorDofs &dofs = col_data.getFieldDofs();
1019 VectorDofs::iterator dit = dofs.begin();
1020 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
1021 if (dAta.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
1022 dAta.forcesOnlyOnEntitiesCol.end()) {
1023 colIndices[ii] = -1;
1024 }
1025 }
1026 }
1027
1028 /*for(int dd1 = 0;dd1<k.size1();dd1++) {
1029 for(int dd2 = 0;dd2<k.size2();dd2++) {
1030 if(k(dd1,dd2)!=k(dd1,dd2)) {
1031 SETERRQ(PETSC_COMM_SELF,1,"Wrong result");
1032 }
1033 }
1034 }*/
1035
1036 CHKERR MatSetValues(getFEMethod()->snes_B, nb_row, row_indices_ptr, nb_col,
1037 col_indices_ptr, &k(0, 0), ADD_VALUES);
1038
1040}
1041
1043 const std::string field_name, BlockData &data, CommonData &common_data,
1044 int tag, bool jacobian, bool ale)
1045 : OpJacobianPiolaKirchhoffStress(field_name, data, common_data, tag,
1046 jacobian, ale, false) {}
1047
1050 const int gg) {
1052
1053 CHKERR dAta.materialAdoublePtr->calculatesIGma_EshelbyStress(
1054 dAta, getNumeredEntFiniteElementPtr());
1055 if (aLe) {
1056 auto &t_sIGma = dAta.materialAdoublePtr->t_sIGma;
1057 auto &t_invH = dAta.materialAdoublePtr->t_invH;
1058 t_sIGma(i, j) = t_sIGma(i, k) * t_invH(j, k);
1059 t_sIGma(i, j) *= dAta.materialAdoublePtr->detH;
1060
1061 }
1062 commonData.sTress[gg].resize(3, 3, false);
1063 for (int dd1 = 0; dd1 < 3; dd1++) {
1064 for (int dd2 = 0; dd2 < 3; dd2++) {
1065 dAta.materialAdoublePtr->sIGma(dd1, dd2) >>=
1066 (commonData.sTress[gg])(dd1, dd2);
1067 }
1068 }
1069
1071}
1072
1074 const std::string field_name, BlockData &data, CommonData &common_data)
1075 : OpRhsPiolaKirchhoff(field_name, data, common_data) {}
1076
1078 const std::string vel_field, const std::string field_name, BlockData &data,
1079 CommonData &common_data)
1080 : OpLhsPiolaKirchhoff_dX(vel_field, field_name, data, common_data) {}
1081
1083 EntitiesFieldData::EntData &col_data, int gg) {
1084 return get_jac<0>(col_data, gg, commonData.jacStress[gg], jac);
1085}
1086
1088 const std::string vel_field, const std::string field_name, BlockData &data,
1089 CommonData &common_data)
1090 : OpLhsPiolaKirchhoff_dx(vel_field, field_name, data, common_data) {}
1091
1093 EntitiesFieldData::EntData &col_data, int gg) {
1094 return get_jac<9>(col_data, gg, commonData.jacStress[gg], jac);
1095}
1096
1099 materialDoublePtr,
1101 materialAdoublePtr) {
1103
1104 if (!materialDoublePtr) {
1106 "Pointer for materialDoublePtr not allocated");
1107 }
1108 if (!materialAdoublePtr) {
1110 "Pointer for materialAdoublePtr not allocated");
1111 }
1112
1114 mField, BLOCKSET | MAT_ELASTICSET, it)) {
1115 Mat_Elastic mydata;
1116 CHKERR it->getAttributeDataStructure(mydata);
1117 int id = it->getMeshsetId();
1118 EntityHandle meshset = it->getMeshset();
1119 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET,
1120 setOfBlocks[id].tEts, true);
1121 setOfBlocks[id].iD = id;
1122 setOfBlocks[id].E = mydata.data.Young;
1123 setOfBlocks[id].PoissonRatio = mydata.data.Poisson;
1124 setOfBlocks[id].materialDoublePtr = materialDoublePtr;
1125 setOfBlocks[id].materialAdoublePtr = materialAdoublePtr;
1126 }
1127
1129}
1130
1132 const std::string element_name,
1133 const std::string spatial_position_field_name,
1134 const std::string material_position_field_name, const bool ale) {
1136
1137 CHKERR mField.add_finite_element(element_name, MF_ZERO);
1139 element_name, spatial_position_field_name);
1141 element_name, spatial_position_field_name);
1143 element_name, spatial_position_field_name);
1144 if (mField.check_field(material_position_field_name)) {
1145 if (ale) {
1147 element_name, material_position_field_name);
1149 element_name, material_position_field_name);
1150 }
1152 element_name, material_position_field_name);
1153 }
1154
1155 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
1156 for (; sit != setOfBlocks.end(); sit++) {
1157 CHKERR mField.add_ents_to_finite_element_by_type(sit->second.tEts, MBTET,
1158 element_name);
1159 }
1160
1162}
1163
1165 const std::string spatial_position_field_name,
1166 const std::string material_position_field_name, const bool ale,
1167 const bool field_disp) {
1169
1170 commonData.spatialPositions = spatial_position_field_name;
1171 commonData.meshPositions = material_position_field_name;
1172
1173 // Rhs
1174 feRhs.getOpPtrVector().push_back(
1175 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
1176 if (mField.check_field(material_position_field_name)) {
1178 material_position_field_name, commonData));
1179 }
1180 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
1181 for (; sit != setOfBlocks.end(); sit++) {
1183 spatial_position_field_name, sit->second, commonData, tAg, false, ale,
1184 field_disp));
1186 spatial_position_field_name, sit->second, commonData));
1187 }
1188
1189 // Energy
1190 feEnergy.getOpPtrVector().push_back(
1191 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
1192 if (mField.check_field(material_position_field_name)) {
1194 material_position_field_name, commonData));
1195 }
1196 sit = setOfBlocks.begin();
1197 for (; sit != setOfBlocks.end(); sit++) {
1198 feEnergy.getOpPtrVector().push_back(
1199 new OpEnergy(spatial_position_field_name, sit->second, commonData,
1200 feEnergy.V, field_disp));
1201 }
1202
1203 // Lhs
1204 feLhs.getOpPtrVector().push_back(
1205 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
1206 if (mField.check_field(material_position_field_name)) {
1208 material_position_field_name, commonData));
1209 }
1210 sit = setOfBlocks.begin();
1211 for (; sit != setOfBlocks.end(); sit++) {
1213 spatial_position_field_name, sit->second, commonData, tAg, true, ale,
1214 field_disp));
1216 spatial_position_field_name, spatial_position_field_name, sit->second,
1217 commonData));
1218 }
1219
1221}
static Index< 'n', 3 > n
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:111
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:172
@ BLOCKSET
Definition: definitions.h:161
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:47
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
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 add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
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_col(const std::string &fe_name, const std::string &name_row)=0
set field col 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.
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:126
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
Definition: Types.hpp:143
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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:1218
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1207
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
const double r
rate factor
double A
double H
Definition: plastic.cpp:109
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 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.
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
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.