]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Written results, modified parameter file to use colorize
authorGiovanni Alzetta <giovannialzetta@hotmail.it>
Wed, 2 May 2018 20:44:46 +0000 (22:44 +0200)
committerMatthias Maier <tamiko@43-1.org>
Tue, 8 May 2018 04:44:55 +0000 (23:44 -0500)
examples/step-60/doc/intro.dox
examples/step-60/doc/results.dox
examples/step-60/step-60.cc

index 221f2b5c5a6e277129cd17864b357b709179cde3..32443edac50fa121d6ace5fcb21f7ec260ff3eb8 100644 (file)
@@ -66,6 +66,10 @@ 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.
 
+As a first example we study the zero Dirichlet boundary condition on
+$\partial\Omega$. The same equations apply if we apply zero Neumann boundary
+conditions on $\partial\Omega$ or a mix of the two.
+
 The variational formulation can be derived by introducing two infinite
 dimensional spaces $V(\Omega)$ and $Q^*(\Gamma)$, respectively for the solution
 $u$ and for the Lagrange multiplier $\lambda$.
@@ -88,7 +92,7 @@ can be taken to be $H^1_0(\Omega)$. The space $Q(\Gamma)$, in the co-dimension
 zero case, should be taken as $H^1(\Gamma)$, while in the co-dimension one case
 should be taken as $H^{1/2}(\Gamma)$.
 
-The function $g$ should, therefore, should be either in $H^1(\Gamma)$ (for the
+The function $g$ should therefore be either in $H^1(\Gamma)$ (for the
 co-dimension zero case) or $H^{1/2}(\Gamma)$ (for the co-dimension one case).
 This leaves us with a Lagrange multiplier $\lambda$ in $Q^*(\Gamma)$, which is
 either $(H^1(\Gamma))^*$ or $H^{-1/2}(\Gamma)$.
@@ -183,17 +187,18 @@ C_{\alpha j} &:=& (v_j, q_\alpha)_\Gamma  &j=1,\dots,n, \alpha = 1,\dots, m \\\\
 G_{\alpha} &:=& (g, q_\alpha)_\Gamma & \alpha = 1,\dots, m.
 @f}
 
-While the matrix $K$ is the standard stiffness matrix for the Poisson problem
-on $\Omega$, and the vector $G$ is a standard right-hand-side vector for a
-finite element problem with forcing term $g$ on $\Gamma$, (see, for example,
-step-3), the matrix $C$ or its transpose $C^T$ are non-standard since they
-couple information on two non-matching grids.
+While the matrix $K$ is the standard stiffness matrix for the Poisson problem on
+$\Omega$, and the vector $G$ is a standard right-hand-side vector for a finite
+element problem with forcing term $g$ on $\Gamma$, (see, for example, step-3),
+the matrix $C$ or its transpose $C^T$ are non-standard since they couple
+information on two non-matching grids.
 
-In particular, the integral that appears in the computation of a single entry of $C$,
-is computed on $\Gamma$. As usual in finite elements, we split this integral on each
-cell of the triangulation used to discretize $\Gamma$, we tranform the integral on $K$ to
-an integral on the reference element $\hat K$, where $F_{K}$ is the corresponding
-shape function, and compute the integral there using a quadrature formula:
+In particular, the integral that appears in the computation of a single entry of
+$C$, is computed on $\Gamma$. As usual in finite elements, we split this
+integral on each cell of the triangulation used to discretize $\Gamma$, we
+tranform the integral on $K$ to an integral on the reference element $\hat K$,
+where $F_{K}$ is the mapping from $\hat K$ to $K$, and compute the integral
+on $\hat K$ using a quadrature formula:
 
 \f[
 C_{\alpha j} := (v_j, q_\alpha)_\Gamma  = \sum_{K\in \Gamma} \int_{\hat K}
@@ -201,19 +206,48 @@ C_{\alpha j} := (v_j, q_\alpha)_\Gamma  = \sum_{K\in \Gamma} \int_{\hat K}
 \sum_{K\in \Gamma} \sum_{i=1}^{n_q}  \big(\hat q_\alpha(\hat x_i)  (v_j \circ F_{K}) (\hat x_i) J_K (\hat x_i) w_i \big)
 \f]
 
