From 25a985bd3e862dd91d80792cf6f8830b1f4f0c7d Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 31 Jul 2009 04:17:33 +0000 Subject: [PATCH] Implement projection onto quads. git-svn-id: https://svn.dealii.org/trunk@19145 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/grid/tria_boundary.cc | 142 ++++++++++++++++- tests/deal.II/project_to_surface_02.cc | 143 ++++++++++++++++++ .../deal.II/project_to_surface_02/cmp/generic | 21 +++ 3 files changed, 299 insertions(+), 7 deletions(-) create mode 100644 tests/deal.II/project_to_surface_02.cc create mode 100644 tests/deal.II/project_to_surface_02/cmp/generic diff --git a/deal.II/deal.II/source/grid/tria_boundary.cc b/deal.II/deal.II/source/grid/tria_boundary.cc index 661397012e..8be8dfdf1a 100644 --- a/deal.II/deal.II/source/grid/tria_boundary.cc +++ b/deal.II/deal.II/source/grid/tria_boundary.cc @@ -11,6 +11,7 @@ // //--------------------------------------------------------------------------- +#include #include #include @@ -536,19 +537,140 @@ project_to_surface (const typename Triangulation::line_iterator & +namespace internal +{ + template + Point + compute_projection (const Iterator &object, + const Point &y, + internal::int2type) + { + // 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|| xi; + for (unsigned int d=0; d x_k; + for (unsigned int i=0; i::vertices_per_cell; ++i) + x_k += object->vertex(i) * + GeometryInfo::d_linear_shape_function (xi, i); + + do + { + Tensor<1,dim> F_k; + for (unsigned int i=0; i::vertices_per_cell; ++i) + F_k += (x_k-y)*object->vertex(i) * + GeometryInfo::d_linear_shape_function_gradient (xi, i); + + Tensor<2,dim> H_k; + for (unsigned int i=0; i::vertices_per_cell; ++i) + for (unsigned int j=0; j::vertices_per_cell; ++j) + { + Tensor<2,dim> tmp; + outer_product (tmp, + GeometryInfo::d_linear_shape_function_gradient (xi, i), + GeometryInfo::d_linear_shape_function_gradient (xi, j)); + H_k += object->vertex(i) * object->vertex(j) * tmp; + } + + const Point delta_xi = - invert(H_k) * F_k; + xi += delta_xi; + + Point old_xk = x_k; + x_k = Point(); + for (unsigned int i=0; i::vertices_per_cell; ++i) + x_k += object->vertex(i) * + GeometryInfo::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 + 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 + Point<2> + compute_projection (const Iterator &, + const Point<2> &y, + /* it's a quad: */internal::int2type<2>) + { + return y; + } +} + + + template Point StraightBoundary:: -project_to_surface (const typename Triangulation::quad_iterator &, - const Point &trial_point) const +project_to_surface (const typename Triangulation::quad_iterator &quad, + const Point &y) const { if (spacedim <= 2) - return trial_point; + return y; else - { - Assert (false, ExcNotImplemented()); - return Point(); - } + return internal::compute_projection (quad, y, + /* it's a quad */internal::int2type<2>()); } @@ -563,6 +685,12 @@ project_to_surface (const typename Triangulation::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(); } diff --git a/tests/deal.II/project_to_surface_02.cc b/tests/deal.II/project_to_surface_02.cc new file mode 100644 index 0000000000..4dd1f4909f --- /dev/null +++ b/tests/deal.II/project_to_surface_02.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + +class Rotate2d +{ + public: + Rotate2d (const double angle) + : + angle(angle) + {} + + template + Point operator() (const Point p) const + { + Point 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 +void do_rotate (Triangulation &tria) +{ + GridTools::transform (Rotate2d(numbers::PI/4), tria); +} + + +void do_rotate (Triangulation<1> &) +{} + + + +template +void create_triangulation(const bool rotate, + Triangulation &tria) +{ + GridGenerator::hyper_cube(tria, 1., 3.); + + if (rotate) + do_rotate (tria); +} + + +template +void test () +{ + deallog << "dim=" << dim << std::endl; + + Triangulation tria; + StraightBoundary 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::active_cell_iterator cell=tria.begin_active(); + Point trial_point; + for (unsigned int i=0; i::quads_per_cell; ++e) + { + const typename Triangulation::quad_iterator + quad = (dim > 2 ? cell->quad(e) : + *reinterpret_cast::quad_iterator*>(&cell)); + + deallog << " Quad " << e << ", projected point="; + + const Point 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 index 0000000000..bf745b93e5 --- /dev/null +++ b/tests/deal.II/project_to_surface_02/cmp/generic @@ -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) -- 2.39.5