Add as_tensor keyword to get_gradient/flux, fixes #56

Author: bleyerj

bindings/python/mgis/fenics/gradient_flux.py

@@ -11,9 +11,9 @@
 import numpy as np
 import ufl
 from .utils import local_project, symmetric_tensor_to_vector, \
-                nonsymmetric_tensor_to_vector, get_quadrature_element
+                nonsymmetric_tensor_to_vector, get_quadrature_element, \
+                vector_to_tensor
 
-<<<<<<< HEAD
 class QuadratureFunction:
     """An abstract class for functions defined at quadrature points"""
     def __init__(self, name, shape):
@@ -31,7 +31,7 @@ def initialize_function(self, mesh, quadrature_degree):
     def update(self, x):
         self.function.vector().set_local(x)
 
-    def project_on(self, space, degree):
+    def project_on(self, space, degree, as_tensor=False):
         """
         Returns the projection on a standard CG/DG space
 
@@ -42,13 +42,20 @@ def project_on(self, space, degree):
             FunctionSpace type ("CG", "DG",...)
         degree: int
             FunctionSpace degree
+        as_tensor: bool
+            Returned as a tensor if True, in vector notation otherwise
         """
-        if self.shape == 1:
+        fun = self.function
+        if as_tensor:
+            fun = vector_to_tensor(self.function)
+            shape = ufl.shape(fun)
+            V = TensorFunctionSpace(self.mesh, space, degree, shape=shape)
+        elif self.shape == 1:
             V = FunctionSpace(self.mesh, space, degree)
         else:
             V = VectorFunctionSpace(self.mesh, space, degree, dim=self.shape)
         v = Function(V, name=self.name)
-        v.assign(project(self.function, V,
+        v.assign(project(fun, V,
                          form_compiler_parameters={"quadrature_degree": self.quadrature_degree}))
         return v
 

bindings/python/mgis/fenics/nonlinear_problem.py

@@ -390,7 +390,7 @@ def update_constitutive_law(self):
         self.update_tangent_blocks()
         self.update_internal_state_variables()
 
-    def get_state_variable(self, name, project_on=None):
+    def get_state_variable(self, name, project_on=None, as_tensor=False):
         """
         Returns the function associated with an internal state variable
 
@@ -411,9 +411,9 @@ def get_state_variable(self, name, project_on=None):
         if project_on is None:
             return self.state_variables["internal"][name].function
         else:
-            return self.state_variables["internal"][name].project_on(*project_on)
+            return self.state_variables["internal"][name].project_on(*project_on, as_tensor)
 
-    def get_flux(self, name, project_on=None):
+    def get_flux(self, name, project_on=None, as_tensor=False):
         """
         Returns the function associated with a flux
 
@@ -434,9 +434,9 @@ def get_flux(self, name, project_on=None):
         if project_on is None:
             return self.fluxes[name].function
         else:
-            return self.fluxes[name].project_on(*project_on)
+            return self.fluxes[name].project_on(*project_on, as_tensor)
 
-    def get_gradient(self, name, project_on=None):
+    def get_gradient(self, name, project_on=None, as_tensor=False):
         """
         Returns the function associated with a gradient
 
@@ -457,7 +457,7 @@ def get_gradient(self, name, project_on=None):
         if project_on is None:
             return self.gradients[name].function
         else:
-            return self.gradients[name].project_on(*project_on)
+            return self.gradients[name].project_on(*project_on, as_tensor)
 
     def get_dissipated_energy(self):
         """Dissipated energy computed from MFront @DissipatedEnergy"""

bindings/python/mgis/fenics/utils.py

@@ -53,6 +53,25 @@ def nonsymmetric_tensor_to_vector(T, T22=0):
     else:
         raise NotImplementedError
 
+def vector_to_tensor(T):
+    """ Return vector following MFront conventions as a tensor """
+    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]]])
+    elif ufl.shape(T)==(5,):
+        return as_tensor([[T[0], T[3]],
+                          T[4], T[1]])
+    elif ufl.shape(T)==(9,):
+        return as_tensor([[T[0], T[3], T[5]],
+                          [T[4], T[1], T[7]],
+                          [T[6], T[8], T[2]]])
+    else:
+        raise NotImplementedError
+
 def axi_grad(r, v):
     """
         Axisymmetric gradient in cylindrical coordinate (er, etheta, ez) for: