v0.15.0
Loading...
Searching...
No Matches
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
10
11#ifndef __FORCES_AND_SOURCES_CORE__HPP__
12#define __FORCES_AND_SOURCES_CORE__HPP__
13
14using namespace boost::numeric;
15
16namespace MoFEM {
17
18/** \brief structure to get information form mofem into EntitiesFieldData
19 * \ingroup mofem_forces_and_sources
20 *
21 */
23
25
27 typedef boost::function<int(int order_row, int order_col, int order_data)>
29
30 typedef boost::function<MoFEMErrorCode(ForcesAndSourcesCore *fe_raw_ptr,
31 int order_row, int order_col,
32 int order_data)>
34
35 /**
36 * \brief Hook to get rule
37 *
38 * \todo check preferred format how works with gcc and clang,
39 * see
40 * <http://www.boost.org/doc/libs/1_64_0/doc/html/function/tutorial.html#idp247873024>
41 */
43
44 /**
45 * @brief Set function to calculate integration rule
46 *
47 */
49
50 /** \brief Data operator to do calculations at integration points.
51 * \ingroup mofem_forces_and_sources
52
53 Is inherited and implemented by user to do calculations. It can be used in
54 many different ways but typically is used to integrate matrices (f.e.
55 stiffness matrix) and the right hand vector (f.e. force vector).
56
57 Note: It is assumed that operator is executed for symmetric problem. That
58 means that is executed for not repeating entities on finite element. For
59 example on triangle we have nodes, 3 edges and one face. Because of symmetry
60 calculations are for: nodes-nodes, nodes-edge0, nodes-edge_1, nodes-edge_2,
61 nodes-face,
62 edges_1-edges_1, edges_1-edge_1, edge_1-edge_2,
63 edge_1-edge_1, edge_1-edge_2,
64 edge_2-edge_2,
65 edge_1-face, edge_1-face, edges_3-face,
66 face - face
67
68 In case of non symmetric problem in addition calculations of lower off
69 diagonal terms. F.e. edge_1-edge_0, esges_3-edge_0, edge_3-edge_1,
70
71 In that case class variable UserDataOperator::sYmm = false;
72
73 NoteL: By default sYmm is set for symmetric problems
74
75 */
76 struct UserDataOperator;
77
78 /** \brief Use to push back operator for row operator
79
80 It can be used to calculate nodal forces or other quantities on the mesh.
81
82 */
83 boost::ptr_deque<UserDataOperator> &getOpPtrVector() { return opPtrVector; }
84
85 /**
86 * @brief Get the Entity Polynomial Base object
87 *
88 * @return boost::shared_ptr<BaseFunction>&&
89 */
91
92 /**
93 * @brief Get the User Polynomial Base object
94 *
95 * @return boost::shared_ptr<BaseFunction>&
96 */
98
99 /**
100 * @brief Matrix of integration points
101 *
102 * Columns is equal to number of integration points, numver of rows
103 * depends on dimension of finite element entity, for example for
104 * tetrahedron rows are x,y,z,weight. Last row is integration weight.
105 *
106 * FIXME: that should be moved to private class data and acessed only by
107 * member function
108 */
110
111 virtual MoFEMErrorCode preProcess();
112 virtual MoFEMErrorCode operator()();
113 virtual MoFEMErrorCode postProcess();
114
115public:
116 /** \brief Get max order of approximation for data fields
117
118 getMeasure getMaxDataOrder () return maximal order on entities, for
119 all data on the element. So for example if finite element is triangle, and
120 triangle base function have order 4 and on edges base function have order
121 2, this function return 4.
122
123 If finite element has for example 2 or more approximated fields, for
124 example Pressure (order 3) and displacement field (order 5), this function
125 returns 5.
126
127 */
128 int getMaxDataOrder() const;
129
130 /// \brief Get max order of approximation for field in rows
131 int getMaxRowOrder() const;
132
133 /// \brief Get max order of approximation for field in columns
134 int getMaxColOrder() const;
135
136 /**
137 * @brief Get the entity data
138 *
139 * @param space
140 * @param type
141 * @param side
142 * @return EntitiesFieldData::EntData&
143 */
144 auto &getEntData(const FieldSpace space, const EntityType type,
145 const int side) {
146 return dataOnElement[space]->dataOnEntities[type][side];
147 }
148
149 /**
150 * @brief Get data on entities and space
151 *
152 * Entities data are stored by space, by entity type, and entity side.
153 *
154 * @return std::array<boost::shared_ptr<EntitiesFieldData>, LASTSPACE>
155 */
157
158 /**
159 * @brief Get derived data on entities and space
160 *
161 * Entities data are stored by space, by entity type, and entity side. Derived
162 * data is used to store data on columns, so it shares information about shape
163 * functions wih rows.
164 *
165 * @return std::array<boost::shared_ptr<EntitiesFieldData>, LASTSPACE>
166 */
168
169protected:
170 /**
171 * \brief get sense (orientation) of entity
172 * @param type type of entity
173 * @param data entity data
174 * @return error code
175 */
177 const EntityType type,
178 boost::ptr_vector<EntitiesFieldData::EntData> &data) const;
179
180 /**
181 * @brief Get the entity data order
182 *
183 * @param type
184 * @param space
185 * @param data
186 * @return MoFEMErrorCode
187 */
189 const EntityType type, const FieldSpace space,
190 boost::ptr_vector<EntitiesFieldData::EntData> &data) const;
191
192 /**
193 * @brief Get the entity sense (orientation)
194 *
195 * @tparam type
196 * @param data
197 * @return MoFEMErrorCode
198 */
199 template <EntityType type>
201 return getEntitySense(type, data.dataOnEntities[type]);
202 }
203
204 /**
205 * @brief Get the entity data order for given space
206 *
207 * @tparam type
208 * @param data
209 * @param space
210 * @return MoFEMErrorCode
211 */
212 template <EntityType type>
214 const FieldSpace space) const {
215 return getEntityDataOrder(type, space, data.dataOnEntities[type]);
216 }
217
218 /** \name Indices */
219
220 /**@{*/
221
222 /// \brief get node indices
223 template <typename EXTRACTOR>
225 getNodesIndices(const int bit_number,
226 FieldEntity_vector_view &ents_field, VectorInt &nodes_indices,
227 VectorInt &local_nodes_indices, EXTRACTOR &&extractor) const;
228
229 /// \brief get row node indices from FENumeredDofEntity_multiIndex
231 const int bit_number) const;
232
233 /// \brief get col node indices from FENumeredDofEntity_multiIndex
235 const int bit_number) const;
236
237 template <typename EXTRACTOR>
238 MoFEMErrorCode getEntityIndices(EntitiesFieldData &data, const int bit_number,
239 FieldEntity_vector_view &ents_field,
240 const EntityType type_lo,
241 const EntityType type_hi,
242 EXTRACTOR &&extractor) const;
243
245 getEntityRowIndices(EntitiesFieldData &data, const int bit_number,
246 const EntityType type_lo = MBVERTEX,
247 const EntityType type_hi = MBPOLYHEDRON) const;
248
250 getEntityColIndices(EntitiesFieldData &data, const int bit_number,
251 const EntityType type_lo = MBVERTEX,
252 const EntityType type_hi = MBPOLYHEDRON) const;
253
254
255 /// \brief get col NoField indices
257 const int bit_number) const;
258
259 /// \brief get col NoField indices
261 const int bit_number) const;
262
263 /**@}*/
264
265 /** \name Data */
266
267 /**@{*/
268
269 /** Get bit ref level in entities, and set it to data
270 */
272
273 /**
274 * \brief Get data on nodes
275 * @param data Data structure
276 * @param field_name Field name
277 * @return Error code
278 */
280 const int bit_number) const;
281
283 getEntityFieldData(EntitiesFieldData &data, const int bit_number,
284 const EntityType type_lo = MBVERTEX,
285 const EntityType type_hi = MBPOLYHEDRON) const;
286
287 /**
288 * @brief Get field data on entities where no field is defined
289 *
290 * This is not standards, since op_data, dataOnElement[NOFIELD], has to be set
291 * for rows and columns for each NOFIELD, also size of sides is not determined.
292 *
293 * @param data @param bit_number @return MoFEMErrorCode
294 */
295 template <typename EXTRACTOR>
297 const int bit_number,
298 EXTRACTOR &&extractor) const;
299
300 /**@}*/
301
302 /// \brief Get nodes on faces
304
305 /// \brief Get field approximation space and base on entities
308
309 /** \name Data form NumeredDofEntity_multiIndex */
310
311 /**@{*/
312
313 /// \brief get indices of nodal indices which are declared for problem but
314 /// not this particular element
316 const NumeredDofEntity_multiIndex &dofs,
317 VectorInt &nodes_indices) const;
318
319 /// \brief get indices by type (generic function) which are declared for
320 /// problem but not this particular element
322 const NumeredDofEntity_multiIndex &dofs,
323 EntityType type, int side_number,
324 VectorInt &indices) const;
325
327 VectorInt &nodes_indices) const;
329 EntityType type, int side_number,
330 VectorInt &indices) const;
332 VectorInt &nodes_indices) const;
334 EntityType type, int side_number,
335 VectorInt &indices) const;
336
337 /**@}*/
338
339 /**
340 * \brief another variant of getRule
341 * @param order_row order of base function on row
342 * @param order_col order of base function on columns
343 * @param order_data order of base function approximating data
344 * @return integration rule
345 *
346 * This function is overloaded by the user. The integration rule
347 * is set such that specific operator implemented by the user is
348 * integrated accurately. For example if user implement bilinear operator
349 * \f[
350 * b(u,v) =
351 * \int_\mathcal{T}
352 * \frac{\partial u_i}{\partial x_j}\frac{\partial v_i}{\partial x_j}
353 * \textrm{d}\mathcal{T}
354 * \f]
355 * then if \f$u\f$ and \f$v\f$ are polynomial of given \em order, then
356 * exact integral would be \code int getRule(int order) { return
357 * 2*(order-1); }; \endcode
358 *
359 * The integration points and weights are set appropriately for given
360 * entity type and integration rule from \ref quad.c
361 *
362 * Method \ref ForcesAndSourcesCore::getRule takes at argument takes
363 * maximal polynomial order set on the element on all fields defined on
364 * the element. If a user likes to have more control, another variant of
365 * this function can be called which distinguishing between field orders
366 * on rows, columns and data, the i.e. first argument of a bilinear form,
367 * the second argument of bilinear form and field coefficients on the
368 * element.
369 *
370 * \note If user set rule to -1 or any other negative integer, then method
371 * \ref ForcesAndSourcesCore::setGaussPts is called. In that method user
372 * can implement own (specific) integration method.
373 *
374 * \bug this function should be const
375 */
376 virtual int getRule(int order_row, int order_col, int order_data);
377
378 /** \brief set user specific integration rule
379
380 This function allows for user defined integration rule. The key is to
381 called matrix gaussPts, which is used by other MoFEM procedures. Matrix
382 has number of rows equal to problem dimension plus one, where last index
383 is used to store weight values. %Number of columns is equal to number of
384 integration points.
385
386 \note This function is called if method \ref
387 ForcesAndSourcesCore::getRule is returning integer -1 or any other
388 negative integer.
389
390 User sets
391 \code
392 MatrixDouble gaussPts;
393 \endcode
394 where
395 \code
396 gaussPts.resize(dim+1,nb_gauss_pts);
397 \endcode
398 number rows represents local coordinates of integration points
399 in reference element, where last index in row is for integration weight.
400
401 */
402 virtual MoFEMErrorCode setGaussPts(int order_row, int order_col,
403 int order_data);
404
405 /**
406 * \brief Calculate base functions
407 * @return Error code
408 */
411
412 /**
413 * \brief Calculate base functions
414 * @return Error code
415 */
417
418 /**
419 * @brief Calculate Bernstein-Bezier base
420 *
421 * @return MoFEMErrorCode
422 */
424
425 /**
426 * @brief Create a entity data on element object
427 *
428 * @return MoFEMErrorCode
429 */
430 MoFEMErrorCode createDataOnElement(EntityType type);
431
432 /**
433 * @brief Iterate user data operators
434 *
435 * @return MoFEMErrorCode
436 */
438
439 /**@{*/
440
441 /** \name Deprecated (do not use) */
442
443 /** \deprecated Use getRule(int row_order, int col_order, int data order)
444 */
445 virtual int getRule(int order);
446
447 /** \deprecated setGaussPts(int row_order, int col_order, int data order);
448 */
449 virtual MoFEMErrorCode setGaussPts(int order);
450
451 /**@}*/
452
453 /**
454 * @brief Entity data on element entity rows fields
455 *
456 */
457 const std::array<boost::shared_ptr<EntitiesFieldData>, LASTSPACE>
459
460 /**
461 * @brief Entity data on element entity columns fields
462 *
463 */
464 const std::array<boost::shared_ptr<EntitiesFieldData>, LASTSPACE>
466
472
473 /**
474 * @brief Vector of finite element users data operators
475 *
476 */
477 boost::ptr_deque<UserDataOperator> opPtrVector;
478
479 friend class UserDataOperator;
480
481protected:
482 /**
483 * @brief Last evaluated type of element entity
484 *
485 */
487
488private:
489 /**
490 * @brief Pointer to entity polynomial base
491 *
492 */
493 boost::shared_ptr<BaseFunction> elementPolynomialBasePtr;
494
495 /**
496 * @brief Pointer to user polynomial base
497 */
498 boost::shared_ptr<BaseFunction> userPolynomialBasePtr;
499
500 /**
501 * @brief Element to integrate on the sides
502 *
503 */
505
506 /**
507 * @brief Set the pointer to face element on the side
508 *
509 * \note Function is is used by face element, while it iterates over
510 * elements on the side
511 *
512 * @param side_fe_ptr
513 * @return MoFEMErrorCode
514 */
516
517 /**u
518 * @brief Element to integrate parent or child
519 *
520 */
522
523 /**
524 * @brief Set the pointer to face element refined
525 *
526 * @param refine_fe_ptr
527 * @return MoFEMErrorCode
528 */
530
536
537 template <int DIM> friend struct OpCopyGeomDataToE;
538 template <typename E> friend struct OpBrokenLoopSide;
539
540protected:
541 MatrixDouble coordsAtGaussPts; ///< coordinated at gauss points
542 double elementMeasure; ///< Depending on dimension of elements, stores length,
543 ///< area or volume of element.
544};
545
547
548 /**
549 * \brief Controls loop over entities on element
550 *
551 * - OPROW is used if row vector is assembled
552 * - OPCOL is usually used if column vector is assembled
553 * - OPROWCOL is usually used for assemble matrices.
554 * - OPSPACE no field is defined for such operator. Is usually used to modify
555 * base
556 *
557 *
558 * For typical problem like Bubnov-Galerkin OPROW and OPCOL are the same. In
559 * more general case for example for non-square matrices columns and rows
560 * could have different numeration and/or different set of fields.
561 *
562 */
563 enum OpType {
564 OPROW = 1 << 0, ///< operator doWork function is executed on FE rows
565 OPCOL = 1 << 1, ///< operator doWork function is executed on FE columns
566 OPROWCOL = 1 << 2, ///< operator doWork is executed on FE rows &columns
567 OPSPACE = 1 << 3, ///< operator do Work is execute on space data
568 OPLAST = 1 << 3 ///< @deprecated would be removed
569 };
570
571 static const char *const OpTypeNames[];
572
573 char opType;
574 std::string rowFieldName;
575 std::string colFieldName;
577
578 /**
579 * This Constructor is used typically when some modification base shape
580 * functions on some approx. space is applied. Operator is run for all data
581 * on space.
582 *
583 * User has no access to field data from this operator.
584 */
585 UserDataOperator(const FieldSpace space, const char type = OPSPACE,
586 const bool symm = true);
587
588 UserDataOperator(const std::string field_name, const char type,
589 const bool symm = true);
590
591 UserDataOperator(const std::string row_field_name,
592 const std::string col_field_name, const char type,
593 const bool symm = true);
594
595 /** \brief Return raw pointer to NumeredEntFiniteElement
596 */
597 inline boost::shared_ptr<const NumeredEntFiniteElement>
599
600 /**
601 * \brief Return finite element entity handle
602 * @return Finite element entity handle
603 */
604 inline EntityHandle getFEEntityHandle() const;
605
606 /**
607 * @brief Get dimension of finite element
608 *
609 * @return int
610 */
611 inline int getFEDim() const;
612
613 /**
614 * @brief Get dimension of finite element
615 *
616 * @return int
617 */
618 inline EntityType getFEType() const;
619
620 /**
621 * @brief Get the side number pointer
622 *
623 * \note For vertex is expectation. Side basses in argument of function doWork
624 * is zero. For other entity types side can be used as argument of this
625 * function.
626 *
627 * @param side_number
628 * @param type
629 * @return boost::weak_ptr<SideNumber>
630 */
631 inline boost::weak_ptr<SideNumber> getSideNumberPtr(const int side_number,
632 const EntityType type);
633
634 /**
635 * @brief Get the side entity
636 *
637 * \note For vertex is expectation. Side basses in argument of function
638 * doWork is zero. For other entity types side can be used as argument of
639 * this function.
640 *
641 * \code
642 * MoFEMErrorCode doWork(int side, EntityType type,
643 * EntitiesFieldData::EntData &data) {
644 * MoFEMFunctionBegin;
645 *
646 * if (type == MBVERTEX) {
647 * for (int n = 0; n != number_of_nodes; ++n)
648 * EntityHandle ent = getSideEntity(n, type);
649 *
650 * // Do somthing
651 *
652 * } else {
653 * EntityHandle ent = getSideEntity(side, type);
654 *
655 * // Do somthing
656 *
657 * }
658 *
659 * MoFEMFunctionReturn(0);
660 * }
661 * \endcode
662 *
663 * @param side_number
664 * @param type
665 */
666 inline EntityHandle getSideEntity(const int side_number,
667 const EntityType type);
668
669 /**
670 * @brief Get the number of nodes on finite element
671 *
672 * @return int
673 */
674 inline int getNumberOfNodesOnElement() const;
675
676 /** \brief Get row indices
677
678 Field could be or not declared for this element but is declared for
679 problem
680
681 \param field_name
682 \param type entity type
683 \param side side number, any number if type is MBVERTEX
684 \return indices
685
686 NOTE: Using those indices to assemble matrix will result in error if new
687 non-zero values need to be created.
688
689 */
690 MoFEMErrorCode getProblemRowIndices(const std::string filed_name,
691 const EntityType type, const int side,
692 VectorInt &indices) const;
693
694 /** \brief Get col indices
695
696 Field could be or not declared for this element but is declared for
697 problem
698
699 \param field_name
700 \param type entity type
701 \param side side number, any number if type is MBVERTEX
702 \return indices
703
704 NOTE: Using those indices to assemble matrix will result in error if new
705 non-zero values need to be created.
706
707 */
708 MoFEMErrorCode getProblemColIndices(const std::string filed_name,
709 const EntityType type, const int side,
710 VectorInt &indices) const;
711
712 /** \brief Return raw pointer to Finite Element Method object
713 */
714 inline const FEMethod *getFEMethod() const;
715
716 /**
717 * \brief Get operator types
718 * @return Return operator type
719 */
720 inline int getOpType() const;
721
722 /**
723 * \brief Set operator type
724 * @param Operator type
725 */
726 inline void setOpType(const OpType type);
727
728 /**
729 * \brief Add operator type
730 */
731 inline void addOpType(const OpType type);
732
733 /**
734 * \brief get number of finite element in the loop
735 * @return number of finite element
736 */
737 inline int getNinTheLoop() const;
738
739 /**
740 * \brief get size of elements in the loop
741 * @return loop size
742 */
743 inline int getLoopSize() const;
744
745 /** \brief Get name of the element
746 */
747 inline std::string getFEName() const;
748
749 /** \name Accessing KSP */
750
751 /**@{*/
752
753 inline const PetscData::Switches &getDataCtx() const;
754
755 inline const KspMethod::KSPContext getKSPCtx() const;
756
757 inline const SnesMethod::SNESContext getSNESCtx() const;
758
759 inline const TSMethod::TSContext getTSCtx() const;
760
761 /**@}*/
762
763 /**@{*/
764
765 inline Vec getKSPf() const;
766
767 inline Mat getKSPA() const;
768
769 inline Mat getKSPB() const;
770
771 /**@}*/
772
773 /** \name Accessing SNES */
774
775 /**@{*/
776
777 inline Vec getSNESf() const;
778
779 inline Vec getSNESx() const;
780
781 inline Mat getSNESA() const;
782
783 inline Mat getSNESB() const;
784
785 /**@}*/
786
787 /** \name Accessing TS */
788
789 /**@{*/
790
791 inline Vec getTSu() const;
792
793 inline Vec getTSu_t() const;
794
795 inline Vec getTSu_tt() const;
796
797 inline Vec getTSf() const;
798
799 inline Mat getTSA() const;
800
801 inline Mat getTSB() const;
802
803 inline int getTSstep() const;
804
805 inline double getTStime() const;
806
807 inline double getTStimeStep() const;
808
809 inline double getTSa() const;
810
811 inline double getTSaa() const;
812
813 /**@}*/
814
815 /**@{*/
816
817 /** \name Base funtions and integration points */
818
819 /** \brief matrix of integration (Gauss) points for Volume Element
820 *
821 * For triangle: columns 0,1 are x,y coordinates respectively and column
822 * 2 is a weight value for example getGaussPts()(1,13) returns y
823 * coordinate of 13th Gauss point on particular volume element
824 *
825 * For tetrahedron: columns 0,1,2 are x,y,z coordinates respectively and
826 * column 3 is a weight value for example getGaussPts()(1,13) returns y
827 * coordinate of 13th Gauss point on particular volume element
828 *
829 */
830 inline MatrixDouble &getGaussPts();
831
832 /**
833 * @brief Get integration weights
834 *
835 * \code
836 * auto t_w = getFTensor0IntegrationWeight();
837 * for(int gg = 0; gg!=getGaussPts.size2(); ++gg) {
838 * // integrate something
839 * ++t_w;
840 * }
841 * \endcode
842 *
843 * @return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>
844 */
845 inline auto getFTensor0IntegrationWeight();
846
847 /**@}*/
848
849 /** \name Coordinates and access to internal data */
850
851 /**@{*/
852
853 /** \brief Gauss points and weight, matrix (nb. of points x 3)
854
855 Column 0-2 integration points coordinate x, y and z, respectively. At rows
856 are integration points.
857
858 */
860
861 /**
862 * \brief Get coordinates at integration points assuming linear geometry
863 *
864 * \code
865 * auto t_coords = getFTensor1CoordsAtGaussPts();
866 * for(int gg = 0;gg!=nb_int_ptrs;gg++) {
867 * // do something
868 * ++t_coords;
869 * }
870 * \endcode
871 *
872 */
873 inline auto getFTensor1CoordsAtGaussPts();
874
875 /**@}*/
876
877 /**@{*/
878
879 /** \name Measures (area, volume, length, etc.) */
880
881 /**
882 * \brief get measure of element
883 * @return volume
884 */
885 inline double getMeasure() const;
886
887 /**
888 * \brief get measure of element
889 * @return volume
890 */
891 inline double &getMeasure();
892
893 /**}*/
894
895 /**@{*/
896
897 /** \name Loops */
898
899 using AdjCache =
900 std::map<EntityHandle,
901 std::vector<boost::weak_ptr<NumeredEntFiniteElement>>>;
902
903 /**
904 * @brief User calls this function to loop over elements on the side of
905 * face. This function calls finite element with its operator to do
906 * calculations.
907 *
908 * @param fe_name name of the side element
909 * @param side_fe pointer to the side element instance
910 * @param dim dimension the of side element
911 * @param ent_for_side entity handle for which adjacent volume or face will
912 * be accessed
913 * @param verb
914 * @param sev
915 * @param adj_cache
916 * @return MoFEMErrorCode
917 */
918 MoFEMErrorCode loopSide(const string &fe_name, ForcesAndSourcesCore *side_fe,
919 const size_t dim, const EntityHandle ent_for_side = 0,
920 boost::shared_ptr<Range> fe_range = nullptr,
921 const int verb = QUIET,
922 const LogManager::SeverityLevel sev = Sev::noisy,
923 AdjCache *adj_cache = nullptr
924
925 );
926
927 /**
928 * @brief User calls this function to loop over the same element using a
929 * different set of integration points. This function calls finite element
930 * with its operator to do calculations.
931 *
932 * @param fe_name
933 * @param this_fe
934 * @param verb
935 * @param sev
936 * @return MoFEMErrorCode
937 */
938 MoFEMErrorCode loopThis(const string &fe_name,
939 ForcesAndSourcesCore *this_fe,
940 const int verb = QUIET,
941 const LogManager::SeverityLevel sev = Sev::noisy);
942
943 /**
944 * @brief User calls this function to loop over parent elements. This function
945 * calls finite element with its operator to do calculations.
946 *
947 * @param fe_name
948 * @param parent_fe
949 * @param verb
950 * @param sev
951 * @return MoFEMErrorCode
952 */
953 MoFEMErrorCode loopParent(const string &fe_name,
954 ForcesAndSourcesCore *parent_fe,
955 const int verb = QUIET,
956 const LogManager::SeverityLevel sev = Sev::noisy);
957
958 /**
959 * @brief User calls this function to loop over parent elements. This function
960 * calls finite element with its operator to do calculations.
961 *
962 * @param fe_name
963 * @param child_fe
964 * @param verb
965 * @param sev
966 * @return MoFEMErrorCode
967 */
968 MoFEMErrorCode loopChildren(const string &fe_name,
969 ForcesAndSourcesCore *child_fe,
970 const int verb = QUIET,
971 const LogManager::SeverityLevel sev = Sev::noisy);
972
973 /**
974 * @brief Iterate over range of elements
975 *
976 * @param fe_name
977 * @param range_fe
978 * @param fe_range
979 * @param verb
980 * @param sev
981 * @return * MoFEMErrorCode
982 */
983 MoFEMErrorCode loopRange(const string &fe_name,
984 ForcesAndSourcesCore *range_fe,
985 boost::shared_ptr<Range> fe_range,
986 const int verb = QUIET,
987 const LogManager::SeverityLevel sev = Sev::noisy);
988
989 /**@}*/
990
991 inline ForcesAndSourcesCore *getPtrFE() const;
992
993 inline ForcesAndSourcesCore *getSidePtrFE() const;
994
995 inline ForcesAndSourcesCore *getRefinePtrFE() const;
996
997
998protected:
1000
1002
1003private:
1008};
1009
1010/// \deprecated Used ForcesAndSourcesCore instead
1012
1013boost::shared_ptr<const NumeredEntFiniteElement>
1014ForcesAndSourcesCore::UserDataOperator::getNumeredEntFiniteElementPtr() const {
1016};
1017
1018EntityHandle ForcesAndSourcesCore::UserDataOperator::getFEEntityHandle() const {
1019 return getNumeredEntFiniteElementPtr()->getEnt();
1020}
1021
1022int ForcesAndSourcesCore::UserDataOperator::getFEDim() const {
1024};
1025
1026EntityType ForcesAndSourcesCore::UserDataOperator::getFEType() const {
1028};
1029
1030boost::weak_ptr<SideNumber>
1031ForcesAndSourcesCore::UserDataOperator::getSideNumberPtr(
1032 const int side_number, const EntityType type) {
1033 auto &side_table_by_side_and_type =
1034 ptrFE->numeredEntFiniteElementPtr->getSideNumberTable().get<1>();
1035 auto side_it =
1036 side_table_by_side_and_type.find(boost::make_tuple(type, side_number));
1037 if (side_it != side_table_by_side_and_type.end())
1038 return *side_it;
1039 else
1040 return boost::weak_ptr<SideNumber>();
1041}
1042
1044ForcesAndSourcesCore::UserDataOperator::getSideEntity(const int side_number,
1045 const EntityType type) {
1046 if (auto side_ptr = getSideNumberPtr(side_number, type).lock())
1047 return side_ptr->ent;
1048 else
1049 return 0;
1050}
1051
1052int ForcesAndSourcesCore::UserDataOperator::getNumberOfNodesOnElement() const {
1053 return ptrFE->getNumberOfNodes();
1054}
1055
1056const FEMethod *ForcesAndSourcesCore::UserDataOperator::getFEMethod() const {
1057 return ptrFE;
1058}
1059
1060int ForcesAndSourcesCore::UserDataOperator::getOpType() const { return opType; }
1061
1062void ForcesAndSourcesCore::UserDataOperator::setOpType(const OpType type) {
1063 opType = type;
1064}
1065
1066void ForcesAndSourcesCore::UserDataOperator::addOpType(const OpType type) {
1067 opType |= type;
1068}
1069
1070int ForcesAndSourcesCore::UserDataOperator::getNinTheLoop() const {
1071 return getFEMethod()->getNinTheLoop();
1072}
1073
1074int ForcesAndSourcesCore::UserDataOperator::getLoopSize() const {
1075 return getFEMethod()->getLoopSize();
1076}
1077
1078std::string ForcesAndSourcesCore::UserDataOperator::getFEName() const {
1079 return getFEMethod()->getFEName();
1080}
1081
1082const PetscData::Switches &
1083ForcesAndSourcesCore::UserDataOperator::getDataCtx() const {
1084 return getFEMethod()->data_ctx;
1085}
1086
1088ForcesAndSourcesCore::UserDataOperator::getKSPCtx() const {
1089 return getFEMethod()->ksp_ctx;
1090}
1091
1093ForcesAndSourcesCore::UserDataOperator::getSNESCtx() const {
1094 return getFEMethod()->snes_ctx;
1095}
1096
1098ForcesAndSourcesCore::UserDataOperator::getTSCtx() const {
1099 return getFEMethod()->ts_ctx;
1100}
1101
1102Vec ForcesAndSourcesCore::UserDataOperator::getKSPf() const {
1103#ifndef NDEBUG
1104 if (getFEMethod()->ksp_f == PETSC_NULLPTR)
1105 THROW_MESSAGE("KSP not set F vector");
1106#endif
1107 return getFEMethod()->ksp_f;
1108}
1109
1110Mat ForcesAndSourcesCore::UserDataOperator::getKSPA() const {
1111#ifndef NDEBUG
1112 if (getFEMethod()->ksp_A == PETSC_NULLPTR)
1113 THROW_MESSAGE("KSP not set A vector");
1114#endif
1115 return getFEMethod()->ksp_A;
1116}
1117
1118Mat ForcesAndSourcesCore::UserDataOperator::getKSPB() const {
1119#ifndef NDEBUG
1120 if (getFEMethod()->ksp_B == PETSC_NULLPTR)
1121 THROW_MESSAGE("KSP not set B vector");
1122#endif
1123 return getFEMethod()->ksp_B;
1124}
1125
1126Vec ForcesAndSourcesCore::UserDataOperator::getSNESf() const {
1127#ifndef NDEBUG
1128 if (getFEMethod()->snes_f == PETSC_NULLPTR)
1129 THROW_MESSAGE("SNES not set F vector");
1130#endif
1131 return getFEMethod()->snes_f;
1132}
1133
1134Vec ForcesAndSourcesCore::UserDataOperator::getSNESx() const {
1135#ifndef NDEBUG
1136 if (getFEMethod()->snes_x == PETSC_NULLPTR)
1137 THROW_MESSAGE("SNESnot set X vector");
1138#endif
1139 return getFEMethod()->snes_x;
1140}
1141
1142Mat ForcesAndSourcesCore::UserDataOperator::getSNESA() const {
1143#ifndef NDEBUG
1144 if (getFEMethod()->snes_A == PETSC_NULLPTR)
1145 THROW_MESSAGE("SNES not set A vector");
1146#endif
1147 return getFEMethod()->snes_A;
1148}
1149
1150Mat ForcesAndSourcesCore::UserDataOperator::getSNESB() const {
1151#ifndef NDEBUG
1152 if (getFEMethod()->snes_B == PETSC_NULLPTR)
1153 THROW_MESSAGE("SNES not set A matrix");
1154#endif
1155 return getFEMethod()->snes_B;
1156}
1157
1158Vec ForcesAndSourcesCore::UserDataOperator::getTSu() const {
1159#ifndef NDEBUG
1160 if (getFEMethod()->ts_u == PETSC_NULLPTR)
1161 THROW_MESSAGE("TS not set U vector");
1162#endif
1163 return getFEMethod()->ts_u;
1164}
1165
1166Vec ForcesAndSourcesCore::UserDataOperator::getTSu_t() const {
1167#ifndef NDEBUG
1168 if (getFEMethod()->ts_u_t == PETSC_NULLPTR)
1169 THROW_MESSAGE("TS not set U_t vector");
1170#endif
1171 return getFEMethod()->ts_u_t;
1172}
1173
1174Vec ForcesAndSourcesCore::UserDataOperator::getTSu_tt() const {
1175#ifndef NDEBUG
1176 if (getFEMethod()->ts_u_tt == PETSC_NULLPTR)
1177 THROW_MESSAGE("TS not set U_tt vector");
1178#endif
1179 return getFEMethod()->ts_u_tt;
1180}
1181
1182Vec ForcesAndSourcesCore::UserDataOperator::getTSf() const {
1183#ifndef NDEBUG
1184 if (getFEMethod()->ts_F == PETSC_NULLPTR)
1185 THROW_MESSAGE("TS not set F vector");
1186#endif
1187 return getFEMethod()->ts_F;
1188}
1189
1190Mat ForcesAndSourcesCore::UserDataOperator::getTSA() const {
1191#ifndef NDEBUG
1192 if (getFEMethod()->ts_A == PETSC_NULLPTR)
1193 THROW_MESSAGE("TS not set A matrix");
1194#endif
1195 return getFEMethod()->ts_A;
1196}
1197
1198Mat ForcesAndSourcesCore::UserDataOperator::getTSB() const {
1199#ifndef NDEBUG
1200 if (getFEMethod()->ts_B == PETSC_NULLPTR)
1201 THROW_MESSAGE("TS not set B matrix");
1202#endif
1203 return getFEMethod()->ts_B;
1204}
1205
1206int ForcesAndSourcesCore::UserDataOperator::getTSstep() const {
1207#ifndef NDEBUG
1208 if ((getFEMethod()->data_ctx & PetscData::PetscData::CtxSetTime).none())
1209 THROW_MESSAGE("TS not set time");
1210#endif
1211 return getFEMethod()->ts_step;
1212}
1213
1214double ForcesAndSourcesCore::UserDataOperator::getTStime() const {
1215#ifndef NDEBUG
1216 if ((getFEMethod()->data_ctx & PetscData::PetscData::CtxSetTime).none())
1217 THROW_MESSAGE("TS not set time");
1218#endif
1219 return getFEMethod()->ts_t;
1220}
1221
1222double ForcesAndSourcesCore::UserDataOperator::getTStimeStep() const {
1223#ifndef NDEBUG
1224 if ((getFEMethod()->data_ctx & PetscData::PetscData::CtxSetTime).none())
1225 THROW_MESSAGE("TS not set time");
1226#endif
1227 return getFEMethod()->ts_dt;
1228}
1229
1230double ForcesAndSourcesCore::UserDataOperator::getTSa() const {
1231#ifndef NDEBUG
1232 if ((getFEMethod()->data_ctx & (PetscData::CtxSetA | PetscData::CtxSetB))
1233 .none() ||
1234 (getFEMethod()->data_ctx & (PetscData::CtxSetX_T)).none())
1235 THROW_MESSAGE("TS not set B matrix");
1236#endif
1237 return getFEMethod()->ts_a;
1238}
1239
1240double ForcesAndSourcesCore::UserDataOperator::getTSaa() const {
1241#ifndef NDEBUG
1242 if ((getFEMethod()->data_ctx & (PetscData::CtxSetA | PetscData::CtxSetB))
1243 .none() ||
1244 (getFEMethod()->data_ctx & (PetscData::CtxSetX_TT)).none())
1245 THROW_MESSAGE("TS not set B matrix");
1246#endif
1247 return getFEMethod()->ts_aa;
1248}
1249
1250MatrixDouble &ForcesAndSourcesCore::UserDataOperator::getGaussPts() {
1251 return static_cast<ForcesAndSourcesCore *>(ptrFE)->gaussPts;
1252}
1253
1254auto ForcesAndSourcesCore::UserDataOperator::getFTensor0IntegrationWeight() {
1256 &(getGaussPts()(getGaussPts().size1() - 1, 0)));
1257}
1258
1259// MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::getPorblemRowIndices(
1260// const std::string filed_name, const EntityType type, const int side,
1261// VectorInt &indices) const {
1262// return getProblemRowIndices(filed_name, type, side, indices);
1263// }
1264
1265ForcesAndSourcesCore *ForcesAndSourcesCore::UserDataOperator::getPtrFE() const {
1266 return ptrFE;
1267}
1268
1270ForcesAndSourcesCore::UserDataOperator::getSidePtrFE() const {
1271 return ptrFE->sidePtrFE;
1272}
1273
1275ForcesAndSourcesCore::UserDataOperator::getRefinePtrFE() const {
1276 return ptrFE->refinePtrFE;
1277}
1278
1279MatrixDouble &ForcesAndSourcesCore::UserDataOperator::getCoordsAtGaussPts() {
1280 return static_cast<ForcesAndSourcesCore *>(ptrFE)->coordsAtGaussPts;
1281}
1282
1283auto ForcesAndSourcesCore::UserDataOperator::getFTensor1CoordsAtGaussPts() {
1285 &getCoordsAtGaussPts()(0, 0), &getCoordsAtGaussPts()(0, 1),
1286 &getCoordsAtGaussPts()(0, 2));
1287}
1288
1289double ForcesAndSourcesCore::UserDataOperator::getMeasure() const {
1290 return static_cast<ForcesAndSourcesCore *>(ptrFE)->elementMeasure;
1291}
1292
1293double &ForcesAndSourcesCore::UserDataOperator::getMeasure() {
1294 return static_cast<ForcesAndSourcesCore *>(ptrFE)->elementMeasure;
1295}
1296
1297/**
1298 * @brief Element used to execute operators on side of the element
1299 *
1300 * @tparam E template for side element type
1301 *
1302 */
1303template <typename E>
1305
1307
1308 OpLoopSide(MoFEM::Interface &m_field, const std::string fe_name,
1309 const int side_dim, boost::shared_ptr<Range> fe_range,
1310 const LogManager::SeverityLevel sev = Sev::noisy,
1311 boost::shared_ptr<AdjCache> adj_cache = nullptr)
1312 : UserDataOperator(NOSPACE, OPSPACE), sideFEPtr(new E(m_field)),
1313 sideFEName(fe_name), sideDim(side_dim), sevLevel(sev),
1314 adjCache(adj_cache), feRange(fe_range) {}
1315
1316 /**
1317 * @brief Construct a new Op Loop Side object
1318 *
1319 * @param m_field
1320 * @param fe_name name of side (domain element)
1321 * @param side_dim dimension
1322 * @param sev severity level
1323 * @param adj_cache if set to null cache of the element is used
1324 */
1325 OpLoopSide(MoFEM::Interface &m_field, const std::string fe_name,
1326 const int side_dim,
1327 const LogManager::SeverityLevel sev = Sev::noisy,
1328 boost::shared_ptr<AdjCache> adj_cache = nullptr)
1329 : OpLoopSide(m_field, fe_name, side_dim, nullptr, sev, adj_cache) {}
1330
1331 MoFEMErrorCode doWork(int side, EntityType type,
1333 return loopSide(sideFEName, sideFEPtr.get(), sideDim, 0, feRange, VERBOSE,
1334 sevLevel, adjCache.get());
1335 };
1336
1337 boost::ptr_deque<UserDataOperator> &getOpPtrVector() {
1338 return sideFEPtr->getOpPtrVector();
1339 }
1340
1341 boost::shared_ptr<E> &getSideFEPtr() { return sideFEPtr; }
1342
1343protected:
1344 const std::string sideFEName;
1345 const int sideDim;
1346 boost::shared_ptr<E> sideFEPtr;
1348 boost::shared_ptr<AdjCache> adjCache;
1349 boost::shared_ptr<Range> feRange;
1350};
1351
1352/**
1353 * @brief Iterate over range of (sub)elements
1354 *
1355 * @tparam E
1356 */
1357template <typename E>
1359
1361
1362 OpLoopRange(MoFEM::Interface &m_field, const std::string fe_name,
1363 boost::shared_ptr<Range> fe_range,
1364 const LogManager::SeverityLevel sev = Sev::noisy)
1365 : UserDataOperator(NOSPACE, OPSPACE), rangeFEPtr(new E(m_field)),
1366 rangeFEName(fe_name), sevLevel(sev), feRange(fe_range) {}
1367
1368 MoFEMErrorCode doWork(int side, EntityType type,
1370 return loopRange(rangeFEName, rangeFEPtr.get(), feRange, VERBOSE, sevLevel);
1371 };
1372
1373 boost::ptr_deque<UserDataOperator> &getOpPtrVector() {
1374 return rangeFEPtr->getOpPtrVector();
1375 }
1376
1377 boost::shared_ptr<E> &getRangeFEPtr() { return rangeFEPtr; }
1378
1379protected:
1380 const std::string rangeFEName;
1381 boost::shared_ptr<E> rangeFEPtr;
1383 boost::shared_ptr<Range> feRange;
1384};
1385
1386/**
1387 * \brief Copy geometry-related data from one element to other
1388 *
1389 * That can be used to copy high order geometry data from coarse element to
1390 * children. That is often a case when higher order geometry is defined only on
1391 * coarse elements.
1392 *
1393 * \note Integration points have to be located at the same geometric positions
1394 *
1395 * FIXME: Write atom test
1396 */
1397template <int DIM>
1399
1400} // namespace MoFEM
1401
1402#endif //__FORCES_AND_SOURCES_CORE__HPP__
1403
1404/**
1405 * \defgroup mofem_forces_and_sources Forces and sources
1406 * \ingroup mofem
1407 *
1408 * \brief Manages complexities related to assembly of vector and matrices at
1409 * single finite element level.
1410 *
1411 **/
@ 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
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
Base face element used to integrate on skeleton.
void setOpType(const OpType type)
Set operator type.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
std::map< EntityHandle, std::vector< boost::weak_ptr< NumeredEntFiniteElement > > > AdjCache
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.
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 form 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()()
function is run for every finite element
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.
boost::function< int(int order_row, int order_col, int order_data)> RuleHookFun
auto & getElementPolynomialBase()
Get the Entity Polynomial Base object.
ForcesAndSourcesCore * sidePtrFE
Element to integrate on the sides.
MatrixDouble coordsAtGaussPts
coordinated at gauss points
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.
boost::function< MoFEMErrorCode(ForcesAndSourcesCore *fe_raw_ptr, int order_row, int order_col, int order_data)> GaussHookFun
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()
function is run at the beginning of loop
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.
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()
function is run at the end of loop
KSPContext
pass information about context of KSP/DM for with finite element is computed
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
static constexpr Switches CtxSetX_TT
std::bitset< 8 > Switches
static constexpr Switches CtxSetX_T
static constexpr Switches CtxSetB
Vec & snes_f
residual
Vec & snes_x
state vector
Mat & snes_B
preconditioner of jacobian matrix
Mat & snes_A
jacobian matrix
Vec & ts_F
residual vector
Vec & ts_u_tt
second time derivative of state vector
Vec & ts_u_t
time derivative of state vector
Vec & ts_u
state vector
Mat & ts_B
Preconditioner for ts_A.
Base volume element used to integrate on contact surface (could be extended to other volume elements)...