v0.8.23
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 SNES */
268 
269  /**@{*/
270 
271  inline Vec getSnesF() const;
272 
273  inline Vec getSnesX() const;
274 
275  inline Mat getSnesA() const;
276 
277  inline Mat getSnesB() const;
278 
279  /**@}*/
280 
281  /** \name Accessing TS */
282 
283  /**@{*/
284 
285  inline Vec getTSu() const;
286 
287  inline Vec getTSu_t() const;
288 
289  inline Vec getTSf() const;
290 
291  inline Mat getTSA() const;
292 
293  inline Mat getTSB() const;
294 
295  inline int getTSstep() const;
296 
297  inline double getTStime() const;
298 
299  inline double getTSa() const;
300 
301  /**@}*/
302 
303  /**@{*/
304 
305  /** \name Base funtions and integration points */
306 
307  /** \brief matrix of integration (Gauss) points for Volume Element
308  *
309  * For triangle: columns 0,1 are x,y coordinates respectively and column
310  * 2 is a weight value for example getGaussPts()(1,13) returns y
311  * coordinate of 13th Gauss point on particular volume element
312  *
313  * For tetrahedron: columns 0,1,2 are x,y,z coordinates respectively and
314  * column 3 is a weight value for example getGaussPts()(1,13) returns y
315  * coordinate of 13th Gauss point on particular volume element
316  *
317  */
318  inline MatrixDouble &getGaussPts();
319 
320  /**
321  * @brief Get integration weights
322  *
323  * \code
324  * auto t_w = getFTensor0IntegrationWeight();
325  * for(int gg = 0; gg!=getGaussPts.size2(); ++gg) {
326  * // integrate something
327  * ++t_w;
328  * }
329  * \endcode
330  *
331  * @return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>
332  */
333  inline auto getFTensor0IntegrationWeight();
334 
335  /**@}*/
336 
337  /**@{*/
338 
339  /** \name Deprecated (do not use) */
340 
341  // \deprecated Deprecated function with spelling mistake
343  getPorblemRowIndices(const std::string filed_name, const EntityType type,
344  const int side, VectorInt &indices) const;
345 
346  /**@}*/
347 
348  protected:
350 
351  inline MoFEMErrorCode setPtrFE(ForcesAndSourcesCore *ptr);
352 
353  inline ForcesAndSourcesCore *getPtrFE() const;
354 
355  inline ForcesAndSourcesCore *getSidePtrFE() const;
356 
357  private:
358  /**
359  * @brief User call this function to loop over elements on the side of
360  * face. This function calls finite element with is operator to do
361  * calculations.
362  *
363  * @param fe_name
364  * @param side_fe
365  * @param dim
366  * @return MoFEMErrorCode
367  */
368  MoFEMErrorCode loopSide(const string &fe_name,
369  ForcesAndSourcesCore *side_fe, const size_t dim);
370 
374  };
375 
376  /** \brief Use to push back operator for row operator
377 
378  It can be used to calculate nodal forces or other quantities on the mesh.
379 
380  */
381  boost::ptr_vector<UserDataOperator> &getOpPtrVector() { return opPtrVector; }
382 
383  /**
384  * @brief Get the Entity Polynomial Base object
385  *
386  * @return boost::shared_ptr<BaseFunction>&&
387  */
388  auto &getElementPolynomialBase() { return elementPolynomialBasePtr; }
389 
390  /**
391  * @brief Get the User Polynomial Base object
392  *
393  * @return boost::shared_ptr<BaseFunction>&
394  */
395  auto &getUserPolynomialBase() { return userPolynomialBasePtr; }
396 
397  /**
398  * @brief Matrix of integration points
399  *
400  * Columns is equal to number of integration points, numver of rows
401  * depends on dimension of finite element entity, for example for
402  * tetrahedron rows are x,y,z,weight. Last row is integration weight.
403  *
404  * FIXME: that should be moved to private class data and acessed only by
405  * member function
406  */
408 
409  virtual MoFEMErrorCode preProcess();
410  virtual MoFEMErrorCode operator()();
411  virtual MoFEMErrorCode postProcess();
412 
413 public:
414  /** \brief Get max order of approximation for data fields
415 
416  Method getMaxDataOrder () return maximal order on entities, for
417  all data on the element. So for example if finite element is triangle, and
418  triangle base function have order 4 and on edges base function have order
419  2, this function return 4.
420 
421  If finite element has for example 2 or more approximated fields, for
422  example Pressure (order 3) and displacement field (order 5), this function
423  returns 5.
424 
425  */
426  int getMaxDataOrder() const;
427 
428  /// \brief Get max order of approximation for field in rows
429  int getMaxRowOrder() const;
430 
431  /// \brief Get max order of approximation for field in columns
432  int getMaxColOrder() const;
433 
434  /// \brief Get number of DOFs on element
435  MoFEMErrorCode getNumberOfNodes(int &num_nodes) const;
436 
437  /**
438  * @brief Get the entity data
439  *
440  * @param space
441  * @param type
442  * @param side
443  * @return const DataForcesAndSourcesCore::EntData&
444  */
446  const EntityType type,
447  const int side) const {
448  return dataOnElement[space]->dataOnEntities[type][side];
449  }
450 
451  /**
452  * @brief Get the entity data
453  *
454  * @param space
455  * @param type
456  * @param side
457  * @return DataForcesAndSourcesCore::EntData&
458  */
460  getEntData(const FieldSpace space, const EntityType type, const int side) {
461  return dataOnElement[space]->dataOnEntities[type][side];
462  }
463 
464 protected:
465  /**
466  * \brief get sense (orientation) of entity
467  * @param type type of entity
468  * @param data entity data
469  * @return error code
470  */
471  MoFEMErrorCode getEntitySense(
472  const EntityType type,
473  boost::ptr_vector<DataForcesAndSourcesCore::EntData> &data) const;
474 
475  /**
476  * @brief Get the entity data order
477  *
478  * @param type
479  * @param space
480  * @param data
481  * @return MoFEMErrorCode
482  */
483  MoFEMErrorCode getEntityDataOrder(
484  const EntityType type, const FieldSpace space,
485  boost::ptr_vector<DataForcesAndSourcesCore::EntData> &data) const;
486 
487  /**
488  * @brief Get the entity sense (orientation)
489  *
490  * @tparam type
491  * @param data
492  * @return MoFEMErrorCode
493  */
494  template <EntityType type>
496  return getEntitySense(type, data.dataOnEntities[type]);
497  }
498 
499  /**
500  * @brief Get the entity data order for given space
501  *
502  * @tparam type
503  * @param data
504  * @param space
505  * @return MoFEMErrorCode
506  */
507  template <EntityType type>
509  const FieldSpace space) const {
510  return getEntityDataOrder(type, space, data.dataOnEntities[type]);
511  }
512 
513  /** \name Indices */
514 
515  /**@{*/
516 
517  /// \brief get node indices
518  MoFEMErrorCode getNodesIndices(const boost::string_ref field_name,
520  VectorInt &nodes_indices,
521  VectorInt &local_nodes_indices) const;
522 
523  /// \brief get row node indices from FENumeredDofEntity_multiIndex
524  MoFEMErrorCode getRowNodesIndices(DataForcesAndSourcesCore &data,
525  const std::string &field_name) const;
526 
527  /// \brief get col node indices from FENumeredDofEntity_multiIndex
528  MoFEMErrorCode getColNodesIndices(DataForcesAndSourcesCore &data,
529  const std::string &field_name) const;
530 
531  MoFEMErrorCode getEntityIndices(
532  DataForcesAndSourcesCore &data, const std::string &field_name,
533  FENumeredDofEntity_multiIndex &dofs, const EntityType type_lo = MBVERTEX,
534  const EntityType type_hi = MBPOLYHEDRON) const;
535 
536  inline MoFEMErrorCode
537  getEntityRowIndices(DataForcesAndSourcesCore &data,
538  const std::string &field_name,
539  const EntityType type_lo = MBVERTEX,
540  const EntityType type_hi = MBPOLYHEDRON) const;
541 
542  inline MoFEMErrorCode
543  getEntityColIndices(DataForcesAndSourcesCore &data,
544  const std::string &field_name,
545  const EntityType type_lo = MBVERTEX,
546  const EntityType type_hi = MBPOLYHEDRON) const;
547 
548  /// \brief get NoField indices
549  MoFEMErrorCode getNoFieldIndices(const std::string &field_name,
551  VectorInt &nodes_indices) const;
552 
553  /// \brief get col NoField indices
554  MoFEMErrorCode getNoFieldRowIndices(DataForcesAndSourcesCore &data,
555  const std::string &field_name) const;
556 
557  /// \brief get col NoField indices
558  MoFEMErrorCode getNoFieldColIndices(DataForcesAndSourcesCore &data,
559  const std::string &field_name) const;
560 
561  /**@}*/
562 
563  /** \name Data */
564 
565  /**@{*/
566 
567  /**
568  * \brief Get field data on nodes
569  * @param field_name Name of field
570  * @param dofs Dofs (element) multi index
571  * @param nodes_data Returned DOFs values
572  * @param nodes_dofs Vector of pointers to DOFs data structure
573  * @param space Get space on nodes (Only H! is valid)
574  * @param base Get base on nodes
575  * @return Error code
576  */
577  MoFEMErrorCode getNodesFieldData(const boost::string_ref field_name,
579  VectorDouble &nodes_data,
580  VectorDofs &nodes_dofs, FieldSpace &space,
581  FieldApproximationBase &base) const;
582 
583  MoFEMErrorCode getNoFieldFieldData(const boost::string_ref field_name,
585  VectorDouble &ent_field_data,
586  VectorDofs &ent_field_dofs) const;
587  MoFEMErrorCode getNoFieldFieldData(DataForcesAndSourcesCore &data,
588  const boost::string_ref field_name) const;
589 
590  /**
591  * \brief Get data on nodes
592  * @param data Data structure
593  * @param field_name Field name
594  * @return Error code
595  */
596  MoFEMErrorCode getNodesFieldData(DataForcesAndSourcesCore &data,
597  const std::string &field_name) const;
598 
600  getEntityFieldData(DataForcesAndSourcesCore &data,
601  const std::string &field_name,
602  const EntityType type_lo = MBVERTEX,
603  const EntityType type_hi = MBPOLYHEDRON) const;
604 
605  /**@}*/
606 
607  /// \brief Get nodes on triangles
608  MoFEMErrorCode getFaceTriNodes(DataForcesAndSourcesCore &data) const;
609 
610  /// \brief Get field approximation space and base on entities
612  getSpacesAndBaseOnEntities(DataForcesAndSourcesCore &data) const;
613 
614  /** \name Data form NumeredDofEntity_multiIndex */
615 
616  /**@{*/
617 
618  /// \brief get indices of nodal indices which are declared for problem but
619  /// not this particular element
620  MoFEMErrorCode getProblemNodesIndices(const std::string &field_name,
621  const NumeredDofEntity_multiIndex &dofs,
622  VectorInt &nodes_indices) const;
623 
624  /// \brief get indices by type (generic function) which are declared for
625  /// problem but not this particular element
626  MoFEMErrorCode getProblemTypeIndices(const std::string &field_name,
627  const NumeredDofEntity_multiIndex &dofs,
628  EntityType type, int side_number,
629  VectorInt &indices) const;
630 
631  MoFEMErrorCode getProblemNodesRowIndices(const std::string &field_name,
632  VectorInt &nodes_indices) const;
633  MoFEMErrorCode getProblemTypeRowIndices(const std::string &field_name,
634  EntityType type, int side_number,
635  VectorInt &indices) const;
636  MoFEMErrorCode getProblemNodesColIndices(const std::string &field_name,
637  VectorInt &nodes_indices) const;
638  MoFEMErrorCode getProblemTypeColIndices(const std::string &field_name,
639  EntityType type, int side_number,
640  VectorInt &indices) const;
641 
642  /**@}*/
643 
644  /**
645  * \brief another variant of getRule
646  * @param order_row order of base function on row
647  * @param order_col order of base function on columns
648  * @param order_data order of base function approximating data
649  * @return integration rule
650  *
651  * This function is overloaded by the user. The integration rule
652  * is set such that specific operator implemented by the user is
653  * integrated accurately. For example if user implement bilinear operator
654  * \f[
655  * b(u,v) =
656  * \int_\mathcal{T}
657  * \frac{\partial u_i}{\partial x_j}\frac{\partial v_i}{\partial x_j}
658  * \textrm{d}\mathcal{T}
659  * \f]
660  * then if \f$u\f$ and \f$v\f$ are polynomial of given \em order, then
661  * exact integral would be \code int getRule(int order) { return
662  * 2*(order-1); }; \endcode
663  *
664  * The integration points and weights are set appropriately for given
665  * entity type and integration rule from \ref quad.c
666  *
667  * Method \ref ForcesAndSourcesCore::getRule takes at argument takes
668  * maximal polynomial order set on the element on all fields defined on
669  * the element. If a user likes to have more control, another variant of
670  * this function can be called which distinguishing between field orders
671  * on rows, columns and data, the i.e. first argument of a bilinear form,
672  * the second argument of bilinear form and field coefficients on the
673  * element.
674  *
675  * \note If user set rule to -1 or any other negative integer, then method
676  * \ref ForcesAndSourcesCore::setGaussPts is called. In that method user
677  * can implement own (specific) integration method.
678  *
679  * \bug this function should be const
680  */
681  virtual int getRule(int order_row, int order_col, int order_data);
682 
683  /** \brief set user specific integration rule
684 
685  This function allows for user defined integration rule. The key is to
686  called matrix gaussPts, which is used by other MoFEM procedures. Matrix
687  has number of rows equal to problem dimension plus one, where last index
688  is used to store weight values. %Number of columns is equal to number of
689  integration points.
690 
691  \note This function is called if method \ref
692  ForcesAndSourcesCore::getRule is returning integer -1 or any other
693  negative integer.
694 
695  User sets
696  \code
697  MatrixDouble gaussPts;
698  \endcode
699  where
700  \code
701  gaussPts.resize(dim+1,nb_gauss_pts);
702  \endcode
703  number rows represents local coordinates of integration points
704  in reference element, where last index in row is for integration weight.
705 
706  */
707  virtual MoFEMErrorCode setGaussPts(int order_row, int order_col,
708  int order_data);
709 
710  /**
711  * \brief Calculate base functions
712  * @return Error code
713  */
715  calculateBaseFunctionsOnElement(const FieldApproximationBase b);
716 
717  /**
718  * \brief Calculate base functions
719  * @return Error code
720  */
721  MoFEMErrorCode calculateBaseFunctionsOnElement();
722 
723  /**
724  * @brief Create a entity data on element object
725  *
726  * @return MoFEMErrorCode
727  */
728  MoFEMErrorCode createDataOnElement();
729 
730  /**
731  * @brief Iterate user data operators
732  *
733  * @return MoFEMErrorCode
734  */
735  MoFEMErrorCode loopOverOperators();
736 
737  /**@{*/
738 
739  /** \name Deprecated (do not use) */
740 
741  /** \deprecated Use getRule(int row_order, int col_order, int data order)
742  */
743  virtual int getRule(int order);
744 
745  /** \deprecated setGaussPts(int row_order, int col_order, int data order);
746  */
747  virtual MoFEMErrorCode setGaussPts(int order);
748 
749  /**@/}*/
750 
751  /**
752  * @brief Entity data on element entity rows fields
753  *
754  */
755  const boost::shared_ptr<DataForcesAndSourcesCore> dataOnElement[LASTSPACE];
756 
757  /**
758  * @brief Entity data on element entity columns fields
759  *
760  */
761  const boost::shared_ptr<DataForcesAndSourcesCore>
762  derivedDataOnElement[LASTSPACE];
763 
769 
770  /**
771  * @brief Vector of finite element users data operators
772  *
773  */
774  boost::ptr_vector<UserDataOperator> opPtrVector;
775 
776  friend class UserDataOperator;
777 
778 private:
779  /**
780  * @brief Last evaluated type of element entity
781  *
782  */
784 
785  /**
786  * @brief Pointer to entity polynomial base
787  *
788  */
789  boost::shared_ptr<BaseFunction> elementPolynomialBasePtr;
790 
791  /**
792  * @brief Pointer to user polynomail base
793  */
794  boost::shared_ptr<BaseFunction> userPolynomialBasePtr;
795 
796  /**
797  * @brief Element to integrate on the sides
798  *
799  */
801 
802  /**
803  * @brief Set the pointer to face element on the side
804  *
805  * \note Function is is used by face element, while it iterates over
806  * elements on the side
807  *
808  * @param side_fe_ptr
809  * @return MoFEMErrorCode
810  */
811  MoFEMErrorCode setSideFEPtr(const ForcesAndSourcesCore *side_fe_ptr);
812 
815 };
816 
817 /// \deprecated Used ForcesAndSourcesCore instead
819 
820 MoFEMErrorCode ForcesAndSourcesCore::getEntityRowIndices(
821  DataForcesAndSourcesCore &data, const std::string &field_name,
822  const EntityType type_lo, const EntityType type_hi) const {
823  return getEntityIndices(data, field_name,
824  const_cast<FENumeredDofEntity_multiIndex &>(
825  numeredEntFiniteElementPtr->getRowsDofs()),
826  type_lo, type_hi);
827 }
828 
829 MoFEMErrorCode ForcesAndSourcesCore::getEntityColIndices(
830  DataForcesAndSourcesCore &data, const std::string &field_name,
831  const EntityType type_lo, const EntityType type_hi) const {
832  return getEntityIndices(data, field_name,
833  const_cast<FENumeredDofEntity_multiIndex &>(
834  numeredEntFiniteElementPtr->getColsDofs()),
835  type_lo, type_hi);
836 }
837 
838 boost::shared_ptr<const NumeredEntFiniteElement>
839 ForcesAndSourcesCore::UserDataOperator::getNumeredEntFiniteElementPtr() const {
840  return ptrFE->numeredEntFiniteElementPtr;
841 };
842 
843 EntityHandle ForcesAndSourcesCore::UserDataOperator::getFEEntityHandle() const {
844  return getNumeredEntFiniteElementPtr()->getEnt();
845 }
846 
847 boost::weak_ptr<SideNumber>
848 ForcesAndSourcesCore::UserDataOperator::getSideNumberPtr(
849  const int side_number, const EntityType type) {
850  auto &side_table_by_side_and_type =
851  ptrFE->numeredEntFiniteElementPtr->getSideNumberTable().get<1>();
852  auto side_it =
853  side_table_by_side_and_type.find(boost::make_tuple(type, side_number));
854  if (side_it != side_table_by_side_and_type.end())
855  return *side_it;
856  else
857  return boost::weak_ptr<SideNumber>();
858 }
859 
861 ForcesAndSourcesCore::UserDataOperator::getSideEntity(const int side_number,
862  const EntityType type) {
863  if (auto side_ptr = getSideNumberPtr(side_number, type).lock())
864  return side_ptr->ent;
865  else
866  return 0;
867 }
868 
869 int ForcesAndSourcesCore::UserDataOperator::getNumberOfNodesOnElement() {
870  int num_nodes;
871  CHKERR ptrFE->getNumberOfNodes(num_nodes);
872  return num_nodes;
873 }
874 
875 const FEMethod *ForcesAndSourcesCore::UserDataOperator::getFEMethod() const {
876  return ptrFE;
877 }
878 
879 int ForcesAndSourcesCore::UserDataOperator::getOpType() const { return opType; }
880 
881 void ForcesAndSourcesCore::UserDataOperator::setOpType(const OpType type) {
882  opType = type;
883 }
884 
885 void ForcesAndSourcesCore::UserDataOperator::addOpType(const OpType type) {
886  opType |= type;
887 }
888 
889 int ForcesAndSourcesCore::UserDataOperator::getNinTheLoop() const {
890  return getFEMethod()->getNinTheLoop();
891 }
892 
893 int ForcesAndSourcesCore::UserDataOperator::getLoopSize() const {
894  return getFEMethod()->getLoopSize();
895 }
896 
897 const std::string &ForcesAndSourcesCore::UserDataOperator::getFEName() const {
898  return getFEMethod()->feName;
899 }
900 
901 Vec ForcesAndSourcesCore::UserDataOperator::getSnesF() const {
902  return getFEMethod()->snes_f;
903 }
904 
905 Vec ForcesAndSourcesCore::UserDataOperator::getSnesX() const {
906  return getFEMethod()->snes_x;
907 }
908 
909 Mat ForcesAndSourcesCore::UserDataOperator::getSnesA() const {
910  return getFEMethod()->snes_A;
911 }
912 
913 Mat ForcesAndSourcesCore::UserDataOperator::getSnesB() const {
914  return getFEMethod()->snes_B;
915 }
916 
917 Vec ForcesAndSourcesCore::UserDataOperator::getTSu() const {
918  return getFEMethod()->ts_u;
919 }
920 
921 Vec ForcesAndSourcesCore::UserDataOperator::getTSu_t() const {
922  return getFEMethod()->ts_u_t;
923 }
924 
925 Vec ForcesAndSourcesCore::UserDataOperator::getTSf() const {
926  return getFEMethod()->ts_F;
927 }
928 
929 Mat ForcesAndSourcesCore::UserDataOperator::getTSA() const {
930  return getFEMethod()->ts_A;
931 }
932 
933 Mat ForcesAndSourcesCore::UserDataOperator::getTSB() const {
934  return getFEMethod()->ts_B;
935 }
936 
937 int ForcesAndSourcesCore::UserDataOperator::getTSstep() const {
938  return getFEMethod()->ts_step;
939 }
940 
941 double ForcesAndSourcesCore::UserDataOperator::getTStime() const {
942  return getFEMethod()->ts_t;
943 }
944 
945 double ForcesAndSourcesCore::UserDataOperator::getTSa() const {
946  return getFEMethod()->ts_a;
947 }
948 
949 MatrixDouble &ForcesAndSourcesCore::UserDataOperator::getGaussPts() {
950  return static_cast<ForcesAndSourcesCore *>(ptrFE)->gaussPts;
951 }
952 
953 auto ForcesAndSourcesCore::UserDataOperator::getFTensor0IntegrationWeight() {
955  &(getGaussPts()(getGaussPts().size1() - 1, 0)));
956 }
957 
958 MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::getPorblemRowIndices(
959  const std::string filed_name, const EntityType type, const int side,
960  VectorInt &indices) const {
961  return getProblemRowIndices(filed_name, type, side, indices);
962 }
963 
965 ForcesAndSourcesCore::UserDataOperator::setPtrFE(ForcesAndSourcesCore *ptr) {
967  ptrFE = ptr;
969 }
970 
971 ForcesAndSourcesCore *ForcesAndSourcesCore::UserDataOperator::getPtrFE() const {
972  return ptrFE;
973 }
974 
976 ForcesAndSourcesCore::UserDataOperator::getSidePtrFE() const {
977  return ptrFE->sidePtrFE;
978 }
979 
980 } // namespace MoFEM
981 
982 #endif //__FORCES_AND_SOURCES_CORE__HPP__
983 
984 /**
985  * \defgroup mofem_forces_and_sources Forces and sources
986  * \ingroup mofem
987  *
988  * \brief Manages complexities related to assembly of vector and matrices at
989  * single finite element level.
990  *
991  **/
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,...
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
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
auto & getUserPolynomialBase()
Get the User Polynomial Base object.
multi_index_container< boost::shared_ptr< FENumeredDofEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< FENumeredDofEntity::interface_type_DofEntity, const UId &, &FENumeredDofEntity::getGlobalUniqueId > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< FENumeredDofEntity::interface_type_DofEntity, EntityHandle, &FENumeredDofEntity::getEnt > >, ordered_non_unique< tag< FieldName_mi_tag >, const_mem_fun< FENumeredDofEntity::interface_type_Field, boost::string_ref, &FENumeredDofEntity::getNameRef > >, ordered_non_unique< tag< Composite_Name_Type_And_Side_Number_mi_tag >, composite_key< FENumeredDofEntity, const_mem_fun< FENumeredDofEntity::interface_type_Field, boost::string_ref, &FENumeredDofEntity::getNameRef >, const_mem_fun< FENumeredDofEntity::interface_type_RefEntity, EntityType, &FENumeredDofEntity::getEntType >, KeyFromKey< member< SideNumber, char, &SideNumber::side_number >, member< FENumeredDofEntity::BaseFEEntity, boost::shared_ptr< SideNumber >, &FENumeredDofEntity::sideNumberPtr > > > >, ordered_non_unique< tag< Composite_Name_And_Type_mi_tag >, composite_key< FENumeredDofEntity, const_mem_fun< FENumeredDofEntity::interface_type_Field, boost::string_ref, &FENumeredDofEntity::getNameRef >, const_mem_fun< FENumeredDofEntity::interface_type_RefEntity, EntityType, &FENumeredDofEntity::getEntType > > >, ordered_non_unique< tag< Composite_Name_And_Ent_mi_tag >, composite_key< FENumeredDofEntity, const_mem_fun< FENumeredDofEntity::interface_type_Field, boost::string_ref, &FENumeredDofEntity::getNameRef >, const_mem_fun< FENumeredDofEntity::interface_type_DofEntity, EntityHandle, &FENumeredDofEntity::getEnt > > > > > FENumeredDofEntity_multiIndex
MultiIndex container keeps FENumeredDofEntity.
boost::shared_ptr< BaseFunction > userPolynomialBasePtr
Pointer to user polynomail base.
int getNinTheLoop() const
get number of evaluated element in the loop
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
OpType
Controls loop over entities on element.
boost::ptr_vector< EntData > dataOnEntities[MBMAXTYPE]
FieldSpace in [ 0, LASTSPACE )
Definition: definitions.h:175
Edge finite elementUser is implementing own operator at Gauss points level, by own object derived fro...
RuleHookFun getRuleHook
Hook to get rule.
DataForcesAndSourcesCore::EntData & getEntData(const FieldSpace space, const EntityType type, const int side)
Get the entity data.
DataForcesAndSourcesCore & dataHcurl
EntityType lastEvaluatedElementEntityType
Last evaluated type of element entity.
FieldApproximationBase
approximation base
Definition: definitions.h:143
boost::shared_ptr< BaseFunction > elementPolynomialBasePtr
Pointer to entity polynomial base.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
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:168
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
#define CHKERR
Inline error check.
Definition: definitions.h:596
ForcesAndSourcesCore::UserDataOperator UserDataOperator
multi_index_container< boost::shared_ptr< FEDofEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< FEDofEntity::interface_type_DofEntity, const UId &, &FEDofEntity::getGlobalUniqueId > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< FEDofEntity::interface_type_DofEntity, EntityHandle, &FEDofEntity::getEnt > >, ordered_non_unique< tag< FieldName_mi_tag >, const_mem_fun< FEDofEntity::interface_type_Field, boost::string_ref, &FEDofEntity::getNameRef > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< FEDofEntity::interface_type_RefEntity, EntityType, &FEDofEntity::getEntType > >, ordered_non_unique< tag< Composite_Name_And_Type_mi_tag >, composite_key< FEDofEntity, const_mem_fun< FEDofEntity::interface_type_Field, boost::string_ref, &FEDofEntity::getNameRef >, const_mem_fun< FEDofEntity::interface_type_RefEntity, EntityType, &FEDofEntity::getEntType > > >, ordered_non_unique< tag< Composite_Name_And_Ent_mi_tag >, composite_key< FEDofEntity, const_mem_fun< FEDofEntity::interface_type_Field, boost::string_ref, &FEDofEntity::getNameRef >, const_mem_fun< FEDofEntity::interface_type_DofEntity, EntityHandle, &FEDofEntity::getEnt > > >, ordered_non_unique< tag< Composite_Name_Type_And_Side_Number_mi_tag >, composite_key< FEDofEntity, const_mem_fun< FEDofEntity::interface_type_Field, boost::string_ref, &FEDofEntity::getNameRef >, const_mem_fun< FEDofEntity::interface_type_RefEntity, EntityType, &FEDofEntity::getEntType >, KeyFromKey< member< SideNumber, char, &SideNumber::side_number >, member< FEDofEntity::BaseFEEntity, boost::shared_ptr< SideNumber >, &FEDofEntity::sideNumberPtr > > > > > > FEDofEntity_multiIndex
MultiIndex container keeps FEDofEntity.
const DataForcesAndSourcesCore::EntData & getEntData(const FieldSpace space, const EntityType type, const int side) const
Get the entity data.
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:71
ublas::vector< boost::shared_ptr< const FEDofEntity >, DofsAllocator > VectorDofs
#define DEPRECATED
Definition: definitions.h:29
DataForcesAndSourcesCore & dataHdiv
Data on single entity (This is passed as argument to DataOperator::doWork)
structure to get information form mofem into DataForcesAndSourcesCore
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...
RuleHookFun setRuleHook
Set function to calculate integration rule.
const int order
approximation order
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