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

Asked by Martin Sandve Alnæs

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:
Answered
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Kristian B. Ølgaard (k.b.oelgaard) said :
#1

On 29 September 2011 12:11, Martin Sandve Alnæs
<email address hidden> wrote:
> New question #172688 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/172688
>
> 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('-')

in your case.
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.
>

Can you help with this problem?

Provide an answer of your own, or ask Martin Sandve Alnæs for more information if necessary.

To post a message you must log in.