from dolfin import * from mfront_wrapper.utils import * from mfront_wrapper.gradient_flux 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 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() 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) self.axi_coeff = 2*pi*abs(x[0]) else: self.axi_coeff = 1 self.solver = NewtonSolver() self.state_variables = [] 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, 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 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) # 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) """ return derivative(self.strain_measure(v), v, v) 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): 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[:, :] = 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); # 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.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] flattened_shapes = [s[0]*s[1] if len(s)==2 else s[0] for s in shapes] for (i,t) in enumerate(self.flux.tangents): s = flattened_shapes[i] # print(s, buff, K.shape, K[137,:]) t.vector().set_local(K[:,buff:buff+s].flatten()) #t.vector().set_local(K.flatten()) buff += s # print(self.flux.tangents[0].vector().get_local()) sizes = self.material.get_state_variable_sizes() for i in range(len(self.state_variables)): size = sizes[i] s_values = self.material.data_manager.s1.internal_state_variables[:, i:(i+size)].flatten() state = self.state_variables[i] if isinstance(state, Flux): state.function.vector().set_local(s_values) shapes = [ufl.shape(t) for t in state.tangents] print(shapes) flattened_shapes = [s[0]*s[1] if len(s)==2 else s[0] if len(s)==1 else 1 for s in shapes] # FIXME: Assumes gradients of state variables are defined after fluxes and in right order for (j,t) in enumerate(state.tangents): s = flattened_shapes[j] t.vector().set_local(K[:,buff:buff+s].flatten()) buff += s else: state.vector().set_local(s_values) 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)