evaluation of normal derivative edge jumps

Hi

It should work though I have not tried a full example. But following code seems to go through

from dolfin import *

mesh=UnitSquare(10,10)

V=FunctionSpace(mesh, "CG", 1)

u = Function(V)
u = interpolate(Expression("x[0]*x[1]"), V)

n = FacetNormal(mesh)
un= dot(grad(u), n)
B = jump(un)*dS
print assemble(B)

Thanks for the quick answer! This seems to work. However, to enable adaptivity, the integral on each edge (instead of the sum over all edges as computed above) is required. How can that be done?

Best,
Martin

On Tue, Nov 01, 2011 at 04:56:25PM -0000, Martin Eigel wrote:
> Question #177108 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/177108
>
> Status: Answered => Open
>
> Martin Eigel is still having a problem:
> Thanks for the quick answer! This seems to work. However, to enable
> adaptivity, the integral on each edge (instead of the sum over all edges
> as computed above) is required. How can that be done?

I only have time for a very quick response: multiply the integral with
a piecewise constant test function and assemble into a vector. Then
each vector entry will contain an indicator for a cell (if the
integral is chosen correctly).

--
Anders

Thanks again. I am afraid I am not experienced enough with dolfin to figure out how to do this. Could someone please elaborate on how to assemble the vector of edge integrals?

Martin

On Wednesday November 2 2011 03:15:44 Martin Eigel wrote:
> Question #177108 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/177108
>
> Status: Answered => Open
>
> Martin Eigel is still having a problem:
> Thanks again. I am afraid I am not experienced enough with dolfin to
> figure out how to do this. Could someone please elaborate on how to
> assemble the vector of edge integrals?

I just modify the code from Praveen here:

from dolfin import *

mesh=UnitSquare(10,10)

V = FunctionSpace(mesh, "CG", 1)
Vdg = FunctionSpace(mesh, "DG", 0)

v = TestFunction(Vdg)
u = interpolate(Expression("x[0]*x[1]"), V)

n = FacetNormal(mesh)
un= dot(grad(u), n)
B = jump(un)*v*dS
b = assemble(B)

Johan

Thanks. This results in an "UFLException: Form argument must be restricted" when assemble(B) is called.
Btw, I am using dolfin 1.0 beta2.

Martin

Try

B = jump(un)*v('+')*dS

On 3 November 2011 06:40, Praveen C
<email address hidden> wrote:
> Question #177108 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/177108
>
>    Status: Open => Answered
>
> Praveen C proposed the following answer:
> Try
>
> B = jump(un)*v('+')*dS

This should compile and run, the question is if the output is as expected?
Since 'v' is defined on a DG0 element, which is constant for each
cell, I expect the result will be the sum of integrals of all adjacent
facets. If the desired output is the integral jump(un) on each
individual facet, then another approach must be taken. Not sure how
though.

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

For residual error estimate, what is required is something like "integral of jump(un) on interior boundary of element". So the sum would be correct in that case.

Yes, one indeed needs the integral of the sum of the two "adjacent functions" (in 2d) along an edge.
Just a couple follow-up questions:
a) I have not come across the syntax v('+') before but now did see it in some undocumented dG-examples. Is it documented anywhere? I am new to dolfin so apologies if this can easily be found in the documentation.
b) In the sample code above, b.array().size is 200, mesh.num_edges() is 320 whereas I'd expect them to be same size. Is this due to boundary edges?
c) Depending on boundary conditions, boundary edges would also have to be evaluated (as jump with regard to prescribed boundary data). How can this be accomplished? As a matter of fact, there should be an integral being evaluated on all edges of the mesh.

Martin

(a) For an edge e, joining cells T+ and T-, v('+') and v('-') refers to restriction to T+ and T- respectively.

(b) The size of b corresponds to the number of cells because it corresponds to the dimension of the function space over which you have defined your test function, in this case DG0.

(c) I'm copy-pasting a complete example for standard residual-based adaptivity for Poisson below. It might not be exactly what you want, but it could be a decent starting point.

__author__ = "Marie E. Rognes (<email address hidden>)"
__date__ = "2009-11-20 -- 2011-11-03"

"""
Demo for adaptive Poisson using a residual-based energy-norm error
estimator

  eta_h**2 = sum_T eta_T**2

with

  eta_T**2 = h_T**2 ||R_T||_T**2 + h_T ||R_dT||_dT**2

where

  R_T = - (f + div grad u_h)
  R_dT = jump(grad u_h * n)

and a 'maximal' marking strategy (refining those cells for
which the error indicator is greater than a certain fraction of the
largest error indicator)
"""

from dolfin import *
from numpy import array, sqrt

# Define Dirichlet boundary (x = 0 or x = 1)
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS

# Error tolerance
tolerance = 0.1
max_iterations = 20

# Refinement fraction
fraction = 0.7

# Create initial mesh
mesh = UnitSquare(4, 4)

for i in range(max_iterations):

    # Define variational problem
    V = FunctionSpace(mesh, "CG", 1)
    v = TestFunction(V)
    u = TrialFunction(V)
    f = Expression("10*exp(-(pow(x[0] - 0.6, 2) + pow(x[1] - 0.4, 2)) / 0.02)",
                   degree=3)
    a = inner(grad(v), grad(u))*dx
    L = v*f*dx

    # Define boundary condition
    bc = DirichletBC(V, 0.0, "near(x[0], 0.0) || near(x[0], 1.0)")

    # Compute solution
    u_h = Function(V)
    solve(a == L, u_h, bc)

    # Define cell and facet residuals
    R_T = - (f + div(grad(u_h)))
    n = FacetNormal(mesh)
    R_dT = dot(grad(u_h), n)

    # Will use space of constants to localize indicator form
    Constants = FunctionSpace(mesh, "DG", 0)
    w = TestFunction(Constants)
    h = CellSize(mesh)

    # Define form for assembling error indicators
    form = (h**2*R_T**2*w*dx + avg(h)*avg(R_dT)**2*2*avg(w)*dS
            + h*R_dT**2*w*ds)

    # Assemble error indicators
    indicators = assemble(form)

    # Calculate error
    error_estimate = sqrt(sum(i for i in indicators.array()))
    print "error_estimate = ", error_estimate

    # Take sqrt of indicators
    indicators = array([sqrt(i) for i in indicators])

    # Stop if error estimates is less than tolerance
    if error_estimate < tolerance:
        print "\nEstimated error < tolerance: %s < %s" % (error_estimate,
                                                        tolerance)
        break

    # Mark cells for refinement based on maximal marking strategy
    largest_error = max(indicators)
    cell_markers = MeshFunction("bool", mesh, mesh.topology().dim())
    for c in cells(mesh):
        cell_markers[c] = indicators[c.index()] > (fraction*largest_error)

    # Refine mesh
    mesh = refine(mesh, cell_markers)

# Plot final solution
plot(u_h)
plot(mesh)

# Hold plots
interactive()

re b) Ok, of course... Just as a remark, since edge jumps dominate the error in case of low-order FEM (and volume terms can be neglected), one would/could prefer marking edges instead of cells. Thus, the evaluation of an edge indicator might be advisable. If there is a direct way to achieve this, please advise.

re c) Splendid! Will try to proceed from there. Would be a nice addition to the documented examples ;)

Thanks to all for the support,
Martin