From: Wolfgang Bangerth Date: Mon, 27 Jul 1998 22:11:06 +0000 (+0000) Subject: Add the criss-cross element. X-Git-Tag: v8.0.0~22797 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6cdde953c580b54ef5806d4dba03d2c25672f372;p=dealii.git Add the criss-cross element. git-svn-id: https://svn.dealii.org/trunk@457 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/fe_lib.criss_cross.h b/deal.II/deal.II/include/fe/fe_lib.criss_cross.h new file mode 100644 index 0000000000..dde7dccce4 --- /dev/null +++ b/deal.II/deal.II/include/fe/fe_lib.criss_cross.h @@ -0,0 +1,355 @@ +/*---------------------------- fe_lib.criss_cross.h ---------------------------*/ +/* $Id$ */ +#ifndef __fe_lib_criss_cross_H +#define __fe_lib_criss_cross_H +/*---------------------------- fe_lib.criss_cross.h ---------------------------*/ + + +#include +#include + + +/** + * 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 +class FECrissCross : public FiniteElement { + 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& p) const; + + /** + * Return the gradient of the #i#th shape + * function at point #p# on the unit cell. + */ + virtual Point shape_grad(const unsigned int i, + const Point& p) const; + + /** + * Refer to the base class for detailed + * information on this function. + */ + virtual void get_unit_ansatz_points (vector > &ansatz_points) const; + + /** + * Refer to the base class for detailed + * information on this function. + */ + virtual void get_ansatz_points (const DoFHandler::cell_iterator &cell, + const Boundary &boundary, + vector > &ansatz_points) const; + + /** + * Refer to the base class for detailed + * information on this function. + */ + virtual void get_face_ansatz_points (const DoFHandler::face_iterator &face, + const Boundary &boundary, + vector > &ansatz_points) const; + + /** + * Refer to the base class for detailed + * information on this function. + */ + virtual void get_local_mass_matrix (const DoFHandler::cell_iterator &cell, + const Boundary &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 &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 shape_grad_transform (const unsigned int i, + const Point &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::face_iterator &face, + const Boundary &boundary, + const vector > &unit_points, + vector &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::face_iterator &face, + const unsigned int subface_no, + const vector > &unit_points, + vector &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::cell_iterator &cell, + const unsigned int face_no, + const Boundary &boundary, + const vector > &unit_points, + vector > &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::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const vector > &unit_points, + vector > &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::cell_iterator &cell, + const vector > &unit_points, + vector &jacobians, + const bool compute_jacobians, + vector > &ansatz_points, + const bool compute_ansatz_points, + vector > &q_points, + const bool compute_q_points, + const dFMatrix &shape_values_transform, + const vector > > &shape_grad_transform, + const Boundary &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 +class QCrissCross1 : public Quadrature { + public: + QCrissCross1 (); + + DeclException0 (ExcNotUseful); +}; + + + +/*---------------------------- fe_lib.criss_cross.h ---------------------------*/ +/* end of #ifndef __fe_lib_criss_cross_H */ +#endif +/*---------------------------- fe_lib.criss_cross.h ---------------------------*/ diff --git a/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc b/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc new file mode 100644 index 0000000000..9604328893 --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc @@ -0,0 +1,685 @@ +/* $Id$ */ + +#include +#include +#include +#include + + + + +/*-----------------------------------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(y1-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 >&) const { + Assert (false, ExcNotUseful()); +}; + + + +template <> +void FECrissCross<1>::get_ansatz_points (const DoFHandler<1>::cell_iterator &, + const Boundary<1> &, + vector > &) const { + Assert (false, ExcNotUseful()); +}; + + + +template <> +void FECrissCross<1>::get_face_ansatz_points (const DoFHandler<1>::face_iterator &, + const Boundary<1> &, + vector > &) 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 > &, + vector &) const { + Assert (false, ExcNotUseful()); +}; + + + +template <> +void FECrissCross<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &, + const unsigned int, + const vector > &, + vector &) const { + Assert (false, ExcNotUseful()); +}; + + + +template <> +void FECrissCross<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, + const unsigned int, + const Boundary<1> &, + const vector > &, + vector > &) const { + Assert (false, ExcNotUseful()); +}; + + + +template <> +void FECrissCross<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, + const unsigned int, + const unsigned int, + const vector > &, + vector > &) const { + Assert (false, ExcNotUseful()); +}; + + + +template <> +void FECrissCross<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &, + const vector > &, + vector &, + const bool, + vector > &, + const bool, + vector > &, + const bool, + const dFMatrix &, + const vector > > &, + 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((i1-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) + + ((y1-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(-1,-1) : Point<2>(0,0)); + case 1: return ((y(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(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 > &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 > &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::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 > &ansatz_points) const { + const unsigned int dim = 2; + + Assert ((ansatz_points.size() == dofs_per_face) && + (ansatz_points.size() == GeometryInfo::vertices_per_face), + ExcWrongFieldDimension (ansatz_points.size(), + GeometryInfo::vertices_per_face)); + + for (unsigned int vertex=0; vertexvertex(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 > &unit_points, + vector &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 > &unit_points, + vector &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 > &unit_points, + vector > &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 > &unit_points, + vector > &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 +void FECrissCross::fill_fe_values (const DoFHandler::cell_iterator &cell, + const vector > &unit_points, + vector &jacobians, + const bool compute_jacobians, + vector > &ansatz_points, + const bool, + vector > &q_points, + const bool compute_q_points, + const dFMatrix &shape_values_transform, + const vector > > &shape_grad_transform, + const Boundary &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 (); + + // 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 gradient = shape_grad_transform[s][l]; + for (unsigned int i=0; i +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; +template class QCrissCross1;