BCs with MixedFunctionSpace

dolfinx
1

Hi,

while changing my code from mixed_element to MixedFunctionSpace (see the previous post here), I ran into convergence problems (solver produces NaN), which I think stem from the bc set up.

This was before:

import numpy as np
from mpi4py import MPI
from dolfinx import fem, mesh
import basix
from dolfinx.fem import locate_dofs_geometrical
from basix.ufl import mixed_element

p = 1
domain = mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [1, 1, 1]], [1, 1, 1], mesh.CellType.hexahedron)

Ue = basix.ufl.element("P", domain.basix_cell(), p, shape=(domain.geometry.dim,))
He = basix.ufl.quadrature_element(domain.basix_cell(), value_shape=(), degree=p)
Te = basix.ufl.blocked_element(He, shape=(domain.geometry.dim, domain.geometry.dim), symmetry=True)

V_el = mixed_element([Ue, Te, He, Te])
V = fem.functionspace(domain, V_el)

Uf = fem.functionspace(domain, Ue)
Tf = fem.functionspace(domain, Te)
Hf = fem.functionspace(domain, He)

V0, _ = V.sub(0).collapse()
u_L = fem.Function(V0)
u_L.x.array[:] = 0.

dofs_L = locate_dofs_geometrical((V.sub(0), V0), lambda x: np.isclose(x[0], 0))
bc_L = fem.dirichletbc(u_L, dofs_L, V.sub(0))

print(dofs_L)

And this is my attempt with MixedFunctionSpace, which results in a different dof vector:

import numpy as np
import ufl
from mpi4py import MPI
from dolfinx import fem, mesh
import basix
from dolfinx import  default_scalar_type

p = 1
domain = mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [1, 1, 1]], [1, 1, 1], mesh.CellType.hexahedron)

Ue = basix.ufl.element("P", domain.basix_cell(), p, shape=(domain.geometry.dim,))
He = basix.ufl.quadrature_element(domain.basix_cell(), value_shape=(), degree=p)
Te = basix.ufl.blocked_element(He, shape=(domain.geometry.dim, domain.geometry.dim), symmetry=True)

Uf = fem.functionspace(domain, Ue)
Tf = fem.functionspace(domain, Te)
Hf = fem.functionspace(domain, He)

V = ufl.MixedFunctionSpace(Uf, Tf, Hf, Tf)

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

dofs_L = fem.locate_dofs_geometrical(Uf, lambda x: np.isclose(x[0], 0))
u_bc = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)

bc_L = fem.dirichletbc(u_bc, dofs_L, Uf)

print(dofs_L)

I am using Fenicx 0.9 on WSL Ubuntu / Windows. Many thanks!

2

You would expect the dof vector to be different; the one references dofs inside the mixed space and the other references dofs of the ‘loose’ functionspace.

At first glance your code looks correct. Also take a look at Multiphysics: Solving PDEs on subdomains — FEniCS Workshop to confirm. Although there locate_dofs_topological is used.

What makes you think it is the BCs?

I’d advise to build up complexity more slowly. First, you could solve a linear elasticity problem on your domain with the MixedFunctionSpace and the BC enforcement per the above. Then solve that linear problem with the nonlinear block solver. And only then move to the nonlinear case that introduces plasticity. That’ll help pinpoint where exactly the issues are introduced.

1 Like
3

Thank you for the quick response! I was able to solve the linear problem with MixedFunctionSpace, but not the nonlinear problem. When reformulating it for the nonlinear block solver, I changed u, = ufl.TrialFunctions(V) to u = fem.Function(Uf) . Is that correct?

MWE:

from dolfinx import log, default_scalar_type
import numpy as np
import ufl
from mpi4py import MPI
from dolfinx import fem, mesh
import basix
from scifem import BlockedNewtonSolver

domain = mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [1, 1, 1]], [1, 1, 1], mesh.CellType.hexahedron)

p = 1
Ue = basix.ufl.element("P", domain.basix_cell(), p, shape=(domain.geometry.dim,))
He = basix.ufl.quadrature_element(domain.basix_cell(), value_shape=(), degree=p)
Te = basix.ufl.blocked_element(He, shape=(domain.geometry.dim, domain.geometry.dim), symmetry=True)

