From: leicht Date: Fri, 9 Feb 2007 06:59:30 +0000 (+0000) Subject: Introduce new flags face_flip and face_rotation for cells in 3D and line_orientation... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f2b96ac1977780e64f7db014641061ea06d61ebc;p=dealii-svn.git Introduce new flags face_flip and face_rotation for cells in 3D and line_orientation for quads in 3D. This is necessary in order to treat previously non-orientable meshes. This is a partial merge from branch_general_meshes. git-svn-id: https://svn.dealii.org/trunk@14444 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/grid/tria_accessor.h b/deal.II/deal.II/include/grid/tria_accessor.h index da1a4559ae..c40c304a32 100644 --- a/deal.II/deal.II/include/grid/tria_accessor.h +++ b/deal.II/deal.II/include/grid/tria_accessor.h @@ -788,13 +788,13 @@ class TriaObjectAccessor : public TriaAccessor * (@p false). Which is the * standard direction is * documented with the - * Triangulation class. In + * GeometryInfo class. In * 1d and 2d, this is always * @p true, but in 3d it may be * different, see the respective * discussion in the * documentation of the - * Triangulation class. + * GeometryInfo class. * * This function is really only * for internal use in the @@ -802,7 +802,124 @@ class TriaObjectAccessor : public TriaAccessor * know what this is all about. */ bool face_orientation (const unsigned int face) const; + + /** + * Return whether the face with index @p + * face is rotated by 180 degrees (@p true) + * or or not (@p false). In 1d and 2d, this + * is always @p false, but in 3d it may be + * different, see the respective discussion + * in the documentation of the + * GeometryInfo class. + * + * This function is really only + * for internal use in the + * library unless you absolutely + * know what this is all about. + */ + bool face_flip (const unsigned int face) const; + + /** + * Return whether the face with index @p + * face is rotated by 90 degrees (@p true) + * or or not (@p false). In 1d and 2d, this + * is always @p false, but in 3d it may be + * different, see the respective discussion + * in the documentation of the + * GeometryInfo class. + * + * This function is really only + * for internal use in the + * library unless you absolutely + * know what this is all about. + */ + bool face_rotation (const unsigned int face) const; + + /** + * Return whether the line with index @p + * line is oriented in standard + * direction. @p true indicates, that the + * line is oriented from vertex 0 to vertex + * 1, whereas it is the other way around + * otherwise. In 1d and 2d, this is always + * @p true, but in 3d it may be different, + * see the respective discussion in the + * documentation of the + * GeometryInfo classe. + * + * This function is really only + * for internal use in the + * library unless you absolutely + * know what this is all about. + */ + bool line_orientation (const unsigned int line) const; + /** + * Set whether the quad with + * index @p face has its normal + * pointing in the standard + * direction (@p true) or + * whether it is the opposite + * (@p false). Which is the + * standard direction is + * documented with the + * GeometryInfo class. + * + * This function is only for + * internal use in the + * library. Setting this flag to + * any other value than the one + * that the triangulation has + * already set is bound to bring + * you desaster. + * + * It is only possible to set the + * face_orientation of cells in 3d + * (i.e. celldim==3 && + * dim==3). + */ + void set_face_orientation (const unsigned int face, + const bool orientation) const; + + /** + * Set the flag indicating, what + * face_flip() will + * return. + * + * It is only possible to set the + * face_flip of cells in 3d + * (i.e. celldim==3 && + * dim==3). + */ + void set_face_flip (const unsigned int face, + const bool flip) const; + + /** + * Set the flag indicating, what + * face_rotation() will + * return. + * + * It is only possible to set the + * face_rotation of cells in 3d + * (i.e. celldim==3 && + * dim==3). + */ + void set_face_rotation (const unsigned int face, + const bool rotation) const; + + /** + * Set the flag indicating, what + * line_orientation() will + * return. + * + * It is only possible to set the + * line_orientation of faces in 3d + * (i.e. celldim==2 && + * dim==3). + */ + void set_line_orientation (const unsigned int line, + const bool orientation) const; + private: /** * Copy operator. This is normally @@ -1307,13 +1424,13 @@ class TriaObjectAccessor<1, dim> : public TriaAccessor * (@p false). Which is the * standard direction is * documented with the - * Triangulation class. In + * GeometryInfo class. In * 1d and 2d, this is always * @p true, but in 3d it may be * different, see the respective * discussion in the * documentation of the - * Triangulation class. + * GeometryInfo class. * * This function is really only * for internal use in the @@ -1322,6 +1439,70 @@ class TriaObjectAccessor<1, dim> : public TriaAccessor */ bool face_orientation (const unsigned int face) const; + /** + * Return whether the face with index @p + * face is rotated by 180 degrees (@p true) + * or or not (@p false). In 1d and 2d, this + * is always @p false, but in 3d it may be + * different, see the respective discussion + * in the documentation of the + * GeometryInfo class. + * + * This function is really only + * for internal use in the + * library unless you absolutely + * know what this is all about. + */ + bool face_flip (const unsigned int face) const; + + /** + * Return whether the face with index @p + * face is rotated by 90 degrees (@p true) + * or or not (@p false). In 1d and 2d, this + * is always @p false, but in 3d it may be + * different, see the respective discussion + * in the documentation of the + * GeometryInfo class. + * + * This function is really only + * for internal use in the + * library unless you absolutely + * know what this is all about. + */ + bool face_rotation (const unsigned int face) const; + + /** + * Return whether the line with index @p + * line is oriented in standard + * direction. @p true indicates, that the + * line is oriented from vertex 0 to vertex + * 1, whereas it is the other way around + * otherwise. In 1d and 2d, this is always + * @p true, but in 3d it may be different, + * see the respective discussion in the + * documentation of the + * GeometryInfo class. + * + * This function is really only + * for internal use in the + * library unless you absolutely + * know what this is all about. + */ + bool line_orientation (const unsigned int line) const; + + /** + * Set the flag indicating, what + * line_orientation() will + * return. + * + * It is only possible to set the + * line_orientation of faces in 3d + * (i.e. celldim==2 && + * dim==3). + */ + void set_line_orientation (const unsigned int line, + const bool orientation) const; + private: /** * Copy operator. This is normally @@ -1857,13 +2038,13 @@ class TriaObjectAccessor<2, dim> : public TriaAccessor * (@p false). Which is the * standard direction is * documented with the - * Triangulation class. In + * GeometryInfo class. In * 1d and 2d, this is always * @p true, but in 3d it may be * different, see the respective * discussion in the * documentation of the - * Triangulation class. + * GeometryInfo class. * * This function is really only * for internal use in the @@ -1872,6 +2053,70 @@ class TriaObjectAccessor<2, dim> : public TriaAccessor */ bool face_orientation (const unsigned int face) const; + /** + * Return whether the face with index @p + * face is rotated by 180 degrees (@p true) + * or or not (@p false). In 1d and 2d, this + * is always @p false, but in 3d it may be + * different, see the respective discussion + * in the documentation of the + * GeometryInfo class. + * + * This function is really only + * for internal use in the + * library unless you absolutely + * know what this is all about. + */ + bool face_flip (const unsigned int face) const; + + /** + * Return whether the face with index @p + * face is rotated by 90 degrees (@p true) + * or or not (@p false). In 1d and 2d, this + * is always @p false, but in 3d it may be + * different, see the respective discussion + * in the documentation of the + * GeometryInfo class. + * + * This function is really only + * for internal use in the + * library unless you absolutely + * know what this is all about. + */ + bool face_rotation (const unsigned int face) const; + + /** + * Return whether the line with index @p + * line is oriented in standard + * direction. @p true indicates, that the + * line is oriented from vertex 0 to vertex + * 1, whereas it is the other way around + * otherwise. In 1d and 2d, this is always + * @p true, but in 3d it may be different, + * see the respective discussion in the + * documentation of the + * GeometryInfo class. + * + * This function is really only + * for internal use in the + * library unless you absolutely + * know what this is all about. + */ + bool line_orientation (const unsigned int line) const; + + /** + * Set the flag indicating, what + * line_orientation() will + * return. + * + * It is only possible to set the + * line_orientation of faces in 3d + * (i.e. celldim==2 && + * dim==3). + */ + void set_line_orientation (const unsigned int line, + const bool orientation) const; + private: /** * Copy operator. This is @@ -2394,13 +2639,13 @@ class TriaObjectAccessor<3, dim> : public TriaAccessor * (@p false). Which is the * standard direction is * documented with the - * Triangulation class. In + * GeometryInfo class. In * 1d and 2d, this is always * @p true, but in 3d it may be * different, see the respective * discussion in the * documentation of the - * Triangulation class. + * GeometryInfo class. * * This function is really only * for internal use in the @@ -2408,7 +2653,57 @@ class TriaObjectAccessor<3, dim> : public TriaAccessor * know what this is all about. */ bool face_orientation (const unsigned int face) const; + /** + * Return whether the face with index @p + * face is rotated by 180 degrees (@p true) + * or or not (@p false). In 1d and 2d, this + * is always @p false, but in 3d it may be + * different, see the respective discussion + * in the documentation of the + * GeometryInfo class. + * + * This function is really only + * for internal use in the + * library unless you absolutely + * know what this is all about. + */ + bool face_flip (const unsigned int face) const; + /** + * Return whether the face with index @p + * face is rotated by 90 degrees (@p true) + * or or not (@p false). In 1d and 2d, this + * is always @p false, but in 3d it may be + * different, see the respective discussion + * in the documentation of the + * GeometryInfo class. + * + * This function is really only + * for internal use in the + * library unless you absolutely + * know what this is all about. + */ + bool face_rotation (const unsigned int face) const; + + /** + * Return whether the line with index @p + * line is oriented in standard + * direction. @p true indicates, that the + * line is oriented from vertex 0 to vertex + * 1, whereas it is the other way around + * otherwise. In 1d and 2d, this is always + * @p true, but in 3d it may be different, + * see the respective discussion in the + * documentation of the + * GeometryInfo class. + * + * This function is really only + * for internal use in the + * library unless you absolutely + * know what this is all about. + */ + bool line_orientation (const unsigned int line) const; + /** * Set whether the quad with * index @p face has its normal @@ -2418,7 +2713,7 @@ class TriaObjectAccessor<3, dim> : public TriaAccessor * (@p false). Which is the * standard direction is * documented with the - * Triangulation class. + * GeometryInfo class. * * This function is only for * internal use in the @@ -2430,6 +2725,32 @@ class TriaObjectAccessor<3, dim> : public TriaAccessor */ void set_face_orientation (const unsigned int face, const bool orientation) const; + + /** + * Set the flag indicating, what + * face_flip() will + * return. + * + * It is only possible to set the + * face_orientation of cells in 3d + * (i.e. celldim==3 && + * dim==3). + */ + void set_face_flip (const unsigned int face, + const bool flip) const; + + /** + * Set the flag indicating, what + * face_rotation() will + * return. + * + * It is only possible to set the + * face_orientation of cells in 3d + * (i.e. celldim==3 && + * dim==3). + */ + void set_face_rotation (const unsigned int face, + const bool rotation) const; private: /** @@ -2747,15 +3068,13 @@ class CellAccessor : public TriaObjectAccessor * for the index of the * child. * - * However, the function - * is more complicated in 3d, - * since there faces may have - * more than one orientation, and - * we have to use - * @p face_orientation for both - * this and the neighbor cell to - * figure out which cell we want - * to have. + * However, the function is more + * complicated in 3d, since there faces may + * have more than one orientation, and we + * have to use @p face_orientation, @p + * face_flip and @p face_rotation for both + * this and the neighbor cell to figure out + * which cell we want to have. * * This can lead to surprising * results: if we are sitting on @@ -2787,6 +3106,12 @@ class CellAccessor : public TriaObjectAccessor * orderings are only identical * if the face orientation is @p * true, and reversed otherwise. + * + * Similarly, effects of + * face_flip()==true and + * face_rotation()==true(), both + * of which indicate a non-standard face + * have to be considered. * * Fortunately, this is only very rarely * of concern, since one is usually only diff --git a/deal.II/deal.II/include/grid/tria_accessor.templates.h b/deal.II/deal.II/include/grid/tria_accessor.templates.h index a34928c5eb..ca3258bd5b 100644 --- a/deal.II/deal.II/include/grid/tria_accessor.templates.h +++ b/deal.II/deal.II/include/grid/tria_accessor.templates.h @@ -302,6 +302,33 @@ TriaObjectAccessor<1,dim>::face_orientation (const unsigned int) const } +template +inline +bool +TriaObjectAccessor<1,dim>::face_flip (const unsigned int) const +{ + return false; +} + + +template +inline +bool +TriaObjectAccessor<1,dim>::face_rotation (const unsigned int) const +{ + return false; +} + + +template +inline +bool +TriaObjectAccessor<1,dim>::line_orientation (const unsigned int) const +{ + return true; +} + + template inline @@ -594,6 +621,53 @@ TriaObjectAccessor<2,dim>::face_orientation (const unsigned int) const } +template +inline +bool +TriaObjectAccessor<2,dim>::face_flip (const unsigned int) const +{ + return false; +} + + +template +inline +bool +TriaObjectAccessor<2,dim>::face_rotation (const unsigned int) const +{ + return false; +} + + +template <> +inline +bool +TriaObjectAccessor<2,3>::line_orientation (const unsigned int line) const +{ + // we cannot usae the quads() function here, + // since it returns a reference to + // TriaObjects, but we need a + // (reference to) TriaObjectsQuad3D + const int dim=3; + Assert (used(), TriaAccessor::ExcCellNotUsed()); + Assert (line::lines_per_cell, + ExcIndexRange (line, 0, GeometryInfo<2>::lines_per_cell)); + Assert (this->present_index * GeometryInfo<2>::lines_per_cell + line + < this->tria->faces->quads.line_orientations.size(), + ExcInternalError()); + return this->tria->faces->quads. + line_orientations[this->present_index * GeometryInfo<2>::lines_per_cell + line]; +} + + +template +inline +bool +TriaObjectAccessor<2,dim>::line_orientation (const unsigned int) const +{ + return true; +} + template inline @@ -785,35 +859,43 @@ TriaObjectAccessor<3,dim>::line_index (const unsigned int i) const Assert (i<12, ExcIndexRange(i,0,12)); // get the line index by asking the - // quads. make sure we handle - // reverted faces correctly + // quads. first assume standard orientation // - // so set up a table that for each + // set up a table that for each // line describes a) from which // quad to take it, b) which line // therein it is if the face is - // oriented correctly, and c) if in - // the opposite direction - static const unsigned int lookup_table[12][3] = - { { 4, 0, 2 }, // take first four lines from bottom face - { 4, 1, 3 }, - { 4, 2, 0 }, - { 4, 3, 1 }, - - { 5, 0, 2 }, // second four lines from top face - { 5, 1, 3 }, - { 5, 2, 0 }, - { 5, 3, 1 }, - - { 0, 0, 2 }, // the rest randomly - { 1, 0, 2 }, - { 0, 1, 3 }, - { 1, 1, 3 }}; + // oriented correctly + static const unsigned int lookup_table[12][2] = + { { 4, 0 }, // take first four lines from bottom face + { 4, 1 }, + { 4, 2 }, + { 4, 3 }, + + { 5, 0 }, // second four lines from top face + { 5, 1 }, + { 5, 2 }, + { 5, 3 }, + + { 0, 0 }, // the rest randomly + { 1, 0 }, + { 0, 1 }, + { 1, 1 }}; + + // respect non-standard faces by calling the + // reordering function from GeometryInfo + + const unsigned int quad_index=lookup_table[i][0]; + const unsigned int std_line_index=lookup_table[i][1]; + + const unsigned int line_index=GeometryInfo::standard_to_real_face_line( + std_line_index, + face_orientation(quad_index), + face_flip(quad_index), + face_rotation(quad_index)); - return (this->quad(lookup_table[i][0]) - ->line_index(face_orientation(lookup_table[i][0]) ? - lookup_table[i][1] : - lookup_table[i][2])); + return (this->quad(quad_index) + ->line_index(line_index)); } @@ -916,15 +998,147 @@ face_orientation (const unsigned int face) const Assert (used(), TriaAccessor::ExcCellNotUsed()); Assert (face::faces_per_cell, ExcIndexRange (face, 0, GeometryInfo<3>::faces_per_cell)); - Assert (this->present_index * GeometryInfo<3>::faces_per_cell + face - < this->tria->levels[this->present_level] - ->hexes.face_orientations.size(), - ExcInternalError()); + Assert (this->present_index * GeometryInfo<3>::faces_per_cell + face + < this->tria->levels[this->present_level] + ->hexes.face_orientations.size(), + ExcInternalError()); + + return (this->tria->levels[this->present_level] + ->hexes.face_orientations[this->present_index * + GeometryInfo<3>::faces_per_cell + + face]); +} + + +template <> +inline +bool +TriaObjectAccessor<3, 3>:: +face_flip (const unsigned int face) const +{ + const int dim=3; + Assert (used(), TriaAccessor::ExcCellNotUsed()); + Assert (face::faces_per_cell, + ExcIndexRange (face, 0, GeometryInfo<3>::faces_per_cell)); + Assert (this->present_index * GeometryInfo<3>::faces_per_cell + face + < this->tria->levels[this->present_level] + ->hexes.face_flips.size(), + ExcInternalError()); + + return (this->tria->levels[this->present_level] + ->hexes.face_flips[this->present_index * + GeometryInfo<3>::faces_per_cell + + face]); +} + + +template <> +inline +bool +TriaObjectAccessor<3, 3>:: +face_rotation (const unsigned int face) const +{ + const int dim=3; + Assert (used(), TriaAccessor::ExcCellNotUsed()); + Assert (face::faces_per_cell, + ExcIndexRange (face, 0, GeometryInfo<3>::faces_per_cell)); + Assert (this->present_index * GeometryInfo<3>::faces_per_cell + face + < this->tria->levels[this->present_level] + ->hexes.face_rotations.size(), + ExcInternalError()); - return (this->tria->levels[this->present_level] - ->hexes.face_orientations[this->present_index * - GeometryInfo<3>::faces_per_cell - + face]); + return (this->tria->levels[this->present_level] + ->hexes.face_rotations[this->present_index * + GeometryInfo<3>::faces_per_cell + + face]); +} + + +template <> +inline +bool +TriaObjectAccessor<3, 3>:: +line_orientation (const unsigned int line) const +{ + const int dim=3; + Assert (used(), TriaAccessor::ExcCellNotUsed()); + Assert (line::lines_per_cell, + ExcIndexRange (line, 0, GeometryInfo<3>::lines_per_cell)); + + // get the line index by asking the + // quads. first assume standard orientation + // + // set up a table that for each + // line describes a) from which + // quad to take it, b) which line + // therein it is if the face is + // oriented correctly + static const unsigned int lookup_table[12][2] = + { { 4, 0 }, // take first four lines from bottom face + { 4, 1 }, + { 4, 2 }, + { 4, 3 }, + + { 5, 0 }, // second four lines from top face + { 5, 1 }, + { 5, 2 }, + { 5, 3 }, + + { 0, 0 }, // the rest randomly + { 1, 0 }, + { 0, 1 }, + { 1, 1 }}; + + const unsigned int quad_index=lookup_table[line][0]; + const unsigned int std_line_index=lookup_table[line][1]; + + const unsigned int line_index=GeometryInfo::standard_to_real_face_line( + std_line_index, + face_orientation(quad_index), + face_flip(quad_index), + face_rotation(quad_index)); + + // now we got to the correct line and ask + // the quad for its line_orientation. however, if + // the face is rotated, it might be possible, + // that a standard orientation of the line + // with respect to the face corrsponds to a + // non-standard orientation for the line with + // respect to the cell. + // + // set up a table indicating if the two + // standard orientations coincide + // + // first index: two pairs of lines 0(lines + // 0/1) and 1(lines 2/3) + // + // second index: face_orientation; 0: + // opposite normal, 1: standard + // + // third index: face_flip; 0: standard, 1: + // face rotated by 180 degrees + // + // forth index: face_rotation: 0: standard, + // 1: face rotated by 90 degrees + + static const bool bool_table[2][2][2][2] = + { { { { true, false }, // lines 0/1, face_orientation=false, face_flip=false, face_rotation=false and true + { false, true }}, // lines 0/1, face_orientation=false, face_flip=true, face_rotation=false and true + { { true, true }, // lines 0/1, face_orientation=true, face_flip=false, face_rotation=false and true + { false, false }}},// linea 0/1, face_orientation=true, face_flip=true, face_rotation=false and true + + { { { true, true }, // lines 2/3 ... + { false, false }}, + { { true, false }, + { false, true }}}}; + + + return (this->quad(quad_index) + ->line_orientation(line_index) + == bool_table[std_line_index/2] + [face_orientation(quad_index)] + [face_flip(quad_index)] + [face_rotation(quad_index)]); } diff --git a/deal.II/deal.II/include/grid/tria_faces.h b/deal.II/deal.II/include/grid/tria_faces.h index a095481006..8981bdbb09 100644 --- a/deal.II/deal.II/include/grid/tria_faces.h +++ b/deal.II/deal.II/include/grid/tria_faces.h @@ -110,7 +110,7 @@ namespace internal * The TriaObject containing the data of quads. */ - TriaObjects quads; + TriaObjectsQuad3D quads; /** * The TriaObject containing the data of lines. diff --git a/deal.II/deal.II/include/grid/tria_objects.h b/deal.II/deal.II/include/grid/tria_objects.h index f3065e9412..ce6ca946b5 100644 --- a/deal.II/deal.II/include/grid/tria_objects.h +++ b/deal.II/deal.II/include/grid/tria_objects.h @@ -23,6 +23,14 @@ DEAL_II_NAMESPACE_OPEN +// TODO: See if we can unify the class hierarchy a bit using partial +// specialization of classes here, e.g. declare a general class +// TriaObjects. +// To consider for this: in that case we would have to duplicate quite a few +// things, e.g. TriaObjectsHex is derived from TriaObjects and +// declares mainly additional data. This would have to be changed in case of a +// partial specialization. + namespace internal { namespace Triangulation @@ -236,6 +244,14 @@ namespace internal * has. */ std::vector face_orientations; + /** + * flip = rotation by 180 degrees + */ + std::vector face_flips; + /** + * rotation by 90 degrees + */ + std::vector face_rotations; /** * Assert that enough space is @@ -270,6 +286,63 @@ namespace internal */ unsigned int memory_consumption () const; }; + + +/** + * For quadrilaterals in 3D the data of TriaObjects needs to be extended, as we + * can obtain faces (quads) with lines in non-standard-orientation, therefore we + * declare a class TriaObjectsQuad3D, which additionaly contains a bool-vector + * of the line-orientations. + * @ingroup grid + */ + + class TriaObjectsQuad3D: public TriaObjects + { + public: + + /** + * In effect, this field has + * 4*n_quads elements, + * being the number of quads + * times the four lines each + * has. + */ + std::vector line_orientations; + + /** + * Assert that enough space is + * allocated to accomodate + * new_objs new objects. + * This function does not only call + * vector::reserve(), but + * does really append the needed + * elements. + */ + void reserve_space (const unsigned int new_objs); + + /** + * Clear all the data contained in this object. + */ + void clear(); + + /** + * Check the memory consistency of the + * different containers. Should only be + * called with the prepro flag @p DEBUG + * set. The function should be called from + * the functions of the higher + * TriaLevel classes. + */ + void monitor_memory (const unsigned int true_dimension) const; + + /** + * Determine an estimate for the + * memory consumption (in bytes) + * of this object. + */ + unsigned int memory_consumption () const; + }; + } } diff --git a/deal.II/deal.II/source/grid/tria_accessor.cc b/deal.II/deal.II/source/grid/tria_accessor.cc index 92605d8832..a554bb5ce8 100644 --- a/deal.II/deal.II/source/grid/tria_accessor.cc +++ b/deal.II/deal.II/source/grid/tria_accessor.cc @@ -282,7 +282,18 @@ int TriaObjectAccessor<2, dim>::vertex_index (const unsigned int corner) const { Assert (corner<4, ExcIndexRange(corner,0,4)); - return line(corner%2)->vertex_index(corner/2); + // table used to switch the vertices, if the + // line orientation is wrong, + // + // first index: line orientation 0: false or + // 1: true=standard + // + // second index: vertex index to be switched + // (or not) + + static const unsigned int switch_table[2][2]={{1,0},{0,1}}; + + return line(corner%2)->vertex_index(switch_table[line_orientation(corner%2)][corner/2]); } @@ -677,6 +688,41 @@ unsigned int TriaObjectAccessor<2, dim>::number_of_children () const } +template +void +TriaObjectAccessor<2, dim>:: +set_line_orientation (const unsigned int , + const bool ) const +{ + Assert(false, + ExcMessage("The line orientation of Quads can only be set in 3D. In 2D the lines are alway in standard orientation."));} + + + +#if deal_II_dimension == 3 + +template <> +void +TriaObjectAccessor<2, 3>:: +set_line_orientation (const unsigned int line, + const bool orientation) const +{ + const int dim=3; + Assert (used(), TriaAccessor::ExcCellNotUsed()); + Assert (line::lines_per_face, + ExcIndexRange (line, 0, GeometryInfo<3>::lines_per_face)); + Assert (this->present_index * GeometryInfo::lines_per_face + line + < this->tria->faces->quads.line_orientations.size(), + ExcInternalError()); + + this->tria->faces + ->quads.line_orientations[this->present_index * + GeometryInfo::lines_per_face + + + line] = orientation; +} + +#endif /*------------------------ Functions: TriaObjectAccessor ---------------------------*/ @@ -704,30 +750,21 @@ int TriaObjectAccessor<3, dim>::vertex_index (const unsigned int corner) const { Assert (corner<8, ExcIndexRange(corner,0,8)); - // get the corner indices by asking - // either the bottom or the top - // face for its vertices. make sure - // we take into account that the - // face might have non-standard - // orientation; if this is not the - // case, transpose vertices 1 and 2 - // of this face - static const unsigned int vertex_translation[4] = - { 0, 2, 1, 3 }; - if (corner<4) - { - if (face_orientation(4) == true) - return quad(4)->vertex_index(corner); - else - return quad(4)->vertex_index(vertex_translation[corner]); - } - else - { - if (face_orientation(5) == true) - return quad(5)->vertex_index(corner-4); - else - return quad(5)->vertex_index(vertex_translation[corner-4]); - } + // get the corner indices by asking either + // the bottom or the top face for its + // vertices. handle non-standard faces by + // calling the vertex reordering function + // from GeometryInfo + + // bottom face (4) for first four vertices, + // top face (5) for the rest + const unsigned int face_index=4+corner/4; + + return quad(face_index)->vertex_index(GeometryInfo + ::standard_to_real_face_vertex(corner%4, + face_orientation(face_index), + face_flip(face_index), + face_rotation(face_index))); } @@ -1649,6 +1686,52 @@ set_face_orientation (const unsigned int face, face] = orientation; } + +template <> +void +TriaObjectAccessor<3, 3>:: +set_face_flip (const unsigned int face, + const bool flip) const +{ + const int dim=3; + Assert (used(), TriaAccessor::ExcCellNotUsed()); + Assert (face::faces_per_cell, + ExcIndexRange (face, 0, GeometryInfo<3>::faces_per_cell)); + Assert (this->present_index * GeometryInfo<3>::faces_per_cell + face + < this->tria->levels[this->present_level] + ->hexes.face_flips.size(), + ExcInternalError()); + + this->tria->levels[this->present_level] + ->hexes.face_flips[this->present_index * + GeometryInfo<3>::faces_per_cell + + + face] = flip; +} + + +template <> +void +TriaObjectAccessor<3, 3>:: +set_face_rotation (const unsigned int face, + const bool rotation) const +{ + const int dim=3; + Assert (used(), TriaAccessor::ExcCellNotUsed()); + Assert (face::faces_per_cell, + ExcIndexRange (face, 0, GeometryInfo<3>::faces_per_cell)); + Assert (this->present_index * GeometryInfo<3>::faces_per_cell + face + < this->tria->levels[this->present_level] + ->hexes.face_rotations.size(), + ExcInternalError()); + + this->tria->levels[this->present_level] + ->hexes.face_rotations[this->present_index * + GeometryInfo<3>::faces_per_cell + + + face] = rotation; +} + #endif // Remark: The following explicit instantiations needed to be moved to // this place here to work around a problem with gcc3.3 on Apple MacOSX. @@ -2183,62 +2266,51 @@ neighbor_child_on_subface (const unsigned int face, ExcMessage ("The neighbor must be on the same level as this cell!")); Assert (this->neighbor(face)->has_children() == true, ExcMessage ("The neighbor must have children!")); - - static const unsigned int subface_translation[4] - = { 0, 2, 1, 3 }; - // see whether face and - // the neighbor's - // counterface share the - // same indexing of - // children. if not so, - // translate child - // indices + // this function returns they neighbor's + // child on a given face and + // subface. However, irrespective of the + // actual face orientation, flip and rotation + // of the cells face it returns that child, + // which is neighboring the subface which + // would be subface number @p subface, if the + // face was in standard orientation. thus, + // our first task is to extract the + // information, which real_subface_no + // corresponds to that for the actual face, + // which might have non-standard orientation, + // flip and rotation. + + // to this end: the subfaces are numbered + // according to the vertices of the + // face. thus we can use the + // standard_to_real_face_vertex function in + // GeometryInfo to translate the given + // subface to the real_subface_no + unsigned int real_subface_no=GeometryInfo<3>:: + standard_to_real_face_vertex(subface, + this->face_orientation(face), + this->face_flip(face), + this->face_rotation(face)); + // get the neighbor's number for the given + // face and the neighbor const unsigned int neighbor_neighbor = this->neighbor_of_neighbor (face); - const bool face_orientations_match - = (this->neighbor(face)->face_orientation(neighbor_neighbor) == - this->face_orientation(face)); + const TriaIterator<3,CellAccessor<3> > neighbor=this->neighbor(face); + + // now use the info provided by GeometryInfo + // to extract the neighbors child number const unsigned int neighbor_child_index - = GeometryInfo<3>::child_cell_on_face(neighbor_neighbor, subface, - face_orientations_match); + = GeometryInfo<3>::child_cell_on_face(neighbor_neighbor, real_subface_no, + neighbor->face_orientation(neighbor_neighbor), + neighbor->face_flip(neighbor_neighbor), + neighbor->face_rotation(neighbor_neighbor)); const TriaIterator<3,CellAccessor<3> > neighbor_child= - this->neighbor(face)->child(neighbor_child_index); + neighbor->child(neighbor_child_index); - // if the face on the side of the present - // cell is in the correct order, then make - // sure that the neighbor child cell we + // make sure that the neighbor child cell we // have found shares the desired subface. - // - // otherwise: if the face is turned - // inside out when viewed from the - // present cell, then the subface we are - // interested in is not - // this->face(face)->child(subface), but - // instead - // this->face(face)->child(subface'), - // where subface' is the translated - // subface number. this is so because we - // only store the face only once, so its - // children are defined in terms of its - // own (circular) orientation, not in - // terms of face_orientation as viewed - // from one of the adjacent cells. in - // that case, cell->face(f)->subface(sf) - // may yield unexpected results; in that - // case, a caller may need to adjust - // according to the face_orientation - // flag, though one in general only wants - // to loop over all subfaces, and not - // pick a particular one - Assert(((this->face_orientation(face) == true) - && - (this->face(face)->child(subface) == - neighbor_child->face(neighbor_neighbor))) - || - ((this->face_orientation(face) == false) - && - (this->face(face)->child(subface_translation[subface]) == - neighbor_child->face(neighbor_neighbor))), + Assert((this->face(face)->child(real_subface_no) == + neighbor_child->face(neighbor_neighbor)), ExcInternalError()); return neighbor_child; diff --git a/deal.II/deal.II/source/grid/tria_objects.cc b/deal.II/deal.II/source/grid/tria_objects.cc index bed78dfc3a..07b1f2c1b7 100644 --- a/deal.II/deal.II/source/grid/tria_objects.cc +++ b/deal.II/deal.II/source/grid/tria_objects.cc @@ -162,10 +162,41 @@ namespace internal new_size * GeometryInfo<3>::faces_per_cell - face_orientations.size(), true); + face_flips.reserve (new_size * GeometryInfo<3>::faces_per_cell); + face_flips.insert (face_flips.end(), + new_size * GeometryInfo<3>::faces_per_cell + - face_flips.size(), + false); + face_rotations.reserve (new_size * GeometryInfo<3>::faces_per_cell); + face_rotations.insert (face_rotations.end(), + new_size * GeometryInfo<3>::faces_per_cell + - face_rotations.size(), + false); }; } + void + TriaObjectsQuad3D::reserve_space (const unsigned int new_quads) + { + const unsigned int new_size = new_quads + + std::count_if (used.begin(), + used.end(), + std::bind2nd (std::equal_to(), true)); + + // see above... + if (new_size>cells.size()) + { + TriaObjects::reserve_space(new_quads); + line_orientations.reserve (new_size * GeometryInfo<2>::lines_per_cell); + line_orientations.insert (line_orientations.end(), + new_size * GeometryInfo<2>::lines_per_cell + - line_orientations.size(), + true); + }; + } + + template<> void TriaObjects::monitor_memory (const unsigned int) const @@ -286,6 +317,33 @@ namespace internal == face_orientations.size(), ExcMemoryInexact (cells.size() * GeometryInfo<3>::faces_per_cell, face_orientations.size())); + Assert (cells.size() * GeometryInfo<3>::faces_per_cell + == face_flips.size(), + ExcMemoryInexact (cells.size() * GeometryInfo<3>::faces_per_cell, + face_flips.size())); + Assert (cells.size() * GeometryInfo<3>::faces_per_cell + == face_rotations.size(), + ExcMemoryInexact (cells.size() * GeometryInfo<3>::faces_per_cell, + face_rotations.size())); + } + + + void + TriaObjectsQuad3D::monitor_memory (const unsigned int) const + { + // check that we have not allocated + // too much memory. note that bool + // vectors allocate their memory in + // chunks of whole integers, so + // they may over-allocate by up to + // as many elements as an integer + // has bits + Assert (cells.size() * GeometryInfo<2>::lines_per_cell + == line_orientations.size(), + ExcMemoryInexact (cells.size() * GeometryInfo<2>::lines_per_cell, + line_orientations.size())); + TriaObjects::monitor_memory (3); + } @@ -307,6 +365,16 @@ namespace internal { TriaObjects::clear(); face_orientations.clear(); + face_flips.clear(); + face_rotations.clear(); + } + + + void + TriaObjectsQuad3D::clear() + { + TriaObjects::clear(); + line_orientations.clear(); } @@ -327,7 +395,17 @@ namespace internal TriaObjectsHex::memory_consumption () const { return (MemoryConsumption::memory_consumption (face_orientations) + - this->TriaObjects::memory_consumption() ); + MemoryConsumption::memory_consumption (face_flips) + + MemoryConsumption::memory_consumption (face_rotations) + + TriaObjects::memory_consumption() ); + } + + + unsigned int + TriaObjectsQuad3D::memory_consumption () const + { + return (MemoryConsumption::memory_consumption (line_orientations) + + this->TriaObjects::memory_consumption() ); } // explicit instantiations