v0.15.0
Loading...
Searching...
No Matches
ForcesAndSourcesCore.hpp
Go to the documentation of this file.
1/** \file ForcesAndSourcesCore.hpp
2 * \brief Generic finite element implementation using operator pipelines
3 *
4 * This file provides the core framework for finite element computations in MoFEM
5 * based on a pipeline of user data operators. Instead of monolithic finite element
6 * implementations, ForcesAndSourcesCore treats finite elements as sequences of
7 * modular operators that can be combined and reused.
8 *
9 * ## Core Concept: Operator Pipeline Architecture
10 *
11 * **Finite Element = Pipeline of Operators**
12 * Rather than traditional finite element classes with hardcoded functionality,
13 * MoFEM implements finite elements as pipelines of small, focused operators.
14 * Each operator performs a specific task (calculate gradients, apply boundary
15 * conditions, assemble matrices, etc.) and can be combined with others to
16 * create complex finite element behaviors.
17 *
18 * **ForcesAndSourcesCore**: Generic base class that manages the operator pipeline
19 * and provides the framework for element assembly. Users inherit from this class
20 * and add specific operators to implement their particular problem.
21 *
22 * ## Operator Types and Usage:
23 *
24 * **OPROW**: Operators that work with test functions (rows of system matrix)
25 * - Used for: Loading vectors, boundary conditions, source terms
26 * - Example: Applying Neumann boundary conditions, body forces
27 * - Operates during: Right-hand side vector assembly
28 *
29 * **OPCOL**: Operators that work with trial functions (columns of system matrix)
30 * - Used for: Geometric operations, coordinate transformations
31 * - Example: Calculating coordinates, setting up geometry
32 * - Operates during: Matrix structure setup and geometric computations
33 *
34 * **OPROWCOL**: Operators that work with both test and trial functions
35 * - Used for: Stiffness matrices, mass matrices, coupling terms
36 * - Example: Linear elasticity stiffness, thermal conductivity
37 * - Operates during: System matrix assembly (most common type)
38 *
39 * **OPSPACE**: Operators that work on the finite element space itself
40 * - Used for: Basis function transformations, coordinate mappings
41 * - Example: Piola transforms for Hdiv/Hcurl spaces, Jacobian operations
42 * - Operates during: Pre-processing before matrix/vector assembly
43 *
44 * ## Usage Pattern:
45 * 1. Inherit from appropriate ForcesAndSourcesCore variant
46 * 2. Add operators to pipeline using getOpPtrVector()
47 * 3. Operators execute automatically during element assembly
48 * 4. Combine different operator types to implement complete finite element
49 *
50 * ## Benefits:
51 * - **Modularity**: Reuse operators across different finite element types
52 * - **Flexibility**: Easy to modify behavior by adding/removing operators
53 * - **Maintainability**: Small, focused operators are easier to debug
54 * - **Extensibility**: Add new physics without modifying existing code
55 *
56 * @note This design enables rapid development of complex multiphysics
57 * finite element applications through operator composition.
58 */
59
60
61
62#ifndef __FORCES_AND_SOURCES_CORE__HPP__
63#define __FORCES_AND_SOURCES_CORE__HPP__
64
65using namespace boost::numeric;
66
67namespace MoFEM {
68
69/** \brief structure to get information from mofem into EntitiesFieldData
70 * \ingroup mofem_forces_and_sources
71 *
72 */
74
76
78 typedef boost::function<int(int order_row, int order_col, int order_data)>
80
81 typedef boost::function<MoFEMErrorCode(ForcesAndSourcesCore *fe_raw_ptr,
82 int order_row, int order_col,
83 int order_data)>
85
86 /**
87 * \brief Hook to get rule
88 *
89 * \todo check preferred format how works with gcc and clang,
90 * see
91 * <http://www.boost.org/doc/libs/1_64_0/doc/html/function/tutorial.html#idp247873024>
92 */
94
95 /**
96 * @brief Set function to calculate integration rule
97 *
98 */
100
101 /** \brief Data operator to do calculations at integration points.
102 * \ingroup mofem_forces_and_sources
103
104 Is inherited and implemented by user to do calculations. It can be used in
105 many different ways but typically is used to integrate matrices (f.e.
106 stiffness matrix) and the right hand vector (f.e. force vector).
107
108 Note: It is assumed that operator is executed for symmetric problem. That
109 means that is executed for not repeating entities on finite element. For
110 example on triangle we have nodes, 3 edges and one face. Because of symmetry
111 calculations are for: nodes-nodes, nodes-edge_0, nodes-edge_1, nodes-edge_2,
112 nodes-face,
113 edge_0-edge_0, edge_0-edge_1, edge_0-edge_2,
114 edge_1-edge_1, edge_1-edge_2,
115 edge_2-edge_2,
116 edge_0-face, edge_1-face, edge_2-face,
117 face - face
118
119 In case of non symmetric problem in addition calculations of lower off
120 diagonal terms. F.e. edge_1-edge_0, edge_2-edge_0, edge_2-edge_1,
121
122 In that case class variable UserDataOperator::sYmm = false;
123
124 Note: By default sYmm is set for symmetric problems
125
126 */
127 struct UserDataOperator;
128
129 /** \brief Use to push back operator for row operator
130
131 It can be used to calculate nodal forces or other quantities on the mesh.
132
133 */
134 boost::ptr_deque<UserDataOperator> &getOpPtrVector() { return opPtrVector; }
135
136 /**
137 * @brief Get the Entity Polynomial Base object
138 *
139 * @return boost::shared_ptr<BaseFunction>&&
140 */
142
143 /**
144 * @brief Get the User Polynomial Base object
145 *
146 * @return boost::shared_ptr<BaseFunction>&
147 */
149
150 /**
151 * @brief Matrix of integration points
152 *
153 * Columns is equal to number of integration points, numver of rows
154 * depends on dimension of finite element entity, for example for
155 * tetrahedron rows are x,y,z,weight. Last row is integration weight.
156 *
157 * FIXME: that should be moved to private class data and acessed only by
158 * member function
159 */
161
162 virtual MoFEMErrorCode preProcess();
163 virtual MoFEMErrorCode operator()();
164 virtual MoFEMErrorCode postProcess();
165
166public:
167 /** \brief Get max order of approximation for data fields
168
169 getMaxDataOrder() returns maximal order on entities, for
170 all data on the element. So for example if finite element is triangle, and
171 triangle base function have order 4 and on edges base function have order
172 2, this function returns 4.
173
174 If finite element has for example 2 or more approximated fields, for
175 example Pressure (order 3) and displacement field (order 5), this function
176 returns 5.
177
178 */
179 int getMaxDataOrder() const;
180
181 /// \brief Get max order of approximation for field in rows
182 int getMaxRowOrder() const;
183
184 /// \brief Get max order of approximation for field in columns
185 int getMaxColOrder() const;
186
187 /**
188 * @brief Get the entity data
189 *
190 * @param space
191 * @param type
192 * @param side
193 * @return EntitiesFieldData::EntData&
194 */
195 auto &getEntData(const FieldSpace space, const EntityType type,
196 const int side) {
197 return dataOnElement[space]->dataOnEntities[type][side];
198 }
199
200 /**
201 * @brief Get data on entities and space
202 *
203 * Entities data are stored by space, by entity type, and entity side.
204 *
205 * @return std::array<boost::shared_ptr<EntitiesFieldData>, LASTSPACE>
206 */
208
209 /**
210 * @brief Get derived data on entities and space
211 *
212 * Entities data are stored by space, by entity type, and entity side. Derived
213 * data is used to store data on columns, so it shares information about shape
214 * functions wih rows.
215 *
216 * @return std::array<boost::shared_ptr<EntitiesFieldData>, LASTSPACE>
217 */
219
220protected:
221 /**
222 * \brief get sense (orientation) of entity
223 * @param type type of entity
224 * @param data entity data
225 * @return error code
226 */
228 const EntityType type,
229 boost::ptr_vector<EntitiesFieldData::EntData> &data) const;
230
231 /**
232 * @brief Get the entity data order
233 *
234 * @param type
235 * @param space
236 * @param data
237 * @return MoFEMErrorCode
238 */
240 const EntityType type, const FieldSpace space,
241 boost::ptr_vector<EntitiesFieldData::EntData> &data) const;
242
243 /**
244 * @brief Get the entity sense (orientation)
245 *
246 * @tparam type
247 * @param data
248 * @return MoFEMErrorCode
249 */
250 template <EntityType type>
252 return getEntitySense(type, data.dataOnEntities[type]);
253 }
254
255 /**
256 * @brief Get the entity data order for given space
257 *
258 * @tparam type
259 * @param data
260 * @param space
261 * @return MoFEMErrorCode
262 */
263 template <EntityType type>
265 const FieldSpace space) const {
266 return getEntityDataOrder(type, space, data.dataOnEntities[type]);
267 }
268
269 /** \name Indices */
270
271 /**@{*/
272
273 /// \brief get node indices
274 template <typename EXTRACTOR>
276 getNodesIndices(const int bit_number,
277 FieldEntity_vector_view &ents_field, VectorInt &nodes_indices,
278 VectorInt &local_nodes_indices, EXTRACTOR &&extractor) const;
279
280 /// \brief get row node indices from FENumeredDofEntity_multiIndex
282 const int bit_number) const;
283
284 /// \brief get col node indices from FENumeredDofEntity_multiIndex
286 const int bit_number) const;
287
288 template <typename EXTRACTOR>
289 MoFEMErrorCode getEntityIndices(EntitiesFieldData &data, const int bit_number,
290 FieldEntity_vector_view &ents_field,
291 const EntityType type_lo,
292 const EntityType type_hi,
293 EXTRACTOR &&extractor) const;
294
296 getEntityRowIndices(EntitiesFieldData &data, const int bit_number,
297 const EntityType type_lo = MBVERTEX,
298 const EntityType type_hi = MBPOLYHEDRON) const;
299
301 getEntityColIndices(EntitiesFieldData &data, const int bit_number,
302 const EntityType type_lo = MBVERTEX,
303 const EntityType type_hi = MBPOLYHEDRON) const;
304
305
306 /// \brief get col NoField indices
308 const int bit_number) const;
309
310 /// \brief get col NoField indices
312 const int bit_number) const;
313
314 /**@}*/
315
316 /** \name Data */
317
318 /**@{*/
319
320 /** Get bit ref level in entities, and set it to data
321 */
323
324 /**
325 * \brief Get data on nodes
326 * @param data Data structure
327 * @param field_name Field name
328 * @return Error code
329 */
331 const int bit_number) const;
332
334 getEntityFieldData(EntitiesFieldData &data, const int bit_number,
335 const EntityType type_lo = MBVERTEX,
336 const EntityType type_hi = MBPOLYHEDRON) const;
337
338 /**
339 * @brief Get field data on entities where no field is defined
340 *
341 * This is not standards, since op_data, dataOnElement[NOFIELD], has to be set
342 * for rows and columns for each NOFIELD, also size of sides is not determined.
343 *
344 * @param data @param bit_number @return MoFEMErrorCode
345 */
346 template <typename EXTRACTOR>
348 const int bit_number,
349 EXTRACTOR &&extractor) const;
350
351 /**@}*/
352
353 /// \brief Get nodes on faces
355
356 /// \brief Get field approximation space and base on entities
359
360 /** \name Data form NumeredDofEntity_multiIndex */
361
362 /**@{*/
363
364 /// \brief get indices of nodal indices which are declared for problem but
365 /// not this particular element
367 const NumeredDofEntity_multiIndex &dofs,
368 VectorInt &nodes_indices) const;
369
370 /// \brief get indices by type (generic function) which are declared for
371 /// problem but not this particular element
373 const NumeredDofEntity_multiIndex &dofs,
374 EntityType type, int side_number,
375 VectorInt &indices) const;
376
378 VectorInt &nodes_indices) const;
380 EntityType type, int side_number,
381 VectorInt &indices) const;
383 VectorInt &nodes_indices) const;
385 EntityType type, int side_number,
386 VectorInt &indices) const;
387
388 /**@}*/
389
390 /**
391 * \brief another variant of getRule
392 * @param order_row order of base function on row
393 * @param order_col order of base function on columns
394 * @param order_data order of base function approximating data
395 * @return integration rule
396 *
397 * This function is overloaded by the user. The integration rule
398 * is set such that specific operator implemented by the user is
399 * integrated accurately. For example if user implement bilinear operator
400 * \f[
401 * b(u,v) =
402 * \int_\mathcal{T}
403 * \frac{\partial u_i}{\partial x_j}\frac{\partial v_i}{\partial x_j}
404 * \textrm{d}\mathcal{T}
405 * \f]
406 * then if \f$u\f$ and \f$v\f$ are polynomial of given \em order, then
407 * exact integral would be \code int getRule(int order) { return
408 * 2*(order-1); }; \endcode
409 *
410 * The integration points and weights are set appropriately for given
411 * entity type and integration rule from \ref quad.c
412 *
413 * Method \ref ForcesAndSourcesCore::getRule takes at argument takes
414 * maximal polynomial order set on the element on all fields defined on
415 * the element. If a user likes to have more control, another variant of
416 * this function can be called which distinguishing between field orders
417 * on rows, columns and data, the i.e. first argument of a bilinear form,
418 * the second argument of bilinear form and field coefficients on the
419 * element.
420 *
421 * \note If user set rule to -1 or any other negative integer, then method
422 * \ref ForcesAndSourcesCore::setGaussPts is called. In that method user
423 * can implement own (specific) integration method.
424 *
425 * \bug this function should be const
426 */
427 virtual int getRule(int order_row, int order_col, int order_data);
428
429 /** \brief set user specific integration rule
430
431 This function allows for user defined integration rule. The key is to
432 called matrix gaussPts, which is used by other MoFEM procedures. Matrix
433 has number of rows equal to problem dimension plus one, where last index
434 is used to store weight values. %Number of columns is equal to number of
435 integration points.
436
437 \note This function is called if method \ref
438 ForcesAndSourcesCore::getRule is returning integer -1 or any other
439 negative integer.
440
441 User sets
442 \code
443 MatrixDouble gaussPts;
444 \endcode
445 where
446 \code
447 gaussPts.resize(dim+1,nb_gauss_pts);
448 \endcode
449 number rows represents local coordinates of integration points
450 in reference element, where last index in row is for integration weight.
451
452 */
453 virtual MoFEMErrorCode setGaussPts(int order_row, int order_col,
454 int order_data);
455
456 /**
457 * \brief Calculate base functions
458 * @return Error code
459 */
462
463 /**
464 * \brief Calculate base functions
465 * @return Error code
466 */
468
469 /**
470 * @brief Calculate Bernstein-Bezier base
471 *
472 * @return MoFEMErrorCode
473 */
475
476 /**
477 * @brief Create a entity data on element object
478 *
479 * @return MoFEMErrorCode
480 */
481 MoFEMErrorCode createDataOnElement(EntityType type);
482
483 /**
484 * @brief Iterate user data operators
485 *
486 * @return MoFEMErrorCode
487 */
489
490 /**@{*/
491
492 /** \name Deprecated (do not use) */
493
494 /** \deprecated Use getRule(int row_order, int col_order, int data order)
495 */
496 virtual int getRule(int order);
497
498 /** \deprecated setGaussPts(int row_order, int col_order, int data order);
499 */
500 virtual MoFEMErrorCode setGaussPts(int order);
501
502 /**@}*/
503
504 /**
505 * @brief Entity data on element entity rows fields
506 *
507 */
508 const std::array<boost::shared_ptr<EntitiesFieldData>, LASTSPACE>
510
511 /**
512 * @brief Entity data on element entity columns fields
513 *
514 */
515 const std::array<boost::shared_ptr<EntitiesFieldData>, LASTSPACE>
517
523
524 /**
525 * @brief Vector of finite element users data operators
526 *
527 */
528 boost::ptr_deque<UserDataOperator> opPtrVector;
529
530 friend class UserDataOperator;
531
532protected:
533 /**
534 * @brief Last evaluated type of element entity
535 *
536 */
538
539private:
540 /**
541 * @brief Pointer to entity polynomial base
542 *
543 */
544 boost::shared_ptr<BaseFunction> elementPolynomialBasePtr;
545
546 /**
547 * @brief Pointer to user polynomial base
548 */
549 boost::shared_ptr<BaseFunction> userPolynomialBasePtr;
550
551 /**
552 * @brief Element to integrate on the sides
553 *
554 */
556
557 /**
558 * @brief Set the pointer to face element on the side
559 *
560 * \note Function is is used by face element, while it iterates over
561 * elements on the side
562 *
563 * @param side_fe_ptr
564 * @return MoFEMErrorCode
565 */
567
568 /**u
569 * @brief Element to integrate parent or child
570 *
571 */
573
574 /**
575 * @brief Set the pointer to face element refined
576 *
577 * @param refine_fe_ptr
578 * @return MoFEMErrorCode
579 */
581
587
588 template <int DIM> friend struct OpCopyGeomDataToE;
589 template <typename E> friend struct OpBrokenLoopSide;
590
591protected:
592 MatrixDouble coordsAtGaussPts; ///< coordinated at gauss points
593 double elementMeasure; ///< Depending on dimension of elements, stores length,
594 ///< area or volume of element.
595};
596
598
599 /**
600 * \brief Controls loop over entities on element
601 *
602 * - OPROW is used if row vector is assembled
603 * - OPCOL is usually used if column vector is assembled
604 * - OPROWCOL is usually used for assemble matrices.
605 * - OPSPACE no field is defined for such operator. Is usually used to modify
606 * base
607 *
608 *
609 * For typical problem like Bubnov-Galerkin OPROW and OPCOL are the same. In
610 * more general case for example for non-square matrices columns and rows
611 * could have different numeration and/or different set of fields.
612 *
613 */
614 enum OpType {
615 OPROW = 1 << 0, ///< operator doWork function is executed on FE rows
616 OPCOL = 1 << 1, ///< operator doWork function is executed on FE columns
617 OPROWCOL = 1 << 2, ///< operator doWork is executed on FE rows &columns
618 OPSPACE = 1 << 3, ///< operator do Work is execute on space data
619 OPLAST = 1 << 3 ///< @deprecated would be removed
620 };
621
622 static const char *const OpTypeNames[];
623
624 char opType;
625 std::string rowFieldName;
626 std::string colFieldName;
628
629 /**
630 * @brief Constructor for operators working on finite element spaces
631 *
632 * This constructor is used when operators modify basis functions on
633 * approximation spaces. The operator runs for all data on the specified
634 * space without access to specific field data. Commonly used for:
635 * - Coordinate transformations (Piola transforms)
636 * - Integration weight modifications
637 * - Basis function scaling or rotation
638 *
639 * @param space Field space to operate on (H1, HDIV, HCURL, L2, etc.)
640 * @param type Operator type (OPSPACE, OPROW, OPCOL, OPROWCOL)
641 * @param symm Whether the operator maintains matrix symmetry (default: true)
642 */
643 UserDataOperator(const FieldSpace space, const char type = OPSPACE,
644 const bool symm = true);
645
646 /**
647 * @brief Constructor for operators working on a single field
648 *
649 * This constructor creates operators that work with data from a specific
650 * field. Used for operations like:
651 * - Calculating field values at integration points
652 * - Computing field gradients or derivatives
653 * - Applying boundary conditions to specific fields
654 * - Post-processing field data
655 *
656 * @param field_name Name of the field to operate on
657 * @param type Operator type:
658 * - OPROW: Works with test functions (RHS assembly)
659 * - OPCOL: Works with trial functions (geometry setup)
660 * - OPROWCOL: Works with both (matrix assembly)
661 * @param symm Whether the operator maintains matrix symmetry (default: true)
662 */
663 UserDataOperator(const std::string field_name, const char type,
664 const bool symm = true);
665
666 /**
667 * @brief Constructor for operators working on two fields (bilinear forms)
668 *
669 * This constructor creates operators that couple two different fields,
670 * essential for:
671 * - Mixed formulations (pressure-velocity coupling)
672 * - Multi-physics problems (thermal-mechanical coupling)
673 * - Contact mechanics (displacement-Lagrange multiplier)
674 * - Interface problems between different field types
675 *
676 * @param row_field_name Name of the field for test functions (rows)
677 * @param col_field_name Name of the field for trial functions (columns)
678 * @param type Operator type (typically OPROWCOL for coupling)
679 * @param symm Whether the operator maintains matrix symmetry (default: true)
680 * Note: Coupling between different fields is often non-symmetric
681 */
682 UserDataOperator(const std::string row_field_name,
683 const std::string col_field_name, const char type,
684 const bool symm = true);
685
686 /** \brief Return raw pointer to NumeredEntFiniteElement
687 */
688 inline boost::shared_ptr<const NumeredEntFiniteElement>
690
691 /**
692 * \brief Return finite element entity handle
693 * @return Finite element entity handle
694 */
695 inline EntityHandle getFEEntityHandle() const;
696
697 /**
698 * @brief Get dimension of finite element
699 *
700 * @return int
701 */
702 inline int getFEDim() const;
703
704 /**
705 * @brief Get dimension of finite element
706 *
707 * @return int
708 */
709 inline EntityType getFEType() const;
710
711 /**
712 * @brief Get the side number pointer
713 *
714 * \note For vertex is expectation. Side basses in argument of function doWork
715 * is zero. For other entity types side can be used as argument of this
716 * function.
717 *
718 * @param side_number
719 * @param type
720 * @return boost::weak_ptr<SideNumber>
721 */
722 inline boost::weak_ptr<SideNumber> getSideNumberPtr(const int side_number,
723 const EntityType type);
724
725 /**
726 * @brief Get the side entity
727 *
728 * \note For vertex is expectation. Side basses in argument of function
729 * doWork is zero. For other entity types side can be used as argument of
730 * this function.
731 *
732 * \code
733 * MoFEMErrorCode doWork(int side, EntityType type,
734 * EntitiesFieldData::EntData &data) {
735 * MoFEMFunctionBegin;
736 *
737 * if (type == MBVERTEX) {
738 * for (int n = 0; n != number_of_nodes; ++n)
739 * EntityHandle ent = getSideEntity(n, type);
740 *
741 * // Do somthing
742 *
743 * } else {
744 * EntityHandle ent = getSideEntity(side, type);
745 *
746 * // Do somthing
747 *
748 * }
749 *
750 * MoFEMFunctionReturn(0);
751 * }
752 * \endcode
753 *
754 * @param side_number
755 * @param type
756 */
757 inline EntityHandle getSideEntity(const int side_number,
758 const EntityType type);
759
760 /**
761 * @brief Get the number of nodes on finite element
762 *
763 * @return int
764 */
765 inline int getNumberOfNodesOnElement() const;
766
767 /** \brief Get row indices
768
769 Field could be or not declared for this element but is declared for
770 problem
771
772 \param field_name
773 \param type entity type
774 \param side side number, any number if type is MBVERTEX
775 \return indices
776
777 NOTE: Using those indices to assemble matrix will result in error if new
778 non-zero values need to be created.
779
780 */
781 MoFEMErrorCode getProblemRowIndices(const std::string filed_name,
782 const EntityType type, const int side,
783 VectorInt &indices) const;
784
785 /** \brief Get col indices
786
787 Field could be or not declared for this element but is declared for
788 problem
789
790 \param field_name
791 \param type entity type
792 \param side side number, any number if type is MBVERTEX
793 \return indices
794
795 NOTE: Using those indices to assemble matrix will result in error if new
796 non-zero values need to be created.
797
798 */
799 MoFEMErrorCode getProblemColIndices(const std::string filed_name,
800 const EntityType type, const int side,
801 VectorInt &indices) const;
802
803 /** \brief Return raw pointer to Finite Element Method object
804 */
805 inline const FEMethod *getFEMethod() const;
806
807 /**
808 * \brief Get operator types
809 * @return Return operator type
810 */
811 inline int getOpType() const;
812
813 /**
814 * \brief Set operator type
815 * @param Operator type
816 */
817 inline void setOpType(const OpType type);
818
819 /**
820 * \brief Add operator type
821 */
822 inline void addOpType(const OpType type);
823
824 /**
825 * \brief get number of finite element in the loop
826 * @return number of finite element
827 */
828 inline int getNinTheLoop() const;
829
830 /**
831 * \brief get size of elements in the loop
832 * @return loop size
833 */
834 inline int getLoopSize() const;
835
836 /** \brief Get name of the element
837 */
838 inline std::string getFEName() const;
839
840 /** \name Accessing KSP */
841
842 /**@{*/
843
844 inline const PetscData::Switches &getDataCtx() const;
845
846 inline const KspMethod::KSPContext getKSPCtx() const;
847
848 inline const SnesMethod::SNESContext getSNESCtx() const;
849
850 inline const TSMethod::TSContext getTSCtx() const;
851
852 /**@}*/
853
854 /**@{*/
855
856 inline Vec getKSPf() const;
857
858 inline Mat getKSPA() const;
859
860 inline Mat getKSPB() const;
861
862 /**@}*/
863
864 /** \name Accessing SNES */
865
866 /**@{*/
867
868 inline Vec getSNESf() const;
869
870 inline Vec getSNESx() const;
871
872 inline Mat getSNESA() const;
873
874 inline Mat getSNESB() const;
875
876 /**@}*/
877
878 /** \name Accessing TS */
879
880 /**@{*/
881
882 inline Vec getTSu() const;
883
884 inline Vec getTSu_t() const;
885
886 inline Vec getTSu_tt() const;
887
888 inline Vec getTSf() const;
889
890 inline Mat getTSA() const;
891
892 inline Mat getTSB() const;
893
894 inline int getTSstep() const;
895
896 inline double getTStime() const;
897
898 inline double getTStimeStep() const;
899
900 inline double getTSa() const;
901
902 inline double getTSaa() const;
903
904 /**@}*/
905
906 /**@{*/
907
908 /**@}*/
909
910 /** \name Base functions and integration points */
911
912 /**@{*/
913
914 /** \brief matrix of integration (Gauss) points for Volume Element
915 *
916 * For triangle: columns 0,1 are x,y coordinates respectively and column
917 * 2 is a weight value for example getGaussPts()(1,13) returns y
918 * coordinate of 13th Gauss point on particular volume element
919 *
920 * For tetrahedron: columns 0,1,2 are x,y,z coordinates respectively and
921 * column 3 is a weight value for example getGaussPts()(1,13) returns y
922 * coordinate of 13th Gauss point on particular volume element
923 *
924 */
925 inline MatrixDouble &getGaussPts();
926
927 /**
928 * @brief Get integration weights
929 *
930 * \code
931 * auto t_w = getFTensor0IntegrationWeight();
932 * for(int gg = 0; gg!=getGaussPts.size2(); ++gg) {
933 * // integrate something
934 * ++t_w;
935 * }
936 * \endcode
937 *
938 * @return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>
939 */
940 inline auto getFTensor0IntegrationWeight();
941
942 /**@}*/
943
944 /** \name Coordinates and access to internal data */
945
946 /**@{*/
947
948 /** \brief Gauss points and weight, matrix (nb. of points x 3)
949
950 Column 0-2 integration points coordinate x, y and z, respectively. At rows
951 are integration points.
952
953 */
955
956 /**
957 * \brief Get coordinates at integration points assuming linear geometry
958 *
959 * \code
960 * auto t_coords = getFTensor1CoordsAtGaussPts();
961 * for(int gg = 0;gg!=nb_int_ptrs;gg++) {
962 * // do something
963 * ++t_coords;
964 * }
965 * \endcode
966 *
967 */
968 inline auto getFTensor1CoordsAtGaussPts();
969
970 /**@}*/
971
972 /**@{*/
973
974 /** \name Measures (area, volume, length, etc.) */
975
976 /**
977 * \brief get measure of element
978 * @return volume
979 */
980 inline double getMeasure() const;
981
982 /**
983 * \brief get measure of element
984 * @return volume
985 */
986 inline double &getMeasure();
987
988 /**}*/
989
990 /**@{*/
991
992 /** \name Loops */
993
994 using AdjCache =
995 std::map<EntityHandle,
996 std::vector<boost::weak_ptr<NumeredEntFiniteElement>>>;
997
998 /**
999 * @brief User calls this function to loop over elements on the side of
1000 * face. This function calls finite element with its operator to do
1001 * calculations.
1002 *
1003 * @param fe_name name of the side element
1004 * @param side_fe pointer to the side element instance
1005 * @param dim dimension the of side element
1006 * @param ent_for_side entity handle for which adjacent volume or face will
1007 * be accessed
1008 * @param verb
1009 * @param sev
1010 * @param adj_cache
1011 * @return MoFEMErrorCode
1012 */
1013 MoFEMErrorCode loopSide(const string &fe_name, ForcesAndSourcesCore *side_fe,
1014 const size_t dim, const EntityHandle ent_for_side = 0,
1015 boost::shared_ptr<Range> fe_range = nullptr,
1016 const int verb = QUIET,
1017 const LogManager::SeverityLevel sev = Sev::noisy,
1018 AdjCache *adj_cache = nullptr
1019
1020 );
1021
1022 /**
1023 * @brief User calls this function to loop over the same element using a
1024 * different set of integration points. This function calls finite element
1025 * with its operator to do calculations.
1026 *
1027 * @param fe_name
1028 * @param this_fe
1029 * @param verb
1030 * @param sev
1031 * @return MoFEMErrorCode
1032 */
1033 MoFEMErrorCode loopThis(const string &fe_name,
1034 ForcesAndSourcesCore *this_fe,
1035 const int verb = QUIET,
1036 const LogManager::SeverityLevel sev = Sev::noisy);
1037
1038 /**
1039 * @brief User calls this function to loop over parent elements. This function
1040 * calls finite element with its operator to do calculations.
1041 *
1042 * @param fe_name
1043 * @param parent_fe
1044 * @param verb
1045 * @param sev
1046 * @return MoFEMErrorCode
1047 */
1048 MoFEMErrorCode loopParent(const string &fe_name,
1049 ForcesAndSourcesCore *parent_fe,
1050 const int verb = QUIET,
1051 const LogManager::SeverityLevel sev = Sev::noisy);
1052
1053 /**
1054 * @brief User calls this function to loop over parent elements. This function
1055 * calls finite element with its operator to do calculations.
1056 *
1057 * @param fe_name
1058 * @param child_fe
1059 * @param verb
1060 * @param sev
1061 * @return MoFEMErrorCode
1062 */
1063 MoFEMErrorCode loopChildren(const string &fe_name,
1064 ForcesAndSourcesCore *child_fe,
1065 const int verb = QUIET,
1066 const LogManager::SeverityLevel sev = Sev::noisy);
1067
1068 /**
1069 * @brief Iterate over range of elements
1070 *
1071 * @param fe_name
1072 * @param range_fe
1073 * @param fe_range
1074 * @param verb
1075 * @param sev
1076 * @return * MoFEMErrorCode
1077 */
1078 MoFEMErrorCode loopRange(const string &fe_name,
1079 ForcesAndSourcesCore *range_fe,
1080 boost::shared_ptr<Range> fe_range,
1081 const int verb = QUIET,
1082 const LogManager::SeverityLevel sev = Sev::noisy);
1083
1084 /**@}*/
1085
1086 inline ForcesAndSourcesCore *getPtrFE() const;
1087
1088 inline ForcesAndSourcesCore *getSidePtrFE() const;
1089
1090 inline ForcesAndSourcesCore *getRefinePtrFE() const;
1091
1092
1093protected:
1095
1097
1098private:
1103};
1104
1105/// \deprecated Used ForcesAndSourcesCore instead
1107
1108boost::shared_ptr<const NumeredEntFiniteElement>
1109ForcesAndSourcesCore::UserDataOperator::getNumeredEntFiniteElementPtr() const {
1111};
1112
1113EntityHandle ForcesAndSourcesCore::UserDataOperator::getFEEntityHandle() const {
1114 return getNumeredEntFiniteElementPtr()->getEnt();
1115}
1116
1117int ForcesAndSourcesCore::UserDataOperator::getFEDim() const {
1119};
1120
1121EntityType ForcesAndSourcesCore::UserDataOperator::getFEType() const {
1123};
1124
1125boost::weak_ptr<SideNumber>
1126ForcesAndSourcesCore::UserDataOperator::getSideNumberPtr(
1127 const int side_number, const EntityType type) {
1128 auto &side_table_by_side_and_type =
1129 ptrFE->numeredEntFiniteElementPtr->getSideNumberTable().get<1>();
1130 auto side_it =
1131 side_table_by_side_and_type.find(boost::make_tuple(type, side_number));
1132 if (side_it != side_table_by_side_and_type.end())
1133 return *side_it;
1134 else
1135 return boost::weak_ptr<SideNumber>();
1136}
1137
1139ForcesAndSourcesCore::UserDataOperator::getSideEntity(const int side_number,
1140 const EntityType type) {
1141 if (auto side_ptr = getSideNumberPtr(side_number, type).lock())
1142 return side_ptr->ent;
1143 else
1144 return 0;
1145}
1146
1147int ForcesAndSourcesCore::UserDataOperator::getNumberOfNodesOnElement() const {
1148 return ptrFE->getNumberOfNodes();
1149}
1150
1151const FEMethod *ForcesAndSourcesCore::UserDataOperator::getFEMethod() const {
1152 return ptrFE;
1153}
1154
1155int ForcesAndSourcesCore::UserDataOperator::getOpType() const { return opType; }
1156
1157void ForcesAndSourcesCore::UserDataOperator::setOpType(const OpType type) {
1158 opType = type;
1159}
1160
1161void ForcesAndSourcesCore::UserDataOperator::addOpType(const OpType type) {
1162 opType |= type;
1163}
1164
1165int ForcesAndSourcesCore::UserDataOperator::getNinTheLoop() const {
1166 return getFEMethod()->getNinTheLoop();
1167}
1168
1169int ForcesAndSourcesCore::UserDataOperator::getLoopSize() const {
1170 return getFEMethod()->getLoopSize();
1171}
1172
1173std::string ForcesAndSourcesCore::UserDataOperator::getFEName() const {
1174 return getFEMethod()->getFEName();
1175}
1176
1177const PetscData::Switches &
1178ForcesAndSourcesCore::UserDataOperator::getDataCtx() const {
1179 return getFEMethod()->data_ctx;
1180}
1181
1183ForcesAndSourcesCore::UserDataOperator::getKSPCtx() const {
1184 return getFEMethod()->ksp_ctx;
1185}
1186
1188ForcesAndSourcesCore::UserDataOperator::getSNESCtx() const {
1189 return getFEMethod()->snes_ctx;
1190}
1191
1193ForcesAndSourcesCore::UserDataOperator::getTSCtx() const {
1194 return getFEMethod()->ts_ctx;
1195}
1196
1197Vec ForcesAndSourcesCore::UserDataOperator::getKSPf() const {
1198#ifndef NDEBUG
1199 if (getFEMethod()->ksp_f == PETSC_NULLPTR)
1200 THROW_MESSAGE("KSP not set F vector");
1201#endif
1202 return getFEMethod()->ksp_f;
1203}
1204
1205Mat ForcesAndSourcesCore::UserDataOperator::getKSPA() const {
1206#ifndef NDEBUG
1207 if (getFEMethod()->ksp_A == PETSC_NULLPTR)
1208 THROW_MESSAGE("KSP not set A vector");
1209#endif
1210 return getFEMethod()->ksp_A;
1211}
1212
1213Mat ForcesAndSourcesCore::UserDataOperator::getKSPB() const {
1214#ifndef NDEBUG
1215 if (getFEMethod()->ksp_B == PETSC_NULLPTR)
1216 THROW_MESSAGE("KSP not set B vector");
1217#endif
1218 return getFEMethod()->ksp_B;
1219}
1220
1221Vec ForcesAndSourcesCore::UserDataOperator::getSNESf() const {
1222#ifndef NDEBUG
1223 if (getFEMethod()->snes_f == PETSC_NULLPTR)
1224 THROW_MESSAGE("SNES not set F vector");
1225#endif
1226 return getFEMethod()->snes_f;
1227}
1228
1229Vec ForcesAndSourcesCore::UserDataOperator::getSNESx() const {
1230#ifndef NDEBUG
1231 if (getFEMethod()->snes_x == PETSC_NULLPTR)
1232 THROW_MESSAGE("SNESnot set X vector");
1233#endif
1234 return getFEMethod()->snes_x;
1235}
1236
1237Mat ForcesAndSourcesCore::UserDataOperator::getSNESA() const {
1238#ifndef NDEBUG
1239 if (getFEMethod()->snes_A == PETSC_NULLPTR)
1240 THROW_MESSAGE("SNES not set A vector");
1241#endif
1242 return getFEMethod()->snes_A;
1243}
1244
1245Mat ForcesAndSourcesCore::UserDataOperator::getSNESB() const {
1246#ifndef NDEBUG
1247 if (getFEMethod()->snes_B == PETSC_NULLPTR)
1248 THROW_MESSAGE("SNES not set A matrix");
1249#endif
1250 return getFEMethod()->snes_B;
1251}
1252
1253Vec ForcesAndSourcesCore::UserDataOperator::getTSu() const {
1254#ifndef NDEBUG
1255 if (getFEMethod()->ts_u == PETSC_NULLPTR)
1256 THROW_MESSAGE("TS not set U vector");
1257#endif
1258 return getFEMethod()->ts_u;
1259}
1260
1261Vec ForcesAndSourcesCore::UserDataOperator::getTSu_t() const {
1262#ifndef NDEBUG
1263 if (getFEMethod()->ts_u_t == PETSC_NULLPTR)
1264 THROW_MESSAGE("TS not set U_t vector");
1265#endif
1266 return getFEMethod()->ts_u_t;
1267}
1268
1269Vec ForcesAndSourcesCore::UserDataOperator::getTSu_tt() const {
1270#ifndef NDEBUG
1271 if (getFEMethod()->ts_u_tt == PETSC_NULLPTR)
1272 THROW_MESSAGE("TS not set U_tt vector");
1273#endif
1274 return getFEMethod()->ts_u_tt;
1275}
1276
1277Vec ForcesAndSourcesCore::UserDataOperator::getTSf() const {
1278#ifndef NDEBUG
1279 if (getFEMethod()->ts_F == PETSC_NULLPTR)
1280 THROW_MESSAGE("TS not set F vector");
1281#endif
1282 return getFEMethod()->ts_F;
1283}
1284
1285Mat ForcesAndSourcesCore::UserDataOperator::getTSA() const {
1286#ifndef NDEBUG
1287 if (getFEMethod()->ts_A == PETSC_NULLPTR)
1288 THROW_MESSAGE("TS not set A matrix");
1289#endif
1290 return getFEMethod()->ts_A;
1291}
1292
1293Mat ForcesAndSourcesCore::UserDataOperator::getTSB() const {
1294#ifndef NDEBUG
1295 if (getFEMethod()->ts_B == PETSC_NULLPTR)
1296 THROW_MESSAGE("TS not set B matrix");
1297#endif
1298 return getFEMethod()->ts_B;
1299}
1300
1301int ForcesAndSourcesCore::UserDataOperator::getTSstep() const {
1302#ifndef NDEBUG
1303 if ((getFEMethod()->data_ctx & PetscData::PetscData::CtxSetTime).none())
1304 THROW_MESSAGE("TS not set time");
1305#endif
1306 return getFEMethod()->ts_step;
1307}
1308
1309double ForcesAndSourcesCore::UserDataOperator::getTStime() const {
1310#ifndef NDEBUG
1311 if ((getFEMethod()->data_ctx & PetscData::PetscData::CtxSetTime).none())
1312 THROW_MESSAGE("TS not set time");
1313#endif
1314 return getFEMethod()->ts_t;
1315}
1316
1317double ForcesAndSourcesCore::UserDataOperator::getTStimeStep() const {
1318#ifndef NDEBUG
1319 if ((getFEMethod()->data_ctx & PetscData::PetscData::CtxSetTime).none())
1320 THROW_MESSAGE("TS not set time");
1321#endif
1322 return getFEMethod()->ts_dt;
1323}
1324
1325double ForcesAndSourcesCore::UserDataOperator::getTSa() const {
1326#ifndef NDEBUG
1327 if ((getFEMethod()->data_ctx & (PetscData::CtxSetA | PetscData::CtxSetB))
1328 .none() ||
1329 (getFEMethod()->data_ctx & (PetscData::CtxSetX_T)).none())
1330 THROW_MESSAGE("TS not set B matrix");
1331#endif
1332 return getFEMethod()->ts_a;
1333}
1334
1335double ForcesAndSourcesCore::UserDataOperator::getTSaa() const {
1336#ifndef NDEBUG
1337 if ((getFEMethod()->data_ctx & (PetscData::CtxSetA | PetscData::CtxSetB))
1338 .none() ||
1339 (getFEMethod()->data_ctx & (PetscData::CtxSetX_TT)).none())
1340 THROW_MESSAGE("TS not set B matrix");
1341#endif
1342 return getFEMethod()->ts_aa;
1343}
1344
1345MatrixDouble &ForcesAndSourcesCore::UserDataOperator::getGaussPts() {
1346 return static_cast<ForcesAndSourcesCore *>(ptrFE)->gaussPts;
1347}
1348
1349auto ForcesAndSourcesCore::UserDataOperator::getFTensor0IntegrationWeight() {
1351 &(getGaussPts()(getGaussPts().size1() - 1, 0)));
1352}
1353
1354// MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::getPorblemRowIndices(
1355// const std::string filed_name, const EntityType type, const int side,
1356// VectorInt &indices) const {
1357// return getProblemRowIndices(filed_name, type, side, indices);
1358// }
1359
1360ForcesAndSourcesCore *ForcesAndSourcesCore::UserDataOperator::getPtrFE() const {
1361 return ptrFE;
1362}
1363
1365ForcesAndSourcesCore::UserDataOperator::getSidePtrFE() const {
1366 return ptrFE->sidePtrFE;
1367}
1368
1370ForcesAndSourcesCore::UserDataOperator::getRefinePtrFE() const {
1371 return ptrFE->refinePtrFE;
1372}
1373
1374MatrixDouble &ForcesAndSourcesCore::UserDataOperator::getCoordsAtGaussPts() {
1375 return static_cast<ForcesAndSourcesCore *>(ptrFE)->coordsAtGaussPts;
1376}
1377
1378auto ForcesAndSourcesCore::UserDataOperator::getFTensor1CoordsAtGaussPts() {
1380 &getCoordsAtGaussPts()(0, 0), &getCoordsAtGaussPts()(0, 1),
1381 &getCoordsAtGaussPts()(0, 2));
1382}
1383
1384double ForcesAndSourcesCore::UserDataOperator::getMeasure() const {
1385 return static_cast<ForcesAndSourcesCore *>(ptrFE)->elementMeasure;
1386}
1387
1388double &ForcesAndSourcesCore::UserDataOperator::getMeasure() {
1389 return static_cast<ForcesAndSourcesCore *>(ptrFE)->elementMeasure;
1390}
1391
1392/**
1393 * @brief Element used to execute operators on side of the element
1394 *
1395 * @tparam E template for side element type
1396 *
1397 */
1398template <typename E>
1400
1402
1403 OpLoopSide(MoFEM::Interface &m_field, const std::string fe_name,
1404 const int side_dim, boost::shared_ptr<Range> fe_range,
1405 const LogManager::SeverityLevel sev = Sev::noisy,
1406 boost::shared_ptr<AdjCache> adj_cache = nullptr)
1407 : UserDataOperator(NOSPACE, OPSPACE), sideFEPtr(new E(m_field)),
1408 sideFEName(fe_name), sideDim(side_dim), sevLevel(sev),
1409 adjCache(adj_cache), feRange(fe_range) {}
1410
1411 /**
1412 * @brief Construct a new Op Loop Side object
1413 *
1414 * @param m_field
1415 * @param fe_name name of side (domain element)
1416 * @param side_dim dimension
1417 * @param sev severity level
1418 * @param adj_cache if set to null cache of the element is used
1419 */
1420 OpLoopSide(MoFEM::Interface &m_field, const std::string fe_name,
1421 const int side_dim,
1422 const LogManager::SeverityLevel sev = Sev::noisy,
1423 boost::shared_ptr<AdjCache> adj_cache = nullptr)
1424 : OpLoopSide(m_field, fe_name, side_dim, nullptr, sev, adj_cache) {}
1425
1426 MoFEMErrorCode doWork(int side, EntityType type,
1428 return loopSide(sideFEName, sideFEPtr.get(), sideDim, 0, feRange, VERBOSE,
1429 sevLevel, adjCache.get());
1430 };
1431
1432 boost::ptr_deque<UserDataOperator> &getOpPtrVector() {
1433 return sideFEPtr->getOpPtrVector();
1434 }
1435
1436 boost::shared_ptr<E> &getSideFEPtr() { return sideFEPtr; }
1437
1438protected:
1439 const std::string sideFEName;
1440 const int sideDim;
1441 boost::shared_ptr<E> sideFEPtr;
1443 boost::shared_ptr<AdjCache> adjCache;
1444 boost::shared_ptr<Range> feRange;
1445};
1446
1447/**
1448 * @brief Iterate over range of (sub)elements
1449 *
1450 * @tparam E
1451 */
1452template <typename E>
1454
1456
1457 OpLoopRange(MoFEM::Interface &m_field, const std::string fe_name,
1458 boost::shared_ptr<Range> fe_range,
1459 const LogManager::SeverityLevel sev = Sev::noisy)
1460 : UserDataOperator(NOSPACE, OPSPACE), rangeFEPtr(new E(m_field)),
1461 rangeFEName(fe_name), sevLevel(sev), feRange(fe_range) {}
1462
1463 MoFEMErrorCode doWork(int side, EntityType type,
1465 return loopRange(rangeFEName, rangeFEPtr.get(), feRange, VERBOSE, sevLevel);
1466 };
1467
1468 boost::ptr_deque<UserDataOperator> &getOpPtrVector() {
1469 return rangeFEPtr->getOpPtrVector();
1470 }
1471
1472 boost::shared_ptr<E> &getRangeFEPtr() { return rangeFEPtr; }
1473
1474protected:
1475 const std::string rangeFEName;
1476 boost::shared_ptr<E> rangeFEPtr;
1478 boost::shared_ptr<Range> feRange;
1479};
1480
1481/**
1482 * \brief Copy geometry-related data from one element to other
1483 *
1484 * That can be used to copy high order geometry data from coarse element to
1485 * children. That is often a case when higher order geometry is defined only on
1486 * coarse elements.
1487 *
1488 * \note Integration points have to be located at the same geometric positions
1489 *
1490 * FIXME: Write atom test
1491 */
1492template <int DIM>
1494
1495} // namespace MoFEM
1496
1497#endif //__FORCES_AND_SOURCES_CORE__HPP__
1498
1499/**
1500 * \defgroup mofem_forces_and_sources Forces and sources
1501 * \ingroup mofem
1502 *
1503 * \brief Manages complexities related to assembly of vector and matrices at
1504 * single finite element level.
1505 *
1506 **/
@ QUIET
@ VERBOSE
FieldApproximationBase
approximation base
Definition definitions.h:58
FieldSpace
approximation spaces
Definition definitions.h:82
@ LASTSPACE
FieldSpace in [ 0, LASTSPACE )
Definition definitions.h:89
@ NOSPACE
Definition definitions.h:83
#define DEPRECATED
Definition definitions.h:17
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
constexpr int order
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.
SeverityLevel
Severity levels.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasVector< int > VectorInt
Definition Types.hpp:67
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
auto type_from_handle(const EntityHandle h)
get type from entity handle
DEPRECATED typedef ForcesAndSourcesCore ForcesAndSurcesCore
auto dimension_from_handle(const EntityHandle h)
get entity dimension form handle
std::vector< boost::weak_ptr< FieldEntity > > FieldEntity_vector_view
constexpr auto field_name
base operator to do operations at Gauss Pt. level
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
data structure for finite element entity
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Structure for user loop methods on finite elements.
EntityHandle getFEEntityHandle() const
Get the entity handle of the current finite element.
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
Shared pointer to finite element database structure.
Base face element used to integrate on skeleton.
void setOpType(const OpType type)
Set operator type.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
MoFEMErrorCode getProblemColIndices(const std::string filed_name, const EntityType type, const int side, VectorInt &indices) const
Get col indices.
int getLoopSize() const
get size of elements in the loop
EntityType getFEType() const
Get dimension of finite element.
MoFEMErrorCode loopSide(const string &fe_name, ForcesAndSourcesCore *side_fe, const size_t dim, const EntityHandle ent_for_side=0, boost::shared_ptr< Range > fe_range=nullptr, const int verb=QUIET, const LogManager::SeverityLevel sev=Sev::noisy, AdjCache *adj_cache=nullptr)
User calls this function to loop over elements on the side of face. This function calls finite elemen...
MoFEMErrorCode loopRange(const string &fe_name, ForcesAndSourcesCore *range_fe, boost::shared_ptr< Range > fe_range, const int verb=QUIET, const LogManager::SeverityLevel sev=Sev::noisy)
Iterate over range of elements.
int getNinTheLoop() const
get number of finite element in the loop
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
int getNumberOfNodesOnElement() const
Get the number of nodes on finite element.
boost::weak_ptr< SideNumber > getSideNumberPtr(const int side_number, const EntityType type)
Get the side number pointer.
const SnesMethod::SNESContext getSNESCtx() const
MoFEMErrorCode loopParent(const string &fe_name, ForcesAndSourcesCore *parent_fe, const int verb=QUIET, const LogManager::SeverityLevel sev=Sev::noisy)
User calls this function to loop over parent elements. This function calls finite element with its op...
EntityHandle getSideEntity(const int side_number, const EntityType type)
Get the side entity.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
std::string getFEName() const
Get name of the element.
MoFEMErrorCode getProblemRowIndices(const std::string filed_name, const EntityType type, const int side, VectorInt &indices) const
Get row indices.
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
std::map< EntityHandle, std::vector< boost::weak_ptr< NumeredEntFiniteElement > > > AdjCache
void addOpType(const OpType type)
Add operator type.
OpType
Controls loop over entities on element.
@ OPCOL
operator doWork function is executed on FE columns
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
@ OPSPACE
operator do Work is execute on space data
MoFEMErrorCode loopChildren(const string &fe_name, ForcesAndSourcesCore *child_fe, const int verb=QUIET, const LogManager::SeverityLevel sev=Sev::noisy)
User calls this function to loop over parent elements. This function calls finite element with its op...
virtual MoFEMErrorCode setPtrFE(ForcesAndSourcesCore *ptr)
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
int getFEDim() const
Get dimension of finite element.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
MoFEMErrorCode loopThis(const string &fe_name, ForcesAndSourcesCore *this_fe, const int verb=QUIET, const LogManager::SeverityLevel sev=Sev::noisy)
User calls this function to loop over the same element using a different set of integration points....
structure to get information from mofem into EntitiesFieldData
MoFEMErrorCode loopOverOperators()
Iterate user data operators.
MoFEMErrorCode getProblemNodesRowIndices(const std::string &field_name, VectorInt &nodes_indices) const
auto & getDataOnElementBySpaceArray()
Get data on entities and space.
MoFEMErrorCode getEntityIndices(EntitiesFieldData &data, const int bit_number, FieldEntity_vector_view &ents_field, const EntityType type_lo, const EntityType type_hi, EXTRACTOR &&extractor) const
MoFEMErrorCode getEntityRowIndices(EntitiesFieldData &data, const int bit_number, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
MoFEMErrorCode getRowNodesIndices(EntitiesFieldData &data, const int bit_number) const
get row node indices from FENumeredDofEntity_multiIndex
int getMaxRowOrder() const
Get max order of approximation for field in rows.
auto & getUserPolynomialBase()
Get the User Polynomial Base object.
MoFEMErrorCode getEntitySense(EntitiesFieldData &data) const
Get the entity sense (orientation)
MoFEMErrorCode getSpacesAndBaseOnEntities(EntitiesFieldData &data) const
Get field approximation space and base on entities.
MoFEMErrorCode getProblemTypeColIndices(const std::string &field_name, EntityType type, int side_number, VectorInt &indices) const
virtual MoFEMErrorCode operator()()
Main operator function executed for each loop iteration.
MoFEMErrorCode getFaceNodes(EntitiesFieldData &data) const
Get nodes on faces.
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > derivedDataOnElement
Entity data on element entity columns fields.
MoFEMErrorCode getNoFieldRowIndices(EntitiesFieldData &data, const int bit_number) const
get col NoField indices
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
MoFEMErrorCode getNodesIndices(const int bit_number, FieldEntity_vector_view &ents_field, VectorInt &nodes_indices, VectorInt &local_nodes_indices, EXTRACTOR &&extractor) const
get node indices
MoFEMErrorCode getEntityFieldData(EntitiesFieldData &data, const int bit_number, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
MoFEMErrorCode setRefineFEPtr(const ForcesAndSourcesCore *refine_fe_ptr)
Set the pointer to face element refined.
auto & getElementPolynomialBase()
Get the Entity Polynomial Base object.
ForcesAndSourcesCore * sidePtrFE
Element to integrate on the sides.
MatrixDouble coordsAtGaussPts
coordinated at gauss points
boost::function< MoFEMErrorCode(ForcesAndSourcesCore *fe_raw_ptr, int order_row, int order_col, int order_data)> GaussHookFun
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
boost::ptr_deque< UserDataOperator > opPtrVector
Vector of finite element users data operators.
MoFEMErrorCode getEntitySense(const EntityType type, boost::ptr_vector< EntitiesFieldData::EntData > &data) const
get sense (orientation) of entity
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > dataOnElement
Entity data on element entity rows fields.
auto & getDerivedDataOnElementBySpaceArray()
Get derived data on entities and space.
MoFEMErrorCode getProblemNodesColIndices(const std::string &field_name, VectorInt &nodes_indices) const
MoFEMErrorCode getNoFieldEntityFieldData(EntitiesFieldData &data, const int bit_number, EXTRACTOR &&extractor) const
Get field data on entities where no field is defined.
MoFEMErrorCode calHierarchicalBaseFunctionsOnElement()
Calculate base functions.
MatrixDouble gaussPts
Matrix of integration points.
MoFEMErrorCode calBernsteinBezierBaseFunctionsOnElement()
Calculate Bernstein-Bezier base.
auto & getEntData(const FieldSpace space, const EntityType type, const int side)
Get the entity data.
MoFEMErrorCode getEntityDataOrder(EntitiesFieldData &data, const FieldSpace space) const
Get the entity data order for given space.
virtual MoFEMErrorCode preProcess()
Pre-processing function executed at loop initialization.
MoFEMErrorCode getColNodesIndices(EntitiesFieldData &data, const int bit_number) const
get col node indices from FENumeredDofEntity_multiIndex
MoFEMErrorCode setSideFEPtr(const ForcesAndSourcesCore *side_fe_ptr)
Set the pointer to face element on the side.
int getMaxColOrder() const
Get max order of approximation for field in columns.
EntityType lastEvaluatedElementEntityType
Last evaluated type of element entity.
boost::function< int(int order_row, int order_col, int order_data)> RuleHookFun
MoFEMErrorCode getProblemNodesIndices(const std::string &field_name, const NumeredDofEntity_multiIndex &dofs, VectorInt &nodes_indices) const
get indices of nodal indices which are declared for problem but not this particular element
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
GaussHookFun setRuleHook
Set function to calculate integration rule.
MoFEMErrorCode getNodesFieldData(EntitiesFieldData &data, const int bit_number) const
Get data on nodes.
MoFEMErrorCode getEntityColIndices(EntitiesFieldData &data, const int bit_number, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
boost::shared_ptr< BaseFunction > elementPolynomialBasePtr
Pointer to entity polynomial base.
MoFEMErrorCode getEntityDataOrder(const EntityType type, const FieldSpace space, boost::ptr_vector< EntitiesFieldData::EntData > &data) const
Get the entity data order.
ForcesAndSourcesCore * refinePtrFE
Element to integrate parent or child.
MoFEMErrorCode getNoFieldColIndices(EntitiesFieldData &data, const int bit_number) const
get col NoField indices
int getMaxDataOrder() const
Get max order of approximation for data fields.
boost::shared_ptr< BaseFunction > userPolynomialBasePtr
Pointer to user polynomial base.
MoFEMErrorCode getProblemTypeRowIndices(const std::string &field_name, EntityType type, int side_number, VectorInt &indices) const
MoFEMErrorCode getProblemTypeIndices(const std::string &field_name, const NumeredDofEntity_multiIndex &dofs, EntityType type, int side_number, VectorInt &indices) const
get indices by type (generic function) which are declared for problem but not this particular element
RuleHookFun getRuleHook
Hook to get rule.
MoFEMErrorCode createDataOnElement(EntityType type)
Create a entity data on element object.
virtual MoFEMErrorCode postProcess()
Post-processing function executed at loop completion.
Mat & ksp_B
Reference to preconditioner matrix in KSP context.
KSPContext
Context enumeration for KSP solver phases.
Mat & ksp_A
Reference to system matrix in KSP context.
Vec & ksp_f
Reference to residual vector in KSP context.
Operator for broken loop side.
Copy geometry-related data from one element to other.
Iterate over range of (sub)elements.
OpLoopRange(MoFEM::Interface &m_field, const std::string fe_name, boost::shared_ptr< Range > fe_range, const LogManager::SeverityLevel sev=Sev::noisy)
boost::shared_ptr< E > & getRangeFEPtr()
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
const std::string rangeFEName
const LogManager::SeverityLevel sevLevel
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< Range > feRange
boost::shared_ptr< E > rangeFEPtr
Element used to execute operators on side of the element.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< Range > feRange
boost::shared_ptr< E > & getSideFEPtr()
const std::string sideFEName
boost::shared_ptr< AdjCache > adjCache
boost::shared_ptr< E > sideFEPtr
OpLoopSide(MoFEM::Interface &m_field, const std::string fe_name, const int side_dim, boost::shared_ptr< Range > fe_range, const LogManager::SeverityLevel sev=Sev::noisy, boost::shared_ptr< AdjCache > adj_cache=nullptr)
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
const LogManager::SeverityLevel sevLevel
OpLoopSide(MoFEM::Interface &m_field, const std::string fe_name, const int side_dim, const LogManager::SeverityLevel sev=Sev::noisy, boost::shared_ptr< AdjCache > adj_cache=nullptr)
Construct a new Op Loop Side object.
static constexpr Switches CtxSetA
Jacobian matrix switch.
static constexpr Switches CtxSetX_TT
Second time derivative switch.
std::bitset< 8 > Switches
Bitset type for context switches.
Switches data_ctx
Current data context switches.
static constexpr Switches CtxSetX_T
First time derivative switch.
static constexpr Switches CtxSetB
Preconditioner matrix switch.
SNESContext
Context enumeration for SNES solver phases.
Vec & snes_f
Reference to residual vector.
Vec & snes_x
Reference to current solution state vector.
Mat & snes_B
Reference to preconditioner of Jacobian matrix.
Mat & snes_A
Reference to Jacobian matrix.
TSContext
Context enumeration for TS solver phases.
Mat & ts_A
Reference to Jacobian matrix: dF/dU + a*dF/dU_t + aa*dF/dU_tt.
Vec & ts_F
Reference to residual vector F(t,U,U_t,U_tt)
Vec & ts_u_tt
Reference to second time derivative of state vector d²U/dt²
Vec & ts_u_t
Reference to first time derivative of state vector dU/dt.
Vec & ts_u
Reference to current state vector U(t)
Mat & ts_B
Reference to preconditioner matrix for ts_A.
Base volume element used to integrate on contact surface (could be extended to other volume elements)...