Merge pull request #45 from luca-heltai/fix-warnings

Fix warnings

include/parsed_tools/grid_generator.h

@@ -16,6 +16,7 @@
 #ifndef parsed_tools_grid_generator_h
 #define parsed_tools_grid_generator_h
 
+#include <deal.II/base/logstream.h>
 #include <deal.II/base/parameter_acceptor.h>
 
 #include <deal.II/distributed/grid_refinement.h>

notebooks/manufactured_solution.ipynb

@@ -0,0 +1,132 @@
+{
+ "cells": [
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Manufactured solutions for Stokes problem\n",
+    "\n",
+    "When solving numerically the stokes problem \n",
+    "\n",
+    "$$\n",
+    "\\begin{aligned}\n",
+    "-\\Delta \\boldsymbol{u} + \\nabla p & = f \\\\\n",
+    "\\nabla \\cdot \\boldsymbol{u} &= 0,\n",
+    "\\end{aligned}\n",
+    "$$\n",
+    "\n",
+    "we can construct manufactured solutions that satisfy the divergence free\n",
+    "constraint, by taking a solution $\\boldsymbol{u}$ that is the curl of a scalar\n",
+    "field (in two-dimensions) or the curl of a vector field in three dimensions.\n",
+    "\n",
+    "For two dimensions: \n",
+    "\n",
+    "$$\n",
+    "\\begin{split}\n",
+    "\\boldsymbol{u}_x & = \\partial_y g \\\\\n",
+    "\\boldsymbol{u}_y & = -\\partial_x g\n",
+    "\\end{split}\n",
+    "$$\n",
+    "\n",
+    "which satisfy by construction $\\nabla \\cdot \\boldsymbol{u} = 0$ for $g\\in\n",
+    "C^2(\\Omega)$.\n",
+    "\n",
+    "We can then build a forcing term $f$ that would force the system to have this\n",
+    "exact solution, by using the equality $f = -\\Delta \\boldsymbol{u} + \\nabla p$.\n",
+    "\n",
+    "If we want to use these in input files for the `stokes` executable of the\n",
+    "`fsi-suite`, we must replace the `**` with `^`, and we must write them in a\n",
+    "format that deal.II understands. This is done in the following cells, where we\n",
+    "first compute symbolically the differential operators, and then produce an\n",
+    "output that is readable by deal.II."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from sympy import *\n",
+    "import textwrap"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "def write_parameters(precursor, p):\n",
+    "    # ux, uy, and p are the exact solution of the velocity and pressure\n",
+    "    ux = diff(precursor, y)\n",
+    "    uy = -diff(precursor, x)\n",
+    "\n",
+    "    # lap_u_x and lap_u_y are minus the laplacian of ux and uy\n",
+    "    lap_u_x = - diff(ux, x, 2) - diff(ux, y, 2)\n",
+    "    lap_u_y = - diff(uy, x, 2) - diff(uy, y, 2)\n",
+    "\n",
+    "    # fx and fy are the resulting forcing terms\n",
+    "    fx = lap_u_x + diff(p, x)\n",
+    "    fy = lap_u_y + diff(p, y)\n",
+    "\n",
+    "    def prm(ux, uy, p, fx, fy):\n",
+    "        def to_prm(x):\n",
+    "            text = str(x).replace('**', '^')\n",
+    "            start_len = len(\"  set Forcing term   =  \")\n",
+    "            return ' \\\\\\n'.join(textwrap.wrap(text, 80-start_len, break_long_words=False, \n",
+    "                                             subsequent_indent=' '*start_len))\n",
+    "\n",
+    "        print(\n",
+    "            \"subsection Functions\\n  set Exact solution = \", to_prm(ux), \";\",\n",
+    "            to_prm(uy), \";\", to_prm(p), \"\\n\",\n",
+    "            \" set Forcing term   = \", to_prm(fx), \";\",\n",
+    "            to_prm(fy), \"; 0\", \"\\nend\\n\"\n",
+    "            )\n",
+    "        \n",
+    "    prm(ux, uy, p, fx, fy)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Spatial variables\n",
+    "x,y = symbols('x y')\n",
+    "\n",
+    "# Precursor of th exact solution. We take the curl of this function to have a \n",
+    "# divergence free velocity field\n",
+    "precursor = sin(pi*x)**2*sin(pi*y)**2\n",
+    "pressure = cos(2*pi*x)*cos(2*pi*y)\n",
+    "\n",
+    "write_parameters(precursor, pressure)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.15"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

source/data_out.cc

@@ -177,7 +177,7 @@ namespace ParsedTools
           {
             static Vector<float> material_ids;
             material_ids.reinit(dh.get_triangulation().n_active_cells());
-            for (const auto cell : dh.active_cell_iterators())
+            for (const auto &cell : dh.active_cell_iterators())
               material_ids(cell->active_cell_index()) = cell->material_id();
             data_out->add_data_vector(
               material_ids,

source/grid_generator.cc

@@ -275,4 +275,4 @@ namespace ParsedTools
   template class GridGenerator<2, 3>;
   template class GridGenerator<3, 3>;
 
-} // namespace ParsedTools
+} // namespace ParsedTools

source/pdes/linear_visco_elasticity.cc

@@ -44,7 +44,7 @@ namespace PDEs
                    "Second Lame coefficient",
                    "Shear viscosity",
                    "Bulk viscosity"})
-    , material_ids_0({{0}})
+    , material_ids_0({0})
     , eulerian_mapping(this->dof_handler, "/LinearViscoElasticity/Mapping")
   {
     this->add_parameter("Material ids of region 0", material_ids_0);

source/pdes/serial/stokes.cc

@@ -52,15 +52,22 @@ namespace PDEs
           ParsedTools::Components::join(std::vector<std::string>(dim + 1, "0"),
                                         ";"),
           "Exact solution")
-      , boundary_conditions("/Stokes/Boundary conditions", component_names)
+      , boundary_conditions(
+          "/Stokes/Boundary conditions",
+          component_names,
+          {{numbers::internal_face_boundary_id}},
+          {"u"},
+          {{ParsedTools::BoundaryConditionType::dirichlet}},
+          {ParsedTools::Components::join(std::vector<std::string>(dim + 1, "0"),
+                                         ";")})
       , error_table(Utilities::split_string_list(component_names),
                     {{VectorTools::H1_norm, VectorTools::L2_norm},
                      {VectorTools::L2_norm}})
       , data_out("/Stokes/Output")
       , velocity(0)
       , pressure(dim)
     {
-      enter_subsection("Error table");
+      enter_subsection("Error table (" + component_names + ")");
       enter_my_subsection(this->prm);
       error_table.add_parameters(this->prm);
       leave_my_subsection(this->prm);