v0.14.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 infromatin 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 /// \brief get NoField indices
256 getNoFieldIndices(const int bit_number,
257 boost::shared_ptr<FENumeredDofEntity_multiIndex> dofs,
258 VectorInt &nodes_indices) const;
259
260 /// \brief get col NoField indices
262 const int bit_number) const;
263
264 /// \brief get col NoField indices
266 const int bit_number) const;
267
268 /**@}*/
269
270 /** \name Data */
271
272 /**@{*/
273
274 /** Get bit ref level in entities, and set it to data
275 */
277
278 /**
279 * \brief Get field data on nodes
280 */
281 MoFEMErrorCode getNoFieldFieldData(const int bit_number,
282 VectorDouble &ent_field_data,
283 VectorDofs &ent_field_dofs,
284 VectorFieldEntities &ent_field) const;
285
287 const int bit_number) const;
288
289 /**
290 * \brief Get data on nodes
291 * @param data Data structure
292 * @param field_name Field name
293 * @return Error code
294 */
296 const int bit_number) const;
297
299 getEntityFieldData(EntitiesFieldData &data, const int bit_number,
300 const EntityType type_lo = MBVERTEX,
301 const EntityType type_hi = MBPOLYHEDRON) const;
302
303 /**@}*/
304
305 /// \brief Get nodes on faces
307
308 /// \brief Get field approximation space and base on entities
311
312 /** \name Data form NumeredDofEntity_multiIndex */
313
314 /**@{*/
315
316 /// \brief get indices of nodal indices which are declared for problem but
317 /// not this particular element
319 const NumeredDofEntity_multiIndex &dofs,
320 VectorInt &nodes_indices) const;
321
322 /// \brief get indices by type (generic function) which are declared for
323 /// problem but not this particular element
325 const NumeredDofEntity_multiIndex &dofs,
326 EntityType type, int side_number,
327 VectorInt &indices) const;
328
330 VectorInt &nodes_indices) const;
332 EntityType type, int side_number,
333 VectorInt &indices) const;
335 VectorInt &nodes_indices) const;
337 EntityType type, int side_number,
338 VectorInt &indices) const;
339
340 /**@}*/
341
342 /**
343 * \brief another variant of getRule
344 * @param order_row order of base function on row
345 * @param order_col order of base function on columns
346 * @param order_data order of base function approximating data
347 * @return integration rule
348 *
349 * This function is overloaded by the user. The integration rule
350 * is set such that specific operator implemented by the user is
351 * integrated accurately. For example if user implement bilinear operator
352 * \f[
353 * b(u,v) =
354 * \int_\mathcal{T}
355 * \frac{\partial u_i}{\partial x_j}\frac{\partial v_i}{\partial x_j}
356 * \textrm{d}\mathcal{T}
357 * \f]
358 * then if \f$u\f$ and \f$v\f$ are polynomial of given \em order, then
359 * exact integral would be \code int getRule(int order) { return
360 * 2*(order-1); }; \endcode
361 *
362 * The integration points and weights are set appropriately for given
363 * entity type and integration rule from \ref quad.c
364 *
365 * Method \ref ForcesAndSourcesCore::getRule takes at argument takes
366 * maximal polynomial order set on the element on all fields defined on
367 * the element. If a user likes to have more control, another variant of
368 * this function can be called which distinguishing between field orders
369 * on rows, columns and data, the i.e. first argument of a bilinear form,
370 * the second argument of bilinear form and field coefficients on the
371 * element.
372 *
373 * \note If user set rule to -1 or any other negative integer, then method
374 * \ref ForcesAndSourcesCore::setGaussPts is called. In that method user
375 * can implement own (specific) integration method.
376 *
377 * \bug this function should be const
378 */
379 virtual int getRule(int order_row, int order_col, int order_data);
380
381 /** \brief set user specific integration rule
382
383 This function allows for user defined integration rule. The key is to
384 called matrix gaussPts, which is used by other MoFEM procedures. Matrix
385 has number of rows equal to problem dimension plus one, where last index
386 is used to store weight values. %Number of columns is equal to number of
387 integration points.
388
389 \note This function is called if method \ref
390 ForcesAndSourcesCore::getRule is returning integer -1 or any other
391 negative integer.
392
393 User sets
394 \code
395 MatrixDouble gaussPts;
396 \endcode
397 where
398 \code
399 gaussPts.resize(dim+1,nb_gauss_pts);
400 \endcode
401 number rows represents local coordinates of integration points
402 in reference element, where last index in row is for integration weight.
403
404 */
405 virtual MoFEMErrorCode setGaussPts(int order_row, int order_col,
406 int order_data);
407
408 /**
409 * \brief Calculate base functions
410 * @return Error code
411 */
414
415 /**
416 * \brief Calculate base functions
417 * @return Error code
418 */
420
421 /**
422 * @brief Calculate Bernstein-Bezier base
423 *
424 * @return MoFEMErrorCode
425 */
427
428 /**
429 * @brief Create a entity data on element object
430 *
431 * @return MoFEMErrorCode
432 */
434
435 /**
436 * @brief Iterate user data operators
437 *
438 * @return MoFEMErrorCode
439 */
441
442 /**@{*/
443
444 /** \name Deprecated (do not use) */
445
446 /** \deprecated Use getRule(int row_order, int col_order, int data order)
447 */
448 virtual int getRule(int order);
449
450 /** \deprecated setGaussPts(int row_order, int col_order, int data order);
451 */
452 virtual MoFEMErrorCode setGaussPts(int order);
453
454 /**@}*/
455
456 /**
457 * @brief Entity data on element entity rows fields
458 *
459 */
460 const std::array<boost::shared_ptr<EntitiesFieldData>, LASTSPACE>
462
463 /**
464 * @brief Entity data on element entity columns fields
465 *
466 */
467 const std::array<boost::shared_ptr<EntitiesFieldData>, LASTSPACE>
469
475
476 /**
477 * @brief Vector of finite element users data operators
478 *
479 */
480 boost::ptr_deque<UserDataOperator> opPtrVector;
481
482 friend class UserDataOperator;
483
484protected:
485 /**
486 * @brief Last evaluated type of element entity
487 *
488 */
490
491private:
492 /**
493 * @brief Pointer to entity polynomial base
494 *
495 */
496 boost::shared_ptr<BaseFunction> elementPolynomialBasePtr;
497
498 /**
499 * @brief Pointer to user polynomail base
500 */
501 boost::shared_ptr<BaseFunction> userPolynomialBasePtr;
502
503 /**
504 * @brief Element to integrate on the sides
505 *
506 */
508
509 /**
510 * @brief Set the pointer to face element on the side
511 *
512 * \note Function is is used by face element, while it iterates over
513 * elements on the side
514 *
515 * @param side_fe_ptr
516 * @return MoFEMErrorCode
517 */
519
520 /**u
521 * @brief Element to integrate parent or child
522 *
523 */
525
526 /**
527 * @brief Set the pointer to face element refined
528 *
529 * @param refine_fe_ptr
530 * @return MoFEMErrorCode
531 */
533
539
540 template <int DIM> friend struct OpCopyGeomDataToE;
541
542
543protected:
544 MatrixDouble coordsAtGaussPts; ///< coordinated at gauss points
545 double elementMeasure; ///< Depending on dimension of elements, stores length,
546 ///< area or volume of element.
547};
548
550
551 /**
552 * \brief Controls loop over entities on element
553 *
554 * - OPROW is used if row vector is assembled
555 * - OPCOL is usually used if column vector is assembled
556 * - OPROWCOL is usually used for assemble matrices.
557 * - OPSPACE no field is defined for such operator. Is usually used to modify
558 * base
559 *
560 *
561 * For typical problem like Bubnov-Galerkin OPROW and OPCOL are the same. In
562 * more general case for example for non-square matrices columns and rows
563 * could have different numeration and/or different set of fields.
564 *
565 */
566 enum OpType {
567 OPROW = 1 << 0, ///< operator doWork function is executed on FE rows
568 OPCOL = 1 << 1, ///< operator doWork function is executed on FE columns
569 OPROWCOL = 1 << 2, ///< operator doWork is executed on FE rows &columns
570 OPSPACE = 1 << 3, ///< operator do Work is execute on space data
571 OPLAST = 1 << 3 ///< @deprecated would be removed
572 };
573
574 static const char *const OpTypeNames[];
575
576 char opType;
577 std::string rowFieldName;
578 std::string colFieldName;
580
581 /**
582 * This Constructor is used typically when some modification base shape
583 * functions on some approx. space is applied. Operator is run for all data
584 * on space.
585 *
586 * User has no access to field data from this operator.
587 */
588 UserDataOperator(const FieldSpace space, const char type = OPSPACE,
589 const bool symm = true);
590
591 UserDataOperator(const std::string field_name, const char type,
592 const bool symm = true);
593
594 UserDataOperator(const std::string row_field_name,
595 const std::string col_field_name, const char type,
596 const bool symm = true);
597
598 /** \brief Return raw pointer to NumeredEntFiniteElement
599 */
600 inline boost::shared_ptr<const NumeredEntFiniteElement>
602
603 /**
604 * \brief Return finite element entity handle
605 * @return Finite element entity handle
606 */
607 inline EntityHandle getFEEntityHandle() const;
608
609 /**
610 * @brief Get dimension of finite element
611 *
612 * @return int
613 */
614 inline int getFEDim() const;
615
616 /**
617 * @brief Get dimension of finite element
618 *
619 * @return int
620 */
621 inline EntityType getFEType() const;
622
623 /**
624 * @brief Get the side number pointer
625 *
626 * \note For vertex is expectation. Side basses in argument of function doWork
627 * is zero. For other entity types side can be used as argument of this
628 * function.
629 *
630 * @param side_number
631 * @param type
632 * @return boost::weak_ptr<SideNumber>
633 */
634 inline boost::weak_ptr<SideNumber> getSideNumberPtr(const int side_number,
635 const EntityType type);
636
637 /**
638 * @brief Get the side entity
639 *
640 * \note For vertex is expectation. Side basses in argument of function
641 * doWork is zero. For other entity types side can be used as argument of
642 * this function.
643 *
644 * \code
645 * MoFEMErrorCode doWork(int side, EntityType type,
646 * EntitiesFieldData::EntData &data) {
647 * MoFEMFunctionBegin;
648 *
649 * if (type == MBVERTEX) {
650 * for (int n = 0; n != number_of_nodes; ++n)
651 * EntityHandle ent = getSideEntity(n, type);
652 *
653 * // Do somthing
654 *
655 * } else {
656 * EntityHandle ent = getSideEntity(side, type);
657 *
658 * // Do somthing
659 *
660 * }
661 *
662 * MoFEMFunctionReturn(0);
663 * }
664 * \endcode
665 *
666 * @param side_number
667 * @param type
668 */
669 inline EntityHandle getSideEntity(const int side_number,
670 const EntityType type);
671
672 /**
673 * @brief Get the number of nodes on finite element
674 *
675 * @return int
676 */
677 inline int getNumberOfNodesOnElement() const;
678
679 /** \brief Get row indices
680
681 Field could be or not declared for this element but is declared for
682 problem
683
684 \param field_name
685 \param type entity type
686 \param side side number, any number if type is MBVERTEX
687 \return indices
688
689 NOTE: Using those indices to assemble matrix will result in error if new
690 non-zero values need to be created.
691
692 */
693 MoFEMErrorCode getProblemRowIndices(const std::string filed_name,
694 const EntityType type, const int side,
695 VectorInt &indices) const;
696
697 /** \brief Get col indices
698
699 Field could be or not declared for this element but is declared for
700 problem
701
702 \param field_name
703 \param type entity type
704 \param side side number, any number if type is MBVERTEX
705 \return indices
706
707 NOTE: Using those indices to assemble matrix will result in error if new
708 non-zero values need to be created.
709
710 */
711 MoFEMErrorCode getProblemColIndices(const std::string filed_name,
712 const EntityType type, const int side,
713 VectorInt &indices) const;
714
715 /** \brief Return raw pointer to Finite Element Method object
716 */
717 inline const FEMethod *getFEMethod() const;
718
719 /**
720 * \brief Get operator types
721 * @return Return operator type
722 */
723 inline int getOpType() const;
724
725 /**
726 * \brief Set operator type
727 * @param Operator type
728 */
729 inline void setOpType(const OpType type);
730
731 /**
732 * \brief Add operator type
733 */
734 inline void addOpType(const OpType type);
735
736 /**
737 * \brief get number of finite element in the loop
738 * @return number of finite element
739 */
740 inline int getNinTheLoop() const;
741
742 /**
743 * \brief get size of elements in the loop
744 * @return loop size
745 */
746 inline int getLoopSize() const;
747
748 /** \brief Get name of the element
749 */
750 inline std::string getFEName() const;
751
752 /** \name Accessing KSP */
753
754 /**@{*/
755
756 inline const PetscData::Switches &getDataCtx() const;
757
758 inline const KspMethod::KSPContext getKSPCtx() const;
759
760 inline const SnesMethod::SNESContext getSNESCtx() const;
761
762 inline const TSMethod::TSContext getTSCtx() const;
763
764 /**@}*/
765
766 /**@{*/
767
768 inline Vec getKSPf() const;
769
770 inline Mat getKSPA() const;
771
772 inline Mat getKSPB() const;
773
774 /**@}*/
775
776 /** \name Accessing SNES */
777
778 /**@{*/
779
780 inline Vec getSNESf() const;
781
782 inline Vec getSNESx() const;
783
784 inline Mat getSNESA() const;
785
786 inline Mat getSNESB() const;
787
788 /**@}*/
789
790 /** \name Accessing TS */
791
792 /**@{*/
793
794 inline Vec getTSu() const;
795
796 inline Vec getTSu_t() const;
797
798 inline Vec getTSu_tt() const;
799
800 inline Vec getTSf() const;
801
802 inline Mat getTSA() const;
803
804 inline Mat getTSB() const;
805
806 inline int getTSstep() const;
807
808 inline double getTStime() const;
809
810 inline double getTStimeStep() const;
811
812 inline double getTSa() const;
813
814 inline double getTSaa() const;
815
816 /**@}*/
817
818 /**@{*/
819
820 /** \name Base funtions and integration points */
821
822 /** \brief matrix of integration (Gauss) points for Volume Element
823 *
824 * For triangle: columns 0,1 are x,y coordinates respectively and column
825 * 2 is a weight value for example getGaussPts()(1,13) returns y
826 * coordinate of 13th Gauss point on particular volume element
827 *
828 * For tetrahedron: columns 0,1,2 are x,y,z coordinates respectively and
829 * column 3 is a weight value for example getGaussPts()(1,13) returns y
830 * coordinate of 13th Gauss point on particular volume element
831 *
832 */
833 inline MatrixDouble &getGaussPts();
834
835 /**
836 * @brief Get integration weights
837 *
838 * \code
839 * auto t_w = getFTensor0IntegrationWeight();
840 * for(int gg = 0; gg!=getGaussPts.size2(); ++gg) {
841 * // integrate something
842 * ++t_w;
843 * }
844 * \endcode
845 *
846 * @return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>
847 */
848 inline auto getFTensor0IntegrationWeight();
849
850 /**@}*/
851
852 /** \name Coordinates and access to internal data */
853
854 /**@{*/
855
856 /** \brief Gauss points and weight, matrix (nb. of points x 3)
857
858 Column 0-2 integration points coordinate x, y and z, respectively. At rows
859 are integration points.
860
861 */
863
864 /**
865 * \brief Get coordinates at integration points assuming linear geometry
866 *
867 * \code
868 * auto t_coords = getFTensor1CoordsAtGaussPts();
869 * for(int gg = 0;gg!=nb_int_ptrs;gg++) {
870 * // do something
871 * ++t_coords;
872 * }
873 * \endcode
874 *
875 */
876 inline auto getFTensor1CoordsAtGaussPts();
877
878 /**@}*/
879
880 /**@{*/
881
882 /** \name Measures (area, volume, length, etc.) */
883
884 /**
885 * \brief get measure of element
886 * @return volume
887 */
888 inline double getMeasure() const;
889
890 /**
891 * \brief get measure of element
892 * @return volume
893 */
894 inline double &getMeasure();
895
896 /**}*/
897
898 /**@{*/
899
900 /** \name Loops */
901
902 using AdjCache =
903 std::map<EntityHandle,
904 std::vector<boost::weak_ptr<NumeredEntFiniteElement>>>;
905
906 /**
907 * @brief User calls this function to loop over elements on the side of
908 * face. This function calls finite element with its operator to do
909 * calculations.
910 *
911 * @param fe_name name of the side element
912 * @param side_fe pointer to the side element instance
913 * @param dim dimension the of side element
914 * @param ent_for_side entity handle for which adjacent volume or face will
915 * be accessed
916 * @param verb
917 * @param sev
918 * @return MoFEMErrorCode
919 */
920 MoFEMErrorCode loopSide(const string &fe_name, ForcesAndSourcesCore *side_fe,
921 const size_t dim, const EntityHandle ent_for_side = 0,
922 const int verb = QUIET,
923 const LogManager::SeverityLevel sev = Sev::noisy,
924 AdjCache *adj_cache = nullptr
925
926 );
927
928 /**
929 * @brief User calls this function to loop over the same element using a
930 * different set of integration points. This function calls finite element
931 * with its operator to do calculations.
932 *
933 * @param fe_name
934 * @param this_fe
935 * @param verb
936 * @param sev
937 * @return MoFEMErrorCode
938 */
939 MoFEMErrorCode loopThis(const string &fe_name,
940 ForcesAndSourcesCore *this_fe,
941 const int verb = QUIET,
942 const LogManager::SeverityLevel sev = Sev::noisy);
943
944 /**
945 * @brief User calls this function to loop over parent elements. This function
946 * calls finite element with its operator to do calculations.
947 *
948 * @param fe_name
949 * @param parent_fe
950 * @param verb
951 * @param sev
952 * @return MoFEMErrorCode
953 */
954 MoFEMErrorCode loopParent(const string &fe_name,
955 ForcesAndSourcesCore *parent_fe,
956 const int verb = QUIET,
957 const LogManager::SeverityLevel sev = Sev::noisy);
958
959 /**
960 * @brief User calls this function to loop over parent elements. This function
961 * calls finite element with its operator to do calculations.
962 *
963 * @param fe_name
964 * @param child_fe
965 * @param verb
966 * @param sev
967 * @return MoFEMErrorCode
968 */
969 MoFEMErrorCode loopChildren(const string &fe_name,
970 ForcesAndSourcesCore *child_fe,
971 const int verb = QUIET,
972 const LogManager::SeverityLevel sev = Sev::noisy);
973
974 /**@}*/
975
976 inline ForcesAndSourcesCore *getPtrFE() const;
977
978 inline ForcesAndSourcesCore *getSidePtrFE() const;
979
980 inline ForcesAndSourcesCore *getRefinePtrFE() const;
981
982
983protected:
985
987
988private:
993};
994
995/// \deprecated Used ForcesAndSourcesCore instead
997
998boost::shared_ptr<const NumeredEntFiniteElement>
999ForcesAndSourcesCore::UserDataOperator::getNumeredEntFiniteElementPtr() const {
1001};
1002
1003EntityHandle ForcesAndSourcesCore::UserDataOperator::getFEEntityHandle() const {
1004 return getNumeredEntFiniteElementPtr()->getEnt();
1005}
1006
1007int ForcesAndSourcesCore::UserDataOperator::getFEDim() const {
1009};
1010
1011EntityType ForcesAndSourcesCore::UserDataOperator::getFEType() const {
1013};
1014
1015boost::weak_ptr<SideNumber>
1016ForcesAndSourcesCore::UserDataOperator::getSideNumberPtr(
1017 const int side_number, const EntityType type) {
1018 auto &side_table_by_side_and_type =
1019 ptrFE->numeredEntFiniteElementPtr->getSideNumberTable().get<1>();
1020 auto side_it =
1021 side_table_by_side_and_type.find(boost::make_tuple(type, side_number));
1022 if (side_it != side_table_by_side_and_type.end())
1023 return *side_it;
1024 else
1025 return boost::weak_ptr<SideNumber>();
1026}
1027
1029ForcesAndSourcesCore::UserDataOperator::getSideEntity(const int side_number,
1030 const EntityType type) {
1031 if (auto side_ptr = getSideNumberPtr(side_number, type).lock())
1032 return side_ptr->ent;
1033 else
1034 return 0;
1035}
1036
1037int ForcesAndSourcesCore::UserDataOperator::getNumberOfNodesOnElement() const {
1038 return ptrFE->getNumberOfNodes();
1039}
1040
1041const FEMethod *ForcesAndSourcesCore::UserDataOperator::getFEMethod() const {
1042 return ptrFE;
1043}
1044
1045int ForcesAndSourcesCore::UserDataOperator::getOpType() const { return opType; }
1046
1047void ForcesAndSourcesCore::UserDataOperator::setOpType(const OpType type) {
1048 opType = type;
1049}
1050
1051void ForcesAndSourcesCore::UserDataOperator::addOpType(const OpType type) {
1052 opType |= type;
1053}
1054
1055int ForcesAndSourcesCore::UserDataOperator::getNinTheLoop() const {
1056 return getFEMethod()->getNinTheLoop();
1057}
1058
1059int ForcesAndSourcesCore::UserDataOperator::getLoopSize() const {
1060 return getFEMethod()->getLoopSize();
1061}
1062
1063std::string ForcesAndSourcesCore::UserDataOperator::getFEName() const {
1064 return getFEMethod()->getFEName();
1065}
1066
1067const PetscData::Switches &
1068ForcesAndSourcesCore::UserDataOperator::getDataCtx() const {
1069 return getFEMethod()->data_ctx;
1070}
1071
1073ForcesAndSourcesCore::UserDataOperator::getKSPCtx() const {
1074 return getFEMethod()->ksp_ctx;
1075}
1076
1078ForcesAndSourcesCore::UserDataOperator::getSNESCtx() const {
1079 return getFEMethod()->snes_ctx;
1080}
1081
1083ForcesAndSourcesCore::UserDataOperator::getTSCtx() const {
1084 return getFEMethod()->ts_ctx;
1085}
1086
1087Vec ForcesAndSourcesCore::UserDataOperator::getKSPf() const {
1088#ifndef NDEBUG
1089 if (getFEMethod()->ksp_f == PETSC_NULL)
1090 THROW_MESSAGE("KSP not set F vector");
1091#endif
1092 return getFEMethod()->ksp_f;
1093}
1094
1095Mat ForcesAndSourcesCore::UserDataOperator::getKSPA() const {
1096#ifndef NDEBUG
1097 if (getFEMethod()->ksp_A == PETSC_NULL)
1098 THROW_MESSAGE("KSP not set A vector");
1099#endif
1100 return getFEMethod()->ksp_A;
1101}
1102
1103Mat ForcesAndSourcesCore::UserDataOperator::getKSPB() const {
1104#ifndef NDEBUG
1105 if (getFEMethod()->ksp_B == PETSC_NULL)
1106 THROW_MESSAGE("KSP not set B vector");
1107#endif
1108 return getFEMethod()->ksp_B;
1109}
1110
1111Vec ForcesAndSourcesCore::UserDataOperator::getSNESf() const {
1112#ifndef NDEBUG
1113 if (getFEMethod()->snes_f == PETSC_NULL)
1114 THROW_MESSAGE("SNES not set F vector");
1115#endif
1116 return getFEMethod()->snes_f;
1117}
1118
1119Vec ForcesAndSourcesCore::UserDataOperator::getSNESx() const {
1120#ifndef NDEBUG
1121 if (getFEMethod()->snes_x == PETSC_NULL)
1122 THROW_MESSAGE("SNESnot set X vector");
1123#endif
1124 return getFEMethod()->snes_x;
1125}
1126
1127Mat ForcesAndSourcesCore::UserDataOperator::getSNESA() const {
1128#ifndef NDEBUG
1129 if (getFEMethod()->snes_A == PETSC_NULL)
1130 THROW_MESSAGE("SNES not set A vector");
1131#endif
1132 return getFEMethod()->snes_A;
1133}
1134
1135Mat ForcesAndSourcesCore::UserDataOperator::getSNESB() const {
1136#ifndef NDEBUG
1137 if (getFEMethod()->snes_B == PETSC_NULL)
1138 THROW_MESSAGE("SNES not set A matrix");
1139#endif
1140 return getFEMethod()->snes_B;
1141}
1142
1143Vec ForcesAndSourcesCore::UserDataOperator::getTSu() const {
1144#ifndef NDEBUG
1145 if (getFEMethod()->ts_u == PETSC_NULL)
1146 THROW_MESSAGE("TS not set U vector");
1147#endif
1148 return getFEMethod()->ts_u;
1149}
1150
1151Vec ForcesAndSourcesCore::UserDataOperator::getTSu_t() const {
1152#ifndef NDEBUG
1153 if (getFEMethod()->ts_u_t == PETSC_NULL)
1154 THROW_MESSAGE("TS not set U_t vector");
1155#endif
1156 return getFEMethod()->ts_u_t;
1157}
1158
1159Vec ForcesAndSourcesCore::UserDataOperator::getTSu_tt() const {
1160#ifndef NDEBUG
1161 if (getFEMethod()->ts_u_tt == PETSC_NULL)
1162 THROW_MESSAGE("TS not set U_tt vector");
1163#endif
1164 return getFEMethod()->ts_u_tt;
1165}
1166
1167Vec ForcesAndSourcesCore::UserDataOperator::getTSf() const {
1168#ifndef NDEBUG
1169 if (getFEMethod()->ts_F == PETSC_NULL)
1170 THROW_MESSAGE("TS not set F vector");
1171#endif
1172 return getFEMethod()->ts_F;
1173}
1174
1175Mat ForcesAndSourcesCore::UserDataOperator::getTSA() const {
1176#ifndef NDEBUG
1177 if (getFEMethod()->ts_A == PETSC_NULL)
1178 THROW_MESSAGE("TS not set A matrix");
1179#endif
1180 return getFEMethod()->ts_A;
1181}
1182
1183Mat ForcesAndSourcesCore::UserDataOperator::getTSB() const {
1184#ifndef NDEBUG
1185 if (getFEMethod()->ts_B == PETSC_NULL)
1186 THROW_MESSAGE("TS not set B matrix");
1187#endif
1188 return getFEMethod()->ts_B;
1189}
1190
1191int ForcesAndSourcesCore::UserDataOperator::getTSstep() const {
1192#ifndef NDEBUG
1193 if ((getFEMethod()->data_ctx & PetscData::PetscData::CtxSetTime).none())
1194 THROW_MESSAGE("TS not set time");
1195#endif
1196 return getFEMethod()->ts_step;
1197}
1198
1199double ForcesAndSourcesCore::UserDataOperator::getTStime() const {
1200#ifndef NDEBUG
1201 if ((getFEMethod()->data_ctx & PetscData::PetscData::CtxSetTime).none())
1202 THROW_MESSAGE("TS not set time");
1203#endif
1204 return getFEMethod()->ts_t;
1205}
1206
1207double ForcesAndSourcesCore::UserDataOperator::getTStimeStep() const {
1208#ifndef NDEBUG
1209 if ((getFEMethod()->data_ctx & PetscData::PetscData::CtxSetTime).none())
1210 THROW_MESSAGE("TS not set time");
1211#endif
1212 return getFEMethod()->ts_dt;
1213}
1214
1215double ForcesAndSourcesCore::UserDataOperator::getTSa() const {
1216#ifndef NDEBUG
1217 if ((getFEMethod()->data_ctx & (PetscData::CtxSetA | PetscData::CtxSetB))
1218 .none() ||
1219 (getFEMethod()->data_ctx & (PetscData::CtxSetX_T)).none())
1220 THROW_MESSAGE("TS not set B matrix");
1221#endif
1222 return getFEMethod()->ts_a;
1223}
1224
1225double ForcesAndSourcesCore::UserDataOperator::getTSaa() const {
1226#ifndef NDEBUG
1227 if ((getFEMethod()->data_ctx & (PetscData::CtxSetA | PetscData::CtxSetB))
1228 .none() ||
1229 (getFEMethod()->data_ctx & (PetscData::CtxSetX_TT)).none())
1230 THROW_MESSAGE("TS not set B matrix");
1231#endif
1232 return getFEMethod()->ts_aa;
1233}
1234
1235MatrixDouble &ForcesAndSourcesCore::UserDataOperator::getGaussPts() {
1236 return static_cast<ForcesAndSourcesCore *>(ptrFE)->gaussPts;
1237}
1238
1239auto ForcesAndSourcesCore::UserDataOperator::getFTensor0IntegrationWeight() {
1241 &(getGaussPts()(getGaussPts().size1() - 1, 0)));
1242}
1243
1244// MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::getPorblemRowIndices(
1245// const std::string filed_name, const EntityType type, const int side,
1246// VectorInt &indices) const {
1247// return getProblemRowIndices(filed_name, type, side, indices);
1248// }
1249
1250ForcesAndSourcesCore *ForcesAndSourcesCore::UserDataOperator::getPtrFE() const {
1251 return ptrFE;
1252}
1253
1255ForcesAndSourcesCore::UserDataOperator::getSidePtrFE() const {
1256 return ptrFE->sidePtrFE;
1257}
1258
1260ForcesAndSourcesCore::UserDataOperator::getRefinePtrFE() const {
1261 return ptrFE->refinePtrFE;
1262}
1263
1264MatrixDouble &ForcesAndSourcesCore::UserDataOperator::getCoordsAtGaussPts() {
1265 return static_cast<ForcesAndSourcesCore *>(ptrFE)->coordsAtGaussPts;
1266}
1267
1268auto ForcesAndSourcesCore::UserDataOperator::getFTensor1CoordsAtGaussPts() {
1270 &getCoordsAtGaussPts()(0, 0), &getCoordsAtGaussPts()(0, 1),
1271 &getCoordsAtGaussPts()(0, 2));
1272}
1273
1274double ForcesAndSourcesCore::UserDataOperator::getMeasure() const {
1275 return static_cast<ForcesAndSourcesCore *>(ptrFE)->elementMeasure;
1276}
1277
1278double &ForcesAndSourcesCore::UserDataOperator::getMeasure() {
1279 return static_cast<ForcesAndSourcesCore *>(ptrFE)->elementMeasure;
1280}
1281
1282/**
1283 * @brief Element used to execute operators on side of the element
1284 *
1285 * @tparam E template for side element type
1286 *
1287 */
1288template <typename E>
1290
1292
1293 /**
1294 * @brief Construct a new Op Loop Side object
1295 *
1296 * @param m_field
1297 * @param fe_name name of side (domain element)
1298 * @param side_dim dimension
1299 */
1300 OpLoopSide(MoFEM::Interface &m_field, const std::string fe_name,
1301 const int side_dim,
1302 const LogManager::SeverityLevel sev = Sev::noisy,
1303 boost::shared_ptr<AdjCache> adj_cache = nullptr)
1304 : UserDataOperator(NOSPACE, OPSPACE), sideFEPtr(new E(m_field)),
1305 fieldName(fe_name), sideDim(side_dim), sevLevel(sev),
1306 adjCache(adj_cache) {}
1307
1311 CHKERR loopSide(fieldName, sideFEPtr.get(), sideDim, 0, VERBOSE, sevLevel,
1312 adjCache.get());
1314 };
1315
1316 boost::ptr_deque<UserDataOperator> &getOpPtrVector() {
1317 return sideFEPtr->getOpPtrVector();
1318 }
1319
1320 boost::shared_ptr<E> &getSideFEPtr() { return sideFEPtr; }
1321
1322protected:
1323 const std::string fieldName;
1324 const int sideDim;
1325 boost::shared_ptr<E> sideFEPtr;
1327 boost::shared_ptr<AdjCache> adjCache;
1328};
1329
1330/**
1331 * \brief Copy geometry-related data from one element to other
1332 *
1333 * That can be used to copy high order geometry data from coarse element to
1334 * children. That is often a case when higher order geometry is defined only on
1335 * coarse elements.
1336 *
1337 * \note Integration points have to be located at the same gometric positions
1338 *
1339 * FIXME: Write atom test
1340 */
1341template <int DIM> struct OpCopyGeomDataToE;
1342
1343} // namespace MoFEM
1344
1345#endif //__FORCES_AND_SOURCES_CORE__HPP__
1346
1347/**
1348 * \defgroup mofem_forces_and_sources Forces and sources
1349 * \ingroup mofem
1350 *
1351 * \brief Manages complexities related to assembly of vector and matrices at
1352 * single finite element level.
1353 *
1354 **/
@ 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 MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define DEPRECATED
Definition definitions.h:17
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
constexpr int order
const int dim
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
ublas::vector< FieldEntity *, FieldEntAllocator > VectorFieldEntities
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
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.
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.
MoFEMErrorCode loopSide(const string &fe_name, ForcesAndSourcesCore *side_fe, const size_t dim, const EntityHandle ent_for_side=0, 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...
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
MoFEMErrorCode getNoFieldIndices(const int bit_number, boost::shared_ptr< FENumeredDofEntity_multiIndex > dofs, VectorInt &nodes_indices) const
get NoField indices
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.
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 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.
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
boost::function< MoFEMErrorCode(ForcesAndSourcesCore *fe_raw_ptr, int order_row, int order_col, int order_data) GaussHookFun)
int getMaxDataOrder() const
Get max order of approximation for data fields.
MoFEMErrorCode getNoFieldFieldData(const int bit_number, VectorDouble &ent_field_data, VectorDofs &ent_field_dofs, VectorFieldEntities &ent_field) const
Get field data on nodes.
boost::shared_ptr< BaseFunction > userPolynomialBasePtr
Pointer to user polynomail 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
Copy geometry-related data from one element to other.
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< E > & getSideFEPtr()
boost::shared_ptr< AdjCache > adjCache
boost::shared_ptr< E > sideFEPtr
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
const std::string fieldName
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)...