update documentation

Author: AlePalu

tutorials/tut_1.rst

@@ -92,14 +92,14 @@ We then define the forcing term as a plain :code:`ScalarField` togheter with the
        return -9*std::pow(p[0], 4) - 12*p[0]*p[0]*p[1]*p[1] + 3*p[0]*p[0] +
                2*p[1]*p[1] - 4*std::pow(p[1], 4) - 10;
    })> f;
-   auto F = integral(unit_square, QS2D6P)(f * v);
+   auto F = integral(unit_square, QS2DP4)(f * v);
 
-Observe that we explicitly require an higher order quadrature specifying the 6 points quadrature formula :code:`QS2D6P` as second argument of the :code:`integral` function. Finally, we define non-homegeneous Dirichlet boundary conditions :math:`g(\boldsymbol{x}) = 3x^2 + 2y^2` on all the boundary of the domain
+Observe that we explicitly require an higher order quadrature specifying the quadrature formula :code:`QS2DP4` for the exact integration of order 4 polynomials as second argument of the :code:`integral` function. Finally, we define non-homegeneous Dirichlet boundary conditions :math:`g(\boldsymbol{x}) = 3x^2 + 2y^2` on all the boundary of the domain
 
 .. code-block:: cpp
 		
    ScalarField<2, decltype([](const PointT& p) { return 3 * p[0] * p[0] + 2 * p[1] * p[1]; })> g;
-   DofHandler<2, 2>& dof_handler = Vh.dof_handler();
+   auto& dof_handler = Vh.dof_handler();
    dof_handler.set_dirichlet_constraint(/* on = */ BoundaryAll, /* data = */ g);
 
 Recall that Dirichlet boundary conditions are implemented as constraints on the degrees of freedom of the linear system :math:`A \boldsymbol{c} = \boldsymbol{b}` deriving form the discretization of the variational problem, and that we must later enforce them on the pair :math:`(A, \boldsymbol{b})` before solving the linear system, using the :code:`enforce_constraints` method.
@@ -172,11 +172,11 @@ The code just assembles :code:`A1` and :code:`A2`, updates the right hand side :
 	 ScalarField<2, decltype([](const PointT& p) {
 	     return -9*std::pow(p[0], 4) - 12*p[0]*p[0]*p[1]*p[1] + 3*p[0]*p[0] + 2*p[1]*p[1] - 4*std::pow(p[1], 4) - 10;
 	 })> f;
-	 auto F = integral(unit_square, QS2D6P)(f * v);
+	 auto F = integral(unit_square, QS2DP4)(f * v);
 
 	 // define dirichlet data
 	 ScalarField<2, decltype([](const PointT& p) { return 3 * p[0] * p[0] + 2 * p[1] * p[1]; })> g;
-	 DofHandler<2, 2>& dof_handler = Vh.dof_handler();
+	 auto& dof_handler = Vh.dof_handler();
 	 dof_handler.set_dirichlet_constraint(/* on = */ BoundaryAll, /* data = */ g);
 
 	 // Newton scheme initialization (solve linearized problem with initial guess u = 0)

tutorials/tut_2.rst

@@ -88,7 +88,7 @@ You can fix the time coordinate calling :code:`f.at(t)`. Subsequent calls of :co
 .. code-block:: cpp
 
    ScalarField<2, decltype([](const PointT& p) { return 0; })> g;
-   DofHandler<2, 2>& dof_handler = Vh.dof_handler();
+   auto& dof_handler = Vh.dof_handler();
    dof_handler.set_dirichlet_constraint(/* on = */ BoundaryAll, /* data = */ g);
 
 Finally, we fix the time step :math:`\Delta t`, set up room for the solution fixing the initial condition to :math:`u_0(\boldsymbol{x}) = 0` and discretizing once and for all the mass matrix :math:`M` and the stiff matrix :math:`\frac{M}{\Delta T} + \frac{1}{2} A`, togheter with the forcing term :math:`F(t)`. Since the matrix :math:`\frac{M}{\Delta T} + \frac{1}{2} A` is SPD and time-invariant, we factorize it outside the time integration loop using a Cholesky factorization:
@@ -145,7 +145,7 @@ Finally, the crank-nicolson time integration loop can start:
          })> f;
 	 // dirichlet data (homogeneous and fixed in time)
 	 ScalarField<2, decltype([](const PointT& p) { return 0; })> g;
-	 DofHandler<2, 2>& dof_handler = Vh.dof_handler();
+	 auto& dof_handler = Vh.dof_handler();
 	 dof_handler.set_dirichlet_constraint(/* on = */ BoundaryAll, /* data = */ g);
 
 	 // crank-nicolson integration
@@ -238,7 +238,7 @@ The script is mostly similar to the Crank-Nicolson time-stepping scheme implemen
 
 	 // dirichlet homoegeneous data (fixed in time)
 	 ScalarField<2, decltype([](const PointT& p) { return 0; })> g_D;
-	 DofHandler<2, 2>& dof_handler = Vh.dof_handler();
+	 auto& dof_handler = Vh.dof_handler();
 	 dof_handler.set_dirichlet_constraint(/* on = */ 0, /* data = */ g_D);
 	 // neumann inflow data
 	 SpaceTimeField<2, decltype([](const PointT& p, double t) { return p[1] * (1 - p[1]) * t * (0.5 - t); })> g_N;

tutorials/tut_5.rst

@@ -93,7 +93,7 @@ In order to apply SUPG, let us first define the stabilization parameter :math:`\
 
    CellDiameterDescriptor h_k(unit_square);
    double delta = 2.85;
-   double b_norm = integral(unit_square, QS2D3P)(b.norm());   // transport field norm
+   double b_norm = integral(unit_square, QS2DP2)(b.norm());   // transport field norm
    
    auto tau_k  = 0.5 * delta * h_k / b_norm;   // stabilization parameter
 
@@ -111,9 +111,9 @@ In order to apply SUPG, let us first define the stabilization parameter :math:`\
 
    .. code-block:: cpp
 
-      double b_norm = integral(unit_square, QS2D3P)(b.norm());
+      double b_norm = integral(unit_square, QS2DP2)(b.norm());
 
-   in the definition of the stablization parameter computes :math:`\int_{[0,1]^2} \| b \|_2` using a two dimensional 3 point simplex rule (:code:`QS2D3P`).
+   in the definition of the stablization parameter computes :math:`\int_{[0,1]^2} \| b \|_2` using a two dimensional 3 point simplex rule (:code:`QS2DP2`).
 
 We have now all the ingredients to assemble the stabilized bilinear forms :math:`\tilde a(u,v)` and :math:`\tilde F(v)`:
 
@@ -163,7 +163,7 @@ Upon discretization and imposition of boundary conditions, the solution of the d
 	 // SUPG correction
 	 CellDiameterDescriptor h_k(unit_square);
 	 double delta = 2.85;
-	 double b_norm = integral(unit_square, QS2D3P)(b.norm());   // transport field norm
+	 double b_norm = integral(unit_square, QS2DP2)(b.norm());   // transport field norm
 	 auto tau_k  = 0.5 * delta * h_k / b_norm;   // stabilization parameter
 	 auto a_supg = a + integral(unit_square)(tau_k * (-mu * (dxx(u) + dyy(u)) + dot(b, grad(u)) + u) * dot(b, grad(v)));
 	 auto F_supg = F + integral(unit_square)(f * v + tau_k * f * dot(b, grad(v)));