Merge pull request thelfer#99 from thelfer/tensor-vec-repr

Fixes thelfer#98

Author: bleyerj

bindings/python/mgis/fenics/gradient_flux.py

@@ -17,13 +17,15 @@
     get_quadrature_element,
     vector_to_tensor,
 )
+import mgis.behaviour as mgis_bv
 
 
 class QuadratureFunction:
     """An abstract class for functions defined at quadrature points"""
-    def __init__(self, name, shape):
+    def __init__(self, name, shape, hypothesis):
         self.shape = shape
         self.name = name
+        self.hypothesis = hypothesis
 
     def initialize_function(self, mesh, quadrature_degree):
         self.quadrature_degree = quadrature_degree
@@ -87,26 +89,32 @@ class Gradient(QuadratureFunction):
     This class is intended for internal use only. Gradient objects must be
     declared by the user using the registration concept.
     """
-    def __init__(self, variable, expression, name, symmetric=None):
+    def __init__(self, variable, expression, name, hypothesis, 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)
-                ])
-            else:
-                self.expression = symmetric_tensor_to_vector(expression)
         else:
-            if ufl.shape(expression) == (2, 2):
-                self.expression = as_vector([
-                    nonsymmetric_tensor_to_vector(expression)[i]
-                    for i in range(5)
-                ])
+            if symmetric:
+                converter = symmetric_tensor_to_vector
+            else:
+                converter = nonsymmetric_tensor_to_vector
+            if hypothesis in [
+                    mgis_bv.Hypothesis.PlaneStrain,
+                    mgis_bv.Hypothesis.Axisymmetrical,
+            ]:
+                if ufl.shape(expression) == (2, 2):
+                    T22 = 0
+                    expression_2d = expression
+                elif ufl.shape(expression) == (3, 3):
+                    T22 = expression[2, 2]
+                    expression_2d = ufl.as_tensor(
+                        [[expression[i, j] for j in range(2)]
+                         for i in range(2)])
+                self.expression = converter(expression_2d, T22)
             else:
-                self.expression = nonsymmetric_tensor_to_vector(expression)
+                self.expression = converter(expression)
+
         shape = ufl.shape(self.expression)
         if len(shape) == 1:
             self.shape = shape[0]
@@ -115,6 +123,7 @@ def __init__(self, variable, expression, name, symmetric=None):
         else:
             self.shape = shape
         self.name = name
+        self.hypothesis = hypothesis
 
     def __call__(self, v):
         return ufl.replace(self.expression, {self.variable: v})
@@ -148,8 +157,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, expression, name, hypothesis):
+        Gradient.__init__(self, variable, expression, name, hypothesis)
 
     def _evaluate_at_quadrature_points(self, x):
         local_project(x, self.function_space, self.dx, self.function)

bindings/python/mgis/fenics/nonlinear_problem.py

@@ -25,7 +25,7 @@
 }
 Id_Grad = {
     mgis_bv.Hypothesis.Tridimensional: lambda u: Identity(3) + grad(u),
-    mgis_bv.Hypothesis.PlaneStrain: lambda u: Identity(2) + grad(u),
+    mgis_bv.Hypothesis.PlaneStrain: lambda u: Identity(3) + grad_3d(u),
     mgis_bv.Hypothesis.Axisymmetrical:
     lambda r, u: Identity(3) + axi_grad(r, u),
 }
@@ -222,8 +222,16 @@ def register_gradient(self, name, expression):
             symmetric = True
         else:
             symmetric = None
-        self.gradients.update(
-            {name: Gradient(self.u, expression, name, symmetric=symmetric)})
+        self.gradients.update({
+            name:
+            Gradient(
+                self.u,
+                expression,
+                name,
+                self.material.hypothesis,
+                symmetric=symmetric,
+            )
+        })
 
     def register_external_state_variable(self, name, expression):
         """
