+++ /dev/null
-##
-# CMake script for the step-47 tutorial program:
-##
-
-# Set the name of the project and target:
-SET(TARGET "step-47")
-
-# Declare all source files the target consists of. Here, this is only
-# the one step-X.cc file, but as you expand your project you may wish
-# to add other source files as well. If your project becomes much larger,
-# you may want to either replace the following statement by something like
-# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
-# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
-# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
-# or switch altogether to the large project CMakeLists.txt file discussed
-# in the "CMake in user projects" page accessible from the "User info"
-# page of the documentation.
-SET(TARGET_SRC
- ${TARGET}.cc
- )
-
-# Usually, you will not need to modify anything beyond this point...
-
-CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
-
-FIND_PACKAGE(deal.II 8.5.0 QUIET
- HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
- )
-IF(NOT ${deal.II_FOUND})
- MESSAGE(FATAL_ERROR "\n"
- "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
- "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
- "or set an environment variable \"DEAL_II_DIR\" that contains this path."
- )
-ENDIF()
-
-DEAL_II_INITIALIZE_CACHED_VARIABLES()
-PROJECT(${TARGET})
-DEAL_II_INVOKE_AUTOPILOT()
+++ /dev/null
-<a name="Intro"></a>
-<h1>Introduction</h1>
-
+++ /dev/null
-unfinished
+++ /dev/null
-#xfem quadrature Points
-0.781998 0.781998
-0.941587 0.72508
-0.72508 0.941587
-0.926335 0.926335
-0.397329 0.926335
-0.583333 0.72508
-0.583333 0.941587
-0.686004 0.781998
-0.72508 0.583333
-0.926335 0.397329
-0.781998 0.686004
-0.941587 0.583333
-0.115331 0.315492
-0.430422 0.60008
-0.069578 0.816587
-0.259669 0.892842
-0.315492 0.115331
-0.816587 0.069578
-0.60008 0.430422
-0.892842 0.259669
+++ /dev/null
-<h1>Results</h1>
+++ /dev/null
-#vertices of xfem subcells
-0 0
-1 0
-1 1
-0 1
-
-0.25 1
-1 0.25
-
-0.625 0.625
-0.75 0.75
-
-0.75 0.75
-0.625 1
-
-0.75 0.75
-1 0.625
-
-0 0
-0.625 0.625
+++ /dev/null
-/* ---------------------------------------------------------------------
- *
- * Copyright (C) 2011 - 2017 by the deal.II authors
- *
- * This file is part of the deal.II library.
- *
- * The deal.II library is free software; you can use it, redistribute
- * it, and/or modify it under the terms of the GNU Lesser General
- * Public License as published by the Free Software Foundation; either
- * version 2.1 of the License, or (at your option) any later version.
- * The full text of the license can be found in the file LICENSE at
- * the top level of the deal.II distribution.
- *
- * ---------------------------------------------------------------------
-
- *
- * Author: Wolfgang Bangerth, University of Heidelberg, 2000
- */
-
-
-#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/sparse_matrix.h>
-#include <deal.II/lac/dynamic_sparsity_pattern.h>
-#include <deal.II/lac/solver_cg.h>
-#include <deal.II/lac/precondition.h>
-#include <deal.II/grid/tria.h>
-#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/tria_accessor.h>
-#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/manifold_lib.h>
-#include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/dofs/dof_tools.h>
-#include <deal.II/fe/fe_q.h>
-#include <deal.II/fe/fe_nothing.h>
-#include <deal.II/fe/fe_system.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/numerics/matrix_tools.h>
-#include <deal.II/numerics/data_out.h>
-
-#include <fstream>
-#include <iostream>
-
-#include <deal.II/fe/fe_q.h>
-#include <deal.II/grid/grid_out.h>
-
-
-#include <deal.II/lac/constraint_matrix.h>
-
-#include <deal.II/grid/grid_refinement.h>
-
-#include <deal.II/numerics/error_estimator.h>
-
-namespace Step47
-{
- using namespace dealii;
-
-
-
- double sign (double d)
- {
- if (d > 0)
- return 1;
- else if (d < 0)
- return -1;
- else
- return 0;
- }
-
-
- template <int dim>
- class LaplaceProblem
- {
- public:
- LaplaceProblem ();
- ~LaplaceProblem ();
-
- void run ();
-
- private:
- bool interface_intersects_cell (const typename Triangulation<dim>::cell_iterator &cell) const;
- std::pair<unsigned int, Quadrature<dim> > compute_quadrature(const Quadrature<dim> &plain_quadrature, const typename hp::DoFHandler<dim>::active_cell_iterator &cell, const std::vector<double> &level_set_values);
- void append_quadrature(const Quadrature<dim> &plain_quadrature,
- const std::vector<Point<dim> > &v ,
- std::vector<Point<dim> > &xfem_points,
- std::vector<double> &xfem_weights);
-
- void setup_system ();
- void assemble_system ();
- void solve ();
- void refine_grid ();
- void output_results (const unsigned int cycle) const;
- void compute_error () const;
-
- Triangulation<dim> triangulation;
-
- hp::DoFHandler<dim> dof_handler;
- hp::FECollection<dim> fe_collection;
-
- ConstraintMatrix constraints;
-
- SparsityPattern sparsity_pattern;
- SparseMatrix<double> system_matrix;
-
- Vector<double> solution;
- Vector<double> system_rhs;
- };
-
-
-
-
- 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;
-
- virtual void value_list (const std::vector<Point<dim> > &points,
- std::vector<double> &values,
- const unsigned int component = 0) const;
- };
-
-
-
- template <int dim>
- double Coefficient<dim>::value (const Point<dim> &p,
- const unsigned int) const
- {
- if (p.square() < 0.5*0.5)
- return 20;
- else
- return 1;
- }
-
-
-
- template <int dim>
- void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points,
- std::vector<double> &values,
- const unsigned int component) const
- {
- const unsigned int n_points = points.size();
-
- Assert (values.size() == n_points,
- ExcDimensionMismatch (values.size(), n_points));
-
- (void) component;
- Assert(component == 0, ExcIndexRange(component, 0, 1));
-
- for (unsigned int i=0; i<n_points; ++i)
- {
- if (points[i].square() < 0.5*0.5)
- values[i] = 20;
- else
- values[i] = 1;
- }
- }
-
-
-
- template <int dim>
- double exact_solution (const Point<dim> &p)
- {
- const double r = p.norm();
-
- return (r < 0.5
- ?
- 1./20 * (-1./4*r*r + 61./16)
- :
- 1./4 * (1-r*r));
- }
-
-
- template <int dim>
- LaplaceProblem<dim>::LaplaceProblem ()
- :
- dof_handler (triangulation)
- {
- fe_collection.push_back (FESystem<dim> (FE_Q<dim>(1), 1,
- FE_Nothing<dim>(), 1));
- fe_collection.push_back (FESystem<dim> (FE_Q<dim>(1), 1,
- FE_Q<dim>(1), 1));
- }
-
-
-
- template <int dim>
- LaplaceProblem<dim>::~LaplaceProblem ()
- {
- dof_handler.clear ();
- }
-
-
-
- template <int dim>
- double
- level_set (const Point<dim> &p)
- {
- return p.norm() - 0.5;
- }
-
-
-
- template <int dim>
- Tensor<1,dim>
- grad_level_set (const Point<dim> &p)
- {
- return p / p.norm();
- }
-
-
-
- template <int dim>
- bool
- LaplaceProblem<dim>::
- interface_intersects_cell (const typename Triangulation<dim>::cell_iterator &cell) const
- {
- for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell-1; ++v)
- if (level_set(cell->vertex(v)) * level_set(cell->vertex(v+1)) < 0)
- return true;
-
- // we get here only if all vertices have the same sign, which means that
- // the cell is not intersected
- return false;
- }
-
-
-
- template <int dim>
- void LaplaceProblem<dim>::setup_system ()
- {
- for (typename hp::DoFHandler<dim>::cell_iterator cell
- = dof_handler.begin_active();
- cell != dof_handler.end(); ++cell)
- if (interface_intersects_cell(cell) == false)
- cell->set_active_fe_index(0);
- else
- cell->set_active_fe_index(1);
-
- dof_handler.distribute_dofs (fe_collection);
-
- solution.reinit (dof_handler.n_dofs());
- system_rhs.reinit (dof_handler.n_dofs());
-
-
- constraints.clear ();
-//TODO: fix this, it currently crashes
- // DoFTools::make_hanging_node_constraints (dof_handler, constraints);
-
-//TODO: component 1 must satisfy zero boundary conditions
- constraints.close();
-
-
- DynamicSparsityPattern dsp(dof_handler.n_dofs());
- DoFTools::make_sparsity_pattern (dof_handler, dsp);
-
- constraints.condense (dsp);
-
- sparsity_pattern.copy_from(dsp);
-
- system_matrix.reinit (sparsity_pattern);
- }
-
-
- template <int dim>
- void LaplaceProblem<dim>::assemble_system ()
- {
- const QGauss<dim> quadrature_formula(3);
-
-
- FEValues<dim> plain_fe_values (fe_collection[0], quadrature_formula,
- update_values | update_gradients |
- update_quadrature_points | update_JxW_values);
-
- const unsigned int n_q_points = quadrature_formula.size();
-
- FullMatrix<double> cell_matrix;
- Vector<double> cell_rhs;
-
- std::vector<types::global_dof_index> local_dof_indices;
-
- const Coefficient<dim> coefficient;
- std::vector<double> coefficient_values (n_q_points);
-
- typename hp::DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active(),
- endc = dof_handler.end();
-
- for (; cell!=endc; ++cell)
- {
- const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
- cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
- cell_rhs.reinit (dofs_per_cell);
-
- cell_matrix = 0;
- cell_rhs = 0;
-
- if (cell->active_fe_index() == 0)
- {
- plain_fe_values.reinit (cell);
-
- coefficient_values.resize (plain_fe_values.n_quadrature_points);
- coefficient.value_list (plain_fe_values.get_quadrature_points(),
- coefficient_values);
-
- for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- cell_matrix(i,j) += (coefficient_values[q_point] *
- plain_fe_values.shape_grad(i,q_point) *
- plain_fe_values.shape_grad(j,q_point) *
- plain_fe_values.JxW(q_point));
-
-
- cell_rhs(i) += (plain_fe_values.shape_value(i,q_point) *
- 1.0 *
- plain_fe_values.JxW(q_point));
- }
- }
- else
- {
-//TODO: verify that the order of support points equals the order of vertices
-//of the cells, as we use below
- Assert (cell->active_fe_index() == 1, ExcInternalError());
- Assert (interface_intersects_cell(cell) == true, ExcInternalError());
-
- std::vector<double> level_set_values (GeometryInfo<dim>::vertices_per_cell);
- for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
- level_set_values[v] = level_set (cell->vertex(v));
-
- FEValues<dim> this_fe_values (fe_collection[1],
- compute_quadrature(quadrature_formula, cell,
- level_set_values).second,
- update_values | update_gradients |
- update_quadrature_points | update_JxW_values );
-
- this_fe_values.reinit (cell);
-
- coefficient_values.resize (this_fe_values.n_quadrature_points);
- coefficient.value_list (this_fe_values.get_quadrature_points(),
- coefficient_values);
-
- for (unsigned int q_point=0; q_point<this_fe_values.n_quadrature_points; ++q_point)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- if (cell->get_fe().system_to_component_index(i).first == 0)
- {
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- if (cell->get_fe().system_to_component_index(j).first == 0)
- cell_matrix(i,j) += (coefficient_values[q_point] *
- this_fe_values.shape_grad(i,q_point) *
- this_fe_values.shape_grad(j,q_point) *
- this_fe_values.JxW(q_point));
- else
- cell_matrix(i,j) += (coefficient_values[q_point] *
- this_fe_values.shape_grad(i,q_point)
- *
- ((std::fabs(level_set(this_fe_values.quadrature_point(q_point)))
- -
- std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(j).second))))*
- this_fe_values.shape_grad(j,q_point)
- +
- grad_level_set(this_fe_values.quadrature_point(q_point)) *
- sign(level_set(this_fe_values.quadrature_point(q_point))) *
- this_fe_values.shape_value(j,q_point)) *
- this_fe_values.JxW(q_point));
-
- cell_rhs(i) += (this_fe_values.shape_value(i,q_point) *
- 1.0 *
- this_fe_values.JxW(q_point));
- }
- else
- {
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- if (cell->get_fe().system_to_component_index(j).first == 0)
- cell_matrix(i,j) += (coefficient_values[q_point] *
- ((std::fabs(level_set(this_fe_values.quadrature_point(q_point)))
- -
- std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(i).second))))*
- this_fe_values.shape_grad(i,q_point)
- +
- grad_level_set(this_fe_values.quadrature_point(q_point)) *
- sign(level_set(this_fe_values.quadrature_point(q_point))) *
- this_fe_values.shape_value(i,q_point)) *
- this_fe_values.shape_grad(j,q_point) *
- this_fe_values.JxW(q_point));
- else
- cell_matrix(i,j) += (coefficient_values[q_point] *
- ((std::fabs(level_set(this_fe_values.quadrature_point(q_point)))
- -
- std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(i).second))))*
- this_fe_values.shape_grad(i,q_point)
- +
- grad_level_set(this_fe_values.quadrature_point(q_point)) *
- sign(level_set(this_fe_values.quadrature_point(q_point))) *
- this_fe_values.shape_value(i,q_point)) *
- ((std::fabs(level_set(this_fe_values.quadrature_point(q_point)))
- -
- std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(j).second))))*
- this_fe_values.shape_grad(j,q_point)
- +
- grad_level_set(this_fe_values.quadrature_point(q_point)) *
- sign(level_set(this_fe_values.quadrature_point(q_point))) *
- this_fe_values.shape_value(j,q_point)) *
- this_fe_values.JxW(q_point));
-
- cell_rhs(i) += ((std::fabs(level_set(this_fe_values.quadrature_point(q_point)))
- -
- std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(i).second))))*
- this_fe_values.shape_value(i,q_point) *
- 1.0 *
- this_fe_values.JxW(q_point));
- }
- }
-
- local_dof_indices.resize (dofs_per_cell);
- cell->get_dof_indices (local_dof_indices);
- constraints.distribute_local_to_global (cell_matrix, cell_rhs,
- local_dof_indices,
- system_matrix, system_rhs);
- }
-
-
- std::map<types::global_dof_index,double> boundary_values;
- VectorTools::interpolate_boundary_values (dof_handler,
- 0,
- ZeroFunction<dim>(2),
- boundary_values);
- MatrixTools::apply_boundary_values (boundary_values,
- system_matrix,
- solution,
- system_rhs);
-
- }
-
-// To integrate the enriched elements we have to find the geometrical
-// decomposition of the original element in subelements. The subelements are
-// used to integrate the elements on both sides of the discontinuity. The
-// discontinuity line is approximated by a piece-wise linear interpolation
-// between the intersection of the discontinuity with the edges of the
-// elements. The vector level_set_values has the values of the level set
-// function at the vertices of the elements. From these values can be found by
-// linear interpolation the intersections. There are three kind of
-// decomposition that are considered. Type 1: there is not cut. Type 2: a
-// corner of the element is cut. Type 3: two corners are cut.
-
- template <int dim>
- std::pair<unsigned int, Quadrature<dim> >
- LaplaceProblem<dim>::compute_quadrature (const Quadrature<dim> &plain_quadrature,
- const typename hp::DoFHandler<dim>::active_cell_iterator &/*cell*/,
- const std::vector<double> &level_set_values)
- {
-
- unsigned int type = 0;
-
- // find the type of cut
- int sign_ls[GeometryInfo<dim>::vertices_per_cell];
- for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
- {
- if (level_set_values[v] > 0) sign_ls[v] = 1;
- else if (level_set_values[v] < 0) sign_ls[v] = -1;
- else sign_ls[v] = 0;
- }
-
- // the sign of the level set function at the 4 nodes of the elements can
- // be positive + or negative - depending on the sign of the level set
- // function we have the following three classes of decomposition type 1:
- // ++++, ---- type 2: -+++, +-++, ++-+, +++-, +---, -+--, --+-, ---+ type
- // 3: +--+, ++--, +-+-, -++-, --++, -+-+
-
- if ( sign_ls[0]==sign_ls[1] &&
- sign_ls[0]==sign_ls[2] &&
- sign_ls[0]==sign_ls[3] )
- type = 1;
- else if ( sign_ls[0]*sign_ls[1]*sign_ls[2]*sign_ls[3] < 0 )
- type = 2;
- else
- type = 3;
-
- unsigned int Pos = 100;
-
- Point<dim> v0(0,0);
- Point<dim> v1(1,0);
- Point<dim> v2(0,1);
- Point<dim> v3(1,1);
-
- Point<dim> A(0,0);
- Point<dim> B(0,0);
- Point<dim> C(0,0);
- Point<dim> D(0,0);
- Point<dim> E(0,0);
- Point<dim> F(0,0);
-
- if (type == 1)
- return std::pair<unsigned int, Quadrature<dim> >(1, plain_quadrature);
-
- if (type==2)
- {
- const unsigned int n_q_points = plain_quadrature.size();
-
- // loop over all subelements for integration in type 2 there are 5
- // subelements
-
- Quadrature<dim> xfem_quadrature(5*n_q_points);
-
- std::vector<Point<dim> > v(GeometryInfo<dim>::vertices_per_cell);
-
- if (sign_ls[0]!=sign_ls[1] && sign_ls[0]!=sign_ls[2] && sign_ls[0]!=sign_ls[3]) Pos = 0;
- else if (sign_ls[1]!=sign_ls[0] && sign_ls[1]!=sign_ls[2] && sign_ls[1]!=sign_ls[3]) Pos = 1;
- else if (sign_ls[2]!=sign_ls[0] && sign_ls[2]!=sign_ls[1] && sign_ls[2]!=sign_ls[3]) Pos = 2;
- else if (sign_ls[3]!=sign_ls[0] && sign_ls[3]!=sign_ls[1] && sign_ls[3]!=sign_ls[2]) Pos = 3;
- else assert(0); // error message
-
- // Find cut coordinates
-
- // deal.ii local coordinates
-
- // 2-------3 | | | | | | 0-------1
-
- if (Pos == 0)
- {
- A[0] = 1. - level_set_values[1]/(level_set_values[1]-level_set_values[0]);
- B[1] = 1. - level_set_values[2]/(level_set_values[2]-level_set_values[0]);
- A(1) = 0.;
- B(0) = 0.;
- C(0) = 0.5*( A(0) + B(0) );
- C(1) = 0.5*( A(1) + B(1) );
- D(0) = 2./3. * C(0);
- D(1) = 2./3. * C(1);
- E(0) = 0.5*A(0);
- E(1) = 0.;
- F(0) = 0.;
- F(1) = 0.5*B(1);
- }
- else if (Pos == 1)
- {
- A[0] = level_set_values[0]/(level_set_values[0]-level_set_values[1]);
- B[1] = 1 - level_set_values[3]/(level_set_values[3]-level_set_values[1]);
- A(1) = 0.;
- B(0) = 1.;
- C(0) = 0.5*( A(0) + B(0) );
- C(1) = 0.5*( A(1) + B(1) );
- D(0) = 1./3. + 2./3. * C(0);
- D(1) = 2./3. * C(1);
- E(0) = 0.5*(1 + A(0));
- E(1) = 0.;
- F(0) = 1.;
- F(1) = 0.5*B(1);
- }
- else if (Pos == 2)
- {
- A[0] = 1 - level_set_values[3]/(level_set_values[3]-level_set_values[2]);
- B[1] = level_set_values[0]/(level_set_values[0]-level_set_values[2]);
- A(1) = 1.;
- B(0) = 0.;
- C(0) = 0.5*( A(0) + B(0) );
- C(1) = 0.5*( A(1) + B(1) );
- D(0) = 2./3. * C(0);
- D(1) = 1./3. + 2./3. * C(1);
- E(0) = 0.5* A(0);
- E(1) = 1.;
- F(0) = 0.;
- F(1) = 0.5*( 1. + B(1) );
- }
- else if (Pos == 3)
- {
- A[0] = level_set_values[2]/(level_set_values[2]-level_set_values[3]);
- B[1] = level_set_values[1]/(level_set_values[1]-level_set_values[3]);
- A(1) = 1.;
- B(0) = 1.;
- C(0) = 0.5*( A(0) + B(0) );
- C(1) = 0.5*( A(1) + B(1) );
- D(0) = 1./3. + 2./3. * C(0);
- D(1) = 1./3. + 2./3. * C(1);
- E(0) = 0.5*( 1. + A(0) );
- E(1) = 1.;
- F(0) = 1.;
- F(1) = 0.5*( 1. + B(1) );
- }
-
- //std::cout << A << std::endl; std::cout << B << std::endl; std::cout
- //<< C << std::endl; std::cout << D << std::endl; std::cout << E <<
- //std::endl; std::cout << F << std::endl;
-
- std::string filename = "vertices.dat";
- std::ofstream output (filename.c_str());
- output << "#vertices of xfem subcells" << std::endl;
- output << v0(0) << " " << v0(1) << std::endl;
- output << v1(0) << " " << v1(1) << std::endl;
- output << v3(0) << " " << v3(1) << std::endl;
- output << v2(0) << " " << v2(1) << std::endl;
- output << std::endl;
- output << A(0) << " " << A(1) << std::endl;
- output << B(0) << " " << B(1) << std::endl;
- output << std::endl;
- output << C(0) << " " << C(1) << std::endl;
- output << D(0) << " " << D(1) << std::endl;
- output << std::endl;
- output << D(0) << " " << D(1) << std::endl;
- output << E(0) << " " << E(1) << std::endl;
- output << std::endl;
- output << D(0) << " " << D(1) << std::endl;
- output << F(0) << " " << F(1) << std::endl;
- output << std::endl;
-
- if (Pos==0)
- output << v3(0) << " " << v3(1) << std::endl;
- else if (Pos==1)
- output << v2(0) << " " << v2(1) << std::endl;
- else if (Pos==2)
- output << v1(0) << " " << v1(1) << std::endl;
- else if (Pos==3)
- output << v0(0) << " " << v0(1) << std::endl;
- output << C(0) << " " << C(1) << std::endl;
-
- Point<dim> subcell_vertices[10];
- subcell_vertices[0] = v0;
- subcell_vertices[1] = v1;
- subcell_vertices[2] = v2;
- subcell_vertices[3] = v3;
- subcell_vertices[4] = A;
- subcell_vertices[5] = B;
- subcell_vertices[6] = C;
- subcell_vertices[7] = D;
- subcell_vertices[8] = E;
- subcell_vertices[9] = F;
-
- std::vector<Point<dim> > xfem_points;
- std::vector<double> xfem_weights;
-
- // lookup table for the decomposition
-
- if (dim==2)
- {
- unsigned int subcell_v_indices[4][5][4] =
- {
- {{0,8,9,7}, {9,7,5,6}, {8,4,7,6}, {5,6,2,3}, {6,4,3,1}},
- {{8,1,7,9}, {4,8,6,7}, {6,7,5,9}, {0,4,2,6}, {2,6,3,5}},
- {{9,7,2,8}, {5,6,9,7}, {6,4,7,8}, {0,1,5,6}, {6,1,4,3}},
- {{7,9,8,3}, {4,6,8,7}, {6,5,7,9}, {0,6,2,4}, {0,1,6,5}}
- };
-
- for (unsigned int subcell = 0; subcell<5; subcell++)
- {
- //std::cout << "subcell : " << subcell << std::endl;
- std::vector<Point<dim> > vertices;
- for (unsigned int i=0; i<4; i++)
- {
- vertices.push_back( subcell_vertices[subcell_v_indices[Pos][subcell][i]] );
- //std::cout << "i : " << i << std::endl; std::cout <<
- //"subcell v : " << subcell_v_indices[Pos][subcell][i] <<
- //std::endl; std::cout << vertices[i](0) << " " <<
- //vertices[i](1) << std::endl;
- }
- //std::cout << std::endl; create quadrature rule
- append_quadrature( plain_quadrature,
- vertices,
- xfem_points,
- xfem_weights);
- //initialize xfem_quadrature with quadrature points of all
- //subelements
- xfem_quadrature.initialize(xfem_points, xfem_weights);
- }
- }
-
- Assert (xfem_quadrature.size() == plain_quadrature.size() * 5, ExcInternalError());
- return std::pair<unsigned int, Quadrature<dim> >(2, xfem_quadrature);
- }
-
- // Type three decomposition (+--+, ++--, +-+-, -++-, --++, -+-+)
-
- if (type==3)
- {
- const unsigned int n_q_points = plain_quadrature.size();
-
- // loop over all subelements for integration in type 2 there are 5
- // subelements
-
- Quadrature<dim> xfem_quadrature(5*n_q_points);
-
- std::vector<Point<dim> > v(GeometryInfo<dim>::vertices_per_cell);
-
- if ( sign_ls[0]==sign_ls[1] && sign_ls[2]==sign_ls[3] )
- {
- Pos = 0;
- A(0) = 0.;
- A(1) = level_set_values[0]/((level_set_values[0]-level_set_values[2]));
- B(0) = 1.;
- B(1) = level_set_values[1]/((level_set_values[1]-level_set_values[3]));
- }
- else if ( sign_ls[0]==sign_ls[2] && sign_ls[1]==sign_ls[3] )
- {
- Pos = 1;
- A(0) = level_set_values[0]/((level_set_values[0]-level_set_values[1]));
- A(1) = 0.;
- B(0) = level_set_values[2]/((level_set_values[2]-level_set_values[3]));
- B(1) = 1.;
- }
- else if ( sign_ls[0]==sign_ls[3] && sign_ls[1]==sign_ls[2] )
- {
- std::cout << "Error: the element has two cut lines and this is not allowed" << std::endl;
- assert(0);
- }
- else
- {
- std::cout << "Error: the level set function has not the right values" << std::endl;
- assert(0);
- }
-
- //std::cout << "Pos " << Pos << std::endl; std::cout << A <<
- //std::endl; std::cout << B << std::endl;
- std::string filename = "vertices.dat";
- std::ofstream output (filename.c_str());
- output << "#vertices of xfem subcells" << std::endl;
- output << A(0) << " " << A(1) << std::endl;
- output << B(0) << " " << B(1) << std::endl;
-
- //fill xfem_quadrature
- Point<dim> subcell_vertices[6];
- subcell_vertices[0] = v0;
- subcell_vertices[1] = v1;
- subcell_vertices[2] = v2;
- subcell_vertices[3] = v3;
- subcell_vertices[4] = A;
- subcell_vertices[5] = B;
-
- std::vector<Point<dim> > xfem_points;
- std::vector<double> xfem_weights;
-
- if (dim==2)
- {
- unsigned int subcell_v_indices[2][2][4] =
- {
- {{0,1,4,5}, {4,5,2,3}},
- {{0,4,2,5}, {4,1,5,3}}
- };
-
- //std::cout << "Pos : " << Pos << std::endl;
- for (unsigned int subcell = 0; subcell<2; subcell++)
- {
- //std::cout << "subcell : " << subcell << std::endl;
- std::vector<Point<dim> > vertices;
- for (unsigned int i=0; i<4; i++)
- {
- vertices.push_back( subcell_vertices[subcell_v_indices[Pos][subcell][i]] );
- //std::cout << "i : " << i << std::endl; std::cout <<
- //"subcell v : " << subcell_v_indices[Pos][subcell][i] <<
- //std::endl; std::cout << vertices[i](0) << " " <<
- //vertices[i](1) << std::endl;
- }
- //std::cout << std::endl; create quadrature rule
- append_quadrature( plain_quadrature,
- vertices,
- xfem_points,
- xfem_weights);
- //initialize xfem_quadrature with quadrature points of all
- //subelements
- xfem_quadrature.initialize(xfem_points, xfem_weights);
- }
- }
- Assert (xfem_quadrature.size() == plain_quadrature.size() * 2, ExcInternalError());
- return std::pair<unsigned int, Quadrature<dim> >(3, xfem_quadrature);
- }
-
- return std::pair<unsigned int, Quadrature<dim> >(0, plain_quadrature);;
-
- }
-
- template <int dim>
- void LaplaceProblem<dim>::append_quadrature ( const Quadrature<dim> &plain_quadrature,
- const std::vector<Point<dim> > &v,
- std::vector<Point<dim> > &xfem_points,
- std::vector<double> &xfem_weights)
-
- {
- // Project integration points into sub-elements. This maps quadrature
- // points from a reference element to a subelement of a reference element.
- // To implement the action of this map the coordinates of the subelements
- // have been calculated (A(0)...F(0),A(1)...F(1)) the coordinates of the
- // quadrature points are given by the bi-linear map defined by the form
- // functions $x^\prime_i = \sum_j v^\prime \phi_j(x^hat_i)$, where the
- // $\phi_j$ are the shape functions of the FEQ.
-
- unsigned int n_v = GeometryInfo<dim>::vertices_per_cell;
-
- std::vector<Point<dim> > q_points = plain_quadrature.get_points();
- std::vector<Point<dim> > q_transf(q_points.size());
- std::vector<double> W = plain_quadrature.get_weights();
- std::vector<double> phi(n_v);
- std::vector<Tensor<1,dim> > grad_phi(n_v);
-
- const unsigned int n_q_points = plain_quadrature.size();
-
- std::vector<double> JxW(n_q_points);
-
- for ( unsigned int i = 0; i < n_q_points; i++)
- {
- switch (dim)
- {
- case 2:
- {
- double xi = q_points[i](0);
- double eta = q_points[i](1);
-
- // Define shape functions on reference element we consider a
- // bi-linear mapping
- phi[0] = (1. - xi) * (1. - eta);
- phi[1] = xi * (1. - eta);
- phi[2] = (1. - xi) * eta;
- phi[3] = xi * eta;
-
- grad_phi[0][0] = (-1. + eta);
- grad_phi[1][0] = (1. - eta);
- grad_phi[2][0] = -eta;
- grad_phi[3][0] = eta;
-
- grad_phi[0][1] = (-1. + xi);
- grad_phi[1][1] = -xi;
- grad_phi[2][1] = 1-xi;
- grad_phi[3][1] = xi;
-
- break;
- }
-
- default:
- Assert (false, ExcNotImplemented());
- }
-
-
- Tensor<2,dim> jacobian;
-
- // Calculate Jacobian of transformation
- for (unsigned int d=0; d<dim; ++d)
- for (unsigned int e=0; e<dim; ++e)
- {
- for (unsigned int j = 0; j<GeometryInfo<dim>::vertices_per_cell; j++)
- {
- jacobian[d][e] += grad_phi[j][e] * v[j](d);
- }
- }
-
- double detJ = determinant(jacobian);
- xfem_weights.push_back (W[i] * detJ);
-
- // Map integration points from reference element to subcell of
- // reference element
- Point<dim> q_prime;
- for (unsigned int d=0; d<dim; ++d)
- for (unsigned int j = 0; j<GeometryInfo<dim>::vertices_per_cell; j++)
- q_prime[d] += v[j](d) * phi[j];
- xfem_points.push_back(q_prime);
- }
-
- }
-
-
- template <int dim>
- void LaplaceProblem<dim>::solve ()
- {
- SolverControl solver_control (1000, 1e-12);
- SolverCG<> solver (solver_control);
-
- PreconditionSSOR<> preconditioner;
- preconditioner.initialize(system_matrix, 1.2);
-
- solver.solve (system_matrix, solution, system_rhs,
- preconditioner);
-
- constraints.distribute (solution);
- }
-
-
-
- template <int dim>
- void LaplaceProblem<dim>::refine_grid ()
- {
- Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
-
- KellyErrorEstimator<dim>::estimate (dof_handler,
- QGauss<dim-1>(3),
- typename FunctionMap<dim>::type(),
- solution,
- estimated_error_per_cell);
-
- GridRefinement::refine_and_coarsen_fixed_number (triangulation,
- estimated_error_per_cell,
- 0.3, 0.03);
-
- triangulation.execute_coarsening_and_refinement ();
- }
-
-
-
- template <int dim>
- class Postprocessor : public DataPostprocessor<dim>
- {
- public:
- virtual
- void
- evaluate_vector_field
- (const dealii::DataPostprocessorInputs::Vector<dim> &inputs,
- std::vector<Vector<double> > &computed_quantities) const;
-
- virtual std::vector<std::string> get_names () const;
-
- virtual
- std::vector<DataComponentInterpretation::DataComponentInterpretation>
- get_data_component_interpretation () const;
-
- virtual UpdateFlags get_needed_update_flags () const;
- };
-
-
- template <int dim>
- std::vector<std::string>
- Postprocessor<dim>::get_names() const
- {
- std::vector<std::string> solution_names (1, "total_solution");
- solution_names.push_back ("error");
- return solution_names;
- }
-
-
- template <int dim>
- std::vector<DataComponentInterpretation::DataComponentInterpretation>
- Postprocessor<dim>::
- get_data_component_interpretation () const
- {
- std::vector<DataComponentInterpretation::DataComponentInterpretation>
- interpretation (2,
- DataComponentInterpretation::component_is_scalar);
- return interpretation;
- }
-
-
- template <int dim>
- UpdateFlags
- Postprocessor<dim>::get_needed_update_flags() const
- {
- return update_values | update_q_points;
- }
-
-
- template <int dim>
- void
- Postprocessor<dim>::
- evaluate_vector_field
- (const dealii::DataPostprocessorInputs::Vector<dim> &inputs,
- std::vector<Vector<double> > &computed_quantities) const
- {
- const unsigned int n_quadrature_points = inputs.solution_values.size();
- Assert (computed_quantities.size() == n_quadrature_points,
- ExcInternalError());
- Assert (inputs.solution_values[0].size() == 2,
- ExcInternalError());
-
- for (unsigned int q=0; q<n_quadrature_points; ++q)
- {
- computed_quantities[q](0)
- = (inputs.solution_values[q](0)
- +
-//TODO: shift in weight function is missing!
- inputs.solution_values[q](1) * std::fabs(level_set(inputs.evaluation_points[q])));
- computed_quantities[q](1)
- = (computed_quantities[q](0)
- -
- exact_solution (inputs.evaluation_points[q]));
- }
- }
-
-
-
- template <int dim>
- void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
- {
- Assert (cycle < 10, ExcNotImplemented());
-
- std::string filename = "solution-";
- filename += ('0' + cycle);
- filename += ".vtk";
-
- std::ofstream output (filename.c_str());
-
- Postprocessor<dim> postprocessor;
- DataOut<dim,hp::DoFHandler<dim> > data_out;
-
- data_out.attach_dof_handler (dof_handler);
- data_out.add_data_vector (solution, "solution");
- data_out.add_data_vector (solution, postprocessor);
- data_out.build_patches (5);
-
- data_out.write_vtk (output);
- }
-
-
-
- template <int dim>
- void LaplaceProblem<dim>::compute_error () const
- {
- hp::QCollection<dim> q_collection;
- q_collection.push_back (QGauss<dim>(2));
- q_collection.push_back (QIterated<dim>(QGauss<1>(2), 4));
-
- hp::FEValues<dim> hp_fe_values (fe_collection, q_collection,
- update_values | update_q_points | update_JxW_values);
-
- double l2_error_square = 0;
-
- std::vector<Vector<double> > solution_values;
-
- typename hp::DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active(),
- endc = dof_handler.end();
-
- for (; cell!=endc; ++cell)
- {
- hp_fe_values.reinit (cell);
-
- const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values ();
-
- solution_values.resize (fe_values.n_quadrature_points,
- Vector<double>(2));
- fe_values.get_function_values (solution,
- solution_values);
-
- for (unsigned int q=0; q<fe_values.n_quadrature_points; ++q)
- {
- const double local_error = (solution_values[q](0)
- +
- std::fabs(level_set(fe_values.quadrature_point(q))) *
- solution_values[q](1)
- -
- exact_solution (fe_values.quadrature_point(q)));
- l2_error_square += local_error * local_error * fe_values.JxW(q);
- }
- }
-
- std::cout << " L2 error = " << std::sqrt (l2_error_square)
- << std::endl;
- }
-
-
-
-
- template <int dim>
- void LaplaceProblem<dim>::run ()
- {
- for (unsigned int cycle=0; cycle<6; ++cycle)
- {
- std::cout << "Cycle " << cycle << ':' << std::endl;
-
- if (cycle == 0)
- {
- GridGenerator::hyper_ball (triangulation);
- //GridGenerator::hyper_cube (triangulation, -1, 1);
-
- static const SphericalManifold<dim> boundary;
- triangulation.set_all_manifold_ids_on_boundary(0);
- triangulation.set_manifold (0, boundary);
-
- triangulation.refine_global (2);
- }
- else
- triangulation.refine_global (1);
-// refine_grid ();
-
-
- std::cout << " Number of active cells: "
- << triangulation.n_active_cells()
- << std::endl;
-
- setup_system ();
-
- std::cout << " Number of degrees of freedom: "
- << dof_handler.n_dofs()
- << std::endl;
-
- assemble_system ();
- solve ();
- compute_error ();
- output_results (cycle);
- }
- }
-}
-
-
-
-int main ()
-{
-
- try
- {
- using namespace dealii;
- using namespace Step47;
-
- LaplaceProblem<2> laplace_problem_2d;
- laplace_problem_2d.run ();
- }
- 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;
-}
+++ /dev/null
-##
-# CMake script for the step-50 tutorial program:
-##
-
-# Set the name of the project and target:
-SET(TARGET "step-50")
-
-# Declare all source files the target consists of. Here, this is only
-# the one step-X.cc file, but as you expand your project you may wish
-# to add other source files as well. If your project becomes much larger,
-# you may want to either replace the following statement by something like
-# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
-# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
-# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
-# or switch altogether to the large project CMakeLists.txt file discussed
-# in the "CMake in user projects" page accessible from the "User info"
-# page of the documentation.
-SET(TARGET_SRC
- ${TARGET}.cc
- )
-
-# Define the output that should be cleaned:
-SET(CLEAN_UP_FILES *.vtu *.pvtu *.visit)
-
-# Usually, you will not need to modify anything beyond this point...
-
-CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
-
-FIND_PACKAGE(deal.II 8.5.0 QUIET
- HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
- )
-IF(NOT ${deal.II_FOUND})
- MESSAGE(FATAL_ERROR "\n"
- "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
- "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
- "or set an environment variable \"DEAL_II_DIR\" that contains this path."
- )
-ENDIF()
-
-#
-# Are all dependencies fulfilled?
-#
-IF(NOT DEAL_II_WITH_MPI OR NOT DEAL_II_WITH_P4EST OR NOT DEAL_II_WITH_TRILINOS) # keep in one line
- MESSAGE(FATAL_ERROR "
-Error! The deal.II library found at ${DEAL_II_PATH} was not configured with
- DEAL_II_WITH_MPI = ON
- DEAL_II_WITH_P4EST = ON
- DEAL_II_WITH_TRILINOS = ON
-One or all of these are OFF in your installation but are required for this tutorial step."
- )
-ENDIF()
-
-DEAL_II_INITIALIZE_CACHED_VARIABLES()
-PROJECT(${TARGET})
-DEAL_II_INVOKE_AUTOPILOT()
+++ /dev/null
-step-16 step-40
+++ /dev/null
-<br>
-
-<i>This program has evolved from a version originally written by Guido
-Kanschat in 2003. It has undergone significant revisions by Bärbel
-Janssen, Guido Kanschat and Wolfgang Bangerth in 2009 and 2010 to demonstrate
-multigrid algorithms on adaptively refined meshes.
-</i>
-
-
-<a name="Intro"></a>
-<h1>Introduction</h1>
-
-
-This example shows the basic usage of the multilevel functions in
-deal.II. It solves the same problem as used in step-6,
-but demonstrating the things one has to provide when using multigrid
-as a preconditioner. In particular, this requires that we define a
-hierarchy of levels, provide transfer operators from one level to the
-next and back, and provide representations of the Laplace operator on
-each level.
-
-In order to allow sufficient flexibility in conjunction with systems of
-differential equations and block preconditioners, quite a few different objects
-have to be created before starting the multilevel method, although
-most of what needs to be done is provided by deal.II itself. These are
-<ul>
-<li>An the object handling transfer between grids; we use the
- MGTransferPrebuilt class for this that does almost all of the work
- inside the library.
-<li>The solver on the coarsest level; here, we use MGCoarseGridHouseholder.
-<li>The smoother on all other levels, which in our case will be the
- MGSmootherRelaxation class using SOR as the underlying method
-<li>And mg::Matrix, a class having a special level multiplication, i.e. we
- basically store one matrix per grid level and allow multiplication
- with it.
-</ul>
-Most of these objects will only be needed inside the function that
-actually solves the linear system. There, these objects are combined
-in an object of type Multigrid, containing the implementation of the
-V-cycle, which is in turn used by the preconditioner PreconditionMG,
-ready for plug-in into a linear solver of the LAC library.
-
-The multilevel method in deal.II follows in many respects the outlines
-of the various publications by James Bramble, Joseph Pasciak and
-Jinchao Xu (i.e. the "BPX" framework). In order to understand many of
-the options, a rough familiarity with their work is quite helpful.
-
-However, in comparison to this framework, the implementation in
-deal.II has to take into account the fact that we want to solve linear
-systems on adaptively refined meshes. This leads to the complication
-that it isn't quite as clear any more what exactly a "level" in a
-multilevel hierarchy of a mesh is. The following image shows what we
-consider to be a "level":
-
-<p align="center">
- @image html "hanging_nodes.png" ""
-</p>
-
-In other words, the fine level in this mesh consists only of the
-degrees of freedom that are defined on the refined cells, but does not
-extend to that part of the domain that is not refined. While this
-guarantees that the overall effort grows as ${\cal O}(N)$ as necessary
-for optimal multigrid complexity, it leads to problems when defining
-where to smooth and what boundary conditions to pose for the operators
-defined on individual levels if the level boundary is not an external
-boundary. These questions are discussed in detail in the
-@ref mg_paper "Multigrid paper by Janssen and Kanschat" that describes
-the implementation in deal.II.
-
-
-
-<h3>The testcase</h3>
-
-The problem we solve here is exactly the same as in
-step-6, the only difference being the solver we use
-here. You may want to look there for a definition of what we solve,
-right hand side and boundary conditions. Obviously, the program would
-also work if we changed the geometry and other pieces of data that
-defines this particular problem.
-
-The things that are new are all those parts that concern the
-multigrid. In particular, this includes the following members of the
-main class:
-- <code>LaplaceProblem::mg_dof_handler</code>
-- <code>LaplaceProblem::mg_sparsity</code>
-- <code>LaplaceProblem::mg_matrices</code>
-- <code>LaplaceProblem::mg_interface_matrices_up</code>
-- <code>LaplaceProblem::assemble_multigrid ()</code>
-- <code>LaplaceProblem::solve ()</code>
-Take a look at these functions.
+++ /dev/null
-unfinished
+++ /dev/null
-<h1>Results</h1>
-
-The output that this program generates is, of course, the same as that
-of step-6, so you may see there for more results. On the
-other hand, since no tutorial program is a good one unless it has at
-least one colorful picture, here is, again, the solution:
-
-
-When run, the output of this program is
-<pre>
-Cycle 0:
- Number of active cells: 20
- Number of degrees of freedom: 25 (by level: 8, 25)
- 7 CG iterations needed to obtain convergence.
-Cycle 1:
- Number of active cells: 44
- Number of degrees of freedom: 57 (by level: 8, 25, 48)
- 8 CG iterations needed to obtain convergence.
-Cycle 2:
- Number of active cells: 92
- Number of degrees of freedom: 117 (by level: 8, 25, 80, 60)
- 9 CG iterations needed to obtain convergence.
-Cycle 3:
- Number of active cells: 188
- Number of degrees of freedom: 221 (by level: 8, 25, 80, 200)
- 12 CG iterations needed to obtain convergence.
-Cycle 4:
- Number of active cells: 416
- Number of degrees of freedom: 485 (by level: 8, 25, 89, 288, 280)
- 13 CG iterations needed to obtain convergence.
-Cycle 5:
- Number of active cells: 800
- Number of degrees of freedom: 925 (by level: 8, 25, 89, 288, 784, 132)
- 14 CG iterations needed to obtain convergence.
-Cycle 6:
- Number of active cells: 1628
- Number of degrees of freedom: 1865 (by level: 8, 25, 89, 304, 1000, 1164, 72)
- 14 CG iterations needed to obtain convergence.
-Cycle 7:
- Number of active cells: 3194
- Number of degrees of freedom: 3603 (by level: 8, 25, 89, 328, 1032, 2200, 1392)
- 16 CG iterations needed to obtain convergence.
-</pre>
-That's not perfect — we would have hoped for a constant number
-of iterations rather than one that increases as we get more and more
-degrees of freedom — but it is also not far away. The reason for
-this is easy enough to understand, however: since we have a strongly
-varying coefficient, the operators that we assembly by quadrature on
-the lower levels become worse and worse approximations of the operator
-on the finest level. Consequently, even if we had perfect solvers on
-the coarser levels, they would not be good preconditioners on the
-finest level. This theory is easily tested by comparing results when
-we use a constant coefficient: in that case, the number of iterations
-remains constant at 9 after the first three or four refinement steps.
-
-We can also compare what this program produces with how @ref step_5
-"step-5" performed. To solve the same problem as in step-5, the only
-two changes that are necessary are (i) to replace the body of the
-function <code>LaplaceProblem::refine_grid</code> by a call to
-<code>triangulation.refine_global(1)</code>, and (ii) to use the same
-SolverControl object and tolerance as in step-5 — the rest of the
-program remains unchanged. In that case, here is how the solvers used
-in step-5 and the multigrid solver used in the current program
-compare:
-<table align="center">
-<tr><th>cells</th><th>step-5</th><th>step-16</th></tr>
-<tr><td>20</td> <td>13</td> <td>6</td> </tr>
-<tr><td>80</td> <td>17</td> <td>7</td> </tr>
-<tr><td>320</td> <td>29</td> <td>9</td> </tr>
-<tr><td>1280</td> <td>51</td> <td>10</td> </tr>
-<tr><td>5120</td> <td>94</td> <td>11</td> </tr>
-<tr><td>20480</td><td>180</td><td>13</td></tr>
-</table>
-This isn't only fewer iterations than in step-5 (each of which
-is, however, much more expensive) but more importantly, the number of
-iterations also grows much more slowly under mesh refinement (again,
-it would be almost constant if the coefficient was constant rather
-than strongly varying as chosen here). This justifies the common
-observation that, whenever possible, multigrid methods should be used
-for second order problems.
-
-
-<h3> Possible extensions </h3>
-
-A close inspection of this program's performance shows that it is mostly
-dominated by matrix-vector operations. step-37 shows one way
-how this can be avoided by working with matrix-free methods.
-
-Another avenue would be to use algebraic multigrid methods. The
-geometric multigrid method used here can at times be a bit awkward to
-implement because it needs all those additional data structures, and
-it becomes even more difficult if the program is to run in %parallel on
-machines coupled through MPI, for example. In that case, it would be
-simpler if one could use a black-box preconditioner that uses some
-sort of multigrid hierarchy for good performance but can figure out
-level matrices and similar things out by itself. Algebraic multigrid
-methods do exactly this, and we will use them in
-step-31 for the solution of a Stokes problem.
+++ /dev/null
-Multigrid on adaptive meshes.
+++ /dev/null
-/* ---------------------------------------------------------------------
- *
- * Copyright (C) 2003 - 2017 by the deal.II authors
- *
- * This file is part of the deal.II library.
- *
- * The deal.II library is free software; you can use it, redistribute
- * it, and/or modify it under the terms of the GNU Lesser General
- * Public License as published by the Free Software Foundation; either
- * version 2.1 of the License, or (at your option) any later version.
- * The full text of the license can be found in the file LICENSE at
- * the top level of the deal.II distribution.
- *
- * ---------------------------------------------------------------------
-
- *
- * Author: Guido Kanschat and Timo Heister
- */
-
-
-// @note: This a work in progress example of parallel geometric
-// multigrid. Some parts are still in heavy development.
-
-// This program is a parallel version of step-16 with a slightly different
-// problem setup.
-
-// @sect3{Include files}
-
-// Again, the first few include files
-// are already known, so we won't
-// comment on them:
-#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/function.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/utilities.h>
-#include <deal.II/base/conditional_ostream.h>
-
-#include <deal.II/lac/constraint_matrix.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/sparse_matrix.h>
-#include <deal.II/lac/solver_cg.h>
-#include <deal.II/lac/solver_gmres.h>
-#include <deal.II/lac/precondition.h>
-
-#include <deal.II/grid/tria.h>
-#include <deal.II/grid/tria_accessor.h>
-#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/grid_out.h>
-#include <deal.II/grid/grid_refinement.h>
-#include <deal.II/grid/tria_boundary_lib.h>
-
-#include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/dofs/dof_tools.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/numerics/data_out.h>
-#include <deal.II/numerics/error_estimator.h>
-
-#include <deal.II/base/index_set.h>
-#include <deal.II/distributed/tria.h>
-#include <deal.II/distributed/grid_refinement.h>
-
-#include <deal.II/multigrid/mg_constrained_dofs.h>
-#include <deal.II/multigrid/multigrid.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/lac/generic_linear_algebra.h>
-
-// #define USE_PETSC_LA PETSc is not quite supported yet
-
-namespace LA
-{
-#ifdef USE_PETSC_LA
- using namespace dealii::LinearAlgebraPETSc;
-#else
- using namespace dealii::LinearAlgebraTrilinos;
-#endif
-}
-
-// This is C++:
-#include <iostream>
-#include <fstream>
-#include <sstream>
-
-// The last step is as in all
-// previous programs:
-namespace Step50
-{
- using namespace dealii;
-
-
- // @sect3{The <code>LaplaceProblem</code> class template}
-
- // This main class is very similar to step-16, except that we are storing a
- // parallel Triangulation and parallel versions of matrices and vectors.
- template <int dim>
- class LaplaceProblem
- {
- public:
- LaplaceProblem (const unsigned int deg);
- void run ();
-
- private:
- void setup_system ();
- void assemble_system ();
- void assemble_multigrid ();
- void solve ();
- void refine_grid ();
- void output_results (const unsigned int cycle) const;
-
- ConditionalOStream pcout;
-
- parallel::distributed::Triangulation<dim> triangulation;
- FE_Q<dim> fe;
- DoFHandler<dim> mg_dof_handler;
-
- typedef LA::MPI::SparseMatrix matrix_t;
- typedef LA::MPI::Vector vector_t;
-
- matrix_t system_matrix;
-
- IndexSet locally_relevant_set;
-
- ConstraintMatrix constraints;
-
- vector_t solution;
- vector_t system_rhs;
-
- const unsigned int degree;
-
- // Finally we are storing the various parallel multigrid matrices. Our
- // problem is self-adjoint, so the interface matrices are the transpose
- // of each other, so we only need to compute/store them once.
- MGLevelObject<matrix_t> mg_matrices;
- MGLevelObject<matrix_t> mg_interface_matrices;
- //
- MGConstrainedDoFs mg_constrained_dofs;
- };
-
-
-
- // @sect3{Nonconstant coefficients}
-
- // The implementation of nonconstant
- // coefficients is copied verbatim
- // from step-5 and step-6:
-
- 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;
-
- virtual void value_list (const std::vector<Point<dim> > &points,
- std::vector<double> &values,
- const unsigned int component = 0) const;
- };
-
-
-
- template <int dim>
- double Coefficient<dim>::value (const Point<dim> &p,
- const unsigned int) const
- {
- if (p.square() < 0.5*0.5)
- return 5;
- else
- return 1;
- }
-
-
-
- template <int dim>
- void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points,
- std::vector<double> &values,
- const unsigned int component) const
- {
- (void)component;
- const unsigned int n_points = points.size();
-
- Assert (values.size() == n_points,
- ExcDimensionMismatch (values.size(), n_points));
-
- Assert (component == 0,
- ExcIndexRange (component, 0, 1));
-
- for (unsigned int i=0; i<n_points; ++i)
- values[i] = Coefficient<dim>::value (points[i]);
- }
-
-
- // @sect3{The <code>LaplaceProblem</code> class implementation}
-
- // @sect4{LaplaceProblem::LaplaceProblem}
-
- // The constructor is left mostly
- // unchanged. We take the polynomial degree
- // of the finite elements to be used as a
- // constructor argument and store it in a
- // member variable.
- //
- // By convention, all adaptively refined
- // triangulations in deal.II never change by
- // more than one level across a face between
- // cells. For our multigrid algorithms,
- // however, we need a slightly stricter
- // guarantee, namely that the mesh also does
- // not change by more than refinement level
- // across vertices that might connect two
- // cells. In other words, we must prevent the
- // following situation:
- //
- // @image html limit_level_difference_at_vertices.png ""
- //
- // This is achieved by passing the
- // Triangulation::limit_level_difference_at_vertices
- // flag to the constructor of the
- // triangulation class.
- template <int dim>
- LaplaceProblem<dim>::LaplaceProblem (const unsigned int degree)
- :
- pcout (std::cout,
- (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)
- == 0)),
- triangulation (MPI_COMM_WORLD,Triangulation<dim>::
- limit_level_difference_at_vertices,
- parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy),
- fe (degree),
- mg_dof_handler (triangulation),
- degree(degree)
- {}
-
-
- // @sect4{LaplaceProblem::setup_system}
-
- // The following function extends what the
- // corresponding one in step-6 did. The top
- // part, apart from the additional output,
- // does the same:
- template <int dim>
- void LaplaceProblem<dim>::setup_system ()
- {
- mg_dof_handler.distribute_dofs (fe);
- mg_dof_handler.distribute_mg_dofs (fe);
-
- DoFTools::extract_locally_relevant_dofs (mg_dof_handler,
- locally_relevant_set);
-
- solution.reinit(mg_dof_handler.locally_owned_dofs(), MPI_COMM_WORLD);
- system_rhs.reinit(mg_dof_handler.locally_owned_dofs(), MPI_COMM_WORLD);
-
- // But it starts to be a wee bit different
- // here, although this still doesn't have
- // anything to do with multigrid
- // methods. step-6 took care of boundary
- // values and hanging nodes in a separate
- // step after assembling the global matrix
- // from local contributions. This works,
- // but the same can be done in a slightly
- // simpler way if we already take care of
- // these constraints at the time of copying
- // local contributions into the global
- // matrix. To this end, we here do not just
- // compute the constraints do to hanging
- // nodes, but also due to zero boundary
- // conditions. We will
- // use this set of constraints later on to
- // help us copy local contributions
- // correctly into the global linear system
- // right away, without the need for a later
- // clean-up stage:
- constraints.reinit (locally_relevant_set);
- DoFTools::make_hanging_node_constraints (mg_dof_handler, constraints);
-
- std::set<types::boundary_id> dirichlet_boundary_ids;
- typename FunctionMap<dim>::type dirichlet_boundary;
- ConstantFunction<dim> homogeneous_dirichlet_bc (1.0);
- dirichlet_boundary_ids.insert(0);
- dirichlet_boundary[0] = &homogeneous_dirichlet_bc;
- VectorTools::interpolate_boundary_values (mg_dof_handler,
- dirichlet_boundary,
- constraints);
- constraints.close ();
-
- DynamicSparsityPattern dsp(mg_dof_handler.n_dofs(), mg_dof_handler.n_dofs());
- DoFTools::make_sparsity_pattern (mg_dof_handler, dsp, constraints);
- system_matrix.reinit (mg_dof_handler.locally_owned_dofs(), dsp, MPI_COMM_WORLD, true);
-
-
- // The multigrid constraints have to be
- // initialized. They need to know about
- // the boundary values as well, so we
- // pass the <code>dirichlet_boundary</code>
- // here as well.
- mg_constrained_dofs.clear();
- mg_constrained_dofs.initialize(mg_dof_handler);
- mg_constrained_dofs.make_zero_boundary_constraints(mg_dof_handler, dirichlet_boundary_ids);
-
-
- // Now for the things that concern the
- // multigrid data structures. First, we
- // resize the multilevel objects to hold
- // matrices and sparsity patterns for every
- // level. The coarse level is zero (this is
- // mandatory right now but may change in a
- // future revision). Note that these
- // functions take a complete, inclusive
- // range here (not a starting index and
- // size), so the finest level is
- // <code>n_levels-1</code>. We first have
- // to resize the container holding the
- // SparseMatrix classes, since they have to
- // release their SparsityPattern before the
- // can be destroyed upon resizing.
- const unsigned int n_levels = triangulation.n_global_levels();
-
- mg_interface_matrices.resize(0, n_levels-1);
- mg_interface_matrices.clear_elements ();
- mg_matrices.resize(0, n_levels-1);
- mg_matrices.clear_elements ();
-
- // Now, we have to provide a matrix on each
- // level. To this end, we first use the
- // MGTools::make_sparsity_pattern function
- // to first generate a preliminary
- // compressed sparsity pattern on each
- // level (see the @ref Sparsity module for
- // more information on this topic) and then
- // copy it over to the one we really
- // want. The next step is to initialize
- // both kinds of level matrices with these
- // sparsity patterns.
- //
- // It may be worth pointing out that the
- // interface matrices only have entries for
- // degrees of freedom that sit at or next
- // to the interface between coarser and
- // finer levels of the mesh. They are
- // therefore even sparser than the matrices
- // on the individual levels of our
- // multigrid hierarchy. If we were more
- // concerned about memory usage (and
- // possibly the speed with which we can
- // multiply with these matrices), we should
- // use separate and different sparsity
- // patterns for these two kinds of
- // matrices.
- for (unsigned int level=0; level<n_levels; ++level)
- {
- DynamicSparsityPattern dsp(mg_dof_handler.n_dofs(level),
- mg_dof_handler.n_dofs(level));
- MGTools::make_sparsity_pattern(mg_dof_handler, dsp, level);
-
- mg_matrices[level].reinit(mg_dof_handler.locally_owned_mg_dofs(level),
- mg_dof_handler.locally_owned_mg_dofs(level),
- dsp,
- MPI_COMM_WORLD, true);
-
- mg_interface_matrices[level].reinit(mg_dof_handler.locally_owned_mg_dofs(level),
- mg_dof_handler.locally_owned_mg_dofs(level),
- dsp,
- MPI_COMM_WORLD, true);
- }
- }
-
-
- // @sect4{LaplaceProblem::assemble_system}
-
- // The following function assembles the
- // linear system on the finest level of the
- // mesh. It is almost exactly the same as in
- // step-6, with the exception that we don't
- // eliminate hanging nodes and boundary
- // values after assembling, but while copying
- // local contributions into the global
- // matrix. This is not only simpler but also
- // more efficient for large problems.
- //
- // This latter trick is something that only
- // found its way into deal.II over time and
- // wasn't used in the initial version of this
- // tutorial program. There is, however, a
- // discussion of this function in the
- // introduction of step-27.
- template <int dim>
- void LaplaceProblem<dim>::assemble_system ()
- {
- const QGauss<dim> quadrature_formula(degree+1);
-
- FEValues<dim> fe_values (fe, quadrature_formula,
- update_values | update_gradients |
- 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();
-
- FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
- Vector<double> cell_rhs (dofs_per_cell);
-
- std::vector<types::global_dof_index> 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)
- if (cell->is_locally_owned())
- {
- cell_matrix = 0;
- cell_rhs = 0;
-
- fe_values.reinit (cell);
-
- coefficient.value_list (fe_values.get_quadrature_points(),
- coefficient_values);
-
- for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- cell_matrix(i,j) += (coefficient_values[q_point] *
- fe_values.shape_grad(i,q_point) *
- fe_values.shape_grad(j,q_point) *
- fe_values.JxW(q_point));
-
- cell_rhs(i) += (fe_values.shape_value(i,q_point) *
- 10.0 *
- fe_values.JxW(q_point));
- }
-
- cell->get_dof_indices (local_dof_indices);
- constraints.distribute_local_to_global (cell_matrix, cell_rhs,
- local_dof_indices,
- system_matrix, system_rhs);
- }
-
- system_matrix.compress(VectorOperation::add);
- system_rhs.compress(VectorOperation::add);
- }
-
-
- // @sect4{LaplaceProblem::assemble_multigrid}
-
- // The next function is the one that builds
- // the linear operators (matrices) that
- // define the multigrid method on each level
- // of the mesh. The integration core is the
- // same as above, but the loop below will go
- // over all existing cells instead of just
- // the active ones, and the results must be
- // entered into the correct matrix. Note also
- // that since we only do multilevel
- // preconditioning, no right-hand side needs
- // to be assembled here.
- //
- // Before we go there, however, we have to
- // take care of a significant amount of book
- // keeping:
- template <int dim>
- void LaplaceProblem<dim>::assemble_multigrid ()
- {
- QGauss<dim> quadrature_formula(1+degree);
-
- FEValues<dim> fe_values (fe, quadrature_formula,
- update_values | update_gradients |
- 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();
-
- FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
-
- std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
-
- const Coefficient<dim> coefficient;
- std::vector<double> coefficient_values (n_q_points);
-
- // Next a few things that are specific to building the multigrid
- // data structures (since we only need them in the current
- // function, rather than also elsewhere, we build them here
- // instead of the <code>setup_system</code> function). Some of the
- // following may be a bit obscure if you're not familiar with the
- // algorithm actually implemented in deal.II to support multilevel
- // algorithms on adaptive meshes; if some of the things below seem
- // strange, take a look at the @ref mg_paper.
- //
- // Our first job is to identify those degrees of freedom on each level
- // that are located on interfaces between adaptively refined levels, and
- // those that lie on the interface but also on the exterior boundary of
- // the domain. The <code>MGConstrainedDoFs</code> already computed the
- // information for us when we called initialize in
-
- // <code>setup_system()</code>.
- // of type IndexSet on each level (get_refinement_edge_indices(),
-
- // The indices just identified will later be used to decide where
- // the assembled value has to be added into on each level. On the
- // other hand, we also have to impose zero boundary conditions on
- // the external boundary of each level. But this the
- // <code>MGConstrainedDoFs</code> knows it. So we simply ask for them
- // by calling <code>get_boundary_indices ()</code>. The third
- // step is to construct constraints on all those degrees of
- // freedom: their value should be zero after each application of
- // the level operators. To this end, we construct ConstraintMatrix
- // objects for each level, and add to each of these constraints
- // for each degree of freedom. Due to the way the ConstraintMatrix
- // stores its data, the function to add a constraint on a single
- // degree of freedom and force it to be zero is called
- // ConstraintMatrix::add_line(); doing so for several degrees of
- // freedom at once can be done using
- // ConstraintMatrix::add_lines():
- std::vector<ConstraintMatrix> boundary_constraints (triangulation.n_global_levels());
- ConstraintMatrix empty_constraints;
- for (unsigned int level=0; level<triangulation.n_global_levels(); ++level)
- {
- IndexSet dofset;
- DoFTools::extract_locally_relevant_level_dofs (mg_dof_handler, level, dofset);
- boundary_constraints[level].reinit(dofset);
- boundary_constraints[level].add_lines (mg_constrained_dofs.get_refinement_edge_indices(level));
- boundary_constraints[level].add_lines (mg_constrained_dofs.get_boundary_indices(level));
-
- boundary_constraints[level].close ();
- }
-
- // Now that we're done with most of our preliminaries, let's start
- // the integration loop. It looks mostly like the loop in
- // <code>assemble_system</code>, with two exceptions: (i) we don't
- // need a right hand side, and more significantly (ii) we don't
- // just loop over all active cells, but in fact all cells, active
- // or not. Consequently, the correct iterator to use is
- // DoFHandler::cell_iterator rather than
- // DoFHandler::active_cell_iterator. Let's go about it:
- typename DoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
- endc = mg_dof_handler.end();
-
- for (; cell!=endc; ++cell)
- if (cell->level_subdomain_id()==triangulation.locally_owned_subdomain())
- {
- cell_matrix = 0;
- fe_values.reinit (cell);
-
- coefficient.value_list (fe_values.get_quadrature_points(),
- coefficient_values);
-
- for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- cell_matrix(i,j) += (coefficient_values[q_point] *
- fe_values.shape_grad(i,q_point) *
- fe_values.shape_grad(j,q_point) *
- fe_values.JxW(q_point));
-
- // The rest of the assembly is again slightly
- // different. This starts with a gotcha that is easily
- // forgotten: The indices of global degrees of freedom we
- // want here are the ones for current level, not for the
- // global matrix. We therefore need the function
- // MGDoFAccessorLLget_mg_dof_indices, not
- // MGDoFAccessor::get_dof_indices as used in the assembly of
- // the global system:
- cell->get_mg_dof_indices (local_dof_indices);
-
- // Next, we need to copy local contributions into the level
- // objects. We can do this in the same way as in the global
- // assembly, using a constraint object that takes care of
- // constrained degrees (which here are only boundary nodes,
- // as the individual levels have no hanging node
- // constraints). Note that the
- // <code>boundary_constraints</code> object makes sure that
- // the level matrices contains no contributions from degrees
- // of freedom at the interface between cells of different
- // refinement level.
- boundary_constraints[cell->level()]
- .distribute_local_to_global (cell_matrix,
- local_dof_indices,
- mg_matrices[cell->level()]);
-
- // The next step is again slightly more obscure (but
- // explained in the @ref mg_paper): We need the remainder of
- // the operator that we just copied into the
- // <code>mg_matrices</code> object, namely the part on the
- // interface between cells at the current level and cells
- // one level coarser. This matrix exists in two directions:
- // for interior DoFs (index $i$) of the current level to
- // those sitting on the interface (index $j$), and the other
- // way around. Of course, since we have a symmetric
- // operator, one of these matrices is the transpose of the
- // other.
- //
- // The way we assemble these matrices is as follows: since
- // the are formed from parts of the local contributions, we
- // first delete all those parts of the local contributions
- // that we are not interested in, namely all those elements
- // of the local matrix for which not $i$ is an interface DoF
- // and $j$ is not. The result is one of the two matrices
- // that we are interested in, and we then copy it into the
- // <code>mg_interface_matrices</code> object. The
- // <code>boundary_interface_constraints</code> object at the
- // same time makes sure that we delete contributions from
- // all degrees of freedom that are not only on the interface
- // but also on the external boundary of the domain.
- //
- // The last part to remember is how to get the other
- // matrix. Since it is only the transpose, we will later (in
- // the <code>solve()</code> function) be able to just pass
- // the transpose matrix where necessary.
-
- const IndexSet &interface_dofs_on_level
- = mg_constrained_dofs.get_refinement_edge_indices(cell->level());
- const unsigned int lvl = cell->level();
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- if (interface_dofs_on_level.is_element(local_dof_indices[i]) // at_refinement_edge(i)
- &&
- !interface_dofs_on_level.is_element(local_dof_indices[j]) // !at_refinement_edge(j)
- &&
- (
- (!mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[i])
- &&
- !mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[j])
- ) // ( !boundary(i) && !boundary(j) )
- ||
- (
- mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[i])
- &&
- local_dof_indices[i]==local_dof_indices[j]
- ) // ( boundary(i) && boundary(j) && i==j )
- )
- )
- {
- // do nothing, so add entries to interface matrix
- }
- else
- {
- cell_matrix(i,j) = 0;
- }
-
-
- empty_constraints
- .distribute_local_to_global (cell_matrix,
- local_dof_indices,
- mg_interface_matrices[cell->level()]);
- }
-
- for (unsigned int i=0; i<triangulation.n_global_levels(); ++i)
- {
- mg_matrices[i].compress(VectorOperation::add);
- mg_interface_matrices[i].compress(VectorOperation::add);
- }
- }
-
-
-
- // @sect4{LaplaceProblem::solve}
-
- // This is the other function that is significantly different in
- // support of the multigrid solver (or, in fact, the preconditioner
- // for which we use the multigrid method).
- //
- // Let us start out by setting up two of the components of
- // multilevel methods: transfer operators between levels, and a
- // solver on the coarsest level. In finite element methods, the
- // transfer operators are derived from the finite element function
- // spaces involved and can often be computed in a generic way
- // independent of the problem under consideration. In that case, we
- // can use the MGTransferPrebuilt class that, given the constraints
- // on the global level and an DoFHandler object computes the
- // matrices corresponding to these transfer operators.
- //
- // The second part of the following lines deals with the coarse grid
- // solver. Since our coarse grid is very coarse indeed, we decide
- // for a direct solver (a Householder decomposition of the coarsest
- // level matrix), even if its implementation is not particularly
- // sophisticated. If our coarse mesh had many more cells than the
- // five we have here, something better suited would obviously be
- // necessary here.
- template <int dim>
- void LaplaceProblem<dim>::solve ()
- {
- // Create the object that deals with the transfer between
- // different refinement levels.
- MGTransferPrebuilt<vector_t> mg_transfer(mg_constrained_dofs);
- // Now the prolongation matrix has to be built.
- mg_transfer.build_matrices(mg_dof_handler);
-
- matrix_t &coarse_matrix = mg_matrices[0];
-
- SolverControl coarse_solver_control (1000, 1e-10, false, false);
- SolverCG<vector_t> coarse_solver(coarse_solver_control);
- PreconditionIdentity id;
- MGCoarseGridLACIteration<SolverCG<vector_t>,vector_t> coarse_grid_solver(coarse_solver,
- coarse_matrix,
- id);
-
- // The next component of a multilevel solver or preconditioner is
- // that we need a smoother on each level. A common choice for this
- // is to use the application of a relaxation method (such as the
- // SOR, Jacobi or Richardson method). The MGSmootherPrecondition
- // class provides support for this kind of smoother. Here, we opt
- // for the application of a single SOR iteration. To this end, we
- // define an appropriate <code>typedef</code> and then setup a
- // smoother object.
- //
- // The last step is to initialize the smoother object with our
- // level matrices and to set some smoothing parameters. The
- // <code>initialize()</code> function can optionally take
- // additional arguments that will be passed to the smoother object
- // on each level. In the current case for the SOR smoother, this
- // could, for example, include a relaxation parameter. However, we
- // here leave these at their default values. The call to
- // <code>set_steps()</code> indicates that we will use two pre-
- // and two post-smoothing steps on each level; to use a variable
- // number of smoother steps on different levels, more options can
- // be set in the constructor call to the <code>mg_smoother</code>
- // object.
- //
- // The last step results from the fact that
- // we use the SOR method as a smoother -
- // which is not symmetric - but we use the
- // conjugate gradient iteration (which
- // requires a symmetric preconditioner)
- // below, we need to let the multilevel
- // preconditioner make sure that we get a
- // symmetric operator even for nonsymmetric
- // smoothers:
- typedef LA::MPI::PreconditionJacobi Smoother;
- MGSmootherPrecondition<matrix_t, Smoother, vector_t> mg_smoother;
- mg_smoother.initialize(mg_matrices, Smoother::AdditionalData(0.5));
- mg_smoother.set_steps(2);
- //mg_smoother.set_symmetric(false);
-
- // The next preparatory step is that we
- // must wrap our level and interface
- // matrices in an object having the
- // required multiplication functions. We
- // will create two objects for the
- // interface objects going from coarse to
- // fine and the other way around; the
- // multigrid algorithm will later use the
- // transpose operator for the latter
- // operation, allowing us to initialize
- // both up and down versions of the
- // operator with the matrices we already
- // built:
- mg::Matrix<vector_t> mg_matrix(mg_matrices);
- mg::Matrix<vector_t> mg_interface_up(mg_interface_matrices);
- mg::Matrix<vector_t> mg_interface_down(mg_interface_matrices);
-
- // Now, we are ready to set up the
- // V-cycle operator and the
- // multilevel preconditioner.
- Multigrid<vector_t > mg(mg_dof_handler,
- mg_matrix,
- coarse_grid_solver,
- mg_transfer,
- mg_smoother,
- mg_smoother);
- //mg.set_debug(6);
- mg.set_edge_matrices(mg_interface_down, mg_interface_up);
-
- PreconditionMG<dim, vector_t, MGTransferPrebuilt<vector_t> >
- preconditioner(mg_dof_handler, mg, mg_transfer);
-
-
- // With all this together, we can finally
- // get about solving the linear system in
- // the usual way:
- SolverControl solver_control (500, 1e-8*system_rhs.l2_norm(), false);
- SolverCG<vector_t> solver (solver_control);
-
- if (false)
- {
- /*
- // code to optionally compare to Trilinos ML
- TrilinosWrappers::PreconditionAMG prec;
-
- TrilinosWrappers::PreconditionAMG::AdditionalData Amg_data;
- // Amg_data.constant_modes = constant_modes;
- Amg_data.elliptic = true;
- Amg_data.higher_order_elements = true;
- Amg_data.smoother_sweeps = 2;
- Amg_data.aggregation_threshold = 0.02;
- // Amg_data.symmetric = true;
-
- prec.initialize (system_matrix,
- Amg_data);
- solver.solve (system_matrix, solution, system_rhs, prec);
- */
- }
- else
- {
- solver.solve (system_matrix, solution, system_rhs,
- preconditioner);
- }
- pcout << " CG converged in " << solver_control.last_step() << " iterations." << std::endl;
-
- constraints.distribute (solution);
- }
-
-
-
- // @sect4{Postprocessing}
-
- // The following two functions postprocess a solution once it is
- // computed. In particular, the first one refines the mesh at the beginning
- // of each cycle while the second one outputs results at the end of each
- // such cycle. The <code>refine_grid()</code> method is almost unchanged
- // from step-6: the only substantial difference is that this method uses a
- // distributed grid refinement function instead of a serial one. The
- // <code>output_results()</code> method is quite different since each
- // processor writes only part of the overall graphical output.
- template <int dim>
- void LaplaceProblem<dim>::refine_grid ()
- {
- Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
-
- LA::MPI::Vector temp_solution;
- temp_solution.reinit(locally_relevant_set, MPI_COMM_WORLD);
- temp_solution = solution;
-
- KellyErrorEstimator<dim>::estimate (static_cast<DoFHandler<dim>&>(mg_dof_handler),
- QGauss<dim-1>(degree+1),
- typename FunctionMap<dim>::type(),
- temp_solution,
- estimated_error_per_cell);
-
- parallel::distributed::GridRefinement::
- refine_and_coarsen_fixed_fraction (triangulation,
- estimated_error_per_cell,
- 0.3, 0.0);
-
- triangulation.execute_coarsening_and_refinement ();
- }
-
-
-
- template <int dim>
- void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
- {
- DataOut<dim> data_out;
-
- LA::MPI::Vector temp_solution;
- temp_solution.reinit(locally_relevant_set, MPI_COMM_WORLD);
- temp_solution = solution;
-
-
- LA::MPI::Vector temp = solution;
- system_matrix.residual(temp,solution,system_rhs);
- LA::MPI::Vector res_ghosted = temp_solution;
- res_ghosted = temp;
-
- data_out.attach_dof_handler (mg_dof_handler);
- data_out.add_data_vector (temp_solution, "solution");
- data_out.add_data_vector (res_ghosted, "res");
- Vector<float> subdomain (triangulation.n_active_cells());
- for (unsigned int i=0; i<subdomain.size(); ++i)
- subdomain(i) = triangulation.locally_owned_subdomain();
- data_out.add_data_vector (subdomain, "subdomain");
-
- data_out.build_patches (0);
-
- const std::string filename = ("solution-" +
- Utilities::int_to_string (cycle, 5) +
- "." +
- Utilities::int_to_string
- (triangulation.locally_owned_subdomain(), 4) +
- ".vtu");
- std::ofstream output (filename.c_str());
- data_out.write_vtu (output);
-
- if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
- {
- std::vector<std::string> filenames;
- for (unsigned int i=0; i<Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); ++i)
- filenames.push_back (std::string("solution-") +
- Utilities::int_to_string (cycle, 5) +
- "." +
- Utilities::int_to_string(i, 4) +
- ".vtu");
- const std::string
- pvtu_master_filename = ("solution-" +
- Utilities::int_to_string (cycle, 5) +
- ".pvtu");
- std::ofstream pvtu_master (pvtu_master_filename.c_str());
- data_out.write_pvtu_record (pvtu_master, filenames);
-
- const std::string
- visit_master_filename = ("solution-" +
- Utilities::int_to_string (cycle, 5) +
- ".visit");
- std::ofstream visit_master (visit_master_filename.c_str());
- DataOutBase::write_visit_record (visit_master, filenames);
-
- std::cout << " wrote " << pvtu_master_filename << std::endl;
-
- }
- }
-
-
- // @sect4{LaplaceProblem::run}
-
- // Like several of the functions above, this
- // is almost exactly a copy of of the
- // corresponding function in step-6. The only
- // difference is the call to
- // <code>assemble_multigrid</code> that takes
- // care of forming the matrices on every
- // level that we need in the multigrid
- // method.
- template <int dim>
- void LaplaceProblem<dim>::run ()
- {
- for (unsigned int cycle=0; cycle<15; ++cycle)
- {
- pcout << "Cycle " << cycle << ':' << std::endl;
-
- if (cycle == 0)
- {
- GridGenerator::hyper_cube (triangulation);
-
- triangulation.refine_global (4);
- }
- else
- refine_grid ();
-
- pcout << " Number of active cells: "
- << triangulation.n_global_active_cells()
- << std::endl;
-
- setup_system ();
-
- pcout << " Number of degrees of freedom: "
- << mg_dof_handler.n_dofs()
- << " (by level: ";
- for (unsigned int level=0; level<triangulation.n_global_levels(); ++level)
- pcout << mg_dof_handler.n_dofs(level)
- << (level == triangulation.n_global_levels()-1
- ? ")" : ", ");
- pcout << std::endl;
-
- assemble_system ();
- assemble_multigrid ();
-
- solve ();
- output_results (cycle);
- }
- }
-}
-
-
-// @sect3{The main() function}
-//
-// This is again the same function as
-// in step-6:
-int main (int argc, char *argv[])
-{
- dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
-
- try
- {
- using namespace dealii;
- using namespace Step50;
-
- LaplaceProblem<2> laplace_problem(1/*degree*/);
- laplace_problem.run ();
- }
- 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;
- }
- catch (...)
- {
- std::cerr << std::endl << std::endl
- << "----------------------------------------------------"
- << std::endl;
- std::cerr << "Unknown exception!" << std::endl
- << "Aborting!" << std::endl
- << "----------------------------------------------------"
- << std::endl;
- throw;
- }
-
- return 0;
-}
projects / dealii.git / commitdiff