From 3baea0d2103530118e0cb04122302a598539df9f Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 17 Mar 1998 16:00:24 +0000 Subject: [PATCH] Final support for hanging nodes. git-svn-id: https://svn.dealii.org/trunk@73 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/Attic/examples/poisson/poisson.cc | 54 +++++++++++--- deal.II/deal.II/include/numerics/base.h | 8 +++ .../deal.II/source/dofs/dof_constraints.cc | 71 +++++++++++++++++++ deal.II/deal.II/source/numerics/base.cc | 15 ++-- tests/big-tests/poisson/poisson.cc | 54 +++++++++++--- 5 files changed, 181 insertions(+), 21 deletions(-) diff --git a/deal.II/deal.II/Attic/examples/poisson/poisson.cc b/deal.II/deal.II/Attic/examples/poisson/poisson.cc index 7b6b860d31..9c4c8b9151 100644 --- a/deal.II/deal.II/Attic/examples/poisson/poisson.cc +++ b/deal.II/deal.II/Attic/examples/poisson/poisson.cc @@ -52,9 +52,13 @@ double PoissonEquation::right_hand_side (const Point &p) const { // cos(2*3.1415926536*p(0))); return p(0)*p(0)*p(0)-3./2.*p(0)*p(0)-6*p(0)+3; case 2: - return ((1-3.1415926536*3.1415926536) * - cos(3.1415926536*p(0)) * - cos(3.1415926536*p(1))); +// return ((1-3.1415926536*3.1415926536) * +// cos(3.1415926536*p(0)) * +// cos(3.1415926536*p(1))); + return (p(0)*p(0)*p(0)+p(1)*p(1)*p(1) + - 3./2.*(p(0)*p(0)+p(1)*p(1)) + - 6*(p(0)+p(1)) + + 6); default: return 0; }; @@ -92,9 +96,9 @@ void PoissonEquation<2>::assemble (dFMatrix &cell_matrix, { for (unsigned int j=0; j(2,3),4); // tria.set_boundary (&boundary); - tria.refine_global (1); +/* tria.refine_global (1); tria.begin_active()->set_refine_flag(); tria.execute_refinement (); + tria.refine_global (3); +*/ + + cout << "Making grid..." << endl; + const unsigned int dim=2; + tria.refine_global (1); + + Triangulation::active_cell_iterator cell, endc; + for (int i=0; i<8; ++i) + { + int n_levels = tria.n_levels(); + cell = tria.begin_active(); + endc = tria.end(); + + for (; cell!=endc; ++cell) + { + double r = rand()*1.0/RAND_MAX, + weight = 1.* + (cell->level()*cell->level()) / + (n_levels*n_levels); + + if (r <= 0.5*weight) + cell->set_refine_flag (); + }; + + tria.execute_refinement (); + }; + tria.refine_global (1); + + cout << "Distributing dofs..." << endl; dof.distribute_dofs (fe); + + cout << "Assembling matrices..." << endl; problem.assemble (equation, quadrature, fe); - -// problem.solve (); + cout << "Solving..." << endl; + problem.solve (); + cout << "Printing..." << endl; DataOut<2> out; ofstream gnuplot("gnuplot.out.5"); problem.fill_data (out); out.write_gnuplot (gnuplot); + gnuplot.close (); return 0; }; diff --git a/deal.II/deal.II/include/numerics/base.h b/deal.II/deal.II/include/numerics/base.h index debc6fd6a5..0a98cc5b5c 100644 --- a/deal.II/deal.II/include/numerics/base.h +++ b/deal.II/deal.II/include/numerics/base.h @@ -7,6 +7,8 @@ #include #include +#include + // forward declaration template class Triangulation; @@ -121,6 +123,12 @@ class ProblemBase { */ dVector solution; + /** + * List of constraints introduced by + * hanging nodes. + */ + ConstraintMatrix constraints; + friend class Assembler; }; diff --git a/deal.II/deal.II/source/dofs/dof_constraints.cc b/deal.II/deal.II/source/dofs/dof_constraints.cc index e250b7a28f..ba588a6a14 100644 --- a/deal.II/deal.II/source/dofs/dof_constraints.cc +++ b/deal.II/deal.II/source/dofs/dof_constraints.cc @@ -482,6 +482,8 @@ void ConstraintMatrix::condense (const dVector &uncondensed, void ConstraintMatrix::condense (dVector &vec) const { + Assert (sorted == true, ExcMatrixNotClosed()); + vector::const_iterator next_constraint = lines.begin(); for (unsigned int row=0; row<(unsigned int)vec.n(); ++row) if (row == (*next_constraint).line) @@ -501,6 +503,75 @@ void ConstraintMatrix::condense (dVector &vec) const { +void ConstraintMatrix::distribute (const dVector &condensed, + dVector &uncondensed) const { + Assert (sorted == true, ExcMatrixNotClosed()); + Assert ((unsigned int)condensed.n()+n_constraints() == (unsigned int)uncondensed.n(), + ExcWrongDimension()); + + // store for each line of the new vector + // its old line number before + // distribution. If the shift is + // -1, this line was condensed away + vector old_line; + + old_line.reserve (uncondensed.n()); + + vector::const_iterator next_constraint = lines.begin(); + unsigned int shift = 0; + unsigned int n_rows = (unsigned int)uncondensed.n(); + for (unsigned int row=0; row!=n_rows; ++row) + if (row == (*next_constraint).line) + { + // this line is constrained + old_line.push_back (-1); + // note that #lines# is ordered + ++next_constraint; + ++shift; + } + else + old_line.push_back (row-shift); + + + next_constraint = lines.begin(); + for (unsigned int line=0; line<(unsigned int)uncondensed.n(); ++line) + if (old_line[line] != -1) + // line was not condensed away + uncondensed(line) = condensed(old_line[line]); + else + { + // line was condensed away, create it newly + // first set it to zero + uncondensed(line) = 0.; + // then add the different contributions + for (unsigned int i=0; i<(*next_constraint).entries.size(); ++i) + uncondensed(line) += (condensed(old_line[(*next_constraint).entries[i].first]) * + (*next_constraint).entries[i].second); + ++next_constraint; + }; +}; + + + +void ConstraintMatrix::distribute (dVector &vec) const { + Assert (sorted == true, ExcMatrixNotClosed()); + + vector::const_iterator next_constraint = lines.begin(); + for (; next_constraint != lines.end(); ++next_constraint) + { + // make entry in line next_constraint.line + // first set it to zero + vec((*next_constraint).line) = 0.; + // then add the different contributions + for (unsigned int i=0; i<(*next_constraint).entries.size(); ++i) + vec((*next_constraint).line) += (vec((*next_constraint).entries[i].first) * + (*next_constraint).entries[i].second); + }; +}; + + + + unsigned int ConstraintMatrix::n_constraints () const { return lines.size(); }; diff --git a/deal.II/deal.II/source/numerics/base.cc b/deal.II/deal.II/source/numerics/base.cc index 755a047850..722129f720 100644 --- a/deal.II/deal.II/source/numerics/base.cc +++ b/deal.II/deal.II/source/numerics/base.cc @@ -45,25 +45,24 @@ void ProblemBase::assemble (const Equation &equation, system_sparsity.reinit (dof_handler->n_dofs(), dof_handler->n_dofs(), dof_handler->max_couplings_between_dofs()); - solution.reinit (dof_handler->n_dofs()); for (unsigned int i=0; ireinit (dof_handler->n_dofs()); // make up sparsity pattern and // compress with constraints - ConstraintMatrix constraints; - + constraints.clear (); dof_handler->make_constraint_matrix (constraints); dof_handler->make_sparsity_pattern (system_sparsity); constraints.condense (system_sparsity); // reinite system matrix system_matrix.reinit (system_sparsity); + solution.reinit (dof_handler->n_dofs()); // create assembler AssemblerData data (*dof_handler, true, true, //assemble matrix and rhs - right_hand_sides.size(), //one rhs + right_hand_sides.size(), *this, quadrature, fe); @@ -81,6 +80,9 @@ void ProblemBase::assemble (const Equation &equation, // condense system matrix in-place constraints.condense (system_matrix); + // condense right hand sides in-place + for (unsigned int i=0; i::assemble (const Equation &equation, template void ProblemBase::solve () { int max_iter = 4000; - double tolerance = 1.e-14; + double tolerance = 1.e-16; Control control1(max_iter,tolerance); PrimitiveVectorMemory memory(right_hand_sides[0]->n()); CG cg(control1,memory); + // solve cg (system_matrix, solution, *right_hand_sides[0]); + // distribute solution + constraints.distribute (solution); }; diff --git a/tests/big-tests/poisson/poisson.cc b/tests/big-tests/poisson/poisson.cc index 7b6b860d31..9c4c8b9151 100644 --- a/tests/big-tests/poisson/poisson.cc +++ b/tests/big-tests/poisson/poisson.cc @@ -52,9 +52,13 @@ double PoissonEquation::right_hand_side (const Point &p) const { // cos(2*3.1415926536*p(0))); return p(0)*p(0)*p(0)-3./2.*p(0)*p(0)-6*p(0)+3; case 2: - return ((1-3.1415926536*3.1415926536) * - cos(3.1415926536*p(0)) * - cos(3.1415926536*p(1))); +// return ((1-3.1415926536*3.1415926536) * +// cos(3.1415926536*p(0)) * +// cos(3.1415926536*p(1))); + return (p(0)*p(0)*p(0)+p(1)*p(1)*p(1) + - 3./2.*(p(0)*p(0)+p(1)*p(1)) + - 6*(p(0)+p(1)) + + 6); default: return 0; }; @@ -92,9 +96,9 @@ void PoissonEquation<2>::assemble (dFMatrix &cell_matrix, { for (unsigned int j=0; j(2,3),4); // tria.set_boundary (&boundary); - tria.refine_global (1); +/* tria.refine_global (1); tria.begin_active()->set_refine_flag(); tria.execute_refinement (); + tria.refine_global (3); +*/ + + cout << "Making grid..." << endl; + const unsigned int dim=2; + tria.refine_global (1); + + Triangulation::active_cell_iterator cell, endc; + for (int i=0; i<8; ++i) + { + int n_levels = tria.n_levels(); + cell = tria.begin_active(); + endc = tria.end(); + + for (; cell!=endc; ++cell) + { + double r = rand()*1.0/RAND_MAX, + weight = 1.* + (cell->level()*cell->level()) / + (n_levels*n_levels); + + if (r <= 0.5*weight) + cell->set_refine_flag (); + }; + + tria.execute_refinement (); + }; + tria.refine_global (1); + + cout << "Distributing dofs..." << endl; dof.distribute_dofs (fe); + + cout << "Assembling matrices..." << endl; problem.assemble (equation, quadrature, fe); - -// problem.solve (); + cout << "Solving..." << endl; + problem.solve (); + cout << "Printing..." << endl; DataOut<2> out; ofstream gnuplot("gnuplot.out.5"); problem.fill_data (out); out.write_gnuplot (gnuplot); + gnuplot.close (); return 0; }; -- 2.39.5