Almost working von Mises plasticity

demos/StationaryHeatTransfer.py

@@ -27,7 +27,7 @@ Tr = 800
 bc = [DirichletBC(V, Constant(Tl), left),
       DirichletBC(V, Constant(Tr), right)]
 
-facets = MeshFunction("size_t", mesh, dim-1)
+facets = MeshFunction("size_t", mesh, 1)
 ds = Measure("ds", subdomain_data=facets)
 
 
@@ -37,17 +37,7 @@ material = mf.MFrontNonlinearMaterial("../materials/src/libBehaviour.so",
                                       hypothesis="plane_strain")
 problem = mf.MFrontNonlinearProblem(T, material)
 problem.bc = bc
-
-j   = Function(Wjs, name="Current heat flux")
-gT  = Function(WgTs, name="Current temperature gradient")
-dj_ddgT = Function(W_dj_ddgT_s, name="Derivative of the Heat Flux with respect to the temperature gradient")
-dj_ddT  = Function(W_dj_ddT_s, name="Derivative of the Heat Flux with respect to the temperature")
-
-T0 = Function(V, name="Temperature at the beginning of the time step")
-T1 = Function(V, name="Temperature estimate at the end of the end step")
-dT = Function(V, name="Iteration correction")
-v = TrialFunction(V)
-T_ = TestFunction(V)
+print(problem.material.data_manager.K.shape)
 
 # definition of the stiffness matrix and residual using standard FEniCS operations
 a_Newton = (inner(grad(v), dot(as_tensor(dj_ddgT), grad(T_)))+
@@ -55,8 +45,8 @@ a_Newton = (inner(grad(v), dot(as_tensor(dj_ddgT), grad(T_)))+
 # a_Newton = inner(grad(v), dot(as_tensor(dj_ddgT), grad(T_)))*dxm
 res = -inner(grad(T_), j)*dxm
 
-problem.a = inner(self.strain_variation(self.du), dot(self.Ct, self.strain_variation(self.u_)))*measure*self.dx
-problam.L = inner(self.strain_variation(self.u_), self.stress)*measure*self.dxl value of the deformation gradient
+#problem.a = inner(self.strain_variation(self.du), dot(self.Ct, self.strain_variation(self.u_)))*measure*self.dx
+#problam.L = inner(self.strain_variation(self.u_), self.stress)*measure*self.dx l value of the deformation gradient
 #  using `MFront` conventions
 T0.assign(Constant(Tl))
 local_project(grad(T0), WgTs, gT)

demos/small_strain_vonMises_plasticity.py

@@ -45,8 +45,9 @@ for hypothesis in ["plane_strain", "axisymmetric"]:
                                           hypothesis=hypothesis,
                                           material_properties=mat_prop)
     problem = mf.MFrontNonlinearProblem(u, material, quadrature_degree=4)
-    problem.set_loading(loading*dot(n, u)*measure*ds(4))
     problem.bc = bc
+    problem.define_form(mf.Stress(u, name="sigma"))
+    problem.set_loading(loading*dot(n, u)*measure*ds(4))
 
     p = problem.get_state_variable(name="EquivalentPlasticStrain")
     assert (ufl.shape(p) == ())

mfront_wrapper/__init__.py

@@ -12,4 +12,5 @@ from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarnin
 warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning)
 
 from .nonlinear_material import MFrontNonlinearMaterial
-from .nonlinear_problem import MFrontNonlinearProblem
\ No newline at end of file
+from .nonlinear_problem import MFrontNonlinearProblem
+from .gradient_flux import *
\ No newline at end of file

mfront_wrapper/gradient_flux.py

@@ -31,6 +31,14 @@ class Gradient:
         """ Directional derivative in direction v """
         return ufl.derivative(self.expression, self.variable, v)
 
+    def initialize_function(self, mesh, quadrature_degree):
+        We = get_quadrature_element(mesh.ufl_cell(), quadrature_degree, self.shape)
+        self.function_space = FunctionSpace(mesh, We)
+        self.function = Function(self.function_space, name=self.name)
+
+    def update(self, x):
+        self.function.vector().set_local(x)
+
 class Grad(Gradient):
     """ Gradient operator : grad(u) """
     def __init__(self, variable, name=None):
@@ -62,7 +70,7 @@ class Var(Gradient):
 
 
 class Flux:
-    def __init__(self, gradients, shape=None, name=None):
+    def __init__(self, gradients, dual_gradient=None, shape=None, name=None):
         if type(gradients)==list:
             self.gradients = gradients
         else:
