dolfin team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

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)

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.

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?

Johan