Level set based domain definition in FEniCS (3D)

[Guangjin_MOU]

Hello FEniCSers,

I want to use level set function to define two different phases in a 1* 1* 1 bloc. My level set function is:
lsf2[x_, y_, z_] =Cos[Pi*(x - 0.5)]^2Cos[Pi(y - 0.5)]^2Cos[Pi(z - 0.5)]^2 - 0.5
when I plot it in Mathematica, that would be a sphere inclusion in a 1* 1* 1 bloc:

but the outcome in FEniCS is different, my whole code is:

from __future__ import print_function
from fenics import *
import numpy as np,sys
import matplotlib.pyplot as plt
import mshr
from matplotlib.pyplot import figure
from matplotlib.ticker import FuncFormatter
from scipy.integrate import dblquad

Nx=Ny=Nz=20
deg=1
lsf=Expression('1/N*(pow((x[0]-0.5),2)+pow(((x[1]-0.5)),2)+pow((x[2]-0.5),2)-0.1)',N=1,degree=1)
mesh=BoxMesh(Point(0.0,0.0,0.0),Point(1.0,1.0,1.0),Nx,Ny,Nz)
materials= MeshFunction('size_t', mesh, 3)  
dx=Measure("dx")(subdomain_data=materials,domain=mesh)
V=FunctionSpace(mesh,'CG',deg)
lsf_project=Function(V)  
lsf_project.assign(project(lsf,V))

class Omega_0(SubDomain):  #inclusion
  def inside(self, x, on_boundary):
        return lsf_project(x[0],x[1],x[2])<=1E-14

materials.set_all(1)
subdomain_0 = Omega_0()
subdomain_0.mark(materials, 0)
xdmffile_materials=XDMFFile('materials_phase.xdmf')
xdmffile_materials.write(materials)
print('phase file saved correctly')

And I illustrate the material_phase.xdmf file in Paraview, which results in:

As we can see, the inclusion is not symmetric as it should be (ex. top right vs top left). Is it a normal phenomena? How to avoid it?

Thanks

[Guangjin_MOU]

By the way, there is no such problem in 2D case

[dokken]

Why should the inclusion be symmetric? You have a grid dependency, as you check if every vertex of your mesh satisfies the criterion in question.
If you look at the levelset in question:

from dolfin import *
Nx=Ny=Nz=20
deg=1
lsf=Expression('1/N*(pow((x[0]-0.5),2)+pow(((x[1]-0.5)),2)+pow((x[2]-0.5),2)-0.1)',N=1,degree=deg)
mesh=BoxMesh(Point(0.0,0.0,0.0),Point(1.0,1.0,1.0),Nx,Ny,Nz)
materials= MeshFunction('size_t', mesh, 3)  
dx=Measure("dx")(subdomain_data=materials,domain=mesh)
V=FunctionSpace(mesh,'CG',deg)
lsf_project=Function(V)  
lsf_project.assign(project(lsf,V))

class Omega_0(SubDomain):  #inclusion
  def inside(self, x, on_boundary):
        return lsf_project(x[0],x[1],x[2])<=1E-14

materials.set_all(1)
subdomain_0 = Omega_0()
subdomain_0.mark(materials, 0)
xdmffile_materials=XDMFFile('materials_phase.xdmf')
xdmffile_materials.write(materials)
print('phase file saved correctly')


xdmffile_lsf=XDMFFile('levelset.xdmf')
xdmffile_lsf.write_checkpoint(lsf_project, "levelset", 0.0, append=False)

You will see that your leveset function has the same behavior. You can for instance increase the tolerance of Omega_0 slightly (from 1e-14 to 1e-10) to get a slightly modified figure.