]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement point_inside in 3d.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 4 Sep 2003 20:14:15 +0000 (20:14 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 4 Sep 2003 20:14:15 +0000 (20:14 +0000)
git-svn-id: https://svn.dealii.org/trunk@7957 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/grid/tria_accessor.cc
deal.II/doc/news/c-4-0.html
tests/bits/Makefile
tests/bits/point_inside_1.cc [new file with mode: 0644]
tests/bits/point_inside_2.cc [new file with mode: 0644]
tests/results/i686-pc-linux-gnu+gcc2.95/bits/point_inside_1.output [new file with mode: 0644]
tests/results/i686-pc-linux-gnu+gcc2.95/bits/point_inside_2.output [new file with mode: 0644]
tests/results/i686-pc-linux-gnu+gcc3.2/bits/point_inside_1.output [new file with mode: 0644]
tests/results/i686-pc-linux-gnu+gcc3.2/bits/point_inside_2.output [new file with mode: 0644]
tests/results/sparc-sun-solaris2.7+gcc2.95/bits/point_inside_1.output [new file with mode: 0644]
tests/results/sparc-sun-solaris2.7+gcc2.95/bits/point_inside_2.output [new file with mode: 0644]

index 84b110f2ab3a607cd9e17731cbfe8d0cdbec20cf..fafd7935df81d43e4f31c3ff8dba361063d65824 100644 (file)
@@ -19,6 +19,7 @@
 #include <grid/tria_accessor.templates.h>
 #include <grid/tria_iterator.templates.h>
 #include <grid/geometry_info.h>
+#include <fe/mapping_q1.h>
 
 #include <cmath>
 
@@ -1832,8 +1833,41 @@ void CellAccessor<3>::recursively_set_material_id (const unsigned char mat_id) c
 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)));
 }
 
 
index 7037cf8cff77a8788b73a8047dae969632542bc8..62a41d7d7e40fdf820734245f5071c5be296e326 100644 (file)
@@ -150,6 +150,13 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK
 <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,
index bf2309fe3c2ee1a61f54c486c4cfed892453c064..14562c8989b863a039e312a355936e8603f0d734 100644 (file)
@@ -90,6 +90,7 @@ gerold_1.exe            : gerold_1.g.$(OBJEXT)             $(libraries)
 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)
