HArD::Core2D
Hybrid Arbitrary Degree::Core 2D - Library to implement 2D schemes with edge and cell polynomials as unknowns
basis.hpp
Go to the documentation of this file.
1 // Core data structures and methods for 2D schemes based on polynomial unknowns in the elements and on the edges
2 //
3 // Provides:
4 // - Full and partial polynomial spaces on the element and edges
5 // - Generic routines to create quadrature nodes over cells and edges of the mesh
6 //
7 // Authors: Daniele Di Pietro (daniele.di-pietro@umontpellier.fr)
8 // Jerome Droniou (jerome.droniou@monash.edu)
9 //
10 
11 /*
12  *
13  * This library was developed around HHO methods, although some parts of it have a more
14  * general purpose. If you use this code or part of it in a scientific publication,
15  * please mention the following book as a reference for the underlying principles
16  * of HHO schemes:
17  *
18  * The Hybrid High-Order Method for Polytopal Meshes: Design, Analysis, and Applications.
19  * D. A. Di Pietro and J. Droniou. Modeling, Simulation and Applications, vol. 19.
20  * Springer International Publishing, 2020, xxxi + 525p. doi: 10.1007/978-3-030-37203-3.
21  * url: https://hal.archives-ouvertes.fr/hal-02151813.
22  *
23  */
24 
25 #ifndef BASIS_HPP
26 #define BASIS_HPP
27 
28 #include <boost/multi_array.hpp>
29 
30 #include <mesh.hpp>
31 //#include <edge.hpp>
32 #include <iostream>
33 
35 
36 #include <quadraturerule.hpp>
37 
38 namespace HArDCore2D
39 {
40 
53  constexpr int dimspace = 2;
54  typedef Eigen::Matrix2d MatrixRd;
55  typedef Eigen::Vector2d VectorRd;
56  typedef Eigen::Vector2i VectorZd;
57 
58  template <typename T>
59  using BasisQuad = boost::multi_array<T, 2>;
60 
61  template <typename T>
62  using FType = std::function<T(const VectorRd &)>;
63 
65  {
66  Scalar = 0,
67  Vector = 1,
68  Matrix = 2
69  };
70 
71 
72  //-----------------------------------------------------------------------------
73  // Powers of monomial basis
74  //-----------------------------------------------------------------------------
75 
77  template<typename GeometricSupport>
79  {
80  // Only specializations are relevant
81  };
82 
83  template<>
85  {
86  static std::vector<VectorZd> complete(size_t degree)
87  {
88  std::vector<VectorZd> powers;
89  powers.reserve( PolynomialSpaceDimension<Cell>::Poly(degree) );
90  for (size_t l = 0; l <= degree; l++)
91  {
92  for (size_t i = 0; i <= l; i++)
93  {
94  powers.push_back(VectorZd(i, l - i));
95  } // for i
96  } // for l
97 
98  return powers;
99  }
100  };
101 
102  //----------------------------------------------------------------------
103  //----------------------------------------------------------------------
104  // STRUCTURE FOR GRADIENT OF MATRIX-VALUED FUNCTIONS
105  //----------------------------------------------------------------------
106  //----------------------------------------------------------------------
107 
109  template<size_t N>
111  {
114  Eigen::Matrix<double, N, N> diff_x,
115  Eigen::Matrix<double, N, N> diff_y
116  ) : dx(diff_x),
117  dy(diff_y)
118  {
119  // do nothing
120  }
121 
123  MatrixGradient() : dx(Eigen::Matrix<double, N, N>::Zero()),
124  dy(Eigen::Matrix<double, N, N>::Zero())
125  {
126  // do nothing
127  }
128 
131  return MatrixGradient<N>(dx + G.dx, dy + G.dy);
132  }
133 
136  this->dx += G.dx;
137  this->dy += G.dy;
138  return *this;
139  }
140 
143  return MatrixGradient<N>(dx - G.dx, dy - G.dy);
144  }
145 
146  // Directional derivatives
147  Eigen::Matrix<double, N, N> dx;
148  Eigen::Matrix<double, N, N> dy;
149 
150  };
151 
153  template <size_t N>
154  MatrixGradient<N> operator*(double scalar, MatrixGradient<N> const & G){
155  return MatrixGradient<N>(scalar * G.dx, scalar * G.dy);
156  };
157 
158  //----------------------------------------------------------------------
159  //----------------------------------------------------------------------
160  // SCALAR MONOMIAL BASIS
161  //----------------------------------------------------------------------
162  //----------------------------------------------------------------------
163 
166  {
167  public:
168  typedef double FunctionValue;
171  typedef void DivergenceValue;
172  typedef void RotorValue;
174 
176 
177  constexpr static const TensorRankE tensorRank = Scalar;
178  constexpr static const bool hasAncestor = false;
179  static const bool hasFunction = true;
180  static const bool hasGradient = true;
181  static const bool hasCurl = true;
182  static const bool hasDivergence = false;
183  static const bool hasRotor = false;
184  static const bool hasHessian = true;
185 
188  const Cell &T,
189  size_t degree
190  );
191 
193  inline size_t dimension() const
194  {
195  return ( m_degree >= 0 ? (m_degree + 1) * (m_degree + 2) / 2 : 0);
196  }
197 
199  FunctionValue function(size_t i, const VectorRd & x) const;
200 
202  GradientValue gradient(size_t i, const VectorRd & x) const;
203 
205  CurlValue curl(size_t i, const VectorRd & x) const;
206 
208  HessianValue hessian(size_t i, const VectorRd & x) const;
209 
211  inline size_t max_degree() const
212  {
213  return m_degree;
214  }
215 
217  inline VectorZd powers(size_t i) const
218  {
219  return m_powers[i];
220  }
221 
222  private:
224  inline VectorRd _coordinate_transform(const VectorRd& x) const
225  {
226  return (x - m_xT) / m_hT;
227  }
228 
229  size_t m_degree;
230  VectorRd m_xT;
231  double m_hT;
232  std::vector<VectorZd> m_powers;
233  Eigen::Matrix2d m_rot; // rotation pi/2
234  };
235 
236  //------------------------------------------------------------------------------
237 
240  {
241  public:
242  typedef double FunctionValue;
245  typedef void DivergenceValue;
246  typedef void RotorValue;
247  typedef void HessianValue;
248 
250 
251  constexpr static const TensorRankE tensorRank = Scalar;
252  constexpr static const bool hasAncestor = false;
253  static const bool hasFunction = true;
254  static const bool hasGradient = true;
255  static const bool hasCurl = false;
256  static const bool hasDivergence = false;
257  static const bool hasRotor = false;
258  static const bool hasHessian = false;
259 
262  const Edge &E,
263  size_t degree
264  );
265 
267  inline size_t dimension() const
268  {
269  return m_degree + 1;
270  }
271 
273  FunctionValue function(size_t i, const VectorRd &x) const;
274 
276  GradientValue gradient(size_t i, const VectorRd &x) const;
277 
279  inline size_t max_degree() const
280  {
281  return m_degree;
282  }
283 
284  private:
285  inline double _coordinate_transform(const VectorRd &x) const
286  {
287  return (x - m_xE).dot(m_tE) / m_hE;
288  }
289 
290  size_t m_degree;
291  VectorRd m_xE;
292  double m_hE;
293  VectorRd m_tE;
294  };
295 
296  //------------------------------------------------------------------------------
297  //------------------------------------------------------------------------------
298  // DERIVED BASIS
299  //------------------------------------------------------------------------------
300  //------------------------------------------------------------------------------
301 
302  //----------------------------FAMILY------------------------------------------
303  //
309  template <typename BasisType>
310  class Family
311  {
312  public:
313  typedef typename BasisType::FunctionValue FunctionValue;
314  typedef typename BasisType::GradientValue GradientValue;
315  typedef typename BasisType::CurlValue CurlValue;
316  typedef typename BasisType::DivergenceValue DivergenceValue;
317  typedef typename BasisType::RotorValue RotorValue;
318  typedef typename BasisType::HessianValue HessianValue;
319 
320  typedef typename BasisType::GeometricSupport GeometricSupport;
321 
322  constexpr static const TensorRankE tensorRank = BasisType::tensorRank;
323  constexpr static const bool hasAncestor = true;
324  static const bool hasFunction = BasisType::hasFunction;
325  static const bool hasGradient = BasisType::hasGradient;
326  static const bool hasCurl = BasisType::hasCurl;
327  static const bool hasDivergence = BasisType::hasDivergence;
328  static const bool hasRotor = BasisType::hasRotor;
329  static const bool hasHessian = BasisType::hasHessian;
330 
331  typedef BasisType AncestorType;
332 
335  const BasisType &basis,
336  const Eigen::MatrixXd &matrix
337  )
338  : m_basis(basis),
339  m_matrix(matrix)
340  {
341  assert((size_t)matrix.cols() == basis.dimension() || "Inconsistent family initialization");
342  }
343 
345  inline size_t dimension() const
346  {
347  return m_matrix.rows();
348  }
349 
351  FunctionValue function(size_t i, const VectorRd &x) const
352  {
353  static_assert(hasFunction, "Call to function() not available");
354 
355  FunctionValue f = m_matrix(i, 0) * m_basis.function(0, x);
356  for (auto j = 1; j < m_matrix.cols(); j++)
357  {
358  f += m_matrix(i, j) * m_basis.function(j, x);
359  } // for j
360  return f;
361  }
362 
364  FunctionValue function(size_t i, size_t iqn, const boost::multi_array<FunctionValue, 2> &ancestor_value_quad) const
365  {
366  static_assert(hasFunction, "Call to function() not available");
367 
368  FunctionValue f = m_matrix(i, 0) * ancestor_value_quad[0][iqn];
369  for (auto j = 1; j < m_matrix.cols(); j++)
370  {
371  f += m_matrix(i, j) * ancestor_value_quad[j][iqn];
372  } // for j
373  return f;
374  }
375 
376 
378  GradientValue gradient(size_t i, const VectorRd &x) const
379  {
380  static_assert(hasGradient, "Call to gradient() not available");
381 
382  GradientValue G = m_matrix(i, 0) * m_basis.gradient(0, x);
383  for (auto j = 1; j < m_matrix.cols(); j++)
384  {
385  G += m_matrix(i, j) * m_basis.gradient(j, x);
386  } // for j
387  return G;
388  }
389 
391  GradientValue gradient(size_t i, size_t iqn, const boost::multi_array<GradientValue, 2> &ancestor_gradient_quad) const
392  {
393  static_assert(hasGradient, "Call to gradient() not available");
394 
395  GradientValue G = m_matrix(i, 0) * ancestor_gradient_quad[0][iqn];
396  for (auto j = 1; j < m_matrix.cols(); j++)
397  {
398  G += m_matrix(i, j) * ancestor_gradient_quad[j][iqn];
399  } // for j
400  return G;
401  }
402 
404  CurlValue curl(size_t i, const VectorRd &x) const
405  {
406  static_assert(hasCurl, "Call to curl() not available");
407 
408  CurlValue C = m_matrix(i, 0) * m_basis.curl(0, x);
409  for (auto j = 1; j < m_matrix.cols(); j++)
410  {
411  C += m_matrix(i, j) * m_basis.curl(j, x);
412  } // for j
413  return C;
414  }
415 
417  CurlValue curl(size_t i, size_t iqn, const boost::multi_array<CurlValue, 2> &ancestor_curl_quad) const
418  {
419  static_assert(hasCurl, "Call to curl() not available");
420 
421  CurlValue C = m_matrix(i, 0) * ancestor_curl_quad[0][iqn];
422  for (auto j = 1; j < m_matrix.cols(); j++)
423  {
424  C += m_matrix(i, j) * ancestor_curl_quad[j][iqn];
425  } // for j
426  return C;
427  }
428 
430  DivergenceValue divergence(size_t i, const VectorRd &x) const
431  {
432  static_assert(hasDivergence, "Call to divergence() not available");
433 
434  DivergenceValue D = m_matrix(i, 0) * m_basis.divergence(0, x);
435  for (auto j = 1; j < m_matrix.cols(); j++)
436  {
437  D += m_matrix(i, j) * m_basis.divergence(j, x);
438  } // for j
439  return D;
440  }
441 
443  DivergenceValue divergence(size_t i, size_t iqn, const boost::multi_array<DivergenceValue, 2> &ancestor_divergence_quad) const
444  {
445  static_assert(hasDivergence, "Call to divergence() not available");
446 
447  DivergenceValue D = m_matrix(i, 0) * ancestor_divergence_quad[0][iqn];
448  for (auto j = 1; j < m_matrix.cols(); j++)
449  {
450  D += m_matrix(i, j) * ancestor_divergence_quad[j][iqn];
451  } // for j
452  return D;
453  }
454 
456  DivergenceValue rotor(size_t i, const VectorRd &x) const
457  {
458  static_assert(hasRotor, "Call to rotor() not available");
459 
460  RotorValue R = m_matrix(i, 0) * m_basis.rotor(0, x);
461  for (auto j = 1; j < m_matrix.cols(); j++)
462  {
463  R += m_matrix(i, j) * m_basis.rotor(j, x);
464  } // for j
465  return R;
466  }
467 
469  RotorValue rotor(size_t i, size_t iqn, const boost::multi_array<RotorValue, 2> &ancestor_rotor_quad) const
470  {
471  static_assert(hasRotor, "Call to rotor() not available");
472 
473  RotorValue R = m_matrix(i, 0) * ancestor_rotor_quad[0][iqn];
474  for (auto j = 1; j < m_matrix.cols(); j++)
475  {
476  R += m_matrix(i, j) * ancestor_rotor_quad[j][iqn];
477  } // for j
478  return R;
479  }
480 
482  inline HessianValue hessian(size_t i, const VectorRd & x) const
483  {
484  static_assert(hasHessian, "Call to hessian() not available");
485 
486  HessianValue H = m_matrix(i, 0) * m_basis.hessian(0, x);
487  for (auto j = 1; j < m_matrix.cols(); j++)
488  {
489  H += m_matrix(i, j) * m_basis.hessian(j, x);
490  } // for j
491  return H;
492  }
493 
495  HessianValue hessian(size_t i, size_t iqn, const boost::multi_array<HessianValue, 2> &ancestor_hessian_quad) const
496  {
497  static_assert(hasHessian, "Call to hessian() not available");
498 
499  HessianValue H = m_matrix(i, 0) * ancestor_hessian_quad[0][iqn];
500  for (auto j = 1; j < m_matrix.cols(); j++)
501  {
502  H += m_matrix(i, j) * ancestor_hessian_quad[j][iqn];
503  } // for j
504  return H;
505  }
506 
508  inline const Eigen::MatrixXd & matrix() const
509  {
510  return m_matrix;
511  }
512 
514  constexpr inline const BasisType &ancestor() const
515  {
516  return m_basis;
517  }
518 
520  inline size_t max_degree() const
521  {
522  return m_basis.max_degree();
523  }
524 
525  private:
526  BasisType m_basis;
527  Eigen::MatrixXd m_matrix;
528  };
529 
530  //----------------------TENSORIZED--------------------------------------------------------
531 
533 
557  template <typename ScalarFamilyType, size_t N>
559  {
560  public:
561  typedef typename Eigen::Matrix<double, N, 1> FunctionValue;
562  typedef typename Eigen::Matrix<double, N, dimspace> GradientValue;
564  typedef double DivergenceValue;
565  typedef double RotorValue;
566  typedef void HessianValue;
567 
568  typedef typename ScalarFamilyType::GeometricSupport GeometricSupport;
569 
570  constexpr static const TensorRankE tensorRank = Vector;
571  constexpr static const bool hasAncestor = true;
572  static const bool hasFunction = ScalarFamilyType::hasFunction;
573  static const bool hasGradient = ScalarFamilyType::hasGradient;
574  // We know how to compute the curl, scalar rotor, and divergence if gradient is available
575  static const bool hasDivergence = ( ScalarFamilyType::hasGradient && N==dimspace );
576  static const bool hasRotor = ( ScalarFamilyType::hasGradient && N==dimspace );
577  static const bool hasCurl = ( ScalarFamilyType::hasGradient && N==dimspace );
578  // In future versions, one could include the computation of second derivatives
579  // checking that BasisType has information on the Hessian
580  static const bool hasHessian = false;
581 
582  typedef ScalarFamilyType AncestorType;
583 
584  TensorizedVectorFamily(const ScalarFamilyType & scalar_family)
585  : m_scalar_family(scalar_family)
586  {
587  static_assert(ScalarFamilyType::tensorRank == Scalar,
588  "Vector family can only be constructed from scalar families");
589  m_rot.row(0) << 0., 1.;
590  m_rot.row(1) << -1., 0.;
591 
592  }
593 
595  inline size_t dimension() const
596  {
597  return m_scalar_family.dimension() * N;
598  }
599 
601  FunctionValue function(size_t i, const VectorRd &x) const
602  {
603  static_assert(hasFunction, "Call to function() not available");
604 
605  FunctionValue ek = Eigen::Matrix<double, N, 1>::Zero();
606  ek(i / m_scalar_family.dimension()) = 1.;
607  return ek * m_scalar_family.function(i % m_scalar_family.dimension(), x);
608  }
609 
611  FunctionValue function(size_t i, size_t iqn, const boost::multi_array<double, 2> &ancestor_value_quad) const
612  {
613  static_assert(hasFunction, "Call to function() not available");
614 
615  FunctionValue ek = Eigen::Matrix<double, N, 1>::Zero();
616  ek(i / m_scalar_family.dimension()) = 1.;
617  return ek * ancestor_value_quad[i % m_scalar_family.dimension()][iqn];
618  }
619 
621  GradientValue gradient(size_t i, const VectorRd &x) const
622  {
623  static_assert(hasGradient, "Call to gradient() not available");
624 
625  GradientValue G = Eigen::Matrix<double, N, dimspace>::Zero();
626  G.row(i / m_scalar_family.dimension()) = m_scalar_family.gradient(i % m_scalar_family.dimension(), x);
627  return G;
628  }
629 
631  GradientValue gradient(size_t i, size_t iqn, const boost::multi_array<VectorRd, 2> &ancestor_gradient_quad) const
632  {
633  static_assert(hasGradient, "Call to gradient() not available");
634 
635  GradientValue G = Eigen::Matrix<double, N, dimspace>::Zero();
636  G.row(i / m_scalar_family.dimension()) = ancestor_gradient_quad[i % m_scalar_family.dimension()][iqn];
637  return G;
638  }
639 
641  CurlValue curl(size_t i, const VectorRd &x) const
642  {
643  static_assert(hasCurl, "Call to curl() not available");
644 
645  return m_rot * gradient(i, x);
646  }
647 
649  CurlValue curl(size_t i, size_t iqn, const boost::multi_array<VectorRd, 2> &ancestor_gradient_quad) const
650  {
651  static_assert(hasCurl, "Call to curl() not available");
652 
653  return m_rot * ancestor_gradient_quad[i][iqn];
654  }
655 
657  DivergenceValue divergence(size_t i, const VectorRd &x) const
658  {
659  static_assert(hasDivergence, "Call to divergence() not available");
660 
661  return m_scalar_family.gradient(i % m_scalar_family.dimension(), x)(i / m_scalar_family.dimension());
662  }
663 
665  DivergenceValue divergence(size_t i, size_t iqn, const boost::multi_array<VectorRd, 2> &ancestor_gradient_quad) const
666  {
667  static_assert(hasDivergence, "Call to divergence() not available");
668 
669  return ancestor_gradient_quad[i % m_scalar_family.dimension()][iqn](i / m_scalar_family.dimension());
670  }
671 
673  RotorValue rotor(size_t i, const VectorRd &x) const
674  {
675  static_assert(hasRotor && hasCurl, "Call to rotor() not available");
676 
677  return curl(i, x).trace();
678  }
679 
681  RotorValue rotor(size_t i, size_t iqn, const boost::multi_array<VectorRd, 2> &ancestor_gradient_quad) const
682  {
683  static_assert(hasRotor, "Call to rotor() not available");
684 
685  return curl(i, iqn, ancestor_gradient_quad).trace();
686  }
687 
689  constexpr inline const ScalarFamilyType &ancestor() const
690  {
691  return m_scalar_family;
692  }
693 
695  inline size_t max_degree() const
696  {
697  return m_scalar_family.max_degree();
698  }
699 
700  private:
701  ScalarFamilyType m_scalar_family;
702  Eigen::Matrix2d m_rot; // Rotation of pi/2
703  };
704 
705  //----------------------MATRIX FAMILY--------------------------------------------------------
706 
708 
713  // Useful formulas to manipulate this basis (all indices starting from 0):
714  // the i-th element in the basis is f_{i%r} E_{i/r}
715  // the matrix E_m has 1 at position (m%N, m/N)
716  // the element f_aE_b is the i-th in the family with i = ( [b/N]N + [b%N] )r + a
717  template<typename ScalarFamilyType, size_t N>
719  {
720  public:
721  typedef typename Eigen::Matrix<double, N, N> FunctionValue;
724  typedef typename Eigen::Matrix<double, N, 1> DivergenceValue;
725  typedef void RotorValue;
726  typedef void HessianValue;
727 
728  typedef typename ScalarFamilyType::GeometricSupport GeometricSupport;
729 
730  constexpr static const TensorRankE tensorRank = Matrix;
731  constexpr static const bool hasAncestor = true;
732  static const bool hasFunction = ScalarFamilyType::hasFunction;
733  static const bool hasGradient = true;
734  static const bool hasDivergence = ( ScalarFamilyType::hasGradient && N==dimspace );
735  static const bool hasRotor = false;
736  static const bool hasCurl = false;
737  static const bool hasHessian = false;
738 
739  typedef ScalarFamilyType AncestorType;
740 
741  MatrixFamily(const ScalarFamilyType & scalar_family)
742  : m_scalar_family(scalar_family),
743  m_E(N*N),
744  m_transposeOperator(Eigen::MatrixXd::Zero(scalar_family.dimension()*N*N, scalar_family.dimension()*N*N))
745  {
746  static_assert(ScalarFamilyType::tensorRank == Scalar,
747  "Vector family can only be constructed from scalar families");
748 
749  // Construct the basis for NxN matrices
750  for (size_t j = 0; j < N; j++){
751  for (size_t i = 0; i < N; i++){
752  m_E[j*N + i] = Eigen::Matrix<double, N, N>::Zero();
753  m_E[j*N + i](i,j) = 1.;
754  }
755  }
756 
757  // Construct the transpose operator
758  for (size_t i = 0; i < dimension(); i++){
759  // The i-th basis element is sent to the j-th basis element s.t. j%r=i%r, (j/r)%N=(i/r)/N, (j/r)/N=(i/r)%N
760  size_t r = m_scalar_family.dimension();
761  size_t j = ( ( (i/r)%N )*N + (i/r)/N ) * r + (i%r);
762  m_transposeOperator(i,j) = 1.;
763  }
764  }
765 
767  inline size_t dimension() const
768  {
769  return m_scalar_family.dimension() * N * N;
770  }
771 
773  FunctionValue function(size_t i, const VectorRd & x) const
774  {
775  static_assert(hasFunction, "Call to function() not available");
776 
777  return m_scalar_family.function(i % m_scalar_family.dimension(), x) * m_E[i / m_scalar_family.dimension()];
778  }
779 
781  FunctionValue function(size_t i, size_t iqn, const boost::multi_array<double, 2> &ancestor_value_quad) const
782  {
783  static_assert(hasFunction, "Call to function() not available");
784 
785  return ancestor_value_quad[i % m_scalar_family.dimension()][iqn] * m_E[i / m_scalar_family.dimension()];
786  }
787 
789  DivergenceValue divergence(size_t i, const VectorRd & x) const
790  {
791  static_assert(hasDivergence, "Call to divergence() not available");
792 
793  Eigen::Matrix<double, N, 1> V = m_E[i/m_scalar_family.dimension()].col( (i/m_scalar_family.dimension()) /N);
794  return m_scalar_family.gradient(i % m_scalar_family.dimension(), x)( (i / m_scalar_family.dimension()) / N) * V;
795  }
796 
798  DivergenceValue divergence(size_t i, size_t iqn, const boost::multi_array<VectorRd, 2> &ancestor_gradient_quad) const
799  {
800  static_assert(hasDivergence, "Call to divergence() not available");
801 
802  Eigen::Matrix<double, N, 1> V = m_E[i/m_scalar_family.dimension()].col( (i/m_scalar_family.dimension()) /N);
803  return ancestor_gradient_quad[i % m_scalar_family.dimension()][iqn]( (i / m_scalar_family.dimension()) / N) * V;
804  }
805 
807  GradientValue gradient(size_t i, const VectorRd & x) const
808  {
809  static_assert(hasGradient, "Call to gradient() not available");
810 
811  const VectorRd gradient_scalar_function = m_scalar_family.gradient(i % m_scalar_family.dimension(), x);
812  const Eigen::Matrix<double, N, N> matrix_basis = m_E[i / m_scalar_family.dimension()];
813  return GradientValue(gradient_scalar_function.x() * matrix_basis, gradient_scalar_function.y() * matrix_basis);
814  }
815 
817  GradientValue gradient(size_t i, size_t iqn, const boost::multi_array<VectorRd, 2> &ancestor_gradient_quad) const
818  {
819  static_assert(hasGradient, "Call to gradient() not available");
820 
821  const VectorRd gradient_scalar_function = ancestor_gradient_quad[i % m_scalar_family.dimension()][iqn];
822  const Eigen::Matrix<double, N, N> matrix_basis = m_E[i / m_scalar_family.dimension()];
823  return GradientValue(gradient_scalar_function.x() * matrix_basis, gradient_scalar_function.y() * matrix_basis);
824  }
825 
827  inline const ScalarFamilyType &ancestor() const
828  {
829  return m_scalar_family;
830  }
831 
833  inline const size_t matrixSize() const
834  {
835  return N;
836  }
837 
839  inline const Eigen::MatrixXd transposeOperator() const
840  {
841  return m_transposeOperator;
842  }
843 
845  inline const Eigen::MatrixXd symmetriseOperator() const
846  {
847  return (Eigen::MatrixXd::Identity(dimension(), dimension()) + m_transposeOperator)/2;
848  }
849 
851  const Eigen::MatrixXd symmetricBasis() const
852  {
853  size_t dim_scalar = m_scalar_family.dimension();
854  Eigen::MatrixXd symOpe = symmetriseOperator();
855 
856  Eigen::MatrixXd SB = Eigen::MatrixXd::Zero(dim_scalar * N*(N+1)/2, dimension());
857 
858  // The matrix consists in selecting certain rows of symmetriseOperator, corresponding to the lower triangle part
859  // To select these rows, we loop over the columns of the underlying matrices, and get the rows in symmetriseOperator()
860  // corresponding to the lower triangular part of each column
861  size_t position = 0;
862  for (size_t i = 0; i<N; i++){
863  // Skip the first i*dim_scalar elements in the column (upper triangle) in symOpe, take the remaining (N-i)*dim_scalar
864  SB.middleRows(position, (N-i)*dim_scalar) = symOpe.middleRows(i*(N+1)*dim_scalar, (N-i)*dim_scalar);
865  // update current position
866  position += (N-i)*dim_scalar;
867  }
868 
869  return SB;
870  }
871 
872  private:
873  ScalarFamilyType m_scalar_family;
874  std::vector<Eigen::Matrix<double, N, N>> m_E;
875  Eigen::MatrixXd m_transposeOperator;
876  };
877 
878 
879  //----------------------TANGENT FAMILY--------------------------------------------------------
880 
882  template <typename ScalarFamilyType>
884  {
885  public:
889  typedef void DivergenceValue;
890  typedef void RotorValue;
891  typedef void HessianValue;
892 
894 
895  constexpr static const TensorRankE tensorRank = Vector;
896  constexpr static const bool hasAncestor = true;
897  static const bool hasFunction = true;
898  static const bool hasGradient = false;
899  static const bool hasCurl = false;
900  static const bool hasDivergence = false;
901  static const bool hasRotor = false;
902  static const bool hasHessian = false;
903 
904  typedef ScalarFamilyType AncestorType;
905 
908  const ScalarFamilyType &scalar_family,
909  const VectorRd &generator
910  )
911  : m_scalar_family(scalar_family),
912  m_generator(generator)
913  {
914  static_assert(ScalarFamilyType::hasFunction, "Call to function() not available");
915  static_assert(std::is_same<typename ScalarFamilyType::GeometricSupport, Edge>::value,
916  "Tangent families can only be defined on edges");
917  }
918 
920  inline size_t dimension() const
921  {
922  return m_scalar_family.dimension();
923  }
924 
926  inline FunctionValue function(size_t i, const VectorRd &x) const
927  {
928  return m_scalar_family.function(i, x) * m_generator;
929  }
930 
932  constexpr inline const ScalarFamilyType &ancestor() const
933  {
934  return m_scalar_family;
935  }
936 
938  inline size_t max_degree() const
939  {
940  return m_scalar_family.max_degree();
941  }
942 
944  inline VectorRd generator() const
945  {
946  return m_generator;
947  }
948 
949  private:
950  ScalarFamilyType m_scalar_family;
951  VectorRd m_generator;
952  };
953 
954 
955  //--------------------SHIFTED BASIS----------------------------------------------------------
956 
958 
959  template <typename BasisType>
961  {
962  public:
963  typedef typename BasisType::FunctionValue FunctionValue;
964  typedef typename BasisType::GradientValue GradientValue;
965  typedef typename BasisType::CurlValue CurlValue;
966  typedef typename BasisType::DivergenceValue DivergenceValue;
967  typedef typename BasisType::RotorValue RotorValue;
968  typedef typename BasisType::HessianValue HessianValue;
969 
970  typedef typename BasisType::GeometricSupport GeometricSupport;
971 
972  constexpr static const TensorRankE tensorRank = BasisType::tensorRank;
973  constexpr static const bool hasAncestor = true;
974  static const bool hasFunction = BasisType::hasFunction;
975  static const bool hasGradient = BasisType::hasGradient;
976  static const bool hasCurl = BasisType::hasCurl;
977  static const bool hasDivergence = BasisType::hasDivergence;
978  static const bool hasRotor = BasisType::hasRotor;
979  static const bool hasHessian = BasisType::hasHessian;
980 
981  typedef BasisType AncestorType;
982 
985  const BasisType &basis,
986  const int shift
987  )
988  : m_basis(basis),
989  m_shift(shift)
990  {
991  // Do nothing
992  }
993 
995  inline size_t dimension() const
996  {
997  return m_basis.dimension() - m_shift;
998  }
999 
1001  constexpr inline const BasisType &ancestor() const
1002  {
1003  return m_basis;
1004  }
1005 
1007  inline size_t max_degree() const
1008  {
1009  return m_basis.max_degree();
1010  }
1011 
1013  inline FunctionValue function(size_t i, const VectorRd &x) const
1014  {
1015  static_assert(hasFunction, "Call to function() not available");
1016 
1017  return m_basis.function(i + m_shift, x);
1018  }
1019 
1021  inline GradientValue gradient(size_t i, const VectorRd &x) const
1022  {
1023  static_assert(hasGradient, "Call to gradient() not available");
1024 
1025  return m_basis.gradient(i + m_shift, x);
1026  }
1027 
1029  inline CurlValue curl(size_t i, const VectorRd &x) const
1030  {
1031  static_assert(hasCurl, "Call to curl() not available");
1032 
1033  return m_basis.curl(i + m_shift, x);
1034  }
1035 
1037  inline DivergenceValue divergence(size_t i, const VectorRd &x) const
1038  {
1039  static_assert(hasDivergence, "Call to divergence() not available");
1040 
1041  return m_basis.divergence(i + m_shift, x);
1042  }
1043 
1045  inline RotorValue rotor(size_t i, const VectorRd &x) const
1046  {
1047  static_assert(hasRotor, "Call to rotor() not available");
1048 
1049  return m_basis.rotor(i + m_shift, x);
1050  }
1051 
1053  inline HessianValue hessian(size_t i, const VectorRd & x) const
1054  {
1055  static_assert(hasHessian, "Call to hessian() not available");
1056 
1057  return m_basis.hessian(i + m_shift, x);
1058  }
1059 
1060  private:
1061  BasisType m_basis;
1062  int m_shift;
1063  };
1064 
1065  //-----------------------RESTRICTED BASIS-------------------------------------------------------
1066 
1068 
1069  template <typename BasisType>
1071  {
1072  public:
1073  typedef typename BasisType::FunctionValue FunctionValue;
1074  typedef typename BasisType::GradientValue GradientValue;
1075  typedef typename BasisType::CurlValue CurlValue;
1076  typedef typename BasisType::DivergenceValue DivergenceValue;
1077  typedef typename BasisType::RotorValue RotorValue;
1078  typedef typename BasisType::HessianValue HessianValue;
1079 
1080  typedef typename BasisType::GeometricSupport GeometricSupport;
1081 
1082  constexpr static const TensorRankE tensorRank = BasisType::tensorRank;
1083  constexpr static const bool hasAncestor = true;
1084  static const bool hasFunction = BasisType::hasFunction;
1085  static const bool hasGradient = BasisType::hasGradient;
1086  static const bool hasCurl = BasisType::hasCurl;
1087  static const bool hasDivergence = BasisType::hasDivergence;
1088  static const bool hasRotor = BasisType::hasRotor;
1089  static const bool hasHessian = BasisType::hasHessian;
1090 
1091  typedef BasisType AncestorType;
1092 
1095  const BasisType &basis,
1096  const size_t &dimension
1097  )
1098  : m_basis(basis),
1099  m_dimension(dimension)
1100  {
1101  // Make sure that the provided dimension is smaller than the one
1102  // of the provided basis
1103  assert(dimension <= basis.dimension());
1104  }
1105 
1107  inline size_t dimension() const
1108  {
1109  return m_dimension;
1110  }
1111 
1113  constexpr inline const BasisType &ancestor() const
1114  {
1115  return m_basis;
1116  }
1117 
1119  inline size_t max_degree() const
1120  {
1121  // We need to find the degree based on the dimension, assumes the basis is hierarchical
1122  size_t deg = 0;
1123  while (PolynomialSpaceDimension<GeometricSupport>::Poly(deg) < m_dimension){
1124  deg++;
1125  }
1126  return deg;
1127  }
1128 
1130  inline FunctionValue function(size_t i, const VectorRd &x) const
1131  {
1132  static_assert(hasFunction, "Call to function() not available");
1133 
1134  return m_basis.function(i, x);
1135  }
1136 
1138  GradientValue gradient(size_t i, const VectorRd &x) const
1139  {
1140  static_assert(hasGradient, "Call to gradient() not available");
1141 
1142  return m_basis.gradient(i, x);
1143  }
1144 
1146  CurlValue curl(size_t i, const VectorRd &x) const
1147  {
1148  static_assert(hasCurl, "Call to curl() not available");
1149 
1150  return m_basis.curl(i, x);
1151  }
1152 
1154  DivergenceValue divergence(size_t i, const VectorRd &x) const
1155  {
1156  static_assert(hasDivergence, "Call to divergence() not available");
1157 
1158  return m_basis.divergence(i, x);
1159  }
1160 
1162  RotorValue rotor(size_t i, const VectorRd &x) const
1163  {
1164  static_assert(hasRotor, "Call to rotor() not available");
1165 
1166  return m_basis.rotor(i, x);
1167  }
1168 
1169  private:
1170  BasisType m_basis;
1171  size_t m_dimension;
1172  };
1173 
1174  //--------------------GRADIENT BASIS----------------------------------------------------------
1175 
1177 
1179  template <typename BasisType>
1181  {
1182  public:
1183  typedef typename BasisType::GradientValue FunctionValue;
1184  typedef Eigen::Matrix<double, dimspace, dimspace> GradientValue;
1186  typedef void DivergenceValue;
1187  typedef void RotorValue;
1188  typedef void HessianValue;
1189 
1190  typedef typename BasisType::GeometricSupport GeometricSupport;
1191 
1192  constexpr static const TensorRankE tensorRank = Vector;
1193  constexpr static const bool hasAncestor = true;
1194  static const bool hasFunction = true;
1195  // In future versions, one could include the computation of derivatives
1196  // checking that BasisType has information on the Hessian
1197  static const bool hasGradient = false;
1198  static const bool hasCurl = false;
1199  static const bool hasDivergence = false;
1200  static const bool hasRotor = false;
1201  static const bool hasHessian = false;
1202 
1203  typedef BasisType AncestorType;
1204 
1206  GradientBasis(const BasisType &basis)
1207  : m_scalar_basis(basis)
1208  {
1209  static_assert(BasisType::hasGradient,
1210  "Gradient basis requires gradient() for the original basis to be available");
1211  // Do nothing
1212  }
1213 
1215  inline size_t dimension() const
1216  {
1217  return m_scalar_basis.dimension();
1218  }
1219 
1221  inline FunctionValue function(size_t i, const VectorRd &x) const
1222  {
1223  return m_scalar_basis.gradient(i, x);
1224  }
1225 
1227  constexpr inline const BasisType &ancestor() const
1228  {
1229  return m_scalar_basis;
1230  }
1231 
1232 
1233  private:
1234  BasisType m_scalar_basis;
1235  };
1236 
1237  //-------------------CURL BASIS-----------------------------------------------------------
1238 
1240 
1241  template <typename BasisType>
1243  {
1244  public:
1246  typedef Eigen::Matrix<double, dimspace, dimspace> GradientValue;
1247  typedef Eigen::Matrix<double, dimspace, dimspace> CurlValue;
1248  typedef void DivergenceValue;
1249  typedef void RotorValue;
1250  typedef void HessianValue;
1251 
1252  typedef typename BasisType::GeometricSupport GeometricSupport;
1253 
1254  constexpr static const TensorRankE tensorRank = Vector;
1255  constexpr static const bool hasAncestor = true;
1256  static const bool hasFunction = true;
1257  static const bool hasGradient = false;
1258  static const bool hasCurl = false;
1259  static const bool hasDivergence = false;
1260  static const bool hasRotor = false;
1261  static const bool hasHessian = false;
1262 
1263  typedef BasisType AncestorType;
1264 
1266  CurlBasis(const BasisType &basis)
1267  : m_basis(basis)
1268  {
1269  static_assert(BasisType::tensorRank == Scalar && std::is_same<typename BasisType::GeometricSupport, Cell>::value,
1270  "Curl basis can only be constructed starting from scalar bases on elements");
1271  static_assert(BasisType::hasCurl,
1272  "Curl basis requires curl() for the original basis to be available");
1273  }
1274 
1276  inline size_t dimension() const
1277  {
1278  return m_basis.dimension();
1279  }
1280 
1282  inline FunctionValue function(size_t i, const VectorRd &x) const
1283  {
1284  return m_basis.curl(i, x);
1285  }
1286 
1288  constexpr inline const BasisType &ancestor() const
1289  {
1290  return m_basis;
1291  }
1292 
1293 
1294  private:
1295  size_t m_degree;
1296  BasisType m_basis;
1297  };
1298 
1299 
1300  //--------------------DIVERGENCE BASIS----------------------------------------------------------
1301 
1303 
1304  template <typename BasisType>
1306  {
1307  public:
1308  typedef double FunctionValue;
1311  typedef double DivergenceValue;
1312  typedef void HessianValue;
1313 
1314  typedef typename BasisType::GeometricSupport GeometricSupport;
1315 
1316  constexpr static const TensorRankE tensorRank = Scalar;
1317  constexpr static const bool hasAncestor = true;
1318  static const bool hasFunction = true;
1319  static const bool hasGradient = false;
1320  static const bool hasCurl = false;
1321  static const bool hasDivergence = false;
1322  static const bool hasHessian = false;
1323 
1324  typedef BasisType AncestorType;
1325 
1327  DivergenceBasis(const BasisType &basis)
1328  : m_vector_basis(basis)
1329  {
1330  static_assert(BasisType::tensorRank == Vector || BasisType::tensorRank == Matrix,
1331  "Divergence basis can only be constructed starting from vector bases");
1332  static_assert(BasisType::hasDivergence,
1333  "Divergence basis requires divergence() for the original basis to be available");
1334  // Do nothing
1335  }
1336 
1338  inline size_t dimension() const
1339  {
1340  return m_vector_basis.dimension();
1341  }
1342 
1344  inline FunctionValue function(size_t i, const VectorRd &x) const
1345  {
1346  return m_vector_basis.divergence(i, x);
1347  }
1348 
1350  constexpr inline const BasisType &ancestor() const
1351  {
1352  return m_vector_basis;
1353  }
1354 
1355  private:
1356  BasisType m_vector_basis;
1357  };
1358 
1359 
1360  //--------------------ROTOR BASIS----------------------------------------------------------
1361 
1363 
1364  template <typename BasisType>
1366  {
1367  public:
1368  typedef double FunctionValue;
1371  typedef void DivergenceValue;
1372  typedef void RotorValue;
1373  typedef void HessianValue;
1374 
1375  typedef typename BasisType::GeometricSupport GeometricSupport;
1376 
1377  constexpr static const TensorRankE tensorRank = Scalar;
1378  constexpr static const bool hasAncestor = true;
1379  static const bool hasFunction = true;
1380  static const bool hasGradient = false;
1381  static const bool hasCurl = false;
1382  static const bool hasDivergence = false;
1383  static const bool hasRotor = false;
1384  static const bool hasHessian = false;
1385 
1386  typedef BasisType AncestorType;
1387 
1389  RotorBasis(const BasisType &basis)
1390  : m_vector_basis(basis)
1391  {
1392  static_assert(BasisType::tensorRank == Vector,
1393  "Rotor basis can only be constructed starting from vector bases");
1394  static_assert(BasisType::hasRotor,
1395  "Rotor basis requires rotor() for the original basis to be available");
1396  // Do nothing
1397  }
1398 
1400  inline size_t dimension() const
1401  {
1402  return m_vector_basis.dimension();
1403  }
1404 
1406  inline FunctionValue function(size_t i, const VectorRd &x) const
1407  {
1408  return m_vector_basis.rotor(i, x);
1409  }
1410 
1412  constexpr inline const BasisType &ancestor() const
1413  {
1414  return m_vector_basis;
1415  }
1416 
1417  private:
1418  BasisType m_vector_basis;
1419  };
1420 
1421 
1422  //--------------------HESSIAN BASIS----------------------------------------------------------
1423 
1425 
1427  template <typename BasisType>
1429  {
1430  public:
1432  // The following types do not make sense in this context and are set to void
1433  typedef void GradientValue;
1434  typedef void CurlValue;
1435  typedef void DivergenceValue;
1436  typedef void RotorValue;
1437  typedef void HessianValue;
1438 
1439  typedef typename BasisType::GeometricSupport GeometricSupport;
1440 
1441  constexpr static const TensorRankE tensorRank = Matrix;
1442  constexpr static const bool hasAncestor = true;
1443  static const bool hasFunction = true;
1444  // In future versions, one could include the computation of derivatives
1445  // checking that BasisType has information on the Hessian
1446  static const bool hasGradient = false;
1447  static const bool hasCurl = false;
1448  static const bool hasDivergence = false;
1449  static const bool hasRotor = false;
1450  static const bool hasHessian = false;
1451 
1452  typedef BasisType AncestorType;
1453 
1455  HessianBasis(const BasisType &basis)
1456  : m_scalar_basis(basis)
1457  {
1458  static_assert(BasisType::tensorRank == Scalar,
1459  "Hessian basis can only be constructed starting from scalar bases");
1460  static_assert(BasisType::hasHessian,
1461  "Hessian basis requires hessian() for the original basis to be available");
1462  // Do nothing
1463  }
1464 
1466  inline size_t dimension() const
1467  {
1468  return m_scalar_basis.dimension();
1469  }
1470 
1472  inline FunctionValue function(size_t i, const VectorRd &x) const
1473  {
1474  return m_scalar_basis.hessian(i, x);
1475  }
1476 
1478  constexpr inline const BasisType &ancestor() const
1479  {
1480  return m_scalar_basis;
1481  }
1482 
1483 
1484  private:
1485  BasisType m_scalar_basis;
1486  };
1487 
1488  //---------------------------------------------------------------------
1489  // BASES FOR THE KOSZUL COMPLEMENTS OF G^k, R^k, H^k
1490  //---------------------------------------------------------------------
1493  {
1494  public:
1496  typedef Eigen::Matrix<double, dimspace, dimspace> GradientValue;
1497  typedef Eigen::Matrix<double, dimspace, dimspace> CurlValue;
1498  typedef double DivergenceValue;
1499  typedef void RotorValue;
1500  typedef void HessianValue;
1501 
1503 
1504  constexpr static const TensorRankE tensorRank = Vector;
1505  constexpr static const bool hasAncestor = false;
1506  static const bool hasFunction = true;
1507  static const bool hasGradient = false;
1508  static const bool hasCurl = false;
1509  static const bool hasDivergence = true;
1510  static const bool hasRotor = false;
1511  static const bool hasHessian = false;
1512 
1515  const Cell &T,
1516  size_t degree
1517  );
1518 
1520  inline size_t dimension() const
1521  {
1523  }
1524 
1526  FunctionValue function(size_t i, const VectorRd &x) const;
1527 
1529  DivergenceValue divergence(size_t i, const VectorRd &x) const;
1530 
1532  inline size_t max_degree() const
1533  {
1534  return m_degree;
1535  }
1536 
1538  inline VectorZd powers(size_t i) const
1539  {
1540  return m_powers[i];
1541  }
1542 
1543  private:
1545  inline VectorRd _coordinate_transform(const VectorRd &x) const
1546  {
1547  return (x - m_xT) / m_hT;
1548  }
1549 
1550  size_t m_degree;
1551  VectorRd m_xT;
1552  double m_hT;
1553  std::vector<VectorZd> m_powers;
1554  };
1555 
1556 
1559  {
1560  public:
1562  typedef Eigen::Matrix<double, dimspace, dimspace> GradientValue;
1563  typedef Eigen::Matrix<double, dimspace, dimspace> CurlValue;
1564  typedef void DivergenceValue;
1565  typedef double RotorValue;
1566  typedef void HessianValue;
1567 
1569 
1570  constexpr static const TensorRankE tensorRank = Vector;
1571  constexpr static const bool hasAncestor = false;
1572  static const bool hasFunction = true;
1573  static const bool hasGradient = false;
1574  static const bool hasCurl = false;
1575  static const bool hasDivergence = false;
1576  static const bool hasRotor = true;
1577  static const bool hasHessian = false;
1578 
1581  const Cell &T,
1582  size_t degree
1583  );
1584 
1586  inline size_t dimension() const
1587  {
1589  }
1590 
1592  FunctionValue function(size_t i, const VectorRd &x) const;
1593 
1595  RotorValue rotor(size_t i, const VectorRd &x) const;
1596 
1598  inline const std::shared_ptr<RolyComplBasisCell> &rck() const
1599  {
1600  return m_Rck_basis;
1601  }
1602 
1603  private:
1604  size_t m_degree;
1605  Eigen::Matrix2d m_rot; // Rotation of pi/2: the basis of Gck is the rotated basis of Rck
1606  std::shared_ptr<RolyComplBasisCell> m_Rck_basis; // A basis of Rck, whose rotation gives a basis of Gck
1607  };
1608 
1611  {
1612  public:
1614  // The following types do not make sense in this context and are set to void
1615  typedef void GradientValue;
1616  typedef void CurlValue;
1617  typedef void DivergenceValue;
1618  typedef void RotorValue;
1619  typedef void HessianValue;
1620 
1622 
1623  constexpr static const TensorRankE tensorRank = Matrix;
1624  constexpr static const bool hasAncestor = false;
1625  static const bool hasFunction = true;
1626  static const bool hasGradient = false;
1627  static const bool hasCurl = false;
1628  static const bool hasDivergence = true;
1629  static const bool hasRotor = false;
1630  static const bool hasHessian = false;
1631 
1634  const Cell &T,
1635  size_t degree
1636  );
1637 
1639  inline size_t dimension() const
1640  {
1642  }
1643 
1645  FunctionValue function(size_t i, const VectorRd &x) const;
1646 
1648  inline size_t max_degree() const
1649  {
1650  return m_degree;
1651  }
1652 
1653  private:
1655  inline VectorRd _coordinate_transform(const VectorRd &x) const
1656  {
1657  return (x - m_xT) / m_hT;
1658  }
1659 
1660  size_t m_degree;
1661  VectorRd m_xT;
1662  double m_hT;
1663  Eigen::Matrix2d m_rot;
1664  TensorizedVectorFamily<MonomialScalarBasisCell, dimspace> m_Pkmo;
1665  };
1666 
1667  //------------------------------------------------------------------------------
1668  // Free functions
1669  //------------------------------------------------------------------------------
1670 
1672  {
1680  Hessian
1681  };
1682 
1684  template <typename outValue, typename inValue, typename FunctionType>
1685  boost::multi_array<outValue, 2> transform_values_quad(
1686  const boost::multi_array<inValue, 2> &B_quad,
1687  const FunctionType &F
1688  )
1689  {
1690  boost::multi_array<outValue, 2> transformed_B_quad( boost::extents[B_quad.shape()[0]][B_quad.shape()[1]] );
1691 
1692  std::transform( B_quad.origin(), B_quad.origin() + B_quad.num_elements(), transformed_B_quad.origin(), F);
1693 
1694  return transformed_B_quad;
1695  };
1696 
1698 
1699  template <typename ScalarBasisType, size_t N>
1701  const ScalarBasisType & B,
1702  const std::vector<Eigen::VectorXd> & v
1703  )
1704  {
1705  size_t r = B.dimension();
1706  size_t k = v.size();
1707 
1708  Eigen::MatrixXd M = Eigen::MatrixXd::Zero(r*k, r*N);
1709 
1710  for (size_t i=0; i < k; i++){
1711  // Check the dimension of vector v[i]
1712  assert(v[i].size() == N);
1713 
1714  // Fill in M
1715  for (size_t j=0; j < N; j++){
1716  M.block(i*r, j*r, r, r) = v[i](j) * Eigen::MatrixXd::Identity(r,r);
1717  }
1718  }
1719 
1722  }
1723 
1725  inline static std::function<Eigen::MatrixXd(const Eigen::MatrixXd &)>
1726  symmetrise_matrix = [](const Eigen::MatrixXd & x)->Eigen::MatrixXd { return 0.5*(x+x.transpose());};
1727 
1729  inline static std::function<Eigen::MatrixXd(const Eigen::MatrixXd &)>
1730  skew_symmetrise_matrix = [](const Eigen::MatrixXd & x)->Eigen::MatrixXd { return 0.5*(x-x.transpose());};
1731 
1733 
1734  template <typename ScalarBasisType, size_t N>
1736  const ScalarBasisType & B
1737  )
1738  {
1739  size_t r = B.dimension();
1740  Eigen::MatrixXd M = Eigen::MatrixXd::Zero(r, r*N*N);
1741 
1742  for (size_t k=0; k < N; k++){
1743  M.block(0, k*r*(N+1), r, r) = Eigen::MatrixXd::Identity(r, r);
1744  }
1745 
1747  return Family<MatrixFamily<ScalarBasisType, N>>(matbasis, M);
1748  }
1749  //------------------------------------------------------------------------------
1750  // BASIS EVALUATION
1751  //------------------------------------------------------------------------------
1752 
1753  namespace detail {
1755 
1756  template<typename BasisType, BasisFunctionE BasisFunction>
1758 
1759  //---------- GENERIC EVALUATION TRAITS -----------------------------------------
1760  // Evaluate the function value at x
1761  template<typename BasisType>
1763  {
1764  static_assert(BasisType::hasFunction, "Call to function not available");
1765  typedef typename BasisType::FunctionValue ReturnValue;
1766  static inline ReturnValue evaluate(const BasisType & basis, size_t i, const VectorRd & x)
1767  {
1768  return basis.function(i, x);
1769  }
1770 
1771  // Computes function value at quadrature node iqn, knowing values of ancestor basis at quadrature nodes
1772  static inline ReturnValue evaluate(
1773  const BasisType &basis,
1774  size_t i,
1775  size_t iqn,
1776  const boost::multi_array<ReturnValue, 2> &ancestor_value_quad
1777  )
1778  {
1779  return basis.function(i, iqn, ancestor_value_quad);
1780  }
1781 
1782  };
1783 
1784  // Evaluate the gradient value at x
1785  template<typename BasisType>
1787  {
1788  static_assert(BasisType::hasGradient, "Call to gradient not available");
1789  typedef typename BasisType::GradientValue ReturnValue;
1790  static inline ReturnValue evaluate(const BasisType & basis, size_t i, const VectorRd & x)
1791  {
1792  return basis.gradient(i, x);
1793  }
1794 
1795  // Computes gradient value at quadrature node iqn, knowing gradients of ancestor basis at quadrature nodes
1796  static inline ReturnValue evaluate(
1797  const BasisType &basis,
1798  size_t i,
1799  size_t iqn,
1800  const boost::multi_array<ReturnValue, 2> &ancestor_gradient_quad
1801  )
1802  {
1803  return basis.gradient(i, iqn, ancestor_gradient_quad);
1804  }
1805  };
1806 
1807  // AC - it's supposed to be used only in case of vector variables (vhhospace)
1808  // Evaluate the symmetric part of gradient at x
1809  template<typename BasisType>
1811  {
1812  static_assert(BasisType::hasGradient, "Call to gradient not available");
1813  typedef typename BasisType::GradientValue ReturnValue;
1814  static inline ReturnValue evaluate(const BasisType & basis, size_t i, const VectorRd & x)
1815  {
1816  return 0.5 * (basis.gradient(i, x) + basis.gradient(i, x).transpose());
1817  }
1818 
1819  // Computes gradient value at quadrature node iqn, knowing gradients of ancestor basis at quadrature nodes
1820  static inline ReturnValue evaluate(
1821  const BasisType &basis,
1822  size_t i,
1823  size_t iqn,
1824  const boost::multi_array<ReturnValue, 2> &ancestor_gradient_quad
1825  )
1826  {
1827  return 0.5 * (basis.gradient(i, iqn, ancestor_gradient_quad) + basis.gradient(i, iqn, ancestor_gradient_quad).transpose());
1828  }
1829  };
1830 
1831  // AC - it's supposed to be used only in case of vector variables (vhhospace)
1832  // Evaluate the symmetric part of gradient at x
1833  template<typename BasisType>
1835  {
1836  static_assert(BasisType::hasGradient, "Call to gradient not available");
1837  typedef typename BasisType::GradientValue ReturnValue;
1838  static inline ReturnValue evaluate(const BasisType & basis, size_t i, const VectorRd & x)
1839  {
1840  return 0.5 * (basis.gradient(i, x) - basis.gradient(i, x).transpose());
1841  }
1842 
1843  // Computes gradient value at quadrature node iqn, knowing gradients of ancestor basis at quadrature nodes
1844  static inline ReturnValue evaluate(
1845  const BasisType &basis,
1846  size_t i,
1847  size_t iqn,
1848  const boost::multi_array<ReturnValue, 2> &ancestor_gradient_quad
1849  )
1850  {
1851  return 0.5 * (basis.gradient(i, iqn, ancestor_gradient_quad) - basis.gradient(i, iqn, ancestor_gradient_quad).transpose());
1852  }
1853  };
1854 
1855  // Evaluate the curl value at x
1856  template<typename BasisType>
1857  struct basis_evaluation_traits<BasisType, Curl>
1858  {
1859  static_assert(BasisType::hasCurl, "Call to curl not available");
1860  typedef typename BasisType::CurlValue ReturnValue;
1861  static inline ReturnValue evaluate(const BasisType & basis, size_t i, const VectorRd & x)
1862  {
1863  return basis.curl(i, x);
1864  }
1865 
1866  // Computes curl value at quadrature node iqn, knowing curls of ancestor basis at quadrature nodes
1867  static inline ReturnValue evaluate(
1868  const BasisType &basis,
1869  size_t i,
1870  size_t iqn,
1871  const boost::multi_array<ReturnValue, 2> &ancestor_curl_quad
1872  )
1873  {
1874  return basis.curl(i, iqn, ancestor_curl_quad);
1875  }
1876  };
1877 
1878  // Evaluate the divergence value at x
1879  template<typename BasisType>
1881  {
1882  static_assert(BasisType::hasDivergence, "Call to divergence not available");
1883  typedef typename BasisType::DivergenceValue ReturnValue;
1884  static inline ReturnValue evaluate(const BasisType & basis, size_t i, const VectorRd & x)
1885  {
1886  return basis.divergence(i, x);
1887  }
1888 
1889  // Computes divergence value at quadrature node iqn, knowing divergences of ancestor basis at quadrature nodes
1890  static inline ReturnValue evaluate(
1891  const BasisType &basis,
1892  size_t i,
1893  size_t iqn,
1894  const boost::multi_array<ReturnValue, 2> &ancestor_divergence_quad
1895  )
1896  {
1897  return basis.divergence(i, iqn, ancestor_divergence_quad);
1898  }
1899  };
1900 
1901  // Evaluate the rotor value at x
1902  template<typename BasisType>
1903  struct basis_evaluation_traits<BasisType, Rotor>
1904  {
1905  static_assert(BasisType::hasRotor, "Call to rotor not available");
1906  typedef typename BasisType::RotorValue ReturnValue;
1907  static inline ReturnValue evaluate(const BasisType & basis, size_t i, const VectorRd & x)
1908  {
1909  return basis.rotor(i, x);
1910  }
1911 
1912  // Computes rotor value at quadrature node iqn, knowing rotors of ancestor basis at quadrature nodes
1913  static inline ReturnValue evaluate(
1914  const BasisType &basis,
1915  size_t i,
1916  size_t iqn,
1917  const boost::multi_array<ReturnValue, 2> &ancestor_rotor_quad
1918  )
1919  {
1920  return basis.rotor(i, iqn, ancestor_rotor_quad);
1921  }
1922  };
1923 
1924  // Evaluate the Hessian value at x
1925  template<typename BasisType>
1927  {
1928  static_assert(BasisType::hasHessian, "Call to Hessian not available");
1929  typedef typename BasisType::HessianValue ReturnValue;
1930  static inline ReturnValue evaluate(const BasisType & basis, size_t i, const VectorRd & x)
1931  {
1932  return basis.hessian(i, x);
1933  }
1934 
1935  // Computes Hessian value at quadrature node iqn, knowing Hessians of ancestor basis at quadrature nodes
1936  static inline ReturnValue evaluate(
1937  const BasisType &basis,
1938  size_t i,
1939  size_t iqn,
1940  const boost::multi_array<ReturnValue, 2> &ancestor_hessian_quad
1941  )
1942  {
1943  return basis.hessian(i, iqn, ancestor_hessian_quad);
1944  }
1945  };
1946 
1947  //---------- TENSORIZED FAMILY EVALUATION TRAITS -----------------------------------------
1948  // Evaluate the value at x of a TensorizedVectorFamily (includes information on the ancestor basis, to optimise eval_quad for tensorized bases)
1949  template <typename ScalarBasisType, size_t N>
1951  {
1952  static_assert(TensorizedVectorFamily<ScalarBasisType, N>::hasFunction, "Call to function not available");
1953 
1955  static const BasisFunctionE AncestorBasisFunction = Function; // Type of values needed from ancestor basis
1957 
1958  // Computes function values at x
1959  static inline ReturnValue evaluate(
1961  size_t i,
1962  const VectorRd &x
1963  )
1964  {
1965  return basis.function(i, x);
1966  }
1967 
1968  // Computes function value at quadrature node iqn, knowing ancestor basis at quadrature nodes
1969  static inline ReturnValue evaluate(
1971  size_t i,
1972  size_t iqn,
1973  const boost::multi_array<double, 2> &ancestor_basis_quad
1974  )
1975  {
1976  return basis.function(i, iqn, ancestor_basis_quad);
1977  }
1978  };
1979 
1980  // Evaluate the gradient at x of a TensorizedVectorFamily (includes information on the ancestor basis, to optimise eval_quad for tensorized bases)
1981  template <typename ScalarBasisType, size_t N>
1983  {
1984  static_assert(TensorizedVectorFamily<ScalarBasisType, N>::hasGradient, "Call to gradient not available");
1985 
1987  static const BasisFunctionE AncestorBasisFunction = Gradient; // Type of values needed from ancestor basis
1989 
1990  // Computes gradient value at x
1991  static inline ReturnValue evaluate(
1993  size_t i,
1994  const VectorRd &x
1995  )
1996  {
1997  return basis.gradient(i, x);
1998  }
1999 
2000  // Computes gradient value at quadrature node iqn, knowing ancestor basis at quadrature nodes
2001  static inline ReturnValue evaluate(
2003  size_t i,
2004  size_t iqn,
2005  const boost::multi_array<VectorRd, 2> &ancestor_basis_quad
2006  )
2007  {
2008  return basis.gradient(i, iqn, ancestor_basis_quad);
2009  }
2010  };
2011 
2012  // AC - supposed to be used only in case of vector variables (vhhospace)
2013  // Evaluate the symmetric gradient at x of a TensorizedVectorFamily (includes information on the ancestor basis, to optimise eval_quad for tensorized bases)
2014  template <typename ScalarBasisType, size_t N>
2016  {
2017  static_assert(TensorizedVectorFamily<ScalarBasisType, N>::hasGradient, "Call to gradient not available");
2018 
2020  static const BasisFunctionE AncestorBasisFunction = Gradient; // Type of values needed from ancestor basis
2022 
2023  // Computes gradient value at x
2024  static inline ReturnValue evaluate(
2026  size_t i,
2027  const VectorRd &x
2028  )
2029  {
2030  return 0.5 * (basis.gradient(i, x) + basis.gradient(i, x).transpose()) ;
2031  }
2032 
2033  // Computes gradient value at quadrature node iqn, knowing ancestor basis at quadrature nodes
2034  static inline ReturnValue evaluate(
2036  size_t i,
2037  size_t iqn,
2038  const boost::multi_array<VectorRd, 2> &ancestor_basis_quad
2039  )
2040  {
2041  return 0.5 * (basis.gradient(i, iqn, ancestor_basis_quad) + basis.gradient(i, iqn, ancestor_basis_quad).transpose());
2042  }
2043  };
2044 
2045  // AC - supposed to be used only in case of vector variables (so gradient is a matrix)
2046  // Evaluate the symmetric gradient at x of a TensorizedVectorFamily (includes information on the ancestor basis, to optimise eval_quad for tensorized bases)
2047  template <typename ScalarBasisType, size_t N>
2049  {
2050  static_assert(TensorizedVectorFamily<ScalarBasisType, N>::hasGradient, "Call to gradient not available");
2051 
2053  static const BasisFunctionE AncestorBasisFunction = Gradient; // Type of values needed from ancestor basis
2055 
2056  // Computes gradient value at x
2057  static inline ReturnValue evaluate(
2059  size_t i,
2060  const VectorRd &x
2061  )
2062  {
2063  return 0.5 * (basis.gradient(i, x) - basis.gradient(i, x).transpose()) ;
2064  }
2065 
2066  // Computes gradient value at quadrature node iqn, knowing ancestor basis at quadrature nodes
2067  static inline ReturnValue evaluate(
2069  size_t i,
2070  size_t iqn,
2071  const boost::multi_array<VectorRd, 2> &ancestor_basis_quad
2072  )
2073  {
2074  return 0.5 * (basis.gradient(i, iqn, ancestor_basis_quad) - basis.gradient(i, iqn, ancestor_basis_quad).transpose());
2075  }
2076  };
2077 
2078  // Evaluate the curl at x of a TensorizedVectorFamily (includes information on the ancestor basis, to optimise eval_quad for tensorized bases)
2079  template <typename ScalarBasisType, size_t N>
2081  {
2082  static_assert(TensorizedVectorFamily<ScalarBasisType, N>::hasCurl, "Call to curl not available");
2083 
2085  static const BasisFunctionE AncestorBasisFunction = Gradient; // Type of values needed from ancestor basis
2087 
2088  // Computes curl values at x
2089  static inline ReturnValue evaluate(
2091  size_t i,
2092  const VectorRd &x
2093  )
2094  {
2095  return basis.curl(i, x);
2096  }
2097 
2098  // Computes curl values at quadrature node iqn, knowing ancestor basis at quadrature nodes
2099  static inline ReturnValue evaluate(
2101  size_t i,
2102  size_t iqn,
2103  const boost::multi_array<VectorRd, 2> &ancestor_basis_quad
2104  )
2105  {
2106  return basis.curl(i, iqn, ancestor_basis_quad);
2107  }
2108 
2109  };
2110 
2111  // Evaluate the divergence at x of a TensorizedVectorFamily (includes information on the ancestor basis, to optimise eval_quad for tensorized bases)
2112  template <typename ScalarBasisType, size_t N>
2114  {
2115  static_assert(TensorizedVectorFamily<ScalarBasisType, N>::hasDivergence, "Call to divergence not available");
2116 
2118  static const BasisFunctionE AncestorBasisFunction = Gradient; // Type of values needed from ancestor basis
2120 
2121  // Computes divergence value at x
2122  static inline ReturnValue evaluate(
2124  size_t i,
2125  const VectorRd &x
2126  )
2127  {
2128  return basis.divergence(i, x);
2129  }
2130 
2131  // Computes divergence value at quadrature node iqn, knowing ancestor basis at quadrature nodes
2132  static inline ReturnValue evaluate(
2134  size_t i,
2135  size_t iqn,
2136  const boost::multi_array<VectorRd, 2> &ancestor_basis_quad
2137  )
2138  {
2139  return basis.divergence(i, iqn, ancestor_basis_quad);
2140  }
2141 
2142  };
2143 
2144  // Evaluate the rotor at x of a TensorizedVectorFamily (includes information on the ancestor basis, to optimise eval_quad for tensorized bases)
2145  template <typename ScalarBasisType, size_t N>
2147  {
2148  static_assert(TensorizedVectorFamily<ScalarBasisType, N>::hasRotor, "Call to rotor not available");
2149 
2151  static const BasisFunctionE AncestorBasisFunction = Gradient; // Type of values needed from ancestor basis
2153 
2154  // Computes divergence value at x
2155  static inline ReturnValue evaluate(
2157  size_t i,
2158  const VectorRd &x
2159  )
2160  {
2161  return basis.rotor(i, x);
2162  }
2163 
2164  // Computes rotor value at quadrature node iqn, knowing ancestor basis at quadrature nodes
2165  static inline ReturnValue evaluate(
2167  size_t i,
2168  size_t iqn,
2169  const boost::multi_array<VectorRd, 2> &ancestor_basis_quad
2170  )
2171  {
2172  return basis.rotor(i, iqn, ancestor_basis_quad);
2173  }
2174 
2175  };
2176 
2177  //---------- MATRIX FAMILY EVALUATION TRAITS -----------------------------------------
2178  // Evaluate the value at x of a MatrixFamily (includes information on the ancestor basis, to optimise eval_quad for matrix families)
2179  template <typename ScalarBasisType, size_t N>
2180  struct basis_evaluation_traits<MatrixFamily<ScalarBasisType, N>, Function>
2181  {
2182  static_assert(MatrixFamily<ScalarBasisType, N>::hasFunction, "Call to function not available");
2183 
2185  static const BasisFunctionE AncestorBasisFunction = Function; // Type of values needed from ancestor basis
2187 
2188  // Computes function values at x
2189  static inline ReturnValue evaluate(
2190  const MatrixFamily<ScalarBasisType, N> &basis,
2191  size_t i,
2192  const VectorRd &x
2193  )
2194  {
2195  return basis.function(i, x);
2196  }
2197 
2198  // Computes function value at quadrature node iqn, knowing ancestor basis at quadrature nodes
2199  static inline ReturnValue evaluate(
2200  const MatrixFamily<ScalarBasisType, N> &basis,
2201  size_t i,
2202  size_t iqn,
2203  const boost::multi_array<double, 2> &ancestor_basis_quad
2204  )
2205  {
2206  return basis.function(i, iqn, ancestor_basis_quad);
2207  }
2208  };
2209 
2210  // Evaluate the divergence at x of a MatrixFamily (includes information on the ancestor basis, to optimise eval_quad for matrix families)
2211  template <typename ScalarBasisType, size_t N>
2212  struct basis_evaluation_traits<MatrixFamily<ScalarBasisType, N>, Divergence>
2213  {
2214  static_assert(MatrixFamily<ScalarBasisType, N>::hasDivergence, "Call to divergence not available");
2215 
2217  static const BasisFunctionE AncestorBasisFunction = Gradient; // Type of values needed from ancestor basis
2219 
2220  // Computes divergence value at x
2221  static inline ReturnValue evaluate(
2222  const MatrixFamily<ScalarBasisType, N> &basis,
2223  size_t i,
2224  const VectorRd &x
2225  )
2226  {
2227  return basis.divergence(i, x);
2228  }
2229 
2230  // Computes divergence value at quadrature node iqn, knowing ancestor basis at quadrature nodes
2231  static inline ReturnValue evaluate(
2232  const MatrixFamily<ScalarBasisType, N> &basis,
2233  size_t i,
2234  size_t iqn,
2235  const boost::multi_array<VectorRd, 2> &ancestor_basis_quad
2236  )
2237  {
2238  return basis.divergence(i, iqn, ancestor_basis_quad);
2239  }
2240 
2241  };
2242 
2243 
2244  // Evaluate the gradient at x of a MatrixFamily (includes information on the ancestor basis, to optimise eval_quad for matrix families)
2245  template <typename ScalarBasisType, size_t N>
2246  struct basis_evaluation_traits<MatrixFamily<ScalarBasisType, N>, Gradient>
2247  {
2248  static_assert(MatrixFamily<ScalarBasisType, N>::hasGradient, "Call to gradient not available");
2249 
2251  static const BasisFunctionE AncestorBasisFunction = Gradient; // Type of values needed from ancestor basis
2253 
2254  // Computes gradient value at x
2255  static inline ReturnValue evaluate(
2256  const MatrixFamily<ScalarBasisType, N> &basis,
2257  size_t i,
2258  const VectorRd &x
2259  )
2260  {
2261  return basis.gradient(i, x);
2262  }
2263 
2264  // Computes gradient value at quadrature node iqn, knowing ancestor basis at quadrature nodes
2265  static inline ReturnValue evaluate(
2266  const MatrixFamily<ScalarBasisType, N> &basis,
2267  size_t i,
2268  size_t iqn,
2269  const boost::multi_array<VectorRd, 2> &ancestor_basis_quad
2270  )
2271  {
2272  return basis.gradient(i, iqn, ancestor_basis_quad);
2273  }
2274 
2275  };
2276 
2277  } // end of namespace detail
2278 
2279 
2280  //-----------------------------------EVALUATE_QUAD--------------------------------------
2281 
2283  template<BasisFunctionE BasisFunction>
2284  struct evaluate_quad {
2286  template<typename BasisType>
2287  static boost::multi_array<typename detail::basis_evaluation_traits<BasisType, BasisFunction>::ReturnValue, 2>
2289  const BasisType & basis,
2290  const QuadratureRule & quad
2291  )
2292  {
2294 
2295  boost::multi_array<typename traits::ReturnValue, 2>
2296  basis_quad( boost::extents[basis.dimension()][quad.size()] );
2297 
2298  for (size_t i = 0; i < basis.dimension(); i++) {
2299  for (size_t iqn = 0; iqn < quad.size(); iqn++) {
2300  basis_quad[i][iqn] = traits::evaluate(basis, i, quad[iqn].vector());
2301  } // for iqn
2302  } // for i
2303 
2304  return basis_quad;
2305  }
2306 
2308  template<typename BasisType>
2309  static boost::multi_array<typename detail::basis_evaluation_traits<Family<BasisType>, BasisFunction>::ReturnValue, 2>
2311  const Family<BasisType> & basis,
2312  const QuadratureRule & quad
2313  )
2314  {
2315  typedef detail::basis_evaluation_traits<Family<BasisType>, BasisFunction> traits;
2316 
2317  boost::multi_array<typename traits::ReturnValue, 2>
2318  basis_quad( boost::extents[basis.dimension()][quad.size()] );
2319 
2320  // Compute values at quadrature note on ancestor basis, and then do a transformation
2321  const auto & ancestor_basis = basis.ancestor();
2322  boost::multi_array<typename traits::ReturnValue, 2>
2323  ancestor_basis_quad = evaluate_quad<BasisFunction>::compute(ancestor_basis, quad);
2324 
2325  Eigen::MatrixXd M = basis.matrix();
2326  for (size_t i = 0; i < size_t(M.rows()); i++) {
2327  for (size_t iqn = 0; iqn < quad.size(); iqn++) {
2328  basis_quad[i][iqn] = M(i,0) * ancestor_basis_quad[0][iqn];
2329  }
2330  for (size_t j = 1; j < size_t(M.cols()); j++) {
2331  for (size_t iqn = 0; iqn < quad.size(); iqn++) {
2332  basis_quad[i][iqn] += M(i,j) * ancestor_basis_quad[j][iqn];
2333  }
2334  } // for j
2335  } // for i
2336 
2337  return basis_quad;
2338  }
2339 
2341  template <typename BasisType, size_t N>
2342  static boost::multi_array<typename detail::basis_evaluation_traits<TensorizedVectorFamily<BasisType, N>, BasisFunction>::ReturnValue, 2>
2345  const QuadratureRule &quad
2346  )
2347  {
2349 
2350  boost::multi_array<typename traits::ReturnValue, 2>
2351  basis_quad(boost::extents[basis.dimension()][quad.size()]);
2352 
2353  const auto &ancestor_basis = basis.ancestor();
2354  const boost::multi_array<typename traits::AncestorBasisFunctionValue, 2> ancestor_basis_quad
2356 
2357  for (size_t i = 0; i < basis.dimension(); i++){
2358  for (size_t iqn = 0; iqn < quad.size(); iqn++){
2359  basis_quad[i][iqn] = traits::evaluate(basis, i, iqn, ancestor_basis_quad);
2360  }
2361  }
2362  return basis_quad;
2363  }
2364 
2366  template <typename BasisType, size_t N>
2367  static boost::multi_array<typename detail::basis_evaluation_traits<MatrixFamily<BasisType, N>, BasisFunction>::ReturnValue, 2>
2369  const MatrixFamily<BasisType, N> &basis,
2370  const QuadratureRule &quad
2371  )
2372  {
2373  typedef detail::basis_evaluation_traits<MatrixFamily<BasisType, N>, BasisFunction> traits;
2374 
2375  boost::multi_array<typename traits::ReturnValue, 2>
2376  basis_quad(boost::extents[basis.dimension()][quad.size()]);
2377 
2378  const auto &ancestor_basis = basis.ancestor();
2379  const boost::multi_array<typename traits::AncestorBasisFunctionValue, 2> ancestor_basis_quad
2381 
2382  for (size_t i = 0; i < basis.dimension(); i++){
2383  for (size_t iqn = 0; iqn < quad.size(); iqn++){
2384  basis_quad[i][iqn] = traits::evaluate(basis, i, iqn, ancestor_basis_quad);
2385  }
2386  }
2387  return basis_quad;
2388  }
2389  };
2390 
2391  //------------------------------------------------------------------------------
2392  // ORTHONORMALISATION
2393  //------------------------------------------------------------------------------
2394 
2403  template<typename T>
2404  Eigen::MatrixXd gram_schmidt(
2405  boost::multi_array<T, 2> & basis_eval,
2406  const std::function<double(size_t, size_t)> & inner_product
2407  )
2408  {
2409  auto induced_norm = [inner_product](size_t i)->double
2410  {
2411  return std::sqrt( inner_product(i, i) );
2412  };
2413 
2414  // Number of basis functions
2415  size_t Nb = basis_eval.shape()[0];
2416  // Number of quadrature nodes
2417  size_t Nqn = basis_eval.shape()[1];
2418 
2419  Eigen::MatrixXd B = Eigen::MatrixXd::Zero(Nb, Nb);
2420 
2421  // Normalise the first element
2422  double norm = induced_norm(0);
2423  for (size_t iqn = 0; iqn < Nqn; iqn++) {
2424  basis_eval[0][iqn] /= norm;
2425  } // for iqn
2426  B(0,0) = 1./norm;
2427 
2428  for (size_t ib = 1; ib < Nb; ib++) {
2429  // 'coeffs' represents the coefficients of the ib-th orthogonal function on the ON basis functions 0 to ib-1
2430  Eigen::RowVectorXd coeffs = Eigen::RowVectorXd::Zero(ib);
2431  for (size_t pb = 0; pb < ib; pb++) {
2432  coeffs(pb) = - inner_product(ib, pb);
2433  } // for pb
2434 
2435  // store the values of the orthogonalised ib-th basis function
2436  for (size_t pb = 0; pb < ib; pb++) {
2437  for (size_t iqn = 0; iqn < Nqn; iqn++) {
2438  basis_eval[ib][iqn] += coeffs(pb) * basis_eval[pb][iqn];
2439  } // for iqn
2440  } // for pb
2441 
2442  // normalise ib-th basis function
2443  double norm = induced_norm(ib);
2444  for (size_t iqn = 0; iqn < Nqn; iqn++) {
2445  basis_eval[ib][iqn] /= norm;
2446  } // for iqn
2447  coeffs /= norm;
2448  // Compute ib-th row of B.
2449  // B.topLeftCorner(ib, ib) contains the rows that represent the ON basis functions 0 to ib-1 on
2450  // the original basis. Multiplying on the left by coeffs gives the coefficients of the ON basis
2451  // functions on the original basis functions from 0 to ib-1.
2452  B.block(ib, 0, 1, ib) = coeffs * B.topLeftCorner(ib, ib);
2453  B(ib, ib) = 1./norm;
2454  }
2455  return B;
2456  }
2457 
2458  //------------------------------------------------------------------------------
2459 
2461  double scalar_product(const double & x, const double & y);
2462 
2464  double scalar_product(const VectorRd & x, const VectorRd & y);
2465 
2467  template<int N>
2468  double scalar_product(const Eigen::Matrix<double, N, N> & x, const Eigen::Matrix<double, N, N> & y)
2469  {
2470  return (x.transpose() * y).trace();
2471  }
2472 
2474  template <typename Value>
2475  boost::multi_array<double, 2> scalar_product(
2476  const boost::multi_array<Value, 2> &basis_quad,
2477  const Value &v
2478  )
2479  {
2480  boost::multi_array<double, 2> basis_dot_v_quad(boost::extents[basis_quad.shape()[0]][basis_quad.shape()[1]]);
2481  std::transform(basis_quad.origin(), basis_quad.origin() + basis_quad.num_elements(),
2482  basis_dot_v_quad.origin(), [&v](const Value &x) -> double { return scalar_product(x,v); });
2483  return basis_dot_v_quad;
2484  }
2485 
2487  template <typename Value>
2488  boost::multi_array<double, 2> scalar_product(
2489  const boost::multi_array<Value, 2> & basis_quad,
2490  const std::vector<Value> & v
2491  )
2492  {
2493  boost::multi_array<double, 2> basis_dot_v_quad(boost::extents[basis_quad.shape()[0]][basis_quad.shape()[1]]);
2494  for (std::size_t i = 0; i < basis_quad.shape()[0]; i++) {
2495  for (std::size_t iqn = 0; iqn < basis_quad.shape()[1]; iqn++) {
2496  basis_dot_v_quad[i][iqn] = scalar_product(basis_quad[i][iqn], v[iqn]);
2497  } // for iqn
2498  } // for i
2499  return basis_dot_v_quad;
2500  }
2501 
2503  boost::multi_array<VectorRd, 2> matrix_vector_product(
2504  const boost::multi_array<MatrixRd, 2> & basis_quad,
2505  const std::vector<VectorRd> & v_quad
2506  );
2507 
2509  boost::multi_array<VectorRd, 2> matrix_vector_product(
2510  const std::vector<MatrixRd> & m_quad,
2511  const boost::multi_array<VectorRd, 2> & basis_quad
2512  );
2513 
2515  boost::multi_array<VectorRd, 2> vector_matrix_product(
2516  const std::vector<VectorRd> & v_quad,
2517  const boost::multi_array<MatrixRd, 2> & basis_quad
2518  );
2519 
2521  template<typename BasisType>
2523  const BasisType & basis,
2524  const QuadratureRule & qr,
2525  boost::multi_array<typename BasisType::FunctionValue, 2> & basis_quad
2526  )
2527  {
2528  // Check that the basis evaluation and quadrature rule are coherent
2529  assert ( basis.dimension() == basis_quad.shape()[0] && qr.size() == basis_quad.shape()[1] );
2530 
2531  // The inner product between the i-th and j-th basis vectors
2532  std::function<double(size_t, size_t)> inner_product
2533  = [&basis_quad, &qr](size_t i, size_t j)->double
2534  {
2535  double r = 0.;
2536  for(size_t iqn = 0; iqn < qr.size(); iqn++) {
2537  r += qr[iqn].w * scalar_product(basis_quad[i][iqn], basis_quad[j][iqn]);
2538  } // for iqn
2539  return r;
2540  };
2541 
2542  Eigen::MatrixXd B = gram_schmidt(basis_quad, inner_product);
2543 
2544  return Family<BasisType>(basis, B);
2545  }
2546 
2548  /* This method is more robust that the previous one based on values at quadrature nodes */
2549  template <typename BasisType>
2551  const BasisType &basis,
2552  const Eigen::MatrixXd & GM
2553  )
2554  {
2555  // Check that the basis and Gram matrix are coherent
2556  assert(basis.dimension() == size_t(GM.rows()) && GM.rows() == GM.cols());
2557 
2558  Eigen::MatrixXd L = GM.llt().matrixL();
2559 
2560  return Family<BasisType>(basis, L.inverse());
2561  }
2562 
2563 
2564  //------------------------------------------------------------------------------
2565  // GRAM MATRICES
2566  //------------------------------------------------------------------------------
2567 
2569  template<typename FunctionValue>
2570  Eigen::MatrixXd compute_gram_matrix(const boost::multi_array<FunctionValue, 2> & B1,
2571  const boost::multi_array<FunctionValue, 2> & B2,
2572  const QuadratureRule & qr,
2573  const size_t nrows,
2574  const size_t ncols,
2575  const std::string sym = "nonsym"
2576  )
2577  {
2578  // Check that the basis evaluation and quadrature rule are coherent
2579  assert ( qr.size() == B1.shape()[1] && qr.size() == B2.shape()[1] );
2580  // Check that we don't ask for more members of family than available
2581  assert ( nrows <= B1.shape()[0] && ncols <= B2.shape()[0] );
2582 
2583  // Re-cast quadrature weights into ArrayXd to make computations faster
2584  Eigen::ArrayXd qr_weights = Eigen::ArrayXd::Zero(qr.size());
2585  for (size_t iqn = 0; iqn < qr.size(); iqn++){
2586  qr_weights(iqn) = qr[iqn].w;
2587  }
2588 
2589  Eigen::MatrixXd M = Eigen::MatrixXd::Zero(nrows, ncols);
2590  for (size_t i = 0; i < nrows; i++) {
2591  size_t jcut = 0;
2592  if ( sym=="sym" ) jcut = i;
2593  for (size_t j = 0; j < jcut; j++){
2594  M(i, j) = M(j, i);
2595  }
2596  for(size_t j = jcut; j < ncols; j++) {
2597  std::vector<double> tmp(B1.shape()[1]);
2598  // Extract values at quadrature nodes for elements i of B1 and j of B2
2599  auto B1i = B1[ boost::indices[i][boost::multi_array_types::index_range(0, B1.shape()[1])] ];
2600  auto B2j = B2[ boost::indices[j][boost::multi_array_types::index_range(0, B1.shape()[1])] ];
2601  // Compute scalar product of B1i and B2j, and recast it as ArrayXd
2602  std::transform(B1i.begin(), B1i.end(), B2j.begin(), tmp.begin(), [](FunctionValue a, FunctionValue b)->double { return scalar_product(a,b);});
2603  Eigen::ArrayXd tmp_array = Eigen::Map<Eigen::ArrayXd, Eigen::Unaligned>(tmp.data(), tmp.size());
2604  // Multiply by quadrature weights and sum (using .sum() of ArrayXd makes this step faster than a loop)
2605  M(i,j) = (qr_weights * tmp_array).sum();
2606 
2607 
2608  // Simple version with loop (replaces everything above after the for on j)
2609  /*
2610  for (size_t iqn = 0; iqn < qr.size(); iqn++) {
2611  M(i,j) += qr[iqn].w * scalar_product(B1[i][iqn], B2[j][iqn]);
2612  } // for iqn
2613  */
2614 
2615  } // for j
2616  } // for i
2617  return M;
2618  }
2619 
2621  template<typename FunctionValue>
2622  Eigen::MatrixXd compute_gram_matrix(const boost::multi_array<FunctionValue, 2> & B1,
2623  const boost::multi_array<FunctionValue, 2> & B2,
2624  const QuadratureRule & qr,
2625  const std::string sym = "nonsym"
2626  )
2627  {
2628  return compute_gram_matrix<FunctionValue>(B1, B2, qr, B1.shape()[0], B2.shape()[0], sym);
2629  }
2630 
2632  template<typename FunctionValue>
2633  inline Eigen::MatrixXd compute_gram_matrix(const boost::multi_array<FunctionValue, 2> & B,
2634  const QuadratureRule & qr
2635  )
2636  {
2637  return compute_gram_matrix<FunctionValue>(B, B, qr, "sym");
2638  }
2639 
2641  Eigen::MatrixXd compute_gram_matrix(const boost::multi_array<VectorRd, 2> & B1,
2642  const boost::multi_array<double, 2> & B2,
2643  const QuadratureRule & qr
2644  );
2645 
2647  Eigen::MatrixXd compute_gram_matrix(const boost::multi_array<double, 2> & B1,
2648  const boost::multi_array<double, 2> & B2,
2649  const QuadratureRule & qr,
2650  const size_t nrows,
2651  const size_t ncols,
2652  const std::string sym = "nonsym"
2653  );
2654 
2656  Eigen::MatrixXd compute_gram_matrix(const boost::multi_array<double, 2> & B1,
2657  const boost::multi_array<double, 2> & B2,
2658  const QuadratureRule & qr,
2659  const std::string sym = "nonsym"
2660  );
2661 
2663  Eigen::MatrixXd compute_gram_matrix(const boost::multi_array<VectorRd, 2> & B1,
2664  const boost::multi_array<VectorRd, 2> & B2,
2665  const QuadratureRule & qr,
2666  const size_t nrows,
2667  const size_t ncols,
2668  const std::string sym = "nonsym"
2669  );
2670 
2672  Eigen::MatrixXd compute_gram_matrix(const boost::multi_array<VectorRd, 2> & B1,
2673  const boost::multi_array<VectorRd, 2> & B2,
2674  const QuadratureRule & qr,
2675  const std::string sym = "nonsym"
2676  );
2677 
2679  template<typename ScalarFamilyType, size_t N>
2681  const boost::multi_array<double, 2> & scalar_family_quad,
2682  const QuadratureRule & qr
2683  )
2684  {
2685  Eigen::MatrixXd Gram = Eigen::MatrixXd::Zero(MatFam.dimension(), MatFam.dimension());
2686 
2687  // Gram matrix of the scalar family
2688  Eigen::MatrixXd scalarGram = compute_gram_matrix<double>(scalar_family_quad, qr);
2689  for (size_t i = 0; i < MatFam.dimension(); i = i+scalarGram.rows()){
2690  Gram.block(i, i, scalarGram.rows(), scalarGram.cols()) = scalarGram;
2691  }
2692 
2693  return Gram;
2694  }
2695 
2696 
2698 
2699  template <typename T>
2700  Eigen::VectorXd integrate(
2701  const FType<T> &f,
2702  const BasisQuad<T> &B,
2703  const QuadratureRule &qr,
2704  size_t n_rows = 0
2705  )
2706  {
2707  // If default, set n_rows to size of family
2708  if (n_rows == 0)
2709  {
2710  n_rows = B.shape()[0];
2711  }
2712 
2713  // Number of quadrature nodes
2714  const size_t num_quads = qr.size();
2715 
2716  // Check number of quadrature nodes is compatible with B
2717  assert(num_quads == B.shape()[1]);
2718 
2719  // Check that we don't ask for more members of family than available
2720  assert(n_rows <= B.shape()[0]);
2721 
2722  Eigen::VectorXd V = Eigen::VectorXd::Zero(n_rows);
2723 
2724  for (size_t iqn = 0; iqn < num_quads; iqn++)
2725  {
2726  double qr_weight = qr[iqn].w;
2727  T f_on_qr = f(qr[iqn].vector());
2728  for (size_t i = 0; i < n_rows; i++)
2729  {
2730  V(i) += qr_weight * scalar_product(B[i][iqn], f_on_qr);
2731  }
2732  }
2733 
2734  return V;
2735  }
2736 
2738  template <typename T, typename U>
2740  const FType<U> &f,
2741  const BasisQuad<T> &B1,
2742  const BasisQuad<T> &B2,
2743  const QuadratureRule &qr,
2744  size_t n_rows = 0,
2745  size_t n_cols = 0,
2746  const std::string sym = "nonsym"
2747  )
2748  {
2749  // If default, set n_rows and n_cols to size of families
2750  if (n_rows == 0 && n_cols == 0)
2751  {
2752  n_rows = B1.shape()[0];
2753  n_cols = B2.shape()[0];
2754  }
2755 
2756  // Number of quadrature nodes
2757  const size_t num_quads = qr.size();
2758  // Check number of quadrature nodes is compatible with B1 and B2
2759  assert(num_quads == B1.shape()[1] && num_quads == B2.shape()[1]);
2760  // Check that we don't ask for more members of family than available
2761  assert(n_rows <= B1.shape()[0] && n_cols <= B2.shape()[0]);
2762 
2763  Eigen::MatrixXd M = Eigen::MatrixXd::Zero(n_rows, n_cols);
2764  for (size_t iqn = 0; iqn < num_quads; iqn++)
2765  {
2766  double qr_weight = qr[iqn].w;
2767  U f_on_qr = f(qr[iqn].vector());
2768  for (size_t i = 0; i < n_rows; i++)
2769  {
2770  T f_B1 = f_on_qr * B1[i][iqn];
2771  size_t jcut = 0;
2772  if (sym == "sym")
2773  jcut = i;
2774  for (size_t j = 0; j < jcut; j++)
2775  {
2776  M(i, j) = M(j, i);
2777  }
2778  for (size_t j = jcut; j < n_cols; j++)
2779  {
2780  M(i, j) += qr_weight * scalar_product(f_B1, B2[j][iqn]);
2781  }
2782  }
2783  }
2784  return M;
2785  }
2786 
2788  template <typename T, typename U>
2790  const FType<U> &f,
2791  const BasisQuad<T> &B1,
2792  const BasisQuad<T> &B2,
2793  const QuadratureRule &qr,
2794  const std::string sym
2795  )
2796  {
2797  return compute_weighted_gram_matrix(f, B1, B2, qr, B1.shape()[0], B2.shape()[0], sym);
2798  }
2799 
2801  Eigen::MatrixXd compute_weighted_gram_matrix(
2802  const FType<VectorRd> &f,
2803  const BasisQuad<VectorRd> &B1,
2804  const BasisQuad<double> &B2,
2805  const QuadratureRule &qr,
2806  size_t n_rows = 0,
2807  size_t n_cols = 0
2808  );
2809 
2810 
2812  Eigen::MatrixXd compute_weighted_gram_matrix(
2813  const FType<VectorRd> &f,
2814  const BasisQuad<double> &B1,
2815  const BasisQuad<VectorRd> &B2,
2816  const QuadratureRule &qr,
2817  size_t n_rows = 0,
2818  size_t n_cols = 0
2819  );
2820 
2822  template <typename Value>
2823  Eigen::MatrixXd compute_closure_matrix(
2824  const boost::multi_array<Value, 2> & f_quad,
2825  const boost::multi_array<Value, 2> & g_quad,
2826  const QuadratureRule & qr_f,
2827  const QuadratureRule & qr_g
2828  )
2829  {
2830  const size_t dimf = f_quad.shape()[0];
2831  const size_t dimg = g_quad.shape()[0];
2832  Eigen::MatrixXd M = Eigen::MatrixXd::Zero(dimf, dimg);
2833 
2834  std::vector<Value> int_f(dimf);
2835  for (size_t i = 0; i < dimf; i++){
2836  int_f[i] = qr_f[0].w * f_quad[i][0];
2837  for (size_t iqn = 1; iqn < qr_f.size(); iqn++){
2838  int_f[i] += qr_f[iqn].w * f_quad[i][iqn];
2839  }
2840  }
2841  std::vector<Value> int_g(dimg);
2842  for (size_t j = 0; j < dimg; j++){
2843  int_g[j] = qr_g[0].w * g_quad[j][0];
2844  for (size_t iqn = 1; iqn < qr_g.size(); iqn++){
2845  int_g[j] += qr_g[iqn].w * g_quad[j][iqn];
2846  }
2847  }
2848 
2849  // Compute (int w_i) and create dot product with (int w_j) for j<=i, then fill in M by symmetry
2850  for (size_t i = 0; i < dimf; i++){
2851  for (size_t j = 0; j < dimg; j++){
2852  M(i,j) = scalar_product(int_f[i], int_g[j]);
2853  }
2854  }
2855 
2856  return M;
2857  }
2858 
2860  template <typename Value>
2861  Eigen::MatrixXd compute_closure_matrix(
2862  const boost::multi_array<Value, 2> & f_quad,
2863  const boost::multi_array<Value, 2> & g_quad,
2864  const QuadratureRule & qr
2865  )
2866  {
2867  return compute_closure_matrix(f_quad, g_quad, qr, qr);
2868  }
2869 
2871  template <typename Value>
2872  Eigen::MatrixXd compute_closure_matrix(
2873  const boost::multi_array<Value, 2> & f_quad,
2874  const QuadratureRule & qr
2875  )
2876  {
2877  return compute_closure_matrix(f_quad, f_quad, qr, qr);
2878  }
2879 
2880  //------------------------------------------------------------------------------
2881  // L2 projection of a function
2882  //------------------------------------------------------------------------------
2883 
2885  template<typename BasisType>
2886  Eigen::VectorXd l2_projection(
2887  const std::function<typename BasisType::FunctionValue(const VectorRd &)> & f,
2888  const BasisType & basis,
2889  QuadratureRule & quad,
2890  const boost::multi_array<typename BasisType::FunctionValue, 2> & basis_quad,
2891  const Eigen::MatrixXd &mass_basis = Eigen::MatrixXd::Zero(1,1)
2892  )
2893  {
2894  Eigen::MatrixXd Mass = mass_basis;
2895  if (Mass.norm() < 1e-13){
2896  Mass = compute_gram_matrix(basis_quad, quad);
2897  }
2898 
2899  Eigen::LDLT<Eigen::MatrixXd> cholesky_mass(Mass);
2900  Eigen::VectorXd b = Eigen::VectorXd::Zero(basis.dimension());
2901  for (size_t i = 0; i < basis.dimension(); i++) {
2902  for (size_t iqn = 0; iqn < quad.size(); iqn++) {
2903  b(i) += quad[iqn].w * scalar_product(f(quad[iqn].vector()), basis_quad[i][iqn]);
2904  } // for iqn
2905  } // for i
2906 
2907  return cholesky_mass.solve(b);
2908  }
2909 
2910  //------------------------------------------------------------------------------
2911  // Decomposition on a polynomial basis
2912  //------------------------------------------------------------------------------
2913 
2915 
2920  template <typename BasisType>
2922  {
2924 
2930  const Edge &E,
2931  const BasisType &basis
2932  ):
2933  m_dim(basis.dimension()),
2934  m_basis(basis),
2935  m_on_basis(basis, Eigen::MatrixXd::Identity(m_dim,m_dim)) // Need to initialise before we compute the real ON basis
2936  {
2937  /* The decomposition will be performed using an orthonormalised version of m_basis, for stability.
2938  The scalar product for which we orthonormalise is a minimal one (much less expensive than integrating
2939  over the edge): \sum_i phi(x_i) psi(x_i) where (x_i)_i are nodes that are sufficient to determine
2940  entirely polynomials up to the considered degree.
2941  The nodes we build here correspond to equidistant P^k nodes on the edge */
2942 
2943  // Nodes
2944  m_nb_nodes = basis.max_degree() + 1;
2945  m_nodes.reserve(m_nb_nodes);
2946  double h = 1./std::max(1., double(basis.max_degree()));
2947  for (size_t i=0; i<m_nb_nodes; i++){
2948  VectorRd xi = E.vertex(0)->coords() + double(i) * h * (E.vertex(1)->coords() - E.vertex(0)->coords());
2949  m_nodes.emplace_back(xi.x(), xi.y(), 1.);
2950  }
2951 
2952  // We then orthonormalise Orthonormalise basis for simple scalar product
2953  m_on_basis_nodes.resize(boost::extents[m_dim][m_nb_nodes]);
2956  };
2957 
2960  const Cell &T,
2961  const BasisType &basis
2962  ):
2963  m_dim(basis.dimension()),
2964  m_basis(basis),
2965  m_on_basis(basis, Eigen::MatrixXd::Identity(m_dim,m_dim)) // Need to initialise before we compute the real ON basis
2966  {
2967  /* The decomposition will be performed using an orthonormalised version of m_basis, for stability.
2968  The scalar product for which we orthonormalise is a minimal one (much less expensive than integrating
2969  over the cell): \sum_i phi(x_i) psi(x_i) where (x_i)_i are nodes that are sufficient to determine
2970  entirely polynomials up to the considered degree.
2971  The nodes we build here correspond to P^k nodes on a tetrahedron */
2972 
2973  // Create nodes
2974  std::vector<Eigen::Vector2i> indices = MonomialPowers<Cell>::complete(basis.max_degree());
2975  VectorRd e1 = VectorRd(1., 0.);
2976  VectorRd e2 = VectorRd(0., 1.);
2977  VectorRd x0 = T.center_mass() - T.diam()*VectorRd(1., 1.)/2.0;
2978 
2979  // Nodes
2980  m_nb_nodes = indices.size();
2981  m_nodes.reserve(m_nb_nodes);
2982  for (size_t i=0; i<m_nb_nodes; i++){
2983  VectorRd xi = x0 + T.diam()*(double(indices[i](0)) * e1 + double(indices[i](1)) * e2)/std::max(1.,double(basis.max_degree()));
2984  m_nodes.emplace_back(xi.x(), xi.y(), xi.z(), 1.);
2985  }
2986 
2987  // We then orthonormalise Orthonormalise basis for simple scalar product
2988  m_on_basis_nodes.resize(boost::extents[m_dim][m_nb_nodes]);
2991  };
2992 
2994  inline QuadratureRule get_nodes() const
2995  {
2996  return m_nodes;
2997  };
2998 
3001  boost::multi_array<typename BasisType::FunctionValue, 2> &values
3002  )
3003  {
3004  // Gram matrix of the polynomials and the ON basis
3005  Eigen::MatrixXd Gram_P_onbasis = compute_gram_matrix(values, m_on_basis_nodes, m_nodes);
3006  // L=matrix of P as a family of the ON basis. In theory, Gram_onbasis is identity, but due to rounding errors it
3007  // can be a bit different. Computing L as if it's not the case improves stability
3008  Eigen::MatrixXd Gram_onbasis = compute_gram_matrix(m_on_basis_nodes, m_nodes);
3009  Eigen::MatrixXd L = Gram_onbasis.ldlt().solve(Gram_P_onbasis.transpose());
3010  return Family<BasisType>(m_basis, L.transpose() * m_on_basis.matrix());
3011  };
3012 
3013 
3014  // Member variables
3015  size_t m_dim; // dimension of basis
3016  BasisType m_basis; // Basis on which we decompose
3017  Family<BasisType> m_on_basis; // Orthonormalised basis (for simple scalar product)
3018  boost::multi_array<typename BasisType::FunctionValue, 2> m_on_basis_nodes; // Values of ON basis at the nodes
3019  size_t m_nb_nodes; // nb of nodes
3021 
3022  };
3023 
3024 
3026 
3027 } // end of namespace HArDCore2D
3028 
3029 #endif
Basis for the space of curls (vectorial rot) of polynomials.
Definition: basis.hpp:1243
Basis (or rather family) of divergence of an existing basis.
Definition: basis.hpp:1306
Definition: basis.hpp:311
Basis for the complement G^{c,k}(T) in P^k(F)^2 of the range of the gradient on a face.
Definition: basis.hpp:1559
Basis for the space of gradients of polynomials.
Definition: basis.hpp:1181
Basis for the space of Hessians of polynomials.
Definition: basis.hpp:1429
Basis for the complement H^{c,k}(T) in P^k(T)^{2x2} of the range of the Hessian on a face.
Definition: basis.hpp:1611
Matrix family obtained from a scalar family.
Definition: basis.hpp:719
Scalar monomial basis on a cell.
Definition: basis.hpp:166
Scalar monomial basis on an edge.
Definition: basis.hpp:240
Generate a basis restricted to the first "dimension" functions.
Definition: basis.hpp:1071
Basis for the complement R^{c,k}(T) in P^k(T)^2 of the range of the vectorial rotational on a face.
Definition: basis.hpp:1493
Basis (or rather family) of scalar rotor of an existing basis.
Definition: basis.hpp:1366
Generate a basis where the function indices are shifted.
Definition: basis.hpp:961
Vector family for polynomial functions that are tangent to an edge (determined by the generator)
Definition: basis.hpp:884
Vector family obtained by tensorization of a scalar family.
Definition: basis.hpp:559
Compute max and min eigenvalues of all matrices for i
Definition: compute_eigs.m:5
function[maxeig, mineig, cond, mat]
Definition: compute_eigs.m:1
static ReturnValue evaluate(const TensorizedVectorFamily< ScalarBasisType, N > &basis, size_t i, const VectorRd &x)
Definition: basis.hpp:1991
void RotorValue
Definition: basis.hpp:1436
VectorRd CurlValue
Definition: basis.hpp:1370
size_t dimension() const
Compute the dimension of the basis.
Definition: basis.hpp:1215
constexpr static const bool hasAncestor
Definition: basis.hpp:1571
static const bool hasRotor
Definition: basis.hpp:1510
Eigen::Matrix< double, N, dimspace > GradientValue
Definition: basis.hpp:562
static const bool hasGradient
Definition: basis.hpp:1507
MatrixGradient(Eigen::Matrix< double, N, N > diff_x, Eigen::Matrix< double, N, N > diff_y)
Constructor.
Definition: basis.hpp:113
GradientValue gradient(size_t i, const VectorRd &x) const
Evaluate the gradient of the i-th function at point x.
Definition: basis.hpp:378
BasisType::GradientValue GradientValue
Definition: basis.hpp:314
Cell GeometricSupport
Definition: basis.hpp:1621
Family< MatrixFamily< ScalarBasisType, N > > IsotropicMatrixFamily(const ScalarBasisType &B)
From a scalar family B, constructs a "Family" of "MatrixFamily" (built on B, of size NxN) that repres...
Definition: basis.hpp:1735
static const bool hasFunction
Definition: basis.hpp:253
size_t max_degree() const
Returns the maximum degree of the basis functions.
Definition: basis.hpp:279
static boost::multi_array< typename detail::basis_evaluation_traits< Family< BasisType >, BasisFunction >::ReturnValue, 2 > compute(const Family< BasisType > &basis, const QuadratureRule &quad)
Evaluate a 'Family' of functions at quadrature nodes (optimised compared the generic basis evaluation...
Definition: basis.hpp:2310
MonomialScalarBasisCell(const Cell &T, size_t degree)
Constructor.
Definition: basis.cpp:13
BasisType AncestorType
Definition: basis.hpp:1386
constexpr static const TensorRankE tensorRank
Definition: basis.hpp:322
constexpr static const TensorRankE tensorRank
Definition: basis.hpp:1504
static ReturnValue evaluate(const TensorizedVectorFamily< ScalarBasisType, N > &basis, size_t i, size_t iqn, const boost::multi_array< VectorRd, 2 > &ancestor_basis_quad)
Definition: basis.hpp:2034
static const bool hasDivergence
Definition: basis.hpp:1087
VectorRd FunctionValue
Definition: basis.hpp:1495
static const bool hasFunction
Definition: basis.hpp:974
size_t dimension() const
Return the dimension of the family.
Definition: basis.hpp:920
double DivergenceValue
Definition: basis.hpp:1311
VectorRd CurlValue
Definition: basis.hpp:244
BasisType::RotorValue ReturnValue
Definition: basis.hpp:1905
HessianValue hessian(size_t i, const VectorRd &x) const
Evaluate the Hessian of the i-th basis function at point x.
Definition: basis.hpp:482
static const bool hasGradient
Definition: basis.hpp:573
static const bool hasFunction
Definition: basis.hpp:572
const size_t matrixSize() const
Return the dimension of the matrices in the family.
Definition: basis.hpp:833
DivergenceValue divergence(size_t i, size_t iqn, const boost::multi_array< VectorRd, 2 > &ancestor_gradient_quad) const
Evaluate the divergence of the i-th basis function at a quadrature point iqn, knowing all the gradien...
Definition: basis.hpp:665
BasisType m_basis
Definition: basis.hpp:3016
BasisType::FunctionValue FunctionValue
Definition: basis.hpp:313
size_t dimension() const
Compute the dimension of the basis.
Definition: basis.hpp:1466
static ReturnValue evaluate(const MatrixFamily< ScalarBasisType, N > &basis, size_t i, const VectorRd &x)
Definition: basis.hpp:2221
BasisType::FunctionValue ReturnValue
Definition: basis.hpp:1764
Family< TensorizedVectorFamily< ScalarBasisType, N > > GenericTensorization(const ScalarBasisType &B, const std::vector< Eigen::VectorXd > &v)
From a scalar family B=(B_1..B_r) and vectors (v_1..v_k) in R^N, constructs a "Family" of "Tensorized...
Definition: basis.hpp:1700
MonomialScalarBasisEdge(const Edge &E, size_t degree)
Constructor.
Definition: basis.cpp:62
static const bool hasHessian
Definition: basis.hpp:1511
CurlValue curl(size_t i, const VectorRd &x) const
Evaluate the curl of the i-th basis function at point x.
Definition: basis.hpp:641
static const bool hasDivergence
Definition: basis.hpp:182
size_t dimension() const
Return the dimension of the family.
Definition: basis.hpp:767
void RotorValue
Definition: basis.hpp:1499
BasisType::GeometricSupport GeometricSupport
Definition: basis.hpp:1252
Eigen::MatrixXd compute_closure_matrix(const boost::multi_array< Value, 2 > &f_quad, const boost::multi_array< Value, 2 > &g_quad, const QuadratureRule &qr_f, const QuadratureRule &qr_g)
Computes closure equation matrices, from evaluations at quadrature nodes. For two families of functio...
Definition: basis.hpp:2823
TensorizedVectorFamily< ScalarBasisType, N >::GradientValue ReturnValue
Definition: basis.hpp:2017
CurlBasis(const BasisType &basis)
Constructor.
Definition: basis.hpp:1266
const ScalarFamilyType & ancestor() const
Return the ancestor (family that has been tensorized)
Definition: basis.hpp:827
void HessianValue
Definition: basis.hpp:566
BasisType::HessianValue HessianValue
Definition: basis.hpp:318
MatrixGradient()
Default constructor.
Definition: basis.hpp:123
VectorRd CurlValue
Definition: basis.hpp:1310
BasisFunctionE
Definition: basis.hpp:1672
double FunctionValue
Definition: basis.hpp:1308
void RotorValue
Definition: basis.hpp:172
static const bool hasHessian
Definition: basis.hpp:1630
static const bool hasCurl
Definition: basis.hpp:1447
static const bool hasFunction
Definition: basis.hpp:1318
static ReturnValue evaluate(const BasisType &basis, size_t i, size_t iqn, const boost::multi_array< ReturnValue, 2 > &ancestor_hessian_quad)
Definition: basis.hpp:1936
static const bool hasDivergence
Definition: basis.hpp:900
BasisType::GeometricSupport GeometricSupport
Definition: basis.hpp:320
static ReturnValue evaluate(const BasisType &basis, size_t i, size_t iqn, const boost::multi_array< ReturnValue, 2 > &ancestor_gradient_quad)
Definition: basis.hpp:1844
static const bool hasRotor
Definition: basis.hpp:1449
DivergenceBasis(const BasisType &basis)
Constructor.
Definition: basis.hpp:1327
TensorizedVectorFamily< ScalarBasisType, N >::CurlValue ReturnValue
Definition: basis.hpp:2082
VectorRd CurlValue
Definition: basis.hpp:888
static const bool hasGradient
Definition: basis.hpp:1257
static ReturnValue evaluate(const BasisType &basis, size_t i, size_t iqn, const boost::multi_array< ReturnValue, 2 > &ancestor_divergence_quad)
Definition: basis.hpp:1890
void HessianValue
Definition: basis.hpp:891
Eigen::Vector2d VectorRd
Definition: basis.hpp:55
Eigen::Vector2i VectorZd
Definition: basis.hpp:56
BasisType::RotorValue RotorValue
Definition: basis.hpp:967
void DivergenceValue
Definition: basis.hpp:1564
BasisType::GeometricSupport GeometricSupport
Definition: basis.hpp:1375
Eigen::Matrix< double, dimspace, dimspace > CurlValue
Definition: basis.hpp:1247
static const bool hasHessian
Definition: basis.hpp:902
static const bool hasGradient
Definition: basis.hpp:898
Eigen::Matrix< double, dimspace, dimspace > GradientValue
Definition: basis.hpp:1246
boost::multi_array< VectorRd, 2 > matrix_vector_product(const boost::multi_array< MatrixRd, 2 > &basis_quad, const std::vector< VectorRd > &v_quad)
Take the product of a matrix-valued basis with a vector.
Definition: basis.cpp:175
static boost::multi_array< typename detail::basis_evaluation_traits< TensorizedVectorFamily< BasisType, N >, BasisFunction >::ReturnValue, 2 > compute(const TensorizedVectorFamily< BasisType, N > &basis, const QuadratureRule &quad)
Evaluate a tensorized family at quadrature nodes (optimised compared the generic basis evaluation,...
Definition: basis.hpp:2343
BasisType::DivergenceValue DivergenceValue
Definition: basis.hpp:316
Family(const BasisType &basis, const Eigen::MatrixXd &matrix)
Constructor.
Definition: basis.hpp:334
BasisType::GradientValue GradientValue
Definition: basis.hpp:1074
static ReturnValue evaluate(const MatrixFamily< ScalarBasisType, N > &basis, size_t i, const VectorRd &x)
Definition: basis.hpp:2255
static ReturnValue evaluate(const BasisType &basis, size_t i, const VectorRd &x)
Definition: basis.hpp:1838
static const bool hasRotor
Definition: basis.hpp:1088
static const bool hasRotor
Definition: basis.hpp:735
BasisType AncestorType
Definition: basis.hpp:1452
void GradientValue
Definition: basis.hpp:1615
static const bool hasFunction
Definition: basis.hpp:1625
static ReturnValue evaluate(const BasisType &basis, size_t i, size_t iqn, const boost::multi_array< ReturnValue, 2 > &ancestor_curl_quad)
Definition: basis.hpp:1867
GradientValue gradient(size_t i, const VectorRd &x) const
Evaluate the gradient of the i-th basis function at point x.
Definition: basis.cpp:31
static const bool hasRotor
Definition: basis.hpp:1260
BasisType AncestorType
Definition: basis.hpp:331
static const bool hasDivergence
Definition: basis.hpp:256
static const bool hasCurl
Definition: basis.hpp:255
void HessianValue
Definition: basis.hpp:1188
static const bool hasCurl
Definition: basis.hpp:899
size_t dimension() const
Dimension of the family. This is actually the number of functions in the family, not necessarily line...
Definition: basis.hpp:345
static const bool hasDivergence
Definition: basis.hpp:1448
Family< BasisType > l2_orthonormalize(const BasisType &basis, const QuadratureRule &qr, boost::multi_array< typename BasisType::FunctionValue, 2 > &basis_quad)
-orthonormalization: simply consists in using gram_schmidt() with the specific l2 inner product
Definition: basis.hpp:2522
static ReturnValue evaluate(const MatrixFamily< ScalarBasisType, N > &basis, size_t i, size_t iqn, const boost::multi_array< double, 2 > &ancestor_basis_quad)
Definition: basis.hpp:2199
static const bool hasCurl
Definition: basis.hpp:1574
void GradientValue
Definition: basis.hpp:1433
RotorValue rotor(size_t i, size_t iqn, const boost::multi_array< VectorRd, 2 > &ancestor_gradient_quad) const
Evaluate the scalar rotor of the i-th basis function at a quadrature point iqn, knowing all the gradi...
Definition: basis.hpp:681
TensorizedVectorFamily< ScalarBasisType, N >::DivergenceValue ReturnValue
Definition: basis.hpp:2115
void HessianValue
Definition: basis.hpp:247
RestrictedBasis(const BasisType &basis, const size_t &dimension)
Constructor.
Definition: basis.hpp:1094
constexpr static const bool hasAncestor
Definition: basis.hpp:731
static ReturnValue evaluate(const TensorizedVectorFamily< ScalarBasisType, N > &basis, size_t i, const VectorRd &x)
Definition: basis.hpp:2155
BasisType::CurlValue CurlValue
Definition: basis.hpp:965
static const bool hasFunction
Definition: basis.hpp:1194
MatrixFamily< ScalarBasisType, N >::GradientValue ReturnValue
Definition: basis.hpp:2248
constexpr static const bool hasAncestor
Definition: basis.hpp:252
BasisType::HessianValue HessianValue
Definition: basis.hpp:1078
Edge GeometricSupport
Definition: basis.hpp:893
static ReturnValue evaluate(const TensorizedVectorFamily< ScalarBasisType, N > &basis, size_t i, const VectorRd &x)
Definition: basis.hpp:2024
const Eigen::MatrixXd symmetriseOperator() const
Return the symmetrisation operator, the rN^2 square matrix that to a given vector of coefficients on ...
Definition: basis.hpp:845
GradientBasis(const BasisType &basis)
Constructor.
Definition: basis.hpp:1206
boost::multi_array< T, 2 > BasisQuad
type for a family of basis functions evaluated on quadrature nodes
Definition: basis.hpp:59
DivergenceValue divergence(size_t i, size_t iqn, const boost::multi_array< VectorRd, 2 > &ancestor_gradient_quad) const
Evaluate the divergence of the i-th basis function at a quadrature point iqn, knowing all the gradien...
Definition: basis.hpp:798
VectorRd generator() const
Returns the generator of the basis functions.
Definition: basis.hpp:944
static const bool hasCurl
Definition: basis.hpp:577
boost::multi_array< outValue, 2 > transform_values_quad(const boost::multi_array< inValue, 2 > &B_quad, const FunctionType &F)
Takes an array B_quad of values at quadrature nodes and applies the function F to all of them....
Definition: basis.hpp:1685
GradientValue gradient(size_t i, const VectorRd &x) const
Evaluate the gradient of the i-th basis function at point x.
Definition: basis.hpp:621
GradientValue gradient(size_t i, const VectorRd &x) const
Evaluate the gradient of the i-th basis function at point x.
Definition: basis.cpp:76
static const bool hasDivergence
Definition: basis.hpp:327
size_t dimension() const
Compute the dimension of the basis.
Definition: basis.hpp:1400
void DivergenceValue
Definition: basis.hpp:245
static const bool hasHessian
Definition: basis.hpp:1322
constexpr static const bool hasAncestor
Definition: basis.hpp:1255
VectorRd FunctionValue
Definition: basis.hpp:1561
static const bool hasHessian
Definition: basis.hpp:1201
Cell GeometricSupport
Definition: basis.hpp:1568
const Eigen::MatrixXd transposeOperator() const
Return the transpose operator, the rN^2 square matrix that to a given vector of coefficients on the M...
Definition: basis.hpp:839
GradientValue gradient(size_t i, const VectorRd &x) const
Evaluate the gradient of the i-th basis function at point x.
Definition: basis.hpp:807
void RotorValue
Definition: basis.hpp:890
RotorBasis(const BasisType &basis)
Constructor.
Definition: basis.hpp:1389
VectorRd GradientValue
Definition: basis.hpp:1369
ScalarFamilyType AncestorType
Definition: basis.hpp:739
size_t dimension() const
Return the dimension of the basis.
Definition: basis.hpp:1107
GradientValue gradient(size_t i, size_t iqn, const boost::multi_array< VectorRd, 2 > &ancestor_gradient_quad) const
Evaluate the gradient of the i-th basis function at a quadrature point iqn, knowing all the gradients...
Definition: basis.hpp:631
Cell GeometricSupport
Definition: basis.hpp:175
GradientValue gradient(size_t i, size_t iqn, const boost::multi_array< GradientValue, 2 > &ancestor_gradient_quad) const
Evaluate the gradient of the i-th function at a quadrature point iqn, knowing all the gradients of an...
Definition: basis.hpp:391
FunctionValue function(size_t i, const VectorRd &x) const
Evaluate the i-th basis function at point x.
Definition: basis.hpp:601
MatrixGradient< N > GradientValue
Definition: basis.hpp:722
double scalar_product(const double &x, const double &y)
Scalar product between two reals.
Definition: basis.cpp:163
void HessianValue
Definition: basis.hpp:1437
size_t dimension() const
Dimension of the basis.
Definition: basis.hpp:1586
Eigen::Matrix< double, N, 1 > FunctionValue
Definition: basis.hpp:561
void DivergenceValue
Definition: basis.hpp:1617
Eigen::MatrixXd compute_gram_matrix(const boost::multi_array< VectorRd, 2 > &B1, const boost::multi_array< double, 2 > &B2, const QuadratureRule &qr)
Compute the Gram-like matrix given a family of vector-valued and one of scalar-valued functions by te...
Definition: basis.cpp:225
void DivergenceValue
Definition: basis.hpp:171
BasisType::FunctionValue FunctionValue
Definition: basis.hpp:1073
BasisType::GeometricSupport GeometricSupport
Definition: basis.hpp:1439
constexpr static const TensorRankE tensorRank
Definition: basis.hpp:1623
QuadratureRule m_nodes
Nodes for the interpolation.
Definition: basis.hpp:3020
BasisType AncestorType
Definition: basis.hpp:1091
BasisType::GradientValue ReturnValue
Definition: basis.hpp:1788
constexpr static const bool hasAncestor
Definition: basis.hpp:178
MatrixRd HessianValue
Definition: basis.hpp:173
Eigen::Matrix< double, dimspace, dimspace > GradientValue
Definition: basis.hpp:1184
static ReturnValue evaluate(const TensorizedVectorFamily< ScalarBasisType, N > &basis, size_t i, size_t iqn, const boost::multi_array< VectorRd, 2 > &ancestor_basis_quad)
Definition: basis.hpp:2132
BasisType::RotorValue RotorValue
Definition: basis.hpp:317
static const bool hasGradient
Definition: basis.hpp:1085
BasisType::DivergenceValue ReturnValue
Definition: basis.hpp:1882
static const bool hasCurl
Definition: basis.hpp:736
size_t max_degree() const
Returns the maximum degree of the basis functions.
Definition: basis.hpp:695
static ReturnValue evaluate(const BasisType &basis, size_t i, size_t iqn, const boost::multi_array< ReturnValue, 2 > &ancestor_gradient_quad)
Definition: basis.hpp:1820
static const bool hasHessian
Definition: basis.hpp:329
BasisType::GeometricSupport GeometricSupport
Definition: basis.hpp:970
static const bool hasDivergence
Definition: basis.hpp:575
static const bool hasDivergence
Definition: basis.hpp:1382
static boost::multi_array< typename detail::basis_evaluation_traits< BasisType, BasisFunction >::ReturnValue, 2 > compute(const BasisType &basis, const QuadratureRule &quad)
Generic basis evaluation.
Definition: basis.hpp:2288
static const bool hasHessian
Definition: basis.hpp:737
void HessianValue
Definition: basis.hpp:1250
static const bool hasHessian
Definition: basis.hpp:1261
HessianValue hessian(size_t i, size_t iqn, const boost::multi_array< HessianValue, 2 > &ancestor_hessian_quad) const
Evaluate the hessian of the i-th function at a quadrature point iqn, knowing all the hessian of ances...
Definition: basis.hpp:495
size_t max_degree() const
Returns the maximum degree of the basis functions.
Definition: basis.hpp:1119
TangentFamily(const ScalarFamilyType &scalar_family, const VectorRd &generator)
Constructor.
Definition: basis.hpp:907
RotorValue rotor(size_t i, const VectorRd &x) const
Evaluate the scalar rotor of the i-th basis function at point x.
Definition: basis.hpp:673
size_t dimension() const
Compute the dimension of the basis.
Definition: basis.hpp:1639
DivergenceValue divergence(size_t i, const VectorRd &x) const
Evaluate the divergence of the i-th basis function at point x.
Definition: basis.hpp:657
size_t dimension() const
Compute the dimension of the basis.
Definition: basis.hpp:1520
constexpr static const bool hasAncestor
Definition: basis.hpp:323
constexpr static const bool hasAncestor
Definition: basis.hpp:1317
static ReturnValue evaluate(const TensorizedVectorFamily< ScalarBasisType, N > &basis, size_t i, size_t iqn, const boost::multi_array< VectorRd, 2 > &ancestor_basis_quad)
Definition: basis.hpp:2001
constexpr const ScalarFamilyType & ancestor() const
Return the ancestor (family used for the tangent)
Definition: basis.hpp:932
Eigen::MatrixXd compute_weighted_gram_matrix(const FType< VectorRd > &f, const BasisQuad< VectorRd > &B1, const BasisQuad< double > &B2, const QuadratureRule &qr, size_t n_rows, size_t n_cols)
Computes the Gram-like matrix of integrals (f dot phi_i, phi_j)
Definition: basis.cpp:354
BasisType::FunctionValue FunctionValue
Definition: basis.hpp:963
void HessianValue
Definition: basis.hpp:1500
static const bool hasHessian
Definition: basis.hpp:1577
static ReturnValue evaluate(const MatrixFamily< ScalarBasisType, N > &basis, size_t i, const VectorRd &x)
Definition: basis.hpp:2189
static const bool hasDivergence
Definition: basis.hpp:1321
size_t dimension() const
Return the dimension of the basis.
Definition: basis.hpp:995
VectorRd FunctionValue
Definition: basis.hpp:886
static const bool hasHessian
Definition: basis.hpp:258
MatrixFamily< ScalarBasisType, N >::DivergenceValue ReturnValue
Definition: basis.hpp:2214
static boost::multi_array< typename detail::basis_evaluation_traits< MatrixFamily< BasisType, N >, BasisFunction >::ReturnValue, 2 > compute(const MatrixFamily< BasisType, N > &basis, const QuadratureRule &quad)
Evaluate a Matrix family at quadrature nodes (optimised compared the generic basis evaluation,...
Definition: basis.hpp:2368
constexpr const BasisType & ancestor() const
Return the ancestor (basis that the gradient was taken of)
Definition: basis.hpp:1478
static const bool hasFunction
Definition: basis.hpp:897
RotorValue rotor(size_t i, const VectorRd &x) const
Evaluate the scalar rotor of the i-th basis function at point x.
Definition: basis.hpp:1045
static const bool hasFunction
Definition: basis.hpp:732
size_t dimension() const
Compute the dimension of the basis.
Definition: basis.hpp:193
size_t max_degree() const
Returns the maximum degree of the basis functions.
Definition: basis.hpp:1007
constexpr const ScalarFamilyType & ancestor() const
Return the ancestor (family that has been tensorized)
Definition: basis.hpp:689
constexpr static const TensorRankE tensorRank
Definition: basis.hpp:1254
GradientValue gradient(size_t i, const VectorRd &x) const
Evaluate the gradient of the i-th basis function at point x.
Definition: basis.hpp:1021
ScalarFamilyType::GeometricSupport GeometricSupport
Definition: basis.hpp:728
CurlValue curl(size_t i, const VectorRd &x) const
Evaluate the curl of the i-th basis function at point x.
Definition: basis.hpp:1029
constexpr const BasisType & ancestor() const
Return the ancestor (basis that the gradient was taken of)
Definition: basis.hpp:1227
const std::shared_ptr< RolyComplBasisCell > & rck() const
Return the Rck basis.
Definition: basis.hpp:1598
static ReturnValue evaluate(const MatrixFamily< ScalarBasisType, N > &basis, size_t i, size_t iqn, const boost::multi_array< VectorRd, 2 > &ancestor_basis_quad)
Definition: basis.hpp:2265
FunctionValue function(size_t i, const VectorRd &x) const
Evaluate the i-th basis function at point x.
Definition: basis.hpp:773
boost::multi_array< typename BasisType::FunctionValue, 2 > m_on_basis_nodes
Definition: basis.hpp:3018
double RotorValue
Definition: basis.hpp:1565
void RotorValue
Definition: basis.hpp:1249
Eigen::Matrix2d MatrixRd
Definition: basis.hpp:54
BasisType::DivergenceValue DivergenceValue
Definition: basis.hpp:966
GradientValue CurlValue
Definition: basis.hpp:563
BasisType::CurlValue ReturnValue
Definition: basis.hpp:1859
RotorValue rotor(size_t i, size_t iqn, const boost::multi_array< RotorValue, 2 > &ancestor_rotor_quad) const
Evaluate the scalar rotor of the i-th function at a quadrature point iqn, knowing all the scalar roto...
Definition: basis.hpp:469
static const bool hasDivergence
Definition: basis.hpp:734
static ReturnValue evaluate(const TensorizedVectorFamily< ScalarBasisType, N > &basis, size_t i, size_t iqn, const boost::multi_array< VectorRd, 2 > &ancestor_basis_quad)
Definition: basis.hpp:2099
CurlValue curl(size_t i, const VectorRd &x) const
Evaluate the curl of the i-th basis function at point x.
Definition: basis.hpp:1146
static const bool hasRotor
Definition: basis.hpp:1200
static ReturnValue evaluate(const TensorizedVectorFamily< ScalarBasisType, N > &basis, size_t i, const VectorRd &x)
Definition: basis.hpp:2057
constexpr static const TensorRankE tensorRank
Definition: basis.hpp:1441
static ReturnValue evaluate(const TensorizedVectorFamily< ScalarBasisType, N > &basis, size_t i, const VectorRd &x)
Definition: basis.hpp:1959
void RotorValue
Definition: basis.hpp:1618
static std::function< Eigen::MatrixXd(const Eigen::MatrixXd &)> symmetrise_matrix
Function to symmetrise a matrix (useful together with transform_values_quad)
Definition: basis.hpp:1726
constexpr const BasisType & ancestor() const
Return the ancestor (basis that the gradient was taken of)
Definition: basis.hpp:1412
DivergenceValue divergence(size_t i, const VectorRd &x) const
Evaluate the divergence of the i-th basis function at point x.
Definition: basis.hpp:1037
VectorZd powers(size_t i) const
Returns the powers of the i-th basis function (its degree can be found using powers(i)....
Definition: basis.hpp:217
DivergenceValue divergence(size_t i, const VectorRd &x) const
Evaluate the divergence of the i-th function at point x.
Definition: basis.hpp:430
constexpr static const TensorRankE tensorRank
Definition: basis.hpp:1192
BasisType::GradientValue FunctionValue
Definition: basis.hpp:1183
static const bool hasGradient
Definition: basis.hpp:1573
Eigen::Matrix< double, dimspace, dimspace > CurlValue
Definition: basis.hpp:1563
size_t max_degree() const
Returns the maximum degree of the basis functions.
Definition: basis.hpp:211
VectorRd FunctionValue
Definition: basis.hpp:1245
void DivergenceValue
Definition: basis.hpp:1186
void DivergenceValue
Definition: basis.hpp:889
BasisType::GradientValue GradientValue
Definition: basis.hpp:964
static ReturnValue evaluate(const TensorizedVectorFamily< ScalarBasisType, N > &basis, size_t i, size_t iqn, const boost::multi_array< VectorRd, 2 > &ancestor_basis_quad)
Definition: basis.hpp:2067
GradientValue gradient(size_t i, const VectorRd &x) const
Evaluate the gradient of the i-th basis function at point x.
Definition: basis.hpp:1138
static const bool hasHessian
Definition: basis.hpp:1384
double DivergenceValue
Definition: basis.hpp:1498
DivergenceValue rotor(size_t i, const VectorRd &x) const
Evaluate the scalar rotor of the i-th function at point x.
Definition: basis.hpp:456
static ReturnValue evaluate(const TensorizedVectorFamily< ScalarBasisType, N > &basis, size_t i, size_t iqn, const boost::multi_array< VectorRd, 2 > &ancestor_basis_quad)
Definition: basis.hpp:2165
VectorRd GradientValue
Definition: basis.hpp:887
constexpr static const TensorRankE tensorRank
Definition: basis.hpp:972
constexpr static const bool hasAncestor
Definition: basis.hpp:1193
Eigen::VectorXd l2_projection(const std::function< typename BasisType::FunctionValue(const VectorRd &)> &f, const BasisType &basis, QuadratureRule &quad, const boost::multi_array< typename BasisType::FunctionValue, 2 > &basis_quad, const Eigen::MatrixXd &mass_basis=Eigen::MatrixXd::Zero(1, 1))
Compute the L2-projection of a function.
Definition: basis.hpp:2886
TensorizedVectorFamily< ScalarBasisType, N >::FunctionValue ReturnValue
Definition: basis.hpp:1952
void CurlValue
Definition: basis.hpp:1616
BasisType AncestorType
Definition: basis.hpp:1263
static const bool hasCurl
Definition: basis.hpp:1086
static const bool hasCurl
Definition: basis.hpp:1198
static const bool hasFunction
Definition: basis.hpp:1506
CurlValue curl(size_t i, const VectorRd &x) const
Evaluate the curl (vector rot) of the i-th basis function at point x.
Definition: basis.cpp:41
static const bool hasRotor
Definition: basis.hpp:1383
void DivergenceValue
Definition: basis.hpp:1435
static const bool hasGradient
Definition: basis.hpp:254
Cell GeometricSupport
Definition: basis.hpp:1502
MatrixFamily(const ScalarFamilyType &scalar_family)
Definition: basis.hpp:741
constexpr int dimspace
Dimension, and generic types for vector in correct dimension (makes it easier to translate a code bet...
Definition: basis.hpp:53
TensorRankE
Definition: basis.hpp:65
RotorValue rotor(size_t i, const VectorRd &x) const
Evaluate the rotor of the i-th basis function at point x.
Definition: basis.cpp:131
double FunctionValue
Definition: basis.hpp:168
constexpr static const bool hasAncestor
Definition: basis.hpp:973
GolyComplBasisCell(const Cell &T, size_t degree)
Constructor.
Definition: basis.cpp:117
Eigen::MatrixXd gram_schmidt(boost::multi_array< T, 2 > &basis_eval, const std::function< double(size_t, size_t)> &inner_product)
Definition: basis.hpp:2404
constexpr static const TensorRankE tensorRank
Definition: basis.hpp:1082
CurlValue curl(size_t i, const VectorRd &x) const
Evaluate the curl of the i-th function at point x.
Definition: basis.hpp:404
static ReturnValue evaluate(const TensorizedVectorFamily< ScalarBasisType, N > &basis, size_t i, const VectorRd &x)
Definition: basis.hpp:2089
VectorRd CurlValue
Definition: basis.hpp:723
static const bool hasRotor
Definition: basis.hpp:183
BasisType::GradientValue ReturnValue
Definition: basis.hpp:1812
DivergenceValue divergence(size_t i, const VectorRd &x) const
Evaluate the divergence of the i-th basis function at point x.
Definition: basis.hpp:1154
void RotorValue
Definition: basis.hpp:1187
MatrixGradient operator+(const MatrixGradient &G)
Addition.
Definition: basis.hpp:130
constexpr const BasisType & ancestor() const
Return the ancestor (basis that the gradient was taken of)
Definition: basis.hpp:1288
static const bool hasCurl
Definition: basis.hpp:326
BasisType::CurlValue CurlValue
Definition: basis.hpp:315
static ReturnValue evaluate(const BasisType &basis, size_t i, const VectorRd &x)
Definition: basis.hpp:1861
static const bool hasGradient
Definition: basis.hpp:1197
MatrixRd FunctionValue
Definition: basis.hpp:1613
static const bool hasGradient
Definition: basis.hpp:1319
Family< BasisType > m_on_basis
Definition: basis.hpp:3017
static std::function< Eigen::MatrixXd(const Eigen::MatrixXd &)> skew_symmetrise_matrix
Function to skew-symmetrise a matrix (useful together with transform_values_quad)
Definition: basis.hpp:1730
static const bool hasCurl
Definition: basis.hpp:1508
void RotorValue
Definition: basis.hpp:725
constexpr const BasisType & ancestor() const
Return the underlying complete basis.
Definition: basis.hpp:1001
TensorizedVectorFamily< ScalarBasisType, N >::GradientValue ReturnValue
Definition: basis.hpp:2050
BasisType::HessianValue HessianValue
Definition: basis.hpp:968
constexpr static const TensorRankE tensorRank
Definition: basis.hpp:1316
static ReturnValue evaluate(const TensorizedVectorFamily< ScalarBasisType, N > &basis, size_t i, const VectorRd &x)
Definition: basis.hpp:2122
const Eigen::MatrixXd & matrix() const
Return the coefficient matrix.
Definition: basis.hpp:508
DivergenceValue divergence(size_t i, size_t iqn, const boost::multi_array< DivergenceValue, 2 > &ancestor_divergence_quad) const
Evaluate the divergence of the i-th function at a quadrature point iqn, knowing all the divergences o...
Definition: basis.hpp:443
constexpr static const TensorRankE tensorRank
Definition: basis.hpp:251
double DivergenceValue
Definition: basis.hpp:564
std::function< T(const VectorRd &)> FType
type for function of point. T is the type of value of the function
Definition: basis.hpp:62
BasisType AncestorType
Definition: basis.hpp:1203
static const bool hasHessian
Definition: basis.hpp:184
DivergenceValue divergence(size_t i, const VectorRd &x) const
Evaluate the divergence of the i-th basis function at point x.
Definition: basis.hpp:789
RotorValue rotor(size_t i, const VectorRd &x) const
Evaluate the rotor of the i-th basis function at point x.
Definition: basis.hpp:1162
CurlValue curl(size_t i, size_t iqn, const boost::multi_array< CurlValue, 2 > &ancestor_curl_quad) const
Evaluate the curl of the i-th function at a quadrature point iqn, knowing all the curls of ancestor b...
Definition: basis.hpp:417
static const bool hasRotor
Definition: basis.hpp:576
Eigen::Matrix< double, dimspace, dimspace > GradientValue
Definition: basis.hpp:1496
Eigen::Matrix< double, N, N > dx
Definition: basis.hpp:147
static const bool hasGradient
Definition: basis.hpp:1380
static const bool hasDivergence
Definition: basis.hpp:1628
MatrixFamily< ScalarBasisType, N >::FunctionValue ReturnValue
Definition: basis.hpp:2182
BasisType::GeometricSupport GeometricSupport
Definition: basis.hpp:1314
static const bool hasDivergence
Definition: basis.hpp:1259
static const bool hasCurl
Definition: basis.hpp:976
Eigen::Matrix< double, N, N > FunctionValue
Definition: basis.hpp:721
VectorRd GradientValue
Definition: basis.hpp:1309
static ReturnValue evaluate(const BasisType &basis, size_t i, const VectorRd &x)
Definition: basis.hpp:1907
static const bool hasRotor
Definition: basis.hpp:328
size_t dimension() const
Compute the dimension of the basis.
Definition: basis.hpp:1276
static const bool hasFunction
Definition: basis.hpp:1256
ScalarFamilyType AncestorType
Definition: basis.hpp:904
size_t dimension() const
Return the dimension of the family.
Definition: basis.hpp:595
static ReturnValue evaluate(const BasisType &basis, size_t i, const VectorRd &x)
Definition: basis.hpp:1930
static const bool hasCurl
Definition: basis.hpp:1258
size_t dimension() const
Compute the dimension of the basis.
Definition: basis.hpp:1338
HolyComplBasisCell(const Cell &T, size_t degree)
Constructor.
Definition: basis.cpp:141
static const bool hasGradient
Definition: basis.hpp:1446
static const bool hasFunction
Definition: basis.hpp:179
DivergenceValue divergence(size_t i, const VectorRd &x) const
Evaluate the divergence of the i-th basis function at point x.
Definition: basis.cpp:106
static ReturnValue evaluate(const BasisType &basis, size_t i, size_t iqn, const boost::multi_array< ReturnValue, 2 > &ancestor_rotor_quad)
Definition: basis.hpp:1913
BasisType::GeometricSupport GeometricSupport
Definition: basis.hpp:1080
VectorRd CurlValue
Definition: basis.hpp:1185
BasisType AncestorType
Definition: basis.hpp:981
void HessianValue
Definition: basis.hpp:1312
GradientValue gradient(size_t i, size_t iqn, const boost::multi_array< VectorRd, 2 > &ancestor_gradient_quad) const
Evaluate the gradient of the i-th basis function at a quadrature point iqn, knowing all the gradients...
Definition: basis.hpp:817
size_t max_degree() const
Returns the maximum degree of the basis functions.
Definition: basis.hpp:1532
MatrixRd FunctionValue
Definition: basis.hpp:1431
Eigen::Matrix< double, dimspace, dimspace > GradientValue
Definition: basis.hpp:1562
VectorRd GradientValue
Definition: basis.hpp:243
size_t max_degree() const
Returns the maximum degree of the basis functions.
Definition: basis.hpp:1648
DecomposePoly(const Cell &T, const BasisType &basis)
Constructor for cell.
Definition: basis.hpp:2959
static const bool hasGradient
Definition: basis.hpp:733
static const bool hasRotor
Definition: basis.hpp:257
static const bool hasCurl
Definition: basis.hpp:1381
static ReturnValue evaluate(const BasisType &basis, size_t i, const VectorRd &x)
Definition: basis.hpp:1884
constexpr static const TensorRankE tensorRank
Definition: basis.hpp:730
static const bool hasRotor
Definition: basis.hpp:1629
void RotorValue
Definition: basis.hpp:1372
static const bool hasDivergence
Definition: basis.hpp:1199
static const bool hasHessian
Definition: basis.hpp:979
void HessianValue
Definition: basis.hpp:1566
constexpr static const TensorRankE tensorRank
Definition: basis.hpp:570
static const bool hasFunction
Definition: basis.hpp:1084
static ReturnValue evaluate(const BasisType &basis, size_t i, size_t iqn, const boost::multi_array< ReturnValue, 2 > &ancestor_gradient_quad)
Definition: basis.hpp:1796
const Eigen::MatrixXd symmetricBasis() const
Returns the matrix that can be used, in a Family of this MatrixBasis, to create a basis of the subspa...
Definition: basis.hpp:851
static const bool hasRotor
Definition: basis.hpp:1576
size_t max_degree() const
Returns the maximum degree of the basis functions.
Definition: basis.hpp:938
VectorRd GradientValue
Definition: basis.hpp:169
HessianValue hessian(size_t i, const VectorRd &x) const
Evaluate the Hessian of the i-th basis function at point x.
Definition: basis.cpp:46
double RotorValue
Definition: basis.hpp:565
BasisType::DivergenceValue DivergenceValue
Definition: basis.hpp:1076
BasisType::GradientValue ReturnValue
Definition: basis.hpp:1836
constexpr const BasisType & ancestor() const
Return the underlying complete basis.
Definition: basis.hpp:1113
static ReturnValue evaluate(const BasisType &basis, size_t i, const VectorRd &x)
Definition: basis.hpp:1766
ScalarFamilyType::GeometricSupport GeometricSupport
Definition: basis.hpp:568
Eigen::VectorXd integrate(const FType< T > &f, const BasisQuad< T > &B, const QuadratureRule &qr, size_t n_rows=0)
Compute the integral of a given function against all functions from a family.
Definition: basis.hpp:2700
static std::vector< VectorZd > complete(size_t degree)
Definition: basis.hpp:86
static ReturnValue evaluate(const BasisType &basis, size_t i, size_t iqn, const boost::multi_array< ReturnValue, 2 > &ancestor_value_quad)
Definition: basis.hpp:1772
VectorZd powers(size_t i) const
Returns the powers of the i-th basis function (not including the vector part)
Definition: basis.hpp:1538
static const bool hasHessian
Definition: basis.hpp:580
size_t max_degree() const
Returns the maximum degree of the basis functions.
Definition: basis.hpp:520
HessianBasis(const BasisType &basis)
Constructor.
Definition: basis.hpp:1455
static const bool hasRotor
Definition: basis.hpp:901
Family< BasisType > family(boost::multi_array< typename BasisType::FunctionValue, 2 > &values)
Returns the decomposed polynomials as a Family of the provided basis.
Definition: basis.hpp:3000
void HessianValue
Definition: basis.hpp:726
Eigen::Matrix< double, N, N > dy
Definition: basis.hpp:148
BasisType::GeometricSupport GeometricSupport
Definition: basis.hpp:1190
constexpr static const bool hasAncestor
Definition: basis.hpp:1378
static ReturnValue evaluate(const BasisType &basis, size_t i, const VectorRd &x)
Definition: basis.hpp:1814
static const bool hasFunction
Definition: basis.hpp:324
static const bool hasGradient
Definition: basis.hpp:180
ShiftedBasis(const BasisType &basis, const int shift)
Constructor.
Definition: basis.hpp:984
constexpr const BasisType & ancestor() const
Return the ancestor (basis that the gradient was taken of)
Definition: basis.hpp:1350
HessianValue hessian(size_t i, const VectorRd &x) const
Evaluate the Hessian of the i-th basis function at point x.
Definition: basis.hpp:1053
ScalarFamilyType AncestorType
Definition: basis.hpp:582
constexpr static const TensorRankE tensorRank
Definition: basis.hpp:895
constexpr static const bool hasAncestor
Definition: basis.hpp:1505
MatrixGradient< N > operator*(double scalar, MatrixGradient< N > const &G)
Multiplication of a MatrixGradient from the left (non-member function to be able to multiply from the...
Definition: basis.hpp:154
void HessianValue
Definition: basis.hpp:1619
constexpr static const bool hasAncestor
Definition: basis.hpp:1083
constexpr static const bool hasAncestor
Definition: basis.hpp:1442
void CurlValue
Definition: basis.hpp:1434
constexpr static const TensorRankE tensorRank
Definition: basis.hpp:1377
static const bool hasHessian
Definition: basis.hpp:1450
CurlValue curl(size_t i, size_t iqn, const boost::multi_array< VectorRd, 2 > &ancestor_gradient_quad) const
Evaluate the curl of the i-th basis function at a quadrature point iqn, knowing all the gradients of ...
Definition: basis.hpp:649
static const bool hasGradient
Definition: basis.hpp:325
constexpr const BasisType & ancestor() const
Return the ancestor.
Definition: basis.hpp:514
size_t dimension() const
Dimension of the basis.
Definition: basis.hpp:267
constexpr static const bool hasAncestor
Definition: basis.hpp:896
void DivergenceValue
Definition: basis.hpp:1248
static const bool hasFunction
Definition: basis.hpp:1443
BasisType AncestorType
Definition: basis.hpp:1324
static const bool hasDivergence
Definition: basis.hpp:1509
constexpr static const TensorRankE tensorRank
Definition: basis.hpp:1570
TensorizedVectorFamily< ScalarBasisType, N >::RotorValue ReturnValue
Definition: basis.hpp:2148
TensorizedVectorFamily< ScalarBasisType, N >::GradientValue ReturnValue
Definition: basis.hpp:1984
MatrixGradient & operator+=(const MatrixGradient &G)
Increment.
Definition: basis.hpp:135
QuadratureRule get_nodes() const
Return the set of nodes (useful to compute value of polynomial to decompose via evaluate_quad)
Definition: basis.hpp:2994
static const bool hasFunction
Definition: basis.hpp:1572
constexpr static const bool hasAncestor
Definition: basis.hpp:571
DecomposePoly(const Edge &E, const BasisType &basis)
Constructor for edge.
Definition: basis.hpp:2929
void HessianValue
Definition: basis.hpp:1373
RolyComplBasisCell(const Cell &T, size_t degree)
Constructor.
Definition: basis.cpp:85
size_t m_dim
Definition: basis.hpp:3011
static const bool hasCurl
Definition: basis.hpp:1627
static const bool hasCurl
Definition: basis.hpp:181
Edge GeometricSupport
Definition: basis.hpp:249
static ReturnValue evaluate(const TensorizedVectorFamily< ScalarBasisType, N > &basis, size_t i, size_t iqn, const boost::multi_array< double, 2 > &ancestor_basis_quad)
Definition: basis.hpp:1969
boost::multi_array< VectorRd, 2 > vector_matrix_product(const std::vector< VectorRd > &v_quad, const boost::multi_array< MatrixRd, 2 > &basis_quad)
Take the product of (the transposed of) a vector with a matrix-valued basis.
Definition: basis.cpp:205
double FunctionValue
Definition: basis.hpp:1368
static ReturnValue evaluate(const MatrixFamily< ScalarBasisType, N > &basis, size_t i, size_t iqn, const boost::multi_array< VectorRd, 2 > &ancestor_basis_quad)
Definition: basis.hpp:2231
BasisType::HessianValue ReturnValue
Definition: basis.hpp:1928
Eigen::Matrix< double, N, 1 > DivergenceValue
Definition: basis.hpp:724
static const bool hasDivergence
Definition: basis.hpp:977
constexpr static const bool hasAncestor
Definition: basis.hpp:1624
size_t m_nb_nodes
Definition: basis.hpp:3019
static const bool hasDivergence
Definition: basis.hpp:1575
static const bool hasRotor
Definition: basis.hpp:978
static const bool hasGradient
Definition: basis.hpp:975
MatrixGradient operator-(const MatrixGradient &G)
Subtraction.
Definition: basis.hpp:142
double FunctionValue
Definition: basis.hpp:242
static const bool hasCurl
Definition: basis.hpp:1320
constexpr static const TensorRankE tensorRank
Definition: basis.hpp:177
static const bool hasHessian
Definition: basis.hpp:1089
static const bool hasGradient
Definition: basis.hpp:1626
static ReturnValue evaluate(const BasisType &basis, size_t i, const VectorRd &x)
Definition: basis.hpp:1790
TensorizedVectorFamily(const ScalarFamilyType &scalar_family)
Definition: basis.hpp:584
static const bool hasFunction
Definition: basis.hpp:1379
BasisType::CurlValue CurlValue
Definition: basis.hpp:1075
Eigen::Matrix< double, dimspace, dimspace > CurlValue
Definition: basis.hpp:1497
void RotorValue
Definition: basis.hpp:246
void DivergenceValue
Definition: basis.hpp:1371
BasisType::RotorValue RotorValue
Definition: basis.hpp:1077
VectorRd CurlValue
Definition: basis.hpp:170
@ Curl
Definition: basis.hpp:1677
@ Function
Definition: basis.hpp:1673
@ Gradient
Definition: basis.hpp:1674
@ SkewsymmetricGradient
Definition: basis.hpp:1676
@ Rotor
Definition: basis.hpp:1679
@ Divergence
Definition: basis.hpp:1678
@ SymmetricGradient
Definition: basis.hpp:1675
@ Hessian
Definition: basis.hpp:1680
@ Scalar
Definition: basis.hpp:66
@ Vector
Definition: basis.hpp:67
@ Matrix
Definition: basis.hpp:68
size_t L
Definition: HHO_DiffAdvecReac.hpp:46
Polytope< DIMENSION > Cell
Definition: Polytope2D.hpp:151
Polytope< 1 > Edge
A Face is a Polytope with object_dim = DIMENSION - 1.
Definition: Polytope2D.hpp:147
std::vector< QuadratureNode > QuadratureRule
Definition: quadraturerule.hpp:55
if(strcmp(field, 'real')) % real valued entries T
Definition: mmread.m:93
for j
Definition: mmread.m:174
Definition: PastixInterface.hpp:16
Definition: ddr-klplate.hpp:27
static auto v
Definition: ddrcore-test.hpp:32
Structure to decompose a set of polynomials on a basis on an edge or a cell.
Definition: basis.hpp:2922
Structure to store the gradient of a matrix-valued function.
Definition: basis.hpp:111
Compute vectors listing the powers of monomial basis functions (for a cell, only this specialization ...
Definition: basis.hpp:79
Basis dimensions for various polynomial spaces on edges/faces/elements (when relevant): Pk,...
Definition: polynomialspacedimension.hpp:55
Basis evaluation traits. Only specialization of 'BasisFunction' (=Function, Gradient,...
Definition: basis.hpp:1757
Evaluate a basis at quadrature nodes. 'BasisFunction' (=Function, Gradient, Curl, Divergence,...
Definition: basis.hpp:2284