]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Final support for hanging nodes.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 17 Mar 1998 16:00:24 +0000 (16:00 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 17 Mar 1998 16:00:24 +0000 (16:00 +0000)
git-svn-id: https://svn.dealii.org/trunk@73 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/Attic/examples/poisson/poisson.cc
deal.II/deal.II/include/numerics/base.h
deal.II/deal.II/source/dofs/dof_constraints.cc
deal.II/deal.II/source/numerics/base.cc
tests/big-tests/poisson/poisson.cc

index 7b6b860d3191643d3c1628fe898dd4ad3cb6f683..9c4c8b9151884ac79b2d46b7e0b1ca3671cf72a4 100644 (file)
@@ -52,9 +52,13 @@ double PoissonEquation<dim>::right_hand_side (const Point<dim> &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<fe_values.total_dofs; ++j)
          cell_matrix(i,j) += (fe_values.shape_grad(i,point) *
-                              fe_values.shape_grad(j,point) /*+
+                              fe_values.shape_grad(j,point) +
                               fe_values.shape_value(i,point) *
-                              fe_values.shape_value(j,point)*/) *
+                              fe_values.shape_value(j,point)) *
                              fe_values.JxW(point);
        rhs[0](i) += fe_values.shape_value(i,point) *
                     right_hand_side(fe_values.quadrature_point(point)) *
@@ -120,20 +124,54 @@ int main () {
 //  tria.create_hyper_ball(Point<2>(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<dim>::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;
 };
index debc6fd6a5aaf627a60350866817deb5c9877b3f..0a98cc5b5c3a36619a39027d486e3bf0991d98c2 100644 (file)
@@ -7,6 +7,8 @@
 
 #include <lac/dsmatrix.h>
 #include <base/exceptions.h>
+#include <grid/dof_constraints.h>
+
 
 // forward declaration
 template <int dim> class Triangulation;
@@ -121,6 +123,12 @@ class ProblemBase {
                                      */
     dVector             solution;
 
+                                    /**
+                                     * List of constraints introduced by
+                                     * hanging nodes.
+                                     */
+    ConstraintMatrix    constraints;
+    
   friend class Assembler<dim>;
 };
 
index e250b7a28f7973c86d3085c7330bc1700f718e0b..ba588a6a1434a56163471f099c5555a2ef58a68b 100644 (file)
@@ -482,6 +482,8 @@ void ConstraintMatrix::condense (const dVector &uncondensed,
 
 
 void ConstraintMatrix::condense (dVector &vec) const {
+  Assert (sorted == true, ExcMatrixNotClosed());
+  
   vector<ConstraintLine>::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<int> old_line;
+
+  old_line.reserve (uncondensed.n());
+
+  vector<ConstraintLine>::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<ConstraintLine>::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();
 };
index 755a0478503e297f02e0fc3e338d47fb9b621f8e..722129f72073fc3bf73f9f72f96084e2f20ce03f 100644 (file)
@@ -45,25 +45,24 @@ void ProblemBase<dim>::assemble (const Equation<dim>   &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; i<right_hand_sides.size(); ++i)
     right_hand_sides[i]->reinit (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<dim> 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<dim>::assemble (const Equation<dim>   &equation,
 
                                   // condense system matrix in-place
   constraints.condense (system_matrix);
+                                  // condense right hand sides in-place
+  for (unsigned int i=0; i<right_hand_sides.size(); ++i)
+    constraints.condense (*right_hand_sides[i]);
 };
 
 
@@ -89,13 +91,16 @@ void ProblemBase<dim>::assemble (const Equation<dim>   &equation,
 template <int dim>
 void ProblemBase<dim>::solve () {
   int    max_iter  = 4000;
-  double tolerance = 1.e-14;
+  double tolerance = 1.e-16;
   
   Control                          control1(max_iter,tolerance);
   PrimitiveVectorMemory<dVector>   memory(right_hand_sides[0]->n());
   CG<dSMatrix,dVector>             cg(control1,memory);
 
+                                  // solve
   cg (system_matrix, solution, *right_hand_sides[0]);
+                                  // distribute solution
+  constraints.distribute (solution);
 };
 
 
index 7b6b860d3191643d3c1628fe898dd4ad3cb6f683..9c4c8b9151884ac79b2d46b7e0b1ca3671cf72a4 100644 (file)
@@ -52,9 +52,13 @@ double PoissonEquation<dim>::right_hand_side (const Point<dim> &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<fe_values.total_dofs; ++j)
          cell_matrix(i,j) += (fe_values.shape_grad(i,point) *
-                              fe_values.shape_grad(j,point) /*+
+                              fe_values.shape_grad(j,point) +
                               fe_values.shape_value(i,point) *
-                              fe_values.shape_value(j,point)*/) *
+                              fe_values.shape_value(j,point)) *
                              fe_values.JxW(point);
        rhs[0](i) += fe_values.shape_value(i,point) *
                     right_hand_side(fe_values.quadrature_point(point)) *
@@ -120,20 +124,54 @@ int main () {
 //  tria.create_hyper_ball(Point<2>(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<dim>::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;
 };

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.