Nonlinear hyperelasticity in dolfin-x: Minimal example, and a few issues

dolfinx
1

Dear community,

While transferring my scripts from dolfin to dolfin-x, a couple of small issues arose for which I haven’t been able to find a satisfactory answer in the demos and the forums yet.
Therefore, I set up a minimal nonlinear hyperelasticity example which works fine in dolfin but has a few issues in dolfin-x.

The example can either be loaded with a uniaxial Dirichlet value or with a uniaxial Neumann load, linearly ramped up from 0 to T (see header of files).

I decided to share the examples (and the mesh files input_domain.xdmf and input_boundary.xdmf) which should work straight away when adapting the basepath variable:

The dolfin-x version:

#!/usr/bin/env python3

# This is a minimal nonlinear quasi-static hyperelasticity example
# with a Neohookean block that is uniaxially loaded, either with a
# prescribed Dirichlet condition (if loading = 'dirichlet')
# or with a prescribed true Neumann traction (if loading = 'neumann');
# loads are ramped up linearly with time t \in [0,1];
# the block is constrained such that a uniaxial stress state (in x-dir)
# is obtained

import time
import sys, os, subprocess, time
import math

import numpy as np
from mpi4py import MPI
from petsc4py import PETSc

from dolfinx import *
from dolfinx.mesh import *
from dolfinx.cpp.mesh import *
from dolfinx.fem import locate_dofs_topological, locate_dofs_geometrical, assemble_matrix, assemble_vector, set_bc, apply_lifting
from dolfinx.io import *
from dolfinx.la import *
from ufl import *

start = time.time()

# base path for reading input from folder in and writing output to folder out
basepath  = '/home/mh/work/fem_scripts/testing'

loading = 'dirichlet' # dirichlet, neumann

T       = 1.0
nstep   = 10
dt      = T/nstep

# in which interval to write results
write_results_every = 1

# encoding and mpi communicator
encoding = dolfinx.cpp.io.XDMFFile.Encoding.ASCII
comm = MPI.COMM_WORLD

# read in xdmf mesh - domain
with XDMFFile(comm,''+basepath+'/in/input_domain.xdmf','r',encoding=encoding) as infile:
    mesh = infile.read_mesh("Grid")

# read in xdmf mesh - boundary
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
with XDMFFile(comm,''+basepath+'/in/input_boundary.xdmf','r',encoding=encoding) as infile:
    mt_s = infile.read_meshtags(mesh, "Grid_surf")

# reference normal on mesh facets (for true Neumann load, or other applications)
n0 = FacetNormal(mesh)

# function spaces for displacement, stresses, strains
# CG: Lagrange, DG: discontinuous Lagrange
V_u          = VectorFunctionSpace(mesh, ("CG", 1))
V_strain     = TensorFunctionSpace(mesh, ("DG", 0))
V_stress     = TensorFunctionSpace(mesh, ("DG", 0))

### define functions
du    = TrialFunction(V_u)            # Incremental displacement
var_u = TestFunction(V_u)             # Test function
u     = Function(V_u, name='Displacement')
# disp increment for Newton
del_u = Function(V_u)
# strain and stress functions for postprocessing
glstrain     = Function(V_strain, name="GreenStrain")
cauchystress = Function(V_stress, name="CauchyStress")

# print infos on problem size
if comm.rank == 0:
    print("Number of degrees of freedom: %i" % (V_u.dim()))
    sys.stdout.flush()


# a zero function
zero = Function(V_u)
zero.vector.set(0.0)

# expression for time-dependent Dirichlet condition
class u_prescr_expr:
    def __init__(self):
        self.t = 0.0

    def eval(self, x):
        u_uni = 0.1
        # Added some spatial/temporal variation here
        return np.full(x.shape[1], u_uni*self.t/T)

u_prescr = u_prescr_expr()
u_prescr.t = 0.0

# interpolate expression into function
u_prescr_func = Function(V_u)
u_prescr_func.interpolate(u_prescr.eval)

# function for time-dependent true follower Neumann load
def trueload(t):
    return Constant(mesh, 10.0*t/T)

# Boundary conditions
dbcs_hom=[]
dbcs_inhom=[]

dbcs_hom.append( DirichletBC(zero, locate_dofs_topological(V_u.sub(0), 2, mt_s.indices[mt_s.values == 1])) ) # x=0
dbcs_hom.append( DirichletBC(zero, locate_dofs_topological(V_u.sub(1), 2, mt_s.indices[mt_s.values == 2])) ) # y=0
dbcs_hom.append( DirichletBC(zero, locate_dofs_topological(V_u.sub(2), 2, mt_s.indices[mt_s.values == 3])) ) # z=0

if loading=='dirichlet':

    dbcs_inhom.append( DirichletBC(u_prescr_func, locate_dofs_topological(V_u.sub(0), 2, mt_s.indices[mt_s.values == 4])) ) # x-load Dirichlet
    dbcs_hom.append( DirichletBC(zero, locate_dofs_topological(V_u.sub(0), 2, mt_s.indices[mt_s.values == 4])) ) # homgenized
    
elif loading=='neumann':

    ds4 = ds(subdomain_data=mt_s, subdomain_id=4) # x-load Neumann
    
else:
    print("Set loading to either dirichlet or neumann!")
    sys.exit()


### nonlinear finite strain kinematics
I = Identity(len(u))  # Identity tensor

