Implementing proper treatment of vector valued data

fenicsadapter/fenicsadapter.py

@@ -9,9 +9,10 @@
 from scipy.interpolate import interp1d
 import numpy as np
 from .config import Config
+from enum import Enum
 
 try:
-    import PySolverInterface
+    import precice
 except ImportError:
     import os
     import sys
@@ -22,8 +23,12 @@
     precice_root = os.getenv('PRECICE_ROOT')
     precice_python_adapter_root = precice_root+"/src/precice/bindings/python"
     sys.path.insert(0, precice_python_adapter_root)
-    import PySolverInterface
+    import precice
 
+class FunctionType(Enum):
+    """ Defines scalar- and vector-valued function """
+    SCALAR = 0 # scalar valued funtion
+    VECTOR = 1 # vector valued function
 
 class CustomExpression(UserExpression):
     """Creates functional representation (for FEniCS) of nodal data
@@ -41,26 +46,62 @@ def update_boundary_data(self, vals, coords_x, coords_y=None, coords_z=None):
             coords_z = np.zeros(self._coords_x.shape)
         self._coords_z = coords_z
 
-        self._vals = vals.flatten()
-        assert (self._vals.shape == self._coords_x.shape)
+        self._vals = vals
 
-    def rbf_interpol(self, x):
-        if x.__len__() == 1:
+        if self.isScalarValues():
+            assert (self._vals.shape == self._coords_x.shape)
+        elif self.isVectorValues():
+            assert (self._vals.shape[0] == self._coords_x.shape[0])
+
+    def rbf_interpol(self, x, dim_no):
+        if x.__len__() == 1: # for 1D only R->R mapping is allowed by preCICE, no need to implement Vector case
             f = Rbf(self._coords_x, self._vals.flatten())
             return f(x)
         if x.__len__() == 2:
-            f = Rbf(self._coords_x, self._coords_y, self._vals.flatten())
+            if self.isScalarValues(): #check if scalar or vector-valued
+                f = Rbf(self._coords_x, self._coords_y, self._vals.flatten())
+            elif self.isVectorValues():
+                f = Rbf(self._coords_x, self._coords_y, self._vals[:,dim_no].flatten()) # extract dim_no element of each vector
+            else: # TODO: this is already checked in isScalarValues()!
+                raise Exception("Dimension of the function is 0 or negative!")
             return f(x[0], x[1])
-        if x.__len__() == 3:
-            f = Rbf(self._coords_x, self._coords_y, self._coords_z, self._vals.flatten())
+        if x.__len__() == 3: # this case has not been tested yet
+            if self.isScalarValues():
+                f = Rbf(self._coords_x, self._coords_y, self._coords_z, self._vals.flatten())
+            elif self.isVectorValues():
+                f = Rbf(self._coords_x, self._coords_y, self._coords_z, self._vals[:,dim_no].flatten())
+            else: # TODO: this is already checked in isScalarValues()!
+                raise Exception("Dimension of the function is 0 or negative!")
             return f(x[0], x[1], x[2])
 
     def lin_interpol(self, x):
         f = interp1d(self._coords_x, self._vals)
         return f(x.x())
 
-    def eval(self, value, x):
-        value[0] = self.rbf_interpol(x)
+    def eval(self, value, x): # overloaded function
+        """
+        Overrides UserExpression.eval(). Called by Expression(x_coord). handles
+        scalar- and vector-valued functions evaluation.
+        """
+        for i in range(self._vals.ndim):
+            value[i] = self.rbf_interpol(x,i)
+
+    def isScalarValues(self):
+        """ Determines if function being interpolated is scalar-valued based on
+        dimension of provided vals vector
+        """
+        if self._vals.ndim == 1:
+            return True
+        elif self._vals.ndim > 1:
+            return False
+        else:
+            raise Exception("Dimension of the function is 0 or negative!")
+
+    def isVectorValues(self):
+        """ Determines if function being interpolated is vector-valued based on
+        dimension of provided vals vector
+        """
+        return not self.isScalarValues()
 
 
 class Adapter:
@@ -75,9 +116,9 @@ def __init__(self, adapter_config_filename='precice-adapter-config.json'):
 
         self._solver_name = self._config.get_solver_name()
 
-        self._interface = PySolverInterface.PySolverInterface(self._solver_name, 0, 1)
+        self._interface = precice.Interface(self._solver_name, 0, 1)
         self._interface.configure(self._config.get_config_file_name())
-        self._dimensions = self._interface.getDimensions()
+        self._dimensions = self._interface.get_dimensions()
 
         self._coupling_subdomain = None # initialized later
         self._mesh_fenics = None # initialized later
@@ -86,18 +127,18 @@ def __init__(self, adapter_config_filename='precice-adapter-config.json'):
         ## coupling mesh related quantities
         self._coupling_mesh_vertices = None # initialized later
         self._mesh_name = self._config.get_coupling_mesh_name()
-        self._mesh_id = self._interface.getMeshID(self._mesh_name)
+        self._mesh_id = self._interface.get_mesh_id(self._mesh_name)
         self._vertex_ids = None # initialized later
         self._n_vertices = None # initialized later
 
         ## write data related quantities (write data is written by this solver to preCICE)
         self._write_data_name = self._config.get_write_data_name()
-        self._write_data_id = self._interface.getDataID(self._write_data_name, self._mesh_id)
+        self._write_data_id = self._interface.get_data_id(self._write_data_name, self._mesh_id)
         self._write_data = None
 
         ## read data related quantities (read data is read by this solver from preCICE)
         self._read_data_name = self._config.get_read_data_name()
-        self._read_data_id = self._interface.getDataID(self._read_data_name, self._mesh_id)
+        self._read_data_id = self._interface.get_data_id(self._read_data_name, self._mesh_id)
         self._read_data = None
 
         ## numerics
@@ -108,6 +149,9 @@ def __init__(self, adapter_config_filename='precice-adapter-config.json'):
         self._t_cp = None  # time of the checkpoint
         self._n_cp = None  # timestep of the checkpoint
 
+        #function space
+        self._function_space = None
+
     def convert_fenics_to_precice(self, data, mesh, subdomain):
         """Converts FEniCS data of type dolfin.Function into
         Numpy array for all x and y coordinates on the boundary.
