]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement find_active_cell_around_point also for hp cells and mapping collections.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 11 Jan 2007 03:39:42 +0000 (03:39 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 11 Jan 2007 03:39:42 +0000 (03:39 +0000)
git-svn-id: https://svn.dealii.org/trunk@14294 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/grid/grid_tools.h
deal.II/deal.II/source/grid/grid_tools.cc

index 72de802a0bb2089a364a0ffb1e715885979af368..99936e5a2dcc88ea084e1a931c4ac0885648a391 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
+//    Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
 DEAL_II_NAMESPACE_OPEN
 
 
+template <int> class DoFHandler;
+template <int> class Mapping;
+namespace hp
+{
+  template <int> class DoFHandler;
+  template <int> class MappingCollection;
+}
+
+
 /**
  * This class is a collection of algorithms working on triangulations,
  * such as shifting or rotating triangulations, but also finding a
@@ -373,6 +382,23 @@ class GridTools
                                    const Container<dim> &container,
                                    const Point<dim>     &p);
 
+                                    /**
+                                     * A version of the previous function
+                                     * where we use that mapping on a given
+                                     * cell that corresponds to the active
+                                     * finite element index of that
+                                     * cell. This is obviously only useful
+                                     * for hp problems, since the active
+                                     * finite element index for all other DoF
+                                     * handlers is always zero.
+                                     */
+    template <int dim>
+    static
+    std::pair<typename hp::DoFHandler<dim>::active_cell_iterator, Point<dim> >
+    find_active_cell_around_point (const hp::MappingCollection<dim>   &mapping,
+                                   const hp::DoFHandler<dim> &container,
+                                   const Point<dim>     &p);
+    
                                      /**
                                       * Use the METIS partitioner to generate
                                       * a partitioning of the active cells
index fb0171b79e54e1717cc96120d6cfcfe15906754d..b076738097e3f1fa922f6866a207acf252930034 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
+//    Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -25,6 +25,7 @@
 #include <dofs/dof_tools.h>
 #include <fe/fe_dgq.h>
 #include <fe/mapping_q1.h>
+#include <fe/mapping_collection.h>
 #include <multigrid/mg_dof_handler.h>
 #include <multigrid/mg_dof_accessor.h>
 
@@ -685,7 +686,8 @@ typename Container::active_cell_iterator
 GridTools::find_active_cell_around_point (const Container  &container,
                                           const Point<dim> &p)
 {
-   return find_active_cell_around_point(StaticMappingQ1<dim>::mapping, container, p).first;
+   return find_active_cell_around_point(StaticMappingQ1<dim>::mapping,
+                                       container, p).first;
 }
 
 
@@ -718,18 +720,18 @@ GridTools::find_active_cell_around_point (const Mapping<dim>   &mapping,
 
    for(; cell != endc; ++cell)
       {
-         Point<dim> p_cell = mapping.transform_real_to_unit_cell(*cell, p);
+         const Point<dim> p_cell = mapping.transform_real_to_unit_cell(*cell, p);
 
                                    // calculate the infinity norm of
                                    // the distance vector to the unit cell.
-         double dist = GeometryInfo<dim>::distance_to_unit_cell(p_cell);
+         const double dist = GeometryInfo<dim>::distance_to_unit_cell(p_cell);
 
                                    // We compare if the point is inside the
                                    // unit cell (or at least not too far
                                    // outside). If it is, it is also checked
                                    // that the cell has a more refined state
-         if(dist < best_distance ||
-            (dist == best_distance && (*cell)->level() > best_level))
+         if (dist < best_distance ||
+            (dist == best_distance && (*cell)->level() > best_level))
             {
                best_distance = dist;
                best_level    = (*cell)->level();
@@ -737,12 +739,71 @@ GridTools::find_active_cell_around_point (const Mapping<dim>   &mapping,
             }
       }
 
-   Assert(best_cell.first.state() == IteratorState::valid, ExcPointNotFound<dim>(p));
+   Assert (best_cell.first.state() == IteratorState::valid,
+          ExcPointNotFound<dim>(p));
 
    return best_cell;
 }
 
 
+
+template <int dim>
+std::pair<typename hp::DoFHandler<dim>::active_cell_iterator, Point<dim> >
+GridTools::find_active_cell_around_point (const hp::MappingCollection<dim>   &mapping,
+                                          const hp::DoFHandler<dim> &container,
+                                          const Point<dim>     &p)
+{
+   typedef typename hp::DoFHandler<dim>::active_cell_iterator cell_iterator;
+
+                                   // The best distance is set to the
+                                   // maximum allowable distance from
+                                   // the unit cell; we assume a
+                                   // max. deviation of 1e-10
+   double best_distance = 1e-10;
+   int    best_level = -1;
+   std::pair<cell_iterator, Point<dim> > best_cell;
+
+                                   // Find closest vertex and determine
+                                   // all adjacent cells
+   unsigned int vertex = find_closest_vertex(container, p);
+   
+   std::vector<cell_iterator> adjacent_cells =
+      find_cells_adjacent_to_vertex(container, vertex);
+
+   typename std::vector<cell_iterator>::const_iterator
+      cell = adjacent_cells.begin(),
+      endc = adjacent_cells.end();
+
+   for(; cell != endc; ++cell)
+      {
+         const Point<dim> p_cell
+          = mapping[(*cell)->active_fe_index()].transform_real_to_unit_cell(*cell, p);
+
+                                   // calculate the infinity norm of
+                                   // the distance vector to the unit cell.
+         const double dist = GeometryInfo<dim>::distance_to_unit_cell(p_cell);
+
+                                   // We compare if the point is inside the
+                                   // unit cell (or at least not too far
+                                   // outside). If it is, it is also checked
+                                   // that the cell has a more refined state
+         if (dist < best_distance ||
+            (dist == best_distance && (*cell)->level() > best_level))
+            {
+               best_distance = dist;
+               best_level    = (*cell)->level();
+               best_cell     = std::make_pair(*cell, p_cell);
+            }
+      }
+
+   Assert (best_cell.first.state() == IteratorState::valid,
+          ExcPointNotFound<dim>(p));
+
+   return best_cell;
+}
+
+
+
 template <int dim>
 void
 GridTools::
@@ -1188,6 +1249,12 @@ GridTools::find_active_cell_around_point (const Mapping<deal_II_dimension> &,
                                           const hp::DoFHandler<deal_II_dimension> &,
                                           const Point<deal_II_dimension> &);
 
+template
+std::pair<hp::DoFHandler<deal_II_dimension>::active_cell_iterator, Point<deal_II_dimension> >
+GridTools::find_active_cell_around_point (const hp::MappingCollection<deal_II_dimension> &,
+                                          const hp::DoFHandler<deal_II_dimension> &,
+                                          const Point<deal_II_dimension> &);
+
 template
 void
 GridTools::partition_triangulation (const unsigned int,

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.