ffc team mailing list archive

["Garth N. Wells" <gnw20@xxxxxxxxx>]

Kent Andre wrote:
> > Not sure what we have to do actually. I suppose I want get any clue from the 
> >
> >   demo/function/restriction
> >
> > demo? 
> >
> > I see that Kent has modified the dolfin.FunctionSpace to reflect the present 
> > (broken) implementation of restriction, anything in that line? Should such an 
> > implementation has been done in FunctionSpaceBase instead of FunctionSpace? 
> > Kent?
> >
> > Johan
>
> This restriction is different in the sense that the restriction is made
> by the standard element on a sub mesh. 
>
> The restriction in focus now is on the complete mesh but for only
> a part of the element. 
>
> Combined, these two features will be very powerful.  I guess the 
> sub-mesh restriction stuff is broken due to the parallel mesh. 


I've attached an example solver. To make it work for P^k with k > 2 we 
need to restrict some functions to cell facets. Using FFC, we do

    scalar  = FiniteElement("Lagrange", "triangle", 3)
    dscalar = FiniteElement("Discontinuous Lagrange", "triangle", 3)

    # This syntax restricts the space 'scalar' to cell facets.
    mixed = scalar["facet"] + dscalar

In PyDOLFIN we could do

    V_dg = FunctionSpace(mesh, "DG", 3)
    V_cg = FunctionSpace(mesh, "CG", 3, "facet")

    mixed = V_cg + V_dg

or

    V_dg = FunctionSpace(mesh, "DG", 3)
    V_cg = FunctionSpace(mesh, "CG", 3)

    mixed = V_cg["scalar"] + V_dg

Garth

> Kent
>
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev
""" Steady state advection-diffusion equation based on Labeur, R.J. 
and Wells, G.N. (2007), http://dx.doi.org/10.1016/j.cma.2007.06.025

"""

__author__ = "Garth N. Wells (gnw20@xxxxxxxxx) and Robert Jan Labeur (r.j.labeur@xxxxxxxxxx)"
__date__ = "2009-08-19"
__copyright__ = "Copyright (C) 2009 Garth N. Wells and Robert Jan Labeur"
__license__  = "GNU LGPL Version 2.1"

from dolfin import *

class Velocity(Function):
    def eval(self, values, x):
        values[0] = 0.0;
        values[1] = 0.0;
        #values[0] = -1.0;
        #values[1] = -0.4;

class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - 1.0) < DOLFIN_EPS and on_boundary

# Load mesh
#mesh = Mesh("unstruct-trimesh-2.xml.gz")
mesh = UnitSquare(32, 32)

# Defining the function spaces
V_dg = FunctionSpace(mesh, "DG", 2)
V_cg = FunctionSpace(mesh, "CG", 2)
V_b  = VectorFunctionSpace(mesh, "CG", 1)

mixed = V_cg + V_dg

# Test and trial functions
vhat, v = TestFunction(mixed)
uhat, u = TrialFunction(mixed)

# Set the the method (-1 = symmetric IP, 1 = Baumann-Oden)
method = -1.0

# Diffusivity
kappa = 1.0

# Source term
#f = 1.0
f = Function(V_cg, "500.0 * exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")

b = Velocity(V_b)

# Penalty term
alpha = 0.0e0

# Mesh-related functions
n = FacetNormal(mesh)
h = CellSize(mesh)

# ( dot(v, n) + |dot(v, n)| )/2.0  = an on outflow,  0 on inflow 
an = (dot(b, n) + abs(dot(b, n)))/2.0

# ( -dot(v, n) + |dot(v, n)| )/2.0  = 0 on outflow,  -an on inflow 
an_min = (-dot(b, n) + abs(dot(b, n)))/2.0

# Standard flux
sigma = -b*u + kappa*grad(u) 

# Interface flux
sigma_a = -dot(b, n)*u + an_min*(uhat - u)
sigma_d = kappa*grad(u) + (alpha/h)*kappa*(uhat -u)*n
sigma_hat = sigma_a + dot(sigma_d, n)

# Bilinear form
a_local = dot(grad(v), sigma)*dx \
         - v('+')*sigma_hat('+')*dS - v('-')*sigma_hat('-')*dS \
         - method*dot(grad(v)('+'), n('+'))*kappa*(uhat('+') - u('+'))*dS \
         - method*dot(grad(v)('-'), n('-'))*kappa*(uhat('-') - u('-'))*dS \
         - v*sigma_hat*ds \
         - method*dot(grad(v), n)*kappa*(uhat - u)*ds

a_global = vhat('+')*sigma_hat('+')*dS + vhat('-')*sigma_hat('-')*dS \
         + vhat*sigma_hat*ds \
         + vhat*an*uhat*ds  # active on outflow boundary
          #+ vhat*hbc*ds  # Flux bc

a = a_local + a_global

# Linear form
L = v*f*dx

# Set up boundary condition (apply strong BCs)
g = Function(V_dg,"sin(pi*5.0*x[1])")
bc = DirichletBC(V_cg, g, DirichletBoundary())

# Solution function
uh = Function(mixed)

# Assemble and apply boundary conditions
A = assemble(a)
b = assemble(L)
bc.apply(A, b)

# Solve system
solve(A, uh.vector(), b)

file = File("u.pvd")
file << uh.sub(0)

# Plot solution
plot(uh.sub(0), interactive=True)