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