Thanks for your response. I refereed to the previous discussion on point source How to construct a sparse matrix used for calculating function values at multiple points - #4 by dokken which is very much helpful. As per the advice, I attempted to apply point load in a simple problem. I created function and obtained the basis values and cells from the function.
After that, it is needed to be inserted in the matrix/vector, which is where I am confused in execution of the same. I tried so far what I understood in the below code on a simple problem of rectangular domain with application of point load on top middle surface, but it is showing Error which is not usual. 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, io
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
from dolfinx import fem, default_scalar_type,mesh
import scipy.sparse as sps
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.fem import (Function, dirichletbc, form, functionspace, locate_dofs_geometrical)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
create_matrix, create_vector, set_bc)
from ufl import TestFunction, TrialFunction, derivative, dx, grad, inner
def construct_force_matrix(V, points):
mesh = V.mesh
sdim = V.dofmap.index_map_bs
num_dofs = V.dofmap.index_map.size_local * sdim
nx, dim = points.shape
_, _, owning_points, cells = dolfinx.cpp.geometry.determine_point_ownership(
mesh._cpp_object, points)
owning_points = np.asarray(owning_points).reshape(-1, 3)
mesh_nodes = mesh.geometry.x
cmap = mesh.geometry.cmaps[0]
ref_x = np.zeros((len(cells), mesh.geometry.dim), dtype=mesh_nodes.dtype)
for i, (point, cell) in enumerate(zip(owning_points, cells)):
geom_dofs = mesh.geometry.dofmap[cell]
ref_x[i] = cmap.pull_back(point.reshape(-1, 3), mesh_nodes[geom_dofs])
u = ufl.TrialFunction(V)
num_dofs = V.dofmap.dof_layout.num_dofs * V.dofmap.bs
if len(cells) > 0:
# Strip out basis function values per cell
expr = dolfinx.fem.Expression(u, ref_x, comm=MPI.COMM_SELF)
values = expr.eval(mesh, cells)
# Strip out basis function values per cell
basis_values = values[:num_dofs:num_dofs*len(cells)]
else:
basis_values = np.zeros(
(0, num_dofs), dtype=dolfinx.default_scalar_type)
print (cells)
print (basis_values)
return cells, basis_values
class NonlinearPDEProblem:
"""Nonlinear problem class for a PDE problem."""
def __init__(self, F, u, bc):
V = u.function_space
du = TrialFunction(V)
self.L = form(F)
self.a = form(derivative(F, u, du))
self.bc = bc
def form(self, x):
x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
def F(self, x, b):
"""Assemble residual vector."""
with b.localForm() as b_local:
b_local.set(0.0)
assemble_vector(b, self.L)
apply_lifting(b, [self.a], bcs=[[self.bc]], x0=[x], scale=-1.0)
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, [self.bc], x, -1.0)
#Constructing the force matrix for the point load
cells, basis_values = construct_force_matrix(x.function_space, points_array)
# Add the point load contribution to the residual vector
with b.localForm() as b_local:
for i in range(len(cells)):
x_vals = np.repeat(x.array[cells[i]], x.function_space.dofmap.index_map_bs)
x_local = np.sum(x_vals * basis_values[i, :])
b_local.array[cells[i]] -= point_load_magnitude * basis_values[i, :] * x_local
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, self.bcs, x, -1.0)
def J(self, x, A):
"""Assemble Jacobian matrix."""
A.zeroEntries()
assemble_matrix(A, self.a, bcs=[self.bc])
A.assemble()
def matrix(self):
return create_matrix(self.a)
def vector(self):
return create_vector(self.L)
#Creating mesh
msh = mesh.create_rectangle (comm=MPI.COMM_WORLD, points=((0.0, 0.0), (6.0,1.0)), n=(96,16), cell_type=mesh.CellType.quadrilateral)
#Creating function space
V = fem.VectorFunctionSpace(msh,("Lagrange",1))
def clamped_boundary(x):
return np.isclose(x[1], 0)
u_zero = np.array((0,) * msh.geometry.dim, dtype=default_scalar_type)
bc = fem.dirichletbc(u_zero, fem.locate_dofs_geometrical(V, clamped_boundary), V)
def left(x):
return np.isclose(x[0],0)
left_facets = mesh.locate_entities_boundary(msh,msh.topology.dim-1,left)
left_dof = fem.locate_dofs_topological(V.sub(0),msh.topology.dim-1,left_facets)
bc_left = fem.dirichletbc(default_scalar_type(0),left_dof,V.sub(0))
bcs = [bc,bc_left]
#Defining point for application of force
point_to_apply_force = np.array([3,1])
points_array = np.array([point_to_apply_force])
# Call the function with the array of points
cells,basis_val = construct_force_matrix(V, points_array)
#Defining traction and body forces
T=fem.Constant(msh,default_scalar_type((0,0)))
B=fem.Constant(msh,default_scalar_type((0,0)))
#Defining test and solution function
v=ufl.TestFunction(V)
u=fem.Function(V)
#Spatial dimension
d=len(u)
#Identity tensor
I=ufl.variable(ufl.Identity(d))
#Deformation gradient
F=ufl.variable(I+grad(u))
# 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 = 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 / 2) * (Ic - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2
# Hyper-elasticity
P = ufl.diff(psi, F)
ds = ufl.Measure('ds')
dx = ufl.Measure('dx')
#Defining the weak form
FF = ufl.inner(ufl.grad(v), P) * dx - ufl.inner(v, B) * dx - ufl.inner(v, T) * ds
# External point load magnitude
point_load_magnitude = 100.0
problem = NonlinearProblem(FF, u, bcs, J=None, form_compiler_options=None, jit_options=None)
#Solving using newton solver
solver=NewtonSolver(msh.comm,problem)
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"
solver.solve(u) #Error when tried to solve.
a=u.x.array
print(a)
The error is may be my proceeding with basis value is wrong on my understanding. If approach of solving is correct, it will be very much helpful for all questions related to application of point load. As per my understanding I started, could you please help me further how to avoid this error . It will be very much helpful to proceed further for my problem. I hope the issues is in the addition of vector to existing vector, please help in further proceeding.
Thanks for all your time and considerations. Any tips or advice on the above problem is greatly appreciated. I will be waiting for your responses.
Error message for your reference:
terminate called after throwing an instance of 'std::runtime_error'
what(): pybind11_object_dealloc(): Tried to deallocate unregistered instance!
[Legion:476503] *** Process received signal ***
[Legion:476503] Signal: Aborted (6)
[Legion:476503] Signal code: (-6)
[Legion:476503] [ 0] /lib/x86_64-linux-gnu/libc.so.6(+0x42520)[0x7fa2ba0d3520]
[Legion:476503] [ 1] /lib/x86_64-linux-gnu/libc.so.6(pthread_kill+0x12c)[0x7fa2ba1279fc]
[Legion:476503] [ 2] /lib/x86_64-linux-gnu/libc.so.6(raise+0x16)[0x7fa2ba0d3476]
................
[Legion:476503] [25] /lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0x80)[0x7fa2ba0bae40]
[Legion:476503] [26] /bin/python3(_start+0x25)[0x55a2d6110f25]
[Legion:476503] *** End of error message ***