# How to applying BC using V..sub().collapse() in dolfinx

Hello all,

I am solving a simple elasticity problem using dolfinx and I am trying to fix one side of a beam and let it deflect under its weight. I am trying the use the “V.sub(2).collapse()[0]” for that. There is no error, but it is not working for some reasons. it shows no deformation at all.

Here is the simple code that I am running :

``````from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import pyvista

from dolfinx.fem import (Constant, Function, FunctionSpace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical,
locate_dofs_topological, assemble, locate_dofs_topological, VectorFunctionSpace)
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc, LinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_unit_square, locate_entities, meshtags, locate_entities_boundary
from dolfinx import mesh as Mesh
from petsc4py.PETSc import ScalarType
import meshio
from dolfinx.plot import create_vtk_mesh
from ufl import (FacetNormal, FiniteElement, TensorElement, Identity, TestFunction, TrialFunction, VectorElement,
div, dot, dx, inner, lhs, nabla_grad, nabla_div, rhs, sym, indices, as_tensor, Measure, MixedElement,
split)

L = 1
W = 0.2
mu = 1
rho = 1
delta = W / L
gamma = 0.4 * delta ** 2
beta = 1.25
lambda_ = beta
g = gamma

mesh = Mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, W])],
[20, 6, 6], cell_type=Mesh.CellType.hexahedron)

element = VectorElement('Lagrange', mesh.ufl_cell(), degree=1, dim=3)
V = FunctionSpace(mesh, element)

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

fdim = mesh.topology.dim - 1

V0, submap = V.sub(0).collapse()

left_dofs = locate_dofs_geometrical((V.sub(2), V0), clamped_boundary)

bc= dirichletbc(ScalarType(0.0), left_dofs[0], V.sub(2))

# ------------------------------------------
T = Constant(mesh, ScalarType((0, 0, 0)))

ds = Measure("ds", domain=mesh)

def epsilon(u):

def sigma(u):
return lambda_ * nabla_div(u) * Identity(len(u)) + 2 * mu * epsilon(u)

u = TrialFunction(V)
v = TestFunction(V)
f = Constant(mesh, ScalarType((0, 0, -rho * g)))
a = inner(sigma(u), epsilon(v)) * dx
L = dot(f, v) * dx + dot(T, v) * ds

problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uf = problem.solve()

# ......................................... Export Results
xdmf = XDMFFile(mesh.comm, "DIS_Beam.xdmf", "w")
xdmf.write_mesh(mesh)
xdmf.write_function(uf
``````

Hi,
Thank you so much for the reply. I solved the problem, but my over all goal is to transform a working code in fenics to fenicsx(0.6.0).
The problem is that I don’t get any error in fenicsx but the results is not correct. Actually, the results is zero all over the domain. Since there is no error, I don’t know how should I find the problem in the fenicsx code: I would appreciate any comments on that.

Here is the code in fenics:

``````from dolfin import *
import fenics as fe
# from fenics import *
# import matplotlib.pyplot as plt
import numpy as np
# from ufl import *
from ufl import indices

# ......................................... Parameters
E = 1 * 152e9
nu = 0.33
la = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
perm = 12.48e-9
epsilon = 11e-9
non_local = 1e-20
mu11 = 1e-6
mu12 = 0.0
mu44 = 0.0
# ................
ESi = 1 * 152e9
laSi = ESi * nu / ((1 + nu) * (1 - 2 * nu))
muSi = ESi / (2 * (1 + nu))
mesh = fe.Mesh("Beam.xml")
Thickness = 20e-9
L = 1000e-9
H = 100e-9
M = 1e-15
Force = 1e-15  # M / (0.5 * L)
Area = 100e-9 ** 2
tr_force = Force / Area
tr_S1 = fe.Constant((0.0, 0, 0.0))
tr_S2 = fe.Constant((0.0, -tr_force, 0.0))
n = FacetNormal(mesh)

# ....................................... Mark boundary regions
tol = 10e-9

boundaries = fe.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)

Mid2 = AutoSubDomain(lambda x: near(x[1], 50e-9, tol))
Mid2.mark(boundaries, 5)

# ......................................