@@ -160,7 +204,7 @@ def set_coupling_mesh(self, mesh, subdomain):
         self._mesh_fenics = mesh
         self._coupling_mesh_vertices, self._n_vertices = self.extract_coupling_boundary_vertices()
         self._vertex_ids = np.zeros(self._n_vertices)
-        self._interface.setMeshVertices(self._mesh_id, self._n_vertices, self._coupling_mesh_vertices.flatten('F'), self._vertex_ids)
+        self._interface.set_mesh_vertices(self._mesh_id, self._n_vertices, self._coupling_mesh_vertices.flatten('F'), self._vertex_ids)
 
     def set_write_field(self, write_function_init):
         """Sets the write field. Called by initalize() function at the
@@ -181,11 +225,10 @@ def set_read_field(self, read_function_init):
     def create_coupling_boundary_condition(self):
         """Creates the coupling boundary conditions using CustomExpression."""
         x_vert, y_vert = self.extract_coupling_boundary_coordinates()
-
         try:  # works with dolfin 1.6.0
-            self._coupling_bc_expression = CustomExpression()
+            self._coupling_bc_expression = CustomExpression(element=self._function_space.ufl_element()) # elemnt information must be provided, else DOLFIN assumes scalar function
         except (TypeError, KeyError):  # works with dolfin 2017.2.0
-            self._coupling_bc_expression = CustomExpression(degree=0)
+            self._coupling_bc_expression = CustomExpression(element=self._function_space.ufl_element(),degree=0)
         self._coupling_bc_expression.set_boundary_data(self._read_data, x_vert, y_vert)
 
     def create_coupling_dirichlet_boundary_condition(self, function_space):
@@ -194,8 +237,9 @@ def create_coupling_dirichlet_boundary_condition(self, function_space):
 
         :return: dolfin.DirichletBC()
         """
