dolfin team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

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

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. 

> > 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?

Johan