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
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