Fixes issues with external state variables

Author: bleyerj

bindings/python/mgis/fenics/gradient_flux.py

@@ -77,6 +77,7 @@ def __init__(self, variable, expression, name, symmetric=None):
         self.variable = variable
         if symmetric is None:
             self.expression = expression
+        # TODO: treat axisymmetric case
         elif symmetric:
             if ufl.shape(expression) == (2, 2):
                 self.expression = as_vector([symmetric_tensor_to_vector(expression)[i] for i in range(4)])
@@ -124,8 +125,8 @@ def _evaluate_at_quadrature_points(self, x):
 
 class Var(Gradient):
     """ A simple variable """
-    # def __init__(self, variable, name):
-    #     return Gradient.__init__(self, variable, variable, name)
+    def __init__(self, variable, expression, name):
+        Gradient.__init__(self, variable, expression, name)
 
     def _evaluate_at_quadrature_points(self, x):
         local_project(x, self.function_space, self.dx, self.function)

bindings/python/mgis/fenics/nonlinear_material.py

@@ -9,6 +9,7 @@
 """
 import mgis.behaviour as mgis_bv
 from .gradient_flux import Var
+from .utils import compute_on_quadrature
 import dolfin
 import subprocess
 import os
@@ -24,7 +25,8 @@ class MFrontNonlinearMaterial:
     This class handles the loading of a MFront behaviour through MGIS.
     """
     def __init__(self, path, name, hypothesis="3d",
-                 material_properties={}, parameters={}):
+                 material_properties={}, parameters={},
+                 rotation_matrix=None):
         """
         Parameters
         -----------
@@ -42,12 +44,17 @@ def __init__(self, path, name, hypothesis="3d",
         parameters : dict
             a dictionary of parameters. The dictionary keys must match the parameter
             names declared in the MFront behaviour. Values must be constants.
+        rotation_matrix : Numpy array, list of list, UFL matrix
+            a 3D rotation matrix expressing the rotation from the global
+            frame to the material frame. The matrix can be spatially variable
+            (either UFL matrix or function of Tensor type)
         """
         self.path = path
         self.name = name
         # Defining the modelling hypothesis
         self.hypothesis = mgis_hypothesis[hypothesis]
         self.material_properties = material_properties
+        self.rotation_matrix = rotation_matrix
         # Loading the behaviour
         try:
             self.load_behaviour()
@@ -74,44 +81,39 @@ def load_behaviour(self):
         else:
             self.behaviour = mgis_bv.load(self.path, self.name, self.hypothesis)
 
-    def set_data_manager(self, ngauss):
+    def set_data_manager(self, degree, ngauss, mesh):
         # Setting the material data manager
         self.data_manager = mgis_bv.MaterialDataManager(self.behaviour, ngauss)
-        self.update_material_properties()
+        self.update_material_properties(degree, mesh)
 
     def update_parameters(self, parameters):
         for (key, value) in parameters.items():
             self.behaviour.setParameter(key, value)
 
-    def update_material_properties(self, material_properties=None):
+    def update_material_properties(self, degree, mesh, material_properties=None):
         if material_properties is not None:
             self.material_properties = material_properties
         for s in [self.data_manager.s0, self.data_manager.s1]:
             for (key, value) in self.material_properties.items():
                 if type(value) in [int, float]:
                     mgis_bv.setMaterialProperty(s, key, value)
                 else:
-                    if isinstance(value, dolfin.Function):
-                        value = value.vector().get_local()
-                    mgis_bv.setMaterialProperty(s, key, value, mgis_bv.MaterialStateManagerStorageMode.LocalStorage)
+                    values = compute_on_quadrature(value, mesh, degree).vector().get_local()
+                    mgis_bv.setMaterialProperty(s, key, values, mgis_bv.MaterialStateManagerStorageMode.LocalStorage)
 
-    def update_external_state_variables(self, external_state_variables):
+    def update_external_state_variables(self, degree, mesh, external_state_variables):
         s = self.data_manager.s1
         for (key, value) in external_state_variables.items():
             if type(value) in [int, float]:
                 mgis_bv.setExternalStateVariable(s, key, value)
             elif isinstance(value, dolfin.Constant):
                 mgis_bv.setExternalStateVariable(s, key, float(value))
             else:
-                if isinstance(value, dolfin.Function):
-                    values = value.vector().get_local()
-                elif isinstance(value, Var):
+                if isinstance(value, Var):
                     value.update()
                     values = value.function.vector().get_local()
                 else:
-                    print(isinstance(value, Var))
-                    print(value)
-                    raise NotImplementedError("{} type is not supported for external state variables".format(type(value)))
+                    values = compute_on_quadrature(value, mesh, degree).vector().get_local()
                 mgis_bv.setExternalStateVariable(s, key, values, mgis_bv.MaterialStateManagerStorageMode.LocalStorage)
 
     def get_parameter(self, name):

bindings/python/mgis/fenics/nonlinear_problem.py

@@ -61,7 +61,7 @@ def __init__(self, u, material, quadrature_degree=2, bcs=None):
         # compute Gauss points numbers
         self.ngauss = Function(self.W0).vector().local_size()
         # Set data manager
