How to condition my modified 3D elasticity FEniCS code?

Linear Algebra
1

Hello,

I am trying to solve a 3D elasticity problem with an additional field for eigenstrains. Geometry is hemisphere with tetrahedral meshing. I have checked the code for syntax: no errors, tuning/relaxing the tolerance parameters and also varied other solver parameters.

When I monitor convergence, the norm either does not decrease or in very rare cases decreases very very slowly. Although the weak form is mathematically correct, I believe there is some problem with the conditioning of stiffness (or global tangent) matrix of the mixed system. I am not sure how to approach this issue.

A MWE is given here: (most of the code is packed in functions)




from copy import deepcopy
from statistics import mode

from numpy import infty, iinfo, int32, argwhere, min, max

from dolfin import *


def sphere_geometry(db_mesh_folder, mesh_file):
     
    #R, mesh_size = geo_parameters

    mymesh = Mesh()
    
    #db_mesh_folder = "/Users/dhariwal/Library/CloudStorage/OneDrive-VirginiaTech/Zirconia/Particle_SIMULATION/Code/single_variant/meshes/"
    #mesh_file = "mesh3D_R_" + str(int(self.R)) + "_size_" + str(int(mesh_size)) + ".xdmf"

    with XDMFFile(db_mesh_folder + mesh_file) as infile:
        infile.read(mymesh)

    print("Domain created, meshing done")

    print("Radius = ", R)
    print("mesh_size = ", mesh_size)
    print("Number of cells = ", mymesh.num_cells())
    print("Number of vertices = ", mymesh.num_vertices())

    # Characteristic size of the elements of the mesh 
    hMESH = CellDiameter(mymesh)
    meshSIZE = 0.5*(mymesh.hmax() + mymesh.hmin())

    print("meshSIZE = ", meshSIZE)

    
    return mymesh, hMESH, meshSIZE


    x1, y1, z1 = geo_parameters[0]
    nx, ny, nz = geo_parameters[1]

    corner1 = Point(0.0, 0.0, 0.0)
    corner2 = Point(x1, y1, z1)

    mymesh = BoxMesh(corner1, corner2, nx, ny, nz)

    print("Number of cells = ", mymesh.num_cells())
    print("Number of vertices = ", mymesh.num_vertices())

    return mymesh

def cuboid_geometry(geo_parameters):

    x1, y1, z1 = geo_parameters[0]
    nx, ny, nz = geo_parameters[1]

    corner1 = Point(0.0, 0.0, 0.0)
    corner2 = Point(x1, y1, z1)

    mymesh = BoxMesh(corner1, corner2, nx, ny, nz)

    print("Number of cells = ", mymesh.num_cells())
    print("Number of vertices = ", mymesh.num_vertices())

    # Characteristic size of the elements of the mesh 
    hMESH = CellDiameter(mymesh)
    meshSIZE = 0.5*(mymesh.hmax() + mymesh.hmin())

    return mymesh, hMESH, meshSIZE


