From: turcksin Date: Fri, 31 May 2013 23:25:45 +0000 (+0000) Subject: Merge from mainline. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=cf32f216a96b19d6225572cf9908d4e18d063b67;p=dealii-svn.git Merge from mainline. git-svn-id: https://svn.dealii.org/branches/branch_bigger_global_dof_indices_4@29705 0785d39b-7218-0410-832d-ea1e28bc413d --- cf32f216a96b19d6225572cf9908d4e18d063b67 diff --cc tests/integrators/elasticity_01.cc index e8129c1ad2,22bf91a8fa..78f25e6c14 --- a/tests/integrators/elasticity_01.cc +++ b/tests/integrators/elasticity_01.cc @@@ -38,23 -43,25 +43,25 @@@ void test_cell(const FEValuesBase std::vector > > ugrad(dim,std::vector >(fev.n_quadrature_points)); - std::vector indices(n); + std::vector indices(n); for (unsigned int i=0;i Vector u2(n2), v2(n2), w2(n2); std::vector > u1val (d,std::vector(fev1.n_quadrature_points)), - u2val (d,std::vector(fev2.n_quadrature_points)); + nullval (d,std::vector(fev2.n_quadrature_points, 0.)); std::vector > > u1grad(d,std::vector >(fev1.n_quadrature_points)), - u2grad(d,std::vector >(fev2.n_quadrature_points)); + nullgrad(d,std::vector >(fev2.n_quadrature_points, Tensor<1,dim>())); - std::vector indices1(n1), indices2(n2); + std::vector indices1(n1), indices2(n2); for (unsigned int i=0;i std::vector > > ugrad(d,std::vector >(fev.n_quadrature_points)); - std::vector indices(n); + std::vector indices(n); for (unsigned int i=0;i Vector u2(n2), v2(n2), w2(n2); std::vector > u1val (d,std::vector(fev1.n_quadrature_points)), - u2val (d,std::vector(fev2.n_quadrature_points)); + nullval (d,std::vector(fev2.n_quadrature_points, 0.)); std::vector > > u1grad(d,std::vector >(fev1.n_quadrature_points)), - u2grad(d,std::vector >(fev2.n_quadrature_points)); + nullgrad(d,std::vector >(fev2.n_quadrature_points, Tensor<1,dim>())); - std::vector indices1(n1), indices2(n2); + std::vector indices1(n1), indices2(n2); for (unsigned int i=0;i + #include + #include + #include + #include + #include + + #include + #include + + #include + #include + #include + #include + #include + + #include + #include + #include + #include + + #include + + #include + #include + #include + #include + #include + #include + + #include + #include + #include + #include + + #include + #include + + namespace Step39 + { + using namespace dealii; + + Functions::SlitSingularityFunction<2> exact_solution; + + + + + template + class MatrixIntegrator : public MeshWorker::LocalIntegrator + { + public: + void cell(MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const; + void boundary(MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const; + void face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const; + }; + + + template + void MatrixIntegrator::cell( + MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const + { + LocalIntegrators::Laplace::cell_matrix(dinfo.matrix(0,false).matrix, info.fe_values()); + } + + + template + void MatrixIntegrator::boundary( + MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const + { + const unsigned int deg = info.fe_values(0).get_fe().tensor_degree(); + LocalIntegrators::Laplace::nitsche_matrix( + dinfo.matrix(0,false).matrix, info.fe_values(0), + LocalIntegrators::Laplace::compute_penalty(dinfo, dinfo, deg, deg)); + } + + template + void MatrixIntegrator::face( + MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const + { + const unsigned int deg = info1.fe_values(0).get_fe().tensor_degree(); + LocalIntegrators::Laplace::ip_matrix( + dinfo1.matrix(0,false).matrix, dinfo1.matrix(0,true).matrix, + dinfo2.matrix(0,true).matrix, dinfo2.matrix(0,false).matrix, + info1.fe_values(0), info2.fe_values(0), + LocalIntegrators::Laplace::compute_penalty(dinfo1, dinfo2, deg, deg)); + } + + template + class RHSIntegrator : public MeshWorker::LocalIntegrator + { + public: + void cell(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const; + }; + + + template + void RHSIntegrator::cell(MeshWorker::DoFInfo &, typename MeshWorker::IntegrationInfo &) const + {} + + + template + void RHSIntegrator::boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const + { + const FEValuesBase &fe = info.fe_values(); + Vector &local_vector = dinfo.vector(0).block(0); + + std::vector boundary_values(fe.n_quadrature_points); + exact_solution.value_list(fe.get_quadrature_points(), boundary_values); + + const unsigned int deg = fe.get_fe().tensor_degree(); + const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure(); + + for (unsigned k=0; k + void RHSIntegrator::face(MeshWorker::DoFInfo &, + MeshWorker::DoFInfo &, + typename MeshWorker::IntegrationInfo &, + typename MeshWorker::IntegrationInfo &) const + {} + + + template + class Estimator : public MeshWorker::LocalIntegrator + { + public: + void cell(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const; + }; + + + template + void Estimator::cell(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const + { + const FEValuesBase &fe = info.fe_values(); + + const std::vector > &DDuh = info.hessians[0][0]; + for (unsigned k=0; kdiameter() * trace(DDuh[k]); + dinfo.value(0) += t*t * fe.JxW(k); + } + dinfo.value(0) = std::sqrt(dinfo.value(0)); + } + + template + void Estimator::boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const + { + const FEValuesBase &fe = info.fe_values(); + + 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) * dinfo.face->measure() / dinfo.cell->measure(); + + for (unsigned k=0; k + void Estimator::face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const + { + const FEValuesBase &fe = info1.fe_values(); + 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) * dinfo1.face->measure() / dinfo1.cell->measure(); + const double penalty2 = deg * (deg+1) * dinfo2.face->measure() / dinfo2.cell->measure(); + const double penalty = penalty1 + penalty2; + const double h = dinfo1.face->measure(); + + for (unsigned k=0; k + class ErrorIntegrator : public MeshWorker::LocalIntegrator + { + public: + void cell(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const; + }; + + + template + void ErrorIntegrator::cell( + MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const + { + const FEValuesBase &fe = info.fe_values(); + std::vector > exact_gradients(fe.n_quadrature_points); + std::vector exact_values(fe.n_quadrature_points); + + exact_solution.gradient_list(fe.get_quadrature_points(), exact_gradients); + exact_solution.value_list(fe.get_quadrature_points(), exact_values); + + const std::vector > &Duh = info.gradients[0][0]; + const std::vector &uh = info.values[0][0]; + + for (unsigned k=0; k + void ErrorIntegrator::boundary( + MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const + { + const FEValuesBase &fe = info.fe_values(); + + std::vector exact_values(fe.n_quadrature_points); + exact_solution.value_list(fe.get_quadrature_points(), exact_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) * dinfo.face->measure() / dinfo.cell->measure(); + + for (unsigned k=0; k + void ErrorIntegrator::face( + MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const + { + const FEValuesBase &fe = info1.fe_values(); + const std::vector &uh1 = info1.values[0][0]; + const std::vector &uh2 = info2.values[0][0]; + + const unsigned int deg = fe.get_fe().tensor_degree(); + const double penalty1 = deg * (deg+1) * dinfo1.face->measure() / dinfo1.cell->measure(); + const double penalty2 = deg * (deg+1) * dinfo2.face->measure() / dinfo2.cell->measure(); + const double penalty = penalty1 + penalty2; + + for (unsigned k=0; k + class InteriorPenaltyProblem + { + public: + typedef MeshWorker::IntegrationInfo CellInfo; + + InteriorPenaltyProblem(const FiniteElement &fe); + + void run(unsigned int n_steps); + + private: + void setup_system (); + void assemble_matrix (); + void assemble_mg_matrix (); + void assemble_right_hand_side (); + void error (); + double estimate (); + void solve (); + void output_results (const unsigned int cycle) const; + + Triangulation triangulation; + const MappingQ1 mapping; + const FiniteElement &fe; + MGDoFHandler mg_dof_handler; + DoFHandler &dof_handler; + + SparsityPattern sparsity; + SparseMatrix matrix; + Vector solution; + Vector right_hand_side; + BlockVector estimates; + + MGLevelObject mg_sparsity; + MGLevelObject > mg_matrix; + + MGLevelObject mg_sparsity_dg_interface; + MGLevelObject > mg_matrix_dg_down; + MGLevelObject > mg_matrix_dg_up; + MGLevelObject > mg_matrix_in_out; + }; + + + template + InteriorPenaltyProblem::InteriorPenaltyProblem(const FiniteElement &fe) + : + mapping(), + fe(fe), + mg_dof_handler(triangulation), + dof_handler(mg_dof_handler), + estimates(1) + { + GridGenerator::hyper_cube_slit(triangulation, -1, 1); + } + + + template + void + InteriorPenaltyProblem::setup_system() + { + dof_handler.distribute_dofs(fe); - unsigned int n_dofs = dof_handler.n_dofs(); ++ types::global_dof_index n_dofs = dof_handler.n_dofs(); + solution.reinit(n_dofs); + right_hand_side.reinit(n_dofs); + + CompressedSparsityPattern c_sparsity(n_dofs); + DoFTools::make_flux_sparsity_pattern(dof_handler, c_sparsity); + sparsity.copy_from(c_sparsity); + matrix.reinit(sparsity); + + const unsigned int n_levels = triangulation.n_levels(); + mg_matrix.resize(0, n_levels-1); + mg_matrix.clear(); + mg_matrix_dg_up.resize(0, n_levels-1); + mg_matrix_dg_up.clear(); + mg_matrix_dg_down.resize(0, n_levels-1); + mg_matrix_dg_down.clear(); + mg_matrix_in_out.resize(0, n_levels-1); + mg_matrix_in_out.clear(); + mg_sparsity.resize(0, n_levels-1); + mg_sparsity_dg_interface.resize(0, n_levels-1); + + for (unsigned int level=mg_sparsity.min_level(); + level<=mg_sparsity.max_level(); ++level) + { + CompressedSparsityPattern c_sparsity(mg_dof_handler.n_dofs(level)); + MGTools::make_flux_sparsity_pattern(mg_dof_handler, c_sparsity, level); + mg_sparsity[level].copy_from(c_sparsity); + mg_matrix[level].reinit(mg_sparsity[level]); + mg_matrix_in_out[level].reinit(mg_sparsity[level]); + + if (level>0) + { + CompressedSparsityPattern ci_sparsity; + ci_sparsity.reinit(mg_dof_handler.n_dofs(level-1), mg_dof_handler.n_dofs(level)); + MGTools::make_flux_sparsity_pattern_edge(mg_dof_handler, ci_sparsity, level); + mg_sparsity_dg_interface[level].copy_from(ci_sparsity); + mg_matrix_dg_up[level].reinit(mg_sparsity_dg_interface[level]); + mg_matrix_dg_down[level].reinit(mg_sparsity_dg_interface[level]); + } + } + } + + + template + void + InteriorPenaltyProblem::assemble_matrix() + { + MeshWorker::IntegrationInfoBox info_box; + UpdateFlags update_flags = update_values | update_gradients; + info_box.add_update_flags_all(update_flags); + info_box.initialize(fe, mapping); + + MeshWorker::DoFInfo dof_info(dof_handler); + + MeshWorker::Assembler::MatrixSimple > assembler; + assembler.initialize(matrix); + + MatrixIntegrator integrator; + MeshWorker::integration_loop( + dof_handler.begin_active(), dof_handler.end(), + dof_info, info_box, + integrator, assembler); + } + + + template + void + InteriorPenaltyProblem::assemble_mg_matrix() + { + MeshWorker::IntegrationInfoBox info_box; + UpdateFlags update_flags = update_values | update_gradients; + info_box.add_update_flags_all(update_flags); + info_box.initialize(fe, mapping); + + MeshWorker::DoFInfo dof_info(mg_dof_handler); + + MeshWorker::Assembler::MGMatrixSimple > assembler; + assembler.initialize(mg_matrix); + assembler.initialize_interfaces(mg_matrix_in_out, mg_matrix_in_out); + assembler.initialize_fluxes(mg_matrix_dg_up, mg_matrix_dg_down); + + MatrixIntegrator integrator; + MeshWorker::integration_loop ( + mg_dof_handler.begin(), mg_dof_handler.end(), + dof_info, info_box, + integrator, assembler); + + for (unsigned int level=mg_matrix_in_out.min_level(); + level<=mg_matrix_in_out.min_level(); ++level) + if (mg_matrix_in_out[level].frobenius_norm() != 0.) + deallog << "Oops!" << std::endl; + } + + + template + void + InteriorPenaltyProblem::assemble_right_hand_side() + { + MeshWorker::IntegrationInfoBox info_box; + UpdateFlags update_flags = update_quadrature_points | update_values | update_gradients; + info_box.add_update_flags_all(update_flags); + info_box.initialize(fe, mapping); + + MeshWorker::DoFInfo dof_info(dof_handler); + + MeshWorker::Assembler::ResidualSimple > assembler; + NamedData* > data; + Vector *rhs = &right_hand_side; + data.add(rhs, "RHS"); + assembler.initialize(data); + + RHSIntegrator integrator; + MeshWorker::integration_loop( + dof_handler.begin_active(), dof_handler.end(), + dof_info, info_box, + integrator, assembler); + + right_hand_side *= -1.; + } + + + template + void + InteriorPenaltyProblem::solve() + { + SolverControl control(1000, 1.e-12); + SolverCG > solver(control); + + MGTransferPrebuilt > mg_transfer; + mg_transfer.build_matrices(mg_dof_handler); + + FullMatrix coarse_matrix; + coarse_matrix.copy_from (mg_matrix[0]); + MGCoarseGridHouseholder > mg_coarse; + mg_coarse.initialize(coarse_matrix); + + GrowingVectorMemory > mem; + typedef PreconditionSOR > RELAXATION; + mg::SmootherRelaxation > + mg_smoother; + RELAXATION::AdditionalData smoother_data(1.); + mg_smoother.initialize(mg_matrix, smoother_data); + + mg_smoother.set_steps(2); + mg_smoother.set_symmetric(true); + mg_smoother.set_variable(false); + + MGMatrix, Vector > mgmatrix(&mg_matrix); + MGMatrix, Vector > mgdown(&mg_matrix_dg_down); + MGMatrix, Vector > mgup(&mg_matrix_dg_up); + MGMatrix, Vector > mgedge(&mg_matrix_in_out); + + Multigrid > mg(mg_dof_handler, mgmatrix, + mg_coarse, mg_transfer, + mg_smoother, mg_smoother); + mg.set_edge_flux_matrices(mgdown, mgup); + mg.set_edge_matrices(mgedge, mgedge); + + PreconditionMG, + MGTransferPrebuilt > > + preconditioner(mg_dof_handler, mg, mg_transfer); + solver.solve(matrix, solution, right_hand_side, preconditioner); + } + + + template + double + InteriorPenaltyProblem::estimate() + { + std::vector old_user_indices; + triangulation.save_user_indices(old_user_indices); + + 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); + + MeshWorker::IntegrationInfoBox info_box; + const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1; + info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points); + + NamedData* > solution_data; + solution_data.add(&solution, "solution"); + + info_box.cell_selector.add("solution", false, false, true); + info_box.boundary_selector.add("solution", true, true, false); + info_box.face_selector.add("solution", true, true, false); + + info_box.add_update_flags_boundary(update_quadrature_points); + info_box.initialize(fe, mapping, solution_data); + + MeshWorker::DoFInfo dof_info(dof_handler); + + MeshWorker::Assembler::CellsAndFaces assembler; + NamedData* > out_data; + BlockVector *est = &estimates; + out_data.add(est, "cells"); + assembler.initialize(out_data, false); + + Estimator integrator; + MeshWorker::integration_loop ( + dof_handler.begin_active(), dof_handler.end(), + dof_info, info_box, + integrator, assembler); + + triangulation.load_user_indices(old_user_indices); + return estimates.block(0).l2_norm(); + } + + + template + void + InteriorPenaltyProblem::error() + { + BlockVector errors(2); + errors.block(0).reinit(triangulation.n_active_cells()); + errors.block(1).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); + + MeshWorker::IntegrationInfoBox info_box; + const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1; + info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points); + + NamedData* > solution_data; + solution_data.add(&solution, "solution"); + + info_box.cell_selector.add("solution", true, true, false); + info_box.boundary_selector.add("solution", true, false, false); + info_box.face_selector.add("solution", true, false, false); + + info_box.add_update_flags_cell(update_quadrature_points); + info_box.add_update_flags_boundary(update_quadrature_points); + info_box.initialize(fe, mapping, solution_data); + + MeshWorker::DoFInfo dof_info(dof_handler); + + MeshWorker::Assembler::CellsAndFaces assembler; + NamedData* > out_data; + BlockVector *est = &errors; + out_data.add(est, "cells"); + assembler.initialize(out_data, false); + + ErrorIntegrator integrator; + MeshWorker::integration_loop ( + dof_handler.begin_active(), dof_handler.end(), + dof_info, info_box, + integrator, assembler); + + deallog << "energy-error: " << errors.block(0).l2_norm() << std::endl; + deallog << "L2-error: " << errors.block(1).l2_norm() << std::endl; + } + + + template + void InteriorPenaltyProblem::output_results (const unsigned int cycle) const + { + char *fn = new char[100]; + sprintf(fn, "sol-%02d", cycle); + + std::string filename(fn); + filename += ".gnuplot"; + deallog << "Writing solution to <" << filename << ">..." + << std::endl << std::endl; + std::ofstream gnuplot_output (filename.c_str()); + + 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 (); + + data_out.write_gnuplot(gnuplot_output); + } + + template + void + InteriorPenaltyProblem::run(unsigned int n_steps) + { + deallog << "Element: " << fe.get_name() << std::endl; + for (unsigned int s=0; s fe1(FE_DGQ<2>(2), 1, FE_Q<2>(2), 1); + InteriorPenaltyProblem<2> test1(fe1); + test1.run(6); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; + }