Uf = fem.functionspace(domain, Ue)
Tf = fem.functionspace(domain, Te)
Hf = fem.functionspace(domain, He)

V = ufl.MixedFunctionSpace(Uf)

u = fem.Function(Uf)
du, = ufl.TestFunctions(V)

mu = fem.Constant(domain, 100.)
la = fem.Constant(domain, 10.)

d = len(u)
Id = ufl.variable(ufl.Identity(d))
P = 2.0 * mu * ufl.sym(ufl.grad(u)) + la * ufl.tr(ufl.sym(ufl.grad(u))) * Id

dx = ufl.Measure("dx", domain=domain)
form = ufl.inner(ufl.grad(du), P) * dx 

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

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

fdim = domain.topology.dim - 1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_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=default_scalar_type)
left_dofs = fem.locate_dofs_topological(Uf, facet_tag.dim, facet_tag.find(1))
bcs = [fem.dirichletbc(u_bc, left_dofs, Uf)]

F_compiled = fem.form(ufl.extract_blocks(form))
J_compiled = fem.form(ufl.extract_blocks(ufl.derivative(form, u, ufl.TrialFunctions(V))))
petsc_options = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"}
log.set_log_level(log.LogLevel.INFO)
blocked_solver = BlockedNewtonSolver(F_compiled, [u], J=J_compiled)
blocked_solver.convergence_criterion = "incremental"
blocked_solver.solve()


What makes you think it is the BCs?

I think I got confused about the meaning of the indices. I expected that locate_dofs_geometrical would give vector with 4x3 entries for the facet in 3D, but now I conclude from it that the indices are vectorial.

4

Strangely, the nonlinear block solver converges for the example above, when I increase the number of elements to

domain = mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [1, 1, 1]], [2, 1, 1], mesh.CellType.hexahedron)

But it diverges with the smallest forcing, which I applied in the following way:

from dolfinx import log, default_scalar_type
import numpy as np
import ufl
from mpi4py import MPI
from dolfinx import fem, mesh
import basix
from scifem import BlockedNewtonSolver

domain = mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [1, 1, 1]], [2, 1, 1], mesh.CellType.hexahedron)
p = 1
Ue = basix.ufl.element("P", domain.basix_cell(), p, shape=(domain.geometry.dim,))
He = basix.ufl.quadrature_element(domain.basix_cell(), value_shape=(), degree=p)
Te = basix.ufl.blocked_element(He, shape=(domain.geometry.dim, domain.geometry.dim), symmetry=True)

Uf = fem.functionspace(domain, Ue)
Tf = fem.functionspace(domain, Te)
Hf = fem.functionspace(domain, He)

V = ufl.MixedFunctionSpace(Uf)

u = fem.Function(Uf)
du, = ufl.TestFunctions(V)

mu = fem.Constant(domain, 100.)
la = fem.Constant(domain, 10.)

d = len(u)
Id = ufl.variable(ufl.Identity(d))
P = 2.0 * mu * ufl.sym(ufl.grad(u)) + la * ufl.tr(ufl.sym(ufl.grad(u))) * Id

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

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

fdim = domain.topology.dim - 1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_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=default_scalar_type)
left_dofs = fem.locate_dofs_topological(Uf, facet_tag.dim, facet_tag.find(1))
bcs = [fem.dirichletbc(u_bc, left_dofs, Uf)]

factor = 0.00001
T = fem.Constant(domain, default_scalar_type((factor, factor, factor)))
metadata = {"quadrature_degree": 4}
dx = ufl.Measure("dx", domain=domain)
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
form = ufl.inner(ufl.grad(du), P) * dx  - ufl.inner(du, T) * ds(2)

F_compiled = fem.form(ufl.extract_blocks(form))
J_compiled = fem.form(ufl.extract_blocks(ufl.derivative(form, u, ufl.TrialFunctions(V))))

petsc_options = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"}
log.set_log_level(log.LogLevel.INFO)
blocked_solver = BlockedNewtonSolver(F_compiled, [u], J=J_compiled)
blocked_solver.convergence_criterion = "incremental"
blocked_solver.solve()