From d28fde1a4767a311d301576b3ef27c5ddb87777f Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Tue, 10 Nov 2009 04:11:40 +0000 Subject: [PATCH] assemble global system without adaptivity git-svn-id: https://svn.dealii.org/trunk@20081 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-39/step-39.cc | 332 ++++++++++++++++++++++++++++ 1 file changed, 332 insertions(+) diff --git a/deal.II/examples/step-39/step-39.cc b/deal.II/examples/step-39/step-39.cc index c5fa40b45b..638037436f 100644 --- a/deal.II/examples/step-39/step-39.cc +++ b/deal.II/examples/step-39/step-39.cc @@ -9,3 +9,335 @@ /* without copyright and license information. Please refer */ /* to the file deal.II/doc/license.html for the text and */ /* further information on this license. */ + + // The include files for the linear + // algebra: A regular SparseMatrix, + // which in turn will include the + // necessary files for + // SparsityPattern and Vector classes. +#include +#include +#include +#include + + // Include files for setting up the + // mesh +#include + + // Include files for FinitElement + // classes and DoFHandler. +#include +#include +#include + + // The include files for using the + // MeshWorker framework +#include +#include + + // Finally, we take our exact + // solution from the library. +#include + +#include +#include + + // All classes of the deal.II library + // are in the namespace dealii. In + // order to save typing, we tell the + // compiler to search names in there + // as well. +using namespace dealii; + + + // @sect3{The local integrators} + + // MeshWorker separates local + // integration from the loops over + // cells and faces. Thus, we have to + // write local integration classes + // for generating matrices, the right + // hand side and the error + // estimator. + + // All these classes have the same + // three functions for integrating + // over cells, boundary faces and + // interior faces, respectively. All + // the information needed for the + // local integration is provided by + // MeshWorker::IntegrationWorker::CellInfo + // and + // MeshWorker::IntegrationWorker::FaceInfo. Note + // that this public interface cannot + // be changed, because it is expected + // by MeshWorker::integration_loop(). +template +class MatrixIntegrator : public Subscriptor +{ + public: + void cell(typename MeshWorker::IntegrationWorker::CellInfo& info) const; + void bdry(typename MeshWorker::IntegrationWorker::FaceInfo& info) const; + void face(typename MeshWorker::IntegrationWorker::FaceInfo& info1, + typename MeshWorker::IntegrationWorker::FaceInfo& info2) const; +}; + + +template +void MatrixIntegrator::cell(typename MeshWorker::IntegrationWorker::CellInfo& info) const +{ + const FEValuesBase& fe = info.fe(); + FullMatrix& local_matrix = info.M1[0].matrix; + + for (unsigned int k=0; k +void MatrixIntegrator::bdry(typename MeshWorker::IntegrationWorker::FaceInfo& info) const +{ + const FEFaceValuesBase& fe = info.fe(); + FullMatrix& local_matrix = info.M1[0].matrix; + + const unsigned int deg = fe.get_fe().tensor_degree(); + const double penalty = 2. * deg * (deg+1) * info.face->measure() / info.cell->measure(); + + for (unsigned k=0;k +void MatrixIntegrator::face(typename MeshWorker::IntegrationWorker::FaceInfo& info1, + typename MeshWorker::IntegrationWorker::FaceInfo& info2) const +{ + const FEFaceValuesBase& fe = info1.fe(); + FullMatrix& matrix_v1u1 = info1.M1[0].matrix; + FullMatrix& matrix_v1u2 = info1.M2[0].matrix; + FullMatrix& matrix_v2u1 = info2.M2[0].matrix; + FullMatrix& matrix_v2u2 = info2.M1[0].matrix; + + const unsigned int deg = fe.get_fe().tensor_degree(); + const double penalty1 = deg * (deg+1) * info1.face->measure() / info1.cell->measure(); + const double penalty2 = deg * (deg+1) * info2.face->measure() / info2.cell->measure(); + const double penalty = penalty1 + penalty2; + + for (unsigned k=0;k +class RHSIntegrator : public Subscriptor +{ + public: + void cell(typename MeshWorker::IntegrationWorker::CellInfo& info) const; + void bdry(typename MeshWorker::IntegrationWorker::FaceInfo& info) const; + void face(typename MeshWorker::IntegrationWorker::FaceInfo& info1, + typename MeshWorker::IntegrationWorker::FaceInfo& info2) const; +}; + + +template +void RHSIntegrator::cell(typename MeshWorker::IntegrationWorker::CellInfo&) const +{} + + +template +void RHSIntegrator::bdry(typename MeshWorker::IntegrationWorker::FaceInfo& info) const +{ + const FEFaceValuesBase& fe = info.fe(); + Vector& local_vector = info.R[0].block(0); + + static Functions::LSingularityFunction lshaped; + std::vector boundary_values(fe.n_quadrature_points); + lshaped.value_list(fe.get_quadrature_points(), boundary_values); + + const unsigned int deg = fe.get_fe().tensor_degree(); + const double penalty = 2. * deg * (deg+1) * info.face->measure() / info.cell->measure(); + + for (unsigned k=0;k +void RHSIntegrator::face(typename MeshWorker::IntegrationWorker::FaceInfo&, + typename MeshWorker::IntegrationWorker::FaceInfo&) const +{} + + + // @sect3{The main class} +template +class Step39 +{ + public: + Step39(const FiniteElement& fe); + + void run(unsigned int n_steps); + + private: + void setup_system (); + void assemble_matrix (); + void assemble_right_hand_side (); + double estimate (); + void solve (); + void refine_grid (); + void output_results (const unsigned int cycle) const; + + Triangulation triangulation; + const MappingQ1 mapping; + const FiniteElement& fe; + MGDoFHandler dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix matrix; + Vector solution; + Vector right_hand_side; +}; + + +template +Step39::Step39(const FiniteElement& fe) + : + fe(fe), + dof_handler(triangulation) +{ + GridGenerator::hyper_L(triangulation); +} + + +template +void +Step39::setup_system() +{ + dof_handler.distribute_dofs(fe); + unsigned int n_dofs = dof_handler.n_dofs(); + + CompressedSparsityPattern c_sparsity(n_dofs); + const DoFHandler& dof = dof_handler; + DoFTools::make_flux_sparsity_pattern(dof, c_sparsity); + sparsity_pattern.copy_from(c_sparsity); + + matrix.reinit(sparsity_pattern); + + solution.reinit(n_dofs); + right_hand_side.reinit(n_dofs); +} + + +template +void +Step39::assemble_matrix() +{ + const MatrixIntegrator local; + MeshWorker::AssemblingIntegrator >, MatrixIntegrator > + integrator(local); + const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1; + integrator.initialize_gauss_quadrature(n_gauss_points, n_gauss_points, n_gauss_points); + UpdateFlags update_flags = update_values | update_gradients; + integrator.add_update_flags(update_flags, true, true, true, true); + + integrator.initialize(matrix); + MeshWorker::IntegrationInfoBox info_box(dof_handler); + info_box.initialize(integrator, fe, mapping); + MeshWorker::integration_loop(dof_handler.begin_active(), dof_handler.end(), info_box, integrator); +} + + +template +void +Step39::assemble_right_hand_side() +{ + const RHSIntegrator local; + MeshWorker::AssemblingIntegrator >, RHSIntegrator > + integrator(local); + const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1; + integrator.initialize_gauss_quadrature(1, n_gauss_points, n_gauss_points); + UpdateFlags update_flags = update_quadrature_points | update_values | update_gradients; + integrator.add_update_flags(update_flags, true, true, true, true); + + NamedData* > data; + Vector* rhs = &right_hand_side; + data.add(rhs, "RHS"); + integrator.initialize(data); + MeshWorker::IntegrationInfoBox info_box(dof_handler); + info_box.initialize(integrator, fe, mapping); + MeshWorker::integration_loop(dof_handler.begin_active(), dof_handler.end(), info_box, integrator); +} + + +template +void +Step39::solve() +{ + SolverControl control(1000, 1.e-12); + SolverCG > cg(control); + cg.solve(matrix, solution, right_hand_side, PreconditionIdentity()); +} + + +template +void +Step39::run(unsigned int n_steps) +{ + for (unsigned int s=0;s dgq1(1); + Step39<2> test1(dgq1); + test1.run(7); +} -- 2.39.5