def func_spaces(mesh, basis):

    if basis == 'linear':
        shape_func = 1
    else:
        shape_func = 2

    # basis for 'vector' displacement
    P1 = VectorElement("Lagrange", mesh.ufl_cell(), shape_func)
    
    # basis for 'scalar' monovariant phase field
    P2 = FiniteElement("Lagrange", mesh.ufl_cell(), shape_func)
    
    W1 = FunctionSpace(mesh, P1)
    W2 = FunctionSpace(mesh, P2)
    
    
    
    from db_variant import Variant
    shape_strains_list = Variant.CuAlNi()

    shape_strains_list = as_tensor(shape_strains_list[0])


    num_variants = 1 

    # single variant
    element = MixedElement([P1, P2])


    # create the space of solutions over the domain 
    W = FunctionSpace(mesh, element)

    # a generic function from W
    w  = Function(W)
    dw = TrialFunction(W)
    wx = TestFunction(W)


    print("Total D.O.F. = ", W.dim())


    return W, W1, W2, w, dw, wx, num_variants, shape_strains_list



    if basis == 'linear':
        shape_func = 1
    else:
        shape_func = 2

    # basis for 'vector' displacement
    P1 = VectorElement("Lagrange", mesh.ufl_cell(), shape_func)
    W = FunctionSpace(mesh, P1)

    w  = Function(W)               # a generic function from W
    dw = TrialFunction(W)
    v  = TestFunction(W)             # test function for variational problem

    num_variants = 0         

    print("This is a pure elastic problem.\n No. of variants = ", num_variants)
    print("Total D.O.F. = ", W.dim())

    return W, w, dw, v, P1





    # Spatial coordinates of mesh
    x = SpatialCoordinate(mesh)

    # Definition of tolerance
    tol = DOLFIN_EPS # FEniCS tolerance - necessary for "boundary_markers" and the correct definition of the boundaries 

    #class Bottom(SubDomain):
    #    def inside(self, x, on_boundary):
    #        return near(x[2], 0.0)

    #bottom = Bottom()

    bT = CompiledSubDomain('near(x[2], R, tol) && on_boundary', R = R, tol = tol)

    bP = CompiledSubDomain('sqrt(x[0]*x[0] + x[1]*x[1]) < meshSIZE \
                            && near(x[2], R, tol) \
                            && on_boundary', meshSIZE = meshSIZE, R=R, tol = tol)

    # Definition of contact zones
    bC = CompiledSubDomain('on_boundary && x[2] < compression + tol', R = R, tol=tol, compression = compression) # Contact search - contact part of the boundary

    boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
    boundaries.set_all(0)
    
    print("Markers for all boundary facets are set")

    bT.mark(boundaries, 2)      # bottom
    bP.mark(boundaries, 4)      # lateral_penalty
    bC.mark(boundaries, 6)      # contact_zone

    print(boundaries.array()[boundaries.array() > 0])



    cells = MeshFunction("size_t", mesh, mesh.topology().dim())
    cells.set_all(1)

    print("Marker for all internal cells is set to 1")

    



    bc1 = DirichletBC(W.sub(2), Constant((-compression)), boundaries, 2)
    bc2 = DirichletBC(W.sub(0), Constant((0)), boundaries, 4)
    bc3 = DirichletBC(W.sub(1), Constant((0)), boundaries, 4)

    bcs = [bc1,bc2,bc3]


    dx = Measure('dx', subdomain_data=cells)
    ds = Measure('ds', subdomain_data=boundaries)

    return dx, ds, boundaries, bcs





    lower_bound = Function(W)
    upper_bound = Function(W)

    #W00 = W.sub(0).collapse()
    #W11 = W.sub(1).collapse()

    #ninfty, pinfty = Function(W00), Function(W00)
    #_lower, c_upper = Function(W11), Function(W11)


    if u_min is None:
        lower_bound.sub(0).vector()[:] = - infty
        lower_bound.sub(1).vector()[:] = - infty
        lower_bound.sub(2).vector()[:] = - infty
        upper_bound.sub(0).vector()[:] = infty
        upper_bound.sub(1).vector()[:] = infty
        upper_bound.sub(2).vector()[:] = infty

        #lower_bound.vector()[:] = -infty
        #upper_bound.vector()[:] = infty
    else:
        lower_bound.vector()[:] = u_min
        upper_bound.vector()[:] = u_max


    #fa = FunctionAssigner(W, [W00, W11])

    #fa.assign(lower_bound, [ninfty, c_lower])
    #fa.assign(upper_bound, [pinfty, c_upper])

    
    return lower_bound, upper_bound






    
    # Form compiler options
    parameters["form_compiler"]["optimize"]     = True
    parameters["form_compiler"]["cpp_optimize"] = True

    from db_compliance import Compliance

    compliance_phase1, compliance_phase2 = Compliance.CuAlNi()

    Compliance_tetragonal = as_tensor(compliance_phase1)
    #Compliance_monoclinic = as_tensor(compliance_phase2)


    # Definition of normal vectors to the elastic body and to the rigid half-space
    nRIGID = Constant((0.0, 0.0, 1.0))
    n = FacetNormal(mesh)




    # matrix-to-vector transform
    def strain3voigt(e):
        return as_vector(
            [e[0, 0], 
            e[1, 1], 
            e[2, 2], 
            e[1, 2] * 2, 
            e[0, 2] * 2, 
            e[0, 1] * 2]
                        )

    # vector-to-matrix transform
    def voigt3stress(s):
        return as_tensor(
            [[s[0], s[5], s[4]], 
            [s[5], s[1], s[3]], 
            [s[4], s[3], s[2]]]
                        )

    def epsilon(u, strain=None):

        if strain is None:
            e = 0.5 * (grad(u) + grad(u).T)
            return e
        
        elif strain=='lagrangian':
            I = Identity(3)
            Fgrad = I + grad(u)
            C = Fgrad.T * Fgrad
            E = 0.5*(C - I)
            return E
        
        elif strain=='hencky_order_22':
            I = Identity(3)
            Fgrad = I + grad(u)
            C = Fgrad.T * Fgrad
            a = C*C - I
            b = (C*C) + 4*C + I
            H = 1.5 * a * inv(b)
            return H
        
    def Vstrain(w):
        return strain3voigt(epsilon(w))

    def Vsigma(w):
        return dot(Compliance_tetragonal, Vstrain(w))

    def sigma(w):
        return voigt3stress(Vsigma(w))

    def gap(w): # Definition of gap function
        x = SpatialCoordinate(mesh)
        return x[2]+w[2]

    def maculay(x): # Definition of Maculay bracket
        return (x+abs(x))/2
    
    def pN(w):
        return dot(dot(nRIGID,sigma(w)),nRIGID) 

    def convert_to_integer(a):
        "Convert to a 32-bit int. Raise exception if this will cause an overflow"
        if abs(a) <= iinfo(int32).max:
            return int32(a)
        else:
            raise OverflowError("Cannot safely convert to a 32-bit int.")


    def elastic(w):
        elastic_energy = 0.5 * inner(sigma(w), epsilon(w))

        return elastic_energy
    
    

   

    # The displacement u must be such that the current configuration x+u
    # remains in the box [xmin = -inf,xmax = inf] x [ymin = 0,ymax = inf]
    # constraint_u = Expression(("xmax - x[0]","ymax - x[1]","zmax - x[2]"),
    #                         xmax=infty,  ymax=infty, zmax=infty, degree=1)
    # constraint_l = Expression(("xmin - x[0]","ymin - x[1]","zmin - x[2]"),
    #                         xmin=-infty, ymin=-infty, zmin=0.0, degree=1)
    # umin = interpolate(constraint_l, W)
    # umax = interpolate(constraint_u, W)

    #w_converged = Function(W)

    time_end = 1.0
    t_curr = 0.0

    t_curr = t_curr + time_step
    d_step = compression * t_curr

   

    # Define the solver parameters

    while t_curr <= time_end:
        
        

        dx, ds, boundaries, bcs = pureElastic_sphere_boundaries(mesh, meshSIZE, R, d_step, W)
        print("Current displacement at boundary = ", d_step, " nm")#, " and i_count = ", i_count)

        #w = w_converged

        # Stored strain energy density (linear elasticity model)
        psi = elastic(w)

        # Total potential energy
        Pi = psi*dx +  (0.5 * penalty * (E/hMESH) * inner(maculay(-gap(w)), maculay(-gap(w)))) * ds(6)

        # Compute first variation of Pi (directional derivative about u in the direction of v)
        F = derivative(Pi, w, v)
        
        # Compute Jacobian of F for iterative "newton_solver"
        J = derivative(F, w, dw)


        problem = NonlinearVariationalProblem(F, w, bcs, J=J)

        umin, umax = PureElastic_func_space_bounds(W)
        problem.set_bounds(umin, umax)

        solver  = NonlinearVariationalSolver(problem)

        snes_solver_parameters = {"nonlinear_solver": "snes",
                            "snes_solver": {"linear_solver": "lu",# "preconditioner": "ilu",
                                            "maximum_iterations": 100,
                                            "absolute_tolerance": 1E-9,
                                            "relative_tolerance": 1E-9,
                                            "report": False,
                                            "method":'vinewtonrsls',
                                            "error_on_nonconvergence": True}}
                                            # ,"krylov_solver": {"divergence_limit": 1e10}
    
        solver.parameters.update(snes_solver_parameters)

        # Solve the nonlinear problem
        (iter, converged) = solver.solve()

        t_curr = t_curr + time_step
        d_step = compression * t_curr

        if not converged:
            t_curr = t_curr - time_step
            time_step = time_step/2.0
            continue


    # WT = TensorFunctionSpace(mesh, "Lagrange", degree=0)
    # # strain field
    # EPS = Function(WT, name="Strain")
    # EPS.assign(project(epsilon(w), WT))

    # # stress field
    # SIGMA = Function(WT, name="Strain")
    # SIGMA.assign(project(sigma(w), WT))


    return w, boundaries