@@ -74,6 +82,10 @@ class Flux:
             self.shape = self.gradient_shapes[0]
         else:
             self.shape = shape
+        if dual_gradient is None:
+            self.dual_gradient = self.gradients[0]
+        else:
+            self.dual_gradient = dual_gradient
         if name is None:
             self.name = "Flux-{}".format(np.random.randint(1e6))
         else:
@@ -87,6 +99,12 @@ class Flux:
                           get_quadrature_element(mesh.ufl_cell(),
                           quadrature_degree, dim=(self.shape, g.shape))), name="d{}_d{}".format(self.name, g.name))
                           for g in self.gradients]
+        for g in self.gradients:
+            g.initialize_function(mesh, quadrature_degree)
+
+    def update(self, x):
+        self.function.vector().set_local(x)
+
 
 
 class Stress(Flux):
@@ -99,16 +117,16 @@ class ThermalFlux(Flux):
         Flux.__init__(self, [Grad(variable, name="grad(T)"), Var(variable, name="T")],
                              variable.geometric_dimension(), name=name)
 
-mesh = UnitSquareMesh(10, 10)
-V = VectorFunctionSpace(mesh, "CG", 1)
-VT = FunctionSpace(mesh, "CG", 1)
-
-u = Function(V, name="u")
-T = Function(VT, name="Temperature")
-v = Function(V, name="v")
-g = Grad_Id(u, name="F")
-sig = Stress(u)
-j = ThermalFlux(T)
-j.initialize_functions(mesh, 2)
-print([ufl.shape(tang) for tang in j.tangents])
-print([tang.name() for tang in j.tangents])
\ No newline at end of file
+# mesh = UnitSquareMesh(10, 10)
+# V = VectorFunctionSpace(mesh, "CG", 1)
+# VT = FunctionSpace(mesh, "CG", 1)
+
+# u = Function(V, name="u")
+# T = Function(VT, name="Temperature")
+# v = Function(V, name="v")
+# g = Grad_Id(u, name="F")
+# sig = Stress(u)
+# j = ThermalFlux(T)
+# j.initialize_functions(mesh, 2)
+# print([ufl.shape(tang) for tang in j.tangents])
+# print([tang.name() for tang in j.tangents])
\ No newline at end of file

mfront_wrapper/nonlinear_problem.py

 from dolfin import *
 from mfront_wrapper.utils import *
+from mfront_wrapper.gradient_flux import *
 import mgis.behaviour as mgis_bv
 
 class MFrontNonlinearProblem(NonlinearProblem):
@@ -11,37 +12,52 @@ class MFrontNonlinearProblem(NonlinearProblem):
         self.du = TrialFunction(self.V)
         self.mesh = self.V.mesh()
         self.material = material
-#        print(self.material.hypothesis)
         self.axisymmetric = self.material.hypothesis==mgis_bv.Hypothesis.Axisymmetrical
         self.integration_type = mgis_bv.IntegrationType.IntegrationWithConsistentTangentOperator
 
         self.quadrature_degree = quadrature_degree
-        self.set_quadrature_function_spaces()
+#        self.set_quadrature_function_spaces()
 
 
+        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 = self.W0.dim()
+        # Set data manager
+        self.material.set_data_manager(self.ngauss)
+        self.finite_strain = self.material.behaviour.getBehaviourType()=="StandardFiniteStrainBehaviour"
+        if self.material.hypothesis == mgis_bv.Hypothesis.Tridimensional:
+            assert self.u.geometric_dimension()==3, "Conflicting geometric dimension and material hypothesis"
+        else:
+            assert self.u.geometric_dimension()==2, "Conflicting geometric dimension and material hypothesis"
 
         self.bc = []
         self.dx = Measure("dx", metadata={"quadrature_degree": self.quadrature_degree,
                                           "quadrature_scheme": "default"})
         if self.axisymmetric:
             x = SpatialCoordinate(self.mesh)
-            measure = 2*pi*abs(x[0])
+            self.axi_coeff = 2*pi*abs(x[0])
         else:
-            measure = 1
-
-        self.initialize_fields()
-
-        # tangent bilinear form
-        self.a = inner(self.strain_variation(self.du), dot(self.Ct, self.strain_variation(self.u_)))*measure*self.dx
-        # residual form (internal forces)
-        self.L = inner(self.strain_variation(self.u_), self.stress)*measure*self.dx
+            self.axi_coeff = 1
 
         self.solver = NewtonSolver()
 
         self.state_variables = []
 
