From 0eb89d842a83e35cba15f2099a1ecf3845bfd430 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 1 Apr 1998 15:52:50 +0000 Subject: [PATCH] Support getting the position of an ansatz function in real space on the real cell rather than on the unit cell. git-svn-id: https://svn.dealii.org/trunk@109 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe.h | 84 +++++++++++++++++--- deal.II/deal.II/include/fe/fe_lib.lagrange.h | 30 ++++--- deal.II/deal.II/source/fe/fe.cc | 48 +++++++++-- deal.II/deal.II/source/fe/fe_lib.linear.cc | 49 ++++++++---- deal.II/deal.II/source/numerics/base.cc | 3 + 5 files changed, 173 insertions(+), 41 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index f3ffc7d85c..6978d7734d 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -178,6 +178,31 @@ class FEValues { * quadrature points. */ const vector > & get_quadrature_points () const; + + /** + * Return the point in real space where + * the #i#th ansatz function is located + * (location is in the sense of where it + * assumes its nominal properties, e.g. at + * the vertex of a cell, at the center of + * a line, etc). + * + * This function is needed for the + * interpolation problem: if we want to + * transfer a continuous function to a + * finite element function by interpolation + * we have to take the continuous + * function's value at the ansatz function + * locations. + */ + const Point & ansatz_point (const unsigned int i) const; + + /** + * Return a pointer to the vector of points + * denoting the location of the ansatz + * functions. + */ + const vector > & get_ansatz_points () const; /** * Return the Jacobi determinant times @@ -288,7 +313,16 @@ class FEValues { * points on the reference element. */ vector > unit_quadrature_points; - + + /** + * Array of points denoting the off-point + * of the ansatz functions. In real space + * (no-one seems to need the off-point + * on the unit cell, so no function is + * provided for this). + */ + vector > ansatz_points; + /** * Store the jacobi matrices at the * different quadrature points. This field @@ -410,7 +444,9 @@ class FiniteElementBase { /** * Compute the Jacobian matrix and the - * quadrature points from the given cell + * quadrature points as well as the ansatz + * function locations on the real cell in + * real space from the given cell * and the given quadrature points on the * unit cell. The Jacobian matrix is to * be computed at every quadrature point. @@ -419,11 +455,11 @@ class FiniteElementBase { * elements need different transformations * of the unit cell to a real cell. * - * The computation of the two fields may + * The computation of the three fields may * share some common code, which is why we * put it in one function. However, it may * not always be necessary to really - * compute both fields, so there are two + * compute all fields, so there are * bool flags which tell the function which * of the fields to actually compute. * @@ -441,8 +477,10 @@ class FiniteElementBase { const vector > &unit_points, vector &jacobians, const bool compute_jacobians, - vector > &points, - const bool compute_points) const; + vector > &ansatz_points, + const bool compute_ansatz_points, + vector > &q_points, + const bool compute_q_points) const; /** * Comparison operator. We also check for @@ -609,7 +647,9 @@ class FiniteElement<1> : public FiniteElementBase<1> { /** * Compute the Jacobian matrix and the - * quadrature points from the given cell + * quadrature points as well as the ansatz + * function locations on the real cell in + * real space from the given cell * and the given quadrature points on the * unit cell. The Jacobian matrix is to * be computed at every quadrature point. @@ -627,14 +667,18 @@ class FiniteElement<1> : public FiniteElementBase<1> { * of higher than first order with * non-equidistant integration points, e.g. * with an exponential dependence from the - * distance to the origin. + * distance to the origin. The standard + * implementation distributes the dofs on + * the line equidistantly. */ virtual void fill_fe_values (const Triangulation<1>::cell_iterator &cell, const vector > &unit_points, vector &jacobians, const bool compute_jacobians, - vector > &points, - const bool compute_points) const; + vector > &ansatz_points, + const bool compute_ansatz_points, + vector > &q_points, + const bool compute_q_points) const; }; @@ -733,7 +777,9 @@ class FiniteElement<2> : public FiniteElementBase<2> { /** * Compute the Jacobian matrix and the - * quadrature points from the given cell + * quadrature points as well as the ansatz + * function locations on the real cell in + * real space from the given cell * and the given quadrature points on the * unit cell. The Jacobian matrix is to * be computed at every quadrature point. @@ -755,8 +801,10 @@ class FiniteElement<2> : public FiniteElementBase<2> { const vector > &unit_points, vector &jacobians, const bool compute_jacobians, - vector > &points, - const bool compute_points) const; + vector > &ansatz_points, + const bool compute_ansatz_points, + vector > &q_points, + const bool compute_q_points) const; }; @@ -796,6 +844,16 @@ FEValues::get_quadrature_points () const { +template +inline +const vector > & +FEValues::get_ansatz_points () const { + Assert (update_flags.update_ansatz_points, ExcAccessToUninitializedField()); + return ansatz_points; +}; + + + template inline const vector & 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 279f872740..c6171f5268 100644 --- a/deal.II/deal.II/include/fe/fe_lib.lagrange.h +++ b/deal.II/deal.II/include/fe/fe_lib.lagrange.h @@ -46,7 +46,9 @@ class FELinear : public FiniteElement { /** * Compute the Jacobian matrix and the - * quadrature points from the given cell + * quadrature points as well as the ansatz + * function locations on the real cell in + * real space from the given cell * and the given quadrature points on the * unit cell. The Jacobian matrix is to * be computed at every quadrature point. @@ -71,8 +73,10 @@ class FELinear : public FiniteElement { const vector > &unit_points, vector &jacobians, const bool compute_jacobians, - vector > &points, - const bool compute_points) const; + vector > &ansatz_points, + const bool compute_ansatz_points, + vector > &q_points, + const bool compute_q_points) const; }; @@ -107,7 +111,9 @@ class FEQuadratic : public FiniteElement { /** * Compute the Jacobian matrix and the - * quadrature points from the given cell + * quadrature points as well as the ansatz + * function locations on the real cell in + * real space from the given cell * and the given quadrature points on the * unit cell. The Jacobian matrix is to * be computed at every quadrature point. @@ -132,8 +138,10 @@ class FEQuadratic : public FiniteElement { const vector > &unit_points, vector &jacobians, const bool compute_jacobians, - vector > &points, - const bool compute_points) const; + vector > &ansatz_points, + const bool compute_ansatz_points, + vector > &q_points, + const bool compute_q_points) const; }; @@ -168,7 +176,9 @@ class FECubic : public FiniteElement { /** * Compute the Jacobian matrix and the - * quadrature points from the given cell + * quadrature points as well as the ansatz + * function locations on the real cell in + * real space from the given cell * and the given quadrature points on the * unit cell. The Jacobian matrix is to * be computed at every quadrature point. @@ -193,8 +203,10 @@ class FECubic : public FiniteElement { const vector > &unit_points, vector &jacobians, const bool compute_jacobians, - vector > &points, - const bool compute_points) const; + vector > &ansatz_points, + const bool compute_ansatz_points, + vector > &q_points, + const bool compute_q_points) const; }; diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index 8e0895ebc0..3a12c2a61c 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -92,6 +92,16 @@ const Point & FEValues::quadrature_point (const unsigned int i) const +template +const Point & FEValues::ansatz_point (const unsigned int i) const { + Assert (i double FEValues::JxW (const unsigned int i) const { Assert (i::reinit (const typename Triangulation::cell_iterator &ce unit_quadrature_points, jacobi_matrices, update_flags.update_jacobians, + ansatz_points, + update_flags.update_ansatz_points, quadrature_points, update_flags.update_q_points); @@ -233,6 +245,8 @@ void FiniteElementBase::fill_fe_values (const typename Triangulation:: vector &, const bool, vector > &, + const bool, + vector > &, const bool) const { Assert (false, ExcPureFunctionCalled()); }; @@ -254,17 +268,39 @@ void FiniteElement<1>::fill_fe_values (const Triangulation<1>::cell_iterator &ce const vector > &unit_points, vector &jacobians, const bool compute_jacobians, - vector > &points, - const bool compute_points) const { + vector > &ansatz_points, + const bool compute_ansatz_points, + vector > &q_points, + const bool compute_q_points) const { // local mesh width const double h=(cell->vertex(1)(0) - cell->vertex(0)(0)); - for (unsigned int i=0; ivertex(0) + h*unit_points[i]; + if (compute_q_points) + q_points[i] = cell->vertex(0) + h*unit_points[i]; + }; + + // compute ansatz points. The first ones + // belong to vertex one, the second ones + // to vertex two, all following are + // equally spaced along the line + if (compute_ansatz_points) + { + ansatz_points.erase (ansatz_points.begin(), ansatz_points.end()); + ansatz_points.reserve (total_dofs); + + // first the dofs in the vertices + for (unsigned int vertex=0; vertex<2; vertex++) + for (unsigned int i=0; ivertex(vertex)); + + // now dofs on line + for (unsigned int i=0; ivertex(0) + + Point<1>((i+1.0)/(total_dofs+1.0)*h)); }; }; @@ -284,6 +320,8 @@ void FiniteElement<2>::fill_fe_values (const Triangulation<2>::cell_iterator &, vector &, const bool, vector > &, + const bool, + vector > &, const bool) const { Assert (false, ExcPureFunctionCalled()); }; 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 a9c6fd16a2..93cb12965e 100644 --- a/deal.II/deal.II/source/fe/fe_lib.linear.cc +++ b/deal.II/deal.II/source/fe/fe_lib.linear.cc @@ -80,12 +80,15 @@ void FELinear<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell, const vector > &unit_points, vector &jacobians, const bool compute_jacobians, - vector > &points, - const bool compute_points) const { + vector > &ansatz_points, + const bool compute_ansatz_points, + vector > &q_points, + const bool compute_q_points) const { // simply pass down FiniteElement<1>::fill_fe_values (cell, unit_points, jacobians, compute_jacobians, - points, compute_points); + ansatz_points, compute_ansatz_points, + q_points, compute_q_points); }; @@ -232,8 +235,10 @@ void FELinear<2>::fill_fe_values (const Triangulation<2>::cell_iterator &cell, const vector > &unit_points, vector &jacobians, const bool compute_jacobians, - vector > &points, - const bool compute_points) const { + vector > &ansatz_points, + const bool compute_ansatz_points, + vector > &q_points, + const bool compute_q_points) const { const unsigned int dim=2; const unsigned int n_vertices=4; unsigned int n_points=unit_points.size(); @@ -243,11 +248,11 @@ void FELinear<2>::fill_fe_values (const Triangulation<2>::cell_iterator &cell, vertices[l] = cell->vertex(l); - if (compute_points) + if (compute_q_points) { // initialize points to zero for (unsigned int i=0; i (); + q_points[i] = Point (); // note: let x_l be the vector of the // lth quadrature point in real space and @@ -259,7 +264,7 @@ void FELinear<2>::fill_fe_values (const Triangulation<2>::cell_iterator &cell, // x_l(xi_l) = sum_j p_j N_j(xi_l) for (unsigned int j=0; j::fill_fe_values (const Triangulation<2>::cell_iterator &cell, jacobians[l].invert(M); }; }; + + // compute ansatz points, which are + // the corners for linear elements + if (compute_ansatz_points) + for (unsigned int vertex=0; vertex<4; ++vertex) + ansatz_points[vertex] = vertices[vertex]; }; @@ -319,12 +330,15 @@ void FEQuadratic<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell const vector > &unit_points, vector &jacobians, const bool compute_jacobians, - vector > &points, - const bool compute_points) const { + vector > &ansatz_points, + const bool compute_ansatz_points, + vector > &q_points, + const bool compute_q_points) const { // simply pass down FiniteElement<1>::fill_fe_values (cell, unit_points, jacobians, compute_jacobians, - points, compute_points); + ansatz_points, compute_ansatz_points, + q_points, compute_q_points); }; @@ -373,6 +387,8 @@ void FEQuadratic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &, vector &, const bool, vector > &, + const bool, + vector > &, const bool) const { Assert (false, typename FiniteElementBase<2>::ExcNotImplemented()); }; @@ -391,12 +407,15 @@ void FECubic<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell, const vector > &unit_points, vector &jacobians, const bool compute_jacobians, - vector > &points, - const bool compute_points) const { + vector > &ansatz_points, + const bool compute_ansatz_points, + vector > &q_points, + const bool compute_q_points) const { // simply pass down FiniteElement<1>::fill_fe_values (cell, unit_points, jacobians, compute_jacobians, - points, compute_points); + ansatz_points, compute_ansatz_points, + q_points, compute_q_points); }; @@ -437,6 +456,8 @@ void FECubic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &, vector &, const bool, vector > &, + const bool, + vector > &, const bool) const { Assert (false, typename FiniteElementBase<2>::ExcNotImplemented()); }; diff --git a/deal.II/deal.II/source/numerics/base.cc b/deal.II/deal.II/source/numerics/base.cc index 64bc08bfc7..1627c1d356 100644 --- a/deal.II/deal.II/source/numerics/base.cc +++ b/deal.II/deal.II/source/numerics/base.cc @@ -400,7 +400,10 @@ ProblemBase::make_boundary_value_list (const DirichletBC &dirichlet_bc, vector dof_values; face->get_dof_indices (face_dofs); + //physical location + Assert (false, ExcNotImplemented()); + (*function_ptr).second->value_list (dof_locations, dof_values); // enter into list -- 2.39.5