]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add another test taken from /deal.II/examples which checks distribution of degrees...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 20 Feb 1999 12:08:34 +0000 (12:08 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 20 Feb 1999 12:08:34 +0000 (12:08 +0000)
git-svn-id: https://svn.dealii.org/trunk@855 0785d39b-7218-0410-832d-ea1e28bc413d

tests/deal.II/dof_test.cc [new file with mode: 0644]

diff --git a/tests/deal.II/dof_test.cc b/tests/deal.II/dof_test.cc
new file mode 100644 (file)
index 0000000..abee91b
--- /dev/null
@@ -0,0 +1,350 @@
+/* $Id$ */
+/* Copyright W. Bangerth, University of Heidelberg, 1998 */
+
+// deal_II_libraries.g=-ldeal_II_2d.g -ldeal_II_3d.g
+// deal_II_libraries=-ldeal_II_2d -ldeal_II_3d
+
+
+
+#include <grid/dof.h>
+#include <grid/tria.h>
+#include <fe/fe_lib.lagrange.h>
+#include <grid/tria_boundary.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_accessor.h>
+#include <lac/dsmatrix.h>
+#include <grid/dof_constraints.h>
+
+
+#include <cmath>
+#include <cstdlib>
+
+
+
+// 1: continuous refinement of the unit square always in the middle
+// 2: refinement of the circle at the boundary
+// 2: refinement of a wiggled area at the boundary
+
+
+
+
+
+
+template <int dim>
+class Ball :
+  public StraightBoundary<dim> {
+  public:
+    virtual Point<dim>
+    get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const {
+      Point<dim> middle = StraightBoundary<dim>::get_new_point_on_line(line);
+      
+      for (int i=0; i<dim; ++i)
+       middle(i) -= .5;
+      middle *= sqrt(dim) / (sqrt(middle.square())*2);
+      for (int i=0; i<dim; ++i)
+       middle(i) += .5;
+      
+      return middle;
+    };
+
+    
+    virtual Point<dim>
+    get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const {
+      Point<dim> middle = StraightBoundary<dim>::get_new_point_on_quad(quad);
+      
+      for (int i=0; i<dim; ++i)
+       middle(i) -= .5;
+      middle *= sqrt(dim) / (sqrt(middle.square())*2);
+      for (int i=0; i<dim; ++i)
+       middle(i) += .5;
+      
+      return middle;
+    };
+};
+
+
+template <int dim>
+class CurvedLine :
+  public StraightBoundary<dim> {
+  public:
+    virtual Point<dim>
+    get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const;
+
+    virtual Point<dim>
+    get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const;
+};
+
+
+
+template <int dim>
+Point<dim>
+CurvedLine<dim>::get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
+{
+  Point<dim> middle = StraightBoundary<dim>::get_new_point_on_line (line);
+
+                                  // if the line is at the top of bottom
+                                  // face: do a special treatment on
+                                  // this line. Note that if the
+                                  // z-value of the midpoint is either
+                                  // 0 or 1, then the z-values of all
+                                  // vertices of the line is like that
+  if (dim>=3)
+    if (((middle(2) == 0) || (middle(2) == 1))
+                                      // find out, if the line is in the
+                                      // interior of the top or bottom face
+                                      // of the domain, or at the edge.
+                                      // lines at the edge need to undergo
+                                      // the usual treatment, while for
+                                      // interior lines taking the midpoint
+                                      // is sufficient
+                                      //
+                                      // note: the trick with the boundary
+                                      // id was invented after the above was
+                                      // written, so we are not very strict
+                                      // here with using these flags
+       && (line->boundary_indicator() == 1))
+      return middle;
+
+
+  double x=middle(0),
+        y=middle(1);
+  
+  if (y<x)
+    if (y<1-x)
+      middle(1) = 0.04*sin(6*3.141592*middle(0));
+    else
+      middle(0) = 1+0.04*sin(6*3.141592*middle(1));
+  
+  else
+    if (y<1-x)
+      middle(0) = 0.04*sin(6*3.141592*middle(1));
+    else
+      middle(1) = 1+0.04*sin(6*3.141592*middle(0));
+  
+  return middle;
+};
+
+
+
+template <int dim>
+Point<dim>
+CurvedLine<dim>::get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const
+{
+  Point<dim> middle = StraightBoundary<dim>::get_new_point_on_quad (quad);
+
+                                  // if the face is at the top of bottom
+                                  // face: do not move the midpoint in
+                                  // x/y direction. Note that if the
+                                  // z-value of the midpoint is either
+                                  // 0 or 1, then the z-values of all
+                                  // vertices of the quad is like that
+  if ((middle(2) == 0) || (middle(2) == 1))
+    return middle;
+  
+  double x=middle(0),
+        y=middle(1);
+  
+  if (y<x)
+    if (y<1-x)
+      middle(1) = 0.04*sin(6*3.141592*middle(0));
+    else
+      middle(0) = 1+0.04*sin(6*3.141592*middle(1));
+  
+  else
+    if (y<1-x)
+      middle(0) = 0.04*sin(6*3.141592*middle(1));
+    else
+      middle(1) = 1+0.04*sin(6*3.141592*middle(0));
+  
+  return middle;
+};
+
+
+
+
+template <int dim>
+class TestCases {
+  public:
+    TestCases ();
+    virtual ~TestCases ();
+    
+    virtual void create_new ();
+    virtual void run (const unsigned int testcase);
+
+  private:
+    Triangulation<dim> *tria;
+    DoFHandler<dim>    *dof;
+};
+
+
+
+template <int dim>
+TestCases<dim>::TestCases () :
+               tria(0), dof(0) {};
+
+
+template <int dim>
+TestCases<dim>::~TestCases () 
+{
+  if (dof)  delete dof;
+  if (tria) delete tria;
+};
+
+
+
+template <int dim>
+void TestCases<dim>::create_new () {
+  if (dof  != 0) delete dof;
+  if (tria != 0) delete tria;
+
+  tria = new Triangulation<dim>();
+  tria->create_hypercube();
+
+  dof = new DoFHandler<dim> (tria);
+};
+
+
+
+
+
+template <int dim>
+void TestCases<dim>::run (const unsigned int test_case) {
+  cout << "Dimension = " << dim
+       << ", Test case = " << test_case << endl
+       << endl;
+  
+  cout << "    Making grid..." << endl;  
+  Boundary<dim> *boundary = 0;
+  
+  switch (test_case) 
+    {
+      case 1: 
+      {
+                                        // refine first cell
+       tria->begin_active()->set_refine_flag();
+       tria->execute_refinement ();
+                                        // refine first active cell
+                                        // on coarsest level
+       tria->begin_active()->set_refine_flag ();
+       tria->execute_refinement ();
+
+       Triangulation<dim>::active_cell_iterator cell;
+       for (int i=0; i<(dim==2 ? 12 : 7); ++i) 
+         {
+                                            // refine the presently
+                                            // second last cell 17
+                                            // times
+           cell = tria->last_active(tria->n_levels()-1);
+           --cell;
+           cell->set_refine_flag ();
+           tria->execute_refinement ();
+         };
+
+       break;
+      }
+      
+      case 2:
+      case 3:
+      {
+       if (dim==3)
+         {
+           tria->begin_active()->face(2)->set_boundary_indicator(1);
+           tria->begin_active()->face(4)->set_boundary_indicator(1);
+         };
+       
+                                        // set the boundary function
+       boundary = (test_case==2 ?
+                   static_cast<Boundary<dim>*>(new Ball<dim>()) :
+                   static_cast<Boundary<dim>*>(new CurvedLine<dim>()));
+       tria->set_boundary (boundary);
+       
+                                        // refine once
+       tria->begin_active()->set_refine_flag();
+       tria->execute_refinement ();
+       
+       Triangulation<dim>::active_cell_iterator cell, endc;
+       for (int i=0; i<5-dim; ++i) 
+         {
+           cell = tria->begin_active();
+           endc = tria->end();
+           
+                                            // refine all
+                                            // boundary cells
+           for (; cell!=endc; ++cell)
+             if (cell->at_boundary())
+               cell->set_refine_flag();
+           
+           tria->execute_refinement();
+         };
+       
+       break;
+      }
+    };
+
+
+  cout << "    Distributing degrees of freedom..." << endl;
+  FELinear<dim> fe;
+  dof->distribute_dofs (fe);
+
+  cout << "    Renumbering degrees of freedom..." << endl;
+  dof->renumber_dofs (Cuthill_McKee, false);
+    
+  dSMatrixStruct sparsity (dof->n_dofs(),
+                          dof->max_couplings_between_dofs());
+  
+  
+  dof->make_sparsity_pattern (sparsity);
+  int unconstrained_bandwidth = sparsity.bandwidth();
+
+  cout << "    Writing sparsity pattern..." << endl;
+  sparsity.print_gnuplot (cout);
+
+
+  
+                                  // computing constraints
+  cout << "    Computing constraints..." << endl;
+  ConstraintMatrix constraints;
+  dof->make_constraint_matrix (constraints);
+  constraints.condense (sparsity);
+  
+  cout << "    Writing condensed sparsity pattern..." << endl;
+  sparsity.print_gnuplot (cout);
+
+
+  cout << endl
+       << "    Total number of cells         = " << tria->n_cells() << endl
+       << "    Total number of active cells  = " << tria->n_active_cells() << endl
+       << "    Number of DoFs                = " << dof->n_dofs() << endl
+       << "    Number of constraints         = " << constraints.n_constraints() << endl
+       << "    Unconstrained matrix bandwidth= " << unconstrained_bandwidth << endl
+       << "    Constrained matrix bandwidth  = " << sparsity.bandwidth()
+       << endl << endl;
+
+                                  // release the lock that dof has to the
+                                  // finite element object
+  dof->clear ();
+  tria->set_boundary (0);
+  if (boundary)
+    delete boundary;
+};
+
+
+
+int main () {
+  for (unsigned int test_case=1; test_case<=3; ++test_case)
+    {
+      TestCases<2> tests;
+      tests.create_new ();
+      tests.run (test_case);
+    };
+
+  for (unsigned int test_case=1; test_case<=3; ++test_case)
+    {
+      TestCases<3> tests;
+      tests.create_new ();
+      tests.run (test_case);
+    };
+  
+  return 0;
+};
+

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.