From: Markus Buerg Date: Thu, 26 Aug 2010 17:14:16 +0000 (+0000) Subject: Parallelization of FETools::compute_embedding_matrices X-Git-Tag: v8.0.0~5659 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=37db1856a6edc9cb579781f9f70906fbfc591ddf;p=dealii.git Parallelization of FETools::compute_embedding_matrices git-svn-id: https://svn.dealii.org/trunk@21739 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/source/fe/fe_tools.cc b/deal.II/deal.II/source/fe/fe_tools.cc index f6567234fe..7a9d4320c1 100644 --- a/deal.II/deal.II/source/fe/fe_tools.cc +++ b/deal.II/deal.II/source/fe/fe_tools.cc @@ -597,45 +597,38 @@ FETools::compute_embedding_matrices(const FiniteElement<2,3>&, #endif -// This function is tested by tests/fe/internals, since it produces the matrices printed there -template -void -FETools::compute_embedding_matrices(const FiniteElement& fe, - std::vector > >& matrices, - const bool isotropic_only) -{ - const unsigned int n = fe.dofs_per_cell; - const unsigned int nd = fe.n_components(); - const unsigned int degree = fe.degree; - - // loop over all possible refinement cases - unsigned int ref_case = (isotropic_only) - ? RefinementCase::isotropic_refinement - : RefinementCase::cut_x; - - for (;ref_case <= RefinementCase::isotropic_refinement; ++ref_case) - { - const unsigned int nc = GeometryInfo::n_children(RefinementCase(ref_case)); - for (unsigned int i=0;i + void + compute_embedding_matrices_for_refinement_case (const FiniteElement& fe, + std::vector >& matrices, + const unsigned int ref_case) + { + const unsigned int n = fe.dofs_per_cell; + const unsigned int nc = GeometryInfo::n_children(RefinementCase(ref_case)); + for (unsigned int i = 0; i < nc; ++i) + { + Assert(matrices[i].n() == n, ExcDimensionMismatch(matrices[i].n (), n)); + Assert(matrices[i].m() == n, ExcDimensionMismatch(matrices[i].m (), n)); + } // Set up meshes, one with a single // reference cell and refine it once - Triangulation tria; - GridGenerator::hyper_cube (tria, 0, 1); - tria.begin_active()->set_refine_flag(RefinementCase(ref_case)); - tria.execute_coarsening_and_refinement(); + Triangulation tria; + GridGenerator::hyper_cube (tria, 0, 1); + tria.begin_active()->set_refine_flag (RefinementCase(ref_case)); + tria.execute_coarsening_and_refinement (); - MappingCartesian mapping; - QGauss q_fine(degree+1); - const unsigned int nq = q_fine.size(); + MappingCartesian mapping; + const unsigned int degree = fe.degree; + QGauss q_fine (degree+1); + const unsigned int nq = q_fine.size(); - FEValues fine (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_values); // We search for the polynomial on // the small cell, being equal to @@ -650,63 +643,70 @@ FETools::compute_embedding_matrices(const FiniteElement& fe, // This matrix is the same for all // children. - fine.reinit(tria.begin_active()); - FullMatrix A(nq*nd, n); - for (unsigned int j=0;j H(A); - - Vector v_coarse(nq*nd); - Vector v_fine(n); - - unsigned int cell_number = 0; - for (typename Triangulation::active_cell_iterator fine_cell - = tria.begin_active(); - fine_cell != tria.end(); ++fine_cell, ++cell_number) - { - fine.reinit(fine_cell); + fine.reinit (tria.begin_active ()); + const unsigned int nd = fe.n_components (); + FullMatrix A (nq*nd, n); + + for (unsigned int j = 0; j < n; ++j) + for (unsigned int d = 0; d < nd; ++d) + for (unsigned int k = 0; k < nq; ++k) + A (k * nd + d, j) = fine.shape_value_component (j, k, d); + + Householder H (A); + static Threads::Mutex mutex; + Vector v_coarse (nq * nd); + Vector v_fine (n); + unsigned int cell_number = 0; + + for (typename Triangulation::active_cell_iterator + fine_cell = tria.begin_active (); fine_cell != tria.end (); + ++fine_cell, ++cell_number) + { + fine.reinit (fine_cell); // evaluate on the coarse cell (which // is the first -- inactive -- cell on // the lowest level of the // triangulation we have created) - const Quadrature q_coarse (fine.get_quadrature_points(), - fine.get_JxW_values()); - FEValues coarse (mapping, fe, q_coarse, update_values); - coarse.reinit(tria.begin(0)); + const Quadrature q_coarse (fine.get_quadrature_points (), + fine.get_JxW_values ()); + FEValues coarse (mapping, fe, q_coarse, update_values); + + coarse.reinit (tria.begin (0)); - FullMatrix &this_matrix = matrices[ref_case-1][cell_number]; - v_coarse = 0; + FullMatrix &this_matrix = matrices[cell_number]; + + v_coarse = 0; // Compute this once for each // coarse grid basis function - for (unsigned int i=0;i& fe, // function to a fine grid // function, the columns // are fine grid. - for (unsigned int j=0;j::n_children(RefinementCase(ref_case)), - ExcInternalError()); - } + for (unsigned int i = 0; i < this_matrix.m (); ++i) + for (unsigned int j = 0; j < this_matrix.n (); ++j) + if (std::fabs (this_matrix (i, j)) < 1e-12) + this_matrix (i, j) = 0.; + + mutex.release (); + } + + Assert (cell_number == GeometryInfo::n_children (RefinementCase (ref_case)), + ExcInternalError ()); + } +} + + +// This function is tested by tests/fe/internals, since it produces the matrices printed there +template +void +FETools::compute_embedding_matrices(const FiniteElement& fe, + std::vector > >& matrices, + const bool isotropic_only) +{ + Threads::TaskGroup task_group; + + // loop over all possible refinement cases + unsigned int ref_case = (isotropic_only) + ? RefinementCase::isotropic_refinement + : RefinementCase::cut_x; + + for (;ref_case <= RefinementCase::isotropic_refinement; ++ref_case) + task_group += Threads::new_task (&compute_embedding_matrices_for_refinement_case, + fe, matrices[ref_case-1], ref_case); + + task_group.join_all (); } diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index f5b3bb9703..b0fab3d28c 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -230,6 +230,14 @@ through DoFHandler::get_tria() and DoFHandler::get_fe(), respectively.
    +
  1. +

    + New: FETools::compute_embedding_matrices now computes the embedding matrix + for every refinement case in its own task. +
    + (Markus Buerg 2010/08/26) +

    +
  2. New: Data structures and a function parent to get the parent of a cell.