from dolfin import * from mfront_wrapper.utils import * import mgis.behaviour as mgis_bv 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.integration_type = mgis_bv.IntegrationType.IntegrationWithConsistentTangentOperator 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, dim=v.geometric_dimension()) 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") mgis_bv.integrate(self.material.data_manager, self.integration_type, 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 mgis_bv.integrate(self.material.data_manager, self.integration_type, 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()) 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: 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.") shape = self.material.get_state_variable_sizes()[position] 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 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)