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