]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Support getting the position of an ansatz function in real space on the real cell...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 1 Apr 1998 15:52:50 +0000 (15:52 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 1 Apr 1998 15:52:50 +0000 (15:52 +0000)
git-svn-id: https://svn.dealii.org/trunk@109 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe.h
deal.II/deal.II/include/fe/fe_lib.lagrange.h
deal.II/deal.II/source/fe/fe.cc
deal.II/deal.II/source/fe/fe_lib.linear.cc
deal.II/deal.II/source/numerics/base.cc

index f3ffc7d85cdbb47daf014c4b8a61b74c501bbbe5..6978d7734d6c343155f021e26cabf128e0053c0f 100644 (file)
@@ -178,6 +178,31 @@ class FEValues {
                                      * quadrature points.
                                      */
     const vector<Point<dim> > & 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<dim> & ansatz_point (const unsigned int i) const;
+
+                                    /**
+                                     * Return a pointer to the vector of points
+                                     * denoting the location of the ansatz
+                                     * functions.
+                                     */
+    const vector<Point<dim> > & get_ansatz_points () const;
     
                                     /**
                                      * Return the Jacobi determinant times
@@ -288,7 +313,16 @@ class FEValues {
                                      * points on the reference element.
                                      */
     vector<Point<dim> >  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<Point<dim> >  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<Point<dim> >               &unit_points,
                                 vector<dFMatrix>    &jacobians,
                                 const bool           compute_jacobians,
-                                vector<Point<dim> > &points,
-                                const bool           compute_points) const;
+                                vector<Point<dim> > &ansatz_points,
+                                const bool           compute_ansatz_points,
+                                vector<Point<dim> > &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<Point<1> >               &unit_points,
                                 vector<dFMatrix>  &jacobians,
                                 const bool         compute_jacobians,
-                                vector<Point<1> > &points,
-                                const bool         compute_points) const;
+                                vector<Point<1> > &ansatz_points,
+                                const bool         compute_ansatz_points,
+                                vector<Point<1> > &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<Point<2> >               &unit_points,
                                 vector<dFMatrix>  &jacobians,
                                 const bool         compute_jacobians,
-                                vector<Point<2> > &points,
-                                const bool         compute_points) const;
+                                vector<Point<2> > &ansatz_points,
+                                const bool         compute_ansatz_points,
+                                vector<Point<2> > &q_points,
+                                const bool         compute_q_points) const;
 };
 
 
