dune-pdelab  2.0.0
laplacedirichletccfv.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_LAPLACEDIRICHLETCCFV_HH
3 #define DUNE_PDELAB_LAPLACEDIRICHLETCCFV_HH
4 
5 #include<dune/common/exceptions.hh>
6 #include<dune/common/fvector.hh>
7 #include<dune/common/static_assert.hh>
8 #include<dune/geometry/referenceelements.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 cell centered finite volumes on axiparallel cube grids
23  // G : grid function for Dirichlet boundary conditions
24  template<typename G>
25  class LaplaceDirichletCCFV : public NumericalJacobianApplySkeleton<LaplaceDirichletCCFV<G> >,
26  public NumericalJacobianApplyBoundary<LaplaceDirichletCCFV<G> >,
27  public NumericalJacobianSkeleton<LaplaceDirichletCCFV<G> >,
28  public NumericalJacobianBoundary<LaplaceDirichletCCFV<G> >,
29  public FullSkeletonPattern,
30  public FullVolumePattern,
32  {
33  public:
34  // pattern assembly flags
35  enum { doPatternVolume = true };
36  enum { doPatternSkeleton = true };
37 
38  // residual assembly flags
39  enum { doAlphaSkeleton = true };
40  enum { doAlphaBoundary = true };
41 
42  LaplaceDirichletCCFV (const G& g_) : g(g_) {}
43 
44  // skeleton integral depending on test and ansatz functions
45  // each face is only visited ONCE!
46  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
47  void alpha_skeleton (const IG& ig,
48  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
49  const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
50  R& r_s, R& r_n) const
51  {
52  // domain and range field type
53  typedef typename LFSU::Traits::FiniteElementType::
54  Traits::LocalBasisType::Traits::DomainFieldType DF;
55  typedef typename LFSU::Traits::FiniteElementType::
56  Traits::LocalBasisType::Traits::RangeFieldType RF;
57 
58  // center in face's reference element
59  const Dune::FieldVector<DF,IG::dimension-1>&
60  face_local = Dune::ReferenceElements<DF,IG::dimension-1>::general(ig.geometry().type()).position(0,0);
61 
62  // face volume for integration
63  RF face_volume = ig.geometry().integrationElement(face_local)
64  *Dune::ReferenceElements<DF,IG::dimension-1>::general(ig.geometry().type()).volume();
65 
66  // cell centers in references elements
67  const Dune::FieldVector<DF,IG::dimension>&
68  inside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
69  const Dune::FieldVector<DF,IG::dimension>&
70  outside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.outside()->type()).position(0,0);
71 
72  // distance between the two cell centers
73  Dune::FieldVector<DF,IG::dimension>
74  inside_global = ig.inside()->geometry().global(inside_local);
75  Dune::FieldVector<DF,IG::dimension>
76  outside_global = ig.outside()->geometry().global(outside_local);
77  inside_global -= outside_global;
78  RF distance = inside_global.two_norm();
79 
80  // contribution to residual on inside element, other residual is computed by symmetric call
81  r_s.accumulate(lfsu_s,0,(x_s(lfsu_s,0)-x_n(lfsu_n,0))*face_volume/distance);
82  r_n.accumulate(lfsu_n,0,-(x_s(lfsu_s,0)-x_n(lfsu_n,0))*face_volume/distance);
83  }
84 
85  // skeleton integral depending on test and ansatz functions
86  // We put the Dirchlet evaluation also in the alpha term to save som e geometry evaluations
87  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
88  void alpha_boundary (const IG& ig,
89  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
90  R& r_s) const
91  {
92  // domain and range field type
93  typedef typename LFSU::Traits::FiniteElementType::
94  Traits::LocalBasisType::Traits::DomainFieldType DF;
95  typedef typename LFSU::Traits::FiniteElementType::
96  Traits::LocalBasisType::Traits::RangeFieldType RF;
97 
98  // center in face's reference element
99  const Dune::FieldVector<DF,IG::dimension-1>&
100  face_local = Dune::ReferenceElements<DF,IG::dimension-1>::general(ig.geometry().type()).position(0,0);
101 
102  // face volume for integration
103  RF face_volume = ig.geometry().integrationElement(face_local)
104  *Dune::ReferenceElements<DF,IG::dimension-1>::general(ig.geometry().type()).volume();
105 
106  // cell center in reference element
107  const Dune::FieldVector<DF,IG::dimension>&
108  inside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
109 
110  // distance between cell center and face center
111  Dune::FieldVector<DF,IG::dimension>
112  inside_global = ig.inside()->geometry().global(inside_local);
113  Dune::FieldVector<DF,IG::dimension>
114  outside_global = ig.geometry().global(face_local);
115  inside_global -= outside_global;
116  RF distance = inside_global.two_norm();
117 
118  // evaluate boundary condition function
119  typename G::Traits::DomainType x = ig.geometryInInside().global(face_local);
120  typename G::Traits::RangeType y;
121  g.evaluate(*(ig.inside()),x,y);
122 
123  // contribution to residual on inside element
124  r_s.accumulate(lfsu_s,0,(x_s(lfsu_s,0)-y[0])*face_volume/distance);
125  }
126 
127  private:
128  const G& g;
129  };
130 
132  } // namespace PDELab
133 } // namespace Dune
134 
135 #endif
Definition: laplacedirichletccfv.hh:40
sparsity pattern generator
Definition: pattern.hh:30
Implement jacobian_apply_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:445
Default flags for all local operators.
Definition: flags.hh:18
Implement jacobian_boundary() based on alpha_boundary()
Definition: defaultimp.hh:251
Definition: laplacedirichletccfv.hh:39
Implement jacobian_apply_boundary() based on alpha_boundary()
Definition: defaultimp.hh:538
const IG & ig
Definition: common/constraints.hh:146
void alpha_skeleton(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &x_n, const LFSV &lfsv_n, R &r_s, R &r_n) const
Definition: laplacedirichletccfv.hh:47
sparsity pattern generator
Definition: pattern.hh:14
Definition: laplacedirichletccfv.hh:25
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:157
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: laplacedirichletccfv.hh:88
Definition: laplacedirichletccfv.hh:35
LaplaceDirichletCCFV(const G &g_)
Definition: laplacedirichletccfv.hh:42
Definition: laplacedirichletccfv.hh:36