From: Martin Kronbichler Date: Fri, 5 Jun 2009 11:10:59 +0000 (+0000) Subject: Let shape_value return a reference to the data, and not only the value. This allows... X-Git-Tag: v8.0.0~7633 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e64c473f75c0c27d23b4b74000eabe8bf3f23ae1;p=dealii.git Let shape_value return a reference to the data, and not only the value. This allows us to use some pointers in performance-critical parts of MatrixCreator and FETools. Moreover, use symmetry in mass and Laplace matrix assembly. git-svn-id: https://svn.dealii.org/trunk@18914 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/fe_values.h b/deal.II/deal.II/include/fe/fe_values.h index 5efbd52402..ee712fc875 100644 --- a/deal.II/deal.II/include/fe/fe_values.h +++ b/deal.II/deal.II/include/fe/fe_values.h @@ -1298,8 +1298,8 @@ class FEValuesBase : protected FEValuesData, * the quadrature point at which * function is to be computed */ - double shape_value (const unsigned int function_no, - const unsigned int point_no) const; + const double & shape_value (const unsigned int function_no, + const unsigned int point_no) const; /** * Compute one vector component of @@ -3789,7 +3789,7 @@ operator[] (const FEValuesExtractors::Vector &vector) const template inline -double +const double & FEValuesBase::shape_value (const unsigned int i, const unsigned int j) const { diff --git a/deal.II/deal.II/source/fe/fe_tools.cc b/deal.II/deal.II/source/fe/fe_tools.cc index 5df260a7ce..9e2c36dbd2 100644 --- a/deal.II/deal.II/source/fe/fe_tools.cc +++ b/deal.II/deal.II/source/fe/fe_tools.cc @@ -637,9 +637,9 @@ FETools::compute_embedding_matrices(const FiniteElement& fe, // children. fine.reinit(tria.begin_active()); FullMatrix A(nq*nd, n); - for (unsigned int d=0;d H(A); @@ -664,6 +664,7 @@ FETools::compute_embedding_matrices(const FiniteElement& fe, coarse.reinit(tria.begin(0)); FullMatrix &this_matrix = matrices[ref_case-1][cell_number]; + v_coarse = 0; // Compute this once for each // coarse grid basis function @@ -675,9 +676,17 @@ FETools::compute_embedding_matrices(const FiniteElement& fe, // function values of the // coarse grid function in // each quadrature point. - for (unsigned int d=0;d& fe, // children. fine.reinit(tria.begin_active()); FullMatrix A(nq*nd, n); - for (unsigned int d=0;d H(A); @@ -849,7 +858,7 @@ FETools::compute_face_embedding_matrices(const FiniteElement& fe, // each quadrature point. for (unsigned int d=0;d& fe, const unsigned int nd = fe.n_components(); const unsigned int degree = fe.degree; + // prepare FEValues, quadrature etc on + // coarse cell + MappingCartesian mapping; + QGauss q_fine(degree+1); + const unsigned int nq = q_fine.size(); + + // create mass matrix on coarse cell. + FullMatrix mass(n, n); + { + // set up a triangulation for coarse cell + Triangulation tr; + GridGenerator::hyper_cube (tr, 0, 1); + + FEValues coarse (mapping, fe, q_fine, + update_JxW_values | update_values); + + typename Triangulation::cell_iterator coarse_cell + = tr.begin(0); + coarse.reinit (coarse_cell); + + const std::vector & JxW = coarse.get_JxW_values(); + for (unsigned int i=0;i::cut_x; - ref_case::isotropic_refinement+1; ++ref_case) + ref_case < RefinementCase::isotropic_refinement+1; ++ref_case) { const unsigned int nc = GeometryInfo::n_children(RefinementCase(ref_case)); @@ -918,45 +973,20 @@ FETools::compute_projection_matrices(const FiniteElement& fe, Assert(matrices[ref_case-1][i].m() == n, ExcDimensionMismatch(matrices[ref_case-1][i].m(),n)); } - + + // create a respective refinement on the + // triangulation Triangulation tr; GridGenerator::hyper_cube (tr, 0, 1); tr.begin_active()->set_refine_flag(RefinementCase(ref_case)); tr.execute_coarsening_and_refinement(); - - MappingCartesian mapping; - QGauss q_fine(degree+1); - const unsigned int nq = q_fine.size(); - - FEValues coarse (mapping, fe, q_fine, - update_quadrature_points | - update_JxW_values | - update_values); + FEValues fine (mapping, fe, q_fine, - update_quadrature_points | - update_JxW_values | + update_quadrature_points | update_JxW_values | update_values); - + typename Triangulation::cell_iterator coarse_cell = tr.begin(0); - // Compute the coarse level mass - // matrix - coarse.reinit(coarse_cell); - FullMatrix A(n, n); - for (unsigned int k=0;k H(A); Vector v_coarse(n); Vector v_fine(n); @@ -977,37 +1007,39 @@ FETools::compute_projection_matrices(const FiniteElement& fe, // Build RHS + const std::vector & JxW = fine.get_JxW_values(); + // Outer loop over all fine // grid shape functions phi_j for (unsigned int j=0;j #include #include -#include -#include +#include +#include #include #include #include -#include -#include -#include #include #include +#include +#include +#include #include #include #include @@ -34,7 +34,6 @@ #include #include #include -#include #include #include @@ -56,10 +55,10 @@ namespace internal { Scratch (const FiniteElement &fe, const UpdateFlags update_flags, - const Function *coefficient, - const Function *rhs_function, + const Function *coefficient, + const Function *rhs_function, const Quadrature &quadrature, - const Mapping &mapping) + const Mapping &mapping) : fe_collection (fe), quadrature_collection (quadrature), @@ -81,10 +80,10 @@ namespace internal Scratch (const ::dealii::hp::FECollection &fe, const UpdateFlags update_flags, - const Function *coefficient, - const Function *rhs_function, - const ::dealii::hp::QCollection &quadrature, - const ::dealii::hp::MappingCollection &mapping) + const Function *coefficient, + const Function *rhs_function, + const ::dealii::hp::QCollection &quadrature, + const ::dealii::hp::MappingCollection &mapping) : fe_collection (fe), quadrature_collection (quadrature), @@ -122,8 +121,8 @@ namespace internal update_flags (data.update_flags) {} - const ::dealii::hp::FECollection fe_collection; - const ::dealii::hp::QCollection quadrature_collection; + const ::dealii::hp::FECollection fe_collection; + const ::dealii::hp::QCollection quadrature_collection; const ::dealii::hp::MappingCollection mapping_collection; ::dealii::hp::FEValues x_fe_values; @@ -133,6 +132,8 @@ namespace internal std::vector rhs_values; std::vector > rhs_vector_values; + std::vector old_JxW; + const Function *coefficient; const Function *rhs_function; @@ -164,7 +165,7 @@ namespace internal const unsigned int dofs_per_cell = fe_values.dofs_per_cell, n_q_points = fe_values.n_quadrature_points; - const FiniteElement &fe = fe_values.get_fe(); + const FiniteElement &fe = fe_values.get_fe(); const unsigned int n_components = fe.n_components(); Assert(data.rhs_function == 0 || @@ -177,15 +178,13 @@ namespace internal ::dealii::MatrixCreator::ExcComponentMismatch()); copy_data.cell_matrix.reinit (dofs_per_cell, dofs_per_cell); - copy_data.cell_matrix = 0; - copy_data.cell_rhs.reinit (dofs_per_cell); - copy_data.cell_rhs = 0; copy_data.dof_indices.resize (dofs_per_cell); cell->get_dof_indices (copy_data.dof_indices); - if (data.rhs_function != 0) + const bool use_rhs_function = data.rhs_function != 0; + if (use_rhs_function) { if (data.rhs_function->n_components==1) { @@ -202,7 +201,8 @@ namespace internal } } - if (data.coefficient != 0) + const bool use_coefficient = data.coefficient != 0; + if (use_coefficient) { if (data.coefficient->n_components==1) { @@ -220,175 +220,116 @@ namespace internal } - if (data.coefficient != 0) - { - if (data.coefficient->n_components==1) - { - for (unsigned int i=0; i & JxW = fe_values.get_JxW_values(); + for (unsigned int i=0; in_components==1) + if (data.coefficient->n_components==1) for (unsigned int point=0; pointn_components==1) + for (unsigned int point=0; pointn_components==1) + if (data.coefficient->n_components==1) for (unsigned int point=0; pointn_components==1) - for (unsigned int point=0; pointn_components==1) for (unsigned int point=0; pointn_components==1) - for (unsigned int point=0; pointget_dof_indices (copy_data.dof_indices); - if (data.rhs_function != 0) + const bool use_rhs_function = data.rhs_function != 0; + if (use_rhs_function) { if (data.rhs_function->n_components==1) { @@ -444,7 +382,8 @@ namespace internal } } - if (data.coefficient != 0) + const bool use_coefficient = data.coefficient != 0; + if (use_coefficient) { if (data.coefficient->n_components==1) { @@ -462,179 +401,120 @@ namespace internal } - if (data.coefficient != 0) - { - if (data.coefficient->n_components==1) - { - for (unsigned int i=0; i & JxW = fe_values.get_JxW_values(); + double add_data; + for (unsigned int i=0; i * grad_phi_i = + &fe_values.shape_grad(i,0); + + // can use symmetry + for (unsigned int j=i; j * grad_phi_j = + & fe_values.shape_grad(j,0); + add_data = 0; + if (use_coefficient) { - if (data.rhs_function->n_components==1) + if (data.coefficient->n_components==1) for (unsigned int point=0; pointn_components==1) + for (unsigned int point=0; pointn_components==1) + if (data.coefficient->n_components==1) for (unsigned int point=0; pointn_components==1) - for (unsigned int point=0; pointn_components==1) for (unsigned int point=0; pointn_components==1) - for (unsigned int point=0; point void copy_local_to_global (const AssemblerData::CopyData &data, @@ -652,8 +532,7 @@ namespace internal (data.cell_rhs.size() == dofs_per_cell), ExcInternalError()); - matrix->add (data.dof_indices, data.cell_matrix); - + matrix->add(data.dof_indices, data.cell_matrix); if (right_hand_side != 0) for (unsigned int i=0; i