dune-pdelab  2.4.1
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 #warning This file is deprecated and will be removed after the Dune-PDELab 2.4 release! Use the ConvectionDiffusionFEM local operator from dune/pdelab/localoperator/convectiondiffusionfem.hh instead!
7 
8 #include <vector>
9 
10 #include <dune/common/fvector.hh>
11 
12 #include <dune/geometry/type.hh>
13 #include <dune/geometry/quadraturerules.hh>
14 
15 #include <dune/localfunctions/common/interfaceswitch.hh>
16 
19 
20 namespace Dune {
21  namespace PDELab {
25 
37  class Laplace
38  : public FullVolumePattern,
40  {
41  public:
42  // pattern assembly flags
43  enum { doPatternVolume = true };
44 
45  // residual assembly flags
46  enum { doAlphaVolume = true };
47 
52  DUNE_DEPRECATED_MSG("Deprecated in Dune-PDELab 2.4, use the local operator ConvectionDiffusionFEM instead!")
53  Laplace (unsigned int quadOrder)
54  : quadOrder_(quadOrder)
55  {}
56 
67  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
68  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
69  {
70  // domain and range field type
71  typedef FiniteElementInterfaceSwitch<
72  typename LFSU::Traits::FiniteElementType
73  > FESwitch;
74  typedef BasisInterfaceSwitch<
75  typename FESwitch::Basis
76  > BasisSwitch;
77  typedef typename BasisSwitch::DomainField DF;
78  typedef typename BasisSwitch::RangeField RF;
79 
80  // dimensions
81  static const int dimLocal = EG::Geometry::mydimension;
82  static const int dimGlobal = EG::Geometry::coorddimension;
83 
84  // select quadrature rule
85  Dune::GeometryType gt = eg.geometry().type();
86  const Dune::QuadratureRule<DF,dimLocal>& rule =
87  Dune::QuadratureRules<DF,dimLocal>::rule(gt,quadOrder_);
88 
89  // loop over quadrature points
90  for(typename Dune::QuadratureRule<DF,dimLocal>::const_iterator it =
91  rule.begin(); it!=rule.end(); ++it)
92  {
93  // evaluate gradient of shape functions
94  // (we assume Galerkin method lfsu=lfsv)
95  std::vector<Dune::FieldMatrix<RF,1,dimGlobal> >
96  gradphiu(lfsu.size());
97  BasisSwitch::gradient(FESwitch::basis(lfsu.finiteElement()),
98  eg.geometry(), it->position(), gradphiu);
99  std::vector<Dune::FieldMatrix<RF,1,dimGlobal> >
100  gradphiv(lfsv.size());
101  BasisSwitch::gradient(FESwitch::basis(lfsv.finiteElement()),
102  eg.geometry(), it->position(), gradphiv);
103 
104  // compute gradient of u
105  Dune::FieldVector<RF,dimGlobal> gradu(0.0);
106  for (size_t i=0; i<lfsu.size(); i++)
107  gradu.axpy(x(lfsu,i),gradphiu[i][0]);
108 
109  // integrate grad u * grad phi_i
110  RF factor = r.weight() * it->weight() * eg.geometry().integrationElement(it->position());
111  for (size_t i=0; i<lfsv.size(); i++)
112  r.rawAccumulate(lfsv,i,(gradu*gradphiv[i][0])*factor);
113  }
114  }
115 
126  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
127  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, M & matrix) const
128  {
129  // Switches between local and global interface
130  typedef FiniteElementInterfaceSwitch<
131  typename LFSU::Traits::FiniteElementType
132  > FESwitch;
133  typedef BasisInterfaceSwitch<
134  typename FESwitch::Basis
135  > BasisSwitch;
136 
137  // domain and range field type
138  typedef typename BasisSwitch::DomainField DF;
139  typedef typename BasisSwitch::RangeField RF;
140  typedef typename LFSU::Traits::SizeType size_type;
141 
142  // dimensions
143  const int dim = EG::Entity::dimension;
144 
145  // select quadrature rule
146  Dune::GeometryType gt = eg.geometry().type();
147  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,quadOrder_);
148 
149  // loop over quadrature points
150  for (typename QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
151  {
152  std::vector<Dune::FieldMatrix<RF,1,dim> > gradphi(lfsu.size());
153  BasisSwitch::gradient(FESwitch::basis(lfsu.finiteElement()),
154  eg.geometry(), it->position(), gradphi);
155 
156  // geometric weight
157  RF factor = it->weight() * eg.geometry().integrationElement(it->position());
158 
159  for (size_type i=0; i<lfsu.size(); i++)
160  {
161  for (size_type j=0; j<lfsv.size(); j++)
162  {
163  // integrate grad u * grad phi
164  matrix.accumulate(lfsv,j,lfsu,i, gradphi[i][0] * gradphi[j][0] * factor);
165  }
166  }
167  }
168  }
169 
170  protected:
171  // Quadrature rule order
172  unsigned int quadOrder_;
173  };
174 
176  } // namespace PDELab
177 } // namespace Dune
178 
179 #endif
sparsity pattern generator
Definition: pattern.hh:13
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 &#39;eg&#39;.
Definition: laplace.hh:127
static const int dim
Definition: adaptivity.hh:83
Definition: adaptivity.hh:27
Default flags for all local operators.
Definition: flags.hh:18
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:68
unsigned int quadOrder_
Definition: laplace.hh:172
Definition: laplace.hh:46
Definition: laplace.hh:37