Merge branch '80-estimate-of-matrix-backend-entry-numbers-is-wrong' into 'master'

Resolve "Estimate of matrix backend entry numbers is wrong"

See merge request !75

Closes #80

CHANGELOG.md

@@ -16,6 +16,9 @@
     cause the corresponding test to fail if the designated output directory did
     not exist.
 * Removed CI artifacts of `test` stage from `deploy` jobs.
+* Correctly estimate the number of row entries in the matrix backend based on
+    the DG discretization and proper static blocking. This can lead to a
+    performance improvement for large problems.
 
 
 ## 1.1.0 (2018-07-27)

dune/dorie/interface/util.hh

@@ -19,28 +19,22 @@
 namespace Dune{
 namespace Dorie{
 
-/// Return the estimation of entries per matrix row for the spatial GridOperator.
+/// 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().
+ *  The values specify the *blocks* per row. DG assembles one block for the
+ *  actual element and one for each of its neighbors.
  *  \param dim Spatial dimension
  *  \param geo Geometry type of grid entities
- *  \return Estimated entries per matrix row
+ *  \return Estimated number of blocks per matrix row
  */
 template<typename R = std::size_t>
-R estimate_mbe_entries (const int dim, const Dune::GeometryType::BasicType geo)
+constexpr R estimate_mbe_entries (const int dim, const Dune::GeometryType::BasicType geo)
 {
-	if(geo==Dune::GeometryType::BasicType::cube){
-		if(dim==2) 		
-			return 5;
-		else if(dim==3)
-			return 7;
+	if (geo==Dune::GeometryType::BasicType::cube){
+		return 2*dim + 1;
 	}
-	else if(geo==Dune::GeometryType::BasicType::simplex){
-		if(dim==2)
-			return 13;
-		else if(dim==3) 
-			return 20;
+	else if (geo==Dune::GeometryType::BasicType::simplex){
+		return dim + 2;
 	}
 	std::cerr << "Cannot provide MBE entry estimation for given dimension and/or GeometryType!" << std::endl;
 	return 1;