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