This source file includes following definitions.
- assemble_and_addScaled_to
- assemble
- evaluate
- assemble_RHS
- massMatrix_for_entity
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16 #include <cmath>
17 #include <config.h>
18 #include <dune/fem/quadrature/cachingquadrature.hh>
19 #include "secOrd_op_rhs_u.h"
20 #include "grid.h"
21 #include "io_parameter.h"
22
23 #ifdef DEBUG
24 #include <iostream>
25 #endif
26
27 using FE_function = Esfem::Grid::Scal_FEfun::Dune_FEfun;
28 using Geometry
29 = FE_function::DiscreteFunctionSpaceType::IteratorType::Entity::Geometry;
30 using Grid_part = FE_function::GridPartType;
31 using Quadrature = Dune::Fem::CachingQuadrature<Grid_part, 0>;
32
33 class RHS_data_u
34 : public Dune::Fem::Function<Esfem::Grid::Grid_and_time::Function_space, RHS_data_u>
35 {
36 public:
37 using Base = Esfem::Grid::Grid_and_time::Function_space;
38 using Domain = Base::DomainType;
39 using Range = Base::RangeType;
40
41 RHS_data_u() = delete;
42 explicit RHS_data_u(const Dune::Fem::TimeProviderBase&,
43 const Esfem::Io::Parameter&);
44 RHS_data_u(const RHS_data_u&) = delete;
45 RHS_data_u& operator=(const RHS_data_u&) = delete;
46
47 void evaluate(const Domain&, Range&) const;
48 Range operator()(const Domain&) const;
49 private:
50 const Dune::Fem::TimeProviderBase& tp;
51 const double a;
52 const double gamma;
53 };
54
55 void assemble_RHS(const RHS_data_u&, FE_function&);
56 void matrixFree_assembly(const Dune::Fem::TimeProviderBase&, const Geometry&,
57 const Quadrature&, FE_function::LocalFunctionType&);
58 void massMatrix_for_entity(const Geometry&, const Quadrature&, const RHS_data_u&,
59 FE_function::LocalFunctionType&);
60
61
62
63
64 struct Esfem::SecOrd_op::Rhs_u::Data{
65 RHS_data_u rhs;
66 const Dune::Fem::TimeProviderBase& tp;
67 const Grid::Grid_and_time::FE_space& fe_space;
68 Data(const Io::Parameter& p, const Grid::Grid_and_time& gt)
69 : rhs {gt.time_provider(), p}, tp {gt.time_provider()},
70 fe_space {gt.fe_space()}
71 {}
72 };
73
74 Esfem::SecOrd_op::Rhs_u::Rhs_u(const Io::Parameter& p, const Grid::Grid_and_time& gt){
75 try{
76 d_ptr = new Data {p, gt};
77 }
78 catch(std::exception&){
79 std::throw_with_nested(std::logic_error
80 {"Error in constructor of Grid_and_time."});
81 }
82 catch(...){
83 throw std::logic_error {"Unknown error in constructor of Grid_and_time."};
84 }
85 }
86 Esfem::SecOrd_op::Rhs_u::~Rhs_u(){
87 delete d_ptr;
88 d_ptr = nullptr;
89 #ifdef DEBUG
90 std::cerr << "~Rhs_u(): delete d_ptr.\n";
91 #endif
92 }
93 void Esfem::SecOrd_op::Rhs_u::assemble_and_addScaled_to(Grid::Scal_FEfun& fef) const{
94 static FE_function tmp {"local_tmp", d_ptr -> fe_space};
95 assemble_RHS(d_ptr -> rhs, tmp);
96 FE_function& dune_fef = fef;
97 dune_fef.axpy(d_ptr -> tp.deltaT(), tmp);
98 }
99 void Esfem::SecOrd_op::Rhs_u::assemble(Grid::Scal_FEfun& fef) const{
100 FE_function& dune_fef = fef;
101 assemble_RHS(d_ptr -> rhs, dune_fef);
102 }
103
104
105
106
107 RHS_data_u::RHS_data_u(const Dune::Fem::TimeProviderBase& tpb,
108 const Esfem::Io::Parameter& p)
109 : tp {tpb}, a {p.tg_a()}, gamma {p.tg_gamma()}
110 {}
111 void RHS_data_u::evaluate(const Domain& d, Range& r) const{
112 static_assert(Domain::dimension == 3, "Bad domain dimension.");
113 static_assert(Range::dimension == 1, "Bad range dimension.");
114 const double x = d[0];
115 const double y = d[1];
116 const double z = d[2];
117 const double t = tp.time();
118 const double g = gamma;
119 const double u = std::exp(-6. * t) * x * y;
120
121 r = u
122
123
124
125
126
127
128
129 ;
130 }
131 RHS_data_u::Range RHS_data_u::operator()(const Domain& d) const{
132 Range r {0};
133 evaluate(d,r);
134 return r;
135 }
136
137 void assemble_RHS(const RHS_data_u& rhs_fun, FE_function& fef){
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 Quadrature quad {entity, 2 * df_space.order() + 1};
143 auto fef_local = fef.localFunction(entity);
144 massMatrix_for_entity(geometry, quad, rhs_fun, fef_local);
145 }
146 fef.communicate();
147 }
148
149 void massMatrix_for_entity(const Geometry& g, const Quadrature& q,
150 const RHS_data_u& rhs_fun,
151 FE_function::LocalFunctionType& f_loc){
152 for(std::size_t pt = 0; pt < q.nop(); ++pt){
153 const auto& x = q.point(pt);
154 RHS_data_u::Range fx {rhs_fun(g.global(x))};
155 fx *= q.weight(pt) * g.integrationElement(x);
156 f_loc.axpy(q[pt], fx);
157 }
158 }
159
160
161
162