+/**
+ * Abstract base class for concrete finite elements which use a
+ * (bi-,tri-)linear mapping from the unit cell to the real cell. Some
+ * functions can be singled out from these elements and are collected
+ * in this one.
+ */
+template <int dim>
+class FELinearMapping : public FiniteElement<dim> {
+ public:
+ /**
+ * Constructor. Simply passes through
+ * its arguments to the base class.
+ */
+ FELinearMapping (const unsigned int dofs_per_vertex,
+ const unsigned int dofs_per_line,
+ const unsigned int dofs_per_quad=0) :
+ FiniteElement<dim> (dofs_per_vertex,
+ dofs_per_line,
+ dofs_per_quad) {};
+
+ /**
+ * 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 Boundary<dim> &boundary) const;
+
+ protected:
+ /**
+ * 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.
+ */
+ double linear_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.
+ * Here, the (bi-)linear basis functions
+ * are meant, which are used for the
+ * computation of the transformation from
+ * unit cell to real space cell.
+ */
+ Point<dim> linear_shape_grad(const unsigned int i,
+ const Point<dim>& p) const;
+};
+
+
+
/*---------------------------- fe.h ---------------------------*/
/* end of #ifndef __fe_H */
* @author Wolfgang Bangerth, 1998
*/
template <int dim>
-class FELinear : public FiniteElement<dim> {
+class FELinear : public FELinearMapping<dim> {
public:
/**
* Constructor
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.
- *
- * 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 straight 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 Boundary<dim> &boundary) const;
-
/**
* Refer to the base class for detailed
* information on this function.
const Boundary<dim> &boundary,
vector<Point<dim> > &ansatz_points) 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.
* to the real cell is implemented.
*/
template <int dim>
-class FEQuadraticSub : public FiniteElement<dim> {
+class FEQuadraticSub : public FELinearMapping<dim> {
public:
/**
* Constructor
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.
- *
- * For one dimensional elements, this
- * function simply passes through to
- * the one implemented in the base class.
- */
- 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 Boundary<dim> &boundary) const;
-
/**
* Refer to the base class for detailed
* information on this function.
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_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.
- */
- 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;
-
- /**
- * 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;
-
- /**
- * 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 subface_no,
- const unsigned int face_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.
virtual void get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &cell,
const Boundary<dim> &boundary,
dFMatrix &local_mass_matrix) const;
-
- private:
- /**
- * 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.
- */
- double linear_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.
- * Here, the (bi-)linear basis functions
- * are meant, which are used for the
- * computation of the transformation from
- * unit cell to real space cell.
- */
- Point<dim> linear_shape_grad(const unsigned int i,
- const Point<dim>& p) const;
};
* freedom.
*/
template <int dim>
-class FECubicSub : public FiniteElement<dim> {
+class FECubicSub : public FELinearMapping<dim> {
public:
/**
* Constructor
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.
- *
- * For one dimensional elements, this
- * function simply passes through to
- * the one implemented in the base class.
- */
- 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 Boundary<dim> &boundary) const;
-
/**
* Refer to the base class for detailed
* information on this function.
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_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.
- */
- 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;
-
- /**
- * 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;
-
- /**
- * 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 subface_no,
- const unsigned int face_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.
virtual void get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &cell,
const Boundary<dim> &boundary,
dFMatrix &local_mass_matrix) const;
-
- private:
- /**
- * 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.
- */
- double linear_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.
- * Here, the (bi-)linear basis functions
- * are meant, which are used for the
- * computation of the transformation from
- * unit cell to real space cell.
- */
- Point<dim> linear_shape_grad(const unsigned int i,
- const Point<dim>& p) const;
};
* freedom.
*/
template <int dim>
-class FEQuarticSub : public FiniteElement<dim> {
+class FEQuarticSub : public FELinearMapping<dim> {
public:
/**
* Constructor
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.
- *
- * For one dimensional elements, this
- * function simply passes through to
- * the one implemented in the base class.
- */
- 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 Boundary<dim> &boundary) const;
-
/**
* Refer to the base class for detailed
* information on this function.
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_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.
- */
- 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;
-
- /**
- * 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;
-
- /**
- * 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 subface_no,
- const unsigned int face_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.
virtual void get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &cell,
const Boundary<dim> &boundary,
dFMatrix &local_mass_matrix) const;
-
- private:
- /**
- * 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.
- */
- double linear_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.
- * Here, the (bi-)linear basis functions
- * are meant, which are used for the
- * computation of the transformation from
- * unit cell to real space cell.
- */
- Point<dim> linear_shape_grad(const unsigned int i,
- const Point<dim>& p) const;
};
+#if deal_II_dimension == 1
+
+template <>
+inline
+double
+FELinearMapping<1>::linear_shape_value(const unsigned int i,
+ const Point<1> &p) const
+{
+ Assert((i<2), ExcInvalidIndex(i));
+ const double xi = p(0);
+ switch (i)
+ {
+ case 0: return 1.-xi;
+ case 1: return xi;
+ }
+ return 0.;
+};
+
+
+
+template <>
+inline
+Point<1>
+FELinearMapping<1>::linear_shape_grad(const unsigned int i,
+ const Point<1>&) const
+{
+ Assert((i<2), ExcInvalidIndex(i));
+ switch (i)
+ {
+ case 0: return Point<1>(-1.);
+ case 1: return Point<1>(1.);
+ }
+ return Point<1>();
+};
+
+
+
+template <>
+void FEQuadraticSub<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
+ const Boundary<1> &,
+ const vector<Point<0> > &,
+ vector<double> &) const {
+ Assert (false, ExcInternalError());
+};
+
+
+
+template <>
+void FEQuadraticSub<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
+ const unsigned int ,
+ const vector<Point<0> > &,
+ vector<double> &) const {
+ Assert (false, ExcInternalError());
+};
+
+
+
+template <>
+void FEQuadraticSub<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, ExcInternalError());
+};
+
+
+
+template <>
+void FEQuadraticSub<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, ExcInternalError());
+};
+
+
+template <>
+void FELinearMapping<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell,
+ const vector<Point<1> > &unit_points,
+ vector<dFMatrix> &jacobians,
+ const bool compute_jacobians,
+ vector<Point<1> > &ansatz_points,
+ const bool compute_ansatz_points,
+ vector<Point<1> > &q_points,
+ const bool compute_q_points,
+ const Boundary<1> &boundary) const {
+ // simply pass down
+ FiniteElement<1>::fill_fe_values (cell, unit_points,
+ jacobians, compute_jacobians,
+ ansatz_points, compute_ansatz_points,
+ q_points, compute_q_points, boundary);
+};
+
+
+
+#endif
+
+
+
+#if deal_II_dimension == 2
+
+template <>
+inline
+double
+FELinearMapping<2>::linear_shape_value (const unsigned int i,
+ const Point<2>& p) const
+{
+ Assert((i<4), ExcInvalidIndex(i));
+ switch (i)
+ {
+ case 0: return (1.-p(0)) * (1.-p(1));
+ case 1: return p(0) * (1.-p(1));
+ case 2: return p(0) * p(1);
+ case 3: return (1.-p(0)) * p(1);
+ }
+ return 0.;
+};
+
+
+
+template <>
+inline
+Point<2>
+FELinearMapping<2>::linear_shape_grad (const unsigned int i,
+ const Point<2>& p) const
+{
+ Assert((i<4), ExcInvalidIndex(i));
+ switch (i)
+ {
+ case 0: return Point<2> (p(1)-1., p(0)-1.);
+ case 1: return Point<2> (1.-p(1), -p(0));
+ case 2: return Point<2> (p(1), p(0));
+ case 3: return Point<2> (-p(1), 1.-p(0));
+ }
+ return Point<2> ();
+};
+
+
+
+template <>
+void FELinearMapping<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 FELinearMapping<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 FELinearMapping<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 FELinearMapping<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())));
+};
+
+#endif
+
+
+
+template <int dim>
+void FELinearMapping<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 compute_ansatz_points,
+ vector<Point<dim> > &q_points,
+ const bool compute_q_points,
+ 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();
+
+ Point<dim> vertices[GeometryInfo<dim>::vertices_per_cell];
+ for (unsigned int l=0; l<GeometryInfo<dim>::vertices_per_cell; ++l)
+ vertices[l] = cell->vertex(l);
+
+
+ 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)
+ //
+ // Here, N_j is the *linear* basis function,
+ // not that of the finite element, since we
+ // use a subparametric mapping
+ for (unsigned int j=0; j<GeometryInfo<dim>::vertices_per_cell; ++j)
+ for (unsigned int l=0; l<n_points; ++l)
+ q_points[l] += vertices[j] * linear_shape_value(j, unit_points[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 vertex 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_vertices
+ 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<GeometryInfo<dim>::vertices_per_cell; ++s)
+ {
+ // we want the linear transform,
+ // so use that function
+ const Point<dim> gradient = linear_shape_grad (s, unit_points[l]);
+ for (unsigned int i=0; i<dim; ++i)
+ for (unsigned int j=0; j<dim; ++j)
+ M(i,j) += vertices[s](i) * gradient(j);
+ };
+ jacobians[l].invert(M);
+ };
+ };
+
+ // compute ansatz points, which are
+ // the corners for linear elements
+ if (compute_ansatz_points)
+ get_ansatz_points (cell, boundary, ansatz_points);
+};
+
+
+
/*------------------------------- Explicit Instantiations -------------*/
template class FiniteElementData<deal_II_dimension>;
template class FiniteElementBase<deal_II_dimension>;
template class FiniteElement<deal_II_dimension>;
+template class FELinearMapping<deal_II_dimension>;
template <>
FECubicSub<1>::FECubicSub () :
- FiniteElement<1> (1, 2) {
+ FELinearMapping<1> (1, 2) {
prolongation[0](0,0) = 1.0;
prolongation[0](0,1) = 0.0;
prolongation[0](0,2) = 0.0;
-template <>
-void FECubicSub<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell,
- const vector<Point<1> > &unit_points,
- vector<dFMatrix> &jacobians,
- const bool compute_jacobians,
- vector<Point<1> > &ansatz_points,
- const bool compute_ansatz_points,
- vector<Point<1> > &q_points,
- const bool compute_q_points,
- const Boundary<1> &boundary) const {
- // simply pass down
- FiniteElement<1>::fill_fe_values (cell, unit_points,
- jacobians, compute_jacobians,
- ansatz_points, compute_ansatz_points,
- q_points, compute_q_points, boundary);
-};
-
-
-
template <>
double
FECubicSub<1>::shape_value(const unsigned int i,
-template <>
-inline
-double
-FECubicSub<1>::linear_shape_value(const unsigned int i,
- const Point<1> &p) const
-{
- Assert((i<2), ExcInvalidIndex(i));
- const double xi = p(0);
- switch (i)
- {
- case 0: return 1.-xi;
- case 1: return xi;
- }
- return 0.;
-};
-
-
-
template <>
Point<1>
FECubicSub<1>::shape_grad(const unsigned int i,
-template <>
-inline
-Point<1>
-FECubicSub<1>::linear_shape_grad(const unsigned int i,
- const Point<1>&) const
-{
- Assert((i<2), ExcInvalidIndex(i));
- switch (i)
- {
- case 0: return Point<1>(-1.);
- case 1: return Point<1>(1.);
- }
- return Point<1>();
-};
-
-
-
template <>
void FECubicSub<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell,
const Boundary<1> &boundary,
-template <>
-void FECubicSub<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
- const Boundary<1> &,
- const vector<Point<0> > &,
- vector<double> &) const {
- Assert (false, ExcInternalError());
-};
-
-
-
-template <>
-void FECubicSub<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
- const unsigned int ,
- const vector<Point<0> > &,
- vector<double> &) const {
- Assert (false, ExcInternalError());
-};
-
-
-
-template <>
-void FECubicSub<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, ExcInternalError());
-};
-
-
-
-template <>
-void FECubicSub<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, ExcInternalError());
-};
-
-
-
template <>
void FECubicSub<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cell,
const Boundary<1> &,
template <>
FECubicSub<2>::FECubicSub () :
- FiniteElement<2> (1, 2, 4)
+ FELinearMapping<2> (1, 2, 4)
{
interface_constraints(0,0) = -1.0/16.0;
interface_constraints(0,1) = -1.0/16.0;
-template <>
-inline
-double
-FECubicSub<2>::linear_shape_value (const unsigned int i,
- const Point<2>& p) const
-{
- Assert((i<4), ExcInvalidIndex(i));
- switch (i)
- {
- case 0: return (1.-p(0)) * (1.-p(1));
- case 1: return p(0) * (1.-p(1));
- case 2: return p(0) * p(1);
- case 3: return (1.-p(0)) * p(1);
- }
- return 0.;
-};
-
-
-
template <>
Point<2>
FECubicSub<2>::shape_grad (const unsigned int i,
-template <>
-inline
-Point<2>
-FECubicSub<2>::linear_shape_grad (const unsigned int i,
- const Point<2>& p) const
-{
- Assert((i<4), ExcInvalidIndex(i));
- switch (i)
- {
- case 0: return Point<2> (p(1)-1., p(0)-1.);
- case 1: return Point<2> (1.-p(1), -p(0));
- case 2: return Point<2> (p(1), p(0));
- case 3: return Point<2> (-p(1), 1.-p(0));
- }
- return Point<2> ();
-};
-
-
-
template <>
void FECubicSub<2>::get_local_mass_matrix (const DoFHandler<2>::cell_iterator &cell,
const Boundary<2> &,
-template <>
-void FECubicSub<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 FECubicSub<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 FECubicSub<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 FECubicSub<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())));
-};
-
#endif
-template <int dim>
-void FECubicSub<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 compute_ansatz_points,
- vector<Point<dim> > &q_points,
- const bool compute_q_points,
- 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();
-
- Point<dim> vertices[GeometryInfo<dim>::vertices_per_cell];
- for (unsigned int l=0; l<GeometryInfo<dim>::vertices_per_cell; ++l)
- vertices[l] = cell->vertex(l);
-
-
- 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)
- //
- // Here, N_j is the *linear* basis function,
- // not that of the finite element, since we
- // use a subparametric mapping
- for (unsigned int j=0; j<GeometryInfo<dim>::vertices_per_cell; ++j)
- for (unsigned int l=0; l<n_points; ++l)
- q_points[l] += vertices[j] * linear_shape_value(j, unit_points[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 vertex shape function at the
- quadrature point l (linear shape functions implied, as these
- are used for the mapping).
-
- We could therefore write:
- l=0..n_points-1
- i=0..dim-1
- j=0..dim-1
- M_{ij}(l) = 0
- s=0..n_vertices
- 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<GeometryInfo<dim>::vertices_per_cell; ++s)
- {
- // we want the linear transform,
- // so use that function
- const Point<dim> gradient = linear_shape_grad (s, unit_points[l]);
- for (unsigned int i=0; i<dim; ++i)
- for (unsigned int j=0; j<dim; ++j)
- M(i,j) += vertices[s](i) * gradient(j);
- };
- jacobians[l].invert(M);
- };
- };
-
- // compute ansatz points, which are
- // the corners for linear elements
- if (compute_ansatz_points)
- get_ansatz_points (cell, boundary, ansatz_points);
-};
-
-
-
// explicit instantiations
template <>
FELinear<1>::FELinear () :
- FiniteElement<1> (1, 0)
+ FELinearMapping<1> (1, 0)
{
// for restriction and prolongation matrices:
// note that we do not add up all the
-template <>
-void FELinear<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell,
- const vector<Point<1> > &unit_points,
- vector<dFMatrix> &jacobians,
- const bool compute_jacobians,
- vector<Point<1> > &ansatz_points,
- const bool compute_ansatz_points,
- vector<Point<1> > &q_points,
- const bool compute_q_points,
- const Boundary<1> &boundary) const {
- // simply pass down
- FiniteElement<1>::fill_fe_values (cell, unit_points,
- jacobians, compute_jacobians,
- ansatz_points, compute_ansatz_points,
- q_points, compute_q_points, boundary);
-};
-
-
-
template <>
void FELinear<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell,
const Boundary<1> &boundary,
-template <>
-void FELinear<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
- const Boundary<1> &,
- const vector<Point<0> > &,
- vector<double> &) const {
- Assert (false, ExcInternalError());
-};
-
-
-
-template <>
-void FELinear<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
- const unsigned int ,
- const vector<Point<0> > &,
- vector<double> &) const {
- Assert (false, ExcInternalError());
-};
-
-
-
-template <>
-void FELinear<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, ExcInternalError());
-};
-
-
-
-template <>
-void FELinear<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, ExcInternalError());
-};
-
-
-
template <>
void FELinear<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cell,
const Boundary<1> &,
template <>
FELinear<2>::FELinear () :
- FiniteElement<2> (1, 0, 0)
+ FELinearMapping<2> (1, 0, 0)
{
interface_constraints(0,0) = 1./2.;
interface_constraints(0,1) = 1./2.;
-template <>
-void FELinear<2>::get_face_jacobians (const DoFHandler<2>::face_iterator &face,
- const Boundary<2> &,
- const vector<Point<1> > &unit_points,
- vector<double> &face_jacobians) const {
- 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 FELinear<2>::get_subface_jacobians (const DoFHandler<2>::face_iterator &face,
- const unsigned int,
- const vector<Point<1> > &unit_points,
- vector<double> &face_jacobians) const {
- 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 FELinear<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 {
- 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 FELinear<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 {
- // 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())));
-};
-
#endif
-template <int dim>
-void FELinear<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 compute_ansatz_points,
- vector<Point<dim> > &q_points,
- const bool compute_q_points,
- 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();
-
- Point<dim> vertices[GeometryInfo<dim>::vertices_per_cell];
- for (unsigned int l=0; l<GeometryInfo<dim>::vertices_per_cell; ++l)
- vertices[l] = cell->vertex(l);
-
-
- 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<GeometryInfo<dim>::vertices_per_cell; ++j)
- for (unsigned int l=0; l<n_points; ++l)
- q_points[l] += vertices[j] *
- FELinear<dim>::shape_value(j, unit_points[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 vertex 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_vertices
- 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<GeometryInfo<dim>::vertices_per_cell; ++s)
- {
- // we want a linear transform and
- // if we prepend the class name in
- // front of the #shape_grad#, we
- // need not use virtual function
- // calls.
- const Point<dim> gradient
- = FELinear<dim>::shape_grad (s, unit_points[l]);
- for (unsigned int i=0; i<dim; ++i)
- for (unsigned int j=0; j<dim; ++j)
- M(i,j) += vertices[s](i) * gradient(j);
- };
- jacobians[l].invert(M);
- };
- };
-
- // compute ansatz points, which are
- // the corners for linear elements
- if (compute_ansatz_points)
- get_ansatz_points (cell, boundary, ansatz_points);
-};
-
-
-
template <int dim>
void FELinear<dim>::get_ansatz_points (const typename DoFHandler<dim>::cell_iterator &cell,
const Boundary<dim> &,
template <>
FEQuadraticSub<1>::FEQuadraticSub () :
- FiniteElement<1> (1, 1) {
+ FELinearMapping<1> (1, 1) {
/*
Get the prolongation matrices by the following little maple script:
-template <>
-void FEQuadraticSub<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell,
- const vector<Point<1> > &unit_points,
- vector<dFMatrix> &jacobians,
- const bool compute_jacobians,
- vector<Point<1> > &ansatz_points,
- const bool compute_ansatz_points,
- vector<Point<1> > &q_points,
- const bool compute_q_points,
- const Boundary<1> &boundary) const {
- // simply pass down
- FiniteElement<1>::fill_fe_values (cell, unit_points,
- jacobians, compute_jacobians,
- ansatz_points, compute_ansatz_points,
- q_points, compute_q_points, boundary);
-};
-
-
-
template <>
double
FEQuadraticSub<1>::shape_value(const unsigned int i,
-template <>
-inline
-double
-FEQuadraticSub<1>::linear_shape_value(const unsigned int i,
- const Point<1> &p) const
-{
- Assert((i<2), ExcInvalidIndex(i));
- const double xi = p(0);
- switch (i)
- {
- case 0: return 1.-xi;
- case 1: return xi;
- }
- return 0.;
-};
-
-
-
template <>
Point<1>
FEQuadraticSub<1>::shape_grad(const unsigned int i,
-template <>
-inline
-Point<1>
-FEQuadraticSub<1>::linear_shape_grad(const unsigned int i,
- const Point<1>&) const
-{
- Assert((i<2), ExcInvalidIndex(i));
- switch (i)
- {
- case 0: return Point<1>(-1.);
- case 1: return Point<1>(1.);
- }
- return Point<1>();
-};
-
-
-
template <>
void FEQuadraticSub<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell,
const Boundary<1> &boundary,
-template <>
-void FEQuadraticSub<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
- const Boundary<1> &,
- const vector<Point<0> > &,
- vector<double> &) const {
- Assert (false, ExcInternalError());
-};
-
-
-
-template <>
-void FEQuadraticSub<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
- const unsigned int ,
- const vector<Point<0> > &,
- vector<double> &) const {
- Assert (false, ExcInternalError());
-};
-
-
-
-template <>
-void FEQuadraticSub<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, ExcInternalError());
-};
-
-
-
-template <>
-void FEQuadraticSub<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, ExcInternalError());
-};
-
-
-
template <>
void FEQuadraticSub<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cell,
const Boundary<1> &,
template <>
FEQuadraticSub<2>::FEQuadraticSub () :
- FiniteElement<2> (1, 1, 1)
+ FELinearMapping<2> (1, 1, 1)
{
interface_constraints(0,2) = 1.0;
interface_constraints(1,0) = 3./8.;
-template <>
-inline
-double
-FEQuadraticSub<2>::linear_shape_value (const unsigned int i,
- const Point<2>& p) const
-{
- Assert((i<4), ExcInvalidIndex(i));
- switch (i)
- {
- case 0: return (1.-p(0)) * (1.-p(1));
- case 1: return p(0) * (1.-p(1));
- case 2: return p(0) * p(1);
- case 3: return (1.-p(0)) * p(1);
- }
- return 0.;
-};
-
-
-
template <>
Point<2>
FEQuadraticSub<2>::shape_grad (const unsigned int i,
-template <>
-inline
-Point<2>
-FEQuadraticSub<2>::linear_shape_grad (const unsigned int i,
- const Point<2>& p) const
-{
- Assert((i<4), ExcInvalidIndex(i));
- switch (i)
- {
- case 0: return Point<2> (p(1)-1., p(0)-1.);
- case 1: return Point<2> (1.-p(1), -p(0));
- case 2: return Point<2> (p(1), p(0));
- case 3: return Point<2> (-p(1), 1.-p(0));
- }
- return Point<2> ();
-};
-
-
-
template <>
void FEQuadraticSub<2>::get_local_mass_matrix (const DoFHandler<2>::cell_iterator &cell,
const Boundary<2> &,
-template <>
-void FEQuadraticSub<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 FEQuadraticSub<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 FEQuadraticSub<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 FEQuadraticSub<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())));
-};
-
#endif
-template <int dim>
-void FEQuadraticSub<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 compute_ansatz_points,
- vector<Point<dim> > &q_points,
- const bool compute_q_points,
- 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();
-
- Point<dim> vertices[GeometryInfo<dim>::vertices_per_cell];
- for (unsigned int l=0; l<GeometryInfo<dim>::vertices_per_cell; ++l)
- vertices[l] = cell->vertex(l);
-
-
- 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)
- //
- // Here, N_j is the *linear* basis function,
- // not that of the finite element, since we
- // use a subparametric mapping
- for (unsigned int j=0; j<GeometryInfo<dim>::vertices_per_cell; ++j)
- for (unsigned int l=0; l<n_points; ++l)
- q_points[l] += vertices[j] * linear_shape_value(j, unit_points[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 vertex 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_vertices
- 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<GeometryInfo<dim>::vertices_per_cell; ++s)
- {
- // we want the linear transform,
- // so use that function
- const Point<dim> gradient = linear_shape_grad (s, unit_points[l]);
- for (unsigned int i=0; i<dim; ++i)
- for (unsigned int j=0; j<dim; ++j)
- M(i,j) += vertices[s](i) * gradient(j);
- };
- jacobians[l].invert(M);
- };
- };
-
- // compute ansatz points, which are
- // the corners for linear elements
- if (compute_ansatz_points)
- get_ansatz_points (cell, boundary, ansatz_points);
-};
-
-
-
-
// explicit instantiations
template class FEQuadraticSub<deal_II_dimension>;
template <>
FEQuarticSub<1>::FEQuarticSub () :
- FiniteElement<1> (1, 3) {
+ FELinearMapping<1> (1, 3) {
prolongation[0](0,0) = 1.0;
prolongation[0](1,3) = 1.0;
prolongation[0](2,0) = 35.0/128.0;
-template <>
-void FEQuarticSub<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell,
- const vector<Point<1> > &unit_points,
- vector<dFMatrix> &jacobians,
- const bool compute_jacobians,
- vector<Point<1> > &ansatz_points,
- const bool compute_ansatz_points,
- vector<Point<1> > &q_points,
- const bool compute_q_points,
- const Boundary<1> &boundary) const {
- // simply pass down
- FiniteElement<1>::fill_fe_values (cell, unit_points,
- jacobians, compute_jacobians,
- ansatz_points, compute_ansatz_points,
- q_points, compute_q_points, boundary);
-};
-
-
-
template <>
double
FEQuarticSub<1>::shape_value(const unsigned int i,
-template <>
-inline
-double
-FEQuarticSub<1>::linear_shape_value(const unsigned int i,
- const Point<1> &p) const
-{
- Assert((i<2), ExcInvalidIndex(i));
- const double xi = p(0);
- switch (i)
- {
- case 0: return 1.-xi;
- case 1: return xi;
- }
- return 0.;
-};
-
-
-
template <>
Point<1>
FEQuarticSub<1>::shape_grad(const unsigned int i,
-template <>
-inline
-Point<1>
-FEQuarticSub<1>::linear_shape_grad(const unsigned int i,
- const Point<1>&) const
-{
- Assert((i<2), ExcInvalidIndex(i));
- switch (i)
- {
- case 0: return Point<1>(-1.);
- case 1: return Point<1>(1.);
- }
- return Point<1>();
-};
-
-
-
template <>
void FEQuarticSub<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell,
const Boundary<1> &boundary,
-template <>
-void FEQuarticSub<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
- const Boundary<1> &,
- const vector<Point<0> > &,
- vector<double> &) const {
- Assert (false, ExcInternalError());
-};
-
-
-
-template <>
-void FEQuarticSub<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
- const unsigned int ,
- const vector<Point<0> > &,
- vector<double> &) const {
- Assert (false, ExcInternalError());
-};
-
-
-
-template <>
-void FEQuarticSub<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, ExcInternalError());
-};
-
-
-
-template <>
-void FEQuarticSub<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, ExcInternalError());
-};
-
-
-
template <>
void FEQuarticSub<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cell,
const Boundary<1> &,
template <>
FEQuarticSub<2>::FEQuarticSub () :
- FiniteElement<2> (1, 3, 9)
+ FELinearMapping<2> (1, 3, 9)
{
interface_constraints(0,3) = 1.0;
interface_constraints(1,0) = 35.0/128.0;
-template <>
-inline
-double
-FEQuarticSub<2>::linear_shape_value (const unsigned int i,
- const Point<2>& p) const
-{
- Assert((i<4), ExcInvalidIndex(i));
- switch (i)
- {
- case 0: return (1.-p(0)) * (1.-p(1));
- case 1: return p(0) * (1.-p(1));
- case 2: return p(0) * p(1);
- case 3: return (1.-p(0)) * p(1);
- }
- return 0.;
-};
-
-
-
template <>
Point<2>
FEQuarticSub<2>::shape_grad (const unsigned int i,
-template <>
-inline
-Point<2>
-FEQuarticSub<2>::linear_shape_grad (const unsigned int i,
- const Point<2>& p) const
-{
- Assert((i<4), ExcInvalidIndex(i));
- switch (i)
- {
- case 0: return Point<2> (p(1)-1., p(0)-1.);
- case 1: return Point<2> (1.-p(1), -p(0));
- case 2: return Point<2> (p(1), p(0));
- case 3: return Point<2> (-p(1), 1.-p(0));
- }
- return Point<2> ();
-};
-
-
-
template <>
void FEQuarticSub<2>::get_local_mass_matrix (const DoFHandler<2>::cell_iterator &cell,
const Boundary<2> &,
-template <>
-void FEQuarticSub<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 FEQuarticSub<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 FEQuarticSub<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 FEQuarticSub<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())));
-};
-
#endif
-template <int dim>
-void FEQuarticSub<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 compute_ansatz_points,
- vector<Point<dim> > &q_points,
- const bool compute_q_points,
- 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();
-
- Point<dim> vertices[GeometryInfo<dim>::vertices_per_cell];
- for (unsigned int l=0; l<GeometryInfo<dim>::vertices_per_cell; ++l)
- vertices[l] = cell->vertex(l);
-
-
- 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)
- //
- // Here, N_j is the *linear* basis function,
- // not that of the finite element, since we
- // use a subparametric mapping
- for (unsigned int j=0; j<GeometryInfo<dim>::vertices_per_cell; ++j)
- for (unsigned int l=0; l<n_points; ++l)
- q_points[l] += vertices[j] * linear_shape_value(j, unit_points[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 vertex shape function at the
- quadrature point l (linear shape functions implied, as these
- are used for the mapping).
-
- We could therefore write:
- l=0..n_points-1
- i=0..dim-1
- j=0..dim-1
- M_{ij}(l) = 0
- s=0..n_vertices
- 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<GeometryInfo<dim>::vertices_per_cell; ++s)
- {
- // we want the linear transform,
- // so use that function
- const Point<dim> gradient = linear_shape_grad (s, unit_points[l]);
- for (unsigned int i=0; i<dim; ++i)
- for (unsigned int j=0; j<dim; ++j)
- M(i,j) += vertices[s](i) * gradient(j);
- };
- jacobians[l].invert(M);
- };
- };
-
- // compute ansatz points, which are
- // the corners for linear elements
- if (compute_ansatz_points)
- get_ansatz_points (cell, boundary, ansatz_points);
-};
-
-
// explicit instantiations