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