@@ -245,7 +253,7 @@ def register_external_state_variable(self, name, expression):
         if type(expression) == float:
             expression = Constant(expression)
         self.state_variables["external"].update(
-            {name: Var(self.u, expression, name)})
+            {name: Var(self.u, expression, name, self.material.hypothesis)})
 
     def set_loading(self, Fext):
         """
@@ -307,7 +315,7 @@ def initialize_fluxes(self):
         fluxes = []
         for (f, size) in zip(self.material.get_flux_names(),
                              self.material.get_flux_sizes()):
-            flux = Flux(f, size)
+            flux = Flux(f, size, self.material.hypothesis)
             flux.initialize_function(self.mesh, self.quadrature_degree)
             fluxes.append(flux)
         self.fluxes = dict(zip(self.material.get_flux_names(), fluxes))
@@ -318,7 +326,8 @@ def initialize_internal_state_variables(self):
                 self.material.get_internal_state_variable_names(),
                 self.material.get_internal_state_variable_sizes(),
         ):
-            state_variable = InternalStateVariable(s, size)
+            state_variable = InternalStateVariable(s, size,
+                                                   self.material.hypothesis)
             state_variable.initialize_function(self.mesh,
                                                self.quadrature_degree)
             state_variables.append(state_variable)

bindings/python/mgis/fenics/utils.py

@@ -32,16 +32,18 @@ def local_project(v, V, dx, u=None):
 
 
 def symmetric_tensor_to_vector(T, T22=0):
-    """ Return symmetric tensor components in vector form notation following MFront conventions
-        T22 can be specified when T is only (2,2)"""
+    """Return symmetric tensor components in vector form notation following MFront conventions
+    T22 can be specified when T is only (2,2)"""
     if ufl.shape(T) == (2, 2):
         return as_vector([T[0, 0], T[1, 1], T22, sqrt(2) * T[0, 1]])
     elif ufl.shape(T) == (3, 3):
         return as_vector([
-            T[0, 0], T[1, 1], T[2, 2],
+            T[0, 0],
+            T[1, 1],
+            T[2, 2],
             sqrt(2) * T[0, 1],
             sqrt(2) * T[0, 2],
-            sqrt(2) * T[1, 2]
+            sqrt(2) * T[1, 2],
         ])
     elif len(ufl.shape(T)) == 1:
         return T
@@ -50,14 +52,21 @@ def symmetric_tensor_to_vector(T, T22=0):
 
 
 def nonsymmetric_tensor_to_vector(T, T22=0):
-    """ Return nonsymmetric tensor components in vector form notation following MFront conventions
-        T22 can be specified when T is only (2,2) """
+    """Return nonsymmetric tensor components in vector form notation following MFront conventions
+    T22 can be specified when T is only (2,2)"""
     if ufl.shape(T) == (2, 2):
         return as_vector([T[0, 0], T[1, 1], T22, T[0, 1], T[1, 0]])
     elif ufl.shape(T) == (3, 3):
         return as_vector([
-            T[0, 0], T[1, 1], T[2, 2], T[0, 1], T[1, 0], T[0, 2], T[2, 0],
-            T[1, 2], T[2, 1]
+            T[0, 0],
+            T[1, 1],
+            T[2, 2],
+            T[0, 1],
+            T[1, 0],
+            T[0, 2],
+            T[2, 0],
+            T[1, 2],
+            T[2, 1],
         ])
     elif len(ufl.shape(T)) == 1:
         return T
@@ -70,9 +79,11 @@ def vector_to_tensor(T):
     if ufl.shape(T) == (4, ):
         return as_tensor([[T[0], T[3] / sqrt(2)], T[3] / sqrt(2), T[1]])
     elif ufl.shape(T) == (6, ):
-        return as_tensor([[T[0], T[3] / sqrt(2), T[4] / sqrt(2)],
-                          [T[3] / sqrt(2), T[1], T[5] / sqrt(2)],
-                          [T[4] / sqrt(2), T[5] / sqrt(2), T[2]]])
+        return as_tensor([
+            [T[0], T[3] / sqrt(2), T[4] / sqrt(2)],
+            [T[3] / sqrt(2), T[1], T[5] / sqrt(2)],
+            [T[4] / sqrt(2), T[5] / sqrt(2), T[2]],
+        ])
     elif ufl.shape(T) == (5, ):
         return as_tensor([[T[0], T[3]], T[4], T[1]])
     elif ufl.shape(T) == (9, ):
@@ -84,15 +95,17 @@ def vector_to_tensor(T):
 
 def axi_grad(r, v):
     """
