util.hh 6.65 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
namespace Dune{
namespace Dorie{

7 8 9 10 11 12 13 14
/// Return the estimation of entries per matrix row for the spatial GridOperator.
/** This supposedly decreases matrix assembly time.
 *  The values are empirical. The actual matrix statistics can be displayed
 *  by assembling the IGO Jacobian and calling patternStatistics().
 *  \param dim Spatial dimension
 *  \param geo Geometry type of grid entities
 *  \return Estimated entries per matrix row
 */
15
template<typename R = std::size_t>
16 17 18 19
R estimate_mbe_entries (const int dim, const Dune::GeometryType::BasicType geo)
{
	if(geo==Dune::GeometryType::BasicType::cube){
		if(dim==2) 		
20
			return 5;
21
		else if(dim==3)
22
			return 7;
23 24 25
	}
	else if(geo==Dune::GeometryType::BasicType::simplex){
		if(dim==2)
26
			return 13;
27
		else if(dim==3) 
28
			return 20;
29 30
	}
	std::cerr << "Cannot provide MBE entry estimation for given dimension and/or GeometryType!" << std::endl;
31
	return 1;
32 33
}

34
/// Provide types and construction of the GridFunctionSpace. Has to be specialized.
35
template<typename GridView, typename RF, int order, Dune::GeometryType::BasicType GeometryType>
36 37 38
struct GridFunctionSpaceHelper
{};

39 40 41 42
template<typename GridView, typename RF>
struct GridFunctionSpaceHelper<GridView,RF,0,Dune::GeometryType::BasicType::simplex>
{
private:
43
	static constexpr int dim = GridView::dimension;
44 45 46 47 48 49
	using DF = typename GridView::Grid::ctype;

public:
	/// Entity set of the GFS
	using ES = Dune::PDELab::OverlappingEntitySet<GridView>;
	/// FiniteElementMap type of GFS
50
	using FEM = typename Dune::PDELab::P0LocalFiniteElementMap<DF,RF,dim>;
51 52 53 54 55 56 57 58 59 60 61 62
	/// Constraints type of the GFS
	using CON = Dune::PDELab::P0ParallelConstraints;
	/// GFS type
	using Type = typename Dune::PDELab::GridFunctionSpace<ES,FEM,CON,
		Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::fixed>>;
	/// 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);
63 64 65
		Dune::GeometryType geo;
		geo.makeSimplex(dim);
		auto fem = std::make_shared<FEM>(geo);
66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82
		auto con = std::make_shared<CON>();
		return Type(es,fem,con);
	}
};

/// Provide types and construction of the GridFunctionSpace
template<typename GridView, typename RF>
struct GridFunctionSpaceHelper<GridView,RF,0,Dune::GeometryType::BasicType::cube>
{
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
83
	using FEM = typename Dune::PDELab::P0LocalFiniteElementMap<DF,RF,dim>;
84 85 86
	/// Constraints type of the GFS
	using CON = Dune::PDELab::P0ParallelConstraints;
	/// GFS type
87
	using Type = typename Dune::PDELab::GridFunctionSpace<ES,FEM,CON,
88 89 90 91 92 93 94
		Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::fixed> >;
	/// GFS Constraints Container type
	using CC = typename Type::template ConstraintsContainer<RF>::Type;

