dolfin team mailing list archive

["Garth N. Wells" <gnw20@xxxxxxxxx>]

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