Merge pull request #156 from precice/fenics-adapter-v1.4.0

Release v1.4.0

CHANGELOG.md

@@ -1,5 +1,10 @@
 # FEniCS-preCICE adapter changelog
 
+## 1.4.0
+
+* Adding CITATION.cff to link the adapter repository to the relevant publication in the journal SoftwareX.
+* Add functionality to define mesh connectivity in 2D cases in the form of triangles.
+
 ## 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).

CITATION.cff

@@ -0,0 +1,45 @@
+# YAML 1.2
+---
+abstract: "The new software FEniCS-preCICE is a middle software layer, sitting in between the existing finite-element library FEniCS and the coupling library preCICE. The middle layer simplifies coupling (existing) FEniCS application codes to other simulation software via preCICE. To this end, FEniCS-preCICE converts between FEniCS and preCICE mesh and data structures, provides easy-to-use coupling conditions, and manages data checkpointing for implicit coupling. The new software is a library itself and follows a FEniCS-native style. Only a few lines of additional code are necessary to prepare a FEniCS application code for coupling. We illustrate the functionality of FEniCS-preCICE by two examples: a FEniCS heat conduction code coupled to OpenFOAM and a FEniCS linear elasticity code coupled to SU2. The results of both scenarios are compared with other simulation software showing good agreement."
+authors: 
+  -
+    affiliation: "Technical University of Munich"
+    family-names: Rodenberg
+    given-names: Benjamin
+    orcid: "https://orcid.org/0000-0002-3116-0133"
+  -
+    affiliation: "University Stuttgart"
+    family-names: Desai
+    given-names: Ishaan
+    orcid: "https://orcid.org/0000-0002-2552-7509"
+  -
+    family-names: Hertrich
+    given-names: Richard
+    orcid: "https://orcid.org/0000-0003-1722-2841"
+  -
+    affiliation: "University of Stuttgart"
+    family-names: Jaust
+    given-names: Alexander
+    orcid: "https://orcid.org/0000-0002-6082-105X"
+  -
+    affiliation: "University of Stuttgart"
+    family-names: Uekermann
+    given-names: Benjamin
+    orcid: "https://orcid.org/0000-0002-1314-9969"
+cff-version: "1.1.0"
+date-released: 2021-01-10
+keywords: 
+  - FEniCS
+  - "Fluid-Structure Interaction"
+  - "Conjugate Heat Transfer"
+  - Multiphysics
+  - "Coupled Problems"
+  - "Finite Element Method"
+  - preCICE
+license: "LGPL-3.0"
+message: "If you use this software, please cite it using these metadata."
+repository-code: "https://github.com/precice/fenics-adapter"
+title: "FEniCS-preCICE: Coupling FEniCS to other Simulation Software"
+version: 1.2.0
+doi: 10.1016/j.softx.2021.100807
+...

README.md

@@ -4,6 +4,10 @@
     <img src="https://img.shields.io/github/license/precice/fenics-adapter.svg" alt="GNU LGPL license">
 </a>
 
+<a style="text-decoration: none" href="https://doi.org/10.1016/j.softx.2021.100807" target="_blank">
+    <img src="https://zenodo.org/badge/DOI/10.1016/j.softx.2021.100807.svg" alt="DOI">
+</a>
+
 <a style="text-decoration: none" href="https://github.com/precice/fenics-adapter/actions/workflows/build-and-test.yml" target="_blank">
     <img src="https://github.com/precice/fenics-adapter/actions/workflows/build-and-test.yml/badge.svg" alt="Build and Test">
 </a>

fenicsprecice/adapter_core.py

@@ -2,7 +2,7 @@
 This module consists of helper functions used in the Adapter class. Names of the functions are self explanatory
 """
 
-from fenics import SubDomain, Point, PointSource, vertices, FunctionSpace, Function, edges
+from fenics import SubDomain, Point, PointSource, vertices, FunctionSpace, Function, edges, cells
 import numpy as np
 from enum import Enum
 import logging
@@ -375,6 +375,8 @@ def get_coupling_boundary_edges(function_space, coupling_subdomain, global_ids,
         Array of first vertex of each edge.
     vertices2_ids : numpy array
         Array of second vertex of each edge.
+    edges_ids : numpy array
+        Array of FEniCS edge local IDs.
     """
 
     def edge_is_on(subdomain, this_edge):
@@ -386,18 +388,58 @@ def edge_is_on(subdomain, this_edge):
 
     vertices1_ids = []
     vertices2_ids = []
+    edges_ids = []
 
     for edge in edges(function_space.mesh()):
         if edge_is_on(coupling_subdomain, edge):
             v1, v2 = list(vertices(edge))
             if v1.global_index() in global_ids and v2.global_index() in global_ids:
                 vertices1_ids.append(id_mapping[v1.global_index()])
                 vertices2_ids.append(id_mapping[v2.global_index()])
+                edges_ids.append(edge.index())
 
     vertices1_ids = np.array(vertices1_ids)
     vertices2_ids = np.array(vertices2_ids)
+    edges_ids = np.array(edges_ids)
 
