From 3bb45d7bd2470386597ed82c87c3d5447dd6f6cc Mon Sep 17 00:00:00 2001 From: heltai Date: Wed, 29 Jan 2014 11:07:50 +0000 Subject: [PATCH] Restored projection to StraightBoundaries git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@32342 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/grid/tria_boundary.h | 9 ++ deal.II/source/grid/tria_boundary.cc | 95 ++++++++++++++++++++ 2 files changed, 104 insertions(+) diff --git a/deal.II/include/deal.II/grid/tria_boundary.h b/deal.II/include/deal.II/grid/tria_boundary.h index 882304fa53..9367c8f3e7 100644 --- a/deal.II/include/deal.II/grid/tria_boundary.h +++ b/deal.II/include/deal.II/grid/tria_boundary.h @@ -164,6 +164,15 @@ public: Point normal_vector(const std::vector > &vertices, const Point &point) const; + /** + * Given a point and a vector of vertices, project the point + * back to the flat boundary. + */ + virtual + Point + project_to_manifold(const std::vector > &vertices, + const Point &point) const; + }; diff --git a/deal.II/source/grid/tria_boundary.cc b/deal.II/source/grid/tria_boundary.cc index d33e1ac01f..9bfca6e8be 100644 --- a/deal.II/source/grid/tria_boundary.cc +++ b/deal.II/source/grid/tria_boundary.cc @@ -260,6 +260,101 @@ StraightBoundary::normal_vector(const std::vector return internal::normalized_alternating_product(grad_F); } +template +Point +StraightBoundary::project_to_manifold(const std::vector > &v, + const Point &trial_point) const { + if(spacedim <= 1) + return trial_point; + else + switch(v.size()) { + case 2: + { + // 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 double s = (trial_point-v[0])*(v[1]-v[0]) / ((v[1]-v[0])*(v[1]-v[0])); + return v[0] + s*(v[1]-v[0]); + } + break; + case 4: + { + // 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 += v[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-trial_point)*v[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 += (v[i] * v[j]) * tmp; + } + + const Point delta_xi = - invert(H_k) * F_k; + xi += delta_xi; + + x_k = Point(); + for (unsigned int i=0; i::vertices_per_cell; ++i) + x_k += v[i] * GeometryInfo::d_linear_shape_function (xi, i); + + if (delta_xi.norm() < 1e-5) + break; + } + while (true); + + return x_k; + } + break; + default: + Assert(false, ExcNotImplemented()); + return trial_point; + } +} + // explicit instantiations #include "tria_boundary.inst" -- 2.39.5