Setting 'pressure' BC in RT/P0 elements for mixed poisson
I want to use RT or BDM elements in a mixedpoisson type approach to solve Darcy flow.
Because the "pressure" is discretised on P0 elements, I cannot define it on the boundary, so I have included a
surface source term in the RHS. Is this the correct approach?
Unfortunately, although it works on the left hand end of my domain, the right hand end gives a slight discontinuity in slope.
Any help much appreciated,
Chris
mesh = UnitSquare(20,20)
n= FacetNormal(mesh)
R = FunctionSpace(
W = FunctionSpace(
V = R+W
(u,p) = TrialFunctions(V)
(v,q) = TestFunctions(V)
# boundary conditions
def boundary0(x):
return x[1]<DOLFIN_EPS
def boundary1(x):
return x[1]>1.0DOLFIN_EPS
u0 = Constant((0.0,0.0))
bc0 = DirichletBC(
bc1 = DirichletBC(
bcs=[bc0,bc1]
f=Expression(
# 'Permeability' k  gaussian in y
k=Expression(
a = (k*dot(v,u)  div(v)*p + q*div(u))*dx
L = dot(v,f*n)*ds
# Compute solution
problem = VariationalProb
(u, p) = problem.
# Plot solution
plot(p, interactive=True)
plot(u, interactive=True)
Question information
 Language:
 English Edit question
 Status:
 Solved
 For:
 DOLFIN Edit question
 Assignee:
 No assignee Edit question
 Solved by:
 Chris Richardson
 Solved:
 20100401
 Last query:
 20100401
 Last reply:
 20100331
Marie Rognes (megsimula) said :  #1 
Chris Richardson wrote:
> New question #105900 on DOLFIN:
> https:/
>
> I want to use RT or BDM elements in a mixedpoisson type approach to solve Darcy flow.
>
> Because the "pressure" is discretised on P0 elements, I cannot define it on the boundary, so I have included a
> surface source term in the RHS. Is this the correct approach?
>
>
Yes, this looks correct.
> Unfortunately, although it works on the left hand end of my domain, the right hand end gives a slight discontinuity in slope.
>
>
I suspect that this is simply an artifact associated with plotting DG_0
functions.
A similar "dip" appears with this simple script:
from dolfin import *
mesh = UnitSquare(20, 20)
V = FunctionSpace(
f = Expression(" 2.0*x[0] + 1.0", degree=1)
g = project(f, V)
plot(g, interactive=True)

Marie
> Any help much appreciated,
> Chris
>
>
> mesh = UnitSquare(20,20)
> n= FacetNormal(mesh)
>
> R = FunctionSpace(
> W = FunctionSpace(
> V = R+W
>
> (u,p) = TrialFunctions(V)
> (v,q) = TestFunctions(V)
>
> # boundary conditions
> def boundary0(x):
> return x[1]<DOLFIN_EPS
> def boundary1(x):
> return x[1]>1.0DOLFIN_EPS
>
> u0 = Constant((0.0,0.0))
> bc0 = DirichletBC(
> bc1 = DirichletBC(
> bcs=[bc0,bc1]
>
> f=Expression(
>
> # 'Permeability' k  gaussian in y
> k=Expression(
>
> a = (k*dot(v,u)  div(v)*p + q*div(u))*dx
> L = dot(v,f*n)*ds
>
> # Compute solution
> problem = VariationalProb
> (u, p) = problem.
>
> # Plot solution
> plot(p, interactive=True)
> plot(u, interactive=True)
>
>
>
Garth Wells (garthwells) said :  #2 
On 30/03/10 02:32, Marie Rognes wrote:
> Chris Richardson wrote:
>> New question #105900 on DOLFIN:
>> https:/
>>
>> I want to use RT or BDM elements in a mixedpoisson type approach to
>> solve Darcy flow.
>>
>> Because the "pressure" is discretised on P0 elements, I cannot define
>> it on the boundary, so I have included a surface source term in the
>> RHS. Is this the correct approach?
>
> Yes, this looks correct.
>
>> Unfortunately, although it works on the left hand end of my domain,
>> the right hand end gives a slight discontinuity in slope.
>>
>
>
> I suspect that this is simply an artifact associated with plotting DG_0
> functions.
Could well be. plot(..) performs a projection onto continuous P1
elements in the background. The VTK file output format supports
cellwise data (i.e. DG0), so try
File(
and view the result with ParaView.
Garth
> A similar "dip" appears with this simple script:
>
> from dolfin import *
> mesh = UnitSquare(20, 20)
> V = FunctionSpace(
> f = Expression(" 2.0*x[0] + 1.0", degree=1)
> g = project(f, V)
> plot(g, interactive=True)
>
>
>
> 
> Marie
>
>
>> Any help much appreciated,
>> Chris
>>
>>
>> mesh = UnitSquare(20,20)
>> n= FacetNormal(mesh)
>>
>> R = FunctionSpace(
>> W = FunctionSpace(
>> V = R+W
>> (u,p) = TrialFunctions(V)
>> (v,q) = TestFunctions(V)
>>
>> # boundary conditions
>> def boundary0(x):
>> return x[1]<DOLFIN_EPS def boundary1(x):
>> return x[1]>1.0DOLFIN_EPS
>>
>> u0 = Constant((0.0,0.0))
>> bc0 = DirichletBC(
>> bc1 = DirichletBC(
>> bcs=[bc0,bc1]
>>
>> f=Expression(
>>
>> # 'Permeability' k  gaussian in y
>> k=Expression(
>>
>> a = (k*dot(v,u)  div(v)*p + q*div(u))*dx
>> L = dot(v,f*n)*ds
>>
>> # Compute solution
>> problem = VariationalProb
>> (u, p) = problem.
>>
>> # Plot solution
>> plot(p, interactive=True)
>> plot(u, interactive=True)
>>
>>
>
>
> _______
> Mailing list: https:/
> Post to : <email address hidden>
> Unsubscribe : https:/
> More help : https:/
Anders Logg (logg) said :  #3 
On Tue, Mar 30, 2010 at 07:21:25AM 0000, Garth Wells wrote:
> Question #105900 on DOLFIN changed:
> https:/
>
> Garth Wells proposed the following answer:
>
> On 30/03/10 02:32, Marie Rognes wrote:
> > Chris Richardson wrote:
> >> New question #105900 on DOLFIN:
> >> https:/
> >>
> >> I want to use RT or BDM elements in a mixedpoisson type approach to
> >> solve Darcy flow.
> >>
> >> Because the "pressure" is discretised on P0 elements, I cannot define
> >> it on the boundary, so I have included a surface source term in the
> >> RHS. Is this the correct approach?
> >
> > Yes, this looks correct.
> >
> >> Unfortunately, although it works on the left hand end of my domain,
> >> the right hand end gives a slight discontinuity in slope.
> >>
> >
> >
> > I suspect that this is simply an artifact associated with plotting DG_0
> > functions.
>
> Could well be. plot(..) performs a projection onto continuous P1
> elements in the background.
I think it does interpolation. Projection should be fine.

Anders
> The VTK file output format supports
> cellwise data (i.e. DG0), so try
>
> File("pressure.
>
> and view the result with ParaView.
>
> Garth
>
> > A similar "dip" appears with this simple script:
> >
> > from dolfin import *
> > mesh = UnitSquare(20, 20)
> > V = FunctionSpace(
> > f = Expression(" 2.0*x[0] + 1.0", degree=1)
> > g = project(f, V)
> > plot(g, interactive=True)
> >
> >
> >
> >
> >
> >> Any help much appreciated,
> >> Chris
> >>
> >>
> >> mesh = UnitSquare(20,20)
> >> n= FacetNormal(mesh)
> >>
> >> R = FunctionSpace(
> >> W = FunctionSpace(
> >> V = R+W
> >> (u,p) = TrialFunctions(V)
> >> (v,q) = TestFunctions(V)
> >>
> >> # boundary conditions
> >> def boundary0(x):
> >> return x[1]<DOLFIN_EPS def boundary1(x):
> >> return x[1]>1.0DOLFIN_EPS
> >>
> >> u0 = Constant((0.0,0.0))
> >> bc0 = DirichletBC(
> >> bc1 = DirichletBC(
> >> bcs=[bc0,bc1]
> >>
> >> f=Expression(
> >>
> >> # 'Permeability' k  gaussian in y
> >> k=Expression(
> >>
> >> a = (k*dot(v,u)  div(v)*p + q*div(u))*dx
> >> L = dot(v,f*n)*ds
> >>
> >> # Compute solution
> >> problem = VariationalProb
> >> (u, p) = problem.
> >>
> >> # Plot solution
> >> plot(p, interactive=True)
> >> plot(u, interactive=True)
> >>
> >>
> >
> >
> > _______
> > Mailing list: https:/
> > Post to : <email address hidden>
> > Unsubscribe : https:/
> > More help : https:/
>
Garth Wells (garthwells) said :  #4 
On 30/03/10 15:39, Anders Logg wrote:
> Question #105900 on DOLFIN changed:
> https:/
>
> Anders Logg proposed the following answer:
> On Tue, Mar 30, 2010 at 07:21:25AM 0000, Garth Wells wrote:
>> Question #105900 on DOLFIN changed:
>> https:/
>>
>> Garth Wells proposed the following answer:
>>
>> On 30/03/10 02:32, Marie Rognes wrote:
>>> Chris Richardson wrote:
>>>> New question #105900 on DOLFIN:
>>>> https:/
>>>>
>>>> I want to use RT or BDM elements in a mixedpoisson type approach to
>>>> solve Darcy flow.
>>>>
>>>> Because the "pressure" is discretised on P0 elements, I cannot define
>>>> it on the boundary, so I have included a surface source term in the
>>>> RHS. Is this the correct approach?
>>>
>>> Yes, this looks correct.
>>>
>>>> Unfortunately, although it works on the left hand end of my domain,
>>>> the right hand end gives a slight discontinuity in slope.
>>>>
>>>
>>>
>>> I suspect that this is simply an artifact associated with plotting DG_0
>>> functions.
>>
>> Could well be. plot(..) performs a projection onto continuous P1
>> elements in the background.
>
> I think it does interpolation. Projection should be fine.
>
Yes, I thought that just after sending the message.
Garth
> 
> Anders
>
>
>> The VTK file output format supports
>> cellwise data (i.e. DG0), so try
>>
>> File("pressure.
>>
>> and view the result with ParaView.
>>
>> Garth
>>
>>> A similar "dip" appears with this simple script:
>>>
>>> from dolfin import *
>>> mesh = UnitSquare(20, 20)
>>> V = FunctionSpace(
>>> f = Expression(" 2.0*x[0] + 1.0", degree=1)
>>> g = project(f, V)
>>> plot(g, interactive=True)
>>>
>>>
>>>
>>>
>>>
>>>> Any help much appreciated,
>>>> Chris
>>>>
>>>>
>>>> mesh = UnitSquare(20,20)
>>>> n= FacetNormal(mesh)
>>>>
>>>> R = FunctionSpace(
>>>> W = FunctionSpace(
>>>> V = R+W
>>>> (u,p) = TrialFunctions(V)
>>>> (v,q) = TestFunctions(V)
>>>>
>>>> # boundary conditions
>>>> def boundary0(x):
>>>> return x[1]<DOLFIN_EPS def boundary1(x):
>>>> return x[1]>1.0DOLFIN_EPS
>>>>
>>>> u0 = Constant((0.0,0.0))
>>>> bc0 = DirichletBC(
>>>> bc1 = DirichletBC(
>>>> bcs=[bc0,bc1]
>>>>
>>>> f=Expression(
>>>>
>>>> # 'Permeability' k  gaussian in y
>>>> k=Expression(
>>>>
>>>> a = (k*dot(v,u)  div(v)*p + q*div(u))*dx
>>>> L = dot(v,f*n)*ds
>>>>
>>>> # Compute solution
>>>> problem = VariationalProb
>>>> (u, p) = problem.
>>>>
>>>> # Plot solution
>>>> plot(p, interactive=True)
>>>> plot(u, interactive=True)
>>>>
>>>>
>>>
>>>
>>> _______
>>> Mailing list: https:/
>>> Post to : <email address hidden>
>>> Unsubscribe : https:/
>>> More help : https:/
>>
>
Chris Richardson (chrisbpi) said :  #5 
OK, thanks that clears up one thing...
My question now becomes how to set the flux BC on the other boundaries when they are not so simple and
straight.
e.g. I thought this would work, but is giving real problems at the corners:
mesh = UnitSquare(20,20)
for x in mesh.coordinates():
x[1] += 0.5*x[0]
n= FacetNormal(mesh)
R = FunctionSpace(
W = FunctionSpace(
V = R+W
(u,p) = TrialFunctions(V)
(v,q) = TestFunctions(V)
# boundary conditions
def boundary0(
return x[0]>DOLFIN_EPS and x[0]<1.0DOLFIN_EPS and on_bound
u0 = Constant((0.0,0.0))
bcs = DirichletBC(
f=Expression(
# f=Constant(0.0)
# 'Permeability' k
k=Expression(
a = (k*dot(v,u)  div(v)*p + q*div(u))*dx
L = dot(v,f*n)*ds
# Compute solution
problem = VariationalProb
(u, p) = problem.
# Plot solution
plot(p, interactive=True)
plot(u, interactive=True)
Marie Rognes (megsimula) said :  #6 
Chris Richardson wrote:
> Question #105900 on DOLFIN changed:
> https:/
>
> Status: Answered => Open
>
> Chris Richardson is still having a problem:
> OK, thanks that clears up one thing...
>
> My question now becomes how to set the flux BC on the other boundaries when they are not so simple and
> straight.
>
> e.g. I thought this would work, but is giving real problems at the
> corners:
>
>
I've attached a revised version. Take a look and see if it makes sense.
Note that with (the previous):
def boundary0(
(*) return x[0]>DOLFIN_EPS and x[0]<1.0DOLFIN_EPS and on_bound
bcs = DirichletBC(
only the facets where both endpoints (and the midpoint) satisfy (*) will
be marked. (Somebody correct me if I'm wrong here.) In particular, the
leftmost "bottom" edge will not be marked.
Further, then v*n will not vanish on this edge, and hence you will get
contributions from this edge in:
L = dot(v,f*n)*ds
Then, the definition of f on the edges starts mattering etc.
The moral is that one should be really careful with the boundary
conditions ;)

Marie
Chris Richardson (chrisbpi) said :  #7 
Thanks  I am still not sure if it is behaving 100% as expected, although maybe my expectations are wrong...
the velocity field does still look a bit odd at the corners.
I suppose I am a bit confused as to the meaning of "x" in the boundary conditions. Is it the midpoint of an edge? or does it have to apply to both ends of the edge?
Anders Logg (logg) said :  #8 
On Tue, Mar 30, 2010 at 02:01:01PM 0000, Chris Richardson wrote:
> Question #105900 on DOLFIN changed:
> https:/
>
> Chris Richardson posted a new comment:
> Thanks  I am still not sure if it is behaving 100% as expected, although maybe my expectations are wrong...
> the velocity field does still look a bit odd at the corners.
>
> I suppose I am a bit confused as to the meaning of "x" in the boundary
> conditions. Is it the midpoint of an edge? or does it have to apply to
> both ends of the edge?
It has to apply to all vertices + the midpoint. Only then is the facet
(edge) marked.
The + midpoint condition was added recently to avoid setting boundary
conditions in cases like this:

/
/ 
/ 
_____/___

Anders
Marie Rognes (megsimula) said :  #9 
Chris Richardson wrote:
> Question #105900 on DOLFIN changed:
> https:/
>
> Chris Richardson posted a new comment:
> Thanks  I am still not sure if it is behaving 100% as expected, although maybe my expectations are wrong...
> the velocity field does still look a bit odd at the corners.
>
>
My best suggestion is to find yourself an exact test case and check
convergence.
If not, I can take a closer look.

Marie
Chris Richardson (chrisbpi) said :  #10 
I still don't really understand why you need the extra qualification in these statements. Surely, the boundary
would be equally well defined without the x[1] clauses. Why do I need to put in the x[1] range so exactly?
left = x[0] < DOLFIN_EPS and (x[1] > DOLFIN_EPS
right = abs(x[0]  1.0) < DOLFIN_EPS and (abs(x[1]  .5) > DOLFIN_EPS
return on_boundary and not (left or right)
Anders Logg (logg) said :  #11 
On Wed, Mar 31, 2010 at 09:21:50AM 0000, Chris Richardson wrote:
> Question #105900 on DOLFIN changed:
> https:/
>
> Status: Answered => Open
>
> Chris Richardson is still having a problem:
>
> I still don't really understand why you need the extra qualification in these statements. Surely, the boundary
> would be equally well defined without the x[1] clauses. Why do I need to put in the x[1] range so exactly?
>
>
> left = x[0] < DOLFIN_EPS and (x[1] > DOLFIN_EPS
> and x[1] < 1.0  DOLFIN_EPS)
If you don't include x[1] here, you will include the points (0, 0) and
(0, 1). That means those point will be in "left" and thus
return on_boundary and not (left or right)
will return false for those points. I assume you want to include those
points in a noslip boundary condition or similar.
> right = abs(x[0]  1.0) < DOLFIN_EPS and (abs(x[1]  .5) > DOLFIN_EPS
> and x[1] < 2.0  DOLFIN_EPS)
This looks strange. What does x[1] range mean here?

Anders
Chris Richardson (chrisbpi) said :  #12 
Thanks, I see the point now.
Is there any other way of setting markers? If I create a MeshFunction, can I access the coordinates of the
facets individually, and iterate through them?
Marie Rognes (megsimula) said :  #13 
Anders Logg wrote:
> On Wed, Mar 31, 2010 at 09:21:50AM 0000, Chris Richardson wrote:
>
>> Question #105900 on DOLFIN changed:
>> https:/
>>
>> Status: Answered => Open
>>
>> Chris Richardson is still having a problem:
>>
>> I still don't really understand why you need the extra qualification in these statements. Surely, the boundary
>> would be equally well defined without the x[1] clauses. Why do I need to put in the x[1] range so exactly?
>>
>>
>> left = x[0] < DOLFIN_EPS and (x[1] > DOLFIN_EPS
>> and x[1] < 1.0  DOLFIN_EPS)
>>
>
> If you don't include x[1] here, you will include the points (0, 0) and
> (0, 1). That means those point will be in "left" and thus
>
> return on_boundary and not (left or right)
>
> will return false for those points. I assume you want to include those
> points in a noslip boundary condition or similar.
>
>
>> right = abs(x[0]  1.0) < DOLFIN_EPS and (abs(x[1]  .5) > DOLFIN_EPS
>> and x[1] < 2.0  DOLFIN_EPS)
>>
>
> This looks strange. What does x[1] range mean here?
>
>
I don't know  I suggested
right = abs(x[0]  1.0) < DOLFIN_EPS and (abs(x[1]  .5) > DOLFIN_EPS

Marie
> 
> Anders
>
>
>
> 
>
> _______
> Mailing list: https:/
> Post to : <email address hidden>
> Unsubscribe : https:/
> More help : https:/
>
Chris Richardson (chrisbpi) said :  #14 
OK, I think I now understand. Thanks