Boundary condition to single direction in the subspace of Mixed function space

[kumar]

Hello everyone,
In my code, I wanted to implement a displacement boundary constraint on the right bottom corner point as fixed only along y direction in a rectangular domain. Since it is a mixed element function space, I am little unsure about the code. I have checked for some previous posts but still remains unclear.

Could anyone please confirm me that the applied boundary condition for the right bottom corner point of the beam is constrained along y direction? . My code is

import os
from mpi4py import MPI
from petsc4py import PETSc
import dolfinx
import numpy as np
import ufl
from basix.ufl import element
from dolfinx import default_real_type, log, plot
from dolfinx.fem import Function, functionspace,assemble_scalar
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.io import XDMFFile
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner
from dolfinx import fem, default_scalar_type,mesh

#Creating mesh
msh = mesh.create_rectangle (comm=MPI.COMM_WORLD, points=((0.0, 0.0), (6.0,1.0)), n=(1,1), cell_type=mesh.CellType.quadrilateral)

#Create function space
V = ufl.VectorElement("Lagrange", msh.ufl_cell(),1) #For deformation map
P1 = element("Lagrange", msh.basix_cell(), 1) #For phase field parametrs

#Creating mixed element space for phi,p,mu,lamda1,lamda2
ME = functionspace(msh, ufl.MixedElement([V, P1, P1, P1, P1]))

#Creating Test functions and function
a, q, v, l1, l2 = ufl.TestFunctions(ME)
u = Function(ME)  # current solution
u0 = Function(ME)  # solution from previous converged step

# Split mixed functions
d, p, mu, lam1, lam2 = ufl.split(u)
d0, p0, mu0, lam01, lam02 = ufl.split(u0)
# Zero u
u.x.array[:] = 0.0

# Interpolate initial condition for the concentration
np.random.seed(30)
u.sub(1).interpolate(lambda x: 0.75 + 0.02 * (0.5 - np.random.rand(x.shape[1])))
u.x.scatter_forward()

#Applying boundary condition to beam
def left_corner(x):
    return np.logical_and(np.isclose(x[1],0.0),np.isclose(x[0],0.0))
    
ME0, _ = ME.sub(0).collapse()
bottom_value1 = fem.Function(ME0)
bottom_value1.x.array[:]= 0
left_corner_vertex = mesh.locate_entities_boundary(msh,0,left_corner)
dofs = fem.locate_dofs_topological((ME.sub(0),ME0),0,left_corner_vertex)
bc_left_corner = fem.dirichletbc(bottom_value1,dofs,ME0)
 
########## 
def right_corner(x):
    return np.logical_and(np.isclose(x[1],0.0),np.isclose(x[0],6.0))
    
bottom_value2 = fem.Function(ME0)
bottom_value2.x.array[:]= 0
right_corner_vertex = mesh.locate_entities_boundary(msh,0,right_corner)
right_dof = fem.locate_dofs_topological ((ME.sub(0),ME0.sub(1)),0,right_corner_vertex) #Having doubt regarding the bc on right corner 
#to constrain displacment only along y direction.
bc_right_corner = fem.dirichletbc(bottom_value2, right_dof, ME.sub(0))
##########

Since I need to include boundary condition to only along y direction to subspace of mixed element function space, it is confusing (right_corner function in code). So it will be great help if someone confirm me regarding the implementation.

Thanks for your time and considerations.

[kumar]

Additionally in the code, when I tried to run, it ends up in showing error, ‘RuntimeError: Failed to successfully call PETSc function ‘KSPSolve’. PETSc error code is: 76, Error in external library’. Any help will be appreciated. Thanks for your time and considerations.
My completed code is,

import os
from mpi4py import MPI
from petsc4py import PETSc
import dolfinx
import numpy as np
import ufl
from basix.ufl import element
from dolfinx import default_real_type, log
from dolfinx.fem import Function, functionspace,assemble_scalar
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.io import XDMFFile
from dolfinx.nls.petsc import NewtonSolver
from ufl import ds, dx, grad, inner
#Additional for mechanical 
from dolfinx import fem, default_scalar_type,mesh

#Creating mesh
msh = mesh.create_rectangle (comm=MPI.COMM_WORLD, points=((0.0, 0.0), (6.0,1.0)), n=(1,1), cell_type=mesh.CellType.quadrilateral)

#Create function space
V = ufl.VectorElement("Lagrange", msh.ufl_cell(),1) #For deformation map
P1 = element("Lagrange", msh.basix_cell(), 1) #For phase field parametrs

#Creating mixed element space for phi,p,mu,lamda1,lamda2
ME = functionspace(msh, ufl.MixedElement([V, P1, P1, P1, P1]))

#Creating Test functions and function
a, q, v, l1, l2 = ufl.TestFunctions(ME)
u = Function(ME)  # current solution
u0 = Function(ME)  # solution from previous converged step

