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  /**@}*/
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  template <typename EXTRACTOR>
559  getNodesIndices(const std::string &field_name,
560  FieldEntity_vector_view &ents_field, VectorInt &nodes_indices,
561  VectorInt &local_nodes_indices, EXTRACTOR &&extractor) const;
562 
563  /// \brief get row node indices from FENumeredDofEntity_multiIndex
564  MoFEMErrorCode getRowNodesIndices(DataForcesAndSourcesCore &data,
565  const std::string &field_name) const;
566 
567  /// \brief get col node indices from FENumeredDofEntity_multiIndex
568  MoFEMErrorCode getColNodesIndices(DataForcesAndSourcesCore &data,
569  const std::string &field_name) const;
570 
571  template <typename EXTRACTOR>
572  MoFEMErrorCode getEntityIndices(DataForcesAndSourcesCore &data,
573  const std::string &field_name,
574  FieldEntity_vector_view &ents_field,
575  const EntityType type_lo,
576  const EntityType type_hi,
577  EXTRACTOR &&extractor) const;
578 
580  getEntityRowIndices(DataForcesAndSourcesCore &data,
581  const std::string &field_name,
582  const EntityType type_lo = MBVERTEX,
583  const EntityType type_hi = MBPOLYHEDRON) const;
584 
586  getEntityColIndices(DataForcesAndSourcesCore &data,
587  const std::string &field_name,
588  const EntityType type_lo = MBVERTEX,
589  const EntityType type_hi = MBPOLYHEDRON) const;
590 
591  /// \brief get NoField indices
593  getNoFieldIndices(const std::string &field_name,
594  boost::shared_ptr<FENumeredDofEntity_multiIndex> dofs,
595  VectorInt &nodes_indices) const;
596 
597  /// \brief get col NoField indices
598  MoFEMErrorCode getNoFieldRowIndices(DataForcesAndSourcesCore &data,
599  const std::string &field_name) const;
600 
601  /// \brief get col NoField indices
602  MoFEMErrorCode getNoFieldColIndices(DataForcesAndSourcesCore &data,
603  const std::string &field_name) const;
604 
605  /**@}*/
606 
607  /** \name Data */
608 
609  /**@{*/
610 
611  /**
612  * \brief Get field data on nodes
613  */
614  MoFEMErrorCode getNoFieldFieldData(const std::string field_name,
615  VectorDouble &ent_field_data,
616  VectorDofs &ent_field_dofs,
617  VectorFieldEntities &ent_field) const;
618 
619  MoFEMErrorCode getNoFieldFieldData(DataForcesAndSourcesCore &data,
620  const std::string field_name) const;
621 
622  /**
623  * \brief Get data on nodes
624  * @param data Data structure
625  * @param field_name Field name
626  * @return Error code
627  */
628  MoFEMErrorCode getNodesFieldData(DataForcesAndSourcesCore &data,
629  const std::string &field_name) const;
630 
632  getEntityFieldData(DataForcesAndSourcesCore &data,
633  const std::string &field_name,
634  const EntityType type_lo = MBVERTEX,
635  const EntityType type_hi = MBPOLYHEDRON) const;
636 
637  /**@}*/
638 
639  /// \brief Get nodes on triangles
640  MoFEMErrorCode getFaceTriNodes(DataForcesAndSourcesCore &data) const;
641 
642  /// \brief Get field approximation space and base on entities
644  getSpacesAndBaseOnEntities(DataForcesAndSourcesCore &data) const;
645 
646  /** \name Data form NumeredDofEntity_multiIndex */
647 
648  /**@{*/
649 
650  /// \brief get indices of nodal indices which are declared for problem but
651  /// not this particular element
652  MoFEMErrorCode getProblemNodesIndices(const std::string &field_name,
653  const NumeredDofEntity_multiIndex &dofs,
654  VectorInt &nodes_indices) const;
655 
656  /// \brief get indices by type (generic function) which are declared for
657  /// problem but not this particular element
658  MoFEMErrorCode getProblemTypeIndices(const std::string &field_name,
659  const NumeredDofEntity_multiIndex &dofs,
660  EntityType type, int side_number,
661  VectorInt &indices) const;
662 
663  MoFEMErrorCode getProblemNodesRowIndices(const std::string &field_name,
664  VectorInt &nodes_indices) const;
665  MoFEMErrorCode getProblemTypeRowIndices(const std::string &field_name,
666  EntityType type, int side_number,
667  VectorInt &indices) const;
668  MoFEMErrorCode getProblemNodesColIndices(const std::string &field_name,
669  VectorInt &nodes_indices) const;
670  MoFEMErrorCode getProblemTypeColIndices(const std::string &field_name,
671  EntityType type, int side_number,
672  VectorInt &indices) const;
673 
674  /**@}*/
675 
676  /**
677  * \brief another variant of getRule
678  * @param order_row order of base function on row
679  * @param order_col order of base function on columns
680  * @param order_data order of base function approximating data
681  * @return integration rule
682  *
683  * This function is overloaded by the user. The integration rule
684  * is set such that specific operator implemented by the user is
685  * integrated accurately. For example if user implement bilinear operator
686  * \f[
687  * b(u,v) =
688  * \int_\mathcal{T}
689  * \frac{\partial u_i}{\partial x_j}\frac{\partial v_i}{\partial x_j}
690  * \textrm{d}\mathcal{T}
691  * \f]
692  * then if \f$u\f$ and \f$v\f$ are polynomial of given \em order, then
693  * exact integral would be \code int getRule(int order) { return
694  * 2*(order-1); }; \endcode
695  *
696  * The integration points and weights are set appropriately for given
697  * entity type and integration rule from \ref quad.c
698  *
699  * Method \ref ForcesAndSourcesCore::getRule takes at argument takes
700  * maximal polynomial order set on the element on all fields defined on
701  * the element. If a user likes to have more control, another variant of
702  * this function can be called which distinguishing between field orders
703  * on rows, columns and data, the i.e. first argument of a bilinear form,
704  * the second argument of bilinear form and field coefficients on the
705  * element.
706  *
707  * \note If user set rule to -1 or any other negative integer, then method
708  * \ref ForcesAndSourcesCore::setGaussPts is called. In that method user
709  * can implement own (specific) integration method.
710  *
711  * \bug this function should be const
712  */
713  virtual int getRule(int order_row, int order_col, int order_data);
714 
715  /** \brief set user specific integration rule
716 
717  This function allows for user defined integration rule. The key is to
718  called matrix gaussPts, which is used by other MoFEM procedures. Matrix
719  has number of rows equal to problem dimension plus one, where last index
720  is used to store weight values. %Number of columns is equal to number of
721  integration points.
722 
723  \note This function is called if method \ref
724  ForcesAndSourcesCore::getRule is returning integer -1 or any other
725  negative integer.
726 
727  User sets
728  \code
729  MatrixDouble gaussPts;
730  \endcode
731  where
732  \code
733  gaussPts.resize(dim+1,nb_gauss_pts);
734  \endcode
735  number rows represents local coordinates of integration points
736  in reference element, where last index in row is for integration weight.
737 
738  */
739  virtual MoFEMErrorCode setGaussPts(int order_row, int order_col,
740  int order_data);
741 
742  /**
743  * \brief Calculate base functions
744  * @return Error code
745  */
747  calHierarchicalBaseFunctionsOnElement(const FieldApproximationBase b);
748 
749  /**
750  * \brief Calculate base functions
751  * @return Error code
752  */
753  MoFEMErrorCode calHierarchicalBaseFunctionsOnElement();
754 
755  /**
756  * @brief Calculate Bernstein-Bezier base
757  *
758  * @return MoFEMErrorCode
759  */
760  MoFEMErrorCode calBernsteinBezierBaseFunctionsOnElement();
761 
762  /**
763  * @brief Create a entity data on element object
764  *
765  * @return MoFEMErrorCode
766  */
767  MoFEMErrorCode createDataOnElement();
768 
769  /**
770  * @brief Iterate user data operators
771  *
772  * @return MoFEMErrorCode
773  */
774  MoFEMErrorCode loopOverOperators();
775 
776  /**@{*/
777 
778  /** \name Deprecated (do not use) */
779 
780  /** \deprecated Use getRule(int row_order, int col_order, int data order)
781  */
782  virtual int getRule(int order);
783 
784  /** \deprecated setGaussPts(int row_order, int col_order, int data order);
785  */
786  virtual MoFEMErrorCode setGaussPts(int order);
787 
788  /**@/}*/
789 
790  /**
791  * @brief Entity data on element entity rows fields
792  *
793  */
794  const std::array<boost::shared_ptr<DataForcesAndSourcesCore>, LASTSPACE>
796 
797  /**
798  * @brief Entity data on element entity columns fields
799  *
800  */
801  const std::array<boost::shared_ptr<DataForcesAndSourcesCore>, LASTSPACE>
803 
809 
810  /**
811  * @brief Vector of finite element users data operators
812  *
813  */
814  boost::ptr_vector<UserDataOperator> opPtrVector;
815 
816  friend class UserDataOperator;
817 
818 protected:
819  /**
820  * @brief Last evaluated type of element entity
821  *
822  */
824 
825 private:
826  /**
827  * @brief Pointer to entity polynomial base
828  *
829  */
830  boost::shared_ptr<BaseFunction> elementPolynomialBasePtr;
831 
832  /**
833  * @brief Pointer to user polynomail base
834  */
835  boost::shared_ptr<BaseFunction> userPolynomialBasePtr;
836 
837  /**
838  * @brief Element to integrate on the sides
839  *
840  */
842 
843  /**
844  * @brief Set the pointer to face element on the side
845  *
846  * \note Function is is used by face element, while it iterates over
847  * elements on the side
848  *
849  * @param side_fe_ptr
850  * @return MoFEMErrorCode
851  */
852  MoFEMErrorCode setSideFEPtr(const ForcesAndSourcesCore *side_fe_ptr);
853 
857 };
858 
859 /// \deprecated Used ForcesAndSourcesCore instead
861 
862 
863 boost::shared_ptr<const NumeredEntFiniteElement>
864 ForcesAndSourcesCore::UserDataOperator::getNumeredEntFiniteElementPtr() const {
865  return ptrFE->numeredEntFiniteElementPtr;
866 };
867 
868 EntityHandle ForcesAndSourcesCore::UserDataOperator::getFEEntityHandle() const {
869  return getNumeredEntFiniteElementPtr()->getEnt();
870 }
871 
872 boost::weak_ptr<SideNumber>
873 ForcesAndSourcesCore::UserDataOperator::getSideNumberPtr(
874  const int side_number, const EntityType type) {
875  auto &side_table_by_side_and_type =
876  ptrFE->numeredEntFiniteElementPtr->getSideNumberTable().get<1>();
877  auto side_it =
878  side_table_by_side_and_type.find(boost::make_tuple(type, side_number));
879  if (side_it != side_table_by_side_and_type.end())
880  return *side_it;
881  else
882  return boost::weak_ptr<SideNumber>();
883 }
884 
886 ForcesAndSourcesCore::UserDataOperator::getSideEntity(const int side_number,
887  const EntityType type) {
888  if (auto side_ptr = getSideNumberPtr(side_number, type).lock())
889  return side_ptr->ent;
890  else
891  return 0;
892 }
893 
894 int ForcesAndSourcesCore::UserDataOperator::getNumberOfNodesOnElement() {
895  int num_nodes;
896  CHKERR ptrFE->getNumberOfNodes(num_nodes);
897  return num_nodes;
898 }
899 
900 const FEMethod *ForcesAndSourcesCore::UserDataOperator::getFEMethod() const {
901  return ptrFE;
902 }
903 
904 int ForcesAndSourcesCore::UserDataOperator::getOpType() const { return opType; }
905 
906 void ForcesAndSourcesCore::UserDataOperator::setOpType(const OpType type) {
907  opType = type;
908 }
909 
910 void ForcesAndSourcesCore::UserDataOperator::addOpType(const OpType type) {
911  opType |= type;
912 }
913 
914 int ForcesAndSourcesCore::UserDataOperator::getNinTheLoop() const {
915  return getFEMethod()->getNinTheLoop();
916 }
917 
918 int ForcesAndSourcesCore::UserDataOperator::getLoopSize() const {
919  return getFEMethod()->getLoopSize();
920 }
921 
922 const std::string &ForcesAndSourcesCore::UserDataOperator::getFEName() const {
923  return getFEMethod()->feName;
924 }
925 
926 const PetscData::Switches &
927 ForcesAndSourcesCore::UserDataOperator::getDataCtx() const {
928  return getFEMethod()->data_ctx;
929 }
930 
932 ForcesAndSourcesCore::UserDataOperator::getKSPCtx() const {
933  return getFEMethod()->ksp_ctx;
934 }
935 
937 ForcesAndSourcesCore::UserDataOperator::getSNESCtx() const {
938  return getFEMethod()->snes_ctx;
939 }
940 
942 ForcesAndSourcesCore::UserDataOperator::getTSCtx() const {
943  return getFEMethod()->ts_ctx;
944 }
945 
946 Vec ForcesAndSourcesCore::UserDataOperator::getKSPf() const {
947  return getFEMethod()->ksp_f;
948 }
949 
950 Mat ForcesAndSourcesCore::UserDataOperator::getKSPA() const {
951  return getFEMethod()->ksp_A;
952 }
953 
954 Mat ForcesAndSourcesCore::UserDataOperator::getKSPB() const {
955  return getFEMethod()->ksp_B;
956 }
957 
958 Vec ForcesAndSourcesCore::UserDataOperator::getSNESf() const {
959  return getFEMethod()->snes_f;
960 }
961 
962 Vec ForcesAndSourcesCore::UserDataOperator::getSNESx() const {
963  return getFEMethod()->snes_x;
964 }
965 
966 Mat ForcesAndSourcesCore::UserDataOperator::getSNESA() const {
967  return getFEMethod()->snes_A;
968 }
969 
970 Mat ForcesAndSourcesCore::UserDataOperator::getSNESB() const {
971  return getFEMethod()->snes_B;
972 }
973 
974 Vec ForcesAndSourcesCore::UserDataOperator::getTSu() const {
975  return getFEMethod()->ts_u;
976 }
977 
978 Vec ForcesAndSourcesCore::UserDataOperator::getTSu_t() const {
979  return getFEMethod()->ts_u_t;
980 }
981 
982 Vec ForcesAndSourcesCore::UserDataOperator::getTSu_tt() const {
983  return getFEMethod()->ts_u_tt;
984 }
985 
986 Vec ForcesAndSourcesCore::UserDataOperator::getTSf() const {
987  return getFEMethod()->ts_F;
988 }
989 
990 Mat ForcesAndSourcesCore::UserDataOperator::getTSA() const {
991  return getFEMethod()->ts_A;
992 }
993 
994 Mat ForcesAndSourcesCore::UserDataOperator::getTSB() const {
995  return getFEMethod()->ts_B;
996 }
997 
998 int ForcesAndSourcesCore::UserDataOperator::getTSstep() const {
999  return getFEMethod()->ts_step;
1000 }
1001 
1002 double ForcesAndSourcesCore::UserDataOperator::getTStime() const {
1003  return getFEMethod()->ts_t;
1004 }
1005 
1006 double ForcesAndSourcesCore::UserDataOperator::getTSa() const {
1007  return getFEMethod()->ts_a;
1008 }
1009 
1010 MatrixDouble &ForcesAndSourcesCore::UserDataOperator::getGaussPts() {
1011  return static_cast<ForcesAndSourcesCore *>(ptrFE)->gaussPts;
1012 }
1013 
1014 auto ForcesAndSourcesCore::UserDataOperator::getFTensor0IntegrationWeight() {
1016  &(getGaussPts()(getGaussPts().size1() - 1, 0)));
1017 }
1018 
1019 MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::getPorblemRowIndices(
1020  const std::string filed_name, const EntityType type, const int side,
1021  VectorInt &indices) const {
1022  return getProblemRowIndices(filed_name, type, side, indices);
1023 }
1024 
1025 ForcesAndSourcesCore *ForcesAndSourcesCore::UserDataOperator::getPtrFE() const {
1026  return ptrFE;
1027 }
1028 
1030 ForcesAndSourcesCore::UserDataOperator::getSidePtrFE() const {
1031  return ptrFE->sidePtrFE;
1032 }
1033 
1034 } // namespace MoFEM
1035 
1036 #endif //__FORCES_AND_SOURCES_CORE__HPP__
1037 
1038 /**
1039  * \defgroup mofem_forces_and_sources Forces and sources
1040  * \ingroup mofem
1041  *
1042  * \brief Manages complexities related to assembly of vector and matrices at
1043  * single finite element level.
1044  *
1045  **/
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
std::vector< boost::weak_ptr< FieldEntity > > FieldEntity_vector_view
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
auto & getUserPolynomialBase()
Get the User Polynomial Base object.
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...
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.
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
RuleHookFun getRuleHook
Hook to get rule.
const int dim
ublas::vector< FieldEntity *, FieldEntAllocator > VectorFieldEntities
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
constexpr int order
FieldApproximationBase
approximation base
Definition: definitions.h:150
boost::shared_ptr< BaseFunction > elementPolynomialBasePtr
Pointer to entity polynomial base.
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
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:604
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:73
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
Base volume element used to integrate on contact surface (could be extended to other volume elements)...
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
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.
Face finite elementUser is implementing own operator at Gauss point level, by own object derived from...
DEPRECATED typedef ForcesAndSourcesCore ForcesAndSurcesCore
DataForcesAndSourcesCore & dataNoField