DUNE PDELab (git)

darcyccfv.hh
1// -*- tab-width: 2; indent-tabs-mode: nil -*-
2#ifndef DUNE_PDELAB_LOCALOPERATOR_DARCYCCFV_HH
3#define DUNE_PDELAB_LOCALOPERATOR_DARCYCCFV_HH
4
8#include<dune/geometry/referenceelements.hh>
10
11#include<dune/pdelab/common/geometrywrapper.hh>
12#include<dune/pdelab/common/function.hh>
13#include<dune/pdelab/localoperator/pattern.hh>
14#include<dune/pdelab/localoperator/flags.hh>
15#include<dune/pdelab/localoperator/idefault.hh>
16#include<dune/pdelab/localoperator/defaultimp.hh>
17#include<dune/pdelab/localoperator/convectiondiffusionparameter.hh>
18
19// A DataHandle class to exchange RT0 coefficients in the overlap
20template<class GV, class V> // mapper type and vector type
21class VectorExchange
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
32public:
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
90template<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;
109 using BCType = typename Dune::PDELab::ConvectionDiffusionBoundaryConditions::Type;
110
111 const T& t;
112 const PL& pl;
114
116 mutable RF determinant;
117 mutable int cachedindex;
118 typename T::Traits::RangeFieldType time;
119
121 GV gv;
122 const IndexSet& is;
123 std::vector<RT0Coeffs> storedcoeffs;
124 mutable std::vector<RT0RangeType> rt0vectors;
125
126public:
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();
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
211 if (bctype==Dune::PDELab::ConvectionDiffusionBoundaryConditions::Dirichlet)
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();
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
223 if (bctype==Dune::PDELab::ConvectionDiffusionBoundaryConditions::Neumann)
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)
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
Provide velocity field for liquid phase.
Definition: darcyccfv.hh:97
CommDataHandleIF describes the features of a data handle for communication in parallel runs using the...
Definition: datahandleif.hh:78
void scatter(MessageBufferImp &buff, const EntityType &e, size_t n)
unpack data from message buffer to user.
Definition: datahandleif.hh:143
bool contains(int dim, int codim) const
returns true if data for given valid codim should be communicated
Definition: datahandleif.hh:94
size_t size(const EntityType &e) const
how many objects of type DataType have to be sent for a given entity
Definition: datahandleif.hh:118
void gather(MessageBufferImp &buff, const EntityType &e) const
pack data from user to message buffer
Definition: datahandleif.hh:129
bool fixedSize(int dim, int codim) const
returns true if size of data per entity of given dim and codim is a constant
Definition: datahandleif.hh:107
void umtv(const X &x, Y &y) const
y += A^T x
Definition: densematrix.hh:419
field_type determinant(bool doPivoting=true) const
calculates the determinant of this matrix
FieldTraits< value_type >::real_type two_norm() const
two norm sqrt(sum over squared values of entries)
Definition: densevector.hh:641
vector space out of a tensor product of fields.
Definition: fvector.hh:91
IndexTypeImp IndexType
The type used for the indices.
Definition: indexidset.hh:92
leaf of a function tree
Definition: function.hh:302
Raviart-Thomas local finite elements for cubes.
Definition: raviartthomascube.hh:40
A few common exception classes.
Traits for type conversions and type information.
Implements a vector constructed from a given type representing a field and a compile-time given size.
@ ForwardCommunication
communicate as given in InterfaceType
Definition: gridenums.hh:171
@ InteriorBorder_All_Interface
send interior and border, receive all entities
Definition: gridenums.hh:88
unspecified value type referenceElement(T &&... t)
Returns a reference element for the objects t....
constexpr Interior interior
PartitionSet for the interior partition.
Definition: partitionset.hh:271
Convenience header that includes all available Raviart-Thomas local finite elements for cubes.
traits class holding the function signature, same as in local function
Definition: function.hh:183
GV::Traits::template Codim< 0 >::Entity ElementType
codim 0 entity
Definition: function.hh:119
GV GridViewType
The type of the grid view the function lives on.
Definition: function.hh:116
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)