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