From: mcbride Date: Thu, 16 Feb 2012 12:45:08 +0000 (+0000) Subject: more small changes X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c77f27d729c90affd9cec6cc7db93d2603f7f963;p=dealii-svn.git more small changes git-svn-id: https://svn.dealii.org/trunk@25098 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-44/doc/intro.dox b/deal.II/examples/step-44/doc/intro.dox index 4a6c347b0a..ad84d560b9 100644 --- a/deal.II/examples/step-44/doc/intro.dox +++ b/deal.II/examples/step-44/doc/intro.dox @@ -201,7 +201,7 @@ If the Helmholtz free energy depends on the left Cauchy-Green tensor $\mathbf{b} Following the multiplicative decomposition of the deformation gradient, the Helmholtz free energy can be decomposed as @f[ - \Psi(\mathbf{b}) = \Psi_{\text{vol}}(\mathbf{J}) + \Psi_{\text{iso}}(\overline{\mathbf{b}}) \, . + \Psi(\mathbf{b}) = \Psi_{\text{vol}}(J) + \Psi_{\text{iso}}(\overline{\mathbf{b}}) \, . @f] Similarly, the Kirchhoff stress can be decomposed into volumetric and isochoric parts as $\boldsymbol{\tau} = \boldsymbol{\tau}_{\text{vol}} + \boldsymbol{\tau}_{\text{iso}}$ where: @f{align*} @@ -216,14 +216,21 @@ Similarly, the Kirchhoff stress can be decomposed into volumetric and isochoric &= \underbrace{( \mathcal{I} - \dfrac{1}{3} \mathbf{I} \otimes \mathbf{I})}_{\mathbb{P}} : \overline{\boldsymbol{\tau}} \, , @f} where -$p = - 1/3 \textrm{tr} \boldsymbol{\sigma} = - 1/3 J^{-1} \textrm{tr} \boldsymbol{\tau}$ -is the hydrostatic pressure and $\mathbb{P}$ is the projection tensor which provides the deviatoric operator in the Eulerian setting. +$p := \dfrac{\partial \Psi_{\text{vol}}(J)}{\partial J}$ is the pressure response. +$\mathbb{P}$ is the projection tensor which provides the deviatoric operator in the Eulerian setting. The fictitious Cauchy stress tensor $\overline{\boldsymbol{\tau}}$ is defined by @f[ \overline{\boldsymbol{\tau}} := 2 \overline{\mathbf{b}} \dfrac{\partial \Psi_{\textrm{iso}}(\overline{\mathbf{b}})}{\partial \overline{\mathbf{b}}} \, . @f] + +@note The pressure response as defined above differs from the widely-used definition of the +pressure in solid mechanics as +$p = - 1/3 \textrm{tr} \boldsymbol{\sigma} = - 1/3 J^{-1} \textrm{tr} \boldsymbol{\tau}$. +Here $p$ is the hydrostatic pressure. +We make use of the pressure response throughout this tut (although we refer to it as the pressure). +

Neo-Hookean materials

