Allow inhomogeneous Dirichlet bcs,rhs and fix results.dox in step-60

diff --git a/doc/news/changes/minor/20211224MarcoFeder b/doc/news/changes/minor/20211224MarcoFeder
new file mode 100644 (file)
index 0000000..6e0357d
--- /dev/null
+++ b/doc/news/changes/minor/20211224MarcoFeder
@@ -0,0 +1,3 @@
+Changed: The step-60 tutorial now works also with non-zero Dirichlet BCs and non-zero right-hand side
+<br>
+(Marco Feder, 2021/12/24)

diff --git a/examples/step-60/doc/intro.dox b/examples/step-60/doc/intro.dox
index 45a2592e5c8c48b8db5b8f52591df9d13c2d2d92..23d3666bbe692ef2a1075ff02e516550b615c26b 100644 (file)
--- a/examples/step-60/doc/intro.dox
+++ b/examples/step-60/doc/intro.dox
@@ -60,19 +60,20 @@ embedded domain (`dim`) is always smaller by one or equal with respect to the
 dimension of the embedding domain $\Omega$ (`spacedim`).
 
 We are going to solve the following differential problem: given a sufficiently
-regular function $g$ on $\Gamma$, find the solution $u$ to
+regular function $g$ on $\Gamma$, a forcing term $f \in L^2(\Omega)$ and a Dirichlet boundary condition
+$u_D$ on $\partial \Omega$, find the solution $u$ to
 
 @f{eqnarray*}{
-- \Delta u + \gamma^T \lambda &=& 0  \text{ in } \Omega\\
+- \Delta u + \gamma^T \lambda &=& f  \text{ in } \Omega\\
 \gamma u &=& g  \text{ in } \Gamma \\
-u & = & 0 \text{ on } \partial\Omega.
+u & = & u_D \text{ on } \partial\Omega.
 @f}
 
-This is a constrained problem, where we are looking for a harmonic function $u$
-that satisfies homogeneous boundary conditions on $\partial\Omega$, subject to
-the constraint $\gamma u = g$ using a Lagrange multiplier.
+This is a constrained problem, where we are looking for a function $u$ that solves the
+Poisson equation and that satisfies Dirichlet boundary conditions $u=u_D$ on $\partial \Omega$,
+subject to the constraint $\gamma u = g$ using a Lagrange multiplier.
 
-This problem has a physical interpretation: harmonic functions, i.e., functions
+When $f=0$ this problem has a physical interpretation: harmonic functions, i.e., functions
 that satisfy the Laplace equation, can be thought of as the displacements of a
 membrane whose boundary values are prescribed. The current situation then
 corresponds to finding the shape of a membrane for which not only the
@@ -98,7 +99,7 @@ conditions on $\partial\Omega$, we obtain the following variational problem:
 
 Given a sufficiently regular function $g$ on $\Gamma$, find the solution $u$ to
 @f{eqnarray*}{
-(\nabla u, \nabla v)_{\Omega} + (\lambda, \gamma v)_{\Gamma} &=& 0 \qquad \forall v \in V(\Omega) \\
+(\nabla u, \nabla v)_{\Omega} + (\lambda, \gamma v)_{\Gamma} &=& (f,v)_{\Omega} \qquad \forall v \in V(\Omega) \\
 (\gamma u, q)_{\Gamma} &=& (g,q)_{\Gamma} \qquad \forall q \in Q(\Gamma),
 @f}
 
@@ -200,7 +201,7 @@ u \\
 \end{pmatrix}
 =
 \begin{pmatrix}
-0 \\
+F \\
 G
 \end{pmatrix}
 \f]
@@ -210,6 +211,7 @@ where
 @f{eqnarray*}{
 K_{ij} &\dealcoloneq& (\nabla v_j, \nabla v_i)_\Omega   \qquad i,j=1,\dots,n \\
 C_{\alpha j} &\dealcoloneq& (v_j, q_\alpha)_\Gamma  \qquad j=1,\dots,n, \alpha = 1,\dots, m \\\\
+F_{i} &\dealcoloneq& (f, v_i)_\Omega   \qquad i=1,\dots,n \\
 G_{\alpha} &\dealcoloneq& (g, q_\alpha)_\Gamma \qquad \alpha = 1,\dots, m.
 @f}
 

