v0.14.0
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
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
203 virtual std::vector<BitRefLevel> &getEntDataBitRefLevel();
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 VectorInt iNdices; ///< Global indices on entity
1071 VectorInt localIndices; ///< Local indices on entity
1072 VectorDofs dOfs; ///< DoFs on entity
1074 VectorDouble fieldData; ///< Field data on entity
1075 std::vector<BitRefLevel> entDataBitRefLevel; ///< Bit ref level in 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
1136 std::vector<BitRefLevel> &getEntDataBitRefLevel();
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
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
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> &
1551
1552boost::shared_ptr<MatrixDouble> &
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 vector
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/Vec/VecSetValues.html>See
1599 * PETSc documentation</a>
1600 *
1601 * @param V
1602 * @param data
1603 * @param vec
1604 * @param iora
1605 * @return MoFEMErrorCode
1606 */
1607template <typename T = EntityStorage>
1609 const EntitiesFieldData::EntData &data,
1610 const VectorDouble &vec, InsertMode iora) {
1611 return VecSetValues<T>(V, data, &*vec.data().begin(), iora);
1612}
1613
1614/**
1615 * @brief Assemble PETSc matrix
1616 *
1617 * Function extract indices from entity data and assemble vector
1618 *
1619 * <a
1620 * href=https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatSetValues.html>See
1621 * PETSc documentation</a>
1622 *
1623 * @param M
1624 * @param row_data
1625 * @param col_data
1626 * @param ptr
1627 * @param iora
1628 * @return MoFEMErrorCode
1629 */
1630template <typename T = EntityStorage>
1632 const EntitiesFieldData::EntData &row_data,
1633 const EntitiesFieldData::EntData &col_data,
1634 const double *ptr, InsertMode iora) {
1635 static_assert(!std::is_same<T, T>::value,
1636 "MatSetValues value for this data storage is not implemented");
1637 return MOFEM_NOT_IMPLEMENTED;
1638}
1639
1640/**
1641 * @brief Assemble PETSc matrix
1642 *
1643 * Function extract indices from entity data and assemble vector
1644 *
1645 * <a
1646 * href=https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatSetValues.html>See
1647 * PETSc documentation</a>
1648 *
1649 * @param M
1650 * @param row_data
1651 * @param col_data
1652 * @param mat
1653 * @param iora
1654 * @return MoFEMErrorCode
1655 */
1656template <typename T = EntityStorage>
1658 const EntitiesFieldData::EntData &row_data,
1659 const EntitiesFieldData::EntData &col_data,
1660 const MatrixDouble &mat, InsertMode iora) {
1661 return MatSetValues<T>(M, row_data, col_data, &*mat.data().begin(), iora);
1662}
1663
1664template <>
1667 const EntitiesFieldData::EntData &col_data,
1668 const double *ptr, InsertMode iora) {
1669 return MatSetValues(
1670 M, row_data.getIndices().size(), &*row_data.getIndices().begin(),
1671 col_data.getIndices().size(), &*col_data.getIndices().begin(), ptr, iora);
1672}
1673
1674/** \name Specializations for tensor base function */
1675
1676/**@{*/
1677
1678template <>
1680EntitiesFieldData::EntData::getFTensor1N<3>(FieldApproximationBase base);
1681
1682template <>
1684EntitiesFieldData::EntData::getFTensor1N<3>(FieldApproximationBase base,
1685 const int gg, const int bb);
1686
1687template <>
1689EntitiesFieldData::EntData::getFTensor2N<3, 3>(FieldApproximationBase base);
1690
1691/**@}*/
1692
1693/** \name Specializations for direcatives of base functions */
1694
1695/**@{*/
1696
1697template <>
1699EntitiesFieldData::EntData::getFTensor1DiffN<3>(
1700 const FieldApproximationBase base);
1701template <>
1703EntitiesFieldData::EntData::getFTensor1DiffN<3>();
1704
1705template <>
1707EntitiesFieldData::EntData::getFTensor1DiffN<2>(
1708 const FieldApproximationBase base);
1709template <>
1711EntitiesFieldData::EntData::getFTensor1DiffN<2>();
1712
1713template <>
1715EntitiesFieldData::EntData::getFTensor2DiffN<3, 2>(FieldApproximationBase base);
1716template <>
1718EntitiesFieldData::EntData::getFTensor2DiffN<3, 2>(FieldApproximationBase base,
1719 const int gg, const int bb);
1720
1721template <>
1724
1725/**@}*/
1726
1727/** \name Specializations for field data */
1728
1729/**@{*/
1730
1731template <>
1733EntitiesFieldData::EntData::getFTensor1FieldData<3>();
1734
1735template <>
1737EntitiesFieldData::EntData::getFTensor1FieldData<2>();
1738
1739template <>
1741EntitiesFieldData::EntData::getFTensor1FieldData<1>();
1742
1743template <>
1745EntitiesFieldData::EntData::getFTensor2FieldData<1, 1>();
1746
1747template <>
1749EntitiesFieldData::EntData::getFTensor2FieldData<1, 2>();
1750
1751template <>
1753EntitiesFieldData::EntData::getFTensor2FieldData<1, 3>();
1754
1755template <>
1757EntitiesFieldData::EntData::getFTensor2FieldData<2, 2>();
1758
1759template <>
1761EntitiesFieldData::EntData::getFTensor2FieldData<3, 3>();
1762
1763template <>
1765EntitiesFieldData::EntData::getFTensor2SymmetricFieldData<3>();
1766
1767template <>
1769EntitiesFieldData::EntData::getFTensor2SymmetricFieldData<2>();
1770
1771/**@}*/
1772
1773/**
1774 * @deprecated Use EntitiesFieldData
1775 */
1777
1778/**
1779 * @deprecated Use DerivedEntitiesFieldData
1780 */
1782
1783} // namespace MoFEM
1784
1785#endif //__ENTITIES_FIELD_DATA_HPP__
1786
1787/**
1788 * \defgroup mofem_forces_and_sources_user_data_operators User data operator
1789 * data structures \ingroup
1790 *
1791 * \brief Users data structures and operator
1792 *
1793 * Data structures passed by argument to MoFEM::DataOperator::doWork and generic
1794 * user operators operating on those structures.
1795 *
1796 */
static Index< 'M', 3 > M
static Index< 'o', 3 > o
#define MOFEM_LOG_C(channel, severity, format,...)
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
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
constexpr int order
const int dim
constexpr int DIM1
Definition level_set.cpp:21
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorShallowArrayAdaptor< int > VectorIntAdaptor
Definition Types.hpp:116
int ApproximationOrder
Approximation on the entity.
Definition Types.hpp:26
UBlasMatrix< int > MatrixInt
Definition Types.hpp:76
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition Types.hpp:115
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
Definition Types.hpp:132
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
ublas::unbounded_array< FEDofEntity *, std::allocator< FEDofEntity * > > DofsAllocator
DEPRECATED typedef DerivedEntitiesFieldData DerivedDataForcesAndSourcesCore
DEPRECATED typedef EntitiesFieldData DataForcesAndSourcesCore
ublas::vector< FieldEntity *, FieldEntAllocator > VectorFieldEntities
EntitiesFieldData::EntData EntData
ublas::unbounded_array< FieldEntity *, std::allocator< FieldEntity * > > FieldEntAllocator
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
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)
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
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)
DerivedEntData(const boost::shared_ptr< EntitiesFieldData::EntData > &ent_data_ptr)
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
DerivedEntitiesFieldData(const boost::shared_ptr< EntitiesFieldData > &data_ptr)
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
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.
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.
std::vector< BitRefLevel > entDataBitRefLevel
Bit ref level in entity.
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
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 std::vector< BitRefLevel > & getEntDataBitRefLevel()
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::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.
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 to child.