The Helmholtz free energy corresponding to a compressible neo-Hookean material is given by diff --git a/deal.II/examples/step-44/step-44.cc b/deal.II/examples/step-44/step-44.cc index 277e32c4f6..bffb02b489 100644 --- a/deal.II/examples/step-44/step-44.cc +++ b/deal.II/examples/step-44/step-44.cc @@ -44,7 +44,7 @@ #include #include #include -#include +#include #include #include @@ -90,10 +90,10 @@ struct FESystem { void FESystem::declare_parameters(ParameterHandler &prm) { prm.enter_subsection("Finite element system"); { - prm.declare_entry("Polynomial degree", "2", Patterns::Integer(), + prm.declare_entry("Polynomial degree", "2", Patterns::Integer(0), "Displacement system polynomial order"); - prm.declare_entry("Quadrature order", "3", Patterns::Integer(), + prm.declare_entry("Quadrature order", "3", Patterns::Integer(0), "Gauss quadrature order"); } prm.leave_subsection(); @@ -127,10 +127,10 @@ struct Geometry { void Geometry::declare_parameters(ParameterHandler &prm) { prm.enter_subsection("Geometry"); { - prm.declare_entry("Global refinement", "2", Patterns::Integer(), + prm.declare_entry("Global refinement", "2", Patterns::Integer(0), "Global refinement level"); - prm.declare_entry("Grid scale", "1e-3", Patterns::Double(), + prm.declare_entry("Grid scale", "1e-3", Patterns::Double(0.0), "Global grid scaling factor"); prm.declare_entry("Pressure ratio p/p0", "100", @@ -186,13 +186,18 @@ void Materials::parse_parameters(ParameterHandler &prm) { } // @sect4{Linear solver} -// Choose both CG solver and SSOR preconditioner settings. -// The default values are optimal for this particular problem. +// Choose both solver and preconditioner settings. +// The use of an effective preconditioner is critical to ensure +// convergence when a large nonlinear motion occurs +// in a Newton increment. +// ToDo: explain +// The default values are optimal for single-thread conditions this particular problem. struct LinearSolver { std::string type_lin; double tol_lin; double max_iterations_lin; - double ssor_relaxation; + std::string preconditioner_type; + double preconditioner_relaxation; static void declare_parameters(ParameterHandler &prm); @@ -206,17 +211,20 @@ void LinearSolver::declare_parameters(ParameterHandler &prm) { prm.declare_entry("Solver type", "CG", Patterns::Selection("CG|Direct"), "Type of solver used to solve the linear system"); - prm.declare_entry("Residual", "1e-6", Patterns::Double(), + prm.declare_entry("Residual", "1e-6", Patterns::Double(0.0), "Linear solver residual (scaled by residual norm)"); prm.declare_entry( "Max iteration multiplier", "1", - Patterns::Double(), + Patterns::Double(0.0), "Linear solver iterations (multiples of the system matrix size)"); - prm.declare_entry("SSOR Relaxation", "0.65", Patterns::Double(), - "SSOR preconditioner relaxation value"); + prm.declare_entry("Preconditioner type", "ssor", Patterns::Selection("jacobi|ssor"), + "Type of preconditioner"); + + prm.declare_entry("Preconditioner relaxation", "0.65", Patterns::Double(0.0), + "Preconditioner relaxation value"); } prm.leave_subsection(); } @@ -227,7 +235,8 @@ void LinearSolver::parse_parameters(ParameterHandler &prm) { 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"); + preconditioner_type = prm.get("Preconditioner type"); + preconditioner_relaxation = prm.get_double("Preconditioner relaxation"); } prm.leave_subsection(); } @@ -252,14 +261,14 @@ void NonlinearSolver::declare_parameters(ParameterHandler &prm) { prm.enter_subsection("Nonlinear solver"); { prm.declare_entry("Max iterations Newton-Raphson", "10", - Patterns::Integer(), + Patterns::Integer(0), "Number of Newton-Raphson iterations allowed"); - prm.declare_entry("Tolerance force", "1.0e-9", Patterns::Double(), + prm.declare_entry("Tolerance force", "1.0e-9", Patterns::Double(0.0), "Force residual tolerance"); prm.declare_entry("Tolerance displacement", "1.0e-6", - Patterns::Double(), "Displacement error tolerance"); + Patterns::Double(0.0), "Displacement error tolerance"); } prm.leave_subsection(); } @@ -642,13 +651,13 @@ public: // Second derivative of the volumetric free energy wrt $\widetilde{J}$ // We need the following computation explicitly in the tangent // so we make it public. - // calculate + // We calculate // $\frac{\partial^2 \Psi_{\textrm{vol}}(\widetilde{J})}{\partial \widetilde{J} \partial \widetilde{J}}$ double get_d2Psi_vol_dJ2(void) const { return kappa * (1.0 + 1.0 / (J_tilde * J_tilde)); } - + // Return various data that we choose to store with the material double get_det_F(void) const { return det_F; } @@ -662,7 +671,7 @@ public: } protected: - // Model properties $\kappa$ and $c_1$ + // Define constitutive model paramaters $\kappa$ and $c_1$ const double kappa; // Bulk modulus const double c_1; // neo-Hookean model parameter @@ -670,7 +679,6 @@ protected: double det_F; double p_tilde; double J_tilde; - SymmetricTensor<2, dim> b_bar; // Determine the volumetric Kirchhoff stress @@ -685,7 +693,7 @@ protected: return AdditionalTools::StandardTensors::dev_P * get_tau_bar(); } - // Determine the fictitious Kirchhoff stress + // Determine the fictitious Kirchhoff stress $\overline{\boldsymbol{\tau}}$ SymmetricTensor<2, dim> get_tau_bar(void) const { return 2.0 * c_1 * b_bar; } @@ -715,7 +723,8 @@ protected: * AdditionalTools::StandardTensors::dev_P; } - // Calculate the fictitious elasticity tensor $\overline{\mathfrak{c}}$ + // Calculate the fictitious elasticity tensor $\overline{\mathfrak{c}}$. + // For the material model chosen this is simply zero. SymmetricTensor<4, dim> get_c_bar() const { SymmetricTensor<4, dim> c_bar; c_bar = 0.0; @@ -726,8 +735,12 @@ protected: // @sect3{Quadrature point history} // As seen in step-18, the PointHistory class offers // a method for storing data at the quadrature points. -// We need to evaluate the Kirchhoff stress $\boldsymbol{\tau}$ and -// the tangent $J\mathfrak{c}$ at the quadrature points. +// Here each quadrature point holds a pointer to a Material. +// Thus, different material models can be used in different regions +// of the domain. +// Among other data, we choose to store the +// Kirchhoff stress $\boldsymbol{\tau}$ and +// the tangent $J\mathfrak{c}$ for the quadrature points. template class PointHistory { @@ -748,7 +761,8 @@ public: // We first create a material object. void setup_lqp(Parameters::AllParameters & parameters) { - // Create an instance of a neo-Hookean material + // Create an instance of a three field + // compressible neo-Hookean material material = new Material_Compressible_Neo_Hook_Three_Field( parameters.mu, parameters.nu); @@ -757,17 +771,14 @@ public: } // Update the stored values and stresses based on the current - // deformation configuration, pressure $p$ and + // deformation measure $\textrm{Grad}\mathbf{u}_{\textrm{n}}$, + // pressure $\widetilde{p}$ and // dilation $\widetilde{J}$ field values. - // The input is the material gradient of the displacement - // $\textrm{Grad}\mathbf{u}_{\textrm{n}}$ void update_values(const Tensor<2, dim> & Grad_u_n, const double p_tilde, const double J_tilde) { - // Store the calculated pressure $p$ - // and dilatation $\widetilde{J}$ - // Various deformation gradient $\mathbf{F}$ from the + // Calculate the deformation gradient $\mathbf{F}$ from the // displacement gradient $\textrm{Grad}\mathbf{u}$, i.e. // $\mathbf{F}(\mathbf{u}) = \mathbf{I} + \textrm{Grad} \mathbf{u}$ static const Tensor<2, dim> I = @@ -778,7 +789,8 @@ public: // We use the inverse of $\mathbf{F}$ frequently so we store it F_inv = invert(F); - // Now we update the material model with the new deformation measures + // Now we update the material model with the new deformation measure, + // pressure and dilatation material->update_material_data(F, p_tilde, J_tilde); // The material has been updated so we now calculate the @@ -805,8 +817,7 @@ public: // and the kinetic variables. // These are used in the material and global - // tangent matrix and residual assembly operations - // so we compute these and store them. + // tangent matrix and residual assembly operations. double get_p_tilde(void) const { return material->get_p_tilde(); } @@ -831,15 +842,12 @@ private: // We specify that each QP has a copy of a material // type in case different materials are used // in different regions of the domain. - // This also - // deals with the issue of preventing data-races during - // multi-threading operations when using shared objects. Material_Compressible_Neo_Hook_Three_Field* material; - // These are all the volume, displacement and strain variables + // The inverse of the deformation gradient Tensor<2, dim> F_inv; - // and the stress-type variables + // the stress-type variables SymmetricTensor<2, dim> tau; double d2Psi_vol_dJ2; double dPsi_vol_dJ; @@ -849,6 +857,8 @@ private: }; // @sect3{Quasi-static quasi-incompressible finite-strain solid} +// The Solid class is the central class in that it represents +// the problem at hand. template class Solid { public: @@ -860,8 +870,10 @@ public: private: - // Threaded building-blocks data structures: - // for the tangent matrix + // Threaded building-blocks data structures + // for the tangent matrix. + // (see the module on @ref distributed + // for a definition, as well as the discussion in step-40) struct PerTaskData_K; struct ScratchData_K; // for the right-hand side @@ -909,11 +921,11 @@ private: ScratchData_SC & scratch, PerTaskData_SC & data); void copy_local_to_global_sc(const PerTaskData_SC & data); - // Apply Dirichlet boundary values + // Apply Dirichlet boundary conditions on the displacement field void make_constraints(const int & it_nr, ConstraintMatrix & constraints); - // Create and update the quadrature points stress and strain values + // Create and update the quadrature points void setup_qph(void); void @@ -925,7 +937,7 @@ private: void copy_local_to_global_UQPH(const PerTaskData_UQPH & data) { } - // Solve for the displacement using a Newton-Rhapson method + // Solve for the displacement using a Newton-Raphson method void solve_nonlinear_timestep(BlockVector & solution_delta); std::pair @@ -953,7 +965,8 @@ private: Time time; TimerOutput timer; - // A storage object for quadrature point information + // A storage object for quadrature point information. + // See step-18 for more on this std::vector > quadrature_point_history; // A description of the finite-element system including the displacement polynomial degree, @@ -969,7 +982,8 @@ private: // Description of how the block-system is arranged // There are 3 blocks, the first contains a vector DOF $\mathbf{u}$ - // while the other two describe scalar DOFs, $p$ and $\widetilde{J}$. + // while the other two describe scalar DOFs, + // $\widetilde{p}$ and $\widetilde{J}$. static const unsigned int n_blocks = 3; static const unsigned int n_components = dim + 2; static const unsigned int first_u_component = 0; @@ -995,6 +1009,7 @@ private: // Objects that store the converged solution and right-hand side vectors, // as well as the tangent matrix. There is a ConstraintMatrix object // used to keep track of constraints. + // We make use of a sparsity pattern designed for a block system. ConstraintMatrix constraints; BlockSparsityPattern sparsity_pattern; BlockSparseMatrix tangent_matrix; @@ -1033,7 +1048,7 @@ private: void get_error_update(const BlockVector & newton_update, Errors & error_update); - double + std::pair get_error_dil(void); // Print information to screen @@ -1092,8 +1107,8 @@ void Solid::run(void) { output_results(); time.increment(); - // Here we define - // $\varDelta \mathbf{\Xi}:= \{\varDelta \mathbf{u},\varDelta p, \varDelta \widetilde{J} \}$. + // Here we define the incremental solution update + // $\varDelta \mathbf{\Xi}:= \{\varDelta \mathbf{u},\varDelta \widetilde{p}, \varDelta \widetilde{J} \}$. BlockVector solution_delta(dofs_per_block); solution_delta.collect_sizes(); @@ -1108,8 +1123,9 @@ void Solid::run(void) { solve_nonlinear_timestep(solution_delta); // $\varDelta \mathbf{\Xi}_{\textrm{n}} = \varDelta \mathbf{\Xi}_{\textrm{n-1}} + \varDelta \mathbf{\Xi}$ solution_n += solution_delta; + // and plot the results output_results(); - + // we then move on happily to the next time step. time.increment(); } } @@ -1119,7 +1135,7 @@ void Solid::run(void) { // @sect4{Threaded-building-blocks structures} // We use TBB to perform as many computationally intensive // distributed tasks as possible. In particular, we assemble the -// tangent matrix and residual vector, the static +// tangent matrix and right hand side vector, the static // condensation contributions, and update data stored // at the quadrature points using TBB. @@ -1141,13 +1157,16 @@ struct Solid::PerTaskData_K { }; // while the ScratchData object stores the larger objects // such as the shape-function values object and a shape function -// gradient and symmetric gradient vector which we will precompute later. +// gradient and symmetric gradient vector which we will compute later. template struct Solid::ScratchData_K { FEValues fe_values_ref; + // interpolation function std::vector > Nx; + // their gradients std::vector > > grad_Nx; + // and their symmetric gradients. std::vector > > symm_grad_Nx; ScratchData_K(const FiniteElement & fe_cell, @@ -1187,7 +1206,7 @@ struct Solid::ScratchData_K { }; -// Next are the same data structures used for the +// Next are the same approach is used for the // right-hand side assembly. // The PerTaskData object again stores local contributions template @@ -1254,7 +1273,7 @@ struct Solid::ScratchData_RHS { // condensed tangent matrix. Recall that we wish to solve // for a displacement-based formulation. // We do the condensation at the element -// level as the $p$ and $\widetilde{J}$ +// level as the $\widetilde{p}$ and $\widetilde{J}$ // fields are element-wise discontinuous. // As these operations are matrix-based, // we need to setup a number of matrices @@ -1266,12 +1285,10 @@ struct Solid::PerTaskData_SC { FullMatrix cell_matrix; std::vector local_dof_indices; - // Calculation matrices (auto resized) FullMatrix k_orig; FullMatrix k_pu; FullMatrix k_pJ; FullMatrix k_JJ; - // Calculation matrices (manual resized) FullMatrix k_pJ_inv; FullMatrix k_bbar; FullMatrix A; @@ -1293,7 +1310,7 @@ struct Solid::PerTaskData_SC { C(n_p, n_u) { } - // Choose not to reset any data as the matrix extraction and + // We choose not to reset any data as the matrix extraction and // replacement tools will take care of this void reset(void) { } @@ -1544,7 +1561,6 @@ void Solid::setup_qph(void) { // only once the grid is refined to its finest level. { triangulation.clear_user_data(); - { std::vector > tmp; tmp.swap(quadrature_point_history); @@ -1601,7 +1617,7 @@ void Solid::update_qph_incremental( PerTaskData_UQPH per_task_data_UQPH; ScratchData_UQPH scratch_data_UQPH(fe, qf_cell, uf_UQPH, solution_total); - // and pass them and the one-cell update function to the workstream to be processed + // and pass them and the one-cell update function to the WorkStream to be processed 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, @@ -1634,7 +1650,6 @@ void Solid::update_qph_incremental_one_cell( // Firstly we need to 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, @@ -1653,10 +1668,10 @@ void Solid::update_qph_incremental_one_cell( } // @sect4{Solid::solve_nonlinear_timestep} +// The driver method for the Newton-Raphson scheme template void Solid::solve_nonlinear_timestep( BlockVector & solution_delta) { - // timer.enter_subsection("Nonlinear solver"); std::cout << std::endl << "Timestep " << time.get_timestep() << " @ " << time.current() << "s" << std::endl; @@ -1695,7 +1710,7 @@ void Solid::solve_nonlinear_timestep( // is an expensive operation and we can potentially avoid an extra // assembly process by not assembling the tangent matrix when convergence // is attained. - assemble_system_rhs(); // Assemble RHS + assemble_system_rhs(); get_error_residual(error_residual); // We store the residual errors after the first iteration @@ -1715,10 +1730,13 @@ void Solid::solve_nonlinear_timestep( return; } - assemble_system_tangent(); // Assemble stiffness matrix - make_constraints(it_nr, constraints); // Make boundary conditions - constraints.condense(tangent_matrix, system_rhs); // Apply BC's + // Now we assemble the tangent + assemble_system_tangent(); + // and make and impose the Dirichlet constraints + make_constraints(it_nr, constraints); + constraints.condense(tangent_matrix, system_rhs); + // Now we actually solve the linearised problem const std::pair lin_solver_output = solve_linear_system(newton_update); @@ -1779,11 +1797,13 @@ void Solid::print_conv_footer(void) { std::cout << "_"; std::cout << std::endl; + const std::pair error_dil = get_error_dil(); + std::cout << "Relative errors:" << std::endl << "Displacement:\t" << error_update.u / error_update_0.u << std::endl << "Force: \t\t" << error_residual.u / error_residual_0.u << std::endl - << "Dilatation:\t" << get_error_dil() << std::endl - << "v / V_0:\t" << vol_current << " / " << vol_reference << " = " << vol_current / vol_reference << std::endl; + << "Dilatation:\t" << error_dil.first << std::endl + << "v / V_0:\t" << vol_current << " / " << vol_reference << " = " << error_dil.second << std::endl; } @@ -1793,8 +1813,12 @@ void Solid::print_conv_footer(void) { // $ \bigl[ \int_{\Omega_0} {[ J - \widetilde{J}]}^{2}\textrm{d}V \bigr]^{1/2}$ // which is then normalised by the current volume // $\int_{\Omega_0} J ~\textrm{d}V = \int_\Omega ~\textrm{d}v$. +// We also return the ratio of the current volume of the domain +// to the reference volume. This is of interest for incompressible media +// where we want to check how well the isochoric constraint has been +// enforced. template -double Solid::get_error_dil(void) { +std::pair Solid::get_error_dil(void) { double dil_L2_error = 0.0; vol_current = 0.0; @@ -1824,7 +1848,10 @@ double Solid::get_error_dil(void) { }Assert(vol_current > 0, ExcInternalError()); } - return (std::sqrt(dil_L2_error)); + std::pair error_dil; + error_dil.first = std::sqrt(dil_L2_error); + error_dil.second = vol_current / vol_reference; + return error_dil; } // Determine the true residual error for the problem. @@ -2377,11 +2404,14 @@ std::pair Solid::solve_linear_system( GrowingVectorMemory > GVM; SolverCG > solver_CG(solver_control, GVM); - // We've chosen a SSOR preconditioner as it appears to provide + // We've chosen by default a SSOR preconditioner as it appears to provide // the fastest solver convergence characteristics for this problem. - PreconditionSSOR > preconditioner; - preconditioner.initialize(tangent_matrix.block(u_dof, u_dof), - parameters.ssor_relaxation); + // However, for multicore computing, the Jacobi preconditioner + // which is multithreaded may converge quicker for larger linear systems. + PreconditionSelector, Vector > preconditioner ( + parameters.preconditioner_type, + parameters.preconditioner_relaxation); + preconditioner.use_matrix(tangent_matrix.block(u_dof, u_dof)); solver_CG.solve(tangent_matrix.block(u_dof, u_dof), newton_update.block(u_dof), system_rhs.block(u_dof), @@ -2502,7 +2532,8 @@ void Solid::assemble_sc_one_cell( // We now extract the contribution of // the dof associated with the current cell // to the global stiffness matrix. - // The discontinuous nature of the p and J + // The discontinuous nature of the $\widetilde{p}$ + // and $\widetilde{J}$ // interpolations mean that their is no // coupling of the local contributions at the // global level. This is not the case with the u dof. @@ -2538,23 +2569,14 @@ void Solid::assemble_sc_one_cell( // added from the surrounding cells, so we need to be careful when we manipulate this block. // We can't just erase the subblocks. // - // 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 | 0 | K_pt^-1 - K_pt | - // | 0 | 0 | 0 | - // // This is the strategy we will employ 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 + // k_store: Since we don't have access to k_{uu}, but we know its contribution is added to the global + // K_{uu} matrix, we just want to add the element wise static-condensation k_bbar. + // k_{pJ}^-1: Similarly, k_pJ exists in the subblock. Since the copy operation is a += operation, we need + // to subtract the existing k_pJ submatrix in addition to "adding" that which we wish to // replace it with. - // K_{tp}^-1: Since the global matrix is symmetric, this block is the same as the one above - // and we can simply use K_pt^-1 as a substitute for this one + // k_{Jp}^-1: Since the global matrix is symmetric, this block is the same as the one above + // and we can simply use k_pJ^-1 as a substitute for this one // We first extract element data from the system matrix. So first // we get the entire subblock for the cell @@ -2624,7 +2646,7 @@ void Solid::output_results(void) const { // can be achieved without physically moving the grid points ourselves. // We first need to copy the solution to a temporary vector and then // create the Eulerian mapping. We also specify the polynomial degree - // to the DataOut object in order to produce a more refined output dataset + // to the DataOut object in order to produce a more refined output data set // when higher order polynomials are used. Vector soln(solution_n.size()); for (unsigned int i = 0; i < soln.size(); ++i)