/usr/include/dune/geometry/genericgeometry/cachedmapping.hh is in libdune-geometry-dev 2.3.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 | // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_GEOMETRY_GENERICGEOMETRY_CACHED_MAPPING_HH
#define DUNE_GEOMETRY_GENERICGEOMETRY_CACHED_MAPPING_HH
#include <dune/geometry/type.hh>
#include <dune/geometry/genericgeometry/topologytypes.hh>
#include <dune/geometry/genericgeometry/referenceelements.hh>
#include <dune/geometry/genericgeometry/matrixhelper.hh>
#include <dune/geometry/genericgeometry/mapping.hh>
namespace Dune
{
namespace GenericGeometry
{
// Internal Forward Declarations
// -----------------------------
template< unsigned int, class >
class CachedJacobianTransposed;
template< unsigned int, class >
class CachedJacobianInverseTransposed;
// CachedStorage
// -------------
template< unsigned int dim, class GeometryTraits >
class CachedStorage
{
friend class CachedJacobianTransposed< dim, GeometryTraits >;
public:
static const unsigned int dimension = dim;
static const unsigned int dimWorld = GeometryTraits::dimWorld;
typedef MappingTraits< typename GeometryTraits::CoordTraits, dimension, dimWorld > Traits;
typedef typename GeometryTraits::Caching Caching;
CachedStorage ()
: affine( false ),
jacobianTransposedComputed( false ),
jacobianInverseTransposedComputed( false ),
integrationElementComputed( false )
{}
typename Traits::JacobianTransposedType jacobianTransposed;
typename Traits::JacobianType jacobianInverseTransposed;
typename Traits::FieldType integrationElement;
bool affine : 1;
bool jacobianTransposedComputed : 1; // = affine, if jacobian transposed was computed
bool jacobianInverseTransposedComputed : 1; // = affine, if jacobian inverse transposed was computed
bool integrationElementComputed : 1; // = affine, if integration element was computed
};
// CachedJacobianTranposed
// -----------------------
template< unsigned int dim, class GeometryTraits >
class CachedJacobianTransposed
{
friend class CachedJacobianInverseTransposed< dim, GeometryTraits >;
typedef CachedStorage< dim, GeometryTraits > Storage;
typedef typename Storage::Traits Traits;
typedef typename Traits::MatrixHelper MatrixHelper;
public:
typedef typename Traits::FieldType ctype;
static const int rows = Traits::dimension;
static const int cols = Traits::dimWorld;
typedef typename Traits::JacobianTransposedType FieldMatrix;
operator bool () const
{
return storage().jacobianTransposedComputed;
}
operator const FieldMatrix & () const
{
return storage().jacobianTransposed;
}
template< class X, class Y >
void mv ( const X &x, Y &y ) const
{
static_cast< const FieldMatrix & >( *this ).mv( x, y );
}
template< class X, class Y >
void mtv ( const X &x, Y &y ) const
{
static_cast< const FieldMatrix & >( *this ).mtv( x, y );
}
template< class X, class Y >
void umv ( const X &x, Y &y ) const
{
static_cast< const FieldMatrix & >( *this ).umv( x, y );
}
template< class X, class Y >
void umtv ( const X &x, Y &y ) const
{
static_cast< const FieldMatrix & >( *this ).umtv( x, y );
}
template< class X, class Y >
void mmv ( const X &x, Y &y ) const
{
static_cast< const FieldMatrix & >( *this ).mmv( x, y );
}
template< class X, class Y >
void mmtv ( const X &x, Y &y ) const
{
static_cast< const FieldMatrix & >( *this ).mmtv( x, y );
}
ctype det () const
{
if( !storage().integrationElementComputed )
{
storage().integrationElement = MatrixHelper::template sqrtDetAAT< rows, cols >( storage().jacobianTransposed );
storage().integrationElementComputed = storage().affine;
}
return storage().integrationElement;
}
private:
Storage &storage () const { return storage_; }
mutable Storage storage_;
};
// CachedJacobianInverseTransposed
// -------------------------------
template< unsigned int dim, class GeometryTraits >
class CachedJacobianInverseTransposed
{
template< class, class > friend class CachedMapping;
typedef CachedJacobianTransposed< dim, GeometryTraits > JacobianTransposed;
typedef typename JacobianTransposed::Storage Storage;
typedef typename JacobianTransposed::Traits Traits;
typedef typename Traits::MatrixHelper MatrixHelper;
public:
typedef typename Traits::FieldType ctype;
static const int rows = Traits::dimWorld;
static const int cols = Traits::dimension;
typedef typename Traits::JacobianType FieldMatrix;
operator bool () const
{
return storage().jacobianInverseTransposedComputed;
}
operator const FieldMatrix & () const
{
return storage().jacobianInverseTransposed;
}
template< class X, class Y >
void mv ( const X &x, Y &y ) const
{
static_cast< const FieldMatrix & >( *this ).mv( x, y );
}
template< class X, class Y >
void mtv ( const X &x, Y &y ) const
{
static_cast< const FieldMatrix & >( *this ).mtv( x, y );
}
template< class X, class Y >
void umv ( const X &x, Y &y ) const
{
static_cast< const FieldMatrix & >( *this ).umv( x, y );
}
template< class X, class Y >
void umtv ( const X &x, Y &y ) const
{
static_cast< const FieldMatrix & >( *this ).umtv( x, y );
}
template< class X, class Y >
void mmv ( const X &x, Y &y ) const
{
static_cast< const FieldMatrix & >( *this ).mmv( x, y );
}
template< class X, class Y >
void mmtv ( const X &x, Y &y ) const
{
static_cast< const FieldMatrix & >( *this ).mmtv( x, y );
}
ctype det () const
{
// integrationElement is always computed with jacobianInverseTransposed
return ctype( 1 ) / storage().integrationElement;
}
private:
JacobianTransposed &jacobianTransposed () { return jacobianTransposed_; }
const JacobianTransposed &jacobianTransposed () const { return jacobianTransposed_; }
Storage &storage () const { return jacobianTransposed().storage(); }
JacobianTransposed jacobianTransposed_;
};
// CachedMapping
// -------------
template< class Topology, class GeometryTraits >
class CachedMapping
{
typedef CachedMapping< Topology, GeometryTraits > This;
typedef typename GeometryTraits::template Mapping< Topology >::type MappingImpl;
public:
typedef MappingTraits< typename GeometryTraits::CoordTraits, Topology::dimension, GeometryTraits::dimWorld > Traits;
typedef GenericGeometry::Mapping< typename GeometryTraits::CoordTraits, Topology, GeometryTraits::dimWorld, MappingImpl > Mapping;
static const unsigned int dimension = Traits::dimension;
static const unsigned int dimWorld = Traits::dimWorld;
typedef typename Traits::FieldType FieldType;
typedef typename Traits::LocalCoordinate LocalCoordinate;
typedef typename Traits::GlobalCoordinate GlobalCoordinate;
typedef CachedStorage< dimension, GeometryTraits > Storage;
typedef CachedJacobianTransposed< dimension, GeometryTraits > JacobianTransposed;
typedef CachedJacobianInverseTransposed< dimension, GeometryTraits > JacobianInverseTransposed;
typedef GenericGeometry::ReferenceElement< Topology, FieldType > ReferenceElement;
// can we safely assume that this mapping is always affine?
static const bool alwaysAffine = Mapping::alwaysAffine;
typedef typename GeometryTraits::Caching Caching;
private:
typedef typename Traits::MatrixHelper MatrixHelper;
public:
template< class CoordVector >
explicit CachedMapping ( const CoordVector &coords )
: mapping_( coords )
{
if( alwaysAffine )
storage().affine = true;
else
computeJacobianTransposed( baryCenter() );
preCompute();
}
template< class CoordVector >
explicit CachedMapping ( const std::pair< const CoordVector &, bool > &coords )
: mapping_( coords.first )
{
storage().affine = coords.second;
preCompute();
}
bool affine () const { return (alwaysAffine || storage().affine); }
Dune::GeometryType type () const { return Dune::GeometryType( Topology() ); }
int numCorners () const { return ReferenceElement::numCorners; }
GlobalCoordinate corner ( int i ) const { return mapping().corner( i ); }
GlobalCoordinate center () const { return global( ReferenceElement::baryCenter() ); }
static bool checkInside ( const LocalCoordinate &x ) { return ReferenceElement::checkInside( x ); }
GlobalCoordinate global ( const LocalCoordinate &x ) const
{
GlobalCoordinate y;
if( jacobianTransposed() )
{
y = corner( 0 );
jacobianTransposed().umtv( x, y );
//MatrixHelper::template ATx< dimension, dimWorld >( jacobianTransposed_, x, y );
//y += corner( 0 );
}
else
mapping().global( x, y );
return y;
}
LocalCoordinate local ( const GlobalCoordinate &y ) const
{
LocalCoordinate x;
if( jacobianInverseTransposed() )
{
GlobalCoordinate z = y - corner( 0 );
jacobianInverseTransposed().mtv( z, x );
// MatrixHelper::template ATx< dimWorld, dimension >( jacobianInverseTransposed(), z, x );
}
else if( affine() )
{
const JacobianTransposed &JT = jacobianTransposed( baryCenter() );
GlobalCoordinate z = y - corner( 0 );
MatrixHelper::template xTRightInvA< dimension, dimWorld >( JT, z, x );
}
else
mapping().local( y, x );
return x;
}
FieldType integrationElement ( const LocalCoordinate &x ) const
{
const EvaluationType evaluateI = Caching::evaluateIntegrationElement;
const EvaluationType evaluateJ = Caching::evaluateJacobianInverseTransposed;
if( ((evaluateI == PreCompute) || (evaluateJ == PreCompute)) && alwaysAffine )
return storage().integrationElement;
else
return jacobianTransposed( x ).det();
}
FieldType volume () const
{
// do we need a quadrature of higher order, here?
const FieldType refVolume = ReferenceElement::volume();
return refVolume * integrationElement( baryCenter() );
}
const JacobianTransposed &jacobianTransposed ( const LocalCoordinate &x ) const
{
const EvaluationType evaluate = Caching::evaluateJacobianTransposed;
if( (evaluate == PreCompute) && alwaysAffine )
return jacobianTransposed();
if( !jacobianTransposed() )
computeJacobianTransposed( x );
return jacobianTransposed();
}
const JacobianInverseTransposed &
jacobianInverseTransposed ( const LocalCoordinate &x ) const
{
const EvaluationType evaluate = Caching::evaluateJacobianInverseTransposed;
if( (evaluate == PreCompute) && alwaysAffine )
return jacobianInverseTransposed();
if( !jacobianInverseTransposed() )
computeJacobianInverseTransposed( x );
return jacobianInverseTransposed();
}
const Mapping &mapping () const { return mapping_; }
private:
static const LocalCoordinate &baryCenter ()
{
return ReferenceElement::baryCenter();
}
Storage &storage () const
{
return jacobianInverseTransposed().storage();
}
const JacobianTransposed &jacobianTransposed () const
{
return jacobianInverseTransposed().jacobianTransposed();
}
const JacobianInverseTransposed &jacobianInverseTransposed () const
{
return jacobianInverseTransposed_;
}
void preCompute ()
{
assert( affine() == mapping().jacobianTransposed( baryCenter(), storage().jacobianTransposed ) );
if( !affine() )
return;
if( (Caching::evaluateJacobianTransposed == PreCompute) && !jacobianTransposed() )
computeJacobianTransposed( baryCenter() );
if( Caching::evaluateJacobianInverseTransposed == PreCompute )
computeJacobianInverseTransposed( baryCenter() );
else if( Caching::evaluateIntegrationElement == PreCompute )
jacobianTransposed().det();
}
void computeJacobianTransposed ( const LocalCoordinate &x ) const
{
storage().affine = mapping().jacobianTransposed( x, storage().jacobianTransposed );
storage().jacobianTransposedComputed = affine();
}
void computeJacobianInverseTransposed ( const LocalCoordinate &x ) const
{
storage().integrationElement
= MatrixHelper::template rightInvA< dimension, dimWorld >( jacobianTransposed( x ), storage().jacobianInverseTransposed );
storage().integrationElementComputed = affine();
storage().jacobianInverseTransposedComputed = affine();
}
private:
Mapping mapping_;
JacobianInverseTransposed jacobianInverseTransposed_;
};
} // namespace GenericGeometry
} // namespace Dune
#endif // #ifndef DUNE_GEOMETRY_GENERICGEOMETRY_CACHED_MAPPING_HH
|