Nonlinear systems of PDE's

Hello there, Could someone please help me with my current fenics problem? Actually I’m working on the PDE problem:

v \in (0, 1), \quad u_z: u(x,y), \quad \mathbb{R}^n \supset \tilde{\omega} \supset \omega\quad \int_{\omega} (v^2 + \eta_{\epsilon})|\nabla u_z|^2 \, dx + G_c \int_{\tilde{\omega}} \left(\frac{1-v^2}{4\epsilon} + \epsilon |\nabla v|^2\right) \, dx

This is the strong form of the problem and the first integral is with respect an elastic energy term and the other one is about fracture energy. Applying the Gateaux derivative and simplifying, I got the following system of equations:

\int_{\omega} (v^2 + \eta_{\epsilon}) \nabla u \cdot \nabla \psi_1 \, dx = 0

\int_{\omega} |\nabla u|^2 \cdot v \cdot \psi_2 \, dx + G_c \int_{\omega} \left(v \cdot \psi_2 - \psi_2\right) \cdot \frac{1}{2\epsilon} + 2\epsilon \cdot \nabla \psi_2 \cdot \nabla v \, dx = 0

And I have the restriction that the fracture term cannot decrease, so I have to start with the function v, by definition of the problem as 1 and the next solution, needs to be lesser than 1 (because there is 1-v² in the fracture energy term). The code that I’ve done, is the following:

from __future__ import print_function
from petsc4py import PETSc
from fenics import *
import matplotlib.pyplot as plt
import numpy as np
import time

combined_file = XDMFFile("combined_solution.xdmf")
combined_file.parameters["flush_output"] = True  # Ensure data is written immediately

T = 25.0            # final time
num_steps = 500     # number of time steps
dt = T / num_steps # time step size
tol = 1E-14; eta = 1*10E-7; eps = 4*10E-2; delta = 0.5; Alternate_tol = 1*10E-2; max_iter = 40
G_c = 1

mesh = UnitSquareMesh(5,5)
mesh_h = UnitSquareMesh(5,5)

V = FunctionSpace(mesh, "Lagrange", 1)
g_t = Expression('t', degree=1, t=0)

def boundary_D(x, on_boundary):
    return on_boundary and (near(x[0], 0, tol) or near(x[0], 5, tol))

bc = DirichletBC(V, g_t, boundary_D)

u = Function(V)
v = Function(V)
u1 = TestFunction(V)
v1 = TestFunction(V)
v_previous = Function(V)
v_previous.interpolate(Constant(1.0))  # Initialize v0 as 1.0

F1 = (v_previous**2 + eta)*u*u1*dx
F2 = ((grad(u)**2)*v*v1 + 2*G_c*(v-1)*v1/(4*eps) - eps*inner(grad(v), grad(v1)))*dx

vtkfile = File('Brittle1.pvd')

for t in range(num_steps):
    # Update the projection of u at the beginning of each time step

    check_tol = False
    i = 0
    while i <= max_iter:
        maxnorm = 1
        solve(F1 == 0, u, bc)
        # Solve for v using the previous time step's solution as an upper bound
        solve(F2 == 0, v)
        u_vec = as_backend_type(u.vector()).vec()
        v_vec = as_backend_type(v.vector()).vec()
        v_previous_vec = as_backend_type(v_previous.vector()).vec()
        diff_vec = v_vec.copy()  
        diff_vec.axpy(-1.0, v_previous_vec)  
        maxnorm = diff_vec.norm(PETSc.NormType.INFINITY)  # Calculate the maximum norm

        print('\n',f"Iteration {i}: MaxNorm(v) = {maxnorm}",'\n')

        v_local = v.vector().get_local()
        v_previous_local = v_previous.vector().get_local()
        condition_1 = np.any(v_local < v_previous_local)
        condition_2 = np.all((v_local >= 0) & (v_local <= 1))
        if i > 0:
            if (condition_1 and condition_2):
                v_previous.vector().set_local(np.minimum(v.vector().get_local(), v_previous.vector().get_local()))  # Adjust values within the v_previous Function
                v_previous.assign(v)  # Update v_previous with the values of v
                i += 1

            if(maxnorm < Alternate_tol):
                check_tol = True
        i = i + 1
    # Update for the next time step
    t += dt 
    g_t.t = t
    combined_file.write(u, t)  # Save solution of u
    combined_file.write(v, t)  # Save solution of v
    vtkfile << (v, t)



My idea was to use interpolation for setting v as 1, and then, mix Picard iteration method with newton’s one (that is already being done in fenics), and at each time step, iterating until a it gets to a convergence criterion (In my case, the max norm between the v and v _previous) and update v only if it’s not bigger than the actual solution, because fractures are not allowed to help.

But, the problem that I’ve been stuck with for straight days, is that the method doesn’t converge, the error only increases, and I never get the right fracture evolution, and I don’t know why. Could someone help me please with this problem?