-Computing this sum is difficult because we have to evaluate $(v_j \circ F_{K}) (\hat x_i)$; this
-complicated process is illustrated in the following picture:
+Computing this sum is non-trivial because we have to evaluate $(v_j \circ F_{K})
+(\hat x_i)$. In general, if $\Gamma$ and $\Omega$ are not aligned, the point
+$F_{K}(\hat x_i)$ is completely arbitrary with respect to $\Omega$, and unless
+we figure out a way to interpolate all basis functions of $V_h(\Omega)$ on an
+arbitrary point on $\Omega$, we cannot compute the integral needed for an entry
+of the matrix $C$.
+
+The process is described by the following algorithm (and by the picture below):
 
-<p align="center">
-  <img src="https://www.dealii.org/images/steps/developer/step-60_C_interpolation.png" alt="">
-</p>
+- For a given cell $K$ in $\Gamma$ compute the real point $y_i \coloneqq F_{K}
+(\hat x_i)$, where $x_i$ is one of the quadrature points used for the integral
+on $K \subseteq \Gamma$.
 
-This translates in the following algorithm:
-- Fix the cell $K$ in $\Gamma$ and compute the real point $y_i \coloneqq F_{K} (\hat x_i)$
 - Find the cell of $\Omega$ in which $y_i$ lies, say it's $T$.
-- To evaluate the basis function use invert shape function $G_{T}$ on $T$:
- $v_j(y_i) = \hat v_j \circ G^{-1}_{T} (y_i)$.
 
-We solve the above saddle point problem by iterating over the Schur complement
-(which is described, for example, in step-20),
-and we construct such complement using LinearOperator classes.
+- To evaluate the basis function use the inverse of the mapping $G_T$ that
+transforms the reference element $\hat T$ into the element $T$: $v_j(y_i) = \hat
+v_j \circ G^{-1}_{T} (y_i)$.
+
+<p align="center"> <img
+  src="https://www.dealii.org/images/steps/developer/step-60_C_interpolation.png"
+  alt=""> </p>
+
+The three steps above can be computed by calling, in turn:
+
+- GridTools::find_active_cell_around_point, followed by
+
+- Mapping::tranform_real_to_unit_cell
+
+- construct a custom Quadrature formula, containing the point in the reference
+cell
+
+- construct an FEValues object, with the given quadrature formula, and initialized with the cell obtained in the first step.
+
+The deal.II library offers a convenient wrapper to perform the first three steps
+on a collection of points: GridTools::compute_point_locations. If one is
+actually interested in computing the full coupling matrix, then it is possible
+to call the method NonMatching::create_coupling_mass_matrix, that performs the
+above steps in an efficient way, reusing all possible data structures, and
+gathering expensive steps together. This is the function we'll be using later
+in this tutorial. 
+
+We solve the final saddle point problem by an iterative solver, applied to the
+Schur complement $S$ (whose construction is described, for example, in step-20), and
+we construct $S$ using LinearOperator classes.
index c48f2a2d22e72bfd64f883f635617e8a2d20387a..ddaafdff0c8193cd1d4cbf470711de406f8b8068 100644 (file)
 <h1>Results</h1>
 
