From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Mon, 18 Aug 2014 21:47:44 +0000 (-0500)
Subject: Patch by k3daevin: Implement at least compute the area of planar faces of 3d cells.
X-Git-Tag: v8.2.0-rc1~197^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a931a7038a2f6c0fab55428e025af0154f589737;p=dealii.git

Patch by k3daevin: Implement at least compute the area of planar faces of 3d cells.
---

diff --git a/doc/news/changes.h b/doc/news/changes.h
index ec80a4f962..d50891f620 100644
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -250,6 +250,12 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> New: TriaAccessor::measure() is now also implemented for faces of
+  3d cells as long as the face is planar.
+  <br>
+  (Kevin Drzycimski, 2014/08/19)
+  </li>
+
   <li> Fixed: Support SLEPc 3.5 by disabling SDFOLD spectrum transformation type
   that has been removed from SLEPc. Therefore, TransformationSpectrumFolding 
   cannot be used with newer SLEPc versions. 
diff --git a/include/deal.II/grid/tria_accessor.h b/include/deal.II/grid/tria_accessor.h
index ace33fe1d1..6aea7debe7 100644
--- a/include/deal.II/grid/tria_accessor.h
+++ b/include/deal.II/grid/tria_accessor.h
@@ -1558,12 +1558,32 @@ public:
   Point<spacedim> barycenter () const;
 
   /**
-   * Volume of the object.  Here, the
-   * volume is defined to be confined by
-   * the (dim-)linear mapping of the unit
-   * cell.  No information about the actual
-   * geometric boundary of the domain is
-   * used.
+   * Compute the dim-dimensional measure of the object. For a
+   * dim-dimensional cell in dim-dimensional space, this equals its
+   * volume. On the other hand, for a 2d cell in 3d space, or if the
+   * current object pointed to is a 2d face of a 3d cell in 3d space,
+   * then the function computes the area the object occupies. For a
+   * one-dimensional object, return its length.
+   *
+   * The function only computes the measure of cells, faces or edges
+   * assumed to be represented by (bi-/tri-)linear mappings. In other
+   * words, it only takes into account the locations of the vertices
+   * that bound the current object but not how the interior of the
+   * object may actually be mapped. In most simple cases, this is
+   * exactly what you want. However, for objects that are not
+   * "straight", e.g. 2d cells embedded in 3d space as part of a
+   * triangulation of a curved domain, two-dimensional faces of 3d
+   * cells that are not just parallelograms, or for faces that are at
+   * the boundary of a domain that is not just bounded by straight
+   * line segments or planes, this function only computes the
+   * dim-dimensional measure of a (bi-/tri-)linear interpolation of
+   * the "real" object as defined by the manifold or boundary object
+   * describing the real geometry of the object in question. If you
+   * want to consider the "real" geometry, you will need to compute
+   * the measure by integrating a function equal to one over the
+   * object, which after applying quadrature equals the summing the
+   * JxW values returned by the FEValues or FEFaceValues object you
+   * will want to use for the integral.
    */
   double measure () const;
 
diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc
index c9f4cb2455..09008f51f8 100644
--- a/source/grid/tria_accessor.cc
+++ b/source/grid/tria_accessor.cc
@@ -970,7 +970,38 @@ namespace
   }
 
 
+  // a 2d face in 3d space
+  double measure (const dealii::TriaAccessor<2,3,3> &accessor)
+  {
+    // If the face is planar, the diagonal from vertex 0 to vertex 3,
+    // v_03, should be in the plane P_012 of vertices 0, 1 and 2.  Get
+    // the normal vector of P_012 and test if v_03 is orthogonal to
+    // that. If so, the face is planar and computing its area is simple.
+    const Point<3> v01 = accessor.vertex(1) - accessor.vertex(0);
+    const Point<3> v02 = accessor.vertex(2) - accessor.vertex(0);
+
+    Point<3> normal;
+    cross_product(normal, v01, v02);
+
+    const Point<3> v03 = accessor.vertex(3) - accessor.vertex(0); //1st diagonal
+
+    if (std::abs((v03 * normal) / (v03.norm() * v01.norm() * v02.norm())) >= 1e-12)
+      {
+	Assert (false,
+		ExcMessage("Computing the measure of a nonplanar face is not implemented!"));
+	return std::numeric_limits<double>::quiet_NaN();
+      }
+
+    // the face is planar. then its area is 1/2 of the norm of the
+    // cross product of the two diagonals
+    const Point<3> v12 = accessor.vertex(2) - accessor.vertex(1);
+    Point<3> twice_area;
+    cross_product(twice_area, v03, v12);
+    return 0.5 * twice_area.norm();
+  }
+
 
+  
   template <int structdim, int dim, int spacedim>
   double
   measure (const TriaAccessor<structdim, dim, spacedim> &)
@@ -978,7 +1009,7 @@ namespace
     // catch-all for all cases not explicitly
     // listed above
     Assert (false, ExcNotImplemented());
-    return 1./0.;
+    return std::numeric_limits<double>::quiet_NaN();
   }
 }
 
