v0.15.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)
207 : MoFEM::VolumeElementForcesAndSourcesCore(m_field),
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
260 EntityType zeroAtType;
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 SETERRQ(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, MOFEM_DATA_INCONSISTENCY, "%s",
451 ss.str().c_str());
452 }
454 }
455
458
459 if (tagS[DAMPERSTRESS] < 0) {
461 }
462
464 for (int gg = 0; gg < nbGaussPts; gg++) {
465
466 MatrixDouble &F =
468 MatrixDouble &F_dot =
470 int nb_active_variables = 0;
471
472 // Activate gradient of defamation
473 for (int dd1 = 0; dd1 < 3; dd1++) {
474 for (int dd2 = 0; dd2 < 3; dd2++) {
475 activeVariables[nb_active_variables++] = F(dd1, dd2);
476 }
477 }
478 // Activate rate of gradient of defamation
479 for (int dd1 = 0; dd1 < 3; dd1++) {
480 for (int dd2 = 0; dd2 < 3; dd2++) {
481 activeVariables[nb_active_variables++] = F_dot(dd1, dd2);
482 }
483 }
484
485 if (nb_active_variables != nbActiveVariables[tagS[DAMPERSTRESS]]) {
486 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
487 "Number of active variables does not much");
488 }
489
491 if (gg == 0) {
493 }
494 commonData.dashpotFirstPiolaKirchhoffStress[gg].resize(3, 3, false);
498 }
499
501 if (gg == 0) {
504 }
507 false);
508 for (int dd = 0; dd < nbActiveResults[tagS[DAMPERSTRESS]]; dd++) {
509 commonData.jacRowPtr[dd] = &commonData.jacStress[gg](dd, 0);
510 }
512 }
513 }
514
516 }
517
518 MoFEMErrorCode doWork(int row_side, EntityType row_type,
519 EntitiesFieldData::EntData &row_data) {
521
522 if (row_type != MBVERTEX)
524 nbGaussPts = row_data.getN().size1();
525
526 commonData.skipThis = false;
527 if (cE.dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
528 cE.dAta.tEts.end()) {
529 commonData.skipThis = true;
531 }
532
533 if (recordOn) {
535 }
537
539 }
540 };
541
545 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
547
548 VectorDouble nF;
549 MoFEMErrorCode aSemble(int row_side, EntityType row_type,
550 EntitiesFieldData::EntData &row_data) {
552 int nb_dofs = row_data.getIndices().size();
553 int *indices_ptr = &row_data.getIndices()[0];
554 CHKERR VecSetValues(getFEMethod()->ts_F, nb_dofs, indices_ptr, &nF[0],
555 ADD_VALUES);
557 }
558 };
559
560 /** \brief Assemble internal force vector
561 \ingroup nonlinear_elastic_elem
562
563 */
564 struct OpRhsStress : public AssembleVector {
567 : AssembleVector(common_data.spatialPositionName),
568 commonData(common_data) {}
569 MoFEMErrorCode doWork(int row_side, EntityType row_type,
570 EntitiesFieldData::EntData &row_data) {
572
573 if (commonData.skipThis) {
575 }
576
577 int nb_dofs = row_data.getIndices().size();
578 if (!nb_dofs) {
580 }
581 nF.resize(nb_dofs, false);
582 nF.clear();
583 int nb_gauss_pts = row_data.getN().size1();
584 for (int gg = 0; gg != nb_gauss_pts; gg++) {
585 const MatrixAdaptor &diffN = row_data.getDiffN(gg, nb_dofs / 3);
586 const MatrixDouble &stress =
588 double val = getVolume() * getGaussPts()(3, gg);
589 for (int dd = 0; dd < nb_dofs / 3; dd++) {
590 for (int rr = 0; rr < 3; rr++) {
591 for (int nn = 0; nn < 3; nn++) {
592 nF[3 * dd + rr] += val * diffN(dd, nn) * stress(rr, nn);
593 }
594 }
595 }
596 }
597 CHKERR aSemble(row_side, row_type, row_data);
599 }
600 };
601
604 AssembleMatrix(string row_name, string col_name)
605 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
606 row_name, col_name, UserDataOperator::OPROWCOL) {}
607
608 MatrixDouble K, transK;
609 MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type,
610 EntityType col_type,
611 EntitiesFieldData::EntData &row_data,
612 EntitiesFieldData::EntData &col_data) {
614 int nb_row = row_data.getIndices().size();
615 int nb_col = col_data.getIndices().size();
616 int *row_indices_ptr = &row_data.getIndices()[0];
617 int *col_indices_ptr = &col_data.getIndices()[0];
618 CHKERR MatSetValues(getFEMethod()->ts_B, nb_row, row_indices_ptr, nb_col,
619 col_indices_ptr, &K(0, 0), ADD_VALUES);
620 if (sYmm) {
621 // Assemble of diagonal terms
622 if (row_side != col_side || row_type != col_type) {
623 transK.resize(nb_col, nb_row, false);
624 noalias(transK) = trans(K);
625 CHKERR MatSetValues(getFEMethod()->ts_B, nb_col, col_indices_ptr,
626 nb_row, row_indices_ptr, &transK(0, 0),
627 ADD_VALUES);
628 }
629 }
631 }
632 };
633
634 /** \brief Assemble matrix
635 */
636 struct OpLhsdxdx : public AssembleMatrix {
638 OpLhsdxdx(CommonData &common_data)
639 : AssembleMatrix(common_data.spatialPositionName,
640 common_data.spatialPositionName),
641 commonData(common_data) {}
642 MatrixDouble dStress_dx;
643 MoFEMErrorCode get_dStress_dx(EntitiesFieldData::EntData &col_data,
644 int gg) {
646 int nb_col = col_data.getIndices().size();
647 dStress_dx.resize(9, nb_col, false);
648 dStress_dx.clear();
649 const MatrixAdaptor diffN = col_data.getDiffN(gg, nb_col / 3);
650 MatrixDouble &jac_stress = commonData.jacStress[gg];
651 for (int dd = 0; dd < nb_col / 3; dd++) { // DoFs in column
652 for (int jj = 0; jj < 3; jj++) { // cont. DoFs in column
653 double a = diffN(dd, jj);
654 for (int rr = 0; rr < 3; rr++) { // Loop over dsigma_ii/dX_rr
655 for (int ii = 0; ii < 9;
656 ii++) { // ii represents components of stress tensor
657 dStress_dx(ii, 3 * dd + rr) += jac_stress(ii, 3 * rr + jj) * a;
658 }
659 }
660 }
661 }
663 }
664 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
665 EntityType col_type,
666 EntitiesFieldData::EntData &row_data,
667 EntitiesFieldData::EntData &col_data) {
669
670 if (commonData.skipThis) {
672 }
673
674 int nb_row = row_data.getIndices().size();
675 int nb_col = col_data.getIndices().size();
676 if (nb_row == 0)
678 if (nb_col == 0)
680 K.resize(nb_row, nb_col, false);
681 K.clear();
682 int nb_gauss_pts = row_data.getN().size1();
683 for (int gg = 0; gg != nb_gauss_pts; gg++) {
684 CHKERR get_dStress_dx(col_data, gg);
685 double val = getVolume() * getGaussPts()(3, gg);
686 // std::cerr << dStress_dx << std::endl;
687 dStress_dx *= val;
688 const MatrixAdaptor &diffN = row_data.getDiffN(gg, nb_row / 3);
689 { // integrate element stiffness matrix
690 for (int dd1 = 0; dd1 < nb_row / 3; dd1++) {
691 for (int rr1 = 0; rr1 < 3; rr1++) {
692 for (int dd2 = 0; dd2 < nb_col / 3; dd2++) {
693 for (int rr2 = 0; rr2 < 3; rr2++) {
694 K(3 * dd1 + rr1, 3 * dd2 + rr2) +=
695 (diffN(dd1, 0) * dStress_dx(3 * rr1 + 0, 3 * dd2 + rr2) +
696 diffN(dd1, 1) * dStress_dx(3 * rr1 + 1, 3 * dd2 + rr2) +
697 diffN(dd1, 2) * dStress_dx(3 * rr1 + 2, 3 * dd2 + rr2));
698 }
699 }
700 }
701 }
702 }
703 }
704 // std::cerr << "G " << getNumeredEntFiniteElementPtr()->getRefEnt() <<
705 // std::endl << K << std::endl;
706 CHKERR aSemble(row_side, col_side, row_type, col_type, row_data,
707 col_data);
708
710 }
711 };
712
713 /** \brief Assemble matrix
714 */
715 struct OpLhsdxdot : public AssembleMatrix {
717 OpLhsdxdot(CommonData &common_data)
718 : AssembleMatrix(common_data.spatialPositionName,
719 common_data.spatialPositionName),
720 commonData(common_data) {}
721 MatrixDouble dStress_dot;
722 MoFEMErrorCode get_dStress_dot(EntitiesFieldData::EntData &col_data,
723 int gg) {
725 int nb_col = col_data.getIndices().size();
726 dStress_dot.resize(9, nb_col, false);
727 dStress_dot.clear();
728 const MatrixAdaptor diffN = col_data.getDiffN(gg, nb_col / 3);
729 MatrixDouble &jac_stress = commonData.jacStress[gg];
730 for (int dd = 0; dd < nb_col / 3; dd++) { // DoFs in column
731 for (int jj = 0; jj < 3; jj++) { // cont. DoFs in column
732 double a = diffN(dd, jj);
733 for (int rr = 0; rr < 3; rr++) { // Loop over dsigma_ii/dX_rr
734 for (int ii = 0; ii < 9;
735 ii++) { // ii represents components of stress tensor
736 dStress_dot(ii, 3 * dd + rr) +=
737 jac_stress(ii, 9 + 3 * rr + jj) * a * getFEMethod()->ts_a;
738 }
739 }
740 }
741 }
743 }
744 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
745 EntityType col_type,
746 EntitiesFieldData::EntData &row_data,
747 EntitiesFieldData::EntData &col_data) {
749
750 if (commonData.skipThis) {
752 }
753
754 int nb_row = row_data.getIndices().size();
755 int nb_col = col_data.getIndices().size();
756 if (nb_row == 0)
758 if (nb_col == 0)
760 K.resize(nb_row, nb_col, false);
761 K.clear();
762 int nb_gauss_pts = row_data.getN().size1();
763 for (int gg = 0; gg != nb_gauss_pts; gg++) {
764 CHKERR get_dStress_dot(col_data, gg);
765 double val = getVolume() * getGaussPts()(3, gg);
766 // std::cerr << dStress_dot << std::endl;
767 dStress_dot *= val;
768 const MatrixAdaptor &diffN = row_data.getDiffN(gg, nb_row / 3);
769 { // integrate element stiffness matrix
770 for (int dd1 = 0; dd1 < nb_row / 3; dd1++) {
771 for (int rr1 = 0; rr1 < 3; rr1++) {
772 for (int dd2 = 0; dd2 < nb_col / 3; dd2++) {
773 for (int rr2 = 0; rr2 < 3; rr2++) {
774 K(3 * dd1 + rr1, 3 * dd2 + rr2) +=
775 (diffN(dd1, 0) * dStress_dot(3 * rr1 + 0, 3 * dd2 + rr2) +
776 diffN(dd1, 1) * dStress_dot(3 * rr1 + 1, 3 * dd2 + rr2) +
777 diffN(dd1, 2) * dStress_dot(3 * rr1 + 2, 3 * dd2 + rr2));
778 }
779 }
780 }
781 }
782 }
783 }
784 // std::cerr << "G " << getNumeredEntFiniteElementPtr()->getRefEnt() <<
785 // std::endl << K << std::endl;
786 CHKERR aSemble(row_side, col_side, row_type, col_type, row_data,
787 col_data);
788
790 }
791 };
792
793 MoFEMErrorCode setBlockDataMap() {
795
797 if (it->getName().compare(0, 6, "DAMPER") == 0) {
798 std::vector<double> data;
799 CHKERR it->getAttributes(data);
800 if (data.size() < 2) {
801 SETERRQ(PETSC_COMM_SELF, 1, "Data inconsistency");
802 }
803 CHKERR mField.get_moab().get_entities_by_type(
804 it->meshset, MBTET, blockMaterialDataMap[it->getMeshsetId()].tEts,
805 true);
806 blockMaterialDataMap[it->getMeshsetId()].gBeta = data[0];
807 blockMaterialDataMap[it->getMeshsetId()].vBeta = data[1];
808 }
809 }
811 }
812
813 MoFEMErrorCode setOperators(const int tag) {
815
816 for (auto &&fe_ptr : {&feRhs, &feLhs}) {
817 // CHKERR addHOOpsVol(commonData.meshNodePositionName, *fe_ptr, true, false,
818 // false, false);
819 fe_ptr->getOpPtrVector().push_back(
820 new OpCalculateVectorFieldGradient<3, 3>(
823 fe_ptr->getOpPtrVector().push_back(new OpGetDataAtGaussPts(
824 commonData.spatialPositionName, commonData, false, true, false));
825 fe_ptr->getOpPtrVector().push_back(
826 new OpCalculateVectorFieldGradientDot<3, 3>(
829 fe_ptr->getOpPtrVector().push_back(new OpGetDataAtGaussPts(
830 commonData.spatialPositionName, commonData, false, true, true));
831 }
832
833 // attach tags for each recorder
834 std::vector<int> tags;
835 tags.push_back(tag);
836
837 ConstitutiveEquationMap::iterator mit = constitutiveEquationMap.begin();
838 for (; mit != constitutiveEquationMap.end(); mit++) {
840 constitutiveEquationMap.at(mit->first);
841 // Right hand side operators
842 feRhs.getOpPtrVector().push_back(new OpJacobian(
843 commonData.spatialPositionName, tags, ce, commonData, true, false));
844 feRhs.getOpPtrVector().push_back(new OpRhsStress(commonData));
845
846 // Left hand side operators
847 feLhs.getOpPtrVector().push_back(new OpJacobian(
848 commonData.spatialPositionName, tags, ce, commonData, false, true));
849 feLhs.getOpPtrVector().push_back(new OpLhsdxdx(commonData));
850 feLhs.getOpPtrVector().push_back(new OpLhsdxdot(commonData));
851 }
852
854 }
855};
856
857#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()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ BLOCKSET
@ MOFEM_IMPOSSIBLE_CASE
Definition definitions.h:35
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
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
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
boost::ptr_map< int, KelvinVoigtDamper::ConstitutiveEquation< adouble > > ConstitutiveEquationMap
KelvinVoigtDamper(MoFEM::Interface &m_field)
MoFEMErrorCode setOperators(const int tag)
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)