2 #ifndef DarcyVelocityFromHeadCCFV_HH
3 #define DarcyVelocityFromHeadCCFV_HH
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>
20 template<
class GV,
class V>
22 :
public Dune::CommDataHandleIF<VectorExchange<GV,V>,
23 typename V::value_type>
25 typedef typename GV::IndexSet IndexSet;
26 typedef typename IndexSet::IndexType IndexType;
30 const IndexSet& indexSet;
38 : gv(gv_), c(c_), indexSet(gv.indexSet())
57 template<
class EntityType>
58 size_t size (EntityType&
e)
const
64 template<
class MessageBuffer,
class EntityType>
65 void gather (MessageBuffer& buff,
const EntityType&
e)
const
67 buff.write(c[indexSet.index(e)]);
74 template<
class MessageBuffer,
class EntityType>
75 void scatter (MessageBuffer& buff,
const EntityType&
e,
size_t n)
79 c[indexSet.index(e)]=x;
90 template<
typename T,
typename PL>
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> >
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;
113 Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0> rt0fe;
115 mutable Dune::FieldMatrix<DF,dim,dim> B;
116 mutable RF determinant;
117 mutable int cachedindex;
118 typename T::Traits::RangeFieldType time;
120 typedef Dune::FieldVector<RF,2*dim> RT0Coeffs;
123 std::vector<RT0Coeffs> storedcoeffs;
124 mutable std::vector<RT0RangeType> rt0vectors;
131 : t(t_), pl(pl_), cachedindex(-1), time(0), gv(pl_.
getGridView()), is(gv.indexSet()), storedcoeffs(is.size(0)),
132 rt0vectors(rt0fe.localBasis().size())
135 typedef typename GV::Traits::template Codim<0>::template Partition<Dune::Interior_Partition>::Iterator ElementIterator;
139 for (ElementIterator it = gv.template begin<0,Dune::Interior_Partition>();
140 it!=gv.template end<0,Dune::Interior_Partition>(); ++it)
143 int index = is.index(*it);
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);
153 typename T::Traits::PermTensorType tensor_inside;
154 tensor_inside = t.A(*it,inside_cell_center_local);
157 typename PL::Traits::RangeType pl_inside;
158 pl.evaluate(*it,inside_cell_center_local,pl_inside);
162 Dune::FieldMatrix<typename Traits::DomainFieldType,dim,dim>
163 B = it->geometry().jacobianInverseTransposed(inside_cell_center_local);
164 RF determinant = B.determinant();
167 IntersectionIterator endit = pl.getGridView().iend(*it);
168 for (IntersectionIterator iit = pl.getGridView().ibegin(*it); iit!=endit; ++iit)
171 vn[iit->indexInInside()] = 0.0;
174 const Dune::FieldVector<DF,dim-1>&
175 face_local = Dune::ReferenceElements<DF,dim-1>::general(iit->geometry().type()).position(0,0);
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);
188 Dune::FieldVector<DF,dim> d(outside_cell_center_global);
189 d -= inside_cell_center_global;
190 RF distance = d.two_norm();
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));
204 typename PL::Traits::RangeType pl_outside;
205 pl.evaluate(outside_cell,outside_cell_center_local,pl_outside);
208 vn[iit->indexInInside()] = k_avg*(pl_inside-pl_outside)/distance;
215 Dune::FieldVector<DF,dim> d = iit->geometry().global(face_local);
216 d -= inside_cell_center_global;
217 RF distance = d.two_norm();
221 bctype = t.bctype(*iit,face_local);
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;
238 typename T::Traits::RangeFieldType j = t.j(*iit,face_local);
239 vn[iit->indexInInside()] = j;
244 Dune::FieldVector<DF,dim> vstar=iit->centerUnitOuterNormal();
245 vstar *= vn[iit->indexInInside()];
246 Dune::FieldVector<RF,dim> normalhat(0);
247 if (iit->indexInInside()%2==0)
248 normalhat[iit->indexInInside()/2] = -1.0;
250 normalhat[iit->indexInInside()/2] = 1.0;
251 Dune::FieldVector<DF,dim> vstarhat(0);
252 B.umtv(vstar,vstarhat);
253 vstarhat *= determinant;
254 storedcoeffs[index][iit->indexInInside()] = vstarhat*normalhat;
260 if (gv.grid().comm().size()>1)
261 gv.grid().communicate(dh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
265 void set_time (
typename T::Traits::RangeFieldType time_)
275 int index = is.index(e);
278 rt0fe.localBasis().evaluateFunction(x,rt0vectors);
280 for (
unsigned int i=0; i<rt0fe.localBasis().size(); i++)
281 yhat.axpy(storedcoeffs[index][i],rt0vectors[i]);
284 if (index != cachedindex)
286 B = e.geometry().jacobianTransposed(x);
287 determinant = B.determinant();
297 return pl.getGridView();
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
R RangeType
range type
Definition: function.hh:60
bool contains(int dim, int codim) const
returns true if data for this codim should be communicated
Definition: darcy_CCFV.hh:42