-<h3> Test case 1: </h3>
-
-For the default problem the value of u on Gamma is 1: this is like imposing a constant
-Dirichlet boundary on the portion of $\Omega$ inside gamma.
-<p align="center">
-  <img src="https://www.dealii.org/images/steps/developer/step-60_1_no_grid.png" alt=""
-   style="width: 50%; height: 50%">
-</p>
-In this second image we can appreciate how, using GridTools::compute_point_locations ,
-we have been able to adaptively refine $\Omega$ exactly where the solution is varying
-the most:
-
-<p align="center">
-  <img src="https://www.dealii.org/images/steps/developer/step-60_1_grid.png" alt=""
-   style="width: 50%; height: 50%">
-</p>
+The first time the code is run, an exception is thrown and nothing is actually
+computed. You'll get the following output:
+@code
+----------------------------------------------------
+Exception on processing: 
+
+--------------------------------------------------------
+An error occurred in line <74> of file <../source/base/parameter_acceptor.cc> in function
+    static void dealii::ParameterAcceptor::initialize(const std::string &, const std::string &, const ParameterHandler::OutputStyle, dealii::ParameterHandler &)
+The violated condition was: 
+    false
+Additional information: 
+    You specified <parameters.prm> as input parameter file, but it does not exist. We created it for you.
+--------------------------------------------------------
+
+Aborting!
+----------------------------------------------------
+@endcode
+
+By inspecting the parameter file, we see the following:
+
+@code
+# Listing of Parameters
+# ---------------------
+subsection Distributed Lagrange<1,2>
+  set Coupling quadrature order                    = 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 embedding space refinement           = 4
+  set Local refinements steps near embedded domain = 3
+  set Use displacement in embedded interface       = false
+  set Verbosity level                              = 10
+
+
+  subsection Embedded configuration
+    # 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  = R=.3, Cx=.4, Cy=.4                 # default: 
+
+    # 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 = R*cos(2*pi*x)+Cx; R*sin(2*pi*x)+Cy # default: 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 Embedded value
+    # 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 = 1     # default: 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
+    set Log result    = true
+    set Max steps     = 1000   # default: 100
+    set Reduction     = 1.e-12 # default: 1.e-2
+    set Tolerance     = 1.e-12 # default: 1.e-10
+  end
+
+end
+@endcode
+
+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-pre
+subsection Distributed Lagrange<1,2>
+  set Coupling quadrature order                    = 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 embedding space refinement           = 4
+  set Local refinements steps near embedded domain = 3
+  set Use displacement in embedded interface       = false
+  set Verbosity level                              = 10
+  subsection Embedded configuration
+    set Function constants  = R=.3, Cx=.4, Cy=.4
+    set Function expression = R*cos(2*pi*x)+Cx; R*sin(2*pi*x)+Cy
+    set Variable names      = x,y,t
+  end
+  subsection Embedded value
+    set Function constants  = 
+    set Function expression = 1
+    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     = 1000
+    set Reduction     = 1.e-12
+    set Tolerance     = 1.e-12
+  end
+end
+@endcode
+
+The rationale behind creating first `parameters.prm` file (the first time the
+program is run) and then a `used_parameters.prm` (every other times you run the
+program), is because you may want to leave most parameters to their default
+values, and only modify a handful of them.
+
+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 embedding space refinement           = 4
+  set Local refinements steps near embedded domain = 3
+  subsection Embedded configuration
+    set Function constants  = R=.3, Cx=.4, Cy=.4
+    set Function expression = R*cos(2*pi*x)+Cx; R*sin(2*pi*x)+Cy
+    set Variable names      = x,y,t
+  end
+  subsection Embedded value
+    set Function constants  =
+    set Function expression = 1
+    set Variable names      = x,y,t
+  end
+end
+@endcode
+
+and you would obtain exactly the same results as in test case 1 below.
+
+<h2> Test case 1: </h2>
+
+For the default problem the value of $u$ on Gamma is set to the constant $1$:
+this is like imposing a constant Dirichlet boundary condition on $\Gamma$, seen
+as boundary of the portion of $\Omega$ inside $\Gamma$. Similarly on $\partial
+\Omega$ we have zero Dirichlet boundary conditions.
+
+<TABLE WIDTH="60%" ALIGN="center">
+  <tr>
+    <td ALIGN="center">
+      <img src="https://www.dealii.org/images/steps/developer/step-60.1_no_grid.png" alt="">
+    </td>
+    <td ALIGN="center">
+      <img src="https://www.dealii.org/images/steps/developer/step-60.1_grid.png" alt="">
+    </td>
+  </tr>
+</table>
+
+The output of the program will look like the following on your pc:
+
+@code
+DEAL::Embedded dofs: 129
+DEAL::Embedding minimal diameter: 0.0110485, embedded maximal diameter: 0.00781250, ratio: 0.707107
+DEAL::Embedding dofs: 2429
+DEAL:cg::Starting value 0.166266
+DEAL:cg::Convergence step 108 value 7.65958e-13
+
+
++---------------------------------------------+------------+------------+
+| Total CPU time elapsed since start          |     0.586s |            |
+|                                             |            |            |
+| 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% |
++---------------------------------+-----------+------------+------------+
+
+
+
++---------------------------------------------+------------+------------+
+| Total wallclock time elapsed since start    |     0.301s |            |
+|                                             |            |            |
+| 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% |
++---------------------------------+-----------+------------+------------+
+
+@endcode
+
+You may notice that, in terms of CPU time, assembling the coupling system is
+twice more expensive than assembling the standard Poisson system, even though
+the matrix is smaller. This is due to the non-matching nature of the
+discretization. Wether 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).
+
+To play around a little bit, we are going to complicate a little the ficticious
+domain as well as the boundary conditions we impose on it.
+
+<h2> Test case 2 and 3: </h2>
+
+If we use the following parameter file:
+@code
+subsection Distributed Lagrange<1,2>
+  set Coupling quadrature order                    = 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            = 8
+  set Initial embedding space refinement           = 4
+  set Local refinements steps near embedded domain = 4
+  set Use displacement in embedded interface       = false
+  set Verbosity level                              = 10
+  subsection Embedded configuration
+    set Function constants  = R=.3, Cx=.5, Cy=.5, r=.1, w=12
+    set Function expression = (R+r*cos(w*pi*x))*cos(2*pi*x)+Cx; (R+r*cos(w*pi*x))*sin(2*pi*x)+Cy
+    set Variable names      = x,y,t
+  end
+  subsection Embedded value
+    set Function constants  = 
+    set Function expression = x-.5
+    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 Reduction     = 1.e-12
+    set Tolerance     = 1.e-12
+  end
+end
+@endcode
+
+we get a "flowery" looking domain, where we impose a linear boundary condition
+$g=x-.5$. This test shows that the method is actually quite accurate in
+recovering an exactly linear function from its boundary conditions, and even
+though the meshes are not aligned, we obtain a pretty good result.
+
+Replacing $x-.5$ with $2(x-.5)^2-2(y-.5)^2$, i.e., modifying the parameter file
+such that we have
+@code
+  ... 
+  subsection Embedded value
+    set Function constants  = 
+    set Function expression = 2*(x-.5)^2-2*(y-.5)^2
+    set Variable names      = x,y,t
+  end
+@endcode
+produce the saddle on the right.
+
+<TABLE WIDTH="60%" ALIGN="center">
+  <tr>
+    <td ALIGN="center">
+      <img src="https://www.dealii.org/images/steps/developer/step-60.3_grid.png" alt="">
+    </td>
+    <td ALIGN="center">
+      <img src="https://www.dealii.org/images/steps/developer/step-60.4_grid.png" alt="">
+    </td>
+  </tr>
+</table>
 
 <a name="extensions"></a>
 <h3>Possibilities for extensions</h3>
 
