From: wolf Date: Tue, 8 Sep 1998 14:46:29 +0000 (+0000) Subject: Change file to at least let it compile properly. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f834ea259afda23b4ea20d59c6d68ce85199f176;p=dealii-svn.git Change file to at least let it compile properly. git-svn-id: https://svn.dealii.org/trunk@589 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/Attic/examples/multigrid/multigrid.cc b/deal.II/deal.II/Attic/examples/multigrid/multigrid.cc index b163635c41..e5e252f876 100644 --- a/deal.II/deal.II/Attic/examples/multigrid/multigrid.cc +++ b/deal.II/deal.II/Attic/examples/multigrid/multigrid.cc @@ -13,9 +13,16 @@ #include #include #include +#include #include #include #include +#include + +#include +#include +#include +#include #include #include @@ -56,10 +63,51 @@ class PoissonEquation : public Equation { template class PoissonProblem { public: + /** + * Declare a data type which denotes a + * mapping between a boundary indicator + * and the function denoting the boundary + * values on this part of the boundary. + * Only one boundary function may be given + * for each boundary indicator, which is + * guaranteed by the #map# data type. + * + * See the general documentation of this + * class for more detail. + */ + typedef map*> FunctionMap; + /** + * Typdedef an iterator which assembles + * matrices and vectors. + */ + typedef TriaActiveIterator > active_assemble_iterator; + PoissonProblem (unsigned int order); void clear (); void create_new (); + void solve (); + + /** + * Initiate the process of assemblage of + * vectors and system matrix. Use the + * given equation object and the given + * quadrature formula. Also use the list + * of dirichlet boundary value functions + * (by default, no dirichlet bc are assumed + * which means that all bc are included + * into the weak formulation). + * + * For what exactly happens here, refer to + * the general doc of this class. + */ + virtual void assemble (const Equation &equation, + const Quadrature &q, + const FiniteElement &fe, + const UpdateFlags update_flags, + const FunctionMap &dirichlet_bc = FunctionMap(), + const Boundary &boundary = StraightBoundary()); + int run (unsigned int level); void print_history (string filename) const; @@ -67,6 +115,33 @@ class PoissonProblem { Triangulation *tria; MGDoFHandler *dof; + + /** + * Sparsity pattern of the system matrix. + */ + dSMatrixStruct system_sparsity; + + /** + * System matrix. + */ + dSMatrix system_matrix; + + /** + * Vector storing the right hand side. + */ + dVector right_hand_side; + + /** + * Solution vector. + */ + dVector solution; + + /** + * List of constraints introduced by + * hanging nodes. + */ + ConstraintMatrix constraints; + Function *rhs; Function *boundary_values; @@ -219,7 +294,11 @@ void PoissonProblem::clear () { boundary_values = 0; }; - ProblemBase::clear (); + system_sparsity.reinit (0,0,1); + system_matrix.clear (); + right_hand_side.reinit (1); + solution.reinit (1); + constraints.clear (); }; @@ -231,11 +310,91 @@ void PoissonProblem::create_new () { tria = new Triangulation(); dof = new MGDoFHandler (tria); - set_tria_and_dof (tria, dof); }; +template +void PoissonProblem::assemble (const Equation &equation, + const Quadrature &quadrature, + const FiniteElement &fe, + const UpdateFlags update_flags, + const FunctionMap &dirichlet_bc, + const Boundary &boundary) { + Assert ((tria!=0) && (dof!=0), ExcNoTriaSelected()); + + system_sparsity.reinit (dof->DoFHandler::n_dofs(), + dof->DoFHandler::n_dofs(), + dof->max_couplings_between_dofs()); + right_hand_side.reinit (dof->DoFHandler::n_dofs()); + + // make up sparsity pattern and + // compress with constraints + constraints.clear (); + dof->DoFHandler::make_constraint_matrix (constraints); + dof->DoFHandler::make_sparsity_pattern (system_sparsity); + constraints.condense (system_sparsity); + + // reinite system matrix + system_matrix.reinit (system_sparsity); + // reinit solution vector, preset + // with zeroes. + solution.reinit (dof->DoFHandler::n_dofs()); + + // create assembler + AssemblerData data (*dof, + true, true, //assemble matrix and rhs + system_matrix, + right_hand_side, + quadrature, + fe, + update_flags, + boundary); + active_assemble_iterator assembler (tria, + tria->begin_active()->level(), + tria->begin_active()->index(), + &data); + // loop over all cells, fill matrix and rhs + do + { + assembler->assemble (equation); + } + while ((++assembler).state() == valid); + + // condense system matrix in-place + constraints.condense (system_matrix); + + // condense right hand side in-place + constraints.condense (right_hand_side); + + // apply Dirichlet bc as described + // in the docs + map boundary_value_list; + VectorTools::interpolate_boundary_values (*dof, + dirichlet_bc, fe, boundary, + boundary_value_list); + MatrixTools::apply_boundary_values (boundary_value_list, + system_matrix, solution, + right_hand_side); +}; + + + +template +void PoissonProblem::solve () { + Assert ((tria!=0) && (dof!=0), ExcNoTriaSelected()); + + SolverControl control(4000, 1e-16); + PrimitiveVectorMemory memory(right_hand_side.size()); + SolverCG cg(control,memory); + + // solve + cg.solve (system_matrix, solution, right_hand_side); + // distribute solution + constraints.distribute (solution); +}; + + template int PoissonProblem::run (const unsigned int level) { @@ -304,14 +463,14 @@ int PoissonProblem::run (const unsigned int level) { cout << " Distributing dofs... "; dof->distribute_dofs (*fe); - cout << dof->n_dofs() << " degrees of freedom." << endl; - n_dofs.push_back (dof->n_dofs()); + cout << dof->DoFHandler::n_dofs() << " degrees of freedom." << endl; + n_dofs.push_back (dof->DoFHandler::n_dofs()); cout << " Assembling matrices..." << endl; UpdateFlags update_flags = UpdateFlags(update_q_points | update_gradients | update_JxW_values); - ProblemBase::FunctionMap dirichlet_bc; + map*> dirichlet_bc; dirichlet_bc[0] = boundary_values; assemble (equation, *quadrature, *fe, update_flags, dirichlet_bc); @@ -323,7 +482,7 @@ int PoissonProblem::run (const unsigned int level) { dVector h1_seminorm_error_per_cell, h1_error_per_cell; cout << " Calculating L1 error... "; - VectorTools::integrate_difference (*dof_handler, + VectorTools::integrate_difference (*dof, solution, sol, l1_error_per_cell, *quadrature, *fe, L1_norm); @@ -331,7 +490,7 @@ int PoissonProblem::run (const unsigned int level) { l1_error.push_back (l1_error_per_cell.l1_norm()); cout << " Calculating L2 error... "; - VectorTools::integrate_difference (*dof_handler, + VectorTools::integrate_difference (*dof, solution, sol, l2_error_per_cell, *quadrature, *fe, L2_norm); @@ -339,7 +498,7 @@ int PoissonProblem::run (const unsigned int level) { l2_error.push_back (l2_error_per_cell.l2_norm()); cout << " Calculating L-infinity error... "; - VectorTools::integrate_difference (*dof_handler, + VectorTools::integrate_difference (*dof, solution, sol, linfty_error_per_cell, *quadrature, *fe, Linfty_norm); @@ -347,7 +506,7 @@ int PoissonProblem::run (const unsigned int level) { linfty_error.push_back (linfty_error_per_cell.linfty_norm()); cout << " Calculating H1-seminorm error... "; - VectorTools::integrate_difference (*dof_handler, + VectorTools::integrate_difference (*dof, solution, sol, h1_seminorm_error_per_cell, *quadrature, *fe, H1_seminorm); @@ -355,14 +514,14 @@ int PoissonProblem::run (const unsigned int level) { h1_seminorm_error.push_back (h1_seminorm_error_per_cell.l2_norm()); cout << " Calculating H1 error... "; - VectorTools::integrate_difference (*dof_handler, + VectorTools::integrate_difference (*dof, solution, sol, h1_error_per_cell, *quadrature, *fe, H1_norm); cout << h1_error_per_cell.l2_norm() << endl; h1_error.push_back (h1_error_per_cell.l2_norm()); - if (dof->n_dofs()<=5000) + if (dof->DoFHandler::n_dofs()<=5000) { dVector l1_error_per_dof, l2_error_per_dof, linfty_error_per_dof; dVector h1_seminorm_error_per_dof, h1_error_per_dof; @@ -382,7 +541,7 @@ int PoissonProblem::run (const unsigned int level) { // sol, projected_solution, false, // *boundary_quadrature); // cout << " Calculating L2 error of projected solution... "; -// VectorTools::integrate_difference (*dof_handler, +// VectorTools::integrate_difference (*dof, // projected_solution, sol, // l2_error_per_cell, // *quadrature, *fe, L2_norm); @@ -398,7 +557,9 @@ int PoissonProblem::run (const unsigned int level) { DataOut out; ofstream o(filename.c_str()); - fill_data (out); + out.attach_dof_handler (*dof); + + out.add_data_vector (solution, "u"); out.add_data_vector (l1_error_per_dof, "L1-Error"); out.add_data_vector (l2_error_per_dof, "L2-Error"); out.add_data_vector (linfty_error_per_dof, "Linfty-Error"); @@ -416,7 +577,7 @@ int PoissonProblem::run (const unsigned int level) { delete quadrature; delete boundary_quadrature; - return dof->n_dofs(); + return dof->DoFHandler::n_dofs(); }; diff --git a/tests/big-tests/multigrid/multigrid.cc b/tests/big-tests/multigrid/multigrid.cc index b163635c41..e5e252f876 100644 --- a/tests/big-tests/multigrid/multigrid.cc +++ b/tests/big-tests/multigrid/multigrid.cc @@ -13,9 +13,16 @@ #include #include #include +#include #include #include #include +#include + +#include +#include +#include +#include #include #include @@ -56,10 +63,51 @@ class PoissonEquation : public Equation { template class PoissonProblem { public: + /** + * Declare a data type which denotes a + * mapping between a boundary indicator + * and the function denoting the boundary + * values on this part of the boundary. + * Only one boundary function may be given + * for each boundary indicator, which is + * guaranteed by the #map# data type. + * + * See the general documentation of this + * class for more detail. + */ + typedef map*> FunctionMap; + /** + * Typdedef an iterator which assembles + * matrices and vectors. + */ + typedef TriaActiveIterator > active_assemble_iterator; + PoissonProblem (unsigned int order); void clear (); void create_new (); + void solve (); + + /** + * Initiate the process of assemblage of + * vectors and system matrix. Use the + * given equation object and the given + * quadrature formula. Also use the list + * of dirichlet boundary value functions + * (by default, no dirichlet bc are assumed + * which means that all bc are included + * into the weak formulation). + * + * For what exactly happens here, refer to + * the general doc of this class. + */ + virtual void assemble (const Equation &equation, + const Quadrature &q, + const FiniteElement &fe, + const UpdateFlags update_flags, + const FunctionMap &dirichlet_bc = FunctionMap(), + const Boundary &boundary = StraightBoundary()); + int run (unsigned int level); void print_history (string filename) const; @@ -67,6 +115,33 @@ class PoissonProblem { Triangulation *tria; MGDoFHandler *dof; + + /** + * Sparsity pattern of the system matrix. + */ + dSMatrixStruct system_sparsity; + + /** + * System matrix. + */ + dSMatrix system_matrix; + + /** + * Vector storing the right hand side. + */ + dVector right_hand_side; + + /** + * Solution vector. + */ + dVector solution; + + /** + * List of constraints introduced by + * hanging nodes. + */ + ConstraintMatrix constraints; + Function *rhs; Function *boundary_values; @@ -219,7 +294,11 @@ void PoissonProblem::clear () { boundary_values = 0; }; - ProblemBase::clear (); + system_sparsity.reinit (0,0,1); + system_matrix.clear (); + right_hand_side.reinit (1); + solution.reinit (1); + constraints.clear (); }; @@ -231,11 +310,91 @@ void PoissonProblem::create_new () { tria = new Triangulation(); dof = new MGDoFHandler (tria); - set_tria_and_dof (tria, dof); }; +template +void PoissonProblem::assemble (const Equation &equation, + const Quadrature &quadrature, + const FiniteElement &fe, + const UpdateFlags update_flags, + const FunctionMap &dirichlet_bc, + const Boundary &boundary) { + Assert ((tria!=0) && (dof!=0), ExcNoTriaSelected()); + + system_sparsity.reinit (dof->DoFHandler::n_dofs(), + dof->DoFHandler::n_dofs(), + dof->max_couplings_between_dofs()); + right_hand_side.reinit (dof->DoFHandler::n_dofs()); + + // make up sparsity pattern and + // compress with constraints + constraints.clear (); + dof->DoFHandler::make_constraint_matrix (constraints); + dof->DoFHandler::make_sparsity_pattern (system_sparsity); + constraints.condense (system_sparsity); + + // reinite system matrix + system_matrix.reinit (system_sparsity); + // reinit solution vector, preset + // with zeroes. + solution.reinit (dof->DoFHandler::n_dofs()); + + // create assembler + AssemblerData data (*dof, + true, true, //assemble matrix and rhs + system_matrix, + right_hand_side, + quadrature, + fe, + update_flags, + boundary); + active_assemble_iterator assembler (tria, + tria->begin_active()->level(), + tria->begin_active()->index(), + &data); + // loop over all cells, fill matrix and rhs + do + { + assembler->assemble (equation); + } + while ((++assembler).state() == valid); + + // condense system matrix in-place + constraints.condense (system_matrix); + + // condense right hand side in-place + constraints.condense (right_hand_side); + + // apply Dirichlet bc as described + // in the docs + map boundary_value_list; + VectorTools::interpolate_boundary_values (*dof, + dirichlet_bc, fe, boundary, + boundary_value_list); + MatrixTools::apply_boundary_values (boundary_value_list, + system_matrix, solution, + right_hand_side); +}; + + + +template +void PoissonProblem::solve () { + Assert ((tria!=0) && (dof!=0), ExcNoTriaSelected()); + + SolverControl control(4000, 1e-16); + PrimitiveVectorMemory memory(right_hand_side.size()); + SolverCG cg(control,memory); + + // solve + cg.solve (system_matrix, solution, right_hand_side); + // distribute solution + constraints.distribute (solution); +}; + + template int PoissonProblem::run (const unsigned int level) { @@ -304,14 +463,14 @@ int PoissonProblem::run (const unsigned int level) { cout << " Distributing dofs... "; dof->distribute_dofs (*fe); - cout << dof->n_dofs() << " degrees of freedom." << endl; - n_dofs.push_back (dof->n_dofs()); + cout << dof->DoFHandler::n_dofs() << " degrees of freedom." << endl; + n_dofs.push_back (dof->DoFHandler::n_dofs()); cout << " Assembling matrices..." << endl; UpdateFlags update_flags = UpdateFlags(update_q_points | update_gradients | update_JxW_values); - ProblemBase::FunctionMap dirichlet_bc; + map*> dirichlet_bc; dirichlet_bc[0] = boundary_values; assemble (equation, *quadrature, *fe, update_flags, dirichlet_bc); @@ -323,7 +482,7 @@ int PoissonProblem::run (const unsigned int level) { dVector h1_seminorm_error_per_cell, h1_error_per_cell; cout << " Calculating L1 error... "; - VectorTools::integrate_difference (*dof_handler, + VectorTools::integrate_difference (*dof, solution, sol, l1_error_per_cell, *quadrature, *fe, L1_norm); @@ -331,7 +490,7 @@ int PoissonProblem::run (const unsigned int level) { l1_error.push_back (l1_error_per_cell.l1_norm()); cout << " Calculating L2 error... "; - VectorTools::integrate_difference (*dof_handler, + VectorTools::integrate_difference (*dof, solution, sol, l2_error_per_cell, *quadrature, *fe, L2_norm); @@ -339,7 +498,7 @@ int PoissonProblem::run (const unsigned int level) { l2_error.push_back (l2_error_per_cell.l2_norm()); cout << " Calculating L-infinity error... "; - VectorTools::integrate_difference (*dof_handler, + VectorTools::integrate_difference (*dof, solution, sol, linfty_error_per_cell, *quadrature, *fe, Linfty_norm); @@ -347,7 +506,7 @@ int PoissonProblem::run (const unsigned int level) { linfty_error.push_back (linfty_error_per_cell.linfty_norm()); cout << " Calculating H1-seminorm error... "; - VectorTools::integrate_difference (*dof_handler, + VectorTools::integrate_difference (*dof, solution, sol, h1_seminorm_error_per_cell, *quadrature, *fe, H1_seminorm); @@ -355,14 +514,14 @@ int PoissonProblem::run (const unsigned int level) { h1_seminorm_error.push_back (h1_seminorm_error_per_cell.l2_norm()); cout << " Calculating H1 error... "; - VectorTools::integrate_difference (*dof_handler, + VectorTools::integrate_difference (*dof, solution, sol, h1_error_per_cell, *quadrature, *fe, H1_norm); cout << h1_error_per_cell.l2_norm() << endl; h1_error.push_back (h1_error_per_cell.l2_norm()); - if (dof->n_dofs()<=5000) + if (dof->DoFHandler::n_dofs()<=5000) { dVector l1_error_per_dof, l2_error_per_dof, linfty_error_per_dof; dVector h1_seminorm_error_per_dof, h1_error_per_dof; @@ -382,7 +541,7 @@ int PoissonProblem::run (const unsigned int level) { // sol, projected_solution, false, // *boundary_quadrature); // cout << " Calculating L2 error of projected solution... "; -// VectorTools::integrate_difference (*dof_handler, +// VectorTools::integrate_difference (*dof, // projected_solution, sol, // l2_error_per_cell, // *quadrature, *fe, L2_norm); @@ -398,7 +557,9 @@ int PoissonProblem::run (const unsigned int level) { DataOut out; ofstream o(filename.c_str()); - fill_data (out); + out.attach_dof_handler (*dof); + + out.add_data_vector (solution, "u"); out.add_data_vector (l1_error_per_dof, "L1-Error"); out.add_data_vector (l2_error_per_dof, "L2-Error"); out.add_data_vector (linfty_error_per_dof, "Linfty-Error"); @@ -416,7 +577,7 @@ int PoissonProblem::run (const unsigned int level) { delete quadrature; delete boundary_quadrature; - return dof->n_dofs(); + return dof->DoFHandler::n_dofs(); };