dune-pdelab  2.5-dev
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 
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  using IndexSet = typename GV::IndexSet;
26  using IndexType = typename IndexSet::IndexType;
27 
28  GV gv;
29  V& c;
30  const IndexSet& indexSet;
31 
32 public:
34  using DataType = typename V::value_type;
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  using GV = typename PL::Traits::GridViewType;
100  using IndexSet = typename GV::IndexSet;
101  using DF = typename GV::Grid::ctype;
102  using RF = typename PL::Traits::RangeFieldType;
103  using RangeType = typename PL::Traits::RangeType;
104  enum { dim = PL::Traits::GridViewType::dimension };
105  using Element = typename GV::Traits::template Codim<0>::Entity;
106  using IntersectionIterator = typename GV::IntersectionIterator;
107  using Intersection = typename IntersectionIterator::Intersection;
108  using RT0RangeType = typename Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0>::Traits::LocalBasisType::Traits::RangeType;
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  using RT0Coeffs = Dune::FieldVector<RF,2*dim>;
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  // compute RT0 coefficients for all interior cells
135  for (const auto& element : elements(gv,Dune::Partitions::interior))
136  {
137  // get local cell number
138  int index = is.index(element);
139 
140  // get geometry
141  auto geo = element.geometry();
142 
143  // cell geometry
144  auto ref_el = referenceElement(geo);
145  auto inside_cell_center_local = ref_el.position(0,0);
146  auto inside_cell_center_global = geo.global(inside_cell_center_local);
147 
148  // absolute permeability in primary cell
149  auto tensor_inside = t.A(element,inside_cell_center_local);
150 
151  // pressure evaluation
152  typename PL::Traits::RangeType pl_inside;
153  pl.evaluate(element,inside_cell_center_local,pl_inside);
154 
155  // for coefficient computation
156  RF vn[2*dim]; // normal velocities
157  auto B = geo.jacobianInverseTransposed(inside_cell_center_local); // the transformation. Assume it is linear
158  auto determinant = B.determinant();
159 
160  // loop over cell neighbors
161  for (const auto& intersection : intersections(pl.getGridView(),element))
162  {
163  // set to zero for processor boundary
164  vn[intersection.indexInInside()] = 0.0;
165 
166  // face geometry
167  auto face_local = referenceElement(intersection.geometry()).position(0,0);
168 
169  // interior face
170  if (intersection.neighbor())
171  {
172  auto outside_cell = intersection.outside();
173  auto outside_cell_center_local = referenceElement(outside_cell.geometry()).position(0,0);
174  auto outside_cell_center_global = outside_cell.geometry().global(outside_cell_center_local);
175 
176  // distance of cell centers
177  auto d(outside_cell_center_global);
178  d -= inside_cell_center_global;
179  auto distance = d.two_norm();
180 
181  // absolute permeability
182  auto tensor_outside = t.A(outside_cell,outside_cell_center_local);
183  auto n_F = intersection.centerUnitOuterNormal();
184  Dune::FieldVector<RF,dim> An_F;
185  tensor_inside.mv(n_F,An_F);
186  auto k_inside = n_F*An_F;
187  tensor_outside.mv(n_F,An_F);
188  auto k_outside = n_F*An_F;
189  auto k_avg = 2.0/(1.0/(k_inside+1E-30) + 1.0/(k_outside+1E-30));
190 
191  // pressure evaluation
192  typename PL::Traits::RangeType pl_outside;
193  pl.evaluate(outside_cell,outside_cell_center_local,pl_outside);
194 
195  // set coefficient
196  vn[intersection.indexInInside()] = k_avg*(pl_inside-pl_outside)/distance;
197  }
198 
199  // boundary face
200  if (intersection.boundary())
201  {
202  // distance of cell center to boundary
203  auto d = intersection.geometry().global(face_local);
204  d -= inside_cell_center_global;
205  auto distance = d.two_norm();
206 
207  // evaluate boundary condition type
208  auto bctype = t.bctype(intersection,face_local);
209 
210  // liquid phase Dirichlet boundary
212  {
213  auto iplocal_s = intersection.geometryInInside().global(face_local);
214  auto g_l = t.g(intersection.inside(),iplocal_s);
215  auto n_F = intersection.centerUnitOuterNormal();
216  Dune::FieldVector<RF,dim> An_F;
217  tensor_inside.mv(n_F,An_F);
218  auto k_inside = n_F*An_F;
219  vn[intersection.indexInInside()] = k_inside * (pl_inside-g_l)/distance;
220  }
221 
222  // liquid phase Neumann boundary
224  {
225  auto j = t.j(intersection,face_local);
226  vn[intersection.indexInInside()] = j;
227  }
228  }
229 
230  // compute coefficient
231  auto vstar=intersection.centerUnitOuterNormal(); // normal on tranformef element
232  vstar *= vn[intersection.indexInInside()];
233  Dune::FieldVector<RF,dim> normalhat(0); // normal on reference element
234  if (intersection.indexInInside()%2==0)
235  normalhat[intersection.indexInInside()/2] = -1.0;
236  else
237  normalhat[intersection.indexInInside()/2] = 1.0;
238  Dune::FieldVector<DF,dim> vstarhat(0);
239  B.umtv(vstar,vstarhat); // Piola backward transformation
240  vstarhat *= determinant;
241  storedcoeffs[index][intersection.indexInInside()] = vstarhat*normalhat;
242  }
243  }
244 
245  // communicate coefficients in overlap
246  VectorExchange<GV,std::vector<RT0Coeffs> > dh(gv,storedcoeffs);
247  if (gv.grid().comm().size()>1)
248  gv.grid().communicate(dh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
249  }
250 
251  // set time where operator is to be evaluated (i.e. end of the time intervall)
252  void set_time (typename T::Traits::RangeFieldType time_)
253  {
254  time = time_;
255  }
256 
257  inline void evaluate (const typename Traits::ElementType& e,
258  const typename Traits::DomainType& x,
259  typename Traits::RangeType& y) const
260  {
261  // local cell number
262  int index = is.index(e);
263 
264  // compute velocity on reference element
265  rt0fe.localBasis().evaluateFunction(x,rt0vectors);
266  typename Traits::RangeType yhat(0);
267  for (unsigned int i=0; i<rt0fe.localBasis().size(); i++)
268  yhat.axpy(storedcoeffs[index][i],rt0vectors[i]);
269 
270  // apply Piola transformation
271  if (index != cachedindex)
272  {
273  B = e.geometry().jacobianTransposed(x); // the transformation. Assume it is linear
274  determinant = B.determinant();
275  cachedindex = index;
276  }
277  y = 0;
278  B.umtv(yhat,y);
279  y *= determinant;
280  }
281 
282  inline const typename Traits::GridViewType& getGridView () const
283  {
284  return pl.getGridView();
285  }
286 };
287 
288 #endif // DUNE_PDELAB_LOCALOPERATOR_DARCYCCFV_HH
DarcyVelocityFromHeadCCFV
Provide velocity field for liquid phase.
Definition: darcyccfv.hh:91
DarcyVelocityFromHeadCCFV::DarcyVelocityFromHeadCCFV
DarcyVelocityFromHeadCCFV(const T &t_, const PL &pl_)
Definition: darcyccfv.hh:130
Dune::PDELab::PowerCompositeGridFunctionTraits::ElementType
GV::Traits::template Codim< 0 >::Entity ElementType
codim 0 entity
Definition: function.hh:118
flags.hh
geometrywrapper.hh
Dune::PDELab::PowerCompositeGridFunctionTraits::GridViewType
GV GridViewType
The type of the grid view the function lives on.
Definition: function.hh:115
VectorExchange
Definition: darcyccfv.hh:21
dim
static const int dim
Definition: adaptivity.hh:84
index
std::size_t index
Definition: interpolate.hh:118
DarcyVelocityFromHeadCCFV::evaluate
void evaluate(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeType &y) const
Definition: darcyccfv.hh:257
Dune::PDELab::FunctionTraits< GV::Grid::ctype, GV::dimension, Dune::FieldVector< GV::Grid::ctype, GV::dimension >, RF, m, R >::RangeType
R RangeType
range type
Definition: function.hh:61
defaultimp.hh
idefault.hh
VectorExchange::DataType
typename V::value_type DataType
export type of data for message buffer
Definition: darcyccfv.hh:34
DarcyVelocityFromHeadCCFV::getGridView
const Traits::GridViewType & getGridView() const
Definition: darcyccfv.hh:282
Dune::PDELab::FunctionTraits< GV::Grid::ctype, GV::dimension, Dune::FieldVector< GV::Grid::ctype, GV::dimension >, RF, m, R >::DomainType
Dune::FieldVector< GV::Grid::ctype, GV::dimension > DomainType
domain type in dim-size coordinates
Definition: function.hh:49
VectorExchange::VectorExchange
VectorExchange(const GV &gv_, V &c_)
constructor
Definition: darcyccfv.hh:37
e
const Entity & e
Definition: localfunctionspace.hh:120
VectorExchange::fixedsize
bool fixedsize(int dim, int codim) const
returns true if size per entity of given dim and codim is a constant
Definition: darcyccfv.hh:48
function.hh
Dune::PDELab::GridFunctionBase
leaf of a function tree
Definition: function.hh:298
VectorExchange::size
size_t size(EntityType &e) const
Definition: darcyccfv.hh:58
Dune::PDELab::ConvectionDiffusionBoundaryConditions::Type
Type
Definition: convectiondiffusionparameter.hh:65
pattern.hh
VectorExchange::contains
bool contains(int dim, int codim) const
returns true if data for this codim should be communicated
Definition: darcyccfv.hh:42
VectorExchange::scatter
void scatter(MessageBuffer &buff, const EntityType &e, size_t n)
Definition: darcyccfv.hh:75
Dune::PDELab::GridFunctionTraits
traits class holding the function signature, same as in local function
Definition: function.hh:176
VectorExchange::gather
void gather(MessageBuffer &buff, const EntityType &e) const
pack data from user to message buffer
Definition: darcyccfv.hh:65
Dune::PDELab::ConvectionDiffusionBoundaryConditions::Dirichlet
@ Dirichlet
Definition: convectiondiffusionparameter.hh:65
convectiondiffusionparameter.hh
DarcyVelocityFromHeadCCFV::set_time
void set_time(typename T::Traits::RangeFieldType time_)
Definition: darcyccfv.hh:252
Dune::PDELab::ConvectionDiffusionBoundaryConditions::Neumann
@ Neumann
Definition: convectiondiffusionparameter.hh:65