From: Timo Heister Date: Mon, 16 Dec 2013 20:41:13 +0000 (+0000) Subject: take over 32032. X-Git-Tag: v8.1.0~15 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=bfbfec1d285334f5cca8799c0b088b7551ca0136;p=dealii.git take over 32032. git-svn-id: https://svn.dealii.org/branches/releases/Branch-8-1@32033 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/developers/porting.html b/deal.II/doc/developers/porting.html index 27fbdd8bb9..b92ca9d935 100644 --- a/deal.II/doc/developers/porting.html +++ b/deal.II/doc/developers/porting.html @@ -47,7 +47,7 @@ DEAL_II_CXX_FLAGS_RELEASE - additional flags for the release library After that try to compile the library with minimal external dependencies (-DDEAL_II_ALLOW_AUTODETECTION=OFF, for - further information see the deal.II CMake + further information see the deal.II CMake documentation).
  • For adding permanent support for the unknown compiler to the diff --git a/deal.II/doc/developers/writing-documentation.html b/deal.II/doc/developers/writing-documentation.html index 249194b52c..5a878e14d2 100644 --- a/deal.II/doc/developers/writing-documentation.html +++ b/deal.II/doc/developers/writing-documentation.html @@ -39,7 +39,7 @@

    In addition to the API documentation, we maintain a series of - + well-documented example programs, which also follow a certain ``literate programming'' style in that the explanations of what is happening are integrated into the program source by means of diff --git a/deal.II/doc/users/cmakelists.html b/deal.II/doc/users/cmakelists.html index b1f13d8d7a..6e5f150f42 100644 --- a/deal.II/doc/users/cmakelists.html +++ b/deal.II/doc/users/cmakelists.html @@ -138,7 +138,7 @@ DEAL_II_SETUP_TARGET(mycode3) If the list gets longer, consider using a loop, possibly -with GLOB. +with GLOB.

    Adding libraries and common source files

    @@ -396,6 +396,7 @@ IF(NOT <expression>) ENDIF() +

    Loops are implemented with the help of WHILE and FOR statements. The former takes the same <expression> as the IF statement:

    @@ -588,7 +589,7 @@ INSTALL(TARGETS mycode DESTINATION bin) INSTALL(DIRECTORY run DESTINATION share/mycode/run) - +

    Autopilot style CMakeLists.txt

    diff --git a/deal.II/examples/step-47/CMakeLists.txt b/deal.II/examples/step-47/CMakeLists.txt deleted file mode 100644 index 0321ee6c15..0000000000 --- a/deal.II/examples/step-47/CMakeLists.txt +++ /dev/null @@ -1,31 +0,0 @@ -## -# 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: -SET(TARGET_SRC - ${TARGET}.cc - # You can specify additional files here! - ) - -# Usually, you will not need to modify anything beyond this point... - -CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8) - -FIND_PACKAGE(deal.II 8.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 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() diff --git a/deal.II/examples/step-47/doc/builds-on b/deal.II/examples/step-47/doc/builds-on deleted file mode 100644 index 9fdd726254..0000000000 --- a/deal.II/examples/step-47/doc/builds-on +++ /dev/null @@ -1 +0,0 @@ -step-5 diff --git a/deal.II/examples/step-47/doc/intro.dox b/deal.II/examples/step-47/doc/intro.dox deleted file mode 100644 index 1518606e43..0000000000 --- a/deal.II/examples/step-47/doc/intro.dox +++ /dev/null @@ -1,3 +0,0 @@ - -

    Introduction

    - diff --git a/deal.II/examples/step-47/doc/kind b/deal.II/examples/step-47/doc/kind deleted file mode 100644 index 15a13db451..0000000000 --- a/deal.II/examples/step-47/doc/kind +++ /dev/null @@ -1 +0,0 @@ -basic diff --git a/deal.II/examples/step-47/doc/points.dat b/deal.II/examples/step-47/doc/points.dat deleted file mode 100644 index 01ce16ad3c..0000000000 --- a/deal.II/examples/step-47/doc/points.dat +++ /dev/null @@ -1,21 +0,0 @@ -#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 diff --git a/deal.II/examples/step-47/doc/results.dox b/deal.II/examples/step-47/doc/results.dox deleted file mode 100644 index f4c6feefb5..0000000000 --- a/deal.II/examples/step-47/doc/results.dox +++ /dev/null @@ -1 +0,0 @@ -

    Results

    diff --git a/deal.II/examples/step-47/doc/tooltip b/deal.II/examples/step-47/doc/tooltip deleted file mode 100644 index 00262bfbfb..0000000000 --- a/deal.II/examples/step-47/doc/tooltip +++ /dev/null @@ -1 +0,0 @@ -Adaptive local refinement. Higher order elements. \ No newline at end of file diff --git a/deal.II/examples/step-47/doc/vertices.dat b/deal.II/examples/step-47/doc/vertices.dat deleted file mode 100644 index 8d1ef87323..0000000000 --- a/deal.II/examples/step-47/doc/vertices.dat +++ /dev/null @@ -1,20 +0,0 @@ -#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 diff --git a/deal.II/examples/step-47/step-47.cc b/deal.II/examples/step-47/step-47.cc deleted file mode 100644 index 87707ba908..0000000000 --- a/deal.II/examples/step-47/step-47.cc +++ /dev/null @@ -1,1137 +0,0 @@ -/* --------------------------------------------------------------------- - * $Id$ - * - * Copyright (C) 2011 - 2013 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 -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include - -#include -#include - - -#include - -#include - -#include - -namespace Step47 -{ - using namespace dealii; - - - - double sign (double d) - { - if (d > 0) - return 1; - else if (d < 0) - return -1; - else - return 0; - } - - - template - class LaplaceProblem - { - public: - LaplaceProblem (); - ~LaplaceProblem (); - - void run (); - - private: - bool interface_intersects_cell (const typename Triangulation::cell_iterator &cell) const; - std::pair > compute_quadrature(const Quadrature &plain_quadrature, const typename hp::DoFHandler::active_cell_iterator &cell, const std::vector &level_set_values); - void append_quadrature(const Quadrature &plain_quadrature, - const std::vector > &v , - std::vector > &xfem_points, - std::vector &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 triangulation; - - hp::DoFHandler dof_handler; - hp::FECollection fe_collection; - - ConstraintMatrix constraints; - - SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; - - Vector solution; - Vector system_rhs; - }; - - - - - template - class Coefficient : public Function - { - public: - Coefficient () : Function() {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; - - virtual void value_list (const std::vector > &points, - std::vector &values, - const unsigned int component = 0) const; - }; - - - - template - double Coefficient::value (const Point &p, - const unsigned int) const - { - if (p.square() < 0.5*0.5) - return 20; - else - return 1; - } - - - - template - void Coefficient::value_list (const std::vector > &points, - std::vector &values, - const unsigned int component) const - { - 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 - double exact_solution (const Point &p) - { - const double r = p.norm(); - - return (r < 0.5 - ? - 1./20 * (-1./4*r*r + 61./16) - : - 1./4 * (1-r*r)); - } - - - template - LaplaceProblem::LaplaceProblem () - : - dof_handler (triangulation) - { - fe_collection.push_back (FESystem (FE_Q(1), 1, - FE_Nothing(), 1)); - fe_collection.push_back (FESystem (FE_Q(1), 1, - FE_Q(1), 1)); - } - - - - template - LaplaceProblem::~LaplaceProblem () - { - dof_handler.clear (); - } - - - - template - double - level_set (const Point &p) - { - return p.norm() - 0.5; - } - - - - template - Tensor<1,dim> - grad_level_set (const Point &p) - { - return p / p.norm(); - } - - - - template - bool - LaplaceProblem:: - interface_intersects_cell (const typename Triangulation::cell_iterator &cell) const - { - for (unsigned int v=0; v::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 - void LaplaceProblem::setup_system () - { - for (typename hp::DoFHandler::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(); - - - CompressedSparsityPattern c_sparsity(dof_handler.n_dofs()); - DoFTools::make_sparsity_pattern (dof_handler, c_sparsity); - - constraints.condense (c_sparsity); - - sparsity_pattern.copy_from(c_sparsity); - - system_matrix.reinit (sparsity_pattern); - } - - - template - void LaplaceProblem::assemble_system () - { - const QGauss quadrature_formula(3); - - - FEValues plain_fe_values (fe_collection[0], quadrature_formula, - update_values | update_gradients | - update_quadrature_points | update_JxW_values); - FEValues enriched_fe_values (fe_collection[1], quadrature_formula, - update_values | update_gradients | - update_quadrature_points | update_JxW_values); - - const unsigned int n_q_points = quadrature_formula.size(); - - FullMatrix cell_matrix; - Vector cell_rhs; - - std::vector local_dof_indices; - - const Coefficient coefficient; - std::vector coefficient_values (n_q_points); - - typename hp::DoFHandler::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_pointactive_fe_index() == 1, ExcInternalError()); - Assert (interface_intersects_cell(cell) == true, ExcInternalError()); - - std::vector level_set_values (GeometryInfo::vertices_per_cell); - for (unsigned int v=0; v::vertices_per_cell; ++v) - level_set_values[v] = level_set (cell->vertex(v)); - - FEValues 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_pointget_fe().system_to_component_index(i).first == 0) - { - for (unsigned int j=0; jget_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; jget_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 boundary_values; - VectorTools::interpolate_boundary_values (dof_handler, - 0, - ZeroFunction(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 - std::pair > - LaplaceProblem::compute_quadrature (const Quadrature &plain_quadrature, - const typename hp::DoFHandler::active_cell_iterator &cell, - const std::vector &level_set_values ) - { - - unsigned int type = 0; - - // find the type of cut - int sign_ls[GeometryInfo::vertices_per_cell]; - for (unsigned int v=0; v::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 v0(0,0); - Point v1(1,0); - Point v2(0,1); - Point v3(1,1); - - Point A(0,0); - Point B(0,0); - Point C(0,0); - Point D(0,0); - Point E(0,0); - Point F(0,0); - - if (type == 1) - return std::pair >(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 xfem_quadrature(5*n_q_points); - - std::vector > v(GeometryInfo::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 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 > xfem_points; - std::vector 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 > 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 >(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 xfem_quadrature(5*n_q_points); - - std::vector > v(GeometryInfo::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 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 > xfem_points; - std::vector 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 > 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 >(3, xfem_quadrature); - } - - return std::pair >(0, plain_quadrature);; - - } - - template - void LaplaceProblem::append_quadrature ( const Quadrature &plain_quadrature, - const std::vector > &v, - std::vector > &xfem_points, - std::vector &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::vertices_per_cell; - - std::vector > q_points = plain_quadrature.get_points(); - std::vector > q_transf(q_points.size()); - std::vector W = plain_quadrature.get_weights(); - std::vector phi(n_v); - std::vector > grad_phi(n_v); - - const unsigned int n_q_points = plain_quadrature.size(); - - std::vector 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::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 q_prime; - for (unsigned int d=0; d::vertices_per_cell; j++) - q_prime[d] += v[j](d) * phi[j]; - xfem_points.push_back(q_prime); - } - - } - - - template - void LaplaceProblem::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 - void LaplaceProblem::refine_grid () - { - Vector estimated_error_per_cell (triangulation.n_active_cells()); - - KellyErrorEstimator::estimate (dof_handler, - QGauss(3), - typename FunctionMap::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 - class Postprocessor : public DataPostprocessor - { - public: - virtual - void - compute_derived_quantities_vector (const std::vector > &uh, - const std::vector > > &duh, - const std::vector > > &dduh, - const std::vector > &normals, - const std::vector > &evaluation_points, - std::vector > &computed_quantities) const; - - virtual std::vector get_names () const; - - virtual - std::vector - get_data_component_interpretation () const; - - virtual UpdateFlags get_needed_update_flags () const; - }; - - - template - std::vector - Postprocessor::get_names() const - { - std::vector solution_names (1, "total_solution"); - solution_names.push_back ("error"); - return solution_names; - } - - - template - std::vector - Postprocessor:: - get_data_component_interpretation () const - { - std::vector - interpretation (2, - DataComponentInterpretation::component_is_scalar); - return interpretation; - } - - - template - UpdateFlags - Postprocessor::get_needed_update_flags() const - { - return update_values | update_q_points; - } - - - template - void - Postprocessor:: - compute_derived_quantities_vector (const std::vector > &uh, - const std::vector > > &/*duh*/, - const std::vector > > &/*dduh*/, - const std::vector > &/*normals*/, - const std::vector > &evaluation_points, - std::vector > &computed_quantities) const - { - const unsigned int n_quadrature_points = uh.size(); - Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError()); - Assert (uh[0].size() == 2, ExcInternalError()); - - for (unsigned int q=0; q - void LaplaceProblem::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 postprocessor; - DataOut > 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 - void LaplaceProblem::compute_error () const - { - hp::QCollection q_collection; - q_collection.push_back (QGauss(2)); - q_collection.push_back (QIterated(QGauss<1>(2), 4)); - - hp::FEValues hp_fe_values (fe_collection, q_collection, - update_values | update_q_points | update_JxW_values); - - double l2_error_square = 0; - - std::vector > solution_values; - - typename hp::DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - - for (; cell!=endc; ++cell) - { - hp_fe_values.reinit (cell); - - const FEValues &fe_values = hp_fe_values.get_present_fe_values (); - - solution_values.resize (fe_values.n_quadrature_points, - Vector(2)); - fe_values.get_function_values (solution, - solution_values); - - for (unsigned int q=0; q boundary; - triangulation.set_boundary (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; - - deallog.depth_console (0); - - 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; -} diff --git a/deal.II/examples/step-50/CMakeLists.txt b/deal.II/examples/step-50/CMakeLists.txt deleted file mode 100644 index 3a16b1326e..0000000000 --- a/deal.II/examples/step-50/CMakeLists.txt +++ /dev/null @@ -1,49 +0,0 @@ -## -# CMake script for the step-32 tutorial program: -## - -# Set the name of the project and target: -SET(TARGET "step-50") - -# Declare all source files the target consists of: -SET(TARGET_SRC - ${TARGET}.cc - # You can specify additional files here! - ) - -# 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.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 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 fullfilled? -# -IF( NOT DEAL_II_WITH_MPI OR - NOT DEAL_II_WITH_P4EST OR - NOT DEAL_II_WITH_TRILINOS ) - 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() diff --git a/deal.II/examples/step-50/doc/builds-on b/deal.II/examples/step-50/doc/builds-on deleted file mode 100644 index 17402734c7..0000000000 --- a/deal.II/examples/step-50/doc/builds-on +++ /dev/null @@ -1 +0,0 @@ -step-6 diff --git a/deal.II/examples/step-50/doc/intro.dox b/deal.II/examples/step-50/doc/intro.dox deleted file mode 100644 index 99666726d3..0000000000 --- a/deal.II/examples/step-50/doc/intro.dox +++ /dev/null @@ -1,90 +0,0 @@ -
    - -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. - - - - -

    Introduction

    - - -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 -
      -
    • An the object handling transfer between grids; we use the - MGTransferPrebuilt class for this that does almost all of the work - inside the library. -
    • The solver on the coarsest level; here, we use MGCoarseGridHouseholder. -
    • The smoother on all other levels, which in our case will be the - MGSmootherRelaxation class using SOR as the underlying method -
    • And MGMatrix, a class having a special level multiplication, i.e. we - basically store one matrix per grid level and allow multiplication - with it. -
    -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": - -

    - @image html "hanging_nodes.png" "" -

    - -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. - - - -

    The testcase

    - -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: -- LaplaceProblem::mg_dof_handler -- LaplaceProblem::mg_sparsity -- LaplaceProblem::mg_matrices -- LaplaceProblem::mg_interface_matrices_up -- LaplaceProblem::assemble_multigrid () -- LaplaceProblem::solve () -Take a look at these functions. diff --git a/deal.II/examples/step-50/doc/kind b/deal.II/examples/step-50/doc/kind deleted file mode 100644 index c1d9154931..0000000000 --- a/deal.II/examples/step-50/doc/kind +++ /dev/null @@ -1 +0,0 @@ -techniques diff --git a/deal.II/examples/step-50/doc/results.dox b/deal.II/examples/step-50/doc/results.dox deleted file mode 100644 index 2d5fe0deb0..0000000000 --- a/deal.II/examples/step-50/doc/results.dox +++ /dev/null @@ -1,98 +0,0 @@ -

    Results

    - -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 -
    -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.
