1 /*! \file secOrd_op_solutionDriven_impl.h
2 \author Christian Power
3 \date 18. March 2016
4
5 \brief Helper classes for `Esfem::SecOrd_op::Solution_driven`
6
7 Revision history
8 --------------------------------------------------
9
10 Revised by Christian Power March 2016
11 Originally written by Christian Power
12 (power22c@gmail.com) March 2016
13
14 Idea
15 --------------------------------------------------
16
17 Implementation of `secOrd_op_solutionDriven.h`
18
19 Created by Christian Power on 18.03.2016
20 Copyright (c) 2016 Christian Power. All rights reserved.
21
22 */
23
24 #ifndef SECORD_OP_SOLUTIONDRIVEN_IMPL_H
25 #define SECORD_OP_SOLUTIONDRIVEN_IMPL_H
26
27 #include <cmath>
28 #include <vector>
29 #include <config.h>
30 #include <dune/fem/operator/common/operator.hh>
31 #include <dune/fem/operator/linear/istloperator.hh>
32 #include <dune/fem/solver/istlsolver.hh>
33 #include <dune/fem/operator/common/differentiableoperator.hh>
34 #include <dune/fem/quadrature/cachingquadrature.hh>
35 #include "io_parameter.h"
36 #include "grid.h"
37
38 namespace Esfem{
39 namespace Impl{
40 class MCF_op
41 : public Dune::Fem::Operator<Esfem::Grid::Vec_FEfun::Dune_FEfun>
42 {
43 public:
44 using Scalar_fef = Esfem::Grid::Scal_FEfun::Dune_FEfun;
45 /*!< \brief \f$ f\colon \R^3 \to \R \f$ */
46 using Vector_fef = Esfem::Grid::Vec_FEfun::Dune_FEfun;
47 /*!< \brief \f$ f\colon \R^3 \to \R^3 \f$ */
48 // Big Matrix
49 // using Linear_operator
50 // = Dune::Fem::ISTLLinearOperator<Vector_fef, Vector_fef>;
51 // Sub matrix for an element
52 // using Local_matrix = Linear_operator::LocalMatrixType;
53 template<typename T>
54 using Local_function = typename T::LocalFunctionType;
55 /*!< \brief Finite element function restricted to a triangle. */
56 template<typename T>
57 using Domain = typename Local_function<T>::DomainType;
58 template<typename T>
59 using Range = typename Local_function<T>::RangeType;
60 using RangeField = typename Local_function<Vector_fef>::RangeFieldType;
61 /*!< \brief Simply \f$ \R \f$ */
62
63 MCF_op(const Io::Parameter&,
64 const Grid::Grid_and_time&,
65 const Scalar_fef& u_input);
66 /*!< \sa Constructor of `Esfem::SecOrd_op::Solution_driven`
67 */
68 MCF_op(const MCF_op&) = delete;
69 MCF_op& operator=(const MCF_op&) = delete;
70
71 void operator()(const Vector_fef& rhs, Vector_fef& lhs) const override;
72 /*!< \brief Applying the mean curvature operator.
73 The cg solver uses this.
74
75 The precise formulation reads as
76 \f{equation*}{
77 \parentheses[\big]{M_3^n + (\alpha + \varepsilon\tau) A_3^n}
78 \nodalValue{X}^{n+1}
79 \f}
80 */
81 void brusselator_rhs(const Vector_fef& rhs, Vector_fef& lhs) const;
82 /*!< \brief Generates rhs for the linear system.
83
84 The new value of the finite element function will be
85 \f{equation*}{
86 (M_3^n + \alpha A_3^n) \nodalValue{X}^n
87 + \tau \delta M_3^n(\nodalValue{u}^n,
88 \nodalValue{\surfaceNormal}).
89 \f}
90 */
91
92 // Assemble matrix for linear system
93 // void jacobian_matrix(const Vector_fef&, Linear_operator&) const;
94
95 /*! \name Variable containing dimensions
96 */
97 //@{
98 static constexpr int dim_vec_domain = Domain<Vector_fef>::dimension;
99 static constexpr int dim_vec_range = Range<Vector_fef>::dimension;
100 static constexpr int dim_scalar_domain = Domain<Scalar_fef>::dimension;
101 static constexpr int dim_scalar_range = Range<Scalar_fef>::dimension;
102 //@}
103
104 // ------------------------------------------------------------
105 // Typedef for helper functions
106 using Geometry
107 = Vector_fef::DiscreteFunctionSpaceType::IteratorType::Entity::Geometry;
108 /*!< \brief Geometry means access to the finite element. */
109 using Grid_part = Vector_fef::GridPartType;
110 /*!< \brief Auxiliary template argument */
111 using Quadrature = Dune::Fem::CachingQuadrature<Grid_part, 0>;
112 template<typename T>
113 using Jacobian_range = typename Local_function<T>::JacobianRangeType;
114 /*!< \brief Jacobian of the local finite element function
115 \warning `Jacobian_range` is used like this: `jac[0][1]`.
116 */
117
118 static_assert(dim_vec_domain == dim_vec_range, "Bad dimension");
119 static_assert(dim_scalar_range == 1, "Bad dimension");
120 private:
121 void mcf_lhs_matrixFree_assembly(const Geometry&,
122 const Quadrature&,
123 const Local_function<Vector_fef>&,
124 Local_function<Vector_fef>&) const;
125 /*!< \brief Used by operator()
126
127 The matrix vector formulation reads as
128 \f{equation*}{
129 \parentheses[\big]{M + (\alpha + \varepsilon \tau) A}
130 \nodalValue{X}
131 \f}
132 */
133 void mcf_rhs_matrixFree_assembly(const Geometry&,
134 const Quadrature&,
135 const Local_function<Vector_fef>&,
136 const Local_function<Scalar_fef>& u_loc,
137 Local_function<Vector_fef>&) const;
138 /*!< \brief Used by rhs()
139
140 The matrix vector formulation reads as
141 \f{equation*}{
142 (M + \alpha A) \nodalValue{X} + \tau \delta
143 M(\nodalValue{u}, \surfaceNormal)
144 \f}
145 */
146 Range<Vector_fef> surface_normal(const Geometry&) const;
147 /*!< \brief Calculates the outwards pointing normal vector. */
148
149 const double alpha;
150 /*!< \brief Velocity Laplace regularization parameter by Lubich */
151 const double delta; /*!< \brief Tumor growth parameter */
152 const double epsilon;
153 /*!< \brief Mean curvature regularization parameter by Elliott */
154 const double eps; /*!< \brief Generic tolerance */
155 const Dune::Fem::TimeProviderBase& tp; /*!< \brief Time step provider */
156 const Scalar_fef& u;
157 /*!< \brief Concentration of the growth promoting chemical substance */
158 };
159 /*!< \sa `Esfem::SecOrd_op::Solution_driven` */
160
161 std::vector<MCF_op::Domain<MCF_op::Vector_fef> >
162 oriented_basis(const MCF_op::Geometry&);
163 /*!< \brief Calculates an oriented basis for the tangent space.
164 Assumes `grid dimension == 2` and
165 `world dimension == 3`.
166 \return Two (numerical) vectors which represent an
167 oriented basis of the tangent space.
168 */
169 MCF_op::Range<MCF_op::Vector_fef> nonUnit_normal
170 (const std::vector<MCF_op::Domain<MCF_op::Vector_fef> >& basis);
171 /*!< \brief Calculates an non-normalized normal vector via
172 the cross product formula.
173 \warning Assumes `basis.size() == 2`
174 and that the basis is correctly oriented.
175 \return Outward pointing non unit normal vector.
176 */
177 inline MCF_op::RangeField
178 euclidean_norm(const MCF_op::Range<MCF_op::Vector_fef>& v){
179 static_assert(MCF_op::dim_vec_range == 3, "Bad dimension");
180 return std::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
181 }
182 /*!< \brief Euclidean norm for a 3-dimensional vector. */
183
184 /*! \name Evaluation helper functions
185 \param pt Quadrature point obtained by `MCF_op::Quadrature`.
186 \param q The quadrature object itself
187 \param cf The local finit element function
188
189 Evaluates the local finite element function
190 on the quadrature points. I have to pass the point and
191 the quadrature
192 object itself since I do not know what its precise type is.
193 I do not consider using template a good practise for this
194 problem.
195
196 \todo Figure out the type of `q[pt]`.
197 */
198 //@{
199 inline MCF_op::Range<MCF_op::Vector_fef>
200 evaluate(const std::size_t pt,
201 const MCF_op::Quadrature& q,
202 const MCF_op::Local_function<MCF_op::Vector_fef>& cf){
203 MCF_op::Range<MCF_op::Vector_fef> X;
204 cf.evaluate(q[pt], X);
205 return X;
206 }
207 /*! \brief Returns `X(p)` */
208 inline MCF_op::Range<MCF_op::Scalar_fef>
209 evaluate(const std::size_t pt,
210 const MCF_op::Quadrature& q,
211 const MCF_op::Local_function<MCF_op::Scalar_fef>& cf){
212 MCF_op::Range<MCF_op::Scalar_fef> u;
213 cf.evaluate(q[pt], u);
214 return u;
215 }
216 /*!< \brief Returns `u(p)` */
217 inline MCF_op::Jacobian_range<MCF_op::Vector_fef>
218 jacobian(const std::size_t pt,
219 const MCF_op::Quadrature& q,
220 const MCF_op::Local_function<MCF_op::Vector_fef>& cf){
221 MCF_op::Jacobian_range<MCF_op::Vector_fef> nabla_X;
222 cf.jacobian(q[pt], nabla_X);
223 return nabla_X;
224 }
225 /*!< \brief Returns `dX(p)`*/
226 //@}
227
228 } // namespace Impl
229 } // namespace Esfem
230
231 #endif // SECORD_OP_SOLUTIONDRIVEN_IMPL_H