diff --git a/tests/deal.II/measure_of_3d_face_01.cc b/tests/deal.II/measure_of_3d_face_01.cc
new file mode 100644
index 0000000000..d718e2266b
--- /dev/null
+++ b/tests/deal.II/measure_of_3d_face_01.cc
@@ -0,0 +1,119 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Test by Kevin Drzycimski: compute the measure of the faces of a 3d cell
+
+#include "../tests.h"
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <fstream>
+#include <iomanip>
+
+
+// move backward two adjacent vertices of the top face up by one
+// unit. all faces remain flat this way
+Point<3> distort_planar(Point<3> p) 
+{
+  if (p(1) > 0.5 && p(2) > 0.5) 
+    {
+      p(1) += 1;
+    }
+  return p;
+}
+
+
+// lift two opposite vertices of the top face up by one unit to create
+// a saddle surface
+Point<3> distort_twisted(Point<3> p) 
+{
+  if (p(2) > 0.5 && (p(0) > 0.5 ^ p(1) > 0.5)) 
+    {
+      p(2) += 1;
+    }
+  return p;
+}
+
+
+
+void test ()
+{
+  Triangulation<3> tria;
+  GridOut gridout;
+
+  deallog << "Planar\n";
+  GridGenerator::hyper_cube(tria);
+  GridTools::transform(&distort_planar, tria);
+  gridout.write_eps(tria, deallog.get_file_stream());
+
+  double measure_planar[] = {1.5, 1.5, 1, ::sqrt(2), 1, 2};
+  deallog << "Face\tExact\tMeasure" << std::endl;
+  Triangulation<3>::active_cell_iterator cell = tria.begin_active();
+  for (int i = 0; i < 6; ++i)
+    {
+      double m = cell->face(i)->measure();
+      deallog << i << '\t' << measure_planar[i] << '\t' << m << std::endl;
+    }
+ 
+  deallog << "Twisted\n";
+  tria.clear();
+  GridGenerator::hyper_cube(tria);
+  GridTools::transform(&distort_twisted, tria);
+  gridout.write_eps(tria, deallog.get_file_stream());
+
+  double measure_twisted[] = {1.5, 1.5, 1.5, 1.5, 1, 5./3};
+  deallog << "Face\tExact\tMeasure" << std::endl;
+  cell = tria.begin_active();
+  for (int i = 0; i < 6; ++i)
+    {
+      double m;
+      try
+	{
+	  m = cell->face(i)->measure();
+	}
+      catch (...)
+	{
+	  m = std::numeric_limits<double>::quiet_NaN();
+	}
+      deallog << i << '\t' << measure_twisted[i] << '\t' << m << std::endl;
+    }
+
+  deallog << std::endl;
+}
+
+
+
+int main()
+{
+  std::ofstream logfile ("output");
+  deallog << std::setprecision (5);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  // run the tests but continue when finding an exception: we will try
+  // out a distorted cell for which TriaAccessor::measure() will error
+  // out, but we'd still like to print and check the data for the
+  // remaining faces
+  deal_II_exceptions::disable_abort_on_exception();
+  
+  test ();
+}
+
diff --git a/tests/deal.II/measure_of_3d_face_01.output b/tests/deal.II/measure_of_3d_face_01.output
new file mode 100644
index 0000000000..3a784e6301
--- /dev/null
+++ b/tests/deal.II/measure_of_3d_face_01.output
@@ -0,0 +1,68 @@
+
+%!PS-Adobe-2.0 EPSF-1.2
+%%Title: deal.II Output
+%%Creator: the deal.II library
+
+%%BoundingBox: 0 0 301 319
+/m {moveto} bind def
+/x {lineto stroke} bind def
+/b {0 0 0 setrgbcolor} def
+/r {1 0 0 setrgbcolor} def
+%%EndProlog
+
+0.500000 setlinewidth
+b 0.00000 40.1924 m 80.3848 109.808 x
+b 139.230 0.00000 m 219.615 69.6152 x
+b 0.00000 40.1924 m 139.230 0.00000 x
+b 80.3848 109.808 m 219.615 69.6152 x
+b 5.15256e-15 179.423 m 160.770 318.653 x
+b 139.230 139.230 m 300.000 278.461 x
+b 5.15256e-15 179.423 m 139.230 139.230 x
+b 160.770 318.653 m 300.000 278.461 x
+b 0.00000 40.1924 m 5.15256e-15 179.423 x
+b 139.230 0.00000 m 139.230 139.230 x
+b 80.3848 109.808 m 160.770 318.653 x
+b 219.615 69.6152 m 300.000 278.461 x
+showpage
+DEAL::Planar
+Face	Exact	Measure
+DEAL::0	1.5000	1.5000
+DEAL::1	1.5000	1.5000
+DEAL::2	1.0000	1.0000
+DEAL::3	1.4142	1.4142
+DEAL::4	1.0000	1.0000
+DEAL::5	2.0000	2.0000
+%!PS-Adobe-2.0 EPSF-1.2
+%%Title: deal.II Output
+%%Creator: the deal.II library
+
+%%BoundingBox: 0 0 301 531
+/m {moveto} bind def
+/x {lineto stroke} bind def
+/b {0 0 0 setrgbcolor} def
+/r {1 0 0 setrgbcolor} def
+%%EndProlog
+
+0.500000 setlinewidth
+b 0.00000 54.9038 m 109.808 150.000 x
+b 190.192 0.00000 m 300.000 95.0962 x
+b 0.00000 54.9038 m 190.192 0.00000 x
+b 109.808 150.000 m 300.000 95.0962 x
+b 7.03853e-15 245.096 m 109.808 530.385 x
+b 190.192 380.385 m 300.000 285.289 x
+b 7.03853e-15 245.096 m 190.192 380.385 x
+b 109.808 530.385 m 300.000 285.289 x
+b 0.00000 54.9038 m 7.03853e-15 245.096 x
+b 190.192 0.00000 m 190.192 380.385 x
+b 109.808 150.000 m 109.808 530.385 x
+b 300.000 95.0962 m 300.000 285.289 x
+showpage
+DEAL::Twisted
+Face	Exact	Measure
+DEAL::0	1.5000	1.5000
+DEAL::1	1.5000	1.5000
+DEAL::2	1.5000	1.5000
+DEAL::3	1.5000	1.5000
+DEAL::4	1.0000	1.0000
+DEAL::5	1.6667	nan
+DEAL::