dune-pdelab  2.4-dev
laplace.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 
4 #ifndef DUNE_PDELAB_LOCALOPERATOR_LAPLACE_HH
5 #define DUNE_PDELAB_LOCALOPERATOR_LAPLACE_HH
6 
7 #include <vector>
8 
9 #include <dune/common/fvector.hh>
10 
11 #include <dune/geometry/type.hh>
12 #include <dune/geometry/quadraturerules.hh>
13 
14 #include <dune/localfunctions/common/interfaceswitch.hh>
15 
18 
19 namespace Dune {
20  namespace PDELab {
24 
36  class Laplace
37  : public FullVolumePattern,
39  {
40  public:
41  // pattern assembly flags
42  enum { doPatternVolume = true };
43 
44  // residual assembly flags
45  enum { doAlphaVolume = true };
46 
51  Laplace (unsigned int quadOrder)
52  : quadOrder_(quadOrder)
53  {}
54 
65  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
66  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
67  {
68  // domain and range field type
69  typedef FiniteElementInterfaceSwitch<
70  typename LFSU::Traits::FiniteElementType
71  > FESwitch;
72  typedef BasisInterfaceSwitch<
73  typename FESwitch::Basis
74  > BasisSwitch;
75  typedef typename BasisSwitch::DomainField DF;
76  typedef typename BasisSwitch::RangeField RF;
77 
78  // dimensions
79  static const int dimLocal = EG::Geometry::mydimension;
80  static const int dimGlobal = EG::Geometry::coorddimension;
81 
82  // select quadrature rule
83  Dune::GeometryType gt = eg.geometry().type();
84  const Dune::QuadratureRule<DF,dimLocal>& rule =
85  Dune::QuadratureRules<DF,dimLocal>::rule(gt,quadOrder_);
86 
87  // loop over quadrature points
88  for(typename Dune::QuadratureRule<DF,dimLocal>::const_iterator it =
89  rule.begin(); it!=rule.end(); ++it)
90  {
91  // evaluate gradient of shape functions
92  // (we assume Galerkin method lfsu=lfsv)
93  std::vector<Dune::FieldMatrix<RF,1,dimGlobal> >
94  gradphiu(lfsu.size());
95  BasisSwitch::gradient(FESwitch::basis(lfsu.finiteElement()),
96  eg.geometry(), it->position(), gradphiu);
97  std::vector<Dune::FieldMatrix<RF,1,dimGlobal> >
98  gradphiv(lfsv.size());
99  BasisSwitch::gradient(FESwitch::basis(lfsv.finiteElement()),
100  eg.geometry(), it->position(), gradphiv);
101 
102  // compute gradient of u
103  Dune::FieldVector<RF,dimGlobal> gradu(0.0);
104  for (size_t i=0; i<lfsu.size(); i++)
105  gradu.axpy(x(lfsu,i),gradphiu[i][0]);
106 
107  // integrate grad u * grad phi_i
108  RF factor = r.weight() * it->weight() * eg.geometry().integrationElement(it->position());
109  for (size_t i=0; i<lfsv.size(); i++)
110  r.rawAccumulate(lfsv,i,(gradu*gradphiv[i][0])*factor);
111  }
112  }
113 
124  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
125  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, M & matrix) const
126  {
127  // Switches between local and global interface
128  typedef FiniteElementInterfaceSwitch<
129  typename LFSU::Traits::FiniteElementType
130  > FESwitch;
131  typedef BasisInterfaceSwitch<
132  typename FESwitch::Basis
133  > BasisSwitch;
134 
135  // domain and range field type
136  typedef typename BasisSwitch::DomainField DF;
137  typedef typename BasisSwitch::RangeField RF;
138  typedef typename LFSU::Traits::SizeType size_type;
139 
140  // dimensions
141  const int dim = EG::Entity::dimension;
142 
143  // select quadrature rule
144  Dune::GeometryType gt = eg.geometry().type();
145  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,quadOrder_);
146 
147  // loop over quadrature points
148  for (typename QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
149  {
150  std::vector<Dune::FieldMatrix<RF,1,dim> > gradphi(lfsu.size());
151  BasisSwitch::gradient(FESwitch::basis(lfsu.finiteElement()),
152  eg.geometry(), it->position(), gradphi);
153 
154  // geometric weight
155  RF factor = it->weight() * eg.geometry().integrationElement(it->position());
156 
157  for (size_type i=0; i<lfsu.size(); i++)
158  {
159  for (size_type j=0; j<lfsv.size(); j++)
160  {
161  // integrate grad u * grad phi
162  matrix.accumulate(lfsv,j,lfsu,i, gradphi[i][0] * gradphi[j][0] * factor);
163  }
164  }
165  }
166  }
167 
168  protected:
169  // Quadrature rule order
170  unsigned int quadOrder_;
171  };
172 
174  } // namespace PDELab
175 } // namespace Dune
176 
177 #endif
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &matrix) const
Compute the Laplace stiffness matrix for the element given in 'eg'.
Definition: laplace.hh:125
Laplace(unsigned int quadOrder)
Constructor.
Definition: laplace.hh:51
Default flags for all local operators.
Definition: flags.hh:18
sparsity pattern generator
Definition: pattern.hh:13
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Compute Laplace matrix times a given vector for one element.
Definition: laplace.hh:66
Definition: adaptivity.hh:27
Definition: laplace.hh:36
static const int dim
Definition: adaptivity.hh:83
unsigned int quadOrder_
Definition: laplace.hh:170
Definition: laplace.hh:45
const EG & eg
Definition: constraints.hh:280