Merge pull request #150 from precice/fenics-adapter-v1.3.0

Release v1.3.0

CHANGELOG.md

@@ -1,5 +1,10 @@
 # FEniCS-preCICE adapter changelog
 
+## 1.3.0
+
+* Adding functionality for 3D cases with PointSource objects at coupling boundaries. See PRs [#133](https://github.com/precice/fenics-adapter/pull/133), [#146](https://github.com/precice/fenics-adapter/pull/146) and [#147](https://github.com/precice/fenics-adapter/pull/147).
+* Fixing an issue of the `config` object not being able to find the adapter configuration file in a Jupyter notebook. [See PR #144](https://github.com/precice/fenics-adapter/pull/144)
+
 ## 1.2.0
 
 * Reduce complexity of initialization to reduce runtime for large cases. [See PR #135](https://github.com/precice/fenics-adapter/pull/135)

README.md

@@ -14,7 +14,7 @@
     <img src="https://github.com/precice/fenics-adapter/actions/workflows/pythonpublish.yml/badge.svg" alt="Upload Python Package">
 </a>
 
-preCICE-adapter for the open source computing platform FEniCS. Note: The adapter **currently only supports 2D simulations in FEniCS.**
+preCICE-adapter for the open source computing platform FEniCS.
 
 ## Installing the package
 
@@ -75,18 +75,14 @@ To create and install the `fenicsprecice` python package the following instructi
 
 ## Citing
 
-If you are using this adapter, please refer to the [citing information on the FEniCS adapter](https://www.precice.org/adapter-fenics.html#how-to-cite).
-
-preCICE is an academic project, developed at the [Technical University of Munich](https://www5.in.tum.de/) and at the [University of Stuttgart](https://www.ipvs.uni-stuttgart.de/). If you use preCICE, please [cite us](https://www.precice.org/publications/):
-
-*H.-J. Bungartz, F. Lindner, B. Gatzhammer, M. Mehl, K. Scheufele, A. Shukaev, and B. Uekermann: preCICE - A Fully Parallel Library for Multi-Physics Surface Coupling. Computers and Fluids, 141, 250–258, 2016.*
-
-If you are using FEniCS, please also consider the information on [the official FEniCS website on citing](https://fenicsproject.org/citing/).
+* FEniCS-preCICE: If you are using this adapter (`fenics-adapter`), please refer to the [citing information on the FEniCS adapter](https://www.precice.org/adapter-fenics.html#how-to-cite).
+* preCICE: preCICE is an academic project, developed at the [Technical University of Munich](https://www5.in.tum.de/) and at the [University of Stuttgart](https://www.ipvs.uni-stuttgart.de/). If you use preCICE, please [cite preCICE](https://precice.org/publications.html#how-to-cite-precice).
+* FEniCS: If you are using FEniCS, please also consider the information on [the official FEniCS website on citing](https://fenicsproject.org/citing/).
 
 ## Development history
 
 The initial version of this adapter was developed by [Benjamin Rodenberg](https://www.in.tum.de/i05/personen/personen/benjamin-rodenberg/) during his research stay at Lund University in the group for [Numerical Analysis](http://www.maths.lu.se/english/research/research-divisions/numerical-analysis/) in close collaboration with [Peter Meisrimel](https://www.lunduniversity.lu.se/lucat/user/09d80f0367a060bcf2a22d7c22e5e504).
 
 [Richard Hertrich](https://github.com/richahert) contributed the possibility to perform FSI simulations using the adapter in his [Bachelor thesis](https://mediatum.ub.tum.de/node?id=1520579).
 
-[Ishaan Desai](https://www.ipvs.uni-stuttgart.de/institute/team/Desai/) improved the user interface and extended the adapter to also allow for parallel FEniCS computations.
+[Ishaan Desai](https://www.ipvs.uni-stuttgart.de/institute/team/Desai/) improved the user interface and extended the adapter to allow for parallel FEniCS computations and 3D cases in certain scenarios.

fenicsprecice/adapter_core.py

@@ -78,24 +78,28 @@ class CouplingMode(Enum):
     UNI_DIRECTIONAL_READ_COUPLING = 6
 
 
-def determine_function_type(input_obj):
+def determine_function_type(input_obj, dims):
     """
     Determines if the function is scalar- or vector-valued based on rank evaluation.
 
     Parameters
     ----------
-    input_obj :
-        A FEniCS function.
+    input_obj : FEniCS Function or FEniCS FunctionSpace
+        A FEniCS Function object or a FEniCS FunctionSpace object
+    dims : int
+        Dimension of problem.
 
     Returns
     -------
     tag : bool
         0 if input_function is SCALAR and 1 if input_function is VECTOR.
     """
-    if isinstance(input_obj, FunctionSpace):  # scalar-valued functions have rank 0 is FEniCS
-        if input_obj.num_sub_spaces() == 0:
+    if isinstance(input_obj, FunctionSpace):
+        obj_dim = input_obj.num_sub_spaces()
+        if obj_dim == 0:
             return FunctionType.SCALAR
-        elif input_obj.num_sub_spaces() == 2:
+        elif obj_dim == 2 or obj_dim == 3:
+            assert obj_dim == dims
             return FunctionType.VECTOR
     elif isinstance(input_obj, Function):
         if input_obj.value_rank() == 0:
@@ -120,7 +124,7 @@ def filter_point_sources(point_sources, filter_out, warn_duplicate=True):
     filter_out: FEniCS domain
         Defines the domain where PointSources should be filtered out.
     warn_duplicate: bool
-        Set False to surpress warnings, if double-boundary points are filtered out.
+        Set False to suppress warnings, if double-boundary points are filtered out.
 
     Returns
     -------
@@ -218,10 +222,7 @@ def get_fenics_vertices(function_space, coupling_subdomain, dims):
         if coupling_subdomain.inside(v.point(), True):
             lids.append(v.index())
             gids.append(v.global_index())
-            if dims == 2:
-                coords.append([v.x(0), v.x(1)])
-            if dims == 3:
-                coords.append([v.x(0), v.x(1), v.x(2)])
+            coords.append([v.x(d) for d in range(dims)])
 
     return np.array(lids), np.array(gids), np.array(coords)
 
@@ -280,13 +281,8 @@ def get_owned_vertices(function_space, coupling_subdomain, dims):
     # Get coordinates and global IDs of all vertices of the mesh  which lie on the coupling boundary.
     # These vertices include shared (owned + unowned) and non-shared vertices in a parallel setting
     gids, lids, coords = [], [], []
-    coord = None
     for v in coupling_vertices:
-        if dims == 2:
-            coord = [v.x(0), v.x(1)]
-        elif dims == 3:
-            coord = [v.x(0), v.x(1), v.x(2)]
-
+        coord = [v.x(d) for d in range(dims)]
         for dof in dofs:
             if (dof == coord).all():
                 gids.append(v.global_index())
@@ -344,14 +340,9 @@ def get_unowned_vertices(function_space, coupling_subdomain, dims):
     # Get coordinates and global IDs of all vertices of the mesh  which lie on the coupling boundary.
     # These vertices include shared (owned + unowned) and non-shared vertices in a parallel setting
     gids = []
-    coord = None
     for v in coupling_verts:
         ownership = False
-        if dims == 2:
-            coord = [v.x(0), v.x(1)]
-        elif dims == 3:
-            coord = [v.x(0), v.x(1), v.x(2)]
-
+        coord = [v.x(d) for d in range(dims)]
         for dof in dofs:
             if (dof == coord).all():
                 ownership = True
@@ -386,12 +377,12 @@ def get_coupling_boundary_edges(function_space, coupling_subdomain, global_ids,
         Array of second vertex of each edge.
     """
 
-    def edge_is_on(subdomain, edge):
+    def edge_is_on(subdomain, this_edge):
         """
         Check whether edge lies within subdomain
         """
-        assert(len(list(vertices(edge))) == 2)
-        return all([subdomain.inside(v.point(), True) for v in vertices(edge)])
+        assert(len(list(vertices(this_edge))) == 2)
+        return all([subdomain.inside(v.point(), True) for v in vertices(this_edge)])
 
     vertices1_ids = []
     vertices2_ids = []
@@ -427,34 +418,34 @@ def get_forces_as_point_sources(fixed_boundary, function_space, data):
 
     Returns
     -------
-    x_forces : list
-        Dictionary carrying X component of forces with reference to each point on the coupling interface.
-    y_forces : list
-        Dictionary carrying Y component of forces with reference to each point on the coupling interface.
+    forces : list
+        D number of lists carrying components of forces with reference to points on the coupling interface.
+        D is dimension of the problem.
     """
-    x_forces = dict()  # dict of PointSources for Forces in x direction
-    y_forces = dict()  # dict of PointSources for Forces in y direction
-
-    fenics_vertices = np.array(list(data.keys()))
+    vertices = np.array(list(data.keys()))
     nodal_data = np.array(list(data.values()))
 
-    # Check for shape of coupling_mesh_vertices and raise Assertion for 3D
-    n_vertices, _ = fenics_vertices.shape
+    n_vertices, dims = vertices.shape
+    assert (dims == 2 or dims == 3), "Provided data does not have the correct dimensions"
 
-    vertices_x = fenics_vertices[:, 0]
-    vertices_y = fenics_vertices[:, 1]
+    forces = []
+    for d in range(dims):
+        forces.append(dict())
 
     for i in range(n_vertices):
-        px, py = vertices_x[i], vertices_y[i]
-        key = (px, py)
-        x_forces[key] = PointSource(function_space.sub(0), Point(px, py), nodal_data[i, 0])
-        y_forces[key] = PointSource(function_space.sub(1), Point(px, py), nodal_data[i, 1])
+        key = []
+        for d in range(dims):
+            key.append(vertices[i, d])
+        key = tuple(key)
+
+        for d in range(dims):
+            forces[d][key] = PointSource(function_space.sub(d), Point(key), nodal_data[i, d])
 
     # Avoid application of PointSource and Dirichlet boundary condition at the same point by filtering
-    x_forces = filter_point_sources(x_forces, fixed_boundary, warn_duplicate=False)
-    y_forces = filter_point_sources(y_forces, fixed_boundary, warn_duplicate=False)
+    for d in range(dims):
+        forces[d] = filter_point_sources(forces[d], fixed_boundary, warn_duplicate=False)
 
-    return x_forces.values(), y_forces.values()  # don't return dictionary, but list of PointSources
+    return (forces[d].values() for d in range(dims))
 
 
 def get_communication_map(comm, function_space, owned_vertices, unowned_vertices):

fenicsprecice/config.py

@@ -33,8 +33,8 @@ def read_json(self, adapter_config_filename):
         adapter_config_filename : string
             Name of the JSON configuration file
         """
-        folder = os.path.dirname(os.path.join(os.getcwd(), os.path.dirname(sys.argv[0]), adapter_config_filename))
-        path = os.path.join(folder, os.path.basename(adapter_config_filename))
+        path = os.path.abspath(adapter_config_filename)
+        folder = os.path.dirname(path)
         read_file = open(path, "r")
         data = json.load(read_file)
         self._config_file_name = os.path.join(folder, data["config_file_name"])

fenicsprecice/expression_core.py

@@ -20,37 +20,39 @@ class CouplingExpression(UserExpression):
     Creates functional representation (for FEniCS) of nodal data provided by preCICE.
     """
 
+    def __init__(self, element, degree):
+        super().__init__(element=element, degree=degree)  # Call constuctor of UserExpression class
+
+        self._function_type = None
+        self._vals = None
+        self._dimension = None
+
+        self._coords_x = self._coords_y = None
+
+        self._f = None
+
     def set_function_type(self, function_type):
         self._function_type = function_type
 
-    def update_boundary_data(self, vals, coords_x, coords_y=None, coords_z=None):
+    def update_boundary_data(self, vals, coords):
         """
         Update object of this class of type FEniCS UserExpression with given point data.
 
         Parameters
         ----------
         vals : double
             Point data to be used to update the Expression.
-        coords_x : double
-            X coordinate of points of which point data is provided.
-        coords_y : double
-            Y coordinate of points of which point data is provided.
-        coords_z : double
-            Z coordinate of points of which point data is provided.
-        """
-        self._coords_x = coords_x
-        self._dimension = 3
-        if coords_y is None:
-            self._dimension -= 1
-            coords_y = np.zeros(self._coords_x.shape)
-        self._coords_y = coords_y
-        if coords_z is None:
-            self._dimension -= 1
-            coords_z = np.zeros(self._coords_x.shape)
-
-        self._coords_y = coords_y
-        self._coords_z = coords_z
+        coords : numpy array
+            The coordinates of fenics vertices in a numpy array [N x D] where
+            N = number of vertices and D = dimensions of geometry.
+        """
         self._vals = vals
+        _, self._dimension = coords.shape
+
+        assert(self._dimension == 2), "Coordinates are of incorrect dimensions"
+
+        self._coords_x = coords[:, 0]
+        self._coords_y = coords[:, 1]
 
         self._f = self.create_interpolant()