[Dune] Help needed with pdelab dof blocking

[Oliver Sander]

Hi Steffen,
thanks a lot.  Changing the OrderingTag does indeed solve the problem.
Best,
Oliver

Am 06.01.2014 17:39, schrieb Steffen Müthing:
> Hi Oliver,
> 
> the problem is that you are using the wrong ordering. The backends are correct in your MWE,
> but the OrderingTags aren’t: You want to block your DOFs by the underlying mesh entity, so you
> need to replace the LexicographicOrderingTag by EntityBlockedOrderingTag.
> 
> It’s unfortunate that you don’t get a better error message, but as the DOF mapping framework is
> so flexible, it’s pretty hard to detect invalid combinations of Orderings and Backends…
> 
> Steffen
> 
> Am 06.01.2014 um 15:33 schrieb Oliver Sander <sander at igpm.rwth-aachen.de>:
> 
>> 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
>>>
>>>
>>
>> <testassembly.cc>_______________________________________________
>> 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: signature.asc
Type: application/pgp-signature
Size: 551 bytes
Desc: OpenPGP digital signature
URL: <https://lists.dune-project.org/pipermail/dune/attachments/20140106/28b36b00/attachment.sig>