dune-pdelab  2.0.0
laplacedirichletp12d.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_LAPLACEDIRICHLETP12D_HH
3 #define DUNE_PDELAB_LAPLACEDIRICHLETP12D_HH
4 
5 #include <dune/common/exceptions.hh>
6 #include <dune/common/fvector.hh>
7 #include <dune/common/fmatrix.hh>
8 #include <dune/common/static_assert.hh>
9 
11 
12 #include"pattern.hh"
13 #include"flags.hh"
14 
15 
16 namespace Dune {
17  namespace PDELab {
18 
19  // a local operator for solving the Laplace equation with Dirichlet boundary conditions
20  // - \Delta u = 0 in \Omega,
21  // u = g on \partial\Omega
22  // with P1 conforming finite elements on triangles
23  class LaplaceDirichletP12D : public NumericalJacobianApplyVolume<LaplaceDirichletP12D>,
24  public NumericalJacobianVolume<LaplaceDirichletP12D>,
25  public FullVolumePattern,
27  {
28  public:
29  // pattern assembly flags
30  enum { doPatternVolume = true };
31 
32  // residual assembly flags
33  enum { doAlphaVolume = true };
34 
35  // volume integral depending on test and ansatz functions
36  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
37  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
38  {
39  // domain and range field type
40  typedef typename LFSU::Traits::FiniteElementType::
41  Traits::LocalBasisType::Traits::DomainFieldType DF;
42  typedef typename LFSU::Traits::FiniteElementType::
43  Traits::LocalBasisType::Traits::RangeFieldType RF;
44 
45  // define integration point (hard coded quadrature)
46  Dune::FieldVector<DF,2> integrationpoint(1.0/3.0);
47 
48  // gradient of shape functions at integration point
49  typedef typename LFSU::Traits::FiniteElementType::
50  Traits::LocalBasisType::Traits::JacobianType JT;
51  std::vector<JT> gradients(lfsu.size());
52  lfsu.finiteElement().localBasis().
53  evaluateJacobian(integrationpoint,gradients);
54 
55  // transformation of gradients to real element
56  const typename EG::Geometry::JacobianInverseTransposed
57  jac = eg.geometry().jacobianInverseTransposed(integrationpoint);
58  Dune::FieldVector<RF,2> gradphi[3];
59  for (int i=0; i<3; i++)
60  {
61  gradphi[i] = 0.0;
62  jac.umv(gradients[i][0],gradphi[i]);
63  }
64 
65  // compute gradient of solution at integration point
66  Dune::FieldVector<RF,2> gradu(0.0);
67  for (int i=0; i<3; i++)
68  gradu.axpy(x[i],gradphi[i]);
69 
70  // integrate grad u * grad phi_i (0.5 is the area of the reference element)
71  RF area = 0.5*eg.geometry().integrationElement(integrationpoint);
72  for (int i=0; i<3; i++)
73  r.accumulate(lfsv, i, (gradu*gradphi[i])*area);
74  }
75  };
76 
78  } // namespace PDELab
79 } // namespace Dune
80 
81 #endif
Default flags for all local operators.
Definition: flags.hh:18
Definition: laplacedirichletp12d.hh:23
Implement jacobian_volume() based on alpha_volume()
Definition: defaultimp.hh:32
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: laplacedirichletp12d.hh:37
Definition: laplacedirichletp12d.hh:33
sparsity pattern generator
Definition: pattern.hh:14
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
const EG & eg
Definition: common/constraints.hh:277
Definition: laplacedirichletp12d.hh:30