[dune-functions] Crash in interpolate

[Carsten Gräser]

Am 20.09.2015 um 20:08 schrieb Oliver Sander:
> Hi Carsten,
>
> [snip]
>>
>> Instead of this it will always interpolate into a vector
>> for the whole global basis using the original global indices.
>> So you would need a
>>
>>   BlockVector<BlockVector<FieldVector<double,1>>>
>>
>> to interpolate the indentity as a vector field, if this is
>> really what you want.
> Thanks for the explanation.  What you propose works as expected.
> Can the code be changed to give an error message instead of the
> segfault?
I was thunking of this, but it's very complicated. The reason
for this is that the nesting structure of your vector is in principle
OK, because it allows to nested indices. The real problem is, that
the inner vector can't be resized. Which currently silently
ommits the resizing. We could (in HierarchicalWrapper.resize())
add a check for this. But then the check would also be run
for each individual entry of large, dynamically sized containers.
Avoiding this would be very complicated because the wrapper would
have know a lot more information about the index structure.
But having checks that are disabled on NDEBUG would be viable.

>>  But then you have the x- and y-coordinates
>> in the first and second velocity components. If you want the
>> Lagrange-points for all velocity dofs, you need the implicite
>> tensor construction, i.e.,
>>
>>   BlockVector<BlockVector<FieldVector<double,dim>>>
> This does not work for me.  It compiles and runs, but the interpolated
> values
> are not correct (each FieldVector has the same value in all entries). 
> However,
> this is not a problem to me.  BV<BV<FV<,dim>>> does not fit the
> multi-indices
> of the TaylorHoodBasis, so it is not surprising that this doesn't work.
It does not fit directly. But it does still work due to the implicit
product space feature I've implemented (in the same way we does it in
dune-fufem). Maybe there's another missunderstanding. I suppose you're
now trying the following:

  BlockVector<BlockVector<FieldVector<double,dim> > > lagrangeNodes;
  
  interpolate(taylorHoodBasis, Dune::TypeTree::hybridTreePath(_0), lagrangeNodes, [](Coordinate x) {return x;});

This will result in the behavior you described. I'll explain how the
subtree interpolation and implicit product spaces work:

Because a subtree corresponds to a product space 'interpolate'
will first map the subtree leaf nodes (which correspond to the
factors of the product space) to the components of the function
range. In the above code the number subtree leaf nodes is dim.
Hence it will associate x[i] to the i-th velocity component.
No it assigns to the dof of the velocity component the associated
function value x[i]. Since the dof is a FV<dim> while x[i] is a scalar,
you get the desired behavior.

For the implicit product feature to drop in you have to blow up
the coefficient vector and the function in the same way. For
the hybridTreePath(_0) the function is already expected to
be a vector field if the dofs are scalar. If the dofs are vector
valued, each component of the vector field should be vector valued,
too. I.e., your function should be matrix valued.

What you really want is to use the same vector valued function
to interpolate for the vector valued coefficients of each velocity
component. I.e.

  BlockVector<BlockVector<FieldVector<double,dim> > > lagrangeNodes;
  for(int i=0; i<dim; ++i)
    interpolate(taylorHoodBasis, Dune::TypeTree::hybridTreePath(_0, i), lagrangeNodes, [](Coordinate x) {return x;});


It takes a moment to understand the mechanism, but in fact
it's quite natural.

BTW: It's even more fun if you consider vector valued basis functions.
Have a look at the documentation of the new DiscreteGlobalBasisFunction.

Best,
Carsten


-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 473 bytes
Desc: OpenPGP digital signature
URL: <https://lists.dune-project.org/pipermail/dune-functions/attachments/20150920/ebb8322e/attachment.sig>