Please help me to check whether boundary condition are right set

Fenicsx is excellent, but I have been tormented by the boundary conditions seting.

I have two equations
- \phi '' -4 A' \phi ' -8 \phi^3 +6 \phi^5 =0
- A '' - 2/3 \phi'^2 =0
the coordinate x \in [0,10] , with the boundary condations
A(0) =0, A'(0) = 0, \phi(0) =0, \phi'(10) =0

The bcs had been set by the following codes but I’m not sure if they are right.
Where \phi \to u1 , \ A \to u2 .

r_0 = 0.0
r_b = 10.0 
 # the coordinate of right boundary
msh = mesh.create_interval(comm=MPI.COMM_WORLD, points=(r_0, r_b), nx=100)
fdim = msh.topology.dim - 1
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
V = FunctionSpace(msh, ufl.MixedElement(P1, P1))

# the left boundary
facets_l = mesh.locate_entities_boundary(msh, dim=fdim, marker=lambda x: np.isclose(x[0], r_0))

dofs_u1 = fem.locate_dofs_topological(V=V.sub(0), entity_dim=fdim, entities=facets_l)
dofs_u2 = fem.locate_dofs_topological(V=V.sub(1), entity_dim=fdim, entities=facets_l)

# the Dirichleb boundary conditions
bc_u1 = fem.dirichletbc(value=ScalarType(0), dofs=dofs_u1, V=V.sub(0))  # u1 = phi
bc_u2 = fem.dirichletbc(value=ScalarType(0), dofs=dofs_u2, V=V.sub(1))  # u2 = A

bcs = [bc_u1, bc_u2]

The values of are dofs_u1 = [0] , dofs_u2 = [2] that confued me.

And the whole code are

import numpy as np

import ufl
from dolfinx import mesh, fem
from dolfinx.fem import Function, FunctionSpace
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, split, TestFunctions

from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType

# parmeters
r_0 = 0.0
r_b = 10.0  # the coordinate of right boundary


msh = mesh.create_interval(comm=MPI.COMM_WORLD, points=(r_0, r_b), nx=100)

fdim = msh.topology.dim - 1
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
V = FunctionSpace(msh, ufl.MixedElement(P1, P1))

# the left boundary
facets_l = mesh.locate_entities_boundary(
    msh, dim=fdim, marker=lambda x: np.isclose(x[0], r_0))

dofs_u1 = fem.locate_dofs_topological(V=V.sub(0), entity_dim=fdim, entities=facets_l)
dofs_u2 = fem.locate_dofs_topological(V=V.sub(1), entity_dim=fdim, entities=facets_l)

# the Dirichleb boundary conditions
bc_u1 = fem.dirichletbc(value=ScalarType(0), dofs=dofs_u1, V=V.sub(0))  # u1 = phi
bc_u2 = fem.dirichletbc(value=ScalarType(0), dofs=dofs_u2, V=V.sub(1))  # u2 = A

bcs = [bc_u1, bc_u2]

# Define functions
u = Function(V)
# Split system functions to access components
#phi, m, sigma = u1, u2, u3
u1, u2 = split(u)

# Define test functions
v_1, v_2 = TestFunctions(V)

# Define variational problem

F1 = v_1.dx(0)*u1.dx(0) + (-4.0 * u2.dx(0)*u1.dx(0) - 8.0*u1**3+6.0*u1**5)*v_1
F2 = v_2.dx(0)*u2.dx(0) - 2.0/3.0*u1.dx(0)**2*v_2
F = (F1+F2)*dx

#  Solve the system.
#
# t1 = data[:,1]
# t2 = data[:,2]
# t3 = data[:,3]
u.x.array[:] = np.random.random(202)

