First working PolarManifold

git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@30391 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/examples/step-1/step-1.cc b/deal.II/examples/step-1/step-1.cc
index eb44c75af1bd009ffe01aff4fac191e6a3063ea8..42def746f606e087cff426dd6cba3135fe984bd0 100644 (file)
--- a/deal.II/examples/step-1/step-1.cc
+++ b/deal.II/examples/step-1/step-1.cc
@@ -38,6 +38,8 @@
 // And this for the declarations of the `sqrt' and `fabs' functions:
 #include <cmath>
 
+#include <deal.II/../../source/grid/manifold_lib.cc>
+
 // The final step in importing deal.II is this: All deal.II functions and
 // classes are in a namespace <code>dealii</code>, to make sure they don't
 // clash with symbols from other libraries you may want to use in conjunction
@@ -100,6 +102,20 @@ void second_grid ()
   GridGenerator::hyper_shell (triangulation,
                               center, inner_radius, outer_radius,
                               10);
+
+   Triangulation<2>::active_cell_iterator
+      cell = triangulation.begin_active(),
+      endc = triangulation.end();
+   for(; cell != endc; ++cell)
+     cell->set_all_manifold_ids(0);
+   
+       
+     // for(unsigned int f=0; f<GeometryInfo<2>::faces_per_cell; ++f)
+     //   if(cell->face(f)->at_boundary())
+     //         cell->face(f)->set_manifold_id(0);
+   
+       
+  
   // By default, the triangulation assumes that all boundaries are straight
   // and given by the cells of the coarse grid (which we just created). It
   // uses this information when cells at the boundary are refined and new
@@ -129,6 +145,9 @@ void second_grid ()
   // In order to demonstrate how to write a loop over all cells, we will
   // refine the grid in five steps towards the inner circle of the domain:
   for (unsigned int step=0; step<5; ++step)
+  // triangulation.refine_global(1);
+  
+  // if(false)
     {
       // Next, we need an iterator which points to a cell and which we will
       // move over all active cells one by one (active cells are those that
@@ -208,6 +227,26 @@ void second_grid ()
   // fine.
 }
 
+void test_case()
+{
+  const unsigned int dim = 2;
+  
+  std::vector<Point<dim> > ps(2);
+  ps[0] = Point<dim>(1.0,1.0);
+  ps[1] = Point<dim>(-1.0,1.0);
+  
+  std::vector<double> ws(2, .5);
+
+  PolarManifold<dim> manifold;
+
+  Point<dim> p = manifold.get_new_point(ps, ws);
+  std::cout << "Point 0 : " << ps[0]
+           << ", Point 1: " << ps[1]
+           << ", computed point: " << p << std::endl;
+  
+  
+}
+
 
 
 // @sect3{The main function}
@@ -216,6 +255,7 @@ void second_grid ()
 // two subfunctions, which produce the two grids.
 int main ()
 {
-  first_grid ();
+  // first_grid ();
   second_grid ();
+  test_case();
 }

diff --git a/deal.II/include/deal.II/grid/manifold.h b/deal.II/include/deal.II/grid/manifold.h
index 922d7d1cab1cf23666c9d59e614884f8d415b5e0..77941d28cd1970899b6970326c5ef4f61e2ec1fe 100644 (file)
--- a/deal.II/include/deal.II/grid/manifold.h
+++ b/deal.II/include/deal.II/grid/manifold.h
@@ -117,8 +117,7 @@ public:
    * of type ExcPureFunctionCalled.
    */
   virtual
-  Point<spacedim>
-  project_to_manifold (const Point<spacedim> &candidate) const;
+  void project_to_manifold (Point<spacedim> &candidate) const;
 };
 
 
@@ -157,6 +156,12 @@ public:
     virtual Point<spacedim>
     get_new_point(const std::vector<Point<spacedim> > &surrounding_points,
                  const std::vector<double> &weights) const;
+
+                                    /**
+                                       Project to FlatManifold: do nothing.
+                                     */
+  virtual
+  void project_to_manifold (Point<spacedim> &candidate) const;
 };
 
 

diff --git a/deal.II/source/grid/manifold.cc b/deal.II/source/grid/manifold.cc
index 07a76418aadf4ad4a907b1dcaa932f8e1bd77a19..10e3a51409e192b2954197d821f844dfbec43512 100644 (file)
--- a/deal.II/source/grid/manifold.cc
+++ b/deal.II/source/grid/manifold.cc
@@ -33,11 +33,10 @@ Manifold<spacedim>::~Manifold ()
 
 
 template <int spacedim>
-Point<spacedim>
-Manifold<spacedim>::project_to_manifold (const Point<spacedim> &candidate) const 
+void
+Manifold<spacedim>::project_to_manifold (Point<spacedim> &candidate) const 
 {
   Assert (false, ExcPureFunctionCalled());
-  return Point<spacedim>();
 }
 
 
@@ -57,14 +56,24 @@ get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
 {
   Assert(surrounding_points.size() == weights.size(),
         ExcDimensionMismatch(surrounding_points.size(), weights.size()));
+
+#ifdef DEBUG
+  double sum=std::accumulate(weights.begin(), weights.end(), 0.0);
+  Assert(std::abs(sum-1.0) < 1e-10, ExcMessage("Weights should sum to 1!"));
+#endif
+  
   
   Point<spacedim> p;
   for(unsigned int i=0; i<surrounding_points.size(); ++i)
     p += surrounding_points[i]*weights[i];
-  
   return p;
 }
 
+template <int spacedim>
+void
+FlatManifold<spacedim>::project_to_manifold (Point<spacedim> &candidate) const 
+{}
+
 
 /* -------------------------- ManifoldChart --------------------- */
 

diff --git a/deal.II/source/grid/manifold_lib.cc b/deal.II/source/grid/manifold_lib.cc
index c7dae42946f3943c82c33c6eb23074452bf3872e..5df357c341b376fe28880d96370f51d44b78c1bb 100644 (file)
--- a/deal.II/source/grid/manifold_lib.cc
+++ b/deal.II/source/grid/manifold_lib.cc
@@ -41,12 +41,12 @@ get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
 
   double radius = 0;
   for(unsigned int i=0; i<surrounding_points.size(); ++i)
-    radius += weights[i]*surrounding_points[i].norm();
+    radius += weights[i]*(surrounding_points[i]-center).norm();
 
   Point<spacedim> p = FlatManifold<spacedim>::get_new_point(surrounding_points, weights);
-
-  p = p*(radius/p.norm());
-  return p;
+  p = p-center;
+  p = p/p.norm()*radius+center;  
+  return p; 
 }
 
 

diff --git a/deal.II/source/grid/manifold_lib.inst.in b/deal.II/source/grid/manifold_lib.inst.in
index 5349d55ec208ea53b2db178d0d3a4a5a5ad446c0..2a5ddfc97891f406f21e2b93f63ca0944e0c77b2 100644 (file)
--- a/deal.II/source/grid/manifold_lib.inst.in
+++ b/deal.II/source/grid/manifold_lib.inst.in
@@ -17,7 +17,7 @@
 
 for (deal_II_dimension : DIMENSIONS)
   {
-    template class PolarManifold<deal_II_dimension>;
+//    template class PolarManifold<deal_II_dimension>;
 //     template class ConeBoundary<deal_II_dimension>;
 //     template class HyperBallBoundary<deal_II_dimension>;
 //     template class HalfHyperBallBoundary<deal_II_dimension>;