Modifs in utils

mfront_wrapper/mfront_wrapper.py

-import mgis.behaviour as mgis_bv
-
-mgis_hypothesis = {"plane_strain": mgis_bv.Hypothesis.PlaneStrain,
-              "plane_stress": mgis_bv.Hypothesis.PlaneStress,
-              "3d": mgis_bv.Hypothesis.Tridimensional,
-              "axisymmetric": mgis_bv.Hypothesis.Axisymmetrical}
-
-
-class NonlinearMaterial:
-    def __init__(self, path, name, hypothesis="3d",
-                 material_properties={}, external_state_variables={"Temperature": 293.15}):
-        self.path = path
-        self.name = name
-        # Defining the modelling hypothesis
-        self.hypothesis = mgis_hypothesis[hypothesis]
-        self.material_properties = material_properties
-        self.external_state_variables = external_state_variables
-        # Loading the behaviour
-        self.behaviour = mgis_bv.load(self.path, self.name, self.hypothesis)
-
-    def set_data_manager(self, ngauss):
-        # Setting the material data manager
-        self.data_manager = mgis_bv.MaterialDataManager(self.behaviour, ngauss)
-        for s in [self.data_manager.s0, self.data_manager.s1]:
-            for (key, value) in self.material_properties.items():
-                mgis_bv.setMaterialProperty(s, key, value)
-            for (key, value) in self.external_state_variables.items():
-                mgis_bv.setExternalStateVariable(s, key, value)
-
-    def get_state_variable_names(self):
-        print(self.behaviour.internal_state_variables[0].getVariableOffsetByString())
-        return [svar.name for svar in self.behaviour.internal_state_variables]
-
-class MFrontNonlinearProblem(NonlinearProblem):
-    def __init__(self, u, material, quadrature_degree=2):
-        NonlinearProblem.__init__(self)
-        self.u = u
-        self.V = self.u.function_space()
-        self.u_ = TestFunction(self.V)
-        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.quadrature_degree = quadrature_degree
-        self.set_quadrature_function_spaces()
-
-        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])
-        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.solver = NewtonSolver()
-
-        self.state_variables = []
-
-    def set_loading(self, Fext):
-        # adds external forces to residual form
-        self.L -= 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 = 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
-        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 strain_measure(self, v):
-        """ Strain measure associated with stress measure:
-            * small strain behaviour: linearized strain tensor epsilon = sym(grad(u))
-            * finite strain behaviour: transformation gradient F = Id + grad(u)
-        """
-        if self.axisymmetric:
-            r = abs(SpatialCoordinate(self.mesh)[0])
-            g = axi_grad(r, v)
-            E = symmetric_tensor_to_vector(sym(g))
-            if v.geometric_dimension()==2:
-                return as_vector([E[i] for i in range(4)])
-        else:
-            g = grad(v)
-            if self.finite_strain:
-                return transformation_gradient(g)
-            else:
-                return symmetric_gradient(g)
-
-    def strain_variation(self, v):
-        """ Variation of strain measure associated with stress measure:
-            * small strain behaviour: linearized strain tensor d_epsilon = sym(grad(du))
-            * finite strain behaviour: displacement gradient d_F = grad(du)
-        """
-        if self.axisymmetric:
-            r = abs(SpatialCoordinate(self.mesh)[0])
-            g = axi_grad(r, v)
-            E = symmetric_tensor_to_vector(sym(g))
-            if v.geometric_dimension()==2:
-                return as_vector([E[i] for i in range(4)])
-        else:
-            g = grad(v)
-            if self.finite_strain:
-                return gradient(g)
-            else:
-                return symmetric_gradient(g)
-
-    def initialize_fields(self):
-        self.stress = Function(self.Wsig, name="Current stress")
-        self.strain = Function(self.Wsig, name="Current strain increment")
-        self.Ct = Function(self.WCt, name="Consistent tangent operator")
-
-        it = mgis_bv.IntegrationType.PredictionWithElasticOperator
-        mgis_bv.integrate(self.material.data_manager, it, 0, 0, self.material.data_manager.n);
-        if self.finite_strain:
-            local_project(self.strain_measure(self.u), self.Wsig, self.dx, self.strain)
-            # copy the strain values to `MGIS`
-            self.material.data_manager.s0.gradients[:, :] = self.strain.vector().get_local().reshape((self.material.data_manager.n, self.strain_dim))
-        else:
-            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)
-        # 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))
-        # integrate the behaviour
-        it = mgis_bv.IntegrationType.IntegrationWithConsistentTangentOperator
-        mgis_bv.integrate(self.material.data_manager, it, 0, 0, self.material.data_manager.n);
-        # getting the stress and consistent tangent operator back to
-        # the FEniCS world.
-        if self.finite_strain:
-            pk1v = self.stress.vector().get_local()
-            mgis_bv.convertFiniteStrainStress(pk1v, self.material.data_manager,
-                                              mgis_bv.FiniteStrainStress.PK1)
-            self.stress.vector().set_local(pk1v)
-            Ctv = self.Ct.vector().get_local()
-            mgis_bv.convertFiniteStrainTangentOperator(Ctv, self.material.data_manager,
-                                                       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())
-            for (s, i) in self.state_variables:
-                s.vector().set_local(self.material.data_manager.s1.internal_state_variables[:, i])
-
-    def register_state_variable(self, name=None, position=None, shape=()):
-        if name is not None:
-            position = self.material.get_state_variable_names().index(name)
-        elif position is not None:
-           name = self.material.get_state_variable_names()[position]
-        else:
-            raise ValueError("Name or position of state variable must be specified.")
-        We = get_quadrature_element(self.mesh.ufl_cell(), self.quadrature_degree, shape)
-        W = FunctionSpace(self.mesh, We)
-        self.state_variables.append([Function(W, name=name), position])
-        return self.state_variables[-1][0]
-
-    def get_state_variable(self, var, name=None, position=None):
-        if name is not None:
-            position = self.material.get_state_variable_names().index(name)
-        var.vector().set_local(self.material.data_manager.s1.internal_state_variables[:, position].flatten())
-
-    def form(self, A, P, b, x):
-        self.update_constitutive_law(self.u)
-        assemble_system(self.a, self.L, A_tensor=A, b_tensor=b, bcs=self.bc, x0=x)
-
-    def F(self,b,x):
-        pass
-
-    def J(self,A,x):
-        pass
-
-    def solve(self, x):
-        self.solver.solve(self, x)
-        mgis_bv.update(self.material.data_manager)
-

