]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement projection onto quads.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 31 Jul 2009 04:17:33 +0000 (04:17 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 31 Jul 2009 04:17:33 +0000 (04:17 +0000)
git-svn-id: https://svn.dealii.org/trunk@19145 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/grid/tria_boundary.cc
tests/deal.II/project_to_surface_02.cc [new file with mode: 0644]
tests/deal.II/project_to_surface_02/cmp/generic [new file with mode: 0644]

index 661397012e623563035b0180d596f4f3955e8ca6..8be8dfdf1a67589509f2707135595b559f1ce92d 100644 (file)
@@ -11,6 +11,7 @@
 //
 //---------------------------------------------------------------------------
 
+#include <iostream>
 
 #include <base/tensor.h>
 #include <grid/tria_boundary.h>
@@ -536,19 +537,140 @@ project_to_surface (const typename Triangulation<dim, spacedim>::line_iterator &
 
 
 
+namespace internal
+{
+  template <typename Iterator, int spacedim, int dim>
+  Point<spacedim>
+  compute_projection (const Iterator        &object,
+                     const Point<spacedim> &y,
+                     internal::int2type<dim>)
+  {
+                                    // let's look at this for
+                                    // simplicity for a quad (dim==2)
+                                    // in a space with spacedim>2:
+
+                                    // all points on the surface are given by
+                                    //   x(\xi) = sum_i v_i phi_x(\xi)
+                                    // where v_i are the vertices of the quad,
+                                    // and \xi=(\xi_1,\xi_2) are the reference
+                                    // coordinates of the quad. so what we are
+                                    // trying to do is find a point x on
+                                    // the surface that is closest to the point
+                                    // y. there are different ways
+                                    // to solve this problem, but in the end
+                                    // it's a nonlinear problem and we have to
+                                    // find reference coordinates \xi so that
+                                    //   J(\xi) = 1/2 || x(\xi)-y ||^2
+                                    // is minimal. x(\xi) is a function that
+                                    // is dim-linear in \xi, so J(\xi) is
+                                    // a polynomial of degree 2*dim that
+                                    // we'd like to minimize. unless dim==1,
+                                    // we'll have to use a Newton
+                                    // method to find the
+                                    // answer. This leads to the
+                                    // following formulation of
+                                    // Newton steps:
+                                    //
+                                    // Given \xi_k, find \delta\xi_k so that
+                                    //   H_k \delta\xi_k = - F_k
+                                    // where H_k is an approximation to the
+                                    // second derivatives of J at \xi_k, and
+                                    // F_k is the first derivative of J.
+                                    // We'll iterate this a number of times
+                                    // until the right hand side is small
+                                    // enough. As a stopping criterion, we
+                                    // terminate if ||\delta\xi||<eps.
+                                    //
+                                    // As for the Hessian, the best choice
+                                    // would be
+                                    //   H_k = J''(\xi_k)
+                                    // but we'll opt for the simpler
+                                    // Gauss-Newton form
+                                    //   H_k = A^T A
+                                    // i.e.
+                                    //   (H_k)_{nm} = \sum_{i,j} v_i*v_j *
+                                    //                   \partial_n phi_i *
+                                    //                   \partial_m phi_j
+                                    // we start at xi=(0.5,0.5).
+    Point<dim> xi;
+    for (unsigned int d=0; d<dim; ++d)
+      xi[d] = 0.5;
+
+    Point<spacedim> x_k;
+    for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+      x_k += object->vertex(i) *
+            GeometryInfo<dim>::d_linear_shape_function (xi, i);
+
+    do 
+      {
+       Tensor<1,dim> F_k;
+       for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+         F_k += (x_k-y)*object->vertex(i) *
+               GeometryInfo<dim>::d_linear_shape_function_gradient (xi, i);
+
+       Tensor<2,dim> H_k;
+       for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+         for (unsigned int j=0; j<GeometryInfo<dim>::vertices_per_cell; ++j)
+           {
+             Tensor<2,dim> tmp;
+             outer_product (tmp,
+                            GeometryInfo<dim>::d_linear_shape_function_gradient (xi, i),
+                            GeometryInfo<dim>::d_linear_shape_function_gradient (xi, j));
+             H_k += object->vertex(i) * object->vertex(j) * tmp;
+           }
+
+       const Point<dim> delta_xi = - invert(H_k) * F_k;
+       xi += delta_xi;
+
+       Point<spacedim> old_xk = x_k;
+       x_k = Point<spacedim>();
+       for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+         x_k += object->vertex(i) *
+                GeometryInfo<dim>::d_linear_shape_function (xi, i);
+
+       if (delta_xi.norm() < 1e-5)
+         break;
+      }
+    while (true);
+
+    return x_k;
+  }
+
+
+                                  // specialization for a quad in 1d
+  template <typename Iterator>
+  Point<1>
+  compute_projection (const Iterator &,
+                     const Point<1> &y,
+                     /* it's a quad: */internal::int2type<2>)
+  {
+    return y;
+  }
+
+                                  // specialization for a quad in 2d
+  template <typename Iterator>
+  Point<2>
+  compute_projection (const Iterator &,
+                     const Point<2> &y,
+                     /* it's a quad: */internal::int2type<2>)
+  {
+    return y;
+  }
+}
+
+
+
 template <int dim, int spacedim>
 Point<spacedim>
 StraightBoundary<dim, spacedim>::
-project_to_surface (const typename Triangulation<dim, spacedim>::quad_iterator &,
-                   const Point<spacedim>                                &trial_point) const
+project_to_surface (const typename Triangulation<dim, spacedim>::quad_iterator &quad,
+                   const Point<spacedim>  &y) const
 {
   if (spacedim <= 2)
-    return trial_point;
+    return y;
   else
-    {
-      Assert (false, ExcNotImplemented());
-      return Point<spacedim>();
-    }
+    return internal::compute_projection (quad, y,
+                                        /* it's a quad */internal::int2type<2>());
 }
 
 
@@ -563,6 +685,12 @@ project_to_surface (const typename Triangulation<dim, spacedim>::hex_iterator &,
     return trial_point;
   else
     {
+                                      // we can presumably call the
+                                      // same function as above (it's
+                                      // written in a generic way)
+                                      // but someone needs to check
+                                      // whether that actually yields
+                                      // the correct result
       Assert (false, ExcNotImplemented());
       return Point<spacedim>();
     }
diff --git a/tests/deal.II/project_to_surface_02.cc b/tests/deal.II/project_to_surface_02.cc
new file mode 100644 (file)
index 0000000..4dd1f49
--- /dev/null
@@ -0,0 +1,143 @@
+//----------------------------  project_to_surface_02.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_02.cc  ---------------------------
+
+
+// test StraightBoundary::project_to_surface for quads
+
+
+
+#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>::quads_per_cell; ++e)
+       {
+         const typename Triangulation<dim>::quad_iterator
+           quad = (dim > 2 ? cell->quad(e) :
+                   *reinterpret_cast<const typename Triangulation<dim>::quad_iterator*>(&cell));
+         
+         deallog << "    Quad " << e << ", projected point=";
+
+         const Point<dim> p = boundary.project_to_surface (quad, trial_point);
+         deallog << p;
+         deallog << "  (quad is from ";
+         deallog << quad->vertex(0);
+         deallog << " to ";
+         deallog << quad->vertex(1);
+         deallog << " to ";
+         deallog << quad->vertex(2);
+         deallog << " to ";
+         deallog << quad->vertex(3);
+         deallog << ")" << std::endl;
+
+                                          // now make sure that p is
+                                          // indeed closer to
+                                          // trial_point than any of
+                                          // the vertices of the quad
+         for (unsigned int v=0; v<4; ++v)
+           Assert (p.distance (trial_point) <=
+                   quad->vertex(v).distance (trial_point),
+                   ExcInternalError());
+       }
+      tria.clear();
+    }
+}
+
+
+
+
+int main ()
+{
+  std::ofstream logfile ("project_to_surface_02/output");
+  deallog << std::setprecision (3);
+  deallog << std::fixed;  
+  deallog.attach(logfile);
+  deallog.depth_console (0);
+
+  test<2>();
+  test<3>();
+}
diff --git a/tests/deal.II/project_to_surface_02/cmp/generic b/tests/deal.II/project_to_surface_02/cmp/generic
new file mode 100644 (file)
index 0000000..bf745b9
--- /dev/null
@@ -0,0 +1,21 @@
+
+DEAL::dim=2
+DEAL::  Case 0
+DEAL::    Quad 0, projected point=1.500 1.500  (quad is from 1.000 1.000 to 3.000 1.000 to 1.000 3.000 to 3.000 3.000)
+DEAL::  Case 1
+DEAL::    Quad 0, projected point=1.500 1.500  (quad is from 0.000 1.414 to 1.414 2.828 to -1.414 2.828 to 0.000 4.243)
+DEAL::dim=3
+DEAL::  Case 0
+DEAL::    Quad 0, projected point=1.000 1.500 1.500  (quad is from 1.000 1.000 1.000 to 1.000 3.000 1.000 to 1.000 1.000 3.000 to 1.000 3.000 3.000)
+DEAL::    Quad 1, projected point=3.000 1.500 1.500  (quad is from 3.000 1.000 1.000 to 3.000 3.000 1.000 to 3.000 1.000 3.000 to 3.000 3.000 3.000)
+DEAL::    Quad 2, projected point=1.500 1.000 1.500  (quad is from 1.000 1.000 1.000 to 1.000 1.000 3.000 to 3.000 1.000 1.000 to 3.000 1.000 3.000)
+DEAL::    Quad 3, projected point=1.500 3.000 1.500  (quad is from 1.000 3.000 1.000 to 1.000 3.000 3.000 to 3.000 3.000 1.000 to 3.000 3.000 3.000)
+DEAL::    Quad 4, projected point=1.500 1.500 1.000  (quad is from 1.000 1.000 1.000 to 3.000 1.000 1.000 to 1.000 3.000 1.000 to 3.000 3.000 1.000)
+DEAL::    Quad 5, projected point=1.500 1.500 3.000  (quad is from 1.000 1.000 3.000 to 3.000 1.000 3.000 to 1.000 3.000 3.000 to 3.000 3.000 3.000)
+DEAL::  Case 1
+DEAL::    Quad 0, projected point=0.707 0.707 1.500  (quad is from 0.000 1.414 1.000 to -1.414 2.828 1.000 to 0.000 1.414 3.000 to -1.414 2.828 3.000)
+DEAL::    Quad 1, projected point=2.121 2.121 1.500  (quad is from 1.414 2.828 1.000 to 0.000 4.243 1.000 to 1.414 2.828 3.000 to 0.000 4.243 3.000)
+DEAL::    Quad 2, projected point=0.793 2.207 1.500  (quad is from 0.000 1.414 1.000 to 0.000 1.414 3.000 to 1.414 2.828 1.000 to 1.414 2.828 3.000)
+DEAL::    Quad 3, projected point=-0.621 3.621 1.500  (quad is from -1.414 2.828 1.000 to -1.414 2.828 3.000 to 0.000 4.243 1.000 to 0.000 4.243 3.000)
+DEAL::    Quad 4, projected point=1.500 1.500 1.000  (quad is from 0.000 1.414 1.000 to 1.414 2.828 1.000 to -1.414 2.828 1.000 to 0.000 4.243 1.000)
+DEAL::    Quad 5, projected point=1.500 1.500 3.000  (quad is from 0.000 1.414 3.000 to 1.414 2.828 3.000 to -1.414 2.828 3.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.