diff --git a/examples/step-60/doc/results.dox b/examples/step-60/doc/results.dox
index ae9c27c51c8ff270d3d7786cd74f675cd1808329..1b4edef822bb9826afb93865b2925d36aa854242 100644 (file)
--- a/examples/step-60/doc/results.dox
+++ b/examples/step-60/doc/results.dox
@@ -32,11 +32,11 @@ parameter file, we see the following:
 # ---------------------
 subsection Distributed Lagrange<1,2>
   set Coupling quadrature order                    = 3
+  set Dirichlet boundary ids                       = 0, 1, 2, 3
   set Embedded configuration finite element degree = 1
   set Embedded space finite element degree         = 1
   set Embedding space finite element degree        = 1
-  set Homogeneous Dirichlet boundary ids           = 0, 1, 2, 3
-  set Initial embedded space refinement            = 7
+  set Initial embedded space refinement            = 8
   set Initial embedding space refinement           = 4
   set Local refinements steps near embedded domain = 3
   set Use displacement in embedded interface       = false
@@ -125,6 +125,88 @@ subsection Distributed Lagrange<1,2>
     set Variable names      = x,y,t
   end
 
+  subsection Embedding Dirichlet boundary conditions
+    # Sometimes it is convenient to use symbolic constants in the expression
+    # that describes the function, rather than having to use its numeric value
+    # everywhere the constant appears. These values can be defined using this
+    # parameter, in the form `var1=value1, var2=value2, ...'.
+    #
+    # A typical example would be to set this runtime parameter to
+    # `pi=3.1415926536' and then use `pi' in the expression of the actual
+    # formula. (That said, for convenience this class actually defines both
+    # `pi' and `Pi' by default, but you get the idea.)
+    set Function constants  =
+
+    # The formula that denotes the function you want to evaluate for
+    # particular values of the independent variables. This expression may
+    # contain any of the usual operations such as addition or multiplication,
+    # as well as all of the common functions such as `sin' or `cos'. In
+    # addition, it may contain expressions like `if(x>0, 1, -1)' where the
+    # expression evaluates to the second argument if the first argument is
+    # true, and to the third argument otherwise. For a full overview of
+    # possible expressions accepted see the documentation of the muparser
+    # library at http://muparser.beltoforion.de/.
+    #
+    # If the function you are describing represents a vector-valued function
+    # with multiple components, then separate the expressions for individual
+    # components by a semicolon.
+    set Function expression = 0
+
+    # The names of the variables as they will be used in the function,
+    # separated by commas. By default, the names of variables at which the
+    # function will be evaluated are `x' (in 1d), `x,y' (in 2d) or `x,y,z' (in
+    # 3d) for spatial coordinates and `t' for time. You can then use these
+    # variable names in your function expression and they will be replaced by
+    # the values of these variables at which the function is currently
+    # evaluated. However, you can also choose a different set of names for the
+    # independent variables at which to evaluate your function expression. For
+    # example, if you work in spherical coordinates, you may wish to set this
+    # input parameter to `r,phi,theta,t' and then use these variable names in
+    # your function expression.
+    set Variable names      = x,y,t
+  end
+
+  subsection Embedding rhs function
+    # Sometimes it is convenient to use symbolic constants in the expression
+    # that describes the function, rather than having to use its numeric value
+    # everywhere the constant appears. These values can be defined using this
+    # parameter, in the form `var1=value1, var2=value2, ...'.
+    #
+    # A typical example would be to set this runtime parameter to
+    # `pi=3.1415926536' and then use `pi' in the expression of the actual
+    # formula. (That said, for convenience this class actually defines both
+    # `pi' and `Pi' by default, but you get the idea.)
+    set Function constants  =
+
+    # The formula that denotes the function you want to evaluate for
+    # particular values of the independent variables. This expression may
+    # contain any of the usual operations such as addition or multiplication,
+    # as well as all of the common functions such as `sin' or `cos'. In
+    # addition, it may contain expressions like `if(x>0, 1, -1)' where the
+    # expression evaluates to the second argument if the first argument is
+    # true, and to the third argument otherwise. For a full overview of
+    # possible expressions accepted see the documentation of the muparser
+    # library at http://muparser.beltoforion.de/.
+    #
+    # If the function you are describing represents a vector-valued function
+    # with multiple components, then separate the expressions for individual
+    # components by a semicolon.
+    set Function expression = 0
+
+    # The names of the variables as they will be used in the function,
+    # separated by commas. By default, the names of variables at which the
+    # function will be evaluated are `x' (in 1d), `x,y' (in 2d) or `x,y,z' (in
+    # 3d) for spatial coordinates and `t' for time. You can then use these
+    # variable names in your function expression and they will be replaced by
+    # the values of these variables at which the function is currently
+    # evaluated. However, you can also choose a different set of names for the
+    # independent variables at which to evaluate your function expression. For
+    # example, if you work in spherical coordinates, you may wish to set this
+    # input parameter to `r,phi,theta,t' and then use these variable names in
+    # your function expression.
+    set Variable names      = x,y,t
+  end
+
   subsection Schur solver control
     set Log frequency = 1
     set Log history   = false