# Split mixed functions
d, p, mu, lam1, lam2 = ufl.split(u)
d0, p0, mu0, lam01, lam02 = ufl.split(u0)
# Zero u
u.x.array[:] = 0.0

# Interpolate initial condition for the concentration
np.random.seed(30)
u.sub(1).interpolate(lambda x: 0.75 + 0.02 * (0.5 - np.random.rand(x.shape[1])))
u.x.scatter_forward()

#Applying boundary condition to beam
def left_corner(x):
    return np.logical_and(np.isclose(x[1],0.0),np.isclose(x[0],0.0))
    
ME0, _ = ME.sub(0).collapse()
bottom_value1 = fem.Function(ME0)
bottom_value1.x.array[:]= 0
left_corner_vertex = mesh.locate_entities_boundary(msh,0,left_corner)
dofs = fem.locate_dofs_topological((ME.sub(0),ME0),0,left_corner_vertex)
bc_left_corner = fem.dirichletbc(bottom_value1,dofs,ME.sub(0))
  
def right_corner(x):
    return np.logical_and(np.isclose(x[1],0.0),np.isclose(x[0],6.0))
    
bottom_value2 = fem.Function(ME0)
bottom_value2.x.array[:]= 0
right_corner_vertex = mesh.locate_entities_boundary(msh,0,right_corner)
right_dof = fem.locate_dofs_topological ((ME.sub(0),ME0.sub(1)),0,right_corner_vertex) 
bc_right_corner = fem.dirichletbc(bottom_value2, right_dof, ME.sub(0))

bcs = [bc_left_corner, bc_right_corner]

                                                    ##### Mechanical solver terms #####
#Defining traction and body forces
T=fem.Constant(msh,default_scalar_type((0,0)))
B=fem.Constant(msh,default_scalar_type((0,0)))

#Defining kinematic quantities
#Spatial dimension
dim=len(d)
#Identity tensor
I=ufl.variable(ufl.Identity(dim))
#Deformation gradient
F=ufl.variable(I+grad(d))
# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)
# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))

# Elasticity parameters
E =default_scalar_type(1e4)
nu = default_scalar_type(0.3)
mu_mech = 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 = (mu_mech / 2) * (Ic - 3) - mu_mech * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2
# Hyper-elasticity
P = ufl.diff(psi, F)

                                                    ##### Phase-field solver terms #####
#Defining parameters
epsilon = 0.625  # interface thickness
areainteng = 21 #Area specific interface energy
dt = 5.0e-06  # time step
k = 0.05 #Dissipation term
theta = 0.5  # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicholson

# Compute the chemical potential df/dc
p = ufl.variable(p)
#f = 100 * p**2 * (1 - p)**2 #Double well potential tutorial
f = areainteng*(((3/2)*epsilon*(grad (p))**2)+(6*(1/epsilon)*p**2*(1-p)**2))
dfdc = ufl.diff(f, p)
# mu_(n+theta)  
mu_mid = (1.0 - theta) * mu0 + theta * mu

#Weak form equations
F0 = inner(grad(a), P) * dx - inner(a, B) * dx - inner(a, T) * ds
F1 = inner(p, q) * dx - inner(p0, q) * dx + dt * k * inner(grad(mu_mid), grad(q)) * dx
F2 = inner(mu, v) * dx - inner(dfdc, v) * dx - 3*epsilon*areainteng*inner(grad(p), grad(v)) * dx
FF = F0+F1+F2

#Solving Nonlinear problem
problem=NonlinearProblem(FF,u,bcs)
#Solver configurations
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"

# We can customize the linear solver used inside the NewtonSolver by
# modifying the PETSc options
ksp = solver.krylov_solver
opts = PETSc.Options()  # type: ignore
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
ksp.setFromOptions()
solver.rtol = np.sqrt(np.finfo(default_real_type).eps) * 1e-2

# Step in time
t = 0.0
#  Reduce run time if on test (CI) server
if "CI" in os.environ.keys() or "GITHUB_ACTIONS" in os.environ.keys():
    T = 1 * dt
else:
    T = 5 * dt

while (t < T):
    t += dt
    r = solver.solve(u)
    print(f"Step {int(t/dt)}: num iterations: {r[0]}")
    u0.x.array[:] = u.x.array

Could anyone please confirm, whether my error is in boundary condition or some where else. Thanks again in advance.

[dokken]

@kumar

left_corner_vertex = mesh.locate_entities_boundary(msh,0,left_corner)
dofs = fem.locate_dofs_topological((ME.sub(0),ME0),0,left_corner_vertex)
bc_left_corner = fem.dirichletbc(bottom_value1,dofs,ME0)

You are taking the first sub-space of the mixed element, which is V, not V.sub(0)

@kumar

