--- /dev/null
+//-----------------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2012 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 for correctness of step-37 (without output and only small sizes)
+
+
+#include "../tests.h"
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/multigrid/multigrid.h>
+#include <deal.II/multigrid/mg_dof_handler.h>
+#include <deal.II/multigrid/mg_dof_accessor.h>
+#include <deal.II/multigrid/mg_transfer.h>
+#include <deal.II/multigrid/mg_tools.h>
+#include <deal.II/multigrid/mg_coarse.h>
+#include <deal.II/multigrid/mg_smoother.h>
+#include <deal.II/multigrid/mg_matrix.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+#include <fstream>
+#include <sstream>
+
+
+namespace Step37
+{
+ using namespace dealii;
+
+ const unsigned int degree_finite_element = 2;
+
+
+
+ template <int dim>
+ class Coefficient : public Function<dim>
+ {
+ public:
+ Coefficient () : Function<dim>() {}
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ template <typename number>
+ number value (const Point<dim,number> &p,
+ const unsigned int component = 0) const;
+
+ virtual void value_list (const std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int component = 0) const;
+ };
+
+
+
+ template <int dim>
+ template <typename number>
+ number Coefficient<dim>::value (const Point<dim,number> &p,
+ const unsigned int /*component*/) const
+ {
+ return 1. / (0.05 + 2.*p.square());
+ }
+
+
+
+ template <int dim>
+ double Coefficient<dim>::value (const Point<dim> &p,
+ const unsigned int component) const
+ {
+ return value<double>(p,component);
+ }
+
+
+
+ template <int dim>
+ void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int component) const
+ {
+ Assert (values.size() == points.size(),
+ ExcDimensionMismatch (values.size(), points.size()));
+ Assert (component == 0,
+ ExcIndexRange (component, 0, 1));
+
+ const unsigned int n_points = points.size();
+ for (unsigned int i=0; i<n_points; ++i)
+ values[i] = value<double>(points[i],component);
+ }
+
+
+
+
+
+ template <int dim, int fe_degree, typename number>
+ class LaplaceOperator : public Subscriptor
+ {
+ public:
+ LaplaceOperator ();
+
+ void clear();
+
+ void reinit (const MGDoFHandler<dim> &dof_handler,
+ const ConstraintMatrix &constraints,
+ const unsigned int level = numbers::invalid_unsigned_int);
+
+ unsigned int m () const;
+ unsigned int n () const;
+
+ void vmult (Vector<double> &dst,
+ const Vector<double> &src) const;
+ void Tvmult (Vector<double> &dst,
+ const Vector<double> &src) const;
+ void vmult_add (Vector<double> &dst,
+ const Vector<double> &src) const;
+ void Tvmult_add (Vector<double> &dst,
+ const Vector<double> &src) const;
+
+ number el (const unsigned int row,
+ const unsigned int col) const;
+ void set_diagonal (const Vector<number> &diagonal);
+
+ std::size_t memory_consumption () const;
+
+ private:
+ void local_apply (const MatrixFree<dim,number> &data,
+ Vector<double> &dst,
+ const Vector<double> &src,
+ const std::pair<unsigned int,unsigned int> &cell_range) const;
+
+ void evaluate_coefficient(const Coefficient<dim> &function);
+
+ MatrixFree<dim,number> data;
+ AlignedVector<VectorizedArray<number> > coefficient;
+
+ Vector<number> diagonal_values;
+ bool diagonal_is_available;
+ };
+
+
+
+ template <int dim, int fe_degree, typename number>
+ LaplaceOperator<dim,fe_degree,number>::LaplaceOperator ()
+ :
+ Subscriptor()
+ {}
+
+
+
+ template <int dim, int fe_degree, typename number>
+ unsigned int
+ LaplaceOperator<dim,fe_degree,number>::m () const
+ {
+ return data.get_vector_partitioner()->size();
+ }
+
+
+
+ template <int dim, int fe_degree, typename number>
+ unsigned int
+ LaplaceOperator<dim,fe_degree,number>::n () const
+ {
+ return data.get_vector_partitioner()->size();
+ }
+
+
+
+ template <int dim, int fe_degree, typename number>
+ void
+ LaplaceOperator<dim,fe_degree,number>::clear ()
+ {
+ data.clear();
+ diagonal_is_available = false;
+ diagonal_values.reinit(0);
+ }
+
+
+
+ template <int dim, int fe_degree, typename number>
+ void
+ LaplaceOperator<dim,fe_degree,number>::reinit (const MGDoFHandler<dim> &dof_handler,
+ const ConstraintMatrix &constraints,
+ const unsigned int level)
+ {
+ typename MatrixFree<dim,number>::AdditionalData additional_data;
+ additional_data.tasks_parallel_scheme =
+ MatrixFree<dim,number>::AdditionalData::partition_color;
+ additional_data.level_mg_handler = level;
+ additional_data.mapping_update_flags = (update_gradients | update_JxW_values |
+ update_quadrature_points);
+ data.reinit (dof_handler, constraints, QGauss<1>(fe_degree+1),
+ additional_data);
+ evaluate_coefficient(Coefficient<dim>());
+ }
+
+
+
+ template <int dim, int fe_degree, typename number>
+ void
+ LaplaceOperator<dim,fe_degree,number>::
+ evaluate_coefficient (const Coefficient<dim> &coefficient_function)
+ {
+ const unsigned int n_cells = data.get_size_info().n_macro_cells;
+ FEEvaluation<dim,fe_degree,fe_degree+1,1,number> phi (data);
+ coefficient.resize (n_cells * phi.n_q_points);
+ for (unsigned int cell=0; cell<n_cells; ++cell)
+ {
+ phi.reinit (cell);
+ for (unsigned int q=0; q<phi.n_q_points; ++q)
+ coefficient[cell*phi.n_q_points+q] =
+ coefficient_function.value(phi.quadrature_point(q));
+ }
+ }
+
+
+
+
+ template <int dim, int fe_degree, typename number>
+ void
+ LaplaceOperator<dim,fe_degree,number>::
+ local_apply (const MatrixFree<dim,number> &data,
+ Vector<double> &dst,
+ const Vector<double> &src,
+ const std::pair<unsigned int,unsigned int> &cell_range) const
+ {
+ FEEvaluation<dim,fe_degree,fe_degree+1,1,number> phi (data);
+ AssertDimension (coefficient.size(),
+ data.get_size_info().n_macro_cells * phi.n_q_points);
+
+ for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
+ {
+ phi.reinit (cell);
+ phi.read_dof_values(src);
+ phi.evaluate (false,true,false);
+ for (unsigned int q=0; q<phi.n_q_points; ++q)
+ phi.submit_gradient (coefficient[cell*phi.n_q_points+q] *
+ phi.get_gradient(q), q);
+ phi.integrate (false,true);
+ phi.distribute_local_to_global (dst);
+ }
+ }
+
+
+
+
+ template <int dim, int fe_degree, typename number>
+ void
+ LaplaceOperator<dim,fe_degree,number>::vmult (Vector<double> &dst,
+ const Vector<double> &src) const
+ {
+ dst = 0;
+ vmult_add (dst, src);
+ }
+
+
+
+ template <int dim, int fe_degree, typename number>
+ void
+ LaplaceOperator<dim,fe_degree,number>::Tvmult (Vector<double> &dst,
+ const Vector<double> &src) const
+ {
+ dst = 0;
+ vmult_add (dst,src);
+ }
+
+
+
+ template <int dim, int fe_degree, typename number>
+ void
+ LaplaceOperator<dim,fe_degree,number>::Tvmult_add (Vector<double> &dst,
+ const Vector<double> &src) const
+ {
+ vmult_add (dst,src);
+ }
+
+
+
+ template <int dim, int fe_degree, typename number>
+ void
+ LaplaceOperator<dim,fe_degree,number>::vmult_add (Vector<double> &dst,
+ const Vector<double> &src) const
+ {
+ data.cell_loop (&LaplaceOperator::local_apply, this, dst, src);
+
+ const std::vector<unsigned int> &
+ constrained_dofs = data.get_constrained_dofs();
+ for (unsigned int i=0; i<constrained_dofs.size(); ++i)
+ dst(constrained_dofs[i]) += src(constrained_dofs[i]);
+ }
+
+
+
+ template <int dim, int fe_degree, typename number>
+ number
+ LaplaceOperator<dim,fe_degree,number>::el (const unsigned int row,
+ const unsigned int col) const
+ {
+ Assert (row == col, ExcNotImplemented());
+ Assert (diagonal_is_available == true, ExcNotInitialized());
+ return diagonal_values(row);
+ }
+
+
+
+ template <int dim, int fe_degree, typename number>
+ void
+ LaplaceOperator<dim,fe_degree,number>::set_diagonal(const Vector<number> &diagonal)
+ {
+ AssertDimension (m(), diagonal.size());
+
+ diagonal_values = diagonal;
+
+ const std::vector<unsigned int> &
+ constrained_dofs = data.get_constrained_dofs();
+ for (unsigned int i=0; i<constrained_dofs.size(); ++i)
+ diagonal_values(constrained_dofs[i]) = 1.0;
+
+ diagonal_is_available = true;
+ }
+
+
+
+ template <int dim, int fe_degree, typename number>
+ std::size_t
+ LaplaceOperator<dim,fe_degree,number>::memory_consumption () const
+ {
+ return (data.memory_consumption () +
+ coefficient.memory_consumption() +
+ diagonal_values.memory_consumption() +
+ MemoryConsumption::memory_consumption(diagonal_is_available));
+ }
+
+
+
+
+ template <int dim>
+ class LaplaceProblem
+ {
+ public:
+ LaplaceProblem ();
+ void run ();
+
+ private:
+ void setup_system ();
+ void assemble_system ();
+ void assemble_multigrid ();
+ void solve ();
+ void output_results (const unsigned int cycle) const;
+
+ typedef LaplaceOperator<dim,degree_finite_element,double> SystemMatrixType;
+ typedef LaplaceOperator<dim,degree_finite_element,float> LevelMatrixType;
+
+ Triangulation<dim> triangulation;
+ FE_Q<dim> fe;
+ MGDoFHandler<dim> mg_dof_handler;
+ ConstraintMatrix constraints;
+
+ SystemMatrixType system_matrix;
+ MGLevelObject<LevelMatrixType> mg_matrices;
+ FullMatrix<float> coarse_matrix;
+ MGLevelObject<ConstraintMatrix> mg_constraints;
+
+ Vector<double> solution;
+ Vector<double> system_rhs;
+ };
+
+
+
+ template <int dim>
+ LaplaceProblem<dim>::LaplaceProblem ()
+ :
+ fe (degree_finite_element),
+ mg_dof_handler (triangulation)
+ {}
+
+
+
+
+ template <int dim>
+ void LaplaceProblem<dim>::setup_system ()
+ {
+ system_matrix.clear();
+ mg_matrices.clear();
+ mg_constraints.clear();
+
+ mg_dof_handler.distribute_dofs (fe);
+
+ deallog << "Number of degrees of freedom: "
+ << mg_dof_handler.n_dofs()
+ << std::endl;
+
+ constraints.clear();
+ VectorTools::interpolate_boundary_values (mg_dof_handler,
+ 0,
+ ZeroFunction<dim>(),
+ constraints);
+ constraints.close();
+
+ system_matrix.reinit (mg_dof_handler, constraints);
+ solution.reinit (mg_dof_handler.n_dofs());
+ system_rhs.reinit (mg_dof_handler.n_dofs());
+
+ const unsigned int nlevels = triangulation.n_levels();
+ mg_matrices.resize(0, nlevels-1);
+ mg_constraints.resize (0, nlevels-1);
+
+ typename FunctionMap<dim>::type dirichlet_boundary;
+ ZeroFunction<dim> homogeneous_dirichlet_bc (1);
+ dirichlet_boundary[0] = &homogeneous_dirichlet_bc;
+ std::vector<std::set<unsigned int> > boundary_indices(triangulation.n_levels());
+ MGTools::make_boundary_list (mg_dof_handler,
+ dirichlet_boundary,
+ boundary_indices);
+ for (unsigned int level=0; level<nlevels; ++level)
+ {
+ std::set<unsigned int>::iterator bc_it = boundary_indices[level].begin();
+ for ( ; bc_it != boundary_indices[level].end(); ++bc_it)
+ mg_constraints[level].add_line(*bc_it);
+
+ mg_constraints[level].close();
+ mg_matrices[level].reinit(mg_dof_handler,
+ mg_constraints[level],
+ level);
+ }
+ coarse_matrix.reinit (mg_dof_handler.n_dofs(0),
+ mg_dof_handler.n_dofs(0));
+ }
+
+
+
+
+ template <int dim>
+ void LaplaceProblem<dim>::assemble_system ()
+ {
+ QGauss<dim> quadrature_formula(fe.degree+1);
+ FEValues<dim> fe_values (fe, quadrature_formula,
+ update_values | update_JxW_values);
+
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ const unsigned int n_q_points = quadrature_formula.size();
+
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+ const Coefficient<dim> coefficient;
+ std::vector<double> coefficient_values (n_q_points);
+
+ typename DoFHandler<dim>::active_cell_iterator cell = mg_dof_handler.begin_active(),
+ endc = mg_dof_handler.end();
+ for (; cell!=endc; ++cell)
+ {
+ cell->get_dof_indices (local_dof_indices);
+ fe_values.reinit (cell);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ double rhs_val = 0;
+ for (unsigned int q=0; q<n_q_points; ++q)
+ rhs_val += (fe_values.shape_value(i,q) * 1.0 *
+ fe_values.JxW(q));
+ system_rhs(local_dof_indices[i]) += rhs_val;
+ }
+ }
+ constraints.condense(system_rhs);
+ }
+
+
+
+ template <int dim>
+ void LaplaceProblem<dim>::assemble_multigrid ()
+ {
+ coarse_matrix = 0;
+ QGauss<dim> quadrature_formula(fe.degree+1);
+ FEValues<dim> fe_values (fe, quadrature_formula,
+ update_gradients | update_inverse_jacobians |
+ update_quadrature_points | update_JxW_values);
+
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ const unsigned int n_q_points = quadrature_formula.size();
+
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+ const Coefficient<dim> coefficient;
+ std::vector<double> coefficient_values (n_q_points);
+ FullMatrix<double> local_matrix (dofs_per_cell, dofs_per_cell);
+ Vector<double> local_diagonal (dofs_per_cell);
+
+ const unsigned int n_levels = triangulation.n_levels();
+ std::vector<Vector<float> > diagonals (n_levels);
+ for (unsigned int level=0; level<n_levels; ++level)
+ diagonals[level].reinit (mg_dof_handler.n_dofs(level));
+
+ std::vector<unsigned int> cell_no(triangulation.n_levels());
+ typename MGDoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
+ endc = mg_dof_handler.end();
+ for (; cell!=endc; ++cell)
+ {
+ const unsigned int level = cell->level();
+ cell->get_mg_dof_indices (local_dof_indices);
+ fe_values.reinit (cell);
+ coefficient.value_list (fe_values.get_quadrature_points(),
+ coefficient_values);
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ double local_diag = 0;
+ for (unsigned int q=0; q<n_q_points; ++q)
+ local_diag += ((fe_values.shape_grad(i,q) *
+ fe_values.shape_grad(i,q)) *
+ coefficient_values[q] * fe_values.JxW(q));
+ local_diagonal(i) = local_diag;
+ }
+ mg_constraints[level].distribute_local_to_global(local_diagonal,
+ local_dof_indices,
+ diagonals[level]);
+
+ if (level == 0)
+ {
+ local_matrix = 0;
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ double add_value = 0;
+ for (unsigned int q=0; q<n_q_points; ++q)
+ add_value += (fe_values.shape_grad(i,q) *
+ fe_values.shape_grad(j,q) *
+ coefficient_values[q] *
+ fe_values.JxW(q));
+ local_matrix(i,j) = add_value;
+ }
+ mg_constraints[0].distribute_local_to_global (local_matrix,
+ local_dof_indices,
+ coarse_matrix);
+ }
+ }
+
+ for (unsigned int level=0; level<n_levels; ++level)
+ mg_matrices[level].set_diagonal (diagonals[level]);
+ }
+
+
+
+ template <int dim>
+ void LaplaceProblem<dim>::solve ()
+ {
+ GrowingVectorMemory<> vector_memory;
+
+ MGTransferPrebuilt<Vector<double> > mg_transfer;
+ mg_transfer.build_matrices(mg_dof_handler);
+
+ MGCoarseGridHouseholder<float, Vector<double> > mg_coarse;
+ mg_coarse.initialize(coarse_matrix);
+
+ typedef PreconditionChebyshev<LevelMatrixType,Vector<double> > SMOOTHER;
+ MGSmootherPrecondition<LevelMatrixType, SMOOTHER, Vector<double> >
+ mg_smoother(vector_memory);
+
+ typename SMOOTHER::AdditionalData smoother_data;
+ smoother_data.smoothing_range = 10.;
+ smoother_data.degree = 6;
+ smoother_data.eig_cg_n_iterations = 10;
+ mg_smoother.initialize(mg_matrices, smoother_data);
+
+ MGMatrix<LevelMatrixType, Vector<double> >
+ mg_matrix(&mg_matrices);
+
+ Multigrid<Vector<double> > mg(mg_dof_handler,
+ mg_matrix,
+ mg_coarse,
+ mg_transfer,
+ mg_smoother,
+ mg_smoother);
+ PreconditionMG<dim, Vector<double>,
+ MGTransferPrebuilt<Vector<double> > >
+ preconditioner(mg_dof_handler, mg, mg_transfer);
+
+ const std::size_t multigrid_memory
+ = (mg_matrices.memory_consumption() +
+ mg_transfer.memory_consumption() +
+ coarse_matrix.memory_consumption());
+
+ SolverControl solver_control (1000, 1e-12*system_rhs.l2_norm());
+ SolverCG<> cg (solver_control);
+
+ cg.solve (system_matrix, solution, system_rhs,
+ preconditioner);
+ }
+
+
+
+ template <int dim>
+ void LaplaceProblem<dim>::run ()
+ {
+ for (unsigned int cycle=0; cycle<3; ++cycle)
+ {
+ deallog << "Cycle " << cycle << std::endl;
+
+ if (cycle == 0)
+ {
+ GridGenerator::hyper_cube (triangulation, 0., 1.);
+ triangulation.refine_global (3-dim);
+ }
+ triangulation.refine_global (1);
+ setup_system ();
+ assemble_system ();
+ assemble_multigrid ();
+ solve ();
+ deallog << std::endl;
+ };
+ }
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile("step-37/output");
+ deallog.attach(logfile);
+ deallog << std::setprecision (3);
+ deallog.threshold_double(1e-10);
+ deallog.depth_console(0);
+
+ {
+ deallog.push("2d");
+ Step37::LaplaceProblem<2> laplace_problem;
+ laplace_problem.run();
+ deallog.pop();
+ }
+ {
+ deallog.push("3d");
+ Step37::LaplaceProblem<3> laplace_problem;
+ laplace_problem.run();
+ deallog.pop();
+ }
+
+
+ return 0;
+}
+
--- /dev/null
+//-----------------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2012 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 for correctness of step-48 (without output and only small sizes)
+
+
+#include "../tests.h"
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/function.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/lac/parallel_vector.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+
+
+namespace Step48
+{
+ using namespace dealii;
+
+ const unsigned int dimension = 2;
+ const unsigned int fe_degree = 4;
+
+
+
+ template <int dim, int fe_degree>
+ class SineGordonOperation
+ {
+ public:
+ SineGordonOperation(const MatrixFree<dim,double> &data_in,
+ const double time_step);
+
+ void apply (parallel::distributed::Vector<double> &dst,
+ const std::vector<parallel::distributed::Vector<double>*> &src) const;
+
+ private:
+ const MatrixFree<dim,double> &data;
+ const VectorizedArray<double> delta_t_sqr;
+ parallel::distributed::Vector<double> inv_mass_matrix;
+
+ void local_apply (const MatrixFree<dim,double> &data,
+ parallel::distributed::Vector<double> &dst,
+ const std::vector<parallel::distributed::Vector<double>*> &src,
+ const std::pair<unsigned int,unsigned int> &cell_range) const;
+ };
+
+
+
+
+ template <int dim, int fe_degree>
+ SineGordonOperation<dim,fe_degree>::
+ SineGordonOperation(const MatrixFree<dim,double> &data_in,
+ const double time_step)
+ :
+ data(data_in),
+ delta_t_sqr(make_vectorized_array(time_step *time_step))
+ {
+ VectorizedArray<double> one = make_vectorized_array (1.);
+
+ data.initialize_dof_vector (inv_mass_matrix);
+
+ FEEvaluationGL<dim,fe_degree> fe_eval(data);
+ const unsigned int n_q_points = fe_eval.n_q_points;
+
+ for (unsigned int cell=0; cell<data.get_size_info().n_macro_cells; ++cell)
+ {
+ fe_eval.reinit(cell);
+ for (unsigned int q=0; q<n_q_points; ++q)
+ fe_eval.submit_value(one,q);
+ fe_eval.integrate (true,false);
+ fe_eval.distribute_local_to_global (inv_mass_matrix);
+ }
+
+ inv_mass_matrix.compress();
+ for (unsigned int k=0; k<inv_mass_matrix.local_size(); ++k)
+ if (inv_mass_matrix.local_element(k)>1e-15)
+ inv_mass_matrix.local_element(k) = 1./inv_mass_matrix.local_element(k);
+ else
+ inv_mass_matrix.local_element(k) = 0;
+ }
+
+
+
+
+
+ template <int dim, int fe_degree>
+ void SineGordonOperation<dim, fe_degree>::
+ local_apply (const MatrixFree<dim> &data,
+ parallel::distributed::Vector<double> &dst,
+ const std::vector<parallel::distributed::Vector<double>*> &src,
+ const std::pair<unsigned int,unsigned int> &cell_range) const
+ {
+ AssertDimension (src.size(), 2);
+ FEEvaluationGL<dim,fe_degree> current (data), old (data);
+ for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
+ {
+ current.reinit (cell);
+ old.reinit (cell);
+
+ current.read_dof_values (*src[0]);
+ old.read_dof_values (*src[1]);
+
+ current.evaluate (true, true, false);
+ old.evaluate (true, false, false);
+
+ for (unsigned int q=0; q<current.n_q_points; ++q)
+ {
+ const VectorizedArray<double> current_value = current.get_value(q);
+ const VectorizedArray<double> old_value = old.get_value(q);
+
+ current.submit_value (2.*current_value - old_value -
+ delta_t_sqr * std::sin(current_value),q);
+ current.submit_gradient (- delta_t_sqr *
+ current.get_gradient(q), q);
+ }
+
+ current.integrate (true,true);
+ current.distribute_local_to_global (dst);
+ }
+ }
+
+
+
+
+ template <int dim, int fe_degree>
+ void SineGordonOperation<dim, fe_degree>::
+ apply (parallel::distributed::Vector<double> &dst,
+ const std::vector<parallel::distributed::Vector<double>*> &src) const
+ {
+ dst = 0;
+ data.cell_loop (&SineGordonOperation<dim,fe_degree>::local_apply,
+ this, dst, src);
+ dst.scale(inv_mass_matrix);
+ }
+
+
+
+ template <int dim>
+ class ExactSolution : public Function<dim>
+ {
+ public:
+ ExactSolution (const unsigned int n_components = 1,
+ const double time = 0.) : Function<dim>(n_components, time) {}
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+ };
+
+ template <int dim>
+ double ExactSolution<dim>::value (const Point<dim> &p,
+ const unsigned int /* component */) const
+ {
+ double t = this->get_time ();
+
+ const double m = 0.5;
+ const double c1 = 0.;
+ const double c2 = 0.;
+ const double factor = (m / std::sqrt(1.-m*m) *
+ std::sin(std::sqrt(1.-m*m)*t+c2));
+ double result = 1.;
+ for (unsigned int d=0; d<dim; ++d)
+ result *= -4. * std::atan (factor / std::cosh(m*p[d]+c1));
+ return result;
+ }
+
+
+
+
+ template <int dim>
+ class SineGordonProblem
+ {
+ public:
+ SineGordonProblem ();
+ void run ();
+
+ private:
+ void make_grid_and_dofs ();
+ void oldstyle_operation ();
+ void assemble_system ();
+ void output_results (const unsigned int timestep_number) const;
+
+ parallel::distributed::Triangulation<dim> triangulation;
+ FE_Q<dim> fe;
+ DoFHandler<dim> dof_handler;
+ ConstraintMatrix constraints;
+ IndexSet locally_relevant_dofs;
+
+ MatrixFree<dim,double> matrix_free_data;
+
+ parallel::distributed::Vector<double> solution, old_solution, old_old_solution;
+
+ const unsigned int n_global_refinements;
+ double time, time_step;
+ const double final_time;
+ const double cfl_number;
+ const unsigned int output_timestep_skip;
+ };
+
+
+
+ template <int dim>
+ SineGordonProblem<dim>::SineGordonProblem ()
+ :
+ triangulation (MPI_COMM_WORLD),
+ fe (QGaussLobatto<1>(fe_degree+1)),
+ dof_handler (triangulation),
+ n_global_refinements (8-2*dim),
+ time (-10),
+ final_time (-9),
+ cfl_number (.1/fe_degree),
+ output_timestep_skip (200)
+ {}
+
+
+ template <int dim>
+ void SineGordonProblem<dim>::make_grid_and_dofs ()
+ {
+ GridGenerator::hyper_cube (triangulation, -15, 15);
+ triangulation.refine_global (n_global_refinements);
+ {
+ typename Triangulation<dim>::active_cell_iterator
+ cell = triangulation.begin_active(),
+ end_cell = triangulation.end();
+ for ( ; cell != end_cell; ++cell)
+ if (cell->is_locally_owned())
+ if (cell->center().norm() < 11)
+ cell->set_refine_flag();
+ triangulation.execute_coarsening_and_refinement();
+
+ cell = triangulation.begin_active();
+ end_cell = triangulation.end();
+ for ( ; cell != end_cell; ++cell)
+ if (cell->is_locally_owned())
+ if (cell->center().norm() < 6)
+ cell->set_refine_flag();
+ triangulation.execute_coarsening_and_refinement();
+ }
+
+ deallog << " Number of global active cells: "
+ << triangulation.n_global_active_cells()
+ << std::endl;
+
+ dof_handler.distribute_dofs (fe);
+
+ deallog << " Number of degrees of freedom: "
+ << dof_handler.n_dofs()
+ << std::endl;
+
+
+ DoFTools::extract_locally_relevant_dofs (dof_handler,
+ locally_relevant_dofs);
+ constraints.clear();
+ constraints.reinit (locally_relevant_dofs);
+ DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+ constraints.close();
+
+ QGaussLobatto<1> quadrature (fe_degree+1);
+ typename MatrixFree<dim>::AdditionalData additional_data;
+ additional_data.mpi_communicator = MPI_COMM_WORLD;
+ additional_data.tasks_parallel_scheme =
+ MatrixFree<dim>::AdditionalData::partition_partition;
+
+ matrix_free_data.reinit (dof_handler, constraints,
+ quadrature, additional_data);
+
+ matrix_free_data.initialize_dof_vector (solution);
+ old_solution.reinit (solution);
+ old_old_solution.reinit (solution);
+ }
+
+
+
+
+ template <int dim>
+ void
+ SineGordonProblem<dim>::output_results (const unsigned int timestep_number) const
+ {
+ parallel::distributed::Vector<double> locally_relevant_solution;
+ locally_relevant_solution.reinit (dof_handler.locally_owned_dofs(),
+ locally_relevant_dofs,
+ MPI_COMM_WORLD);
+ locally_relevant_solution.copy_from (solution);
+ locally_relevant_solution.update_ghost_values ();
+ constraints.distribute (locally_relevant_solution);
+
+ Vector<float> norm_per_cell (triangulation.n_active_cells());
+ VectorTools::integrate_difference (dof_handler,
+ locally_relevant_solution,
+ ZeroFunction<dim>(),
+ norm_per_cell,
+ QGauss<dim>(fe_degree+1),
+ VectorTools::L2_norm);
+ const double solution_norm =
+ std::sqrt(Utilities::MPI::sum (norm_per_cell.norm_sqr(), MPI_COMM_WORLD));
+
+ deallog << " Time:"
+ << std::setw(8) << std::setprecision(3) << time
+ << ", solution norm: "
+ << std::setprecision(5) << std::setw(7) << solution_norm
+ << std::endl;
+ }
+
+
+
+ template <int dim>
+ void
+ SineGordonProblem<dim>::run ()
+ {
+ make_grid_and_dofs();
+
+ const double local_min_cell_diameter =
+ triangulation.last()->diameter()/std::sqrt(dim);
+ const double global_min_cell_diameter
+ = -Utilities::MPI::max(-local_min_cell_diameter, MPI_COMM_WORLD);
+ time_step = cfl_number * global_min_cell_diameter;
+ time_step = (final_time-time)/(int((final_time-time)/time_step));
+ deallog << " Time step size: " << time_step << ", finest cell: "
+ << global_min_cell_diameter << std::endl << std::endl;
+
+
+ VectorTools::interpolate (dof_handler,
+ ExactSolution<dim> (1, time),
+ solution);
+ VectorTools::interpolate (dof_handler,
+ ExactSolution<dim> (1, time-time_step),
+ old_solution);
+ output_results (0);
+
+ std::vector<parallel::distributed::Vector<double>*> previous_solutions;
+ previous_solutions.push_back(&old_solution);
+ previous_solutions.push_back(&old_old_solution);
+
+ SineGordonOperation<dim,fe_degree> sine_gordon_op (matrix_free_data,
+ time_step);
+
+ unsigned int timestep_number = 1;
+
+ for (time+=time_step; time<=final_time; time+=time_step, ++timestep_number)
+ {
+ old_old_solution.swap (old_solution);
+ old_solution.swap (solution);
+ sine_gordon_op.apply (solution, previous_solutions);
+
+ if (timestep_number % output_timestep_skip == 0)
+ output_results(timestep_number / output_timestep_skip);
+ }
+ output_results(timestep_number / output_timestep_skip + 1);
+
+ deallog << std::endl
+ << " Performed " << timestep_number << " time steps."
+ << std::endl << std::endl;
+ }
+}
+
+
+
+int main (int argc, char ** argv)
+{
+ Utilities::System::MPI_InitFinalize mpi_initialization(argc, argv);
+
+ unsigned int myid = Utilities::System::get_this_mpi_process (MPI_COMM_WORLD);
+ deallog.push(Utilities::int_to_string(myid));
+
+ if (myid == 0)
+ {
+ std::ofstream logfile(output_file_for_mpi("step-48").c_str());
+ deallog.attach(logfile);
+ deallog << std::setprecision(4);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ {
+ deallog.push("2d");
+ Step48::SineGordonProblem<2> sg_problem;
+ sg_problem.run ();
+ deallog.pop();
+ }
+ {
+ deallog.push("3d");
+ Step48::SineGordonProblem<3> sg_problem;
+ sg_problem.run ();
+ deallog.pop();
+ }
+ }
+ else
+ {
+ deallog.depth_console(0);
+ {
+ Step48::SineGordonProblem<2> sg_problem;
+ sg_problem.run ();
+ }
+ {
+ Step48::SineGordonProblem<3> sg_problem;
+ sg_problem.run ();
+ }
+ }
+}
+
projects / dealii-svn.git / commitdiff