From 4d8555b5509f0eb6304fbe693ef0eb8b17a355fa Mon Sep 17 00:00:00 2001 From: wolf Date: Sat, 20 Feb 1999 12:08:34 +0000 Subject: [PATCH] Add another test taken from /deal.II/examples which checks distribution of degrees of freedom and handling of hanging nodes. git-svn-id: https://svn.dealii.org/trunk@855 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/deal.II/dof_test.cc | 350 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 350 insertions(+) create mode 100644 tests/deal.II/dof_test.cc diff --git a/tests/deal.II/dof_test.cc b/tests/deal.II/dof_test.cc new file mode 100644 index 0000000000..abee91b85f --- /dev/null +++ b/tests/deal.II/dof_test.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include + + +#include +#include + + + +// 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 +class Ball : + public StraightBoundary { + public: + virtual Point + get_new_point_on_line (const typename Triangulation::line_iterator &line) const { + Point middle = StraightBoundary::get_new_point_on_line(line); + + for (int i=0; i + get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const { + Point middle = StraightBoundary::get_new_point_on_quad(quad); + + for (int i=0; i +class CurvedLine : + public StraightBoundary { + public: + virtual Point + get_new_point_on_line (const typename Triangulation::line_iterator &line) const; + + virtual Point + get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; +}; + + + +template +Point +CurvedLine::get_new_point_on_line (const typename Triangulation::line_iterator &line) const +{ + Point middle = StraightBoundary::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 +Point +CurvedLine::get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const +{ + Point middle = StraightBoundary::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 +class TestCases { + public: + TestCases (); + virtual ~TestCases (); + + virtual void create_new (); + virtual void run (const unsigned int testcase); + + private: + Triangulation *tria; + DoFHandler *dof; +}; + + + +template +TestCases::TestCases () : + tria(0), dof(0) {}; + + +template +TestCases::~TestCases () +{ + if (dof) delete dof; + if (tria) delete tria; +}; + + + +template +void TestCases::create_new () { + if (dof != 0) delete dof; + if (tria != 0) delete tria; + + tria = new Triangulation(); + tria->create_hypercube(); + + dof = new DoFHandler (tria); +}; + + + + + +template +void TestCases::run (const unsigned int test_case) { + cout << "Dimension = " << dim + << ", Test case = " << test_case << endl + << endl; + + cout << " Making grid..." << endl; + Boundary *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::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*>(new Ball()) : + static_cast*>(new CurvedLine())); + tria->set_boundary (boundary); + + // refine once + tria->begin_active()->set_refine_flag(); + tria->execute_refinement (); + + Triangulation::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 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; +}; + -- 2.39.5