v0.13.1
Loading...
Searching...
No Matches
EntitiesFieldData.hpp
Go to the documentation of this file.
1/** \file EntitiesFieldData.hpp
2
3\brief Data structures for accessing information about finite element and its
4degrees of freedom.
5
6*/
7
8
9
10#ifndef __ENTITIES_FIELD_DATA_HPP__
11#define __ENTITIES_FIELD_DATA_HPP__
12
13using namespace boost::numeric;
14
15namespace MoFEM {
16
17using DofsAllocator = ublas::unbounded_array<
18
19 FEDofEntity *, std::allocator<FEDofEntity *>
20
21 >;
22
23using VectorDofs = ublas::vector<FEDofEntity *, DofsAllocator>;
24
25using FieldEntAllocator = ublas::unbounded_array<
26
27 FieldEntity *, std::allocator<FieldEntity *>
28
29 >;
30
31using VectorFieldEntities = ublas::vector<FieldEntity *, FieldEntAllocator>;
32
33/** \brief data structure for finite element entity
34 * \ingroup mofem_forces_and_sources_user_data_operators
35 *
36 * It keeps that about indices of degrees of freedom, dofs data, base functions
37 * functions, entity side number, type of entities, approximation order, etc.
38 *
39 */
41
42 struct EntData;
43
44 std::bitset<LASTSPACE> sPace; ///< spaces on element
45 std::bitset<LASTBASE> bAse; ///< bases on element
46 MatrixInt facesNodes; ///< nodes on finite element faces
47 MatrixInt facesNodesOrder; ///< order of face nodes on element
48
49 std::array<std::bitset<LASTSPACE>, MBMAXTYPE>
50 spacesOnEntities; ///< spaces on entity types
51 std::array<std::bitset<LASTBASE>, MBMAXTYPE>
52 basesOnEntities; ///< bases on entity types
53 std::array<std::bitset<LASTBASE>, LASTSPACE>
54 basesOnSpaces; ///< base on spaces
55 std::array<boost::ptr_vector<EntData>, MBMAXTYPE>
56 dataOnEntities; ///< data on nodes, base
57 ///< function, dofs
58 ///< values, etc.
59
60 EntitiesFieldData(const EntityType type);
61 virtual ~EntitiesFieldData() = default;
62
63 virtual MoFEMErrorCode setElementType(const EntityType type);
64
65 /**
66 * Reset data associated with particular field name
67 * @return error code
68 */
70
71 /**
72 * @brief Swap approximation base
73 *
74 * Bernstein-Bezier (BB) base is not hierarchical, and is calculated for
75 * particular field, since it all shape functions change with the order. BB
76 * base is precalculated for every field, and when user push operator with
77 * paricular field using BB base, pointers to shape funtions and
78 * BaseDerivatives of shape functions are set to particular location, once
79 * operator is executed, pointers are switch back to its oroginal position.
80 *
81 * getNSharedPtr(base) <=== getBBNSharedPtr(field_name);
82 * // DO OPERATOR WORK
83 * getNSharedPtr(base) ==> getBBNSharedPtr(field_name);
84 *
85 * @param field_name
86 * @param base
87 * @return MoFEMErrorCode
88 */
89 virtual MoFEMErrorCode baseSwap(const std::string &field_name,
90 const FieldApproximationBase base);
91
92 friend std::ostream &operator<<(std::ostream &os, const EntitiesFieldData &e);
93
94protected:
96};
97
98/** \brief this class derive data form other data structure
99 * \ingroup mofem_forces_and_sources_user_data_operators
100 *
101 *
102 * It behaves like normal data structure it is used to share base functions with
103 * other data structures. Dofs values, approx. order and
104 * indices are not shared.
105 *
106 * Shape functions, senses are shared with other data structure.
107 *
108 */
110
111 struct DerivedEntData;
112
114 const boost::shared_ptr<EntitiesFieldData> &data_ptr);
116
117private:
118 const boost::shared_ptr<EntitiesFieldData> dataPtr;
119};
120
121/** \brief Data on single entity (This is passed as argument to
122 * DataOperator::doWork) \ingroup mofem_forces_and_sources_user_data_operators
123 * \nosubgrouping
124 *
125 * \todo Hdiv and Hcurl functions should be accessed through common interface.
126 */
128
136 };
137
138 /** \name Constructor and destructor */
139
140 /**@{*/
141
142 EntData(const bool allocate_base_matrices = true);
143 virtual ~EntData() = default;
144
145 /**@}*/
146
147 /** \name Sense, order and indices */
148
149 /**@{*/
150
151 /// \brief get entity sense, need to calculate base functions with
152 /// conforming approximation fields
153 virtual int getSense() const;
154
155 /// \brief get approximation order
156 inline ApproximationOrder getOrder() const;
157
158 /// \brief Get global indices of dofs on entity
159 inline const VectorInt &getIndices() const;
160
161 /// \brief get global indices of dofs on entity up to given order
163
164 /// \brief get local indices of dofs on entity
165 inline const VectorInt &getLocalIndices() const;
166
167 /// \brief get local indices of dofs on entity up to given order
169
170 inline int &getSense();
171
173
174 inline VectorInt &getIndices();
175
176 inline VectorInt &getLocalIndices();
177
178 /**@}*/
179
180 /** \name Data on entity */
181
182 /**@{*/
183
184 /// \brief get dofs values
185 inline const VectorDouble &getFieldData() const;
186
187 /// \brief get dofs values up to given order
188 inline const VectorAdaptor getFieldDataUpToOrder(int order);
189
190 /// \brief get dofs data stature FEDofEntity
191 inline const VectorDofs &getFieldDofs() const;
192
193 /// \brief get dofs data stature FEDofEntity
194 inline VectorDouble &getFieldData();
195
196 /// \brief get field entities
197 inline const VectorFieldEntities &getFieldEntities() const;
198
199 /// \brief get field entities
201
202 //// \brief get entity bit ref level
204
205 /**
206 * @brief Return FTensor of rank 1, i.e. vector from filed data coeffinects
207 *
208 * \code
209 * auto t_vec = data.getFTensor1FieldData<3>();
210 * \endcode
211 *
212 * @tparam Tensor_Dim size of vector
213 * @return FTensor::Tensor1<FTensor::PackPtr<double *, Tensor_Dim>,
214 * Tensor_Dim>
215 */
216 template <int Tensor_Dim>
219
220 /**
221 * @brief Return FTensor rank 2, i.e. matrix from filed data coeffinects
222 *
223 * \code
224 * auto t_mat = data.getFTensor2FieldData<3,3>();
225 * \endcode
226 *
227 * @tparam Tensor_Dim0
228 * @tparam Tensor_Dim1
229 * @return FTensor::Tensor2<FTensor::PackPtr<double *, Tensor_Dim0 *
230 * Tensor_Dim1>, Tensor_Dim0, Tensor_Dim1>
231 */
232 template <int Tensor_Dim0, int Tensor_Dim1>
234 Tensor_Dim0, Tensor_Dim1>
236
237 /**
238 * @brief Return symmetric FTensor rank 2, i.e. matrix from filed data
239 * coeffinects
240 *
241 * \code
242 * auto t_mat = data.getFTensor2SymmetricFieldData<3>();
243 * \endcode
244 *
245 * @tparam Tensor_Dim dimension of the tensor
246 * @return FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, (Tensor_Dim
247 * * (Tensor_Dim + 1)) / 2>, Tensor_Dim>
248 */
249 template <int Tensor_Dim>
251 FTensor::PackPtr<double *, (Tensor_Dim * (Tensor_Dim + 1)) / 2>,
252 Tensor_Dim>
254
255 /**
256 * @brief Resturn scalar files as a FTensor of rank 0
257 *
258 * @return FTensor::Tensor0<FTensor::PackPtr<double *,1> >
259 */
261
262 inline VectorDofs &getFieldDofs();
263
264 /**@}*/
265
266 /** \name Base and space */
267
268 /**@{*/
269
270 /**
271 * \brief Get approximation base
272 * @return Approximation base
273 */
275
276 /**
277 * \brief Get field space
278 * @return Field space
279 */
280 inline FieldSpace &getSpace();
281
282 /**
283 * Get shared pointer to base base functions
284 */
285 virtual boost::shared_ptr<MatrixDouble> &
287 const BaseDerivatives direvatie);
288
289 /**
290 * Get shared pointer to base base functions
291 */
292 virtual boost::shared_ptr<MatrixDouble> &
294
295 /**
296 * Get shared pointer to derivatives of base base functions
297 */
298 virtual boost::shared_ptr<MatrixDouble> &
300
301 /**@}*/
302
303 /** \name Get base functions for H1/L2 */
304
305 /**@{*/
306
307 /** \brief get base functions
308 * this return matrix (nb. of rows is equal to nb. of Gauss pts, nb. of
309 * columns is equal to number of base functions on this entity.
310 *
311 * \note Note that for vectorial base, like Hdiv or Hcurl, in columns are
312 * vectorial base functions. For tonsorial would be tonsorial base
313 * functions. Interpretation depends on type of base, scalar, vectorial or
314 * tonsorial and dimension fo problem.
315 *
316 */
317 inline MatrixDouble &getN(const FieldApproximationBase base);
318
319 /**
320 * @copydoc MoFEM::EntitiesFieldData::EntData::getN
321 */
322 inline MatrixDouble &getN(const std::string &field_name);
323
324 /**
325 * @copydoc MoFEM::EntitiesFieldData::EntData::getN
326 */
327 inline MatrixDouble &getN();
328
329 /** \brief get derivatives of base functions
330 *
331 * Matrix at rows has nb. of Gauss pts, at columns it has derivative of
332 * base functions. Columns are structured as follows, [ dN1/dx, dN1/dy,
333 * dN1/dz, dN2/dx, dN2/dy, dN2/dz, ... ]
334 *
335 * Scalar base functions:
336 * Note that base functions are calculated in file H1.c
337 * Above description not apply for derivatives of nodal functions, since
338 * derivative of nodal functions in case of simplexes, EDGES, TRIANGLES and
339 * TETS are constant. So that matrix rows represents nb. of base
340 * functions, columns are derivatives. Nb. of columns depend on element
341 * dimension, for EDGES is one, for TRIS is 2 and TETS is 3.
342 *
343 * \note Note that for node element this function make no sense.
344 *
345 * Tonsorial base functions:
346 * \note Note: In rows ale integration pts, columns are formatted that that
347 * components of vectors and then derivatives, for example row for given
348 * integration points is formatted in array
349 * \f[
350 * t_{0,0}, t_{1,0}, t_{1,0}, t_{0,1}, t_{1,1}, t_{1,1}, t_{0,2}, t_{1,2},
351 * t_{1,2} \f] where comma express derivative, i.e. \f$t_{2,1} =
352 * \frac{\partial t_2}{\partial \xi_1}\f$
353 *
354 */
355 inline MatrixDouble &getDiffN(const FieldApproximationBase base);
356
357 /**
358 * @brief Get base function derivative
359 *
360 * @param base base
361 * @param derivative derivative
362 * @return MatrixDouble&
363 */
364 inline MatrixDouble &getN(const FieldApproximationBase base,
365 const BaseDerivatives derivative);
366
367 /**
368 * @copydoc MoFEM::EntitiesFieldData::EntData::getDiffN
369 */
370 inline MatrixDouble &getDiffN(const std::string &field_name);
371
372 /**
373 * @copydoc MoFEM::EntitiesFieldData::EntData::getDiffN
374 */
375 inline MatrixDouble &getDiffN();
376
377 /**
378 * @brief Get base function derivative
379 *
380 * @param derivative
381 * @return MatrixDouble&
382 */
383 inline MatrixDouble &getN(const BaseDerivatives derivative);
384
385 /// \brief get base functions at Gauss pts
386 inline const VectorAdaptor getN(const FieldApproximationBase base,
387 const int gg);
388
389 /// \brief get base functions at Gauss pts
390 inline const VectorAdaptor getN(const int gg);
391
392 /** \brief get derivative of base functions at Gauss pts
393
394 * returned matrix on rows has base functions, in column its derivatives.
395 *
396 * \param base Approximation base
397 * \param gg Nb. of Gauss pts.
398 *
399 */
400 inline const MatrixAdaptor getDiffN(const FieldApproximationBase base,
401 const int gg);
402
403 /** \brief get derivative of base functions at Gauss pts
404
405 * returned matrix on rows has base functions, in column its derivatives.
406 *
407 * \param gg nb. of Gauss pts.
408 *
409 */
410 inline const MatrixAdaptor getDiffN(const int gg);
411
412 /** \brief get base functions at Gauss pts
413
414 * Note that multi field element, two different field can have different
415 * approximation orders. Since we use hierarchical approximation basis,
416 * base functions are calculated once for element, using maximal
417 * approximation order on given entity.
418 *
419 * \param base Approximation base
420 * \param gg number of Gauss point
421 * \param nb_base_functions number of of base functions returned
422
423 */
424 inline const VectorAdaptor getN(const FieldApproximationBase base,
425 const int gg, const int nb_base_functions);
426
427 /** \brief get base functions at Gauss pts
428
429 * Note that multi field element, two different field can have different
430 * approximation orders. Since we use hierarchical approximation basis,
431 * base functions are calculated once for element, using maximal
432 * approximation order on given entity.
433 *
434 * \param gg number of Gauss point
435 * \param nb_base_functions number of of base functions returned
436
437 */
438 inline const VectorAdaptor getN(const int gg, const int nb_base_functions);
439
440 /** \brief get derivatives of base functions at Gauss pts
441 *
442 * Note that multi field element, two different field can have different
443 * approximation orders. Since we use hierarchical approximation basis,
444 * base functions are calculated once for element, using maximal
445 * approximation order on given entity.
446 *
447 * \param base Approximation base
448 * \param gg nb. of Gauss point
449 * \param nb_base_functions number of of base functions
450 *
451 */
452 inline const MatrixAdaptor getDiffN(const FieldApproximationBase base,
453 const int gg,
454 const int nb_base_functions);
455
456 /** \brief get derivatives of base functions at Gauss pts
457 *
458 * Note that multi field element, two different field can have different
459 * approximation orders. Since we use hierarchical approximation basis,
460 * base functions are calculated once for element, using maximal
461 * approximation order on given entity.
462 *
463 * \param gg nb. of Gauss point
464 * \param nb_base_functions number of of base functions
465 *
466 */
467 inline const MatrixAdaptor getDiffN(const int gg,
468 const int nb_base_functions);
469
470 /**@}*/
471
472 /** \name Get base functions for vectorial approximation basese, i.e.
473 * Hdiv/Hcurl */
474
475 /**@{*/
476
477 /** \brief get Hdiv of base functions at Gauss pts
478 *
479 * \param base Approximation base
480 * \param gg nb. of Gauss point
481 *
482 */
483 template <int DIM>
484 inline const MatrixAdaptor getVectorN(const FieldApproximationBase base,
485 const int gg);
486
487 /** \brief get Hdiv of base functions at Gauss pts
488 *
489 * \param gg nb. of Gauss point
490 * \param number of of base functions
491 *
492 */
493 template <int DIM> inline const MatrixAdaptor getVectorN(const int gg);
494
495 /** \brief get DiffHdiv of base functions at Gauss pts
496 *
497 * \param base Approximation base
498 * \param gg nb. of Gauss point
499 * \param number of of base functions
500 *
501 */
502 template <int DIM0, int DIM1>
504 const int gg);
505
506 /** \brief get DiffHdiv of base functions at Gauss pts
507 *
508 * \param gg nb. of Gauss point
509 * \param number of of base functions
510 *
511 */
512 template <int DIM0, int DIM1>
513 inline const MatrixAdaptor getVectorDiffN(const int gg);
514
515 /** \brief get DiffHdiv of base functions at Gauss pts
516 *
517 * \param base Approximation base
518 * \param gg nb. of Gauss point
519 * \param number of of base functions
520 *
521 */
522 template <int DIM0, int DIM1>
524 const int dof, const int gg);
525
526 /** \brief get DiffHdiv of base functions at Gauss pts
527 *
528 * \param gg nb. of Gauss point
529 * \param number of of base functions
530 *
531 */
532 template <int DIM0, int DIM1>
533 inline const MatrixAdaptor getVectorDiffN(const int dof, const int gg);
534
535 /**@}*/
536
537 /** \name Get base functions with FTensor */
538
539 /**@{*/
540
541 /**
542 * \brief Get base function as Tensor0
543 *
544 * \param base
545 * \return Tensor0
546 *
547 */
550
551 /**
552 * \brief Get base function as Tensor0
553 *
554 * Return base functions for field base
555 *
556 * \return Tensor0
557 *
558 */
560
561 /**
562 * \brief Get base function as Tensor0 (Loop by integration points)
563 *
564 * \param base
565 * \param gg integration points
566 * \param bb base function
567 * \return Tensor0
568
569 Note that:
570 \code
571 t0 = data.getFTensor0N(base,bb);
572 ++t0
573 \endcode
574 Increment in above code will move pointer to base function in next
575 integration point.
576
577 *
578 */
580 getFTensor0N(const FieldApproximationBase base, const int gg, const int bb);
581
582 /**
583 * \brief Get base function as Tensor0 (Loop by integration points)
584 *
585 * Return base functions for field base
586 *
587 * \param bb base function
588 * \return Tensor0
589 *
590 */
592 getFTensor0N(const int gg, const int bb);
593
594 /**
595 * \brief Get derivatives of base functions
596 *
597 * For volume element like tetrahedral or prism,
598 * \code
599 * Tensor1<double*,3> diff_base = data.getFTensor1DiffN<3>();
600 * \endcode
601 *
602 * For face element like triangle or quad
603 * \code
604 * Tensor1<double*,2> diff_base = data.getFTensor1DiffN<2>();
605 * \endcode
606 *
607 * \param base functions
608 * \return Tensor rank 1 (vector)
609 *
610 */
611 template <int Tensor_Dim>
614
615 /**
616 * \brief Get derivatives of base functions
617 *
618 * For volume element like tetrahedral or prism,
619 * \code
620 * Tensor1<double*,3> diff_base = data.getFTensor1DiffN<3>();
621 * \endcode
622 *
623 * For face element like triangle or quad
624 * \code
625 * Tensor1<double*,2> diff_base = data.getFTensor1DiffN<2>();
626 * \endcode
627 *
628 * \return Tensor rank 1 (vector)
629 *
630 */
631 template <int Tensor_Dim>
634
635 /**
636 * \brief Get derivatives of base functions (Loop by integration points)
637 *
638 * For volume element like tetrahedral or prism,
639 * \code
640 * Tensor1<double*,3> diff_base = data.getFTensor1DiffN<3>(base,gg,bb);
641 * \endcode
642 * where bb is base function and gg is integration pt. Operator ++diff_base
643 * will move tensor pointer to next integration point.
644 *
645 * For face element like triangle or quad
646 * \code
647 * Tensor1<double*,2> diff_base = data.getFTensor1DiffN<2>(base,gg,bb);
648 * \endcode
649 *
650 * \return Tensor rank 1 (vector)
651 *
652 */
653 template <int Tensor_Dim>
655 getFTensor1DiffN(const FieldApproximationBase base, const int gg,
656 const int bb);
657
658 /**
659 * \brief Get derivatives of base functions (Loop by integration points)
660 *
661 * For volume element like tetrahedral or prism,
662 * \code
663 * Tensor1<double*,3> diff_base = data.getFTensor1DiffN<3>(gg,bb);
664 * \endcode
665 * where bb is base function and gg is integration pt. Operator ++diff_base
666 * will move tensor pointer to next base function.
667 *
668 * For face element like triangle or quad
669 * \code
670 * Tensor1<double*,2> diff_base = data.getFTensor1DiffN<2>(gg,bb);
671 * \endcode
672 *
673 * \return Tensor rank 1 (vector)
674 *
675 */
676 template <int Tensor_Dim>
678 getFTensor1DiffN(const int gg, const int bb);
679
680 /** \brief Get base functions for Hdiv/Hcurl spaces
681
682 \note You probably like to use getFTensor1N(), in typical use base is
683 set automatically based on base set to field.
684
685 * @param base Approximation base
686
687 Example:
688 \code
689 FTensor::Index<'i',3> i;
690 int nb_dofs = data.getFieldData().size();
691 auto t_n_hdiv = data.getFTensor1N<3>();
692 for(int gg = 0;gg!=nb_gauss_pts;gg++) {
693 int ll = 0;
694 for(;ll!=nb_dofs;ll++) {
695 double dot_product = t_n_hdiv(i)*t_n_hdiv(i);
696 ++t_n_hdiv;
697 }
698 for(;ll!=data.getVectorN().size2()/3;ll++) {
699 ++t_n_hdiv;
700 }
701 }
702 \endcode
703
704 */
705 template <int Tensor_Dim>
708
709 /** \brief Get base functions for Hdiv space
710
711 Example:
712 \code
713 FTensor::Index<'i',3> i;
714 int nb_dofs = data.getFieldData().size();
715 auto t_n_hdiv = data.getFTensor1N<3>();
716 for(int gg = 0;gg!=nb_gauss_pts;gg++) {
717 int ll = 0;
718 for(;ll!=nb_dofs;ll++) {
719 double dot_product = t_n_hdiv(i)*t_n_hdiv(i);
720 ++t_n_hdiv;
721 }
722 for(;ll!=data.getVectorN().size2()/3;ll++) {
723 ++t_n_hdiv;
724 }
725 }
726 \endcode
727
728 */
729 template <int Tensor_Dim> auto getFTensor1N();
730
731 /** \brief Get derivatives of base functions for Hdiv space
732 */
733 template <int Tensor_Dim0, int Tensor_Dim1>
735 Tensor_Dim0, Tensor_Dim1>
737
738 /** \brief Get derivatives of base functions for Hdiv space at integration
739 * pts
740 */
741 template <int Tensor_Dim0, int Tensor_Dim1>
743 Tensor_Dim0, Tensor_Dim1>
744 getFTensor2DiffN(FieldApproximationBase base, const int gg, const int bb);
745
746 /** \brief Get derivatives of base functions for Hdiv space
747 */
748 template <int Tensor_Dim0, int Tensor_Dim1>
750 Tensor_Dim0, Tensor_Dim1>
752 return getFTensor2DiffN<Tensor_Dim0, Tensor_Dim1>(bAse);
753 }
754
755 /** \brief Get derivatives of base functions for Hdiv space at integration
756 * pts
757 */
758 template <int Tensor_Dim0, int Tensor_Dim1>
760 Tensor_Dim0, Tensor_Dim1>
761 getFTensor2DiffN(const int gg, const int bb) {
762 return getFTensor2DiffN<Tensor_Dim0, Tensor_Dim1>(bAse, gg, bb);
763 }
764
765 /** \brief Get second derivatives of base functions for Hvec space
766 */
767 template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2>
770 Tensor_Dim0, Tensor_Dim1, Tensor_Dim2>
772
773 /** \brief Get second derivatives of base functions for Hvec space
774 */
775 template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2>
776 inline FTensor::Tensor3<
778 Tensor_Dim0, Tensor_Dim1, Tensor_Dim2>
780 return getFTensor3Diff2N<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2>(bAse);
781 }
782
783 /**
784 * \brief Get Hdiv base functions at integration point
785
786 \code
787 FTensor::Index<'i',3> i;
788 for(int gg = 0;gg!=nb_gauss_pts;gg++) {
789 auto t_base = data.getFTensor1N(base,gg,bb);
790 for(int bb = 0;bb!=nb_base_functions;bb++) {
791 auto dot = t_base(i)*t_base(i);
792 }
793 }
794 \endcode
795
796 */
797 template <int Tensor_Dim>
799 getFTensor1N(FieldApproximationBase base, const int gg, const int bb);
800
801 /**
802 * \brief Get Hdiv base functions at integration point
803
804 \code
805 FTensor::Index<'i',3> i;
806 for(int gg = 0;gg!=nb_gauss_pts;gg++) {
807 auto t_base = data.getFTensor1N(gg,0);
808 for(int bb = 0;bb!=nb_base_functions;bb++) {
809 double dot = t_base(i)*t_base(i);
810 }
811 }
812 \endcode
813
814 */
815 template <int Tensor_Dim>
816 inline auto getFTensor1N(const int gg, const int bb);
817
818 /** \brief Get base functions for Hdiv/Hcurl spaces
819
820 \note You probably like to use getFTensor1N(), in typical use base is
821 set automatically based on base set to field.
822
823 * @param base Approximation base
824
825 Example:
826 \code
827 FTensor::Index<'i',3> i;
828 FTensor::Index<'i',3> j;
829 int nb_dofs = data.getFieldData().size();
830 auto t_n_hdiv = data.getFTensor2N<3,3>();
831 for(int gg = 0;gg!=nb_gauss_pts;gg++) {
832 int ll = 0;
833 for(;ll!=nb_dofs;ll++) {
834 double dot_product = t_n_hdiv(i,j)*t_n_hdiv(i,j);
835 ++t_n_hdiv;
836 }
837 for(;ll!=data.getVectorN().size2()/3;ll++) {
838 ++t_n_hdiv;
839 }
840 }
841 \endcode
842
843 */
844 template <int Tensor_Dim0, int Tensor_Dim1>
846 Tensor_Dim0, Tensor_Dim1>
848
849 /** \brief Get base functions for Hdiv space
850
851 Example:
852 \code
853 FTensor::Index<'i',3> i;
854 FTensor::Index<'j',3> j;
855
856 int nb_dofs = data.getFieldData().size();
857 auto t_n_hdiv = data.getFTensor2N<3,3>();
858 for(int gg = 0;gg!=nb_gauss_pts;gg++) {
859 int ll = 0;
860 for(;ll!=nb_dofs;ll++) {
861 double dot_product = t_n_hdiv(i,j)*t_n_hdiv(i,j);
862 ++t_n_hdiv;
863 }
864 for(;ll!=data.getVectorN().size2()/3;ll++) {
865 ++t_n_hdiv;
866 }
867 }
868 \endcode
869
870 */
871 template <int Tensor_Dim0, int Tensor_Dim1> auto getFTensor2N();
872
873 /** \brief Get base functions for tensor Hdiv/Hcurl spaces
874
875 \note You probably like to use getFTensor2N(), in typical use base is
876 set automatically based on base set to field.
877
878 @param base Approximation base
879
880 Example:
881 \code
882 FTensor::Index<'i',3> i;
883 FTensor::Index<'j',3> i;
884 int nb_dofs = data.getFieldData().size();
885 for(int gg = 0;gg!=nb_gauss_pts;gg++) {
886 auto t_n_hdiv = data.getFTensor2N<3>(base,gg,bb);
887 int ll = 0;
888 for(;ll!=nb_dofs;ll++) {
889 double dot_product = t_n_hdiv(i,j)*t_n_hdiv(i,j);
890 ++t_n_hdiv;
891 }
892 for(;ll!=data.getVectorN().size2()/3;ll++) {
893 ++t_n_hdiv;
894 }
895 }
896 \endcode
897
898 */
899 template <int Tensor_Dim0, int Tensor_Dim1>
901 Tensor_Dim0, Tensor_Dim1>
902 getFTensor2N(FieldApproximationBase base, const int gg, const int bb);
903
904 /** \brief Get base functions for Hdiv space
905
906 Example:
907 \code
908 FTensor::Index<'i',3> i;
909 FTensor::Index<'j',3> j;
910 int nb_dofs = data.getFieldData().size();
911 for(int gg = 0;gg!=nb_gauss_pts;++gg) {
912 int ll = 0;
913 auto t_n_hdiv = data.getFTensor2N<3,3>(gg,0);
914 for(;ll!=nb_dofs;ll++) {
915 double dot_product = t_n_hdiv(i)*t_n_hdiv(i);
916 ++t_n_hdiv;
917 }
918 for(;ll!=data.getVectorN().size2()/3;ll++) {
919 ++t_n_hdiv;
920 }
921 }
922 \endcode
923
924 */
925 template <int Tensor_Dim0, int Tensor_Dim1>
926 auto getFTensor2N(const int gg, const int bb);
927
928 /**@}*/
929
930 /** \name Auxiliary functions */
931
932 /**@{*/
933
934 friend std::ostream &operator<<(std::ostream &os,
936
937 /**
938 * Reset data associated with particular field name
939 * @return error code
940 */
942
943 /**@}*/
944
945 /** \name Bernstein-Bezier base only functions */
946
947 /**@{*/
948
949 /**
950 * @brief Get orders at the nodes
951 *
952 * @return VectorInt&
953 */
954 inline VectorInt &getBBNodeOrder();
955
956 /**
957 * @brief Get file BB indices
958 *
959 * @return MatrixInt&
960 */
962
963 virtual boost::shared_ptr<MatrixInt> &
964 getBBAlphaIndicesSharedPtr(const std::string &field_name);
965
966 /**
967 * Get shared pointer to BB base base functions
968 */
969 virtual boost::shared_ptr<MatrixDouble> &
970 getBBNSharedPtr(const std::string &field_name);
971
972 /**
973 * Get shared pointer to BB base base functions
974 */
975 virtual const boost::shared_ptr<MatrixDouble> &
976 getBBNSharedPtr(const std::string &field_name) const;
977
978 /**
979 * Get shared pointer to BB derivatives of base base functions
980 */
981 virtual boost::shared_ptr<MatrixDouble> &
982 getBBDiffNSharedPtr(const std::string &field_name);
983
984 /**
985 * Get shared pointer to derivatives of BB base base functions
986 */
987 virtual const boost::shared_ptr<MatrixDouble> &
988 getBBDiffNSharedPtr(const std::string &field_name) const;
989
990 virtual std::map<std::string, boost::shared_ptr<MatrixInt>> &
992
993 /**
994 * @brief get hash map of base function for BB base, key is a field name
995 *
996 * @return std::map<std::string, boost::shared_ptr<MatrixDouble>>&
997 */
998 virtual std::map<std::string, boost::shared_ptr<MatrixDouble>> &getBBNMap();
999
1000 /**
1001 * @brief get hash map of direvarives base function for BB base, key is a
1002 * field name
1003 *
1004 * @return std::map<std::string, boost::shared_ptr<MatrixDouble>>&
1005 */
1006 virtual std::map<std::string, boost::shared_ptr<MatrixDouble>> &
1007 getBBDiffNMap();
1008
1009 /**
1010 * @brief get ALpha indices for BB base by order
1011 *
1012 * @param o approximation order
1013 * @return boost::shared_ptr<MatrixInt>&
1014 */
1015 virtual boost::shared_ptr<MatrixInt> &
1017
1018 /**
1019 * @brief get BB base by order
1020 *
1021 * @param o
1022 * @return boost::shared_ptr<MatrixDouble>&
1023 */
1024 virtual boost::shared_ptr<MatrixDouble> &
1025 getBBNByOrderSharedPtr(const size_t o);
1026
1027 /**
1028 * @brief get BB base derivative by order
1029 *
1030 * @param o
1031 * @return boost::shared_ptr<MatrixDouble>&
1032 */
1033 virtual boost::shared_ptr<MatrixDouble> &
1034 getBBDiffNByOrderSharedPtr(const size_t o);
1035
1036 static constexpr size_t MaxBernsteinBezierOrder = BITFEID_SIZE;
1037
1038 virtual std::array<boost::shared_ptr<MatrixInt>, MaxBernsteinBezierOrder> &
1040
1041 virtual std::array<boost::shared_ptr<MatrixDouble>, MaxBernsteinBezierOrder> &
1043
1044 virtual std::array<boost::shared_ptr<MatrixDouble>, MaxBernsteinBezierOrder> &
1046
1047 /**
1048 * @brief Swap bases functions
1049 *
1050 * Some base are not hierarchical and depene on approximation order. Such case
1051 * demand special handling, that appropiate base order is set depending on
1052 * field, such that is accessible in operator.
1053 *
1054 * @note Base is not swap on meshsets
1055 *
1056 * @param field_name
1057 * @param base
1058 * @return MoFEMErrorCode
1059 */
1060 virtual MoFEMErrorCode baseSwap(const std::string &field_name,
1061 const FieldApproximationBase base);
1062
1063 /**@}*/
1064
1065protected:
1066 int sEnse; ///< Entity sense (orientation)
1067 ApproximationOrder oRder; ///< Entity order
1068 FieldSpace sPace; ///< Entity space
1069 FieldApproximationBase bAse; ///< Field approximation base
1070 BitRefLevel entDataBitRefLevel; ///< Bit ref level in entity
1071 VectorInt iNdices; ///< Global indices on entity
1072 VectorInt localIndices; ///< Local indices on entity
1073 VectorDofs dOfs; ///< DoFs on entity
1075 VectorDouble fieldData; ///< Field data on entity
1076
1077 std::array<std::array<boost::shared_ptr<MatrixDouble>, LASTBASE>,
1080
1081 std::array<boost::shared_ptr<MatrixDouble>, LASTBASE> &N; ///< Base functions
1082 std::array<boost::shared_ptr<MatrixDouble>, LASTBASE>
1083 &diffN; ///< Derivatives of base functions
1084
1085 std::string bbFieldName; ///< field name
1086 VectorInt bbNodeOrder; ///< order of nodes
1087 std::map<std::string, boost::shared_ptr<MatrixDouble>> bbN;
1088 std::map<std::string, boost::shared_ptr<MatrixDouble>> bbDiffN;
1089 std::map<std::string, boost::shared_ptr<MatrixInt>>
1090 bbAlphaIndices; ///< Indices for Bernstein-Bezier (BB) base
1091
1092 std::array<boost::shared_ptr<MatrixDouble>, MaxBernsteinBezierOrder>
1093 bbNByOrder; ///< BB base functions by order
1094 std::array<boost::shared_ptr<MatrixDouble>, MaxBernsteinBezierOrder>
1095 bbDiffNByOrder; ///< BB base functions direvatives by order
1096 std::array<boost::shared_ptr<MatrixInt>, MaxBernsteinBezierOrder>
1097 bbAlphaIndicesByOrder; ///< BB alpha indices by order
1098
1099protected:
1100 /**
1101 * @brief Used by Bernstein base to keep temporally pointer
1102 *
1103 * @copydoc MoFEM::EntitiesFieldData::baseSwap
1104 */
1105 boost::shared_ptr<MatrixDouble> swapBaseNPtr;
1106
1107 /**
1108 * @brief Used by Bernstein base to keep temporally pointer
1109 *
1110 * @copydoc MoFEM::EntitiesFieldData::baseSwap
1111 */
1112 boost::shared_ptr<MatrixDouble> swapBaseDiffNPtr;
1113
1114 friend struct OpAddParentEntData;
1115};
1116
1118
1119/** \brief Derived ata on single entity (This is passed as argument to
1120 * DataOperator::doWork) \ingroup mofem_forces_and_sources_user_data_operators
1121 * \nosubgrouping
1122 *
1123 * DerivedEntData share part information with EntData except infomation about
1124 * base functions.
1125 *
1126 */
1129
1131 const boost::shared_ptr<EntitiesFieldData::EntData> &ent_data_ptr);
1132
1133 int getSense() const;
1134
1135 //// \brief get entity bit ref level
1137
1138 boost::shared_ptr<MatrixDouble> &
1140 const BaseDerivatives derivative);
1141
1142 boost::shared_ptr<MatrixDouble> &
1144
1145 boost::shared_ptr<MatrixDouble> &
1147
1148 const boost::shared_ptr<MatrixDouble> &
1149 getNSharedPtr(const FieldApproximationBase base) const;
1150
1151 const boost::shared_ptr<MatrixDouble> &
1152 getDiffNSharedPtr(const FieldApproximationBase base) const;
1153
1154 inline boost::shared_ptr<MatrixDouble> &
1156
1157 inline boost::shared_ptr<MatrixDouble> &
1159
1160 boost::shared_ptr<MatrixInt> &
1161 getBBAlphaIndicesSharedPtr(const std::string &field_name);
1162
1163 /**
1164 * Get shared pointer to BB base base functions
1165 */
1166 boost::shared_ptr<MatrixDouble> &
1167 getBBNSharedPtr(const std::string &field_name);
1168
1169 /**
1170 * Get shared pointer to BB base base functions
1171 */
1172 const boost::shared_ptr<MatrixDouble> &
1173 getBBNSharedPtr(const std::string &field_name) const;
1174
1175 /**
1176 * Get shared pointer to BB derivatives of base base functions
1177 */
1178 boost::shared_ptr<MatrixDouble> &
1179 getBBDiffNSharedPtr(const std::string &field_name);
1180
1181 /**
1182 * Get shared pointer to derivatives of BB base base functions
1183 */
1184 const boost::shared_ptr<MatrixDouble> &
1185 getBBDiffNSharedPtr(const std::string &field_name) const;
1186
1187 /**
1188 * @copydoc MoFEM::EntitiesFieldData::EntData::swapBaseNPtr
1189 */
1190 MoFEMErrorCode baseSwap(const std::string &field_name,
1191 const FieldApproximationBase base);
1192
1193protected:
1194 const boost::shared_ptr<EntitiesFieldData::EntData> entDataPtr;
1195};
1196
1198 return oRder;
1199}
1200
1202 return iNdices;
1203}
1204
1205const VectorIntAdaptor
1207 unsigned int size = 0;
1208 if (auto dof = dOfs[0]) {
1209 size = dof->getOrderNbDofs(order) * dof->getNbOfCoeffs();
1210 size = size < iNdices.size() ? size : iNdices.size();
1211 }
1212 int *data = &*iNdices.data().begin();
1213 return VectorIntAdaptor(size, ublas::shallow_array_adaptor<int>(size, data));
1214}
1215
1217 return localIndices;
1218}
1219
1220const VectorIntAdaptor
1222 unsigned int size = 0;
1223 if (auto dof = dOfs[0]) {
1224 size = dof->getOrderNbDofs(order) * dof->getNbOfCoeffs();
1225 size = size < localIndices.size() ? size : localIndices.size();
1226 }
1227 int *data = &*localIndices.data().begin();
1228 return VectorIntAdaptor(size, ublas::shallow_array_adaptor<int>(size, data));
1229}
1230
1232
1234
1236
1238 return localIndices;
1239}
1240
1242 return fieldData;
1243}
1244
1245const VectorAdaptor
1247 unsigned int size = 0;
1248 if (auto dof = dOfs[0]) {
1249 size = dof->getOrderNbDofs(order) * dof->getNbOfCoeffs();
1250 size = size < fieldData.size() ? size : fieldData.size();
1251 }
1252 double *data = &*fieldData.data().begin();
1253 return getVectorAdaptor(data, size);
1254}
1255
1257 return dOfs;
1258}
1259
1261
1263
1265 return fieldEntities;
1266}
1267
1268const VectorFieldEntities &
1270 return fieldEntities;
1271}
1272
1273template <int Tensor_Dim>
1276 std::stringstream s;
1277 s << "Not implemented for this dimension dim = " << Tensor_Dim;
1278 THROW_MESSAGE(s.str());
1279}
1280
1281template <int Tensor_Dim0, int Tensor_Dim1>
1283 Tensor_Dim0, Tensor_Dim1>
1285 std::stringstream s;
1286 s << "Not implemented for this dimension dim0 = " << Tensor_Dim0;
1287 s << " and dim1 " << Tensor_Dim1;
1288 THROW_MESSAGE(s.str());
1289}
1290
1291template <int Tensor_Dim>
1293 FTensor::PackPtr<double *, (Tensor_Dim * (Tensor_Dim + 1)) / 2>, Tensor_Dim>
1295 std::stringstream s;
1296 s << "Not implemented for this dimension dim = " << Tensor_Dim;
1297 THROW_MESSAGE(s.str());
1298}
1299
1301
1303
1306 return *(getNSharedPtr(base));
1307}
1308
1310 return *(getBBNSharedPtr(field_name));
1311}
1312
1314
1317 return *(getDiffNSharedPtr(base));
1318}
1319
1322 const BaseDerivatives derivative) {
1323#ifndef NDEBUG
1324 if (!getNSharedPtr(base, derivative)) {
1325 MOFEM_LOG_C("SELF", Sev::error,
1326 "Ptr to base %s functions derivative %d is null",
1327 ApproximationBaseNames[base], derivative);
1328 THROW_MESSAGE("Null pointer");
1329 }
1330#endif
1331 return *(getNSharedPtr(base, derivative));
1332}
1333
1336 return *(getBBDiffNSharedPtr(field_name));
1337}
1338
1340
1343 return getN(bAse, derivative);
1344}
1345
1346const VectorAdaptor
1348 const int gg) {
1349 int size = getN(base).size2();
1350 double *data = &getN(base)(gg, 0);
1351 return VectorAdaptor(size, ublas::shallow_array_adaptor<double>(size, data));
1352}
1353
1355 return getN(bAse, gg);
1356}
1357
1358const MatrixAdaptor
1360 const int gg) {
1361 // FIXME: That is bug, it will not work if number of integration pts is
1362 // equal to number of nodes on entity. User who not implementing low
1363 // level DataOperator will not experience this.
1364 if (getN(base).size1() == getDiffN(base).size1()) {
1365 int size = getN(base).size2();
1366 int dim = getDiffN(base).size2() / size;
1367 double *data = &getDiffN(base)(gg, 0);
1368 return MatrixAdaptor(
1369 getN(base).size2(), dim,
1370 ublas::shallow_array_adaptor<double>(getDiffN(base).size2(), data));
1371 } else {
1372 // in some cases, f.e. for derivatives of nodal base functions at only
1373 // one gauss point is needed
1374 return MatrixAdaptor(
1375 getN(base).size1(), getN(base).size2(),
1376 ublas::shallow_array_adaptor<double>(getDiffN(base).data().size(),
1377 &getDiffN(base).data()[0]));
1378 }
1379}
1380
1382 return getDiffN(bAse, gg);
1383}
1384
1385const VectorAdaptor
1387 const int gg, const int nb_base_functions) {
1388 (void)getN()(gg, nb_base_functions -
1389 1); // throw error if nb_base_functions is to big
1390 double *data = &getN(base)(gg, 0);
1391 return VectorAdaptor(nb_base_functions, ublas::shallow_array_adaptor<double>(
1392 nb_base_functions, data));
1393}
1394
1395const VectorAdaptor
1396EntitiesFieldData::EntData::getN(const int gg, const int nb_base_functions) {
1397 return getN(bAse, gg, nb_base_functions);
1398}
1399
1400const MatrixAdaptor
1402 const int gg,
1403 const int nb_base_functions) {
1404 // FIXME: That is bug, it will not work if number of integration pts is
1405 // equal to number of nodes on entity. User who not implementing low
1406 // level DataOperator will not experience this.
1407 if (getN(base).size1() == getDiffN(base).size1()) {
1408 (void)getN(base)(gg,
1409 nb_base_functions -
1410 1); // throw error if nb_base_functions is to big
1411 int dim = getDiffN(base).size2() / getN(base).size2();
1412 double *data = &getDiffN(base)(gg, 0);
1413 return MatrixAdaptor(
1414 nb_base_functions, dim,
1415 ublas::shallow_array_adaptor<double>(dim * nb_base_functions, data));
1416 } else {
1417 // in some cases, f.e. for derivatives of nodal base functions only one
1418 // gauss point is needed
1419 return MatrixAdaptor(
1420 getN(base).size1(), getN(base).size2(),
1421 ublas::shallow_array_adaptor<double>(getDiffN(base).data().size(),
1422 &getDiffN(base).data()[0]));
1423 }
1424}
1425
1426const MatrixAdaptor
1428 const int nb_base_functions) {
1429 return getDiffN(bAse, gg, nb_base_functions);
1430}
1431
1432template <int DIM>
1433const MatrixAdaptor
1435 const int gg) {
1436 if (PetscUnlikely(getN(base).size2() % DIM)) {
1437 THROW_MESSAGE("Wrong dimension");
1438 }
1439
1440 const int nb_base_functions = getN(base).size2() / DIM;
1441 double *data = &getN(base)(gg, 0);
1442 return MatrixAdaptor(
1443 nb_base_functions, DIM,
1444 ublas::shallow_array_adaptor<double>(DIM * nb_base_functions, data));
1445}
1446
1447template <int DIM>
1449 return getVectorN<DIM>(bAse, gg);
1450}
1451
1452template <int DIM0, int DIM1>
1453const MatrixAdaptor
1455 const int gg) {
1456 if (PetscUnlikely(getDiffN(base).size2() % (DIM0 * DIM1))) {
1457 THROW_MESSAGE("Wrong dimension");
1458 }
1459
1460 const int nb_base_functions = getN(base).size2() / (DIM0 * DIM1);
1461 double *data = &getN(base)(gg, 0);
1462 return MatrixAdaptor(nb_base_functions, DIM0 * DIM1,
1463 ublas::shallow_array_adaptor<double>(
1464 DIM0 * DIM1 * nb_base_functions, data));
1465}
1466
1467template <int DIM0, int DIM1>
1469 return getVectorDiffN<DIM0, DIM1>(bAse, gg);
1470}
1471
1472template <int DIM0, int DIM1>
1473const MatrixAdaptor
1475 const int dof, const int gg) {
1476 double *data =
1477 &EntitiesFieldData::EntData::getDiffN(base)(gg, DIM0 * DIM1 * dof);
1478 return MatrixAdaptor(DIM0, DIM1,
1479 ublas::shallow_array_adaptor<double>(DIM0 * DIM1, data));
1480}
1481
1482template <int DIM0, int DIM1>
1484 const int gg) {
1485 return getVectorDiffN<DIM0, DIM1>(bAse, dof, gg);
1486}
1487
1490 double *ptr = &*getN(base).data().begin();
1492};
1493
1496 return getFTensor0N(bAse);
1497};
1498
1501 const int gg, const int bb) {
1502 double *ptr = &getN(base)(gg, bb);
1504};
1505
1507EntitiesFieldData::EntData::getFTensor0N(const int gg, const int bb) {
1508 return getFTensor0N(bAse, gg, bb);
1509};
1510
1511template <int Tensor_Dim> auto EntitiesFieldData::EntData::getFTensor1N() {
1512 return getFTensor1N<Tensor_Dim>(bAse);
1513}
1514
1515template <int Tensor_Dim>
1516auto EntitiesFieldData::EntData::getFTensor1N(const int gg, const int bb) {
1517 return getFTensor1N<Tensor_Dim>(bAse, gg, bb);
1518}
1519
1520template <int Tensor_Dim0, int Tensor_Dim1>
1522 return getFTensor2N<Tensor_Dim0, Tensor_Dim1>(bAse);
1523}
1524
1525template <int Tensor_Dim0, int Tensor_Dim1>
1526auto EntitiesFieldData::EntData::getFTensor2N(const int gg, const int bb) {
1527 return getFTensor2N<Tensor_Dim0, Tensor_Dim1>(bAse, gg, bb);
1528}
1529
1530/** \name Bernstein-Bezier base only functions */
1531
1532/**@{*/
1533
1535
1537 return *getBBAlphaIndicesSharedPtr(bbFieldName);
1538}
1539
1540/**@}*/
1541
1542/** \name DerivedEntData */
1543
1544/**@{*/
1545
1546boost::shared_ptr<MatrixDouble> &
1548 const FieldApproximationBase base) {
1549 return N[base];
1550}
1551
1552boost::shared_ptr<MatrixDouble> &
1554 const FieldApproximationBase base) {
1555 return diffN[base];
1556}
1557
1558/**@}*/
1559
1560/**
1561 * @brief Assemble PETSc vector
1562 *
1563 * Function extract indices from entity data and assemble vector
1564 *
1565 * <a
1566 * href=https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecSetValues.html>See
1567 * PETSc documentation</a>
1568 *
1569 * @param V
1570 * @param data
1571 * @param ptr
1572 * @param iora
1573 * @return MoFEMErrorCode
1574 */
1575template <typename T = EntityStorage>
1577 const EntitiesFieldData::EntData &data,
1578 const double *ptr, InsertMode iora) {
1579 static_assert(!std::is_same<T, T>::value,
1580 "VecSetValues value for this data storage is not implemented");
1581 return MOFEM_NOT_IMPLEMENTED;
1582}
1583
1584template <>
1587 const double *ptr, InsertMode iora) {
1588 return VecSetValues(V, data.getIndices().size(), &*data.getIndices().begin(),
1589 ptr, iora);
1590}
1591
1592/**
1593 * @brief Assemble PETSc matrix
1594 *
1595 * Function extract indices from entity data and assemble vector
1596 *
1597 * <a
1598 * href=https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatSetValues.html>See
1599 * PETSc documentation</a>
1600 *
1601 * @param M
1602 * @param row_data
1603 * @param col_data
1604 * @param ptr
1605 * @param iora
1606 * @return MoFEMErrorCode
1607 */
1608template <typename T = EntityStorage>
1610 const EntitiesFieldData::EntData &row_data,
1611 const EntitiesFieldData::EntData &col_data,
1612 const double *ptr, InsertMode iora) {
1613 static_assert(!std::is_same<T, T>::value,
1614 "MatSetValues value for this data storage is not implemented");
1615 return MOFEM_NOT_IMPLEMENTED;
1616}
1617
1618template <>
1621 const EntitiesFieldData::EntData &col_data,
1622 const double *ptr, InsertMode iora) {
1623 return MatSetValues(
1624 M, row_data.getIndices().size(), &*row_data.getIndices().begin(),
1625 col_data.getIndices().size(), &*col_data.getIndices().begin(), ptr, iora);
1626}
1627
1628/** \name Specializations for tensor base function */
1629
1630/**@{*/
1631
1632template <>
1634EntitiesFieldData::EntData::getFTensor1N<3>(FieldApproximationBase base);
1635
1636template <>
1638EntitiesFieldData::EntData::getFTensor1N<3>(FieldApproximationBase base,
1639 const int gg, const int bb);
1640
1641template <>
1643EntitiesFieldData::EntData::getFTensor2N<3, 3>(FieldApproximationBase base);
1644
1645/**@}*/
1646
1647/** \name Specializations for direcatives of base functions */
1648
1649/**@{*/
1650
1651template <>
1653EntitiesFieldData::EntData::getFTensor1DiffN<3>(
1654 const FieldApproximationBase base);
1655template <>
1657EntitiesFieldData::EntData::getFTensor1DiffN<3>();
1658
1659template <>
1661EntitiesFieldData::EntData::getFTensor1DiffN<2>(
1662 const FieldApproximationBase base);
1663template <>
1665EntitiesFieldData::EntData::getFTensor1DiffN<2>();
1666
1667template <>
1669EntitiesFieldData::EntData::getFTensor2DiffN<3, 2>(FieldApproximationBase base);
1670template <>
1672EntitiesFieldData::EntData::getFTensor2DiffN<3, 2>(FieldApproximationBase base,
1673 const int gg, const int bb);
1674
1675template <>
1678
1679/**@}*/
1680
1681/** \name Specializations for field data */
1682
1683/**@{*/
1684
1685template <>
1687EntitiesFieldData::EntData::getFTensor1FieldData<3>();
1688
1689template <>
1691EntitiesFieldData::EntData::getFTensor1FieldData<2>();
1692
1693template <>
1695EntitiesFieldData::EntData::getFTensor1FieldData<1>();
1696
1697template <>
1699EntitiesFieldData::EntData::getFTensor2FieldData<3, 3>();
1700
1701template <>
1703EntitiesFieldData::EntData::getFTensor2SymmetricFieldData<3>();
1704
1705template <>
1707EntitiesFieldData::EntData::getFTensor2SymmetricFieldData<2>();
1708
1709/**@}*/
1710
1711/**
1712 * @deprecated Use EntitiesFieldData
1713 */
1715
1716/**
1717 * @deprecated Use DerivedEntitiesFieldData
1718 */
1720
1721} // namespace MoFEM
1722
1723#endif //__ENTITIES_FIELD_DATA_HPP__
1724
1725/**
1726 * \defgroup mofem_forces_and_sources_user_data_operators User data operator
1727 * data structures \ingroup
1728 *
1729 * \brief Users data structures and operator
1730 *
1731 * Data structures passed by argument to MoFEM::DataOperator::doWork and generic
1732 * user operators operating on those structures.
1733 *
1734 */
static Index< 'M', 3 > M
static Index< 'o', 3 > o
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:304
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ LASTBASE
Definition: definitions.h:69
FieldSpace
approximation spaces
Definition: definitions.h:82
@ LASTSPACE
FieldSpace in [ 0, LASTSPACE )
Definition: definitions.h:89
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
static const char *const ApproximationBaseNames[]
Definition: definitions.h:72
#define DEPRECATED
Definition: definitions.h:17
#define BITFEID_SIZE
max number of finite elements
Definition: definitions.h:221
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:561
const int dim
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
VectorShallowArrayAdaptor< int > VectorIntAdaptor
Definition: Types.hpp:116
int ApproximationOrder
Approximation on the entity.
Definition: Types.hpp:26
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:115
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
Definition: Types.hpp:132
UBlasMatrix< int > MatrixInt
Definition: Types.hpp:76
UBlasVector< int > VectorInt
Definition: Types.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:31
DEPRECATED typedef DerivedEntitiesFieldData DerivedDataForcesAndSourcesCore
ublas::vector< FieldEntity *, FieldEntAllocator > VectorFieldEntities
DEPRECATED typedef EntitiesFieldData DataForcesAndSourcesCore
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
ublas::unbounded_array< FieldEntity *, std::allocator< FieldEntity * > > FieldEntAllocator
MoFEMErrorCode VecSetValues< EntityStorage >(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
MoFEMErrorCode MatSetValues< EntityStorage >(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
ublas::unbounded_array< FEDofEntity *, std::allocator< FEDofEntity * > > DofsAllocator
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
constexpr auto field_name
const int N
Definition: speed_test.cpp:3
Derived ata on single entity (This is passed as argument to DataOperator::doWork)
boost::shared_ptr< MatrixDouble > & getBBNSharedPtr(const std::string &field_name)
boost::shared_ptr< MatrixDouble > & getDerivedDiffNSharedPtr(const FieldApproximationBase base)
const boost::shared_ptr< EntitiesFieldData::EntData > entDataPtr
MoFEMErrorCode baseSwap(const std::string &field_name, const FieldApproximationBase base)
Used by Bernstein base to keep temporally pointer.
int getSense() const
get entity sense, need to calculate base functions with conforming approximation fields
boost::shared_ptr< MatrixInt > & getBBAlphaIndicesSharedPtr(const std::string &field_name)
boost::shared_ptr< MatrixDouble > & getBBDiffNSharedPtr(const std::string &field_name)
boost::shared_ptr< MatrixDouble > & getDiffNSharedPtr(const FieldApproximationBase base)
boost::shared_ptr< MatrixDouble > & getNSharedPtr(const FieldApproximationBase base, const BaseDerivatives derivative)
boost::shared_ptr< MatrixDouble > & getDerivedNSharedPtr(const FieldApproximationBase base)
this class derive data form other data structure
MoFEMErrorCode setElementType(const EntityType type)
const boost::shared_ptr< EntitiesFieldData > dataPtr
Data on single entity (This is passed as argument to DataOperator::doWork)
std::array< boost::shared_ptr< MatrixDouble >, MaxBernsteinBezierOrder > bbDiffNByOrder
BB base functions direvatives by order.
boost::shared_ptr< MatrixDouble > swapBaseDiffNPtr
Used by Bernstein base to keep temporally pointer.
virtual std::array< boost::shared_ptr< MatrixDouble >, MaxBernsteinBezierOrder > & getBBDiffNByOrderArray()
const VectorIntAdaptor getIndicesUpToOrder(int order)
get global indices of dofs on entity up to given order
virtual std::array< boost::shared_ptr< MatrixDouble >, MaxBernsteinBezierOrder > & getBBNByOrderArray()
int sEnse
Entity sense (orientation)
auto getFTensor2N()
Get base functions for Hdiv space.
VectorDouble fieldData
Field data on entity.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const int gg, const int bb)
Get derivatives of base functions (Loop by integration points)
MatrixDouble & getN()
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
FieldApproximationBase & getBase()
Get approximation base.
boost::shared_ptr< MatrixDouble > swapBaseNPtr
Used by Bernstein base to keep temporally pointer.
ApproximationOrder getOrder() const
get approximation order
std::map< std::string, boost::shared_ptr< MatrixInt > > bbAlphaIndices
Indices for Bernstein-Bezier (BB) base.
friend std::ostream & operator<<(std::ostream &os, const EntitiesFieldData::EntData &e)
VectorInt localIndices
Local indices on entity.
const MatrixAdaptor getVectorDiffN(FieldApproximationBase base, const int gg)
get DiffHdiv of base functions at Gauss pts
const VectorFieldEntities & getFieldEntities() const
get field entities
virtual std::array< boost::shared_ptr< MatrixInt >, MaxBernsteinBezierOrder > & getBBAlphaIndicesByOrderArray()
virtual std::map< std::string, boost::shared_ptr< MatrixInt > > & getBBAlphaIndicesMap()
const VectorIntAdaptor getLocalIndicesUpToOrder(int order)
get local indices of dofs on entity up to given order
const VectorAdaptor getFieldDataUpToOrder(int order)
get dofs values up to given order
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN()
Get derivatives of base functions.
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(const int gg, const int bb)
Get derivatives of base functions for Hdiv space at integration pts.
static constexpr size_t MaxBernsteinBezierOrder
VectorInt iNdices
Global indices on entity.
virtual boost::shared_ptr< MatrixInt > & getBBAlphaIndicesSharedPtr(const std::string &field_name)
const VectorDouble & getFieldData() const
get dofs values
virtual boost::shared_ptr< MatrixDouble > & getBBDiffNByOrderSharedPtr(const size_t o)
get BB base derivative by order
FTensor::Tensor3< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 *Tensor_Dim2 >, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > getFTensor3Diff2N(FieldApproximationBase base)
Get second derivatives of base functions for Hvec space.
virtual std::map< std::string, boost::shared_ptr< MatrixDouble > > & getBBNMap()
get hash map of base function for BB base, key is a field name
BitRefLevel entDataBitRefLevel
Bit ref level in entity.
const MatrixAdaptor getVectorN(const FieldApproximationBase base, const int gg)
get Hdiv of base functions at Gauss pts
virtual MoFEMErrorCode baseSwap(const std::string &field_name, const FieldApproximationBase base)
Swap bases functions.
const VectorInt & getLocalIndices() const
get local indices of dofs on entity
ApproximationOrder oRder
Entity order.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from filed data coeffinects.
virtual std::map< std::string, boost::shared_ptr< MatrixDouble > > & getBBDiffNMap()
get hash map of direvarives base function for BB base, key is a field name
virtual int getSense() const
get entity sense, need to calculate base functions with conforming approximation fields
VectorFieldEntities fieldEntities
Field entities.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N()
Get base function as Tensor0.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
MatrixInt & getBBAlphaIndices()
Get file BB indices.
std::array< boost::shared_ptr< MatrixInt >, MaxBernsteinBezierOrder > bbAlphaIndicesByOrder
BB alpha indices by order.
MatrixDouble & getDiffN()
get derivatives of base functions
VectorInt & getBBNodeOrder()
Get orders at the nodes.
virtual boost::shared_ptr< MatrixInt > & getBBAlphaIndicesByOrderSharedPtr(const size_t o)
get ALpha indices for BB base by order
std::array< boost::shared_ptr< MatrixDouble >, LASTBASE > & diffN
Derivatives of base functions.
std::map< std::string, boost::shared_ptr< MatrixDouble > > bbDiffN
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *,(Tensor_Dim *(Tensor_Dim+1))/2 >, Tensor_Dim > getFTensor2SymmetricFieldData()
Return symmetric FTensor rank 2, i.e. matrix from filed data coeffinects.
FTensor::Tensor3< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 *Tensor_Dim2 >, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > getFTensor3Diff2N()
Get second derivatives of base functions for Hvec space.
FieldSpace & getSpace()
Get field space.
auto getFTensor1N()
Get base functions for Hdiv space.
virtual boost::shared_ptr< MatrixDouble > & getBBNSharedPtr(const std::string &field_name)
std::array< boost::shared_ptr< MatrixDouble >, MaxBernsteinBezierOrder > bbNByOrder
BB base functions by order.
std::array< std::array< boost::shared_ptr< MatrixDouble >, LASTBASE >, LastDerivative > baseFunctionsAndBaseDerivatives
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN()
Get derivatives of base functions for Hdiv space.
std::map< std::string, boost::shared_ptr< MatrixDouble > > bbN
virtual boost::shared_ptr< MatrixDouble > & getDiffNSharedPtr(const FieldApproximationBase base)
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2FieldData()
Return FTensor rank 2, i.e. matrix from filed data coeffinects.
virtual boost::shared_ptr< MatrixDouble > & getBBNByOrderSharedPtr(const size_t o)
get BB base by order
virtual boost::shared_ptr< MatrixDouble > & getNSharedPtr(const FieldApproximationBase base, const BaseDerivatives direvatie)
const VectorInt & getIndices() const
Get global indices of dofs on entity.
FieldApproximationBase bAse
Field approximation base.
virtual boost::shared_ptr< MatrixDouble > & getBBDiffNSharedPtr(const std::string &field_name)
std::array< boost::shared_ptr< MatrixDouble >, LASTBASE > & N
Base functions.
virtual BitRefLevel & getEntDataBitRefLevel()
data structure for finite element entity
friend std::ostream & operator<<(std::ostream &os, const EntitiesFieldData &e)
MoFEMErrorCode resetFieldDependentData()
std::bitset< LASTSPACE > sPace
spaces on element
MatrixInt facesNodes
nodes on finite element faces
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
std::bitset< LASTBASE > bAse
bases on element
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
std::array< std::bitset< LASTBASE >, LASTSPACE > basesOnSpaces
base on spaces
virtual ~EntitiesFieldData()=default
virtual MoFEMErrorCode baseSwap(const std::string &field_name, const FieldApproximationBase base)
Swap approximation base.
virtual MoFEMErrorCode setElementType(const EntityType type)
MatrixInt facesNodesOrder
order of face nodes on element
std::array< std::bitset< LASTBASE >, MBMAXTYPE > basesOnEntities
bases on entity types
keeps information about indexed dofs for the finite element
Struct keeps handle to entity in the field.
Operator to project base functions from parent entity.