This source file includes following definitions.
- massMatrix_loc
- stiffnessMatrix_loc
- quad_massMatrix_loc
- jacobian_quadMass_loc
- vec_jacRange
- mass_matrix
- massMatrix_constOne
- jacobian
- heat_model
- quad_massMatrix_model
- jacobian_matrix_heat
- jacobian_matrix_quadMass
- calculate_matrix_row_size
- prepare_linearOperator_matrix
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16 #include <config.h>
17 #include <dune/fem/quadrature/cachingquadrature.hh>
18 #include "secOrd_op_brusselator_impl.h"
19 #include "io_parameter.h"
20 #include "grid.h"
21
22 #ifdef DEBUG
23 #include <iostream>
24 #endif
25
26 using namespace std;
27 using FE_function = Esfem::Grid::Scal_FEfun::Dune_FEfun;
28 using FE_space = FE_function::DiscreteFunctionSpaceType;
29 using Entity = FE_space::IteratorType::Entity;
30 using Geometry = Entity::Geometry;
31 using Grid_part = FE_function::GridPartType;
32 using Quadrature = Dune::Fem::CachingQuadrature<Grid_part, 0>;
33 using Local_function = FE_function::LocalFunctionType;
34 using Domain = Local_function::DomainType;
35 using Range = Local_function::RangeType;
36 using Jacobian_range = Local_function::JacobianRangeType;
37 using Range_field = Local_function::RangeFieldType;
38
39 static constexpr int dim_domain
40 = Esfem::Grid::Grid_and_time::Function_space::dimDomain;
41 static constexpr int dim_range
42 = Esfem::Grid::Grid_and_time::Function_space::dimRange;
43
44 using Linear_operator = Dune::Fem::ISTLLinearOperator<FE_function, FE_function>;
45 using Local_matrix = Linear_operator::LocalMatrixType;
46 using Diffusion_tensor = Dune::FieldMatrix<Range_field, dim_domain, dim_domain>;
47
48 inline Range massMatrix_loc(const std::size_t pt, const Quadrature& q,
49 const Local_function& cf){
50 Range u;
51 cf.evaluate(q[pt], u);
52 return u;
53 }
54 inline Jacobian_range stiffnessMatrix_loc(const std::size_t pt, const Quadrature& q,
55 const Local_function& cf){
56 Jacobian_range nabla_u;
57 cf.jacobian(q[pt], nabla_u);
58 return nabla_u;
59 }
60 inline Range quad_massMatrix_loc(const std::size_t pt, const Quadrature& q,
61 const Local_function& rhs_loc,
62 const Local_function& arg1_loc,
63 const Local_function& arg2_loc){
64 Range u;
65 Range tmp;
66 rhs_loc.evaluate(q[pt], u);
67 arg1_loc.evaluate(q[pt], tmp);
68 u *= tmp;
69 arg2_loc.evaluate(q[pt], tmp);
70 u *= tmp;
71 return u;
72 }
73 inline Range jacobian_quadMass_loc(const std::size_t pt, const Quadrature& q,
74 const Local_function& arg1_loc,
75 const Local_function& arg2_loc){
76 Range u;
77 Range tmp;
78 arg1_loc.evaluate(q[pt], u);
79 arg2_loc.evaluate(q[pt], tmp);
80 u *= tmp;
81 return u;
82 }
83
84 static void prepare_linearOperator_matrix(const FE_space&, Linear_operator&);
85 static std::size_t calculate_matrix_row_size(const FE_function&);
86
87
88
89
90 struct Brusselator_op::Data{
91 const FE_function& first_arg;
92 const FE_function& second_arg;
93 std::vector<Range> vec_range {};
94 std::vector<Jacobian_range> vec_jacRange {};
95 double massMatrix_lhs {0.};
96 double stiffnessMatrix_lhs {0.};
97 double quadMassMatrix_lhs {0.};
98 double massMatrix_rhs {0.};
99 Data(const Esfem::Io::Parameter&, const Esfem::Grid::Grid_and_time&,
100 const Esfem::Growth, const FE_function&, const FE_function&);
101 };
102
103 Brusselator_op::Data::Data(const Esfem::Io::Parameter& p,
104 const Esfem::Grid::Grid_and_time& gt,
105 const Esfem::Growth type,
106 const FE_function& quadMassMatrix_firstArg,
107 const FE_function& quadMassMatrix_secondArg)
108 :first_arg {quadMassMatrix_firstArg},
109 second_arg {quadMassMatrix_secondArg},
110 vec_range(calculate_matrix_row_size(first_arg)),
111 vec_jacRange(calculate_matrix_row_size(first_arg))
112 {
113 const auto dT = gt.time_provider().deltaT();
114 switch(type){
115 case Esfem::Growth::promoting:
116 massMatrix_lhs = 1 + dT * p.tg_gamma();
117 stiffnessMatrix_lhs = dT;
118 quadMassMatrix_lhs = (-1) * dT * p.tg_gamma();
119 massMatrix_rhs = dT * p.tg_gamma() * p.tg_a();
120 break;
121 case Esfem::Growth::inhibiting:
122 massMatrix_lhs = 1;
123 stiffnessMatrix_lhs = dT * p.tg_Dc();
124 quadMassMatrix_lhs = dT * p.tg_gamma();
125 massMatrix_rhs = dT * p.tg_gamma() * p.tg_b();
126 break;
127 default:
128 throw std::logic_error {"Error in constructor of Data. "
129 "Unkown growth type."};
130 break;
131 };
132 }
133
134 Brusselator_op::Brusselator_op(const Esfem::Io::Parameter& p,
135 const Esfem::Grid::Grid_and_time& gt,
136 const Esfem::Growth type,
137 const FE_function& quadMassMatrix_firstArg,
138 const FE_function& quadMassMatrix_secondArg)
139 try :d_ptr {make_unique<Data>(p, gt, type,
140 quadMassMatrix_firstArg, quadMassMatrix_secondArg)}
141 {}
142 catch(const std::exception&){
143 std::throw_with_nested(std::logic_error {"Error in constructor of "
144 "Brusselator_op."});
145 }
146 catch(...){
147 throw std::logic_error {"Unkown error in constructor of Brusselator_op."};
148 }
149 Brusselator_op::~Brusselator_op() = default;
150
151 void Brusselator_op::operator()(const FE_function& rhs, FE_function& lhs) const{
152 lhs.clear();
153 const auto& df_space = lhs.space();
154 for(const auto& entity : df_space){
155 const auto& rhs_loc = rhs.localFunction(entity);
156 auto lhs_loc = lhs.localFunction(entity);
157 heat_model(entity, rhs_loc, lhs_loc);
158 quad_massMatrix_model(entity, rhs_loc, lhs_loc);
159 }
160 }
161 void Brusselator_op::mass_matrix(const FE_function& rhs, FE_function& lhs) const{
162 lhs.clear();
163 const auto& df_space = lhs.space();
164 for(const auto& entity : df_space){
165 const auto& rhs_loc = rhs.localFunction(entity);
166 auto lhs_loc = lhs.localFunction(entity);
167 const auto& g = entity.geometry();
168 Quadrature q {entity, lhs_loc.order() + rhs_loc.order()};
169 for(std::size_t pt = 0; pt < q.nop(); ++pt){
170 const auto& x = q.point(pt);
171 const auto weight = q.weight(pt) * g.integrationElement(x);
172 auto Mu = massMatrix_loc(pt, q, rhs_loc);
173
174 Mu *= weight;
175 lhs_loc.axpy(q[pt], Mu);
176 }
177 }
178 }
179 void Brusselator_op::massMatrix_constOne(FE_function& lhs) const{
180 lhs.clear();
181 const auto& df_space = lhs.space();
182 for(const auto& entity : df_space){
183 auto lhs_loc = lhs.localFunction(entity);
184 Quadrature q {entity, lhs_loc.order()};
185 const auto& g = entity.geometry();
186 for(std::size_t pt = 0; pt < q.nop(); ++pt){
187 const auto& x = q.point(pt);
188 const auto weight = q.weight(pt) * g.integrationElement(x);
189
190 Range Mu {d_ptr -> massMatrix_rhs * weight};
191 lhs_loc.axpy(q[pt], Mu);
192 }
193 }
194 }
195 void Brusselator_op::jacobian(const FE_function& fef, Linear_operator& matrix) const{
196 const auto& df_space = fef.space();
197
198 prepare_linearOperator_matrix(df_space, matrix);
199
200 for(const auto& entity : df_space){
201 const auto fef_loc = fef.localFunction(entity);
202 auto matrix_loc = matrix.localMatrix(entity, entity);
203
204 jacobian_matrix_heat(entity, fef_loc, matrix_loc);
205 jacobian_matrix_quadMass(entity, fef_loc, matrix_loc);
206 }
207 matrix.communicate();
208 }
209
210
211
212
213 void Brusselator_op::heat_model(const Entity& e, const Local_function& rhs_loc,
214 Local_function& lhs_loc) const{
215 const auto& g = e.geometry();
216 Quadrature q {e, lhs_loc.order() + rhs_loc.order()};
217
218 for(std::size_t pt = 0; pt < q.nop(); ++pt){
219
220 const auto& x = q.point(pt);
221 const auto weight = q.weight(pt) * g.integrationElement(x);
222
223 auto Mu = massMatrix_loc(pt, q, rhs_loc);
224 auto Au = stiffnessMatrix_loc(pt, q, rhs_loc);
225
226 Mu *= d_ptr -> massMatrix_lhs * weight;
227 Au *= d_ptr -> stiffnessMatrix_lhs * weight;
228
229 lhs_loc.axpy(q[pt], Mu, Au);
230 }
231 }
232 void Brusselator_op::
233 quad_massMatrix_model(const Entity& e, const Local_function& rhs_loc,
234 Local_function& lhs_loc) const{
235 const auto& g = e.geometry();
236 const auto& arg1_loc = d_ptr -> first_arg.localFunction(e);
237 const auto& arg2_loc = d_ptr -> second_arg.localFunction(e);
238 Quadrature q
239 {e, lhs_loc.order() + arg1_loc.order() + arg2_loc.order() + rhs_loc.order()};
240
241 for(std::size_t pt = 0; pt < q.nop(); ++pt){
242 const auto& x = q.point(pt);
243 const auto integral_factor
244 = q.weight(pt) * g.integrationElement(x)
245 * d_ptr -> quadMassMatrix_lhs;
246
247 auto Mu = quad_massMatrix_loc(pt, q, rhs_loc, arg1_loc, arg2_loc);
248 Mu *= integral_factor;
249
250 lhs_loc.axpy(q[pt], Mu);
251 }
252 }
253 void Brusselator_op::
254 jacobian_matrix_heat(const Entity& entity,
255 const Local_function& fef_loc, Local_matrix& matrix_loc) const{
256 auto& vec_range = d_ptr -> vec_range;
257 auto& vec_jacRange = d_ptr -> vec_jacRange;
258
259 const auto& geometry = entity.geometry();
260 const auto& base_set = matrix_loc.domainBasisFunctionSet();
261
262 Quadrature q {entity, 2 * fef_loc.order()};
263 for(std::size_t pt = 0; pt < q.nop(); ++pt){
264 const auto& x = q.point(pt);
265 const auto integral_factor = q.weight(pt) * geometry.integrationElement(x);
266
267 base_set.evaluateAll(q[pt], vec_range);
268 base_set.jacobianAll(q[pt], vec_jacRange);
269
270
271
272
273 for(std::size_t localCol = 0; localCol < base_set.size(); ++localCol){
274 Range basisFunction = vec_range[localCol];
275 basisFunction *= d_ptr -> massMatrix_lhs;
276 Jacobian_range grad_basisFunction = vec_jacRange[localCol];
277 grad_basisFunction *= d_ptr -> stiffnessMatrix_lhs;
278 matrix_loc
279 .column(localCol)
280 .axpy(vec_range, vec_jacRange,
281 basisFunction, grad_basisFunction,
282 integral_factor);
283 }
284 }
285 }
286 void Brusselator_op::
287 jacobian_matrix_quadMass(const Entity& entity,
288 const Local_function& fef_loc, Local_matrix& matrix_loc) const{
289 auto& vec_range = d_ptr -> vec_range;
290 const auto& arg1_loc = d_ptr -> first_arg.localFunction(entity);
291 const auto& arg2_loc = d_ptr -> second_arg.localFunction(entity);
292
293 const auto& geometry = entity.geometry();
294 const auto& base_set = matrix_loc.domainBasisFunctionSet();
295
296 Quadrature q {entity, 2 * fef_loc.order() + arg1_loc.order() + arg2_loc.order()};
297 for(std::size_t pt = 0; pt < q.nop(); ++pt){
298 const auto& x = q.point(pt);
299 const auto integral_factor
300 = q.weight(pt) * geometry.integrationElement(x)
301 * d_ptr -> quadMassMatrix_lhs;
302
303 base_set.evaluateAll(q[pt], vec_range);
304 const auto u0 = jacobian_quadMass_loc(pt, q, arg1_loc, arg2_loc);
305
306 for(std::size_t localCol = 0; localCol < base_set.size(); ++localCol){
307 Range basisFunction = vec_range[localCol] * u0;
308 matrix_loc
309 .column(localCol)
310 .axpy(vec_range, basisFunction, integral_factor);
311 }
312 }
313 }
314
315
316
317
318 std::size_t calculate_matrix_row_size(const FE_function& fef){
319 const auto& df_space = fef.space();
320 return df_space.blockMapper().maxNumDofs() * df_space.localBlockSize;
321
322 }
323 void prepare_linearOperator_matrix(const FE_space& df_space,
324 Linear_operator& matrix){
325 Dune::Fem::DiagonalStencil<FE_space, FE_space> stencil {df_space, df_space};
326 matrix.reserve(stencil);
327 matrix.clear();
328 }
329
330
331
332