dune-pdelab  2.4.1
darcyccfv.hh
Go to the documentation of this file.
1 // -*- tab-width: 2; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_LOCALOPERATOR_DARCYCCFV_HH
3 #define DUNE_PDELAB_LOCALOPERATOR_DARCYCCFV_HH
4 
5 #include<dune/common/exceptions.hh>
6 #include<dune/common/fvector.hh>
7 #include<dune/common/typetraits.hh>
8 #include<dune/geometry/referenceelements.hh>
9 #include<dune/localfunctions/raviartthomas/raviartthomascube.hh>
10 
19 
20 // A DataHandle class to exchange RT0 coefficients in the overlap
21 template<class GV, class V> // mapper type and vector type
23  : public Dune::CommDataHandleIF<VectorExchange<GV,V>,
24  typename V::value_type>
25 {
26  using IndexSet = typename GV::IndexSet;
27  using IndexType = typename IndexSet::IndexType;
28 
29  GV gv;
30  V& c;
31  const IndexSet& indexSet;
32 
33 public:
35  using DataType = typename V::value_type;
36 
38  VectorExchange (const GV& gv_, V& c_)
39  : gv(gv_), c(c_), indexSet(gv.indexSet())
40  {}
41 
43  bool contains (int dim, int codim) const
44  {
45  return (codim==0);
46  }
47 
49  bool fixedsize (int dim, int codim) const
50  {
51  return true;
52  }
53 
58  template<class EntityType>
59  size_t size (EntityType& e) const
60  {
61  return 1;
62  }
63 
65  template<class MessageBuffer, class EntityType>
66  void gather (MessageBuffer& buff, const EntityType& e) const
67  {
68  buff.write(c[indexSet.index(e)]);
69  }
70 
75  template<class MessageBuffer, class EntityType>
76  void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
77  {
78  DataType x;
79  buff.read(x);
80  c[indexSet.index(e)]=x;
81  }
82 };
83 
91 template<typename T, typename PL>
93  : public Dune::PDELab::GridFunctionBase<Dune::PDELab::GridFunctionTraits<typename PL::Traits::GridViewType,
94  typename PL::Traits::RangeFieldType,
95  PL::Traits::GridViewType::dimension,
96  Dune::FieldVector<typename PL::Traits::RangeFieldType,PL::Traits::GridViewType::dimension> >,
97  DarcyVelocityFromHeadCCFV<T,PL> >
98 {
99  // extract useful types
100  using GV = typename PL::Traits::GridViewType;
101  using IndexSet = typename GV::IndexSet;
102  using DF = typename GV::Grid::ctype;
103  using RF = typename PL::Traits::RangeFieldType;
104  using RangeType = typename PL::Traits::RangeType;
105  enum { dim = PL::Traits::GridViewType::dimension };
106  using Element = typename GV::Traits::template Codim<0>::Entity;
107  using IntersectionIterator = typename GV::IntersectionIterator;
108  using Intersection = typename IntersectionIterator::Intersection;
109  using RT0RangeType = typename Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0>::Traits::LocalBasisType::Traits::RangeType;
111 
112  const T& t;
113  const PL& pl;
114  Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0> rt0fe;
115 
116  mutable Dune::FieldMatrix<DF,dim,dim> B;
117  mutable RF determinant;
118  mutable int cachedindex;
119  typename T::Traits::RangeFieldType time;
120 
121  using RT0Coeffs = Dune::FieldVector<RF,2*dim>;
122  GV gv;
123  const IndexSet& is;
124  std::vector<RT0Coeffs> storedcoeffs;
125  mutable std::vector<RT0RangeType> rt0vectors;
126 
127 public:
130 
131  DarcyVelocityFromHeadCCFV (const T& t_, const PL& pl_)
132  : t(t_), pl(pl_), cachedindex(-1), time(0), gv(pl_.getGridView()), is(gv.indexSet()), storedcoeffs(is.size(0)),
133  rt0vectors(rt0fe.localBasis().size())
134  {
135  // compute RT0 coefficients for all interior cells
136  for (const auto& element : elements(gv,Dune::Partitions::interior))
137  {
138  // get local cell number
139  int index = is.index(element);
140 
141  // get geometry
142  auto geo = element.geometry();
143 
144  // cell geometry
145  auto ref_el = referenceElement(geo);
146  auto inside_cell_center_local = ref_el.position(0,0);
147  auto inside_cell_center_global = geo.global(inside_cell_center_local);
148 
149  // absolute permeability in primary cell
150  auto tensor_inside = t.A(element,inside_cell_center_local);
151 
152  // pressure evaluation
153  typename PL::Traits::RangeType pl_inside;
154  pl.evaluate(element,inside_cell_center_local,pl_inside);
155 
156  // for coefficient computation
157  RF vn[2*dim]; // normal velocities
158  auto B = geo.jacobianInverseTransposed(inside_cell_center_local); // the transformation. Assume it is linear
159  auto determinant = B.determinant();
160 
161  // loop over cell neighbors
162  for (const auto& intersection : intersections(pl.getGridView(),element))
163  {
164  // set to zero for processor boundary
165  vn[intersection.indexInInside()] = 0.0;
166 
167  // face geometry
168  auto face_local = referenceElement(intersection.geometry()).position(0,0);
169 
170  // interior face
171  if (intersection.neighbor())
172  {
173  auto outside_cell = intersection.outside();
174  auto outside_cell_center_local = referenceElement(outside_cell.geometry()).position(0,0);
175  auto outside_cell_center_global = outside_cell.geometry().global(outside_cell_center_local);
176 
177  // distance of cell centers
178  auto d(outside_cell_center_global);
179  d -= inside_cell_center_global;
180  auto distance = d.two_norm();
181 
182  // absolute permeability
183  auto tensor_outside = t.A(outside_cell,outside_cell_center_local);
184  auto n_F = intersection.centerUnitOuterNormal();
185  Dune::FieldVector<RF,dim> An_F;
186  tensor_inside.mv(n_F,An_F);
187  auto k_inside = n_F*An_F;
188  tensor_outside.mv(n_F,An_F);
189  auto k_outside = n_F*An_F;
190  auto k_avg = 2.0/(1.0/(k_inside+1E-30) + 1.0/(k_outside+1E-30));
191 
192  // pressure evaluation
193  typename PL::Traits::RangeType pl_outside;
194  pl.evaluate(outside_cell,outside_cell_center_local,pl_outside);
195 
196  // set coefficient
197  vn[intersection.indexInInside()] = k_avg*(pl_inside-pl_outside)/distance;
198  }
199 
200  // boundary face
201  if (intersection.boundary())
202  {
203  // distance of cell center to boundary
204  auto d = intersection.geometry().global(face_local);
205  d -= inside_cell_center_global;
206  auto distance = d.two_norm();
207 
208  // evaluate boundary condition type
209  auto bctype = t.bctype(intersection,face_local);
210 
211  // liquid phase Dirichlet boundary
213  {
214  auto iplocal_s = intersection.geometryInInside().global(face_local);
215  auto g_l = t.g(intersection.inside(),iplocal_s);
216  auto n_F = intersection.centerUnitOuterNormal();
217  Dune::FieldVector<RF,dim> An_F;
218  tensor_inside.mv(n_F,An_F);
219  auto k_inside = n_F*An_F;
220  vn[intersection.indexInInside()] = k_inside * (pl_inside-g_l)/distance;
221  }
222 
223  // liquid phase Neumann boundary
225  {
226  auto j = t.j(intersection,face_local);
227  vn[intersection.indexInInside()] = j;
228  }
229  }
230 
231  // compute coefficient
232  auto vstar=intersection.centerUnitOuterNormal(); // normal on tranformef element
233  vstar *= vn[intersection.indexInInside()];
234  Dune::FieldVector<RF,dim> normalhat(0); // normal on reference element
235  if (intersection.indexInInside()%2==0)
236  normalhat[intersection.indexInInside()/2] = -1.0;
237  else
238  normalhat[intersection.indexInInside()/2] = 1.0;
239  Dune::FieldVector<DF,dim> vstarhat(0);
240  B.umtv(vstar,vstarhat); // Piola backward transformation
241  vstarhat *= determinant;
242  storedcoeffs[index][intersection.indexInInside()] = vstarhat*normalhat;
243  }
244  }
245 
246  // communicate coefficients in overlap
247  VectorExchange<GV,std::vector<RT0Coeffs> > dh(gv,storedcoeffs);
248  if (gv.grid().comm().size()>1)
249  gv.grid().communicate(dh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
250  }
251 
252  // set time where operator is to be evaluated (i.e. end of the time intervall)
253  void set_time (typename T::Traits::RangeFieldType time_)
254  {
255  time = time_;
256  }
257 
258  inline void evaluate (const typename Traits::ElementType& e,
259  const typename Traits::DomainType& x,
260  typename Traits::RangeType& y) const
261  {
262  // local cell number
263  int index = is.index(e);
264 
265  // compute velocity on reference element
266  rt0fe.localBasis().evaluateFunction(x,rt0vectors);
267  typename Traits::RangeType yhat(0);
268  for (unsigned int i=0; i<rt0fe.localBasis().size(); i++)
269  yhat.axpy(storedcoeffs[index][i],rt0vectors[i]);
270 
271  // apply Piola transformation
272  if (index != cachedindex)
273  {
274  B = e.geometry().jacobianTransposed(x); // the transformation. Assume it is linear
275  determinant = B.determinant();
276  cachedindex = index;
277  }
278  y = 0;
279  B.umtv(yhat,y);
280  y *= determinant;
281  }
282 
283  inline const typename Traits::GridViewType& getGridView () const
284  {
285  return pl.getGridView();
286  }
287 };
288 
289 #endif
size_t size(EntityType &e) const
Definition: darcyccfv.hh:59
Definition: darcyccfv.hh:22
Provide velocity field for liquid phase.
Definition: darcyccfv.hh:92
void gather(MessageBuffer &buff, const EntityType &e) const
pack data from user to message buffer
Definition: darcyccfv.hh:66
bool contains(int dim, int codim) const
returns true if data for this codim should be communicated
Definition: darcyccfv.hh:43
static const int dim
Definition: adaptivity.hh:83
VectorExchange(const GV &gv_, V &c_)
constructor
Definition: darcyccfv.hh:38
GV GridViewType
The type of the grid view the function lives on.
Definition: function.hh:114
ReferenceElementWrapper< ReferenceElement< typename Geometry::ctype, Geometry::mydimension > > referenceElement(const Geometry &geo)
Returns the reference element for the given geometry.
Definition: referenceelements.hh:144
GV::Traits::template Codim< 0 >::Entity ElementType
codim 0 entity
Definition: function.hh:117
typename V::value_type DataType
export type of data for message buffer
Definition: darcyccfv.hh:35
Definition: convectiondiffusionparameter.hh:65
Type
Definition: convectiondiffusionparameter.hh:65
const Entity & e
Definition: localfunctionspace.hh:111
Dune::FieldVector< GV::Grid::ctype, GV::dimension > DomainType
domain type in dim-size coordinates
Definition: function.hh:48
Definition: convectiondiffusionparameter.hh:65
void scatter(MessageBuffer &buff, const EntityType &e, size_t n)
Definition: darcyccfv.hh:76
void evaluate(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeType &y) const
Definition: darcyccfv.hh:258
bool fixedsize(int dim, int codim) const
returns true if size per entity of given dim and codim is a constant
Definition: darcyccfv.hh:49
const Traits::GridViewType & getGridView() const
Definition: darcyccfv.hh:283
void set_time(typename T::Traits::RangeFieldType time_)
Definition: darcyccfv.hh:253
DarcyVelocityFromHeadCCFV(const T &t_, const PL &pl_)
Definition: darcyccfv.hh:131
leaf of a function tree
Definition: function.hh:576
traits class holding the function signature, same as in local function
Definition: function.hh:175