def func_space_bounds(W, u_min=None, u_max=None, c_min=None, c_max=None):

    lower_bound = Function(W)
    upper_bound = Function(W)


    S = [W.sub(i).collapse() for i in range(num_variants + 1)]


    ninfty, pinfty = Function(S[0]), Function(S[0])

    c_lower = [Function(S[i+1]) for i in range(num_variants)]
    c_upper = [Function(S[i+1]) for i in range(num_variants)]


    if u_min is None:

        ninfty.vector()[:] = - 50
        pinfty.vector()[:] = 50

        if c_min is None:

            for i in range(num_variants):
                c_lower[i].vector()[:] = - 0.1
                c_upper[i].vector()[:] = 1

    else: 
        ninfty.vector()[:] = u_min
        pinfty.vector()[:] = u_max

        for i in range(num_variants):
                c_lower[i].vector()[:] = c_min
                c_upper[i].vector()[:] = c_max



    fa = FunctionAssigner(W, S)

    lb = [ninfty]
    ub = [pinfty]

    for i in c_lower: lb.append(i)
    for i in c_upper: ub.append(i)

    fa.assign(lower_bound, lb)
    fa.assign(upper_bound, ub)

    print(lower_bound.vector()[:], upper_bound.vector()[:])

    return lower_bound, upper_bound


