Multiplying scalar and vector functions on components of mixed elements

General
1

Hello everyone!

I am working up to implementing some nonlinear beam modeling in FEniCSx, but as a stepping stone to debug some of the implementation issues, I am playing with a 1D catenary problem. For a 1D cable in 2D space, parameterized by its arclength s, the position vector is \bar{r}(s)=(x(s),y(s)), and the governing equations become

\frac{\partial}{\partial s} \left( N(s) \frac{\bar{t}(s)}{|\bar{t}(s)|} \right) + \bar{f}(s) = 0
N(s) = k(|\bar{t}(s)| - 1)

where N(s) is the scalar tension in the cable, \bar{t}(s) = \nabla_s\bar{r} is the tangent vector, \bar{f} is a distributed force, and k is a stiffness. For this problem, I am tracking both the position (\bar{r}) and the tension (N) as state variables (while not strictly necessary here, this will make my life much easier in future problems).

In this problem (and the problems I will be moving to in the future), the governing PDEs include both scalar and vector components that are multiplied together. I have seen some other examples on the discussion boards, but not which deal with functions that are part of the same mixed element formulation. When I attempt to implement a version of the suggestions from other posts, I hit a wall.

I have included the code I have for this problem up to the formation of the weak form. When I run this code, I get an error with the “normal_force_vector” function saying Invalid ranks 1 and 2 in product. (I have also attempted to apply the multiplication to the vector of node values, but I get errors saying that “Grad” doesn’t have a node vector like the Functions do). Any suggestions on how to move forward with problems involving scalar-vector multiplication like this would be greatly appreciated!

# static 1D hung cable (caternary)
import numpy as np
from mpi4py import MPI
import basix
from dolfinx.fem import (Function, dirichletbc,functionspace,
                         locate_dofs_topological,Constant,Expression,assemble_vector)
from dolfinx.io import VTKFile
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import locate_entities_boundary,create_interval
from dolfinx import default_scalar_type, geometry
from ufl import (SpatialCoordinate,inner, TestFunction, TrialFunction, div, grad, dx, dot,split, sqrt)
import matplotlib.pyplot as plt


## Constants for the problem
L = 1.0 # length of the cable
k = 1.0 # stiffness of the cable
fy = -0.1 # force per unit length in y direction

## Create the domain over the arclength of the cable
NUM_ELEM = 10
domain = create_interval(MPI.COMM_WORLD, NUM_ELEM, [0, L])

S = SpatialCoordinate(domain) # array of arclength parameter

# convert model constants to the domain
L = Constant(domain,L)
k = Constant(domain,k)
f = Constant(domain,[0,fy])

## Create the elements -- using Hermite here for compatibility with beam/plate equations
# Create Hermite order 3 on a interval (for more informations see:
#    https://defelement.com/elements/examples/interval-Hermite-3.html )
beam_vector_element=basix.ufl.element(  basix.ElementFamily.Hermite,     # family of element (Hermite polynomial so is C^1)
                                        basix.CellType.interval,         # cell type that the element is defined over (here a line)
                                        3,                               # degree of the element (3 for C^1 Hermite)
                                        shape=(2,))                      # create a vector-valued element (and it derivative bc Hermite)

beam_scalar_element=basix.ufl.element(  basix.ElementFamily.Hermite,     # family of element (Hermite polynomial so is C^1)
                                        basix.CellType.interval,         # cell type that the element is defined over (here a line)
                                        3,                               # degree of the element (3 for C^1 Hermite)
                                        shape=(1,))                      # create a scalar-valued element (and it derivative bc Hermite)

## We will need 1 vectors (r: position) and 1 scalar (T: tension)
GEB_beam_element=basix.ufl.mixed_element([  beam_vector_element,    # position
                                            beam_scalar_element])  # tension

V = functionspace(domain,GEB_beam_element) # create a vector function space over the domain
V_vec = functionspace(domain,beam_vector_element)
V_sca = functionspace(domain,beam_scalar_element)

### defining the functions for the weak form
state = Function(V, name='state')
r,T = split(state)  ## splitting the state into the vector (position) and scalar (tension) component

d_state = TestFunction(V)
dr,dT = split(d_state) ## splitting the test function


# initial condition
state.x.array[:] = 0
state.sub(0).interpolate(lambda x: np.vstack((x[0], 0*x[0])) ) # interpolating the initial position of the beam 



## defining all of the relevant functions
t = grad(r) ## get the tangent vector

def eps(t):
    # calculate the strain
    return sqrt(inner(t,t)) - 1

def normal_force_vector(T,t):
    # returns the normal force vector (scalar tension times unit tangent vector)
    T_force = Function(V_vec)

    T_force.interpolate(Expression( T * t / sqrt(inner(t,t))  ,V_vec.element.interpolation_points() ))
    return T_force



## Weak form Equation 1: T(s)t(s) * dv/ds * dx + f * v * dx = 0

F1 = inner( normal_force_vector(T,t), grad(dr) )*dx + inner(f,dr)*dx

## Weak form Equation 2: T(s) * v(s) * dx - k eps(s) * v(s) * dx = 0

F2 = inner( T, dT )*dx - k*inner( eps(t), dT )*dx

## combined problem 
F = F1 + F2
2

You have:

(Pdb++) T.ufl_shape
(1,)
(Pdb++) t.ufl_shape
(2, 1)

So maybe you mean computing t*T ?
additionally, since you are in 1-D, maybe you want to use dr.dx(0) rather than grad(dr), or use grad(dr)[0]?

You can easily inspect the shape of the various objects, as I showed above, and think about what makes mathematical sense.

3

Thanks so much! The .ufl_shape tool was very useful. Using this tool, I found that the bigger issue was that I was making a 1-by-1 vector instead of a scalar when I was creating the Hermite elements. I fixed this by changing

beam_scalar_element=basix.ufl.element(  basix.ElementFamily.Hermite,     # family of element (Hermite polynomial so is C^1)
                                        basix.CellType.interval,         # cell type that the element is defined over (here a line)
                                        3,                               # degree of the element (3 for C^1 Hermite)
                                        shape=(1,))                      # create a scalar-valued element (and it derivative bc Hermite)

to

beam_scalar_element=basix.ufl.element(  basix.ElementFamily.Hermite,     # family of element (Hermite polynomial so is C^1)
                                        basix.CellType.interval,         # cell type that the element is defined over (here a line)
                                        3,                               # degree of the element (3 for C^1 Hermite)
                                        )                      # create a scalar-valued element (and it derivative bc Hermite)

(Removing the shape argument). This actually allowed either T*t and t*T to produce the desired result.

Additionally, thank you for pointing out the .dx(0) function – I was not aware of it and it is definitely the better tool to use here.