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  /**
474  * @brief Get the entity data
475  *
476  * @param space
477  * @param type
478  * @param side
479  * @return const DataForcesAndSourcesCore::EntData&
480  */
482  const EntityType type,
483  const int side) const {
484  return dataOnElement[space]->dataOnEntities[type][side];
485  }
486 
487  /**
488  * @brief Get the entity data
489  *
490  * @param space
491  * @param type
492  * @param side
493  * @return DataForcesAndSourcesCore::EntData&
494  */
496  getEntData(const FieldSpace space, const EntityType type, const int side) {
497  return dataOnElement[space]->dataOnEntities[type][side];
498  }
499 
500 protected:
501  /**
502  * \brief get sense (orientation) of entity
503  * @param type type of entity
504  * @param data entity data
505  * @return error code
506  */
507  MoFEMErrorCode getEntitySense(
508  const EntityType type,
509  boost::ptr_vector<DataForcesAndSourcesCore::EntData> &data) const;
510 
511  /**
512  * @brief Get the entity data order
513  *
514  * @param type
515  * @param space
516  * @param data
517  * @return MoFEMErrorCode
518  */
519  MoFEMErrorCode getEntityDataOrder(
520  const EntityType type, const FieldSpace space,
521  boost::ptr_vector<DataForcesAndSourcesCore::EntData> &data) const;
522 
523  /**
524  * @brief Get the entity sense (orientation)
525  *
526  * @tparam type
527  * @param data
528  * @return MoFEMErrorCode
529  */
530  template <EntityType type>
532  return getEntitySense(type, data.dataOnEntities[type]);
533  }
534 
535  /**
536  * @brief Get the entity data order for given space
537  *
538  * @tparam type
539  * @param data
540  * @param space
541  * @return MoFEMErrorCode
542  */
543  template <EntityType type>
545  const FieldSpace space) const {
546  return getEntityDataOrder(type, space, data.dataOnEntities[type]);
547  }
548 
549  /** \name Indices */
550 
551  /**@{*/
552 
553  /// \brief get node indices
554  MoFEMErrorCode getNodesIndices(const boost::string_ref field_name,
556  VectorInt &nodes_indices,
557  VectorInt &local_nodes_indices) const;
558 
559  /// \brief get row node indices from FENumeredDofEntity_multiIndex
560  MoFEMErrorCode getRowNodesIndices(DataForcesAndSourcesCore &data,
561  const std::string &field_name) const;
562 
563  /// \brief get col node indices from FENumeredDofEntity_multiIndex
564  MoFEMErrorCode getColNodesIndices(DataForcesAndSourcesCore &data,
565  const std::string &field_name) const;
566 
567  MoFEMErrorCode getEntityIndices(
568  DataForcesAndSourcesCore &data, const std::string &field_name,
569  FENumeredDofEntity_multiIndex &dofs, const EntityType type_lo = MBVERTEX,
570  const EntityType type_hi = MBPOLYHEDRON) const;
571 
572  inline MoFEMErrorCode
573  getEntityRowIndices(DataForcesAndSourcesCore &data,
574  const std::string &field_name,
575  const EntityType type_lo = MBVERTEX,
576  const EntityType type_hi = MBPOLYHEDRON) const;
577 
578  inline MoFEMErrorCode
579  getEntityColIndices(DataForcesAndSourcesCore &data,
580  const std::string &field_name,
581  const EntityType type_lo = MBVERTEX,
582  const EntityType type_hi = MBPOLYHEDRON) const;
583 
584  /// \brief get NoField indices
585  MoFEMErrorCode getNoFieldIndices(const std::string &field_name,
587  VectorInt &nodes_indices) const;
588 
589  /// \brief get col NoField indices
590  MoFEMErrorCode getNoFieldRowIndices(DataForcesAndSourcesCore &data,
591  const std::string &field_name) const;
592 
593  /// \brief get col NoField indices
594  MoFEMErrorCode getNoFieldColIndices(DataForcesAndSourcesCore &data,
595  const std::string &field_name) const;
596 
597  /**@}*/
598 
599  /** \name Data */
600 
601  /**@{*/
602 
603  /**
604  * \brief Get field data on nodes
605  */
606  MoFEMErrorCode getNoFieldFieldData(const boost::string_ref field_name,
608  VectorDouble &ent_field_data,
609  VectorDofs &ent_field_dofs) const;
610 
611  MoFEMErrorCode getNoFieldFieldData(DataForcesAndSourcesCore &data,
612  const boost::string_ref field_name) const;
613 
614  /**
615  * \brief Get data on nodes
616  * @param data Data structure
617  * @param field_name Field name
618  * @return Error code
619  */
620  MoFEMErrorCode getNodesFieldData(DataForcesAndSourcesCore &data,
621  const std::string &field_name) const;
622 
624  getEntityFieldData(DataForcesAndSourcesCore &data,
625  const std::string &field_name,
626  const EntityType type_lo = MBVERTEX,
627  const EntityType type_hi = MBPOLYHEDRON) const;
628 
629  /**@}*/
630 
631  /// \brief Get nodes on triangles
632  MoFEMErrorCode getFaceTriNodes(DataForcesAndSourcesCore &data) const;
633 
634  /// \brief Get field approximation space and base on entities
636  getSpacesAndBaseOnEntities(DataForcesAndSourcesCore &data) const;
637 
638  /** \name Data form NumeredDofEntity_multiIndex */
639 
640  /**@{*/
641 
642  /// \brief get indices of nodal indices which are declared for problem but
643  /// not this particular element
644  MoFEMErrorCode getProblemNodesIndices(const std::string &field_name,
645  const NumeredDofEntity_multiIndex &dofs,
646  VectorInt &nodes_indices) const;
647 
648  /// \brief get indices by type (generic function) which are declared for
649  /// problem but not this particular element
650  MoFEMErrorCode getProblemTypeIndices(const std::string &field_name,
651  const NumeredDofEntity_multiIndex &dofs,
652  EntityType type, int side_number,
653  VectorInt &indices) const;
654 
655  MoFEMErrorCode getProblemNodesRowIndices(const std::string &field_name,
656  VectorInt &nodes_indices) const;
657  MoFEMErrorCode getProblemTypeRowIndices(const std::string &field_name,
658  EntityType type, int side_number,
659  VectorInt &indices) const;
660  MoFEMErrorCode getProblemNodesColIndices(const std::string &field_name,
661  VectorInt &nodes_indices) const;
662  MoFEMErrorCode getProblemTypeColIndices(const std::string &field_name,
663  EntityType type, int side_number,
664  VectorInt &indices) const;
665 
666  /**@}*/
667 
668  /**
669  * \brief another variant of getRule
670  * @param order_row order of base function on row
671  * @param order_col order of base function on columns
672  * @param order_data order of base function approximating data
673  * @return integration rule
674  *
675  * This function is overloaded by the user. The integration rule
676  * is set such that specific operator implemented by the user is
677  * integrated accurately. For example if user implement bilinear operator
678  * \f[
679  * b(u,v) =
680  * \int_\mathcal{T}
681  * \frac{\partial u_i}{\partial x_j}\frac{\partial v_i}{\partial x_j}
682  * \textrm{d}\mathcal{T}
683  * \f]
684  * then if \f$u\f$ and \f$v\f$ are polynomial of given \em order, then
685  * exact integral would be \code int getRule(int order) { return
686  * 2*(order-1); }; \endcode
687  *
688  * The integration points and weights are set appropriately for given
689  * entity type and integration rule from \ref quad.c
690  *
691  * Method \ref ForcesAndSourcesCore::getRule takes at argument takes
692  * maximal polynomial order set on the element on all fields defined on
693  * the element. If a user likes to have more control, another variant of
694  * this function can be called which distinguishing between field orders
695  * on rows, columns and data, the i.e. first argument of a bilinear form,
696  * the second argument of bilinear form and field coefficients on the
697  * element.
698  *
699  * \note If user set rule to -1 or any other negative integer, then method
700  * \ref ForcesAndSourcesCore::setGaussPts is called. In that method user
701  * can implement own (specific) integration method.
702  *
703  * \bug this function should be const
704  */
705  virtual int getRule(int order_row, int order_col, int order_data);
706 
707  /** \brief set user specific integration rule
708 
709  This function allows for user defined integration rule. The key is to
710  called matrix gaussPts, which is used by other MoFEM procedures. Matrix
711  has number of rows equal to problem dimension plus one, where last index
712  is used to store weight values. %Number of columns is equal to number of
713  integration points.
714 
715  \note This function is called if method \ref
716  ForcesAndSourcesCore::getRule is returning integer -1 or any other
717  negative integer.
718 
719  User sets
720  \code
721  MatrixDouble gaussPts;
722  \endcode
723  where
724  \code
725  gaussPts.resize(dim+1,nb_gauss_pts);
726  \endcode
727  number rows represents local coordinates of integration points
728  in reference element, where last index in row is for integration weight.
729 
730  */
731  virtual MoFEMErrorCode setGaussPts(int order_row, int order_col,
732  int order_data);
733 
734  /**
735  * \brief Calculate base functions
736  * @return Error code
737  */
739  calHierarchicalBaseFunctionsOnElement(const FieldApproximationBase b);
740 
741  /**
742  * \brief Calculate base functions
743  * @return Error code
744  */
745  MoFEMErrorCode calHierarchicalBaseFunctionsOnElement();
746 
747  /**
748  * @brief Calculate Bernstein-Bezier base
749  *
750  * @return MoFEMErrorCode
751  */
752  MoFEMErrorCode calBernsteinBezierBaseFunctionsOnElement();
753 
754  /**
755  * @brief Create a entity data on element object
756  *
757  * @return MoFEMErrorCode
758  */
759  MoFEMErrorCode createDataOnElement();
760 
761  /**
762  * @brief Iterate user data operators
763  *
764  * @return MoFEMErrorCode
765  */
766  MoFEMErrorCode loopOverOperators();
767 
768  /**@{*/
769 
770  /** \name Deprecated (do not use) */
771 
772  /** \deprecated Use getRule(int row_order, int col_order, int data order)
773  */
774  virtual int getRule(int order);
775 
776  /** \deprecated setGaussPts(int row_order, int col_order, int data order);
777  */
778  virtual MoFEMErrorCode setGaussPts(int order);
779 
780  /**@/}*/
781 
782  /**
783  * @brief Entity data on element entity rows fields
784  *
785  */
786  const std::array<boost::shared_ptr<DataForcesAndSourcesCore>, LASTSPACE>
788 
789  /**
790  * @brief Entity data on element entity columns fields
791  *
792  */
793  const std::array<boost::shared_ptr<DataForcesAndSourcesCore>, LASTSPACE>
795 
801 
802  /**
803  * @brief Vector of finite element users data operators
804  *
805  */
806  boost::ptr_vector<UserDataOperator> opPtrVector;
807 
808  friend class UserDataOperator;
809 
810 protected:
811  /**
812  * @brief Last evaluated type of element entity
813  *
814  */
816 
817 private:
818  /**
819  * @brief Pointer to entity polynomial base
820  *
821  */
822  boost::shared_ptr<BaseFunction> elementPolynomialBasePtr;
823 
824  /**
825  * @brief Pointer to user polynomail base
826  */
827  boost::shared_ptr<BaseFunction> userPolynomialBasePtr;
828 
829  /**
830  * @brief Element to integrate on the sides
831  *
832  */
834 
835  /**
836  * @brief Set the pointer to face element on the side
837  *
838  * \note Function is is used by face element, while it iterates over
839  * elements on the side
840  *
841  * @param side_fe_ptr
842  * @return MoFEMErrorCode
843  */
844  MoFEMErrorCode setSideFEPtr(const ForcesAndSourcesCore *side_fe_ptr);
845 
848 };
849 
850 /// \deprecated Used ForcesAndSourcesCore instead
852 
853 MoFEMErrorCode ForcesAndSourcesCore::getEntityRowIndices(
854  DataForcesAndSourcesCore &data, const std::string &field_name,
855  const EntityType type_lo, const EntityType type_hi) const {
856  return getEntityIndices(data, field_name,
857  const_cast<FENumeredDofEntity_multiIndex &>(
858  numeredEntFiniteElementPtr->getRowsDofs()),
859  type_lo, type_hi);
860 }
861 
862 MoFEMErrorCode ForcesAndSourcesCore::getEntityColIndices(
863  DataForcesAndSourcesCore &data, const std::string &field_name,
864  const EntityType type_lo, const EntityType type_hi) const {
865  return getEntityIndices(data, field_name,
866  const_cast<FENumeredDofEntity_multiIndex &>(
867  numeredEntFiniteElementPtr->getColsDofs()),
868  type_lo, type_hi);
869 }
870 
871 boost::shared_ptr<const NumeredEntFiniteElement>
872 ForcesAndSourcesCore::UserDataOperator::getNumeredEntFiniteElementPtr() const {
873  return ptrFE->numeredEntFiniteElementPtr;
874 };
875 
876 EntityHandle ForcesAndSourcesCore::UserDataOperator::getFEEntityHandle() const {
877  return getNumeredEntFiniteElementPtr()->getEnt();
878 }
879 
880 boost::weak_ptr<SideNumber>
881 ForcesAndSourcesCore::UserDataOperator::getSideNumberPtr(
882  const int side_number, const EntityType type) {
883  auto &side_table_by_side_and_type =
884  ptrFE->numeredEntFiniteElementPtr->getSideNumberTable().get<1>();
885  auto side_it =
886  side_table_by_side_and_type.find(boost::make_tuple(type, side_number));
887  if (side_it != side_table_by_side_and_type.end())
888  return *side_it;
889  else
890  return boost::weak_ptr<SideNumber>();
891 }
892 
894 ForcesAndSourcesCore::UserDataOperator::getSideEntity(const int side_number,
895  const EntityType type) {
896  if (auto side_ptr = getSideNumberPtr(side_number, type).lock())
897  return side_ptr->ent;
898  else
899  return 0;
900 }
901 
902 int ForcesAndSourcesCore::UserDataOperator::getNumberOfNodesOnElement() {
903  int num_nodes;
904  CHKERR ptrFE->getNumberOfNodes(num_nodes);
905  return num_nodes;
906 }
907 
908 const FEMethod *ForcesAndSourcesCore::UserDataOperator::getFEMethod() const {
909  return ptrFE;
910 }
911 
912 int ForcesAndSourcesCore::UserDataOperator::getOpType() const { return opType; }
913 
914 void ForcesAndSourcesCore::UserDataOperator::setOpType(const OpType type) {
915  opType = type;
916 }
917 
918 void ForcesAndSourcesCore::UserDataOperator::addOpType(const OpType type) {
919  opType |= type;
920 }
921 
922 int ForcesAndSourcesCore::UserDataOperator::getNinTheLoop() const {
923  return getFEMethod()->getNinTheLoop();
924 }
925 
926 int ForcesAndSourcesCore::UserDataOperator::getLoopSize() const {
927  return getFEMethod()->getLoopSize();
928 }
929 
930 const std::string &ForcesAndSourcesCore::UserDataOperator::getFEName() const {
931  return getFEMethod()->feName;
932 }
933 
934 const PetscData::Switches &
935 ForcesAndSourcesCore::UserDataOperator::getDataCtx() const {
936  return getFEMethod()->data_ctx;
937 }
938 
940 ForcesAndSourcesCore::UserDataOperator::getKSPCtx() const {
941  return getFEMethod()->ksp_ctx;
942 }
943 
945 ForcesAndSourcesCore::UserDataOperator::getSNESCtx() const {
946  return getFEMethod()->snes_ctx;
947 }
948 
950 ForcesAndSourcesCore::UserDataOperator::getTSCtx() const {
951  return getFEMethod()->ts_ctx;
952 }
953 
954 Vec ForcesAndSourcesCore::UserDataOperator::getKSPf() const {
955  return getFEMethod()->ksp_f;
956 }
957 
958 Mat ForcesAndSourcesCore::UserDataOperator::getKSPA() const {
959  return getFEMethod()->ksp_A;
960 }
961 
962 Mat ForcesAndSourcesCore::UserDataOperator::getKSPB() const {
963  return getFEMethod()->ksp_B;
964 }
965 
966 Vec ForcesAndSourcesCore::UserDataOperator::getSNESf() const {
967  return getFEMethod()->snes_f;
968 }
969 
970 Vec ForcesAndSourcesCore::UserDataOperator::getSNESx() const {
971  return getFEMethod()->snes_x;
972 }
973 
974 Mat ForcesAndSourcesCore::UserDataOperator::getSNESA() const {
975  return getFEMethod()->snes_A;
976 }
977 
978 Mat ForcesAndSourcesCore::UserDataOperator::getSNESB() const {
979  return getFEMethod()->snes_B;
980 }
981 
982 Vec ForcesAndSourcesCore::UserDataOperator::getTSu() const {
983  return getFEMethod()->ts_u;
984 }
985 
986 Vec ForcesAndSourcesCore::UserDataOperator::getTSu_t() const {
987  return getFEMethod()->ts_u_t;
988 }
989 
990 Vec ForcesAndSourcesCore::UserDataOperator::getTSu_tt() const {
991  return getFEMethod()->ts_u_tt;
992 }
993 
994 Vec ForcesAndSourcesCore::UserDataOperator::getTSf() const {
995  return getFEMethod()->ts_F;
996 }
997 
998 Mat ForcesAndSourcesCore::UserDataOperator::getTSA() const {
999  return getFEMethod()->ts_A;
1000 }
1001 
1002 Mat ForcesAndSourcesCore::UserDataOperator::getTSB() const {
1003  return getFEMethod()->ts_B;
1004 }
1005 
1006 int ForcesAndSourcesCore::UserDataOperator::getTSstep() const {
1007  return getFEMethod()->ts_step;
1008 }
1009 
1010 double ForcesAndSourcesCore::UserDataOperator::getTStime() const {
1011  return getFEMethod()->ts_t;
1012 }
1013 
1014 double ForcesAndSourcesCore::UserDataOperator::getTSa() const {
1015  return getFEMethod()->ts_a;
1016 }
1017 
1018 MatrixDouble &ForcesAndSourcesCore::UserDataOperator::getGaussPts() {
1019  return static_cast<ForcesAndSourcesCore *>(ptrFE)->gaussPts;
1020 }
1021 
1022 auto ForcesAndSourcesCore::UserDataOperator::getFTensor0IntegrationWeight() {
1024  &(getGaussPts()(getGaussPts().size1() - 1, 0)));
1025 }
1026 
1027 MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::getPorblemRowIndices(
1028  const std::string filed_name, const EntityType type, const int side,
1029  VectorInt &indices) const {
1030  return getProblemRowIndices(filed_name, type, side, indices);
1031 }
1032 
1033 ForcesAndSourcesCore *ForcesAndSourcesCore::UserDataOperator::getPtrFE() const {
1034  return ptrFE;
1035 }
1036 
1038 ForcesAndSourcesCore::UserDataOperator::getSidePtrFE() const {
1039  return ptrFE->sidePtrFE;
1040 }
1041 
1042 } // namespace MoFEM
1043 
1044 #endif //__FORCES_AND_SOURCES_CORE__HPP__
1045 
1046 /**
1047  * \defgroup mofem_forces_and_sources Forces and sources
1048  * \ingroup mofem
1049  *
1050  * \brief Manages complexities related to assembly of vector and matrices at
1051  * single finite element level.
1052  *
1053  **/
MoFEMErrorCode getEntityDataOrder(DataForcesAndSourcesCore &data, const FieldSpace space) const
Get the entity data order for given space.
DataForcesAndSourcesCore & dataL2
structure for User Loop Methods on finite elementsIt can be used to calculate stiffness matrices,...
Deprecated interface functions.
boost::ptr_vector< UserDataOperator > opPtrVector
Vector of finite element users data operators.
MatrixDouble gaussPts
Matrix of integration points.
base operator to do operations at Gauss Pt. level
ForcesAndSourcesCore * sidePtrFE
Element to integrate on the sides.
data structure for finite element entityIt keeps that about indices of degrees of freedom,...
DataForcesAndSourcesCore & dataH1
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
std::bitset< 8 > Switches
Definition: LoopMethods.hpp:59
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
auto & getUserPolynomialBase()
Get the User Polynomial Base object.
multi_index_container< boost::shared_ptr< FENumeredDofEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< FENumeredDofEntity::interface_type_DofEntity, const UId &, &FENumeredDofEntity::getGlobalUniqueId > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< FENumeredDofEntity::interface_type_DofEntity, EntityHandle, &FENumeredDofEntity::getEnt > >, ordered_non_unique< tag< FieldName_mi_tag >, const_mem_fun< FENumeredDofEntity::interface_type_Field, boost::string_ref, &FENumeredDofEntity::getNameRef > >, ordered_non_unique< tag< Composite_Name_Type_And_Side_Number_mi_tag >, composite_key< FENumeredDofEntity, const_mem_fun< FENumeredDofEntity::interface_type_Field, boost::string_ref, &FENumeredDofEntity::getNameRef >, const_mem_fun< FENumeredDofEntity::interface_type_RefEntity, EntityType, &FENumeredDofEntity::getEntType >, KeyFromKey< member< SideNumber, char, &SideNumber::side_number >, member< FENumeredDofEntity::BaseFEEntity, boost::shared_ptr< SideNumber >, &FENumeredDofEntity::sideNumberPtr > > > >, ordered_non_unique< tag< Composite_Name_And_Type_mi_tag >, composite_key< FENumeredDofEntity, const_mem_fun< FENumeredDofEntity::interface_type_Field, boost::string_ref, &FENumeredDofEntity::getNameRef >, const_mem_fun< FENumeredDofEntity::interface_type_RefEntity, EntityType, &FENumeredDofEntity::getEntType > > >, ordered_non_unique< tag< Composite_Name_And_Ent_mi_tag >, composite_key< FENumeredDofEntity, const_mem_fun< FENumeredDofEntity::interface_type_Field, boost::string_ref, &FENumeredDofEntity::getNameRef >, const_mem_fun< FENumeredDofEntity::interface_type_DofEntity, EntityHandle, &FENumeredDofEntity::getEnt > > > > > FENumeredDofEntity_multiIndex
MultiIndex container keeps FENumeredDofEntity.
boost::shared_ptr< BaseFunction > userPolynomialBasePtr
Pointer to user polynomail base.
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
OpType
Controls loop over entities on element.
FieldSpace in [ 0, LASTSPACE )
Definition: definitions.h:181
Edge finite elementUser is implementing own operator at Gauss points level, by own object derived fro...
RuleHookFun getRuleHook
Hook to get rule.
const int dim
DataForcesAndSourcesCore::EntData & getEntData(const FieldSpace space, const EntityType type, const int side)
Get the entity data.
DataForcesAndSourcesCore & dataHcurl
EntityType lastEvaluatedElementEntityType
Last evaluated type of element entity.
const std::array< boost::shared_ptr< DataForcesAndSourcesCore >, LASTSPACE > derivedDataOnElement
Entity data on element entity columns fields.
ForcesAndSourcesCore::UserDataOperator UserDataOperator
FieldApproximationBase
approximation base
Definition: definitions.h:150
boost::shared_ptr< BaseFunction > elementPolynomialBasePtr
Pointer to entity polynomial base.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp: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:174
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
#define CHKERR
Inline error check.
Definition: definitions.h:602
constexpr int order
multi_index_container< boost::shared_ptr< FEDofEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< FEDofEntity::interface_type_DofEntity, const UId &, &FEDofEntity::getGlobalUniqueId > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< FEDofEntity::interface_type_DofEntity, EntityHandle, &FEDofEntity::getEnt > >, ordered_non_unique< tag< FieldName_mi_tag >, const_mem_fun< FEDofEntity::interface_type_Field, boost::string_ref, &FEDofEntity::getNameRef > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< FEDofEntity::interface_type_RefEntity, EntityType, &FEDofEntity::getEntType > >, ordered_non_unique< tag< Composite_Name_And_Type_mi_tag >, composite_key< FEDofEntity, const_mem_fun< FEDofEntity::interface_type_Field, boost::string_ref, &FEDofEntity::getNameRef >, const_mem_fun< FEDofEntity::interface_type_RefEntity, EntityType, &FEDofEntity::getEntType > > >, ordered_non_unique< tag< Composite_Name_And_Ent_mi_tag >, composite_key< FEDofEntity, const_mem_fun< FEDofEntity::interface_type_Field, boost::string_ref, &FEDofEntity::getNameRef >, const_mem_fun< FEDofEntity::interface_type_DofEntity, EntityHandle, &FEDofEntity::getEnt > > >, ordered_non_unique< tag< Composite_Name_Type_And_Side_Number_mi_tag >, composite_key< FEDofEntity, const_mem_fun< FEDofEntity::interface_type_Field, boost::string_ref, &FEDofEntity::getNameRef >, const_mem_fun< FEDofEntity::interface_type_RefEntity, EntityType, &FEDofEntity::getEntType >, KeyFromKey< member< SideNumber, char, &SideNumber::side_number >, member< FEDofEntity::BaseFEEntity, boost::shared_ptr< SideNumber >, &FEDofEntity::sideNumberPtr > > > > > > FEDofEntity_multiIndex
MultiIndex container keeps FEDofEntity.
const DataForcesAndSourcesCore::EntData & getEntData(const FieldSpace space, const EntityType type, const int side) const
Get the entity data.
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:72
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:73
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...
const std::array< boost::shared_ptr< DataForcesAndSourcesCore >, LASTSPACE > dataOnElement
Entity data on element entity rows fields.
RuleHookFun setRuleHook
Set function to calculate integration rule.
auto & getElementPolynomialBase()
Get the Entity Polynomial Base object.
multi_index_container< boost::shared_ptr< NumeredDofEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, const UId &, &NumeredDofEntity::getGlobalUniqueId > >, ordered_non_unique< tag< Unique_Ent_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, const UId &, &NumeredDofEntity::getEntGlobalUniqueId > >, ordered_non_unique< tag< Part_mi_tag >, member< NumeredDofEntity, unsigned int, &NumeredDofEntity::pArt > >, ordered_non_unique< tag< Idx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::dofIdx > >, ordered_non_unique< tag< FieldName_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_Field, boost::string_ref, &NumeredDofEntity::getNameRef > >, ordered_non_unique< tag< PetscGlobalIdx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::petscGloablDofIdx > >, ordered_non_unique< tag< PetscLocalIdx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::petscLocalDofIdx > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, EntityHandle, &NumeredDofEntity::getEnt > >, ordered_non_unique< tag< Composite_Name_And_Ent_And_EntDofIdx_mi_tag >, composite_key< NumeredDofEntity, const_mem_fun< NumeredDofEntity::interface_type_Field, boost::string_ref, &NumeredDofEntity::getNameRef >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, EntityHandle, &NumeredDofEntity::getEnt >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, DofIdx, &NumeredDofEntity::getEntDofIdx > > >, ordered_non_unique< tag< Composite_Name_And_Part_mi_tag >, composite_key< NumeredDofEntity, const_mem_fun< NumeredDofEntity::interface_type_Field, boost::string_ref, &NumeredDofEntity::getNameRef >, member< NumeredDofEntity, unsigned int, &NumeredDofEntity::pArt > > >, ordered_non_unique< tag< Composite_Name_Ent_And_Part_mi_tag >, composite_key< NumeredDofEntity, const_mem_fun< NumeredDofEntity::interface_type_Field, boost::string_ref, &NumeredDofEntity::getNameRef >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, EntityHandle, &NumeredDofEntity::getEnt >, member< NumeredDofEntity, unsigned int, &NumeredDofEntity::pArt > > > > > NumeredDofEntity_multiIndex
MultiIndex container keeps NumeredDofEntity.
Face finite elementUser is implementing own operator at Gauss point level, by own object derived from...
DEPRECATED typedef ForcesAndSourcesCore ForcesAndSurcesCore
DataForcesAndSourcesCore & dataNoField