def sphere_boundaries(mesh, meshSIZE, R, compression, W):

    # Spatial coordinates of mesh
    x = SpatialCoordinate(mesh)

    # Definition of tolerance
    tol = DOLFIN_EPS # FEniCS tolerance - necessary for "boundary_markers" and the correct definition of the boundaries 


    bT = CompiledSubDomain('near(x[2], R, tol) && on_boundary', R = R, tol = tol)

    bP = CompiledSubDomain('sqrt(x[0]*x[0] + x[1]*x[1]) < meshSIZE \
                            && near(x[2], R, tol) \
                            && on_boundary', meshSIZE = meshSIZE, R=R, tol = tol)

    # Definition of contact zones
    bC = CompiledSubDomain('on_boundary && x[2] < compression + tol', R = R, tol=tol, compression = compression) # Contact search - contact part of the boundary

    boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
    boundaries.set_all(0)
    

    bT.mark(boundaries, 2)      # bottom
    bP.mark(boundaries, 4)      # lateral_penalty
    bC.mark(boundaries, 6)      # contact_zone

    #print(boundaries.array()[boundaries.array() == 2])
    print(boundaries.array()[boundaries.array() == 4])
    print(boundaries.array()[boundaries.array() == 6])



    cells = MeshFunction("size_t", mesh, mesh.topology().dim())
    cells.set_all(1)



    bc1 = DirichletBC(W.sub(0).sub(2), Constant((-compression)), boundaries, 2)
    bc2 = DirichletBC(W.sub(0).sub(0), Constant((0)), boundaries, 4)
    bc3 = DirichletBC(W.sub(0).sub(1), Constant((0)), boundaries, 4)

    bc4 = DirichletBC(W.sub(1), Constant((0)), boundaries, 4)
    
    # for more than one variants
    #bc5 = DirichletBC(W.sub(2), Constant((0)), boundaries, 4)
    #bc6 = DirichletBC(W.sub(3), Constant((0)), boundaries, 4)
    #bc7 = DirichletBC(W.sub(4), Constant((0)), boundaries, 4)
    #bc8 = DirichletBC(W.sub(5), Constant((0)), boundaries, 4)
    #bc9 = DirichletBC(W.sub(6), Constant((0)), boundaries, 4)
    
    # bcs = [bc1,bc2,bc3]
    # bcs = [bc1,bc2,bc3, bc4, bc5, bc6, bc7, bc8, bc9]
    bcs = [bc1,bc2,bc3, bc4] #, bc5, bc6, bc7, bc8, bc9]


    dx = Measure('dx', subdomain_data=cells)
    ds = Measure('ds', subdomain_data=boundaries)

    return dx, ds, boundaries, bcs


