v0.10.0
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 
491  static MoFEMErrorCode
492  setBlocks(MoFEM::Interface &m_field,
493  boost::shared_ptr<map<int, BlockData>> &block_sets_ptr);
494 
496  string element_name, string velocity_field_name,
497  string spatial_position_field_name,
498  string material_position_field_name = "MESH_NODE_POSITIONS",
499  bool ale = false, BitRefLevel bit = BitRefLevel());
500 
502  string element_name, string velocity_field_name,
503  string spatial_position_field_name,
504  string material_position_field_name = "MESH_NODE_POSITIONS",
505  bool ale = false, BitRefLevel bit = BitRefLevel());
506 
508  string element_name, string velocity_field_name,
509  string spatial_position_field_name,
510  string material_position_field_name = "MESH_NODE_POSITIONS",
511  bool ale = false, BitRefLevel bit = BitRefLevel(),
512  Range *intersected = NULL);
513 
515  string velocity_field_name, string spatial_position_field_name,
516  string material_position_field_name = "MESH_NODE_POSITIONS",
517  bool ale = false, bool linear = false);
518 
520  string velocity_field_name, string spatial_position_field_name,
521  string material_position_field_name = "MESH_NODE_POSITIONS",
522  bool ale = false);
523 
525  string velocity_field_name, string spatial_position_field_name,
526  string material_position_field_name = "MESH_NODE_POSITIONS",
527  Range *forces_on_entities_ptr = NULL);
528 
530  string velocity_field_name, string spatial_position_field_name,
531  string material_position_field_name = "MESH_NODE_POSITIONS",
532  bool linear = false);
533 
534  struct MatShellCtx {
535 
536  Mat K, M;
537  VecScatter scatterU, scatterV;
538  double ts_a; //,scale;
539 
541  MatShellCtx();
542  virtual ~MatShellCtx();
543 
544  Mat barK;
545  Vec u, v, Ku, Mv;
547 
549 
550  friend MoFEMErrorCode MultOpA(Mat A, Vec x, Vec f);
551  friend MoFEMErrorCode ZeroEntriesOp(Mat A);
552  };
553 
554  /** \brief Mult operator for shell matrix
555  *
556  * \f[
557  \left[
558  \begin{array}{cc}
559  \mathbf{M} & \mathbf{K} \\
560  \mathbf{I} & -\mathbf{I}a
561  \end{array}
562  \right]
563  \left[
564  \begin{array}{c}
565  \mathbf{v} \\
566  \mathbf{u}
567  \end{array}
568  \right] =
569  \left[
570  \begin{array}{c}
571  \mathbf{r}_u \\
572  \mathbf{r}_v
573  \end{array}
574  \right]
575  * \f]
576  *
577  */
578  static MoFEMErrorCode MultOpA(Mat A, Vec x, Vec f) {
580 
581  void *void_ctx;
582  CHKERR MatShellGetContext(A, &void_ctx);
583  MatShellCtx *ctx = (MatShellCtx *)void_ctx;
584  if (!ctx->iNitialized) {
585  CHKERR ctx->iNit();
586  }
587  CHKERR VecZeroEntries(f);
588  // Mult Ku
589  CHKERR VecScatterBegin(ctx->scatterU, x, ctx->u, INSERT_VALUES,
590  SCATTER_FORWARD);
591  CHKERR VecScatterEnd(ctx->scatterU, x, ctx->u, INSERT_VALUES,
592  SCATTER_FORWARD);
593  CHKERR MatMult(ctx->K, ctx->u, ctx->Ku);
594  CHKERR VecScatterBegin(ctx->scatterU, ctx->Ku, f, INSERT_VALUES,
595  SCATTER_REVERSE);
596  CHKERR VecScatterEnd(ctx->scatterU, ctx->Ku, f, INSERT_VALUES,
597  SCATTER_REVERSE);
598  // Mult Mv
599  CHKERR VecScatterBegin(ctx->scatterV, x, ctx->v, INSERT_VALUES,
600  SCATTER_FORWARD);
601  CHKERR VecScatterEnd(ctx->scatterV, x, ctx->v, INSERT_VALUES,
602  SCATTER_FORWARD);
603  CHKERR MatMult(ctx->M, ctx->v, ctx->Mv);
604  CHKERR VecScatterBegin(ctx->scatterU, ctx->Mv, f, ADD_VALUES,
605  SCATTER_REVERSE);
606  CHKERR VecScatterEnd(ctx->scatterU, ctx->Mv, f, ADD_VALUES,
607  SCATTER_REVERSE);
608  // Velocities
609  CHKERR VecAXPY(ctx->v, -ctx->ts_a, ctx->u);
610  // CHKERR VecScale(ctx->v,ctx->scale);
611  CHKERR VecScatterBegin(ctx->scatterV, ctx->v, f, INSERT_VALUES,
612  SCATTER_REVERSE);
613  CHKERR VecScatterEnd(ctx->scatterV, ctx->v, f, INSERT_VALUES,
614  SCATTER_REVERSE);
615  // Assemble
616  CHKERR VecAssemblyBegin(f);
617  CHKERR VecAssemblyEnd(f);
619  }
620 
623 
624  void *void_ctx;
625  CHKERR MatShellGetContext(A, &void_ctx);
626  MatShellCtx *ctx = (MatShellCtx *)void_ctx;
627  CHKERR MatZeroEntries(ctx->K);
628  CHKERR MatZeroEntries(ctx->M);
630  }
631 
632  struct PCShellCtx {
633 
634  Mat shellMat;
635  bool initPC; ///< check if PC is initialized
636 
637  PCShellCtx(Mat shell_mat) : shellMat(shell_mat), initPC(false) {}
638 
639  PC pC;
640 
642 
644 
645  friend MoFEMErrorCode PCShellSetUpOp(PC pc);
646  friend MoFEMErrorCode PCShellDestroy(PC pc);
647  friend MoFEMErrorCode PCShellApplyOp(PC pc, Vec f, Vec x);
648  };
649 
652 
653  void *void_ctx;
654  CHKERR PCShellGetContext(pc, &void_ctx);
655  PCShellCtx *ctx = (PCShellCtx *)void_ctx;
656  CHKERR ctx->iNit();
657  MatShellCtx *shell_mat_ctx;
658  CHKERR MatShellGetContext(ctx->shellMat, &shell_mat_ctx);
659  CHKERR PCSetFromOptions(ctx->pC);
660  CHKERR PCSetOperators(ctx->pC, shell_mat_ctx->barK, shell_mat_ctx->barK);
661  CHKERR PCSetUp(ctx->pC);
663  }
664 
667 
668  void *void_ctx;
669  CHKERR PCShellGetContext(pc, &void_ctx);
670  PCShellCtx *ctx = (PCShellCtx *)void_ctx;
671  CHKERR ctx->dEstroy();
673  }
674 
675  /** \brief apply pre-conditioner for shell matrix
676  *
677  * \f[
678  \left[
679  \begin{array}{cc}
680  \mathbf{M} & \mathbf{K} \\
681  \mathbf{I} & -\mathbf{I}a
682  \end{array}
683  \right]
684  \left[
685  \begin{array}{c}
686  \mathbf{v} \\
687  \mathbf{u}
688  \end{array}
689  \right] =
690  \left[
691  \begin{array}{c}
692  \mathbf{r}_u \\
693  \mathbf{r}_v
694  \end{array}
695  \right]
696  * \f]
697  *
698  * where \f$\mathbf{v} = \mathbf{r}_v + a\mathbf{u}\f$ and
699  \f$\mathbf{u}=(a\mathbf{M}+\mathbf{K})^{-1}(\mathbf{r}_u -
700  \mathbf{M}\mathbf{r}_v\f$.
701  *
702  */
703  static MoFEMErrorCode PCShellApplyOp(PC pc, Vec f, Vec x) {
705 
706  void *void_ctx;
707  CHKERR PCShellGetContext(pc, &void_ctx);
708  PCShellCtx *ctx = (PCShellCtx *)void_ctx;
709  MatShellCtx *shell_mat_ctx;
710  CHKERR MatShellGetContext(ctx->shellMat, &shell_mat_ctx);
711  // forward
712  CHKERR VecScatterBegin(shell_mat_ctx->scatterU, f, shell_mat_ctx->Ku,
713  INSERT_VALUES, SCATTER_FORWARD);
714  CHKERR VecScatterEnd(shell_mat_ctx->scatterU, f, shell_mat_ctx->Ku,
715  INSERT_VALUES, SCATTER_FORWARD);
716  CHKERR VecScatterBegin(shell_mat_ctx->scatterV, f, shell_mat_ctx->v,
717  INSERT_VALUES, SCATTER_FORWARD);
718  CHKERR VecScatterEnd(shell_mat_ctx->scatterV, f, shell_mat_ctx->v,
719  INSERT_VALUES, SCATTER_FORWARD);
720  // CHKERR VecScale(shell_mat_ctx->v,1/shell_mat_ctx->scale);
721  // apply pre-conditioner and calculate u
722  CHKERR MatMult(shell_mat_ctx->M, shell_mat_ctx->v,
723  shell_mat_ctx->Mv); // Mrv
724  CHKERR VecAXPY(shell_mat_ctx->Ku, -1, shell_mat_ctx->Mv); // f-Mrv
725  CHKERR PCApply(ctx->pC, shell_mat_ctx->Ku,
726  shell_mat_ctx->u); // u = (aM+K)^(-1)(ru-Mrv)
727  // VecView(shell_mat_ctx->u,PETSC_VIEWER_STDOUT_WORLD);
728  // calculate velocities
729  CHKERR VecAXPY(shell_mat_ctx->v, shell_mat_ctx->ts_a,
730  shell_mat_ctx->u); // v = v + a*u
731  // VecView(shell_mat_ctx->v,PETSC_VIEWER_STDOUT_WORLD);
732  // reverse
733  CHKERR VecZeroEntries(x);
734  CHKERR VecScatterBegin(shell_mat_ctx->scatterU, shell_mat_ctx->u, x,
735  INSERT_VALUES, SCATTER_REVERSE);
736  CHKERR VecScatterEnd(shell_mat_ctx->scatterU, shell_mat_ctx->u, x,
737  INSERT_VALUES, SCATTER_REVERSE);
738  CHKERR VecScatterBegin(shell_mat_ctx->scatterV, shell_mat_ctx->v, x,
739  INSERT_VALUES, SCATTER_REVERSE);
740  CHKERR VecScatterEnd(shell_mat_ctx->scatterV, shell_mat_ctx->v, x,
741  INSERT_VALUES, SCATTER_REVERSE);
742  CHKERR VecAssemblyBegin(x);
743  CHKERR VecAssemblyEnd(x);
745  }
746 
747  struct ShellResidualElement : public FEMethod {
750 
751  // variables bellow need to be set by user
752  MatShellCtx *shellMatCtx; ///< pointer to shell matrix
753 
755 
757  };
758 
759 #ifdef __DIRICHLET_HPP__
760 
761  /** \brief blocked element/problem
762  *
763  * Blocked element run loops for different problem than TS problem. It is
764  * used to calculate matrices of shell matrix.
765  *
766  */
767  struct ShellMatrixElement : public FEMethod {
768 
770  ShellMatrixElement(MoFEM::Interface &m_field);
771 
772  typedef std::pair<std::string, FEMethod *> PairNameFEMethodPtr;
773  typedef std::vector<PairNameFEMethodPtr> LoopsToDoType;
774  LoopsToDoType loopK; ///< methods to calculate K shell matrix
775  LoopsToDoType loopM; ///< methods to calculate M shell matrix
776  LoopsToDoType loopAuxM; ///< methods to calculate derivatives of inertia
777  ///< forces over displacements shell matrix
778 
779  // variables bellow need to be set by user
780  string problemName; ///< name of shell problem
781  MatShellCtx *shellMatCtx; ///< pointer to shell matrix
782  DirichletDisplacementBc *DirichletBcPtr; ///< boundary conditions
783 
784  MoFEMErrorCode preProcess();
785  };
786 
787 #endif //__DIRICHLET_HPP__
788 };
789 
790 #endif //__CONVECTIVE_MASS_ELEMENT_HPP
791 
792 /**
793 * \defgroup convective_mass_elem Mass Element
794 * \ingroup user_modules
795 **/
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:509
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:76
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:516
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)
ForcesAndSourcesCore::UserDataOperator UserDataOperator
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)
DataForcesAndSourcesCore::EntData EntData
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
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:69
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:102
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:604
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
constexpr int order
std::vector< std::vector< double * > > jacVelRowPtr
ConvectiveMassElement(MoFEM::Interface &m_field, short int tag)
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
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:74
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