v0.10.0
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  inline double getTSaa() const;
340 
341  /**@}*/
342 
343  /**@{*/
344 
345  /** \name Base funtions and integration points */
346 
347  /** \brief matrix of integration (Gauss) points for Volume Element
348  *
349  * For triangle: columns 0,1 are x,y coordinates respectively and column
350  * 2 is a weight value for example getGaussPts()(1,13) returns y
351  * coordinate of 13th Gauss point on particular volume element
352  *
353  * For tetrahedron: columns 0,1,2 are x,y,z coordinates respectively and
354  * column 3 is a weight value for example getGaussPts()(1,13) returns y
355  * coordinate of 13th Gauss point on particular volume element
356  *
357  */
358  inline MatrixDouble &getGaussPts();
359 
360  /**
361  * @brief Get integration weights
362  *
363  * \code
364  * auto t_w = getFTensor0IntegrationWeight();
365  * for(int gg = 0; gg!=getGaussPts.size2(); ++gg) {
366  * // integrate something
367  * ++t_w;
368  * }
369  * \endcode
370  *
371  * @return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>
372  */
373  inline auto getFTensor0IntegrationWeight();
374 
375  /**@}*/
376 
377  /**@{*/
378 
379  /** \name Deprecated (do not use) */
380 
381  // \deprecated Deprecated function with spelling mistake
383  getPorblemRowIndices(const std::string filed_name, const EntityType type,
384  const int side, VectorInt &indices) const;
385 
386  /**@}*/
387 
388  protected:
390 
391  virtual MoFEMErrorCode setPtrFE(ForcesAndSourcesCore *ptr);
392 
393  inline ForcesAndSourcesCore *getPtrFE() const;
394 
395  inline ForcesAndSourcesCore *getSidePtrFE() const;
396 
397  private:
398  /**
399  * @brief User call this function to loop over elements on the side of
400  * face. This function calls finite element with is operator to do
401  * calculations.
402  *
403  * @param fe_name name of the side element
404  * @param side_fe pointer to the side element instance
405  * @param dim dimension the of side element
406  * @param ent_for_side entity handle for which adjacent volume or face will
407  * be accessed
408  * @return MoFEMErrorCode
409  */
410  MoFEMErrorCode loopSide(const string &fe_name,
411  ForcesAndSourcesCore *side_fe, const size_t dim,
412  const EntityHandle ent_for_side = 0);
413 
414  friend class ForcesAndSourcesCore;
418  };
419 
420  /** \brief Use to push back operator for row operator
421 
422  It can be used to calculate nodal forces or other quantities on the mesh.
423 
424  */
425  boost::ptr_vector<UserDataOperator> &getOpPtrVector() { return opPtrVector; }
426 
427  /**
428  * @brief Get the Entity Polynomial Base object
429  *
430  * @return boost::shared_ptr<BaseFunction>&&
431  */
432  auto &getElementPolynomialBase() { return elementPolynomialBasePtr; }
433 
434  /**
435  * @brief Get the User Polynomial Base object
436  *
437  * @return boost::shared_ptr<BaseFunction>&
438  */
439  auto &getUserPolynomialBase() { return userPolynomialBasePtr; }
440 
441  /**
442  * @brief Matrix of integration points
443  *
444  * Columns is equal to number of integration points, numver of rows
445  * depends on dimension of finite element entity, for example for
446  * tetrahedron rows are x,y,z,weight. Last row is integration weight.
447  *
448  * FIXME: that should be moved to private class data and acessed only by
449  * member function
450  */
452 
453  virtual MoFEMErrorCode preProcess();
454  virtual MoFEMErrorCode operator()();
455  virtual MoFEMErrorCode postProcess();
456 
457 public:
458  /** \brief Get max order of approximation for data fields
459 
460  Method getMaxDataOrder () return maximal order on entities, for
461  all data on the element. So for example if finite element is triangle, and
462  triangle base function have order 4 and on edges base function have order
463  2, this function return 4.
464 
465  If finite element has for example 2 or more approximated fields, for
466  example Pressure (order 3) and displacement field (order 5), this function
467  returns 5.
468 
469  */
470  int getMaxDataOrder() const;
471 
472  /// \brief Get max order of approximation for field in rows
473  int getMaxRowOrder() const;
474 
475  /// \brief Get max order of approximation for field in columns
476  int getMaxColOrder() const;
477 
478  /**
479  * @brief Get the entity data
480  *
481  * @param space
482  * @param type
483  * @param side
484  * @return const DataForcesAndSourcesCore::EntData&
485  */
487  const EntityType type,
488  const int side) const {
489  return dataOnElement[space]->dataOnEntities[type][side];
490  }
491 
492  /**
493  * @brief Get the entity data
494  *
495  * @param space
496  * @param type
497  * @param side
498  * @return DataForcesAndSourcesCore::EntData&
499  */
501  getEntData(const FieldSpace space, const EntityType type, const int side) {
502  return dataOnElement[space]->dataOnEntities[type][side];
503  }
504 
505 protected:
506  /**
507  * \brief get sense (orientation) of entity
508  * @param type type of entity
509  * @param data entity data
510  * @return error code
511  */
512  MoFEMErrorCode getEntitySense(
513  const EntityType type,
514  boost::ptr_vector<DataForcesAndSourcesCore::EntData> &data) const;
515 
516  /**
517  * @brief Get the entity data order
518  *
519  * @param type
520  * @param space
521  * @param data
522  * @return MoFEMErrorCode
523  */
524  MoFEMErrorCode getEntityDataOrder(
525  const EntityType type, const FieldSpace space,
526  boost::ptr_vector<DataForcesAndSourcesCore::EntData> &data) const;
527 
528  /**
529  * @brief Get the entity sense (orientation)
530  *
531  * @tparam type
532  * @param data
533  * @return MoFEMErrorCode
534  */
535  template <EntityType type>
537  return getEntitySense(type, data.dataOnEntities[type]);
538  }
539 
540  /**
541  * @brief Get the entity data order for given space
542  *
543  * @tparam type
544  * @param data
545  * @param space
546  * @return MoFEMErrorCode
547  */
548  template <EntityType type>
550  const FieldSpace space) const {
551  return getEntityDataOrder(type, space, data.dataOnEntities[type]);
552  }
553 
554  /** \name Indices */
555 
556  /**@{*/
557 
558  /// \brief get node indices
559  template <typename EXTRACTOR>
561  getNodesIndices(const std::string &field_name,
562  FieldEntity_vector_view &ents_field, VectorInt &nodes_indices,
563  VectorInt &local_nodes_indices, EXTRACTOR &&extractor) const;
564 
565  /// \brief get row node indices from FENumeredDofEntity_multiIndex
566  MoFEMErrorCode getRowNodesIndices(DataForcesAndSourcesCore &data,
567  const std::string &field_name) const;
568 
569  /// \brief get col node indices from FENumeredDofEntity_multiIndex
570  MoFEMErrorCode getColNodesIndices(DataForcesAndSourcesCore &data,
571  const std::string &field_name) const;
572 
573  template <typename EXTRACTOR>
574  MoFEMErrorCode getEntityIndices(DataForcesAndSourcesCore &data,
575  const std::string &field_name,
576  FieldEntity_vector_view &ents_field,
577  const EntityType type_lo,
578  const EntityType type_hi,
579  EXTRACTOR &&extractor) const;
580 
582  getEntityRowIndices(DataForcesAndSourcesCore &data,
583  const std::string &field_name,
584  const EntityType type_lo = MBVERTEX,
585  const EntityType type_hi = MBPOLYHEDRON) const;
586 
588  getEntityColIndices(DataForcesAndSourcesCore &data,
589  const std::string &field_name,
590  const EntityType type_lo = MBVERTEX,
591  const EntityType type_hi = MBPOLYHEDRON) const;
592 
593  /// \brief get NoField indices
595  getNoFieldIndices(const std::string &field_name,
596  boost::shared_ptr<FENumeredDofEntity_multiIndex> dofs,
597  VectorInt &nodes_indices) const;
598 
599  /// \brief get col NoField indices
600  MoFEMErrorCode getNoFieldRowIndices(DataForcesAndSourcesCore &data,
601  const std::string &field_name) const;
602 
603  /// \brief get col NoField indices
604  MoFEMErrorCode getNoFieldColIndices(DataForcesAndSourcesCore &data,
605  const std::string &field_name) const;
606 
607  /**@}*/
608 
609  /** \name Data */
610 
611  /**@{*/
612 
613  /**
614  * \brief Get field data on nodes
615  */
616  MoFEMErrorCode getNoFieldFieldData(const std::string field_name,
617  VectorDouble &ent_field_data,
618  VectorDofs &ent_field_dofs,
619  VectorFieldEntities &ent_field) const;
620 
621  MoFEMErrorCode getNoFieldFieldData(DataForcesAndSourcesCore &data,
622  const std::string field_name) const;
623 
624  /**
625  * \brief Get data on nodes
626  * @param data Data structure
627  * @param field_name Field name
628  * @return Error code
629  */
630  MoFEMErrorCode getNodesFieldData(DataForcesAndSourcesCore &data,
631  const std::string &field_name) const;
632 
634  getEntityFieldData(DataForcesAndSourcesCore &data,
635  const std::string &field_name,
636  const EntityType type_lo = MBVERTEX,
637  const EntityType type_hi = MBPOLYHEDRON) const;
638 
639  /**@}*/
640 
641  /// \brief Get nodes on triangles
642  MoFEMErrorCode getFaceTriNodes(DataForcesAndSourcesCore &data) const;
643 
644  /// \brief Get field approximation space and base on entities
646  getSpacesAndBaseOnEntities(DataForcesAndSourcesCore &data) const;
647 
648  /** \name Data form NumeredDofEntity_multiIndex */
649 
650  /**@{*/
651 
652  /// \brief get indices of nodal indices which are declared for problem but
653  /// not this particular element
654  MoFEMErrorCode getProblemNodesIndices(const std::string &field_name,
655  const NumeredDofEntity_multiIndex &dofs,
656  VectorInt &nodes_indices) const;
657 
658  /// \brief get indices by type (generic function) which are declared for
659  /// problem but not this particular element
660  MoFEMErrorCode getProblemTypeIndices(const std::string &field_name,
661  const NumeredDofEntity_multiIndex &dofs,
662  EntityType type, int side_number,
663  VectorInt &indices) const;
664 
665  MoFEMErrorCode getProblemNodesRowIndices(const std::string &field_name,
666  VectorInt &nodes_indices) const;
667  MoFEMErrorCode getProblemTypeRowIndices(const std::string &field_name,
668  EntityType type, int side_number,
669  VectorInt &indices) const;
670  MoFEMErrorCode getProblemNodesColIndices(const std::string &field_name,
671  VectorInt &nodes_indices) const;
672  MoFEMErrorCode getProblemTypeColIndices(const std::string &field_name,
673  EntityType type, int side_number,
674  VectorInt &indices) const;
675 
676  /**@}*/
677 
678  /**
679  * \brief another variant of getRule
680  * @param order_row order of base function on row
681  * @param order_col order of base function on columns
682  * @param order_data order of base function approximating data
683  * @return integration rule
684  *
685  * This function is overloaded by the user. The integration rule
686  * is set such that specific operator implemented by the user is
687  * integrated accurately. For example if user implement bilinear operator
688  * \f[
689  * b(u,v) =
690  * \int_\mathcal{T}
691  * \frac{\partial u_i}{\partial x_j}\frac{\partial v_i}{\partial x_j}
692  * \textrm{d}\mathcal{T}
693  * \f]
694  * then if \f$u\f$ and \f$v\f$ are polynomial of given \em order, then
695  * exact integral would be \code int getRule(int order) { return
696  * 2*(order-1); }; \endcode
697  *
698  * The integration points and weights are set appropriately for given
699  * entity type and integration rule from \ref quad.c
700  *
701  * Method \ref ForcesAndSourcesCore::getRule takes at argument takes
702  * maximal polynomial order set on the element on all fields defined on
703  * the element. If a user likes to have more control, another variant of
704  * this function can be called which distinguishing between field orders
705  * on rows, columns and data, the i.e. first argument of a bilinear form,
706  * the second argument of bilinear form and field coefficients on the
707  * element.
708  *
709  * \note If user set rule to -1 or any other negative integer, then method
710  * \ref ForcesAndSourcesCore::setGaussPts is called. In that method user
711  * can implement own (specific) integration method.
712  *
713  * \bug this function should be const
714  */
715  virtual int getRule(int order_row, int order_col, int order_data);
716 
717  /** \brief set user specific integration rule
718 
719  This function allows for user defined integration rule. The key is to
720  called matrix gaussPts, which is used by other MoFEM procedures. Matrix
721  has number of rows equal to problem dimension plus one, where last index
722  is used to store weight values. %Number of columns is equal to number of
723  integration points.
724 
725  \note This function is called if method \ref
726  ForcesAndSourcesCore::getRule is returning integer -1 or any other
727  negative integer.
728 
729  User sets
730  \code
731  MatrixDouble gaussPts;
732  \endcode
733  where
734  \code
735  gaussPts.resize(dim+1,nb_gauss_pts);
736  \endcode
737  number rows represents local coordinates of integration points
738  in reference element, where last index in row is for integration weight.
739 
740  */
741  virtual MoFEMErrorCode setGaussPts(int order_row, int order_col,
742  int order_data);
743 
744  /**
745  * \brief Calculate base functions
746  * @return Error code
747  */
749  calHierarchicalBaseFunctionsOnElement(const FieldApproximationBase b);
750 
751  /**
752  * \brief Calculate base functions
753  * @return Error code
754  */
755  MoFEMErrorCode calHierarchicalBaseFunctionsOnElement();
756 
757  /**
758  * @brief Calculate Bernstein-Bezier base
759  *
760  * @return MoFEMErrorCode
761  */
762  MoFEMErrorCode calBernsteinBezierBaseFunctionsOnElement();
763 
764  /**
765  * @brief Create a entity data on element object
766  *
767  * @return MoFEMErrorCode
768  */
769  MoFEMErrorCode createDataOnElement();
770 
771  /**
772  * @brief Iterate user data operators
773  *
774  * @return MoFEMErrorCode
775  */
776  MoFEMErrorCode loopOverOperators();
777 
778  /**@{*/
779 
780  /** \name Deprecated (do not use) */
781 
782  /** \deprecated Use getRule(int row_order, int col_order, int data order)
783  */
784  virtual int getRule(int order);
785 
786  /** \deprecated setGaussPts(int row_order, int col_order, int data order);
787  */
788  virtual MoFEMErrorCode setGaussPts(int order);
789 
790  /**@/}*/
791 
792  /**
793  * @brief Entity data on element entity rows fields
794  *
795  */
796  const std::array<boost::shared_ptr<DataForcesAndSourcesCore>, LASTSPACE>
798 
799  /**
800  * @brief Entity data on element entity columns fields
801  *
802  */
803  const std::array<boost::shared_ptr<DataForcesAndSourcesCore>, LASTSPACE>
805 
811 
812  /**
813  * @brief Vector of finite element users data operators
814  *
815  */
816  boost::ptr_vector<UserDataOperator> opPtrVector;
817 
818  friend class UserDataOperator;
819 
820 protected:
821  /**
822  * @brief Last evaluated type of element entity
823  *
824  */
826 
827 private:
828  /**
829  * @brief Pointer to entity polynomial base
830  *
831  */
832  boost::shared_ptr<BaseFunction> elementPolynomialBasePtr;
833 
834  /**
835  * @brief Pointer to user polynomail base
836  */
837  boost::shared_ptr<BaseFunction> userPolynomialBasePtr;
838 
839  /**
840  * @brief Element to integrate on the sides
841  *
842  */
844 
845  /**
846  * @brief Set the pointer to face element on the side
847  *
848  * \note Function is is used by face element, while it iterates over
849  * elements on the side
850  *
851  * @param side_fe_ptr
852  * @return MoFEMErrorCode
853  */
854  MoFEMErrorCode setSideFEPtr(const ForcesAndSourcesCore *side_fe_ptr);
855 
859 };
860 
861 /// \deprecated Used ForcesAndSourcesCore instead
863 
864 
865 boost::shared_ptr<const NumeredEntFiniteElement>
866 ForcesAndSourcesCore::UserDataOperator::getNumeredEntFiniteElementPtr() const {
867  return ptrFE->numeredEntFiniteElementPtr;
868 };
869 
870 EntityHandle ForcesAndSourcesCore::UserDataOperator::getFEEntityHandle() const {
871  return getNumeredEntFiniteElementPtr()->getEnt();
872 }
873 
874 boost::weak_ptr<SideNumber>
875 ForcesAndSourcesCore::UserDataOperator::getSideNumberPtr(
876  const int side_number, const EntityType type) {
877  auto &side_table_by_side_and_type =
878  ptrFE->numeredEntFiniteElementPtr->getSideNumberTable().get<1>();
879  auto side_it =
880  side_table_by_side_and_type.find(boost::make_tuple(type, side_number));
881  if (side_it != side_table_by_side_and_type.end())
882  return *side_it;
883  else
884  return boost::weak_ptr<SideNumber>();
885 }
886 
887 EntityHandle
888 ForcesAndSourcesCore::UserDataOperator::getSideEntity(const int side_number,
889  const EntityType type) {
890  if (auto side_ptr = getSideNumberPtr(side_number, type).lock())
891  return side_ptr->ent;
892  else
893  return 0;
894 }
895 
896 int ForcesAndSourcesCore::UserDataOperator::getNumberOfNodesOnElement() {
897  int num_nodes;
898  CHKERR ptrFE->getNumberOfNodes(num_nodes);
899  return num_nodes;
900 }
901 
902 const FEMethod *ForcesAndSourcesCore::UserDataOperator::getFEMethod() const {
903  return ptrFE;
904 }
905 
906 int ForcesAndSourcesCore::UserDataOperator::getOpType() const { return opType; }
907 
908 void ForcesAndSourcesCore::UserDataOperator::setOpType(const OpType type) {
909  opType = type;
910 }
911 
912 void ForcesAndSourcesCore::UserDataOperator::addOpType(const OpType type) {
913  opType |= type;
914 }
915 
916 int ForcesAndSourcesCore::UserDataOperator::getNinTheLoop() const {
917  return getFEMethod()->getNinTheLoop();
918 }
919 
920 int ForcesAndSourcesCore::UserDataOperator::getLoopSize() const {
921  return getFEMethod()->getLoopSize();
922 }
923 
924 const std::string &ForcesAndSourcesCore::UserDataOperator::getFEName() const {
925  return getFEMethod()->feName;
926 }
927 
928 const PetscData::Switches &
929 ForcesAndSourcesCore::UserDataOperator::getDataCtx() const {
930  return getFEMethod()->data_ctx;
931 }
932 
934 ForcesAndSourcesCore::UserDataOperator::getKSPCtx() const {
935  return getFEMethod()->ksp_ctx;
936 }
937 
939 ForcesAndSourcesCore::UserDataOperator::getSNESCtx() const {
940  return getFEMethod()->snes_ctx;
941 }
942 
944 ForcesAndSourcesCore::UserDataOperator::getTSCtx() const {
945  return getFEMethod()->ts_ctx;
946 }
947 
948 Vec ForcesAndSourcesCore::UserDataOperator::getKSPf() const {
949  return getFEMethod()->ksp_f;
950 }
951 
952 Mat ForcesAndSourcesCore::UserDataOperator::getKSPA() const {
953  return getFEMethod()->ksp_A;
954 }
955 
956 Mat ForcesAndSourcesCore::UserDataOperator::getKSPB() const {
957  return getFEMethod()->ksp_B;
958 }
959 
960 Vec ForcesAndSourcesCore::UserDataOperator::getSNESf() const {
961  return getFEMethod()->snes_f;
962 }
963 
964 Vec ForcesAndSourcesCore::UserDataOperator::getSNESx() const {
965  return getFEMethod()->snes_x;
966 }
967 
968 Mat ForcesAndSourcesCore::UserDataOperator::getSNESA() const {
969  return getFEMethod()->snes_A;
970 }
971 
972 Mat ForcesAndSourcesCore::UserDataOperator::getSNESB() const {
973  return getFEMethod()->snes_B;
974 }
975 
976 Vec ForcesAndSourcesCore::UserDataOperator::getTSu() const {
977  return getFEMethod()->ts_u;
978 }
979 
980 Vec ForcesAndSourcesCore::UserDataOperator::getTSu_t() const {
981  return getFEMethod()->ts_u_t;
982 }
983 
984 Vec ForcesAndSourcesCore::UserDataOperator::getTSu_tt() const {
985  return getFEMethod()->ts_u_tt;
986 }
987 
988 Vec ForcesAndSourcesCore::UserDataOperator::getTSf() const {
989  return getFEMethod()->ts_F;
990 }
991 
992 Mat ForcesAndSourcesCore::UserDataOperator::getTSA() const {
993  return getFEMethod()->ts_A;
994 }
995 
996 Mat ForcesAndSourcesCore::UserDataOperator::getTSB() const {
997  return getFEMethod()->ts_B;
998 }
999 
1000 int ForcesAndSourcesCore::UserDataOperator::getTSstep() const {
1001  return getFEMethod()->ts_step;
1002 }
1003 
1004 double ForcesAndSourcesCore::UserDataOperator::getTStime() const {
1005  return getFEMethod()->ts_t;
1006 }
1007 
1008 double ForcesAndSourcesCore::UserDataOperator::getTSa() const {
1009  return getFEMethod()->ts_a;
1010 }
1011 
1012 double ForcesAndSourcesCore::UserDataOperator::getTSaa() const {
1013  return getFEMethod()->ts_aa;
1014 }
1015 
1016 MatrixDouble &ForcesAndSourcesCore::UserDataOperator::getGaussPts() {
1017  return static_cast<ForcesAndSourcesCore *>(ptrFE)->gaussPts;
1018 }
1019 
1020 auto ForcesAndSourcesCore::UserDataOperator::getFTensor0IntegrationWeight() {
1022  &(getGaussPts()(getGaussPts().size1() - 1, 0)));
1023 }
1024 
1025 MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::getPorblemRowIndices(
1026  const std::string filed_name, const EntityType type, const int side,
1027  VectorInt &indices) const {
1028  return getProblemRowIndices(filed_name, type, side, indices);
1029 }
1030 
1031 ForcesAndSourcesCore *ForcesAndSourcesCore::UserDataOperator::getPtrFE() const {
1032  return ptrFE;
1033 }
1034 
1036 ForcesAndSourcesCore::UserDataOperator::getSidePtrFE() const {
1037  return ptrFE->sidePtrFE;
1038 }
1039 
1040 } // namespace MoFEM
1041 
1042 #endif //__FORCES_AND_SOURCES_CORE__HPP__
1043 
1044 /**
1045  * \defgroup mofem_forces_and_sources Forces and sources
1046  * \ingroup mofem
1047  *
1048  * \brief Manages complexities related to assembly of vector and matrices at
1049  * single finite element level.
1050  *
1051  **/
ForcesAndSourcesCore::UserDataOperator UserDataOperator
FieldApproximationBase
approximation base
Definition: definitions.h:150
FieldSpace
approximation spaces
Definition: definitions.h:174
@ LASTSPACE
FieldSpace in [ 0, LASTSPACE )
Definition: definitions.h:181
#define DEPRECATED
Definition: definitions.h:29
#define CHKERR
Inline error check.
Definition: definitions.h:604
const int dim
multi_index_container< boost::shared_ptr< NumeredDofEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, UId, &NumeredDofEntity::getLocalUniqueId > >, 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< 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 > > > > NumeredDofEntity_multiIndex
MultiIndex container keeps NumeredDofEntity.
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:73
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
ublas::vector< FieldEntity *, FieldEntAllocator > VectorFieldEntities
DEPRECATED typedef ForcesAndSourcesCore ForcesAndSurcesCore
std::vector< boost::weak_ptr< FieldEntity > > FieldEntity_vector_view
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
Data on single entity (This is passed as argument to DataOperator::doWork)
data structure for finite element entity
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
base operator to do operations at Gauss Pt. level
Deprecated interface functions.
structure for User Loop Methods on finite elements
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
Data operator to do calculations at integration points.
OpType
Controls loop over entities on element.
structure to get information form mofem into DataForcesAndSourcesCore
auto & getUserPolynomialBase()
Get the User Polynomial Base object.
DataForcesAndSourcesCore & dataNoField
const std::array< boost::shared_ptr< DataForcesAndSourcesCore >, LASTSPACE > derivedDataOnElement
Entity data on element entity columns fields.
const std::array< boost::shared_ptr< DataForcesAndSourcesCore >, LASTSPACE > dataOnElement
Entity data on element entity rows fields.
auto & getElementPolynomialBase()
Get the Entity Polynomial Base object.
MoFEMErrorCode getEntitySense(DataForcesAndSourcesCore &data) const
Get the entity sense (orientation)
ForcesAndSourcesCore * sidePtrFE
Element to integrate on the sides.
const DataForcesAndSourcesCore::EntData & getEntData(const FieldSpace space, const EntityType type, const int side) const
Get the entity data.
DataForcesAndSourcesCore & dataH1
DataForcesAndSourcesCore & dataHdiv
DataForcesAndSourcesCore & dataHcurl
MatrixDouble gaussPts
Matrix of integration points.
DataForcesAndSourcesCore::EntData & getEntData(const FieldSpace space, const EntityType type, const int side)
Get the entity data.
MoFEMErrorCode getEntityDataOrder(DataForcesAndSourcesCore &data, const FieldSpace space) const
Get the entity data order for given space.
RuleHookFun setRuleHook
Set function to calculate integration rule.
EntityType lastEvaluatedElementEntityType
Last evaluated type of element entity.
boost::function< int(int order_row, int order_col, int order_data)> RuleHookFun
boost::ptr_vector< UserDataOperator > opPtrVector
Vector of finite element users data operators.
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
DataForcesAndSourcesCore & dataL2
boost::shared_ptr< BaseFunction > elementPolynomialBasePtr
Pointer to entity polynomial base.
boost::shared_ptr< BaseFunction > userPolynomialBasePtr
Pointer to user polynomail base.
RuleHookFun getRuleHook
Hook to get rule.
KSPContext
pass information about context of KSP/DM for with finite element is computed
Definition: LoopMethods.hpp:96
std::bitset< 8 > Switches
Definition: LoopMethods.hpp:59
Base volume element used to integrate on contact surface (could be extended to other volume elements)...