[dune-pdelab] [dune-pdelab-commit] dune-pdelab r471 - trunk/dune/pdelab/gridfunctionspace
Bernd Flemisch
bernd at iws.uni-stuttgart.de
Mon May 17 16:10:38 CEST 2010
Hey Christian,
I am still not able to come up with a stripped down example, but with
some more details:
I use
typedef Dune::PDELab::Q1LocalFiniteElementMap<Scalar,Scalar,dim> FEM;
typedef Dune::PDELab::GridFunctionSpace<GridView, FEM, Constraints,
Dune::PDELab::ISTLVectorBackend<2> > ScalarGridFunctionSpace;
typedef Dune::PDELab::PowerGridFunctionSpace<ScalarGridFunctionSpace, 2,
Dune::PDELab::GridFunctionSpaceBlockwiseMapper> GridFunctionSpace;
The problem is that it does not enter the specialization
template<typename T, int k, int s>
class
PowerGridFunctionSpaceBase<T,k,GridFunctionSpaceComponentBlockwiseMapper<s>
>
at line 1340 in gridfunctionspace.hh, but the general
template<typename T, int k, typename P>
class PowerGridFunctionSpaceBase
at line 1159. But I guess that it is intended to enter the
specialization if GridFunctionSpaceBlockwiseMapper is used, right?
If I use
typedef Dune::PDELab::PowerGridFunctionSpace<ScalarGridFunctionSpace, 2,
Dune::PDELab::GridFunctionSpaceComponentBlockwiseMapper<1> >
GridFunctionSpace;
instead, it enters the specialization and everything works as before.
If it does not enter the specialization, the indexing gets mixed up
between elementwise and phys.wise.
I tried to build a specializiation for
template<typename T, int k>
class PowerGridFunctionSpaceBase<T,k,GridFunctionSpaceBlockwiseMapper>
My first efforts failed miserably. In principle, would this be the
correct way to go?
Kind regards
Bernd
On 05/17/2010 02:13 PM, Christian Engwer wrote:
> Up to recently the you had two options regarding the blockin of you
> dofs. Assume two cells A,B with two physical values x,y and two DOFs
> 1,2.
>
> you can have Elementwise blocking of physical value wise
> blocking. This results two different numberings:
>
> Elementwise : [Ax1, Ay1, Ax2, Ay2, Bx1, By1, By2, By2]
> Phys. wise : [Ax1, Ax2, Bx1, Bx2, Ay1, Ay2, By1, By2]
>
> For elementwise blocking it was necessary that all physical values
> have the same number of dofs. This rendered certain setups
> impossible. Think of Stokes equation.
>
> For a Stoeks DG discretitzation you have a function space like this:
> (P_k)^dim x P_{k-1}
>
> In order to benefit from you DG structure you want to do element wise
> blocking and you want a structure in the following way:
> [Av1_1, .. Av1_n, Av2_1, .. Av2_n, Av3_1, .. Av3_n, Ap_1, .. Ap_m, Bv1_1, ...]
>
> The changes in the code are intended to allow exactly this
> behaviour. I tried to keep them backwards compatible.
>
> You please describe which exact setup you have and how the indices
> changed. Please give a boiled down test case ;-)
>
> Christian
>
> On Mon, May 17, 2010 at 11:17:59AM +0200, Bernd Flemisch wrote:
>
>> Dear Dune-PDELab, hey Christian,
>>
>> this commit unfortunately breaks my code in the sense that I obtain
>> nonsense solutions. I guess that some assumption I put in on the
>> indexing is not valid any more. Could you maybe explain your commit?
>>
>> Thank you. Kind regards
>> Bernd
>>
>> On 05/07/2010 08:05 PM, christi at conan.iwr.uni-heidelberg.de wrote:
>>
>>> Author: christi
>>> Date: Fri May 7 20:05:07 2010
>>> New Revision: 471
>>> URL: http://svn.dune-project.org/websvn/listing.php?repname=dune-pdelab&path=/&rev=471&sc=1
>>>
>>> Log:
>>> allow to change the blocking structure
>>>
>>> Modified:
>>> trunk/dune/pdelab/gridfunctionspace/gridfunctionspace.hh
>>>
>>> Modified: trunk/dune/pdelab/gridfunctionspace/gridfunctionspace.hh
>>> ==============================================================================
>>> --- trunk/dune/pdelab/gridfunctionspace/gridfunctionspace.hh Fri May 7 11:34:26 2010 (r470)
>>> +++ trunk/dune/pdelab/gridfunctionspace/gridfunctionspace.hh Fri May 7 20:05:07 2010 (r471)
>>> @@ -4,22 +4,22 @@
>>> #define DUNE_PDELAB_GRIDFUNCTIONSPACE_HH
>>>
>>> #include <cstddef>
>>> -#include<map>
>>> +#include <map>
>>> #include <ostream>
>>> -#include<set>
>>> -#include<vector>
>>> +#include <set>
>>> +#include <vector>
>>>
>>> -#include<dune/common/exceptions.hh>
>>> -#include<dune/common/geometrytype.hh>
>>> -#include<dune/common/stdstreams.hh>
>>> -#include<dune/grid/common/genericreferenceelements.hh>
>>> +#include <dune/common/exceptions.hh>
>>> +#include <dune/common/geometrytype.hh>
>>> +#include <dune/common/stdstreams.hh>
>>> +#include <dune/grid/common/genericreferenceelements.hh>
>>>
>>> #include <dune/localfunctions/common/localkey.hh>
>>>
>>> -#include"../common/countingptr.hh"
>>> -#include"../common/multitypetree.hh"
>>> -#include"../common/cpstoragepolicy.hh"
>>> -#include"../common/geometrywrapper.hh"
>>> +#include "../common/countingptr.hh"
>>> +#include "../common/multitypetree.hh"
>>> +#include "../common/cpstoragepolicy.hh"
>>> +#include "../common/geometrywrapper.hh"
>>>
>>> #include"localfunctionspace.hh"
>>>
>>> @@ -175,12 +175,8 @@
>>> /** \brief Tag indicating an arbitrary number of unkowns per entity.
>>> *
>>> * class used to pass compile-time parameter to the GridFunctionSpace.
>>> - *
>>> */
>>> - struct GridFunctionGeneralMapper
>>> - {
>>> - enum {dummy=0} ;
>>> - };
>>> + struct GridFunctionGeneralMapper {};
>>>
>>>
>>> // Empty constraints assembler class
>>> @@ -494,17 +490,12 @@
>>> std::set<unsigned int> codimUsed;
>>> };
>>>
>>> -
>>> -
>>> /** \brief Tag indicating a fixed number of unkowns per entity (known at compile time).
>>> *
>>> * class used to pass compile-time parameter to the GridFunctionSpace.
>>> *
>>> */
>>> - struct GridFunctionRestrictedMapper
>>> - {
>>> - enum {dummy=1} ;
>>> - };
>>> + struct GridFunctionRestrictedMapper {};
>>>
>>> //! \}
>>>
>>> @@ -782,7 +773,6 @@
>>> template<typename IIS>
>>> struct GridFunctionStaticSize
>>> {
>>> - enum {dummy=2} ;
>>> typedef IIS IntersectionIndexSet;
>>> };
>>>
>>> @@ -1114,35 +1104,57 @@
>>> //! \brief Indicates lexicographics ordering of the unknowns for composite
>>> //! grid function spaces.
>>> //!
>>> - //! this class may be used to pass compile-time
>>> - //! parameters to the implementation of
>>> + //! this class may be used to pass compile-time
>>> + //! parameters to the implementation of
>>> //! \link PowerGridFunctionSpace PowerGridFunctionSpace \endlink or
>>> //! \link CompositeGridFunctionSpace CompositeGridFunctionSpace \endlink
>>> - struct GridFunctionSpaceLexicographicMapper
>>> - {
>>> - enum {dummy=0} ;
>>> - };
>>> + struct GridFunctionSpaceLexicographicMapper {};
>>>
>>> //! \brief Indicates using block-wise ordering of the unknowns for composite
>>> //! grid function spaces.
>>> //!
>>> - //! this class may be used to pass compile-time
>>> - //! parameters to the implementation of
>>> + //! The exact blocking structure can be passed as template parameters
>>> + //!
>>> + //! this class may be used to pass compile-time
>>> + //! parameters to the implementation of
>>> //! \link PowerGridFunctionSpace PowerGridFunctionSpace \endlink or
>>> //! \link CompositeGridFunctionSpace CompositeGridFunctionSpace \endlink
>>> - struct GridFunctionSpaceBlockwiseMapper
>>> - {
>>> - enum {dummy=0} ;
>>> - };
>>> + template<int s0 = 1, int s1 = 1, int s2 = 1, int s3 = 1, int s4 = 1, int s5 = 1, int s6 = 1, int s7 = 1, int s8 = 1, int s9 = 1>
>>> + struct GridFunctionSpaceComponentBlockwiseMapper
>>> + {
>>> + static const int size[];
>>> + static const int offset[];
>>> + };
>>> + template<int s0, int s1, int s2, int s3, int s4, int s5, int s6, int s7, int s8, int s9>
>>> + const int GridFunctionSpaceComponentBlockwiseMapper<s0,s1,s2,s3,s4,s5,s6,s7,s8,s9>::
>>> + size[] = { s0, s1, s2, s3, s4, s5, s6, s7, s8, s9 };
>>> + template<int s0, int s1, int s2, int s3, int s4, int s5, int s6, int s7, int s8, int s9>
>>> + const int GridFunctionSpaceComponentBlockwiseMapper<s0,s1,s2,s3,s4,s5,s6,s7,s8,s9>::
>>> + offset[] = { 0, s0, s0+s1, s0+s1+s2, s0+s1+s2+s3, s0+s1+s2+s3+s4,
>>> + s0+s1+s2+s3+s4+s5, s0+s1+s2+s3+s4+s5+s6, s0+s1+s2+s3+s4+s5+s6+s7,
>>> + s0+s1+s2+s3+s4+s5+s6+s7+s8, s0+s1+s2+s3+s4+s5+s6+s7+s8+s9 };
>>> +
>>> + //! \brief Indicates using block-wise ordering of the unknowns for composite
>>> + //! grid function spaces.
>>> + //!
>>> + //! this class may be used to pass compile-time
>>> + //! parameters to the implementation of
>>> + //! \link PowerGridFunctionSpace PowerGridFunctionSpace \endlink or
>>> + //! \link CompositeGridFunctionSpace CompositeGridFunctionSpace \endlink
>>> + struct GridFunctionSpaceBlockwiseMapper : GridFunctionSpaceComponentBlockwiseMapper<> {};
>>>
>>> //! \}
>>>
>>> template<typename T, int k, typename P>
>>> class PowerGridFunctionSpace;
>>>
>>> - // product of identical grid function spaces
>>> - // base class that holds implementation of the methods
>>> - // this is the default version with lexicographic ordering
>>> + //! product of identical grid function spaces
>>> + //! base class that holds implementation of the methods
>>> + //! this is the default version with lexicographic ordering
>>> + //!
>>> + //! \tparam T PLEASE DOCUMENT
>>> + //! \tparam k PLEASE DOCUMENT
>>> + //! \tparam P PLEASE DOCUMENT
>>> template<typename T, int k, typename P>
>>> class PowerGridFunctionSpaceBase
>>> : public PowerNode<T,k,CountingPointerStoragePolicy>,
>>> @@ -1297,7 +1309,7 @@
>>> private:
>>> void setup ()
>>> {
>>> - Dune::dinfo << "power grid function space(lexicographic version):"
>>> + Dune::dinfo << "PowerGridFunctionSpace(lexicographic version):"
>>> << std::endl;
>>> Dune::dinfo << "( ";
>>> offset[0] = 0;
>>> @@ -1325,18 +1337,17 @@
>>> // product of identical grid function spaces
>>> // base class that holds implementation of the methods
>>> // specialization for blockwise ordering
>>> - template<typename T, int k>
>>> - class PowerGridFunctionSpaceBase<T,k,GridFunctionSpaceBlockwiseMapper>
>>> + template<typename T, int k, int s>
>>> + class PowerGridFunctionSpaceBase<T,k,GridFunctionSpaceComponentBlockwiseMapper<s> >
>>> : public PowerNode<T,k,CountingPointerStoragePolicy>,
>>> public Countable
>>> {
>>> - friend class PowerGridFunctionSpace<T,k,GridFunctionSpaceBlockwiseMapper>;
>>> -
>>> + friend class PowerGridFunctionSpace<T,k,GridFunctionSpaceComponentBlockwiseMapper<s> >;
>>> public:
>>> //! export traits class
>>> typedef PowerCompositeGridFunctionSpaceTraits<typename T::Traits::GridViewType,
>>> typename T::Traits::BackendType,
>>> - GridFunctionSpaceBlockwiseMapper, k>
>>> + GridFunctionSpaceComponentBlockwiseMapper<s>, k>
>>> Traits;
>>>
>>> //! extract type of container storing Es
>>> @@ -1408,7 +1419,7 @@
>>> template<int i>
>>> typename Traits::SizeType subMap (typename Traits::SizeType j) const
>>> {
>>> - return j*k+i;
>>> + return (j%s)+(j/s)*k*s+i*s;
>>> }
>>>
>>> //------------------------------
>>> @@ -1478,21 +1489,26 @@
>>> private:
>>> void setup ()
>>> {
>>> - Dune::dinfo << "power grid function space(blockwise version):"
>>> + Dune::dinfo << "PowerGridFunctionSpace(blockwise version):"
>>> << std::endl;
>>> Dune::dinfo << "( ";
>>> offset[0] = 0;
>>> maxlocalsize = 0;
>>> for (int i=0; i<k; i++)
>>> - {
>>> - childSize[i] = this->getChild(i).globalSize();
>>> - Dune::dinfo << childSize[i] << " ";
>>> - offset[i+1] = offset[i]+childSize[i];
>>> - maxlocalsize += this->getChild(i).maxLocalSize();
>>> - }
>>> + {
>>> + childSize[i] = this->getChild(i).globalSize();
>>> + offset[i+1] = offset[i]+childSize[i];
>>> + Dune::dinfo << childSize[i] << "[" << offset[i] << "] ";
>>> + maxlocalsize += this->getChild(i).maxLocalSize();
>>> + }
>>> Dune::dinfo << ") total size = " << offset[k]
>>> << " max local size = " << maxlocalsize
>>> << std::endl;
>>> + /* check the local block size */
>>> + if (this->getChild(0).maxLocalSize()%s != 0)
>>> + DUNE_THROW(Exception,
>>> + "number of DOFs (" << this->getChild(0).maxLocalSize() << ") per component "
>>> + "must be a multiple of the BlockSize (" << s << ")");
>>> for (int i=1; i<k; i++)
>>> if (childSize[i]!=childSize[0])
>>> DUNE_THROW(Exception, "components must be of equal size");
>>> @@ -2181,7 +2197,7 @@
>>> private:
>>> void setup ()
>>> {
>>> - Dune::dinfo << "composite grid function space(lexicographic version):"
>>> + Dune::dinfo << "CompositeGridFunctionSpace(lexicographic version):"
>>> << std::endl;
>>>
>>> CompositeGridFunctionSpaceBaseVisitChildMetaProgram<CompositeGridFunctionSpaceBase,BaseT::CHILDREN,0>::
>>> @@ -2216,14 +2232,15 @@
>>> // P is the ordering parameter
>>> // Ti are all grid function spaces
>>> template<typename T0, typename T1, typename T2, typename T3,
>>> - typename T4, typename T5, typename T6,
>>> - typename T7, typename T8>
>>> - class CompositeGridFunctionSpaceBase<GridFunctionSpaceBlockwiseMapper,T0,T1,T2,T3,T4,T5,T6,T7,T8>
>>> + typename T4, typename T5, typename T6, typename T7, typename T8,
>>> + int s0, int s1, int s2, int s3, int s4, int s5, int s6, int s7, int s8, int s9>
>>> + class CompositeGridFunctionSpaceBase<GridFunctionSpaceComponentBlockwiseMapper<s0,s1,s2,s3,s4,s5,s6,s7,s8,s9>,
>>> + T0,T1,T2,T3,T4,T5,T6,T7,T8>
>>> : public CompositeNode<CountingPointerStoragePolicy,T0,T1,T2,T3,T4,T5,T6,T7,T8>,
>>> public Countable
>>> {
>>> - friend class CompositeGridFunctionSpace<GridFunctionSpaceBlockwiseMapper,
>>> - T0,T1,T2,T3,T4,T5,T6,T7,T8>; // for setup
>>> + typedef GridFunctionSpaceComponentBlockwiseMapper<s0,s1,s2,s3,s4,s5,s6,s7,s8,s9> BlockwiseMapper;
>>> + friend class CompositeGridFunctionSpace<BlockwiseMapper, T0,T1,T2,T3,T4,T5,T6,T7,T8>; // for setup
>>>
>>> typedef CompositeNode<CountingPointerStoragePolicy,T0,T1,T2,T3,T4,T5,T6,T7,T8> BaseT;
>>>
>>> @@ -2231,7 +2248,7 @@
>>> //! export traits class
>>> typedef PowerCompositeGridFunctionSpaceTraits<typename T0::Traits::GridViewType,
>>> typename T0::Traits::BackendType,
>>> - GridFunctionSpaceBlockwiseMapper,
>>> + BlockwiseMapper,
>>> NonEmptyChilds<T0,T1,T2,T3,T4,T5,
>>> T6,T7,T8>::value>
>>> Traits;
>>> @@ -2301,7 +2318,15 @@
>>> template<int i>
>>> typename Traits::SizeType subMap (typename Traits::SizeType j) const
>>> {
>>> - return j*BaseT::CHILDREN+i;
>>> + // make the block sizes and offsets available in an array
>>> + static const int blockSize[] = { s0, s1, s2, s3, s4, s5, s6, s7, s8, s9 };
>>> + static const int blockOffset[] = { 0, s0, s0+s1, s0+s1+s2, s0+s1+s2+s3, s0+s1+s2+s3+s4,
>>> + s0+s1+s2+s3+s4+s5, s0+s1+s2+s3+s4+s5+s6, s0+s1+s2+s3+s4+s5+s6+s7,
>>> + s0+s1+s2+s3+s4+s5+s6+s7+s8, s0+s1+s2+s3+s4+s5+s6+s7+s8+s9 };
>>> + return (j%BlockwiseMapper::size[i])
>>> + +(j/BlockwiseMapper::size[i])*BlockwiseMapper::offset[BaseT::CHILDREN]
>>> + +BlockwiseMapper::offset[i];
>>> + return (j%blockSize[i])+(j/blockSize[i])*blockOffset[BaseT::CHILDREN]+blockOffset[i];
>>> }
>>>
>>> //------------------------------
>>> @@ -2357,15 +2382,24 @@
>>> private:
>>> void setup ()
>>> {
>>> - Dune::dinfo << "composite grid function space(blockwise version):"
>>> + Dune::dinfo << "CompositeGridFunctionSpace(blockwise version):"
>>> << std::endl;
>>>
>>> CompositeGridFunctionSpaceBaseVisitChildMetaProgram<CompositeGridFunctionSpaceBase,BaseT::CHILDREN,0>::
>>> setup(*this,childGlobalSize,childLocalSize);
>>>
>>> + // make the block sizes available in an array
>>> + static const int blockSize[] = { s0, s1, s2, s3, s4, s5, s6, s7, s8, s9 };
>>> + // check for compatible sizes
>>> for (int i=1; i<BaseT::CHILDREN; i++)
>>> - if (childGlobalSize[i]!=childGlobalSize[0])
>>> + {
>>> + if (childLocalSize[i]%blockSize[i]!=0)
>>> + DUNE_THROW(Exception,
>>> + "number of DOFs (" << childLocalSize[i] << ") per component "
>>> + "must be a multiple of the BlockSize (" << blockSize[i] << ")");
>>> + if (childGlobalSize[i]/blockSize[i]!=childGlobalSize[0]/blockSize[0])
>>> DUNE_THROW(Exception, "components must be of equal size");
>>> + }
>>>
>>> Dune::dinfo << "( ";
>>> offset[0] = 0;
>>> @@ -2389,7 +2423,6 @@
>>> mutable std::vector<typename Traits::SizeType> childglobal;
>>> };
>>>
>>> -
>>> //! \addtogroup GridFunctionSpace \{
>>>
>>> //! \brief grid function space composed of other grid function spaces
>>> @@ -3199,7 +3232,7 @@
>>> GridFunctionSubSpace (const GFS& gfs)
>>> : GridFunctionSubSpaceIntermediateBase<GFS,k,GFS::isLeaf>(gfs)
>>> {
>>> - Dune::dinfo << "grid function subspace:" << std::endl;
>>> + Dune::dinfo << "GridFunctionSubSpace:" << std::endl;
>>> Dune::dinfo << "root space size = " << gfs.globalSize()
>>> << " max local size = " << this->maxLocalSize()
>>> << std::endl;
>>>
>>> _______________________________________________
>>> dune-pdelab-commit mailing list
>>> dune-pdelab-commit at dune-project.org
>>> http://lists.dune-project.org/mailman/listinfo/dune-pdelab-commit
>>>
>>>
>>>
>>
>
--
_____________________________________________________________________
Bernd Flemisch phone: +49 711 685 69162
IWS, Universität Stuttgart fax: +49 711 685 60430
Pfaffenwaldring 61 email: bernd at iws.uni-stuttgart.de
D-70569 Stuttgart url: www.hydrosys.uni-stuttgart.de
_____________________________________________________________________
More information about the dune-pdelab
mailing list