From 0d9504dbae78a4042e4eef674634480a5f977305 Mon Sep 17 00:00:00 2001 From: guido Date: Fri, 4 Sep 1998 15:05:30 +0000 Subject: [PATCH] ansatz replaced by support git-svn-id: https://svn.dealii.org/trunk@572 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/Todo | 22 +------ deal.II/deal.II/include/fe/fe.h | 56 ++++++++--------- deal.II/deal.II/include/numerics/matrices.h | 10 +++ deal.II/deal.II/source/fe/fe.cc | 62 +++++++++---------- .../deal.II/source/fe/fe_linear_mapping.cc | 18 +++--- deal.II/deal.II/source/numerics/matrices.cc | 17 ++++- deal.II/deal.II/source/numerics/vectors.cc | 10 +-- 7 files changed, 100 insertions(+), 95 deletions(-) diff --git a/deal.II/deal.II/Todo b/deal.II/deal.II/Todo index 38ac430e50..20797c0285 100644 --- a/deal.II/deal.II/Todo +++ b/deal.II/deal.II/Todo @@ -137,29 +137,12 @@ I suppose the [mg_]get_dof_values really belongs to the [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) @@ -171,4 +154,5 @@ FeValues: add flexibility for update flags 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? diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index cbaa74528d..29d3f03701 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -445,7 +445,7 @@ struct FiniteElementBase : * 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). * @@ -454,9 +454,9 @@ struct FiniteElementBase : * 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 @@ -533,7 +533,7 @@ class FiniteElement : public FiniteElementBase { * 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. */ @@ -550,7 +550,7 @@ class FiniteElement : public FiniteElementBase { /** * 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 @@ -608,8 +608,8 @@ class FiniteElement : public FiniteElementBase { const vector > &unit_points, vector &jacobians, const bool compute_jacobians, - vector > &ansatz_points, - const bool compute_ansatz_points, + vector > &support_points, + const bool compute_support_points, vector > &q_points, const bool compute_q_points, const dFMatrix &shape_values_transform, @@ -658,13 +658,13 @@ class FiniteElement : public FiniteElementBase { * 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 @@ -694,11 +694,11 @@ class FiniteElement : public FiniteElementBase { * 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 @@ -726,8 +726,8 @@ class FiniteElement : public FiniteElementBase { const vector > &global_unit_points, vector &jacobians, const bool compute_jacobians, - vector > &ansatz_points, - const bool compute_ansatz_points, + vector > &support_points, + const bool compute_support_points, vector > &q_points, const bool compute_q_points, vector &face_jacobi_determinants, @@ -749,7 +749,7 @@ class FiniteElement : public FiniteElementBase { * 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. * @@ -784,8 +784,8 @@ class FiniteElement : public FiniteElementBase { const Boundary &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 @@ -799,7 +799,7 @@ class FiniteElement : public FiniteElementBase { * line. For all other dimensions, an * overwritten function has to be provided. */ - virtual void get_unit_ansatz_points (vector > &unit_points) const; + virtual void get_unit_support_points (vector > &unit_points) const; /** * Compute the off-points of the finite @@ -818,7 +818,7 @@ class FiniteElement : public FiniteElementBase { * 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. @@ -829,18 +829,18 @@ class FiniteElement : public FiniteElementBase { * line. For all other dimensions, an * overwritten function has to be provided. */ - virtual void get_ansatz_points (const DoFHandler::cell_iterator &cell, + virtual void get_support_points (const DoFHandler::cell_iterator &cell, const Boundary &boundary, - vector > &ansatz_points) const; + vector > &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 @@ -863,7 +863,7 @@ class FiniteElement : public FiniteElementBase { * 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. @@ -873,9 +873,9 @@ class FiniteElement : public FiniteElementBase { * derived classes should throw an error * when called with #dim==1#. */ - virtual void get_face_ansatz_points (const DoFHandler::face_iterator &face, + virtual void get_face_support_points (const DoFHandler::face_iterator &face, const Boundary &boundary, - vector > &ansatz_points) const =0; + vector > &support_points) const =0; /** * This is the second separated function @@ -980,7 +980,7 @@ class FiniteElement : public FiniteElementBase { * 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 diff --git a/deal.II/deal.II/include/numerics/matrices.h b/deal.II/deal.II/include/numerics/matrices.h index 2a6bcaa4f7..b417f8179b 100644 --- a/deal.II/deal.II/include/numerics/matrices.h +++ b/deal.II/deal.II/include/numerics/matrices.h @@ -286,6 +286,16 @@ class MatrixCreator { dVector &rhs_vector, const Function *a = 0); + /** + * Build Lagrange interpolation + matrix of different finite + elements. + */ + static void create_interpolation_matrix(const FiniteElement &high, + const FiniteElement &low, + dFMatrix& result); + + /** * Exception */ diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index 7ec449c73b..4160b83abb 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -169,9 +169,9 @@ bool FiniteElementBase::operator == (const FiniteElementBase &f) const #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 > &ansatz_points) const; +// vector > &support_points) const; template <> @@ -179,8 +179,8 @@ void FiniteElement<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell, const vector > &unit_points, vector &jacobians, const bool compute_jacobians, - vector > &ansatz_points, - const bool compute_ansatz_points, + vector > &support_points, + const bool compute_support_points, vector > &q_points, const bool compute_q_points, const dFMatrix &, @@ -190,8 +190,8 @@ void FiniteElement<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell, 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 @@ -208,12 +208,12 @@ void FiniteElement<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell, 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); }; @@ -263,34 +263,34 @@ void FiniteElement<1>::fill_fe_subface_values (const DoFHandler<1>::cell_iterato template <> -void FiniteElement<1>::get_unit_ansatz_points (vector > &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 > &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(0); + support_points[next++] = Point<1>(0); for (unsigned int i=0; i(1); + support_points[next++] = Point<1>(1); // now dofs on line for (unsigned int i=0; i((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 > &ansatz_points) const { - Assert (ansatz_points.size() == total_dofs, - ExcWrongFieldDimension(ansatz_points.size(), total_dofs)); - // compute ansatz points. The first ones + vector > &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 @@ -300,11 +300,11 @@ void FiniteElement<1>::get_ansatz_points (const DoFHandler<1>::cell_iterator &ce // first the dofs in the vertices for (unsigned int vertex=0; vertex<2; vertex++) for (unsigned int i=0; ivertex(vertex); + support_points[next++] = cell->vertex(vertex); // now dofs on line for (unsigned int i=0; ivertex(0) + + support_points[next++] = cell->vertex(0) + Point<1>((i+1.0)/(dofs_per_line+1.0)*h); }; @@ -336,8 +336,8 @@ void FiniteElement::fill_fe_face_values (const DoFHandler::cell_iterat const vector > &global_unit_points, vector &jacobians, const bool compute_jacobians, - vector > &ansatz_points, - const bool compute_ansatz_points, + vector > &support_points, + const bool compute_support_points, vector > &q_points, const bool compute_q_points, vector &face_jacobi_determinants, @@ -353,8 +353,8 @@ void FiniteElement::fill_fe_face_values (const DoFHandler::cell_iterat 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 > dummy(total_dofs); fill_fe_values (cell, global_unit_points, @@ -364,8 +364,8 @@ void FiniteElement::fill_fe_face_values (const DoFHandler::cell_iterat 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, @@ -423,14 +423,14 @@ void FiniteElement::fill_fe_subface_values (const DoFHandler::cell_ite template -void FiniteElement::get_unit_ansatz_points (vector > &) const { +void FiniteElement::get_unit_support_points (vector > &) const { Assert (false, ExcPureFunctionCalled()); }; template -void FiniteElement::get_ansatz_points (const DoFHandler::cell_iterator &, +void FiniteElement::get_support_points (const DoFHandler::cell_iterator &, const Boundary &, vector > &) const { Assert (false, ExcPureFunctionCalled()); diff --git a/deal.II/deal.II/source/fe/fe_linear_mapping.cc b/deal.II/deal.II/source/fe/fe_linear_mapping.cc index dc44a59619..5e1b2cb8b8 100644 --- a/deal.II/deal.II/source/fe/fe_linear_mapping.cc +++ b/deal.II/deal.II/source/fe/fe_linear_mapping.cc @@ -97,8 +97,8 @@ void FELinearMapping<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cel const vector > &unit_points, vector &jacobians, const bool compute_jacobians, - vector > &ansatz_points, - const bool compute_ansatz_points, + vector > &support_points, + const bool compute_support_points, vector > &q_points, const bool compute_q_points, const dFMatrix &shape_values_transform, @@ -107,7 +107,7 @@ void FELinearMapping<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cel // 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); @@ -281,8 +281,8 @@ void FELinearMapping::fill_fe_values (const DoFHandler::cell_iterator const vector > &unit_points, vector &jacobians, const bool compute_jacobians, - vector > &ansatz_points, - const bool compute_ansatz_points, + vector > &support_points, + const bool compute_support_points, vector > &q_points, const bool compute_q_points, const dFMatrix &shape_values_transform, @@ -292,8 +292,8 @@ void FELinearMapping::fill_fe_values (const DoFHandler::cell_iterator 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(); @@ -364,8 +364,8 @@ void FELinearMapping::fill_fe_values (const DoFHandler::cell_iterator }; }; - if (compute_ansatz_points) - get_ansatz_points (cell, boundary, ansatz_points); + if (compute_support_points) + get_support_points (cell, boundary, support_points); }; diff --git a/deal.II/deal.II/source/numerics/matrices.cc b/deal.II/deal.II/source/numerics/matrices.cc index 482ae8e962..425f34dae9 100644 --- a/deal.II/deal.II/source/numerics/matrices.cc +++ b/deal.II/deal.II/source/numerics/matrices.cc @@ -797,11 +797,22 @@ void LaplaceMatrix::assemble (dVector &rhs, rhs(i) += values(i,point) * rhs_values[point] * weights[point]; -}; - - +} +template void +MatrixCreator::create_interpolation_matrix(const FiniteElement &high, + const FiniteElement &low, + dFMatrix& result) +{ + result.reinit (high.total_dofs, low.total_dofs); + vector > unit_support_points (high.total_dofs); + high.get_unit_support_points (unit_support_points); + + for (unsigned int i=0; i::interpolate (const DoFHandler &dof, endc = dof.end(); vector dofs_on_cell (fe.total_dofs); vector dof_values_on_cell (fe.total_dofs); - vector > ansatz_points (fe.total_dofs); + vector > 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 @@ -266,7 +266,7 @@ VectorTools::interpolate_boundary_values (const DoFHandler &dof, // 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 @@ -396,7 +396,7 @@ void VectorTools::integrate_difference (const DoFHandler &dof, // \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 -- 2.39.5