From 38f17e38f1b591a73e995ffa003a7dc397e2c680 Mon Sep 17 00:00:00 2001 From: Ralf Hartmann Date: Wed, 2 Nov 2005 07:45:31 +0000 Subject: [PATCH] Merge changes of branch_renumbering to main branch. git-svn-id: https://svn.dealii.org/trunk@11683 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/geometry_info.h | 201 +-- deal.II/base/source/data_out_base.cc | 3 +- deal.II/base/source/geometry_info.cc | 204 ++- deal.II/base/source/quadrature.cc | 112 +- deal.II/deal.II/include/fe/fe.h | 10 +- deal.II/deal.II/include/fe/fe_poly.h | 4 +- deal.II/deal.II/include/fe/fe_q.h | 122 +- .../deal.II/include/fe/fe_q_hierarchical.h | 121 +- deal.II/deal.II/include/fe/mapping.h | 35 - .../deal.II/include/fe/mapping_cartesian.h | 2 +- deal.II/deal.II/include/fe/mapping_q.h | 2 +- deal.II/deal.II/include/fe/mapping_q1.h | 16 +- .../deal.II/include/grid/persistent_tria.h | 12 +- deal.II/deal.II/include/grid/tria.h | 362 +++--- .../include/grid/tria_accessor.templates.h | 52 +- deal.II/deal.II/source/dofs/dof_tools.cc | 20 +- deal.II/deal.II/source/fe/fe_dgq.cc | 26 +- deal.II/deal.II/source/fe/fe_nedelec.cc | 149 +-- deal.II/deal.II/source/fe/fe_nedelec_2d.cc | 32 +- deal.II/deal.II/source/fe/fe_nedelec_3d.cc | 205 ++- deal.II/deal.II/source/fe/fe_q.cc | 32 +- .../deal.II/source/fe/fe_q_hierarchical.cc | 187 +-- .../source/fe/fe_tools.all_dimensions.cc | 99 +- deal.II/deal.II/source/fe/mapping.cc | 31 - .../deal.II/source/fe/mapping_cartesian.cc | 46 +- deal.II/deal.II/source/fe/mapping_q.cc | 181 +-- deal.II/deal.II/source/fe/mapping_q1.cc | 84 +- .../deal.II/source/fe/mapping_q1_eulerian.cc | 4 +- deal.II/deal.II/source/grid/grid_generator.cc | 191 ++- deal.II/deal.II/source/grid/grid_in.cc | 14 +- deal.II/deal.II/source/grid/grid_out.cc | 42 +- .../deal.II/source/grid/grid_reordering.cc | 42 +- deal.II/deal.II/source/grid/grid_tools.cc | 137 +- .../deal.II/source/grid/persistent_tria.cc | 11 + deal.II/deal.II/source/grid/tria.cc | 1123 ++++++++++------- deal.II/deal.II/source/grid/tria_accessor.cc | 139 +- deal.II/deal.II/source/grid/tria_boundary.cc | 12 +- .../deal.II/source/grid/tria_boundary_lib.cc | 12 +- deal.II/deal.II/source/numerics/data_out.cc | 44 +- .../deal.II/source/numerics/data_out_stack.cc | 20 +- tests/bits/fe_q_3d_01.cc | 42 +- tests/bits/mesh_3d_13.cc | 2 +- tests/bits/mesh_3d_14.cc | 2 +- tests/bits/mesh_3d_15.cc | 2 +- tests/bits/mesh_3d_2.cc | 2 +- tests/bits/mesh_3d_3.cc | 8 +- tests/bits/q_points.cc | 2 +- tests/deal.II/constraints.cc | 12 +- tests/deal.II/dof_test.cc | 2 +- tests/deal.II/gradients.cc | 14 +- tests/deal.II/grid_test.cc | 10 +- tests/deal.II/measure_et_al.cc | 2 +- tests/deal.II/second_derivatives.cc | 6 +- tests/deal.II/subcelldata.cc | 14 +- tests/fe/mapping.cc | 22 +- tests/fe/mapping_q1_eulerian.cc | 2 +- tests/fe/numbering.cc | 129 +- 57 files changed, 2350 insertions(+), 2064 deletions(-) diff --git a/deal.II/base/include/base/geometry_info.h b/deal.II/base/include/base/geometry_info.h index 606fa6c2b5..f3cb9349e3 100644 --- a/deal.II/base/include/base/geometry_info.h +++ b/deal.II/base/include/base/geometry_info.h @@ -108,16 +108,92 @@ struct GeometryInfo<0> +//! Topological description of cells. +/** + * Dimension independent base class for the GeometryInfo + * classes with dim=1,2,3,4. Includes all data and methods + * which can be defined in an dimension indendent way. + * + * @author Ralf Hartmann, 2005 + */ +class GeometryInfoBase +{ + private: + /** + * The maximal dimension for + * which data in this class is + * implemented. + */ + static const unsigned int max_dim = 4; + + /** + * Number of faces of a cell for + * dim=max_dim. + */ + static const unsigned int faces_per_cell_max_dim = 2*max_dim; + + public: + + /** + * For each face of the reference + * cell, this field stores the + * coordinate direction in which + * its normal vector points. In + * dim dimension these + * are the 2*dim first + * entries of + * {0,0,1,1,2,2,3,3}. + * + * Remark that this is only the + * coordinate number. The actual + * direction of the normal vector + * is obtained by multiplying the + * unit vector in this direction + * with #unit_normal_orientation. + */ + static const unsigned int unit_normal_direction[faces_per_cell_max_dim]; + + /** + * Orientation of the unit normal + * vector of a face of the + * reference cell. In + * dim dimension these + * are the 2*dim first + * entries of + * {-1,1,-1,1,-1,1,-1,1}. + * + * Each value is either + * 1 or -1, + * corresponding to a normal + * vector pointing in the + * positive or negative + * coordinate direction, + * respectively. + */ + static const int unit_normal_orientation[faces_per_cell_max_dim]; + + /** + * List of numbers which denotes + * which face is opposite to a + * given face. Its entries are + * { 1, 0, 3, 2, 5, 4, 7, 6}. + */ + static const unsigned int opposite_face[faces_per_cell_max_dim]; +}; + + + /** * This template specifies the interface to all topological structure * of the mesh cells. + * @ref Instantiations: few (dim = 1,2,3,4 and a specialization for dim=0) * * @ingroup grid * @ref Instantiations few (dim = 1,2,3,4 and a specialization for dim=0) * @author Wolfgang Bangerth, 1998, Ralf Hartmann, 2005 */ template -struct GeometryInfo +struct GeometryInfo: public GeometryInfoBase { /** @@ -198,14 +274,32 @@ struct GeometryInfo GeometryInfo::quads_per_cell); /** - * List of numbers which denotes which - * face is opposite to a given face. In - * 1d, this list is {1,0}, in 2d - * {2, 3, 0, 1}, in 3d {1, - * 0, 4, 5, 2, 3}. + * Rearrange vertices for UCD + * output. For a cell being + * written in UCD format, each + * entry in this field contains + * the number of a vertex in + * deal.II that corresponds + * to the UCD numbering at this + * location. + * + * Typical example: write a cell + * and arrange the vertices, such + * that UCD understands them. + * + * \begin{verbatim} + * for (i=0; i< n_vertices; ++i) + * out << cell->vertex(ucd_to_deal[i]); + * \end{verbatim} + * + * As the vertex numbering in + * deal.II versions <= 5.1 + * happened to coincide with the + * UCD numbering, this field can + * also be used like a + * old_to_lexicographic mapping. */ - static const unsigned int opposite_face[faces_per_cell]; - + static const unsigned int ucd_to_deal[vertices_per_cell]; /** * Rearrange vertices for OpenDX @@ -228,36 +322,6 @@ struct GeometryInfo */ static const unsigned int dx_to_deal[vertices_per_cell]; - /** - * For each face of the reference - * cell, this field stores the - * coordinate direction in which - * its normal vector points. - * - * Remark that this is only the - * coordinate number. The acual - * direction of the normal vector - * is obtained by multiplying the - * unit vector in this direction - * with #unit_normal_orientation. - */ - static const unsigned int unit_normal_direction[faces_per_cell]; - - /** - * Orientation of the unit normal - * vector of a face of the - * reference cell. - * - * Each value is either - * 1 or -1, - * corresponding to a normal - * vector pointing in the - * positive or negative - * coordinate direction, - * respectively. - */ - static const int unit_normal_orientation[faces_per_cell]; - /** * This field stores which child * cells are adjacent to a @@ -267,13 +331,13 @@ struct GeometryInfo * For example, in 2D the layout of * a cell is as follows: * @verbatim - * . 2 - * . 3-->--2 + * . 3 + * . 2-->--3 * . | | - * . 3 ^ ^ 1 + * . 0 ^ ^ 1 * . | | * . 0-->--1 - * . 0 + * . 2 * @endverbatim * Vertices and faces are indicated * with their numbers, faces also with @@ -283,16 +347,17 @@ struct GeometryInfo * like this: * @verbatim * *--*--* - * | 3|2 | + * | 2|3 | * *--*--* * | 0|1 | * *--*--* * @endverbatim * - * Thus, the child cells on face zero - * are (ordered in the direction of the - * face) 0 and 1, on face 2 they are - * 3 and 2, etc. + * Thus, the child cells on face + * 0 are (ordered in the + * direction of the face) 0 and + * 2, on face 3 they are 2 and 3, + * etc. * * For three spatial dimensions, * the exact order of the @@ -319,7 +384,7 @@ struct GeometryInfo * cell vertex number of the * vertexth vertex of * line line, e.g. - * GeometryInfo<2>::line_to_cell_vertices(2,0)=3. + * GeometryInfo<2>::line_to_cell_vertices(3,0)=2. * * The order of the lines, as * well as their direction (which @@ -345,7 +410,7 @@ struct GeometryInfo * cell vertex number of the * vertexth vertex of * face face, e.g. - * GeometryInfo<2>::face_to_cell_vertices(2,0)=3. + * GeometryInfo<2>::face_to_cell_vertices(3,0)=2. * * Through the * face_orientation @@ -378,7 +443,7 @@ struct GeometryInfo * cell line number of the * lineth line of face * face, e.g. - * GeometryInfo<3>::face_to_cell_lines(3,1)=5. + * GeometryInfo<3>::face_to_cell_lines(5,0)=4. * * Through the * face_orientation @@ -388,11 +453,8 @@ struct GeometryInfo * orientation. * face_orientation * defaults to true - * (standard orientation). - * - * This function is useful and - * implemented for - * dim=3, only. + * (standard orientation) and has + * no effect in 2d. */ static unsigned int face_to_cell_lines (const unsigned int face, const unsigned int line, @@ -484,9 +546,7 @@ GeometryInfo<1>::unit_cell_vertex (const unsigned int vertex) Assert (vertex < vertices_per_cell, ExcIndexRange (vertex, 0, vertices_per_cell)); - static const Point<1> vertices[vertices_per_cell] = - { Point<1>(0.), Point<1>(1.) }; - return vertices[vertex]; + return Point<1>(static_cast(vertex)); } @@ -499,10 +559,7 @@ GeometryInfo<2>::unit_cell_vertex (const unsigned int vertex) Assert (vertex < vertices_per_cell, ExcIndexRange (vertex, 0, vertices_per_cell)); - static const Point<2> vertices[vertices_per_cell] = - { Point<2>(0., 0.), Point<2>(1., 0.), - Point<2>(1.,1.), Point<2>(0.,1.) }; - return vertices[vertex]; + return Point<2>(vertex%2, vertex/2); } @@ -515,15 +572,11 @@ GeometryInfo<3>::unit_cell_vertex (const unsigned int vertex) Assert (vertex < vertices_per_cell, ExcIndexRange (vertex, 0, vertices_per_cell)); - static const Point<3> vertices[vertices_per_cell] = - { Point<3>(0., 0., 0.), Point<3>(1., 0., 0.), - Point<3>(1., 0., 1.), Point<3>(0., 0., 1.), - Point<3>(0., 1., 0.), Point<3>(1., 1., 0.), - Point<3>(1., 1., 1.), Point<3>(0., 1., 1.) }; - return vertices[vertex]; + return Point<3>(vertex%2, vertex/2%2, vertex/4); } + template inline Point @@ -557,8 +610,8 @@ GeometryInfo<2>::child_cell_from_point (const Point<2> &p) Assert ((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1])); return (p[0] <= 0.5 ? - (p[1] <= 0.5 ? 0 : 3) : - (p[1] <= 0.5 ? 1 : 2)); + (p[1] <= 0.5 ? 0 : 2) : + (p[1] <= 0.5 ? 1 : 3)); } @@ -574,11 +627,11 @@ GeometryInfo<3>::child_cell_from_point (const Point<3> &p) return (p[0] <= 0.5 ? (p[1] <= 0.5 ? - (p[2] <= 0.5 ? 0 : 3) : - (p[2] <= 0.5 ? 4 : 7)) : + (p[2] <= 0.5 ? 0 : 4) : + (p[2] <= 0.5 ? 2 : 6)) : (p[1] <= 0.5 ? - (p[2] <= 0.5 ? 1 : 2) : - (p[2] <= 0.5 ? 5 : 6))); + (p[2] <= 0.5 ? 1 : 5) : + (p[2] <= 0.5 ? 3 : 7))); } diff --git a/deal.II/base/source/data_out_base.cc b/deal.II/base/source/data_out_base.cc index dd2c0a0fae..298099df6a 100644 --- a/deal.II/base/source/data_out_base.cc +++ b/deal.II/base/source/data_out_base.cc @@ -72,8 +72,9 @@ template bool DataOutBase::Patch::operator == (const Patch &patch) const { + const double eps=3e-16; for (unsigned int i=0; i::vertices_per_cell; ++i) - if (vertices[i] != patch.vertices[i]) + if (vertices[i].distance(patch.vertices[i]) > eps) return false; for (unsigned int i=0; i::faces_per_cell; ++i) diff --git a/deal.II/base/source/geometry_info.cc b/deal.II/base/source/geometry_info.cc index 6993809624..52d3ead7e1 100644 --- a/deal.II/base/source/geometry_info.cc +++ b/deal.II/base/source/geometry_info.cc @@ -49,30 +49,49 @@ namespace internal +const unsigned int +GeometryInfoBase::unit_normal_direction[GeometryInfoBase::faces_per_cell_max_dim] += { 0, 0, 1, 1, 2, 2, 3, 3 }; + +const int +GeometryInfoBase::unit_normal_orientation[GeometryInfoBase::faces_per_cell_max_dim] += { -1, 1, -1, 1, -1, 1, -1, 1 }; + +const unsigned int +GeometryInfoBase::opposite_face[GeometryInfoBase::faces_per_cell_max_dim] += { 1, 0, 3, 2, 5, 4, 7, 6 }; -template <> -const unsigned int GeometryInfo<1>::opposite_face[GeometryInfo<1>::faces_per_cell] -= { 0, 1 }; template <> -const unsigned int GeometryInfo<2>::opposite_face[GeometryInfo<2>::faces_per_cell] -= { 2, 3, 0, 1 }; +const unsigned int GeometryInfo<1>::ucd_to_deal[GeometryInfo<1>::vertices_per_cell] += { 0, 1}; template <> -const unsigned int GeometryInfo<3>::opposite_face[GeometryInfo<3>::faces_per_cell] -= { 1, 0, 4, 5, 2, 3 }; +const unsigned int GeometryInfo<2>::ucd_to_deal[GeometryInfo<2>::vertices_per_cell] += { 0, 1, 3, 2}; template <> -const unsigned int GeometryInfo<4>::opposite_face[GeometryInfo<4>::faces_per_cell] -= { invalid_unsigned_int, - invalid_unsigned_int, - invalid_unsigned_int, - invalid_unsigned_int, - invalid_unsigned_int, - invalid_unsigned_int, - invalid_unsigned_int, - invalid_unsigned_int }; +const unsigned int GeometryInfo<3>::ucd_to_deal[GeometryInfo<3>::vertices_per_cell] += { 0, 1, 5, 4, 2, 3, 7, 6}; +template <> +const unsigned int GeometryInfo<4>::ucd_to_deal[GeometryInfo<4>::vertices_per_cell] += { invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int}; template <> @@ -81,70 +100,30 @@ const unsigned int GeometryInfo<1>::dx_to_deal[GeometryInfo<1>::vertices_per_cel template <> const unsigned int GeometryInfo<2>::dx_to_deal[GeometryInfo<2>::vertices_per_cell] -= { 0, 3, 1, 2}; += { 0, 2, 1, 3}; template <> const unsigned int GeometryInfo<3>::dx_to_deal[GeometryInfo<3>::vertices_per_cell] -= { 0, 3, 4, 7, 1, 2, 5, 6}; += { 0, 4, 2, 6, 1, 5, 3, 7}; template <> const unsigned int GeometryInfo<4>::dx_to_deal[GeometryInfo<4>::vertices_per_cell] -= { invalid_unsigned_int, - invalid_unsigned_int, - invalid_unsigned_int, - invalid_unsigned_int, - invalid_unsigned_int, - invalid_unsigned_int, - invalid_unsigned_int, - invalid_unsigned_int, - invalid_unsigned_int, - invalid_unsigned_int, - invalid_unsigned_int, - invalid_unsigned_int, - invalid_unsigned_int, - invalid_unsigned_int, - invalid_unsigned_int, - invalid_unsigned_int }; - - -//TODO: Use these values in mappings. -template <> -const unsigned int -GeometryInfo<1>::unit_normal_direction[GeometryInfo<1>::faces_per_cell] -= { 0, 0 }; -template <> -const int -GeometryInfo<1>::unit_normal_orientation[GeometryInfo<1>::faces_per_cell] -= { -1, 1 }; - -template <> -const unsigned int -GeometryInfo<2>::unit_normal_direction[GeometryInfo<2>::faces_per_cell] -= { 1, 0, 1, 0 }; -template <> -const int -GeometryInfo<2>::unit_normal_orientation[GeometryInfo<2>::faces_per_cell] -= { -1, 1, 1, -1 }; - -template <> -const unsigned int -GeometryInfo<3>::unit_normal_direction[GeometryInfo<3>::faces_per_cell] -= { 1, 1, 2, 0, 2, 0 }; -template <> -const int -GeometryInfo<3>::unit_normal_orientation[GeometryInfo<3>::faces_per_cell] -= { -1, 1, -1, 1, 1, -1 }; - -// These are already for the new, regular numbering -template <> -const unsigned int -GeometryInfo<4>::unit_normal_direction[GeometryInfo<4>::faces_per_cell] -= { 0, 0, 1, 1, 2, 2, 3, 3 }; -template <> -const int -GeometryInfo<4>::unit_normal_orientation[GeometryInfo<4>::faces_per_cell] -= { -1, 1, -1, 1, -1, 1, -1, 1 }; - += { invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int, + invalid_unsigned_int}; template <> @@ -171,10 +150,11 @@ GeometryInfo<2>::child_cell_on_face (const unsigned int face, Assert (subface::child_cell_on_face (const unsigned int face, Assert (face::child_cell_on_face (const unsigned int face, template <> unsigned int GeometryInfo<4>::child_cell_on_face (const unsigned int, - const unsigned int, - const bool) + const unsigned int, + const bool) { Assert(false, ExcNotImplemented()); return invalid_unsigned_int; @@ -243,19 +223,19 @@ GeometryInfo<3>::line_to_cell_vertices (const unsigned int line, const unsigned int vertex) { Assert (line::line_to_cell_vertices (const unsigned int line, template <> unsigned int GeometryInfo<4>::line_to_cell_vertices (const unsigned int, - const unsigned int) + const unsigned int) { Assert(false, ExcNotImplemented()); return invalid_unsigned_int; @@ -312,24 +292,26 @@ GeometryInfo<3>::face_to_cell_lines (const unsigned int face, Assert (face +template unsigned int -GeometryInfo<4>::face_to_cell_lines (const unsigned int, - const unsigned int, - const bool) +GeometryInfo::face_to_cell_lines (const unsigned int, + const unsigned int, + const bool) { Assert(false, ExcNotImplemented()); return invalid_unsigned_int; diff --git a/deal.II/base/source/quadrature.cc b/deal.II/base/source/quadrature.cc index 7d194347d8..82ea10678c 100644 --- a/deal.II/base/source/quadrature.cc +++ b/deal.II/base/source/quadrature.cc @@ -391,16 +391,16 @@ QProjector<2>::project_to_face (const Quadrature<1> &quadrature, switch (face_no) { case 0: - q_points[p] = Point(quadrature.point(p)(0),0); + q_points[p] = Point(0,quadrature.point(p)(0)); break; case 1: q_points[p] = Point(1,quadrature.point(p)(0)); break; case 2: - q_points[p] = Point(quadrature.point(p)(0),1); + q_points[p] = Point(quadrature.point(p)(0),0); break; case 3: - q_points[p] = Point(0,quadrature.point(p)(0)); + q_points[p] = Point(quadrature.point(p)(0),1); break; default: Assert (false, ExcInternalError()); @@ -424,35 +424,35 @@ QProjector<3>::project_to_face (const Quadrature<2> &quadrature, switch (face_no) { case 0: - q_points[p] = Point(quadrature.point(p)(0), - 0, + q_points[p] = Point(0, + quadrature.point(p)(0), quadrature.point(p)(1)); break; case 1: - q_points[p] = Point(quadrature.point(p)(0), - 1, + q_points[p] = Point(1, + quadrature.point(p)(0), quadrature.point(p)(1)); break; case 2: - q_points[p] = Point(quadrature.point(p)(0), - quadrature.point(p)(1), - 0); + q_points[p] = Point(quadrature.point(p)(1), + 0, + quadrature.point(p)(0)); break; case 3: - q_points[p] = Point(1, - quadrature.point(p)(0), - quadrature.point(p)(1)); + q_points[p] = Point(quadrature.point(p)(1), + 1, + quadrature.point(p)(0)); break; case 4: q_points[p] = Point(quadrature.point(p)(0), quadrature.point(p)(1), - 1); + 0); break; case 5: - q_points[p] = Point(0, - quadrature.point(p)(0), - quadrature.point(p)(1)); - break; + q_points[p] = Point(quadrature.point(p)(0), + quadrature.point(p)(1), + 1); + break; default: Assert (false, ExcInternalError()); @@ -493,12 +493,10 @@ QProjector<2>::project_to_subface (const Quadrature<1> &quadrature, switch (subface_no) { case 0: - q_points[p] - = Point(quadrature.point(p)(0)/2,0); + q_points[p] = Point(0,quadrature.point(p)(0)/2); break; case 1: - q_points[p] - = Point(quadrature.point(p)(0)/2+0.5,0); + q_points[p] = Point(0,quadrature.point(p)(0)/2+0.5); break; default: Assert (false, ExcInternalError()); @@ -521,10 +519,12 @@ QProjector<2>::project_to_subface (const Quadrature<1> &quadrature, switch (subface_no) { case 0: - q_points[p] = Point(quadrature.point(p)(0)/2,1); + q_points[p] + = Point(quadrature.point(p)(0)/2,0); break; case 1: - q_points[p] = Point(quadrature.point(p)(0)/2+0.5,1); + q_points[p] + = Point(quadrature.point(p)(0)/2+0.5,0); break; default: Assert (false, ExcInternalError()); @@ -534,15 +534,16 @@ QProjector<2>::project_to_subface (const Quadrature<1> &quadrature, switch (subface_no) { case 0: - q_points[p] = Point(0,quadrature.point(p)(0)/2); + q_points[p] = Point(quadrature.point(p)(0)/2,1); break; case 1: - q_points[p] = Point(0,quadrature.point(p)(0)/2+0.5); + q_points[p] = Point(quadrature.point(p)(0)/2+0.5,1); break; default: Assert (false, ExcInternalError()); }; break; + default: Assert (false, ExcInternalError()); }; @@ -572,44 +573,43 @@ QProjector<3>::project_to_subface (const Quadrature<2> &quadrature, switch (face_no) { case 0: - q_points[p] = Point(quadrature.point(p)(0)/2, - 0, + q_points[p] = Point(0, + quadrature.point(p)(0)/2, quadrature.point(p)(1)/2); switch (subface_no) { case 0: break; case 1: - q_points[p][0] += 1./2.; + q_points[p][1] += 1./2.; break; case 2: - q_points[p][0] += 1./2.; q_points[p][2] += 1./2.; break; case 3: + q_points[p][1] += 1./2.; q_points[p][2] += 1./2.; break; default: Assert (false, ExcInternalError()); }; - break; case 1: - q_points[p] = Point(quadrature.point(p)(0)/2, - 1, + q_points[p] = Point(1, + quadrature.point(p)(0)/2, quadrature.point(p)(1)/2); switch (subface_no) { case 0: break; case 1: - q_points[p][0] += 1./2.; + q_points[p][1] += 1./2.; break; case 2: - q_points[p][0] += 1./2.; q_points[p][2] += 1./2.; break; case 3: + q_points[p][1] += 1./2.; q_points[p][2] += 1./2.; break; default: @@ -617,43 +617,43 @@ QProjector<3>::project_to_subface (const Quadrature<2> &quadrature, }; break; case 2: - q_points[p] = Point(quadrature.point(p)(0)/2, - quadrature.point(p)(1)/2, - 0); + q_points[p] = Point(quadrature.point(p)(1)/2, + 0, + quadrature.point(p)(0)/2); switch (subface_no) { case 0: break; case 1: - q_points[p][0] += 1./2.; + q_points[p][2] += 1./2.; break; case 2: q_points[p][0] += 1./2.; - q_points[p][1] += 1./2.; break; case 3: - q_points[p][1] += 1./2.; + q_points[p][0] += 1./2.; + q_points[p][2] += 1./2.; break; default: Assert (false, ExcInternalError()); }; break; case 3: - q_points[p] = Point(1, - quadrature.point(p)(0)/2, - quadrature.point(p)(1)/2); + q_points[p] = Point(quadrature.point(p)(1)/2, + 1, + quadrature.point(p)(0)/2); switch (subface_no) { case 0: break; case 1: - q_points[p][1] += 1./2.; + q_points[p][2] += 1./2.; break; case 2: - q_points[p][1] += 1./2.; - q_points[p][2] += 1./2.; + q_points[p][0] += 1./2.; break; case 3: + q_points[p][0] += 1./2.; q_points[p][2] += 1./2.; break; default: @@ -663,7 +663,7 @@ QProjector<3>::project_to_subface (const Quadrature<2> &quadrature, case 4: q_points[p] = Point(quadrature.point(p)(0)/2, quadrature.point(p)(1)/2, - 1); + 0); switch (subface_no) { case 0: @@ -672,10 +672,10 @@ QProjector<3>::project_to_subface (const Quadrature<2> &quadrature, q_points[p][0] += 1./2.; break; case 2: - q_points[p][0] += 1./2.; q_points[p][1] += 1./2.; break; case 3: + q_points[p][0] += 1./2.; q_points[p][1] += 1./2.; break; default: @@ -683,27 +683,27 @@ QProjector<3>::project_to_subface (const Quadrature<2> &quadrature, }; break; case 5: - q_points[p] = Point(0, - quadrature.point(p)(0)/2, - quadrature.point(p)(1)/2); + q_points[p] = Point(quadrature.point(p)(0)/2, + quadrature.point(p)(1)/2, + 1); switch (subface_no) { case 0: break; case 1: - q_points[p][1] += 1./2.; + q_points[p][0] += 1./2.; break; case 2: q_points[p][1] += 1./2.; - q_points[p][2] += 1./2.; break; case 3: - q_points[p][2] += 1./2.; + q_points[p][0] += 1./2.; + q_points[p][1] += 1./2.; break; default: Assert (false, ExcInternalError()); }; - break; + break; default: Assert (false, ExcInternalError()); }; diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index 1d33fe5d29..0a1bd5e55d 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -160,15 +160,15 @@ template class FECollection; * from the following sketch, where the overall order of the different * dof groups is depicted: * @verbatim - * *--13--3--14--* + * *--15--4--16--* * | | | - * 16 20 7 19 12 + * 10 19 6 20 12 * | | | - * 4--8---0--6---2 + * 1--7---0--8---2 * | | | - * 15 17 5 18 11 + * 9 17 5 18 11 * | | | - * *--9---1--10--* + * *--13--3--14--* * @endverbatim * The numbering of vertices and lines, as well as the numbering of * children within a line is consistent with the one described in diff --git a/deal.II/deal.II/include/fe/fe_poly.h b/deal.II/deal.II/include/fe/fe_poly.h index bffcdc241d..075e284855 100644 --- a/deal.II/deal.II/include/fe/fe_poly.h +++ b/deal.II/deal.II/include/fe/fe_poly.h @@ -49,8 +49,8 @@ * Furthermore, this class assumes that shape functions of the * FiniteElement under consideration do not depend on the * actual shape of the cells in real space, i.e. update_once() - * includes update_values. For FiniteElements whose shape functions - * depend on the cells in real space, the update_once() and + * includes update_values. For FiniteElements whose shape + * functions depend on the cells in real space, the update_once() and * update_each() functions must be overloaded. * * Todos: diff --git a/deal.II/deal.II/include/fe/fe_q.h b/deal.II/deal.II/include/fe/fe_q.h index c0ad6aec70..5d05291c94 100644 --- a/deal.II/deal.II/include/fe/fe_q.h +++ b/deal.II/deal.II/include/fe/fe_q.h @@ -82,7 +82,7 @@ * *
  • 2D case: * @verbatim - * 3-------2 + * 2-------3 * | | * | | * | | @@ -91,12 +91,12 @@ * *
  • 3D case: * @verbatim - * 7-------6 7-------6 + * 6-------7 6-------7 * /| | / /| * / | | / / | * / | | / / | - * 3 | | 3-------2 | - * | 4-------5 | | 5 + * 4 | | 4-------5 | + * | 2-------3 | | 3 * | / / | | / * | / / | | / * |/ / | |/ @@ -108,12 +108,12 @@ *
      *
    • Index 0: [0, 0, 0]; *
    • Index 1: [1, 0, 0]; - *
    • Index 2: [1, 0, 1]; - *
    • Index 3: [0, 0, 1]; - *
    • Index 4: [0, 1, 0]; - *
    • Index 5: [1, 1, 0]; - *
    • Index 6: [1, 1, 1]; - *
    • Index 7: [0, 1, 1]; + *
    • Index 2: [0, 1, 0]; + *
    • Index 3: [1, 1, 0]; + *
    • Index 4: [0, 0, 1]; + *
    • Index 5: [1, 0, 1]; + *
    • Index 6: [0, 1, 1]; + *
    • Index 7: [1, 1, 1]; *
    * *

    Q2 elements

    @@ -125,34 +125,34 @@ * *
  • 2D case: * @verbatim - * 3---6---2 + * 2---7---3 * | | - * 7 8 5 + * 4 8 5 * | | - * 0---4---1 + * 0---6---1 * @endverbatim * *
  • 3D case: * @verbatim - * 7--14---6 7--14---6 + * 6--15---7 6--15---7 * /| | / /| - * 19 | 13 19 1813 - * / 15 | / / | - * 3 | | 3---10--2 | - * | 4--12---5 | | 5 - * | / / | 9 / - * 11 16 17 11 | 17 + * 12 | 19 12 1319 + * / 18 | / / | + * 4 | | 4---14--5 | + * | 2---11--3 | | 3 + * | / / | 17 / + * 16 8 9 16 | 9 * |/ / | |/ - * 0---8---1 0---8---1 + * 0---10--1 0---8---1 * * *-------* *-------* * /| | / /| - * / | 21 | / 24 / | + * / | 23 | / 25 / | * / | | / / | * * | | *-------* | - * |25 *-------* | |23 * - * | / / | 20 | / - * | / 22 / | | / + * |20 *-------* | |21 * + * | / / | 22 | / + * | / 24 / | | / * |/ / | |/ * *-------* *-------* * @endverbatim @@ -163,30 +163,30 @@ *
      *
    • Index 0: [0, 0, 0]; *
    • Index 1: [1, 0, 0]; - *
    • Index 2: [1, 0, 1]; - *
    • Index 3: [0, 0, 1]; - *
    • Index 4: [0, 1, 0]; - *
    • Index 5: [1, 1, 0]; - *
    • Index 6: [1, 1, 1]; - *
    • Index 7: [0, 1, 1]; - *
    • Index 8: [1/2, 0, 0]; - *
    • Index 9: [1, 0, 1/2]; - *
    • Index 10: [1/2, 0, 1]; - *
    • Index 11: [0, 0, 1/2]; - *
    • Index 12: [1/2, 1, 0]; - *
    • Index 13: [1, 1, 1/2]; - *
    • Index 14: [1/2, 1, 1]; - *
    • Index 15: [0, 1, 1/2]; - *
    • Index 16: [0, 1/2, 0]; - *
    • Index 17: [1, 1/2, 0]; - *
    • Index 18: [1, 1/2, 1]; - *
    • Index 19: [0, 1/2, 1]; - *
    • Index 20: [1/2, 0, 1/2]; - *
    • Index 21: [1/2, 1, 1/2]; - *
    • Index 22: [1/2, 1/2, 0]; - *
    • Index 23: [1, 1/2, 1/2]; - *
    • Index 24: [1/2, 1/2, 1]; - *
    • Index 25: [0, 1/2, 1/2]; + *
    • Index 2: [0, 1, 0]; + *
    • Index 3: [1, 1, 0]; + *
    • Index 4: [0, 0, 1]; + *
    • Index 5: [1, 0, 1]; + *
    • Index 6: [0, 1, 1]; + *
    • Index 7: [1, 1, 1]; + *
    • Index 8: [0, 1/2, 0]; + *
    • Index 9: [1, 1/2, 0]; + *
    • Index 10: [1/2, 0, 0]; + *
    • Index 11: [1/2, 1, 0]; + *
    • Index 12: [0, 1/2, 1]; + *
    • Index 13: [1, 1/2, 1]; + *
    • Index 14: [1/2, 0, 1]; + *
    • Index 15: [1/2, 1, 1]; + *
    • Index 16: [0, 0, 1/2]; + *
    • Index 17: [1, 0, 1/2]; + *
    • Index 18: [0, 1, 1/2]; + *
    • Index 19: [1, 1, 1/2]; + *
    • Index 20: [0, 1/2, 1/2]; + *
    • Index 21: [1, 1/2, 1/2]; + *
    • Index 22: [1/2, 0, 1/2]; + *
    • Index 23: [1/2, 1, 1/2]; + *
    • Index 24: [1/2, 1/2, 0]; + *
    • Index 25: [1/2, 1/2, 1]; *
    • Index 26: [1/2, 1/2, 1/2]; *
    * @@ -199,16 +199,14 @@ * *
  • 2D case: * @verbatim - * 3--8--9--2 + * 2--10-11-3 * | | - * 11 14 15 7 + * 5 14 15 7 * | | - * 10 12 13 6 + * 4 12 13 6 * | | - * 0--4--5--1 + * 0--8--9--1 * @endverbatim - * Note the reverse ordering of degrees of freedom on the left and - * upper line. * *

    Q4 elements

    *
      @@ -219,21 +217,19 @@ * *
    • 2D case: * @verbatim - * 3--10-11-12-2 + * 2--13-14-15-3 * | | - * 15 22 23 24 9 + * 6 22 23 24 9 * | | - * 14 19 20 21 8 + * 5 19 20 21 8 * | | - * 13 16 17 18 7 + * 4 16 17 18 7 * | | - * 0--4--5--6--1 + * 0--10-11-12-1 * @endverbatim *
    - * Note the reverse ordering of degrees of freedom on the left and upper - * line. * - * @author Wolfgang Bangerth, 1998, 2003; Guido Kanschat, 2001; Ralf Hartmann, 2001, 2004; Oliver Kayser-Herold, 2004 + * @author Wolfgang Bangerth, 1998, 2003; Guido Kanschat, 2001; Ralf Hartmann, 2001, 2004, 2005; Oliver Kayser-Herold, 2004 */ template class FE_Q : public FE_Poly,dim> diff --git a/deal.II/deal.II/include/fe/fe_q_hierarchical.h b/deal.II/deal.II/include/fe/fe_q_hierarchical.h index 6e43e9f023..af16fac6d1 100644 --- a/deal.II/deal.II/include/fe/fe_q_hierarchical.h +++ b/deal.II/deal.II/include/fe/fe_q_hierarchical.h @@ -78,7 +78,7 @@ template class MappingQ; * *
  • 2D case: * @verbatim - * 3-------2 + * 2-------3 * | | * | | * | | @@ -87,28 +87,29 @@ template class MappingQ; * *
  • 3D case: * @verbatim - * 7-------6 7-------6 + * 6-------7 6-------7 * /| | / /| * / | | / / | * / | | / / | - * 3 | | 3-------2 | - * | 4-------5 | | 5 + * 4 | | 4-------5 | + * | 2-------3 | | 3 * | / / | | / * | / / | | / * |/ / | |/ * 0-------1 0-------1 + * @endverbatim * * The respective coordinate values of the support points of the degrees * of freedom are as follows: *
      *
    • Index 0: [0, 0, 0]; *
    • Index 1: [1, 0, 0]; - *
    • Index 2: [1, 0, 1]; - *
    • Index 3: [0, 0, 1]; - *
    • Index 4: [0, 1, 0]; - *
    • Index 5: [1, 1, 0]; - *
    • Index 6: [1, 1, 1]; - *
    • Index 7: [0, 1, 1]; + *
    • Index 2: [0, 1, 0]; + *
    • Index 3: [1, 1, 0]; + *
    • Index 4: [0, 0, 1]; + *
    • Index 5: [1, 0, 1]; + *
    • Index 6: [0, 1, 1]; + *
    • Index 7: [1, 1, 1]; *
    * *

    Q2 elements

    @@ -120,34 +121,34 @@ template class MappingQ; * *
  • 2D case: * @verbatim - * 3---6---2 + * 2---7---3 * | | - * 7 8 5 + * 4 8 5 * | | - * 0---4---1 + * 0---6---1 * @endverbatim * *
  • 3D case: * @verbatim - * 7--14---6 7--14---6 + * 6--15---7 6--15---7 * /| | / /| - * 19 | 13 19 1813 - * / 15 | / / | - * 3 | | 3---10--2 | - * | 4--12---5 | | 5 - * | / / | 9 / - * 11 16 17 11 | 17 + * 12 | 19 12 1319 + * / 18 | / / | + * 4 | | 4---14--5 | + * | 2---11--3 | | 3 + * | / / | 17 / + * 16 8 9 16 | 9 * |/ / | |/ - * 0---8---1 0---8---1 + * 0---10--1 0---8---1 * * *-------* *-------* * /| | / /| - * / | 21 | / 24 / | + * / | 23 | / 25 / | * / | | / / | * * | | *-------* | - * |25 *-------* | |23 * - * | / / | 20 | / - * | / 22 / | | / + * |20 *-------* | |21 * + * | / / | 22 | / + * | / 24 / | | / * |/ / | |/ * *-------* *-------* * @endverbatim @@ -158,30 +159,30 @@ template class MappingQ; *
      *
    • Index 0: [0, 0, 0]; *
    • Index 1: [1, 0, 0]; - *
    • Index 2: [1, 0, 1]; - *
    • Index 3: [0, 0, 1]; - *
    • Index 4: [0, 1, 0]; - *
    • Index 5: [1, 1, 0]; - *
    • Index 6: [1, 1, 1]; - *
    • Index 7: [0, 1, 1]; - *
    • Index 8: [1/2, 0, 0]; - *
    • Index 9: [1, 0, 1/2]; - *
    • Index 10: [1/2, 0, 1]; - *
    • Index 11: [0, 0, 1/2]; - *
    • Index 12: [1/2, 1, 0]; - *
    • Index 13: [1, 1, 1/2]; - *
    • Index 14: [1/2, 1, 1]; - *
    • Index 15: [0, 1, 1/2]; - *
    • Index 16: [0, 1/2, 0]; - *
    • Index 17: [1, 1/2, 0]; - *
    • Index 18: [1, 1/2, 1]; - *
    • Index 19: [0, 1/2, 1]; - *
    • Index 20: [1/2, 0, 1/2]; - *
    • Index 21: [1/2, 1, 1/2]; - *
    • Index 22: [1/2, 1/2, 0]; - *
    • Index 23: [1, 1/2, 1/2]; - *
    • Index 24: [1/2, 1/2, 1]; - *
    • Index 25: [0, 1/2, 1/2]; + *
    • Index 2: [0, 1, 0]; + *
    • Index 3: [1, 1, 0]; + *
    • Index 4: [0, 0, 1]; + *
    • Index 5: [1, 0, 1]; + *
    • Index 6: [0, 1, 1]; + *
    • Index 7: [1, 1, 1]; + *
    • Index 8: [0, 1/2, 0]; + *
    • Index 9: [1, 1/2, 0]; + *
    • Index 10: [1/2, 0, 0]; + *
    • Index 11: [1/2, 1, 0]; + *
    • Index 12: [0, 1/2, 1]; + *
    • Index 13: [1, 1/2, 1]; + *
    • Index 14: [1/2, 0, 1]; + *
    • Index 15: [1/2, 1, 1]; + *
    • Index 16: [0, 0, 1/2]; + *
    • Index 17: [1, 0, 1/2]; + *
    • Index 18: [0, 1, 1/2]; + *
    • Index 19: [1, 1, 1/2]; + *
    • Index 20: [0, 1/2, 1/2]; + *
    • Index 21: [1, 1/2, 1/2]; + *
    • Index 22: [1/2, 0, 1/2]; + *
    • Index 23: [1/2, 1, 1/2]; + *
    • Index 24: [1/2, 1/2, 0]; + *
    • Index 25: [1/2, 1/2, 1]; *
    • Index 26: [1/2, 1/2, 1/2]; *
    * @@ -194,16 +195,14 @@ template class MappingQ; * *
  • 2D case: * @verbatim - * 3--8--9--2 + * 2--10-11-3 * | | - * 11 14 15 7 + * 5 14 15 7 * | | - * 10 12 13 6 + * 4 12 13 6 * | | - * 0--4--5--1 + * 0--8--9--1 * @endverbatim - * Note the reverse ordering of degrees of freedom on the left and - * upper line. * *

    Q4 elements

    *
      @@ -214,19 +213,17 @@ template class MappingQ; * *
    • 2D case: * @verbatim - * 3--10-11-12-2 + * 2--13-14-15-3 * | | - * 15 22 23 24 9 + * 6 22 23 24 9 * | | - * 14 19 20 21 8 + * 5 19 20 21 8 * | | - * 13 16 17 18 7 + * 4 16 17 18 7 * | | - * 0--4--5--6--1 + * 0--10-11-12-1 * @endverbatim *
    - * Note the reverse ordering of degrees of freedom on the left and upper - * line. * * @author Brian Carnes, 2002, Ralf Hartmann 2004, 2005 */ diff --git a/deal.II/deal.II/include/fe/mapping.h b/deal.II/deal.II/include/fe/mapping.h index 144c3290be..157d7c74b4 100644 --- a/deal.II/deal.II/include/fe/mapping.h +++ b/deal.II/deal.II/include/fe/mapping.h @@ -360,30 +360,6 @@ class Mapping : public Subscriptor */ DeclException0 (ExcInvalidData); - - - protected: - - /** - * Vector of unit normal - * directions. The entry divided by - * 2 determines the non-zero - * component of the normal vector: - * 0 means x, 1 means y and 2 means - * z. The entry modulo 2 determines - * the orientation of the first - * tangential vector in the - * cross-product. This has to be - * chosen such that the normal - * vector points outwards. - * - * This variable is purely for - * internal use and its values are - * determined by its usage in the - * source code. - */ - static const unsigned int normal_directions[2*dim]; - private: /** @@ -510,17 +486,6 @@ class Mapping : public Subscriptor /*@}*/ -/* -------------- declaration of explicit specializations ------------- */ - - -// declaration of explicit specializations of member variables, if the -// compiler allows us to do that (the standard says we must) -#ifndef DEAL_II_MEMBER_VAR_SPECIALIZATION_BUG -template<> const unsigned int Mapping<1>::normal_directions[2]; -template<> const unsigned int Mapping<2>::normal_directions[4]; -template<> const unsigned int Mapping<3>::normal_directions[6]; -#endif - /* ------------------------- inline functions ------------------------- */ diff --git a/deal.II/deal.II/include/fe/mapping_cartesian.h b/deal.II/deal.II/include/fe/mapping_cartesian.h index 716ca732a8..dde71233a3 100644 --- a/deal.II/deal.II/include/fe/mapping_cartesian.h +++ b/deal.II/deal.II/include/fe/mapping_cartesian.h @@ -32,7 +32,7 @@ * cartesian meshes. Apply this mapping to a general mesh to get * strange results. * - * @author Guido Kanschat, 2001 + * @author Guido Kanschat, 2001; Ralf Hartmann, 2005 */ template class MappingCartesian : public Mapping diff --git a/deal.II/deal.II/include/fe/mapping_q.h b/deal.II/deal.II/include/fe/mapping_q.h index 8de37b52b1..135e855400 100644 --- a/deal.II/deal.II/include/fe/mapping_q.h +++ b/deal.II/deal.II/include/fe/mapping_q.h @@ -33,7 +33,7 @@ template class TensorProductPolynomials; * deal.II/doc/reports/mapping_q/index.html in the `Reports' * section of `Documentation'. * - * @author Ralf Hartmann, Guido Kanschat 2000, 2001 + * @author Ralf Hartmann, 2000, 2001, 2005; Guido Kanschat 2000, 2001 */ template class MappingQ : public MappingQ1 diff --git a/deal.II/deal.II/include/fe/mapping_q1.h b/deal.II/deal.II/include/fe/mapping_q1.h index 286a826e7b..d02932a4ec 100644 --- a/deal.II/deal.II/include/fe/mapping_q1.h +++ b/deal.II/deal.II/include/fe/mapping_q1.h @@ -41,7 +41,7 @@ * element @p FE_Q of order 1. Therefore, coupling these two yields * an isoparametric element. * - * @author Guido Kanschat, Ralf Hartmann, 2000, 2001 + * @author Guido Kanschat, 2000, 2001; Ralf Hartmann, 2000, 2001, 2005 */ template class MappingQ1 : public Mapping @@ -550,13 +550,6 @@ class MappingQ1 : public Mapping const Point &p, InternalData &mdata, Point &p_unit) const; - - /** - * Mapping between tensor product - * ordering and real ordering of - * the vertices. - */ - static const unsigned int vertex_mapping[1<::InternalData::derivative (const unsigned int qpoint, #ifndef DOXYGEN -// declaration of explicit specializations of member variables, if the -// compiler allows us to do that (the standard says we must) -#ifndef DEAL_II_MEMBER_VAR_SPECIALIZATION_BUG -template<> const unsigned int MappingQ1<1>::vertex_mapping[2]; -template<> const unsigned int MappingQ1<2>::vertex_mapping[4]; -template<> const unsigned int MappingQ1<3>::vertex_mapping[8]; -#endif template<> void MappingQ1<1>::compute_shapes_virtual ( const std::vector > &unit_points, InternalData& data) const; diff --git a/deal.II/deal.II/include/grid/persistent_tria.h b/deal.II/deal.II/include/grid/persistent_tria.h index 7f07f46ca3..9d3db9b630 100644 --- a/deal.II/deal.II/include/grid/persistent_tria.h +++ b/deal.II/deal.II/include/grid/persistent_tria.h @@ -217,7 +217,17 @@ class PersistentTriangulation : public Triangulation */ virtual void create_triangulation (const std::vector > &vertices, const std::vector > &cells, - const SubCellData &subcelldata); + const SubCellData &subcelldata); + + /** + * Throw an error, since this function + * is not useful in the context of this + * class. + */ + virtual void create_triangulation_compatibility ( + const std::vector > &vertices, + const std::vector > &cells, + const SubCellData &subcelldata); /** * Writes all refine and coarsen diff --git a/deal.II/deal.II/include/grid/tria.h b/deal.II/deal.II/include/grid/tria.h index e310389848..b693d613a0 100644 --- a/deal.II/deal.II/include/grid/tria.h +++ b/deal.II/deal.II/include/grid/tria.h @@ -1105,25 +1105,45 @@ struct TriaNumberCache<3> : public TriaNumberCache<2> * apply some smoothing for multigrid algorithms, but this has to be decided * upon later. * + * N4/ face lines in 3d are ordered, such that the induced 2d local + * coordinate system (x,y) implies (right hand rule) a normal in + * face normal direction, see N2/ * *

    Implementation conventions for two spatial dimensions

    + * + * From version 5.2 onwards deal.II is based on a numbering scheme, + * which uses a lexicographic ordering (with x running fastest) + * whereever possible, hence trying to adopt a kind of 'canonical' + * ordering. + * + * The ordering of vertices and faces (lines) in 2d is defined by + * + * N1) vertices are numbered in lexicographic ordering * - * There is a convention about the direction of the bounding lines of quads in - * 2D. The direction of a line is the direction of point 0 towards point 1. We - * define, that allowed cells contain of lines of which the direction is - * as follows: + * N2) faces (lines in 2d): first the two faces with normals in x- + * and then y-direction. For each two faces: first the face with + * normal in negative coordinate direction, then the one with normal + * in positive direction, i.e. the faces are ordered according to + * their normals pointing in -x, x, -y, y direction. + * + * N3) the direction of a line is represented by the direction of + * point 0 towards point 1 and is always in one of the coordinate + * directions + * + * The resulting numbering of vertices and faces (lines) in 2d as + * well as the directions of lines is shown in the following. * @verbatim + * 3 + * 2-->--3 + * | | + * 0^ ^1 + * | | + * 0-->--1 * 2 - * 3--->---2 - * | | - * 3^ ^1 - * | | - * 0--->---1 - * 0 * @endverbatim - * The number of the vertices and lines is also indicated. This orientation of - * lines has to be checked/generated upon construction of a grid and is - * preserved upon refinement. + * We note, that the orientation of lines has to be + * checked/generated upon construction of a grid and is preserved + * upon refinement. * * Further we define, that child lines have the same direction as their parent, * i.e. that subline(0).vertex(0)==line.vertex(0) and @@ -1144,15 +1164,16 @@ struct TriaNumberCache<3> : public TriaNumberCache<2> * quadrature formulae or the point of definition of trial functions), we * define the following coordinate system for the unit cell: * @verbatim - * y^ 3-------2 - * | | | - * | | | - * | | | - * | 0-------1 - * *-------------->x + * y^ 2-----3 + * | | | + * | | | + * | | | + * | 0-----1 + * *------------>x * @endverbatim - * with vertex 0 being the origin of the coordinate system, vertex 1 having - * coordinates (1,0), vertex 2 at (1,1) and vertex 3 at (0,1). + * with vertex 0 being the origin of the coordinate system, vertex 1 + * having coordinates (1,0), vertex 2 at (0,1) and + * vertex 3 at (1,1). * * *

    Implementation conventions for three spatial dimensions

    @@ -1161,18 +1182,18 @@ struct TriaNumberCache<3> : public TriaNumberCache<2> * for vertices, lines and faces of hexahedra in three space * dimensions. Before giving these conventions we declare the * following sketch to be the standard way of drawing 3d pictures of - * hexahedra: + * hexahedra: * @verbatim - * *-------* *-------* - * /| | / /| - * / | | / / | - * / | | / / | - * * | | *-------* | - * | *-------* | | * - * | / / | | / - * | / / | | / - * |/ / | |/ - * *-------* *-------* + * *-------* *-------* + * /| | / /| + * / | | / / | + * z / | | / / | + * ^ * | | *-------* | + * | ^y | *-------* | | * + * | / | / / | | / + * | / | / / | | / + * |/ |/ / | |/ + * *------>x *-------* *-------* * @endverbatim * The left part of the picture shows the left, bottom and back face of the * cube, while the right one shall be the top, right and front face. You may @@ -1182,45 +1203,55 @@ struct TriaNumberCache<3> : public TriaNumberCache<2> * conventions can be extracted at run- or compile-time from the * member functions and variables of the GeometryInfo classes. * - * *

    Vertices

    + * + * The ordering of vertices in 3d is defined by the same rules as in + * the 2d case, i.e. + * + * N1) vertices are numbered in lexicographic ordering * - * The vertices on the front face are numbered exactly the same way as are - * the vertices on a quadrilateral. The vertices on the back face are numbered - * similarly by moving the front face to the back (no turning, no twisting, - * just a shift): + * Hence, the vertices are numbered as follows * @verbatim - * 7-------6 7-------6 + * 6-------7 6-------7 * /| | / /| * / | | / / | * / | | / / | - * 3 | | 3-------2 | - * | 4-------5 | | 5 + * 4 | | 4-------5 | + * | 2-------3 | | 3 * | / / | | / * | / / | | / * |/ / | |/ * 0-------1 0-------1 * @endverbatim + * + * We note, that first the vertices on the bottom face (z=0) are + * numbered exactly the same way as are the vertices on a + * quadrilateral. Then the vertices on the top face (z=1) are + * numbered similarly by moving the bottom face to the top. * *

    Lines

    * - * Here, the same holds as for the vertices: the lines of the front face are - * numbered as for the quadrilateral, for the back face they are just shifted. - * Finally, the four lines connecting front and back face are numbered: + * Here, the same holds as for the vertices: + * + * N4) line ordering in 3d: + *
      + *
    • first the lines of face (z=0) in 2d line ordering, + *
    • then the lines of face (z=1) in 2d line ordering, + *
    • finally the lines in z direction in lexicographic ordering + *
    * @verbatim - * *---6---* *---6---* + * *---7---* *---7---* * /| | / /| - * 11 | 5 11 10 5 - * / 7 | / / | - * * | | *---2---* | - * | *---4---* | | * - * | / / | 1 / - * 3 8 9 3 | 9 + * 4 | 11 4 5 11 + * / 10 | / / | + * * | | *---6---* | + * | *---3---* | | * + * | / / | 9 / + * 8 0 1 8 | 1 * |/ / | |/ - * *---0---* *---0---* + * *---2---* *---2---* * @endverbatim - * The directions of the front and back lines is as for the respective faces, while - * the connecting lines always point to the back: + * As in 2d lines are directed in coordinate directions, see N3. * @verbatim * *--->---* *--->---* * /| | / /| @@ -1250,93 +1281,99 @@ struct TriaNumberCache<3> : public TriaNumberCache<2> * *

    Faces

    * - * The faces are numbered in the same order as the lines were numbered: front - * face, back face, then the four side faces: + * The numbering of faces in 3d is defined by a rule analogous to 2d: + * + * N2a) faces (quads in 3d): first the two faces with normals in x-, + * then y- and z-direction. For each two faces: first the face with + * normal in negative coordinate direction, then the one with normal + * in positive direction, i.e. the faces are ordered according to + * their normals pointing in -x, x, -y, y, -z, z direction. + * + * Therefore, the faces are numbered in the ordering: left, right, + * front, back, bottom and top face: * @verbatim * *-------* *-------* * /| | / /| - * / | 1 | / 4 / | + * / | 3 | / 5 / | * / | | / / | * * | | *-------* | - * | 5 *-------* | | 3 * + * | 0 *-------* | | 1 * * | / / | | / - * | / 2 / | 0 | / + * | / 4 / | 2 | / * |/ / | |/ * *-------* *-------* * @endverbatim * - * The standard direction of the faces is determined by the - * numbers the lines have within a given face. This is like follows: + * The standard direction of the faces is such, that the + * induced 2d local coordinate system (x,y) implies (right hand + * rule) a normal in face normal direction, see N2a). In the + * following we show the local coordinate system and the numbering + * of face lines: *
      *
    • Faces 0 and 1: * @verbatim - * *---2---* *-------* + * Face 0 Face 1 + * *-------* *-------* * /| | / /| - * / | 1 / / | - * / 3 | / / | - * * | | *---2---* | - * | *---0---* | | * - * | / / | 1 / - * | / / 3 | / + * 3 1 | / 3 1 + * y/ | | / y/ | + * * |x | *-------* |x + * | *-------* | | * + * 0 / / | 0 / + * | 2 / | | 2 * |/ / | |/ - * *-------* *---0---* + * *-------* *-------* * @endverbatim * - *
    • Faces 2 and 4: + *
    • Faces 2 and 3: * @verbatim - * *-------* *---2---* + * x Face 3 Face 2 + * *---1---* *-------* * /| | / /| - * / | | 3 1 | - * / | | / / | - * * | | *---0---* | - * | *---2---* | | * - * | / / | | / - * | 3 1 | | / + * / | 3 / / | + * / 2 | x/ / | + * * | | *---1---* | + * | *---0---*y | | * + * | / / | 3 / + * | / / 2 | / * |/ / | |/ - * *---0---* *-------* + * *-------* *---0---*y * @endverbatim * - *
    • Faces 3 and 5: + *
    • Faces 4 and 5: * @verbatim - * *-------* *-------* + * Face 4 y Face 5 + * *-------* *---3---* * /| | / /| - * 2 1 | / 2 1 + * / | | 0 1 | * / | | / / | - * * | | *-------* | - * | *-------* | | * - * 3 / / | 3 / - * | 0 / | | 0 + * * |y | *---2---* x | + * | *---3---* | | * + * | / / | | / + * | 0 1 | | / * |/ / | |/ - * *-------* *-------* + * *---2---* x *-------* * @endverbatim *
    - * - * Due to this numbering, the following lines are identical in the - * standard orientation: + * + * The face line numbers (0,1,2,3) correspond to following cell line + * numbers. *
      - *
    • Line 0 of face 0, and line 0 of face 2; - *
    • Line 1 of face 0, and line 3 of face 3; - *
    • Line 2 of face 0, and line 0 of face 4; - *
    • Line 3 of face 0, and line 3 of face 5; - *
    • Line 0 of face 1, and line 2 of face 2; - *
    • Line 1 of face 1, and line 1 of face 3; - *
    • Line 2 of face 1, and line 2 of face 4; - *
    • Line 3 of face 1, and line 1 of face 5; - *
    • Line 3 of face 2, and line 0 of face 5; - *
    • Line 1 of face 2, and line 0 of face 3; - *
    • Line 1 of face 4, and line 2 of face 3; - *
    • Line 3 of face 4, and line 2 of face 5. + *
    • Face 0: lines 8, 10, 0, 4; + *
    • Face 1: lines 9, 11, 1, 5; + *
    • Face 2: lines 2, 6, 8, 9; + *
    • Face 3: lines 3, 7, 10, 11; + *
    • Face 4: lines 0, 1, 2, 3; + *
    • Face 5: lines 4, 5, 6, 7; *
    + * You can get these numbers using the + * GeometryInfo<3>::face_to_cell_lines() function. * - * This standard orientation of faces in 3d can also be depicted by - * assigning a normal vector to each face. The direction of this - * vector (pointing into or out of the cell) is implied by the - * direction of its bounding lines: if you look onto a face and the - * lines are numbered in counter-clockwise sense, then the normal is - * pointing towards you. Thus, faces 1, 2, and 5 have normals that - * point into the cell in standard face orientation, while the - * normals of faces 0, 3, and 4 point outward. Note that opposite - * faces have parallel normal vectors. + * The face normals can be deduced from the face orientation by + * applying the right hand side rule (x,y -> normal). We note, that + * in the standard orientation of faces in 3d, faces 0, 2, and 4 + * have normals that point into the cell, while the normals of faces + * 1, 3, and 5 point outward. * * However, it turns out that a significant number of meshes cannot * satisfy this convention. This is due to the fact that the face @@ -1362,15 +1399,15 @@ struct TriaNumberCache<3> : public TriaNumberCache<2> * *

    Children

    * - * The eight children of a cell are numbered as follows: + * The eight children of a cell are numbered according to the vertices: * @verbatim * *-------* *-------* - * /| 7 6 | / 7 6 /| - * /7| | / /6| - * / | | / 3 2 / | - * * | 4 5 | *-------*2 5| - * |3 4*-------* | 3 2 | * - * | / 4 5 / | | / + * /| 6 7 | / 6 7 /| + * /6| | / /7| + * / | | / 4 5 / | + * * | 2 3 | *-------*5 3| + * |4 2*-------* | 4 5 | * + * | / 2 3 / | | / * |0/ / | |1/ * |/0 1 / | 0 1 |/ * *-------* *-------* @@ -1379,71 +1416,69 @@ struct TriaNumberCache<3> : public TriaNumberCache<2> * Taking into account the orientation of the faces, the following * children are adjacent to the respective faces: *
      - *
    • Face 0: children 0, 1, 2, 3; - *
    • Face 1: children 4, 5, 6, 7; - *
    • Face 2: children 0, 1, 5, 4; - *
    • Face 3: children 1, 5, 6, 2; - *
    • Face 4: children 3, 2, 6, 7; - *
    • Face 5: children 0, 4, 7, 3. + *
    • Face 0: children 0, 2, 4, 6; + *
    • Face 1: children 1, 3, 5, 7; + *
    • Face 2: children 0, 4, 1, 5; + *
    • Face 3: children 2, 6, 3, 7; + *
    • Face 4: children 0, 1, 2, 3; + *
    • Face 5: children 4, 5, 6, 7. *
    - * You can get these numbers using the GeometryInfo<3>@p ::child_cell_on_face - * function. Each child is adjacent to the vertex with the same number. + * You can get these numbers using the + * GeometryInfo<3>::child_cell_on_face() function. As each child is + * adjacent to the vertex with the same number these numbers are + * also given by the GeometryInfo<3>::face_to_cell_vertices() + * function. * * Note that, again, the above list only holds for faces in their * standard orientation. If a face is not in standard orientation, - * then the children at positions 1 and 3 (counting from 0 to 3) - * would be swapped. + * then the children at positions 1 and 2 (counting from 0 to 3) + * would be swapped. In fact, this is what the child_cell_on_face + * and the face_to_cell_vertices functions of GeometryInfo<3> do, + * when invoked with a face_orientation=false argument. * * The information which child cell is at which position of which * face is most often used when computing jump terms across faces * with hanging nodes, using objects of type - * FESubfaceValues. Sitting on one cell, you would look at - * face and figure out which child of the neighbor is sitting on a - * given subface between the present and the neighboring cell. To - * avoid having to query the standard orientation of the faces of - * the two cells every time in such cases, you should use a function - * call like - * cell->neighbor_child_on_subface(face_no,subface_no), which - * returns the correct result both in 2d (where face orientations - * are immaterial) and 3d (where it is necessary to query the face - * orientation and possibly swap the result of - * GeometryInfo<3>::child_cell_on_face). In general, the use of - * GeometryInfo<3>::child_cell_on_face is best avoided due to - * these problems. + * FESubfaceValues. Sitting on one cell, you would look at face and + * figure out which child of the neighbor is sitting on a given + * subface between the present and the neighboring cell. To avoid + * having to query the standard orientation of the faces of the two + * cells every time in such cases, you should use a function call + * like + * cell->neighbor_child_on_subface(face_no,subface_no), + * which returns the correct result both in 2d (where face + * orientations are immaterial) and 3d (where it is necessary to use + * the face orientation as additional argument to + * GeometryInfo<3>::child_cell_on_face). * *

    Coordinate systems

    * * We define the following coordinate system for the explicit coordinates of * the vertices of the unit cell: * @verbatim - * 7-------6 7-------6 + * 6-------7 6-------7 * /| | / /| * / | | / / | * z / | | / / | - * ^ 3 | | 3-------2 | - * | ^y | 4-------5 | | 5 + * ^ 4 | | 4-------5 | + * | ^y | 2-------3 | | 3 * | / | / / | | / * | / | / / | | / * |/ |/ / | |/ * *------>x 0-------1 0-------1 * @endverbatim - * This convention in conjunction with the numbering of the vertices is a bit - * unfortunate, since the vertices 0 through 3 have the coordinates (x,0,z) - * with @p x and @p z being the same as the @p x and @p y coordinates of a quad - * in the plane; more intuitive would have been if they had the coordinates - * (x,y,0). However, the vertex numbering was historically chosen as shown. * * By the convention laid down as above, the vertices have the following - * coordinates: + * coordinates (lexicographic, with x running fastest): *
      *
    • Vertex 0: (0,0,0); *
    • Vertex 1: (1,0,0); - *
    • Vertex 2: (1,0,1); - *
    • Vertex 3: (0,0,1); - *
    • Vertex 4: (0,1,0); - *
    • Vertex 5: (1,1,0); - *
    • Vertex 6: (1,1,1); - *
    • Vertex 7: (0,1,1). + *
    • Vertex 2: (0,1,0); + *
    • Vertex 3: (1,1,0); + *
    • Vertex 4: (0,0,1); + *
    • Vertex 5: (1,0,1); + *
    • Vertex 6: (0,1,1); + *
    • Vertex 7: (1,1,1). *
    * * @@ -1455,7 +1490,7 @@ struct TriaNumberCache<3> : public TriaNumberCache<2> * data stored in the triangulation. * * @ingroup grid - * @author Wolfgang Bangerth, 1998 + * @author Wolfgang Bangerth, 1998; Ralf Hartmann, 2005 */ template class Triangulation : public Subscriptor @@ -1788,7 +1823,22 @@ class Triangulation : public Subscriptor */ virtual void create_triangulation (const std::vector > &vertices, const std::vector > &cells, - const SubCellData &subcelldata); + const SubCellData &subcelldata); + + /** + * For backward compatibility, + * only. This function takes the + * cell data in the ordering as + * requested by deal.II versions + * up to 5.1, converts it to the + * new (lexicographic) ordering + * and calls + * create_triangulation(). + */ + virtual void create_triangulation_compatibility ( + const std::vector > &vertices, + const std::vector > &cells, + const SubCellData &subcelldata); /** * Distort the grid by randomly 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 9adff0ac21..8d54a2a9fa 100644 --- a/deal.II/deal.II/include/grid/tria_accessor.templates.h +++ b/deal.II/deal.II/include/grid/tria_accessor.templates.h @@ -590,20 +590,20 @@ TriaObjectAccessor<3,dim>::line (const unsigned int i) const // oriented correctly, and c) if in // the opposite direction static const unsigned int lookup_table[12][3] = - { { 0, 0, 3 }, // take first four lines from front face - { 0, 1, 2 }, - { 0, 2, 1 }, - { 0, 3, 0 }, + { { 4, 0, 2 }, // take first four lines from bottom face + { 4, 1, 3 }, + { 4, 2, 0 }, + { 4, 3, 1 }, - { 1, 0, 3 }, // second four lines from back face - { 1, 1, 2 }, - { 1, 2, 1 }, - { 1, 3, 0 }, + { 5, 0, 2 }, // second four lines from top face + { 5, 1, 3 }, + { 5, 2, 0 }, + { 5, 3, 1 }, - { 2, 3, 0 }, // the rest randomly - { 2, 1, 2 }, - { 4, 1, 2 }, - { 4, 3, 0 }}; + { 0, 0, 2 }, // the rest randomly + { 1, 0, 2 }, + { 0, 1, 3 }, + { 1, 1, 3 }}; return (this->quad(lookup_table[i][0]) ->line(face_orientation(lookup_table[i][0]) ? @@ -648,20 +648,20 @@ TriaObjectAccessor<3,dim>::line_index (const unsigned int i) const // oriented correctly, and c) if in // the opposite direction static const unsigned int lookup_table[12][3] = - { { 0, 0, 3 }, // take first four lines from front face - { 0, 1, 2 }, - { 0, 2, 1 }, - { 0, 3, 0 }, - - { 1, 0, 3 }, // second four lines from back face - { 1, 1, 2 }, - { 1, 2, 1 }, - { 1, 3, 0 }, - - { 2, 3, 0 }, // the rest randomly - { 2, 1, 2 }, - { 4, 1, 2 }, - { 4, 3, 0 }}; + { { 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 }}; return (this->quad(lookup_table[i][0]) ->line_index(face_orientation(lookup_table[i][0]) ? diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index 432cc8a96e..55c071e6bc 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -1421,16 +1421,16 @@ void DoFTools::make_hanging_node_constraints ( // assert some consistency // assumptions - Assert ((face->child(0)->vertex_index(2) == - face->child(1)->vertex_index(3)) && - (face->child(0)->vertex_index(2) == - face->child(2)->vertex_index(0)) && - (face->child(0)->vertex_index(2) == - face->child(3)->vertex_index(1)), + Assert ((face->child(0)->vertex_index(3) == + face->child(1)->vertex_index(2)) && + (face->child(0)->vertex_index(3) == + face->child(2)->vertex_index(1)) && + (face->child(0)->vertex_index(3) == + face->child(3)->vertex_index(0)), ExcInternalError()); for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) dofs_on_children[next_index++] - = face->child(0)->vertex_dof_index(2,dof); + = face->child(0)->vertex_dof_index(3,dof); // dof numbers on the centers of // the lines bounding this face @@ -1448,13 +1448,13 @@ void DoFTools::make_hanging_node_constraints ( = face->child(0)->line(1)->dof_index(dof); for (unsigned int dof=0; dofchild(1)->line(2)->dof_index(dof); + = face->child(2)->line(1)->dof_index(dof); for (unsigned int dof=0; dofchild(2)->line(3)->dof_index(dof); + = face->child(0)->line(3)->dof_index(dof); for (unsigned int dof=0; dofchild(3)->line(0)->dof_index(dof); + = face->child(1)->line(3)->dof_index(dof); // dofs on the bordering lines for (unsigned int line=0; line<4; ++line) diff --git a/deal.II/deal.II/source/fe/fe_dgq.cc b/deal.II/deal.II/source/fe/fe_dgq.cc index 02bf00a1c2..d88db932ce 100644 --- a/deal.II/deal.II/source/fe/fe_dgq.cc +++ b/deal.II/deal.II/source/fe/fe_dgq.cc @@ -477,13 +477,13 @@ FE_DGQ::has_support_on_face (const unsigned int shape_index, case 2: { - if (face_index==0 && shape_index < n) + if (face_index==0 && (shape_index % n) == 0) return true; if (face_index==1 && (shape_index % n) == this->degree) return true; - if (face_index==2 && shape_index >= this->dofs_per_cell-n) + if (face_index==2 && shape_index < n) return true; - if (face_index==3 && (shape_index % n) == 0) + if (face_index==3 && shape_index >= this->dofs_per_cell-n) return true; return false; }; @@ -492,23 +492,23 @@ FE_DGQ::has_support_on_face (const unsigned int shape_index, { const unsigned int in2 = shape_index % n2; + // x=0 + if (face_index==0 && (shape_index % n) == 0) + return true; + // x=1 + if (face_index==1 && (shape_index % n) == n-1) + return true; // y=0 - if (face_index==0 && in2 < n ) + if (face_index==2 && in2 < n ) return true; // y=1 - if (face_index==1 && in2 >= n2-n) + if (face_index==3 && in2 >= n2-n) return true; // z=0 - if (face_index==2 && shape_index < n2) - return true; - // x=1 - if (face_index==3 && (shape_index % n) == n-1) + if (face_index==4 && shape_index < n2) return true; // z=1 - if (face_index==4 && shape_index >= this->dofs_per_cell - n2) - return true; - // x=0 - if (face_index==5 && (shape_index % n) == 0) + if (face_index==5 && shape_index >= this->dofs_per_cell - n2) return true; return false; }; diff --git a/deal.II/deal.II/source/fe/fe_nedelec.cc b/deal.II/deal.II/source/fe/fe_nedelec.cc index bf1c410234..781252c9a1 100644 --- a/deal.II/deal.II/source/fe/fe_nedelec.cc +++ b/deal.II/deal.II/source/fe/fe_nedelec.cc @@ -135,14 +135,14 @@ FE_Nedelec<2>::shape_value_component (const unsigned int i, { switch (i) { - // (1-y, 0) - case 0: return (component == 0 ? 1-p(1) : 0); + // (0, 1-x) + case 0: return (component == 0 ? 0 : 1-p(0)); // (0,x) case 1: return (component == 0 ? 0 : p(0)); + // (1-y, 0) + case 2: return (component == 0 ? 1-p(1) : 0); // (y, 0) - case 2: return (component == 0 ? p(1) : 0); - // (0, 1-x) - case 3: return (component == 0 ? 0 : 1-p(0)); + case 3: return (component == 0 ? p(1) : 0); // there are only // four shape @@ -190,12 +190,12 @@ FE_Nedelec<3>::shape_value_component (const unsigned int i, // a little. these directions // are: // - // for lines 0, 2, 4, 6: + // for lines 2, 3, 6, 7: // (1,0,0) - // for lines 1, 3, 5, 7: - // (0,0,1) - // for lines 8, 9, 10, 11: + // for lines 0, 1, 4, 5: // (0,1,0) + // for lines 8, 9, 10, 11: + // (0,0,1) // // thus, sort out all those // cases where the component @@ -204,9 +204,9 @@ FE_Nedelec<3>::shape_value_component (const unsigned int i, // spatially dependent part // which is then also the // return value - if (((i<8) && (((i%2==0) && (component!=0)) || - ((i%2==1) && (component!=2)))) || - ((i>=8) && (component != 1))) + if (((i<8) && (((i%4>=2) && (component!=0)) || + ((i%4<2) && (component!=1)))) || + ((i>=8) && (component != 2))) return 0; // now we know that the @@ -218,20 +218,21 @@ FE_Nedelec<3>::shape_value_component (const unsigned int i, z = p(2); switch (i) { - case 0: return (1-y)*(1-z); - case 2: return (1-y)*z; - case 1: return x*(1-y); - case 3: return (1-x)*(1-y); - - case 4: return y*(1-z); - case 6: return y*z; - case 5: return x*y; - case 7: return (1-x)*y; + case 0: return (1-x)*(1-z); + case 1: return x*(1-z); + case 2: return (1-y)*(1-z); + case 3: return y*(1-z); + + case 4: return (1-x)*z; + case 5: return x*z; + case 6: return (1-y)*z; + case 7: return y*z; + + case 8: return (1-x)*(1-y); + case 9: return x*(1-y); + case 10: return (1-x)*y; + case 11: return x*y; - case 8: return (1-x)*(1-z); - case 9: return x*(1-z); - case 10: return x*z; - case 11: return (1-x)*z; default: Assert (false, ExcInternalError()); return 0; @@ -293,10 +294,10 @@ FE_Nedelec<2>::shape_grad_component (const unsigned int i, // third index=component // within gradient static const double unit_gradients[4][2][2] - = { { {0.,-1.}, {0.,0.} }, - { {0.,0.}, {1.,0.} }, - { {0.,+1.}, {0.,0.} }, - { {0.,0.}, {-1.,0.} } }; + = { { {0., 0.}, {-1.,0.} }, + { {0., 0.}, {+1.,0.} }, + { {0.,-1.}, { 0.,0.} }, + { {0.,+1.}, { 0.,0.} } }; return Tensor<1,dim>(unit_gradients[i][component]); }; @@ -352,20 +353,20 @@ FE_Nedelec<3>::shape_grad_component (const unsigned int i, y = p(1), z = p(2); const double unit_gradients[12][3][3] - = { { {0,-(1-z), -(1-y)}, {0,0,0}, { 0, 0, 0} }, - { {0, 0, 0}, {0,0,0}, { (1-y), -x, 0} }, - { {0, -z, (1-y)}, {0,0,0}, { 0, 0, 0} }, - { {0, 0, 0}, {0,0,0}, {-(1-y), -(1-x), 0} }, - - { {0, (1-z), -y}, {0,0,0}, { 0, 0, 0} }, - { {0, 0, 0}, {0,0,0}, { y, x, 0} }, - { {0, z, y}, {0,0,0}, { 0, 0, 0} }, - { {0, 0, 0}, {0,0,0}, { -y, (1-x), 0} }, - - { {0, 0, 0}, {-(1-z), 0, -(1-x)}, {0, 0, 0} }, - { {0, 0, 0}, { (1-z), 0, -x}, {0, 0, 0} }, - { {0, 0, 0}, { z, 0, x}, {0, 0, 0} }, - { {0, 0, 0}, { -z, 0, (1-x)}, {0, 0, 0} } }; + = { { {0, 0, 0}, {-(1-z), 0, -(1-x)}, {0, 0, 0} }, + { {0, 0, 0}, { (1-z), 0, -x}, {0, 0, 0} }, + { {0, -(1-z), -(1-y)}, {0, 0, 0}, {0, 0, 0} }, + { {0, (1-z), -y}, {0, 0, 0}, {0, 0, 0} }, + + { {0, 0, 0}, {-z, 0, (1-x)}, {0, 0, 0} }, + { {0, 0, 0}, { z, 0, x}, {0, 0, 0} }, + { {0, -z, (1-y)}, { 0, 0, 0}, {0, 0, 0} }, + { {0, z, y}, { 0, 0, 0}, {0, 0, 0} }, + + { {0, 0, 0}, {0, 0, 0}, {-(1-y), -(1-x), 0} }, + { {0, 0, 0}, {0, 0, 0}, { (1-y), -x, 0} }, + { {0, 0, 0}, {0, 0, 0}, { -y, (1-x), 0} }, + { {0, 0, 0}, {0, 0, 0}, { y, x, 0} } }; // note: simple check whether // this can at all be: build // the sum over all these @@ -469,29 +470,33 @@ FE_Nedelec<3>::shape_grad_grad_component (const unsigned int i, // second derivative static const double unit_grad_grads[12][3][3][3] = { - { { {0, 0, 0}, {0, 0, 1}, {0, 1, 0} }, - { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} }, + { { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} }, + { {0, 0, 1}, {0, 0, 0}, {1, 0, 0} }, { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} } }, { { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} }, + { {0, 0,-1}, {0, 0, 0}, {-1, 0, 0} }, + { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} } }, + + { { {0, 0, 0}, {0, 0, 1}, {0, 1, 0} }, { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} }, - { {0,-1, 0}, {-1, 0, 0}, {0, 0, 0} } }, + { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} } }, { { {0, 0, 0}, {0, 0,-1}, {0,-1, 0} }, { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} }, { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} } }, { { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} }, - { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} }, - { {0, 1, 0}, { 1, 0, 0}, {0, 0, 0} } }, - - { { {0, 0, 0}, {0, 0,-1}, {0,-1, 0} }, - { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} }, + { {0, 0,-1}, {0, 0, 0}, {-1, 0, 0} }, { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} } }, { { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} }, + { {0, 0, 1}, {0, 0, 0}, {1, 0, 0} }, + { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} } }, + + { { {0, 0, 0}, {0, 0,-1}, {0,-1, 0} }, { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} }, - { {0, 1, 0}, { 1, 0, 0}, {0, 0, 0} } }, + { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} } }, { { {0, 0, 0}, {0, 0, 1}, {0, 1, 0} }, { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} }, @@ -499,23 +504,20 @@ FE_Nedelec<3>::shape_grad_grad_component (const unsigned int i, { { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} }, { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} }, - { {0,-1, 0}, {-1, 0, 0}, {0, 0, 0} } }, + { {0, 1, 0}, {1, 0, 0}, {0, 0, 0} } }, { { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} }, - { {0, 0, 1}, {0, 0, 0}, {1, 0, 0} }, - { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} } }, - - { { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} }, - { {0, 0,-1}, {0, 0, 0}, {-1, 0, 0} }, - { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} } }, + { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} }, + { {0,-1, 0}, {-1, 0, 0}, {0, 0, 0} } }, { { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} }, - { {0, 0, 1}, {0, 0, 0}, {1, 0, 0} }, - { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} } }, + { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} }, + { {0,-1, 0}, {-1, 0, 0}, {0, 0, 0} } }, { { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} }, - { {0, 0,-1}, {0, 0, 0}, {-1, 0, 0} }, - { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} } } }; + { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} }, + { {0, 1, 0}, {1, 0, 0}, {0, 0, 0} } } + }; return Tensor<2,dim>(unit_grad_grads[i][component]); }; @@ -684,11 +686,11 @@ FE_Nedelec::initialize_restriction () for (unsigned int c=0; c::children_per_cell; ++c) this->restriction[c].reinit (this->dofs_per_cell, this->dofs_per_cell); - + this->restriction[0](0,0) = 2.; this->restriction[1](1,1) = 2.; - this->restriction[3](2,2) = 2.; - this->restriction[0](3,3) = 2.; + this->restriction[0](2,2) = 2.; + this->restriction[2](3,3) = 2.; break; }; @@ -726,14 +728,14 @@ FE_Nedelec::initialize_restriction () this->restriction[c].reinit (this->dofs_per_cell, this->dofs_per_cell); this->restriction[0](0,0) = 2.; - this->restriction[0](3,3) = 2.; this->restriction[1](1,1) = 2.; - this->restriction[3](2,2) = 2.; + this->restriction[0](2,2) = 2.; + this->restriction[2](3,3) = 2.; this->restriction[4](4,4) = 2.; - this->restriction[4](7,7) = 2.; this->restriction[5](5,5) = 2.; - this->restriction[7](6,6) = 2.; + this->restriction[4](6,6) = 2.; + this->restriction[6](7,7) = 2.; this->restriction[0](8,8) = 2.; this->restriction[1](9,9) = 2.; @@ -1444,7 +1446,7 @@ FE_Nedelec::has_support_on_face (const unsigned int shape_index, // these in a table const unsigned int opposite_faces[GeometryInfo<2>::faces_per_cell] - = { 2, 3, 0, 1}; + = { 1, 0, 3, 2}; return (face_index != opposite_faces[shape_index]); }; @@ -1460,8 +1462,9 @@ FE_Nedelec::has_support_on_face (const unsigned int shape_index, // defined on const unsigned int opposite_faces[GeometryInfo<3>::lines_per_cell][2] - = { {1,4}, {1,5}, {1,2}, {1,3}, {0,4}, {0,5}, - {0,2}, {0,3}, {3,4}, {4,5}, {2,5}, {2,3}}; + = { {1,5}, {0,5}, {3,5}, {2,5}, + {1,4}, {0,4}, {3,4}, {2,4}, + {1,3}, {0,3}, {1,2}, {0,2}}; return ((face_index != opposite_faces[shape_index][0]) && diff --git a/deal.II/deal.II/source/fe/fe_nedelec_2d.cc b/deal.II/deal.II/source/fe/fe_nedelec_2d.cc index 5fe81b5cf0..9eecf7abe5 100644 --- a/deal.II/deal.II/source/fe/fe_nedelec_2d.cc +++ b/deal.II/deal.II/source/fe/fe_nedelec_2d.cc @@ -35,34 +35,34 @@ namespace FE_Nedelec_2d { static const double q1_into_q1_refined_0[] = { - .5, 0, 0 , 0, - 0, 0.25,0, 0.25, - 0.25, 0, 0.25,0, - 0, 0, 0, .5 + 0.5, 0, 0, 0, + 0.25, 0.25, 0, 0, + 0, 0, 0.5, 0, + 0, 0, 0.25, 0.25 }; static const double q1_into_q1_refined_1[] = { - .5, 0., 0., 0., - 0., .5, 0., 0., - 0.25, 0., 0.25, 0., - 0., 0.25, 0., 0.25, + 0.25, 0.25, 0, 0, + 0, 0.5, 0, 0, + 0, 0, 0.5, 0, + 0, 0, 0.25, 0.25 }; static const double q1_into_q1_refined_2[] = { - 0.25, 0., 0.25, 0., - 0., .5, 0., 0., - 0., 0., .5, 0., - 0., 0.25, 0., 0.25, + 0.5, 0, 0, 0, + 0.25, 0.25, 0, 0, + 0, 0, 0.25, 0.25, + 0, 0, 0, 0.5 }; static const double q1_into_q1_refined_3[] = { - 0.25, 0., 0.25, 0., - 0., 0.25, 0., 0.25, - 0., 0., .5, 0., - 0., 0., 0., .5, + 0.25, 0.25, 0, 0, + 0, 0.5, 0, 0, + 0, 0, 0.25, 0.25, + 0, 0, 0, 0.5 }; } // namespace FE_Nedelec_2d diff --git a/deal.II/deal.II/source/fe/fe_nedelec_3d.cc b/deal.II/deal.II/source/fe/fe_nedelec_3d.cc index 7337e6e58f..ee8d0f9f57 100644 --- a/deal.II/deal.II/source/fe/fe_nedelec_3d.cc +++ b/deal.II/deal.II/source/fe/fe_nedelec_3d.cc @@ -36,133 +36,130 @@ namespace FE_Nedelec_3d { static const double q1_into_q1_refined_0[] = { - .5, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0.25, 0., 0.25, 0., 0., 0., 0.,0.,0.,0.,0., - 0.25, 0., 0.25, 0., 0., 0., 0., 0.,0.,0.,0.,0., - 0., 0., 0., .5, 0., 0., 0., 0.,0.,0.,0.,0., - 0.25, 0., 0., 0., 0.25, 0., 0., 0.,0.,0.,0.,0., - 0., 0.125, 0., 0.125, 0., 0.125, 0., 0.125,0.,0.,0.,0., - 0.125, 0., 0.125, 0., 0.125, 0., 0.125, 0.,0.,0.,0.,0., - 0., 0., 0., 0.25, 0., 0., 0., 0.25,0.,0.,0.,0., - 0., 0., 0., 0., 0., 0., 0., 0., .5, 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0.125, 0.125, 0.125, 0.125, - 0., 0., 0., 0., 0., 0., 0., 0., 0.25, .0, 0., 0.25, + 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0.25, 0.25, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0.25, 0.25, 0, 0, 0, 0, 0, 0, 0, 0, + 0.25, 0, 0, 0, 0.25, 0, 0, 0, 0, 0, 0, 0, + 0.125,0.125,0, 0, 0.125,0.125,0, 0, 0, 0, 0, 0, + 0, 0, 0.25, 0, 0, 0, 0.25, 0, 0, 0, 0, 0, + 0, 0, 0.125,0.125,0, 0, 0.125,0.125,0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0.25, 0.25, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0.25, 0, 0.25, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0.125,0.125,0.125,0.125 }; static const double q1_into_q1_refined_1[] = { - .5, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., .5, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0.25, 0., 0.25, 0., 0., 0., 0., 0.,0.,0.,0.,0., - 0., 0.25, 0., 0.25, 0., 0., 0., 0.,0.,0.,0.,0., - 0.25, 0., 0., 0., 0.25, 0., 0., 0.,0.,0.,0.,0., - 0., 0.25, 0., 0., 0., 0.25, 0., 0.,0.,0.,0.,0., - 0.125, 0., 0.125, 0., 0.125, 0., 0.125, 0.,0.,0.,0.,0., - 0., 0.125, 0., 0.125, 0., 0.125, 0., 0.125,0.,0.,0.,0., - 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., .5, 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0.125, 0.125, 0.125, 0.125, + 0.25, 0.25, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0.25, 0.25, 0, 0, 0, 0, 0, 0, 0, 0, + 0.125,0.125,0, 0, 0.125,0.125,0, 0, 0, 0, 0, 0, + 0, 0.25, 0, 0, 0, 0.25, 0, 0, 0, 0, 0, 0, + 0, 0, 0.25, 0, 0, 0, 0.25, 0, 0, 0, 0, 0, + 0, 0, 0.125,0.125,0, 0, 0.125,0.125,0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0.25, 0.25, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0.125,0.125,0.125,0.125, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.25, 0, 0.25 }; static const double q1_into_q1_refined_2[] = { - 0.25, 0., 0.25, 0., 0., 0., 0., 0.,0.,0.,0.,0., - 0., .5, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., .5, 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0.25, 0., 0.25, 0., 0., 0., 0.,0.,0.,0.,0., - 0.125, 0., 0.125, 0., 0.125, 0., 0.125, 0.,0.,0.,0.,0., - 0., 0.25, 0., 0., 0., 0.25, 0., 0.,0.,0.,0.,0., - 0., 0., 0.25, 0., 0., 0., 0.25, 0.,0.,0.,0.,0., - 0., 0.125, 0., 0.125, 0., 0.125, 0., 0.125,0.,0.,0.,0., - 0., 0., 0., 0., 0., 0., 0., 0., 0.125, 0.125, 0.125, 0.125, - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .5, 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, + 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0.25, 0.25, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0.25, 0.25, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, + 0.25, 0, 0, 0, 0.25, 0, 0, 0, 0, 0, 0, 0, + 0.125,0.125,0, 0, 0.125,0.125,0, 0, 0, 0, 0, 0, + 0, 0, 0.125,0.125,0, 0, 0.125,0.125,0, 0, 0, 0, + 0, 0, 0, 0.25, 0, 0, 0, 0.25, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0.25, 0, 0.25, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0.125,0.125,0.125,0.125, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.25, 0.25 }; static const double q1_into_q1_refined_3[] = { - 0.25, 0., 0.25, 0., 0., 0., 0., 0.,0.,0.,0.,0., - 0., 0.25, 0., 0.25, 0., 0., 0., 0.,0.,0.,0.,0., - 0., 0., .5, 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., .5, 0., 0., 0., 0., 0., 0., 0., 0., - 0.125, 0., 0.125, 0., 0.125, 0., 0.125, 0.,0.,0.,0.,0., - 0., 0.125, 0., 0.125, 0., 0.125, 0., 0.125,0.,0.,0.,0., - 0., 0., 0.25, 0., 0., 0., 0.25, 0.,0.,0.,0.,0., - 0., 0., 0., 0.25, 0., 0., 0., 0.25,0.,0.,0.,0., - 0., 0., 0., 0., 0., 0., 0., 0.,0.25,0.,0.,0.25, - 0., 0., 0., 0., 0., 0., 0., 0., 0.125, 0.125, 0.125, 0.125, - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .5, + 0.25, 0.25, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0.25, 0.25, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, + 0.125,0.125,0, 0, 0.125,0.125,0, 0, 0, 0, 0, 0, + 0, 0.25, 0, 0, 0, 0.25, 0, 0, 0, 0, 0, 0, + 0, 0, 0.125,0.125,0, 0, 0.125,0.125,0, 0, 0, 0, + 0, 0, 0, 0.25, 0, 0, 0, 0.25, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0.125,0.125,0.125,0.125, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.25, 0, 0.25, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.25, 0.25, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5 }; static const double q1_into_q1_refined_4[] = { - 0.25, 0., 0., 0., 0.25, 0., 0., 0.,0.,0.,0.,0., - 0., 0.125, 0., 0.125, 0., 0.125, 0., 0.125,0.,0.,0.,0., - 0.125, 0., 0.125, 0., 0.125, 0., 0.125, 0.,0.,0.,0.,0., - 0., 0., 0., 0.25, 0., 0., 0., 0.25, 0., 0., 0., 0., - 0., 0., 0., 0., .5, 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0.25, 0., 0.25,0.,0.,0.,0., - 0., 0., 0., 0., 0.25, 0., 0.25, 0.,0.,0.,0.,0., - 0., 0., 0., 0., 0., 0., 0., .5, 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., .5, 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0.125, 0.125, 0.125, 0.125, - 0., 0., 0., 0., 0., 0., 0., 0.,0.25,0.,0.,0.25, + 0.25, 0, 0, 0, 0.25, 0, 0, 0, 0, 0, 0, 0, + 0.125,0.125,0, 0, 0.125,0.125,0, 0, 0, 0, 0, 0, + 0, 0, 0.25, 0, 0, 0, 0.25, 0, 0, 0, 0, 0, + 0, 0, 0.125,0.125,0, 0, 0.125,0.125,0, 0, 0, 0, + 0, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0.25, 0.25, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0.25, 0.25, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0.25, 0.25, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0.25, 0, 0.25, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0.125,0.125,0.125,0.125 }; static const double q1_into_q1_refined_5[] = - { - 0.25, 0., 0., 0., 0.25, 0., 0., 0.,0.,0.,0.,0., - 0., 0.25, 0., 0., 0., 0.25, 0., 0.,0.,0.,0.,0., - 0.125, 0., 0.125, 0., 0.125, 0., 0.125, 0.,0.,0.,0.,0., - 0., 0.125, 0., 0.125, 0., 0.125, 0., 0.125,0.,0.,0.,0., - 0., 0., 0., 0., .5, 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., .5, 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0.25, 0., 0.25, 0.,0.,0.,0.,0., - 0., 0., 0., 0., 0., 0.25, 0., 0.25,0.,0.,0.,0., - 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., .5, 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0.125, 0.125, 0.125, 0.125, - + { + 0.125,0.125,0, 0, 0.125,0.125,0, 0, 0, 0, 0, 0, + 0, 0.25, 0, 0, 0, 0.25, 0, 0, 0, 0, 0, 0, + 0, 0, 0.25, 0, 0, 0, 0.25, 0, 0, 0, 0, 0, + 0, 0, 0.125,0.125,0, 0, 0.125,0.125,0, 0, 0, 0, + 0, 0, 0, 0, 0.25, 0.25, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0.25, 0.25, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0.25, 0.25, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0.125,0.125,0.125,0.125, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.25, 0, 0.25 }; static const double q1_into_q1_refined_6[] = { - 0.125, 0., 0.125, 0., 0.125, 0., 0.125, 0.,0.,0.,0.,0., - 0., 0.25, 0., 0., 0., 0.25, 0., 0.,0.,0.,0.,0., - 0., 0., 0.25, 0., 0., 0., 0.25, 0.,0.,0.,0.,0., - 0., 0.125, 0., 0.125, 0., 0.125, 0., 0.125,0.,0.,0.,0., - 0., 0., 0., 0., 0.25, 0., 0.25, 0.,0.,0.,0.,0., - 0., 0., 0., 0., 0., .5, 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., .5, 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0.25, 0., 0.25, 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0.125, 0.125, 0.125, 0.125, - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .5, 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, - - + 0.25, 0, 0, 0, 0.25, 0, 0, 0, 0, 0, 0, 0, + 0.125,0.125,0, 0, 0.125,0.125,0, 0, 0, 0, 0, 0, + 0, 0, 0.125,0.125,0, 0, 0.125,0.125,0, 0, 0, 0, + 0, 0, 0, 0.25, 0, 0, 0, 0.25, 0, 0, 0, 0, + 0, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0.25, 0.25, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0.25, 0.25, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0.25, 0, 0.25, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0.125,0.125,0.125,0.125, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.25, 0.25 }; static const double q1_into_q1_refined_7[] = { - 0.125, 0., 0.125, 0., 0.125, 0., 0.125, 0.,0.,0.,0.,0., - 0., 0.125, 0., 0.125, 0., 0.125, 0., 0.125,0.,0.,0.,0., - 0., 0., 0.25, 0., 0., 0., 0.25, 0.,0.,0.,0.,0., - 0., 0., 0., 0.25, 0., 0., 0., 0.25,0.,0.,0.,0., - 0., 0., 0., 0., 0.25, 0., 0.25, 0.,0.,0.,0.,0., - 0., 0., 0., 0., 0., 0.25, 0., 0.25,0.,0.,0.,0., - 0., 0., 0., 0., 0., 0., .5, 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., .5, 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0., 0., 0.25, - 0., 0., 0., 0., 0., 0., 0., 0., 0.125, 0.125, 0.125, 0.125, - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5, + 0.125,0.125,0, 0, 0.125,0.125,0, 0, 0, 0, 0, 0, + 0, 0.25, 0, 0, 0, 0.25, 0, 0, 0, 0, 0, 0, + 0, 0, 0.125,0.125,0, 0, 0.125,0.125,0, 0, 0, 0, + 0, 0, 0, 0.25, 0, 0, 0, 0.25, 0, 0, 0, 0, + 0, 0, 0, 0, 0.25, 0.25, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0.25, 0.25, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0.125,0.125,0.125,0.125, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.25, 0, 0.25, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.25, 0.25, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5 }; } // namespace FE_Nedelec_3d @@ -201,10 +198,10 @@ namespace FE_Nedelec_3d { static const double constraint_q1[] = { - 0, .25, 0, .25, // first the four interior lines - .25, 0, .25, 0, - 0, .25, 0, .25, - .25, 0, .25, 0, + .25, .25, 0, 0, // first the four interior lines + .25, .25, 0, 0, + 0, 0, .25, .25, + 0, 0, .25, .25, .5, 0, 0, 0, // then the two child lines of each of the four outer .5, 0, 0, 0, // ones. since the shape functions are constant on each 0, .5, 0, 0, // line, the two child lines get the same weights, modulo diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index 47abd0e0db..4febfb4e0d 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -613,15 +613,15 @@ FE_Q<3>::initialize_constraints () // described in the documentation // of the FiniteElement class // (fe_base.h), i.e. - // *--13--3--14--* + // *--15--4--16--* // | | | - // 16 20 7 19 12 + // 10 19 6 20 12 // | | | - // 4--8---0--6---2 + // 1--7---0--8---2 // | | | - // 15 17 5 18 11 + // 9 17 5 18 11 // | | | - // *--9---1--10--* + // *--13--3--14--* std::vector > constraint_points; // Add midpoint @@ -629,10 +629,10 @@ FE_Q<3>::initialize_constraints () // Add midpoints of lines of // "mother-face" - constraint_points.push_back (Point (0.5, 0)); + constraint_points.push_back (Point (0, 0.5)); constraint_points.push_back (Point (1, 0.5)); + constraint_points.push_back (Point (0.5, 0)); constraint_points.push_back (Point (0.5, 1)); - constraint_points.push_back (Point (0, 0.5)); if (this->degree>1) { @@ -649,22 +649,22 @@ FE_Q<3>::initialize_constraints () // Add nodes of lines interior // in the "mother-face" - // line 5: use line 15 - QProjector::project_to_subface(qline, 3, 0, p_line); + // line 5: use line 9 + QProjector::project_to_subface(qline, 0, 0, p_line); for (unsigned int i=0; i (0.5, 0)); // line 6: use line 10 QProjector::project_to_subface(qline, 0, 1, p_line); - for (unsigned int i=0; i (0, 0.5)); - // line 7: use line 16 - QProjector::project_to_subface(qline, 3, 1, p_line); for (unsigned int i=0; i (0.5, 0)); - // line 8: use line 9 - QProjector::project_to_subface(qline, 0, 0, p_line); + // line 7: use line 13 + QProjector::project_to_subface(qline, 2, 0, p_line); for (unsigned int i=0; i (0, 0.5)); + constraint_points.push_back (p_line[i] + Point (0, 0.5)); + // line 8: use line 14 + QProjector::project_to_subface(qline, 2, 1, p_line); + for (unsigned int i=0; i (0, 0.5)); // DoFs on bordering lines // lines 9-16 diff --git a/deal.II/deal.II/source/fe/fe_q_hierarchical.cc b/deal.II/deal.II/source/fe/fe_q_hierarchical.cc index 0c0a459d90..4e998dcaff 100644 --- a/deal.II/deal.II/source/fe/fe_q_hierarchical.cc +++ b/deal.II/deal.II/source/fe/fe_q_hierarchical.cc @@ -235,17 +235,17 @@ initialize_constraints (const std::vector > &dofs_subcell) // boundary vertex nodes this->interface_constraints(1,face_renumber[i]) = - dofs_subcell[0](1, i % dofs_1d) * - dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d); + dofs_subcell[0](0, i % dofs_1d) * + dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d); this->interface_constraints(2,face_renumber[i]) = dofs_subcell[1](1, i % dofs_1d) * dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d); this->interface_constraints(3,face_renumber[i]) = + dofs_subcell[0](1, i % dofs_1d) * + dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d); + this->interface_constraints(4,face_renumber[i]) = dofs_subcell[1](0, i % dofs_1d) * dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d); - this->interface_constraints(4,face_renumber[i]) = - dofs_subcell[0](0, i % dofs_1d) * - dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d); // interior edges for (unsigned int j=0; j<(degree-1); j++) @@ -254,26 +254,26 @@ initialize_constraints (const std::vector > &dofs_subcell) dofs_subcell[0](1, i % dofs_1d) * dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d); this->interface_constraints(5 + (degree-1) + j,face_renumber[i]) = - dofs_subcell[1](2 + j, i % dofs_1d) * - dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d); - this->interface_constraints(5 + 2*(degree-1) + j,face_renumber[i]) = dofs_subcell[0](1,i % dofs_1d) * dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d); - this->interface_constraints(5 + 3*(degree-1) + j,face_renumber[i]) = + this->interface_constraints(5 + 2*(degree-1) + j,face_renumber[i]) = dofs_subcell[0](2 + j,i % dofs_1d) * dofs_subcell[1](0, (i - (i % dofs_1d)) / dofs_1d); + this->interface_constraints(5 + 3*(degree-1) + j,face_renumber[i]) = + dofs_subcell[1](2 + j, i % dofs_1d) * + dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d); } // boundary edges for (unsigned int j=0; j<(degree-1); j++) { - // bottom edge + // left edge this->interface_constraints(5 + 4*(degree-1) + j,face_renumber[i]) = - dofs_subcell[0](2 + j, i % dofs_1d) * - dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d); + dofs_subcell[0](0, i % dofs_1d) * + dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d); this->interface_constraints(5 + 4*(degree-1) + (degree-1) + j,face_renumber[i]) = - dofs_subcell[1](2 + j, i % dofs_1d) * - dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d); + dofs_subcell[0](0, i % dofs_1d) * + dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d); // right edge this->interface_constraints(5 + 4*(degree-1) + 2*(degree-1) + j,face_renumber[i]) = dofs_subcell[1](1, i % dofs_1d) * @@ -281,20 +281,20 @@ initialize_constraints (const std::vector > &dofs_subcell) this->interface_constraints(5 + 4*(degree-1) + 3*(degree-1) + j,face_renumber[i]) = dofs_subcell[1](1, i % dofs_1d) * dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d); - // top edge + // bottom edge this->interface_constraints(5 + 4*(degree-1) + 4*(degree-1) + j,face_renumber[i]) = dofs_subcell[0](2 + j, i % dofs_1d) * - dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d); + dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d); this->interface_constraints(5 + 4*(degree-1) + 5*(degree-1) + j,face_renumber[i]) = dofs_subcell[1](2 + j, i % dofs_1d) * - dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d); - // left edge + dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d); + // top edge this->interface_constraints(5 + 4*(degree-1) + 6*(degree-1) + j,face_renumber[i]) = - dofs_subcell[0](0, i % dofs_1d) * - dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d); + dofs_subcell[0](2 + j, i % dofs_1d) * + dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d); this->interface_constraints(5 + 4*(degree-1) + 7*(degree-1) + j,face_renumber[i]) = - dofs_subcell[0](0, i % dofs_1d) * - dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d); + dofs_subcell[1](2 + j, i % dofs_1d) * + dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d); } // interior faces @@ -311,11 +311,11 @@ initialize_constraints (const std::vector > &dofs_subcell) dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d); // subcell 2 this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1) + 2*(degree-1)*(degree-1),face_renumber[i]) = - dofs_subcell[1](2 + j, i % dofs_1d) * + dofs_subcell[0](2 + j, i % dofs_1d) * dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d); // subcell 3 this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1) + 3*(degree-1)*(degree-1),face_renumber[i]) = - dofs_subcell[0](2 + j, i % dofs_1d) * + dofs_subcell[1](2 + j, i % dofs_1d) * dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d); } } @@ -375,8 +375,10 @@ initialize_embedding_and_restriction (const std::vector > &do { for (unsigned int c=0; c::children_per_cell; ++c) { - unsigned int c0 = ((c==1) || (c==2)) ? 1 : 0; - unsigned int c1 = ((c==2) || (c==3)) ? 1 : 0; + // left/right line: 0/1 + const unsigned int c0 = c%2; + // bottom/top line: 0/1 + const unsigned int c1 = c/2; this->prolongation[c](j,i) = dofs_subcell[c0](renumber[j] % dofs_1d, @@ -397,9 +399,12 @@ initialize_embedding_and_restriction (const std::vector > &do { for (unsigned int c=0; c::children_per_cell; ++c) { - unsigned int c0 = ((c==1) || (c==2) || (c==5) || (c==6)) ? 1 : 0; - unsigned int c1 = ((c==4) || (c==5) || (c==6) || (c==7)) ? 1 : 0; - unsigned int c2 = ((c==2) || (c==3) || (c==6) || (c==7)) ? 1 : 0; + // left/right face: 0/1 + const unsigned int c0 = c%2; + // front/back face: 0/1 + const unsigned int c1 = (c%4)/2; + // bottom/top face: 0/1 + const unsigned int c2 = c/4; this->prolongation[c](j,i) = dofs_subcell[c0](renumber[j] % dofs_1d, @@ -612,10 +617,10 @@ hierarchic_to_fe_q_hierarchical_numbering (const FiniteElementData &fe) // Example: degree=3 // // hierarchical numbering: - // 3 8 9 2 - // 11 14 15 7 - // 10 12 13 6 - // 0 4 5 1 + // 2 10 11 3 + // 5 14 15 7 + // 4 12 13 6 + // 0 8 9 1 // // fe_q_hierarchical numbering: // 4 6 7 5 @@ -626,20 +631,20 @@ hierarchic_to_fe_q_hierarchical_numbering (const FiniteElementData &fe) // first the four vertices h2l[next_index++] = 0; h2l[next_index++] = 1; - h2l[next_index++] = n+1; h2l[next_index++] = n; - // bottom line + h2l[next_index++] = n+1; + // left line for (unsigned int i=0; i &fe) unsigned int next_index = 0; const unsigned int n2=n*n; // first the eight vertices - // front face, counterclock wise + // bottom face, lexicographic h2l[next_index++] = 0; h2l[next_index++] = 1; - h2l[next_index++] = n2+1; - h2l[next_index++] = n2; - // back face, counterclock wise h2l[next_index++] = n; h2l[next_index++] = n+1; - h2l[next_index++] = n2+n+1; + // top face, lexicographic + h2l[next_index++] = n2; + h2l[next_index++] = n2+1; h2l[next_index++] = n2+n; + h2l[next_index++] = n2+n+1; // now the lines - // front face, counterclock wise - for (unsigned int i=0; i &fe) for (unsigned int i=0; i::face_fe_q_hierarchical_to_hierarchic_numbering (const unsi template <> bool FE_Q_Hierarchical<1>::has_support_on_face (const unsigned int shape_index, - const unsigned int face_index) const + const unsigned int face_index) const { Assert (shape_index < this->dofs_per_cell, - ExcIndexRange (shape_index, 0, this->dofs_per_cell)); + ExcIndexRange (shape_index, 0, this->dofs_per_cell)); Assert (face_index < GeometryInfo<1>::faces_per_cell, - ExcIndexRange (face_index, 0, GeometryInfo<1>::faces_per_cell)); + ExcIndexRange (face_index, 0, GeometryInfo<1>::faces_per_cell)); - // in 1d, things are simple. since - // there is only one degree of - // freedom per vertex in this - // class, the first is on vertex 0 - // (==face 0 in some sense), the - // second on face 1: + // in 1d, things are simple. since + // there is only one degree of + // freedom per vertex in this + // class, the first is on vertex 0 + // (==face 0 in some sense), the + // second on face 1: return (((shape_index == 0) && (face_index == 0)) || - ((shape_index == 1) && (face_index == 1))); + ((shape_index == 1) && (face_index == 1))); } #endif + template bool FE_Q_Hierarchical::has_support_on_face (const unsigned int shape_index, @@ -801,10 +806,10 @@ FE_Q_Hierarchical::has_support_on_face (const unsigned int shape_index, Assert (face_index < GeometryInfo::faces_per_cell, ExcIndexRange (face_index, 0, GeometryInfo::faces_per_cell)); - // first, special-case interior - // shape functions, since they - // have no support no-where on - // the boundary + // first, special-case interior + // shape functions, since they + // have no support no-where on + // the boundary if (((dim==2) && (shape_index>=this->first_quad_index)) || ((dim==3) && (shape_index>=this->first_hex_index))) @@ -814,11 +819,11 @@ FE_Q_Hierarchical::has_support_on_face (const unsigned int shape_index, // vertex if (shape_index < this->first_line_index) { - // for Q elements, there is one - // dof per vertex, so + // for Q elements, there is + // one dof per vertex, so // shape_index==vertex_number. check - // whether this vertex is on - // the given face + // whether this vertex is + // on the given face. const unsigned int vertex_no = shape_index; Assert (vertex_no < GeometryInfo::vertices_per_cell, ExcInternalError()); @@ -849,10 +854,10 @@ FE_Q_Hierarchical::has_support_on_face (const unsigned int shape_index, static_cast(GeometryInfo::quads_per_cell), ExcInternalError()); - // in 2d, cell bubble are zero - // on all faces. but we have - // treated this case above - // already + // in 2d, cell bubble are + // zero on all faces. but + // we have treated this + // case above already Assert (dim != 2, ExcInternalError()); // in 3d, @@ -860,7 +865,7 @@ FE_Q_Hierarchical::has_support_on_face (const unsigned int shape_index, if (dim == 3) return (quad_index == face_index); else - Assert (false, ExcNotImplemented()); + Assert (false, ExcNotImplemented()); } else // dof on hex diff --git a/deal.II/deal.II/source/fe/fe_tools.all_dimensions.cc b/deal.II/deal.II/source/fe/fe_tools.all_dimensions.cc index 501198b33a..c46d112199 100644 --- a/deal.II/deal.II/source/fe/fe_tools.all_dimensions.cc +++ b/deal.II/deal.II/source/fe/fe_tools.all_dimensions.cc @@ -63,20 +63,25 @@ FETools::hierarchic_to_lexicographic_numbering (const FiniteElementData &fe // first the four vertices h2l[next_index++] = 0; h2l[next_index++] = n-1; - h2l[next_index++] = n*n-1; h2l[next_index++] = n*(n-1); - // first line + h2l[next_index++] = n*n-1; + + // left line for (unsigned int i=0; i &fe { unsigned int next_index = 0; // first the eight vertices - h2l[next_index++] = 0; - h2l[next_index++] = n-1; - h2l[next_index++] = (n-1)*(n*n+1); - h2l[next_index++] = (n-1)*n*n; - h2l[next_index++] = n*(n-1); - h2l[next_index++] = n*n-1; - h2l[next_index++] = n*n*n-1; - h2l[next_index++] = (n-1)*(n*n+n); + h2l[next_index++] = 0; // 0 + h2l[next_index++] = ( 1)*degree; // 1 + h2l[next_index++] = ( n )*degree; // 2 + h2l[next_index++] = ( n+1)*degree; // 3 + h2l[next_index++] = (n*n )*degree; // 4 + h2l[next_index++] = (n*n +1)*degree; // 5 + h2l[next_index++] = (n*n+n )*degree; // 6 + h2l[next_index++] = (n*n+n+1)*degree; // 7 + // line 0 for (unsigned int i=0; i &fe - template void FETools::lexicographic_to_hierarchic_numbering (const FiniteElementData &fe, diff --git a/deal.II/deal.II/source/fe/mapping.cc b/deal.II/deal.II/source/fe/mapping.cc index ab61598620..9a8b8d4d11 100644 --- a/deal.II/deal.II/source/fe/mapping.cc +++ b/deal.II/deal.II/source/fe/mapping.cc @@ -15,37 +15,6 @@ #include -#if (deal_II_dimension == 1) - -template<> -const unsigned int Mapping<1>::normal_directions[2] = -{ - 1, 0 -}; - -#endif - -#if (deal_II_dimension == 2) - -template<> -const unsigned int Mapping<2>::normal_directions[4] = -{ - 2, 0, 3, 1 -}; - -#endif - -#if (deal_II_dimension == 3) - -template<> -const unsigned int Mapping<3>::normal_directions[6] = -{ - 3, 2, 5, 0, 4, 1 -}; - -#endif - - template Mapping::~Mapping () {} diff --git a/deal.II/deal.II/source/fe/mapping_cartesian.cc b/deal.II/deal.II/source/fe/mapping_cartesian.cc index 5a5a2150a2..06f734d7ba 100644 --- a/deal.II/deal.II/source/fe/mapping_cartesian.cc +++ b/deal.II/deal.II/source/fe/mapping_cartesian.cc @@ -176,12 +176,12 @@ MappingCartesian::compute_fill (const typename Triangulation::cell_ite break; case 2: data.length[0] = cell->vertex(1)(0) - start(0); - data.length[1] = cell->vertex(3)(1) - start(1); + data.length[1] = cell->vertex(2)(1) - start(1); break; case 3: data.length[0] = cell->vertex(1)(0) - start(0); - data.length[1] = cell->vertex(4)(1) - start(1); - data.length[2] = cell->vertex(3)(2) - start(2); + data.length[1] = cell->vertex(2)(1) - start(1); + data.length[2] = cell->vertex(4)(2) - start(2); break; default: Assert(false, ExcNotImplemented()); @@ -240,10 +240,10 @@ MappingCartesian::compute_fill (const typename Triangulation::cell_ite { static const Point normals[GeometryInfo<2>::faces_per_cell] - = { Point(0, -1), - Point(1, 0), - Point(0, 1), - Point(-1, 0) }; + = { Point(-1, 0), + Point( 1, 0), + Point( 0,-1), + Point( 0, 1) }; std::fill (normal_vectors.begin(), normal_vectors.end(), normals[face_no]); @@ -254,12 +254,12 @@ MappingCartesian::compute_fill (const typename Triangulation::cell_ite { static const Point normals[GeometryInfo<3>::faces_per_cell] - = { Point(0, -1, 0), - Point(0, 1, 0), - Point(0, 0, -1), - Point(1, 0, 0), - Point(0, 0, 1), - Point(-1, 0, 0) }; + = { Point(-1, 0, 0), + Point( 1, 0, 0), + Point( 0,-1, 0), + Point( 0, 1, 0), + Point( 0, 0,-1), + Point( 0, 0, 1) }; std::fill (normal_vectors.begin(), normal_vectors.end(), normals[face_no]); @@ -340,8 +340,8 @@ MappingCartesian::fill_fe_face_values (const typename Triangulation::c // product of the local lengths // since the jacobian is diagonal double J = 1.; - for (unsigned int d=0;dnormal_directions[face_no]/2)) + for (unsigned int d=0; d::unit_normal_direction[face_no]) J *= data.length[d]; if (data.current_update_flags() & update_JxW_values) @@ -383,8 +383,8 @@ MappingCartesian::fill_fe_subface_values (const typename Triangulation // product of the local lengths // since the jacobian is diagonal double J = 1.; - for (unsigned int d=0;dnormal_directions[face_no]/2)) + for (unsigned int d=0; d::unit_normal_direction[face_no]) J *= data.length[d]; if (data.current_update_flags() & update_JxW_values) @@ -552,12 +552,12 @@ MappingCartesian::transform_unit_to_real_cell ( break; case 2: length[0] = cell->vertex(1)(0) - start(0); - length[1] = cell->vertex(3)(1) - start(1); + length[1] = cell->vertex(2)(1) - start(1); break; case 3: length[0] = cell->vertex(1)(0) - start(0); - length[1] = cell->vertex(4)(1) - start(1); - length[2] = cell->vertex(3)(2) - start(2); + length[1] = cell->vertex(2)(1) - start(1); + length[2] = cell->vertex(4)(2) - start(2); break; default: Assert(false, ExcNotImplemented()); @@ -589,12 +589,12 @@ MappingCartesian::transform_real_to_unit_cell ( break; case 2: real(0) /= cell->vertex(1)(0) - start(0); - real(1) /= cell->vertex(3)(1) - start(1); + real(1) /= cell->vertex(2)(1) - start(1); break; case 3: real(0) /= cell->vertex(1)(0) - start(0); - real(1) /= cell->vertex(4)(1) - start(1); - real(2) /= cell->vertex(3)(2) - start(2); + real(1) /= cell->vertex(2)(1) - start(1); + real(2) /= cell->vertex(4)(2) - start(2); break; default: Assert(false, ExcNotImplemented()); diff --git a/deal.II/deal.II/source/fe/mapping_q.cc b/deal.II/deal.II/source/fe/mapping_q.cc index b8b8205a33..5e53b1efc3 100644 --- a/deal.II/deal.II/source/fe/mapping_q.cc +++ b/deal.II/deal.II/source/fe/mapping_q.cc @@ -459,14 +459,14 @@ MappingQ::set_laplace_on_quad_vector(Table<2,double> &loqvs) const { // (same as above) static const double loqv3[4*12] - ={80/1053., 1/81., 11/1053., 1/81., 25/117., 44/351., - 7/117., 16/351., 7/117., 16/351., 25/117., 44/351., - 1/81., 80/1053., 1/81., 11/1053., 44/351., 25/117., - 25/117., 44/351., 16/351., 7/117., 7/117., 16/351., - 1/81., 11/1053., 1/81., 80/1053., 7/117., 16/351., - 16/351., 7/117., 25/117., 44/351., 44/351., 25/117., - 11/1053., 1/81., 80/1053., 1/81., 16/351., 7/117., - 44/351., 25/117., 44/351., 25/117., 16/351., 7/117.}; + ={80/1053., 1/81., 1/81., 11/1053., 25/117., 44/351., + 7/117., 16/351., 25/117., 44/351., 7/117., 16/351., + 1/81., 80/1053., 11/1053., 1/81., 7/117., 16/351., + 25/117., 44/351., 44/351., 25/117., 16/351., 7/117., + 1/81., 11/1053., 80/1053., 1/81., 44/351., 25/117., + 16/351., 7/117., 7/117., 16/351., 25/117., 44/351., + 11/1053., 1/81., 1/81., 80/1053., 16/351., 7/117., + 44/351., 25/117., 16/351., 7/117., 44/351., 25/117.}; Assert (sizeof(loqv3)/sizeof(loqv3[0]) == n_inner_2d * n_outer_2d, ExcInternalError()); @@ -480,86 +480,86 @@ MappingQ::set_laplace_on_quad_vector(Table<2,double> &loqvs) const case 4: { static const double loqv4[9*16] - ={0.074059218503115712, -0.0010757446289059922, - 0.0019142922390714627, -0.0010757446289060069, - 0.22312738654318912, 0.13468513060151868, - 0.03812914216116723, 0.029131600026332517, - 0.02200737428129396, 0.016008355644312244, - 0.029131600026332527, 0.022007374281293915, - 0.016008355644312224, 0.22312738654318917, - 0.13468513060151874, 0.038129142161167237, - - 0.0066480315133420603, 0.0066480315133427186, - 0.0028734528616576258, 0.0028734528616574584, - 0.15277168188202372, 0.23481527607092728, - 0.15277168188202397, 0.079035726825841868, - 0.059692382812499903, 0.036198648174158153, - 0.024962693117977778, 0.040819489554071289, - 0.024962693117977889, 0.079035726825843783, - 0.059692382812500312, 0.036198648174158236, - - -0.0010757446289069229, 0.074059218503115892, - -0.0010757446289058842, 0.0019142922390713395, - 0.038129142161167293, 0.13468513060151846, - 0.22312738654318981, 0.22312738654318776, - 0.13468513060151835, 0.038129142161167202, - 0.016008355644312171, 0.022007374281293943, - 0.029131600026332617, 0.029131600026335094, - 0.02200737428129395, 0.016008355644312293, - - 0.0066480315133420733, 0.0028734528616574727, - 0.0028734528616576362, 0.0066480315133427255, - 0.079035726825843755, 0.059692382812500257, - 0.036198648174158243, 0.024962693117977792, - 0.04081948955407131, 0.024962693117977903, - 0.079035726825841868, 0.059692382812499799, - 0.036198648174158098, 0.15277168188202384, - 0.23481527607092734, 0.15277168188202397, - - 0.011067708333333018, 0.011067708333333358, - 0.011067708333333736, 0.01106770833333337, - 0.067708333333334217, 0.1035156250000009, - 0.067708333333334356, 0.067708333333333759, - 0.10351562499999903, 0.067708333333333995, - 0.067708333333333703, 0.10351562499999885, - 0.067708333333333898, 0.067708333333334245, - 0.10351562500000108, 0.067708333333334397, - - 0.0028734528616571847, 0.0066480315133423621, - 0.0066480315133430378, 0.0028734528616573343, - 0.036198648174158195, 0.059692382812500278, - 0.079035726825844074, 0.15277168188202336, - 0.23481527607092598, 0.15277168188202367, - 0.036198648174158042, 0.059692382812499861, - 0.079035726825842201, 0.024962693117977788, - 0.040819489554074009, 0.024962693117977879, - - -0.0010757446289069133, 0.0019142922390713401, - -0.0010757446289058651, 0.07405921850311592, - 0.029131600026335087, 0.022007374281293911, - 0.016008355644312279, 0.016008355644312175, - 0.022007374281293957, 0.029131600026332641, - 0.22312738654318776, 0.13468513060151827, - 0.038129142161167182, 0.038129142161167286, - 0.13468513060151852, 0.22312738654318992, - - 0.0028734528616571756, 0.0028734528616573209, - 0.0066480315133430369, 0.0066480315133423742, - 0.024962693117977757, 0.040819489554073919, - 0.024962693117977847, 0.036198648174158049, - 0.059692382812499924, 0.079035726825842215, - 0.15277168188202328, 0.23481527607092575, - 0.15277168188202356, 0.036198648174158167, - 0.059692382812500319, 0.079035726825844088, - - 0.0019142922390712369, -0.0010757446289068027, - 0.07405921850311617, -0.0010757446289067784, - 0.016008355644312283, 0.022007374281293967, - 0.02913160002633523, 0.038129142161167258, - 0.13468513060151821, 0.22312738654318864, - 0.038129142161167265, 0.13468513060151813, - 0.22312738654318859, 0.016008355644312276, - 0.022007374281294009, 0.029131600026335244 + ={0.07405921850311571, -0.001075744628905992, + -0.001075744628906007, 0.001914292239071463, + 0.2231273865431892, 0.1346851306015187, + 0.03812914216116724, 0.02913160002633252, + 0.02200737428129396, 0.01600835564431224, + 0.2231273865431891, 0.1346851306015187, + 0.03812914216116723, 0.02913160002633253, + 0.02200737428129391, 0.01600835564431222, + + 0.00664803151334206, 0.006648031513342719, + 0.002873452861657458, 0.002873452861657626, + 0.07903572682584378, 0.05969238281250031, + 0.03619864817415824, 0.07903572682584187, + 0.0596923828124999, 0.03619864817415815, + 0.1527716818820237, 0.2348152760709273, + 0.152771681882024, 0.02496269311797778, + 0.04081948955407129, 0.02496269311797789, + + -0.001075744628906923, 0.07405921850311589, + 0.001914292239071339, -0.001075744628905884, + 0.02913160002633509, 0.02200737428129395, + 0.01600835564431229, 0.2231273865431878, + 0.1346851306015183, 0.0381291421611672, + 0.03812914216116729, 0.1346851306015185, + 0.2231273865431898, 0.01600835564431217, + 0.02200737428129394, 0.02913160002633262, + + 0.006648031513342073, 0.002873452861657473, + 0.006648031513342726, 0.002873452861657636, + 0.1527716818820238, 0.2348152760709273, + 0.152771681882024, 0.02496269311797779, + 0.04081948955407131, 0.0249626931179779, + 0.07903572682584376, 0.05969238281250026, + 0.03619864817415824, 0.07903572682584187, + 0.0596923828124998, 0.0361986481741581, + + 0.01106770833333302, 0.01106770833333336, + 0.01106770833333337, 0.01106770833333374, + 0.06770833333333424, 0.1035156250000011, + 0.0677083333333344, 0.06770833333333376, + 0.103515624999999, 0.06770833333333399, + 0.06770833333333422, 0.1035156250000009, + 0.06770833333333436, 0.0677083333333337, + 0.1035156249999988, 0.0677083333333339, + + 0.002873452861657185, 0.006648031513342362, + 0.002873452861657334, 0.006648031513343038, + 0.02496269311797779, 0.04081948955407401, + 0.02496269311797788, 0.1527716818820234, + 0.234815276070926, 0.1527716818820237, + 0.03619864817415819, 0.05969238281250028, + 0.07903572682584407, 0.03619864817415804, + 0.05969238281249986, 0.0790357268258422, + + -0.001075744628906913, 0.00191429223907134, + 0.07405921850311592, -0.001075744628905865, + 0.03812914216116729, 0.1346851306015185, + 0.2231273865431899, 0.01600835564431217, + 0.02200737428129396, 0.02913160002633264, + 0.02913160002633509, 0.02200737428129391, + 0.01600835564431228, 0.2231273865431878, + 0.1346851306015183, 0.03812914216116718, + + 0.002873452861657176, 0.002873452861657321, + 0.006648031513342374, 0.006648031513343037, + 0.03619864817415817, 0.05969238281250032, + 0.07903572682584409, 0.03619864817415805, + 0.05969238281249992, 0.07903572682584221, + 0.02496269311797776, 0.04081948955407392, + 0.02496269311797785, 0.1527716818820233, + 0.2348152760709258, 0.1527716818820236, + + 0.001914292239071237, -0.001075744628906803, + -0.001075744628906778, 0.07405921850311617, + 0.01600835564431228, 0.02200737428129401, + 0.02913160002633524, 0.03812914216116726, + 0.1346851306015182, 0.2231273865431886, + 0.01600835564431228, 0.02200737428129397, + 0.02913160002633523, 0.03812914216116726, + 0.1346851306015181, 0.2231273865431886, }; Assert (sizeof(loqv4)/sizeof(loqv4[0]) == @@ -823,7 +823,7 @@ MappingQ::compute_mapping_support_points( a.resize(GeometryInfo::vertices_per_cell); for (unsigned int i=0; i::vertices_per_cell; ++i) - a[i] = cell->vertex(this->vertex_mapping[i]); + a[i] = cell->vertex(i); } } @@ -853,6 +853,7 @@ MappingQ::compute_support_points_laplace(const typename Triangulation: break; case 3: + { // in 3d also add the // points located on // the boundary faces @@ -860,7 +861,7 @@ MappingQ::compute_support_points_laplace(const typename Triangulation: add_quad_support_points (cell, a); apply_laplace_vector (laplace_on_hex_vector, a); break; - + } default: Assert(false, ExcNotImplemented()); break; diff --git a/deal.II/deal.II/source/fe/mapping_q1.cc b/deal.II/deal.II/source/fe/mapping_q1.cc index 554166b344..b867bd7fb2 100644 --- a/deal.II/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/deal.II/source/fe/mapping_q1.cc @@ -82,10 +82,6 @@ MappingQ1::compute_shapes (const std::vector > &unit_points, #if (deal_II_dimension == 1) -template<> -const unsigned int MappingQ1::vertex_mapping[2] = -{ 0, 1 }; - template<> void @@ -118,10 +114,6 @@ MappingQ1<1>::compute_shapes_virtual (const std::vector > &unit_points, #if (deal_II_dimension == 2) -template<> const unsigned int -MappingQ1<2>::vertex_mapping[4] = -{ 0, 1, 3, 2 }; - template<> void @@ -163,10 +155,6 @@ MappingQ1<2>::compute_shapes_virtual (const std::vector > &unit_points, #if (deal_II_dimension == 3) -template<> -const unsigned int MappingQ1::vertex_mapping[8] = -{ 0, 1, 4, 5, 3, 2, 7, 6 }; - template<> void @@ -362,40 +350,46 @@ MappingQ1::compute_face_data (const UpdateFlags update_flags, const unsigned int nfaces = GeometryInfo::faces_per_cell; data.unit_tangentials.resize (nfaces*(dim-1), std::vector > (n_original_q_points)); - for (unsigned int i=0; inormal_directions[i]/2 + dim; - - // First tangential has a - // non-zero in component - // (i+1)%dim, if normal is - // non-zero in i. - Tensor<1,dim> tangential; - tangential[(nindex+1)%dim] = (this->normal_directions[i]%2) ? -1 : 1; - std::fill (data.unit_tangentials[i].begin(), - data.unit_tangentials[i].end(), - tangential); - - if (dim>2) + // ensure a counterclock wise + // orientation of tangentials + static const int tangential_orientation[4]={-1,1,1,-1}; + for (unsigned int i=0; i tang; + tang[1-i/2]=tangential_orientation[i]; + std::fill (data.unit_tangentials[i].begin(), + data.unit_tangentials[i].end(), tang); + } + } + else if (dim==3) + { + for (unsigned int i=0; i tangential; - tangential[(nindex-1)%dim] = 1.; - std::fill (data.unit_tangentials[i+nfaces].begin(), - data.unit_tangentials[i+nfaces].end(), - tangential); + Tensor<1,dim> tang1, tang2; + + const unsigned int nd= + GeometryInfo::unit_normal_direction[i]; + + // first tangential + // vector in direction + // of the (nd+1)%3 axis + // and inverted in case + // of unit inward normal + tang1[(nd+1)%dim]=GeometryInfo::unit_normal_orientation[i]; + // second tangential + // vector in direction + // of the (nd+2)%3 axis + tang2[(nd+2)%dim]=1.; + + // same unit tangents + // for all quadrature + // points on this face + std::fill (data.unit_tangentials[i].begin(), + data.unit_tangentials[i].end(), tang1); + std::fill (data.unit_tangentials[nfaces+i].begin(), + data.unit_tangentials[nfaces+i].end(), tang2); } } } @@ -536,7 +530,7 @@ MappingQ1::compute_mapping_support_points( a.resize(GeometryInfo::vertices_per_cell); for (unsigned int i=0; i::vertices_per_cell; ++i) - a[i] = cell->vertex(vertex_mapping[i]); + a[i] = cell->vertex(i); } diff --git a/deal.II/deal.II/source/fe/mapping_q1_eulerian.cc b/deal.II/deal.II/source/fe/mapping_q1_eulerian.cc index f58396141e..e305c1471e 100644 --- a/deal.II/deal.II/source/fe/mapping_q1_eulerian.cc +++ b/deal.II/deal.II/source/fe/mapping_q1_eulerian.cc @@ -94,12 +94,12 @@ compute_mapping_support_points(const typename Triangulation::cell_iterator // are always numbered first, // we can access them easily for (unsigned int j=0; jvertex_mapping[i]*dim+j); + shift_vector[j] = mapping_values(i*dim+j); // compute new support point by // old (reference) value and // added shift - a[i] = cell->vertex(this->vertex_mapping[i]) + shift_vector; + a[i] = cell->vertex(i) + shift_vector; } } diff --git a/deal.II/deal.II/source/grid/grid_generator.cc b/deal.II/deal.II/source/grid/grid_generator.cc index b128bb451f..8643d592cf 100644 --- a/deal.II/deal.II/source/grid/grid_generator.cc +++ b/deal.II/deal.II/source/grid/grid_generator.cc @@ -92,7 +92,8 @@ GridGenerator::hyper_rectangle (Triangulation &tria, cells[0].vertices[i] = i; cells[0].material_id = 0; - tria.create_triangulation (vertices, cells, SubCellData()); + SubCellData subcelldata; + tria.create_triangulation_compatibility (vertices, cells, subcelldata); // Assign boundary indicators if (colorize) @@ -119,25 +120,8 @@ GridGenerator::colorize_hyper_rectangle (Triangulation &tria) // there is only one cell, so // simple task const typename Triangulation::cell_iterator cell = tria.begin(); - switch(dim) - { - case 2: - cell->face(0)->set_boundary_indicator (2); - cell->face(1)->set_boundary_indicator (1); - cell->face(2)->set_boundary_indicator (3); - cell->face(3)->set_boundary_indicator (0); - break; - case 3: - cell->face(0)->set_boundary_indicator (2); - cell->face(1)->set_boundary_indicator (3); - cell->face(2)->set_boundary_indicator (4); - cell->face(3)->set_boundary_indicator (1); - cell->face(4)->set_boundary_indicator (5); - cell->face(5)->set_boundary_indicator (0); - break; - default: - Assert(false, ExcNotImplemented()); - }; + for (unsigned int f=0; f::faces_per_cell; ++f) + cell->face(f)->set_boundary_indicator (f); } #endif @@ -171,9 +155,9 @@ GridGenerator::parallelogram ( std::vector > vertices (GeometryInfo::vertices_per_cell); vertices[1] = corners[0]; - vertices[3] = corners[1]; - vertices[2] = vertices[1]; - vertices[2] += vertices[3]; + vertices[2] = corners[1]; + vertices[3] = vertices[1]; + vertices[3] += vertices[2]; // Prepare cell data std::vector > cells (1); for (unsigned int i=0;i::vertices_per_cell;++i) @@ -251,8 +235,8 @@ GridGenerator::subdivided_hyper_cube (Triangulation &tria, const unsigned int c = x+y*repetitions; cells[c].vertices[0] = y*(repetitions+1)+x; cells[c].vertices[1] = y*(repetitions+1)+x+1; - cells[c].vertices[2] = (y+1)*(repetitions+1)+x+1; - cells[c].vertices[3] = (y+1)*(repetitions+1)+x; + cells[c].vertices[2] = (y+1)*(repetitions+1)+x; + cells[c].vertices[3] = (y+1)*(repetitions+1)+x+1; cells[c].material_id = 0; } break; @@ -365,8 +349,8 @@ subdivided_hyper_rectangle (Triangulation &tria, const unsigned int c = x+y*repetitions[0]; cells[c].vertices[0] = y*(repetitions[0]+1)+x; cells[c].vertices[1] = y*(repetitions[0]+1)+x+1; - cells[c].vertices[2] = (y+1)*(repetitions[0]+1)+x+1; - cells[c].vertices[3] = (y+1)*(repetitions[0]+1)+x; + cells[c].vertices[2] = (y+1)*(repetitions[0]+1)+x; + cells[c].vertices[3] = (y+1)*(repetitions[0]+1)+x+1; cells[c].material_id = 0; } break; @@ -386,12 +370,12 @@ subdivided_hyper_rectangle (Triangulation &tria, z*repetitions[0]*repetitions[1]; cells[c].vertices[0] = z*n_xy + y*n_x + x; cells[c].vertices[1] = z*n_xy + y*n_x + x+1; - cells[c].vertices[2] = (z+1)*n_xy + y*n_x + x+1; - cells[c].vertices[3] = (z+1)*n_xy + y*n_x + x; - cells[c].vertices[4] = z*n_xy + (y+1)*n_x + x; - cells[c].vertices[5] = z*n_xy + (y+1)*n_x + x+1; - cells[c].vertices[6] = (z+1)*n_xy + (y+1)*n_x + x+1; - cells[c].vertices[7] = (z+1)*n_xy + (y+1)*n_x + x; + cells[c].vertices[2] = z*n_xy + (y+1)*n_x + x; + cells[c].vertices[3] = z*n_xy + (y+1)*n_x + x+1; + cells[c].vertices[4] = (z+1)*n_xy + y*n_x + x; + cells[c].vertices[5] = (z+1)*n_xy + y*n_x + x+1; + cells[c].vertices[6] = (z+1)*n_xy + (y+1)*n_x + x; + cells[c].vertices[7] = (z+1)*n_xy + (y+1)*n_x + x+1; cells[c].material_id = 0; } break; @@ -572,8 +556,8 @@ void GridGenerator::enclosed_hyper_cube (Triangulation<2> &tria, { cells[k].vertices[0] = i1+4*i0; cells[k].vertices[1] = i1+4*i0+1; - cells[k].vertices[2] = i1+4*i0+5; - cells[k].vertices[3] = i1+4*i0+4; + cells[k].vertices[2] = i1+4*i0+4; + cells[k].vertices[3] = i1+4*i0+5; if (colorize) cells[k].material_id = materials[k]; ++k; @@ -602,10 +586,10 @@ GridGenerator::hyper_cube_slit (Triangulation<2> &tria, Point<2>(left, right), Point<2>(right,right), Point<2>(rl2, left ) }; - const int cell_vertices[4][4] = { { 0,1,2,3 }, - { 9,4,5,2 }, - { 3,2,6,7 }, - { 2,5,8,6 } }; + const int cell_vertices[4][4] = { { 0,1,3,2 }, + { 9,4,2,5 }, + { 3,2,7,6 }, + { 2,5,6,8 } }; std::vector > cells (4, CellData<2>()); for (unsigned int i=0; i<4; ++i) { @@ -613,9 +597,10 @@ GridGenerator::hyper_cube_slit (Triangulation<2> &tria, cells[i].vertices[j] = cell_vertices[i][j]; cells[i].material_id = 0; }; - tria.create_triangulation (std::vector >(&vertices[0], &vertices[10]), - cells, - SubCellData()); // no boundary information + tria.create_triangulation ( + std::vector >(&vertices[0], &vertices[10]), + cells, + SubCellData()); // no boundary information if (colorize) { @@ -642,9 +627,9 @@ GridGenerator::hyper_L (Triangulation<2> &tria, Point (b,(a+b)/2), Point (a,b), Point ((a+b)/2,b) }; - const int cell_vertices[3][4] = {{0, 1, 4, 3}, - {1, 2, 5, 4}, - {3, 4, 7, 6}}; + const int cell_vertices[3][4] = {{0, 1, 3, 4}, + {1, 2, 4, 5}, + {3, 4, 6, 7}}; std::vector > cells (3, CellData<2>()); @@ -655,9 +640,10 @@ GridGenerator::hyper_L (Triangulation<2> &tria, cells[i].material_id = 0; }; - tria.create_triangulation (std::vector >(&vertices[0], &vertices[8]), - cells, - SubCellData()); // no boundary information + tria.create_triangulation ( + std::vector >(&vertices[0], &vertices[8]), + cells, + SubCellData()); // no boundary information } @@ -680,11 +666,11 @@ GridGenerator::hyper_ball (Triangulation<2> &tria, p+Point<2>(-1,+1)*(radius/std::sqrt(2.0)), p+Point<2>(+1,+1)*(radius/std::sqrt(2.0)) }; - const int cell_vertices[5][4] = {{0, 1, 3, 2}, - {0, 2, 4, 6}, - {2, 3, 5, 4}, - {1, 7, 5, 3}, - {6, 4, 5, 7}}; + const int cell_vertices[5][4] = {{0, 1, 2, 3}, + {0, 2, 6, 4}, + {2, 3, 4, 5}, + {1, 7, 3, 5}, + {6, 4, 7, 5}}; std::vector > cells (5, CellData<2>()); @@ -695,9 +681,10 @@ GridGenerator::hyper_ball (Triangulation<2> &tria, cells[i].material_id = 0; }; - tria.create_triangulation (std::vector >(&vertices[0], &vertices[8]), - cells, - SubCellData()); // no boundary information + tria.create_triangulation ( + std::vector >(&vertices[0], &vertices[8]), + cells, + SubCellData()); // no boundary information } @@ -751,13 +738,14 @@ void GridGenerator::hyper_shell (Triangulation<2> &tria, { cells[i].vertices[0] = i; cells[i].vertices[1] = (i+1)%N; - cells[i].vertices[2] = N+((i+1)%N); - cells[i].vertices[3] = N+i; + cells[i].vertices[2] = N+i; + cells[i].vertices[3] = N+((i+1)%N); cells[i].material_id = 0; }; - tria.create_triangulation (vertices, cells, SubCellData()); + tria.create_triangulation ( + vertices, cells, SubCellData()); } @@ -813,10 +801,10 @@ GridGenerator::half_hyper_ball (Triangulation<2> &tria, p+Point<2>(0,+1)*radius, p+Point<2>(+1,+1)*(radius/std::sqrt(2.0)) }; - const int cell_vertices[5][4] = {{0, 1, 3, 2}, - {2, 3, 5, 4}, - {1, 7, 5, 3}, - {6, 4, 5, 7}}; + const int cell_vertices[5][4] = {{0, 1, 2, 3}, + {2, 3, 4, 5}, + {1, 7, 3, 5}, + {6, 4, 7, 5}}; std::vector > cells (4, CellData<2>()); @@ -827,9 +815,10 @@ GridGenerator::half_hyper_ball (Triangulation<2> &tria, cells[i].material_id = 0; }; - tria.create_triangulation (std::vector >(&vertices[0], &vertices[8]), - cells, - SubCellData()); // no boundary information + tria.create_triangulation ( + std::vector >(&vertices[0], &vertices[8]), + cells, + SubCellData()); // no boundary information } @@ -892,8 +881,8 @@ GridGenerator::half_hyper_shell (Triangulation<2> &tria, { cells[i].vertices[0] = i; cells[i].vertices[1] = (i+1)%(N+1); - cells[i].vertices[2] = N+1+((i+1)%(N+1)); - cells[i].vertices[3] = N+1+i; + cells[i].vertices[2] = N+1+i; + cells[i].vertices[3] = N+1+((i+1)%(N+1)); cells[i].material_id = 0; }; @@ -958,19 +947,20 @@ void GridGenerator::enclosed_hyper_cube (Triangulation<3> &tria, { cells[k].vertices[0] = i2+4*i1+16*i0; cells[k].vertices[1] = i2+4*i1+16*i0+1; - cells[k].vertices[2] = i2+4*i1+16*i0+5; - cells[k].vertices[3] = i2+4*i1+16*i0+4; - cells[k].vertices[4] = i2+4*i1+16*i0+16; - cells[k].vertices[5] = i2+4*i1+16*i0+17; - cells[k].vertices[6] = i2+4*i1+16*i0+21; - cells[k].vertices[7] = i2+4*i1+16*i0+20; + cells[k].vertices[2] = i2+4*i1+16*i0+16; + cells[k].vertices[3] = i2+4*i1+16*i0+17; + cells[k].vertices[4] = i2+4*i1+16*i0+4; + cells[k].vertices[5] = i2+4*i1+16*i0+5; + cells[k].vertices[6] = i2+4*i1+16*i0+20; + cells[k].vertices[7] = i2+4*i1+16*i0+21; if (colorize) cells[k].material_id = materials[k]; ++k; } - tria.create_triangulation (vertices, - cells, - SubCellData()); // no boundary information + tria.create_triangulation ( + vertices, + cells, + SubCellData()); // no boundary information } @@ -1032,10 +1022,12 @@ GridGenerator::hyper_L (Triangulation<3> &tria, cells[i].vertices[j] = cell_vertices[i][j]; cells[i].material_id = 0; }; - - tria.create_triangulation (std::vector >(&vertices[0], &vertices[26]), - cells, - SubCellData()); // no boundary information + + SubCellData subcelldata; + tria.create_triangulation_compatibility ( + std::vector >(&vertices[0], &vertices[26]), + cells, + subcelldata); // no boundary information } @@ -1094,10 +1086,12 @@ GridGenerator::hyper_ball (Triangulation<3> &tria, }; GridReordering<3>::reorder_cells (cells); - - tria.create_triangulation (std::vector >(&vertices[0], &vertices[n_vertices]), - cells, - SubCellData()); // no boundary information + + SubCellData subcelldata; + tria.create_triangulation_compatibility ( + std::vector >(&vertices[0], &vertices[n_vertices]), + cells, + subcelldata); // no boundary information } @@ -1166,10 +1160,12 @@ GridGenerator::cylinder (Triangulation<3> &tria, }; GridReordering<3>::reorder_cells (cells); - - tria.create_triangulation (std::vector >(&vertices[0], &vertices[24]), - cells, - SubCellData()); // no boundary information + + SubCellData subcelldata; + tria.create_triangulation_compatibility ( + std::vector >(&vertices[0], &vertices[24]), + cells, + subcelldata); // no boundary information Triangulation<3>::cell_iterator cell = tria.begin(); Triangulation<3>::cell_iterator end = tria.end(); @@ -1267,18 +1263,19 @@ void GridGenerator::cylinder_shell (Triangulation<3> &tria, { cells[i+j*N_r].vertices[0] = i + (j+1)*2*N_r; cells[i+j*N_r].vertices[1] = (i+1)%N_r + (j+1)*2*N_r; - cells[i+j*N_r].vertices[2] = N_r+((i+1)%N_r) + (j+1)*2*N_r; - cells[i+j*N_r].vertices[3] = N_r+i + (j+1)*2*N_r; - - cells[i+j*N_r].vertices[4] = i + j*2*N_r; - cells[i+j*N_r].vertices[5] = (i+1)%N_r + j*2*N_r; - cells[i+j*N_r].vertices[6] = N_r+((i+1)%N_r) + j*2*N_r; - cells[i+j*N_r].vertices[7] = N_r+i + j*2*N_r; + cells[i+j*N_r].vertices[2] = i + j*2*N_r; + cells[i+j*N_r].vertices[3] = (i+1)%N_r + j*2*N_r; + + cells[i+j*N_r].vertices[4] = N_r+i + (j+1)*2*N_r; + cells[i+j*N_r].vertices[5] = N_r+((i+1)%N_r) + (j+1)*2*N_r; + cells[i+j*N_r].vertices[6] = N_r+i + j*2*N_r; + cells[i+j*N_r].vertices[7] = N_r+((i+1)%N_r) + j*2*N_r; cells[i+j*N_r].material_id = 0; } - tria.create_triangulation (vertices_3d, cells, SubCellData()); + tria.create_triangulation ( + vertices_3d, cells, SubCellData()); } diff --git a/deal.II/deal.II/source/grid/grid_in.cc b/deal.II/deal.II/source/grid/grid_in.cc index 6582f69b0d..e7b84e2aec 100644 --- a/deal.II/deal.II/source/grid/grid_in.cc +++ b/deal.II/deal.II/source/grid/grid_in.cc @@ -201,7 +201,7 @@ void GridIn::read_ucd (std::istream &in) // ... and cells GridReordering::invert_all_cells_of_negative_grid (vertices, cells); GridReordering::reorder_cells (cells); - tria->create_triangulation (vertices, cells, subcelldata); + tria->create_triangulation_compatibility (vertices, cells, subcelldata); } @@ -366,7 +366,7 @@ void GridIn::read_dbmesh (std::istream &in) // ...and cells GridReordering::invert_all_cells_of_negative_grid (vertices, cells); GridReordering::reorder_cells (cells); - tria->create_triangulation (vertices, cells, subcelldata); + tria->create_triangulation_compatibility (vertices, cells, subcelldata); } @@ -448,7 +448,7 @@ void GridIn<2>::read_xda (std::istream &in) // ... and cells GridReordering<2>::invert_all_cells_of_negative_grid (vertices, cells); GridReordering<2>::reorder_cells (cells); - tria->create_triangulation (vertices, cells, subcelldata); + tria->create_triangulation_compatibility (vertices, cells, subcelldata); } #endif // #if deal_II_dimension == 2 @@ -531,7 +531,7 @@ void GridIn<3>::read_xda (std::istream &in) // ... and cells GridReordering<3>::invert_all_cells_of_negative_grid (vertices, cells); GridReordering<3>::reorder_cells (cells); - tria->create_triangulation (vertices, cells, subcelldata); + tria->create_triangulation_compatibility (vertices, cells, subcelldata); } #endif // #if deal_II_dimension == 3 @@ -735,7 +735,7 @@ void GridIn::read_msh (std::istream &in) // ... and cells GridReordering::invert_all_cells_of_negative_grid (vertices, cells); GridReordering::reorder_cells (cells); - tria->create_triangulation (vertices, cells, subcelldata); + tria->create_triangulation_compatibility (vertices, cells, subcelldata); } @@ -1002,7 +1002,7 @@ void GridIn<2>::read_netcdf (const std::string &filename) SubCellData subcelldata; GridTools::delete_unused_vertices(vertices, cells, subcelldata); GridReordering::reorder_cells (cells); - tria->create_triangulation (vertices, cells, subcelldata); + tria->create_triangulation_compatibility (vertices, cells, subcelldata); #endif } @@ -1212,7 +1212,7 @@ void GridIn<3>::read_netcdf (const std::string &filename) if (switch_y_z) GridReordering::invert_all_cells_of_negative_grid (vertices, cells); GridReordering::reorder_cells (cells); - tria->create_triangulation (vertices, cells, subcelldata); + tria->create_triangulation_compatibility (vertices, cells, subcelldata); #endif } diff --git a/deal.II/deal.II/source/grid/grid_out.cc b/deal.II/deal.II/source/grid/grid_out.cc index c2983f7d4e..56423878fa 100644 --- a/deal.II/deal.II/source/grid/grid_out.cc +++ b/deal.II/deal.II/source/grid/grid_out.cc @@ -352,10 +352,10 @@ void GridOut::write_ucd (const Triangulation &tria, // AVS Developer's Guide, Release 4, // May, 1992, p. E6 // - // note: vertex numbers are 1-base + // note: vertex numbers are 1-base for (unsigned int vertex=0; vertex::vertices_per_cell; ++vertex) - out << cell->vertex_index(vertex)+1 << ' '; + out << cell->vertex_index(GeometryInfo::ucd_to_deal[vertex])+1 << ' '; out << std::endl; }; @@ -443,7 +443,8 @@ void GridOut::write_xfig (const Triangulation& tria, // (dots/inch) for (unsigned int k=0;k<=nv;++k) { - const Point& p = cell->vertex(k % nv); + const Point& p = cell->vertex( + GeometryInfo::ucd_to_deal[k % nv]); for (unsigned int d=0;d& tria, out << std::endl; } // Now write boundary edges + static const unsigned int face_reorder[4]={2,1,3,0}; if (xfig_flags.draw_boundary) for (unsigned int f=0;f::face_iterator - face = cell->face(f); + face = cell->face(face_reorder[f]); const unsigned char bi = face->boundary_indicator(); if (bi != 255) { @@ -571,7 +573,7 @@ void GridOut::write_ucd_faces (const Triangulation &tria, }; // note: vertex numbers are 1-base for (unsigned int vertex=0; vertex::vertices_per_face; ++vertex) - out << face->vertex_index(vertex)+1 << ' '; + out << face->vertex_index(GeometryInfo::ucd_to_deal[vertex])+1 << ' '; out << std::endl; ++index; @@ -662,7 +664,7 @@ void GridOut::write_gnuplot (const Triangulation<2> &tria, // drawing pencil at the // end for (unsigned int i=0; i::vertices_per_cell; ++i) - out << cell->vertex(i) + out << cell->vertex(GeometryInfo::ucd_to_deal[i]) << ' ' << cell->level() << ' ' << static_cast(cell->material_id()) << std::endl; out << cell->vertex(0) @@ -779,7 +781,7 @@ void GridOut::write_gnuplot (const Triangulation<3> &tria, out << "# cell " << cell << std::endl; if (mapping==0 || n_points==2 || !cell->has_boundary_lines()) - { + { // front face out << cell->vertex(0) << ' ' << cell->level() @@ -787,10 +789,10 @@ void GridOut::write_gnuplot (const Triangulation<3> &tria, << cell->vertex(1) << ' ' << cell->level() << ' ' << static_cast(cell->material_id()) << std::endl - << cell->vertex(2) + << cell->vertex(5) << ' ' << cell->level() << ' ' << static_cast(cell->material_id()) << std::endl - << cell->vertex(3) + << cell->vertex(4) << ' ' << cell->level() << ' ' << static_cast(cell->material_id()) << std::endl << cell->vertex(0) @@ -798,19 +800,19 @@ void GridOut::write_gnuplot (const Triangulation<3> &tria, << ' ' << static_cast(cell->material_id()) << std::endl << std::endl; // back face - out << cell->vertex(4) + out << cell->vertex(2) << ' ' << cell->level() << ' ' << static_cast(cell->material_id()) << std::endl - << cell->vertex(5) + << cell->vertex(3) << ' ' << cell->level() << ' ' << static_cast(cell->material_id()) << std::endl - << cell->vertex(6) + << cell->vertex(7) << ' ' << cell->level() << ' ' << static_cast(cell->material_id()) << std::endl - << cell->vertex(7) + << cell->vertex(6) << ' ' << cell->level() << ' ' << static_cast(cell->material_id()) << std::endl - << cell->vertex(4) + << cell->vertex(2) << ' ' << cell->level() << ' ' << static_cast(cell->material_id()) << std::endl << std::endl; @@ -819,28 +821,28 @@ void GridOut::write_gnuplot (const Triangulation<3> &tria, out << cell->vertex(0) << ' ' << cell->level() << ' ' << static_cast(cell->material_id()) << std::endl - << cell->vertex(4) + << cell->vertex(2) << ' ' << cell->level() << ' ' << static_cast(cell->material_id()) << std::endl << std::endl; out << cell->vertex(1) << ' ' << cell->level() << ' ' << static_cast(cell->material_id()) << std::endl - << cell->vertex(5) + << cell->vertex(3) << ' ' << cell->level() << ' ' << static_cast(cell->material_id()) << std::endl << std::endl; - out << cell->vertex(2) + out << cell->vertex(5) << ' ' << cell->level() << ' ' << static_cast(cell->material_id()) << std::endl - << cell->vertex(6) + << cell->vertex(7) << ' ' << cell->level() << ' ' << static_cast(cell->material_id()) << std::endl << std::endl; - out << cell->vertex(3) + out << cell->vertex(4) << ' ' << cell->level() << ' ' << static_cast(cell->material_id()) << std::endl - << cell->vertex(7) + << cell->vertex(6) << ' ' << cell->level() << ' ' << static_cast(cell->material_id()) << std::endl << std::endl; diff --git a/deal.II/deal.II/source/grid/grid_reordering.cc b/deal.II/deal.II/source/grid/grid_reordering.cc index ec33913760..33c137c9be 100644 --- a/deal.II/deal.II/source/grid/grid_reordering.cc +++ b/deal.II/deal.II/source/grid/grid_reordering.cc @@ -1451,24 +1451,34 @@ GridReordering<3>::invert_all_cells_of_negative_grid( const std::vector > &all_vertices, std::vector > &cells) { + int vertices_lex[GeometryInfo<3>::vertices_per_cell]; unsigned int n_negative_cells=0; for (unsigned int cell_no=0; cell_no 0, - ExcInternalError()); - } - + { + // GridTools::cell_measure + // requires the vertices to be + // in lexicographic ordering + for (unsigned int i=0; i::vertices_per_cell; ++i) + vertices_lex[GeometryInfo<3>::ucd_to_deal[i]]=cells[cell_no].vertices[i]; + if (GridTools::cell_measure(all_vertices, vertices_lex) < 0) + { + ++n_negative_cells; + for (unsigned int i=0; i<4; ++i) + std::swap(cells[cell_no].vertices[i], cells[cell_no].vertices[i+4]); + + // check whether the + // resulting cell is now ok. + // if not, then the grid is + // seriously broken and + // should be sticked into the + // bin + for (unsigned int i=0; i::vertices_per_cell; ++i) + vertices_lex[GeometryInfo<3>::ucd_to_deal[i]]=cells[cell_no].vertices[i]; + AssertThrow(GridTools::cell_measure(all_vertices, vertices_lex) > 0, + ExcInternalError()); + } + } + // We assuming that all cells of a // grid have either positive or // negative volumes but not both diff --git a/deal.II/deal.II/source/grid/grid_tools.cc b/deal.II/deal.II/source/grid/grid_tools.cc index d9350c3d75..1666ef9c80 100644 --- a/deal.II/deal.II/source/grid/grid_tools.cc +++ b/deal.II/deal.II/source/grid/grid_tools.cc @@ -109,6 +109,12 @@ double GridTools::cell_measure(const std::vector > &all_vertices, const int vertex_indices[GeometryInfo<3>::vertices_per_cell]) { + // note that this is the + // cell_measure based on the new + // deal.II numbering. When called + // from inside GridReordering make + // sure that you reorder the + // vertex_indices before const double x[8] = { all_vertices[vertex_indices[0]](0), all_vertices[vertex_indices[1]](0), all_vertices[vertex_indices[2]](0), @@ -135,60 +141,85 @@ GridTools::cell_measure(const std::vector > &all_vertices, all_vertices[vertex_indices[7]](2) }; /* - Get the computation of the measure by the same little Maple script - as in TriaObjectAccessor<3, 3>::barycenter(), see tria_accessor.cc. + This is the same Maple script as in the barycenter method above + except of that here the shape functions tphi[0]-tphi[7] are ordered + according to the lexicographic numbering. + + x := array(0..7): + y := array(0..7): + z := array(0..7): + tphi[0] := (1-xi)*(1-eta)*(1-zeta): + tphi[1] := xi*(1-eta)*(1-zeta): + tphi[2] := (1-xi)* eta*(1-zeta): + tphi[3] := xi* eta*(1-zeta): + tphi[4] := (1-xi)*(1-eta)*zeta: + tphi[5] := xi*(1-eta)*zeta: + tphi[6] := (1-xi)* eta*zeta: + tphi[7] := xi* eta*zeta: + x_real := sum(x[s]*tphi[s], s=0..7): + y_real := sum(y[s]*tphi[s], s=0..7): + z_real := sum(z[s]*tphi[s], s=0..7): + with (linalg): + J := matrix(3,3, [[diff(x_real, xi), diff(x_real, eta), diff(x_real, zeta)], + [diff(y_real, xi), diff(y_real, eta), diff(y_real, zeta)], + [diff(z_real, xi), diff(z_real, eta), diff(z_real, zeta)]]): + detJ := det (J): + + measure := simplify ( int ( int ( int (detJ, xi=0..1), eta=0..1), zeta=0..1)): + + readlib(C): + + C(measure, optimized); + + The C code produced by this maple script is further optimized by + hand. In particular, division by 12 is performed only once, not + hundred of times. */ - double s1, s2, s3; - - s3 = -x[2]*y[7]*z[3]/12.0+x[1]*y[4]*z[0]/12.0-x[0]*y[3]*z[7]/12.0-z[1]*x - [2]*y[0]/12.0+x[5]*y[4]*z[0]/12.0+x[1]*y[5]*z[0]/12.0-x[1]*y[2]*z[0]/12.0-y[2]* - x[3]*z[7]/12.0-x[1]*y[0]*z[4]/12.0-y[1]*x[0]*z[2]/12.0-y[0]*x[7]*z[3]/12.0+y[5] - *x[4]*z[7]/12.0-x[0]*y[3]*z[4]/12.0-y[1]*x[0]*z[3]/12.0-z[0]*x[7]*y[4]/12.0+x - [1]*y[0]*z[3]/12.0-x[2]*y[3]*z[0]/12.0-y[1]*x[2]*z[5]/12.0; - s2 = s3+y[0]*x[3]*z[4]/12.0+y[1]*x[2]*z[0]/12.0+y[2]*x[7]*z[3]/12.0+z[1]* - x[4]*y[0]/12.0-y[5]*x[4]*z[0]/12.0-y[1]*x[5]*z[0]/12.0+x[5]*y[7]*z[4]/12.0+y[1] - *x[3]*z[0]/12.0+x[0]*y[4]*z[7]/12.0+z[0]*x[7]*y[3]/12.0-z[2]*x[3]*y[0]/12.0+y - [1]*x[0]*z[5]/12.0-x[1]*y[5]*z[2]/12.0+y[1]*x[0]*z[4]/12.0+z[1]*x[2]*y[5]/12.0- - x[5]*y[0]*z[4]/12.0+y[2]*x[3]*z[1]/12.0+x[1]*y[0]*z[2]/12.0; - s3 = -z[1]*x[0]*y[4]/12.0+z[1]*x[0]*y[2]/12.0-x[0]*y[7]*z[4]/12.0+z[0]*x - [4]*y[7]/12.0+x[2]*y[3]*z[7]/12.0+y[2]*x[6]*z[3]/12.0+y[1]*x[6]*z[2]/12.0-y[2]* - x[1]*z[3]/12.0+x[1]*y[2]*z[5]/12.0+y[0]*x[3]*z[7]/12.0+y[2]*x[3]*z[0]/12.0-y[0] - *x[4]*z[3]/12.0-y[2]*x[0]*z[3]/12.0-y[1]*x[4]*z[0]/12.0+z[1]*x[5]*y[0]/12.0+x - [1]*y[2]*z[6]/12.0-x[1]*y[0]*z[5]/12.0-x[2]*y[3]*z[1]/12.0; - s1 = s3-z[0]*x[3]*y[4]/12.0+z[1]*x[0]*y[3]/12.0-z[1]*x[5]*y[2]/12.0+z[1]* - x[2]*y[6]/12.0-z[2]*x[3]*y[1]/12.0+z[2]*x[3]*y[6]/12.0+x[2]*y[1]*z[3]/12.0-z[1] - *x[3]*y[0]/12.0+x[2]*y[3]*z[6]/12.0-z[0]*x[3]*y[7]/12.0+y[5]*x[0]*z[4]/12.0-z - [1]*x[6]*y[2]/12.0-z[2]*x[6]*y[3]/12.0+z[2]*x[0]*y[3]/12.0+z[2]*x[3]*y[7]/12.0+ - z[2]*x[1]*y[3]/12.0+y[1]*x[5]*z[2]/12.0-z[2]*x[7]*y[3]/12.0+s2; - s3 = x[0]*y[7]*z[3]/12.0-z[1]*x[0]*y[5]/12.0-x[1]*y[3]*z[0]/12.0-x[1]*y - [6]*z[2]/12.0-x[2]*y[6]*z[3]/12.0-x[5]*y[4]*z[7]/12.0+z[0]*x[4]*y[3]/12.0+x[0]* - y[4]*z[3]/12.0-x[6]*y[5]*z[7]/12.0-x[6]*y[7]*z[2]/12.0+z[1]*x[6]*y[5]/12.0+y[6] - *x[4]*z[7]/12.0+y[1]*x[5]*z[6]/12.0+x[1]*y[6]*z[5]/12.0-y[5]*x[7]*z[4]/12.0+y - [0]*x[7]*z[4]/12.0-y[0]*x[4]*z[7]/12.0-z[5]*x[0]*y[4]/12.0; - s2 = s3+z[6]*x[3]*y[7]/12.0+z[6]*x[7]*y[4]/12.0-z[6]*x[4]*y[7]/12.0-z[6]* - x[7]*y[3]/12.0+y[6]*x[5]*z[7]/12.0-y[6]*x[7]*z[5]/12.0-x[2]*y[5]*z[6]/12.0+z[6] - *x[2]*y[7]/12.0+z[6]*x[7]*y[5]/12.0-z[6]*x[5]*y[7]/12.0-z[6]*x[7]*y[2]/12.0+x - [6]*y[7]*z[5]/12.0+z[5]*x[7]*y[4]/12.0-z[5]*x[4]*y[7]/12.0+y[5]*x[1]*z[4]/12.0- - y[5]*x[6]*z[4]/12.0-y[5]*x[4]*z[1]/12.0+y[5]*x[4]*z[6]/12.0; - s3 = y[6]*x[7]*z[3]/12.0+x[6]*y[5]*z[2]/12.0-x[5]*y[6]*z[2]/12.0-x[6]*y - [2]*z[5]/12.0+x[5]*y[2]*z[6]/12.0+x[6]*y[2]*z[7]/12.0+z[5]*x[4]*y[0]/12.0-y[6]* - x[2]*z[7]/12.0+x[2]*y[6]*z[5]/12.0+y[6]*x[7]*z[2]/12.0-y[6]*x[3]*z[7]/12.0-x[7] - *y[4]*z[3]/12.0-x[6]*y[4]*z[7]/12.0-x[5]*y[1]*z[4]/12.0+x[6]*y[7]*z[4]/12.0+x - [7]*y[3]*z[4]/12.0-z[1]*x[5]*y[6]/12.0+x[3]*y[4]*z[7]/12.0+x[5]*y[6]*z[4]/12.0; - - // the vertices entered into the - // above field and the vertices - // used by the maple script use - // different orderings (i.e. one is - // the inside-out orientation of - // the other one). the measure is - // thus negative. Thus take the - // negative value of the result - return -(s3+x[5]*y[4]*z[1]/12.0-x[5]*y[4]*z[6]/12.0-x[3]*y[7]*z[4]/12.0+x[6]* - y[3]*z[7]/12.0-y[1]*x[6]*z[5]/12.0-z[5]*x[1]*y[4]/12.0+z[5]*x[6]*y[4]/12.0+z[5] - *x[4]*y[1]/12.0-z[5]*x[4]*y[6]/12.0-x[6]*y[7]*z[3]/12.0-x[4]*y[3]*z[7]/12.0+x - [4]*y[7]*z[3]/12.0-x[1]*y[5]*z[6]/12.0-y[6]*x[7]*z[4]/12.0-y[1]*x[2]*z[6]/12.0- - y[2]*x[3]*z[6]/12.0+s2+s1+x[2]*y[0]*z[3]/12.0); + + const double t3 = y[3]*x[2]; + const double t5 = z[1]*x[5]; + const double t9 = z[3]*x[2]; + const double t11 = x[1]*y[0]; + const double t14 = x[4]*y[0]; + const double t18 = x[5]*y[7]; + const double t20 = y[1]*x[3]; + const double t22 = y[5]*x[4]; + const double t26 = z[7]*x[6]; + const double t28 = x[0]*y[4]; + const double t34 = z[3]*x[1]*y[2]+t3*z[1]-t5*y[7]+y[7]*x[4]*z[6]+t9*y[6]-t11*z[4]-t5*y[3]-t14*z[2]+z[1]*x[4]*y[0]-t18*z[3]+t20*z[0]-t22*z[0]-y[0]*x[5]*z[4]-t26*y[3]+t28*z[2]-t9*y[1]-y[1]*x[4]*z[0]-t11*z[5]; + const double t37 = y[1]*x[0]; + const double t44 = x[1]*y[5]; + const double t46 = z[1]*x[0]; + const double t49 = x[0]*y[2]; + const double t52 = y[5]*x[7]; + const double t54 = x[3]*y[7]; + const double t56 = x[2]*z[0]; + const double t58 = x[3]*y[2]; + const double t64 = -x[6]*y[4]*z[2]-t37*z[2]+t18*z[6]-x[3]*y[6]*z[2]+t11*z[2]+t5*y[0]+t44*z[4]-t46*y[4]-t20*z[7]-t49*z[6]-t22*z[1]+t52*z[3]-t54*z[2]-t56*y[4]-t58*z[0]+y[1]*x[2]*z[0]+t9*y[7]+t37*z[4]; + const double t66 = x[1]*y[7]; + const double t68 = y[0]*x[6]; + const double t70 = x[7]*y[6]; + const double t73 = z[5]*x[4]; + const double t76 = x[6]*y[7]; + const double t90 = x[4]*z[0]; + const double t92 = x[1]*y[3]; + const double t95 = -t66*z[3]-t68*z[2]-t70*z[2]+t26*y[5]-t73*y[6]-t14*z[6]+t76*z[2]-t3*z[6]+x[6]*y[2]*z[4]-z[3]*x[6]*y[2]+t26*y[4]-t44*z[3]-x[1]*y[2]*z[0]+x[5]*y[6]*z[4]+t54*z[5]+t90*y[2]-t92*z[2]+t46*y[2]; + const double t102 = x[2]*y[0]; + const double t107 = y[3]*x[7]; + const double t114 = x[0]*y[6]; + const double t125 = y[0]*x[3]*z[2]-z[7]*x[5]*y[6]-x[2]*y[6]*z[4]+t102*z[6]-t52*z[6]+x[2]*y[4]*z[6]-t107*z[5]-t54*z[6]+t58*z[6]-x[7]*y[4]*z[6]+t37*z[5]-t114*z[4]+t102*z[4]-z[1]*x[2]*y[0]+t28*z[6]-y[5]*x[6]*z[4]-z[5]*x[1]*y[4]-t73*y[7]; + const double t129 = z[0]*x[6]; + const double t133 = y[1]*x[7]; + const double t145 = y[1]*x[5]; + const double t156 = t90*y[6]-t129*y[4]+z[7]*x[2]*y[6]-t133*z[5]+x[5]*y[3]*z[7]-t26*y[2]-t70*z[3]+t46*y[3]+z[5]*x[7]*y[4]+z[7]*x[3]*y[6]-t49*z[4]+t145*z[7]-x[2]*y[7]*z[6]+t70*z[5]+t66*z[5]-z[7]*x[4]*y[6]+t18*z[4]+x[1]*y[4]*z[0]; + const double t160 = x[5]*y[4]; + const double t165 = z[1]*x[7]; + const double t178 = z[1]*x[3]; + const double t181 = t107*z[6]+t22*z[7]+t76*z[3]+t160*z[1]-x[4]*y[2]*z[6]+t70*z[4]+t165*y[5]+x[7]*y[2]*z[6]-t76*z[5]-t76*z[4]+t133*z[3]-t58*z[1]+y[5]*x[0]*z[4]+t114*z[2]-t3*z[7]+t20*z[2]+t178*y[7]+t129*y[2]; + const double t207 = t92*z[7]+t22*z[6]+z[3]*x[0]*y[2]-x[0]*y[3]*z[2]-z[3]*x[7]*y[2]-t165*y[3]-t9*y[0]+t58*z[7]+y[3]*x[6]*z[2]+t107*z[2]+t73*y[0]-x[3]*y[5]*z[7]+t3*z[0]-t56*y[6]-z[5]*x[0]*y[4]+t73*y[1]-t160*z[6]+t160*z[0]; + const double t228 = -t44*z[7]+z[5]*x[6]*y[4]-t52*z[4]-t145*z[4]+t68*z[4]+t92*z[5]-t92*z[0]+t11*z[3]+t44*z[0]+t178*y[5]-t46*y[5]-t178*y[0]-t145*z[0]-t20*z[5]-t37*z[3]-t160*z[7]+t145*z[3]+x[4]*y[6]*z[2]; + + return (t34+t64+t95+t125+t156+t181+t207+t228)/12.; } #else diff --git a/deal.II/deal.II/source/grid/persistent_tria.cc b/deal.II/deal.II/source/grid/persistent_tria.cc index 0aaf528474..716b38cc90 100644 --- a/deal.II/deal.II/source/grid/persistent_tria.cc +++ b/deal.II/deal.II/source/grid/persistent_tria.cc @@ -127,6 +127,17 @@ PersistentTriangulation::create_triangulation (const std::vector } +template +void +PersistentTriangulation::create_triangulation_compatibility ( + const std::vector > &, + const std::vector > &, + const SubCellData &) +{ + Assert (false, ExcImpossibleInDim(dim)); +} + + template void PersistentTriangulation::write_flags(std::ostream &out) const diff --git a/deal.II/deal.II/source/grid/tria.cc b/deal.II/deal.II/source/grid/tria.cc index 3f9eb96d73..6b5073737f 100644 --- a/deal.II/deal.II/source/grid/tria.cc +++ b/deal.II/deal.II/source/grid/tria.cc @@ -246,6 +246,17 @@ void Triangulation::copy_triangulation (const Triangulation &old_tria) #if deal_II_dimension == 1 +template <> +void Triangulation<1>::create_triangulation_compatibility ( + const std::vector > &v, + const std::vector > &cells, + const SubCellData &subcelldata) +{ + // no reordering required + create_triangulation(v, cells, subcelldata); +} + + template <> void Triangulation<1>::create_triangulation (const std::vector > &v, const std::vector > &cells, @@ -377,6 +388,46 @@ void Triangulation<1>::create_triangulation (const std::vector > &v, #if deal_II_dimension == 2 +template <> +void Triangulation<2>::create_triangulation_compatibility ( + const std::vector > &v, + const std::vector > &c, + const SubCellData &subcelldata) +{ + const bool output=false; + + // copy cell data as it will be + // modified in the following + std::vector > cells(c); + + for (unsigned int cell=0; cell void Triangulation<2>::create_triangulation (const std::vector > &v, const std::vector > &cells, @@ -716,6 +767,84 @@ struct QuadComparator }; +template <> +void +Triangulation<3>::create_triangulation_compatibility ( + const std::vector > &v, + const std::vector > &c, + const SubCellData &subc) +{ + const bool output1=false; + const bool output2=false; + + // copy cell and subcell data as + // they will be modified in the + // following + std::vector > cells(c); + SubCellData subcelldata(subc); + + unsigned int tmp[GeometryInfo<3>::vertices_per_cell]; + for (unsigned int cell=0; cell::vertices_per_cell; ++vertex) + std::cout << cells[cell].vertices[vertex] << " "; + } + + for (unsigned int i=0; i::vertices_per_cell; ++i) + tmp[i]=cells[cell].vertices[i]; + for (unsigned int i=0; i::vertices_per_cell; ++i) + cells[cell].vertices[GeometryInfo<3>::ucd_to_deal[i]]=tmp[i]; + + if (output1) + { + std::cout << " after swap: "; + for (unsigned int vertex=0; + vertex::vertices_per_cell; ++vertex) + std::cout << cells[cell].vertices[vertex] << " "; + std::cout << std::endl; + std::cout << "vertices:" << std::endl; + for (unsigned int vertex=0; + vertex::vertices_per_cell; ++vertex) + std::cout << v[cells[cell].vertices[vertex]] << std::endl; + } + } + + if (output2) + { + std::cout << "vertices:" << std::endl; + for (unsigned int i=0; i >::iterator boundary_quad + = subcelldata.boundary_quads.begin(); + std::vector >::iterator end_quad + = subcelldata.boundary_quads.end(); + for (unsigned int quad_no=0; boundary_quad!=end_quad; ++boundary_quad, ++quad_no) + { + if (output2) + std::cout << quad_no << ": " << boundary_quad->vertices[0] + << ", " << boundary_quad->vertices[1] + << ", " << boundary_quad->vertices[2] + << ", " << boundary_quad->vertices[3]; + + std::swap(boundary_quad->vertices[2], boundary_quad->vertices[3]); + + if (output2) + std::cout << ", after swap: " << boundary_quad->vertices[0] + << ", " << boundary_quad->vertices[1] + << ", " << boundary_quad->vertices[2] + << ", " << boundary_quad->vertices[3] << std::endl; + } + + create_triangulation(v, cells, subcelldata); +} + + template <> void Triangulation<3>::create_triangulation (const std::vector > &v, @@ -966,12 +1095,12 @@ Triangulation<3>::create_triangulation (const std::vector > &v, // set. thus, vertex 0 is // the one from which two // edges originate, and - // vertex 2 is the one to + // vertex 3 is the one to // which they converge. we // are then left with // orientations 0-1-2-3 and - // 0-3-2-1 for the order of - // vertices. the + // 2-3-0-1 for the order of + // lines. the // corresponding quad can // be easily constructed by // exchanging lines. we do @@ -983,8 +1112,8 @@ Triangulation<3>::create_triangulation (const std::vector > &v, // new one and instead // later set the // face_orientation flag - const Quad test_quad (quad.line(3), quad.line(2), - quad.line(1), quad.line(0)); + const Quad test_quad (quad.line(2), quad.line(3), + quad.line(0), quad.line(1)); if (needed_quads.find (test_quad) == needed_quads.end()) needed_quads[quad] = end_quad(); } @@ -1073,8 +1202,8 @@ Triangulation<3>::create_triangulation (const std::vector > &v, // face must be // available in // reverse order then - const Quad test_quad (quad.line(3), quad.line(2), - quad.line(1), quad.line(0)); + const Quad test_quad (quad.line(2), quad.line(3), + quad.line(0), quad.line(1)); Assert (needed_quads.find (test_quad) != needed_quads.end(), ExcInternalError()); face_iterator[face] = needed_quads[test_quad]; @@ -1155,11 +1284,11 @@ Triangulation<3>::create_triangulation (const std::vector > &v, // that line_no in // standard face // orientation is - // 3-line_no in + // (line_no+2)%4 in // non-standard face // orientation - Assert (face_iterator[face1]->line(face_orientation[face1] ? line1 : 3-line1) == - face_iterator[face2]->line(face_orientation[face2] ? line2 : 3-line2), + Assert (face_iterator[face1]->line(face_orientation[face1] ? line1 : (line1+2)%4) == + face_iterator[face2]->line(face_orientation[face2] ? line2 : (line2+2)%4), ExcInternalErrorOnCell(c)); } #endif @@ -1298,6 +1427,7 @@ Triangulation<3>::create_triangulation (const std::vector > &v, line->set_boundary_indicator (boundary_line->material_id); } + // now go on with boundary faces std::vector >::const_iterator boundary_quad = subcelldata.boundary_quads.begin(); @@ -1358,7 +1488,7 @@ Triangulation<3>::create_triangulation (const std::vector > &v, // orientation does not carry // any information. Quad quad_compare_1(line[0]->index(), line[1]->index(), line[2]->index(), line[3]->index()); - Quad quad_compare_2(line[3]->index(), line[2]->index(), line[1]->index(), line[0]->index()); + Quad quad_compare_2(line[2]->index(), line[3]->index(), line[0]->index(), line[1]->index()); // try to find the quad with // lines situated as @@ -1367,19 +1497,35 @@ Triangulation<3>::create_triangulation (const std::vector > &v, // the boundary lines 3 times // until it is found or it does // not exist. + + // mapping from counterclock to + // lexicographic ordering of + // quad lines + static const unsigned int lex2cclock[4]={3,1,0,2}; + // copy lines from + // lexicographic to + // counterclock ordering, as + // rotation is much simpler in + // counterclock ordering + line_iterator line_counterclock[4]; + for (unsigned int i=0; i<4; ++i) + line_counterclock[lex2cclock[i]]=line[i]; unsigned int n_rotations=0; bool not_found_quad_1; while ( (not_found_quad_1=(needed_quads.find(quad_compare_1) == needed_quads.end())) && ( needed_quads.find(quad_compare_2) == needed_quads.end()) && (n_rotations<4)) { - rotate(line, line+1, line+4); + // use the rotate defined + // in + rotate(line_counterclock, line_counterclock+1, line_counterclock+4); // update the quads with - // rotated lines + // rotated lines (i runs in + // lexicographic ordering) for (unsigned int i=0; i<4; ++i) { - quad_compare_1.set_line(i, line[i]->index()); - quad_compare_2.set_line(3-i, line[i]->index()); + quad_compare_1.set_line(i, line_counterclock[lex2cclock[i]]->index()); + quad_compare_2.set_line((i+2)%4, line_counterclock[lex2cclock[i]]->index()); } ++n_rotations; @@ -4730,33 +4876,33 @@ Triangulation<2>::execute_refinement () Set up an array of the 3x3 vertices, which are distributed on the cell (the array consists of indices into the @p{vertices} std::vector - 6--5--4 + 2--7--3 | | | - 7--8--3 + 4--8--5 | | | - 0--1--2 + 0--6--1 Second: Set up an array of the new lines (the array consists of iterator pointers into the lines arrays) - .-4-.-5-. The directions are: .->-.->-. - 7 9 3 ^ ^ ^ - .-10.11- . .->-.->-. - 6 8 2 ^ ^ ^ - .-0-.-1-. .->-.->-. + .-6-.-7-. The directions are: .->-.->-. + 1 9 3 ^ ^ ^ + .-10.11-. .->-.->-. + 0 8 2 ^ ^ ^ + .-4-.-5-. .->-.->-. Third: Set up an array of neighbors: - 4 5 + 6 7 .--.--. - 7| | |3 + 1| | |3 .--.--. - 6| | |2 + 0| | |2 .--.--. - 0 1 + 4 5 We need this array for two reasons: first to get the lines which will bound the four subcells (if the neighboring cell is refined, these @@ -4783,7 +4929,7 @@ Triangulation<2>::execute_refinement () The created children are numbered like this: .--.--. - |3 . 2| + |2 . 3| .--.--. |0 | 1| .--.--. @@ -4801,12 +4947,12 @@ Triangulation<2>::execute_refinement () // collect the indices // all vertices int new_vertices[9] = {cell->vertex_index(0), - cell->line(0)->child(0)->vertex_index(1), cell->vertex_index(1), - cell->line(1)->child(0)->vertex_index(1), cell->vertex_index(2), - cell->line(2)->child(0)->vertex_index(1), cell->vertex_index(3), + cell->line(0)->child(0)->vertex_index(1), + cell->line(1)->child(0)->vertex_index(1), + cell->line(2)->child(0)->vertex_index(1), cell->line(3)->child(0)->vertex_index(1), next_unused_vertex}; @@ -4847,10 +4993,10 @@ Triangulation<2>::execute_refinement () ExcCellShouldBeUnused()); } - new_lines[8] ->set (Line(new_vertices[1], new_vertices[8])); - new_lines[9] ->set (Line(new_vertices[8], new_vertices[5])); - new_lines[10]->set (Line(new_vertices[7], new_vertices[8])); - new_lines[11]->set (Line(new_vertices[8], new_vertices[3])); + new_lines[8] ->set (Line(new_vertices[6], new_vertices[8])); + new_lines[9] ->set (Line(new_vertices[8], new_vertices[7])); + new_lines[10]->set (Line(new_vertices[4], new_vertices[8])); + new_lines[11]->set (Line(new_vertices[8], new_vertices[5])); for (l=8; l<12; ++l) { @@ -4889,19 +5035,19 @@ Triangulation<2>::execute_refinement () Assert(n_children==4, ExcNotImplemented()); subcells[0]->set (Quad(new_lines[0]->index(), new_lines[8]->index(), - new_lines[10]->index(), - new_lines[6]->index())); - subcells[1]->set (Quad(new_lines[1]->index(), + new_lines[4]->index(), + new_lines[10]->index())); + subcells[1]->set (Quad(new_lines[8]->index(), new_lines[2]->index(), - new_lines[11]->index(), - new_lines[8]->index())); - subcells[2]->set (Quad(new_lines[11]->index(), - new_lines[3]->index(), new_lines[5]->index(), - new_lines[9]->index())); - subcells[3]->set (Quad(new_lines[10]->index(), + new_lines[11]->index())); + subcells[2]->set (Quad(new_lines[1]->index(), new_lines[9]->index(), - new_lines[4]->index(), + new_lines[10]->index(), + new_lines[6]->index())); + subcells[3]->set (Quad(new_lines[9]->index(), + new_lines[3]->index(), + new_lines[11]->index(), new_lines[7]->index())); for (unsigned int i=0; i::execute_refinement () // in the new cells subcells[0]->set_neighbor (0, neighbor_cells[0][0]); subcells[0]->set_neighbor (1, subcells[1]); - subcells[0]->set_neighbor (2, subcells[3]); - subcells[0]->set_neighbor (3, neighbor_cells[3][0]); + subcells[0]->set_neighbor (2, neighbor_cells[2][0]); + subcells[0]->set_neighbor (3, subcells[2]); - subcells[1]->set_neighbor (0, neighbor_cells[0][1]); + subcells[1]->set_neighbor (0, subcells[0]); subcells[1]->set_neighbor (1, neighbor_cells[1][0]); - subcells[1]->set_neighbor (2, subcells[2]); - subcells[1]->set_neighbor (3, subcells[0]); + subcells[1]->set_neighbor (2, neighbor_cells[2][1]); + subcells[1]->set_neighbor (3, subcells[3]); - subcells[2]->set_neighbor (0, subcells[1]); - subcells[2]->set_neighbor (1, neighbor_cells[1][1]); - subcells[2]->set_neighbor (2, neighbor_cells[2][1]); - subcells[2]->set_neighbor (3, subcells[3]); + subcells[2]->set_neighbor (0, neighbor_cells[0][1]); + subcells[2]->set_neighbor (1, subcells[3]); + subcells[2]->set_neighbor (2, subcells[0]); + subcells[2]->set_neighbor (3, neighbor_cells[3][0]); - subcells[3]->set_neighbor (0, subcells[0]); - subcells[3]->set_neighbor (1, subcells[2]); - subcells[3]->set_neighbor (2, neighbor_cells[2][0]); + subcells[3]->set_neighbor (0, subcells[2]); + subcells[3]->set_neighbor (1, neighbor_cells[1][1]); + subcells[3]->set_neighbor (2, subcells[1]); subcells[3]->set_neighbor (3, neighbor_cells[3][1]); // now we need to set the @@ -5544,16 +5690,16 @@ Triangulation<3>::execute_refinement () // first collect the // indices of the five // vertices: - // *--2--* + // *--3--* // | | | - // 3--4--1 + // 0--4--1 // | | | - // *--0--* + // *--2--* // the lines are numbered // as follows: // *--*--* - // | 2 | - // *3-*-1* + // | 1 | + // *2-*-3* // | 0 | // *--*--* const unsigned int vertex_indices[5] @@ -5564,10 +5710,10 @@ Triangulation<3>::execute_refinement () next_unused_vertex }; - new_lines[0]->set (Line(vertex_indices[0], vertex_indices[4])); - new_lines[1]->set (Line(vertex_indices[4], vertex_indices[1])); - new_lines[2]->set (Line(vertex_indices[4], vertex_indices[2])); - new_lines[3]->set (Line(vertex_indices[3], vertex_indices[4])); + new_lines[0]->set (Line(vertex_indices[2], vertex_indices[4])); + new_lines[1]->set (Line(vertex_indices[4], vertex_indices[3])); + new_lines[2]->set (Line(vertex_indices[0], vertex_indices[4])); + new_lines[3]->set (Line(vertex_indices[4], vertex_indices[1])); for (unsigned int i=0; i<4; ++i) { @@ -5585,20 +5731,20 @@ Triangulation<3>::execute_refinement () // about the indices of // the lines, with the // following numbering: - // *5-*-4* - // 6 10 3 - // *11*-9* - // 7 8 2 - // *0-*-1* + // .-6-.-7-. + // 1 9 3 + // .-10.11-. + // 0 8 2 + // .-4-.-5-. const unsigned int line_indices[12] = { quad->line(0)->child(0)->index(), quad->line(0)->child(1)->index(), quad->line(1)->child(0)->index(), quad->line(1)->child(1)->index(), - quad->line(2)->child(1)->index(), quad->line(2)->child(0)->index(), - quad->line(3)->child(1)->index(), + quad->line(2)->child(1)->index(), quad->line(3)->child(0)->index(), + quad->line(3)->child(1)->index(), new_lines[0]->index(), new_lines[1]->index(), new_lines[2]->index(), @@ -5639,20 +5785,20 @@ Triangulation<3>::execute_refinement () new_quads[0]->set (Quad(line_indices[0], line_indices[8], - line_indices[11], - line_indices[7])); - new_quads[1]->set (Quad(line_indices[1], - line_indices[2], - line_indices[9], - line_indices[8])); - new_quads[2]->set (Quad(line_indices[9], - line_indices[3], line_indices[4], line_indices[10])); - new_quads[3]->set (Quad(line_indices[11], - line_indices[10], + new_quads[1]->set (Quad(line_indices[8], + line_indices[2], line_indices[5], + line_indices[11])); + new_quads[2]->set (Quad(line_indices[1], + line_indices[9], + line_indices[10], line_indices[6])); + new_quads[3]->set (Quad(line_indices[9], + line_indices[3], + line_indices[11], + line_indices[7])); for (unsigned int i=0; i<4; ++i) { new_quads[i]->set_used_flag(); @@ -5756,7 +5902,7 @@ Triangulation<3>::execute_refinement () // orientation for (unsigned int face=0; face::faces_per_cell; ++face) - vertices[next_unused_vertex] += hex->face(face)->child(0)->vertex(2) * + vertices[next_unused_vertex] += hex->face(face)->child(0)->vertex(3) * 1./12.; // now that we created @@ -5788,40 +5934,40 @@ Triangulation<3>::execute_refinement () // cutting the hex in two // vertically and // horizontally) - // *--2--* *--5--* + // *--3--* *--5--* // / / / | | | - // 3--6--1 3--6--1 + // 0--6--1 0--6--1 // / / / | | | - // *--0--* *--4--* + // *--2--* *--4--* // the lines are numbered // as follows: // *--*--* *--*--* - // / 2 / | 5 | - // *3-*-1* *3-*-1* + // / 1 / | 5 | + // *2-*-3* *2-*-3* // / 0 / | 4 | // *--*--* *--*--* // // note that both asking // for child 0 and for - // vertex 2 within that + // vertex 3 within that // is invariant with // respect to the face // orientation, so we do // not have to ask here const unsigned int vertex_indices[7] - = { hex->face(0)->child(0)->vertex_index(2), - hex->face(3)->child(0)->vertex_index(2), - hex->face(1)->child(0)->vertex_index(2), - hex->face(5)->child(0)->vertex_index(2), - hex->face(2)->child(0)->vertex_index(2), - hex->face(4)->child(0)->vertex_index(2), + = { hex->face(0)->child(0)->vertex_index(3), + hex->face(1)->child(0)->vertex_index(3), + hex->face(2)->child(0)->vertex_index(3), + hex->face(3)->child(0)->vertex_index(3), + hex->face(4)->child(0)->vertex_index(3), + hex->face(5)->child(0)->vertex_index(3), next_unused_vertex }; - new_lines[0]->set (Line(vertex_indices[0], vertex_indices[6])); - new_lines[1]->set (Line(vertex_indices[6], vertex_indices[1])); - new_lines[2]->set (Line(vertex_indices[6], vertex_indices[2])); - new_lines[3]->set (Line(vertex_indices[3], vertex_indices[6])); + new_lines[0]->set (Line(vertex_indices[2], vertex_indices[6])); + new_lines[1]->set (Line(vertex_indices[6], vertex_indices[3])); + new_lines[2]->set (Line(vertex_indices[0], vertex_indices[6])); + new_lines[3]->set (Line(vertex_indices[6], vertex_indices[1])); new_lines[4]->set (Line(vertex_indices[4], vertex_indices[6])); new_lines[5]->set (Line(vertex_indices[6], vertex_indices[5])); @@ -5841,71 +5987,68 @@ Triangulation<3>::execute_refinement () // collect some data // about the indices of // the lines, with the - // following numbering: - // front plane *---*---* - // | 2 | - // *3--*--1* - // | 0 | - // *---*---* - // - // middle plane *9--*--8* - // 10 14 7 - // *15-*-13* - // 11 12 6 - // *4--*--5* - // - // back plane *---*---* - // | 18 | - // *19-*-17* - // | 16 | - // *---*---* - // - // left plane (the - // left-to-right planes - // are displayed twice - // each, for better - // readability of the - // indices; the left part - // is already determined - // by the pictures above) - // * * - // /| /| - // * | * | - // /| * /| * - // *10/| * |21 - // | * | | * | - // |/| * |20 * - // *11/ * |/ - // | * | * - // |/ |/ - // * * - // - // middle plane - // * * - // /| 23| - // * 18 * | - // /| * 22| * - // *14/16 * |25 - // | * | | * | - // 2/| * |24 * - // *12/ * |27 - // | * | * - // 0/ |26 - // * * - // - // - // right plane - // * * - // /| /| - // * | * | - // /| * /| * - // * 6/| * |29 - // | * | | * | - // |/| * |28 * - // * 7/ * |/ - // | * | * - // |/ |/ - // * * + // following numbering: + // (note that face 0 and + // 1 each are shown twice + // for better + // readability) + + // face 0: left plane + // * * + // /| /| + // * | * | + // /| * /| * + // * 1/| * |3| + // | * | | * | + // |/| * |2| * + // * 0/ * |/ + // | * | * + // |/ |/ + // * * + // face 1: right plane + // * * + // /| /| + // * | * | + // /| * /| * + // * 5/| * |7| + // | * | | * | + // |/| * |6| * + // * 4/ * |/ + // | * | * + // |/ |/ + // * * + // face 2: front plane + // (note: x,y exchanged) + // *---*---* + // | 11 | + // *-8-*-9-* + // | 10 | + // *---*---* + // face 3: back plane + // (note: x,y exchanged) + // *---*---* + // | 15 | + // *12-*-13* + // | 14 | + // *---*---* + // face 4: bottom plane + // *---*---* + // / 17 / + // *18-*-19* + // / 16 / + // *---*---* + // face 5: top plane + // *---*---* + // / 21 / + // *22-*-23* + // / 20 / + // *---*---* + // middle planes + // *---*---* *---*---* + // / 25 / | 29 | + // *26-*-27* *26-*-27* + // / 24 / | 28 | + // *---*---* *---*---* // // this time we have to // take into account @@ -5914,6 +6057,8 @@ Triangulation<3>::execute_refinement () // correctly or in the // opposite direction, so // store that up front + + // TODO: shorten this const bool face_orientation[6] = { hex->face_orientation (0), hex->face_orientation (1), @@ -5924,67 +6069,71 @@ Triangulation<3>::execute_refinement () const unsigned int line_indices[30] = { - hex->face(0)->child(0 ) - ->line_index(face_orientation[0] ? 1 : 2), //0 - hex->face(0)->child(face_orientation[0] ? 1 : 3) - ->line_index(face_orientation[0] ? 2 : 1), //1 - hex->face(0)->child(2 ) - ->line_index(face_orientation[0] ? 3 : 0), //2 - hex->face(0)->child(face_orientation[0] ? 3 : 1) - ->line_index(face_orientation[0] ? 0 : 3), //3 - - hex->face(2)->child(0 ) - ->line_index(face_orientation[2] ? 2 : 1), //4 - hex->face(2)->child(face_orientation[2] ? 1 : 3) - ->line_index(face_orientation[2] ? 2 : 1), //5 - hex->face(3)->child(0 ) - ->line_index(face_orientation[3] ? 1 : 2), //6 - hex->face(3)->child(face_orientation[3] ? 3 : 1) - ->line_index(face_orientation[3] ? 1 : 2), //7 - - hex->face(4)->child(face_orientation[4] ? 1 : 3) - ->line_index(face_orientation[4] ? 2 : 1), //8 - hex->face(4)->child(0 ) - ->line_index(face_orientation[4] ? 2 : 1), //9 - hex->face(5)->child(face_orientation[5] ? 3 : 1) - ->line_index(face_orientation[5] ? 1 : 2), //10 - hex->face(5)->child(0 ) - ->line_index(face_orientation[5] ? 1 : 2), //11 - - new_lines[4]->index(), //12 - new_lines[1]->index(), //13 - new_lines[5]->index(), //14 - new_lines[3]->index(), //15 - - hex->face(1)->child(0 ) - ->line_index(face_orientation[1] ? 1 : 2), //16 - hex->face(1)->child(face_orientation[1] ? 1 : 3) - ->line_index(face_orientation[1] ? 2 : 1), //17 - hex->face(1)->child(2 ) - ->line_index(face_orientation[1] ? 3 : 0), //18 - hex->face(1)->child(face_orientation[1] ? 3 : 1) - ->line_index(face_orientation[1] ? 0 : 3), //19 - - hex->face(5)->child(0 ) - ->line_index(face_orientation[5] ? 2 : 1), //20 - hex->face(5)->child(face_orientation[5] ? 1 : 3) - ->line_index(face_orientation[5] ? 2 : 1), //21 - hex->face(4)->child(0 ) - ->line_index(face_orientation[4] ? 1 : 2), //22 - hex->face(4)->child(face_orientation[4] ? 3 : 1) - ->line_index(face_orientation[4] ? 1 : 2), //23 - - new_lines[0]->index(), //24 - new_lines[2]->index(), //25 - hex->face(2)->child(0 ) - ->line_index(face_orientation[2] ? 1 : 2), //26 - hex->face(2)->child(face_orientation[2] ? 3 : 1) - ->line_index(face_orientation[2] ? 1 : 2), //27 - - hex->face(3)->child(0 ) - ->line_index(face_orientation[3] ? 2 : 1), //28 - hex->face(3)->child(face_orientation[3] ? 1 : 3) - ->line_index(face_orientation[3] ? 2 : 1) //29 + // ask for + // children 0 and + // 3 as they are + // invariant + // wrt. face_orientation + hex->face(0)->child(0) + ->line_index(face_orientation[0] ? 1 : 3), //0 + hex->face(0)->child(3) + ->line_index(face_orientation[0] ? 0 : 2), //1 + hex->face(0)->child(0) + ->line_index(face_orientation[0] ? 3 : 1), //2 + hex->face(0)->child(3) + ->line_index(face_orientation[0] ? 2 : 0), //3 + + hex->face(1)->child(0) + ->line_index(face_orientation[1] ? 1 : 3), //4 + hex->face(1)->child(3) + ->line_index(face_orientation[1] ? 0 : 2), //5 + hex->face(1)->child(0) + ->line_index(face_orientation[1] ? 3 : 1), //6 + hex->face(1)->child(3) + ->line_index(face_orientation[1] ? 2 : 0), //7 + + hex->face(2)->child(0) + ->line_index(face_orientation[2] ? 1 : 3), //8 + hex->face(2)->child(3) + ->line_index(face_orientation[2] ? 0 : 2), //9 + hex->face(2)->child(0) + ->line_index(face_orientation[2] ? 3 : 1), //10 + hex->face(2)->child(3) + ->line_index(face_orientation[2] ? 2 : 0), //11 + + hex->face(3)->child(0) + ->line_index(face_orientation[3] ? 1 : 3), //12 + hex->face(3)->child(3) + ->line_index(face_orientation[3] ? 0 : 2), //13 + hex->face(3)->child(0) + ->line_index(face_orientation[3] ? 3 : 1), //14 + hex->face(3)->child(3) + ->line_index(face_orientation[3] ? 2 : 0), //15 + + hex->face(4)->child(0) + ->line_index(face_orientation[4] ? 1 : 3), //16 + hex->face(4)->child(3) + ->line_index(face_orientation[4] ? 0 : 2), //17 + hex->face(4)->child(0) + ->line_index(face_orientation[4] ? 3 : 1), //18 + hex->face(4)->child(3) + ->line_index(face_orientation[4] ? 2 : 0), //19 + + hex->face(5)->child(0) + ->line_index(face_orientation[5] ? 1 : 3), //20 + hex->face(5)->child(3) + ->line_index(face_orientation[5] ? 0 : 2), //21 + hex->face(5)->child(0) + ->line_index(face_orientation[5] ? 3 : 1), //22 + hex->face(5)->child(3) + ->line_index(face_orientation[5] ? 2 : 0), //23 + + new_lines[0]->index(), //24 + new_lines[1]->index(), //25 + new_lines[2]->index(), //26 + new_lines[3]->index(), //27 + new_lines[4]->index(), //28 + new_lines[5]->index() //29 }; // find some space for @@ -6005,88 +6154,83 @@ Triangulation<3>::execute_refinement () // set up the 12 quads, // numbered as follows - // (shown are the three - // planes cutting the hex - // in two): - // - // *-----*-----* - // | 3 | 2 | - // | | | - // *-----*-----* - // | | | - // | 0 | 1 | - // *-----*-----* + // (left quad numbering, + // right line numbering + // extracted from above) // - // *----*----* - // / 7 / 6 / - // *----*----* - // / 4 / 5 / - // *----*----* + // * * + // /| 21| + // * | * 15 + // y/|3* 20| * + // * |/| * |/| + // |2* |x 11 * 14 + // |/|1* |/| * + // * |/ * |17 + // |0* 10 * + // |/ |16 + // * * // + // x + // *---*---* *22-*-23* + // | 5 | 7 | 1 29 5 + // *---*---* *26-*-27* + // | 4 | 6 | 0 28 4 + // *---*---*y *18-*-19* // - // * - // /| - // / | - // *10| - // /| * - // / | /| - // * |/ | - // |11* | - // | /| 9* - // |/ | / - // * |/ - // |8 * - // | / - // |/ - // * - new_quads[0]->set (Quad(line_indices[4], - line_indices[12], - line_indices[15], - line_indices[11])); - new_quads[1]->set (Quad(line_indices[5], - line_indices[6], - line_indices[13], - line_indices[12])); - new_quads[2]->set (Quad(line_indices[13], - line_indices[7], - line_indices[8], - line_indices[14])); - new_quads[3]->set (Quad(line_indices[15], - line_indices[14], - line_indices[9], - line_indices[10])); - new_quads[4]->set (Quad(line_indices[3], - line_indices[24], - line_indices[15], - line_indices[20])); - new_quads[5]->set (Quad(line_indices[1], + // y + // *----*----* *-12-*-13-* + // / 10 / 11 / 3 25 7 + // *----*----* *-26-*-27-* + // / 8 / 9 / 2 24 6 + // *----*----*x *--8-*--9-* + new_quads[0]->set (Quad(line_indices[10], line_indices[28], - line_indices[13], + line_indices[16], line_indices[24])); - new_quads[6]->set (Quad(line_indices[13], - line_indices[29], + new_quads[1]->set (Quad(line_indices[28], + line_indices[14], line_indices[17], line_indices[25])); - new_quads[7]->set (Quad(line_indices[15], + new_quads[2]->set (Quad(line_indices[11], + line_indices[29], + line_indices[24], + line_indices[20])); + new_quads[3]->set (Quad(line_indices[29], + line_indices[15], line_indices[25], - line_indices[19], line_indices[21])); - new_quads[8]->set (Quad(line_indices[26], - line_indices[12], + new_quads[4]->set (Quad(line_indices[18], + line_indices[26], + line_indices[0], + line_indices[28])); + new_quads[5]->set (Quad(line_indices[26], + line_indices[22], + line_indices[1], + line_indices[29])); + new_quads[6]->set (Quad(line_indices[19], + line_indices[27], + line_indices[28], + line_indices[4])); + new_quads[7]->set (Quad(line_indices[27], + line_indices[23], + line_indices[29], + line_indices[5])); + new_quads[8]->set (Quad(line_indices[2], line_indices[24], - line_indices[0])); - new_quads[9]->set (Quad(line_indices[27], - line_indices[16], - line_indices[25], - line_indices[12])); - new_quads[10]->set (Quad(line_indices[25], - line_indices[18], - line_indices[23], - line_indices[14])); - new_quads[11]->set (Quad(line_indices[24], - line_indices[14], - line_indices[22], - line_indices[2])); + line_indices[8], + line_indices[26])); + new_quads[9]->set (Quad(line_indices[24], + line_indices[6], + line_indices[9], + line_indices[27])); + new_quads[10]->set (Quad(line_indices[3], + line_indices[25], + line_indices[26], + line_indices[12])); + new_quads[11]->set (Quad(line_indices[25], + line_indices[7], + line_indices[27], + line_indices[13])); for (unsigned int i=0; i<12; ++i) { new_quads[i]->set_used_flag(); @@ -6110,48 +6254,29 @@ Triangulation<3>::execute_refinement () // // planes in the interior // of the old hex: - // *-----*-----* - // | 3 | 2 | - // | | | - // *-----*-----* - // | | | - // | 0 | 1 | - // *-----*-----* - // - // *----*----* - // / 7 / 6 / - // *----*----* - // / 4 / 5 / - // *----*----* - // - // - // * - // /| - // / | - // *10| - // /| * - // / | /| - // * |/ | - // |11* | - // | /| 9* - // |/ | / - // * |/ - // |8 * - // | / - // |/ - // * + // * + // /| + // * | + // /|3* *---*---* *----*----* + // * |/| | 5 | 7 | / 10 / 11 / + // |2* | *---*---* *----*----* + // |/|1* | 4 | 6 | / 8 / 9 / + // * |/ *---*---*y *----*----*x + // |0* + // |/ + // * // // children of the faces // of the old hex // *-------* *-------* - // /|19 18| /31 30/| - // 34| | / /26 - // / | | /28 29/ | - // * |16 17| *-------*27 | - // 3533*-------* |15 14| 25* - // | /23 22/ | | / - // 32/ / | |24 - // |/20 21/ |12 13|/ + // /|25 27| /34 35/| + // 15| | / /19 + // / | | /32 33/ | + // * |24 26| *-------*18 | + // 1413*-------* |21 23| 17* + // | /30 31/ | | / + // 12/ / | |16 + // |/28 29/ |20 22|/ // *-------* *-------* // // note that we have to @@ -6159,11 +6284,13 @@ Triangulation<3>::execute_refinement () // orientation of // faces. as an // optimization: asking - // for child 0 or 2 of a + // for child 0 or 3 of a // face is invariant // under the orientation, // so we don't have to // ask for it then + // + // TODO: simplify this const unsigned int quad_indices[36] = { new_quads[0]->index(), //0 @@ -6180,34 +6307,34 @@ Triangulation<3>::execute_refinement () new_quads[11]->index(), //11 hex->face(0)->child_index(0), //12 - hex->face(0)->child_index(face_orientation[0] ? 1 : 3), - hex->face(0)->child_index(2), - hex->face(0)->child_index(face_orientation[0] ? 3 : 1), + hex->face(0)->child_index(face_orientation[0] ? 1 : 2), + hex->face(0)->child_index(face_orientation[0] ? 2 : 1), + hex->face(0)->child_index(3), hex->face(1)->child_index(0), //16 - hex->face(1)->child_index(face_orientation[1] ? 1 : 3), - hex->face(1)->child_index(2), - hex->face(1)->child_index(face_orientation[1] ? 3 : 1), + hex->face(1)->child_index(face_orientation[1] ? 1 : 2), + hex->face(1)->child_index(face_orientation[1] ? 2 : 1), + hex->face(1)->child_index(3), hex->face(2)->child_index(0), //20 - hex->face(2)->child_index(face_orientation[2] ? 1 : 3), - hex->face(2)->child_index(2), - hex->face(2)->child_index(face_orientation[2] ? 3 : 1), + hex->face(2)->child_index(face_orientation[2] ? 1 : 2), + hex->face(2)->child_index(face_orientation[2] ? 2 : 1), + hex->face(2)->child_index(3), hex->face(3)->child_index(0), //24 - hex->face(3)->child_index(face_orientation[3] ? 1 : 3), - hex->face(3)->child_index(2), - hex->face(3)->child_index(face_orientation[3] ? 3 : 1), + hex->face(3)->child_index(face_orientation[3] ? 1 : 2), + hex->face(3)->child_index(face_orientation[3] ? 2 : 1), + hex->face(3)->child_index(3), hex->face(4)->child_index(0), //28 - hex->face(4)->child_index(face_orientation[4] ? 1 : 3), - hex->face(4)->child_index(2), - hex->face(4)->child_index(face_orientation[4] ? 3 : 1), + hex->face(4)->child_index(face_orientation[4] ? 1 : 2), + hex->face(4)->child_index(face_orientation[4] ? 2 : 1), + hex->face(4)->child_index(3), hex->face(5)->child_index(0), //32 - hex->face(5)->child_index(face_orientation[5] ? 1 : 3), - hex->face(5)->child_index(2), - hex->face(5)->child_index(face_orientation[5] ? 3 : 1) + hex->face(5)->child_index(face_orientation[5] ? 1 : 2), + hex->face(5)->child_index(face_orientation[5] ? 2 : 1), + hex->face(5)->child_index(3) }; @@ -6236,57 +6363,57 @@ Triangulation<3>::execute_refinement () // present cell hex->set_children (new_hexes[0]->index()); - // front children + // bottom children new_hexes[0]->set (Hexahedron(quad_indices[12], quad_indices[0], quad_indices[20], - quad_indices[8], quad_indices[4], - quad_indices[32])); - new_hexes[1]->set (Hexahedron(quad_indices[13], - quad_indices[1], - quad_indices[21], - quad_indices[24], - quad_indices[5], + quad_indices[28], quad_indices[8])); - new_hexes[2]->set (Hexahedron(quad_indices[14], - quad_indices[2], - quad_indices[5], - quad_indices[27], + new_hexes[1]->set (Hexahedron(quad_indices[0], + quad_indices[16], + quad_indices[22], + quad_indices[6], quad_indices[29], - quad_indices[11])); - new_hexes[3]->set (Hexahedron(quad_indices[15], - quad_indices[3], + quad_indices[9])); + new_hexes[2]->set (Hexahedron(quad_indices[13], + quad_indices[1], quad_indices[4], - quad_indices[11], - quad_indices[28], - quad_indices[35])); + quad_indices[24], + quad_indices[30], + quad_indices[10])); + new_hexes[3]->set (Hexahedron(quad_indices[1], + quad_indices[17], + quad_indices[6], + quad_indices[26], + quad_indices[31], + quad_indices[11])); - // back children - new_hexes[4]->set (Hexahedron(quad_indices[0], - quad_indices[16], + // top children + new_hexes[4]->set (Hexahedron(quad_indices[14], + quad_indices[2], + quad_indices[21], + quad_indices[5], + quad_indices[8], + quad_indices[32])); + new_hexes[5]->set (Hexahedron(quad_indices[2], + quad_indices[18], quad_indices[23], - quad_indices[9], quad_indices[7], + quad_indices[9], quad_indices[33])); - new_hexes[5]->set (Hexahedron(quad_indices[1], - quad_indices[17], - quad_indices[22], + new_hexes[6]->set (Hexahedron(quad_indices[15], + quad_indices[3], + quad_indices[5], quad_indices[25], - quad_indices[6], - quad_indices[9])); - new_hexes[6]->set (Hexahedron(quad_indices[2], - quad_indices[18], - quad_indices[6], - quad_indices[26], - quad_indices[30], - quad_indices[10])); + quad_indices[10], + quad_indices[34])); new_hexes[7]->set (Hexahedron(quad_indices[3], quad_indices[19], quad_indices[7], - quad_indices[10], - quad_indices[31], - quad_indices[34])); + quad_indices[27], + quad_indices[11], + quad_indices[35])); for (unsigned int i=0; i<8; ++i) @@ -6509,61 +6636,114 @@ Triangulation<3>::execute_refinement () // now we've got all // neighbors, so set them // in the new cells + + // the subface ordering + // is as follows + // face 0 and 1: + // * + // /| + // * | + // y/|3* + // * |/| + // |2* |x + // |/|1* + // * |/ + // |0* + // |/ + // * + // face 2 and 3: + // x + // *---*---* + // | 1 | 3 | + // *---*---* + // | 0 | 2 | + // *---*---*y + // face 4 and 5: + // y + // *---*---* + // / 2 / 3 / + // *---*---* + // / 0 / 1 / + // *---*---*x new_hexes[0]->set_neighbor (0, neighbor_cells[0][0]); - new_hexes[0]->set_neighbor (1, new_hexes[4]); + new_hexes[0]->set_neighbor (1, new_hexes[1]); new_hexes[0]->set_neighbor (2, neighbor_cells[2][0]); - new_hexes[0]->set_neighbor (3, new_hexes[1]); - new_hexes[0]->set_neighbor (4, new_hexes[3]); - new_hexes[0]->set_neighbor (5, neighbor_cells[5][0]); - - new_hexes[1]->set_neighbor (0, neighbor_cells[0][1]); - new_hexes[1]->set_neighbor (1, new_hexes[5]); - new_hexes[1]->set_neighbor (2, neighbor_cells[2][1]); - new_hexes[1]->set_neighbor (3, neighbor_cells[3][0]); - new_hexes[1]->set_neighbor (4, new_hexes[2]); - new_hexes[1]->set_neighbor (5, new_hexes[0]); - - new_hexes[2]->set_neighbor (0, neighbor_cells[0][2]); - new_hexes[2]->set_neighbor (1, new_hexes[6]); - new_hexes[2]->set_neighbor (2, new_hexes[1]); - new_hexes[2]->set_neighbor (3, neighbor_cells[3][3]); - new_hexes[2]->set_neighbor (4, neighbor_cells[4][1]); - new_hexes[2]->set_neighbor (5, new_hexes[3]); - - new_hexes[3]->set_neighbor (0, neighbor_cells[0][3]); - new_hexes[3]->set_neighbor (1, new_hexes[7]); - new_hexes[3]->set_neighbor (2, new_hexes[0]); - new_hexes[3]->set_neighbor (3, new_hexes[2]); - new_hexes[3]->set_neighbor (4, neighbor_cells[4][0]); - new_hexes[3]->set_neighbor (5, neighbor_cells[5][3]); - - new_hexes[4]->set_neighbor (0, new_hexes[0]); - new_hexes[4]->set_neighbor (1, neighbor_cells[1][0]); - new_hexes[4]->set_neighbor (2, neighbor_cells[2][3]); - new_hexes[4]->set_neighbor (3, new_hexes[5]); - new_hexes[4]->set_neighbor (4, new_hexes[7]); - new_hexes[4]->set_neighbor (5, neighbor_cells[5][1]); - - new_hexes[5]->set_neighbor (0, new_hexes[1]); - new_hexes[5]->set_neighbor (1, neighbor_cells[1][1]); - new_hexes[5]->set_neighbor (2, neighbor_cells[2][2]); - new_hexes[5]->set_neighbor (3, neighbor_cells[3][1]); - new_hexes[5]->set_neighbor (4, new_hexes[6]); - new_hexes[5]->set_neighbor (5, new_hexes[4]); - - new_hexes[6]->set_neighbor (0, new_hexes[2]); - new_hexes[6]->set_neighbor (1, neighbor_cells[1][2]); - new_hexes[6]->set_neighbor (2, new_hexes[5]); - new_hexes[6]->set_neighbor (3, neighbor_cells[3][2]); - new_hexes[6]->set_neighbor (4, neighbor_cells[4][2]); - new_hexes[6]->set_neighbor (5, new_hexes[7]); - - new_hexes[7]->set_neighbor (0, new_hexes[3]); + new_hexes[0]->set_neighbor (3, new_hexes[2]); + new_hexes[0]->set_neighbor (4, neighbor_cells[4][0]); + new_hexes[0]->set_neighbor (5, new_hexes[4]); + Assert(GeometryInfo::child_cell_on_face(0,0)==0, ExcInternalError()); + Assert(GeometryInfo::child_cell_on_face(2,0)==0, ExcInternalError()); + Assert(GeometryInfo::child_cell_on_face(4,0)==0, ExcInternalError()); + + new_hexes[1]->set_neighbor (0, new_hexes[0]); + new_hexes[1]->set_neighbor (1, neighbor_cells[1][0]); + new_hexes[1]->set_neighbor (2, neighbor_cells[2][2]); + new_hexes[1]->set_neighbor (3, new_hexes[3]); + new_hexes[1]->set_neighbor (4, neighbor_cells[4][1]); + new_hexes[1]->set_neighbor (5, new_hexes[5]); + Assert(GeometryInfo::child_cell_on_face(1,0)==1, ExcInternalError()); + Assert(GeometryInfo::child_cell_on_face(2,2)==1, ExcInternalError()); + Assert(GeometryInfo::child_cell_on_face(4,1)==1, ExcInternalError()); + + new_hexes[2]->set_neighbor (0, neighbor_cells[0][1]); + new_hexes[2]->set_neighbor (1, new_hexes[3]); + new_hexes[2]->set_neighbor (2, new_hexes[0]); + new_hexes[2]->set_neighbor (3, neighbor_cells[3][0]); + new_hexes[2]->set_neighbor (4, neighbor_cells[4][2]); + new_hexes[2]->set_neighbor (5, new_hexes[6]); + Assert(GeometryInfo::child_cell_on_face(0,1)==2, ExcInternalError()); + Assert(GeometryInfo::child_cell_on_face(3,0)==2, ExcInternalError()); + Assert(GeometryInfo::child_cell_on_face(4,2)==2, ExcInternalError()); + + new_hexes[3]->set_neighbor (0, new_hexes[2]); + new_hexes[3]->set_neighbor (1, neighbor_cells[1][1]); + new_hexes[3]->set_neighbor (2, new_hexes[1]); + new_hexes[3]->set_neighbor (3, neighbor_cells[3][2]); + new_hexes[3]->set_neighbor (4, neighbor_cells[4][3]); + new_hexes[3]->set_neighbor (5, new_hexes[7]); + Assert(GeometryInfo::child_cell_on_face(1,1)==3, ExcInternalError()); + Assert(GeometryInfo::child_cell_on_face(3,2)==3, ExcInternalError()); + Assert(GeometryInfo::child_cell_on_face(4,3)==3, ExcInternalError()); + + new_hexes[4]->set_neighbor (0, neighbor_cells[0][2]); + new_hexes[4]->set_neighbor (1, new_hexes[5]); + new_hexes[4]->set_neighbor (2, neighbor_cells[2][1]); + new_hexes[4]->set_neighbor (3, new_hexes[6]); + new_hexes[4]->set_neighbor (4, new_hexes[0]); + new_hexes[4]->set_neighbor (5, neighbor_cells[5][0]); + Assert(GeometryInfo::child_cell_on_face(0,2)==4, ExcInternalError()); + Assert(GeometryInfo::child_cell_on_face(2,1)==4, ExcInternalError()); + Assert(GeometryInfo::child_cell_on_face(5,0)==4, ExcInternalError()); + + new_hexes[5]->set_neighbor (0, new_hexes[4]); + new_hexes[5]->set_neighbor (1, neighbor_cells[1][2]); + new_hexes[5]->set_neighbor (2, neighbor_cells[2][3]); + new_hexes[5]->set_neighbor (3, new_hexes[7]); + new_hexes[5]->set_neighbor (4, new_hexes[1]); + new_hexes[5]->set_neighbor (5, neighbor_cells[5][1]); + Assert(GeometryInfo::child_cell_on_face(1,2)==5, ExcInternalError()); + Assert(GeometryInfo::child_cell_on_face(2,3)==5, ExcInternalError()); + Assert(GeometryInfo::child_cell_on_face(5,1)==5, ExcInternalError()); + + new_hexes[6]->set_neighbor (0, neighbor_cells[0][3]); + new_hexes[6]->set_neighbor (1, new_hexes[7]); + new_hexes[6]->set_neighbor (2, new_hexes[4]); + new_hexes[6]->set_neighbor (3, neighbor_cells[3][1]); + new_hexes[6]->set_neighbor (4, new_hexes[2]); + new_hexes[6]->set_neighbor (5, neighbor_cells[5][2]); + Assert(GeometryInfo::child_cell_on_face(0,3)==6, ExcInternalError()); + Assert(GeometryInfo::child_cell_on_face(3,1)==6, ExcInternalError()); + Assert(GeometryInfo::child_cell_on_face(5,2)==6, ExcInternalError()); + + new_hexes[7]->set_neighbor (0, new_hexes[6]); new_hexes[7]->set_neighbor (1, neighbor_cells[1][3]); - new_hexes[7]->set_neighbor (2, new_hexes[4]); - new_hexes[7]->set_neighbor (3, new_hexes[6]); - new_hexes[7]->set_neighbor (4, neighbor_cells[4][3]); - new_hexes[7]->set_neighbor (5, neighbor_cells[5][2]); + new_hexes[7]->set_neighbor (2, new_hexes[5]); + new_hexes[7]->set_neighbor (3, neighbor_cells[3][3]); + new_hexes[7]->set_neighbor (4, new_hexes[3]); + new_hexes[7]->set_neighbor (5, neighbor_cells[5][3]); + Assert(GeometryInfo::child_cell_on_face(1,3)==7, ExcInternalError()); + Assert(GeometryInfo::child_cell_on_face(3,3)==7, ExcInternalError()); + Assert(GeometryInfo::child_cell_on_face(5,3)==7, ExcInternalError()); // now we need to set the @@ -7986,21 +8166,21 @@ void Triangulation<2>::delete_children (cell_iterator &cell) { // that they may appear at unwanted // places later on... cell->child(0)->line(1)->clear_user_pointer(); - cell->child(0)->line(2)->clear_user_pointer(); - cell->child(2)->line(0)->clear_user_pointer(); - cell->child(2)->line(3)->clear_user_pointer(); + cell->child(0)->line(3)->clear_user_pointer(); + cell->child(3)->line(0)->clear_user_pointer(); + cell->child(3)->line(2)->clear_user_pointer(); // same for user flags cell->child(0)->line(1)->clear_user_flag(); - cell->child(0)->line(2)->clear_user_flag(); - cell->child(2)->line(0)->clear_user_flag(); - cell->child(2)->line(3)->clear_user_flag(); + cell->child(0)->line(3)->clear_user_flag(); + cell->child(3)->line(0)->clear_user_flag(); + cell->child(3)->line(2)->clear_user_flag(); // delete the four interior lines cell->child(0)->line(1)->clear_used_flag(); - cell->child(0)->line(2)->clear_used_flag(); - cell->child(2)->line(0)->clear_used_flag(); - cell->child(2)->line(3)->clear_used_flag(); + cell->child(0)->line(3)->clear_used_flag(); + cell->child(3)->line(0)->clear_used_flag(); + cell->child(3)->line(2)->clear_used_flag(); // for the four faces: if the // neighbor does not itself need @@ -8114,7 +8294,7 @@ void Triangulation<3>::delete_children (cell_iterator &cell) { // is the vertex at the heart of // this cell, which is the sixth of // the first child - vertices_used[cell->child(0)->vertex_index(6)] = false; + vertices_used[cell->child(0)->vertex_index(7)] = false; /////////////////////////////////////// // delete interior quads and lines @@ -8123,30 +8303,29 @@ void Triangulation<3>::delete_children (cell_iterator &cell) { // line's and quad's indices const quad_iterator interior_quads[12] = { cell->child(0)->face(1), - cell->child(1)->face(1), cell->child(2)->face(1), - cell->child(3)->face(1), + cell->child(4)->face(1), + cell->child(6)->face(1), cell->child(0)->face(3), - cell->child(3)->face(3), + cell->child(1)->face(3), cell->child(4)->face(3), - cell->child(7)->face(3), + cell->child(5)->face(3), - cell->child(0)->face(4), - cell->child(1)->face(4), - cell->child(4)->face(4), - cell->child(5)->face(4) + cell->child(0)->face(5), + cell->child(1)->face(5), + cell->child(2)->face(5), + cell->child(3)->face(5) }; const line_iterator interior_lines[6] - = { cell->child(0)->line(10), - cell->child(4)->line(10), - - cell->child(0)->line(6), - cell->child(1)->line(6), + = { cell->child(0)->line(5), + cell->child(0)->line(7), + cell->child(0)->line(11), - cell->child(0)->line(5), - cell->child(3)->line(5) + cell->child(7)->line(0), + cell->child(7)->line(2), + cell->child(7)->line(8), }; // clear user pointers, to avoid that @@ -8183,13 +8362,13 @@ void Triangulation<3>::delete_children (cell_iterator &cell) { quad_iterator quad = cell->face(face); // delete middle vertex - vertices_used[quad->child(0)->vertex_index(2)] = false; + vertices_used[quad->child(0)->vertex_index(3)] = false; const line_iterator interior_lines[4] = { quad->child(0)->line(1), - quad->child(0)->line(2), - quad->child(2)->line(0), - quad->child(2)->line(3) + quad->child(0)->line(3), + quad->child(3)->line(0), + quad->child(3)->line(2) }; // delete interior lines diff --git a/deal.II/deal.II/source/grid/tria_accessor.cc b/deal.II/deal.II/source/grid/tria_accessor.cc index d5b62abf5a..598c4cfc94 100644 --- a/deal.II/deal.II/source/grid/tria_accessor.cc +++ b/deal.II/deal.II/source/grid/tria_accessor.cc @@ -267,8 +267,7 @@ int TriaObjectAccessor<2, dim>::vertex_index (const unsigned int corner) const { Assert (corner<4, ExcIndexRange(corner,0,4)); - static const int corner_convention[4] = { 0,0,1,1 }; - return line(corner)->vertex_index(corner_convention[corner]); + return line(corner%2)->vertex_index(corner/2); } @@ -437,8 +436,8 @@ bool TriaObjectAccessor<2, dim>::at_boundary () const template double TriaObjectAccessor<2, dim>::diameter () const { - return std::sqrt(std::max((vertex(2)-vertex(0)).square(), - (vertex(3)-vertex(1)).square())); + return std::sqrt(std::max((vertex(3)-vertex(0)).square(), + (vertex(2)-vertex(1)).square())); } @@ -476,9 +475,9 @@ Point<2> TriaObjectAccessor<2, 2>::barycenter () const x := array(0..3); y := array(0..3); tphi[0] := (1-xi)*(1-eta): - tphi[1] := xi*(1-eta): - tphi[2] := xi*eta: - tphi[3] := (1-xi)*eta: + tphi[1] := xi*(1-eta): + tphi[2] := (1-xi)*eta: + tphi[3] := xi*eta: x_real := sum(x[s]*tphi[s], s=0..3): y_real := sum(y[s]*tphi[s], s=0..3): detJ := diff(x_real,xi)*diff(y_real,eta) - diff(x_real,eta)*diff(y_real,xi): @@ -501,30 +500,22 @@ Point<2> TriaObjectAccessor<2, 2>::barycenter () const vertex(2)(1), vertex(3)(1) }; const double t1 = x[0]*x[1]; - const double t4 = x[2]*x[3]; - const double t7 = x[0]*x[3]; - const double t9 = x[1]*x[2]; - const double t11 = x[0]*x[0]; - const double t13 = x[1]*x[1]; - const double t15 = x[2]*x[2]; - const double t17 = x[3]*x[3]; - const double t25 = -t1*y[0]+t1*y[1]-t4*y[2]+t4*y[3]+t7*y[0]-t9*y[1] - +t11*y[1]-t13*y[0]+t15*y[3]-t17*y[2]-t11*y[3] - +t13*y[2]-t15*y[1]+t17*y[0]+t9*y[2]-t7*y[3]; - const double t35 = 1/(-x[0]*y[3]-x[1]*y[0]+x[2]*y[3]+x[3]*y[0] - +x[0]*y[1]-x[3]*y[2]+x[1]*y[2]-x[2]*y[1]); - const double t38 = y[3]*x[3]; - const double t40 = y[0]*x[0]; - const double t42 = y[1]*x[1]; - const double t44 = y[2]*x[2]; - const double t50 = y[0]*y[0]; - const double t52 = y[1]*y[1]; - const double t54 = y[2]*y[2]; - const double t56 = y[3]*y[3]; - const double t62 = -t38*y[2]+t40*y[1]-t42*y[0]+t44*y[3]-t40*y[3] - +t42*y[2]-t44*y[1]+t38*y[0]-x[1]*t50+x[0]*t52 - -x[3]*t54+x[2]*t56+x[3]*t50-x[2]*t52+x[1]*t54-x[0]*t56; - return Point<2> (t25*t35/3.0, t62*t35/3.0); + const double t3 = x[0]*x[0]; + const double t5 = x[1]*x[1]; + const double t9 = y[0]*x[0]; + const double t11 = y[1]*x[1]; + const double t14 = x[2]*x[2]; + const double t16 = x[3]*x[3]; + const double t20 = x[2]*x[3]; + const double t27 = t1*y[1]+t3*y[1]-t5*y[0]-t3*y[2]+t5*y[3]+t9*x[2]-t11*x[3]-t1*y[0]-t14*y[3]+t16*y[2]-t16*y[1]+t14*y[0]-t20*y[3]-x[0]*x[2]*y[2]+x[1]*x[3]*y[3]+t20*y[2]; + const double t37 = 1/(-x[1]*y[0]+x[1]*y[3]+y[0]*x[2]+x[0]*y[1]-x[0]*y[2]-y[1]*x[3]-x[2]*y[3]+x[3]*y[2]); + const double t39 = y[2]*y[2]; + const double t51 = y[0]*y[0]; + const double t53 = y[1]*y[1]; + const double t59 = y[3]*y[3]; + const double t63 = t39*x[3]+y[2]*y[0]*x[2]+y[3]*x[3]*y[2]-y[2]*x[2]*y[3]-y[3]*y[1]*x[3]-t9*y[2]+t11*y[3]+t51*x[2]-t53*x[3]-x[1]*t51+t9*y[1]-t11*y[0]+x[0]*t53-t59*x[2]+t59*x[1]-t39*x[0]; + + return Point<2> (t27*t37/3, t63*t37/3); } @@ -550,9 +541,9 @@ double TriaObjectAccessor<2, 2>::measure () const x := array(0..3); y := array(0..3); tphi[0] := (1-xi)*(1-eta): - tphi[1] := xi*(1-eta): - tphi[2] := xi*eta: - tphi[3] := (1-xi)*eta: + tphi[1] := xi*(1-eta): + tphi[2] := (1-xi)*eta: + tphi[3] := xi*eta: x_real := sum(x[s]*tphi[s], s=0..3): y_real := sum(y[s]*tphi[s], s=0..3): detJ := diff(x_real,xi)*diff(y_real,eta) - diff(x_real,eta)*diff(y_real,xi): @@ -572,8 +563,7 @@ double TriaObjectAccessor<2, 2>::measure () const vertex(2)(1), vertex(3)(1) }; - return (-x[0]*y[3]/2.0-x[1]*y[0]/2.0+x[2]*y[3]/2.0+x[3]*y[0]/2.0+ - x[0]*y[1]/2.0-x[3]*y[2]/2.0+x[1]*y[2]/2.0-x[2]*y[1]/2.0); + return (-x[1]*y[0]/2+x[1]*y[3]/2+y[0]*x[2]/2+x[0]*y[1]/2-x[0]*y[2]/2-y[1]*x[3]/2-x[2]*y[3]/2+x[3]*y[2]/2); } @@ -688,28 +678,28 @@ 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 front or the back + // 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 3 + // case, transpose vertices 1 and 2 // of this face static const unsigned int vertex_translation[4] = - { 0, 3, 2, 1 }; + { 0, 2, 1, 3 }; if (corner<4) { - if (face_orientation(0) == true) - return quad(0)->vertex_index(corner); + if (face_orientation(4) == true) + return quad(4)->vertex_index(corner); else - return quad(0)->vertex_index(vertex_translation[corner]); + return quad(4)->vertex_index(vertex_translation[corner]); } else { - if (face_orientation(1) == true) - return quad(1)->vertex_index(corner-4); + if (face_orientation(5) == true) + return quad(5)->vertex_index(corner-4); else - return quad(1)->vertex_index(vertex_translation[corner-4]); + return quad(5)->vertex_index(vertex_translation[corner-4]); } } @@ -877,10 +867,10 @@ bool TriaObjectAccessor<3, dim>::at_boundary () const template double TriaObjectAccessor<3, dim>::diameter () const { - return std::sqrt(std::max( std::max((vertex(6)-vertex(0)).square(), - (vertex(7)-vertex(1)).square()), - std::max((vertex(4)-vertex(2)).square(), - (vertex(5)-vertex(3)).square()) )); + return std::sqrt(std::max( std::max((vertex(7)-vertex(0)).square(), + (vertex(6)-vertex(1)).square()), + std::max((vertex(2)-vertex(5)).square(), + (vertex(3)-vertex(4)).square()) )); } @@ -946,32 +936,55 @@ Point<3> TriaObjectAccessor<3, 3>::barycenter () const readlib(C): C(array(1..3, [xs, ys, zs])); + + + This script takes more than several hours when using an old version + of maple on an old and slow computer. Therefore, when changing to + the new deal.II numbering scheme (lexicographic numbering) the code + lines below have not been reproduced with maple but only the + ordering of points in the definitions of x[], y[] and z[] have been + changed. + + For the case, someone is willing to rerun the maple script, he/she + should use following ordering of shape functions: + + tphi[0] := (1-xi)*(1-eta)*(1-zeta): + tphi[1] := xi*(1-eta)*(1-zeta): + tphi[2] := (1-xi)* eta*(1-zeta): + tphi[3] := xi* eta*(1-zeta): + tphi[4] := (1-xi)*(1-eta)*zeta: + tphi[5] := xi*(1-eta)*zeta: + tphi[6] := (1-xi)* eta*zeta: + tphi[7] := xi* eta*zeta: + + and change the ordering of points in the definitions of x[], y[] and + z[] back to the standard ordering. */ const double x[8] = { vertex(0)(0), vertex(1)(0), + vertex(5)(0), + vertex(4)(0), vertex(2)(0), vertex(3)(0), - vertex(4)(0), - vertex(5)(0), - vertex(6)(0), - vertex(7)(0) }; + vertex(7)(0), + vertex(6)(0) }; const double y[8] = { vertex(0)(1), vertex(1)(1), + vertex(5)(1), + vertex(4)(1), vertex(2)(1), vertex(3)(1), - vertex(4)(1), - vertex(5)(1), - vertex(6)(1), - vertex(7)(1) }; + vertex(7)(1), + vertex(6)(1) }; const double z[8] = { vertex(0)(2), vertex(1)(2), + vertex(5)(2), + vertex(4)(2), vertex(2)(2), vertex(3)(2), - vertex(4)(2), - vertex(5)(2), - vertex(6)(2), - vertex(7)(2) }; + vertex(7)(2), + vertex(6)(2) }; double s1, s2, s3, s4, s5, s6, s7, s8; @@ -1758,7 +1771,7 @@ bool CellAccessor<2>::point_inside (const Point<2> &p) const // we want the faces in counter // clockwise orientation - static const int direction[4]={1,1,-1,-1}; + static const int direction[4]={-1,1,1,-1}; for (unsigned int f=0; f<4; ++f) { // vector from the first vertex @@ -2145,7 +2158,7 @@ neighbor_child_on_subface (const unsigned int face, ExcMessage ("The neighbor must have children!")); static const unsigned int subface_translation[4] - = { 0, 3, 2, 1 }; + = { 0, 2, 1, 3 }; // see whether face and // the neighbor's // counterface share the diff --git a/deal.II/deal.II/source/grid/tria_boundary.cc b/deal.II/deal.II/source/grid/tria_boundary.cc index dd917dd1be..78e5b89276 100644 --- a/deal.II/deal.II/source/grid/tria_boundary.cc +++ b/deal.II/deal.II/source/grid/tria_boundary.cc @@ -237,8 +237,8 @@ get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &quad, for (unsigned int j=0; j::face_iterator &face, Boundary<3>::FaceVertexNormals &face_vertex_normals) const { const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face; - + + static const unsigned int neighboring_vertices[4][2]= + { {1,2},{3,0},{0,3},{2,1}}; for (unsigned int vertex=0; vertex::face_iterator &face, // radiating away from this // vertex const Tensor<1,3> tangents[2] - = { face->vertex((vertex+1) % vertices_per_face) + = { face->vertex(neighboring_vertices[vertex][0]) - face->vertex(vertex), - face->vertex((vertex+vertices_per_face-1) % vertices_per_face) + face->vertex(neighboring_vertices[vertex][1]) - face->vertex(vertex) }; // then compute the normal by diff --git a/deal.II/deal.II/source/grid/tria_boundary_lib.cc b/deal.II/deal.II/source/grid/tria_boundary_lib.cc index fe583730fd..f65247d4e7 100644 --- a/deal.II/deal.II/source/grid/tria_boundary_lib.cc +++ b/deal.II/deal.II/source/grid/tria_boundary_lib.cc @@ -152,16 +152,16 @@ CylinderBoundary<3>::get_intermediate_points_on_quad ( unsigned int m=static_cast (std::sqrt(static_cast(points.size()))); Assert(points.size()==m*m, ExcInternalError()); - std::vector > lp3(m); + std::vector > lp0(m); std::vector > lp1(m); - get_intermediate_points_on_line(quad->line(3), lp3); + get_intermediate_points_on_line(quad->line(0), lp0); get_intermediate_points_on_line(quad->line(1), lp1); std::vector > lps(m); for (unsigned int i=0; i::get_intermediate_points_on_quad ( unsigned int m=static_cast (std::sqrt(static_cast(points.size()))); Assert(points.size()==m*m, ExcInternalError()); - std::vector > lp3(m); + std::vector > lp0(m); std::vector > lp1(m); - get_intermediate_points_on_line(quad->line(3), lp3); + get_intermediate_points_on_line(quad->line(0), lp0); get_intermediate_points_on_line(quad->line(1), lp1); std::vector > lps(m); for (unsigned int i=0; i::build_some_patches (Data data) { Assert (patch != this->patches.end(), ExcInternalError()); + // use ucd_to_deal map as patch + // vertices are in the old, + // unnatural ordering for (unsigned int vertex=0; vertex::vertices_per_cell; ++vertex) - patch->vertices[vertex] = data.mapping->transform_unit_to_real_cell - (cell, GeometryInfo::unit_cell_vertex (vertex)); + patch->vertices[vertex] = data.mapping->transform_unit_to_real_cell + (cell, GeometryInfo::unit_cell_vertex ( + GeometryInfo::ucd_to_deal[vertex])); if (data.n_datasets > 0) { @@ -515,40 +519,8 @@ void DataOut::build_some_patches (Data data) // neighbor, so get its // patch number and set it // for the neighbor index - const unsigned int neighbor_patch_index - = this->patches[(*data.cell_to_patch_index_map) - [neighbor->level()][neighbor->index()]].patch_index; - - switch (dim) - { - case 1: - patch->neighbors[f] = neighbor_patch_index; - break; - - case 2: - switch (f) - { - case 0: patch->neighbors[2] = neighbor_patch_index; break; - case 1: patch->neighbors[1] = neighbor_patch_index; break; - case 2: patch->neighbors[3] = neighbor_patch_index; break; - case 3: patch->neighbors[0] = neighbor_patch_index; break; - } - break; - case 3: - switch (f) - { - case 0: patch->neighbors[2] = neighbor_patch_index; break; - case 1: patch->neighbors[3] = neighbor_patch_index; break; - case 2: patch->neighbors[4] = neighbor_patch_index; break; - case 3: patch->neighbors[1] = neighbor_patch_index; break; - case 4: patch->neighbors[5] = neighbor_patch_index; break; - case 5: patch->neighbors[0] = neighbor_patch_index; break; - } - break; - - default: - Assert(false, ExcNotImplemented()); - } + patch->neighbors[f] = this->patches[(*data.cell_to_patch_index_map) + [neighbor->level()][neighbor->index()]].patch_index; }; // next cell (patch) in this diff --git a/deal.II/deal.II/source/numerics/data_out_stack.cc b/deal.II/deal.II/source/numerics/data_out_stack.cc index 8cd39eff8d..aad6f5c855 100644 --- a/deal.II/deal.II/source/numerics/data_out_stack.cc +++ b/deal.II/deal.II/source/numerics/data_out_stack.cc @@ -345,7 +345,7 @@ void DataOutStack::build_patches (const unsigned int n_subdivisions) // cell to both the bottom // and top surfaces: // - // 3-------2 + // 2-------3 // / / // / / // / / @@ -354,7 +354,7 @@ void DataOutStack::build_patches (const unsigned int n_subdivisions) // i.e. the 3d vertex indices // 0, 1, ..., 7 are mapped to // the 2d indices 0, 1, 1, 0, - // 3, 2, 2, 3, but with + // 2, 3, 3, 2, but with // different z-coordinates patch->vertices[0] = Point(cell->vertex(0)(0), cell->vertex(0)(1), @@ -369,17 +369,17 @@ void DataOutStack::build_patches (const unsigned int n_subdivisions) cell->vertex(0)(1), parameter); - patch->vertices[4] = Point(cell->vertex(3)(0), - cell->vertex(3)(1), - parameter-parameter_step); - patch->vertices[5] = Point(cell->vertex(2)(0), + patch->vertices[4] = Point(cell->vertex(2)(0), cell->vertex(2)(1), parameter-parameter_step); - patch->vertices[6] = Point(cell->vertex(2)(0), - cell->vertex(2)(1), - parameter); - patch->vertices[7] = Point(cell->vertex(3)(0), + patch->vertices[5] = Point(cell->vertex(3)(0), cell->vertex(3)(1), + parameter-parameter_step); + patch->vertices[6] = Point(cell->vertex(3)(0), + cell->vertex(3)(1), + parameter); + patch->vertices[7] = Point(cell->vertex(2)(0), + cell->vertex(2)(1), parameter); break; diff --git a/tests/bits/fe_q_3d_01.cc b/tests/bits/fe_q_3d_01.cc index 024dcb7959..cf6610b62b 100644 --- a/tests/bits/fe_q_3d_01.cc +++ b/tests/bits/fe_q_3d_01.cc @@ -29,11 +29,11 @@ namespace FE_Q_3d { static const double constraint_q1[] = { - .25,.25,.25,.25, - .5,.5,0.,0., - 0.,.5,.5,0., - 0.,0.,.5,.5, - .5,0.,0.,.5 + 0.25,0.25,0.25,0.25, + 0.5,0,0.5,0, + 0,0.5,0,0.5, + 0.5,0.5,0,0, + 0,0,0.5,0.5 }; static const double constraint_q2[] = @@ -43,22 +43,22 @@ namespace FE_Q_3d 0,0,0,0,0,1,0,0,0, 0,0,0,0,0,0,1,0,0, 0,0,0,0,0,0,0,1,0, - 0,0,0,0,.375,0,-.125,0,.75, - 0,0,0,0,0,.375,0,-.125,.75, - 0,0,0,0,-.125,0,.375,0,.75, - 0,0,0,0,0,-.125,0,.375,.75, - .375,-.125,0,0,.75,0,0,0,0, - -.125,.375,0,0,.75,0,0,0,0, - 0,.375,-.125,0,0,.75,0,0,0, - 0,-.125,.375,0,0,.75,0,0,0, - 0,0,-.125,.375,0,0,.75,0,0, - 0,0,.375,-.125,0,0,.75,0,0, - .375,0,0,-.125,0,0,0,.75,0, - -.125,0,0,.375,0,0,0,.75,0, - .140625,-.046875,.015625,-.046875,.28125,-.09375,-.09375,.28125,.5625, - -.046875,.140625,-.046875,.015625,.28125,.28125,-.09375,-.09375,.5625, - .015625,-.046875,.140625,-.046875,-.09375,.28125,.28125,-.09375,.5625, - -.046875,.015625,-.046875,.140625,-.09375,-.09375,.28125,.28125,.5625 + 0,0,0,0,0,0,0.375,-0.125,0.75, + 0,0,0,0,0,0,-0.125,0.375,0.75, + 0,0,0,0,0.375,-0.125,0,0,0.75, + 0,0,0,0,-0.125,0.375,0,0,0.75, + 0.375,0,-0.125,0,0.75,0,0,0,0, + -0.125,0,0.375,0,0.75,0,0,0,0, + 0,0.375,0,-0.125,0,0.75,0,0,0, + 0,-0.125,0,0.375,0,0.75,0,0,0, + 0.375,-0.125,0,0,0,0,0.75,0,0, + -0.125,0.375,0,0,0,0,0.75,0,0, + 0,0,0.375,-0.125,0,0,0,0.75,0, + 0,0,-0.125,0.375,0,0,0,0.75,0, + 0.140625,-0.046875,-0.046875,0.015625,0.28125,-0.09375,0.28125,-0.09375,0.5625, + -0.046875,0.140625,0.015625,-0.046875,-0.09375,0.28125,0.28125,-0.09375,0.5625, + -0.046875,0.015625,0.140625,-0.046875,0.28125,-0.09375,-0.09375,0.28125,0.5625, + 0.015625,-0.046875,-0.046875,0.140625,-0.09375,0.28125,-0.09375,0.28125,0.5625 }; } diff --git a/tests/bits/mesh_3d_13.cc b/tests/bits/mesh_3d_13.cc index a08deabc3e..5202661742 100644 --- a/tests/bits/mesh_3d_13.cc +++ b/tests/bits/mesh_3d_13.cc @@ -58,7 +58,7 @@ void check_this (Triangulation<3> &tria) const bool orientation_flag = (neighbor->face_orientation(neighbor_neighbor) == true); static const unsigned int subface_translation[4] - = { 0, 3, 2, 1 }; + = { 0, 2, 1, 3 }; const unsigned int neighbor_child_index = (GeometryInfo<3>:: child_cell_on_face(neighbor_neighbor, diff --git a/tests/bits/mesh_3d_14.cc b/tests/bits/mesh_3d_14.cc index de5e7e12c4..2f1f8dd04b 100644 --- a/tests/bits/mesh_3d_14.cc +++ b/tests/bits/mesh_3d_14.cc @@ -80,7 +80,7 @@ void check_this (Triangulation<3> &tria) = (neighbor->face_orientation(neighbor_neighbor) == cell->face_orientation(face_no)); static const unsigned int subface_translation[4] - = { 0, 3, 2, 1 }; + = { 0, 2, 1, 3 }; const unsigned int neighbor_child_index = (GeometryInfo<3>:: child_cell_on_face(neighbor_neighbor, diff --git a/tests/bits/mesh_3d_15.cc b/tests/bits/mesh_3d_15.cc index 8b5a82d3de..69388e7d4e 100644 --- a/tests/bits/mesh_3d_15.cc +++ b/tests/bits/mesh_3d_15.cc @@ -73,7 +73,7 @@ void check_this (Triangulation<3> &tria) = (neighbor->face_orientation(neighbor_neighbor) == cell->face_orientation(face_no)); static const unsigned int subface_translation[4] - = { 0, 3, 2, 1 }; + = { 0, 2, 1, 3 }; const unsigned int neighbor_child_index = (GeometryInfo<3>:: child_cell_on_face(neighbor_neighbor, diff --git a/tests/bits/mesh_3d_2.cc b/tests/bits/mesh_3d_2.cc index 16ffc300ad..51b3f9ac1a 100644 --- a/tests/bits/mesh_3d_2.cc +++ b/tests/bits/mesh_3d_2.cc @@ -68,7 +68,7 @@ int main () // we know that from the second // cell, the common face must have // wrong orientation. check this - Assert ((++coarse_grid.begin_active())->face_orientation(4) + Assert ((++coarse_grid.begin_active())->face_orientation(5) == false, ExcInternalError()); } diff --git a/tests/bits/mesh_3d_3.cc b/tests/bits/mesh_3d_3.cc index 93d838ced4..8d6ea27cf2 100644 --- a/tests/bits/mesh_3d_3.cc +++ b/tests/bits/mesh_3d_3.cc @@ -55,13 +55,13 @@ int main () // make sure by hand that certain // vertices match up - Assert (cells[0]->vertex(1) == cells[1]->vertex(3), + Assert (cells[0]->vertex(1) == cells[1]->vertex(4), ExcInternalError()); - Assert (cells[0]->vertex(2) == cells[1]->vertex(2), + Assert (cells[0]->vertex(3) == cells[1]->vertex(6), ExcInternalError()); - Assert (cells[0]->vertex(5) == cells[1]->vertex(7), + Assert (cells[0]->vertex(5) == cells[1]->vertex(5), ExcInternalError()); - Assert (cells[0]->vertex(6) == cells[1]->vertex(6), + Assert (cells[0]->vertex(7) == cells[1]->vertex(7), ExcInternalError()); } diff --git a/tests/bits/q_points.cc b/tests/bits/q_points.cc index 1795458e44..096e61052b 100644 --- a/tests/bits/q_points.cc +++ b/tests/bits/q_points.cc @@ -66,7 +66,7 @@ void create_two_cubes (Triangulation<3> &coarse_grid) // finally generate a triangulation // out of this GridReordering<3>::reorder_cells (cells); - coarse_grid.create_triangulation (vertices, cells, SubCellData()); + coarse_grid.create_triangulation_compatibility (vertices, cells, SubCellData()); } diff --git a/tests/deal.II/constraints.cc b/tests/deal.II/constraints.cc index cbe4d7d864..a856d721a3 100644 --- a/tests/deal.II/constraints.cc +++ b/tests/deal.II/constraints.cc @@ -68,7 +68,8 @@ void make_tria (Triangulation<3> &tria, int step) cells[0].material_id = 0; cells[1].material_id = 0; - tria.create_triangulation (std::vector >(&vertices[0], &vertices[12]), + tria.create_triangulation_compatibility ( + std::vector >(&vertices[0], &vertices[12]), cells, SubCellData()); // no boundary information @@ -109,7 +110,8 @@ void make_tria (Triangulation<3> &tria, int step) cells[0].material_id = 0; cells[1].material_id = 0; - tria.create_triangulation (std::vector >(&vertices[0], &vertices[12]), + tria.create_triangulation_compatibility ( + std::vector >(&vertices[0], &vertices[12]), cells, SubCellData()); // no boundary information @@ -150,7 +152,8 @@ void make_tria (Triangulation<3> &tria, int step) cells[0].material_id = 0; cells[1].material_id = 0; - tria.create_triangulation (std::vector >(&vertices[0], &vertices[12]), + tria.create_triangulation_compatibility ( + std::vector >(&vertices[0], &vertices[12]), cells, SubCellData()); // no boundary information @@ -206,7 +209,8 @@ void make_tria (Triangulation<3> &tria, int step) cells[2].material_id = 0; cells[3].material_id = 0; - tria.create_triangulation (std::vector >(&vertices[0], &vertices[18]), + tria.create_triangulation_compatibility ( + std::vector >(&vertices[0], &vertices[18]), cells, SubCellData()); // no boundary information diff --git a/tests/deal.II/dof_test.cc b/tests/deal.II/dof_test.cc index a8a81c9fa6..46bd59cc17 100644 --- a/tests/deal.II/dof_test.cc +++ b/tests/deal.II/dof_test.cc @@ -264,8 +264,8 @@ void TestCases::run (const unsigned int test_case) { if (dim==3) { - tria->begin_active()->face(2)->set_boundary_indicator(1); tria->begin_active()->face(4)->set_boundary_indicator(1); + tria->begin_active()->face(5)->set_boundary_indicator(1); }; // set the boundary function diff --git a/tests/deal.II/gradients.cc b/tests/deal.II/gradients.cc index 2ad6af4234..27455c75c4 100644 --- a/tests/deal.II/gradients.cc +++ b/tests/deal.II/gradients.cc @@ -40,7 +40,7 @@ int main () Triangulation<2> tria; GridGenerator::hyper_cube (tria,0,1); - tria.begin_active()->vertex(2)(0) = 2; + tria.begin_active()->vertex(3)(0) = 2; FE_Q<2> fe(1); DoFHandler<2> dof(tria); @@ -84,17 +84,17 @@ int main () (grads[3] == Point<2>(0,-1))); break; case 2: - ok = ((grads[0] == Point<2>(0,0)) && - (grads[1] == Point<2>(0.5,0)) && - (grads[2] == Point<2>(0,1)) && - (grads[3] == Point<2>(0.5,0.5))); - break; - case 3: ok = ((grads[0] == Point<2>(0,1)) && (grads[1] == Point<2>(-0.5,1)) && (grads[2] == Point<2>(0,0)) && (grads[3] == Point<2>(-0.5,0.5))); break; + case 3: + ok = ((grads[0] == Point<2>(0,0)) && + (grads[1] == Point<2>(0.5,0)) && + (grads[2] == Point<2>(0,1)) && + (grads[3] == Point<2>(0.5,0.5))); + break; }; deallog << " Shape function " << vertex diff --git a/tests/deal.II/grid_test.cc b/tests/deal.II/grid_test.cc index e0ba4a3959..4c11c9122f 100644 --- a/tests/deal.II/grid_test.cc +++ b/tests/deal.II/grid_test.cc @@ -202,10 +202,9 @@ void test (const int test_case) for (int i=0; i<(dim==2 ? 3 : 2); ++i) { // refine the presently - // second last cell 17 + // last cell several // times cell = tria.last_active(tria.n_levels()-1); - --cell; cell->set_refine_flag (); tria.execute_coarsening_and_refinement (); }; @@ -218,8 +217,8 @@ void test (const int test_case) { if (dim==3) { - tria.begin_active()->face(2)->set_boundary_indicator(1); tria.begin_active()->face(4)->set_boundary_indicator(1); + tria.begin_active()->face(5)->set_boundary_indicator(1); }; @@ -254,11 +253,10 @@ void test (const int test_case) tria.set_boundary (1); break; } - }; - + } GridOut().write_ucd (tria, logfile); - + deallog << " Total number of cells = " << tria.n_cells() << std::endl << " Total number of active cells = " << tria.n_active_cells() << std::endl; diff --git a/tests/deal.II/measure_et_al.cc b/tests/deal.II/measure_et_al.cc index 6842dffa5c..362210c57d 100644 --- a/tests/deal.II/measure_et_al.cc +++ b/tests/deal.II/measure_et_al.cc @@ -48,7 +48,7 @@ void create_triangulation(const unsigned int case_no, GridGenerator::hyper_cube(tria, 1., 3.); Point<2> &v0=tria.begin_active()->vertex(0); v0(0) = 0.; - Point<2> &v2=tria.begin_active()->vertex(2); + Point<2> &v2=tria.begin_active()->vertex(3); v2(0) = 5.; v2(1) = 4.; // exact_areas.push_back(7.); diff --git a/tests/deal.II/second_derivatives.cc b/tests/deal.II/second_derivatives.cc index ba1ff4f61f..4735a766eb 100644 --- a/tests/deal.II/second_derivatives.cc +++ b/tests/deal.II/second_derivatives.cc @@ -64,12 +64,12 @@ int main () if (loop==1) { tria.begin_active()->vertex(1)(0) = 2; - tria.begin_active()->vertex(2)(0) = 2; + tria.begin_active()->vertex(3)(0) = 2; } if (loop==2) { - tria.begin_active()->vertex(2)(0) = 3; - tria.begin_active()->vertex(2)(1) = 3; + tria.begin_active()->vertex(3)(0) = 3; + tria.begin_active()->vertex(3)(1) = 3; } FEValues<2> fevalues(fe,q,update_second_derivatives); diff --git a/tests/deal.II/subcelldata.cc b/tests/deal.II/subcelldata.cc index e0c5e5def1..f34557f622 100644 --- a/tests/deal.II/subcelldata.cc +++ b/tests/deal.II/subcelldata.cc @@ -21,6 +21,15 @@ #include +static const unsigned subcells[6][4] = {{0, 1, 2, 3}, + {4, 5, 6, 7}, + {0, 1, 5, 4}, + {1, 5, 6, 2}, + {3, 2, 6, 7}, + {0, 4, 7, 3}}; + + + template void test() { @@ -75,14 +84,13 @@ void test() for(unsigned int f=0; f::faces_per_cell; ++f) { for(unsigned int i=0; i::vertices_per_face; ++i) - subcelldata.boundary_quads[f].vertices[i]= - GeometryInfo::face_to_cell_vertices(f,i); + subcelldata.boundary_quads[f].vertices[i]=subcells[f][i]; subcelldata.boundary_quads[f].material_id=10*f+1; } } Triangulation tria; - tria.create_triangulation (vertices, cells, subcelldata); + tria.create_triangulation_compatibility (vertices, cells, subcelldata); GridOutFlags::Ucd ucd_flags(true,true); GridOut grid_out; diff --git a/tests/fe/mapping.cc b/tests/fe/mapping.cc index 7346a60d95..043b38bc33 100644 --- a/tests/fe/mapping.cc +++ b/tests/fe/mapping.cc @@ -257,15 +257,15 @@ void create_triangulations(std::vector *> &tria_ptr, Point<2>(right,left), Point<2>(right,right), Point<2>(left,right) }; - const int cell_vertices[1][4] = { { 1,2,3,0 } }; + const int cell_vertices[1][4] = { { 1,2,0,3 } }; std::vector > cells (1, CellData<2>()); for (unsigned int j=0; j<4; ++j) cells[0].vertices[j] = cell_vertices[0][j]; cells[0].material_id = 0; tria->create_triangulation (std::vector >(&vertices[0], &vertices[4]), - cells, - SubCellData()); + cells, + SubCellData()); exact_areas.push_back(9.); } @@ -277,9 +277,9 @@ void create_triangulations(std::vector *> &tria_ptr, GridGenerator::hyper_cube(*tria, 1., 3.); Point<2> &v0=tria->begin_quad()->vertex(0); v0(0) = 0.; - Point<2> &v2=tria->begin_quad()->vertex(2); - v2(0) = 5.; - v2(1) = 4.; + Point<2> &v3=tria->begin_quad()->vertex(3); + v3(0) = 5.; + v3(1) = 4.; exact_areas.push_back(7.); show[1][0] = 1; } @@ -296,14 +296,14 @@ void create_triangulations(std::vector *> &tria_ptr, GridGenerator::hyper_cube(*tria, 1., 5.); Point<2> &v2 = tria->begin_active()->vertex(2); Point<2> &v3 = tria->begin_active()->vertex(3); - v2(0) = 3.; + v2(0) = 1.; v2(1) = 3.; - v3(0) = 1.; + v3(0) = 3.; v3(1) = 3.; tria->set_boundary(1,*boundary1); tria->set_boundary(2,*boundary2); tria->begin_active()->face(1)->set_boundary_indicator(1); - tria->begin_active()->face(2)->set_boundary_indicator(2); + tria->begin_active()->face(3)->set_boundary_indicator(2); double pi=std::acos(-1.); double alpha=2*std::atan(0.5); exact_areas.push_back(4+pi-2.5*(alpha-std::sin(alpha))); @@ -374,7 +374,7 @@ void create_triangulations(std::vector *> &tria_ptr, tria=new Triangulation<3>(); tria_ptr.push_back(tria); GridGenerator::hyper_cube(*tria, 1., 3.); - Point<3> &v=tria->begin()->vertex(6); + Point<3> &v=tria->begin()->vertex(7); v(0) = 5.; v(1) = 4.; v(2) = 4.5; @@ -396,7 +396,7 @@ void create_triangulations(std::vector *> &tria_ptr, tria_ptr.push_back(tria); GridGenerator::hyper_cube(*tria, 1., 3.); tria->set_boundary(1,*boundary1); - tria->begin_active()->face(3)->set_boundary_indicator(1); + tria->begin_active()->face(1)->set_boundary_indicator(1); exact_areas.push_back(8.+pi/3*h*h*(3*r-h)); } diff --git a/tests/fe/mapping_q1_eulerian.cc b/tests/fe/mapping_q1_eulerian.cc index 29d36f59f2..4c49dd1dc5 100644 --- a/tests/fe/mapping_q1_eulerian.cc +++ b/tests/fe/mapping_q1_eulerian.cc @@ -32,7 +32,7 @@ show_values(FiniteElement& fe, // shift one point of the cell // somehow if (dim > 1) - tr.begin_active()->vertex(2)(dim-1) += 1./std::sqrt(2.); + tr.begin_active()->vertex(dim==2 ? 3 : 5)(dim-1) += 1./std::sqrt(2.); DoFHandler dof(tr); dof.distribute_dofs(fe); diff --git a/tests/fe/numbering.cc b/tests/fe/numbering.cc index 179b4b353b..69023fa3d3 100644 --- a/tests/fe/numbering.cc +++ b/tests/fe/numbering.cc @@ -26,9 +26,22 @@ void check (const FE_Q &fe) const unsigned int dofs_per_cell = fe.dofs_per_cell; // polynomial degree const unsigned int degree = fe.dofs_per_line+1; - // number of grid points in each direction + // number of grid points in each + // direction const unsigned int n = degree+1; + // the following lines of code are + // somewhat odd, due to the way the + // hierarchic numbering is + // organized. if someone would + // really want to understand these + // lines, you better draw some + // pictures where you indicate the + // indices and orders of vertices, + // lines, etc, along with the + // numbers of the degrees of + // freedom in hierarchical and + // lexicographical order switch (dim) { case 1: @@ -39,7 +52,7 @@ void check (const FE_Q &fe) hierarchic_to_lexicographic_numbering[i] = i-1; break; - }; + } case 2: { @@ -47,20 +60,25 @@ void check (const FE_Q &fe) // first the four vertices hierarchic_to_lexicographic_numbering[next_index++] = 0; hierarchic_to_lexicographic_numbering[next_index++] = n-1; - hierarchic_to_lexicographic_numbering[next_index++] = n*n-1; hierarchic_to_lexicographic_numbering[next_index++] = n*(n-1); - // first line + hierarchic_to_lexicographic_numbering[next_index++] = n*n-1; + + // left line for (unsigned int i=0; i &fe) for (unsigned int j=0; j &fe) for (unsigned int i=0; i