From 2f18bc3295cf1c54fd0fff68c15447a123fba6b6 Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 9 Sep 2011 04:02:31 +0000 Subject: [PATCH] Move things local to each program into a local namespace. git-svn-id: https://svn.dealii.org/trunk@24295 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-40/step-40.cc | 1438 ++++++++--------- deal.II/examples/step-45/step-45.cc | 836 +++++----- deal.II/examples/step-46/step-46.cc | 2203 ++++++++++++++------------- deal.II/examples/step-47/step-47.cc | 1736 ++++++++++----------- 4 files changed, 3119 insertions(+), 3094 deletions(-) diff --git a/deal.II/examples/step-40/step-40.cc b/deal.II/examples/step-40/step-40.cc index fcc14db60c..732097b0d9 100644 --- a/deal.II/examples/step-40/step-40.cc +++ b/deal.II/examples/step-40/step-40.cc @@ -4,7 +4,7 @@ /* $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 */ @@ -125,743 +125,746 @@ #include #include -using namespace dealii; - - // @sect3{The LaplaceProblem 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 mpi_communicator - // 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 pcout - // 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 -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 triangulation; - - DoFHandler dof_handler; - FE_Q fe; - - IndexSet locally_owned_dofs; - IndexSet locally_relevant_dofs; - - ConstraintMatrix constraints; + using namespace dealii; + + // @sect3{The LaplaceProblem 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 mpi_communicator + // 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 pcout + // 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 + 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 triangulation; + + DoFHandler dof_handler; + FE_Q 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 LaplaceProblem 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 pcout variable + // to only allow processor zero to output + // anything: + template + LaplaceProblem::LaplaceProblem () + : + mpi_communicator (MPI_COMM_WORLD), + triangulation (mpi_communicator, + typename Triangulation::MeshSmoothing + (Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening)), + dof_handler (triangulation), + fe (2), + pcout (std::cout, + (Utilities::System:: + get_this_mpi_process(mpi_communicator) + == 0)) + {} + + + + template + LaplaceProblem::~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 + void LaplaceProblem::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(), + 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 cell->is_ghost() + // || cell->is_artificial() 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 + void LaplaceProblem::assemble_system () + { + const QGauss quadrature_formula(3); + + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | + update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::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 + 0.5+0.25*std::sin(4.0 * numbers::PI * + fe_values.quadrature_point(q_point)[0]) + ? 1 : -1); + + for (unsigned int i=0; iget_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 + void LaplaceProblem::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 LaplaceProblem 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 pcout variable - // to only allow processor zero to output - // anything: -template -LaplaceProblem::LaplaceProblem () - : - mpi_communicator (MPI_COMM_WORLD), - triangulation (mpi_communicator, - typename Triangulation::MeshSmoothing - (Triangulation::smoothing_on_refinement | - Triangulation::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 -LaplaceProblem::~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 -void LaplaceProblem::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(), - 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 + void LaplaceProblem::refine_grid () + { + Vector estimated_error_per_cell (triangulation.n_active_cells()); + KellyErrorEstimator::estimate (dof_handler, + QGauss(3), + typename FunctionMap::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 cell->is_ghost() - // || cell->is_artificial() 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 -void LaplaceProblem::assemble_system () -{ - const QGauss quadrature_formula(3); - - FEValues fe_values (fe, quadrature_formula, - update_values | update_gradients | - update_quadrature_points | - update_JxW_values); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); - - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); - Vector cell_rhs (dofs_per_cell); - - std::vector local_dof_indices (dofs_per_cell); - - typename DoFHandler::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 + void LaplaceProblem::output_results (const unsigned int cycle) const + { + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (locally_relevant_solution, "u"); + + Vector subdomain (triangulation.n_active_cells()); + for (unsigned int i=0; isolution-XX-PPPP.vtu + // where XX indicates + // the refinement cycle, + // PPPP 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 + // .vtu 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 filenames; + for (unsigned int i=0; + ipcout instead of + // std::cout 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 + void LaplaceProblem::run () + { + const unsigned int n_cycles = 8; + for (unsigned int cycle=0; cycle - 0.5+0.25*std::sin(4.0 * numbers::PI * - fe_values.quadrature_point(q_point)[0]) - ? 1 : -1); - - for (unsigned int i=0; iget_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 -void LaplaceProblem::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 -void LaplaceProblem::refine_grid () -{ - Vector estimated_error_per_cell (triangulation.n_active_cells()); - KellyErrorEstimator::estimate (dof_handler, - QGauss(3), - typename FunctionMap::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 -void LaplaceProblem::output_results (const unsigned int cycle) const -{ - DataOut data_out; - data_out.attach_dof_handler (dof_handler); - data_out.add_data_vector (locally_relevant_solution, "u"); - - Vector subdomain (triangulation.n_active_cells()); - for (unsigned int i=0; isolution-XX-PPPP.vtu - // where XX indicates - // the refinement cycle, - // PPPP 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 - // .vtu 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 filenames; - for (unsigned int i=0; - ipcout instead of - // std::cout 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 -void LaplaceProblem::run () -{ - const unsigned int n_cycles = 8; - for (unsigned int cycle=0; cycle laplace_problem_2d; laplace_problem_2d.run (); } - + PetscFinalize(); } catch (std::exception &exc) diff --git a/deal.II/examples/step-45/step-45.cc b/deal.II/examples/step-45/step-45.cc index 6a7a6b086d..56494a2423 100644 --- a/deal.II/examples/step-45/step-45.cc +++ b/deal.II/examples/step-45/step-45.cc @@ -3,7 +3,7 @@ /* $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 */ @@ -44,443 +44,446 @@ #include -using namespace dealii; - - // @sect3{The LaplaceProblem class} - - // The class LaplaceProblem 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 - // constraints variables that - // will hold the constraints from the - // periodic boundary condition. We will - // initialize it in the - // make_periodicity_constraints() - // function which we call from - // make_grid_and_dofs(). -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 system_matrix; - Vector system_rhs; - Vector solution; - - void assemble_system (); - void output_results (); - void make_grid_and_dofs (); - void make_periodicity_constraints (); - void solve (); -}; - - - // @sect3{The RightHandSide 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 LaplaceProblem class} + + // The class LaplaceProblem 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 + // constraints variables that + // will hold the constraints from the + // periodic boundary condition. We will + // initialize it in the + // make_periodicity_constraints() + // function which we call from + // make_grid_and_dofs(). + class LaplaceProblem + { + public: + LaplaceProblem (); + void run (); + + private: + Triangulation<2> triangulation; + + FE_Q<2> fe; + DoFHandler<2> dof_handler; + + ConstraintMatrix constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + Vector system_rhs; + Vector solution; + + void assemble_system (); + void output_results (); + void make_grid_and_dofs (); + void make_periodicity_constraints (); + void solve (); + }; + + + // @sect3{The RightHandSide 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 LaplaceProblem 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 run(). 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 + // ConstraintMatrix 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 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 + // p of this map, + // p-@>first corresponds + // to the "key" of an entry (the + // global number of the degree of + // freedom), whereas + // p-@>second 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::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 cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); - // @sect3{Implementation of the LaplaceProblem 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 run(). 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 - // ConstraintMatrix 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 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 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 - // p of this map, - // p-@>first corresponds - // to the "key" of an entry (the - // global number of the degree of - // freedom), whereas - // p-@>second 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::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_pointget_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 cell_matrix (dofs_per_cell, dofs_per_cell); - Vector 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 > precondition; - std::vector 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_pointget_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 > 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 main function} @@ -491,6 +494,9 @@ int main () { try { + using namespace dealii; + using namespace Step45; + deallog.depth_console (0); LaplaceProblem laplace_problem; @@ -509,7 +515,7 @@ int main () return 1; } - catch (...) + catch (...) { std::cerr << std::endl << std::endl << "----------------------------------------------------" diff --git a/deal.II/examples/step-46/step-46.cc b/deal.II/examples/step-46/step-46.cc index 52a0882115..75dc091cfb 100644 --- a/deal.II/examples/step-46/step-46.cc +++ b/deal.II/examples/step-46/step-46.cc @@ -56,1196 +56,1200 @@ #include #include -using namespace dealii; - - // @sect3{The FluidStructureProblem 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 - // (cell_is_in_fluid_domain, - // cell_is_in_solid_domain) 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 - // (make_grid, set_active_fe_indices, - // assemble_interface_terms) 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 - // (viscosity, lambda, eta) - // describes the material properties used for - // the two physics models. -template -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::cell_iterator &cell); - - static bool - cell_is_in_solid_domain (const typename hp::DoFHandler::cell_iterator &cell); - - - void make_grid (); - void set_active_fe_indices (); - void setup_dofs (); - void assemble_system (); - void assemble_interface_term (const FEFaceValuesBase &elasticity_fe_face_values, - const FEFaceValuesBase &stokes_fe_face_values, - std::vector > &elasticity_phi, - std::vector > &stokes_phi_grads_u, - std::vector &stokes_phi_p, - FullMatrix &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 triangulation; - FESystem stokes_fe; - FESystem elasticity_fe; - hp::FECollection fe_collection; - hp::DoFHandler dof_handler; - - ConstraintMatrix constraints; - - SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; - - Vector solution; - Vector 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 -class StokesBoundaryValues : public Function -{ - public: - StokesBoundaryValues () : Function(dim+1+dim) {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; - - virtual void vector_value (const Point &p, - Vector &value) const; -}; - -template -double -StokesBoundaryValues::value (const Point &p, - const unsigned int component) const +namespace Step46 { - Assert (component < this->n_components, - ExcIndexRange (component, 0, this->n_components)); + using namespace dealii; + + // @sect3{The FluidStructureProblem 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 + // (cell_is_in_fluid_domain, + // cell_is_in_solid_domain) 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 + // (make_grid, set_active_fe_indices, + // assemble_interface_terms) 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 + // (viscosity, lambda, eta) + // describes the material properties used for + // the two physics models. + template + 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::cell_iterator &cell); + + static bool + cell_is_in_solid_domain (const typename hp::DoFHandler::cell_iterator &cell); + + + void make_grid (); + void set_active_fe_indices (); + void setup_dofs (); + void assemble_system (); + void assemble_interface_term (const FEFaceValuesBase &elasticity_fe_face_values, + const FEFaceValuesBase &stokes_fe_face_values, + std::vector > &elasticity_phi, + std::vector > &stokes_phi_grads_u, + std::vector &stokes_phi_p, + FullMatrix &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 triangulation; + FESystem stokes_fe; + FESystem elasticity_fe; + hp::FECollection fe_collection; + hp::DoFHandler dof_handler; + + ConstraintMatrix constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector 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 + class StokesBoundaryValues : public Function + { + public: + StokesBoundaryValues () : Function(dim+1+dim) {} - return 0; -} + virtual double value (const Point &p, + const unsigned int component = 0) const; + virtual void vector_value (const Point &p, + Vector &value) const; + }; -template -void -StokesBoundaryValues::vector_value (const Point &p, - Vector &values) const -{ - for (unsigned int c=0; cn_components; ++c) - values(c) = StokesBoundaryValues::value (p, c); -} + template + double + StokesBoundaryValues::value (const Point &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 -class RightHandSide : public Function -{ - public: - RightHandSide () : Function(dim+1) {} + return 0; + } - virtual double value (const Point &p, - const unsigned int component = 0) const; - virtual void vector_value (const Point &p, - Vector &value) const; + template + void + StokesBoundaryValues::vector_value (const Point &p, + Vector &values) const + { + for (unsigned int c=0; cn_components; ++c) + values(c) = StokesBoundaryValues::value (p, c); + } -}; -template -double -RightHandSide::value (const Point &/*p*/, - const unsigned int /*component*/) const -{ - return 0; -} + template + class RightHandSide : public Function + { + public: + RightHandSide () : Function(dim+1) {} + virtual double value (const Point &p, + const unsigned int component = 0) const; -template -void -RightHandSide::vector_value (const Point &p, - Vector &values) const -{ - for (unsigned int c=0; cn_components; ++c) - values(c) = RightHandSide::value (p, c); -} + virtual void vector_value (const Point &p, + Vector &value) const; + }; - // @sect3{The FluidStructureProblem 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 -FluidStructureProblem:: -FluidStructureProblem (const unsigned int stokes_degree, - const unsigned int elasticity_degree) - : - stokes_degree (stokes_degree), - elasticity_degree (elasticity_degree), - triangulation (Triangulation::maximum_smoothing), - stokes_fe (FE_Q(stokes_degree+1), dim, - FE_Q(stokes_degree), 1, - FE_Nothing(), dim), - elasticity_fe (FE_Nothing(), dim, - FE_Nothing(), 1, - FE_Q(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 + double + RightHandSide::value (const Point &/*p*/, + const unsigned int /*component*/) const + { + return 0; + } + template + void + RightHandSide::vector_value (const Point &p, + Vector &values) const + { + for (unsigned int c=0; cn_components; ++c) + values(c) = RightHandSide::value (p, c); + } -template -bool -FluidStructureProblem:: -cell_is_in_fluid_domain (const typename hp::DoFHandler::cell_iterator &cell) -{ - return (cell->material_id() == fluid_domain_id); -} + // @sect3{The FluidStructureProblem 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 + FluidStructureProblem:: + FluidStructureProblem (const unsigned int stokes_degree, + const unsigned int elasticity_degree) + : + stokes_degree (stokes_degree), + elasticity_degree (elasticity_degree), + triangulation (Triangulation::maximum_smoothing), + stokes_fe (FE_Q(stokes_degree+1), dim, + FE_Q(stokes_degree), 1, + FE_Nothing(), dim), + elasticity_fe (FE_Nothing(), dim, + FE_Nothing(), 1, + FE_Q(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 -bool -FluidStructureProblem:: -cell_is_in_solid_domain (const typename hp::DoFHandler::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. make_grid, 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 -void -FluidStructureProblem::make_grid () -{ - GridGenerator::subdivided_hyper_cube (triangulation, 8, -1, 1); - - for (typename Triangulation::active_cell_iterator - cell = triangulation.begin_active(); - cell != triangulation.end(); ++cell) - for (unsigned int f=0; f::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::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 + bool + FluidStructureProblem:: + cell_is_in_fluid_domain (const typename hp::DoFHandler::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 -void -FluidStructureProblem::set_active_fe_indices () -{ - for (typename hp::DoFHandler::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 + bool + FluidStructureProblem:: + cell_is_in_solid_domain (const typename hp::DoFHandler::cell_iterator &cell) + { + return (cell->material_id() == solid_domain_id); + } - // @sect4{FluidStructureProblem::setup_dofs} - - // 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 -void -FluidStructureProblem::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. make_grid, 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 + void + FluidStructureProblem::make_grid () { - constraints.clear (); - DoFTools::make_hanging_node_constraints (dof_handler, - constraints); + GridGenerator::subdivided_hyper_cube (triangulation, 8, -1, 1); - std::vector velocity_mask (dim+1+dim, false); - for (unsigned int d=0; d(), - constraints, - velocity_mask); - - std::vector elasticity_mask (dim+1+dim, false); - for (unsigned int d=dim+1; d(dim+1+dim), - constraints, - elasticity_mask); + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + for (unsigned int f=0; f::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::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 + void + FluidStructureProblem::set_active_fe_indices () { - std::vector local_face_dof_indices (stokes_fe.dofs_per_face); for (typename hp::DoFHandler::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::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; sfface(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; iset_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{FluidStructureProblem::setup_dofs} + + // 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 + void + FluidStructureProblem::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 velocity_mask (dim+1+dim, false); + for (unsigned int d=0; d(), + constraints, + velocity_mask); + + std::vector elasticity_mask (dim+1+dim, false); + for (unsigned int d=dim+1; d(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 local_face_dof_indices (stokes_fe.dofs_per_face); + for (typename hp::DoFHandler::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::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; sfface(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=dim+1) && (d>=dim+1))) - cell_coupling[c][d] = DoFTools::always; - - if ((c>=dim+1) && (dface(f)->get_dof_indices (local_face_dof_indices, 0); + for (unsigned int i=0; i 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=dim+1) && (d>=dim+1))) + cell_coupling[c][d] = DoFTools::always; + + if ((c>=dim+1) && (dFluidStructureProblem::assemble_system} + // 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 + void FluidStructureProblem::assemble_system () + { + system_matrix=0; + system_rhs=0; + + const QGauss stokes_quadrature(stokes_degree+2); + const QGauss elasticity_quadrature(elasticity_degree+2); + + hp::QCollection q_collection; + q_collection.push_back (stokes_quadrature); + q_collection.push_back (elasticity_quadrature); + + hp::FEValues hp_fe_values (fe_collection, q_collection, + update_values | + update_quadrature_points | + update_JxW_values | + update_gradients); + + const QGauss common_face_quadrature(std::max (stokes_degree+2, + elasticity_degree+2)); + + FEFaceValues stokes_fe_face_values (stokes_fe, + common_face_quadrature, + update_JxW_values | + update_normal_vectors | + update_gradients); + FEFaceValues elasticity_fe_face_values (elasticity_fe, + common_face_quadrature, + update_values); + FESubfaceValues stokes_fe_subface_values (stokes_fe, + common_face_quadrature, + update_JxW_values | + update_normal_vectors | + update_gradients); + FESubfaceValues 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 local_matrix; + FullMatrix local_interface_matrix (elasticity_dofs_per_cell, + stokes_dofs_per_cell); + Vector local_rhs; + + std::vector local_dof_indices; + std::vector neighbor_dof_indices (stokes_dofs_per_cell); + + const RightHandSide 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 > stokes_phi_grads_u (stokes_dofs_per_cell); + std::vector stokes_div_phi_u (stokes_dofs_per_cell); + std::vector stokes_phi_p (stokes_dofs_per_cell); + + std::vector > elasticity_phi_grad (elasticity_dofs_per_cell); + std::vector elasticity_phi_div (elasticity_dofs_per_cell); + std::vector > 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::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + hp_fe_values.reinit (cell); + + const FEValues &fe_values = hp_fe_values.get_present_fe_values(); + + 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 + // if 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{FluidStructureProblem::assemble_system} + for (unsigned int q=0; q -void FluidStructureProblem::assemble_system () -{ - system_matrix=0; - system_rhs=0; - - const QGauss stokes_quadrature(stokes_degree+2); - const QGauss elasticity_quadrature(elasticity_degree+2); - - hp::QCollection q_collection; - q_collection.push_back (stokes_quadrature); - q_collection.push_back (elasticity_quadrature); - - hp::FEValues hp_fe_values (fe_collection, q_collection, - update_values | - update_quadrature_points | - update_JxW_values | - update_gradients); - - const QGauss common_face_quadrature(std::max (stokes_degree+2, - elasticity_degree+2)); - - FEFaceValues stokes_fe_face_values (stokes_fe, - common_face_quadrature, - update_JxW_values | - update_normal_vectors | - update_gradients); - FEFaceValues elasticity_fe_face_values (elasticity_fe, - common_face_quadrature, - update_values); - FESubfaceValues stokes_fe_subface_values (stokes_fe, - common_face_quadrature, - update_JxW_values | - update_normal_vectors | - update_gradients); - FESubfaceValues 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 local_matrix; - FullMatrix local_interface_matrix (elasticity_dofs_per_cell, - stokes_dofs_per_cell); - Vector local_rhs; - - std::vector local_dof_indices; - std::vector neighbor_dof_indices (stokes_dofs_per_cell); - - const RightHandSide 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 > stokes_phi_grads_u (stokes_dofs_per_cell); - std::vector stokes_div_phi_u (stokes_dofs_per_cell); - std::vector stokes_phi_p (stokes_dofs_per_cell); - - std::vector > elasticity_phi_grad (elasticity_dofs_per_cell); - std::vector elasticity_phi_div (elasticity_dofs_per_cell); - std::vector > 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::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - { - hp_fe_values.reinit (cell); - - const FEValues &fe_values = hp_fe_values.get_present_fe_values(); - - 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 - // if 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; iget_fe().dofs_per_cell; + Assert (dofs_per_cell == elasticity_dofs_per_cell, + ExcInternalError()); - for (unsigned int q=0; qget_fe().dofs_per_cell; - Assert (dofs_per_cell == elasticity_dofs_per_cell, - ExcInternalError()); + for (unsigned int q=0; qlocal_rhs 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::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 f + // 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 - // local_rhs 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::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 f - // 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; - subfaceface(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; + subfaceface(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 -void -FluidStructureProblem:: -assemble_interface_term (const FEFaceValuesBase &elasticity_fe_face_values, - const FEFaceValuesBase &stokes_fe_face_values, - std::vector > &elasticity_phi, - std::vector > &stokes_phi_grads_u, - std::vector &stokes_phi_p, - FullMatrix &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 normal_vector = stokes_fe_face_values.normal_vector(q); - - for (unsigned int k=0; kFluidStructureProblem::solve} + // 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 + void + FluidStructureProblem:: + assemble_interface_term (const FEFaceValuesBase &elasticity_fe_face_values, + const FEFaceValuesBase &stokes_fe_face_values, + std::vector > &elasticity_phi, + std::vector > &stokes_phi_grads_u, + std::vector &stokes_phi_p, + FullMatrix &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 normal_vector = stokes_fe_face_values.normal_vector(q); + + for (unsigned int k=0; k -void -FluidStructureProblem::solve () -{ - SparseDirectUMFPACK direct_solver; - direct_solver.initialize (system_matrix); - direct_solver.vmult (solution, system_rhs); - constraints.distribute (solution); -} + // @sect4{FluidStructureProblem::solve} + + // 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 + void + FluidStructureProblem::solve () + { + SparseDirectUMFPACK direct_solver; + direct_solver.initialize (system_matrix); + direct_solver.vmult (solution, system_rhs); + constraints.distribute (solution); + } - // @sect4{FluidStructureProblem::output_results} - // 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 -void -FluidStructureProblem:: -output_results (const unsigned int refinement_cycle) const -{ - std::vector solution_names (dim, "velocity"); - solution_names.push_back ("pressure"); - for (unsigned int d=0; dFluidStructureProblem::output_results} - std::vector - 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 + void + FluidStructureProblem:: + output_results (const unsigned int refinement_cycle) const + { + std::vector solution_names (dim, "velocity"); + solution_names.push_back ("pressure"); + for (unsigned int d=0; d + 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 > data_out; - data_out.attach_dof_handler (dof_handler); + DataOut > data_out; + data_out.attach_dof_handler (dof_handler); - data_out.add_data_vector (solution, solution_names, - DataOut >::type_dof_data, - data_component_interpretation); - data_out.build_patches (); + data_out.add_data_vector (solution, solution_names, + DataOut >::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{FluidStructureProblem::refine_mesh} - - // 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 -void -FluidStructureProblem::refine_mesh () -{ - Vector - stokes_estimated_error_per_cell (triangulation.n_active_cells()); - Vector - elasticity_estimated_error_per_cell (triangulation.n_active_cells()); - - const QGauss stokes_face_quadrature(stokes_degree+2); - const QGauss elasticity_face_quadrature(elasticity_degree+2); - - hp::QCollection face_q_collection; - face_q_collection.push_back (stokes_face_quadrature); - face_q_collection.push_back (elasticity_face_quadrature); - - std::vector stokes_component_mask (dim+1+dim, false); - for (unsigned int d=0; d::estimate (dof_handler, - face_q_collection, - typename FunctionMap::type(), - solution, - stokes_estimated_error_per_cell, - stokes_component_mask); - - std::vector elasticity_component_mask (dim+1+dim, false); - for (unsigned int d=0; d::estimate (dof_handler, - face_q_collection, - typename FunctionMap::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 - 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 assemble_system. + // @sect4{FluidStructureProblem::refine_mesh} + + // 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 + void + FluidStructureProblem::refine_mesh () { - unsigned int cell_index = 0; - for (typename hp::DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(); - cell != dof_handler.end(); ++cell, ++cell_index) - for (unsigned int f=0; f::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 + stokes_estimated_error_per_cell (triangulation.n_active_cells()); + Vector + elasticity_estimated_error_per_cell (triangulation.n_active_cells()); + + const QGauss stokes_face_quadrature(stokes_degree+2); + const QGauss elasticity_face_quadrature(elasticity_degree+2); + + hp::QCollection face_q_collection; + face_q_collection.push_back (stokes_face_quadrature); + face_q_collection.push_back (elasticity_face_quadrature); + + std::vector stokes_component_mask (dim+1+dim, false); + for (unsigned int d=0; d::estimate (dof_handler, + face_q_collection, + typename FunctionMap::type(), + solution, + stokes_estimated_error_per_cell, + stokes_component_mask); + + std::vector elasticity_component_mask (dim+1+dim, false); + for (unsigned int d=0; d::estimate (dof_handler, + face_q_collection, + typename FunctionMap::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 + 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 assemble_system. + { + unsigned int cell_index = 0; + for (typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + cell != dof_handler.end(); ++cell, ++cell_index) + for (unsigned int f=0; f::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{FluidStructureProblem::run} + // @sect4{FluidStructureProblem::run} - // 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 -void FluidStructureProblem::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 + void FluidStructureProblem::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; + } + } } @@ -1259,6 +1263,9 @@ int main () { try { + using namespace dealii; + using namespace Step46; + deallog.depth_console (0); FluidStructureProblem<2> flow_problem(1, 1); diff --git a/deal.II/examples/step-47/step-47.cc b/deal.II/examples/step-47/step-47.cc index a442fb28ba..0580e15af9 100644 --- a/deal.II/examples/step-47/step-47.cc +++ b/deal.II/examples/step-47/step-47.cc @@ -49,393 +49,395 @@ #include -using namespace dealii; - - - -double sign (double d) -{ - if (d > 0) - return 1; - else if (d < 0) - return -1; - else - return 0; -} - - -template -class LaplaceProblem +namespace Step47 { - public: - LaplaceProblem (); - ~LaplaceProblem (); - - void run (); + using namespace dealii; - private: - bool interface_intersects_cell (const typename Triangulation::cell_iterator &cell) const; - std::pair > compute_quadrature(const Quadrature &plain_quadrature, const typename hp::DoFHandler::active_cell_iterator &cell, const std::vector &level_set_values); - void append_quadrature(const Quadrature &plain_quadrature, - const std::vector > &v , - std::vector > &xfem_points, - std::vector &xfem_weights); - void setup_system (); - void assemble_system (); - void solve (); - void refine_grid (); - void output_results (const unsigned int cycle) const; - void compute_error () const; - Triangulation triangulation; + double sign (double d) + { + if (d > 0) + return 1; + else if (d < 0) + return -1; + else + return 0; + } - hp::DoFHandler dof_handler; - hp::FECollection fe_collection; - ConstraintMatrix constraints; + template + class LaplaceProblem + { + public: + LaplaceProblem (); + ~LaplaceProblem (); - SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; + void run (); - Vector solution; - Vector system_rhs; -}; + private: + bool interface_intersects_cell (const typename Triangulation::cell_iterator &cell) const; + std::pair > compute_quadrature(const Quadrature &plain_quadrature, const typename hp::DoFHandler::active_cell_iterator &cell, const std::vector &level_set_values); + void append_quadrature(const Quadrature &plain_quadrature, + const std::vector > &v , + std::vector > &xfem_points, + std::vector &xfem_weights); + void setup_system (); + void assemble_system (); + void solve (); + void refine_grid (); + void output_results (const unsigned int cycle) const; + void compute_error () const; + Triangulation triangulation; + hp::DoFHandler dof_handler; + hp::FECollection fe_collection; -template -class Coefficient : public Function -{ - public: - Coefficient () : Function() {} + ConstraintMatrix constraints; - virtual double value (const Point &p, - const unsigned int component = 0) const; + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; - virtual void value_list (const std::vector > &points, - std::vector &values, - const unsigned int component = 0) const; -}; + Vector solution; + Vector system_rhs; + }; -template -double Coefficient::value (const Point &p, - const unsigned int) const -{ - if (p.square() < 0.5*0.5) - return 20; - else - return 1; -} + template + class Coefficient : public Function + { + public: + Coefficient () : Function() {} + virtual double value (const Point &p, + const unsigned int component = 0) const; -template -void Coefficient::value_list (const std::vector > &points, - std::vector &values, - const unsigned int component) const -{ - const unsigned int n_points = points.size(); + virtual void value_list (const std::vector > &points, + std::vector &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 + double Coefficient::value (const Point &p, + const unsigned int) const + { + if (p.square() < 0.5*0.5) + return 20; + else + return 1; + } -template -double exact_solution (const Point &p) -{ - const double r = p.norm(); + template + void Coefficient::value_list (const std::vector > &points, + std::vector &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 -LaplaceProblem::LaplaceProblem () - : - dof_handler (triangulation) -{ - fe_collection.push_back (FESystem (FE_Q(1), 1, - FE_Nothing(), 1)); - fe_collection.push_back (FESystem (FE_Q(1), 1, - FE_Q(1), 1)); -} + for (unsigned int i=0; i -LaplaceProblem::~LaplaceProblem () -{ - dof_handler.clear (); -} + template + double exact_solution (const Point &p) + { + const double r = p.norm(); + return (r < 0.5 + ? + 1./20 * (-1./4*r*r + 61./16) + : + 1./4 * (1-r*r)); + } -template -double -level_set (const Point &p) -{ - return p.norm() - 0.5; -} + template + LaplaceProblem::LaplaceProblem () + : + dof_handler (triangulation) + { + fe_collection.push_back (FESystem (FE_Q(1), 1, + FE_Nothing(), 1)); + fe_collection.push_back (FESystem (FE_Q(1), 1, + FE_Q(1), 1)); + } -template -Tensor<1,dim> -grad_level_set (const Point &p) -{ - return p / p.norm(); -} + template + LaplaceProblem::~LaplaceProblem () + { + dof_handler.clear (); + } -template -bool -LaplaceProblem:: -interface_intersects_cell (const typename Triangulation::cell_iterator &cell) const -{ - for (unsigned int v=0; v::vertices_per_cell-1; ++v) - if (level_set(cell->vertex(v)) * level_set(cell->vertex(v+1)) < 0) - return true; - - // we get here only if all vertices - // have the same sign, which means - // that the cell is not intersected - return false; -} + template + double + level_set (const Point &p) + { + return p.norm() - 0.5; + } -template -void LaplaceProblem::setup_system () -{ - for (typename hp::DoFHandler::cell_iterator cell - = dof_handler.begin_active(); - cell != dof_handler.end(); ++cell) - if (interface_intersects_cell(cell) == false) - cell->set_active_fe_index(0); - else - cell->set_active_fe_index(1); + template + Tensor<1,dim> + grad_level_set (const Point &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 + bool + LaplaceProblem:: + interface_intersects_cell (const typename Triangulation::cell_iterator &cell) const + { + for (unsigned int v=0; v::vertices_per_cell-1; ++v) + if (level_set(cell->vertex(v)) * level_set(cell->vertex(v+1)) < 0) + return true; - 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 + void LaplaceProblem::setup_system () + { + for (typename hp::DoFHandler::cell_iterator cell + = dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + if (interface_intersects_cell(cell) == false) + cell->set_active_fe_index(0); + else + cell->set_active_fe_index(1); - 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 -void LaplaceProblem::assemble_system () -{ - const QGauss quadrature_formula(3); +//TODO: component 1 must satisfy zero boundary conditions + constraints.close(); - FEValues plain_fe_values (fe_collection[0], quadrature_formula, - update_values | update_gradients | - update_quadrature_points | update_JxW_values); - FEValues enriched_fe_values (fe_collection[1], quadrature_formula, - update_values | update_gradients | - update_quadrature_points | update_JxW_values); + 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 cell_matrix; - Vector cell_rhs; + sparsity_pattern.copy_from(c_sparsity); - std::vector local_dof_indices; + system_matrix.reinit (sparsity_pattern); + } - const Coefficient coefficient; - std::vector coefficient_values (n_q_points); - typename hp::DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); + template + void LaplaceProblem::assemble_system () + { + const QGauss 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 plain_fe_values (fe_collection[0], quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + FEValues enriched_fe_values (fe_collection[1], quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); - 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 cell_matrix; + Vector cell_rhs; - for (unsigned int q_point=0; q_point local_dof_indices; + const Coefficient coefficient; + std::vector 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::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); - std::vector level_set_values (GeometryInfo::vertices_per_cell); - for (unsigned int v=0; v::vertices_per_cell; ++v) - level_set_values[v] = level_set (cell->vertex(v)); + 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 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_pointget_fe().system_to_component_index(i).first == 0) - { - for (unsigned int j=0; jget_fe().system_to_component_index(j).first == 0) - cell_matrix(i,j) += (coefficient_values[q_point] * - this_fe_values.shape_grad(i,q_point) * - this_fe_values.shape_grad(j,q_point) * - this_fe_values.JxW(q_point)); - else - cell_matrix(i,j) += (coefficient_values[q_point] * - this_fe_values.shape_grad(i,q_point) - * - ((std::fabs(level_set(this_fe_values.quadrature_point(q_point))) - - - std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(j).second))))* - this_fe_values.shape_grad(j,q_point) - + - grad_level_set(this_fe_values.quadrature_point(q_point)) * - sign(level_set(this_fe_values.quadrature_point(q_point))) * - this_fe_values.shape_value(j,q_point)) * - this_fe_values.JxW(q_point)); - - cell_rhs(i) += (this_fe_values.shape_value(i,q_point) * - 1.0 * - this_fe_values.JxW(q_point)); - } - else + for (unsigned int q_point=0; q_pointget_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 level_set_values (GeometryInfo::vertices_per_cell); + for (unsigned int v=0; v::vertices_per_cell; ++v) + level_set_values[v] = level_set (cell->vertex(v)); + + FEValues this_fe_values (fe_collection[1], + compute_quadrature(quadrature_formula, cell, + level_set_values).second, + update_values | update_gradients | + update_quadrature_points | update_JxW_values ); + + this_fe_values.reinit (cell); + + coefficient_values.resize (this_fe_values.n_quadrature_points); + coefficient.value_list (this_fe_values.get_quadrature_points(), + coefficient_values); + + for (unsigned int q_point=0; q_pointget_fe().system_to_component_index(i).first == 0) + { + for (unsigned int j=0; jget_fe().system_to_component_index(j).first == 0) + cell_matrix(i,j) += (coefficient_values[q_point] * + this_fe_values.shape_grad(i,q_point) * + this_fe_values.shape_grad(j,q_point) * + this_fe_values.JxW(q_point)); + else + cell_matrix(i,j) += (coefficient_values[q_point] * + this_fe_values.shape_grad(i,q_point) + * + ((std::fabs(level_set(this_fe_values.quadrature_point(q_point))) + - + std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(j).second))))* + this_fe_values.shape_grad(j,q_point) + + + grad_level_set(this_fe_values.quadrature_point(q_point)) * + sign(level_set(this_fe_values.quadrature_point(q_point))) * + this_fe_values.shape_value(j,q_point)) * + this_fe_values.JxW(q_point)); + + cell_rhs(i) += (this_fe_values.shape_value(i,q_point) * + 1.0 * + this_fe_values.JxW(q_point)); + } + else + { + for (unsigned int j=0; jget_fe().system_to_component_index(j).first == 0) + cell_matrix(i,j) += (coefficient_values[q_point] * + ((std::fabs(level_set(this_fe_values.quadrature_point(q_point))) + - + std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(i).second))))* + this_fe_values.shape_grad(i,q_point) + + + grad_level_set(this_fe_values.quadrature_point(q_point)) * + sign(level_set(this_fe_values.quadrature_point(q_point))) * + this_fe_values.shape_value(i,q_point)) * + this_fe_values.shape_grad(j,q_point) * + this_fe_values.JxW(q_point)); + else + cell_matrix(i,j) += (coefficient_values[q_point] * + ((std::fabs(level_set(this_fe_values.quadrature_point(q_point))) + - + std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(i).second))))* + this_fe_values.shape_grad(i,q_point) + + + grad_level_set(this_fe_values.quadrature_point(q_point)) * + sign(level_set(this_fe_values.quadrature_point(q_point))) * + this_fe_values.shape_value(i,q_point)) * + ((std::fabs(level_set(this_fe_values.quadrature_point(q_point))) + - + std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(j).second))))* + this_fe_values.shape_grad(j,q_point) + + + grad_level_set(this_fe_values.quadrature_point(q_point)) * + sign(level_set(this_fe_values.quadrature_point(q_point))) * + this_fe_values.shape_value(j,q_point)) * + this_fe_values.JxW(q_point)); + + cell_rhs(i) += ((std::fabs(level_set(this_fe_values.quadrature_point(q_point))) + - + std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(i).second))))* + this_fe_values.shape_value(i,q_point) * + 1.0 * + this_fe_values.JxW(q_point)); + } + } - local_dof_indices.resize (dofs_per_cell); - cell->get_dof_indices (local_dof_indices); - constraints.distribute_local_to_global (cell_matrix, cell_rhs, - local_dof_indices, - system_matrix, system_rhs); - } + local_dof_indices.resize (dofs_per_cell); + cell->get_dof_indices (local_dof_indices); + constraints.distribute_local_to_global (cell_matrix, cell_rhs, + local_dof_indices, + system_matrix, system_rhs); + } - std::map boundary_values; - VectorTools::interpolate_boundary_values (dof_handler, - 0, - ZeroFunction(2), - boundary_values); - MatrixTools::apply_boundary_values (boundary_values, - system_matrix, - solution, - system_rhs); + std::map boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction(2), + boundary_values); + MatrixTools::apply_boundary_values (boundary_values, + system_matrix, + solution, + system_rhs); -} + } // To integrate the enriched elements we have to find the geometrical decomposition // of the original element in subelements. The subelements are used to integrate @@ -447,654 +449,655 @@ void LaplaceProblem::assemble_system () // are considered. // Type 1: there is not cut. Type 2: a corner of the element is cut. Type 3: two corners are cut. -template -std::pair > -LaplaceProblem::compute_quadrature (const Quadrature &plain_quadrature, - const typename hp::DoFHandler::active_cell_iterator &cell, - const std::vector &level_set_values ) -{ - - unsigned int type = 0; - - // find the type of cut - int sign_ls[GeometryInfo::vertices_per_cell]; - for (unsigned int v=0; v::vertices_per_cell; ++v) - { - if (level_set_values[v] > 0) sign_ls[v] = 1; - else if (level_set_values[v] < 0) sign_ls[v] = -1; - else sign_ls[v] = 0; - } + template + std::pair > + LaplaceProblem::compute_quadrature (const Quadrature &plain_quadrature, + const typename hp::DoFHandler::active_cell_iterator &cell, + const std::vector &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::vertices_per_cell]; + for (unsigned int v=0; v::vertices_per_cell; ++v) + { + if (level_set_values[v] > 0) sign_ls[v] = 1; + else if (level_set_values[v] < 0) sign_ls[v] = -1; + else sign_ls[v] = 0; + } - 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 v0(0,0); - Point v1(1,0); - Point v2(0,1); - Point 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 A(0,0); - Point B(0,0); - Point C(0,0); - Point D(0,0); - Point E(0,0); - Point F(0,0); - - if (type == 1) - return std::pair >(1, plain_quadrature); - - if (type==2) - { - const unsigned int n_q_points = plain_quadrature.size(); - - // loop over all subelements for integration - // in type 2 there are 5 subelements - - Quadrature xfem_quadrature(5*n_q_points); - - std::vector > v(GeometryInfo::vertices_per_cell); - - if (sign_ls[0]!=sign_ls[1] && sign_ls[0]!=sign_ls[2] && sign_ls[0]!=sign_ls[3]) Pos = 0; - else if (sign_ls[1]!=sign_ls[0] && sign_ls[1]!=sign_ls[2] && sign_ls[1]!=sign_ls[3]) Pos = 1; - else if (sign_ls[2]!=sign_ls[0] && sign_ls[2]!=sign_ls[1] && sign_ls[2]!=sign_ls[3]) Pos = 2; - else if (sign_ls[3]!=sign_ls[0] && sign_ls[3]!=sign_ls[1] && sign_ls[3]!=sign_ls[2]) Pos = 3; - else assert(0); // error message - - // Find cut coordinates - - // deal.ii local coordinates - - // 2-------3 - // | | - // | | - // | | - // 0-------1 - - if (Pos == 0) - { - A[0] = 1. - level_set_values[1]/(level_set_values[1]-level_set_values[0]); - B[1] = 1. - level_set_values[2]/(level_set_values[2]-level_set_values[0]); - A(1) = 0.; - B(0) = 0.; - C(0) = 0.5*( A(0) + B(0) ); - C(1) = 0.5*( A(1) + B(1) ); - D(0) = 2./3. * C(0); - D(1) = 2./3. * C(1); - E(0) = 0.5*A(0); - E(1) = 0.; - F(0) = 0.; - F(1) = 0.5*B(1); - } - else if (Pos == 1) - { - A[0] = level_set_values[0]/(level_set_values[0]-level_set_values[1]); - B[1] = 1 - level_set_values[3]/(level_set_values[3]-level_set_values[1]); - A(1) = 0.; - B(0) = 1.; - C(0) = 0.5*( A(0) + B(0) ); - C(1) = 0.5*( A(1) + B(1) ); - D(0) = 1./3. + 2./3. * C(0); - D(1) = 2./3. * C(1); - E(0) = 0.5*(1 + A(0)); - E(1) = 0.; - F(0) = 1.; - F(1) = 0.5*B(1); - } - else if (Pos == 2) - { - A[0] = 1 - level_set_values[3]/(level_set_values[3]-level_set_values[2]); - B[1] = level_set_values[0]/(level_set_values[0]-level_set_values[2]); - A(1) = 1.; - B(0) = 0.; - C(0) = 0.5*( A(0) + B(0) ); - C(1) = 0.5*( A(1) + B(1) ); - D(0) = 2./3. * C(0); - D(1) = 1./3. + 2./3. * C(1); - E(0) = 0.5* A(0); - E(1) = 1.; - F(0) = 0.; - F(1) = 0.5*( 1. + B(1) ); - } - else if (Pos == 3) - { - A[0] = level_set_values[2]/(level_set_values[2]-level_set_values[3]); - B[1] = level_set_values[1]/(level_set_values[1]-level_set_values[3]); - A(1) = 1.; - B(0) = 1.; - C(0) = 0.5*( A(0) + B(0) ); - C(1) = 0.5*( A(1) + B(1) ); - D(0) = 1./3. + 2./3. * C(0); - D(1) = 1./3. + 2./3. * C(1); - E(0) = 0.5*( 1. + A(0) ); - E(1) = 1.; - F(0) = 1.; - F(1) = 0.5*( 1. + B(1) ); - } - - //std::cout << A << std::endl; - //std::cout << B << std::endl; - //std::cout << C << std::endl; - //std::cout << D << std::endl; - //std::cout << E << std::endl; - //std::cout << F << std::endl; - - std::string filename = "vertices.dat"; - std::ofstream output (filename.c_str()); - output << "#vertices of xfem subcells" << std::endl; - output << v0(0) << " " << v0(1) << std::endl; - output << v1(0) << " " << v1(1) << std::endl; - output << v3(0) << " " << v3(1) << std::endl; - output << v2(0) << " " << v2(1) << std::endl; - output << std::endl; - output << A(0) << " " << A(1) << std::endl; - output << B(0) << " " << B(1) << std::endl; - output << std::endl; - output << C(0) << " " << C(1) << std::endl; - output << D(0) << " " << D(1) << std::endl; - output << std::endl; - output << D(0) << " " << D(1) << std::endl; - output << E(0) << " " << E(1) << std::endl; - output << std::endl; - output << D(0) << " " << D(1) << std::endl; - output << F(0) << " " << F(1) << std::endl; - output << std::endl; - - if (Pos==0) - output << v3(0) << " " << v3(1) << std::endl; - else if (Pos==1) - output << v2(0) << " " << v2(1) << std::endl; - else if (Pos==2) - output << v1(0) << " " << v1(1) << std::endl; - else if (Pos==3) - output << v0(0) << " " << v0(1) << std::endl; - output << C(0) << " " << C(1) << std::endl; - - Point subcell_vertices[10]; - subcell_vertices[0] = v0; - subcell_vertices[1] = v1; - subcell_vertices[2] = v2; - subcell_vertices[3] = v3; - subcell_vertices[4] = A; - subcell_vertices[5] = B; - subcell_vertices[6] = C; - subcell_vertices[7] = D; - subcell_vertices[8] = E; - subcell_vertices[9] = F; - - std::vector > xfem_points; - std::vector xfem_weights; - - // lookup table for the decomposition - - if (dim==2) - { - unsigned int subcell_v_indices[4][5][4] = { - {{0,8,9,7}, {9,7,5,6}, {8,4,7,6}, {5,6,2,3}, {6,4,3,1}}, - {{8,1,7,9}, {4,8,6,7}, {6,7,5,9}, {0,4,2,6}, {2,6,3,5}}, - {{9,7,2,8}, {5,6,9,7}, {6,4,7,8}, {0,1,5,6}, {6,1,4,3}}, - {{7,9,8,3}, {4,6,8,7}, {6,5,7,9}, {0,6,2,4}, {0,1,6,5}} - }; - - for (unsigned int subcell = 0; subcell<5; subcell++) - { - //std::cout << "subcell : " << subcell << std::endl; - std::vector > vertices; - for (unsigned int i=0; i<4; i++) - { - vertices.push_back( subcell_vertices[subcell_v_indices[Pos][subcell][i]] ); - //std::cout << "i : " << i << std::endl; - //std::cout << "subcell v : " << subcell_v_indices[Pos][subcell][i] << std::endl; - //std::cout << vertices[i](0) << " " << vertices[i](1) << std::endl; - } - //std::cout << std::endl; - // create quadrature rule - append_quadrature( plain_quadrature, - vertices, - xfem_points, - xfem_weights); - //initialize xfem_quadrature with quadrature points of all subelements - xfem_quadrature.initialize(xfem_points, xfem_weights); - } - } + unsigned int Pos = 100; - Assert (xfem_quadrature.size() == plain_quadrature.size() * 5, ExcInternalError()); - return std::pair >(2, xfem_quadrature); - } + Point v0(0,0); + Point v1(1,0); + Point v2(0,1); + Point v3(1,1); - // Type three decomposition - // (+--+, ++--, +-+-, -++-, --++, -+-+) + Point A(0,0); + Point B(0,0); + Point C(0,0); + Point D(0,0); + Point E(0,0); + Point F(0,0); - if (type==3) - { - const unsigned int n_q_points = plain_quadrature.size(); - - // loop over all subelements for integration - // in type 2 there are 5 subelements - - Quadrature xfem_quadrature(5*n_q_points); - - std::vector > v(GeometryInfo::vertices_per_cell); - - if ( sign_ls[0]==sign_ls[1] && sign_ls[2]==sign_ls[3] ) - { - Pos = 0; - A(0) = 0.; - A(1) = level_set_values[0]/((level_set_values[0]-level_set_values[2])); - B(0) = 1.; - B(1) = level_set_values[1]/((level_set_values[1]-level_set_values[3])); - } - else if ( sign_ls[0]==sign_ls[2] && sign_ls[1]==sign_ls[3] ) - { - Pos = 1; - A(0) = level_set_values[0]/((level_set_values[0]-level_set_values[1])); - A(1) = 0.; - B(0) = level_set_values[2]/((level_set_values[2]-level_set_values[3])); - B(1) = 1.; - } - else if ( sign_ls[0]==sign_ls[3] && sign_ls[1]==sign_ls[2] ) - { - std::cout << "Error: the element has two cut lines and this is not allowed" << std::endl; - assert(0); - } - else - { - std::cout << "Error: the level set function has not the right values" << std::endl; - assert(0); - } - - //std::cout << "Pos " << Pos << std::endl; - //std::cout << A << std::endl; - //std::cout << B << std::endl; - std::string filename = "vertices.dat"; - std::ofstream output (filename.c_str()); - output << "#vertices of xfem subcells" << std::endl; - output << A(0) << " " << A(1) << std::endl; - output << B(0) << " " << B(1) << std::endl; - - //fill xfem_quadrature - Point subcell_vertices[6]; - subcell_vertices[0] = v0; - subcell_vertices[1] = v1; - subcell_vertices[2] = v2; - subcell_vertices[3] = v3; - subcell_vertices[4] = A; - subcell_vertices[5] = B; - - std::vector > xfem_points; - std::vector xfem_weights; - - if (dim==2) - { - unsigned int subcell_v_indices[2][2][4] = { - {{0,1,4,5}, {4,5,2,3}}, - {{0,4,2,5}, {4,1,5,3}} - }; - - //std::cout << "Pos : " << Pos << std::endl; - for (unsigned int subcell = 0; subcell<2; subcell++) - { - //std::cout << "subcell : " << subcell << std::endl; - std::vector > vertices; - for (unsigned int i=0; i<4; i++) - { - vertices.push_back( subcell_vertices[subcell_v_indices[Pos][subcell][i]] ); - //std::cout << "i : " << i << std::endl; - //std::cout << "subcell v : " << subcell_v_indices[Pos][subcell][i] << std::endl; - //std::cout << vertices[i](0) << " " << vertices[i](1) << std::endl; - } - //std::cout << std::endl; - // create quadrature rule - append_quadrature( plain_quadrature, - vertices, - xfem_points, - xfem_weights); - //initialize xfem_quadrature with quadrature points of all subelements - xfem_quadrature.initialize(xfem_points, xfem_weights); - } - } - Assert (xfem_quadrature.size() == plain_quadrature.size() * 2, ExcInternalError()); - return std::pair >(3, xfem_quadrature); - } + if (type == 1) + return std::pair >(1, plain_quadrature); - return std::pair >(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 -void LaplaceProblem::append_quadrature ( const Quadrature &plain_quadrature, - const std::vector > &v, - std::vector > &xfem_points, - std::vector &xfem_weights) + Quadrature 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 > v(GeometryInfo::vertices_per_cell); - unsigned int n_v = GeometryInfo::vertices_per_cell; + if (sign_ls[0]!=sign_ls[1] && sign_ls[0]!=sign_ls[2] && sign_ls[0]!=sign_ls[3]) Pos = 0; + else if (sign_ls[1]!=sign_ls[0] && sign_ls[1]!=sign_ls[2] && sign_ls[1]!=sign_ls[3]) Pos = 1; + else if (sign_ls[2]!=sign_ls[0] && sign_ls[2]!=sign_ls[1] && sign_ls[2]!=sign_ls[3]) Pos = 2; + else if (sign_ls[3]!=sign_ls[0] && sign_ls[3]!=sign_ls[1] && sign_ls[3]!=sign_ls[2]) Pos = 3; + else assert(0); // error message - std::vector > q_points = plain_quadrature.get_points(); - std::vector > q_transf(q_points.size()); - std::vector W = plain_quadrature.get_weights(); - std::vector phi(n_v); - std::vector > grad_phi(n_v); + // Find cut coordinates - const unsigned int n_q_points = plain_quadrature.size(); + // deal.ii local coordinates - std::vector 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 subcell_vertices[10]; + subcell_vertices[0] = v0; + subcell_vertices[1] = v1; + subcell_vertices[2] = v2; + subcell_vertices[3] = v3; + subcell_vertices[4] = A; + subcell_vertices[5] = B; + subcell_vertices[6] = C; + subcell_vertices[7] = D; + subcell_vertices[8] = E; + subcell_vertices[9] = F; + + std::vector > xfem_points; + std::vector xfem_weights; + + // lookup table for the decomposition + + if (dim==2) { - for (unsigned int j = 0; j::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 > 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 >(2, xfem_quadrature); + } - // Map integration points from reference element to subcell of reference element - Point q_prime; - for (unsigned int d=0; d::vertices_per_cell; j++) - q_prime[d] += v[j](d) * phi[j]; - xfem_points.push_back(q_prime); - } + // 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 -void LaplaceProblem::solve () -{ - SolverControl solver_control (1000, 1e-12); - SolverCG<> solver (solver_control); + Quadrature xfem_quadrature(5*n_q_points); + + std::vector > v(GeometryInfo::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 subcell_vertices[6]; + subcell_vertices[0] = v0; + subcell_vertices[1] = v1; + subcell_vertices[2] = v2; + subcell_vertices[3] = v3; + subcell_vertices[4] = A; + subcell_vertices[5] = B; + + std::vector > xfem_points; + std::vector xfem_weights; + + if (dim==2) + { + unsigned int subcell_v_indices[2][2][4] = { + {{0,1,4,5}, {4,5,2,3}}, + {{0,4,2,5}, {4,1,5,3}} + }; - 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 > vertices; + for (unsigned int i=0; i<4; i++) + { + vertices.push_back( subcell_vertices[subcell_v_indices[Pos][subcell][i]] ); + //std::cout << "i : " << i << std::endl; + //std::cout << "subcell v : " << subcell_v_indices[Pos][subcell][i] << std::endl; + //std::cout << vertices[i](0) << " " << vertices[i](1) << std::endl; + } + //std::cout << std::endl; + // create quadrature rule + append_quadrature( plain_quadrature, + vertices, + xfem_points, + xfem_weights); + //initialize xfem_quadrature with quadrature points of all subelements + xfem_quadrature.initialize(xfem_points, xfem_weights); + } + } + Assert (xfem_quadrature.size() == plain_quadrature.size() * 2, ExcInternalError()); + return std::pair >(3, xfem_quadrature); + } + return std::pair >(0, plain_quadrature);; + } -template -void LaplaceProblem::refine_grid () -{ - Vector estimated_error_per_cell (triangulation.n_active_cells()); + template + void LaplaceProblem::append_quadrature ( const Quadrature &plain_quadrature, + const std::vector > &v, + std::vector > &xfem_points, + std::vector &xfem_weights) - KellyErrorEstimator::estimate (dof_handler, - QGauss(3), - typename FunctionMap::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::vertices_per_cell; - triangulation.execute_coarsening_and_refinement (); -} + std::vector > q_points = plain_quadrature.get_points(); + std::vector > q_transf(q_points.size()); + std::vector W = plain_quadrature.get_weights(); + std::vector phi(n_v); + std::vector > grad_phi(n_v); + const unsigned int n_q_points = plain_quadrature.size(); + std::vector JxW(n_q_points); -template -class Postprocessor : public DataPostprocessor -{ - public: - virtual - void - compute_derived_quantities_vector (const std::vector > &uh, - const std::vector > > &duh, - const std::vector > > &dduh, - const std::vector > &normals, - const std::vector > &evaluation_points, - std::vector > &computed_quantities) const; + 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 get_names () const; + default: + Assert (false, ExcNotImplemented()); + } - virtual unsigned int n_output_variables() const; - virtual - std::vector - 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::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 -std::vector -Postprocessor::get_names() const -{ - std::vector 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 q_prime; + for (unsigned int d=0; d::vertices_per_cell; j++) + q_prime[d] += v[j](d) * phi[j]; + xfem_points.push_back(q_prime); + } + } -template -unsigned int -Postprocessor::n_output_variables() const -{ - return get_names().size(); -} + template + void LaplaceProblem::solve () + { + SolverControl solver_control (1000, 1e-12); + SolverCG<> solver (solver_control); -template -std::vector -Postprocessor:: -get_data_component_interpretation () const -{ - std::vector - 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 -UpdateFlags -Postprocessor::get_needed_update_flags() const -{ - return update_values | update_q_points; -} + constraints.distribute (solution); + } -template -void -Postprocessor:: -compute_derived_quantities_vector (const std::vector > &uh, - const std::vector > > &/*duh*/, - const std::vector > > &/*dduh*/, - const std::vector > &/*normals*/, - const std::vector > &evaluation_points, - std::vector > &computed_quantities) const -{ - const unsigned int n_quadrature_points = uh.size(); - Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError()); - Assert (uh[0].size() == 2, ExcInternalError()); - Assert (computed_quantities[0].size()==n_output_variables(),ExcInternalError()); - for (unsigned int q=0; q + void LaplaceProblem::refine_grid () + { + Vector estimated_error_per_cell (triangulation.n_active_cells()); + KellyErrorEstimator::estimate (dof_handler, + QGauss(3), + typename FunctionMap::type(), + solution, + estimated_error_per_cell); + GridRefinement::refine_and_coarsen_fixed_number (triangulation, + estimated_error_per_cell, + 0.3, 0.03); -template -void LaplaceProblem::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 postprocessor; - DataOut > data_out; + template + class Postprocessor : public DataPostprocessor + { + public: + virtual + void + compute_derived_quantities_vector (const std::vector > &uh, + const std::vector > > &duh, + const std::vector > > &dduh, + const std::vector > &normals, + const std::vector > &evaluation_points, + std::vector > &computed_quantities) const; - 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 get_names () const; - //data_out.write_vtk (output); - data_out.write_gmv (output); -} + virtual unsigned int n_output_variables() const; + virtual + std::vector + get_data_component_interpretation () const; + virtual UpdateFlags get_needed_update_flags () const; + }; -template -void LaplaceProblem::compute_error () const -{ - hp::QCollection q_collection; - q_collection.push_back (QGauss(2)); - q_collection.push_back (QIterated(QGauss<1>(2), 4)); - hp::FEValues hp_fe_values (fe_collection, q_collection, - update_values | update_q_points | update_JxW_values); + template + std::vector + Postprocessor::get_names() const + { + std::vector solution_names (1, "total_solution"); + solution_names.push_back ("error"); + return solution_names; + } - double l2_error_square = 0; - std::vector > solution_values; + template + unsigned int + Postprocessor::n_output_variables() const + { + return get_names().size(); + } - typename hp::DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - { - hp_fe_values.reinit (cell); - - const FEValues &fe_values = hp_fe_values.get_present_fe_values (); - - solution_values.resize (fe_values.n_quadrature_points, - Vector(2)); - fe_values.get_function_values (solution, - solution_values); - - for (unsigned int q=0; q + std::vector + Postprocessor:: + get_data_component_interpretation () const + { + std::vector + interpretation (2, + DataComponentInterpretation::component_is_scalar); + return interpretation; + } + + + template + UpdateFlags + Postprocessor::get_needed_update_flags() const + { + return update_values | update_q_points; + } + + + template + void + Postprocessor:: + compute_derived_quantities_vector (const std::vector > &uh, + const std::vector > > &/*duh*/, + const std::vector > > &/*dduh*/, + const std::vector > &/*normals*/, + const std::vector > &evaluation_points, + std::vector > &computed_quantities) const + { + const unsigned int n_quadrature_points = uh.size(); + Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError()); + Assert (uh[0].size() == 2, ExcInternalError()); + Assert (computed_quantities[0].size()==n_output_variables(),ExcInternalError()); + + for (unsigned int q=0; q + void LaplaceProblem::output_results (const unsigned int cycle) const + { + Assert (cycle < 10, ExcNotImplemented()); + + std::string filename = "solution-"; + filename += ('0' + cycle); + //filename += ".vtk"; + filename += ".gmv"; + + std::ofstream output (filename.c_str()); + + Postprocessor postprocessor; + DataOut > data_out; + + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, "solution"); + data_out.add_data_vector (solution, postprocessor); + data_out.build_patches (5); + + //data_out.write_vtk (output); + data_out.write_gmv (output); + } + + + + template + void LaplaceProblem::compute_error () const + { + hp::QCollection q_collection; + q_collection.push_back (QGauss(2)); + q_collection.push_back (QIterated(QGauss<1>(2), 4)); + + hp::FEValues hp_fe_values (fe_collection, q_collection, + update_values | update_q_points | update_JxW_values); + + double l2_error_square = 0; + + std::vector > solution_values; + + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + for (; cell!=endc; ++cell) + { + hp_fe_values.reinit (cell); + + const FEValues &fe_values = hp_fe_values.get_present_fe_values (); + + solution_values.resize (fe_values.n_quadrature_points, + Vector(2)); + fe_values.get_function_values (solution, + solution_values); + + for (unsigned int q=0; q + void LaplaceProblem::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 boundary; - triangulation.set_boundary (0, boundary); + static const HyperBallBoundary 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); + } + } } @@ -1104,6 +1107,9 @@ int main () try { + using namespace dealii; + using namespace Step47; + deallog.depth_console (0); LaplaceProblem<2> laplace_problem_2d; -- 2.39.5