+        self._function_space = function_space
         self.create_coupling_boundary_condition()
-        return dolfin.DirichletBC(function_space, self._coupling_bc_expression, self._coupling_subdomain)
+        return dolfin.DirichletBC(self._function_space, self._coupling_bc_expression, self._coupling_subdomain)
 
     def create_coupling_neumann_boundary_condition(self, test_functions):
         """Creates the coupling Neumann boundary conditions using
@@ -205,59 +249,60 @@ def create_coupling_neumann_boundary_condition(self, test_functions):
          Langtangen, Hans Petter, and Anders Logg. "Solving PDEs in Python The
          FEniCS Tutorial Volume I." (2016).)
         """
+        self._function_space = test_functions.function_space()
         self.create_coupling_boundary_condition()
         return self._coupling_bc_expression * test_functions * 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).)
 
     def advance(self, write_function, u_np1, u_n, t, dt, n):
-        """Calls preCICE advance function using PySolverInterface and manages checkpointing.
+        """Calls preCICE advance function using precice and manages checkpointing.
         The solution u_n is updated by this function via call-by-reference. The corresponding values for t and n are returned.
-        
+
         This means:
         * either, the checkpoint self._u_cp is assigned to u_n to repeat the iteration,
         * or u_n+1 is assigned to u_n and the checkpoint is updated correspondingly.
-        
+
         :param write_function: a FEniCS function being sent to the other participant as boundary condition at the coupling interface
         :param u_np1: new value of FEniCS solution u_n+1 at time t_n+1 = t+dt
         :param u_n: old value of FEniCS solution u_n at time t_n = t; updated via call-by-reference
         :param t: current time t_n for timestep n
         :param dt: timestep size dt = t_n+1 - t_n
         :param n: current timestep
-        :return: return starting time t and timestep n for next FEniCS solver iteration. u_n is updated by advance correspondingly. 
+        :return: return starting time t and timestep n for next FEniCS solver iteration. u_n is updated by advance correspondingly.
         """
 
         # sample write data at interface
         x_vert, y_vert = self.extract_coupling_boundary_coordinates()
         self._write_data = self.convert_fenics_to_precice(write_function, self._mesh_fenics, self._coupling_subdomain)
 
         # communication
-        self._interface.writeBlockScalarData(self._write_data_id, self._n_vertices, self._vertex_ids, self._write_data)
+        self._interface.write_block_scalar_data(self._write_data_id, self._n_vertices, self._vertex_ids, self._write_data)
         max_dt = self._interface.advance(dt)
-        self._interface.readBlockScalarData(self._read_data_id, self._n_vertices, self._vertex_ids, self._read_data)
+        self._interface.read_block_scalar_data(self._read_data_id, self._n_vertices, self._vertex_ids, self._read_data)
 
         # update boundary condition with read data
         self._coupling_bc_expression.update_boundary_data(self._read_data, x_vert, y_vert)
 
         precice_step_complete = False
 
         # checkpointing
-        if self._interface.isActionRequired(PySolverInterface.PyActionReadIterationCheckpoint()):
+        if self._interface.is_action_required(precice.action_read_iteration_checkpoint()):
             # continue FEniCS computation from checkpoint
             u_n.assign(self._u_cp)  # set u_n to value of checkpoint
             t = self._t_cp
             n = self._n_cp
-            self._interface.fulfilledAction(PySolverInterface.PyActionReadIterationCheckpoint())
+            self._interface.fulfilled_action(precice.action_read_iteration_checkpoint())
         else:
             u_n.assign(u_np1)
             t = new_t = t + dt  # todo the variables new_t, new_n could be saved, by just using t and n below, however I think it improved readability.
             n = new_n = n + 1
 
-        if self._interface.isActionRequired(PySolverInterface.PyActionWriteIterationCheckpoint()):
+        if self._interface.is_action_required(precice.action_write_iteration_checkpoint()):
             # continue FEniCS computation with u_np1
             # update checkpoint
             self._u_cp.assign(u_np1)
             self._t_cp = new_t
             self._n_cp = new_n
-            self._interface.fulfilledAction(PySolverInterface.PyActionWriteIterationCheckpoint())
+            self._interface.fulfilled_action(precice.action_write_iteration_checkpoint())
             precice_step_complete = True
 
         return t, n, precice_step_complete, max_dt
@@ -273,20 +318,31 @@ def initialize(self, coupling_subdomain, mesh, read_field, write_field, u_n, t=0
         self.set_write_field(write_field)
         self._precice_tau = self._interface.initialize()
 
-        if self._interface.isActionRequired(PySolverInterface.PyActionWriteInitialData()):
-            self._interface.writeBlockScalarData(self._write_data_id, self._n_vertices, self._vertex_ids, self._write_data)
-            self._interface.fulfilledAction(PySolverInterface.PyActionWriteInitialData())
+        if self._interface.is_action_required(precice.action_write_initial_data()):
+            if  self.function_type(write_field) is FunctionType.SCALAR: #write_field.value_rank() == 0:
+                self._interface.write_block_scalar_data(self._write_data_id, self._n_vertices, self._vertex_ids, self._write_data)
+            elif self.function_type(write_field) is FunctionType.VECTOR:
+                self._interface.write_block_vector_data(self._write_data_id, self._n_vertices, self._vertex_ids, self._write_data.ravel())
+            else:
+                raise Exception("Rank of function space is neither 0 nor 1")
+            self._interface.fulfilled_action(precice.action_write_initial_data())
 
-        self._interface.initializeData()
+        self._interface.initialize_data()
 
-        if self._interface.isReadDataAvailable():
-            self._interface.readBlockScalarData(self._read_data_id, self._n_vertices, self._vertex_ids, self._read_data)
+        if self._interface.is_read_data_available():
+            if self.function_type(read_field) is FunctionType.SCALAR:
+                self._interface.read_block_scalar_data(self._read_data_id, self._n_vertices, self._vertex_ids, self._read_data)
+            elif self.function_type(read_field) is FunctionType.VECTOR:
+                self._interface.read_block_vector_data(self._read_data_id, self._n_vertices, self._vertex_ids, self._read_data.ravel())
+            else:
+                raise Exception("Rank of function space is neither 0 nor 1")
 
-        if self._interface.isActionRequired(PySolverInterface.PyActionWriteIterationCheckpoint()):
+
+        if self._interface.is_action_required(precice.action_write_iteration_checkpoint()):
             self._u_cp = u_n.copy(deepcopy=True)
             self._t_cp = t
             self._n_cp = n
-            self._interface.fulfilledAction(PySolverInterface.PyActionWriteIterationCheckpoint())
+            self._interface.fulfilled_action(precice.action_write_iteration_checkpoint())
 
         return self._precice_tau
 
@@ -296,7 +352,7 @@ def is_coupling_ongoing(self):
 
         :return: True if the coupling is ongoing, False otherwise
         """
-        return self._interface.isCouplingOngoing()
+        return self._interface.is_coupling_ongoing()
 
     def extract_coupling_boundary_coordinates(self):
         """Extracts the coordinates of vertices that lay on the boundary. 3D
@@ -316,6 +372,18 @@ def extract_coupling_boundary_coordinates(self):
             # todo this has to be fixed for "proper" 3D coupling. Currently this is a workaround for the coupling of 2D fenics with pseudo 3D openfoam
             return vertices_x, vertices_y
 
+    def function_type(self, function):
+        """ Determines if the function is scalar- or vector-valued based on
+        rank evaluation.
+        """
+        # TODO: is is better to use FunctionType.SCALAR.value here ?
+        if function.value_rank() == 0: # scalar-valued functions have rank 0 is FEniCS
+            return FunctionType.SCALAR
+        elif function.value_rank() == 1: # vector-valued functions have rank 1 in FEniCS
+            return FunctionType.VECTOR
+        else:
+            raise Exception("Error determining function type")
+
     def finalize(self):
         """Finalizes the coupling interface."""
         self._interface.finalize()