# ......................................... Define function spaces
V = fe.VectorElement("CG", tetrahedron, 2)  # Displacement
S = fe.TensorElement("CG", tetrahedron, 1)  # Relaxed strain
LM = fe.TensorElement("CG", tetrahedron, 1)  # Lagrange multipliers
P = fe.FiniteElement("CG", tetrahedron, 1)  # Electric Potential
# ......................................... Define mixed element
MIXEL = fe.MixedElement([V, S, LM, P])
Space = FunctionSpace(mesh, MIXEL)
# ......................................... Define Variational Variables
v_tr = TrialFunction(Space)
v_test = TestFunction(Space)
Result = Function(Space)
u, gradu, lamu, p = fe.split(v_tr)
del_u, del_gradu, del_lamu, del_p = fe.split(v_test)
# ......................................... Define measure parameter & BCs
ds = fe.Measure("ds", domain=mesh, subdomain_data=boundaries)

# ......................................... Boundary conditions

b1 = fe.DirichletBC(Space.sub(0).sub(0), fe.Constant(0.0), Left)
b2 = fe.DirichletBC(Space.sub(0).sub(1), fe.Constant(0.0), Left)
b3 = fe.DirichletBC(Space.sub(0).sub(2), fe.Constant(0.0), Left)

bp = fe.DirichletBC(Space.sub(3), fe.Constant(0.0), Mid2)
bc = [b1, b2, b3, bp]
# ......................................... Define tensor
delta = Identity(V.cell().geometric_dimension())
i, j, k, l, m, n = indices(6)

# Strain Tensor
def StrainT(u):
return as_tensor((1. / 2. * (u[i].dx(j) + u[j].dx(i))), [i, j])

# Stress Tensor
def Elas(la, mu):
return as_tensor(la * delta[i, j] * delta[k, l] + 2. * mu * delta[i, k] * delta[j, l],
(i, j, k, l))  # Elasticity tensor

def StressT(u, la, mu):
return as_tensor((Elas(la, mu)[i, j, k, l] * StrainT(u)[k, l]), [i, j])

# Relaxed Strain Tensor
def RStrain(u):
return as_tensor((u[j].dx(i)), [i, j])

# Strain Gradient Tensor using relaxed strain
return as_tensor((1. / 2. * (gradu[k, j].dx(i) + gradu[j, k].dx(i))), [i, j, k])

# Higher Order Stress
def SGElas(la, non_local, mu):
return as_tensor((la * (non_local ** 2) * delta[i, j] * delta[l, m] * delta[k, n] +
2 * mu * (non_local ** 2) * delta[i, l] * delta[j, m] * delta[k, n]),
(i, j, k, l, m, n))  # Strain Gradient Elasticity Tensor

return as_tensor(
(SGElas(la, non_local, mu)[i, j, k, l, m, n] * RStrianGrad(gradu)[l, m, n]),
(i, j, k))

# Flexoelectricity Tensor
fl1 = np.zeros((3, 3, 3, 3))
fl1[0, 0, 0, 0] = 1
fl1[1, 1, 1, 1] = 1
fl1[2, 2, 2, 2] = 1
flexo1 = fe.Constant(fl1)
flexo = as_tensor(((mu11 - mu12 - 2 * mu44) * flexo1[i, j, k, l] + mu12 * delta[i, j] * delta[
k, l] + mu12 * delta[i, k] * delta[j, l] + mu12 * delta[i, l] * delta[j, k]), (i, j, k, l))

def HStressElec(p):
return as_tensor((-1 * flexo[l, i, j, k] * (-1 * grad(p))[l]), (i, j, k))

return as_tensor(perm * delta[i, j] * (-1 * grad(p))[j] + (flexo[i, j, k, l] * RStrianGrad(gradu)[j, k, l]), i)

# ......................................... Dimensionless Form
Fl = 1  # 10 ** (-6)
Fs = 1  # E
NONDIM = fe.Constant(1.0)
NONDIM2 = fe.Constant(1.0)
# ......................................... Define Variational Formulation
F1 = NONDIM * StressT(u, la, mu)[j, k] * StrainT(del_u)[j, k] * dx - NONDIM * lamu[j, k] * del_u[k].dx(
j) * dx - NONDIM * \
tr_S1[i] * del_u[i] * ds(1) - NONDIM * tr_S2[i] * del_u[i] * ds(2)
F2 = NONDIM2 * HStress(gradu)[i, j, k] * RStrianGrad(del_gradu)[i, j, k] * dx + NONDIM * lamu[j, k] * \
j, k] * dx

