2 #ifndef DUNE_PDELAB_LOCALOPERATOR_DARCYCCFV_HH 3 #define DUNE_PDELAB_LOCALOPERATOR_DARCYCCFV_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> 21 template<
class GV,
class V>
23 :
public Dune::CommDataHandleIF<VectorExchange<GV,V>,
24 typename V::value_type>
26 using IndexSet =
typename GV::IndexSet;
27 using IndexType =
typename IndexSet::IndexType;
31 const IndexSet& indexSet;
39 : gv(gv_), c(c_), indexSet(gv.indexSet())
58 template<
class EntityType>
59 size_t size (EntityType&
e)
const 65 template<
class MessageBuffer,
class EntityType>
66 void gather (MessageBuffer& buff,
const EntityType&
e)
const 68 buff.write(c[indexSet.index(e)]);
75 template<
class MessageBuffer,
class EntityType>
76 void scatter (MessageBuffer& buff,
const EntityType&
e,
size_t n)
80 c[indexSet.index(e)]=x;
91 template<
typename T,
typename PL>
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> >
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;
114 Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0> rt0fe;
116 mutable Dune::FieldMatrix<DF,dim,dim> B;
117 mutable RF determinant;
118 mutable int cachedindex;
119 typename T::Traits::RangeFieldType time;
121 using RT0Coeffs = Dune::FieldVector<RF,2*dim>;
124 std::vector<RT0Coeffs> storedcoeffs;
125 mutable std::vector<RT0RangeType> rt0vectors;
132 : t(t_), pl(pl_), cachedindex(-1), time(0), gv(pl_.getGridView()), is(gv.indexSet()), storedcoeffs(is.
size(0)),
133 rt0vectors(rt0fe.localBasis().
size())
136 for (
const auto& element : elements(gv,Dune::Partitions::interior))
139 int index = is.index(element);
142 auto geo = element.geometry();
146 auto inside_cell_center_local = ref_el.position(0,0);
147 auto inside_cell_center_global = geo.global(inside_cell_center_local);
150 auto tensor_inside = t.A(element,inside_cell_center_local);
153 typename PL::Traits::RangeType pl_inside;
154 pl.evaluate(element,inside_cell_center_local,pl_inside);
158 auto B = geo.jacobianInverseTransposed(inside_cell_center_local);
159 auto determinant = B.determinant();
162 for (
const auto& intersection : intersections(pl.getGridView(),element))
165 vn[intersection.indexInInside()] = 0.0;
171 if (intersection.neighbor())
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);
178 auto d(outside_cell_center_global);
179 d -= inside_cell_center_global;
180 auto distance = d.two_norm();
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));
193 typename PL::Traits::RangeType pl_outside;
194 pl.evaluate(outside_cell,outside_cell_center_local,pl_outside);
197 vn[intersection.indexInInside()] = k_avg*(pl_inside-pl_outside)/distance;
201 if (intersection.boundary())
204 auto d = intersection.geometry().global(face_local);
205 d -= inside_cell_center_global;
206 auto distance = d.two_norm();
209 auto bctype = t.bctype(intersection,face_local);
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;
226 auto j = t.j(intersection,face_local);
227 vn[intersection.indexInInside()] = j;
232 auto vstar=intersection.centerUnitOuterNormal();
233 vstar *= vn[intersection.indexInInside()];
234 Dune::FieldVector<RF,dim> normalhat(0);
235 if (intersection.indexInInside()%2==0)
236 normalhat[intersection.indexInInside()/2] = -1.0;
238 normalhat[intersection.indexInInside()/2] = 1.0;
239 Dune::FieldVector<DF,dim> vstarhat(0);
240 B.umtv(vstar,vstarhat);
241 vstarhat *= determinant;
242 storedcoeffs[index][intersection.indexInInside()] = vstarhat*normalhat;
248 if (gv.grid().comm().size()>1)
249 gv.grid().communicate(dh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
253 void set_time (
typename T::Traits::RangeFieldType time_)
263 int index = is.index(e);
266 rt0fe.localBasis().evaluateFunction(x,rt0vectors);
268 for (
unsigned int i=0; i<rt0fe.localBasis().size(); i++)
269 yhat.axpy(storedcoeffs[index][i],rt0vectors[i]);
272 if (index != cachedindex)
274 B = e.geometry().jacobianTransposed(x);
275 determinant = B.determinant();
285 return pl.getGridView();
289 #endif // DUNE_PDELAB_LOCALOPERATOR_DARCYCCFV_HH void evaluate(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeType &y) const
Definition: darcyccfv.hh:258
VectorExchange(const GV &gv_, V &c_)
constructor
Definition: darcyccfv.hh:38
static const int dim
Definition: adaptivity.hh:83
bool contains(int dim, int codim) const
returns true if data for this codim should be communicated
Definition: darcyccfv.hh:43
Type
Definition: convectiondiffusionparameter.hh:65
void gather(MessageBuffer &buff, const EntityType &e) const
pack data from user to message buffer
Definition: darcyccfv.hh:66
Dune::FieldVector< GV::Grid::ctype, GV::dimension > DomainType
domain type in dim-size coordinates
Definition: function.hh:48
typename V::value_type DataType
export type of data for message buffer
Definition: darcyccfv.hh:35
ReferenceElementWrapper< ReferenceElement< typename Geometry::ctype, Geometry::mydimension > > referenceElement(const Geometry &geo)
Returns the reference element for the given geometry.
Definition: referenceelements.hh:144
const Traits::GridViewType & getGridView() const
Definition: darcyccfv.hh:283
void scatter(MessageBuffer &buff, const EntityType &e, size_t n)
Definition: darcyccfv.hh:76
traits class holding the function signature, same as in local function
Definition: function.hh:175
Definition: convectiondiffusionparameter.hh:65
Provide velocity field for liquid phase.
Definition: darcyccfv.hh:92
const Entity & e
Definition: localfunctionspace.hh:111
Definition: convectiondiffusionparameter.hh:65
Definition: darcyccfv.hh:22
DarcyVelocityFromHeadCCFV(const T &t_, const PL &pl_)
Definition: darcyccfv.hh:131
GV::Traits::template Codim< 0 >::Entity ElementType
codim 0 entity
Definition: function.hh:117
size_t size(EntityType &e) const
Definition: darcyccfv.hh:59
GV GridViewType
The type of the grid view the function lives on.
Definition: function.hh:114
leaf of a function tree
Definition: function.hh:575
void set_time(typename T::Traits::RangeFieldType time_)
Definition: darcyccfv.hh:253
R RangeType
range type
Definition: function.hh:60
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