Using position-dependent body forces

Hi!

I read many FEniCS examples, where body forces used in the variational formulation are constants.
e.g. in the linear elasticity solver Solving PDEs in Python - <br> The FEniCS Tutorial Volume I with weight force : f = Constant((0, 0, -rho*g)).

I would like to test body forces of which the expression depends on the location in the mesh.
Here are my questions:

  • Would you have some examples of implementation of position-dependent body forces?
  • It seems that body forces are usually expressed with the Cartesian coordinates. Is there any FEniCS function that enable to express easily the position-dependence with spherical coordinates? Would you also have some examples?

Many thanks,
Anne

DOLFIN (and DOLFINx) has the function

x = SpatialCoordinate(mesh)
f = as_vector((x[0], x[1], cos(x[2]))

It has for instance been done for cylindrical coordinates by @CastriMik in Scattering from a sphere in the axisymmetric formulation — GSoC 2022 - Expanding FEniCSx electromagnetic demos

or
https://fenicsproject.org/qa/10476/the-weak-formulation-of-spherical-coordinates/
or

Thank you for the lightening

Just one precision, does ‘x = SpatialCoordinate(mesh)’ enable to compute spherical coordinates? Many thanks.

The mesh is agnostic to the coordinate system, although it typically assumes an orthogonal spatial basis. Your implementation of integral/differential operators can be used to define the coordinate system.

see e.g. Del in cylindrical and spherical coordinates - Wikipedia

See also the excellent (although in old FEniCS) example here: Axisymmetric formulation for elastic structures of revolution — Numerical tours of continuum mechanics using FEniCS master documentation.

Thank you nate.

Actually, in the FEniCS example, I don’t really get how the solver can “understand” that in the bilinear form “a = inner(sigma(du), eps(u_))*x[0]*dx” for instance, x[0] refers to r. Or that x[1] refers to z. Could you bring me some further explanation?

In my experiments, the weak form I use does not specify the coordinates system. But I would like to express the body force in the spherical basis.
So, I first started to use spherical coordinates as follow:
x = SpatialCoordinate(mesh)
radius = safe_sqrt( x[0]*x[0] + x[1]*x[1] + x[2]*x[2])
theta = atan((safe_sqrt( x[0]*x[0] + x[1]*x[1])/x[2])
phi = atan(x[1]/x[0])
But then, I used instead:
r, phi, theta = SpatialCoordinate(mesh)
Which formulation should I prefer?

Thank you for the help.

Hi!

Actually, I would like to compute for each node of the mesh the spherical coordinates and associated local spherical basis. As I am in dynamic case, I would like them to be updated at each solver iteration.

Here are my questions:

  • I understood, but maybe I am wrong, that I needed to define the nodal spherical coordinates and basis as finite element functions to be sure they would be updated. Is it correct? Is there another way?
  • In the below code, I cannot associate a ‘fe.as_vector’ (dimension 3) to ‘e_r_vect’ (np.array([0.,0.,0…]) ) and get a ‘setting an array element with a sequence’ ValueError. As e_r_vect is supposed to contain the 3d components for all node (I thought shape would be : (number of node, 3) ), how is it possible to access them precisely from the shape(3 * number of nodes,)-structure of ‘e_r_vect’?

For now, I tried the following code:

Define vector space and associated spherical basis vectors functions

VectorSpace_DG1 = VectorFunctionSpace(mesh, “Lagrange”, 1)
e_r = Function(VectorSpace_DG1)
e_theta = Function(VectorSpace_DG1)
e_phi = Function(VectorSpace_DG1)

Define scalar space and associated functions

ScalarSpace_DG1 = FunctionSpace(mesh, “Lagrange”, 1)
r = Function(ScalarSpace_DG1)
theta = Function(ScalarSpace_DG1)
phi = Function(ScalarSpace_DG1)

coordinates_all_nodes = np.array([ np.array([node.point().x(), node.point().y(), node.point().z()]) for dof, node in enumerate(fe.vertices(mesh))])

Dynamic spherical coordinates?

r_vec = r.vector()
theta_vec = theta.vector()
phi_vec = phi.vector()

for dof, node in enumerate(vertices(fenics_mesh)):

r_vec[dof] = sqrt( coordinates_all_nodes[dof][0]*coordinates_all_nodes[dof][0] +
coordinates_all_nodes[dof][1]*coordinates_all_nodes[dof][1] +
coordinates_all_nodes[dof][2]*coordinates_all_nodes[dof][2] )
theta_vec[dof] = acos(coordinates_all_nodes[dof][2]/r_vec[dof])
phi_vec[dof] = atan(coordinates_all_nodes[dof][1]/coordinates_all_nodes[dof][0])

r.vector()[:] = r_vec
theta.vector()[:] = theta_vec
phi.vector()[:] = phi_vec

Dynamic spherical basis?

e_r_vect = e_r.vector()[:]
e_theta_vect = e_theta.vector()[:]
e_phi_vect = e_phi.vector()[:]

for dof, node in enumerate(vertices(fenics_mesh)):

e_r_vect[dof] = as_vector( (sin(theta_vec[dof])*cos(phi_vec[dof]), sin(theta_vec[dof])*sin(phi_vec[dof]), cos(theta_vec[dof])) ) → setting an array element with a sequence’ ValueError
e_theta_vect[dof] = as_vector( (cos(theta_vec[dof])*cos(phi_vec[dof]), cos(theta_vec[dof])*sin(phi_vec[dof]), sin(theta_vec[dof])) )
e_phi_vect[dof] = as_vector( (-sin(phi_vec[dof]), cos(phi_vec[dof]), 0.) )

e_r.vector()[:] = e_r_vect
e_theta.vector()[:] = e_theta_vect
e_phi.vector()[:] = e_phi_vect

Many thanks!