diff --git a/partitioned-heat-conduction-complex/README.md b/partitioned-heat-conduction-complex/README.md
new file mode 100644
index 000000000..fef453fd4
--- /dev/null
+++ b/partitioned-heat-conduction-complex/README.md
@@ -0,0 +1,27 @@
+---
+title: Partitioned heat conduction (complex setup)
+permalink: tutorials-partitioned-heat-conduction-complex.html
+keywords: FEniCS, Heat conduction
+summary: This tutorial is the advanced version of the `partitioned-heat-conduction` tutorial. More advanced features and geometries may be used.
+---
+
+{% include important.html content="We have not yet ported the documentation of the preCICE tutorials from the preCICE wiki to here. Please go to the [preCICE wiki](https://github.com/precice/precice/wiki#2-getting-started---tutorials)" %}
+
+## Setup
+
+This case is an advanced version of `partitioned-heat-conduction`. Some advanced features offered by this case:
+
+* Geometries may be chosen arbitrarily. One possibility is to use a circle and a rectangular plate with a hole, but you can also provide your own geometry, if you want.
+* You may combine arbitrary mesh resolutions at the coupling interface.
+* Nearest projection mapping is used.
+* The Dirichlet and Neumann participants may be swapped arbitrarily.
+* The exchanged temperature is still scalar valued, but the heat flux is vector valued.
+* You can decide to use a time dependent heat flux and right-hand side to make the problem more challenging.
+
+## Available solvers and dependencies
+
+See `partitioned-heat-conduction`, only `fenics` is provided as a solver.
+
+## Running the simulation
+
+See `partitioned-heat-conduction`. The additional featured mentioned above can be activated via command line arguments. Please run `python3 fenics/heat.py --help` for a full list of provided arguments.
diff --git a/partitioned-heat-conduction/fenics-fenics/.gitignore b/partitioned-heat-conduction-complex/fenics/.gitignore
similarity index 100%
rename from partitioned-heat-conduction/fenics-fenics/.gitignore
rename to partitioned-heat-conduction-complex/fenics/.gitignore
diff --git a/partitioned-heat-conduction/fenics-fenics/Allclean b/partitioned-heat-conduction-complex/fenics/Allclean
similarity index 100%
rename from partitioned-heat-conduction/fenics-fenics/Allclean
rename to partitioned-heat-conduction-complex/fenics/Allclean
diff --git a/partitioned-heat-conduction/fenics-fenics/errorcomputation.py b/partitioned-heat-conduction-complex/fenics/errorcomputation.py
similarity index 100%
rename from partitioned-heat-conduction/fenics-fenics/errorcomputation.py
rename to partitioned-heat-conduction-complex/fenics/errorcomputation.py
diff --git a/partitioned-heat-conduction/fenics-fenics/heat.py b/partitioned-heat-conduction-complex/fenics/heat.py
similarity index 76%
rename from partitioned-heat-conduction/fenics-fenics/heat.py
rename to partitioned-heat-conduction-complex/fenics/heat.py
index e528dd9d9..5af963950 100644
--- a/partitioned-heat-conduction/fenics-fenics/heat.py
+++ b/partitioned-heat-conduction-complex/fenics/heat.py
@@ -29,7 +29,7 @@
File, solve, lhs, rhs, grad, inner, dot, dx, ds, interpolate, VectorFunctionSpace, MeshFunction, MPI
from fenicsprecice import Adapter
from errorcomputation import compute_errors
-from my_enums import ProblemType, Subcycling, DomainPart
+from my_enums import ProblemType, DomainPart
import argparse
import numpy as np
from problem_setup import get_geometry, get_problem_setup
@@ -53,9 +53,12 @@ def determine_gradient(V_g, u, flux):
solve(a == L, flux)
-parser = argparse.ArgumentParser(description='Solving heat equation for simple or complex interface case')
-parser.add_argument("-d", "--dirichlet", help="create a dirichlet problem", dest='dirichlet', action='store_true')
-parser.add_argument("-n", "--neumann", help="create a neumann problem", dest='neumann', action='store_true')
+parser = argparse.ArgumentParser(
+ description='Solving heat equation for simple or complex interface case')
+parser.add_argument("-d", "--dirichlet", help="create a dirichlet problem",
+ dest='dirichlet', action='store_true')
+parser.add_argument("-n", "--neumann", help="create a neumann problem",
+ dest='neumann', action='store_true')
parser.add_argument("-g", "--gamma", help="parameter gamma to set temporal dependence of heat flux", default=0.0,
type=float)
parser.add_argument("-a", "--arbitrary-coupling-interface",
@@ -69,32 +72,9 @@ def determine_gradient(V_g, u, flux):
args = parser.parse_args()
-subcycle = Subcycling.NONE
-
-fenics_dt = None
-error_tol = None
-
-# for all scenarios, we assume precice_dt == .1
-if subcycle is Subcycling.NONE and not args.arbitrary_coupling_interface:
- fenics_dt = .1 # time step size
- error_tol = 10 ** -6 # Error is bounded by coupling accuracy. In theory we would obtain the analytical solution.
- print("Subcycling = NO, Arbitrary coupling interface = NO, error tolerance = {}".format(error_tol))
-elif subcycle is Subcycling.NONE and args.arbitrary_coupling_interface:
- fenics_dt = .1 # time step size
- error_tol = 10 ** -3 # error low, if we do not subcycle. In theory we would obtain the analytical solution.
- # TODO For reasons, why we currently still have a relatively high error, see milestone https://github.com/precice/fenics-adapter/milestone/1
- print("Subcycling = NO, Arbitrary coupling interface = YES, error tolerance = {}".format(error_tol))
-elif subcycle is Subcycling.MATCHING:
- fenics_dt = .01 # time step size
- error_tol = 10 ** -2 # error increases. If we use subcycling, we cannot assume that we still get the exact solution.
- # TODO Using waveform relaxation, we should be able to obtain the exact solution here, as well.
- print("Subcycling = YES, Matching. error tolerance = {}".format(error_tol))
-elif subcycle is Subcycling.NONMATCHING:
- fenics_dt = .03 # time step size
- error_tol = 10 ** -1 # error increases. If we use subcycling, we cannot assume that we still get the exact solution.
- # TODO Using waveform relaxation, we should be able to obtain the exact solution here, as well.
- print("Subcycling = YES, Non-matching. error tolerance = {}".format(error_tol))
-
+fenics_dt = .1
+# Error is bounded by coupling accuracy. In theory we can obtain the analytical solution.
+error_tol = 10 ** -6
alpha = 3 # parameter alpha
beta = 1.3 # parameter beta
gamma = args.gamma # parameter gamma, dependence of heat flux on time
@@ -112,7 +92,8 @@ def determine_gradient(V_g, u, flux):
beta=beta, gamma=gamma, t=0)
u_D_function = interpolate(u_D, V)
# Define flux in x direction on coupling interface (grad(u_D) in normal direction)
-f_N = Expression(("2 * gamma*t*x[0] + 2 * (1-gamma)*x[0]", "2 * alpha*x[1]"), degree=1, gamma=gamma, alpha=alpha, t=0)
+f_N = Expression(("2 * gamma*t*x[0] + 2 * (1-gamma)*x[0]",
+ "2 * alpha*x[1]"), degree=1, gamma=gamma, alpha=alpha, t=0)
f_N_function = interpolate(f_N, V_g)
# Define initial value
@@ -124,10 +105,12 @@ def determine_gradient(V_g, u, flux):
# Initialize the adapter according to the specific participant
if problem is ProblemType.DIRICHLET:
precice = Adapter(adapter_config_filename="precice-adapter-config-D.json")
- precice_dt = precice.initialize(coupling_boundary, read_function_space=V, write_object=f_N_function)
+ precice_dt = precice.initialize(
+ coupling_boundary, read_function_space=V, write_object=f_N_function)
elif problem is ProblemType.NEUMANN:
precice = Adapter(adapter_config_filename="precice-adapter-config-N.json")
- precice_dt = precice.initialize(coupling_boundary, read_function_space=V_g, write_object=u_D_function)
+ precice_dt = precice.initialize(
+ coupling_boundary, read_function_space=V_g, write_object=u_D_function)
boundary_marker = False
@@ -152,7 +135,8 @@ def determine_gradient(V_g, u, flux):
# modify Neumann boundary condition on coupling interface, modify weak form correspondingly
if not boundary_marker: # there is only 1 Neumann-BC which is at the coupling boundary -> integration over whole boundary
if coupling_expression.is_scalar_valued():
- F += v * coupling_expression * dolfin.ds # this term has to be added to weak form to add a Neumann BC (see e.g. p. 83ff Langtangen, Hans Petter, and Anders Logg. "Solving PDEs in Python The FEniCS Tutorial Volume I." (2016).)
+ # this term has to be added to weak form to add a Neumann BC (see e.g. p. 83ff Langtangen, Hans Petter, and Anders Logg. "Solving PDEs in Python The FEniCS Tutorial Volume I." (2016).)
+ F += v * coupling_expression * dolfin.ds
elif coupling_expression.is_vector_valued():
normal = FacetNormal(mesh)
F += -v * dot(normal, coupling_expression) * dolfin.ds
@@ -199,7 +183,8 @@ def determine_gradient(V_g, u, flux):
error_out << error_pointwise
# set t_1 = t_0 + dt, this gives u_D^1
-u_D.t = t + dt(0) # call dt(0) to evaluate FEniCS Constant. Todo: is there a better way?
+# call dt(0) to evaluate FEniCS Constant. Todo: is there a better way?
+u_D.t = t + dt(0)
f.t = t + dt(0)
flux = Function(V_g)
@@ -207,7 +192,8 @@ def determine_gradient(V_g, u, flux):
while precice.is_coupling_ongoing():
- if precice.is_action_required(precice.action_write_iteration_checkpoint()): # write checkpoint
+ # write checkpoint
+ if precice.is_action_required(precice.action_write_iteration_checkpoint()):
precice.store_checkpoint(u_n, t, n)
read_data = precice.read_data()
@@ -235,7 +221,8 @@ def determine_gradient(V_g, u, flux):
precice_dt = precice.advance(dt(0))
- if precice.is_action_required(precice.action_read_iteration_checkpoint()): # roll back to checkpoint
+ # roll back to checkpoint
+ if precice.is_action_required(precice.action_read_iteration_checkpoint()):
u_cp, t_cp, n_cp = precice.retrieve_checkpoint()
u_n.assign(u_cp)
t = t_cp
@@ -248,7 +235,8 @@ def determine_gradient(V_g, u, flux):
if precice.is_time_window_complete():
u_ref = interpolate(u_D, V)
u_ref.rename("reference", " ")
- error, error_pointwise = compute_errors(u_n, u_ref, V, total_error_tol=error_tol)
+ error, error_pointwise = compute_errors(
+ u_n, u_ref, V, total_error_tol=error_tol)
print('n = %d, t = %.2f: L2 error on domain = %.3g' % (n, t, error))
# output solution and reference solution at t_n+1
print('output u^%d and u_ref^%d' % (n, n))
diff --git a/partitioned-heat-conduction/fenics-fenics/my_enums.py b/partitioned-heat-conduction-complex/fenics/my_enums.py
similarity index 51%
rename from partitioned-heat-conduction/fenics-fenics/my_enums.py
rename to partitioned-heat-conduction-complex/fenics/my_enums.py
index e94988681..27fc5d9ae 100644
--- a/partitioned-heat-conduction/fenics-fenics/my_enums.py
+++ b/partitioned-heat-conduction-complex/fenics/my_enums.py
@@ -1,5 +1,6 @@
from enum import Enum
+
class ProblemType(Enum):
"""
Enum defines problem type. Details see above.
@@ -16,15 +17,3 @@ class DomainPart(Enum):
RIGHT = 2 # right part of domain in simple interface case
CIRCULAR = 3 # circular part of domain in complex interface case
RECTANGLE = 4 # domain excluding circular part of complex interface case
-
-
-class Subcycling(Enum):
- """
- Enum defines which kind of subcycling is used
- """
- NONE = 0 # no subcycling, precice_dt == fenics_dt
- MATCHING = 1 # subcycling, where fenics_dt fits into precice_dt, mod(precice_dt, fenics_dt) == 0
- NONMATCHING = 2 # subcycling, where fenics_dt does not fit into precice_dt, mod(precice_dt, fenics_dt) != 0
-
- # note: the modulo expressions above should be understood in an exact way (no floating point round off problems. For
- # details, see https://stackoverflow.com/questions/14763722/python-modulo-on-floats)
\ No newline at end of file
diff --git a/partitioned-heat-conduction-complex/fenics/precice-adapter-config-D.json b/partitioned-heat-conduction-complex/fenics/precice-adapter-config-D.json
new file mode 100644
index 000000000..9430d6e1e
--- /dev/null
+++ b/partitioned-heat-conduction-complex/fenics/precice-adapter-config-D.json
@@ -0,0 +1,9 @@
+{
+ "participant_name": "Dirichlet",
+ "config_file_name": "../precice-config.xml",
+ "interface": {
+ "coupling_mesh_name": "Dirichlet-Mesh",
+ "write_data_name": "Flux",
+ "read_data_name": "Temperature"
+ }
+}
diff --git a/partitioned-heat-conduction-complex/fenics/precice-adapter-config-N.json b/partitioned-heat-conduction-complex/fenics/precice-adapter-config-N.json
new file mode 100644
index 000000000..b21673b23
--- /dev/null
+++ b/partitioned-heat-conduction-complex/fenics/precice-adapter-config-N.json
@@ -0,0 +1,9 @@
+{
+ "participant_name": "Neumann",
+ "config_file_name": "../precice-config.xml",
+ "interface": {
+ "coupling_mesh_name": "Neumann-Mesh",
+ "write_data_name": "Temperature",
+ "read_data_name": "Flux"
+ }
+}
diff --git a/partitioned-heat-conduction/fenics-fenics/problem_setup.py b/partitioned-heat-conduction-complex/fenics/problem_setup.py
similarity index 100%
rename from partitioned-heat-conduction/fenics-fenics/problem_setup.py
rename to partitioned-heat-conduction-complex/fenics/problem_setup.py
diff --git a/partitioned-heat-conduction-complex/precice-config.xml b/partitioned-heat-conduction-complex/precice-config.xml
new file mode 100644
index 000000000..547a4f207
--- /dev/null
+++ b/partitioned-heat-conduction-complex/precice-config.xml
@@ -0,0 +1,61 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/partitioned-heat-conduction/README.md b/partitioned-heat-conduction/README.md
index 945bc0cf5..a78c6ff99 100644
--- a/partitioned-heat-conduction/README.md
+++ b/partitioned-heat-conduction/README.md
@@ -1,8 +1,70 @@
---
title: Partitioned heat conduction
permalink: tutorials-partitioned-heat-conduction.html
-keywords:
-summary:
+keywords: FEniCS, Nutils, Heat conduction
+summary: In this tutorial we solve a simple heat equation. The domain is partitioned and the coupling is established in a Dirichlet-Neumann fashion.
---
{% include important.html content="We have not yet ported the documentation of the preCICE tutorials from the preCICE wiki to here. Please go to the [preCICE wiki](https://github.com/precice/precice/wiki#2-getting-started---tutorials)" %}
+
+## Setup
+
+In this tutorial we solve a partitioned heat equation. For information on the non-partitioned case, please refer to [1, p.37ff]. In this tutorial the computational domain is partitioned and coupled via preCICE. The coupling roughly follows the approach described in [2].
+
+
+
+Case setup from [3]. `D` denotes the Dirichlet participant and `N` denotes the Neumann participant.
+
+The heat equation is solved on a rectangular domain `Omega = [0,2] x [0,1]` with given Dirichlet boundary conditions. We split the domain at `x_c = 1` using a straight vertical line, the coupling interface. The left part of the domain will be referred to as the Dirichlet partition and the right part as the Neumann partition. To couple the two participant we use Dirichlet-Neumann coupling. Here, the Dirichlet participant receives Dirichlet boundary conditions (`Temperature`) at the coupling interface and solves the heat equation using these boundary conditions on the left part of the domain. Then the Dirichlet participant computes the resulting heat flux (`Flux`) from the solution and sends it to the Neumann participant. The Neumann participant uses the flux as a Neumann boundary condition to solve the heat equation on the right part of the domain. We then extract the temperature from the solution and send it back to the Dirichlet participant. This establishes the coupling between the two participants.
+
+This simple case allows us to compare the solution for the partitioned case to a known analytical solution (method of manufactures solutions, see [1, p.37ff]). For more usage examples and details, please refer to [3, sect. 4.1].
+
+## Available solvers and dependencies
+
+You can either couple a solver with itself or different solvers with each other. In any case you will need to have preCICE and the python bindings installed on your system.
+
+* `fenics`, requires you to install [FEniCS](https://fenicsproject.org/download/) and the [FEniCS-adapter](https://github.com/precice/fenics-adapter). The code is largely based on this [fenics-tutorial](https://github.com/hplgit/fenics-tutorial/blob/master/pub/python/vol1/ft03_heat.py) from [1].
+
+
+* :construction: This case is still under construction. See https://github.com/precice/tutorials/issues/152. :construction: `nutils`, requires you to install [Nutils](http://www.nutils.org/en/latest/).
+
+## Running the simulation
+
+For choosing whether you want to run the Dirichlet-kind and a Neumann-kind participant, please provide the following commandline input:
+
+* `-d` flag will enforce Dirichlet boundary conditions on the coupling interface.
+* `-n` flag will enforce Neumann boundary conditions on the coupling interface.
+
+For running the case, open two terminals and:
+
+```
+cd fenics
+python3 heat.py -d
+```
+
+and
+
+```
+cd fenics
+python3 heat.py -n
+```
+
+If you want to use Nutils for one or both sides of the setup, just `cd nutils`. The FEniCS case also supports parallel runs. Here, simply execute
+
+```
+mpirun -n python3 heat.py -d
+```
+
+## Visualization
+
+Output is written into the folder `fenics/out`. You can visualize the content with paraview by opening the `*.pvd` files. The files `Dirichlet.pvd` and `Neumann.pvd` correspond to the numerical solution of the Dirichlet, respectively Neumann, problem, while the files with the prefix `ref` correspond to the analytical reference solution, the files with `error` show the error and the files with `ranks` the ranks of the solvers (if executed in parallel).
+
+
+
+Visualization in paraview for `x_c = 1.5`.
+
+## References
+
+[1] Hans Petter Langtangen and Anders Logg. "Solving PDEs in Minutes-The FEniCS Tutorial Volume I." (2016). [pdf](https://fenicsproject.org/pub/tutorial/pdf/fenics-tutorial-vol1.pdf)
+[2] Azahar Monge and Philipp Birken. "Convergence Analysis of the Dirichlet-Neumann Iteration for Finite Element Discretizations." (2016). Proceedings in Applied Mathematics and Mechanics. [doi](https://doi.org/10.1002/pamm.201610355)
+[3] Benjamin Rüth, Benjamin Uekermann, Miriam Mehl, Philipp Birken, Azahar Monge, and Hans Joachim Bungartz. "Quasi-Newton waveform iteration for partitioned surface-coupled multiphysics applications." (2020). International Journal for Numerical Methods in Engineering. [doi](https://doi.org/10.1002/nme.6443)
diff --git a/partitioned-heat-conduction/fenics-fenics/README.md b/partitioned-heat-conduction/fenics-fenics/README.md
deleted file mode 100644
index 44a8987c7..000000000
--- a/partitioned-heat-conduction/fenics-fenics/README.md
+++ /dev/null
@@ -1,3 +0,0 @@
-# Tutorial for a partitioned heat equation using FEniCS
-
-This tutorial is described in the [preCICE wiki](https://github.com/precice/precice/wiki/Tutorial-for-solving-the-heat-equation-in-a-partitioned-fashion-using-FEniCS).
diff --git a/partitioned-heat-conduction/fenics-fenics/precice-adapter-config-D.json b/partitioned-heat-conduction/fenics-fenics/precice-adapter-config-D.json
deleted file mode 100644
index 46addb18c..000000000
--- a/partitioned-heat-conduction/fenics-fenics/precice-adapter-config-D.json
+++ /dev/null
@@ -1,9 +0,0 @@
-{
- "participant_name": "HeatDirichlet",
- "config_file_name": "precice-config.xml",
- "interface": {
- "coupling_mesh_name": "DirichletNodes",
- "write_data_name": "Flux",
- "read_data_name": "Temperature"
- }
-}
diff --git a/partitioned-heat-conduction/fenics-fenics/precice-adapter-config-N.json b/partitioned-heat-conduction/fenics-fenics/precice-adapter-config-N.json
deleted file mode 100644
index 335a2b9b8..000000000
--- a/partitioned-heat-conduction/fenics-fenics/precice-adapter-config-N.json
+++ /dev/null
@@ -1,9 +0,0 @@
-{
- "participant_name": "HeatNeumann",
- "config_file_name": "precice-config.xml",
- "interface": {
- "coupling_mesh_name": "NeumannNodes",
- "write_data_name": "Temperature",
- "read_data_name": "Flux"
- }
-}
diff --git a/partitioned-heat-conduction/fenics-fenics/precice-config.xml b/partitioned-heat-conduction/fenics-fenics/precice-config.xml
deleted file mode 100644
index 800f6ade5..000000000
--- a/partitioned-heat-conduction/fenics-fenics/precice-config.xml
+++ /dev/null
@@ -1,61 +0,0 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/partitioned-heat-conduction/fenics/.gitignore b/partitioned-heat-conduction/fenics/.gitignore
new file mode 100644
index 000000000..3cf3e6819
--- /dev/null
+++ b/partitioned-heat-conduction/fenics/.gitignore
@@ -0,0 +1,7 @@
+*.txt
+venv
+.idea
+*.pyc
+out
+*.log
+*events.json
diff --git a/partitioned-heat-conduction/fenics/Allclean b/partitioned-heat-conduction/fenics/Allclean
new file mode 100755
index 000000000..1dbf98f13
--- /dev/null
+++ b/partitioned-heat-conduction/fenics/Allclean
@@ -0,0 +1,11 @@
+#!/bin/bash
+
+echo "Cleaning logs and output files from previous runs..."
+
+rm -r precice-run/
+
+rm out/*
+
+rm -f *.log
+
+echo "Cleaning successful!"
diff --git a/partitioned-heat-conduction/fenics/errorcomputation.py b/partitioned-heat-conduction/fenics/errorcomputation.py
new file mode 100644
index 000000000..8964cd884
--- /dev/null
+++ b/partitioned-heat-conduction/fenics/errorcomputation.py
@@ -0,0 +1,13 @@
+from fenics import inner, assemble, dx, project, sqrt
+
+
+def compute_errors(u_approx, u_ref, V, total_error_tol=10 ** -4):
+ error_normalized = (u_ref - u_approx) / u_ref # compute pointwise L2 error
+ error_pointwise = project(abs(error_normalized), V) # project onto function space
+ error_total = sqrt(
+ assemble(inner(error_pointwise, error_pointwise) * dx)) # determine L2 norm to estimate total error
+ error_pointwise.rename("error", " ")
+
+ assert (error_total < total_error_tol)
+
+ return error_total, error_pointwise
diff --git a/partitioned-heat-conduction/fenics/heat.py b/partitioned-heat-conduction/fenics/heat.py
new file mode 100644
index 000000000..20bfcd130
--- /dev/null
+++ b/partitioned-heat-conduction/fenics/heat.py
@@ -0,0 +1,234 @@
+"""
+The basic example is taken from "Langtangen, Hans Petter, and Anders Logg. Solving PDEs in Python: The FEniCS
+Tutorial I. Springer International Publishing, 2016."
+
+The example code has been extended with preCICE API calls and mixed boundary conditions to allow for a Dirichlet-Neumann
+coupling of two separate heat equations.
+
+The original source code can be found on https://github.com/hplgit/fenics-tutorial/blob/master/pub/python/vol1/ft03_heat.py.
+
+Heat equation with Dirichlet conditions. (Dirichlet problem)
+ u'= Laplace(u) + f in the unit square [0,1] x [0,1]
+ u = u_C on the coupling boundary at x = 1
+ u = u_D on the remaining boundary
+ u = u_0 at t = 0
+ u = 1 + x^2 + alpha*y^2 + \beta*t
+ f = beta - 2 - 2*alpha
+
+Heat equation with mixed boundary conditions. (Neumann problem)
+ u'= Laplace(u) + f in the shifted unit square [1,2] x [0,1]
+ du/dn = f_N on the coupling boundary at x = 1
+ u = u_D on the remaining boundary
+ u = u_0 at t = 0
+ u = 1 + x^2 + alpha*y^2 + \beta*t
+ f = beta - 2 - 2*alpha
+"""
+
+from __future__ import print_function, division
+from fenics import Function, FunctionSpace, Expression, Constant, DirichletBC, TrialFunction, TestFunction, \
+ File, solve, lhs, rhs, grad, inner, dot, dx, ds, interpolate, VectorFunctionSpace, MeshFunction, MPI
+from fenicsprecice import Adapter
+from errorcomputation import compute_errors
+from my_enums import ProblemType, DomainPart
+import argparse
+import numpy as np
+from problem_setup import get_geometry
+import dolfin
+from dolfin import FacetNormal, dot
+
+
+def determine_gradient(V_g, u, flux):
+ """
+ compute flux following http://hplgit.github.io/INF5620/doc/pub/fenics_tutorial1.1/tu2.html#tut-poisson-gradu
+ :param V_g: Vector function space
+ :param u: solution where gradient is to be determined
+ :param flux: returns calculated flux into this value
+ """
+
+ w = TrialFunction(V_g)
+ v = TestFunction(V_g)
+
+ a = inner(w, v) * dx
+ L = inner(grad(u), v) * dx
+ solve(a == L, flux)
+
+
+parser = argparse.ArgumentParser(
+ description='Solving heat equation for simple or complex interface case')
+parser.add_argument("-d", "--dirichlet", help="create a dirichlet problem",
+ dest='dirichlet', action='store_true')
+parser.add_argument("-n", "--neumann", help="create a neumann problem",
+ dest='neumann', action='store_true')
+
+args = parser.parse_args()
+
+fenics_dt = .1 # time step size
+# Error is bounded by coupling accuracy. In theory we would obtain the analytical solution.
+error_tol = 10 ** -6
+
+alpha = 3 # parameter alpha
+beta = 1.3 # parameter beta
+
+if args.dirichlet and not args.neumann:
+ problem = ProblemType.DIRICHLET
+ domain_part = DomainPart.LEFT
+elif args.neumann and not args.dirichlet:
+ problem = ProblemType.NEUMANN
+ domain_part = DomainPart.RIGHT
+
+mesh, coupling_boundary, remaining_boundary = get_geometry(domain_part)
+
+# Define function space using mesh
+V = FunctionSpace(mesh, 'P', 2)
+V_g = VectorFunctionSpace(mesh, 'P', 1)
+W = V_g.sub(0).collapse()
+
+# Define boundary conditions
+u_D = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t',
+ degree=2, alpha=alpha, beta=beta, t=0)
+u_D_function = interpolate(u_D, V)
+
+if problem is ProblemType.DIRICHLET:
+ # Define flux in x direction
+ f_N = Expression("2 * x[0]", degree=1, alpha=alpha, t=0)
+ f_N_function = interpolate(f_N, W)
+
+# Define initial value
+u_n = interpolate(u_D, V)
+u_n.rename("Temperature", "")
+
+precice, precice_dt, initial_data = None, 0.0, None
+
+# Initialize the adapter according to the specific participant
+if problem is ProblemType.DIRICHLET:
+ precice = Adapter(adapter_config_filename="precice-adapter-config-D.json")
+ precice_dt = precice.initialize(
+ coupling_boundary, read_function_space=V, write_object=f_N_function)
+elif problem is ProblemType.NEUMANN:
+ precice = Adapter(adapter_config_filename="precice-adapter-config-N.json")
+ precice_dt = precice.initialize(
+ coupling_boundary, read_function_space=W, write_object=u_D_function)
+
+dt = Constant(0)
+dt.assign(np.min([fenics_dt, precice_dt]))
+
+# Define variational problem
+u = TrialFunction(V)
+v = TestFunction(V)
+f = Expression('beta - 2 - 2*alpha', degree=2, alpha=alpha, beta=beta, t=0)
+F = u * v / dt * dx + dot(grad(u), grad(v)) * dx - (u_n / dt + f) * v * dx
+
+bcs = [DirichletBC(V, u_D, remaining_boundary)]
+
+# Set boundary conditions at coupling interface once wrt to the coupling
+# expression
+coupling_expression = precice.create_coupling_expression()
+if problem is ProblemType.DIRICHLET:
+ # modify Dirichlet boundary condition on coupling interface
+ bcs.append(DirichletBC(V, coupling_expression, coupling_boundary))
+if problem is ProblemType.NEUMANN:
+ # modify Neumann boundary condition on coupling interface, modify weak
+ # form correspondingly
+ F += v * coupling_expression * dolfin.ds
+
+a, L = lhs(F), rhs(F)
+
+# Time-stepping
+u_np1 = Function(V)
+u_np1.rename("Temperature", "")
+t = 0
+
+# reference solution at t=0
+u_ref = interpolate(u_D, V)
+u_ref.rename("reference", " ")
+
+# mark mesh w.r.t ranks
+mesh_rank = MeshFunction("size_t", mesh, mesh.topology().dim())
+if problem is ProblemType.NEUMANN:
+ mesh_rank.set_all(MPI.rank(MPI.comm_world) + 4)
+else:
+ mesh_rank.set_all(MPI.rank(MPI.comm_world) + 0)
+mesh_rank.rename("myRank", "")
+
+# Generating output files
+temperature_out = File("out/%s.pvd" % precice.get_participant_name())
+ref_out = File("out/ref%s.pvd" % precice.get_participant_name())
+error_out = File("out/error%s.pvd" % precice.get_participant_name())
+ranks = File("out/ranks%s.pvd" % precice.get_participant_name())
+
+# output solution and reference solution at t=0, n=0
+n = 0
+print('output u^%d and u_ref^%d' % (n, n))
+temperature_out << u_n
+ref_out << u_ref
+ranks << mesh_rank
+
+error_total, error_pointwise = compute_errors(u_n, u_ref, V)
+error_out << error_pointwise
+
+# set t_1 = t_0 + dt, this gives u_D^1
+# call dt(0) to evaluate FEniCS Constant. Todo: is there a better way?
+u_D.t = t + dt(0)
+f.t = t + dt(0)
+
+if problem is ProblemType.DIRICHLET:
+ flux = Function(V_g)
+ flux.rename("Flux", "")
+
+while precice.is_coupling_ongoing():
+
+ # write checkpoint
+ if precice.is_action_required(precice.action_write_iteration_checkpoint()):
+ precice.store_checkpoint(u_n, t, n)
+
+ read_data = precice.read_data()
+
+ # Update the coupling expression with the new read data
+ precice.update_coupling_expression(coupling_expression, read_data)
+
+ dt.assign(np.min([fenics_dt, precice_dt]))
+
+ # Compute solution u^n+1, use bcs u_D^n+1, u^n and coupling bcs
+ solve(a == L, u_np1, bcs)
+
+ # Write data to preCICE according to which problem is being solved
+ if problem is ProblemType.DIRICHLET:
+ # Dirichlet problem reads temperature and writes flux on boundary to Neumann problem
+ determine_gradient(V_g, u_np1, flux)
+ flux_x = interpolate(flux.sub(0), W)
+ precice.write_data(flux_x)
+ elif problem is ProblemType.NEUMANN:
+ # Neumann problem reads flux and writes temperature on boundary to Dirichlet problem
+ precice.write_data(u_np1)
+
+ precice_dt = precice.advance(dt(0))
+
+ # roll back to checkpoint
+ if precice.is_action_required(precice.action_read_iteration_checkpoint()):
+ u_cp, t_cp, n_cp = precice.retrieve_checkpoint()
+ u_n.assign(u_cp)
+ t = t_cp
+ n = n_cp
+ else: # update solution
+ u_n.assign(u_np1)
+ t += float(dt)
+ n += 1
+
+ if precice.is_time_window_complete():
+ u_ref = interpolate(u_D, V)
+ u_ref.rename("reference", " ")
+ error, error_pointwise = compute_errors(
+ u_n, u_ref, V, total_error_tol=error_tol)
+ print('n = %d, t = %.2f: L2 error on domain = %.3g' % (n, t, error))
+ # output solution and reference solution at t_n+1
+ print('output u^%d and u_ref^%d' % (n, n))
+ temperature_out << u_n
+ ref_out << u_ref
+ error_out << error_pointwise
+
+ # Update Dirichlet BC
+ u_D.t = t + float(dt)
+ f.t = t + float(dt)
+
+# Hold plot
+precice.finalize()
diff --git a/partitioned-heat-conduction/fenics/my_enums.py b/partitioned-heat-conduction/fenics/my_enums.py
new file mode 100644
index 000000000..aa558178b
--- /dev/null
+++ b/partitioned-heat-conduction/fenics/my_enums.py
@@ -0,0 +1,18 @@
+from enum import Enum
+
+class ProblemType(Enum):
+ """
+ Enum defines problem type. Details see above.
+ """
+ DIRICHLET = 1 # Dirichlet problem
+ NEUMANN = 2 # Neumann problem
+
+
+class DomainPart(Enum):
+ """
+ Enum defines which part of the domain [x_left, x_right] x [y_bottom, y_top] we compute.
+ """
+ LEFT = 1 # left part of domain in simple interface case
+ RIGHT = 2 # right part of domain in simple interface case
+ CIRCULAR = 3 # circular part of domain in complex interface case
+ RECTANGLE = 4 # domain excluding circular part of complex interface case
diff --git a/partitioned-heat-conduction/fenics/precice-adapter-config-D.json b/partitioned-heat-conduction/fenics/precice-adapter-config-D.json
new file mode 100644
index 000000000..c2a7788c0
--- /dev/null
+++ b/partitioned-heat-conduction/fenics/precice-adapter-config-D.json
@@ -0,0 +1,9 @@
+{
+ "participant_name": "Dirichlet",
+ "config_file_name": "../precice-config.xml",
+ "interface": {
+ "coupling_mesh_name": "Dirichlet-Mesh",
+ "write_data_name": "Flux",
+ "read_data_name": "Temperature"
+ }
+}
diff --git a/partitioned-heat-conduction/fenics/precice-adapter-config-N.json b/partitioned-heat-conduction/fenics/precice-adapter-config-N.json
new file mode 100644
index 000000000..36ea15d19
--- /dev/null
+++ b/partitioned-heat-conduction/fenics/precice-adapter-config-N.json
@@ -0,0 +1,9 @@
+{
+ "participant_name": "Neumann",
+ "config_file_name": "../precice-config.xml",
+ "interface": {
+ "coupling_mesh_name": "Neumann-Mesh",
+ "write_data_name": "Temperature",
+ "read_data_name": "Flux"
+ }
+}
diff --git a/partitioned-heat-conduction/fenics/problem_setup.py b/partitioned-heat-conduction/fenics/problem_setup.py
new file mode 100644
index 000000000..97610a7d2
--- /dev/null
+++ b/partitioned-heat-conduction/fenics/problem_setup.py
@@ -0,0 +1,59 @@
+"""
+Problem setup for partitioned-heat-conduction/fenics-fenics tutorial
+"""
+
+from fenics import SubDomain, Point, RectangleMesh, near, Function, VectorFunctionSpace, Expression
+from my_enums import DomainPart
+
+
+y_bottom, y_top = 0, 1
+x_left, x_right = 0, 2
+x_coupling = 1.0 # x coordinate of coupling interface
+radius = 0.2
+midpoint = Point(0.5, 0.5)
+
+
+class ExcludeStraightBoundary(SubDomain):
+ def get_user_input_args(self, args):
+ self._interface = args.interface
+
+ def inside(self, x, on_boundary):
+ tol = 1E-14
+ if on_boundary and not near(x[0], x_coupling, tol) or near(x[1], y_top, tol) or near(x[1], y_bottom, tol):
+ return True
+ else:
+ return False
+
+
+class StraightBoundary(SubDomain):
+ def inside(self, x, on_boundary):
+ tol = 1E-14
+ if on_boundary and near(x[0], x_coupling, tol):
+ return True
+ else:
+ return False
+
+
+def get_geometry(domain_part):
+ nx = 5
+ ny = 10
+ low_resolution = 5
+ high_resolution = 5
+ n_vertices = 20
+
+ if domain_part is DomainPart.LEFT:
+ nx = nx*3
+ p0 = Point(x_left, y_bottom)
+ p1 = Point(x_coupling, y_top)
+ elif domain_part is DomainPart.RIGHT:
+ ny = ny*2
+ p0 = Point(x_coupling, y_bottom)
+ p1 = Point(x_right, y_top)
+ else:
+ raise Exception("invalid domain_part: {}".format(domain_part))
+
+ mesh = RectangleMesh(p0, p1, nx, ny)
+ coupling_boundary = StraightBoundary()
+ remaining_boundary = ExcludeStraightBoundary()
+
+ return mesh, coupling_boundary, remaining_boundary
diff --git a/partitioned-heat-conduction/images/HT_FEniCS_movie.gif b/partitioned-heat-conduction/images/HT_FEniCS_movie.gif
new file mode 100644
index 000000000..55a7fc8c3
Binary files /dev/null and b/partitioned-heat-conduction/images/HT_FEniCS_movie.gif differ
diff --git a/partitioned-heat-conduction/images/setup.png b/partitioned-heat-conduction/images/setup.png
new file mode 100644
index 000000000..e442f20c4
Binary files /dev/null and b/partitioned-heat-conduction/images/setup.png differ
diff --git a/partitioned-heat-conduction/precice-config.xml b/partitioned-heat-conduction/precice-config.xml
new file mode 100644
index 000000000..400214c2f
--- /dev/null
+++ b/partitioned-heat-conduction/precice-config.xml
@@ -0,0 +1,61 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+