# j=0 and k=0
C00 = NONDIM * gradu[0, 0] * del_lamu[0, 0] * dx - NONDIM * u[0].dx(0) * del_lamu[0, 0] * dx
# j=0 and k=1
C01 = NONDIM * gradu[0, 1] * del_lamu[0, 1] * dx - NONDIM * u[1].dx(0) * del_lamu[0, 1] * dx
# j=0 and k=2
C02 = NONDIM * gradu[0, 2] * del_lamu[0, 2] * dx - NONDIM * u[2].dx(0) * del_lamu[0, 2] * dx
# j=1 and k=0
C10 = NONDIM * gradu[1, 0] * del_lamu[1, 0] * dx - NONDIM * u[0].dx(1) * del_lamu[1, 0] * dx
# j=1 and k=1
C11 = NONDIM * gradu[1, 1] * del_lamu[1, 1] * dx - NONDIM * u[1].dx(1) * del_lamu[1, 1] * dx
# j=1 and k=2
C12 = NONDIM * gradu[1, 2] * del_lamu[1, 2] * dx - NONDIM * u[2].dx(1) * del_lamu[1, 2] * dx
# j=2 and k=0
C20 = NONDIM * gradu[2, 0] * del_lamu[2, 0] * dx - NONDIM * u[0].dx(2) * del_lamu[2, 0] * dx
# j=2 and k=1
C21 = NONDIM * gradu[2, 1] * del_lamu[2, 1] * dx - NONDIM * u[1].dx(2) * del_lamu[2, 1] * dx
# j=2 and k=2
C22 = NONDIM * gradu[2, 2] * del_lamu[2, 2] * dx - NONDIM * u[2].dx(2) * del_lamu[2, 2] * dx

VarForm = F1 + F2 + F3 + C00 + C01 + C02 + C10 + C11 + C12 + C20 + C21 + C22  # + F4
# ......................................... Preparing the form
F = VarForm
a = lhs(F)
L = rhs(F)
# ......................................... Solving
(A, b) = assemble_system(a, L, bc)
# ---------------------------Parallel Run
sol = LUSolver("mumps")
sol.solve(A, Result.vector(), b)
uf, Gu, LaMA, pot = Result.split()

# ......................................... Export Results

filename1 = "DIS_Beam.xdmf"
hdf5_file = XDMFFile(MPI.comm_world, filename1)
hdf5_file.write_checkpoint(uf, "F", XDMFFile.Encoding.HDF5)

filename2 = "pot_Beam.xdmf"
hdf5_file = XDMFFile(MPI.comm_world, filename2)
hdf5_file.write_checkpoint(pot, "F", XDMFFile.Encoding.HDF5)
``````

And here is the code in FEnicsx:

``````from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import pyvista

from dolfinx.fem import (Constant, Function, FunctionSpace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical,
locate_dofs_topological, assemble, locate_dofs_topological)
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc, LinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_unit_square, locate_entities, meshtags, locate_entities_boundary
from dolfinx import mesh as Mesh
from petsc4py.PETSc import ScalarType
import meshio
from dolfinx.plot import create_vtk_mesh
from ufl import (FacetNormal, FiniteElement, TensorElement, Identity, TestFunctions, TrialFunctions, VectorElement,
div, dot, dx, inner, lhs, nabla_grad, rhs, sym, indices, as_tensor, Measure, MixedElement, split)

# ......................................... Parameters
qe = 1.6021 * 1e-19
non_local = 1e-20
# -------------------------Paper
# E_1 = 1 * 152e9
# E_2 = 1 * 152e9
# E_3 = 1 * 152e9
E = 1 * 152e9
nu = 0.33
la = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
e0 = 8.854 * 1e-12
perm = 11.8 * e0  # 141.6e-9

f11 = 1 * 1e-6
f12 = 1 * 1e-6
f44 = 0.0

