+//! Topological description of cells.
+/**
+ * Dimension independent base class for the <tt>GeometryInfo<dim></tt>
+ * classes with <tt>dim=1,2,3,4</tt>. 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
+ * <tt>dim=max_dim</tt>.
+ */
+ 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
+ * <tt>dim</tt> dimension these
+ * are the <tt>2*dim</tt> first
+ * entries of
+ * <tt>{0,0,1,1,2,2,3,3}</tt>.
+ *
+ * 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
+ * <tt>dim</tt> dimension these
+ * are the <tt>2*dim<tt> first
+ * entries of
+ * <tt>{-1,1,-1,1,-1,1,-1,1}</tt>.
+ *
+ * Each value is either
+ * <tt>1</tt> or <tt>-1</tt>,
+ * 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
+ * <tt>{ 1, 0, 3, 2, 5, 4, 7, 6}</tt>.
+ */
+ 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 <int dim>
-struct GeometryInfo
+struct GeometryInfo: public GeometryInfoBase
{
/**
GeometryInfo<dim-1>::quads_per_cell);
/**
- * List of numbers which denotes which
- * face is opposite to a given face. In
- * 1d, this list is <tt>{1,0}</tt>, in 2d
- * <tt>{2, 3, 0, 1}</tt>, in 3d <tt>{1,
- * 0, 4, 5, 2, 3}</tt>.
+ * 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
+ * <tt>deal.II</tt> 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
*/
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
- * <tt>1</tt> or <tt>-1</tt>,
- * 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
* 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
* 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
* cell vertex number of the
* <tt>vertex</tt>th vertex of
* line <tt>line</tt>, e.g.
- * <tt>GeometryInfo<2>::line_to_cell_vertices(2,0)=3</tt>.
+ * <tt>GeometryInfo<2>::line_to_cell_vertices(3,0)=2</tt>.
*
* The order of the lines, as
* well as their direction (which
* cell vertex number of the
* <tt>vertex</tt>th vertex of
* face <tt>face</tt>, e.g.
- * <tt>GeometryInfo<2>::face_to_cell_vertices(2,0)=3</tt>.
+ * <tt>GeometryInfo<2>::face_to_cell_vertices(3,0)=2</tt>.
*
* Through the
* <tt>face_orientation</tt>
* cell line number of the
* <tt>line</tt>th line of face
* <tt>face</tt>, e.g.
- * <tt>GeometryInfo<3>::face_to_cell_lines(3,1)=5</tt>.
+ * <tt>GeometryInfo<3>::face_to_cell_lines(5,0)=4</tt>.
*
* Through the
* <tt>face_orientation</tt>
* orientation.
* <tt>face_orientation</tt>
* defaults to <tt>true</tt>
- * (standard orientation).
- *
- * This function is useful and
- * implemented for
- * <tt>dim=3</tt>, only.
+ * (standard orientation) and has
+ * no effect in 2d.
*/
static unsigned int face_to_cell_lines (const unsigned int face,
const unsigned int line,
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<double>(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);
}
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 <int dim>
inline
Point<dim>
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));
}
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)));
}
bool
DataOutBase::Patch<dim,spacedim>::operator == (const Patch &patch) const
{
+ const double eps=3e-16;
for (unsigned int i=0; i<GeometryInfo<dim>::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<GeometryInfo<dim>::faces_per_cell; ++i)
+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 <>
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 <>
Assert (subface<subfaces_per_face, ExcIndexRange(subface, 0, subfaces_per_face));
static const unsigned
- subcells[faces_per_cell][subfaces_per_face] = {{0,1},
- {1,2},
- {3,2},
- {0,3}};
+ subcells[faces_per_cell][subfaces_per_face] = {{0,2},
+ {1,3},
+ {0,1},
+ {2,3}};
+
return subcells[face][subface];
}
Assert (face<faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
Assert (subface<subfaces_per_face, ExcIndexRange(subface, 0, subfaces_per_face));
- static const unsigned int flip[subfaces_per_face] = { 0, 3, 2, 1 };
+ static const unsigned int flip[subfaces_per_face] = { 0, 2, 1, 3 };
static const unsigned
- subcells[faces_per_cell][subfaces_per_face] = {{0, 1, 2, 3},
- {4, 5, 6, 7},
- {0, 1, 5, 4},
- {1, 5, 6, 2},
- {3, 2, 6, 7},
- {0, 4, 7, 3}};
+ subcells[faces_per_cell][subfaces_per_face] = {{0, 2, 4, 6},
+ {1, 3, 5, 7},
+ {0, 4, 1, 5},
+ {2, 6, 3, 7},
+ {0, 1, 2, 3},
+ {4, 5, 6, 7}};
return face_orientation ? subcells[face][subface] : subcells[face][flip[subface]];
}
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;
const unsigned int vertex)
{
Assert (line<lines_per_cell, ExcIndexRange(line, 0, lines_per_cell));
- Assert (vertex<2, ExcIndexRange(vertex, 0, 2));
+ Assert (vertex<2, ExcIndexRange(vertex, 0, 2));
static const unsigned
- vertices[lines_per_cell][2] = {{0, 1}, // front face
- {1, 2},
- {3, 2},
- {0, 3},
- {4, 5}, // back face
- {5, 6},
- {7, 6},
- {4, 7},
- {0, 4}, // connects of front and back face
- {1, 5},
+ vertices[lines_per_cell][2] = {{0, 2}, // bottom face
+ {1, 3},
+ {0, 1},
+ {2, 3},
+ {4, 6}, // top face
+ {5, 7},
+ {4, 5},
+ {6, 7},
+ {0, 4}, // connects of bottom
+ {1, 5}, // top face
{2, 6},
{3, 7}};
return vertices[line][vertex];
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;
Assert (face<faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
Assert (line<lines_per_face, ExcIndexRange(line, 0, lines_per_face));
+ // also represented by (line+2)%4
+ static const unsigned int flip[lines_per_face] = { 2, 3, 0, 1 };
+
static const unsigned
- lines[faces_per_cell][lines_per_face] = {{0, 1, 2, 3},
- {4, 5, 6, 7},
- {0, 9, 4, 8},
- {9, 5,10, 1},
- {2,10, 6,11},
- {8, 7,11, 3}};
-
- return face_orientation ? lines[face][line] : lines[face][3-line];
+ lines[faces_per_cell][lines_per_face] = {{8,10, 0, 4}, // left face
+ {9,11, 1, 5}, // right face
+ {2, 6, 8, 9}, // front face
+ {3, 7,10,11}, // back face
+ {0, 1, 2, 3}, // bottom face
+ {4, 5, 6, 7}};// top face
+ return face_orientation ? lines[face][line] : lines[face][flip[line]];
}
-template <>
+template<int dim>
unsigned int
-GeometryInfo<4>::face_to_cell_lines (const unsigned int,
- const unsigned int,
- const bool)
+GeometryInfo<dim>::face_to_cell_lines (const unsigned int,
+ const unsigned int,
+ const bool)
{
Assert(false, ExcNotImplemented());
return invalid_unsigned_int;
switch (face_no)
{
case 0:
- q_points[p] = Point<dim>(quadrature.point(p)(0),0);
+ q_points[p] = Point<dim>(0,quadrature.point(p)(0));
break;
case 1:
q_points[p] = Point<dim>(1,quadrature.point(p)(0));
break;
case 2:
- q_points[p] = Point<dim>(quadrature.point(p)(0),1);
+ q_points[p] = Point<dim>(quadrature.point(p)(0),0);
break;
case 3:
- q_points[p] = Point<dim>(0,quadrature.point(p)(0));
+ q_points[p] = Point<dim>(quadrature.point(p)(0),1);
break;
default:
Assert (false, ExcInternalError());
switch (face_no)
{
case 0:
- q_points[p] = Point<dim>(quadrature.point(p)(0),
- 0,
+ q_points[p] = Point<dim>(0,
+ quadrature.point(p)(0),
quadrature.point(p)(1));
break;
case 1:
- q_points[p] = Point<dim>(quadrature.point(p)(0),
- 1,
+ q_points[p] = Point<dim>(1,
+ quadrature.point(p)(0),
quadrature.point(p)(1));
break;
case 2:
- q_points[p] = Point<dim>(quadrature.point(p)(0),
- quadrature.point(p)(1),
- 0);
+ q_points[p] = Point<dim>(quadrature.point(p)(1),
+ 0,
+ quadrature.point(p)(0));
break;
case 3:
- q_points[p] = Point<dim>(1,
- quadrature.point(p)(0),
- quadrature.point(p)(1));
+ q_points[p] = Point<dim>(quadrature.point(p)(1),
+ 1,
+ quadrature.point(p)(0));
break;
case 4:
q_points[p] = Point<dim>(quadrature.point(p)(0),
quadrature.point(p)(1),
- 1);
+ 0);
break;
case 5:
- q_points[p] = Point<dim>(0,
- quadrature.point(p)(0),
- quadrature.point(p)(1));
- break;
+ q_points[p] = Point<dim>(quadrature.point(p)(0),
+ quadrature.point(p)(1),
+ 1);
+ break;
default:
Assert (false, ExcInternalError());
switch (subface_no)
{
case 0:
- q_points[p]
- = Point<dim>(quadrature.point(p)(0)/2,0);
+ q_points[p] = Point<dim>(0,quadrature.point(p)(0)/2);
break;
case 1:
- q_points[p]
- = Point<dim>(quadrature.point(p)(0)/2+0.5,0);
+ q_points[p] = Point<dim>(0,quadrature.point(p)(0)/2+0.5);
break;
default:
Assert (false, ExcInternalError());
switch (subface_no)
{
case 0:
- q_points[p] = Point<dim>(quadrature.point(p)(0)/2,1);
+ q_points[p]
+ = Point<dim>(quadrature.point(p)(0)/2,0);
break;
case 1:
- q_points[p] = Point<dim>(quadrature.point(p)(0)/2+0.5,1);
+ q_points[p]
+ = Point<dim>(quadrature.point(p)(0)/2+0.5,0);
break;
default:
Assert (false, ExcInternalError());
switch (subface_no)
{
case 0:
- q_points[p] = Point<dim>(0,quadrature.point(p)(0)/2);
+ q_points[p] = Point<dim>(quadrature.point(p)(0)/2,1);
break;
case 1:
- q_points[p] = Point<dim>(0,quadrature.point(p)(0)/2+0.5);
+ q_points[p] = Point<dim>(quadrature.point(p)(0)/2+0.5,1);
break;
default:
Assert (false, ExcInternalError());
};
break;
+
default:
Assert (false, ExcInternalError());
};
switch (face_no)
{
case 0:
- q_points[p] = Point<dim>(quadrature.point(p)(0)/2,
- 0,
+ q_points[p] = Point<dim>(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<dim>(quadrature.point(p)(0)/2,
- 1,
+ q_points[p] = Point<dim>(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:
};
break;
case 2:
- q_points[p] = Point<dim>(quadrature.point(p)(0)/2,
- quadrature.point(p)(1)/2,
- 0);
+ q_points[p] = Point<dim>(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<dim>(1,
- quadrature.point(p)(0)/2,
- quadrature.point(p)(1)/2);
+ q_points[p] = Point<dim>(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:
case 4:
q_points[p] = Point<dim>(quadrature.point(p)(0)/2,
quadrature.point(p)(1)/2,
- 1);
+ 0);
switch (subface_no)
{
case 0:
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:
};
break;
case 5:
- q_points[p] = Point<dim>(0,
- quadrature.point(p)(0)/2,
- quadrature.point(p)(1)/2);
+ q_points[p] = Point<dim>(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());
};
* 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
* Furthermore, this class assumes that shape functions of the
* FiniteElement under consideration do <em>not</em> 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 <tt>update_values</tt>. For FiniteElements whose shape
+ * functions depend on the cells in real space, the update_once() and
* update_each() functions must be overloaded.
*
* Todos:
*
* <li> 2D case:
* @verbatim
- * 3-------2
+ * 2-------3
* | |
* | |
* | |
*
* <li> 3D case:
* @verbatim
- * 7-------6 7-------6
+ * 6-------7 6-------7
* /| | / /|
* / | | / / |
* / | | / / |
- * 3 | | 3-------2 |
- * | 4-------5 | | 5
+ * 4 | | 4-------5 |
+ * | 2-------3 | | 3
* | / / | | /
* | / / | | /
* |/ / | |/
* <ul>
* <li> Index 0: <tt>[0, 0, 0]</tt>;
* <li> Index 1: <tt>[1, 0, 0]</tt>;
- * <li> Index 2: <tt>[1, 0, 1]</tt>;
- * <li> Index 3: <tt>[0, 0, 1]</tt>;
- * <li> Index 4: <tt>[0, 1, 0]</tt>;
- * <li> Index 5: <tt>[1, 1, 0]</tt>;
- * <li> Index 6: <tt>[1, 1, 1]</tt>;
- * <li> Index 7: <tt>[0, 1, 1]</tt>;
+ * <li> Index 2: <tt>[0, 1, 0]</tt>;
+ * <li> Index 3: <tt>[1, 1, 0]</tt>;
+ * <li> Index 4: <tt>[0, 0, 1]</tt>;
+ * <li> Index 5: <tt>[1, 0, 1]</tt>;
+ * <li> Index 6: <tt>[0, 1, 1]</tt>;
+ * <li> Index 7: <tt>[1, 1, 1]</tt>;
* </ul>
* </ul>
* <h4>Q2 elements</h4>
*
* <li> 2D case:
* @verbatim
- * 3---6---2
+ * 2---7---3
* | |
- * 7 8 5
+ * 4 8 5
* | |
- * 0---4---1
+ * 0---6---1
* @endverbatim
*
* <li> 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
* <ul>
* <li> Index 0: <tt>[0, 0, 0]</tt>;
* <li> Index 1: <tt>[1, 0, 0]</tt>;
- * <li> Index 2: <tt>[1, 0, 1]</tt>;
- * <li> Index 3: <tt>[0, 0, 1]</tt>;
- * <li> Index 4: <tt>[0, 1, 0]</tt>;
- * <li> Index 5: <tt>[1, 1, 0]</tt>;
- * <li> Index 6: <tt>[1, 1, 1]</tt>;
- * <li> Index 7: <tt>[0, 1, 1]</tt>;
- * <li> Index 8: <tt>[1/2, 0, 0]</tt>;
- * <li> Index 9: <tt>[1, 0, 1/2]</tt>;
- * <li> Index 10: <tt>[1/2, 0, 1]</tt>;
- * <li> Index 11: <tt>[0, 0, 1/2]</tt>;
- * <li> Index 12: <tt>[1/2, 1, 0]</tt>;
- * <li> Index 13: <tt>[1, 1, 1/2]</tt>;
- * <li> Index 14: <tt>[1/2, 1, 1]</tt>;
- * <li> Index 15: <tt>[0, 1, 1/2]</tt>;
- * <li> Index 16: <tt>[0, 1/2, 0]</tt>;
- * <li> Index 17: <tt>[1, 1/2, 0]</tt>;
- * <li> Index 18: <tt>[1, 1/2, 1]</tt>;
- * <li> Index 19: <tt>[0, 1/2, 1]</tt>;
- * <li> Index 20: <tt>[1/2, 0, 1/2]</tt>;
- * <li> Index 21: <tt>[1/2, 1, 1/2]</tt>;
- * <li> Index 22: <tt>[1/2, 1/2, 0]</tt>;
- * <li> Index 23: <tt>[1, 1/2, 1/2]</tt>;
- * <li> Index 24: <tt>[1/2, 1/2, 1]</tt>;
- * <li> Index 25: <tt>[0, 1/2, 1/2]</tt>;
+ * <li> Index 2: <tt>[0, 1, 0]</tt>;
+ * <li> Index 3: <tt>[1, 1, 0]</tt>;
+ * <li> Index 4: <tt>[0, 0, 1]</tt>;
+ * <li> Index 5: <tt>[1, 0, 1]</tt>;
+ * <li> Index 6: <tt>[0, 1, 1]</tt>;
+ * <li> Index 7: <tt>[1, 1, 1]</tt>;
+ * <li> Index 8: <tt>[0, 1/2, 0]</tt>;
+ * <li> Index 9: <tt>[1, 1/2, 0]</tt>;
+ * <li> Index 10: <tt>[1/2, 0, 0]</tt>;
+ * <li> Index 11: <tt>[1/2, 1, 0]</tt>;
+ * <li> Index 12: <tt>[0, 1/2, 1]</tt>;
+ * <li> Index 13: <tt>[1, 1/2, 1]</tt>;
+ * <li> Index 14: <tt>[1/2, 0, 1]</tt>;
+ * <li> Index 15: <tt>[1/2, 1, 1]</tt>;
+ * <li> Index 16: <tt>[0, 0, 1/2]</tt>;
+ * <li> Index 17: <tt>[1, 0, 1/2]</tt>;
+ * <li> Index 18: <tt>[0, 1, 1/2]</tt>;
+ * <li> Index 19: <tt>[1, 1, 1/2]</tt>;
+ * <li> Index 20: <tt>[0, 1/2, 1/2]</tt>;
+ * <li> Index 21: <tt>[1, 1/2, 1/2]</tt>;
+ * <li> Index 22: <tt>[1/2, 0, 1/2]</tt>;
+ * <li> Index 23: <tt>[1/2, 1, 1/2]</tt>;
+ * <li> Index 24: <tt>[1/2, 1/2, 0]</tt>;
+ * <li> Index 25: <tt>[1/2, 1/2, 1]</tt>;
* <li> Index 26: <tt>[1/2, 1/2, 1/2]</tt>;
* </ul>
* </ul>
*
* <li> 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.
* </ul>
* <h4>Q4 elements</h4>
* <ul>
*
* <li> 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
* </ul>
- * 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 <int dim>
class FE_Q : public FE_Poly<TensorProductPolynomials<dim>,dim>
*
* <li> 2D case:
* @verbatim
- * 3-------2
+ * 2-------3
* | |
* | |
* | |
*
* <li> 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:
* <ul>
* <li> Index 0: <tt>[0, 0, 0]</tt>;
* <li> Index 1: <tt>[1, 0, 0]</tt>;
- * <li> Index 2: <tt>[1, 0, 1]</tt>;
- * <li> Index 3: <tt>[0, 0, 1]</tt>;
- * <li> Index 4: <tt>[0, 1, 0]</tt>;
- * <li> Index 5: <tt>[1, 1, 0]</tt>;
- * <li> Index 6: <tt>[1, 1, 1]</tt>;
- * <li> Index 7: <tt>[0, 1, 1]</tt>;
+ * <li> Index 2: <tt>[0, 1, 0]</tt>;
+ * <li> Index 3: <tt>[1, 1, 0]</tt>;
+ * <li> Index 4: <tt>[0, 0, 1]</tt>;
+ * <li> Index 5: <tt>[1, 0, 1]</tt>;
+ * <li> Index 6: <tt>[0, 1, 1]</tt>;
+ * <li> Index 7: <tt>[1, 1, 1]</tt>;
* </ul>
* </ul>
* <h4>Q2 elements</h4>
*
* <li> 2D case:
* @verbatim
- * 3---6---2
+ * 2---7---3
* | |
- * 7 8 5
+ * 4 8 5
* | |
- * 0---4---1
+ * 0---6---1
* @endverbatim
*
* <li> 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
* <ul>
* <li> Index 0: <tt>[0, 0, 0]</tt>;
* <li> Index 1: <tt>[1, 0, 0]</tt>;
- * <li> Index 2: <tt>[1, 0, 1]</tt>;
- * <li> Index 3: <tt>[0, 0, 1]</tt>;
- * <li> Index 4: <tt>[0, 1, 0]</tt>;
- * <li> Index 5: <tt>[1, 1, 0]</tt>;
- * <li> Index 6: <tt>[1, 1, 1]</tt>;
- * <li> Index 7: <tt>[0, 1, 1]</tt>;
- * <li> Index 8: <tt>[1/2, 0, 0]</tt>;
- * <li> Index 9: <tt>[1, 0, 1/2]</tt>;
- * <li> Index 10: <tt>[1/2, 0, 1]</tt>;
- * <li> Index 11: <tt>[0, 0, 1/2]</tt>;
- * <li> Index 12: <tt>[1/2, 1, 0]</tt>;
- * <li> Index 13: <tt>[1, 1, 1/2]</tt>;
- * <li> Index 14: <tt>[1/2, 1, 1]</tt>;
- * <li> Index 15: <tt>[0, 1, 1/2]</tt>;
- * <li> Index 16: <tt>[0, 1/2, 0]</tt>;
- * <li> Index 17: <tt>[1, 1/2, 0]</tt>;
- * <li> Index 18: <tt>[1, 1/2, 1]</tt>;
- * <li> Index 19: <tt>[0, 1/2, 1]</tt>;
- * <li> Index 20: <tt>[1/2, 0, 1/2]</tt>;
- * <li> Index 21: <tt>[1/2, 1, 1/2]</tt>;
- * <li> Index 22: <tt>[1/2, 1/2, 0]</tt>;
- * <li> Index 23: <tt>[1, 1/2, 1/2]</tt>;
- * <li> Index 24: <tt>[1/2, 1/2, 1]</tt>;
- * <li> Index 25: <tt>[0, 1/2, 1/2]</tt>;
+ * <li> Index 2: <tt>[0, 1, 0]</tt>;
+ * <li> Index 3: <tt>[1, 1, 0]</tt>;
+ * <li> Index 4: <tt>[0, 0, 1]</tt>;
+ * <li> Index 5: <tt>[1, 0, 1]</tt>;
+ * <li> Index 6: <tt>[0, 1, 1]</tt>;
+ * <li> Index 7: <tt>[1, 1, 1]</tt>;
+ * <li> Index 8: <tt>[0, 1/2, 0]</tt>;
+ * <li> Index 9: <tt>[1, 1/2, 0]</tt>;
+ * <li> Index 10: <tt>[1/2, 0, 0]</tt>;
+ * <li> Index 11: <tt>[1/2, 1, 0]</tt>;
+ * <li> Index 12: <tt>[0, 1/2, 1]</tt>;
+ * <li> Index 13: <tt>[1, 1/2, 1]</tt>;
+ * <li> Index 14: <tt>[1/2, 0, 1]</tt>;
+ * <li> Index 15: <tt>[1/2, 1, 1]</tt>;
+ * <li> Index 16: <tt>[0, 0, 1/2]</tt>;
+ * <li> Index 17: <tt>[1, 0, 1/2]</tt>;
+ * <li> Index 18: <tt>[0, 1, 1/2]</tt>;
+ * <li> Index 19: <tt>[1, 1, 1/2]</tt>;
+ * <li> Index 20: <tt>[0, 1/2, 1/2]</tt>;
+ * <li> Index 21: <tt>[1, 1/2, 1/2]</tt>;
+ * <li> Index 22: <tt>[1/2, 0, 1/2]</tt>;
+ * <li> Index 23: <tt>[1/2, 1, 1/2]</tt>;
+ * <li> Index 24: <tt>[1/2, 1/2, 0]</tt>;
+ * <li> Index 25: <tt>[1/2, 1/2, 1]</tt>;
* <li> Index 26: <tt>[1/2, 1/2, 1/2]</tt>;
* </ul>
* </ul>
*
* <li> 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.
* </ul>
* <h4>Q4 elements</h4>
* <ul>
*
* <li> 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
* </ul>
- * Note the reverse ordering of degrees of freedom on the left and upper
- * line.
*
* @author Brian Carnes, 2002, Ralf Hartmann 2004, 2005
*/
*/
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:
/**
/*@}*/
-/* -------------- 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 ------------------------- */
* 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 <int dim>
class MappingCartesian : public Mapping<dim>
* <tt>deal.II/doc/reports/mapping_q/index.html</tt> in the `Reports'
* section of `Documentation'.
*
- * @author Ralf Hartmann, Guido Kanschat 2000, 2001
+ * @author Ralf Hartmann, 2000, 2001, 2005; Guido Kanschat 2000, 2001
*/
template <int dim>
class MappingQ : public MappingQ1<dim>
* 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 <int dim>
class MappingQ1 : public Mapping<dim>
const Point<dim> &p,
InternalData &mdata,
Point<dim> &p_unit) const;
-
- /**
- * Mapping between tensor product
- * ordering and real ordering of
- * the vertices.
- */
- static const unsigned int vertex_mapping[1<<dim];
private:
/**
#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<Point<1> > &unit_points,
InternalData& data) const;
*/
virtual void create_triangulation (const std::vector<Point<dim> > &vertices,
const std::vector<CellData<dim> > &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<Point<dim> > &vertices,
+ const std::vector<CellData<dim> > &cells,
+ const SubCellData &subcelldata);
/**
* Writes all refine and coarsen
* 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/
*
* <h4>Implementation conventions for two spatial dimensions</h4>
+ *
+ * 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 <tt>subline(0).vertex(0)==line.vertex(0)</tt> and
* 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 <tt>(1,0)</tt>, vertex 2 at <tt>(1,1)</tt> and vertex 3 at <tt>(0,1)</tt>.
+ * with vertex 0 being the origin of the coordinate system, vertex 1
+ * having coordinates <tt>(1,0)</tt>, vertex 2 at <tt>(0,1)</tt> and
+ * vertex 3 at <tt>(1,1)</tt>.
*
*
* <h3>Implementation conventions for three spatial dimensions</h3>
* 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
* conventions can be extracted at run- or compile-time from the
* member functions and variables of the GeometryInfo classes.
*
- *
* <h4>Vertices</h4>
+ *
+ * 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.
*
* <h4>Lines</h4>
*
- * 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:
+ * <ul>
+ * <li>first the lines of face (z=0) in 2d line ordering,
+ * <li>then the lines of face (z=1) in 2d line ordering,
+ * <li>finally the lines in z direction in lexicographic ordering
+ * </ul>
* @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
* *--->---* *--->---*
* /| | / /|
*
* <h4>Faces</h4>
*
- * 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 <em>standard</em> direction of the faces is determined by the
- * numbers the lines have within a given face. This is like follows:
+ * The <em>standard</em> 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:
* <ul>
* <li> 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
*
- * <li> Faces 2 and 4:
+ * <li> 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
*
- * <li> Faces 3 and 5:
+ * <li> 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
* </ul>
- *
- * 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.
* <ul>
- * <li> Line 0 of face 0, and line 0 of face 2;
- * <li> Line 1 of face 0, and line 3 of face 3;
- * <li> Line 2 of face 0, and line 0 of face 4;
- * <li> Line 3 of face 0, and line 3 of face 5;
- * <li> Line 0 of face 1, and line 2 of face 2;
- * <li> Line 1 of face 1, and line 1 of face 3;
- * <li> Line 2 of face 1, and line 2 of face 4;
- * <li> Line 3 of face 1, and line 1 of face 5;
- * <li> Line 3 of face 2, and line 0 of face 5;
- * <li> Line 1 of face 2, and line 0 of face 3;
- * <li> Line 1 of face 4, and line 2 of face 3;
- * <li> Line 3 of face 4, and line 2 of face 5.
+ * <li> Face 0: lines 8, 10, 0, 4;
+ * <li> Face 1: lines 9, 11, 1, 5;
+ * <li> Face 2: lines 2, 6, 8, 9;
+ * <li> Face 3: lines 3, 7, 10, 11;
+ * <li> Face 4: lines 0, 1, 2, 3;
+ * <li> Face 5: lines 4, 5, 6, 7;
* </ul>
+ * 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
*
* <h4>Children</h4>
*
- * 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 |/
* *-------* *-------*
* Taking into account the orientation of the faces, the following
* children are adjacent to the respective faces:
* <ul>
- * <li> Face 0: children 0, 1, 2, 3;
- * <li> Face 1: children 4, 5, 6, 7;
- * <li> Face 2: children 0, 1, 5, 4;
- * <li> Face 3: children 1, 5, 6, 2;
- * <li> Face 4: children 3, 2, 6, 7;
- * <li> Face 5: children 0, 4, 7, 3.
+ * <li> Face 0: children 0, 2, 4, 6;
+ * <li> Face 1: children 1, 3, 5, 7;
+ * <li> Face 2: children 0, 4, 1, 5;
+ * <li> Face 3: children 2, 6, 3, 7;
+ * <li> Face 4: children 0, 1, 2, 3;
+ * <li> Face 5: children 4, 5, 6, 7.
* </ul>
- * 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 <tt>face_orientation=false</tt> 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
- * <tt>cell->neighbor_child_on_subface(face_no,subface_no)</tt>, 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
- * <tt>GeometryInfo<3>::child_cell_on_face</tt>). In general, the use of
- * <tt>GeometryInfo<3>::child_cell_on_face</tt> 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
+ * <tt>cell->neighbor_child_on_subface(face_no,subface_no)</tt>,
+ * 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
+ * <tt>GeometryInfo<3>::child_cell_on_face</tt>).
*
* <h4>Coordinate systems</h4>
*
* 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 <tt>(x,0,z)</tt>
- * 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
- * <tt>(x,y,0)</tt>. 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):
* <ul>
* <li> Vertex 0: <tt>(0,0,0)</tt>;
* <li> Vertex 1: <tt>(1,0,0)</tt>;
- * <li> Vertex 2: <tt>(1,0,1)</tt>;
- * <li> Vertex 3: <tt>(0,0,1)</tt>;
- * <li> Vertex 4: <tt>(0,1,0)</tt>;
- * <li> Vertex 5: <tt>(1,1,0)</tt>;
- * <li> Vertex 6: <tt>(1,1,1)</tt>;
- * <li> Vertex 7: <tt>(0,1,1)</tt>.
+ * <li> Vertex 2: <tt>(0,1,0)</tt>;
+ * <li> Vertex 3: <tt>(1,1,0)</tt>;
+ * <li> Vertex 4: <tt>(0,0,1)</tt>;
+ * <li> Vertex 5: <tt>(1,0,1)</tt>;
+ * <li> Vertex 6: <tt>(0,1,1)</tt>;
+ * <li> Vertex 7: <tt>(1,1,1)</tt>.
* </ul>
*
*
* data stored in the triangulation.
*
* @ingroup grid
- * @author Wolfgang Bangerth, 1998
+ * @author Wolfgang Bangerth, 1998; Ralf Hartmann, 2005
*/
template <int dim>
class Triangulation : public Subscriptor
*/
virtual void create_triangulation (const std::vector<Point<dim> > &vertices,
const std::vector<CellData<dim> > &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<Point<dim> > &vertices,
+ const std::vector<CellData<dim> > &cells,
+ const SubCellData &subcelldata);
/**
* Distort the grid by randomly
// 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]) ?
// 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]) ?
// 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
= face->child(0)->line(1)->dof_index(dof);
for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
dofs_on_children[next_index++]
- = face->child(1)->line(2)->dof_index(dof);
+ = face->child(2)->line(1)->dof_index(dof);
for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
dofs_on_children[next_index++]
- = face->child(2)->line(3)->dof_index(dof);
+ = face->child(0)->line(3)->dof_index(dof);
for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
dofs_on_children[next_index++]
- = face->child(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)
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;
};
{
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;
};
{
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
// 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
// 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
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;
// 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]);
};
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
// 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} },
{ { {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]);
};
for (unsigned int c=0; c<GeometryInfo<dim>::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;
};
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.;
// 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]);
};
// 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])
&&
{
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
{
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
{
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
// 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<Point<dim-1> > constraint_points;
// Add midpoint
// Add midpoints of lines of
// "mother-face"
- constraint_points.push_back (Point<dim-1> (0.5, 0));
+ constraint_points.push_back (Point<dim-1> (0, 0.5));
constraint_points.push_back (Point<dim-1> (1, 0.5));
+ constraint_points.push_back (Point<dim-1> (0.5, 0));
constraint_points.push_back (Point<dim-1> (0.5, 1));
- constraint_points.push_back (Point<dim-1> (0, 0.5));
if (this->degree>1)
{
// Add nodes of lines interior
// in the "mother-face"
- // line 5: use line 15
- QProjector<dim-1>::project_to_subface(qline, 3, 0, p_line);
+ // line 5: use line 9
+ QProjector<dim-1>::project_to_subface(qline, 0, 0, p_line);
for (unsigned int i=0; i<n; ++i)
constraint_points.push_back (p_line[i] + Point<dim-1> (0.5, 0));
// line 6: use line 10
QProjector<dim-1>::project_to_subface(qline, 0, 1, p_line);
- for (unsigned int i=0; i<n; ++i)
- constraint_points.push_back (p_line[i] + Point<dim-1> (0, 0.5));
- // line 7: use line 16
- QProjector<dim-1>::project_to_subface(qline, 3, 1, p_line);
for (unsigned int i=0; i<n; ++i)
constraint_points.push_back (p_line[i] + Point<dim-1> (0.5, 0));
- // line 8: use line 9
- QProjector<dim-1>::project_to_subface(qline, 0, 0, p_line);
+ // line 7: use line 13
+ QProjector<dim-1>::project_to_subface(qline, 2, 0, p_line);
for (unsigned int i=0; i<n; ++i)
- constraint_points.push_back (p_line[i] + Point<dim-1> (0, 0.5));
+ constraint_points.push_back (p_line[i] + Point<dim-1> (0, 0.5));
+ // line 8: use line 14
+ QProjector<dim-1>::project_to_subface(qline, 2, 1, p_line);
+ for (unsigned int i=0; i<n; ++i)
+ constraint_points.push_back (p_line[i] + Point<dim-1> (0, 0.5));
// DoFs on bordering lines
// lines 9-16
// 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++)
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) *
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
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);
}
}
{
for (unsigned int c=0; c<GeometryInfo<2>::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,
{
for (unsigned int c=0; c<GeometryInfo<3>::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,
// 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
// 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.dofs_per_line; ++i)
- h2l[next_index++] = 2+i;
+ h2l[next_index++] = (2+i)*n;
// right line
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
h2l[next_index++] = (2+i)*n+1;
+ // bottom line
+ for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+ h2l[next_index++] = 2+i;
// top line
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
h2l[next_index++] = n+2+i;
- // left line
- for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- h2l[next_index++] = (2+i)*n;
// inside quad
Assert (fe.dofs_per_quad == fe.dofs_per_line*fe.dofs_per_line,
ExcInternalError());
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.dofs_per_line; ++i)
- h2l[next_index++] = 2+i;
+ // bottom face
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- h2l[next_index++] = (2+i)*n2+1;
+ h2l[next_index++] = (2+i)*n;
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- h2l[next_index++] = n2+2+i;
+ h2l[next_index++] = (2+i)*n+1;
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- h2l[next_index++] = (2+i)*n2;
- // back face, counterclock wise
+ h2l[next_index++] = 2+i;
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
h2l[next_index++] = n+2+i;
+ // top face
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- h2l[next_index++] = (2+i)*n2+n+1;
+ h2l[next_index++] = n2+(2+i)*n;
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- h2l[next_index++] = n2+n+2+i;
+ h2l[next_index++] = n2+(2+i)*n+1;
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- h2l[next_index++] = (2+i)*n2+n;
- // lines in y-direction,
- // counterclock wise
+ h2l[next_index++] = n2+2+i;
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- h2l[next_index++] = (2+i)*n;
+ h2l[next_index++] = n2+n+2+i;
+ // lines in z-direction
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- h2l[next_index++] = (2+i)*n+1;
+ h2l[next_index++] = (2+i)*n2;
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- h2l[next_index++] = n2+(2+i)*n+1;
+ h2l[next_index++] = (2+i)*n2+1;
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- h2l[next_index++] = n2+(2+i)*n;
+ h2l[next_index++] = (2+i)*n2+n;
+ for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+ h2l[next_index++] = (2+i)*n2+n+1;
// inside quads
Assert (fe.dofs_per_quad == fe.dofs_per_line*fe.dofs_per_line,
ExcInternalError());
+ // left face
+ for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+ for (unsigned int j=0; j<fe.dofs_per_line; ++j)
+ h2l[next_index++] = (2+i)*n2+(2+j)*n;
+ // right face
+ for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+ for (unsigned int j=0; j<fe.dofs_per_line; ++j)
+ h2l[next_index++] = (2+i)*n2+(2+j)*n+1;
// front face
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
for (unsigned int j=0; j<fe.dofs_per_line; ++j)
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
for (unsigned int j=0; j<fe.dofs_per_line; ++j)
h2l[next_index++] = (2+i)*n+2+j;
- // right face
- for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- for (unsigned int j=0; j<fe.dofs_per_line; ++j)
- h2l[next_index++] = (2+i)*n2+(2+j)*n+1;
// top face
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
for (unsigned int j=0; j<fe.dofs_per_line; ++j)
h2l[next_index++] = n2+(2+i)*n+2+j;
- // left face
- for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- for (unsigned int j=0; j<fe.dofs_per_line; ++j)
- h2l[next_index++] = (2+i)*n2+(2+j)*n;
// inside hex
Assert (fe.dofs_per_hex == fe.dofs_per_quad*fe.dofs_per_line,
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 <int dim>
bool
FE_Q_Hierarchical<dim>::has_support_on_face (const unsigned int shape_index,
Assert (face_index < GeometryInfo<dim>::faces_per_cell,
ExcIndexRange (face_index, 0, GeometryInfo<dim>::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)))
// 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<dim>::vertices_per_cell,
ExcInternalError());
static_cast<signed int>(GeometryInfo<dim>::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,
if (dim == 3)
return (quad_index == face_index);
else
- Assert (false, ExcNotImplemented());
+ Assert (false, ExcNotImplemented());
}
else
// dof on hex
// 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.dofs_per_line; ++i)
- h2l[next_index++] = 1+i;
- // second line
+ h2l[next_index++] = (1+i)*n;
+
+ // right line
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
h2l[next_index++] = (2+i)*n-1;
- // third line
+
+ // bottom line
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- h2l[next_index++] = n*(n-1)+i+1;
- // fourth line
+ h2l[next_index++] = 1+i;
+
+ // top line
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- h2l[next_index++] = (1+i)*n;
+ h2l[next_index++] = n*(n-1)+i+1;
+
// inside quad
Assert (fe.dofs_per_quad == fe.dofs_per_line*fe.dofs_per_line,
ExcInternalError());
{
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.dofs_per_line; ++i)
- h2l[next_index++] = 1+i;
+ h2l[next_index++] = (i+1)*n;
+ // line 1
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- h2l[next_index++] = n-1+(i+1)*n*n;
+ h2l[next_index++] = n-1+(i+1)*n;
+ // line 2
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- h2l[next_index++] = n*n*(n-1)+i+1;
+ h2l[next_index++] = 1+i;
+ // line 3
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- h2l[next_index++] = (i+1)*n*n;
+ h2l[next_index++] = 1+i+n*(n-1);
+ // line 4
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- h2l[next_index++] = 1+i+n*(n-1);
+ h2l[next_index++] = (n-1)*n*n+(i+1)*n;
+ // line 5
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- h2l[next_index++] = n-1+(i+1)*n*n+n*(n-1);
+ h2l[next_index++] = (n-1)*(n*n+1)+(i+1)*n;
+ // line 6
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- h2l[next_index++] = n*n*(n-1)+i+1+n*(n-1);
+ h2l[next_index++] = n*n*(n-1)+i+1;
+ // line 7
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- h2l[next_index++] = (i+1)*n*n+n*(n-1);
+ h2l[next_index++] = n*n*(n-1)+i+1+n*(n-1);
+ // line 8
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- h2l[next_index++] = (i+1)*n;
+ h2l[next_index++] = (i+1)*n*n;
+ // line 9
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- h2l[next_index++] = n-1+(i+1)*n;
+ h2l[next_index++] = n-1+(i+1)*n*n;
+ // line 10
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- h2l[next_index++] = (n-1)*(n*n+1)+(i+1)*n;
+ h2l[next_index++] = (i+1)*n*n+n*(n-1);
+ // line 11
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- h2l[next_index++] = (n-1)*n*n+(i+1)*n;
+ h2l[next_index++] = n-1+(i+1)*n*n+n*(n-1);
+
// inside quads
Assert (fe.dofs_per_quad == fe.dofs_per_line*fe.dofs_per_line,
ExcInternalError());
- // quad 1
+ // face 0
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
for (unsigned int j=0; j<fe.dofs_per_line; ++j)
- h2l[next_index++] = (i+1)*n*n+j+1;
- // quad 2
+ h2l[next_index++] = (i+1)*n*n+n*(j+1);
+ // face 1
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
for (unsigned int j=0; j<fe.dofs_per_line; ++j)
- h2l[next_index++] = (i+1)*n*n+n*(n-1)+j+1;
- // quad 3
+ h2l[next_index++] = (i+1)*n*n+n-1+n*(j+1);
+ // face 2, note the orientation!
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
for (unsigned int j=0; j<fe.dofs_per_line; ++j)
- h2l[next_index++] = n*(i+1)+j+1;
- // quad 4
+ h2l[next_index++] = (j+1)*n*n+i+1;
+ // face 3, note the orientation!
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
for (unsigned int j=0; j<fe.dofs_per_line; ++j)
- h2l[next_index++] = (i+1)*n*n+n-1+n*(j+1);
- // quad 5
+ h2l[next_index++] = (j+1)*n*n+n*(n-1)+i+1;
+ // face 4
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
for (unsigned int j=0; j<fe.dofs_per_line; ++j)
- h2l[next_index++] = (n-1)*n*n+n*(i+1)+j+1;
- // quad 6
+ h2l[next_index++] = n*(i+1)+j+1;
+ // face 5
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
for (unsigned int j=0; j<fe.dofs_per_line; ++j)
- h2l[next_index++] = (i+1)*n*n+n*(j+1);
+ h2l[next_index++] = (n-1)*n*n+n*(i+1)+j+1;
// inside hex
Assert (fe.dofs_per_hex == fe.dofs_per_quad*fe.dofs_per_line,
-
template <int dim>
void
FETools::lexicographic_to_hierarchic_numbering (const FiniteElementData<dim> &fe,
#include <fe/mapping.h>
-#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 <int dim>
Mapping<dim>::~Mapping ()
{}
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());
{
static const Point<dim>
normals[GeometryInfo<2>::faces_per_cell]
- = { Point<dim>(0, -1),
- Point<dim>(1, 0),
- Point<dim>(0, 1),
- Point<dim>(-1, 0) };
+ = { Point<dim>(-1, 0),
+ Point<dim>( 1, 0),
+ Point<dim>( 0,-1),
+ Point<dim>( 0, 1) };
std::fill (normal_vectors.begin(),
normal_vectors.end(),
normals[face_no]);
{
static const Point<dim>
normals[GeometryInfo<3>::faces_per_cell]
- = { Point<dim>(0, -1, 0),
- Point<dim>(0, 1, 0),
- Point<dim>(0, 0, -1),
- Point<dim>(1, 0, 0),
- Point<dim>(0, 0, 1),
- Point<dim>(-1, 0, 0) };
+ = { Point<dim>(-1, 0, 0),
+ Point<dim>( 1, 0, 0),
+ Point<dim>( 0,-1, 0),
+ Point<dim>( 0, 1, 0),
+ Point<dim>( 0, 0,-1),
+ Point<dim>( 0, 0, 1) };
std::fill (normal_vectors.begin(),
normal_vectors.end(),
normals[face_no]);
// product of the local lengths
// since the jacobian is diagonal
double J = 1.;
- for (unsigned int d=0;d<dim;++d)
- if (d != (this->normal_directions[face_no]/2))
+ for (unsigned int d=0; d<dim; ++d)
+ if (d != GeometryInfo<dim>::unit_normal_direction[face_no])
J *= data.length[d];
if (data.current_update_flags() & update_JxW_values)
// product of the local lengths
// since the jacobian is diagonal
double J = 1.;
- for (unsigned int d=0;d<dim;++d)
- if (d != (this->normal_directions[face_no]/2))
+ for (unsigned int d=0; d<dim; ++d)
+ if (d != GeometryInfo<dim>::unit_normal_direction[face_no])
J *= data.length[d];
if (data.current_update_flags() & update_JxW_values)
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());
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());
{
// (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());
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]) ==
a.resize(GeometryInfo<dim>::vertices_per_cell);
for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
- a[i] = cell->vertex(this->vertex_mapping[i]);
+ a[i] = cell->vertex(i);
}
}
break;
case 3:
+ {
// in 3d also add the
// points located on
// the boundary faces
add_quad_support_points (cell, a);
apply_laplace_vector (laplace_on_hex_vector, a);
break;
-
+ }
default:
Assert(false, ExcNotImplemented());
break;
#if (deal_II_dimension == 1)
-template<>
-const unsigned int MappingQ1<deal_II_dimension>::vertex_mapping[2] =
-{ 0, 1 };
-
template<>
void
#if (deal_II_dimension == 2)
-template<> const unsigned int
-MappingQ1<2>::vertex_mapping[4] =
-{ 0, 1, 3, 2 };
-
template<>
void
#if (deal_II_dimension == 3)
-template<>
-const unsigned int MappingQ1<deal_II_dimension>::vertex_mapping[8] =
-{ 0, 1, 4, 5, 3, 2, 7, 6 };
-
template<>
void
const unsigned int nfaces = GeometryInfo<dim>::faces_per_cell;
data.unit_tangentials.resize (nfaces*(dim-1),
std::vector<Tensor<1,dim> > (n_original_q_points));
- for (unsigned int i=0; i<nfaces; ++i)
+ if (dim==2)
{
- // Base index of the
- // non-zero entry of the
- // tangential vector. Add
- // dim so we can subtract 1
- // without getting negative
- // values.
- unsigned int nindex = this->normal_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<nfaces; ++i)
+ {
+ Tensor<1,dim> 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<nfaces; ++i)
{
- // Second tangential
- // has a non-zero in
- // component (i-1)%dim,
- // if normal is
- // non-zero in
- // i. Creates a
- // right-handed system.
- Tensor<1,dim> 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<dim>::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<dim>::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);
}
}
}
a.resize(GeometryInfo<dim>::vertices_per_cell);
for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
- a[i] = cell->vertex(vertex_mapping[i]);
+ a[i] = cell->vertex(i);
}
// are always numbered first,
// we can access them easily
for (unsigned int j=0; j<dim; ++j)
- shift_vector[j] = mapping_values(this->vertex_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;
}
}
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)
// there is only one cell, so
// simple task
const typename Triangulation<dim>::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<GeometryInfo<dim>::faces_per_cell; ++f)
+ cell->face(f)->set_boundary_indicator (f);
}
#endif
std::vector<Point<dim> > vertices (GeometryInfo<dim>::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<CellData<dim> > cells (1);
for (unsigned int i=0;i<GeometryInfo<dim>::vertices_per_cell;++i)
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;
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;
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;
{
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;
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<CellData<2> > cells (4, CellData<2>());
for (unsigned int i=0; i<4; ++i)
{
cells[i].vertices[j] = cell_vertices[i][j];
cells[i].material_id = 0;
};
- tria.create_triangulation (std::vector<Point<2> >(&vertices[0], &vertices[10]),
- cells,
- SubCellData()); // no boundary information
+ tria.create_triangulation (
+ std::vector<Point<2> >(&vertices[0], &vertices[10]),
+ cells,
+ SubCellData()); // no boundary information
if (colorize)
{
Point<dim> (b,(a+b)/2),
Point<dim> (a,b),
Point<dim> ((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<CellData<2> > cells (3, CellData<2>());
cells[i].material_id = 0;
};
- tria.create_triangulation (std::vector<Point<dim> >(&vertices[0], &vertices[8]),
- cells,
- SubCellData()); // no boundary information
+ tria.create_triangulation (
+ std::vector<Point<dim> >(&vertices[0], &vertices[8]),
+ cells,
+ SubCellData()); // no boundary information
}
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<CellData<2> > cells (5, CellData<2>());
cells[i].material_id = 0;
};
- tria.create_triangulation (std::vector<Point<2> >(&vertices[0], &vertices[8]),
- cells,
- SubCellData()); // no boundary information
+ tria.create_triangulation (
+ std::vector<Point<2> >(&vertices[0], &vertices[8]),
+ cells,
+ SubCellData()); // no boundary information
}
{
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());
}
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<CellData<2> > cells (4, CellData<2>());
cells[i].material_id = 0;
};
- tria.create_triangulation (std::vector<Point<2> >(&vertices[0], &vertices[8]),
- cells,
- SubCellData()); // no boundary information
+ tria.create_triangulation (
+ std::vector<Point<2> >(&vertices[0], &vertices[8]),
+ cells,
+ SubCellData()); // no boundary information
}
{
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;
};
{
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
}
cells[i].vertices[j] = cell_vertices[i][j];
cells[i].material_id = 0;
};
-
- tria.create_triangulation (std::vector<Point<dim> >(&vertices[0], &vertices[26]),
- cells,
- SubCellData()); // no boundary information
+
+ SubCellData subcelldata;
+ tria.create_triangulation_compatibility (
+ std::vector<Point<dim> >(&vertices[0], &vertices[26]),
+ cells,
+ subcelldata); // no boundary information
}
};
GridReordering<3>::reorder_cells (cells);
-
- tria.create_triangulation (std::vector<Point<3> >(&vertices[0], &vertices[n_vertices]),
- cells,
- SubCellData()); // no boundary information
+
+ SubCellData subcelldata;
+ tria.create_triangulation_compatibility (
+ std::vector<Point<3> >(&vertices[0], &vertices[n_vertices]),
+ cells,
+ subcelldata); // no boundary information
}
};
GridReordering<3>::reorder_cells (cells);
-
- tria.create_triangulation (std::vector<Point<3> >(&vertices[0], &vertices[24]),
- cells,
- SubCellData()); // no boundary information
+
+ SubCellData subcelldata;
+ tria.create_triangulation_compatibility (
+ std::vector<Point<3> >(&vertices[0], &vertices[24]),
+ cells,
+ subcelldata); // no boundary information
Triangulation<3>::cell_iterator cell = tria.begin();
Triangulation<3>::cell_iterator end = tria.end();
{
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());
}
// ... and cells
GridReordering<dim>::invert_all_cells_of_negative_grid (vertices, cells);
GridReordering<dim>::reorder_cells (cells);
- tria->create_triangulation (vertices, cells, subcelldata);
+ tria->create_triangulation_compatibility (vertices, cells, subcelldata);
}
// ...and cells
GridReordering<dim>::invert_all_cells_of_negative_grid (vertices, cells);
GridReordering<dim>::reorder_cells (cells);
- tria->create_triangulation (vertices, cells, subcelldata);
+ tria->create_triangulation_compatibility (vertices, cells, subcelldata);
}
// ... 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
// ... 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
// ... and cells
GridReordering<dim>::invert_all_cells_of_negative_grid (vertices, cells);
GridReordering<dim>::reorder_cells (cells);
- tria->create_triangulation (vertices, cells, subcelldata);
+ tria->create_triangulation_compatibility (vertices, cells, subcelldata);
}
SubCellData subcelldata;
GridTools::delete_unused_vertices(vertices, cells, subcelldata);
GridReordering<dim>::reorder_cells (cells);
- tria->create_triangulation (vertices, cells, subcelldata);
+ tria->create_triangulation_compatibility (vertices, cells, subcelldata);
#endif
}
if (switch_y_z)
GridReordering<dim>::invert_all_cells_of_negative_grid (vertices, cells);
GridReordering<dim>::reorder_cells (cells);
- tria->create_triangulation (vertices, cells, subcelldata);
+ tria->create_triangulation_compatibility (vertices, cells, subcelldata);
#endif
}
// 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<GeometryInfo<dim>::vertices_per_cell;
++vertex)
- out << cell->vertex_index(vertex)+1 << ' ';
+ out << cell->vertex_index(GeometryInfo<dim>::ucd_to_deal[vertex])+1 << ' ';
out << std::endl;
};
// (dots/inch)
for (unsigned int k=0;k<=nv;++k)
{
- const Point<dim>& p = cell->vertex(k % nv);
+ const Point<dim>& p = cell->vertex(
+ GeometryInfo<dim>::ucd_to_deal[k % nv]);
for (unsigned int d=0;d<dim;++d)
{
int val = (int)(1200 * xfig_flags.scaling(d) *
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<nf;++f)
{
typename Triangulation<dim>::face_iterator
- face = cell->face(f);
+ face = cell->face(face_reorder[f]);
const unsigned char bi = face->boundary_indicator();
if (bi != 255)
{
};
// note: vertex numbers are 1-base
for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
- out << face->vertex_index(vertex)+1 << ' ';
+ out << face->vertex_index(GeometryInfo<dim-1>::ucd_to_deal[vertex])+1 << ' ';
out << std::endl;
++index;
// drawing pencil at the
// end
for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
- out << cell->vertex(i)
+ out << cell->vertex(GeometryInfo<dim>::ucd_to_deal[i])
<< ' ' << cell->level()
<< ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl;
out << cell->vertex(0)
out << "# cell " << cell << std::endl;
if (mapping==0 || n_points==2 || !cell->has_boundary_lines())
- {
+ {
// front face
out << cell->vertex(0)
<< ' ' << cell->level()
<< cell->vertex(1)
<< ' ' << cell->level()
<< ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << cell->vertex(2)
+ << cell->vertex(5)
<< ' ' << cell->level()
<< ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << cell->vertex(3)
+ << cell->vertex(4)
<< ' ' << cell->level()
<< ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
<< cell->vertex(0)
<< ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
<< std::endl;
// back face
- out << cell->vertex(4)
+ out << cell->vertex(2)
<< ' ' << cell->level()
<< ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << cell->vertex(5)
+ << cell->vertex(3)
<< ' ' << cell->level()
<< ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << cell->vertex(6)
+ << cell->vertex(7)
<< ' ' << cell->level()
<< ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << cell->vertex(7)
+ << cell->vertex(6)
<< ' ' << cell->level()
<< ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << cell->vertex(4)
+ << cell->vertex(2)
<< ' ' << cell->level()
<< ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
<< std::endl;
out << cell->vertex(0)
<< ' ' << cell->level()
<< ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << cell->vertex(4)
+ << cell->vertex(2)
<< ' ' << cell->level()
<< ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
<< std::endl;
out << cell->vertex(1)
<< ' ' << cell->level()
<< ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << cell->vertex(5)
+ << cell->vertex(3)
<< ' ' << cell->level()
<< ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
<< std::endl;
- out << cell->vertex(2)
+ out << cell->vertex(5)
<< ' ' << cell->level()
<< ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << cell->vertex(6)
+ << cell->vertex(7)
<< ' ' << cell->level()
<< ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
<< std::endl;
- out << cell->vertex(3)
+ out << cell->vertex(4)
<< ' ' << cell->level()
<< ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << cell->vertex(7)
+ << cell->vertex(6)
<< ' ' << cell->level()
<< ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
<< std::endl;
const std::vector<Point<3> > &all_vertices,
std::vector<CellData<3> > &cells)
{
+ int vertices_lex[GeometryInfo<3>::vertices_per_cell];
unsigned int n_negative_cells=0;
for (unsigned int cell_no=0; cell_no<cells.size(); ++cell_no)
- if (GridTools::cell_measure(all_vertices, cells[cell_no].vertices) < 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
- AssertThrow(GridTools::cell_measure(all_vertices, cells[cell_no].vertices) > 0,
- ExcInternalError());
- }
-
+ {
+ // GridTools::cell_measure
+ // requires the vertices to be
+ // in lexicographic ordering
+ for (unsigned int i=0; i<GeometryInfo<3>::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<GeometryInfo<3>::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
GridTools::cell_measure(const std::vector<Point<3> > &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),
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
}
+template <int dim>
+void
+PersistentTriangulation<dim>::create_triangulation_compatibility (
+ const std::vector<Point<dim> > &,
+ const std::vector<CellData<dim> > &,
+ const SubCellData &)
+{
+ Assert (false, ExcImpossibleInDim(dim));
+}
+
+
template <int dim>
void
PersistentTriangulation<dim>::write_flags(std::ostream &out) const
#if deal_II_dimension == 1
+template <>
+void Triangulation<1>::create_triangulation_compatibility (
+ const std::vector<Point<1> > &v,
+ const std::vector<CellData<1> > &cells,
+ const SubCellData &subcelldata)
+{
+ // no reordering required
+ create_triangulation(v, cells, subcelldata);
+}
+
+
template <>
void Triangulation<1>::create_triangulation (const std::vector<Point<1> > &v,
const std::vector<CellData<1> > &cells,
#if deal_II_dimension == 2
+template <>
+void Triangulation<2>::create_triangulation_compatibility (
+ const std::vector<Point<2> > &v,
+ const std::vector<CellData<2> > &c,
+ const SubCellData &subcelldata)
+{
+ const bool output=false;
+
+ // copy cell data as it will be
+ // modified in the following
+ std::vector<CellData<2> > cells(c);
+
+ for (unsigned int cell=0; cell<cells.size(); ++cell)
+ {
+ if (output)
+ {
+ std::cout << "cell" << cell << ": ";
+ for (unsigned int vertex=0; vertex<4; ++vertex)
+ std::cout << cells[cell].vertices[vertex] << " ";
+ }
+
+ std::swap(cells[cell].vertices[2],cells[cell].vertices[3]);
+
+
+ if (output)
+ {
+ std::cout << " after swap: ";
+ for (unsigned int vertex=0; vertex<4; ++vertex)
+ std::cout << cells[cell].vertices[vertex] << " ";
+ std::cout << std::endl;
+ }
+ }
+
+ // subcelldata.boundary_lines refer
+ // to the vertices in v. they don't
+ // need to be changed
+
+ create_triangulation(v, cells, subcelldata);
+}
+
template <>
void Triangulation<2>::create_triangulation (const std::vector<Point<2> > &v,
const std::vector<CellData<2> > &cells,
};
+template <>
+void
+Triangulation<3>::create_triangulation_compatibility (
+ const std::vector<Point<3> > &v,
+ const std::vector<CellData<3> > &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<CellData<3> > cells(c);
+ SubCellData subcelldata(subc);
+
+ unsigned int tmp[GeometryInfo<3>::vertices_per_cell];
+ for (unsigned int cell=0; cell<cells.size(); ++cell)
+ {
+ if (output1)
+ {
+ std::cout << "cell" << cell << ": ";
+ for (unsigned int vertex=0;
+ vertex<GeometryInfo<3>::vertices_per_cell; ++vertex)
+ std::cout << cells[cell].vertices[vertex] << " ";
+ }
+
+ for (unsigned int i=0; i<GeometryInfo<3>::vertices_per_cell; ++i)
+ tmp[i]=cells[cell].vertices[i];
+ for (unsigned int i=0; i<GeometryInfo<3>::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<GeometryInfo<3>::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<GeometryInfo<3>::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<vertices.size(); ++i)
+ std::cout << vertices[i] << std::endl;
+ }
+ // now points in boundary quads
+ std::vector<CellData<2> >::iterator boundary_quad
+ = subcelldata.boundary_quads.begin();
+ std::vector<CellData<2> >::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<Point<3> > &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
// 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();
}
// 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];
// 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
line->set_boundary_indicator (boundary_line->material_id);
}
+
// now go on with boundary faces
std::vector<CellData<2> >::const_iterator boundary_quad
= subcelldata.boundary_quads.begin();
// 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
// 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 <algorithms>
+ 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;
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
The created children are numbered like this:
.--.--.
- |3 . 2|
+ |2 . 3|
.--.--.
|0 | 1|
.--.--.
// 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};
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)
{
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<n_children; ++i)
// 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
// 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]
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)
{
// 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(),
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();
// orientation
for (unsigned int face=0;
face<GeometryInfo<dim>::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
// 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]));
// 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
// 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),
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
// 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();
//
// 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
// 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
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)
};
// 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)
// 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<dim>::child_cell_on_face(0,0)==0, ExcInternalError());
+ Assert(GeometryInfo<dim>::child_cell_on_face(2,0)==0, ExcInternalError());
+ Assert(GeometryInfo<dim>::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<dim>::child_cell_on_face(1,0)==1, ExcInternalError());
+ Assert(GeometryInfo<dim>::child_cell_on_face(2,2)==1, ExcInternalError());
+ Assert(GeometryInfo<dim>::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<dim>::child_cell_on_face(0,1)==2, ExcInternalError());
+ Assert(GeometryInfo<dim>::child_cell_on_face(3,0)==2, ExcInternalError());
+ Assert(GeometryInfo<dim>::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<dim>::child_cell_on_face(1,1)==3, ExcInternalError());
+ Assert(GeometryInfo<dim>::child_cell_on_face(3,2)==3, ExcInternalError());
+ Assert(GeometryInfo<dim>::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<dim>::child_cell_on_face(0,2)==4, ExcInternalError());
+ Assert(GeometryInfo<dim>::child_cell_on_face(2,1)==4, ExcInternalError());
+ Assert(GeometryInfo<dim>::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<dim>::child_cell_on_face(1,2)==5, ExcInternalError());
+ Assert(GeometryInfo<dim>::child_cell_on_face(2,3)==5, ExcInternalError());
+ Assert(GeometryInfo<dim>::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<dim>::child_cell_on_face(0,3)==6, ExcInternalError());
+ Assert(GeometryInfo<dim>::child_cell_on_face(3,1)==6, ExcInternalError());
+ Assert(GeometryInfo<dim>::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<dim>::child_cell_on_face(1,3)==7, ExcInternalError());
+ Assert(GeometryInfo<dim>::child_cell_on_face(3,3)==7, ExcInternalError());
+ Assert(GeometryInfo<dim>::child_cell_on_face(5,3)==7, ExcInternalError());
// now we need to set the
// 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
// 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
// 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
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
{
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);
}
template <int dim>
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()));
}
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):
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);
}
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):
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);
}
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]);
}
}
template <int dim>
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()) ));
}
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;
// 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
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
for (unsigned int j=0; j<m; ++j, x+=ds)
points[i*m+j]=((1-x) * vertices[0] +
x * vertices[1]) * (1-y) +
- ((1-x) * vertices[3] +
- x * vertices[2]) * y;
+ ((1-x) * vertices[2] +
+ x * vertices[3]) * y;
}
}
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<vertices_per_face; ++vertex)
{
// first define the two tangent
// 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
unsigned int m=static_cast<unsigned int> (std::sqrt(static_cast<double>(points.size())));
Assert(points.size()==m*m, ExcInternalError());
- std::vector<Point<3> > lp3(m);
+ std::vector<Point<3> > lp0(m);
std::vector<Point<3> > 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<Point<3> > lps(m);
for (unsigned int i=0; i<m; ++i)
{
- get_intermediate_points_between_points(lp3[i], lp1[i], lps);
+ get_intermediate_points_between_points(lp0[i], lp1[i], lps);
for (unsigned int j=0; j<m; ++j)
points[i*m+j]=lps[j];
unsigned int m=static_cast<unsigned int> (std::sqrt(static_cast<double>(points.size())));
Assert(points.size()==m*m, ExcInternalError());
- std::vector<Point<3> > lp3(m);
+ std::vector<Point<3> > lp0(m);
std::vector<Point<3> > 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<Point<3> > lps(m);
for (unsigned int i=0; i<m; ++i)
{
- get_intermediate_points_between_points(lp3[i], lp1[i], lps);
+ get_intermediate_points_between_points(lp0[i], lp1[i], lps);
for (unsigned int j=0; j<m; ++j)
points[i*m+j]=lps[j];
{
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<GeometryInfo<dim>::vertices_per_cell; ++vertex)
- patch->vertices[vertex] = data.mapping->transform_unit_to_real_cell
- (cell, GeometryInfo<dim>::unit_cell_vertex (vertex));
+ patch->vertices[vertex] = data.mapping->transform_unit_to_real_cell
+ (cell, GeometryInfo<dim>::unit_cell_vertex (
+ GeometryInfo<dim>::ucd_to_deal[vertex]));
if (data.n_datasets > 0)
{
// 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
// cell to both the bottom
// and top surfaces:
//
- // 3-------2
+ // 2-------3
// / /
// / /
// / /
// 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<dim+1>(cell->vertex(0)(0),
cell->vertex(0)(1),
cell->vertex(0)(1),
parameter);
- patch->vertices[4] = Point<dim+1>(cell->vertex(3)(0),
- cell->vertex(3)(1),
- parameter-parameter_step);
- patch->vertices[5] = Point<dim+1>(cell->vertex(2)(0),
+ patch->vertices[4] = Point<dim+1>(cell->vertex(2)(0),
cell->vertex(2)(1),
parameter-parameter_step);
- patch->vertices[6] = Point<dim+1>(cell->vertex(2)(0),
- cell->vertex(2)(1),
- parameter);
- patch->vertices[7] = Point<dim+1>(cell->vertex(3)(0),
+ patch->vertices[5] = Point<dim+1>(cell->vertex(3)(0),
cell->vertex(3)(1),
+ parameter-parameter_step);
+ patch->vertices[6] = Point<dim+1>(cell->vertex(3)(0),
+ cell->vertex(3)(1),
+ parameter);
+ patch->vertices[7] = Point<dim+1>(cell->vertex(2)(0),
+ cell->vertex(2)(1),
parameter);
break;
{
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[] =
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
};
}
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,
= (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,
= (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,
// 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());
}
// 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());
}
// 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());
}
cells[0].material_id = 0;
cells[1].material_id = 0;
- tria.create_triangulation (std::vector<Point<3> >(&vertices[0], &vertices[12]),
+ tria.create_triangulation_compatibility (
+ std::vector<Point<3> >(&vertices[0], &vertices[12]),
cells,
SubCellData()); // no boundary information
cells[0].material_id = 0;
cells[1].material_id = 0;
- tria.create_triangulation (std::vector<Point<3> >(&vertices[0], &vertices[12]),
+ tria.create_triangulation_compatibility (
+ std::vector<Point<3> >(&vertices[0], &vertices[12]),
cells,
SubCellData()); // no boundary information
cells[0].material_id = 0;
cells[1].material_id = 0;
- tria.create_triangulation (std::vector<Point<3> >(&vertices[0], &vertices[12]),
+ tria.create_triangulation_compatibility (
+ std::vector<Point<3> >(&vertices[0], &vertices[12]),
cells,
SubCellData()); // no boundary information
cells[2].material_id = 0;
cells[3].material_id = 0;
- tria.create_triangulation (std::vector<Point<3> >(&vertices[0], &vertices[18]),
+ tria.create_triangulation_compatibility (
+ std::vector<Point<3> >(&vertices[0], &vertices[18]),
cells,
SubCellData()); // no boundary information
{
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
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);
(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
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 ();
};
{
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);
};
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;
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.);
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);
#include <fstream>
+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 <int dim>
void test()
{
for(unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
{
for(unsigned int i=0; i<GeometryInfo<dim>::vertices_per_face; ++i)
- subcelldata.boundary_quads[f].vertices[i]=
- GeometryInfo<dim>::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<dim> tria;
- tria.create_triangulation (vertices, cells, subcelldata);
+ tria.create_triangulation_compatibility (vertices, cells, subcelldata);
GridOutFlags::Ucd ucd_flags(true,true);
GridOut grid_out;
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<CellData<2> > 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<Point<2> >(&vertices[0], &vertices[4]),
- cells,
- SubCellData());
+ cells,
+ SubCellData());
exact_areas.push_back(9.);
}
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;
}
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)));
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;
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));
}
// 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<dim> dof(tr);
dof.distribute_dofs(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:
hierarchic_to_lexicographic_numbering[i] = i-1;
break;
- };
+ }
case 2:
{
// 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.dofs_per_line; ++i)
- hierarchic_to_lexicographic_numbering[next_index++] = 1+i;
- // second line
+ hierarchic_to_lexicographic_numbering[next_index++] = (1+i)*n;
+
+ // right line
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
hierarchic_to_lexicographic_numbering[next_index++] = (2+i)*n-1;
- // third line
+
+ // bottom line
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- hierarchic_to_lexicographic_numbering[next_index++] = n*(n-1)+i+1;
- // fourth line
+ hierarchic_to_lexicographic_numbering[next_index++] = 1+i;
+
+ // top line
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- hierarchic_to_lexicographic_numbering[next_index++] = (1+i)*n;
+ hierarchic_to_lexicographic_numbering[next_index++] = n*(n-1)+i+1;
+
// inside quad
Assert (fe.dofs_per_quad == fe.dofs_per_line*fe.dofs_per_line,
ExcInternalError());
for (unsigned int j=0; j<fe.dofs_per_line; ++j)
hierarchic_to_lexicographic_numbering[next_index++] = n*(i+1)+j+1;
+ Assert (next_index == fe.dofs_per_cell, ExcInternalError());
+
break;
- };
+ }
case 3:
{
unsigned int next_index = 0;
// first the eight vertices
- hierarchic_to_lexicographic_numbering[next_index++] = 0;
- hierarchic_to_lexicographic_numbering[next_index++] = n-1;
- hierarchic_to_lexicographic_numbering[next_index++] = (n-1)*(n*n+1);
- hierarchic_to_lexicographic_numbering[next_index++] = (n-1)*n*n;
- hierarchic_to_lexicographic_numbering[next_index++] = n*(n-1);
- hierarchic_to_lexicographic_numbering[next_index++] = n*n-1;
- hierarchic_to_lexicographic_numbering[next_index++] = n*n*n-1;
- hierarchic_to_lexicographic_numbering[next_index++] = (n-1)*(n*n+n);
+ hierarchic_to_lexicographic_numbering[next_index++] = 0; // 0
+ hierarchic_to_lexicographic_numbering[next_index++] = ( 1)*degree; // 1
+ hierarchic_to_lexicographic_numbering[next_index++] = ( n )*degree; // 2
+ hierarchic_to_lexicographic_numbering[next_index++] = ( n+1)*degree; // 3
+ hierarchic_to_lexicographic_numbering[next_index++] = (n*n )*degree; // 4
+ hierarchic_to_lexicographic_numbering[next_index++] = (n*n +1)*degree; // 5
+ hierarchic_to_lexicographic_numbering[next_index++] = (n*n+n )*degree; // 6
+ hierarchic_to_lexicographic_numbering[next_index++] = (n*n+n+1)*degree; // 7
+ // line 0
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- hierarchic_to_lexicographic_numbering[next_index++] = 1+i;
+ hierarchic_to_lexicographic_numbering[next_index++] = (i+1)*n;
+ // line 1
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- hierarchic_to_lexicographic_numbering[next_index++] = n-1+(i+1)*n*n;
+ hierarchic_to_lexicographic_numbering[next_index++] = n-1+(i+1)*n;
+ // line 2
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- hierarchic_to_lexicographic_numbering[next_index++] = n*n*(n-1)+i+1;
+ hierarchic_to_lexicographic_numbering[next_index++] = 1+i;
+ // line 3
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- hierarchic_to_lexicographic_numbering[next_index++] = (i+1)*n*n;
+ hierarchic_to_lexicographic_numbering[next_index++] = 1+i+n*(n-1);
+ // line 4
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- hierarchic_to_lexicographic_numbering[next_index++] = 1+i+n*(n-1);
+ hierarchic_to_lexicographic_numbering[next_index++] = (n-1)*n*n+(i+1)*n;
+ // line 5
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- hierarchic_to_lexicographic_numbering[next_index++] = n-1+(i+1)*n*n+n*(n-1);
+ hierarchic_to_lexicographic_numbering[next_index++] = (n-1)*(n*n+1)+(i+1)*n;
+ // line 6
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- hierarchic_to_lexicographic_numbering[next_index++] = n*n*(n-1)+i+1+n*(n-1);
+ hierarchic_to_lexicographic_numbering[next_index++] = n*n*(n-1)+i+1;
+ // line 7
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- hierarchic_to_lexicographic_numbering[next_index++] = (i+1)*n*n+n*(n-1);
+ hierarchic_to_lexicographic_numbering[next_index++] = n*n*(n-1)+i+1+n*(n-1);
+ // line 8
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- hierarchic_to_lexicographic_numbering[next_index++] = (i+1)*n;
+ hierarchic_to_lexicographic_numbering[next_index++] = (i+1)*n*n;
+ // line 9
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- hierarchic_to_lexicographic_numbering[next_index++] = n-1+(i+1)*n;
+ hierarchic_to_lexicographic_numbering[next_index++] = n-1+(i+1)*n*n;
+ // line 10
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- hierarchic_to_lexicographic_numbering[next_index++] = (n-1)*(n*n+1)+(i+1)*n;
+ hierarchic_to_lexicographic_numbering[next_index++] = (i+1)*n*n+n*(n-1);
+ // line 11
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- hierarchic_to_lexicographic_numbering[next_index++] = (n-1)*n*n+(i+1)*n;
+ hierarchic_to_lexicographic_numbering[next_index++] = n-1+(i+1)*n*n+n*(n-1);
+
// inside quads
Assert (fe.dofs_per_quad == fe.dofs_per_line*fe.dofs_per_line,
ExcInternalError());
- // quad 1
+ // face 0
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
for (unsigned int j=0; j<fe.dofs_per_line; ++j)
- hierarchic_to_lexicographic_numbering[next_index++] = (i+1)*n*n+j+1;
- // quad 2
+ hierarchic_to_lexicographic_numbering[next_index++] = (i+1)*n*n+n*(j+1);
+ // face 1
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
for (unsigned int j=0; j<fe.dofs_per_line; ++j)
- hierarchic_to_lexicographic_numbering[next_index++] = (i+1)*n*n+n*(n-1)+j+1;
- // quad 3
+ hierarchic_to_lexicographic_numbering[next_index++] = (i+1)*n*n+n-1+n*(j+1);
+ // face 2, note the orientation!
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
for (unsigned int j=0; j<fe.dofs_per_line; ++j)
- hierarchic_to_lexicographic_numbering[next_index++] = n*(i+1)+j+1;
- // quad 4
+ hierarchic_to_lexicographic_numbering[next_index++] = (j+1)*n*n+i+1;
+ // face 3, note the orientation!
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
for (unsigned int j=0; j<fe.dofs_per_line; ++j)
- hierarchic_to_lexicographic_numbering[next_index++] = (i+1)*n*n+n-1+n*(j+1);
- // quad 5
+ hierarchic_to_lexicographic_numbering[next_index++] = (j+1)*n*n+n*(n-1)+i+1;
+ // face 4
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
for (unsigned int j=0; j<fe.dofs_per_line; ++j)
- hierarchic_to_lexicographic_numbering[next_index++] = (n-1)*n*n+n*(i+1)+j+1;
- // quad 6
+ hierarchic_to_lexicographic_numbering[next_index++] = n*(i+1)+j+1;
+ // face 5
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
for (unsigned int j=0; j<fe.dofs_per_line; ++j)
- hierarchic_to_lexicographic_numbering[next_index++] = (i+1)*n*n+n*(j+1);
+ hierarchic_to_lexicographic_numbering[next_index++] = (n-1)*n*n+n*(i+1)+j+1;
// inside hex
Assert (fe.dofs_per_hex == fe.dofs_per_quad*fe.dofs_per_line,
for (unsigned int i=0; i<fe.dofs_per_line; ++i)
for (unsigned int j=0; j<fe.dofs_per_line; ++j)
for (unsigned int k=0; k<fe.dofs_per_line; ++k)
- hierarchic_to_lexicographic_numbering[next_index++]
- = n*n*(i+1)+n*(j+1)+k+1;
+ hierarchic_to_lexicographic_numbering[next_index++] = n*n*(i+1)+n*(j+1)+k+1;
+
+ Assert (next_index == fe.dofs_per_cell, ExcInternalError());
break;
- };
-
+ }
default:
Assert (false, ExcNotImplemented());
- };
+ }
// now check with data from the
// lib: there we have the mapping