F = I + nabla_grad(u) # deformation gradient
J = det(F)            # determinant of deformation gradient
# use variable expression since later we diff w.r.t these
C = variable(F.T*F)   # right Cauchy-Green tensor

# Cauchy-Green invariants
Ic   = tr(C)
IIc  = 0.5*(tr(C)**2. - tr(C*C))
IIIc = det(C)
# isochoric Cauchy-Green invariants
Ic_bar   = IIIc**(-1./3.) * Ic
IIc_bar  = IIIc**(-2./3.) * IIc


### material law - Neohookean, isochoric-volumetric split
mu, nu = 10., 0.49
kappa = mu/(1.-2.*nu)
# isochoric strain-energy function
psi_iso = (mu/2.) * (Ic_bar - 3.) # isochoric strain energy function
# volumetric strain-energy function
psi_vol = (kappa/4.) * (IIIc - 2.*ln(sqrt(IIIc)) - 1.) # Ogden-type volumetric part
# total strain-energy function
psi = psi_iso + psi_vol


#2nd Piola-Kirchhoff stress 
def PK2_stress(r):
    S = 2.*diff(psi,C)
    return S

# Green-Lagrange strain (only for post-processing)
def GL_strain(r):
    F = I + grad(r)
    E = 0.5*(F.T*F - I)    
    return E

#Cauchy stress (only for post-processing)
def Cauchy_stress(r):
    F = I + grad(r)
    J = det(F)
    S = PK2_stress(r)
    sigma = (1./J) * F*S*F.T
    return sigma


# internal virtual work
def dW_int(u_):
    
    return 0.5*inner(PK2_stress(u_),derivative(C, u_, var_u))*dx


# external virtual work
def dW_ext(u_):
    
    if loading=='dirichlet':
        
        return 0.0
    
    elif loading=='neumann':
        
        return trueload(t)*J*dot(dot(inv(F).T,n0), var_u)*ds4
    
    else:
        print("Set loading to either dirichlet or neumann!")
        sys.exit()


### LocalSolver seems to have been removed in Dolfin-x
def local_project(var_u, V, u=None):
    """Element-wise projection using LocalSolver"""
    du_    = TrialFunction(V)
    var_u_ = TestFunction(V)
    a_proj = inner(du_, var_u_)*dx
    b_proj = inner(var_u, var_u_)*dx
    solver = LocalSolver(a_proj, b_proj)
    solver.factorize()
    if u is None:
        u = Function(V)
        solver.solve_local_rhs(u)
        return u
    else:
        solver.solve_local_rhs(u)
        return

results_xdmf_file = XDMFFile(comm,''+basepath+'/out/results.xdmf','w',encoding=encoding)
results_xdmf_file.write_mesh(mesh)


# load/time stepping
interval = np.linspace(0, T, nstep+1)

for (N, dt) in enumerate(np.diff(interval)):
    
    t = interval[N+1]

    # update time-dependent expressions
    u_prescr.t = t
    u_prescr_func.interpolate(u_prescr.eval)
   
    # nonlinear variational form: internal minus external virtual work
    varform = dW_int(u) - dW_ext(u)

    # Jacobian of nonlinear variational form
    jac = derivative(varform, u, du)

    # Newton loop
    struct_res_norm = 1.0
    struct_inc_norm = 1.0
    tol_struct_res = 1.0e-8
    tol_struct_inc = 1.0e-8
    it = 0
    maxiter = 25
    
    # apply inhomogeneous bcs to solution vector
    set_bc(u.vector, dbcs_inhom)
    
    K = assemble_matrix(jac, dbcs_hom)
    K.assemble()
    
    r = assemble_vector(-varform)
    
    #apply_lifting(r, [-varform], [dbcs_hom])
    r.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
    set_bc(r, dbcs_hom)

    struct_disp_norm_0 = u.vector.norm()
    struct_res_norm_0 = r.norm()
    
    
    if comm.rank == 0:
        print("Predictor (disp 2-norm %.4e) yields residual 2-norm:\n      %.4e" % (struct_disp_norm_0,struct_res_norm_0))
        print("iter  struct res 2-norm  struct incr 2-norm")
        sys.stdout.flush()

    while it < maxiter:
        
        it += 1
        
        # solve linearized system
        solve(K, del_u.vector, r)
    
        # update solution
        u.vector.axpy(1.0, del_u.vector)
        
        K = assemble_matrix(jac, dbcs_hom)
        K.assemble()
        
        r = assemble_vector(-varform)
        #apply_lifting(r, [-varform], [dbcs])
        r.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        set_bc(r, dbcs_hom)
        
        struct_res_norm = r.norm()
        struct_inc_norm = del_u.vector.norm()
        
        if comm.rank == 0:
            print("%i     %.4e         %.4e" % (it, struct_res_norm, struct_inc_norm))
            sys.stdout.flush()
        
        # check if converged
        if struct_res_norm <= tol_struct_res and struct_inc_norm <= tol_struct_inc:
            break
    
    else:

        if comm.rank == 0:
            print("Newton did not converge after %i iterations!" % (it))
            sys.stdout.flush()
        
        sys.exit()    


    # write results every write_results_every steps
    if (N+1) % write_results_every == 0:
        
        # Save solution to XDMF format
        #results_xdmf_file.write(u, t)
        results_xdmf_file.write_function(u, t)
        # Compute strains and save to file
        #local_project(GL_strain(u), V_strain, glstrain)
        ##results_xdmf_file.write(glstrain, t)
        #results_xdmf_file.write_function(glstrain, t)
        ## Compute Cauchy stresses and save to file
        #local_project(Cauchy_stress(u), V_stress, cauchystress)
        ##results_xdmf_file.write(cauchystress, t)
        #results_xdmf_file.write_function(cauchystress, t)
        

    if comm.rank == 0: # only proc 0 should print this
        print("### TIME STEP %i / %i successfully completed, TIME: %.4f" % (N+1,nstep,t))
        print("-----------------------------------------------------------")
        sys.stdout.flush()
        

