/usr/include/dune/grid/albertagrid/gridfactory.hh is in libdune-grid-dev 2.4.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 | // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ALBERTA_GRIDFACTORY_HH
#define DUNE_ALBERTA_GRIDFACTORY_HH
/** \file
* \author Martin Nolte
* \brief specialization of the generic GridFactory for AlbertaGrid
*/
#include <algorithm>
#include <limits>
#include <map>
#include <dune/common/array.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/grid/common/gridfactory.hh>
#include <dune/grid/utility/grapedataioformattypes.hh>
#include <dune/grid/albertagrid/agrid.hh>
#if HAVE_ALBERTA
namespace Dune
{
/** \brief specialization of the generic GridFactory for AlbertaGrid
*
* \ingroup GridFactory
*
* The GridFactory for AlbertaGrid adds some extensions to the standard
* GridFactoryInterface. It provides the following additional features:
* - It allows to set boundary ids via insertBoundary. For ALBERTA 1.2,
* these boundary ids are ignored, though.
* - For ALBERTA 3.0 and above, you can add face transformation to identify
* faces. This allows the construction of periodic grids.
* - The grid can be written in ALBERTA's native format for macro
* triangulations via write (both ASCII and XDR format are supported).
* - On grid creation, a name can be assigned to the grid.
* - On grid creation, the grid can be reordered such that ALBERTA uses
* the longest edge as refinement edge during recursive bisection.
* .
*/
template< int dim, int dimworld >
class GridFactory< AlbertaGrid< dim, dimworld > >
: public GridFactoryInterface< AlbertaGrid< dim, dimworld > >
{
typedef GridFactory< AlbertaGrid< dim, dimworld > > This;
public:
//! type of grid this factory is for
typedef AlbertaGrid< dim, dimworld > Grid;
//! type of (scalar) coordinates
typedef typename Grid::ctype ctype;
//! dimension of the grid
static const int dimension = Grid::dimension;
//! dimension of the world
static const int dimensionworld = Grid::dimensionworld;
//! type of vector for world coordinates
typedef FieldVector< ctype, dimensionworld > WorldVector;
//! type of matrix from world coordinates to world coordinates
typedef FieldMatrix< ctype, dimensionworld, dimensionworld > WorldMatrix;
typedef DuneBoundaryProjection< dimensionworld > DuneProjection;
typedef Dune::shared_ptr< const DuneProjection > DuneProjectionPtr;
typedef Dune::BoundarySegment< dimension, dimensionworld > BoundarySegment;
template< int codim >
struct Codim
{
typedef typename Grid::template Codim< codim >::Entity Entity;
};
private:
typedef Dune::BoundarySegmentWrapper< dimension, dimensionworld > BoundarySegmentWrapper;
static const int numVertices
= Alberta::NumSubEntities< dimension, dimension >::value;
typedef Alberta::MacroElement< dimension > MacroElement;
typedef Alberta::ElementInfo< dimension > ElementInfo;
typedef Alberta::MacroData< dimension > MacroData;
typedef Alberta::NumberingMap< dimension, Alberta::Dune2AlbertaNumbering > NumberingMap;
typedef Alberta::DuneBoundaryProjection< dimension > Projection;
typedef array< unsigned int, dimension > FaceId;
typedef std::map< FaceId, size_t > BoundaryMap;
class ProjectionFactory;
public:
//! are boundary ids supported by this factory?
static const bool supportsBoundaryIds = true;
//! is the factory able to create periodic meshes?
static const bool supportPeriodicity = MacroData::supportPeriodicity;
/** default constructor */
GridFactory ()
: globalProjection_( (const DuneProjection *) 0 )
{
macroData_.create();
}
virtual ~GridFactory ();
/** \brief insert a vertex into the macro grid
*
* \param[in] pos position of the vertex (in world coordinates)
*/
virtual void insertVertex ( const WorldVector &pos )
{
macroData_.insertVertex( pos );
}
/** \brief insert an element into the macro grid
*
* \param[in] type GeometryType of the new element
* \param[in] vertices indices of the element vertices (in DUNE numbering)
*/
virtual void insertElement ( const GeometryType &type,
const std::vector< unsigned int > &vertices )
{
if( (int)type.dim() != dimension )
DUNE_THROW( AlbertaError, "Inserting element of wrong dimension: " << type.dim() );
if( !type.isSimplex() )
DUNE_THROW( AlbertaError, "Alberta supports only simplices." );
if( vertices.size() != (size_t)numVertices )
DUNE_THROW( AlbertaError, "Wrong number of vertices passed: " << vertices.size() << "." );
int array[ numVertices ];
for( int i = 0; i < numVertices; ++i )
array[ i ] = vertices[ numberingMap_.alberta2dune( dimension, i ) ];
macroData_.insertElement( array );
}
/** \brief mark a face as boundary (and assign a boundary id)
*
* \internal
*
* \param[in] element index of the element, the face belongs to
* \param[in] face local number of the face within the element
* \param[in] id boundary id to assign to the face
*
* \note ALBERTA supports only boundary id in the range 1,...,127.
*/
virtual void insertBoundary ( int element, int face, int id )
{
if( (id <= 0) || (id > 127) )
DUNE_THROW( AlbertaError, "Invalid boundary id: " << id << "." );
macroData_.boundaryId( element, numberingMap_.dune2alberta( 1, face ) ) = id;
}
/** \brief insert a boundary projection into the macro grid
*
* \internal
*
* \param[in] type geometry type of boundary face
* \param[in] vertices vertices of the boundary face
* \param[in] projection boundary projection
*
* \note The grid takes control of the projection object.
*/
virtual void
insertBoundaryProjection ( const GeometryType &type,
const std::vector< unsigned int > &vertices,
const DuneProjection *projection )
{
if( (int)type.dim() != dimension-1 )
DUNE_THROW( AlbertaError, "Inserting boundary face of wrong dimension: " << type.dim() );
if( !type.isSimplex() )
DUNE_THROW( AlbertaError, "Alberta supports only simplices." );
FaceId faceId;
if( vertices.size() != faceId.size() )
DUNE_THROW( AlbertaError, "Wrong number of face vertices passed: " << vertices.size() << "." );
for( size_t i = 0; i < faceId.size(); ++i )
faceId[ i ] = vertices[ i ];
std::sort( faceId.begin(), faceId.end() );
typedef std::pair< typename BoundaryMap::iterator, bool > InsertResult;
const InsertResult result = boundaryMap_.insert( std::make_pair( faceId, boundaryProjections_.size() ) );
if( !result.second )
DUNE_THROW( GridError, "Only one boundary projection can be attached to a face." );
boundaryProjections_.push_back( DuneProjectionPtr( projection ) );
}
/** \brief insert a global (boundary) projection into the macro grid
*
* \internal
*
* \param[in] projection global (boundary) projection
*
* \note The grid takes control of the projection object.
*/
virtual void insertBoundaryProjection ( const DuneProjection *projection )
{
if( globalProjection_ )
DUNE_THROW( GridError, "Only one global boundary projection can be attached to a grid." );
globalProjection_ = DuneProjectionPtr( projection );
}
/** \brief insert a boundary segment into the macro grid
*
* Only influences the ordering of the boundary segments
* \param[in] vertices vertex indices of boundary face
*/
virtual void
insertBoundarySegment ( const std::vector< unsigned int >& vertices )
{
typedef typename GenericGeometry::SimplexTopology< dimension-1 >::type Topology;
insertBoundaryProjection( GeometryType( Topology() ), vertices, 0 );
}
/** \brief insert a shaped boundary segment into the macro grid
*
* \param[in] vertices vertex indices of boundary face
* \param[in] boundarySegment geometric realization of shaped boundary
*/
virtual void
insertBoundarySegment ( const std::vector< unsigned int > &vertices,
const shared_ptr< BoundarySegment > &boundarySegment )
{
const ReferenceElement< ctype, dimension-1 > &refSimplex
= ReferenceElements< ctype, dimension-1 >::simplex();
if( !boundarySegment )
DUNE_THROW( GridError, "Trying to insert null as a boundary segment." );
if( (int)vertices.size() != refSimplex.size( dimension-1 ) )
DUNE_THROW( GridError, "Wrong number of face vertices passed: " << vertices.size() << "." );
std::vector< WorldVector > coords( refSimplex.size( dimension-1 ) );
for( int i = 0; i < dimension; ++i )
{
Alberta::GlobalVector &x = macroData_.vertex( vertices[ i ] );
for( int j = 0; j < dimensionworld; ++j )
coords[ i ][ j ] = x[ j ];
if( ((*boundarySegment)( refSimplex.position( i, dimension-1 ) ) - coords[ i ]).two_norm() > 1e-6 )
DUNE_THROW( GridError, "Boundary segment does not interpolate the corners." );
}
const GeometryType gt = refSimplex.type( 0, 0 );
const DuneProjection *prj = new BoundarySegmentWrapper( gt, coords, boundarySegment );
insertBoundaryProjection( gt, vertices, prj );
}
/** \brief add a face transformation (for periodic identification)
*
* A face transformation is an affine mapping T from world coordinates
* to world coordinates. ALBERTA periodically identifies two faces f and g
* if T( f ) = g or T( g ) = f.
*
* \param[in] matrix matrix describing the linear part of T
* \param[in] shift vector describing T( 0 )
*
* \note ALBERTA requires the matrix to be orthogonal.
*
* \note ALBERTA automatically adds the inverse transformation.
*/
void insertFaceTransformation ( const WorldMatrix &matrix, const WorldVector &shift );
/** \brief mark the longest edge as refinemet edge
*
* Marking the longest edge avoids cycles in the recursive bisection
* algorithm, if the longest edge of each element is unique. It also
* makes sure the angles degenerate least. It can, hoowever, produce
* more nonlocal refinements than necessary. Therefore this feature is
* disabled by default.
*/
void markLongestEdge ()
{
macroData_.markLongestEdge();
}
/** \brief finalize grid creation and hand over the grid
*
* This version of createGrid is original to the AlbertaGrid grid factroy,
* allowing to specity a grid name.
*
* \returns a pointer to the newly created grid
*
* \note The caller takes responsibility of freeing the memory allocated
* for the grid.
* \note ALBERTA's grid factory provides a static method for freeing the
* grid (destroyGrid).
*/
Grid *createGrid ()
{
macroData_.finalize();
if( macroData_.elementCount() == 0 )
DUNE_THROW( GridError, "Cannot create empty AlbertaGrid." );
if( dimension < 3 )
macroData_.setOrientation( Alberta::Real( 1 ) );
assert( macroData_.checkNeighbors() );
macroData_.checkCycles();
ProjectionFactory projectionFactory( *this );
return new Grid( macroData_, projectionFactory );
}
/** \brief destroy a grid previously obtain from this factory
*
* \param[in] grid pointer to the grid to destroy
*/
static void destroyGrid ( Grid *grid )
{
delete grid;
}
/** \brief write out the macro triangulation in native grid file format
*
* \tparam type type of file to write (either ascii or xdr)
*
* \param[in] filename name of the file to write to
*
* \returns \c true on success
*/
template< GrapeIOFileFormatType type >
bool write ( const std::string &filename )
{
static_assert( type != pgm, "AlbertaGridFactory: writing pgm format is not supported." );
macroData_.finalize();
if( dimension < 3 )
macroData_.setOrientation( Alberta::Real( 1 ) );
assert( macroData_.checkNeighbors() );
return macroData_.write( filename, (type == xdr) );
}
/** \brief write out the macro triangulation in native grid file format
*
* The grid is written in human readable form (ascii).
*
* \param[in] filename name of the file to write to
*
* \returns \c true on success
*/
virtual bool write ( const std::string &filename )
{
return write< ascii >( filename );
}
virtual unsigned int
insertionIndex ( const typename Codim< 0 >::Entity &entity ) const
{
return insertionIndex( Grid::getRealImplementation( entity ).elementInfo() );
}
virtual unsigned int
insertionIndex ( const typename Codim< dimension >::Entity &entity ) const
{
const int elIndex = insertionIndex( Grid::getRealImplementation( entity ).elementInfo() );
const typename MacroData::ElementId &elementId = macroData_.element( elIndex );
return elementId[ Grid::getRealImplementation( entity ).subEntity() ];
}
virtual unsigned int
insertionIndex ( const typename Grid::LeafIntersection &intersection ) const
{
const Grid &grid = Grid::getRealImplementation( intersection ).grid();
const ElementInfo &elementInfo = Grid::getRealImplementation( intersection ).elementInfo();
const int face = grid.generic2alberta( 1, intersection.indexInInside() );
return insertionIndex( elementInfo, face );
}
virtual bool
wasInserted ( const typename Grid::LeafIntersection &intersection ) const
{
return (insertionIndex( intersection ) < std::numeric_limits< unsigned int >::max());
}
private:
unsigned int insertionIndex ( const ElementInfo &elementInfo ) const;
unsigned int insertionIndex ( const ElementInfo &elementInfo, const int face ) const;
FaceId faceId ( const ElementInfo &elementInfo, const int face ) const;
MacroData macroData_;
NumberingMap numberingMap_;
DuneProjectionPtr globalProjection_;
BoundaryMap boundaryMap_;
std::vector< DuneProjectionPtr > boundaryProjections_;
};
template< int dim, int dimworld >
GridFactory< AlbertaGrid< dim, dimworld > >::~GridFactory ()
{
macroData_.release();
}
template< int dim, int dimworld >
inline void
GridFactory< AlbertaGrid< dim, dimworld > >
::insertFaceTransformation ( const WorldMatrix &matrix, const WorldVector &shift )
{
// make sure the matrix is orthogonal
for( int i = 0; i < dimworld; ++i )
for( int j = 0; j < dimworld; ++j )
{
const ctype delta = (i == j ? ctype( 1 ) : ctype( 0 ));
const ctype epsilon = (8*dimworld)*std::numeric_limits< ctype >::epsilon();
if( std::abs( matrix[ i ] * matrix[ j ] - delta ) > epsilon )
{
DUNE_THROW( AlbertaError,
"Matrix of face transformation is not orthogonal." );
}
}
// copy matrix
Alberta::GlobalMatrix M;
for( int i = 0; i < dimworld; ++i )
for( int j = 0; j < dimworld; ++j )
M[ i ][ j ] = matrix[ i ][ j ];
// copy shift
Alberta::GlobalVector t;
for( int i = 0; i < dimworld; ++i )
t[ i ] = shift[ i ];
// insert into ALBERTA macro data
macroData_.insertWallTrafo( M, t );
}
template< int dim, int dimworld >
inline unsigned int
GridFactory< AlbertaGrid< dim, dimworld > >
::insertionIndex ( const ElementInfo &elementInfo ) const
{
const MacroElement ¯oElement = elementInfo.macroElement();
const unsigned int index = macroElement.index;
#ifndef NDEBUG
const typename MacroData::ElementId &elementId = macroData_.element( index );
for( int i = 0; i <= dimension; ++i )
{
const Alberta::GlobalVector &x = macroData_.vertex( elementId[ i ] );
const Alberta::GlobalVector &y = macroElement.coordinate( i );
for( int j = 0; j < dimensionworld; ++j )
{
if( x[ j ] != y[ j ] )
DUNE_THROW( GridError, "Vertex in macro element does not coincide with same vertex in macro data structure." );
}
}
#endif // #ifndef NDEBUG
return index;
}
template< int dim, int dimworld >
inline unsigned int
GridFactory< AlbertaGrid< dim, dimworld > >
::insertionIndex ( const ElementInfo &elementInfo, const int face ) const
{
typedef typename BoundaryMap::const_iterator Iterator;
const Iterator it = boundaryMap_.find( faceId( elementInfo, face ) );
if( it != boundaryMap_.end() )
return it->second;
else
return std::numeric_limits< unsigned int >::max();
}
template< int dim, int dimworld >
inline typename GridFactory< AlbertaGrid< dim, dimworld > >::FaceId
GridFactory< AlbertaGrid< dim, dimworld > >
::faceId ( const ElementInfo &elementInfo, const int face ) const
{
const unsigned int index = insertionIndex( elementInfo );
const typename MacroData::ElementId &elementId = macroData_.element( index );
FaceId faceId;
for( size_t i = 0; i < faceId.size(); ++i )
{
const int k = Alberta::MapVertices< dimension, 1 >::apply( face, i );
faceId[ i ] = elementId[ k ];
}
std::sort( faceId.begin(), faceId.end() );
return faceId;
}
// GridFactory::ProjectionFactory
// ------------------------------
template< int dim, int dimworld >
class GridFactory< AlbertaGrid< dim, dimworld > >::ProjectionFactory
: public Alberta::ProjectionFactory< Alberta::DuneBoundaryProjection< dim >, ProjectionFactory >
{
typedef ProjectionFactory This;
typedef Alberta::ProjectionFactory< Alberta::DuneBoundaryProjection< dim >, ProjectionFactory > Base;
typedef typename Dune::GridFactory< AlbertaGrid< dim, dimworld > > Factory;
public:
typedef typename Base::Projection Projection;
typedef typename Base::ElementInfo ElementInfo;
typedef typename Projection::Projection DuneProjection;
ProjectionFactory( const GridFactory &gridFactory )
: gridFactory_( gridFactory )
{}
bool hasProjection ( const ElementInfo &elementInfo, const int face ) const
{
if( gridFactory().globalProjection_ )
return true;
const unsigned int index = gridFactory().insertionIndex( elementInfo, face );
if( index < std::numeric_limits< unsigned int >::max() )
return bool( gridFactory().boundaryProjections_[ index ] );
else
return false;
}
bool hasProjection ( const ElementInfo &elementInfo ) const
{
return bool( gridFactory().globalProjection_ );
}
Projection projection ( const ElementInfo &elementInfo, const int face ) const
{
const unsigned int index = gridFactory().insertionIndex( elementInfo, face );
if( index < std::numeric_limits< unsigned int >::max() )
{
const DuneProjectionPtr &projection = gridFactory().boundaryProjections_[ index ];
if( projection )
return Projection( projection );
}
assert( gridFactory().globalProjection_ );
return Projection( gridFactory().globalProjection_ );
};
Projection projection ( const ElementInfo &elementInfo ) const
{
assert( gridFactory().globalProjection_ );
return Projection( gridFactory().globalProjection_ );
};
const GridFactory &gridFactory () const
{
return gridFactory_;
}
private:
const GridFactory &gridFactory_;
};
}
#endif // #if HAVE_ALBERTA
#endif // #ifndef DUNE_ALBERTA_GRIDFACTORY_HH
|