From: heltai Date: Mon, 10 Feb 2014 17:50:30 +0000 (+0000) Subject: Fixed codimension one bug in manifold case. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=08dd8675f61c00269c1a7131f04158fe25b6d639;p=dealii-svn.git Fixed codimension one bug in manifold case. git-svn-id: https://svn.dealii.org/branches/branch_manifold_id_intermediate@32449 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/source/grid/tria.cc b/deal.II/source/grid/tria.cc index ac23fb319f..11188685fa 100644 --- a/deal.II/source/grid/tria.cc +++ b/deal.II/source/grid/tria.cc @@ -4054,8 +4054,13 @@ namespace internal // to the information // stored in the // boundary class + const Boundary &manifold = + (cell->manifold_id() != numbers::invalid_manifold_id ? + cell->get_manifold() : + triangulation.get_boundary(static_cast(cell->material_id()))); + triangulation.vertices[next_unused_vertex] = - cell->get_boundary().get_new_point_on_quad(cell); + manifold.get_new_point_on_quad (cell); } } @@ -4441,14 +4446,20 @@ namespace internal // midpoint; otherwise we // have to ask the manifold // object - if (dim == spacedim) - triangulation.vertices[next_unused_vertex] = - (cell->vertex(0) + cell->vertex(1)) / 2; - else - triangulation.vertices[next_unused_vertex] = - cell->get_boundary().get_new_point_on_line(cell); - triangulation.vertices_used[next_unused_vertex] = true; - + + // Now we always ask the manifold where to put the + // new point. The get_boundary function will return + // a flat boundary if invalid_manifold_id is set, + // behaving as before in the flat case. Backward + // compatibility requires us to use the material id + // of the cell in the codimension one case, however + // we only do this if the manifold_id is the invalid + // one, otherwise use the material id. This is done + // internally in the get_boundary() function. + triangulation.vertices[next_unused_vertex] = + cell->get_boundary().get_new_point_on_line(cell); + triangulation.vertices_used[next_unused_vertex] = true; + // search for next two // unused cell (++ takes // care of the end of the @@ -4865,21 +4876,23 @@ namespace internal // lines we can compute the // midpoint as the mean of // the two vertices: - if (line->at_boundary()) + // if (line->at_boundary()) triangulation.vertices[next_unused_vertex] = line->get_boundary().get_new_point_on_line (line); - else - triangulation.vertices[next_unused_vertex] - = (line->vertex(0) + line->vertex(1)) / 2; } else - // however, if spacedim>dim, we - // always have to ask the - // boundary object for its - // answer - triangulation.vertices[next_unused_vertex] - = triangulation.get_boundary(line->user_index()).get_new_point_on_line (line); - + // however, if spacedim>dim, we always have to ask + // the boundary object for its answer. We use the + // same object of the cell (which was stored in + // line->user_index() before) unless a manifold_id + // has been set on this very line. + if(line->manifold_id() == numbers::invalid_manifold_id) + triangulation.vertices[next_unused_vertex] + = triangulation.get_boundary(line->user_index()).get_new_point_on_line (line); + else + triangulation.vertices[next_unused_vertex] + = line->get_boundary().get_new_point_on_line (line); + // now that we created // the right point, make // up the two child @@ -6014,7 +6027,7 @@ namespace internal // care about boundary quads // here triangulation.vertices[next_unused_vertex] - = (middle_line->vertex(0) + middle_line->vertex(1)) / 2; + = middle_line->get_boundary().get_new_point_on_line(middle_line); triangulation.vertices_used[next_unused_vertex] = true; // now search a slot for the two @@ -6080,7 +6093,8 @@ namespace internal // set the middle vertex // appropriately - if (quad->at_boundary()) + if (quad->at_boundary() || + (quad->manifold_id() != numbers::invalid_manifold_id) ) triangulation.vertices[next_unused_vertex] = quad->get_boundary().get_new_point_on_quad (quad); else diff --git a/tests/manifold/flat_manifold_01.cc b/tests/manifold/flat_manifold_01.cc deleted file mode 100644 index 9d463df4e4..0000000000 --- a/tests/manifold/flat_manifold_01.cc +++ /dev/null @@ -1,55 +0,0 @@ -//---------------------------- manifold_id_01.cc --------------------------- -// Copyright (C) 2011, 2013 by the mathLab team. -// -// This file is subject to LGPL 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. -// -//---------------------------- flat_manifold_01.cc --------------------------- - - -// Test periodicity conditions using FlatManifold - -#include "../tests.h" -#include -#include - - -// all include files you need here -#include -#include - -template -void test() { - Point period; - for(unsigned int i=0; i manifold(period); - deallog << "Testing spacedim = " << spacedim << std::endl; - deallog << "Period: " << period << std::endl; - std::vector > p(2); - std::vector w(2, 0.5); - double eps = .1; - p[0] = period*eps*2; - p[1] = period*(1.0-eps); - Point middle = manifold.get_new_point(p, w); - deallog << "P(0): " << p[0] << std::endl; - deallog << "P(1): " << p[1] << std::endl; - deallog << "Middle: " << middle << std::endl; -} - -int main () -{ - std::ofstream logfile("output"); - deallog.attach(logfile); - deallog.depth_console(0); - deallog.threshold_double(1.e-5); - - test<1>(); - test<2>(); - test<3>(); - - return 0; -} - diff --git a/tests/manifold/flat_manifold_01.output b/tests/manifold/flat_manifold_01.output deleted file mode 100644 index f6f885fc5b..0000000000 --- a/tests/manifold/flat_manifold_01.output +++ /dev/null @@ -1,16 +0,0 @@ - -DEAL::Testing spacedim = 1 -DEAL::Period: 1.00000 -DEAL::P(0): 0.200000 -DEAL::P(1): 0.900000 -DEAL::Middle: 0.0500000 -DEAL::Testing spacedim = 2 -DEAL::Period: 1.00000 1.00000 -DEAL::P(0): 0.200000 0.200000 -DEAL::P(1): 0.900000 0.900000 -DEAL::Middle: 0.0500000 0.0500000 -DEAL::Testing spacedim = 3 -DEAL::Period: 1.00000 1.00000 1.00000 -DEAL::P(0): 0.200000 0.200000 0.200000 -DEAL::P(1): 0.900000 0.900000 0.900000 -DEAL::Middle: 0.0500000 0.0500000 0.0500000 diff --git a/tests/manifold/manifold_id_05.cc b/tests/manifold/manifold_id_05.cc index 3d94f71390..e6a5a52b55 100644 --- a/tests/manifold/manifold_id_05.cc +++ b/tests/manifold/manifold_id_05.cc @@ -22,7 +22,7 @@ #include #include #include -#include +#include #include // Helper function @@ -35,13 +35,13 @@ void test(unsigned int ref=1){ for(unsigned int i=0; i boundary(center); + HyperBallBoundary boundary(center, .5*std::sqrt((double)dim)); Triangulation tria; GridGenerator::hyper_cube (tria); typename Triangulation::active_cell_iterator cell; tria.begin_active()->face(0)->set_manifold_id(1); - tria.set_manifold(1,boundary); + tria.set_boundary(1,boundary); tria.refine_global(2); diff --git a/tests/manifold/manifold_id_06.cc b/tests/manifold/manifold_id_06.cc index ccf108b4ff..8fa4e892ca 100644 --- a/tests/manifold/manifold_id_06.cc +++ b/tests/manifold/manifold_id_06.cc @@ -22,8 +22,8 @@ #include #include #include -#include #include +#include // Helper function template @@ -37,7 +37,7 @@ void test(unsigned int ref=1){ double radius=center.norm(); - PolarManifold boundary(center); + HyperBallBoundary boundary(center, .25*std::sqrt((double)spacedim)); Triangulation tria; GridGenerator::hyper_cube (tria); typename Triangulation::active_cell_iterator cell; @@ -52,7 +52,7 @@ void test(unsigned int ref=1){ if(cell->face(f)->center().distance(center)< radius) cell->face(f)->set_all_manifold_ids(1); - tria.set_manifold(1,boundary); + tria.set_boundary(1,boundary); tria.refine_global(2); GridOut gridout; diff --git a/tests/manifold/manifold_id_07.cc b/tests/manifold/manifold_id_07.cc index 9430331025..9e0662f45a 100644 --- a/tests/manifold/manifold_id_07.cc +++ b/tests/manifold/manifold_id_07.cc @@ -24,7 +24,7 @@ #include #include #include -#include +#include #include // Helper function @@ -38,12 +38,12 @@ void test(unsigned int ref=1){ center[i] = .5; Triangulation tria; - PolarManifold boundary(center); + HyperBallBoundary boundary(center,.5*std::sqrt((double)dim)); GridGenerator::hyper_cube (tria); typename Triangulation::active_cell_iterator cell; tria.begin_active()->face(0)->set_manifold_id(1); - tria.set_manifold(1,boundary); + tria.set_boundary(1,boundary); tria.refine_global(2);