Subject: Re: [FFC-dev] Test functions which are zero on the Dirichlet boundary?

From: Chong Luo <[email-masked]>
Date: Sat, 7 Jun 2008 17:24:24 -0700 (PDT)

Thank you for your clarification. I think Dirichlet boundary won't be a problem as I thought.

I'm glad to see that there's a Mixed Poisson demo in FEniCS. However I found that it cannot treat Neumann boundary conditions. So I still stick to the mixed formulation that I proposed in the first email. I seemed to have proved that for b(q,v) = (grad(q), v), choosing the space of q to be M=P(k+1) and the space of v to be X=P(k)^2 satisfies the LBB condition. However, this choice is unstable in FEniCS, while choosing M=P(k) and X to be P(k)^2 gives stability in FEniCS (at least in all the model problems I tested). Is there any reason for that? (Sorry, maybe this is more a numerical analysis question than a software question.)

Best,
Chong Luo

----- Original Message ----
From: Anders Logg <[email-masked]>
To: [email-masked]
Sent: Saturday, June 7, 2008 3:19:33 AM
Subject: Re: [FFC-dev] Test functions which are zero on the Dirichlet boundary?

On Fri, Jun 06, 2008 at 07:15:33PM -0700, Chong Luo wrote:
> Is it possible to specify test functions which are zero on the Dirichlet
> boundary?

No, you have to apply such constraints afterwards to the linear
system, using the DirichletBC class in DOLFIN.

> For example, let's say we want to use mixed finite element to solve Poisson's
> equation with both Neumann and Dirichlet boundary conditions:
> - Laplace (p) = f  in Omega
> p = pD on Gamma0
> dp/dn = uN on Gamma1
> 
> So we introduce velocity u = - grad(p), and get the following mixed
> formulation:
> (u, v) + (v, grad(p)) = 0                           for all v in X
> (u, grad(q))  = -(f,q) + <q, uN>              for all q in M
> where (,) is integration on Omega, while <,> is integration on the boundary of
> Omega.
> 
> I found that I need to take test function q such that q is zero on the
> Dirichlet boundary Gamma0 (which has to be nonempty), for this mixed
> formulation to satisfy LBB condition.  Is it possible to specify this
> constraint in FFC/FEniCS?
> 
> Thank you!

There is a demo that demonstrates how to implement the mixed Poisson
system in DOLFIN. Take a look in

  demo/pde/mixed-poisson/

This demo uses BDM1 elements for the flux and DG0 elements for the
pressure to get stability.

-- 
Anders

---

From: Marie Rognes <[email-masked]>
Date: Mon, 09 Jun 2008 11:47:17 +0200

Chong Luo wrote:
> I'm glad to see that there's a Mixed Poisson demo in FEniCS. However I 
> found that it cannot treat Neumann boundary conditions.

This demo uses the H(div) x L^2 formulation of the Mixed Poisson 
problem. With this formulation and *your* notation, p = 0 is a natural 
boundary condition, while conditions of the form u * n = g are essential 
boundary conditions.

The demo as it is assumes pure natural boundary conditions, but it is 
easy to extend. Again with your notation, simply replace

    a = (dot(u, v) - div(v)*p + q*div(u))*dx

by
     
    a = (dot(u, v) - div(v)*p + q*div(u))*dx + dot(v, n)*p*ds

where for instance
   
    n = FacetNormal("triangle", mesh)

and add the condition v*n = g using the DirichletBC.

> So I still stick to the mixed formulation that I proposed in the first 
> email. I seemed to have proved that for b(q,v) = (grad(q), v), 
> choosing the space of q to be M=P(k+1) and the space of v to be 
> X=P(k)^2 satisfies the LBB condition.

Continuous polynomials of degree k+1 for M and *discontinuous* 
polynomials of degree k for X, satisfy the stability conditions for the  
L^2 x H^1 mixed formulation.

> However, this choice is unstable in FEniCS,

If you were talking about continuous polynomials of degree k+1 for M and 
continuous polynomials of degree k for X, I think this choice would be 
unstable independently of software ;)

> while choosing M=P(k) and X to be P(k)^2 gives stability in FEniCS (at 
> least in all the model problems I tested). Is there any reason for 
> that? (Sorry, maybe this is more a numerical analysis question than a 
> software question.)
>

Really? With M = P(1) and X = P(1)^2 (both continuous), I get rather 
clear oscillations for the DOLFIN mixed-poisson demo problem using the 
L^2 - H^1 formulation, p = 0 on the x = 0 boundary and u * n = 0 on the 
rest.

--
Marie

> ------------------------------------------------------------------------
>
> _______________________________________________
> FFC-dev mailing list
> [email-masked]
> http://www.fenics.org/mailman/listinfo/ffc-dev
>   

-- 
Marie E. Rognes
Ph.D Fellow, 
Centre of Mathematics for Applications, 
University of Oslo
http://folk.uio.no/meg