Setting 'pressure' BC in RT/P0 elements for mixed poisson

Chris Richardson wrote:
> New question #105900 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/105900
>
> I want to use RT or BDM elements in a mixed-poisson 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(mesh,"DG",0)
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(mesh,"BDM",1)
> W = FunctionSpace(mesh,"DG",0)
> 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.0-DOLFIN_EPS
>
> u0 = Constant((0.0,0.0))
> bc0 = DirichletBC(V.sub(0), u0, boundary0)
> bc1 = DirichletBC(V.sub(0), u0, boundary1)
> bcs=[bc0,bc1]
>
> f=Expression("-1.0*(x[0]==0.0)+1.0*(x[0]==1.0)")
>
> # 'Permeability' k - gaussian in y
> k=Expression("1.0+exp(-40*pow(x[1]-0.5,2))")
>
> a = (k*dot(v,u) - div(v)*p + q*div(u))*dx
> L = dot(v,f*n)*ds
>
> # Compute solution
> problem = VariationalProblem(a, L,bcs)
> (u, p) = problem.solve().split()
>
> # Plot solution
> plot(p, interactive=True)
> plot(u, interactive=True)
>
>
>

On 30/03/10 02:32, Marie Rognes wrote:
> Chris Richardson wrote:
>> New question #105900 on DOLFIN:
>> https://answers.launchpad.net/dolfin/+question/105900
>>
>> I want to use RT or BDM elements in a mixed-poisson 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
cell-wise data (i.e. DG0), so try

     File("pressure.pvd") << p

and view the result with ParaView.

Garth

> A similar "dip" appears with this simple script:
>
> from dolfin import *
> mesh = UnitSquare(20, 20)
> V = FunctionSpace(mesh,"DG",0)
> 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(mesh,"BDM",1)
>> W = FunctionSpace(mesh,"DG",0)
>> 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.0-DOLFIN_EPS
>>
>> u0 = Constant((0.0,0.0))
>> bc0 = DirichletBC(V.sub(0), u0, boundary0)
>> bc1 = DirichletBC(V.sub(0), u0, boundary1)
>> bcs=[bc0,bc1]
>>
>> f=Expression("-1.0*(x[0]==0.0)+1.0*(x[0]==1.0)")
>>
>> # 'Permeability' k - gaussian in y
>> k=Expression("1.0+exp(-40*pow(x[1]-0.5,2))")
>>
>> a = (k*dot(v,u) - div(v)*p + q*div(u))*dx
>> L = dot(v,f*n)*ds
>>
>> # Compute solution
>> problem = VariationalProblem(a, L,bcs)
>> (u, p) = problem.solve().split()
>>
>> # Plot solution
>> plot(p, interactive=True)
>> plot(u, interactive=True)
>>
>>
>
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp

On Tue, Mar 30, 2010 at 07:21:25AM -0000, Garth Wells wrote:
> Question #105900 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/105900
>
> Garth Wells proposed the following answer:
>
> On 30/03/10 02:32, Marie Rognes wrote:
> > Chris Richardson wrote:
> >> New question #105900 on DOLFIN:
> >> https://answers.launchpad.net/dolfin/+question/105900
> >>
> >> I want to use RT or BDM elements in a mixed-poisson 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
> cell-wise data (i.e. DG0), so try
>
> File("pressure.pvd") << p
>
> and view the result with ParaView.
>
> Garth
>
> > A similar "dip" appears with this simple script:
> >
> > from dolfin import *
> > mesh = UnitSquare(20, 20)
> > V = FunctionSpace(mesh,"DG",0)
> > 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(mesh,"BDM",1)
> >> W = FunctionSpace(mesh,"DG",0)
> >> 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.0-DOLFIN_EPS
> >>
> >> u0 = Constant((0.0,0.0))
> >> bc0 = DirichletBC(V.sub(0), u0, boundary0)
> >> bc1 = DirichletBC(V.sub(0), u0, boundary1)
> >> bcs=[bc0,bc1]
> >>
> >> f=Expression("-1.0*(x[0]==0.0)+1.0*(x[0]==1.0)")
> >>
> >> # 'Permeability' k - gaussian in y
> >> k=Expression("1.0+exp(-40*pow(x[1]-0.5,2))")
> >>
> >> a = (k*dot(v,u) - div(v)*p + q*div(u))*dx
> >> L = dot(v,f*n)*ds
> >>
> >> # Compute solution
> >> problem = VariationalProblem(a, L,bcs)
> >> (u, p) = problem.solve().split()
> >>
> >> # Plot solution
> >> plot(p, interactive=True)
> >> plot(u, interactive=True)
> >>
> >>
> >
> >
> > _______________________________________________
> > Mailing list: https://launchpad.net/~dolfin
> > Post to : <email address hidden>
> > Unsubscribe : https://launchpad.net/~dolfin
> > More help : https://help.launchpad.net/ListHelp
>

