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