[examples] Use constraints to setup Dirichlet BC in Stokes example

examples/stokes-taylorhood.cc

@@ -30,6 +30,8 @@
 #include <dune/functions/functionspacebases/lagrangebasis.hh>
 #include <dune/functions/functionspacebases/subspacebasis.hh>
 #include <dune/functions/functionspacebases/boundarydofs.hh>
+#include <dune/fufem/constraints/affineconstraints.hh>
+#include <dune/fufem/constraints/boundaryconstraints.hh>
 
 #include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
 #include <dune/functions/gridfunctions/gridviewfunction.hh>
@@ -200,6 +202,8 @@ int main (int argc, char *argv[]) try
   //  Assemble the system
   /////////////////////////////////////////////////////////
 
+  using namespace Indices;
+
   // { rhs_assembly_begin }
   Vector rhs;
 
@@ -216,7 +220,6 @@ int main (int argc, char *argv[]) try
 
   {
     Dune::Timer timer;
-    using namespace Indices;
     using namespace ::Dune::Fufem::Forms;
     auto velocityBasis = Functions::subspaceBasis(taylorHoodBasis, _0);
     auto pressureBasis = Functions::subspaceBasis(taylorHoodBasis, _1);
@@ -251,25 +254,6 @@ int main (int argc, char *argv[]) try
   // Only velocity components have Dirichlet boundary values
   /////////////////////////////////////////////////////////
 
-  // { initialize_boundary_dofs_vector_begin }
-
-  BitVector isBoundary;
-
-  auto isBoundaryBackend = Dune::Functions::istlVectorBackend(isBoundary);
-  isBoundaryBackend.resize(taylorHoodBasis);
-  isBoundary = false;
-  // { initialize_boundary_dofs_vector_end }
-
-  // { determine_boundary_dofs_begin }
-  using namespace Indices;
-  Functions::forEachBoundaryDOF(
-          Functions::subspaceBasis(taylorHoodBasis, _0),
-          [&] (auto&& index) {
-            isBoundaryBackend[index] = true;
-          });
-  // { determine_boundary_dofs_end }
-
-  // { interpolate_dirichlet_values_begin }
   using Coordinate = GridView::Codim<0> ::Geometry::GlobalCoordinate;
   using VelocityRange = FieldVector<double,dim>;
   auto&& velocityDirichletData = [](Coordinate x)
@@ -277,11 +261,8 @@ int main (int argc, char *argv[]) try
     return VelocityRange{0.0, double(x[0] < 1e-8)};
   };
 
-  Functions::interpolate(
-          Functions::subspaceBasis(taylorHoodBasis, _0), rhs,
-          velocityDirichletData,
-          isBoundary);
-  // { interpolate_dirichlet_values_end }
+  auto constraints = Dune::Fufem::makeAffineConstraints<BitVector, Vector>(taylorHoodBasis);
+  computeBoundaryConstraints(constraints, Functions::subspaceBasis(taylorHoodBasis, _0), velocityDirichletData);
 
   ////////////////////////////////////////////
   //   Modify Dirichlet rows
@@ -289,7 +270,7 @@ int main (int argc, char *argv[]) try
 
   // loop over the matrix rows
   // { set_dirichlet_matrix_begin }
-  incorporateEssentialConstraints(stiffnessMatrix, rhs, isBoundary, rhs);
+  constraints.constrainLinearSystem(stiffnessMatrix, rhs);
   // { set_dirichlet_matrix_end }
 
   ////////////////////////////