* 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
Point<spacedim>
project_to_surface (const typename Triangulation<dim,spacedim>::quad_iterator &quad,
const Point<spacedim> &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<spacedim>
+ project_to_surface (const typename Triangulation<dim,spacedim>::hex_iterator &hex,
+ const Point<spacedim> &candidate) const;
};
* 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
Point<spacedim>
project_to_surface (const typename Triangulation<dim,spacedim>::quad_iterator &quad,
const Point<spacedim> &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<spacedim>
+ project_to_surface (const typename Triangulation<dim,spacedim>::hex_iterator &hex,
+ const Point<spacedim> &candidate) const;
};
+template <int dim, int spacedim>
+Point<spacedim>
+Boundary<dim, spacedim>::
+project_to_surface (const typename Triangulation<dim, spacedim>::hex_iterator &,
+ const Point<spacedim> &trial_point) const
+{
+ if (spacedim <= 3)
+ return trial_point;
+ else
+ {
+ Assert (false, ExcPureFunctionCalled());
+ return Point<spacedim>();
+ }
+}
+
+
/* -------------------------- StraightBoundary --------------------- */
template <int dim, int spacedim>
Point<spacedim>
StraightBoundary<dim, spacedim>::
-project_to_surface (const typename Triangulation<dim, spacedim>::line_iterator &,
+project_to_surface (const typename Triangulation<dim, spacedim>::line_iterator &line,
const Point<spacedim> &trial_point) const
{
if (spacedim <= 1)
return trial_point;
else
{
- Assert (false, ExcPureFunctionCalled());
- return Point<spacedim>();
+ // 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<spacedim> 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);
}
}
return trial_point;
else
{
- Assert (false, ExcPureFunctionCalled());
+ Assert (false, ExcNotImplemented());
+ return Point<spacedim>();
+ }
+}
+
+
+
+template <int dim, int spacedim>
+Point<spacedim>
+StraightBoundary<dim, spacedim>::
+project_to_surface (const typename Triangulation<dim, spacedim>::hex_iterator &,
+ const Point<spacedim> &trial_point) const
+{
+ if (spacedim <= 3)
+ return trial_point;
+ else
+ {
+ Assert (false, ExcNotImplemented());
return Point<spacedim>();
}
}
get_finest_common_cells_* \
project_*[0-9] \
project_boundary_* \
+ project_to_surface_* \
minimal_cell_diameter \
maximal_cell_diameter \
union_triangulation \
--- /dev/null
+//---------------------------- 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 <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>::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>();
+}
--- /dev/null
+
+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)