v0.14.0
Loading...
Searching...
No Matches
KelvinVoigtDamper.hpp
Go to the documentation of this file.
1/** \file KelvinVoigtDamper.hpp
2 * \brief Implementation dashpot, i.e. damper
3 * \ingroup nonlinear_elastic_elem
4 *
5 */
6
7
8
9#ifndef __KELVIN_VOIGT_DAMPER_HPP__
10#define __KELVIN_VOIGT_DAMPER_HPP__
11
12#ifndef WITH_ADOL_C
13#error "MoFEM need to be compiled with ADOL-C"
14#endif
15
16/** \brief Implementation of Kelvin Voigt Damper
17\ingroup nonlinear_elastic_elem
18
19*/
21
23
25
26 /** \brief Dumper material parameters
27 \ingroup nonlinear_elastic_elem
28 */
31 double vBeta; ///< Poisson ration spring alpha
32 double gBeta; ///< Sheer modulus spring alpha
33 bool lInear;
34 BlockMaterialData() : vBeta(0), gBeta(1), lInear(false) {}
35 };
36
37 std::map<int, BlockMaterialData> blockMaterialDataMap;
38
39 /** \brief Constitutive model functions
40 \ingroup nonlinear_elastic_elem
41
42 */
43 template <typename TYPE> struct ConstitutiveEquation {
44
48
51
52 ConstitutiveEquation(BlockMaterialData &data, bool is_displacement = true)
53 : dAta(data), isDisplacement(is_displacement) {}
54 virtual ~ConstitutiveEquation() = default;
55
56 MatrixBoundedArray<TYPE, 9> F; ///< Gradient of deformation
57 MatrixBoundedArray<TYPE, 9> FDot; ///< Rate of gradient of deformation
58 MatrixBoundedArray<TYPE, 9>
59 gradientUDot; ///< Rate of gradient of displacements
60 MatrixBoundedArray<TYPE, 9> engineringStrainDot;
61 MatrixBoundedArray<TYPE, 9>
62 dashpotCauchyStress; ///< Stress generated by spring beta
63 MatrixBoundedArray<TYPE, 9>
64 dashpotFirstPiolaKirchhoffStress; ///< Stress generated by spring beta
65 MatrixBoundedArray<TYPE, 9> invF; ///< Inverse of gradient of deformation
66
68 TYPE J; ///< Jacobian of gradient of deformation
69
70 /** \brief Calculate strain rate
71
72 \f[
73 \dot{\varepsilon}_{ij} = \frac{1}{2}
74 \left(
75 \frac{\partial v_i}{\partial X_j}
76 +
77 \frac{\partial v_j}{\partial X_i}
78 \right)
79 \f]
80
81 */
82 virtual MoFEMErrorCode calculateEngineeringStrainDot() {
84 gradientUDot.resize(3, 3, false);
85 noalias(gradientUDot) = FDot;
86
87 // for (int ii = 0; ii < 3; ii++)
88 // gradientUDot(ii, ii) -= 1;
89
91 for (int ii = 0; ii < 3; ii++) {
93 }
94 engineringStrainDot.resize(3, 3, false);
98 }
99
100 /** \brief Calculate Cauchy dashpot stress
101
102 Calculate dashpot Cauchy stress. It has to be pull back to reference
103 configuration before use in total Lagrangian formulation.
104
105 \f[
106 \sigma^\beta_{ij} = 2G^\beta\left[
107 \dot{\varepsilon}_{ij}
108 + \frac{v^\beta}{1-2v^\beta}\dot{\varepsilon}_{kk}\delta_{ij}
109 \right]
110 \f]
111
112 */
113 virtual MoFEMErrorCode calculateDashpotCauchyStress() {
115 dashpotCauchyStress.resize(3, 3, false);
116 double a = 2.0 * dAta.gBeta;
117 double b = a * (dAta.vBeta / (1.0 - 2.0 * dAta.vBeta));
119 for (int ii = 0; ii < 3; ii++) {
121 }
123 }
124
125 /** \brief Calculate First Piola-Kirchhoff Stress Dashpot stress
126
127 \f[
128 P^\beta_{ij} = J \sigma^\beta_{ik} F^{-1}_{jk}
129 \f]
130
131 */
132 virtual MoFEMErrorCode calculateFirstPiolaKirchhoffStress() {
134 dashpotFirstPiolaKirchhoffStress.resize(3, 3, false);
135 if (dAta.lInear) {
137 } else {
138
139 invF.resize(3, 3, false);
140
141 auto t_dashpotFirstPiolaKirchhoffStress = getFTensor2FromArray3by3(
143 auto t_dashpotCauchyStress = getFTensor2FromArray3by3(
145 auto t_F = getFTensor2FromArray3by3(F, FTensor::Number<0>(), 0);
146 auto t_invF = getFTensor2FromArray3by3(invF, FTensor::Number<0>(), 0);
147 if (isDisplacement) {
148 t_F(0, 0) += 1;
149 t_F(1, 1) += 1;
150 t_F(2, 2) += 1;
151 }
152
153 J = determinantTensor3by3(t_F);
154 CHKERR invertTensor3by3(t_F, J, t_invF);
155 t_dashpotFirstPiolaKirchhoffStress(i, j) =
156 J * (t_dashpotCauchyStress(i, k) * t_invF(j, k));
157 }
159 }
160 };
161
162 typedef boost::ptr_map<int, KelvinVoigtDamper::ConstitutiveEquation<adouble>>
165
166 /** \brief Common data for nonlinear_elastic_elem model
167 \ingroup nonlinear_elastic_elem
168 */
169 struct CommonData {
170
174
175 std::map<std::string, std::vector<VectorDouble>> dataAtGaussPts;
176 std::map<std::string, std::vector<MatrixDouble>> gradAtGaussPts;
177
178 boost::shared_ptr<MatrixDouble> dataAtGaussTmpPtr;
179 boost::shared_ptr<MatrixDouble> gradDataAtGaussTmpPtr;
180
181 std::vector<MatrixDouble> dashpotFirstPiolaKirchhoffStress;
182
183 std::vector<double *> jacRowPtr;
184 std::vector<MatrixDouble> jacStress;
185
189
190 CommonData() : recordOn(true), skipThis(true) {
191 dataAtGaussTmpPtr = boost::make_shared<MatrixDouble>();
192 dataAtGaussTmpPtr->resize(3, 1);
193 gradDataAtGaussTmpPtr = boost::make_shared<MatrixDouble>();
194 gradDataAtGaussTmpPtr->resize(9, 1);
195 meshNodePositionName = "MESH_NODE_POSITIONS";
196 }
197 };
199
200 /// \brief definition of volume element
202
204 int addToRule; ///< Takes into account HO geometry
205
206 DamperFE(MoFEM::Interface &m_field, CommonData &common_data)
208 commonData(common_data), addToRule(1) {}
209
210 int getRule(int order) { return order + addToRule; }
211
212 MoFEMErrorCode preProcess() {
213
216
217 // if (ts_ctx == CTX_TSSETIFUNCTION) {
218
219 // CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
220 // problemPtr, commonData.spatialPositionName,
221 // commonData.spatialPositionNameDot, COL, ts_u_t, INSERT_VALUES,
222 // SCATTER_REVERSE);
223 // }
224
226 }
227
228 MoFEMErrorCode postProcess() {
229
231
232 // CHKERR MoFEM::VolumeElementForcesAndSourcesCore::postProcess();
233
234 // if (ts_ctx == CTX_TSSETIFUNCTION) {
235 // CHKERR VecAssemblyBegin(ts_F);
236 // CHKERR VecAssemblyEnd(ts_F);
237 // }
238 // if (ts_ctx == CTX_TSSETIJACOBIAN) {
239 // CHKERR MatAssemblyBegin(ts_B, MAT_FLUSH_ASSEMBLY);
240 // CHKERR MatAssemblyEnd(ts_B, MAT_FLUSH_ASSEMBLY);
241 // }
242
244 }
245 };
246
248
250 : mField(m_field), feRhs(m_field, commonData),
251 feLhs(m_field, commonData) {}
252
255
261
262 OpGetDataAtGaussPts(const std::string field_name, CommonData &common_data,
263 bool calc_val, bool calc_grad, bool calc_dot = false,
264 EntityType zero_at_type = MBVERTEX)
265 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
267 commonData(common_data), calcVal(calc_val), calcGrad(calc_grad),
268 calcDot(calc_dot), zeroAtType(zero_at_type) {}
269
270 /** \brief Operator field value
271 *
272 */
273 MoFEMErrorCode doWork(int side, EntityType type,
274 EntitiesFieldData::EntData &data) {
276
277 int nb_dofs = data.getFieldData().size();
278 if (nb_dofs == 0) {
280 }
281 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
282 int nb_gauss_pts = data.getN().size1();
283
285 if (calcDot)
287 // Initialize
288 if (calcVal) {
289 commonData.dataAtGaussPts[field_name].resize(nb_gauss_pts);
290 for (int gg = 0; gg < nb_gauss_pts; gg++) {
291 commonData.dataAtGaussPts[field_name][gg].resize(rank, false);
293 }
294 }
295 if (calcGrad) {
296 commonData.gradAtGaussPts[field_name].resize(nb_gauss_pts);
297 for (int gg = 0; gg < nb_gauss_pts; gg++) {
298 commonData.gradAtGaussPts[field_name][gg].resize(rank, 3, false);
300 }
301 }
302
303 // Zero values
304 // if (type == zeroAtType) {
305 // for (int gg = 0; gg < nb_gauss_pts; gg++) {
306 // if (calcVal) {
307 // commonData.dataAtGaussPts[rowFieldName][gg].clear();
308 // }
309 // if (calcGrad) {
310 // commonData.gradAtGaussPts[rowFieldName][gg].clear();
311 // }
312 // }
313 // }
314
315 auto t_disp = getFTensor1FromMat<3>(*commonData.dataAtGaussTmpPtr);
316 auto t_diff_disp =
317 getFTensor2FromMat<3, 3>(*commonData.gradDataAtGaussTmpPtr);
318
319 // Calculate values at integration points
320 for (int gg = 0; gg < nb_gauss_pts; gg++) {
321 for (int rr1 = 0; rr1 < rank; rr1++) {
322 if (calcVal) {
323 commonData.dataAtGaussPts[field_name][gg][rr1] = t_disp(rr1);
324 }
325 if (calcGrad) {
326 for (int rr2 = 0; rr2 < 3; rr2++) {
327 commonData.gradAtGaussPts[field_name][gg](rr1, rr2) =
328 t_diff_disp(rr1, rr2);
329 }
330 }
331 }
332 ++t_disp;
333 ++t_diff_disp;
334 }
335
337 }
338 };
339
342
343 std::vector<int> tagS;
346
349 bool &recordOn;
350 std::map<int, int> &nbActiveVariables;
351 std::map<int, int> &nbActiveResults;
352
353 OpJacobian(const std::string field_name, std::vector<int> tags,
355 CommonData &common_data, bool calculate_residual,
356 bool calculate_jacobian)
357 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
359 tagS(tags), cE(ce), commonData(common_data),
360 calculateResidualBool(calculate_residual),
361 calculateJacobianBool(calculate_jacobian),
362 recordOn(common_data.recordOn),
364 nbActiveResults(common_data.nbActiveResults) {}
365
367 VectorDouble activeVariables;
368
369 MoFEMErrorCode recordDamperStress() {
371
372 if (tagS[DAMPERSTRESS] < 0) {
374 }
375
376 cE.F.resize(3, 3, false);
377 cE.FDot.resize(3, 3, false);
378 MatrixDouble &F =
380 MatrixDouble &F_dot =
382 trace_on(tagS[DAMPERSTRESS]);
383 {
384 // Activate gradient of defamation
386 for (int dd1 = 0; dd1 < 3; dd1++) {
387 for (int dd2 = 0; dd2 < 3; dd2++) {
388 cE.F(dd1, dd2) <<= F(dd1, dd2);
390 }
391 }
392 for (int dd1 = 0; dd1 < 3; dd1++) {
393 for (int dd2 = 0; dd2 < 3; dd2++) {
394 cE.FDot(dd1, dd2) <<= F_dot(dd1, dd2);
396 }
397 }
398
399 // Do calculations
403
404 // Results
407 commonData.dashpotFirstPiolaKirchhoffStress[0].resize(3, 3, false);
408 for (int d1 = 0; d1 < 3; d1++) {
409 for (int d2 = 0; d2 < 3; d2++) {
413 }
414 }
415 }
416 trace_off();
418 }
419
420 MoFEMErrorCode calculateFunction(TagEvaluate te, double *ptr) {
422
423 int r;
424 // play recorder for values
425 r = ::function(tagS[te], nbActiveResults[tagS[te]],
426 nbActiveVariables[tagS[te]], &activeVariables[0], ptr);
427 if (r < 3) { // function is locally analytic
428 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
429 "ADOL-C function evaluation with error r = %d", r);
430 }
431
433 }
434
435 MoFEMErrorCode calculateJacobian(TagEvaluate te) {
437
438 try {
439 int r;
440 r = jacobian(tagS[te], nbActiveResults[tagS[te]],
442 &(commonData.jacRowPtr[0]));
443 if (r < 3) {
444 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
445 "ADOL-C function evaluation with error");
446 }
447 } catch (const std::exception &ex) {
448 std::ostringstream ss;
449 ss << "throw in method: " << ex.what() << std::endl;
450 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
451 }
453 }
454
457
458 if (tagS[DAMPERSTRESS] < 0) {
460 }
461
463 for (int gg = 0; gg < nbGaussPts; gg++) {
464
465 MatrixDouble &F =
467 MatrixDouble &F_dot =
469 int nb_active_variables = 0;
470
471 // Activate gradient of defamation
472 for (int dd1 = 0; dd1 < 3; dd1++) {
473 for (int dd2 = 0; dd2 < 3; dd2++) {
474 activeVariables[nb_active_variables++] = F(dd1, dd2);
475 }
476 }
477 // Activate rate of gradient of defamation
478 for (int dd1 = 0; dd1 < 3; dd1++) {
479 for (int dd2 = 0; dd2 < 3; dd2++) {
480 activeVariables[nb_active_variables++] = F_dot(dd1, dd2);
481 }
482 }
483
484 if (nb_active_variables != nbActiveVariables[tagS[DAMPERSTRESS]]) {
485 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
486 "Number of active variables does not much");
487 }
488
490 if (gg == 0) {
492 }
493 commonData.dashpotFirstPiolaKirchhoffStress[gg].resize(3, 3, false);
497 }
498
500 if (gg == 0) {
503 }
506 false);
507 for (int dd = 0; dd < nbActiveResults[tagS[DAMPERSTRESS]]; dd++) {
508 commonData.jacRowPtr[dd] = &commonData.jacStress[gg](dd, 0);
509 }
511 }
512 }
513
515 }
516
517 MoFEMErrorCode doWork(int row_side, EntityType row_type,
518 EntitiesFieldData::EntData &row_data) {
520
521 if (row_type != MBVERTEX)
523 nbGaussPts = row_data.getN().size1();
524
525 commonData.skipThis = false;
526 if (cE.dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
527 cE.dAta.tEts.end()) {
528 commonData.skipThis = true;
530 }
531
532 if (recordOn) {
534 }
536
538 }
539 };
540
544 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
546
547 VectorDouble nF;
548 MoFEMErrorCode aSemble(int row_side, EntityType row_type,
549 EntitiesFieldData::EntData &row_data) {
551 int nb_dofs = row_data.getIndices().size();
552 int *indices_ptr = &row_data.getIndices()[0];
553 CHKERR VecSetValues(getFEMethod()->ts_F, nb_dofs, indices_ptr, &nF[0],
554 ADD_VALUES);
556 }
557 };
558
559 /** \brief Assemble internal force vector
560 \ingroup nonlinear_elastic_elem
561
562 */
563 struct OpRhsStress : public AssembleVector {
566 : AssembleVector(common_data.spatialPositionName),
567 commonData(common_data) {}
568 MoFEMErrorCode doWork(int row_side, EntityType row_type,
569 EntitiesFieldData::EntData &row_data) {
571
572 if (commonData.skipThis) {
574 }
575
576 int nb_dofs = row_data.getIndices().size();
577 if (!nb_dofs) {
579 }
580 nF.resize(nb_dofs, false);
581 nF.clear();
582 int nb_gauss_pts = row_data.getN().size1();
583 for (int gg = 0; gg != nb_gauss_pts; gg++) {
584 const MatrixAdaptor &diffN = row_data.getDiffN(gg, nb_dofs / 3);
585 const MatrixDouble &stress =
587 double val = getVolume() * getGaussPts()(3, gg);
588 for (int dd = 0; dd < nb_dofs / 3; dd++) {
589 for (int rr = 0; rr < 3; rr++) {
590 for (int nn = 0; nn < 3; nn++) {
591 nF[3 * dd + rr] += val * diffN(dd, nn) * stress(rr, nn);
592 }
593 }
594 }
595 }
596 CHKERR aSemble(row_side, row_type, row_data);
598 }
599 };
600
603 AssembleMatrix(string row_name, string col_name)
604 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
605 row_name, col_name, UserDataOperator::OPROWCOL) {}
606
607 MatrixDouble K, transK;
608 MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type,
609 EntityType col_type,
610 EntitiesFieldData::EntData &row_data,
611 EntitiesFieldData::EntData &col_data) {
613 int nb_row = row_data.getIndices().size();
614 int nb_col = col_data.getIndices().size();
615 int *row_indices_ptr = &row_data.getIndices()[0];
616 int *col_indices_ptr = &col_data.getIndices()[0];
617 CHKERR MatSetValues(getFEMethod()->ts_B, nb_row, row_indices_ptr, nb_col,
618 col_indices_ptr, &K(0, 0), ADD_VALUES);
619 if (sYmm) {
620 // Assemble of diagonal terms
621 if (row_side != col_side || row_type != col_type) {
622 transK.resize(nb_col, nb_row, false);
623 noalias(transK) = trans(K);
624 CHKERR MatSetValues(getFEMethod()->ts_B, nb_col, col_indices_ptr,
625 nb_row, row_indices_ptr, &transK(0, 0),
626 ADD_VALUES);
627 }
628 }
630 }
631 };
632
633 /** \brief Assemble matrix
634 */
635 struct OpLhsdxdx : public AssembleMatrix {
637 OpLhsdxdx(CommonData &common_data)
638 : AssembleMatrix(common_data.spatialPositionName,
639 common_data.spatialPositionName),
640 commonData(common_data) {}
641 MatrixDouble dStress_dx;
642 MoFEMErrorCode get_dStress_dx(EntitiesFieldData::EntData &col_data,
643 int gg) {
645 int nb_col = col_data.getIndices().size();
646 dStress_dx.resize(9, nb_col, false);
647 dStress_dx.clear();
648 const MatrixAdaptor diffN = col_data.getDiffN(gg, nb_col / 3);
649 MatrixDouble &jac_stress = commonData.jacStress[gg];
650 for (int dd = 0; dd < nb_col / 3; dd++) { // DoFs in column
651 for (int jj = 0; jj < 3; jj++) { // cont. DoFs in column
652 double a = diffN(dd, jj);
653 for (int rr = 0; rr < 3; rr++) { // Loop over dsigma_ii/dX_rr
654 for (int ii = 0; ii < 9;
655 ii++) { // ii represents components of stress tensor
656 dStress_dx(ii, 3 * dd + rr) += jac_stress(ii, 3 * rr + jj) * a;
657 }
658 }
659 }
660 }
662 }
663 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
664 EntityType col_type,
665 EntitiesFieldData::EntData &row_data,
666 EntitiesFieldData::EntData &col_data) {
668
669 if (commonData.skipThis) {
671 }
672
673 int nb_row = row_data.getIndices().size();
674 int nb_col = col_data.getIndices().size();
675 if (nb_row == 0)
677 if (nb_col == 0)
679 K.resize(nb_row, nb_col, false);
680 K.clear();
681 int nb_gauss_pts = row_data.getN().size1();
682 for (int gg = 0; gg != nb_gauss_pts; gg++) {
683 CHKERR get_dStress_dx(col_data, gg);
684 double val = getVolume() * getGaussPts()(3, gg);
685 // std::cerr << dStress_dx << std::endl;
686 dStress_dx *= val;
687 const MatrixAdaptor &diffN = row_data.getDiffN(gg, nb_row / 3);
688 { // integrate element stiffness matrix
689 for (int dd1 = 0; dd1 < nb_row / 3; dd1++) {
690 for (int rr1 = 0; rr1 < 3; rr1++) {
691 for (int dd2 = 0; dd2 < nb_col / 3; dd2++) {
692 for (int rr2 = 0; rr2 < 3; rr2++) {
693 K(3 * dd1 + rr1, 3 * dd2 + rr2) +=
694 (diffN(dd1, 0) * dStress_dx(3 * rr1 + 0, 3 * dd2 + rr2) +
695 diffN(dd1, 1) * dStress_dx(3 * rr1 + 1, 3 * dd2 + rr2) +
696 diffN(dd1, 2) * dStress_dx(3 * rr1 + 2, 3 * dd2 + rr2));
697 }
698 }
699 }
700 }
701 }
702 }
703 // std::cerr << "G " << getNumeredEntFiniteElementPtr()->getRefEnt() <<
704 // std::endl << K << std::endl;
705 CHKERR aSemble(row_side, col_side, row_type, col_type, row_data,
706 col_data);
707
709 }
710 };
711
712 /** \brief Assemble matrix
713 */
714 struct OpLhsdxdot : public AssembleMatrix {
716 OpLhsdxdot(CommonData &common_data)
717 : AssembleMatrix(common_data.spatialPositionName,
718 common_data.spatialPositionName),
719 commonData(common_data) {}
720 MatrixDouble dStress_dot;
721 MoFEMErrorCode get_dStress_dot(EntitiesFieldData::EntData &col_data,
722 int gg) {
724 int nb_col = col_data.getIndices().size();
725 dStress_dot.resize(9, nb_col, false);
726 dStress_dot.clear();
727 const MatrixAdaptor diffN = col_data.getDiffN(gg, nb_col / 3);
728 MatrixDouble &jac_stress = commonData.jacStress[gg];
729 for (int dd = 0; dd < nb_col / 3; dd++) { // DoFs in column
730 for (int jj = 0; jj < 3; jj++) { // cont. DoFs in column
731 double a = diffN(dd, jj);
732 for (int rr = 0; rr < 3; rr++) { // Loop over dsigma_ii/dX_rr
733 for (int ii = 0; ii < 9;
734 ii++) { // ii represents components of stress tensor
735 dStress_dot(ii, 3 * dd + rr) +=
736 jac_stress(ii, 9 + 3 * rr + jj) * a * getFEMethod()->ts_a;
737 }
738 }
739 }
740 }
742 }
743 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
744 EntityType col_type,
745 EntitiesFieldData::EntData &row_data,
746 EntitiesFieldData::EntData &col_data) {
748
749 if (commonData.skipThis) {
751 }
752
753 int nb_row = row_data.getIndices().size();
754 int nb_col = col_data.getIndices().size();
755 if (nb_row == 0)
757 if (nb_col == 0)
759 K.resize(nb_row, nb_col, false);
760 K.clear();
761 int nb_gauss_pts = row_data.getN().size1();
762 for (int gg = 0; gg != nb_gauss_pts; gg++) {
763 CHKERR get_dStress_dot(col_data, gg);
764 double val = getVolume() * getGaussPts()(3, gg);
765 // std::cerr << dStress_dot << std::endl;
766 dStress_dot *= val;
767 const MatrixAdaptor &diffN = row_data.getDiffN(gg, nb_row / 3);
768 { // integrate element stiffness matrix
769 for (int dd1 = 0; dd1 < nb_row / 3; dd1++) {
770 for (int rr1 = 0; rr1 < 3; rr1++) {
771 for (int dd2 = 0; dd2 < nb_col / 3; dd2++) {
772 for (int rr2 = 0; rr2 < 3; rr2++) {
773 K(3 * dd1 + rr1, 3 * dd2 + rr2) +=
774 (diffN(dd1, 0) * dStress_dot(3 * rr1 + 0, 3 * dd2 + rr2) +
775 diffN(dd1, 1) * dStress_dot(3 * rr1 + 1, 3 * dd2 + rr2) +
776 diffN(dd1, 2) * dStress_dot(3 * rr1 + 2, 3 * dd2 + rr2));
777 }
778 }
779 }
780 }
781 }
782 }
783 // std::cerr << "G " << getNumeredEntFiniteElementPtr()->getRefEnt() <<
784 // std::endl << K << std::endl;
785 CHKERR aSemble(row_side, col_side, row_type, col_type, row_data,
786 col_data);
787
789 }
790 };
791
792 MoFEMErrorCode setBlockDataMap() {
794
796 if (it->getName().compare(0, 6, "DAMPER") == 0) {
797 std::vector<double> data;
798 CHKERR it->getAttributes(data);
799 if (data.size() < 2) {
800 SETERRQ(PETSC_COMM_SELF, 1, "Data inconsistency");
801 }
802 CHKERR mField.get_moab().get_entities_by_type(
803 it->meshset, MBTET, blockMaterialDataMap[it->getMeshsetId()].tEts,
804 true);
805 blockMaterialDataMap[it->getMeshsetId()].gBeta = data[0];
806 blockMaterialDataMap[it->getMeshsetId()].vBeta = data[1];
807 }
808 }
810 }
811
812 MoFEMErrorCode setOperators(const int tag) {
814
815 for (auto &&fe_ptr : {&feRhs, &feLhs}) {
816 // CHKERR addHOOpsVol(commonData.meshNodePositionName, *fe_ptr, true, false,
817 // false, false);
818 fe_ptr->getOpPtrVector().push_back(
822 fe_ptr->getOpPtrVector().push_back(new OpGetDataAtGaussPts(
823 commonData.spatialPositionName, commonData, false, true, false));
824 fe_ptr->getOpPtrVector().push_back(
825 new OpCalculateVectorFieldGradientDot<3, 3>(
828 fe_ptr->getOpPtrVector().push_back(new OpGetDataAtGaussPts(
829 commonData.spatialPositionName, commonData, false, true, true));
830 }
831
832 // attach tags for each recorder
833 std::vector<int> tags;
834 tags.push_back(tag);
835
836 ConstitutiveEquationMap::iterator mit = constitutiveEquationMap.begin();
837 for (; mit != constitutiveEquationMap.end(); mit++) {
839 constitutiveEquationMap.at(mit->first);
840 // Right hand side operators
841 feRhs.getOpPtrVector().push_back(new OpJacobian(
842 commonData.spatialPositionName, tags, ce, commonData, true, false));
843 feRhs.getOpPtrVector().push_back(new OpRhsStress(commonData));
844
845 // Left hand side operators
846 feLhs.getOpPtrVector().push_back(new OpJacobian(
847 commonData.spatialPositionName, tags, ce, commonData, false, true));
848 feLhs.getOpPtrVector().push_back(new OpLhsdxdx(commonData));
849 feLhs.getOpPtrVector().push_back(new OpLhsdxdot(commonData));
850 }
851
853 }
854};
855
856#endif //__KELVIN_VOIGT_DAMPER_HPP__
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
#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
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_IMPOSSIBLE_CASE
Definition: definitions.h:35
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
#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
constexpr int order
@ F
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
constexpr auto field_name
MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
AssembleMatrix(string row_name, string col_name)
MoFEMErrorCode aSemble(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
double vBeta
Poisson ration spring alpha.
double gBeta
Sheer modulus spring alpha.
Common data for nonlinear_elastic_elem model.
std::map< int, int > nbActiveVariables
boost::shared_ptr< MatrixDouble > dataAtGaussTmpPtr
std::vector< MatrixDouble > dashpotFirstPiolaKirchhoffStress
std::vector< double * > jacRowPtr
std::vector< MatrixDouble > jacStress
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
std::map< std::string, std::vector< VectorDouble > > dataAtGaussPts
std::map< int, int > nbActiveResults
boost::shared_ptr< MatrixDouble > gradDataAtGaussTmpPtr
MatrixBoundedArray< TYPE, 9 > FDot
Rate of gradient of deformation.
MatrixBoundedArray< TYPE, 9 > invF
Inverse of gradient of deformation.
MatrixBoundedArray< TYPE, 9 > F
Gradient of deformation.
virtual MoFEMErrorCode calculateFirstPiolaKirchhoffStress()
Calculate First Piola-Kirchhoff Stress Dashpot stress.
MatrixBoundedArray< TYPE, 9 > gradientUDot
Rate of gradient of displacements.
virtual MoFEMErrorCode calculateDashpotCauchyStress()
Calculate Cauchy dashpot stress.
MatrixBoundedArray< TYPE, 9 > engineringStrainDot
virtual MoFEMErrorCode calculateEngineeringStrainDot()
Calculate strain rate.
MatrixBoundedArray< TYPE, 9 > dashpotFirstPiolaKirchhoffStress
Stress generated by spring beta.
ConstitutiveEquation(BlockMaterialData &data, bool is_displacement=true)
MatrixBoundedArray< TYPE, 9 > dashpotCauchyStress
Stress generated by spring beta.
TYPE J
Jacobian of gradient of deformation.
definition of volume element
int addToRule
Takes into account HO geometry.
MoFEMErrorCode postProcess()
function is run at the end of loop
DamperFE(MoFEM::Interface &m_field, CommonData &common_data)
MoFEMErrorCode preProcess()
function is run at the beginning of loop
OpGetDataAtGaussPts(const std::string field_name, CommonData &common_data, bool calc_val, bool calc_grad, bool calc_dot=false, EntityType zero_at_type=MBVERTEX)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator field value.
std::map< int, int > & nbActiveVariables
MoFEMErrorCode calculateJacobian(TagEvaluate te)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
std::map< int, int > & nbActiveResults
KelvinVoigtDamper::ConstitutiveEquation< adouble > & cE
MoFEMErrorCode calculateAtIntPtsDamperStress()
MoFEMErrorCode calculateFunction(TagEvaluate te, double *ptr)
OpJacobian(const std::string field_name, std::vector< int > tags, KelvinVoigtDamper::ConstitutiveEquation< adouble > &ce, CommonData &common_data, bool calculate_residual, bool calculate_jacobian)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
MoFEMErrorCode get_dStress_dot(EntitiesFieldData::EntData &col_data, int gg)
OpLhsdxdot(CommonData &common_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpLhsdxdx(CommonData &common_data)
MoFEMErrorCode get_dStress_dx(EntitiesFieldData::EntData &col_data, int gg)
Assemble internal force vector.
OpRhsStress(CommonData &common_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Implementation of Kelvin Voigt Damper.
MoFEMErrorCode setBlockDataMap()
ConstitutiveEquationMap constitutiveEquationMap
std::map< int, BlockMaterialData > blockMaterialDataMap
KelvinVoigtDamper(MoFEM::Interface &m_field)
MoFEMErrorCode setOperators(const int tag)
boost::ptr_map< int, KelvinVoigtDamper::ConstitutiveEquation< adouble > > ConstitutiveEquationMap
MoFEM::Interface & mField
virtual moab::Interface & get_moab()=0
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
@ OPCOL
operator doWork function is executed on FE columns
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
virtual MoFEMErrorCode preProcess()
function is run at the beginning of loop
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
VolumeElementForcesAndSourcesCore(Interface &m_field, const EntityType type=MBTET)