Variational formulation of elasticity problem when using different function types

variational formulation
1

Hi there,

I am trying to define regions with different material properties, on which to solve an elasticity problem with displacement defined as a CG2 function.
When mu and lambda are also set using CG2 functions, the model runs fine, but I feel material properties should be defined element-wise. (MWE1)

MWE1:

import numpy as np
import ufl

from mpi4py import MPI
from petsc4py.PETSc import ScalarType

from dolfinx import mesh, fem, plot, io, nls,log
from dolfinx.fem import FunctionSpace,VectorFunctionSpace, Function,TensorFunctionSpace, Constant
from ufl import TestFunction, TrialFunction, dx, inner,grad, Measure,derivative,system
from petsc4py import PETSc
import time

startT = time.time()

L = 15.0
l= 10.0
d= 5.0
domain = mesh.create_box(MPI.COMM_WORLD,[[0.0,0.0,0.0], [L, l, d]], [15,20,5], mesh.CellType.hexahedron)
V = fem.VectorFunctionSpace(domain, ("CG", 2))
Vdim1=fem.FunctionSpace(domain, ("CG", 2))

xdmf = io.XDMFFile(domain.comm, "YoungsModulus.xdmf", "w")
xdmf.write_mesh(domain)



def bottom(x):
    return np.isclose(x[1], 0)

def top(x):
    return np.isclose(x[1], l)


