SubDomain as a function of the solution

mesh
1

Hi everyone,

I’m working on a problem using the Cahn-Hilliard equation and I want to keep track of the mass of each pure phase on each time step, that is to say the domains where c = 0.0 and c=1.0. Since my transition region is really small compared to the pure phase, my idea was to define a SubDomain where c < 0.5 and another where c > 0.5 .

My working example interpolates a circle of radius 0.25 and c = 1.0 as the solution u and creates a SubDomain wherever u > 0.0 by looping through all the cells. If there is a more efficient way to do this, I’d really appreciate it.

Cheers

from dolfin import *
import math
import time

tic     = time.time()

class ExpessionTest(UserExpression):
    def __init__(self, **kwargs):
        super().__init__(**kwargs)
    def eval_cell(self, value, x, ufc_cell):
        if ( (x[0] - .5)**2 + (x[1] - .5)**2 ) < 0.25**2:
            value[0] = 1.0
        else:
            value[0] = 0.0
    def value_shape(self):
        return ()

N                   = 500
mesh                = UnitSquareMesh(N, N)
V                   = FunctionSpace(mesh, "CG", 1)
u                   = Function(V)
u                   = interpolate(ExpessionTest(degree = 1), V)
File("Value.pvd")   << u

sub_domains         = MeshFunction("size_t", mesh, mesh.topology().dim())
sub_domains.set_all(0)

for c in cells(mesh):
    md      = c.midpoint()
    val     = u( md.x(), md.y() )
    if val > 0.0:
        sub_domains[c]  = 1

File("Subdomains.pvd") << sub_domains

dv = Measure('dx', domain = mesh, subdomain_data = sub_domains)

Vol_n   = assemble(u*dv(1))
Vol_a   = math.pi * .25**2
error   = abs(Vol_a - Vol_n) / Vol_a

toc     = time.time()

print("Numerical volume: %.5e   Analytical volume: %.5e Error: %.5e" % (Vol_n, Vol_a, error))
print("Simulation time in minutes: %.2f" % ((toc-tic)/60.0) )
2

Hi, consider

from dolfin import *
import numpy as np
import time


f = Expression('pow(x[0]-0.5, 2) + pow(x[1]-0.5, 2) < pow(0.25, 2) ? 1: 0', degree=1)


def fill0(subdomains, u):
    mesh = subdomains.mesh()
    for c in cells(mesh):
        md = c.midpoint()
        val = u( md.x(), md.y() )
        if val > 0.0:
            sub_domains[c]  = 1

    return subdomains


def fill1(subdomains, u):
    mesh = subdomains.mesh()
    Q = FunctionSpace(mesh, 'DG', 0)
    q = TestFunction(Q)
    # Project to P0 space. If u linear integration mounts to eval at
    # cell midpoint
    b = assemble((1/CellVolume(mesh))*inner(u, q)*dx)

    subdomains.array()[:] = np.where(b.get_local() > 0.0, 1.0, 0.0)
    
    return subdomains

# --------------------------------------------------------------------

if __name__ == '__main__':
    N = 1024
    mesh = UnitSquareMesh(N, N)
    V = FunctionSpace(mesh, "CG", 1)
    u = interpolate(f, V)
    

    sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
    
    tic = time.time()
    sub_domains = fill0(sub_domains, u)
    toc = time.time()
    dv = Measure('dx', domain = mesh, subdomain_data = sub_domains)

    Vol_n   = assemble(u*dv(1))
    Vol_a   = pi * .25**2
    error   = abs(Vol_a - Vol_n) / Vol_a

    print("Numerical volume: %.5e   Analytical volume: %.5e Error: %.5e" % (Vol_n, Vol_a, error))
    print("Simulation time in s: %.2f" % (toc-tic) )

    sub_domains.set_all(0)

    tic = time.time()
    sub_domains = fill1(sub_domains, u)
    toc = time.time()
    dv = Measure('dx', domain = mesh, subdomain_data = sub_domains)

    Vol_n   = assemble(u*dv(1))
    Vol_a   = pi * .25**2
    error   = abs(Vol_a - Vol_n) / Vol_a

    print("Numerical volume: %.5e   Analytical volume: %.5e Error: %.5e" % (Vol_n, Vol_a, error))
    print("Simulation time in s: %.2f" % (toc-tic) )
3

Thanks! That is a lot faster than my method. For the sake of generality, if I where to use a higher order element, would another faster method exists?

EDIT: Assuming that the subdomain is a growing circle with center in the upper edge of a rectangle domain and that I wanted to output the radius of said circle, is there a way to do that without a CellIterator?

The snippet of a working code would be

size            = np.sum(sub_domains.array())
y_coord         = np.zeros(size)
itr             = 0
for cell in SubsetIterator(sub_domains, 1):
    y_coord[itr] = cell.midpoint().y()
    itr          += 1
Radius          = height - np.amin(y_coord)