if comm.rank == 0: # only proc 0 should print this
    print('Time for computation: %.1f min' % ( (time.time()-start)/60. ))
    sys.stdout.flush()

The dolfin version (latest development dolfin):

#!/usr/bin/env python3

# This is a minimal nonlinear quasi-static hyperelasticity example
# with a Neohookean block that is uniaxially loaded, either with a
# prescribed Dirichlet condition (if loading = 'dirichlet')
# or with a prescribed true Neumann traction (if loading = 'neumann');
# loads are ramped up linearly with time t \in [0,1];
# the block is constrained such that a uniaxial stress state (in x-dir)
# is obtained

import time
import sys, os, subprocess, time
import math

import matplotlib.pyplot as plt
from dolfin import *
import numpy as np

start = time.time()

# base path for reading input from folder in and writing output to folder out
basepath  = '/home/mh/work/fem_scripts/testing'

loading = 'dirichlet' # dirichlet, neumann

# time/load stepping
T       = 1.0
nstep   = 10
dt      = T/nstep

# in which interval to write results
write_results_every = 1

# mpi communicator
comm = MPI.comm_world

# read in xdmf mesh - domain
mesh = Mesh()
with XDMFFile(''+basepath+'/in/input_domain.xdmf') as infile:
    infile.read(mesh)

# read in xdmf mesh - boundary
mvc_s = MeshValueCollection("size_t", mesh, 2)
with XDMFFile(''+basepath+'/in/input_boundary.xdmf') as infile:
    infile.read(mvc_s, "boundaries_surf")

mf_s = cpp.mesh.MeshFunctionSizet(mesh, mvc_s)

# reference normal on mesh facets (for true Neumann load, or other applications)
n0 = FacetNormal(mesh)

# function spaces for displacement, stresses, strains
# CG: Lagrange, DG: discontinuous Lagrange
V_u          = VectorFunctionSpace(mesh, "CG", 1)
V_strain     = TensorFunctionSpace(mesh, "DG", 0)
V_stress     = TensorFunctionSpace(mesh, "DG", 0)

### define functions
du    = TrialFunction(V_u)            # Incremental displacement
var_u = TestFunction(V_u)             # Test function
u     = Function(V_u, name='Displacement')
# disp increment for Newton
del_u = Function(V_u)
# strain and stress functions for postprocessing
glstrain     = Function(V_strain, name="GreenStrain")
cauchystress = Function(V_stress, name="CauchyStress")

# print infos on problem size
if comm.rank == 0:
    print("Number of degrees of freedom: %i" % (V_u.dim()))
    sys.stdout.flush()


# a zero constant
zeroval = Constant(0.0)

# expression for time-dependent Dirichlet condition
u_prescr_expr = Expression(('u_uni*t/'+str(T)+''), u_uni=0.1, t=0, degree=0)

# function for time-dependent true follower Neumann load
def trueload(t):
    return Constant(10.0*t/T)

# Boundary conditions
dbcs_hom=[]
dbcs_inhom=[]

dbcs_hom.append( DirichletBC(V_u.sub(0), zeroval, mf_s, 1) ) # x=0
dbcs_hom.append( DirichletBC(V_u.sub(1), zeroval, mf_s, 2) ) # y=0
dbcs_hom.append( DirichletBC(V_u.sub(2), zeroval, mf_s, 3) ) # z=0

if loading=='dirichlet':

    dbcs_inhom.append( DirichletBC(V_u.sub(0), u_prescr_expr, mf_s, 4) ) # x-load Dirichlet
    dbcs_hom.append( DirichletBC(V_u.sub(0), zeroval, mf_s, 4) ) # homogenized
    
elif loading=='neumann':

    ds4 = ds(subdomain_data=mf_s, subdomain_id=4) # x-load Neumann
    
else:
    print("Set loading to either dirichlet or neumann!")
    sys.exit()


### nonlinear finite strain kinematics
I = Identity(len(u))  # Identity tensor

F = I + nabla_grad(u) # deformation gradient
J = det(F)            # determinant of deformation gradient
# use variable expression since later we diff w.r.t these
C = variable(F.T*F)   # right Cauchy-Green tensor

# Cauchy-Green invariants
Ic   = tr(C)
IIc  = 0.5*(tr(C)**2. - tr(C*C))
IIIc = det(C)
# isochoric Cauchy-Green invariants
Ic_bar   = IIIc**(-1./3.) * Ic
IIc_bar  = IIIc**(-2./3.) * IIc


### material law - Neohookean, isochoric-volumetric split
mu, nu = 10., 0.49
kappa = mu/(1.-2.*nu)
# isochoric strain-energy function
psi_iso = (mu/2.) * (Ic_bar - 3.) # isochoric strain energy function
# volumetric strain-energy function
psi_vol = (kappa/4.) * (IIIc - 2.*ln(sqrt(IIIc)) - 1.) # Ogden-type volumetric part
# total strain-energy function
psi = psi_iso + psi_vol


