# Spheres Don't Looks Spherical

Hi

All. I am using fenics/dolfin version 2019.1.0.

I am trying to model the waves due to a moving spherical body. However, I am not sure how to correctly define a spherical object using dolfin expressions. The code below creates a sphere, but as you can see, when I project the expression to a function space and plot it in 2D, the spherical body is far from spherical and highly dependent on mesh resolution. Does anyone know how to properly define a spherical object to avoid this issue? Thanks!

``````import dolfin as d
import matplotlib.pyplot as plt
import numpy as np

## Creating a simple square mesh:
mesh = d.UnitSquareMesh(50, 50)

## An expression for a spherical object centered at (0.5, 0.5) with radius 0.1:
sph = d.Expression("pow(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2), 0.5) <= 0.1 ? 100 : 0", degree=2)

## Projecting into a function space and plotting:
V = d.FunctionSpace(mesh, "CG", 2)
sph_proj = d.project(sph, V)
max = np.round(sph_proj.vector().max(),2) ## maximum value of the projection

plt.figure()
col = d.plot(sph_proj)
d.plot(mesh, alpha = 0.5)
plt.title(f'Sph. source [max val = {max}]')
plt.colorbar(col)
plt.show()
``````

The most accurate representation would be to use `SpatialCoordinate(mesh)`, as this will be evaluated at the quadrature points of the mesh.
You can control the quadrature degree and thus the resolution of the sphere.

This is expected. Please note that matplotlib plotting

interpolates the solution into P-1. you should use `XDMFFile.write_checkpoint` and paraview to retain the second order polynomials.

Hi, thanks, I was not familiar with the `SpatialCoordinate` function. A quick follow up question – I am not entirely sure how to use the coordinates obtained from `SpatialCoordinate` as a part of an expression. Here is my naive attempt to do it:

``````import dolfin as d
import matplotlib.pyplot as plt
import numpy as np

## Creating a simple square mesh:
mesh = d.UnitSquareMesh(512, 512)
x_phys = d.SpatialCoordinate(mesh)
x = x_phys

class sph(d.UserExpression):
def eval(self, value, x):
value[0] = 100 if (pow(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2), 0.5) <= 0.1) else 0
def value_shape(self):
return ()

sph0 = sph()
V = d.FunctionSpace(mesh, "DG", 0)
dens = d.interpolate(sph0, V)
max = np.round(dens.vector().max(),2) ## maximum value of the projection

plt.figure()
col = d.plot(dens)
plt.title(f'Sph. source [max val = {max}]')
plt.colorbar(col)
plt.show()
``````

Is this the correct way of using the x’s from `SpatialCoordinate`? Also I have noticed that the results obtained by `dens = d.interpolate(sph0, V)` sometimes (depending on the resolution) differ from `dens = d.project(sph0, V)`. Which of these should I be using in this case? Thanks!

Consider the following MWE:

``````import dolfin as d
import matplotlib.pyplot as plt
import numpy as np
try:
from ufl import conditional
except ModuleNotFoundError:
from ufl_legacy import conditional

## Creating a simple square mesh:
mesh = d.UnitSquareMesh(512, 512)
x_phys = d.SpatialCoordinate(mesh)

sph0 = conditional(d.lt(d.sqrt((x_phys[0] - 0.5)**2 + (x_phys[1] - 0.5)**2), 0.1), 100, 0)

V = d.FunctionSpace(mesh, "DG", 0)
dens = d.project(sph0, V)
max = np.round(dens.vector().max(),2) ## maximum value of the projection

plt.figure()
col = d.plot(dens)
plt.title(f'Sph. source [max val = {max}]')
plt.colorbar(col)
plt.savefig('sph.png')
``````

which yields

Project and interpolation are two different operations. See for instance:
https://jsdokken.com/FEniCS23-tutorial/src/approximations.html
that explains some of the things happening behind the scenes.