@@ -103,7 +104,7 @@ cylinder.exe            : cylinder.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 \
diff --git a/tests/bits/point_inside_1.cc b/tests/bits/point_inside_1.cc
new file mode 100644 (file)
index 0000000..feaf94e
--- /dev/null
@@ -0,0 +1,71 @@
+//----------------------------  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 ();
+}
diff --git a/tests/bits/point_inside_2.cc b/tests/bits/point_inside_2.cc
new file mode 100644 (file)
index 0000000..abf9365
--- /dev/null
@@ -0,0 +1,111 @@
+//----------------------------  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 ();
+}
diff --git a/tests/results/i686-pc-linux-gnu+gcc2.95/bits/point_inside_1.output b/tests/results/i686-pc-linux-gnu+gcc2.95/bits/point_inside_1.output
new file mode 100644 (file)
index 0000000..fba48bc
--- /dev/null
@@ -0,0 +1,12 @@
+
+DEAL::0.500000 0.500000 0.500000 inside 1
+DEAL::2.00000 0.500000 0.500000 inside 0
+DEAL::0.500000 2.00000 0.500000 inside 0
+DEAL::0.500000 0.500000 2.00000 inside 0
+DEAL::-2.00000 0.500000 0.500000 inside 0
+DEAL::0.500000 -2.00000 0.500000 inside 0
+DEAL::0.500000 0.500000 -2.00000 inside 0
+DEAL::0.900000 0.900000 0.900000 inside 1
+DEAL::1.00000 0.500000 0.500000 inside 1
+DEAL::1.000000 0.500000 0.500000 inside 1
+DEAL::1.00000 0.500000 0.500000 inside 0
diff --git a/tests/results/i686-pc-linux-gnu+gcc2.95/bits/point_inside_2.output b/tests/results/i686-pc-linux-gnu+gcc2.95/bits/point_inside_2.output
new file mode 100644 (file)
index 0000000..8facc9b
--- /dev/null
@@ -0,0 +1,15 @@
+
+DEAL::0.500000 0.500000 0.500000        inside 1 expected 1
+DEAL::2.00000 0.500000 0.500000         inside 0 expected 0
+DEAL::0.500000 2.00000 0.500000         inside 0 expected 0
+DEAL::0.500000 0.500000 2.00000         inside 0 expected 0
+DEAL::-2.00000 0.500000 0.500000        inside 0 expected 0
+DEAL::0.500000 -2.00000 0.500000        inside 0 expected 0
+DEAL::0.500000 0.500000 -2.00000        inside 0 expected 0
+DEAL::0.900000 0.900000 0.900000        inside 1 expected 1
+DEAL::1.00000 0.500000 0.500000         inside 1 expected 1
+DEAL::1.000000 0.500000 0.500000        inside 1 expected 1
+DEAL::1.00000 0.500000 0.500000         inside 0 expected 0
+DEAL::-0.100000 0.100000 0.100000       inside 1 expected 1
+DEAL::-0.240000 0.500000 0.500000       inside 1 expected 1
+DEAL::-0.260000 0.500000 0.500000       inside 0 expected 0
diff --git a/tests/results/i686-pc-linux-gnu+gcc3.2/bits/point_inside_1.output b/tests/results/i686-pc-linux-gnu+gcc3.2/bits/point_inside_1.output
new file mode 100644 (file)
index 0000000..fba48bc
--- /dev/null
@@ -0,0 +1,12 @@
+
+DEAL::0.500000 0.500000 0.500000 inside 1
+DEAL::2.00000 0.500000 0.500000 inside 0
+DEAL::0.500000 2.00000 0.500000 inside 0
+DEAL::0.500000 0.500000 2.00000 inside 0
+DEAL::-2.00000 0.500000 0.500000 inside 0
+DEAL::0.500000 -2.00000 0.500000 inside 0
+DEAL::0.500000 0.500000 -2.00000 inside 0
+DEAL::0.900000 0.900000 0.900000 inside 1
+DEAL::1.00000 0.500000 0.500000 inside 1
+DEAL::1.000000 0.500000 0.500000 inside 1
+DEAL::1.00000 0.500000 0.500000 inside 0
diff --git a/tests/results/i686-pc-linux-gnu+gcc3.2/bits/point_inside_2.output b/tests/results/i686-pc-linux-gnu+gcc3.2/bits/point_inside_2.output
new file mode 100644 (file)
index 0000000..8facc9b
--- /dev/null
@@ -0,0 +1,15 @@
+
+DEAL::0.500000 0.500000 0.500000        inside 1 expected 1
+DEAL::2.00000 0.500000 0.500000         inside 0 expected 0
+DEAL::0.500000 2.00000 0.500000         inside 0 expected 0
+DEAL::0.500000 0.500000 2.00000         inside 0 expected 0
+DEAL::-2.00000 0.500000 0.500000        inside 0 expected 0
+DEAL::0.500000 -2.00000 0.500000        inside 0 expected 0
+DEAL::0.500000 0.500000 -2.00000        inside 0 expected 0
+DEAL::0.900000 0.900000 0.900000        inside 1 expected 1
+DEAL::1.00000 0.500000 0.500000         inside 1 expected 1
+DEAL::1.000000 0.500000 0.500000        inside 1 expected 1
+DEAL::1.00000 0.500000 0.500000         inside 0 expected 0
+DEAL::-0.100000 0.100000 0.100000       inside 1 expected 1
+DEAL::-0.240000 0.500000 0.500000       inside 1 expected 1
+DEAL::-0.260000 0.500000 0.500000       inside 0 expected 0
diff --git a/tests/results/sparc-sun-solaris2.7+gcc2.95/bits/point_inside_1.output b/tests/results/sparc-sun-solaris2.7+gcc2.95/bits/point_inside_1.output
new file mode 100644 (file)
index 0000000..fba48bc
--- /dev/null
@@ -0,0 +1,12 @@
+
+DEAL::0.500000 0.500000 0.500000 inside 1
+DEAL::2.00000 0.500000 0.500000 inside 0
+DEAL::0.500000 2.00000 0.500000 inside 0
+DEAL::0.500000 0.500000 2.00000 inside 0
+DEAL::-2.00000 0.500000 0.500000 inside 0
+DEAL::0.500000 -2.00000 0.500000 inside 0
+DEAL::0.500000 0.500000 -2.00000 inside 0
+DEAL::0.900000 0.900000 0.900000 inside 1
+DEAL::1.00000 0.500000 0.500000 inside 1
+DEAL::1.000000 0.500000 0.500000 inside 1
+DEAL::1.00000 0.500000 0.500000 inside 0
diff --git a/tests/results/sparc-sun-solaris2.7+gcc2.95/bits/point_inside_2.output b/tests/results/sparc-sun-solaris2.7+gcc2.95/bits/point_inside_2.output
new file mode 100644 (file)
index 0000000..8facc9b
--- /dev/null
@@ -0,0 +1,15 @@
+
+DEAL::0.500000 0.500000 0.500000        inside 1 expected 1
+DEAL::2.00000 0.500000 0.500000         inside 0 expected 0
+DEAL::0.500000 2.00000 0.500000         inside 0 expected 0
+DEAL::0.500000 0.500000 2.00000         inside 0 expected 0
+DEAL::-2.00000 0.500000 0.500000        inside 0 expected 0
+DEAL::0.500000 -2.00000 0.500000        inside 0 expected 0
+DEAL::0.500000 0.500000 -2.00000        inside 0 expected 0
+DEAL::0.900000 0.900000 0.900000        inside 1 expected 1
+DEAL::1.00000 0.500000 0.500000         inside 1 expected 1
+DEAL::1.000000 0.500000 0.500000        inside 1 expected 1
+DEAL::1.00000 0.500000 0.500000         inside 0 expected 0
+DEAL::-0.100000 0.100000 0.100000       inside 1 expected 1
+DEAL::-0.240000 0.500000 0.500000       inside 1 expected 1
+DEAL::-0.260000 0.500000 0.500000       inside 0 expected 0

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.