diff --git a/packages/core.cmake b/packages/core.cmake index a75192465..98c369df4 100644 --- a/packages/core.cmake +++ b/packages/core.cmake @@ -1,480 +1,485 @@ #=============================================================================== # @file 00_core.cmake # # @author Guillaume Anciaux # @author Nicolas Richart # # @date creation: Mon Nov 21 2011 # @date last modification: Fri Sep 19 2014 # # @brief package description for core # # @section LICENSE # # Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # #=============================================================================== package_declare(core NOT_OPTIONAL DESCRIPTION "core package for Akantu") package_declare_sources(core common/aka_array.cc common/aka_array.hh common/aka_array_tmpl.hh common/aka_blas_lapack.hh common/aka_circular_array.hh common/aka_circular_array_inline_impl.cc common/aka_common.cc common/aka_common.hh common/aka_common_inline_impl.cc common/aka_csr.hh common/aka_element_classes_info_inline_impl.cc common/aka_error.cc common/aka_error.hh common/aka_event_handler_manager.hh common/aka_extern.cc common/aka_fwd.hh common/aka_grid_dynamic.hh common/aka_math.cc common/aka_math.hh common/aka_math_tmpl.hh common/aka_memory.cc common/aka_memory.hh common/aka_memory_inline_impl.cc common/aka_random_generator.hh common/aka_safe_enum.hh common/aka_static_memory.cc common/aka_static_memory.hh common/aka_static_memory_inline_impl.cc common/aka_static_memory_tmpl.hh common/aka_typelist.hh common/aka_types.hh common/aka_visitor.hh common/aka_voigthelper.hh common/aka_voigthelper.cc fe_engine/element_class.cc fe_engine/element_class.hh fe_engine/element_class_tmpl.hh fe_engine/element_classes/element_class_hexahedron_8_inline_impl.cc fe_engine/element_classes/element_class_hexahedron_20_inline_impl.cc fe_engine/element_classes/element_class_pentahedron_6_inline_impl.cc fe_engine/element_classes/element_class_pentahedron_15_inline_impl.cc fe_engine/element_classes/element_class_point_1_inline_impl.cc fe_engine/element_classes/element_class_quadrangle_4_inline_impl.cc fe_engine/element_classes/element_class_quadrangle_8_inline_impl.cc fe_engine/element_classes/element_class_segment_2_inline_impl.cc fe_engine/element_classes/element_class_segment_3_inline_impl.cc fe_engine/element_classes/element_class_tetrahedron_10_inline_impl.cc fe_engine/element_classes/element_class_tetrahedron_4_inline_impl.cc fe_engine/element_classes/element_class_triangle_3_inline_impl.cc fe_engine/element_classes/element_class_triangle_6_inline_impl.cc fe_engine/fe_engine.cc fe_engine/fe_engine.hh fe_engine/fe_engine_inline_impl.cc fe_engine/fe_engine_template.hh fe_engine/fe_engine_template_tmpl.hh fe_engine/geometrical_element.cc fe_engine/integration_element.cc fe_engine/integrator.hh fe_engine/integrator_gauss.hh fe_engine/integrator_gauss_inline_impl.cc fe_engine/interpolation_element.cc fe_engine/interpolation_element_tmpl.hh fe_engine/quadrature_point.hh fe_engine/shape_functions.hh fe_engine/shape_functions_inline_impl.cc fe_engine/shape_lagrange.cc fe_engine/shape_lagrange.hh fe_engine/shape_lagrange_inline_impl.cc fe_engine/shape_linked.cc fe_engine/shape_linked.hh fe_engine/shape_linked_inline_impl.cc fe_engine/element.hh fe_engine/quadrature_point.hh io/dumper/dumpable.hh io/dumper/dumpable.cc io/dumper/dumpable_dummy.hh io/dumper/dumpable_inline_impl.hh io/dumper/dumper_field.hh io/dumper/dumper_material_padders.hh io/dumper/dumper_filtered_connectivity.hh io/dumper/dumper_element_type.hh io/dumper/dumper_element_partition.hh io/mesh_io.cc io/mesh_io.hh io/mesh_io/mesh_io_abaqus.cc io/mesh_io/mesh_io_abaqus.hh io/mesh_io/mesh_io_diana.cc io/mesh_io/mesh_io_diana.hh io/mesh_io/mesh_io_msh.cc io/mesh_io/mesh_io_msh.hh io/model_io.cc io/model_io.hh io/parser/algebraic_parser.hh io/parser/input_file_parser.hh io/parser/parsable.cc io/parser/parsable.hh io/parser/parsable_tmpl.hh io/parser/parser.cc io/parser/parser.hh io/parser/parser_tmpl.hh io/parser/cppargparse/cppargparse.hh io/parser/cppargparse/cppargparse.cc io/parser/cppargparse/cppargparse_tmpl.hh mesh/element_group.cc mesh/element_group.hh mesh/element_group_inline_impl.cc mesh/element_type_map.hh mesh/element_type_map_tmpl.hh mesh/element_type_map_filter.hh mesh/group_manager.cc mesh/group_manager.hh mesh/group_manager_inline_impl.cc mesh/mesh.cc mesh/mesh.hh mesh/mesh_accessor.hh mesh/mesh_events.hh mesh/mesh_filter.hh mesh/mesh_data.cc mesh/mesh_data.hh mesh/mesh_data_tmpl.hh mesh/mesh_inline_impl.cc mesh/node_group.cc mesh/node_group.hh mesh/node_group_inline_impl.cc mesh_utils/mesh_partition.cc mesh_utils/mesh_partition.hh mesh_utils/mesh_partition/mesh_partition_mesh_data.cc mesh_utils/mesh_partition/mesh_partition_mesh_data.hh mesh_utils/mesh_partition/mesh_partition_scotch.hh mesh_utils/mesh_pbc.cc mesh_utils/mesh_utils.cc mesh_utils/mesh_utils.hh mesh_utils/mesh_utils_inline_impl.cc model/boundary_condition.hh model/boundary_condition_functor.hh model/boundary_condition_functor_inline_impl.cc model/boundary_condition_tmpl.hh model/dof_manager.cc model/dof_manager.hh model/dof_manager_inline_impl.cc model/model_solver.cc model/model_solver.hh model/non_linear_solver.hh + model/non_linear_solver.cc model/non_linear_solver_default.cc model/non_linear_solver_default.hh model/solver_callback.hh model/time_step_solver.hh + model/time_step_solvers/time_step_solver.cc model/time_step_solvers/time_step_solver_default.cc model/time_step_solvers/time_step_solver_default.hh model/time_step_solvers/time_step_solver_default_explicit.hh model/non_linear_solver_callback.hh - model/integration_scheme/integration_scheme.hh model/time_step_solvers/time_step_solver_default_solver_callback.hh model/integration_scheme/generalized_trapezoidal.hh - model/integration_scheme/generalized_trapezoidal_inline_impl.cc + model/integration_scheme/generalized_trapezoidal.cc + model/integration_scheme/integration_scheme.hh + model/integration_scheme/integration_scheme.cc model/integration_scheme/integration_scheme_1st_order.hh + model/integration_scheme/integration_scheme_1st_order.cc model/integration_scheme/integration_scheme_2nd_order.hh + model/integration_scheme/integration_scheme_2nd_order.cc model/integration_scheme/newmark-beta.hh - model/integration_scheme/newmark-beta_inline_impl.cc + model/integration_scheme/newmark-beta.cc model/model.cc model/model.hh model/model_inline_impl.cc model/solid_mechanics/material.cc model/solid_mechanics/material.hh model/solid_mechanics/material_inline_impl.cc model/solid_mechanics/material_selector.hh model/solid_mechanics/material_selector_tmpl.hh model/solid_mechanics/materials/internal_field.hh model/solid_mechanics/materials/internal_field_tmpl.hh model/solid_mechanics/materials/random_internal_field.hh model/solid_mechanics/materials/random_internal_field_tmpl.hh model/solid_mechanics/solid_mechanics_model.cc model/solid_mechanics/solid_mechanics_model.hh model/solid_mechanics/solid_mechanics_model_inline_impl.cc model/solid_mechanics/solid_mechanics_model_mass.cc model/solid_mechanics/solid_mechanics_model_material.cc model/solid_mechanics/solid_mechanics_model_tmpl.hh model/solid_mechanics/solid_mechanics_model_event_handler.hh model/solid_mechanics/materials/plane_stress_toolbox.hh model/solid_mechanics/materials/plane_stress_toolbox_tmpl.hh model/solid_mechanics/materials/material_core_includes.hh model/solid_mechanics/materials/material_elastic.cc model/solid_mechanics/materials/material_elastic.hh model/solid_mechanics/materials/material_elastic_inline_impl.cc model/solid_mechanics/materials/material_thermal.cc model/solid_mechanics/materials/material_thermal.hh model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc model/solid_mechanics/materials/material_elastic_linear_anisotropic.hh model/solid_mechanics/materials/material_elastic_orthotropic.cc model/solid_mechanics/materials/material_elastic_orthotropic.hh model/solid_mechanics/materials/material_damage/material_damage.hh model/solid_mechanics/materials/material_damage/material_damage_tmpl.hh model/solid_mechanics/materials/material_damage/material_marigo.cc model/solid_mechanics/materials/material_damage/material_marigo.hh model/solid_mechanics/materials/material_damage/material_marigo_inline_impl.cc model/solid_mechanics/materials/material_damage/material_mazars.cc model/solid_mechanics/materials/material_damage/material_mazars.hh model/solid_mechanics/materials/material_damage/material_mazars_inline_impl.cc model/solid_mechanics/materials/material_finite_deformation/material_neohookean.cc model/solid_mechanics/materials/material_finite_deformation/material_neohookean.hh model/solid_mechanics/materials/material_finite_deformation/material_neohookean_inline_impl.cc model/solid_mechanics/materials/material_plastic/material_plastic.cc model/solid_mechanics/materials/material_plastic/material_plastic.hh model/solid_mechanics/materials/material_plastic/material_plastic_inline_impl.cc model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening.cc model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening.hh model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening_inline_impl.cc model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.cc model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.hh solver/sparse_solver.cc solver/sparse_solver.hh solver/sparse_solver_inline_impl.cc solver/sparse_matrix.cc solver/sparse_matrix.hh solver/sparse_matrix_inline_impl.cc synchronizer/communication_buffer.hh synchronizer/communication_buffer_inline_impl.cc synchronizer/data_accessor.cc synchronizer/data_accessor.hh synchronizer/data_accessor_inline_impl.cc synchronizer/distributed_synchronizer.cc synchronizer/distributed_synchronizer.hh synchronizer/distributed_synchronizer_tmpl.hh synchronizer/dof_synchronizer.cc synchronizer/dof_synchronizer.hh synchronizer/dof_synchronizer_inline_impl.cc synchronizer/filtered_synchronizer.cc synchronizer/filtered_synchronizer.hh synchronizer/mpi_type_wrapper.hh synchronizer/pbc_synchronizer.cc synchronizer/pbc_synchronizer.hh synchronizer/real_static_communicator.hh synchronizer/static_communicator.cc synchronizer/static_communicator.hh synchronizer/static_communicator_dummy.hh synchronizer/static_communicator_inline_impl.hh synchronizer/synchronizer.cc synchronizer/synchronizer.hh synchronizer/synchronizer_registry.cc synchronizer/synchronizer_registry.hh ) package_declare_elements(core ELEMENT_TYPES _point_1 _segment_2 _segment_3 _triangle_3 _triangle_6 _quadrangle_4 _quadrangle_8 _tetrahedron_4 _tetrahedron_10 _pentahedron_6 _pentahedron_15 _hexahedron_8 _hexahedron_20 KIND regular GEOMETRICAL_TYPES _gt_point _gt_segment_2 _gt_segment_3 _gt_triangle_3 _gt_triangle_6 _gt_quadrangle_4 _gt_quadrangle_8 _gt_tetrahedron_4 _gt_tetrahedron_10 _gt_hexahedron_8 _gt_hexahedron_20 _gt_pentahedron_6 _gt_pentahedron_15 INTERPOLATION_TYPES _itp_lagrange_point_1 _itp_lagrange_segment_2 _itp_lagrange_segment_3 _itp_lagrange_triangle_3 _itp_lagrange_triangle_6 _itp_lagrange_quadrangle_4 _itp_serendip_quadrangle_8 _itp_lagrange_tetrahedron_4 _itp_lagrange_tetrahedron_10 _itp_lagrange_hexahedron_8 _itp_serendip_hexahedron_20 _itp_lagrange_pentahedron_6 _itp_lagrange_pentahedron_15 GEOMETRICAL_SHAPES _gst_point _gst_triangle _gst_square _gst_prism GAUSS_INTEGRATION_TYPES _git_point _git_segment _git_triangle _git_tetrahedron _git_pentahedron INTERPOLATION_KIND _itk_lagrangian FE_ENGINE_LISTS gradient_on_quadrature_points interpolate_on_quadrature_points interpolate compute_normals_on_control_points inverse_map contains compute_shapes compute_shapes_derivatives get_shapes_derivatives ) package_declare_material_infos(core LIST AKANTU_CORE_MATERIAL_LIST INCLUDE material_core_includes.hh ) package_declare_documentation_files(core manual.sty manual.cls manual.tex manual-macros.sty manual-titlepages.tex manual-introduction.tex manual-gettingstarted.tex manual-io.tex manual-solidmechanicsmodel.tex manual-constitutive-laws.tex manual-lumping.tex manual-elements.tex manual-appendix-elements.tex manual-appendix-materials.tex manual-appendix-packages.tex manual-backmatter.tex manual-bibliography.bib manual-bibliographystyle.bst figures/bc_and_ic_example.pdf figures/boundary.pdf figures/boundary.svg figures/dirichlet.pdf figures/dirichlet.svg figures/doc_wheel.pdf figures/doc_wheel.svg figures/dynamic_analysis.png figures/explicit_dynamic.pdf figures/explicit_dynamic.svg figures/static.pdf figures/static.svg figures/hooke_law.pdf figures/hot-point-1.png figures/hot-point-2.png figures/implicit_dynamic.pdf figures/implicit_dynamic.svg figures/insertion.pdf figures/interpolate.pdf figures/interpolate.svg figures/problemDomain.pdf_tex figures/problemDomain.pdf figures/static_analysis.png figures/stress_strain_el.pdf figures/tangent.pdf figures/tangent.svg figures/vectors.pdf figures/vectors.svg figures/stress_strain_neo.pdf figures/visco_elastic_law.pdf figures/isotropic_hardening_plasticity.pdf figures/stress_strain_visco.pdf figures/elements/hexahedron_8.pdf figures/elements/hexahedron_8.svg figures/elements/quadrangle_4.pdf figures/elements/quadrangle_4.svg figures/elements/quadrangle_8.pdf figures/elements/quadrangle_8.svg figures/elements/segment_2.pdf figures/elements/segment_2.svg figures/elements/segment_3.pdf figures/elements/segment_3.svg figures/elements/tetrahedron_10.pdf figures/elements/tetrahedron_10.svg figures/elements/tetrahedron_4.pdf figures/elements/tetrahedron_4.svg figures/elements/triangle_3.pdf figures/elements/triangle_3.svg figures/elements/triangle_6.pdf figures/elements/triangle_6.svg figures/elements/xtemp.pdf ) package_declare_documentation(core "This package is the core engine of \\akantu. It depends on:" "\\begin{itemize}" "\\item A C++ compiler (\\href{http://gcc.gnu.org/}{GCC} >= 4, or \\href{https://software.intel.com/en-us/intel-compilers}{Intel})." "\\item The cross-platform, open-source \\href{http://www.cmake.org/}{CMake} build system." "\\item The \\href{http://www.boost.org/}{Boost} C++ portable libraries." "\\item The \\href{http://www.zlib.net/}{zlib} compression library." "\\end{itemize}" "" "Under Ubuntu (14.04 LTS) the installation can be performed using the commands:" "\\begin{command}" " > sudo apt-get install cmake libboost-dev zlib1g-dev g++" "\\end{command}" "" "Under Mac OS X the installation requires the following steps:" "\\begin{itemize}" "\\item Install Xcode" "\\item Install the command line tools." "\\item Install the MacPorts project which allows to automatically" "download and install opensource packages." "\\end{itemize}" "Then the following commands should be typed in a terminal:" "\\begin{command}" " > sudo port install cmake gcc48 boost" "\\end{command}" ) find_program(READLINK_COMMAND readlink) find_program(ADDR2LINE_COMMAND addr2line) mark_as_advanced(READLINK_COMMAND) mark_as_advanced(ADDR2LINE_COMMAND) include(CheckFunctionExists) check_function_exists(clock_gettime _clock_gettime) if(NOT _clock_gettime) set(AKANTU_USE_OBSOLETE_GETTIMEOFDAY ON CACHE INTERNAL "" FORCE) else() set(AKANTU_USE_OBSOLETE_GETTIMEOFDAY OFF CACHE INTERNAL "" FORCE) endif() diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh index 632fdf86f..3451e72b7 100644 --- a/src/common/aka_common.hh +++ b/src/common/aka_common.hh @@ -1,464 +1,466 @@ /** * @file aka_common.hh * * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 * @date last modification: Mon Sep 15 2014 * * @brief common type descriptions for akantu * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @section DESCRIPTION * * All common things to be included in the projects files * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COMMON_HH__ #define __AKANTU_COMMON_HH__ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #define __BEGIN_AKANTU__ namespace akantu { #define __END_AKANTU__ }; /* -------------------------------------------------------------------------- */ #define __BEGIN_AKANTU_DUMPER__ namespace dumper { #define __END_AKANTU_DUMPER__ } /* -------------------------------------------------------------------------- */ #if defined(WIN32) # define __attribute__(x) #endif /* -------------------------------------------------------------------------- */ #include "aka_config.hh" #include "aka_error.hh" #include "aka_safe_enum.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* Common types */ /* -------------------------------------------------------------------------- */ typedef std::string ID; static const Real UINT_INIT_VALUE = Real(0.); #ifdef AKANTU_NDEBUG static const Real REAL_INIT_VALUE = Real(0.); #else static const Real REAL_INIT_VALUE = std::numeric_limits::quiet_NaN(); #endif /* -------------------------------------------------------------------------- */ /* Memory types */ /* -------------------------------------------------------------------------- */ typedef UInt MemoryID; typedef std::string Surface; typedef std::pair SurfacePair; typedef std::list< SurfacePair > SurfacePairList; /* -------------------------------------------------------------------------- */ extern const UInt _all_dimensions; /* -------------------------------------------------------------------------- */ /* Mesh/FEM/Model types */ /* -------------------------------------------------------------------------- */ __END_AKANTU__ #include "aka_element_classes_info.hh" __BEGIN_AKANTU__ /// small help to use names for directions enum SpacialDirection { _x = 0, _y = 1, _z = 2 }; /// enum MeshIOType type of mesh reader/writer enum MeshIOType { _miot_auto, ///< Auto guess of the reader to use based on the extension _miot_gmsh, ///< Gmsh files _miot_gmsh_struct, ///< Gsmh reader with reintpretation of elements has structures elements _miot_diana, ///< TNO Diana mesh format _miot_abaqus ///< Abaqus mesh format }; /// enum AnalysisMethod type of solving method used to solve the equation of motion enum AnalysisMethod { _static, _implicit_dynamic, _explicit_lumped_mass, _explicit_lumped_capacity, _explicit_consistent_mass }; /// enum DOFSupportType defines which kind of dof that can exists enum DOFSupportType { _dst_nodal, _dst_generic }; /// Type of non linear resolution available in akantu enum NonLinearSolverType { _nls_linear, ///< No non linear convergence loop _nls_newton_raphson, ///< Regular Newton-Raphson _nls_newton_raphson_modified, ///< Newton-Raphson with initial tangent }; /// Type of time stepping solver enum TimeStepSolverType { - _tsst_forward_euler, ///< GeneralizedTrapezoidal(0) - _tsst_trapezoidal_rule_1, ///< GeneralizedTrapezoidal(1/2) - _tsst_backward_euler, ///< GeneralizedTrapezoidal(1) - _tsst_central_difference, ///< NewmarkBeta(0, 1/2) - _tsst_fox_goodwin, ///< NewmarkBeta(1/6, 1/2) - _tsst_trapezoidal_rule_2, ///< NewmarkBeta(1/2, 1/2) - _tsst_linear_acceleration, ///< NewmarkBeta(1/3, 1/2) - _tsst_newmark_beta, ///< generic NewmarkBeta with user defined - /// alpha and beta + _tsst_forward_euler, ///< GeneralizedTrapezoidal(0) + _tsst_forward_euler_lumped, ///< GeneralizedTrapezoidal(0) with lumped capacity + _tsst_trapezoidal_rule_1, ///< GeneralizedTrapezoidal(1/2) + _tsst_backward_euler, ///< GeneralizedTrapezoidal(1) + _tsst_central_difference, ///< NewmarkBeta(0, 1/2) + _tsst_central_difference_lumped, ///< NewmarkBeta(0, 1/2) with lumped mass + _tsst_fox_goodwin, ///< NewmarkBeta(1/6, 1/2) + _tsst_trapezoidal_rule_2, ///< NewmarkBeta(1/2, 1/2) + _tsst_linear_acceleration, ///< NewmarkBeta(1/3, 1/2) + _tsst_newmark_beta, ///< generic NewmarkBeta with user defined + /// alpha and beta _tsst_generalized_trapezoidal, ///< generic GeneralizedTrapezoidal with user ///defined alpha }; /// enum SolveConvergenceMethod different resolution algorithms enum SolveConvergenceMethod { _scm_newton_raphson_tangent, ///< Newton-Raphson with tangent matrix _scm_newton_raphson_tangent_modified, ///< Newton-Raphson with constant tangent matrix _scm_newton_raphson_tangent_not_computed ///< Newton-Raphson with constant tangent matrix (user has to assemble it) }; /// enum SolveConvergenceCriteria different convergence criteria enum SolveConvergenceCriteria { _scc_residual, ///< Use residual to test the convergence _scc_solution, ///< Use solution to test the convergence _scc_residual_mass_wgh ///< Use residual weighted by inv. nodal mass to testb }; /// enum CohesiveMethod type of insertion of cohesive elements enum CohesiveMethod { _intrinsic, _extrinsic }; //! enum ContactResolutionMethod types of solving for the contact enum ContactResolutionMethod { _penalty, _lagrangian, _augmented_lagrangian, _nitsche, _mortar }; //! enum ContactImplementationMethod types for different contact implementations enum ContactImplementationMethod { _none, _uzawa, _generalized_newton }; /// myfunction(double * position, double * stress/force, double * normal, unsigned int material_id) typedef void (*BoundaryFunction)(double *,double *, double*, unsigned int); /// @enum BoundaryFunctionType type of function passed for boundary conditions enum BoundaryFunctionType { _bft_stress, _bft_traction, _bft_traction_local }; /// @enum SparseMatrixType type of sparse matrix used enum MatrixType { _unsymmetric, _symmetric }; /* -------------------------------------------------------------------------- */ /* Contact */ /* -------------------------------------------------------------------------- */ typedef ID ContactID; typedef ID ContactSearchID; typedef ID ContactNeighborStructureID; enum ContactType { _ct_not_defined = 0, _ct_2d_expli = 1, _ct_3d_expli = 2, _ct_rigid = 3 }; enum ContactSearchType { _cst_not_defined = 0, _cst_2d_expli = 1, _cst_expli = 2 }; enum ContactNeighborStructureType { _cnst_not_defined = 0, _cnst_regular_grid = 1, _cnst_2d_grid = 2 }; /* -------------------------------------------------------------------------- */ /* Friction */ /* -------------------------------------------------------------------------- */ typedef ID FrictionID; /* -------------------------------------------------------------------------- */ /* Ghosts handling */ /* -------------------------------------------------------------------------- */ typedef ID SynchronizerID; /// @enum CommunicatorType type of communication method to use enum CommunicatorType { _communicator_mpi, _communicator_dummy }; /// @enum SynchronizationTag type of synchronizations enum SynchronizationTag { //--- SolidMechanicsModel tags --- _gst_smm_mass, //< synchronization of the SolidMechanicsModel.mass _gst_smm_for_gradu, //< synchronization of the SolidMechanicsModel.displacement _gst_smm_boundary, //< synchronization of the boundary, forces, velocities and displacement _gst_smm_uv, //< synchronization of the nodal velocities and displacement _gst_smm_res, //< synchronization of the nodal residual _gst_smm_init_mat, //< synchronization of the data to initialize materials _gst_smm_stress, //< synchronization of the stresses to compute the internal forces _gst_smmc_facets, //< synchronization of facet data to setup facet synch _gst_smmc_facets_conn, //< synchronization of facet global connectivity _gst_smmc_facets_stress, //< synchronization of facets' stress to setup facet synch _gst_smmc_damage, //< synchronization of damage //--- CohesiveElementInserter tags --- _gst_ce_inserter, //< synchronization of global nodes id of newly inserted cohesive elements _gst_ce_groups, //< synchronization of cohesive element insertion depending on facet groups //--- GroupManager tags --- _gst_gm_clusters, //< synchronization of clusters //--- HeatTransfer tags --- _gst_htm_capacity, //< synchronization of the nodal heat capacity _gst_htm_temperature, //< synchronization of the nodal temperature _gst_htm_gradient_temperature, //< synchronization of the element gradient temperature //--- LevelSet tags --- /// synchronization of the nodal level set value phi _gst_htm_phi, /// synchronization of the element gradient phi _gst_htm_gradient_phi, //--- Material non local --- _gst_mnl_for_average, //< synchronization of data to average in non local material _gst_mnl_weight, //< synchronization of data for the weight computations //--- General tags --- _gst_test, //< Test tag _gst_user_1, //< tag for user simulations _gst_user_2, //< tag for user simulations _gst_material_id, //< synchronization of the material ids _gst_for_dump, //< everything that needs to be synch before dump //--- Contact & Friction --- _gst_cf_nodal, //< synchronization of disp, velo, and current position _gst_cf_incr, //< synchronization of increment ///--- Solver tags --- _gst_solver_solution //< synchronization of the solution obained with the PETSc solver }; /// standard output stream operator for SynchronizationTag inline std::ostream & operator <<(std::ostream & stream, SynchronizationTag type); /// @enum GhostType type of ghost enum GhostType { _not_ghost, _ghost, _casper // not used but a real cute ghost }; /* -------------------------------------------------------------------------- */ struct GhostType_def { typedef GhostType type; static const type _begin_ = _not_ghost; static const type _end_ = _casper; }; typedef safe_enum ghost_type_t; /// standard output stream operator for GhostType inline std::ostream & operator <<(std::ostream & stream, GhostType type); /// @enum SynchronizerOperation reduce operation that the synchronizer can perform enum SynchronizerOperation { _so_sum, _so_min, _so_max, _so_null }; /* -------------------------------------------------------------------------- */ /* Global defines */ /* -------------------------------------------------------------------------- */ #define AKANTU_MIN_ALLOCATION 2000 #define AKANTU_INDENT " " #define AKANTU_INCLUDE_INLINE_IMPL /* -------------------------------------------------------------------------- */ // int 2 type construct template struct Int2Type { enum { result = d }; }; // type 2 type construct template class Type2Type { typedef T OriginalType; }; /* -------------------------------------------------------------------------- */ template struct is_scalar { enum{ value = false }; }; #define AKANTU_SPECIFY_IS_SCALAR(type) \ template<> \ struct is_scalar { \ enum { value = true }; \ } AKANTU_SPECIFY_IS_SCALAR(Real); AKANTU_SPECIFY_IS_SCALAR(UInt); AKANTU_SPECIFY_IS_SCALAR(Int); AKANTU_SPECIFY_IS_SCALAR(bool); template < typename T1, typename T2 > struct is_same { enum { value = false }; // is_same represents a bool. }; template < typename T > struct is_same { enum { value = true }; }; /* -------------------------------------------------------------------------- */ #define AKANTU_SET_MACRO(name, variable, type) \ inline void set##name (type variable) { \ this->variable = variable; \ } #define AKANTU_GET_MACRO(name, variable, type) \ inline type get##name () const { \ return variable; \ } #define AKANTU_GET_MACRO_NOT_CONST(name, variable, type) \ inline type get##name () { \ return variable; \ } #define AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, \ support, con) \ inline con Array & \ get##name (const support & el_type, \ const GhostType & ghost_type = _not_ghost) con { \ return variable(el_type, ghost_type); \ } #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType,) #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, const) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType,) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, const) /* -------------------------------------------------------------------------- */ /// initialize the static part of akantu void initialize(int & argc, char ** & argv); /// initialize the static part of akantu and read the global input_file void initialize(const std::string & input_file, int & argc, char ** & argv); /* -------------------------------------------------------------------------- */ /// finilize correctly akantu and clean the memory void finalize (); /* -------------------------------------------------------------------------- */ /* * For intel compiler annoying remark */ #if defined(__INTEL_COMPILER) /// remark #981: operands are evaluated in unspecified order #pragma warning ( disable : 981 ) /// remark #383: value copied to temporary, reference to temporary used #pragma warning ( disable : 383 ) #endif //defined(__INTEL_COMPILER) /* -------------------------------------------------------------------------- */ /* string manipulation */ /* -------------------------------------------------------------------------- */ inline std::string to_lower(const std::string & str); /* -------------------------------------------------------------------------- */ inline std::string trim(const std::string & to_trim); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /// give a string representation of the a human readable size in bit template std::string printMemorySize(UInt size); /* -------------------------------------------------------------------------- */ __END_AKANTU__ #include "aka_fwd.hh" __BEGIN_AKANTU__ /// get access to the internal argument parser cppargparse::ArgumentParser & getStaticArgumentParser(); /// get access to the internal input file parser Parser & getStaticParser(); /// get access to the user part of the internal input file parser const ParserSection & getUserParser(); __END_AKANTU__ #include "aka_common_inline_impl.cc" #endif /* __AKANTU_COMMON_HH__ */ diff --git a/src/model/dof_manager.hh b/src/model/dof_manager.hh index 109d83abd..b314e20fb 100644 --- a/src/model/dof_manager.hh +++ b/src/model/dof_manager.hh @@ -1,304 +1,304 @@ /** * @file dof_manager.hh * * @author Nicolas Richart * * @date Wed Jul 22 11:43:43 2015 * * @brief Class handling the different types of dofs * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_memory.hh" #include "non_linear_solver.hh" #include "time_step_solver.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_DOF_MANAGER_HH__ #define __AKANTU_DOF_MANAGER_HH__ namespace akantu { class NonLinearSolverCallback; class SparseMatrixAIJ; } __BEGIN_AKANTU__ class DOFManager : protected Memory { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: DOFManager(const Mesh & mesh, const ID & id = "dof_manager", const MemoryID & memory_id = 0); virtual ~DOFManager(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// register an array of degree of freedom void registerDOFs(const ID & dof_id, Array & dofs_array, const DOFSupportType & support_type); /// register an array of derivatives for a particular dof array void registerDOFsDerivative(const ID & dof_id, UInt order, Array & dofs_derivative); /// register array representing the blocked degree of freedoms void registerBlockedDOFs(const ID & dof_id, Array & blocked_dofs); /// Assemble an array to the global residual array virtual void assembleToResidual(const ID & dof_id, const Array & array_to_assemble, Real scale_factor = 1.) = 0; /** * Assemble elementary values to a local array of the size nb_nodes * * nb_dof_per_node. The dof number is implicitly considered as * conn(el, n) * nb_nodes_per_element + d. * With 0 < n < nb_nodes_per_element and 0 < d < nb_dof_per_node **/ virtual void assembleElementalArrayLocalArray( const Array & elementary_vect, Array & array_assembeled, const ElementType & type, const GhostType & ghost_type, Real scale_factor = 1., const Array & filter_elements = empty_filter); /** * Assemble elementary values to the global residual array. The dof number is * implicitly considered as conn(el, n) * nb_nodes_per_element + d. * With 0 < n < nb_nodes_per_element and 0 < d < nb_dof_per_node **/ virtual void assembleElementalArrayResidual( const ID & dof_id, const Array & elementary_vect, const ElementType & type, const GhostType & ghost_type, Real scale_factor = 1., const Array & filter_elements = empty_filter); /** * Assemble elementary values to the global residual array. The dof number is * implicitly considered as conn(el, n) * nb_nodes_per_element + d. With 0 < * n < nb_nodes_per_element and 0 < d < nb_dof_per_node **/ virtual void assembleElementalMatricesToMatrix( const ID & matrix_id, const ID & dof_id, const Array & elementary_mat, const ElementType & type, const GhostType & ghost_type, const MatrixType & elemental_matrix_type = _symmetric, const Array & filter_elements = empty_filter) = 0; /// notation fully defined yet... // virtual void assemblePreassembledMatrix(const ID & matrix_id, // const ID & dof_id_m, // const ID & dof_id_n, // const Matrix & matrix) = 0; /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: /// splits the solution storage from a global view to the per dof storages void splitSolutionPerDOFs(); /// minimum fonctionality to implement per derived version of the DOFManager /// to allow the splitSolutionPerDOFs function to work virtual void getSolutionPerDOFs(const ID & dof_id, Array & solution_array) = 0; /// fill a Vector with the equation numbers corresponding to the given /// connectivity inline void extractElementEquationNumber( const Array & equation_numbers, const Vector & connectivity, UInt nb_degree_of_freedom, Vector & local_equation_number); /// register a matrix void registerSparseMatrix(const ID & matrix_id, SparseMatrix & matrix); /// register a time step solver instantiated by a derived class void registerTimeStepSolver(const ID & time_step_solver_id, TimeStepSolver & time_step_solver); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the equation numbers (in local numbering) corresponding to a dof ID Array & getLocalEquationNumbers(const ID & dof_id); /// get the equation numbers (in global numbering) corresponding to a dof ID Array & getGlobalEquationNumbers(const ID & dof_id); /// get the local equation number corresponding to the global equation number /// get the equation numbers (in local numbering) corresponding to a dof ID const Array & getLocalEquationNumbers(const ID & dof_id) const; /// get the equation numbers (in global numbering) corresponding to a dof ID const Array & getGlobalEquationNumbers(const ID & dof_id) const; /// Global number of dofs AKANTU_GET_MACRO(SystemSize, this->system_size, UInt); /// Local number of dofs AKANTU_GET_MACRO(LocalSystemSize, this->local_system_size, UInt); /// Get the solver callback stored in the dof_manager AKANTU_GET_MACRO(SolverCallback, *solver_callback, NonLinearSolverCallback &); /* ------------------------------------------------------------------------ */ /* DOFs and derivatives accessors */ /* ------------------------------------------------------------------------ */ /// Get a reference to the registered dof array for a given id Array & getDOFs(const ID & dofs_id); /// Get a reference to the registered dof derivatives array for a given id Array & getDOFsDerivatives(const ID & dofs_id, UInt order); // /// Get a reference to the blocked dofs array registered for the given id // Array & getBlockedDOFs(const ID & dofs_id); /// Get a reference to the blocked dofs array registered for the given id const Array & getBlockedDOFs(const ID & dofs_id) const; /// Get a reference to the solution array registered for the given id const Array & getSolution(const ID & dofs_id) const; // /// Get a reference to the right hand side array registered for the given // dof // /// id // const Array & getRHS() const; /* ------------------------------------------------------------------------ */ /* Matrices accessors */ /* ------------------------------------------------------------------------ */ /// Get an instance of a new SparseMatrix virtual SparseMatrix & getNewMatrix(const ID & matrix_id, const MatrixType & matrix_type); /// Get an instance of a new SparseMatrix as a copy of the SparseMatrix /// matrix_to_copy_id virtual SparseMatrix & getNewMatrix(const ID & matrix_id, const ID & matrix_to_copy_id); /// Get the reference of an existing matrix - SparseMatrixAIJ & getMatrix(const ID & matrix_id); + SparseMatrix & getMatrix(const ID & matrix_id); /* ------------------------------------------------------------------------ */ /* Non linear system solver */ /* ------------------------------------------------------------------------ */ /// Get instance of a non linear solver virtual NonLinearSolver & getNewNonLinearSolver(const ID & nls_solver_id, const NonLinearSolverType & _non_linear_solver_type); /// get instance of a non linear solver virtual NonLinearSolver & getNonLinearSolver(const ID & nls_solver_id); /* ------------------------------------------------------------------------ */ /* TimeStepSolver */ /* ------------------------------------------------------------------------ */ /// Get instance of a time step solver virtual TimeStepSolver & getNewTimeStepSolver(const ID & time_step_solver_id, const ID & dof_id, const TimeStepSolverType & type); /// get instance of a time step solver virtual TimeStepSolver & getTimeStepSolver(const ID & time_step_solver_id); /* ------------------------------------------------------------------------*/ /* Class Members */ /* ------------------------------------------------------------------------*/ protected: /// dof representations in the dof manager struct DOFData { /// DOF support type (nodal, general) this is needed to determine how the /// dof are shared among processors DOFSupportType support_type; /// Degree of freedom array Array * dof; /// Blocked degree of freedoms array Array * blocked_dofs; /// Solution associated to the dof Array * solution; /* ---------------------------------------------------------------------- */ /* data for dynamic simulations */ /* ---------------------------------------------------------------------- */ /// Degree of freedom derivatives arrays std::vector *> dof_derivatives; /// local numbering equation numbers Array local_equation_number; }; /// equation number in global numbering Array global_equation_number; /// type to store dofs informations typedef std::map DOFStorage; /// type to store all the matrices typedef std::map SparseMatricesMap; /// type to store all the non linear solver typedef std::map NonLinearSolversMap; /// type to store all the time step solver typedef std::map TimeStepSolversMap; /// store a reference to the dof arrays DOFStorage dofs; /// list of sparse matrices that where created SparseMatricesMap matrices; /// non linear solvers storage NonLinearSolversMap non_linear_solvers; /// time step solvers storage TimeStepSolversMap time_step_solvers; /// reference to the underlying mesh const Mesh & mesh; /// Solver callback to assemble residual and Jacobian matrix NonLinearSolverCallback * solver_callback; /// Total number of degrees of freedom (size with the ghosts) UInt local_system_size; /// Number of purely local dofs (size without the ghosts) UInt pure_local_system_size; /// Total number of degrees of freedom UInt system_size; }; __END_AKANTU__ #include "dof_manager_inline_impl.cc" #endif /* __AKANTU_DOF_MANAGER_HH__ */ diff --git a/src/model/dof_manager_default.cc b/src/model/dof_manager_default.cc index 2cde8b29d..4b7516a5b 100644 --- a/src/model/dof_manager_default.cc +++ b/src/model/dof_manager_default.cc @@ -1,326 +1,326 @@ /** * @file dof_manager_default.cc * * @author Nicolas Richart * * @date Tue Aug 11 16:21:01 2015 * * @brief Implementation of the default DOFManager * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "dof_manager_default.hh" #include "sparse_matrix_aij.hh" #include "time_step_solver_default.hh" #include "static_communicator.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ inline void DOFManagerDefault::addSymmetricElementalMatrixToSymmetric( SparseMatrixAIJ & matrix, const Matrix & elementary_mat, const Vector & equation_numbers, UInt max_size) { for (UInt i = 0; i < elementary_mat.rows(); ++i) { UInt c_irn = equation_numbers(i); if (c_irn < max_size) { for (UInt j = i; j < elementary_mat.cols(); ++j) { UInt c_jcn = equation_numbers(j); if (c_jcn < max_size) { matrix(c_irn, c_jcn) += elementary_mat(i, j); } } } } } /* -------------------------------------------------------------------------- */ inline void DOFManagerDefault::addUnsymmetricElementalMatrixToSymmetric( SparseMatrixAIJ & matrix, const Matrix & elementary_mat, const Vector & equation_numbers, UInt max_size) { for (UInt i = 0; i < elementary_mat.rows(); ++i) { UInt c_irn = equation_numbers(i); if (c_irn < max_size) { for (UInt j = 0; j < elementary_mat.cols(); ++j) { UInt c_jcn = equation_numbers(j); if (c_jcn < max_size) { if (c_jcn >= c_irn) { matrix(c_irn, c_jcn) += elementary_mat(i, j); } } } } } } /* -------------------------------------------------------------------------- */ inline void DOFManagerDefault::addElementalMatrixToUnsymmetric( SparseMatrixAIJ & matrix, const Matrix & elementary_mat, const Vector & equation_numbers, UInt max_size) { for (UInt i = 0; i < elementary_mat.rows(); ++i) { UInt c_irn = equation_numbers(i); if (c_irn < max_size) { for (UInt j = 0; j < elementary_mat.cols(); ++j) { UInt c_jcn = equation_numbers(j); if (c_jcn < max_size) { matrix(c_irn, c_jcn) += elementary_mat(i, j); } } } } } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ DOFManagerDefault::DOFManagerDefault(const Mesh & mesh, const ID & id, const MemoryID & memory_id) : DOFManager(mesh, id, memory_id) {} /* -------------------------------------------------------------------------- */ DOFManagerDefault::~DOFManagerDefault() {} /* -------------------------------------------------------------------------- */ void DOFManagerDefault::registerDOFs(const ID & dof_id, Array & dofs_array, DOFSupportType & support_type) { // stores the current numbers of dofs UInt local_nb_dofs = this->local_system_size; UInt pure_local_nb_dofs = this->pure_local_system_size; // update or create the dof_data DOFManager::registerDOFs(dof_id, dofs_array, support_type); // Count the number of pure local dofs per proc StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); UInt prank = comm.whoAmI(); UInt psize = comm.getNbProc(); Array nb_dofs_per_proc(psize); nb_dofs_per_proc(prank) = this->pure_local_system_size - pure_local_nb_dofs; comm.allGather(nb_dofs_per_proc); UInt first_global_dofs_id = std::accumulate( nb_dofs_per_proc.begin(), nb_dofs_per_proc.begin() + prank, 0); // nb local dofs to account for UInt nb_dofs = this->local_system_size - local_nb_dofs; DOFData & dof_data = *dofs[dof_id]; this->global_equation_number.resize(this->local_system_size); // set the equation numbers if (support_type == _dst_nodal) { UInt first_dof_id = local_nb_dofs; dof_data.local_equation_number.resize(nb_dofs); for (UInt d = 0; d < nb_dofs; ++d) { UInt local_eq_num = first_dof_id + d; dof_data.local_equation_number(d) = local_eq_num; this->global_equation_number(local_eq_num) = first_global_dofs_id + d; } } } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManagerDefault::getNewMatrix(const ID & matrix_id, const MatrixType & matrix_type) { std::stringstream sstr; sstr << this->id << ":" << matrix_id; SparseMatrix * sm = new SparseMatrixAIJ(*this, matrix_type, sstr.str(), this->memory_id); this->registerSparseMatrix(matrix_id, *sm); return *sm; } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManagerDefault::getNewMatrix(const ID & matrix_id, const ID & matrix_to_copy_id) { std::stringstream sstr; sstr << this->id << ":" << matrix_id; SparseMatrixAIJ & sm_to_copy = this->getMatrix(matrix_to_copy_id); SparseMatrix * sm = new SparseMatrixAIJ(sm_to_copy, sstr.str(), this->memory_id); this->registerSparseMatrix(matrix_id, *sm); return *sm; } /* -------------------------------------------------------------------------- */ SparseMatrixAIJ & DOFManagerDefault::getMatrix(const ID & matrix_id) { AIJMatrixMap::iterator it = this->aij_matrices.find(matrix_id); if (it == this->aij_matrices.end()) AKANTU_EXCEPTION("The matrix " << matrix_id << " does not exists in the DOFManager " << this->id); return *(it->second); } /* -------------------------------------------------------------------------- */ TimeStepSolver & DOFManagerDefault::getNewTimeStepSolver( - const ID & dof_id, const ID & time_step_solver_id, + const ID & time_step_solver_id, const TimeStepSolverType & time_step_solver_type) { std::stringstream sstr; sstr << this->id << ":" << time_step_solver_id; - TimeStepSolver * tss = new TimeStepSolverDefault( - *this, dof_id, time_step_solver_type, sstr.str(), this->memory_id); + TimeStepSolver * tss = new TimeStepSolverDefault(*this, time_step_solver_type, + sstr.str(), this->memory_id); this->registerTimeStepSolver(time_step_solver_id, *tss); return *tss; } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::getSolutionPerDOFs(const ID & dof_id, Array & solution_array) { AKANTU_DEBUG_IN(); const Array & equation_number = this->getLocalEquationNumbers(dof_id); UInt nb_degree_of_freedoms = solution_array.getSize() * solution_array.getNbComponent(); AKANTU_DEBUG_ASSERT( equation_number.getSize() == nb_degree_of_freedoms, "The array to get the solution does not have a correct size." << " (" << solution_array.getID() << ")"); Real * sol_it = solution_array.storage(); UInt * equ_it = equation_number.storage(); Array solution; for (UInt d = 0; d < nb_degree_of_freedoms; ++d, ++sol_it, ++equ_it) { (*sol_it) = solution(*equ_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleToResidual(const ID & dof_id, const Array & array_to_assemble, Real scale_factor) { AKANTU_DEBUG_IN(); const Array & equation_number = this->getLocalEquationNumbers(dof_id); UInt nb_degree_of_freedoms = array_to_assemble.getSize() * array_to_assemble.getNbComponent(); AKANTU_DEBUG_ASSERT(equation_number.getSize() == nb_degree_of_freedoms, "The array to assemble does not have a correct size." << " (" << array_to_assemble.getID() << ")"); Real * arr_it = array_to_assemble.storage(); UInt * equ_it = equation_number.storage(); for (UInt d = 0; d < nb_degree_of_freedoms; ++d, ++arr_it, ++equ_it) { residual(*equ_it) += scale_factor * (*arr_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleElementalMatricesToMatrix( const ID & matrix_id, const ID & dof_id, const Array & elementary_mat, const ElementType & type, const GhostType & ghost_type, const MatrixType & elemental_matrix_type, const Array & filter_elements) { AKANTU_DEBUG_IN(); const Array & equation_number = this->getLocalEquationNumbers(dof_id); SparseMatrixAIJ & A = this->getMatrix(matrix_id); UInt nb_element; if (ghost_type == _not_ghost) { nb_element = mesh.getNbElement(type); } else { AKANTU_DEBUG_TO_IMPLEMENT(); } UInt * filter_it = NULL; if (filter_elements != empty_filter) { nb_element = filter_elements.getSize(); filter_it = filter_elements.storage(); } else { nb_element = mesh.getNbElement(type, ghost_type); } AKANTU_DEBUG_ASSERT(elementary_mat.getSize() == nb_element, "The vector elementary_mat(" << elementary_mat.getID() << ") has not the good size."); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_degree_of_freedom = elementary_mat.getNbComponent() / (nb_nodes_per_element * nb_nodes_per_element); const Array connectivity = this->mesh.getConnectivity(type, ghost_type); Array::const_vector_iterator conn_begin = connectivity.begin(nb_nodes_per_element); Array::const_vector_iterator conn_it = conn_begin; UInt size_mat = nb_nodes_per_element * nb_degree_of_freedom; Vector local_eq_nb(nb_degree_of_freedom * nb_nodes_per_element); Array::const_matrix_iterator el_mat_it = elementary_mat.begin(size_mat, size_mat); for (UInt e = 0; e < nb_element; ++e, ++el_mat_it) { if (filter_it != NULL) conn_it = conn_begin + *filter_it; this->extractElementEquationNumber(equation_number, *conn_it, nb_degree_of_freedom, local_eq_nb); if (filter_it != NULL) ++filter_it; else ++conn_it; if (A.getMatrixType() == _symmetric) if (elemental_matrix_type == _symmetric) this->addSymmetricElementalMatrixToSymmetric(A, *el_mat_it, local_eq_nb, A.getSize()); else this->addUnsymmetricElementalMatrixToSymmetric( A, *el_mat_it, local_eq_nb, A.getSize()); else this->addElementalMatrixToUnsymmetric(A, *el_mat_it, local_eq_nb, A.getSize()); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/dof_manager_default.hh b/src/model/dof_manager_default.hh index b1b217e95..b3cf5ce28 100644 --- a/src/model/dof_manager_default.hh +++ b/src/model/dof_manager_default.hh @@ -1,174 +1,175 @@ /** * @file dof_manager_default.hh * * @author Nicolas Richart * * @date Tue Aug 11 14:06:18 2015 * * @brief Default implementation of the dof manager * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "dof_manager.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_DOF_MANAGER_DEFAULT_HH__ #define __AKANTU_DOF_MANAGER_DEFAULT_HH__ namespace akantu { class SparseMatrixAIJ; class NonLinearSolverDefault; class TimeStepSolverDefault; } __BEGIN_AKANTU__ class DOFManagerDefault : public DOFManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: DOFManagerDefault(const Mesh & mesh, const ID & id = "dof_manager_default", const MemoryID & memory_id = 0); virtual ~DOFManagerDefault(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// register an array of degree of freedom void registerDOFs(const ID & dof_id, Array & dofs_array, DOFSupportType & support_type); /// Assemble an array to the global residual array virtual void assembleToResidual(const ID & dof_id, const Array & array_to_assemble, Real scale_factor = 1.); /** * Assemble elementary values to the global residual array. The dof number is * implicitly considered as conn(el, n) * nb_nodes_per_element + d. * With 0 < n < nb_nodes_per_element and 0 < d < nb_dof_per_node **/ virtual void assembleElementalArrayResidual( const ID & dof_id, const Array & array_to_assemble, const ElementType & type, const GhostType & ghost_type, Real scale_factor = 1.); /** * Assemble elementary values to the global matrix. The dof number is * implicitly considered as conn(el, n) * nb_nodes_per_element + d. * With 0 < n < nb_nodes_per_element and 0 < d < nb_dof_per_node **/ virtual void assembleElementalMatricesToMatrix( const ID & matrix_id, const ID & dof_id, const Array & elementary_mat, const ElementType & type, const GhostType & ghost_type, const MatrixType & elemental_matrix_type, const Array & filter_elements); protected: /// Get the part of the solution corresponding to the dof_id virtual void getSolutionPerDOFs(const ID & dof_id, Array & solution_array); private: /// Add a symmetric matrices to a symmetric sparse matrix void addSymmetricElementalMatrixToSymmetric( SparseMatrixAIJ & matrix, const Matrix & element_mat, const Vector & equation_numbers, UInt max_size); /// Add a unsymmetric matrices to a symmetric sparse matrix (i.e. cohesive /// elements) void addUnsymmetricElementalMatrixToSymmetric( SparseMatrixAIJ & matrix, const Matrix & element_mat, const Vector & equation_numbers, UInt max_size); /// Add a matrices to a unsymmetric sparse matrix void addElementalMatrixToUnsymmetric(SparseMatrixAIJ & matrix, const Matrix & element_mat, const Vector & equation_numbers, UInt max_size); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// Get an instance of a new SparseMatrix virtual SparseMatrix & getNewMatrix(const ID & matrix_id, const MatrixType & matrix_type); /// Get an instance of a new SparseMatrix as a copy of the SparseMatrix /// matrix_to_copy_id virtual SparseMatrix & getNewMatrix(const ID & matrix_id, const ID & matrix_to_copy_id); /// Get the reference of an existing matrix SparseMatrixAIJ & getMatrix(const ID & matrix_id); /*------------------------------------------------------------------------- */ /// Get an instance of a time step solver virtual TimeStepSolver & - getNewTimeStepSolver(const ID & time_step_solver, const ID & dof_id, + getNewTimeStepSolver(const ID & time_step_solver, const TimeStepSolverType & time_step_solver_type); /* ------------------------------------------------------------------------ */ /// Get the solution array AKANTU_GET_MACRO_NOT_CONST(GlobalSolution, global_solution, Array &); /// Get the residual array AKANTU_GET_MACRO_NOT_CONST(Residual, residual, Array &); /// Get the blocked dofs array AKANTU_GET_MACRO(GlobalBlockedDOFs, global_blocked_dofs, const Array &); + /// Get the location type of a given dof bool isLocalOrMasterDOF(UInt dof_num); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: typedef std::map AIJMatrixMap; typedef std::map DefaultNonLinearSolversMap; typedef std::map DefaultTimeStepSolversMap; /// rhs to the system of equation corresponding to the residual linked to the /// different dofs Array residual; /// solution of the system of equation corresponding to the different dofs Array global_solution; /// blocked degree of freedom in the system equation corresponding to the /// different dofs Array global_blocked_dofs; /// define the dofs type, local, shared, ghost Array dofs_type; /// Map of the different matrices stored in the dof manager AIJMatrixMap aij_matrices; /// Map of the different time step solvers stored with there real type DefaultTimeStepSolversMap default_time_step_solver_map; }; __END_AKANTU__ #include "dof_manager_default_inline_impl.cc" #endif /* __AKANTU_DOF_MANAGER_DEFAULT_HH__ */ diff --git a/src/model/integration_scheme/generalized_trapezoidal.cc b/src/model/integration_scheme/generalized_trapezoidal.cc new file mode 100644 index 000000000..afd1db163 --- /dev/null +++ b/src/model/integration_scheme/generalized_trapezoidal.cc @@ -0,0 +1,172 @@ +/** + * @file generalized_trapezoidal.cc + * + * @author Nicolas Richart + * + * @date creation: Mon Jul 04 2011 + * @date last modification: Thu Jun 05 2014 + * + * @brief implementation of inline functions + * + * @section LICENSE + * + * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) any + * later version. + * + * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR + * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with Akantu. If not, see . + * + */ + +/* -------------------------------------------------------------------------- */ +#include "generalized_trapezoidal.hh" +#include "aka_array.hh" +#include "dof_manager.hh" +#include "sparse_matrix.hh" +/* -------------------------------------------------------------------------- */ + +__BEGIN_AKANTU__ + +/* -------------------------------------------------------------------------- */ +void GeneralizedTrapezoidal::predictor(Real delta_t, Array & u, + Array & u_dot, + const Array & blocked_dofs) const { + AKANTU_DEBUG_IN(); + + UInt nb_nodes = u.getSize(); + UInt nb_degree_of_freedom = u.getNbComponent() * nb_nodes; + + Real * u_val = u.storage(); + Real * u_dot_val = u_dot.storage(); + bool * blocked_dofs_val = blocked_dofs.storage(); + + for (UInt d = 0; d < nb_degree_of_freedom; d++) { + if (!(*blocked_dofs_val)) { + *u_val += (1. - alpha) * delta_t * *u_dot_val; + } + u_val++; + u_dot_val++; + blocked_dofs_val++; + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void GeneralizedTrapezoidal::corrector(const SolutionType & type, Real delta_t, + Array & u, Array & u_dot, + const Array & blocked_dofs, + const Array & delta) const { + AKANTU_DEBUG_IN(); + + switch (type) { + case _temperature: + this->allCorrector<_temperature>(delta_t, u, u_dot, blocked_dofs, delta); + break; + case _temperature_rate: + this->allCorrector<_temperature_rate>(delta_t, u, u_dot, blocked_dofs, + delta); + break; + default: + AKANTU_EXCEPTION("The corrector type : " + << type + << " is not supported by this type of integration scheme"); + } + + AKANTU_DEBUG_OUT(); +} +/* -------------------------------------------------------------------------- */ +Real +GeneralizedTrapezoidal::getTemperatureCoefficient(const SolutionType & type, + Real delta_t) const { + switch (type) { + case _temperature: + return 1.; + case _temperature_rate: + return alpha * delta_t; + default: + AKANTU_EXCEPTION("The corrector type : " + << type + << " is not supported by this type of integration scheme"); + } +} + +/* -------------------------------------------------------------------------- */ +Real +GeneralizedTrapezoidal::getTemperatureRateCoefficient(const SolutionType & type, + Real delta_t) const { + switch (type) { + case _temperature: + return 1. / (alpha * delta_t); + case _temperature_rate: + return 1.; + default: + AKANTU_EXCEPTION("The corrector type : " + << type + << " is not supported by this type of integration scheme"); + } +} + +/* -------------------------------------------------------------------------- */ +template +void GeneralizedTrapezoidal::allCorrector(Real delta_t, Array & u, + Array & u_dot, + const Array & blocked_dofs, + const Array & delta) const { + AKANTU_DEBUG_IN(); + + UInt nb_nodes = u.getSize(); + UInt nb_degree_of_freedom = u.getNbComponent() * nb_nodes; + + Real e = getTemperatureCoefficient(type, delta_t); + Real d = getTemperatureRateCoefficient(type, delta_t); + + Real * u_val = u.storage(); + Real * u_dot_val = u_dot.storage(); + Real * delta_val = delta.storage(); + bool * blocked_dofs_val = blocked_dofs.storage(); + + for (UInt dof = 0; dof < nb_degree_of_freedom; dof++) { + if (!(*blocked_dofs_val)) { + *u_val += e * *delta_val; + *u_dot_val += d * *delta_val; + } + u_val++; + u_dot_val++; + delta_val++; + blocked_dofs_val++; + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void GeneralizedTrapezoidal::assembleJacobian(const SolutionType & type, + Real delta_t) { + AKANTU_DEBUG_IN(); + + SparseMatrix & J = this->dof_manager.getMatrix("J"); + + const SparseMatrix & M = this->dof_manager.getMatrix("M"); + const SparseMatrix & K = this->dof_manager.getMatrix("K"); + + Real c = this->getTemperatureRateCoefficient(type, delta_t); + Real e = this->getTemperatureCoefficient(type, delta_t); + + J.add(M, e); + J.add(K, c); + + AKANTU_DEBUG_OUT(); +} + +__END_AKANTU__ diff --git a/src/model/integration_scheme/generalized_trapezoidal.hh b/src/model/integration_scheme/generalized_trapezoidal.hh index 7c1fe0bce..0f647ae29 100644 --- a/src/model/integration_scheme/generalized_trapezoidal.hh +++ b/src/model/integration_scheme/generalized_trapezoidal.hh @@ -1,161 +1,155 @@ /** * @file generalized_trapezoidal.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Mon Jul 04 2011 * @date last modification: Thu Jun 05 2014 * * @brief Generalized Trapezoidal Method. This implementation is taken from * Méthodes numériques en mécanique des solides by Alain Curnier \note{ISBN: * 2-88074-247-1} * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_GENERALIZED_TRAPEZOIDAL_HH__ #define __AKANTU_GENERALIZED_TRAPEZOIDAL_HH__ #include "integration_scheme_1st_order.hh" __BEGIN_AKANTU__ /** * The two differentiate equation (thermal and kinematic) are : * @f{eqnarray*}{ * C\dot{u}_{n+1} + Ku_{n+1} = q_{n+1}\\ * u_{n+1} = u_{n} + (1-\alpha) \Delta t \dot{u}_{n} + \alpha \Delta t *\dot{u}_{n+1} * @f} * * To solve it : * Predictor : * @f{eqnarray*}{ * u^0_{n+1} &=& u_{n} + (1-\alpha) \Delta t v_{n} \\ * \dot{u}^0_{n+1} &=& \dot{u}_{n} * @f} * * Solve : * @f[ (a C + b K^i_{n+1}) w = q_{n+1} - f^i_{n+1} - C \dot{u}^i_{n+1} @f] * * Corrector : * @f{eqnarray*}{ * \dot{u}^{i+1}_{n+1} &=& \dot{u}^{i}_{n+1} + a w \\ * u^{i+1}_{n+1} &=& u^{i}_{n+1} + b w * @f} * * a and b depends on the resolution method : temperature (u) or temperature *rate (@f$\dot{u}@f$) * * For temperature : @f$ w = \delta u, a = 1 / (\alpha \Delta t) , b = 1 @f$ @n * For temperature rate : @f$ w = \delta \dot{u}, a = 1, b = \alpha \Delta t @f$ */ class GeneralizedTrapezoidal : public IntegrationScheme1stOrder { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - GeneralizedTrapezoidal(Real alpha) : alpha(alpha){}; + GeneralizedTrapezoidal(DOFManager & dof_manager, Real alpha) : IntegrationScheme1stOrder(dof_manager), alpha(alpha){}; virtual ~GeneralizedTrapezoidal(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: - virtual void integrationSchemePred(Real delta_t, Array & u, - Array & u_dot, - const Array & blocked_dofs) const; + virtual void predictor(Real delta_t, Array & u, Array & u_dot, + const Array & blocked_dofs) const; - virtual void integrationSchemeCorrTemp(Real delta_t, Array & u, - Array & u_dot, - const Array & blocked_dofs, - const Array & delta) const; + virtual void corrector(const SolutionType & type, Real delta_t, + Array & u, Array & u_dot, + const Array & blocked_dofs, + const Array & delta) const; - virtual void integrationSchemeCorrTempRate(Real delta_t, Array & u, - Array & u_dot, - const Array & blocked_dofs, - const Array & delta) const; + virtual void assembleJacobian(const SolutionType & type, + Real time_step); public: /// the coeffichent @f{b@f} in the description - template - Real getTemperatureCoefficient(Real delta_t) const; + Real getTemperatureCoefficient(const SolutionType & type, Real delta_t) const; /// the coeffichent @f{a@f} in the description - template - Real getTemperatureRateCoefficient(Real delta_t) const; + Real getTemperatureRateCoefficient(const SolutionType & type, + Real delta_t) const; private: - template - void integrationSchemeCorr(Real delta_t, Array & u, Array & u_dot, - const Array & blocked_dofs, const Array & delta) const; + template + void allCorrector(Real delta_t, Array & u, Array & u_dot, + const Array & blocked_dofs, + const Array & delta) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Alpha, alpha, Real); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// the @f$\alpha@f$ parameter const Real alpha; }; /* -------------------------------------------------------------------------- */ -/* inline functions */ /* -------------------------------------------------------------------------- */ -#if defined(AKANTU_INCLUDE_INLINE_IMPL) -#include "generalized_trapezoidal_inline_impl.cc" -#endif - /** * Forward Euler (explicit) -> condition on delta_t */ class ForwardEuler : public GeneralizedTrapezoidal { public: - ForwardEuler() : GeneralizedTrapezoidal(0.){}; + ForwardEuler(DOFManager & dof_manager) : GeneralizedTrapezoidal(dof_manager, 0.){}; }; /** * Trapezoidal rule (implicit), midpoint rule or Crank-Nicolson */ class TrapezoidalRule1 : public GeneralizedTrapezoidal { public: - TrapezoidalRule1() : GeneralizedTrapezoidal(.5){}; + TrapezoidalRule1(DOFManager & dof_manager) : GeneralizedTrapezoidal(dof_manager, .5){}; }; /** * Backward Euler (implicit) */ class BackwardEuler : public GeneralizedTrapezoidal { public: - BackwardEuler() : GeneralizedTrapezoidal(1.){}; + BackwardEuler(DOFManager & dof_manager) : GeneralizedTrapezoidal(dof_manager, 1.){}; }; +/* -------------------------------------------------------------------------- */ + __END_AKANTU__ #endif /* __AKANTU_GENERALIZED_TRAPEZOIDAL_HH__ */ diff --git a/src/model/integration_scheme/generalized_trapezoidal_inline_impl.cc b/src/model/integration_scheme/generalized_trapezoidal_inline_impl.cc deleted file mode 100644 index ba822e26a..000000000 --- a/src/model/integration_scheme/generalized_trapezoidal_inline_impl.cc +++ /dev/null @@ -1,146 +0,0 @@ -/** - * @file generalized_trapezoidal_inline_impl.cc - * - * @author Nicolas Richart - * - * @date creation: Mon Jul 04 2011 - * @date last modification: Thu Jun 05 2014 - * - * @brief implementation of inline functions - * - * @section LICENSE - * - * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) - * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) - * - * Akantu is free software: you can redistribute it and/or modify it under the - * terms of the GNU Lesser General Public License as published by the Free - * Software Foundation, either version 3 of the License, or (at your option) any - * later version. - * - * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY - * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR - * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more - * details. - * - * You should have received a copy of the GNU Lesser General Public License - * along with Akantu. If not, see . - * - */ - -/* -------------------------------------------------------------------------- */ - -inline void GeneralizedTrapezoidal::integrationSchemePred(Real delta_t, - Array & u, - Array & u_dot, - const Array & blocked_dofs) const { - AKANTU_DEBUG_IN(); - - UInt nb_nodes = u.getSize(); - UInt nb_degree_of_freedom = u.getNbComponent() * nb_nodes; - - Real * u_val = u.storage(); - Real * u_dot_val = u_dot.storage(); - bool * blocked_dofs_val = blocked_dofs.storage(); - - for (UInt d = 0; d < nb_degree_of_freedom; d++) { - if(!(*blocked_dofs_val)) { - *u_val += (1. - alpha) * delta_t * *u_dot_val; - } - u_val++; - u_dot_val++; - blocked_dofs_val++; - } - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -inline void GeneralizedTrapezoidal::integrationSchemeCorrTemp(Real delta_t, - Array & u, - Array & u_dot, - const Array & blocked_dofs, - const Array & delta) const { - AKANTU_DEBUG_IN(); - - integrationSchemeCorr(delta_t, - u, - u_dot, - blocked_dofs, - delta); - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -inline void GeneralizedTrapezoidal::integrationSchemeCorrTempRate(Real delta_t, - Array & u, - Array & u_dot, - const Array & blocked_dofs, - const Array & delta) const { - AKANTU_DEBUG_IN(); - - integrationSchemeCorr(delta_t, - u, - u_dot, - blocked_dofs, - delta); - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -template<> -inline Real GeneralizedTrapezoidal::getTemperatureCoefficient(__attribute__ ((unused)) Real delta_t) const { - return 1.; -} - -template<> -inline Real GeneralizedTrapezoidal::getTemperatureRateCoefficient(Real delta_t) const { - return 1./(alpha * delta_t); -} - -/* -------------------------------------------------------------------------- */ -template<> -inline Real GeneralizedTrapezoidal::getTemperatureCoefficient(Real delta_t) const { - return alpha * delta_t; -} - -template<> -inline Real GeneralizedTrapezoidal::getTemperatureRateCoefficient(__attribute__ ((unused)) Real delta_t) const { - return 1.; -} - -/* -------------------------------------------------------------------------- */ -template -inline void GeneralizedTrapezoidal::integrationSchemeCorr(Real delta_t, - Array & u, - Array & u_dot, - const Array & blocked_dofs, - const Array & delta) const { - AKANTU_DEBUG_IN(); - - UInt nb_nodes = u.getSize(); - UInt nb_degree_of_freedom = u.getNbComponent() * nb_nodes; - - Real e = getTemperatureCoefficient(delta_t); - Real d = getTemperatureRateCoefficient(delta_t); - - Real * u_val = u.storage(); - Real * u_dot_val = u_dot.storage(); - Real * delta_val = delta.storage(); - bool * blocked_dofs_val = blocked_dofs.storage(); - - for (UInt dof = 0; dof < nb_degree_of_freedom; dof++) { - if(!(*blocked_dofs_val)) { - *u_val += e * *delta_val; - *u_dot_val += d * *delta_val; - } - u_val++; - u_dot_val++; - delta_val++; - blocked_dofs_val++; - } - - AKANTU_DEBUG_OUT(); -} diff --git a/src/model/time_step_solvers/time_step_solver.cc b/src/model/integration_scheme/integration_scheme.cc similarity index 71% copy from src/model/time_step_solvers/time_step_solver.cc copy to src/model/integration_scheme/integration_scheme.cc index 8616755a7..b665cfc1a 100644 --- a/src/model/time_step_solvers/time_step_solver.cc +++ b/src/model/integration_scheme/integration_scheme.cc @@ -1,42 +1,40 @@ /** - * @file time_step_solver.cc + * @file integration_scheme.cc * * @author Nicolas Richart * - * @date Mon Oct 12 16:56:43 2015 + * @date Fri Oct 23 15:33:36 2015 * - * @brief Implementation of common part of TimeStepSolvers + * @brief Common interface to all interface schemes * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ -#include "time_step_solver.hh" +#include "integration_scheme.hh" +#include "dof_manager.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ -TimeStepSolver::TimeStepSolver(DOFManager & dof_manager, const ID & dof_id, - const TimeStepSolverType & type, const ID & id, - UInt memory_id) - : Memory(id, memory_id), _dof_manager(dof_manager), dof_id(dof_id), - type(type) {} +IntegrationScheme::IntegrationScheme(DOFManager & dof_manager, UInt order) + : dof_manager(dof_manager), order(order) {} __END_AKANTU__ diff --git a/src/model/integration_scheme/integration_scheme.hh b/src/model/integration_scheme/integration_scheme.hh index 1f385f269..5f4a3fc0e 100644 --- a/src/model/integration_scheme/integration_scheme.hh +++ b/src/model/integration_scheme/integration_scheme.hh @@ -1,61 +1,91 @@ /** * @file integration_scheme.hh * * @author Nicolas Richart * * @date Mon Sep 28 10:43:18 2015 * * @brief This class is just a base class for the integration schemes * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ +#include "aka_common.hh" +/* -------------------------------------------------------------------------- */ #ifndef __AKANTU_INTEGRATION_SCHEME_HH__ #define __AKANTU_INTEGRATION_SCHEME_HH__ +namespace akantu { +class DOFManager; +} + __BEGIN_AKANTU__ class IntegrationScheme { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - IntegrationScheme(UInt order) : order(order) {} + enum SolutionType { + _displacement = 0, + _temperature = 0, + _velocity = 1, + _temperature_rate = 1, + _acceleration = 2, + }; + + IntegrationScheme(DOFManager & dof_manager, UInt order); virtual ~IntegrationScheme() {} + /* ------------------------------------------------------------------------ */ + /* Methods */ + /* ------------------------------------------------------------------------ */ +public: + /// generic interface of a predictor + virtual void predictor(const ID & dof_id, Real delta_t) = 0; + + /// generic interface of a corrector + virtual void corrector(const SolutionType & type, const ID & dof_id, + Real delta_t) = 0; + + /// assemble the jacobian matrix + virtual void assembleJacobian(const SolutionType & type, Real delta_t) = 0; + /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the order of the integration scheme UInt getOrder(); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: + DOFManager & dof_manager; + /// The order of the integrator UInt order; }; __END_AKANTU__ #endif /* __AKANTU_INTEGRATION_SCHEME_HH__ */ diff --git a/src/model/non_linear_solver.cc b/src/model/integration_scheme/integration_scheme_1st_order.cc similarity index 52% copy from src/model/non_linear_solver.cc copy to src/model/integration_scheme/integration_scheme_1st_order.cc index 09df956c0..fedfc49c4 100644 --- a/src/model/non_linear_solver.cc +++ b/src/model/integration_scheme/integration_scheme_1st_order.cc @@ -1,78 +1,69 @@ /** - * @file non_linear_solver.cc + * @file integration_scheme_1st_order.cc * * @author Nicolas Richart * - * @date Tue Oct 13 15:34:43 2015 + * @date Fri Oct 23 12:31:32 2015 * - * @brief Implementation of the base class NonLinearSolver + * @brief Implementation of the common functions for 1st order time + *integrations * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ -#include "non_linear_solver.hh" +#include "integration_scheme_1st_order.hh" +#include "dof_manager.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ -void NonLinearSolver::callbackPredictors() { - NonLinearSolverCallbackMap it = solver_callbacks.begin(); - NonLinearSolverCallbackMap end = solver_callbacks.end(); +void IntegrationScheme1stOrder::predictor(const ID & dof_id, Real delta_t) { + AKANTU_DEBUG_IN(); - for (; it != end; ++it) { - it->second->predictor(); - } -} + Array & u = this->dof_manager.getDOFs(dof_id); + Array & u_dot = this->dof_manager.getDOFsDerivatives(dof_id, 1); + const Array & blocked_dofs = this->dof_manager.getBlockedDOFs(dof_id); -/* -------------------------------------------------------------------------- */ -void NonLinearSolver::callbackCorrectors() { - NonLinearSolverCallbackMap it = solver_callbacks.begin(); - NonLinearSolverCallbackMap end = solver_callbacks.end(); + this->predictor(delta_t, u, u_dot, blocked_dofs); - for (; it != end; ++it) { - it->second->corrector(); - } + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void NonLinearSolver::callbackAssembleJacobian() { - NonLinearSolverCallbackMap it = solver_callbacks.begin(); - NonLinearSolverCallbackMap end = solver_callbacks.end(); +void IntegrationScheme1stOrder::corrector(const SolutionType & type, + const ID & dof_id, Real delta_t) { + AKANTU_DEBUG_IN(); - for (; it != end; ++it) { - it->second->assembleJacobian(); - } -} + Array & u = this->dof_manager.getDOFs(dof_id); + Array & u_dot = this->dof_manager.getDOFsDerivatives(dof_id, 1); -/* -------------------------------------------------------------------------- */ -void NonLinearSolver::callbackAssembleResidual() { - NonLinearSolverCallbackMap it = solver_callbacks.begin(); - NonLinearSolverCallbackMap end = solver_callbacks.end(); + const Array & solution = this->dof_manager.getSolution(dof_id); + const Array & blocked_dofs = this->dof_manager.getBlockedDOFs(dof_id); + + this->corrector(type, delta_t, u, u_dot, blocked_dofs, solution); - for (; it != end; ++it) { - it->second->assembleResidual(); - } + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/integration_scheme/integration_scheme_1st_order.hh b/src/model/integration_scheme/integration_scheme_1st_order.hh index eaba46c09..a3d6635d6 100644 --- a/src/model/integration_scheme/integration_scheme_1st_order.hh +++ b/src/model/integration_scheme/integration_scheme_1st_order.hh @@ -1,88 +1,87 @@ /** * @file integration_scheme_1st_order.hh * * @author Nicolas Richart * * @date creation: Mon Jul 04 2011 * @date last modification: Wed Mar 13 2013 * * @brief Interface of the time integrator of first order * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "integration_scheme.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_INTEGRATION_SCHEME_1ST_ORDER_HH__ #define __AKANTU_INTEGRATION_SCHEME_1ST_ORDER_HH__ __BEGIN_AKANTU__ class IntegrationScheme1stOrder : public IntegrationScheme { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - enum IntegrationSchemeCorrectorType { - _temperature_corrector, - _temperature_rate_corrector - }; + IntegrationScheme1stOrder(DOFManager & dof_manager) + : IntegrationScheme(dof_manager, 1){}; - IntegrationScheme1stOrder() : IntegrationScheme(1){ }; - - virtual ~IntegrationScheme1stOrder() {}; + virtual ~IntegrationScheme1stOrder(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: - virtual void integrationSchemePred(Real delta_t, - Array & u, - Array & u_dot, - const Array & boundary) const = 0; + virtual void predictor(const ID & dof_id, Real delta_t); + virtual void corrector(const SolutionType & type, const ID & dof_id, Real delta_t); - virtual void integrationSchemeCorrTemp(Real delta_t, Array & u, - Array & u_dot, - const Array & boundary, - const Array & delta) const = 0; + /// generic interface of a predictor of 1st order + virtual void predictor(Real delta_t, Array & u, Array & u_dot, + const Array & boundary) const = 0; - virtual void integrationSchemeCorrTempRate(Real delta_t, Array & u, - Array & u_dot, - const Array & boundary, - const Array & delta) const = 0; + /// generic interface of a corrector of 1st order + virtual void corrector(const SolutionType & type, Real delta_t, + Array & u, Array & u_dot, + const Array & boundary, + const Array & delta) const = 0; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ -public: +protected: + virtual Real getTemperatureCoefficient(const SolutionType & type, + Real delta_t) const = 0; + virtual Real getTemperatureRateCoefficient(const SolutionType & type, + Real delta_t) const = 0; + /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: }; __END_AKANTU__ #include "generalized_trapezoidal.hh" #endif /* __AKANTU_INTEGRATION_SCHEME_1ST_ORDER_HH__ */ diff --git a/src/model/integration_scheme/integration_scheme_2nd_order.cc b/src/model/integration_scheme/integration_scheme_2nd_order.cc new file mode 100644 index 000000000..1048933af --- /dev/null +++ b/src/model/integration_scheme/integration_scheme_2nd_order.cc @@ -0,0 +1,72 @@ +/** + * @file integration_scheme_2nd_order.cc + * + * @author Nicolas Richart + * + * @date Tue Oct 20 10:41:12 2015 + * + * @brief Implementation of the common part of 2nd order integration schemes + * + * @section LICENSE + * + * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) any + * later version. + * + * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR + * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with Akantu. If not, see . + * + */ + +/* -------------------------------------------------------------------------- */ +#include "integration_scheme_2nd_order.hh" +#include "dof_manager.hh" +/* -------------------------------------------------------------------------- */ + +__BEGIN_AKANTU__ + +/* -------------------------------------------------------------------------- */ +void IntegrationScheme2ndOrder::predictor(const ID & dof_id, + Real delta_t) { + AKANTU_DEBUG_IN(); + + Array & u = this->dof_manager.getDOFs(dof_id); + Array & u_dot = this->dof_manager.getDOFsDerivatives(dof_id, 1); + Array & u_dot_dot = this->dof_manager.getDOFsDerivatives(dof_id, 2); + const Array & blocked_dofs = this->dof_manager.getBlockedDOFs(dof_id); + + this->predictor(delta_t, u, u_dot, u_dot_dot, blocked_dofs); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void IntegrationScheme2ndOrder::corrector(const SolutionType & type, + const ID & dof_id, + Real delta_t) { + AKANTU_DEBUG_IN(); + + Array & u = this->dof_manager.getDOFs(dof_id); + Array & u_dot = this->dof_manager.getDOFsDerivatives(dof_id, 1); + Array & u_dot_dot = this->dof_manager.getDOFsDerivatives(dof_id, 2); + + const Array & solution = this->dof_manager.getSolution(dof_id); + const Array & blocked_dofs = this->dof_manager.getBlockedDOFs(dof_id); + + this->corrector(type, delta_t, u, u_dot, u_dot_dot, blocked_dofs, solution); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ + +__END_AKANTU__ diff --git a/src/model/integration_scheme/integration_scheme_2nd_order.hh b/src/model/integration_scheme/integration_scheme_2nd_order.hh index 46aae47c2..766ee1c8a 100644 --- a/src/model/integration_scheme/integration_scheme_2nd_order.hh +++ b/src/model/integration_scheme/integration_scheme_2nd_order.hh @@ -1,100 +1,103 @@ /** * @file integration_scheme_2nd_order.hh * * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Wed Aug 04 2010 * @date last modification: Thu Jun 05 2014 * * @brief Interface of the integrator of second order * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ +#include "aka_array.hh" +#include "integration_scheme.hh" +/* -------------------------------------------------------------------------- */ #ifndef __AKANTU_INTEGRATION_SCHEME_2ND_ORDER_HH__ #define __AKANTU_INTEGRATION_SCHEME_2ND_ORDER_HH__ -/* -------------------------------------------------------------------------- */ -#include "aka_common.hh" -#include "aka_array.hh" -#include "integration_scheme.hh" -/* -------------------------------------------------------------------------- */ +namespace akantu { + class SparseMatrix; +} __BEGIN_AKANTU__ class IntegrationScheme2ndOrder : public IntegrationScheme { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - enum IntegrationSchemeCorrectorType { - _acceleration_corrector, - _velocity_corrector, - _displacement_corrector - }; - - IntegrationScheme2ndOrder() : IntegrationScheme(2){}; + IntegrationScheme2ndOrder(DOFManager & dof_manager) : IntegrationScheme(dof_manager, 2){}; virtual ~IntegrationScheme2ndOrder(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: - virtual void integrationSchemePred(Real delta_t, Array & u, - Array & u_dot, - Array & u_dot_dot, - const Array & blocked_dofs) const = 0; - - virtual void integrationSchemeCorrDispl(Real delta_t, Array & u, - Array & u_dot, - Array & u_dot_dot, - const Array & blocked_dofs, - const Array & delta_u) const = 0; - - virtual void integrationSchemeCorrVeloc(Real delta_t, Array & u, - Array & u_dot, - Array & u_dot_dot, - const Array & blocked_dofs, - const Array & delta_u_dot) const = 0; - - virtual void integrationSchemeCorrAccel(Real delta_t, Array & u, - Array & u_dot, - Array & u_dot_dot, - const Array & blocked_dofs, - const Array & delta_u_dot_dot) const = 0; + /// generic interface of a predictor + virtual void predictor(const ID & dof_id, Real delta_t); + + /// generic interface of a corrector + virtual void corrector(const SolutionType & type, const ID & dof_id, + Real delta_t); + + /// generic interface of a predictor of 2nd order + virtual void predictor(Real delta_t, Array & u, Array & u_dot, + Array & u_dot_dot, + const Array & blocked_dofs) const = 0; + + /// generic interface of a corrector of 2nd order + virtual void corrector(const SolutionType & type, + Real delta_t, Array & u, Array & u_dot, + Array & u_dot_dot, + const Array & blocked_dofs, + const Array & delta) const = 0; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ -public: +protected: + virtual Real + getAccelerationCoefficient(const SolutionType & type, + Real delta_t) const; + + virtual Real + getVelocityCoefficient(const SolutionType & type, + Real delta_t) const; + + virtual Real + getDisplacementCoefficient(const SolutionType & type, + Real delta_t) const; + /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: }; __END_AKANTU__ #include "newmark-beta.hh" #endif /* __AKANTU_INTEGRATION_SCHEME_2ND_ORDER_HH__ */ diff --git a/src/model/integration_scheme/newmark-beta.cc b/src/model/integration_scheme/newmark-beta.cc new file mode 100644 index 000000000..56e675631 --- /dev/null +++ b/src/model/integration_scheme/newmark-beta.cc @@ -0,0 +1,227 @@ +/** + * @file newmark-beta.cc + * + * @author David Simon Kammer + * @author Nicolas Richart + * + * @date creation: Tue Oct 05 2010 + * @date last modification: Thu Jun 05 2014 + * + * @brief implementation of the newmark-@f$\beta@f$ integration scheme. This + * implementation is taken from Méthodes numériques en mécanique des solides by + * Alain Curnier \note{ISBN: 2-88074-247-1} + * + * @section LICENSE + * + * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) any + * later version. + * + * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR + * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with Akantu. If not, see . + * + */ + +/* -------------------------------------------------------------------------- */ +#include "newmark-beta.hh" +#include "dof_manager.hh" +#include "sparse_matrix.hh" +/* -------------------------------------------------------------------------- */ + +__BEGIN_AKANTU__ + +/* + * @f$ \tilde{u_{n+1}} = u_{n} + \Delta t \dot{u}_n + \frac{\Delta t^2}{2} + * \ddot{u}_n @f$ + * @f$ \tilde{\dot{u}_{n+1}} = \dot{u}_{n} + \Delta t \ddot{u}_{n} @f$ + * @f$ \tilde{\ddot{u}_{n}} = \ddot{u}_{n} @f$ + */ +void NewmarkBeta::predictor(Real delta_t, Array & u, Array & u_dot, + Array & u_dot_dot, + const Array & blocked_dofs) const { + AKANTU_DEBUG_IN(); + + UInt nb_nodes = u.getSize(); + UInt nb_degree_of_freedom = u.getNbComponent() * nb_nodes; + + Real * u_val = u.storage(); + Real * u_dot_val = u_dot.storage(); + Real * u_dot_dot_val = u_dot_dot.storage(); + bool * blocked_dofs_val = blocked_dofs.storage(); + + for (UInt d = 0; d < nb_degree_of_freedom; d++) { + if (!(*blocked_dofs_val)) { + Real dt_a_n = delta_t * *u_dot_dot_val; + + *u_val += (1 - k * alpha) * delta_t * *u_dot_val + + (.5 - h * alpha * beta) * delta_t * dt_a_n; + *u_dot_val = (1 - k) * *u_dot_val + (1 - h * beta) * dt_a_n; + *u_dot_dot_val = (1 - h) * *u_dot_dot_val; + } + u_val++; + u_dot_val++; + u_dot_dot_val++; + blocked_dofs_val++; + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void NewmarkBeta::corrector(const SolutionType & type, Real delta_t, + Array & u, Array & u_dot, + Array & u_dot_dot, + const Array & blocked_dofs, + const Array & delta) const { + AKANTU_DEBUG_IN(); + + switch (type) { + case _acceleration: { + this->allCorrector<_acceleration>(delta_t, u, u_dot, u_dot_dot, + blocked_dofs, delta); + break; + } + case _velocity: { + this->allCorrector<_velocity>(delta_t, u, u_dot, u_dot_dot, blocked_dofs, + delta); + break; + } + case _displacement: { + this->allCorrector<_displacement>(delta_t, u, u_dot, u_dot_dot, + blocked_dofs, delta); + break; + } + default: + AKANTU_EXCEPTION("The corrector type : " + << type + << " is not supported by this type of integration scheme"); + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +Real NewmarkBeta::getAccelerationCoefficient(const SolutionType & type, + Real delta_t) const { + switch (type) { + case _acceleration: + return 1.; + case _velocity: + return 1. / (beta * delta_t); + case _displacement: + return 1. / (alpha * beta * delta_t * delta_t); + default: + AKANTU_EXCEPTION("The corrector type : " + << type + << " is not supported by this type of integration scheme"); + } +} + +/* -------------------------------------------------------------------------- */ +Real NewmarkBeta::getVelocityCoefficient(const SolutionType & type, + Real delta_t) const { + switch (type) { + case _acceleration: + return beta * delta_t; + case _velocity: + return 1.; + case _displacement: + return 1. / (alpha * delta_t); + default: + AKANTU_EXCEPTION("The corrector type : " + << type + << " is not supported by this type of integration scheme"); + } +} + +/* -------------------------------------------------------------------------- */ +Real NewmarkBeta::getDisplacementCoefficient(const SolutionType & type, + Real delta_t) const { + switch (type) { + case _acceleration: + return alpha * beta * delta_t * delta_t; + case _velocity: + return alpha * delta_t; + case _displacement: + return 1.; + default: + AKANTU_EXCEPTION("The corrector type : " + << type + << " is not supported by this type of integration scheme"); + } +} + +/* -------------------------------------------------------------------------- */ +template +void NewmarkBeta::allCorrector(Real delta_t, Array & u, + Array & u_dot, Array & u_dot_dot, + const Array & blocked_dofs, + const Array & delta) const { + AKANTU_DEBUG_IN(); + + UInt nb_nodes = u.getSize(); + UInt nb_degree_of_freedom = u.getNbComponent() * nb_nodes; + + Real c = getAccelerationCoefficient(type, delta_t); + Real d = getVelocityCoefficient(type, delta_t); + Real e = getDisplacementCoefficient(type, delta_t); + + Real * u_val = u.storage(); + Real * u_dot_val = u_dot.storage(); + Real * u_dot_dot_val = u_dot_dot.storage(); + Real * delta_val = delta.storage(); + bool * blocked_dofs_val = blocked_dofs.storage(); + + for (UInt dof = 0; dof < nb_degree_of_freedom; dof++) { + if (!(*blocked_dofs_val)) { + *u_val += e * *delta_val; + *u_dot_val += d * *delta_val; + *u_dot_dot_val += c * *delta_val; + } + u_val++; + u_dot_val++; + u_dot_dot_val++; + delta_val++; + blocked_dofs_val++; + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void NewmarkBeta::assembleJacobian(const SolutionType & type, + Real delta_t) { + AKANTU_DEBUG_IN(); + + SparseMatrix & J = this->dof_manager.getMatrix("J"); + + const SparseMatrix & M = this->dof_manager.getMatrix("M"); + const SparseMatrix & K = this->dof_manager.getMatrix("K"); + + Real c = this->getAccelerationCoefficient(type, delta_t); + Real e = this->getDisplacementCoefficient(type, delta_t); + + J.add(K, e); + J.add(M, c); + + try { + Real d = this->getVelocityCoefficient(type, delta_t); + const SparseMatrix & C = this->dof_manager.getMatrix("C"); + J.add(C, d); + } catch (debug::Exception &) { + } + + AKANTU_DEBUG_OUT(); +} +/* -------------------------------------------------------------------------- */ + +__END_AKANTU__ diff --git a/src/model/integration_scheme/newmark-beta.hh b/src/model/integration_scheme/newmark-beta.hh index 5bc9f6875..8d4771a87 100644 --- a/src/model/integration_scheme/newmark-beta.hh +++ b/src/model/integration_scheme/newmark-beta.hh @@ -1,200 +1,186 @@ /** * @file newmark-beta.hh * * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Tue Oct 05 2010 * @date last modification: Thu Jun 05 2014 * * @brief implementation of the newmark-@f$\beta@f$ integration scheme. This * implementation is taken from Méthodes numériques en mécanique des solides by * Alain Curnier \note{ISBN: 2-88074-247-1} * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_NEWMARK_BETA_HH__ #define __AKANTU_NEWMARK_BETA_HH__ /* -------------------------------------------------------------------------- */ #include "integration_scheme_2nd_order.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /** * The three differentiate equations (dynamic and cinematic) are : * @f{eqnarray*}{ * M \ddot{u}_{n+1} + C \dot{u}_{n+1} + K u_{n+1} &=& q_{n+1} \\ * u_{n+1} &=& u_{n} + (1 - \alpha) \Delta t \dot{u}_{n} + \alpha \Delta t *\dot{u}_{n+1} + (1/2 - \alpha) \Delta t^2 \ddot{u}_n \\ * \dot{u}_{n+1} &=& \dot{u}_{n} + (1 - \beta) \Delta t \ddot{u}_{n} + \beta *\Delta t \ddot{u}_{n+1} * @f} * * Predictor: * @f{eqnarray*}{ * u^{0}_{n+1} &=& u_{n} + \Delta t \dot{u}_n + \frac{\Delta t^2}{2} *\ddot{u}_n \\ * \dot{u}^{0}_{n+1} &=& \dot{u}_{n} + \Delta t \ddot{u}_{n} \\ * \ddot{u}^{0}_{n+1} &=& \ddot{u}_{n} * @f} * * Solve : * @f[ (c M + d C + e K^i_{n+1}) w = = q_{n+1} - f^i_{n+1} - C \dot{u}^i_{n+1} *- M \ddot{u}^i_{n+1} @f] * * Corrector : * @f{eqnarray*}{ * \ddot{u}^{i+1}_{n+1} &=& \ddot{u}^{i}_{n+1} + c w \\ * \dot{u}^{i+1}_{n+1} &=& \dot{u}^{i}_{n+1} + d w \\ * u^{i+1}_{n+1} &=& u^{i}_{n+1} + e w * @f} * * c, d and e are parameters depending on the method used to solve the equations *@n * For acceleration : @f$ w = \delta \ddot{u}, e = \alpha \beta \Delta t^2, d = *\beta \Delta t, c = 1 @f$ @n * For velocity : @f$ w = \delta \dot{u}, e = 1/\beta \Delta t, d = *1, c = \alpha \Delta t @f$ @n * For displacement : @f$ w = \delta u, e = 1, d = *1/\alpha \Delta t, c = 1/\alpha \beta \Delta t^2 @f$ */ class NewmarkBeta : public IntegrationScheme2ndOrder { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - NewmarkBeta(Real alpha, Real beta) : beta(beta), alpha(alpha), k(0.), h(0.){}; + NewmarkBeta(DOFManager & dof_manager, Real alpha, Real beta) + : IntegrationScheme2ndOrder(dof_manager), beta(beta), alpha(alpha), k(0.), + h(0.){}; ~NewmarkBeta(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: - inline void integrationSchemePred(Real delta_t, Array & u, - Array & u_dot, - Array & u_dot_dot, - const Array & blocked_dofs) const; - - inline void integrationSchemeCorrAccel(Real delta_t, Array & u, - Array & u_dot, - Array & u_dot_dot, - const Array & blocked_dofs, - const Array & delta) const; - - inline void integrationSchemeCorrVeloc(Real delta_t, Array & u, - Array & u_dot, - Array & u_dot_dot, - const Array & blocked_dofs, - const Array & delta) const; - - inline void integrationSchemeCorrDispl(Real delta_t, Array & u, - Array & u_dot, - Array & u_dot_dot, - const Array & blocked_dofs, - const Array & delta) const; + void predictor(Real delta_t, Array & u, Array & u_dot, + Array & u_dot_dot, + const Array & blocked_dofs) const; + + void corrector(const SolutionType & type, Real delta_t, Array & u, + Array & u_dot, Array & u_dot_dot, + const Array & blocked_dofs, + const Array & delta) const; + + void assembleJacobian(const SolutionType & type, Real delta_t); public: - template - Real getAccelerationCoefficient(Real delta_t) const; + Real getAccelerationCoefficient(const SolutionType & type, + Real delta_t) const; - template - Real getVelocityCoefficient(Real delta_t) const; + Real getVelocityCoefficient(const SolutionType & type, Real delta_t) const; - template - Real getDisplacementCoefficient(Real delta_t) const; + Real getDisplacementCoefficient(const SolutionType & type, + Real delta_t) const; private: - template - void integrationSchemeCorr(Real delta_t, Array & u, Array & u_dot, - Array & u_dot_dot, - const Array & blocked_dofs, - const Array & delta) const; + template + void allCorrector(Real delta_t, Array & u, Array & u_dot, + Array & u_dot_dot, const Array & blocked_dofs, + const Array & delta) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Beta, beta, Real); AKANTU_GET_MACRO(Alpha, alpha, Real); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// the @f$\beta@f$ parameter const Real beta; /// the @f$\alpha@f$ parameter const Real alpha; const Real k; const Real h; }; -/* -------------------------------------------------------------------------- */ -/* inline functions */ -/* -------------------------------------------------------------------------- */ - -#if defined(AKANTU_INCLUDE_INLINE_IMPL) -#include "newmark-beta_inline_impl.cc" -#endif - /** * central difference method (explicit) * undamped stability condition : * @f$ \Delta t = \alpha \Delta t_{crit} = \frac{2}{\omega_{max}} \leq \min_{e} *\frac{l_e}{c_e} * */ class CentralDifference : public NewmarkBeta { public: - CentralDifference() : NewmarkBeta(0., 1./2.){}; + CentralDifference(DOFManager & dof_manager) + : NewmarkBeta(dof_manager, 0., 1. / 2.){}; }; //#include "integration_scheme/central_difference.hh" /// undamped trapezoidal rule (implicit) class TrapezoidalRule2 : public NewmarkBeta { public: - TrapezoidalRule2() : NewmarkBeta(1./2., 1./2.){}; + TrapezoidalRule2(DOFManager & dof_manager) + : NewmarkBeta(dof_manager, 1. / 2., 1. / 2.){}; }; /// Fox-Goodwin rule (implicit) class FoxGoodwin : public NewmarkBeta { public: - FoxGoodwin() : NewmarkBeta(1./6., 1./2.){}; + FoxGoodwin(DOFManager & dof_manager) + : NewmarkBeta(dof_manager, 1. / 6., 1. / 2.){}; }; /// Linear acceleration (implicit) class LinearAceleration : public NewmarkBeta { public: - LinearAceleration() : NewmarkBeta(1./3., 1./2.){}; + LinearAceleration(DOFManager & dof_manager) + : NewmarkBeta(dof_manager, 1. / 3., 1. / 2.){}; }; +/* -------------------------------------------------------------------------- */ + __END_AKANTU__ #endif /* __AKANTU_NEWMARK_BETA_HH__ */ diff --git a/src/model/integration_scheme/newmark-beta_inline_impl.cc b/src/model/integration_scheme/newmark-beta_inline_impl.cc deleted file mode 100644 index 98eb1b432..000000000 --- a/src/model/integration_scheme/newmark-beta_inline_impl.cc +++ /dev/null @@ -1,230 +0,0 @@ -/** - * @file newmark-beta_inline_impl.cc - * - * @author David Simon Kammer - * @author Nicolas Richart - * - * @date creation: Tue Oct 05 2010 - * @date last modification: Thu Jun 05 2014 - * - * @brief implementation of the newmark-@f$\beta@f$ integration scheme. This - * implementation is taken from Méthodes numériques en mécanique des solides by - * Alain Curnier \note{ISBN: 2-88074-247-1} - * - * @section LICENSE - * - * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) - * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) - * - * Akantu is free software: you can redistribute it and/or modify it under the - * terms of the GNU Lesser General Public License as published by the Free - * Software Foundation, either version 3 of the License, or (at your option) any - * later version. - * - * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY - * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR - * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more - * details. - * - * You should have received a copy of the GNU Lesser General Public License - * along with Akantu. If not, see . - * - */ - -/* -------------------------------------------------------------------------- */ -/* - * @f$ \tilde{u_{n+1}} = u_{n} + \Delta t \dot{u}_n + \frac{\Delta t^2}{2} \ddot{u}_n @f$ - * @f$ \tilde{\dot{u}_{n+1}} = \dot{u}_{n} + \Delta t \ddot{u}_{n} @f$ - * @f$ \tilde{\ddot{u}_{n}} = \ddot{u}_{n} @f$ - */ -inline void NewmarkBeta::integrationSchemePred(Real delta_t, - Array & u, - Array & u_dot, - Array & u_dot_dot, - const Array & blocked_dofs) const { - AKANTU_DEBUG_IN(); - - UInt nb_nodes = u.getSize(); - UInt nb_degree_of_freedom = u.getNbComponent() * nb_nodes; - - Real * u_val = u.storage(); - Real * u_dot_val = u_dot.storage(); - Real * u_dot_dot_val = u_dot_dot.storage(); - bool * blocked_dofs_val = blocked_dofs.storage(); - - for (UInt d = 0; d < nb_degree_of_freedom; d++) { - if(!(*blocked_dofs_val)) { - Real dt_a_n = delta_t * *u_dot_dot_val; - - *u_val += (1 - k*alpha) * delta_t * *u_dot_val + (.5 - h*alpha*beta) * delta_t * dt_a_n; - *u_dot_val = (1 - k) * *u_dot_val + (1 - h*beta) * dt_a_n; - *u_dot_dot_val = (1 - h) * *u_dot_dot_val; - } - u_val++; - u_dot_val++; - u_dot_dot_val++; - blocked_dofs_val++; - } - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -/* - * @f$ u_{n+1} = \tilde{u_{n+1}} + \alpha \beta \Delta t^2 \delta \ddot{u}_n @f$ - * @f$ \dot{u}_{n+1} = \tilde{\dot{u}_{n+1}} + \beta \Delta t * \delta \ddot{u}_{n+1} @f$ - * @f$ \ddot{u}_{n+1} = \tilde{\ddot{u}_{n+1}} + \delta \ddot{u}_{n+1} @f$ - */ -inline void NewmarkBeta::integrationSchemeCorrAccel(Real delta_t, - Array & u, - Array & u_dot, - Array & u_dot_dot, - const Array & blocked_dofs, - const Array & delta - ) const { - AKANTU_DEBUG_IN(); - - integrationSchemeCorr<_acceleration_corrector>(delta_t, - u, - u_dot, - u_dot_dot, - blocked_dofs, - delta); - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -/* - * @f$ u_{n+1} = \tilde{u_{n+1}} + \alpha \Delta t \delta \dot{u}_n @f$ - * @f$ \dot{u}_{n+1} = \tilde{\dot{u}_{n+1}} + \delta \dot{u}_{n+1} @f$ - * @f$ \ddot{u}_{n+1} = \tilde{\ddot{u}_{n+1}} + \frac{1}{\beta \Delta t} \delta \dot{u}_{n+1} @f$ - */ -inline void NewmarkBeta::integrationSchemeCorrVeloc(Real delta_t, - Array & u, - Array & u_dot, - Array & u_dot_dot, - const Array & blocked_dofs, - const Array & delta - ) const { - AKANTU_DEBUG_IN(); - - integrationSchemeCorr<_velocity_corrector>(delta_t, - u, - u_dot, - u_dot_dot, - blocked_dofs, - delta); - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -/* - * @f$ u_{n+1} = \tilde{u_{n+1}} + \delta u_n @f$ - * @f$ \dot{u}_{n+1} = \tilde{\dot{u}_{n+1}} + \frac{1}{\alpha \Delta t} \delta u_{n+1} @f$ - * @f$ \ddot{u}_{n+1} = \tilde{\ddot{u}_{n+1}} + \frac{1}{\alpha \beta \Delta t^2} \delta u_{n+1} @f$ - */ -inline void NewmarkBeta::integrationSchemeCorrDispl(Real delta_t, - Array & u, - Array & u_dot, - Array & u_dot_dot, - const Array & blocked_dofs, - const Array & delta - ) const { - AKANTU_DEBUG_IN(); - - integrationSchemeCorr<_displacement_corrector>(delta_t, - u, - u_dot, - u_dot_dot, - blocked_dofs, - delta); - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -template<> -inline Real NewmarkBeta::getAccelerationCoefficient(__attribute__((unused)) Real delta_t) const { - return 1.; -} -template<> -inline Real NewmarkBeta::getAccelerationCoefficient(Real delta_t) const { - return 1. / (beta * delta_t); -} -template<> -inline Real NewmarkBeta::getAccelerationCoefficient(Real delta_t) const { - return 1. / (alpha * beta * delta_t * delta_t); -} - -/* -------------------------------------------------------------------------- */ -template<> -inline Real NewmarkBeta::getVelocityCoefficient(Real delta_t) const { - return beta * delta_t; -} -template<> -inline Real NewmarkBeta::getVelocityCoefficient(__attribute__((unused)) Real delta_t) const { - return 1.; -} -template<> -inline Real NewmarkBeta::getVelocityCoefficient(Real delta_t) const { - return 1. / (alpha * delta_t); -} - -/* -------------------------------------------------------------------------- */ -template<> -inline Real NewmarkBeta::getDisplacementCoefficient(Real delta_t) const { - return alpha * beta * delta_t * delta_t; -} -template<> -inline Real NewmarkBeta::getDisplacementCoefficient(Real delta_t) const { - return alpha * delta_t; -} -template<> -inline Real NewmarkBeta::getDisplacementCoefficient(__attribute__((unused)) Real delta_t) const { - return 1.; -} - - - - -/* -------------------------------------------------------------------------- */ -template -void NewmarkBeta::integrationSchemeCorr(Real delta_t, - Array & u, - Array & u_dot, - Array & u_dot_dot, - const Array & blocked_dofs, - const Array & delta - ) const { - AKANTU_DEBUG_IN(); - - UInt nb_nodes = u.getSize(); - UInt nb_degree_of_freedom = u.getNbComponent() * nb_nodes; - - Real c = getAccelerationCoefficient(delta_t); - Real d = getVelocityCoefficient(delta_t); - Real e = getDisplacementCoefficient(delta_t); - - Real * u_val = u.storage(); - Real * u_dot_val = u_dot.storage(); - Real * u_dot_dot_val = u_dot_dot.storage(); - Real * delta_val = delta.storage(); - bool * blocked_dofs_val = blocked_dofs.storage(); - - for (UInt dof = 0; dof < nb_degree_of_freedom; dof++) { - if(!(*blocked_dofs_val)) { - *u_val += e * *delta_val; - *u_dot_val += d * *delta_val; - *u_dot_dot_val += c * *delta_val; - } - u_val++; - u_dot_val++; - u_dot_dot_val++; - delta_val++; - blocked_dofs_val++; - } - - AKANTU_DEBUG_OUT(); -} diff --git a/src/model/model_solver.cc b/src/model/model_solver.cc index 0d5551a25..c0c3859ce 100644 --- a/src/model/model_solver.cc +++ b/src/model/model_solver.cc @@ -1,384 +1,437 @@ /** * @file model_solver.cc * * @author Nicolas Richart * * @date Wed Aug 12 13:31:56 2015 * * @brief Implementation of ModelSolver * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "model_solver.hh" #include "dof_manager.hh" #if defined(AKANTU_USE_MPI) #include "mpi_type_wrapper.hh" #endif #if defined(AKANTU_USE_MUMPS) #include "dof_manager_default.hh" #endif #if defined(AKANTU_USE_PETSC) #include "dof_manager_petsc.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ ModelSolver::ModelSolver(const Mesh & mesh, const ID & id, UInt memory_id) : Parsable(_st_solver, id), dof_manager(NULL) { std::pair sub_sect = getStaticParser().getSubSections(_st_solver); if (sub_sect.first != sub_sect.second) { AKANTU_EXCEPTION("More than on solver section present in the input file"); } const ParserSection & section = *sub_sect.first; std::string solver_type = section.getName(); if (solver_type == "petsc") { #if defined(AKANTU_USE_PETSC) - this->dof_manager = new DOFManagerPETSc(mesh, id + ":dof_manager_petsc", memory_id); + this->dof_manager = + new DOFManagerPETSc(mesh, id + ":dof_manager_petsc", memory_id); #else AKANTU_EXCEPTION( "To use PETSc you have to activate it in the compilations options"); #endif } else if (solver_type == "mumps") { #if defined(AKANTU_USE_MUMPS) this->dof_manager = new DOFManagerDefault(mesh, id + ":dof_manager_default", memory_id); #else AKANTU_EXCEPTION( "To use MUMPS you have to activate it in the compilations options"); #endif } else { AKANTU_EXCEPTION( "To use the solver " << solver_type << " you will have to code it. This is an unknown solver type."); } } +/* -------------------------------------------------------------------------- */ +void ModelSolver::solveStep() { + AKANTU_DEBUG_IN(); + + this->solveStep(this->default_time_step_solver_id, + this->default_non_linear_solver_id); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void ModelSolver::solveStep(const ID & time_step_solver_id, + const ID & non_linear_solver_id) { + AKANTU_DEBUG_IN(); + + std::set::const_iterator tss_it = + time_step_solvers.find(time_step_solver_id); + if (tss_it == time_step_solvers.end()) { + AKANTU_EXCEPTION("No time step solver was instantiated with this id"); + } + + std::set::const_iterator nls_it = + time_step_solvers.find(time_step_solver_id); + if (tss_it == time_step_solvers.end()) { + AKANTU_EXCEPTION("No non linear solver was instantiated with this id"); + } + + NonLinearSolver & nls = + this->dof_manager->getNonLinearSolver(non_linear_solver_id); + TimeStepSolver & tss = + this->dof_manager->getTimeStepSolver(time_step_solver_id); + + NonLinearSolverCallback & previous_nls_callback = nls.getCallbacks(); + NonLinearSolverCallback & previous_tss_callback = nls.getCallbacks(); + + // set the callbacks of the time step solver to the current model + // solver + tss.registerCallback(*this); + + // set the callbacks of the non linear solver to the chosen time step solver + nls.registerCallback(tss); + + // make one non linear solve + nls.solve(); + + /// restores the callbacks + nls.registerCallback(previous_nls_callback); + tss.registerCallback(previous_tss_callback); + + AKANTU_DEBUG_OUT(); +} + /* -------------------------------------------------------------------------- */ void ModelSolver::initTimeStepSolver( const ID & time_step_solver_id, const ID & dof_id, const TimeStepSolverType & time_step_solver_type) { - if(this->default_time_step_solver == "") { - this->default_time_step_solver = time_step_solver_id; + if (this->default_time_step_solver_id == "") { + this->default_time_step_solver_id = time_step_solver_id; } this->dof_manager->getNewTimeStepSolver(time_step_solver_id, dof_id, time_step_solver_type); } /* -------------------------------------------------------------------------- */ // /* -------------------------------------------------------------------------- // */ // template // void SolidMechanicsModel::solve(Array &increment, Real block_val, // bool need_factorize, bool // has_profile_changed) { // if(has_profile_changed) { // this->initJacobianMatrix(); // need_factorize = true; // } // updateResidualInternal(); //doesn't do anything for static // if(need_factorize) { // Real c = 0.,d = 0.,e = 0.; // if(method == _static) { // AKANTU_DEBUG_INFO("Solving K inc = r"); // e = 1.; // } else { // AKANTU_DEBUG_INFO("Solving (c M + d C + e K) inc = r"); // NewmarkBeta * nmb_int = dynamic_cast(integrator); // c = nmb_int->getAccelerationCoefficient(time_step); // d = nmb_int->getVelocityCoefficient(time_step); // e = nmb_int->getDisplacementCoefficient(time_step); // } // jacobian_matrix->clear(); // // J = c M + d C + e K // if(stiffness_matrix) // jacobian_matrix->add(*stiffness_matrix, e); // if(mass_matrix) // jacobian_matrix->add(*mass_matrix, c); // #if !defined(AKANTU_NDEBUG) // if(mass_matrix && AKANTU_DEBUG_TEST(dblDump)) // mass_matrix->saveMatrix("M.mtx"); // #endif // if(velocity_damping_matrix) // jacobian_matrix->add(*velocity_damping_matrix, d); // jacobian_matrix->applyBoundary(*blocked_dofs, block_val); // #if !defined(AKANTU_NDEBUG) // if(AKANTU_DEBUG_TEST(dblDump)) // jacobian_matrix->saveMatrix("J.mtx"); // #endif // solver->factorize(); // } // // if (rhs.getSize() != 0) // // solver->setRHS(rhs); // // else // solver->setOperators(); // solver->setRHS(*residual); // // solve @f[ J \delta w = r @f] // solver->solve(increment); // UInt nb_nodes = displacement-> getSize(); // UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes; // bool * blocked_dofs_val = blocked_dofs->storage(); // Real * increment_val = increment.storage(); // for (UInt j = 0; j < nb_degree_of_freedom; // ++j,++increment_val, ++blocked_dofs_val) { // if ((*blocked_dofs_val)) // *increment_val = 0.0; // } // } // /* -------------------------------------------------------------------------- // */ // template // bool SolidMechanicsModel::solveStatic(Real tolerance, UInt max_iteration, // bool do_not_factorize) { // AKANTU_DEBUG_INFO("Solving Ku = f"); // AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL, // "You should first initialize the implicit solver and // assemble the stiffness matrix by calling // initImplicit"); // AnalysisMethod analysis_method=method; // Real error = 0.; // method=_static; // bool converged = this->template solveStep(tolerance, // error, max_iteration, do_not_factorize); // method=analysis_method; // return converged; // } // /* -------------------------------------------------------------------------- // */ // template // bool SolidMechanicsModel::solveStep(Real tolerance, // UInt max_iteration) { // Real error = 0.; // return this->template solveStep(tolerance, // error, // max_iteration); // } // /* -------------------------------------------------------------------------- // */ // template // bool SolidMechanicsModel::solveStep(Real tolerance, Real & error, UInt // max_iteration, // bool do_not_factorize) { // EventManager::sendEvent(SolidMechanicsModelEvent::BeforeSolveStepEvent(method)); // this->implicitPred(); // this->updateResidual(); // AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL, // "You should first initialize the implicit solver and // assemble the stiffness matrix"); // bool need_factorize = !do_not_factorize; // if (method==_implicit_dynamic) { // AKANTU_DEBUG_ASSERT(mass_matrix != NULL, // "You should first initialize the implicit solver and // assemble the mass matrix"); // } // switch (cmethod) { // case _scm_newton_raphson_tangent: // case _scm_newton_raphson_tangent_not_computed: // break; // case _scm_newton_raphson_tangent_modified: // this->assembleStiffnessMatrix(); // break; // default: // AKANTU_DEBUG_ERROR("The resolution method " << cmethod << " has not been // implemented!"); // } // this->n_iter = 0; // bool converged = false; // error = 0.; // if(criteria == _scc_residual) { // converged = this->testConvergence (tolerance, error); // if(converged) return converged; // } // do { // if (cmethod == _scm_newton_raphson_tangent) // this->assembleStiffnessMatrix(); // solve (*increment, 1., // need_factorize); // this->implicitCorr(); // if(criteria == _scc_residual) this->updateResidual(); // converged = this->testConvergence (tolerance, error); // if(criteria == _scc_increment && !converged) this->updateResidual(); // //this->dump(); // this->n_iter++; // AKANTU_DEBUG_INFO("[" << criteria << "] Convergence iteration " // << std::setw(std::log10(max_iteration)) << this->n_iter // << ": error " << error << (converged ? " < " : " > ") // << tolerance); // switch (cmethod) { // case _scm_newton_raphson_tangent: // need_factorize = true; // break; // case _scm_newton_raphson_tangent_not_computed: // case _scm_newton_raphson_tangent_modified: // need_factorize = false; // break; // default: // AKANTU_DEBUG_ERROR("The resolution method " << cmethod << " has not // been implemented!"); // } // } while (!converged && this->n_iter < max_iteration); // // this makes sure that you have correct strains and stresses after the // solveStep function (e.g., for dumping) // if(criteria == _scc_increment) this->updateResidual(); // if (converged) { // EventManager::sendEvent(SolidMechanicsModelEvent::AfterSolveStepEvent(method)); // } else if(this->n_iter == max_iteration) { // AKANTU_DEBUG_WARNING("[" << criteria << "] Convergence not reached after // " // << std::setw(std::log10(max_iteration)) << // this->n_iter << // " iteration" << (this->n_iter == 1 ? "" : "s") << // "!" << std::endl); // } // return converged; // } // void SolidMechanicsModel::updateResidualInternal() { // AKANTU_DEBUG_IN(); // AKANTU_DEBUG_INFO("Update the residual"); // // f = f_ext - f_int - Ma - Cv = r - Ma - Cv; // if(method != _static) { // // f -= Ma // if(mass_matrix) { // // if full mass_matrix // Array * Ma = new Array(*acceleration, true, "Ma"); // *Ma *= *mass_matrix; // /// \todo check unit conversion for implicit dynamics // // *Ma /= f_m2a // *residual -= *Ma; // delete Ma; // } else if (mass) { // // else lumped mass // UInt nb_nodes = acceleration->getSize(); // UInt nb_degree_of_freedom = acceleration->getNbComponent(); // Real * mass_val = mass->storage(); // Real * accel_val = acceleration->storage(); // Real * res_val = residual->storage(); // bool * blocked_dofs_val = blocked_dofs->storage(); // for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) { // if(!(*blocked_dofs_val)) { // *res_val -= *accel_val * *mass_val /f_m2a; // } // blocked_dofs_val++; // res_val++; // mass_val++; // accel_val++; // } // } else { // AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix."); // } // // f -= Cv // if(velocity_damping_matrix) { // Array * Cv = new Array(*velocity); // *Cv *= *velocity_damping_matrix; // *residual -= *Cv; // delete Cv; // } // } // AKANTU_DEBUG_OUT(); // } // /* -------------------------------------------------------------------------- // */ // void SolidMechanicsModel::solveLumped(Array & x, // const Array & A, // const Array & b, // const Array & blocked_dofs, // Real alpha) { // Real * A_val = A.storage(); // Real * b_val = b.storage(); // Real * x_val = x.storage(); // bool * blocked_dofs_val = blocked_dofs.storage(); // UInt nb_degrees_of_freedom = x.getSize() * x.getNbComponent(); // for (UInt n = 0; n < nb_degrees_of_freedom; ++n) { // if(!(*blocked_dofs_val)) { // *x_val = alpha * (*b_val / *A_val); // } // x_val++; // A_val++; // b_val++; // blocked_dofs_val++; // } // } __END_AKANTU__ diff --git a/src/model/model_solver.hh b/src/model/model_solver.hh index f34febd8c..180c55065 100644 --- a/src/model/model_solver.hh +++ b/src/model/model_solver.hh @@ -1,87 +1,113 @@ /** * @file model_solver.hh * * @author Nicolas Richart * * @date Wed Jul 22 10:53:10 2015 * * @brief Class regrouping the common solve interface to the different models * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "parsable.hh" #include "dof_manager.hh" #include "non_linear_solver_callback.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MODEL_SOLVER_HH__ #define __AKANTU_MODEL_SOLVER_HH__ __BEGIN_AKANTU__ class ModelSolver : public Parsable, public NonLinearSolverCallback { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ModelSolver(const Mesh & mesh, const ID & id, UInt memory_id); virtual ~ModelSolver(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// Solve a step, implicit, explicit, static, ... any things based on the /// current configuration void solveStep(); + /// solve a step using a given pre instantiated time step solver and + /// nondynamic linear solver + void solveStep(const ID & time_step_solver_id, + const ID & non_linear_solver_id); + /// Initialize a time solver that can be used afterwards with its id void initTimeStepSolver(const ID & time_step_solver_id, const ID & dof_id, const TimeStepSolverType & time_step_solver_type); -protected: + /* ------------------------------------------------------------------------ */ + /* NonLinearSolverCallback interface */ + /* ------------------------------------------------------------------------ */ +public: + /// Predictor interface for the callback + void predictor(); + + /// Corrector interface for the callback + void corrector(); + + /// AssembleResidual interface for the callback + virtual void assembleResidual() = 0; + + /// AssembleJacobian interface for the callback + virtual void assembleJacobian() = 0; + /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ protected: DOFManager & getDOFManager() { return *this->dof_manager; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// Underlying dof_manager (the brain...) DOFManager * dof_manager; /// List of instantiated time step solvers std::set time_step_solvers; /// Default time step solver to use - ID default_time_step_solver; + ID default_time_step_solver_id; + + /// List of instantiated non_linear_solvers + std::set non_linear_solvers; + + /// Default nondynamic linear solver + ID default_non_linear_solver_id; }; __END_AKANTU__ #endif /* __AKANTU_MODEL_SOLVER_HH__ */ diff --git a/src/model/non_linear_solver.cc b/src/model/non_linear_solver.cc index 09df956c0..195435857 100644 --- a/src/model/non_linear_solver.cc +++ b/src/model/non_linear_solver.cc @@ -1,78 +1,49 @@ /** * @file non_linear_solver.cc * * @author Nicolas Richart * * @date Tue Oct 13 15:34:43 2015 * * @brief Implementation of the base class NonLinearSolver * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "non_linear_solver.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ -void NonLinearSolver::callbackPredictors() { - NonLinearSolverCallbackMap it = solver_callbacks.begin(); - NonLinearSolverCallbackMap end = solver_callbacks.end(); - - for (; it != end; ++it) { - it->second->predictor(); - } -} - -/* -------------------------------------------------------------------------- */ -void NonLinearSolver::callbackCorrectors() { - NonLinearSolverCallbackMap it = solver_callbacks.begin(); - NonLinearSolverCallbackMap end = solver_callbacks.end(); - - for (; it != end; ++it) { - it->second->corrector(); - } -} +NonLinearSolver::NonLinearSolver( + DOFManager & dof_manager, + const NonLinearSolverType & non_linear_solver_type, const ID & id, + UInt memory_id) + : Memory(id, memory_id), _dof_manager(dof_manager), + non_linear_solver_type(non_linear_solver_type) {} /* -------------------------------------------------------------------------- */ -void NonLinearSolver::callbackAssembleJacobian() { - NonLinearSolverCallbackMap it = solver_callbacks.begin(); - NonLinearSolverCallbackMap end = solver_callbacks.end(); - - for (; it != end; ++it) { - it->second->assembleJacobian(); - } -} - -/* -------------------------------------------------------------------------- */ -void NonLinearSolver::callbackAssembleResidual() { - NonLinearSolverCallbackMap it = solver_callbacks.begin(); - NonLinearSolverCallbackMap end = solver_callbacks.end(); - - for (; it != end; ++it) { - it->second->assembleResidual(); - } -} +NonLinearSolver::~NonLinearSolver() {} /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/non_linear_solver.hh b/src/model/non_linear_solver.hh index 2505f0290..55db2c8e9 100644 --- a/src/model/non_linear_solver.hh +++ b/src/model/non_linear_solver.hh @@ -1,100 +1,96 @@ /** * @file non_linear_solver.hh * * @author Nicolas Richart * * @date Mon Aug 24 23:48:41 2015 * * @brief Non linear solver interface * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_NON_LINEAR_SOLVER_HH__ #define __AKANTU_NON_LINEAR_SOLVER_HH__ namespace akantu { class DOFManager; class NonLinearSolverCallback; } __BEGIN_AKANTU__ class NonLinearSolver : Memory { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: NonLinearSolver(DOFManager & dof_manager, const NonLinearSolverType & non_linear_solver_type, const ID & id = "non_linear_solver", UInt memory_id = 0); virtual ~NonLinearSolver(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// solve the system described by the jacobian matrix, and rhs contained in /// the dof manager virtual void solve() = 0; -protected: - /// call all the registered predictors - virtual void callbackPredictors(); - - /// call all the registered correctors - virtual void callbackCorrectors(); + /// register a callback object for a given dof_id + virtual void registerCallback(NonLinearSolverCallback & callbacks); - /// call all the registered predictors - virtual void callbackAssembleJacobian(); + /// register a callback object for a given dof_id + virtual NonLinearSolverCallback & getCallbacks(); - /// call all the registered predictors - virtual void callbackAssembleResidual(); +protected: /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// configure the solver from ParserSection virtual void setParameters(const ParserSection & parameters_section); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ -protected: - typedef std::map NonLinearSolverCallbackMap; +private: + DOFManager & _dof_manager; +protected: /// Set of callbacks to use in the solver for jacobian assembly, residual /// assembly, corrector & predictor if needed - NonLinearSolverCallbackMap solver_callbacks; + NonLinearSolverCallback * solver_callback; /// type of non linear solver NonLinearSolverType non_linear_solver_type; }; __END_AKANTU__ #endif /* __AKANTU_NON_LINEAR_SOLVER_HH__ */ diff --git a/src/model/non_linear_solver_default.cc b/src/model/non_linear_solver_default.cc index f7d578911..6ca98f431 100644 --- a/src/model/non_linear_solver_default.cc +++ b/src/model/non_linear_solver_default.cc @@ -1,168 +1,168 @@ /** * @file non_linear_solver_default.cc * * @author Nicolas Richart * * @date Tue Aug 25 00:57:00 2015 * * @brief Implementation of the default NonLinearSolver * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "non_linear_solver_default.hh" #include "dof_manager_default.hh" #include "non_linear_solver_callback.hh" #include "static_communicator.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ NonLinearSolverDefault::NonLinearSolverDefault( DOFManagerDefault & dof_manager, const NonLinearSolverType & non_linear_solver_type, const ID & id, UInt memory_id) : NonLinearSolver(dof_manager, non_linear_solver_type, id, memory_id), dof_manager(dof_manager), solver(dof_manager, "jacobian", id + ":sparse_solver", memory_id), convergence_criteria_type(_scc_solution), convergence_criteria(1e-10), max_iterations(10), n_iter(0), error(0.), converged(false) {} /* -------------------------------------------------------------------------- */ NonLinearSolverDefault::~NonLinearSolverDefault() {} /* ------------------------------------------------------------------------ */ void NonLinearSolverDefault::solve() { // EventManager::sendEvent(NonLinearSolver::BeforeNonLinearSolverSolve(method)); - this->callbackPredictors(); + this->solver_callback->predictor(); switch (this->non_linear_solver_type) { case _nls_linear: case _nls_newton_raphson: break; case _nls_newton_raphson_modified: - this->callbackAssembleJacobian(); + this->solver_callback->assembleJacobian(); break; default: AKANTU_DEBUG_ERROR("The resolution method " << this->non_linear_solver_type << " has not been implemented!"); } this->n_iter = 0; this->converged = false; if (this->convergence_criteria_type == _scc_residual) { this->converged = this->testConvergence(this->dof_manager.getResidual()); if (this->converged) return; } do { if (this->non_linear_solver_type == _nls_newton_raphson) - this->callbackAssembleJacobian(); + this->solver_callback->assembleJacobian(); this->solver.solve(); - this->callbackCorrectors(); + this->solver_callback->corrector(); // EventManager::sendEvent(NonLinearSolver::AfterSparseSolve(method)); if (this->convergence_criteria_type == _scc_residual) { - this->callbackAssembleResidual(); + this->solver_callback->assembleResidual(); this->converged = this->testConvergence(this->dof_manager.getResidual()); } else { this->converged = this->testConvergence(this->dof_manager.getGlobalSolution()); } if (this->convergence_criteria_type == _scc_solution && !this->converged) - this->callbackAssembleResidual(); + this->solver_callback->assembleResidual(); // this->dump(); this->n_iter++; AKANTU_DEBUG_INFO( "[" << this->convergence_criteria_type << "] Convergence iteration " << std::setw(std::log10(this->max_iterations)) << this->n_iter << ": error " << this->error << (this->converged ? " < " : " > ") << this->convergence_criteria); } while (!this->converged && this->n_iter < this->max_iterations); // this makes sure that you have correct strains and stresses after the // solveStep function (e.g., for dumping) if (this->convergence_criteria_type == _scc_solution) - this->callbackAssembleResidual(); + this->solver_callback->assembleResidual(); if (this->converged) { // EventManager::sendEvent( // SolidMechanicsModelEvent::AfterNonLinearSolverSolves(method)); } else if (this->n_iter == this->max_iterations) { AKANTU_DEBUG_WARNING("[" << this->convergence_criteria_type << "] Convergence not reached after " << std::setw(std::log10(this->max_iterations)) << this->n_iter << " iteration" << (this->n_iter == 1 ? "" : "s") << "!" << std::endl); } return; } /* -------------------------------------------------------------------------- */ void NonLinearSolverDefault::setParameters( const ParserSection & parameters_section) {} /* -------------------------------------------------------------------------- */ bool NonLinearSolverDefault::testConvergence(const Array & array) { AKANTU_DEBUG_IN(); const Array & blocked_dofs = this->dof_manager.getGlobalBlockedDOFs(); UInt nb_degree_of_freedoms = array.getSize(); Array::const_scalar_iterator arr_it = array.begin(); Array::const_scalar_iterator bld_it = blocked_dofs.begin(); Real norm = 0.; for (UInt n = 0; n < nb_degree_of_freedoms; ++n, ++arr_it, ++bld_it) { bool is_local_node = this->dof_manager.isLocalOrMasterDOF(n); if(!(*bld_it) && is_local_node) { norm += *arr_it * *arr_it; } } StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum); norm = std::sqrt(norm); AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase"); this->error = norm; return (error < this->convergence_criteria); } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/time_step_solver.hh b/src/model/time_step_solver.hh index ccc0dc157..3e086c26c 100644 --- a/src/model/time_step_solver.hh +++ b/src/model/time_step_solver.hh @@ -1,92 +1,92 @@ /** * @file time_step_solver.hh * * @author Nicolas Richart * * @date Mon Aug 24 12:42:04 2015 * * @brief This corresponding to the time step evolution solver * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_memory.hh" #include "aka_array.hh" +#include "non_linear_solver_callback.hh" + /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_TIME_STEP_SOLVER_HH__ #define __AKANTU_TIME_STEP_SOLVER_HH__ namespace akantu { class DOFManager; } __BEGIN_AKANTU__ -class TimeStepSolver : public Memory { +class TimeStepSolver : public Memory, public NonLinearSolverCallback { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - TimeStepSolver(DOFManager & dof_manager, const ID & dof_id, - const TimeStepSolverType & type, const ID & id, - UInt memory_id); + TimeStepSolver(DOFManager & dof_manager, const TimeStepSolverType & type, + const ID & id, UInt memory_id); virtual ~TimeStepSolver(); /* ------------------------------------------------------------------------ */ - /* Methods */ + /* Accessors */ /* ------------------------------------------------------------------------ */ public: - void solveStep(); -private: - virtual void predictor() = 0; + /// register a callback object for a given dof_id + virtual void registerCallback(NonLinearSolverCallback & callbacks); - virtual void corrector() = 0; + /// register a callback object for a given dof_id + virtual NonLinearSolverCallback & getCallbacks(); - /* ------------------------------------------------------------------------ */ - /* Accessors */ - /* ------------------------------------------------------------------------ */ -public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// Underlying dof manager containing the dof to treat DOFManager & _dof_manager; /// The dof to make evolve ID dof_id; /// Type of solver TimeStepSolverType type; /// The increment of the dof from @f[ u_{n+1} = u_{n} + inc @f] Array increment; /// The time step for this solver Real time_step; + + /// NonLinearSolverCallback to use to get the dofs updates + NonLinearSolverCallback * callbacks; }; __END_AKANTU__ #endif /* __AKANTU_TIME_STEP_SOLVER_HH__ */ diff --git a/src/model/time_step_solvers/time_step_solver.cc b/src/model/time_step_solvers/time_step_solver.cc index 8616755a7..496911e26 100644 --- a/src/model/time_step_solvers/time_step_solver.cc +++ b/src/model/time_step_solvers/time_step_solver.cc @@ -1,42 +1,50 @@ /** * @file time_step_solver.cc * * @author Nicolas Richart * * @date Mon Oct 12 16:56:43 2015 * * @brief Implementation of common part of TimeStepSolvers * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "time_step_solver.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ -TimeStepSolver::TimeStepSolver(DOFManager & dof_manager, const ID & dof_id, +/* -------------------------------------------------------------------------- */ +TimeStepSolver::TimeStepSolver(DOFManager & dof_manager, const TimeStepSolverType & type, const ID & id, UInt memory_id) - : Memory(id, memory_id), _dof_manager(dof_manager), dof_id(dof_id), - type(type) {} + : Memory(id, memory_id), _dof_manager(dof_manager), type(type) {} + +/* -------------------------------------------------------------------------- */ +TimeStepSolver::~TimeStepSolver() {} + +/* -------------------------------------------------------------------------- */ + + +/* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/time_step_solvers/time_step_solver_default.cc b/src/model/time_step_solvers/time_step_solver_default.cc index 7620d8e33..1451bebb8 100644 --- a/src/model/time_step_solvers/time_step_solver_default.cc +++ b/src/model/time_step_solvers/time_step_solver_default.cc @@ -1,218 +1,214 @@ /** * @file time_step_solver_default.cc * * @author Nicolas Richart * * @date Wed Sep 16 10:20:55 2015 * * @brief Default implementation of the time step solver * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "time_step_solver_default.hh" #include "dof_manager_default.hh" +#include "sparse_matrix_aij.hh" + #include "integration_scheme_1st_order.hh" #include "integration_scheme_2nd_order.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ // void TimeStepSolverDefault::updateAcceleration() { // AKANTU_DEBUG_IN(); // updateResidualInternal(); // if (method == _explicit_lumped_mass) { // /* residual = residual_{n+1} - M * acceleration_n therefore // solution = increment acceleration not acceleration */ // solveLumped(*increment_acceleration, *mass, *residual, *blocked_dofs, // f_m2a); // } else if (method == _explicit_consistent_mass) { // solve(*increment_acceleration); // } // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ TimeStepSolverDefault::TimeStepSolverDefault(DOFManagerDefault & dof_manager, - const ID & dof_id, const TimeStepSolverType & type, const ID & id, UInt memory_id) - : TimeStepSolver(dof_manager, dof_id, type, id, memory_id), - dof_manager(dof_manager) { + : TimeStepSolver(dof_manager, type, id, memory_id), + dof_manager(dof_manager), is_mass_lumped(false) { + switch (type) { + case _tsst_forward_euler_lumped: + this->is_mass_lumped = true; case _tsst_forward_euler: { - this->integration_scheme = new ForwardEuler(); + + this->integration_scheme = new ForwardEuler(dof_manager); break; } case _tsst_trapezoidal_rule_1: { - this->integration_scheme = new TrapezoidalRule1(); + this->integration_scheme = new TrapezoidalRule1(dof_manager); break; } case _tsst_backward_euler: { - this->integration_scheme = new BackwardEuler(); + this->integration_scheme = new BackwardEuler(dof_manager); break; } + case _tsst_central_difference_lumped: + this->is_mass_lumped = true; case _tsst_central_difference: { - this->integration_scheme = new CentralDifference(); + this->integration_scheme = new CentralDifference(dof_manager); break; } case _tsst_fox_goodwin: { - this->integration_scheme = new FoxGoodwin(); + this->integration_scheme = new FoxGoodwin(dof_manager); break; } case _tsst_trapezoidal_rule_2: { - this->integration_scheme = new TrapezoidalRule2(); + this->integration_scheme = new TrapezoidalRule2(dof_manager); break; } case _tsst_linear_acceleration: { - this->integration_scheme = new LinearAceleration(); + this->integration_scheme = new LinearAceleration(dof_manager); break; } // Write a c++11 version of the constructor with initializer list that // contains the arguments for the integration scheme case _tsst_generalized_trapezoidal: case _tsst_newmark_beta: AKANTU_EXCEPTION( "This time step solvers cannot be created with this constructor"); } } /* -------------------------------------------------------------------------- */ TimeStepSolverDefault::~TimeStepSolverDefault() { delete this->integration_scheme; } /* -------------------------------------------------------------------------- */ void TimeStepSolverDefault::predictor() { AKANTU_DEBUG_IN(); - Array & u = this->dof_manager.getDOFs(this->dof_id); - const Array & blocked_dofs = - this->dof_manager.getBlockedDOFs(this->dof_id); - - // increment.copy(u); - - if (this->integration_scheme->getOrder() == 1) { - Array & u_dot = dof_manager.getDOFsDerivatives(this->dof_id, 1); - IntegrationScheme1stOrder & int_scheme = - *dynamic_cast(this->integration_scheme); - int_scheme.integrationSchemePred(this->time_step, u, u_dot, blocked_dofs); - } else if (this->integration_scheme->getOrder() == 2) { + TimeStepSolver::predictor(); - Array & u_dot = dof_manager.getDOFsDerivatives(this->dof_id, 1); - Array & u_dot_dot = dof_manager.getDOFsDerivatives(this->dof_id, 2); - - IntegrationScheme2ndOrder & int_scheme = - *dynamic_cast(this->integration_scheme); - int_scheme.integrationSchemePred(this->time_step, u, u_dot, u_dot_dot, - blocked_dofs); - } + this->integration_scheme->predictor(this->dof_id, this->time_step); // UInt nb_degree_of_freedom = u.getSize() * u.getNbComponent(); // Array::scalar_iterator incr_it = // increment.begin_reinterpret(nb_degree_of_freedom); // Array::const_scalar_iterator u_it = // u.begin_reinterpret(nb_degree_of_freedom); // Array::const_scalar_iterator u_end = // u.end_reinterpret(nb_degree_of_freedom); // for (; u_it != u_end; ++u_it, ++incr_it) { // *incr_it = *u_it - *incr_it; // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void TimeStepSolverDefault::corrector() { AKANTU_DEBUG_IN(); - Array & u = this->dof_manager.getDOFs(this->dof_id); - const Array & solution = this->dof_manager.getSolution(this->dof_id); - const Array & blocked_dofs = - this->dof_manager.getBlockedDOFs(this->dof_id); - - // increment.copy(u); - - if (this->integration_scheme->getOrder() == 1) { - Array & u_dot = dof_manager.getDOFsDerivatives(this->dof_id, 1); - IntegrationScheme1stOrder & int_scheme = - *dynamic_cast(this->integration_scheme); - - switch (this->corrector_type) { - case IntegrationScheme1stOrder::_temperature_corrector: - int_scheme.integrationSchemeCorrTemp(this->time_step, u, u_dot, - blocked_dofs, solution); - break; - case IntegrationScheme1stOrder::_temperature_rate_corrector: - int_scheme.integrationSchemeCorrTempRate(this->time_step, u, u_dot, - blocked_dofs, solution); - break; - } - } else if (this->integration_scheme->getOrder() == 2) { - - Array & u_dot = dof_manager.getDOFsDerivatives(this->dof_id, 1); - Array & u_dot_dot = dof_manager.getDOFsDerivatives(this->dof_id, 2); - - IntegrationScheme2ndOrder & int_scheme = - *dynamic_cast(this->integration_scheme); - - switch (this->corrector_type) { - case IntegrationScheme2ndOrder::_displacement_corrector: - int_scheme.integrationSchemeCorrDispl(this->time_step, u, u_dot, - u_dot_dot, blocked_dofs, solution); - break; - case IntegrationScheme2ndOrder::_velocity_corrector: - int_scheme.integrationSchemeCorrVeloc(this->time_step, u, u_dot, - u_dot_dot, blocked_dofs, solution); - break; - case IntegrationScheme2ndOrder::_acceleration_corrector: - int_scheme.integrationSchemeCorrAccel(this->time_step, u, u_dot, - u_dot_dot, blocked_dofs, solution); - break; - } - } + TimeStepSolver::corrector(); + + this->integration_scheme->corrector( + IntegrationScheme::SolutionType(this->solution_type), this->dof_id, + this->time_step); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void TimeStepSolverDefault::solveStep() { +void TimeStepSolverDefault::assembleJacobian() { AKANTU_DEBUG_IN(); - // EventManager::sendEvent( - // SolidMechanicsModelEvent::BeforeSolveStepEvent(method)); + TimeStepSolver::assembleJacobian(); - // this->predictor(); - // this->solver->solve(); - // this->corrector(); + this->integration_scheme->assembleJacobian( + IntegrationScheme::SolutionType(this->solution_type), this->time_step); - // EventManager::sendEvent( - // SolidMechanicsModelEvent::AfterSolveStepEvent(method)); + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void TimeStepSolverDefault::assembleResidual() { + AKANTU_DEBUG_IN(); + + TimeStepSolver::assembleResidual(); + + // if (!this->is_mass_lumped) { + + // Array * Ma = new Array(*acceleration, true, "Ma"); + // *Ma *= *mass_matrix; + // /// \todo check unit conversion for implicit dynamics + // // *Ma /= f_m2a + // *residual -= *Ma; + // delete Ma; + // } else if (mass) { + + // // else lumped mass + // UInt nb_nodes = acceleration->getSize(); + // UInt nb_degree_of_freedom = acceleration->getNbComponent(); + + // Real * mass_val = mass->storage(); + // Real * accel_val = acceleration->storage(); + // Real * res_val = residual->storage(); + // bool * blocked_dofs_val = blocked_dofs->storage(); + + // for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) { + // if (!(*blocked_dofs_val)) { + // *res_val -= *accel_val * *mass_val / f_m2a; + // } + // blocked_dofs_val++; + // res_val++; + // mass_val++; + // accel_val++; + // } + // } else { + // AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix."); + // } + + // // f -= Cv + // if (velocity_damping_matrix) { + // Array * Cv = new Array(*velocity); + // *Cv *= *velocity_damping_matrix; + // *residual -= *Cv; + // delete Cv; + // } + // } AKANTU_DEBUG_OUT(); } +/* -------------------------------------------------------------------------- */ + __END_AKANTU__ diff --git a/src/model/time_step_solvers/time_step_solver_default.hh b/src/model/time_step_solvers/time_step_solver_default.hh index 8b310740d..3b05cdb8c 100644 --- a/src/model/time_step_solvers/time_step_solver_default.hh +++ b/src/model/time_step_solvers/time_step_solver_default.hh @@ -1,87 +1,94 @@ /** * @file time_step_solver_default.hh * * @author Nicolas Richart * * @date Mon Aug 24 17:10:29 2015 * * @brief Default implementation for the time stepper * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "time_step_solver.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_TIME_STEP_SOLVER_DEFAULT_HH__ #define __AKANTU_TIME_STEP_SOLVER_DEFAULT_HH__ namespace akantu { -class IntegrationScheme; -class DOFManagerDefault; + class IntegrationScheme; + class DOFManagerDefault; } __BEGIN_AKANTU__ class TimeStepSolverDefault : public TimeStepSolver { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - TimeStepSolverDefault(DOFManagerDefault & dof_manager, const ID & dof_id, + TimeStepSolverDefault(DOFManagerDefault & dof_manager, const TimeStepSolverType & type, const ID & id, UInt memory_id); virtual ~TimeStepSolverDefault(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// implementation of the TimeStepSolver::predictor() virtual void predictor(); /// implementation of the TimeStepSolver::corrector() virtual void corrector(); + /// implementation of the TimeStepSolver::assembleJacobian() + virtual void assembleJacobian(); + /// implementation of the TimeStepSolver::assembleResidual() + virtual void assembleResidual(); /// implementation of the generic TimeStepSolver::solveStep() virtual void solveStep(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// DOFManager with its real type DOFManagerDefault & dof_manager; /// Underlying integration scheme IntegrationScheme * integration_scheme; /// Type of corrector to use - UInt corrector_type; + UInt solution_type; + + /// define if the mass matrix is lumped or not + bool is_mass_lumped; }; __END_AKANTU__ #endif /* __AKANTU_TIME_STEP_SOLVER_DEFAULT_HH__ */