@@ -141,15 +223,13 @@ If you now run the program, you will get a file called `used_parameters.prm`,
 containing a shorter version of the above parameters (without comments and
 documentation), documenting all parameters that were used to run your program:
 @code
-# Parameter file generated with
-# DEAL_II_PACKAGE_VERSION = 9.0.0
 subsection Distributed Lagrange<1,2>
   set Coupling quadrature order                    = 3
+  set Dirichlet boundary ids                       = 0, 1, 2, 3
   set Embedded configuration finite element degree = 1
   set Embedded space finite element degree         = 1
   set Embedding space finite element degree        = 1
-  set Homogeneous Dirichlet boundary ids           = 0, 1, 2, 3
-  set Initial embedded space refinement            = 7
+  set Initial embedded space refinement            = 8
   set Initial embedding space refinement           = 4
   set Local refinements steps near embedded domain = 3
   set Use displacement in embedded interface       = false
@@ -164,6 +244,16 @@ subsection Distributed Lagrange<1,2>
     set Function expression = 1
     set Variable names      = x,y,t
   end
+  subsection Embedding Dirichlet boundary conditions
+    set Function constants  =
+    set Function expression = 0
+    set Variable names      = x,y,t
+  end
+  subsection Embedding rhs function
+    set Function constants  =
+    set Function expression = 0
+    set Variable names      = x,y,t
+  end
   subsection Schur solver control
     set Log frequency = 1
     set Log history   = false
@@ -184,7 +274,7 @@ For example, you could use the following (perfectly valid) parameter file with
 this tutorial program:
 @code
 subsection Distributed Lagrange<1,2>
-  set Initial embedded space refinement            = 7
+  set Initial embedded space refinement            = 8
   set Initial embedding space refinement           = 4
   set Local refinements steps near embedded domain = 3
   subsection Embedded configuration
@@ -197,6 +287,16 @@ subsection Distributed Lagrange<1,2>
     set Function expression = 1
     set Variable names      = x,y,t
   end
+  subsection Embedding Dirichlet boundary conditions
+    set Function constants  =
+    set Function expression = 0
+    set Variable names      = x,y,t
+  end
+  subsection Embedding rhs function
+    set Function constants  =
+    set Function expression = 0
+    set Variable names      = x,y,t
+  end
 end
 @endcode
 
@@ -230,55 +330,50 @@ as boundary of the portion of $\Omega$ inside $\Gamma$. Similarly on $\partial
 The output of the program will look like the following:
 
 @code