#2nd Piola-Kirchhoff stress 
def PK2_stress(r):
    S = 2.*diff(psi,C)
    return S

# Green-Lagrange strain (only for post-processing)
def GL_strain(r):
    F = I + grad(r)
    E = 0.5*(F.T*F - I)    
    return E

#Cauchy stress (only for post-processing)
def Cauchy_stress(r):
    F = I + grad(r)
    J = det(F)
    S = PK2_stress(r)
    sigma = (1./J) * F*S*F.T
    return sigma


# internal virtual work
def dW_int(u_):
    
    return 0.5*inner(PK2_stress(u_),derivative(C, u_, var_u))*dx


# external virtual work
def dW_ext(u_):
    
    if loading=='dirichlet':
        
        return 0.0
    
    elif loading=='neumann':
        
        return trueload(t)*J*dot(dot(inv(F).T,n0), var_u)*ds4
    
    else:
        print("Set loading to either dirichlet or neumann!")
        sys.exit()
    

def local_project(var_u, V, u=None):
    """Element-wise projection using LocalSolver"""
    du_    = TrialFunction(V)
    var_u_ = TestFunction(V)
    a_proj = inner(du_, var_u_)*dx
    b_proj = inner(var_u, var_u_)*dx
    solver = LocalSolver(a_proj, b_proj)
    solver.factorize()
    if u is None:
        u = Function(V)
        solver.solve_local_rhs(u)
        return u
    else:
        solver.solve_local_rhs(u)
        return


results_xdmf_file = XDMFFile(''+basepath+'/out/results.xdmf')
results_xdmf_file.parameters["flush_output"] = True
results_xdmf_file.parameters["functions_share_mesh"] = True
results_xdmf_file.parameters["rewrite_function_mesh"] = False


# load/time stepping
interval = np.linspace(0, T, nstep+1)

for (N, dt) in enumerate(np.diff(interval)):
    
    t = interval[N+1]

    # update time-dependent expressions
    u_prescr_expr.t = t

    # nonlinear variational form: internal minus external virtual work
    varform = dW_int(u) - dW_ext(u)

    # Jacobian of nonlinear variational form
    jac = derivative(varform, u, du)

    # Newton loop
    struct_res_norm = 1.0
    struct_inc_norm = 1.0
    tol_struct_res = 1.0e-8
    tol_struct_inc = 1.0e-8
    it = 0
    maxiter = 25

    # apply inhomogeneous bcs to solution vector
    for bc in dbcs_inhom: bc.apply(u.vector())
    
    # assemble system
    K, r = assemble_system(jac, -varform, dbcs_hom)

    # get initial norms
    struct_disp_norm_0 = norm(u.vector(),'l2')
    struct_res_norm_0 = norm(r,'l2')
    

    if comm.rank == 0:
        print("Predictor (disp 2-norm %.4e) yields residual 2-norm:\n      %.4e" % (struct_disp_norm_0,struct_res_norm_0))
        print("iter  struct res 2-norm  struct incr 2-norm")
        sys.stdout.flush()
    
    while it < maxiter:
        
        it += 1
        
        # solve linearized system
        solve(K, del_u.vector(), r, 'superlu_dist', 'none')
        
        # update solution
        u.vector().axpy(1.0, del_u.vector())
        
        # assemble system
        K, r = assemble_system(jac, -varform, dbcs_hom)

        struct_res_norm = norm(r,'l2')
        struct_inc_norm = norm(del_u.vector(),'l2')

        if comm.rank == 0:
            print("%i     %.4e         %.4e" % (it, struct_res_norm, struct_inc_norm))
            sys.stdout.flush()
        
        # check if converged
        if struct_res_norm <= tol_struct_res and struct_inc_norm <= tol_struct_inc:
            break
    
    else:

        if comm.rank == 0:
            print("Newton did not converge after %i iterations!" % (it))
            sys.stdout.flush()
        
        sys.exit()


    # write results every write_results_every steps
    if (N+1) % write_results_every == 0:
        
        # Save solution to XDMF format
        results_xdmf_file.write(u, t)
        # Compute strains and save to file
        local_project(GL_strain(u), V_strain, glstrain)
        results_xdmf_file.write(glstrain, t)
        # Compute Cauchy stresses and save to file
        local_project(Cauchy_stress(u), V_stress, cauchystress)
        results_xdmf_file.write(cauchystress, t)


    if comm.rank == 0: # only proc 0 should print this
        print("### TIME STEP %i / %i successfully completed, TIME: %.4f" % (N+1,nstep,t))
        print("-----------------------------------------------------------")
        sys.stdout.flush()


if comm.rank == 0: # only proc 0 should print this
    print('Time for computation: %.1f min' % ( (time.time()-start)/60. ))
    sys.stdout.flush()

The mesh file input_domain.xdmf:

