From f0bc90dee63b1c6ed9e7eebb788e45f5f7ca87d3 Mon Sep 17 00:00:00 2001 From: "jeremy.bleyer" Date: Mon, 24 Feb 2020 15:15:44 +0100 Subject: [PATCH] Modifs in utils --- mfront_wrapper/mfront_wrapper.py | 201 ---------------------------- mfront_wrapper/nonlinear_problem.py | 18 +-- mfront_wrapper/utils.py | 10 +- 3 files changed, 15 insertions(+), 214 deletions(-) delete mode 100644 mfront_wrapper/mfront_wrapper.py diff --git a/mfront_wrapper/mfront_wrapper.py b/mfront_wrapper/mfront_wrapper.py deleted file mode 100644 index 281bb60..0000000 --- a/mfront_wrapper/mfront_wrapper.py +++ /dev/null @@ -1,201 +0,0 @@ -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) - diff --git a/mfront_wrapper/nonlinear_problem.py b/mfront_wrapper/nonlinear_problem.py index bdfc6e0..f3da3d5 100644 --- a/mfront_wrapper/nonlinear_problem.py +++ b/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") diff --git a/mfront_wrapper/utils.py b/mfront_wrapper/utils.py index cbde703..bff808a 100644 --- a/mfront_wrapper/utils.py +++ b/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: -- GitLab