f11_1 = 1 * 1e-6
f11_3 = 1 * 1e-6
#
f12_1 = 1 * 1e-6
f12_3 = 1 * 1e-6
f44_1 = 0.0
f44_3 = 0.0

mu11 = 1e-6
mu12 = 1e-6
mu44 = 0.0

L = 1000e-9
H = 100e-9
h = H / 2

def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:, :2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read": [cell_data]})
return out_mesh

Vol_mesh = create_mesh(msh, "tetra", prune_z=False)
Surf_mesh = create_mesh(msh, "triangle", prune_z=False)
meshio.write("mesh.xdmf", Vol_mesh)
meshio.write("mt.xdmf", Surf_mesh)

with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)
with XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "r") as xdmf:

Force = 1e-1  # M / (0.5 * L)
Area = H ** 2
tr_force = Force / Area
tr_S1 = Constant(mesh, ScalarType((0, 0, 0.0)))
# tr_S2 = fe.Constant((0.0, -tr_force, 0.0))
tr_S2 = Constant(mesh, ScalarType((0, 0, -tr_force)))
n = FacetNormal(mesh)
# ......................................... Define Elements
v = VectorElement("CG", mesh.ufl_cell(), 2)  # Displacement Element
s = TensorElement("CG", mesh.ufl_cell(), 1)  # Relaxed strain Element
lm = TensorElement("CG", mesh.ufl_cell(), 1)  # Lagrange multipliers Element
p = FiniteElement("CG", mesh.ufl_cell(), 1)  # Electric Potential Element
# ......................................... Define mixed element
MIXEL = MixedElement([v, s, lm, p])
Space = FunctionSpace(mesh, MIXEL)
# ......................................... Define Trail and Test Functions
(u, gradu, lamu, p) = TrialFunctions(Space)
(del_u, del_gradu, del_lamu, del_p) = TestFunctions(Space)
# ....................................... Mark boundary regions
tol = 0.1 * 10e-9

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

def MidZ(x):
return np.isclose(x[2], H / 2, tol)

# # ....................................... ds
boundaries = [(1, lambda x: np.isclose(x[0], 0, tol)),
(2, lambda x: np.isclose(x[0], L, tol)),
(3, lambda x: np.isclose(x[2], H / 2, tol))]

facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)

# ....................................... BCs
S0, S0_to_Space = Space.sub(0).collapse()
# Clamp BC for left side
clamp = Function(S0)
clamp.x.set(0.0)
clamp_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, Left)
clamp_dofs = locate_dofs_topological((Space.sub(0), S0), mesh.topology.dim - 1, clamp_facets)

bc0 = dirichletbc(clamp, clamp_dofs, Space.sub(0))
bc = [bc0]

# ......................................... Define tensor
delta = Identity(len(u))
i, j, k, l, m, n = indices(6)

# Strain Tensor
def StrainT(u):
return as_tensor((1. / 2. * (u[i].dx(j) + u[j].dx(i))), [i, j])

# Stress Tensor
def Elas(la, mu):
return as_tensor(la * delta[i, j] * delta[k, l] + 2. * mu * delta[i, k] * delta[j, l],
(i, j, k, l))  # Elasticity tensor

def StressT(u):
return as_tensor((Elas(la, mu)[i, j, k, l] * StrainT(u)[k, l]), [i, j])

# Relaxed Strain Tensor
def RStrain(u):
return as_tensor((u[j].dx(i)), [i, j])

# Strain Gradient Tensor using relaxed strain
return as_tensor((1. / 2. * (gradu[k, j].dx(i) + gradu[j, k].dx(i))), [i, j, k])

# Higher Order Stress
def SGElas(la, non_local, mu):
return as_tensor((la * (non_local ** 2) * delta[i, j] * delta[l, m] * delta[k, n] +
2 * mu * (non_local ** 2) * delta[i, l] * delta[j, m] * delta[k, n]),
(i, j, k, l, m, n))  # Strain Gradient Elasticity Tensor

return as_tensor(
(SGElas(la, non_local, mu)[i, j, k, l, m, n] * RStrianGrad(gradu)[l, m, n]),
(i, j, k))