On 30/03/10 15:39, Anders Logg wrote:
> Question #105900 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/105900
>
> 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://answers.launchpad.net/dolfin/+question/105900
>>
>> Garth Wells proposed the following answer:
>>
>> On 30/03/10 02:32, Marie Rognes wrote:
>>> Chris Richardson wrote:
>>>> New question #105900 on DOLFIN:
>>>> https://answers.launchpad.net/dolfin/+question/105900
>>>>
>>>> I want to use RT or BDM elements in a mixed-poisson 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
>> cell-wise data (i.e. DG0), so try
>>
>> File("pressure.pvd")<< p
>>
>> and view the result with ParaView.
>>
>> Garth
>>
>>> A similar "dip" appears with this simple script:
>>>
>>> from dolfin import *
>>> mesh = UnitSquare(20, 20)
>>> V = FunctionSpace(mesh,"DG",0)
>>> 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(mesh,"BDM",1)
>>>> W = FunctionSpace(mesh,"DG",0)
>>>> 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.0-DOLFIN_EPS
>>>>
>>>> u0 = Constant((0.0,0.0))
>>>> bc0 = DirichletBC(V.sub(0), u0, boundary0)
>>>> bc1 = DirichletBC(V.sub(0), u0, boundary1)
>>>> bcs=[bc0,bc1]
>>>>
>>>> f=Expression("-1.0*(x[0]==0.0)+1.0*(x[0]==1.0)")
>>>>
>>>> # 'Permeability' k - gaussian in y
>>>> k=Expression("1.0+exp(-40*pow(x[1]-0.5,2))")
>>>>
>>>> a = (k*dot(v,u) - div(v)*p + q*div(u))*dx
>>>> L = dot(v,f*n)*ds
>>>>
>>>> # Compute solution
>>>> problem = VariationalProblem(a, L,bcs)
>>>> (u, p) = problem.solve().split()
>>>>
>>>> # Plot solution
>>>> plot(p, interactive=True)
>>>> plot(u, interactive=True)
>>>>
>>>>
>>>
>>>
>>> _______________________________________________
>>> Mailing list: https://launchpad.net/~dolfin
>>> Post to : <email address hidden>
>>> Unsubscribe : https://launchpad.net/~dolfin
>>> More help : https://help.launchpad.net/ListHelp
>>
>

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(mesh,"BDM",1)
W = FunctionSpace(mesh,"DG",0)
V = R+W

(u,p) = TrialFunctions(V)
(v,q) = TestFunctions(V)

# boundary conditions
def boundary0(x,on_bound):
    return x[0]>DOLFIN_EPS and x[0]<1.0-DOLFIN_EPS and on_bound

u0 = Constant((0.0,0.0))
bcs = DirichletBC(V.sub(0), u0, boundary0)

f=Expression("-1.0*(x[0]==0.0)+1.0*(x[0]==1.0)")
# f=Constant(0.0)

# 'Permeability' k
k=Expression("1.0+exp(-40*pow(x[1]-0.7,2))")

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

# Compute solution
problem = VariationalProblem(a, L,bcs)
(u, p) = problem.solve().split()

# Plot solution
plot(p, interactive=True)
plot(u, interactive=True)

Chris Richardson wrote:
> Question #105900 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/105900
>
> 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(x,on_bound):
    (*) return x[0]>DOLFIN_EPS and x[0]<1.0-DOLFIN_EPS and on_bound
    bcs = DirichletBC(V.sub(0), u0, boundary0)

only the facets where both endpoints (and the midpoint) satisfy (*) will
be marked. (Somebody correct me if I'm wrong here.) In particular, the
left-most "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

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 mid-point of an edge? or does it have to apply to both ends of the edge?

On Tue, Mar 30, 2010 at 02:01:01PM -0000, Chris Richardson wrote:
> Question #105900 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/105900
>
> 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 mid-point 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

Chris Richardson wrote:
> Question #105900 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/105900
>
> 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

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)
        right = abs(x[0] - 1.0) < DOLFIN_EPS and (abs(x[1] - .5) > DOLFIN_EPS
                                                  and x[1] < 2.0 - DOLFIN_EPS)
        return on_boundary and not (left or right)

On Wed, Mar 31, 2010 at 09:21:50AM -0000, Chris Richardson wrote:
> Question #105900 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/105900
>
> 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 no-slip 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

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?

Anders Logg wrote:
> On Wed, Mar 31, 2010 at 09:21:50AM -0000, Chris Richardson wrote:
>
>> Question #105900 on DOLFIN changed:
>> https://answers.launchpad.net/dolfin/+question/105900
>>
>> 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 no-slip 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
                                                  and x[1] < 1.5 - DOLFIN_EPS)

--
Marie

> --
> Anders
>
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp
>

OK, I think I now understand. Thanks