* 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
* 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
/**
* 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.
* 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.
*
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
/**
* 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.
* 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;
};
/**
* 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.
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;
};
+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> &
/**
* 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.
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;
};
/**
* 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.
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;
};
/**
* 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.
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;
};
+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));
unit_quadrature_points,
jacobi_matrices,
update_flags.update_jacobians,
+ ansatz_points,
+ update_flags.update_ansatz_points,
quadrature_points,
update_flags.update_q_points);
vector<dFMatrix> &,
const bool,
vector<Point<dim> > &,
+ const bool,
+ vector<Point<dim> > &,
const bool) const {
Assert (false, ExcPureFunctionCalled());
};
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));
};
};
vector<dFMatrix> &,
const bool,
vector<Point<2> > &,
+ const bool,
+ vector<Point<2> > &,
const bool) const {
Assert (false, ExcPureFunctionCalled());
};
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);
};
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();
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
// 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]);
};
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];
};
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);
};
vector<dFMatrix> &,
const bool,
vector<Point<2> > &,
+ const bool,
+ vector<Point<2> > &,
const bool) const {
Assert (false, typename FiniteElementBase<2>::ExcNotImplemented());
};
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);
};
vector<dFMatrix> &,
const bool,
vector<Point<2> > &,
+ const bool,
+ vector<Point<2> > &,
const bool) const {
Assert (false, typename FiniteElementBase<2>::ExcNotImplemented());
};
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