dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Sun, Dec 07, 2008 at 11:09:11PM +0100, Johan Hake wrote:
> On Sunday 07 December 2008 22:34:42 Anders Logg wrote:
> > On Sun, Dec 07, 2008 at 10:04:51PM +0100, Johan Hake wrote:
> > > On Sunday 07 December 2008 20:21:05 Anders Logg wrote:
> > > > On Sun, Dec 07, 2008 at 07:59:40PM +0100, Johan Hake wrote:
> > > > > On Sunday 07 December 2008 16:42:02 Anders Logg wrote:
> > > > > > On Sun, Dec 07, 2008 at 02:57:01PM +0100, Martin Sandve Alnæs wrote:
> > > > > > > 2008/12/6 Johan Hake <hake@xxxxxxxxx>:
> > > > > > > > On Saturday 06 December 2008 20:54:34 Anders Logg wrote:
> > > > > > > >> On Sat, Dec 06, 2008 at 04:03:25PM +0100, Johan Hake wrote:
> > > > > > > >> > On Saturday 06 December 2008 15:43:13 Anders Logg wrote:
> > > > > > > >> > > On Sat, Dec 06, 2008 at 03:34:40PM +0100, Johan Hake wrote:
> > > > > > > >> > > > On Saturday 06 December 2008 14:56:59 Anders Logg wrote:
> > > > > > > >> > > > > On Sat, Dec 06, 2008 at 12:05:41PM +0100, DOLFIN wrote:
> > > > > > > >> > > > > > One or more new changesets pushed to the primary
> > > > > > > >> > > > > > dolfin repository. A short summary of the last three
> > > > > > > >> > > > > > changesets is included below.
> > > > > > > >> > > > > >
> > > > > > > >> > > > > > changeset:
> > > > > > > >> > > > > > 5262:0e349fbe09ce4179252652fe5c1d58725951bc3f tag:
> > > > > > > >> > > > > >   tip
> > > > > > > >> > > > > > user:        "Johan Hake <hake@xxxxxxxxx>"
> > > > > > > >> > > > > > date:        Sat Dec 06 12:05:42 2008 +0100
> > > > > > > >> > > > > > files:      
> > > > > > > >> > > > > > demo/pde/stokes/taylor-hood/python/demo.py
> > > > > > > >> > > > > > dolfin/function/SpecialFunctions.h
> > > > > > > >> > > > > > dolfin/swig/dolfin_function_pre.i
> > > > > > > >> > > > > > site-packages/dolfin/function.py description:
> > > > > > > >> > > > > > Added support for f.split()
> > > > > > > >> > > > > >  - Renamed operator[] to f._sub instead of f.sub
> > > > > > > >> > > > > >  - f.sub(i) now returns an instantiated sub function
> > > > > > > >> > > > > >  - f.split() uses f.sub() to return a tuple of all
> > > > > > > >> > > > > > sub functions - stoke/taylor-hood demo now runs.
> > > > > > > >> > > > >
> > > > > > > >> > > > > Excellent!
> > > > > > > >> > > > >
> > > > > > > >> > > > > It runs now but the solution looks completely crazy.
> > > > > > > >> > > > > The problem is that the boundary conditions are not
> > > > > > > >> > > > > set correctly since we use V and Q to set the boundary
> > > > > > > >> > > > > conditions for the sub systems and they don't know the
> > > > > > > >> > > > > offsets
> > > > > > > >> > > > > (DofMap::offset()).
> > > > > > > >> > > >
> > > > > > > >> > > > The problem here is that these spaces need to be
> > > > > > > >> > > > SubSpaces?
> > > > > > > >> > >
> > > > > > > >> > > Yes.
> > > > > > > >> > >
> > > > > > > >> > > > If that is the case, we could extract the subspaces
> > > > > > > >> > > > after a MixedFunctionSpace is created and then store
> > > > > > > >> > > > these in the spaces attribute, either as pure
> > > > > > > >> > > > cpp.SubSpaces or add another python class, SubSpace,
> > > > > > > >> > > > which is a cpp.SubSpace and stores the original
> > > > > > > >> > > > ffc.element too, or something?
> > > > > > > >> > >
> > > > > > > >> > > I was thinking something like this:
> > > > > > > >> > >
> > > > > > > >> > > W = VectorFunctionSpace(mesh, "triangle", 2) +
> > > > > > > >> > > FunctionSpace(mesh, "triangle", 1) V, Q = W.split()
> > > > > > > >> >
> > > > > > > >> > Doesn't look too nice, and it will create a non intuitive
> > > > > > > >> > work flow. Combining two interfaces aren't easy!
> > > > > > > >> >
> > > > > > > >> > > The split function needs to both create SubSpaces (which
> > > > > > > >> > > will lead to DofMaps with correct offsets) and set the
> > > > > > > >> > > element correctly, which can be done by looking at the
> > > > > > > >> > > spaces attribute in MixedFunctionSpace.
> > > > > > > >> > >
> > > > > > > >> > > It's a bit weird since what we really do is
> > > > > > > >> > >
> > > > > > > >> > > V = VectorFunctionSpace(mesh, "triangle", 2)
> > > > > > > >> > > Q = FunctionSpace(mesh, "triangle", 1)
> > > > > > > >> > > W = V + Q
> > > > > > > >> > > V, Q = W.split()
> > > > > > > >> > >
> > > > > > > >> > > First, we put V and Q together to create a mixed function
> > > > > > > >> > > space. Then we split W again into V and Q. After the split
> > > > > > > >> > > V and Q know that they are part of the bigger space (at
> > > > > > > >> > > least they know the offset into the bigger space).
> > > > > > > >> >
> > > > > > > >> > Yes I see this. Instead, or in addition, of using split, we
> > > > > > > >> > could implement W.sub(i) that returns a subspace, which is
> > > > > > > >> > then used when setting the bc, and other relevant places.
> > > > > > > >> >
> > > > > > > >> > > The other option would be to let + have a side effect on V
> > > > > > > >> > > and Q but that does not seem to be a good solution.
> > > > > > > >> >
> > > > > > > >> > I thought of this too. But I cannot se how we could do it
> > > > > > > >> > with the present implementation of FunctionSpace/Subspace.
> > > > > > > >>
> > > > > > > >> If we can get make_subspace working as you suggest below, then
> > > > > > > >> it would be easy to modify __add__ for FunctionSpaceBase:
> > > > > > > >>
> > > > > > > >>   W = MixedFunctionSpace([self, other])
> > > > > > > >>   self.make_subspace(W, 0)
> > > > > > > >>   other.make_subspace(W, 1)
> > > > > > > >>   return W
> > > > > > > >
> > > > > > > > Yes something like that.
> > > > > > > >
> > > > > > > >> This would make the interface look nice but someone (you know
> > > > > > > >> who you are) might think we are insane... :-)
> > > > > > > >
> > > > > > > > *laugh*
> > > > > > > >
> > > > > > > > I am not very happy with it either, as the __add__ operator
> > > > > > > > should not do such a thing. What we are saying is that V and Q
> > > > > > > > are individual FunctionSpaces before we make make the add.
> > > > > > > > After we have made it they are suddenly subspaces of the result
> > > > > > > > of an addition. Not very intuitive.
> > > > > > > >
> > > > > > > > One slution could be to remove the __add__ operator and instead
> > > > > > > > force the user to instantiate the Mixed space with a
> > > > > > > > constructor, e.g.
> > > > > > > >
> > > > > > > >   W = MixedFunctionSpace(V,Q)
> > > > > > > >
> > > > > > > > With this it is easier to justify that we are doing somthing
> > > > > > > > with V and Q.
> > > > > > >
> > > > > > > No it isn't.
> > > > > > >
> > > > > > > This isn't a matter of taste, it just doesn't scale outside a
> > > > > > > minimal example.
> > > > > > >
> > > > > > > What will happen here?
> > > > > > >   W0 = MixedFunctionSpace(V0,V1) # or V0 + V1, same thing
> > > > > > >   W1 = MixedFunctionSpace(V1,V2)
> > > > > > >   W2 = MixedFunctionSpace(V2,V3)
> > > > > > >
> > > > > > > >> > If we could just call a member function in a FunctionSpace,
> > > > > > > >> > let say V.make_subspace(W,i), and this would then create the
> > > > > > > >> > needed stuff in V, to make it SubSpaceable. But now SubSpace
> > > > > > > >> > is a subclass which effectively prevents this approach.
> > > > > > > >> >
> > > > > > > >> > This approach is abit intrucive too. But what difference
> > > > > > > >> > would it make for the user, i.e., if V in addition to be a
> > > > > > > >> > FunctionSpace now also is a> >> > subspace of W?
> > > > > > > >>
> > > > > > > >> The only reason we need this is to be able to set boundary
> > > > > > > >> conditions for sub systems. For that we need two things: the
> > > > > > > >> offset into the global system vector and the DofMap for the
> > > > > > > >> subsystem in question. There might be other ways to define the
> > > > > > > >> interface for DirichletBC, for example:
> > > > > > > >>
> > > > > > > >>   bc = DirichletBC(W, 0, ...)
> > > > > > > >
> > > > > > > > This is maybee the best one? The DirichletBC operates at the
> > > > > > > > global vector from the mixed space, and then it is more
> > > > > > > > intuitive to actually use the mixed space to set the BC.
> > > > > > >
> > > > > > > Agree.
> > > > > >
> > > > > > I'm still not convinced. The thing is that the C++ implementation
> > > > > > of DirichletBC assumes that it gets a FunctionSpace. Sometimes,
> > > > > > this may be a SubSpace (which is a subclass of FunctionSpace) but
> > > > > > DirichletBC doesn't care, it just sets the bcs for the
> > > > > > FunctionSpace. Before, in 0.8.1, there were 6 different
> > > > > > constructors and an extra SubSystem argument was required, but the
> > > > > > introduction of SubSpace simplified the implementation.
> > > > > >
> > > > > > This all works very fine in C++, but only because we don't see the
> > > > > > initial FunctionSpaces V and Q that are used to create the
> > > > > > MixedFunctionSpace W. We just get the MixedFunctionSpace W directly
> > > > > > and then it's natural to define V and Q from W. The problem in
> > > > > > Python is that we see the original V and Q, so it becomes awkward
> > > > > > to put them together and then split them:
> > > > > >
> > > > > >   W = V + Q
> > > > > >   (V, Q) = W.split()
> > > > > >
> > > > > > How about the following (just renaming the variables):
> > > > > >
> > > > > >   W = V + Q
> > > > > >   (W0, W1) = W.split()
> > > > > >
> > > > > >   bc0 = DirichletBC(W0, domain, u0)
> > > > > >   bc1 = DirichletBC(W1, domain, p0)
> > > > >
> > > > > If we implement this as we did with the Function.split(), we could
> > > > > also support:
> > > > >
> > > > >    W = V + Q
> > > > >
> > > > >    bc0 = DirichletBC(u0, W.sub(0), domain)
> > > > >    bc1 = DirichletBC(p0, W.sub(1), domain)
> > > >
> > > > Looks good.
> > > >
> > > > What about the order of the arguments? Should the value be first, or
> > > > the FunctionSpace?
> > >
> > > No opinion. Just changed them to illustrate the present implementation.
> > >
> > > > > It is a bit fragile though, as the user can potentially send in the
> > > > > wrong FunctionSpace. The only difference between V and W.sub(0), as I
> > > > > understand it, is that the latter includes information about the
> > > > > offset into the global vector.
> > > >
> > > > Yes.
> > > >
> > > > > Would it be possible to add a check in DirichletBC.apply (or/and in
> > > > > BoundaryCondition.apply) that check the number of dofs in the vector
> > > > > it is acting on and compare that with the number of dofs in the
> > > > > FunctionSpace?
> > > > >
> > > > > If the number of dofs in the vector is larger, then the FuncitonSpace
> > > > > have to be a SubSpace, and issue an error or arning if it is not?
> > > >
> > > > Yes, that would be good to add. Using dynamic_cast?
> > >
> > > Or, we could just check for V.dof_map().global_dimension() == x.size().
> >
> > But we also need to check that if the above does not hold, then V must
> > be a subspace. Something like this (in pseudo-code):
> >
> >   if V.dim() < x.size() and not V is SubSpace:
> >       error("Illegal dimension for function space.")
> 
> Yes, of course, I forgot the last part :P 
> 
> What with:
> 
>   V.dof_map().global_dimension() == x.size() and V.dof_map().offset() != 0

We need to check more: the dimension cannot be larger than x.size().
Also, the offset may be zero also for a SubSpace (the first one).

How about this:

    if (V.dof_map().global_dimension() > x.size())
      error("...");
    else if (V.dof_map().global_dimension() < x.size() && not isinstance(V, SubSpace))
      error("...");

> I have added sub to MixedFunctionSpace. But it does not seem to work. There 
> also seems to be an issue with swig garbage collecting the created subspace, 
> so wee need to keep a reference. 

Strange. I've added some printing of offsets and it looks right.

> > > I also see that there are no dimension check for the vectors and matrices
> > > in the apply function.
> > >
> > >   A.size(0) == A.size(1) == b.size() == x.size()
> > >
> > > We should assume the user throws anything at it. This probably goes for a
> > > lot of functions.
> >
> > Yes! We need lots of more checks in many places...
> 
> Is it a registered bug?

I don't know. Everything is in such a flux right now so I rarely look
at the bug tracking system.

-- 
Anders

Attachment: signature.asc
Description: Digital signature