-DEAL::Embedded dofs: 129
-DEAL::Embedding minimal diameter: 0.0110485, embedded maximal diameter: 0.00781250, ratio: 0.707107
+DEAL::Embedded dofs: 257
+DEAL::Embedding minimal diameter: 0.0110485, embedded maximal diameter: 0.00736292, ratio: 0.666416
 DEAL::Embedding dofs: 2429
-DEAL:cg::Starting value 0.166266
-DEAL:cg::Convergence step 108 value 7.65958e-13
+DEAL:cg::Starting value 0.117692
+DEAL:cg::Convergence step 594 value 8.06558e-13
 
 
 +---------------------------------------------+------------+------------+
-| Total CPU time elapsed since start          |     0.586s |            |
+| Total CPU time elapsed since start          |      1.48s |            |
 |                                             |            |            |
 | Section                         | no. calls |  CPU time  | % of total |
 +---------------------------------+-----------+------------+------------+
-| Assemble coupling system        |         1 |     0.132s |        23% |
-| Assemble system                 |         1 |    0.0733s |        12% |
-| Output results                  |         1 |     0.087s |        15% |
-| Setup coupling                  |         1 |    0.0244s |       4.2% |
-| Setup grids and dofs            |         1 |    0.0907s |        15% |
-| Solve system                    |         1 |     0.178s |        30% |
+| Assemble coupling system        |         1 |    0.0381s |       2.6% |
+| Assemble system                 |         1 |     0.153s |        10% |
+| Output results                  |         1 |      0.11s |       7.4% |
+| Setup coupling                  |         1 |    0.0364s |       2.5% |
+| Setup grids and dofs            |         1 |     0.168s |        11% |
+| Solve system                    |         1 |     0.974s |        66% |
 +---------------------------------+-----------+------------+------------+
 
 
 
 +---------------------------------------------+------------+------------+
-| Total wallclock time elapsed since start    |     0.301s |            |
+| Total wallclock time elapsed since start    |     0.798s |            |
 |                                             |            |            |
 | Section                         | no. calls |  wall time | % of total |
 +---------------------------------+-----------+------------+------------+
-| Assemble coupling system        |         1 |    0.0385s |        13% |
-| Assemble system                 |         1 |    0.0131s |       4.3% |
-| Output results                  |         1 |    0.0736s |        24% |
-| Setup coupling                  |         1 |    0.0234s |       7.7% |
-| Setup grids and dofs            |         1 |    0.0679s |        23% |
-| Solve system                    |         1 |    0.0832s |        28% |
+| Assemble coupling system        |         1 |    0.0469s |       5.9% |
+| Assemble system                 |         1 |    0.0348s |       4.4% |
+| Output results                  |         1 |    0.0821s |        10% |
+| Setup coupling                  |         1 |    0.0371s |       4.7% |
+| Setup grids and dofs            |         1 |     0.157s |        20% |
+| Solve system                    |         1 |     0.436s |        55% |
 +---------------------------------+-----------+------------+------------+
 
 @endcode
 
-You may notice that, in terms of CPU time, assembling the coupling system is
-twice as expensive as assembling the standard Poisson system, even though the
-matrix is smaller. This is due to the non-matching nature of the discretization.
-Whether this is acceptable or not, depends on the applications.
-
 If the problem was set in a three-dimensional setting, and the immersed mesh was
 time dependent, it would be much more expensive to recreate the mesh at each
 step rather than use the technique we present here. Moreover, you may be able to
 create a very fast and optimized solver on a uniformly refined square or cubic
 grid, and embed the domain where you want to perform your computation using the
 technique presented here. This would require you to only have a surface
-representatio of your domain (a much cheaper and easier mesh to produce).
+representation of your domain (a much cheaper and easier mesh to produce).
 
 To play around a little bit, we are going to complicate a little the fictitious
 domain as well as the boundary conditions we impose on it.
