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