Please consider reading Phase field fracture implementation in FEniCS.

Hello there, sorry for taking a while to response you. I tried to implement the code from this papper:

from dolfin import *

mesh = Mesh('mesh.xml')
V = FunctionSpace(mesh, 'CG', 1)
W = VectorFunctionSpace(mesh, 'CG', 1)
WW = FunctionSpace(mesh, 'DG', 0)
p, q = TrialFunction(V), TestFunction(V)
u, v = TrialFunction(W), TestFunction(W)
Gc = 2.7; l = 0.015; lmbda = 121.1538e3; mu = 80.7692e3

def epsilon(u):
    return sym(grad(u))

def sigma(u):
    return 2.0*mu*epsilon(u)+lmbda*tr(epsilon(u))*Identity(len(u))

def psi(u):
    return 0.5*(lmbda+mu)*(0.5*(tr(epsilon(u))+abs(tr(epsilon(u)))))**2+\

def H(uold,unew,Hold):
    return conditional(lt(psi(uold),psi(unew)),psi(unew),Hold)

top = CompiledSubDomain("near(x[1], 0.5) && on_boundary")
bot = CompiledSubDomain("near(x[1], -0.5) && on_boundary")
def Crack(x):
    return abs(x[1]) < 1e-03 and x[0] <= 0.0

load = Expression("t", t = 0.0, degree=1)
bcbot= DirichletBC(W, Constant((0.0,0.0)), bot)
bctop = DirichletBC(W.sub(1), load, top)
bc_u = [bcbot, bctop]
bc_phi = [DirichletBC(V, Constant(1.0), Crack)]
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
ds = Measure("ds")(subdomain_data="domainBoundaries")

n = FacetNormal(mesh)
unew, uold = Function(W), Function(W)
pnew, pold, Hold = Function(V), Function(V), Function(V)
pnew, pold, Hold = Function(V), Function(V), Function(V)
E_du = ((1.0-pold)**2)*inner(grad(v),sigma(u))*dx
E_phi = (Gc*l*inner(grad(p),grad(q))+((Gc/l)+2.0*H(uold,unew,Hold))\
p_disp = LinearVariationalProblem(lhs(E_du), rhs(E_du), unew, bc_u)
p_phi = LinearVariationalProblem(lhs(E_phi), rhs(E_phi), pnew, bc_phi)
solver_disp = LinearVariationalSolver(p_disp)
solver_phi = LinearVariationalSolver(p_phi)
t = 0
u_r = 0.007
deltaT = 0.1
tol = 1e-3
conc_f = File ("./ResultsDir/phi.pvd")
fname = open('ForcevsDisp.txt', 'w')
while t<=1.0:
    t += deltaT
    if t >=0.7:
        deltaT = 0.0001
    iter = 0
    err = 1
    while err > tol:
        iter +=1
    err_u = errornorm(unew,uold,norm_type = 'l2',mesh = None)
    err_phi = errornorm(pnew,pold,norm_type = 'l2',mesh = None)
    err = max(err_u,err_phi)
    Hold.assign(project(psi(unew), WW))
    if err < tol:
        print ('Iterations:', iter, ', Total time', t)
    if round(t*1e4) % 10 == 0:
        conc_f << pnew

    Traction = dot(sigma(unew),n)
    fy = Traction[1]*ds(1)
    fname.write(str(t*u_r) + "\t")
    fname.write(str(assemble(fy)) + "\n")

print ("Simulation completed")

And i got the following error:

Traceback (most recent call last):
  File "/home/jandui/Documents/MATH code/First FEM approx/", line 2, in <module>
    mesh = Mesh('mesh.xml')

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
*** -------------------------------------------------------------------------
*** Error:   Unable to read data from XML file.
*** Reason:  Unable to open file "mesh.xml".
*** Where:   This error was encountered inside XMLFile.cpp.
*** Process: 0
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  ubuntu
*** -------------------------------------------------------------------------

And, i’m working with a simple version of this one, without considering any tensor equation, only the module of the gradient of the function u. Another question is why it’s not correct my approach, cause i think that there is a similar logic in both of them. Thanks.

You can deal with this by 2 approaches.

  1. Use the below piece of code to generate the mesh using mshr (not maintained any longer, therefore can use UnitSquareMesh(n, n) if mshr doesn’t work).
from dolfin import *
from mshr import *

and refine mesh at the crack tip using the code given here.
Reference : Phase Field Fracture

  1. You can use Gmsh and take advantage of increasing mesh field at a point. After doing meshing, you can export file in .msh format and then refer to Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio on how to import this mesh.

Which version are you talking about ? If I remember the paper has only one code version that you are using.

I’ll look into it when am able to, but one thing I can tell that this form doesn’t match with what you wrote in mathematics above.

1 Like