-
    -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 LaplaceProblem::refine_grid by a call to -triangulation.refine_global(1), 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: - - - - - - - - -
    cellsstep-5step-16
    20 13 6
    80 17 7
    320 29 9
    1280 51 10
    5120 94 11
    2048018013
    -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. - - -

    Possible extensions

    - -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. diff --git a/deal.II/examples/step-50/doc/tooltip b/deal.II/examples/step-50/doc/tooltip deleted file mode 100644 index 2fd65590a0..0000000000 --- a/deal.II/examples/step-50/doc/tooltip +++ /dev/null @@ -1 +0,0 @@ -Multigrid on adaptive meshes. diff --git a/deal.II/examples/step-50/step-50.cc b/deal.II/examples/step-50/step-50.cc deleted file mode 100644 index 768281e4e7..0000000000 --- a/deal.II/examples/step-50/step-50.cc +++ /dev/null @@ -1,1075 +0,0 @@ -/* --------------------------------------------------------------------- - * $Id$ - * - * Copyright (C) 2003 - 2013 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, University of Heidelberg, 2003 - * Baerbel Janssen, University of Heidelberg, 2010 - * Wolfgang Bangerth, Texas A&M University, 2010 - */ - - -// parallel geometric multigrid. work in progress! - -// As discussed in the introduction, most of -// this program is copied almost verbatim -// from step-6, which itself is only a slight -// modification of step-5. Consequently, a -// significant part of this program is not -// new if you've read all the material up to -// step-6, and we won't comment on that part -// of the functionality that is -// unchanged. Rather, we will focus on those -// aspects of the program that have to do -// with the multigrid functionality which -// forms the new aspect of this tutorial -// program. - -// @sect3{Include files} - -// Again, the first few include files -// are already known, so we won't -// comment on them: -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include - -#include -#include - -#include -#include - -#include -#include -#include - -#include -#include -#include - -#include -#include -#include - -// These, now, are the include necessary for -// the multilevel methods. The first two -// declare classes that allow us to enumerate -// degrees of freedom not only on the finest -// mesh level, but also on intermediate -// levels (that's what the MGDoFHandler class -// does) as well as allow to access this -// information (iterators and accessors over -// these cells). -// -// The rest of the include files deals with -// the mechanics of multigrid as a linear -// operator (solver or preconditioner). -#include -#include -#include -#include -#include -#include -#include -#include - -// This is C++: -#include -#include - -// The last step is as in all -// previous programs: -namespace Step50 -{ - using namespace dealii; - - - // @sect3{The LaplaceProblem class template} - - // This main class is basically the same - // class as in step-6. As far as member - // functions is concerned, the only addition - // is the assemble_multigrid - // function that assembles the matrices that - // correspond to the discrete operators on - // intermediate levels: - template - 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; - - parallel::distributed::Triangulation triangulation; - FE_Q fe; - DoFHandler mg_dof_handler; - - typedef TrilinosWrappers::SparseMatrix matrix_t; - typedef TrilinosWrappers::MPI::Vector vector_t; - - matrix_t system_matrix; - - IndexSet locally_relevant_set; - - // We need an additional object for the - // hanging nodes constraints. They are - // handed to the transfer object in the - // multigrid. Since we call a compress - // inside the multigrid these constraints - // are not allowed to be inhomogeneous so - // we store them in different ConstraintMatrix - // objects. - ConstraintMatrix hanging_node_constraints; - ConstraintMatrix constraints; - - vector_t solution; - vector_t system_rhs; - - const unsigned int degree; - - // The following four objects are the - // only additional member variables, - // compared to step-6. They first three - // represent the - // operators that act on individual - // levels of the multilevel hierarchy, - // rather than on the finest mesh as do - // the objects above while the last object - // stores information about the boundary - // indices on each level and information - // about indices lying on a refinement - // edge between two different refinement - // levels. - // - // To facilitate having objects on each - // level of a multilevel hierarchy, - // deal.II has the MGLevelObject class - // template that provides storage for - // objects on each level. What we need - // here are matrices on each level, which - // implies that we also need sparsity - // patterns on each level. As outlined in - // the @ref mg_paper, the operators - // (matrices) that we need are actually - // twofold: one on the interior of each - // level, and one at the interface - // between each level and that part of - // the domain where the mesh is - // coarser. In fact, we will need the - // latter in two versions: for the - // direction from coarse to fine mesh and - // from fine to coarse. Fortunately, - // however, we here have a self-adjoint - // problem for which one of these is the - // transpose of the other, and so we only - // have to build one; we choose the one - // from coarse to fine. - MGLevelObject mg_matrices; - MGLevelObject 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 - class Coefficient : public Function - { - public: - Coefficient () : Function() {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; - - virtual void value_list (const std::vector > &points, - std::vector &values, - const unsigned int component = 0) const; - }; - - - - template - double Coefficient::value (const Point &p, - const unsigned int) const - { - if (p.square() < 0.5*0.5) - return 20; - else - return 1; - } - - - - template - void Coefficient::value_list (const std::vector > &points, - std::vector &values, - const unsigned int component) const - { - 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::value (points[i]); - } - - - // @sect3{The LaplaceProblem 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 - LaplaceProblem::LaplaceProblem (const unsigned int degree) - : - triangulation (MPI_COMM_WORLD,Triangulation:: - limit_level_difference_at_vertices, - parallel::distributed::Triangulation::construct_multigrid_hierarchy), - fe (degree), - mg_dof_handler (triangulation), - degree(degree) - { - if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)!=0) - deallog.depth_console(0); - } - - - - // @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 - void LaplaceProblem::setup_system () - { - mg_dof_handler.distribute_dofs (fe); - mg_dof_handler.distribute_mg_dofs (fe); - - - // Here we output not only the - // degrees of freedom on the finest - // level, but also in the - // multilevel structure - deallog << "Number of degrees of freedom: " - << mg_dof_handler.n_dofs(); - - for (unsigned int l=0; l::type dirichlet_boundary; - ZeroFunction homogeneous_dirichlet_bc (1); - dirichlet_boundary[0] = &homogeneous_dirichlet_bc; - VectorTools::interpolate_boundary_values (mg_dof_handler, - dirichlet_boundary, - constraints); - constraints.close (); - hanging_node_constraints.close (); - - CompressedSimpleSparsityPattern csp(mg_dof_handler.n_dofs(), mg_dof_handler.n_dofs()); - DoFTools::make_sparsity_pattern (mg_dof_handler, csp, constraints); - system_matrix.reinit (mg_dof_handler.locally_owned_dofs(), csp, 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 dirichlet_boundary - // here as well. - mg_constrained_dofs.clear(); - mg_constrained_dofs.initialize(mg_dof_handler, dirichlet_boundary); - - - // 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 - // n_levels-1. 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 (); - mg_matrices.resize(0, n_levels-1); - mg_matrices.clear (); - - // 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 - void LaplaceProblem::assemble_system () - { - const QGauss quadrature_formula(degree+1); - - FEValues 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 cell_matrix (dofs_per_cell, dofs_per_cell); - Vector cell_rhs (dofs_per_cell); - - std::vector local_dof_indices (dofs_per_cell); - - const Coefficient coefficient; - std::vector coefficient_values (n_q_points); - - typename DoFHandler::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_pointget_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 - void LaplaceProblem::assemble_multigrid () - { - QGauss quadrature_formula(1+degree); - - FEValues 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 cell_matrix (dofs_per_cell, dofs_per_cell); - - std::vector local_dof_indices (dofs_per_cell); - - const Coefficient coefficient; - std::vector 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 setup_system 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. As in many other parts of the - // library, we do this by using Boolean masks, i.e. vectors of - // Booleans each element of which indicates whether the - // corresponding degree of freedom index is an interface DoF or - // not. The MGConstraints already computed the - // information for us when we called initialize in - // setup_system(). - std::vector > interface_dofs - = mg_constrained_dofs.get_refinement_edge_indices (); - std::vector > boundary_interface_dofs - = mg_constrained_dofs.get_refinement_edge_boundary_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 - // MGConstraints knows it. So we simply ask for them - // by calling get_boundary_indices (). 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 boundary_constraints (triangulation.n_global_levels()); - std::vector boundary_interface_constraints (triangulation.n_global_levels()); - for (unsigned int level=0; levelassemble_system, 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 - // MGDoFHandler::cell_iterator rather than - // MGDoFHandler::active_cell_iterator. Let's go about it: - typename DoFHandler::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_pointget_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 - // boundary_constraints 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 - // mg_matrices 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 - // mg_interface_matrices object. The - // boundary_interface_constraints 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 solve() function) be able to just pass - // the transpose matrix where necessary. - for (unsigned int i=0; ilevel()][local_dof_indices[i]]==true && - interface_dofs[cell->level()][local_dof_indices[j]]==false)) - cell_matrix(i,j) = 0; - - boundary_interface_constraints[cell->level()] - .distribute_local_to_global (cell_matrix, - local_dof_indices, - mg_interface_matrices[cell->level()]); - } - - for (unsigned int i=0; i - void LaplaceProblem::solve () - { - - // Create the object that deals with the transfer between - // different refinement levels. We need to pass it the hanging - // node constraints. - MGTransferPrebuilt mg_transfer(hanging_node_constraints, mg_constrained_dofs); - // Now the prolongation matrix has to be built. This matrix needs - // to take the boundary values on each level into account and - // needs to know about the indices at the refinement edges. The - // MGConstraints knows about that so pass it as an - // argument. - mg_transfer.build_matrices(mg_dof_handler); - - matrix_t &coarse_matrix = mg_matrices[0]; - //coarse_matrix.copy_from (mg_matrices[0]); - //MGCoarseGridHouseholder coarse_grid_solver; - //coarse_grid_solver.initialize (coarse_matrix); - - SolverControl coarse_solver_control (1000, 1e-10, false, false); - SolverGMRES coarse_solver(coarse_solver_control); - PreconditionIdentity id; - MGCoarseGridLACIteration,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 typedef 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 - // initialize() 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 - // set_steps() 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 mg_smoother - // 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 TrilinosWrappers::PreconditionJacobi Smoother; - MGSmootherPrecondition mg_smoother; - mg_smoother.initialize(mg_matrices); - 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: - MGMatrix mg_matrix(&mg_matrices); - MGMatrix mg_interface_up(&mg_interface_matrices); - MGMatrix mg_interface_down(&mg_interface_matrices); - - // Now, we are ready to set up the - // V-cycle operator and the - // multilevel preconditioner. - Multigrid 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 > - 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 cg (solver_control); - - solution = 0; - - 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); - cg.solve (system_matrix, solution, system_rhs, - prec); - } - else - { - cg.solve (system_matrix, solution, system_rhs, - preconditioner); - } - - 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 functions are almost - // unchanged from those in step-6, with the - // exception of two minor differences: The - // KellyErrorEstimator::estimate function - // wants an argument of type DoFHandler, not - // MGDoFHandler, and so we have to cast from - // derived to base class; and we generate - // output in VTK format, to use the more - // modern visualization programs available - // today compared to those that were - // available when step-6 was written. - template - void LaplaceProblem::refine_grid () - { - Vector estimated_error_per_cell (triangulation.n_active_cells()); - - TrilinosWrappers::MPI::Vector temp_solution; - temp_solution.reinit(locally_relevant_set, MPI_COMM_WORLD); - temp_solution = solution; - - KellyErrorEstimator::estimate (static_cast&>(mg_dof_handler), - QGauss(3), - typename FunctionMap::type(), - temp_solution, - estimated_error_per_cell); - parallel::distributed::GridRefinement:: - refine_and_coarsen_fixed_fraction (triangulation, - estimated_error_per_cell, - 0.3, 0.03); - triangulation.execute_coarsening_and_refinement (); - } - - - - template - void LaplaceProblem::output_results (const unsigned int cycle) const - { - DataOut data_out; - - TrilinosWrappers::MPI::Vector temp_solution; - temp_solution.reinit(locally_relevant_set, MPI_COMM_WORLD); - temp_solution = solution; - - - TrilinosWrappers::MPI::Vector temp = solution; - system_matrix.residual(temp,solution,system_rhs); - TrilinosWrappers::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 subdomain (triangulation.n_active_cells()); - for (unsigned int i=0; i filenames; - for (unsigned int i=0; iassemble_multigrid that takes - // care of forming the matrices on every - // level that we need in the multigrid - // method. - template - void LaplaceProblem::run () - { - for (unsigned int cycle=0; cycle<13; ++cycle) - { - deallog << "Cycle " << cycle << ':' << std::endl; - - if (cycle == 0) - { - GridGenerator::hyper_cube (triangulation); - - triangulation.refine_global (3); - } - else - refine_grid (); - - deallog << " Number of active cells: " - << triangulation.n_global_active_cells() - << std::endl; - - setup_system (); - - deallog << " Number of degrees of freedom: " - << mg_dof_handler.n_dofs() - << " (by level: "; - for (unsigned int level=0; level laplace_problem(1); - 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; - - 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; -}