<h3>Specific improvements</h3>
<ol>
+<li> New: There is now a function GridTools::volume computing the volume
+of a triangulation.
+<br>
+(Wolfgang Bangerth, 2011/09/02)
+
<li> New: Code like <code>cell-@>face(1)-@>set_boundary_indicator(42);</code>
now also works in 1d.
<br>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/fe/mapping_q1.h>
#include <list>
* individual functions for more information.
*
* @ingroup grid
- * @author Wolfgang Bangerth, 2001, 2003, 2004, Ralf Hartmann, 2005
*/
class GridTools
{
template <int dim, int spacedim>
static
double diameter (const Triangulation<dim, spacedim> &tria);
+
+ /**
+ * Compute the volume (i.e. the dim-dimensional measure) of the
+ * triangulation. We compute the measure using the integral
+ * $\int 1 \; dx$. The integral approximated is approximated
+ * via quadrature for which we need the mapping argument.
+ *
+ * This function also works for objects of type
+ * parallel::distributed::Triangulation, in which case the
+ * function is a collective operation.
+ */
+ template <int dim, int spacedim>
+ static
+ double volume (const Triangulation<dim,spacedim> &tria,
+ const Mapping<dim,spacedim> &mapping = StaticMappingQ1<dim,spacedim>::mapping);
/**
* Given a list of vertices (typically
//---------------------------------------------------------------------------
-
+#include <deal.II/base/quadrature_lib.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/fe/fe_values.h>
#include <deal.II/hp/mapping_collection.h>
#include <deal.II/multigrid/mg_dof_handler.h>
#include <deal.II/multigrid/mg_dof_accessor.h>
double
GridTools::diameter (const Triangulation<dim, spacedim> &tria)
{
+ Assert ((dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&tria)
+ == 0),
+ ExcNotImplemented());
+
// the algorithm used simply
// traverses all cells and picks
// out the boundary vertices. it
+template <int dim, int spacedim>
+double
+GridTools::volume (const Triangulation<dim, spacedim> &triangulation,
+ const Mapping<dim,spacedim> &mapping)
+{
+ // get the degree of the mapping if possible. if not, just assume 1
+ const unsigned int mapping_degree
+ = (dynamic_cast<const MappingQ<dim,spacedim>*>(&mapping) != 0 ?
+ dynamic_cast<const MappingQ<dim,spacedim>*>(&mapping)->get_degree() :
+ 1);
+
+ // then initialize an appropriate quadrature formula
+ const QGauss<dim> quadrature_formula (mapping_degree + 1);
+ const unsigned int n_q_points = quadrature_formula.size();
+
+ FE_DGQ<dim,spacedim> dummy_fe(0);
+ FEValues<dim,spacedim> fe_values (mapping, dummy_fe, quadrature_formula,
+ update_JxW_values);
+
+ typename Triangulation<dim,spacedim>::active_cell_iterator
+ cell = triangulation.begin_active(),
+ endc = triangulation.end();
+
+ double local_volume = 0;
+
+ // compute the integral quantities by quadrature
+ for (; cell!=endc; ++cell)
+ if (!cell->is_ghost() && !cell->is_artificial())
+ {
+ fe_values.reinit (cell);
+ for (unsigned int q=0; q<n_q_points; ++q)
+ local_volume += fe_values.JxW(q);
+ }
+
+ double global_volume = 0;
+
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+ if (const parallel::distributed::Triangulation<dim,spacedim>* p_tria
+ = dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&triangulation))
+ MPI_Allreduce (&local_volume, &global_volume, 1, MPI_DOUBLE,
+ MPI_SUM,
+ p_tria->get_communicator());
+ else
+ global_volume = local_volume;
+#else
+ global_volume = local_volume;
+#endif
+
+ return global_volume;
+}
+
+
template <>
double
GridTools::cell_measure<3>(const std::vector<Point<3> > &all_vertices,
GridTools::diameter<deal_II_dimension, deal_II_dimension+1> (const Triangulation<deal_II_dimension, deal_II_dimension+1> &);
#endif
+ template
+ double
+ GridTools::volume<deal_II_dimension> (const Triangulation<deal_II_dimension> &,
+ const Mapping<deal_II_dimension> &);
+
+#if deal_II_dimension < 3
+ template
+ double
+ GridTools::volume<deal_II_dimension, deal_II_dimension+1> (const Triangulation<deal_II_dimension, deal_II_dimension+1> &,
+ const Mapping<deal_II_dimension, deal_II_dimension+1> &);
+#endif
+
template
void GridTools::delete_unused_vertices (std::vector<Point<deal_II_dimension> > &,
std::vector<CellData<deal_II_dimension> > &,
--- /dev/null
+//---------------------------- grid_tools.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2008, 2011 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- grid_tools.cc ---------------------------
+
+// check GridTools::volume
+
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_out.h>
+
+#include <fstream>
+#include <iomanip>
+
+std::ofstream logfile("grid_tools_03/output");
+
+
+
+template <int dim>
+void test1 ()
+{
+ // test 1: hypercube
+ if (true)
+ {
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria);
+
+ for (unsigned int i=0; i<2; ++i)
+ {
+ tria.refine_global(2);
+ deallog << dim << "d, "
+ << "hypercube volume, "
+ << i*2
+ << " refinements: "
+ << GridTools::volume (tria)
+ << std::endl;
+ };
+ };
+
+ // test 2: hyperball
+ if (dim >= 2)
+ {
+ Triangulation<dim> tria;
+ GridGenerator::hyper_ball(tria, Point<dim>(), 1);
+
+ static const HyperBallBoundary<dim> boundary;
+ tria.set_boundary (0, boundary);
+
+ for (unsigned int i=0; i<4; ++i)
+ {
+ tria.refine_global(1);
+ deallog << dim << "d, "
+ << "hyperball volume, "
+ << i
+ << " refinements: "
+ << GridTools::volume (tria)
+ << std::endl;
+ }
+ deallog << "exact value="
+ << (dim==2 ? numbers::PI :
+ 4./3.*numbers::PI)
+ << std::endl;
+ }
+}
+
+
+int main ()
+{
+ deallog << std::setprecision(4);
+ logfile << std::setprecision(4);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ test1<1> ();
+ test1<2> ();
+ test1<3> ();
+
+ return 0;
+}
+
--- /dev/null
+
+DEAL::1d, hypercube volume, 0 refinements: 1.000
+DEAL::1d, hypercube volume, 2 refinements: 1.000
+DEAL::2d, hypercube volume, 0 refinements: 1.000
+DEAL::2d, hypercube volume, 2 refinements: 1.000
+DEAL::2d, hyperball volume, 0 refinements: 2.828
+DEAL::2d, hyperball volume, 1 refinements: 3.061
+DEAL::2d, hyperball volume, 2 refinements: 3.121
+DEAL::2d, hyperball volume, 3 refinements: 3.137
+DEAL::exact value=3.142
+DEAL::3d, hypercube volume, 0 refinements: 1.000
+DEAL::3d, hypercube volume, 2 refinements: 1.000
+DEAL::3d, hyperball volume, 0 refinements: 3.210
+DEAL::3d, hyperball volume, 1 refinements: 3.917
+DEAL::3d, hyperball volume, 2 refinements: 4.119
+DEAL::3d, hyperball volume, 3 refinements: 4.171
+DEAL::exact value=4.189
--- /dev/null
+//---------------------------- grid_tools.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2008, 2011 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- grid_tools.cc ---------------------------
+
+// check GridTools::volume for codim-one
+
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_out.h>
+
+#include <fstream>
+#include <iomanip>
+
+std::ofstream logfile("grid_tools_04/output");
+
+
+
+template <int dim>
+void test1 ()
+{
+ // test 1: hypercube
+ if (true)
+ {
+ Triangulation<dim,dim+1> tria;
+ GridGenerator::hyper_cube(tria);
+
+ for (unsigned int i=0; i<2; ++i)
+ {
+ tria.refine_global(2);
+ deallog << dim << "d, "
+ << "hypercube volume, "
+ << i*2
+ << " refinements: "
+ << GridTools::volume (tria)
+ << std::endl;
+ };
+ };
+}
+
+
+int main ()
+{
+ deallog << std::setprecision(4);
+ logfile << std::setprecision(4);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ test1<1> ();
+ test1<2> ();
+
+ return 0;
+}
+
--- /dev/null
+
+DEAL::1d, hypercube volume, 0 refinements: 1.000
+DEAL::1d, hypercube volume, 2 refinements: 1.000
+DEAL::2d, hypercube volume, 0 refinements: 1.000
+DEAL::2d, hypercube volume, 2 refinements: 1.000
projects / dealii-svn.git / commitdiff