-        Axisymmetric gradient in cylindrical coordinate (er, etheta, ez) for:
-        * a scalar v(r, z)
-        * a 2d-vectorial (vr(r,z), vz(r, z))
-        * a 3d-vectorial (vr(r,z), 0, vz(r, z))
+    Axisymmetric gradient in cylindrical coordinate (er, etheta, ez) for:
+    * a scalar v(r, z)
+    * a 2d-vectorial (vr(r,z), vz(r, z))
+    * a 3d-vectorial (vr(r,z), 0, vz(r, z))
     """
     if ufl.shape(v) == (3, ):
-        return as_matrix([[v[0].dx(0), -v[1] / r, v[0].dx(1)],
-                          [v[1].dx(0), v[0] / r, v[1].dx(1)],
-                          [v[2].dx(0), 0, v[2].dx(1)]])
+        return as_matrix([
+            [v[0].dx(0), -v[1] / r, v[0].dx(1)],
+            [v[1].dx(0), v[0] / r, v[1].dx(1)],
+            [v[2].dx(0), 0, v[2].dx(1)],
+        ])
     elif ufl.shape(v) == (2, ):
         return as_matrix([[v[0].dx(0), v[0].dx(1), 0],
                           [v[1].dx(0), v[1].dx(1), 0], [0, 0, v[0] / r]])
@@ -102,6 +115,11 @@ def axi_grad(r, v):
         raise NotImplementedError
 
 
+def grad_3d(u):
+    return as_matrix([[u[0].dx(0), u[0].dx(1), 0], [u[1].dx(0), u[1].dx(1), 0],
+                      [0, 0, 0]])
+
+
 def symmetric_gradient(g):
     """ Return symmetric gradient components in vector form"""
     return symmetric_tensor_to_vector(sym(g))
@@ -122,26 +140,26 @@ def get_quadrature_element(cell, degree, dim=0):
         return FiniteElement("Quadrature",
                              cell,
                              degree=degree,
-                             quad_scheme='default')
+                             quad_scheme="default")
     elif type(dim) == int:
         return VectorElement("Quadrature",
                              cell,
                              degree=degree,
                              dim=dim,
-                             quad_scheme='default')
+                             quad_scheme="default")
     elif len(dim) == 1 or (len(dim) == 2 and 0 in dim):
         d = [dd for dd in list(dim) if dd != 0][0]
         return VectorElement("Quadrature",
                              cell,
                              degree=degree,
                              dim=d,
-                             quad_scheme='default')
+                             quad_scheme="default")
     elif type(dim) == tuple:
         return TensorElement("Quadrature",
                              cell,
                              degree=degree,
                              shape=dim,
-                             quad_scheme='default')
+                             quad_scheme="default")
     else:
         raise ValueError("Wrong shape for dim=", dim)
 
@@ -170,12 +188,14 @@ 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
-                 })
+    dx = Measure(
+        "dx",
+        domain=mesh,
+        metadata={
+            "quadrature_scheme": "default",
+            "quadrature_degree": degree
+        },
+    )
     fg = Function(Vg)
     fg.assign(local_project(f, Vg, dx))
     return fg
@@ -185,12 +205,14 @@ 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
-                 })
+    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)