Pointwise Dirichlet BC

Try removing 'on_boundary'.

Same thing: singular matrix.

I have also tried:

class PinPoint(SubDomain):
    ...

but with no success.

Do you have any working example where pointwise Dirichlet BC are imposed?

I tried in a simpler code than the one I am working with. The code below works:

from dolfin import *

# Create the mesh
xMin, yMin = 0.0, 0.0
xMax, yMax = DOLFIN_PI, DOLFIN_PI

mesh = Rectangle(xMin, yMin, xMax, yMax, 32, 32)

# Define the function spaces
V = FunctionSpace(mesh, "CG", 1)

# Define the trial and test functions
u = TrialFunction(V)
v = TestFunction(V)

f = Expression("2*cos(x[0])*cos(x[1])")

# Define the variational forms
a = inner(grad(u), grad(v))*dx
L = f*v*dx

def pinpoint(x):
    return abs(x[0] - xMin) < DOLFIN_EPS and abs(x[1] - yMin) < DOLFIN_EPS

bc = DirichletBC(V, Constant(1.0), pinpoint, "pointwise")

problem = VariationalProblem(a, L, bc)

u = problem.solve()

plot(u)
interactive()
# ============================

But I am using mixed finite element and I am having a bad time to set pointwise Dirichlet BC.

- Allan

So, when I use mixed finite element with BDM and DG function space, is that correct to impose Dirichlet BC pointwise for the scalar field in DG? Because the value is not on a point anymore, but in a cell.

Is there a different procedure in this case?

If you run the code below for a pure Neumann problem with Mixed FEM, you you see that the gradient field sigma can be "determined", but the scalar field u assumes exorbitant high numbers, probably because the pointwise Dirichlet BC is not being processed.

from dolfin import *

# Create the mesh
xMin, yMin = 0.0, 0.0
xMax, yMax = DOLFIN_PI, DOLFIN_PI

mesh = Rectangle(xMin, yMin, xMax, yMax, 32, 32)

# Define function spaces and mixed (product) space
BDM = FunctionSpace(mesh, "BDM", 1)
DG = FunctionSpace(mesh, "DG", 0)
W = BDM * DG

# Define trial and test functions
(sigma, u) = TrialFunctions(W)
(tau, v) = TestFunctions(W)

# Define source function
f = Expression("2*cos(x[0])*cos(x[1])")

# Define variational form
a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx
L = - f*v*dx

# Define function u on the boundary
ub = Expression(("-sin(x[0])*cos(x[1])", "-cos(x[0])*sin(x[1])"))

# Define p on corner
pc = Constant(1.0)

# Define essential boundary
def boundary(x):
    return abs(x[0] - xMin) < DOLFIN_EPS or \
           abs(x[0] - xMax) < DOLFIN_EPS or \
           abs(x[1] - yMin) < DOLFIN_EPS or \
           abs(x[1] - yMax) < DOLFIN_EPS

def pinpoint(x):
    return abs(x[0] - xMin) < DOLFIN_EPS and abs(x[1] - yMin) < DOLFIN_EPS

bc = [DirichletBC(W.sub(0), ub, boundary), \
      DirichletBC(W.sub(1), pc, pinpoint, "pointwise")]

# Compute solution
w = Function(W)

A = assemble(a)
b = assemble(L)

[condition.apply(A, b) for condition in bc]

solve(A, w.vector(), b)

(sigma, u) = w.split()

# Plot sigma and u
plot(sigma)
plot(u)
interactive()

Applying an essential boundary condition for the scalar field
in a mixed Poisson formulation is very seldomly the right thing to do.

The issue at hand is that the scalar field is only determined
up to a constant when only Neumann conditions are applied. You can
remedy this by adding some additional constraint, for instance
that the average value (integral over the domain) of u is zero.
Here is a code snippet that does this (minor modification of
demo/pde/mixed-poisson/python/demo_mixed-poisson.py):

------
from dolfin import *

# Create mesh
mesh = UnitSquare(32, 32)

# Define function spaces and mixed (product) space
BDM = FunctionSpace(mesh, "BDM", 1)
DG = FunctionSpace(mesh, "DG", 0)
R = FunctionSpace(mesh, "R", 0)
W = MixedFunctionSpace([BDM, DG, R])

# Define trial and test functions
(sigma, u, c) = TrialFunctions(W)
(tau, v, d) = TestFunctions(W)

