Compute function maximum/minimum on subdomains

1

Hi all!
I’m trying to compute the maximum/minimum of a function on a specific subdomain (ie. not over the whole geometry).

I’m currently doing this:

from fenics import *

def calculate_minimum_volume(u, cell_function, subdomain_id):
    V = u.function_space()
    dofmap = V.dofmap()
mesh = V.mesh()
mini = u.vector().max()
for cell in cells(mesh):
    if cell_function[cell.index()] == subdomain_id:
        dofs = dofmap.cell_dofs(cell.index())
        for dof in dofs:
            try:
                [dof][0]
                if u.vector()[dof][0] < mini:
                    mini = u.vector()[dof][0]
            except:
                if u.vector()[dof] < mini:
                    mini = u.vector()[dof]
return mini

mesh = UnitSquareMesh(5, 5)
V = FunctionSpace(mesh, 'P', 1)
u = Function(V)

mf = MeshFunction("size_t", mesh, 2, 0)


domain = CompiledSubDomain('x[1] > 0.5')
domain.mark(mf, 1)

print(calculate_minimum_volume(u, mf, 1))

Though on big meshes it can take ages. Do you guys know a better way of doing this ?

Cheers,

Rem

2

The following seems snappy enough

from dolfin import *
import numpy as np


def minimum(f, subdomains, subd_id):
    '''Minimum of f over subdomains cells marked with subd_id'''
    V = f.function_space()
    
    dm = V.dofmap()

    subd_dofs = np.unique(np.hstack(
        [dm.cell_dofs(c.index()) for c in SubsetIterator(subdomains, subd_id)]))
    
    return np.min(f.vector().get_local()[subd_dofs])

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

if __name__ == '__main__':
    
    domain = CompiledSubDomain('x[1] > 0.5 - DOLFIN_EPS')
    f = Expression('x[0]+x[1]', degree=1)
    
    # Check scaling
    for n in (8, 16, 32, 64, 128, 256):
        
        mesh = UnitSquareMesh(n, n)
        V = FunctionSpace(mesh, 'P', 1)
        u = interpolate(f, V)

        mf = MeshFunction('size_t', mesh, 2, 0)
        domain.mark(mf, 1)

        t = Timer('minimum')
        min_value = minimum(u, mf, 1)
        dt = t.stop()

        print dt, min_value
2 Likes
3

Works like 100 times faster. Thanks!

4

I’ve tried to change this function to get the maximum but it doesn’t seem to be working :confused:
Do you know what I’ve done wrong ?

def maximum(f, subdomains, subd_id):
    '''Maximum of f over subdomains cells marked with subd_id'''
    V = f.function_space()
    dm = V.dofmap()
    subd_dofs = np.unique(np.hstack(
        [dm.cell_dofs(c.index()) for c in SubsetIterator(subdomains, subd_id)]))
    return np.max(f.vector().get_local()[subd_dofs])

The error message is:
index -2147483648 is out of bounds for axis 0 with size 10201

5

I cannot reproduce your behavior. Just by changing Mirok’s code to maximum:

from dolfin import *
import numpy as np


def maximum(f, subdomains, subd_id):
    '''Minimum of f over subdomains cells marked with subd_id'''
    V = f.function_space()

    dm = V.dofmap()

    subd_dofs = np.unique(np.hstack(
        [dm.cell_dofs(c.index()) for c in SubsetIterator(subdomains, subd_id)]))

    return np.max(f.vector().get_local()[subd_dofs])

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

if __name__ == '__main__':

    domain = CompiledSubDomain('x[1] > 0.5 - DOLFIN_EPS')
    f = Expression('x[0]+x[1]', degree=1)

    # Check scaling
    for n in (8, 16, 32, 64, 128, 256):

        mesh = UnitSquareMesh(n, n)
        V = FunctionSpace(mesh, 'P', 1)
        u = interpolate(f, V)

        mf = MeshFunction('size_t', mesh, 2, 0)
        domain.mark(mf, 1)

        t = Timer('max')
        max_value = maximum(u, mf, 1)
        dt = t.stop()

        print( dt, max_value)

I get the following output

Iterating over subset, found 64 entities out of 128.
0.00044637500000000005 2.0
Iterating over subset, found 256 entities out of 512.
0.000759978 2.0
Iterating over subset, found 1024 entities out of 2048.
0.0024261440000000003 2.0
Iterating over subset, found 4096 entities out of 8192.
0.009333566 2.0
Iterating over subset, found 16384 entities out of 32768.
0.035768256000000005 2.0
Iterating over subset, found 65536 entities out of 131072.
0.13794973700000002 2.0