#include <grid/tria_accessor.templates.h>
#include <grid/tria_iterator.templates.h>
#include <grid/geometry_info.h>
+#include <fe/mapping_q1.h>
#include <cmath>
template <>
bool CellAccessor<3>::point_inside (const Point<3> &p) const
{
- Assert (false, ExcNotImplemented() );
- return false;
+ // original implementation by Joerg
+ // Weimar
+
+ // we first eliminate points based
+ // on the maximum and minumum of
+ // the corner coordinates, then
+ // transform to the unit cell, and
+ // check there.
+ const unsigned int dim = 3;
+ Point<dim> maxp = this->vertex(0);
+ Point<dim> minp = this->vertex(0);
+
+ for (unsigned int v=1; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+ for (unsigned int d=0; d<dim; ++d)
+ {
+ maxp[d] = std::max (maxp[d],this->vertex(v)[d]);
+ minp[d] = std::min (minp[d],this->vertex(v)[d]);
+ }
+
+ // rule out points outside the
+ // bounding box of this cell
+ for (unsigned int d=0; d<dim; d++)
+ if ((p[d] < minp[d]) || (p[d] > maxp[d]))
+ return false;
+
+ // now we need to check more
+ // carefully: transform to the
+ // unit cube
+ // and check there.
+ static const MappingQ1<dim> mapping;
+
+ const TriaRawIterator<dim, CellAccessor<dim> > cell_iterator (*this);
+ return (GeometryInfo<dim>::
+ is_inside_unit_cell (mapping.transform_real_to_unit_cell(cell_iterator,
+ p)));
}
<h3>deal.II</h3>
<ol>
+ <li> <p> New: The <code class="class">TriaAccessor</code>::<code
+ class="member">point_inside</code> function is now also
+ implemented in 3d.
+ <br>
+ (Joerg Weimar, WB 2002/09/04)
+ </p>
+
<li> <p> New: The <code class="class">TriaAccessor</code>::<code
class="member">recursively_set_material_id</code> function sets
the material id of the present cell and of all its children,
roy_1.exe : roy_1.g.$(OBJEXT) $(libraries)
denis_1.exe : denis_1.g.$(OBJEXT) $(libraries)
point_inside_1.exe : point_inside_1.g.$(OBJEXT) $(libraries)
+point_inside_2.exe : point_inside_2.g.$(OBJEXT) $(libraries)
unit_support_points.exe : unit_support_points.g.$(OBJEXT) $(libraries)
parameter_handler_1.exe : parameter_handler_1.g.$(OBJEXT) $(libraries)
parameter_handler_2.exe : parameter_handler_2.g.$(OBJEXT) $(libraries)
tests = anna_1 anna_2 anna_3 anna_4 anna_5 anna_6 \
- geometry_info_1 point_inside_1 \
+ geometry_info_1 point_inside_1 point_inside_2 \
data_out_01 data_out_02 data_out_faces_01 data_out_faces_02 \
data_out_rotation_01 data_out_rotation_02 data_out_stack_01 \
dof_tools_00a \
--- /dev/null
+//---------------------------- point_inside_1.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2002, 2003 by the deal.II authors and Anna Schneebeli
+//
+// 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.
+//
+//---------------------------- point_inside_1.cc ---------------------------
+
+
+// check TriaAccessor<3>::point_inside for a cell that is aligned with
+// the coordinate axes
+//
+// this program is a modified version of one by Joerg Weimar,
+// TU Braunschweig
+
+#include <base/logstream.h>
+#include <grid/tria.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_generator.h>
+#include <fstream>
+
+
+void check ()
+{
+ Triangulation<3> triangulation;
+
+ GridGenerator::hyper_cube (triangulation);
+
+ // Now get the cell
+ const Triangulation<3>::cell_iterator cell = triangulation.begin();
+
+ double testcoord[11][3] = {{0.5,0.5,0.5},
+ {2,0.5,0.5},
+ {0.5,2,0.5},
+ {0.5,0.5,2},
+ {-2,0.5,0.5},
+ {0.5,-2,0.5},
+ {0.5,0.5,-2},
+ {0.9,0.9,0.9},
+ {1.0,0.5,0.5},
+ {0.9999999,0.5,0.5},
+ {1.0000001,0.5,0.5} };
+
+ int expected[] = {1,0,0,0,0,0,0,1,1,1,0};
+
+ for (int i=0; i<11; i++)
+ {
+ const Point<3> testpoint(testcoord[i][0],
+ testcoord[i][1],
+ testcoord[i][2]);
+ bool res = cell->point_inside(testpoint);
+ deallog << testpoint << " inside " << res <<std::endl;
+ Assert (res == expected[i], ExcInternalError());
+ }
+}
+
+
+int main ()
+{
+ std::ofstream logfile("point_inside_1.output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ check ();
+}
--- /dev/null
+//---------------------------- point_inside_2.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2002, 2003 by the deal.II authors and Anna Schneebeli
+//
+// 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.
+//
+//---------------------------- point_inside_2.cc ---------------------------
+
+
+// check TriaAccessor<3>::point_inside for a cell that is _not_
+// aligned with the coordinate axes
+//
+// this program is a modified version of one by Joerg Weimar,
+// TU Braunschweig
+
+#include <base/logstream.h>
+#include <grid/tria.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_generator.h>
+#include <fstream>
+
+
+void check ()
+{
+ // use a rather complicated way to
+ // generate a new mesh with one
+ // moved vertex -- one could move
+ // it by hand in the original
+ // triangulation, but this is how
+ // the original testcase came in,
+ // and it tests some more stuff, so
+ // why change things
+ Triangulation<3> triangulation;
+
+ // we generate a hyper_cube and
+ // modify one vertex.
+ GridGenerator::hyper_cube (triangulation);
+
+ // Now get the cell
+ Triangulation<3>::cell_iterator cell = triangulation.begin();
+
+ std::vector<CellData<3> > cells (1, CellData<3>());
+ // get the indices
+ for (unsigned int j=0; j<8; ++j) {
+ cells[0].vertices[j] = cell->vertex_index(j);
+ }
+
+ // get the vertex coordinates
+ Point<3> vertices[8];
+ for (unsigned int j=0; j<8; ++j) {
+ vertices[j] = cell->vertex(j);
+ }
+ // modify one vertex.
+ vertices[0] = Point<3>(-1,0,0);
+
+ SubCellData boundary_info;
+ Triangulation<3> triangulation2;
+
+ // create a new triangulation.
+ triangulation2.create_triangulation (
+ std::vector<Point<3> >(&vertices[0],&vertices[8]),
+ cells , boundary_info
+ );
+
+ Triangulation<3>::cell_iterator cell2 = triangulation2.begin();
+
+ // and test it.
+ double testcoord[14][3] = {{0.5,0.5,0.5},
+ {2,0.5,0.5},
+ {0.5,2,0.5},
+ {0.5,0.5,2},
+ {-2,0.5,0.5},
+ {0.5,-2,0.5},
+ {0.5,0.5,-2},
+ {0.9,0.9,0.9},
+ {1.0,0.5,0.5},
+ {0.9999999,0.5,0.5},
+ {1.0000001,0.5,0.5},
+ {-0.1,0.1,0.1},
+ {-0.24,0.5,0.5},
+ {-0.26,0.5,0.5}
+ };
+
+ int expected[] = {1,0,0,0,0,0,0,1,1,1,0,1,1,0};
+ for (int i=0; i<14; i++)
+ {
+ const Point<3> testpoint(testcoord[i][0],
+ testcoord[i][1],
+ testcoord[i][2]);
+ bool res = cell2->point_inside(testpoint);
+ deallog << testpoint << " \t inside " << res
+ << " expected " << expected[i] << std::endl;
+ Assert (res == expected[i], ExcInternalError());
+ }
+}
+
+
+int main ()
+{
+ std::ofstream logfile("point_inside_2.output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ check ();
+}