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!
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')
