//
//---------------------------------------------------------------------------
+#include <iostream>
#include <base/tensor.h>
#include <grid/tria_boundary.h>
+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>());
}
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>();
}
--- /dev/null
+//---------------------------- 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>();
+}
--- /dev/null
+
+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)
projects / dealii.git / commitdiff