v0.13.1
ForcesAndSourcesCore.hpp
Go to the documentation of this file.
1/** \file ForcesAndSourcesCore.hpp
2 \brief Implementation of elements on entities.
3
4 Those element are inherited by user to implement specific implementation of
5 particular problem.
6
7*/
8
9/* This file is part of MoFEM.
10 * MoFEM is free software: you can redistribute it and/or modify it under
11 * the terms of the GNU Lesser General Public License as published by the
12 * Free Software Foundation, either version 3 of the License, or (at your
13 * option) any later version.
14 *
15 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18 * License for more details.
19 *
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
22
23#ifndef __FORCES_AND_SOURCES_CORE__HPP__
24#define __FORCES_AND_SOURCES_CORE__HPP__
25
26using namespace boost::numeric;
27
28namespace MoFEM {
29
30/** \brief structure to get information form mofem into EntitiesFieldData
31 * \ingroup mofem_forces_and_sources
32 *
33 */
35
37
39 typedef boost::function<int(int order_row, int order_col, int order_data)>
41
42 typedef boost::function<MoFEMErrorCode(ForcesAndSourcesCore *fe_raw_ptr,
43 int order_row, int order_col,
44 int order_data)>
46
47 /**
48 * \brief Hook to get rule
49 *
50 * \todo check preferred format how works with gcc and clang,
51 * see
52 * <http://www.boost.org/doc/libs/1_64_0/doc/html/function/tutorial.html#idp247873024>
53 */
55
56 /**
57 * @brief Set function to calculate integration rule
58 *
59 */
61
62 /** \brief Data operator to do calculations at integration points.
63 * \ingroup mofem_forces_and_sources
64
65 Is inherited and implemented by user to do calculations. It can be used in
66 many different ways but typically is used to integrate matrices (f.e.
67 stiffness matrix) and the right hand vector (f.e. force vector).
68
69 Note: It is assumed that operator is executed for symmetric problem. That
70 means that is executed for not repeating entities on finite element. For
71 example on triangle we have nodes, 3 edges and one face. Because of symmetry
72 calculations are for: nodes-nodes, nodes-edge0, nodes-edge_1, nodes-edge_2,
73 nodes-face,
74 edges_1-edges_1, edges_1-edge_1, edge_1-edge_2,
75 edge_1-edge_1, edge_1-edge_2,
76 edge_2-edge_2,
77 edge_1-face, edge_1-face, edges_3-face,
78 face - face
79
80 In case of non symmetric problem in addition calculations of lower off
81 diagonal terms. F.e. edge_1-edge_0, esges_3-edge_0, edge_3-edge_1,
82
83 In that case class variable UserDataOperator::sYmm = false;
84
85 NoteL: By default sYmm is set for symmetric problems
86
87 */
88 struct UserDataOperator;
89
90 /** \brief Use to push back operator for row operator
91
92 It can be used to calculate nodal forces or other quantities on the mesh.
93
94 */
95 boost::ptr_vector<UserDataOperator> &getOpPtrVector() { return opPtrVector; }
96
97 /**
98 * @brief Get the Entity Polynomial Base object
99 *
100 * @return boost::shared_ptr<BaseFunction>&&
101 */
103
104 /**
105 * @brief Get the User Polynomial Base object
106 *
107 * @return boost::shared_ptr<BaseFunction>&
108 */
110
111 /**
112 * @brief Matrix of integration points
113 *
114 * Columns is equal to number of integration points, numver of rows
115 * depends on dimension of finite element entity, for example for
116 * tetrahedron rows are x,y,z,weight. Last row is integration weight.
117 *
118 * FIXME: that should be moved to private class data and acessed only by
119 * member function
120 */
122
123 virtual MoFEMErrorCode preProcess();
124 virtual MoFEMErrorCode operator()();
125 virtual MoFEMErrorCode postProcess();
126
127public:
128 /** \brief Get max order of approximation for data fields
129
130 getMeasure getMaxDataOrder () return maximal order on entities, for
131 all data on the element. So for example if finite element is triangle, and
132 triangle base function have order 4 and on edges base function have order
133 2, this function return 4.
134
135 If finite element has for example 2 or more approximated fields, for
136 example Pressure (order 3) and displacement field (order 5), this function
137 returns 5.
138
139 */
140 int getMaxDataOrder() const;
141
142 /// \brief Get max order of approximation for field in rows
143 int getMaxRowOrder() const;
144
145 /// \brief Get max order of approximation for field in columns
146 int getMaxColOrder() const;
147
148 /**
149 * @brief Get the entity data
150 *
151 * @param space
152 * @param type
153 * @param side
154 * @return EntitiesFieldData::EntData&
155 */
156 auto &getEntData(const FieldSpace space, const EntityType type,
157 const int side) {
158 return dataOnElement[space]->dataOnEntities[type][side];
159 }
160
161 /**
162 * @brief Get data on entities and space
163 *
164 * Entities data are stored by space, by entity type, and entity side.
165 *
166 * @return std::array<boost::shared_ptr<EntitiesFieldData>, LASTSPACE>
167 */
169
170 /**
171 * @brief Get derived data on entities and space
172 *
173 * Entities data are stored by space, by entity type, and entity side. Derived
174 * data is used to store data on columns, so it shares infromatin about shape
175 * functions wih rows.
176 *
177 * @return std::array<boost::shared_ptr<EntitiesFieldData>, LASTSPACE>
178 */
180
181protected:
182 /**
183 * \brief get sense (orientation) of entity
184 * @param type type of entity
185 * @param data entity data
186 * @return error code
187 */
189 const EntityType type,
190 boost::ptr_vector<EntitiesFieldData::EntData> &data) const;
191
192 /**
193 * @brief Get the entity data order
194 *
195 * @param type
196 * @param space
197 * @param data
198 * @return MoFEMErrorCode
199 */
201 const EntityType type, const FieldSpace space,
202 boost::ptr_vector<EntitiesFieldData::EntData> &data) const;
203
204 /**
205 * @brief Get the entity sense (orientation)
206 *
207 * @tparam type
208 * @param data
209 * @return MoFEMErrorCode
210 */
211 template <EntityType type>
213 return getEntitySense(type, data.dataOnEntities[type]);
214 }
215
216 /**
217 * @brief Get the entity data order for given space
218 *
219 * @tparam type
220 * @param data
221 * @param space
222 * @return MoFEMErrorCode
223 */
224 template <EntityType type>
226 const FieldSpace space) const {
227 return getEntityDataOrder(type, space, data.dataOnEntities[type]);
228 }
229
230 /** \name Indices */
231
232 /**@{*/
233
234 /// \brief get node indices
235 template <typename EXTRACTOR>
237 getNodesIndices(const int bit_number,
238 FieldEntity_vector_view &ents_field, VectorInt &nodes_indices,
239 VectorInt &local_nodes_indices, EXTRACTOR &&extractor) const;
240
241 /// \brief get row node indices from FENumeredDofEntity_multiIndex
243 const int bit_number) const;
244
245 /// \brief get col node indices from FENumeredDofEntity_multiIndex
247 const int bit_number) const;
248
249 template <typename EXTRACTOR>
250 MoFEMErrorCode getEntityIndices(EntitiesFieldData &data, const int bit_number,
251 FieldEntity_vector_view &ents_field,
252 const EntityType type_lo,
253 const EntityType type_hi,
254 EXTRACTOR &&extractor) const;
255
257 getEntityRowIndices(EntitiesFieldData &data, const int bit_number,
258 const EntityType type_lo = MBVERTEX,
259 const EntityType type_hi = MBPOLYHEDRON) const;
260
262 getEntityColIndices(EntitiesFieldData &data, const int bit_number,
263 const EntityType type_lo = MBVERTEX,
264 const EntityType type_hi = MBPOLYHEDRON) const;
265
266 /// \brief get NoField indices
268 getNoFieldIndices(const int bit_number,
269 boost::shared_ptr<FENumeredDofEntity_multiIndex> dofs,
270 VectorInt &nodes_indices) const;
271
272 /// \brief get col NoField indices
274 const int bit_number) const;
275
276 /// \brief get col NoField indices
278 const int bit_number) const;
279
280 /**@}*/
281
282 /** \name Data */
283
284 /**@{*/
285
286 /** Get bit ref level in entities, and set it to data
287 */
289
290 /**
291 * \brief Get field data on nodes
292 */
293 MoFEMErrorCode getNoFieldFieldData(const int bit_number,
294 VectorDouble &ent_field_data,
295 VectorDofs &ent_field_dofs,
296 VectorFieldEntities &ent_field) const;
297
299 const int bit_number) const;
300
301 /**
302 * \brief Get data on nodes
303 * @param data Data structure
304 * @param field_name Field name
305 * @return Error code
306 */
308 const int bit_number) const;
309
311 getEntityFieldData(EntitiesFieldData &data, const int bit_number,
312 const EntityType type_lo = MBVERTEX,
313 const EntityType type_hi = MBPOLYHEDRON) const;
314
315 /**@}*/
316
317 /// \brief Get nodes on faces
319
320 /// \brief Get field approximation space and base on entities
323
324 /** \name Data form NumeredDofEntity_multiIndex */
325
326 /**@{*/
327
328 /// \brief get indices of nodal indices which are declared for problem but
329 /// not this particular element
331 const NumeredDofEntity_multiIndex &dofs,
332 VectorInt &nodes_indices) const;
333
334 /// \brief get indices by type (generic function) which are declared for
335 /// problem but not this particular element
337 const NumeredDofEntity_multiIndex &dofs,
338 EntityType type, int side_number,
339 VectorInt &indices) const;
340
342 VectorInt &nodes_indices) const;
344 EntityType type, int side_number,
345 VectorInt &indices) const;
347 VectorInt &nodes_indices) const;
349 EntityType type, int side_number,
350 VectorInt &indices) const;
351
352 /**@}*/
353
354 /**
355 * \brief another variant of getRule
356 * @param order_row order of base function on row
357 * @param order_col order of base function on columns
358 * @param order_data order of base function approximating data
359 * @return integration rule
360 *
361 * This function is overloaded by the user. The integration rule
362 * is set such that specific operator implemented by the user is
363 * integrated accurately. For example if user implement bilinear operator
364 * \f[
365 * b(u,v) =
366 * \int_\mathcal{T}
367 * \frac{\partial u_i}{\partial x_j}\frac{\partial v_i}{\partial x_j}
368 * \textrm{d}\mathcal{T}
369 * \f]
370 * then if \f$u\f$ and \f$v\f$ are polynomial of given \em order, then
371 * exact integral would be \code int getRule(int order) { return
372 * 2*(order-1); }; \endcode
373 *
374 * The integration points and weights are set appropriately for given
375 * entity type and integration rule from \ref quad.c
376 *
377 * Method \ref ForcesAndSourcesCore::getRule takes at argument takes
378 * maximal polynomial order set on the element on all fields defined on
379 * the element. If a user likes to have more control, another variant of
380 * this function can be called which distinguishing between field orders
381 * on rows, columns and data, the i.e. first argument of a bilinear form,
382 * the second argument of bilinear form and field coefficients on the
383 * element.
384 *
385 * \note If user set rule to -1 or any other negative integer, then method
386 * \ref ForcesAndSourcesCore::setGaussPts is called. In that method user
387 * can implement own (specific) integration method.
388 *
389 * \bug this function should be const
390 */
391 virtual int getRule(int order_row, int order_col, int order_data);
392
393 /** \brief set user specific integration rule
394
395 This function allows for user defined integration rule. The key is to
396 called matrix gaussPts, which is used by other MoFEM procedures. Matrix
397 has number of rows equal to problem dimension plus one, where last index
398 is used to store weight values. %Number of columns is equal to number of
399 integration points.
400
401 \note This function is called if method \ref
402 ForcesAndSourcesCore::getRule is returning integer -1 or any other
403 negative integer.
404
405 User sets
406 \code
407 MatrixDouble gaussPts;
408 \endcode
409 where
410 \code
411 gaussPts.resize(dim+1,nb_gauss_pts);
412 \endcode
413 number rows represents local coordinates of integration points
414 in reference element, where last index in row is for integration weight.
415
416 */
417 virtual MoFEMErrorCode setGaussPts(int order_row, int order_col,
418 int order_data);
419
420 /**
421 * \brief Calculate base functions
422 * @return Error code
423 */
426
427 /**
428 * \brief Calculate base functions
429 * @return Error code
430 */
432
433 /**
434 * @brief Calculate Bernstein-Bezier base
435 *
436 * @return MoFEMErrorCode
437 */
439
440 /**
441 * @brief Create a entity data on element object
442 *
443 * @return MoFEMErrorCode
444 */
446
447 /**
448 * @brief Iterate user data operators
449 *
450 * @return MoFEMErrorCode
451 */
453
454 /**@{*/
455
456 /** \name Deprecated (do not use) */
457
458 /** \deprecated Use getRule(int row_order, int col_order, int data order)
459 */
460 virtual int getRule(int order);
461
462 /** \deprecated setGaussPts(int row_order, int col_order, int data order);
463 */
464 virtual MoFEMErrorCode setGaussPts(int order);
465
466 /**@}*/
467
468 /**
469 * @brief Entity data on element entity rows fields
470 *
471 */
472 const std::array<boost::shared_ptr<EntitiesFieldData>, LASTSPACE>
474
475 /**
476 * @brief Entity data on element entity columns fields
477 *
478 */
479 const std::array<boost::shared_ptr<EntitiesFieldData>, LASTSPACE>
481
487
488 /**
489 * @brief Vector of finite element users data operators
490 *
491 */
492 boost::ptr_vector<UserDataOperator> opPtrVector;
493
494 friend class UserDataOperator;
495
496protected:
497 /**
498 * @brief Last evaluated type of element entity
499 *
500 */
502
503private:
504 /**
505 * @brief Pointer to entity polynomial base
506 *
507 */
508 boost::shared_ptr<BaseFunction> elementPolynomialBasePtr;
509
510 /**
511 * @brief Pointer to user polynomail base
512 */
513 boost::shared_ptr<BaseFunction> userPolynomialBasePtr;
514
515 /**
516 * @brief Element to integrate on the sides
517 *
518 */
520
521 /**
522 * @brief Set the pointer to face element on the side
523 *
524 * \note Function is is used by face element, while it iterates over
525 * elements on the side
526 *
527 * @param side_fe_ptr
528 * @return MoFEMErrorCode
529 */
531
532 /**u
533 * @brief Element to integrate parent or child
534 *
535 */
537
538 /**
539 * @brief Set the pointer to face element refined
540 *
541 * @param refine_fe_ptr
542 * @return MoFEMErrorCode
543 */
545
551
552protected:
553 MatrixDouble coordsAtGaussPts; ///< coordinated at gauss points
554 double elementMeasure; ///< Depending on dimension of elements, stores length,
555 ///< area or volume of element.
556};
557
559
560 /**
561 * \brief Controls loop over entities on element
562 *
563 * - OPROW is used if row vector is assembled
564 * - OPCOL is usually used if column vector is assembled
565 * - OPROWCOL is usually used for assemble matrices.
566 * - OPSPACE no field is defined for such operator. Is usually used to modify
567 * base
568 *
569 *
570 * For typical problem like Bubnov-Galerkin OPROW and OPCOL are the same. In
571 * more general case for example for non-square matrices columns and rows
572 * could have different numeration and/or different set of fields.
573 *
574 */
575 enum OpType {
576 OPROW = 1 << 0, ///< operator doWork function is executed on FE rows
577 OPCOL = 1 << 1, ///< operator doWork function is executed on FE columns
578 OPROWCOL = 1 << 2, ///< operator doWork is executed on FE rows &columns
579 OPSPACE = 1 << 3, ///< operator do Work is execute on space data
580 OPLAST = 1 << 3 ///< @deprecated would be removed
581 };
582
583 static const char *const OpTypeNames[];
584
585 char opType;
586 std::string rowFieldName;
587 std::string colFieldName;
589
590 /**
591 * This Constructor is used typically when some modification base shape
592 * functions on some approx. space is applied. Operator is run for all data
593 * on space.
594 *
595 * User has no access to field data from this operator.
596 */
597 UserDataOperator(const FieldSpace space, const char type = OPSPACE,
598 const bool symm = true);
599
600 UserDataOperator(const std::string &field_name, const char type,
601 const bool symm = true);
602
603 UserDataOperator(const std::string &row_field_name,
604 const std::string &col_field_name, const char type,
605 const bool symm = true);
606
607 /** \brief Return raw pointer to NumeredEntFiniteElement
608 */
609 inline boost::shared_ptr<const NumeredEntFiniteElement>
611
612 /**
613 * \brief Return finite element entity handle
614 * @return Finite element entity handle
615 */
616 inline EntityHandle getFEEntityHandle() const;
617
618 /**
619 * @brief Get dimension of finite element
620 *
621 * @return int
622 */
623 inline int getFEDim() const;
624
625 /**
626 * @brief Get dimension of finite element
627 *
628 * @return int
629 */
630 inline EntityType getFEType() const;
631
632 /**
633 * @brief Get the side number pointer
634 *
635 * \note For vertex is expection. Side basses in argument of function doWork
636 * is zero. For other entity types side can be used as argument of this
637 * function.
638 *
639 * @param side_number
640 * @param type
641 * @return boost::weak_ptr<SideNumber>
642 */
643 inline boost::weak_ptr<SideNumber> getSideNumberPtr(const int side_number,
644 const EntityType type);
645
646 /**
647 * @brief Get the side entity
648 *
649 * \note For vertex is expection. Side basses in argument of function
650 * doWork is zero. For other entity types side can be used as argument of
651 * this function.
652 *
653 * \code
654 * MoFEMErrorCode doWork(int side, EntityType type,
655 * EntitiesFieldData::EntData &data) {
656 * MoFEMFunctionBegin;
657 *
658 * if (type == MBVERTEX) {
659 * for (int n = 0; n != number_of_nodes; ++n)
660 * EntityHandle ent = getSideEntity(n, type);
661 *
662 * // Do somthing
663 *
664 * } else {
665 * EntityHandle ent = getSideEntity(side, type);
666 *
667 * // Do somthing
668 *
669 * }
670 *
671 * MoFEMFunctionReturn(0);
672 * }
673 * \endcode
674 *
675 * @param side_number
676 * @param type
677 */
678 inline EntityHandle getSideEntity(const int side_number,
679 const EntityType type);
680
681 /**
682 * @brief Get the number of nodes on finite element
683 *
684 * @return int
685 */
686 inline int getNumberOfNodesOnElement() const;
687
688 /** \brief Get row indices
689
690 Field could be or not declared for this element but is declared for
691 problem
692
693 \param field_name
694 \param type entity type
695 \param side side number, any number if type is MBVERTEX
696 \return indices
697
698 NOTE: Using those indices to assemble matrix will result in error if new
699 non-zero values need to be created.
700
701 */
702 MoFEMErrorCode getProblemRowIndices(const std::string filed_name,
703 const EntityType type, const int side,
704 VectorInt &indices) const;
705
706 /** \brief Get col indices
707
708 Field could be or not declared for this element but is declared for
709 problem
710
711 \param field_name
712 \param type entity type
713 \param side side number, any number if type is MBVERTEX
714 \return indices
715
716 NOTE: Using those indices to assemble matrix will result in error if new
717 non-zero values need to be created.
718
719 */
720 MoFEMErrorCode getProblemColIndices(const std::string filed_name,
721 const EntityType type, const int side,
722 VectorInt &indices) const;
723
724 /** \brief Return raw pointer to Finite Element Method object
725 */
726 inline const FEMethod *getFEMethod() const;
727
728 /**
729 * \brief Get operator types
730 * @return Return operator type
731 */
732 inline int getOpType() const;
733
734 /**
735 * \brief Set operator type
736 * @param Operator type
737 */
738 inline void setOpType(const OpType type);
739
740 /**
741 * \brief Add operator type
742 */
743 inline void addOpType(const OpType type);
744
745 /**
746 * \brief get number of finite element in the loop
747 * @return number of finite element
748 */
749 inline int getNinTheLoop() const;
750
751 /**
752 * \brief get size of elements in the loop
753 * @return loop size
754 */
755 inline int getLoopSize() const;
756
757 /** \brief Get name of the element
758 */
759 inline std::string getFEName() const;
760
761 /** \name Accessing KSP */
762
763 /**@{*/
764
765 inline const PetscData::Switches &getDataCtx() const;
766
767 inline const KspMethod::KSPContext getKSPCtx() const;
768
769 inline const SnesMethod::SNESContext getSNESCtx() const;
770
771 inline const TSMethod::TSContext getTSCtx() const;
772
773 /**@}*/
774
775 /**@{*/
776
777 inline Vec getKSPf() const;
778
779 inline Mat getKSPA() const;
780
781 inline Mat getKSPB() const;
782
783 /**@}*/
784
785 /** \name Accessing SNES */
786
787 /**@{*/
788
789 inline Vec getSNESf() const;
790
791 inline Vec getSNESx() const;
792
793 inline Mat getSNESA() const;
794
795 inline Mat getSNESB() const;
796
797 /**@}*/
798
799 /** \name Accessing TS */
800
801 /**@{*/
802
803 inline Vec getTSu() const;
804
805 inline Vec getTSu_t() const;
806
807 inline Vec getTSu_tt() const;
808
809 inline Vec getTSf() const;
810
811 inline Mat getTSA() const;
812
813 inline Mat getTSB() const;
814
815 inline int getTSstep() const;
816
817 inline double getTStime() const;
818
819 inline double getTSa() const;
820
821 inline double getTSaa() const;
822
823 /**@}*/
824
825 /**@{*/
826
827 /** \name Base funtions and integration points */
828
829 /** \brief matrix of integration (Gauss) points for Volume Element
830 *
831 * For triangle: columns 0,1 are x,y coordinates respectively and column
832 * 2 is a weight value for example getGaussPts()(1,13) returns y
833 * coordinate of 13th Gauss point on particular volume element
834 *
835 * For tetrahedron: columns 0,1,2 are x,y,z coordinates respectively and
836 * column 3 is a weight value for example getGaussPts()(1,13) returns y
837 * coordinate of 13th Gauss point on particular volume element
838 *
839 */
840 inline MatrixDouble &getGaussPts();
841
842 /**
843 * @brief Get integration weights
844 *
845 * \code
846 * auto t_w = getFTensor0IntegrationWeight();
847 * for(int gg = 0; gg!=getGaussPts.size2(); ++gg) {
848 * // integrate something
849 * ++t_w;
850 * }
851 * \endcode
852 *
853 * @return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>
854 */
855 inline auto getFTensor0IntegrationWeight();
856
857 /**@}*/
858
859 /** \name Coordinates and access to internal data */
860
861 /**@{*/
862
863 /** \brief Gauss points and weight, matrix (nb. of points x 3)
864
865 Column 0-2 integration points coordinate x, y and z, respectively. At rows
866 are integration points.
867
868 */
870
871 /**
872 * \brief Get coordinates at integration points assuming linear geometry
873 *
874 * \code
875 * auto t_coords = getFTensor1CoordsAtGaussPts();
876 * for(int gg = 0;gg!=nb_int_ptrs;gg++) {
877 * // do something
878 * ++t_coords;
879 * }
880 * \endcode
881 *
882 */
883 inline auto getFTensor1CoordsAtGaussPts();
884
885 /**@}*/
886
887 /**@{*/
888
889 /** \name Measures (area, volume, length, etc.) */
890
891 /**
892 * \brief get measure of element
893 * @return volume
894 */
895 inline double getMeasure() const;
896
897 /**
898 * \brief get measure of element
899 * @return volume
900 */
901 inline double &getMeasure();
902
903 /**}*/
904
905 /**@{*/
906
907 /** \name Deprecated (do not use) */
908
909 // //! \deprecated Use getSNESF intead
910 // DEPRECATED inline Vec getSnesF() const { return getSNESf(); }
911
912 // //! \deprecated Use getSNESX intead
913 // DEPRECATED inline Vec getSnesX() const { return getSNESx(); }
914
915 // //! \deprecated Use getSNESA intead
916 // DEPRECATED inline Mat getSnesA() const { return getSNESA(); }
917
918 // //! \deprecated Use getSNESB intead
919 // DEPRECATED inline Mat getSnesB() const { return getSNESB(); }
920
921 // // \deprecated Deprecated function with spelling mistake
922 // DEPRECATED inline MoFEMErrorCode
923 // getPorblemRowIndices(const std::string filed_name, const EntityType type,
924 // const int side, VectorInt &indices) const;
925
926 /**@}*/
927
928 /**@{*/
929
930 /** \name Loops */
931
932 /**
933 * @brief User call this function to loop over elements on the side of
934 * face. This function calls finite element with is operator to do
935 * calculations.
936 *
937 * @param fe_name name of the side element
938 * @param side_fe pointer to the side element instance
939 * @param dim dimension the of side element
940 * @param ent_for_side entity handle for which adjacent volume or face will
941 * be accessed
942 * @param verb
943 * @param sev
944 * @return MoFEMErrorCode
945 */
946 MoFEMErrorCode loopSide(const string &fe_name, ForcesAndSourcesCore *side_fe,
947 const size_t dim, const EntityHandle ent_for_side = 0,
948 const int verb = QUIET,
949 const LogManager::SeverityLevel sev = Sev::noisy);
950
951 /**
952 * @brief User call this function to loop over parent elements. This function
953 * calls finite element with is operator to do calculations.
954 *
955 * @param fe_name
956 * @param parent_fe
957 * @param verb
958 * @param sev
959 * @return MoFEMErrorCode
960 */
961 MoFEMErrorCode loopThis(const string &fe_name,
962 ForcesAndSourcesCore *parent_fe,
963 const int verb = QUIET,
964 const LogManager::SeverityLevel sev = Sev::noisy);
965
966 /**
967 * @brief User call this function to loop over parent elements. This function
968 * calls finite element with is operator to do calculations.
969 *
970 * @param fe_name
971 * @param parent_fe
972 * @param verb
973 * @param sev
974 * @return MoFEMErrorCode
975 */
976 MoFEMErrorCode loopParent(const string &fe_name,
977 ForcesAndSourcesCore *parent_fe,
978 const int verb = QUIET,
979 const LogManager::SeverityLevel sev = Sev::noisy);
980
981 /**
982 * @brief User call this function to loop over parent elements. This function
983 * calls finite element with is operator to do calculations.
984 *
985 * @param fe_name
986 * @param child_fe
987 * @param verb
988 * @param sev
989 * @return MoFEMErrorCode
990 */
991 MoFEMErrorCode loopChildren(const string &fe_name,
992 ForcesAndSourcesCore *child_fe,
993 const int verb = QUIET,
994 const LogManager::SeverityLevel sev = Sev::noisy);
995
996 /**@}*/
997
998 inline ForcesAndSourcesCore *getPtrFE() const;
999
1000 inline ForcesAndSourcesCore *getSidePtrFE() const;
1001
1002 inline ForcesAndSourcesCore *getRefinePtrFE() const;
1003
1004
1005protected:
1007
1009
1010private:
1015};
1016
1017/// \deprecated Used ForcesAndSourcesCore instead
1019
1020boost::shared_ptr<const NumeredEntFiniteElement>
1023};
1024
1026 return getNumeredEntFiniteElementPtr()->getEnt();
1027}
1028
1031};
1032
1035};
1036
1037boost::weak_ptr<SideNumber>
1039 const int side_number, const EntityType type) {
1040 auto &side_table_by_side_and_type =
1041 ptrFE->numeredEntFiniteElementPtr->getSideNumberTable().get<1>();
1042 auto side_it =
1043 side_table_by_side_and_type.find(boost::make_tuple(type, side_number));
1044 if (side_it != side_table_by_side_and_type.end())
1045 return *side_it;
1046 else
1047 return boost::weak_ptr<SideNumber>();
1048}
1049
1050EntityHandle
1052 const EntityType type) {
1053 if (auto side_ptr = getSideNumberPtr(side_number, type).lock())
1054 return side_ptr->ent;
1055 else
1056 return 0;
1057}
1058
1060 return ptrFE->getNumberOfNodes();
1061}
1062
1064 return ptrFE;
1065}
1066
1068
1070 opType = type;
1071}
1072
1074 opType |= type;
1075}
1076
1078 return getFEMethod()->getNinTheLoop();
1079}
1080
1082 return getFEMethod()->getLoopSize();
1083}
1084
1086 return getFEMethod()->getFEName();
1087}
1088
1089const PetscData::Switches &
1091 return getFEMethod()->data_ctx;
1092}
1093
1096 return getFEMethod()->ksp_ctx;
1097}
1098
1101 return getFEMethod()->snes_ctx;
1102}
1103
1106 return getFEMethod()->ts_ctx;
1107}
1108
1110#ifndef NDEBUG
1111 if (getFEMethod()->ksp_f == PETSC_NULL)
1112 THROW_MESSAGE("KSP not set F vector");
1113#endif
1114 return getFEMethod()->ksp_f;
1115}
1116
1118#ifndef NDEBUG
1119 if (getFEMethod()->ksp_A == PETSC_NULL)
1120 THROW_MESSAGE("KSP not set A vector");
1121#endif
1122 return getFEMethod()->ksp_A;
1123}
1124
1126#ifndef NDEBUG
1127 if (getFEMethod()->ksp_B == PETSC_NULL)
1128 THROW_MESSAGE("KSP not set B vector");
1129#endif
1130 return getFEMethod()->ksp_B;
1131}
1132
1134#ifndef NDEBUG
1135 if (getFEMethod()->snes_f == PETSC_NULL)
1136 THROW_MESSAGE("SNES not set F vector");
1137#endif
1138 return getFEMethod()->snes_f;
1139}
1140
1142#ifndef NDEBUG
1143 if (getFEMethod()->snes_x == PETSC_NULL)
1144 THROW_MESSAGE("SNESnot set X vector");
1145#endif
1146 return getFEMethod()->snes_x;
1147}
1148
1150#ifndef NDEBUG
1151 if (getFEMethod()->snes_A == PETSC_NULL)
1152 THROW_MESSAGE("SNES not set A vector");
1153#endif
1154 return getFEMethod()->snes_A;
1155}
1156
1158#ifndef NDEBUG
1159 if (getFEMethod()->snes_B == PETSC_NULL)
1160 THROW_MESSAGE("SNES not set A matrix");
1161#endif
1162 return getFEMethod()->snes_B;
1163}
1164
1166#ifndef NDEBUG
1167 if (getFEMethod()->ts_u == PETSC_NULL)
1168 THROW_MESSAGE("TS not set U vector");
1169#endif
1170 return getFEMethod()->ts_u;
1171}
1172
1174#ifndef NDEBUG
1175 if (getFEMethod()->ts_u_t == PETSC_NULL)
1176 THROW_MESSAGE("TS not set U_t vector");
1177#endif
1178 return getFEMethod()->ts_u_t;
1179}
1180
1182#ifndef NDEBUG
1183 if (getFEMethod()->ts_u_tt == PETSC_NULL)
1184 THROW_MESSAGE("TS not set U_tt vector");
1185#endif
1186 return getFEMethod()->ts_u_tt;
1187}
1188
1190#ifndef NDEBUG
1191 if (getFEMethod()->ts_F == PETSC_NULL)
1192 THROW_MESSAGE("TS not set F vector");
1193#endif
1194 return getFEMethod()->ts_F;
1195}
1196
1198#ifndef NDEBUG
1199 if (getFEMethod()->ts_A == PETSC_NULL)
1200 THROW_MESSAGE("TS not set A matrix");
1201#endif
1202 return getFEMethod()->ts_A;
1203}
1204
1206#ifndef NDEBUG
1207 if (getFEMethod()->ts_B == PETSC_NULL)
1208 THROW_MESSAGE("TS not set B matrix");
1209#endif
1210 return getFEMethod()->ts_B;
1211}
1212
1214#ifndef NDEBUG
1215 if ((getFEMethod()->data_ctx & PetscData::PetscData::CtxSetTime).none())
1216 THROW_MESSAGE("TS not set time");
1217#endif
1218 return getFEMethod()->ts_step;
1219}
1220
1222#ifndef NDEBUG
1223 if ((getFEMethod()->data_ctx & PetscData::PetscData::CtxSetTime).none())
1224 THROW_MESSAGE("TS not set time");
1225#endif
1226 return getFEMethod()->ts_t;
1227}
1228
1230#ifndef NDEBUG
1231 if ((getFEMethod()->data_ctx & (PetscData::CtxSetA | PetscData::CtxSetB))
1232 .none() ||
1233 (getFEMethod()->data_ctx & (PetscData::CtxSetX_T)).none())
1234 THROW_MESSAGE("TS not set B matrix");
1235#endif
1236 return getFEMethod()->ts_a;
1237}
1238
1240#ifndef NDEBUG
1241 if ((getFEMethod()->data_ctx & (PetscData::CtxSetA | PetscData::CtxSetB))
1242 .none() ||
1243 (getFEMethod()->data_ctx & (PetscData::CtxSetX_TT)).none())
1244 THROW_MESSAGE("TS not set B matrix");
1245#endif
1246 return getFEMethod()->ts_aa;
1247}
1248
1250 return static_cast<ForcesAndSourcesCore *>(ptrFE)->gaussPts;
1251}
1252
1255 &(getGaussPts()(getGaussPts().size1() - 1, 0)));
1256}
1257
1258// MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::getPorblemRowIndices(
1259// const std::string filed_name, const EntityType type, const int side,
1260// VectorInt &indices) const {
1261// return getProblemRowIndices(filed_name, type, side, indices);
1262// }
1263
1265 return ptrFE;
1266}
1267
1270 return ptrFE->sidePtrFE;
1271}
1272
1275 return ptrFE->refinePtrFE;
1276}
1277
1279 return static_cast<ForcesAndSourcesCore *>(ptrFE)->coordsAtGaussPts;
1280}
1281
1284 &getCoordsAtGaussPts()(0, 0), &getCoordsAtGaussPts()(0, 1),
1285 &getCoordsAtGaussPts()(0, 2));
1286}
1287
1289 return static_cast<ForcesAndSourcesCore *>(ptrFE)->elementMeasure;
1290}
1291
1293 return static_cast<ForcesAndSourcesCore *>(ptrFE)->elementMeasure;
1294}
1295
1296} // namespace MoFEM
1297
1298#endif //__FORCES_AND_SOURCES_CORE__HPP__
1299
1300/**
1301 * \defgroup mofem_forces_and_sources Forces and sources
1302 * \ingroup mofem
1303 *
1304 * \brief Manages complexities related to assembly of vector and matrices at
1305 * single finite element level.
1306 *
1307 **/
@ QUIET
Definition: definitions.h:221
FieldApproximationBase
approximation base
Definition: definitions.h:71
FieldSpace
approximation spaces
Definition: definitions.h:95
@ LASTSPACE
FieldSpace in [ 0, LASTSPACE )
Definition: definitions.h:102
#define DEPRECATED
Definition: definitions.h:30
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:574
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.
Definition: LogManager.hpp:43
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasVector< int > VectorInt
Definition: Types.hpp:78
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1487
ublas::vector< FieldEntity *, FieldEntAllocator > VectorFieldEntities
DEPRECATED typedef ForcesAndSourcesCore ForcesAndSurcesCore
auto dimension_from_handle(const EntityHandle h)
get entity dimension form handle
Definition: Templates.hpp:1495
std::vector< boost::weak_ptr< FieldEntity > > FieldEntity_vector_view
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
constexpr auto field_name
base operator to do operations at Gauss Pt. level
Deprecated interface functions.
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.
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.
const KspMethod::KSPContext getKSPCtx() const
UserDataOperator(const FieldSpace space, const char type=OPSPACE, const bool symm=true)
int getNinTheLoop() const
get number of finite element in the loop
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
MoFEMErrorCode loopThis(const string &fe_name, ForcesAndSourcesCore *parent_fe, const int verb=QUIET, const LogManager::SeverityLevel sev=Sev::noisy)
User call this function to loop over parent elements. This function calls finite element with is oper...
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 call this function to loop over parent elements. This function calls finite element with is oper...
const PetscData::Switches & getDataCtx() const
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 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)
User call this function to loop over elements on the side of face. This function calls finite 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 call this function to loop over parent elements. This function calls finite element with is oper...
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
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.
ForcesAndSourcesCore(Interface &m_field)
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
MatrixDouble coordsAtGaussPts
coordinated at gauss points
boost::function< MoFEMErrorCode(ForcesAndSourcesCore *fe_raw_ptr, int order_row, int order_col, int order_data)> GaussHookFun
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
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
boost::ptr_vector< UserDataOperator > opPtrVector
Vector of finite element users data operators.
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
GaussHookFun setRuleHook
Set function to calculate integration rule.
MoFEMErrorCode getBitRefLevelOnData()
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.
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
Definition: LoopMethods.hpp:83
static constexpr Switches CtxSetA
Definition: LoopMethods.hpp:47
static constexpr Switches CtxSetX_TT
Definition: LoopMethods.hpp:51
std::bitset< 8 > Switches
Definition: LoopMethods.hpp:43
static constexpr Switches CtxSetX_T
Definition: LoopMethods.hpp:50
static constexpr Switches CtxSetB
Definition: LoopMethods.hpp:48
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)...
Base volume element used to integrate on skeleton.