util.hh 7.25 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
#include <dune/dorie/common/finite_element/pk_dg_fem.hh>
18

19 20 21
namespace Dune{
namespace Dorie{

22
/**
23
 * @brief      Enum class for output policy.
24
 * @details    It defines when a model class should produce an output.
25
 */
26
enum class OutputPolicy {None,EndOfRichardsStep,EndOfTransportStep};
27 28

/**
29 30
 * @brief      Enum class for output policy.
 * @details    It defines which variable is the target in order to mark the grid
31
 */
32
enum class AdaptivityPolicy {None,WaterFlux,SoluteFlux};
33

34
/// Return the estimation of entries per matrix row for the spatial GridOperator
35
/** This supposedly decreases matrix assembly time.
36 37
 *  The values specify the *blocks* per row. DG assembles one block for the
 *  actual element and one for each of its neighbors.
38 39
 *  \param dim Spatial dimension
 *  \param geo Geometry type of grid entities
40
 *  \return Estimated number of blocks per matrix row
41
 */
42
template<typename R = std::size_t>
43
constexpr R estimate_mbe_entries (const int dim, const Dune::GeometryType::BasicType geo)
44
{
45 46
	if (geo==Dune::GeometryType::BasicType::cube){
		return 2*dim + 1;
47
	}
48 49
	else if (geo==Dune::GeometryType::BasicType::simplex){
		return dim + 2;
50
	}
51 52 53

	throw std::runtime_error("Cannot provide matrix backend entry estimate "
							 "for given geometry type!");
54 55
}

56 57 58
/// Provide types and construction of the GFS for the linear solver
template<typename GridView, typename RF, Dune::GeometryType::BasicType GeometryType>
struct LinearSolverGridFunctionSpaceHelper
59 60
{};

61
/// Provide types and construction of the GridFunctionSpace
62
template<typename GridView, typename RF>
63
struct LinearSolverGridFunctionSpaceHelper<GridView,RF,Dune::GeometryType::BasicType::simplex>
64 65
{
private:
66
	static constexpr int dim = GridView::dimension;
67 68 69 70 71 72
	using DF = typename GridView::Grid::ctype;

public:
	/// Entity set of the GFS
	using ES = Dune::PDELab::OverlappingEntitySet<GridView>;
	/// FiniteElementMap type of GFS
73
	using FEM = typename Dune::PDELab::P0LocalFiniteElementMap<DF,RF,dim>;
74 75 76 77
	/// Constraints type of the GFS
	using CON = Dune::PDELab::P0ParallelConstraints;
	/// GFS type
	using Type = typename Dune::PDELab::GridFunctionSpace<ES,FEM,CON,
78
		Dune::PDELab::ISTL::VectorBackend<>>;
79 80 81 82 83 84 85
	/// 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
86
		auto fem = std::make_shared<FEM>(Dune::GeometryTypes::simplex(dim));
87 88 89 90 91 92 93
		auto con = std::make_shared<CON>();
		return Type(es,fem,con);
	}
};

/// Provide types and construction of the GridFunctionSpace
template<typename GridView, typename RF>
94
struct LinearSolverGridFunctionSpaceHelper<GridView,RF,Dune::GeometryType::BasicType::cube>
95 96 97 98 99 100 101 102 103
{
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
104
	using FEM = typename Dune::PDELab::P0LocalFiniteElementMap<DF,RF,dim>;
105 106 107
	/// Constraints type of the GFS
	using CON = Dune::PDELab::P0ParallelConstraints;
	/// GFS type
108
	using Type = typename Dune::PDELab::GridFunctionSpace<ES,FEM,CON,
109
		Dune::PDELab::ISTL::VectorBackend<>>;
110 111 112 113 114 115
	/// GFS Constraints Container type
	using CC = typename Type::template ConstraintsContainer<RF>::Type;

	/// create GFS from GridView
	static Type create (const GridView& gv)
	{
116
		ES es(gv);
Lukas Riedel's avatar
Lukas Riedel committed
117
		auto fem = std::make_shared<FEM>(Dune::GeometryTypes::cube(dim));
118
		auto con = std::make_shared<CON>();
119
		return Type(es,fem,con);
120 121 122
	}
};

123 124 125 126 127
/// Provide types and construction of the GridFunctionSpace. Has to be specialized.
template<typename GridView, typename RF, int order, Dune::GeometryType::BasicType GeometryType>
struct GridFunctionSpaceHelper
{};

128 129 130
/// Provide types and construction of the GridFunctionSpace
/*  \todo add constraints specialization depending on grid
 */
131
template<typename GridView, typename RF, int order>
132
struct GridFunctionSpaceHelper<GridView,RF,order,Dune::GeometryType::BasicType::simplex>
133
{
134
private:
135
	static constexpr int dim = GridView::dimension;
136
	using DF = typename GridView::Grid::ctype;
137 138

public:
139
	/// Entity set of the GFS
140
	using ES = Dune::PDELab::OverlappingEntitySet<GridView>;
141
	/// FiniteElementMap type of GFS
142
	using FEM = typename Dune::Dorie::PkDGLocalFiniteElementMap<DF, RF, dim, order>;
143
	/// Constraints type of the GFS
144
	using CON = Dune::PDELab::P0ParallelConstraints;
145
	/// GFS type
146
	using Type = typename Dune::PDELab::GridFunctionSpace<ES,FEM,CON,
Lukas Riedel's avatar
Lukas Riedel committed
147 148 149 150 151 152
		Dune::PDELab::ISTL::VectorBackend<
			Dune::PDELab::ISTL::Blocking::fixed,
			Dune::PB::pk_size(order, dim)
		>
	>;

153
	/// GFS Constraints Container type
154 155
	using CC = typename Type::template ConstraintsContainer<RF>::Type;

156
	/// create GFS from GridView
157
	static Type create (const GridView& gv)
158
	{
159
		ES es(gv);
160
		auto fem = std::make_shared<FEM>();
161 162
		auto con = std::make_shared<CON>();
		return Type(es,fem,con);
163 164 165
	}
};

166
/// Provide types and construction of the GridFunctionSpace
167 168 169
template<typename GridView, typename RF, int order>
struct GridFunctionSpaceHelper<GridView,RF,order,Dune::GeometryType::BasicType::cube>
{
170
private:
171 172
	static constexpr int dim = GridView::dimension;
	using DF = typename GridView::Grid::ctype;
173 174

public:
175
	/// Entity set of the GFS
176
	using ES = Dune::PDELab::OverlappingEntitySet<GridView>;
177
	/// FiniteElementMap type of GFS
178
	using FEM = typename Dune::PDELab::QkDGLocalFiniteElementMap<DF,RF,order,dim>;
179
	/// Constraints type of the GFS
180
	using CON = Dune::PDELab::P0ParallelConstraints;
181
	/// GFS type
182
	using Type = typename Dune::PDELab::GridFunctionSpace<ES,FEM,CON,
183
		Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed>
Lukas Riedel's avatar
Lukas Riedel committed
184
	>;
185
	/// GFS Constraints Container type
186 187
	using CC = typename Type::template ConstraintsContainer<RF>::Type;

188
	/// create GFS from GridView
189
	static Type create (const GridView& gv)
190
	{
191
		ES es(gv);
192
		auto fem = std::make_shared<FEM>();
193 194
		auto con = std::make_shared<CON>();
		return Type(es,fem,con);
195
	}
196 197
};

198
/// Traits struct defining basic types based on grid an geometry
199
template<class GridType, Dune::GeometryType::BasicType GeometryType>
200
struct BaseTraits
201
{
202
	static constexpr int dim = GridType::dimension;
Lukas Riedel's avatar
Lukas Riedel committed
203
	static constexpr Dune::GeometryType::BasicType GridGeometryType = GeometryType;
204

205 206 207
	using TF = double;
	using TimeField = TF;

208
	using RF = double;
209 210
	using RangeField = RF;
	using Array = std::vector<RangeField>;
211 212
	using Scalar = Dune::FieldVector<RF,1>;
	using Vector = Dune::FieldVector<RF,dim>;
213 214 215
	using Tensor = Dune::FieldMatrix<RF,dim,dim>;
	using Index = unsigned int;
	using IndexArray = Dune::FieldVector<Index,dim>;
216

217
	using Grid = GridType;
218 219 220 221 222
	using DomainField = typename Grid::ctype;
	using DF = DomainField;
	using Domain = Dune::FieldVector<DF,dim>;
	using IntersectionDomain = Dune::FieldVector<DF,dim-1>;

223
	using GV = typename Grid::LeafGridView;
224
	using GridView = GV;
225 226
	using Element = typename GV::Traits::template Codim<0>::Entity;
	using Intersection = typename GV::Traits::Intersection;
227 228 229 230 231
};

}
}

232
#endif // DUNE_DORIE_UTIL_HH