/* $Id$ */
/* */
-/* Copyright (C) 2009, 2010 by Timo Heister and the deal.II authors */
+/* Copyright (C) 2009, 2010, 2011 by Timo Heister and the deal.II authors */
/* */
/* This file is subject to QPL and may not be distributed */
/* without copyright and license information. Please refer */
#include <fstream>
#include <iostream>
-using namespace dealii;
-
- // @sect3{The <code>LaplaceProblem</code> class template}
-
- // Next let's declare the main class of this
- // program. Its structure is almost exactly
- // that of the step-6 tutorial program. The
- // only significant differences are:
- // - The <code>mpi_communicator</code>
- // variable that describes the set of
- // processors we want this code to run
- // on. In practice, this will be
- // MPI_COMM_WORLD, i.e. all processors the
- // batch scheduling system has assigned to
- // this particular job.
- // - The presence of the <code>pcout</code>
- // variable of type ConditionOStream.
- // - The obvious use of
- // parallel::distributed::Triangulation
- // instead of Triangulation.
- // - The presence of two IndexSet objects
- // that denote which sets of degrees of
- // freedom (and associated elements of
- // solution and right hand side vectors) we
- // own on the current processor and which
- // we need (as ghost elements) for the
- // algorithms in this program to work.
- // - The fact that all matrices and
- // vectors are now distributed. We
- // use their PETScWrapper versions
- // for this since deal.II's own
- // classes do not provide %parallel
- // functionality. Note that as part
- // of this class, we store a
- // solution vector that does not
- // only contain the degrees of
- // freedom the current processor
- // owns, but also (as ghost
- // elements) all those vector
- // elements that correspond to
- // "locally relevant" degrees of
- // freedom (i.e. all those that
- // live on locally owned cells or
- // the layer of ghost cells that
- // surround it).
-template <int dim>
-class LaplaceProblem
+namespace Step40
{
- public:
- LaplaceProblem ();
- ~LaplaceProblem ();
-
- void run ();
-
- private:
- void setup_system ();
- void assemble_system ();
- void solve ();
- void refine_grid ();
- void output_results (const unsigned int cycle) const;
-
- MPI_Comm mpi_communicator;
-
- parallel::distributed::Triangulation<dim> triangulation;
-
- DoFHandler<dim> dof_handler;
- FE_Q<dim> fe;
-
- IndexSet locally_owned_dofs;
- IndexSet locally_relevant_dofs;
-
- ConstraintMatrix constraints;
+ using namespace dealii;
+
+ // @sect3{The <code>LaplaceProblem</code> class template}
+
+ // Next let's declare the main class of this
+ // program. Its structure is almost exactly
+ // that of the step-6 tutorial program. The
+ // only significant differences are:
+ // - The <code>mpi_communicator</code>
+ // variable that describes the set of
+ // processors we want this code to run
+ // on. In practice, this will be
+ // MPI_COMM_WORLD, i.e. all processors the
+ // batch scheduling system has assigned to
+ // this particular job.
+ // - The presence of the <code>pcout</code>
+ // variable of type ConditionOStream.
+ // - The obvious use of
+ // parallel::distributed::Triangulation
+ // instead of Triangulation.
+ // - The presence of two IndexSet objects
+ // that denote which sets of degrees of
+ // freedom (and associated elements of
+ // solution and right hand side vectors) we
+ // own on the current processor and which
+ // we need (as ghost elements) for the
+ // algorithms in this program to work.
+ // - The fact that all matrices and
+ // vectors are now distributed. We
+ // use their PETScWrapper versions
+ // for this since deal.II's own
+ // classes do not provide %parallel
+ // functionality. Note that as part
+ // of this class, we store a
+ // solution vector that does not
+ // only contain the degrees of
+ // freedom the current processor
+ // owns, but also (as ghost
+ // elements) all those vector
+ // elements that correspond to
+ // "locally relevant" degrees of
+ // freedom (i.e. all those that
+ // live on locally owned cells or
+ // the layer of ghost cells that
+ // surround it).
+ template <int dim>
+ class LaplaceProblem
+ {
+ public:
+ LaplaceProblem ();
+ ~LaplaceProblem ();
+
+ void run ();
+
+ private:
+ void setup_system ();
+ void assemble_system ();
+ void solve ();
+ void refine_grid ();
+ void output_results (const unsigned int cycle) const;
+
+ MPI_Comm mpi_communicator;
+
+ parallel::distributed::Triangulation<dim> triangulation;
+
+ DoFHandler<dim> dof_handler;
+ FE_Q<dim> fe;
+
+ IndexSet locally_owned_dofs;
+ IndexSet locally_relevant_dofs;
+
+ ConstraintMatrix constraints;
+
+ PETScWrappers::MPI::SparseMatrix system_matrix;
+ PETScWrappers::MPI::Vector locally_relevant_solution;
+ PETScWrappers::MPI::Vector system_rhs;
+
+ ConditionalOStream pcout;
+ };
+
+
+ // @sect3{The <code>LaplaceProblem</code> class implementation}
+
+ // @sect4{Constructors and destructors}
+
+ // Constructors and destructors are rather
+ // trivial. In addition to what we do in
+ // step-6, we set the set of processors we
+ // want to work on to all machines available
+ // (MPI_COMM_WORLD); ask the triangulation to
+ // ensure that the mesh remains smooth and
+ // free to refined islands, for example; and
+ // initialize the <code>pcout</code> variable
+ // to only allow processor zero to output
+ // anything:
+ template <int dim>
+ LaplaceProblem<dim>::LaplaceProblem ()
+ :
+ mpi_communicator (MPI_COMM_WORLD),
+ triangulation (mpi_communicator,
+ typename Triangulation<dim>::MeshSmoothing
+ (Triangulation<dim>::smoothing_on_refinement |
+ Triangulation<dim>::smoothing_on_coarsening)),
+ dof_handler (triangulation),
+ fe (2),
+ pcout (std::cout,
+ (Utilities::System::
+ get_this_mpi_process(mpi_communicator)
+ == 0))
+ {}
+
+
+
+ template <int dim>
+ LaplaceProblem<dim>::~LaplaceProblem ()
+ {
+ dof_handler.clear ();
+ }
+
+
+ // @sect4{LaplaceProblem::setup_system}
+
+ // The following function is, arguably, the
+ // most interesting one in the entire program
+ // since it goes to the heart of what
+ // distinguishes %parallel step-40 from
+ // sequential step-6.
+ //
+ // At the top we do what we always do: tell
+ // the DoFHandler object to distribute
+ // degrees of freedom. Since the
+ // triangulation we use here is distributed,
+ // the DoFHandler object is smart enough to
+ // recognize that on each processor it can
+ // only distribute degrees of freedom on
+ // cells it owns; this is followed by an
+ // exchange step in which processors tell
+ // each other about degrees of freedom on
+ // ghost cell. The result is a DoFHandler
+ // that knows about the degrees of freedom on
+ // locally owned cells and ghost cells
+ // (i.e. cells adjacent to locally owned
+ // cells) but nothing about cells that are
+ // further away, consistent with the basic
+ // philosophy of distributed computing that
+ // no processor can know everything.
+ template <int dim>
+ void LaplaceProblem<dim>::setup_system ()
+ {
+ dof_handler.distribute_dofs (fe);
+
+ // The next two lines extract some
+ // informatino we will need later
+ // on, namely two index sets that
+ // provide information about which
+ // degrees of freedom are owned by
+ // the current processor (this
+ // information will be used to
+ // initialize solution and right
+ // hand side vectors, and the
+ // system matrix, indicating which
+ // elements to store on the current
+ // processor and which to expect to
+ // be stored somewhere else); and
+ // an index set that indicates
+ // which degrees of freedom are
+ // locally relevant (i.e. live on
+ // cells that the current processor
+ // owns or on the layer of ghost
+ // cells around the locally owned
+ // cells; we need all of these
+ // degrees of freedom, for example,
+ // to estimate the error on the
+ // local cells).
+ locally_owned_dofs = dof_handler.locally_owned_dofs ();
+ DoFTools::extract_locally_relevant_dofs (dof_handler,
+ locally_relevant_dofs);
+
+ // Next, let us initialize the
+ // solution and right hand side
+ // vectors. As mentioned above, the
+ // solution vector we seek does not
+ // only store elements we own, but
+ // also ghost entries; on the other
+ // hand, the right hand side vector
+ // only needs to have the entries
+ // the current processor owns since
+ // all we will ever do is write
+ // into it, never read from it on
+ // locally owned cells (of course
+ // the linear solvers will read
+ // from it, but they do not care
+ // about the geometric location of
+ // degrees of freedom).
+ locally_relevant_solution.reinit (mpi_communicator,
+ locally_owned_dofs,
+ locally_relevant_dofs);
+ locally_relevant_solution = 0;
+ system_rhs.reinit (mpi_communicator,
+ dof_handler.n_dofs(),
+ dof_handler.n_locally_owned_dofs());
+ system_rhs = 0;
+
+ // The next step is to compute hanging node
+ // and boundary value constraints, which we
+ // combine into a single object storing all
+ // constraints.
+ //
+ // As with all other things in %parallel,
+ // the mantra must be that no processor can
+ // store all information about the entire
+ // universe. As a consequence, we need to
+ // tell the constraints object for which
+ // degrees of freedom it can store
+ // constraints and for which it may not
+ // expect any information to store. In our
+ // case, as explained in the @ref
+ // distributed module, the degrees of
+ // freedom we need to care about on each
+ // processor are the locally relevant ones,
+ // so we pass this to the
+ // ConstraintMatrix::reinit function. As a
+ // side note, if you forget to pass this
+ // argument, the ConstraintMatrix class
+ // will allocate an array with length equal
+ // to the largest DoF index it has seen so
+ // far. For processors with high MPI
+ // process number, this may be very large
+ // -- maybe on the order of billions. The
+ // program would then allocate more memory
+ // than for likely all other operations
+ // combined for this single array.
+ constraints.clear ();
+ constraints.reinit (locally_relevant_dofs);
+ DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 0,
+ ZeroFunction<dim>(),
+ constraints);
+ constraints.close ();
+
+ // The last part of this function deals
+ // with initializing the matrix with
+ // accompanying sparsity pattern. As in
+ // previous tutorial programs, we use the
+ // CompressedSimpleSparsityPattern as an
+ // intermediate with which we then
+ // initialize the PETSc matrix. To do so we
+ // have to tell the sparsity pattern its
+ // size but as above there is no way the
+ // resulting object will be able to store
+ // even a single pointer for each global
+ // degree of freedom; the best we can hope
+ // for is that it stores information about
+ // each locally relevant degree of freedom,
+ // i.e. all those that we may ever touch in
+ // the process of assembling the matrix
+ // (the @ref distributed_paper
+ // "distributed computing paper" has a long
+ // discussion why one really needs the
+ // locally relevant, and not the small set
+ // of locally active degrees of freedom in
+ // this context).
+ //
+ // So we tell the sparsity pattern its size
+ // and what DoFs to store anything for and
+ // then ask DoFTools::make_sparsity_pattern
+ // to fill it (this function ignores all
+ // cells that are not locally owned,
+ // mimicking what we will do below in the
+ // assembly process). After this, we call a
+ // function that exchanges entries in these
+ // sparsity pattern between processors so
+ // that in the end each processor really
+ // knows about all the entries that will
+ // exist in that part of the finite element
+ // matrix that it will own. The final step
+ // is to initialize the matrix with the
+ // sparsity pattern.
+ CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
+ dof_handler.n_dofs(),
+ locally_relevant_dofs);
+ DoFTools::make_sparsity_pattern (dof_handler,
+ csp,
+ constraints, false);
+ SparsityTools::distribute_sparsity_pattern (csp,
+ dof_handler.n_locally_owned_dofs_per_processor(),
+ mpi_communicator,
+ locally_relevant_dofs);
+ system_matrix.reinit (mpi_communicator,
+ csp,
+ dof_handler.n_locally_owned_dofs_per_processor(),
+ dof_handler.n_locally_owned_dofs_per_processor(),
+ Utilities::System::get_this_mpi_process(mpi_communicator));
+ }
+
+
+
+ // @sect4{LaplaceProblem::assemble_system}
+
+ // The function that then assembles the
+ // linear system is comparatively boring,
+ // being almost exactly what we've seen
+ // before. The points to watch out for are:
+ // - Assembly must only loop over locally
+ // owned cells. We test this by comparing
+ // a cell's subdomain_id against
+ // information from the triangulation
+ // but an equally valid condition would
+ // have been to skip all cells for which
+ // the condition <code>cell->is_ghost()
+ // || cell->is_artificial()</code> is
+ // true.
+ // - Copying local contributions into the
+ // global matrix must include distributing
+ // constraints and boundary values. In
+ // other words, we can now (as we did in
+ // step-6) first copy every local
+ // contribution into the global matrix and
+ // only in a later step take care of
+ // hanging node constraints and boundary
+ // values. The reason is, as discussed in
+ // step-17, that PETSc does not provide
+ // access to arbitrary elements of the
+ // matrix once they have been assembled
+ // into it -- in parts because they may
+ // simple no longer reside on the current
+ // processor but have instead been shipped
+ // to a different machine.
+ // - The way we compute the right hand side
+ // (given the formula stated in the
+ // introduction) may not be the most
+ // elegant but will do for a program whose
+ // focus lies somewhere entirely different.
+ template <int dim>
+ void LaplaceProblem<dim>::assemble_system ()
+ {
+ const QGauss<dim> quadrature_formula(3);
+
+ 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<unsigned int> local_dof_indices (dofs_per_cell);
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
+ if (cell->subdomain_id() == triangulation.locally_owned_subdomain())
+ {
+ cell_matrix = 0;
+ cell_rhs = 0;
+
+ fe_values.reinit (cell);
+
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ {
+ const double
+ rhs_value
+ = (fe_values.quadrature_point(q_point)[1]
+ >
+ 0.5+0.25*std::sin(4.0 * numbers::PI *
+ fe_values.quadrature_point(q_point)[0])
+ ? 1 : -1);
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
+ fe_values.shape_grad(j,q_point) *
+ fe_values.JxW(q_point));
+
+ cell_rhs(i) += (rhs_value *
+ fe_values.shape_value(i,q_point) *
+ 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);
+ }
- PETScWrappers::MPI::SparseMatrix system_matrix;
- PETScWrappers::MPI::Vector locally_relevant_solution;
- PETScWrappers::MPI::Vector system_rhs;
+ system_matrix.compress ();
+ system_rhs.compress ();
+ }
+
+
+
+ // @sect4{LaplaceProblem::solve}
+
+ // Even though solving linear systems
+ // on potentially tens of thousands
+ // of processors is by far not a
+ // trivial job, the function that
+ // does this is -- at least at the
+ // outside -- relatively simple. Most
+ // of the parts you've seen
+ // before. There are really only two
+ // things worth mentioning:
+ // - Solvers and preconditioners are
+ // built on the deal.II wrappers of
+ // PETSc functionality. It is
+ // relatively well known that the
+ // primary bottleneck of massively
+ // %parallel linear solvers is not
+ // actually the communication
+ // between processors, but the fact
+ // that it is difficult to produce
+ // preconditioners that scale well
+ // to large numbers of
+ // processors. Over the second half
+ // of the first decade of the 21st
+ // century, it has become clear
+ // that algebraic multigrid (AMG)
+ // methods turn out to be extremely
+ // efficient in this context, and
+ // we will use one of them -- the
+ // BoomerAMG implementation of the
+ // Hypre package that can be
+ // interfaced to through PETSc --
+ // for the current program. The
+ // rest of the solver itself is
+ // boilerplate and has been shown
+ // before. Since the linear system
+ // is symmetric and positive
+ // definite, we can use the CG
+ // method as the outer solver.
+ // - Ultimately, we want a vector
+ // that stores not only the
+ // elements of the solution for
+ // degrees of freedom the current
+ // processor owns, but also all
+ // other locally relevant degrees
+ // of freedom. On the other hand,
+ // the solver itself needs a vector
+ // that is uniquely split between
+ // processors, without any
+ // overlap. We therefore create a
+ // vector at the beginning of this
+ // function that has these
+ // properties, use it to solve the
+ // linear system, and only assign
+ // it to the vector we want at the
+ // very end. This last step ensures
+ // that all ghost elements are also
+ // copied as necessary.
+ template <int dim>
+ void LaplaceProblem<dim>::solve ()
+ {
+ PETScWrappers::MPI::Vector
+ completely_distributed_solution (mpi_communicator,
+ dof_handler.n_dofs(),
+ dof_handler.n_locally_owned_dofs());
- ConditionalOStream pcout;
-};
+ SolverControl solver_control (dof_handler.n_dofs(), 1e-12);
+ PETScWrappers::SolverCG solver(solver_control, mpi_communicator);
- // @sect3{The <code>LaplaceProblem</code> class implementation}
+ // Ask for a symmetric preconditioner by
+ // setting the first parameter in
+ // AdditionalData to true.
+ PETScWrappers::PreconditionBoomerAMG
+ preconditioner(system_matrix,
+ PETScWrappers::PreconditionBoomerAMG::AdditionalData(true));
- // @sect4{Constructors and destructors}
+ solver.solve (system_matrix, completely_distributed_solution, system_rhs,
+ preconditioner);
- // Constructors and destructors are rather
- // trivial. In addition to what we do in
- // step-6, we set the set of processors we
- // want to work on to all machines available
- // (MPI_COMM_WORLD); ask the triangulation to
- // ensure that the mesh remains smooth and
- // free to refined islands, for example; and
- // initialize the <code>pcout</code> variable
- // to only allow processor zero to output
- // anything:
-template <int dim>
-LaplaceProblem<dim>::LaplaceProblem ()
- :
- mpi_communicator (MPI_COMM_WORLD),
- triangulation (mpi_communicator,
- typename Triangulation<dim>::MeshSmoothing
- (Triangulation<dim>::smoothing_on_refinement |
- Triangulation<dim>::smoothing_on_coarsening)),
- dof_handler (triangulation),
- fe (2),
- pcout (std::cout,
- (Utilities::System::
- get_this_mpi_process(mpi_communicator)
- == 0))
-{}
+ pcout << " Solved in " << solver_control.last_step()
+ << " iterations." << std::endl;
+ constraints.distribute (completely_distributed_solution);
+ locally_relevant_solution = completely_distributed_solution;
+ locally_relevant_solution.update_ghost_values();
+ }
-template <int dim>
-LaplaceProblem<dim>::~LaplaceProblem ()
-{
- dof_handler.clear ();
-}
- // @sect4{LaplaceProblem::setup_system}
+ // @sect4{LaplaceProblem::refine_grid}
- // The following function is, arguably, the
- // most interesting one in the entire program
- // since it goes to the heart of what
- // distinguishes %parallel step-40 from
- // sequential step-6.
- //
- // At the top we do what we always do: tell
- // the DoFHandler object to distribute
- // degrees of freedom. Since the
- // triangulation we use here is distributed,
- // the DoFHandler object is smart enough to
- // recognize that on each processor it can
- // only distribute degrees of freedom on
- // cells it owns; this is followed by an
- // exchange step in which processors tell
- // each other about degrees of freedom on
- // ghost cell. The result is a DoFHandler
- // that knows about the degrees of freedom on
- // locally owned cells and ghost cells
- // (i.e. cells adjacent to locally owned
- // cells) but nothing about cells that are
- // further away, consistent with the basic
- // philosophy of distributed computing that
- // no processor can know everything.
-template <int dim>
-void LaplaceProblem<dim>::setup_system ()
-{
- dof_handler.distribute_dofs (fe);
-
- // The next two lines extract some
- // informatino we will need later
- // on, namely two index sets that
- // provide information about which
- // degrees of freedom are owned by
- // the current processor (this
- // information will be used to
- // initialize solution and right
- // hand side vectors, and the
- // system matrix, indicating which
- // elements to store on the current
- // processor and which to expect to
- // be stored somewhere else); and
- // an index set that indicates
- // which degrees of freedom are
- // locally relevant (i.e. live on
- // cells that the current processor
- // owns or on the layer of ghost
- // cells around the locally owned
- // cells; we need all of these
- // degrees of freedom, for example,
- // to estimate the error on the
- // local cells).
- locally_owned_dofs = dof_handler.locally_owned_dofs ();
- DoFTools::extract_locally_relevant_dofs (dof_handler,
- locally_relevant_dofs);
-
- // Next, let us initialize the
- // solution and right hand side
- // vectors. As mentioned above, the
- // solution vector we seek does not
- // only store elements we own, but
- // also ghost entries; on the other
- // hand, the right hand side vector
- // only needs to have the entries
- // the current processor owns since
- // all we will ever do is write
- // into it, never read from it on
- // locally owned cells (of course
- // the linear solvers will read
- // from it, but they do not care
- // about the geometric location of
- // degrees of freedom).
- locally_relevant_solution.reinit (mpi_communicator,
- locally_owned_dofs,
- locally_relevant_dofs);
- locally_relevant_solution = 0;
- system_rhs.reinit (mpi_communicator,
- dof_handler.n_dofs(),
- dof_handler.n_locally_owned_dofs());
- system_rhs = 0;
-
- // The next step is to compute hanging node
- // and boundary value constraints, which we
- // combine into a single object storing all
- // constraints.
+ // The function that estimates the
+ // error and refines the grid is
+ // again almost exactly like the one
+ // in step-6. The only difference is
+ // that the function that flags cells
+ // to be refined is now in namespace
+ // parallel::distributed::GridRefinement
+ // -- a namespace that has functions
+ // that can communicate between all
+ // involved processors and determine
+ // global thresholds to use in
+ // deciding which cells to refine and
+ // which to coarsen.
//
- // As with all other things in %parallel,
- // the mantra must be that no processor can
- // store all information about the entire
- // universe. As a consequence, we need to
- // tell the constraints object for which
- // degrees of freedom it can store
- // constraints and for which it may not
- // expect any information to store. In our
- // case, as explained in the @ref
- // distributed module, the degrees of
- // freedom we need to care about on each
- // processor are the locally relevant ones,
- // so we pass this to the
- // ConstraintMatrix::reinit function. As a
- // side note, if you forget to pass this
- // argument, the ConstraintMatrix class
- // will allocate an array with length equal
- // to the largest DoF index it has seen so
- // far. For processors with high MPI
- // process number, this may be very large
- // -- maybe on the order of billions. The
- // program would then allocate more memory
- // than for likely all other operations
- // combined for this single array.
- constraints.clear ();
- constraints.reinit (locally_relevant_dofs);
- DoFTools::make_hanging_node_constraints (dof_handler, constraints);
- VectorTools::interpolate_boundary_values (dof_handler,
- 0,
- ZeroFunction<dim>(),
- constraints);
- constraints.close ();
-
- // The last part of this function deals
- // with initializing the matrix with
- // accompanying sparsity pattern. As in
- // previous tutorial programs, we use the
- // CompressedSimpleSparsityPattern as an
- // intermediate with which we then
- // initialize the PETSc matrix. To do so we
- // have to tell the sparsity pattern its
- // size but as above there is no way the
- // resulting object will be able to store
- // even a single pointer for each global
- // degree of freedom; the best we can hope
- // for is that it stores information about
- // each locally relevant degree of freedom,
- // i.e. all those that we may ever touch in
- // the process of assembling the matrix
- // (the @ref distributed_paper
- // "distributed computing paper" has a long
- // discussion why one really needs the
- // locally relevant, and not the small set
- // of locally active degrees of freedom in
- // this context).
+ // Note that we didn't have to do
+ // anything special about the
+ // KellyErrorEstimator class: we just
+ // give it a vector with as many
+ // elements as the local
+ // triangulation has cells (locally
+ // owned cells, ghost cells, and
+ // artificial ones), but it only
+ // fills those entries that
+ // correspond to cells that are
+ // locally owned.
+ 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(),
+ locally_relevant_solution,
+ estimated_error_per_cell);
+ parallel::distributed::GridRefinement::
+ refine_and_coarsen_fixed_number (triangulation,
+ estimated_error_per_cell,
+ 0.3, 0.03);
+ triangulation.execute_coarsening_and_refinement ();
+ }
+
+
+
+ // @sect4{LaplaceProblem::output_results}
+
+ // Compared to the corresponding
+ // function in step-6, the one here
+ // is a tad more complicated. There
+ // are two reasons: the first one is
+ // that we do not just want to output
+ // the solution but also for each
+ // cell which processor owns it
+ // (i.e. which "subdomain" it is
+ // in). Secondly, as discussed at
+ // length in step-17 and step-18,
+ // generating graphical data can be a
+ // bottleneck in parallelizing. In
+ // step-18, we have moved this step
+ // out of the actual computation but
+ // shifted it into a separate program
+ // that later combined the output
+ // from various processors into a
+ // single file. But this doesn't
+ // scale: if the number of processors
+ // is large, this may mean that the
+ // step of combining data on a single
+ // processor later becomes the
+ // longest running part of the
+ // program, or it may produce a file
+ // that's so large that it can't be
+ // visualized any more. We here
+ // follow a more sensible approach,
+ // namely creating individual files
+ // for each MPI process and leaving
+ // it to the visualization program to
+ // make sense of that.
//
- // So we tell the sparsity pattern its size
- // and what DoFs to store anything for and
- // then ask DoFTools::make_sparsity_pattern
- // to fill it (this function ignores all
- // cells that are not locally owned,
- // mimicking what we will do below in the
- // assembly process). After this, we call a
- // function that exchanges entries in these
- // sparsity pattern between processors so
- // that in the end each processor really
- // knows about all the entries that will
- // exist in that part of the finite element
- // matrix that it will own. The final step
- // is to initialize the matrix with the
- // sparsity pattern.
- CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
- dof_handler.n_dofs(),
- locally_relevant_dofs);
- DoFTools::make_sparsity_pattern (dof_handler,
- csp,
- constraints, false);
- SparsityTools::distribute_sparsity_pattern (csp,
- dof_handler.n_locally_owned_dofs_per_processor(),
- mpi_communicator,
- locally_relevant_dofs);
- system_matrix.reinit (mpi_communicator,
- csp,
- dof_handler.n_locally_owned_dofs_per_processor(),
- dof_handler.n_locally_owned_dofs_per_processor(),
- Utilities::System::get_this_mpi_process(mpi_communicator));
-}
-
-
-
- // @sect4{LaplaceProblem::assemble_system}
-
- // The function that then assembles the
- // linear system is comparatively boring,
- // being almost exactly what we've seen
- // before. The points to watch out for are:
- // - Assembly must only loop over locally
- // owned cells. We test this by comparing
- // a cell's subdomain_id against
- // information from the triangulation
- // but an equally valid condition would
- // have been to skip all cells for which
- // the condition <code>cell->is_ghost()
- // || cell->is_artificial()</code> is
- // true.
- // - Copying local contributions into the
- // global matrix must include distributing
- // constraints and boundary values. In
- // other words, we can now (as we did in
- // step-6) first copy every local
- // contribution into the global matrix and
- // only in a later step take care of
- // hanging node constraints and boundary
- // values. The reason is, as discussed in
- // step-17, that PETSc does not provide
- // access to arbitrary elements of the
- // matrix once they have been assembled
- // into it -- in parts because they may
- // simple no longer reside on the current
- // processor but have instead been shipped
- // to a different machine.
- // - The way we compute the right hand side
- // (given the formula stated in the
- // introduction) may not be the most
- // elegant but will do for a program whose
- // focus lies somewhere entirely different.
-template <int dim>
-void LaplaceProblem<dim>::assemble_system ()
-{
- const QGauss<dim> quadrature_formula(3);
-
- 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<unsigned int> local_dof_indices (dofs_per_cell);
-
- typename DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active(),
- endc = dof_handler.end();
- for (; cell!=endc; ++cell)
- if (cell->subdomain_id() == triangulation.locally_owned_subdomain())
+ // To start, the top of the function
+ // looks like always. In addition to
+ // attaching the solution vector (the
+ // one that has entries for all
+ // locally relevant, not only the
+ // locally owned, elements), we
+ // attach a data vector that stores,
+ // for each cell, the subdomain the
+ // cell belongs to. This is slightly
+ // tricky, because of course not
+ // every processor knows about every
+ // cell. The vector we attach
+ // therefore has an entry for every
+ // cell that the current processor
+ // has in its mesh (locally owned
+ // onces, ghost cells, and artificial
+ // cells), but the DataOut class will
+ // ignore all entries that correspond
+ // to cells that are not owned by the
+ // current processor. As a
+ // consequence, it doesn't actually
+ // matter what values we write into
+ // these vector entries: we simply
+ // fill the entire vector with the
+ // number of the current MPI process
+ // (i.e. the subdomain_id of the
+ // current process); this correctly
+ // sets the values we care for,
+ // i.e. the entries that correspond
+ // to locally owned cells, while
+ // providing the wrong value for all
+ // other elements -- but these are
+ // then ignored anyway.
+ template <int dim>
+ void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
+ {
+ DataOut<dim> data_out;
+ data_out.attach_dof_handler (dof_handler);
+ data_out.add_data_vector (locally_relevant_solution, "u");
+
+ 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 ();
+
+ // The next step is to write this
+ // data to disk. We choose file
+ // names of the form
+ // <code>solution-XX-PPPP.vtu</code>
+ // where <code>XX</code> indicates
+ // the refinement cycle,
+ // <code>PPPP</code> refers to the
+ // processor number (enough for up
+ // to 10,000 processors, though we
+ // hope that nobody ever tries to
+ // generate this much data -- you
+ // would likely overflow all file
+ // system quotas), and
+ // <code>.vtu</code> indicates the
+ // XML-based Visualization Toolkit
+ // (VTK) file format.
+ const std::string filename = ("solution-" +
+ Utilities::int_to_string (cycle, 2) +
+ "." +
+ Utilities::int_to_string
+ (triangulation.locally_owned_subdomain(), 4));
+ std::ofstream output ((filename + ".vtu").c_str());
+ data_out.write_vtu (output);
+
+ // The last step is to write a
+ // "master record" that lists for
+ // the visualization program the
+ // names of the various files that
+ // combined represents the
+ // graphical data for the entire
+ // domain. The
+ // DataOutBase::write_pvtu_record
+ // does this, and it needs a list
+ // of filenames that we create
+ // first. Note that only one
+ // processor needs to generate this
+ // file; we arbitrarily choose
+ // processor zero to take over this
+ // job.
+ if (Utilities::System::get_this_mpi_process(mpi_communicator) == 0)
{
- cell_matrix = 0;
- cell_rhs = 0;
-
- fe_values.reinit (cell);
+ std::vector<std::string> filenames;
+ for (unsigned int i=0;
+ i<Utilities::System::get_n_mpi_processes(mpi_communicator);
+ ++i)
+ filenames.push_back ("solution-" +
+ Utilities::int_to_string (cycle, 2) +
+ "." +
+ Utilities::int_to_string (i, 4) +
+ ".vtu");
+
+ std::ofstream master_output ((filename + ".pvtu").c_str());
+ data_out.write_pvtu_record (master_output, filenames);
+ }
+ }
+
+
+
+ // @sect4{LaplaceProblem::run}
+
+ // The function that controls the
+ // overall behavior of the program is
+ // again like the one in step-6. The
+ // minor difference are the use of
+ // <code>pcout</code> instead of
+ // <code>std::cout</code> for output
+ // to the console (see also step-17)
+ // and that we only generate
+ // graphical output if at most 32
+ // processors are involved. Without
+ // this limit, it would be just too
+ // easy for people carelessly running
+ // this program without reading it
+ // first to bring down the cluster
+ // interconnect and fill any file
+ // system available :-)
+ //
+ // A functional difference to step-6
+ // is the use of a square domain and
+ // that we start with a slightly
+ // finer mesh (5 global refinement
+ // cycles) -- there just isn't much
+ // of a point showing a massively
+ // %parallel program starting on 4
+ // cells (although admittedly the
+ // point is only slightly stronger
+ // starting on 1024).
+ template <int dim>
+ void LaplaceProblem<dim>::run ()
+ {
+ const unsigned int n_cycles = 8;
+ for (unsigned int cycle=0; cycle<n_cycles; ++cycle)
+ {
+ pcout << "Cycle " << cycle << ':' << std::endl;
- for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ if (cycle == 0)
{
- const double
- rhs_value
- = (fe_values.quadrature_point(q_point)[1]
- >
- 0.5+0.25*std::sin(4.0 * numbers::PI *
- fe_values.quadrature_point(q_point)[0])
- ? 1 : -1);
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
- fe_values.shape_grad(j,q_point) *
- fe_values.JxW(q_point));
-
- cell_rhs(i) += (rhs_value *
- fe_values.shape_value(i,q_point) *
- fe_values.JxW(q_point));
- }
+ GridGenerator::hyper_cube (triangulation);
+ triangulation.refine_global (5);
}
+ else
+ refine_grid ();
- 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 ();
- system_rhs.compress ();
-}
-
-
-
- // @sect4{LaplaceProblem::solve}
-
- // Even though solving linear systems
- // on potentially tens of thousands
- // of processors is by far not a
- // trivial job, the function that
- // does this is -- at least at the
- // outside -- relatively simple. Most
- // of the parts you've seen
- // before. There are really only two
- // things worth mentioning:
- // - Solvers and preconditioners are
- // built on the deal.II wrappers of
- // PETSc functionality. It is
- // relatively well known that the
- // primary bottleneck of massively
- // %parallel linear solvers is not
- // actually the communication
- // between processors, but the fact
- // that it is difficult to produce
- // preconditioners that scale well
- // to large numbers of
- // processors. Over the second half
- // of the first decade of the 21st
- // century, it has become clear
- // that algebraic multigrid (AMG)
- // methods turn out to be extremely
- // efficient in this context, and
- // we will use one of them -- the
- // BoomerAMG implementation of the
- // Hypre package that can be
- // interfaced to through PETSc --
- // for the current program. The
- // rest of the solver itself is
- // boilerplate and has been shown
- // before. Since the linear system
- // is symmetric and positive
- // definite, we can use the CG
- // method as the outer solver.
- // - Ultimately, we want a vector
- // that stores not only the
- // elements of the solution for
- // degrees of freedom the current
- // processor owns, but also all
- // other locally relevant degrees
- // of freedom. On the other hand,
- // the solver itself needs a vector
- // that is uniquely split between
- // processors, without any
- // overlap. We therefore create a
- // vector at the beginning of this
- // function that has these
- // properties, use it to solve the
- // linear system, and only assign
- // it to the vector we want at the
- // very end. This last step ensures
- // that all ghost elements are also
- // copied as necessary.
-template <int dim>
-void LaplaceProblem<dim>::solve ()
-{
- PETScWrappers::MPI::Vector
- completely_distributed_solution (mpi_communicator,
- dof_handler.n_dofs(),
- dof_handler.n_locally_owned_dofs());
-
- SolverControl solver_control (dof_handler.n_dofs(), 1e-12);
+ setup_system ();
- PETScWrappers::SolverCG solver(solver_control, mpi_communicator);
+ pcout << " Number of active cells: "
+ << triangulation.n_global_active_cells()
+ << std::endl
+ << " Number of degrees of freedom: "
+ << dof_handler.n_dofs()
+ << std::endl;
- // Ask for a symmetric preconditioner by
- // setting the first parameter in
- // AdditionalData to true.
- PETScWrappers::PreconditionBoomerAMG
- preconditioner(system_matrix,
- PETScWrappers::PreconditionBoomerAMG::AdditionalData(true));
+ assemble_system ();
+ solve ();
- solver.solve (system_matrix, completely_distributed_solution, system_rhs,
- preconditioner);
+ if (Utilities::System::get_n_mpi_processes(mpi_communicator) <= 32)
+ output_results (cycle);
- pcout << " Solved in " << solver_control.last_step()
- << " iterations." << std::endl;
-
- constraints.distribute (completely_distributed_solution);
-
- locally_relevant_solution = completely_distributed_solution;
- locally_relevant_solution.update_ghost_values();
-}
-
-
-
- // @sect4{LaplaceProblem::refine_grid}
-
- // The function that estimates the
- // error and refines the grid is
- // again almost exactly like the one
- // in step-6. The only difference is
- // that the function that flags cells
- // to be refined is now in namespace
- // parallel::distributed::GridRefinement
- // -- a namespace that has functions
- // that can communicate between all
- // involved processors and determine
- // global thresholds to use in
- // deciding which cells to refine and
- // which to coarsen.
- //
- // Note that we didn't have to do
- // anything special about the
- // KellyErrorEstimator class: we just
- // give it a vector with as many
- // elements as the local
- // triangulation has cells (locally
- // owned cells, ghost cells, and
- // artificial ones), but it only
- // fills those entries that
- // correspond to cells that are
- // locally owned.
-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(),
- locally_relevant_solution,
- estimated_error_per_cell);
- parallel::distributed::GridRefinement::
- refine_and_coarsen_fixed_number (triangulation,
- estimated_error_per_cell,
- 0.3, 0.03);
- triangulation.execute_coarsening_and_refinement ();
-}
-
-
-
- // @sect4{LaplaceProblem::output_results}
-
- // Compared to the corresponding
- // function in step-6, the one here
- // is a tad more complicated. There
- // are two reasons: the first one is
- // that we do not just want to output
- // the solution but also for each
- // cell which processor owns it
- // (i.e. which "subdomain" it is
- // in). Secondly, as discussed at
- // length in step-17 and step-18,
- // generating graphical data can be a
- // bottleneck in parallelizing. In
- // step-18, we have moved this step
- // out of the actual computation but
- // shifted it into a separate program
- // that later combined the output
- // from various processors into a
- // single file. But this doesn't
- // scale: if the number of processors
- // is large, this may mean that the
- // step of combining data on a single
- // processor later becomes the
- // longest running part of the
- // program, or it may produce a file
- // that's so large that it can't be
- // visualized any more. We here
- // follow a more sensible approach,
- // namely creating individual files
- // for each MPI process and leaving
- // it to the visualization program to
- // make sense of that.
- //
- // To start, the top of the function
- // looks like always. In addition to
- // attaching the solution vector (the
- // one that has entries for all
- // locally relevant, not only the
- // locally owned, elements), we
- // attach a data vector that stores,
- // for each cell, the subdomain the
- // cell belongs to. This is slightly
- // tricky, because of course not
- // every processor knows about every
- // cell. The vector we attach
- // therefore has an entry for every
- // cell that the current processor
- // has in its mesh (locally owned
- // onces, ghost cells, and artificial
- // cells), but the DataOut class will
- // ignore all entries that correspond
- // to cells that are not owned by the
- // current processor. As a
- // consequence, it doesn't actually
- // matter what values we write into
- // these vector entries: we simply
- // fill the entire vector with the
- // number of the current MPI process
- // (i.e. the subdomain_id of the
- // current process); this correctly
- // sets the values we care for,
- // i.e. the entries that correspond
- // to locally owned cells, while
- // providing the wrong value for all
- // other elements -- but these are
- // then ignored anyway.
-template <int dim>
-void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
-{
- DataOut<dim> data_out;
- data_out.attach_dof_handler (dof_handler);
- data_out.add_data_vector (locally_relevant_solution, "u");
-
- 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 ();
-
- // The next step is to write this
- // data to disk. We choose file
- // names of the form
- // <code>solution-XX-PPPP.vtu</code>
- // where <code>XX</code> indicates
- // the refinement cycle,
- // <code>PPPP</code> refers to the
- // processor number (enough for up
- // to 10,000 processors, though we
- // hope that nobody ever tries to
- // generate this much data -- you
- // would likely overflow all file
- // system quotas), and
- // <code>.vtu</code> indicates the
- // XML-based Visualization Toolkit
- // (VTK) file format.
- const std::string filename = ("solution-" +
- Utilities::int_to_string (cycle, 2) +
- "." +
- Utilities::int_to_string
- (triangulation.locally_owned_subdomain(), 4));
- std::ofstream output ((filename + ".vtu").c_str());
- data_out.write_vtu (output);
-
- // The last step is to write a
- // "master record" that lists for
- // the visualization program the
- // names of the various files that
- // combined represents the
- // graphical data for the entire
- // domain. The
- // DataOutBase::write_pvtu_record
- // does this, and it needs a list
- // of filenames that we create
- // first. Note that only one
- // processor needs to generate this
- // file; we arbitrarily choose
- // processor zero to take over this
- // job.
- if (Utilities::System::get_this_mpi_process(mpi_communicator) == 0)
- {
- std::vector<std::string> filenames;
- for (unsigned int i=0;
- i<Utilities::System::get_n_mpi_processes(mpi_communicator);
- ++i)
- filenames.push_back ("solution-" +
- Utilities::int_to_string (cycle, 2) +
- "." +
- Utilities::int_to_string (i, 4) +
- ".vtu");
-
- std::ofstream master_output ((filename + ".pvtu").c_str());
- data_out.write_pvtu_record (master_output, filenames);
- }
-}
-
-
-
- // @sect4{LaplaceProblem::run}
-
- // The function that controls the
- // overall behavior of the program is
- // again like the one in step-6. The
- // minor difference are the use of
- // <code>pcout</code> instead of
- // <code>std::cout</code> for output
- // to the console (see also step-17)
- // and that we only generate
- // graphical output if at most 32
- // processors are involved. Without
- // this limit, it would be just too
- // easy for people carelessly running
- // this program without reading it
- // first to bring down the cluster
- // interconnect and fill any file
- // system available :-)
- //
- // A functional difference to step-6
- // is the use of a square domain and
- // that we start with a slightly
- // finer mesh (5 global refinement
- // cycles) -- there just isn't much
- // of a point showing a massively
- // %parallel program starting on 4
- // cells (although admittedly the
- // point is only slightly stronger
- // starting on 1024).
-template <int dim>
-void LaplaceProblem<dim>::run ()
-{
- const unsigned int n_cycles = 8;
- for (unsigned int cycle=0; cycle<n_cycles; ++cycle)
- {
- pcout << "Cycle " << cycle << ':' << std::endl;
-
- if (cycle == 0)
- {
- GridGenerator::hyper_cube (triangulation);
- triangulation.refine_global (5);
- }
- else
- refine_grid ();
-
- setup_system ();
-
- pcout << " Number of active cells: "
- << triangulation.n_global_active_cells()
- << std::endl
- << " Number of degrees of freedom: "
- << dof_handler.n_dofs()
- << std::endl;
-
- assemble_system ();
- solve ();
-
- if (Utilities::System::get_n_mpi_processes(mpi_communicator) <= 32)
- output_results (cycle);
-
- pcout << std::endl;
- }
+ pcout << std::endl;
+ }
+ }
}
{
try
{
+ using namespace dealii;
+ using namespace Step40;
+
PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
deallog.depth_console (0);
-
+
{
LaplaceProblem<2> laplace_problem_2d;
laplace_problem_2d.run ();
}
-
+
PetscFinalize();
}
catch (std::exception &exc)
/* $Id$ */
/* */
-/* Copyright (C) 2010 by the deal.II authors */
+/* Copyright (C) 2010, 2011 by the deal.II authors */
/* */
/* This file is subject to QPL and may not be distributed */
/* without copyright and license information. Please refer */
#include <fstream>
-using namespace dealii;
-
- // @sect3{The <code>LaplaceProblem</code> class}
-
- // The class <code>LaplaceProblem</code> is
- // the main class of this problem. As
- // mentioned in the introduction, it is
- // fashioned after the corresponding class in
- // step-3. Correspondingly, the documentation
- // from that tutorial program applies here as
- // well. The only new member variable is the
- // <code>constraints</code> variables that
- // will hold the constraints from the
- // periodic boundary condition. We will
- // initialize it in the
- // <code>make_periodicity_constraints()</code>
- // function which we call from
- // <code>make_grid_and_dofs()</code>.
-class LaplaceProblem
+namespace Step45
{
- public:
- LaplaceProblem ();
- void run ();
-
- private:
- Triangulation<2> triangulation;
-
- FE_Q<2> fe;
- DoFHandler<2> dof_handler;
-
- ConstraintMatrix constraints;
-
- SparsityPattern sparsity_pattern;
- SparseMatrix<double> system_matrix;
- Vector<double> system_rhs;
- Vector<double> solution;
-
- void assemble_system ();
- void output_results ();
- void make_grid_and_dofs ();
- void make_periodicity_constraints ();
- void solve ();
-};
-
-
- // @sect3{The <code>RightHandSide</code> class}
-
- // The following implements the right hand
- // side function discussed in the
- // introduction. Its implementation is
- // obvious given what has been shown in
- // step-4 before:
-class RightHandSide: public Function<2>
-{
- public:
- RightHandSide ();
-
- virtual double value (const Point<2>& p,
- const unsigned int component = 0) const;
-};
-
-
-RightHandSide::RightHandSide ()
- :
- Function<2> ()
-{}
+ using namespace dealii;
+
+ // @sect3{The <code>LaplaceProblem</code> class}
+
+ // The class <code>LaplaceProblem</code> is
+ // the main class of this problem. As
+ // mentioned in the introduction, it is
+ // fashioned after the corresponding class in
+ // step-3. Correspondingly, the documentation
+ // from that tutorial program applies here as
+ // well. The only new member variable is the
+ // <code>constraints</code> variables that
+ // will hold the constraints from the
+ // periodic boundary condition. We will
+ // initialize it in the
+ // <code>make_periodicity_constraints()</code>
+ // function which we call from
+ // <code>make_grid_and_dofs()</code>.
+ class LaplaceProblem
+ {
+ public:
+ LaplaceProblem ();
+ void run ();
+
+ private:
+ Triangulation<2> triangulation;
+
+ FE_Q<2> fe;
+ DoFHandler<2> dof_handler;
+
+ ConstraintMatrix constraints;
+
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+ Vector<double> system_rhs;
+ Vector<double> solution;
+
+ void assemble_system ();
+ void output_results ();
+ void make_grid_and_dofs ();
+ void make_periodicity_constraints ();
+ void solve ();
+ };
+
+
+ // @sect3{The <code>RightHandSide</code> class}
+
+ // The following implements the right hand
+ // side function discussed in the
+ // introduction. Its implementation is
+ // obvious given what has been shown in
+ // step-4 before:
+ class RightHandSide: public Function<2>
+ {
+ public:
+ RightHandSide ();
+
+ virtual double value (const Point<2>& p,
+ const unsigned int component = 0) const;
+ };
+
+
+ RightHandSide::RightHandSide ()
+ :
+ Function<2> ()
+ {}
+
+
+ double
+ RightHandSide::value (const Point<2>&p,
+ const unsigned int) const
+ {
+ return (std::cos (2 * numbers::PI * p(0)) *
+ std::exp (- 2 * p(0)) *
+ std::cos (2 * numbers::PI * p(1)) *
+ std::exp (- 2 * p(1)));
+ }
+
+ // @sect3{Implementation of the <code>LaplaceProblem</code> class}
+
+ // The first part of implementing the main
+ // class is the constructor. It is unchanged
+ // from step-3 and step-4:
+ LaplaceProblem::LaplaceProblem ()
+ :
+ fe (1),
+ dof_handler (triangulation)
+ {}
+
+
+ // @sect4{LaplaceProblem::make_grid_and_dofs}
+
+ // The following is the first function to be
+ // called in <code>run()</code>. It sets up
+ // the mesh and degrees of freedom.
+ //
+ // We start by creating the usual square mesh
+ // and changing the boundary indicator on the
+ // parts of the boundary where we have
+ // Dirichlet boundary conditions (top and
+ // bottom, i.e. faces two and three of the
+ // reference cell as defined by
+ // GeometryInfo), so that we can distinguish
+ // between the parts of the boundary where
+ // periodic and where Dirichlet boundary
+ // conditions hold. We then refine the mesh a
+ // fixed number of times, with child faces
+ // inheriting the boundary indicators
+ // previously set on the coarse mesh from
+ // their parents.
+ void LaplaceProblem::make_grid_and_dofs ()
+ {
+ GridGenerator::hyper_cube (triangulation);
+ triangulation.begin_active ()->face (2)->set_boundary_indicator (1);
+ triangulation.begin_active ()->face (3)->set_boundary_indicator (1);
+ triangulation.refine_global (5);
+
+ // The next step is to distribute the
+ // degrees of freedom and produce a little
+ // bit of graphical output:
+ dof_handler.distribute_dofs (fe);
+ std::cout << "Number of active cells: "
+ << triangulation.n_active_cells ()
+ << std::endl
+ << "Degrees of freedom: " << dof_handler.n_dofs ()
+ << std::endl;
+
+ // Now it is the time for the constraints
+ // that come from the periodicity
+ // constraints. We do this in the
+ // following, separate function, after
+ // clearing any possible prior content from
+ // the constraints object:
+ constraints.clear ();
+ make_periodicity_constraints ();
+
+ // We also incorporate the homogeneous
+ // Dirichlet boundary conditions on the
+ // upper and lower parts of the boundary
+ // (i.e. the ones with boundary indicator
+ // 1) and close the
+ // <code>ConstraintMatrix</code> object:
+ VectorTools::interpolate_boundary_values (dof_handler, 1,
+ ZeroFunction<2> (),
+ constraints);
+ constraints.close ();
+
+ // Then we create the sparsity pattern and
+ // the system matrix and initialize the
+ // solution and right-hand side
+ // vectors. This is again as in step-3 or
+ // step-6, for example:
+ CompressedSparsityPattern c_sparsity_pattern (dof_handler.n_dofs(),
+ dof_handler.n_dofs());
+ DoFTools::make_sparsity_pattern (dof_handler,
+ c_sparsity_pattern,
+ constraints,
+ false);
+ c_sparsity_pattern.compress ();
+ sparsity_pattern.copy_from (c_sparsity_pattern);
+
+ system_matrix.reinit (sparsity_pattern);
+ system_rhs.reinit (dof_handler.n_dofs());
+ solution.reinit (dof_handler.n_dofs());
+ }
+
+
+
+ // @sect4{LaplaceProblem::make_periodicity_constraints}
+
+ // This is the function that provides the new
+ // material of this tutorial program. The
+ // general outline of the algorithm is as
+ // follows: we first loop over all the
+ // degrees of freedom on the right boundary
+ // and record their $y$-locations in a map
+ // together with their global indices. Then
+ // we go along the left boundary, find
+ // matching $y$-locations for each degree of
+ // freedom, and then add constraints that
+ // identify these matched degrees of freedom.
+ //
+ // In this function, we make use of the fact
+ // that we have a scalar element (i.e. the
+ // only valid vector component that can be
+ // passed to DoFAccessor::vertex_dof_index is
+ // zero) and that we have a $Q_1$ element for
+ // which all degrees of freedom live in the
+ // vertices of the cell. Furthermore, we have
+ // assumed that we are in 2d and that meshes
+ // were not refined adaptively — the
+ // latter assumption would imply that there
+ // may be vertices that aren't matched
+ // one-to-one and for which we won't be able
+ // to compute constraints this easily. We
+ // will discuss in the "outlook" part of the
+ // results section below other strategies to
+ // write the current function that can work
+ // in cases like this as well.
+ void LaplaceProblem::make_periodicity_constraints ()
+ {
+ // To start with the actual implementation,
+ // we loop over all active cells and check
+ // whether the cell is located at the right
+ // boundary (i.e. face 1 — the one at
+ // the right end of the cell — is at
+ // the boundary). If that is so, then we
+ // use that for the currently used finite
+ // element, each degree of freedom of the
+ // face is located on one vertex, and store
+ // their $y$-coordinate along with the
+ // global number of this degree of freedom
+ // in the following map:
+ std::map<unsigned int, double> dof_locations;
+
+ for (DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active ();
+ cell != dof_handler.end (); ++cell)
+ if (cell->at_boundary ()
+ &&
+ cell->face(1)->at_boundary ())
+ {
+ dof_locations[cell->face(1)->vertex_dof_index(0, 0)]
+ = cell->face(1)->vertex(0)[1];
+ dof_locations[cell->face(1)->vertex_dof_index(1, 0)]
+ = cell->face(1)->vertex(1)[1];
+ }
+ // Note that in the above block, we add
+ // vertices zero and one of the affected
+ // face to the map. This means that we will
+ // add each vertex twice, once from each of
+ // the two adjacent cells (unless the
+ // vertex is a corner of the domain). Since
+ // the coordinates of the vertex are the
+ // same both times of course, there is no
+ // harm: we replace one value in the map
+ // with itself the second time we visit an
+ // entry.
+ //
+ // The same will be true below where we add
+ // the same constraint twice to the
+ // ConstraintMatrix — again, we will
+ // overwrite the constraint with itself,
+ // and no harm is done.
+
+ // Now we have to find the corresponding
+ // degrees of freedom on the left part of
+ // the boundary. Therefore we loop over all
+ // cells again and choose the ones where
+ // face 0 is at the boundary:
+ for (DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active ();
+ cell != dof_handler.end (); ++cell)
+ if (cell->at_boundary ()
+ &&
+ cell->face (0)->at_boundary ())
+ {
+ // Every degree of freedom on this
+ // face needs to have a corresponding
+ // one on the right side of the face,
+ // and our goal is to add a
+ // constraint for the one on the left
+ // in terms of the one on the
+ // right. To this end we first add a
+ // new line to the constraint matrix
+ // for this one degree of
+ // freedom. Then we identify it with
+ // the corresponding degree of
+ // freedom on the right part of the
+ // boundary by constraining the
+ // degree of freedom on the left with
+ // the one on the right times a
+ // weight of 1.0.
+ //
+ // Consequently, we loop over the two
+ // vertices of each face we find and
+ // then loop over all the
+ // $y$-locations we've previously
+ // recorded to find which degree of
+ // freedom on the right boundary
+ // corresponds to the one we
+ // currently look at. Note that we
+ // have entered these into a map, and
+ // when looping over the iterators
+ // <code>p</code> of this map,
+ // <code>p-@>first</code> corresponds
+ // to the "key" of an entry (the
+ // global number of the degree of
+ // freedom), whereas
+ // <code>p-@>second</code> is the
+ // "value" (the $y$-location we have
+ // entered above).
+ //
+ // We are quite sure here that we
+ // should be finding such a
+ // corresponding degree of
+ // freedom. However, sometimes stuff
+ // happens and so the bottom of the
+ // block contains an assertion that
+ // our assumption was indeed correct
+ // and that a vertex was found.
+ for (unsigned int face_vertex = 0; face_vertex<2; ++face_vertex)
+ {
+ constraints.add_line (cell->face(0)->vertex_dof_index (face_vertex, 0));
+
+ std::map<unsigned int, double>::const_iterator p = dof_locations.begin();
+ for (; p != dof_locations.end(); ++p)
+ if (std::fabs(p->second - cell->face(0)->vertex(face_vertex)[1]) < 1e-8)
+ {
+ constraints.add_entry (cell->face(0)->vertex_dof_index (face_vertex, 0),
+ p->first, 1.0);
+ break;
+ }
+ Assert (p != dof_locations.end(),
+ ExcMessage ("No corresponding degree of freedom was found!"));
+ }
+ }
+ }
+
+
+
+ // @sect4{LaplaceProblem::assemble_system}
+
+ // Assembling the system matrix and the
+ // right-hand side vector is done as in other
+ // tutorials before.
+ //
+ // The only difference here is that we don't
+ // copy elements from local contributions
+ // into the global matrix and later fix up
+ // constrained degrees of freedom, but that
+ // we let the ConstraintMatrix do this job in
+ // one swoop for us using the
+ // ConstraintMatrix::distribute_local_to_global
+ // function(). This was previously already
+ // demonstrated in step-16, step-22, for
+ // example, along with a discussion in the
+ // introduction of step-27.
+ void LaplaceProblem::assemble_system ()
+ {
+ QGauss<2> quadrature_formula(2);
+ FEValues<2> 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();
-double
-RightHandSide::value (const Point<2>&p,
- const unsigned int) const
-{
- return (std::cos (2 * numbers::PI * p(0)) *
- std::exp (- 2 * p(0)) *
- std::cos (2 * numbers::PI * p(1)) *
- std::exp (- 2 * p(1)));
-}
+ FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_rhs (dofs_per_cell);
- // @sect3{Implementation of the <code>LaplaceProblem</code> class}
-
- // The first part of implementing the main
- // class is the constructor. It is unchanged
- // from step-3 and step-4:
-LaplaceProblem::LaplaceProblem ()
- :
- fe (1),
- dof_handler (triangulation)
-{}
-
-
- // @sect4{LaplaceProblem::make_grid_and_dofs}
-
- // The following is the first function to be
- // called in <code>run()</code>. It sets up
- // the mesh and degrees of freedom.
- //
- // We start by creating the usual square mesh
- // and changing the boundary indicator on the
- // parts of the boundary where we have
- // Dirichlet boundary conditions (top and
- // bottom, i.e. faces two and three of the
- // reference cell as defined by
- // GeometryInfo), so that we can distinguish
- // between the parts of the boundary where
- // periodic and where Dirichlet boundary
- // conditions hold. We then refine the mesh a
- // fixed number of times, with child faces
- // inheriting the boundary indicators
- // previously set on the coarse mesh from
- // their parents.
-void LaplaceProblem::make_grid_and_dofs ()
-{
- GridGenerator::hyper_cube (triangulation);
- triangulation.begin_active ()->face (2)->set_boundary_indicator (1);
- triangulation.begin_active ()->face (3)->set_boundary_indicator (1);
- triangulation.refine_global (5);
-
- // The next step is to distribute the
- // degrees of freedom and produce a little
- // bit of graphical output:
- dof_handler.distribute_dofs (fe);
- std::cout << "Number of active cells: "
- << triangulation.n_active_cells ()
- << std::endl
- << "Degrees of freedom: " << dof_handler.n_dofs ()
- << std::endl;
-
- // Now it is the time for the constraints
- // that come from the periodicity
- // constraints. We do this in the
- // following, separate function, after
- // clearing any possible prior content from
- // the constraints object:
- constraints.clear ();
- make_periodicity_constraints ();
-
- // We also incorporate the homogeneous
- // Dirichlet boundary conditions on the
- // upper and lower parts of the boundary
- // (i.e. the ones with boundary indicator
- // 1) and close the
- // <code>ConstraintMatrix</code> object:
- VectorTools::interpolate_boundary_values (dof_handler, 1,
- ZeroFunction<2> (),
- constraints);
- constraints.close ();
-
- // Then we create the sparsity pattern and
- // the system matrix and initialize the
- // solution and right-hand side
- // vectors. This is again as in step-3 or
- // step-6, for example:
- CompressedSparsityPattern c_sparsity_pattern (dof_handler.n_dofs(),
- dof_handler.n_dofs());
- DoFTools::make_sparsity_pattern (dof_handler,
- c_sparsity_pattern,
- constraints,
- false);
- c_sparsity_pattern.compress ();
- sparsity_pattern.copy_from (c_sparsity_pattern);
-
- system_matrix.reinit (sparsity_pattern);
- system_rhs.reinit (dof_handler.n_dofs());
- solution.reinit (dof_handler.n_dofs());
-}
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+ const RightHandSide right_hand_side;
-
- // @sect4{LaplaceProblem::make_periodicity_constraints}
-
- // This is the function that provides the new
- // material of this tutorial program. The
- // general outline of the algorithm is as
- // follows: we first loop over all the
- // degrees of freedom on the right boundary
- // and record their $y$-locations in a map
- // together with their global indices. Then
- // we go along the left boundary, find
- // matching $y$-locations for each degree of
- // freedom, and then add constraints that
- // identify these matched degrees of freedom.
- //
- // In this function, we make use of the fact
- // that we have a scalar element (i.e. the
- // only valid vector component that can be
- // passed to DoFAccessor::vertex_dof_index is
- // zero) and that we have a $Q_1$ element for
- // which all degrees of freedom live in the
- // vertices of the cell. Furthermore, we have
- // assumed that we are in 2d and that meshes
- // were not refined adaptively — the
- // latter assumption would imply that there
- // may be vertices that aren't matched
- // one-to-one and for which we won't be able
- // to compute constraints this easily. We
- // will discuss in the "outlook" part of the
- // results section below other strategies to
- // write the current function that can work
- // in cases like this as well.
-void LaplaceProblem::make_periodicity_constraints ()
-{
- // To start with the actual implementation,
- // we loop over all active cells and check
- // whether the cell is located at the right
- // boundary (i.e. face 1 — the one at
- // the right end of the cell — is at
- // the boundary). If that is so, then we
- // use that for the currently used finite
- // element, each degree of freedom of the
- // face is located on one vertex, and store
- // their $y$-coordinate along with the
- // global number of this degree of freedom
- // in the following map:
- std::map<unsigned int, double> dof_locations;
-
- for (DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active ();
- cell != dof_handler.end (); ++cell)
- if (cell->at_boundary ()
- &&
- cell->face(1)->at_boundary ())
- {
- dof_locations[cell->face(1)->vertex_dof_index(0, 0)]
- = cell->face(1)->vertex(0)[1];
- dof_locations[cell->face(1)->vertex_dof_index(1, 0)]
- = cell->face(1)->vertex(1)[1];
- }
- // Note that in the above block, we add
- // vertices zero and one of the affected
- // face to the map. This means that we will
- // add each vertex twice, once from each of
- // the two adjacent cells (unless the
- // vertex is a corner of the domain). Since
- // the coordinates of the vertex are the
- // same both times of course, there is no
- // harm: we replace one value in the map
- // with itself the second time we visit an
- // entry.
- //
- // The same will be true below where we add
- // the same constraint twice to the
- // ConstraintMatrix — again, we will
- // overwrite the constraint with itself,
- // and no harm is done.
-
- // Now we have to find the corresponding
- // degrees of freedom on the left part of
- // the boundary. Therefore we loop over all
- // cells again and choose the ones where
- // face 0 is at the boundary:
- for (DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active ();
- cell != dof_handler.end (); ++cell)
- if (cell->at_boundary ()
- &&
- cell->face (0)->at_boundary ())
+ DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
{
- // Every degree of freedom on this
- // face needs to have a corresponding
- // one on the right side of the face,
- // and our goal is to add a
- // constraint for the one on the left
- // in terms of the one on the
- // right. To this end we first add a
- // new line to the constraint matrix
- // for this one degree of
- // freedom. Then we identify it with
- // the corresponding degree of
- // freedom on the right part of the
- // boundary by constraining the
- // degree of freedom on the left with
- // the one on the right times a
- // weight of 1.0.
- //
- // Consequently, we loop over the two
- // vertices of each face we find and
- // then loop over all the
- // $y$-locations we've previously
- // recorded to find which degree of
- // freedom on the right boundary
- // corresponds to the one we
- // currently look at. Note that we
- // have entered these into a map, and
- // when looping over the iterators
- // <code>p</code> of this map,
- // <code>p-@>first</code> corresponds
- // to the "key" of an entry (the
- // global number of the degree of
- // freedom), whereas
- // <code>p-@>second</code> is the
- // "value" (the $y$-location we have
- // entered above).
- //
- // We are quite sure here that we
- // should be finding such a
- // corresponding degree of
- // freedom. However, sometimes stuff
- // happens and so the bottom of the
- // block contains an assertion that
- // our assumption was indeed correct
- // and that a vertex was found.
- for (unsigned int face_vertex = 0; face_vertex<2; ++face_vertex)
- {
- constraints.add_line (cell->face(0)->vertex_dof_index (face_vertex, 0));
-
- std::map<unsigned int, double>::const_iterator p = dof_locations.begin();
- for (; p != dof_locations.end(); ++p)
- if (std::fabs(p->second - cell->face(0)->vertex(face_vertex)[1]) < 1e-8)
- {
- constraints.add_entry (cell->face(0)->vertex_dof_index (face_vertex, 0),
- p->first, 1.0);
- break;
- }
- Assert (p != dof_locations.end(),
- ExcMessage ("No corresponding degree of freedom was found!"));
- }
+ fe_values.reinit (cell);
+ cell_matrix = 0;
+ cell_rhs = 0;
+
+ 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) += (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) *
+ right_hand_side.value (fe_values.quadrature_point (q_point)) *
+ 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);
}
-}
-
-
+ }
- // @sect4{LaplaceProblem::assemble_system}
-
- // Assembling the system matrix and the
- // right-hand side vector is done as in other
- // tutorials before.
- //
- // The only difference here is that we don't
- // copy elements from local contributions
- // into the global matrix and later fix up
- // constrained degrees of freedom, but that
- // we let the ConstraintMatrix do this job in
- // one swoop for us using the
- // ConstraintMatrix::distribute_local_to_global
- // function(). This was previously already
- // demonstrated in step-16, step-22, for
- // example, along with a discussion in the
- // introduction of step-27.
-void LaplaceProblem::assemble_system ()
-{
- QGauss<2> quadrature_formula(2);
- FEValues<2> 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();
+ // @sect4{LaplaceProblem::solve}
- FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
- Vector<double> cell_rhs (dofs_per_cell);
+ // To solve the linear system of equations
+ // $Au=b$ we use the CG solver with an
+ // SSOR-preconditioner. This is, again,
+ // copied almost verbatim from step-4, with
+ // the exception of the preconditioner. As in
+ // step-6, we need to make sure that
+ // constrained degrees of freedom get their
+ // correct values after solving by calling
+ // the ConstraintMatrix::distribute function:
+ void LaplaceProblem::solve ()
+ {
+ SolverControl solver_control (dof_handler.n_dofs (), 1e-12);
+ PreconditionSSOR<SparseMatrix<double> > precondition;
- std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+ precondition.initialize (system_matrix);
- const RightHandSide right_hand_side;
+ SolverCG<> cg (solver_control);
- DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active(),
- endc = dof_handler.end();
- for (; cell!=endc; ++cell)
- {
- fe_values.reinit (cell);
- cell_matrix = 0;
- cell_rhs = 0;
-
- 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) += (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) *
- right_hand_side.value (fe_values.quadrature_point (q_point)) *
- 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);
- }
-}
+ cg.solve (system_matrix, solution, system_rhs, precondition);
+ constraints.distribute (solution);
+ }
- // @sect4{LaplaceProblem::solve}
+ // @sect4{LaplaceProblem::output_results}
- // To solve the linear system of equations
- // $Au=b$ we use the CG solver with an
- // SSOR-preconditioner. This is, again,
- // copied almost verbatim from step-4, with
- // the exception of the preconditioner. As in
- // step-6, we need to make sure that
- // constrained degrees of freedom get their
- // correct values after solving by calling
- // the ConstraintMatrix::distribute function:
-void LaplaceProblem::solve ()
-{
- SolverControl solver_control (dof_handler.n_dofs (), 1e-12);
- PreconditionSSOR<SparseMatrix<double> > precondition;
-
- precondition.initialize (system_matrix);
-
- SolverCG<> cg (solver_control);
-
- cg.solve (system_matrix, solution, system_rhs, precondition);
- constraints.distribute (solution);
-}
+ // This is another function copied from
+ // previous tutorial programs. It generates
+ // graphical output in VTK format:
+ void LaplaceProblem::output_results ()
+ {
+ DataOut<2> data_out;
+ data_out.attach_dof_handler (dof_handler);
+ data_out.add_data_vector (solution, "u");
+ data_out.build_patches ();
- // @sect4{LaplaceProblem::output_results}
+ std::ofstream output ("solution.vtk");
- // This is another function copied from
- // previous tutorial programs. It generates
- // graphical output in VTK format:
-void LaplaceProblem::output_results ()
-{
- DataOut<2> data_out;
-
- data_out.attach_dof_handler (dof_handler);
- data_out.add_data_vector (solution, "u");
- data_out.build_patches ();
-
- std::ofstream output ("solution.vtk");
-
- data_out.write_vtk (output);
-}
+ data_out.write_vtk (output);
+ }
- // @sect4{LaplaceProblem::run}
+ // @sect4{LaplaceProblem::run}
- // And another function copied from previous
- // programs:
-void LaplaceProblem::run ()
-{
- make_grid_and_dofs();
- assemble_system ();
- solve ();
- output_results ();
+ // And another function copied from previous
+ // programs:
+ void LaplaceProblem::run ()
+ {
+ make_grid_and_dofs();
+ assemble_system ();
+ solve ();
+ output_results ();
+ }
}
// @sect3{The <code>main</code> function}
{
try
{
+ using namespace dealii;
+ using namespace Step45;
+
deallog.depth_console (0);
LaplaceProblem laplace_problem;
return 1;
}
- catch (...)
+ catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
#include <fstream>
#include <sstream>
-using namespace dealii;
-
- // @sect3{The <code>FluidStructureProblem</code> class template}
-
- // This is the main class. It is, if you
- // want, a combination of step-8 and step-22
- // in that it has member variables that
- // either address the global problem (the
- // Triangulation and hp::DoFHandler objects,
- // as well as the hp::FECollection and
- // various linear algebra objects) or that
- // pertain to either the elasticity or Stokes
- // sub-problems. The general structure of the
- // class, however, is like that of most of
- // the other programs implementing stationary
- // problems.
- //
- // There are a few helper functions
- // (<code>cell_is_in_fluid_domain,
- // cell_is_in_solid_domain</code>) of
- // self-explanatory nature (operating on the
- // symbolic names for the two subdomains that
- // will be used as material_ids for cells
- // belonging to the subdomains, as explained
- // in the introduction) and a few functions
- // (<code>make_grid, set_active_fe_indices,
- // assemble_interface_terms</code>) that have
- // been broken out of other functions that
- // can be found in many of the other tutorial
- // programs and that will be discussed as we
- // get to their implementation.
- //
- // The final set of variables
- // (<code>viscosity, lambda, eta</code>)
- // describes the material properties used for
- // the two physics models.
-template <int dim>
-class FluidStructureProblem
-{
- public:
- FluidStructureProblem (const unsigned int stokes_degree,
- const unsigned int elasticity_degree);
- void run ();
-
- private:
- enum
- {
- fluid_domain_id,
- solid_domain_id
- };
-
- static bool
- cell_is_in_fluid_domain (const typename hp::DoFHandler<dim>::cell_iterator &cell);
-
- static bool
- cell_is_in_solid_domain (const typename hp::DoFHandler<dim>::cell_iterator &cell);
-
-
- void make_grid ();
- void set_active_fe_indices ();
- void setup_dofs ();
- void assemble_system ();
- void assemble_interface_term (const FEFaceValuesBase<dim> &elasticity_fe_face_values,
- const FEFaceValuesBase<dim> &stokes_fe_face_values,
- std::vector<Tensor<1,dim> > &elasticity_phi,
- std::vector<SymmetricTensor<2,dim> > &stokes_phi_grads_u,
- std::vector<double> &stokes_phi_p,
- FullMatrix<double> &local_interface_matrix) const;
- void solve ();
- void output_results (const unsigned int refinement_cycle) const;
- void refine_mesh ();
-
- const unsigned int stokes_degree;
- const unsigned int elasticity_degree;
-
- Triangulation<dim> triangulation;
- FESystem<dim> stokes_fe;
- FESystem<dim> elasticity_fe;
- hp::FECollection<dim> fe_collection;
- hp::DoFHandler<dim> dof_handler;
-
- ConstraintMatrix constraints;
-
- SparsityPattern sparsity_pattern;
- SparseMatrix<double> system_matrix;
-
- Vector<double> solution;
- Vector<double> system_rhs;
-
- const double viscosity;
- const double lambda;
- const double mu;
-};
-
-
- // @sect3{Boundary values and right hand side}
-
- // The following classes do as their names
- // suggest. The boundary values for the
- // velocity are $\mathbf u=(0, \sin(\pi
- // x))^T$ in 2d and $\mathbf u=(0, 0,
- // \sin(\pi x)\sin(\pi y))^T$ in 3d,
- // respectively. The remaining boundary
- // conditions for this problem are all
- // homogenous and have been discussed in the
- // introduction. The right hand side forcing
- // term is zero for both the fluid and the
- // solid.
-template <int dim>
-class StokesBoundaryValues : public Function<dim>
-{
- public:
- StokesBoundaryValues () : Function<dim>(dim+1+dim) {}
-
- virtual double value (const Point<dim> &p,
- const unsigned int component = 0) const;
-
- virtual void vector_value (const Point<dim> &p,
- Vector<double> &value) const;
-};
-
-template <int dim>
-double
-StokesBoundaryValues<dim>::value (const Point<dim> &p,
- const unsigned int component) const
+namespace Step46
{
- Assert (component < this->n_components,
- ExcIndexRange (component, 0, this->n_components));
+ using namespace dealii;
+
+ // @sect3{The <code>FluidStructureProblem</code> class template}
+
+ // This is the main class. It is, if you
+ // want, a combination of step-8 and step-22
+ // in that it has member variables that
+ // either address the global problem (the
+ // Triangulation and hp::DoFHandler objects,
+ // as well as the hp::FECollection and
+ // various linear algebra objects) or that
+ // pertain to either the elasticity or Stokes
+ // sub-problems. The general structure of the
+ // class, however, is like that of most of
+ // the other programs implementing stationary
+ // problems.
+ //
+ // There are a few helper functions
+ // (<code>cell_is_in_fluid_domain,
+ // cell_is_in_solid_domain</code>) of
+ // self-explanatory nature (operating on the
+ // symbolic names for the two subdomains that
+ // will be used as material_ids for cells
+ // belonging to the subdomains, as explained
+ // in the introduction) and a few functions
+ // (<code>make_grid, set_active_fe_indices,
+ // assemble_interface_terms</code>) that have
+ // been broken out of other functions that
+ // can be found in many of the other tutorial
+ // programs and that will be discussed as we
+ // get to their implementation.
+ //
+ // The final set of variables
+ // (<code>viscosity, lambda, eta</code>)
+ // describes the material properties used for
+ // the two physics models.
+ template <int dim>
+ class FluidStructureProblem
+ {
+ public:
+ FluidStructureProblem (const unsigned int stokes_degree,
+ const unsigned int elasticity_degree);
+ void run ();
- if (component == dim-1)
- switch (dim)
+ private:
+ enum
{
- case 2:
- return std::sin(numbers::PI*p[0]);
- case 3:
- return std::sin(numbers::PI*p[0]) * std::sin(numbers::PI*p[1]);
- default:
- Assert (false, ExcNotImplemented());
- }
+ fluid_domain_id,
+ solid_domain_id
+ };
+
+ static bool
+ cell_is_in_fluid_domain (const typename hp::DoFHandler<dim>::cell_iterator &cell);
+
+ static bool
+ cell_is_in_solid_domain (const typename hp::DoFHandler<dim>::cell_iterator &cell);
+
+
+ void make_grid ();
+ void set_active_fe_indices ();
+ void setup_dofs ();
+ void assemble_system ();
+ void assemble_interface_term (const FEFaceValuesBase<dim> &elasticity_fe_face_values,
+ const FEFaceValuesBase<dim> &stokes_fe_face_values,
+ std::vector<Tensor<1,dim> > &elasticity_phi,
+ std::vector<SymmetricTensor<2,dim> > &stokes_phi_grads_u,
+ std::vector<double> &stokes_phi_p,
+ FullMatrix<double> &local_interface_matrix) const;
+ void solve ();
+ void output_results (const unsigned int refinement_cycle) const;
+ void refine_mesh ();
+
+ const unsigned int stokes_degree;
+ const unsigned int elasticity_degree;
+
+ Triangulation<dim> triangulation;
+ FESystem<dim> stokes_fe;
+ FESystem<dim> elasticity_fe;
+ hp::FECollection<dim> fe_collection;
+ hp::DoFHandler<dim> dof_handler;
+
+ ConstraintMatrix constraints;
+
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+
+ Vector<double> solution;
+ Vector<double> system_rhs;
+
+ const double viscosity;
+ const double lambda;
+ const double mu;
+ };
+
+
+ // @sect3{Boundary values and right hand side}
+
+ // The following classes do as their names
+ // suggest. The boundary values for the
+ // velocity are $\mathbf u=(0, \sin(\pi
+ // x))^T$ in 2d and $\mathbf u=(0, 0,
+ // \sin(\pi x)\sin(\pi y))^T$ in 3d,
+ // respectively. The remaining boundary
+ // conditions for this problem are all
+ // homogenous and have been discussed in the
+ // introduction. The right hand side forcing
+ // term is zero for both the fluid and the
+ // solid.
+ template <int dim>
+ class StokesBoundaryValues : public Function<dim>
+ {
+ public:
+ StokesBoundaryValues () : Function<dim>(dim+1+dim) {}
- return 0;
-}
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+ virtual void vector_value (const Point<dim> &p,
+ Vector<double> &value) const;
+ };
-template <int dim>
-void
-StokesBoundaryValues<dim>::vector_value (const Point<dim> &p,
- Vector<double> &values) const
-{
- for (unsigned int c=0; c<this->n_components; ++c)
- values(c) = StokesBoundaryValues<dim>::value (p, c);
-}
+ template <int dim>
+ double
+ StokesBoundaryValues<dim>::value (const Point<dim> &p,
+ const unsigned int component) const
+ {
+ Assert (component < this->n_components,
+ ExcIndexRange (component, 0, this->n_components));
+ if (component == dim-1)
+ switch (dim)
+ {
+ case 2:
+ return std::sin(numbers::PI*p[0]);
+ case 3:
+ return std::sin(numbers::PI*p[0]) * std::sin(numbers::PI*p[1]);
+ default:
+ Assert (false, ExcNotImplemented());
+ }
-template <int dim>
-class RightHandSide : public Function<dim>
-{
- public:
- RightHandSide () : Function<dim>(dim+1) {}
+ return 0;
+ }
- virtual double value (const Point<dim> &p,
- const unsigned int component = 0) const;
- virtual void vector_value (const Point<dim> &p,
- Vector<double> &value) const;
+ template <int dim>
+ void
+ StokesBoundaryValues<dim>::vector_value (const Point<dim> &p,
+ Vector<double> &values) const
+ {
+ for (unsigned int c=0; c<this->n_components; ++c)
+ values(c) = StokesBoundaryValues<dim>::value (p, c);
+ }
-};
-template <int dim>
-double
-RightHandSide<dim>::value (const Point<dim> &/*p*/,
- const unsigned int /*component*/) const
-{
- return 0;
-}
+ template <int dim>
+ class RightHandSide : public Function<dim>
+ {
+ public:
+ RightHandSide () : Function<dim>(dim+1) {}
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
-template <int dim>
-void
-RightHandSide<dim>::vector_value (const Point<dim> &p,
- Vector<double> &values) const
-{
- for (unsigned int c=0; c<this->n_components; ++c)
- values(c) = RightHandSide<dim>::value (p, c);
-}
+ virtual void vector_value (const Point<dim> &p,
+ Vector<double> &value) const;
+ };
- // @sect3{The <code>FluidStructureProblem</code> implementation}
-
- // @sect4{Constructors and helper functions}
-
- // Let's now get to the implementation of the
- // primary class of this program. The first
- // few functions are the constructor and the
- // helper functions that can be used to
- // determine which part of the domain a cell
- // is in. Given the discussion of these
- // topics in the introduction, their
- // implementation is rather obvious. In the
- // constructor, note that we have to
- // construct the hp::FECollection object from
- // the base elements for Stokes and
- // elasticity; using the
- // hp::FECollection::push_back function
- // assigns them spots zero and one in this
- // collection, an order that we have to
- // remember and use consistently in the rest
- // of the program.
-template <int dim>
-FluidStructureProblem<dim>::
-FluidStructureProblem (const unsigned int stokes_degree,
- const unsigned int elasticity_degree)
- :
- stokes_degree (stokes_degree),
- elasticity_degree (elasticity_degree),
- triangulation (Triangulation<dim>::maximum_smoothing),
- stokes_fe (FE_Q<dim>(stokes_degree+1), dim,
- FE_Q<dim>(stokes_degree), 1,
- FE_Nothing<dim>(), dim),
- elasticity_fe (FE_Nothing<dim>(), dim,
- FE_Nothing<dim>(), 1,
- FE_Q<dim>(elasticity_degree), dim),
- dof_handler (triangulation),
- viscosity (2),
- lambda (1),
- mu (1)
-{
- fe_collection.push_back (stokes_fe);
- fe_collection.push_back (elasticity_fe);
-}
+ template <int dim>
+ double
+ RightHandSide<dim>::value (const Point<dim> &/*p*/,
+ const unsigned int /*component*/) const
+ {
+ return 0;
+ }
+ template <int dim>
+ void
+ RightHandSide<dim>::vector_value (const Point<dim> &p,
+ Vector<double> &values) const
+ {
+ for (unsigned int c=0; c<this->n_components; ++c)
+ values(c) = RightHandSide<dim>::value (p, c);
+ }
-template <int dim>
-bool
-FluidStructureProblem<dim>::
-cell_is_in_fluid_domain (const typename hp::DoFHandler<dim>::cell_iterator &cell)
-{
- return (cell->material_id() == fluid_domain_id);
-}
+ // @sect3{The <code>FluidStructureProblem</code> implementation}
+
+ // @sect4{Constructors and helper functions}
+
+ // Let's now get to the implementation of the
+ // primary class of this program. The first
+ // few functions are the constructor and the
+ // helper functions that can be used to
+ // determine which part of the domain a cell
+ // is in. Given the discussion of these
+ // topics in the introduction, their
+ // implementation is rather obvious. In the
+ // constructor, note that we have to
+ // construct the hp::FECollection object from
+ // the base elements for Stokes and
+ // elasticity; using the
+ // hp::FECollection::push_back function
+ // assigns them spots zero and one in this
+ // collection, an order that we have to
+ // remember and use consistently in the rest
+ // of the program.
+ template <int dim>
+ FluidStructureProblem<dim>::
+ FluidStructureProblem (const unsigned int stokes_degree,
+ const unsigned int elasticity_degree)
+ :
+ stokes_degree (stokes_degree),
+ elasticity_degree (elasticity_degree),
+ triangulation (Triangulation<dim>::maximum_smoothing),
+ stokes_fe (FE_Q<dim>(stokes_degree+1), dim,
+ FE_Q<dim>(stokes_degree), 1,
+ FE_Nothing<dim>(), dim),
+ elasticity_fe (FE_Nothing<dim>(), dim,
+ FE_Nothing<dim>(), 1,
+ FE_Q<dim>(elasticity_degree), dim),
+ dof_handler (triangulation),
+ viscosity (2),
+ lambda (1),
+ mu (1)
+ {
+ fe_collection.push_back (stokes_fe);
+ fe_collection.push_back (elasticity_fe);
+ }
-template <int dim>
-bool
-FluidStructureProblem<dim>::
-cell_is_in_solid_domain (const typename hp::DoFHandler<dim>::cell_iterator &cell)
-{
- return (cell->material_id() == solid_domain_id);
-}
- // @sect4{Meshes and assigning subdomains}
-
- // The next pair of functions deals with
- // generating a mesh and making sure all
- // flags that denote subdomains are
- // correct. <code>make_grid</code>, as
- // discussed in the introduction, generates
- // an $8\times 8$ mesh (or an $8\times
- // 8\times 8$ mesh in 3d) to make sure that
- // each coarse mesh cell is completely within
- // one of the subdomains. After generating
- // this mesh, we loop over its boundary and
- // set the boundary indicator to one at the
- // top boundary, the only place where we set
- // nonzero Dirichlet boundary
- // conditions. After this, we loop again over
- // all cells to set the material indicator
- // — used to denote which part of the
- // domain we are in, to either the fluid or
- // solid indicator.
-template <int dim>
-void
-FluidStructureProblem<dim>::make_grid ()
-{
- GridGenerator::subdivided_hyper_cube (triangulation, 8, -1, 1);
-
- for (typename Triangulation<dim>::active_cell_iterator
- cell = triangulation.begin_active();
- cell != triangulation.end(); ++cell)
- for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
- if (cell->face(f)->at_boundary()
- &&
- (cell->face(f)->center()[dim-1] == 1))
- cell->face(f)->set_all_boundary_indicators(1);
-
-
- for (typename Triangulation<dim>::active_cell_iterator
- cell = dof_handler.begin_active();
- cell != dof_handler.end(); ++cell)
- if (((std::fabs(cell->center()[0]) < 0.25)
- &&
- (cell->center()[dim-1] > 0.5))
- ||
- ((std::fabs(cell->center()[0]) >= 0.25)
- &&
- (cell->center()[dim-1] > -0.5)))
- cell->set_material_id (fluid_domain_id);
- else
- cell->set_material_id (solid_domain_id);
-}
+ template <int dim>
+ bool
+ FluidStructureProblem<dim>::
+ cell_is_in_fluid_domain (const typename hp::DoFHandler<dim>::cell_iterator &cell)
+ {
+ return (cell->material_id() == fluid_domain_id);
+ }
- // The second part of this pair of functions
- // determines which finite element to use on
- // each cell. Above we have set the material
- // indicator for each coarse mesh cell, and
- // as mentioned in the introduction, this
- // information is inherited from mother to
- // child cell upon mesh refinement.
- //
- // In other words, whenever we have refined
- // (or created) the mesh, we can rely on the
- // material indicators to be a correct
- // description of which part of the domain a
- // cell is in. We then use this to set the
- // active FE index of the cell to the
- // corresponding element of the
- // hp::FECollection member variable of this
- // class: zero for fluid cells, one for solid
- // cells.
-template <int dim>
-void
-FluidStructureProblem<dim>::set_active_fe_indices ()
-{
- for (typename hp::DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active();
- cell != dof_handler.end(); ++cell)
- {
- if (cell_is_in_fluid_domain(cell))
- cell->set_active_fe_index (0);
- else if (cell_is_in_solid_domain(cell))
- cell->set_active_fe_index (1);
- else
- Assert (false, ExcNotImplemented());
- }
-}
+ template <int dim>
+ bool
+ FluidStructureProblem<dim>::
+ cell_is_in_solid_domain (const typename hp::DoFHandler<dim>::cell_iterator &cell)
+ {
+ return (cell->material_id() == solid_domain_id);
+ }
- // @sect4{<code>FluidStructureProblem::setup_dofs</code>}
-
- // The next step is to setup the data
- // structures for the linear system. To this
- // end, we first have to set the active FE
- // indices with the function immediately
- // above, then distribute degrees of freedom,
- // and then determine constraints on the
- // linear system. The latter includes hanging
- // node constraints as usual, but also the
- // inhomogenous boundary values at the top
- // fluid boundary, and zero boundary values
- // along the perimeter of the solid
- // subdomain.
-template <int dim>
-void
-FluidStructureProblem<dim>::setup_dofs ()
-{
- set_active_fe_indices ();
- dof_handler.distribute_dofs (fe_collection);
+ // @sect4{Meshes and assigning subdomains}
+
+ // The next pair of functions deals with
+ // generating a mesh and making sure all
+ // flags that denote subdomains are
+ // correct. <code>make_grid</code>, as
+ // discussed in the introduction, generates
+ // an $8\times 8$ mesh (or an $8\times
+ // 8\times 8$ mesh in 3d) to make sure that
+ // each coarse mesh cell is completely within
+ // one of the subdomains. After generating
+ // this mesh, we loop over its boundary and
+ // set the boundary indicator to one at the
+ // top boundary, the only place where we set
+ // nonzero Dirichlet boundary
+ // conditions. After this, we loop again over
+ // all cells to set the material indicator
+ // — used to denote which part of the
+ // domain we are in, to either the fluid or
+ // solid indicator.
+ template <int dim>
+ void
+ FluidStructureProblem<dim>::make_grid ()
{
- constraints.clear ();
- DoFTools::make_hanging_node_constraints (dof_handler,
- constraints);
+ GridGenerator::subdivided_hyper_cube (triangulation, 8, -1, 1);
- std::vector<bool> velocity_mask (dim+1+dim, false);
- for (unsigned int d=0; d<dim; ++d)
- velocity_mask[d] = true;
- VectorTools::interpolate_boundary_values (dof_handler,
- 1,
- StokesBoundaryValues<dim>(),
- constraints,
- velocity_mask);
-
- std::vector<bool> elasticity_mask (dim+1+dim, false);
- for (unsigned int d=dim+1; d<dim+1+dim; ++d)
- elasticity_mask[d] = true;
- VectorTools::interpolate_boundary_values (dof_handler,
- 0,
- ZeroFunction<dim>(dim+1+dim),
- constraints,
- elasticity_mask);
+ for (typename Triangulation<dim>::active_cell_iterator
+ cell = triangulation.begin_active();
+ cell != triangulation.end(); ++cell)
+ for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ if (cell->face(f)->at_boundary()
+ &&
+ (cell->face(f)->center()[dim-1] == 1))
+ cell->face(f)->set_all_boundary_indicators(1);
+
+
+ for (typename Triangulation<dim>::active_cell_iterator
+ cell = dof_handler.begin_active();
+ cell != dof_handler.end(); ++cell)
+ if (((std::fabs(cell->center()[0]) < 0.25)
+ &&
+ (cell->center()[dim-1] > 0.5))
+ ||
+ ((std::fabs(cell->center()[0]) >= 0.25)
+ &&
+ (cell->center()[dim-1] > -0.5)))
+ cell->set_material_id (fluid_domain_id);
+ else
+ cell->set_material_id (solid_domain_id);
}
- // There are more constraints we have to
- // handle, though: we have to make sure
- // that the velocity is zero at the
- // interface between fluid and solid. The
- // following piece of code was already
- // presented in the introduction:
+
+ // The second part of this pair of functions
+ // determines which finite element to use on
+ // each cell. Above we have set the material
+ // indicator for each coarse mesh cell, and
+ // as mentioned in the introduction, this
+ // information is inherited from mother to
+ // child cell upon mesh refinement.
+ //
+ // In other words, whenever we have refined
+ // (or created) the mesh, we can rely on the
+ // material indicators to be a correct
+ // description of which part of the domain a
+ // cell is in. We then use this to set the
+ // active FE index of the cell to the
+ // corresponding element of the
+ // hp::FECollection member variable of this
+ // class: zero for fluid cells, one for solid
+ // cells.
+ template <int dim>
+ void
+ FluidStructureProblem<dim>::set_active_fe_indices ()
{
- std::vector<unsigned int> local_face_dof_indices (stokes_fe.dofs_per_face);
for (typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active();
cell != dof_handler.end(); ++cell)
- if (cell_is_in_fluid_domain (cell))
- for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
- if (!cell->at_boundary(f))
- {
- bool face_is_on_interface = false;
-
- if ((cell->neighbor(f)->has_children() == false)
- &&
- (cell_is_in_solid_domain (cell->neighbor(f))))
- face_is_on_interface = true;
- else if (cell->neighbor(f)->has_children() == true)
- {
- for (unsigned int sf=0; sf<cell->face(f)->n_children(); ++sf)
- if (cell_is_in_solid_domain (cell->neighbor_child_on_subface
- (f, sf)))
- {
- face_is_on_interface = true;
- break;
- }
- }
-
- if (face_is_on_interface)
- {
- cell->face(f)->get_dof_indices (local_face_dof_indices, 0);
- for (unsigned int i=0; i<local_face_dof_indices.size(); ++i)
- if (stokes_fe.face_system_to_component_index(i).first < dim)
- constraints.add_line (local_face_dof_indices[i]);
- }
- }
+ {
+ if (cell_is_in_fluid_domain(cell))
+ cell->set_active_fe_index (0);
+ else if (cell_is_in_solid_domain(cell))
+ cell->set_active_fe_index (1);
+ else
+ Assert (false, ExcNotImplemented());
+ }
}
- // At the end of all this, we can declare
- // to the constraints object that we now
- // have all constraints ready to go and
- // that the object can rebuild its internal
- // data structures for better efficiency:
- constraints.close ();
-
- std::cout << " Number of active cells: "
- << triangulation.n_active_cells()
- << std::endl
- << " Number of degrees of freedom: "
- << dof_handler.n_dofs()
- << std::endl;
-
- // In the rest of this function we create a
- // sparsity pattern as discussed
- // extensively in the introduction, and use
- // it to initialize the matrix; then also
- // set vectors to their correct sizes:
+
+ // @sect4{<code>FluidStructureProblem::setup_dofs</code>}
+
+ // The next step is to setup the data
+ // structures for the linear system. To this
+ // end, we first have to set the active FE
+ // indices with the function immediately
+ // above, then distribute degrees of freedom,
+ // and then determine constraints on the
+ // linear system. The latter includes hanging
+ // node constraints as usual, but also the
+ // inhomogenous boundary values at the top
+ // fluid boundary, and zero boundary values
+ // along the perimeter of the solid
+ // subdomain.
+ template <int dim>
+ void
+ FluidStructureProblem<dim>::setup_dofs ()
{
- CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
- dof_handler.n_dofs());
+ set_active_fe_indices ();
+ dof_handler.distribute_dofs (fe_collection);
+
+ {
+ constraints.clear ();
+ DoFTools::make_hanging_node_constraints (dof_handler,
+ constraints);
+
+ std::vector<bool> velocity_mask (dim+1+dim, false);
+ for (unsigned int d=0; d<dim; ++d)
+ velocity_mask[d] = true;
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 1,
+ StokesBoundaryValues<dim>(),
+ constraints,
+ velocity_mask);
+
+ std::vector<bool> elasticity_mask (dim+1+dim, false);
+ for (unsigned int d=dim+1; d<dim+1+dim; ++d)
+ elasticity_mask[d] = true;
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 0,
+ ZeroFunction<dim>(dim+1+dim),
+ constraints,
+ elasticity_mask);
+ }
- Table<2,DoFTools::Coupling> cell_coupling (fe_collection.n_components(),
- fe_collection.n_components());
- Table<2,DoFTools::Coupling> face_coupling (fe_collection.n_components(),
- fe_collection.n_components());
+ // There are more constraints we have to
+ // handle, though: we have to make sure
+ // that the velocity is zero at the
+ // interface between fluid and solid. The
+ // following piece of code was already
+ // presented in the introduction:
+ {
+ std::vector<unsigned int> local_face_dof_indices (stokes_fe.dofs_per_face);
+ for (typename hp::DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active();
+ cell != dof_handler.end(); ++cell)
+ if (cell_is_in_fluid_domain (cell))
+ for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ if (!cell->at_boundary(f))
+ {
+ bool face_is_on_interface = false;
+
+ if ((cell->neighbor(f)->has_children() == false)
+ &&
+ (cell_is_in_solid_domain (cell->neighbor(f))))
+ face_is_on_interface = true;
+ else if (cell->neighbor(f)->has_children() == true)
+ {
+ for (unsigned int sf=0; sf<cell->face(f)->n_children(); ++sf)
+ if (cell_is_in_solid_domain (cell->neighbor_child_on_subface
+ (f, sf)))
+ {
+ face_is_on_interface = true;
+ break;
+ }
+ }
- for (unsigned int c=0; c<fe_collection.n_components(); ++c)
- for (unsigned int d=0; d<fe_collection.n_components(); ++d)
- {
- if (((c<dim+1) && (d<dim+1)
- && !((c==dim) && (d==dim)))
- ||
- ((c>=dim+1) && (d>=dim+1)))
- cell_coupling[c][d] = DoFTools::always;
-
- if ((c>=dim+1) && (d<dim+1))
- face_coupling[c][d] = DoFTools::always;
- }
-
- DoFTools::make_flux_sparsity_pattern (dof_handler, csp,
- cell_coupling, face_coupling);
- constraints.condense (csp);
- sparsity_pattern.copy_from (csp);
+ if (face_is_on_interface)
+ {
+ cell->face(f)->get_dof_indices (local_face_dof_indices, 0);
+ for (unsigned int i=0; i<local_face_dof_indices.size(); ++i)
+ if (stokes_fe.face_system_to_component_index(i).first < dim)
+ constraints.add_line (local_face_dof_indices[i]);
+ }
+ }
+ }
+
+ // At the end of all this, we can declare
+ // to the constraints object that we now
+ // have all constraints ready to go and
+ // that the object can rebuild its internal
+ // data structures for better efficiency:
+ constraints.close ();
+
+ std::cout << " Number of active cells: "
+ << triangulation.n_active_cells()
+ << std::endl
+ << " Number of degrees of freedom: "
+ << dof_handler.n_dofs()
+ << std::endl;
+
+ // In the rest of this function we create a
+ // sparsity pattern as discussed
+ // extensively in the introduction, and use
+ // it to initialize the matrix; then also
+ // set vectors to their correct sizes:
+ {
+ CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
+ dof_handler.n_dofs());
+
+ Table<2,DoFTools::Coupling> cell_coupling (fe_collection.n_components(),
+ fe_collection.n_components());
+ Table<2,DoFTools::Coupling> face_coupling (fe_collection.n_components(),
+ fe_collection.n_components());
+
+ for (unsigned int c=0; c<fe_collection.n_components(); ++c)
+ for (unsigned int d=0; d<fe_collection.n_components(); ++d)
+ {
+ if (((c<dim+1) && (d<dim+1)
+ && !((c==dim) && (d==dim)))
+ ||
+ ((c>=dim+1) && (d>=dim+1)))
+ cell_coupling[c][d] = DoFTools::always;
+
+ if ((c>=dim+1) && (d<dim+1))
+ face_coupling[c][d] = DoFTools::always;
+ }
+
+ DoFTools::make_flux_sparsity_pattern (dof_handler, csp,
+ cell_coupling, face_coupling);
+ constraints.condense (csp);
+ sparsity_pattern.copy_from (csp);
+ }
+
+ system_matrix.reinit (sparsity_pattern);
+
+ solution.reinit (dof_handler.n_dofs());
+ system_rhs.reinit (dof_handler.n_dofs());
}
- system_matrix.reinit (sparsity_pattern);
- solution.reinit (dof_handler.n_dofs());
- system_rhs.reinit (dof_handler.n_dofs());
-}
+ // @sect4{<code>FluidStructureProblem::assemble_system</code>}
+ // Following is the central function of this
+ // program: the one that assembles the linear
+ // system. It has a long section of setting
+ // up auxiliary functions at the beginning:
+ // from creating the quadrature formulas and
+ // setting up the FEValues, FEFaceValues and
+ // FESubfaceValues objects necessary to
+ // integrate the cell terms as well as the
+ // interface terms for the case where cells
+ // along the interface come together at same
+ // size or with differing levels of
+ // refinement...
+ template <int dim>
+ void FluidStructureProblem<dim>::assemble_system ()
+ {
+ system_matrix=0;
+ system_rhs=0;
+
+ const QGauss<dim> stokes_quadrature(stokes_degree+2);
+ const QGauss<dim> elasticity_quadrature(elasticity_degree+2);
+
+ hp::QCollection<dim> q_collection;
+ q_collection.push_back (stokes_quadrature);
+ q_collection.push_back (elasticity_quadrature);
+
+ hp::FEValues<dim> hp_fe_values (fe_collection, q_collection,
+ update_values |
+ update_quadrature_points |
+ update_JxW_values |
+ update_gradients);
+
+ const QGauss<dim-1> common_face_quadrature(std::max (stokes_degree+2,
+ elasticity_degree+2));
+
+ FEFaceValues<dim> stokes_fe_face_values (stokes_fe,
+ common_face_quadrature,
+ update_JxW_values |
+ update_normal_vectors |
+ update_gradients);
+ FEFaceValues<dim> elasticity_fe_face_values (elasticity_fe,
+ common_face_quadrature,
+ update_values);
+ FESubfaceValues<dim> stokes_fe_subface_values (stokes_fe,
+ common_face_quadrature,
+ update_JxW_values |
+ update_normal_vectors |
+ update_gradients);
+ FESubfaceValues<dim> elasticity_fe_subface_values (elasticity_fe,
+ common_face_quadrature,
+ update_values);
+
+ // ...to objects that are needed to
+ // describe the local contributions to the
+ // global linear system...
+ const unsigned int stokes_dofs_per_cell = stokes_fe.dofs_per_cell;
+ const unsigned int elasticity_dofs_per_cell = elasticity_fe.dofs_per_cell;
+
+ FullMatrix<double> local_matrix;
+ FullMatrix<double> local_interface_matrix (elasticity_dofs_per_cell,
+ stokes_dofs_per_cell);
+ Vector<double> local_rhs;
+
+ std::vector<unsigned int> local_dof_indices;
+ std::vector<unsigned int> neighbor_dof_indices (stokes_dofs_per_cell);
+
+ const RightHandSide<dim> right_hand_side;
+
+ // ...to variables that allow us to extract
+ // certain components of the shape
+ // functions and cache their values rather
+ // than having to recompute them at every
+ // quadrature point:
+ const FEValuesExtractors::Vector velocities (0);
+ const FEValuesExtractors::Scalar pressure (dim);
+ const FEValuesExtractors::Vector displacements (dim+1);
+
+ std::vector<SymmetricTensor<2,dim> > stokes_phi_grads_u (stokes_dofs_per_cell);
+ std::vector<double> stokes_div_phi_u (stokes_dofs_per_cell);
+ std::vector<double> stokes_phi_p (stokes_dofs_per_cell);
+
+ std::vector<Tensor<2,dim> > elasticity_phi_grad (elasticity_dofs_per_cell);
+ std::vector<double> elasticity_phi_div (elasticity_dofs_per_cell);
+ std::vector<Tensor<1,dim> > elasticity_phi (elasticity_dofs_per_cell);
+
+ // Then comes the main loop over all cells
+ // and, as in step-27, the initialization
+ // of the hp::FEValues object for the
+ // current cell and the extraction of a
+ // FEValues object that is appropriate for
+ // the current cell:
+ 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();
+
+ local_matrix.reinit (cell->get_fe().dofs_per_cell,
+ cell->get_fe().dofs_per_cell);
+ local_rhs.reinit (cell->get_fe().dofs_per_cell);
+
+ // With all of this done, we continue
+ // to assemble the cell terms for cells
+ // that are part of the Stokes and
+ // elastic regions. While we could in
+ // principle do this in one formula, in
+ // effect implementing the one bilinear
+ // form stated in the introduction, we
+ // realize that our finite element
+ // spaces are chosen in such a way that
+ // on each cell, one set of variables
+ // (either velocities and pressure, or
+ // displacements) are always zero, and
+ // consequently a more efficient way of
+ // computing local integrals is to do
+ // only what's necessary based on an
+ // <code>if</code> clause that tests
+ // which part of the domain we are in.
+ //
+ // The actual computation of the local
+ // matrix is the same as in step-22 as
+ // well as that given in the @ref
+ // vector_valued documentation module
+ // for the elasticity equations:
+ if (cell_is_in_fluid_domain (cell))
+ {
+ const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
+ Assert (dofs_per_cell == stokes_dofs_per_cell,
+ ExcInternalError());
- // @sect4{<code>FluidStructureProblem::assemble_system</code>}
+ for (unsigned int q=0; q<fe_values.n_quadrature_points; ++q)
+ {
+ for (unsigned int k=0; k<dofs_per_cell; ++k)
+ {
+ stokes_phi_grads_u[k] = fe_values[velocities].symmetric_gradient (k, q);
+ stokes_div_phi_u[k] = fe_values[velocities].divergence (k, q);
+ stokes_phi_p[k] = fe_values[pressure].value (k, q);
+ }
- // Following is the central function of this
- // program: the one that assembles the linear
- // system. It has a long section of setting
- // up auxiliary functions at the beginning:
- // from creating the quadrature formulas and
- // setting up the FEValues, FEFaceValues and
- // FESubfaceValues objects necessary to
- // integrate the cell terms as well as the
- // interface terms for the case where cells
- // along the interface come together at same
- // size or with differing levels of
- // refinement...
-template <int dim>
-void FluidStructureProblem<dim>::assemble_system ()
-{
- system_matrix=0;
- system_rhs=0;
-
- const QGauss<dim> stokes_quadrature(stokes_degree+2);
- const QGauss<dim> elasticity_quadrature(elasticity_degree+2);
-
- hp::QCollection<dim> q_collection;
- q_collection.push_back (stokes_quadrature);
- q_collection.push_back (elasticity_quadrature);
-
- hp::FEValues<dim> hp_fe_values (fe_collection, q_collection,
- update_values |
- update_quadrature_points |
- update_JxW_values |
- update_gradients);
-
- const QGauss<dim-1> common_face_quadrature(std::max (stokes_degree+2,
- elasticity_degree+2));
-
- FEFaceValues<dim> stokes_fe_face_values (stokes_fe,
- common_face_quadrature,
- update_JxW_values |
- update_normal_vectors |
- update_gradients);
- FEFaceValues<dim> elasticity_fe_face_values (elasticity_fe,
- common_face_quadrature,
- update_values);
- FESubfaceValues<dim> stokes_fe_subface_values (stokes_fe,
- common_face_quadrature,
- update_JxW_values |
- update_normal_vectors |
- update_gradients);
- FESubfaceValues<dim> elasticity_fe_subface_values (elasticity_fe,
- common_face_quadrature,
- update_values);
-
- // ...to objects that are needed to
- // describe the local contributions to the
- // global linear system...
- const unsigned int stokes_dofs_per_cell = stokes_fe.dofs_per_cell;
- const unsigned int elasticity_dofs_per_cell = elasticity_fe.dofs_per_cell;
-
- FullMatrix<double> local_matrix;
- FullMatrix<double> local_interface_matrix (elasticity_dofs_per_cell,
- stokes_dofs_per_cell);
- Vector<double> local_rhs;
-
- std::vector<unsigned int> local_dof_indices;
- std::vector<unsigned int> neighbor_dof_indices (stokes_dofs_per_cell);
-
- const RightHandSide<dim> right_hand_side;
-
- // ...to variables that allow us to extract
- // certain components of the shape
- // functions and cache their values rather
- // than having to recompute them at every
- // quadrature point:
- const FEValuesExtractors::Vector velocities (0);
- const FEValuesExtractors::Scalar pressure (dim);
- const FEValuesExtractors::Vector displacements (dim+1);
-
- std::vector<SymmetricTensor<2,dim> > stokes_phi_grads_u (stokes_dofs_per_cell);
- std::vector<double> stokes_div_phi_u (stokes_dofs_per_cell);
- std::vector<double> stokes_phi_p (stokes_dofs_per_cell);
-
- std::vector<Tensor<2,dim> > elasticity_phi_grad (elasticity_dofs_per_cell);
- std::vector<double> elasticity_phi_div (elasticity_dofs_per_cell);
- std::vector<Tensor<1,dim> > elasticity_phi (elasticity_dofs_per_cell);
-
- // Then comes the main loop over all cells
- // and, as in step-27, the initialization
- // of the hp::FEValues object for the
- // current cell and the extraction of a
- // FEValues object that is appropriate for
- // the current cell:
- 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();
-
- local_matrix.reinit (cell->get_fe().dofs_per_cell,
- cell->get_fe().dofs_per_cell);
- local_rhs.reinit (cell->get_fe().dofs_per_cell);
-
- // With all of this done, we continue
- // to assemble the cell terms for cells
- // that are part of the Stokes and
- // elastic regions. While we could in
- // principle do this in one formula, in
- // effect implementing the one bilinear
- // form stated in the introduction, we
- // realize that our finite element
- // spaces are chosen in such a way that
- // on each cell, one set of variables
- // (either velocities and pressure, or
- // displacements) are always zero, and
- // consequently a more efficient way of
- // computing local integrals is to do
- // only what's necessary based on an
- // <code>if</code> clause that tests
- // which part of the domain we are in.
- //
- // The actual computation of the local
- // matrix is the same as in step-22 as
- // well as that given in the @ref
- // vector_valued documentation module
- // for the elasticity equations:
- if (cell_is_in_fluid_domain (cell))
- {
- const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
- Assert (dofs_per_cell == stokes_dofs_per_cell,
- ExcInternalError());
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ local_matrix(i,j) += (2 * viscosity * stokes_phi_grads_u[i] * stokes_phi_grads_u[j]
+ - stokes_div_phi_u[i] * stokes_phi_p[j]
+ - stokes_phi_p[i] * stokes_div_phi_u[j])
+ * fe_values.JxW(q);
+ }
+ }
+ else
+ {
+ const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
+ Assert (dofs_per_cell == elasticity_dofs_per_cell,
+ ExcInternalError());
- for (unsigned int q=0; q<fe_values.n_quadrature_points; ++q)
- {
- for (unsigned int k=0; k<dofs_per_cell; ++k)
- {
- stokes_phi_grads_u[k] = fe_values[velocities].symmetric_gradient (k, q);
- stokes_div_phi_u[k] = fe_values[velocities].divergence (k, q);
- stokes_phi_p[k] = fe_values[pressure].value (k, q);
- }
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- local_matrix(i,j) += (2 * viscosity * stokes_phi_grads_u[i] * stokes_phi_grads_u[j]
- - stokes_div_phi_u[i] * stokes_phi_p[j]
- - stokes_phi_p[i] * stokes_div_phi_u[j])
- * fe_values.JxW(q);
- }
- }
- else
- {
- const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
- Assert (dofs_per_cell == elasticity_dofs_per_cell,
- ExcInternalError());
+ for (unsigned int q=0; q<fe_values.n_quadrature_points; ++q)
+ {
+ for (unsigned int k=0; k<dofs_per_cell; ++k)
+ {
+ elasticity_phi_grad[k] = fe_values[displacements].gradient (k, q);
+ elasticity_phi_div[k] = fe_values[displacements].divergence (k, q);
+ }
- for (unsigned int q=0; q<fe_values.n_quadrature_points; ++q)
- {
- for (unsigned int k=0; k<dofs_per_cell; ++k)
- {
- elasticity_phi_grad[k] = fe_values[displacements].gradient (k, q);
- elasticity_phi_div[k] = fe_values[displacements].divergence (k, q);
- }
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ local_matrix(i,j)
+ += (lambda *
+ elasticity_phi_div[i] * elasticity_phi_div[j]
+ +
+ mu *
+ scalar_product(elasticity_phi_grad[i], elasticity_phi_grad[j])
+ +
+ mu *
+ scalar_product(elasticity_phi_grad[i], transpose(elasticity_phi_grad[j]))
+ )
+ *
+ fe_values.JxW(q);
+ }
+ }
+ }
+
+ // Once we have the contributions from
+ // cell integrals, we copy them into
+ // the global matrix (taking care of
+ // constraints right away, through the
+ // ConstraintMatrix::distribute_local_to_global
+ // function). Note that we have not
+ // written anything into the
+ // <code>local_rhs</code> variable,
+ // though we still need to pass it
+ // along since the elimination of
+ // nonzero boundary values requires the
+ // modification of local and
+ // consequently also global right hand
+ // side values:
+ local_dof_indices.resize (cell->get_fe().dofs_per_cell);
+ cell->get_dof_indices (local_dof_indices);
+ constraints.distribute_local_to_global (local_matrix, local_rhs,
+ local_dof_indices,
+ system_matrix, system_rhs);
+
+ // The more interesting part of this
+ // function is where we see about face
+ // terms along the interface between
+ // the two subdomains. To this end, we
+ // first have to make sure that we only
+ // assemble them once even though a
+ // loop over all faces of all cells
+ // would encounter each part of the
+ // interface twice. We arbitrarily make
+ // the decision that we will only
+ // evaluate interface terms if the
+ // current cell is part of the solid
+ // subdomain and if, consequently, a
+ // face is not at the boundary and the
+ // potential neighbor behind it is part
+ // of the fluid domain. Let's start
+ // with these conditions:
+ if (cell_is_in_solid_domain (cell))
+ for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ if (cell->at_boundary(f) == false)
+ {
+ // At this point we know that
+ // the current cell is a
+ // candidate for integration
+ // and that a neighbor behind
+ // face <code>f</code>
+ // exists. There are now three
+ // possibilities:
+ //
+ // - The neighbor is at the
+ // same refinement level and
+ // has no children.
+ // - The neighbor has children.
+ // - The neighbor is coarser.
+ //
+ // In all three cases, we are
+ // only interested in it if it
+ // is part of the fluid
+ // subdomain. So let us start
+ // with the first and simplest
+ // case: if the neighbor is at
+ // the same level, has no
+ // children, and is a fluid
+ // cell, then the two cells
+ // share a boundary that is
+ // part of the interface along
+ // which we want to integrate
+ // interface terms. All we have
+ // to do is initialize two
+ // FEFaceValues object with the
+ // current face and the face of
+ // the neighboring cell (note
+ // how we find out which face
+ // of the neighboring cell
+ // borders on the current cell)
+ // and pass things off to the
+ // function that evaluates the
+ // interface terms (the third
+ // through fifth arguments to
+ // this function provide it
+ // with scratch arrays). The
+ // result is then again copied
+ // into the global matrix,
+ // using a function that knows
+ // that the DoF indices of rows
+ // and columns of the local
+ // matrix result from different
+ // cells:
+ if ((cell->neighbor(f)->level() == cell->level())
+ &&
+ (cell->neighbor(f)->has_children() == false)
+ &&
+ cell_is_in_fluid_domain (cell->neighbor(f)))
{
- local_matrix(i,j)
- += (lambda *
- elasticity_phi_div[i] * elasticity_phi_div[j]
- +
- mu *
- scalar_product(elasticity_phi_grad[i], elasticity_phi_grad[j])
- +
- mu *
- scalar_product(elasticity_phi_grad[i], transpose(elasticity_phi_grad[j]))
- )
- *
- fe_values.JxW(q);
+ elasticity_fe_face_values.reinit (cell, f);
+ stokes_fe_face_values.reinit (cell->neighbor(f),
+ cell->neighbor_of_neighbor(f));
+
+ assemble_interface_term (elasticity_fe_face_values, stokes_fe_face_values,
+ elasticity_phi, stokes_phi_grads_u, stokes_phi_p,
+ local_interface_matrix);
+
+ cell->neighbor(f)->get_dof_indices (neighbor_dof_indices);
+ constraints.distribute_local_to_global(local_interface_matrix,
+ local_dof_indices,
+ neighbor_dof_indices,
+ system_matrix);
}
- }
- }
- // Once we have the contributions from
- // cell integrals, we copy them into
- // the global matrix (taking care of
- // constraints right away, through the
- // ConstraintMatrix::distribute_local_to_global
- // function). Note that we have not
- // written anything into the
- // <code>local_rhs</code> variable,
- // though we still need to pass it
- // along since the elimination of
- // nonzero boundary values requires the
- // modification of local and
- // consequently also global right hand
- // side values:
- local_dof_indices.resize (cell->get_fe().dofs_per_cell);
- cell->get_dof_indices (local_dof_indices);
- constraints.distribute_local_to_global (local_matrix, local_rhs,
- local_dof_indices,
- system_matrix, system_rhs);
-
- // The more interesting part of this
- // function is where we see about face
- // terms along the interface between
- // the two subdomains. To this end, we
- // first have to make sure that we only
- // assemble them once even though a
- // loop over all faces of all cells
- // would encounter each part of the
- // interface twice. We arbitrarily make
- // the decision that we will only
- // evaluate interface terms if the
- // current cell is part of the solid
- // subdomain and if, consequently, a
- // face is not at the boundary and the
- // potential neighbor behind it is part
- // of the fluid domain. Let's start
- // with these conditions:
- if (cell_is_in_solid_domain (cell))
- for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
- if (cell->at_boundary(f) == false)
- {
- // At this point we know that
- // the current cell is a
- // candidate for integration
- // and that a neighbor behind
- // face <code>f</code>
- // exists. There are now three
- // possibilities:
- //
- // - The neighbor is at the
- // same refinement level and
- // has no children.
- // - The neighbor has children.
- // - The neighbor is coarser.
- //
- // In all three cases, we are
- // only interested in it if it
- // is part of the fluid
- // subdomain. So let us start
- // with the first and simplest
- // case: if the neighbor is at
- // the same level, has no
- // children, and is a fluid
- // cell, then the two cells
- // share a boundary that is
- // part of the interface along
- // which we want to integrate
- // interface terms. All we have
- // to do is initialize two
- // FEFaceValues object with the
- // current face and the face of
- // the neighboring cell (note
- // how we find out which face
- // of the neighboring cell
- // borders on the current cell)
- // and pass things off to the
- // function that evaluates the
- // interface terms (the third
- // through fifth arguments to
- // this function provide it
- // with scratch arrays). The
- // result is then again copied
- // into the global matrix,
- // using a function that knows
- // that the DoF indices of rows
- // and columns of the local
- // matrix result from different
- // cells:
- if ((cell->neighbor(f)->level() == cell->level())
- &&
- (cell->neighbor(f)->has_children() == false)
- &&
- cell_is_in_fluid_domain (cell->neighbor(f)))
- {
- elasticity_fe_face_values.reinit (cell, f);
- stokes_fe_face_values.reinit (cell->neighbor(f),
- cell->neighbor_of_neighbor(f));
-
- assemble_interface_term (elasticity_fe_face_values, stokes_fe_face_values,
- elasticity_phi, stokes_phi_grads_u, stokes_phi_p,
- local_interface_matrix);
-
- cell->neighbor(f)->get_dof_indices (neighbor_dof_indices);
- constraints.distribute_local_to_global(local_interface_matrix,
- local_dof_indices,
- neighbor_dof_indices,
- system_matrix);
- }
-
- // The second case is if the
- // neighbor has further
- // children. In that case, we
- // have to loop over all the
- // children of the neighbor to
- // see if they are part of the
- // fluid subdomain. If they
- // are, then we integrate over
- // the common interface, which
- // is a face for the neighbor
- // and a subface of the current
- // cell, requiring us to use an
- // FEFaceValues for the
- // neighbor and an
- // FESubfaceValues for the
- // current cell:
- else if ((cell->neighbor(f)->level() == cell->level())
- &&
- (cell->neighbor(f)->has_children() == true))
- {
- for (unsigned int subface=0;
- subface<cell->face(f)->n_children();
- ++subface)
- if (cell_is_in_fluid_domain (cell->neighbor_child_on_subface
- (f, subface)))
- {
- elasticity_fe_subface_values.reinit (cell,
- f,
- subface);
- stokes_fe_face_values.reinit (cell->neighbor_child_on_subface (f, subface),
- cell->neighbor_of_neighbor(f));
-
- assemble_interface_term (elasticity_fe_subface_values,
- stokes_fe_face_values,
- elasticity_phi,
- stokes_phi_grads_u, stokes_phi_p,
- local_interface_matrix);
-
- cell->neighbor_child_on_subface (f, subface)
- ->get_dof_indices (neighbor_dof_indices);
- constraints.distribute_local_to_global(local_interface_matrix,
- local_dof_indices,
- neighbor_dof_indices,
- system_matrix);
- }
- }
-
- // The last option is that the
- // neighbor is coarser. In that
- // case we have to use an
- // FESubfaceValues object for
- // the neighbor and a
- // FEFaceValues for the current
- // cell; the rest is the same
- // as before:
- else if (cell->neighbor_is_coarser(f)
- &&
- cell_is_in_fluid_domain(cell->neighbor(f)))
- {
- elasticity_fe_face_values.reinit (cell, f);
- stokes_fe_subface_values.reinit (cell->neighbor(f),
- cell->neighbor_of_coarser_neighbor(f).first,
- cell->neighbor_of_coarser_neighbor(f).second);
-
- assemble_interface_term (elasticity_fe_face_values,
- stokes_fe_subface_values,
- elasticity_phi,
- stokes_phi_grads_u, stokes_phi_p,
- local_interface_matrix);
-
- cell->neighbor(f)->get_dof_indices (neighbor_dof_indices);
- constraints.distribute_local_to_global(local_interface_matrix,
- local_dof_indices,
- neighbor_dof_indices,
- system_matrix);
-
- }
- }
- }
-}
+ // The second case is if the
+ // neighbor has further
+ // children. In that case, we
+ // have to loop over all the
+ // children of the neighbor to
+ // see if they are part of the
+ // fluid subdomain. If they
+ // are, then we integrate over
+ // the common interface, which
+ // is a face for the neighbor
+ // and a subface of the current
+ // cell, requiring us to use an
+ // FEFaceValues for the
+ // neighbor and an
+ // FESubfaceValues for the
+ // current cell:
+ else if ((cell->neighbor(f)->level() == cell->level())
+ &&
+ (cell->neighbor(f)->has_children() == true))
+ {
+ for (unsigned int subface=0;
+ subface<cell->face(f)->n_children();
+ ++subface)
+ if (cell_is_in_fluid_domain (cell->neighbor_child_on_subface
+ (f, subface)))
+ {
+ elasticity_fe_subface_values.reinit (cell,
+ f,
+ subface);
+ stokes_fe_face_values.reinit (cell->neighbor_child_on_subface (f, subface),
+ cell->neighbor_of_neighbor(f));
+
+ assemble_interface_term (elasticity_fe_subface_values,
+ stokes_fe_face_values,
+ elasticity_phi,
+ stokes_phi_grads_u, stokes_phi_p,
+ local_interface_matrix);
+
+ cell->neighbor_child_on_subface (f, subface)
+ ->get_dof_indices (neighbor_dof_indices);
+ constraints.distribute_local_to_global(local_interface_matrix,
+ local_dof_indices,
+ neighbor_dof_indices,
+ system_matrix);
+ }
+ }
+ // The last option is that the
+ // neighbor is coarser. In that
+ // case we have to use an
+ // FESubfaceValues object for
+ // the neighbor and a
+ // FEFaceValues for the current
+ // cell; the rest is the same
+ // as before:
+ else if (cell->neighbor_is_coarser(f)
+ &&
+ cell_is_in_fluid_domain(cell->neighbor(f)))
+ {
+ elasticity_fe_face_values.reinit (cell, f);
+ stokes_fe_subface_values.reinit (cell->neighbor(f),
+ cell->neighbor_of_coarser_neighbor(f).first,
+ cell->neighbor_of_coarser_neighbor(f).second);
+
+ assemble_interface_term (elasticity_fe_face_values,
+ stokes_fe_subface_values,
+ elasticity_phi,
+ stokes_phi_grads_u, stokes_phi_p,
+ local_interface_matrix);
+
+ cell->neighbor(f)->get_dof_indices (neighbor_dof_indices);
+ constraints.distribute_local_to_global(local_interface_matrix,
+ local_dof_indices,
+ neighbor_dof_indices,
+ system_matrix);
+ }
+ }
+ }
+ }
- // In the function that assembles the global
- // system, we passed computing interface
- // terms to a separate function we discuss
- // here. The key is that even though we can't
- // predict the combination of FEFaceValues
- // and FESubfaceValues objects, they are both
- // derived from the FEFaceValuesBase class
- // and consequently we don't have to care:
- // the function is simply called with two
- // such objects denoting the values of the
- // shape functions on the quadrature points
- // of the two sides of the face. We then do
- // what we always do: we fill the scratch
- // arrays with the values of shape functions
- // and their derivatives, and then loop over
- // all entries of the matrix to compute the
- // local integrals. The details of the
- // bilinear form we evaluate here are given
- // in the introduction.
-template <int dim>
-void
-FluidStructureProblem<dim>::
-assemble_interface_term (const FEFaceValuesBase<dim> &elasticity_fe_face_values,
- const FEFaceValuesBase<dim> &stokes_fe_face_values,
- std::vector<Tensor<1,dim> > &elasticity_phi,
- std::vector<SymmetricTensor<2,dim> > &stokes_phi_grads_u,
- std::vector<double> &stokes_phi_p,
- FullMatrix<double> &local_interface_matrix) const
-{
- Assert (stokes_fe_face_values.n_quadrature_points ==
- elasticity_fe_face_values.n_quadrature_points,
- ExcInternalError());
- const unsigned int n_face_quadrature_points
- = elasticity_fe_face_values.n_quadrature_points;
-
- const FEValuesExtractors::Vector velocities (0);
- const FEValuesExtractors::Scalar pressure (dim);
- const FEValuesExtractors::Vector displacements (dim+1);
-
- local_interface_matrix = 0;
- for (unsigned int q=0; q<n_face_quadrature_points; ++q)
- {
- const Tensor<1,dim> normal_vector = stokes_fe_face_values.normal_vector(q);
-
- for (unsigned int k=0; k<stokes_fe_face_values.dofs_per_cell; ++k)
- stokes_phi_grads_u[k] = stokes_fe_face_values[velocities].symmetric_gradient (k, q);
- for (unsigned int k=0; k<elasticity_fe_face_values.dofs_per_cell; ++k)
- elasticity_phi[k] = elasticity_fe_face_values[displacements].value (k,q);
-
- for (unsigned int i=0; i<elasticity_fe_face_values.dofs_per_cell; ++i)
- for (unsigned int j=0; j<stokes_fe_face_values.dofs_per_cell; ++j)
- local_interface_matrix(i,j) += -((2 * viscosity *
- (stokes_phi_grads_u[j] *
- normal_vector)
- +
- stokes_phi_p[j] *
- normal_vector) *
- elasticity_phi[i] *
- stokes_fe_face_values.JxW(q));
- }
-}
- // @sect4{<code>FluidStructureProblem::solve</code>}
+ // In the function that assembles the global
+ // system, we passed computing interface
+ // terms to a separate function we discuss
+ // here. The key is that even though we can't
+ // predict the combination of FEFaceValues
+ // and FESubfaceValues objects, they are both
+ // derived from the FEFaceValuesBase class
+ // and consequently we don't have to care:
+ // the function is simply called with two
+ // such objects denoting the values of the
+ // shape functions on the quadrature points
+ // of the two sides of the face. We then do
+ // what we always do: we fill the scratch
+ // arrays with the values of shape functions
+ // and their derivatives, and then loop over
+ // all entries of the matrix to compute the
+ // local integrals. The details of the
+ // bilinear form we evaluate here are given
+ // in the introduction.
+ template <int dim>
+ void
+ FluidStructureProblem<dim>::
+ assemble_interface_term (const FEFaceValuesBase<dim> &elasticity_fe_face_values,
+ const FEFaceValuesBase<dim> &stokes_fe_face_values,
+ std::vector<Tensor<1,dim> > &elasticity_phi,
+ std::vector<SymmetricTensor<2,dim> > &stokes_phi_grads_u,
+ std::vector<double> &stokes_phi_p,
+ FullMatrix<double> &local_interface_matrix) const
+ {
+ Assert (stokes_fe_face_values.n_quadrature_points ==
+ elasticity_fe_face_values.n_quadrature_points,
+ ExcInternalError());
+ const unsigned int n_face_quadrature_points
+ = elasticity_fe_face_values.n_quadrature_points;
+
+ const FEValuesExtractors::Vector velocities (0);
+ const FEValuesExtractors::Scalar pressure (dim);
+ const FEValuesExtractors::Vector displacements (dim+1);
+
+ local_interface_matrix = 0;
+ for (unsigned int q=0; q<n_face_quadrature_points; ++q)
+ {
+ const Tensor<1,dim> normal_vector = stokes_fe_face_values.normal_vector(q);
+
+ for (unsigned int k=0; k<stokes_fe_face_values.dofs_per_cell; ++k)
+ stokes_phi_grads_u[k] = stokes_fe_face_values[velocities].symmetric_gradient (k, q);
+ for (unsigned int k=0; k<elasticity_fe_face_values.dofs_per_cell; ++k)
+ elasticity_phi[k] = elasticity_fe_face_values[displacements].value (k,q);
+
+ for (unsigned int i=0; i<elasticity_fe_face_values.dofs_per_cell; ++i)
+ for (unsigned int j=0; j<stokes_fe_face_values.dofs_per_cell; ++j)
+ local_interface_matrix(i,j) += -((2 * viscosity *
+ (stokes_phi_grads_u[j] *
+ normal_vector)
+ +
+ stokes_phi_p[j] *
+ normal_vector) *
+ elasticity_phi[i] *
+ stokes_fe_face_values.JxW(q));
+ }
+ }
- // As discussed in the introduction, we use a
- // rather trivial solver here: we just pass
- // the linear system off to the
- // SparseDirectUMFPACK direct solver (see,
- // for example, step-29). The only thing we
- // have to do after solving is ensure that
- // hanging node and boundary value
- // constraints are correct.
-template <int dim>
-void
-FluidStructureProblem<dim>::solve ()
-{
- SparseDirectUMFPACK direct_solver;
- direct_solver.initialize (system_matrix);
- direct_solver.vmult (solution, system_rhs);
- constraints.distribute (solution);
-}
+ // @sect4{<code>FluidStructureProblem::solve</code>}
+
+ // As discussed in the introduction, we use a
+ // rather trivial solver here: we just pass
+ // the linear system off to the
+ // SparseDirectUMFPACK direct solver (see,
+ // for example, step-29). The only thing we
+ // have to do after solving is ensure that
+ // hanging node and boundary value
+ // constraints are correct.
+ template <int dim>
+ void
+ FluidStructureProblem<dim>::solve ()
+ {
+ SparseDirectUMFPACK direct_solver;
+ direct_solver.initialize (system_matrix);
+ direct_solver.vmult (solution, system_rhs);
+ constraints.distribute (solution);
+ }
- // @sect4{<code>FluidStructureProblem::output_results</code>}
- // Generating graphical output is rather
- // trivial here: all we have to do is
- // identify which components of the solution
- // vector belong to scalars and/or vectors
- // (see, for example, step-22 for a previous
- // example), and then pass it all on to the
- // DataOut class (with the second template
- // argument equal to hp::DoFHandler instead
- // of the usual default DoFHandler):
-template <int dim>
-void
-FluidStructureProblem<dim>::
-output_results (const unsigned int refinement_cycle) const
-{
- std::vector<std::string> solution_names (dim, "velocity");
- solution_names.push_back ("pressure");
- for (unsigned int d=0; d<dim; ++d)
- solution_names.push_back ("displacement");
+ // @sect4{<code>FluidStructureProblem::output_results</code>}
- std::vector<DataComponentInterpretation::DataComponentInterpretation>
- data_component_interpretation
- (dim, DataComponentInterpretation::component_is_part_of_vector);
- data_component_interpretation
- .push_back (DataComponentInterpretation::component_is_scalar);
- for (unsigned int d=0; d<dim; ++d)
+ // Generating graphical output is rather
+ // trivial here: all we have to do is
+ // identify which components of the solution
+ // vector belong to scalars and/or vectors
+ // (see, for example, step-22 for a previous
+ // example), and then pass it all on to the
+ // DataOut class (with the second template
+ // argument equal to hp::DoFHandler instead
+ // of the usual default DoFHandler):
+ template <int dim>
+ void
+ FluidStructureProblem<dim>::
+ output_results (const unsigned int refinement_cycle) const
+ {
+ std::vector<std::string> solution_names (dim, "velocity");
+ solution_names.push_back ("pressure");
+ for (unsigned int d=0; d<dim; ++d)
+ solution_names.push_back ("displacement");
+
+ std::vector<DataComponentInterpretation::DataComponentInterpretation>
+ data_component_interpretation
+ (dim, DataComponentInterpretation::component_is_part_of_vector);
data_component_interpretation
- .push_back (DataComponentInterpretation::component_is_part_of_vector);
+ .push_back (DataComponentInterpretation::component_is_scalar);
+ for (unsigned int d=0; d<dim; ++d)
+ data_component_interpretation
+ .push_back (DataComponentInterpretation::component_is_part_of_vector);
- DataOut<dim,hp::DoFHandler<dim> > data_out;
- data_out.attach_dof_handler (dof_handler);
+ DataOut<dim,hp::DoFHandler<dim> > data_out;
+ data_out.attach_dof_handler (dof_handler);
- data_out.add_data_vector (solution, solution_names,
- DataOut<dim,hp::DoFHandler<dim> >::type_dof_data,
- data_component_interpretation);
- data_out.build_patches ();
+ data_out.add_data_vector (solution, solution_names,
+ DataOut<dim,hp::DoFHandler<dim> >::type_dof_data,
+ data_component_interpretation);
+ data_out.build_patches ();
- std::ostringstream filename;
- filename << "solution-"
- << Utilities::int_to_string (refinement_cycle, 2)
- << ".vtk";
+ std::ostringstream filename;
+ filename << "solution-"
+ << Utilities::int_to_string (refinement_cycle, 2)
+ << ".vtk";
- std::ofstream output (filename.str().c_str());
- data_out.write_vtk (output);
-}
+ std::ofstream output (filename.str().c_str());
+ data_out.write_vtk (output);
+ }
- // @sect4{<code>FluidStructureProblem::refine_mesh</code>}
-
- // The next step is to refine the mesh. As
- // was discussed in the introduction, this is
- // a bit tricky primarily because the fluid
- // and the solid subdomains use variables
- // that have different physical dimensions
- // and for which the absolute magnitude of
- // error estimates is consequently not
- // directly comparable. We will therefore
- // have to scale them. At the top of the
- // function, we therefore first compute error
- // estimates for the different variables
- // separately (using the velocities but not
- // the pressure for the fluid domain, and the
- // displacements in the solid domain):
-template <int dim>
-void
-FluidStructureProblem<dim>::refine_mesh ()
-{
- Vector<float>
- stokes_estimated_error_per_cell (triangulation.n_active_cells());
- Vector<float>
- elasticity_estimated_error_per_cell (triangulation.n_active_cells());
-
- const QGauss<dim-1> stokes_face_quadrature(stokes_degree+2);
- const QGauss<dim-1> elasticity_face_quadrature(elasticity_degree+2);
-
- hp::QCollection<dim-1> face_q_collection;
- face_q_collection.push_back (stokes_face_quadrature);
- face_q_collection.push_back (elasticity_face_quadrature);
-
- std::vector<bool> stokes_component_mask (dim+1+dim, false);
- for (unsigned int d=0; d<dim; ++d)
- stokes_component_mask[d] = true;
- KellyErrorEstimator<dim>::estimate (dof_handler,
- face_q_collection,
- typename FunctionMap<dim>::type(),
- solution,
- stokes_estimated_error_per_cell,
- stokes_component_mask);
-
- std::vector<bool> elasticity_component_mask (dim+1+dim, false);
- for (unsigned int d=0; d<dim; ++d)
- elasticity_component_mask[dim+1+d] = true;
- KellyErrorEstimator<dim>::estimate (dof_handler,
- face_q_collection,
- typename FunctionMap<dim>::type(),
- solution,
- elasticity_estimated_error_per_cell,
- elasticity_component_mask);
-
- // We then normalize error estimates by
- // dividing by their norm and scale the
- // fluid error indicators by a factor of 4
- // as discussed in the introduction. The
- // results are then added together into a
- // vector that contains error indicators
- // for all cells:
- stokes_estimated_error_per_cell
- *= 4. / stokes_estimated_error_per_cell.l2_norm();
- elasticity_estimated_error_per_cell
- *= 1. / elasticity_estimated_error_per_cell.l2_norm();
-
- Vector<float>
- estimated_error_per_cell (triangulation.n_active_cells());
-
- estimated_error_per_cell += stokes_estimated_error_per_cell;
- estimated_error_per_cell += elasticity_estimated_error_per_cell;
-
- // The second to last part of the function,
- // before actually refining the mesh,
- // involves a heuristic that we have
- // already mentioned in the introduction:
- // because the solution is discontinuous,
- // the KellyErrorEstimator class gets all
- // confused about cells that sit at the
- // boundary between subdomains: it believes
- // that the error is large there because
- // the jump in the gradient is large, even
- // though this is entirely expected and a
- // feature that is in fact present in the
- // exact solution as well and therefore not
- // indicative of any numerical error.
- //
- // Consequently, we set the error
- // indicators to zero for all cells at the
- // interface; the conditions determining
- // which cells this affects are slightly
- // awkward because we have to account for
- // the possibility of adaptively refined
- // meshes, meaning that the neighboring
- // cell can be coarser than the current
- // one, or could in fact be refined some
- // more. The structure of these nested
- // conditions is much the same as we
- // encountered when assembling interface
- // terms in <code>assemble_system</code>.
+ // @sect4{<code>FluidStructureProblem::refine_mesh</code>}
+
+ // The next step is to refine the mesh. As
+ // was discussed in the introduction, this is
+ // a bit tricky primarily because the fluid
+ // and the solid subdomains use variables
+ // that have different physical dimensions
+ // and for which the absolute magnitude of
+ // error estimates is consequently not
+ // directly comparable. We will therefore
+ // have to scale them. At the top of the
+ // function, we therefore first compute error
+ // estimates for the different variables
+ // separately (using the velocities but not
+ // the pressure for the fluid domain, and the
+ // displacements in the solid domain):
+ template <int dim>
+ void
+ FluidStructureProblem<dim>::refine_mesh ()
{
- unsigned int cell_index = 0;
- for (typename hp::DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active();
- cell != dof_handler.end(); ++cell, ++cell_index)
- for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
- if (cell_is_in_solid_domain (cell))
- {
- if ((cell->at_boundary(f) == false)
- &&
- (((cell->neighbor(f)->level() == cell->level())
- &&
- (cell->neighbor(f)->has_children() == false)
- &&
- cell_is_in_fluid_domain (cell->neighbor(f)))
- ||
- ((cell->neighbor(f)->level() == cell->level())
- &&
- (cell->neighbor(f)->has_children() == true)
- &&
- (cell_is_in_fluid_domain (cell->neighbor_child_on_subface
- (f, 0))))
- ||
- (cell->neighbor_is_coarser(f)
- &&
- cell_is_in_fluid_domain(cell->neighbor(f)))
- ))
- estimated_error_per_cell(cell_index) = 0;
- }
- else
- {
- if ((cell->at_boundary(f) == false)
- &&
- (((cell->neighbor(f)->level() == cell->level())
- &&
- (cell->neighbor(f)->has_children() == false)
- &&
- cell_is_in_solid_domain (cell->neighbor(f)))
- ||
- ((cell->neighbor(f)->level() == cell->level())
- &&
- (cell->neighbor(f)->has_children() == true)
+ Vector<float>
+ stokes_estimated_error_per_cell (triangulation.n_active_cells());
+ Vector<float>
+ elasticity_estimated_error_per_cell (triangulation.n_active_cells());
+
+ const QGauss<dim-1> stokes_face_quadrature(stokes_degree+2);
+ const QGauss<dim-1> elasticity_face_quadrature(elasticity_degree+2);
+
+ hp::QCollection<dim-1> face_q_collection;
+ face_q_collection.push_back (stokes_face_quadrature);
+ face_q_collection.push_back (elasticity_face_quadrature);
+
+ std::vector<bool> stokes_component_mask (dim+1+dim, false);
+ for (unsigned int d=0; d<dim; ++d)
+ stokes_component_mask[d] = true;
+ KellyErrorEstimator<dim>::estimate (dof_handler,
+ face_q_collection,
+ typename FunctionMap<dim>::type(),
+ solution,
+ stokes_estimated_error_per_cell,
+ stokes_component_mask);
+
+ std::vector<bool> elasticity_component_mask (dim+1+dim, false);
+ for (unsigned int d=0; d<dim; ++d)
+ elasticity_component_mask[dim+1+d] = true;
+ KellyErrorEstimator<dim>::estimate (dof_handler,
+ face_q_collection,
+ typename FunctionMap<dim>::type(),
+ solution,
+ elasticity_estimated_error_per_cell,
+ elasticity_component_mask);
+
+ // We then normalize error estimates by
+ // dividing by their norm and scale the
+ // fluid error indicators by a factor of 4
+ // as discussed in the introduction. The
+ // results are then added together into a
+ // vector that contains error indicators
+ // for all cells:
+ stokes_estimated_error_per_cell
+ *= 4. / stokes_estimated_error_per_cell.l2_norm();
+ elasticity_estimated_error_per_cell
+ *= 1. / elasticity_estimated_error_per_cell.l2_norm();
+
+ Vector<float>
+ estimated_error_per_cell (triangulation.n_active_cells());
+
+ estimated_error_per_cell += stokes_estimated_error_per_cell;
+ estimated_error_per_cell += elasticity_estimated_error_per_cell;
+
+ // The second to last part of the function,
+ // before actually refining the mesh,
+ // involves a heuristic that we have
+ // already mentioned in the introduction:
+ // because the solution is discontinuous,
+ // the KellyErrorEstimator class gets all
+ // confused about cells that sit at the
+ // boundary between subdomains: it believes
+ // that the error is large there because
+ // the jump in the gradient is large, even
+ // though this is entirely expected and a
+ // feature that is in fact present in the
+ // exact solution as well and therefore not
+ // indicative of any numerical error.
+ //
+ // Consequently, we set the error
+ // indicators to zero for all cells at the
+ // interface; the conditions determining
+ // which cells this affects are slightly
+ // awkward because we have to account for
+ // the possibility of adaptively refined
+ // meshes, meaning that the neighboring
+ // cell can be coarser than the current
+ // one, or could in fact be refined some
+ // more. The structure of these nested
+ // conditions is much the same as we
+ // encountered when assembling interface
+ // terms in <code>assemble_system</code>.
+ {
+ unsigned int cell_index = 0;
+ for (typename hp::DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active();
+ cell != dof_handler.end(); ++cell, ++cell_index)
+ for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ if (cell_is_in_solid_domain (cell))
+ {
+ if ((cell->at_boundary(f) == false)
&&
- (cell_is_in_solid_domain (cell->neighbor_child_on_subface
- (f, 0))))
- ||
- (cell->neighbor_is_coarser(f)
+ (((cell->neighbor(f)->level() == cell->level())
+ &&
+ (cell->neighbor(f)->has_children() == false)
+ &&
+ cell_is_in_fluid_domain (cell->neighbor(f)))
+ ||
+ ((cell->neighbor(f)->level() == cell->level())
+ &&
+ (cell->neighbor(f)->has_children() == true)
+ &&
+ (cell_is_in_fluid_domain (cell->neighbor_child_on_subface
+ (f, 0))))
+ ||
+ (cell->neighbor_is_coarser(f)
+ &&
+ cell_is_in_fluid_domain(cell->neighbor(f)))
+ ))
+ estimated_error_per_cell(cell_index) = 0;
+ }
+ else
+ {
+ if ((cell->at_boundary(f) == false)
&&
- cell_is_in_solid_domain(cell->neighbor(f)))
- ))
- estimated_error_per_cell(cell_index) = 0;
- }
- }
+ (((cell->neighbor(f)->level() == cell->level())
+ &&
+ (cell->neighbor(f)->has_children() == false)
+ &&
+ cell_is_in_solid_domain (cell->neighbor(f)))
+ ||
+ ((cell->neighbor(f)->level() == cell->level())
+ &&
+ (cell->neighbor(f)->has_children() == true)
+ &&
+ (cell_is_in_solid_domain (cell->neighbor_child_on_subface
+ (f, 0))))
+ ||
+ (cell->neighbor_is_coarser(f)
+ &&
+ cell_is_in_solid_domain(cell->neighbor(f)))
+ ))
+ estimated_error_per_cell(cell_index) = 0;
+ }
+ }
- GridRefinement::refine_and_coarsen_fixed_number (triangulation,
- estimated_error_per_cell,
- 0.3, 0.0);
- triangulation.execute_coarsening_and_refinement ();
-}
+ GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+ estimated_error_per_cell,
+ 0.3, 0.0);
+ triangulation.execute_coarsening_and_refinement ();
+ }
- // @sect4{<code>FluidStructureProblem::run</code>}
+ // @sect4{<code>FluidStructureProblem::run</code>}
- // This is, as usual, the function that
- // controls the overall flow of operation. If
- // you've read through tutorial programs
- // step-1 through step-6, for example, then
- // you are already quite familiar with the
- // following structure:
-template <int dim>
-void FluidStructureProblem<dim>::run ()
-{
- make_grid ();
+ // This is, as usual, the function that
+ // controls the overall flow of operation. If
+ // you've read through tutorial programs
+ // step-1 through step-6, for example, then
+ // you are already quite familiar with the
+ // following structure:
+ template <int dim>
+ void FluidStructureProblem<dim>::run ()
+ {
+ make_grid ();
- for (unsigned int refinement_cycle = 0; refinement_cycle<10-2*dim;
- ++refinement_cycle)
- {
- std::cout << "Refinement cycle " << refinement_cycle << std::endl;
+ for (unsigned int refinement_cycle = 0; refinement_cycle<10-2*dim;
+ ++refinement_cycle)
+ {
+ std::cout << "Refinement cycle " << refinement_cycle << std::endl;
- if (refinement_cycle > 0)
- refine_mesh ();
+ if (refinement_cycle > 0)
+ refine_mesh ();
- setup_dofs ();
+ setup_dofs ();
- std::cout << " Assembling..." << std::endl;
- assemble_system ();
+ std::cout << " Assembling..." << std::endl;
+ assemble_system ();
- std::cout << " Solving..." << std::endl;
- solve ();
+ std::cout << " Solving..." << std::endl;
+ solve ();
- std::cout << " Writing output..." << std::endl;
- output_results (refinement_cycle);
+ std::cout << " Writing output..." << std::endl;
+ output_results (refinement_cycle);
- std::cout << std::endl;
- }
+ std::cout << std::endl;
+ }
+ }
}
{
try
{
+ using namespace dealii;
+ using namespace Step46;
+
deallog.depth_console (0);
FluidStructureProblem<2> flow_problem(1, 1);
#include <deal.II/numerics/error_estimator.h>
-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
+namespace Step47
{
- public:
- LaplaceProblem ();
- ~LaplaceProblem ();
-
- void run ();
+ using namespace dealii;
- 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;
+ double sign (double d)
+ {
+ if (d > 0)
+ return 1;
+ else if (d < 0)
+ return -1;
+ else
+ return 0;
+ }
- hp::DoFHandler<dim> dof_handler;
- hp::FECollection<dim> fe_collection;
- ConstraintMatrix constraints;
+ template <int dim>
+ class LaplaceProblem
+ {
+ public:
+ LaplaceProblem ();
+ ~LaplaceProblem ();
- SparsityPattern sparsity_pattern;
- SparseMatrix<double> system_matrix;
+ void run ();
- Vector<double> solution;
- Vector<double> system_rhs;
-};
+ 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;
-template <int dim>
-class Coefficient : public Function<dim>
-{
- public:
- Coefficient () : Function<dim>() {}
+ ConstraintMatrix constraints;
- virtual double value (const Point<dim> &p,
- const unsigned int component = 0) const;
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
- virtual void value_list (const std::vector<Point<dim> > &points,
- std::vector<double> &values,
- const unsigned int component = 0) const;
-};
+ Vector<double> solution;
+ Vector<double> system_rhs;
+ };
-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>
+ class Coefficient : public Function<dim>
+ {
+ public:
+ Coefficient () : Function<dim>() {}
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
-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();
+ virtual void value_list (const std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int component = 0) const;
+ };
- 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)
- {
- if (points[i].square() < 0.5*0.5)
- values[i] = 20;
- else
- values[i] = 1;
- }
-}
+ 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>
-double exact_solution (const Point<dim> &p)
-{
- const double r = p.norm();
+ 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();
- return (r < 0.5
- ?
- 1./20 * (-1./4*r*r + 61./16)
- :
- 1./4 * (1-r*r));
-}
+ Assert (values.size() == n_points,
+ ExcDimensionMismatch (values.size(), n_points));
+ Assert (component == 0,
+ ExcIndexRange (component, 0, 1));
-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));
-}
+ 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>
-LaplaceProblem<dim>::~LaplaceProblem ()
-{
- dof_handler.clear ();
-}
+ 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>
-double
-level_set (const Point<dim> &p)
-{
- return p.norm() - 0.5;
-}
+ 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>
-Tensor<1,dim>
-grad_level_set (const Point<dim> &p)
-{
- return p / p.norm();
-}
+ template <int dim>
+ LaplaceProblem<dim>::~LaplaceProblem ()
+ {
+ dof_handler.clear ();
+ }
-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>
+ double
+ level_set (const Point<dim> &p)
+ {
+ return p.norm() - 0.5;
+ }
-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);
+ template <int dim>
+ Tensor<1,dim>
+ grad_level_set (const Point<dim> &p)
+ {
+ return p / p.norm();
+ }
- dof_handler.distribute_dofs (fe_collection);
- solution.reinit (dof_handler.n_dofs());
- system_rhs.reinit (dof_handler.n_dofs());
+ 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;
- constraints.clear ();
-//TODO: fix this, it currently crashes
- // DoFTools::make_hanging_node_constraints (dof_handler,
- // constraints);
+ // we get here only if all vertices
+ // have the same sign, which means
+ // that the cell is not intersected
+ return false;
+ }
-//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);
+ 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);
- constraints.condense (c_sparsity);
+ dof_handler.distribute_dofs (fe_collection);
- sparsity_pattern.copy_from(c_sparsity);
+ solution.reinit (dof_handler.n_dofs());
+ system_rhs.reinit (dof_handler.n_dofs());
- system_matrix.reinit (sparsity_pattern);
-}
+ constraints.clear ();
+//TODO: fix this, it currently crashes
+ // DoFTools::make_hanging_node_constraints (dof_handler,
+ // constraints);
-template <int dim>
-void LaplaceProblem<dim>::assemble_system ()
-{
- const QGauss<dim> quadrature_formula(3);
+//TODO: component 1 must satisfy zero boundary conditions
+ constraints.close();
- FEValues<dim> plain_fe_values (fe_collection[0], quadrature_formula,
- update_values | update_gradients |
- update_quadrature_points | update_JxW_values);
- FEValues<dim> enriched_fe_values (fe_collection[1], quadrature_formula,
- update_values | update_gradients |
- update_quadrature_points | update_JxW_values);
+ CompressedSparsityPattern c_sparsity(dof_handler.n_dofs());
+ DoFTools::make_sparsity_pattern (dof_handler, c_sparsity);
- const unsigned int n_q_points = quadrature_formula.size();
+ constraints.condense (c_sparsity);
- FullMatrix<double> cell_matrix;
- Vector<double> cell_rhs;
+ sparsity_pattern.copy_from(c_sparsity);
- std::vector<unsigned int> local_dof_indices;
+ system_matrix.reinit (sparsity_pattern);
+ }
- 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();
+ template <int dim>
+ void LaplaceProblem<dim>::assemble_system ()
+ {
+ const QGauss<dim> quadrature_formula(3);
- 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;
+ FEValues<dim> plain_fe_values (fe_collection[0], quadrature_formula,
+ update_values | update_gradients |
+ update_quadrature_points | update_JxW_values);
+ FEValues<dim> enriched_fe_values (fe_collection[1], quadrature_formula,
+ update_values | update_gradients |
+ update_quadrature_points | update_JxW_values);
- if (cell->active_fe_index() == 0)
- {
- plain_fe_values.reinit (cell);
+ const unsigned int n_q_points = quadrature_formula.size();
- coefficient_values.resize (plain_fe_values.n_quadrature_points);
- coefficient.value_list (plain_fe_values.get_quadrature_points(),
- coefficient_values);
+ FullMatrix<double> cell_matrix;
+ Vector<double> cell_rhs;
- 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));
+ std::vector<unsigned int> local_dof_indices;
+ const Coefficient<dim> coefficient;
+ std::vector<double> coefficient_values (n_q_points);
- 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
-//TODO: remove update_support_points and friends, since they aren't implemented anyway
- Assert (cell->active_fe_index() == 1, ExcInternalError());
- Assert (interface_intersects_cell(cell) == true, ExcInternalError());
+ typename hp::DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
- 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));
+ 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);
- 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 );
+ cell_matrix = 0;
+ cell_rhs = 0;
- this_fe_values.reinit (cell);
+ if (cell->active_fe_index() == 0)
+ {
+ plain_fe_values.reinit (cell);
- coefficient_values.resize (this_fe_values.n_quadrature_points);
- coefficient.value_list (this_fe_values.get_quadrature_points(),
- coefficient_values);
+ 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<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 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)
- 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) *
+ 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 *
- this_fe_values.JxW(q_point));
+ 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
+//TODO: remove update_support_points and friends, since they aren't implemented anyway
+ 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);
- }
+ 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<unsigned int,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);
+ std::map<unsigned int,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
// 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;
- }
+ 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 )
+ {
- // 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 folloing three classes of decomposition
- // type 1: ++++, ----
- // type 2: -+++, +-++, ++-+, +++-, +---, -+--, --+-, ---+
- // type 3: +--+, ++--, +-+-, -++-, --++, -+-+
+ unsigned int type = 0;
- 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;
+ // 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;
+ }
- unsigned int Pos = 100;
+ // 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 folloing three classes of decomposition
+ // type 1: ++++, ----
+ // type 2: -+++, +-++, ++-+, +++-, +---, -+--, --+-, ---+
+ // type 3: +--+, ++--, +-+-, -++-, --++, -+-+
- Point<dim> v0(0,0);
- Point<dim> v1(1,0);
- Point<dim> v2(0,1);
- Point<dim> v3(1,1);
+ 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;
- 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);
- }
- }
+ unsigned int Pos = 100;
- Assert (xfem_quadrature.size() == plain_quadrature.size() * 5, ExcInternalError());
- return std::pair<unsigned int, Quadrature<dim> >(2, xfem_quadrature);
- }
+ Point<dim> v0(0,0);
+ Point<dim> v1(1,0);
+ Point<dim> v2(0,1);
+ Point<dim> v3(1,1);
- // Type three decomposition
- // (+--+, ++--, +-+-, -++-, --++, -+-+)
+ 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==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);
- }
+ if (type == 1)
+ return std::pair<unsigned int, Quadrature<dim> >(1, plain_quadrature);
- return std::pair<unsigned int, Quadrature<dim> >(0, 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
-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)
+ Quadrature<dim> xfem_quadrature(5*n_q_points);
-{
- // 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.
+ std::vector<Point<dim> > v(GeometryInfo<dim>::vertices_per_cell);
- unsigned int n_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
- 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);
+ // Find cut coordinates
- const unsigned int n_q_points = plain_quadrature.size();
+ // deal.ii local coordinates
- std::vector<double> JxW(n_q_points);
+ // 2-------3
+ // | |
+ // | |
+ // | |
+ // 0-------1
- for ( unsigned int i = 0; i < n_q_points; i++)
- {
- switch (dim)
- {
- case 2:
+ if (Pos == 0)
{
- 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;
+ 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) );
}
- default:
- Assert (false, ExcNotImplemented());
- }
-
-
- Tensor<2,dim> jacobian;
+ //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;
- // Calculate Jacobian of transformation
- for (unsigned int d=0; d<dim; ++d)
- for (unsigned int e=0; e<dim; ++e)
+ 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)
{
- for (unsigned int j = 0; j<GeometryInfo<dim>::vertices_per_cell; j++)
+ 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++)
{
- jacobian[d][e] += grad_phi[j][e] * v[j](d);
+ //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);
}
}
- double detJ = determinant(jacobian);
- xfem_weights.push_back (W[i] * detJ);
+ Assert (xfem_quadrature.size() == plain_quadrature.size() * 5, ExcInternalError());
+ return std::pair<unsigned int, Quadrature<dim> >(2, xfem_quadrature);
+ }
- // 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);
- }
+ // 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
-template <int dim>
-void LaplaceProblem<dim>::solve ()
-{
- SolverControl solver_control (1000, 1e-12);
- SolverCG<> solver (solver_control);
+ Quadrature<dim> xfem_quadrature(5*n_q_points);
+
+ std::vector<Point<dim> > v(GeometryInfo<dim>::vertices_per_cell);
- PreconditionSSOR<> preconditioner;
- preconditioner.initialize(system_matrix, 1.2);
+ 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);
+ }
- solver.solve (system_matrix, solution, system_rhs,
- preconditioner);
+ //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}}
+ };
- constraints.distribute (solution);
-}
+ //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>::refine_grid ()
-{
- Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+ 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)
- KellyErrorEstimator<dim>::estimate (dof_handler,
- QGauss<dim-1>(3),
- typename FunctionMap<dim>::type(),
- solution,
- estimated_error_per_cell);
+ {
+ // 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.
- GridRefinement::refine_and_coarsen_fixed_number (triangulation,
- estimated_error_per_cell,
- 0.3, 0.03);
+ unsigned int n_v = GeometryInfo<dim>::vertices_per_cell;
- triangulation.execute_coarsening_and_refinement ();
-}
+ 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);
-template <int dim>
-class Postprocessor : public DataPostprocessor<dim>
-{
- public:
- virtual
- void
- compute_derived_quantities_vector (const std::vector<Vector<double> > &uh,
- const std::vector<std::vector<Tensor<1,dim> > > &duh,
- const std::vector<std::vector<Tensor<2,dim> > > &dduh,
- const std::vector<Point<dim> > &normals,
- const std::vector<Point<dim> > &evaluation_points,
- std::vector<Vector<double> > &computed_quantities) const;
+ 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;
+ }
- virtual std::vector<std::string> get_names () const;
+ default:
+ Assert (false, ExcNotImplemented());
+ }
- virtual unsigned int n_output_variables() const;
- virtual
- std::vector<DataComponentInterpretation::DataComponentInterpretation>
- get_data_component_interpretation () const;
+ Tensor<2,dim> jacobian;
- virtual UpdateFlags get_needed_update_flags () const;
-};
+ // 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);
-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;
-}
+ // 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>
-unsigned int
-Postprocessor<dim>::n_output_variables() const
-{
- return get_names().size();
-}
+ template <int dim>
+ void LaplaceProblem<dim>::solve ()
+ {
+ SolverControl solver_control (1000, 1e-12);
+ SolverCG<> solver (solver_control);
-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;
-}
+ PreconditionSSOR<> preconditioner;
+ preconditioner.initialize(system_matrix, 1.2);
+ solver.solve (system_matrix, solution, system_rhs,
+ preconditioner);
-template <int dim>
-UpdateFlags
-Postprocessor<dim>::get_needed_update_flags() const
-{
- return update_values | update_q_points;
-}
+ constraints.distribute (solution);
+ }
-template <int dim>
-void
-Postprocessor<dim>::
-compute_derived_quantities_vector (const std::vector<Vector<double> > &uh,
- const std::vector<std::vector<Tensor<1,dim> > > &/*duh*/,
- const std::vector<std::vector<Tensor<2,dim> > > &/*dduh*/,
- const std::vector<Point<dim> > &/*normals*/,
- const std::vector<Point<dim> > &evaluation_points,
- std::vector<Vector<double> > &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());
- Assert (computed_quantities[0].size()==n_output_variables(),ExcInternalError());
- for (unsigned int q=0; q<n_quadrature_points; ++q)
- {
- computed_quantities[q](0)
- = (uh[q](0)
- +
- uh[q](1) * std::fabs(level_set(evaluation_points[q])));
- computed_quantities[q](1)
- = (computed_quantities[q](0)
- -
- exact_solution (evaluation_points[q]));
- }
-}
+ 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);
-template <int dim>
-void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
-{
- Assert (cycle < 10, ExcNotImplemented());
+ triangulation.execute_coarsening_and_refinement ();
+ }
- std::string filename = "solution-";
- filename += ('0' + cycle);
- //filename += ".vtk";
- filename += ".gmv";
- std::ofstream output (filename.c_str());
- Postprocessor<dim> postprocessor;
- DataOut<dim,hp::DoFHandler<dim> > data_out;
+ template <int dim>
+ class Postprocessor : public DataPostprocessor<dim>
+ {
+ public:
+ virtual
+ void
+ compute_derived_quantities_vector (const std::vector<Vector<double> > &uh,
+ const std::vector<std::vector<Tensor<1,dim> > > &duh,
+ const std::vector<std::vector<Tensor<2,dim> > > &dduh,
+ const std::vector<Point<dim> > &normals,
+ const std::vector<Point<dim> > &evaluation_points,
+ std::vector<Vector<double> > &computed_quantities) const;
- 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);
+ virtual std::vector<std::string> get_names () const;
- //data_out.write_vtk (output);
- data_out.write_gmv (output);
-}
+ virtual unsigned int n_output_variables() const;
+ virtual
+ std::vector<DataComponentInterpretation::DataComponentInterpretation>
+ get_data_component_interpretation () const;
+ virtual UpdateFlags get_needed_update_flags () const;
+ };
-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);
+ 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;
+ }
- double l2_error_square = 0;
- std::vector<Vector<double> > solution_values;
+ template <int dim>
+ unsigned int
+ Postprocessor<dim>::n_output_variables() const
+ {
+ return get_names().size();
+ }
- 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);
- }
- }
+ 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>::
+ compute_derived_quantities_vector (const std::vector<Vector<double> > &uh,
+ const std::vector<std::vector<Tensor<1,dim> > > &/*duh*/,
+ const std::vector<std::vector<Tensor<2,dim> > > &/*dduh*/,
+ const std::vector<Point<dim> > &/*normals*/,
+ const std::vector<Point<dim> > &evaluation_points,
+ std::vector<Vector<double> > &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());
+ Assert (computed_quantities[0].size()==n_output_variables(),ExcInternalError());
+
+ for (unsigned int q=0; q<n_quadrature_points; ++q)
+ {
+ computed_quantities[q](0)
+ = (uh[q](0)
+ +
+ uh[q](1) * std::fabs(level_set(evaluation_points[q])));
+ computed_quantities[q](1)
+ = (computed_quantities[q](0)
+ -
+ exact_solution (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";
+ filename += ".gmv";
+
+ 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);
+ data_out.write_gmv (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;
-}
+ 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;
+ 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);
+ if (cycle == 0)
+ {
+ GridGenerator::hyper_ball (triangulation);
+ //GridGenerator::hyper_cube (triangulation, -1, 1);
- static const HyperBallBoundary<dim> boundary;
- triangulation.set_boundary (0, boundary);
+ static const HyperBallBoundary<dim> boundary;
+ triangulation.set_boundary (0, boundary);
- triangulation.refine_global (2);
- }
- else
- triangulation.refine_global (1);
+ triangulation.refine_global (2);
+ }
+ else
+ triangulation.refine_global (1);
// refine_grid ();
- std::cout << " Number of active cells: "
- << triangulation.n_active_cells()
- << std::endl;
+ std::cout << " Number of active cells: "
+ << triangulation.n_active_cells()
+ << std::endl;
- setup_system ();
+ setup_system ();
- std::cout << " Number of degrees of freedom: "
- << dof_handler.n_dofs()
- << std::endl;
+ std::cout << " Number of degrees of freedom: "
+ << dof_handler.n_dofs()
+ << std::endl;
- assemble_system ();
- solve ();
- compute_error ();
- output_results (cycle);
- }
+ assemble_system ();
+ solve ();
+ compute_error ();
+ output_results (cycle);
+ }
+ }
}
try
{
+ using namespace dealii;
+ using namespace Step47;
+
deallog.depth_console (0);
LaplaceProblem<2> laplace_problem_2d;
projects / dealii-svn.git / commitdiff