util.hh 7.35 KB
Newer Older
Lukas Riedel's avatar
Lukas Riedel committed
1 2
#ifndef DUNE_DORIE_UTIL_HH
#define DUNE_DORIE_UTIL_HH
3

4 5 6 7 8 9 10 11 12
#include <iostream>
#include <memory>
#include <vector>

#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>

#include <dune/geometry/type.hh>

Santiago Ospina's avatar
Santiago Ospina committed
13 14 15 16
#include <dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
#include <dune/pdelab/constraints/p0.hh>
#include <dune/pdelab/finiteelementmap/p0fem.hh>
#include <dune/pdelab/finiteelementmap/qkdg.hh>
17 18
#include <dune/pdelab/finiteelementmap/opbfem.hh>

19 20 21
namespace Dune{
namespace Dorie{

22
/// Return the estimation of entries per matrix row for the spatial GridOperator
23
/** This supposedly decreases matrix assembly time.
24 25
 *  The values specify the *blocks* per row. DG assembles one block for the
 *  actual element and one for each of its neighbors.
26 27
 *  \param dim Spatial dimension
 *  \param geo Geometry type of grid entities
28
 *  \return Estimated number of blocks per matrix row
29
 */
30
template<typename R = std::size_t>
31
constexpr R estimate_mbe_entries (const int dim, const Dune::GeometryType::BasicType geo)
32
{
33 34
	if (geo==Dune::GeometryType::BasicType::cube){
		return 2*dim + 1;
35
	}
36 37
	else if (geo==Dune::GeometryType::BasicType::simplex){
		return dim + 2;
38 39
	}
	std::cerr << "Cannot provide MBE entry estimation for given dimension and/or GeometryType!" << std::endl;
40
	return 1;
41 42
}

43 44 45
/// Provide types and construction of the GFS for the linear solver
template<typename GridView, typename RF, Dune::GeometryType::BasicType GeometryType>
struct LinearSolverGridFunctionSpaceHelper
46 47
{};

48
/// Provide types and construction of the GridFunctionSpace
49
template<typename GridView, typename RF>
50
struct LinearSolverGridFunctionSpaceHelper<GridView,RF,Dune::GeometryType::BasicType::simplex>
51 52
{
private:
53
	static constexpr int dim = GridView::dimension;
54 55 56 57 58 59
	using DF = typename GridView::Grid::ctype;

public:
	/// Entity set of the GFS
	using ES = Dune::PDELab::OverlappingEntitySet<GridView>;
	/// FiniteElementMap type of GFS
60
	using FEM = typename Dune::PDELab::P0LocalFiniteElementMap<DF,RF,dim>;
61 62 63 64
	/// Constraints type of the GFS
	using CON = Dune::PDELab::P0ParallelConstraints;
	/// GFS type
	using Type = typename Dune::PDELab::GridFunctionSpace<ES,FEM,CON,
Lukas Riedel's avatar
Lukas Riedel committed
65
		Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed, 1>>;
66 67 68 69 70 71 72
	/// GFS Constraints Container type
	using CC = typename Type::template ConstraintsContainer<RF>::Type;

	/// create GFS from GridView
	static Type create (const GridView& gv)
	{
		ES es(gv);
Lukas Riedel's avatar
Lukas Riedel committed
73
		auto fem = std::make_shared<FEM>(Dune::GeometryTypes::simplex(dim));
74 75 76 77 78 79 80
		auto con = std::make_shared<CON>();
		return Type(es,fem,con);
	}
};

/// Provide types and construction of the GridFunctionSpace
template<typename GridView, typename RF>
81
struct LinearSolverGridFunctionSpaceHelper<GridView,RF,Dune::GeometryType::BasicType::cube>
82 83 84 85 86 87 88 89 90
{
private:
	static constexpr int dim = GridView::dimension;
	using DF = typename GridView::Grid::ctype;

public:
	/// Entity set of the GFS
	using ES = Dune::PDELab::OverlappingEntitySet<GridView>;
	/// FiniteElementMap type of GFS
91
	using FEM = typename Dune::PDELab::P0LocalFiniteElementMap<DF,RF,dim>;
92 93 94
	/// Constraints type of the GFS
	using CON = Dune::PDELab::P0ParallelConstraints;
	/// GFS type
95
	using Type = typename Dune::PDELab::GridFunctionSpace<ES,FEM,CON,
Lukas Riedel's avatar
Lukas Riedel committed
96
		Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed, 1>>;
97 98 99 100 101 102
	/// GFS Constraints Container type
	using CC = typename Type::template ConstraintsContainer<RF>::Type;

	/// create GFS from GridView
	static Type create (const GridView& gv)
	{
103
		ES es(gv);
Lukas Riedel's avatar
Lukas Riedel committed
104
		auto fem = std::make_shared<FEM>(Dune::GeometryTypes::cube(dim));
105
		auto con = std::make_shared<CON>();
106
		return Type(es,fem,con);
107 108 109
	}
};

110 111 112 113 114
/// Provide types and construction of the GridFunctionSpace. Has to be specialized.
template<typename GridView, typename RF, int order, Dune::GeometryType::BasicType GeometryType>
struct GridFunctionSpaceHelper
{};

115 116 117
/// Provide types and construction of the GridFunctionSpace
/*  \todo add constraints specialization depending on grid
 */
118
template<typename GridView, typename RF, int order>
119
struct GridFunctionSpaceHelper<GridView,RF,order,Dune::GeometryType::BasicType::simplex>
120
{
121
private:
122
	static constexpr int dim = GridView::dimension;
123
	using DF = typename GridView::Grid::ctype;
124 125

public:
126
	/// Entity set of the GFS
127
	using ES = Dune::PDELab::OverlappingEntitySet<GridView>;
128
	/// FiniteElementMap type of GFS
129
	using FEM = typename Dune::PDELab::OPBLocalFiniteElementMap<DF, RF, order, dim, Dune::GeometryType::BasicType::simplex>;
130
	/// Constraints type of the GFS
131
	using CON = Dune::PDELab::P0ParallelConstraints;
132
	/// GFS type
133
	using Type = typename Dune::PDELab::GridFunctionSpace<ES,FEM,CON,
Lukas Riedel's avatar
Lukas Riedel committed
134 135 136 137 138 139
		Dune::PDELab::ISTL::VectorBackend<
			Dune::PDELab::ISTL::Blocking::fixed,
			Dune::PB::pk_size(order, dim)
		>
	>;

140
	/// GFS Constraints Container type
141 142
	using CC = typename Type::template ConstraintsContainer<RF>::Type;

143
	/// create GFS from GridView
144
	static Type create (const GridView& gv)
145
	{
146
		ES es(gv);
147
		auto fem = std::make_shared<FEM>();
148 149
		auto con = std::make_shared<CON>();
		return Type(es,fem,con);
150 151 152
	}
};

153
/// Provide types and construction of the GridFunctionSpace
154 155 156
template<typename GridView, typename RF, int order>
struct GridFunctionSpaceHelper<GridView,RF,order,Dune::GeometryType::BasicType::cube>
{
157
private:
158 159
	static constexpr int dim = GridView::dimension;
	using DF = typename GridView::Grid::ctype;
160 161

public:
162
	/// Entity set of the GFS
163
	using ES = Dune::PDELab::OverlappingEntitySet<GridView>;
164
	/// FiniteElementMap type of GFS
165
	using FEM = typename Dune::PDELab::QkDGLocalFiniteElementMap<DF,RF,order,dim>;
166
	/// Constraints type of the GFS
167
	using CON = Dune::PDELab::P0ParallelConstraints;
168
	/// GFS type
169
	using Type = typename Dune::PDELab::GridFunctionSpace<ES,FEM,CON,
170
		Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed>
Lukas Riedel's avatar
Lukas Riedel committed
171
	>;
172
	/// GFS Constraints Container type
173 174
	using CC = typename Type::template ConstraintsContainer<RF>::Type;

175
	/// create GFS from GridView
176
	static Type create (const GridView& gv)
177
	{
178
		ES es(gv);
179
		auto fem = std::make_shared<FEM>();
180 181
		auto con = std::make_shared<CON>();
		return Type(es,fem,con);
182
	}
183 184
};

185
/// Traits struct defining basic types based on grid an geometry
186
template< template<int> class GridType, Dune::GeometryType::BasicType GeometryType, int dimension, int order , bool output, bool adaptivity>
187
struct BaseTraits
188
{
189
	static constexpr int dim = dimension;
190
	static constexpr int fem_order = order;
191
	static constexpr bool write_output = output;
192
	static constexpr bool adapt_grid = adaptivity;
Lukas Riedel's avatar
Lukas Riedel committed
193
	static constexpr Dune::GeometryType::BasicType GridGeometryType = GeometryType;
194

195
	using RF = double;
196 197
	using RangeField = RF;
	using Array = std::vector<RangeField>;
198 199
	using Scalar = Dune::FieldVector<RF,1>;
	using Vector = Dune::FieldVector<RF,dim>;
200 201 202
	using Tensor = Dune::FieldMatrix<RF,dim,dim>;
	using Index = unsigned int;
	using IndexArray = Dune::FieldVector<Index,dim>;
203 204 205 206 207 208 209

	using Grid = GridType<dim>;
	using DomainField = typename Grid::ctype;
	using DF = DomainField;
	using Domain = Dune::FieldVector<DF,dim>;
	using IntersectionDomain = Dune::FieldVector<DF,dim-1>;

210
	using GV = typename Grid::LeafGridView;
211
	using GridView = GV;
212 213
	using Element = typename GV::Traits::template Codim<0>::Entity;
	using Intersection = typename GV::Traits::Intersection;
Lukas Riedel's avatar
Lukas Riedel committed
214
	//using GFS = GridFunctionSpaceHelper<GV,RF,order,GeometryType>;
215 216
	//using AdaptivityGFS = GridFunctionSpaceHelper<GV,RF,0,GeometryType>;
	//using LinearSolverGFS = GridFunctionSpaceHelper<GV,RF,0,GeometryType>;
217

218 219
};

220 221


222 223 224
}
}

Lukas Riedel's avatar
Lukas Riedel committed
225
#endif // DUNE_DORIE_UTIL_HH