-/* Author: Guido Kanschat, Texas A&M University, 2009 */
-
-/* $Id: step-39.cc 28552 2013-02-25 15:18:15Z kanschat $ */
-/* */
-/* Copyright (C) 2010-2013 by the deal.II authors */
-/* */
-/* This file is subject to QPL and may not be distributed */
-/* without copyright and license information. Please refer */
-/* to the file deal.II/doc/license.html for the text and */
-/* further information on this license. */
+//----------------------------------------------------------------------
+// $Id: mesh_worker_01.cc 28614 2013-02-27 15:20:34Z kanschat $
+//
+// Copyright (C) 2013 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//----------------------------------------------------------------------
+
+// Test problem in MGTools::make_flux_sparsity_pattern_edge
#include "../tests.h"
#include <deal.II/lac/sparse_matrix.h>
{
using namespace dealii;
- Functions::SlitSingularityFunction<2> exact_solution;
-
-
-
-
- template <int dim>
- class MatrixIntegrator : public MeshWorker::LocalIntegrator<dim>
- {
- public:
- void cell(MeshWorker::DoFInfo<dim> &dinfo,
- typename MeshWorker::IntegrationInfo<dim> &info) const;
- void boundary(MeshWorker::DoFInfo<dim> &dinfo,
- typename MeshWorker::IntegrationInfo<dim> &info) const;
- void face(MeshWorker::DoFInfo<dim> &dinfo1,
- MeshWorker::DoFInfo<dim> &dinfo2,
- typename MeshWorker::IntegrationInfo<dim> &info1,
- typename MeshWorker::IntegrationInfo<dim> &info2) const;
- };
-
-
- template <int dim>
- void MatrixIntegrator<dim>::cell(
- MeshWorker::DoFInfo<dim> &dinfo,
- typename MeshWorker::IntegrationInfo<dim> &info) const
- {
- LocalIntegrators::Laplace::cell_matrix(dinfo.matrix(0,false).matrix, info.fe_values());
- }
-
-
- template <int dim>
- void MatrixIntegrator<dim>::boundary(
- MeshWorker::DoFInfo<dim> &dinfo,
- typename MeshWorker::IntegrationInfo<dim> &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 <int dim>
- void MatrixIntegrator<dim>::face(
- MeshWorker::DoFInfo<dim> &dinfo1,
- MeshWorker::DoFInfo<dim> &dinfo2,
- typename MeshWorker::IntegrationInfo<dim> &info1,
- typename MeshWorker::IntegrationInfo<dim> &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 <int dim>
- class RHSIntegrator : public MeshWorker::LocalIntegrator<dim>
- {
- public:
- void cell(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const;
- void boundary(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const;
- void face(MeshWorker::DoFInfo<dim> &dinfo1,
- MeshWorker::DoFInfo<dim> &dinfo2,
- typename MeshWorker::IntegrationInfo<dim> &info1,
- typename MeshWorker::IntegrationInfo<dim> &info2) const;
- };
-
-
- template <int dim>
- void RHSIntegrator<dim>::cell(MeshWorker::DoFInfo<dim> &, typename MeshWorker::IntegrationInfo<dim> &) const
- {}
-
-
- template <int dim>
- void RHSIntegrator<dim>::boundary(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const
- {
- const FEValuesBase<dim> &fe = info.fe_values();
- Vector<double> &local_vector = dinfo.vector(0).block(0);
-
- std::vector<double> 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<fe.n_quadrature_points; ++k)
- for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- local_vector(i) += (- fe.shape_value(i,k) * penalty * boundary_values[k]
- + (fe.normal_vector(k) * fe.shape_grad(i,k)) * boundary_values[k])
- * fe.JxW(k);
- }
-
-
- template <int dim>
- void RHSIntegrator<dim>::face(MeshWorker::DoFInfo<dim> &,
- MeshWorker::DoFInfo<dim> &,
- typename MeshWorker::IntegrationInfo<dim> &,
- typename MeshWorker::IntegrationInfo<dim> &) const
- {}
-
-
- template <int dim>
- class Estimator : public MeshWorker::LocalIntegrator<dim>
- {
- public:
- void cell(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const;
- void boundary(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const;
- void face(MeshWorker::DoFInfo<dim> &dinfo1,
- MeshWorker::DoFInfo<dim> &dinfo2,
- typename MeshWorker::IntegrationInfo<dim> &info1,
- typename MeshWorker::IntegrationInfo<dim> &info2) const;
- };
-
-
- template <int dim>
- void Estimator<dim>::cell(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const
- {
- const FEValuesBase<dim> &fe = info.fe_values();
-
- const std::vector<Tensor<2,dim> > &DDuh = info.hessians[0][0];
- for (unsigned k=0; k<fe.n_quadrature_points; ++k)
- {
- const double t = dinfo.cell->diameter() * trace(DDuh[k]);
- dinfo.value(0) += t*t * fe.JxW(k);
- }
- dinfo.value(0) = std::sqrt(dinfo.value(0));
- }
-
- template <int dim>
- void Estimator<dim>::boundary(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const
- {
- const FEValuesBase<dim> &fe = info.fe_values();
-
- std::vector<double> boundary_values(fe.n_quadrature_points);
- exact_solution.value_list(fe.get_quadrature_points(), boundary_values);
-
- const std::vector<double> &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<fe.n_quadrature_points; ++k)
- dinfo.value(0) += penalty * (boundary_values[k] - uh[k]) * (boundary_values[k] - uh[k])
- * fe.JxW(k);
- dinfo.value(0) = std::sqrt(dinfo.value(0));
- }
-
-
- template <int dim>
- void Estimator<dim>::face(MeshWorker::DoFInfo<dim> &dinfo1,
- MeshWorker::DoFInfo<dim> &dinfo2,
- typename MeshWorker::IntegrationInfo<dim> &info1,
- typename MeshWorker::IntegrationInfo<dim> &info2) const
- {
- const FEValuesBase<dim> &fe = info1.fe_values();
- const std::vector<double> &uh1 = info1.values[0][0];
- const std::vector<double> &uh2 = info2.values[0][0];
- const std::vector<Tensor<1,dim> > &Duh1 = info1.gradients[0][0];
- const std::vector<Tensor<1,dim> > &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<fe.n_quadrature_points; ++k)
- {
- double diff1 = uh1[k] - uh2[k];
- double diff2 = fe.normal_vector(k) * Duh1[k] - fe.normal_vector(k) * Duh2[k];
- dinfo1.value(0) += (penalty * diff1*diff1 + h * diff2*diff2)
- * fe.JxW(k);
- }
- dinfo1.value(0) = std::sqrt(dinfo1.value(0));
- dinfo2.value(0) = dinfo1.value(0);
- }
-
-
-
- template <int dim>
- class ErrorIntegrator : public MeshWorker::LocalIntegrator<dim>
- {
- public:
- void cell(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const;
- void boundary(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const;
- void face(MeshWorker::DoFInfo<dim> &dinfo1,
- MeshWorker::DoFInfo<dim> &dinfo2,
- typename MeshWorker::IntegrationInfo<dim> &info1,
- typename MeshWorker::IntegrationInfo<dim> &info2) const;
- };
-
-
- template <int dim>
- void ErrorIntegrator<dim>::cell(
- MeshWorker::DoFInfo<dim> &dinfo,
- typename MeshWorker::IntegrationInfo<dim> &info) const
- {
- const FEValuesBase<dim> &fe = info.fe_values();
- std::vector<Tensor<1,dim> > exact_gradients(fe.n_quadrature_points);
- std::vector<double> 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<Tensor<1,dim> > &Duh = info.gradients[0][0];
- const std::vector<double> &uh = info.values[0][0];
-
- for (unsigned k=0; k<fe.n_quadrature_points; ++k)
- {
- double sum = 0;
- for (unsigned int d=0; d<dim; ++d)
- {
- const double diff = exact_gradients[k][d] - Duh[k][d];
- sum += diff*diff;
- }
- const double diff = exact_values[k] - uh[k];
- dinfo.value(0) += sum * fe.JxW(k);
- dinfo.value(1) += diff*diff * fe.JxW(k);
- }
- dinfo.value(0) = std::sqrt(dinfo.value(0));
- dinfo.value(1) = std::sqrt(dinfo.value(1));
- }
-
-
- template <int dim>
- void ErrorIntegrator<dim>::boundary(
- MeshWorker::DoFInfo<dim> &dinfo,
- typename MeshWorker::IntegrationInfo<dim> &info) const
- {
- const FEValuesBase<dim> &fe = info.fe_values();
-
- std::vector<double> exact_values(fe.n_quadrature_points);
- exact_solution.value_list(fe.get_quadrature_points(), exact_values);
-
- const std::vector<double> &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<fe.n_quadrature_points; ++k)
- {
- const double diff = exact_values[k] - uh[k];
- dinfo.value(0) += penalty * diff * diff * fe.JxW(k);
- }
- dinfo.value(0) = std::sqrt(dinfo.value(0));
- }
-
-
- template <int dim>
- void ErrorIntegrator<dim>::face(
- MeshWorker::DoFInfo<dim> &dinfo1,
- MeshWorker::DoFInfo<dim> &dinfo2,
- typename MeshWorker::IntegrationInfo<dim> &info1,
- typename MeshWorker::IntegrationInfo<dim> &info2) const
- {
- const FEValuesBase<dim> &fe = info1.fe_values();
- const std::vector<double> &uh1 = info1.values[0][0];
- const std::vector<double> &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<fe.n_quadrature_points; ++k)
- {
- double diff = uh1[k] - uh2[k];
- dinfo1.value(0) += (penalty * diff*diff)
- * fe.JxW(k);
- }
- dinfo1.value(0) = std::sqrt(dinfo1.value(0));
- dinfo2.value(0) = dinfo1.value(0);
- }
-
-
-
-
template <int dim>
class InteriorPenaltyProblem
{
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;
parallel::distributed::Triangulation<dim> triangulation;
const MappingQ1<dim> mapping;
IndexSet locally_relevant_set;
TrilinosWrappers::SparseMatrix matrix;
- TrilinosWrappers::MPI::Vector solution;
- TrilinosWrappers::MPI::Vector right_hand_side;
- BlockVector<double> estimates;
MGLevelObject<TrilinosWrappers::SparseMatrix> mg_matrix;
dof_handler.distribute_mg_dofs (fe);
DoFTools::extract_locally_relevant_dofs (dof_handler, locally_relevant_set);
- solution.reinit(dof_handler.locally_owned_dofs(), MPI_COMM_WORLD);
- right_hand_side.reinit(dof_handler.locally_owned_dofs(), MPI_COMM_WORLD);
-
+
CompressedSimpleSparsityPattern c_sparsity(dof_handler.n_dofs(), dof_handler.n_dofs());
DoFTools::make_flux_sparsity_pattern(dof_handler, c_sparsity);
matrix.reinit(dof_handler.locally_owned_dofs(), c_sparsity, MPI_COMM_WORLD, true);
}
- template <int dim>
- void
- InteriorPenaltyProblem<dim>::assemble_matrix()
- {
- MeshWorker::IntegrationInfoBox<dim> info_box;
- UpdateFlags update_flags = update_values | update_gradients;
- info_box.add_update_flags_all(update_flags);
- info_box.initialize(fe, mapping);
-
- MeshWorker::DoFInfo<dim> dof_info(dof_handler);
-
- MeshWorker::Assembler::MatrixSimple<TrilinosWrappers::SparseMatrix> assembler;
- assembler.initialize(matrix);
-
- FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
- begin(IteratorFilters::LocallyOwnedCell(), dof_handler.begin_active());
- FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
- end(IteratorFilters::LocallyOwnedCell(), dof_handler.end());
-
- MatrixIntegrator<dim> integrator;
- MeshWorker::integration_loop<dim, dim>(
- begin, end,
- dof_info, info_box,
- integrator, assembler);
-
- matrix.compress(VectorOperation::add);
- }
-
-
- template <int dim>
- void
- InteriorPenaltyProblem<dim>::assemble_mg_matrix()
- {
- MeshWorker::IntegrationInfoBox<dim> info_box;
- UpdateFlags update_flags = update_values | update_gradients;
- info_box.add_update_flags_all(update_flags);
- info_box.initialize(fe, mapping);
-
- MeshWorker::DoFInfo<dim> dof_info(dof_handler);
-
- MeshWorker::Assembler::MGMatrixSimple<TrilinosWrappers::SparseMatrix> assembler;
- assembler.initialize(mg_matrix);
- assembler.initialize_fluxes(mg_matrix_dg_up, mg_matrix_dg_down);
-
- FilteredIterator<typename DoFHandler<dim>::level_cell_iterator>
- begin(IteratorFilters::LocallyOwnedLevelCell(), dof_handler.begin());
- FilteredIterator<typename DoFHandler<dim>::level_cell_iterator>
- end(IteratorFilters::LocallyOwnedLevelCell(), dof_handler.end());
-
- MatrixIntegrator<dim> integrator;
- MeshWorker::integration_loop<dim, dim> (
- begin, end,
- dof_info, info_box,
- integrator, assembler);
-
- for (unsigned int level=mg_matrix.min_level(); level <= mg_matrix.max_level();++level)
- {
- mg_matrix[level].compress(VectorOperation::add);
- if (level > mg_matrix.min_level())
- {
- mg_matrix_dg_up[level].compress(VectorOperation::add);
- mg_matrix_dg_down[level].compress(VectorOperation::add);
- }
- }
- }
-
-
- template <int dim>
- void
- InteriorPenaltyProblem<dim>::assemble_right_hand_side()
- {
- MeshWorker::IntegrationInfoBox<dim> 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<dim> dof_info(dof_handler);
-
- MeshWorker::Assembler::ResidualSimple<TrilinosWrappers::MPI::Vector> assembler;
- NamedData<TrilinosWrappers::MPI::Vector*> data;
- TrilinosWrappers::MPI::Vector *rhs = &right_hand_side;
- data.add(rhs, "RHS");
- assembler.initialize(data);
-
- FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
- begin(IteratorFilters::LocallyOwnedCell(), dof_handler.begin_active());
- FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
- end(IteratorFilters::LocallyOwnedCell(), dof_handler.end());
-
-
- RHSIntegrator<dim> integrator;
- MeshWorker::integration_loop<dim, dim>(begin, end,
- dof_info, info_box,
- integrator, assembler);
-
- right_hand_side.compress(VectorOperation::add);
- right_hand_side *= -1.;
- }
-
-
- template <int dim>
- void
- InteriorPenaltyProblem<dim>::solve()
- {
- SolverControl control(1000, 1.e-12);
- SolverCG<TrilinosWrappers::MPI::Vector > solver(control);
-
- MGTransferPrebuilt<TrilinosWrappers::MPI::Vector > mg_transfer;
- mg_transfer.build_matrices(dof_handler);
-
- SolverControl coarse_solver_control (1000, 1e-10, false, false);
- SolverCG<TrilinosWrappers::MPI::Vector> coarse_solver(coarse_solver_control);
- PreconditionIdentity identity;
- TrilinosWrappers::SparseMatrix coarse_matrix = mg_matrix[0];
- MGCoarseGridLACIteration<SolverCG<TrilinosWrappers::MPI::Vector>,TrilinosWrappers::MPI::Vector>
- coarse_grid_solver(coarse_solver, coarse_matrix, identity);
-
- typedef TrilinosWrappers::PreconditionJacobi Smoother;
- MGSmootherPrecondition<TrilinosWrappers::SparseMatrix, Smoother, TrilinosWrappers::MPI::Vector> mg_smoother;
- mg_smoother.initialize(mg_matrix, .7);
- mg_smoother.set_steps(2);
-
- MGMatrix<TrilinosWrappers::SparseMatrix, TrilinosWrappers::MPI::Vector> mgmatrix(&mg_matrix);
- MGMatrix<TrilinosWrappers::SparseMatrix, TrilinosWrappers::MPI::Vector> mgdown(&mg_matrix_dg_down);
- MGMatrix<TrilinosWrappers::SparseMatrix, TrilinosWrappers::MPI::Vector> mgup(&mg_matrix_dg_up);
-
- Multigrid<TrilinosWrappers::MPI::Vector > mg(dof_handler, mgmatrix,
- coarse_grid_solver, mg_transfer,
- mg_smoother, mg_smoother);
- mg.set_edge_flux_matrices(mgdown, mgup);
-
- PreconditionMG<dim, TrilinosWrappers::MPI::Vector,
- MGTransferPrebuilt<TrilinosWrappers::MPI::Vector > >
- preconditioner(dof_handler, mg, mg_transfer);
- solver.solve(matrix, solution, right_hand_side, preconditioner);
- }
-
-
- template <int dim>
- double
- InteriorPenaltyProblem<dim>::estimate()
- {
- TrilinosWrappers::MPI::Vector ghost;
- ghost.reinit(locally_relevant_set, MPI_COMM_WORLD);
- ghost = solution;
-
- std::vector<unsigned int> 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<dim>::active_cell_iterator cell = triangulation.begin_active();
- cell != triangulation.end(); ++cell,++i)
- cell->set_user_index(i);
-
- MeshWorker::IntegrationInfoBox<dim> 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<TrilinosWrappers::MPI::Vector* > solution_data;
- solution_data.add(&ghost, "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<dim> dof_info(dof_handler);
-
- MeshWorker::Assembler::CellsAndFaces<double> assembler;
- NamedData<BlockVector<double>* > out_data;
- BlockVector<double> *est = &estimates;
- out_data.add(est, "cells");
- assembler.initialize(out_data, false);
-
- FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
- begin(IteratorFilters::LocallyOwnedCell(), dof_handler.begin_active());
- FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
- end(IteratorFilters::LocallyOwnedCell(), dof_handler.end());
-
- Estimator<dim> integrator;
- MeshWorker::integration_loop<dim, dim> (
- begin, end,
- dof_info, info_box,
- integrator, assembler);
-
- triangulation.load_user_indices(old_user_indices);
- return estimates.block(0).l2_norm();
- }
-
-
- // template <int dim>
- // void
- // InteriorPenaltyProblem<dim>::error()
- // {
- // BlockVector<double> 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<dim>::active_cell_iterator cell = triangulation.begin_active();
- // cell != triangulation.end(); ++cell,++i)
- // cell->set_user_index(i);
-
- // MeshWorker::IntegrationInfoBox<dim> 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<TrilinosWrappers::MPI::Vector* > 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<dim> dof_info(dof_handler);
-
- // MeshWorker::Assembler::CellsAndFaces<double> assembler;
- // NamedData<BlockVector<double>* > out_data;
- // BlockVector<double> *est = &errors;
- // out_data.add(est, "cells");
- // assembler.initialize(out_data, false);
-
- // ErrorIntegrator<dim> integrator;
- // MeshWorker::integration_loop<dim, dim> (
- // 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 <int dim>
- // void InteriorPenaltyProblem<dim>::output_results (const unsigned int cycle) const
- // {
- // char *fn = new char[100];
- // sprintf(fn, "step-39/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<dim> 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 <int dim>
void
InteriorPenaltyProblem<dim>::run(unsigned int n_steps)
for (unsigned int s=0; s<n_steps; ++s)
{
deallog << "Step " << s << std::endl;
- if (estimates.block(0).size() == 0)
+ if (s==0)
triangulation.refine_global(1);
else
{
if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
triangulation.begin_active(triangulation.n_levels()-1)->set_refine_flag();
- // parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction (triangulation,
- // estimates.block(0),
- // 0.5, 0.0);
triangulation.execute_coarsening_and_refinement ();
}
for (unsigned int l=0; l<triangulation.n_levels(); ++l)
deallog << ' ' << dof_handler.n_dofs(l);
deallog << std::endl;
-
- deallog << "Assemble matrix" << std::endl;
- assemble_matrix();
- deallog << "Assemble multilevel matrix" << std::endl;
- assemble_mg_matrix();
- deallog << "Assemble right hand side" << std::endl;
- assemble_right_hand_side();
- deallog << "Solve" << std::endl;
- solve();
- //error();
- deallog << "Estimate " << estimate() << std::endl;
- //output_results(s);
}
}
}
try
{
- FE_DGQ<2> fe1(2);
+ FE_DGQ<2> fe1(0);
InteriorPenaltyProblem<2> test1(fe1);
test1.run(20);
}
projects / dealii-svn.git / commitdiff