ME = functionspace(msh, ufl.MixedElement([V, P1, P1, P1, P1]))

You can check this by printing ME.element.num_sub_elements, which gives you 5, not 6.
Thus, you would have to change:

@kumar

right_dof = fem.locate_dofs_topological ((ME.sub(0),ME0.sub(1)),0,right_corner_vertex) #Having doubt regarding the bc on right corner 
#to constrain displacment only along y direction.
bc_right_corner = fem.dirichletbc(bottom_value2, right_dof, ME.sub(0))

To

ME0y, _ = ME.sub(0).sub(1).collapse()
bottom_value2 = fem.Function(ME0y)
bottom_value2.x.array[:]= 0
right_corner_vertex = mesh.locate_entities_boundary(msh,0,right_corner)
right_dof = fem.locate_dofs_topological ((ME.sub(0).sub(1),ME0y),0,right_corner_vertex)
bc_right_corner = fem.dirichletbc(bottom_value2, right_dof, ME.sub(0).sub(1))

[kumar]

Thanks for your immediate response Mr.Dokken. I have a general question which is very much necessary for my problem. In my case of coupled problem of Cahn-Hilliard and Mechanics, I need to solve total of 5 set of equations. Out of it, three were needed to be solved for the whole domain and the final two equations were solved only at the nodal points whereas which need not to be interpolated to the whole domain.

As a part of solution, I need to assemble my elemental stiffness matrix where the one part needed to calculate at the whole domain level where interpolation is done (shown in grey in Figure) and the other part, need not to be interpolated which in turn solved only at the nodes, so it is coupled directly to the assembly. Because integral is mentioned for the expression mentioned in grey, but the other are just scalar expressions of the constraints, needed to add in the final assembly as per the instruction. The elemental stiffness matrix is shown in Figure below.

. As shown in the above figure the components not in the grey, are not needed to be interpolated and are scalar expression values directly added to the matrix explicitly. Like this it is needed to be added explicitly for all the elements in domain.

So could you please confirm me, whether such an assembly is possible in FEniCSx or any other possible approaches to assemble like above and solve the problem.

Thanks for your time and considerations. Any help on the above is truly appreciated and thanks everyone in advance. I am expecting for all your help and suggestions.

[dokken]

I don’t quite understand what the two different parts of your matrix is supposed to be here. However, if you can express each of them in DOLFINx, you can use the PETSc nest functionality to make a block matrix, see for instance; Stokes equations using Taylor-Hood elements — DOLFINx 0.7.3 documentation

[kumar]

I hope the equations below makes a bit clear understanding of my problem. The grey part of the assembly were obtained from equations,

So these are the ones can be obtained by our normal procedure by solving in FEniCSx since these are procedure involved normally in solving weak form equations but the others are not in integral and just expressions only at nodes and not interpolated to whole domain. Those expression are,

.
These are only only at nodes and not for the whole domain. So these are directly needed to be added in the final assembled element matrix along with the grey ones. These are expanded as ,

I hope these makes a clear understanding. I did not getting an idea for solving these elemental matrix with half at domain level and other at only 4 nodal points without interpolation at gauss points.

Could you please suggest some ideas for solving these. Thanks for your time and considerations.

[dokken]

You can create the remaining matrices with petsc4py, and use them with petsc-nest

[kumar]

Thanks for your immediate response. So from my understanding, my weak form equations (grey box) were solved normally using FEniCsx mixed function space and the rest two constraint equations (Fischer-Burmeister condition) were created using petsc4py and then both are them are added using the petsc-nest ?
Such that after assembly, I can further proceed to solve normally in FEniCSx. Is my understanding of the solution procedure is correct? Could you please confirm?.

Many thanks for your time and considerations.

[dokken]

Yes, that would be the idea

[kumar]

I can clarify my question with a simple example of the Cahn-Hilliard tutorial. In the following code, I am solving the weak form equations of the phase field as F0 + F1. The solver part of the code looks like,

# Weak statement of the equations
F0 = inner(c, q) * dx - inner(c0, q) * dx + dt * inner(grad(mu_mid), grad(q)) * dx
F1 = inner(mu, v) * dx - inner(dfdc, v) * dx - lmbda * inner(grad(c), grad(v)) * dx
F = F0 + F1   

# Create nonlinear problem and Newton solver
problem = NonlinearProblem(F, u)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = np.sqrt(np.finfo(default_real_type).eps) * 1e-2

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

Confusing part for me is that, can I be able to add my new petsc matrix to the assembled elemental stiffness matrix of this F ? Because usually we directly give solver.solve(), once the solver has been modified. So in between can I able to control the matrix of elemental stiffness matrix modification?

It will be very much helpful, if I get some hint, so that I can proceed my problem further. Any hints for proceeding my problem is greatly appreciated. Many thanks for all your time and considerations.