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 protected]
> > Unsubscribe : https://launchpad.net/~dolfin
> > More help : https://help.launchpad.net/ListHelp
>
signature.asc
Description: Digital signature
_______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : [email protected] Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp
