How to create a source term that is a spatial function?

variational formulation
1

This is the original code I took from The FEniCSx tutorial:

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

I want to change the constant -rho*g in the source term f to a function of the spatial coordinates. How do I do that?

2

See for instance Implementation — FEniCSx tutorial

3

I can’t believe I missed that.
Thank you very much!

4

I tried to substitude the constant term in the

with a spatial function and it gives the error:

ValueError: setting an array element with a sequence.

and when I tried to using this spatial function directly in the linear elasticity problem like the following:

L = ufl.dot(p,v) * ds(7) #+ ufl.dot(T, v) * ds

it gives the error:

UFLException: Dot product requires non-scalar arguments, got arguments with ranks 0 and 1.

Can somebody please tell me what is wrong with the code?

EDIT: I tried writing source term as a vector and insert into the linear elasticity problem. And it seems worked: (p is a spatial function) I just don’t know if it is the right way to do it.

f = ufl.as_vector([0,0,-p])
5

You would have to make a reproducible example for anyone to help you.

6

Thank you for the reply, dokken. Here is my MWE:

from dolfinx import fem,mesh
import ufl
from mpi4py import MPI
import numpy as np


N = 30
domain = mesh.create_box(MPI.COMM_WORLD,[[0.0,0.0,0.0],[1.0,1.0,-1.0]],[N,N,N//2],mesh.CellType.hexahedron)
x = ufl.SpatialCoordinate(domain)
d = 0.02
R = 0.5
a = np.sqrt(R*d)

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

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

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

def bottom(x):
    return np.isclose(x[2],-1) 

def bound_xy(x):
    return np.logical_and(np.logical_and(np.greater(x[0],0) , np.less(x[0],a)),np.logical_and(np.greater(x[1],0), np.less(x[1],a)))

fdim = domain.topology.dim - 1
symx_facets = mesh.locate_entities_boundary(domain, fdim, symmetry_x)
symy_facets = mesh.locate_entities_boundary(domain, fdim, symmetry_y)
bottom_facets = mesh.locate_entities_boundary(domain, fdim, bottom)
top_facets = mesh.locate_entities_boundary(domain, fdim, top)
xy_facets = mesh.locate_entities_boundary(domain, fdim, bound_xy)

marked_facets = np.hstack([symx_facets,symy_facets,bottom_facets,top_facets,x_facets,y_facets,xy_facets])
marked_values = np.hstack([np.full_like(symx_facets,1), np.full_like(symy_facets,2), np.full_like(bottom_facets,3),np.full_like(top_facets,4),np.full_like(xy_facets,5)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets],marked_values[sorted_facets])

V = fem.VectorFunctionSpace(domain,("CG",1))

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# metrics
E = 10
nu = 0.3
lmbda = E*nu/((1+nu)*(1-2*nu))
mu = E/(2*(1+nu))

# load
F = 4/3.*E/(1-nu**2)*a*d
p0 = 3*F/(2*math.pi*a**2)
p = p0*(1-(x[0]**2+x[1]**2)/a**2)**0.5
f = ufl.as_vector([0,0,-p])

sub_id = 5
ds = ufl.Measure("ds",domain=domain,subdomain_data=facet_tag,metadata = metadata, subdomain_id=sub_id)

# BCs
symx_dofs = fem.locate_dofs_topological(V.sub(0),fdim, facet_tag.find(1))
symy_dofs = fem.locate_dofs_topological(V.sub(1),fdim, facet_tag.find(2))
bottom_dofs = fem.locate_dofs_topological(V,fdim, facet_tag.find(3))

u0 = np.array([0.0,0.0,0.0], dtype=ScalarType)
bc_bottom = fem.dirichletbc(u0,bottom_dofs,V)
bc_x = fem.dirichletbc(ScalarType(0),symx_dofs,V.sub(0))
bc_y = fem.dirichletbc(ScalarType(0),symy_dofs,V.sub(1))

bcs = [bc_bottom,bc_x,bc_y]

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 lmbda * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2*mu*epsilon(u)

a_ = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f,v) * ds #+ ufl.dot(T, v) * ds
problem = fem.petsc.LinearProblem(a_, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

I Want to implement this code as the load control variant of the Hertzian contact problem. Which has a load p(r) as a spatial function in the contact area within the a

7

What you are doing with

Is correct. This adds a spatially varying load in the z direction.

2 Likes