From: bangerth Date: Fri, 9 Sep 2011 04:12:49 +0000 (+0000) Subject: Move things local to each program into a local namespace. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c22c3fff7d6217a955c6bf2fc0bb83563123f1e7;p=dealii-svn.git Move things local to each program into a local namespace. git-svn-id: https://svn.dealii.org/trunk@24297 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-42/step-42.cc b/deal.II/examples/step-42/step-42.cc index b195fa35e1..2431fec161 100644 --- a/deal.II/examples/step-42/step-42.cc +++ b/deal.II/examples/step-42/step-42.cc @@ -3,7 +3,7 @@ /* $Id$ */ /* */ -/* Copyright (C) 2008, 2009, 2010 by the deal.II authors */ +/* Copyright (C) 2008, 2009, 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 */ @@ -71,1076 +71,1079 @@ #include -using namespace dealii; - +namespace Step42 +{ + using namespace dealii; -template -struct InnerPreconditioner; + template + struct InnerPreconditioner; -template <> -struct InnerPreconditioner<2> -{ - typedef SparseDirectUMFPACK type; -}; + template <> + struct InnerPreconditioner<2> + { + typedef SparseDirectUMFPACK type; + }; -template <> -struct InnerPreconditioner<3> -{ - typedef SparseILU type; -}; -template -void copy(const MATRIX &matrix, - FullMatrix &full_matrix) -{ - const unsigned int m = matrix.m(); - const unsigned int n = matrix.n(); - full_matrix.reinit(n,m); + template <> + struct InnerPreconditioner<3> + { + typedef SparseILU type; + }; - Vector unit (n); - Vector result (m); - for(unsigned int i=0; i + void copy(const MATRIX &matrix, + FullMatrix &full_matrix) { - unit(i) = 1; - for(unsigned int j=0; j unit (n); + Vector result (m); + for(unsigned int i=0; i -class StokesProblem -{ - public: - StokesProblem (const unsigned int degree); - void run (); + template + class StokesProblem + { + public: + StokesProblem (const unsigned int degree); + void run (); - private: - void setup_dofs (); - void assemble_system (); - void assemble_multigrid (); - void solve (); - void solve_block (); + private: + void setup_dofs (); + void assemble_system (); + void assemble_multigrid (); + void solve (); + void solve_block (); - void find_dofs_on_lower_level (std::vector > &lower_dofs, - std::vector > &boundary_dofs); + void find_dofs_on_lower_level (std::vector > &lower_dofs, + std::vector > &boundary_dofs); - void output_results (const unsigned int refinement_cycle) const; - void refine_mesh (); + void output_results (const unsigned int refinement_cycle) const; + void refine_mesh (); - const unsigned int degree; + const unsigned int degree; - Triangulation triangulation; - FESystem fe; - MGDoFHandler dof_handler; + Triangulation triangulation; + FESystem fe; + MGDoFHandler dof_handler; - ConstraintMatrix constraints; + ConstraintMatrix constraints; - BlockSparsityPattern sparsity_pattern; - BlockSparseMatrix system_matrix; + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_matrix; - BlockVector solution; - BlockVector system_rhs; + BlockVector solution; + BlockVector system_rhs; - MGLevelObject mg_constraints; - MGLevelObject mg_sparsity; - MGLevelObject > mg_matrices; + MGLevelObject mg_constraints; + MGLevelObject mg_sparsity; + MGLevelObject > mg_matrices; - MGLevelObject > mg_interface_matrices; - MGConstrainedDoFs mg_constrained_dofs; - std::vector > mg_dofs_per_component; + MGLevelObject > mg_interface_matrices; + MGConstrainedDoFs mg_constrained_dofs; + std::vector > mg_dofs_per_component; - std::vector::type> > mg_A_preconditioner; - std_cxx1x::shared_ptr::type> A_preconditioner; -}; + std::vector::type> > mg_A_preconditioner; + std_cxx1x::shared_ptr::type> A_preconditioner; + }; -template -class BoundaryValues : public Function -{ - public: - BoundaryValues () : Function(dim+1) {} + template + class BoundaryValues : public Function + { + public: + BoundaryValues () : Function(dim+1) {} - virtual double value (const Point &p, - const unsigned int component = 0) const; + virtual double value (const Point &p, + const unsigned int component = 0) const; - virtual void vector_value (const Point &p, - Vector &value) const; -}; + virtual void vector_value (const Point &p, + Vector &value) const; + }; -template -double -BoundaryValues::value (const Point &p, - const unsigned int component) const -{ - Assert (component < this->n_components, - ExcIndexRange (component, 0, this->n_components)); + template + double + BoundaryValues::value (const Point &p, + const unsigned int component) const + { + Assert (component < this->n_components, + ExcIndexRange (component, 0, this->n_components)); - if (component == 0 && p[0] == 0) - return (dim == 2 ? - p[1]*(p[1]-1.) : p[1]*(p[1]-1.) * p[2]*(p[2]-1.)); - return 0; -} + if (component == 0 && p[0] == 0) + return (dim == 2 ? - p[1]*(p[1]-1.) : p[1]*(p[1]-1.) * p[2]*(p[2]-1.)); + return 0; + } -template -void -BoundaryValues::vector_value (const Point &p, - Vector &values) const -{ - for (unsigned int c=0; cn_components; ++c) - values(c) = BoundaryValues::value (p, c); -} + template + void + BoundaryValues::vector_value (const Point &p, + Vector &values) const + { + for (unsigned int c=0; cn_components; ++c) + values(c) = BoundaryValues::value (p, c); + } -template -class RightHandSide : public Function -{ - public: - RightHandSide () : Function(dim+1) {} + template + class RightHandSide : public Function + { + public: + RightHandSide () : Function(dim+1) {} - virtual double value (const Point &p, - const unsigned int component = 0) const; + virtual double value (const Point &p, + const unsigned int component = 0) const; - virtual void vector_value (const Point &p, - Vector &value) const; + virtual void vector_value (const Point &p, + Vector &value) const; -}; + }; -template -double -RightHandSide::value (const Point &/*p*/, - const unsigned int component) const -{ - return (component == 1 ? 1 : 0); -} + template + double + RightHandSide::value (const Point &/*p*/, + const unsigned int component) const + { + return (component == 1 ? 1 : 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 + 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 -class InverseMatrix : public Subscriptor -{ - public: - InverseMatrix (const Matrix &m, - const Preconditioner &preconditioner); + template + class InverseMatrix : public Subscriptor + { + public: + InverseMatrix (const Matrix &m, + const Preconditioner &preconditioner); - void vmult (Vector &dst, - const Vector &src) const; + void vmult (Vector &dst, + const Vector &src) const; - mutable std::string name; - private: - const SmartPointer matrix; - const SmartPointer preconditioner; -}; + mutable std::string name; + private: + const SmartPointer matrix; + const SmartPointer preconditioner; + }; -template -InverseMatrix::InverseMatrix (const Matrix &m, - const Preconditioner &preconditioner) - : - matrix (&m), - preconditioner (&preconditioner) -{} + template + InverseMatrix::InverseMatrix (const Matrix &m, + const Preconditioner &preconditioner) + : + matrix (&m), + preconditioner (&preconditioner) + {} -template -void InverseMatrix::vmult (Vector &dst, - const Vector &src) const -{ - SolverControl solver_control (src.size(), 1.0e-12*src.l2_norm()); - SolverCG<> cg (solver_control); + template + void InverseMatrix::vmult (Vector &dst, + const Vector &src) const + { + SolverControl solver_control (src.size(), 1.0e-12*src.l2_norm()); + SolverCG<> cg (solver_control); - dst = 0; + dst = 0; - try - { - cg.solve (*matrix, dst, src, *preconditioner); - } - catch (...) - { - std::cout << "Failure in " << __PRETTY_FUNCTION__ << std::endl; - abort (); - } + try + { + cg.solve (*matrix, dst, src, *preconditioner); + } + catch (...) + { + std::cout << "Failure in " << __PRETTY_FUNCTION__ << std::endl; + abort (); + } #ifdef STEP_42_TEST - if (name == "in schur") - std::cout << " " << solver_control.last_step() - << " inner CG steps inside the Schur complement "; - else if (name == "top left") - std::cout << " " << solver_control.last_step() - << " CG steps on the top left block "; - else if (name == "rhs") - std::cout << " " << solver_control.last_step() - << " CG steps for computing the r.h.s. "; - else - abort (); - - std::cout << solver_control.initial_value() << "->" << solver_control.last_value() - << std::endl; + if (name == "in schur") + std::cout << " " << solver_control.last_step() + << " inner CG steps inside the Schur complement "; + else if (name == "top left") + std::cout << " " << solver_control.last_step() + << " CG steps on the top left block "; + else if (name == "rhs") + std::cout << " " << solver_control.last_step() + << " CG steps for computing the r.h.s. "; + else + abort (); + + std::cout << solver_control.initial_value() << "->" << solver_control.last_value() + << std::endl; #endif -} + } -template -class BlockSchurPreconditioner : public Subscriptor -{ - public: - BlockSchurPreconditioner (const BlockSparseMatrix &S, - const InverseMatrix,PreconditionerMp> &Mpinv, - const PreconditionerA &Apreconditioner); - - void vmult (BlockVector &dst, - const BlockVector &src) const; - - private: - const SmartPointer > system_matrix; - const SmartPointer, - PreconditionerMp > > m_inverse; - const PreconditionerA &a_preconditioner; - - mutable Vector tmp; - -}; - -template -BlockSchurPreconditioner::BlockSchurPreconditioner( - const BlockSparseMatrix &S, - const InverseMatrix,PreconditionerMp> &Mpinv, - const PreconditionerA &Apreconditioner - ) - : - system_matrix (&S), - m_inverse (&Mpinv), - a_preconditioner (Apreconditioner), - tmp (S.block(1,1).m()) -{} - - // Now the interesting function, the multiplication of - // the preconditioner with a BlockVector. -template -void BlockSchurPreconditioner::vmult ( - BlockVector &dst, - const BlockVector &src) const -{ - // Form u_new = A^{-1} u - a_preconditioner.vmult (dst.block(0), src.block(0)); - // Form tmp = - B u_new + p - // (SparseMatrix::residual - // does precisely this) - system_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1)); - // Change sign in tmp - tmp *= -1; - // Multiply by approximate Schur complement - // (i.e. a pressure mass matrix) - m_inverse->vmult (dst.block(1), tmp); -} + template + class BlockSchurPreconditioner : public Subscriptor + { + public: + BlockSchurPreconditioner (const BlockSparseMatrix &S, + const InverseMatrix,PreconditionerMp> &Mpinv, + const PreconditionerA &Apreconditioner); + + void vmult (BlockVector &dst, + const BlockVector &src) const; + + private: + const SmartPointer > system_matrix; + const SmartPointer, + PreconditionerMp > > m_inverse; + const PreconditionerA &a_preconditioner; + + mutable Vector tmp; + + }; + + template + BlockSchurPreconditioner::BlockSchurPreconditioner( + const BlockSparseMatrix &S, + const InverseMatrix,PreconditionerMp> &Mpinv, + const PreconditionerA &Apreconditioner + ) + : + system_matrix (&S), + m_inverse (&Mpinv), + a_preconditioner (Apreconditioner), + tmp (S.block(1,1).m()) + {} + + // Now the interesting function, the multiplication of + // the preconditioner with a BlockVector. + template + void BlockSchurPreconditioner::vmult ( + BlockVector &dst, + const BlockVector &src) const + { + // Form u_new = A^{-1} u + a_preconditioner.vmult (dst.block(0), src.block(0)); + // Form tmp = - B u_new + p + // (SparseMatrix::residual + // does precisely this) + system_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1)); + // Change sign in tmp + tmp *= -1; + // Multiply by approximate Schur complement + // (i.e. a pressure mass matrix) + m_inverse->vmult (dst.block(1), tmp); + } -template -class SchurComplement : public Subscriptor -{ - public: - SchurComplement (const BlockSparseMatrix &system_matrix, - const InverseMatrix, Preconditioner> &A_inverse); + template + class SchurComplement : public Subscriptor + { + public: + SchurComplement (const BlockSparseMatrix &system_matrix, + const InverseMatrix, Preconditioner> &A_inverse); - void vmult (Vector &dst, - const Vector &src) const; + void vmult (Vector &dst, + const Vector &src) const; - unsigned int m() const - { - return system_matrix->block(1,1).m(); - } + unsigned int m() const + { + return system_matrix->block(1,1).m(); + } - unsigned int n() const - { - return system_matrix->block(1,1).n(); - } + unsigned int n() const + { + return system_matrix->block(1,1).n(); + } - private: - const SmartPointer > system_matrix; - const SmartPointer, Preconditioner> > A_inverse; + private: + const SmartPointer > system_matrix; + const SmartPointer, Preconditioner> > A_inverse; - mutable Vector tmp1, tmp2; -}; + mutable Vector tmp1, tmp2; + }; -template -SchurComplement:: -SchurComplement (const BlockSparseMatrix &system_matrix, - const InverseMatrix,Preconditioner> &A_inverse) - : - system_matrix (&system_matrix), - A_inverse (&A_inverse), - tmp1 (system_matrix.block(0,0).m()), - tmp2 (system_matrix.block(0,0).m()) -{} + template + SchurComplement:: + SchurComplement (const BlockSparseMatrix &system_matrix, + const InverseMatrix,Preconditioner> &A_inverse) + : + system_matrix (&system_matrix), + A_inverse (&A_inverse), + tmp1 (system_matrix.block(0,0).m()), + tmp2 (system_matrix.block(0,0).m()) + {} -template -void SchurComplement::vmult (Vector &dst, - const Vector &src) const -{ - system_matrix->block(0,1).vmult (tmp1, src); - A_inverse->name = "in schur"; - A_inverse->vmult (tmp2, tmp1); - system_matrix->block(1,0).vmult (dst, tmp2); - dst *= -1; - system_matrix->block(1,1).vmult_add (dst, src); - dst *= -1; -} + template + void SchurComplement::vmult (Vector &dst, + const Vector &src) const + { + system_matrix->block(0,1).vmult (tmp1, src); + A_inverse->name = "in schur"; + A_inverse->vmult (tmp2, tmp1); + system_matrix->block(1,0).vmult (dst, tmp2); + dst *= -1; + system_matrix->block(1,1).vmult_add (dst, src); + dst *= -1; + } -template -StokesProblem::StokesProblem (const unsigned int degree) - : - degree (degree), - triangulation (Triangulation::limit_level_difference_at_vertices), - fe (FE_Q(degree+1), dim, - FE_Q(degree), 1), - dof_handler (triangulation) -{} + template + StokesProblem::StokesProblem (const unsigned int degree) + : + degree (degree), + triangulation (Triangulation::limit_level_difference_at_vertices), + fe (FE_Q(degree+1), dim, + FE_Q(degree), 1), + dof_handler (triangulation) + {} -template -void StokesProblem::setup_dofs () -{ - A_preconditioner.reset (); - mg_A_preconditioner.resize (0); - system_matrix.clear (); + template + void StokesProblem::setup_dofs () + { + A_preconditioner.reset (); + mg_A_preconditioner.resize (0); + system_matrix.clear (); - dof_handler.distribute_dofs (fe); + dof_handler.distribute_dofs (fe); // DoFRenumbering::Cuthill_McKee (dof_handler); - std::vector block_component (dim+1,0); - block_component[dim] = 1; - DoFRenumbering::component_wise (dof_handler, block_component); + std::vector block_component (dim+1,0); + block_component[dim] = 1; + DoFRenumbering::component_wise (dof_handler, block_component); - { - constraints.clear (); - typename FunctionMap::type dirichlet_boundary; - ZeroFunction homogeneous_dirichlet_bc (dim+1); //TODO: go back to BoundaryValues - - dirichlet_boundary[0] = &homogeneous_dirichlet_bc; - MappingQ1 mapping; - - std::vector component_mask (dim+1, true); - component_mask[dim] = false; - VectorTools::interpolate_boundary_values (mapping, - dof_handler, - dirichlet_boundary, - constraints, - component_mask); - - DoFTools::make_hanging_node_constraints (dof_handler, - constraints); - - mg_constrained_dofs.clear(); - mg_constrained_dofs.initialize(dof_handler, dirichlet_boundary); - } + { + constraints.clear (); + typename FunctionMap::type dirichlet_boundary; + ZeroFunction homogeneous_dirichlet_bc (dim+1); //TODO: go back to BoundaryValues + + dirichlet_boundary[0] = &homogeneous_dirichlet_bc; + MappingQ1 mapping; + + std::vector component_mask (dim+1, true); + component_mask[dim] = false; + VectorTools::interpolate_boundary_values (mapping, + dof_handler, + dirichlet_boundary, + constraints, + component_mask); + + DoFTools::make_hanging_node_constraints (dof_handler, + constraints); + + mg_constrained_dofs.clear(); + mg_constrained_dofs.initialize(dof_handler, dirichlet_boundary); + } - constraints.close (); + constraints.close (); - std::vector dofs_per_block (2); - DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, - block_component); - const unsigned int n_u = dofs_per_block[0], - n_p = dofs_per_block[1]; + std::vector dofs_per_block (2); + DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, + block_component); + const unsigned int n_u = dofs_per_block[0], + n_p = dofs_per_block[1]; - std::cout << " Number of active cells: " - << triangulation.n_active_cells() - << std::endl - << " Number of degrees of freedom: " - << dof_handler.n_dofs() - << " (" << n_u << '+' << n_p << ')' - << std::endl; + std::cout << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl + << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << " (" << n_u << '+' << n_p << ')' + << std::endl; - { - BlockCompressedSimpleSparsityPattern csp (2,2); + { + BlockCompressedSimpleSparsityPattern csp (2,2); - csp.block(0,0).reinit (n_u, n_u); - csp.block(1,0).reinit (n_p, n_u); - csp.block(0,1).reinit (n_u, n_p); - csp.block(1,1).reinit (n_p, n_p); + csp.block(0,0).reinit (n_u, n_u); + csp.block(1,0).reinit (n_p, n_u); + csp.block(0,1).reinit (n_u, n_p); + csp.block(1,1).reinit (n_p, n_p); - csp.collect_sizes(); + csp.collect_sizes(); - DoFTools::make_sparsity_pattern ( + DoFTools::make_sparsity_pattern ( static_cast&>(dof_handler), csp, constraints, false); - sparsity_pattern.copy_from (csp); - } + sparsity_pattern.copy_from (csp); + } - system_matrix.reinit (sparsity_pattern); + system_matrix.reinit (sparsity_pattern); - solution.reinit (2); - solution.block(0).reinit (n_u); - solution.block(1).reinit (n_p); - solution.collect_sizes (); + solution.reinit (2); + solution.block(0).reinit (n_u); + solution.block(1).reinit (n_p); + solution.collect_sizes (); - system_rhs.reinit (2); - system_rhs.block(0).reinit (n_u); - system_rhs.block(1).reinit (n_p); - system_rhs.collect_sizes (); + system_rhs.reinit (2); + system_rhs.block(0).reinit (n_u); + system_rhs.block(1).reinit (n_p); + system_rhs.collect_sizes (); - //now setup stuff for mg - const unsigned int nlevels = triangulation.n_levels(); + //now setup stuff for mg + const unsigned int nlevels = triangulation.n_levels(); - mg_matrices.resize(0, nlevels-1); - mg_matrices.clear (); - mg_interface_matrices.resize(0, nlevels-1); - mg_interface_matrices.clear (); - mg_sparsity.resize(0, nlevels-1); + mg_matrices.resize(0, nlevels-1); + mg_matrices.clear (); + mg_interface_matrices.resize(0, nlevels-1); + mg_interface_matrices.clear (); + mg_sparsity.resize(0, nlevels-1); - mg_dofs_per_component.resize (nlevels); - for (unsigned int level=0; level -void StokesProblem::assemble_system () -{ - system_matrix=0; - system_rhs=0; - - QGauss quadrature_formula(degree+2); - - FEValues fe_values (fe, quadrature_formula, - update_values | - update_quadrature_points | - update_JxW_values | - update_gradients); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); - - FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); - Vector local_rhs (dofs_per_cell); - - std::vector local_dof_indices (dofs_per_cell); + template + void StokesProblem::assemble_system () + { + system_matrix=0; + system_rhs=0; - const RightHandSide right_hand_side; - std::vector > rhs_values (n_q_points, - Vector(dim+1)); + QGauss quadrature_formula(degree+2); + FEValues fe_values (fe, quadrature_formula, + update_values | + update_quadrature_points | + update_JxW_values | + update_gradients); - const FEValuesExtractors::Vector velocities (0); - const FEValuesExtractors::Scalar pressure (dim); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); + Vector local_rhs (dofs_per_cell); - std::vector > phi_grads_u (dofs_per_cell); - std::vector div_phi_u (dofs_per_cell); - std::vector phi_p (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); - typename MGDoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - { - fe_values.reinit (cell); - local_matrix = 0; - local_rhs = 0; + const RightHandSide right_hand_side; + std::vector > rhs_values (n_q_points, + Vector(dim+1)); - right_hand_side.vector_value_list(fe_values.get_quadrature_points(), - rhs_values); - for (unsigned int q=0; q > phi_grads_u (dofs_per_cell); + std::vector div_phi_u (dofs_per_cell); + std::vector phi_p (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); - } -} + typename MGDoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + local_matrix = 0; + local_rhs = 0; + + right_hand_side.vector_value_list(fe_values.get_quadrature_points(), + rhs_values); + + for (unsigned int q=0; qget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global (local_matrix, local_rhs, + local_dof_indices, + system_matrix, system_rhs); + } + } -template -void StokesProblem::assemble_multigrid () -{ - QGauss quadrature_formula(degree+2); - FEValues fe_values (fe, quadrature_formula, - update_values | - update_quadrature_points | - update_JxW_values | - update_gradients); + template + void StokesProblem::assemble_multigrid () + { + QGauss quadrature_formula(degree+2); + FEValues fe_values (fe, quadrature_formula, + update_values | + update_quadrature_points | + update_JxW_values | + update_gradients); - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); - FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); - std::vector local_dof_indices (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); - const FEValuesExtractors::Vector velocities (0); - const FEValuesExtractors::Scalar pressure (dim); + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); - std::vector > phi_grads_u (dofs_per_cell); - std::vector div_phi_u (dofs_per_cell); - std::vector phi_p (dofs_per_cell); + std::vector > phi_grads_u (dofs_per_cell); + std::vector div_phi_u (dofs_per_cell); + std::vector phi_p (dofs_per_cell); - std::vector > interface_dofs - = mg_constrained_dofs.get_refinement_edge_indices (); - std::vector > boundary_interface_dofs - = mg_constrained_dofs.get_refinement_edge_boundary_indices (); + std::vector > interface_dofs + = mg_constrained_dofs.get_refinement_edge_indices (); + std::vector > boundary_interface_dofs + = mg_constrained_dofs.get_refinement_edge_boundary_indices (); - std::vector boundary_constraints (triangulation.n_levels()); - std::vector boundary_interface_constraints (triangulation.n_levels()); - for (unsigned int level=0; level boundary_constraints (triangulation.n_levels()); + std::vector boundary_interface_constraints (triangulation.n_levels()); + for (unsigned int level=0; level::cell_iterator - cell = dof_handler.begin(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - { - // Remember the level of the - // current cell. - const unsigned int level = cell->level(); - // Compute the values specified - // by update flags above. - fe_values.reinit (cell); - local_matrix = 0; - - for (unsigned int q=0; q::cell_iterator + cell = dof_handler.begin(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) { - for (unsigned int k=0; klevel(); + // Compute the values specified + // by update flags above. + fe_values.reinit (cell); + local_matrix = 0; + + for (unsigned int q=0; qget_mg_dof_indices (local_dof_indices); + boundary_constraints[level] + .distribute_local_to_global (local_matrix, + local_dof_indices, + mg_matrices[level]); + + for (unsigned int i=0; iget_mg_dof_indices (local_dof_indices); - boundary_constraints[level] - .distribute_local_to_global (local_matrix, - local_dof_indices, - mg_matrices[level]); - - for (unsigned int i=0; i::type>(new typename InnerPreconditioner::type()); - mg_A_preconditioner[level] - ->initialize (mg_matrices[level].block(0,0), - typename InnerPreconditioner::type::AdditionalData()); - } -} + mg_A_preconditioner.resize (triangulation.n_levels()); + for (unsigned int level=0; level::type>(new typename InnerPreconditioner::type()); + mg_A_preconditioner[level] + ->initialize (mg_matrices[level].block(0,0), + typename InnerPreconditioner::type::AdditionalData()); + } + } -template -class SchurComplementSmoother -{ - public: - struct AdditionalData - { - const InnerPreconditioner *A_preconditioner; - }; + template + class SchurComplementSmoother + { + public: + struct AdditionalData + { + const InnerPreconditioner *A_preconditioner; + }; - void initialize (const BlockSparseMatrix &system_matrix, - const AdditionalData &data); + void initialize (const BlockSparseMatrix &system_matrix, + const AdditionalData &data); - void vmult (BlockVector &dst, - const BlockVector &src) const; + void vmult (BlockVector &dst, + const BlockVector &src) const; - void Tvmult (BlockVector &dst, - const BlockVector &src) const; + void Tvmult (BlockVector &dst, + const BlockVector &src) const; - void clear (); + void clear (); - private: - SmartPointer > system_matrix; - SmartPointer A_preconditioner; -}; + private: + SmartPointer > system_matrix; + SmartPointer A_preconditioner; + }; -template -void -SchurComplementSmoother:: -initialize (const BlockSparseMatrix &system_matrix, - const AdditionalData &data) -{ - this->system_matrix = &system_matrix; - this->A_preconditioner = data.A_preconditioner; -} + template + void + SchurComplementSmoother:: + initialize (const BlockSparseMatrix &system_matrix, + const AdditionalData &data) + { + this->system_matrix = &system_matrix; + this->A_preconditioner = data.A_preconditioner; + } -template -void -SchurComplementSmoother:: -vmult (BlockVector &dst, - const BlockVector &src) const -{ + template + void + SchurComplementSmoother:: + vmult (BlockVector &dst, + const BlockVector &src) const + { #ifdef STEP_42_TEST - std::cout << "Entering smoother with " << dst.size() << " unknowns" << std::endl; + std::cout << "Entering smoother with " << dst.size() << " unknowns" << std::endl; #endif - SparseDirectUMFPACK direct_solver; - direct_solver.initialize(*system_matrix); - Vector solution, rhs; - solution = dst; - rhs = src; - direct_solver.vmult(solution, rhs); - dst = solution; + SparseDirectUMFPACK direct_solver; + direct_solver.initialize(*system_matrix); + Vector solution, rhs; + solution = dst; + rhs = src; + direct_solver.vmult(solution, rhs); + dst = solution; /* const InverseMatrix,InnerPreconditioner> - A_inverse (system_matrix->block(0,0), *A_preconditioner); + A_inverse (system_matrix->block(0,0), *A_preconditioner); Vector tmp (dst.block(0).size()); { - Vector schur_rhs (dst.block(1).size()); - A_inverse.name = "rhs"; - A_inverse.vmult (tmp, src.block(0)); + Vector schur_rhs (dst.block(1).size()); + A_inverse.name = "rhs"; + A_inverse.vmult (tmp, src.block(0)); // std::cout << " TMP " << tmp.l2_norm() << std::endl; - system_matrix->block(1,0).vmult (schur_rhs, tmp); - schur_rhs -= src.block(1); +system_matrix->block(1,0).vmult (schur_rhs, tmp); +schur_rhs -= src.block(1); // std::cout << " BLOCK 1 " << src.block(1).l2_norm() << std::endl; // std::cout << " SCHUR RHS " << schur_rhs.l2_norm() << std::endl; - SchurComplement - schur_complement (*system_matrix, A_inverse); - - // The usual control structures for - // the solver call are created... - SolverControl solver_control (dst.block(1).size(), - 1e-1*schur_rhs.l2_norm()); - SolverGMRES<> cg (solver_control); - -#ifdef STEP_42_TEST - std::cout << " Starting Schur complement solver -- " - << schur_complement.m() << " unknowns" - << std::endl; -#endif - try - { - cg.solve (schur_complement, dst.block(1), schur_rhs, - PreconditionIdentity()); - } - catch (...) - { - std::cout << "Failure in " << __PRETTY_FUNCTION__ << std::endl; - std::cout << schur_rhs.l2_norm () << std::endl; - abort (); - } +SchurComplement +schur_complement (*system_matrix, A_inverse); + + // The usual control structures for + // the solver call are created... + SolverControl solver_control (dst.block(1).size(), + 1e-1*schur_rhs.l2_norm()); + SolverGMRES<> cg (solver_control); + + #ifdef STEP_42_TEST + std::cout << " Starting Schur complement solver -- " + << schur_complement.m() << " unknowns" + << std::endl; + #endif + try + { + cg.solve (schur_complement, dst.block(1), schur_rhs, + PreconditionIdentity()); + } + catch (...) + { + std::cout << "Failure in " << __PRETTY_FUNCTION__ << std::endl; + std::cout << schur_rhs.l2_norm () << std::endl; + abort (); + } // no constraints to be taken care of here #ifdef STEP_42_TEST - std::cout << " " - << solver_control.last_step() - << " CG Schur complement iterations in smoother " - << solver_control.initial_value() << "->" << solver_control.last_value() - << std::endl; +std::cout << " " +<< solver_control.last_step() +<< " CG Schur complement iterations in smoother " +<< solver_control.initial_value() << "->" << solver_control.last_value() +<< std::endl; #endif - } +} - { - system_matrix->block(0,1).vmult (tmp, dst.block(1)); - tmp *= -1; - tmp += src.block(0); +{ +system_matrix->block(0,1).vmult (tmp, dst.block(1)); +tmp *= -1; +tmp += src.block(0); - A_inverse.name = "top left"; - A_inverse.vmult (dst.block(0), tmp); +A_inverse.name = "top left"; +A_inverse.vmult (dst.block(0), tmp); // no constraints here either - } +} #ifdef STEP_42_TEST - std::cout << "Exiting smoother with " << dst.size() << " unknowns" << std::endl; +std::cout << "Exiting smoother with " << dst.size() << " unknowns" << std::endl; #endif */ -} + } -template -void -SchurComplementSmoother::clear () -{} + template + void + SchurComplementSmoother::clear () + {} -template -void -SchurComplementSmoother:: -Tvmult (BlockVector &, - const BlockVector &) const -{ - Assert (false, ExcNotImplemented()); -} + template + void + SchurComplementSmoother:: + Tvmult (BlockVector &, + const BlockVector &) const + { + Assert (false, ExcNotImplemented()); + } -template -void StokesProblem::solve () -{ - system_matrix.block(1,1) = 0; - assemble_multigrid (); - typedef PreconditionMG, MGTransferPrebuilt > > - MGPREC; - - GrowingVectorMemory > mg_vector_memory; - - MGTransferPrebuilt > mg_transfer(constraints, mg_constrained_dofs); - std::vector block_component (dim+1,0); - block_component[dim] = 1; - mg_transfer.set_component_to_block_map (block_component); - mg_transfer.build_matrices(dof_handler); - - FullMatrix mg_coarse_matrix; - mg_coarse_matrix.copy_from (mg_matrices[0]); - MGCoarseGridHouseholder > mg_coarse; - mg_coarse.initialize(mg_coarse_matrix); - - MGMatrix, BlockVector > - mg_matrix(&mg_matrices); - MGMatrix, BlockVector > - mg_interface_up(&mg_interface_matrices); - MGMatrix, BlockVector > - mg_interface_down(&mg_interface_matrices); - - typedef - SchurComplementSmoother::type> - Smoother; - - MGSmootherPrecondition, - Smoother, - BlockVector > - mg_smoother(mg_vector_memory); - - MGLevelObject - smoother_data (0, triangulation.n_levels()-1); - - for (unsigned int level=0; level > mg(dof_handler, - mg_matrix, - mg_coarse, - mg_transfer, - mg_smoother, - mg_smoother); - mg.set_debug(3); - mg.set_edge_matrices(mg_interface_down, mg_interface_up); - - MGPREC preconditioner(dof_handler, mg, mg_transfer); - - SolverControl solver_control (system_matrix.m(), - 1e-6*system_rhs.l2_norm()); - GrowingVectorMemory > vector_memory; - SolverGMRES >::AdditionalData gmres_data; - gmres_data.max_n_tmp_vectors = 100; - - SolverGMRES > gmres(solver_control, vector_memory, - gmres_data); + template + void StokesProblem::solve () + { + system_matrix.block(1,1) = 0; + assemble_multigrid (); + typedef PreconditionMG, MGTransferPrebuilt > > + MGPREC; + + GrowingVectorMemory > mg_vector_memory; + + MGTransferPrebuilt > mg_transfer(constraints, mg_constrained_dofs); + std::vector block_component (dim+1,0); + block_component[dim] = 1; + mg_transfer.set_component_to_block_map (block_component); + mg_transfer.build_matrices(dof_handler); + + FullMatrix mg_coarse_matrix; + mg_coarse_matrix.copy_from (mg_matrices[0]); + MGCoarseGridHouseholder > mg_coarse; + mg_coarse.initialize(mg_coarse_matrix); + + MGMatrix, BlockVector > + mg_matrix(&mg_matrices); + MGMatrix, BlockVector > + mg_interface_up(&mg_interface_matrices); + MGMatrix, BlockVector > + mg_interface_down(&mg_interface_matrices); + + typedef + SchurComplementSmoother::type> + Smoother; + + MGSmootherPrecondition, + Smoother, + BlockVector > + mg_smoother(mg_vector_memory); + + MGLevelObject + smoother_data (0, triangulation.n_levels()-1); + + for (unsigned int level=0; level > mg(dof_handler, + mg_matrix, + mg_coarse, + mg_transfer, + mg_smoother, + mg_smoother); + mg.set_debug(3); + mg.set_edge_matrices(mg_interface_down, mg_interface_up); + + MGPREC preconditioner(dof_handler, mg, mg_transfer); + + SolverControl solver_control (system_matrix.m(), + 1e-6*system_rhs.l2_norm()); + GrowingVectorMemory > vector_memory; + SolverGMRES >::AdditionalData gmres_data; + gmres_data.max_n_tmp_vectors = 100; + + SolverGMRES > gmres(solver_control, vector_memory, + gmres_data); // PreconditionIdentity precondition_identity; #ifdef STEP_42_TEST - std::cout << "Starting outer GMRES complement solver" << std::endl; + std::cout << "Starting outer GMRES complement solver" << std::endl; #endif - try - { - gmres.solve(system_matrix, solution, system_rhs, - preconditioner); - } - catch (...) - { - std::cout << "Failure in " << __PRETTY_FUNCTION__ << std::endl; - abort (); - } + try + { + gmres.solve(system_matrix, solution, system_rhs, + preconditioner); + } + catch (...) + { + std::cout << "Failure in " << __PRETTY_FUNCTION__ << std::endl; + abort (); + } - constraints.distribute (solution); + constraints.distribute (solution); - std::cout << solver_control.last_step() - << " outer GMRES iterations "; -} + std::cout << solver_control.last_step() + << " outer GMRES iterations "; + } template -void StokesProblem::solve_block () -{ - std::cout << " Computing preconditioner..." << std::endl << std::flush; - - A_preconditioner - = std_cxx1x::shared_ptr::type>(new typename InnerPreconditioner::type()); - A_preconditioner->initialize (system_matrix.block(0,0), - typename InnerPreconditioner::type::AdditionalData()); - - SparseMatrix pressure_mass_matrix; - pressure_mass_matrix.reinit(sparsity_pattern.block(1,1)); - pressure_mass_matrix.copy_from(system_matrix.block(1,1)); - system_matrix.block(1,1) = 0; - - SparseILU pmass_preconditioner; - pmass_preconditioner.initialize (pressure_mass_matrix, - SparseILU::AdditionalData()); - - InverseMatrix,SparseILU > - m_inverse (pressure_mass_matrix, pmass_preconditioner); - - BlockSchurPreconditioner::type, - SparseILU > + void StokesProblem::solve_block () + { + std::cout << " Computing preconditioner..." << std::endl << std::flush; + + A_preconditioner + = std_cxx1x::shared_ptr::type>(new typename InnerPreconditioner::type()); + A_preconditioner->initialize (system_matrix.block(0,0), + typename InnerPreconditioner::type::AdditionalData()); + + SparseMatrix pressure_mass_matrix; + pressure_mass_matrix.reinit(sparsity_pattern.block(1,1)); + pressure_mass_matrix.copy_from(system_matrix.block(1,1)); + system_matrix.block(1,1) = 0; + + SparseILU pmass_preconditioner; + pmass_preconditioner.initialize (pressure_mass_matrix, + SparseILU::AdditionalData()); + + InverseMatrix,SparseILU > + m_inverse (pressure_mass_matrix, pmass_preconditioner); + + BlockSchurPreconditioner::type, + SparseILU > preconditioner (system_matrix, m_inverse, *A_preconditioner); - SolverControl solver_control (system_matrix.m(), - 1e-6*system_rhs.l2_norm()); - GrowingVectorMemory > vector_memory; - SolverGMRES >::AdditionalData gmres_data; - gmres_data.max_n_tmp_vectors = 100; + SolverControl solver_control (system_matrix.m(), + 1e-6*system_rhs.l2_norm()); + GrowingVectorMemory > vector_memory; + SolverGMRES >::AdditionalData gmres_data; + gmres_data.max_n_tmp_vectors = 100; - SolverGMRES > gmres(solver_control, vector_memory, - gmres_data); + SolverGMRES > gmres(solver_control, vector_memory, + gmres_data); - gmres.solve(system_matrix, solution, system_rhs, - preconditioner); + gmres.solve(system_matrix, solution, system_rhs, + preconditioner); - constraints.distribute (solution); + constraints.distribute (solution); - std::cout << " " - << solver_control.last_step() - << " block GMRES iterations "; -} + std::cout << " " + << solver_control.last_step() + << " block GMRES iterations "; + } -template -void -StokesProblem::output_results (const unsigned int refinement_cycle) const -{ - std::vector solution_names (dim, "velocity"); - solution_names.push_back ("pressure"); + template + void + StokesProblem::output_results (const unsigned int refinement_cycle) const + { + std::vector solution_names (dim, "velocity"); + solution_names.push_back ("pressure"); - std::vector + std::vector + data_component_interpretation + (dim, DataComponentInterpretation::component_is_part_of_vector); data_component_interpretation - (dim, DataComponentInterpretation::component_is_part_of_vector); - data_component_interpretation - .push_back (DataComponentInterpretation::component_is_scalar); - - 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 (); - - 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); -} + .push_back (DataComponentInterpretation::component_is_scalar); + + 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 (); + + 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); + } -template -void -StokesProblem::refine_mesh () -{ - Vector estimated_error_per_cell (triangulation.n_active_cells()); - - std::vector component_mask (dim+1, false); - component_mask[dim] = true; - KellyErrorEstimator::estimate (static_cast&>(dof_handler), - QGauss(degree+1), - typename FunctionMap::type(), - solution, - estimated_error_per_cell, - component_mask); - - GridRefinement::refine_and_coarsen_fixed_number (triangulation, - estimated_error_per_cell, - 0.3, 0.0); - triangulation.execute_coarsening_and_refinement (); -} + template + void + StokesProblem::refine_mesh () + { + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + std::vector component_mask (dim+1, false); + component_mask[dim] = true; + KellyErrorEstimator::estimate (static_cast&>(dof_handler), + QGauss(degree+1), + typename FunctionMap::type(), + solution, + estimated_error_per_cell, + component_mask); + + GridRefinement::refine_and_coarsen_fixed_number (triangulation, + estimated_error_per_cell, + 0.3, 0.0); + triangulation.execute_coarsening_and_refinement (); + } -template -void StokesProblem::run () -{ + template + void StokesProblem::run () { - std::vector subdivisions (dim, 1); - subdivisions[0] = 1; - - const Point bottom_left = (dim == 2 ? - Point(0,0) : - Point(0,0,0)); - const Point top_right = (dim == 2 ? - Point(1,1) : - Point(1,1,1)); - - GridGenerator::subdivided_hyper_rectangle (triangulation, - subdivisions, - bottom_left, - top_right); - } + { + std::vector subdivisions (dim, 1); + subdivisions[0] = 1; + + const Point bottom_left = (dim == 2 ? + Point(0,0) : + Point(0,0,0)); + const Point top_right = (dim == 2 ? + Point(1,1) : + Point(1,1,1)); + + GridGenerator::subdivided_hyper_rectangle (triangulation, + subdivisions, + bottom_left, + top_right); + } - 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)->center()[0] == 1) - cell->face(f)->set_all_boundary_indicators(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)->center()[0] == 1) + cell->face(f)->set_all_boundary_indicators(1); - triangulation.refine_global (1); + triangulation.refine_global (1); - for (unsigned int refinement_cycle = 0; refinement_cycle<10; - ++refinement_cycle) - { - std::cout << "Refinement cycle " << refinement_cycle << std::endl; + for (unsigned int refinement_cycle = 0; refinement_cycle<10; + ++refinement_cycle) + { + std::cout << "Refinement cycle " << refinement_cycle << std::endl; - if (refinement_cycle > 0) - refine_mesh (); + if (refinement_cycle > 0) + refine_mesh (); - std::ostringstream out_filename; - out_filename << "gitter" - << refinement_cycle - << ".eps"; + std::ostringstream out_filename; + out_filename << "gitter" + << refinement_cycle + << ".eps"; - std::ofstream grid_output (out_filename.str().c_str()); - GridOut grid_out; - grid_out.write_eps (triangulation, grid_output); + std::ofstream grid_output (out_filename.str().c_str()); + GridOut grid_out; + grid_out.write_eps (triangulation, grid_output); - setup_dofs (); + setup_dofs (); - std::cout << " Assembling..." << std::endl << std::flush; - assemble_system (); + std::cout << " Assembling..." << std::endl << std::flush; + assemble_system (); - std::cout << " Solving..." << std::flush; + std::cout << " Solving..." << std::flush; - solve_block (); - output_results (refinement_cycle); - system ("mv solution-* block"); + solve_block (); + output_results (refinement_cycle); + system ("mv solution-* block"); - solution = 0; + solution = 0; - solve (); - output_results (refinement_cycle); - system ("mv solution-* mg"); + solve (); + output_results (refinement_cycle); + system ("mv solution-* mg"); - std::cout << std::endl; - } + std::cout << std::endl; + } + } } @@ -1149,6 +1152,9 @@ int main () { try { + using namespace dealii; + using namespace Step42; + deallog.depth_console (0); StokesProblem<2> flow_problem(1); diff --git a/deal.II/examples/step-43/step-43.cc b/deal.II/examples/step-43/step-43.cc index b3f5d3a521..feb9f2a8cf 100644 --- a/deal.II/examples/step-43/step-43.cc +++ b/deal.II/examples/step-43/step-43.cc @@ -4,7 +4,7 @@ /* $Id$ */ /* */ -/* Copyright (C) 2010 by Chih-Che Chueh and the deal.II authors */ +/* Copyright (C) 2010, 2011 by Chih-Che Chueh and the deal.II authors */ /* */ /* This file is subject to QPL and may not be distributed */ /* without copyright and license information. Please refer */ @@ -18,14 +18,14 @@ // the functionality of these well-known // deal.II library files and some C++ header // files. - // + // // In this program, we use a tensor-valued // coefficient. Since it may have a spatial // dependence, we consider it a tensor-valued // function. The following include file // provides the TensorFunction class that // offers such functionality: - // + // // Then we need to include some header files // that provide vector, matrix, and // preconditioner classes that implement @@ -35,7 +35,7 @@ // interfaces to the matrix and vector // classes based on Trilinos as well as // Trilinos preconditioners: - // + // // At the end of this top-matter, we import // all deal.II names into the global // namespace: @@ -79,2434 +79,2441 @@ #include #include -using namespace dealii; +namespace Step43 +{ + using namespace dealii; - // @sect3{The InverseMatrix class template} + // @sect3{The InverseMatrix class template} - // This part is exactly the same as that used in step-31. + // This part is exactly the same as that used in step-31. - // @sect3{Schur complement preconditioner} + // @sect3{Schur complement preconditioner} - // This part for the Schur complement - // preconditioner is almost the same as that - // used in step-31. The only difference is - // that the original variable name - // stokes_matrix is replaced by another name - // darcy_matrix to satisfy our problem. -namespace LinearSolvers -{ - template - class InverseMatrix : public Subscriptor + // This part for the Schur complement + // preconditioner is almost the same as that + // used in step-31. The only difference is + // that the original variable name + // stokes_matrix is replaced by another name + // darcy_matrix to satisfy our problem. + namespace LinearSolvers { - public: - InverseMatrix (const Matrix &m, - const Preconditioner &preconditioner); + template + class InverseMatrix : public Subscriptor + { + public: + InverseMatrix (const Matrix &m, + const Preconditioner &preconditioner); - template - void vmult (VectorType &dst, - const VectorType &src) const; + template + void vmult (VectorType &dst, + const VectorType &src) const; - private: - const SmartPointer matrix; - const Preconditioner &preconditioner; - }; + private: + const SmartPointer matrix; + const Preconditioner &preconditioner; + }; - template - InverseMatrix:: - InverseMatrix (const Matrix &m, - const Preconditioner &preconditioner) - : - matrix (&m), - preconditioner (preconditioner) - {} + template + InverseMatrix:: + InverseMatrix (const Matrix &m, + const Preconditioner &preconditioner) + : + matrix (&m), + preconditioner (preconditioner) + {} - template - template - void - InverseMatrix:: - vmult (VectorType &dst, - const VectorType &src) const - { - SolverControl solver_control (src.size(), 1e-7*src.l2_norm()); - SolverCG cg (solver_control); + template + template + void + InverseMatrix:: + vmult (VectorType &dst, + const VectorType &src) const + { + SolverControl solver_control (src.size(), 1e-7*src.l2_norm()); + SolverCG cg (solver_control); - dst = 0; + dst = 0; - try - { - cg.solve (*matrix, dst, src, preconditioner); - } - catch (std::exception &e) - { - Assert (false, ExcMessage(e.what())); - } + try + { + cg.solve (*matrix, dst, src, preconditioner); + } + catch (std::exception &e) + { + Assert (false, ExcMessage(e.what())); + } + } + + template + class BlockSchurPreconditioner : public Subscriptor + { + public: + BlockSchurPreconditioner ( + const TrilinosWrappers::BlockSparseMatrix &S, + const InverseMatrix &Mpinv, + const PreconditionerA &Apreconditioner); + + void vmult (TrilinosWrappers::BlockVector &dst, + const TrilinosWrappers::BlockVector &src) const; + + private: + const SmartPointer darcy_matrix; + const SmartPointer > m_inverse; + const PreconditionerA &a_preconditioner; + + mutable TrilinosWrappers::Vector tmp; + }; + + + + template + BlockSchurPreconditioner:: + BlockSchurPreconditioner(const TrilinosWrappers::BlockSparseMatrix &S, + const InverseMatrix &Mpinv, + const PreconditionerA &Apreconditioner) + : + darcy_matrix (&S), + m_inverse (&Mpinv), + a_preconditioner (Apreconditioner), + tmp (darcy_matrix->block(1,1).m()) + {} + + + template + void BlockSchurPreconditioner::vmult ( + TrilinosWrappers::BlockVector &dst, + const TrilinosWrappers::BlockVector &src) const + { + a_preconditioner.vmult (dst.block(0), src.block(0)); + darcy_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1)); + tmp *= -1; + m_inverse->vmult (dst.block(1), tmp); + } } - template - class BlockSchurPreconditioner : public Subscriptor + + // @sect3{The TwoPhaseFlowProblem class} + + // The definition of the class that defines + // the top-level logic of solving the + // time-dependent advection-dominated + // two-phase flow problem (or + // Buckley-Leverett problem + // [Buckley 1942]) is mainly based on + // three tutorial programs (step-21, step-31, + // step-33). The main difference is that, + // since adaptive operator splitting is + // considered, we need a bool-type variable + // solve_pressure_velocity_part to tell us + // when we need to solve the pressure and + // velocity part, need another bool-type + // variable + // previous_solve_pressure_velocity_part to + // determine if we have to cumulate + // micro-time steps that we need them to do + // extrapolation for the total velocity, and + // some solution vectors + // (e.g. nth_darcy_solution_after_solving_pressure_part + // and + // n_minus_oneth_darcy_solution_after_solving_pressure_part) + // to store some solutions in previous time + // steps after the solution of the pressure + // and velocity part. + // + // The member functions within this class + // have been named so properly so that + // readers can easily understand what they + // are doing. + // + // Like step-31, this tutorial uses two + // DoFHandler objects for the darcy system + // (presure and velocity) and + // saturation. This is because we want it to + // run faster, which reasons have been + // described in step-31. + // + // There is yet another important thing: + // unlike step-31. this step uses one more + // ConstraintMatrix object called + // darcy_preconditioner_constraints. This + // constraint object only for assembling the + // matrix for darcy preconditioner includes + // hanging node constrants as well as + // Dirichlet boundary value + // constraints. Without this constraint + // object for the preconditioner, we cannot + // get the convergence results when we solve + // darcy linear system. + // + // The last one variable indicates whether + // the matrix needs to be rebuilt the next + // time the corresponding build functions are + // called. This allows us to move the + // corresponding if into the function and + // thereby keeping our main run() function + // clean and easy to read. + template + class TwoPhaseFlowProblem { public: - BlockSchurPreconditioner ( - const TrilinosWrappers::BlockSparseMatrix &S, - const InverseMatrix &Mpinv, - const PreconditionerA &Apreconditioner); - - void vmult (TrilinosWrappers::BlockVector &dst, - const TrilinosWrappers::BlockVector &src) const; + TwoPhaseFlowProblem (const unsigned int degree); + void run (); private: - const SmartPointer darcy_matrix; - const SmartPointer > m_inverse; - const PreconditionerA &a_preconditioner; - - mutable TrilinosWrappers::Vector tmp; + void setup_dofs (); + void assemble_darcy_preconditioner (); + void build_darcy_preconditioner (); + void assemble_darcy_system (); + void assemble_saturation_system (); + void assemble_saturation_matrix (); + void assemble_saturation_rhs (); + void assemble_saturation_rhs_cell_term (const FEValues &saturation_fe_values, + const FEValues &darcy_fe_values, + const std::vector &local_dof_indices, + const double global_u_infty, + const double global_S_variation, + const double global_Omega_diameter); + void assemble_saturation_rhs_boundary_term (const FEFaceValues &saturation_fe_face_values, + const FEFaceValues &darcy_fe_face_values, + const std::vector &local_dof_indices); + double get_maximal_velocity () const; + std::pair get_extrapolated_saturation_range () const; + void solve (); + bool determine_whether_to_solve_pressure_velocity_part () const; + void compute_refinement_indicators (Vector &indicator) const; + void refine_grid (const Vector &indicator); + void project_back_saturation (); + void output_results () const; + + static + double + compute_viscosity(const std::vector &old_saturation, + const std::vector &old_old_saturation, + const std::vector > &old_saturation_grads, + const std::vector > &old_old_saturation_grads, + const std::vector > &present_darcy_values, + const double global_u_infty, + const double global_S_variation, + const double global_Omega_diameter, + const double cell_diameter, + const double old_time_step, + const double viscosity); + + + const unsigned int degree; + + Triangulation triangulation; + + const unsigned int darcy_degree; + FESystem darcy_fe; + DoFHandler darcy_dof_handler; + ConstraintMatrix darcy_constraints; + + ConstraintMatrix darcy_preconditioner_constraints; + + TrilinosWrappers::BlockSparseMatrix darcy_matrix; + TrilinosWrappers::BlockSparseMatrix darcy_preconditioner_matrix; + + TrilinosWrappers::BlockVector darcy_solution; + TrilinosWrappers::BlockVector darcy_rhs; + + TrilinosWrappers::BlockVector nth_darcy_solution_after_solving_pressure_part; + TrilinosWrappers::BlockVector n_minus_oneth_darcy_solution_after_solving_pressure_part; + + const unsigned int saturation_degree; + FE_Q saturation_fe; + DoFHandler saturation_dof_handler; + ConstraintMatrix saturation_constraints; + + TrilinosWrappers::SparseMatrix saturation_matrix; + + TrilinosWrappers::Vector predictor_saturation_solution; + TrilinosWrappers::Vector saturation_solution; + TrilinosWrappers::Vector old_saturation_solution; + TrilinosWrappers::Vector old_old_saturation_solution; + TrilinosWrappers::Vector saturation_rhs; + + TrilinosWrappers::Vector nth_saturation_solution_after_solving_pressure_part; + + const unsigned int n_refinement_steps; + bool solve_pressure_velocity_part; + bool previous_solve_pressure_velocity_part; + + const double saturation_level; + const double saturation_value; + + double n_minus_oneth_time_step; + double cumulative_nth_time_step; + + double time_step; + double old_time_step; + unsigned int timestep_number; + double viscosity; + + std_cxx1x::shared_ptr Amg_preconditioner; + std_cxx1x::shared_ptr Mp_preconditioner; + + bool rebuild_saturation_matrix; }; + // @sect3{Pressure right hand side, Pressure boundary values and saturation initial value classes} - template - BlockSchurPreconditioner:: - BlockSchurPreconditioner(const TrilinosWrappers::BlockSparseMatrix &S, - const InverseMatrix &Mpinv, - const PreconditionerA &Apreconditioner) - : - darcy_matrix (&S), - m_inverse (&Mpinv), - a_preconditioner (Apreconditioner), - tmp (darcy_matrix->block(1,1).m()) - {} - - - template - void BlockSchurPreconditioner::vmult ( - TrilinosWrappers::BlockVector &dst, - const TrilinosWrappers::BlockVector &src) const + // This part is directly taken from step-21 + // so there is no need to repeat the same + // descriptions. + template + class PressureRightHandSide : public Function { - a_preconditioner.vmult (dst.block(0), src.block(0)); - darcy_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1)); - tmp *= -1; - m_inverse->vmult (dst.block(1), tmp); - } -} - - - // @sect3{The TwoPhaseFlowProblem class} - - // The definition of the class that defines - // the top-level logic of solving the - // time-dependent advection-dominated - // two-phase flow problem (or - // Buckley-Leverett problem - // [Buckley 1942]) is mainly based on - // three tutorial programs (step-21, step-31, - // step-33). The main difference is that, - // since adaptive operator splitting is - // considered, we need a bool-type variable - // solve_pressure_velocity_part to tell us - // when we need to solve the pressure and - // velocity part, need another bool-type - // variable - // previous_solve_pressure_velocity_part to - // determine if we have to cumulate - // micro-time steps that we need them to do - // extrapolation for the total velocity, and - // some solution vectors - // (e.g. nth_darcy_solution_after_solving_pressure_part - // and - // n_minus_oneth_darcy_solution_after_solving_pressure_part) - // to store some solutions in previous time - // steps after the solution of the pressure - // and velocity part. - // - // The member functions within this class - // have been named so properly so that - // readers can easily understand what they - // are doing. - // - // Like step-31, this tutorial uses two - // DoFHandler objects for the darcy system - // (presure and velocity) and - // saturation. This is because we want it to - // run faster, which reasons have been - // described in step-31. - // - // There is yet another important thing: - // unlike step-31. this step uses one more - // ConstraintMatrix object called - // darcy_preconditioner_constraints. This - // constraint object only for assembling the - // matrix for darcy preconditioner includes - // hanging node constrants as well as - // Dirichlet boundary value - // constraints. Without this constraint - // object for the preconditioner, we cannot - // get the convergence results when we solve - // darcy linear system. - // - // The last one variable indicates whether - // the matrix needs to be rebuilt the next - // time the corresponding build functions are - // called. This allows us to move the - // corresponding if into the function and - // thereby keeping our main run() function - // clean and easy to read. -template -class TwoPhaseFlowProblem -{ - public: - TwoPhaseFlowProblem (const unsigned int degree); - void run (); - - private: - void setup_dofs (); - void assemble_darcy_preconditioner (); - void build_darcy_preconditioner (); - void assemble_darcy_system (); - void assemble_saturation_system (); - void assemble_saturation_matrix (); - void assemble_saturation_rhs (); - void assemble_saturation_rhs_cell_term (const FEValues &saturation_fe_values, - const FEValues &darcy_fe_values, - const std::vector &local_dof_indices, - const double global_u_infty, - const double global_S_variation, - const double global_Omega_diameter); - void assemble_saturation_rhs_boundary_term (const FEFaceValues &saturation_fe_face_values, - const FEFaceValues &darcy_fe_face_values, - const std::vector &local_dof_indices); - double get_maximal_velocity () const; - std::pair get_extrapolated_saturation_range () const; - void solve (); - bool determine_whether_to_solve_pressure_velocity_part () const; - void compute_refinement_indicators (Vector &indicator) const; - void refine_grid (const Vector &indicator); - void project_back_saturation (); - void output_results () const; - - static - double - compute_viscosity(const std::vector &old_saturation, - const std::vector &old_old_saturation, - const std::vector > &old_saturation_grads, - const std::vector > &old_old_saturation_grads, - const std::vector > &present_darcy_values, - const double global_u_infty, - const double global_S_variation, - const double global_Omega_diameter, - const double cell_diameter, - const double old_time_step, - const double viscosity); - - - const unsigned int degree; - - Triangulation triangulation; - - const unsigned int darcy_degree; - FESystem darcy_fe; - DoFHandler darcy_dof_handler; - ConstraintMatrix darcy_constraints; - - ConstraintMatrix darcy_preconditioner_constraints; - - TrilinosWrappers::BlockSparseMatrix darcy_matrix; - TrilinosWrappers::BlockSparseMatrix darcy_preconditioner_matrix; - - TrilinosWrappers::BlockVector darcy_solution; - TrilinosWrappers::BlockVector darcy_rhs; - - TrilinosWrappers::BlockVector nth_darcy_solution_after_solving_pressure_part; - TrilinosWrappers::BlockVector n_minus_oneth_darcy_solution_after_solving_pressure_part; - - const unsigned int saturation_degree; - FE_Q saturation_fe; - DoFHandler saturation_dof_handler; - ConstraintMatrix saturation_constraints; - - TrilinosWrappers::SparseMatrix saturation_matrix; - - TrilinosWrappers::Vector predictor_saturation_solution; - TrilinosWrappers::Vector saturation_solution; - TrilinosWrappers::Vector old_saturation_solution; - TrilinosWrappers::Vector old_old_saturation_solution; - TrilinosWrappers::Vector saturation_rhs; - - TrilinosWrappers::Vector nth_saturation_solution_after_solving_pressure_part; - - const unsigned int n_refinement_steps; - bool solve_pressure_velocity_part; - bool previous_solve_pressure_velocity_part; - - const double saturation_level; - const double saturation_value; - - double n_minus_oneth_time_step; - double cumulative_nth_time_step; - - double time_step; - double old_time_step; - unsigned int timestep_number; - double viscosity; - - std_cxx1x::shared_ptr Amg_preconditioner; - std_cxx1x::shared_ptr Mp_preconditioner; - - bool rebuild_saturation_matrix; -}; - - - // @sect3{Pressure right hand side, Pressure boundary values and saturation initial value classes} - - // This part is directly taken from step-21 - // so there is no need to repeat the same - // descriptions. -template -class PressureRightHandSide : public Function -{ - public: - PressureRightHandSide () : Function(1) {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; -}; - - - -template -double -PressureRightHandSide::value (const Point &/*p*/, - const unsigned int /*component*/) const -{ - return 0; -} - - -template -class PressureBoundaryValues : public Function -{ - public: - PressureBoundaryValues () : Function(1) {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; -}; - - -template -double -PressureBoundaryValues::value (const Point &p, - const unsigned int /*component*/) const -{ - return 1-p[0]; -} - - -template -class SaturationBoundaryValues : public Function -{ - public: - SaturationBoundaryValues () : Function(1) {} + public: + PressureRightHandSide () : Function(1) {} - virtual double value (const Point &p, - const unsigned int component = 0) const; -}; + virtual double value (const Point &p, + const unsigned int component = 0) const; + }; -template -double -SaturationBoundaryValues::value (const Point &p, - const unsigned int /*component*/) const -{ - if (p[0] == 0) - return 1; - else + template + double + PressureRightHandSide::value (const Point &/*p*/, + const unsigned int /*component*/) const + { return 0; -} - - -template -class SaturationInitialValues : public Function -{ - public: - SaturationInitialValues () : Function(1) {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; - - virtual void vector_value (const Point &p, - Vector &value) const; - -}; + } -template -double -SaturationInitialValues::value (const Point &/*p*/, - const unsigned int /*component*/) const -{ - return 0; -} + template + class PressureBoundaryValues : public Function + { + public: + PressureBoundaryValues () : Function(1) {} + virtual double value (const Point &p, + const unsigned int component = 0) const; + }; -template -void -SaturationInitialValues::vector_value (const Point &p, - Vector &values) const -{ - for (unsigned int c=0; cn_components; ++c) - values(c) = SaturationInitialValues::value (p,c); -} + template + double + PressureBoundaryValues::value (const Point &p, + const unsigned int /*component*/) const + { + return 1-p[0]; + } - // @sect3{Permeability models} - // In this tutorial, we still use two - // permeability models previous used in - // step-21 so we refrain from excessive - // comments about them. But we want to note - // that if ones use the Random Medium model, - // they can change one parameter called the - // number of high-permeability regions/points - // to increase the amount of permeability in - // the computational domain. -namespace SingleCurvingCrack -{ template - class KInverse : public TensorFunction<2,dim> + class SaturationBoundaryValues : public Function { public: - KInverse () - : - TensorFunction<2,dim> () - {} + SaturationBoundaryValues () : Function(1) {} - virtual void value_list (const std::vector > &points, - std::vector > &values) const; + virtual double value (const Point &p, + const unsigned int component = 0) const; }; + template - void - KInverse::value_list (const std::vector > &points, - std::vector > &values) const + double + SaturationBoundaryValues::value (const Point &p, + const unsigned int /*component*/) const { - Assert (points.size() == values.size(), - ExcDimensionMismatch (points.size(), values.size())); - - for (unsigned int p=0; p - class KInverse : public TensorFunction<2,dim> + class SaturationInitialValues : public Function { public: - KInverse () - : - TensorFunction<2,dim> () - {} + SaturationInitialValues () : Function(1) {} - virtual void value_list (const std::vector > &points, - std::vector > &values) const; + virtual double value (const Point &p, + const unsigned int component = 0) const; - private: - static std::vector > centers; + virtual void vector_value (const Point &p, + Vector &value) const; - static std::vector > get_centers (); }; - template - std::vector > - KInverse::centers = KInverse::get_centers(); + double + SaturationInitialValues::value (const Point &/*p*/, + const unsigned int /*component*/) const + { + return 0; + } template - std::vector > - KInverse::get_centers () + void + SaturationInitialValues::vector_value (const Point &p, + Vector &values) const { - const unsigned int N = (dim == 2 ? - 40 : - (dim == 3 ? - 100 : - throw ExcNotImplemented())); - - std::vector > centers_list (N); - for (unsigned int i=0; i(rand())/RAND_MAX; - - return centers_list; + for (unsigned int c=0; cn_components; ++c) + values(c) = SaturationInitialValues::value (p,c); } + // @sect3{Permeability models} - template - void - KInverse::value_list (const std::vector > &points, - std::vector > &values) const + // In this tutorial, we still use two + // permeability models previous used in + // step-21 so we refrain from excessive + // comments about them. But we want to note + // that if ones use the Random Medium model, + // they can change one parameter called the + // number of high-permeability regions/points + // to increase the amount of permeability in + // the computational domain. + namespace SingleCurvingCrack { - Assert (points.size() == values.size(), - ExcDimensionMismatch (points.size(), values.size())); + template + class KInverse : public TensorFunction<2,dim> + { + public: + KInverse () + : + TensorFunction<2,dim> () + {} - for (unsigned int p=0; p > &points, + std::vector > &values) const; + }; - double permeability = 0; - for (unsigned int i=0; i + void + KInverse::value_list (const std::vector > &points, + std::vector > &values) const + { + Assert (points.size() == values.size(), + ExcDimensionMismatch (points.size(), values.size())); - for (unsigned int d=0; d + class KInverse : public TensorFunction<2,dim> + { + public: + KInverse () + : + TensorFunction<2,dim> () + {} - // The implementations of all the physical - // quantities such as total mobility - // $\lambda_t$ and fractional flow of water - // $F$ are taken from step-21 so again we - // don't have do any comment about them. -double mobility_inverse (const double S, - const double viscosity) -{ - return 1.0 /(1.0/viscosity * S * S + (1-S) * (1-S)); -} + virtual void value_list (const std::vector > &points, + std::vector > &values) const; -double f_saturation (const double S, - const double viscosity) -{ - return S*S /( S * S +viscosity * (1-S) * (1-S)); -} + private: + static std::vector > centers; -double get_fractional_flow_derivative (const double S, - const double viscosity) -{ - const double temp = ( S * S + viscosity * (1-S) * (1-S) ); + static std::vector > get_centers (); + }; - const double numerator = 2.0 * S * temp - - - S * S * - ( 2.0 * S - 2.0 * viscosity * (1-S) ); - const double denomerator = std::pow(temp, 2.0 ); - return numerator / denomerator; -} + template + std::vector > + KInverse::centers = KInverse::get_centers(); + + + template + std::vector > + KInverse::get_centers () + { + const unsigned int N = (dim == 2 ? + 40 : + (dim == 3 ? + 100 : + throw ExcNotImplemented())); + + std::vector > centers_list (N); + for (unsigned int i=0; i(rand())/RAND_MAX; + + return centers_list; + } - // @sect3{TwoPhaseFlowProblem::TwoPhaseFlowProblem} - - // The constructor of this class is an - // extension of the constructor in step-21 - // and step-31. We need to add the various - // variables that concern the saturation. As - // discussed in the introduction, we are - // going to use $Q_2 \times Q_1$ - // (Taylor-Hood) elements again for the darcy - // system, which element combination fulfills - // the Ladyzhenskaya-Babuska-Brezzi (LBB) - // conditions - // [Brezzi and Fortin 1991, Chen 2005], and $Q_1$ - // elements for the saturation. However, by - // using variables that store the polynomial - // degree of the darcy and temperature finite - // elements, it is easy to consistently - // modify the degree of the elements as well - // as all quadrature formulas used on them - // downstream. Moreover, we initialize the - // time stepping, variables related to - // operator splitting as well as the option - // for matrix assembly and preconditioning: -template -TwoPhaseFlowProblem::TwoPhaseFlowProblem (const unsigned int degree) - : - degree (degree), - darcy_degree (degree), - darcy_fe (FE_Q(darcy_degree+1), dim, - FE_Q(darcy_degree), 1), - darcy_dof_handler (triangulation), - - saturation_degree (degree), - saturation_fe (saturation_degree), - saturation_dof_handler (triangulation), - - n_refinement_steps (4), - solve_pressure_velocity_part (false), - previous_solve_pressure_velocity_part (false), - - saturation_level (2), - saturation_value (0.5), - - time_step (0), - old_time_step (0), - viscosity (0.2), - - rebuild_saturation_matrix (true) -{} - - - // @sect3{TwoPhaseFlowProblem::setup_dofs} - - // This is the function that sets up the - // DoFHandler objects we have here (one for - // the darcy part and one for the saturation - // part) as well as set to the right sizes - // the various objects required for the - // linear algebra in this program. Its basic - // operations are similar to what authors in - // step-31 did. - // - // The body of the function first enumerates - // all degrees of freedom for the darcy and - // saturation systems. For the darcy part, - // degrees of freedom are then sorted to - // ensure that velocities precede pressure - // DoFs so that we can partition the darcy - // matrix into a $2 \times 2$ matrix. Like - // step-31, the present step does not perform - // any additional DoF renumbering. - // - // Then, we need to incorporate hanging node - // constraints and Dirichlet boundary value - // constraints into - // darcy_preconditioner_constraints. However, - // this constraints are only set to the - // pressure component since the Schur - // complement preconditioner that corresponds - // to the porous media flow operator in - // non-mixed form, $-\nabla \cdot [\mathbf K - // \lambda_t(S)]\nabla$. Therefore, we use a - // component_mask that filters out the - // velocity component, so that the - // condensation is performed on pressure - // degrees of freedom only. - // - // After having done so, we count the number - // of degrees of freedom in the various - // blocks: - // - // The next step is to create the sparsity - // pattern for the darcy and saturation - // system matrices as well as the - // preconditioner matrix from which we build - // the darcy preconditioner. As in step-31, - // we choose to create the pattern not as in - // the first few tutorial programs, but by - // using the blocked version of - // CompressedSimpleSparsityPattern. The - // reason for doing this is mainly memory, - // that is, the SparsityPattern class would - // consume too much memory when used in three - // spatial dimensions as we intend to do for - // this program. So, for this, we follow the - // same way as step-31 did and we don't have - // to repeat descriptions again for the rest - // of the member function. -template -void TwoPhaseFlowProblem::setup_dofs () -{ - std::vector darcy_block_component (dim+1,0); - darcy_block_component[dim] = 1; - { - darcy_dof_handler.distribute_dofs (darcy_fe); - DoFRenumbering::Cuthill_McKee (darcy_dof_handler); - DoFRenumbering::component_wise (darcy_dof_handler, darcy_block_component); - darcy_constraints.clear (); - DoFTools::make_hanging_node_constraints (darcy_dof_handler, darcy_constraints); - darcy_constraints.close (); + template + void + KInverse::value_list (const std::vector > &points, + std::vector > &values) const + { + Assert (points.size() == values.size(), + ExcDimensionMismatch (points.size(), values.size())); + + for (unsigned int p=0; p component_mask (dim+1, false); - component_mask[dim] = true; + double get_fractional_flow_derivative (const double S, + const double viscosity) + { + const double temp = ( S * S + viscosity * (1-S) * (1-S) ); + const double numerator = 2.0 * S * temp + - + S * S * + ( 2.0 * S - 2.0 * viscosity * (1-S) ); - DoFTools::make_hanging_node_constraints (darcy_dof_handler, darcy_preconditioner_constraints); - DoFTools::make_zero_boundary_constraints (darcy_dof_handler, darcy_preconditioner_constraints, component_mask); + const double denomerator = std::pow(temp, 2.0 ); - darcy_preconditioner_constraints.close (); + return numerator / denomerator; } - std::vector darcy_dofs_per_block (2); - DoFTools::count_dofs_per_block (darcy_dof_handler, darcy_dofs_per_block, darcy_block_component); - const unsigned int n_u = darcy_dofs_per_block[0], - n_p = darcy_dofs_per_block[1], - n_s = saturation_dof_handler.n_dofs(); - - std::cout << "Number of active cells: " - << triangulation.n_active_cells() - << " (on " - << triangulation.n_levels() - << " levels)" - << std::endl - << "Number of degrees of freedom: " - << n_u + n_p + n_s - << " (" << n_u << '+' << n_p << '+'<< n_s <<')' - << std::endl - << std::endl; + // @sect3{TwoPhaseFlowProblem::TwoPhaseFlowProblem} + + // The constructor of this class is an + // extension of the constructor in step-21 + // and step-31. We need to add the various + // variables that concern the saturation. As + // discussed in the introduction, we are + // going to use $Q_2 \times Q_1$ + // (Taylor-Hood) elements again for the darcy + // system, which element combination fulfills + // the Ladyzhenskaya-Babuska-Brezzi (LBB) + // conditions + // [Brezzi and Fortin 1991, Chen 2005], and $Q_1$ + // elements for the saturation. However, by + // using variables that store the polynomial + // degree of the darcy and temperature finite + // elements, it is easy to consistently + // modify the degree of the elements as well + // as all quadrature formulas used on them + // downstream. Moreover, we initialize the + // time stepping, variables related to + // operator splitting as well as the option + // for matrix assembly and preconditioning: + template + TwoPhaseFlowProblem::TwoPhaseFlowProblem (const unsigned int degree) + : + degree (degree), + darcy_degree (degree), + darcy_fe (FE_Q(darcy_degree+1), dim, + FE_Q(darcy_degree), 1), + darcy_dof_handler (triangulation), - { - darcy_matrix.clear (); + saturation_degree (degree), + saturation_fe (saturation_degree), + saturation_dof_handler (triangulation), - BlockCompressedSimpleSparsityPattern csp (2,2); + n_refinement_steps (4), + solve_pressure_velocity_part (false), + previous_solve_pressure_velocity_part (false), - csp.block(0,0).reinit (n_u, n_u); - csp.block(0,1).reinit (n_u, n_p); - csp.block(1,0).reinit (n_p, n_u); - csp.block(1,1).reinit (n_p, n_p); + saturation_level (2), + saturation_value (0.5), - csp.collect_sizes (); + time_step (0), + old_time_step (0), + viscosity (0.2), - Table<2,DoFTools::Coupling> coupling (dim+1, dim+1); + rebuild_saturation_matrix (true) + {} - for (unsigned int c=0; c::setup_dofs} + + // This is the function that sets up the + // DoFHandler objects we have here (one for + // the darcy part and one for the saturation + // part) as well as set to the right sizes + // the various objects required for the + // linear algebra in this program. Its basic + // operations are similar to what authors in + // step-31 did. + // + // The body of the function first enumerates + // all degrees of freedom for the darcy and + // saturation systems. For the darcy part, + // degrees of freedom are then sorted to + // ensure that velocities precede pressure + // DoFs so that we can partition the darcy + // matrix into a $2 \times 2$ matrix. Like + // step-31, the present step does not perform + // any additional DoF renumbering. + // + // Then, we need to incorporate hanging node + // constraints and Dirichlet boundary value + // constraints into + // darcy_preconditioner_constraints. However, + // this constraints are only set to the + // pressure component since the Schur + // complement preconditioner that corresponds + // to the porous media flow operator in + // non-mixed form, $-\nabla \cdot [\mathbf K + // \lambda_t(S)]\nabla$. Therefore, we use a + // component_mask that filters out the + // velocity component, so that the + // condensation is performed on pressure + // degrees of freedom only. + // + // After having done so, we count the number + // of degrees of freedom in the various + // blocks: + // + // The next step is to create the sparsity + // pattern for the darcy and saturation + // system matrices as well as the + // preconditioner matrix from which we build + // the darcy preconditioner. As in step-31, + // we choose to create the pattern not as in + // the first few tutorial programs, but by + // using the blocked version of + // CompressedSimpleSparsityPattern. The + // reason for doing this is mainly memory, + // that is, the SparsityPattern class would + // consume too much memory when used in three + // spatial dimensions as we intend to do for + // this program. So, for this, we follow the + // same way as step-31 did and we don't have + // to repeat descriptions again for the rest + // of the member function. + template + void TwoPhaseFlowProblem::setup_dofs () + { + std::vector darcy_block_component (dim+1,0); + darcy_block_component[dim] = 1; + { + darcy_dof_handler.distribute_dofs (darcy_fe); + DoFRenumbering::Cuthill_McKee (darcy_dof_handler); + DoFRenumbering::component_wise (darcy_dof_handler, darcy_block_component); - DoFTools::make_sparsity_pattern (darcy_dof_handler, coupling, csp, - darcy_constraints, false); + darcy_constraints.clear (); + DoFTools::make_hanging_node_constraints (darcy_dof_handler, darcy_constraints); + darcy_constraints.close (); + } + { + saturation_dof_handler.distribute_dofs (saturation_fe); - darcy_matrix.reinit (csp); - } + saturation_constraints.clear (); + DoFTools::make_hanging_node_constraints (saturation_dof_handler, saturation_constraints); + saturation_constraints.close (); + } + { + darcy_preconditioner_constraints.clear (); - { - Amg_preconditioner.reset (); - Mp_preconditioner.reset (); - darcy_preconditioner_matrix.clear (); + std::vector component_mask (dim+1, false); + component_mask[dim] = true; - BlockCompressedSimpleSparsityPattern csp (2,2); - csp.block(0,0).reinit (n_u, n_u); - csp.block(0,1).reinit (n_u, n_p); - csp.block(1,0).reinit (n_p, n_u); - csp.block(1,1).reinit (n_p, n_p); + DoFTools::make_hanging_node_constraints (darcy_dof_handler, darcy_preconditioner_constraints); + DoFTools::make_zero_boundary_constraints (darcy_dof_handler, darcy_preconditioner_constraints, component_mask); - csp.collect_sizes (); + darcy_preconditioner_constraints.close (); + } - Table<2,DoFTools::Coupling> coupling (dim+1, dim+1); - for (unsigned int c=0; c darcy_dofs_per_block (2); + DoFTools::count_dofs_per_block (darcy_dof_handler, darcy_dofs_per_block, darcy_block_component); + const unsigned int n_u = darcy_dofs_per_block[0], + n_p = darcy_dofs_per_block[1], + n_s = saturation_dof_handler.n_dofs(); + + std::cout << "Number of active cells: " + << triangulation.n_active_cells() + << " (on " + << triangulation.n_levels() + << " levels)" + << std::endl + << "Number of degrees of freedom: " + << n_u + n_p + n_s + << " (" << n_u << '+' << n_p << '+'<< n_s <<')' + << std::endl + << std::endl; - darcy_preconditioner_matrix.reinit (csp); - } + { + darcy_matrix.clear (); + BlockCompressedSimpleSparsityPattern csp (2,2); - { - saturation_matrix.clear (); + csp.block(0,0).reinit (n_u, n_u); + csp.block(0,1).reinit (n_u, n_p); + csp.block(1,0).reinit (n_p, n_u); + csp.block(1,1).reinit (n_p, n_p); - CompressedSimpleSparsityPattern csp (n_s, n_s); + csp.collect_sizes (); - DoFTools::make_sparsity_pattern (saturation_dof_handler, csp, - saturation_constraints, false); + Table<2,DoFTools::Coupling> coupling (dim+1, dim+1); + for (unsigned int c=0; c coupling (dim+1, dim+1); + for (unsigned int c=0; c::assemble_darcy_preconditioner} - - // This function assembles the matrix we use - // for preconditioning the darcy system. What - // we need are a vector matrix weighted by - // $\left(\mathbf{K} \lambda_t\right)^{-1}$ - // on the velocity components and a mass - // matrix weighted by $\left(\mathbf{K} - // \lambda_t\right)$ on the pressure - // component. We start by generating a - // quadrature object of appropriate order, - // the FEValues object that can give values - // and gradients at the quadrature points - // (together with quadrature weights). Next - // we create data structures for the cell - // matrix and the relation between local and - // global DoFs. The vectors phi_u and - // grad_phi_p are going to hold the values of - // the basis functions in order to faster - // build up the local matrices, as was - // already done in step-22. Before we start - // the loop over all active cells, we have to - // specify which components are pressure and - // which are velocity. - // - // The creation of the local matrix is rather - // simple. There are only a term weighted by - // $\left(\mathbf{K} \lambda_t\right)^{-1}$ - // (on the velocity) and a mass matrix - // weighted by $\left(\mathbf{K} - // \lambda_t\right)$ to be generated, so the - // creation of the local matrix is done in - // two lines. Once the local matrix is ready - // (loop over rows and columns in the local - // matrix on each quadrature point), we get - // the local DoF indices and write the local - // information into the global matrix. We do - // this by directly applying the constraints - // (i.e. darcy_preconditioner_constraints) - // from hanging nodes locally and Dirichlet - // boundary conditions with zero values. By - // doing so, we don't have to do that - // afterwards, and we don't also write into - // entries of the matrix that will actually - // be set to zero again later when - // eliminating constraints. -template -void -TwoPhaseFlowProblem::assemble_darcy_preconditioner () -{ - std::cout << " Rebuilding darcy preconditioner..." << std::endl; - - darcy_preconditioner_matrix = 0; - - const QGauss quadrature_formula(darcy_degree+2); - FEValues darcy_fe_values (darcy_fe, quadrature_formula, - update_JxW_values | - update_values | - update_gradients | - update_quadrature_points); - FEValues saturation_fe_values (saturation_fe, quadrature_formula, - update_values); - - const unsigned int dofs_per_cell = darcy_fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); - - const RandomMedium::KInverse k_inverse; -// const SingleCurvingCrack::KInverse k_inverse; - std::vector > k_inverse_values (n_q_points); - Tensor<2,dim> k_value; - - std::vector old_saturation_values (n_q_points); - - FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); - std::vector local_dof_indices (dofs_per_cell); - - std::vector > phi_u (dofs_per_cell); - std::vector > grad_phi_p (dofs_per_cell); - - const FEValuesExtractors::Vector velocities (0); - const FEValuesExtractors::Scalar pressure (dim); - - typename DoFHandler::active_cell_iterator - cell = darcy_dof_handler.begin_active(), - endc = darcy_dof_handler.end(); - typename DoFHandler::active_cell_iterator - saturation_cell = saturation_dof_handler.begin_active(); - - for (; cell!=endc; ++cell, ++saturation_cell) { - darcy_fe_values.reinit (cell); - saturation_fe_values.reinit (saturation_cell); - - local_matrix = 0; - - saturation_fe_values.get_function_values (old_saturation_solution, old_saturation_values); - - k_inverse.value_list (darcy_fe_values.get_quadrature_points(), - k_inverse_values); - - for (unsigned int q=0; qget_dof_indices (local_dof_indices); - darcy_preconditioner_constraints.distribute_local_to_global (local_matrix, - local_dof_indices, - darcy_preconditioner_matrix); - } -} + saturation_matrix.clear (); + CompressedSimpleSparsityPattern csp (n_s, n_s); - // @sect3{TwoPhaseFlowProblem::build_darcy_preconditioner} - - // This function generates the inner - // preconditioners that are going to be used - // for the Schur complement block - // preconditioner. The preconditioners need - // to be regenerated at every saturation time - // step since they contain the independent - // variables saturation $S$ with time. - // - // Next, we set up the preconditioner for the - // velocity-velocity matrix - // $\mathbf{M}^{\mathbf{u}}$ and the Schur - // complement $\mathbf{S}$. As explained in - // the introduction, we are going to use an - // IC preconditioner based on a vector matrix - // (which is spectrally close to the darcy - // matrix $\mathbf{M}^{\mathbf{u}}$) and - // another based on a Laplace vector matrix - // (which is spectrally close to the - // non-mixed pressure matrix - // $\mathbf{S}$). Usually, the - // TrilinosWrappers::PreconditionIC class can - // be seen as a good black-box preconditioner - // which does not need any special knowledge. -template -void -TwoPhaseFlowProblem::build_darcy_preconditioner () -{ - assemble_darcy_preconditioner (); + DoFTools::make_sparsity_pattern (saturation_dof_handler, csp, + saturation_constraints, false); - Amg_preconditioner = std_cxx1x::shared_ptr - (new TrilinosWrappers::PreconditionIC()); - Amg_preconditioner->initialize(darcy_preconditioner_matrix.block(0,0)); - Mp_preconditioner = std_cxx1x::shared_ptr - (new TrilinosWrappers::PreconditionIC()); - Mp_preconditioner->initialize(darcy_preconditioner_matrix.block(1,1)); + saturation_matrix.reinit (csp); + } -} + darcy_solution.reinit (2); + darcy_solution.block(0).reinit (n_u); + darcy_solution.block(1).reinit (n_p); + darcy_solution.collect_sizes (); + nth_darcy_solution_after_solving_pressure_part.reinit (2); + nth_darcy_solution_after_solving_pressure_part.block(0).reinit (n_u); + nth_darcy_solution_after_solving_pressure_part.block(1).reinit (n_p); + nth_darcy_solution_after_solving_pressure_part.collect_sizes (); - // @sect3{TwoPhaseFlowProblem::assemble_darcy_system} - - // This is the function that assembles the - // linear system for the darcy system. - // - // Regarding the technical details of - // implementation, the procedures are similar - // to those in step-22 and step-31 we reset - // matrix and vector, create a quadrature - // formula on the cells, and then create the - // respective FEValues object. For the update - // flags, we require basis function - // derivatives only in case of a full - // assembly, since they are not needed for - // the right hand side; as always, choosing - // the minimal set of flags depending on what - // is currently needed makes the call to - // FEValues::reinit further down in the - // program more efficient. - // - // There is one thing that needs to be - // commented ¡V since we have a separate - // finite element and DoFHandler for the - // saturation, we need to generate a second - // FEValues object for the proper evaluation - // of the saturation solution. This isn't too - // complicated to realize here: just use the - // saturation structures and set an update - // flag for the basis function values which - // we need for evaluation of the saturation - // solution. The only important part to - // remember here is that the same quadrature - // formula is used for both FEValues objects - // to ensure that we get matching information - // when we loop over the quadrature points of - // the two objects. - // - // The declarations proceed with some - // shortcuts for array sizes, the creation of - // the local matrix, right hand side as well - // as the vector for the indices of the local - // dofs compared to the global system. - // - // Note that in its present form, the - // function uses the permeability implemented - // in the RandomMedium::KInverse - // class. Switching to the single curved - // crack permeability function is as simple - // as just changing the namespace name. - // - // Here's the an important step: we have to - // get the values of the saturation function - // of the previous time step at the - // quadrature points. To this end, we can use - // the FEValues::get_function_values - // (previously already used in step-9, - // step-14 and step-15), a function that - // takes a solution vector and returns a list - // of function values at the quadrature - // points of the present cell. In fact, it - // returns the complete vector-valued - // solution at each quadrature point, - // i.e. not only the saturation but also the - // velocities and pressure: - // - // Next we need a vector that will contain - // the values of the saturation solution at - // the previous time level at the quadrature - // points to assemble the source term in the - // right hand side of the momentum - // equation. Let's call this vector - // old_saturation_values. - // - // The set of vectors we create next hold the - // evaluations of the basis functions as well - // as their gradients and symmetrized - // gradients that will be used for creating - // the matrices. Putting these into their own - // arrays rather than asking the FEValues - // object for this information each time it - // is needed is an optimization to accelerate - // the assembly process, see step-22 for - // details. - // - // The last two declarations are used to - // extract the individual blocks (velocity, - // pressure, saturation) from the total FE - // system. - // - // Now start the loop over all cells in the - // problem. We are working on two different - // DoFHandlers for this assembly routine, so - // we must have two different cell iterators - // for the two objects in use. This might - // seem a bit peculiar, since both the darcy - // system and the saturation system use the - // same grid, but that's the only way to keep - // degrees of freedom in sync. The first - // statements within the loop are again all - // very familiar, doing the update of the - // finite element data as specified by the - // update flags, zeroing out the local arrays - // and getting the values of the old solution - // at the quadrature points. Then we are - // ready to loop over the quadrature points - // on the cell. - // - // Once this is done, we start the loop over - // the rows and columns of the local matrix - // and feed the matrix with the relevant - // products. - // - // The last step in the loop over all cells - // is to enter the local contributions into - // the global matrix and vector structures to - // the positions specified in - // local_dof_indices. Again, we let the - // ConstraintMatrix class do the insertion of - // the cell matrix elements to the global - // matrix, which already condenses the - // hanging node constraints. -template -void TwoPhaseFlowProblem::assemble_darcy_system () -{ - darcy_matrix = 0; - darcy_rhs = 0; + n_minus_oneth_darcy_solution_after_solving_pressure_part.reinit (2); + n_minus_oneth_darcy_solution_after_solving_pressure_part.block(0).reinit (n_u); + n_minus_oneth_darcy_solution_after_solving_pressure_part.block(1).reinit (n_p); + n_minus_oneth_darcy_solution_after_solving_pressure_part.collect_sizes (); - QGauss quadrature_formula(darcy_degree+2); - QGauss face_quadrature_formula(darcy_degree+2); + darcy_rhs.reinit (2); + darcy_rhs.block(0).reinit (n_u); + darcy_rhs.block(1).reinit (n_p); + darcy_rhs.collect_sizes (); - FEValues darcy_fe_values (darcy_fe, quadrature_formula, - update_values | update_gradients | - update_quadrature_points | update_JxW_values); + predictor_saturation_solution.reinit (n_s); + saturation_solution.reinit (n_s); + old_saturation_solution.reinit (n_s); + old_old_saturation_solution.reinit (n_s); - FEValues saturation_fe_values (saturation_fe, quadrature_formula, - update_values); + nth_saturation_solution_after_solving_pressure_part.reinit (n_s); + + saturation_rhs.reinit (n_s); + } - FEFaceValues darcy_fe_face_values (darcy_fe, face_quadrature_formula, - update_values | update_normal_vectors | - update_quadrature_points | update_JxW_values); - const unsigned int dofs_per_cell = darcy_fe.dofs_per_cell; + // @sect3{TwoPhaseFlowProblem::assemble_darcy_preconditioner} + + // This function assembles the matrix we use + // for preconditioning the darcy system. What + // we need are a vector matrix weighted by + // $\left(\mathbf{K} \lambda_t\right)^{-1}$ + // on the velocity components and a mass + // matrix weighted by $\left(\mathbf{K} + // \lambda_t\right)$ on the pressure + // component. We start by generating a + // quadrature object of appropriate order, + // the FEValues object that can give values + // and gradients at the quadrature points + // (together with quadrature weights). Next + // we create data structures for the cell + // matrix and the relation between local and + // global DoFs. The vectors phi_u and + // grad_phi_p are going to hold the values of + // the basis functions in order to faster + // build up the local matrices, as was + // already done in step-22. Before we start + // the loop over all active cells, we have to + // specify which components are pressure and + // which are velocity. + // + // The creation of the local matrix is rather + // simple. There are only a term weighted by + // $\left(\mathbf{K} \lambda_t\right)^{-1}$ + // (on the velocity) and a mass matrix + // weighted by $\left(\mathbf{K} + // \lambda_t\right)$ to be generated, so the + // creation of the local matrix is done in + // two lines. Once the local matrix is ready + // (loop over rows and columns in the local + // matrix on each quadrature point), we get + // the local DoF indices and write the local + // information into the global matrix. We do + // this by directly applying the constraints + // (i.e. darcy_preconditioner_constraints) + // from hanging nodes locally and Dirichlet + // boundary conditions with zero values. By + // doing so, we don't have to do that + // afterwards, and we don't also write into + // entries of the matrix that will actually + // be set to zero again later when + // eliminating constraints. + template + void + TwoPhaseFlowProblem::assemble_darcy_preconditioner () + { + std::cout << " Rebuilding darcy preconditioner..." << std::endl; - const unsigned int n_q_points = quadrature_formula.size(); - const unsigned int n_face_q_points = face_quadrature_formula.size(); + darcy_preconditioner_matrix = 0; - FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); - Vector local_rhs (dofs_per_cell); + const QGauss quadrature_formula(darcy_degree+2); + FEValues darcy_fe_values (darcy_fe, quadrature_formula, + update_JxW_values | + update_values | + update_gradients | + update_quadrature_points); + FEValues saturation_fe_values (saturation_fe, quadrature_formula, + update_values); - std::vector local_dof_indices (dofs_per_cell); + const unsigned int dofs_per_cell = darcy_fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); - const PressureRightHandSide pressure_right_hand_side; - const PressureBoundaryValues pressure_boundary_values; - const RandomMedium::KInverse k_inverse; + const RandomMedium::KInverse k_inverse; // const SingleCurvingCrack::KInverse k_inverse; - std::vector pressure_rhs_values (n_q_points); - std::vector boundary_values (n_face_q_points); - std::vector > k_inverse_values (n_q_points); + std::vector > k_inverse_values (n_q_points); + Tensor<2,dim> k_value; - std::vector old_saturation_values (n_q_points); + std::vector old_saturation_values (n_q_points); - std::vector > phi_u (dofs_per_cell); - std::vector div_phi_u (dofs_per_cell); - std::vector phi_p (dofs_per_cell); + FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); - const FEValuesExtractors::Vector velocities (0); - const FEValuesExtractors::Scalar pressure (dim); + std::vector > phi_u (dofs_per_cell); + std::vector > grad_phi_p (dofs_per_cell); - typename DoFHandler::active_cell_iterator - cell = darcy_dof_handler.begin_active(), - endc = darcy_dof_handler.end(); - typename DoFHandler::active_cell_iterator - saturation_cell = saturation_dof_handler.begin_active(); + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); - for (; cell!=endc; ++cell, ++saturation_cell) - { - darcy_fe_values.reinit (cell); - saturation_fe_values.reinit (saturation_cell); - - local_matrix = 0; - local_rhs = 0; - - saturation_fe_values.get_function_values (old_saturation_solution, old_saturation_values); - - pressure_right_hand_side.value_list (darcy_fe_values.get_quadrature_points(), - pressure_rhs_values); - k_inverse.value_list (darcy_fe_values.get_quadrature_points(), - k_inverse_values); - - for (unsigned int q=0; q::active_cell_iterator + cell = darcy_dof_handler.begin_active(), + endc = darcy_dof_handler.end(); + typename DoFHandler::active_cell_iterator + saturation_cell = saturation_dof_handler.begin_active(); - for (unsigned int face_no=0; - face_no::faces_per_cell; - ++face_no) - if (cell->at_boundary(face_no)) - { - darcy_fe_face_values.reinit (cell, face_no); - - pressure_boundary_values - .value_list (darcy_fe_face_values.get_quadrature_points(), - boundary_values); - - for (unsigned int q=0; q - phi_i_u = darcy_fe_face_values[velocities].value (i, q); - - local_rhs(i) += -(phi_i_u * - darcy_fe_face_values.normal_vector(q) * - boundary_values[q] * - darcy_fe_face_values.JxW(q)); - } - } + for (; cell!=endc; ++cell, ++saturation_cell) + { + darcy_fe_values.reinit (cell); + saturation_fe_values.reinit (saturation_cell); + + local_matrix = 0; + + saturation_fe_values.get_function_values (old_saturation_solution, old_saturation_values); + + k_inverse.value_list (darcy_fe_values.get_quadrature_points(), + k_inverse_values); + + for (unsigned int q=0; qget_dof_indices (local_dof_indices); + darcy_preconditioner_constraints.distribute_local_to_global (local_matrix, + local_dof_indices, + darcy_preconditioner_matrix); + } + } - for (unsigned int i=0; iget_dof_indices (local_dof_indices); + // @sect3{TwoPhaseFlowProblem::build_darcy_preconditioner} + + // This function generates the inner + // preconditioners that are going to be used + // for the Schur complement block + // preconditioner. The preconditioners need + // to be regenerated at every saturation time + // step since they contain the independent + // variables saturation $S$ with time. + // + // Next, we set up the preconditioner for the + // velocity-velocity matrix + // $\mathbf{M}^{\mathbf{u}}$ and the Schur + // complement $\mathbf{S}$. As explained in + // the introduction, we are going to use an + // IC preconditioner based on a vector matrix + // (which is spectrally close to the darcy + // matrix $\mathbf{M}^{\mathbf{u}}$) and + // another based on a Laplace vector matrix + // (which is spectrally close to the + // non-mixed pressure matrix + // $\mathbf{S}$). Usually, the + // TrilinosWrappers::PreconditionIC class can + // be seen as a good black-box preconditioner + // which does not need any special knowledge. + template + void + TwoPhaseFlowProblem::build_darcy_preconditioner () + { + assemble_darcy_preconditioner (); - darcy_constraints.distribute_local_to_global (local_matrix, - local_rhs, - local_dof_indices, - darcy_matrix, - darcy_rhs); + Amg_preconditioner = std_cxx1x::shared_ptr + (new TrilinosWrappers::PreconditionIC()); + Amg_preconditioner->initialize(darcy_preconditioner_matrix.block(0,0)); - } -} + Mp_preconditioner = std_cxx1x::shared_ptr + (new TrilinosWrappers::PreconditionIC()); + Mp_preconditioner->initialize(darcy_preconditioner_matrix.block(1,1)); + } - // @sect3{TwoPhaseFlowProblem::assemble_saturation_system} - - // This function is to assemble the linear - // system for the saturation transport - // equation. It includes two member - // functions: assemble_saturation_matrix () - // and assemble_saturation_rhs (). The former - // function that assembles the saturation - // left hand side needs to be changed only - // when grids have been changed since the - // matrix is filled only with basis - // functions. However, the latter that - // assembles the right hand side must be - // changed at every saturation time step - // since it depends on an unknown variable - // saturation. -template -void TwoPhaseFlowProblem::assemble_saturation_system () -{ - if ( rebuild_saturation_matrix == true ) - { - saturation_matrix = 0; - assemble_saturation_matrix (); - } - saturation_rhs = 0; - assemble_saturation_rhs (); -} + // @sect3{TwoPhaseFlowProblem::assemble_darcy_system} + + // This is the function that assembles the + // linear system for the darcy system. + // + // Regarding the technical details of + // implementation, the procedures are similar + // to those in step-22 and step-31 we reset + // matrix and vector, create a quadrature + // formula on the cells, and then create the + // respective FEValues object. For the update + // flags, we require basis function + // derivatives only in case of a full + // assembly, since they are not needed for + // the right hand side; as always, choosing + // the minimal set of flags depending on what + // is currently needed makes the call to + // FEValues::reinit further down in the + // program more efficient. + // + // There is one thing that needs to be + // commented ¡V since we have a separate + // finite element and DoFHandler for the + // saturation, we need to generate a second + // FEValues object for the proper evaluation + // of the saturation solution. This isn't too + // complicated to realize here: just use the + // saturation structures and set an update + // flag for the basis function values which + // we need for evaluation of the saturation + // solution. The only important part to + // remember here is that the same quadrature + // formula is used for both FEValues objects + // to ensure that we get matching information + // when we loop over the quadrature points of + // the two objects. + // + // The declarations proceed with some + // shortcuts for array sizes, the creation of + // the local matrix, right hand side as well + // as the vector for the indices of the local + // dofs compared to the global system. + // + // Note that in its present form, the + // function uses the permeability implemented + // in the RandomMedium::KInverse + // class. Switching to the single curved + // crack permeability function is as simple + // as just changing the namespace name. + // + // Here's the an important step: we have to + // get the values of the saturation function + // of the previous time step at the + // quadrature points. To this end, we can use + // the FEValues::get_function_values + // (previously already used in step-9, + // step-14 and step-15), a function that + // takes a solution vector and returns a list + // of function values at the quadrature + // points of the present cell. In fact, it + // returns the complete vector-valued + // solution at each quadrature point, + // i.e. not only the saturation but also the + // velocities and pressure: + // + // Next we need a vector that will contain + // the values of the saturation solution at + // the previous time level at the quadrature + // points to assemble the source term in the + // right hand side of the momentum + // equation. Let's call this vector + // old_saturation_values. + // + // The set of vectors we create next hold the + // evaluations of the basis functions as well + // as their gradients and symmetrized + // gradients that will be used for creating + // the matrices. Putting these into their own + // arrays rather than asking the FEValues + // object for this information each time it + // is needed is an optimization to accelerate + // the assembly process, see step-22 for + // details. + // + // The last two declarations are used to + // extract the individual blocks (velocity, + // pressure, saturation) from the total FE + // system. + // + // Now start the loop over all cells in the + // problem. We are working on two different + // DoFHandlers for this assembly routine, so + // we must have two different cell iterators + // for the two objects in use. This might + // seem a bit peculiar, since both the darcy + // system and the saturation system use the + // same grid, but that's the only way to keep + // degrees of freedom in sync. The first + // statements within the loop are again all + // very familiar, doing the update of the + // finite element data as specified by the + // update flags, zeroing out the local arrays + // and getting the values of the old solution + // at the quadrature points. Then we are + // ready to loop over the quadrature points + // on the cell. + // + // Once this is done, we start the loop over + // the rows and columns of the local matrix + // and feed the matrix with the relevant + // products. + // + // The last step in the loop over all cells + // is to enter the local contributions into + // the global matrix and vector structures to + // the positions specified in + // local_dof_indices. Again, we let the + // ConstraintMatrix class do the insertion of + // the cell matrix elements to the global + // matrix, which already condenses the + // hanging node constraints. + template + void TwoPhaseFlowProblem::assemble_darcy_system () + { + darcy_matrix = 0; + darcy_rhs = 0; + QGauss quadrature_formula(darcy_degree+2); + QGauss face_quadrature_formula(darcy_degree+2); + FEValues darcy_fe_values (darcy_fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); - // @sect3{TwoPhaseFlowProblem::assemble_saturation_matrix} + FEValues saturation_fe_values (saturation_fe, quadrature_formula, + update_values); - // This function is easily understood since - // it only forms a simple mass matrix for the - // left hand side of the saturation linear - // system by basis functions phi_i_s and - // phi_j_s only. Finally, as usual, we enter - // the local contribution into the global - // matrix by specifying the position in - // local_dof_indices. This is done by letting - // the ConstraintMatrix class do the - // insertion of the cell matrix elements to - // the global matrix, which already condenses - // the hanging node constraints. -template -void TwoPhaseFlowProblem::assemble_saturation_matrix () -{ - QGauss quadrature_formula(saturation_degree+2); + FEFaceValues darcy_fe_face_values (darcy_fe, face_quadrature_formula, + update_values | update_normal_vectors | + update_quadrature_points | update_JxW_values); - FEValues saturation_fe_values (saturation_fe, quadrature_formula, - update_values | update_JxW_values); + const unsigned int dofs_per_cell = darcy_fe.dofs_per_cell; - const unsigned int dofs_per_cell = saturation_fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + const unsigned int n_face_q_points = face_quadrature_formula.size(); - const unsigned int n_q_points = quadrature_formula.size(); + FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); + Vector local_rhs (dofs_per_cell); - FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); - Vector local_rhs (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); - std::vector local_dof_indices (dofs_per_cell); + const PressureRightHandSide pressure_right_hand_side; + const PressureBoundaryValues pressure_boundary_values; + const RandomMedium::KInverse k_inverse; +// const SingleCurvingCrack::KInverse k_inverse; - typename DoFHandler::active_cell_iterator - cell = saturation_dof_handler.begin_active(), - endc = saturation_dof_handler.end(); - for (; cell!=endc; ++cell) - { - saturation_fe_values.reinit (cell); - local_matrix = 0; - local_rhs = 0; - - for (unsigned int q=0; qget_dof_indices (local_dof_indices); - - saturation_constraints.distribute_local_to_global (local_matrix, - local_dof_indices, - saturation_matrix); + std::vector pressure_rhs_values (n_q_points); + std::vector boundary_values (n_face_q_points); + std::vector > k_inverse_values (n_q_points); - } -} + std::vector old_saturation_values (n_q_points); + std::vector > phi_u (dofs_per_cell); + std::vector div_phi_u (dofs_per_cell); + std::vector phi_p (dofs_per_cell); + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); - // @sect3{TwoPhaseFlowProblem::assemble_saturation_rhs} - - // This function is to assemble the right - // hand side of the saturation transport - // equation. Before assembling it, we have to - // call two FEValues objects for the darcy - // and saturation systems respectively and, - // even more, two FEFaceValues objects for - // the both systems because we have a - // boundary integral term in the weak form of - // saturation equation. For the FEFaceValues - // object of the saturation system, we also - // enter the normal vectors with an update - // flag update_normal_vectors. - // - // Next, before looping over all the cells, - // we have to compute some parameters - // (e.g. global_u_infty, global_S_variasion, - // and global_Omega_diameter) that the - // artificial viscosity $\nu$ needs, which - // desriptions have been appearing in - // step-31. - // - // Next, we start to loop over all the - // saturation and darcy cells to put the - // local contributions into the global - // vector. In this loop, in order to simplify - // the implementation in this function, we - // generate two more functions: one is - // assemble_saturation_rhs_cell_term and the - // other is - // assemble_saturation_rhs_boundary_term, - // which is contained in an inner boudary - // loop. The former is to assemble the - // integral cell term with neccessary - // arguments and the latter is to assemble - // the integral global boundary $\Omega$ - // terms. It should be noted that we achieve - // the insertion of the cell or boundary - // vector elements to the global vector in - // the two functions rather than in this - // present function by giving these two - // functions with a common argument - // local_dof_indices, and two arguments - // saturation_fe_values darcy_fe_values for - // assemble_saturation_rhs_cell_term and - // another two arguments - // saturation_fe_face_values - // darcy_fe_face_values for - // assemble_saturation_rhs_boundary_term. -template -void TwoPhaseFlowProblem::assemble_saturation_rhs () -{ - QGauss quadrature_formula(saturation_degree+2); - QGauss face_quadrature_formula(saturation_degree+2); - - FEValues saturation_fe_values (saturation_fe, quadrature_formula, - update_values | update_gradients | - update_quadrature_points | update_JxW_values); - FEValues darcy_fe_values (darcy_fe, quadrature_formula, - update_values); - FEFaceValues saturation_fe_face_values (saturation_fe, face_quadrature_formula, - update_values | update_normal_vectors | - update_quadrature_points | update_JxW_values); - FEFaceValues darcy_fe_face_values (darcy_fe, face_quadrature_formula, - update_values); - FEFaceValues saturation_fe_face_values_neighbor (saturation_fe, face_quadrature_formula, - update_values); - - const unsigned int dofs_per_cell = saturation_dof_handler.get_fe().dofs_per_cell; - std::vector local_dof_indices (dofs_per_cell); - - const double global_u_infty = get_maximal_velocity (); - const std::pair - global_S_range = get_extrapolated_saturation_range (); - const double global_S_variasion = global_S_range.second - global_S_range.first; - const double global_Omega_diameter = GridTools::diameter (triangulation); - - typename DoFHandler::active_cell_iterator - cell = saturation_dof_handler.begin_active(), - endc = saturation_dof_handler.end(); - typename DoFHandler::active_cell_iterator - darcy_cell = darcy_dof_handler.begin_active(); - for (; cell!=endc; ++cell, ++darcy_cell) - { - saturation_fe_values.reinit (cell); - darcy_fe_values.reinit (darcy_cell); - - cell->get_dof_indices (local_dof_indices); - - assemble_saturation_rhs_cell_term(saturation_fe_values, - darcy_fe_values, - local_dof_indices, - global_u_infty, - global_S_variasion, - global_Omega_diameter); - - for (unsigned int face_no=0; face_no::faces_per_cell; - ++face_no) - { - - if (cell->at_boundary(face_no)) - { - darcy_fe_face_values.reinit (darcy_cell, face_no); - saturation_fe_face_values.reinit (cell, face_no); - assemble_saturation_rhs_boundary_term (saturation_fe_face_values, - darcy_fe_face_values, - local_dof_indices); - } - } - } -} + typename DoFHandler::active_cell_iterator + cell = darcy_dof_handler.begin_active(), + endc = darcy_dof_handler.end(); + typename DoFHandler::active_cell_iterator + saturation_cell = saturation_dof_handler.begin_active(); + for (; cell!=endc; ++cell, ++saturation_cell) + { + darcy_fe_values.reinit (cell); + saturation_fe_values.reinit (saturation_cell); + + local_matrix = 0; + local_rhs = 0; + + saturation_fe_values.get_function_values (old_saturation_solution, old_saturation_values); + + pressure_right_hand_side.value_list (darcy_fe_values.get_quadrature_points(), + pressure_rhs_values); + k_inverse.value_list (darcy_fe_values.get_quadrature_points(), + k_inverse_values); + + for (unsigned int q=0; q::faces_per_cell; + ++face_no) + if (cell->at_boundary(face_no)) + { + darcy_fe_face_values.reinit (cell, face_no); + + pressure_boundary_values + .value_list (darcy_fe_face_values.get_quadrature_points(), + boundary_values); + + for (unsigned int q=0; q + phi_i_u = darcy_fe_face_values[velocities].value (i, q); + + local_rhs(i) += -(phi_i_u * + darcy_fe_face_values.normal_vector(q) * + boundary_values[q] * + darcy_fe_face_values.JxW(q)); + } + } + + for (unsigned int i=0; iget_dof_indices (local_dof_indices); + + darcy_constraints.distribute_local_to_global (local_matrix, + local_rhs, + local_dof_indices, + darcy_matrix, + darcy_rhs); + } + } - // @sect3{TwoPhaseFlowProblem::assemble_saturation_rhs_cell_term} - - // In this function, we actually compute - // every artificial viscosity for every - // element. Then, with the artificial value, - // we can finish assembling the saturation - // right hand side cell integral - // terms. Finally, we can pass the local - // contributions on to the global vector with - // the position specified in - // local_dof_indices. -template -void -TwoPhaseFlowProblem:: -assemble_saturation_rhs_cell_term (const FEValues &saturation_fe_values, - const FEValues &darcy_fe_values, - const std::vector &local_dof_indices, - const double global_u_infty, - const double global_S_variation, - const double global_Omega_diameter) -{ - const unsigned int dofs_per_cell = saturation_fe_values.dofs_per_cell; - const unsigned int n_q_points = saturation_fe_values.n_quadrature_points; - - Vector local_rhs (dofs_per_cell); - - std::vector old_saturation_solution_values(n_q_points); - std::vector old_old_saturation_solution_values(n_q_points); - std::vector > old_grad_saturation_solution_values(n_q_points); - std::vector > old_old_grad_saturation_solution_values(n_q_points); - std::vector > present_darcy_solution_values(n_q_points, Vector(dim+1)); - - saturation_fe_values.get_function_values (old_saturation_solution, old_saturation_solution_values); - saturation_fe_values.get_function_values (old_old_saturation_solution, old_old_saturation_solution_values); - saturation_fe_values.get_function_grads (old_saturation_solution, old_grad_saturation_solution_values); - saturation_fe_values.get_function_grads (old_old_saturation_solution, old_old_grad_saturation_solution_values); - darcy_fe_values.get_function_values (darcy_solution, present_darcy_solution_values); - - const double nu - = compute_viscosity (old_saturation_solution_values, - old_old_saturation_solution_values, - old_grad_saturation_solution_values, - old_old_grad_saturation_solution_values, - present_darcy_solution_values, - global_u_infty, - global_S_variation, - global_Omega_diameter, - saturation_fe_values.get_cell()->diameter(), - old_time_step, - viscosity); - - for (unsigned int q=0; q::assemble_saturation_system} + + // This function is to assemble the linear + // system for the saturation transport + // equation. It includes two member + // functions: assemble_saturation_matrix () + // and assemble_saturation_rhs (). The former + // function that assembles the saturation + // left hand side needs to be changed only + // when grids have been changed since the + // matrix is filled only with basis + // functions. However, the latter that + // assembles the right hand side must be + // changed at every saturation time step + // since it depends on an unknown variable + // saturation. + template + void TwoPhaseFlowProblem::assemble_saturation_system () + { + if ( rebuild_saturation_matrix == true ) { - const double old_s = old_saturation_solution_values[q]; - Tensor<1,dim> present_u; - for (unsigned int d=0; d grad_phi_i_s = saturation_fe_values.shape_grad (i, q); - - local_rhs(i) += (time_step * - f_saturation(old_s,viscosity) * - present_u * - grad_phi_i_s - - - time_step * - nu * - old_grad_saturation_solution_values[q] * grad_phi_i_s - + - old_s * phi_i_s) - * - saturation_fe_values.JxW(q); + saturation_matrix = 0; + assemble_saturation_matrix (); } - saturation_constraints.distribute_local_to_global (local_rhs, - local_dof_indices, - saturation_rhs); -} + saturation_rhs = 0; + assemble_saturation_rhs (); + } - // @sect3{TwoPhaseFlowProblem::assemble_saturation_rhs_boundary_term} - - // In this function, we have to give - // upwinding in the global boundary faces, - // i.e. we impose the Dirichlet boundary - // conditions only on inflow parts of global - // boundary, which has been described in - // step-21 so we refrain from giving more - // descriptions about that. -template -void -TwoPhaseFlowProblem:: -assemble_saturation_rhs_boundary_term (const FEFaceValues &saturation_fe_face_values, - const FEFaceValues &darcy_fe_face_values, - const std::vector &local_dof_indices) -{ - const unsigned int dofs_per_cell = saturation_fe_face_values.dofs_per_cell; - const unsigned int n_face_q_points = saturation_fe_face_values.n_quadrature_points; - Vector local_rhs (dofs_per_cell); + // @sect3{TwoPhaseFlowProblem::assemble_saturation_matrix} - std::vector old_saturation_solution_values_face(n_face_q_points); - std::vector > present_darcy_solution_values_face(n_face_q_points, Vector(dim+1)); - std::vector neighbor_saturation (n_face_q_points); + // This function is easily understood since + // it only forms a simple mass matrix for the + // left hand side of the saturation linear + // system by basis functions phi_i_s and + // phi_j_s only. Finally, as usual, we enter + // the local contribution into the global + // matrix by specifying the position in + // local_dof_indices. This is done by letting + // the ConstraintMatrix class do the + // insertion of the cell matrix elements to + // the global matrix, which already condenses + // the hanging node constraints. + template + void TwoPhaseFlowProblem::assemble_saturation_matrix () + { + QGauss quadrature_formula(saturation_degree+2); - saturation_fe_face_values.get_function_values (old_saturation_solution, old_saturation_solution_values_face); - darcy_fe_face_values.get_function_values (darcy_solution, present_darcy_solution_values_face); + FEValues saturation_fe_values (saturation_fe, quadrature_formula, + update_values | update_JxW_values); - SaturationBoundaryValues saturation_boundary_values; - saturation_boundary_values - .value_list (saturation_fe_face_values.get_quadrature_points(), - neighbor_saturation); + const unsigned int dofs_per_cell = saturation_fe.dofs_per_cell; - for (unsigned int q=0; q present_u_face; - for (unsigned int d=0; d local_matrix (dofs_per_cell, dofs_per_cell); + Vector local_rhs (dofs_per_cell); - const bool is_outflow_q_point = (normal_flux >= 0); + std::vector local_dof_indices (dofs_per_cell); - for (unsigned int i=0; i::active_cell_iterator + cell = saturation_dof_handler.begin_active(), + endc = saturation_dof_handler.end(); + for (; cell!=endc; ++cell) + { + saturation_fe_values.reinit (cell); + local_matrix = 0; + local_rhs = 0; + + for (unsigned int q=0; qget_dof_indices (local_dof_indices); + + saturation_constraints.distribute_local_to_global (local_matrix, + local_dof_indices, + saturation_matrix); + } + } - // @sect3{TwoPhaseFlowProblem::solve} - - // This function is to implement the operator - // splitting algorithm. At the beginning of - // the implementation, we decide whther to - // solve the pressure-velocity part by - // running an a posteriori criterion, which - // will be described in the following - // function. If we get the bool variable true - // from that function, we will solve the - // pressure-velocity part for updated - // velocity. Then, we use GMRES with the - // Schur complement preconditioner to solve - // this linear system, as is described in the - // Introduction. After solving the velocity - // and pressure, we need to keep the - // solutions for linear extrapolations in the - // future. It is noted that we always solve - // the pressure-velocity part in the first - // three micro time steps to ensure accuracy - // at the beginning of computation, and to - // provide starting data to linearly - // extrapolate previously computed velocities - // to the current time step. - // - // On the other hand, if we get a false - // variable from the criterion, we will - // directly use linear extrapolation to - // compute the updated velocity for the - // solution of saturation later. - // - // Next, like step-21, this program need to - // compute the present time step. - // - // Next, we need to use two bool variables - // solve_pressure_velocity_part and - // previous_solve_pressure_velocity_part to - // decide whether we stop or continue - // cumulating the micro time steps for linear - // extropolations in the next iteration. With - // the reason, we need one variable - // cumulative_nth_time_step for keeping the - // present aggregated micro time steps and - // anther one n_minus_oneth_time_step for - // retaining the previous micro time steps. - // - // Finally, we start to calculate the - // saturation part with the use of the - // incomplete Cholesky decomposition for - // preconditioning. -template -void TwoPhaseFlowProblem::solve () -{ - solve_pressure_velocity_part = determine_whether_to_solve_pressure_velocity_part (); - if ( timestep_number <= 3 || solve_pressure_velocity_part == true ) - { - std::cout << " Solving darcy system (pressure-velocity part)..." << std::endl; - assemble_darcy_system (); - build_darcy_preconditioner (); + // @sect3{TwoPhaseFlowProblem::assemble_saturation_rhs} + + // This function is to assemble the right + // hand side of the saturation transport + // equation. Before assembling it, we have to + // call two FEValues objects for the darcy + // and saturation systems respectively and, + // even more, two FEFaceValues objects for + // the both systems because we have a + // boundary integral term in the weak form of + // saturation equation. For the FEFaceValues + // object of the saturation system, we also + // enter the normal vectors with an update + // flag update_normal_vectors. + // + // Next, before looping over all the cells, + // we have to compute some parameters + // (e.g. global_u_infty, global_S_variasion, + // and global_Omega_diameter) that the + // artificial viscosity $\nu$ needs, which + // desriptions have been appearing in + // step-31. + // + // Next, we start to loop over all the + // saturation and darcy cells to put the + // local contributions into the global + // vector. In this loop, in order to simplify + // the implementation in this function, we + // generate two more functions: one is + // assemble_saturation_rhs_cell_term and the + // other is + // assemble_saturation_rhs_boundary_term, + // which is contained in an inner boudary + // loop. The former is to assemble the + // integral cell term with neccessary + // arguments and the latter is to assemble + // the integral global boundary $\Omega$ + // terms. It should be noted that we achieve + // the insertion of the cell or boundary + // vector elements to the global vector in + // the two functions rather than in this + // present function by giving these two + // functions with a common argument + // local_dof_indices, and two arguments + // saturation_fe_values darcy_fe_values for + // assemble_saturation_rhs_cell_term and + // another two arguments + // saturation_fe_face_values + // darcy_fe_face_values for + // assemble_saturation_rhs_boundary_term. + template + void TwoPhaseFlowProblem::assemble_saturation_rhs () + { + QGauss quadrature_formula(saturation_degree+2); + QGauss face_quadrature_formula(saturation_degree+2); + + FEValues saturation_fe_values (saturation_fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + FEValues darcy_fe_values (darcy_fe, quadrature_formula, + update_values); + FEFaceValues saturation_fe_face_values (saturation_fe, face_quadrature_formula, + update_values | update_normal_vectors | + update_quadrature_points | update_JxW_values); + FEFaceValues darcy_fe_face_values (darcy_fe, face_quadrature_formula, + update_values); + FEFaceValues saturation_fe_face_values_neighbor (saturation_fe, face_quadrature_formula, + update_values); + + const unsigned int dofs_per_cell = saturation_dof_handler.get_fe().dofs_per_cell; + std::vector local_dof_indices (dofs_per_cell); + + const double global_u_infty = get_maximal_velocity (); + const std::pair + global_S_range = get_extrapolated_saturation_range (); + const double global_S_variasion = global_S_range.second - global_S_range.first; + const double global_Omega_diameter = GridTools::diameter (triangulation); + typename DoFHandler::active_cell_iterator + cell = saturation_dof_handler.begin_active(), + endc = saturation_dof_handler.end(); + typename DoFHandler::active_cell_iterator + darcy_cell = darcy_dof_handler.begin_active(); + for (; cell!=endc; ++cell, ++darcy_cell) { - const LinearSolvers::InverseMatrix - mp_inverse (darcy_preconditioner_matrix.block(1,1), *Mp_preconditioner); + saturation_fe_values.reinit (cell); + darcy_fe_values.reinit (darcy_cell); + + cell->get_dof_indices (local_dof_indices); + + assemble_saturation_rhs_cell_term(saturation_fe_values, + darcy_fe_values, + local_dof_indices, + global_u_infty, + global_S_variasion, + global_Omega_diameter); + + for (unsigned int face_no=0; face_no::faces_per_cell; + ++face_no) + { + + if (cell->at_boundary(face_no)) + { + darcy_fe_face_values.reinit (darcy_cell, face_no); + saturation_fe_face_values.reinit (cell, face_no); + assemble_saturation_rhs_boundary_term (saturation_fe_face_values, + darcy_fe_face_values, + local_dof_indices); + } + } + } + } - const LinearSolvers::BlockSchurPreconditioner - preconditioner (darcy_matrix, mp_inverse, *Amg_preconditioner); - SolverControl solver_control (darcy_matrix.m(), - 1e-6*darcy_rhs.l2_norm()); - SolverGMRES - gmres (solver_control, - SolverGMRES::AdditionalData(100)); + // @sect3{TwoPhaseFlowProblem::assemble_saturation_rhs_cell_term} - for (unsigned int i=0; i + void + TwoPhaseFlowProblem:: + assemble_saturation_rhs_cell_term (const FEValues &saturation_fe_values, + const FEValues &darcy_fe_values, + const std::vector &local_dof_indices, + const double global_u_infty, + const double global_S_variation, + const double global_Omega_diameter) + { + const unsigned int dofs_per_cell = saturation_fe_values.dofs_per_cell; + const unsigned int n_q_points = saturation_fe_values.n_quadrature_points; + + Vector local_rhs (dofs_per_cell); + + std::vector old_saturation_solution_values(n_q_points); + std::vector old_old_saturation_solution_values(n_q_points); + std::vector > old_grad_saturation_solution_values(n_q_points); + std::vector > old_old_grad_saturation_solution_values(n_q_points); + std::vector > present_darcy_solution_values(n_q_points, Vector(dim+1)); + + saturation_fe_values.get_function_values (old_saturation_solution, old_saturation_solution_values); + saturation_fe_values.get_function_values (old_old_saturation_solution, old_old_saturation_solution_values); + saturation_fe_values.get_function_grads (old_saturation_solution, old_grad_saturation_solution_values); + saturation_fe_values.get_function_grads (old_old_saturation_solution, old_old_grad_saturation_solution_values); + darcy_fe_values.get_function_values (darcy_solution, present_darcy_solution_values); + + const double nu + = compute_viscosity (old_saturation_solution_values, + old_old_saturation_solution_values, + old_grad_saturation_solution_values, + old_old_grad_saturation_solution_values, + present_darcy_solution_values, + global_u_infty, + global_S_variation, + global_Omega_diameter, + saturation_fe_values.get_cell()->diameter(), + old_time_step, + viscosity); + + for (unsigned int q=0; q present_u; + for (unsigned int d=0; d grad_phi_i_s = saturation_fe_values.shape_grad (i, q); + + local_rhs(i) += (time_step * + f_saturation(old_s,viscosity) * + present_u * + grad_phi_i_s + - + time_step * + nu * + old_grad_saturation_solution_values[q] * grad_phi_i_s + + + old_s * phi_i_s) + * + saturation_fe_values.JxW(q); + } - gmres.solve(darcy_matrix, darcy_solution, darcy_rhs, preconditioner); + saturation_constraints.distribute_local_to_global (local_rhs, + local_dof_indices, + saturation_rhs); + } - darcy_constraints.distribute (darcy_solution); - std::cout << " " - << solver_control.last_step() - << " GMRES iterations for darcy system (pressure-velocity part)." - << std::endl; + // @sect3{TwoPhaseFlowProblem::assemble_saturation_rhs_boundary_term} - } + // In this function, we have to give + // upwinding in the global boundary faces, + // i.e. we impose the Dirichlet boundary + // conditions only on inflow parts of global + // boundary, which has been described in + // step-21 so we refrain from giving more + // descriptions about that. + template + void + TwoPhaseFlowProblem:: + assemble_saturation_rhs_boundary_term (const FEFaceValues &saturation_fe_face_values, + const FEFaceValues &darcy_fe_face_values, + const std::vector &local_dof_indices) + { + const unsigned int dofs_per_cell = saturation_fe_face_values.dofs_per_cell; + const unsigned int n_face_q_points = saturation_fe_face_values.n_quadrature_points; - { - n_minus_oneth_darcy_solution_after_solving_pressure_part = nth_darcy_solution_after_solving_pressure_part; - nth_darcy_solution_after_solving_pressure_part = darcy_solution; + Vector local_rhs (dofs_per_cell); - nth_saturation_solution_after_solving_pressure_part = saturation_solution; - } - } - else - { - darcy_solution.block(0) = nth_darcy_solution_after_solving_pressure_part.block(0); - darcy_solution.block(0).sadd (2.0, -1.0, n_minus_oneth_darcy_solution_after_solving_pressure_part.block(0) ); + std::vector old_saturation_solution_values_face(n_face_q_points); + std::vector > present_darcy_solution_values_face(n_face_q_points, Vector(dim+1)); + std::vector neighbor_saturation (n_face_q_points); - double extrapolated_time_step = GridTools::minimal_cell_diameter(triangulation) / - get_maximal_velocity() / 8.0; + saturation_fe_face_values.get_function_values (old_saturation_solution, old_saturation_solution_values_face); + darcy_fe_face_values.get_function_values (darcy_solution, present_darcy_solution_values_face); - double local_cumulative_time_step = cumulative_nth_time_step + extrapolated_time_step; - double coef_1 = local_cumulative_time_step / n_minus_oneth_time_step; - double coef_2 = ( 1.0 + coef_1 ); + SaturationBoundaryValues saturation_boundary_values; + saturation_boundary_values + .value_list (saturation_fe_face_values.get_quadrature_points(), + neighbor_saturation); - TrilinosWrappers::Vector tmp (darcy_solution.block(0).size()); - tmp = nth_darcy_solution_after_solving_pressure_part.block(0); + for (unsigned int q=0; q present_u_face; + for (unsigned int d=0; d= 0); + + for (unsigned int i=0; i::solve} + + // This function is to implement the operator + // splitting algorithm. At the beginning of + // the implementation, we decide whther to + // solve the pressure-velocity part by + // running an a posteriori criterion, which + // will be described in the following + // function. If we get the bool variable true + // from that function, we will solve the + // pressure-velocity part for updated + // velocity. Then, we use GMRES with the + // Schur complement preconditioner to solve + // this linear system, as is described in the + // Introduction. After solving the velocity + // and pressure, we need to keep the + // solutions for linear extrapolations in the + // future. It is noted that we always solve + // the pressure-velocity part in the first + // three micro time steps to ensure accuracy + // at the beginning of computation, and to + // provide starting data to linearly + // extrapolate previously computed velocities + // to the current time step. + // + // On the other hand, if we get a false + // variable from the criterion, we will + // directly use linear extrapolation to + // compute the updated velocity for the + // solution of saturation later. + // + // Next, like step-21, this program need to + // compute the present time step. + // + // Next, we need to use two bool variables + // solve_pressure_velocity_part and + // previous_solve_pressure_velocity_part to + // decide whether we stop or continue + // cumulating the micro time steps for linear + // extropolations in the next iteration. With + // the reason, we need one variable + // cumulative_nth_time_step for keeping the + // present aggregated micro time steps and + // anther one n_minus_oneth_time_step for + // retaining the previous micro time steps. + // + // Finally, we start to calculate the + // saturation part with the use of the + // incomplete Cholesky decomposition for + // preconditioning. + template + void TwoPhaseFlowProblem::solve () + { + solve_pressure_velocity_part = determine_whether_to_solve_pressure_velocity_part (); + if ( timestep_number <= 3 || solve_pressure_velocity_part == true ) + { + std::cout << " Solving darcy system (pressure-velocity part)..." << std::endl; - old_time_step = time_step; - time_step = GridTools::minimal_cell_diameter(triangulation) / - get_maximal_velocity() / 8.0; + assemble_darcy_system (); + build_darcy_preconditioner (); - if ( timestep_number <= 3 || ( solve_pressure_velocity_part == true && previous_solve_pressure_velocity_part == true ) ) - { - n_minus_oneth_time_step = time_step; - cumulative_nth_time_step = 0.0; - } - else if ( solve_pressure_velocity_part == true && previous_solve_pressure_velocity_part == false ) - { - n_minus_oneth_time_step = cumulative_nth_time_step; - cumulative_nth_time_step = 0.0; - } - else - { - cumulative_nth_time_step += time_step; - } + { + const LinearSolvers::InverseMatrix + mp_inverse (darcy_preconditioner_matrix.block(1,1), *Mp_preconditioner); - previous_solve_pressure_velocity_part = solve_pressure_velocity_part; + const LinearSolvers::BlockSchurPreconditioner + preconditioner (darcy_matrix, mp_inverse, *Amg_preconditioner); - std::cout << " Solving saturation transport equation..." << std::endl; + SolverControl solver_control (darcy_matrix.m(), + 1e-6*darcy_rhs.l2_norm()); - assemble_saturation_system (); + SolverGMRES + gmres (solver_control, + SolverGMRES::AdditionalData(100)); - { - SolverControl solver_control (saturation_matrix.m(), - 1e-8*saturation_rhs.l2_norm()); - SolverCG cg (solver_control); + for (unsigned int i=0; i::determine_whether_to_solve_pressure_velocity_part} - - // This function is to implement the a - // posteriori criterion for - // adaptive operator splitting. As mentioned - // in step-31, we use two FEValues objects - // initialized with two cell iterators that - // we walk in parallel through the two - // DoFHandler objects associated with the - // same Triangulation object; for these two - // FEValues objects, we use of course the - // same quadrature objects so that we can - // iterate over the same set of quadrature - // points, but each FEValues object will get - // update flags only according to what it - // actually needs to compute. - // - // In addition to this, if someone doesn't - // want to perform their simulation with - // operator splitting, they can lower the - // criterion value (default value is $5.0$) - // down to zero ad therefore numerical - // algorithm becomes the original IMPES - // method. -template -bool -TwoPhaseFlowProblem::determine_whether_to_solve_pressure_velocity_part () const -{ - if (timestep_number <= 3) - return true; + darcy_solution.block(0).sadd (0.5, 0.5, tmp); + } - const QGauss quadrature_formula(saturation_degree+2); - const unsigned int n_q_points = quadrature_formula.size(); - FEValues fe_values (saturation_fe, quadrature_formula, - update_values | update_quadrature_points); + old_time_step = time_step; + time_step = GridTools::minimal_cell_diameter(triangulation) / + get_maximal_velocity() / 8.0; - std::vector old_saturation_after_solving_pressure (n_q_points); - std::vector present_saturation (n_q_points); + if ( timestep_number <= 3 || ( solve_pressure_velocity_part == true && previous_solve_pressure_velocity_part == true ) ) + { + n_minus_oneth_time_step = time_step; + cumulative_nth_time_step = 0.0; + } + else if ( solve_pressure_velocity_part == true && previous_solve_pressure_velocity_part == false ) + { + n_minus_oneth_time_step = cumulative_nth_time_step; + cumulative_nth_time_step = 0.0; + } + else + { + cumulative_nth_time_step += time_step; + } - const RandomMedium::KInverse k_inverse; -// const SingleCurvingCrack::KInverse k_inverse; + previous_solve_pressure_velocity_part = solve_pressure_velocity_part; - std::vector > k_inverse_values (n_q_points); + std::cout << " Solving saturation transport equation..." << std::endl; - double max_global_aop_indicator = 0.0; + assemble_saturation_system (); - typename DoFHandler::active_cell_iterator - cell = saturation_dof_handler.begin_active(), - endc = saturation_dof_handler.end(); - for (; cell!=endc; ++cell) { - double max_local_mobility_reciprocal_difference = 0.0; - double max_local_permeability_inverse_l1_norm = 0.0; + SolverControl solver_control (saturation_matrix.m(), + 1e-8*saturation_rhs.l2_norm()); + SolverCG cg (solver_control); - fe_values.reinit(cell); - fe_values.get_function_values (nth_saturation_solution_after_solving_pressure_part, - old_saturation_after_solving_pressure); - fe_values.get_function_values (saturation_solution, - present_saturation); + TrilinosWrappers::PreconditionIC preconditioner; + preconditioner.initialize (saturation_matrix); - k_inverse.value_list (fe_values.get_quadrature_points(), - k_inverse_values); + cg.solve (saturation_matrix, saturation_solution, + saturation_rhs, preconditioner); - for (unsigned int q=0; q 5.0 ) - { - return true; } - else - { - std::cout << " Activating adaptive operating splitting" << std::endl; - return false; - } -} - + } - // @sect3{TwoPhaseFlowProblem::compute_refinement_indicators} - - // This function is to to compute the - // refinement indicator discussed in the - // introduction for each cell and its - // implementation is similar to that - // contained in step-33. There is no need to - // repeat descriptions about it. -template -void -TwoPhaseFlowProblem:: -compute_refinement_indicators (Vector &refinement_indicators) const -{ - const QMidpoint quadrature_formula; - FEValues fe_values (saturation_fe, quadrature_formula, update_gradients); - std::vector > grad_saturation (1); - double max_refinement_indicator = 0.0; + // @sect3{TwoPhaseFlowProblem::determine_whether_to_solve_pressure_velocity_part} + + // This function is to implement the a + // posteriori criterion for + // adaptive operator splitting. As mentioned + // in step-31, we use two FEValues objects + // initialized with two cell iterators that + // we walk in parallel through the two + // DoFHandler objects associated with the + // same Triangulation object; for these two + // FEValues objects, we use of course the + // same quadrature objects so that we can + // iterate over the same set of quadrature + // points, but each FEValues object will get + // update flags only according to what it + // actually needs to compute. + // + // In addition to this, if someone doesn't + // want to perform their simulation with + // operator splitting, they can lower the + // criterion value (default value is $5.0$) + // down to zero ad therefore numerical + // algorithm becomes the original IMPES + // method. + template + bool + TwoPhaseFlowProblem::determine_whether_to_solve_pressure_velocity_part () const + { + if (timestep_number <= 3) + return true; - typename DoFHandler::active_cell_iterator - cell = saturation_dof_handler.begin_active(), - endc = saturation_dof_handler.end(); - for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) - { - fe_values.reinit(cell); - fe_values.get_function_grads (predictor_saturation_solution, - grad_saturation); - - refinement_indicators(cell_no) - = std::log( 1.0 + std::sqrt( grad_saturation[0] * - grad_saturation[0] ) ); - max_refinement_indicator = std::max(max_refinement_indicator, - refinement_indicators(cell_no)); - } + const QGauss quadrature_formula(saturation_degree+2); + const unsigned int n_q_points = quadrature_formula.size(); -// std::cout << "max_refinement_indicator =" << max_refinement_indicator << std::endl; -} + FEValues fe_values (saturation_fe, quadrature_formula, + update_values | update_quadrature_points); + std::vector old_saturation_after_solving_pressure (n_q_points); + std::vector present_saturation (n_q_points); + const RandomMedium::KInverse k_inverse; +// const SingleCurvingCrack::KInverse k_inverse; - // @sect3{TwoPhaseFlowProblem::refine_grid} + std::vector > k_inverse_values (n_q_points); - // This function is to decide if every cell - // is refined or coarsened with computed - // refinement indicators in the previous - // function and do the interpolations of the - // solution vectors. The main difference from - // the previous time-dependent tutorials is - // that there is no need to do the solution - // interpolations if we don't have any cell - // that is refined or coarsend, saving some - // additional computing time. -template -void -TwoPhaseFlowProblem:: -refine_grid (const Vector &refinement_indicators) -{ - const double current_saturation_level = saturation_level + - n_refinement_steps; + double max_global_aop_indicator = 0.0; - { typename DoFHandler::active_cell_iterator cell = saturation_dof_handler.begin_active(), endc = saturation_dof_handler.end(); - - for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) + for (; cell!=endc; ++cell) { - cell->clear_coarsen_flag(); - cell->clear_refine_flag(); + double max_local_mobility_reciprocal_difference = 0.0; + double max_local_permeability_inverse_l1_norm = 0.0; - if ((cell->level() < current_saturation_level) && - (std::fabs(refinement_indicators(cell_no)) > saturation_value)) - cell->set_refine_flag(); - else - if ((cell->level() > double(n_refinement_steps)) && - (std::fabs(refinement_indicators(cell_no)) < 0.75 * saturation_value)) - cell->set_coarsen_flag(); + fe_values.reinit(cell); + fe_values.get_function_values (nth_saturation_solution_after_solving_pressure_part, + old_saturation_after_solving_pressure); + fe_values.get_function_values (saturation_solution, + present_saturation); + + k_inverse.value_list (fe_values.get_quadrature_points(), + k_inverse_values); + + for (unsigned int q=0; q 5.0 ) + { + return true; + } + else + { + std::cout << " Activating adaptive operating splitting" << std::endl; + return false; } } - triangulation.prepare_coarsening_and_refinement (); - unsigned int number_of_cells_refine = 0; - unsigned int number_of_cells_coarsen = 0; + // @sect3{TwoPhaseFlowProblem::compute_refinement_indicators} + + // This function is to to compute the + // refinement indicator discussed in the + // introduction for each cell and its + // implementation is similar to that + // contained in step-33. There is no need to + // repeat descriptions about it. + template + void + TwoPhaseFlowProblem:: + compute_refinement_indicators (Vector &refinement_indicators) const { + + const QMidpoint quadrature_formula; + FEValues fe_values (saturation_fe, quadrature_formula, update_gradients); + std::vector > grad_saturation (1); + + double max_refinement_indicator = 0.0; + typename DoFHandler::active_cell_iterator cell = saturation_dof_handler.begin_active(), endc = saturation_dof_handler.end(); + for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) + { + fe_values.reinit(cell); + fe_values.get_function_grads (predictor_saturation_solution, + grad_saturation); + + refinement_indicators(cell_no) + = std::log( 1.0 + std::sqrt( grad_saturation[0] * + grad_saturation[0] ) ); + max_refinement_indicator = std::max(max_refinement_indicator, + refinement_indicators(cell_no)); + } - for (; cell!=endc; ++cell) - if (cell->refine_flag_set()) - ++number_of_cells_refine; - else - if (cell->coarsen_flag_set()) - ++number_of_cells_coarsen; +// std::cout << "max_refinement_indicator =" << max_refinement_indicator << std::endl; } - std::cout << " " - << number_of_cells_refine - << " cell(s) are going to be refined." - << std::endl; - std::cout << " " - << number_of_cells_coarsen - << " cell(s) are going to be coarsened." - << std::endl; - std::cout << std::endl; - if ( number_of_cells_refine > 0 || number_of_cells_coarsen > 0 ) - { - std::vector x_saturation (3); - x_saturation[0] = saturation_solution; - x_saturation[1] = old_saturation_solution; - x_saturation[2] = nth_saturation_solution_after_solving_pressure_part; - - std::vector x_darcy (2); - x_darcy[0] = nth_darcy_solution_after_solving_pressure_part; - x_darcy[1] = n_minus_oneth_darcy_solution_after_solving_pressure_part; + // @sect3{TwoPhaseFlowProblem::refine_grid} - SolutionTransfer saturation_soltrans(saturation_dof_handler); + // This function is to decide if every cell + // is refined or coarsened with computed + // refinement indicators in the previous + // function and do the interpolations of the + // solution vectors. The main difference from + // the previous time-dependent tutorials is + // that there is no need to do the solution + // interpolations if we don't have any cell + // that is refined or coarsend, saving some + // additional computing time. + template + void + TwoPhaseFlowProblem:: + refine_grid (const Vector &refinement_indicators) + { + const double current_saturation_level = saturation_level + + n_refinement_steps; - SolutionTransfer darcy_soltrans(darcy_dof_handler); + { + typename DoFHandler::active_cell_iterator + cell = saturation_dof_handler.begin_active(), + endc = saturation_dof_handler.end(); + + for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) + { + cell->clear_coarsen_flag(); + cell->clear_refine_flag(); + + if ((cell->level() < current_saturation_level) && + (std::fabs(refinement_indicators(cell_no)) > saturation_value)) + cell->set_refine_flag(); + else + if ((cell->level() > double(n_refinement_steps)) && + (std::fabs(refinement_indicators(cell_no)) < 0.75 * saturation_value)) + cell->set_coarsen_flag(); + } + } + triangulation.prepare_coarsening_and_refinement (); - triangulation.prepare_coarsening_and_refinement(); - saturation_soltrans.prepare_for_coarsening_and_refinement(x_saturation); + unsigned int number_of_cells_refine = 0; + unsigned int number_of_cells_coarsen = 0; - darcy_soltrans.prepare_for_coarsening_and_refinement(x_darcy); + { + typename DoFHandler::active_cell_iterator + cell = saturation_dof_handler.begin_active(), + endc = saturation_dof_handler.end(); - triangulation.execute_coarsening_and_refinement (); - setup_dofs (); + for (; cell!=endc; ++cell) + if (cell->refine_flag_set()) + ++number_of_cells_refine; + else + if (cell->coarsen_flag_set()) + ++number_of_cells_coarsen; + } - std::vector tmp_saturation (3); - tmp_saturation[0].reinit (saturation_solution); - tmp_saturation[1].reinit (saturation_solution); - tmp_saturation[2].reinit (saturation_solution); - saturation_soltrans.interpolate(x_saturation, tmp_saturation); + std::cout << " " + << number_of_cells_refine + << " cell(s) are going to be refined." + << std::endl; + std::cout << " " + << number_of_cells_coarsen + << " cell(s) are going to be coarsened." + << std::endl; - saturation_solution = tmp_saturation[0]; - old_saturation_solution = tmp_saturation[1]; - nth_saturation_solution_after_solving_pressure_part = tmp_saturation[2]; + std::cout << std::endl; - std::vector tmp_darcy (2); - tmp_darcy[0].reinit (darcy_solution); - tmp_darcy[1].reinit (darcy_solution); - darcy_soltrans.interpolate(x_darcy, tmp_darcy); + if ( number_of_cells_refine > 0 || number_of_cells_coarsen > 0 ) + { + std::vector x_saturation (3); + x_saturation[0] = saturation_solution; + x_saturation[1] = old_saturation_solution; + x_saturation[2] = nth_saturation_solution_after_solving_pressure_part; - nth_darcy_solution_after_solving_pressure_part = tmp_darcy[0]; - n_minus_oneth_darcy_solution_after_solving_pressure_part = tmp_darcy[1]; + std::vector x_darcy (2); + x_darcy[0] = nth_darcy_solution_after_solving_pressure_part; + x_darcy[1] = n_minus_oneth_darcy_solution_after_solving_pressure_part; - rebuild_saturation_matrix = true; - } - else - { - rebuild_saturation_matrix = false; - - std::vector darcy_block_component (dim+1,0); - darcy_block_component[dim] = 1; - - std::vector darcy_dofs_per_block (2); - DoFTools::count_dofs_per_block (darcy_dof_handler, darcy_dofs_per_block, darcy_block_component); - const unsigned int n_u = darcy_dofs_per_block[0], - n_p = darcy_dofs_per_block[1], - n_s = saturation_dof_handler.n_dofs(); - - std::cout << "Number of active cells: " - << triangulation.n_active_cells() - << " (on " - << triangulation.n_levels() - << " levels)" - << std::endl - << "Number of degrees of freedom: " - << n_u + n_p + n_s - << " (" << n_u << '+' << n_p << '+'<< n_s <<')' - << std::endl - << std::endl; - } + SolutionTransfer saturation_soltrans(saturation_dof_handler); -} + SolutionTransfer darcy_soltrans(darcy_dof_handler); + triangulation.prepare_coarsening_and_refinement(); + saturation_soltrans.prepare_for_coarsening_and_refinement(x_saturation); - // @sect3{TwoPhaseFlowProblem::output_results} + darcy_soltrans.prepare_for_coarsening_and_refinement(x_darcy); - // This function to process the output - // data. We only store the results when we - // actually solve the pressure and velocity - // part at the present time step. The rest of - // the implementation is similar to that - // output function in step-31, which - // implementations has been explained in that - // tutorial. -template -void TwoPhaseFlowProblem::output_results () const -{ - if ( solve_pressure_velocity_part == false ) - return; + triangulation.execute_coarsening_and_refinement (); + setup_dofs (); - const FESystem joint_fe (darcy_fe, 1, - saturation_fe, 1); - DoFHandler joint_dof_handler (triangulation); - joint_dof_handler.distribute_dofs (joint_fe); - Assert (joint_dof_handler.n_dofs() == - darcy_dof_handler.n_dofs() + saturation_dof_handler.n_dofs(), - ExcInternalError()); + std::vector tmp_saturation (3); + tmp_saturation[0].reinit (saturation_solution); + tmp_saturation[1].reinit (saturation_solution); + tmp_saturation[2].reinit (saturation_solution); + saturation_soltrans.interpolate(x_saturation, tmp_saturation); - Vector joint_solution (joint_dof_handler.n_dofs()); + saturation_solution = tmp_saturation[0]; + old_saturation_solution = tmp_saturation[1]; + nth_saturation_solution_after_solving_pressure_part = tmp_saturation[2]; - { - std::vector local_joint_dof_indices (joint_fe.dofs_per_cell); - std::vector local_darcy_dof_indices (darcy_fe.dofs_per_cell); - std::vector local_saturation_dof_indices (saturation_fe.dofs_per_cell); + std::vector tmp_darcy (2); + tmp_darcy[0].reinit (darcy_solution); + tmp_darcy[1].reinit (darcy_solution); + darcy_soltrans.interpolate(x_darcy, tmp_darcy); - typename DoFHandler::active_cell_iterator - joint_cell = joint_dof_handler.begin_active(), - joint_endc = joint_dof_handler.end(), - darcy_cell = darcy_dof_handler.begin_active(), - saturation_cell = saturation_dof_handler.begin_active(); + nth_darcy_solution_after_solving_pressure_part = tmp_darcy[0]; + n_minus_oneth_darcy_solution_after_solving_pressure_part = tmp_darcy[1]; - for (; joint_cell!=joint_endc; ++joint_cell, ++darcy_cell, ++saturation_cell) + rebuild_saturation_matrix = true; + } + else { - joint_cell->get_dof_indices (local_joint_dof_indices); - darcy_cell->get_dof_indices (local_darcy_dof_indices); - saturation_cell->get_dof_indices (local_saturation_dof_indices); - - for (unsigned int i=0; i darcy_block_component (dim+1,0); + darcy_block_component[dim] = 1; + + std::vector darcy_dofs_per_block (2); + DoFTools::count_dofs_per_block (darcy_dof_handler, darcy_dofs_per_block, darcy_block_component); + const unsigned int n_u = darcy_dofs_per_block[0], + n_p = darcy_dofs_per_block[1], + n_s = saturation_dof_handler.n_dofs(); + + std::cout << "Number of active cells: " + << triangulation.n_active_cells() + << " (on " + << triangulation.n_levels() + << " levels)" + << std::endl + << "Number of degrees of freedom: " + << n_u + n_p + n_s + << " (" << n_u << '+' << n_p << '+'<< n_s <<')' + << std::endl + << std::endl; } + } - std::vector joint_solution_names; - switch (dim) - { - case 2: - joint_solution_names.push_back ("u"); - joint_solution_names.push_back ("v"); - break; - - case 3: - joint_solution_names.push_back ("u"); - joint_solution_names.push_back ("v"); - joint_solution_names.push_back ("w"); - break; - - default: - Assert (false, ExcNotImplemented()); - } - joint_solution_names.push_back ("pressure"); - joint_solution_names.push_back ("saturation"); - std::vector - data_component_interpretation - (dim, DataComponentInterpretation::component_is_part_of_vector); - data_component_interpretation - .push_back (DataComponentInterpretation::component_is_scalar); - data_component_interpretation - .push_back (DataComponentInterpretation::component_is_scalar); - DataOut data_out; - data_out.attach_dof_handler (joint_dof_handler); - data_out.add_data_vector (joint_solution, joint_solution_names, - DataOut::type_dof_data, - data_component_interpretation); + // @sect3{TwoPhaseFlowProblem::output_results} - data_out.build_patches (); + // This function to process the output + // data. We only store the results when we + // actually solve the pressure and velocity + // part at the present time step. The rest of + // the implementation is similar to that + // output function in step-31, which + // implementations has been explained in that + // tutorial. + template + void TwoPhaseFlowProblem::output_results () const + { + if ( solve_pressure_velocity_part == false ) + return; - std::string filename = "solution-" + - Utilities::int_to_string (timestep_number, 5) + ".tec"; - std::ofstream output (filename.c_str()); - data_out.write_tecplot (output); -} + const FESystem joint_fe (darcy_fe, 1, + saturation_fe, 1); + DoFHandler joint_dof_handler (triangulation); + joint_dof_handler.distribute_dofs (joint_fe); + Assert (joint_dof_handler.n_dofs() == + darcy_dof_handler.n_dofs() + saturation_dof_handler.n_dofs(), + ExcInternalError()); + Vector joint_solution (joint_dof_handler.n_dofs()); + { + std::vector local_joint_dof_indices (joint_fe.dofs_per_cell); + std::vector local_darcy_dof_indices (darcy_fe.dofs_per_cell); + std::vector local_saturation_dof_indices (saturation_fe.dofs_per_cell); - // @sect3{TwoPhaseFlowProblem::THE_REMAINING_FUNCTIONS} + typename DoFHandler::active_cell_iterator + joint_cell = joint_dof_handler.begin_active(), + joint_endc = joint_dof_handler.end(), + darcy_cell = darcy_dof_handler.begin_active(), + saturation_cell = saturation_dof_handler.begin_active(); + + for (; joint_cell!=joint_endc; ++joint_cell, ++darcy_cell, ++saturation_cell) + { + joint_cell->get_dof_indices (local_joint_dof_indices); + darcy_cell->get_dof_indices (local_darcy_dof_indices); + saturation_cell->get_dof_indices (local_saturation_dof_indices); + + for (unsigned int i=0; i -void -TwoPhaseFlowProblem::project_back_saturation () -{ - for (unsigned int i=0; i 1) - saturation_solution(i) = 1; -} + } + } + std::vector joint_solution_names; + switch (dim) + { + case 2: + joint_solution_names.push_back ("u"); + joint_solution_names.push_back ("v"); + break; + + case 3: + joint_solution_names.push_back ("u"); + joint_solution_names.push_back ("v"); + joint_solution_names.push_back ("w"); + break; + + default: + Assert (false, ExcNotImplemented()); + } + joint_solution_names.push_back ("pressure"); + joint_solution_names.push_back ("saturation"); + std::vector + data_component_interpretation + (dim, DataComponentInterpretation::component_is_part_of_vector); + data_component_interpretation + .push_back (DataComponentInterpretation::component_is_scalar); + data_component_interpretation + .push_back (DataComponentInterpretation::component_is_scalar); -template -double -TwoPhaseFlowProblem::get_maximal_velocity () const -{ - QGauss quadrature_formula(darcy_degree+2); - const unsigned int n_q_points - = quadrature_formula.size(); - - FEValues darcy_fe_values (darcy_fe, quadrature_formula, - update_values); - std::vector > darcy_solution_values(n_q_points, - Vector(dim+1)); - double max_velocity = 0; - - typename DoFHandler::active_cell_iterator - cell = darcy_dof_handler.begin_active(), - endc = darcy_dof_handler.end(); - for (; cell!=endc; ++cell) - { - darcy_fe_values.reinit (cell); - darcy_fe_values.get_function_values (darcy_solution, darcy_solution_values); - - for (unsigned int q=0; q velocity; - for (unsigned int i=0; i data_out; - return max_velocity; -} + data_out.attach_dof_handler (joint_dof_handler); + data_out.add_data_vector (joint_solution, joint_solution_names, + DataOut::type_dof_data, + data_component_interpretation); + data_out.build_patches (); -template -std::pair -TwoPhaseFlowProblem::get_extrapolated_saturation_range () const -{ - const QGauss quadrature_formula(saturation_degree+2); - const unsigned int n_q_points = quadrature_formula.size(); + std::string filename = "solution-" + + Utilities::int_to_string (timestep_number, 5) + ".tec"; + std::ofstream output (filename.c_str()); + data_out.write_tecplot (output); + } - FEValues fe_values (saturation_fe, quadrature_formula, - update_values); - std::vector old_saturation_values(n_q_points); - std::vector old_old_saturation_values(n_q_points); - if (timestep_number != 0) - { - double min_saturation = (1. + time_step/old_time_step) * - old_saturation_solution.linfty_norm() - + - time_step/old_time_step * - old_old_saturation_solution.linfty_norm(), - max_saturation = -min_saturation; - typename DoFHandler::active_cell_iterator - cell = saturation_dof_handler.begin_active(), - endc = saturation_dof_handler.end(); - for (; cell!=endc; ++cell) - { - fe_values.reinit (cell); - fe_values.get_function_values (old_saturation_solution, - old_saturation_values); - fe_values.get_function_values (old_old_saturation_solution, - old_old_saturation_values); - - for (unsigned int q=0; q::THE_REMAINING_FUNCTIONS} - typename DoFHandler::active_cell_iterator - cell = saturation_dof_handler.begin_active(), - endc = saturation_dof_handler.end(); - for (; cell!=endc; ++cell) - { - fe_values.reinit (cell); - fe_values.get_function_values (old_saturation_solution, - old_saturation_values); + // The remaining functions that have been + // used in step-31 so we don't have to + // describe their implementations. + template + void + TwoPhaseFlowProblem::project_back_saturation () + { + for (unsigned int i=0; i 1) + saturation_solution(i) = 1; + } - for (unsigned int q=0; q + double + TwoPhaseFlowProblem::get_maximal_velocity () const + { + QGauss quadrature_formula(darcy_degree+2); + const unsigned int n_q_points + = quadrature_formula.size(); - return std::make_pair(min_saturation, max_saturation); - } -} + FEValues darcy_fe_values (darcy_fe, quadrature_formula, + update_values); + std::vector > darcy_solution_values(n_q_points, + Vector(dim+1)); + double max_velocity = 0; -template -double -TwoPhaseFlowProblem:: -compute_viscosity (const std::vector &old_saturation, - const std::vector &old_old_saturation, - const std::vector > &old_saturation_grads, - const std::vector > &old_old_saturation_grads, - const std::vector > &present_darcy_values, - const double global_u_infty, - const double global_S_variation, - const double global_Omega_diameter, - const double cell_diameter, - const double old_time_step, - const double viscosity) -{ - const double beta = 0.08 * dim; - const double alpha = 1; + typename DoFHandler::active_cell_iterator + cell = darcy_dof_handler.begin_active(), + endc = darcy_dof_handler.end(); + for (; cell!=endc; ++cell) + { + darcy_fe_values.reinit (cell); + darcy_fe_values.get_function_values (darcy_solution, darcy_solution_values); + + for (unsigned int q=0; q velocity; + for (unsigned int i=0; i + std::pair + TwoPhaseFlowProblem::get_extrapolated_saturation_range () const + { + const QGauss quadrature_formula(saturation_degree+2); + const unsigned int n_q_points = quadrature_formula.size(); - double max_residual = 0; - double max_velocity = 0; + FEValues fe_values (saturation_fe, quadrature_formula, + update_values); + std::vector old_saturation_values(n_q_points); + std::vector old_old_saturation_values(n_q_points); - for (unsigned int q=0; q < n_q_points; ++q) - { - Tensor<1,dim> u; - for (unsigned int d=0; d::active_cell_iterator + cell = saturation_dof_handler.begin_active(), + endc = saturation_dof_handler.end(); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + fe_values.get_function_values (old_saturation_solution, + old_saturation_values); + fe_values.get_function_values (old_old_saturation_solution, + old_old_saturation_values); + + for (unsigned int q=0; q::active_cell_iterator + cell = saturation_dof_handler.begin_active(), + endc = saturation_dof_handler.end(); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + fe_values.get_function_values (old_saturation_solution, + old_saturation_values); + + for (unsigned int q=0; q + double + TwoPhaseFlowProblem:: + compute_viscosity (const std::vector &old_saturation, + const std::vector &old_old_saturation, + const std::vector > &old_saturation_grads, + const std::vector > &old_old_saturation_grads, + const std::vector > &present_darcy_values, + const double global_u_infty, + const double global_S_variation, + const double global_Omega_diameter, + const double cell_diameter, + const double old_time_step, + const double viscosity) + { + const double beta = 0.08 * dim; + const double alpha = 1; - const double dF_dS = get_fractional_flow_derivative ((old_saturation[q] + old_old_saturation[q]) / 2.0, - viscosity); + if (global_u_infty == 0) + return 5e-3 * cell_diameter; - const double u_grad_S = u * dF_dS * - (old_saturation_grads[q] + old_old_saturation_grads[q]) / 2.0; + const unsigned int n_q_points = old_saturation.size(); - const double residual - = std::abs((dS_dt + u_grad_S) * - std::pow((old_saturation[q]+old_old_saturation[q]) / 2, - alpha-1.)); + double max_residual = 0; + double max_velocity = 0; - max_residual = std::max (residual, max_residual); - max_velocity = std::max (std::sqrt (u*u), max_velocity); - } + for (unsigned int q=0; q < n_q_points; ++q) + { + Tensor<1,dim> u; + for (unsigned int d=0; d::run} - - // In this function, we follow the structure - // of the same function partly in step-21 and - // partly in step-31 so again there is no - // need to repeat it. However, since we - // consider the simulation with grid - // adaptivity, we need to compute a - // saturation predictor, which implementation - // was first used in step-33, for the - // function that computes the refinement - // indicators. -template -void TwoPhaseFlowProblem::run () -{ - unsigned int pre_refinement_step = 0; + const double residual + = std::abs((dS_dt + u_grad_S) * + std::pow((old_saturation[q]+old_old_saturation[q]) / 2, + alpha-1.)); - GridGenerator::hyper_cube (triangulation, 0, 1); - triangulation.refine_global (n_refinement_steps); + max_residual = std::max (residual, max_residual); + max_velocity = std::max (std::sqrt (u*u), max_velocity); + } - setup_dofs (); + const double global_scaling = global_u_infty * global_S_variation / + std::pow(global_Omega_diameter, alpha - 2.); - start_time_iteration: + return (beta * + max_velocity * + std::min (cell_diameter, + std::pow(cell_diameter,alpha) * + max_residual / global_scaling)); + } - VectorTools::project (saturation_dof_handler, - saturation_constraints, - QGauss(saturation_degree+2), - SaturationInitialValues(), - old_saturation_solution); - timestep_number = 0; - double time = 0; + // @sect3{TwoPhaseFlowProblem::run} - do - { - std::cout << "Timestep " << timestep_number - << ": t=" << time - << ", dt=" << time_step - << std::endl; + // In this function, we follow the structure + // of the same function partly in step-21 and + // partly in step-31 so again there is no + // need to repeat it. However, since we + // consider the simulation with grid + // adaptivity, we need to compute a + // saturation predictor, which implementation + // was first used in step-33, for the + // function that computes the refinement + // indicators. + template + void TwoPhaseFlowProblem::run () + { + unsigned int pre_refinement_step = 0; - solve (); + GridGenerator::hyper_cube (triangulation, 0, 1); + triangulation.refine_global (n_refinement_steps); - output_results (); + setup_dofs (); - solve_pressure_velocity_part = false; + start_time_iteration: - if ((timestep_number == 0) && - (pre_refinement_step < saturation_level)) - { - predictor_saturation_solution = saturation_solution; - predictor_saturation_solution.sadd (2.0, -1.0, old_saturation_solution); - Vector refinement_indicators (triangulation.n_active_cells()); - compute_refinement_indicators(refinement_indicators); - refine_grid(refinement_indicators); - ++pre_refinement_step; - goto start_time_iteration; - } - else - { - predictor_saturation_solution = saturation_solution; - predictor_saturation_solution.sadd (2.0, -1.0, old_saturation_solution); - Vector refinement_indicators (triangulation.n_active_cells()); - compute_refinement_indicators(refinement_indicators); - refine_grid(refinement_indicators); - } + VectorTools::project (saturation_dof_handler, + saturation_constraints, + QGauss(saturation_degree+2), + SaturationInitialValues(), + old_saturation_solution); - time += time_step; - ++timestep_number; + timestep_number = 0; + double time = 0; - old_old_saturation_solution = old_saturation_solution; - old_saturation_solution = saturation_solution; + do + { + std::cout << "Timestep " << timestep_number + << ": t=" << time + << ", dt=" << time_step + << std::endl; + + solve (); + + output_results (); + + solve_pressure_velocity_part = false; + + if ((timestep_number == 0) && + (pre_refinement_step < saturation_level)) + { + predictor_saturation_solution = saturation_solution; + predictor_saturation_solution.sadd (2.0, -1.0, old_saturation_solution); + Vector refinement_indicators (triangulation.n_active_cells()); + compute_refinement_indicators(refinement_indicators); + refine_grid(refinement_indicators); + ++pre_refinement_step; + goto start_time_iteration; + } + else + { + predictor_saturation_solution = saturation_solution; + predictor_saturation_solution.sadd (2.0, -1.0, old_saturation_solution); + Vector refinement_indicators (triangulation.n_active_cells()); + compute_refinement_indicators(refinement_indicators); + refine_grid(refinement_indicators); + } - } - while (time <= 250); + time += time_step; + ++timestep_number; + + old_old_saturation_solution = old_saturation_solution; + old_saturation_solution = saturation_solution; + + } + while (time <= 250); + } } + int main () { try { + using namespace dealii; + using namespace Step43; + deallog.depth_console (0); TwoPhaseFlowProblem<3> two_phase_flow_problem(1); diff --git a/deal.II/examples/step-44/step-44.cc b/deal.II/examples/step-44/step-44.cc index e6b2b68d73..426a571184 100644 --- a/deal.II/examples/step-44/step-44.cc +++ b/deal.II/examples/step-44/step-44.cc @@ -53,25 +53,28 @@ #include #include -using namespace dealii; -// @sect3{Run-time parameters} -namespace Parameters +namespace Step44 { + using namespace dealii; + +// @sect3{Run-time parameters} + namespace Parameters + { // Finite Element system -struct FESystem -{ - int poly_degree; - int quad_order; + struct FESystem + { + int poly_degree; + int quad_order; - static void declare_parameters (ParameterHandler &prm); - void parse_parameters (ParameterHandler &prm); -}; + static void declare_parameters (ParameterHandler &prm); + void parse_parameters (ParameterHandler &prm); + }; -void FESystem::declare_parameters (ParameterHandler &prm) -{ - prm.enter_subsection("Finite element system"); + void FESystem::declare_parameters (ParameterHandler &prm) { + prm.enter_subsection("Finite element system"); + { prm.declare_entry("Polynomial degree", "1", Patterns::Integer(), @@ -81,35 +84,35 @@ void FESystem::declare_parameters (ParameterHandler &prm) "2", Patterns::Integer(), "Gauss quadrature order"); + } + prm.leave_subsection(); } - prm.leave_subsection(); -} -void FESystem::parse_parameters (ParameterHandler &prm) -{ - prm.enter_subsection("Finite element system"); + void FESystem::parse_parameters (ParameterHandler &prm) { + prm.enter_subsection("Finite element system"); + { poly_degree = prm.get_integer("Polynomial degree"); quad_order = prm.get_integer("Quadrature order"); + } + prm.leave_subsection(); } - prm.leave_subsection(); -} // Geometry -struct Geometry -{ - int global_refinement; - double scale; - double p_p0; + struct Geometry + { + int global_refinement; + double scale; + double p_p0; - static void declare_parameters (ParameterHandler &prm); - void parse_parameters (ParameterHandler &prm); -}; + static void declare_parameters (ParameterHandler &prm); + void parse_parameters (ParameterHandler &prm); + }; -void Geometry::declare_parameters (ParameterHandler &prm) -{ - prm.enter_subsection("Geometry"); + void Geometry::declare_parameters (ParameterHandler &prm) { + prm.enter_subsection("Geometry"); + { prm.declare_entry("Global refinement", "2", Patterns::Integer(), @@ -124,35 +127,35 @@ void Geometry::declare_parameters (ParameterHandler &prm) "40", Patterns::Selection("20|40|60|80|100"), "Ratio of applied pressure to reference pressure"); + } + prm.leave_subsection(); } - prm.leave_subsection(); -} -void Geometry::parse_parameters (ParameterHandler &prm) -{ - prm.enter_subsection("Geometry"); + void Geometry::parse_parameters (ParameterHandler &prm) { + prm.enter_subsection("Geometry"); + { global_refinement = prm.get_integer("Global refinement"); scale = prm.get_double("Grid scale"); p_p0= prm.get_double("Pressure ratio p/p0"); + } + prm.leave_subsection(); } - prm.leave_subsection(); -} // Materials -struct Materials -{ - double nu; - double mu; + struct Materials + { + double nu; + double mu; - static void declare_parameters (ParameterHandler &prm); - void parse_parameters (ParameterHandler &prm); -}; + static void declare_parameters (ParameterHandler &prm); + void parse_parameters (ParameterHandler &prm); + }; -void Materials::declare_parameters (ParameterHandler &prm) -{ - prm.enter_subsection("Material properties"); + void Materials::declare_parameters (ParameterHandler &prm) { + prm.enter_subsection("Material properties"); + { prm.declare_entry("Poisson's ratio", "0.49", Patterns::Double(), @@ -162,36 +165,36 @@ void Materials::declare_parameters (ParameterHandler &prm) "1.0e6", Patterns::Double(), "Shear modulus"); + } + prm.leave_subsection(); } - prm.leave_subsection(); -} -void Materials::parse_parameters (ParameterHandler &prm) -{ - prm.enter_subsection("Material properties"); + void Materials::parse_parameters (ParameterHandler &prm) { + prm.enter_subsection("Material properties"); + { nu = prm.get_double("Poisson's ratio"); mu = prm.get_double("Shear modulus"); + } + prm.leave_subsection(); } - prm.leave_subsection(); -} // Linear solver -struct LinearSolver -{ - std::string type_lin; - double tol_lin; - double max_iterations_lin; - double ssor_relaxation; + struct LinearSolver + { + std::string type_lin; + double tol_lin; + double max_iterations_lin; + double ssor_relaxation; - static void declare_parameters (ParameterHandler &prm); - void parse_parameters (ParameterHandler &prm); -}; + static void declare_parameters (ParameterHandler &prm); + void parse_parameters (ParameterHandler &prm); + }; -void LinearSolver::declare_parameters (ParameterHandler &prm) -{ - prm.enter_subsection("Linear solver"); + void LinearSolver::declare_parameters (ParameterHandler &prm) { + prm.enter_subsection("Linear solver"); + { prm.declare_entry("Solver type", "CG", Patterns::Selection("CG|Direct"), @@ -211,37 +214,37 @@ void LinearSolver::declare_parameters (ParameterHandler &prm) "0.6", Patterns::Double(), "SSOR preconditioner relaxation value"); + } + prm.leave_subsection(); } - prm.leave_subsection(); -} -void LinearSolver::parse_parameters (ParameterHandler &prm) -{ - prm.enter_subsection("Linear solver"); + void LinearSolver::parse_parameters (ParameterHandler &prm) { + prm.enter_subsection("Linear solver"); + { type_lin = prm.get("Solver type"); tol_lin = prm.get_double("Residual"); max_iterations_lin = prm.get_double("Max iteration multiplier"); ssor_relaxation = prm.get_double("SSOR Relaxation"); + } + prm.leave_subsection(); } - prm.leave_subsection(); -} // Nonlinear solver -struct NonlinearSolver -{ - unsigned int max_iterations_NR; - double tol_f; - double tol_u; + struct NonlinearSolver + { + unsigned int max_iterations_NR; + double tol_f; + double tol_u; - static void declare_parameters (ParameterHandler &prm); - void parse_parameters (ParameterHandler &prm); -}; + static void declare_parameters (ParameterHandler &prm); + void parse_parameters (ParameterHandler &prm); + }; -void NonlinearSolver::declare_parameters (ParameterHandler &prm) -{ - prm.enter_subsection("Nonlinear solver"); + void NonlinearSolver::declare_parameters (ParameterHandler &prm) { + prm.enter_subsection("Nonlinear solver"); + { prm.declare_entry("Max iterations Newton-Raphson", "10", Patterns::Integer(), @@ -256,35 +259,35 @@ void NonlinearSolver::declare_parameters (ParameterHandler &prm) "1.0e-3", Patterns::Double(), "Displacement error tolerance"); + } + prm.leave_subsection(); } - prm.leave_subsection(); -} -void NonlinearSolver::parse_parameters (ParameterHandler &prm) -{ - prm.enter_subsection("Nonlinear solver"); + void NonlinearSolver::parse_parameters (ParameterHandler &prm) { + prm.enter_subsection("Nonlinear solver"); + { max_iterations_NR = prm.get_integer("Max iterations Newton-Raphson"); tol_f = prm.get_double("Tolerance force"); tol_u = prm.get_double("Tolerance displacement"); + } + prm.leave_subsection(); } - prm.leave_subsection(); -} // Time -struct Time -{ - double end_time; - double delta_t; + struct Time + { + double end_time; + double delta_t; - static void declare_parameters (ParameterHandler &prm); - void parse_parameters (ParameterHandler &prm); -}; + static void declare_parameters (ParameterHandler &prm); + void parse_parameters (ParameterHandler &prm); + }; -void Time::declare_parameters (ParameterHandler &prm) -{ - prm.enter_subsection("Time"); + void Time::declare_parameters (ParameterHandler &prm) { + prm.enter_subsection("Time"); + { prm.declare_entry("End time", "1", Patterns::Double(), @@ -294,23 +297,23 @@ void Time::declare_parameters (ParameterHandler &prm) "0.1", Patterns::Double(), "Time step size"); + } + prm.leave_subsection(); } - prm.leave_subsection(); -} -void Time::parse_parameters (ParameterHandler &prm) -{ - prm.enter_subsection("Time"); + void Time::parse_parameters (ParameterHandler &prm) { + prm.enter_subsection("Time"); + { end_time = prm.get_double("End time"); delta_t = prm.get_double("Time step size"); + } + prm.leave_subsection(); } - prm.leave_subsection(); -} // All parameters -struct AllParameters - : + struct AllParameters + : public FESystem, public Geometry, public Materials, @@ -318,683 +321,683 @@ struct AllParameters public NonlinearSolver, public Time -{ - AllParameters (const std::string & input_file); + { + AllParameters (const std::string & input_file); - static void declare_parameters (ParameterHandler &prm); - void parse_parameters (ParameterHandler &prm); -}; + static void declare_parameters (ParameterHandler &prm); + void parse_parameters (ParameterHandler &prm); + }; -AllParameters::AllParameters (const std::string & input_file) -{ - ParameterHandler prm; - declare_parameters(prm); - prm.read_input (input_file); - parse_parameters(prm); -} + AllParameters::AllParameters (const std::string & input_file) + { + ParameterHandler prm; + declare_parameters(prm); + prm.read_input (input_file); + parse_parameters(prm); + } -void AllParameters::declare_parameters (ParameterHandler &prm) -{ - FESystem::declare_parameters(prm); - Geometry::declare_parameters(prm); - Materials::declare_parameters(prm); - LinearSolver::declare_parameters(prm); - NonlinearSolver::declare_parameters(prm); - Time::declare_parameters(prm); -} + void AllParameters::declare_parameters (ParameterHandler &prm) + { + FESystem::declare_parameters(prm); + Geometry::declare_parameters(prm); + Materials::declare_parameters(prm); + LinearSolver::declare_parameters(prm); + NonlinearSolver::declare_parameters(prm); + Time::declare_parameters(prm); + } -void AllParameters::parse_parameters (ParameterHandler &prm) -{ - FESystem::parse_parameters(prm); - Geometry::parse_parameters(prm); - Materials::parse_parameters(prm); - LinearSolver::parse_parameters(prm); - NonlinearSolver::parse_parameters(prm); - Time::parse_parameters(prm); -} -} + void AllParameters::parse_parameters (ParameterHandler &prm) + { + FESystem::parse_parameters(prm); + Geometry::parse_parameters(prm); + Materials::parse_parameters(prm); + LinearSolver::parse_parameters(prm); + NonlinearSolver::parse_parameters(prm); + Time::parse_parameters(prm); + } + } // @sect3{General tools} -namespace AdditionalTools -{ -template -void extract_submatrix(const std::vector< unsigned int > &row_index_set, - const std::vector< unsigned int > &column_index_set, - const MatrixType &matrix, - FullMatrix< double > &sub_matrix ) -{ + namespace AdditionalTools + { + template + void extract_submatrix(const std::vector< unsigned int > &row_index_set, + const std::vector< unsigned int > &column_index_set, + const MatrixType &matrix, + FullMatrix< double > &sub_matrix ) + { - const unsigned int n_rows_submatrix = row_index_set.size(); - const unsigned int n_cols_submatrix = column_index_set.size(); + const unsigned int n_rows_submatrix = row_index_set.size(); + const unsigned int n_cols_submatrix = column_index_set.size(); - sub_matrix.reinit(n_rows_submatrix, n_cols_submatrix); + sub_matrix.reinit(n_rows_submatrix, n_cols_submatrix); - for (unsigned int sub_row = 0; sub_row < n_rows_submatrix; ++sub_row) { + for (unsigned int sub_row = 0; sub_row < n_rows_submatrix; ++sub_row) { const unsigned int row = row_index_set[sub_row]; Assert (row<=matrix.m(), ExcIndexRange(row, 0, matrix.m())); for (unsigned int sub_col = 0; sub_col < n_cols_submatrix; ++sub_col) { - const unsigned int col = column_index_set[sub_col]; - Assert (col<=matrix.n(), ExcIndexRange(col, 0, matrix.n())); + const unsigned int col = column_index_set[sub_col]; + Assert (col<=matrix.n(), ExcIndexRange(col, 0, matrix.n())); - sub_matrix(sub_row,sub_col) = matrix(row, col); + sub_matrix(sub_row,sub_col) = matrix(row, col); } + } } -} -template -void replace_submatrix(const std::vector< unsigned int > &row_index_set, - const std::vector< unsigned int > &column_index_set, - const MatrixType &sub_matrix, - FullMatrix< double > &matrix) -{ - const unsigned int n_rows_submatrix = row_index_set.size(); - Assert (n_rows_submatrix<=sub_matrix.m(), ExcIndexRange(n_rows_submatrix, 0, sub_matrix.m())); - const unsigned int n_cols_submatrix = column_index_set.size(); - Assert (n_cols_submatrix<=sub_matrix.n(), ExcIndexRange(n_cols_submatrix, 0, sub_matrix.n())); + template + void replace_submatrix(const std::vector< unsigned int > &row_index_set, + const std::vector< unsigned int > &column_index_set, + const MatrixType &sub_matrix, + FullMatrix< double > &matrix) + { + const unsigned int n_rows_submatrix = row_index_set.size(); + Assert (n_rows_submatrix<=sub_matrix.m(), ExcIndexRange(n_rows_submatrix, 0, sub_matrix.m())); + const unsigned int n_cols_submatrix = column_index_set.size(); + Assert (n_cols_submatrix<=sub_matrix.n(), ExcIndexRange(n_cols_submatrix, 0, sub_matrix.n())); - for (unsigned int sub_row = 0; sub_row < n_rows_submatrix; ++sub_row) { + for (unsigned int sub_row = 0; sub_row < n_rows_submatrix; ++sub_row) { const unsigned int row = row_index_set[sub_row]; Assert (row<=matrix.m(), ExcIndexRange(row, 0, matrix.m())); for (unsigned int sub_col = 0; sub_col < n_cols_submatrix; ++sub_col) { - const unsigned int col = column_index_set[sub_col]; - Assert (col<=matrix.n(), ExcIndexRange(col, 0, matrix.n())); + const unsigned int col = column_index_set[sub_col]; + Assert (col<=matrix.n(), ExcIndexRange(col, 0, matrix.n())); - matrix(row, col) = sub_matrix(sub_row, sub_col); + matrix(row, col) = sub_matrix(sub_row, sub_col); } + } } -} -} + } // @sect3{Time class} -class Time { -public: - Time (const double & time_end, - const double & delta_t) - : - timestep (0), - time_current (0.0), - time_end (time_end), - delta_t (delta_t) - {} - virtual ~Time (void) {} - - const double & current (void) const {return time_current;} - const double & end (void) const {return time_end;} - const double & get_delta_t (void) const {return delta_t;} - const unsigned int & get_timestep (void) const {return timestep;} - void increment (void) {time_current += delta_t; ++timestep;} - -private: - unsigned int timestep; - double time_current; - const double time_end; - const double delta_t; -}; + class Time { + public: + Time (const double & time_end, + const double & delta_t) + : + timestep (0), + time_current (0.0), + time_end (time_end), + delta_t (delta_t) + {} + virtual ~Time (void) {} + + const double & current (void) const {return time_current;} + const double & end (void) const {return time_end;} + const double & get_delta_t (void) const {return delta_t;} + const unsigned int & get_timestep (void) const {return timestep;} + void increment (void) {time_current += delta_t; ++timestep;} + + private: + unsigned int timestep; + double time_current; + const double time_end; + const double delta_t; + }; // @sect3{Neo-Hookean material} -template -class Material_NH -{ -public: - /// \brief Class constructor - Material_NH (const double & lambda, - const double & mu) - : - lambda_0 (lambda), - mu_0 (mu), - kappa_0 (lambda + 2.0/3.0*mu) - { } - virtual ~Material_NH (void) {}; - - // Stress and constitutive tensors - virtual SymmetricTensor<2, dim> get_T (const double & J, - const SymmetricTensor <2, dim> & B) - { - const double dW_dJ = get_dU_dtheta (J); - return mu_0*B + dW_dJ*J*I; - } + template + class Material_NH + { + public: + /// \brief Class constructor + Material_NH (const double & lambda, + const double & mu) + : + lambda_0 (lambda), + mu_0 (mu), + kappa_0 (lambda + 2.0/3.0*mu) + { } + virtual ~Material_NH (void) {}; - virtual SymmetricTensor<4, dim> get_JC (const double & J, - const SymmetricTensor <2, dim> & B) - { - const double dW_dJ = get_dU_dtheta (J); - const double d2W_dJ2 = get_d2U_dtheta2 (J); - return J*( (dW_dJ + J*d2W_dJ2)*IxI - (2.0*dW_dJ)*II ); - } + // Stress and constitutive tensors + virtual SymmetricTensor<2, dim> get_T (const double & J, + const SymmetricTensor <2, dim> & B) + { + const double dW_dJ = get_dU_dtheta (J); + return mu_0*B + dW_dJ*J*I; + } + + virtual SymmetricTensor<4, dim> get_JC (const double & J, + const SymmetricTensor <2, dim> & B) + { + const double dW_dJ = get_dU_dtheta (J); + const double d2W_dJ2 = get_d2U_dtheta2 (J); + return J*( (dW_dJ + J*d2W_dJ2)*IxI - (2.0*dW_dJ)*II ); + } - // Volumetric quantities methods - double get_dU_dtheta (const double & d) {return kappa_0*(d - 1.0/d);} - double get_d2U_dtheta2 (const double & d) {return kappa_0*(1.0 + 1.0/(d*d));} + // Volumetric quantities methods + double get_dU_dtheta (const double & d) {return kappa_0*(d - 1.0/d);} + double get_d2U_dtheta2 (const double & d) {return kappa_0*(1.0 + 1.0/(d*d));} -protected: - // Material properties - const double lambda_0; // Lame modulus - const double mu_0; // Shear modulus - const double kappa_0; // Bulk modulus + protected: + // Material properties + const double lambda_0; // Lame modulus + const double mu_0; // Shear modulus + const double kappa_0; // Bulk modulus - static SymmetricTensor<2, dim> const I; - static SymmetricTensor<4, dim> const IxI; - static SymmetricTensor<4, dim> const II; -}; + static SymmetricTensor<2, dim> const I; + static SymmetricTensor<4, dim> const IxI; + static SymmetricTensor<4, dim> const II; + }; -template SymmetricTensor<2, dim> const Material_NH::I = SymmetricTensor<2, dim> (unit_symmetric_tensor ()); -template SymmetricTensor<4, dim> const Material_NH::IxI = SymmetricTensor<4, dim> (outer_product (I, I)); -template SymmetricTensor<4, dim> const Material_NH::II = SymmetricTensor<4, dim> (identity_tensor ()); + template SymmetricTensor<2, dim> const Material_NH::I = SymmetricTensor<2, dim> (unit_symmetric_tensor ()); + template SymmetricTensor<4, dim> const Material_NH::IxI = SymmetricTensor<4, dim> (outer_product (I, I)); + template SymmetricTensor<4, dim> const Material_NH::II = SymmetricTensor<4, dim> (identity_tensor ()); // @sect3{Quadrature point history} -template -class PointHistory -{ -public: - PointHistory (void) - : - material (NULL), - dilatation_n (1.0), - pressure_n (0.0) - { } - virtual ~PointHistory (void) {delete material;} - - void setup_lqp ( Parameters::AllParameters & parameters ) - { - const double lambda = 2.0*parameters.mu*parameters.nu / (1.0-2.0*parameters.nu); - material = new Material_NH (lambda, - parameters.mu); + template + class PointHistory + { + public: + PointHistory (void) + : + material (NULL), + dilatation_n (1.0), + pressure_n (0.0) + { } + virtual ~PointHistory (void) {delete material;} - // Initialise all tensors correctly - update_values (Tensor <2,dim> (), 0.0, 1.0); - } + void setup_lqp ( Parameters::AllParameters & parameters ) + { + const double lambda = 2.0*parameters.mu*parameters.nu / (1.0-2.0*parameters.nu); + material = new Material_NH (lambda, + parameters.mu); - // Total Variables - void update_values (const Tensor<2, dim> & grad_u_n, - const double & pressure, - const double & dilatation) - { - // Calculated variables from displacement, displacement gradients - const Tensor <2,dim> F = static_cast > (unit_symmetric_tensor ()) + grad_u_n; - J = determinant(F); - F_inv = invert(F); - B_bar = std::pow(get_J(), -2.0/3.0) * symmetrize ( F* transpose (F) ); - - // Precalculated pressure, dilatation - pressure_n = pressure; - dilatation_n = dilatation; - - // Now that all the necessary variables are set, we can update the stress tensors - // Stress update can only update the stresses once the - // dilatation has been set as p = p(d) - T_bar = material->get_T (get_J(), get_B_bar()); - T_iso = dev_P*get_T_bar(); // Note: T_iso depends on T_bar - T_vol = get_pressure()*get_J()*I; - } + // Initialise all tensors correctly + update_values (Tensor <2,dim> (), 0.0, 1.0); + } - // Displacement and strain - const double & get_dilatation(void) const {return dilatation_n;} - const double & get_J (void) const {return J;} - const Tensor <2,dim> & get_F_inv (void) const {return F_inv;} - const SymmetricTensor <2,dim> & get_B_bar (void) const {return B_bar;} + // Total Variables + void update_values (const Tensor<2, dim> & grad_u_n, + const double & pressure, + const double & dilatation) + { + // Calculated variables from displacement, displacement gradients + const Tensor <2,dim> F = static_cast > (unit_symmetric_tensor ()) + grad_u_n; + J = determinant(F); + F_inv = invert(F); + B_bar = std::pow(get_J(), -2.0/3.0) * symmetrize ( F* transpose (F) ); + + // Precalculated pressure, dilatation + pressure_n = pressure; + dilatation_n = dilatation; + + // Now that all the necessary variables are set, we can update the stress tensors + // Stress update can only update the stresses once the + // dilatation has been set as p = p(d) + T_bar = material->get_T (get_J(), get_B_bar()); + T_iso = dev_P*get_T_bar(); // Note: T_iso depends on T_bar + T_vol = get_pressure()*get_J()*I; + } + + // Displacement and strain + const double & get_dilatation(void) const {return dilatation_n;} + const double & get_J (void) const {return J;} + const Tensor <2,dim> & get_F_inv (void) const {return F_inv;} + const SymmetricTensor <2,dim> & get_B_bar (void) const {return B_bar;} - // Volumetric terms - double get_dU_dtheta (void) { + // Volumetric terms + double get_dU_dtheta (void) { return material->get_dU_dtheta(get_dilatation()); - } + } - double get_d2U_dtheta2 (void) { + double get_d2U_dtheta2 (void) { return material->get_d2U_dtheta2(get_dilatation()); - } + } - // Stress - double get_pressure(void) {return pressure_n;} - const SymmetricTensor<2, dim> & get_T_iso (void) const {return T_iso;} - const SymmetricTensor<2, dim> & get_T_vol (void) const {return T_vol;}; + // Stress + double get_pressure(void) {return pressure_n;} + const SymmetricTensor<2, dim> & get_T_iso (void) const {return T_iso;} + const SymmetricTensor<2, dim> & get_T_vol (void) const {return T_vol;}; - // Tangent matrices - SymmetricTensor <4,dim> get_C_iso(void) - { - const double & J = get_J(); - const SymmetricTensor<2, dim> & B_bar = get_B_bar(); - const SymmetricTensor<2, dim> & T_iso = get_T_iso(); + // Tangent matrices + SymmetricTensor <4,dim> get_C_iso(void) + { + const double & J = get_J(); + const SymmetricTensor<2, dim> & B_bar = get_B_bar(); + const SymmetricTensor<2, dim> & T_iso = get_T_iso(); - const SymmetricTensor <4,dim> T_iso_x_I = outer_product(T_iso, I); - const SymmetricTensor <4,dim> I_x_T_iso = outer_product(I, T_iso); - const SymmetricTensor <4,dim> CC_bar = material->get_JC (J, B_bar); + const SymmetricTensor <4,dim> T_iso_x_I = outer_product(T_iso, I); + const SymmetricTensor <4,dim> I_x_T_iso = outer_product(I, T_iso); + const SymmetricTensor <4,dim> CC_bar = material->get_JC (J, B_bar); - return 2.0/3.0*trace(get_T_bar())*dev_P - - 2.0/3.0*(T_iso_x_I + I_x_T_iso) - + dev_P*CC_bar*dev_P; - } + return 2.0/3.0*trace(get_T_bar())*dev_P + - 2.0/3.0*(T_iso_x_I + I_x_T_iso) + + dev_P*CC_bar*dev_P; + } - SymmetricTensor <4,dim> get_C_vol(void) - { - const double & p = get_pressure(); - const double & J = get_J(); - return p*J*(IxI - 2.0*II); - } + SymmetricTensor <4,dim> get_C_vol(void) + { + const double & p = get_pressure(); + const double & J = get_J(); + return p*J*(IxI - 2.0*II); + } -private: - // === MATERIAL === - Material_NH * material; - - // ==== VOLUME, DISPLACEMENT AND STRAIN VARIABLES ==== - double dilatation_n; // Current dilatation - double J; - Tensor <2,dim> F_inv; - SymmetricTensor <2,dim> B_bar; - SymmetricTensor <2,dim> E; - - // ==== STRESS VARIABLES ==== - double pressure_n; // Current pressure - SymmetricTensor<2, dim> T_bar; - SymmetricTensor<2, dim> T_iso; - SymmetricTensor<2, dim> T_vol; - const SymmetricTensor<2, dim> & get_T_bar (void) const {return T_bar;} - - // Basis tensors - static SymmetricTensor<2, dim> const I; - static SymmetricTensor<4, dim> const IxI; - static SymmetricTensor<4, dim> const II; - static SymmetricTensor<4, dim> const dev_P; -}; - -template SymmetricTensor<2,dim> const PointHistory::I -= SymmetricTensor<2,dim> (unit_symmetric_tensor ()); -template SymmetricTensor<4,dim> const PointHistory::IxI -= SymmetricTensor<4,dim> (outer_product (I, I)); -template SymmetricTensor<4,dim> const PointHistory::II -= SymmetricTensor<4,dim> (identity_tensor ()); -template SymmetricTensor<4,dim> const PointHistory::dev_P -= SymmetricTensor<4,dim> (II - 1.0/3.0*IxI); + private: + // === MATERIAL === + Material_NH * material; + + // ==== VOLUME, DISPLACEMENT AND STRAIN VARIABLES ==== + double dilatation_n; // Current dilatation + double J; + Tensor <2,dim> F_inv; + SymmetricTensor <2,dim> B_bar; + SymmetricTensor <2,dim> E; + + // ==== STRESS VARIABLES ==== + double pressure_n; // Current pressure + SymmetricTensor<2, dim> T_bar; + SymmetricTensor<2, dim> T_iso; + SymmetricTensor<2, dim> T_vol; + const SymmetricTensor<2, dim> & get_T_bar (void) const {return T_bar;} + + // Basis tensors + static SymmetricTensor<2, dim> const I; + static SymmetricTensor<4, dim> const IxI; + static SymmetricTensor<4, dim> const II; + static SymmetricTensor<4, dim> const dev_P; + }; + + template SymmetricTensor<2,dim> const PointHistory::I + = SymmetricTensor<2,dim> (unit_symmetric_tensor ()); + template SymmetricTensor<4,dim> const PointHistory::IxI + = SymmetricTensor<4,dim> (outer_product (I, I)); + template SymmetricTensor<4,dim> const PointHistory::II + = SymmetricTensor<4,dim> (identity_tensor ()); + template SymmetricTensor<4,dim> const PointHistory::dev_P + = SymmetricTensor<4,dim> (II - 1.0/3.0*IxI); // @sect3{Quasi-static quasi-incompressible finite-strain solid} -template -class Solid -{ -public: - Solid (const std::string & input_file); - virtual ~Solid (void); - void run (void); + template + class Solid + { + public: + Solid (const std::string & input_file); + virtual ~Solid (void); + void run (void); + + private: + + // === DATA STRUCTS === + + struct PerTaskData_K + { + FullMatrix cell_matrix; + std::vector local_dof_indices; + + PerTaskData_K (const unsigned int dofs_per_cell) + : + cell_matrix (dofs_per_cell, + dofs_per_cell), + local_dof_indices (dofs_per_cell) + { } + + void reset (void) { + cell_matrix = 0.0; + } + }; -private: + struct ScratchData_K + { + FEValues fe_values_ref; - // === DATA STRUCTS === + std::vector < std::vector< double > > Nx; + std::vector < std::vector< Tensor<2, dim> > > grad_Nx; + std::vector < std::vector< SymmetricTensor<2, dim> > > symm_grad_Nx; - struct PerTaskData_K - { - FullMatrix cell_matrix; - std::vector local_dof_indices; - - PerTaskData_K (const unsigned int dofs_per_cell) - : - cell_matrix (dofs_per_cell, - dofs_per_cell), - local_dof_indices (dofs_per_cell) - { } - - void reset (void) { - cell_matrix = 0.0; - } - }; - - struct ScratchData_K - { - FEValues fe_values_ref; - - std::vector < std::vector< double > > Nx; - std::vector < std::vector< Tensor<2, dim> > > grad_Nx; - std::vector < std::vector< SymmetricTensor<2, dim> > > symm_grad_Nx; - - ScratchData_K ( const FiniteElement & fe_cell, - const QGauss & qf_cell, - const UpdateFlags uf_cell) - : - fe_values_ref (fe_cell, - qf_cell, - uf_cell), - Nx (qf_cell.size(), - std::vector< double >(fe_cell.dofs_per_cell)), - grad_Nx (qf_cell.size(), - std::vector< Tensor<2, dim> >(fe_cell.dofs_per_cell)), - symm_grad_Nx (qf_cell.size(), - std::vector< SymmetricTensor<2, dim> >(fe_cell.dofs_per_cell)) - { } - - ScratchData_K ( const ScratchData_K & rhs ) : - fe_values_ref ( rhs.fe_values_ref.get_fe(), - rhs.fe_values_ref.get_quadrature(), - rhs.fe_values_ref.get_update_flags() ), - Nx (rhs.Nx), - grad_Nx (rhs.grad_Nx), - symm_grad_Nx (rhs.symm_grad_Nx) - { } - - void reset (void) { + ScratchData_K ( const FiniteElement & fe_cell, + const QGauss & qf_cell, + const UpdateFlags uf_cell) + : + fe_values_ref (fe_cell, + qf_cell, + uf_cell), + Nx (qf_cell.size(), + std::vector< double >(fe_cell.dofs_per_cell)), + grad_Nx (qf_cell.size(), + std::vector< Tensor<2, dim> >(fe_cell.dofs_per_cell)), + symm_grad_Nx (qf_cell.size(), + std::vector< SymmetricTensor<2, dim> >(fe_cell.dofs_per_cell)) + { } + + ScratchData_K ( const ScratchData_K & rhs ) : + fe_values_ref ( rhs.fe_values_ref.get_fe(), + rhs.fe_values_ref.get_quadrature(), + rhs.fe_values_ref.get_update_flags() ), + Nx (rhs.Nx), + grad_Nx (rhs.grad_Nx), + symm_grad_Nx (rhs.symm_grad_Nx) + { } + + void reset (void) { for (unsigned int q_point=0; q_point < grad_Nx.size(); ++q_point) { - for (unsigned int k=0; k < Nx.size(); ++k) { - Nx[q_point][k] = 0.0; - grad_Nx[q_point][k] = 0.0; - symm_grad_Nx[q_point][k] = 0.0; - } + for (unsigned int k=0; k < Nx.size(); ++k) { + Nx[q_point][k] = 0.0; + grad_Nx[q_point][k] = 0.0; + symm_grad_Nx[q_point][k] = 0.0; + } } - } + } - }; - - struct PerTaskData_F - { - Vector cell_rhs; - std::vector local_dof_indices; - - PerTaskData_F (const unsigned int dofs_per_cell) - : - cell_rhs (dofs_per_cell), - local_dof_indices (dofs_per_cell) - { } - - void reset (void) { cell_rhs = 0.0; } - }; + }; - struct ScratchData_F - { - FEValues fe_values_ref; - FEFaceValues fe_face_values_ref; - - std::vector < std::vector< double > > Nx; - std::vector < std::vector< SymmetricTensor<2, dim> > > symm_grad_Nx; - std::vector< Vector > rhs_values; - - // Solution data - std::vector< std::vector > > solution_grads; - - ScratchData_F ( const FiniteElement & fe_cell, - const QGauss & qf_cell, - const UpdateFlags uf_cell, - const QGauss & qf_face, - const UpdateFlags uf_face) - : - fe_values_ref (fe_cell, - qf_cell, - uf_cell), - fe_face_values_ref (fe_cell, - qf_face, - uf_face), - Nx (qf_cell.size(), - std::vector< double >(fe_cell.dofs_per_cell)), - symm_grad_Nx (qf_cell.size(), - std::vector< SymmetricTensor<2, dim> >(fe_cell.dofs_per_cell)), - rhs_values (qf_cell.size(), - Vector(dim)) - { } - - ScratchData_F ( const ScratchData_F & rhs ) - : - fe_values_ref ( rhs.fe_values_ref.get_fe(), - rhs.fe_values_ref.get_quadrature(), - rhs.fe_values_ref.get_update_flags() ), - fe_face_values_ref ( rhs.fe_face_values_ref.get_fe(), - rhs.fe_face_values_ref.get_quadrature(), - rhs.fe_face_values_ref.get_update_flags() ), - Nx (rhs.Nx), - symm_grad_Nx (rhs.symm_grad_Nx), - rhs_values (rhs.rhs_values) - { } - - void reset (void) { - for (unsigned int q_point=0; q_point < symm_grad_Nx.size(); ++q_point) { - for (unsigned int k=0; k < symm_grad_Nx[q_point].size(); ++k) { - Nx[q_point][k] = 0.0; - symm_grad_Nx[q_point][k] = 0.0; - rhs_values[q_point] = 0.0; - } - } - } + struct PerTaskData_F + { + Vector cell_rhs; + std::vector local_dof_indices; - }; + PerTaskData_F (const unsigned int dofs_per_cell) + : + cell_rhs (dofs_per_cell), + local_dof_indices (dofs_per_cell) + { } - struct PerTaskData_SC - { - FullMatrix cell_matrix; - std::vector local_dof_indices; - - // Calculation matrices (auto resized) - FullMatrix K_orig; - FullMatrix K_pu; - FullMatrix K_pt; - FullMatrix K_tt; - // Calculation matrices (manual resized) - FullMatrix K_pt_inv; - FullMatrix K_tt_inv; - FullMatrix K_con; - FullMatrix A; - FullMatrix B; - FullMatrix C; - - PerTaskData_SC (const unsigned int & dofs_per_cell, - const unsigned int & n_u, - const unsigned int & n_p, - const unsigned int & n_t) - : - cell_matrix (dofs_per_cell, - dofs_per_cell), - local_dof_indices (dofs_per_cell), - K_pt_inv (n_t, n_p), - K_tt_inv (n_t, n_t), - K_con (n_u, n_u), - A (n_t, n_u), - B (n_t, n_u), - C (n_p, n_u) - { } - - // Choose not to reset any data - // The matrix extraction and replacement tools will take care of this - void reset(void) { } - }; + void reset (void) { cell_rhs = 0.0; } + }; - // Dummy struct for TBB - struct ScratchData_SC - { - ScratchData_SC (void) { } - ScratchData_SC (const ScratchData_SC & rhs) { } - void reset (void) { } - }; + struct ScratchData_F + { + FEValues fe_values_ref; + FEFaceValues fe_face_values_ref; - // Dummy struct for TBB - struct PerTaskData_UQPH - { - PerTaskData_UQPH (void) { } - void reset(void) { } - }; + std::vector < std::vector< double > > Nx; + std::vector < std::vector< SymmetricTensor<2, dim> > > symm_grad_Nx; + std::vector< Vector > rhs_values; - struct ScratchData_UQPH - { - FEValues fe_values_ref; - std::vector< Tensor< 2, dim> > solution_grads_u_total; - std::vector solution_values_p_total; - std::vector solution_values_t_total; - const BlockVector & solution_total; + // Solution data + std::vector< std::vector > > solution_grads; - ScratchData_UQPH (const FiniteElement & fe_cell, + ScratchData_F ( const FiniteElement & fe_cell, const QGauss & qf_cell, const UpdateFlags uf_cell, - const BlockVector & solution_total) - : - fe_values_ref (fe_cell, - qf_cell, - uf_cell), - solution_grads_u_total (qf_cell.size()), - solution_values_p_total (qf_cell.size()), - solution_values_t_total (qf_cell.size()), - solution_total (solution_total) - { } - - ScratchData_UQPH (const ScratchData_UQPH & rhs) - : - fe_values_ref (rhs.fe_values_ref.get_fe(), - rhs.fe_values_ref.get_quadrature(), - rhs.fe_values_ref.get_update_flags()), - solution_grads_u_total (rhs.solution_grads_u_total), - solution_values_p_total (rhs.solution_values_p_total), - solution_values_t_total (rhs.solution_values_t_total), - solution_total (rhs.solution_total) - { } - - void reset (void) - { - // Is this necessary? Won't the call to fe_values.get_gradient overwrite this data? - for (unsigned int q=0; q < qf_cell.size(); ++q) + const QGauss & qf_face, + const UpdateFlags uf_face) + : + fe_values_ref (fe_cell, + qf_cell, + uf_cell), + fe_face_values_ref (fe_cell, + qf_face, + uf_face), + Nx (qf_cell.size(), + std::vector< double >(fe_cell.dofs_per_cell)), + symm_grad_Nx (qf_cell.size(), + std::vector< SymmetricTensor<2, dim> >(fe_cell.dofs_per_cell)), + rhs_values (qf_cell.size(), + Vector(dim)) + { } + + ScratchData_F ( const ScratchData_F & rhs ) + : + fe_values_ref ( rhs.fe_values_ref.get_fe(), + rhs.fe_values_ref.get_quadrature(), + rhs.fe_values_ref.get_update_flags() ), + fe_face_values_ref ( rhs.fe_face_values_ref.get_fe(), + rhs.fe_face_values_ref.get_quadrature(), + rhs.fe_face_values_ref.get_update_flags() ), + Nx (rhs.Nx), + symm_grad_Nx (rhs.symm_grad_Nx), + rhs_values (rhs.rhs_values) + { } + + void reset (void) { + for (unsigned int q_point=0; q_point < symm_grad_Nx.size(); ++q_point) { + for (unsigned int k=0; k < symm_grad_Nx[q_point].size(); ++k) { + Nx[q_point][k] = 0.0; + symm_grad_Nx[q_point][k] = 0.0; + rhs_values[q_point] = 0.0; + } + } + } + + }; + + struct PerTaskData_SC + { + FullMatrix cell_matrix; + std::vector local_dof_indices; + + // Calculation matrices (auto resized) + FullMatrix K_orig; + FullMatrix K_pu; + FullMatrix K_pt; + FullMatrix K_tt; + // Calculation matrices (manual resized) + FullMatrix K_pt_inv; + FullMatrix K_tt_inv; + FullMatrix K_con; + FullMatrix A; + FullMatrix B; + FullMatrix C; + + PerTaskData_SC (const unsigned int & dofs_per_cell, + const unsigned int & n_u, + const unsigned int & n_p, + const unsigned int & n_t) + : + cell_matrix (dofs_per_cell, + dofs_per_cell), + local_dof_indices (dofs_per_cell), + K_pt_inv (n_t, n_p), + K_tt_inv (n_t, n_t), + K_con (n_u, n_u), + A (n_t, n_u), + B (n_t, n_u), + C (n_p, n_u) + { } + + // Choose not to reset any data + // The matrix extraction and replacement tools will take care of this + void reset(void) { } + }; + + // Dummy struct for TBB + struct ScratchData_SC + { + ScratchData_SC (void) { } + ScratchData_SC (const ScratchData_SC & rhs) { } + void reset (void) { } + }; + + // Dummy struct for TBB + struct PerTaskData_UQPH + { + PerTaskData_UQPH (void) { } + void reset(void) { } + }; + + struct ScratchData_UQPH + { + FEValues fe_values_ref; + std::vector< Tensor< 2, dim> > solution_grads_u_total; + std::vector solution_values_p_total; + std::vector solution_values_t_total; + const BlockVector & solution_total; + + ScratchData_UQPH (const FiniteElement & fe_cell, + const QGauss & qf_cell, + const UpdateFlags uf_cell, + const BlockVector & solution_total) + : + fe_values_ref (fe_cell, + qf_cell, + uf_cell), + solution_grads_u_total (qf_cell.size()), + solution_values_p_total (qf_cell.size()), + solution_values_t_total (qf_cell.size()), + solution_total (solution_total) + { } + + ScratchData_UQPH (const ScratchData_UQPH & rhs) + : + fe_values_ref (rhs.fe_values_ref.get_fe(), + rhs.fe_values_ref.get_quadrature(), + rhs.fe_values_ref.get_update_flags()), + solution_grads_u_total (rhs.solution_grads_u_total), + solution_values_p_total (rhs.solution_values_p_total), + solution_values_t_total (rhs.solution_values_t_total), + solution_total (rhs.solution_total) + { } + + void reset (void) { - solution_grads_u_total[q] = 0.0; - solution_values_p_total[q] = 0.0; - solution_values_t_total[q] = 0.0; + // Is this necessary? Won't the call to fe_values.get_gradient overwrite this data? + for (unsigned int q=0; q < qf_cell.size(); ++q) + { + solution_grads_u_total[q] = 0.0; + solution_values_p_total[q] = 0.0; + solution_values_t_total[q] = 0.0; + } } - } - }; - - // === METHODS === - - /// \brief Print out a greeting for the user - void make_grid (void); - /// \brief Setup the Finite Element system to be solved - void system_setup (void); - void determine_component_extractors(void); - - /// \brief Assemble the system and right hand side matrices using multi-threading - void assemble_system_K (void); - void assemble_system_K_one_cell (const typename DoFHandler::active_cell_iterator & cell, - ScratchData_K & scratch, - PerTaskData_K & data); - void copy_local_to_global_K (const PerTaskData_K & data); - void assemble_system_F (void); - void assemble_system_F_one_cell (const typename DoFHandler::active_cell_iterator & cell, - ScratchData_F & scratch, - PerTaskData_F & data); - void copy_local_to_global_F (const PerTaskData_F & data); - void assemble_SC (void); - void assemble_SC_one_cell (const typename DoFHandler::active_cell_iterator & cell, - ScratchData_SC & scratch, - PerTaskData_SC & data); - void copy_local_to_global_SC (const PerTaskData_SC & data); - /// \brief Apply Dirichlet boundary values - void make_constraints (const int & it_nr, - ConstraintMatrix & constraints); - - // /// \brief Setup the quadrature point history for each cell - void setup_qph(void); - // /// \brief Update the quadrature points stress and strain values, and fibre directions - void update_qph_incremental ( const BlockVector & solution_delta ); - void update_qph_incremental_one_cell (const typename DoFHandler::active_cell_iterator & cell, - ScratchData_UQPH & scratch, - PerTaskData_UQPH & data); - void copy_local_to_global_UQPH (const PerTaskData_UQPH & data) {} - /// \brief Solve for the displacement using a Newton-Rhapson method - void solve_nonlinear_timestep (BlockVector & solution_delta); - void solve_linear_system (BlockVector & newton_update); - - /// \brief Error measurement - void get_error_res (const BlockVector & residual, BlockVector & error_res); - void get_error_update (const BlockVector & newton_update, BlockVector & error_update); - double get_error_dil (void); - - // Solution - BlockVector get_solution_total (const BlockVector & solution_delta); - - // Postprocessing - void output_results(void); - - // === ATTRIBUTES === - // Parameters - Parameters::AllParameters parameters; - - // Geometry - Triangulation triangulation; // Describes the triangulation - - // Time - Time time; - TimerOutput timer; - - // === Quadrature points === - std::vector< PointHistory > quadrature_point_history; // Quadrature point history - - // === Finite element system === - DoFHandler dof_handler_ref; // Describes the degrees of freedom - const unsigned int degree; - const FESystem fe; // Describes the global FE system - - unsigned int dofs_per_cell; // Number of degrees of freedom on each cell - const FEValuesExtractors::Vector u_fe; - const FEValuesExtractors::Scalar p_fe; - const FEValuesExtractors::Scalar t_fe; - - // Block description - static const unsigned int n_blocks = 3; - static const unsigned int n_components = dim + 2; - static const unsigned int first_u_component = 0; - static const unsigned int p_component = dim; - static const unsigned int t_component = dim + 1; - - enum {u_dof=0 , p_dof, t_dof}; - std::vector dofs_per_block; - std::vector element_indices_u; - std::vector element_indices_p; - std::vector element_indices_t; - - // === Quadrature === - QGauss qf_cell; // Cell quadrature formula - QGauss qf_face; // Face quadrature formula - unsigned int n_q_points; // Number of quadrature points in a cell - unsigned int n_q_points_f; // Number of quadrature points in a face - - // === Stiffness matrix setup ==== - ConstraintMatrix constraints; // Matrix to keep track of all constraints - BlockSparsityPattern sparsity_pattern; // Sparsity pattern for the stiffness matrix - BlockSparseMatrix tangent_matrix; // Global stiffness matrix - BlockVector residual; // Holds the residual vector - BlockVector solution_n; // Holds the solution vector: Total displacement over all time-steps -}; + }; + + // === METHODS === + + /// \brief Print out a greeting for the user + void make_grid (void); + /// \brief Setup the Finite Element system to be solved + void system_setup (void); + void determine_component_extractors(void); + + /// \brief Assemble the system and right hand side matrices using multi-threading + void assemble_system_K (void); + void assemble_system_K_one_cell (const typename DoFHandler::active_cell_iterator & cell, + ScratchData_K & scratch, + PerTaskData_K & data); + void copy_local_to_global_K (const PerTaskData_K & data); + void assemble_system_F (void); + void assemble_system_F_one_cell (const typename DoFHandler::active_cell_iterator & cell, + ScratchData_F & scratch, + PerTaskData_F & data); + void copy_local_to_global_F (const PerTaskData_F & data); + void assemble_SC (void); + void assemble_SC_one_cell (const typename DoFHandler::active_cell_iterator & cell, + ScratchData_SC & scratch, + PerTaskData_SC & data); + void copy_local_to_global_SC (const PerTaskData_SC & data); + /// \brief Apply Dirichlet boundary values + void make_constraints (const int & it_nr, + ConstraintMatrix & constraints); + + // /// \brief Setup the quadrature point history for each cell + void setup_qph(void); + // /// \brief Update the quadrature points stress and strain values, and fibre directions + void update_qph_incremental ( const BlockVector & solution_delta ); + void update_qph_incremental_one_cell (const typename DoFHandler::active_cell_iterator & cell, + ScratchData_UQPH & scratch, + PerTaskData_UQPH & data); + void copy_local_to_global_UQPH (const PerTaskData_UQPH & data) {} + /// \brief Solve for the displacement using a Newton-Rhapson method + void solve_nonlinear_timestep (BlockVector & solution_delta); + void solve_linear_system (BlockVector & newton_update); + + /// \brief Error measurement + void get_error_res (const BlockVector & residual, BlockVector & error_res); + void get_error_update (const BlockVector & newton_update, BlockVector & error_update); + double get_error_dil (void); + + // Solution + BlockVector get_solution_total (const BlockVector & solution_delta); + + // Postprocessing + void output_results(void); + + // === ATTRIBUTES === + // Parameters + Parameters::AllParameters parameters; + + // Geometry + Triangulation triangulation; // Describes the triangulation + + // Time + Time time; + TimerOutput timer; + + // === Quadrature points === + std::vector< PointHistory > quadrature_point_history; // Quadrature point history + + // === Finite element system === + DoFHandler dof_handler_ref; // Describes the degrees of freedom + const unsigned int degree; + const FESystem fe; // Describes the global FE system + + unsigned int dofs_per_cell; // Number of degrees of freedom on each cell + const FEValuesExtractors::Vector u_fe; + const FEValuesExtractors::Scalar p_fe; + const FEValuesExtractors::Scalar t_fe; + + // Block description + static const unsigned int n_blocks = 3; + static const unsigned int n_components = dim + 2; + static const unsigned int first_u_component = 0; + static const unsigned int p_component = dim; + static const unsigned int t_component = dim + 1; + + enum {u_dof=0 , p_dof, t_dof}; + std::vector dofs_per_block; + std::vector element_indices_u; + std::vector element_indices_p; + std::vector element_indices_t; + + // === Quadrature === + QGauss qf_cell; // Cell quadrature formula + QGauss qf_face; // Face quadrature formula + unsigned int n_q_points; // Number of quadrature points in a cell + unsigned int n_q_points_f; // Number of quadrature points in a face + + // === Stiffness matrix setup ==== + ConstraintMatrix constraints; // Matrix to keep track of all constraints + BlockSparsityPattern sparsity_pattern; // Sparsity pattern for the stiffness matrix + BlockSparseMatrix tangent_matrix; // Global stiffness matrix + BlockVector residual; // Holds the residual vector + BlockVector solution_n; // Holds the solution vector: Total displacement over all time-steps + }; // @sect3{Implementation of the Solid class} // @sect4{Public interface} -template -Solid::Solid (const std::string & input_file) - : - parameters (input_file), - triangulation (Triangulation::maximum_smoothing), - time (parameters.end_time, parameters.delta_t), - timer (std::cout, - TimerOutput::summary, - TimerOutput::wall_times), - dof_handler_ref (triangulation), - degree (parameters.poly_degree), - fe (FE_Q(parameters.poly_degree), dim, // displacement - FE_DGPMonomial(parameters.poly_degree-1), 1, // pressure - FE_DGPMonomial(parameters.poly_degree-1), 1), // dilatation - u_fe (first_u_component), - p_fe (p_component), - t_fe (t_component), - dofs_per_block (n_blocks), - qf_cell (parameters.quad_order), - qf_face (parameters.quad_order) -{ + template + Solid::Solid (const std::string & input_file) + : + parameters (input_file), + triangulation (Triangulation::maximum_smoothing), + time (parameters.end_time, parameters.delta_t), + timer (std::cout, + TimerOutput::summary, + TimerOutput::wall_times), + dof_handler_ref (triangulation), + degree (parameters.poly_degree), + fe (FE_Q(parameters.poly_degree), dim, // displacement + FE_DGPMonomial(parameters.poly_degree-1), 1, // pressure + FE_DGPMonomial(parameters.poly_degree-1), 1), // dilatation + u_fe (first_u_component), + p_fe (p_component), + t_fe (t_component), + dofs_per_block (n_blocks), + qf_cell (parameters.quad_order), + qf_face (parameters.quad_order) + { n_q_points = qf_cell.size(); n_q_points_f = qf_face.size(); dofs_per_cell = fe.dofs_per_cell; determine_component_extractors(); -} + } -template -Solid::~Solid (void) -{ + template + Solid::~Solid (void) + { dof_handler_ref.clear (); -} + } -template -void Solid::run (void) -{ - // Pre-processing + template + void Solid::run (void) + { + // Pre-processing make_grid (); system_setup (); output_results (); // Output initial grid position @@ -1004,129 +1007,129 @@ void Solid::run (void) solution_delta.collect_sizes (); while (time.current() <= time.end()) { - solution_delta = 0.0; + solution_delta = 0.0; - // Solve step and update total solution vector - solve_nonlinear_timestep (solution_delta); - solution_n += solution_delta; + // Solve step and update total solution vector + solve_nonlinear_timestep (solution_delta); + solution_n += solution_delta; - output_results (); - time.increment(); + output_results (); + time.increment(); } -} + } // @sect4{Solid::make_grid} -template -void Solid::make_grid (void) -{ + template + void Solid::make_grid (void) + { GridGenerator::hyper_rectangle ( triangulation, - Point (0.0, 0.0, 0.0), - Point (1.0, 1.0, 1.0), - true ); + Point (0.0, 0.0, 0.0), + Point (1.0, 1.0, 1.0), + true ); GridTools::scale (parameters.scale, triangulation); - // Need to refine at least once for the indentation problem + // Need to refine at least once for the indentation problem if (parameters.global_refinement == 0) triangulation.refine_global (1); else triangulation.refine_global (parameters.global_refinement); - // Apply different BC's to a patch on the top surface + // Apply different BC's to a patch on the top surface typename Triangulation::active_cell_iterator - cell = triangulation.begin_active(), - endc = triangulation.end(); + cell = triangulation.begin_active(), + endc = triangulation.end(); for (; cell!=endc; ++cell) - { + { if (cell->at_boundary() == true) { - for (unsigned int face=0; face < GeometryInfo::faces_per_cell; ++face) { - // Find faces on the +y surface - if ( cell->face(face)->at_boundary() == true - && cell->face(face)->center()[2] == 1.0*parameters.scale) - { - if ( cell->face(face)->center()[0] < 0.5*parameters.scale - && cell->face(face)->center()[1] < 0.5*parameters.scale) - { - cell->face(face)->set_boundary_indicator (6); // Set a new boundary id on a patch - } - } - } + for (unsigned int face=0; face < GeometryInfo::faces_per_cell; ++face) { + // Find faces on the +y surface + if ( cell->face(face)->at_boundary() == true + && cell->face(face)->center()[2] == 1.0*parameters.scale) + { + if ( cell->face(face)->center()[0] < 0.5*parameters.scale + && cell->face(face)->center()[1] < 0.5*parameters.scale) + { + cell->face(face)->set_boundary_indicator (6); // Set a new boundary id on a patch + } + } + } } - } -} + } + } // @sect4{Solid::system_setup} -template -void Solid::system_setup (void) -{ + template + void Solid::system_setup (void) + { timer.enter_subsection ("Setup system"); - // Number of components per block + // Number of components per block std::vector block_component (n_components, u_dof); // Displacement block_component[p_component] = p_dof; // Pressure block_component[t_component] = t_dof; // Dilatation - // Setup DOF handler + // Setup DOF handler dof_handler_ref.distribute_dofs (fe); DoFRenumbering::Cuthill_McKee (dof_handler_ref); DoFRenumbering::component_wise (dof_handler_ref, block_component); - // Count number of dofs per block + // Count number of dofs per block DoFTools::count_dofs_per_block (dof_handler_ref, dofs_per_block, block_component); std::cout - << "Triangulation:" - << "\n\t Number of active cells: " << triangulation.n_active_cells() - << "\n\t Number of degrees of freedom: " << dof_handler_ref.n_dofs() - << std::endl; - - // the global system matrix will have the following structure - // | K'_uu | K_up | 0 | | dU_u | | dR_u | - // K = | K_pu | K_tt^-1 | K_pt^-1 | , dU = | dU_p | , dR = | dR_p | - // | 0 | K_tp | K_tt | | dU_t | | dR_t | - // reflect this structure in the sparsity pattern + << "Triangulation:" + << "\n\t Number of active cells: " << triangulation.n_active_cells() + << "\n\t Number of degrees of freedom: " << dof_handler_ref.n_dofs() + << std::endl; + + // the global system matrix will have the following structure + // | K'_uu | K_up | 0 | | dU_u | | dR_u | + // K = | K_pu | K_tt^-1 | K_pt^-1 | , dU = | dU_p | , dR = | dR_p | + // | 0 | K_tp | K_tt | | dU_t | | dR_t | + // reflect this structure in the sparsity pattern Table<2,DoFTools::Coupling> coupling (n_components, n_components); for (unsigned int ii = 0; ii < n_components; ++ii) { - for (unsigned int jj = ii; jj < n_components; ++jj) { - if ((ii < p_component) && (jj == t_component)) { - coupling[jj][ii] = DoFTools::none; - coupling[ii][jj] = DoFTools::none; - } - else { - coupling[ii][jj] = DoFTools::always; - coupling[jj][ii] = DoFTools::always; - } - } + for (unsigned int jj = ii; jj < n_components; ++jj) { + if ((ii < p_component) && (jj == t_component)) { + coupling[jj][ii] = DoFTools::none; + coupling[ii][jj] = DoFTools::none; + } + else { + coupling[ii][jj] = DoFTools::always; + coupling[jj][ii] = DoFTools::always; + } + } } - // Setup system matrix + // Setup system matrix tangent_matrix.clear (); { - const unsigned int n_dofs_u = dofs_per_block[u_dof]; - const unsigned int n_dofs_p = dofs_per_block[p_dof]; - const unsigned int n_dofs_t = dofs_per_block[t_dof]; + const unsigned int n_dofs_u = dofs_per_block[u_dof]; + const unsigned int n_dofs_p = dofs_per_block[p_dof]; + const unsigned int n_dofs_t = dofs_per_block[t_dof]; - BlockCompressedSimpleSparsityPattern csp (n_blocks, n_blocks); + BlockCompressedSimpleSparsityPattern csp (n_blocks, n_blocks); - csp.block(u_dof,u_dof).reinit (n_dofs_u, n_dofs_u); - csp.block(u_dof,p_dof).reinit (n_dofs_u, n_dofs_p); - csp.block(u_dof,t_dof).reinit (n_dofs_u, n_dofs_t); + csp.block(u_dof,u_dof).reinit (n_dofs_u, n_dofs_u); + csp.block(u_dof,p_dof).reinit (n_dofs_u, n_dofs_p); + csp.block(u_dof,t_dof).reinit (n_dofs_u, n_dofs_t); - csp.block(p_dof,u_dof).reinit (n_dofs_p, n_dofs_u); - csp.block(p_dof,p_dof).reinit (n_dofs_p, n_dofs_p); - csp.block(p_dof,t_dof).reinit (n_dofs_p, n_dofs_t); + csp.block(p_dof,u_dof).reinit (n_dofs_p, n_dofs_u); + csp.block(p_dof,p_dof).reinit (n_dofs_p, n_dofs_p); + csp.block(p_dof,t_dof).reinit (n_dofs_p, n_dofs_t); - csp.block(t_dof,u_dof).reinit (n_dofs_t, n_dofs_u); - csp.block(t_dof,p_dof).reinit (n_dofs_t, n_dofs_p); - csp.block(t_dof,t_dof).reinit (n_dofs_t, n_dofs_t); - csp.collect_sizes(); + csp.block(t_dof,u_dof).reinit (n_dofs_t, n_dofs_u); + csp.block(t_dof,p_dof).reinit (n_dofs_t, n_dofs_p); + csp.block(t_dof,t_dof).reinit (n_dofs_t, n_dofs_t); + csp.collect_sizes(); - DoFTools::make_sparsity_pattern (dof_handler_ref, csp); - // DoFTools::make_sparsity_pattern (dof_handler_ref, csp, constraints, false); - // DoFTools::make_sparsity_pattern (dof_handler_ref, coupling, csp, constraints, false); - sparsity_pattern.copy_from (csp); + DoFTools::make_sparsity_pattern (dof_handler_ref, csp); + // DoFTools::make_sparsity_pattern (dof_handler_ref, csp, constraints, false); + // DoFTools::make_sparsity_pattern (dof_handler_ref, coupling, csp, constraints, false); + sparsity_pattern.copy_from (csp); } tangent_matrix.reinit (sparsity_pattern); - // Setup storage vectors + // Setup storage vectors residual.reinit (dofs_per_block); residual.collect_sizes (); @@ -1134,95 +1137,95 @@ void Solid::system_setup (void) solution_n.collect_sizes (); solution_n.block(t_dof) = 1.0; // Dilatation is 1 in the initial configuration - // Set up the quadrature point history + // Set up the quadrature point history setup_qph (); timer.leave_subsection(); -} + } // A way to extract subblocks from the matrix -template -void Solid::determine_component_extractors(void) -{ + template + void Solid::determine_component_extractors(void) + { element_indices_u.clear(); element_indices_p.clear(); element_indices_t.clear(); for (unsigned int k=0; k < fe.dofs_per_cell; ++k) { - // 0 = u, 1 = p, 2 = dilatation interpolation fields - const unsigned int k_group = fe.system_to_base_index(k).first.first; - if (k_group == u_dof) { - element_indices_u.push_back(k); - } - else if (k_group == p_dof) { - element_indices_p.push_back(k); - } - else if (k_group == t_dof) { - element_indices_t.push_back(k); - } - else { - Assert (k_group <= t_dof, ExcInternalError()); - } + // 0 = u, 1 = p, 2 = dilatation interpolation fields + const unsigned int k_group = fe.system_to_base_index(k).first.first; + if (k_group == u_dof) { + element_indices_u.push_back(k); + } + else if (k_group == p_dof) { + element_indices_p.push_back(k); + } + else if (k_group == t_dof) { + element_indices_t.push_back(k); + } + else { + Assert (k_group <= t_dof, ExcInternalError()); + } } -} + } // @sect4{Solid::setup_qph} -template -void Solid::setup_qph (void) -{ + template + void Solid::setup_qph (void) + { std::cout << " Setting up quadrature point data..." << std::endl; { - typename Triangulation::active_cell_iterator - cell = triangulation.begin_active(), - endc = triangulation.end(); - - unsigned int our_cells = 0; - for (; cell != endc; ++cell) { - cell->clear_user_pointer(); - ++our_cells; - } - - { - std::vector > tmp; - tmp.swap(quadrature_point_history); - } - - quadrature_point_history.resize(our_cells * n_q_points); - - unsigned int history_index = 0; - for (cell = triangulation.begin_active(); cell != endc; ++cell) { - cell->set_user_pointer(&quadrature_point_history[history_index]); - history_index += n_q_points; - } - - Assert(history_index == quadrature_point_history.size(), ExcInternalError()); + typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(), + endc = triangulation.end(); + + unsigned int our_cells = 0; + for (; cell != endc; ++cell) { + cell->clear_user_pointer(); + ++our_cells; + } + + { + std::vector > tmp; + tmp.swap(quadrature_point_history); + } + + quadrature_point_history.resize(our_cells * n_q_points); + + unsigned int history_index = 0; + for (cell = triangulation.begin_active(); cell != endc; ++cell) { + cell->set_user_pointer(&quadrature_point_history[history_index]); + history_index += n_q_points; + } + + Assert(history_index == quadrature_point_history.size(), ExcInternalError()); } - // Setup initial data + // Setup initial data typename DoFHandler::active_cell_iterator - cell = dof_handler_ref.begin_active(), - endc = dof_handler_ref.end(); + cell = dof_handler_ref.begin_active(), + endc = dof_handler_ref.end(); for (; cell != endc; ++cell) { - PointHistory* lqph = reinterpret_cast*> (cell->user_pointer()); - Assert(lqph >= &quadrature_point_history.front(), ExcInternalError()); - Assert(lqph < &quadrature_point_history.back(), ExcInternalError()); - - // Setup any initial information at displacement gauss points - for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) { - lqph[q_point].setup_lqp( parameters ); - } + PointHistory* lqph = reinterpret_cast*> (cell->user_pointer()); + Assert(lqph >= &quadrature_point_history.front(), ExcInternalError()); + Assert(lqph < &quadrature_point_history.back(), ExcInternalError()); + + // Setup any initial information at displacement gauss points + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) { + lqph[q_point].setup_lqp( parameters ); + } } -} + } // @sect4{Solid::update_qph_incremental} -template -void Solid::update_qph_incremental (const BlockVector & solution_delta) -{ + template + void Solid::update_qph_incremental (const BlockVector & solution_delta) + { timer.enter_subsection("Update QPH data"); std::cout << "Update QPH data..."<< std::endl; - // Get total solution as it stands at this update increment + // Get total solution as it stands at this update increment const BlockVector solution_total = get_solution_total(solution_delta); const UpdateFlags uf_UQPH ( update_values | update_gradients ); PerTaskData_UQPH per_task_data_UQPH; @@ -1232,21 +1235,21 @@ void Solid::update_qph_incremental (const BlockVector & solution_d solution_total); WorkStream::run ( dof_handler_ref.begin_active(), - dof_handler_ref.end(), - *this, - &Solid::update_qph_incremental_one_cell, - &Solid::copy_local_to_global_UQPH, - scratch_data_UQPH, - per_task_data_UQPH); + dof_handler_ref.end(), + *this, + &Solid::update_qph_incremental_one_cell, + &Solid::copy_local_to_global_UQPH, + scratch_data_UQPH, + per_task_data_UQPH); timer.leave_subsection(); -} + } -template -void Solid::update_qph_incremental_one_cell (const typename DoFHandler::active_cell_iterator & cell, - ScratchData_UQPH & scratch, - PerTaskData_UQPH & data) -{ + template + void Solid::update_qph_incremental_one_cell (const typename DoFHandler::active_cell_iterator & cell, + ScratchData_UQPH & scratch, + PerTaskData_UQPH & data) + { PointHistory* lqph = reinterpret_cast*> (cell->user_pointer()); Assert(lqph >= &quadrature_point_history.front(), ExcInternalError()); Assert(lqph < &quadrature_point_history.back(), ExcInternalError()); @@ -1255,35 +1258,35 @@ void Solid::update_qph_incremental_one_cell (const typename DoFHandler Assert(scratch.solution_values_p_total.size() == n_q_points, ExcInternalError()); Assert(scratch.solution_values_t_total.size() == n_q_points, ExcInternalError()); - // Find the values and gradients at quadrature points inside the current cell + // Find the values and gradients at quadrature points inside the current cell scratch.fe_values_ref.reinit(cell); scratch.fe_values_ref[u_fe].get_function_gradients (scratch.solution_total, scratch.solution_grads_u_total); scratch.fe_values_ref[p_fe].get_function_values (scratch.solution_total, scratch.solution_values_p_total); scratch.fe_values_ref[t_fe].get_function_values (scratch.solution_total,scratch. solution_values_t_total); - // === UPDATE DATA AT EACH GAUSS POINT === - // Update displacement and deformation gradient at all quadrature points + // === UPDATE DATA AT EACH GAUSS POINT === + // Update displacement and deformation gradient at all quadrature points for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) { - lqph[q_point].update_values (scratch.solution_grads_u_total [q_point], - scratch.solution_values_p_total[q_point], - scratch.solution_values_t_total[q_point]); + lqph[q_point].update_values (scratch.solution_grads_u_total [q_point], + scratch.solution_values_p_total[q_point], + scratch.solution_values_t_total[q_point]); } -} + } // @sect4{Solid::solve_nonlinear_timestep} -template -void Solid::solve_nonlinear_timestep (BlockVector & solution_delta) -{ - // timer.enter_subsection("Nonlinear solver"); + template + void Solid::solve_nonlinear_timestep (BlockVector & solution_delta) + { + // timer.enter_subsection("Nonlinear solver"); std::cout - << "Timestep " << time.get_timestep() - << std::endl; + << "Timestep " << time.get_timestep() + << std::endl; - // Newton update vector + // Newton update vector BlockVector newton_update (dofs_per_block); newton_update.collect_sizes (); - // Solution error vectors + // Solution error vectors BlockVector soln_error_res (dofs_per_block); // Holds the true residual vector BlockVector soln_error_update (dofs_per_block); // Holds the update error vector soln_error_res.collect_sizes (); @@ -1292,73 +1295,73 @@ void Solid::solve_nonlinear_timestep (BlockVector & solution_delta double res_u = 0.0, res_f = 0.0; double res_u_0 = 1.0, res_f_0 = 1.0; for (unsigned int it_nr=0; it_nr < parameters.max_iterations_NR; ++ it_nr) - { + { std::cout - << std::endl - << "Newton iteration: " << it_nr - << std::endl; + << std::endl + << "Newton iteration: " << it_nr + << std::endl; tangent_matrix = 0.0; residual = 0.0; - // Check residual + // Check residual make_constraints (it_nr, constraints); // Make boundary conditions assemble_system_F (); // Assemble RHS get_error_res(residual, soln_error_res); - // Residual scaling factors + // Residual scaling factors res_f = soln_error_res.block(u_dof).l2_norm(); if (it_nr == 0) res_f_0 = res_f; - // Check for solution convergence + // Check for solution convergence if ( it_nr > 0 - && res_u/res_u_0 <= parameters.tol_u - && res_f/res_f_0 <= parameters.tol_f) - { + && res_u/res_u_0 <= parameters.tol_u + && res_f/res_f_0 <= parameters.tol_f) + { std::cout - << std::endl - << "Solution for timestep " << time.get_timestep() - << " converged on Newton iteration " << it_nr-1 << "." - << std::endl - << "Relative displacement error: " << res_u/res_u_0 - << "\t Relative force error: " << res_f/res_f_0 - << "\t Dilatation error: " << get_error_dil() - << std::endl << std::endl; - - // timer.leave_subsection(); + << std::endl + << "Solution for timestep " << time.get_timestep() + << " converged on Newton iteration " << it_nr-1 << "." + << std::endl + << "Relative displacement error: " << res_u/res_u_0 + << "\t Relative force error: " << res_f/res_f_0 + << "\t Dilatation error: " << get_error_dil() + << std::endl << std::endl; + + // timer.leave_subsection(); return; - } + } - // No convergence -> continue with calculations - // Assemble stiffness matrix + // No convergence -> continue with calculations + // Assemble stiffness matrix assemble_system_K (); - // Do the static condensation to make K'_uu, and put K_pt^{-1} - // in the K_pt block and K_tt^{-1} in the K_pp block + // Do the static condensation to make K'_uu, and put K_pt^{-1} + // in the K_pt block and K_tt^{-1} in the K_pp block assemble_SC(); - // Do the static condensation to make K'_uu, and put K_pt^{-1} - // in the K_pt block and K_tt^{-1} in the K_pp block + // Do the static condensation to make K'_uu, and put K_pt^{-1} + // in the K_pt block and K_tt^{-1} in the K_pp block assemble_SC(); constraints.condense (tangent_matrix, residual); // Apply BC's solve_linear_system (newton_update); constraints.distribute(newton_update); // Populate the constrained DOF's with their values - // Newton update error + // Newton update error get_error_update(newton_update, soln_error_update); res_u = soln_error_update.block(u_dof).l2_norm(); - // Residual scaling factors + // Residual scaling factors if (it_nr == 0) res_u_0 = res_u; std::cout - << "Nonlinear system error: " - << std::endl << std::scientific - << " Solution update \t ||dU||: " << soln_error_update.l2_norm() - << "\t ||dU_u||: " << soln_error_update.block(u_dof).l2_norm() - << "\t ||dU_p||: " << soln_error_update.block(p_dof).l2_norm() - << "\t ||dU_t||: " << soln_error_update.block(t_dof).l2_norm() - << std::endl; + << "Nonlinear system error: " + << std::endl << std::scientific + << " Solution update \t ||dU||: " << soln_error_update.l2_norm() + << "\t ||dU_u||: " << soln_error_update.block(u_dof).l2_norm() + << "\t ||dU_p||: " << soln_error_update.block(p_dof).l2_norm() + << "\t ||dU_t||: " << soln_error_update.block(t_dof).l2_norm() + << std::endl; std::cout << std::scientific << " Residual \t ||dF||: " << soln_error_res.l2_norm() << "\t ||dR_u||: " << soln_error_res.block(u_dof).l2_norm() @@ -1371,70 +1374,70 @@ void Solid::solve_nonlinear_timestep (BlockVector & solution_delta << "\t Dilatation error: " << get_error_dil() << std::endl; - // Update and continue iterating + // Update and continue iterating solution_delta += newton_update; // Update current solution update_qph_incremental (solution_delta); // Update quadrature point information - } + } throw(ExcMessage("No convergence in nonlinear solver!")); -} + } -template -void Solid::get_error_res (const BlockVector & residual, BlockVector & error_res) -{ + template + void Solid::get_error_res (const BlockVector & residual, BlockVector & error_res) + { for (unsigned int i=0; i < dof_handler_ref.n_dofs(); ++i) - if (!constraints.is_constrained(i)) - error_res(i) = residual(i); -} + if (!constraints.is_constrained(i)) + error_res(i) = residual(i); + } -template -void Solid::get_error_update (const BlockVector & newton_update, BlockVector & error_update) -{ + template + void Solid::get_error_update (const BlockVector & newton_update, BlockVector & error_update) + { for (unsigned int i=0; i < dof_handler_ref.n_dofs(); ++i) - if (!constraints.is_constrained(i)) - error_update(i) = newton_update(i); -} + if (!constraints.is_constrained(i)) + error_update(i) = newton_update(i); + } -template -double Solid::get_error_dil (void) -{ + template + double Solid::get_error_dil (void) + { double v_e = 0.0; // Volume in current configuration double V_e = 0.0; // Volume in reference configuration FEValues fe_values_ref (fe, qf_cell, update_JxW_values); typename DoFHandler::active_cell_iterator - cell = dof_handler_ref.begin_active(), - endc = dof_handler_ref.end(); + cell = dof_handler_ref.begin_active(), + endc = dof_handler_ref.end(); for (; cell != endc; ++cell) { - fe_values_ref.reinit (cell); - PointHistory* lqph = reinterpret_cast*> (cell->user_pointer()); - Assert(lqph >= &quadrature_point_history.front(), ExcInternalError()); - Assert(lqph < &quadrature_point_history.back(), ExcInternalError()); - - for (unsigned int q_point=0; q_point < n_q_points; ++q_point) { - v_e += lqph[q_point].get_dilatation() * fe_values_ref.JxW(q_point); - V_e += fe_values_ref.JxW(q_point); - } + fe_values_ref.reinit (cell); + PointHistory* lqph = reinterpret_cast*> (cell->user_pointer()); + Assert(lqph >= &quadrature_point_history.front(), ExcInternalError()); + Assert(lqph < &quadrature_point_history.back(), ExcInternalError()); + + for (unsigned int q_point=0; q_point < n_q_points; ++q_point) { + v_e += lqph[q_point].get_dilatation() * fe_values_ref.JxW(q_point); + V_e += fe_values_ref.JxW(q_point); + } } return std::abs((v_e - V_e)/V_e); // Difference between initial and current volume -} + } // Solution (valid at any Newton step) -template -BlockVector Solid::get_solution_total (const BlockVector & solution_delta) -{ + template + BlockVector Solid::get_solution_total (const BlockVector & solution_delta) + { BlockVector solution_total (solution_n); solution_total += solution_delta; return solution_total; -} + } // @sect4{Solid::solve_linear_system} -template -void Solid::solve_linear_system (BlockVector & newton_update) -{ + template + void Solid::solve_linear_system (BlockVector & newton_update) + { std::cout << "Solve linear system..." << std::endl; BlockVector A (dofs_per_block); @@ -1442,83 +1445,83 @@ void Solid::solve_linear_system (BlockVector & newton_update) A.collect_sizes (); B.collect_sizes (); - // | K'_uu | K_up | 0 | | dU_u | | dR_u | - // K = | K_pu | K_tt^-1 | K_pt^-1 | , dU = | dU_p | , dR = | dR_p | - // | 0 | K_tp | K_tt | | dU_t | | dR_t | + // | K'_uu | K_up | 0 | | dU_u | | dR_u | + // K = | K_pu | K_tt^-1 | K_pt^-1 | , dU = | dU_p | , dR = | dR_p | + // | 0 | K_tp | K_tt | | dU_t | | dR_t | - // Solve for du + // Solve for du { - // K'uu du = Ru − Kup Ktp^-1 (Rt − Ktt Kpt^{-1} Rp) - tangent_matrix.block(p_dof, t_dof).vmult(A.block(t_dof), residual.block(p_dof)); - tangent_matrix.block(t_dof, t_dof).vmult (B.block(t_dof), A.block(t_dof)); - A.block(t_dof).equ(1.0, residual.block(t_dof), -1.0, B.block(t_dof)); - tangent_matrix.block(p_dof, t_dof).Tvmult(A.block(p_dof), A.block(t_dof)); - tangent_matrix.block(u_dof, p_dof).vmult(A.block(u_dof), A.block(p_dof)); - residual.block(u_dof) -= A.block(u_dof); + // K'uu du = Ru − Kup Ktp^-1 (Rt − Ktt Kpt^{-1} Rp) + tangent_matrix.block(p_dof, t_dof).vmult(A.block(t_dof), residual.block(p_dof)); + tangent_matrix.block(t_dof, t_dof).vmult (B.block(t_dof), A.block(t_dof)); + A.block(t_dof).equ(1.0, residual.block(t_dof), -1.0, B.block(t_dof)); + tangent_matrix.block(p_dof, t_dof).Tvmult(A.block(p_dof), A.block(t_dof)); + tangent_matrix.block(u_dof, p_dof).vmult(A.block(u_dof), A.block(p_dof)); + residual.block(u_dof) -= A.block(u_dof); - timer.enter_subsection("Linear solver"); - if (parameters.type_lin == "CG") + timer.enter_subsection("Linear solver"); + if (parameters.type_lin == "CG") { - const int solver_its = tangent_matrix.block(u_dof, u_dof).m() * parameters.max_iterations_lin; - const double tol_sol = parameters.tol_lin * residual.block(u_dof).l2_norm(); + const int solver_its = tangent_matrix.block(u_dof, u_dof).m() * parameters.max_iterations_lin; + const double tol_sol = parameters.tol_lin * residual.block(u_dof).l2_norm(); - SolverControl solver_control (solver_its , tol_sol); + SolverControl solver_control (solver_its , tol_sol); - GrowingVectorMemory < Vector > GVM; - SolverCG < Vector > solver_CG (solver_control, GVM); + GrowingVectorMemory < Vector > GVM; + SolverCG < Vector > solver_CG (solver_control, GVM); - // SSOR -> much better than Jacobi for symmetric systems - PreconditionSSOR > preconditioner; - preconditioner.initialize (tangent_matrix.block(u_dof, u_dof), parameters.ssor_relaxation); + // SSOR -> much better than Jacobi for symmetric systems + PreconditionSSOR > preconditioner; + preconditioner.initialize (tangent_matrix.block(u_dof, u_dof), parameters.ssor_relaxation); - solver_CG.solve (tangent_matrix.block(u_dof, u_dof), - newton_update.block(u_dof), - residual.block(u_dof), - preconditioner); + solver_CG.solve (tangent_matrix.block(u_dof, u_dof), + newton_update.block(u_dof), + residual.block(u_dof), + preconditioner); - std::cout - << "\t Iterations: " << solver_control.last_step() - << "\n\t Residual: " << solver_control.last_value() - << std::endl; + std::cout + << "\t Iterations: " << solver_control.last_step() + << "\n\t Residual: " << solver_control.last_value() + << std::endl; } - else if (parameters.type_lin == "Direct") + else if (parameters.type_lin == "Direct") { - SparseDirectUMFPACK A_direct; - A_direct.initialize(tangent_matrix.block(u_dof, u_dof)); - A_direct.vmult (newton_update.block(u_dof), - residual.block(u_dof)); + SparseDirectUMFPACK A_direct; + A_direct.initialize(tangent_matrix.block(u_dof, u_dof)); + A_direct.vmult (newton_update.block(u_dof), + residual.block(u_dof)); } - else throw (ExcMessage("Linear solver type not implemented")); - timer.leave_subsection(); + else throw (ExcMessage("Linear solver type not implemented")); + timer.leave_subsection(); } timer.enter_subsection("Linear solver postprocessing"); - // Postprocess for dp + // Postprocess for dp { - // dp = Ktp^{-1} ( Rt − Ktt Kpt^{-1} (Rp − Kpu du) ) - tangent_matrix.block(p_dof, u_dof).vmult (A.block(p_dof), newton_update.block(u_dof)); - B.block(p_dof).equ(1.0, residual.block(p_dof), -1.0, A.block(p_dof)); - tangent_matrix.block(p_dof, t_dof).vmult(A.block(t_dof), B.block(p_dof)); - tangent_matrix.block(t_dof, t_dof).vmult(B.block(t_dof), A.block(t_dof)); - A.block(t_dof).equ (1.0, residual.block(t_dof), -1.0, B.block(t_dof)); - tangent_matrix.block(p_dof, t_dof).Tvmult (newton_update.block(p_dof), A.block(t_dof)); + // dp = Ktp^{-1} ( Rt − Ktt Kpt^{-1} (Rp − Kpu du) ) + tangent_matrix.block(p_dof, u_dof).vmult (A.block(p_dof), newton_update.block(u_dof)); + B.block(p_dof).equ(1.0, residual.block(p_dof), -1.0, A.block(p_dof)); + tangent_matrix.block(p_dof, t_dof).vmult(A.block(t_dof), B.block(p_dof)); + tangent_matrix.block(t_dof, t_dof).vmult(B.block(t_dof), A.block(t_dof)); + A.block(t_dof).equ (1.0, residual.block(t_dof), -1.0, B.block(t_dof)); + tangent_matrix.block(p_dof, t_dof).Tvmult (newton_update.block(p_dof), A.block(t_dof)); } - // Postprocess for dt + // Postprocess for dt { - // dt = Ktt^{-1} (Rt − Ktp dp) - tangent_matrix.block(t_dof, p_dof).vmult (A.block(t_dof), newton_update.block(p_dof)); - residual.block(t_dof) -= A.block(t_dof); - tangent_matrix.block(p_dof, p_dof).vmult (newton_update.block(t_dof), residual.block(t_dof)); + // dt = Ktt^{-1} (Rt − Ktp dp) + tangent_matrix.block(t_dof, p_dof).vmult (A.block(t_dof), newton_update.block(p_dof)); + residual.block(t_dof) -= A.block(t_dof); + tangent_matrix.block(p_dof, p_dof).vmult (newton_update.block(t_dof), residual.block(t_dof)); } timer.leave_subsection(); -} + } // @sect4{Solid::assemble_system_K} -template -void Solid::assemble_system_K (void) -{ + template + void Solid::assemble_system_K (void) + { timer.enter_subsection("Assemble system matrix"); std::cout << "Assemble system matrix..."<< std::endl; @@ -1530,124 +1533,124 @@ void Solid::assemble_system_K (void) ScratchData_K scratch_data (fe, qf_cell, uf_cell); WorkStream::run ( dof_handler_ref.begin_active(), - dof_handler_ref.end(), - *this, - &Solid::assemble_system_K_one_cell, - &Solid::copy_local_to_global_K, - scratch_data, - per_task_data); + dof_handler_ref.end(), + *this, + &Solid::assemble_system_K_one_cell, + &Solid::copy_local_to_global_K, + scratch_data, + per_task_data); timer.leave_subsection(); -} + } -template -void Solid::copy_local_to_global_K (const PerTaskData_K & data) -{ - // Add the local contribution to the system matrix + template + void Solid::copy_local_to_global_K (const PerTaskData_K & data) + { + // Add the local contribution to the system matrix for (unsigned int i=0; i -void Solid::assemble_system_K_one_cell (const typename DoFHandler::active_cell_iterator & cell, - ScratchData_K & scratch, - PerTaskData_K & data) -{ + for (unsigned int j=0; j + void Solid::assemble_system_K_one_cell (const typename DoFHandler::active_cell_iterator & cell, + ScratchData_K & scratch, + PerTaskData_K & data) + { data.reset(); // Reset data in the PerTaskData_K storage unit scratch.reset(); // Reset data in the Scratch storage unit scratch.fe_values_ref.reinit (cell); cell->get_dof_indices (data.local_dof_indices); // Find out which global numbers the degrees of freedom on this cell have PointHistory *lqph = reinterpret_cast*>(cell->user_pointer()); - // Set up cell shape function gradients + // Set up cell shape function gradients static const SymmetricTensor<2, dim> I = unit_symmetric_tensor (); for (unsigned int q_point=0; q_point < n_q_points; ++q_point) { - const Tensor<2, dim> F_inv = lqph[q_point].get_F_inv(); + const Tensor<2, dim> F_inv = lqph[q_point].get_F_inv(); - for (unsigned int k=0; k< dofs_per_cell; ++k) { - const unsigned int k_group = fe.system_to_base_index(k).first.first; + for (unsigned int k=0; k< dofs_per_cell; ++k) { + const unsigned int k_group = fe.system_to_base_index(k).first.first; - if (k_group == u_dof) { - scratch.grad_Nx[q_point][k] = scratch.fe_values_ref[u_fe].gradient(k, q_point) * F_inv; - scratch.symm_grad_Nx[q_point][k] = symmetrize(scratch.grad_Nx[q_point][k]); - } - else if (k_group == p_dof) { - scratch.Nx[q_point][k] = scratch.fe_values_ref[p_fe].value(k, q_point); - } - else if (k_group == t_dof) { - scratch.Nx[q_point][k] = scratch.fe_values_ref[t_fe].value(k, q_point); - } - else { - Assert (k_group <= t_dof, ExcInternalError()); - } + if (k_group == u_dof) { + scratch.grad_Nx[q_point][k] = scratch.fe_values_ref[u_fe].gradient(k, q_point) * F_inv; + scratch.symm_grad_Nx[q_point][k] = symmetrize(scratch.grad_Nx[q_point][k]); + } + else if (k_group == p_dof) { + scratch.Nx[q_point][k] = scratch.fe_values_ref[p_fe].value(k, q_point); + } + else if (k_group == t_dof) { + scratch.Nx[q_point][k] = scratch.fe_values_ref[t_fe].value(k, q_point); } + else { + Assert (k_group <= t_dof, ExcInternalError()); + } + } } - // Build cell stiffness matrix - // Global and local system matrices are symmetric - // => Take advantage of this: Build only the lower half of the local matrix - // Only assemble 1/2 of the K_uu, K_pp = 0, K_tt blocks and the whole K_pt, K_ut, K_up blocks + // Build cell stiffness matrix + // Global and local system matrices are symmetric + // => Take advantage of this: Build only the lower half of the local matrix + // Only assemble 1/2 of the K_uu, K_pp = 0, K_tt blocks and the whole K_pt, K_ut, K_up blocks for (unsigned int q_point=0; q_point < n_q_points; ++q_point) { - const Tensor <2,dim> T = static_cast < Tensor<2, dim> > (lqph[q_point].get_T_iso() + lqph[q_point].get_T_vol()); - const SymmetricTensor <4,dim> C = lqph[q_point].get_C_iso() + lqph[q_point].get_C_vol(); - const double C_v = lqph[q_point].get_d2U_dtheta2(); - const double J = lqph[q_point].get_J(); - - const std::vector & N = scratch.Nx[q_point]; - const std::vector< SymmetricTensor <2,dim> > & symm_B = scratch.symm_grad_Nx[q_point]; - const std::vector< Tensor <2,dim> > & B = scratch.grad_Nx[q_point]; - const double & JxW = scratch.fe_values_ref.JxW(q_point); - - for (unsigned int i=0; i < dofs_per_cell; ++i) { - - const unsigned int component_i = fe.system_to_component_index(i).first; - const unsigned int i_group = fe.system_to_base_index(i).first.first; - - // Only assemble the lower diagonal part of the local matrix - for (unsigned int j=0; j <= i; ++j) { - - const unsigned int component_j = fe.system_to_component_index(j).first; - const unsigned int j_group = fe.system_to_base_index(j).first.first; - - if ( (i_group == j_group) && (i_group == u_dof ) ) { - data.cell_matrix(i,j) - += ( symm_B[i] * C * symm_B[j] // Material stiffness - + ( component_i == component_j ? - B[i][component_i] * T * B[j][component_j] : - 0.0 ) // Geometric stiffness. Only add this along local diagonals - ) * JxW; // K_uu - } - else if ( (i_group == p_dof) && (j_group == u_dof) ) { - data.cell_matrix(i,j) += N[i]*J*(symm_B[j]*I)*JxW; // K_pu - } - else if ( (i_group == t_dof) && (j_group == p_dof) ) { - data.cell_matrix(i,j) -= N[i]*N[j]*JxW; // K_tp - } - else if ( (i_group == j_group) && (i_group == t_dof) ) { - data.cell_matrix(i,j) += N[i]*C_v*N[j]*JxW; // K_tt - } - else Assert ((i_group <= t_dof) && (j_group <= t_dof), ExcInternalError()); - } // END j LOOP - } // END i LOOP + const Tensor <2,dim> T = static_cast < Tensor<2, dim> > (lqph[q_point].get_T_iso() + lqph[q_point].get_T_vol()); + const SymmetricTensor <4,dim> C = lqph[q_point].get_C_iso() + lqph[q_point].get_C_vol(); + const double C_v = lqph[q_point].get_d2U_dtheta2(); + const double J = lqph[q_point].get_J(); + + const std::vector & N = scratch.Nx[q_point]; + const std::vector< SymmetricTensor <2,dim> > & symm_B = scratch.symm_grad_Nx[q_point]; + const std::vector< Tensor <2,dim> > & B = scratch.grad_Nx[q_point]; + const double & JxW = scratch.fe_values_ref.JxW(q_point); + + for (unsigned int i=0; i < dofs_per_cell; ++i) { + + const unsigned int component_i = fe.system_to_component_index(i).first; + const unsigned int i_group = fe.system_to_base_index(i).first.first; + + // Only assemble the lower diagonal part of the local matrix + for (unsigned int j=0; j <= i; ++j) { + + const unsigned int component_j = fe.system_to_component_index(j).first; + const unsigned int j_group = fe.system_to_base_index(j).first.first; + + if ( (i_group == j_group) && (i_group == u_dof ) ) { + data.cell_matrix(i,j) + += ( symm_B[i] * C * symm_B[j] // Material stiffness + + ( component_i == component_j ? + B[i][component_i] * T * B[j][component_j] : + 0.0 ) // Geometric stiffness. Only add this along local diagonals + ) * JxW; // K_uu + } + else if ( (i_group == p_dof) && (j_group == u_dof) ) { + data.cell_matrix(i,j) += N[i]*J*(symm_B[j]*I)*JxW; // K_pu + } + else if ( (i_group == t_dof) && (j_group == p_dof) ) { + data.cell_matrix(i,j) -= N[i]*N[j]*JxW; // K_tp + } + else if ( (i_group == j_group) && (i_group == t_dof) ) { + data.cell_matrix(i,j) += N[i]*C_v*N[j]*JxW; // K_tt + } + else Assert ((i_group <= t_dof) && (j_group <= t_dof), ExcInternalError()); + } // END j LOOP + } // END i LOOP } // END q_point LOOP - // Global and local system matrices are symmetric - // => Copy the upper half of the local matrix in the bottom half of the local matrix + // Global and local system matrices are symmetric + // => Copy the upper half of the local matrix in the bottom half of the local matrix for (unsigned int i=0; i -void Solid::assemble_system_F (void) -{ + template + void Solid::assemble_system_F (void) + { timer.enter_subsection("Assemble system RHS"); std::cout << "Assemble system RHS..."<< std::endl; @@ -1664,132 +1667,132 @@ void Solid::assemble_system_F (void) uf_face); WorkStream::run ( dof_handler_ref.begin_active(), - dof_handler_ref.end(), - *this, - &Solid::assemble_system_F_one_cell, - &Solid::copy_local_to_global_F, - scratch_data, - per_task_data ); + dof_handler_ref.end(), + *this, + &Solid::assemble_system_F_one_cell, + &Solid::copy_local_to_global_F, + scratch_data, + per_task_data ); timer.leave_subsection(); -} + } -template -void Solid::copy_local_to_global_F (const PerTaskData_F & data) -{ - // Add the local contribution to the system RHS vector + template + void Solid::copy_local_to_global_F (const PerTaskData_F & data) + { + // Add the local contribution to the system RHS vector for (unsigned int i=0; i -void Solid::assemble_system_F_one_cell (const typename DoFHandler::active_cell_iterator & cell, - ScratchData_F & scratch, - PerTaskData_F & data) -{ + template + void Solid::assemble_system_F_one_cell (const typename DoFHandler::active_cell_iterator & cell, + ScratchData_F & scratch, + PerTaskData_F & data) + { data.reset(); // Reset data in the PerTaskData_K storage unit scratch.reset(); // Reset data in the ScratchData_F storage unit scratch.fe_values_ref.reinit (cell); cell->get_dof_indices (data.local_dof_indices); // Find out which global numbers the degrees of freedom on this cell have PointHistory *lqph = reinterpret_cast*>(cell->user_pointer()); - // Precompute some data + // Precompute some data for (unsigned int q_point=0; q_point < n_q_points; ++q_point) { - const Tensor<2, dim> F_inv = lqph[q_point].get_F_inv(); + const Tensor<2, dim> F_inv = lqph[q_point].get_F_inv(); - for (unsigned int k=0; k T = lqph[q_point].get_T_iso() + lqph[q_point].get_T_vol(); - const double J = lqph[q_point].get_J(); - const double D = lqph[q_point].get_dilatation(); - const double p = lqph[q_point].get_pressure(); - const double p_star = lqph[q_point].get_dU_dtheta(); + const SymmetricTensor <2,dim> T = lqph[q_point].get_T_iso() + lqph[q_point].get_T_vol(); + const double J = lqph[q_point].get_J(); + const double D = lqph[q_point].get_dilatation(); + const double p = lqph[q_point].get_pressure(); + const double p_star = lqph[q_point].get_dU_dtheta(); - const std::vector< double > & N = scratch.Nx[q_point]; - const std::vector< SymmetricTensor <2,dim> > & symm_B = scratch.symm_grad_Nx[q_point]; - const double JxW = scratch.fe_values_ref.JxW(q_point); + const std::vector< double > & N = scratch.Nx[q_point]; + const std::vector< SymmetricTensor <2,dim> > & symm_B = scratch.symm_grad_Nx[q_point]; + const double JxW = scratch.fe_values_ref.JxW(q_point); - for (unsigned int i=0; iat_boundary() == true) - { + { static const Tensor <2, dim> I = static_cast < Tensor <2, dim> > ( unit_symmetric_tensor () ); for (unsigned int face=0; face < GeometryInfo::faces_per_cell; ++face) - { + { if ( cell->face(face)->at_boundary() == true && cell->face(face)->boundary_indicator() == 6 ) - { + { scratch.fe_face_values_ref.reinit (cell, face); for (unsigned int f_q_point=0; f_q_point < n_q_points_f; ++f_q_point) - { + { const Tensor <1, dim> & N = scratch.fe_face_values_ref.normal_vector(f_q_point); - // Traction in reference configuration - // t_0 = p*N + // Traction in reference configuration + // t_0 = p*N static const double p0 = -4.0/(parameters.scale*parameters.scale); // Reference pressure of 4 Pa const double time_ramp = (time.current() / time.end()); // Linearly ramp up the pressure with time const double pressure = p0 * parameters.p_p0 * time_ramp; const Tensor <1,dim> traction = pressure * N; for (unsigned int i=0; i < dofs_per_cell; ++i) { - // Determine the dimensional component that matches the dof component (i.e. i % dim) - const unsigned int i_group = fe.system_to_base_index(i).first.first; - - if (i_group == u_dof) { - const unsigned int component_i = fe.system_to_component_index(i).first; - const double & Ni = scratch.fe_face_values_ref.shape_value(i,f_q_point); - const double & JxW = scratch.fe_face_values_ref.JxW(f_q_point); - - // Add traction vector contribution to the local RHS vector (displacement dofs only) - data.cell_rhs(i) += (Ni * traction[component_i]) // Contribution from external forces - * JxW; - } + // Determine the dimensional component that matches the dof component (i.e. i % dim) + const unsigned int i_group = fe.system_to_base_index(i).first.first; + + if (i_group == u_dof) { + const unsigned int component_i = fe.system_to_component_index(i).first; + const double & Ni = scratch.fe_face_values_ref.shape_value(i,f_q_point); + const double & JxW = scratch.fe_face_values_ref.JxW(f_q_point); + + // Add traction vector contribution to the local RHS vector (displacement dofs only) + data.cell_rhs(i) += (Ni * traction[component_i]) // Contribution from external forces + * JxW; + } } // END i LOOP - } // END face q_point LOOP - } // END at boundary check LOOP + } // END face q_point LOOP + } // END at boundary check LOOP - } // END face LOOP - } -} + } // END face LOOP + } + } // @sect4{Solid::assemble_system_SC} -template -void Solid::assemble_SC (void) -{ + template + void Solid::assemble_SC (void) + { timer.enter_subsection("Perform static condensation"); PerTaskData_SC per_task_data (dofs_per_cell, @@ -1799,75 +1802,75 @@ void Solid::assemble_SC (void) ScratchData_SC scratch_data; WorkStream::run ( dof_handler_ref.begin_active(), - dof_handler_ref.end(), - *this, - &Solid::assemble_SC_one_cell, - &Solid::copy_local_to_global_SC, - scratch_data, - per_task_data ); + dof_handler_ref.end(), + *this, + &Solid::assemble_SC_one_cell, + &Solid::copy_local_to_global_SC, + scratch_data, + per_task_data ); timer.leave_subsection(); -} + } -template -void Solid::copy_local_to_global_SC (const PerTaskData_SC & data) -{ - // Add the local contribution to the system matrix + template + void Solid::copy_local_to_global_SC (const PerTaskData_SC & data) + { + // Add the local contribution to the system matrix for (unsigned int i=0; i -void Solid::assemble_SC_one_cell (const typename DoFHandler::active_cell_iterator & cell, - ScratchData_SC & scratch, - PerTaskData_SC & data) -{ + for (unsigned int j=0; j + void Solid::assemble_SC_one_cell (const typename DoFHandler::active_cell_iterator & cell, + ScratchData_SC & scratch, + PerTaskData_SC & data) + { data.reset(); scratch.reset(); cell->get_dof_indices (data.local_dof_indices); // Find out which global numbers the degrees of freedom on this cell have - // The local stifness matrix K_e is: - // | K_uu | K_up | 0 | - // | K_pu | 0 | K_pt | - // | 0 | K_tp | K_tt | - // - // We are going to exploit the zeros for post-processing as: - // | K'_uu | K_up | 0 | - // | K_pu | K_tt^-1 | K_pt^-1 | - // | 0 | K_tp | K_tt | - // with K'_uu = K_uu + Kup Ktp^{-1} Ktt Kpt^{-1} Kpu - - // NOTE: - // GLOBAL Data already exists in the K_uu, K_pt, K_tp subblocks - // - // For the K_uu block in particular, this means that contributions have been - // added from the surrounding cells, so we need to be careful when we manipulate this block. - // We can't just erase the subblocks and - // Additionally the copy_local_to_global operation is a "+=" operation -> need to take this - // into account - // - // So the intermediate matrix that we need to get from what we have in K_uu and what we - // are actually wanting is: - // | K'_uu - K_uu | 0 | 0 | - // | 0 | K_tt^-1 | K_pt^-1 - K_pt | - // | 0 | 0 | 0 | - // - // Strategy to get the subblocks we want: - // K'_uu: Since we don't have access to K_uu^h, but we know its contribution is added to the global - // K_uu matrix, we just want to add the element wise static-condensation - // K'_uu^h = K_uu^h + K_up^h K_tp^{-1}^h K_tt^h K_pt^{-1}^h K_pu^h - // Since we already have K_uu^h in the system matrix, we just need to do the following - // K'_uu^h == (K_uu^h += K_up^h K_tp^{-1}^h K_tt^h K_pt^{-1}^h K_pu^h) - // K_pt^-1: Similarly, K_pt exists in the subblock. Since the copy operation is a += operation, we need - // to subtract the existing K_pt submatrix in addition to "adding" that which we wish to - // replace it with. - // K_tp^-1: Same as above - // K_tt^-1: Nothing exists in the original K_pp subblock, so we can just add this contribution as is. - - // Extract element data from the system matrix + // The local stifness matrix K_e is: + // | K_uu | K_up | 0 | + // | K_pu | 0 | K_pt | + // | 0 | K_tp | K_tt | + // + // We are going to exploit the zeros for post-processing as: + // | K'_uu | K_up | 0 | + // | K_pu | K_tt^-1 | K_pt^-1 | + // | 0 | K_tp | K_tt | + // with K'_uu = K_uu + Kup Ktp^{-1} Ktt Kpt^{-1} Kpu + + // NOTE: + // GLOBAL Data already exists in the K_uu, K_pt, K_tp subblocks + // + // For the K_uu block in particular, this means that contributions have been + // added from the surrounding cells, so we need to be careful when we manipulate this block. + // We can't just erase the subblocks and + // Additionally the copy_local_to_global operation is a "+=" operation -> need to take this + // into account + // + // So the intermediate matrix that we need to get from what we have in K_uu and what we + // are actually wanting is: + // | K'_uu - K_uu | 0 | 0 | + // | 0 | K_tt^-1 | K_pt^-1 - K_pt | + // | 0 | 0 | 0 | + // + // Strategy to get the subblocks we want: + // K'_uu: Since we don't have access to K_uu^h, but we know its contribution is added to the global + // K_uu matrix, we just want to add the element wise static-condensation + // K'_uu^h = K_uu^h + K_up^h K_tp^{-1}^h K_tt^h K_pt^{-1}^h K_pu^h + // Since we already have K_uu^h in the system matrix, we just need to do the following + // K'_uu^h == (K_uu^h += K_up^h K_tp^{-1}^h K_tt^h K_pt^{-1}^h K_pu^h) + // K_pt^-1: Similarly, K_pt exists in the subblock. Since the copy operation is a += operation, we need + // to subtract the existing K_pt submatrix in addition to "adding" that which we wish to + // replace it with. + // K_tp^-1: Same as above + // K_tt^-1: Nothing exists in the original K_pp subblock, so we can just add this contribution as is. + + // Extract element data from the system matrix AdditionalTools::extract_submatrix(data.local_dof_indices, data.local_dof_indices, @@ -1886,7 +1889,7 @@ void Solid::assemble_SC_one_cell (const typename DoFHandler::active_ce data.K_orig, data.K_tt); - // Place K_pt^-1 in the K_pt block + // Place K_pt^-1 in the K_pt block data.K_pt_inv.invert(data.K_pt); data.K_pt_inv.add (-1.0, data.K_pt); AdditionalTools::replace_submatrix(element_indices_p, @@ -1894,14 +1897,14 @@ void Solid::assemble_SC_one_cell (const typename DoFHandler::active_ce data.K_pt_inv, data.cell_matrix); - // Place K_tt^-1 in the K_pp block + // Place K_tt^-1 in the K_pp block data.K_tt_inv.invert(data.K_tt); AdditionalTools::replace_submatrix(element_indices_p, element_indices_p, data.K_tt_inv, data.cell_matrix); - // Make condensation terms to add to the K_uu block + // Make condensation terms to add to the K_uu block data.K_pt_inv.mmult(data.A, data.K_pu); data.K_tt.mmult(data.B, data.A); data.K_pt_inv.Tmmult(data.C, data.B); // Symmetric matrix @@ -1910,97 +1913,97 @@ void Solid::assemble_SC_one_cell (const typename DoFHandler::active_ce element_indices_u, data.K_con, data.cell_matrix); -} + } // @sect4{Solid::make_constraints} -template -void Solid::make_constraints (const int & it_nr, - ConstraintMatrix & constraints) -{ + template + void Solid::make_constraints (const int & it_nr, + ConstraintMatrix & constraints) + { std::cout << "Make constraints..."<< std::endl; constraints.clear(); const bool apply_dirichlet_bc = (it_nr == 0); - // Boundary conditions: - // b_id 0: -x face: Zero x-component of displacement : Symmetry plane - // b_id 2: -y face: Zero y-component of displacement : Symmetry plane - // b_id 4: -z face: Zero z-component of displacement : Symmetry plane + // Boundary conditions: + // b_id 0: -x face: Zero x-component of displacement : Symmetry plane + // b_id 2: -y face: Zero y-component of displacement : Symmetry plane + // b_id 4: -z face: Zero z-component of displacement : Symmetry plane - // b_id 5: +z face: Zero x-component and Zero y-component - // b_id 6: Applied pressure face: Zero x-component and Zero y-component - // b_id 1: +x face: Traction free - // b_id 3: +y face: Traction free + // b_id 5: +z face: Zero x-component and Zero y-component + // b_id 6: Applied pressure face: Zero x-component and Zero y-component + // b_id 1: +x face: Traction free + // b_id 3: +y face: Traction free { - const int boundary_id = 0; + const int boundary_id = 0; - std::vector< bool > components (n_components, false); - components[0] = true; + std::vector< bool > components (n_components, false); + components[0] = true; - if (apply_dirichlet_bc == true) { - VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); - } - else { - VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); - } + if (apply_dirichlet_bc == true) { + VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); + } + else { + VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); + } } { - const int boundary_id = 2; + const int boundary_id = 2; - std::vector< bool > components (n_components, false); - components[1] = true; + std::vector< bool > components (n_components, false); + components[1] = true; - if (apply_dirichlet_bc == true) { - VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); - } - else { - VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); - } + if (apply_dirichlet_bc == true) { + VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); + } + else { + VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); + } } { - const int boundary_id = 4; - std::vector< bool > components (n_components, false); - components[2] = true; - - if (apply_dirichlet_bc == true) { - VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); - } - else { - VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); - } + const int boundary_id = 4; + std::vector< bool > components (n_components, false); + components[2] = true; + + if (apply_dirichlet_bc == true) { + VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); + } + else { + VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); + } } { - const int boundary_id = 5; - std::vector< bool > components (n_components, true); - components[2] = false; - - if (apply_dirichlet_bc == true) { - VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); - } - else { - VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); - } + const int boundary_id = 5; + std::vector< bool > components (n_components, true); + components[2] = false; + + if (apply_dirichlet_bc == true) { + VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); + } + else { + VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); + } } { - const int boundary_id = 6; - std::vector< bool > components (n_components, true); - components[2] = false; - - if (apply_dirichlet_bc == true) { - VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); - } - else { - VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); - } + const int boundary_id = 6; + std::vector< bool > components (n_components, true); + components[2] = false; + + if (apply_dirichlet_bc == true) { + VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); + } + else { + VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); + } } constraints.close(); -} + } // @sect4{Solid::output_results} -template -void Solid::output_results(void) -{ + template + void Solid::output_results(void) + { DataOut data_out; std::vector data_component_interpretation (dim, DataComponentInterpretation::component_is_part_of_vector); @@ -2015,8 +2018,8 @@ void Solid::output_results(void) data_out.add_data_vector (solution_n, solution_name, DataOut::type_dof_data, data_component_interpretation); - // MappingQEulerian q_mapping (degree, solution_n.block(u_dof), dof_handler_ref); - // MappingQEulerian q_mapping (degree, solution_n, dof_handler_ref); + // MappingQEulerian q_mapping (degree, solution_n.block(u_dof), dof_handler_ref); + // MappingQEulerian q_mapping (degree, solution_n, dof_handler_ref); Vector soln; soln.reinit(solution_n.size()); for (unsigned int i=0; i < soln.size(); ++i) soln(i) = solution_n(i); @@ -2030,43 +2033,48 @@ void Solid::output_results(void) std::ofstream output (filename.str().c_str()); data_out.write_vtk (output); + } } + // @sect3{Main function} int main () { - try + try { - deallog.depth_console (0); + using namespace dealii; + using namespace Step44; - Solid<3> solid_3d ("parameters.prm"); - solid_3d.run(); + deallog.depth_console (0); + + Solid<3> solid_3d ("parameters.prm"); + solid_3d.run(); } - catch (std::exception &exc) + catch (std::exception &exc) { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Exception on processing: " << std::endl - << exc.what() << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; - return 1; + return 1; } - catch (...) + catch (...) { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; } - return 0; + return 0; }