cell index for use in a function form & ('-') or ('+') for the exterior facets

Question

Hi. I have two questions (both for DG elements):

1) I need to create a function form to evaluate $$\frac{A}{V}$$ for each cell. I use this for an Interior-Penalty term. Below is what I've written, but I'd like to make it as general as possible. So, I think, instead of cell_idx below, I should use the local cell index number as I do not know the cell numbering on each processors. Any thoughts what is the best approach?

def AreaOverVolume(mesh):
    V0 = FunctionSpace(mesh, "DG", 0)
    AoV=Function(V0)
    D = mesh.topology().dim()
    for cell_idx in xrange(mesh.num_cells()):
        cell = Cell(mesh, cell_idx)
        facet = cell.entities(D-1)
        area = 0.0
        for f in xrange(cell.num_entities(D-1)):
            area += cell.facet_area(f) # perimeter in 2D
        AoV.vector()[cell_idx] = area/cell.volume()
    return AoV

2) In the DG formulation, we have a jump across facets. For interior facets, it is not really important which cell is ('+') or ('-'). For the exterior facets (ds), however, I think this distinction is important. In my little experience with FEniCS, it appears that using either one or even the average of the two returns the same identical answer. I need help understanding what is the appropriate choice; e.g., AoV('-')*ds or AoV('+')*ds

Regards,

Alireza

Comments

I don't fully understand your question, but are you looking for something different from:

h = CellVolume(mesh)/FacetArea(mesh)

Answer 1

1) Your function works the same on one or many processors. The num_cells() and cell_indices etc will refer to the local processor only. The function AoV is global as always and since all cells belong on some processor then all cells will get a value. The only thing is that you need to not use cell_index as the dof number when you call AoV.vector()[dof_number] = area/cell.volume(). You get the proper dof_number from V0.dofmap().cell_dofs(cell_index)

2) There is only one side with any value on an external ds facet, so you must write only AoV and not include any restriction + or -. It is only in dS integrals where you need to put restrictions on the integrand.

Comments

I am not 100% sure but you may also need to run AoV.vector().apply('insert') at the end before the return statement

---

Thank you Nate and Tormod. I guess I was doing it the hard way;)

I could use FacetArea(mesh), but FEniCS was not happy because I was restricting it with +/- for the ds integrals, which I thought it's needed. For some reasons, I could restrict my AoV function with +/- and the results appeared to be correct. In any case, FEnicS is now happy again, thanks to your comments.

A side note: It appears that the AoV value computed with my function is slightly different than the one obtained with the simple CellVolume(mesh)/FacetArea(mesh).