<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.
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;
}
+ // 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> &)
// catch-all for all cases not explicitly
// listed above
Assert (false, ExcNotImplemented());
- return 1./0.;
+ return std::numeric_limits<double>::quiet_NaN();
}
}
--- /dev/null
+// ---------------------------------------------------------------------
+// $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 ();
+}
+
--- /dev/null
+
+%!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::
projects / dealii.git / commitdiff