From ba7721b611cb0ac8ab352e6cbd96f6a2eccf5308 Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 22 Apr 2011 21:41:27 +0000 Subject: [PATCH] Fix accidental checkin of a merge commit. git-svn-id: https://svn.dealii.org/trunk@23625 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-44/step-44.cc | 484 ++++++++++++++-------------- 1 file changed, 238 insertions(+), 246 deletions(-) diff --git a/deal.II/examples/step-44/step-44.cc b/deal.II/examples/step-44/step-44.cc index e8c203822e..55cf8cdbce 100644 --- a/deal.II/examples/step-44/step-44.cc +++ b/deal.II/examples/step-44/step-44.cc @@ -2,7 +2,7 @@ /* Authors: Jean-Paul Pelteret, University of Cape Town, */ /* Andrew McBride, University of Erlangen-Nuremberg, 2010 */ /* */ -/* Copyright (C) 2010 by the deal.II authors */ +/* Copyright (C) 2010, 2011 by the deal.II authors */ /* & Jean-Paul Pelteret and Andrew McBride */ /* */ /* This file is subject to QPL and may not be distributed */ @@ -63,7 +63,7 @@ struct FESystem { int poly_degree; int quad_order; - + static void declare_parameters (ParameterHandler &prm); void parse_parameters (ParameterHandler &prm); }; @@ -76,7 +76,7 @@ void FESystem::declare_parameters (ParameterHandler &prm) "1", Patterns::Integer(), "Displacement system polynomial order"); - + prm.declare_entry("Quadrature order", "2", Patterns::Integer(), @@ -101,7 +101,7 @@ struct Geometry int global_refinement; double scale; double p_p0; - + static void declare_parameters (ParameterHandler &prm); void parse_parameters (ParameterHandler &prm); }; @@ -114,12 +114,12 @@ void Geometry::declare_parameters (ParameterHandler &prm) "2", Patterns::Integer(), "Global refinement level"); - + prm.declare_entry("Grid scale", "1.0", Patterns::Double(), "Global grid scaling factor"); - + prm.declare_entry("Pressure ratio p/p0", "40", Patterns::Selection("20|40|60|80|100"), @@ -144,7 +144,7 @@ struct Materials { double nu; double mu; - + static void declare_parameters (ParameterHandler &prm); void parse_parameters (ParameterHandler &prm); }; @@ -157,7 +157,7 @@ void Materials::declare_parameters (ParameterHandler &prm) "0.49", Patterns::Double(), "Poisson's ratio"); - + prm.declare_entry("Shear modulus", "1.0e6", Patterns::Double(), @@ -183,7 +183,7 @@ struct LinearSolver double tol_lin; double max_iterations_lin; double ssor_relaxation; - + static void declare_parameters (ParameterHandler &prm); void parse_parameters (ParameterHandler &prm); }; @@ -196,17 +196,17 @@ void LinearSolver::declare_parameters (ParameterHandler &prm) "CG", Patterns::Selection("CG|Direct"), "Type of solver used to solve the linear system"); - + prm.declare_entry("Residual", "1e-6", Patterns::Double(), "Linear solver residual (scaled by residual norm)"); - + prm.declare_entry("Max iteration multiplier", "2", Patterns::Double(), "Linear solver iterations (multiples of the system matrix size)"); - + prm.declare_entry("SSOR Relaxation", "0.6", Patterns::Double(), @@ -233,7 +233,7 @@ struct NonlinearSolver unsigned int max_iterations_NR; double tol_f; double tol_u; - + static void declare_parameters (ParameterHandler &prm); void parse_parameters (ParameterHandler &prm); }; @@ -246,12 +246,12 @@ void NonlinearSolver::declare_parameters (ParameterHandler &prm) "10", Patterns::Integer(), "Number of Newton-Raphson iterations allowed"); - + prm.declare_entry("Tolerance force", "1.0e-9", Patterns::Double(), "Force residual tolerance"); - + prm.declare_entry("Tolerance displacement", "1.0e-3", Patterns::Double(), @@ -276,7 +276,7 @@ struct Time { double end_time; double delta_t; - + static void declare_parameters (ParameterHandler &prm); void parse_parameters (ParameterHandler &prm); }; @@ -289,7 +289,7 @@ void Time::declare_parameters (ParameterHandler &prm) "1", Patterns::Double(), "End time"); - + prm.declare_entry("Time step size", "0.1", Patterns::Double(), @@ -317,10 +317,10 @@ struct AllParameters public LinearSolver, public NonlinearSolver, public Time - + { AllParameters (const std::string & input_file); - + static void declare_parameters (ParameterHandler &prm); void parse_parameters (ParameterHandler &prm); }; @@ -363,20 +363,20 @@ void extract_submatrix(const std::vector< unsigned int > &row_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(); - + sub_matrix.reinit(n_rows_submatrix, n_cols_submatrix); - + 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())); - + sub_matrix(sub_row,sub_col) = matrix(row, col); } } @@ -392,17 +392,17 @@ void replace_submatrix(const std::vector< unsigned int > &row_index_set, 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) { 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())); - + matrix(row, col) = sub_matrix(sub_row, sub_col); - + } } } @@ -421,13 +421,13 @@ public: 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; @@ -449,7 +449,7 @@ public: 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) @@ -457,7 +457,7 @@ public: 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) { @@ -465,17 +465,17 @@ public: 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));} - + 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; @@ -497,17 +497,17 @@ public: 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); - + // Initialise all tensors correctly update_values (Tensor <2,dim> (), 0.0, 1.0); } - + // Total Variables void update_values (const Tensor<2, dim> & grad_u_n, const double & pressure, @@ -518,11 +518,11 @@ public: 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) @@ -530,68 +530,68 @@ public: 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) { return material->get_dU_dtheta(get_dilatation()); } - + 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;}; - + // 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); - + 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); } - + 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; @@ -617,36 +617,36 @@ 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; } }; - + 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) @@ -661,7 +661,7 @@ private: 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(), @@ -670,7 +670,7 @@ private: 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) { @@ -680,35 +680,35 @@ private: } } } - + }; - + 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, @@ -728,7 +728,7 @@ private: rhs_values (qf_cell.size(), Vector(dim)) { } - + ScratchData_F ( const ScratchData_F & rhs ) : fe_values_ref ( rhs.fe_values_ref.get_fe(), @@ -741,7 +741,7 @@ private: 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) { @@ -751,14 +751,14 @@ private: } } } - + }; - + struct PerTaskData_SC { FullMatrix cell_matrix; std::vector local_dof_indices; - + // Calculation matrices (auto resized) FullMatrix K_orig; FullMatrix K_pu; @@ -771,7 +771,7 @@ private: FullMatrix A; FullMatrix B; FullMatrix C; - + PerTaskData_SC (const unsigned int & dofs_per_cell, const unsigned int & n_u, const unsigned int & n_p, @@ -787,12 +787,12 @@ private: 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 { @@ -800,14 +800,14 @@ private: 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; @@ -815,7 +815,7 @@ private: 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, @@ -829,7 +829,7 @@ private: 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(), @@ -840,7 +840,7 @@ private: 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? @@ -852,15 +852,15 @@ private: } } }; - + // === 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, @@ -880,7 +880,7 @@ private: /// \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 @@ -892,61 +892,61 @@ private: /// \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 @@ -999,22 +999,17 @@ void Solid::run (void) system_setup (); output_results (); // Output initial grid position time.increment(); - + BlockVector solution_delta (dofs_per_block); solution_delta.collect_sizes (); -<<<<<<< .mine while (time.current() <= time.end()) { -======= - - while (time.current() <= time.end()) { ->>>>>>> .r23597 solution_delta = 0.0; - + // Solve step and update total solution vector solve_nonlinear_timestep (solution_delta); solution_n += solution_delta; - + output_results (); time.increment(); } @@ -1029,11 +1024,11 @@ void Solid::make_grid (void) Point (1.0, 1.0, 1.0), true ); GridTools::scale (parameters.scale, triangulation); - + // 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 typename Triangulation::active_cell_iterator cell = triangulation.begin_active(), @@ -1062,25 +1057,25 @@ template void Solid::system_setup (void) { timer.enter_subsection ("Setup system"); - + // 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 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 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 | @@ -1096,52 +1091,52 @@ void Solid::system_setup (void) else { coupling[ii][jj] = DoFTools::always; coupling[jj][ii] = DoFTools::always; - } + } } } - + // 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]; - + 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(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(); - + 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 residual.reinit (dofs_per_block); residual.collect_sizes (); - + solution_n.reinit (dofs_per_block); solution_n.collect_sizes (); solution_n.block(t_dof) = 1.0; // Dilatation is 1 in the initial configuration - + // Set up the quadrature point history setup_qph (); - + timer.leave_subsection(); } @@ -1152,7 +1147,7 @@ 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; @@ -1176,34 +1171,34 @@ 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()); } - + // Setup initial data typename DoFHandler::active_cell_iterator cell = dof_handler_ref.begin_active(), @@ -1212,7 +1207,7 @@ void Solid::setup_qph (void) 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 ); @@ -1226,7 +1221,7 @@ void Solid::update_qph_incremental (const BlockVector & solution_d { timer.enter_subsection("Update QPH data"); std::cout << "Update QPH data..."<< std::endl; - + // 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 ); @@ -1235,7 +1230,7 @@ void Solid::update_qph_incremental (const BlockVector & solution_d qf_cell, uf_UQPH, solution_total); - + WorkStream::run ( dof_handler_ref.begin_active(), dof_handler_ref.end(), *this, @@ -1243,7 +1238,7 @@ void Solid::update_qph_incremental (const BlockVector & solution_d &Solid::copy_local_to_global_UQPH, scratch_data_UQPH, per_task_data_UQPH); - + timer.leave_subsection(); } @@ -1255,17 +1250,17 @@ void Solid::update_qph_incremental_one_cell (const typename DoFHandler PointHistory* lqph = reinterpret_cast*> (cell->user_pointer()); Assert(lqph >= &quadrature_point_history.front(), ExcInternalError()); Assert(lqph < &quadrature_point_history.back(), ExcInternalError()); - + Assert(scratch.solution_grads_u_total.size() == n_q_points, ExcInternalError()); 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 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 for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) { @@ -1283,17 +1278,17 @@ void Solid::solve_nonlinear_timestep (BlockVector & solution_delta std::cout << "Timestep " << time.get_timestep() << std::endl; - + // Newton update vector BlockVector newton_update (dofs_per_block); newton_update.collect_sizes (); - + // 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 (); soln_error_update .collect_sizes (); - + 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) @@ -1302,10 +1297,10 @@ void Solid::solve_nonlinear_timestep (BlockVector & solution_delta << std::endl << "Newton iteration: " << it_nr << std::endl; - + tangent_matrix = 0.0; residual = 0.0; - + // Check residual make_constraints (it_nr, constraints); // Make boundary conditions assemble_system_F (); // Assemble RHS @@ -1313,7 +1308,7 @@ void Solid::solve_nonlinear_timestep (BlockVector & solution_delta // 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 if ( it_nr > 0 && res_u/res_u_0 <= parameters.tol_u @@ -1328,19 +1323,19 @@ void Solid::solve_nonlinear_timestep (BlockVector & solution_delta << "\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 - assemble_system_K (); - + 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 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 @@ -1349,11 +1344,11 @@ void Solid::solve_nonlinear_timestep (BlockVector & solution_delta 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 get_error_update(newton_update, soln_error_update); res_u = soln_error_update.block(u_dof).l2_norm(); - + // Residual scaling factors if (it_nr == 0) res_u_0 = res_u; std::cout @@ -1375,28 +1370,28 @@ void Solid::solve_nonlinear_timestep (BlockVector & solution_delta << "\t Relative force error: " << res_f/res_f_0 << "\t Dilatation error: " << get_error_dil() << std::endl; - + // 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) -{ +{ for (unsigned int i=0; i < dof_handler_ref.n_dofs(); ++i) - if (!constraints.is_constrained(i)) + if (!constraints.is_constrained(i)) error_res(i) = residual(i); } 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)) + if (!constraints.is_constrained(i)) error_update(i) = newton_update(i); } @@ -1405,9 +1400,9 @@ 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(); @@ -1416,13 +1411,13 @@ double Solid::get_error_dil (void) 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 } @@ -1432,7 +1427,7 @@ BlockVector Solid::get_solution_total (const BlockVector { BlockVector solution_total (solution_n); solution_total += solution_delta; - + return solution_total; } @@ -1441,22 +1436,19 @@ template void Solid::solve_linear_system (BlockVector & newton_update) { std::cout << "Solve linear system..." << std::endl; - + BlockVector A (dofs_per_block); BlockVector B (dofs_per_block); 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 | - + // Solve for du { -<<<<<<< .mine -======= ->>>>>>> .r23597 // 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)); @@ -1464,27 +1456,27 @@ void Solid::solve_linear_system (BlockVector & newton_update) 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") { 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); - + 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); - + 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() @@ -1500,7 +1492,7 @@ void Solid::solve_linear_system (BlockVector & newton_update) else throw (ExcMessage("Linear solver type not implemented")); timer.leave_subsection(); } - + timer.enter_subsection("Linear solver postprocessing"); // Postprocess for dp { @@ -1512,7 +1504,7 @@ void Solid::solve_linear_system (BlockVector & newton_update) 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 { // dt = Ktt^{-1} (Rt − Ktp dp) @@ -1529,14 +1521,14 @@ void Solid::assemble_system_K (void) { timer.enter_subsection("Assemble system matrix"); std::cout << "Assemble system matrix..."<< std::endl; - + tangent_matrix = 0.0; // Clear the matrix - + const UpdateFlags uf_cell ( update_values | update_gradients | update_JxW_values ); - + PerTaskData_K per_task_data (dofs_per_cell); // Initialise members of per_task_data to the correct sizes. ScratchData_K scratch_data (fe, qf_cell, uf_cell); - + WorkStream::run ( dof_handler_ref.begin_active(), dof_handler_ref.end(), *this, @@ -1544,7 +1536,7 @@ void Solid::assemble_system_K (void) &Solid::copy_local_to_global_K, scratch_data, per_task_data); - + timer.leave_subsection(); } @@ -1569,15 +1561,15 @@ void Solid::assemble_system_K_one_cell (const typename DoFHandler::act 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 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(); - + 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]); @@ -1593,7 +1585,7 @@ void Solid::assemble_system_K_one_cell (const typename DoFHandler::act } } } - + // Build cell stiffness matrix // Global and local system matrices are symmetric // => Take advantage of this: Build only the lower half of the local matrix @@ -1603,23 +1595,23 @@ void Solid::assemble_system_K_one_cell (const typename DoFHandler::act 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 @@ -1640,9 +1632,9 @@ void Solid::assemble_system_K_one_cell (const typename DoFHandler::act 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 for (unsigned int i=0; i::assemble_system_F (void) { timer.enter_subsection("Assemble system RHS"); std::cout << "Assemble system RHS..."<< std::endl; - + residual = 0.0; // Clear the vector - + const UpdateFlags uf_cell ( update_values | update_gradients | update_JxW_values ); const UpdateFlags uf_face ( update_values | update_normal_vectors | update_JxW_values); - + PerTaskData_F per_task_data (dofs_per_cell); // Initialise members of per_task_data to the correct sizes. ScratchData_F scratch_data (fe, qf_cell, uf_cell, qf_face, uf_face); - + WorkStream::run ( dof_handler_ref.begin_active(), dof_handler_ref.end(), *this, @@ -1678,7 +1670,7 @@ void Solid::assemble_system_F (void) &Solid::copy_local_to_global_F, scratch_data, per_task_data ); - + timer.leave_subsection(); } @@ -1701,14 +1693,14 @@ void Solid::assemble_system_F_one_cell (const typename DoFHandler::act 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 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(); - + for (unsigned int k=0; k::assemble_system_F_one_cell (const typename DoFHandler::act else Assert (k_group <= t_dof, ExcInternalError()); } } - + // Assembly for residual contribution for (unsigned int q_point=0; q_point < n_q_points; ++q_point) { const SymmetricTensor <2,dim> T = lqph[q_point].get_T_iso() + lqph[q_point].get_T_vol(); @@ -1729,14 +1721,14 @@ void Solid::assemble_system_F_one_cell (const typename DoFHandler::act 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); - + for (unsigned int i=0; i::assemble_system_F_one_cell (const typename DoFHandler::act else Assert (i_group <= t_dof, ExcInternalError()); } // END i LOOP } // END q_point LOOP - + // Assembly for Neumann RHS contribution if (cell->at_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 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; @@ -1789,7 +1781,7 @@ void Solid::assemble_system_F_one_cell (const typename DoFHandler::act } // END i LOOP } // END face q_point LOOP } // END at boundary check LOOP - + } // END face LOOP } } @@ -1799,13 +1791,13 @@ template void Solid::assemble_SC (void) { timer.enter_subsection("Perform static condensation"); - + PerTaskData_SC per_task_data (dofs_per_cell, element_indices_u.size(), element_indices_p.size(), element_indices_t.size()); // Initialise members of per_task_data to the correct sizes. ScratchData_SC scratch_data; - + WorkStream::run ( dof_handler_ref.begin_active(), dof_handler_ref.end(), *this, @@ -1813,7 +1805,7 @@ void Solid::assemble_SC (void) &Solid::copy_local_to_global_SC, scratch_data, per_task_data ); - + timer.leave_subsection(); } @@ -1836,7 +1828,7 @@ void Solid::assemble_SC_one_cell (const typename DoFHandler::active_ce 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 | @@ -1847,7 +1839,7 @@ void Solid::assemble_SC_one_cell (const typename DoFHandler::active_ce // | 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 // @@ -1874,9 +1866,9 @@ void Solid::assemble_SC_one_cell (const typename DoFHandler::active_ce // 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, tangent_matrix, @@ -1893,7 +1885,7 @@ void Solid::assemble_SC_one_cell (const typename DoFHandler::active_ce element_indices_t, data.K_orig, data.K_tt); - + // 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); @@ -1901,14 +1893,14 @@ void Solid::assemble_SC_one_cell (const typename DoFHandler::active_ce element_indices_t, data.K_pt_inv, data.cell_matrix); - + // 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 data.K_pt_inv.mmult(data.A, data.K_pu); data.K_tt.mmult(data.B, data.A); @@ -1926,25 +1918,25 @@ 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 - + // 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; - + 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 ); } @@ -1954,10 +1946,10 @@ void Solid::make_constraints (const int & it_nr, } { const int boundary_id = 2; - + 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 ); } @@ -1969,7 +1961,7 @@ void Solid::make_constraints (const int & it_nr, 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 ); } @@ -1981,7 +1973,7 @@ void Solid::make_constraints (const int & it_nr, 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 ); } @@ -1993,7 +1985,7 @@ void Solid::make_constraints (const int & it_nr, 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 ); } @@ -2001,7 +1993,7 @@ void Solid::make_constraints (const int & it_nr, VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); } } - + constraints.close(); } @@ -2010,15 +2002,15 @@ template void Solid::output_results(void) { DataOut data_out; - + 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); - + std::vector solution_name (dim, "displacement"); solution_name.push_back ("pressure"); solution_name.push_back ("dilatation"); - + data_out.attach_dof_handler (dof_handler_ref); data_out.add_data_vector (solution_n, solution_name, @@ -2030,12 +2022,12 @@ void Solid::output_results(void) for (unsigned int i=0; i < soln.size(); ++i) soln(i) = solution_n(i); MappingQEulerian q_mapping (degree, soln, dof_handler_ref); data_out.build_patches (q_mapping,degree); - + std::ostringstream filename; filename << "solution-" << time.get_timestep() << ".vtk"; - + std::ofstream output (filename.str().c_str()); data_out.write_vtk (output); } @@ -2046,7 +2038,7 @@ int main () try { deallog.depth_console (0); - + Solid<3> solid_3d ("parameters.prm"); solid_3d.run(); } @@ -2060,7 +2052,7 @@ int main () << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; - + return 1; } catch (...) @@ -2074,7 +2066,7 @@ int main () << std::endl; return 1; } - + return 0; } -- 2.39.5