@@ -308,11 +403,21 @@ subsection Distributed Lagrange<1,2>
     set Function expression = x-.5
     set Variable names      = x,y,t
   end
+  subsection Embedding Dirichlet boundary conditions
+    set Function constants  =
+    set Function expression = 0
+    set Variable names      = x,y,t
+  end
+  subsection Embedding rhs function
+    set Function constants  =
+    set Function expression = 0.0
+    set Variable names      = x,y,t
+  end
   subsection Schur solver control
     set Log frequency = 1
     set Log history   = false
     set Log result    = true
-    set Max steps     = 100000
+    set Max steps     = 1000
     set Reduction     = 1.e-12
     set Tolerance     = 1.e-12
   end

diff --git a/examples/step-60/step-60.cc b/examples/step-60/step-60.cc
index 162f876f5cdca3a06ef502ce467e4092341f3d8a..7e286064486e476f65c5fe9bff0b2fea735f3952 100644 (file)
--- a/examples/step-60/step-60.cc
+++ b/examples/step-60/step-60.cc
@@ -247,12 +247,12 @@ namespace Step60
       // $\Gamma$.
       unsigned int initial_embedded_refinement = 8;
 
-      // The list of boundary ids where we impose homogeneous Dirichlet boundary
-      // conditions. On the remaining boundary ids (if any), we impose
-      // homogeneous Neumann boundary conditions.
-      // As a default problem we have zero Dirichlet boundary conditions on
+      // The list of boundary ids where we impose (possibly inhomogeneous)
+      // Dirichlet boundary conditions. On the remaining boundary ids (if any),
+      // we impose homogeneous Neumann boundary conditions. As a default problem
+      // we have zero Dirichlet boundary conditions on
       // $\partial \Omega$
-      std::list<types::boundary_id> homogeneous_dirichlet_ids{0, 1, 2, 3};
+      std::list<types::boundary_id> dirichlet_ids{0, 1, 2, 3};
 
       // FiniteElement degree of the embedding space: $V_h(\Omega)$
       unsigned int embedding_space_finite_element_degree = 1;
@@ -364,12 +364,22 @@ namespace Step60
 
     std::unique_ptr<Mapping<dim, spacedim>> embedded_mapping;
 
-    // We do the same thing to specify the value of the function $g$,
-    // which is what we want our solution to be in the embedded space.
+    // We do the same thing to specify the value of the forcing term $f$.
     // In this case the Function is a scalar one.
+    ParameterAcceptorProxy<Functions::ParsedFunction<spacedim>>
+      embedding_rhs_function;
+
+    // Then, we specify the value of the function $g$, which is what we want
+    // our solution to be in the embedded space.
+    // Again, in this case the Function is a scalar one.
     ParameterAcceptorProxy<Functions::ParsedFunction<spacedim>>
       embedded_value_function;
 
+    // Finally, the value of the Dirichlet boundary conditions on $\partial
+    // \Omega$ is specified.
+    ParameterAcceptorProxy<Functions::ParsedFunction<spacedim>>
+      embedding_dirichlet_boundary_function;
+
     // Similarly to what we have done with the Functions::ParsedFunction class,
     // we repeat the same for the ReductionControl class, allowing us to
     // specify all possible stopping criteria for the Schur complement
@@ -387,7 +397,7 @@ namespace Step60
     AffineConstraints<double> constraints;
 
     Vector<double> solution;
-    Vector<double> rhs;
+    Vector<double> embedding_rhs;
 
     Vector<double> lambda;
     Vector<double> embedded_rhs;
@@ -526,8 +536,7 @@ namespace Step60
     add_parameter("Local refinements steps near embedded domain",
                   delta_refinement);
 
-    add_parameter("Homogeneous Dirichlet boundary ids",
-                  homogeneous_dirichlet_ids);
+    add_parameter("Dirichlet boundary ids", dirichlet_ids);
 
     add_parameter("Use displacement in embedded interface", use_displacement);
 
