diff --git a/packages/core.cmake b/packages/core.cmake
index 20c0fda86..c8a07a8e0 100644
--- a/packages/core.cmake
+++ b/packages/core.cmake
@@ -1,514 +1,516 @@
#===============================================================================
# @file core.cmake
#
# @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
# @author Nicolas Richart <nicolas.richart@epfl.ch>
#
# @date creation: Mon Nov 21 2011
# @date last modification: Mon Jan 18 2016
#
# @brief package description for core
#
# @section LICENSE
#
# Copyright (©) 2010-2012, 2014, 2015 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 <http://www.gnu.org/licenses/>.
#
#===============================================================================
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/integration_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
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_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_real.cc
io/parser/parser_random.cc
io/parser/parser_types.cc
io/parser/parser_input_files.cc
io/parser/parser.hh
io/parser/parser_tmpl.hh
io/parser/parser_grammar_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_utils_pbc.cc
mesh_utils/mesh_utils.cc
mesh_utils/mesh_utils.hh
mesh_utils/mesh_utils_inline_impl.cc
mesh_utils/global_ids_updater.hh
mesh_utils/global_ids_updater.cc
mesh_utils/global_ids_updater_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.cc
model/non_linear_solver.hh
model/non_linear_solver_default.hh
model/non_linear_solver_lumped.cc
model/non_linear_solver_lumped.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/time_step_solvers/time_step_solver_default_solver_callback.hh
- model/integration_scheme/generalized_trapezoidal.hh
model/integration_scheme/generalized_trapezoidal.cc
- model/integration_scheme/integration_scheme.hh
+ model/integration_scheme/generalized_trapezoidal.hh
model/integration_scheme/integration_scheme.cc
- model/integration_scheme/integration_scheme_1st_order.hh
+ model/integration_scheme/integration_scheme.hh
model/integration_scheme/integration_scheme_1st_order.cc
- model/integration_scheme/integration_scheme_2nd_order.hh
+ model/integration_scheme/integration_scheme_1st_order.hh
model/integration_scheme/integration_scheme_2nd_order.cc
- model/integration_scheme/newmark-beta.hh
+ model/integration_scheme/integration_scheme_2nd_order.hh
model/integration_scheme/newmark-beta.cc
+ model/integration_scheme/newmark-beta.hh
+ model/integration_scheme/pseudo_time.cc
+ model/integration_scheme/pseudo_time.hh
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
model/common/neighborhood_base.hh
model/common/neighborhood_base.cc
model/common/neighborhood_base_inline_impl.cc
model/common/neighborhoods_criterion_evaluation/neighborhood_max_criterion.hh
model/common/neighborhoods_criterion_evaluation/neighborhood_max_criterion.cc
model/common/neighborhoods_criterion_evaluation/neighborhood_max_criterion_inline_impl.cc
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/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
synchronizer/grid_synchronizer.cc
synchronizer/grid_synchronizer.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_integration_points
interpolate_on_integration_points
interpolate
compute_normals_on_integration_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-feengine.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)
find_program(PATCH_COMMAND patch)
mark_as_advanced(READLINK_COMMAND)
mark_as_advanced(ADDR2LINE_COMMAND)
include(CheckFunctionExists)
check_function_exists(clock_gettime _clock_gettime)
include(CheckCXXSymbolExists)
check_cxx_symbol_exists(strdup cstring AKANTU_HAS_STRDUP)
if(NOT _clock_gettime)
set(AKANTU_USE_OBSOLETE_GETTIMEOFDAY ON CACHE INTERNAL "" FORCE)
else()
set(AKANTU_USE_OBSOLETE_GETTIMEOFDAY OFF CACHE INTERNAL "" FORCE)
endif()
package_declare_extra_files_to_package(core
SOURCES
common/aka_element_classes_info.hh.in
common/aka_config.hh.in
model/solid_mechanics/material_list.hh.in
)
diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh
index dcb3d0b9e..631a6d501 100644
--- a/src/common/aka_common.hh
+++ b/src/common/aka_common.hh
@@ -1,429 +1,430 @@
/**
* @file aka_common.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Mon Jun 14 2010
* @date last modification: Thu Jan 21 2016
*
* @brief common type descriptions for akantu
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 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 <http://www.gnu.org/licenses/>.
*
* @section DESCRIPTION
*
* All common things to be included in the projects files
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_COMMON_HH__
#define __AKANTU_COMMON_HH__
/* -------------------------------------------------------------------------- */
#include <list>
#include <limits>
/* -------------------------------------------------------------------------- */
#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;
#ifdef AKANTU_NDEBUG
static const Real REAL_INIT_VALUE = Real(0.);
#else
static const Real REAL_INIT_VALUE = std::numeric_limits<Real>::quiet_NaN();
#endif
/* -------------------------------------------------------------------------- */
/* Memory types */
/* -------------------------------------------------------------------------- */
typedef UInt MemoryID;
typedef std::string Surface;
typedef std::pair<Surface, Surface> 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 = 0,
_implicit_dynamic = 1,
_explicit_lumped_mass = 2,
_explicit_lumped_capacity = 2,
_explicit_consistent_mass = 3
};
/// 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
_nls_lumped, ///< Case of lumped mass or equivalent matrix
_nls_auto, ///< This will take a default value that make sense in case of
///model::getNewSolver
};
/// Type of time stepping solver
enum TimeStepSolverType {
_tsst_static, ///< Static solution
_tsst_dynamic, ///< Dynamic solver
_tsst_dynamic_lumped, ///< Dynamic solver with lumped mass
};
/// Type of integration scheme
enum IntegrationSchemeType {
+ _ist_pseudo_time, ///< Pseudo Time
_ist_forward_euler, ///< GeneralizedTrapezoidal(0)
_ist_trapezoidal_rule_1, ///< GeneralizedTrapezoidal(1/2)
_ist_backward_euler, ///< GeneralizedTrapezoidal(1)
_ist_central_difference, ///< NewmarkBeta(0, 1/2)
_ist_fox_goodwin, ///< NewmarkBeta(1/6, 1/2)
_ist_trapezoidal_rule_2, ///< NewmarkBeta(1/2, 1/2)
_ist_linear_acceleration, ///< NewmarkBeta(1/3, 1/2)
_ist_newmark_beta, ///< generic NewmarkBeta with user defined
/// alpha and beta
_ist_generalized_trapezoidal, ///< generic GeneralizedTrapezoidal with user
/// defined alpha
};
/// 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 SparseMatrixType type of sparse matrix used
enum MatrixType { _unsymmetric, _symmetric };
/* -------------------------------------------------------------------------- */
/* 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
// --- GlobalIdsUpdater tags ---
_gst_giu_global_conn, ///< synchronization of global connectivities
// --- CohesiveElementInserter tags ---
_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 ---
_gst_htm_phi, ///< synchronization of the nodal level set value phi
_gst_htm_gradient_phi, ///< synchronization of the element 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
// --- NeighborhoodSynchronization tags ---
_gst_nh_criterion,
// --- 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<GhostType_def> 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_prod,
_so_land,
_so_band,
_so_lor,
_so_bor,
_so_lxor,
_so_bxor,
_so_min_loc,
_so_max_loc,
_so_null
};
/* -------------------------------------------------------------------------- */
/* Global defines */
/* -------------------------------------------------------------------------- */
#define AKANTU_MIN_ALLOCATION 2000
#define AKANTU_INDENT " "
#define AKANTU_INCLUDE_INLINE_IMPL
/* -------------------------------------------------------------------------- */
template <class T> struct is_scalar {
enum { value = false };
};
#define AKANTU_SPECIFY_IS_SCALAR(type) \
template <> struct is_scalar<type> { \
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<T, T> {
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<type> & 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();
/* -------------------------------------------------------------------------- */
/// Read an new input file
void readInputFile(const std::string & input_file);
/* -------------------------------------------------------------------------- */
/*
* 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 <typename T> 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"
/* -------------------------------------------------------------------------- */
#if defined(AKANTU_UNORDERED_MAP_IS_CXX11)
__BEGIN_AKANTU_UNORDERED_MAP__
#if AKANTU_INTEGER_SIZE == 4
#define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b9
#elif AKANTU_INTEGER_SIZE == 8
#define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b97f4a7c13LL
#endif
/**
* Hashing function for pairs based on hash_combine from boost The magic number
* is coming from the golden number @f[\phi = \frac{1 + \sqrt5}{2}@f]
* @f[\frac{2^32}{\phi} = 0x9e3779b9@f]
* http://stackoverflow.com/questions/4948780/magic-number-in-boosthash-combine
* http://burtleburtle.net/bob/hash/doobs.html
*/
template <typename a, typename b> struct hash<std::pair<a, b> > {
public:
hash() : ah(), bh() {}
size_t operator()(const std::pair<a, b> & p) const {
size_t seed = ah(p.first);
return bh(p.second) + AKANTU_HASH_COMBINE_MAGIC_NUMBER + (seed << 6) +
(seed >> 2);
}
private:
const hash<a> ah;
const hash<b> bh;
};
__END_AKANTU_UNORDERED_MAP__
#endif
#endif /* __AKANTU_COMMON_HH__ */
diff --git a/src/model/dof_manager.hh b/src/model/dof_manager.hh
index 68c17fbbd..ba899b11d 100644
--- a/src/model/dof_manager.hh
+++ b/src/model/dof_manager.hh
@@ -1,306 +1,307 @@
/**
* @file dof_manager.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_memory.hh"
#include "non_linear_solver.hh"
#include "time_step_solver.hh"
+#include "mesh_events.hh"
/* -------------------------------------------------------------------------- */
#include <map>
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_DOF_MANAGER_HH__
#define __AKANTU_DOF_MANAGER_HH__
__BEGIN_AKANTU__
-class DOFManager : protected Memory {
+class DOFManager : protected Memory, protected MeshEventHandler {
/* ------------------------------------------------------------------------ */
/* 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<Real> & 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<Real> & dofs_derivative);
/// register array representing the blocked degree of freedoms
void registerBlockedDOFs(const ID & dof_id, Array<bool> & blocked_dofs);
/// Assemble an array to the global residual array
virtual void assembleToResidual(const ID & dof_id,
const Array<Real> & 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<Real> & elementary_vect, Array<Real> & array_assembeled,
const ElementType & type, const GhostType & ghost_type,
Real scale_factor = 1.,
const Array<UInt> & 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<Real> & elementary_vect,
const ElementType & type, const GhostType & ghost_type,
Real scale_factor = 1.,
const Array<UInt> & 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<Real> & elementary_mat, const ElementType & type,
const GhostType & ghost_type,
const MatrixType & elemental_matrix_type = _symmetric,
const Array<UInt> & 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<Real> & matrix) = 0;
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<Real> & solution_array) = 0;
/// fill a Vector with the equation numbers corresponding to the given
/// connectivity
inline void extractElementEquationNumber(
const Array<UInt> & equation_numbers, const Vector<UInt> & connectivity,
UInt nb_degree_of_freedom, Vector<UInt> & local_equation_number);
/// converts local equation numbers to global equation numbers;
template <class S> inline void localToGlobalEquationNumber(S & inout);
/* ------------------------------------------------------------------------ */
/// register a matrix
void registerSparseMatrix(const ID & matrix_id, SparseMatrix & matrix);
/// register a non linear solver instantiated by a derived class
void registerNonLinearSolver(const ID & non_linear_solver_id,
NonLinearSolver & non_linear_solver);
/// register a time step solver instantiated by a derived class
void registerTimeStepSolver(const ID & time_step_solver_id,
TimeStepSolver & time_step_solver);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
protected:
struct DOFData;
inline DOFData & getDOFData(const ID & dof_id);
inline const DOFData & getDOFData(const ID & dof_id) const;
public:
/// get the equation numbers (in local numbering) corresponding to a dof ID
inline const Array<UInt> & getLocalEquationNumbers(const ID & dof_id) const;
/// return the local index of the global equation number
inline UInt globalToLocalEquationNumber(UInt global) 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);
/* ------------------------------------------------------------------------ */
/* DOFs and derivatives accessors */
/* ------------------------------------------------------------------------ */
/// Get a reference to the registered dof array for a given id
inline Array<Real> & getDOFs(const ID & dofs_id);
/// Get a reference to the registered dof derivatives array for a given id
inline Array<Real> & getDOFsDerivatives(const ID & dofs_id, UInt order);
/// Get a reference to the blocked dofs array registered for the given id
inline const Array<bool> & getBlockedDOFs(const ID & dofs_id) const;
/// Get a reference to the solution array registered for the given id
inline const Array<Real> & getSolution(const ID & dofs_id) const;
/* ------------------------------------------------------------------------ */
/* Matrices accessors */
/* ------------------------------------------------------------------------ */
/// Get an instance of a new SparseMatrix
virtual SparseMatrix & getNewMatrix(const ID & matrix_id,
const MatrixType & matrix_type) = 0;
/// 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) = 0;
/// Get the reference of an existing matrix
SparseMatrix & getMatrix(const ID & matrix_id);
/// Get an instance of a new lumped matrix
virtual Array<Real> & getNewLumpedMatrix(const ID & matrix_id);
/// Get the lumped version of a given matrix
const Array<Real> & getLumpedMatrix(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) = 0;
/// get instance of a non linear solver
virtual NonLinearSolver & getNonLinearSolver(const ID & nls_solver_id);
/* ------------------------------------------------------------------------ */
/* Time-Step Solver */
/* ------------------------------------------------------------------------ */
/// Get instance of a time step solver
virtual TimeStepSolver &
getNewTimeStepSolver(const ID & time_step_solver_id,
const TimeStepSolverType & type,
NonLinearSolver & non_linear_solver) = 0;
/// 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<Real> * dof;
/// Blocked degree of freedoms array
Array<bool> * blocked_dofs;
/// Solution associated to the dof
Array<Real> * solution;
/* ---------------------------------------------------------------------- */
/* data for dynamic simulations */
/* ---------------------------------------------------------------------- */
/// Degree of freedom derivatives arrays
std::vector<Array<Real> *> dof_derivatives;
/// local numbering equation numbers
Array<UInt> local_equation_number;
};
/// equation number in global numbering
Array<UInt> global_equation_number;
typedef unordered_map<UInt, UInt>::type equation_numbers_map;
/// dual information of global_equation_number
equation_numbers_map global_to_local_mapping;
/// type to store dofs informations
typedef std::map<ID, DOFData *> DOFStorage;
/// type to store all the matrices
typedef std::map<ID, SparseMatrix *> SparseMatricesMap;
/// type to store all the lumped matrices
typedef std::map<ID, Array<Real> *> LumpedMatricesMap;
/// type to store all the non linear solver
typedef std::map<ID, NonLinearSolver *> NonLinearSolversMap;
/// type to store all the time step solver
typedef std::map<ID, TimeStepSolver *> TimeStepSolversMap;
/// store a reference to the dof arrays
DOFStorage dofs;
/// list of sparse matrices that where created
SparseMatricesMap matrices;
/// list of lumped matrices
LumpedMatricesMap lumped_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;
/// 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/integration_scheme/pseudo_time.cc b/src/model/integration_scheme/pseudo_time.cc
new file mode 100644
index 000000000..017cd34c3
--- /dev/null
+++ b/src/model/integration_scheme/pseudo_time.cc
@@ -0,0 +1,85 @@
+/**
+ * @file pseudo_time.cc
+ *
+ * @author Nicolas Richart <nicolas.richart@epfl.ch>
+ *
+ * @date Wed Feb 17 09:49:10 2016
+ *
+ * @brief Implementation of a really simple integration scheme
+ *
+ * @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 <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+#include "pseudo_time.hh"
+#include "dof_manager.hh"
+#include "sparse_matrix.hh"
+/* -------------------------------------------------------------------------- */
+
+__BEGIN_AKANTU__
+
+/* -------------------------------------------------------------------------- */
+PseudoTime::PseudoTime(DOFManager & dof_manager, const ID & dof_id)
+ : IntegrationScheme(dof_manager, dof_id, 0) {}
+
+/* -------------------------------------------------------------------------- */
+void PseudoTime::predictor(Real delta_t) {}
+
+/* -------------------------------------------------------------------------- */
+void PseudoTime::corrector(const SolutionType & type, Real delta_t) {
+ Array<Real> & u = this->dof_manager.getDOFs(this->dof_id);
+ const Array<Real> & delta = this->dof_manager.getSolution(this->dof_id);
+ const Array<bool> & blocked_dofs =
+ this->dof_manager.getBlockedDOFs(this->dof_id);
+
+ UInt nb_nodes = u.getSize();
+ UInt nb_degree_of_freedom = u.getNbComponent() * nb_nodes;
+
+ Array<Real>::scalar_iterator u_it = u.begin_reinterpret(nb_degree_of_freedom);
+ Array<Real>::scalar_iterator u_end = u.end_reinterpret(nb_degree_of_freedom);
+ Array<Real>::const_scalar_iterator delta_it =
+ delta.begin_reinterpret(nb_degree_of_freedom);
+ Array<bool>::const_scalar_iterator blocked_dofs_it =
+ blocked_dofs.begin_reinterpret(nb_degree_of_freedom);
+
+ for (; u_it != u_end; ++u_it) {
+ if (!(*blocked_dofs_it)) {
+ *u_it += *delta_it;
+ }
+ ++u_it;
+ ++delta_it;
+ ++blocked_dofs_it;
+ }
+}
+
+/* -------------------------------------------------------------------------- */
+void PseudoTime::assembleJacobian(const SolutionType & type, Real delta_t) {
+ SparseMatrix & J = this->dof_manager.getMatrix("J");
+ const SparseMatrix & K = this->dof_manager.getMatrix("K");
+
+ J.add(K);
+}
+
+/* -------------------------------------------------------------------------- */
+void PseudoTime::assembleResidual(bool is_lumped) {}
+
+/* -------------------------------------------------------------------------- */
+
+__END_AKANTU__
diff --git a/src/model/integration_scheme/pseudo_time.hh b/src/model/integration_scheme/pseudo_time.hh
new file mode 100644
index 000000000..c409ab342
--- /dev/null
+++ b/src/model/integration_scheme/pseudo_time.hh
@@ -0,0 +1,68 @@
+/**
+ * @file pseudo_time.hh
+ *
+ * @author Nicolas Richart <nicolas.richart@epfl.ch>
+ *
+ * @date Wed Feb 17 09:46:05 2016
+ *
+ * @brief Pseudo time integration scheme
+ *
+ * @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 <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+#include "integration_scheme.hh"
+/* -------------------------------------------------------------------------- */
+
+
+#ifndef __AKANTU_PSEUDO_TIME_HH__
+#define __AKANTU_PSEUDO_TIME_HH__
+
+__BEGIN_AKANTU__
+
+class PseudoTime : public IntegrationScheme {
+ /* ------------------------------------------------------------------------ */
+ /* Constructors/Destructors */
+ /* ------------------------------------------------------------------------ */
+public:
+ PseudoTime(DOFManager & dof_manager, const ID & dof_id);
+ virtual ~PseudoTime() {}
+
+ /* ------------------------------------------------------------------------ */
+ /* Methods */
+ /* ------------------------------------------------------------------------ */
+public:
+ /// generic interface of a predictor
+ virtual void predictor(Real delta_t);
+
+ /// generic interface of a corrector
+ virtual void corrector(const SolutionType & type, Real delta_t);
+
+ /// assemble the jacobian matrix
+ virtual void assembleJacobian(const SolutionType & type,
+ Real delta_t);
+
+ /// assemble the residual
+ virtual void assembleResidual(bool is_lumped);
+};
+
+__END_AKANTU__
+
+#endif /* __AKANTU_PSEUDO_TIME_HH__ */
diff --git a/src/model/model_solver.cc b/src/model/model_solver.cc
index e5474d278..76586bd97 100644
--- a/src/model/model_solver.cc
+++ b/src/model/model_solver.cc
@@ -1,457 +1,469 @@
/**
* @file model_solver.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "model_solver.hh"
#include "mesh.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), parent_id(id), parent_memory_id(memory_id), mesh(mesh),
- dof_manager(NULL) {}
+ : Parsable(_st_solver, id), parent_id(id), parent_memory_id(memory_id),
+ mesh(mesh), dof_manager(NULL) {}
/* -------------------------------------------------------------------------- */
-ModelSolver::~ModelSolver() {
- delete this->dof_manager;
-}
+ModelSolver::~ModelSolver() { delete this->dof_manager; }
/* -------------------------------------------------------------------------- */
void ModelSolver::initDOFManager() {
std::pair<Parser::const_section_iterator, Parser::const_section_iterator>
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();
this->initDOFManager(solver_type);
}
/* -------------------------------------------------------------------------- */
void ModelSolver::initDOFManager(const ID & solver_type) {
if (solver_type == "petsc") {
#if defined(AKANTU_USE_PETSC)
ID id = this->parent_id + ":dof_manager_petsc";
- this->dof_manager = new DOFManagerPETSc(this->mesh, id, this->parent_memory_id);
+ this->dof_manager =
+ new DOFManagerPETSc(this->mesh, id, this->parent_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)
ID id = this->parent_id + ":dof_manager_default";
- this->dof_manager = new DOFManagerDefault(this->mesh, id, this->parent_memory_id);
+ this->dof_manager =
+ new DOFManagerDefault(this->mesh, id, this->parent_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(ID solver_id) {
AKANTU_DEBUG_IN();
if (solver_id == "")
solver_id = default_solver_id;
TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(solver_id);
// make one non linear solve
tss.solveStep(*this);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
-void ModelSolver::getNewSolver(
- const ID & solver_id, const TimeStepSolverType & time_step_solver_type,
- const NonLinearSolverType & non_linear_solver_type) {
+void ModelSolver::getNewSolver(const ID & solver_id,
+ TimeStepSolverType time_step_solver_type,
+ NonLinearSolverType non_linear_solver_type) {
if (this->default_solver_id == "") {
this->default_solver_id = solver_id;
}
+ if (non_linear_solver_type == _nls_auto) {
+ switch (time_step_solver_type) {
+ case _tsst_dynamic:
+ case _tsst_static:
+ non_linear_solver_type = _nls_newton_raphson;
+ break;
+ case _tsst_dynamic_lumped:
+ non_linear_solver_type = _nls_lumped;
+ break;
+ }
+ }
+
NonLinearSolver & nls = this->dof_manager->getNewNonLinearSolver(
solver_id, non_linear_solver_type);
this->dof_manager->getNewTimeStepSolver(solver_id, time_step_solver_type,
nls);
}
/* -------------------------------------------------------------------------- */
void ModelSolver::setIntegrationScheme(
const ID & solver_id, const ID & dof_id,
const IntegrationSchemeType & integration_scheme_type) {
TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(solver_id);
tss.setIntegrationScheme(dof_id, integration_scheme_type);
}
/* -------------------------------------------------------------------------- */
void ModelSolver::predictor() {}
/* -------------------------------------------------------------------------- */
void ModelSolver::corrector() {}
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
// void ModelSolver::setIntegrationScheme(
// const ID & solver_id, const ID & dof_id,
// const IntegrationSchemeType & integration_scheme_type) {}
/* -------------------------------------------------------------------------- */
// /* --------------------------------------------------------------------------
// */
// template <NewmarkBeta::IntegrationSchemeCorrectorType type>
// void SolidMechanicsModel::solve(Array<Real> &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<NewmarkBeta *>(integrator);
// c = nmb_int->getAccelerationCoefficient<type>(time_step);
// d = nmb_int->getVelocityCoefficient<type>(time_step);
// e = nmb_int->getDisplacementCoefficient<type>(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<SolveConvergenceMethod cmethod, SolveConvergenceCriteria criteria>
// 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<cmethod, criteria>(tolerance,
// error, max_iteration, do_not_factorize);
// method=analysis_method;
// return converged;
// }
// /* --------------------------------------------------------------------------
// */
// template<SolveConvergenceMethod cmethod, SolveConvergenceCriteria criteria>
// bool SolidMechanicsModel::solveStep(Real tolerance,
// UInt max_iteration) {
// Real error = 0.;
// return this->template solveStep<cmethod,criteria>(tolerance,
// error,
// max_iteration);
// }
// /* --------------------------------------------------------------------------
// */
// template<SolveConvergenceMethod cmethod, SolveConvergenceCriteria criteria>
// 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<criteria> (tolerance, error);
// if(converged) return converged;
// }
// do {
// if (cmethod == _scm_newton_raphson_tangent)
// this->assembleStiffnessMatrix();
// solve<NewmarkBeta::_displacement_corrector> (*increment, 1.,
// need_factorize);
// this->implicitCorr();
// if(criteria == _scc_residual) this->updateResidual();
// converged = this->testConvergence<criteria> (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<Real> * Ma = new Array<Real>(*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<Real> * Cv = new Array<Real>(*velocity);
// *Cv *= *velocity_damping_matrix;
// *residual -= *Cv;
// delete Cv;
// }
// }
// AKANTU_DEBUG_OUT();
// }
// /* --------------------------------------------------------------------------
// */
// void SolidMechanicsModel::solveLumped(Array<Real> & x,
// const Array<Real> & A,
// const Array<Real> & b,
// const Array<bool> & 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++;
// }
// }
/* -------------------------------------------------------------------------- */
// 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<NewmarkBeta::_acceleration_corrector>(*increment_acceleration);
// }
// AKANTU_DEBUG_OUT();
// }
__END_AKANTU__
diff --git a/src/model/model_solver.hh b/src/model/model_solver.hh
index 3323273c6..9d2e9f647 100644
--- a/src/model/model_solver.hh
+++ b/src/model/model_solver.hh
@@ -1,126 +1,126 @@
/**
* @file model_solver.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "parsable.hh"
#include "solver_callback.hh"
/* -------------------------------------------------------------------------- */
#include <set>
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MODEL_SOLVER_HH__
#define __AKANTU_MODEL_SOLVER_HH__
namespace akantu {
class Mesh;
class DOFManager;
class IntegrationScheme;
}
__BEGIN_AKANTU__
class ModelSolver : public Parsable, public SolverCallback {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
ModelSolver(const Mesh & mesh, const ID & id, UInt memory_id);
virtual ~ModelSolver();
/// initialize the dof manager based on solver type passed in the input file
void initDOFManager();
/// initialize the dof manager based on the used chosen solver type
void initDOFManager(const ID & solver_type);
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// solve a step using a given pre instantiated time step solver and
/// nondynamic linear solver
void solveStep(ID time_step_solver_id = "");
/// Initialize a time solver that can be used afterwards with its id
void
getNewSolver(const ID & solver_id,
- const TimeStepSolverType & time_step_solver_type,
- const NonLinearSolverType & non_linear_solver_type = _nls_auto);
+ TimeStepSolverType time_step_solver_type,
+ NonLinearSolverType non_linear_solver_type = _nls_auto);
/// set an integration scheme for a given dof and a given solver
void
setIntegrationScheme(const ID & solver_id, const ID & dof_id,
const IntegrationSchemeType & integration_scheme_type);
/// set an externally instantiated integration scheme
void setIntegrationScheme(const ID & solver_id, const ID & dof_id,
IntegrationScheme & integration_scheme);
/* ------------------------------------------------------------------------ */
/* SolverCallback interface */
/* ------------------------------------------------------------------------ */
public:
/// Predictor interface for the callback
virtual void predictor();
/// Corrector interface for the callback
virtual void corrector();
/// AssembleResidual interface for the callback
virtual void assembleResidual() = 0;
/// AssembleJacobian interface for the callback
virtual void assembleJacobian() = 0;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
DOFManager & getDOFManager() { return *this->dof_manager; }
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
ID parent_id;
UInt parent_memory_id;
/// Underlying mesh
const Mesh & mesh;
/// Underlying dof_manager (the brain...)
DOFManager * dof_manager;
/// Default time step solver to use
ID default_solver_id;
};
__END_AKANTU__
#endif /* __AKANTU_MODEL_SOLVER_HH__ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc
index 9f9cd9eb2..f212d4fef 100644
--- a/src/model/solid_mechanics/solid_mechanics_model.cc
+++ b/src/model/solid_mechanics/solid_mechanics_model.cc
@@ -1,1462 +1,1463 @@
/**
* @file solid_mechanics_model.cc
*
* @author Ramin Aghababaei <ramin.aghababaei@epfl.ch>
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Clement Roux <clement.roux@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Tue Jul 27 2010
* @date last modification: Tue Jan 19 2016
*
* @brief Implementation of the SolidMechanicsModel class
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_math.hh"
#include "aka_common.hh"
#include "solid_mechanics_model.hh"
#include "group_manager_inline_impl.cc"
#include "element_group.hh"
#include "static_communicator.hh"
#include "dumpable_inline_impl.hh"
#ifdef AKANTU_USE_IOHELPER
#include "dumper_field.hh"
#include "dumper_paraview.hh"
#include "dumper_homogenizing_field.hh"
#include "dumper_internal_material_field.hh"
#include "dumper_elemental_field.hh"
#include "dumper_material_padders.hh"
#include "dumper_element_partition.hh"
#include "dumper_iohelper.hh"
#endif
#ifdef AKANTU_DAMAGE_NON_LOCAL
#include "non_local_manager.hh"
#endif
/* -------------------------------------------------------------------------- */
#include <cmath>
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
const SolidMechanicsModelOptions
default_solid_mechanics_model_options(_explicit_lumped_mass, false);
/* -------------------------------------------------------------------------- */
/**
* A solid mechanics model need a mesh and a dimension to be created. the model
* by it self can not do a lot, the good init functions should be called in
* order to configure the model depending on what we want to do.
*
* @param mesh mesh representing the model we want to simulate
* @param dim spatial dimension of the problem, if dim = 0 (default value) the
* dimension of the problem is assumed to be the on of the mesh
* @param id an id to identify the model
*/
SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh, UInt dim, const ID & id,
const MemoryID & memory_id)
: Model(mesh, dim, id, memory_id), BoundaryCondition<SolidMechanicsModel>(),
time_step(NAN), f_m2a(1.0), displacement(NULL),
previous_displacement(NULL), increment(NULL), mass(NULL), velocity(NULL),
acceleration(NULL), increment_acceleration(NULL), external_force(NULL),
internal_force(NULL), blocked_dofs(NULL), current_position(NULL),
material_index("material index", id),
material_local_numbering("material local numbering", id), materials(0),
material_selector(new DefaultMaterialSelector(material_index)),
is_default_material_selector(true), increment_flag(false),
synch_parallel(NULL), are_materials_instantiated(false) {
AKANTU_DEBUG_IN();
this->createSynchronizerRegistry(this);
this->registerFEEngineObject<MyFEEngineType>("SolidMechanicsFEEngine", mesh,
spatial_dimension);
this->mesh.registerEventHandler(*this);
#if defined(AKANTU_USE_IOHELPER)
this->mesh.registerDumper<DumperParaview>("paraview_all", id, true);
this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular);
#endif
+
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
SolidMechanicsModel::~SolidMechanicsModel() {
AKANTU_DEBUG_IN();
std::vector<Material *>::iterator mat_it;
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
delete *mat_it;
}
materials.clear();
delete synch_parallel;
if (is_default_material_selector) {
delete material_selector;
material_selector = NULL;
}
flatten_internal_map::iterator fl_it = this->registered_internals.begin();
flatten_internal_map::iterator fl_end = this->registered_internals.end();
for (; fl_it != fl_end; ++fl_it) {
delete fl_it->second;
}
#ifdef AKANTU_DAMAGE_NON_LOCAL
delete non_local_manager;
non_local_manager = NULL;
#endif
delete pbc_synch;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::setTimeStep(Real time_step) {
this->time_step = time_step;
#if defined(AKANTU_USE_IOHELPER)
this->mesh.getDumper().setTimeStep(time_step);
#endif
}
/* -------------------------------------------------------------------------- */
/* Initialization */
/* -------------------------------------------------------------------------- */
/**
* This function groups many of the initialization in on function. For most of
* basics case the function should be enough. The functions initialize the
* model, the internal vectors, set them to 0, and depending on the parameters
* it also initialize the explicit or implicit solver.
*
* @param material_file the file containing the materials to use
* @param method the analysis method wanted. See the akantu::AnalysisMethod for
* the different possibilities
*/
void SolidMechanicsModel::initFull(const ModelOptions & options) {
Model::initFull(options);
const SolidMechanicsModelOptions & smm_options =
dynamic_cast<const SolidMechanicsModelOptions &>(options);
this->method = smm_options.analysis_method;
// initialize the vectors
this->initArrays();
// set the initial condition to 0
external_force->clear();
velocity->clear();
acceleration->clear();
displacement->clear();
// initialize pbc
if (this->pbc_pair.size() != 0)
this->initPBC();
// initialize the time integration schemes
// switch (this->method) {
// case _explicit_lumped_mass:
// this->initExplicit();
// break;
// case _explicit_consistent_mass:
// this->initSolver();
// this->initExplicit();
// break;
// case _implicit_dynamic:
// this->initImplicit(true);
// break;
// case _static:
// this->initImplicit(false);
// this->initArraysPreviousDisplacment();
// break;
// default:
// AKANTU_EXCEPTION("analysis method not recognised by SolidMechanicsModel");
// break;
// }
#ifdef AKANTU_DAMAGE_NON_LOCAL
/// create the non-local manager object for non-local damage computations
this->non_local_manager = new NonLocalManager(*this);
#endif
// initialize the materials
if (this->parser->getLastParsedFile() != "") {
this->instantiateMaterials();
}
if (!smm_options.no_init_materials) {
this->initMaterials();
}
if (increment_flag)
this->initBC(*this, *displacement, *increment, *external_force);
else
this->initBC(*this, *displacement, *external_force);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initParallel(MeshPartition * partition,
DataAccessor * data_accessor) {
AKANTU_DEBUG_IN();
if (data_accessor == NULL)
data_accessor = this;
synch_parallel = &createParallelSynch(partition, data_accessor);
synch_registry->registerSynchronizer(*synch_parallel, _gst_material_id);
synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_mass);
synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_stress);
synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_boundary);
synch_registry->registerSynchronizer(*synch_parallel, _gst_for_dump);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initFEEngineBoundary() {
FEEngine & fem_boundary = getFEEngineBoundary();
fem_boundary.initShapeFunctions(_not_ghost);
fem_boundary.initShapeFunctions(_ghost);
fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost);
fem_boundary.computeNormalsOnIntegrationPoints(_ghost);
}
/* -------------------------------------------------------------------------- */
// void SolidMechanicsModel::initExplicit(AnalysisMethod analysis_method) {
// AKANTU_DEBUG_IN();
// // in case of switch from implicit to explicit
// if (!this->isExplicit())
// method = analysis_method;
// // if (integrator) delete integrator;
// // integrator = new CentralDifference();
// UInt nb_nodes = acceleration->getSize();
// UInt nb_degree_of_freedom = acceleration->getNbComponent();
// std::stringstream sstr;
// sstr << id << ":increment_acceleration";
// increment_acceleration =
// &(alloc<Real>(sstr.str(), nb_nodes, nb_degree_of_freedom, Real()));
// AKANTU_DEBUG_OUT();
// }
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initArraysPreviousDisplacment() {
AKANTU_DEBUG_IN();
this->setIncrementFlagOn();
if (not this->previous_displacement) {
UInt nb_nodes = this->mesh.getNbNodes();
std::stringstream sstr_disp_t;
sstr_disp_t << this->id << ":previous_displacement";
this->previous_displacement = &(this->alloc<Real>(
sstr_disp_t.str(), nb_nodes, this->spatial_dimension, 0.));
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Allocate all the needed vectors. By default their are not necessarily set to
* 0
*/
void SolidMechanicsModel::initArrays() {
AKANTU_DEBUG_IN();
UInt nb_nodes = mesh.getNbNodes();
std::stringstream sstr_disp;
sstr_disp << id << ":displacement";
std::stringstream sstr_forc;
sstr_forc << id << ":force";
std::stringstream sstr_ifor;
sstr_forc << id << ":internal_force";
std::stringstream sstr_boun;
sstr_boun << id << ":blocked_dofs";
/* ------------------------------------------------------------------------ */
// for static
this->displacement = &(alloc<Real>(sstr_disp.str(), nb_nodes,
spatial_dimension, REAL_INIT_VALUE));
this->internal_force = &(alloc<Real>(sstr_ifor.str(), nb_nodes,
spatial_dimension, REAL_INIT_VALUE));
this->external_force = &(alloc<Real>(sstr_forc.str(), nb_nodes,
spatial_dimension, REAL_INIT_VALUE));
this->blocked_dofs =
&(alloc<bool>(sstr_boun.str(), nb_nodes, spatial_dimension, false));
this->getDOFManager().registerDOFs("displacements", *this->displacement,
_dst_nodal);
this->getDOFManager().registerBlockedDOFs("displacements",
*this->blocked_dofs);
std::stringstream sstr_curp;
sstr_curp << id << ":current_position";
this->current_position =
&(alloc<Real>(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE));
/* ------------------------------------------------------------------------ */
// for dynamic
std::stringstream sstr_velo;
sstr_velo << id << ":velocity";
std::stringstream sstr_acce;
sstr_acce << id << ":acceleration";
this->velocity = &(alloc<Real>(sstr_velo.str(), nb_nodes, spatial_dimension,
REAL_INIT_VALUE));
this->acceleration = &(alloc<Real>(sstr_acce.str(), nb_nodes,
spatial_dimension, REAL_INIT_VALUE));
this->getDOFManager().registerDOFsDerivative("displacements", 1,
*this->velocity);
this->getDOFManager().registerDOFsDerivative("displacements", 2,
*this->acceleration);
for (UInt g = _not_ghost; g <= _ghost; ++g) {
GhostType gt = (GhostType)g;
Mesh::type_iterator it =
mesh.firstType(spatial_dimension, gt, _ek_not_defined);
Mesh::type_iterator end =
mesh.lastType(spatial_dimension, gt, _ek_not_defined);
for (; it != end; ++it) {
UInt nb_element = mesh.getNbElement(*it, gt);
this->material_index.alloc(nb_element, 1, *it, gt);
this->material_local_numbering.alloc(nb_element, 1, *it, gt);
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Initialize the model,basically it pre-compute the shapes, shapes derivatives
* and jacobian
*
*/
void SolidMechanicsModel::initModel() {
/// \todo add the current position as a parameter to initShapeFunctions for
/// large deformation
getFEEngine().initShapeFunctions(_not_ghost);
getFEEngine().initShapeFunctions(_ghost);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initPBC() {
Model::initPBC();
registerPBCSynchronizer();
// as long as there are ones on the diagonal of the matrix, we can put
// boudandary true for slaves
std::map<UInt, UInt>::iterator it = pbc_pair.begin();
std::map<UInt, UInt>::iterator end = pbc_pair.end();
UInt dim = mesh.getSpatialDimension();
while (it != end) {
for (UInt i = 0; i < dim; ++i)
(*blocked_dofs)((*it).first, i) = true;
++it;
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::registerPBCSynchronizer() {
pbc_synch = new PBCSynchronizer(pbc_pair);
synch_registry->registerSynchronizer(*pbc_synch, _gst_smm_uv);
synch_registry->registerSynchronizer(*pbc_synch, _gst_smm_mass);
synch_registry->registerSynchronizer(*pbc_synch, _gst_smm_res);
synch_registry->registerSynchronizer(*pbc_synch, _gst_for_dump);
// changeLocalEquationNumberForPBC(pbc_pair, mesh.getSpatialDimension());
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::assembleResidual() {
AKANTU_DEBUG_IN();
this->assembleInternalForces();
this->getDOFManager().assembleToResidual("displacements",
*this->external_force, 1);
this->getDOFManager().assembleToResidual("displacements",
*this->internal_force, -1);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* This function computes the internal forces as F_{int} = \int_{\Omega} N
* \sigma d\Omega@f$
*/
void SolidMechanicsModel::assembleInternalForces() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Assemble the internal forces");
this->internal_force->clear();
AKANTU_DEBUG_INFO("Compute local stresses");
std::vector<Material *>::iterator mat_it;
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllStresses(_not_ghost);
}
#ifdef AKANTU_DAMAGE_NON_LOCAL
/* ------------------------------------------------------------------------ */
/* Computation of the non local part */
this->non_local_manager->computeAllNonLocalStresses();
#endif
/* ------------------------------------------------------------------------ */
/* assembling the forces internal */
// communicate the stress
AKANTU_DEBUG_INFO("Send data for residual assembly");
synch_registry->asynchronousSynchronize(_gst_smm_stress);
AKANTU_DEBUG_INFO("Assemble residual for local elements");
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.assembleInternalForces(_not_ghost);
}
AKANTU_DEBUG_INFO("Wait distant stresses");
// finalize communications
synch_registry->waitEndSynchronize(_gst_smm_stress);
AKANTU_DEBUG_INFO("Assemble residual for ghost elements");
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.assembleInternalForces(_ghost);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::assembleJacobian() {
this->assembleStiffnessMatrix();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::assembleStiffnessMatrix() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Assemble the new stiffness matrix.");
// call compute stiffness matrix on each local elements
std::vector<Material *>::iterator mat_it;
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
(*mat_it)->assembleStiffnessMatrix(_not_ghost);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateCurrentPosition() {
AKANTU_DEBUG_IN();
this->current_position->copy(this->mesh.getNodes());
Array<Real>::vector_iterator cpos_it =
this->current_position->begin(spatial_dimension);
Array<Real>::vector_iterator cpos_end =
this->current_position->end(spatial_dimension);
Array<Real>::const_vector_iterator disp_it =
this->displacement->begin(spatial_dimension);
for (; cpos_it != cpos_end; ++cpos_it, ++disp_it) {
*cpos_it += *disp_it;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
// void SolidMechanicsModel::initializeUpdateResidualData() {
// AKANTU_DEBUG_IN();
// UInt nb_nodes = mesh.getNbNodes();
// internal_force->resize(nb_nodes);
// /// copy the forces in residual for boundary conditions
// this->getDOFManager().assembleToResidual("displacements",
// *this->external_force);
// // start synchronization
// synch_registry->asynchronousSynchronize(_gst_smm_uv);
// synch_registry->waitEndSynchronize(_gst_smm_uv);
// this->updateCurrentPosition();
// AKANTU_DEBUG_OUT();
// }
/* -------------------------------------------------------------------------- */
/* Explicit scheme */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
// void SolidMechanicsModel::computeStresses() {
// if (isExplicit()) {
// // start synchronization
// synch_registry->asynchronousSynchronize(_gst_smm_uv);
// synch_registry->waitEndSynchronize(_gst_smm_uv);
// // compute stresses on all local elements for each materials
// std::vector<Material *>::iterator mat_it;
// for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
// Material & mat = **mat_it;
// mat.computeAllStresses(_not_ghost);
// }
// #ifdef AKANTU_DAMAGE_NON_LOCAL
// /* Computation of the non local part */
// this->non_local_manager->computeAllNonLocalStresses();
// #endif
// } else {
// std::vector<Material *>::iterator mat_it;
// for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
// Material & mat = **mat_it;
// mat.computeAllStressesFromTangentModuli(_not_ghost);
// }
// }
// }
/* -------------------------------------------------------------------------- */
/* Implicit scheme */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
// /**
// * Initialize the solver and create the sparse matrices needed.
// *
// */
// void SolidMechanicsModel::initSolver(__attribute__((unused))
// SolverOptions & options) {
// UInt nb_global_nodes = mesh.getNbGlobalNodes();
// jacobian_matrix = &(this->getDOFManager().getNewMatrix("jacobian",
// _symmetric));
// // jacobian_matrix->buildProfile(mesh, *dof_synchronizer,
// spatial_dimension);
// if (!isExplicit()) {
// delete stiffness_matrix;
// std::stringstream sstr_sti;
// sstr_sti << id << ":stiffness_matrix";
// stiffness_matrix = &(this->getDOFManager().getNewMatrix("stiffness",
// _symmetric));
// }
// if (solver) solver->initialize(options);
// }
// /* --------------------------------------------------------------------------
// */
// void SolidMechanicsModel::initJacobianMatrix() {
// // @todo make it more flexible: this is an ugly patch to treat the case of
// non
// // fix profile of the K matrix
// delete jacobian_matrix;
// jacobian_matrix = &(this->getDOFManager().getNewMatrix("jacobian",
// "stiffness"));
// std::stringstream sstr_solv; sstr_solv << id << ":solver";
// delete solver;
// solver = new SolverMumps(*jacobian_matrix, sstr_solv.str());
// if(solver)
// solver->initialize(_solver_no_options);
// }
/* -------------------------------------------------------------------------- */
/**
* Initialize the implicit solver, either for dynamic or static cases,
*
* @param dynamic
*/
// void SolidMechanicsModel::initImplicit(bool dynamic) {
// AKANTU_DEBUG_IN();
// method = dynamic ? _implicit_dynamic : _static;
// if (!increment)
// setIncrementFlagOn();
// initSolver();
// // if(method == _implicit_dynamic) {
// // if(integrator) delete integrator;
// // integrator = new TrapezoidalRule2();
// // }
// AKANTU_DEBUG_OUT();
// }
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
// SparseMatrix & SolidMechanicsModel::initVelocityDampingMatrix() {
// return this->getDOFManager().getNewMatrix("velocity_damping", "jacobian");
// }
// /* --------------------------------------------------------------------------
// */
// void SolidMechanicsModel::implicitPred() {
// AKANTU_DEBUG_IN();
// if(previous_displacement) previous_displacement->copy(*displacement);
// if(method == _implicit_dynamic)
// integrator->integrationSchemePred(time_step,
// *displacement,
// *velocity,
// *acceleration,
// *blocked_dofs);
// AKANTU_DEBUG_OUT();
// }
// /* --------------------------------------------------------------------------
// */
// void SolidMechanicsModel::implicitCorr() {
// AKANTU_DEBUG_IN();
// if(method == _implicit_dynamic) {
// integrator->integrationSchemeCorrDispl(time_step,
// *displacement,
// *velocity,
// *acceleration,
// *blocked_dofs,
// *increment);
// } else {
// UInt nb_nodes = displacement->getSize();
// UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes;
// Real * incr_val = increment->storage();
// Real * disp_val = displacement->storage();
// bool * boun_val = blocked_dofs->storage();
// for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++disp_val, ++incr_val,
// ++boun_val){
// *incr_val *= (1. - *boun_val);
// *disp_val += *incr_val;
// }
// }
// AKANTU_DEBUG_OUT();
// }
/* -------------------------------------------------------------------------- */
// void SolidMechanicsModel::updateIncrement() {
// AKANTU_DEBUG_IN();
// AKANTU_DEBUG_ASSERT(
// previous_displacement,
// "The previous displacement has to be initialized."
// << " Are you working with Finite or Ineslactic deformations?");
// UInt nb_nodes = displacement->getSize();
// UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes;
// Real * incr_val = increment->storage();
// Real * disp_val = displacement->storage();
// Real * prev_disp_val = previous_displacement->storage();
// for (UInt j = 0; j < nb_degree_of_freedom;
// ++j, ++disp_val, ++incr_val, ++prev_disp_val)
// *incr_val = (*disp_val - *prev_disp_val);
// AKANTU_DEBUG_OUT();
// }
/* -------------------------------------------------------------------------- */
// void SolidMechanicsModel::updatePreviousDisplacement() {
// AKANTU_DEBUG_IN();
// AKANTU_DEBUG_ASSERT(
// previous_displacement,
// "The previous displacement has to be initialized."
// << " Are you working with Finite or Ineslactic deformations?");
// previous_displacement->copy(*displacement);
// AKANTU_DEBUG_OUT();
// }
/* -------------------------------------------------------------------------- */
/* Information */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::synchronizeBoundaries() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(synch_registry,
"Synchronizer registry was not initialized."
<< " Did you call initParallel?");
synch_registry->synchronize(_gst_smm_boundary);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::synchronizeResidual() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(synch_registry,
"Synchronizer registry was not initialized."
<< " Did you call initPBC?");
synch_registry->synchronize(_gst_smm_res);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::setIncrementFlagOn() {
AKANTU_DEBUG_IN();
if (!increment) {
UInt nb_nodes = mesh.getNbNodes();
std::stringstream sstr_inc;
sstr_inc << id << ":increment";
increment = &(alloc<Real>(sstr_inc.str(), nb_nodes, spatial_dimension, 0.));
}
increment_flag = true;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getStableTimeStep() {
AKANTU_DEBUG_IN();
Real min_dt = getStableTimeStep(_not_ghost);
/// reduction min over all processors
StaticCommunicator::getStaticCommunicator().allReduce(&min_dt, 1, _so_min);
AKANTU_DEBUG_OUT();
return min_dt;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getStableTimeStep(const GhostType & ghost_type) {
AKANTU_DEBUG_IN();
Material ** mat_val = &(materials.at(0));
Real min_dt = std::numeric_limits<Real>::max();
this->updateCurrentPosition();
Element elem;
elem.ghost_type = ghost_type;
elem.kind = _ek_regular;
Mesh::type_iterator it =
mesh.firstType(spatial_dimension, ghost_type, _ek_regular);
Mesh::type_iterator end =
mesh.lastType(spatial_dimension, ghost_type, _ek_regular);
for (; it != end; ++it) {
elem.type = *it;
UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*it);
UInt nb_element = mesh.getNbElement(*it);
Array<UInt>::const_scalar_iterator mat_indexes =
material_index(*it, ghost_type).begin();
Array<UInt>::const_scalar_iterator mat_loc_num =
material_local_numbering(*it, ghost_type).begin();
Array<Real> X(0, nb_nodes_per_element * spatial_dimension);
FEEngine::extractNodalToElementField(mesh, *current_position, X, *it,
_not_ghost);
Array<Real>::matrix_iterator X_el =
X.begin(spatial_dimension, nb_nodes_per_element);
for (UInt el = 0; el < nb_element;
++el, ++X_el, ++mat_indexes, ++mat_loc_num) {
elem.element = *mat_loc_num;
Real el_h = getFEEngine().getElementInradius(*X_el, *it);
Real el_c = mat_val[*mat_indexes]->getCelerity(elem);
Real el_dt = el_h / el_c;
min_dt = std::min(min_dt, el_dt);
}
}
AKANTU_DEBUG_OUT();
return min_dt;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getKineticEnergy() {
AKANTU_DEBUG_IN();
if (!mass && !mass_matrix)
AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix.");
Real ekin = 0.;
UInt nb_nodes = mesh.getNbNodes();
Real * vel_val = velocity->storage();
Real * mass_val = mass->storage();
for (UInt n = 0; n < nb_nodes; ++n) {
Real mv2 = 0;
bool is_local_node = mesh.isLocalOrMasterNode(n);
bool is_not_pbc_slave_node = !isPBCSlaveNode(n);
bool count_node = is_local_node && is_not_pbc_slave_node;
for (UInt i = 0; i < spatial_dimension; ++i) {
if (count_node)
mv2 += *vel_val * *vel_val * *mass_val;
vel_val++;
mass_val++;
}
ekin += mv2;
}
StaticCommunicator::getStaticCommunicator().allReduce(&ekin, 1, _so_sum);
AKANTU_DEBUG_OUT();
return ekin * .5;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getKineticEnergy(const ElementType & type,
UInt index) {
AKANTU_DEBUG_IN();
UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type);
Array<Real> vel_on_quad(nb_quadrature_points, spatial_dimension);
Array<UInt> filter_element(1, 1, index);
getFEEngine().interpolateOnIntegrationPoints(*velocity, vel_on_quad,
spatial_dimension, type,
_not_ghost, filter_element);
Array<Real>::vector_iterator vit = vel_on_quad.begin(spatial_dimension);
Array<Real>::vector_iterator vend = vel_on_quad.end(spatial_dimension);
Vector<Real> rho_v2(nb_quadrature_points);
Real rho = materials[material_index(type)(index)]->getRho();
for (UInt q = 0; vit != vend; ++vit, ++q) {
rho_v2(q) = rho * vit->dot(*vit);
}
AKANTU_DEBUG_OUT();
return .5 * getFEEngine().integrate(rho_v2, type, index);
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getExternalWork() {
AKANTU_DEBUG_IN();
Real * incr_or_velo = NULL;
if (this->method == _static) {
incr_or_velo = this->increment->storage();
} else {
incr_or_velo = this->velocity->storage();
}
Real * forc = external_force->storage();
Real * resi = internal_force->storage();
bool * boun = blocked_dofs->storage();
Real work = 0.;
UInt nb_nodes = this->mesh.getNbNodes();
for (UInt n = 0; n < nb_nodes; ++n) {
bool is_local_node = this->mesh.isLocalOrMasterNode(n);
bool is_not_pbc_slave_node = !this->isPBCSlaveNode(n);
bool count_node = is_local_node && is_not_pbc_slave_node;
for (UInt i = 0; i < this->spatial_dimension; ++i) {
if (count_node) {
if (*boun)
work -= *resi * *incr_or_velo;
else
work += *forc * *incr_or_velo;
}
++incr_or_velo;
++forc;
++resi;
++boun;
}
}
StaticCommunicator::getStaticCommunicator().allReduce(&work, 1, _so_sum);
if (this->method != _static)
work *= this->time_step;
AKANTU_DEBUG_OUT();
return work;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getEnergy(const std::string & energy_id) {
AKANTU_DEBUG_IN();
if (energy_id == "kinetic") {
return getKineticEnergy();
} else if (energy_id == "external work") {
return getExternalWork();
}
Real energy = 0.;
std::vector<Material *>::iterator mat_it;
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
energy += (*mat_it)->getEnergy(energy_id);
}
/// reduction sum over all processors
StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum);
AKANTU_DEBUG_OUT();
return energy;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getEnergy(const std::string & energy_id,
const ElementType & type, UInt index) {
AKANTU_DEBUG_IN();
if (energy_id == "kinetic") {
return getKineticEnergy(type, index);
}
std::vector<Material *>::iterator mat_it;
UInt mat_index = this->material_index(type, _not_ghost)(index);
UInt mat_loc_num = this->material_local_numbering(type, _not_ghost)(index);
Real energy =
this->materials[mat_index]->getEnergy(energy_id, type, mat_loc_num);
AKANTU_DEBUG_OUT();
return energy;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onElementsAdded(const Array<Element> & element_list,
const NewElementsEvent & event) {
AKANTU_DEBUG_IN();
this->getFEEngine().initShapeFunctions(_not_ghost);
this->getFEEngine().initShapeFunctions(_ghost);
for (UInt g = _not_ghost; g <= _ghost; ++g) {
GhostType gt = (GhostType)g;
Mesh::type_iterator it =
this->mesh.firstType(spatial_dimension, gt, _ek_not_defined);
Mesh::type_iterator end =
this->mesh.lastType(spatial_dimension, gt, _ek_not_defined);
for (; it != end; ++it) {
UInt nb_element = this->mesh.getNbElement(*it, gt);
if (!material_index.exists(*it, gt)) {
this->material_index.alloc(nb_element, 1, *it, gt);
this->material_local_numbering.alloc(nb_element, 1, *it, gt);
} else {
this->material_index(*it, gt).resize(nb_element);
this->material_local_numbering(*it, gt).resize(nb_element);
}
}
}
Array<Element>::const_iterator<Element> it = element_list.begin();
Array<Element>::const_iterator<Element> end = element_list.end();
ElementTypeMapArray<UInt> filter("new_element_filter", this->getID());
for (UInt el = 0; it != end; ++it, ++el) {
const Element & elem = *it;
if (!filter.exists(elem.type, elem.ghost_type))
filter.alloc(0, 1, elem.type, elem.ghost_type);
filter(elem.type, elem.ghost_type).push_back(elem.element);
}
this->assignMaterialToElements(&filter);
std::vector<Material *>::iterator mat_it;
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
(*mat_it)->onElementsAdded(element_list, event);
}
if (method == _explicit_lumped_mass)
this->assembleMassLumped();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onElementsRemoved(
__attribute__((unused)) const Array<Element> & element_list,
const ElementTypeMapArray<UInt> & new_numbering,
const RemovedElementsEvent & event) {
this->getFEEngine().initShapeFunctions(_not_ghost);
this->getFEEngine().initShapeFunctions(_ghost);
std::vector<Material *>::iterator mat_it;
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
(*mat_it)->onElementsRemoved(element_list, new_numbering, event);
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onNodesAdded(const Array<UInt> & nodes_list,
__attribute__((unused))
const NewNodesEvent & event) {
AKANTU_DEBUG_IN();
UInt nb_nodes = mesh.getNbNodes();
if (displacement)
displacement->resize(nb_nodes);
if (mass)
mass->resize(nb_nodes);
if (velocity)
velocity->resize(nb_nodes);
if (acceleration)
acceleration->resize(nb_nodes);
if (external_force)
external_force->resize(nb_nodes);
if (internal_force)
internal_force->resize(nb_nodes);
if (blocked_dofs)
blocked_dofs->resize(nb_nodes);
if (previous_displacement)
previous_displacement->resize(nb_nodes);
if (increment_acceleration)
increment_acceleration->resize(nb_nodes);
if (increment)
increment->resize(nb_nodes);
std::vector<Material *>::iterator mat_it;
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
(*mat_it)->onNodesAdded(nodes_list, event);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onNodesRemoved(__attribute__((unused))
const Array<UInt> & element_list,
const Array<UInt> & new_numbering,
__attribute__((unused))
const RemovedNodesEvent & event) {
if (displacement)
mesh.removeNodesFromArray(*displacement, new_numbering);
if (mass)
mesh.removeNodesFromArray(*mass, new_numbering);
if (velocity)
mesh.removeNodesFromArray(*velocity, new_numbering);
if (acceleration)
mesh.removeNodesFromArray(*acceleration, new_numbering);
if (internal_force)
mesh.removeNodesFromArray(*internal_force, new_numbering);
if (external_force)
mesh.removeNodesFromArray(*external_force, new_numbering);
if (blocked_dofs)
mesh.removeNodesFromArray(*blocked_dofs, new_numbering);
// if (increment_acceleration)
// mesh.removeNodesFromArray(*increment_acceleration, new_numbering);
// if (increment)
// mesh.removeNodesFromArray(*increment, new_numbering);
// if (method != _explicit_lumped_mass) {
// this->initSolver();
// }
}
/* -------------------------------------------------------------------------- */
bool SolidMechanicsModel::isInternal(const std::string & field_name,
const ElementKind & element_kind) {
bool is_internal = false;
/// check if at least one material contains field_id as an internal
for (UInt m = 0; m < materials.size() && !is_internal; ++m) {
is_internal |= materials[m]->isInternal<Real>(field_name, element_kind);
}
return is_internal;
}
/* -------------------------------------------------------------------------- */
ElementTypeMap<UInt>
SolidMechanicsModel::getInternalDataPerElem(const std::string & field_name,
const ElementKind & element_kind) {
if (!(this->isInternal(field_name, element_kind)))
AKANTU_EXCEPTION("unknown internal " << field_name);
for (UInt m = 0; m < materials.size(); ++m) {
if (materials[m]->isInternal<Real>(field_name, element_kind))
return materials[m]->getInternalDataPerElem<Real>(field_name,
element_kind);
}
return ElementTypeMap<UInt>();
}
/* -------------------------------------------------------------------------- */
ElementTypeMapArray<Real> &
SolidMechanicsModel::flattenInternal(const std::string & field_name,
const ElementKind & kind,
const GhostType ghost_type) {
std::pair<std::string, ElementKind> key(field_name, kind);
if (this->registered_internals.count(key) == 0) {
this->registered_internals[key] =
new ElementTypeMapArray<Real>(field_name, this->id);
}
ElementTypeMapArray<Real> * internal_flat = this->registered_internals[key];
typedef ElementTypeMapArray<Real>::type_iterator iterator;
iterator tit = internal_flat->firstType(spatial_dimension, ghost_type, kind);
iterator end = internal_flat->lastType(spatial_dimension, ghost_type, kind);
for (; tit != end; ++tit) {
ElementType type = *tit;
(*internal_flat)(type, ghost_type).clear();
}
for (UInt m = 0; m < materials.size(); ++m) {
if (materials[m]->isInternal<Real>(field_name, kind))
materials[m]->flattenInternal(field_name, *internal_flat, ghost_type,
kind);
}
return *internal_flat;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::flattenAllRegisteredInternals(
const ElementKind & kind) {
typedef std::map<std::pair<std::string, ElementKind>,
ElementTypeMapArray<Real> *>::iterator iterator;
iterator it = this->registered_internals.begin();
iterator end = this->registered_internals.end();
while (it != end) {
if (kind == it->first.second)
this->flattenInternal(it->first.first, kind);
++it;
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onDump() {
this->flattenAllRegisteredInternals(_ek_regular);
}
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
dumper::Field * SolidMechanicsModel::createElementalField(
const std::string & field_name, const std::string & group_name,
bool padding_flag, const UInt & spatial_dimension,
const ElementKind & kind) {
dumper::Field * field = NULL;
if (field_name == "partitions")
field = mesh.createElementalField<UInt, dumper::ElementPartitionField>(
mesh.getConnectivities(), group_name, spatial_dimension, kind);
else if (field_name == "material_index")
field = mesh.createElementalField<UInt, Vector, dumper::ElementalField>(
material_index, group_name, spatial_dimension, kind);
else {
// this copy of field_name is used to compute derivated data such as
// strain and von mises stress that are based on grad_u and stress
std::string field_name_copy(field_name);
if (field_name == "strain" || field_name == "Green strain" ||
field_name == "principal strain" ||
field_name == "principal Green strain")
field_name_copy = "grad_u";
else if (field_name == "Von Mises stress")
field_name_copy = "stress";
bool is_internal = this->isInternal(field_name_copy, kind);
if (is_internal) {
ElementTypeMap<UInt> nb_data_per_elem =
this->getInternalDataPerElem(field_name_copy, kind);
ElementTypeMapArray<Real> & internal_flat =
this->flattenInternal(field_name_copy, kind);
field = mesh.createElementalField<Real, dumper::InternalMaterialField>(
internal_flat, group_name, spatial_dimension, kind, nb_data_per_elem);
if (field_name == "strain") {
dumper::ComputeStrain<false> * foo =
new dumper::ComputeStrain<false>(*this);
field = dumper::FieldComputeProxy::createFieldCompute(field, *foo);
} else if (field_name == "Von Mises stress") {
dumper::ComputeVonMisesStress * foo =
new dumper::ComputeVonMisesStress(*this);
field = dumper::FieldComputeProxy::createFieldCompute(field, *foo);
} else if (field_name == "Green strain") {
dumper::ComputeStrain<true> * foo =
new dumper::ComputeStrain<true>(*this);
field = dumper::FieldComputeProxy::createFieldCompute(field, *foo);
} else if (field_name == "principal strain") {
dumper::ComputePrincipalStrain<false> * foo =
new dumper::ComputePrincipalStrain<false>(*this);
field = dumper::FieldComputeProxy::createFieldCompute(field, *foo);
} else if (field_name == "principal Green strain") {
dumper::ComputePrincipalStrain<true> * foo =
new dumper::ComputePrincipalStrain<true>(*this);
field = dumper::FieldComputeProxy::createFieldCompute(field, *foo);
}
// treat the paddings
if (padding_flag) {
if (field_name == "stress") {
if (spatial_dimension == 2) {
dumper::StressPadder<2> * foo = new dumper::StressPadder<2>(*this);
field = dumper::FieldComputeProxy::createFieldCompute(field, *foo);
}
} else if (field_name == "strain" || field_name == "Green strain") {
if (spatial_dimension == 2) {
dumper::StrainPadder<2> * foo = new dumper::StrainPadder<2>(*this);
field = dumper::FieldComputeProxy::createFieldCompute(field, *foo);
}
}
}
// homogenize the field
dumper::ComputeFunctorInterface * foo =
dumper::HomogenizerProxy::createHomogenizer(*field);
field = dumper::FieldComputeProxy::createFieldCompute(field, *foo);
}
}
return field;
}
/* -------------------------------------------------------------------------- */
dumper::Field *
SolidMechanicsModel::createNodalFieldReal(const std::string & field_name,
const std::string & group_name,
bool padding_flag) {
std::map<std::string, Array<Real> *> real_nodal_fields;
real_nodal_fields["displacement"] = this->displacement;
real_nodal_fields["mass"] = this->mass;
real_nodal_fields["velocity"] = this->velocity;
real_nodal_fields["acceleration"] = this->acceleration;
real_nodal_fields["force"] = this->external_force;
real_nodal_fields["residual"] = this->internal_force;
real_nodal_fields["increment"] = this->increment;
dumper::Field * field = NULL;
if (padding_flag)
field = this->mesh.createNodalField(real_nodal_fields[field_name],
group_name, 3);
else
field =
this->mesh.createNodalField(real_nodal_fields[field_name], group_name);
return field;
}
/* -------------------------------------------------------------------------- */
dumper::Field *
SolidMechanicsModel::createNodalFieldBool(const std::string & field_name,
const std::string & group_name,
bool padding_flag) {
std::map<std::string, Array<bool> *> uint_nodal_fields;
uint_nodal_fields["blocked_dofs"] = blocked_dofs;
dumper::Field * field = NULL;
field = mesh.createNodalField(uint_nodal_fields[field_name], group_name);
return field;
}
/* -------------------------------------------------------------------------- */
#else
/* -------------------------------------------------------------------------- */
dumper::Field * SolidMechanicsModel::createElementalField(
__attribute__((unused)) const std::string & field_name,
__attribute__((unused)) const std::string & group_name,
__attribute__((unused)) bool padding_flag,
__attribute__((unused)) const UInt & spatial_dimension,
__attribute__((unused)) const ElementKind & kind) {
return NULL;
}
/* -------------------------------------------------------------------------- */
dumper::Field * SolidMechanicsModel::createNodalFieldReal(
__attribute__((unused)) const std::string & field_name,
__attribute__((unused)) const std::string & group_name,
__attribute__((unused)) bool padding_flag) {
return NULL;
}
/* -------------------------------------------------------------------------- */
dumper::Field * SolidMechanicsModel::createNodalFieldBool(
__attribute__((unused)) const std::string & field_name,
__attribute__((unused)) const std::string & group_name,
__attribute__((unused)) bool padding_flag) {
return NULL;
}
#endif
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(const std::string & dumper_name) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
mesh.dump(dumper_name);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(const std::string & dumper_name, UInt step) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
mesh.dump(dumper_name, step);
}
/* ------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(const std::string & dumper_name, Real time,
UInt step) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
mesh.dump(dumper_name, time, step);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump() {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
mesh.dump();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(UInt step) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
mesh.dump(step);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(Real time, UInt step) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
mesh.dump(time, step);
}
/* -------------------------------------------------------------------------- */
// void SolidMechanicsModel::computeCauchyStresses() {
// AKANTU_DEBUG_IN();
// // call compute stiffness matrix on each local elements
// std::vector<Material *>::iterator mat_it;
// for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
// Material & mat = **mat_it;
// if (mat.isFiniteDeformation())
// mat.computeAllCauchyStresses(_not_ghost);
// }
// AKANTU_DEBUG_OUT();
// }
/* -------------------------------------------------------------------------- */
// void SolidMechanicsModel::saveStressAndStrainBeforeDamage() {
// EventManager::sendEvent(
// SolidMechanicsModelEvent::BeginningOfDamageIterationEvent());
// }
/* -------------------------------------------------------------------------- */
// void SolidMechanicsModel::updateEnergiesAfterDamage() {
// EventManager::sendEvent(SolidMechanicsModelEvent::AfterDamageEvent());
// }
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::printself(std::ostream & stream, int indent) const {
std::string space;
for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
;
stream << space << "Solid Mechanics Model [" << std::endl;
stream << space << " + id : " << id << std::endl;
stream << space << " + spatial dimension : " << spatial_dimension
<< std::endl;
stream << space << " + fem [" << std::endl;
getFEEngine().printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << " + nodals information [" << std::endl;
displacement->printself(stream, indent + 2);
mass->printself(stream, indent + 2);
velocity->printself(stream, indent + 2);
acceleration->printself(stream, indent + 2);
external_force->printself(stream, indent + 2);
internal_force->printself(stream, indent + 2);
blocked_dofs->printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << " + material information [" << std::endl;
material_index.printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << " + materials [" << std::endl;
std::vector<Material *>::const_iterator mat_it;
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
const Material & mat = *(*mat_it);
mat.printself(stream, indent + 1);
}
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
__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 7466219b0..fecfd5c4f 100644
--- a/src/model/time_step_solvers/time_step_solver_default.cc
+++ b/src/model/time_step_solvers/time_step_solver_default.cc
@@ -1,242 +1,247 @@
/**
* @file time_step_solver_default.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "time_step_solver_default.hh"
#include "dof_manager_default.hh"
#include "sparse_matrix_aij.hh"
+#include "pseudo_time.hh"
#include "integration_scheme_1st_order.hh"
#include "integration_scheme_2nd_order.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
TimeStepSolverDefault::TimeStepSolverDefault(
DOFManagerDefault & dof_manager, const TimeStepSolverType & type,
NonLinearSolver & non_linear_solver, const ID & id, UInt memory_id)
: TimeStepSolver(dof_manager, type, non_linear_solver, id, memory_id),
dof_manager(dof_manager), is_mass_lumped(false) {
switch (type) {
case _tsst_dynamic:
break;
case _tsst_dynamic_lumped:
this->is_mass_lumped = true;
break;
case _tsst_static:
/// initialize a static time solver for allback dofs
break;
}
}
/* -------------------------------------------------------------------------- */
void TimeStepSolverDefault::setIntegrationScheme(
const ID & dof_id, const IntegrationSchemeType & type) {
if (this->integration_schemes.find(dof_id) !=
this->integration_schemes.end()) {
AKANTU_EXCEPTION("Their DOFs "
<< dof_id
<< " have already an integration scheme associated");
}
IntegrationScheme * integration_scheme = NULL;
if (this->is_mass_lumped) {
switch (type) {
case _ist_forward_euler: {
integration_scheme = new ForwardEuler(dof_manager, dof_id);
break;
}
case _ist_central_difference: {
integration_scheme = new CentralDifference(dof_manager, dof_id);
break;
}
default:
AKANTU_EXCEPTION(
"This integration scheme cannot be used in lumped dynamic");
}
} else {
switch (type) {
+ case _ist_pseudo_time: {
+ integration_scheme = new PseudoTime(dof_manager, dof_id);
+ break;
+ }
case _ist_forward_euler: {
integration_scheme = new ForwardEuler(dof_manager, dof_id);
break;
}
case _ist_trapezoidal_rule_1: {
integration_scheme = new TrapezoidalRule1(dof_manager, dof_id);
break;
}
case _ist_backward_euler: {
integration_scheme = new BackwardEuler(dof_manager, dof_id);
break;
}
case _ist_central_difference: {
integration_scheme = new CentralDifference(dof_manager, dof_id);
break;
}
case _ist_fox_goodwin: {
integration_scheme = new FoxGoodwin(dof_manager, dof_id);
break;
}
case _ist_trapezoidal_rule_2: {
integration_scheme = new TrapezoidalRule2(dof_manager, dof_id);
break;
}
case _ist_linear_acceleration: {
integration_scheme = new LinearAceleration(dof_manager, dof_id);
break;
}
// Write a c++11 version of the constructor with initializer list that
// contains the arguments for the integration scheme
case _ist_generalized_trapezoidal:
case _ist_newmark_beta:
AKANTU_EXCEPTION(
"This time step solvers cannot be created with this constructor");
}
}
this->integration_schemes[dof_id] = integration_scheme;
this->integration_schemes_owner.insert(dof_id);
}
/* -------------------------------------------------------------------------- */
TimeStepSolverDefault::~TimeStepSolverDefault() {
DOFsIntegrationSchemesOwner::iterator it =
this->integration_schemes_owner.begin();
DOFsIntegrationSchemesOwner::iterator end =
this->integration_schemes_owner.end();
for (; it != end; ++it) {
delete this->integration_schemes[*it];
}
this->integration_schemes.clear();
}
/* -------------------------------------------------------------------------- */
void TimeStepSolverDefault::solveStep(SolverCallback & solver_callback) {
this->solver_callback = &solver_callback;
this->non_linear_solver.solve(*this);
this->solver_callback = NULL;
}
/* -------------------------------------------------------------------------- */
void TimeStepSolverDefault::predictor() {
AKANTU_DEBUG_IN();
TimeStepSolver::predictor();
DOFsIntegrationSchemes::iterator integration_scheme_it =
this->integration_schemes.begin();
DOFsIntegrationSchemes::iterator integration_scheme_end =
this->integration_schemes.end();
for (; integration_scheme_it != integration_scheme_end;
++integration_scheme_it) {
integration_scheme_it->second->predictor(this->time_step);
}
// UInt nb_degree_of_freedom = u.getSize() * u.getNbComponent();
// Array<Real>::scalar_iterator incr_it =
// increment.begin_reinterpret(nb_degree_of_freedom);
// Array<Real>::const_scalar_iterator u_it =
// u.begin_reinterpret(nb_degree_of_freedom);
// Array<Real>::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();
TimeStepSolver::corrector();
DOFsIntegrationSchemes::iterator integration_scheme_it =
this->integration_schemes.begin();
DOFsIntegrationSchemes::iterator integration_scheme_end =
this->integration_schemes.end();
for (; integration_scheme_it != integration_scheme_end;
++integration_scheme_it) {
integration_scheme_it->second->corrector(
IntegrationScheme::SolutionType(this->solution_type), this->time_step);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void TimeStepSolverDefault::assembleJacobian() {
AKANTU_DEBUG_IN();
TimeStepSolver::assembleJacobian();
DOFsIntegrationSchemes::iterator integration_scheme_it =
this->integration_schemes.begin();
DOFsIntegrationSchemes::iterator integration_scheme_end =
this->integration_schemes.end();
for (; integration_scheme_it != integration_scheme_end;
++integration_scheme_it) {
integration_scheme_it->second->assembleJacobian(
IntegrationScheme::SolutionType(this->solution_type), this->time_step);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void TimeStepSolverDefault::assembleResidual() {
AKANTU_DEBUG_IN();
TimeStepSolver::assembleResidual();
DOFsIntegrationSchemes::iterator integration_scheme_it =
this->integration_schemes.begin();
DOFsIntegrationSchemes::iterator integration_scheme_end =
this->integration_schemes.end();
for (; integration_scheme_it != integration_scheme_end;
++integration_scheme_it) {
integration_scheme_it->second->assembleResidual(this->is_mass_lumped);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
__END_AKANTU__