dune-pdelab  2.0.0
diffusionmfd.hh
Go to the documentation of this file.
1 #ifndef DUNE_PDELAB_DIFFUSIONMFD_HH
2 #define DUNE_PDELAB_DIFFUSIONMFD_HH
3 
4 #include "pattern.hh"
5 #include "flags.hh"
6 #include "mfdcommon.hh"
7 #include "diffusionparam.hh"
8 
9 namespace Dune
10 {
11  namespace PDELab
12  {
13 
17 
30  template<class Data, class WBuilder = MimeticBrezziW<typename Data::ctype,Data::dimension> >
32  : public FullVolumePattern
34  {
35  static const unsigned int dim = Data::dimension;
36  typedef typename Data::ctype ctype;
37  typedef typename Data::rtype rtype;
38 
39  public:
40  // pattern assembly flags
41  enum { doPatternVolume = true };
42 
43  // residual assembly flags
44  enum { doAlphaVolume = true };
45  enum { doSkeletonTwoSided = true };
46  enum { doAlphaSkeleton = true };
47  enum { doAlphaBoundary = true };
48  enum { doAlphaVolumePostSkeleton = true };
49 
50  enum { doLambdaVolume = true };
51  enum { doLambdaBoundary = true };
52 
53  DiffusionMFD(const Data& data_, const WBuilder& wbuilder_ = WBuilder())
54  : data(data_), wbuilder(wbuilder_)
55  {}
56 
57  // alpha ***********************************************
58 
59  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
60  void alpha_volume(const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
61  {
62  cell.init(eg.entity());
63  }
64 
65  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
66  void alpha_skeleton(const IG& ig,
67  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
68  const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
69  R& r_s, R& r_n) const
70  {
71  cell.add_face(ig.intersection());
72  }
73 
74  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
75  void alpha_boundary(const IG& ig, const LFSU& lfsu_s, const X& x_s,
76  const LFSV& lfsv_s, R& r_s) const
77  {
78  cell.add_face(ig.intersection());
79  }
80 
81  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
82  void alpha_volume_post_skeleton(const EG& eg, const LFSU& lfsu, const X& x,
83  const LFSV& lfsv, R& r) const
84  {
85  // extract subspaces
86  typedef typename LFSU::template Child<0>::Type CellUnknowns;
87  const CellUnknowns& cell_space = lfsu.template child<0>();
88  typedef typename LFSU::template Child<1>::Type FaceUnknowns;
89  const FaceUnknowns& face_space = lfsu.template child<1>();
90 
91  // get permeability for current cell
92  GeometryType gt = eg.geometry().type();
93  FieldVector<ctype,dim> localcenter = ReferenceElements<ctype,dim>::general(gt).position(0,0);
94  FieldMatrix<rtype,dim,dim> K = data.K(eg.entity(), localcenter);
95 
96  // build matrix W
97  wbuilder.build_W(cell, K, W);
98 
99  // Compute residual
100  for(unsigned int e = 0, m = 0; e < cell.num_faces; ++e)
101  for(unsigned int f = 0; f < cell.num_faces; ++f, ++m)
102  {
103  r.accumulate(cell_space,0,W[m]
104  * (x(cell_space,0) - x(face_space,f)));
105  r.accumulate(face_space,e,-W[m]
106  * (x(cell_space,0) - x(face_space,f)));
107  }
108 
109  // Helmholtz term
110  r.accumulate(cell_space,0,data.a_0(eg.entity(), localcenter)
111  * x(cell_space,0));
112  }
113 
114  // jacobian ********************************************
115 
116  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
117  void jacobian_volume(const EG& eg, const LFSU& lfsu, const X& x,
118  const LFSV& lfsv, M& mat) const
119  {
120  cell.init(eg.entity());
121  }
122 
123  template<typename IG, typename LFSU, typename X, typename LFSV, typename M>
124  void jacobian_skeleton(const IG& ig,
125  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
126  const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
127  M& mat_ss, M& mat_sn,
128  M& mat_ns, M& mat_nn) const
129  {
130  cell.add_face(ig.intersection());
131  }
132 
133  template<typename IG, typename LFSU, typename X, typename LFSV, typename M>
134  void jacobian_boundary(const IG& ig, const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
135  M& mat_ss) const
136  {
137  cell.add_face(ig.intersection());
138  }
139 
140  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
141  void jacobian_volume_post_skeleton(const EG& eg, const LFSU& lfsu, const X& x,
142  const LFSV& lfsv, M& mat) const
143  {
144  // extract subspaces
145  typedef typename LFSU::template Child<0>::Type CellUnknowns;
146  const CellUnknowns& cell_space = lfsu.template child<0>();
147  typedef typename LFSU::template Child<1>::Type FaceUnknowns;
148  const FaceUnknowns& face_space = lfsu.template child<1>();
149 
150  // get permeability for current cell
151  GeometryType gt = eg.geometry().type();
152  FieldVector<ctype,dim> localcenter = ReferenceElements<ctype,dim>::general(gt).position(0,0);
153  FieldMatrix<rtype,dim,dim> K = data.K(eg.entity(), localcenter);
154 
155  // build matrix W
156  wbuilder.build_W(cell, K, W);
157 
158  // Compute residual
159  for(unsigned int e = 0, m = 0; e < cell.num_faces; ++e)
160  for(unsigned int f = 0; f < cell.num_faces; ++f, ++m)
161  {
162  mat.accumulate(cell_space,0,cell_space,0,W[m]);
163  mat.accumulate(cell_space,0,face_space,f,-W[m]);
164  mat.accumulate(face_space,f,cell_space,0,-W[m]);
165  mat.accumulate(face_space,f,face_space,f,W[m]);
166  }
167 
168  // Helmholtz term
169  mat.accumulate(cell_space,0,cell_space,0,
170  data.a_0(eg.entity(), localcenter));
171  }
172 
173  // jacobian_apply **************************************
174 
175  template<typename EG, typename LFSU, typename X, typename LFSV, typename Y>
176  void jacobian_apply_volume(const EG& eg, const LFSU& lfsu, const X& x,
177  const LFSV& lfsv, Y& y) const
178  {
179  cell.init(eg.entity());
180  }
181 
182  template<typename IG, typename LFSU, typename X, typename LFSV, typename Y>
183  void jacobian_apply_skeleton(const IG& ig,
184  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
185  const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
186  Y& y_s, Y& y_n) const
187  {
188  cell.add_face(ig.intersection());
189  }
190 
191  template<typename IG, typename LFSU, typename X, typename LFSV, typename Y>
192  void jacobian_apply_boundary(const IG& ig, const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
193  Y& y_s) const
194  {
195  cell.add_face(ig.intersection());
196  }
197 
198  template<typename EG, typename LFSU, typename X, typename LFSV, typename Y>
199  void jacobian_apply_volume_post_skeleton(const EG& eg, const LFSU& lfsu, const X& x,
200  const LFSV& lfsv, Y& y) const
201  {
202  // extract subspaces
203  typedef typename LFSU::template Child<0>::Type CellUnknowns;
204  const CellUnknowns& cell_space = lfsu.template child<0>();
205  typedef typename LFSU::template Child<1>::Type FaceUnknowns;
206  const FaceUnknowns& face_space = lfsu.template child<1>();
207 
208  // get permeability for current cell
209  GeometryType gt = eg.geometry().type();
210  FieldVector<ctype,dim> localcenter = ReferenceElements<ctype,dim>::general(gt).position(0,0);
211  FieldMatrix<rtype,dim,dim> K = data.K(eg.entity(), localcenter);
212 
213  // build matrix W
214  wbuilder.build_W(cell, K, W);
215 
216  // Compute residual
217  for(int e = 0, m = 0; e < cell.num_faces; ++e)
218  for(int f = 0; f < cell.num_faces; ++f, ++m)
219  {
220  y[cell_space.localIndex(0)] += W[m] * x[cell_space.localIndex(0)];
221  y[cell_space.localIndex(0)] -= W[m] * x[face_space.localIndex(f)];
222  y[face_space.localIndex(f)] -= W[m] * x[cell_space.localIndex(0)];
223  y[face_space.localIndex(f)] += W[m] * x[face_space.localIndex(f)];
224  }
225 
226  // Helmholtz term
227  y[cell_space.localIndex(0)] += data.a_0(eg.entity(), localcenter)
228  * x[cell_space.localIndex(0)];
229  }
230 
231  // lambda **********************************************
232 
233  template<typename EG, typename LFSV, typename R>
234  void lambda_volume(const EG& eg, const LFSV& lfsv, R& r) const
235  {
236  // extract subspaces
237  typedef typename LFSV::template Child<0>::Type CellUnknowns;
238  const CellUnknowns& cell_space = lfsv.template child<0>();
239 
240  GeometryType gt = eg.geometry().type();
241  FieldVector<ctype,dim> localcenter = ReferenceElements<ctype,dim>::general(gt).position(0,0);
242  r.accumulate(cell_space,0,-cell.volume * data.f(eg.entity(), localcenter));
243  }
244 
245  template<typename IG, typename LFSV, typename R>
246  void lambda_boundary(const IG& ig, const LFSV& lfsv, R& r) const
247  {
248  // local index of current face
249  unsigned int e = cell.num_faces - 1;
250 
251  GeometryType gt = ig.intersection().type();
252  FieldVector<ctype,dim-1> center = ReferenceElements<ctype,dim-1>::general(gt).position(0,0);
253  FieldVector<ctype,dim> local_face_center = ig.geometryInInside().global(center);
254 
255  if (data.bcType(ig, center) == Data::bcNeumann)
256  {
257  typedef typename LFSV::template Child<1>::Type FaceUnknowns;
258  const FaceUnknowns& face_space = lfsv.template child<1>();
259 
260  r.accumulate(face_space,e,cell.face_areas[e]
261  * data.j(*(ig.inside()), local_face_center));
262  }
263  }
264 
265  private:
266  const Data& data;
267  const WBuilder wbuilder;
269  mutable std::vector<rtype> W;
270  };
271 
273  }
274 }
275 
276 #endif
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: diffusionmfd.hh:75
Default flags for all local operators.
Definition: flags.hh:18
void jacobian_apply_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, Y &y_s) const
Definition: diffusionmfd.hh:192
void jacobian_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, M &mat_ss, M &mat_sn, M &mat_ns, M &mat_nn) const
Definition: diffusionmfd.hh:124
void jacobian_apply_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, Y &y) const
Definition: diffusionmfd.hh:176
Definition: diffusionmfd.hh:41
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: diffusionmfd.hh:246
void jacobian_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, M &mat_ss) const
Definition: diffusionmfd.hh:134
void jacobian_apply_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, Y &y_s, Y &y_n) const
Definition: diffusionmfd.hh:183
Definition: diffusionmfd.hh:44
void jacobian_apply_volume_post_skeleton(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, Y &y) const
Definition: diffusionmfd.hh:199
Definition: diffusionmfd.hh:50
Definition: diffusionmfd.hh:51
void jacobian_volume_post_skeleton(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: diffusionmfd.hh:141
Definition: diffusionmfd.hh:31
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: diffusionmfd.hh:117
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: diffusionmfd.hh:60
Definition: mfdcommon.hh:15
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: diffusionmfd.hh:234
const IG & ig
Definition: common/constraints.hh:146
void alpha_volume_post_skeleton(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: diffusionmfd.hh:82
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: diffusionmfd.hh:66
Definition: diffusionmfd.hh:47
Definition: diffusionmfd.hh:46
DiffusionMFD(const Data &data_, const WBuilder &wbuilder_=WBuilder())
Definition: diffusionmfd.hh:53
const F & f
Definition: common/constraints.hh:145
sparsity pattern generator
Definition: pattern.hh:14
const EG & eg
Definition: common/constraints.hh:277
const E & e
Definition: interpolate.hh:172