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

dolfinx
1

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.

image
image827×427 3.16 KB

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)

# ......................................... Reading mesh file
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):
    return sym(nabla_grad(u))  # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)


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
2

See for instance Component Wise Dirichlet BC - #2 by dokken or Component-wise Dirichlet BC — FEniCSx tutorial

3

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))
# ......................................... Reading mesh file
mesh = fe.Mesh("Beam.xml")
Thickness = 20e-9
L = 1000e-9
H = 100e-9
# ......................................... Define Loading
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
def RStrianGrad(gradu):
    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


def HStress(gradu):
    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))


def ElecDis(gradu):
    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] * \
     del_gradu[
         j, k] * dx

F3 = HStressElec(p)[i, j, k] * RStrianGrad(del_gradu)[i, j, k] * dx - ElecDis(gradu)[i] * (-1 * grad(del_p)[i]) * 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

# ......................................... Reading mesh file
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



msh = meshio.read("BeamPartitioned.msh")

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 = xdmf.read_mesh(name="Grid")
    ct = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)
with XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "r") as xdmf:
    ft = xdmf.read_meshtags(mesh, name="Grid")

# ......................................... Define Loading

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
def RStrianGrad(gradu):
    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


def HStress(gradu):
    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))


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


def Polarization(gradu):
    return as_tensor((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] * \
     del_gradu[
         j, k] * dx

F3 = HStressElec(p)[i, j, k] * RStrianGrad(del_gradu)[i, j, k] * dx - ElecDis(gradu)[i] * (
        -1 * nabla_grad(del_p)[i]) * 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()
gmsh.model.add("Beam")
# =====================================    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, 0, lc, 1)
gmsh.model.geo.addPoint(L, 0, 0, lc, 2)
gmsh.model.geo.addPoint(L, H, 0, lc, 3)
gmsh.model.geo.addPoint(0, H, 0, lc, 4)

gmsh.model.geo.addPoint(0, 0, h, lc, 5)
gmsh.model.geo.addPoint(L, 0, h, lc, 6)
gmsh.model.geo.addPoint(L, H, h, lc, 7)
gmsh.model.geo.addPoint(0, H, h, lc, 8)

# 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)

gmsh.model.geo.addPoint(0, 0, W, lc, 9)
gmsh.model.geo.addPoint(L, 0, W, lc, 10)
gmsh.model.geo.addPoint(L, H, W, lc, 11)
gmsh.model.geo.addPoint(0, H, W, lc, 12)
# CREATING LINES & PLANES
gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(2, 3, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(4, 1, 4)
gmsh.model.geo.addCurveLoop([1, 2, 3, 4], 1)
gmsh.model.geo.addPlaneSurface([1], 1)

gmsh.model.geo.addLine(5, 6, 5)
gmsh.model.geo.addLine(6, 7, 6)
gmsh.model.geo.addLine(7, 8, 7)
gmsh.model.geo.addLine(8, 5, 8)

gmsh.model.geo.addLine(9, 10, 9)
gmsh.model.geo.addLine(10, 11, 10)
gmsh.model.geo.addLine(11, 12, 11)
gmsh.model.geo.addLine(12, 9, 12)

gmsh.model.geo.addCurveLoop([9, 10, 11, 12], 2)
gmsh.model.geo.addPlaneSurface([2], 2)

gmsh.model.geo.addLine(1, 5, 13)
gmsh.model.geo.addLine(2, 6, 14)
gmsh.model.geo.addCurveLoop([1, 14, -5, -13], 3)
gmsh.model.geo.addPlaneSurface([3], 3)

gmsh.model.geo.addLine(3, 7, 15)
gmsh.model.geo.addCurveLoop([2, 15, -6, -14], 4)
gmsh.model.geo.addPlaneSurface([4], 4)

gmsh.model.geo.addLine(4, 8, 16)
gmsh.model.geo.addCurveLoop([13, -8, -16, 4], 5)
gmsh.model.geo.addPlaneSurface([5], 5)

gmsh.model.geo.addCurveLoop([16, -7, -15, 3], 6)
gmsh.model.geo.addPlaneSurface([6], 6)

gmsh.model.geo.addLine(5, 9, 17)
gmsh.model.geo.addLine(6, 10, 18)
gmsh.model.geo.addCurveLoop([5, 18, -9, -17], 7)
gmsh.model.geo.addPlaneSurface([7], 7)

gmsh.model.geo.addLine(7, 11, 19)
gmsh.model.geo.addLine(8, 12, 20)
gmsh.model.geo.addCurveLoop([-7, 19, 11, -20], 8)
gmsh.model.geo.addPlaneSurface([8], 8)

gmsh.model.geo.addCurveLoop([6, 19, -10, -18], 9)
gmsh.model.geo.addPlaneSurface([9], 9)

gmsh.model.geo.addCurveLoop([8, 17, -12, -20], 10)
gmsh.model.geo.addPlaneSurface([10], 10)

# -----------------  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.addVolume([1], 1)
gmsh.model.geo.synchronize()
gmsh.model.geo.addCurveLoop([5, 6, 7, 8], 11)
gmsh.model.geo.addPlaneSurface([11], 11)

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

#############
gmsh.model.addPhysicalGroup(2,[5,10], name="Left")
gmsh.model.addPhysicalGroup(2,[4,9], name="Right")
gmsh.model.addPhysicalGroup(2,[2], name="Top")
gmsh.model.addPhysicalGroup(2,[11], name="Mid")
gmsh.model.addPhysicalGroup(3,[1], name="Beam")
###################

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()
4

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
5

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

6

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

7

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