<Xdmf Version="3.0">
  <Domain>
    <Grid Name="Grid">
      <Geometry GeometryType="XYZ">
        <DataItem DataType="Float" Dimensions="14 3" Format="XML" Precision="8">
          0.0 0.0 1.0
          0.0 0.0 0.0
          0.0 1.0 1.0
          0.0 1.0 0.0
          1.0 0.0 1.0
          1.0 0.0 0.0
          1.0 1.0 1.0
          1.0 1.0 0.0
          0.0 0.5 0.5
          1.0 0.5 0.5
          0.5 0.0 0.5
          0.5 1.0 0.5
          0.5 0.5 0.0
          0.5 0.5 1.0
        </DataItem>
      </Geometry>
      <Topology NodesPerElement="4" NumberOfElements="24" TopologyType="tetrahedron">
        <DataItem DataType="Int" Dimensions="24 4" Format="XML" Precision="4">
          9 11 10 12
          8 13 11 10
          11 10 13 9
          8 10 11 12
          1 0 8 10
          0 2 8 13
          10 0 13 4
          3 11 8 2
          1 8 3 12
          11 13 2 6
          4 13 9 6
          6 11 9 7
          10 4 9 5
          11 7 3 12
          12 9 7 5
          12 1 10 5
          13 8 0 10
          8 2 11 13
          10 8 1 12
          11 3 8 12
          7 11 9 12
          13 9 10 4
          11 9 13 6
          12 10 9 5
        </DataItem>
      </Topology>
    </Grid>
  </Domain>
</Xdmf>

The mesh file input_boundary.xdmf:

<Xdmf Version="3.0">
  <Domain>
    <Grid Name="Grid_surf">
      <Topology NodesPerElement="3" NumberOfElements="24" TopologyType="triangle">
        <DataItem DataType="Int" Dimensions="24 3" Format="XML" Precision="4">
          1 0 8
          0 2 8
          3 1 8
          2 3 8
          0 1 10
          4 0 10
          1 5 10
          5 4 10
          1 3 12
          5 1 12
          3 7 12
          7 5 12
          5 9 4
          4 9 6
          7 9 5
          6 9 7
          2 11 3
          6 11 2
          3 11 7
          7 11 6
          0 13 2
          4 13 0
          2 13 6
          6 13 4
        </DataItem>
      </Topology>
      <Attribute AttributeType="Scalar" Center="Face" Name="boundaries_surf">
        <DataItem DataType="Int" Dimensions="24" Format="XML" Precision="4">
          1
          1
          1
          1
          2
          2
          2
          2
          3
          3
          3
          3
          4
          4
          4
          4
          5
          5
          5
          5
          6
          6
          6
          6
        </DataItem>
      </Attribute>
    </Grid>
    <Grid Name="Grid_edge">
      <Topology NodesPerElement="2" NumberOfElements="2" TopologyType="Polyline">
        <DataItem DataType="Int" Dimensions="2 2" Format="XML" Precision="4">
          1 3
          1 0
        </DataItem>
      </Topology>
      <Attribute AttributeType="Scalar" Center="Edge" Name="boundaries_edge">
        <DataItem DataType="Int" Dimensions="2" Format="XML" Precision="4">
          1
          2
        </DataItem>
      </Attribute>
    </Grid>
    <Grid Name="Grid_point">
      <Topology NodesPerElement="1" NumberOfElements="1" TopologyType="Polyvertex">
        <DataItem DataType="Int" Dimensions="1 1" Format="XML" Precision="4">
          1
        </DataItem>
      </Topology>
      <Attribute AttributeType="Scalar" Center="Node" Name="boundaries_point">
        <DataItem DataType="Int" Dimensions="1" Format="XML" Precision="4">
          1
        </DataItem>
      </Attribute>
    </Grid>
  </Domain>
</Xdmf>

The dolfin-x and dolfin versions yield the same results in serial. However, the main issue I’m having with the dolfin-x version is that it does not work in parallel, and I think it might have to do with how Dirichlet conditions are applied. From the demos, I deduced that a routine called
apply_lifting(r, [-varform], [dbcs_hom]) (line 250 dolfin-x example) has to be called after rhs assemble, however I had to comment this function out since it produced the following error on my end:

Traceback (most recent call last):
File “./dolfinx_minimal.py”, line 250, in
apply_lifting(r, [-varform], [dbcs_hom])
File “/home/mh/.local/lib/python3.6/site-packages/dolfinx/fem/assemble.py”, line 314, in apply_lifting
cpp.fem.apply_lifting(b, _create_cpp_form(a), bcs, x0, scale)
IndexError: vector::_M_range_check: __n (which is 1) >= this->size() (which is 1)

I have no idea why this is happening, and if this function is the cause for parallel non-functionality (since in serial it works and gives the same results as the dolfin code).
Am I missing something fundamental here?

A second issue is: When setting loading=‘neumann’ in my example, my time dependent function trueload(t) is not accepted by ufl, I get:

File “/home/mh/.local/lib/python3.6/site-packages/ufl/exproperators.py”, line 182, in _mul
return _mult(self, o)
File “/home/mh/.local/lib/python3.6/site-packages/ufl/exproperators.py”, line 114, in _mult
r1, r2 = len(s1), len(s2)
TypeError: object of type ‘numpy.float64’ has no len()

while in dolfin, this works fine.

A third smaller question would be on the LocalSover in the local_project function for projecting of derived quantities like stresses or strains for postprocessing. It seems that the LocalSolver has gone in dolfin-x, so I do not know how I would perform the local solve for my stresses and strains. Has this been renamed or replaced?

I know that I addressed rather different topics here, but since they all relate to the same minimal example, I decided to group these three issues since at least they all deal with how things translate from dolfin to dolfin-x.
Maybe, these minimal examples also help others how to tansfer and address specific problems in the latest dolfin-x version.

Thanks!

Best wishes,
Marc

2 Likes
2

As this is a rather lengthy example, I do not have time to go through it all.
One thing I note is that you are not applying bcs to the sub-spaces correctly.
See:


on how to apply bcs to subspaces.

