Derivate of Expression

1

Hello everybody,
I have a Function f : \Omega \subset \mathbb{R}^3 \to \mathbb{R} with

f \in H^1_0(\Omega) and i want to calculate \nabla f \in H_0(curl)
My function is
f(x)=sin(\pi \cdot x_1)\cdot sin(\pi \cdot x_2)\cdot sin(\pi\cdot x_3)
and I know that
\nabla f(x)= \left( \begin{array}{c} \pi\cdot cos(\pi\cdot _1)\cdot sin(\pi \cdot x_2)\cdot sin(\pi \cdot x_3) \\ \pi \cdot sin(\pi \cdot x_1) \cdot cos(\pi \cdot x_2)\cdot sin(\pi \cdot x_3) \\ \pi \cdot sin(\pi \cdot x_1) \cdot sin(\pi \cdot x_2)\cdot cos(\pi \cdot x_3) \end{array}\right)

I have some questions:
-how I have to define f? I tried

pi = math.pi

eta = Expression('sin(pi*x[0])+sin(pi*x[1])+sin(pi*x[2])', degree=1)`

Is this so ok?

-how can i set \nabla f manually and what I have to do so that Fenics computes it automatically?

from dolfin import*
mesh = UnitCubeMesh(10, 10, 10)

V = FiniteElement("N1curl", mesh.ufl_cell(), 1)
S = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
Y = VectorFunctionSpace(mesh, 'P', 1)
Z = V * S
W = FunctionSpace(mesh, Z)

Thanks for your help.
Noya

2

Consider the following:

from dolfin import *

mesh = UnitCubeMesh(10, 10, 10)

n1curl = FiniteElement("N1curl", mesh.ufl_cell(), 1)
cg1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
Z = MixedElement([n1curl, cg1])

W = FunctionSpace(mesh, Z)
w = Function(W)

W0 = W.sub(0).collapse()

x = SpatialCoordinate(mesh)
eta = sin(pi*x[0])+sin(pi*x[1])+sin(pi*x[2])
grad_eta = project(grad(eta), W0)
assigner = FunctionAssigner(W.sub(0), W0)
assigner.assign(w.sub(0), grad_eta)
print(w.vector().get_local())
2 Likes
3

Thanks for your answer and your help. This works very well.
But in the end I want to calculate a norm like \parallel eta - p_h\parallel_{H^1_0} where p_h is solution from my code and eta is the analytical solution.

(yh, ph) = TrialFunctions(W)

But it doesn’t work and the error is:

Traceback (most recent call last):
File “/home/noya/test.py”, line 91, in
error_p_exact = errornorm(eta, ph, ‘H10’)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/norms.py”, line 246, in errornorm
degree_u = u.ufl_element().degree()
AttributeError: ‘Sum’ object has no attribute ‘ufl_element’

And I want to set the analytical solution y as zero but when I try something like

y = Constant((0, 0, 0)

I’ve got the following error when I try to compute the norm \parallel y - y_h\parallel_{H_0(curl)}

Traceback (most recent call last):
File “/home/noya/test.py”, line 89, in
error_y_exact = errornorm(y, yh, ‘Hcurl’)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/norms.py”, line 233, in errornorm
raise RuntimeError(“Cannot compute error norm. Value shapes do not match.”)
RuntimeError: Cannot compute error norm. Value shapes do not match.

I’d would be nice if you could help me to understand and correct my mistakes.
Noya

4

As you are not supplying the code reproducing the error, I cannot really help you any further.
Please note that you cannot use a TrialFunction in your error-calculation (as you have indicated that you do with your snippets).

6

You have redefined y

as

in your code, and thus y has 4 components, (3 in first space and 1 in the second space), while yh has 3 components.

8

I would suggest the following:

from dolfin import *
mesh = UnitCubeMesh(6, 6, 6)
print('dim mesh:', mesh.geometry().dim())  # Dimension mesh


V = FiniteElement("N1curl", mesh.ufl_cell(), 1)
S = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
Z = MixedElement([V, S])
W = FunctionSpace(mesh, Z)

print('Dim = DoF: ', W.dim())  # Dimension W


def boundary(x, on_boundary):
    return on_boundary


w = Function(W)
bc = DirichletBC(W, w, boundary)

(vh, qh) = TestFunctions(W)
(yh, ph) = TrialFunctions(W)


y = Constant((0, 0, 0))

W0 = W.sub(0).collapse()
x = SpatialCoordinate(mesh)
eta = sin(pi*x[0])+sin(pi*x[1])+sin(pi*x[2])
eta_expr = project(eta, FunctionSpace(mesh, "CG", 4))
grad_eta = project(grad(eta), W0)
assigner = FunctionAssigner(W.sub(0), W0)
assigner.assign(w.sub(0), grad_eta)
print(w.vector().get_local())
p = eta
# define f, nu

f = grad_eta
nu = Expression('1-(1/(2*(s*s+1)))', degree=2, s=0)

# weak formulation
# left side
a = nu * inner(curl(yh), curl(vh)) * dx
b = inner(vh, grad(ph)) * dx
c = inner(yh, grad(qh)) * dx
#A = a + b + c
A = a+b+c
# right side
L = inner(f, vh) * dx

# initial value
z0 = Function(W)
y0, p0 = z0.split()
nu.s = norm(y0, 'Hcurl0')


maxiter = 100
tol = 1e-8

# Iteration
m = Function(W)
for n in range(0, maxiter):
    solve(A == L, m, bc)
    yh, ph = m.split()
    error1 = errornorm(yh, y0, 'Hcurl')

    error2 = errornorm(ph, p0, 'H10')

    error = error1 + error2

    error_y_exact = errornorm(y, yh, 'Hcurl')
    error_p_exact = errornorm(eta_expr, ph, 'H10')

    error_exact = error_y_exact + error_p_exact
    print(n, error_exact)
    if error < tol:
        # if error & error_exact < tol
        break
    nu.s = norm(yh, 'Hcurl0')
    z0.assign(m)

However, as your variational form is not dependent of z0 (or y0 and p0), I do not see how you are expecting anything to change from iteration to iteration.