	/// create GFS from GridView
	static Type create (const GridView& gv)
	{
95
		ES es(gv);
96 97 98
		Dune::GeometryType geo;
		geo.makeCube(dim);
		auto fem = std::make_shared<FEM>(geo);
99
		auto con = std::make_shared<CON>();
100
		return Type(es,fem,con);
101 102 103
	}
};

104 105 106
/// Provide types and construction of the GridFunctionSpace
/*  \todo add constraints specialization depending on grid
 */
107
template<typename GridView, typename RF, int order>
108
struct GridFunctionSpaceHelper<GridView,RF,order,Dune::GeometryType::BasicType::simplex>
109
{
110
private:
111
	using DF = typename GridView::Grid::ctype;
112 113

public:
114
	/// Entity set of the GFS
115
	using ES = Dune::PDELab::OverlappingEntitySet<GridView>;
116
	/// FiniteElementMap type of GFS
117 118
	using FEM = typename Dune::PDELab::PkLocalFiniteElementMap<ES,DF,RF,order>;
	/// Constraints type of the GFS
119
	using CON = Dune::PDELab::P0ParallelConstraints;
120
	/// GFS type
121 122
	using Type = typename Dune::PDELab::GridFunctionSpace<ES,FEM,CON,
		Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::fixed>>;
123
	/// GFS Constraints Container type
124 125
	using CC = typename Type::template ConstraintsContainer<RF>::Type;

126
	/// create GFS from GridView
127
	static Type create (const GridView& gv)
128
	{
129 130 131 132
		ES es(gv);
		auto fem = std::make_shared<FEM>(es);
		auto con = std::make_shared<CON>();
		return Type(es,fem,con);
133 134 135
	}
};

136
/// Provide types and construction of the GridFunctionSpace
137 138 139
template<typename GridView, typename RF, int order>
struct GridFunctionSpaceHelper<GridView,RF,order,Dune::GeometryType::BasicType::cube>
{
140
private:
141 142
	static constexpr int dim = GridView::dimension;
	using DF = typename GridView::Grid::ctype;
143 144

public:
145
	/// Entity set of the GFS
146
	using ES = Dune::PDELab::OverlappingEntitySet<GridView>;
147
	/// FiniteElementMap type of GFS
148
	using FEM = typename Dune::PDELab::QkDGLocalFiniteElementMap<DF,RF,order,dim>;
149
	/// Constraints type of the GFS
150
	using CON = Dune::PDELab::P0ParallelConstraints;
151
	/// GFS type
152 153
	using Type = typename Dune::PDELab::GridFunctionSpace<ES,FEM,CON,
		Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::fixed> >;
154
	/// GFS Constraints Container type
155 156
	using CC = typename Type::template ConstraintsContainer<RF>::Type;

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

167
/// Traits struct defining basic types based on grid an geometry
168
template< template<int> class GridType, Dune::GeometryType::BasicType GeometryType, int dimension, int order , bool output, bool adaptivity>
169
struct BaseTraits
170
{
171
	static constexpr int dim = dimension;
172
	static constexpr int fem_order = order;
173
	static constexpr bool write_output = output;
174
	static constexpr bool adapt_grid = adaptivity;
Lukas Riedel's avatar
Lukas Riedel committed
175
	static constexpr Dune::GeometryType::BasicType GridGeometryType = GeometryType;
176

177
	using RF = double;
178 179
	using RangeField = RF;
	using Array = std::vector<RangeField>;
180 181
	using Scalar = Dune::FieldVector<RF,1>;
	using Vector = Dune::FieldVector<RF,dim>;
182 183 184
	using Tensor = Dune::FieldMatrix<RF,dim,dim>;
	using Index = unsigned int;
	using IndexArray = Dune::FieldVector<Index,dim>;
185 186 187 188 189 190 191

	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>;

192
	using GV = typename Grid::LeafGridView;
193
	using GridView = GV;
194 195
	using Element = typename GV::Traits::template Codim<0>::Entity;
	using Intersection = typename GV::Traits::Intersection;
Lukas Riedel's avatar
Lukas Riedel committed
196
	//using GFS = GridFunctionSpaceHelper<GV,RF,order,GeometryType>;
197 198
	//using AdaptivityGFS = GridFunctionSpaceHelper<GV,RF,0,GeometryType>;
	//using LinearSolverGFS = GridFunctionSpaceHelper<GV,RF,0,GeometryType>;
199

200 201
};

202 203


204 205 206
}
}

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