A doubt about imposing bc on mixed element

General
,
1

Hi all! I have a question about imposing the dirichletbc on the boundaries.
I consider a mixed element on a rectangular domain with length L and height d. I clamp the left end and impose a Delta at the right end. Please see below

import dolfinx
from dolfinx import nls
import ufl 
import numpy as np
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import dolfinx.nls.petsc
from dolfinx import geometry
import matplotlib.pyplot as plt
from dolfinx import fem, io, mesh
import math
from ufl import (Measure, SpatialCoordinate, TestFunctions, TrialFunctions,
                 div, exp, inner)

# Define domain parameters
L = 0.166  # total length
delta = 2.0  # aspect ratio
d=L *delta  # thickness
Delta = 0.01  # displacement at the right end
f0 =3.1415926-0.2

# Create domain
my_domain = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0.0, 0.0]), np.array([L, d])], [10,20], cell_type=dolfinx.mesh.CellType.quadrilateral)

# Define mixed mesh
u_el = ufl.VectorElement('Lagrange', my_domain.ufl_cell(), degree=1, dim=2)
phi_el = ufl.FiniteElement('Lagrange', my_domain.ufl_cell(), degree=1)
mixed_el = ufl.MixedElement([u_el, phi_el]) # mixed element
V = dolfinx.fem.FunctionSpace(my_domain, mixed_el) 
#Define the trial and testfunctions on the space V
u_gen = fem.Function(V)
u, f = ufl.split(u_gen)
u_trial_gen = ufl.TrialFunction(V)
u_trial, f_trial = ufl.split(u_trial_gen)
u_test_gen = ufl.TestFunction(V)
u_test, f_test = ufl.split(u_test_gen)

I = ufl.Identity(my_domain.topology.dim)
F = I + ufl.grad(u)
C = F.T * F
E = ufl.variable(1/2*(C-I)) # Green-Lagrange strain
J = np.array([[[1.0, 0.0], [1.0, 0.0]],
              [[0.0, 1.0], [0.0, 1.0]]])
Edag=-ufl.sin(f)**2*I /2

i,j,k = ufl.indices(3)
J_ufl = ufl.as_tensor(J)
E_deviation=ufl.as_tensor(J_ufl[i,j,k]*E[j,k],i) 

#potensial energy
psi=ufl.tr(E-Edag) ** 2+E_deviation**2
dx = Measure("dx", my_domain)
potential_energy = psi * dx

#Initialization:
u_gen.sub(1).x.array[:] =f0


# Define variational problem
residual = ufl.derivative(potential_energy,u_gen,u_test_gen)

#find the facets for the left and right boundaries
fdim = my_domain.topology.dim - 1
facets_left = mesh.locate_entities_boundary(my_domain, fdim, lambda x: np.isclose(x[0], 0.0))
facets_right = mesh.locate_entities_boundary(my_domain, fdim, lambda x: np.isclose(x[0],L))
Q, _ = V.sub(0).collapse()

#find the dofs
dofs_left_u = fem.locate_dofs_topological((V.sub(0), Q), fdim, facets_left)
dofs_left_f = fem.locate_dofs_topological(V.sub(1), fdim, facets_left)
dofs_right_u = fem.locate_dofs_topological((V.sub(0), Q), fdim, facets_right)
dofs_right_f = fem.locate_dofs_topological(V.sub(1), fdim, facets_right)

#boudary conditions
def left(x):
    vals = np.zeros((2, x.shape[1]))
    return vals
    
def right(x):
    vals= np.zeros((2, x.shape[1]))
    vals[0, :] = Delta
    return vals

#set the left bcs
uDL = fem.Function(Q)
uDL.interpolate(left)
bc_u_left = fem.dirichletbc(uDL, dofs_left_u, V.sub(0))
bc_f_left =fem.dirichletbc(ScalarType(f0), dofs_left_f, V.sub(1))

#set the right bcs
uDR = fem.Function(Q)
uDR.interpolate(right)
bc_u_right = fem.dirichletbc(uDR, dofs_right_u, V.sub(0))
bc_f_right =fem.dirichletbc(ScalarType(f0), dofs_right_f, V.sub(1))

#combine the bcs
bc = [bc_u_left, bc_f_left,bc_u_right,bc_f_right]


# Create nonlinear problem
problem = dolfinx.fem.petsc.NonlinearProblem(residual,u_gen,  bcs=bc)
solver = nls.petsc.NewtonSolver(my_domain.comm, problem)
# Set Newton solver options
solver.atol = 1e-12
solver.rtol = 1e-12
solver.convergence_criterion = "incremental"
from dolfinx import log
log.set_log_level(log.LogLevel.INFO)
niterations, reason = solver.solve(u_gen)
print('number of iterations in Newton solver:',niterations)

u_h, f_h=u_gen.split()
with dolfinx.io.XDMFFile(my_domain.comm, "out_test/solution.xdmf", "w") as file:
    file.write_mesh(my_domain)
    file.write_function(u_h)

However, the results in Paraview shows that the boundary conditions at the right end is not fully imposed

image
image436×756 74 KB

Could you please inform me where the errors come from?
Many thanks for your help!

2

First try to change this to:

with dolfinx.io.XDMFFile(my_domain.comm, "out_test/solution.xdmf", "w") as file:
    file.write_mesh(my_domain)
    file.write_function(u_h.collapse())

Do you get the same result?

3

Thanks dokken!

Now the boundary conditions are correct. Could you explain why we need …collapse()?

4

When just using Function.split, the underlying array containing degrees of freedom is not split into a continuous chunk of data, which is required by write_function.

1 Like