From d2e64817135269c6f1fe116193eb354015e63ce7 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 7 Feb 2011 16:31:37 +0000 Subject: [PATCH] Make VectorTools::compute_mean_value and VectorTools::integrate_difference also work in parallel. git-svn-id: https://svn.dealii.org/trunk@23302 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 7 + .../deal.II/numerics/vectors.templates.h | 565 +++++++++--------- 2 files changed, 304 insertions(+), 268 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 883f4931dc..5871044557 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -78,6 +78,13 @@ should be fixed now.

Specific improvements

    +
  1. Fixed: The functions VectorTools::compute_mean_value and +VectorTools::integrate_difference now also work for distributed +triangulations of type parallel::distributed::Triangulation. +
    +(Wolfgang Bangerth, 2011/02/07) +
  2. +
  3. Changed: If the libz library was detected during library configuration, the function DataOutBase::write_vtu now writes data in compressed format, saving a good fraction of disk space (80-90% for big output files). diff --git a/deal.II/include/deal.II/numerics/vectors.templates.h b/deal.II/include/deal.II/numerics/vectors.templates.h index 5e74334edf..cc2abca4bd 100644 --- a/deal.II/include/deal.II/numerics/vectors.templates.h +++ b/deal.II/include/deal.II/numerics/vectors.templates.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name: $ // -// Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors +// Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010, 2011 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -4194,272 +4194,280 @@ namespace internal typename DH::active_cell_iterator cell = dof.begin_active(), endc = dof.end(); for (unsigned int index=0; cell != endc; ++cell, ++index) - { - double diff=0; - // initialize for this cell - x_fe_values.reinit (cell); - - const dealii::FEValues &fe_values = x_fe_values.get_present_fe_values (); - const unsigned int n_q_points = fe_values.n_quadrature_points; - - // resize all out scratch - // arrays to the number of - // quadrature points we use - // for the present cell - function_values.resize (n_q_points, - dealii::Vector(n_components)); - function_grads.resize (n_q_points, - std::vector >(n_components)); - - weight_values.resize (n_q_points); - weight_vectors.resize (n_q_points, - dealii::Vector(n_components)); - - psi_values.resize (n_q_points, - dealii::Vector(n_components)); - psi_grads.resize (n_q_points, - std::vector >(n_components)); - psi_scalar.resize (n_q_points); - - tmp_values.resize (n_q_points); - tmp_gradients.resize (n_q_points); - - if (weight!=0) - { - if (weight->n_components>1) - weight->vector_value_list (fe_values.get_quadrature_points(), - weight_vectors); - else - { - weight->value_list (fe_values.get_quadrature_points(), - weight_values); - for (unsigned int k=0;kis_artificial() && !cell->is_ghost()) + { + double diff=0; + // initialize for this cell + x_fe_values.reinit (cell); + + const dealii::FEValues &fe_values = x_fe_values.get_present_fe_values (); + const unsigned int n_q_points = fe_values.n_quadrature_points; + + // resize all out scratch + // arrays to the number of + // quadrature points we use + // for the present cell + function_values.resize (n_q_points, + dealii::Vector(n_components)); + function_grads.resize (n_q_points, + std::vector >(n_components)); + + weight_values.resize (n_q_points); + weight_vectors.resize (n_q_points, + dealii::Vector(n_components)); + + psi_values.resize (n_q_points, + dealii::Vector(n_components)); + psi_grads.resize (n_q_points, + std::vector >(n_components)); + psi_scalar.resize (n_q_points); + + tmp_values.resize (n_q_points); + tmp_gradients.resize (n_q_points); + + if (weight!=0) + { + if (weight->n_components>1) + weight->vector_value_list (fe_values.get_quadrature_points(), + weight_vectors); + else + { + weight->value_list (fe_values.get_quadrature_points(), + weight_values); + for (unsigned int k=0;k &mapping, UpdateFlags(update_JxW_values | update_values)); - typename DoFHandler::active_cell_iterator c; + typename DoFHandler::active_cell_iterator cell; std::vector > values(quadrature.size(), Vector (dof.get_fe().n_components())); double mean = 0.; double area = 0.; // Compute mean value - for (c = dof.begin_active(); c != dof.end(); ++c) + for (cell = dof.begin_active(); cell != dof.end(); ++cell) + if (!cell->is_artificial() && !cell->is_ghost()) + { + fe.reinit (cell); + fe.get_function_values(v, values); + for (unsigned int k=0; k< quadrature.size(); ++k) + { + mean += fe.JxW(k) * values[k](component); + area += fe.JxW(k); + } + } + +#if DEAL_II_USE_P4EST + // if this was a distributed + // DoFHandler, we need to do the + // reduction over the entire domain + if (const parallel::distributed::Triangulation * + p_d_triangulation + = dynamic_cast *>(&dof.get_tria())) { - fe.reinit (c); - fe.get_function_values(v, values); - for (unsigned int k=0; k< quadrature.size(); ++k) - { - mean += fe.JxW(k) * values[k](component); - area += fe.JxW(k); - }; - }; + double my_values[2] = { mean, area }; + double global_values[2]; + + MPI_Reduce (&my_values, &global_values, 2, MPI_DOUBLE, + MPI_SUM, 0, + p_d_triangulation->get_communicator()); + + mean = global_values[0]; + area = global_values[1]; + } +#endif return (mean/area); } -- 2.39.5