v0.9.1
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 /* Implementation of convective mass element
8  *
9  * This file is part of MoFEM.
10  * MoFEM is free software: you can redistribute it and/or modify it under
11  * the terms of the GNU Lesser General Public License as published by the
12  * Free Software Foundation, either version 3 of the License, or (at your
13  * option) any later version.
14  *
15  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
16  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18  * License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public
21  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
22 
23 #ifndef __KELVIN_VOIGT_DAMPER_HPP__
24 #define __KELVIN_VOIGT_DAMPER_HPP__
25 
26 #ifndef WITH_ADOL_C
27 #error "MoFEM need to be compiled with ADOL-C"
28 #endif
29 
30 /** \brief Implementation of Kelvin Voigt Damper
31 \ingroup nonlinear_elastic_elem
32 
33 */
35 
37 
39 
40  /** \brief Dumper material parameters
41  \ingroup nonlinear_elastic_elem
42  */
44  Range tEts;
45  double vBeta; ///< Poisson ration spring alpha
46  double gBeta; ///< Sheer modulus spring alpha
47  bool lInear;
48  BlockMaterialData() : vBeta(0), gBeta(1), lInear(false) {}
49  };
50 
51  std::map<int, BlockMaterialData> blockMaterialDataMap;
52 
53  /** \brief Constitutive model functions
54  \ingroup nonlinear_elastic_elem
55 
56  */
57  template <typename TYPE> struct ConstitutiveEquation {
58 
60 
62  virtual ~ConstitutiveEquation() {}
63 
64  ublas::matrix<TYPE> F; ///< Gradient of deformation
65  ublas::matrix<TYPE> FDot; ///< Rate of gradient of deformation
66  ublas::matrix<TYPE> gradientUDot; ///< Rate of gradient of displacements
67  ublas::matrix<TYPE> engineringStrainDot;
68  ublas::matrix<TYPE>
69  dashpotCauchyStress; ///< Stress generated by spring beta
70  ublas::matrix<TYPE>
71  dashpotFirstPiolaKirchhoffStress; ///< Stress generated by spring beta
72 
74  TYPE J; ///< Jacobian of gradient of deformation
75  ublas::matrix<TYPE> invF; ///< Inverse of gradient of deformation
76 
77  /** \brief Calculate determinant of 3x3 matrix
78  */
79  MoFEMErrorCode dEterminant(ublas::matrix<TYPE> a, TYPE &det) {
81  // a11a22a33
82  //+a21a32a13
83  //+a31a12a23
84  //-a11a32a23
85  //-a31a22a13
86  //-a21a12a33
87  // http://www.cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html
88  // http://mathworld.wolfram.com/MatrixInverse.html
89  det = a(0, 0) * a(1, 1) * a(2, 2) + a(1, 0) * a(2, 1) * a(0, 2) +
90  a(2, 0) * a(0, 1) * a(1, 2) - a(0, 0) * a(2, 1) * a(1, 2) -
91  a(2, 0) * a(1, 1) * a(0, 2) - a(1, 0) * a(0, 1) * a(2, 2);
93  }
94 
95  /** \brief Calculate inverse of 3x3 matrix
96  */
97  MoFEMErrorCode iNvert(TYPE det, ublas::matrix<TYPE> a,
98  ublas::matrix<TYPE> &inv_a) {
100  //
101  inv_a.resize(3, 3);
102  // http://www.cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html
103  // http://mathworld.wolfram.com/MatrixInverse.html
104  inv_a(0, 0) = a(1, 1) * a(2, 2) - a(1, 2) * a(2, 1);
105  inv_a(0, 1) = a(0, 2) * a(2, 1) - a(0, 1) * a(2, 2);
106  inv_a(0, 2) = a(0, 1) * a(1, 2) - a(0, 2) * a(1, 1);
107  inv_a(1, 0) = a(1, 2) * a(2, 0) - a(1, 0) * a(2, 2);
108  inv_a(1, 1) = a(0, 0) * a(2, 2) - a(0, 2) * a(2, 0);
109  inv_a(1, 2) = a(0, 2) * a(1, 0) - a(0, 0) * a(1, 2);
110  inv_a(2, 0) = a(1, 0) * a(2, 1) - a(1, 1) * a(2, 0);
111  inv_a(2, 1) = a(0, 1) * a(2, 0) - a(0, 0) * a(2, 1);
112  inv_a(2, 2) = a(0, 0) * a(1, 1) - a(0, 1) * a(1, 0);
113  inv_a /= det;
115  }
116 
117  /** \brief Calculate strain rate
118 
119  \f[
120  \dot{\varepsilon}_{ij} = \frac{1}{2}
121  \left(
122  \frac{\partial v_i}{\partial X_j}
123  +
124  \frac{\partial v_j}{\partial X_i}
125  \right)
126  \f]
127 
128  */
131  gradientUDot.resize(3, 3, false);
132  noalias(gradientUDot) = FDot;
133  for (int ii = 0; ii < 3; ii++) {
134  gradientUDot(ii, ii) -= 1;
135  }
137  for (int ii = 0; ii < 3; ii++) {
139  }
140  engineringStrainDot.resize(3, 3, false);
141  noalias(engineringStrainDot) = gradientUDot + trans(gradientUDot);
142  engineringStrainDot *= 0.5;
144  }
145 
146  /** \brief Calculate Cauchy dashpot stress
147 
148  Calculate dashpot Cauchy stress. It has to be pull back to reference
149  configuration before use in total Lagrangian formulation.
150 
151  \f[
152  \sigma^\beta_{ij} = 2G^\beta\left[
153  \dot{\varepsilon}_{ij}
154  + \frac{v^\beta}{1-2v^\beta}\dot{\varepsilon}_{kk}\delta_{ij}
155  \right]
156  \f]
157 
158  */
161  dashpotCauchyStress.resize(3, 3, false);
162  double a = 2.0 * dAta.gBeta;
163  double b = a * (dAta.vBeta / (1.0 - 2.0 * dAta.vBeta));
165  for (int ii = 0; ii < 3; ii++) {
167  }
169  }
170 
171  /** \brief Calculate First Piola-Kirchhoff Stress Dashpot stress
172 
173  \f[
174  P^\beta_{ij} = J \sigma^\beta_{ik} F^{-1}_{jk}
175  \f]
176 
177  */
180  dashpotFirstPiolaKirchhoffStress.resize(3, 3, false);
181  if (dAta.lInear) {
183  } else {
184 
185  invF.resize(3, 3, false);
186  CHKERR dEterminant(F, J);
187  CHKERR iNvert(J, F, invF);
189  J * prod(dashpotCauchyStress, trans(invF));
190  }
192  }
193  };
194 
195  typedef boost::ptr_map<int, KelvinVoigtDamper::ConstitutiveEquation<adouble>>
198 
199  /** \brief Common data for nonlinear_elastic_elem model
200  \ingroup nonlinear_elastic_elem
201  */
202  struct CommonData {
203 
206 
207  std::map<std::string, std::vector<VectorDouble>> dataAtGaussPts;
208  std::map<std::string, std::vector<MatrixDouble>> gradAtGaussPts;
209 
210  std::vector<MatrixDouble> dashpotFirstPiolaKirchhoffStress;
211 
212  std::vector<double *> jacRowPtr;
213  std::vector<MatrixDouble> jacStress;
214 
215  bool recordOn;
216  bool skipThis;
217  std::map<int, int> nbActiveVariables, nbActiveResults;
218 
219  CommonData() : recordOn(true), skipThis(true) {}
220  };
222 
223  /// \brief definition of volume element
225 
227  int addToRule; ///< Takes into account HO geometry
228 
229  DamperFE(MoFEM::Interface &m_field, CommonData &common_data)
231  commonData(common_data), addToRule(1) {}
232 
233  int getRule(int order) { return order + addToRule; }
234 
236 
239 
240  if (ts_ctx == CTX_TSSETIFUNCTION) {
241 
242  CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
244  commonData.spatialPositionNameDot, COL, ts_u_t, INSERT_VALUES,
245  SCATTER_REVERSE);
246  }
247 
249  }
250 
252 
254 
256 
257  if (ts_ctx == CTX_TSSETIFUNCTION) {
258  CHKERR VecAssemblyBegin(ts_F);
259  CHKERR VecAssemblyEnd(ts_F);
260  }
261  if (ts_ctx == CTX_TSSETIJACOBIAN) {
262  CHKERR MatAssemblyBegin(ts_B, MAT_FLUSH_ASSEMBLY);
263  CHKERR MatAssemblyEnd(ts_B, MAT_FLUSH_ASSEMBLY);
264  }
265 
267  }
268  };
269 
271 
273  : mField(m_field), feRhs(m_field, commonData),
274  feLhs(m_field, commonData) {}
275 
278 
280  bool calcVal;
281  bool calcGrad;
282  EntityType zeroAtType;
283 
284  OpGetDataAtGaussPts(const std::string field_name, CommonData &common_data,
285  bool calc_val, bool calc_grad,
286  EntityType zero_at_type = MBVERTEX)
288  field_name, UserDataOperator::OPCOL),
289  commonData(common_data), calcVal(calc_val), calcGrad(calc_grad),
290  zeroAtType(zero_at_type) {}
291 
292  /** \brief Operator field value
293  *
294  */
295  MoFEMErrorCode doWork(int side, EntityType type,
298 
299  int nb_dofs = data.getFieldData().size();
300  if (nb_dofs == 0) {
302  }
303  int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
304  int nb_gauss_pts = data.getN().size1();
305 
306  // Initialize
307  if (calcVal) {
308  commonData.dataAtGaussPts[rowFieldName].resize(nb_gauss_pts);
309  for (int gg = 0; gg < nb_gauss_pts; gg++) {
310  commonData.dataAtGaussPts[rowFieldName][gg].resize(rank, false);
311  }
312  }
313  if (calcGrad) {
314  commonData.gradAtGaussPts[rowFieldName].resize(nb_gauss_pts);
315  for (int gg = 0; gg < nb_gauss_pts; gg++) {
316  commonData.gradAtGaussPts[rowFieldName][gg].resize(rank, 3, false);
317  }
318  }
319 
320  // Zero values
321  if (type == zeroAtType) {
322  for (int gg = 0; gg < nb_gauss_pts; gg++) {
323  if (calcVal) {
325  }
326  if (calcGrad) {
328  }
329  }
330  }
331 
332  VectorDouble &values = data.getFieldData();
333 
334  // Calculate values at integration points
335  for (int gg = 0; gg < nb_gauss_pts; gg++) {
336  VectorDouble N = data.getN(gg, nb_dofs / rank);
337  MatrixDouble diffN = data.getDiffN(gg, nb_dofs / rank);
338  for (int dd = 0; dd < nb_dofs / rank; dd++) {
339  for (int rr1 = 0; rr1 < rank; rr1++) {
340  if (calcVal) {
342  N[dd] * values[rank * dd + rr1];
343  }
344  if (calcGrad) {
345  for (int rr2 = 0; rr2 < 3; rr2++) {
346  commonData.gradAtGaussPts[rowFieldName][gg](rr1, rr2) +=
347  diffN(dd, rr2) * values[rank * dd + rr1];
348  }
349  }
350  }
351  }
352  }
353 
355  }
356  };
357 
358  struct OpJacobian
360 
361  std::vector<int> tagS;
364 
367  bool &recordOn;
368  std::map<int, int> &nbActiveVariables;
369  std::map<int, int> &nbActiveResults;
370 
371  OpJacobian(const std::string field_name, std::vector<int> tags,
373  CommonData &common_data, bool calculate_residual,
374  bool calculate_jacobian)
376  field_name, UserDataOperator::OPROW),
377  tagS(tags), cE(ce), commonData(common_data),
378  calculateResidualBool(calculate_residual),
379  calculateJacobianBool(calculate_jacobian),
380  recordOn(common_data.recordOn),
382  nbActiveResults(common_data.nbActiveResults) {}
383 
386 
389 
390  if (tagS[DAMPERSTRESS] < 0) {
392  }
393 
394  cE.F.resize(3, 3, false);
395  cE.FDot.resize(3, 3, false);
396  MatrixDouble &F =
398  MatrixDouble &F_dot =
400  trace_on(tagS[DAMPERSTRESS]);
401  {
402  // Activate gradient of defamation
404  for (int dd1 = 0; dd1 < 3; dd1++) {
405  for (int dd2 = 0; dd2 < 3; dd2++) {
406  cE.F(dd1, dd2) <<= F(dd1, dd2);
408  }
409  }
410  for (int dd1 = 0; dd1 < 3; dd1++) {
411  for (int dd2 = 0; dd2 < 3; dd2++) {
412  cE.FDot(dd1, dd2) <<= F_dot(dd1, dd2);
414  }
415  }
416 
417  // Do calculations
421 
422  // Results
425  commonData.dashpotFirstPiolaKirchhoffStress[0].resize(3, 3, false);
426  for (int d1 = 0; d1 < 3; d1++) {
427  for (int d2 = 0; d2 < 3; d2++) {
431  }
432  }
433  }
434  trace_off();
436  }
437 
440 
441  int r;
442  // play recorder for values
443  r = ::function(tagS[te], nbActiveResults[tagS[te]],
444  nbActiveVariables[tagS[te]], &activeVariables[0], ptr);
445  if (r < 3) { // function is locally analytic
446  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
447  "ADOL-C function evaluation with error r = %d", r);
448  }
449 
451  }
452 
455 
456  try {
457  int r;
458  r = jacobian(tagS[te], nbActiveResults[tagS[te]],
460  &(commonData.jacRowPtr[0]));
461  if (r < 3) {
462  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
463  "ADOL-C function evaluation with error");
464  }
465  } catch (const std::exception &ex) {
466  std::ostringstream ss;
467  ss << "throw in method: " << ex.what() << std::endl;
468  SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
469  }
471  }
472 
475 
476  if (tagS[DAMPERSTRESS] < 0) {
478  }
479 
481  for (int gg = 0; gg < nbGaussPts; gg++) {
482 
483  MatrixDouble &F =
485  MatrixDouble &F_dot =
487  int nb_active_variables = 0;
488 
489  // Activate gradient of defamation
490  for (int dd1 = 0; dd1 < 3; dd1++) {
491  for (int dd2 = 0; dd2 < 3; dd2++) {
492  activeVariables[nb_active_variables++] = F(dd1, dd2);
493  }
494  }
495  // Activate rate of gradient of defamation
496  for (int dd1 = 0; dd1 < 3; dd1++) {
497  for (int dd2 = 0; dd2 < 3; dd2++) {
498  activeVariables[nb_active_variables++] = F_dot(dd1, dd2);
499  }
500  }
501 
502  if (nb_active_variables != nbActiveVariables[tagS[DAMPERSTRESS]]) {
503  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE,
504  "Number of active variables does not much");
505  }
506 
507  if (calculateResidualBool) {
508  if (gg == 0) {
510  }
511  commonData.dashpotFirstPiolaKirchhoffStress[gg].resize(3, 3, false);
513  DAMPERSTRESS,
515  }
516 
517  if (calculateJacobianBool) {
518  if (gg == 0) {
521  }
524  false);
525  for (int dd = 0; dd < nbActiveResults[tagS[DAMPERSTRESS]]; dd++) {
527  }
529  }
530  }
531 
533  }
534 
535  MoFEMErrorCode doWork(int row_side, EntityType row_type,
538 
539  if (row_type != MBVERTEX)
541  nbGaussPts = row_data.getN().size1();
542 
543  commonData.skipThis = false;
544  if (cE.dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
545  cE.dAta.tEts.end()) {
546  commonData.skipThis = true;
548  }
549 
550  if (recordOn) {
552  }
554 
556  }
557  };
558 
561  AssembleVector(string field_name)
563  field_name, UserDataOperator::OPROW) {}
564 
566  MoFEMErrorCode aSemble(int row_side, EntityType row_type,
569  int nb_dofs = row_data.getIndices().size();
570  int *indices_ptr = &row_data.getIndices()[0];
571  CHKERR VecSetValues(getFEMethod()->ts_F, nb_dofs, indices_ptr, &nF[0],
572  ADD_VALUES);
574  }
575  };
576 
577  /** \brief Assemble internal force vector
578  \ingroup nonlinear_elastic_elem
579 
580  */
581  struct OpRhsStress : public AssembleVector {
583  OpRhsStress(CommonData &common_data)
584  : AssembleVector(common_data.spatialPositionName),
585  commonData(common_data) {}
586  MoFEMErrorCode doWork(int row_side, EntityType row_type,
589 
590  if (commonData.skipThis) {
592  }
593 
594  int nb_dofs = row_data.getIndices().size();
595  if (!nb_dofs) {
597  }
598  nF.resize(nb_dofs, false);
599  nF.clear();
600  int nb_gauss_pts = row_data.getN().size1();
601  for (int gg = 0; gg != nb_gauss_pts; gg++) {
602  const MatrixAdaptor &diffN = row_data.getDiffN(gg, nb_dofs / 3);
603  const MatrixDouble &stress =
605  double val = getVolume() * getGaussPts()(3, gg);
606  if (getHoGaussPtsDetJac().size() > 0) {
607  val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
608  }
609  for (int dd = 0; dd < nb_dofs / 3; dd++) {
610  for (int rr = 0; rr < 3; rr++) {
611  for (int nn = 0; nn < 3; nn++) {
612  nF[3 * dd + rr] += val * diffN(dd, nn) * stress(rr, nn);
613  }
614  }
615  }
616  }
617  CHKERR aSemble(row_side, row_type, row_data);
619  }
620  };
621 
624  AssembleMatrix(string row_name, string col_name)
626  row_name, col_name, UserDataOperator::OPROWCOL) {}
627 
629  MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type,
630  EntityType col_type,
634  int nb_row = row_data.getIndices().size();
635  int nb_col = col_data.getIndices().size();
636  int *row_indices_ptr = &row_data.getIndices()[0];
637  int *col_indices_ptr = &col_data.getIndices()[0];
638  CHKERR MatSetValues(getFEMethod()->ts_B, nb_row, row_indices_ptr, nb_col,
639  col_indices_ptr, &K(0, 0), ADD_VALUES);
640  if (sYmm) {
641  // Assemble of diagonal terms
642  if (row_side != col_side || row_type != col_type) {
643  transK.resize(nb_col, nb_row, false);
644  noalias(transK) = trans(K);
645  CHKERR MatSetValues(getFEMethod()->ts_B, nb_col, col_indices_ptr,
646  nb_row, row_indices_ptr, &transK(0, 0),
647  ADD_VALUES);
648  }
649  }
651  }
652  };
653 
654  /** \brief Assemble matrix
655  */
656  struct OpLhsdxdx : public AssembleMatrix {
658  OpLhsdxdx(CommonData &common_data)
659  : AssembleMatrix(common_data.spatialPositionName,
660  common_data.spatialPositionName),
661  commonData(common_data) {}
664  int gg) {
666  int nb_col = col_data.getIndices().size();
667  dStress_dx.resize(9, nb_col, false);
668  dStress_dx.clear();
669  const MatrixAdaptor diffN = col_data.getDiffN(gg, nb_col / 3);
670  MatrixDouble &jac_stress = commonData.jacStress[gg];
671  for (int dd = 0; dd < nb_col / 3; dd++) { // DoFs in column
672  for (int jj = 0; jj < 3; jj++) { // cont. DoFs in column
673  double a = diffN(dd, jj);
674  for (int rr = 0; rr < 3; rr++) { // Loop over dsigma_ii/dX_rr
675  for (int ii = 0; ii < 9;
676  ii++) { // ii represents components of stress tensor
677  dStress_dx(ii, 3 * dd + rr) += jac_stress(ii, 3 * rr + jj) * a;
678  }
679  }
680  }
681  }
683  }
684  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
685  EntityType col_type,
689 
690  if (commonData.skipThis) {
692  }
693 
694  int nb_row = row_data.getIndices().size();
695  int nb_col = col_data.getIndices().size();
696  if (nb_row == 0)
698  if (nb_col == 0)
700  K.resize(nb_row, nb_col, false);
701  K.clear();
702  int nb_gauss_pts = row_data.getN().size1();
703  for (int gg = 0; gg != nb_gauss_pts; gg++) {
704  CHKERR get_dStress_dx(col_data, gg);
705  double val = getVolume() * getGaussPts()(3, gg);
706  if (getHoGaussPtsDetJac().size() > 0) {
707  val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
708  }
709  // std::cerr << dStress_dx << std::endl;
710  dStress_dx *= val;
711  const MatrixAdaptor &diffN = row_data.getDiffN(gg, nb_row / 3);
712  { // integrate element stiffness matrix
713  for (int dd1 = 0; dd1 < nb_row / 3; dd1++) {
714  for (int rr1 = 0; rr1 < 3; rr1++) {
715  for (int dd2 = 0; dd2 < nb_col / 3; dd2++) {
716  for (int rr2 = 0; rr2 < 3; rr2++) {
717  K(3 * dd1 + rr1, 3 * dd2 + rr2) +=
718  (diffN(dd1, 0) * dStress_dx(3 * rr1 + 0, 3 * dd2 + rr2) +
719  diffN(dd1, 1) * dStress_dx(3 * rr1 + 1, 3 * dd2 + rr2) +
720  diffN(dd1, 2) * dStress_dx(3 * rr1 + 2, 3 * dd2 + rr2));
721  }
722  }
723  }
724  }
725  }
726  }
727  // std::cerr << "G " << getNumeredEntFiniteElementPtr()->getRefEnt() <<
728  // std::endl << K << std::endl;
729  CHKERR aSemble(row_side, col_side, row_type, col_type, row_data,
730  col_data);
731 
733  }
734  };
735 
736  /** \brief Assemble matrix
737  */
738  struct OpLhsdxdot : public AssembleMatrix {
740  OpLhsdxdot(CommonData &common_data)
741  : AssembleMatrix(common_data.spatialPositionName,
742  common_data.spatialPositionName),
743  commonData(common_data) {}
746  int gg) {
748  int nb_col = col_data.getIndices().size();
749  dStress_dot.resize(9, nb_col, false);
750  dStress_dot.clear();
751  const MatrixAdaptor diffN = col_data.getDiffN(gg, nb_col / 3);
752  MatrixDouble &jac_stress = commonData.jacStress[gg];
753  for (int dd = 0; dd < nb_col / 3; dd++) { // DoFs in column
754  for (int jj = 0; jj < 3; jj++) { // cont. DoFs in column
755  double a = diffN(dd, jj);
756  for (int rr = 0; rr < 3; rr++) { // Loop over dsigma_ii/dX_rr
757  for (int ii = 0; ii < 9;
758  ii++) { // ii represents components of stress tensor
759  dStress_dot(ii, 3 * dd + rr) +=
760  jac_stress(ii, 9 + 3 * rr + jj) * a * getFEMethod()->ts_a;
761  }
762  }
763  }
764  }
766  }
767  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
768  EntityType col_type,
772 
773  if (commonData.skipThis) {
775  }
776 
777  int nb_row = row_data.getIndices().size();
778  int nb_col = col_data.getIndices().size();
779  if (nb_row == 0)
781  if (nb_col == 0)
783  K.resize(nb_row, nb_col, false);
784  K.clear();
785  int nb_gauss_pts = row_data.getN().size1();
786  for (int gg = 0; gg != nb_gauss_pts; gg++) {
787  CHKERR get_dStress_dot(col_data, gg);
788  double val = getVolume() * getGaussPts()(3, gg);
789  if (getHoGaussPtsDetJac().size() > 0) {
790  val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
791  }
792  // std::cerr << dStress_dot << std::endl;
793  dStress_dot *= val;
794  const MatrixAdaptor &diffN = row_data.getDiffN(gg, nb_row / 3);
795  { // integrate element stiffness matrix
796  for (int dd1 = 0; dd1 < nb_row / 3; dd1++) {
797  for (int rr1 = 0; rr1 < 3; rr1++) {
798  for (int dd2 = 0; dd2 < nb_col / 3; dd2++) {
799  for (int rr2 = 0; rr2 < 3; rr2++) {
800  K(3 * dd1 + rr1, 3 * dd2 + rr2) +=
801  (diffN(dd1, 0) * dStress_dot(3 * rr1 + 0, 3 * dd2 + rr2) +
802  diffN(dd1, 1) * dStress_dot(3 * rr1 + 1, 3 * dd2 + rr2) +
803  diffN(dd1, 2) * dStress_dot(3 * rr1 + 2, 3 * dd2 + rr2));
804  }
805  }
806  }
807  }
808  }
809  }
810  // std::cerr << "G " << getNumeredEntFiniteElementPtr()->getRefEnt() <<
811  // std::endl << K << std::endl;
812  CHKERR aSemble(row_side, col_side, row_type, col_type, row_data,
813  col_data);
814 
816  }
817  };
818 
821 
823  if (it->getName().compare(0, 6, "DAMPER") == 0) {
824  std::vector<double> data;
825  CHKERR it->getAttributes(data);
826  if (data.size() < 2) {
827  SETERRQ(PETSC_COMM_SELF, 1, "Data inconsistency");
828  }
829  CHKERR mField.get_moab().get_entities_by_type(
830  it->meshset, MBTET, blockMaterialDataMap[it->getMeshsetId()].tEts,
831  true);
832  blockMaterialDataMap[it->getMeshsetId()].gBeta = data[0];
833  blockMaterialDataMap[it->getMeshsetId()].vBeta = data[1];
834  }
835  }
837  }
838 
839  MoFEMErrorCode setOperators(const int tag) {
841 
842  DamperFE *fe_ptr[] = {&feRhs, &feLhs};
843  for (int ss = 0; ss < 2; ss++) {
844  fe_ptr[ss]->getOpPtrVector().push_back(new OpGetDataAtGaussPts(
845  commonData.spatialPositionName, commonData, false, true));
846  fe_ptr[ss]->getOpPtrVector().push_back(new OpGetDataAtGaussPts(
848  }
849 
850  // attach tags for each recorder
851  std::vector<int> tags;
852  tags.push_back(tag);
853 
854  ConstitutiveEquationMap::iterator mit = constitutiveEquationMap.begin();
855  for (; mit != constitutiveEquationMap.end(); mit++) {
857  constitutiveEquationMap.at(mit->first);
858  // Right hand side operators
859  feRhs.getOpPtrVector().push_back(new OpJacobian(
860  commonData.spatialPositionName, tags, ce, commonData, true, false));
861  feRhs.getOpPtrVector().push_back(new OpRhsStress(commonData));
862 
863  // Left hand side operators
864  feLhs.getOpPtrVector().push_back(new OpJacobian(
865  commonData.spatialPositionName, tags, ce, commonData, false, true));
866  feLhs.getOpPtrVector().push_back(new OpLhsdxdx(commonData));
867  feLhs.getOpPtrVector().push_back(new OpLhsdxdot(commonData));
868  }
869 
871  }
872 };
873 
874 #endif //__KELVIN_VOIGT_DAMPER_HPP__
ConstitutiveEquation(BlockMaterialData &data)
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
MoFEMErrorCode setBlockDataMap()
Deprecated interface functions.
std::map< std::string, std::vector< VectorDouble > > dataAtGaussPts
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
virtual moab::Interface & get_moab()=0
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
OpGetDataAtGaussPts(const std::string field_name, CommonData &common_data, bool calc_val, bool calc_grad, EntityType zero_at_type=MBVERTEX)
MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
std::vector< MatrixDouble > jacStress
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
int addToRule
Takes into account HO geometry.
Mat & ts_B
Preconditioner for ts_A.
std::vector< MatrixDouble > dashpotFirstPiolaKirchhoffStress
ublas::matrix< TYPE > gradientUDot
Rate of gradient of displacements.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
std::map< int, int > & nbActiveVariables
boost::ptr_map< int, KelvinVoigtDamper::ConstitutiveEquation< adouble > > ConstitutiveEquationMap
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513
AssembleMatrix(string row_name, string col_name)
virtual MoFEMErrorCode calculateEngineeringStrainDot()
Calculate strain rate.
UserDataOperator(const FieldSpace space, const char type=OPLAST, const bool symm=true)
TYPE J
Jacobian of gradient of deformation.
Assemble internal force vector.
KelvinVoigtDamper::ConstitutiveEquation< adouble > & cE
MoFEMErrorCode calculateAtIntPtsDamperStress()
MoFEMErrorCode setOperators(const int tag)
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
std::map< int, int > & nbActiveResults
MoFEMErrorCode preProcess()
function is run at the beginning of loop
ublas::matrix< TYPE > dashpotFirstPiolaKirchhoffStress
Stress generated by spring beta.
MoFEM::Interface & mField
const Problem * problemPtr
raw pointer to problem
ublas::matrix< TYPE > dashpotCauchyStress
Stress generated by spring beta.
OpRhsStress(CommonData &common_data)
OpLhsdxdot(CommonData &common_data)
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
MoFEMErrorCode calculateJacobian(TagEvaluate te)
Vec & ts_u_t
time derivative of state vector
bool sYmm
If true assume that matrix is symmetric structure.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
virtual MoFEMErrorCode calculateFirstPiolaKirchhoffStress()
Calculate First Piola-Kirchhoff Stress Dashpot stress.
PetscReal ts_a
shift for U_tt (see PETSc Time Solver)
ublas::matrix< TYPE > F
Gradient of deformation.
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
MoFEMErrorCode get_dStress_dx(DataForcesAndSourcesCore::EntData &col_data, int gg)
ublas::matrix< TYPE > FDot
Rate of gradient of deformation.
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
Definition: Types.hpp:121
constexpr int order
DataForcesAndSourcesCore::EntData EntData
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
OpJacobian(const std::string field_name, std::vector< int > tags, KelvinVoigtDamper::ConstitutiveEquation< adouble > &ce, CommonData &common_data, bool calculate_residual, bool calculate_jacobian)
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
MoFEMErrorCode aSemble(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
std::map< int, BlockMaterialData > blockMaterialDataMap
double vBeta
Poisson ration spring alpha.
std::vector< double * > jacRowPtr
OpLhsdxdx(CommonData &common_data)
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
#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.
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
#define CHKERR
Inline error check.
Definition: definitions.h:601
virtual MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator field value.
Vec & ts_F
residual vector
Implementation of Kelvin Voigt Damper.
DamperFE(MoFEM::Interface &m_field, CommonData &common_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
virtual MoFEMErrorCode calculateDashpotCauchyStress()
Calculate Cauchy dashpot stress.
ConstitutiveEquationMap constitutiveEquationMap
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
std::map< int, int > nbActiveVariables
MoFEMErrorCode calculateFunction(TagEvaluate te, double *ptr)
TSContext ts_ctx
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
std::map< int, int > nbActiveResults
virtual MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode dEterminant(ublas::matrix< TYPE > a, TYPE &det)
Calculate determinant of 3x3 matrix.
const int N
Definition: speed_test.cpp:3
MoFEMErrorCode get_dStress_dot(DataForcesAndSourcesCore::EntData &col_data, int gg)
ublas::matrix< TYPE > invF
Inverse of gradient of deformation.
Common data for nonlinear_elastic_elem model.
double gBeta
Sheer modulus spring alpha.
KelvinVoigtDamper(MoFEM::Interface &m_field)
definition of volume element
MoFEMErrorCode iNvert(TYPE det, ublas::matrix< TYPE > a, ublas::matrix< TYPE > &inv_a)
Calculate inverse of 3x3 matrix.