Instead of LocalSolver, you can project your solution to an appropriate DG space.

2 Likes
3

Ok, thanks for your help. Sorry the example is a bit lengthy, however I could not shorten it further (a lot of length comes from the time/load stepping and the Newton scheme).

I changed my DBCs according to the demo. I collapsed the three subspaces
V_u0 = V_u.sub(0).collapse()
…
interpolated the prescribed functions (zeros and u_prescr)
u_prescr_func = Function(V_u0)

and changed all DBC lines
dbcs_inhom.append( DirichletBC(u_prescr_func, locate_dofs_topological((V_u.sub(0), V_u0), 2, mt_s.indices[mt_s.values == 4]),V_u.sub(0)) )

which yields the same serial results. The problem with not being able to call apply_lifting remains with the same error.

I searched for the project function (https://fenicsproject.org/docs/dolfin/1.6.0/python/programmers-reference/fem/projection/project.html) but it seems to have gone…
Is this replaced by the interpolate functionality?

Thanks!

Bw,
Marc

4

A projection is simply solving a problem with a mass matrix on the LHS, see:
https://bitbucket.org/fenics-project/dolfin/src/f0a90b8fffc958ed612484cd3465e59994b695f9/python/dolfin/fem/projection.py?at=master#lines-126:138
This can easily be implemented manually in dolfinx (or dolfin).

There are a lot of things you can do to shorten the example. As you cannot apply the lifting (even at the first time step), the time loop can go. You should also experiment with different linear forms, to see what distinguish you form from those used in the examples demos.

A minimal working example does not have to solve the complete problem, it should rather be the shortest code possible to get to the error message.

5

Thanks, I figured out how to apply the lifting and solved the parallel issue with DBCs.

Seems to me that the projection in this example assembles the whole system and then solves. There was this local solving in the elastodynamics example

def local_project(var_u, V, u=None):
    """Element-wise projection using LocalSolver"""
    du_    = TrialFunction(V)
    var_u_ = TestFunction(V)
    a_proj = inner(du_, var_u_)*dx
    b_proj = inner(var_u, var_u_)*dx
    solver = LocalSolver(a_proj, b_proj)
    solver.factorize()
    if u is None:
        u = Function(V)
        solver.solve_local_rhs(u)
        return u
    else:
        solver.solve_local_rhs(u)
        return

which efficiently solved at integration point level. Is there any chance to do this that efficiently in dolfinx or do I have to do another global assemble and solve?

Thx!
Marc

6

I cant say how things will progress in dolfinx. It might be worth raising an issue requesting the local solver feature, or at least someone will tell you if it is faster than the DG projection. For the sake of other readers, please share your solution to the apply lifting problem.

1 Like
7

Ok, I’ll do that!

Yes, of course I will share my solution. The system matrix and rhs assemble have to look like this:

# assemble system matrix
K = assemble_matrix(jac, dbcs)
K.assemble()

# assemble rhs vector
r = assemble_vector(varform)
apply_lifting(r, [jac], [dbcs], x0=[u.vector], scale=-1.0)
r.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(r, dbcs, x0=u.vector, scale=-1.0)

The apply_lifting and the way to call set_bc like this also work if there are non-homogeneous DBCs present. One does not have to manually create a homogeneous version or call dbc.homogenize.
Solving the linear system then should be called like this:

# solve linearized system
solve(K, del_u.vector, -r)

# update solution
u.vector.axpy(1.0, del_u.vector)
u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

where the ghostUpdate on u is essential for parallel functionality.

Bw,
Marc

8

@marchirschvogel Thank you so much for sharing! Your post was very helpful.

I know it is not the main issue of your post, but I was wondering, why do you use:

Instead of

F = I + grad(u)

as in the hyperelasticity demo? From my understanding, it should be I + grad(u).

9

Hi Alberto,

Yes, you are absolutely right! I’ve already bumped into this “mistake”: nabla_grad actually computes the gradient transposed.
So it does not matter if you’re using conservative loading and since the material is formulated w.r.t. a symmetric measure C = F.T*F, it’s ok.
But, if you’re using a follower load or you have F isolated (not in a symmetric version) in your virtual work expression, this leads to wrong results.

So, since there were a couple of other smaller mistakes, I just share the new and corrected dolfin-x example (compatible with dolfin-x from today). Feel free to test it for the two loading cases “dirichlet” or “neumann” if it behaves well. No guarantee for correctness, so would be great if you let me know if there are still mistakes (works with the above posted input files - you may have to adapt the paths for reading in…):

#!/usr/bin/env python3

# This is a minimal nonlinear quasi-static hyperelasticity example
# with a Neohookean block that is uniaxially loaded, either with a
# prescribed Dirichlet condition (if loading = 'dirichlet')
# or with a prescribed true Neumann traction (if loading = 'neumann');
# loads are ramped up with time t \in [0,1];
# the block is constrained such that a uniaxial stress state (in x-dir)
# is obtained

import time
import sys, os, subprocess, time
import math

import numpy as np
from mpi4py import MPI
from petsc4py import PETSc

from dolfinx import Function, FunctionSpace, VectorFunctionSpace, TensorFunctionSpace, DirichletBC
from dolfinx.fem import locate_dofs_topological, locate_dofs_geometrical, assemble_matrix, assemble_vector, set_bc, apply_lifting
from dolfinx.io import XDMFFile
from dolfinx.la import solve
from ufl import TrialFunction, TestFunction, FacetNormal, inner, dx, ds, Identity, grad, det, tr, dot, inv, derivative, variable, diff, ln, sqrt, cos, pi

start = time.time()

# base path for reading input from folder in and writing output to folder out
basepath  = '/home/mh/work/fem_scripts/testing'

loading = 'dirichlet' # dirichlet, neumann

T       = 1.0
nstep   = 10
dt      = T/nstep

# in which interval to write results
write_results_every = -999

# encoding and mpi communicator
encoding = XDMFFile.Encoding.ASCII
comm = MPI.COMM_WORLD

# read in xdmf mesh - domain
with XDMFFile(comm,''+basepath+'/input/block_domain.xdmf','r',encoding=encoding) as infile:
    mesh = infile.read_mesh(name="Grid")

# read in xdmf mesh - boundary
mesh.topology.create_connectivity(2, mesh.topology.dim)
with XDMFFile(comm,''+basepath+'/input/block_boundary.xdmf','r',encoding=encoding) as infile:
    mt_s = infile.read_meshtags(mesh, name="Grid_surf")

# reference normal on mesh facets (for true Neumann load, or other applications)
n0 = FacetNormal(mesh)

# function spaces for displacement, stresses, strains
# CG: Lagrange, DG: discontinuous Lagrange
V_u          = VectorFunctionSpace(mesh, ("CG", 1))
V_tensor     = TensorFunctionSpace(mesh, ("DG", 0))
V_real       = FunctionSpace(mesh, ("R", 0))

### define functions
du    = TrialFunction(V_u)            # Incremental displacement
var_u = TestFunction(V_u)             # Test function
u     = Function(V_u, name='Displacement')
# disp increment for Newton
del_u = Function(V_u)
# strain and stress functions for postprocessing
glstrain     = Function(V_tensor, name="GreenStrain")
cauchystress = Function(V_tensor, name="CauchyStress")

# print infos on problem size
if comm.rank == 0:
    print("Number of degrees of freedom: %i" % (V_u.dim))
    sys.stdout.flush()


# a zero function
zero = Function(V_u)
zero.vector.set(0.0)

# expression to interpolate time-dependent loads into
class template_expression:
    def __init__(self):
        self.val = 0.0

    def evaluate(self, x):
        return np.full(x.shape[1], self.val)

# expression for time-dependent Dirichlet condition
def dbc_uni(t):
    umax = 0.1
    return umax*t/T

# expression for time-dependent Neumann load
def neumann_true(t):
    pmax = 10.
    return 0.5*pmax*(1.-cos(pi*t/T))

loadneu = template_expression()
loadneu_func = Function(V_real)
loadneu_func.interpolate(loadneu.evaluate)

loaddbc = template_expression()
loaddbc_func = Function(V_u)
loaddbc_func.interpolate(loaddbc.evaluate)


# Boundary conditions
dbcs=[]

dbcs.append( DirichletBC(zero, locate_dofs_topological(V_u.sub(0), 2, mt_s.indices[mt_s.values == 1])) ) # x=0
dbcs.append( DirichletBC(zero, locate_dofs_topological(V_u.sub(1), 2, mt_s.indices[mt_s.values == 2])) ) # y=0
dbcs.append( DirichletBC(zero, locate_dofs_topological(V_u.sub(2), 2, mt_s.indices[mt_s.values == 3])) ) # z=0

if loading=='dirichlet':

    dbcs.append( DirichletBC(loaddbc_func, locate_dofs_topological(V_u.sub(0), 2, mt_s.indices[mt_s.values == 4])) ) # x-load Dirichlet
    
elif loading=='neumann':

    ds4 = ds(subdomain_data=mt_s, subdomain_id=4) # x-load Neumann
    
else:
    raise NameError("Set loading to either dirichlet or neumann!")


### nonlinear finite strain kinematics
I = Identity(len(u))  # Identity tensor

F = I + grad(u) # deformation gradient
J = det(F)            # determinant of deformation gradient
# use variable expression since later we diff w.r.t these
C = variable(F.T*F)   # right Cauchy-Green tensor

# Cauchy-Green invariants
Ic   = tr(C)
IIc  = 0.5*(tr(C)**2. - tr(C*C))
IIIc = det(C)
# isochoric Cauchy-Green invariants
Ic_bar   = IIIc**(-1./3.) * Ic
IIc_bar  = IIIc**(-2./3.) * IIc


### material law - Neohookean, isochoric-volumetric split
mu, nu = 10., 0.49
kappa = mu/(1.-2.*nu)
# isochoric strain-energy function
psi_iso = (mu/2.) * (Ic_bar - 3.) # isochoric strain energy function
# volumetric strain-energy function
psi_vol = (kappa/4.) * (IIIc - 2.*ln(sqrt(IIIc)) - 1.) # Ogden-type volumetric part
# total strain-energy function
psi = psi_iso + psi_vol


#2nd Piola-Kirchhoff stress 
def S():
    return 2.*diff(psi,C)

# Green-Lagrange strain (only for post-processing)
def E():
    return 0.5*(F.T*F - I)    

#Cauchy stress (only for post-processing)
def sig():
    return (1./J) * F*S()*F.T


# internal virtual work
def dW_int():
    
    return 0.5*inner(S(),derivative(C, u, var_u))*dx


# external virtual work
def dW_ext():
    
    if loading=='dirichlet':
        
        return 0.0
    
    elif loading=='neumann':
        
        return loadneu_func*J*dot(dot(inv(F).T,n0), var_u)*ds4
    
    else:
        raise NameError("Set loading to either dirichlet or neumann!")


# nonlinear variational form: internal minus external virtual work
varform = dW_int() - dW_ext()

# Jacobian of nonlinear variational form
jac = derivative(varform, u, du)


if write_results_every > -1:
    results_xdmf_file = XDMFFile(comm,''+basepath+'/out/results.xdmf','w',encoding=encoding)
    results_xdmf_file.write_mesh(mesh)


# load/time stepping
interval = np.linspace(0, T, nstep+1)

for (N, dt) in enumerate(np.diff(interval)):
    
    t = interval[N+1]

    # update time-dependent expressions
    loadneu.val = neumann_true(t)
    loadneu_func.interpolate(loadneu.evaluate)

    loaddbc.val = dbc_uni(t)
    loaddbc_func.interpolate(loaddbc.evaluate)

    # Newton loop
    struct_res_norm = 1.0
    struct_inc_norm = 1.0
    tolres = 1.0e-8
    tolinc = 1.0e-8
    it = 0
    maxiter = 25

    # create solver
    ksp = PETSc.KSP().create(comm)
    ksp.setType("preonly")
    ksp.getPC().setType("lu")
    ksp.getPC().setFactorSolverType("superlu_dist")

    while it < maxiter:
        
        # assemble rhs vector
        r = assemble_vector(varform)
        apply_lifting(r, [jac], [dbcs], x0=[u.vector], scale=-1.0)
        r.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        set_bc(r, dbcs, x0=u.vector, scale=-1.0)
        
        # assemble system matrix
        K = assemble_matrix(jac, dbcs)
        K.assemble()

        # solve linear system
        ksp.setOperators(K)
        ksp.solve(-r, del_u.vector)
        
        # get residual and increment norm
        struct_res_norm = r.norm()
        struct_inc_norm = del_u.vector.norm()
        
        # update solution
        u.vector.axpy(1.0, del_u.vector)
        u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
        
        if comm.rank == 0:
            print("%i     %.4e         %.4e" % (it, struct_res_norm, struct_inc_norm))
            sys.stdout.flush()
        
        # increase iteration index
        it += 1

        # check if converged
        if struct_res_norm <= tolres and struct_inc_norm <= tolinc:
            break
    
    else:

        if comm.rank == 0:
            print("Newton did not converge after %i iterations!" % (it))
            sys.stdout.flush()
        
        sys.exit()    


    # write results every write_results_every steps
    if (N+1) % write_results_every == 0:
        
        # Save solution to XDMF format
        results_xdmf_file.write_function(u, t)
        
    if comm.rank == 0: # only proc 0 should print this
        print("### TIME STEP %i / %i successfully completed, TIME: %.4f" % (N+1,nstep,t))
        print("-----------------------------------------------------------")
        sys.stdout.flush()
        

if comm.rank == 0: # only proc 0 should print this
    print('Time for computation: %.1f min' % ( (time.time()-start)/60. ))
    sys.stdout.flush()
1 Like
10

Thanks, Marc. As far as I can tell, the code is correct. I have examined this line carefully

and it appears correct to me. Confirmation came when I realized that it does match the corresponding syntax in the elasticity hands-on tutorial.

Thanks again!

11

Ah cool, I did not know this tutorial. So, the load is a follower pressure/tension load acting on the current configuration, so the true traction would be t = -p n (pressure if p positive, tension if negative). Pull-back then should yield t0 = -p J F^(-T) n0. In the loaded configuration, the axial Cauchy stress should equal the value of p.

Bw,
Marc

12

Having a closer look at the demo, I am not sure if the loading is of the same nature as mine.
The line J T F^(-T) n = g in the tutorial is just another way of saying the boundary traction is t0 = P n0 = g, while I assume their n is my n0 (reference normal) and the T (I prefer calling it sigma) is the Cauchy stress. (My t0 is first Piola-Kirchhoff traction and P first Piola-Kirchhoff stress.)
In my case, the value of g actually is deformation-dependent, so in tutorial way this would be J sigma F^(-T) n0 = -p J F^(-T) n0 which would just be sigma n = -p n (with sigma as Cauchy stress, n as true normal).
So guess this is different as far as I see. Which is not wrong of course, just a different kind of load. The virtual work expression for a “standard” traction in my terminology would be
dot(loadneu_func_vec, var_u)*ds4 where now loadneu_func_vec can be some prescribed vector function.

Let me know if this makes sense or maybe if I got something wrong.

Bw,
Marc

1 Like
13

Hi !
I have a question. You expressed the second Piola-Kirchhoff stress in your code. Similarly, can I express the first Piola-Kirchhoff stress as

FF = variable(F)
P = diff(psi,FF)

Thank you very much!

14

Hi, sure, this should yield equivalent results!
You only need to change the internal virtual work expression accordingly, from 0.5*inner(S(),derivative(C, u, var_u))*dx to inner(P(),derivative(F, u, var_u))*dx. And guess you can also express derivative(F, u, var_u) with grad(var_u).

Bw,
Marc

15

Hi dokken, is there any chance for the local solver to be added in dolfinx? I’m solving a plasticity problem where the global projection to internal variables (defined on quad points) is really slow. The local projection feature is therefore highly desired.

16

You would have to create an issue on the DOLFINx issue tracker. As I’ve never used the LocalSolver, im not aware of all the details of what it does, and what would be required to recreate it.