def cuboid_boundaries(mesh, compression, W):

    x = SpatialCoordinate(mesh)

    tol = DOLFIN_EPS

    bLEFT = CompiledSubDomain("near(x[0], 0.0, tol) && on_boundary", tol = tol)
    bRIGHT = CompiledSubDomain("near(x[0], 200.0, tol) && on_boundary", tol = tol)


    boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
    boundaries.set_all(0)
    
    print("Markers for all boundary facets are set")

    bLEFT.mark(boundaries, 2)       # left
    bRIGHT.mark(boundaries, 4)      # right
    

    #print(boundaries.array()[boundaries.array() == 2])
    print(boundaries.array()[boundaries.array() == 2])
    print(boundaries.array()[boundaries.array() == 4])



    cells = MeshFunction("size_t", mesh, mesh.topology().dim())
    cells.set_all(1)

    print("Marker for all internal cells is set to 1")

    bc1 = DirichletBC(W.sub(0).sub(0), Constant((compression)), boundaries, 4)
    bc2 = DirichletBC(W.sub(0).sub(1), Constant((0)), boundaries, 4)
    bc3 = DirichletBC(W.sub(0).sub(2), Constant((0)), boundaries, 4)

    bc4 = DirichletBC(W.sub(0).sub(0), Constant((0)), boundaries, 2)
    bc5 = DirichletBC(W.sub(0).sub(1), Constant((0)), boundaries, 2)
    bc6 = DirichletBC(W.sub(0).sub(2), Constant((0)), boundaries, 2)


    bcs = [bc1,bc2,bc3, bc4, bc5, bc6]


    dx = Measure('dx', subdomain_data=cells)
    ds = Measure('ds', subdomain_data=boundaries)

    return dx, ds, boundaries, bcs


