[MG]Line/QuadAccessor classes, but now it is in the
[MG]CellAccessor. Correct this some time.
-
-
-
-
-
-DEAL:
-Is dvector::operator= (double) really useful or does it make more
- confusion than it helps?
-
-
-Give mia::State a much better name, same for Control ! :-]]]
-
-
Let all the reinit functions in /lac free their memory, if reinit
- is called with less requirements. Maybe give free all memory
- if the given dimension is zero. If so, check all deal.II files
+ is called with the given dimension is zero. If so, check all deal.II files
for use of reinit. (At present, new dimension==0 is not allowed.)
[Done for dFMatrix and dVector; still to be done for the other
classes.]
-Why are there all these unsafe casts from VectorBase to dVector
- in dvector.cc? What is VectorBase there for anyway
-
Use unsigned integers for the colnums array in dSMatrixStruct. This
would enhance safety since colnum=-1 would no longer point to a
valid address. How do you mark non-used columns? (gk)
Remove all fe& in vectors.h
-Collection of transfer matrices: how to design?
+No support points for non-Lagrangian elements? Check for Langrange in
+interpolation?
* function value (as opposed to derivatives or the like, as used in the
* Hermitean finite element class or in the quintic Argyris element). It is
* further assumed that a basis function takes its nominal value at a
- * certain point (e.g. linear ansatz functions take their value in the
+ * certain point (e.g. linear trial functions take their value in the
* corners of the element; this restriction rules out spectral elements for
* the present library).
*
* places where this is used in the library come to mind to the author,
* namely the treating of boundary values in the #ProblemBase# class and
* the interpolation in the #VectorTools# collection. You should also
- * look out for other places where explicit use of the ansatz points is
+ * look out for other places where explicit use of the support points is
* made if you want to use elements of other classes. A hint may be the
- * use of the #get_ansatz_points# and #get_face_ansatz_points# functions
+ * use of the #get_support_points# and #get_face_support_points# functions
* of this class.
*
* This also is used in some sense in the
* function of the transformation mapping
* from unit cell to real cell. For
* isoparametric elements, this function
- * is the same as the ansatz functions,
+ * is the same as the trial functions,
* but for sublinear or other mappings,
* they differ.
*/
/**
* Compute the Jacobian matrix and the
- * quadrature points as well as the ansatz
+ * quadrature points as well as the trial
* function locations on the real cell in
* real space from the given cell
* and the given quadrature points on the
const vector<Point<dim> > &unit_points,
vector<dFMatrix> &jacobians,
const bool compute_jacobians,
- vector<Point<dim> > &ansatz_points,
- const bool compute_ansatz_points,
+ vector<Point<dim> > &support_points,
+ const bool compute_support_points,
vector<Point<dim> > &q_points,
const bool compute_q_points,
const dFMatrix &shape_values_transform,
* jacobi matrices as explained in the
* documentation of the #FEValues# class.
*
- * The ansatz points are the
- * off-points of those ansatz functions
+ * The support points are the
+ * off-points of those trial functions
* located on the given face; this
* information is taken over from the
- * #get_face_ansatz_points# function.
+ * #get_face_support_points# function.
*
- * The order of ansatz functions is the
+ * The order of trial functions is the
* same as if it were a cell of dimension
* one less than the present. E.g. in
* two dimensions, the order is first
* specific in this standard implementation:
* The jacobi determinants of the
* transformation from the unit face to the
- * real face, the ansatz points
+ * real face, the support points
* and the outward normal vectors. For
* these fields, there exist pure
* virtual functions, #get_face_jacobians#,
- * #get_face_ansatz_points# and
+ * #get_face_support_points# and
* #get_normal_vectors#.
*
* Though there is a standard
const vector<Point<dim> > &global_unit_points,
vector<dFMatrix> &jacobians,
const bool compute_jacobians,
- vector<Point<dim> > &ansatz_points,
- const bool compute_ansatz_points,
+ vector<Point<dim> > &support_points,
+ const bool compute_support_points,
vector<Point<dim> > &q_points,
const bool compute_q_points,
vector<double> &face_jacobi_determinants,
* of the other parameters is the same as
* for the #fill_fe_face_values# function.
*
- * Since the usage of ansatz points on
+ * Since the usage of support points on
* subfaces is not useful, it is excluded
* from the interface to this function.
*
const Boundary<dim> &boundary) const;
/**
- * Return the ansatz points of the
- * ansatz functions on the unit cell.
+ * Return the support points of the
+ * trial functions on the unit cell.
*
* The function assumes that the
* #unit_points# array already has the
* line. For all other dimensions, an
* overwritten function has to be provided.
*/
- virtual void get_unit_ansatz_points (vector<Point<dim> > &unit_points) const;
+ virtual void get_unit_support_points (vector<Point<dim> > &unit_points) const;
/**
* Compute the off-points of the finite
* abovementioned one.
*
* The function assumes that the
- * #ansatz_points# array already has the
+ * #support_points# array already has the
* right size. The order of points in
* the array matches that returned by
* the #cell->get_dof_indices# function.
* line. For all other dimensions, an
* overwritten function has to be provided.
*/
- virtual void get_ansatz_points (const DoFHandler<dim>::cell_iterator &cell,
+ virtual void get_support_points (const DoFHandler<dim>::cell_iterator &cell,
const Boundary<dim> &boundary,
- vector<Point<dim> > &ansatz_points) const;
+ vector<Point<dim> > &support_points) const;
/**
* Compute the off-points of the finite
* element basis functions located on the
* face. It only returns the off-points
- * of the ansatz functions which are
+ * of the trial functions which are
* located on the face, rather than of
* all basis functions, which is done by
- * the #get_ansatz_points# function.
+ * the #get_support_points# function.
*
* This function produces a subset of
* the information provided by the
* the #fill_fe_face_values# function.
*
* The function assumes that the
- * #ansatz_points# array already has the
+ * #support_points# array already has the
* right size. The order of points in
* the array matches that returned by
* the #face->get_dof_indices# function.
* derived classes should throw an error
* when called with #dim==1#.
*/
- virtual void get_face_ansatz_points (const DoFHandler<dim>::face_iterator &face,
+ virtual void get_face_support_points (const DoFHandler<dim>::face_iterator &face,
const Boundary<dim> &boundary,
- vector<Point<dim> > &ansatz_points) const =0;
+ vector<Point<dim> > &support_points) const =0;
/**
* This is the second separated function
* the computation of the local mass matrix
* is reduced to the computation of a
* weighted evaluation of a polynom in
- * the coordinates of the ansatz points
+ * the coordinates of the support points
* in real space (for linear mappings,
* these are the corner points, for
* quadratic mappings also the center of
dVector &rhs_vector,
const Function<dim> *a = 0);
+ /**
+ * Build Lagrange interpolation
+ matrix of different finite
+ elements.
+ */
+ static void create_interpolation_matrix(const FiniteElement<dim> &high,
+ const FiniteElement<dim> &low,
+ dFMatrix& result);
+
+
/**
* Exception
*/
#if deal_II_dimension == 1
//template <>
-//void FiniteElement<1>::get_ansatz_points (const DoFHandler<1>::cell_iterator &cell,
+//void FiniteElement<1>::get_support_points (const DoFHandler<1>::cell_iterator &cell,
// const Boundary<1> &,
-// vector<Point<1> > &ansatz_points) const;
+// vector<Point<1> > &support_points) const;
template <>
const vector<Point<1> > &unit_points,
vector<dFMatrix> &jacobians,
const bool compute_jacobians,
- vector<Point<1> > &ansatz_points,
- const bool compute_ansatz_points,
+ vector<Point<1> > &support_points,
+ const bool compute_support_points,
vector<Point<1> > &q_points,
const bool compute_q_points,
const dFMatrix &,
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 (support_points.size() == total_dofs,
+ ExcWrongFieldDimension(support_points.size(), total_dofs));
// local mesh width
q_points[i] = cell->vertex(0) + h*unit_points[i];
};
- // compute ansatz points. The first ones
+ // compute support 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)
- get_ansatz_points (cell, boundary, ansatz_points);
+ if (compute_support_points)
+ get_support_points (cell, boundary, support_points);
};
template <>
-void FiniteElement<1>::get_unit_ansatz_points (vector<Point<1> > &ansatz_points) const {
- Assert (ansatz_points.size() == total_dofs,
- ExcWrongFieldDimension(ansatz_points.size(), total_dofs));
- // compute ansatz points. The first ones
+void FiniteElement<1>::get_unit_support_points (vector<Point<1> > &support_points) const {
+ Assert (support_points.size() == total_dofs,
+ ExcWrongFieldDimension(support_points.size(), total_dofs));
+ // compute support points. The first ones
// belong to vertex one, the second ones
// to vertex two, all following are
// equally spaced along the line
unsigned int next = 0;
// first the dofs in the vertices
for (unsigned int i=0; i<dofs_per_vertex; ++i)
- ansatz_points[next++] = Point<1>(0);
+ support_points[next++] = Point<1>(0);
for (unsigned int i=0; i<dofs_per_vertex; ++i)
- ansatz_points[next++] = Point<1>(1);
+ support_points[next++] = Point<1>(1);
// now dofs on line
for (unsigned int i=0; i<dofs_per_line; ++i)
- ansatz_points[next++] = Point<1>((i+1.0)/(dofs_per_line+1.0));
+ support_points[next++] = Point<1>((i+1.0)/(dofs_per_line+1.0));
};
template <>
-void FiniteElement<1>::get_ansatz_points (const DoFHandler<1>::cell_iterator &cell,
+void FiniteElement<1>::get_support_points (const DoFHandler<1>::cell_iterator &cell,
const Boundary<1> &,
- vector<Point<1> > &ansatz_points) const {
- Assert (ansatz_points.size() == total_dofs,
- ExcWrongFieldDimension(ansatz_points.size(), total_dofs));
- // compute ansatz points. The first ones
+ vector<Point<1> > &support_points) const {
+ Assert (support_points.size() == total_dofs,
+ ExcWrongFieldDimension(support_points.size(), total_dofs));
+ // compute support points. The first ones
// belong to vertex one, the second ones
// to vertex two, all following are
// equally spaced along the line
// 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[next++] = cell->vertex(vertex);
+ support_points[next++] = cell->vertex(vertex);
// now dofs on line
for (unsigned int i=0; i<dofs_per_line; ++i)
- ansatz_points[next++] = cell->vertex(0) +
+ support_points[next++] = cell->vertex(0) +
Point<1>((i+1.0)/(dofs_per_line+1.0)*h);
};
const vector<Point<dim> > &global_unit_points,
vector<dFMatrix> &jacobians,
const bool compute_jacobians,
- vector<Point<dim> > &ansatz_points,
- const bool compute_ansatz_points,
+ vector<Point<dim> > &support_points,
+ const bool compute_support_points,
vector<Point<dim> > &q_points,
const bool compute_q_points,
vector<double> &face_jacobi_determinants,
ExcWrongFieldDimension(q_points.size(), unit_points.size()));
Assert (global_unit_points.size() == unit_points.size(),
ExcWrongFieldDimension(global_unit_points.size(), unit_points.size()));
- Assert (ansatz_points.size() == dofs_per_face,
- ExcWrongFieldDimension(ansatz_points.size(), dofs_per_face));
+ Assert (support_points.size() == dofs_per_face,
+ ExcWrongFieldDimension(support_points.size(), dofs_per_face));
vector<Point<dim> > dummy(total_dofs);
fill_fe_values (cell, global_unit_points,
shape_values_transform, shape_gradients_transform,
boundary);
- if (compute_ansatz_points)
- get_face_ansatz_points (cell->face(face_no), boundary, ansatz_points);
+ if (compute_support_points)
+ get_face_support_points (cell->face(face_no), boundary, support_points);
if (compute_face_jacobians)
get_face_jacobians (cell->face(face_no), boundary,
template <int dim>
-void FiniteElement<dim>::get_unit_ansatz_points (vector<Point<dim> > &) const {
+void FiniteElement<dim>::get_unit_support_points (vector<Point<dim> > &) const {
Assert (false, ExcPureFunctionCalled());
};
template <int dim>
-void FiniteElement<dim>::get_ansatz_points (const DoFHandler<dim>::cell_iterator &,
+void FiniteElement<dim>::get_support_points (const DoFHandler<dim>::cell_iterator &,
const Boundary<dim> &,
vector<Point<dim> > &) const {
Assert (false, ExcPureFunctionCalled());
const vector<Point<1> > &unit_points,
vector<dFMatrix> &jacobians,
const bool compute_jacobians,
- vector<Point<1> > &ansatz_points,
- const bool compute_ansatz_points,
+ vector<Point<1> > &support_points,
+ const bool compute_support_points,
vector<Point<1> > &q_points,
const bool compute_q_points,
const dFMatrix &shape_values_transform,
// simply pass down
FiniteElement<1>::fill_fe_values (cell, unit_points,
jacobians, compute_jacobians,
- ansatz_points, compute_ansatz_points,
+ support_points, compute_support_points,
q_points, compute_q_points,
shape_values_transform, shape_gradients_transform,
boundary);
const vector<Point<dim> > &unit_points,
vector<dFMatrix> &jacobians,
const bool compute_jacobians,
- vector<Point<dim> > &ansatz_points,
- const bool compute_ansatz_points,
+ vector<Point<dim> > &support_points,
+ const bool compute_support_points,
vector<Point<dim> > &q_points,
const bool compute_q_points,
const dFMatrix &shape_values_transform,
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 (support_points.size() == total_dofs,
+ ExcWrongFieldDimension(support_points.size(), total_dofs));
unsigned int n_points=unit_points.size();
};
};
- if (compute_ansatz_points)
- get_ansatz_points (cell, boundary, ansatz_points);
+ if (compute_support_points)
+ get_support_points (cell, boundary, support_points);
};
rhs(i) += values(i,point) *
rhs_values[point] *
weights[point];
-};
-
-
+}
+template<int dim> void
+MatrixCreator<dim>::create_interpolation_matrix(const FiniteElement<dim> &high,
+ const FiniteElement<dim> &low,
+ dFMatrix& result)
+{
+ result.reinit (high.total_dofs, low.total_dofs);
+ vector<Point<dim> > unit_support_points (high.total_dofs);
+ high.get_unit_support_points (unit_support_points);
+
+ for (unsigned int i=0; i<high.total_dofs; ++i)
+ for (unsigned int j=0; j<low.total_dofs; ++j)
+ result(i,j) = low.shape_value (j, unit_support_points[i]);
+}
endc = dof.end();
vector<int> dofs_on_cell (fe.total_dofs);
vector<double> dof_values_on_cell (fe.total_dofs);
- vector<Point<dim> > ansatz_points (fe.total_dofs);
+ vector<Point<dim> > support_points (fe.total_dofs);
for (; cell!=endc; ++cell)
{
// for each cell:
// get location of finite element
// off-points
- fe.get_ansatz_points (cell, boundary, ansatz_points);
+ fe.get_support_points (cell, boundary, support_points);
// get function values at these points
- function.value_list (ansatz_points, dof_values_on_cell);
+ function.value_list (support_points, dof_values_on_cell);
// get indices of the dofs on this cell
cell->get_dof_indices (dofs_on_cell);
// distribute function values to the
// boundary values of dofs on this
// face
face->get_dof_indices (face_dofs);
- fe.get_face_ansatz_points (face, boundary, dof_locations);
+ fe.get_face_support_points (face, boundary, dof_locations);
function_ptr->second->value_list (dof_locations, dof_values);
// enter into list
// \psi(x_j)=\sum_i v_i \phi_i(x_j)
// with v_i the nodal values of the
// fe_function and \phi_i(x_j) the
- // matrix of the ansatz function
+ // matrix of the trial function
// values at the integration point
// x_j. Then the vector
// of the \psi(x_j) is v*Phi with