4 #ifndef DUNE_GEOMETRY_AXISALIGNED_CUBE_GEOMETRY_HH
5 #define DUNE_GEOMETRY_AXISALIGNED_CUBE_GEOMETRY_HH
13 #include <dune/common/fvector.hh>
14 #include <dune/common/fmatrix.hh>
15 #include <dune/common/diagonalmatrix.hh>
16 #include <dune/common/unused.hh>
48 template <
class CoordType,
unsigned int dim,
unsigned int coorddim>
76 typedef typename conditional<dim==coorddim,
77 DiagonalMatrix<ctype,dim>,
86 typedef typename conditional<dim==coorddim,
87 DiagonalMatrix<ctype,dim>,
95 const Dune::FieldVector<ctype,coorddim> upper)
99 jacobianTransposed_(0),
100 jacobianInverseTransposed_(0)
103 axes_ = (1<<coorddim)-1;
114 const Dune::FieldVector<ctype,coorddim> upper,
115 const std::bitset<coorddim>& axes)
119 jacobianTransposed_(0),
120 jacobianInverseTransposed_(0)
122 assert(axes.count()==dim);
123 for (
size_t i=0; i<coorddim; i++)
125 upper_[i] = lower_[i];
134 jacobianTransposed_(0),
135 jacobianInverseTransposed_(0)
141 lower_ = other.lower_;
142 upper_ = other.upper_;
144 jacobianTransposed_ = other.jacobianTransposed_;
145 jacobianInverseTransposed_ = other.jacobianInverseTransposed_;
159 if (dim == coorddim) {
160 for (
size_t i=0; i<coorddim; i++)
161 result[i] = lower_[i] + local[i]*(upper_[i] - lower_[i]);
166 for (
size_t i=0; i<coorddim; i++)
167 result[i] = (axes_[i])
168 ? lower_[i] + local[lc++]*(upper_[i] - lower_[i])
178 if (dim == coorddim) {
179 for (
size_t i=0; i<dim; i++)
180 result[i] = (global[i] - lower_[i]) / (upper_[i] - lower_[i]);
181 }
else if (dim != 0) {
183 for (
size_t i=0; i<coorddim; i++)
185 result[lc++] = (global[i] - lower_[i]) / (upper_[i] - lower_[i]);
194 return jacobianTransposed_;
201 return jacobianInverseTransposed_;
220 for (
size_t i=0; i<coorddim; i++)
221 result[i] = 0.5 * (lower_[i] + upper_[i]);
236 if (dim == coorddim) {
237 for (
size_t i=0; i<coorddim; i++)
238 result[i] = (k & (1<<i)) ? upper_[i] : lower_[i];
242 unsigned int mask = 1;
244 for (
size_t i=0; i<coorddim; i++) {
246 result[i] = lower_[i];
248 result[i] = (k & mask) ? upper_[i] : lower_[i];
262 if (dim == coorddim) {
263 for (
size_t i=0; i<dim; i++)
264 vol *= upper_[i] - lower_[i];
266 }
else if (dim != 0) {
267 for (
size_t i=0; i<coorddim; i++)
269 vol *= upper_[i] - lower_[i];
284 for (
size_t i=0; i<dim; i++)
285 jacobianTransposed.diagonal()[i] = upper_[i] - lower_[i];
295 for (
size_t i=0; i<coorddim; i++)
297 jacobianTransposed[lc++][i] = upper_[i] - lower_[i];
303 for (
size_t i=0; i<dim; i++)
304 jacobianInverseTransposed.diagonal()[i] = 1.0 / (upper_[i] - lower_[i]);
314 for (
size_t i=0; i<coorddim; i++)
316 jacobianInverseTransposed[i][lc++] = 1.0 / (upper_[i] - lower_[i]);
319 Dune::FieldVector<ctype,coorddim> lower_;
321 Dune::FieldVector<ctype,coorddim> upper_;
323 std::bitset<coorddim> axes_;