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>
20 template<
class GV,
class V>
22 :
public Dune::CommDataHandleIF<VectorExchange<GV,V>,
23 typename V::value_type>
25 using IndexSet =
typename GV::IndexSet;
26 using IndexType =
typename IndexSet::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 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;
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 using RT0Coeffs = Dune::FieldVector<RF,2*dim>;
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 for (
const auto& element : elements(gv,Dune::Partitions::interior))
138 int index = is.index(element);
141 auto geo = element.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);
149 auto tensor_inside = t.A(element,inside_cell_center_local);
152 typename PL::Traits::RangeType pl_inside;
153 pl.evaluate(element,inside_cell_center_local,pl_inside);
157 auto B = geo.jacobianInverseTransposed(inside_cell_center_local);
158 auto determinant = B.determinant();
161 for (
const auto& intersection : intersections(pl.getGridView(),element))
164 vn[intersection.indexInInside()] = 0.0;
167 auto face_local = referenceElement(intersection.geometry()).position(0,0);
170 if (intersection.neighbor())
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);
177 auto d(outside_cell_center_global);
178 d -= inside_cell_center_global;
179 auto distance = d.two_norm();
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));
192 typename PL::Traits::RangeType pl_outside;
193 pl.evaluate(outside_cell,outside_cell_center_local,pl_outside);
196 vn[intersection.indexInInside()] = k_avg*(pl_inside-pl_outside)/distance;
200 if (intersection.boundary())
203 auto d = intersection.geometry().global(face_local);
204 d -= inside_cell_center_global;
205 auto distance = d.two_norm();
208 auto bctype = t.bctype(intersection,face_local);
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;
225 auto j = t.j(intersection,face_local);
226 vn[intersection.indexInInside()] = j;
231 auto vstar=intersection.centerUnitOuterNormal();
232 vstar *= vn[intersection.indexInInside()];
233 Dune::FieldVector<RF,dim> normalhat(0);
234 if (intersection.indexInInside()%2==0)
235 normalhat[intersection.indexInInside()/2] = -1.0;
237 normalhat[intersection.indexInInside()/2] = 1.0;
238 Dune::FieldVector<DF,dim> vstarhat(0);
239 B.umtv(vstar,vstarhat);
240 vstarhat *= determinant;
241 storedcoeffs[
index][intersection.indexInInside()] = vstarhat*normalhat;
247 if (gv.grid().comm().size()>1)
248 gv.grid().communicate(dh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
252 void set_time (
typename T::Traits::RangeFieldType time_)
265 rt0fe.localBasis().evaluateFunction(x,rt0vectors);
267 for (
unsigned int i=0; i<rt0fe.localBasis().size(); i++)
268 yhat.axpy(storedcoeffs[
index][i],rt0vectors[i]);
271 if (
index != cachedindex)
273 B =
e.geometry().jacobianTransposed(x);
274 determinant = B.determinant();
284 return pl.getGridView();
288 #endif // DUNE_PDELAB_LOCALOPERATOR_DARCYCCFV_HH