From: Wolfgang Bangerth 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.

Specific improvements

    +
  1. New: TriaAccessor::measure() is now also implemented for faces of + 3d cells as long as the face is planar. +
    + (Kevin Drzycimski, 2014/08/19) +
    2. +
  2. 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 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::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 double measure (const TriaAccessor &) @@ -978,7 +1009,7 @@ namespace // catch-all for all cases not explicitly // listed above Assert (false, ExcNotImplemented()); - return 1./0.; + return std::numeric_limits::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 +#include +#include +#include +#include +#include +#include +#include + + +// 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::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::