-    def add_stress(self, gradients=[]):
-        pass
+    def define_form(self, flux):
+        self.flux = flux
+        self.flux.initialize_functions(self.mesh, self.quadrature_degree)
+
+        # tangent bilinear form
+        dg = self.flux.dual_gradient.variation(self.u_)
+        tangent_terms = [inner(dg, dot(t, g.variation(self.du)))
+                         for (t, g) in zip(self.flux.tangents, self.flux.gradients)]
+        self.a = sum(tangent_terms)*self.axi_coeff*self.dx
+        # residual form (internal forces)
+        self.L = inner(dg, self.flux.function)*self.axi_coeff*self.dx
+
 
     def set_loading(self, Fext):
         # adds external forces to residual form
@@ -56,11 +72,7 @@ class MFrontNonlinearProblem(NonlinearProblem):
         self.ngauss = self.W0.dim()
         # Set data manager
         self.material.set_data_manager(self.ngauss)
-        self.finite_strain = self.material.behaviour.getBehaviourType()=="StandardFiniteStrainBehaviour"
-        if self.material.hypothesis == mgis_bv.Hypothesis.Tridimensional:
-            assert self.u.geometric_dimension()==3, "Conflicting geometric dimension and material hypothesis"
-        else:
-            assert self.u.geometric_dimension()==2, "Conflicting geometric dimension and material hypothesis"
+
         # Get strain measure dimension
         self.strain_dim = ufl.shape(self.strain_measure(self.u))[0]
         # Define quadrature spaces for stress/strain and tangent matrix
@@ -111,9 +123,10 @@ class MFrontNonlinearProblem(NonlinearProblem):
             self.Ct.vector().set_local(self.material.data_manager.K.flatten())
 
     def update_constitutive_law(self, u):
-        local_project(self.strain_measure(u), self.Wsig, self.dx, self.strain)
+        g0 = self.flux.gradients[0]
+        local_project(g0(self.u), g0.function_space, self.dx, g0.function)
         # copy the strain values to `MGIS`
-        self.material.data_manager.s1.gradients[:, :] = self.strain.vector().get_local().reshape((self.material.data_manager.n, self.strain_dim))
+        self.material.data_manager.s1.gradients[:, :] = g0.function.vector().get_local().reshape((self.material.data_manager.n, g0.shape))
         # integrate the behaviour
         mgis_bv.integrate(self.material.data_manager, self.integration_type,
                           0, 0, self.material.data_manager.n);
@@ -129,13 +142,26 @@ class MFrontNonlinearProblem(NonlinearProblem):
                                                        mgis_bv.FiniteStrainTangentOperator.DPK1_DF)
             self.Ct.vector().set_local(Ctv)
         else:
-            self.stress.vector().set_local(self.material.data_manager.s1.thermodynamic_forces.flatten())
-            self.Ct.vector().set_local(self.material.data_manager.K.flatten())
-
-            sizes = self.material.get_state_variable_sizes()
-            for (s, i) in self.state_variables:
-                size = sizes[i]
-                s.vector().set_local(self.material.data_manager.s1.internal_state_variables[:, i:(i+size)].flatten())
+            self.flux.update(self.material.data_manager.s1.thermodynamic_forces.flatten())
+            K = self.material.data_manager.K
+            buff = 0
+            shapes = [ufl.shape(t) for t in self.flux.tangents]
+            ntot = sum([s[0]*s[1] if len(s)==2 else s[0] for s in shapes])
+            for t in self.flux.tangents:
+                s1, s2 = ufl.shape(t)
+                print(s1, s2)
+                t.vector().set_local(np.lib.stride_tricks.as_strided(K[buff:],
+                                                          shape = (K.size // ntot, s1, s2),
+                                                          strides = (K.strides[1] * ntot,
+                                                                     K.strides[1] * s1,
+                                                                     K.strides[1])).flatten())
+                t.vector().set_local(K.flatten())
+                buff += s1*s2
+            print(self.flux.tangents[0].vector().get_local())
+            # sizes = self.material.get_state_variable_sizes()
+            # for (s, i) in self.state_variables:
+            #     size = sizes[i]
+            #     s.vector().set_local(self.material.data_manager.s1.internal_state_variables[:, i:(i+size)].flatten())
 
     def get_state_variable(self, name=None, position=None):
         if name is not None: