How to apply a force load at a known node number

General
9

Hi Jorgen,

I was able to write the code using the point source example and installed scifem but when I run the code I get that petsc4py.typing is not installed, even when the line “from petsc4py import PETSc” works.

Traceback (most recent call last):

File ~/.local/lib/python3.10/site-packages/spyder_kernels/customize/utils.py:209 in exec_encapsulate_locals
exec_fun(compile(code_ast, filename, “exec”), globals)

File ~/files/fenics_point_loading_example_2.py:19
import scifem

File ~/.local/lib/python3.10/site-packages/scifem/init.py:11
from .solvers import NewtonSolver

File ~/.local/lib/python3.10/site-packages/scifem/solvers.py:11
import petsc4py.typing

ModuleNotFoundError: No module named ‘petsc4py.typing’

The code is the following:

from petsc4py import PETSc
from mpi4py import MPI
import ufl
from dolfinx import mesh, fem, default_scalar_type, io, la
from dolfinx.fem import Function, form, petsc
from dolfinx.fem.petsc import assemble_matrix, apply_lifting
import numpy as np

import scifem



mu = 1
lambda_ = 1.25


domain = mesh.create_unit_cube(MPI.COMM_WORLD, nx=1, ny=1, nz=1, cell_type=mesh.CellType.tetrahedron)
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))

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

fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)

u_D = np.array([0, 0, 0], dtype=default_scalar_type)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)

T = fem.Constant(domain, default_scalar_type((0, 0, 0)))

ds = ufl.Measure("ds", domain=domain)


def epsilon(u):
    return ufl.sym(ufl.grad(u))  # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)

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


u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx

L = ufl.dot(T, v) * ds






b = Function(V)
b.x.array[:] = 0


geom_dtype = domain.geometry.x.dtype
if domain.comm.rank == 0:
    points = np.array([[1.0, 1.0, 1.0]], dtype=geom_dtype)
else:
    points = np.zeros((0, 3), dtype=geom_dtype)



point_source = scifem.PointSource(V, points, magnitude=np.array([[1.0, 1.0, 1.0]], dtype=geom_dtype))
point_source.apply_to_vector(b)


a_compiled = form(a)
apply_lifting(b.vector, [a_compiled], [[bc]])
b.x.scatter_reverse(la.InsertMode.add)
petsc.set_bc(b.vector, [bc])
b.x.scatter_forward()


A = assemble_matrix(a_compiled, bcs=[bc])
A.assemble()


ksp = PETSc.KSP().create(domain.comm)
ksp.setOperators(A)
ksp.setType(PETSc.KSP.Type.PREONLY)
ksp.getPC().setType(PETSc.PC.Type.LU)
ksp.getPC().setFactorSolverType("mumps")


# Compute solution
uh = Function(V)
ksp.solve(b.vector, uh.vector)
uh.x.scatter_forward()




with io.XDMFFile(domain.comm, "deformation.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    uh.name = "Deformation"
    xdmf.write_function(uh)