From b34603879c1b86d4a3808230bfdfdbc594db3c98 Mon Sep 17 00:00:00 2001 From: heltai Date: Wed, 21 Aug 2013 22:30:46 +0000 Subject: [PATCH] First working PolarManifold git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@30391 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-1/step-1.cc | 42 +++++++++++++++++++++++- deal.II/include/deal.II/grid/manifold.h | 9 +++-- deal.II/source/grid/manifold.cc | 17 +++++++--- deal.II/source/grid/manifold_lib.cc | 8 ++--- deal.II/source/grid/manifold_lib.inst.in | 2 +- 5 files changed, 66 insertions(+), 12 deletions(-) diff --git a/deal.II/examples/step-1/step-1.cc b/deal.II/examples/step-1/step-1.cc index eb44c75af1..42def746f6 100644 --- 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 +#include + // The final step in importing deal.II is this: All deal.II functions and // classes are in a namespace dealii, 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::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 > ps(2); + ps[0] = Point(1.0,1.0); + ps[1] = Point(-1.0,1.0); + + std::vector ws(2, .5); + + PolarManifold manifold; + + Point 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 922d7d1cab..77941d28cd 100644 --- 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 - project_to_manifold (const Point &candidate) const; + void project_to_manifold (Point &candidate) const; }; @@ -157,6 +156,12 @@ public: virtual Point get_new_point(const std::vector > &surrounding_points, const std::vector &weights) const; + + /** + Project to FlatManifold: do nothing. + */ + virtual + void project_to_manifold (Point &candidate) const; }; diff --git a/deal.II/source/grid/manifold.cc b/deal.II/source/grid/manifold.cc index 07a76418aa..10e3a51409 100644 --- a/deal.II/source/grid/manifold.cc +++ b/deal.II/source/grid/manifold.cc @@ -33,11 +33,10 @@ Manifold::~Manifold () template -Point -Manifold::project_to_manifold (const Point &candidate) const +void +Manifold::project_to_manifold (Point &candidate) const { Assert (false, ExcPureFunctionCalled()); - return Point(); } @@ -57,14 +56,24 @@ get_new_point (const std::vector > &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 p; for(unsigned int i=0; i +void +FlatManifold::project_to_manifold (Point &candidate) const +{} + /* -------------------------- ManifoldChart --------------------- */ diff --git a/deal.II/source/grid/manifold_lib.cc b/deal.II/source/grid/manifold_lib.cc index c7dae42946..5df357c341 100644 --- 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 > &surrounding_points, double radius = 0; for(unsigned int i=0; i p = FlatManifold::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 5349d55ec2..2a5ddfc978 100644 --- 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; +// template class PolarManifold; // template class ConeBoundary; // template class HyperBallBoundary; // template class HalfHyperBallBoundary; -- 2.39.5