-        self.material.set_data_manager(self.ngauss)
+        self.material.set_data_manager(self.quadrature_degree, self.ngauss, self.mesh)
         # dummy function used for computing global quantity
         self._dummy_function = Function(self.W0)
 
@@ -90,6 +90,14 @@ def __init__(self, u, material, quadrature_degree=2, bcs=None):
 
         self.dt = 0
 
+        if self.material.rotation_matrix is not None:
+            if isinstance(self.material.rotation_matrix, (list, tuple, np.ndarray)):
+                self.rotation_values = np.asarray(self.material.rotation_matrix).ravel()
+            else:
+                self.rotation_values = compute_on_quadrature(self.material.rotation_matrix,
+                                                             self.mesh, self.quadrature_degree)
+                self.rotation_values = self.rotation_values.vector().get_local()
+
         self.state_variables =  {"internal": None,
                                  "external": dict.fromkeys(self.material.get_external_state_variable_names(), None)}
         self.gradients = dict.fromkeys(self.material.get_gradient_names(), None)
@@ -194,12 +202,9 @@ def register_external_state_variable(self, name, expression):
         vtype = self.material.behaviour.external_state_variables[pos].type
         if vtype != mgis_bv.VariableType.Scalar:
             raise NotImplementedError("Only scalar external state variables are handled")
-        if isinstance(expression, Constant):
-            self.state_variables["external"].update({name: expression})
-        elif type(expression) == float:
-            self.state_variables["external"].update({name: Constant(expression)})
-        else:
-            self.state_variables["external"].update({name: Var(self.u, expression, name)})
+        if type(expression) == float:
+            expression = Constant(expression)
+        self.state_variables["external"].update({name: Var(self.u, expression, name)})
 
     def set_loading(self, Fext):
         """
@@ -212,38 +217,20 @@ def set_loading(self, Fext):
         """
         self._Fext = ufl.replace(Fext, {self.u: self.u_})
 
-    def set_quadrature_function_spaces(self):
-        cell = self.mesh.ufl_cell()
-        W0e = get_quadrature_element(cell, self.quadrature_degree)
-        # scalar quadrature space
-        self.W0 = FunctionSpace(self.mesh, W0e)
-        # compute Gauss points numbers
-        self.ngauss = Function(self.W0).vector().local_size()
-        # Set data manager
-        self.material.set_data_manager(self.ngauss)
-
-        # Get strain measure dimension
-        self.strain_dim = ufl.shape(self.strain_measure(self.u))[0]
-        # Define quadrature spaces for stress/strain and tangent matrix
-        Wsige = get_quadrature_element(cell, self.quadrature_degree, self.strain_dim)
-        # stress/strain quadrature space
-        self.Wsig = FunctionSpace(self.mesh, Wsige)
-        Wce = get_quadrature_element(cell, self.quadrature_degree, (self.strain_dim, self.strain_dim))
-        # tangent matrix quadrature space
-        self.WCt = FunctionSpace(self.mesh, Wce)
-
     def initialize_external_state_variables(self):
         for (s, size) in zip(self.material.get_external_state_variable_names(), self.material.get_external_state_variable_sizes()):
             state_var = self.state_variables["external"][s]
-            if isinstance(state_var, Gradient):
-                state_var.initialize_function(self.mesh, self.quadrature_degree)
-                values = state_var.function.vector().get_local()
-                mgis_bv.setExternalStateVariable(self.material.data_manager.s0, s,
-                                                 values, mgis_bv.MaterialStateManagerStorageMode.LocalStorage)
-            elif isinstance(state_var, Constant):
+            if isinstance(state_var, Constant):
                 mgis_bv.setExternalStateVariable(self.material.data_manager.s0, s, float(state_var))
             else:
-                raise ValueError("External state variable '{}' has not been registered.".format(s))
+                if isinstance(state_var, Var):
+                    state_var.initialize_function(self.mesh, self.quadrature_degree)
+                    values = state_var.function.vector().get_local()
+                else:
+                    values = compute_on_quadrature(state_var, self.mesh,
+                                               self.quadrature_degree).vector().get_local()
+                mgis_bv.setExternalStateVariable(self.material.data_manager.s0, s,
+                                                 values, mgis_bv.MaterialStateManagerStorageMode.LocalStorage)
 
     def initialize_gradients(self):
         buff=0
@@ -287,7 +274,7 @@ def initialize_tangent_blocks(self):
                     except:
                         value = self.state_variables["external"].get(t[1], None)
                         if value is not None and isinstance(value, Var):
-                                flux_gradients.append(value)
+                            flux_gradients.append(value)
                         else:
                             raise ValueError("'{}' could not be associated with a registered gradient or state variable.".format(t[1]))
             self.fluxes[f].initialize_tangent_blocks(flux_gradients)
@@ -347,15 +334,23 @@ def update_tangent_blocks(self):
                 except KeyError:
                     raise KeyError("'{}' could not be found as a flux or an internal state variable.")
             block_shape = self.flattened_block_shapes[i]