@@ -556,7 +565,10 @@ namespace Step60
     const Parameters &parameters)
     : parameters(parameters)
     , embedded_configuration_function("Embedded configuration", spacedim)
+    , embedding_rhs_function("Embedding rhs function")
     , embedded_value_function("Embedded value")
+    , embedding_dirichlet_boundary_function(
+        "Embedding Dirichlet boundary conditions")
     , schur_solver_control("Schur solver control")
     , monitor(std::cout, TimerOutput::summary, TimerOutput::cpu_and_wall_times)
   {
@@ -581,9 +593,15 @@ namespace Step60
                                    "R*cos(2*pi*x)+Cx; R*sin(2*pi*x)+Cy");
       });
 
+    embedding_rhs_function.declare_parameters_call_back.connect(
+      []() -> void { ParameterAcceptor::prm.set("Function expression", "0"); });
+
     embedded_value_function.declare_parameters_call_back.connect(
       []() -> void { ParameterAcceptor::prm.set("Function expression", "1"); });
 
+    embedding_dirichlet_boundary_function.declare_parameters_call_back.connect(
+      []() -> void { ParameterAcceptor::prm.set("Function expression", "0"); });
+
     schur_solver_control.declare_parameters_call_back.connect([]() -> void {
       ParameterAcceptor::prm.set("Max steps", "1000");
       ParameterAcceptor::prm.set("Reduction", "1.e-12");
@@ -828,8 +846,7 @@ namespace Step60
 
   // We now set up the DoFs of $\Omega$ and $\Gamma$: since they are
   // fundamentally independent (except for the fact that $\Omega$'s mesh is more
-  // refined "around"
-  // $\Gamma$) the procedure is standard.
+  // refined "around" $\Gamma$) the procedure is standard.
   template <int dim, int spacedim>
   void DistributedLagrangeProblem<dim, spacedim>::setup_embedding_dofs()
   {
@@ -839,10 +856,10 @@ namespace Step60
     space_dh->distribute_dofs(*space_fe);
 
     DoFTools::make_hanging_node_constraints(*space_dh, constraints);
-    for (auto id : parameters.homogeneous_dirichlet_ids)
+    for (auto id : parameters.dirichlet_ids)
       {
         VectorTools::interpolate_boundary_values(
-          *space_dh, id, Functions::ZeroFunction<spacedim>(), constraints);
+          *space_dh, id, embedding_dirichlet_boundary_function, constraints);
       }
     constraints.close();
 
@@ -852,7 +869,7 @@ namespace Step60
     stiffness_sparsity.copy_from(dsp);
     stiffness_matrix.reinit(stiffness_sparsity);
     solution.reinit(space_dh->n_dofs());
-    rhs.reinit(space_dh->n_dofs());
+    embedding_rhs.reinit(space_dh->n_dofs());
 
     deallog << "Embedding dofs: " << space_dh->n_dofs() << std::endl;
   }
@@ -914,11 +931,13 @@ namespace Step60
     {
       TimerOutput::Scope timer_section(monitor, "Assemble system");
 
-      // Embedding stiffness matrix $K$, and the right hand side $G$.
+      // Embedding stiffness matrix $K$ and right hand sides $F$ and $G$.
       MatrixTools::create_laplace_matrix(
         *space_dh,
         QGauss<spacedim>(2 * space_fe->degree + 1),
         stiffness_matrix,
+        embedding_rhs_function,
+        embedding_rhs,
         static_cast<const Function<spacedim> *>(nullptr),
         constraints);
 
@@ -978,9 +997,9 @@ namespace Step60
     SolverCG<Vector<double>> solver_cg(schur_solver_control);
     auto S_inv = inverse_operator(S, solver_cg, PreconditionIdentity());
 
-    lambda = S_inv * embedded_rhs;
+    lambda = S_inv * (C * K_inv * embedding_rhs - embedded_rhs);
 
-    solution = K_inv * Ct * lambda;
+    solution = K_inv * (embedding_rhs - Ct * lambda);
 
     constraints.distribute(solution);
   }