-Add something
+The `parameters.prm` file is a perfect tool to play with: we can easily impose a
+Neumann boundary, change the parametrization, the function values and the
+refinements.
+
+However, non aligned meshes are complex to handle, as you shall see after some
+changes. For example set the following parameters:
+
+@code
+subsection Distributed Lagrange<1,2>
+  set Homogeneous Dirichlet boundary ids           = 0, 1, 3
+  subsection Embedded configuration
+    set Function constants  = R=.2, Cx=.5, Cy=.4
+    set Function expression = R*sin(2*pi*x)+0.5*R*sin(4*pi*x)+Cx; 1.5*R*cos(2*pi*x)+Cy
+  end
+end
+@endcode
+
+Keeping the default refinement options you shall notice some glitch in the upper part
+of this hot air balloon shape. Why? What's the difference with the following?
+
+<TABLE WIDTH="60%" ALIGN="center">
+  <tr>
+    <td ALIGN="center">
+      <img src="https://www.dealii.org/images/steps/developer/step-60.2_no_grid.png" alt="">
+    </td>
+    <td ALIGN="center">
+      <img src="https://www.dealii.org/images/steps/developer/step-60.2_grid.png" alt="">
+    </td>
+  </tr>
+</table>
+
+To solve this problem we need to improve the way in which adaptive refinement is
+performed. Currently the vector of support points contains the vertices of
+$\Gamma$, if the finite dimensional space for the configuration is of degree
+one: these are too far away to make $\Omega$ refine properly and simply
+increasing the refinements of $\Gamma$ will not be enough, the conjugate
+gradient won't converge.
+
+As you modify the code you shall notice that, for instance, with complex
+geometries there shall be other convergence problems: parameters are handy, but do
+not solve all problems!
 
 <h4> Parallel Code </h4>
 
