Calculating an exact solution on mesh vertices with a function defined

1

Hi,

I have an exact solution to a problem which I must calculate on mesh vertices, and I cannot simply use an Expression to define it. So I have to define it as a function. I think my problem when I solve it the solution is kind of “discretized” and I cannot “interpolate” on my domain. I do not understand what is causing it. I’d appreciate any help. Here is a simply example code to what I am trying to achieve:

from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from dolfin import *

# Exact Solution Defined as a Function

def Analytic_Solution(x,y):
    
    #this is just to show what I am trying to do.
    u_ex = x**2+y**2
    
    return u_ex

# Define geometry
radius_c = 1
c_points = 30
xcyl = 0
ycyl = 0
domain = Circle(Point(xcyl, ycyl), radius_c, c_points)

# Mesh the domain
mesh = generate_mesh(domain,10)
plot(mesh)

# I need help here

F = FunctionSpace(mesh, 'CG', 1)

n = F.dim()                                                                      
d = mesh.geometry().dim()

f = Function(F)

F_dof_coordinates = F.tabulate_dof_coordinates()
F_dof_coordinates.resize((n, d))
F_dof_x = F_dof_coordinates[:, 0]
F_dof_y = F_dof_coordinates[:, 1]

# This part does not give a smooth solution
u_temp = Analytic_Solution(F_dof_x,F_dof_y)

f.vector().set_local(u_temp)

plot(f)
2

Please Edit your post and encapsulate the code with ``` to make it readable

1 Like
3

Which version of fenics are you using? The resulting plot looks like this for me (using the docker image quay.io/fenicsproject/stable:latest). This does not seem wrong to me?

circle_int
circle_int798×721 179 KB

1 Like
4

This is the result I get too. Well, to my surprise I understand my problem is something else now.

I checked if I will get any errors with an exact solution defined with an expression as below:

# Exact solution defined with an expression
V = FunctionSpace(mesh, 'P', 1)
u = interpolate(Expression('x[0]*x[0] + x[1]*x[1]', degree=2), V)
plot(u, interactive=True)

# Check if there is any error
E = errornorm(u, f)
print(E)

Looks like either I was not able to find the source of my error or something I do here is different and I get the correct answer by chance. I suspect about the change of ordering of nodal values and vertex values, I am confused about that part at the moment. Anyway, I hope this helps other people. I will add my question here later if I can find the “real problem” in my actual code. I am not having the right answer at the moment in my code but this one here looks alright. Thank you dokken anyway.