Thank you @dokken for the reference, i followed the link and changed my snes_problem.py file and code accordingly but now this error comes:
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
<ipython-input-21-e028e09a1a81> in <module>
2 E_alpha_alpha = ufl.derivative(E_alpha,alpha,ufl.TrialFunction(V_alpha))
3
----> 4 damage_problem = SNESProblem(E_alpha, E_alpha_alpha, alpha, bcs_alpha)
5
6 b = dolfinx.cpp.la.create_vector(V_alpha.dofmap.index_map, V_alpha.dofmap.index_map_bs)
~/Documents/FEM_FEniCS_tuts/newfrac-fenicsx-training/notebooks/python/snes_problem.py in __init__(self, F, J, u, bcs)
28 # Create matrix and vector to be used for assembly
29 # of the non-linear problem
---> 30 self.A = dolfinx.fem.create_matrix(self.a)
31 self.b = dolfinx.fem.create_vector(self.L)
32
/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py in form(form, dtype, form_compiler_params, jit_params)
137 return form
138
--> 139 return _create_form(form)
140
141
/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py in _create_form(form)
132 return form argument"""
133 if isinstance(form, ufl.Form):
--> 134 return _form(form)
135 elif isinstance(form, collections.abc.Iterable):
136 return list(map(lambda sub_form: _create_form(sub_form), form))
/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py in _form(form)
126 _cpp.fem.IntegralType.vertex: subdomains.get("vertex")}
127
--> 128 return formcls(ufcx_form, V, coeffs, constants, subdomains, mesh, code)
129
130 def _create_form(form):
/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py in __init__(self, form, V, coeffs, constants, subdomains, mesh, code)
46 self._ufcx_form = form
47 ffi = cffi.FFI()
---> 48 super().__init__(ffi.cast("uintptr_t", ffi.addressof(self._ufcx_form)),
49 V, coeffs, constants, subdomains, mesh)
50
RuntimeError: Incompatible mesh
Code in snes_problem.py:
from typing import List
from petsc4py import PETSc
import dolfinx
import ufl
class SNESProblem:
"""Nonlinear problem class compatible with PETSC.SNES solver.
"""
def __init__(self, F: ufl.form.Form, J: ufl.form.Form, u: dolfinx.Function, bcs: List[dolfinx.DirichletBC]):
"""This class set up structures for solving a non-linear problem using Newton's method.
Parameters
==========
F: Residual.
J: Jacobian.
u: Solution.
bcs: Dirichlet boundary conditions.
"""
self.L = F
self.a = J
self.bcs = bcs
self.u = u
# Create matrix and vector to be used for assembly
# of the non-linear problem
self.A = dolfinx.fem.create_matrix(self.a)
self.b = dolfinx.fem.create_vector(self.L)
def F(self, snes: PETSc.SNES, x: PETSc.Vec, b: PETSc.Vec):
"""Assemble the residual F into the vector b.
Parameters
==========
snes: the snes object
x: Vector containing the latest solution.
b: Vector to assemble the residual into.
"""
# We need to assign the vector to the function
x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
x.copy(self.u.vector)
self.u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
# Zero the residual vector
with b.localForm() as b_local:
b_local.set(0.0)
dolfinx.fem.assemble_vector(b, self.L)
# Apply boundary conditions
dolfinx.fem.apply_lifting(b, [self.a], [self.bcs], [x], -1.0)
b.ghostUpdate(addv=PETSc.InsertMode.ADD,mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.set_bc(b, self.bcs, x, -1.0)
def J(self, snes, x: PETSc.Vec, A: PETSc.Mat, P: PETSc.Mat):
"""Assemble the Jacobian matrix.
Parameters
==========
x: Vector containing the latest solution.
A: Matrix to assemble the Jacobian into.
"""
A.zeroEntries()
dolfinx.fem.assemble_matrix(A, self.a, self.bcs)
A.assemble()
Code now -
import matplotlib.pyplot as plt
import numpy as np
import dolfinx
import dolfinx.plot
import ufl
from mpi4py import MPI
from petsc4py import PETSc
import pyvista
from pyvista.utilities.xvfb import start_xvfb
start_xvfb(wait=0.5)
import sys
sys.path.append("/home/ratan/Documents/FEM_FEniCS_tuts/newfrac-fenicsx-training/notebooks/python/")
from snes_problem import SNESProblem
from plots import plot_damage_state
L = 1.; H = 0.3;
ell_ = 0.1
cell_size = ell_/6;
nx = int(L/cell_size)
ny = int(H/cell_size)
mesh = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD, [(0.0, 0.0), (L, H)], [nx, ny])
ndim = mesh.geometry.dim
pyvista.OFF_SCREEN
dolfinx.plot.create_vtk_mesh(mesh, ndim)
grid = pyvista.UnstructuredGrid(*dolfinx.plot.create_vtk_mesh(mesh, ndim))
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True, show_scalar_bar=True)
plotter.view_xy()
plotter.add_axes()
plotter.set_scale(5,5)
#plotter.reset_camera(render=True, bounds=(-L/2, L/2, -H/2, H/2, 0, 0))
# if not pyvista.OFF_SCREEN:
# plotter.show()
# from pathlib import Path
# Path("output").mkdir(parents=True, exist_ok=True)
# figure = plotter.screenshot("output/mesh.png")
element_u = ufl.VectorElement('Lagrange',mesh.ufl_cell(),degree=1,dim=2)
V_u = dolfinx.fem.FunctionSpace(mesh, element_u)
element_alpha = ufl.FiniteElement('Lagrange',mesh.ufl_cell(),degree=1)
V_alpha = dolfinx.fem.FunctionSpace(mesh, element_alpha)
# Define the state
u = dolfinx.fem.Function(V_u, name="Displacement")
alpha = dolfinx.fem.Function(V_alpha, name="Damage")
state = {"u": u, "alpha": alpha}
# need upper/lower bound for the damage field
alpha_lb = dolfinx.fem.Function(V_alpha, name="Lower bound")
alpha_ub = dolfinx.fem.Function(V_alpha, name="Upper bound")
# Measures
dx = ufl.Measure("dx",domain=mesh)
ds = ufl.Measure("ds",domain=mesh)
def bottom(x):
return np.isclose(x[1], 0.0)
def right(x):
return np.isclose(x[0], L)
def top(x):
return np.isclose(x[1], H)
def left(x):
return np.isclose(x[0], 0.0)
blocked_dofs_left_u = dolfinx.fem.locate_dofs_geometrical((V_u.sub(0), V_u.sub(0).collapse()[0]), left)
blocked_dofs_right_u = dolfinx.fem.locate_dofs_geometrical((V_u.sub(0), V_u.sub(0).collapse()[0]), right)
blocked_dofs_bottom_u = dolfinx.fem.locate_dofs_geometrical((V_u.sub(1), V_u.sub(1).collapse()[0]), bottom)
blocked_dofs_left_alpha = dolfinx.fem.locate_dofs_geometrical(V_alpha, left)
blocked_dofs_right_alpha = dolfinx.fem.locate_dofs_geometrical(V_alpha, right)
zero_u = dolfinx.fem.Function(V_u.sub(0).collapse()[0])
with zero_u.vector.localForm() as bc_local:
bc_local.set(0.0)
nonzero_u = dolfinx.fem.Function(V_u.sub(0).collapse()[0])
with nonzero_u.vector.localForm() as bc_local:
bc_local.set(1.0)
one_alpha = dolfinx.fem.Function(V_alpha)
with one_alpha.vector.localForm() as bc_local:
bc_local.set(1.0)
zero_alpha = dolfinx.fem.Function(V_alpha)
with zero_alpha.vector.localForm() as bc_local:
bc_local.set(0.0)
bc_u0 = dolfinx.fem.dirichletbc(zero_u, blocked_dofs_left_u, V_u.sub(1))
bc_u1 = dolfinx.fem.dirichletbc(nonzero_u, blocked_dofs_right_u, V_u.sub(1))
bc_u2 = dolfinx.fem.dirichletbc(zero_u, blocked_dofs_bottom_u, V_u.sub(0))
bc_alpha0 = dolfinx.fem.dirichletbc(zero_alpha, blocked_dofs_left_alpha)
bc_alpha1 = dolfinx.fem.dirichletbc(zero_alpha, blocked_dofs_right_alpha)
bcs_u = [bc_u0,bc_u1,bc_u2]
bcs_alpha = [bc_alpha0,bc_alpha1]
# setting the upper bound to 0 where BCs are applied
alpha_ub.interpolate(one_alpha)
dolfinx.fem.petsc.set_bc(alpha_ub.vector, bcs_alpha)
E, nu = dolfinx.fem.Constant(mesh, 100.0), dolfinx.fem.Constant(mesh, 0.3)
Gc = dolfinx.fem.Constant(mesh, 1.0)
ell = dolfinx.fem.Constant(mesh, ell_)
def w(alpha):
"""Dissipated energy function as a function of the damage """
return alpha
def a(alpha, k_ell=1.e-6):
"""Stiffness modulation as a function of the damage """
return (1 - alpha) ** 2 + k_ell
def eps(u):
"""Strain tensor as a function of the displacement"""
return ufl.sym(ufl.grad(u))
def sigma_0(u):
"""Stress tensor of the undamaged material as a function of the displacement"""
mu = E / (2.0 * (1.0 + nu))
lmbda = E * nu / (1.0 - nu ** 2)
return 2.0 * mu * eps(u) + lmbda * ufl.tr(eps(u)) * ufl.Identity(ndim)
def sigma(u,alpha):
"""Stress tensor of the damaged material as a function of the displacement and the damage"""
return a(alpha) * sigma_0(u)
import sympy
z = sympy.Symbol("z")
c_w = 4*sympy.integrate(sympy.sqrt(w(z)),(z,0,1))
print("c_w = ",c_w)
c_1w = sympy.integrate(sympy.sqrt(1/w(z)),(z,0,1))
print("c_1/w = ",c_1w)
tmp = 2*(sympy.diff(w(z),z)/sympy.diff(1/a(z),z)).subs({"z":0})
sigma_c = sympy.sqrt(tmp * Gc.value * E.value / (c_w * ell.value))
print("sigma_c = %2.3f"%sigma_c)
eps_c = float(sigma_c/E.value)
print("eps_c = %2.3f"%eps_c)
f = dolfinx.fem.Constant(mesh,(0.,0.))
elastic_energy = 0.5 * ufl.inner(sigma(u,alpha), eps(u)) * dx
dissipated_energy = Gc / float(c_w) * (w(alpha) / ell + ell * ufl.dot(ufl.grad(alpha), ufl.grad(alpha))) * dx
external_work = ufl.dot(f, u) * dx
total_energy = elastic_energy + dissipated_energy - external_work
E_u = ufl.derivative(total_energy,u,ufl.TestFunction(V_u))
E_du = ufl.replace(E_u,{u: ufl.TrialFunction(V_u)})
problem_u = dolfinx.fem.petsc.LinearProblem(a=ufl.lhs(E_du), L=ufl.rhs(E_du), bcs=bcs_u, u=u,
petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
problem_u.solve()
load = .6
with nonzero_u.vector.localForm() as bc_local:
bc_local.set(load)
problem_u.solve()
plot_damage_state(state,load=load)
EE_alpha = ufl.derivative(total_energy,alpha,ufl.TestFunction(V_alpha))
E_alpha_alpha = ufl.derivative(E_alpha,alpha,ufl.TrialFunction(V_alpha))
damage_problem = SNESProblem(E_alpha, E_alpha_alpha, alpha, bcs_alpha)
b = dolfinx.cpp.la.create_vector(V_alpha.dofmap.index_map, V_alpha.dofmap.index_map_bs)
J = dolfinx.fem.create_matrix(damage_problem.a)
# Create Newton solver and solve
solver_alpha_snes = PETSc.SNES().create()
solver_alpha_snes.setType("vinewtonrsls")
solver_alpha_snes.setFunction(damage_problem.F, b)
solver_alpha_snes.setJacobian(damage_problem.J, J)
solver_alpha_snes.setTolerances(rtol=1.0e-9, max_it=50)
solver_alpha_snes.getKSP().setType("preonly")
solver_alpha_snes.getKSP().setTolerances(rtol=1.0e-9)
solver_alpha_snes.getKSP().getPC().setType("lu")
# We set the bound (Note: they are passed as reference and not as values)
solver_alpha_snes.setVariableBounds(alpha_lb.vector,alpha_ub.vector)
Thanks for your patience and time.