# Flexoelectricity Tensor
fl1 = np.zeros((3, 3, 3, 3))
fl1[0, 0, 0, 0] = 1
fl1[1, 1, 1, 1] = 1
fl1[2, 2, 2, 2] = 1

flexo1 = Constant(mesh, PETSc.ScalarType(fl1))
flexo = as_tensor(((mu11 - mu12 - 2 * mu44) * flexo1[i, j, k, l] + mu12 * delta[i, j] * delta[
k, l] + mu12 * delta[i, k] * delta[j, l] + mu12 * delta[i, l] * delta[j, k]), (i, j, k, l))

def HStressElec(p):
return as_tensor((-1 * flexo[l, i, j, k] * (-1 * nabla_grad(p))[l]), (i, j, k))

return as_tensor(perm * delta[i, j] * (-1 * nabla_grad(p))[j] + (flexo[i, j, k, l] * RStrianGrad(gradu)[j, k, l]),
i)

# ......................................... Dimensionless Form
# ......................................... Dimensionless Form
Fl = 1  # 10 ** (-6)
Fs = 1  # E
NONDIM = Constant(mesh, PETSc.ScalarType(1.0))
NONDIM2 = Constant(mesh, PETSc.ScalarType(1.0))
# ......................................... Define Variational Formulation
F1 = NONDIM * StressT(u)[j, k] * StrainT(del_u)[j, k] * dx - NONDIM * lamu[j, k] * del_u[k].dx(j) * dx - NONDIM * \
tr_S1[i] * del_u[i] * ds(1) - NONDIM * tr_S2[i] * del_u[i] * ds(2)
F2 = NONDIM2 * HStress(gradu)[i, j, k] * RStrianGrad(del_gradu)[i, j, k] * dx + NONDIM * lamu[j, k] * \
j, k] * dx

# F4 = 1e-20 * del_p * ds(1)  # Electrod BC
F4 = p.dx(i).dx(i) * ds(1)  # Electrod BC

# CONSTRAINT=psi[j,k]*del_lamu[j,k]*dx-u[k].dx(j)*del_lamu[j,k]
# j=0 and k=0
C00 = NONDIM * gradu[0, 0] * del_lamu[0, 0] * dx - NONDIM * u[0].dx(0) * del_lamu[0, 0] * dx
# j=0 and k=1
C01 = NONDIM * gradu[0, 1] * del_lamu[0, 1] * dx - NONDIM * u[1].dx(0) * del_lamu[0, 1] * dx
# j=0 and k=2
C02 = NONDIM * gradu[0, 2] * del_lamu[0, 2] * dx - NONDIM * u[2].dx(0) * del_lamu[0, 2] * dx
# j=1 and k=0
C10 = NONDIM * gradu[1, 0] * del_lamu[1, 0] * dx - NONDIM * u[0].dx(1) * del_lamu[1, 0] * dx
# j=1 and k=1
C11 = NONDIM * gradu[1, 1] * del_lamu[1, 1] * dx - NONDIM * u[1].dx(1) * del_lamu[1, 1] * dx
# j=1 and k=2
C12 = NONDIM * gradu[1, 2] * del_lamu[1, 2] * dx - NONDIM * u[2].dx(1) * del_lamu[1, 2] * dx
# j=2 and k=0
C20 = NONDIM * gradu[2, 0] * del_lamu[2, 0] * dx - NONDIM * u[0].dx(2) * del_lamu[2, 0] * dx
# j=2 and k=1
C21 = NONDIM * gradu[2, 1] * del_lamu[2, 1] * dx - NONDIM * u[1].dx(2) * del_lamu[2, 1] * dx
# j=2 and k=2
C22 = NONDIM * gradu[2, 2] * del_lamu[2, 2] * dx - NONDIM * u[2].dx(2) * del_lamu[2, 2] * dx

VarForm = F1 + F2 + F3 + C00 + C01 + C02 + C10 + C11 + C12 + C20 + C21 + C22 + F4

# ......................................... Preparing the form
F = VarForm
a = form(lhs(F))
L = form(rhs(F))

A = assemble_matrix(a, bcs=bc)
A.assemble()
b = create_vector(L)

# ......................................... create a linear solver