# Create nonlinear problem and Newton solver
problem = NonlinearProblem(F, u, bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-5


# We can customize the linear solver used inside the NewtonSolver by
# modifying the PETSc options
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
# -

r = solver.solve(u)

The above codes can’t reach the a convegence result.

Hi

Regarding the boundary condition assignments of a 2-field problem I suggest you have a look at https://fenicsproject.discourse.group/t/dirichletbcs-assignment-for-coupled-vector-field-problem/9050.

Hi

Despite my above suggestion, I actually have the same problem as @zhaozhh that the problem does not converge.

I have simply copied the hyperelastic demo problem https://jorgensd.github.io/dolfinx-tutorial/chapter2/hyperelasticity.html and simply expanded the function space making it a two-field problem. For testing, I fixed ALL degrees of freedom of the second unknown field, wf, set the work-conjugate director stress, df, to a constant zero vector and effectively solve only for the unknown displacement field, u (wf is completely prescribed). In this sense, the two-field problem boils down to the original demo problem and should converge without any issues. However, it does not converge.

from contextlib import ExitStack

import numpy as np

import ufl
from dolfinx import fem, nls, mesh, plot, log
from mpi4py import MPI
from petsc4py import PETSc

import pdb


# ----------------------------------------------------------------------------------------------------------------------
# Create a box mesh of a beam
L = 20.0
H = 1.0
msh = mesh.create_box(MPI.COMM_WORLD, [np.array([0.0, 0.0, 0.0]), np.array([L, H, H])], [10, 1, 1],
                      mesh.CellType.hexahedron)
P1 = ufl.VectorElement("CG", msh.ufl_cell(), 1, dim=3)  # Lagrange family, 1st order
P2 = ufl.VectorElement("CG", msh.ufl_cell(), 1, dim=3)  # Lagrange family, 1st order
VY = fem.FunctionSpace(msh, ufl.MixedElement([P1, P2]))


# ----------------------------------------------------------------------------------------------------------------------
# Boundary conditions and body loads
def left(x):
    return np.isclose(x[0], 0)

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

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

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


def fixed_displacement_expression(x):
    return np.stack((np.zeros(x.shape[1]), np.zeros(x.shape[1]), np.zeros(x.shape[1])))

def incremented_displacement_expression(x):
    return np.stack((np.full(x.shape[1], 1.0e-02), np.zeros(x.shape[1]), np.zeros(x.shape[1])))

def fixed_director_expression(x):
    return np.stack((np.zeros(x.shape[1]), np.zeros(x.shape[1]), np.zeros(x.shape[1])))

# ----------------------------------------------
# Dirichlet boundary application to affected displacement dofs

# Fix displacement DOF on left face and apply axial displacement on right face
V0, submap = VY.sub(0).collapse()
fixed_displacement = fem.Function(V0)
fixed_displacement.interpolate(fixed_displacement_expression)
incremented_displacement = fem.Function(V0)
incremented_displacement.interpolate(incremented_displacement_expression)
left_u_dofs = fem.locate_dofs_geometrical((VY.sub(0),V0), left)
right_u_dofs = fem.locate_dofs_geometrical((VY.sub(0),V0), right)
bc1 = fem.dirichletbc(fixed_displacement, left_u_dofs, V0)
bc2 = fem.dirichletbc(incremented_displacement, right_u_dofs, V0)

# Fix ALL change of director degrees of freedom in the domain
V1, submap = VY.sub(1).collapse()
fixed_director = fem.Function(V1)
fixed_director.interpolate(fixed_director_expression)
top_wf_dofs = fem.locate_dofs_geometrical((VY.sub(1),V1), top)
bottom_wf_dofs = fem.locate_dofs_geometrical((VY.sub(1),V1), bottom)
bc3 = fem.dirichletbc(fixed_director, top_wf_dofs, V1)
bc4 = fem.dirichletbc(fixed_director, bottom_wf_dofs,V1)

bcs = [bc1,  bc2, bc3, bc4]
pdb.set_trace()

# ----------------------------------------------------------------------------------------------------------------------
# Function space
v, yf = ufl.TestFunctions(VY)     # test functions
uwf = fem.Function(VY)            # current displacement and change of director
u, wf = ufl.split(uwf)

# ----------------------------------------------------------------------------------------------------------------------
# Kinematics

# Spatial dimension
d = len(u)

# Identity tensor
I = ufl.variable(ufl.Identity(d))

# Deformation gradient
F = ufl.variable(I + ufl.grad(u))

# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)

# Invariants
Ic = ufl.variable(ufl.tr(C))
J  = ufl.variable(ufl.det(F))

# ----------------------------------------------------------------------------------------------------------------------
# Constitutive law

# Set the elasticity parameters
E = PETSc.ScalarType(1.0e4)
nu = PETSc.ScalarType(0.3)
mu = fem.Constant(msh, E/(2*(1 + nu)))
lmbda = fem.Constant(msh, E*nu/((1 + nu)*(1 - 2*nu)))

# Stored strain energy density (compressible neo-Hookean model)
psi_m = (mu / 2.0) * (Ic - 3.0) - mu * ufl.ln(J) + (lmbda / 2.0) * (ufl.ln(J))**2

# 1st Piola Kirchhoff stress tensor
P = ufl.diff(psi_m, F)                                     # 1st Piola Kirchhoff stress tensor
df = fem.Constant(msh, PETSc.ScalarType((0.0, 0.0, 0.0)))  # dummy director stress

# ----------------------------------------------------------------------------------------------------------------------
# Variational formulation
metadata = {"quadrature_degree": 3}
dx = ufl.Measure("dx", domain=msh, metadata=metadata)

# Define form Pi (we want to find u such that Pi(u) = 0)
Pi = ufl.inner(P, ufl.grad(v))*dx + ufl.inner(df, yf)*dx

# Initialize nonlinear problem
problem = fem.petsc.NonlinearProblem(Pi, uwf, bcs)

# ----------------------------------------------------------------------------------------------------------------------
# Initialise Newton Rhapson solver
solver = nls.petsc.NewtonSolver(msh.comm, problem)

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

# ----------------------------------------------------------------------------------------------------------------------
# Incremental application of the traction loading
log.set_log_level(log.LogLevel.INFO)
for n in range(1, 2):
    num_its, converged = solver.solve(uwf)
    assert(converged)
1 Like

Thank you for your response, @skatulla. I am a very new user of fenicsx, so I can not understand your example. In my opinion, my above problem is quit simple but I also can not to confirm the boundary conditions is set right or no . Unfortunately, I can’t find out a simliar example in the documents about fenicsx or fecics.

After I changed the initial guess solution the codes converged, but gave me the wrong results. So I doubt the BCS in my codes are wrong. The following are the new codes,


'''
coordinate:
     x->[0,10]
equations:
    -u1'' -4* u2'*u1' - 8*u1**3+6*u1**5 =0 , where the ' = d/dx,
    -u2''- 2/3 *u1'**2
boundary conditions
    u1(0) =0, u1'(10) =0
    u2(0) =0, u2'(0) =0


'''

import numpy as np

import ufl
from dolfinx import mesh, fem
from dolfinx.fem import Function, FunctionSpace
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, split, TestFunctions

from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType

# parmeters
r_0 = 0.0
r_b = 10.0  # the coordinate of right boundary


a=0.4
b=2.0
c=1.0
beta=0.2
msh = mesh.create_interval(comm=MPI.COMM_WORLD, points=(r_0, r_b), nx=100)

fdim = msh.topology.dim - 1
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
V = FunctionSpace(msh, ufl.MixedElement(P1, P1))

# the left boundary
facets_l = mesh.locate_entities_boundary(
    msh, dim=fdim, marker=lambda x: np.isclose(x[0], r_0))

dofs_u1 = fem.locate_dofs_topological(V=V.sub(0), entity_dim=fdim, entities=facets_l)
dofs_u2 = fem.locate_dofs_topological(V=V.sub(1), entity_dim=fdim, entities=facets_l)

# the Dirichleb boundary conditions
bc_u1 = fem.dirichletbc(value=ScalarType(0), dofs=dofs_u1, V=V.sub(0))  # u1 = phi
bc_u2 = fem.dirichletbc(value=ScalarType(0), dofs=dofs_u2, V=V.sub(1))  # u2 = A

bcs = [bc_u1, bc_u2]

# Define functions
u = Function(V)
# Split system functions to access components
#phi, m, sigma = u1, u2, u3
u1, u2 = split(u)

# Define test functions
v_1, v_2 = TestFunctions(V)

# Define variational problem

F1 = v_1.dx(0)*u1.dx(0) + (-4.0 * u2.dx(0)*u1.dx(0) +2.0*a*u1- 4.0*b*u1**3 + 6.0*c*u1**5)*v_1
F2 = v_2.dx(0)*u2.dx(0) - 2.0/3.0*u1.dx(0)**2*v_2
F = (F1+F2)*dx

#  Solve the system.
#
# t1 = data[:,1]
# t2 = data[:,2]
# t3 = data[:,3]



x = np.linspace(r_0,r_b,101)
def guess(x):
      
    phi0 = np.sqrt((np.sqrt(b**2-3*a*c)+b)/3*c)
 	#print phi0
 	#print a,b,c
    t1=2
    y1_ = np.zeros(101)
    y2_ = np.zeros(101)
    j=0
    for i in x: 
        if i <= t1:
            y1_[j] = 0.0
            y2_[j] = 0.0	
            
        else:
            y1_[j]=phi0
            y2_[j]=-0.3*i	
        j=j+1
    return y1_ ,y2_ 
yguess=guess(x)
u.x.array[:] = np.r_[yguess[0],yguess[1]]

# Create nonlinear problem and Newton solver
problem = NonlinearProblem(F, u, bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-5
solver.nonzero_initial_guess = True
solver.max_it = 100
solver.report = True

# We can customize the linear solver used inside the NewtonSolver by
# modifying the PETSc options
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
# -

r = solver.solve(u)


import matplotlib.pyplot as plt
u1,u2=split(u)
plt.figure(1,dpi=300)  

plt.plot ( u.x.array[0:100])
plt.plot ( u.x.array[101:201])
filename = 'bvp_02_solution.png'
plt.savefig ( filename )
print ( '  Graphics saved as "%s"' % ( filename ) )
plt.show()
plt.close ( )








The right results shoud be as

@zhaozhh I also only started to use FEniCSx a week ago and didn’t get my nonlinear 2-field problem running. Can’t help you unfortunately. I normally have extensive log file output during the code development stage which helps me getting things working in a conventional FEM development environment. FEniCSx is very high-level which is its strength but it also seems to make effective debugging difficult.

Dear @skatulla, thank you for your advice, I will try to use the log file in my codes.

I think this solution is not as wrong as you think. When you plot the solution, you are making an incorrect assumption about the dof layout (i.e. you are assuming that the dofs corresponding of u1 are stored in positions 0 through 99 of u.x, and the dofs of u2 are stored in entries 101 through 200 of u.x):

However, as you can see by using

V0, dofs0 = V.sub(0).collapse()
V1, dofs1 = V.sub(1).collapse()
print(dofs0)
print(dofs1)

this is not how the dofs are arranged. For example, I obtain

>>> V0, dofs0 = V.sub(0).collapse()
>>> dofs0
[0, 1, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42,
    44, 46, 48, 50, 52, 54, 56, 58, 60, 62, 64, 66, 68, 70, 72, 74, 76, 78, 80, 82,
    84, 86, 88, 90, 92, 94, 96, 98, 100, 102, 104, 106, 108, 110, 112, 114, 116, 118,
    120, 122, 124, 126, 128, 130, 132, 134, 136, 138, 140, 142, 144, 146, 148, 150,
    152, 154, 156, 158, 160, 162, 164, 166, 168, 170, 172, 174, 176, 178, 180, 182,
    184, 186, 188, 190, 192, 194, 196, 198, 200]
>>> V1, dofs1 = V.sub(1).collapse()
>>> dofs1
[2, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 35, 37, 39, 41, 43,
    45, 47, 49, 51, 53, 55, 57, 59, 61, 63, 65, 67, 69, 71, 73, 75, 77, 79, 81, 83,
    85, 87, 89, 91, 93, 95, 97, 99, 101, 103, 105, 107, 109, 111, 113, 115, 117, 119,
    121, 123, 125, 127, 129, 131, 133, 135, 137, 139, 141, 143, 145, 147, 149,
    151, 153, 155, 157, 159, 161, 163, 165, 167, 169, 171, 173, 175, 177, 179,
    181, 183, 185, 187, 189, 191, 193, 195, 197, 199, 201]

You should use something like the following for plotting.

V0, dofs0 = V.sub(0).collapse()
V1, dofs1 = V.sub(1).collapse()
plt.plot ( u.x.array[dofs0])
plt.plot ( u.x.array[dofs1])

When I export the result to XDMF and plot with ParaView, I obtain:

from dolfinx import io
u1, u2 = u.split()
with (io.XDMFFile(MPI.COMM_WORLD, "u1.xdmf", "w") as xdmf,
        io.XDMFFile(MPI.COMM_WORLD, "u2.xdmf", "w") as xdmf2):
    xdmf.write_mesh(msh)
    xdmf.write_function(u1, 0.0)
    xdmf2.write_mesh(msh)
    xdmf2.write_function(u2, 0.0)

2 Likes

It appears to me that this only fixes wf on the top and bottom faces, rather than fixing all the dofs as you suggest. To fix all the dofs, I think you should try

# Fix ALL change of director degrees of freedom in the domain
def all_dofs(x):
    return np.full((x.shape[1],), True)
V1, submap = VY.sub(1).collapse()
fixed_director = fem.Function(V1)
fixed_director.interpolate(fixed_director_expression)
all_wf_dofs = fem.locate_dofs_geometrical((VY.sub(1),V1), all_dofs)
bc3 = fem.dirichletbc(fixed_director, all_wf_dofs, V1)

bcs = [bc1,  bc2, bc3]
2 Likes

@conpierce8 , thank you very much! Your response give me a right understanding about the concepts of dofs and collapse(). Now I’m sure that my BCS is right. And I think the wrong results maybe relate to the problem itself.

u2 holds the Dirichlet and Neumman conditions at the same end point.

u2(0) =0, u2'(0) =0

My Dirichlet conditions is u2(0) =0, which maybe let fenicsx understant the BCs are

u2(0) =0, u2'(10) =0

So fenicsx gives the result

But scipy.integrate.solve_bvp gives the right one

You are correct, FEniCSx (and the finite element method in general) implicitly assumes a zero-flux (i.e. u_2' = 0) boundary condition on all non-Dirichlet boundaries where the flux is not specified. FEniCSx interprets your problem in this way since you have not included a boundary integral (ds) in your variational form.

Offhand, I don’t know how to specify multiple boundary conditions (i.e. Dirichlet and Neumann) at a single point.

Dear @conpierce8 , I have solved the problem

just now by adding a function and equation with following method

u3  = u2'
u3' = -2/3 u1'**2

So

u2(0)=0,
u3(0)=0.

Namely, the u2’s Neumann BC has been changed into the Dirichlet BC of u3.

Following are the changed codes,


'''
coordinate:
     x->[0,10]
equations:
    -u1'' -4* u2'*u1' - 8*u1**3+6*u1**5 =0 , where the ' = d/dx,
   # -u2''- 2/3 *u1'**2 =0 splited into
     u3' +  2/3 *u1'**2 =0
     u3  =  u2'
boundary conditions
    u1(0) =0, # u1'(10) =0
    u2(0) =0, #u2'(0) =0
    u3(0) =0


'''

import numpy as np

import ufl
from dolfinx import mesh, fem,log
from dolfinx.fem import Function, FunctionSpace
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, split, TestFunctions

from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
import matplotlib.pyplot as plt
# parmeters
r_0 = 0.0
r_b = 10.0  # the coordinate of right boundary


a=0.2
b=2.0
c=1.0
beta=0.2
msh = mesh.create_interval(comm=MPI.COMM_WORLD, points=(r_0, r_b), nx=100)

fdim = msh.topology.dim - 1
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
V = FunctionSpace(msh, ufl.MixedElement(P1, P1,P1))

# the left boundary
facets_l = mesh.locate_entities_boundary(msh, dim=fdim, marker=lambda x: np.isclose(x[0], r_0))

dofs_u1 = fem.locate_dofs_topological(V=V.sub(0), entity_dim=fdim, entities=facets_l)
dofs_u2 = fem.locate_dofs_topological(V=V.sub(1), entity_dim=fdim, entities=facets_l)
dofs_u3 = fem.locate_dofs_topological(V=V.sub(2), entity_dim=fdim, entities=facets_l)

# the Dirichleb boundary conditions
bc_u1 = fem.dirichletbc(value=ScalarType(0), dofs=dofs_u1, V=V.sub(0))  # u1 = phi
bc_u2 = fem.dirichletbc(value=ScalarType(0), dofs=dofs_u2, V=V.sub(1))  # u2 = A
bc_u3 = fem.dirichletbc(value=ScalarType(0), dofs=dofs_u3, V=V.sub(2))  # u2 = A

bcs = [bc_u1, bc_u2, bc_u3]

# Define functions
u = Function(V)
# Split system functions to access components
#phi, m, sigma = u1, u2, u3
u1, u2,u3 = split(u)

# Define test functions
v1, v2, v3 = TestFunctions(V)

# Define variational problem

F1 = v1.dx(0)*u1.dx(0) + (-4.0 * u2.dx(0)*u1.dx(0) +2.0*a*u1- 4.0*b*u1**3 + 6.0*c*u1**5)*v1
F2 = (u3.dx(0)+ 2.0/3.0*u1.dx(0)**2)*v2
F3 = (u3-u2.dx(0))*v3
F = (F1+F2+F3)*dx

#  Solve the system.
#
# t1 = data[:,1]
# t2 = data[:,2]
# t3 = data[:,3]



x = np.linspace(r_0,r_b,101)
def guess(x):
      
    t1=2
    y1_ = np.zeros(101)
    y2_ = np.zeros(101)
    y3_ = np.zeros(101)
    j=0
    for i in x: 
        if i <= t1:
            y1_[j] = 0.0
            y2_[j] = 0.0	
            y3_[j] = 0.0	
            
        else:
            y1_[j] = 1
            y2_[j]=-0.2*i	
            y3_[j]=-0.2
        j=j+1
    return y1_ ,y2_ ,y3_
yguess=guess(x)
V0, dofs0 = V.sub(0).collapse()
V1, dofs1 = V.sub(1).collapse()
V2, dofs2 = V.sub(2).collapse()
u.x.array[dofs0]=yguess[0]
u.x.array[dofs1]=yguess[1]
u.x.array[dofs2]=yguess[2]

#u.x.array[:] = np.r_[yguess[0],yguess[1]]
plt.figure(2,dpi=300) 
plt.plot(yguess[0])
plt.plot(yguess[1])
plt.plot(yguess[2])
plt.show()
# Create nonlinear problem and Newton solver
problem = NonlinearProblem(F, u, bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-5
solver.nonzero_initial_guess = True
solver.max_it = 1000
solver.report = True

# We can customize the linear solver used inside the NewtonSolver by
# modifying the PETSc options
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
# -
#log.set_log_level(log.LogLevel.INFO)
r = solver.solve(u)



#1_,u2_, u3_=u.split()


plt.figure(2,dpi=300)

plt.plot (x, u.x.array[dofs0],label="u1")
plt.plot (x, u.x.array[dofs1],label="u2")
plt.legend()
plt.show()

plt.close ( )

The right result is obtained at last.

1 Like

@conpierce8 Thanks for the suggestion. In principle, you are correct. However, in this special case of choosing a discretization of a single 8-node brick element in both thickness directions of this beam problem, fixing the DOF of all nodes on top and bottom surfaces has the same effect.

@skatulla Ah! You caught me not paying close attention. :grimacing:

It looks like you should not apply your Dirichlet bcs to the collapsed subspaces, i.e. replacing

bc1 = fem.dirichletbc(fixed_displacement, left_u_dofs, V0)
bc2 = fem.dirichletbc(incremented_displacement, right_u_dofs, V0)

bc3 = fem.dirichletbc(fixed_director, top_wf_dofs, V1)
bc4 = fem.dirichletbc(fixed_director, bottom_wf_dofs,V1)

with

bc1 = fem.dirichletbc(fixed_displacement, left_u_dofs, VY.sub(0))
bc2 = fem.dirichletbc(incremented_displacement, right_u_dofs, VY.sub(0))

bc3 = fem.dirichletbc(fixed_director, top_wf_dofs, VY.sub(1))
bc4 = fem.dirichletbc(fixed_director, bottom_wf_dofs,VY.sub(1))

achieves convergence for me.

Complete code:

from contextlib import ExitStack

import numpy as np

import ufl
from dolfinx import fem, nls, mesh, plot, log, io
from mpi4py import MPI
from petsc4py import PETSc

import pdb


# ----------------------------------------------------------------------------------------------------------------------
# Create a box mesh of a beam
L = 20.0
H = 1.0
msh = mesh.create_box(MPI.COMM_WORLD, [np.array([0.0, 0.0, 0.0]), np.array([L, H, H])], [10, 1, 1],
                      mesh.CellType.hexahedron)
P1 = ufl.VectorElement("CG", msh.ufl_cell(), 1, dim=3)  # Lagrange family, 1st order
P2 = ufl.VectorElement("CG", msh.ufl_cell(), 1, dim=3)  # Lagrange family, 1st order
VY = fem.FunctionSpace(msh, ufl.MixedElement([P1, P2]))


# ----------------------------------------------------------------------------------------------------------------------
# Boundary conditions and body loads
def left(x):
    return np.isclose(x[0], 0)

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

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

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


def fixed_displacement_expression(x):
    return np.stack((np.zeros(x.shape[1]), np.zeros(x.shape[1]), np.zeros(x.shape[1])))

def incremented_displacement_expression(x):
    return np.stack((np.full(x.shape[1], 1.0e-02), np.zeros(x.shape[1]), np.zeros(x.shape[1])))

def fixed_director_expression(x):
    return np.stack((np.zeros(x.shape[1]), np.zeros(x.shape[1]), np.zeros(x.shape[1])))

# ----------------------------------------------
# Dirichlet boundary application to affected displacement dofs

# Fix displacement DOF on left face and apply axial displacement on right face
V0, submap = VY.sub(0).collapse()
fixed_displacement = fem.Function(V0)
fixed_displacement.interpolate(fixed_displacement_expression)
incremented_displacement = fem.Function(V0)
incremented_displacement.interpolate(incremented_displacement_expression)
left_u_dofs = fem.locate_dofs_geometrical((VY.sub(0),V0), left)
right_u_dofs = fem.locate_dofs_geometrical((VY.sub(0),V0), right)
bc1 = fem.dirichletbc(fixed_displacement, left_u_dofs, VY.sub(0))
bc2 = fem.dirichletbc(incremented_displacement, right_u_dofs, VY.sub(0))

# Fix ALL change of director degrees of freedom in the domain
V1, submap = VY.sub(1).collapse()
fixed_director = fem.Function(V1)
fixed_director.interpolate(fixed_director_expression)
top_wf_dofs = fem.locate_dofs_geometrical((VY.sub(1),V1), top)
bottom_wf_dofs = fem.locate_dofs_geometrical((VY.sub(1),V1), bottom)
bc3 = fem.dirichletbc(fixed_director, top_wf_dofs, VY.sub(1))
bc4 = fem.dirichletbc(fixed_director, bottom_wf_dofs,VY.sub(1))

bcs = [bc1,  bc2, bc3, bc4]
pdb.set_trace()

# ----------------------------------------------------------------------------------------------------------------------
# Function space
v, yf = ufl.TestFunctions(VY)     # test functions
uwf = fem.Function(VY)            # current displacement and change of director
u, wf = ufl.split(uwf)

# ----------------------------------------------------------------------------------------------------------------------
# Kinematics

# Spatial dimension
d = len(u)

# Identity tensor
I = ufl.variable(ufl.Identity(d))

# Deformation gradient
F = ufl.variable(I + ufl.grad(u))

# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)

# Invariants
Ic = ufl.variable(ufl.tr(C))
J  = ufl.variable(ufl.det(F))

# ----------------------------------------------------------------------------------------------------------------------
# Constitutive law

# Set the elasticity parameters
E = PETSc.ScalarType(1.0e4)
nu = PETSc.ScalarType(0.3)
mu = fem.Constant(msh, E/(2*(1 + nu)))
lmbda = fem.Constant(msh, E*nu/((1 + nu)*(1 - 2*nu)))

# Stored strain energy density (compressible neo-Hookean model)
psi_m = (mu / 2.0) * (Ic - 3.0) - mu * ufl.ln(J) + (lmbda / 2.0) * (ufl.ln(J))**2

# 1st Piola Kirchhoff stress tensor
P = ufl.diff(psi_m, F)                                     # 1st Piola Kirchhoff stress tensor
df = fem.Constant(msh, PETSc.ScalarType((0.0, 0.0, 0.0)))  # dummy director stress

# ----------------------------------------------------------------------------------------------------------------------
# Variational formulation
metadata = {"quadrature_degree": 3}
dx = ufl.Measure("dx", domain=msh, metadata=metadata)

# Define form Pi (we want to find u such that Pi(u) = 0)
Pi = ufl.inner(P, ufl.grad(v))*dx + ufl.inner(df, yf)*dx

# Initialize nonlinear problem
problem = fem.petsc.NonlinearProblem(Pi, uwf, bcs)

# ----------------------------------------------------------------------------------------------------------------------
# Initialise Newton Rhapson solver
solver = nls.petsc.NewtonSolver(msh.comm, problem)

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

# ----------------------------------------------------------------------------------------------------------------------
# Incremental application of the traction loading
log.set_log_level(log.LogLevel.INFO)
for n in range(1, 2):
    num_its, converged = solver.solve(uwf)
    assert(converged)

@conpierce8 Ah, very nice! I see your small changes defining the Dirichlet boundary conditions, instead of

P1 = ufl.VectorElement("CG", msh.ufl_cell(), 1, dim=3)  # Lagrange family, 1st order
P2 = ufl.VectorElement("CG", msh.ufl_cell(), 1, dim=3)  # Lagrange family, 1st order
VY = fem.FunctionSpace(msh, ufl.MixedElement([P1, P2]))
bc1 = fem.dirichletbc(fixed_displacement, left_u_dofs, V0)
bc2 = fem.dirichletbc(incremented_displacement, right_u_dofs, V0)
bc3 = fem.dirichletbc(fixed_director, top_wf_dofs, V1)
bc4 = fem.dirichletbc(fixed_director, bottom_wf_dofs,V1)

you changed them to

P1 = ufl.VectorElement("CG", msh.ufl_cell(), 1, dim=3)  # Lagrange family, 1st order
P2 = ufl.VectorElement("CG", msh.ufl_cell(), 1, dim=3)  # Lagrange family, 1st order
VY = fem.FunctionSpace(msh, ufl.MixedElement([P1, P2]))
bc1 = fem.dirichletbc(fixed_displacement, left_u_dofs, VY.sub(0))
bc2 = fem.dirichletbc(incremented_displacement, right_u_dofs, VY.sub(0))
bc3 = fem.dirichletbc(fixed_director, top_wf_dofs, VY.sub(1))
bc4 = fem.dirichletbc(fixed_director, bottom_wf_dofs,VY.sub(1))

Thanks a lot!

1 Like