def constitutive_law(mesh, num_variants, R, W, w, dw, wx, compression, penalty, time_step, strain_approx=None):

    #set_log_level(10)

    # Form compiler options
    parameters["form_compiler"]["optimize"]     = True
    parameters["form_compiler"]["cpp_optimize"] = True

    parameters.parse()



    from db_compliance import Compliance

    compliance_phase1, compliance_phase2 = Compliance.CuAlNi()

    Compliance_tetragonal = as_tensor(compliance_phase1)
    Compliance_monoclinic = as_tensor(compliance_phase2)

    # print(Compliance_monoclinic)
    # print(Compliance_tetragonal)


    # Definition of normal vectors to the elastic body and to the rigid half-space
    nRIGID = Constant((0.0, 0.0, 1.0))
    n = FacetNormal(mesh)




    ##################################
    # Phase field definitions
    ##################################


    def mar(c):

        # tot_mar = sum(c)
        
        # for one variant
        tot_mar = c

        return tot_mar
    

    def aus(c):
        return 1.0 - mar(c)


    def trans_strain(c):

        # strain = [c[i] * shape_strains_list[i] for i in range(num_variants)]
        #tr_strain = sum(strain)

        # for one variant
        strain = c * shape_strains_list

        #return tr_strain
        return strain      
  


    def mixture_compliance(c):
        return (aus(c) * Compliance_tetragonal) + (mar(c) * Compliance_monoclinic)



    ##################################
    # Displacement field definitions
    ##################################


    # matrix-to-vector transform
    def strain3voigt(e):
        return as_vector(
            [e[0, 0], 
            e[1, 1], 
            e[2, 2], 
            e[1, 2] * 2, 
            e[0, 2] * 2, 
            e[0, 1] * 2]
                        )

    # vector-to-matrix transform
    def voigt3stress(s):
        return as_tensor(
            [[s[0], s[5], s[4]], 
            [s[5], s[1], s[3]], 
            [s[4], s[3], s[2]]]
                        )


    def epsilon(u):

        if strain_approx is None:
            e = sym(grad(u))
            return e 
        
        elif strain_approx =='lagrangian':
            I = Identity(3)
            Fgrad = I + grad(u)
            C = Fgrad.T * Fgrad
            E = 0.5*(C - I)
            return E
        
        elif strain_approx =='hencky_order_22':
            I = Identity(3)
            Fgrad = I + grad(u)
            C = Fgrad.T * Fgrad
            a = C*C - I
            b = (C*C) + 4*C + I
            H = 1.5 * a * inv(b)
            return H

    def epsilon_elastic(u, c):
        #u, c = split(w)
        return epsilon(u) - trans_strain(c)
  
    def Vstrain(u, c):
        return strain3voigt(epsilon_elastic(u, c))

    def Vsigma(u, c):
        #u, c = split(w)
        return dot(mixture_compliance(c), Vstrain(u, c))

    def sigma(u, c):
        return voigt3stress(Vsigma(u, c))

    def dilatation(c):
        return det(trans_strain(c))


    def gap(u): # Definition of gap function
        x = SpatialCoordinate(mesh)
        return x[2] + u[2]

    def maculay(x): # Definition of Maculay bracket
        return (x+abs(x))/2
    
    def pN(u, c):
        return dot(dot(nRIGID, sigma(u, c)), nRIGID) 

    def convert_to_integer(a):
        "Convert to a 32-bit int. Raise exception if this will cause an overflow"
        if abs(a) <= iinfo(int32).max:
            return int32(a)
        else:
            raise OverflowError("Cannot safely convert to a 32-bit int.")


    def psi(w):

        u = split(w)[0]
        c = split(w)[1]

        #c = split(w)[1:]
        

        def elastic(u, c):
            elastic_energy = 0.5 * dilatation(c) * inner(sigma(u, c), epsilon_elastic(u, c))
            
            # approximate energy without volume dilatation (det F = 1)
            # elastic_energy = 0.5 * inner(sigma(u, c), epsilon_elastic(u, c))
            
            return elastic_energy
    
        def interface(c):

            AM_surface_energy = 0.0020     # J/m^2 = Pa.m = Pa'.nm (no change in magnitude in new scaling)
            AM_length_scale   = 50.0     # nm 

            MM_surface_energy = 0.00     # J/m^2 = Pa.m (no change in magnitude in new scaling)
            MM_length_scale   = 1.00     # nm


            AM_interface = (4.0/DOLFIN_PI) * (AM_surface_energy/AM_length_scale) * (aus(c)*mar(c))

            #MM_interface = (4.0/DOLFIN_PI) * (MM_surface_energy/MM_length_scale) * diffuse

            energy_cost = AM_interface #+ MM_interface
            
            return energy_cost

        def chemical(c):
            
            phi = 5.0e-3  # MPa = 1e-3 * Pa'
            chemical_energy = phi*mar(c)
            
            return chemical_energy


        def contact_pen(u):
            pen = (0.5 * penalty *inner(maculay(-gap(u)),maculay(-gap(u))))
            return pen

        def phase_pen(c):


            # penp1 = [inner(maculay(c[i] - 1), maculay(c[i] - 1)) for i in range(num_variants)]
            # penp0 = [inner(maculay(-c[i]), maculay(-c[i])) for i in range(num_variants)]

            #penp = sum(penp1) + sum(penp0)


            # for one variant
            penp1 = inner(maculay(c - 1), maculay(c - 1))
            penp0 = inner(maculay(-c), maculay(-c))

            penp = penp1 + penp0
                        
            return 0.5 * (10*penalty) * penp


        # return elastic(u, c), contact_pen(u), phase_pen(c)
        return elastic(u, c) + interface(c) + chemical(c), contact_pen(u), phase_pen(c)




    time_end = 1.0
    t_curr = 0.0

    t_curr = t_curr + time_step
    d_step = compression * t_curr




    while t_curr <= time_end:


        dx, ds, boundaries, bcs = sphere_boundaries(mesh, meshSIZE, R, d_step, W)

        print("Current displacement at boundary = ", d_step, " nm")



        # Total potential energy (linear elasticity, gradient-interface)
        Pi = psi(w)[0]*dx + psi(w)[1]*ds(6) + psi(w)[2]*dx

        # Compute first variation of Pi (directional derivative about u in the direction of v)
        F = derivative(Pi, w, wx)

        # Compute Jacobian of F for iterative "newton_solver"
        J = derivative(F, w, dw)



        wmin, wmax = func_space_bounds(W)

        for bc in bcs:
            bc.apply(wmin.vector())
            bc.apply(wmax.vector())


        # initial guess
        w.vector()[:] = d_step




        problem = NonlinearVariationalProblem(F, w, bcs, J=J)
        problem.set_bounds(wmin, wmax)
        solver  = NonlinearVariationalSolver(problem)


        snes_solver_parameters = {"nonlinear_solver": "snes",              
                                "snes_solver": {"linear_solver": 'lu', #"lu",
                                                 "preconditioner": "petsc_amg", #"bjacobi", 
                                                 "maximum_iterations": 10,
                                                 "absolute_tolerance": 1E-3,
                                                 "relative_tolerance": 1E-3,
                                                 "report": True,
                                                 "method":'vinewtonrsls',
                                                 "error_on_nonconvergence": True,
                                                 }}
        
        solver.parameters.update(snes_solver_parameters)
        

        #Solve the nonlinear problem
        iter, converged = solver.solve()



        class PhaseProblem(OptimisationProblem):

            def __init__(self):
                OptimisationProblem.__init__(self)

            # Objective function
            def f(self, x):
                w.vector()[:] = x
                return assemble(Pi)

            # Gradient of the objective function
            def F(self, b, x):
                w.vector()[:] = x
                assemble(F, tensor=b)

            # Hessian of the objective function
            def J(self, A, x):
                w.vector()[:] = x
                assemble(J, tensor=A)


        # solver = PETScTAOSolver()  

        # tao_solver_parameters = {"method": "bqpip",
        #                         #"line_search": "gpcg",
        #                         #"linear_solver": "lu",
        #                         #"preconditioner": "ilu",
        #                         "gradient_absolute_tol": 1e-3,
        #                         "gradient_relative_tol": 1e-3,
        #                         "maximum_iterations": 150,
        #                         "monitor_convergence": True,
        #                         "report": True
        # }

        # solver.parameters.update(tao_solver_parameters)

        # info(solver.parameters, True)

        # iter, converged = solver.solve(PhaseProblem(), w.vector(), wmin.vector(), wmax.vector())




        t_curr = t_curr + time_step
        d_step = compression * t_curr

    
    
        _u = w.split(deepcopy=True)[0]
        _c = w.split(deepcopy=True)[1]
        
        
        print("c_max = ", max(_c.vector()[:]), "c_min = ", min(_c.vector()[:]))

        print("u_max = ", max(_u.vector()[:]), "u_min = ", min(_u.vector()[:]))

        #from vedo.dolfin import plot
        #plot(_u.sub(2)) 






    return w, _u, _c, mar(_c), boundaries

    

   

    
 
    
