Accuracy on subdomain area

variational formulation
1

Hello,

I am working with a level set formulation, in which the domain of interest is the zero level set of a function Phi. I need to calculate the area of this domain. However, since I’m using a binary strategy for marking the inside and outside of the domain, the area value is not accurate unless I use really refined meshes. I wonder if there is a simple way of improving this result, for exemple using intermediate values of density to each element instead of a binary approach.

Here I have a simple exemple in which I calculate the area of a circle increasing the mesh size and comparing to the expected value.

from dolfin import *
import numpy as np
from matplotlib import pyplot as plt
from tabulate import tabulate

r = 0.2
def phi(x):
    return (x[0] - 0.5)**2 + (x[1]-0.5)**2 - r**2 #Zero level set is a circle centered in (0.5,0.5) with radius 0.2
class Omega(SubDomain):
    def inside(self, x, on_boundary):
        return phi(x) > 0

#Total domain area: 1
#Interior expected area: pi*r^2
#Exterior expected area: 1 - pi*r^2

a = 10; q = 2; length = 6
size = [a *  q** (n - 1) for n in range(1, length + 1)]
a_int = []; a_ext = []; e_int = []; e_ext=[]; k = 0

for N in size:
    print('Mesh size:', N)
    mesh = UnitSquareMesh(N, N,'crossed')
    domains = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
    omega = Omega() 
    omega.mark(domains, 1)  
    dx_int = Measure("dx", domain=mesh, subdomain_data=domains)
    a_int.append(assemble(Constant(1)*dx_int(0)))
    a_ext.append(assemble(Constant(1)*dx_int(1)))
    e_int.append(pi*r**2 - a_int[k])
    e_ext.append(1 - pi*r**2 - a_ext[k])
    print('Interior area:', a_int[k])
    print('Error on interior area:', e_int[k])
    print('Exterior area:', a_ext[k])
    print('Error on extrior area:', e_ext[k])
    k+=1

e_size = np.ones(6)/size
headers=['mesh size','element size', 'interior area', 'exterior area', 'interior error', 'exterior error']
table = [size, e_size,a_int,a_ext, e_int, e_ext]
print(tabulate(np.transpose(table),headers=headers, floatfmt=".4f"))

Thank you!

1 Like
2

I’m far from an expert on level set methods, but I thought the idea was to use the function itself as a definition of the geometry of the interface \vec{x}_\text{interface} = \{\vec{x} : \phi(\vec{x}) = 0\}.

To this end I just compute (please forgive the haphazard notation)
A_\text{interior} = \int_\Omega \begin{cases} 1 & \phi < 0 \\ 0 & \text{otherwise} \end{cases} \mathrm{d}x
and
A_\text{exterior} = \int_\Omega \begin{cases} 1 & \phi \geq 0 \\ 0 & \text{otherwise} \end{cases} \mathrm{d}x.

See the following modifications:

from dolfin import *
import numpy as np
from matplotlib import pyplot as plt
from tabulate import tabulate

r = 0.2
def phi(x):
    return (x[0] - 0.5)**2 + (x[1]-0.5)**2 - r**2 #Zero level set is a circle centered in (0.5,0.5) with radius 0.2
class Omega(SubDomain):
    def inside(self, x, on_boundary):
        return phi(x) > 0


class PhiExpr(UserExpression):
    def eval(self, value, x):
        value[0] = phi(x)
    def value_shape(self):
        return ()

#Total domain area: 1
#Interior expected area: pi*r^2
#Exterior expected area: 1 - pi*r^2

a = 10; q = 2; length = 6
size = [a *  q** (n - 1) for n in range(1, length + 1)]
a_int = []; a_ext = []; e_int = []; e_ext=[]; k = 0
algebraic_int, algebraic_ext, e_algebraic_int, e_algebraic_ext = [], [], [], []

