v0.9.1
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
402  * @param side_fe
403  * @param dim
404  * @return MoFEMErrorCode
405  */
406  MoFEMErrorCode loopSide(const string &fe_name,
407  ForcesAndSourcesCore *side_fe, const size_t dim);
408 
409  friend class ForcesAndSourcesCore;
413  };
414 
415  /** \brief Use to push back operator for row operator
416 
417  It can be used to calculate nodal forces or other quantities on the mesh.
418 
419  */
420  boost::ptr_vector<UserDataOperator> &getOpPtrVector() { return opPtrVector; }
421 
422  /**
423  * @brief Get the Entity Polynomial Base object
424  *
425  * @return boost::shared_ptr<BaseFunction>&&
426  */
427  auto &getElementPolynomialBase() { return elementPolynomialBasePtr; }
428 
429  /**
430  * @brief Get the User Polynomial Base object
431  *
432  * @return boost::shared_ptr<BaseFunction>&
433  */
434  auto &getUserPolynomialBase() { return userPolynomialBasePtr; }
435 
436  /**
437  * @brief Matrix of integration points
438  *
439  * Columns is equal to number of integration points, numver of rows
440  * depends on dimension of finite element entity, for example for
441  * tetrahedron rows are x,y,z,weight. Last row is integration weight.
442  *
443  * FIXME: that should be moved to private class data and acessed only by
444  * member function
445  */
447 
448  virtual MoFEMErrorCode preProcess();
449  virtual MoFEMErrorCode operator()();
450  virtual MoFEMErrorCode postProcess();
451 
452 public:
453  /** \brief Get max order of approximation for data fields
454 
455  Method getMaxDataOrder () return maximal order on entities, for
456  all data on the element. So for example if finite element is triangle, and
457  triangle base function have order 4 and on edges base function have order
458  2, this function return 4.
459 
460  If finite element has for example 2 or more approximated fields, for
461  example Pressure (order 3) and displacement field (order 5), this function
462  returns 5.
463 
464  */
465  int getMaxDataOrder() const;
466 
467  /// \brief Get max order of approximation for field in rows
468  int getMaxRowOrder() const;
469 
470  /// \brief Get max order of approximation for field in columns
471  int getMaxColOrder() const;
472 
473  /// \brief Get number of DOFs on element
474  MoFEMErrorCode getNumberOfNodes(int &num_nodes) 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 
851 };
852 
853 /// \deprecated Used ForcesAndSourcesCore instead
855 
856 MoFEMErrorCode ForcesAndSourcesCore::getEntityRowIndices(
857  DataForcesAndSourcesCore &data, const std::string &field_name,
858  const EntityType type_lo, const EntityType type_hi) const {
859  return getEntityIndices(data, field_name,
860  const_cast<FENumeredDofEntity_multiIndex &>(
861  numeredEntFiniteElementPtr->getRowsDofs()),
862  type_lo, type_hi);
863 }
864 
865 MoFEMErrorCode ForcesAndSourcesCore::getEntityColIndices(
866  DataForcesAndSourcesCore &data, const std::string &field_name,
867  const EntityType type_lo, const EntityType type_hi) const {
868  return getEntityIndices(data, field_name,
869  const_cast<FENumeredDofEntity_multiIndex &>(
870  numeredEntFiniteElementPtr->getColsDofs()),
871  type_lo, type_hi);
872 }
873 
874 boost::shared_ptr<const NumeredEntFiniteElement>
875 ForcesAndSourcesCore::UserDataOperator::getNumeredEntFiniteElementPtr() const {
876  return ptrFE->numeredEntFiniteElementPtr;
877 };
878 
879 EntityHandle ForcesAndSourcesCore::UserDataOperator::getFEEntityHandle() const {
880  return getNumeredEntFiniteElementPtr()->getEnt();
881 }
882 
883 boost::weak_ptr<SideNumber>
884 ForcesAndSourcesCore::UserDataOperator::getSideNumberPtr(
885  const int side_number, const EntityType type) {
886  auto &side_table_by_side_and_type =
887  ptrFE->numeredEntFiniteElementPtr->getSideNumberTable().get<1>();
888  auto side_it =
889  side_table_by_side_and_type.find(boost::make_tuple(type, side_number));
890  if (side_it != side_table_by_side_and_type.end())
891  return *side_it;
892  else
893  return boost::weak_ptr<SideNumber>();
894 }
895 
897 ForcesAndSourcesCore::UserDataOperator::getSideEntity(const int side_number,
898  const EntityType type) {
899  if (auto side_ptr = getSideNumberPtr(side_number, type).lock())
900  return side_ptr->ent;
901  else
902  return 0;
903 }
904 
905 int ForcesAndSourcesCore::UserDataOperator::getNumberOfNodesOnElement() {
906  int num_nodes;
907  CHKERR ptrFE->getNumberOfNodes(num_nodes);
908  return num_nodes;
909 }
910 
911 const FEMethod *ForcesAndSourcesCore::UserDataOperator::getFEMethod() const {
912  return ptrFE;
913 }
914 
915 int ForcesAndSourcesCore::UserDataOperator::getOpType() const { return opType; }
916 
917 void ForcesAndSourcesCore::UserDataOperator::setOpType(const OpType type) {
918  opType = type;
919 }
920 
921 void ForcesAndSourcesCore::UserDataOperator::addOpType(const OpType type) {
922  opType |= type;
923 }
924 
925 int ForcesAndSourcesCore::UserDataOperator::getNinTheLoop() const {
926  return getFEMethod()->getNinTheLoop();
927 }
928 
929 int ForcesAndSourcesCore::UserDataOperator::getLoopSize() const {
930  return getFEMethod()->getLoopSize();
931 }
932 
933 const std::string &ForcesAndSourcesCore::UserDataOperator::getFEName() const {
934  return getFEMethod()->feName;
935 }
936 
937 const PetscData::Switches &
938 ForcesAndSourcesCore::UserDataOperator::getDataCtx() const {
939  return getFEMethod()->data_ctx;
940 }
941 
943 ForcesAndSourcesCore::UserDataOperator::getKSPCtx() const {
944  return getFEMethod()->ksp_ctx;
945 }
946 
948 ForcesAndSourcesCore::UserDataOperator::getSNESCtx() const {
949  return getFEMethod()->snes_ctx;
950 }
951 
953 ForcesAndSourcesCore::UserDataOperator::getTSCtx() const {
954  return getFEMethod()->ts_ctx;
955 }
956 
957 Vec ForcesAndSourcesCore::UserDataOperator::getKSPf() const {
958  return getFEMethod()->ksp_f;
959 }
960 
961 Mat ForcesAndSourcesCore::UserDataOperator::getKSPA() const {
962  return getFEMethod()->ksp_A;
963 }
964 
965 Mat ForcesAndSourcesCore::UserDataOperator::getKSPB() const {
966  return getFEMethod()->ksp_B;
967 }
968 
969 Vec ForcesAndSourcesCore::UserDataOperator::getSNESf() const {
970  return getFEMethod()->snes_f;
971 }
972 
973 Vec ForcesAndSourcesCore::UserDataOperator::getSNESx() const {
974  return getFEMethod()->snes_x;
975 }
976 
977 Mat ForcesAndSourcesCore::UserDataOperator::getSNESA() const {
978  return getFEMethod()->snes_A;
979 }
980 
981 Mat ForcesAndSourcesCore::UserDataOperator::getSNESB() const {
982  return getFEMethod()->snes_B;
983 }
984 
985 Vec ForcesAndSourcesCore::UserDataOperator::getTSu() const {
986  return getFEMethod()->ts_u;
987 }
988 
989 Vec ForcesAndSourcesCore::UserDataOperator::getTSu_t() const {
990  return getFEMethod()->ts_u_t;
991 }
992 
993 Vec ForcesAndSourcesCore::UserDataOperator::getTSu_tt() const {
994  return getFEMethod()->ts_u_tt;
995 }
996 
997 Vec ForcesAndSourcesCore::UserDataOperator::getTSf() const {
998  return getFEMethod()->ts_F;
999 }
1000 
1001 Mat ForcesAndSourcesCore::UserDataOperator::getTSA() const {
1002  return getFEMethod()->ts_A;
1003 }
1004 
1005 Mat ForcesAndSourcesCore::UserDataOperator::getTSB() const {
1006  return getFEMethod()->ts_B;
1007 }
1008 
1009 int ForcesAndSourcesCore::UserDataOperator::getTSstep() const {
1010  return getFEMethod()->ts_step;
1011 }
1012 
1013 double ForcesAndSourcesCore::UserDataOperator::getTStime() const {
1014  return getFEMethod()->ts_t;
1015 }
1016 
1017 double ForcesAndSourcesCore::UserDataOperator::getTSa() const {
1018  return getFEMethod()->ts_a;
1019 }
1020 
1021 MatrixDouble &ForcesAndSourcesCore::UserDataOperator::getGaussPts() {
1022  return static_cast<ForcesAndSourcesCore *>(ptrFE)->gaussPts;
1023 }
1024 
1025 auto ForcesAndSourcesCore::UserDataOperator::getFTensor0IntegrationWeight() {
1027  &(getGaussPts()(getGaussPts().size1() - 1, 0)));
1028 }
1029 
1030 MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::getPorblemRowIndices(
1031  const std::string filed_name, const EntityType type, const int side,
1032  VectorInt &indices) const {
1033  return getProblemRowIndices(filed_name, type, side, indices);
1034 }
1035 
1036 ForcesAndSourcesCore *ForcesAndSourcesCore::UserDataOperator::getPtrFE() const {
1037  return ptrFE;
1038 }
1039 
1041 ForcesAndSourcesCore::UserDataOperator::getSidePtrFE() const {
1042  return ptrFE->sidePtrFE;
1043 }
1044 
1045 } // namespace MoFEM
1046 
1047 #endif //__FORCES_AND_SOURCES_CORE__HPP__
1048 
1049 /**
1050  * \defgroup mofem_forces_and_sources Forces and sources
1051  * \ingroup mofem
1052  *
1053  * \brief Manages complexities related to assembly of vector and matrices at
1054  * single finite element level.
1055  *
1056  **/
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:74
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:180
Edge finite elementUser is implementing own operator at Gauss points level, by own object derived fro...
RuleHookFun getRuleHook
Hook to get rule.
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.
FieldApproximationBase
approximation base
Definition: definitions.h:149
boost::shared_ptr< BaseFunction > elementPolynomialBasePtr
Pointer to entity polynomial base.
constexpr int order
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
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:173
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
#define CHKERR
Inline error check.
Definition: definitions.h:601
ForcesAndSourcesCore::UserDataOperator UserDataOperator
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:71
ublas::vector< boost::shared_ptr< const FEDofEntity >, DofsAllocator > VectorDofs
ContactPrism finite elementUser is implementing own operator at Gauss points level,...
#define DEPRECATED
Definition: definitions.h:29
DataForcesAndSourcesCore & dataHdiv
Data on single entity (This is passed as argument to DataOperator::doWork)
structure to get information form mofem into DataForcesAndSourcesCore
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
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