This source file includes following definitions.
- euclidean_norm
- evaluate
- evaluate
- jacobian
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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
46 using Vector_fef = Esfem::Grid::Vec_FEfun::Dune_FEfun;
47
48
49
50
51
52
53 template<typename T>
54 using Local_function = typename T::LocalFunctionType;
55
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
62
63 MCF_op(const Io::Parameter&,
64 const Grid::Grid_and_time&,
65 const Scalar_fef& u_input);
66
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
73
74
75
76
77
78
79
80
81 void brusselator_rhs(const Vector_fef& rhs, Vector_fef& lhs) const;
82
83
84
85
86
87
88
89
90
91
92
93
94
95
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
106 using Geometry
107 = Vector_fef::DiscreteFunctionSpaceType::IteratorType::Entity::Geometry;
108
109 using Grid_part = Vector_fef::GridPartType;
110
111 using Quadrature = Dune::Fem::CachingQuadrature<Grid_part, 0>;
112 template<typename T>
113 using Jacobian_range = typename Local_function<T>::JacobianRangeType;
114
115
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
126
127
128
129
130
131
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
139
140
141
142
143
144
145
146 Range<Vector_fef> surface_normal(const Geometry&) const;
147
148
149 const double alpha;
150
151 const double delta;
152 const double epsilon;
153
154 const double eps;
155 const Dune::Fem::TimeProviderBase& tp;
156 const Scalar_fef& u;
157
158 };
159
160
161 std::vector<MCF_op::Domain<MCF_op::Vector_fef> >
162 oriented_basis(const MCF_op::Geometry&);
163
164
165
166
167
168
169 MCF_op::Range<MCF_op::Vector_fef> nonUnit_normal
170 (const std::vector<MCF_op::Domain<MCF_op::Vector_fef> >& basis);
171
172
173
174
175
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
183
184
185
186
187
188
189
190
191
192
193
194
195
196
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
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
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
226
227
228 }
229 }
230
231 #endif