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