--- /dev/null
+/*---------------------------- fe_lib.criss_cross.h ---------------------------*/
+/* $Id$ */
+#ifndef __fe_lib_criss_cross_H
+#define __fe_lib_criss_cross_H
+/*---------------------------- fe_lib.criss_cross.h ---------------------------*/
+
+
+#include <fe/fe.h>
+#include <fe/quadrature.h>
+
+
+/**
+ * This class implements a rather unusual macro element, the so-called
+ * criss-cross element. Its purpose is mostly to demonstrate the absence
+ * of superconvergence effects on triangular meshes where at each vertex
+ * more or less than six elements meet, but never exactly six.
+ *
+ * The construction of the element is best explained in 2d. Consider a
+ * quadrilateral with basis functions at each vertex and one at the
+ * crossing-point of the two diagonals. The element is divided by the
+ * diagonals into four triangles and assume that each vertex basis
+ * function has support only on the two triangles adjacent to the
+ * respective vertex and is constant zero on the other two triangles;
+ * they are linear on each of the triangles and globally continuous.
+ * The center basis function lives on each of the four triangles, is
+ * linear on each triangles and vanishes at the faces of the quadrilateral.
+ *
+ * Now, on the unit element, these basis functions are the same as for
+ * a triangular ansatz space, namely the class of ${\cal P}_1$ Lagrange
+ * elements. Due to the arrangement of the four triangles on the
+ * quadrilateral, it is clear that considering the whole triangulation
+ * of the domain, always four triangles meet at the points which
+ * correspond with the centers of the quadrilaterals and $2*s$ triangles
+ * meet at the vertices of the quadrilaterals, if $s$ is the number of
+ * quadrilaterals meeting there. Thus, in most cases the number of
+ * triangles meeting are four or eight, which effectively destroys
+ * superconvergence at nodes.
+ *
+ * This element is not quite equivalent from beginning to the linear
+ * triangular elements. The reason for this is that if we use a bilinear
+ * mapping of the unit quadrilateral to the real cell, the diagonals will
+ * in general not be straight lines. Therefore, the shape functions will
+ * in general not be linear on the real cell, unlike for the linear
+ * triangular element, which uses a linear mapping. The missing linearity
+ * makes assemblage of matrices a bit more complicated, since the gradient
+ * is not constant and we need more than one quadrature point, as well
+ * as some other dubtle difficulties. This problem can, however, be cured
+ * using a trick: the usual transformation from unit coordinates $\vec\xi$
+ * to real coordinates $\vec x(\vec\xi)$ looks like
+ * $$
+ * \vec x(\vec\xi) = \sum_{i=0}^3 \phi_i^L(\vec\xi) \vec x_i
+ * $$
+ * with $\phi_i^L$ being the bilinear basis functions associated with the
+ * vertices and $\vec x_i$ being the coordinates of the vertices in real
+ * space. Now, we could also choose
+ * $$
+ * \vec x(\vec\xi) = \sum_{i=0}^4 \phi_i(\vec\xi) \vec x_i
+ * $$
+ * with the basis functions $\phi_i$ of this element, the four vertices
+ * in real space $\vec x_0..\vec x_3$ and an interior point in real space
+ * $\vec x_4$. We can choose the interior point quite arbitrarily and it
+ * will become clear in a moment how we have to do so. First let us note
+ * that because the vertex basis functions are linear on the faces,
+ * because they vanish on the two faces not adjacent to the associated
+ * vertex and because the center basis function vanishes at the four
+ * faces, the four sides of the unit cell are mapped to straight lines
+ * in real space, just like for the bilinear mapping.
+ *
+ * Now, to ensure that the mapping of each of the four triangles to the
+ * real space is linear, we have to require that the two diagonals are
+ * mapped to straight lines. One necessary condition for this is, that the
+ * center point of the unit cell is mapped to the crossing point of the
+ * two diagonals in real space. Therefore, we choose $\vec x_4$ to be
+ * this point. Then we note, that because the vertex basis functions vanish
+ * on the diagonal not through the vertex and are constant zero beyond that,
+ * the mapping of the line from the center to a vertex is described entirely
+ * by the basis function of that vertex and the center basis function; but
+ * because they both are linear on that line, the line is also a straight
+ * one in real space. This proves that by this construction of the mapping
+ * between unit and real cell, the mapping of each of the four triangles
+ * is linear (note that this does not hold for the whole element; the
+ * mapping of the quadrilaterals is only piecewise linear and globally
+ * continuous). It also proves that the ansatz space using this element
+ * is equivalent to the ansatz space using triangles and linear elements.
+ *
+ * Since in one space dimension, this element equals two linear elements,
+ * i.e. a linear ansatz on a mesh once more refined than the present one,
+ * this element is not implemented for one dimension. There may be an
+ * analogue to the criss-cross element in more than two space dimensions,
+ * but it is not implemented at present.
+ *
+ * As stated above, the element is not really a good one. It may, however,
+ * serve to study superconvergence effects and also to satisfy the author's
+ * curiosity. At least for the first of these two reasons, it is better
+ * suited than using a genuine triangulation of the domain (i.e. using real
+ * triangles rather than subdividing quadrilaterals as shown above), since
+ * the construction of triangulations with four or eight cells meeting at
+ * each vertex is certainly not feasible other than by hand, while the
+ * decomposition of a domain using quadrilaterals is easier.
+ *
+ *
+ * \section{Hanging nodes}
+ *
+ * Hanging nodes are handled exactly like for any other element. It should
+ * however be noted that the support of basis functions get quite
+ * complicated in the presence of hanging nodes, as the following figure
+ * depicts:
+ * \begin{verbatim}
+ * *-----------------*--------*----
+ * | /|\ |
+ * | /..|.\ |
+ * | /....|...\ |
+ * | /......|.....\ |
+ * | /.......|.......\|
+ * | /.........*--------*----
+ * | /..........|......./|
+ * | /............|....../ |
+ * | /.............|..../ |
+ * | /...............|.....\ |
+ * |/................|.......\|
+ * *-----------------o--------*-----
+ * \end{verbatim}
+ * The dotted area is the support of the basis function associated with the
+ * bottom middle vertex (denoted by #o#) after the hanging node in the center
+ * of the `picture' was eliminated. This strange structure of the support
+ * should not pose too many problems in practice, it is only note here for
+ * completeness and for curiosity.
+ *
+ *
+ * \section{Experience with the criss-cross element}
+ *
+ * Experience is that the error for the criss cross element shows
+ * the same convergence rate as the linear Lagrange element ($h^2$ for the
+ * $L^2$ error, $h$ for the $H^1$ error). The $L^2$ error is about the same
+ * size for the same number of elements as for the linear element; since
+ * the criss-cross elements has about twice as many degrees of freedom as
+ * the linear element for the same triangulation, the $L^2$ error really
+ * is about twice as large as a function of the number of degrees of freedom.
+ *
+ * Converse to that, the $H^1$ error is about a factor of 1.2 smaller for
+ * the same number of degrees of freedoms.
+ *
+ * Apart from all this data, it must not be forgotten that we cannot
+ * expect superconvergence neither in the Gauss points not in the vertices.
+ * Thus we may improve the accuracy of the solution obtained with the linear
+ * element by a postprocess, while we can't do so for the criss-cross element.
+ *
+ * All given data refer to a Poisson equation with nonhomogeneous boundary
+ * values on the unit disk (resp. a triangulation of that) and hanging nodes.
+ *
+ *
+ * \section{Using quadrature formulae for this element}
+ *
+ * When using one of the usual quadrature formulae, a common problem is
+ * that some of the quadrature points lie on the interfaces of the
+ * triangles. For this reason, there is a family of quadrature formulae
+ * defined below, names \ref{QCrissCross1} and higher order, which
+ * resemble the quadrature formulae used on triangular domains, but
+ * taken four-fold, i.e. for each of the four subtriangles.
+ *
+ *
+ * @author Wolfgang Bangerth, 1998
+ */
+template <int dim>
+class FECrissCross : public FiniteElement<dim> {
+ public:
+ /**
+ * Constructor
+ */
+ FECrissCross ();
+
+ /**
+ * Return the value of the #i#th shape
+ * function at point #p# on the unit cell.
+ */
+ virtual double shape_value(const unsigned int i,
+ const Point<dim>& p) const;
+
+ /**
+ * Return the gradient of the #i#th shape
+ * function at point #p# on the unit cell.
+ */
+ virtual Point<dim> shape_grad(const unsigned int i,
+ const Point<dim>& p) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ */
+ virtual void get_unit_ansatz_points (vector<Point<dim> > &ansatz_points) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ */
+ virtual void get_ansatz_points (const DoFHandler<dim>::cell_iterator &cell,
+ const Boundary<dim> &boundary,
+ vector<Point<dim> > &ansatz_points) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ */
+ virtual void get_face_ansatz_points (const DoFHandler<dim>::face_iterator &face,
+ const Boundary<dim> &boundary,
+ vector<Point<dim> > &ansatz_points) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ */
+ virtual void get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &cell,
+ const Boundary<dim> &boundary,
+ dFMatrix &local_mass_matrix) const;
+
+ /**
+ * Return the value of the #i#th shape
+ * function at point #p# on the unit cell.
+ * Here, the (bi-)linear basis functions
+ * are meant, which are used for the
+ * computation of the transformation from
+ * unit cell to real space cell.
+ */
+ virtual double shape_value_transform (const unsigned int i,
+ const Point<dim> &p) const;
+
+ /**
+ * Return the gradient of the #i#th shape
+ * function at point #p# on the unit cell.
+ * Here, the (bi-)linear basis functions
+ * are meant, which are used for the
+ * computation of the transformation from
+ * unit cell to real space cell.
+ */
+ virtual Point<dim> shape_grad_transform (const unsigned int i,
+ const Point<dim> &p) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ *
+ * In two spatial dimensions, this function
+ * simply returns the length of the face.
+ */
+ virtual void get_face_jacobians (const DoFHandler<dim>::face_iterator &face,
+ const Boundary<dim> &boundary,
+ const vector<Point<dim-1> > &unit_points,
+ vector<double> &face_jacobi_determinants) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ *
+ * In two spatial dimensions, this function
+ * simply returns half the length of the
+ * whole face.
+ */
+ virtual void get_subface_jacobians (const DoFHandler<dim>::face_iterator &face,
+ const unsigned int subface_no,
+ const vector<Point<dim-1> > &unit_points,
+ vector<double> &face_jacobi_determinants) const;
+
+ /**
+ * Return the normal vectors to the
+ * face with number #face_no# of #cell#.
+ *
+ * For linear finite elements, this function
+ * is particularly simple since all normal
+ * vectors are equal and can easiliy be
+ * computed from the direction of the face
+ * without using the transformation (Jacobi)
+ * matrix, at least for two dimensions.
+ *
+ * Refer to the base class for detailed
+ * information on this function.
+ */
+ virtual void get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const Boundary<dim> &boundary,
+ const vector<Point<dim-1> > &unit_points,
+ vector<Point<dim> > &normal_vectors) const;
+
+ /**
+ * Return the normal vectors to the
+ * subface with number #subface_no# of
+ * the face with number #face_no# of #cell#.
+ *
+ * For linear finite elements, this function
+ * is particularly simple since all normal
+ * vectors are equal and can easiliy be
+ * computed from the direction of the face
+ * without using the transformation (Jacobi)
+ * matrix, at least for two dimensions.
+ *
+ * Refer to the base class for detailed
+ * information on this function.
+ */
+ virtual void get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const unsigned int subface_no,
+ const vector<Point<dim-1> > &unit_points,
+ vector<Point<dim> > &normal_vectors) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ *
+ * For one dimensional elements, this
+ * function simply passes through to
+ * the one implemented in the base class.
+ * For higher dimensional finite elements
+ * we use linear mappings and therefore
+ * the boundary object is ignored since
+ * the boundary is approximated using
+ * piecewise multilinear boundary segments.
+ */
+ virtual void fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
+ const vector<Point<dim> > &unit_points,
+ vector<dFMatrix> &jacobians,
+ const bool compute_jacobians,
+ vector<Point<dim> > &ansatz_points,
+ const bool compute_ansatz_points,
+ vector<Point<dim> > &q_points,
+ const bool compute_q_points,
+ const dFMatrix &shape_values_transform,
+ const vector<vector<Point<dim> > > &shape_grad_transform,
+ const Boundary<dim> &boundary) const;
+
+ DeclException0 (ExcNotUseful);
+};
+
+
+
+
+/**
+ * Quadrature formula for the criss-cross element. This quadrature
+ * formula uses one point at the barycenter of each of the four subtriangles.
+ *
+ * For the same reason as for the criss-cross element itself, this
+ * formula is not implemented for one space dimension.
+ */
+template <int dim>
+class QCrissCross1 : public Quadrature<dim> {
+ public:
+ QCrissCross1 ();
+
+ DeclException0 (ExcNotUseful);
+};
+
+
+
+/*---------------------------- fe_lib.criss_cross.h ---------------------------*/
+/* end of #ifndef __fe_lib_criss_cross_H */
+#endif
+/*---------------------------- fe_lib.criss_cross.h ---------------------------*/
--- /dev/null
+/* $Id$ */
+
+#include <fe/fe_lib.criss_cross.h>
+#include <grid/tria_iterator.h>
+#include <grid/dof_accessor.h>
+#include <grid/geometry_info.h>
+
+
+
+
+/*-----------------------------------2d------------------------------------
+ Maple script to automate some of the error-prone computations on
+ this strange sort of element.
+
+ n_functions := 5:
+
+ # note: ansatz_points[i] is a vector which is indexed from
+ # one and not from zero!
+ ansatz_points[0] := [0,0]:
+ ansatz_points[1] := [1,0]:
+ ansatz_points[2] := [1,1]:
+ ansatz_points[3] := [0,1]:
+ ansatz_points[4] := [1/2,1/2]:
+
+ phi[0] := proc(x,y) if(y<1-x) then 1-x-y; else 0; fi; end:
+ phi[1] := proc(x,y) if(y<x) then x-y; else 0; fi; end:
+ phi[2] := proc(x,y) if(y>1-x) then x+y-1; else 0; fi; end:
+ phi[3] := proc(x,y) if(y>x) then y-x; else 0; fi; end:
+ phi[4] := proc(x,y) 1 - phi[0](x,y) - phi[1](x,y)
+ - phi[2](x,y) - phi[3](x,y) ; end:
+
+ #points on children: let them be indexed one-based, as are
+ #the ansatz_points
+ points[0] := array(0..n_functions-1, 1..2):
+ points[1] := array(0..n_functions-1, 1..2):
+ points[2] := array(0..n_functions-1, 1..2):
+ points[3] := array(0..n_functions-1, 1..2):
+ for i from 0 to n_functions-1 do
+ points[0][i,1] := ansatz_points[i][1]/2:
+ points[0][i,2] := ansatz_points[i][2]/2:
+
+ points[1][i,1] := ansatz_points[i][1]/2+1/2:
+ points[1][i,2] := ansatz_points[i][2]/2:
+
+ points[2][i,1] := ansatz_points[i][1]/2+1/2:
+ points[2][i,2] := ansatz_points[i][2]/2+1/2:
+
+ points[3][i,1] := ansatz_points[i][1]/2:
+ points[3][i,2] := ansatz_points[i][2]/2+1/2:
+ od:
+
+ prolongation := array(0..3,0..n_functions-1, 0..n_functions-1):
+ print ("Computing prolongation matrices"):
+ for i from 0 to 3 do
+ print ("child", i):
+ for j from 0 to n_functions-1 do
+ for k from 0 to n_functions-1 do
+ prolongation[i,j,k] := phi[k](points[i][j,1], points[i][j,2]):
+ od:
+ od:
+ od:
+
+ eq_sys := {(1-t)*x0 + t*x2 = (1-s)*x1 + s*x3,
+ (1-t)*y0 + t*y2 = (1-s)*y1 + s*y3}:
+ solution := solve (eq_sys, {s,t});
+
+ xs := subs (solution, (1-t)*x0 + t*x2):
+ ys := subs (solution, (1-t)*y0 + t*y2):
+ ps := array(1..2, [xs, ys]):
+
+ print ("writing data to files"):
+ readlib(C):
+ C(prolongation, filename=prolongation_2d):
+ C(ps, filename=crosspoint):
+
+ --------------------------------------------------------------------
+
+ Postprocess the prolongation matrix by the commands
+
+ perl -pi -e 's/\[(\d+)\]\[(\d+)\]\[(\d+)\]/[$1]($2,$3)/g;' prolongation_2d
+ perl -pi -e 's/.*= 0.0;\n//g;' prolongation_2d
+
+-----------------------------------------------------------------------------*/
+
+
+#if deal_II_dimension == 1
+
+template <>
+FECrissCross<1>::FECrissCross () :
+ FiniteElement<1> (0,0,0,0)
+{
+ Assert (false, ExcNotUseful());
+};
+
+
+
+template <>
+double FECrissCross<1>::shape_value (const unsigned int, const Point<1> &) const {
+ Assert (false, ExcNotUseful());
+ return 0;
+};
+
+
+
+template <>
+Point<1> FECrissCross<1>::shape_grad (const unsigned int, const Point<1> &) const {
+ Assert (false, ExcNotUseful());
+ return Point<1>();
+};
+
+
+
+template <>
+void FECrissCross<1>::get_unit_ansatz_points (vector<Point<1> >&) const {
+ Assert (false, ExcNotUseful());
+};
+
+
+
+template <>
+void FECrissCross<1>::get_ansatz_points (const DoFHandler<1>::cell_iterator &,
+ const Boundary<1> &,
+ vector<Point<1> > &) const {
+ Assert (false, ExcNotUseful());
+};
+
+
+
+template <>
+void FECrissCross<1>::get_face_ansatz_points (const DoFHandler<1>::face_iterator &,
+ const Boundary<1> &,
+ vector<Point<1> > &) const {
+ Assert (false, ExcNotUseful());
+};
+
+
+
+template <>
+void FECrissCross<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &,
+ const Boundary<1> &,
+ dFMatrix &) const {
+ Assert (false, ExcNotUseful());
+};
+
+
+
+template <>
+double FECrissCross<1>::shape_value_transform (const unsigned int,
+ const Point<1> &) const {
+ Assert (false, ExcNotUseful());
+ return 0;
+};
+
+
+
+template <>
+Point<1> FECrissCross<1>::shape_grad_transform (const unsigned int,
+ const Point<1> &) const {
+ Assert (false, ExcNotUseful());
+ return Point<1>();
+};
+
+
+
+template <>
+void FECrissCross<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
+ const Boundary<1> &,
+ const vector<Point<0> > &,
+ vector<double> &) const {
+ Assert (false, ExcNotUseful());
+};
+
+
+
+template <>
+void FECrissCross<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
+ const unsigned int,
+ const vector<Point<0> > &,
+ vector<double> &) const {
+ Assert (false, ExcNotUseful());
+};
+
+
+
+template <>
+void FECrissCross<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
+ const unsigned int,
+ const Boundary<1> &,
+ const vector<Point<0> > &,
+ vector<Point<1> > &) const {
+ Assert (false, ExcNotUseful());
+};
+
+
+
+template <>
+void FECrissCross<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
+ const unsigned int,
+ const unsigned int,
+ const vector<Point<0> > &,
+ vector<Point<1> > &) const {
+ Assert (false, ExcNotUseful());
+};
+
+
+
+template <>
+void FECrissCross<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &,
+ const vector<Point<1> > &,
+ vector<dFMatrix> &,
+ const bool,
+ vector<Point<1> > &,
+ const bool,
+ vector<Point<1> > &,
+ const bool,
+ const dFMatrix &,
+ const vector<vector<Point<1> > > &,
+ const Boundary<1> &) const {
+ Assert (false, ExcNotUseful());
+};
+
+#endif
+
+
+
+
+
+
+#if deal_II_dimension == 2
+
+template <>
+FECrissCross<2>::FECrissCross () :
+ FiniteElement<2> (1,0,1,5)
+{
+ interface_constraints(0,0) = 1./2.;
+ interface_constraints(0,1) = 1./2.;
+
+ prolongation[0](0,0) = 1.0;
+ prolongation[0](1,0) = 1.0/2.0;
+ prolongation[0](1,1) = 1.0/2.0;
+ prolongation[0](2,4) = 1.0;
+ prolongation[0](3,0) = 1.0/2.0;
+ prolongation[0](3,3) = 1.0/2.0;
+ prolongation[0](4,0) = 1.0/2.0;
+ prolongation[0](4,4) = 1.0/2.0;
+ prolongation[1](0,0) = 1.0/2.0;
+ prolongation[1](0,1) = 1.0/2.0;
+ prolongation[1](1,1) = 1.0;
+ prolongation[1](2,1) = 1.0/2.0;
+ prolongation[1](2,2) = 1.0/2.0;
+ prolongation[1](3,4) = 1.0;
+ prolongation[1](4,1) = 1.0/2.0;
+ prolongation[1](4,4) = 1.0/2.0;
+ prolongation[2](0,4) = 1.0;
+ prolongation[2](1,1) = 1.0/2.0;
+ prolongation[2](1,2) = 1.0/2.0;
+ prolongation[2](2,2) = 1.0;
+ prolongation[2](3,2) = 1.0/2.0;
+ prolongation[2](3,3) = 1.0/2.0;
+ prolongation[2](4,2) = 1.0/2.0;
+ prolongation[2](4,4) = 1.0/2.0;
+ prolongation[3](0,0) = 1.0/2.0;
+ prolongation[3](0,3) = 1.0/2.0;
+ prolongation[3](1,4) = 1.0;
+ prolongation[3](2,2) = 1.0/2.0;
+ prolongation[3](2,3) = 1.0/2.0;
+ prolongation[3](3,3) = 1.0;
+ prolongation[3](4,3) = 1.0/2.0;
+ prolongation[3](4,4) = 1.0/2.0;
+};
+
+
+
+template <>
+inline
+double FECrissCross<2>::shape_value (const unsigned int i,
+ const Point<2> &p) const {
+ Assert((i<total_dofs), ExcInvalidIndex(i));
+
+ const double x = p(0),
+ y = p(1);
+ switch (i)
+ {
+ case 0: return ((y<1-x) ? 1-x-y : 0);
+ case 1: return ((y<x) ? x-y : 0);
+ case 2: return ((y>1-x) ? x+y-1 : 0);
+ case 3: return ((y>x) ? y-x : 0);
+
+ // I am too lazy to optimize the
+ // following myself. Let the
+ // compiler do this
+ case 4: return (1-(((y<1-x) ? 1-x-y : 0) +
+ ((y<x) ? x-y : 0) +
+ ((y>1-x) ? x+y-1 : 0) +
+ ((y>x) ? y-x : 0)));
+ }
+ return 0.;
+};
+
+
+
+template <>
+inline
+Point<2> FECrissCross<2>::shape_grad (const unsigned int i, const Point<2> &p) const {
+ Assert((i<total_dofs), ExcInvalidIndex(i));
+
+ const double x = p(0),
+ y = p(1);
+ switch (i)
+ {
+ case 0: return ((y<1-x) ? Point<2>(-1,-1) : Point<2>(0,0));
+ case 1: return ((y<x) ? Point<2>(1,-1) : Point<2>(0,0));
+ case 2: return ((y>1-x) ? Point<2>(1,1) : Point<2>(0,0));
+ case 3: return ((y>x) ? Point<2>(-1,1) : Point<2>(0,0));
+
+ // I am too lazy to optimize the
+ // following myself. Let the
+ // compiler do this
+ case 4: return -1.*(((y<1-x) ? Point<2>(-1,-1) : Point<2>(0,0)) +
+ ((y<x) ? Point<2>(1,-1) : Point<2>(0,0)) +
+ ((y>1-x) ? Point<2>(1,1) : Point<2>(0,0)) +
+ ((y>x) ? Point<2>(-1,1) : Point<2>(0,0)));
+ }
+ return Point<2>();
+};
+
+
+
+template <>
+void FECrissCross<2>::get_unit_ansatz_points (vector<Point<2> > &unit_points) const {
+ Assert(unit_points.size()==total_dofs,
+ ExcWrongFieldDimension (unit_points.size(), total_dofs));
+
+ unit_points[0] = Point<2> (0,0);
+ unit_points[1] = Point<2> (1,0);
+ unit_points[2] = Point<2> (1,1);
+ unit_points[3] = Point<2> (0,1);
+ unit_points[4] = Point<2> (0.5, 0.5);
+};
+
+
+
+template <>
+void FECrissCross<2>::get_ansatz_points (const DoFHandler<2>::cell_iterator &cell,
+ const Boundary<2> &,
+ vector<Point<2> > &ansatz_points) const {
+ const unsigned int dim = 2;
+
+ Assert (ansatz_points.size() == total_dofs,
+ ExcWrongFieldDimension (ansatz_points.size(), total_dofs));
+
+ // copy vertices
+ for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell; ++vertex)
+ ansatz_points[vertex] = cell->vertex(vertex);
+
+/*
+ last ansatz point is the common point of the two diagonals; the formula for
+ the computation is a bit lengthy but straightforward. You can get it with
+ the small Maple script printed at the beginning of this file.
+*/
+ const double x0 = cell->vertex(0)(0),
+ y0 = cell->vertex(0)(1),
+ x1 = cell->vertex(1)(0),
+ y1 = cell->vertex(1)(1),
+ x2 = cell->vertex(2)(0),
+ y2 = cell->vertex(2)(1),
+ x3 = cell->vertex(3)(0),
+ y3 = cell->vertex(3)(1);
+ const double t1 = x0*y1;
+ const double t2 = x0*y3;
+ const double t4 = x1*y0;
+ const double t5 = x3*y0;
+ const double t14 = (t1-t2+x1*y3-t4+t5-x3*y1)/(t1-t2-x2*y1+x2*y3-t4+x1*y2+t5-x3*y2);
+ const double t15 = 1.0-t14;
+ ansatz_points[4](0) = t15*x0+t14*x2;
+ ansatz_points[4](1) = t15*y0+t14*y2;
+};
+
+
+
+template <>
+void FECrissCross<2>::get_face_ansatz_points (const DoFHandler<2>::face_iterator &face,
+ const Boundary<2> &,
+ vector<Point<2> > &ansatz_points) const {
+ const unsigned int dim = 2;
+
+ Assert ((ansatz_points.size() == dofs_per_face) &&
+ (ansatz_points.size() == GeometryInfo<dim>::vertices_per_face),
+ ExcWrongFieldDimension (ansatz_points.size(),
+ GeometryInfo<dim>::vertices_per_face));
+
+ for (unsigned int vertex=0; vertex<dofs_per_face; ++vertex)
+ ansatz_points[vertex] = face->vertex(vertex);
+};
+
+
+
+template <>
+void FECrissCross<2>::get_local_mass_matrix (const DoFHandler<2>::cell_iterator &,
+ const Boundary<2> &,
+ dFMatrix &local_mass_matrix) const {
+ Assert (local_mass_matrix.n() == total_dofs,
+ ExcWrongFieldDimension(local_mass_matrix.n(),total_dofs));
+ Assert (local_mass_matrix.m() == total_dofs,
+ ExcWrongFieldDimension(local_mass_matrix.m(),total_dofs));
+
+ // in this special element, some of the
+ // entries are zero (which is not the
+ // case for most other elements, so
+ // we first reset all elements and only
+ // fill in those that are nonzero
+ local_mass_matrix.clear ();
+
+ Assert (false, ExcNotUseful());
+};
+
+
+
+template <>
+inline
+double FECrissCross<2>::shape_value_transform (const unsigned int i,
+ const Point<2> &p) const {
+ // use an isoparametric ansatz
+ return shape_value(i,p);
+};
+
+
+
+template <>
+Point<2> FECrissCross<2>::shape_grad_transform (const unsigned int i,
+ const Point<2> &p) const {
+ // use an isoparametric ansatz
+ return shape_grad(i,p);
+};
+
+
+
+template <>
+void FECrissCross<2>::get_face_jacobians (const DoFHandler<2>::face_iterator &face,
+ const Boundary<2> &,
+ const vector<Point<1> > &unit_points,
+ vector<double> &face_jacobians) const {
+ // more or less copied from the linear
+ // finite element
+ Assert (unit_points.size() == face_jacobians.size(),
+ ExcWrongFieldDimension (unit_points.size(), face_jacobians.size()));
+
+ // a linear mapping for a single line
+ // produces particularly simple
+ // expressions for the jacobi
+ // determinant :-)
+ const double h = sqrt((face->vertex(1) - face->vertex(0)).square());
+ fill_n (face_jacobians.begin(),
+ unit_points.size(),
+ h);
+};
+
+
+
+template <>
+void FECrissCross<2>::get_subface_jacobians (const DoFHandler<2>::face_iterator &face,
+ const unsigned int,
+ const vector<Point<1> > &unit_points,
+ vector<double> &face_jacobians) const {
+ // more or less copied from the linear
+ // finite element
+ Assert (unit_points.size() == face_jacobians.size(),
+ ExcWrongFieldDimension (unit_points.size(), face_jacobians.size()));
+ Assert (face->at_boundary() == false,
+ ExcBoundaryFaceUsed ());
+
+ // a linear mapping for a single line
+ // produces particularly simple
+ // expressions for the jacobi
+ // determinant :-)
+ const double h = sqrt((face->vertex(1) - face->vertex(0)).square());
+ fill_n (face_jacobians.begin(),
+ unit_points.size(),
+ h/2);
+};
+
+
+
+template <>
+void FECrissCross<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell,
+ const unsigned int face_no,
+ const Boundary<2> &,
+ const vector<Point<1> > &unit_points,
+ vector<Point<2> > &normal_vectors) const {
+ // more or less copied from the linear
+ // finite element
+ Assert (unit_points.size() == normal_vectors.size(),
+ ExcWrongFieldDimension (unit_points.size(), normal_vectors.size()));
+
+ const DoFHandler<2>::face_iterator face = cell->face(face_no);
+ // compute direction of line
+ const Point<2> line_direction = (face->vertex(1) - face->vertex(0));
+ // rotate to the right by 90 degrees
+ const Point<2> normal_direction(line_direction(1),
+ -line_direction(0));
+
+ if (face_no <= 1)
+ // for sides 0 and 1: return the correctly
+ // scaled vector
+ fill (normal_vectors.begin(), normal_vectors.end(),
+ normal_direction / sqrt(normal_direction.square()));
+ else
+ // for sides 2 and 3: scale and invert
+ // vector
+ fill (normal_vectors.begin(), normal_vectors.end(),
+ normal_direction / (-sqrt(normal_direction.square())));
+};
+
+
+
+template <>
+void FECrissCross<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell,
+ const unsigned int face_no,
+ const unsigned int,
+ const vector<Point<1> > &unit_points,
+ vector<Point<2> > &normal_vectors) const {
+ // more or less copied from the linear
+ // finite element
+ // note, that in 2D the normal vectors to the
+ // subface have the same direction as that
+ // for the face
+ Assert (unit_points.size() == normal_vectors.size(),
+ ExcWrongFieldDimension (unit_points.size(), normal_vectors.size()));
+ Assert (cell->face(face_no)->at_boundary() == false,
+ ExcBoundaryFaceUsed ());
+
+ const DoFHandler<2>::face_iterator face = cell->face(face_no);
+ // compute direction of line
+ const Point<2> line_direction = (face->vertex(1) - face->vertex(0));
+ // rotate to the right by 90 degrees
+ const Point<2> normal_direction(line_direction(1),
+ -line_direction(0));
+
+ if (face_no <= 1)
+ // for sides 0 and 1: return the correctly
+ // scaled vector
+ fill (normal_vectors.begin(), normal_vectors.end(),
+ normal_direction / sqrt(normal_direction.square()));
+ else
+ // for sides 2 and 3: scale and invert
+ // vector
+ fill (normal_vectors.begin(), normal_vectors.end(),
+ normal_direction / (-sqrt(normal_direction.square())));
+};
+
+
+
+template <int dim>
+void FECrissCross<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
+ const vector<Point<dim> > &unit_points,
+ vector<dFMatrix> &jacobians,
+ const bool compute_jacobians,
+ vector<Point<dim> > &ansatz_points,
+ const bool,
+ vector<Point<dim> > &q_points,
+ const bool compute_q_points,
+ const dFMatrix &shape_values_transform,
+ const vector<vector<Point<dim> > > &shape_grad_transform,
+ const Boundary<dim> &boundary) const {
+ Assert (jacobians.size() == unit_points.size(),
+ ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
+ Assert (q_points.size() == unit_points.size(),
+ ExcWrongFieldDimension(q_points.size(), unit_points.size()));
+ Assert (ansatz_points.size() == total_dofs,
+ ExcWrongFieldDimension(ansatz_points.size(), total_dofs));
+
+
+ unsigned int n_points=unit_points.size();
+
+ // we need the ansatz points in any
+ // way, wanted or not by the user
+ get_ansatz_points (cell, boundary, ansatz_points);
+
+ if (compute_q_points)
+ {
+ // initialize points to zero
+ for (unsigned int i=0; i<n_points; ++i)
+ q_points[i] = Point<dim> ();
+
+ // note: let x_l be the vector of the
+ // lth quadrature point in real space and
+ // xi_l that on the unit cell, let further
+ // p_j be the vector of the jth vertex
+ // of the cell in real space and
+ // N_j(xi_l) be the value of the associated
+ // basis function at xi_l, then
+ // x_l(xi_l) = sum_j p_j N_j(xi_l)
+ for (unsigned int j=0; j<n_transform_functions; ++j)
+ for (unsigned int l=0; l<n_points; ++l)
+ q_points[l] += ansatz_points[j] * shape_values_transform(j, l);
+ };
+
+
+/* jacobi matrices: compute d(x)/d(xi) and invert this
+ Let M(l) be the inverse of J at the quadrature point l, then
+ M_{ij}(l) = sum_s p_i(s) d(N_s(l))/d(xi_j)
+ where p_i(s) is the i-th coordinate of the s-th vertex vector,
+ N_s(l) is the value of the s-th shape function at the
+ quadrature point l.
+
+ We could therefore write:
+ l=0..n_points-1
+ i=0..dim-1
+ j=0..dim-1
+ M_{ij}(l) = 0
+ s=0..n_ansatz_points
+ M_{ij}(l) += p_i(s) d(N_s(l))/d(xi_j)
+
+ However, we rewrite the loops to only compute the gradient once for
+ each integration point and basis function.
+*/
+ if (compute_jacobians)
+ {
+ dFMatrix M(dim,dim);
+ for (unsigned int l=0; l<n_points; ++l)
+ {
+ M.clear ();
+ for (unsigned int s=0; s<n_transform_functions; ++s)
+ {
+ // we want the linear transform,
+ // so use that function
+ const Point<dim> gradient = shape_grad_transform[s][l];
+ for (unsigned int i=0; i<dim; ++i)
+ for (unsigned int j=0; j<dim; ++j)
+ M(i,j) += ansatz_points[s](i) * gradient(j);
+ };
+ jacobians[l].invert(M);
+ };
+ };
+};
+
+#endif
+
+
+
+
+
+/*--------------------------- QCrissCross* ------------------------------------*/
+
+
+#if deal_II_dimension == 1
+
+template <>
+QCrissCross1<1>::QCrissCross1 () :
+ Quadrature<1> (1)
+{
+ Assert (false, ExcNotUseful());
+};
+
+#endif
+
+
+
+#if deal_II_dimension == 2
+
+template <>
+QCrissCross1<2>::QCrissCross1 () :
+ Quadrature<2> (4)
+{
+ // let quadrature points be the
+ // barycenters of the four triangles
+ quadrature_points[0] = Point<2>(1./2., 1./6.);
+ quadrature_points[1] = Point<2>(5./6., 1./2.);
+ quadrature_points[2] = Point<2>(1./2., 5./6.);
+ quadrature_points[3] = Point<2>(1./6., 1./2.);
+
+ weights[0] = 1./4.;
+ weights[1] = 1./4.;
+ weights[2] = 1./4.;
+ weights[3] = 1./4.;
+};
+
+#endif
+
+
+
+// explicit instantiations
+
+template class FECrissCross<deal_II_dimension>;
+template class QCrissCross1<deal_II_dimension>;