Frequency Domain Example

1

Esteemed,

I am studying about frequency domain and found the following page describing an example of a frequency domain problem. But in the code available there are some points that I’m not able to solve, such as the term i,j,k,l=indices(4), where I’ve never seen this function for fenics.

I made some modifications so I don’t need to use the .xml file and updated the mixed functions, but I still can’t run the code.

Can this code still be “updated” or is there already a better way to solve the frequency domain problem?

Follow the website: http://mypages.iit.edu/~asriva13/?page_id=779
Link with theory: http://mypages.iit.edu/~asriva13/wp-content/uploads/2015/09/Frequency-Domain.pdf
Code:

from dolfin import *
import numpy
from ufl import indices

frequency=50000
omega=2*numpy.pi*frequency

#Mesh information

#mesh = Mesh("geom.xml");
#cd=MeshFunction('size_t',mesh,"geom_physical_region.xml");
#fd=MeshFunction('size_t',mesh,"geom_facet_region.xml");

L = 0.2
H = 0.01
Nx = 200
Ny = 10
mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny)

#Function space definition over the mesh
#V = VectorFunctionSpace(mesh, "CG", 1) ##Continuous Galerkin for the displacement field
#Vcomplex = V*V

Element1 = FiniteElement("CG", mesh.ufl_cell(), 1)
Element2 = FiniteElement("CG", mesh.ufl_cell(), 1)

# Define the mixed element
V_elem = MixedElement([Element1, Element2])

# Define the mixed function space
Vcomplex = FunctionSpace(mesh, V_elem)

V0=FunctionSpace(mesh, 'DG', 0) ##Discontinuous Galerkin for the material property fields
M0=TensorFunctionSpace(mesh, 'DG', 0, shape=(2,2,2,2))

#Material properties
E, rho, nu = 300E+8, 8000, 0.3

lam=E*nu/((1+nu)*(1-2*nu))
mu=E/(2*(1+nu))
i,j,k,l=indices(4)
delta=Identity(2)

C=as_tensor((lam*(delta[i,j]*delta[k,l])+mu*(delta[i,k]*delta[j,l]+delta[i,l]*delta[j,k])),(i,j,k,l))

#Applying the Boundary conditions
zero=Constant((0.0, 0.0))
one=Constant((1.0, 0.0))

def left(x, on_boundary):
    return near(x[0], 0.)

boundary = [DirichletBC(Vcomplex.sub(0), zero, left), DirichletBC(Vcomplex.sub(1), zero, left)]

#boundary = [DirichletBC(Vcomplex.sub(0), zero, fd, 12), DirichletBC(Vcomplex.sub(1), zero, fd, 12)]
#source = [DirichletBC(Vcomplex.sub(0), one, fd, 11), DirichletBC(Vcomplex.sub(1), one, fd, 11)]

#bc=boundary+source
bc=boundary

#Real and Imaginary parts of the trial and test functions
uR, uI = TrialFunctions(Vcomplex)
wR, wI = TestFunctions(Vcomplex)

strainR=as_tensor(0.5*(Dx(uR[i],j)+Dx(uR[j],i)),(i,j))
strainI=as_tensor(0.5*(Dx(uI[i],j)+Dx(uI[j],i)),(i,j))

stressR=as_tensor((C[i,j,k,l]*strainR[k,l]),[i,j])
stressI=as_tensor((C[i,j,k,l]*strainI[k,l]),[i,j])

F=omega**2*rho*(dot(uR,wR)+dot(uI,wI))*dx-(inner(stressR, grad(wR))+inner(stressI, grad(wI)))*dx
a, L = lhs(F), rhs(F)

# Set up the PDE
solnU = Function(Vcomplex)
problem = LinearVariationalProblem(a, L, solnU, bcs=bc)
solver = LinearVariationalSolver(problem)
solver.solve()
solnUR, solnUI=solnU.split()
plot(solnUI)
interactive()
2

indices is a UFL function that can be import by calling from ufl import indices.

In general, you could switch to DOLFINx, as you would not have to split your problem into a real and imaginary component as there is support for complex numbers.
See for instance: https://github.com/FEniCS/dolfinx/blob/main/python/demo/helmholtz2D/demo_helmholtz_2d.py

1 Like
3

Dear Dokken,

Thank you for your reply. But I’m still working with Dolfin-Adjoint, which doesn’t support DOLFINx yet.
I made a small edit to the code by getting the UFL import right. And I also changed the boundary condition to a cantilever beam.

The following error still occurs:

*** Error:   Unable to create Dirichlet boundary condition.
*** Reason:  Expecting a scalar boundary value but given function is vector-valued.
*** Where:   This error was encountered inside DirichletBC.cpp.
*** Process: 0
4

You cannot enforce a 2D vector on the subspace 0:

as

is defined by

scalar elements/

1 Like
5

Dear Dokken,

I made a function space change for elements 1 and 2. Adjusted the material properties, entered the unit charge on the right side of the mesh, and it worked. However, as the result was not the same as that of the site of origin posted above (it is not possible to know if the result is exactly from the variables in question), I would like to know if any specialist could verify if the results make sense.

Here’s the new code:

from dolfin import *
import numpy
from ufl import indices

frequency=50000
omega=2*numpy.pi*frequency

L = 0.2
H = 0.01
Nx = 200
Ny = 10
mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny)

#Function space definition over the mesh
Element1 = VectorElement("CG", mesh.ufl_cell(), 1)
Element2 = VectorElement("CG", mesh.ufl_cell(), 1)

# Define the mixed element
V_elem = MixedElement([Element1, Element2])

# Define the mixed function space
Vcomplex = FunctionSpace(mesh, V_elem)

#Material properties
E, rho, nu = 300E+8, 8000, 0.3

# Lame coefficient for constitutive relation
mu = E/2./(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

def eps(v):
    return sym(grad(v))
def sigma(v):
    dim = v.geometric_dimension()
    return 2.0*mu*eps(v) + lmbda*tr(eps(v))*Identity(dim)


#Applying the Boundary conditions
zero=Constant((0.0, 0.0))
one=Constant((0.0, -1.0))

def left(x, on_boundary):
    return near(x[0], 0.)

def rigth(x, on_boundary):
    return near(x[0], 0.2)

boundary = [DirichletBC(Vcomplex.sub(0), zero, left), DirichletBC(Vcomplex.sub(1), zero, left)]
source = [DirichletBC(Vcomplex.sub(0), one, rigth), DirichletBC(Vcomplex.sub(1), one, rigth)]

bc=boundary+source

#Real and Imaginary parts of the trial and test functions
uR, uI = TrialFunctions(Vcomplex)
wR, wI = TestFunctions(Vcomplex)

F=omega**2*rho*(dot(uR,wR)+dot(uI,wI))*dx-(inner(sigma(uR), grad(wR))+inner(sigma(uI), grad(wI)))*dx

a, L = lhs(F), rhs(F)

# Set up the PDE
solnU = Function(Vcomplex)
problem = LinearVariationalProblem(a, L, solnU, bcs=bc)
solver = LinearVariationalSolver(problem)
solver.solve()
solnUR, solnUI=solnU.split()
plot(solnUI)
plot(solnUR)

file_results = XDMFFile(MPI.comm_world, "solution_frquency_300E+8.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.write(solnUI, 0.)
file_results.write(solnUR, 0.)

Result on paraview:

frequency0
frequency0790×632 36.1 KB

frequency
frequency790×632 47.6 KB