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 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 : 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.