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