10 #include <dune/common/fvector.hh>
11 #include <dune/common/fmatrix.hh>
12 #include <dune/common/bigunsignedint.hh>
13 #include <dune/common/parallel/collectivecommunication.hh>
14 #include <dune/common/reservedvector.hh>
15 #include <dune/geometry/genericgeometry/topologytypes.hh>
16 #include <dune/geometry/axisalignedcubegeometry.hh>
19 #include <dune/grid/sgrid/numbering.hh>
43 template<
int dim,
int dimworld,
class Gr
idImp>
class SGeometry;
44 template<
int codim,
int dim,
class Gr
idImp>
class SEntity;
46 template<
int codim,
class Gr
idImp>
class SEntitySeed;
47 template<
int codim, PartitionIteratorType,
class Gr
idImp>
class SLevelIterator;
48 template<
int dim,
int dimworld,
class ctype>
class SGrid;
53 namespace FacadeOptions
56 template<
int dim,
int dimworld,
class ctype,
int mydim,
int cdim>
60 static const bool v =
false;
63 template<
int dim,
int dimworld,
class ctype,
int mydim,
int cdim>
67 static const bool v =
false;
103 template<
int mydim,
int cdim,
class Gr
idImp>
105 :
public AxisAlignedCubeGeometry<typename GridImp::ctype,mydim,cdim>
109 typedef typename GridImp::ctype
ctype;
118 void make (
const FieldVector<ctype,cdim>& lower,
119 const FieldMatrix<ctype,mydim,cdim>& A)
123 static_cast< AxisAlignedCubeGeometry<ctype,mydim,cdim> &
>( *this ) = AxisAlignedCubeGeometry<ctype,mydim,cdim>(lower);
128 FieldVector<ctype, cdim> upper = lower;
129 for (
int i=0; i<mydim; i++)
133 std::bitset<cdim> axes;
135 for (
size_t i=0; i<cdim; i++)
136 if ((upper[i] - lower[i]) > 1e-10)
140 static_cast< AxisAlignedCubeGeometry<ctype,mydim,cdim> &
>( *this ) = AxisAlignedCubeGeometry<ctype,mydim,cdim>(lower, upper, axes);
145 : AxisAlignedCubeGeometry<
ctype,mydim,cdim>(FieldVector<
ctype,cdim>(0),FieldVector<
ctype,cdim>(0))
155 template<
int codim,
int dim,
class Gr
idImp,
template<
int,
int,
class>
class EntityImp>
161 enum { dimworld = GridImp::dimensionworld };
163 typedef typename GridImp::Traits::template Codim< codim >::GeometryImpl GeometryImpl;
166 typedef typename GridImp::ctype
ctype;
167 typedef typename GridImp::template Codim<codim>::Geometry
Geometry;
228 void make (GridImp* _grid,
int _l,
int _id);
231 void make (
int _l,
int _id);
239 return grid->persistentIndex(
l, codim,
z);
253 if (codim<dim || l==grid->maxLevel())
258 array<int,dim> coord;
259 for (
int k=0; k<dim; k++)
260 coord[k] =
z[k]*(1<<(
grid->maxLevel()-
l));
263 return grid->n(
grid->maxLevel(),coord);
269 DUNE_THROW(NotImplemented,
"subIndex for entities with codimension > 0 is not implemented");
276 DUNE_THROW(NotImplemented,
"subIndex for entities with codimension > 0 is not implemented");
296 template<
int codim,
int dim,
class Gr
idImp>
334 template<
int dim,
class Gr
idImp>
337 enum { dimworld = GridImp::dimensionworld };
344 typedef typename GridImp::Traits::template Codim< 0 >::GeometryImpl GeometryImpl;
345 typedef typename GridImp::Traits::template Codim< 0 >::LocalGeometryImpl LocalGeometryImpl;
351 typedef typename GridImp::ctype
ctype;
352 typedef typename GridImp::template Codim<0>::Geometry
Geometry;
371 template<
int cc>
int count ()
const;
384 return (this->
grid)->n(this->
l, this->
grid->subz(this->z,i,codim));
394 assert(this->
l == this->
grid->maxLevel());
396 return (this->
grid)->n(this->
l, this->
grid->subz(this->z,i,codim));
404 return this->
grid->persistentIndex(this->
l, codim, this->
grid->subz(this->z,i,codim));
430 bool hasFather ()
const
432 return (this->
level()>0);
438 return ( this->
grid->maxLevel() == this->
level() );
452 LocalGeometry geometryInFather ()
const;
467 SEntity (GridImp* _grid,
int _l,
int _index) :
478 void make (GridImp* _grid,
int _l,
int _id)
481 built_father =
false;
488 built_father =
false;
495 mutable bool built_father;
496 mutable int father_index;
497 mutable LocalGeometryImpl in_father_local;
498 void make_father()
const;
519 template<
class Gr
idImp>
524 enum { dim = GridImp::dimension };
525 enum { dimworld = GridImp::dimensionworld };
532 typedef typename GridImp::template Codim<0>::Entity
Entity;
533 typedef typename GridImp::ctype
ctype;
546 int _maxLevel,
bool makeend) :
554 orig_l = this->
entity().level();
555 orig_index = _grid->getRealImplementation(this->
entity()).compressedIndex();
559 stack.push(originalElement);
565 push_sons(orig_l,orig_index);
573 int orig_l, orig_index;
576 std::stack<SHierarchicStackElem, Dune::ReservedVector<SHierarchicStackElem,GridImp::MAXL> > stack;
578 void push_sons (
int level,
int fatherid);
588 template<
class Gr
idImp>
591 enum { dim=GridImp::dimension };
592 enum { dimworld=GridImp::dimensionworld };
594 typedef typename GridImp::Traits::template Codim< 1 >::GeometryImpl GeometryImpl;
595 typedef typename GridImp::Traits::template Codim< 1 >::LocalGeometryImpl LocalGeometryImpl;
600 typedef typename GridImp::template Codim<0>::Entity
Entity;
602 typedef typename GridImp::template Codim<1>::Geometry
Geometry;
611 typedef typename GridImp::ctype
ctype;
642 return grid->boundarySegmentIndex(
self.level(), count, zred);
676 self(*_self), ne(self), grid(_grid),
677 partition(_grid->partition(grid->getRealImplementation(ne).l,_self->z)),
678 zred(_grid->compress(grid->getRealImplementation(ne).l,_self->z)),
686 self(other.self), ne(other.ne), grid(other.grid),
687 partition(other.partition), zred(other.zred),
688 count(other.count), valid_count(other.valid_count),
689 valid_nb(other.valid_nb), is_on_boundary(other.is_on_boundary),
690 built_intersections(false),
698 assert(grid == other.grid);
703 partition = other.partition;
706 valid_count = other.valid_count;
707 valid_nb = other.valid_nb;
708 is_on_boundary = other.is_on_boundary;
711 built_intersections =
false;
717 void make (
int _count)
const;
718 void makeintersections ()
const;
721 const GridImp * grid;
725 mutable bool valid_count;
726 mutable bool valid_nb;
727 mutable bool is_on_boundary;
728 mutable bool built_intersections;
729 mutable LocalGeometryImpl is_self_local;
730 mutable GeometryImpl is_global;
731 mutable LocalGeometryImpl is_nb_local;
735 template<
class Gr
idImp>
738 enum { dim=GridImp::dimension };
739 enum { dimworld=GridImp::dimensionworld };
741 typedef typename GridImp::template Codim<0>::Entity
Entity;
743 typedef typename GridImp::template Codim<1>::Geometry
Geometry;
753 typedef typename GridImp::ctype
ctype;
757 return is.boundary();
763 return is.boundaryId();
769 return is.boundarySegmentIndex();
775 return is.neighbor();
799 return is.geometryInInside();
805 return is.geometryInOutside();
811 return is.geometry();
823 return is.indexInInside();
829 return is.indexInOutside();
842 n *= is.geometry().integrationElement(local);
855 FieldVector<ctype, dimworld> normal(0.0);
856 normal[is.count/2] = (is.count%2) ? 1.0 : -1.0;
864 #ifndef DOXYGEN // doxygen can't handle this recursive usage
880 while(! this->empty())
891 template<
int codim,
class Gr
idImp>
894 enum { dim = GridImp::dimension };
898 typedef typename GridImp::template Codim<codim>::Entity
Entity;
987 grid->getRealImplementation(*e).make(_grid, _l,_id);
1000 template<
int codim,
class Gr
idImp>
1003 enum { dim = GridImp::dimension };
1014 _l(l), _index(index)
1024 int index ()
const {
return this->_index; }
1036 template<
int codim, PartitionIteratorType pitype,
class Gr
idImp>
1041 enum { dim = GridImp::dimension };
1047 typedef typename GridImp::template Codim<codim>::Entity
Entity;
1065 template<
class Gr
idImp>
1071 enum { dim = GridImp::dimension };
1088 int index (
const typename GridImp::Traits::template Codim<cd>::Entity& e)
const
1090 return grid.getRealImplementation(e).compressedIndex();
1094 int subIndex (
const typename GridImp::Traits::template Codim< cc >::Entity &e,
1095 int i,
unsigned int codim )
const
1098 return grid.getRealImplementation(e).subCompressedIndex(codim, i);
1100 DUNE_THROW( NotImplemented,
"subIndex for higher codimension entity not implemented for SGrid." );
1104 template<
class EntityType >
1107 return (e.level() == level);
1113 return grid.size( level, type );
1119 return grid.size( level, codim );
1123 const std::vector<GeometryType>&
geomTypes (
int codim)
const
1125 return mytypes[codim];
1129 const GridImp& grid;
1131 std::vector<GeometryType> mytypes[GridImp::dimension+1];
1143 template<
class Gr
idImp>
1145 public IdSet<GridImp,SGridGlobalIdSet<GridImp>, typename remove_const<GridImp>::type::PersistentIndexType>
1160 typedef typename remove_const<GridImp>::type::PersistentIndexType
IdType;
1168 IdType id (
const typename remove_const<GridImp>::type::Traits::template Codim<cd>::Entity& e)
const
1170 return GridImp::getRealImplementation(e).persistentIndex();
1178 IdType subId (
const typename remove_const< GridImp >::type::Traits::template Codim< 0 >::Entity &e,
1179 int i,
unsigned int codim )
const
1181 return GridImp::getRealImplementation(e).subPersistentIndex(codim, i);
1186 template<
int dim,
int dimworld,
class ctype>
1201 bigunsignedint<dim*sgrid_dim_bits+sgrid_level_bits+sgrid_codim_bits>,
1203 bigunsignedint<dim*sgrid_dim_bits+sgrid_level_bits+sgrid_codim_bits>,
1204 CollectiveCommunication<Dune::SGrid<dim,dimworld,ctype> >,
1266 template<
int dim,
int dimworld,
typename _ctype = sgr
id_ctype>
1295 SGrid (
const int *
const N_,
const ctype *
const H_);
1304 SGrid (
const int *
const N_,
const ctype *
const L_,
const ctype *
const H_);
1315 SGrid (FieldVector<int,dim> N_, FieldVector<ctype,dim> L_, FieldVector<ctype,dim> H_);
1328 template<
int cd, PartitionIteratorType pitype>
1332 template<
int cd, PartitionIteratorType pitype>
1339 return lbegin<cd,All_Partition>(level);
1346 return lend<cd,All_Partition>(level);
1350 template<
int cd, PartitionIteratorType pitype>
1351 typename Traits::template Codim<cd>::template Partition<pitype>::LeafIterator
leafbegin ()
const;
1354 template<
int cd, PartitionIteratorType pitype>
1355 typename Traits::template Codim<cd>::template Partition<pitype>::LeafIterator
leafend ()
const;
1361 return leafbegin<cd,All_Partition>();
1368 return leafend<cd,All_Partition>();
1372 template <
typename Seed>
1373 typename Traits::template Codim<Seed::codimension>::EntityPointer
1376 enum { codim = Seed::codimension };
1395 template<
class T,
template<
class>
class P,
int codim>
1403 int size (
int level,
int codim)
const;
1414 return (type.isCube()) ?
size(level,dim-type.dim()) : 0;
1426 return boundarysize;
1450 const array<int, dim>&
dims(
int level)
const {
1474 return theglobalidset;
1479 return theglobalidset;
1484 assert(level>=0 && level<=
maxLevel());
1485 return *(indexsets[level]);
1490 return *indexsets.back();
1497 template<
class DataHandle>
1501 template<
class DataHandle>
1505 const CollectiveCommunication<SGrid>&
comm ()
const
1554 template<int codim_, int dim_, class GridImp_, template<int,int,class> class EntityImp_>
1557 template<int codim_, class GridImp_>
1560 template<int codim_, int dim_, class GridImp_, template<int,int,class> class EntityImp_>
1564 FieldVector<ctype, dimworld> pos (int level, array<int,dim>& z) const;
1567 int calc_codim (int level, const array<int,dim>& z) const;
1570 int n (int level, const array<int,dim>& z) const;
1573 array<int,dim> z (int level, int i, int codim) const;
1576 array<int,dim> subz (const array<int,dim> & z, int i, int codim) const;
1579 array<int,dim> compress (int level, const array<int,dim>& z) const;
1582 array<int,dim> expand (int level, const array<int,dim>& r, int b) const;
1587 int partition (int level, const array<int,dim>& z) const;
1590 bool exists (int level, const array<int,dim>& zred) const;
1593 int boundarySegmentIndex (int l, int face, const array<int,dim> & zentity) const
1595 array<int,dim-1> zface;
1599 for (
int i=0; i<dir; i++) zface[i] = zentity[i]/(1<<l);
1600 for (
int i=dir+1; i<dim; i++) zface[i-1] = zentity[i]/(1<<l);
1601 zface = boundarymapper[dir].expand(zface, 0);
1603 int index = boundarymapper[dir].n(zface);
1605 for (
int i=0; i<dir; i++)
1606 index += 2*boundarymapper[i].elements(0);
1607 index += side*boundarymapper[dir].elements(0);
1612 PersistentIndexType persistentIndex (
int l,
int codim,
const array<int,dim> & zentity)
const
1625 for (
int i=dim-1; i>=0; i--)
1637 int trailing = 1000;
1638 for (
int i=0; i<dim; i++)
1642 for (
int j=0; j<l; j++)
1643 if (zentity[i]&(1<<(j+1)))
1647 trailing =
std::min(trailing,zeros);
1651 int level = l-trailing;
1661 for (
int i=dim-1; i>=0; i--)
1673 SGrid & operator = (
const SGrid &) {
return *
this; }
1675 void makeSGrid (
const array<int,dim>& N_,
const FieldVector<ctype, dim>& L_,
const FieldVector<ctype, dim>& H_);
1680 CollectiveCommunication<SGrid> ccobj;
1682 ReservedVector<SGridLevelIndexSet<const SGrid<dim,dimworld> >*,
MAXL> indexsets;
1683 SGridGlobalIdSet<const SGrid<dim,dimworld> > theglobalidset;
1686 FieldVector<ctype, dim> low;
1687 FieldVector<ctype, dim> H;
1688 std::vector<array<int,dim> > N;
1689 std::vector<FieldVector<ctype, dim> > h;
1690 mutable CubeMapper<dim> *mapper;
1693 array<CubeMapper<dim-1>, dim> boundarymapper;
1697 namespace Capabilities
1711 template<
int dim,
int dimw>
1714 static const bool v =
true;
1715 static const unsigned int topologyId = GenericGeometry :: CubeTopology< dim > :: type :: id ;
1721 template<
int dim,
int dimw>
1724 static const bool v =
true;
1730 template<
int dim,
int dimw,
int cdim>
1733 static const bool v =
true;
1739 template<
int dim,
int dimw>
1742 static const bool v =
true;
1748 template<
int dim,
int dimw>
1751 static const bool v =
true;
1758 #include "sgrid/sgrid.cc"