# Define source function
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")

# Define variational form
a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v + c*v + d*u)*dx
L = - f*v*dx

bc = DirichletBC(W.sub(0), (0.0, 0.0), "on_boundary")

# Compute solution
w = Function(W)
solve(a == L, w, bc)
(sigma, u, c) = w.split()
----

--
Marie

On 08/27/11 19:41, Allan Leal wrote:
> Allan Moreira Mulin Leal gave more information on the question:
> So, when I use mixed finite element with BDM and DG function space, is
> that correct to impose Dirichlet BC pointwise for the scalar field in
> DG? Because the value is not on a point anymore, but in a cell.
>
> Is there a different procedure in this case?
>
> If you run the code below for a pure Neumann problem with Mixed FEM, you
> you see that the gradient field sigma can be "determined", but the
> scalar field u assumes exorbitant high numbers, probably because the
> pointwise Dirichlet BC is not being processed.
>
> from dolfin import *
>
> # Create the mesh =
>
> xMin, yMin =3D 0.0, 0.0
> xMax, yMax =3D DOLFIN_PI, DOLFIN_PI
>
> mesh =3D Rectangle(xMin, yMin, xMax, yMax, 32, 32)
>
> # Define function spaces and mixed (product) space
> BDM =3D FunctionSpace(mesh, "BDM", 1)
> DG =3D FunctionSpace(mesh, "DG", 0)
> W =3D BDM * DG
>
> # Define trial and test functions
> (sigma, u) =3D TrialFunctions(W)
> (tau, v) =3D TestFunctions(W)
>
> # Define source function
> f =3D Expression("2*cos(x[0])*cos(x[1])")
>
> # Define variational form
> a =3D (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx
> L =3D - f*v*dx
>
> # Define function u on the boundary
> ub =3D Expression(("-sin(x[0])*cos(x[1])", "-cos(x[0])*sin(x[1])"))
>
> # Define p on corner
> pc =3D Constant(1.0)
>
> # Define essential boundary
> def boundary(x):
> return abs(x[0] - xMin)< DOLFIN_EPS or \
> abs(x[0] - xMax)< DOLFIN_EPS or \
> abs(x[1] - yMin)< DOLFIN_EPS or \
> abs(x[1] - yMax)< DOLFIN_EPS
>
> def pinpoint(x):
> return abs(x[0] - xMin)< DOLFIN_EPS and abs(x[1] - yMin)< DOLFIN_EPS
>
> bc =3D [DirichletBC(W.sub(0), ub, boundary), \
> DirichletBC(W.sub(1), pc, pinpoint, "pointwise")]
>
> # Compute solution
> w =3D Function(W)
>
> A =3D assemble(a)
> b =3D assemble(L)
>
> [condition.apply(A, b) for condition in bc]
>
> solve(A, w.vector(), b)
>
> (sigma, u) =3D w.split()
>
> # Plot sigmaDelivered-To: <email address hidden>
> Received: by 10.231.37.13 with SMTP id v13cs198772ibd;
> Sat, 27 Aug 2011 10:41:04 -0700 (PDT)
> Received: by 10.227.201.144 with SMTP id fa16mr2127753wbb.70.1314466863299;
> Sat, 27 Aug 2011 10:41:03 -0700 (PDT)
> Return-Path:<email address hidden>
> Received: from indium.canonical.com (indium.canonical.com [91.189.90.7])
> by mx.google.com with ESMTPS id et9si7131638wbb.44.2011.08.27.10.41.02
> (version=TLSv1/SSLv3 cipher=OTHER);
> Sat, 27 Aug 2011 10:41:03 -0700 (PDT)
> Received-SPF: pass (google.com: best guess record for domain of <email address hidden> designates 91.189.90.7 as permitted sender) client-ip=91.189.90.7;
> Authentication-Results: mx.google.com; spf=pass (google.com: best guess record for domain of <email address hidden> designates 91.189.90.7 as permitted sender) <email address hidden>
> Received: from loganberry.canonical.com ([91.189.90.37])
> by indium.canonical.com with esmtp (Exim 4.71 #1 (Debian))
> id 1QxMsP-0007LM-OE
> for<email address hidden>; Sat, 27 Aug 2011 17:41:01 +0000
> Received: from loganberry.canonical.com (localhost [127.0.0.1])
> by loganberry.canonical.com (Postfix) with ESMTP id 7D4952E810F
> for<email address hidden>; Sat, 27 Aug 2011 17:41:00 +0000 (UTC)
> Content-Type: text/plain; charset="utf-8"
> MIME-Version: 1.0
> Content-Transfer-Encoding: quoted-printable
> X-Launchpad-Question: product=dolfin; status=Open; assignee=None;
> priority=Normal; language=en
> X-Launchpad-Message-Rationale: Answer Contact (dolfin) @dolfin
> Reply-To: <email address hidden>
> References:<email address hidden>
> To: <email address hidden>
> From: Allan Leal<email address hidden>
> Subject: Re: [Question #169325]: Pointwise Dirichlet BC
> Message-Id:<email address hidden>
> Date: Sat, 27 Aug 2011 17:41:00 -0000
> Sender: <email address hidden>
> Errors-To: <email address hidden>
> Precedence: bulk
> X-Generated-By: Launchpad (canonical.com); Revision="13794";
> Instance="initZopeless config overlay"
> X-Launchpad-Hash: e1b0d24c821071565c8dc97e9d0ed4f5cde97dc8
>
> Question #169325 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/169325
>
> Allan Moreira Mulin Leal gave more information on the question:
> So, when I use mixed finite element with BDM and DG function space, is
> that correct to impose Dirichlet BC pointwise for the scalar field in
> DG? Because the value is not on a point anymore, but in a cell.
>
> Is there a different procedure in this case?
>
> If you run the code below for a pure Neumann problem with Mixed FEM, you
> you see that the gradient field sigma can be "determined", but the
> scalar field u assumes exorbitant high numbers, probably because the
> pointwise Dirichlet BC is not being processed.
>
> from dolfin import *
>
> # Create the mesh =
>
> xMin, yMin =3D 0.0, 0.0
> xMax, yMax =3D DOLFIN_PI, DOLFIN_PI
>
> mesh =3D Rectangle(xMin, yMin, xMax, yMax, 32, 32)
>
> # Define function spaces and mixed (product) space
> BDM =3D FunctionSpace(mesh, "BDM", 1)
> DG =3D FunctionSpace(mesh, "DG", 0)
> W =3D BDM * DG
>
> # Define trial and test functions
> (sigma, u) =3D TrialFunctions(W)
> (tau, v) =3D TestFunctions(W)
>
> # Define source function
> f =3D Expression("2*cos(x[0])*cos(x[1])")
>
> # Define variational form
> a =3D (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx
> L =3D - f*v*dx
>
> # Define function u on the boundary
> ub =3D Expression(("-sin(x[0])*cos(x[1])", "-cos(x[0])*sin(x[1])"))
>
> # Define p on corner
> pc =3D Constant(1.0)
>
> # Define essential boundary
> def boundary(x):
> return abs(x[0] - xMin)< DOLFIN_EPS or \
> abs(x[0] - xMax)< DOLFIN_EPS or \
> abs(x[1] - yMin)< DOLFIN_EPS or \
> abs(x[1] - yMax)< DOLFIN_EPS
>
> def pinpoint(x):
> return abs(x[0] - xMin)< DOLFIN_EPS and abs(x[1] - yMin)< DOLFIN_EPS
>
> bc =3D [DirichletBC(W.sub(0), ub, boundary), \
> DirichletBC(W.sub(1), pc, pinpoint, "pointwise")]
>
> # Compute solution
> w =3D Function(W)
>
> A =3D assemble(a)
> b =3D assemble(L)
>
> [condition.apply(A, b) for condition in bc]
>
> solve(A, w.vector(), b)
>
> (sigma, u) =3D w.split()
>
> # Plot sigma and u
> plot(sigma)
> plot(u)
> interactive()
>
> -- =
>
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

Thanks Marie Rognes, that solved my question.

Hi Marie,

I have already done something quite similar to what you wrote. I just took the demo for pure Neumann problem that you wrote (demo_neumann-poisson.py) and adapted to the mixed finite element case.

At first, I was having problem with the notation:

W = BDM * DG * R

Then, I found that the MixedFunctionSpace could be constructed as:

W = MixedFunctionSpace([BDM, DG, R])

So, now I have another question: why the operator overloading '*' does not work for more than two function spaces?

Try:

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

W = V * V * V

(u, v, w) = TrialFunctions(W)

- Allan

Thanks Kristian B. Ølgaard, that solved my question.