fdim = domain.topology.dim -1
bottom_facets = mesh.locate_entities_boundary(domain, fdim, bottom)
top_facets = mesh.locate_entities_boundary(domain, fdim, top)
marked_facets = np.hstack([bottom_facets, top_facets])
marked_values = np.hstack([ np.full_like(bottom_facets, 1), np.full_like(top_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

u_bc = np.array((0,) * domain.geometry.dim, dtype=PETSc.ScalarType)
bottom_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(1))
bc1=fem.dirichletbc(u_bc, bottom_dofs, V)
loaded_boundary_dofs_z = fem.locate_dofs_topological(V.sub(1), facet_tag.dim, facet_tag.find(2)) 
u_AD = ScalarType(-0.1)
bc2 = fem.dirichletbc(u_AD, loaded_boundary_dofs_z, V.sub(1)) 

bcs = [bc1,bc2]


B = fem.Constant(domain, PETSc.ScalarType((0, 0, 0)))
T = fem.Constant(domain, PETSc.ScalarType((0, 0, 0)))

v = ufl.TestFunction(V)
u = fem.Function(V)

d = len(u)

I = ufl.variable(ufl.Identity(d))
F = ufl.variable(I + ufl.grad(u))
C = ufl.variable(F.T * F)
Ic = ufl.variable(ufl.tr(C))
J  = ufl.variable(ufl.det(F))

E_gel = fem.Constant(domain,PETSc.ScalarType(150))
pr_gel = fem.Constant(domain,PETSc.ScalarType(0.3))
E_miner = fem.Constant(domain,PETSc.ScalarType(15000))

E_field=Function(Vdim1)
mu_field=Function(Vdim1)
lmbda_field=Function(Vdim1)

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

coord=Vdim1.tabulate_dof_coordinates()
Evalues=E_field.x.array 
for i, (x,y,z) in enumerate(coord):
    if y<=5.2:
        Evalues[i]=E_miner.value
    else:
        Evalues[i]=E_gel.value
E_field.x.array[:]=Evalues

lmbda_field.x.array[:] = E_field.x.array*pr_gel.value/((1+pr_gel.value)*(1-2*pr_gel.value))
mu_field.x.array[:] = E_field.x.array/(2*(1+pr_gel.value))  

# Stored strain energy density (compressible neo-Hookean model)
psi = (mu_field / 2) * (Ic - 3) - mu_field * ufl.ln(J) + (lmbda_field / 2) * (ufl.ln(J))**2

P = ufl.diff(psi, F)

metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
F = ufl.inner(ufl.grad(v), P)*dx - ufl.inner(v, B)*dx - ufl.inner(v, T)*ds(2) 

problem = fem.petsc.NonlinearProblem(F, u, bcs)
solver = nls.petsc.NewtonSolver(domain.comm, problem)

# Set Newton solver options
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"


# Solve problem in time steps
log.set_log_level(log.LogLevel.INFO)
tval0 = 0
T.value[2] = tval0
num_its, converged = solver.solve(u)
assert(converged)
u.x.scatter_forward()
xdmf.write_function(E_field, 1)          
xdmf.close()


endT = time.time()
print("The time of execution of above program is :",(endT-startT),flush=True)

When mu and lambda are defined using DG0 functions, and displacement with (Lagrange,1) the model also runs but plotting material properties shows a layer of elements with arbitrary low Youngs modulus (MWE2).

image
image1141×651 49.6 KB

MWE2:

import numpy as np
import ufl

from mpi4py import MPI
from petsc4py.PETSc import ScalarType

from dolfinx import mesh, fem, plot, io, nls,log
from dolfinx.fem import FunctionSpace,VectorFunctionSpace, Function,TensorFunctionSpace, Constant
from ufl import TestFunction, TrialFunction, dx, inner,grad, Measure,derivative,system
from petsc4py import PETSc
import time

startT = time.time()

L = 15.0
l= 10.0
d= 5.0

domain = mesh.create_box(MPI.COMM_WORLD,[[0.0,0.0,0.0], [L, l, d]], [15,20,5], mesh.CellType.hexahedron)

E_gel = fem.Constant(domain,PETSc.ScalarType(150))
pr_gel = fem.Constant(domain,PETSc.ScalarType(0.3))
E_miner = fem.Constant(domain,PETSc.ScalarType(15000))

"""V = fem.VectorFunctionSpace(domain, ("CG", 2))
Vdim1=fem.FunctionSpace(domain, ("CG", 2))"""
V = fem.VectorFunctionSpace(domain, ("Lagrange", 1))
Vdim1=fem.FunctionSpace(domain, ("Lagrange", 1))
Q = FunctionSpace(domain, ("DG", 0))


xdmf = io.XDMFFile(domain.comm, "YoungsModulusDG.xdmf", "w")
xdmf.write_mesh(domain)


def Omega_0(x):
    return x[1] <= 5.2

def Omega_1(x):
    return x[1] > 5.2

Young=Function(Q)
mu_field=Function(Q)
lmbda_field=Function(Q)


cells_0 = mesh.locate_entities(domain, domain.topology.dim, Omega_0)
cells_1 = mesh.locate_entities(domain, domain.topology.dim, Omega_1)

Young.x.array[cells_0] = np.full_like(cells_0, E_miner.value, dtype=PETSc.ScalarType)
Young.x.array[cells_1] = np.full_like(cells_1, E_gel.value, dtype=PETSc.ScalarType)

xdmf.write_function(Young, 1)          
xdmf.close()


def bottom(x):
    return np.isclose(x[1], 0)

def top(x):
    return np.isclose(x[1], l)


fdim = domain.topology.dim -1
bottom_facets = mesh.locate_entities_boundary(domain, fdim, bottom)
top_facets = mesh.locate_entities_boundary(domain, fdim, top)
marked_facets = np.hstack([bottom_facets, top_facets])
marked_values = np.hstack([ np.full_like(bottom_facets, 1), np.full_like(top_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

u_bc = np.array((0,) * domain.geometry.dim, dtype=PETSc.ScalarType)
bottom_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(1))
bc1=fem.dirichletbc(u_bc, bottom_dofs, V)
loaded_boundary_dofs_z = fem.locate_dofs_topological(V.sub(1), facet_tag.dim, facet_tag.find(2)) 
u_AD = ScalarType(-0.1)
bc2 = fem.dirichletbc(u_AD, loaded_boundary_dofs_z, V.sub(1)) 

bcs = [bc1,bc2]


B = fem.Constant(domain, PETSc.ScalarType((0, 0, 0)))
T = fem.Constant(domain, PETSc.ScalarType((0, 0, 0)))

v = ufl.TestFunction(V)
u = fem.Function(V)

d = len(u)

I = ufl.variable(ufl.Identity(d))
F = ufl.variable(I + ufl.grad(u))
C = ufl.variable(F.T * F)
Ic = ufl.variable(ufl.tr(C))
J  = ufl.variable(ufl.det(F))



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

coord=Vdim1.tabulate_dof_coordinates()


lmbda_field.x.array[:] = Young.x.array*pr_gel.value/((1+pr_gel.value)*(1-2*pr_gel.value))
mu_field.x.array[:] = Young.x.array/(2*(1+pr_gel.value))  

# Stored strain energy density (compressible neo-Hookean model)
psi = (mu_field / 2) * (Ic - 3) - mu_field * ufl.ln(J) + (lmbda_field / 2) * (ufl.ln(J))**2

P = ufl.diff(psi, F)

metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
F = ufl.inner(ufl.grad(v), P)*dx - ufl.inner(v, B)*dx - ufl.inner(v, T)*ds(2) 

problem = fem.petsc.NonlinearProblem(F, u, bcs)
solver = nls.petsc.NewtonSolver(domain.comm, problem)

# Set Newton solver options
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"


# Solve problem in time steps
log.set_log_level(log.LogLevel.INFO)
tval0 = 0
T.value[2] = tval0
num_its, converged = solver.solve(u)
assert(converged)
u.x.scatter_forward()



endT = time.time()
pritype or paste code here

Finally, when mu and lambda are defined using DG0 functions, and displacement with CG2, the solver throws error: RuntimeError: Failed to successfully call PETSc function ‘KSPSolve’. PETSc error code is: 76,

I realise I am most likely facing dimension inconsistencies but I am really confused as to how to best address this problem.

Thank you in advance for your help.

2

This only gives you all cells where every vertex satisfies the constraint.
This means that in the transfer region where some vertices are in omega 0 and others in omega 1 you have cells that fall outside each category.
See A question about Function numerical setting - #3 by hherlyng

2 Likes
3

Oh I see, thank you very much, that solved my issue.