Merge pull request #68 from thelfer/revert-67-fix-esv-tangent-blocks

Revert "Fix update of state variables as UFL expression"

Author: bleyerj

CMakeLists.txt

@@ -3,9 +3,6 @@ cmake_minimum_required(VERSION 3.0)
 include(cmake/modules/mgis.cmake)
 project("mfront-generic-interface")
 
-set(CMAKE_CXX_STANDARD 17)
-set(CXX_STANDARD_REQUIRED ON)
-
 # portable-build
 option(enable-portable-build "produce binary that can be shared between various machine (same architecture, same gcc version, different processors" OFF)
 

bindings/python/mgis/fenics/gradient_flux.py

@@ -77,7 +77,6 @@ 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)])
@@ -125,8 +124,8 @@ def _evaluate_at_quadrature_points(self, x):
 
 class Var(Gradient):
     """ A simple variable """
-    def __init__(self, variable, expression, name):
-        Gradient.__init__(self, variable, expression, name)
+    # def __init__(self, variable, name):
+    #     return Gradient.__init__(self, variable, variable, 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,7 +9,6 @@
 """
 import mgis.behaviour as mgis_bv
 from .gradient_flux import Var
-from .utils import compute_on_quadrature
 import dolfin
 import subprocess
 import os
@@ -25,8 +24,7 @@ class MFrontNonlinearMaterial:
     This class handles the loading of a MFront behaviour through MGIS.
     """
     def __init__(self, path, name, hypothesis="3d",
-                 material_properties={}, parameters={},
-                 rotation_matrix=None):
+                 material_properties={}, parameters={}):
         """
         Parameters
         -----------
@@ -44,17 +42,12 @@ 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()
@@ -81,39 +74,44 @@ def load_behaviour(self):
         else:
             self.behaviour = mgis_bv.load(self.path, self.name, self.hypothesis)
 
-    def set_data_manager(self, degree, ngauss, mesh):
+    def set_data_manager(self, ngauss):
         # Setting the material data manager
         self.data_manager = mgis_bv.MaterialDataManager(self.behaviour, ngauss)
-        self.update_material_properties(degree, mesh)
+        self.update_material_properties()
 
     def update_parameters(self, parameters):
         for (key, value) in parameters.items():
             self.behaviour.setParameter(key, value)
 
-    def update_material_properties(self, degree, mesh, material_properties=None):
+    def update_material_properties(self, 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:
-                    values = compute_on_quadrature(value, mesh, degree).vector().get_local()
-                    mgis_bv.setMaterialProperty(s, key, values, mgis_bv.MaterialStateManagerStorageMode.LocalStorage)
+                    if isinstance(value, dolfin.Function):
+                        value = value.vector().get_local()
+                    mgis_bv.setMaterialProperty(s, key, value, mgis_bv.MaterialStateManagerStorageMode.LocalStorage)
 
-    def update_external_state_variables(self, degree, mesh, external_state_variables):
+    def update_external_state_variables(self, 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, Var):
+                if isinstance(value, dolfin.Function):
+                    values = value.vector().get_local()
+                elif isinstance(value, Var):
                     value.update()
                     values = value.function.vector().get_local()
                 else:
-                    values = compute_on_quadrature(value, mesh, degree).vector().get_local()
+                    print(isinstance(value, Var))
+                    print(value)
+                    raise NotImplementedError("{} type is not supported for external state variables".format(type(value)))
                 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.quadrature_degree, self.ngauss, self.mesh)
+        self.material.set_data_manager(self.ngauss)
         # dummy function used for computing global quantity
         self._dummy_function = Function(self.W0)
 
@@ -90,14 +90,6 @@ 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)
@@ -202,9 +194,12 @@ 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 type(expression) == float:
-            expression = Constant(expression)
-        self.state_variables["external"].update({name: Var(self.u, expression, name)})
+        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)})
 
     def set_loading(self, Fext):
         """
@@ -217,20 +212,38 @@ 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, Constant):
-                mgis_bv.setExternalStateVariable(self.material.data_manager.s0, s, float(state_var))
-            else:
-                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()
+            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):
+                mgis_bv.setExternalStateVariable(self.material.data_manager.s0, s, float(state_var))
+            else:
+                raise ValueError("External state variable '{}' has not been registered.".format(s))
 
     def initialize_gradients(self):
         buff=0
@@ -274,7 +287,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)
@@ -334,23 +347,15 @@ 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]
-            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)
+            t.vector().set_local(self.material.data_manager.K[:,buff:buff+block_shape].flatten())
             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_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)
+            flux.function.vector().set_local(self.material.data_manager.s1.thermodynamic_forces[:,buff:buff+block_shape].flatten())
             buff += block_shape
 
     def update_gradients(self):
@@ -360,9 +365,6 @@ 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:
@@ -382,8 +384,7 @@ 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.quadrature_degree, self.mesh,
-                                                      self.state_variables["external"])
+        self.material.update_external_state_variables(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,48 +118,4 @@ 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)
-
-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)
+        raise ValueError("Wrong shape for dim=", dim)