From 755e8ea8daf73eb13c7a0f275b273932dcb5e1b0 Mon Sep 17 00:00:00 2001 From: wolf Date: Sat, 7 Nov 1998 02:15:57 +0000 Subject: [PATCH] More work with regard to the transition of gradients from data type Point to Tensor<1,dim>. git-svn-id: https://svn.dealii.org/trunk@655 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe.h | 10 +- .../deal.II/include/fe/fe_lib.criss_cross.h | 14 +-- .../deal.II/include/fe/fe_linear_mapping.h | 6 +- deal.II/deal.II/include/fe/fe_values.h | 10 +- deal.II/deal.II/include/grid/tria.h | 4 +- deal.II/deal.II/source/dofs/dof_handler.cc | 10 +- deal.II/deal.II/source/fe/fe.cc | 12 +- .../deal.II/source/fe/fe_lib.criss_cross.cc | 111 ++++++++++------- deal.II/deal.II/source/fe/fe_lib.linear.cc | 14 ++- .../deal.II/source/fe/fe_linear_mapping.cc | 115 ++++++++++-------- deal.II/deal.II/source/fe/fe_values.cc | 25 ++-- .../source/numerics/error_estimator.cc | 12 +- deal.II/deal.II/source/numerics/vectors.cc | 10 +- 13 files changed, 196 insertions(+), 157 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index 89f9270928..69a3bc511d 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -578,8 +578,8 @@ class FiniteElement : public FiniteElementBase { * #i#th shape function for the mapping * from unit to real cell. */ - virtual Point shape_grad_transform (const unsigned int i, - const Point &p) const = 0; + virtual Tensor<1,dim> shape_grad_transform (const unsigned int i, + const Point &p) const = 0; /** * Compute the Jacobian matrix and the @@ -646,7 +646,7 @@ class FiniteElement : public FiniteElementBase { vector > &q_points, const bool compute_q_points, const dFMatrix &shape_values_transform, - const vector > > &shape_grads_transform, + const vector > > &shape_grads_transform, const Boundary &boundary) const; /** @@ -768,7 +768,7 @@ class FiniteElement : public FiniteElementBase { vector > &normal_vectors, const bool compute_normal_vectors, const dFMatrix &shape_values_transform, - const vector > > &shape_grads_transform, + const vector > > &shape_grads_transform, const Boundary &boundary) const; /** @@ -813,7 +813,7 @@ class FiniteElement : public FiniteElementBase { vector > &normal_vectors, const bool compute_normal_vectors, const dFMatrix &shape_values_transform, - const vector > > &shape_grads_transform, + const vector > > &shape_grads_transform, const Boundary &boundary) const; /** diff --git a/deal.II/deal.II/include/fe/fe_lib.criss_cross.h b/deal.II/deal.II/include/fe/fe_lib.criss_cross.h index d27931ed96..68d41b8ea3 100644 --- a/deal.II/deal.II/include/fe/fe_lib.criss_cross.h +++ b/deal.II/deal.II/include/fe/fe_lib.criss_cross.h @@ -210,16 +210,16 @@ class FECrissCross : public FiniteElement { * information on this function. */ virtual void get_support_points (const DoFHandler::cell_iterator &cell, - const Boundary &boundary, - vector > &support_points) const; + const Boundary &boundary, + vector > &support_points) const; /** * Refer to the base class for detailed * information on this function. */ virtual void get_face_support_points (const DoFHandler::face_iterator &face, - const Boundary &boundary, - vector > &support_points) const; + const Boundary &boundary, + vector > &support_points) const; /** * Refer to the base class for detailed @@ -248,8 +248,8 @@ class FECrissCross : public FiniteElement { * computation of the transformation from * unit cell to real space cell. */ - virtual Point shape_grad_transform (const unsigned int i, - const Point &p) const; + virtual Tensor<1,dim> shape_grad_transform (const unsigned int i, + const Point &p) const; /** * Refer to the base class for detailed @@ -339,7 +339,7 @@ class FECrissCross : public FiniteElement { vector > &q_points, const bool compute_q_points, const dFMatrix &shape_values_transform, - const vector > > &shape_grad_transform, + const vector > > &shape_grad_transform, const Boundary &boundary) const; DeclException0 (ExcNotUseful); diff --git a/deal.II/deal.II/include/fe/fe_linear_mapping.h b/deal.II/deal.II/include/fe/fe_linear_mapping.h index 735023da01..9be39b199f 100644 --- a/deal.II/deal.II/include/fe/fe_linear_mapping.h +++ b/deal.II/deal.II/include/fe/fe_linear_mapping.h @@ -51,8 +51,8 @@ class FELinearMapping : public FiniteElement { * computation of the transformation from * unit cell to real space cell. */ - virtual Point shape_grad_transform (const unsigned int i, - const Point &p) const; + virtual Tensor<1,dim> shape_grad_transform (const unsigned int i, + const Point &p) const; /** * Refer to the base class for detailed @@ -142,7 +142,7 @@ class FELinearMapping : public FiniteElement { vector > &q_points, const bool compute_q_points, const dFMatrix &shape_values_transform, - const vector > > &shape_grad_transform, + const vector > > &shape_grad_transform, const Boundary &boundary) const; }; diff --git a/deal.II/deal.II/include/fe/fe_values.h b/deal.II/deal.II/include/fe/fe_values.h index 12fa99453f..303c126a24 100644 --- a/deal.II/deal.II/include/fe/fe_values.h +++ b/deal.II/deal.II/include/fe/fe_values.h @@ -297,7 +297,7 @@ class FEValuesBase { * function at the #j# quadrature point. * If you want to get the derivative in * one of the coordinate directions, use - * the appropriate function of the #Point# + * the appropriate function of the #Tensor# * class to extract one component. Since * only a reference to the gradient's value * is returned, there should be no major @@ -327,8 +327,8 @@ class FEValuesBase { * #gradients# object already has the * right size. */ - void get_function_grads (const dVector &fe_function, - vector > &gradients) const; + void get_function_grads (const dVector &fe_function, + vector > &gradients) const; /** * Return the position of the #i#th @@ -704,7 +704,7 @@ class FEValues : public FEValuesBase { * Analogous to the #shape_values_transform# * array of the base class. */ - vector > > unit_shape_gradients_transform; + vector > > unit_shape_gradients_transform; /** * Array of quadrature points in the unit @@ -861,7 +861,7 @@ class FEFaceValuesBase : public FEValuesBase { * Analogous to the #shape_values_transform# * array of the base class. */ - vector > > > unit_shape_gradients_transform; + vector > > > unit_shape_gradients_transform; /** * Array of quadrature points on the diff --git a/deal.II/deal.II/include/grid/tria.h b/deal.II/deal.II/include/grid/tria.h index f670670333..1cb011c71f 100644 --- a/deal.II/deal.II/include/grid/tria.h +++ b/deal.II/deal.II/include/grid/tria.h @@ -8,7 +8,7 @@ #include #include #include - +#include //forward declaration needed @@ -1144,7 +1144,7 @@ class TriaDimensionInfo<2> { * @author Wolfgang Bangerth, 1998 */ template -class Triangulation : public TriaDimensionInfo { +class Triangulation : public TriaDimensionInfo, public Subscriptor { public: // insert these definitions for gcc2.8, // since it can't inherit typedefs (I diff --git a/deal.II/deal.II/source/dofs/dof_handler.cc b/deal.II/deal.II/source/dofs/dof_handler.cc index fb090012be..5a79b35868 100644 --- a/deal.II/deal.II/source/dofs/dof_handler.cc +++ b/deal.II/deal.II/source/dofs/dof_handler.cc @@ -25,12 +25,16 @@ template DoFHandler::DoFHandler (Triangulation *tria) : tria(tria), - used_dofs (0) {}; + used_dofs (0) +{ + tria->subscribe (); +}; template -DoFHandler::~DoFHandler () -{} +DoFHandler::~DoFHandler () { + tria->unsubscribe (); +} diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index 0833482a4a..57151da0e9 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -185,7 +185,7 @@ void FiniteElement<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell, vector > &q_points, const bool compute_q_points, const dFMatrix &, - const vector > > &, + const vector > > &, const Boundary<1> &boundary) const { Assert (jacobians.size() == unit_points.size(), ExcWrongFieldDimension(jacobians.size(), unit_points.size())); @@ -235,7 +235,7 @@ void FiniteElement<1>::fill_fe_face_values (const DoFHandler<1>::cell_iterator & vector > &, const bool, const dFMatrix &, - const vector > > &, + const vector > > &, const Boundary<1> &) const { Assert (false, ExcNotImplemented()); }; @@ -256,7 +256,7 @@ void FiniteElement<1>::fill_fe_subface_values (const DoFHandler<1>::cell_iterato vector > &, const bool, const dFMatrix &, - const vector > > &, + const vector > > &, const Boundary<1> &) const { Assert (false, ExcNotImplemented()); }; @@ -323,7 +323,7 @@ void FiniteElement::fill_fe_values (const DoFHandler::cell_iterator &, vector > &, const bool, const dFMatrix &, - const vector > > &, + const vector > > &, const Boundary &) const { Assert (false, ExcPureFunctionCalled()); }; @@ -346,7 +346,7 @@ void FiniteElement::fill_fe_face_values (const DoFHandler::cell_iterat vector > &normal_vectors, const bool compute_normal_vectors, const dFMatrix &shape_values_transform, - const vector > > &shape_gradients_transform, + const vector > > &shape_gradients_transform, const Boundary &boundary) const { Assert (jacobians.size() == unit_points.size(), ExcWrongFieldDimension(jacobians.size(), unit_points.size())); @@ -395,7 +395,7 @@ void FiniteElement::fill_fe_subface_values (const DoFHandler::cell_ite vector > &normal_vectors, const bool compute_normal_vectors, const dFMatrix &shape_values_transform, - const vector > > &shape_gradients_transform, + const vector > > &shape_gradients_transform, const Boundary &boundary) const { Assert (jacobians.size() == unit_points.size(), ExcWrongFieldDimension(jacobians.size(), unit_points.size())); diff --git a/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc b/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc index 5cace7e1d7..0af836516f 100644 --- a/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc +++ b/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc @@ -347,8 +347,8 @@ void FECrissCross<1>::get_unit_support_points (vector >&) const { template <> void FECrissCross<1>::get_support_points (const DoFHandler<1>::cell_iterator &, - const Boundary<1> &, - vector > &) const { + const Boundary<1> &, + vector > &) const { Assert (false, ExcNotUseful()); }; @@ -356,8 +356,8 @@ void FECrissCross<1>::get_support_points (const DoFHandler<1>::cell_iterator &, template <> void FECrissCross<1>::get_face_support_points (const DoFHandler<1>::face_iterator &, - const Boundary<1> &, - vector > &) const { + const Boundary<1> &, + vector > &) const { Assert (false, ExcNotUseful()); }; @@ -382,8 +382,8 @@ double FECrissCross<1>::shape_value_transform (const unsigned int, template <> -Point<1> FECrissCross<1>::shape_grad_transform (const unsigned int, - const Point<1> &) const { +Tensor<1,1> FECrissCross<1>::shape_grad_transform (const unsigned int, + const Point<1> &) const { Assert (false, ExcNotUseful()); return Point<1>(); }; @@ -442,7 +442,7 @@ void FECrissCross<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &, vector > &, const bool, const dFMatrix &, - const vector > > &, + const vector > > &, const Boundary<1> &) const { Assert (false, ExcNotUseful()); }; @@ -594,8 +594,8 @@ void FECrissCross<2>::get_unit_support_points (vector > &unit_points) c template <> void FECrissCross<2>::get_support_points (const DoFHandler<2>::cell_iterator &cell, - const Boundary<2> &, - vector > &support_points) const { + const Boundary<2> &, + vector > &support_points) const { const unsigned int dim = 2; Assert (support_points.size() == total_dofs, @@ -632,8 +632,8 @@ void FECrissCross<2>::get_support_points (const DoFHandler<2>::cell_iterator &ce template <> void FECrissCross<2>::get_face_support_points (const DoFHandler<2>::face_iterator &face, - const Boundary<2> &, - vector > &support_points) const { + const Boundary<2> &, + vector > &support_points) const { const unsigned int dim = 2; Assert ((support_points.size() == dofs_per_face) && @@ -793,7 +793,7 @@ void FECrissCross<2>::get_local_mass_matrix (const DoFHandler<2>::cell_iterator template <> inline double FECrissCross<2>::shape_value_transform (const unsigned int i, - const Point<2> &p) const { + const Point<2> &p) const { // use an isoparametric ansatz return shape_value(i,p); }; @@ -801,8 +801,8 @@ double FECrissCross<2>::shape_value_transform (const unsigned int i, template <> -Point<2> FECrissCross<2>::shape_grad_transform (const unsigned int i, - const Point<2> &p) const { +Tensor<1,2> FECrissCross<2>::shape_grad_transform (const unsigned int i, + const Point<2> &p) const { // use an isoparametric ansatz return shape_grad(i,p); }; @@ -934,7 +934,7 @@ void FECrissCross::fill_fe_values (const DoFHandler::cell_iterator &ce vector > &q_points, const bool compute_q_points, const dFMatrix &shape_values_transform, - const vector > > &/*shape_grad_transform*/, + const vector > > &/*shape_grad_transform*/, const Boundary &boundary) const { Assert (jacobians.size() == unit_points.size(), ExcWrongFieldDimension(jacobians.size(), unit_points.size())); @@ -999,40 +999,57 @@ void FECrissCross::fill_fe_values (const DoFHandler::cell_iterator &ce vertices[l] = cell->vertex(l); if (compute_jacobians) - for (unsigned int point=0; point::get_unit_support_points (vector > &unit_points) const template -void FELinear::get_support_points (const typename DoFHandler::cell_iterator &cell, - const Boundary &, - vector > &support_points) const { +void +FELinear::get_support_points (const typename DoFHandler::cell_iterator &cell, + const Boundary &, + vector > &support_points) const { Assert (support_points.size() == total_dofs, ExcWrongFieldDimension (support_points.size(), total_dofs)); @@ -392,9 +393,10 @@ void FELinear::get_support_points (const typename DoFHandler::cell_ite template -void FELinear::get_face_support_points (const typename DoFHandler::face_iterator &face, - const Boundary &, - vector > &support_points) const { +void +FELinear::get_face_support_points (const typename DoFHandler::face_iterator &face, + const Boundary &, + vector > &support_points) const { Assert ((support_points.size() == dofs_per_face) && (support_points.size() == GeometryInfo::vertices_per_face), ExcWrongFieldDimension (support_points.size(), 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 e2728ad763..f9694d653f 100644 --- a/deal.II/deal.II/source/fe/fe_linear_mapping.cc +++ b/deal.II/deal.II/source/fe/fe_linear_mapping.cc @@ -36,7 +36,7 @@ FELinearMapping<1>::shape_value_transform (const unsigned int i, template <> inline -Point<1> +Tensor<1,1> FELinearMapping<1>::shape_grad_transform(const unsigned int i, const Point<1>&) const { @@ -102,7 +102,7 @@ void FELinearMapping<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cel vector > &q_points, const bool compute_q_points, const dFMatrix &shape_values_transform, - const vector > > &shape_gradients_transform, + const vector > > &shape_gradients_transform, const Boundary<1> &boundary) const { // simply pass down FiniteElement<1>::fill_fe_values (cell, unit_points, @@ -142,7 +142,7 @@ FELinearMapping<2>::shape_value_transform (const unsigned int i, template <> inline -Point<2> +Tensor<1,2> FELinearMapping<2>::shape_grad_transform (const unsigned int i, const Point<2>& p) const { @@ -286,7 +286,7 @@ void FELinearMapping::fill_fe_values (const DoFHandler::cell_iterator vector > &q_points, const bool compute_q_points, const dFMatrix &shape_values_transform, - const vector > > &/*shape_grad_transform*/, + const vector > > &/*shape_grad_transform*/, const Boundary &boundary) const { Assert (jacobians.size() == unit_points.size(), ExcWrongFieldDimension(jacobians.size(), unit_points.size())); @@ -379,11 +379,19 @@ void FELinearMapping::fill_fe_values (const DoFHandler::cell_iterator is multiplied to the unit cell gradient *from the right*! be very careful with these things. - The following little program tests the correct behaviour (you have to find - out the right include files, I tested it within a whole project with far - more include files than necessary): + The following little program tests the correct behaviour: ------------------------------------------- + #include + #include + #include + #include + #include + #include + #include + #include + #include + int main () { Triangulation<2> tria; tria.create_hypercube (0,1); @@ -428,53 +436,58 @@ void FELinearMapping::fill_fe_values (const DoFHandler::cell_iterator */ if (compute_jacobians) - for (unsigned int point=0; point::shape_grad (const unsigned int i, template -void FEValuesBase::get_function_grads (const dVector &fe_function, - vector > &gradients) const { +void FEValuesBase::get_function_grads (const dVector &fe_function, + vector > &gradients) const { Assert (gradients.size() == n_quadrature_points, ExcWrongVectorSize(gradients.size(), n_quadrature_points)); @@ -112,7 +112,7 @@ void FEValuesBase::get_function_grads (const dVector &fe_function, present_cell->get_dof_values (fe_function, dof_values); // initialize with zero - fill_n (gradients.begin(), n_quadrature_points, Point()); + fill_n (gradients.begin(), n_quadrature_points, Tensor<1,dim>()); // add up contributions of trial // functions @@ -128,7 +128,8 @@ void FEValuesBase::get_function_grads (const dVector &fe_function, template -const Point & FEValuesBase::quadrature_point (const unsigned int i) const { +const Point & +FEValuesBase::quadrature_point (const unsigned int i) const { Assert (i & FEValuesBase::quadrature_point (const unsigned int i) co template -const Point & FEValuesBase::support_point (const unsigned int i) const { +const Point & +FEValuesBase::support_point (const unsigned int i) const { Assert (i::FEValues (const FiniteElement &fe, unit_shape_2nd_derivatives(fe.total_dofs, vector >(quadrature.n_quadrature_points)), unit_shape_gradients_transform(fe.n_transform_functions, - vector >(quadrature.n_quadrature_points)), + vector >(quadrature.n_quadrature_points)), unit_quadrature_points(quadrature.get_quad_points()) { Assert ((update_flags & update_normal_vectors) == false, @@ -320,8 +322,8 @@ FEFaceValuesBase::FEFaceValuesBase (const unsigned int n_q_points, vector > >(n_dofs, vector >(n_q_points))), unit_shape_gradients_transform (n_faces_or_subfaces, - vector > >(n_transform_functions, - vector >(n_q_points))), + vector > >(n_transform_functions, + vector >(n_q_points))), unit_face_quadrature_points (n_q_points, Point()), unit_quadrature_points (n_faces_or_subfaces, vector >(n_q_points, Point())), @@ -332,7 +334,8 @@ FEFaceValuesBase::FEFaceValuesBase (const unsigned int n_q_points, template -const Point & FEFaceValuesBase::normal_vector (const unsigned int i) const { +const Point & +FEFaceValuesBase::normal_vector (const unsigned int i) const { Assert (i::reinit (const typename DoFHandler::cell_iterator &c { fill_n (shape_gradients[i].begin(), n_quadrature_points, - Point()); + Tensor<1,dim>()); for (unsigned int j=0; j::reinit (const typename DoFHandler::cell_iterator { fill_n (shape_gradients[i].begin(), n_quadrature_points, - Point()); + Tensor<1,dim>()); for (unsigned int j=0; j > psi(n_q_points); + vector > psi(n_q_points); fe_face_values_cell.get_function_grads (solution, psi); @@ -280,14 +280,14 @@ integrate_over_regular_face (const active_cell_iterator &cell, // get a list of the gradients of // the finite element solution // restricted to the neighbor cell - vector > neighbor_psi (n_q_points); + vector > neighbor_psi (n_q_points); fe_face_values_neighbor.get_function_grads (solution, neighbor_psi); // compute the jump in the gradients transform (psi.begin(), psi.end(), neighbor_psi.begin(), psi.begin(), - minus >()); + minus >()); }; @@ -397,7 +397,7 @@ integrate_over_irregular_face (const active_cell_iterator &cell, // // let psi be a short name for // [a grad u_h] - vector > psi(n_q_points); + vector > psi(n_q_points); // store which number #cell# has in the // list of neighbors of #neighbor# @@ -434,7 +434,7 @@ integrate_over_irregular_face (const active_cell_iterator &cell, // restrict the finite element on the // neighbor cell to the common #subface#. // store the gradient in #neighbor_psi# - vector > neighbor_psi (n_q_points); + vector > neighbor_psi (n_q_points); fe_face_values.reinit (neighbor_child, neighbor_neighbor, boundary); fe_face_values.get_function_grads (solution, neighbor_psi); @@ -442,7 +442,7 @@ integrate_over_irregular_face (const active_cell_iterator &cell, transform (psi.begin(), psi.end(), neighbor_psi.begin(), psi.begin(), - minus >()); + minus >()); // next we have to multiply this with // the normal vector. Since we have diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc index 9f21f6718d..3e7ce4b779 100644 --- a/deal.II/deal.II/source/numerics/vectors.cc +++ b/deal.II/deal.II/source/numerics/vectors.cc @@ -32,8 +32,8 @@ inline double sqr (const double x) { template -inline double sqr_point (const Point &p) { - return p.square(); +inline double sqr_point (const Tensor<1,dim> &p) { + return p * p; }; @@ -507,7 +507,7 @@ void VectorTools::integrate_difference (const DoFHandler &dof, // same procedure as above, but now // psi is a vector of gradients const unsigned int n_q_points = q.n_quadrature_points; - vector > psi (n_q_points); + vector > psi (n_q_points); // in praxi: first compute // exact fe_function vector @@ -518,13 +518,13 @@ void VectorTools::integrate_difference (const DoFHandler &dof, // fe_function if (true) { - vector > function_grads (n_q_points, Point()); + vector > function_grads (n_q_points, Tensor<1,dim>()); fe_values.get_function_grads (fe_function, function_grads); transform (psi.begin(), psi.end(), function_grads.begin(), psi.begin(), - minus >()); + minus >()); }; // take square of integrand vector psi_square (psi.size(), 0.0); -- 2.39.5