v0.9.1
ConvectiveMassElement.hpp
Go to the documentation of this file.
1 /** \file ConvectiveMassElement.hpp
2  * \brief Operators and data structures for mass and convective mass element
3  * \ingroup convective_mass_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 __CONVECTIVE_MASS_ELEMENT_HPP
24 #define __CONVECTIVE_MASS_ELEMENT_HPP
25 
26 #ifndef WITH_ADOL_C
27 #error "MoFEM need to be compiled with ADOL-C"
28 #endif
29 
30 /** \brief structure grouping operators and data used for calculation of mass
31  * (convective) element \ingroup convective_mass_elem \ingroup
32  * nonlinear_elastic_elem
33  *
34  * In order to assemble matrices and right hand vectors, the loops over
35  * elements, entities over that elements and finally loop over integration
36  * points are executed.
37  *
38  * Following implementation separate those three celeries of loops and to each
39  * loop attach operator.
40  *
41  */
43 
44  /// \brief definition of volume element
46 
47  Mat A;
48  Vec F;
49  bool initV; ///< check if ghost vector used to accumalte Kinetin energy is
50  ///< created
51 
52  MyVolumeFE(MoFEM::Interface &m_field);
53 
54  /** \brief it is used to calculate nb. of Gauss integration points
55  *
56  * for more details pleas look
57  * Reference:
58  *
59  * Albert Nijenhuis, Herbert Wilf,
60  * Combinatorial Algorithms for Computers and Calculators,
61  * Second Edition,
62  * Academic Press, 1978,
63  * ISBN: 0-12-519260-6,
64  * LC: QA164.N54.
65  *
66  * More details about algorithm
67  * http://people.sc.fsu.edu/~jburkardt/cpp_src/gm_rule/gm_rule.html
68  **/
69  int getRule(int order);
70 
71  Vec V;
72  double eNergy;
73 
76  };
77 
78  MyVolumeFE feMassRhs; ///< calculate right hand side for tetrahedral elements
80  return feMassRhs;
81  } ///< get rhs volume element
82  MyVolumeFE feMassLhs; ///< calculate left hand side for tetrahedral
83  ///< elements,i.e. mass element
85  return feMassLhs;
86  } ///< get lhs volume element
87  MyVolumeFE feMassAuxLhs; ///< calculate left hand side for tetrahedral
88  ///< elements for Kuu shell matrix
90  return feMassAuxLhs;
91  } ///< get lhs volume element for Kuu shell matrix
92 
93  MyVolumeFE feVelRhs; ///< calculate right hand side for tetrahedral elements
94  MyVolumeFE &getLoopFeVelRhs() { return feVelRhs; } ///< get rhs volume element
95  MyVolumeFE feVelLhs; ///< calculate left hand side for tetrahedral elements
96  MyVolumeFE &getLoopFeVelLhs() { return feVelLhs; } ///< get lhs volume element
97 
98  MyVolumeFE feTRhs; ///< calculate right hand side for tetrahedral elements
99  MyVolumeFE &getLoopFeTRhs() { return feTRhs; } ///< get rhs volume element
100  MyVolumeFE feTLhs; ///< calculate left hand side for tetrahedral elements
101  MyVolumeFE &getLoopFeTLhs() { return feTLhs; } ///< get lhs volume element
102 
103  MyVolumeFE feEnergy; ///< calculate kinetic energy
105  return feEnergy;
106  } ///< get kinetic energy element
107 
109  short int tAg;
110 
111  ConvectiveMassElement(MoFEM::Interface &m_field, short int tag);
112 
113  /** \brief data for calculation inertia forces
114  * \ingroup user_modules
115  */
116  struct BlockData {
117  double rho0; ///< reference density
118  VectorDouble a0; ///< constant acceleration
119  Range tEts; ///< elements in block set
120  };
121  std::map<int, BlockData>
122  setOfBlocks; ///< maps block set id with appropriate BlockData
123 
124  /** \brief common data used by volume elements
125  * \ingroup user_modules
126  */
127  struct CommonData {
128 
129  bool lInear;
131  CommonData() : lInear(false), staticOnly(false) {}
132 
133  std::map<std::string, std::vector<VectorDouble>> dataAtGaussPts;
134  std::map<std::string, std::vector<MatrixDouble>> gradAtGaussPts;
138  std::vector<VectorDouble> valVel;
139  std::vector<std::vector<double *>> jacVelRowPtr;
140  std::vector<MatrixDouble> jacVel;
141  std::vector<VectorDouble> valMass;
142  std::vector<std::vector<double *>> jacMassRowPtr;
143  std::vector<MatrixDouble> jacMass;
144  std::vector<VectorDouble> valT;
145  std::vector<std::vector<double *>> jacTRowPtr;
146  std::vector<MatrixDouble> jacT;
147  };
149 
150  boost::ptr_vector<MethodForForceScaling> methodsOp;
151 
154 
155  std::vector<VectorDouble> &valuesAtGaussPts;
156  std::vector<MatrixDouble> &gradientAtGaussPts;
157  const EntityType zeroAtType;
158 
159  OpGetDataAtGaussPts(const std::string field_name,
160  std::vector<VectorDouble> &values_at_gauss_pts,
161  std::vector<MatrixDouble> &gardient_at_gauss_pts);
162 
163  /** \brief operator calculating deformation gradient
164  *
165  */
166  MoFEMErrorCode doWork(int side, EntityType type,
168  };
169 
171  OpGetCommonDataAtGaussPts(const std::string field_name,
172  CommonData &common_data);
173  };
174 
176 
177  template <typename TYPE>
179  ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>> &a,
180  TYPE &det) {
182  // a11a22a33
183  //+a21a32a13
184  //+a31a12a23
185  //-a11a32a23
186  //-a31a22a13
187  //-a21a12a33
188  // http://www.cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html
189  // http://mathworld.wolfram.com/MatrixInverse.html
190  det = a(0, 0) * a(1, 1) * a(2, 2) + a(1, 0) * a(2, 1) * a(0, 2) +
191  a(2, 0) * a(0, 1) * a(1, 2) - a(0, 0) * a(2, 1) * a(1, 2) -
192  a(2, 0) * a(1, 1) * a(0, 2) - a(1, 0) * a(0, 1) * a(2, 2);
194  }
195 
196  template <typename TYPE>
198  TYPE det,
199  ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>> &a,
200  ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>>
201  &inv_a) {
203  //
204  inv_a.resize(3, 3);
205  // http://www.cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html
206  // http://mathworld.wolfram.com/MatrixInverse.html
207  inv_a(0, 0) = a(1, 1) * a(2, 2) - a(1, 2) * a(2, 1);
208  inv_a(0, 1) = a(0, 2) * a(2, 1) - a(0, 1) * a(2, 2);
209  inv_a(0, 2) = a(0, 1) * a(1, 2) - a(0, 2) * a(1, 1);
210  inv_a(1, 0) = a(1, 2) * a(2, 0) - a(1, 0) * a(2, 2);
211  inv_a(1, 1) = a(0, 0) * a(2, 2) - a(0, 2) * a(2, 0);
212  inv_a(1, 2) = a(0, 2) * a(1, 0) - a(0, 0) * a(1, 2);
213  inv_a(2, 0) = a(1, 0) * a(2, 1) - a(1, 1) * a(2, 0);
214  inv_a(2, 1) = a(0, 1) * a(2, 0) - a(0, 0) * a(2, 1);
215  inv_a(2, 2) = a(0, 0) * a(1, 1) - a(0, 1) * a(1, 0);
216  inv_a /= det;
218  }
219  };
220 
224 
227  int tAg;
228  bool jAcobian;
229  bool &lInear;
230  bool fieldDisp;
231 
232  boost::ptr_vector<MethodForForceScaling> &methodsOp;
233 
234  OpMassJacobian(const std::string field_name, BlockData &data,
235  CommonData &common_data,
236  boost::ptr_vector<MethodForForceScaling> &methods_op,
237  int tag, bool linear = false);
238 
239  ublas::vector<adouble, ublas::bounded_array<adouble, 3>> a, dot_W, dp_dt,
240  a_res;
241  ublas::matrix<adouble, ublas::row_major, ublas::bounded_array<adouble, 9>>
242  h, H, invH, F, g, G;
243  std::vector<double> active;
244 
245  MoFEMErrorCode doWork(int row_side, EntityType row_type,
247  };
248 
251 
254 
255  OpMassRhs(const std::string field_name, BlockData &data,
256  CommonData &common_data);
257 
259 
260  MoFEMErrorCode doWork(int row_side, EntityType row_type,
262  };
263 
267 
271 
272  OpMassLhs_dM_dv(const std::string vel_field, const std::string field_name,
273  BlockData &data, CommonData &common_data,
274  Range *forcesonlyonentities_ptr = NULL);
275 
277 
279  int gg);
280 
281  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
282  EntityType col_type,
285  };
286 
288 
289  OpMassLhs_dM_dx(const std::string field_name, const std::string col_field,
290  BlockData &data, CommonData &common_data);
291 
293  };
294 
296 
297  OpMassLhs_dM_dX(const std::string field_name, const std::string col_field,
298  BlockData &data, CommonData &common_data);
299 
301  };
302 
305 
308  Vec *Vptr;
309  bool &lInear;
310 
311  OpEnergy(const std::string field_name, BlockData &data,
312  CommonData &common_data, Vec *v_ptr);
313 
316 
317  MoFEMErrorCode doWork(int row_side, EntityType row_type,
319  };
320 
324 
327  int tAg;
329 
330  OpVelocityJacobian(const std::string field_name, BlockData &data,
331  CommonData &common_data, int tag, bool jacobian = true);
332 
333  ublas::vector<adouble, ublas::bounded_array<adouble, 3>> a_res;
334  ublas::vector<adouble, ublas::bounded_array<adouble, 3>> v, dot_w, dot_W;
335  ublas::matrix<adouble, ublas::row_major, ublas::bounded_array<adouble, 9>>
336  h, H, invH, F;
337  ublas::vector<adouble, ublas::bounded_array<adouble, 3>> dot_u;
339 
340  std::vector<double> active;
341 
342  MoFEMErrorCode doWork(int row_side, EntityType row_type,
344  };
345 
349 
352 
353  OpVelocityRhs(const std::string field_name, BlockData &data,
354  CommonData &common_data);
355 
357 
358  MoFEMErrorCode doWork(int row_side, EntityType row_type,
360  };
361 
363 
364  OpVelocityLhs_dV_dv(const std::string vel_field,
365  const std::string field_name, BlockData &data,
366  CommonData &common_data);
367 
369  int gg);
370  };
371 
373 
374  OpVelocityLhs_dV_dx(const std::string vel_field,
375  const std::string field_name, BlockData &data,
376  CommonData &common_data);
377 
379  int gg);
380  };
381 
383 
384  OpVelocityLhs_dV_dX(const std::string vel_field,
385  const std::string field_name, BlockData &data,
386  CommonData &common_data);
387 
389  int gg);
390  };
391 
395 
398  int tAg;
399  bool jAcobian;
400  bool fieldDisp;
401 
402  OpEshelbyDynamicMaterialMomentumJacobian(const std::string field_name,
403  BlockData &data,
404  CommonData &common_data, int tag,
405  bool jacobian = true);
406 
407  ublas::vector<adouble, ublas::bounded_array<adouble, 3>> a, v, a_T;
408  ublas::matrix<adouble, ublas::row_major, ublas::bounded_array<adouble, 9>>
409  g, H, invH, h, F, G;
411 
412  MoFEMErrorCode doWork(int row_side, EntityType row_type,
414  };
415 
419 
423 
424  OpEshelbyDynamicMaterialMomentumRhs(const std::string field_name,
425  BlockData &data,
426  CommonData &common_data,
427  Range *forcesonlyonentities_ptr);
428 
430 
431  MoFEMErrorCode doWork(int row_side, EntityType row_type,
433  };
434 
436 
437  OpEshelbyDynamicMaterialMomentumLhs_dv(const std::string vel_field,
438  const std::string field_name,
439  BlockData &data,
440  CommonData &common_data,
441  Range *forcesonlyonentities_ptr);
442 
444  int gg);
445  };
446 
449 
450  OpEshelbyDynamicMaterialMomentumLhs_dx(const std::string vel_field,
451  const std::string field_name,
452  BlockData &data,
453  CommonData &common_data,
454  Range *forcesonlyonentities_ptr);
455 
457  int gg);
458  };
459 
462 
463  OpEshelbyDynamicMaterialMomentumLhs_dX(const std::string vel_field,
464  const std::string field_name,
465  BlockData &data,
466  CommonData &common_data,
467  Range *forcesonlyonentities_ptr);
468 
470  int gg);
471  };
472 
473  struct UpdateAndControl : public FEMethod {
474 
476  TS tS;
477  const std::string velocityField;
478  const std::string spatialPositionField;
479 
481  UpdateAndControl(MoFEM::Interface &m_field, TS _ts,
482  const std::string velocity_field,
483  const std::string spatial_position_field);
484 
487  };
488 
490 
492  string element_name, string velocity_field_name,
493  string spatial_position_field_name,
494  string material_position_field_name = "MESH_NODE_POSITIONS",
495  bool ale = false, BitRefLevel bit = BitRefLevel());
496 
498  string element_name, string velocity_field_name,
499  string spatial_position_field_name,
500  string material_position_field_name = "MESH_NODE_POSITIONS",
501  bool ale = false, BitRefLevel bit = BitRefLevel());
502 
504  string element_name, string velocity_field_name,
505  string spatial_position_field_name,
506  string material_position_field_name = "MESH_NODE_POSITIONS",
507  bool ale = false, BitRefLevel bit = BitRefLevel(),
508  Range *intersected = NULL);
509 
511  string velocity_field_name, string spatial_position_field_name,
512  string material_position_field_name = "MESH_NODE_POSITIONS",
513  bool ale = false, bool linear = false);
514 
516  string velocity_field_name, string spatial_position_field_name,
517  string material_position_field_name = "MESH_NODE_POSITIONS",
518  bool ale = false);
519 
521  string velocity_field_name, string spatial_position_field_name,
522  string material_position_field_name = "MESH_NODE_POSITIONS",
523  Range *forces_on_entities_ptr = NULL);
524 
526  string velocity_field_name, string spatial_position_field_name,
527  string material_position_field_name = "MESH_NODE_POSITIONS",
528  bool linear = false);
529 
530  struct MatShellCtx {
531 
532  Mat K, M;
533  VecScatter scatterU, scatterV;
534  double ts_a; //,scale;
535 
537  MatShellCtx();
538  virtual ~MatShellCtx();
539 
540  Mat barK;
541  Vec u, v, Ku, Mv;
543 
545 
546  friend MoFEMErrorCode MultOpA(Mat A, Vec x, Vec f);
547  friend MoFEMErrorCode ZeroEntriesOp(Mat A);
548  };
549 
550  /** \brief Mult operator for shell matrix
551  *
552  * \f[
553  \left[
554  \begin{array}{cc}
555  \mathbf{M} & \mathbf{K} \\
556  \mathbf{I} & -\mathbf{I}a
557  \end{array}
558  \right]
559  \left[
560  \begin{array}{c}
561  \mathbf{v} \\
562  \mathbf{u}
563  \end{array}
564  \right] =
565  \left[
566  \begin{array}{c}
567  \mathbf{r}_u \\
568  \mathbf{r}_v
569  \end{array}
570  \right]
571  * \f]
572  *
573  */
574  static MoFEMErrorCode MultOpA(Mat A, Vec x, Vec f) {
576 
577  void *void_ctx;
578  CHKERR MatShellGetContext(A, &void_ctx);
579  MatShellCtx *ctx = (MatShellCtx *)void_ctx;
580  if (!ctx->iNitialized) {
581  CHKERR ctx->iNit();
582  }
583  CHKERR VecZeroEntries(f);
584  // Mult Ku
585  CHKERR VecScatterBegin(ctx->scatterU, x, ctx->u, INSERT_VALUES,
586  SCATTER_FORWARD);
587  CHKERR VecScatterEnd(ctx->scatterU, x, ctx->u, INSERT_VALUES,
588  SCATTER_FORWARD);
589  CHKERR MatMult(ctx->K, ctx->u, ctx->Ku);
590  CHKERR VecScatterBegin(ctx->scatterU, ctx->Ku, f, INSERT_VALUES,
591  SCATTER_REVERSE);
592  CHKERR VecScatterEnd(ctx->scatterU, ctx->Ku, f, INSERT_VALUES,
593  SCATTER_REVERSE);
594  // Mult Mv
595  CHKERR VecScatterBegin(ctx->scatterV, x, ctx->v, INSERT_VALUES,
596  SCATTER_FORWARD);
597  CHKERR VecScatterEnd(ctx->scatterV, x, ctx->v, INSERT_VALUES,
598  SCATTER_FORWARD);
599  CHKERR MatMult(ctx->M, ctx->v, ctx->Mv);
600  CHKERR VecScatterBegin(ctx->scatterU, ctx->Mv, f, ADD_VALUES,
601  SCATTER_REVERSE);
602  CHKERR VecScatterEnd(ctx->scatterU, ctx->Mv, f, ADD_VALUES,
603  SCATTER_REVERSE);
604  // Velocities
605  CHKERR VecAXPY(ctx->v, -ctx->ts_a, ctx->u);
606  // CHKERR VecScale(ctx->v,ctx->scale);
607  CHKERR VecScatterBegin(ctx->scatterV, ctx->v, f, INSERT_VALUES,
608  SCATTER_REVERSE);
609  CHKERR VecScatterEnd(ctx->scatterV, ctx->v, f, INSERT_VALUES,
610  SCATTER_REVERSE);
611  // Assemble
612  CHKERR VecAssemblyBegin(f);
613  CHKERR VecAssemblyEnd(f);
615  }
616 
619 
620  void *void_ctx;
621  CHKERR MatShellGetContext(A, &void_ctx);
622  MatShellCtx *ctx = (MatShellCtx *)void_ctx;
623  CHKERR MatZeroEntries(ctx->K);
624  CHKERR MatZeroEntries(ctx->M);
626  }
627 
628  struct PCShellCtx {
629 
630  Mat shellMat;
631  bool initPC; ///< check if PC is initialized
632 
633  PCShellCtx(Mat shell_mat) : shellMat(shell_mat), initPC(false) {}
634 
635  PC pC;
636 
638 
640 
641  friend MoFEMErrorCode PCShellSetUpOp(PC pc);
642  friend MoFEMErrorCode PCShellDestroy(PC pc);
643  friend MoFEMErrorCode PCShellApplyOp(PC pc, Vec f, Vec x);
644  };
645 
648 
649  void *void_ctx;
650  CHKERR PCShellGetContext(pc, &void_ctx);
651  PCShellCtx *ctx = (PCShellCtx *)void_ctx;
652  CHKERR ctx->iNit();
653  MatShellCtx *shell_mat_ctx;
654  CHKERR MatShellGetContext(ctx->shellMat, &shell_mat_ctx);
655  CHKERR PCSetFromOptions(ctx->pC);
656  CHKERR PCSetOperators(ctx->pC, shell_mat_ctx->barK, shell_mat_ctx->barK);
657  CHKERR PCSetUp(ctx->pC);
659  }
660 
663 
664  void *void_ctx;
665  CHKERR PCShellGetContext(pc, &void_ctx);
666  PCShellCtx *ctx = (PCShellCtx *)void_ctx;
667  CHKERR ctx->dEstroy();
669  }
670 
671  /** \brief apply pre-conditioner for shell matrix
672  *
673  * \f[
674  \left[
675  \begin{array}{cc}
676  \mathbf{M} & \mathbf{K} \\
677  \mathbf{I} & -\mathbf{I}a
678  \end{array}
679  \right]
680  \left[
681  \begin{array}{c}
682  \mathbf{v} \\
683  \mathbf{u}
684  \end{array}
685  \right] =
686  \left[
687  \begin{array}{c}
688  \mathbf{r}_u \\
689  \mathbf{r}_v
690  \end{array}
691  \right]
692  * \f]
693  *
694  * where \f$\mathbf{v} = \mathbf{r}_v + a\mathbf{u}\f$ and
695  \f$\mathbf{u}=(a\mathbf{M}+\mathbf{K})^{-1}(\mathbf{r}_u -
696  \mathbf{M}\mathbf{r}_v\f$.
697  *
698  */
699  static MoFEMErrorCode PCShellApplyOp(PC pc, Vec f, Vec x) {
701 
702  void *void_ctx;
703  CHKERR PCShellGetContext(pc, &void_ctx);
704  PCShellCtx *ctx = (PCShellCtx *)void_ctx;
705  MatShellCtx *shell_mat_ctx;
706  CHKERR MatShellGetContext(ctx->shellMat, &shell_mat_ctx);
707  // forward
708  CHKERR VecScatterBegin(shell_mat_ctx->scatterU, f, shell_mat_ctx->Ku,
709  INSERT_VALUES, SCATTER_FORWARD);
710  CHKERR VecScatterEnd(shell_mat_ctx->scatterU, f, shell_mat_ctx->Ku,
711  INSERT_VALUES, SCATTER_FORWARD);
712  CHKERR VecScatterBegin(shell_mat_ctx->scatterV, f, shell_mat_ctx->v,
713  INSERT_VALUES, SCATTER_FORWARD);
714  CHKERR VecScatterEnd(shell_mat_ctx->scatterV, f, shell_mat_ctx->v,
715  INSERT_VALUES, SCATTER_FORWARD);
716  // CHKERR VecScale(shell_mat_ctx->v,1/shell_mat_ctx->scale);
717  // apply pre-conditioner and calculate u
718  CHKERR MatMult(shell_mat_ctx->M, shell_mat_ctx->v,
719  shell_mat_ctx->Mv); // Mrv
720  CHKERR VecAXPY(shell_mat_ctx->Ku, -1, shell_mat_ctx->Mv); // f-Mrv
721  CHKERR PCApply(ctx->pC, shell_mat_ctx->Ku,
722  shell_mat_ctx->u); // u = (aM+K)^(-1)(ru-Mrv)
723  // VecView(shell_mat_ctx->u,PETSC_VIEWER_STDOUT_WORLD);
724  // calculate velocities
725  CHKERR VecAXPY(shell_mat_ctx->v, shell_mat_ctx->ts_a,
726  shell_mat_ctx->u); // v = v + a*u
727  // VecView(shell_mat_ctx->v,PETSC_VIEWER_STDOUT_WORLD);
728  // reverse
729  CHKERR VecZeroEntries(x);
730  CHKERR VecScatterBegin(shell_mat_ctx->scatterU, shell_mat_ctx->u, x,
731  INSERT_VALUES, SCATTER_REVERSE);
732  CHKERR VecScatterEnd(shell_mat_ctx->scatterU, shell_mat_ctx->u, x,
733  INSERT_VALUES, SCATTER_REVERSE);
734  CHKERR VecScatterBegin(shell_mat_ctx->scatterV, shell_mat_ctx->v, x,
735  INSERT_VALUES, SCATTER_REVERSE);
736  CHKERR VecScatterEnd(shell_mat_ctx->scatterV, shell_mat_ctx->v, x,
737  INSERT_VALUES, SCATTER_REVERSE);
738  CHKERR VecAssemblyBegin(x);
739  CHKERR VecAssemblyEnd(x);
741  }
742 
743  struct ShellResidualElement : public FEMethod {
746 
747  // variables bellow need to be set by user
748  MatShellCtx *shellMatCtx; ///< pointer to shell matrix
749 
751 
753  };
754 
755 #ifdef __DIRICHLET_HPP__
756 
757  /** \brief blocked element/problem
758  *
759  * Blocked element run loops for different problem than TS problem. It is
760  * used to calculate matrices of shell matrix.
761  *
762  */
763  struct ShellMatrixElement : public FEMethod {
764 
766  ShellMatrixElement(MoFEM::Interface &m_field);
767 
768  typedef std::pair<std::string, FEMethod *> PairNameFEMethodPtr;
769  typedef std::vector<PairNameFEMethodPtr> LoopsToDoType;
770  LoopsToDoType loopK; ///< methods to calculate K shell matrix
771  LoopsToDoType loopM; ///< methods to calculate M shell matrix
772  LoopsToDoType loopAuxM; ///< methods to calculate derivatives of inertia
773  ///< forces over displacements shell matrix
774 
775  // variables bellow need to be set by user
776  string problemName; ///< name of shell problem
777  MatShellCtx *shellMatCtx; ///< pointer to shell matrix
778  DirichletDisplacementBc *DirichletBcPtr; ///< boundary conditions
779 
780  MoFEMErrorCode preProcess();
781  };
782 
783 #endif //__DIRICHLET_HPP__
784 };
785 
786 #endif //__CONVECTIVE_MASS_ELEMENT_HPP
787 
788 /**
789 * \defgroup convective_mass_elem Mass Element
790 * \ingroup user_modules
791 **/
static MoFEMErrorCode PCShellDestroy(PC pc)
ublas::vector< adouble, ublas::bounded_array< adouble, 3 > > dot_w
ublas::matrix< adouble, ublas::row_major, ublas::bounded_array< adouble, 9 > > H
ublas::matrix< adouble, ublas::row_major, ublas::bounded_array< adouble, 9 > > F
Deprecated interface functions.
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
MoFEMErrorCode getJac(DataForcesAndSourcesCore::EntData &col_data, int gg)
MyVolumeFE & getLoopFeTLhs()
get lhs volume element
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
ublas::vector< adouble, ublas::bounded_array< adouble, 3 > > dot_W
virtual MoFEMErrorCode getJac(DataForcesAndSourcesCore::EntData &col_data, int gg)
friend MoFEMErrorCode PCShellDestroy(PC pc)
MoFEMErrorCode setShellMatrixMassOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool linear=false)
MyVolumeFE feVelLhs
calculate left hand side for tetrahedral elements
MyVolumeFE feEnergy
calculate kinetic energy
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:506
MyVolumeFE & getLoopFeEnergy()
get kinetic energy element
OpGetDataAtGaussPts(const std::string field_name, std::vector< VectorDouble > &values_at_gauss_pts, std::vector< MatrixDouble > &gardient_at_gauss_pts)
static MoFEMErrorCode PCShellSetUpOp(PC pc)
OpEnergy(const std::string field_name, BlockData &data, CommonData &common_data, Vec *v_ptr)
ublas::vector< adouble, ublas::bounded_array< adouble, 3 > > a_res
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
structure grouping operators and data used for calculation of mass (convective) element \ nonlinear_e...
OpEshelbyDynamicMaterialMomentumJacobian(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian=true)
MoFEMErrorCode setKinematicEshelbyOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", Range *forces_on_entities_ptr=NULL)
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
ublas::vector< adouble, ublas::bounded_array< adouble, 3 > > dp_dt
friend MoFEMErrorCode PCShellApplyOp(PC pc, Vec f, Vec x)
MyVolumeFE & getLoopFeTRhs()
get rhs volume element
MoFEMErrorCode setVelocityOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false)
UpdateAndControl(MoFEM::Interface &m_field, TS _ts, const std::string velocity_field, const std::string spatial_position_field)
static MoFEMErrorCode MultOpA(Mat A, Vec x, Vec f)
Mult operator for shell matrix.
ublas::vector< adouble, ublas::bounded_array< adouble, 3 > > dot_u
virtual MoFEMErrorCode getJac(DataForcesAndSourcesCore::EntData &col_data, int gg)
boost::ptr_vector< MethodForForceScaling > methodsOp
MyVolumeFE & getLoopFeVelLhs()
get lhs volume element
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513
ublas::matrix< adouble, ublas::row_major, ublas::bounded_array< adouble, 9 > > H
MoFEMErrorCode addConvectiveMassElement(string element_name, string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, BitRefLevel bit=BitRefLevel())
ublas::matrix< adouble, ublas::row_major, ublas::bounded_array< adouble, 9 > > invH
common data used by volume elements
std::map< std::string, std::vector< VectorDouble > > dataAtGaussPts
OpMassJacobian(const std::string field_name, BlockData &data, CommonData &common_data, boost::ptr_vector< MethodForForceScaling > &methods_op, int tag, bool linear=false)
ublas::vector< adouble, ublas::bounded_array< adouble, 3 > > dot_W
boost::ptr_vector< MethodForForceScaling > & methodsOp
ublas::matrix< adouble, ublas::row_major, ublas::bounded_array< adouble, 9 > > g
MatShellCtx * shellMatCtx
pointer to shell matrix
OpMassRhs(const std::string field_name, BlockData &data, CommonData &common_data)
OpVelocityRhs(const std::string field_name, BlockData &data, CommonData &common_data)
OpVelocityLhs_dV_dv(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
std::vector< std::vector< double * > > jacTRowPtr
ublas::matrix< adouble, ublas::row_major, ublas::bounded_array< adouble, 9 > > F
ublas::matrix< adouble, ublas::row_major, ublas::bounded_array< adouble, 9 > > G
ublas::matrix< adouble, ublas::row_major, ublas::bounded_array< adouble, 9 > > h
OpMassLhs_dM_dX(const std::string field_name, const std::string col_field, BlockData &data, CommonData &common_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
ublas::vector< adouble, ublas::bounded_array< adouble, 3 > > a_res
int getRule(int order)
it is used to calculate nb. of Gauss integration points
ublas::matrix< adouble, ublas::row_major, ublas::bounded_array< adouble, 9 > > h
MoFEMErrorCode dEterminant(ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 >> &a, TYPE &det)
constexpr int order
DataForcesAndSourcesCore::EntData EntData
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
VectorDouble a0
constant acceleration
ublas::vector< adouble, ublas::bounded_array< adouble, 3 > > a
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
operator calculating deformation gradient
friend MoFEMErrorCode MultOpA(Mat A, Vec x, Vec f)
ublas::vector< adouble, ublas::bounded_array< adouble, 3 > > a
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:43
ublas::vector< adouble, ublas::bounded_array< adouble, 3 > > v
MyVolumeFE & getLoopFeVelRhs()
get rhs volume element
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
MyVolumeFE feTLhs
calculate left hand side for tetrahedral elements
OpGetCommonDataAtGaussPts(const std::string field_name, CommonData &common_data)
MoFEMErrorCode setConvectiveMassOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, bool linear=false)
ublas::matrix< adouble, ublas::row_major, ublas::bounded_array< adouble, 9 > > F
OpMassLhs_dM_dv(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr=NULL)
MoFEMErrorCode getJac(DataForcesAndSourcesCore::EntData &col_data, int gg)
ublas::matrix< adouble, ublas::row_major, ublas::bounded_array< adouble, 9 > > invH
virtual MoFEMErrorCode getJac(DataForcesAndSourcesCore::EntData &col_data, int gg)
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:97
OpVelocityLhs_dV_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
OpEshelbyDynamicMaterialMomentumRhs(const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
#define CHKERR
Inline error check.
Definition: definitions.h:601
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
static MoFEMErrorCode PCShellApplyOp(PC pc, Vec f, Vec x)
apply pre-conditioner for shell matrix
MyVolumeFE & getLoopFeMassAuxLhs()
get lhs volume element for Kuu shell matrix
MyVolumeFE & getLoopFeMassLhs()
get lhs volume element
virtual MoFEMErrorCode getJac(DataForcesAndSourcesCore::EntData &col_data, int gg)
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
MyVolumeFE(MoFEM::Interface &m_field)
ublas::matrix< adouble, ublas::row_major, ublas::bounded_array< adouble, 9 > > g
std::vector< std::vector< double * > > jacVelRowPtr
ConvectiveMassElement(MoFEM::Interface &m_field, short int tag)
ForcesAndSourcesCore::UserDataOperator UserDataOperator
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:87
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
OpVelocityJacobian(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian=true)
MoFEMErrorCode iNvert(TYPE det, ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 >> &a, ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 >> &inv_a)
ublas::vector< adouble, ublas::bounded_array< adouble, 3 > > v
static MoFEMErrorCode ZeroEntriesOp(Mat A)
virtual MoFEMErrorCode getJac(DataForcesAndSourcesCore::EntData &col_data, int gg)
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
virtual MoFEMErrorCode getJac(DataForcesAndSourcesCore::EntData &col_data, int gg)
OpEshelbyDynamicMaterialMomentumLhs_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
MoFEMErrorCode addVelocityElement(string element_name, string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, BitRefLevel bit=BitRefLevel())
friend MoFEMErrorCode PCShellSetUpOp(PC pc)
friend MoFEMErrorCode ZeroEntriesOp(Mat A)
ublas::vector< adouble, ublas::bounded_array< adouble, 3 > > a_T
OpEshelbyDynamicMaterialMomentumLhs_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
ublas::matrix< adouble, ublas::row_major, ublas::bounded_array< adouble, 9 > > G
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
virtual MoFEMErrorCode getJac(DataForcesAndSourcesCore::EntData &col_data, int gg)
ublas::matrix< adouble, ublas::row_major, ublas::bounded_array< adouble, 9 > > H
ublas::matrix< adouble, ublas::row_major, ublas::bounded_array< adouble, 9 > > h
std::vector< std::vector< double * > > jacMassRowPtr
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
OpEshelbyDynamicMaterialMomentumLhs_dv(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
data for calculation inertia forces
MoFEMErrorCode addEshelbyDynamicMaterialMomentum(string element_name, string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, BitRefLevel bit=BitRefLevel(), Range *intersected=NULL)
OpVelocityLhs_dV_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
OpMassLhs_dM_dx(const std::string field_name, const std::string col_field, BlockData &data, CommonData &common_data)
MyVolumeFE feVelRhs
calculate right hand side for tetrahedral elements
MyVolumeFE feMassRhs
calculate right hand side for tetrahedral elements
MyVolumeFE & getLoopFeMassRhs()
get rhs volume element
ublas::matrix< adouble, ublas::row_major, ublas::bounded_array< adouble, 9 > > invH
MyVolumeFE feTRhs
calculate right hand side for tetrahedral elements
bool initPC
check if PC is initialized