From 276c6cadca2d2da28db8bb9a31c2733d6b63ae51 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Mon, 25 Jan 2010 16:38:19 +0000 Subject: [PATCH] step 39 compiles, but crashes in estimate, because hessians are not computed git-svn-id: https://svn.dealii.org/trunk@20447 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/numerics/mesh_worker_info.h | 8 +- deal.II/examples/step-39/step-39.cc | 152 +++++++++++------- 2 files changed, 100 insertions(+), 60 deletions(-) diff --git a/deal.II/deal.II/include/numerics/mesh_worker_info.h b/deal.II/deal.II/include/numerics/mesh_worker_info.h index 4a432bf09f..36fedf501c 100644 --- a/deal.II/deal.II/include/numerics/mesh_worker_info.h +++ b/deal.II/deal.II/include/numerics/mesh_worker_info.h @@ -560,8 +560,9 @@ namespace MeshWorker const FiniteElement& el, const Mapping& mapping); - template + template void initialize(const WORKER&, + ASSEMBLER &assembler, const FiniteElement& el, const Mapping& mapping, const NamedData& data); @@ -877,15 +878,16 @@ namespace MeshWorker template - template + template void IntegrationInfoBox::initialize( const WORKER& integrator, + ASSEMBLER &assembler, const FiniteElement& el, const Mapping& mapping, const NamedData& data) { - initialize(integrator, el, mapping); + initialize(integrator, assembler, el, mapping); boost::shared_ptr > p; p = boost::shared_ptr >(new VectorData (integrator.cell_selector)); diff --git a/deal.II/examples/step-39/step-39.cc b/deal.II/examples/step-39/step-39.cc index 98f87b5099..54097ff146 100644 --- a/deal.II/examples/step-39/step-39.cc +++ b/deal.II/examples/step-39/step-39.cc @@ -97,10 +97,10 @@ 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; + static void cell(typename MeshWorker::IntegrationWorker::CellInfo& info); + static void bdry(typename MeshWorker::IntegrationWorker::FaceInfo& info); + static void face(typename MeshWorker::IntegrationWorker::FaceInfo& info1, + typename MeshWorker::IntegrationWorker::FaceInfo& info2); }; @@ -112,7 +112,7 @@ class MatrixIntegrator : public Subscriptor // degrees of freedom associated // with the shape functions. template -void MatrixIntegrator::cell(typename MeshWorker::IntegrationWorker::CellInfo& info) const +void MatrixIntegrator::cell(typename MeshWorker::IntegrationWorker::CellInfo& info) { const FEValuesBase& fe = info.fe(); FullMatrix& local_matrix = info.M1[0].matrix; @@ -127,7 +127,7 @@ void MatrixIntegrator::cell(typename MeshWorker::IntegrationWorker::Ce // On boundary faces, we use the // Nitsche boundary condition template -void MatrixIntegrator::bdry(typename MeshWorker::IntegrationWorker::FaceInfo& info) const +void MatrixIntegrator::bdry(typename MeshWorker::IntegrationWorker::FaceInfo& info) { const FEFaceValuesBase& fe = info.fe(); FullMatrix& local_matrix = info.M1[0].matrix; @@ -147,7 +147,7 @@ void MatrixIntegrator::bdry(typename MeshWorker::IntegrationWorker::Fa template void MatrixIntegrator::face(typename MeshWorker::IntegrationWorker::FaceInfo& info1, - typename MeshWorker::IntegrationWorker::FaceInfo& info2) const + typename MeshWorker::IntegrationWorker::FaceInfo& info2) { const FEFaceValuesBase& fe1 = info1.fe(); const FEFaceValuesBase& fe2 = info2.fe(); @@ -196,20 +196,20 @@ template 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; + static void cell(typename MeshWorker::IntegrationWorker::CellInfo& info); + static void bdry(typename MeshWorker::IntegrationWorker::FaceInfo& info); + static void face(typename MeshWorker::IntegrationWorker::FaceInfo& info1, + typename MeshWorker::IntegrationWorker::FaceInfo& info2); }; template -void RHSIntegrator::cell(typename MeshWorker::IntegrationWorker::CellInfo&) const +void RHSIntegrator::cell(typename MeshWorker::IntegrationWorker::CellInfo&) {} template -void RHSIntegrator::bdry(typename MeshWorker::IntegrationWorker::FaceInfo& info) const +void RHSIntegrator::bdry(typename MeshWorker::IntegrationWorker::FaceInfo& info) { const FEFaceValuesBase& fe = info.fe(); Vector& local_vector = info.R[0].block(0); @@ -230,7 +230,7 @@ void RHSIntegrator::bdry(typename MeshWorker::IntegrationWorker::FaceI template void RHSIntegrator::face(typename MeshWorker::IntegrationWorker::FaceInfo&, - typename MeshWorker::IntegrationWorker::FaceInfo&) const + typename MeshWorker::IntegrationWorker::FaceInfo&) {} @@ -238,15 +238,15 @@ template class Estimator : 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; + static void cell(typename MeshWorker::IntegrationWorker::CellInfo& info); + static void bdry(typename MeshWorker::IntegrationWorker::FaceInfo& info); + static void face(typename MeshWorker::IntegrationWorker::FaceInfo& info1, + typename MeshWorker::IntegrationWorker::FaceInfo& info2); }; template -void Estimator::cell(typename MeshWorker::IntegrationWorker::CellInfo& info) const +void Estimator::cell(typename MeshWorker::IntegrationWorker::CellInfo& info) { const FEValuesBase& fe = info.fe(); @@ -260,7 +260,7 @@ void Estimator::cell(typename MeshWorker::IntegrationWorker::CellInfo& template -void Estimator::bdry(typename MeshWorker::IntegrationWorker::FaceInfo& info) const +void Estimator::bdry(typename MeshWorker::IntegrationWorker::FaceInfo& info) { const FEFaceValuesBase& fe = info.fe(); @@ -280,7 +280,7 @@ void Estimator::bdry(typename MeshWorker::IntegrationWorker::FaceInfo& template void Estimator::face(typename MeshWorker::IntegrationWorker::FaceInfo& info1, - typename MeshWorker::IntegrationWorker::FaceInfo& info2) const + typename MeshWorker::IntegrationWorker::FaceInfo& info2) { const FEFaceValuesBase& fe = info1.fe(); const std::vector& uh1 = info1.values[0][0]; @@ -310,6 +310,9 @@ template class Step39 { public: + typedef typename MeshWorker::IntegrationWorker::CellInfo CellInfo; + typedef typename MeshWorker::IntegrationWorker::FaceInfo FaceInfo; + Step39(const FiniteElement& fe); void run(unsigned int n_steps); @@ -412,18 +415,24 @@ template void Step39::assemble_matrix() { - const MatrixIntegrator local; - MeshWorker::AssemblingIntegrator >, MatrixIntegrator > - integrator(local); + MeshWorker::IntegrationWorker integration_worker; + MeshWorker::Assembler::MatrixSimple > assembler; + 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); + integration_worker.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); + integration_worker.add_update_flags(update_flags, true, true, true, true); - integrator.initialize(matrix); + assembler.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); + info_box.initialize(integration_worker, assembler, fe, mapping); + MeshWorker::integration_loop( + dof_handler.begin_active(), dof_handler.end(), + info_box, + &MatrixIntegrator::cell, + &MatrixIntegrator::bdry, + &MatrixIntegrator::face, + assembler); } @@ -431,19 +440,25 @@ template void Step39::assemble_mg_matrix() { - const MatrixIntegrator local; - MeshWorker::AssemblingIntegrator >, MatrixIntegrator > - integrator(local); + MeshWorker::IntegrationWorker integration_worker; + MeshWorker::Assembler::MGMatrixSimple > assembler; + const unsigned int n_gauss_points = mg_dof_handler.get_fe().tensor_degree()+1; - integrator.initialize_gauss_quadrature(n_gauss_points, n_gauss_points, n_gauss_points); + integration_worker.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); + integration_worker.add_update_flags(update_flags, true, true, true, true); - integrator.initialize(mg_matrix); - integrator.initialize_fluxes(mg_matrix_dg_up, mg_matrix_dg_down); + assembler.initialize(mg_matrix); + assembler.initialize_fluxes(mg_matrix_dg_up, mg_matrix_dg_down); MeshWorker::IntegrationInfoBox info_box(mg_dof_handler); - info_box.initialize(integrator, fe, mapping); - MeshWorker::integration_loop(mg_dof_handler.begin(), mg_dof_handler.end(), info_box, integrator); + info_box.initialize(integration_worker, assembler, fe, mapping); + MeshWorker::integration_loop( + dof_handler.begin_active(), dof_handler.end(), + info_box, + &MatrixIntegrator::cell, + &MatrixIntegrator::bdry, + &MatrixIntegrator::face, + assembler); } @@ -451,21 +466,28 @@ template void Step39::assemble_right_hand_side() { - const RHSIntegrator local; - MeshWorker::AssemblingIntegrator >, RHSIntegrator > - integrator(local); + MeshWorker::IntegrationWorker integration_worker; + MeshWorker::Assembler::ResidualSimple > assembler; + const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1; - integrator.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points); + integration_worker.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points); UpdateFlags update_flags = update_quadrature_points | update_values | update_gradients; - integrator.add_update_flags(update_flags, true, true, true, true); + integration_worker.add_update_flags(update_flags, true, true, true, true); NamedData* > data; Vector* rhs = &right_hand_side; data.add(rhs, "RHS"); - integrator.initialize(data); + assembler.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); + info_box.initialize(integration_worker, assembler, fe, mapping); + + MeshWorker::integration_loop( + dof_handler.begin_active(), dof_handler.end(), + info_box, + &RHSIntegrator::cell, + &RHSIntegrator::bdry, + &RHSIntegrator::face, + assembler); right_hand_side *= -1.; } @@ -553,23 +575,39 @@ Step39::estimate() cell != triangulation.end();++cell,++i) cell->set_user_index(i); - const Estimator local; - MeshWorker::AssemblingIntegrator, Estimator > - integrator(local); + MeshWorker::IntegrationWorker integration_worker; + MeshWorker::Assembler::CellsAndFaces assembler; + const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1; - integrator.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points); - UpdateFlags update_flags = update_values | update_gradients; - integrator.add_update_flags(update_flags | update_hessians, true, true, true, true); - integrator.add_update_flags(update_quadrature_points, false, true, false, false); + integration_worker.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points); + + NamedData* > solution_data; + solution_data.add(&solution, "solution"); + + MeshWorker::VectorSelector cs; + MeshWorker::VectorSelector fs; + cs.add("solution", true, true, true); + fs.add("solution", true, true, false); + + integration_worker.cell_selector = cs; + integration_worker.boundary_selector = fs; + integration_worker.face_selector = fs; + integration_worker.initialize_update_flags(); + integration_worker.add_update_flags(update_quadrature_points, false, true, false, false); NamedData* > out_data; BlockVector* est = &estimates; out_data.add(est, "cells"); - integrator.initialize(out_data, false); + assembler.initialize(out_data, false); MeshWorker::IntegrationInfoBox info_box(dof_handler); - info_box.initialize_data(&solution, std::string("solution"), true, true, true); - info_box.initialize(integrator, fe, mapping); - MeshWorker::integration_loop(dof_handler.begin_active(), dof_handler.end(), info_box, integrator); + info_box.initialize(integration_worker, assembler, fe, mapping); + MeshWorker::integration_loop ( + dof_handler.begin_active(), dof_handler.end(), + info_box, + &Estimator::cell, + &Estimator::bdry, + &Estimator::face, + assembler); return estimates.block(0).l2_norm(); } -- 2.39.5