1 /*! \file secOrd_op_rhs_impl.h
2 \brief Helper classes for `secOrd_op_rhs.cpp`
3
4 Revision history
5 --------------------------------------------------
6
7 Revised by Christian Power April 2016
8 Originally written by Christian Power
9 (power22c@gmail.com) April 2016
10
11 Idea
12 --------------------------------------------------
13
14 Implement `Vec_rhs::Data`
15
16 \author Christian Power
17 \date 16. April 2016
18 \copyright Copyright (c) 2016 Christian Power. All rights reserved.
19 */
20
21 #ifndef SECORD_OP_RHS_IMPL_H
22 #define SECORD_OP_RHS_IMPL_H
23
24 #include <functional>
25 #include <config.h>
26 #include <dune/fem/quadrature/cachingquadrature.hh>
27 #include "grid.h"
28 #include "secOrd_op_rhs.h"
29
30
31 namespace Esfem{
32 //! Implementation details
33 namespace Impl{
34 //! Scalar valued right-hand side for the surface PDE
35 class Rhs_fun
36 :public Dune::Fem::Function
37 <Esfem::Grid::Grid_and_time::Function_space, Rhs_fun>
38 {
39 public:
40 //! Template argument
41 using Base = Esfem::Grid::Grid_and_time::Function_space;
42 //! \f$ \R^3\f$
43 using Domain = Base::DomainType;
44 //! \f$ \R^1\f$
45 using Range = Base::RangeType;
46
47 static_assert(Domain::dimension == 3, "Bad Domain dimension");
48 static_assert(Range::dimension == 1, "Bad Range dimension");
49
50 //! Get time, time step and indicate right-hand side for `u` or `w`
51 /*! \warning If time provider outlives `this`, then you are in problem. */
52 Rhs_fun(const Dune::Fem::TimeProviderBase&,
53 const Growth);
54 //! No copy constructor
55 Rhs_fun(const Rhs_fun&) = delete;
56 //! No copy assignment
57 Rhs_fun& operator=(const Rhs_fun&) = delete;
58
59 //! y = f(x)
60 void evaluate(const Domain&, Range&) const;
61 //! \copybrief evaluate()
62 Range operator()(const Domain&) const;
63 private:
64 //! Current time and time step
65 const Dune::Fem::TimeProviderBase& tp;
66 //! Right-hand side for \f$ u\f$ or \f$ w\f$
67 std::function<void(const Domain&,Range&)> fun_impl;
68 //! Dynamic assert
69 void dassert(const bool assertion, const std::string& msg);
70 };
71
72 //! Vector valued right-hand side for the mean curvature equation
73 /*! This is my right-hand function:
74 \f[
75 \biggl[ k \Bigl( 1 - \frac{r(t)}{r_{end}}\Bigr) \lvert X \rvert
76 + 2 \Bigl(\alpha k \Bigl( 1 - \frac{r(t)}{r_{end}}\Bigr)
77 + \varepsilon \Bigr) \frac{1}{\lvert X \rvert}
78 - \delta xy e^{-6t}\biggr]
79 \f]
80 Currently I am testing. So, I use an easier right-hand side.
81 I am doing classic mean curvature flow on the sphere.
82 Right-hand side must be zero for this example. Also I expect
83 unit sphere as starting value, alpha = 0 and epsilon = 1.
84
85 Old test, which I want to check again later:
86 - \f$ \surface_0 = S^2\f$
87 - Exact flow: \f$ \Phi(x,t) = e^{-t} x\f$
88 - Velocity: \f$ v(x,t) = -x\f$
89 - Normal: \f$ n = x/ |x|\f$
90 - Mean curvature: \f$ H = 2/ |x|\f$
91 - \f$\delta = \varepsilon = 0\f$, hence the operator is
92 \f[
93 v - \alpha \Delta v = g(x,t) = -(1 + \alpha 2/ |x|^2) x.
94 \f]
95 \todo Finish test and clean up comments.
96 */
97 class Vec_rhs_fun
98 :public Dune::Fem::Function
99 <Esfem::Grid::Grid_and_time::Vec_Function_space, Vec_rhs_fun>
100 {
101 public:
102 //! Template argument
103 using Base = Esfem::Grid::Grid_and_time::Vec_Function_space;
104 //! \f$ \R^3 \f$
105 using Domain = Base::DomainType;
106 //! \f$ \R^3 \f$
107 using Range = Base::RangeType;
108
109 static_assert(Domain::dimension == 3, "Bad Domain dimension");
110 static_assert(Range::dimension == 3, "Bad Range dimension");
111
112 //! Get time and time step
113 /*! \copydetails Rhs_fun::Rhs_fun() */
114 explicit Vec_rhs_fun(const Dune::Fem::TimeProviderBase&);
115 //! No copy construct
116 Vec_rhs_fun(const Vec_rhs_fun&) = delete;
117 //! No copy assignment
118 Vec_rhs_fun& operator=(const Vec_rhs_fun&) = delete;
119
120 //! \copybrief Rhs_fun::evaluate()
121 void evaluate(const Domain&, Range&) const;
122 //! \copybrief evaluate()
123 Range operator()(const Domain&) const;
124 private:
125 //! Time and time step
126 const Dune::Fem::TimeProviderBase& tp;
127 //! \f$\alpha \Delta v\f$
128 double alpha;
129 //! \f$ \varepsilon \Delta X\f$
130 double epsilon;
131 //! Initial size of sphere
132 double r_start;
133 //! Final size of sphere
134 double r_end;
135 //! Steepness of logistic growh
136 double k;
137 //! \f$\delta u\f$
138 double delta;
139 //! Helper variable, which save computations
140 mutable double cache[4];
141 //! Conditional update member `cache`
142 /*! `cache[0]` holds the last time. If the current time does not equals `cache[0]`,
143 then update the cache. The content of cache is
144 - `cache[0] = tp.time()`,
145 - \f$ k( 1 - \frac{r(t)}{r_{end}})\f$,
146 - \f$2 ( \alpha k( 1 - \frac{r(t)}{r_{end}}) + \varepsilon)\f$,
147 - \f$\delta e^{-6t}\f$,
148
149 where \f$r(t) = \frac{r_{end}r_0}{r_{end} e^{-kt} + r_0 (1 - e^{-kt})}\f$
150 */
151 void update_cache() const;
152 };
153
154 //! Right-hand side for surface logistic sphere experiment
155 /*! \pre Use evaluate() on the exact surface.
156 I assume that it is a sphere, such that \f$ r(t) = |x|\f$.
157 \sa Esfem::Brusselator_scheme::eoc_sls() */
158 struct sls_rhs
159 : Dune::Fem::Function
160 <Esfem::Grid::Grid_and_time::Vec_Function_space, sls_rhs>,
161 SecOrd_op::vRhs{
162 //! Dune function
163 using dBase = Esfem::Grid::Grid_and_time::Vec_Function_space;
164 //! \f$ \R^3\f$
165 using dom = dBase::DomainType;
166 //! \f$ \R^3\f$
167 using ran = dBase::RangeType;
168 //! For my generic algorithm
169 using Range = ran;
170
171 //! Get time and finit element space
172 /*! \post Grid and time outlive this object. */
173 sls_rhs(const Grid::Grid_and_time&);
174 sls_rhs* clone() override{ return new sls_rhs {*this}; }
175 void addScaled_to(Grid::Vec_FEfun& rhs) override;
176 //! Needed for interpolation
177 ran operator()(const dom&) const;
178 private:
179 //! Manifold dimension
180 static constexpr int dim {2};
181 //! Time step
182 const Dune::Fem::TimeProviderBase& tp;
183 //! Load vector
184 Esfem::Grid::Vec_FEfun::Dune_FEfun lvec;
185 //! Initial radius (population)
186 double rA;
187 //! Carrying capacity
188 double rE;
189 //! \f$\alpha\f$
190 double a;
191 //! \f$\varepsilon\f$
192 double e;
193 //! Growth rate
194 double k;
195 //! Growth rate scalar function
196 double delta;
197 };
198
199 //! Right-hand side for surface Dalquist test equation
200 /*! \pre Grid_and_time must outlive this object.
201 \sa Esfem::Brusselator_scheme::sd() */
202 struct sd_rhs
203 : Dune::Fem::Function
204 <Esfem::Grid::Grid_and_time::Vec_Function_space, sd_rhs>,
205 SecOrd_op::vRhs{
206 //! Dune function
207 using dBase = Esfem::Grid::Grid_and_time::Vec_Function_space;
208 //! \f$ \R^3\f$
209 using dom = dBase::DomainType;
210 //! \f$ \R^3\f$
211 using ran = dBase::RangeType;
212 //! For my generic algorithm
213 using Range = ran;
214
215 //! Get time and finit element space
216 /*! \post Grid and time outlive this object. */
217 sd_rhs(const Grid::Grid_and_time&);
218 sd_rhs* clone() override{ return new sd_rhs {*this}; }
219 void addScaled_to(Grid::Vec_FEfun& rhs) override;
220 //! Needed for interpolation
221 ran operator()(const dom&) const;
222 private:
223 //! Manifold dimension
224 static constexpr int dim {2};
225 //! Time step
226 const Dune::Fem::TimeProviderBase& tp;
227 //! Load vector
228 Esfem::Grid::Vec_FEfun::Dune_FEfun lvec;
229 //! \f$\alpha\f$
230 double a;
231 //! \f$\varepsilon\f$
232 double e;
233 };
234
235 //! Right-hand side for surface logistic sphere experiment
236 /*! \pre Use evaluate() on the exact surface.
237 I assume that it is a sphere, such that \f$ r(t) = |x|\f$.
238 \sa Esfem::Brusselator_scheme::eoc_sls() */
239 struct sdp_u_rhs
240 : Dune::Fem::Function
241 <Esfem::Grid::Grid_and_time::Function_space, sdp_u_rhs>,
242 SecOrd_op::sRhs{
243 //! Dune function
244 using dBase = Esfem::Grid::Grid_and_time::Function_space;
245 //! \f$ \R^3\f$
246 using dom = dBase::DomainType;
247 //! \f$ \R\f$
248 using ran = dBase::RangeType;
249 //! For my generic algorithm
250 using Range = ran;
251
252 //! Get time and finit element space
253 /*! \post Grid and time outlive this object. */
254 sdp_u_rhs(const Grid::Grid_and_time&);
255 sdp_u_rhs* clone() override{ return new sdp_u_rhs {*this}; }
256 void addScaled_to(Grid::Scal_FEfun& rhs) override;
257 //! Needed for interpolation
258 ran operator()(const dom&) const;
259 private:
260 //! Manifold dimension
261 static constexpr int dim {2};
262 //! Time step
263 const Dune::Fem::TimeProviderBase& tp;
264 //! Load vector
265 Esfem::Grid::Scal_FEfun::Dune_FEfun lscal;
266 //! Initial radius (population)
267 double rA;
268 //! Carrying capacity
269 double rE;
270 //! Growth rate
271 double k;
272 };
273
274 //! Assemble load vector
275 /*! \tparam Rhs Deduce Rhs_fun or Vec_rhs_fun
276 \tparam Fef Deduce
277 Scal_FEfun::Dune_FEfun or Vec_FEfun::Dune_FEfun
278 */
279 template<typename Rhs, typename Fef>
280 void assemble_RHS(const Rhs& rhs, Fef& fef);
281 } // namespace Impl
282
283 //! %Data details of Rhs
284 struct SecOrd_op::Rhs::Data{
285 //! Time and time step
286 const Dune::Fem::TimeProviderBase& tp;
287 //! Scalar valued right-hand side function
288 Impl::Rhs_fun rhs;
289 //! Finite element load vector
290 Esfem::Grid::Scal_FEfun::Dune_FEfun load_vector;
291
292 //! Get space-time manifold and finite element space
293 Data(const Grid::Grid_and_time&, const Growth);
294 };
295
296 //! %Data details of Vec_rhs
297 struct SecOrd_op::Vec_rhs::Data{
298 //! \copybrief Rhs::Data::tp
299 const Dune::Fem::TimeProviderBase& tp;
300 //! Vector valued right-hand side function
301 Impl::Vec_rhs_fun rhs;
302 //! \copybrief Rhs::Data::load_vector
303 Esfem::Grid::Vec_FEfun::Dune_FEfun load_vector;
304
305 //! \copybrief Rhs::Data::Data()
306 Data(const Grid::Grid_and_time&);
307 };
308 } // namespace Esfem
309
310 // ----------------------------------------------------------------------
311 // Template implementation
312
313 template<typename Rhs, typename Fef>
314 void Esfem::Impl::assemble_RHS(const Rhs& rhs, Fef& fef){
315 using Range = typename Rhs::Range;
316 using Grid_part = typename Fef::GridPartType;
317 using Quadrature = Dune::Fem::CachingQuadrature<Grid_part, 0>;
318
319 fef.clear();
320 const auto& df_space = fef.space();
321 for(const auto& entity : df_space){
322 const auto& geometry = entity.geometry();
323 const Quadrature quad {entity, 2 * df_space.order() + 1};
324 auto fef_local = fef.localFunction(entity);
325 for(std::size_t pt = 0; pt < quad.nop(); ++pt){
326 const auto& x = quad.point(pt);
327 Range fx {rhs(geometry.global(x))};
328 // rhs.evaluate(geometry.global(x), fx);
329 fx *= quad.weight(pt) * geometry.integrationElement(x);
330 fef_local.axpy(quad[pt], fx);
331 }
332 }
333 fef.communicate();
334 }
335
336 // ----------------------------------------------------------------------
337 // Inline implementation
338
339 inline void Esfem::Impl::Rhs_fun::evaluate(const Domain& d, Range& r) const{
340 fun_impl(d,r);
341 }
342 inline Esfem::Impl::Rhs_fun::Range
343 Esfem::Impl::Rhs_fun::operator()(const Domain& d) const{
344 Range r {0};
345 evaluate(d,r);
346 return r;
347 }
348
349 #endif // SECORD_OP_RHS_IMPL_H