UPDATED: Constructing nodal force vector involving shape function derivatives

1

[UPDATED POST BELOW WITH SHORTER MWE, HOPEFULLY CLEARER - LEAVING THIS ONE FOR CONTEXT]

For a shell problem, I need to implement a source term \mathbf{f} defined at the dofs of a mixed function space Z. In this MWE we have two unknowns: u and theta.
My source term is the componentwise summation of a global vector expression, evaluated at several points in the domain: [p_1,p_2,...p_n], and in the most general and simplified case, takes on the form:

\mathbf{f} = [...,f_{ui},f_{\theta i},...]^T =\Sigma_{p=1}^{n} [..., \frac{\partial N_i}{\partial \xi}\frac{\partial \xi}{\partial x} , N_i,...] . The reasons for doing this are complicated to explain and for this MWE this is sufficient. Here i is the ith node and N and \xi refer to the element shape function and parent coordinate, respectively. Following @kamensky 's answer to this post: I have the following MWE. However, I am not sure I handle the global dof numbering correctly because I am not getting the results I expected when I apply this to my problem. Also, I would really appreciate a cleaner way of implementing this if it is possible. I’m still using dolfin 2019 because I need doflin-adjoint. thank you in advance for your help

from dolfin import *
import numpy as np
from ufl import JacobianInverse

mesh = UnitSquareMesh(1,1)
J_inv = project(JacobianInverse(mesh))

tree = mesh.bounding_box_tree()
P2 = FiniteElement("P", triangle, degree=2)
bubble = FiniteElement("B", triangle, degree=3)
Z = FunctionSpace(mesh, MixedElement([P2 + bubble] + [P2]))
z = Function(Z)
u, theta = split(z)

# Evaluate the i-th basis function at point x:
def basis_func(i,x):
    cell_index = tree.compute_first_entity_collision(Point(*x))
    cell_global_dofs = Z.dofmap().cell_dofs(cell_index)
    for local_dof in range(0, len(cell_global_dofs)):
        if(i==cell_global_dofs[local_dof]):
            cell = Cell(mesh, cell_index)
            values = np.array([0,])
            return Z.element().evaluate_basis(local_dof, x,
                                             cell.get_vertex_coordinates(),
                                             cell.orientation())[0]

    # If none of this cell's shape functions map to the i-th basis function,
    # then the i-th basis function is zero at x.
    return 0.0

# Evaluate the i-th basis function n-th derivative in direction j at point x:
def basis_func_derivs(i, j, x, n=1):
    cell_index = tree.compute_first_entity_collision(Point(*x))
    cell_global_dofs = Z.dofmap().cell_dofs(cell_index)
    for local_dof in range(0, len(cell_global_dofs)):
        if(i==cell_global_dofs[local_dof]):
            cell = Cell(mesh, cell_index)
            values = np.array([0,])

            return Z.element().evaluate_basis_derivatives(local_dof, n, x,
                                              cell.get_vertex_coordinates(),
                                              cell.orientation())[j]

    # If none of this cell's shape functions map to the i-th basis function,
    # then the i-th basis function is zero at x.
    return 0.0

# Build a Python list of all basis functions:
list_of_basis_funcs = []
list_of_basis_funcs_derivs = []

for i in range(0,Z.dim()):
    list_of_basis_funcs += [(lambda i: lambda x : basis_func(i,x))(i),]

for i in range(0,Z.dim()):
    list_of_basis_funcs_derivs += [(lambda i: lambda x : basis_func_derivs(i,0,x,n=1))(i),]

def f_u(i,x):
    # the term corresponding to u
    return list_of_basis_funcs_derivs[i](x)*J_inv(x).reshape(2,2)[0,0]

def f_theta(i,x):
    # term corresponding to theta
    return list_of_basis_funcs[i](x)

# Some points for the source term:
pts = [[0.1, 0.75], [0.78, 0.18]]
    
f_vec = np.zeros(Z.dim())
for pt in pts:
    f_vec += np.array([f_u(i,pt) for i in range(Z.dim())])
    f_vec += np.array([f_theta(i,pt) for i in range(Z.dim())])
    
f=Function(Z)
f.vector()[:] = f_vec
2

[POST UPDATED FOR CLARITY, SHORTER MWE]

I construct a vector of nodal forces, which are functions of the coordinates of some points, given as a list. If a cell contains one of these points, it is to be assigned nodal forces. To simplify this post, without a loss of generality, these nodal dof forces are the element shape function derivatives (“strains”), evaluated at the point’s coordinates. Example code follows. I am looking for a simpler way to do this, preferably at a higher level using ufl directly (the ideal solution would play nice with dolfin-adjoint “out of the box”).

(For more context: Low level basis function operations in UFL)

from dolfin import *
import numpy as np
from ufl import JacobianInverse

mesh = UnitSquareMesh(1,1)
tree = mesh.bounding_box_tree()
V = FunctionSpace(mesh, 'CG', 1)

# Evaluate the i-th basis function n-th derivative in direction j at point x:
def basis_func_derivs(i, j, x, n=1):
    cell_index = tree.compute_first_entity_collision(Point(*x))
    cell_global_dofs = V.dofmap().cell_dofs(cell_index)
    for local_dof in range(len(cell_global_dofs)):
        if(i==cell_global_dofs[local_dof]):
            cell = Cell(mesh, cell_index)
            values = np.array([0,])

            return V.element().evaluate_basis_derivatives(local_dof, n, x,
                                              cell.get_vertex_coordinates(),
                                              cell.orientation())[j]

    # If none of this cell's shape functions map to the i-th basis function,
    # then the i-th basis function is zero at x.
    return 0.0

# Build a Python list of basis function derivatives dN/dxi:
list_of_basis_funcs_derivs = []
for i in range(0,V.dim()):
    list_of_basis_funcs_derivs += [(lambda i: lambda x : basis_func_derivs(i,0,x,n=1))(i),]

J_inv = project(JacobianInverse(mesh))

def f_i(i,x):      # the nodal force term
    return list_of_basis_funcs_derivs[i](x)*J_inv(x).reshape(2,2)[0,0]

# Some points to evaluate the source term:
pts = [[1,0],[0,1]]
f_vec = np.zeros(V.dim())

for pt in pts:
    f_vec += np.array([f_i(i,pt) for i in range(V.dim())])

print(f_vec)
f = Function(V)
f.vector()[:] = f_vec
3

Not sure why I didn’t think of this at first, but the terms involving the basis functions can be trivially obtained using PointSource. I’m posting an updated code with a test below. Also, I believe I was mistaken about including the inverse jacobian term for the derivatives as it seems from looking at the source code the derivatives are already transformed. I’d be very grateful if someone could offer suggestions for how to simplify obtaining the source terms involving the derivatives.

from dolfin import *
import numpy as np

mesh = UnitSquareMesh(1,1)
tree = mesh.bounding_box_tree()
P2 = FiniteElement("P", triangle, degree=2)
bubble = FiniteElement("B", triangle, degree=3)
Z = FunctionSpace(mesh, MixedElement([P2 + bubble] + [P2]))
subspace = 0  # or 1, for testing

# Evaluate the i-th basis function at point x:
def basis_func(i,x,subspace=0):
    cell_index = tree.compute_first_entity_collision(Point(*x))
    cell_global_dofs = Z.sub(subspace).dofmap().cell_dofs(cell_index)
    for local_dof in range(0, len(cell_global_dofs)):
        if(i==cell_global_dofs[local_dof]):
            cell = Cell(mesh, cell_index)
            values = np.array([0,])
            return Z.sub(subspace).element().evaluate_basis(local_dof, x,
                                             cell.get_vertex_coordinates(),
                                             cell.orientation())[0]
    # If none of this cell's shape functions map to the i-th basis function,
    # then the i-th basis function is zero at x.
    return 0.0

# Evaluate the i-th basis function n-th derivative in direction j at point x:
def basis_func_derivs(i, j, x, n=1):
    cell_index = tree.compute_first_entity_collision(Point(*x))
    cell_global_dofs = Z.dofmap().cell_dofs(cell_index)
    for local_dof in range(0, len(cell_global_dofs)):
        if(i==cell_global_dofs[local_dof]):
            cell = Cell(mesh, cell_index)
            values = np.array([0,])
            return Z.element().evaluate_basis_derivatives(local_dof, n, x,
                                              cell.get_vertex_coordinates(),
                                              cell.orientation())[j]
    # If none of this cell's shape functions map to the i-th basis function,
    # then the i-th basis function is zero at x.
    return 0.0

# Build a Python list of basis functions and derivatives:
N_list = []
dNdx_list = []

for i in range(0,Z.dim()):
    N_list += [(lambda i: lambda x : basis_func(i,x, subspace=subspace))(i),]

for i in range(0,Z.dim()):
    dNdx_list += [(lambda i: lambda x : basis_func_derivs(i,0,x,n=1))(i),]

def f_u(i,x):  # term containing basis function derivatives
    return dNdx_list[i](x)

def f_theta(i,x):  # term containing basis functions
    return N_list[i](x)

# Some points for the source term for testing:
pts = [[0.1, 0.75], [0.78, 0.18]]

f_vec_u = np.zeros(Z.dim())
f_vec_theta = np.zeros(Z.dim())
for pt in pts:
    f_vec_u += np.array([f_u(i,pt) for i in range(Z.dim())])
    f_vec_theta += np.array([f_theta(i,pt) for i in range(Z.dim())])

ftheta = Function(Z)
ftheta.vector()[:] = f_vec_theta

#### now replicate the functionality of f_theta using PointSource (works)
p = Function(Z)
for pt in pts: 
    pf = PointSource(Z.sub(subspace), Point(pt), 1)
    pf.apply(p.vector())
np.testing.assert_almost_equal(f_vec_theta, p.vector()[:])

# any way to do the same with dNdx?
1 Like