Rheolef  7.1
an efficient C++ finite element environment
form_element.h
Go to the documentation of this file.
1 #ifndef _RHEO_FORM_ELEMENT_H
2 #define _RHEO_FORM_ELEMENT_H
23 //
24 // for checking the spectral convergence exp(-ak)
25 // on the reference element
26 //
27 // author: Pierre.Saramito@imag.fr
28 //
29 // date: 1 october 2017
30 //
31 #include "rheolef/quadrature.h"
32 #include "rheolef/basis.h"
33 
34 namespace rheolef {
35 
36 // -----------------------------------------------------------------------------
37 // mass
38 // -----------------------------------------------------------------------------
39 template<class Value>
41 product(const Value& a, const Value& b) { return a*b; }
42 template<class T>
43 T
44 product(const point_basic<T>& a, const point_basic<T>& b) { return dot(a,b); }
45 
46 template<class T, class Value, class Basis>
47 void
49  const Basis& b,
50  reference_element hat_K,
51  Eigen::SparseMatrix<T,Eigen::RowMajor>& mass)
52 {
54  const T eps = 1e3*std::numeric_limits<T>::epsilon();
55  quadrature_option qopt;
56  qopt.set_order(2*b.degree());
57  quadrature<Float> quad (qopt);
58  size_type loc_ndof = b.ndof (hat_K);
59  Eigen::Matrix<Value,Eigen::Dynamic,1> phi (loc_ndof);
60  mass.resize (loc_ndof, loc_ndof);
61  for (typename quadrature<T>::const_iterator iter_q = quad.begin(hat_K),
62  last_q = quad.end(hat_K); iter_q != last_q; iter_q++) {
63  b.evaluate (hat_K, (*iter_q).x, phi);
64  for (size_type loc_idof = 0; loc_idof < loc_ndof; ++loc_idof) {
65  for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
66  T coeff_w = (*iter_q).w*product(phi[loc_idof],phi[loc_jdof]);
67  if (fabs(coeff_w) < eps) continue;
68  mass.coeffRef(loc_idof,loc_jdof) += coeff_w;
69  }}
70  }
71  // TODO: some zeros could appear as sum upon quadrature (e.g. Pk_sherwin):
72  // => do a second pass to eliminate it
73  mass.makeCompressed();
74 }
75 template<class T, class Basis>
76 void
78  const Basis& b,
79  reference_element hat_K,
80  Eigen::SparseMatrix<T,Eigen::RowMajor>& mass)
81 {
82  switch (b.valued_tag()) {
83  case space_constant::scalar: build_mass_generic<T,T,Basis> (b, hat_K, mass); break;
84  case space_constant::vector: build_mass_generic<T,point_basic<T>,Basis> (b, hat_K, mass); break;
85  default: error_macro("tensor-like mass_element test: not yet");
86  }
87 }
88 // -----------------------------------------------------------------------------
89 // mass summed on element boundary, for "Pkd[sides]" tests
90 // -----------------------------------------------------------------------------
91 template<class T, class Basis>
92 void
94  const Basis& b,
95  reference_element tilde_K,
96  Eigen::SparseMatrix<T,Eigen::RowMajor>& mass_bdr)
97 {
99  const T eps = 1e3*std::numeric_limits<T>::epsilon();
100  quadrature_option qopt;
101  qopt.set_order(2*b.degree());
102  quadrature<Float> quad (qopt);
103  size_type loc_ndof = b.ndof (tilde_K);
104  Eigen::Matrix<T,Eigen::Dynamic,1> phi (loc_ndof);
105  mass_bdr.resize (loc_ndof, loc_ndof);
106  for (size_type loc_isid = 0, loc_nsid = tilde_K.n_side(); loc_isid < loc_nsid; ++loc_isid) {
107  side_information_type sid (tilde_K, loc_isid);
108  reference_element hat_S = sid.hat;
109  T tilde_Js = tilde_K.side_measure(loc_isid)/measure(hat_S);
110  for (typename quadrature<T>::const_iterator iter_q = quad.begin(hat_S),
111  last_q = quad.end(hat_S); iter_q != last_q; iter_q++) {
112  T tilde_wq = tilde_Js*(*iter_q).w;
113  if (b.option().is_restricted_to_sides()) {
114  b.evaluate_on_side (tilde_K, sid, (*iter_q).x, phi);
115  } else {
116  b.evaluate (tilde_K, (*iter_q).x, phi);
117  }
118  for (size_type loc_idof = 0; loc_idof < loc_ndof; ++loc_idof) {
119  for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
120  T coeff_w = tilde_wq*product(phi[loc_idof],phi[loc_jdof]);
121  if (fabs(coeff_w) < eps) continue;
122  mass_bdr.coeffRef(loc_idof,loc_jdof) += coeff_w;
123  }}
124  }
125  }
126 }
127 // -----------------------------------------------------------------------------
128 // grad_grad
129 // -----------------------------------------------------------------------------
130 template<class T, class Basis>
131 void
133  const Basis& b,
134  reference_element hat_K,
135  Eigen::SparseMatrix<T,Eigen::RowMajor>& grad_grad)
136 {
138  const T eps = 1e3*std::numeric_limits<T>::epsilon();
139  quadrature_option qopt;
140  qopt.set_order(2*b.degree());
141  quadrature<Float> quad (qopt);
142  size_type loc_ndof = b.ndof (hat_K);
143  Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> grad_phi (loc_ndof);
144  grad_grad.resize (loc_ndof, loc_ndof);
145  for (typename quadrature<T>::const_iterator iter_q = quad.begin(hat_K),
146  last_q = quad.end(hat_K); iter_q != last_q; iter_q++) {
147  b.grad_evaluate (hat_K, (*iter_q).x, grad_phi);
148  for (size_type loc_idof = 0; loc_idof < loc_ndof; ++loc_idof) {
149  for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
150  T coeff_w = (*iter_q).w*dot(grad_phi[loc_idof],grad_phi[loc_jdof]);
151  if (fabs(coeff_w) < eps) continue;
152  grad_grad.coeffRef (loc_idof,loc_jdof) += coeff_w;
153  }}
154  }
155 }
156 
157 } // namespace rheolef
158 #endif // _RHEO_FORM_ELEMENT_H
field::size_type size_type
Definition: branch.cc:425
Float phi(const point &nu, Float a, Float b)
rep::const_iterator const_iterator
Definition: quadrature.h:195
const_iterator end(reference_element hat_K) const
Definition: quadrature.h:219
const_iterator begin(reference_element hat_K) const
Definition: quadrature.h:218
see the reference_element page for the full documentation
Float side_measure(size_type loc_isid) const
std::vector< int >::size_type size_type
integrate_option quadrature_option
size_t size_type
Definition: basis_get.cc:76
Expr1::float_type T
Definition: field_expr.h:261
This file is part of Rheolef.
void build_mass_bdr(const Basis &b, reference_element tilde_K, Eigen::SparseMatrix< T, Eigen::RowMajor > &mass_bdr)
Definition: form_element.h:93
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
dot(x,y): see the expression page for the full documentation
Definition: vec_expr_v2.h:415
Float measure(reference_element hat_K)
void build_mass(const Basis &b, reference_element hat_K, Eigen::SparseMatrix< T, Eigen::RowMajor > &mass)
Definition: form_element.h:77
void build_mass_generic(const Basis &b, reference_element hat_K, Eigen::SparseMatrix< T, Eigen::RowMajor > &mass)
Definition: form_element.h:48
void build_grad_grad(const Basis &b, reference_element hat_K, Eigen::SparseMatrix< T, Eigen::RowMajor > &grad_grad)
Definition: form_element.h:132
scalar_traits< Value >::type product(const Value &a, const Value &b)
Definition: form_element.h:41
Definition: phi.h:25
Float epsilon