hex = 0


R = 100  # nm
mesh_size = 5.0


mesh_folder = "./meshes/"
if hex == 1:
    subfolder = 'hexahedral/'
else:
    subfolder = "tetrahedral/"

mesh_folder = mesh_folder + subfolder
mesh_file = "mesh3D_R_" + str(int(R)) + "_size_" + str(int(mesh_size)) + ".xdmf"


# density = 25

# coordinates = (200.0, 100.0, 100.0)
# num_points = (density, density, density)

#mymesh, hMESH, meshSIZE, = cuboid_geometry((coordinates, num_points))


mymesh, hMESH, meshSIZE, = sphere_geometry(mesh_folder, mesh_file)


basis = 'linear'
W, W1, W2, w, dw, wx, num_variants, shape_strains_list = func_spaces(mymesh, basis)


compression = R/100.0

#dx, ds, boundaries, bcs = sphere_boundaries(mymesh, meshSIZE, R, compression, W)


penalty = 100
#E = 200  # GPa = Pa'

time_step = 5e-1

print("hMESH = ", hMESH)




_w, _u, _c, mar, boundaries = constitutive_law(mymesh, num_variants, R, W, w, dw, wx, compression, penalty, time_step)


from vedo.dolfin import plot
plot(_u)



file = File("./displacement-pen-phase.pvd")
file << project(_u, W1)

print(mar)

#file = File("./phase-pen-phase.pvd")
#file << project(mar, W2)




#w, boundaries = pureElastic_constitutive_law(mymesh, meshSIZE, hMESH, R, W, w, dw, v, compression, penalty, E, time_step)

#from vedo.dolfin import plot

#plot(w)

#with XDMFFile("u_PEN_" + ".xdmf") as xdmf:
#        xdmf.write_checkpoint(w, "u", 0.0, XDMFFile.Encoding.HDF5, append=False)


# Savu u as .pvd 
# file = File("./displacement-pen.pvd")
# file << w
2

Skimming through the code, it looks like you’re using a direct solver (LU decomposition) in each Newton step, which is usually insensitive to conditioning of the tangent matrix, except in really extreme cases. More likely is that the time/load steps in the outer loop are too large, so the initial guess of the Newton iteration is not close enough to the solution for it to converge. Some common strategies would be to first try a line search method, e.g., adding the key/value pair

"line_search":"bt"

to the value for key "snes_solver" in the dictionary snes_solver_parameters (which could be written on its own line as

snes_solver_parameters["snes_solver"]["line_search"] = "bt"

if the nesting of dictionaries is unclear from my wording above), then, if that doesn’t work, try shrinking the time/load step size. However, debugging nonlinear convergence can sometimes involve a lot of problem-specific troubleshooting.

3 Likes
3

Thanks for your response. I am now getting convergence by adaptively setting the time-step. Nonlinearity seems to be the likely cause behind this. I have been asked to implement hexahedral mesh to reduce the necessity of fine refinement and stability.

