v0.9.2
ForcesAndSourcesCore.hpp
Go to the documentation of this file.
1 /** \file ForcesAndSourcesCore.hpp
2  \brief Implementation of elements on entities.
3 
4  Those element are inherited by user to implement specific implementation of
5  particular problem.
6 
7 */
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 __FORCES_AND_SOURCES_CORE__HPP__
24 #define __FORCES_AND_SOURCES_CORE__HPP__
25 
26 using namespace boost::numeric;
27 
28 namespace MoFEM {
29 
30 /** \brief structure to get information form mofem into DataForcesAndSourcesCore
31  * \ingroup mofem_forces_and_sources
32  *
33  */
34 struct ForcesAndSourcesCore : public FEMethod {
35 
37 
39 
40  typedef boost::function<int(int order_row, int order_col, int order_data)>
42 
43  /**
44  * \brief Hook to get rule
45  *
46  * \todo check preferred format how works with gcc and clang,
47  * see
48  * <http://www.boost.org/doc/libs/1_64_0/doc/html/function/tutorial.html#idp247873024>
49  */
51 
52  /**
53  * @brief Set function to calculate integration rule
54  *
55  */
57 
58  /** \brief Data operator to do calculations at integration points.
59  * \ingroup mofem_forces_and_sources
60 
61  Is inherited and implemented by user to do calculations. It can be used in
62  many different ways but typically is used to integrate matrices (f.e.
63  stiffness matrix) and the right hand vector (f.e. force vector).
64 
65  Note: It is assumed that operator is executed for symmetric problem. That
66  means that is executed for not repeating entities on finite element. For
67  example on triangle we have nodes, 3 edges and one face. Because of symmetry
68  calculations are for: nodes-nodes, nodes-edge0, nodes-edge_1, nodes-edge_2,
69  nodes-face,
70  edges_1-edges_1, edges_1-edge_1, edge_1-edge_2,
71  edge_1-edge_1, edge_1-edge_2,
72  edge_2-edge_2,
73  edge_1-face, edge_1-face, edges_3-face,
74  face - face
75 
76  In case of non symmetric problem in addition calculations of lower off
77  diagonal terms. F.e. edge_1-edge_0, esges_3-edge_0, edge_3-edge_1,
78 
79  In that case class variable UserDataOperator::sYmm = false;
80 
81  NoteL: By default sYmm is set for symmetric problems
82 
83  */
84  struct UserDataOperator : public DataOperator {
85 
86  /**
87  * \brief Controls loop over entities on element
88  *
89  * OPROW is used if row vector is assembled
90  * OPCOL is usually used if column vector is assembled
91  * OPROWCOL is usually used for assemble matrices.
92  *
93  * For typical problem like Bubnov-Galerkin OPROW and OPCOL are the same. In
94  * more general case for example for non-square matrices columns and rows
95  * could have different numeration and/or different set of fields.
96  *
97  */
98  enum OpType {
99  OPROW = 1 << 0,
100  OPCOL = 1 << 1,
101  OPROWCOL = 1 << 2,
102  OPLAST = 1 << 3
103  };
104 
105  char opType;
106  std::string rowFieldName;
107  std::string colFieldName;
109 
110  /**
111  * This Constructor is used typically when some modification base shape
112  * functions on some approx. space is applied. Operator is run for all data
113  * on space.
114  *
115  * User has no access to field data from this operator.
116  */
117  UserDataOperator(const FieldSpace space, const char type = OPLAST,
118  const bool symm = true);
119 
120  UserDataOperator(const std::string &field_name, const char type,
121  const bool symm = true);
122 
123  UserDataOperator(const std::string &row_field_name,
124  const std::string &col_field_name, const char type,
125  const bool symm = true);
126 
127  /** \brief Return raw pointer to NumeredEntFiniteElement
128  */
129  inline boost::shared_ptr<const NumeredEntFiniteElement>
130  getNumeredEntFiniteElementPtr() const;
131 
132  /**
133  * \brief Return finite element entity handle
134  * @return Finite element entity handle
135  */
136  inline EntityHandle getFEEntityHandle() const;
137 
138  /**
139  * @brief Get the side number pointer
140  *
141  * \note For vertex is expection. Side basses in argument of function doWork
142  * is zero. For other entity types side can be used as argument of this
143  * function.
144  *
145  * @param side_number
146  * @param type
147  * @return boost::weak_ptr<SideNumber>
148  */
149  inline boost::weak_ptr<SideNumber> getSideNumberPtr(const int side_number,
150  const EntityType type);
151 
152  /**
153  * @brief Get the side entity
154  *
155  * \note For vertex is expection. Side basses in argument of function
156  * doWork is zero. For other entity types side can be used as argument of
157  * this function.
158  *
159  * \code
160  * MoFEMErrorCode doWork(int side, EntityType type,
161  * DataForcesAndSourcesCore::EntData &data) {
162  * MoFEMFunctionBegin;
163  *
164  * if (type == MBVERTEX) {
165  * for (int n = 0; n != number_of_nodes; ++n)
166  * EntityHandle ent = getSideEntity(n, type);
167  *
168  * // Do somthing
169  *
170  * } else {
171  * EntityHandle ent = getSideEntity(side, type);
172  *
173  * // Do somthing
174  *
175  * }
176  *
177  * MoFEMFunctionReturn(0);
178  * }
179  * \endcode
180  *
181  * @param side_number
182  * @param type
183  */
184  inline EntityHandle getSideEntity(const int side_number,
185  const EntityType type);
186 
187  /**
188  * @brief Get the number of nodes on finite element
189  *
190  * @return int
191  */
192  inline int getNumberOfNodesOnElement();
193 
194  /** \brief Get row indices
195 
196  Field could be or not declared for this element but is declared for
197  problem
198 
199  \param field_name
200  \param type entity type
201  \param side side number, any number if type is MBVERTEX
202  \return indices
203 
204  NOTE: Using those indices to assemble matrix will result in error if new
205  non-zero values need to be created.
206 
207  */
208  MoFEMErrorCode getProblemRowIndices(const std::string filed_name,
209  const EntityType type, const int side,
210  VectorInt &indices) const;
211 
212  /** \brief Get col indices
213 
214  Field could be or not declared for this element but is declared for
215  problem
216 
217  \param field_name
218  \param type entity type
219  \param side side number, any number if type is MBVERTEX
220  \return indices
221 
222  NOTE: Using those indices to assemble matrix will result in error if new
223  non-zero values need to be created.
224 
225  */
226  MoFEMErrorCode getProblemColIndices(const std::string filed_name,
227  const EntityType type, const int side,
228  VectorInt &indices) const;
229 
230  /** \brief Return raw pointer to Finite Element Method object
231  */
232  inline const FEMethod *getFEMethod() const;
233 
234  /**
235  * \brief Get operator types
236  * @return Return operator type
237  */
238  inline int getOpType() const;
239 
240  /**
241  * \brief Set operator type
242  * @param Operator type
243  */
244  inline void setOpType(const OpType type);
245 
246  /**
247  * \brief Add operator type
248  */
249  inline void addOpType(const OpType type);
250 
251  /**
252  * \brief get number of finite element in the loop
253  * @return number of finite element
254  */
255  inline int getNinTheLoop() const;
256 
257  /**
258  * \brief get size of elements in the loop
259  * @return loop size
260  */
261  inline int getLoopSize() const;
262 
263  /** \brief Get name of the element
264  */
265  inline const std::string &getFEName() const;
266 
267  /** \name Accessing KSP */
268 
269  /**@{*/
270 
271  inline const PetscData::Switches &getDataCtx() const;
272 
273  inline const KspMethod::KSPContext getKSPCtx() const;
274 
275  inline const SnesMethod::SNESContext getSNESCtx() const;
276 
277  inline const TSMethod::TSContext getTSCtx() const;
278 
279  /**@}*/
280 
281  /**@{*/
282 
283  inline Vec getKSPf() const;
284 
285  inline Mat getKSPA() const;
286 
287  inline Mat getKSPB() const;
288 
289  /**@}*/
290 
291  /** \name Accessing SNES */
292 
293  /**@{*/
294 
295  inline Vec getSNESf() const;
296 
297  inline Vec getSNESx() const;
298 
299  inline Mat getSNESA() const;
300 
301  inline Mat getSNESB() const;
302 
303  //! \deprecated Use getSNESF intead
304  DEPRECATED inline Vec getSnesF() const { return getSNESf(); }
305 
306  //! \deprecated Use getSNESX intead
307  DEPRECATED inline Vec getSnesX() const { return getSNESx(); }
308 
309  //! \deprecated Use getSNESA intead
310  DEPRECATED inline Mat getSnesA() const { return getSNESA(); }
311 
312  //! \deprecated Use getSNESB intead
313  DEPRECATED inline Mat getSnesB() const { return getSNESB(); }
314 
315  /**@}*/
316 
317  /** \name Accessing TS */
318 
319  /**@{*/
320 
321  inline Vec getTSu() const;
322 
323  inline Vec getTSu_t() const;
324 
325  inline Vec getTSu_tt() const;
326 
327  inline Vec getTSf() const;
328 
329  inline Mat getTSA() const;
330 
331  inline Mat getTSB() const;
332 
333  inline int getTSstep() const;
334 
335  inline double getTStime() const;
336 
337  inline double getTSa() const;
338 
339  /**@}*/
340 
341  /**@{*/
342 
343  /** \name Base funtions and integration points */
344 
345  /** \brief matrix of integration (Gauss) points for Volume Element
346  *
347  * For triangle: columns 0,1 are x,y coordinates respectively and column
348  * 2 is a weight value for example getGaussPts()(1,13) returns y
349  * coordinate of 13th Gauss point on particular volume element
350  *
351  * For tetrahedron: columns 0,1,2 are x,y,z coordinates respectively and
352  * column 3 is a weight value for example getGaussPts()(1,13) returns y
353  * coordinate of 13th Gauss point on particular volume element
354  *
355  */
356  inline MatrixDouble &getGaussPts();
357 
358  /**
359  * @brief Get integration weights
360  *
361  * \code
362  * auto t_w = getFTensor0IntegrationWeight();
363  * for(int gg = 0; gg!=getGaussPts.size2(); ++gg) {
364  * // integrate something
365  * ++t_w;
366  * }
367  * \endcode
368  *
369  * @return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>
370  */
371  inline auto getFTensor0IntegrationWeight();
372 
373  /**@}*/
374 
375  /**@{*/
376 
377  /** \name Deprecated (do not use) */
378 
379  // \deprecated Deprecated function with spelling mistake
381  getPorblemRowIndices(const std::string filed_name, const EntityType type,
382  const int side, VectorInt &indices) const;
383 
384  /**@}*/
385 
386  protected:
388 
389  virtual MoFEMErrorCode setPtrFE(ForcesAndSourcesCore *ptr);
390 
391  inline ForcesAndSourcesCore *getPtrFE() const;
392 
393  inline ForcesAndSourcesCore *getSidePtrFE() const;
394 
395  private:
396  /**
397  * @brief User call this function to loop over elements on the side of
398  * face. This function calls finite element with is operator to do
399  * calculations.
400  *
401  * @param fe_name name of the side element
402  * @param side_fe pointer to the side element instance
403  * @param dim dimension the of side element
404  * @param ent_for_side entity handle for which adjacent volume or face will
405  * be accessed
406  * @return MoFEMErrorCode
407  */
408  MoFEMErrorCode loopSide(const string &fe_name,
409  ForcesAndSourcesCore *side_fe, const size_t dim,
410  const EntityHandle ent_for_side = 0);
411 
412  friend class ForcesAndSourcesCore;
416  };
417 
418  /** \brief Use to push back operator for row operator
419 
420  It can be used to calculate nodal forces or other quantities on the mesh.
421 
422  */
423  boost::ptr_vector<UserDataOperator> &getOpPtrVector() { return opPtrVector; }
424 
425  /**
426  * @brief Get the Entity Polynomial Base object
427  *
428  * @return boost::shared_ptr<BaseFunction>&&
429  */
430  auto &getElementPolynomialBase() { return elementPolynomialBasePtr; }
431 
432  /**
433  * @brief Get the User Polynomial Base object
434  *
435  * @return boost::shared_ptr<BaseFunction>&
436  */
437  auto &getUserPolynomialBase() { return userPolynomialBasePtr; }
438 
439  /**
440  * @brief Matrix of integration points
441  *
442  * Columns is equal to number of integration points, numver of rows
443  * depends on dimension of finite element entity, for example for
444  * tetrahedron rows are x,y,z,weight. Last row is integration weight.
445  *
446  * FIXME: that should be moved to private class data and acessed only by
447  * member function
448  */
450 
451  virtual MoFEMErrorCode preProcess();
452  virtual MoFEMErrorCode operator()();
453  virtual MoFEMErrorCode postProcess();
454 
455 public:
456  /** \brief Get max order of approximation for data fields
457 
458  Method getMaxDataOrder () return maximal order on entities, for
459  all data on the element. So for example if finite element is triangle, and
460  triangle base function have order 4 and on edges base function have order
461  2, this function return 4.
462 
463  If finite element has for example 2 or more approximated fields, for
464  example Pressure (order 3) and displacement field (order 5), this function
465  returns 5.
466 
467  */
468  int getMaxDataOrder() const;
469 
470  /// \brief Get max order of approximation for field in rows
471  int getMaxRowOrder() const;
472 
473  /// \brief Get max order of approximation for field in columns
474  int getMaxColOrder() const;
475 
476  /**
477  * @brief Get the entity data
478  *
479  * @param space
480  * @param type
481  * @param side
482  * @return const DataForcesAndSourcesCore::EntData&
483  */
485  const EntityType type,
486  const int side) const {
487  return dataOnElement[space]->dataOnEntities[type][side];
488  }
489 
490  /**
491  * @brief Get the entity data
492  *
493  * @param space
494  * @param type
495  * @param side
496  * @return DataForcesAndSourcesCore::EntData&
497  */
499  getEntData(const FieldSpace space, const EntityType type, const int side) {
500  return dataOnElement[space]->dataOnEntities[type][side];
501  }
502 
503 protected:
504  /**
505  * \brief get sense (orientation) of entity
506  * @param type type of entity
507  * @param data entity data
508  * @return error code
509  */
510  MoFEMErrorCode getEntitySense(
511  const EntityType type,
512  boost::ptr_vector<DataForcesAndSourcesCore::EntData> &data) const;
513 
514  /**
515  * @brief Get the entity data order
516  *
517  * @param type
518  * @param space
519  * @param data
520  * @return MoFEMErrorCode
521  */
522  MoFEMErrorCode getEntityDataOrder(
523  const EntityType type, const FieldSpace space,
524  boost::ptr_vector<DataForcesAndSourcesCore::EntData> &data) const;
525 
526  /**
527  * @brief Get the entity sense (orientation)
528  *
529  * @tparam type
530  * @param data
531  * @return MoFEMErrorCode
532  */
533  template <EntityType type>
535  return getEntitySense(type, data.dataOnEntities[type]);
536  }
537 
538  /**
539  * @brief Get the entity data order for given space
540  *
541  * @tparam type
542  * @param data
543  * @param space
544  * @return MoFEMErrorCode
545  */
546  template <EntityType type>
548  const FieldSpace space) const {
549  return getEntityDataOrder(type, space, data.dataOnEntities[type]);
550  }
551 
552  /** \name Indices */
553 
554  /**@{*/
555 
556  /// \brief get node indices
557  MoFEMErrorCode getNodesIndices(const boost::string_ref field_name,
559  VectorInt &nodes_indices,
560  VectorInt &local_nodes_indices) const;
561 
562  /// \brief get row node indices from FENumeredDofEntity_multiIndex
563  MoFEMErrorCode getRowNodesIndices(DataForcesAndSourcesCore &data,
564  const std::string &field_name) const;
565 
566  /// \brief get col node indices from FENumeredDofEntity_multiIndex
567  MoFEMErrorCode getColNodesIndices(DataForcesAndSourcesCore &data,
568  const std::string &field_name) const;
569 
570  MoFEMErrorCode getEntityIndices(
571  DataForcesAndSourcesCore &data, const std::string &field_name,
572  FENumeredDofEntity_multiIndex &dofs, const EntityType type_lo = MBVERTEX,
573  const EntityType type_hi = MBPOLYHEDRON) const;
574 
575  inline MoFEMErrorCode
576  getEntityRowIndices(DataForcesAndSourcesCore &data,
577  const std::string &field_name,
578  const EntityType type_lo = MBVERTEX,
579  const EntityType type_hi = MBPOLYHEDRON) const;
580 
581  inline MoFEMErrorCode
582  getEntityColIndices(DataForcesAndSourcesCore &data,
583  const std::string &field_name,
584  const EntityType type_lo = MBVERTEX,
585  const EntityType type_hi = MBPOLYHEDRON) const;
586 
587  /// \brief get NoField indices
588  MoFEMErrorCode getNoFieldIndices(const std::string &field_name,
590  VectorInt &nodes_indices) const;
591 
592  /// \brief get col NoField indices
593  MoFEMErrorCode getNoFieldRowIndices(DataForcesAndSourcesCore &data,
594  const std::string &field_name) const;
595 
596  /// \brief get col NoField indices
597  MoFEMErrorCode getNoFieldColIndices(DataForcesAndSourcesCore &data,
598  const std::string &field_name) const;
599 
600  /**@}*/
601 
602  /** \name Data */
603 
604  /**@{*/
605 
606  /**
607  * \brief Get field data on nodes
608  */
609  MoFEMErrorCode getNoFieldFieldData(const boost::string_ref field_name,
611  VectorDouble &ent_field_data,
612  VectorDofs &ent_field_dofs) const;
613 
614  MoFEMErrorCode getNoFieldFieldData(DataForcesAndSourcesCore &data,
615  const boost::string_ref field_name) const;
616 
617  /**
618  * \brief Get data on nodes
619  * @param data Data structure
620  * @param field_name Field name
621  * @return Error code
622  */
623  MoFEMErrorCode getNodesFieldData(DataForcesAndSourcesCore &data,
624  const std::string &field_name) const;
625 
627  getEntityFieldData(DataForcesAndSourcesCore &data,
628  const std::string &field_name,
629  const EntityType type_lo = MBVERTEX,
630  const EntityType type_hi = MBPOLYHEDRON) const;
631 
632  /**@}*/
633 
634  /// \brief Get nodes on triangles
635  MoFEMErrorCode getFaceTriNodes(DataForcesAndSourcesCore &data) const;
636 
637  /// \brief Get field approximation space and base on entities
639  getSpacesAndBaseOnEntities(DataForcesAndSourcesCore &data) const;
640 
641  /** \name Data form NumeredDofEntity_multiIndex */
642 
643  /**@{*/
644 
645  /// \brief get indices of nodal indices which are declared for problem but
646  /// not this particular element
647  MoFEMErrorCode getProblemNodesIndices(const std::string &field_name,
648  const NumeredDofEntity_multiIndex &dofs,
649  VectorInt &nodes_indices) const;
650 
651  /// \brief get indices by type (generic function) which are declared for
652  /// problem but not this particular element
653  MoFEMErrorCode getProblemTypeIndices(const std::string &field_name,
654  const NumeredDofEntity_multiIndex &dofs,
655  EntityType type, int side_number,
656  VectorInt &indices) const;
657 
658  MoFEMErrorCode getProblemNodesRowIndices(const std::string &field_name,
659  VectorInt &nodes_indices) const;
660  MoFEMErrorCode getProblemTypeRowIndices(const std::string &field_name,
661  EntityType type, int side_number,
662  VectorInt &indices) const;
663  MoFEMErrorCode getProblemNodesColIndices(const std::string &field_name,
664  VectorInt &nodes_indices) const;
665  MoFEMErrorCode getProblemTypeColIndices(const std::string &field_name,
666  EntityType type, int side_number,
667  VectorInt &indices) const;
668 
669  /**@}*/
670 
671  /**
672  * \brief another variant of getRule
673  * @param order_row order of base function on row
674  * @param order_col order of base function on columns
675  * @param order_data order of base function approximating data
676  * @return integration rule
677  *
678  * This function is overloaded by the user. The integration rule
679  * is set such that specific operator implemented by the user is
680  * integrated accurately. For example if user implement bilinear operator
681  * \f[
682  * b(u,v) =
683  * \int_\mathcal{T}
684  * \frac{\partial u_i}{\partial x_j}\frac{\partial v_i}{\partial x_j}
685  * \textrm{d}\mathcal{T}
686  * \f]
687  * then if \f$u\f$ and \f$v\f$ are polynomial of given \em order, then
688  * exact integral would be \code int getRule(int order) { return
689  * 2*(order-1); }; \endcode
690  *
691  * The integration points and weights are set appropriately for given
692  * entity type and integration rule from \ref quad.c
693  *
694  * Method \ref ForcesAndSourcesCore::getRule takes at argument takes
695  * maximal polynomial order set on the element on all fields defined on
696  * the element. If a user likes to have more control, another variant of
697  * this function can be called which distinguishing between field orders
698  * on rows, columns and data, the i.e. first argument of a bilinear form,
699  * the second argument of bilinear form and field coefficients on the
700  * element.
701  *
702  * \note If user set rule to -1 or any other negative integer, then method
703  * \ref ForcesAndSourcesCore::setGaussPts is called. In that method user
704  * can implement own (specific) integration method.
705  *
706  * \bug this function should be const
707  */
708  virtual int getRule(int order_row, int order_col, int order_data);
709 
710  /** \brief set user specific integration rule
711 
712  This function allows for user defined integration rule. The key is to
713  called matrix gaussPts, which is used by other MoFEM procedures. Matrix
714  has number of rows equal to problem dimension plus one, where last index
715  is used to store weight values. %Number of columns is equal to number of
716  integration points.
717 
718  \note This function is called if method \ref
719  ForcesAndSourcesCore::getRule is returning integer -1 or any other
720  negative integer.
721 
722  User sets
723  \code
724  MatrixDouble gaussPts;
725  \endcode
726  where
727  \code
728  gaussPts.resize(dim+1,nb_gauss_pts);
729  \endcode
730  number rows represents local coordinates of integration points
731  in reference element, where last index in row is for integration weight.
732 
733  */
734  virtual MoFEMErrorCode setGaussPts(int order_row, int order_col,
735  int order_data);
736 
737  /**
738  * \brief Calculate base functions
739  * @return Error code
740  */
742  calHierarchicalBaseFunctionsOnElement(const FieldApproximationBase b);
743 
744  /**
745  * \brief Calculate base functions
746  * @return Error code
747  */
748  MoFEMErrorCode calHierarchicalBaseFunctionsOnElement();
749 
750  /**
751  * @brief Calculate Bernstein-Bezier base
752  *
753  * @return MoFEMErrorCode
754  */
755  MoFEMErrorCode calBernsteinBezierBaseFunctionsOnElement();
756 
757  /**
758  * @brief Create a entity data on element object
759  *
760  * @return MoFEMErrorCode
761  */
762  MoFEMErrorCode createDataOnElement();
763 
764  /**
765  * @brief Iterate user data operators
766  *
767  * @return MoFEMErrorCode
768  */
769  MoFEMErrorCode loopOverOperators();
770 
771  /**@{*/
772 
773  /** \name Deprecated (do not use) */
774 
775  /** \deprecated Use getRule(int row_order, int col_order, int data order)
776  */
777  virtual int getRule(int order);
778 
779  /** \deprecated setGaussPts(int row_order, int col_order, int data order);
780  */
781  virtual MoFEMErrorCode setGaussPts(int order);
782 
783  /**@/}*/
784 
785  /**
786  * @brief Entity data on element entity rows fields
787  *
788  */
789  const std::array<boost::shared_ptr<DataForcesAndSourcesCore>, LASTSPACE>
791 
792  /**
793  * @brief Entity data on element entity columns fields
794  *
795  */
796  const std::array<boost::shared_ptr<DataForcesAndSourcesCore>, LASTSPACE>
798 
804 
805  /**
806  * @brief Vector of finite element users data operators
807  *
808  */
809  boost::ptr_vector<UserDataOperator> opPtrVector;
810 
811  friend class UserDataOperator;
812 
813 protected:
814  /**
815  * @brief Last evaluated type of element entity
816  *
817  */
819 
820 private:
821  /**
822  * @brief Pointer to entity polynomial base
823  *
824  */
825  boost::shared_ptr<BaseFunction> elementPolynomialBasePtr;
826 
827  /**
828  * @brief Pointer to user polynomail base
829  */
830  boost::shared_ptr<BaseFunction> userPolynomialBasePtr;
831 
832  /**
833  * @brief Element to integrate on the sides
834  *
835  */
837 
838  /**
839  * @brief Set the pointer to face element on the side
840  *
841  * \note Function is is used by face element, while it iterates over
842  * elements on the side
843  *
844  * @param side_fe_ptr
845  * @return MoFEMErrorCode
846  */
847  MoFEMErrorCode setSideFEPtr(const ForcesAndSourcesCore *side_fe_ptr);
848 
852 };
853 
854 /// \deprecated Used ForcesAndSourcesCore instead
856 
857 MoFEMErrorCode ForcesAndSourcesCore::getEntityRowIndices(
858  DataForcesAndSourcesCore &data, const std::string &field_name,
859  const EntityType type_lo, const EntityType type_hi) const {
860  return getEntityIndices(data, field_name,
861  const_cast<FENumeredDofEntity_multiIndex &>(
862  numeredEntFiniteElementPtr->getRowsDofs()),
863  type_lo, type_hi);
864 }
865 
866 MoFEMErrorCode ForcesAndSourcesCore::getEntityColIndices(
867  DataForcesAndSourcesCore &data, const std::string &field_name,
868  const EntityType type_lo, const EntityType type_hi) const {
869  return getEntityIndices(data, field_name,
870  const_cast<FENumeredDofEntity_multiIndex &>(
871  numeredEntFiniteElementPtr->getColsDofs()),
872  type_lo, type_hi);
873 }
874 
875 boost::shared_ptr<const NumeredEntFiniteElement>
876 ForcesAndSourcesCore::UserDataOperator::getNumeredEntFiniteElementPtr() const {
877  return ptrFE->numeredEntFiniteElementPtr;
878 };
879 
880 EntityHandle ForcesAndSourcesCore::UserDataOperator::getFEEntityHandle() const {
881  return getNumeredEntFiniteElementPtr()->getEnt();
882 }
883 
884 boost::weak_ptr<SideNumber>
885 ForcesAndSourcesCore::UserDataOperator::getSideNumberPtr(
886  const int side_number, const EntityType type) {
887  auto &side_table_by_side_and_type =
888  ptrFE->numeredEntFiniteElementPtr->getSideNumberTable().get<1>();
889  auto side_it =
890  side_table_by_side_and_type.find(boost::make_tuple(type, side_number));
891  if (side_it != side_table_by_side_and_type.end())
892  return *side_it;
893  else
894  return boost::weak_ptr<SideNumber>();
895 }
896 
898 ForcesAndSourcesCore::UserDataOperator::getSideEntity(const int side_number,
899  const EntityType type) {
900  if (auto side_ptr = getSideNumberPtr(side_number, type).lock())
901  return side_ptr->ent;
902  else
903  return 0;
904 }
905 
906 int ForcesAndSourcesCore::UserDataOperator::getNumberOfNodesOnElement() {
907  int num_nodes;
908  CHKERR ptrFE->getNumberOfNodes(num_nodes);
909  return num_nodes;
910 }
911 
912 const FEMethod *ForcesAndSourcesCore::UserDataOperator::getFEMethod() const {
913  return ptrFE;
914 }
915 
916 int ForcesAndSourcesCore::UserDataOperator::getOpType() const { return opType; }
917 
918 void ForcesAndSourcesCore::UserDataOperator::setOpType(const OpType type) {
919  opType = type;
920 }
921 
922 void ForcesAndSourcesCore::UserDataOperator::addOpType(const OpType type) {
923  opType |= type;
924 }
925 
926 int ForcesAndSourcesCore::UserDataOperator::getNinTheLoop() const {
927  return getFEMethod()->getNinTheLoop();
928 }
929 
930 int ForcesAndSourcesCore::UserDataOperator::getLoopSize() const {
931  return getFEMethod()->getLoopSize();
932 }
933 
934 const std::string &ForcesAndSourcesCore::UserDataOperator::getFEName() const {
935  return getFEMethod()->feName;
936 }
937 
938 const PetscData::Switches &
939 ForcesAndSourcesCore::UserDataOperator::getDataCtx() const {
940  return getFEMethod()->data_ctx;
941 }
942 
944 ForcesAndSourcesCore::UserDataOperator::getKSPCtx() const {
945  return getFEMethod()->ksp_ctx;
946 }
947 
949 ForcesAndSourcesCore::UserDataOperator::getSNESCtx() const {
950  return getFEMethod()->snes_ctx;
951 }
952 
954 ForcesAndSourcesCore::UserDataOperator::getTSCtx() const {
955  return getFEMethod()->ts_ctx;
956 }
957 
958 Vec ForcesAndSourcesCore::UserDataOperator::getKSPf() const {
959  return getFEMethod()->ksp_f;
960 }
961 
962 Mat ForcesAndSourcesCore::UserDataOperator::getKSPA() const {
963  return getFEMethod()->ksp_A;
964 }
965 
966 Mat ForcesAndSourcesCore::UserDataOperator::getKSPB() const {
967  return getFEMethod()->ksp_B;
968 }
969 
970 Vec ForcesAndSourcesCore::UserDataOperator::getSNESf() const {
971  return getFEMethod()->snes_f;
972 }
973 
974 Vec ForcesAndSourcesCore::UserDataOperator::getSNESx() const {
975  return getFEMethod()->snes_x;
976 }
977 
978 Mat ForcesAndSourcesCore::UserDataOperator::getSNESA() const {
979  return getFEMethod()->snes_A;
980 }
981 
982 Mat ForcesAndSourcesCore::UserDataOperator::getSNESB() const {
983  return getFEMethod()->snes_B;
984 }
985 
986 Vec ForcesAndSourcesCore::UserDataOperator::getTSu() const {
987  return getFEMethod()->ts_u;
988 }
989 
990 Vec ForcesAndSourcesCore::UserDataOperator::getTSu_t() const {
991  return getFEMethod()->ts_u_t;
992 }
993 
994 Vec ForcesAndSourcesCore::UserDataOperator::getTSu_tt() const {
995  return getFEMethod()->ts_u_tt;
996 }
997 
998 Vec ForcesAndSourcesCore::UserDataOperator::getTSf() const {
999  return getFEMethod()->ts_F;
1000 }
1001 
1002 Mat ForcesAndSourcesCore::UserDataOperator::getTSA() const {
1003  return getFEMethod()->ts_A;
1004 }
1005 
1006 Mat ForcesAndSourcesCore::UserDataOperator::getTSB() const {
1007  return getFEMethod()->ts_B;
1008 }
1009 
1010 int ForcesAndSourcesCore::UserDataOperator::getTSstep() const {
1011  return getFEMethod()->ts_step;
1012 }
1013 
1014 double ForcesAndSourcesCore::UserDataOperator::getTStime() const {
1015  return getFEMethod()->ts_t;
1016 }
1017 
1018 double ForcesAndSourcesCore::UserDataOperator::getTSa() const {
1019  return getFEMethod()->ts_a;
1020 }
1021 
1022 MatrixDouble &ForcesAndSourcesCore::UserDataOperator::getGaussPts() {
1023  return static_cast<ForcesAndSourcesCore *>(ptrFE)->gaussPts;
1024 }
1025 
1026 auto ForcesAndSourcesCore::UserDataOperator::getFTensor0IntegrationWeight() {
1028  &(getGaussPts()(getGaussPts().size1() - 1, 0)));
1029 }
1030 
1031 MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::getPorblemRowIndices(
1032  const std::string filed_name, const EntityType type, const int side,
1033  VectorInt &indices) const {
1034  return getProblemRowIndices(filed_name, type, side, indices);
1035 }
1036 
1037 ForcesAndSourcesCore *ForcesAndSourcesCore::UserDataOperator::getPtrFE() const {
1038  return ptrFE;
1039 }
1040 
1042 ForcesAndSourcesCore::UserDataOperator::getSidePtrFE() const {
1043  return ptrFE->sidePtrFE;
1044 }
1045 
1046 } // namespace MoFEM
1047 
1048 #endif //__FORCES_AND_SOURCES_CORE__HPP__
1049 
1050 /**
1051  * \defgroup mofem_forces_and_sources Forces and sources
1052  * \ingroup mofem
1053  *
1054  * \brief Manages complexities related to assembly of vector and matrices at
1055  * single finite element level.
1056  *
1057  **/
MoFEMErrorCode getEntityDataOrder(DataForcesAndSourcesCore &data, const FieldSpace space) const
Get the entity data order for given space.
DataForcesAndSourcesCore & dataL2
structure for User Loop Methods on finite elementsIt can be used to calculate stiffness matrices,...
Deprecated interface functions.
boost::ptr_vector< UserDataOperator > opPtrVector
Vector of finite element users data operators.
MatrixDouble gaussPts
Matrix of integration points.
base operator to do operations at Gauss Pt. level
ForcesAndSourcesCore * sidePtrFE
Element to integrate on the sides.
data structure for finite element entityIt keeps that about indices of degrees of freedom,...
DataForcesAndSourcesCore & dataH1
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
std::bitset< 8 > Switches
Definition: LoopMethods.hpp:59
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
auto & getUserPolynomialBase()
Get the User Polynomial Base object.
multi_index_container< boost::shared_ptr< FENumeredDofEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< FENumeredDofEntity::interface_type_DofEntity, const UId &, &FENumeredDofEntity::getGlobalUniqueId > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< FENumeredDofEntity::interface_type_DofEntity, EntityHandle, &FENumeredDofEntity::getEnt > >, ordered_non_unique< tag< FieldName_mi_tag >, const_mem_fun< FENumeredDofEntity::interface_type_Field, boost::string_ref, &FENumeredDofEntity::getNameRef > >, ordered_non_unique< tag< Composite_Name_Type_And_Side_Number_mi_tag >, composite_key< FENumeredDofEntity, const_mem_fun< FENumeredDofEntity::interface_type_Field, boost::string_ref, &FENumeredDofEntity::getNameRef >, const_mem_fun< FENumeredDofEntity::interface_type_RefEntity, EntityType, &FENumeredDofEntity::getEntType >, KeyFromKey< member< SideNumber, char, &SideNumber::side_number >, member< FENumeredDofEntity::BaseFEEntity, boost::shared_ptr< SideNumber >, &FENumeredDofEntity::sideNumberPtr > > > >, ordered_non_unique< tag< Composite_Name_And_Type_mi_tag >, composite_key< FENumeredDofEntity, const_mem_fun< FENumeredDofEntity::interface_type_Field, boost::string_ref, &FENumeredDofEntity::getNameRef >, const_mem_fun< FENumeredDofEntity::interface_type_RefEntity, EntityType, &FENumeredDofEntity::getEntType > > >, ordered_non_unique< tag< Composite_Name_And_Ent_mi_tag >, composite_key< FENumeredDofEntity, const_mem_fun< FENumeredDofEntity::interface_type_Field, boost::string_ref, &FENumeredDofEntity::getNameRef >, const_mem_fun< FENumeredDofEntity::interface_type_DofEntity, EntityHandle, &FENumeredDofEntity::getEnt > > > > > FENumeredDofEntity_multiIndex
MultiIndex container keeps FENumeredDofEntity.
boost::shared_ptr< BaseFunction > userPolynomialBasePtr
Pointer to user polynomail base.
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
OpType
Controls loop over entities on element.
FieldSpace in [ 0, LASTSPACE )
Definition: definitions.h:181
Edge finite elementUser is implementing own operator at Gauss points level, by own object derived fro...
RuleHookFun getRuleHook
Hook to get rule.
const int dim
DataForcesAndSourcesCore::EntData & getEntData(const FieldSpace space, const EntityType type, const int side)
Get the entity data.
DataForcesAndSourcesCore & dataHcurl
EntityType lastEvaluatedElementEntityType
Last evaluated type of element entity.
const std::array< boost::shared_ptr< DataForcesAndSourcesCore >, LASTSPACE > derivedDataOnElement
Entity data on element entity columns fields.
ForcesAndSourcesCore::UserDataOperator UserDataOperator
FieldApproximationBase
approximation base
Definition: definitions.h:150
boost::shared_ptr< BaseFunction > elementPolynomialBasePtr
Pointer to entity polynomial base.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
KSPContext
pass information about context of KSP/DM for with finite element is computed
Definition: LoopMethods.hpp:96
boost::function< int(int order_row, int order_col, int order_data)> RuleHookFun
MoFEMErrorCode getEntitySense(DataForcesAndSourcesCore &data) const
Get the entity sense (orientation)
FieldSpace
approximation spaces
Definition: definitions.h:174
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
#define CHKERR
Inline error check.
Definition: definitions.h:602
constexpr int order
multi_index_container< boost::shared_ptr< FEDofEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< FEDofEntity::interface_type_DofEntity, const UId &, &FEDofEntity::getGlobalUniqueId > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< FEDofEntity::interface_type_DofEntity, EntityHandle, &FEDofEntity::getEnt > >, ordered_non_unique< tag< FieldName_mi_tag >, const_mem_fun< FEDofEntity::interface_type_Field, boost::string_ref, &FEDofEntity::getNameRef > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< FEDofEntity::interface_type_RefEntity, EntityType, &FEDofEntity::getEntType > >, ordered_non_unique< tag< Composite_Name_And_Type_mi_tag >, composite_key< FEDofEntity, const_mem_fun< FEDofEntity::interface_type_Field, boost::string_ref, &FEDofEntity::getNameRef >, const_mem_fun< FEDofEntity::interface_type_RefEntity, EntityType, &FEDofEntity::getEntType > > >, ordered_non_unique< tag< Composite_Name_And_Ent_mi_tag >, composite_key< FEDofEntity, const_mem_fun< FEDofEntity::interface_type_Field, boost::string_ref, &FEDofEntity::getNameRef >, const_mem_fun< FEDofEntity::interface_type_DofEntity, EntityHandle, &FEDofEntity::getEnt > > >, ordered_non_unique< tag< Composite_Name_Type_And_Side_Number_mi_tag >, composite_key< FEDofEntity, const_mem_fun< FEDofEntity::interface_type_Field, boost::string_ref, &FEDofEntity::getNameRef >, const_mem_fun< FEDofEntity::interface_type_RefEntity, EntityType, &FEDofEntity::getEntType >, KeyFromKey< member< SideNumber, char, &SideNumber::side_number >, member< FEDofEntity::BaseFEEntity, boost::shared_ptr< SideNumber >, &FEDofEntity::sideNumberPtr > > > > > > FEDofEntity_multiIndex
MultiIndex container keeps FEDofEntity.
const DataForcesAndSourcesCore::EntData & getEntData(const FieldSpace space, const EntityType type, const int side) const
Get the entity data.
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:72
ContactPrism finite elementUser is implementing own operator at Gauss points level,...
#define DEPRECATED
Definition: definitions.h:29
DataForcesAndSourcesCore & dataHdiv
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
Data on single entity (This is passed as argument to DataOperator::doWork)
structure to get information form mofem into DataForcesAndSourcesCore
Base volume element used to integrate on contact surface (could be extended to other volume elements)...
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:73
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...
const std::array< boost::shared_ptr< DataForcesAndSourcesCore >, LASTSPACE > dataOnElement
Entity data on element entity rows fields.
RuleHookFun setRuleHook
Set function to calculate integration rule.
auto & getElementPolynomialBase()
Get the Entity Polynomial Base object.
multi_index_container< boost::shared_ptr< NumeredDofEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, const UId &, &NumeredDofEntity::getGlobalUniqueId > >, ordered_non_unique< tag< Unique_Ent_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, const UId &, &NumeredDofEntity::getEntGlobalUniqueId > >, ordered_non_unique< tag< Part_mi_tag >, member< NumeredDofEntity, unsigned int, &NumeredDofEntity::pArt > >, ordered_non_unique< tag< Idx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::dofIdx > >, ordered_non_unique< tag< FieldName_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_Field, boost::string_ref, &NumeredDofEntity::getNameRef > >, ordered_non_unique< tag< PetscGlobalIdx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::petscGloablDofIdx > >, ordered_non_unique< tag< PetscLocalIdx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::petscLocalDofIdx > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, EntityHandle, &NumeredDofEntity::getEnt > >, ordered_non_unique< tag< Composite_Name_And_Ent_And_EntDofIdx_mi_tag >, composite_key< NumeredDofEntity, const_mem_fun< NumeredDofEntity::interface_type_Field, boost::string_ref, &NumeredDofEntity::getNameRef >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, EntityHandle, &NumeredDofEntity::getEnt >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, DofIdx, &NumeredDofEntity::getEntDofIdx > > >, ordered_non_unique< tag< Composite_Name_And_Part_mi_tag >, composite_key< NumeredDofEntity, const_mem_fun< NumeredDofEntity::interface_type_Field, boost::string_ref, &NumeredDofEntity::getNameRef >, member< NumeredDofEntity, unsigned int, &NumeredDofEntity::pArt > > >, ordered_non_unique< tag< Composite_Name_Ent_And_Part_mi_tag >, composite_key< NumeredDofEntity, const_mem_fun< NumeredDofEntity::interface_type_Field, boost::string_ref, &NumeredDofEntity::getNameRef >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, EntityHandle, &NumeredDofEntity::getEnt >, member< NumeredDofEntity, unsigned int, &NumeredDofEntity::pArt > > > > > NumeredDofEntity_multiIndex
MultiIndex container keeps NumeredDofEntity.
Face finite elementUser is implementing own operator at Gauss point level, by own object derived from...
DEPRECATED typedef ForcesAndSourcesCore ForcesAndSurcesCore
DataForcesAndSourcesCore & dataNoField