]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Allow projection onto hex. Add testcase.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 30 Jul 2009 20:12:08 +0000 (20:12 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 30 Jul 2009 20:12:08 +0000 (20:12 +0000)
git-svn-id: https://svn.dealii.org/trunk@19144 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/grid/tria_boundary.h
deal.II/deal.II/source/grid/tria_boundary.cc
tests/deal.II/Makefile
tests/deal.II/project_to_surface_01.cc [new file with mode: 0644]
tests/deal.II/project_to_surface_01/cmp/generic [new file with mode: 0644]

index ef8b5ee517d1fd0522ce21d333a7b561009403e1..1cf7b52a37492d1c9e0a7127ee1a788540c69050 100644 (file)
@@ -324,7 +324,7 @@ class Boundary : public Subscriptor
                                      * characterized by the given
                                      * quad.
                                      *
-                                     * If spacedim==2, then the line
+                                     * If spacedim<=2, then the surface
                                      * represented by the quad
                                      * iterator is the entire space
                                      * (i.e. it is a cell, not a part
@@ -336,6 +336,26 @@ class Boundary : public Subscriptor
     Point<spacedim>
     project_to_surface (const typename Triangulation<dim,spacedim>::quad_iterator &quad,
                        const Point<spacedim> &candidate) const;
+
+                                    /**
+                                     * Same function as above but for
+                                     * a point that is to be
+                                     * projected onto the area
+                                     * characterized by the given
+                                     * quad.
+                                     *
+                                     * If spacedim<=3, then the manifold
+                                     * represented by the hex
+                                     * iterator is the entire space
+                                     * (i.e. it is a cell, not a part
+                                     * of the boundary), and the
+                                     * returned point equals the
+                                     * given input point.
+                                     */
+    virtual
+    Point<spacedim>
+    project_to_surface (const typename Triangulation<dim,spacedim>::hex_iterator &hex,
+                       const Point<spacedim> &candidate) const;
 };
 
 
@@ -491,7 +511,7 @@ class StraightBoundary : public Boundary<dim,spacedim>
                                      * vertices of the given quad
                                      * iterator.
                                      *
-                                     * If spacedim==2, then the line
+                                     * If spacedim<=2, then the surface
                                      * represented by the quad
                                      * iterator is the entire space
                                      * (i.e. it is a cell, not a part
@@ -503,6 +523,33 @@ class StraightBoundary : public Boundary<dim,spacedim>
     Point<spacedim>
     project_to_surface (const typename Triangulation<dim,spacedim>::quad_iterator &quad,
                        const Point<spacedim> &candidate) const;
+
+                                    /**
+                                     * Same function as above but for
+                                     * a point that is to be
+                                     * projected onto the area
+                                     * characterized by the given
+                                     * quad.
+                                     *
+                                     * The point returned is the
+                                     * projection of the candidate
+                                     * point onto the trilinear
+                                     * manifold spanned by the eight
+                                     * vertices of the given hex
+                                     * iterator.
+                                     *
+                                     * If spacedim<=3, then the manifold
+                                     * represented by the hex
+                                     * iterator is the entire space
+                                     * (i.e. it is a cell, not a part
+                                     * of the boundary), and the
+                                     * returned point equals the
+                                     * given input point.
+                                     */
+    virtual
+    Point<spacedim>
+    project_to_surface (const typename Triangulation<dim,spacedim>::hex_iterator &hex,
+                       const Point<spacedim> &candidate) const;
 };
 
 
index dd6ba2ef5d3b8c60fc5c0092ee74084c067b66d2..661397012e623563035b0180d596f4f3955e8ca6 100644 (file)
@@ -197,6 +197,22 @@ project_to_surface (const typename Triangulation<dim, spacedim>::quad_iterator &
 
 
 
+template <int dim, int spacedim>
+Point<spacedim>
+Boundary<dim, spacedim>::
+project_to_surface (const typename Triangulation<dim, spacedim>::hex_iterator &,
+                   const Point<spacedim>                                &trial_point) const
+{
+  if (spacedim <= 3)
+    return trial_point;
+  else
+    {
+      Assert (false, ExcPureFunctionCalled());
+      return Point<spacedim>();
+    }
+}
+
+
 /* -------------------------- StraightBoundary --------------------- */
 
 
@@ -498,15 +514,23 @@ get_normals_at_vertices (const Triangulation<3>::face_iterator &face,
 template <int dim, int spacedim>
 Point<spacedim>
 StraightBoundary<dim, spacedim>::
-project_to_surface (const typename Triangulation<dim, spacedim>::line_iterator &,
+project_to_surface (const typename Triangulation<dim, spacedim>::line_iterator &line,
                    const Point<spacedim>                                &trial_point) const
 {
   if (spacedim <= 1)
     return trial_point;
   else
     {
-      Assert (false, ExcPureFunctionCalled());
-      return Point<spacedim>();
+                                      // find the point that lies on
+                                      // the line p1--p2. the
+                                      // formulas pan out to
+                                      // something rather simple
+                                      // because the mapping to the
+                                      // line is linear
+      const Point<spacedim> p1 = line->vertex(0),
+                           p2 = line->vertex(1);
+      const double s = (trial_point-p1)*(p2-p1) / ((p2-p1)*(p2-p1));
+      return p1 + s*(p2-p1);
     }
 }
 
@@ -522,7 +546,24 @@ project_to_surface (const typename Triangulation<dim, spacedim>::quad_iterator &
     return trial_point;
   else
     {
-      Assert (false, ExcPureFunctionCalled());
+      Assert (false, ExcNotImplemented());
+      return Point<spacedim>();
+    }
+}
+
+
+
+template <int dim, int spacedim>
+Point<spacedim>
+StraightBoundary<dim, spacedim>::
+project_to_surface (const typename Triangulation<dim, spacedim>::hex_iterator &,
+                   const Point<spacedim>                                &trial_point) const
+{
+  if (spacedim <= 3)
+    return trial_point;
+  else
+    {
+      Assert (false, ExcNotImplemented());
       return Point<spacedim>();
     }
 }
index fb49128eee3f61f80616ba815cda33f7e88aceec..091e969279b8a6f288ba8978fa8f9d0cac1787c1 100644 (file)
@@ -68,6 +68,7 @@ tests_x = block_matrices \
        get_finest_common_cells_* \
        project_*[0-9] \
        project_boundary_* \
+       project_to_surface_* \
        minimal_cell_diameter \
        maximal_cell_diameter \
        union_triangulation \
diff --git a/tests/deal.II/project_to_surface_01.cc b/tests/deal.II/project_to_surface_01.cc
new file mode 100644 (file)
index 0000000..9df3928
--- /dev/null
@@ -0,0 +1,136 @@
+//----------------------------  project_to_surface_01.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2005, 2008, 2009 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.
+//
+//----------------------------  project_to_surface_01.cc  ---------------------------
+
+
+// test StraightBoundary::project_to_surface for lines
+
+
+
+#include "../tests.h"
+#include <base/quadrature_lib.h>
+#include <base/logstream.h>
+#include <grid/tria.h>
+#include <grid/tria_boundary.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_tools.h>
+
+#include <fstream>
+#include <iomanip>
+
+
+class Rotate2d
+{
+  public:
+    Rotate2d (const double angle)
+                   :
+                   angle(angle)
+      {}
+
+    template <int spacedim>
+    Point<spacedim> operator() (const Point<spacedim> p) const
+      {
+       Point<spacedim> q;
+       q[0] = std::cos(angle)*p(0) - std::sin(angle) * p(1);
+       q[1] = std::sin(angle)*p(0) + std::cos(angle) * p(1);
+       for (unsigned d=2; d<spacedim; ++d)
+         q[d] = p[d];
+       return q;
+      }
+  private:
+    const double angle;
+};
+
+
+template <int dim>
+void do_rotate (Triangulation<dim> &tria)
+{
+  GridTools::transform (Rotate2d(numbers::PI/4), tria);
+}
+
+
+void do_rotate (Triangulation<1> &)
+{}
+
+
+
+template <int dim>
+void create_triangulation(const bool rotate,
+                         Triangulation<dim> &tria)
+{
+  GridGenerator::hyper_cube(tria, 1., 3.);
+
+  if (rotate)
+    do_rotate (tria);
+}
+
+
+template <int dim>
+void test ()
+{
+  deallog << "dim=" << dim << std::endl;
+  
+  Triangulation<dim> tria;
+  StraightBoundary<dim> boundary;
+
+  for (unsigned int case_no=0; case_no<2; ++case_no)
+    {
+      deallog << "  Case " << case_no << std::endl;
+      create_triangulation((case_no==1), tria);
+
+      const typename Triangulation<dim>::active_cell_iterator cell=tria.begin_active();
+      Point<dim> trial_point;
+      for (unsigned int i=0; i<dim; ++i)
+       trial_point[i] = 1.5;
+      
+      for (unsigned int e=0; e<GeometryInfo<dim>::lines_per_cell; ++e)
+       {
+         deallog << "    Line " << e << ", projected point=";
+         if (dim > 1)
+           deallog << boundary.project_to_surface (cell->line(e), trial_point);
+         else
+           deallog << boundary.project_to_surface (cell, trial_point);
+
+         deallog << "  (line is from ";
+         if (dim > 1)
+           deallog << cell->line(e)->vertex(0);
+         else
+           deallog << cell->vertex(0);
+         deallog << " to ";
+         if (dim > 1)
+           deallog << cell->line(e)->vertex(1);
+         else
+           deallog << cell->vertex(1);
+         
+         deallog << ")" << std::endl;
+       }
+      tria.clear();
+    }
+}
+
+
+
+
+int main ()
+{
+  std::ofstream logfile ("project_to_surface_01/output");
+  deallog << std::setprecision (3);
+  deallog << std::fixed;  
+  deallog.attach(logfile);
+  deallog.depth_console (0);
+
+  test<1>();
+  test<2>();
+  test<3>();
+}
diff --git a/tests/deal.II/project_to_surface_01/cmp/generic b/tests/deal.II/project_to_surface_01/cmp/generic
new file mode 100644 (file)
index 0000000..6d08e77
--- /dev/null
@@ -0,0 +1,44 @@
+
+DEAL::dim=1
+DEAL::  Case 0
+DEAL::    Line 0, projected point=1.500  (line is from 1.000 to 3.000)
+DEAL::  Case 1
+DEAL::    Line 0, projected point=1.500  (line is from 1.000 to 3.000)
+DEAL::dim=2
+DEAL::  Case 0
+DEAL::    Line 0, projected point=1.000 1.500  (line is from 1.000 1.000 to 1.000 3.000)
+DEAL::    Line 1, projected point=3.000 1.500  (line is from 3.000 1.000 to 3.000 3.000)
+DEAL::    Line 2, projected point=1.500 1.000  (line is from 1.000 1.000 to 3.000 1.000)
+DEAL::    Line 3, projected point=1.500 3.000  (line is from 1.000 3.000 to 3.000 3.000)
+DEAL::  Case 1
+DEAL::    Line 0, projected point=0.707 0.707  (line is from 0.000 1.414 to -1.414 2.828)
+DEAL::    Line 1, projected point=2.121 2.121  (line is from 1.414 2.828 to 0.000 4.243)
+DEAL::    Line 2, projected point=0.793 2.207  (line is from 0.000 1.414 to 1.414 2.828)
+DEAL::    Line 3, projected point=-0.621 3.621  (line is from -1.414 2.828 to 0.000 4.243)
+DEAL::dim=3
+DEAL::  Case 0
+DEAL::    Line 0, projected point=1.000 1.500 1.000  (line is from 1.000 1.000 1.000 to 1.000 3.000 1.000)
+DEAL::    Line 1, projected point=3.000 1.500 1.000  (line is from 3.000 1.000 1.000 to 3.000 3.000 1.000)
+DEAL::    Line 2, projected point=1.500 1.000 1.000  (line is from 1.000 1.000 1.000 to 3.000 1.000 1.000)
+DEAL::    Line 3, projected point=1.500 3.000 1.000  (line is from 1.000 3.000 1.000 to 3.000 3.000 1.000)
+DEAL::    Line 4, projected point=1.000 1.500 3.000  (line is from 1.000 1.000 3.000 to 1.000 3.000 3.000)
+DEAL::    Line 5, projected point=3.000 1.500 3.000  (line is from 3.000 1.000 3.000 to 3.000 3.000 3.000)
+DEAL::    Line 6, projected point=1.500 1.000 3.000  (line is from 1.000 1.000 3.000 to 3.000 1.000 3.000)
+DEAL::    Line 7, projected point=1.500 3.000 3.000  (line is from 1.000 3.000 3.000 to 3.000 3.000 3.000)
+DEAL::    Line 8, projected point=1.000 1.000 1.500  (line is from 1.000 1.000 1.000 to 1.000 1.000 3.000)
+DEAL::    Line 9, projected point=3.000 1.000 1.500  (line is from 3.000 1.000 1.000 to 3.000 1.000 3.000)
+DEAL::    Line 10, projected point=1.000 3.000 1.500  (line is from 1.000 3.000 1.000 to 1.000 3.000 3.000)
+DEAL::    Line 11, projected point=3.000 3.000 1.500  (line is from 3.000 3.000 1.000 to 3.000 3.000 3.000)
+DEAL::  Case 1
+DEAL::    Line 0, projected point=0.707 0.707 1.000  (line is from 0.000 1.414 1.000 to -1.414 2.828 1.000)
+DEAL::    Line 1, projected point=2.121 2.121 1.000  (line is from 1.414 2.828 1.000 to 0.000 4.243 1.000)
+DEAL::    Line 2, projected point=0.793 2.207 1.000  (line is from 0.000 1.414 1.000 to 1.414 2.828 1.000)
+DEAL::    Line 3, projected point=-0.621 3.621 1.000  (line is from -1.414 2.828 1.000 to 0.000 4.243 1.000)
+DEAL::    Line 4, projected point=0.707 0.707 3.000  (line is from 0.000 1.414 3.000 to -1.414 2.828 3.000)
+DEAL::    Line 5, projected point=2.121 2.121 3.000  (line is from 1.414 2.828 3.000 to 0.000 4.243 3.000)
+DEAL::    Line 6, projected point=0.793 2.207 3.000  (line is from 0.000 1.414 3.000 to 1.414 2.828 3.000)
+DEAL::    Line 7, projected point=-0.621 3.621 3.000  (line is from -1.414 2.828 3.000 to 0.000 4.243 3.000)
+DEAL::    Line 8, projected point=0.000 1.414 1.500  (line is from 0.000 1.414 1.000 to 0.000 1.414 3.000)
+DEAL::    Line 9, projected point=1.414 2.828 1.500  (line is from 1.414 2.828 1.000 to 1.414 2.828 3.000)
+DEAL::    Line 10, projected point=-1.414 2.828 1.500  (line is from -1.414 2.828 1.000 to -1.414 2.828 3.000)
+DEAL::    Line 11, projected point=0.000 4.243 1.500  (line is from 0.000 4.243 1.000 to 0.000 4.243 3.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.