-    return vertices1_ids, vertices2_ids
+    return vertices1_ids, vertices2_ids, edges_ids
+
+
+def get_coupling_triangles(function_space, coupling_subdomain, precice_edge_dict):
+    """
+    Extracts triangles of mesh which lie on the coupling region.
+
+    Parameters
+    ----------
+    function_space : FEniCS function space
+        Function space on which the finite element problem definition lives.
+    coupling_subdomain : FEniCS Domain
+        FEniCS domain of the coupling interface region.
+    precice_edge_dict: dict
+        Dictionary with FEniCS IDs of coupling mesh edges as keys and preCICE IDs of the edges as values
+
+    Returns
+    -------
+    edges : numpy array
+        Array of edges indices (3 per triangle)
+    """
+
+    def cell_is_in(subdomain, this_cell):
+        """
+        Check whether edge lies within subdomain
+        """
+        assert(len(list(vertices(this_cell))) == 3), "Only triangular meshes are supported"
+        return all([subdomain.inside(v.point(), True) for v in vertices(this_cell)])
+
+    edges_ids = []
+
+    for cell in cells(function_space.mesh()):
+        if cell_is_in(coupling_subdomain, cell):
+            e1, e2, e3 = list(edges(cell))
+            if all(edge in precice_edge_dict.keys() for edge in [e1.index(), e2.index(), e3.index()]):
+                edges_ids.append([e1.index(), e2.index(), e3.index()])
+
+    return np.array(edges_ids)
 
 
 def get_forces_as_point_sources(fixed_boundary, function_space, data):

fenicsprecice/fenicsprecice.py

@@ -8,7 +8,7 @@
 import precice
 from .adapter_core import FunctionType, determine_function_type, convert_fenics_to_precice, get_fenics_vertices, \
     get_owned_vertices, get_unowned_vertices, get_coupling_boundary_edges, get_forces_as_point_sources, \
-    get_communication_map, communicate_shared_vertices, CouplingMode, Vertices, VertexType, filter_point_sources
+    get_communication_map, communicate_shared_vertices, CouplingMode, Vertices, VertexType, filter_point_sources, get_coupling_triangles
 from .expression_core import SegregatedRBFInterpolationExpression, EmptyExpression
 from .solverstate import SolverState
 from fenics import Function, FunctionSpace
@@ -65,6 +65,7 @@ def __init__(self, adapter_config_filename='precice-adapter-config.json'):
         self._unowned_vertices = Vertices(VertexType.UNOWNED)
         self._fenics_vertices = Vertices(VertexType.FENICS)
         self._precice_vertex_ids = None  # initialized later
+        self._precice_edge_dict = dict()
 
         # read data related quantities (read data is read from preCICE and applied in FEniCS)
         self._read_function_type = None  # stores whether read function is scalar or vector valued
@@ -410,14 +411,25 @@ def initialize(self, coupling_subdomain, read_function_space=None, write_object=
         # Define a mapping between coupling vertices and their IDs in preCICE
         id_mapping = {key: value for key, value in zip(self._owned_vertices.get_global_ids(), self._precice_vertex_ids)}
 
-        edge_vertex_ids1, edge_vertex_ids2 = get_coupling_boundary_edges(function_space, coupling_subdomain,
-                                                                         self._owned_vertices.get_global_ids(),
-                                                                         id_mapping)
+        edge_vertex_ids1, edge_vertex_ids2, edges_ids = get_coupling_boundary_edges(
+            function_space, coupling_subdomain, self._owned_vertices.get_global_ids(), id_mapping)
 
         for i in range(len(edge_vertex_ids1)):
             assert (edge_vertex_ids1[i] != edge_vertex_ids2[i])
-            self._interface.set_mesh_edge(self._interface.get_mesh_id(self._config.get_coupling_mesh_name()),
-                                          edge_vertex_ids1[i], edge_vertex_ids2[i])
+            self._precice_edge_dict[edges_ids[i]] = self._interface.set_mesh_edge(
+                self._interface.get_mesh_id(self._config.get_coupling_mesh_name()),
+                edge_vertex_ids1[i], edge_vertex_ids2[i])
+
+        # Configure mesh connectivity (triangles from edges) for 2D simulations
+        if self._fenics_dims == 2:
+            edges = get_coupling_triangles(function_space, coupling_subdomain, self._precice_edge_dict)
+            for edges_ids in edges:
+                self._interface.set_mesh_triangle(self._interface.get_mesh_id(self._config.get_coupling_mesh_name()),
+                                                  self._precice_edge_dict[edges_ids[0]],
+                                                  self._precice_edge_dict[edges_ids[1]],
+                                                  self._precice_edge_dict[edges_ids[2]])
+        else:
+            print("Mesh connectivity information is not written for 3D cases.")
 
         precice_dt = self._interface.initialize()
 

tests/unit/test_adapter_core.py

@@ -29,7 +29,7 @@ def inside(self, x, on_boundary):
             if right_edge.inside(v.point(), True):
                 global_ids.append(v.global_index())
 
-        edge_vertex_ids1, edge_vertex_ids2 = get_coupling_boundary_edges(V, right_edge, global_ids, id_mapping)
+        edge_vertex_ids1, edge_vertex_ids2, _ = get_coupling_boundary_edges(V, right_edge, global_ids, id_mapping)
 
         self.assertEqual(len(edge_vertex_ids1), 10)
         self.assertEqual(len(edge_vertex_ids2), 10)