@@ -796,6 +844,16 @@ FEValues<dim>::get_quadrature_points () const {
 
 
 
+template <int dim>
+inline
+const vector<Point<dim> > &
+FEValues<dim>::get_ansatz_points () const {
+  Assert (update_flags.update_ansatz_points, ExcAccessToUninitializedField());
+  return ansatz_points;
+};
+
+
+
 template <int dim>
 inline
 const vector<double> &
index 279f872740ad03f8e343380bd81443ba7e4bb08b..c6171f526864fc3af4b9932996a135e55910ab48 100644 (file)
@@ -46,7 +46,9 @@ class FELinear : public FiniteElement<dim> {
 
                                     /**
                                      * 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<dim> {
                                 const vector<Point<dim> >               &unit_points,
                                 vector<dFMatrix>    &jacobians,
                                 const bool           compute_jacobians,
-                                vector<Point<dim> > &points,
-                                const bool           compute_points) const;
+                                vector<Point<dim> > &ansatz_points,
+                                const bool           compute_ansatz_points,
+                                vector<Point<dim> > &q_points,
+                                const bool           compute_q_points) const;
 };
 
 
@@ -107,7 +111,9 @@ class FEQuadratic : public FiniteElement<dim> {
 
                                     /**
                                      * 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<dim> {
                                 const vector<Point<dim> >               &unit_points,
                                 vector<dFMatrix>    &jacobians,
                                 const bool           compute_jacobians,
-                                vector<Point<dim> > &points,
-                                const bool           compute_points) const;
+                                vector<Point<dim> > &ansatz_points,
+                                const bool           compute_ansatz_points,
+                                vector<Point<dim> > &q_points,
+                                const bool           compute_q_points) const;
 };
 
 
@@ -168,7 +176,9 @@ class FECubic : public FiniteElement<dim> {
 
                                     /**
                                      * 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<dim> {
                                 const vector<Point<dim> >               &unit_points,
                                 vector<dFMatrix>    &jacobians,
                                 const bool           compute_jacobians,
-                                vector<Point<dim> > &points,
-                                const bool           compute_points) const;
+                                vector<Point<dim> > &ansatz_points,
+                                const bool           compute_ansatz_points,
+                                vector<Point<dim> > &q_points,
+                                const bool           compute_q_points) const;
 };
 
 
index 8e0895ebc073db3179313269d7db00371545bee3..3a12c2a61cf1a6dacca4bbfd98044892ccc5ce1c 100644 (file)
@@ -92,6 +92,16 @@ const Point<dim> & FEValues<dim>::quadrature_point (const unsigned int i) const
 
 
 
+template <int dim>
+const Point<dim> & FEValues<dim>::ansatz_point (const unsigned int i) const {
+  Assert (i<ansatz_points.size(), ExcInvalidIndex(i, ansatz_points.size()));
+  Assert (update_flags.update_ansatz_points, ExcAccessToUninitializedField());
+  
+  return ansatz_points[i];
+};
+
+
+
 template <int dim>
 double FEValues<dim>::JxW (const unsigned int i) const {
   Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
@@ -112,6 +122,8 @@ void FEValues<dim>::reinit (const typename Triangulation<dim>::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<dim>::fill_fe_values (const typename Triangulation<dim>::
                                             vector<dFMatrix> &,
                                             const bool,
                                             vector<Point<dim> > &,
+                                            const bool,
+                                            vector<Point<dim> > &,
                                             const bool) const {
   Assert (false, ExcPureFunctionCalled());
 };
@@ -254,17 +268,39 @@ void FiniteElement<1>::fill_fe_values (const Triangulation<1>::cell_iterator &ce
                                       const vector<Point<1> > &unit_points,
                                       vector<dFMatrix>  &jacobians,
                                       const bool         compute_jacobians,
-                                      vector<Point<1> > &points,
-                                      const bool         compute_points) const {
+                                      vector<Point<1> > &ansatz_points,
+                                      const bool         compute_ansatz_points,
+                                      vector<Point<1> > &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; i<points.size(); ++i) 
+  for (unsigned int i=0; i<q_points.size(); ++i) 
     {
       if (compute_jacobians)
        jacobians[i](0,0) = 1./h;
-      if (compute_points)
-       points[i] = cell->vertex(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; i<dofs_per_vertex; ++i)
+         ansatz_points.push_back (cell->vertex(vertex));
+      
+                                      // now dofs on line
+      for (unsigned int i=0; i<dofs_per_line; ++i) 
+       ansatz_points.push_back (cell->vertex(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<dFMatrix> &,
                                       const bool,
                                       vector<Point<2> > &,
+                                      const bool,
+                                      vector<Point<2> > &,
                                       const bool) const {
   Assert (false, ExcPureFunctionCalled());
 };
index a9c6fd16a27668f8f97e652336b490e14e9fc008..93cb12965e89915a5bfc17d22d7b1bb50ca3327f 100644 (file)
@@ -80,12 +80,15 @@ void FELinear<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell,
                                  const vector<Point<1> >               &unit_points,
                                  vector<dFMatrix>  &jacobians,
                                  const bool         compute_jacobians,
-                                 vector<Point<1> > &points,
-                                 const bool         compute_points) const {
+                                 vector<Point<1> > &ansatz_points,
+                                 const bool         compute_ansatz_points,
+                                 vector<Point<1> > &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<Point<2> >               &unit_points,
                                  vector<dFMatrix>  &jacobians,
                                  const bool         compute_jacobians,
-                                 vector<Point<2> > &points,
-                                 const bool         compute_points) const {
+                                 vector<Point<2> > &ansatz_points,
+                                 const bool         compute_ansatz_points,
+                                 vector<Point<2> > &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<n_points; ++i)
-       points[i] = Point<dim> ();
+       q_points[i] = Point<dim> ();
       
                                       // 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<n_vertices; ++j) 
        for (unsigned int l=0; l<n_points; ++l) 
-         points[l] += vertices[j] * shape_value(j, unit_points[l]);
+         q_points[l] += vertices[j] * shape_value(j, unit_points[l]);
     };
   
 
@@ -303,6 +308,12 @@ void FELinear<2>::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<Point<1> >               &unit_points,
                                     vector<dFMatrix>  &jacobians,
                                     const bool         compute_jacobians,
-                                    vector<Point<1> > &points,
-                                    const bool         compute_points) const {
+                                    vector<Point<1> > &ansatz_points,
+                                    const bool         compute_ansatz_points,
+                                    vector<Point<1> > &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<dFMatrix>  &,
                                     const bool,
                                     vector<Point<2> > &,
+                                    const bool,
+                                    vector<Point<2> > &,
                                     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<Point<1> >               &unit_points,
                                 vector<dFMatrix>  &jacobians,
                                 const bool         compute_jacobians,
-                                vector<Point<1> > &points,
-                                const bool          compute_points) const {
+                                vector<Point<1> > &ansatz_points,
+                                const bool         compute_ansatz_points,
+                                vector<Point<1> > &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<dFMatrix>  &,
                                 const bool,
                                 vector<Point<2> > &,
+                                const bool,
+                                vector<Point<2> > &,
                                 const bool) const {
   Assert (false, typename FiniteElementBase<2>::ExcNotImplemented());
 };
index 64bc08bfc7404e4ad6a8eeba218cd26007e9a3d2..1627c1d3566bce78978966d237c58b42539b1a63 100644 (file)
@@ -400,7 +400,10 @@ ProblemBase<dim>::make_boundary_value_list (const DirichletBC &dirichlet_bc,
        vector<double>      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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.