]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement GridTools::volume.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 3 Sep 2011 00:23:29 +0000 (00:23 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 3 Sep 2011 00:23:29 +0000 (00:23 +0000)
git-svn-id: https://svn.dealii.org/trunk@24245 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/grid/grid_tools.h
deal.II/source/grid/grid_tools.cc
deal.II/source/grid/grid_tools.inst.in
tests/deal.II/grid_tools_03.cc [new file with mode: 0644]
tests/deal.II/grid_tools_03/cmp/generic [new file with mode: 0644]
tests/deal.II/grid_tools_04.cc [new file with mode: 0644]
tests/deal.II/grid_tools_04/cmp/generic [new file with mode: 0644]

index 7e18a31931e0d6e25fe1f530221ca490372565cb..c23460535577073ff112e07c2edc7d25e078df58 100644 (file)
@@ -283,6 +283,11 @@ and DoF handlers embedded in higher dimensional space have been added.
 <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>
index c57c608c9c7b13b1cc1a7e87e269eb8e26e9c749..3e32d1f13242da4f1ed4f7309ed8bf34afa0c6b4 100644 (file)
@@ -18,6 +18,7 @@
 #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>
 
@@ -42,7 +43,6 @@ class SparsityPattern;
  * individual functions for more information.
  *
  * @ingroup grid
- * @author Wolfgang Bangerth, 2001, 2003, 2004, Ralf Hartmann, 2005
  */
 class GridTools
 {
@@ -62,6 +62,21 @@ 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
index 08131de96a51f06f4e23aa66277854ef6816d7e1..79ea68df191cfe65df18a5817bc64e2d66aea90e 100644 (file)
@@ -12,7 +12,7 @@
 //---------------------------------------------------------------------------
 
 
-
+#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>
@@ -27,6 +27,8 @@
 #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>
@@ -78,6 +80,10 @@ template <int dim, int spacedim>
 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
@@ -126,6 +132,58 @@ GridTools::diameter (const Triangulation<dim, spacedim> &tria)
 
 
 
+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,
index 9d0533cd67e5e6eae22ecda108d74b778e281314..a6dcc68690eb1a9d594053af244861d2dd966364 100644 (file)
@@ -62,6 +62,18 @@ for (deal_II_dimension : DIMENSIONS)
       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> > &,
diff --git a/tests/deal.II/grid_tools_03.cc b/tests/deal.II/grid_tools_03.cc
new file mode 100644 (file)
index 0000000..0659ae4
--- /dev/null
@@ -0,0 +1,94 @@
+//----------------------------  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;
+}
+
diff --git a/tests/deal.II/grid_tools_03/cmp/generic b/tests/deal.II/grid_tools_03/cmp/generic
new file mode 100644 (file)
index 0000000..f785507
--- /dev/null
@@ -0,0 +1,17 @@
+
+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
diff --git a/tests/deal.II/grid_tools_04.cc b/tests/deal.II/grid_tools_04.cc
new file mode 100644 (file)
index 0000000..3ea6d46
--- /dev/null
@@ -0,0 +1,68 @@
+//----------------------------  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;
+}
+
diff --git a/tests/deal.II/grid_tools_04/cmp/generic b/tests/deal.II/grid_tools_04/cmp/generic
new file mode 100644 (file)
index 0000000..3df88e4
--- /dev/null
@@ -0,0 +1,5 @@
+
+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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.