From dfd37a810e4617263eb6b2439806791fd88d96f9 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 16 May 2000 13:34:23 +0000 Subject: [PATCH] Use the old scheme to compute the Jacobian matrices again. The one which uses Maple output was horribly slow (about 40 times, or so) and also overly complicated. A simple heat equation program ran 4 times faster after this change (the whole program, not only matrix assembly...), so this change is probably justified. git-svn-id: https://svn.dealii.org/trunk@2869 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/q1_mapping.h | 10 - deal.II/deal.II/source/fe/q1_mapping.cc | 57 +- .../deal.II/source/fe/q1_mapping.jacobians.cc | 1203 ----------------- 3 files changed, 29 insertions(+), 1241 deletions(-) diff --git a/deal.II/deal.II/include/fe/q1_mapping.h b/deal.II/deal.II/include/fe/q1_mapping.h index 298edb09b4..9867e1e8ac 100644 --- a/deal.II/deal.II/include/fe/q1_mapping.h +++ b/deal.II/deal.II/include/fe/q1_mapping.h @@ -156,16 +156,6 @@ class FEQ1Mapping : public FiniteElement const FullMatrix &shape_values_transform, const vector > > &shape_grad_transform) const; - /** - * Compute the jacobian matrices - * of the mapping between unit - * and real cell at the given - * points on the unit cell. - */ - static void compute_jacobian_matrices (const DoFHandler::cell_iterator &cell, - const vector > &unit_points, - vector > &jacobians); - /** * Compute the gradients of the jacobian * matrices of the mapping between unit diff --git a/deal.II/deal.II/source/fe/q1_mapping.cc b/deal.II/deal.II/source/fe/q1_mapping.cc index 97689773cc..4545fbf7fc 100644 --- a/deal.II/deal.II/source/fe/q1_mapping.cc +++ b/deal.II/deal.II/source/fe/q1_mapping.cc @@ -543,7 +543,7 @@ void FEQ1Mapping::fill_fe_values (const DoFHandler::cell_iterator &cel vector > &q_points, const bool compute_q_points, const FullMatrix &shape_values_transform, - const vector > > &/*shape_grad_transform*/) const + const vector > > &shape_grad_transform) const { Assert ((!compute_jacobians) || (jacobians.size() == unit_points.size()), ExcWrongFieldDimension(jacobians.size(), unit_points.size())); @@ -604,32 +604,6 @@ void FEQ1Mapping::fill_fe_values (const DoFHandler::cell_iterator &cel However, we rewrite the loops to only compute the gradient once for each integration point and basis function. - Indeed, this was the old way we did it; the code is below. However, there - is a more efficient way, namely setting up M analytically, doing the - inversion analyically and only then doing the evaluation; a small Maple - script does this (it is part of the script in the - subdirectory). - - if (compute_jacobians) - { - FullMatrix M(dim,dim); - for (unsigned int l=0; l::vertices_per_cell; ++s) - { - // we want the linear transform, - // so use that function - const Point gradient = shape_grad_transform[s][l]; - for (unsigned int i=0; i::fill_fe_values (const DoFHandler::cell_iterator &cel */ if (compute_jacobians) - compute_jacobian_matrices (cell, unit_points, jacobians); + { + if (dim == 1) + { + const double h = (cell->vertex(1)(0)-cell->vertex(0)(0)); + + for (unsigned int point=0; point M; + for (unsigned int l=0; l::vertices_per_cell; ++s) + { + // we want the linear transform, + // so use that function + const Point gradient = shape_grad_transform[s][l]; + for (unsigned int i=0; i -void FEQ1Mapping<1>::compute_jacobian_matrices (const DoFHandler<1>::cell_iterator &cell, - const vector > &unit_points, - vector > &jacobians) -{ - Assert (unit_points.size() == jacobians.size(), - ExcWrongFieldDimension(jacobians.size(), unit_points.size())); - - const double h = (cell->vertex(1)(0)-cell->vertex(0)(0)); - - for (unsigned int point=0; point void FEQ1Mapping<1>::compute_jacobian_gradients (const DoFHandler<1>::cell_iterator &, const vector > &unit_points, @@ -67,58 +52,6 @@ void FEQ1Mapping<1>::compute_jacobian_gradients (const DoFHandler<1>::cell_itera #if deal_II_dimension == 2 -template <> -void FEQ1Mapping<2>::compute_jacobian_matrices (const DoFHandler<2>::cell_iterator &cell, - const vector > &unit_points, - vector > &jacobians) -{ - const unsigned int dim = 2; - - Assert (unit_points.size() == jacobians.size(), - ExcWrongFieldDimension(jacobians.size(), unit_points.size())); - - Point vertices[GeometryInfo::vertices_per_cell]; - for (unsigned int l=0; l::vertices_per_cell; ++l) - vertices[l] = cell->vertex(l); - - for (unsigned int point=0; point void FEQ1Mapping<2>::compute_jacobian_gradients (const DoFHandler<2>::cell_iterator &cell, const vector > &unit_points, @@ -202,1142 +135,6 @@ void FEQ1Mapping<2>::compute_jacobian_gradients (const DoFHandler<2>::cell_itera #if deal_II_dimension == 3 -template <> -void FEQ1Mapping<3>::compute_jacobian_matrices (const DoFHandler<3>::cell_iterator &cell, - const vector > &unit_points, - vector > &jacobians) -{ - const unsigned int dim = 3; - - Assert (unit_points.size() == jacobians.size(), - ExcWrongFieldDimension(jacobians.size(), unit_points.size())); - - Point vertices[GeometryInfo::vertices_per_cell]; - for (unsigned int l=0; l::vertices_per_cell; ++l) - vertices[l] = cell->vertex(l); - - for (unsigned int point=0; point void FEQ1Mapping<3>::compute_jacobian_gradients (const DoFHandler<3>::cell_iterator &cell, const vector > &unit_points, -- 2.39.5