dolfin team mailing list archive

[Johan Hake <johan.hake@xxxxxxxxx>]

On Monday November 15 2010 09:57:36 David Maxwell wrote:
> Question #133939 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/133939
> 
>     Status: Answered => Open
> 
> David Maxwell is still having a problem:
> This would be too easy. I'm looking for  a function defined on (a part
> of) the boundary, not just its average value.

Ok.

Restricting your FunctionSpace to a geometrical subdomain is not supported 
(yet) in DOLFIN.

You should, however, be able to hack the solution (not pretty, and maybe not 
totaly correct wrt FEM analysis) by solving a reduces system if you are using 
'CG' 1 elements.

 lhs = assemble(uu*vv*ds(bottom), exterior_facet_domains=mesh_function)
 rhs = assemble(inner(grad(u),grad(vv))*dx - vv*f*dx)

 verts = VertexFunction('uint', mesh)

 # Mark the vertices at the boundary with the bottom marker
 inds = (verts.array()==bottom)

 # Get the reduces system
 rhs0 = rhs[inds] # Approximation of the rhs contribution 
 lhs0 = lhs[inds,inds] # Just getting rid of zeros
 x = lhs0.copy()

 solve(lhs0,x,rhs0)

 # Put reduced solution on the full solution function
 u.vector()[:] = 0.0
 u.vector()[inds] = x

 plot(u)
 interactive()

Johan

> -David
> 
> On Nov 15, 2010, at 8:33 AM, Johan Hake wrote:
> > Your question #133939 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/133939
> > 
> >    Status: Open => Answered
> > 
> > Johan Hake proposed the following answer:
> > Couldn't you just assemble the boundary condition directly from the
> > solution?
> > 
> >  neuman_vector = assemble(u*vv*ds(bottom),\
> >  
> >                  exterior_facet_domains=mesh_function)
> > 
> > If this is what you want you can preassemble a boundary mass matrix:
> >  Mb = assemble(uu*vv*ds(bottom), exterior_facet_domains=mesh_function)
> > 
> > and then each times step:
> >  neuman_vector = Mb*u.vector()
> > 
> > Johan
> > 
> > On Saturday November 13 2010 12:58:06 David Maxwell wrote:
> >> New question #133939 on DOLFIN:
> >> https://answers.launchpad.net/dolfin/+question/133939
> >> 
> >> I'd like to extract Neumann data from a solution to an elliptic PDE.
> >> 
> >> For example, if u is a weak solution of -Laplacian(u) = f, then its
> >> Neumann data (\partial_n u) on the boundary is well defined as follows:
> >> 
> >> Given sufficiently regular v defined on the boundary, extend v to a
> >> function in H^1 in the domain.  Then
> >> 
> >> <\partial_n u, v> = a_\Omega \nabla u \nabla v - f v
> >> 
> >> Since v belongs to L^2 of the boundary, I can interpret \partial_n as an
> >> element of L^2 of the boundary.
> >> 
> >> In truth, I really only want the Neumann data on part of the boundary.
> >> 
> >> I thought I might proceed as follows.  In the following, 'u' is the
> >> solution of the PDE in function space V, and mesh_function is a mesh
> >> function on the edges that equals 1 on the 'top' boundary, 2 on the
> >> remainder of the boundary, and 0 on all other edges.
> >> 
> >> My first attempt went something like:
> >> 
> >> uu = TestFunction(V)
> >> vv = TrialFunction(V)
> >> 
> >> inside=0; top=1; bottom=2
> >> 
> >> a = uu*vv*ds(bottom)
> >> rhs = inner(grad(u),grad(vv))*dx - vv*f*dx
> >> 
> >> top_bc = DirichletBC(V,Constant(0),mesh_function,top)
> >> interior_bc = DirichletBC(V,Constant(0),mesh_function,inside) # I
> >> thought I was being clever here!
> >> 
> >> u_n =
> >> VariationalProblem(a,rhs,bcs=[top_bc,interior_bc],exterior_facet_domains
> >> =m esh_function).solve()
> >> 
> >> This crashes on applying the boundary conditions using 'ident', so I
> >> tried setting the 'use_ident' parameter to False on the boundary
> >> conditions.
> >> 
> >> This still doesn't work because the interior_bc is (rightfully) too
> >> enthusiastic -- every vertex is joined to some interior edge and so I
> >> pick up the zero solution.
> >> 
> >> I tried some other tacts, but they are perhaps too foolish to describe
> >> here.
> >> 
> >> And I have the sense that I really ought to be solving a problem on a
> >> BoundaryMesh (rather than a highly constrained problem on the original
> >> mesh) but I don't know how to extend Functions on a boundary mesh to
> >> Functions on the original mesh.
> >> 
> >> Help!
> >> 
> >> David Maxwell
> >> 
> >> 
> >> 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@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.
> 
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~dolfin
> More help   : https://help.launchpad.net/ListHelp