Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
// This file is part of the dune-gdt project:
// https://github.com/dune-community/dune-gdt
// Copyright 2010-2018 dune-gdt developers and contributors. All rights reserved.
// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
// or GPL-2.0+ (http://opensource.org/licenses/gpl-license)
// with "runtime exception" (http://www.dune-project.org/license.html)
// Authors:
// Felix Schindler (2019)
#ifndef DUNE_GDT_TOOLS_local_mass_matrices_HH
#define DUNE_GDT_TOOLS_local_mass_matrices_HH
#include <mutex>
#include <vector>
#include <dune/xt/common/parallel/threadstorage.hh>
#include <dune/xt/la/container/common.hh>
#include <dune/xt/la/container/conversion.hh>
#include <dune/xt/la/matrix-inverter.hh>
#include <dune/xt/grid/functors/interfaces.hh>
#include <dune/xt/grid/type_traits.hh>
#include <dune/gdt/exceptions.hh>
#include <dune/gdt/local/bilinear-forms/integrals.hh>
#include <dune/gdt/local/integrands/product.hh>
#include <dune/gdt/spaces/interface.hh>
#include <dune/gdt/spaces/mapper/finite-volume.hh>
namespace Dune {
namespace GDT {
template <class GV, size_t r = 1, size_t rC = 1, class F = double, class AGV = GV>
class LocalMassMatrixProvider
: public XT::Grid::ElementFunctor<GV>
, public XT::Common::ThreadResultPropagator<
LocalMassMatrixProvider<GV, r, rC, F>,
std::map<size_t, std::pair<XT::LA::CommonDenseMatrix<F>, XT::LA::CommonDenseMatrix<F>>>,
XT::Common::merge_container<
std::map<size_t, std::pair<XT::LA::CommonDenseMatrix<F>, XT::LA::CommonDenseMatrix<F>>>>>
{
static_assert(XT::Grid::is_view<AGV>::value, "");
using ThisType = LocalMassMatrixProvider<GV, r, rC, F>;
using BaseType = XT::Grid::ElementFunctor<GV>;
using Propagator = XT::Common::ThreadResultPropagator<
LocalMassMatrixProvider<GV, r, rC, F>,
std::map<size_t, std::pair<XT::LA::CommonDenseMatrix<F>, XT::LA::CommonDenseMatrix<F>>>,
XT::Common::merge_container<
std::map<size_t, std::pair<XT::LA::CommonDenseMatrix<F>, XT::LA::CommonDenseMatrix<F>>>>>;
friend Propagator;
public:
using AssemblyGridView = AGV;
using SpaceType = SpaceInterface<GV, r, rC, F>;
using typename BaseType::E;
using typename BaseType::ElementType;
LocalMassMatrixProvider(const AssemblyGridView& grid_view, const SpaceType& space)
: Propagator(this)
, grid_view_(grid_view)
, space_(space)
, element_mapper_(grid_view_)
, instance_counter_(0)
, local_basis_(space_.basis().localize())
{}
LocalMassMatrixProvider(const ThisType& other)
: Propagator(this)
, grid_view_(other.grid_view_)
, space_(other.space_)
, element_mapper_(grid_view_)
, instance_counter_(other.instance_counter_ + 1)
, local_basis_(space_.basis().localize())
{}
void apply_local(const ElementType& element)
{
local_basis_->bind(element);
const size_t id = element_mapper_.global_index(element, 0);
const LocalElementIntegralBilinearForm<E, r, rC, F, F> local_l2_bilinear_form(
LocalElementProductIntegrand<E, r, rC, F, F>(1.));
auto matrix =
XT::LA::convert_to<XT::LA::CommonDenseMatrix<F>>(local_l2_bilinear_form.apply2(*local_basis_, *local_basis_));
auto inverse = XT::LA::invert_matrix(matrix);
local_mass_matrices_.insert(std::make_pair(id, std::make_pair(std::move(matrix), std::move(inverse))));
} // ... apply_local(...)
BaseType* copy() override final
{
return Propagator::copy_imp();
}
void finalize() override final
{
this->finalize_imp();
}
/**
* \note This is only used to merge the thread-local results after the grid walk, you are probably interested in
* local_mass_matrix() and local_mass_matrix_inverse().
*/
std::map<size_t, std::pair<XT::LA::CommonDenseMatrix<F>, XT::LA::CommonDenseMatrix<F>>> result() const
{
return local_mass_matrices_;
}
void set_result(std::map<size_t, std::pair<XT::LA::CommonDenseMatrix<F>, XT::LA::CommonDenseMatrix<F>>> res)
{
local_mass_matrices_ = res;
}
const XT::LA::CommonDenseMatrix<F>& local_mass_matrix(const ElementType& element) const
{
const size_t id = element_mapper_.global_index(element, 0);
DUNE_THROW_IF(local_mass_matrices_.count(id) == 0,
XT::Common::Exceptions::this_should_not_happen,
"Missing local mass matrix for id " << id << "!");
return local_mass_matrices_.at(id).first;
}
const XT::LA::CommonDenseMatrix<F>& local_mass_matrix_inverse(const ElementType& element) const
{
const size_t id = element_mapper_.global_index(element, 0);
DUNE_THROW_IF(local_mass_matrices_.count(id) == 0,
XT::Common::Exceptions::this_should_not_happen,
"Missing local mass matrix for id " << id << "!");
return local_mass_matrices_.at(id).second;
}
private:
const AssemblyGridView& grid_view_;
const SpaceType& space_;
const FiniteVolumeMapper<GV> element_mapper_;
size_t instance_counter_;
std::unique_ptr<typename SpaceType::GlobalBasisType::LocalizedBasisType> local_basis_;
std::map<size_t, std::pair<XT::LA::CommonDenseMatrix<F>, XT::LA::CommonDenseMatrix<F>>> local_mass_matrices_;
}; // class LocalMassMatrixProvider
} // namespace GDT
} // namespace Dune
#endif // DUNE_GDT_TOOLS_local_mass_matrices_HH