const bool,
const Boundary<1> &) const {
Assert (false, ExcNotImplemented());
-}
+};
+
+void FiniteElement<1>::fill_fe_subface_values (const DoFHandler<1>::cell_iterator &,
+ const unsigned int ,
+ const unsigned int ,
+ const vector<Point<0> > &,
+ const vector<Point<1> > &,
+ vector<dFMatrix> &,
+ const bool ,
+ vector<Point<1> > &,
+ const bool ,
+ vector<double> &,
+ const bool ,
+ vector<Point<1> > &,
+ const bool,
+ const Boundary<1> &) const {
+ Assert (false, ExcNotImplemented());
+};
+
+
+
+template <int dim>
+void FiniteElement<dim>::fill_fe_subface_values (const DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const unsigned int subface_no,
+ const vector<Point<dim-1> > &unit_points,
+ const vector<Point<dim> > &global_unit_points,
+ vector<dFMatrix> &jacobians,
+ const bool compute_jacobians,
+ vector<Point<dim> > &q_points,
+ 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()));
+ Assert (q_points.size() == unit_points.size(),
+ ExcWrongFieldDimension(q_points.size(), unit_points.size()));
+ Assert (global_unit_points.size() == unit_points.size(),
+ ExcWrongFieldDimension(global_unit_points.size(), unit_points.size()));
+
+ static vector<Point<dim> > dummy(total_dofs);
+ fill_fe_values (cell, global_unit_points,
+ jacobians, compute_jacobians,
+ dummy, false,
+ q_points, compute_q_points,
+ boundary);
+
+ if (compute_face_jacobians)
+ get_subface_jacobians (cell->face(face_no), subface_no,
+ unit_points, face_jacobi_determinants);
+
+ if (compute_normal_vectors)
+ get_normal_vectors (cell, face_no, subface_no,
+ unit_points, normal_vectors);
+};
+
+
/*------------------------------- Explicit Instantiations -------------*/
template class FiniteElementData<1>;
+void FELinear<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
+ const unsigned int ,
+ const vector<Point<0> > &,
+ vector<double> &) const {
+ Assert (false, ExcInternalError());
+};
+
+
+
void FELinear<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
const unsigned int,
const Boundary<1> &,
+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());
+};
+
+
+
FELinear<2>::FELinear () :
+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);
+};
+
+
+
void FELinear<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell,
const unsigned int face_no,
const Boundary<2> &,
+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())));
+};
+
+
+
+void FEQuadratic<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
+ const unsigned int ,
+ const vector<Point<0> > &,
+ vector<double> &) const {
+ Assert (false, ExcInternalError());
+};
+
+
+
void FEQuadratic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
const unsigned int,
const Boundary<1> &,
+void FEQuadratic<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());
+};
+
+
FEQuadratic<2>::FEQuadratic () :
template <int dim>
void FEQuadratic<dim>::get_face_jacobians (const DoFHandler<dim>::face_iterator &,
- const Boundary<dim> &,
- const vector<Point<dim-1> > &,
- vector<double> &) const {
+ const Boundary<dim> &,
+ const vector<Point<dim-1> > &,
+ vector<double> &) const {
+ Assert (false, ExcNotImplemented());
+};
+
+
+
+template <int dim>
+void FEQuadratic<dim>::get_subface_jacobians (const DoFHandler<dim>::face_iterator &face,
+ const unsigned int ,
+ const vector<Point<dim-1> > &,
+ vector<double> &) const {
+ Assert (face->at_boundary() == false,
+ ExcBoundaryFaceUsed ());
+
Assert (false, ExcNotImplemented());
};
const Boundary<dim> &,
const vector<Point<dim-1> > &,
vector<Point<dim> > &) const {
- Assert (false, ExcInternalError());
+ Assert (false, ExcNotImplemented());
+};
+
+
+
+template <int dim>
+void FEQuadratic<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const unsigned int,
+ const vector<Point<dim-1> > &,
+ vector<Point<dim> > &) const {
+ Assert (cell->face(face_no)->at_boundary() == false,
+ ExcBoundaryFaceUsed ());
+ Assert (false, ExcNotImplemented());
};
+void FECubic<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
+ const unsigned int ,
+ const vector<Point<0> > &,
+ vector<double> &) const {
+ Assert (false, ExcInternalError());
+};
+
+
+
void FECubic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
const unsigned int,
const Boundary<1> &,
+void FECubic<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());
+};
+
+
FECubic<2>::FECubic () :
FiniteElement<2> (1, 2, 4) {};
+template <int dim>
+void FECubic<dim>::get_subface_jacobians (const DoFHandler<dim>::face_iterator &face,
+ const unsigned int ,
+ const vector<Point<dim-1> > &,
+ vector<double> &) const {
+ Assert (face->at_boundary() == false,
+ ExcBoundaryFaceUsed ());
+
+ Assert (false, ExcNotImplemented());
+};
+
+
+
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());
+ Assert (false, ExcNotImplemented());
+};
+
+
+
+template <int dim>
+void FECubic<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const unsigned int ,
+ const vector<Point<dim-1> > &,
+ vector<Point<dim> > &) const {
+ Assert (cell->face(face_no)->at_boundary() == false,
+ ExcBoundaryFaceUsed ());
+
+ Assert (false, ExcNotImplemented());
};
vector<Point<dim> >(quadrature.n_quadrature_points)),
unit_quadrature_points(quadrature.get_quad_points())
{
- Assert ((update_flags | update_normal_vectors) == false,
+ Assert ((update_flags & update_normal_vectors) == false,
ExcInvalidUpdateFlag());
for (unsigned int i=0; i<fe.total_dofs; ++i)
n_dofs,
n_faces_or_subfaces,
update_flags),
+ unit_shape_gradients (n_faces_or_subfaces,
+ vector<vector<Point<dim> > >(n_dofs,
+ vector<Point<dim> >(n_q_points))),
unit_face_quadrature_points (n_q_points, Point<dim-1>()),
unit_quadrature_points (n_faces_or_subfaces,
vector<Point<dim> >(n_q_points, Point<dim>())),
face_jacobi_determinants (n_q_points, 0),
normal_vectors (n_q_points)
-{
- for (unsigned int i=0; i<n_faces_or_subfaces; ++i)
- {
- unit_shape_gradients[i].resize (n_dofs,
- vector<Point<dim> >(n_q_points));
- unit_quadrature_points[i].resize (n_q_points,
- Point<dim>());
- };
-};
+{};
Assert (update_flags & update_jacobians, ExcCannotInitializeField());
for (unsigned int i=0; i<fe.total_dofs; ++i)
- for (unsigned int j=0; j<n_quadrature_points; ++j)
- {
- shape_gradients[i][j] = Point<dim>();
-
+ {
+ fill_n (shape_gradients[i].begin(),
+ n_quadrature_points,
+ Point<dim>());
+ for (unsigned int j=0; j<n_quadrature_points; ++j)
for (unsigned int s=0; s<dim; ++s)
// (grad psi)_s =
// (grad_{\xi\eta})_b J_{bs}
shape_gradients[i][j](s)
+= (unit_shape_gradients[face_no][i][j](b) *
jacobi_matrices[j](b,s));
- };
+ };
};
const unsigned int subface_no,
const FiniteElement<dim> &fe,
const Boundary<dim> &boundary) {
+ Assert (cell->face(face_no)->at_boundary() == false,
+ ExcReinitCalledWithBoundaryFace());
+
present_cell = cell;
selected_dataset = face_no*(1<<(dim-1)) + subface_no;
// fill jacobi matrices and real
{
Assert (update_flags & update_jacobians, ExcCannotInitializeField());
- for (unsigned int i=0; i<fe.total_dofs; ++i)
- for (unsigned int j=0; j<n_quadrature_points; ++j)
- {
- shape_gradients[i][j] = Point<dim>();
-
+ for (unsigned int i=0; i<fe.total_dofs; ++i)
+ {
+ fill_n (shape_gradients[i].begin(),
+ n_quadrature_points,
+ Point<dim>());
+ for (unsigned int j=0; j<n_quadrature_points; ++j)
for (unsigned int s=0; s<dim; ++s)
// (grad psi)_s =
// (grad_{\xi\eta})_b J_{bs}
shape_gradients[i][j](s)
+= (unit_shape_gradients[selected_dataset][i][j](b) *
jacobi_matrices[j](b,s));
- };
+ };
};