From cd6c54f8026be271940ffe0978ee91fba40dfb58 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 25 Jun 1998 11:20:23 +0000 Subject: [PATCH] Start of implementation of (bi-)quadratic elements. git-svn-id: https://svn.dealii.org/trunk@411 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_lib.lagrange.h | 6 +- deal.II/deal.II/source/fe/fe_lib.linear.cc | 365 ++++++++++++++----- 2 files changed, 271 insertions(+), 100 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe_lib.lagrange.h b/deal.II/deal.II/include/fe/fe_lib.lagrange.h index dda7a46759..1b3f403463 100644 --- a/deal.II/deal.II/include/fe/fe_lib.lagrange.h +++ b/deal.II/deal.II/include/fe/fe_lib.lagrange.h @@ -165,16 +165,16 @@ class FELinear : public FiniteElement { /** * Define a (bi-, tri-, etc)quadratic finite element in #dim# space dimensions. - * In one space dimension, a linear (subparametric) mapping from the unit cell + * A linear (subparametric) mapping from the unit cell * to the real cell is implemented. */ template -class FEQuadratic : public FiniteElement { +class FEQuadraticSub : public FiniteElement { public: /** * Constructor */ - FEQuadratic (); + FEQuadraticSub (); /** * Return the value of the #i#th shape diff --git a/deal.II/deal.II/source/fe/fe_lib.linear.cc b/deal.II/deal.II/source/fe/fe_lib.linear.cc index 0dc692538b..cfdf0bb204 100644 --- a/deal.II/deal.II/source/fe/fe_lib.linear.cc +++ b/deal.II/deal.II/source/fe/fe_lib.linear.cc @@ -47,7 +47,7 @@ FELinear<1>::FELinear () : template <> double FELinear<1>::shape_value(const unsigned int i, - const Point<1>& p) const + const Point<1> &p) const { Assert((i::get_face_ansatz_points (const typename DoFHandler::face #if deal_II_dimension == 1 template <> -FEQuadratic<1>::FEQuadratic () : +FEQuadraticSub<1>::FEQuadraticSub () : FiniteElement<1> (1, 1) {}; template <> -void FEQuadratic<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell, +void FEQuadraticSub<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell, const vector > &unit_points, vector &jacobians, const bool compute_jacobians, @@ -672,16 +672,53 @@ void FEQuadratic<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell, template <> -void FEQuadratic<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell, - const Boundary<1> &boundary, - vector > &ansatz_points) const { +double +FEQuadraticSub<1>::shape_value(const unsigned int i, + const Point<1> &p) const +{ + Assert((i +inline +Point<1> +FEQuadraticSub<1>::shape_grad(const unsigned int i, + const Point<1> &p) const +{ + Assert((i(-3+4*xi); + case 1: return Point<1>(4*xi-1); + case 2: return Point<1>(4-8*xi); + } + return Point<1>(); +}; + + + +template <> +void FEQuadraticSub<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell, + const Boundary<1> &boundary, + vector > &ansatz_points) const { FiniteElement<1>::get_ansatz_points (cell, boundary, ansatz_points); }; template <> -void FEQuadratic<1>::get_face_ansatz_points (const typename DoFHandler<1>::face_iterator &, +void FEQuadraticSub<1>::get_face_ansatz_points (const typename DoFHandler<1>::face_iterator &, const Boundary<1> &, vector > &) const { Assert (false, ExcInternalError()); @@ -690,42 +727,42 @@ void FEQuadratic<1>::get_face_ansatz_points (const typename DoFHandler<1>::face_ template <> -void FEQuadratic<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &, - const Boundary<1> &, - const vector > &, - vector &) const { +void FEQuadraticSub<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &, + const Boundary<1> &, + const vector > &, + vector &) const { Assert (false, ExcInternalError()); }; template <> -void FEQuadratic<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &, - const unsigned int , - const vector > &, - vector &) const { +void FEQuadraticSub<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &, + const unsigned int , + const vector > &, + vector &) const { Assert (false, ExcInternalError()); }; template <> -void FEQuadratic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, - const unsigned int, - const Boundary<1> &, - const vector > &, - vector > &) const { +void FEQuadraticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, + const unsigned int, + const Boundary<1> &, + const vector > &, + vector > &) const { Assert (false, ExcInternalError()); }; template <> -void FEQuadratic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, - const unsigned int, - const unsigned int, - const vector > &, - vector > &) const { +void FEQuadraticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, + const unsigned int, + const unsigned int, + const vector > &, + vector > &) const { Assert (false, ExcInternalError()); }; @@ -735,7 +772,7 @@ void FEQuadratic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, #if deal_II_dimension == 2 template <> -FEQuadratic<2>::FEQuadratic () : +FEQuadraticSub<2>::FEQuadraticSub () : FiniteElement<2> (1, 1, 1) { interface_constraints(0,2) = 1.0; @@ -751,123 +788,257 @@ FEQuadratic<2>::FEQuadratic () : Assert (false, ExcNotImplemented()); }; -#endif - - -template +template <> double -FEQuadratic::shape_value (const unsigned int i, - const Point &) const +FEQuadraticSub<2>::shape_value (const unsigned int i, + const Point<2> &p) const { - Assert (i::ExcInvalidIndex(i)); - Assert (false, ExcNotImplemented()); - return 0.; + Assert (i -Point -FEQuadratic::shape_grad (const unsigned int i, - const Point &) const +template <> +Point<2> +FEQuadraticSub<2>::shape_grad (const unsigned int i, + const Point<2> &p) const { - Assert (i::ExcInvalidIndex(i)); - Assert (false, ExcNotImplemented()); - return Point (); -}; + Assert (i((-4*xi+3) * (1.-eta)*( 2*eta-1), + (1.-xi)*( 2*xi-1) * (-4*eta+3)); + case 1: return Point<2>((-4*xi+1) * (1.-eta)*( 2*eta-1) , + xi *(-2*xi+1) * (-4*eta+3)); + case 2: return Point<2>((-4*xi+1) * eta *(-2*eta+1), + xi *(-2*xi+1) * (-4*eta+1)); + case 3: return Point<2>((-4*xi+3) * eta *(-2*eta+1), + (1.-xi)*( 2*xi-1) * (-4*eta+1)); + case 4: return Point<2>(4 * (1-2*xi) * (1-eta)*(1-2*eta), + 4 * (1-xi)*xi * (4*eta-3)); + case 5: return Point<2>(4 * (4*xi-1) * (1-eta)*eta, + 4 * xi *(-1+2*xi) * (1-2*eta)); + case 6: return Point<2>(4 * (1-2*xi) * eta *(-1+2*eta), + 4 * (1-xi)*xi * (4*eta-1)); + case 7: return Point<2>(4 * (4*xi-3) * (1-eta)*eta, + 4 * (1.-xi)*(1-2*xi) * (1-2*eta)); + case 8: return Point<2>(16 * (1-2*xi) * eta*(1-eta), + 16 * xi*(1-xi) * (1-2*eta)); + }; + return Point<2> (); +}; -template -void FEQuadratic::fill_fe_values (const DoFHandler::cell_iterator &, - const vector > &unit_points, - vector &jacobians, - const bool, - vector > &ansatz_points, - const bool, - vector > &q_points, - const bool, - const 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())); +template <> +void FEQuadraticSub<2>::get_ansatz_points (const typename DoFHandler<2>::cell_iterator &cell, + const Boundary<2>&, + vector > &ansatz_points) const { Assert (ansatz_points.size() == total_dofs, - ExcWrongFieldDimension(ansatz_points.size(), total_dofs)); + ExcWrongFieldDimension (ansatz_points.size(), total_dofs)); + + for (unsigned int vertex=0; vertex<4; ++vertex) + ansatz_points[vertex] = cell->vertex(vertex); - Assert (false, ExcNotImplemented()); + // for the bilinear mapping, the centers + // of the face on the unit cell are mapped + // to the mean coordinates of the vertices + for (unsigned int line=0; line<4; ++line) + ansatz_points[4+line] = (cell->line(line)->vertex(0) + + cell->line(line)->vertex(1)) / 2; + // same for the center of the square: + // since all four linear basis functions + // take on the value 1/4 at the center, + // the center is mapped to the mean + // coordinates of the four vertices + ansatz_points[8] = (ansatz_points[0] + + ansatz_points[1] + + ansatz_points[2] + + ansatz_points[3]) / 4; }; -template -void FEQuadratic::get_ansatz_points (const typename DoFHandler::cell_iterator &, - const Boundary &, - vector > &) const { - Assert (false, ExcNotImplemented()); -}; - - +template <> +void FEQuadraticSub<2>::get_face_ansatz_points (const typename DoFHandler<2>::face_iterator &face, + const Boundary<2> &, + vector > &ansatz_points) const { + Assert (ansatz_points.size() == dofs_per_face, + ExcWrongFieldDimension (ansatz_points.size(), dofs_per_face)); -template -void FEQuadratic::get_face_ansatz_points (const typename DoFHandler::face_iterator &, - const Boundary &, - vector > &) const { - Assert (false, ExcNotImplemented()); + for (unsigned int vertex=0; vertex<2; ++vertex) + ansatz_points[vertex] = face->vertex(vertex); + ansatz_points[2] = (ansatz_points[0] + ansatz_points[1]) / 1; }; -template -void FEQuadratic::get_face_jacobians (const DoFHandler::face_iterator &, - const Boundary &, - const vector > &, - vector &) const { - Assert (false, ExcNotImplemented()); +template <> +void FEQuadraticSub<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 FEQuadratic::get_subface_jacobians (const DoFHandler::face_iterator &face, +template <> +void FEQuadraticSub<2>::get_subface_jacobians (const DoFHandler<2>::face_iterator &face, const unsigned int , - const vector > &, - vector &) const { + 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 ()); - Assert (false, ExcNotImplemented()); + // 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 FEQuadratic::get_normal_vectors (const DoFHandler::cell_iterator &, - const unsigned int, - const Boundary &, - const vector > &, - vector > &) const { - Assert (false, ExcNotImplemented()); +template <> +void FEQuadraticSub<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 FEQuadratic::get_normal_vectors (const DoFHandler::cell_iterator &cell, - const unsigned int face_no, +template <> +void FEQuadraticSub<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell, + const unsigned int face_no, const unsigned int, - const vector > &, - vector > &) const { + 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()))); +}; + +#endif + + + + + +template +void FEQuadraticSub::fill_fe_values (const DoFHandler::cell_iterator &, + const vector > &unit_points, + vector &jacobians, + const bool, + vector > &ansatz_points, + const bool, + vector > &q_points, + const bool, + const 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)); + Assert (false, ExcNotImplemented()); }; + + + + template -void FEQuadratic::get_local_mass_matrix (const DoFHandler::cell_iterator &, +void FEQuadraticSub::get_local_mass_matrix (const DoFHandler::cell_iterator &, const Boundary &, dFMatrix &) const { Assert (false, ExcNotImplemented()); @@ -1102,6 +1273,6 @@ void FECubic::get_local_mass_matrix (const DoFHandler::cell_iterator & // explicit instantiations template class FELinear; -template class FEQuadratic; +template class FEQuadraticSub; template class FECubic; -- 2.39.5