# How do I get the area of internal facets for DG methods?

I'm trying to implement a DG scheme that has a penalty term

h_F = min_{T \in (T+,T-)} |T| / |dT|

I've tried to implement this in PyDOLFIN as:

h_T = CellSize(mesh)
h_dT = FacetArea(mesh)
h_p = h_T('+') / h_dT('+')
h_m = h_T('-') / h_dT('-')
# h_F = min(h_p, h_m) = h_p < h_m ? h_p: h_m
h_F = conditional(lt(h_p, h_m), h_p, h_m)

but FacetArea is only defined for exterior facets.
Is there something corresponding for interior facets?

## Question information

Language:
English Edit question
Status:
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Last query:
 Revision history for this message Kristian B. Ølgaard (k.b.oelgaard) said on 2011-09-29: #1

On 29 September 2011 12:11, Martin Sandve Alnæs
> New question #172688 on DOLFIN:
>
> I'm trying to implement a DG scheme that has a penalty term
>
>  h_F = min_{T \in (T+,T-)}  |T| / |dT|
>
> I've tried to implement this in PyDOLFIN as:
>
>  h_T = CellSize(mesh)
>  h_dT = FacetArea(mesh)
>  h_p = h_T('+') / h_dT('+')
>  h_m = h_T('-') / h_dT('-')
>  # h_F = min(h_p, h_m) = h_p < h_m ? h_p: h_m
>  h_F = conditional(lt(h_p, h_m), h_p, h_m)
>
> but FacetArea is only defined for exterior facets.
> Is there something corresponding for interior facets?

I don't understand.

The form file:

element = FiniteElement("Discontinuous Lagrange", "triangle", 1)
u = TrialFunction(element)
v = TestFunction(element)
h = Constant(triangle)
h_avg = (h('+') + h('-'))/2
a = h_avg*u('+')*v('-')*dS

and main.cpp file:

// Create mesh
UnitSquare mesh(2, 2);

// Create functions
FacetArea h(mesh);

// Create funtion space
Poisson::FunctionSpace V(mesh);

// Define variational problem
Poisson::BilinearForm a(V, V);
a.h = h;

Matrix A;
assemble(A, a);

info(A, true);

works fine here.

Is there a typo in SpecialFunctions.cpp FacetArea::eval():

dolfin_error("SpecialFunctions.cpp",
"evaluate FacetArea expression",
"Facet area is only defined on mesh boundaries");

should it be 'cell boundaries' instead?

FacetArea only returns the area of the facet currently being
integrated which means that

h_dT = FacetArea(mesh)
h_dT('+') = h_dT('-')

Perhaps what you want is the sum of facet areas of a cell?
If that is the case, we can attach this ('surface_area'?) to Cell in
UFL like, cell.volume, cell.n, etc.

Kristian

> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>