Compatibility update for v3.0.0

fenicsprecice/adapter_core.py

@@ -383,29 +383,29 @@ def edge_is_on(subdomain, this_edge):
         """
         Check whether edge lies within subdomain
         """
-        assert(len(list(vertices(this_edge))) == 2)
+        assert (len(list(vertices(this_edge))) == 2)
         return all([subdomain.inside(v.point(), True) for v in vertices(this_edge)])
 
     vertices1_ids = []
     vertices2_ids = []
-    edges_ids = []
+    fenics_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())
+                fenics_edges_ids.append(edge.index())
 
     vertices1_ids = np.array(vertices1_ids)
     vertices2_ids = np.array(vertices2_ids)
-    edges_ids = np.array(edges_ids)
+    fenics_edges_ids = np.array(fenics_edges_ids)
 
-    return vertices1_ids, vertices2_ids, edges_ids
+    return vertices1_ids, vertices2_ids, fenics_edges_ids
 
 
-def get_coupling_triangles(function_space, coupling_subdomain, precice_edge_dict):
+def get_coupling_triangles(function_space, coupling_subdomain, fenics_edge_ids, id_mapping):
     """
     Extracts triangles of mesh which lie on the coupling region.
 
@@ -415,31 +415,36 @@ def get_coupling_triangles(function_space, coupling_subdomain, precice_edge_dict
         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
+    fenics_edge_ids: numpy array
+        Array with FEniCS IDs of coupling mesh edges
 
     Returns
     -------
-    edges : numpy array
-        Array of edges indices (3 per triangle)
+    vertex_ids : numpy array
+        Array of indices of vertices which make up triangles (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"
+        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 = []
-
+    vertex_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)
+            if all(edge_ids in fenics_edge_ids for edge_ids in [e1.index(), e2.index(), e3.index()]):
+                v1, v2 = vertices(e1)
+                _, v3 = vertices(e2)
+                assert (v3 != v1)
+                assert (v2 != v2)
+                vertex_ids.append([id_mapping[v1.global_index()],
+                                   id_mapping[v2.global_index()],
+                                   id_mapping[v3.global_index()]])
+
+    return np.array(vertex_ids)
 
 
 def get_forces_as_point_sources(fixed_boundary, function_space, data):

fenicsprecice/expression_core.py

@@ -49,7 +49,7 @@ def update_boundary_data(self, vals, coords):
         self._vals = vals
         _, self._dimension = coords.shape
 
-        assert(self._dimension == 2), "Coordinates are of incorrect dimensions"
+        assert (self._dimension == 2), "Coordinates are of incorrect dimensions"
 
         self._coords_x = coords[:, 0]
         self._coords_y = coords[:, 1]
@@ -128,10 +128,10 @@ def is_scalar_valued(self):
         """
         try:
             if self._vals.ndim == 1:
-                assert(self._function_type is FunctionType.SCALAR)
+                assert (self._function_type is FunctionType.SCALAR)
                 return True
             elif self._vals.ndim > 1:
-                assert(self._function_type is FunctionType.VECTOR)
+                assert (self._function_type is FunctionType.VECTOR)
                 return False
             else:
                 raise Exception("Dimension of the function is 0 or negative!")
@@ -149,10 +149,10 @@ def is_vector_valued(self):
         """
         try:
             if self._vals.ndim > 1:
-                assert(self._function_type is FunctionType.VECTOR)
+                assert (self._function_type is FunctionType.VECTOR)
                 return True
             elif self._vals.ndim == 1:
-                assert(self._function_type is FunctionType.SCALAR)
+                assert (self._function_type is FunctionType.SCALAR)
                 return False
             else:
                 raise Exception("Dimension of the function is 0 or negative!")
@@ -170,7 +170,7 @@ class SegregatedRBFInterpolationExpression(CouplingExpression):
     """
 
     def segregated_interpolant_2d(self, coords_x, coords_y, data):
-        assert(coords_x.shape == coords_y.shape)
+        assert (coords_x.shape == coords_y.shape)
         # create least squares system to approximate a * x ** 2 + b * x + c ~= y
 
         def lstsq_interp(x, y, w): return w[0] * x ** 2 + w[1] * y ** 2 + w[2] * x * y + w[3] * x + w[4] * y + w[5]
@@ -242,7 +242,7 @@ def eval(self, value, x):
         :param x: coordinate where expression has to be evaluated
         :param value: buffer where result has to be returned to
         """
-        assert(MPI.COMM_WORLD.Get_size() > 1)
+        assert (MPI.COMM_WORLD.Get_size() > 1)
         for i in range(self._vals.ndim):
             value[i] = 0