mfront_wrapper/nonlinear_problem.py

@@ -18,6 +18,8 @@ class MFrontNonlinearProblem(NonlinearProblem):
         self.quadrature_degree = quadrature_degree
         self.set_quadrature_function_spaces()
 
+
+
         self.bc = []
         self.dx = Measure("dx", metadata={"quadrature_degree": self.quadrature_degree,
                                           "quadrature_scheme": "default"})
@@ -38,6 +40,9 @@ class MFrontNonlinearProblem(NonlinearProblem):
 
         self.state_variables = []
 
+    def add_stress(self, gradients=[]):
+        pass
+
     def set_loading(self, Fext):
         # adds external forces to residual form
         self.L -= ufl.replace(Fext, {self.u: self.u_})
@@ -89,18 +94,7 @@ class MFrontNonlinearProblem(NonlinearProblem):
             * small strain behaviour: linearized strain tensor d_epsilon = sym(grad(du))
             * finite strain behaviour: displacement gradient d_F = grad(du)
         """
-        if self.axisymmetric:
-            r = abs(SpatialCoordinate(self.mesh)[0])
-            g = axi_grad(r, v)
-            E = symmetric_tensor_to_vector(sym(g))
-            if v.geometric_dimension()==2:
-                return as_vector([E[i] for i in range(4)])
-        else:
-            g = grad(v)
-            if self.finite_strain:
-                return gradient(g)
-            else:
-                return symmetric_gradient(g)
+        return derivative(self.strain_measure(v), v, v)
 
     def initialize_fields(self):
         self.stress = Function(self.Wsig, name="Current stress")

mfront_wrapper/utils.py

@@ -36,6 +36,8 @@ def symmetric_tensor_to_vector(T, T22=0):
         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], sqrt(2)*T[0, 1], sqrt(2)*T[0, 2], sqrt(2)*T[1, 2]])
+    elif len(ufl.shape(T)) == 1:
+        return T
     else:
         raise NotImplementedError
 
@@ -46,6 +48,8 @@ def nonsymmetric_tensor_to_vector(T, T22=0):
          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]])
+    elif len(ufl.shape(T)) == 1:
+        return T
     else:
         raise NotImplementedError
 
@@ -85,9 +89,13 @@ def get_quadrature_element(cell, degree, dim=0):
     if dim in [0, 1, (), (0,), (1,)]:
         return FiniteElement("Quadrature", cell, degree=degree,
                             quad_scheme='default')
-    elif type(dim) == int or len(dim)==1:
+    elif type(dim) == int:
         return VectorElement("Quadrature", cell, degree=degree,
                               dim=dim, 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')
     elif type(dim) == tuple:
         return TensorElement("Quadrature", cell, degree=degree, shape=dim, quad_scheme='default')
     else: