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