[Dune] Help needed with pdelab dof blocking

[Oliver Sander]

Hi all,
I have a follow-up to my question.  While Andreas' suggestion gets the types correct,
I nevertheless cannot assemble the matrix.  Instead, I get an ISTLError thrown, because
some non-existent entry is accessed in the BCRSMatrix.  Supposed I need to tweak another
parameter?  I'll attach a minimal test case that triggers the exception.  Thanks for
your help.
Best,
Oliver

Am 05.01.2014 21:03, schrieb Andreas Nüßing:
> Hi Oliver,
> 
> without me being a PDELab-expert and assuming you are using the PDELab trunk, the following GFS-type does the trick for me. The Blocking template parameter of the VectorBackend for the
> VectorGridFunctionSpace was missing and so defaulted to "no_blocking".
> 
> typedef PDELab::VectorGridFunctionSpace<
>         GridType::LeafGridView,
>         FEM,
>         dim,
> PDELab::ISTLVectorBackend<PDELab::ISTLParameters::static_blocking>,
>         PDELab::ISTLVectorBackend<>,
>         PDELab::ConformingDirichletConstraints,
>         PDELab::LexicographicOrderingTag,
>         PDELab::DefaultLeafOrderingTag
>         > GFS;
> 
> gives the desired output
> 
> Vector type: Dune::BlockVector<Dune::FieldVector<double, 3>, std::allocator<Dune::FieldVector<double, 3> > >
> Matrix type: Dune::BCRSMatrix<Dune::FieldMatrix<double, 3, 3>, std::allocator<Dune::FieldMatrix<double, 3, 3> > >
> 
> Andreas
> 
> 
> 
> Am 05.01.2014 20:02, schrieb Oliver Sander:
>> Dear Dune users,
>> [I'm writing this to the general Dune list instead of the pdelab list,
>> because it is rare (for me) to receive an answer on latter list.]
>>
>> I need help setting up an FEM space.  I want to solve a standard linear
>> elasticity problem on a structured grid using Q1 elements.  ISTL is to
>> be used as the linear algebra backend, and I expect the matrix and vector
>> types pdelab gives me to be
>>
>> BCRSMatrix<FieldMatrix<double,3,3> > and BlockVector<FieldVector<double,3> >
>>
>> Unfortunately, the little test program I attach below gives me
>>
>> BCRSMatrix<FieldMatrix<double,1,1> > and BlockVector<FieldVector<double,1> >
>>
>> instead.  Surely I am using the wrong parameters, but as pdelab is slightly
>> under-documented I am stuck now.  Can somebody please tell me how I need
>> to adjust my code to get the matrix and vector types I want?
>>
>> Many thanks,
>> Oliver
>>
>>
>> #include "config.h"
>>
>> #include<iostream>
>> #include<dune/common/fvector.hh>
>>
>> #include<dune/grid/yaspgrid.hh>
>>
>> #include<dune/pdelab/finiteelementmap/q1fem.hh>
>> #include<dune/pdelab/constraints/conforming.hh>
>> #include<dune/pdelab/gridfunctionspace/vectorgridfunctionspace.hh>
>> #include<dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh>
>> #include<dune/pdelab/gridfunctionspace/interpolate.hh>
>> #include <dune/pdelab/gridfunctionspace/vtk.hh>
>> #include<dune/pdelab/constraints/common/constraints.hh>
>> #include<dune/pdelab/common/function.hh>
>> #include<dune/pdelab/backend/istlvectorbackend.hh>
>> #include<dune/pdelab/backend/istlmatrixbackend.hh>
>> #include<dune/pdelab/localoperator/linearelasticity.hh>
>> #include<dune/pdelab/constraints/constraintsparameters.hh>
>> #include<dune/pdelab/gridoperator/gridoperator.hh>
>>
>> using namespace Dune;
>>
>> const int dim = 3;
>>
>> typedef YaspGrid<dim> GridType;
>>
>> typedef PDELab::Q1LocalFiniteElementMap<double,double,dim> FEM;
>> typedef PDELab::VectorGridFunctionSpace<
>>    GridType::LeafGridView,
>>    FEM,
>>    dim,
>>    PDELab::ISTLVectorBackend<>,
>>    PDELab::ISTLVectorBackend<>,
>>    PDELab::ConformingDirichletConstraints,
>>    PDELab::LexicographicOrderingTag,
>>    PDELab::DefaultLeafOrderingTag
>>    > GFS;
>>
>> typedef GFS::ConstraintsContainer<double>::Type C;
>>
>> typedef PDELab::GridOperator<GFS,GFS,PDELab::LinearElasticity,PDELab::ISTLMatrixBackend,double,double,double,C,C> GO;
>> typedef GO::Traits::Domain V;
>> typedef GO::Jacobian M;
>> typedef M::BaseT ISTL_M;
>> typedef V::BaseT ISTL_V;
>>
>>
>> int main(int argc, char** argv)
>> {
>>      std::cout << "Vector type: " << className<ISTL_V>() << std::endl;
>>      std::cout << "Matrix type: " << className<ISTL_M>() << std::endl;
>>
>> }
>>
>>
>>
>>
>>
>> _______________________________________________
>> Dune mailing list
>> Dune at dune-project.org
>> http://lists.dune-project.org/mailman/listinfo/dune
> 
> 

-------------- next part --------------
A non-text attachment was scrubbed...
Name: testassembly.cc
Type: text/x-c++src
Size: 2684 bytes
Desc: not available
URL: <https://lists.dune-project.org/pipermail/dune/attachments/20140106/0923e5a6/attachment.cc>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 551 bytes
Desc: OpenPGP digital signature
URL: <https://lists.dune-project.org/pipermail/dune/attachments/20140106/0923e5a6/attachment.sig>