Minimization over functions of different shapes

variational formulation
1

Hello,
I successfully minimized the energy of a strained sheet for the displacement field (2d vector). However, when I add the tensorial field q_{\alpha\beta} to the energy, A^{\alpha\beta\gamma\delta}q_{\alpha\beta}q_{\gamma\delta}, the resulting energy becomes 300 orders of magnitude higher, 2.3046e+302 instead of 0.02325. Of course, this doesn’t make sense because the tensorial field could have vanished, so the energy would not have changed. What did I do wrong?

Here is a minimal working example:

import ufl
from ufl import *
import dolfinx
from dolfinx import *
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
import numpy as np
domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

#The trial functions are the 2d vector u and the degree 2 tensor q
VFS = ufl.VectorElement('CG',domain.ufl_cell(),1, dim=2)
TFS = ufl.TensorElement('CG', domain.ufl_cell(),degree=2)
mel = MixedElement([VFS,TFS])
V = dolfinx.fem.FunctionSpace(domain, VFS)
Q = dolfinx.fem.FunctionSpace(domain, TFS)
Y = dolfinx.fem.FunctionSpace(domain, mel)
y = fem.Function(Y)
u, q = ufl.split(y)
yt = ufl.TestFunction(Y)
v,p = ufl.split(yt)

#impose clamped boundary conditions
def lower_boundary(x):
    return np.isclose(x[1],0);
def upper_boundary(x):
    return np.isclose(x[1],1);

fdim = domain.topology.dim - 1
boundary_facets_up = mesh.locate_entities_boundary(domain, fdim, upper_boundary)
dofs_up = fem.locate_dofs_topological(V, fdim, boundary_facets_up)
u_up = ScalarType((0,0.1))
bc_up = fem.dirichletbc(u_up, dofs_up, V)

boundary_facets_down = mesh.locate_entities_boundary(domain, fdim, lower_boundary)
dofs_down = fem.locate_dofs_topological(V, fdim, boundary_facets_down)
u_down = ScalarType((0,-0.1))
bc_down = fem.dirichletbc(u_down, dofs_down, V)

bc = [bc_up, bc_down]


#The energy
mu = 0.4
Lambda = 0.5
tdim=domain.topology.dim
# define strain and stress

def epsilon(u):
    return sym(grad(u))
def sigma(u):
    return 2*mu*epsilon(u) + Lambda*tr(epsilon(u))*Identity(tdim)
def sigmaq(q):
    return 2*mu*q + Lambda*tr(q)*Identity(q.ufl_shape[0])
Energy = 1/2*inner(sigma(u),epsilon(u))*dx+1/8*inner(sigmaq(q),q)*dx


# Minimization
Res = derivative(Energy, y, yt)
problem = fem.petsc.NonlinearProblem(Res, y, bcs=bc)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.report = True
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(y)
assert(converged)
print(f"Number of interations: {n:d}")
print("Energy = ", fem.assemble_scalar(fem.form(Energy)))
2

Applying Dirichlet bcs to mixed spaces requires modified syntax.

This says that the bc should be applied to the function space V. However, in the mixed formulation you are solving for a function y in the mixed space Y, so these bcs have no effect; i.e. you are solving an ill-posed problem.

Modifying the bcs as follows yields the original energy:

Number of interations: 3
Energy =  0.023253452428557055

boundary_facets_up = mesh.locate_entities_boundary(domain, fdim, upper_boundary)
dofs_up = fem.locate_dofs_topological((Y.sub(0), V), fdim, boundary_facets_up)
u_up = fem.Function(V)
u_up.interpolate(lambda x: np.broadcast_to(ScalarType([[0],[0.1]]), (2, x.shape[1])))
bc_up = fem.dirichletbc(u_up, dofs_up, Y)

boundary_facets_down = mesh.locate_entities_boundary(domain, fdim, lower_boundary)
dofs_down = fem.locate_dofs_topological((Y.sub(0), V), fdim, boundary_facets_down)
u_down = fem.Function(V)
u_down.interpolate(lambda x: np.broadcast_to(ScalarType([[0],[-0.1]]), (2, x.shape[1])))
bc_down = fem.dirichletbc(u_down, dofs_down, Y)
2 Likes