From: Wolfgang Bangerth Date: Thu, 30 Jul 2009 20:12:08 +0000 (+0000) Subject: Allow projection onto hex. Add testcase. X-Git-Tag: v8.0.0~7421 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3491cf33c55673aa9bbc7e20e928a1ca2e331745;p=dealii.git Allow projection onto hex. Add testcase. git-svn-id: https://svn.dealii.org/trunk@19144 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/grid/tria_boundary.h b/deal.II/deal.II/include/grid/tria_boundary.h index ef8b5ee517..1cf7b52a37 100644 --- a/deal.II/deal.II/include/grid/tria_boundary.h +++ b/deal.II/deal.II/include/grid/tria_boundary.h @@ -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 project_to_surface (const typename Triangulation::quad_iterator &quad, const Point &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 + project_to_surface (const typename Triangulation::hex_iterator &hex, + const Point &candidate) const; }; @@ -491,7 +511,7 @@ class StraightBoundary : public Boundary * 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 Point project_to_surface (const typename Triangulation::quad_iterator &quad, const Point &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 + project_to_surface (const typename Triangulation::hex_iterator &hex, + const Point &candidate) const; }; diff --git a/deal.II/deal.II/source/grid/tria_boundary.cc b/deal.II/deal.II/source/grid/tria_boundary.cc index dd6ba2ef5d..661397012e 100644 --- a/deal.II/deal.II/source/grid/tria_boundary.cc +++ b/deal.II/deal.II/source/grid/tria_boundary.cc @@ -197,6 +197,22 @@ project_to_surface (const typename Triangulation::quad_iterator & +template +Point +Boundary:: +project_to_surface (const typename Triangulation::hex_iterator &, + const Point &trial_point) const +{ + if (spacedim <= 3) + return trial_point; + else + { + Assert (false, ExcPureFunctionCalled()); + return Point(); + } +} + + /* -------------------------- StraightBoundary --------------------- */ @@ -498,15 +514,23 @@ get_normals_at_vertices (const Triangulation<3>::face_iterator &face, template Point StraightBoundary:: -project_to_surface (const typename Triangulation::line_iterator &, +project_to_surface (const typename Triangulation::line_iterator &line, const Point &trial_point) const { if (spacedim <= 1) return trial_point; else { - Assert (false, ExcPureFunctionCalled()); - return Point(); + // 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 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::quad_iterator & return trial_point; else { - Assert (false, ExcPureFunctionCalled()); + Assert (false, ExcNotImplemented()); + return Point(); + } +} + + + +template +Point +StraightBoundary:: +project_to_surface (const typename Triangulation::hex_iterator &, + const Point &trial_point) const +{ + if (spacedim <= 3) + return trial_point; + else + { + Assert (false, ExcNotImplemented()); return Point(); } } diff --git a/tests/deal.II/Makefile b/tests/deal.II/Makefile index fb49128eee..091e969279 100644 --- a/tests/deal.II/Makefile +++ b/tests/deal.II/Makefile @@ -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 index 0000000000..9df39288c9 --- /dev/null +++ b/tests/deal.II/project_to_surface_01.cc @@ -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 +#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::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 index 0000000000..6d08e77c78 --- /dev/null +++ b/tests/deal.II/project_to_surface_01/cmp/generic @@ -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)