[Dune] Help needed with pdelab dof blocking

[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: 495 bytes
Desc: Message signed with OpenPGP using GPGMail
URL: <https://lists.dune-project.org/pipermail/dune/attachments/20140106/5c71c820/attachment.sig>