From: Guido Kanschat Date: Wed, 2 Dec 2009 04:12:15 +0000 (+0000) Subject: step 39 running X-Git-Tag: v8.0.0~6767 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=710a1127ec97a953593913faec69b56ba1bafc1b;p=dealii.git step 39 running git-svn-id: https://svn.dealii.org/trunk@20195 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/numerics/mesh_worker_assembler.h b/deal.II/deal.II/include/numerics/mesh_worker_assembler.h index 33af1b0de0..b0077fd09a 100644 --- a/deal.II/deal.II/include/numerics/mesh_worker_assembler.h +++ b/deal.II/deal.II/include/numerics/mesh_worker_assembler.h @@ -160,7 +160,7 @@ namespace MeshWorker * communicating the cell and * face vectors. */ - typedef NamedData,CellsAndFaces > > DataVectors; + typedef NamedData*> DataVectors; /** * The initialization @@ -219,7 +219,7 @@ namespace MeshWorker */ template void initialize_info(DoFInfo& info, - bool interior_face); + bool interior_face) const; /** * Assemble the local values @@ -962,6 +962,116 @@ namespace MeshWorker }; +//----------------------------------------------------------------------// + + template + inline void + Functional::initialize(unsigned int n) + { + results.resize(n); + } + + + template + template + inline void + Functional::assemble(const DoFInfo& info) + { + for (unsigned int i=0;i + template + inline void + Functional::assemble(const DoFInfo& info1, + const DoFInfo& info2) + { + for (unsigned int i=0;i + inline number + Functional::operator() (unsigned int i) const + { + AssertIndexRange(i, results.size()); + return results[i]; + } + +//----------------------------------------------------------------------// + + template + inline void + CellsAndFaces::initialize(DataVectors& r, bool sep) + { + Assert(r.name(0) == "cells", typename DataVectors::ExcNameMismatch(0, "cells")); + if (sep) + { + Assert(r.name(1) == "faces", typename DataVectors::ExcNameMismatch(1, "faces")); + AssertDimension(r(0)->n_blocks(), r(1)->n_blocks()); + } + + results = r; + separate_faces = sep; + } + + + template + template + inline void + CellsAndFaces::initialize_info(DoFInfo& info, bool) const + { + info.initialize_numbers(results(0)->n_blocks()); + } + + + template + template + inline void + CellsAndFaces::assemble(const DoFInfo& info) + { + for (unsigned int i=0;iblock(i)(info.face->user_index()) += info.J[i]; + else + results(0)->block(i)(info.cell->user_index()) += info.J[i]; + } + } + + + template + template + inline void + CellsAndFaces::assemble(const DoFInfo& info1, + const DoFInfo& info2) + { + for (unsigned int i=0;iblock(i)(info1.face->user_index()) += J; + if (info2.face != info1.face) + results(1)->block(i)(info2.face->user_index()) += J; + } + else + { + results(0)->block(i)(info1.cell->user_index()) += .5*info1.J[i]; + results(0)->block(i)(info2.cell->user_index()) += .5*info2.J[i]; + } + } + } + + + //----------------------------------------------------------------------// template 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 947a66b23a..cef55748d3 100644 --- a/deal.II/deal.II/include/numerics/mesh_worker_info.h +++ b/deal.II/deal.II/include/numerics/mesh_worker_info.h @@ -47,7 +47,8 @@ namespace MeshWorker MatrixBlock >& M, unsigned int row, unsigned int col); public: - void initialize_vectors(const unsigned int n_vectors); + void initialize_numbers(const unsigned int n); + void initialize_vectors(const unsigned int n); /** * Allocate @p n local * matrices. Additionally, @@ -542,19 +543,20 @@ namespace MeshWorker */ template IntegrationInfoBox(const T&); - + + template + void initialize_data(const VECTOR*, const std::string& name, + bool values, bool gradients, bool hessians); + template void initialize(const WORKER&, const FiniteElement& el, const Mapping& mapping); - template - void initialize(const WORKER&, - const FiniteElement& el, - const Mapping& mapping, - const NamedData >& data); // private: - + + boost::shared_ptr > data; + CellInfo cell_info; FaceInfo bdry_info; FaceInfo face_info; @@ -565,6 +567,14 @@ namespace MeshWorker //----------------------------------------------------------------------// + template + inline void + LocalResults::initialize_numbers(unsigned int n) + { + J.resize(n); + } + + template inline void LocalResults::initialize_vectors(unsigned int n) @@ -623,6 +633,8 @@ namespace MeshWorker inline void LocalResults::reinit(const BlockIndices& bi) { + for (unsigned int i=0;i + template + void + IntegrationInfoBox::initialize_data( + const VECTOR* v, const std::string& name, + bool values, bool gradients, bool hessians) + { + boost::shared_ptr > + p = boost::shared_ptr >(new VectorData ()); + p->add(name, values, gradients, hessians); + p->initialize(v, name); + data = p; + cell_info.initialize_data(data); + bdry_info.initialize_data(data); + face_info.initialize_data(data); + subface_info.initialize_data(data); + neighbor_info.initialize_data(data); + } + + template template void @@ -836,33 +870,14 @@ namespace MeshWorker cell_info.initialize >(el, mapping, integrator.cell_quadrature, integrator.cell_flags); bdry_info.initialize >(el, mapping, integrator.bdry_quadrature, - integrator.face_flags); + integrator.bdry_flags); face_info.initialize >(el, mapping, integrator.face_quadrature, integrator.face_flags); subface_info.initialize >(el, mapping, integrator.face_quadrature, integrator.face_flags); neighbor_info.initialize >(el, mapping, integrator.face_quadrature, integrator.ngbr_flags); - } - - - template - template - void - IntegrationInfoBox::initialize( - const WORKER& integrator, - const FiniteElement& el, - const Mapping& mapping, - const NamedData >& data) - { - cell_info.initialize_data(data); - bdry_info.initialize_data(data); - face_info.initialize_data(data); - subface_info.initialize_data(data); - neighbor_info.initialize_data(data); - - initialize(integrator, el, mapping); - } + } } diff --git a/deal.II/deal.II/include/numerics/mesh_worker_info.templates.h b/deal.II/deal.II/include/numerics/mesh_worker_info.templates.h index c92dbc15e5..95c5274699 100644 --- a/deal.II/deal.II/include/numerics/mesh_worker_info.templates.h +++ b/deal.II/deal.II/include/numerics/mesh_worker_info.templates.h @@ -106,7 +106,7 @@ namespace MeshWorker // For all components for (unsigned int j=0;j void print (STREAM& s, const NamedData& v) const; + /** + * Print the number of selections to the stream. + */ + template + void print (STREAM& s) const; + protected: /** * Selection of the vectors @@ -238,9 +244,11 @@ namespace MeshWorker unsigned int size) const; }; + /** * Based on VectorSelector, this is the class that implements the - * function VectorDataBase::fill() for a certain type of vector. + * function VectorDataBase::fill() for a certain type of vector, using + * NamedData to identify vectors by name. * * @author Guido Kanschat, 2009 */ @@ -250,6 +258,7 @@ namespace MeshWorker { public: void initialize(const NamedData&); + void initialize(const VECTOR*, const std::string& name); virtual void fill( std::vector > >& values, @@ -262,8 +271,9 @@ namespace MeshWorker unsigned int start, unsigned int size) const; private: - NamedData > > data; + NamedData > > data; }; + //----------------------------------------------------------------------// @@ -357,6 +367,17 @@ namespace MeshWorker } + template + inline void + VectorSelector::print(STREAM& s) const + { + s << "values: " << n_values() + << " gradients: " << n_gradients() + << " hessians: " << n_hessians() + << std::endl; + } + + template inline void VectorSelector::print(STREAM& s, const NamedData& v) const diff --git a/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.templates.h b/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.templates.h index d0d3e6a6df..e691467949 100644 --- a/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.templates.h +++ b/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.templates.h @@ -47,6 +47,16 @@ namespace MeshWorker } + template + void + VectorData::initialize(const VECTOR* v, const std::string& name) + { + SmartPointer > p = v; + data.add(p, name); + VectorSelector::initialize(data); + } + + template void VectorData::fill( diff --git a/deal.II/deal.II/source/numerics/mesh_worker_assembler.all_dimensions.cc b/deal.II/deal.II/source/numerics/mesh_worker_assembler.all_dimensions.cc index c57593172c..b6c27e8d8d 100644 --- a/deal.II/deal.II/source/numerics/mesh_worker_assembler.all_dimensions.cc +++ b/deal.II/deal.II/source/numerics/mesh_worker_assembler.all_dimensions.cc @@ -22,107 +22,9 @@ namespace MeshWorker namespace Assembler { //----------------------------------------------------------------------// - - template - inline void - Functional::initialize(unsigned int n) - { - results.resize(n); - } - - - template - template - inline void - Functional::assemble(const DoFInfo& info) - { - for (unsigned int i=0;i - template - inline void - Functional::assemble(const DoFInfo& info1, - const DoFInfo& info2) - { - for (unsigned int i=0;i - inline number - Functional::operator() (unsigned int i) const - { - AssertIndexRange(i, results.size()); - return results[i]; - } - //----------------------------------------------------------------------// - template - inline void - CellsAndFaces::initialize(DataVectors& r, bool sep) - { - Assert(r.name(0) == "cells", DataVectors::ExcNameMismatch(0, "cells")); - if (sep) - Assert(r.name(1) == "faces", DataVectors::ExcNameMismatch(1, "faces")); - AssertDimension(r(0).n_blocks(), r(1).n_blocks()); - - results = r; - separate_faces = sep; - } - - - template - template - inline void - CellsAndFaces::assemble(const DoFInfo& info) - { - for (unsigned int i=0;iuser_index()) += info.J[i]; - else - results.vector(0).block(i)(info.cell->user_index()) += info.J[i]; - } - } - - - template - template - inline void - CellsAndFaces::assemble(const DoFInfo& info1, - const DoFInfo& info2) - { - for (unsigned int i=0;iuser_index()) += J; - if (info2.face != info1.face) - results.vector(1).block(i)(info2.face->user_index()) += J; - } - else - { - results.vector(0).block(i)(info1.cell->user_index()) += info1.J[i]; - results.vector(0).block(i)(info2.cell->user_index()) += info2.J[i]; - } - } - } - - -//----------------------------------------------------------------------// - - - } } diff --git a/deal.II/examples/step-39/step-39.cc b/deal.II/examples/step-39/step-39.cc index f17955777c..2961469b46 100644 --- a/deal.II/examples/step-39/step-39.cc +++ b/deal.II/examples/step-39/step-39.cc @@ -20,10 +20,12 @@ #include #include #include +#include // Include files for setting up the // mesh #include +#include // Include files for FiniteElement // classes and DoFHandler. @@ -102,6 +104,13 @@ class MatrixIntegrator : public Subscriptor }; + // On each cell, we integrate the + // Dirichlet form. All local + // integrations consist of nested + // loops, first over all quadrature + // points and the iner loops over the + // degrees of freedom associated + // with the shape functions. template void MatrixIntegrator::cell(typename MeshWorker::IntegrationWorker::CellInfo& info) const { @@ -115,7 +124,8 @@ void MatrixIntegrator::cell(typename MeshWorker::IntegrationWorker::Ce * fe.JxW(k); } - + // On boundary faces, we use the + // Nitsche boundary condition template void MatrixIntegrator::bdry(typename MeshWorker::IntegrationWorker::FaceInfo& info) const { @@ -153,7 +163,7 @@ void MatrixIntegrator::face(typename MeshWorker::IntegrationWorker::Fa { Assert (info1.face == info2.face, ExcInternalError()); Assert (info1.face->has_children(), ExcInternalError()); - Assert (info1.cell->has_children(), ExcInternalError()); +// Assert (info1.cell->has_children(), ExcInternalError()); penalty1 *= 2; } const double penalty = penalty1 + penalty2; @@ -224,6 +234,77 @@ void RHSIntegrator::face(typename MeshWorker::IntegrationWorker::FaceI {} +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; +}; + + +template +void Estimator::cell(typename MeshWorker::IntegrationWorker::CellInfo& info) const +{ + const FEValuesBase& fe = info.fe(); + + const std::vector >& DDuh = info.hessians[0][0]; + for (unsigned k=0;kdiameter() * trace(DDuh[k]); + info.J[0] += t*t * fe.JxW(k); + } +} + + +template +void Estimator::bdry(typename MeshWorker::IntegrationWorker::FaceInfo& info) const +{ + const FEFaceValuesBase& fe = info.fe(); + + std::vector boundary_values(fe.n_quadrature_points); + exact_solution.value_list(fe.get_quadrature_points(), boundary_values); + + const std::vector& uh = info.values[0][0]; + + 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 Estimator::face(typename MeshWorker::IntegrationWorker::FaceInfo& info1, + typename MeshWorker::IntegrationWorker::FaceInfo& info2) const +{ + const FEFaceValuesBase& fe = info1.fe(); + const std::vector& uh1 = info1.values[0][0]; + const std::vector& uh2 = info2.values[0][0]; + const std::vector >& Duh1 = info1.gradients[0][0]; + const std::vector >& Duh2 = info2.gradients[0][0]; + + 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; + const double h = info1.face->measure(); + + for (unsigned k=0;k class Step39 @@ -239,7 +320,7 @@ class Step39 void assemble_mg_matrix (); void assemble_right_hand_side (); void error (); - void estimate (); + double estimate (); void solve (); void refine_grid (); void output_results (const unsigned int cycle) const; @@ -254,7 +335,8 @@ class Step39 SparseMatrix matrix; Vector solution; Vector right_hand_side; - + BlockVector estimates; + MGLevelObject mg_sparsity; MGLevelObject mg_sparsity_dg_interface; MGLevelObject > mg_matrix; @@ -268,11 +350,10 @@ Step39::Step39(const FiniteElement& fe) : fe(fe), mg_dof_handler(triangulation), - dof_handler(mg_dof_handler) + dof_handler(mg_dof_handler), + estimates(1) { GridGenerator::hyper_L(triangulation, -1, 1); - triangulation.begin()->set_refine_flag(); - triangulation.execute_coarsening_and_refinement(); } @@ -474,7 +555,39 @@ Step39::error() QGauss quadrature(n_gauss_points); VectorTools::integrate_difference(mapping, dof_handler, solution, exact_solution, cell_errors, quadrature, VectorTools::L2_norm); - deallog << "Error " << cell_errors.l2_norm() << std::endl; + deallog << "Error " << cell_errors.l2_norm() << std::endl; +} + + +template +double +Step39::estimate() +{ + estimates.block(0).reinit(triangulation.n_active_cells()); + unsigned int i=0; + for (typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); + cell != triangulation.end();++cell,++i) + cell->set_user_index(i); + + const Estimator local; + MeshWorker::AssemblingIntegrator, Estimator > + 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+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_hessians, true, false, false, false); + integrator.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); + 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); + return estimates.block(0).l2_norm(); } @@ -483,10 +596,10 @@ void Step39::output_results (const unsigned int cycle) const { // Output of the solution in // gnuplot format. - std::string filename = "sol-"; - filename += ('0' + cycle); - Assert (cycle < 10, ExcInternalError()); - + char * fn = new char[100]; + sprintf(fn, "sol-%02d", cycle); + + std::string filename(fn); filename += ".gnuplot"; std::cout << "Writing solution to <" << filename << ">..." << std::endl << std::endl; @@ -495,6 +608,7 @@ void Step39::output_results (const unsigned int cycle) const DataOut data_out; data_out.attach_dof_handler (dof_handler); data_out.add_data_vector (solution, "u"); + data_out.add_data_vector (estimates.block(0), "est"); data_out.build_patches (); @@ -511,7 +625,17 @@ Step39::run(unsigned int n_steps) { deallog << "Step " << s << std::endl; if (s != 0) - triangulation.refine_global(1); + { + if (estimates.block(0).size() == 0) + triangulation.refine_global(1); + else + { + GridRefinement::refine_and_coarsen_fixed_fraction (triangulation, + estimates.block(0), + 0.5, 0.0); + triangulation.execute_coarsening_and_refinement (); + } + } deallog << "Triangulation " << triangulation.n_active_cells() << " cells, " @@ -531,8 +655,8 @@ Step39::run(unsigned int n_steps) assemble_right_hand_side(); deallog << "Solve" << std::endl; solve(); - deallog << "Error" << std::endl; error(); + deallog << "Estimate " << estimate() << std::endl; output_results(s); } } @@ -540,7 +664,7 @@ Step39::run(unsigned int n_steps) int main() { - FE_DGQ<2> fe1(1); + FE_DGQ<2> fe1(3); Step39<2> test1(fe1); - test1.run(7); + test1.run(13); }