* performed by the specific finite element
* class.
*
- * Two fields remain to be finite element
+ * Three fields remain to be finite element
* specific in this standard implementation:
* The jacobi determinants of the
* transformation from the unit face to the
- * real face and the ansatz points. For
- * these two fields, there exist pure
- * virtual functions, #get_face_jacobians#
- * and the #get_face_ansatz_points#
- * function.
+ * real face, the ansatz points
+ * and the outward normal vectors. For
+ * these fields, there exist pure
+ * virtual functions, #get_face_jacobians#,
+ * #get_face_ansatz_points# and
+ * #get_normal_vectors#.
*
* Though there is a standard
* implementation, there
const bool compute_q_points,
vector<double> &face_jacobi_determinants,
const bool compute_face_jacobians,
+ vector<Point<dim> > &normal_vectors,
+ const bool compute_normal_vectors,
const Boundary<dim> &boundary) const;
/**
const Boundary<dim> &boundary,
const vector<Point<dim-1> > &unit_points,
vector<double> &face_jacobi_determinants) const =0;
+
+ /**
+ * Compute the normal vectors to the cell
+ * at the quadrature points. See the
+ * documentation for the #fill_fe_face_values#
+ * function for more details. The function
+ * must guarantee that the length of the
+ * vectors be one.
+ *
+ * Since any implementation for one
+ * dimension would be senseless, all
+ * derived classes should throw an error
+ * when called with #dim==1#.
+ */
+ 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 =0;
/**
* Exception
const Boundary<dim> &boundary,
const vector<Point<dim-1> > &unit_points,
vector<double> &face_jacobi_determinants) const;
+
+ /**
+ * 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;
};
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_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;
};
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_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;
};
* on which the ansatz functions are
* located.
*/
- update_ansatz_points = 16
+ update_ansatz_points = 16,
+ /**
+ * Update the outward normal vectors
+ * to the face relative to this cell.
+ * This flag is only evaluated by
+ * the #FEFaceValues# class.
+ */
+ update_normal_vectors = 32
};
-
/*---------------------------- fe_update_flags.h ---------------------------*/
/* end of #ifndef __fe_update_flags_H */
#endif
* segments, but higher order elements
* may use other ways.)
*/
- void reinit (const DoFHandler<dim>::cell_iterator &,
+ void reinit (const typename DoFHandler<dim>::cell_iterator &,
const FiniteElement<dim> &,
const Boundary<dim> &);
for the transformation of the unit face to the real face (needed to
compute the weight factors for integration along faces). These two
concepts have to be carefully separated.
+
+ Finally, we will often need the outward normal to a cell at the quadrature
+ points. While this could in principle be easily done using the Jacobi
+ matrices at the quadrature points and the normal vectors to the unit cell
+ (also easily derived, since they have an appealingly easy form for the unit
+ cell ;-), it is more efficiently done by the finite element class itself.
+ For example for (bi-, tri-)linear mappings the normal vector is readily
+ available without compicated matrix-vector-multiplications.
*/
template <int dim>
class FEFaceValues {
* alike cells.
*/
const vector<double> & get_JxW_values () const;
+
+ /**
+ * Return the outward normal vector to
+ * the cell at the #i#th quadrature
+ * point. The length of the vector
+ * is normalized to one.
+ */
+ const Point<dim> & normal_vector (const unsigned int i) const;
+
+ /**
+ * Return the list of outward normal
+ * vectors to the cell at the
+ * quadrature points.
+ */
+ const vector<Point<dim> > & get_normal_vectors () const;
/**
* Reinitialize the gradients, Jacobi
* segments, but higher order elements
* may use other ways.)
*/
- void reinit (const DoFHandler<dim>::cell_iterator &cell,
+ void reinit (const typename DoFHandler<dim>::cell_iterator &cell,
const unsigned int face_no,
const FiniteElement<dim> &fe,
const Boundary<dim> &boundary);
* compute the JxW values.
*/
vector<double> face_jacobi_determinants;
+
+ /**
+ * List of outward normal vectors at the
+ * quadrature points. This field is filled
+ * in by the finite element class.
+ */
+ vector<Point<dim> > normal_vectors;
/**
* Store which fields are to be updated by
+template <int dim>
+inline
+const vector<Point<dim> > &
+FEFaceValues<dim>::get_normal_vectors () const {
+ Assert (update_flags & update_normal_vectors, ExcAccessToUninitializedField());
+ return normal_vectors;
+};
+
+
+
/*---------------------------- fe_values.h ---------------------------*/
/* end of #ifndef __fe_values_H */
const bool ,
vector<double> &,
const bool ,
+ vector<Point<1> > &,
+ const bool,
const Boundary<1> &) const {
Assert (false, ExcNotImplemented());
}
const bool compute_q_points,
vector<double> &face_jacobi_determinants,
const bool compute_face_jacobians,
+ vector<Point<dim> > &normal_vectors,
+ const bool compute_normal_vectors,
const Boundary<dim> &boundary) const {
Assert (jacobians.size() == unit_points.size(),
ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
q_points, compute_q_points,
boundary);
- cout << "Global unit points:\n";
- for (unsigned int p=0; p<unit_points.size(); ++p)
- cout << " " << global_unit_points[p] << endl;
-
if (compute_ansatz_points)
get_face_ansatz_points (cell->face(face_no), boundary, ansatz_points);
if (compute_face_jacobians)
get_face_jacobians (cell->face(face_no), boundary,
unit_points, face_jacobi_determinants);
+
+ if (compute_normal_vectors)
+ get_normal_vectors (cell, face_no, boundary,
+ unit_points, normal_vectors);
+ cout << "Global unit points:\n";
+ for (unsigned int p=0; p<unit_points.size(); ++p)
+ cout << " " << global_unit_points[p] << endl;
+
cout << "Global ansatz points:\n";
for (unsigned int p=0; p<unit_points.size(); ++p)
cout << " " << ansatz_points[p] << endl;
+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());
+};
+
+
+
+
+
FELinear<2>::FELinear () :
FiniteElement<2> (1, 0, 0)
{
+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())));
+};
+
+
+
+void FEQuadratic<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());
+};
+
+
+
+
+
FEQuadratic<2>::FEQuadratic () :
FiniteElement<2> (1, 1, 1)
{
+template <int dim>
+void FEQuadratic<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &,
+ const unsigned int,
+ const Boundary<dim> &,
+ const vector<Point<dim-1> > &,
+ vector<Point<dim> > &) const {
+ Assert (false, ExcInternalError());
+};
+
+
+
+void FECubic<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());
+};
+
+
+
FECubic<2>::FECubic () :
FiniteElement<2> (1, 2, 4) {};
+template <int dim>
+void FECubic<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &,
+ const unsigned int,
+ const Boundary<dim> &,
+ const vector<Point<dim-1> > &,
+ vector<Point<dim> > &) const {
+ Assert (false, ExcInternalError());
+};
+
+
// explicit instantiations
ansatz_points (fe.total_dofs, Point<dim>()),
jacobi_matrices (quadrature.n_quadrature_points,dFMatrix(dim,dim)),
face_jacobi_determinants (quadrature.n_quadrature_points,0),
+ normal_vectors (quadrature.n_quadrature_points,Point<dim>()),
update_flags (update_flags),
selected_face(0)
{
+template <int dim>
+const Point<dim> & FEFaceValues<dim>::normal_vector (const unsigned int i) const {
+ Assert (i<normal_vectors.size(), ExcInvalidIndex(i, normal_vectors.size()));
+ Assert (update_flags & update_normal_vectors, ExcAccessToUninitializedField());
+
+ return normal_vectors[i];
+};
+
+
+
template <int dim>
double FEFaceValues<dim>::JxW (const unsigned int i) const {
Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
update_flags & update_q_points,
face_jacobi_determinants,
update_flags & update_JxW_values,
+ normal_vectors,
+ update_flags & update_normal_vectors,
boundary);
// compute gradients on real element if