]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Patch by k3daevin: Implement at least compute the area of planar faces of 3d cells.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 18 Aug 2014 21:47:44 +0000 (16:47 -0500)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 19 Aug 2014 15:31:37 +0000 (10:31 -0500)
doc/news/changes.h
include/deal.II/grid/tria_accessor.h
source/grid/tria_accessor.cc
tests/deal.II/measure_of_3d_face_01.cc [new file with mode: 0644]
tests/deal.II/measure_of_3d_face_01.output [new file with mode: 0644]

index ec80a4f962eae51aa8b1c508b8ffde769f4246e5..d50891f620706a31c183af7d0a6112c83089dfd6 100644 (file)
@@ -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. 
index ace33fe1d19b5ea84771835470e70a295ddb1632..6aea7debe7df7514defdf5a0c0a4a4635025b635 100644 (file)
@@ -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;
 
index c9f4cb2455a7c3a1476519e524d0361dca33a4a6..09008f51f8a122aef920c095498eaba246e7db35 100644 (file)
@@ -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 (file)
index 0000000..d718e22
--- /dev/null
@@ -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 (file)
index 0000000..3a784e6
--- /dev/null
@@ -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::

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.