dolfin team mailing list archive

[Garth Wells <question105900@xxxxxxxxxxxxxxxxxxxxx>]

Question #105900 on DOLFIN changed:
https://answers.launchpad.net/dolfin/+question/105900

Garth Wells proposed the following answer:

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 : dolfin@xxxxxxxxxxxxxxxxxxx
>>> Unsubscribe : https://launchpad.net/~dolfin
>>> More help : https://help.launchpad.net/ListHelp
>>
>

-- 
You received this question notification because you are a member of
DOLFIN Team, which is an answer contact for DOLFIN.