@@ -33,14 +390,18 @@ parallel::distributed::Triangulation ).
 When using non-matching grids in parallel a problem arises:
 to compute the matrix $C$ a process needs information about both
 meshes on the same portion of real space but, when working with
-distributed meshes, this might not be the case.
+distributed meshes, this information may not be available.
+
 Various strategies can be implemented to tackle this problem:
 
 - distribute the two meshes so that this constraint is satisfied
-- use communication for the parts of real space where the constraint is
- not satisfied
-- make use of a shared triangulation and a distributed triangulation
 
-The latter strategy is clearly the easier to implement, as all
-the function used in this tutorial program can work letting $\Omega$
-be distributed and $\Gamma$ be a shared triangulation.
+- use communication for the parts of real space where the constraint is not
+  satisfied
+
+- use a shared triangulation for the embedding space, and a distributed
+  triangulation for the emdedded configuration
+
+The latter strategy is clearly the easiest to implement, as most of the
+functions used in this tutorial program will work unchanged also in the parallel
+case.
index 2c05c8c39cc77288c924f4593acf368ee7f9fe31..c3c24f016c9095f6999371d7bf4e29f7f6599e6b 100644 (file)
@@ -13,7 +13,7 @@
  *
  * ---------------------------------------------------------------------
  *
- * Author: Luca Heltai, Giovanni Alzetta,
+ * Authors: Luca Heltai, Giovanni Alzetta,
  * International School for Advanced Studies, Trieste, 2018
  */
 
@@ -243,10 +243,12 @@ namespace Step60
       // $\Gamma$.
       unsigned int initial_embedded_refinement                  = 7;
 
-      // A list of boundary ids where we impose homogeneous Dirichlet boundary
+      // 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
-      std::list<types::boundary_id> homogeneous_dirichlet_ids   {0};
+      // 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};
 
       // FiniteElement degree of the embedding space: $V_h(\Omega)$
       unsigned int embedding_space_finite_element_degree        = 1;
@@ -648,7 +650,10 @@ namespace Step60
     // Initializing $\Omega$:
     // constructing the Triangulation and wrapping it into a unique_ptr
     space_grid = std_cxx14::make_unique<Triangulation<spacedim> >();
-    GridGenerator::hyper_cube(*space_grid);
+    // The last argument is set to true: this activates colorize, which
+    // we use to assign the Dirichlet and Neumann conditions. See
+    // GridGenerator::hyper_rectangle for the details on the ids used.
+    GridGenerator::hyper_cube(*space_grid,0,1,true);
     // Requesting the varius values to the parameters object, which is
     // of type DistributedLagrangeProblemParameters
     space_grid->refine_global(parameters.initial_refinement);
@@ -739,6 +744,9 @@ namespace Step60
                                            *embedded_dh,
                                            support_points);
 
+    for (unsigned int i=1; i<embedded_dh->n_dofs(); ++i)
+      support_points.emplace_back((support_points[i-1]+support_points[i])/2);
+
     // Once we have the support points of the embedded finite element space, we
     // would like to identify what cells of the embedding space contain what
     // support point, to get a chance at refining the embedding grid where it is
@@ -958,9 +966,6 @@ namespace Step60
     SparseDirectUMFPACK K_inv_umfpack;
     K_inv_umfpack.initialize(stiffness_matrix);
 
-    // Same thing, for the embedded space
-//    SparseDirectUMFPACK A_inv_umfpack;
-//    A_inv_umfpack.initialize(embedded_stiffness_matrix);
     // Initializing the operators, as described in the introduction
     auto K = linear_operator(stiffness_matrix);
     auto A = linear_operator(embedded_stiffness_matrix);
@@ -968,12 +973,11 @@ namespace Step60
     auto C = transpose_operator(Ct);
 
     auto K_inv = linear_operator(K, K_inv_umfpack);
-//    auto A_inv = linear_operator(A, A_inv_umfpack);
 
-    auto S = C*K_inv*Ct;
     // Using the Schur complement method
+    auto S = C*K_inv*Ct;
     SolverCG<Vector<double> > solver_cg(schur_solver_control);
-    auto S_inv = inverse_operator(S, solver_cg, PreconditionIdentity() );//A_inv);
+    auto S_inv = inverse_operator(S, solver_cg, PreconditionIdentity() );
 
     lambda = S_inv * embedded_rhs;
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.