for N in size:
    print('Mesh size:', N)
    mesh = UnitSquareMesh(N, N,'crossed')

    V = FunctionSpace(mesh, "CG", 1)
    phi_function = Function(V)
    phi_function.interpolate(PhiExpr())

    algebraic_int.append(assemble(conditional(lt(phi_function, 0.0), 1, 0)*dx))
    algebraic_ext.append(assemble(conditional(gt(phi_function, 0.0), 1, 0)*dx))
    e_algebraic_int.append(pi*r**2 - algebraic_int[k])
    e_algebraic_ext.append(1 - pi*r**2 - algebraic_ext[k])

    domains = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
    omega = Omega() 
    omega.mark(domains, 1)  
    dx_int = Measure("dx", domain=mesh, subdomain_data=domains)
    a_int.append(assemble(Constant(1)*dx_int(0)))
    a_ext.append(assemble(Constant(1)*dx_int(1)))
    e_int.append(pi*r**2 - a_int[k])
    e_ext.append(1 - pi*r**2 - a_ext[k])
    print('Interior area:', a_int[k])
    print('Error on interior area:', e_int[k])
    print('Exterior area:', a_ext[k])
    print('Error on extrior area:', e_ext[k])

    print('Interior algebraic area:', algebraic_int[k])
    print('Error on algebraic interior area:', e_algebraic_int[k])
    print('Exterior algebraic area:', algebraic_ext[k])
    print('Error on algebraic exterior area:', e_algebraic_ext[k])

    k+=1

e_size = np.ones(length)/size
headers=['mesh size','element size', 'interior area', 'exterior area', 'interior error', 'exterior error']
table = [size, e_size,a_int,a_ext, e_int, e_ext]
print(tabulate(np.transpose(table),headers=headers, floatfmt=".3e"))

headers2=['mesh size','element size', 'algebraic interior area', 'algebraic exterior area', 'algebraic interior error', 'algebraic exterior error']
table2 = [size, e_size,algebraic_int,algebraic_ext, e_algebraic_int, e_algebraic_ext]
print(tabulate(np.transpose(table2),headers=headers2, floatfmt=".3e"))

where table2 gives a slightly better approximation

  mesh size    element size    interior area    exterior area    interior error    exterior error
-----------  --------------  ---------------  ---------------  ----------------  ----------------
  1.000e+01       1.000e-01        1.800e-01        8.200e-01        -5.434e-02         5.434e-02
  2.000e+01       5.000e-02        1.500e-01        8.500e-01        -2.434e-02         2.434e-02
  4.000e+01       2.500e-02        1.375e-01        8.625e-01        -1.184e-02         1.184e-02
  8.000e+01       1.250e-02        1.313e-01        8.687e-01        -5.586e-03         5.586e-03
  1.600e+02       6.250e-03        1.286e-01        8.714e-01        -2.930e-03         2.930e-03
  3.200e+02       3.125e-03        1.271e-01        8.729e-01        -1.465e-03         1.465e-03
  mesh size    element size    algebraic interior area    algebraic exterior area    algebraic interior error    algebraic exterior error
-----------  --------------  -------------------------  -------------------------  --------------------------  --------------------------
  1.000e+01       1.000e-01                  1.400e-01                  8.600e-01                  -1.434e-02                   1.434e-02
  2.000e+01       5.000e-02                  1.250e-01                  8.750e-01                   6.637e-04                  -6.637e-04
  4.000e+01       2.500e-02                  1.250e-01                  8.750e-01                   6.637e-04                  -6.637e-04
  8.000e+01       1.250e-02                  1.253e-01                  8.747e-01                   3.512e-04                  -3.512e-04
  1.600e+02       6.250e-03                  1.256e-01                  8.744e-01                   3.871e-05                  -3.871e-05
  3.200e+02       3.125e-03                  1.257e-01                  8.743e-01                  -3.564e-07                   3.563e-07

You can probably get even better results by modifying the space in which \phi lives or some fancy quadrature scheme.

I’m pretty sure you can do wonderful things with splines to get even more precise and smooth representations of level sets using splines. See for example @kamensky’s tIGAr.

1 Like
3

Hi Nate,

Thank you so much for your reply. This computation you mention is exactly what I was trying to do with the class Omega. But I was getting a binary choice for each element: 1 for the interior and 0 for the exterior. Your approach seems to be much wiser, although I’m not sure if I totally get it (I’m sorry, I am still a bit new to FEniCS).

Is this line phi_function.interpolate(PhiExpr()) actually interpolating Phi function inside each element as if I was getting a density value for each element? Or does it only consider nodal values of Phi on the mesh?

Thank you much!