solver = PETSc.KSP().create(mesh.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

# ......................................... Solving
Result = Function(Space)
solver.solve(b, Result.vector)

# uf, Gu, LaMA, pot = Result.split()
uf = Result.sub(0).collapse()

# ......................................... Export Results
xdmf = XDMFFile(mesh.comm, "DIS_Beam.xdmf", "w")
xdmf.write_mesh(mesh)
xdmf.write_function(uf)

one = Constant(mesh, ScalarType(1.))
f = form(one * ds(2))
area = assemble_scalar(f)
print(area)

TR = Constant(mesh, ScalarType(tr_force))
f2 = form(TR * ds(2))
Traction = assemble_scalar(f2)
print(Traction)
``````

And the code I use to creat mesh is:

``````# -*- coding: utf-8 -*-
"""
Created on Fri Jun 17 11:13:55 2022

@author: Arash
"""

import gmsh
import sys

##################
gmsh.initialize()
# =====================================    Goccmetry      ======================
import gmsh
import sys

gmsh.initialize()

lc = 0

H = 100e-9
L = 1000e-9
W = 100e-9
h = W / 2

# H = 1e-2
# L = 5 * 1e-2
# W = 500e-6
# h = W / 2
# CREATING POINTS

# gmsh.model.geo.addPoint(0, 0, W - h, lc, 9)
# gmsh.model.geo.addPoint(L, 0, W - h, lc, 10)
# gmsh.model.geo.addPoint(L, H, W - h, lc, 11)
# gmsh.model.geo.addPoint(0, H, W - h, lc, 12)

# CREATING LINES & PLANES

# -----------------  Creat list
NX = 100
NY = 10
NZ = 3
type = "Progression "  # Bump";
PRD = 1
XEdgeList = [1, 5, 9, 3, 7, 11]
YEdgeList = [2, 6, 10, 4, 8, 12]
ZEdgeList = [3, 17, 16, 20, 15, 19, 14, 18]

gmsh.model.geo.synchronize()
gmsh.model.geo.addSurfaceLoop([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], 1)
gmsh.model.geo.synchronize()

gmsh.model.geo.synchronize()
gmsh.model.mesh.embed(2, [11], 3, 1)
gmsh.model.geo.synchronize()
# gmsh.model.geo.mesh.setTransfiniteSurface(11)

#############
###################

gmsh.model.geo.synchronize()
gmsh.option.setNumber('Mesh.MeshSizeMax', 3 * 10e-9)
gmsh.model.mesh.generate(3)
##############
gmsh.write("BeamPartitioned.msh")
######################
if '-nopopup' not in sys.argv:
gmsh.fltk.run()
#################################
gmsh.finalize()

``````

Your code is far from a minimal working example, and requires a lot of time for people to familiarize themself with.

What I would start with is to check that all boundary conditions are set correctly, then next up check that all subdomains (that being facets or cell tags) are marked correctly.

I would also check the convergence status of the ksp solver

``````solver.getConvergedReason()
``````
1 Like

Thank you so much and I appreciate your time. I worked on it and by changing the solver from

``````solver = PETSc.KSP().create(mesh.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
``````

TO

``````ksp = PETSc.KSP().create(mesh.comm)
ksp.setOperators(A)
ksp.setType(PETSc.KSP.Type.MINRES)
pc = ksp.getPC()
pc.setType("hypre")
pc.setHYPREType("boomeramg")
pc.setUp()
ksp.solve(b, Result.vector)
``````

Part of the problem is now solved. In this code , I am using a mixed formulation and I use Lagrange multiplier. Now, the solid mechanics part is working well but the results for the electrical domain is zero and I think the problem is with the way I set up the solver. Do you have any comments that might help me?

Thanks

Thank you so much @dokken for your comments. After checking my code, I solved the issue by modifying the solver as

``````problem = LinearProblem(a, L, bcs=bc,
petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
Result = problem.solve()
``````

However, I don’t know exactly what the problem was and why only this set up works for my problem. Generally, I would appreciate it if you could mention any source that can help me know about how to set up an efficient solver for a problem.

Thanks

See KSP: Linear System Solvers — PETSc v3.19.5-1107-g8cfeadfb465 documentation and each individual solver’s documentation.