From c7e52a2cbf759590c39d2556d9ba567c5621309c Mon Sep 17 00:00:00 2001 From: bangerth Date: Tue, 26 Jun 2012 14:03:47 +0000 Subject: [PATCH] New test, currently failing. git-svn-id: https://svn.dealii.org/trunk@25646 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/deal.II/no_flux_10.cc | 265 +++++++++++++++++++++++++++ tests/deal.II/no_flux_10/cmp/generic | 22 +++ 2 files changed, 287 insertions(+) create mode 100644 tests/deal.II/no_flux_10.cc create mode 100644 tests/deal.II/no_flux_10/cmp/generic diff --git a/tests/deal.II/no_flux_10.cc b/tests/deal.II/no_flux_10.cc new file mode 100644 index 0000000000..d06a73394b --- /dev/null +++ b/tests/deal.II/no_flux_10.cc @@ -0,0 +1,265 @@ +//---------------------------- vectors_rhs_hp_02.cc --------------------------- +// $Id: vectors_rhs_hp_02.cc 23710 2011-05-17 04:50:10Z bangerth $ +// Version: $Name$ +// +// Copyright (C) 2012 by the deal.II authors +// +// 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. +// +//---------------------------- --------------------------- + + +// test by Jennifer Worthen -- we get an error of the following kind: +// -------------------------------------------------------- +// An error occurred in line <304> of file <.../constraint_matrix.cc> in function +// void dealii::ConstraintMatrix::close() +// The violated condition was: +// dof_index != line->line +// The name and call sequence of the exception was: +// ExcMessage ("Cycle in constraints detected!") +// Additional Information: +// Cycle in constraints detected! +// -------------------------------------------------------- + + + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +template +static void +colorize_sixty_deg_hyper_shell(Triangulation & tria, + const Point& center, + const double inner_radius, + const double outer_radius); + +template +static void sixty_deg_hyper_shell (Triangulation &tria, + const Point ¢er, + const double inner_radius, + const double outer_radius, + const unsigned int n_cells = 0, + const bool colorize = false); + +template +class SixtyDegHyperShellBoundary : public HyperShellBoundary +{ + public: + + SixtyDegHyperShellBoundary (const Point ¢er, + const double inner_radius, + const double outer_radius); + + private: + const double inner_radius; + const double outer_radius; +}; + +// Implementation for 3D only +template <> +void +colorize_sixty_deg_hyper_shell(Triangulation<3> & tria, + const Point<3>& center, + const double inner_radius, + const double outer_radius) +{ + + // if (tria.n_cells() != 4) + // AssertThrow (false, ExcNotImplemented()); + + double middle = (outer_radius-inner_radius)/2e0 + inner_radius; + double eps = 1e-3*middle; + Triangulation<3>::cell_iterator cell = tria.begin(); + + for (;cell!=tria.end();++cell) + for(unsigned int f=0; f::faces_per_cell; ++f) + { + if(!cell->face(f)->at_boundary()) + continue; + + double radius = cell->face(f)->center().norm() - center.norm(); + if (std::fabs(cell->face(f)->center()(2) - sqrt(3.)*cell->face(f)->center()(0)) < eps ) // z = sqrt(3)x set boundary 2 + { + cell->face(f)->set_boundary_indicator(2); + for (unsigned int j=0;j::lines_per_face;++j) + if(cell->face(f)->line(j)->at_boundary()) + if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) > eps) + cell->face(f)->line(j)->set_boundary_indicator(2); + } + else if (std::fabs(cell->face(f)->center()(2) + sqrt(3.)*cell->face(f)->center()(0)) < eps) // z = -sqrt(3)x set boundary 3 + { + cell->face(f)->set_boundary_indicator(3); + for (unsigned int j=0;j::lines_per_face;++j) + if(cell->face(f)->line(j)->at_boundary()) + if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) > eps) + cell->face(f)->line(j)->set_boundary_indicator(3); + } + else if (std::fabs(cell->face(f)->center()(2) - sqrt(3.)*cell->face(f)->center()(1)) < eps ) // z = sqrt(3)y set boundary 4 + { + cell->face(f)->set_boundary_indicator(4); + for (unsigned int j=0;j::lines_per_face;++j) + if(cell->face(f)->line(j)->at_boundary()) + if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) > eps) + cell->face(f)->line(j)->set_boundary_indicator(4); + } + else if (std::fabs(cell->face(f)->center()(2) + sqrt(3.)*cell->face(f)->center()(1)) < eps ) // z = -sqrt(3)y set boundary 5 + { + cell->face(f)->set_boundary_indicator(5); + for (unsigned int j=0;j::lines_per_face;++j) + if(cell->face(f)->line(j)->at_boundary()) + if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) > eps) + cell->face(f)->line(j)->set_boundary_indicator(5); + } + else if (radius < middle) // inner radius set boundary 0 + { + cell->face(f)->set_boundary_indicator(0); + for (unsigned int j=0;j::lines_per_face;++j) + if(cell->face(f)->line(j)->at_boundary()) + if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) < eps) + cell->face(f)->line(j)->set_boundary_indicator(0); + } + else if (radius > middle) // outer radius set boundary 1 + { + cell->face(f)->set_boundary_indicator(1); + for (unsigned int j=0;j::lines_per_face;++j) + if(cell->face(f)->line(j)->at_boundary()) + if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) < eps) + cell->face(f)->line(j)->set_boundary_indicator(1); + } + else + AssertThrow (false, ExcInternalError()); + } +} + +// Implementation for 3D only +template <> +void sixty_deg_hyper_shell (Triangulation<3> & tria, + const Point<3>& center, + const double inner_radius, + const double outer_radius, + const unsigned int n, + const bool colorize) +{ + // Assert ((inner_radius > 0) && (inner_radius < outer_radius), + // ExcInvalidRadii ()); + if (n == 0 || n == 2) + { + const double r0 = inner_radius; + const double r1 = outer_radius; + + std::vector > vertices; + + vertices.push_back (center+Point<3>( 1.0/sqrt(5.)*r0, 1.0/sqrt(5.)*r0, sqrt(3./5.)*r0)); //8 -> 0 + vertices.push_back (center+Point<3>( 1.0/sqrt(5.)*r1, 1.0/sqrt(5.)*r1, sqrt(3./5.)*r1)); //9 -> 1 + vertices.push_back (center+Point<3>( 1.0/sqrt(5.)*r0, -1.0/sqrt(5.)*r0, sqrt(3./5.)*r0)); //10 -> 2 + vertices.push_back (center+Point<3>( 1.0/sqrt(5.)*r1, -1.0/sqrt(5.)*r1, sqrt(3./5.)*r1)); //11 -> 3 + vertices.push_back (center+Point<3>( -1.0/sqrt(5.)*r0, 1.0/sqrt(5.)*r0, sqrt(3./5.)*r0)); //14 -> 4 + vertices.push_back (center+Point<3>( -1.0/sqrt(5.)*r1, 1.0/sqrt(5.)*r1, sqrt(3./5.)*r1)); //15 -> 5 + vertices.push_back (center+Point<3>( -1.0/sqrt(5.)*r0, -1.0/sqrt(5.)*r0, sqrt(3./5.)*r0)); //16 -> 6 + vertices.push_back (center+Point<3>( -1.0/sqrt(5.)*r1, -1.0/sqrt(5.)*r1, sqrt(3./5.)*r1)); //17 -> 7 + + const int cell_vertices[1][8] = { + {6, 2, 4, 0, 7, 3, 5, 1}, + }; + + std::vector > cells(1); + + for (unsigned int i=0; i<1; ++i) + { + for (unsigned int j=0; j<8; ++j) + cells[i].vertices[j] = cell_vertices[i][j]; + cells[i].material_id = 0; + } + + tria.create_triangulation ( vertices, cells, SubCellData()); // no boundary information + } + else + { + AssertThrow(false, ExcNotImplemented()); + } + + if (colorize) + colorize_sixty_deg_hyper_shell(tria, center, inner_radius, outer_radius); +} + +template +SixtyDegHyperShellBoundary::SixtyDegHyperShellBoundary (const Point ¢er, + const double inner_radius, + const double outer_radius) + : + HyperShellBoundary (center), + inner_radius (inner_radius), + outer_radius (outer_radius) +{ + if (dim > 2) + Assert ((inner_radius >= 0) && + (outer_radius > 0) && + (outer_radius > inner_radius), + ExcMessage ("Inner and outer radii must be specified explicitly in 3d.")); +} + +template +void run() +{ + Triangulation triangulation; + FESystem fe(FE_Q(1), dim); + DoFHandler dof_handler (triangulation); + ConstraintMatrix constraints; + + sixty_deg_hyper_shell (triangulation, + Point(), + 0.5, + 1.0, + 2, + true); + + static SixtyDegHyperShellBoundary boundary(Point(), + 0.5, + 1.0); + triangulation.set_boundary (0, boundary); + triangulation.set_boundary (1, boundary); + + triangulation.refine_global(2); + + dof_handler.distribute_dofs (fe); + + std::set no_normal_flux_boundaries; + no_normal_flux_boundaries.insert (0); + no_normal_flux_boundaries.insert (2); + VectorTools::compute_no_normal_flux_constraints + (dof_handler, 0, + no_normal_flux_boundaries, + constraints); + + constraints.close(); + constraints.print(deallog.get_file_stream()); + + + deallog << "OK" << std::endl; +} + + +int main () +{ + std::ofstream logfile ("no_flux_10/output"); + logfile.precision (4); + logfile.setf(std::ios::fixed); + deallog.attach(logfile); + deallog.depth_console (0); + + run<3> (); +} diff --git a/tests/deal.II/no_flux_10/cmp/generic b/tests/deal.II/no_flux_10/cmp/generic new file mode 100644 index 0000000000..1ac20abe9d --- /dev/null +++ b/tests/deal.II/no_flux_10/cmp/generic @@ -0,0 +1,22 @@ + + 2 = 0 + 5 = 0 + 8 = 0 + 7 6: -0.4142 + 11 9: 1.0000 + 11 10: 1.0000 + 19 18: -0.4889 + 19 20: -0.9881 + 22 21: -1.0000 + 22 23: -0.8648 + 25 = 0 + 26 = 0 + 27 = 0 + 28 = 0 + 29 = 0 + 31 = 0 + 34 = 0 + 35 33: -1.0000 + 37 = 0 + 40 = 0 + 41 39: -0.4142 -- 2.39.5