I generate both hexahedron-element and tetra-element mesh using pygmsh. However, using the hex mesh in the same above code makes dolfin crash.

FYI the MWE code to generate hex elements is below:

mesh_size = 10.0
R = 50.0

import pygmsh, gmsh
import meshio
from math import sqrt

for hex in range(2):

    geom = pygmsh.occ.Geometry()

    model3D = geom.__enter__()

    if hex == 1:

        algo2D = 8 # Frontal-Delaunay for Quads
        algo3D = 1 # Delaunay
        algoRecom = 0 # 1. Blossom, 3. Blossom full-quad
        algoSubDiv = 2 # All hexa

        gmsh.option.setNumber("Mesh.Algorithm", algo2D)
        gmsh.option.setNumber("Mesh.Algorithm3D", algo3D)
        gmsh.option.setNumber("Mesh.RecombinationAlgorithm", algoRecom)
        gmsh.option.setNumber("Mesh.RecombineAll", 1)
        gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", algoSubDiv)


    #model3D.characteristic_length_min = mesh_size
    model3D.characteristic_length_max = mesh_size

    box =  model3D.add_box(x0=(-R, -R, R), extents=(2*R, 2*R, 1.1*R))
    #box1 =  model3D.add_box([0.5, 0.5, 1], [0.5, 0.5, 1])
    ball = model3D.add_ball([0.0, 0.0, R], R)

    union_minus_ball = model3D.boolean_difference([ball], [box])

    model3D.synchronize()
    model3D.add_physical(union_minus_ball, "union_minus_ball")

    gmesh = model3D.generate_mesh()

    if hex == 1:
        gmesh.write("./meshes/hexahedral/mesh3D_R_" + str(int(R)) + "_size_" + str(int(mesh_size)) + ".vtu")
    else:
        gmesh.write("./meshes/tetrahedral/mesh3D_R_" + str(int(R)) + "_size_" + str(int(mesh_size)) + ".vtu")




    points, cells, cell_data_dict = gmesh.points, gmesh.cells, gmesh.cell_data_dict

    print(cells)
    print(cell_data_dict)

    import numpy as np

    tetra_cells = None
    hexa_cells = None

    for cell in cells:
        
        print(cell.type)
        
        if cell.type == 'hexahedron':
            if hexa_cells is None:
                hexa_cells = cell.data
            else:
                hexa_cells = np.vstack((hexa_cells, cell.data))
            
        if cell.type == "tetra":
            # Collect the individual meshes
            if tetra_cells is None:
                tetra_cells = cell.data
            else:
                tetra_cells = np.vstack((tetra_cells, cell.data))

    tetra_data = None
    hexa_data = None

    for key in cell_data_dict["union_minus_ball"].keys():
        print(key)

        if key == 'hexahedron':
            hexa_data = cell_data_dict["union_minus_ball"][key]
            hexa_mesh = meshio.Mesh(points=points, cells={"hexahedron": hexa_cells},
                                    cell_data={"name_to_read":[hexa_data]})
        
        if key == "tetra":
            tetra_data = cell_data_dict["union_minus_ball"][key]
        
            # Create tetra mesh with cell data attached
            tetra_mesh = meshio.Mesh(points=points, cells={"tetra": tetra_cells},
                                    cell_data={"name_to_read":[tetra_data]})

    if hex == 1:
        meshio.write("./meshes/hexahedral/mesh3D_R_" + str(int(R)) + "_size_" + str(int(mesh_size)) + ".xdmf", hexa_mesh)

        print("Hex meshing DONE, Radius = ", R, " max_edge_length = ", mesh_size, "\n")

    else:
        meshio.write("./meshes/tetrahedral/mesh3D_R_" + str(int(R)) + "_size_" + str(int(mesh_size)) + ".xdmf", tetra_mesh)

        print("Tet meshing DONE, Radius = ", R, " max_edge_length = ", mesh_size, "\n")

Where could the issue be? in my meshing or in dolphin’s support?

Appreciate your help! Thanks.

4

The legacy dolfin code has very limited support for hex meshes. This has been vastly improved in DOLFINx

5

Thanks @dokken. Will shift to DOLFINx in a short period of time.

The above code generates a correct tetrahedral mesh, however when I open the xdmf file of hex mesh, paraview presents element type error. Not sure what is the issue here.

I have written this code referring the method explained on this link