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