Geometric dimension error while compiling biharmonic pde

1

Hi all ,
I am working on biharmonic pde . I have neumann boundary condition . while compiling below code getting the error prescribed below
Cannot determine geometric dimension from expression.

Traceback (most recent call last):
  File "biharmonic3.py", line 87, in <module>
    + alpha/h_avg*inner(jump(z,n), jump(grad(v),n))*dS
  File "/usr/lib/python3/dist-packages/ufl/operators.py", line 381, in grad
    return Grad(f)
  File "/usr/lib/python3/dist-packages/ufl/differentiation.py", line 159, in __init__
    self._dim = find_geometric_dimension(f)
  File "/usr/lib/python3/dist-packages/ufl/domain.py", line 384, in find_geometric_dimension
    error("Cannot determine geometric dimension from expression.")
  File "/usr/lib/python3/dist-packages/ufl/log.py", line 172, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Cannot determine geometric dimension from expression.

My code is given below:

import matplotlib.pyplot as plt
import numpy as np
from dolfin import *

# Next, some parameters for the form compiler are set::

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True

# A mesh is created, and a quadratic finite element function space::

# Make mesh ghosted for evaluation of DG terms
parameters["ghost_mode"] = "shared_facet"

# Create mesh and define function space
mesh = UnitSquareMesh(8,8)
V = FunctionSpace(mesh, "CG", 2)


class Source(UserExpression):
    def eval(self, values, x):
        values[0] = 112*pi**4*(sin(pi*x[0]))**4*(sin(pi*x[1]))**4 + 24*pi**4*(cos(pi*x[0]))**4*(sin(pi*x[1]))**4 + 24*pi**4*(cos(pi*x[1]))**4*(sin(pi*x[0]))**4 - 288*pi**4*(cos(pi*x[0]))**2*(sin(pi*x[0]))**2*(sin(pi*x[1]))**4 -288*pi**4*(cos(pi*x[1]))**2*(sin(pi*x[0]))**4*(sin(pi*x[1]))**2 + 288*pi**4*(cos(pi*x[0]))**2*(cos(pi*x[1]))**2*(sin(pi*x[0]))**2*(sin(pi*x[1]))**2


# On the finite element space ``V``, trial and test functions are
# created::

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)

# Define normal component, mesh size and right-hand side
h = CellDiameter(mesh)
h_avg = (h('+') + h('-'))/2.0
n = FacetNormal(mesh)
g = Source(degree=2)

# Penalty parameter
alpha = Constant(6.0)

class Source(UserExpression):
    def eval(self, value, x):
        value[0] = (sin(pi*x[0]))*(sin(pi*x[0]))*(sin(pi*x[0]))*(sin(pi*x[0]))*(sin(pi*x[1]))*(sin(pi*x[1]))*(sin(pi*x[1]))*(sin(pi*x[1]))

f = Source(degree=2)
 
class Source(UserExpression):

    def eval(self, value, x):
        values[0] = 4*pi*cos(pi*x[0])*(sin(pi*x[0])**3)*(sin(pi*x[1])**4)
        values[1] = 4*pi*cos(pi*x[1])*(sin(pi*x[0])**4)*(sin(pi*x[1])**3)
    def value_shape(self):
         return (2,)

z = Source(degree=2)

b = inner(div(grad(f)), div(grad(v)))*dx \
  - inner(avg(div(grad(f))), jump(grad(v), n))*dS \
  - inner(jump(z, n), avg(div(grad(v))))*dS \
  + alpha/h_avg*inner(jump(z,n), jump(grad(v),n))*dS

a = inner(div(grad(u)), div(grad(v)))*dx \
  - inner(avg(div(grad(u))), jump(grad(v), n))*dS \
  - inner(jump(grad(u), n), avg(div(grad(v))))*dS \
  + alpha/h_avg*dot(jump(grad(u),n), jump(grad(v),n))*dS



# Define linear form
L = b - g*v*dx

# A :py:class:`Function <dolfin.functions.function.Function>` is created
# to store the solution and the variational problem is solved::

# Solve variational problem
u = Function(V)
solve(a == L, u)

x = SpatialCoordinate(mesh)
u_ex = (cos(pi*x[0]))*(cos(pi*x[1]))
error_L2 = sqrt(assemble(((u-u_ex)**2)*dx))
print('error_L2  =', error_L2)

# The computed solution is written to a file in VTK format and plotted to
# the screen. ::

# Save solution to file
file = File("biharmonic.pvd")
file << u

# Plot solution
plot(u)
plt.show()
2

You should be able to solve this by adding the following class to your userexpression:

def value_shape(self):
    return ()
3

Thanku for reply. Same error I am getting again. I want to know what is geometric dimension error?

4

Instead of using UserExpression, you can use ufl-expressions:

x = SpatialCoordinate(mesh)
f =(sin(pi*x[0]))*(sin(pi*x[0]))*(sin(pi*x[0]))*(sin(pi*x[0]))*(sin(pi*x[1]))*(sin(pi*x[1]))*(sin(pi*x[1]))*(sin(pi*x[1]))

z = as_vector((4*pi*cos(pi*x[0])*(sin(pi*x[0])**3)*(sin(pi*x[1])**4)
               ,4*pi*cos(pi*x[1])*(sin(pi*x[0])**4)*(sin(pi*x[1])**3)))
5

It work but L2 error is too large while I am using h=1/128 , I am getting 0.5000000952… . Am I doing something wrong? Kindly help me.

6

That means that you haven’t provided the correct sources for the analytical solution. I suggest debugging this by visualizing the exact solution (project it to an appropriate function space) and the approximate solution in i.e. Paraview

7

You can also consider this question on the old forum:
https://fenicsproject.org/qa/11686/biharmonic-equation-with-dirichlet-and-neumann-bcs/