Fixed the call to function spaces on functions with ufl_function_space in the elasticity module.

src/bvpy/bvp.py

@@ -420,6 +420,7 @@ class BVP(object):
                 logger.warn(str(key) + ' is not in the list of parameter')
 
     def solve(self, **kwargs):
+        print(f' trial function in the solve method: {type(self.trialFunction)}')
         lhs, rhs = self.vform.get_form(self.trialFunction,
                                        self.testFunction,
                                        self.solution)

src/bvpy/vforms/elasticity.py

@@ -21,8 +21,7 @@
 # -----------------------------------------------------------------------
 
 # import fenics as fe
-
-from ufl import nabla_grad, nabla_div, ln
+from typing import Optional
 
 from bvpy.vforms.abstract import AbstractVform, Element
 from bvpy.utils.tensor import applyElementwise
@@ -30,6 +29,7 @@ from bvpy.domains.geometry import normal
 
 import ufl
 from dolfinx import fem
+from ufl import grad, sym, tr
 
 __all__ = ['ElasticForm', 'LinearElasticForm', 'HyperElasticForm']
 
@@ -133,7 +133,7 @@ class ElasticForm(AbstractVform):
 
         """
         tnsr_func_space =\
-            fem.TensorFunctionSpace(displacement.function_space().mesh(),
+            fem.TensorFunctionSpace(displacement.ufl_function_space().mesh,
                                    'DG', 0)
 
         return fem.project(self._strain(displacement), tnsr_func_space)
@@ -154,7 +154,7 @@ class ElasticForm(AbstractVform):
 
         """
         tnsr_func_space =\
-            fem.TensorFunctionSpace(displacement.function_space().mesh(),
+            fem.TensorFunctionSpace(displacement.ufl_function_space().mesh,
                                    'DG', 0)
 
         return fem.project(self._stress(displacement), tnsr_func_space)
@@ -164,47 +164,47 @@ class LinearElasticForm(ElasticForm):
     """Linear elastic vform in the small deformation approximation.
     """
 
-    def _strain(self, u=None):
+    def _strain(self, u:Optional[fem.Function]=None) -> fem.Function:
         """Computes the linear strain defined as the symmetric part of the
         deformation gradient.
 
         Parameters
         ----------
-        u : :class:`dolfin.cpp.function.Function<dolfin.cpp.function.Function>`
-            Displacement field.
+        u :
+            The displacement field to consider.
 
         Returns
         -------
-        :class:`dolfin.cpp.function.Function<dolfin.cpp.function.Function>`
-            The computed strain.
+        The computed strain.
 
         """
-        return 0.5 * (nabla_grad(u) + nabla_grad(u).T)
+        # return 0.5 * (nabla_grad(u) + nabla_grad(u).T)
+        return .5 * sym(grad(u))
 
-    def _stress(self, u=None):
+    def _stress(self, u:Optional[fem.Function]=None) -> fem.Function:
         """Computes the stress as the energetic conjugate of the linear strain
         tensor field.
 
         Parameters
         ----------
-        u : :class:`dolfin.cpp.function.Function<dolfin.cpp.function.Function>`
-            Displacement field.
+        u : 
+            The displacement field to consider.
 
         Returns
         -------
-        :class:`dolfin.cpp.function.Function<dolfin.cpp.function.Function>`
-            The computed stress.
+        The computed stress.
 
         """
         d = u.geometric_dimension()
-
-        fe_stress = self._parameters['lmbda'] * nabla_div(u) * ufl.Identity(d) +\
-                    2 * self._parameters['mu'] * self._strain(u)
+        lmbd = self._parameters['lmbda']
+        mu = self._parameters['mu']
+        fe_stress = lmbd * tr(self._strain(u)) * ufl.Identity(d) +\
+                    2 * mu * self._strain(u)
 
         if self._plane_stress and d==3:
             # Use instead FacetNormal from ufl
             # n = ufl.FacetNormal(mesh)
-            nrl = normal(u.function_space().mesh())
+            nrl = normal(u.ufl_function_space().mesh)
             tgt_proj = ufl.Identity(d) - ufl.outer(nrl, nrl)
 
             return ufl.dot(tgt_proj, ufl.dot(fe_stress, tgt_proj))
@@ -217,7 +217,7 @@ class LinearElasticForm(ElasticForm):
                    * ufl.inner(self._stress(u), self._strain(v)) * self.dx
 
         self.rhs = ufl.dot(self._parameters['source'], v)\
-            * self.dx(domain=v.function_space().mesh())
+            * self.dx(domain=v.ufl_function_space().mesh)
 
     def set_expression(self):
         self._expression = "thickness*stress(u)*strain(v)*dx "