-            t.vector().set_local(self.material.data_manager.K[:,buff:buff+block_shape].flatten())
+            tang_block_vals = self.material.data_manager.K[:,buff:buff+block_shape].ravel()
+            if self.material.rotation_matrix is not None:
+                mgis_bv.rotateTangentOperatorBlocks(tang_block_vals, self.material.behaviour,
+                                                    self.rotation_values)
+            t.vector().set_local(tang_block_vals)
             buff += block_shape
 
     def update_fluxes(self):
         buff = 0
         for (i, f) in enumerate(self.material.get_flux_names()):
             flux = self.fluxes[f]
             block_shape = self.material.get_flux_sizes()[i]
-            flux.function.vector().set_local(self.material.data_manager.s1.thermodynamic_forces[:,buff:buff+block_shape].flatten())
+            flux_vals = self.material.data_manager.s1.thermodynamic_forces[:,buff:buff+block_shape].ravel()
+            if self.material.rotation_matrix is not None:
+                mgis_bv.rotateThermodynamicForces(flux_vals, self.material.behaviour,
+                                                  self.rotation_values)
+            flux.function.vector().set_local(flux_vals)
             buff += block_shape
 
     def update_gradients(self):
@@ -365,6 +360,9 @@ def update_gradients(self):
             gradient.update()
             block_shape = self.material.get_gradient_sizes()[i]
             grad_vals = gradient.function.vector().get_local()
+            if self.material.rotation_matrix is not None:
+                mgis_bv.rotateGradients(grad_vals, self.material.behaviour,
+                                        self.rotation_values)
             if gradient.shape > 0:
                 grad_vals = grad_vals.reshape((self.material.data_manager.n, gradient.shape))
             else:
@@ -384,7 +382,8 @@ def update_internal_state_variables(self):
     def update_constitutive_law(self):
         """Performs the consitutive law update"""
         self.update_gradients()
-        self.material.update_external_state_variables(self.state_variables["external"])
+        self.material.update_external_state_variables(self.quadrature_degree, self.mesh,
+                                                      self.state_variables["external"])
         # integrate the behaviour
         mgis_bv.integrate(self.material.data_manager, self.integration_type,
                           self.dt, 0, self.material.data_manager.n);

bindings/python/mgis/fenics/utils.py

@@ -118,4 +118,48 @@ def get_quadrature_element(cell, degree, dim=0):
     elif type(dim) == tuple:
         return TensorElement("Quadrature", cell, degree=degree, shape=dim, quad_scheme='default')
     else:
-        raise ValueError("Wrong shape for dim=", dim)
+        raise ValueError("Wrong shape for dim=", dim)
+
+def interpolate_on_quadrature(f, degree):
+    V = f.function_space()
+    Ve = V.ufl_element()
+    if Ve.family()=="Quadrature":
+        return f
+    else:
+        shape = Ve.value_shape()
+        if len(shape) == 1:
+            shape = shape[0]
+        Vge = Ve.__class__("Quadrature", Ve.cell(), degree, shape,
+                           quad_scheme="default")
+        Vg = FunctionSpace(V.mesh(), Vge)
+        fg = Function(Vg)
+        fg.interpolate(f)
+        return fg
+
+def project_on_quadrature(f, mesh, degree):
+    shape = f.ufl_shape
+    Vge = get_quadrature_element(mesh.ufl_cell(), degree, shape)
+    Vg = FunctionSpace(mesh, Vge)
+    dx = Measure("dx", domain=mesh, metadata={"quadrature_scheme": "default",
+                                              "quadrature_degree": degree})
+    fg = Function(Vg)
+    fg.assign(local_project(f, Vg, dx))
+    return fg
+
+def compute_on_quadrature(f, mesh, degree):
+    shape = f.ufl_shape
+    Vge = get_quadrature_element(mesh.ufl_cell(), degree, shape)
+    Vg = FunctionSpace(mesh, Vge)
+    dx = Measure("dx", domain=mesh, metadata={"quadrature_scheme": "default",
+                                              "quadrature_degree": degree})
+    fg = Function(Vg)
+    if (isinstance(f, function.function.Function) or
+        isinstance(f, function.function.Constant) or
+        isinstance(f, function.function.Expression)):
+        fg.interpolate(f)
+        # return interpolate_on_quadrature(f, degree)
+    else:
+        fg.assign(local_project(f, Vg, dx))
+    return fg
+        # fg.interpolate(f)
+        # return project_on_quadrature(f, degree)

docs/web/release-notes-1.2.md

@@ -127,6 +127,10 @@ module:
 
 # Issues solved
 
+## Issue #66: `mgis.fenics`: Missing update of external state variables when seen as UFL expression
+
+Declaring external state variables as UFL expressions were not properly updated between time steps.
+
 ## Issues #54: Inform the calling code about `@DissipatedEnergy` and/or `@InternalEnergy`
 
 The `Behaviour` class now exposes two new boolean data members:
@@ -241,4 +245,4 @@ The previous instruction imports the `mgis::MFrontGenericInterface` target, whic
 
 For details, see <https://github.com/thelfer/MFrontGenericInterfaceSupport/issues/25>
 
-# References
+# References