]> https://gitweb.dealii.org/ - dealii.git/commitdiff
New test, currently failing.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 26 Jun 2012 14:03:47 +0000 (14:03 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 26 Jun 2012 14:03:47 +0000 (14:03 +0000)
git-svn-id: https://svn.dealii.org/trunk@25646 0785d39b-7218-0410-832d-ea1e28bc413d

tests/deal.II/no_flux_10.cc [new file with mode: 0644]
tests/deal.II/no_flux_10/cmp/generic [new file with mode: 0644]

diff --git a/tests/deal.II/no_flux_10.cc b/tests/deal.II/no_flux_10.cc
new file mode 100644 (file)
index 0000000..d06a733
--- /dev/null
@@ -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 <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+#include <dofs/dof_handler.h>
+#include <fe/fe_system.h>
+#include <fe/mapping_q1.h>
+#include <numerics/vectors.h>
+#include <numerics/data_out.h>
+#include <base/exceptions.h>
+#include <base/function.h>
+
+
+template <int dim>
+static void
+colorize_sixty_deg_hyper_shell(Triangulation<dim> & tria,
+                              const Point<dim>& center,
+                              const double inner_radius,
+                              const double outer_radius);
+
+template<int dim>
+static void sixty_deg_hyper_shell (Triangulation<dim>   &tria,
+                                  const Point<dim>     &center,
+                                  const double        inner_radius,
+                                  const double        outer_radius,
+                                  const unsigned int  n_cells = 0,
+                                  const bool colorize = false);
+
+template <int dim>
+class SixtyDegHyperShellBoundary : public HyperShellBoundary<dim>
+{
+  public:
+
+    SixtyDegHyperShellBoundary (const Point<dim> &center,
+                               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<GeometryInfo<3>::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<GeometryInfo<3>::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<GeometryInfo<3>::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<GeometryInfo<3>::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<GeometryInfo<3>::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<GeometryInfo<3>::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<GeometryInfo<3>::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<Point<3> > 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<CellData<3> > 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<int dim>
+SixtyDegHyperShellBoundary<dim>::SixtyDegHyperShellBoundary (const Point<dim> &center,
+                                                            const double inner_radius,
+                                                            const double outer_radius)
+                :
+                HyperShellBoundary<dim> (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 <int dim>
+void run()
+{
+  Triangulation<dim> triangulation;
+  FESystem<dim> fe(FE_Q<dim>(1), dim);
+  DoFHandler<dim> dof_handler (triangulation);
+  ConstraintMatrix constraints;
+
+  sixty_deg_hyper_shell (triangulation,
+                        Point<dim>(),
+                        0.5,
+                        1.0,
+                        2,
+                        true);
+
+  static SixtyDegHyperShellBoundary<dim> boundary(Point<dim>(),
+                                                 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<unsigned char> 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 (file)
index 0000000..1ac20ab
--- /dev/null
@@ -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

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.