integrate along single facet

Question

Hi,

I would like to integrate the quantity avg(u) of a discontinuous solution 'u' across a single facet in my mesh, but cant seem to figure out how to do this. How would I go about this?

Best,
Stein

Comments

Is it a single facet of an element in particular or is it a particular face of, say, a square?

---

Essentially I just want to loop over all edges/facets in the mesh, and do a certain integration

Answer 1

What about this:

from dolfin import *
import numpy as np

V = FunctionSpace(UnitSquareMesh(4,4, 'crossed'), 'DG', 0)
f = Expression("x[0]", degree=1)
u = project(f, V)

mf = MeshFunctionSizet(V.mesh(), 1, 0)  # Default value for facets is set to 0
facet_ids = np.zeros(V.mesh().num_facets(), dtype=np.uintp)
facet_ids[::6] = 1   # <<---- Mark those you want
mf.set_values(facet_ids)  # Mark your facet ids with 1
# Alternatively:
#mf.set_value(one_facet_id, 1)

assemble(avg(u)*dS(subdomain_data=mf, subdomain_id=1)) 

And (bonus! ;) in case you want to check what facets you are marking, you can use this:

import matplotlib.pyplot as pl
def plot_facets(mesh:Mesh, ids:list):
    """ Plots the mesh and highlights the facets given by ids."""
    plot(mesh, alpha=0.5)
    coords = mesh.coordinates()
    for fid in ids:
        f = Facet(mesh, fid)
        v1, v2 = f.entities(0)
        cc = np.vstack(coords[[v1,v2]])
        pl.plot(cc[:,0], cc[:,1], c='red')

plot_facets(V.mesh(), np.where(facet_ids == 1)[0])

Comments

Thanks for having a look at this. Your approach seems very promising. But i'm getting the error "Unable to extract all indices". I'm not very confident about my own understanding of the internal workings of MeshFunctions. However, shouldnt it be:

facet_ids[6] = 1

Instead of

facet_ids[::6] = 1

And why aren't we using a FacetFunction instead of a MeshFunction?
I still cant get it to work though. Right now I am interested in comparing integrals across facets of an element, so:

mf = MeshFunctionSizet(V.mesh(), 1, 0)
for c in cells(V.mesh()):
    print "Cell:", c.index()
    for fct in facets(c):
        facet_ids = np.zeros(V.mesh().num_facets(), dtype=np.uintp)
        facet_ids[fct.index()] = 1
        mf.set_values(facet_ids)

        print assemble((avg(u))*dS(subdomain_data=mf, subdomain_id=1))
    print ''

---

I suppose there will be some more efficient way of doing it, but your code does work for me. (But if I may suggest a micro-optimization: you can declare facet_ids before the loops to avoid constructing it in every iteration if you reset the only value that changed after the assembly: facet_ids[fct.index()] = 0. Of course the cost is negligible but hey... ;)

Also: facet_ids[::6] is just taking every sixth element of the array, it was just an example.

And yes FacetFunctionSizet should work exactly the same (you just don't need to specify the dimension anymore, so the call would be mf = FacetFunctionSizet(V.mesh(), 0)), I just wanted to demo the more general concept.

Good luck!

---

In my implementation, which failed, I was accidentally using a trialfunction instead of a regular function for 'u'. Hence the error message. Now its working. I'm not too worried about optimization as I just mean to do some simple analysis on a very coarse mesh for some benchmark problems.
Thanks for your help!