This source file includes following definitions.
- mass_matrix
- stiffness_matrix
- solve
- apply_massMatrix_to
- mass_matrix
- mass_matrix
- mass_matrix
- matrixFree_assembly
- massMatrixFree_assembly
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16 #include <stdexcept>
17 #include <config.h>
18 #include <dune/fem/operator/common/operator.hh>
19 #include <dune/fem/solver/cginverseoperator.hh>
20 #include <dune/fem/solver/oemsolver.hh>
21 #include <dune/fem/quadrature/cachingquadrature.hh>
22 #include "secOrd_op_linearHeat.h"
23 #include "io_parameter.h"
24 #include "grid.h"
25
26 using namespace std;
27 using FE_function = Esfem::Grid::Scal_FEfun::Dune_FEfun;
28 using Solver = Dune::Fem::CGInverseOperator<FE_function>;
29 using Geometry
30 = FE_function::DiscreteFunctionSpaceType::IteratorType::Entity::Geometry;
31 using Grid_part = FE_function::GridPartType;
32 using Quadrature = Dune::Fem::CachingQuadrature<Grid_part, 0>;
33 using Domain = FE_function::LocalFunctionType::DomainType;
34 using Range = FE_function::LocalFunctionType::RangeType;
35 using Jacobian_range = FE_function::LocalFunctionType::JacobianRangeType;
36
37 class Linear_heat_op : public Dune::Fem::Operator<FE_function>{
38 public:
39 explicit Linear_heat_op(const Esfem::Io::Parameter&,
40 const Esfem::Grid::Grid_and_time&);
41 Linear_heat_op(const Linear_heat_op&) = delete;
42 Linear_heat_op& operator=(const Linear_heat_op&) = delete;
43
44 void operator()(const FE_function& rhs, FE_function& lhs) const override;
45 void mass_matrix(FE_function&);
46 void mass_matrix(const FE_function& rhs, FE_function& lhs) const;
47 private:
48 double bdf_alpha_lead_coeff {1.};
49 const Dune::Fem::TimeProviderBase& tp;
50 FE_function tmp_fef;
51 };
52
53 void matrixFree_assembly(const double dT, const Geometry&, const Quadrature&,
54 const FE_function::LocalFunctionType&,
55 FE_function::LocalFunctionType&);
56 void massMatrixFree_assembly(const Geometry&, const Quadrature&,
57 const FE_function::LocalFunctionType&,
58 FE_function::LocalFunctionType&);
59 inline
60 Range mass_matrix(const std::size_t pt, const Quadrature& q,
61 const FE_function::LocalFunctionType& cf){
62 Range u;
63 cf.evaluate(q[pt], u);
64 return u;
65 }
66 inline
67 Jacobian_range stiffness_matrix(const std::size_t pt, const Quadrature& q,
68 const FE_function::LocalFunctionType& cf){
69 Jacobian_range nabla_u;
70 cf.jacobian(q[pt], nabla_u);
71 return nabla_u;
72 }
73
74
75
76
77 struct Esfem::SecOrd_op::Linear_heat::Data{
78 Linear_heat_op heat_op;
79 Solver heat_solver;
80 Data(const Io::Parameter& p, const Grid::Grid_and_time& gt)
81 : heat_op {p, gt}, heat_solver {heat_op, p.eps(), p.eps()}
82 {}
83 };
84 Esfem::SecOrd_op::Linear_heat::Linear_heat(const Io::Parameter& p,
85 const Grid::Grid_and_time& gt)
86 try : d_ptr {make_unique<Data>(p, gt)}
87 {}
88 catch(const std::exception&){
89 std::throw_with_nested(std::logic_error{"Error in constructor of Linear_heat."});
90 }
91 catch(...){
92 throw std::logic_error{"Unkown error in constructor of Linear_heat."};
93 }
94
95 Esfem::SecOrd_op::Linear_heat::~Linear_heat() = default;
96
97 void Esfem::SecOrd_op::Linear_heat::
98 solve(const Grid::Scal_FEfun& rhs, Grid::Scal_FEfun& lhs) const{
99 const FE_function& fef1 = rhs;
100 FE_function& fef2 = lhs;
101 d_ptr -> heat_solver(fef1, fef2);
102 }
103 void Esfem::SecOrd_op::Linear_heat::
104 apply_massMatrix_to(Grid::Scal_FEfun& sfef) const{
105 FE_function& fef = sfef;
106 d_ptr -> heat_op.mass_matrix(fef);
107 }
108 void Esfem::SecOrd_op::Linear_heat::
109 mass_matrix(const Grid::Scal_FEfun& rhs, Grid::Scal_FEfun& lhs) const{
110 const FE_function& rhs_ref = rhs;
111 FE_function& lhs_ref = lhs;
112 d_ptr -> heat_op.mass_matrix(rhs_ref, lhs_ref);
113 }
114
115
116
117 Linear_heat_op::Linear_heat_op(const Esfem::Io::Parameter& p,
118 const Esfem::Grid::Grid_and_time& gt)
119 : bdf_alpha_lead_coeff {p.bdf_alphas().back()},
120 tp {gt.time_provider()}, tmp_fef {"tmp_fef", gt.fe_space()}
121 {}
122 void Linear_heat_op::operator()(const FE_function& cfef, FE_function& fef) const{
123 fef.clear();
124 const auto& df_space = fef.space();
125 for(const auto& entity : df_space){
126 const auto& geometry = entity.geometry();
127 const auto cfef_loc = cfef.localFunction(entity);
128 auto fef_loc = fef.localFunction(entity);
129 Quadrature quad {entity, cfef_loc.order() + fef_loc.order()};
130 matrixFree_assembly(tp.deltaT(), geometry, quad,
131 cfef_loc, fef_loc);
132 }
133 fef.communicate();
134 }
135 void Linear_heat_op::mass_matrix(FE_function& fef){
136 tmp_fef.assign(fef);
137 try{
138 fef.clear();
139 const auto& df_space = fef.space();
140 for(const auto& entity : df_space){
141 const auto& geometry = entity.geometry();
142 const auto& cfef_loc = tmp_fef.localFunction(entity);
143 auto fef_loc = fef.localFunction(entity);
144 Quadrature quad {entity, fef_loc.order()};
145 massMatrixFree_assembly(geometry, quad, cfef_loc, fef_loc);
146 }
147 }
148 catch(...){
149 fef.assign(tmp_fef);
150 throw std::runtime_error{"Error in Linear_heat::mass_matrix()"};
151 }
152 }
153 void Linear_heat_op::mass_matrix(const FE_function& rhs, FE_function& lhs) const{
154 lhs.clear();
155 const auto& df_space = lhs.space();
156 for(const auto& entity : df_space){
157 const auto& geometry = entity.geometry();
158 const auto& cfef_loc = rhs.localFunction(entity);
159 auto fef_loc = lhs.localFunction(entity);
160 Quadrature quad {entity, fef_loc.order() + cfef_loc.order()};
161 massMatrixFree_assembly(geometry, quad, cfef_loc, fef_loc);
162 }
163 }
164
165 void matrixFree_assembly(const double dT, const Geometry& g,
166 const Quadrature& q,
167 const FE_function::LocalFunctionType& cf,
168 FE_function::LocalFunctionType& f){
169 static_assert(Domain::dimension == 3, "Bad domain dimension.");
170 static_assert(Range::dimension == 1, "Bad range dimension.");
171
172 for(std::size_t pt = 0; pt < q.nop(); ++pt){
173
174 auto Mu = mass_matrix(pt, q, cf);
175 auto Au = stiffness_matrix(pt, q, cf);
176
177 const auto& x = q.point(pt);
178 const auto weight = q.weight(pt) * g.integrationElement(x);
179
180 Mu *= weight;
181 Au *= dT * weight;
182 f.axpy(q[pt], Mu, Au);
183 }
184 }
185 void massMatrixFree_assembly(const Geometry& g,
186 const Quadrature& q,
187 const FE_function::LocalFunctionType& cf,
188 FE_function::LocalFunctionType& f){
189 static_assert(Domain::dimension == 3, "Bad domain dimension.");
190 static_assert(Range::dimension == 1, "Bad range dimension.");
191
192 for(std::size_t pt = 0; pt < q.nop(); ++pt){
193
194 auto Mu = mass_matrix(pt, q, cf);
195
196 const auto& x = q.point(pt);
197 const auto weight = q.weight(pt) * g.integrationElement(x);
198
199 Mu *= weight;
200 f.axpy(q[pt], Mu);
201 }
202 }
203
204
205