diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh index 0c33fabbf..7bda2cc7b 100644 --- a/src/common/aka_common.hh +++ b/src/common/aka_common.hh @@ -1,648 +1,662 @@ /** * @file aka_common.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 * @date last modification: Mon Feb 12 2018 * * @brief common type descriptions for akantu * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @section DESCRIPTION * * All common things to be included in the projects files * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COMMON_HH__ #define __AKANTU_COMMON_HH__ #include "aka_compatibilty_with_cpp_standard.hh" /* -------------------------------------------------------------------------- */ #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" /* -------------------------------------------------------------------------- */ #include #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /* Constants */ /* -------------------------------------------------------------------------- */ namespace { __attribute__((unused)) constexpr UInt _all_dimensions{ std::numeric_limits::max()}; #ifdef AKANTU_NDEBUG __attribute__((unused)) constexpr Real REAL_INIT_VALUE{0.}; #else __attribute__((unused)) constexpr Real REAL_INIT_VALUE{ std::numeric_limits::quiet_NaN()}; #endif } // namespace /* -------------------------------------------------------------------------- */ /* Common types */ /* -------------------------------------------------------------------------- */ using ID = std::string; using MemoryID = UInt; } // namespace akantu /* -------------------------------------------------------------------------- */ #include "aka_enum_macros.hh" /* -------------------------------------------------------------------------- */ #include "aka_element_classes_info.hh" namespace akantu { /* -------------------------------------------------------------------------- */ /* Mesh/FEM/Model types */ /* -------------------------------------------------------------------------- */ /// small help to use names for directions enum SpatialDirection { _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 MeshEventHandlerPriority defines relative order of execution of /// events enum EventHandlerPriority { _ehp_highest = 0, _ehp_mesh = 5, _ehp_fe_engine = 9, _ehp_synchronizer = 10, _ehp_dof_manager = 20, _ehp_model = 94, _ehp_non_local_manager = 100, _ehp_lowest = 100 }; // clang-format off #define AKANTU_MODEL_TYPES \ (model) \ (solid_mechanics_model) \ (solid_mechanics_model_cohesive) \ (heat_transfer_model) \ (structural_mechanics_model) \ (embedded_model) \ (phase_field_model) \ (coupler_solid_phasefield) // clang-format on /// enum ModelType defines which type of physics is solved AKANTU_CLASS_ENUM_DECLARE(ModelType, AKANTU_MODEL_TYPES) AKANTU_CLASS_ENUM_OUTPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) AKANTU_CLASS_ENUM_INPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) /// 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 }; #if !defined(DOXYGEN) // clang-format off #define AKANTU_NON_LINEAR_SOLVER_TYPES \ (linear) \ (newton_raphson) \ (newton_raphson_modified) \ (lumped) \ (gmres) \ (bfgs) \ (cg) \ (auto) // clang-format on AKANTU_CLASS_ENUM_DECLARE(NonLinearSolverType, AKANTU_NON_LINEAR_SOLVER_TYPES) AKANTU_CLASS_ENUM_OUTPUT_STREAM(NonLinearSolverType, AKANTU_NON_LINEAR_SOLVER_TYPES) AKANTU_CLASS_ENUM_INPUT_STREAM(NonLinearSolverType, AKANTU_NON_LINEAR_SOLVER_TYPES) #else /// Type of non linear resolution available in akantu enum class NonLinearSolverType { _linear, ///< No non linear convergence loop _newton_raphson, ///< Regular Newton-Raphson _newton_raphson_modified, ///< Newton-Raphson with initial tangent _lumped, ///< Case of lumped mass or equivalent matrix _gmres, _bfgs, _cg, _auto ///< This will take a default value that make sense in case of /// model::getNewSolver }; #endif #if !defined(DOXYGEN) // clang-format off #define AKANTU_TIME_STEP_SOLVER_TYPE \ (static) \ (dynamic) \ (dynamic_lumped) \ (not_defined) // clang-format on AKANTU_CLASS_ENUM_DECLARE(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE) AKANTU_CLASS_ENUM_OUTPUT_STREAM(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE) AKANTU_CLASS_ENUM_INPUT_STREAM(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE) #else /// Type of time stepping solver enum class TimeStepSolverType { _static, ///< Static solution _dynamic, ///< Dynamic solver _dynamic_lumped, ///< Dynamic solver with lumped mass _not_defined, ///< For not defined cases }; #endif #if !defined(DOXYGEN) // clang-format off #define AKANTU_INTEGRATION_SCHEME_TYPE \ (pseudo_time) \ (forward_euler) \ (trapezoidal_rule_1) \ (backward_euler) \ (central_difference) \ (fox_goodwin) \ (trapezoidal_rule_2) \ (linear_acceleration) \ (newmark_beta) \ (generalized_trapezoidal) // clang-format on AKANTU_CLASS_ENUM_DECLARE(IntegrationSchemeType, AKANTU_INTEGRATION_SCHEME_TYPE) AKANTU_CLASS_ENUM_OUTPUT_STREAM(IntegrationSchemeType, AKANTU_INTEGRATION_SCHEME_TYPE) AKANTU_CLASS_ENUM_INPUT_STREAM(IntegrationSchemeType, AKANTU_INTEGRATION_SCHEME_TYPE) #else /// Type of integration scheme enum class IntegrationSchemeType { _pseudo_time, ///< Pseudo Time _forward_euler, ///< GeneralizedTrapezoidal(0) _trapezoidal_rule_1, ///< GeneralizedTrapezoidal(1/2) _backward_euler, ///< GeneralizedTrapezoidal(1) _central_difference, ///< NewmarkBeta(0, 1/2) _fox_goodwin, ///< NewmarkBeta(1/6, 1/2) _trapezoidal_rule_2, ///< NewmarkBeta(1/2, 1/2) _linear_acceleration, ///< NewmarkBeta(1/3, 1/2) _newmark_beta, ///< generic NewmarkBeta with user defined /// alpha and beta _generalized_trapezoidal ///< generic GeneralizedTrapezoidal with user /// defined alpha }; #endif #if !defined(DOXYGEN) // clang-format off #define AKANTU_SOLVE_CONVERGENCE_CRITERIA \ (residual) \ (solution) \ (residual_mass_wgh) // clang-format on AKANTU_CLASS_ENUM_DECLARE(SolveConvergenceCriteria, AKANTU_SOLVE_CONVERGENCE_CRITERIA) AKANTU_CLASS_ENUM_OUTPUT_STREAM(SolveConvergenceCriteria, AKANTU_SOLVE_CONVERGENCE_CRITERIA) AKANTU_CLASS_ENUM_INPUT_STREAM(SolveConvergenceCriteria, AKANTU_SOLVE_CONVERGENCE_CRITERIA) #else /// enum SolveConvergenceCriteria different convergence criteria enum class SolveConvergenceCriteria { _residual, ///< Use residual to test the convergence _solution, ///< Use solution to test the convergence _residual_mass_wgh ///< Use residual weighted by inv. nodal mass to ///< testb }; #endif /// enum CohesiveMethod type of insertion of cohesive elements enum CohesiveMethod { _intrinsic, _extrinsic }; /// @enum SparseMatrixType type of sparse matrix used enum MatrixType { _unsymmetric, _symmetric, _mt_not_defined }; /* -------------------------------------------------------------------------- */ /* Ghosts handling */ /* -------------------------------------------------------------------------- */ /// @enum CommunicatorType type of communication method to use enum CommunicatorType { _communicator_mpi, _communicator_dummy }; #if !defined(DOXYGEN) // clang-format off #define AKANTU_SYNCHRONIZATION_TAG \ (whatever) \ (update) \ (ask_nodes) \ (size) \ (smm_mass) \ (smm_for_gradu) \ (smm_boundary) \ (smm_uv) \ (smm_res) \ (smm_init_mat) \ (smm_stress) \ (smmc_facets) \ (smmc_facets_conn) \ (smmc_facets_stress) \ (smmc_damage) \ (giu_global_conn) \ (ce_groups) \ (gm_clusters) \ (htm_temperature) \ (htm_gradient_temperature) \ (htm_phi) \ (htm_gradient_phi) \ (pfm_damage) \ - (pfm_gradient_damage) \ + (pfm_driving) \ + (pfm_history) \ + (pfm_energy) \ + (csp_damage) \ + (csp_strain) \ (mnl_for_average) \ (mnl_weight) \ (nh_criterion) \ (test) \ (user_1) \ (user_2) \ (material_id) \ (for_dump) \ (cf_nodal) \ (cf_incr) \ (solver_solution) // clang-format on AKANTU_CLASS_ENUM_DECLARE(SynchronizationTag, AKANTU_SYNCHRONIZATION_TAG) AKANTU_CLASS_ENUM_OUTPUT_STREAM(SynchronizationTag, AKANTU_SYNCHRONIZATION_TAG) #else /// @enum SynchronizationTag type of synchronizations enum class SynchronizationTag { //--- Generic tags --- _whatever, _update, _ask_nodes, _size, //--- SolidMechanicsModel tags --- _smm_mass, ///< synchronization of the SolidMechanicsModel.mass _smm_for_gradu, ///< synchronization of the /// SolidMechanicsModel.displacement _smm_boundary, ///< synchronization of the boundary, forces, velocities /// and displacement _smm_uv, ///< synchronization of the nodal velocities and displacement _smm_res, ///< synchronization of the nodal residual _smm_init_mat, ///< synchronization of the data to initialize materials _smm_stress, ///< synchronization of the stresses to compute the ///< internal /// forces _smmc_facets, ///< synchronization of facet data to setup facet synch _smmc_facets_conn, ///< synchronization of facet global connectivity _smmc_facets_stress, ///< synchronization of facets' stress to setup ///< facet /// synch _smmc_damage, ///< synchronization of damage // --- GlobalIdsUpdater tags --- _giu_global_conn, ///< synchronization of global connectivities // --- CohesiveElementInserter tags --- _ce_groups, ///< synchronization of cohesive element insertion depending /// on facet groups // --- GroupManager tags --- _gm_clusters, ///< synchronization of clusters // --- HeatTransfer tags --- _htm_temperature, ///< synchronization of the nodal temperature _htm_gradient_temperature, ///< synchronization of the element gradient /// temperature // --- PhaseFieldModel tags --- _pfm_damage, ///< synchronization of the nodal damage - _pfm_gradient_damage, ///< synchronization of the element - /// gradient damage - + _pfm_driving, ///< synchronization of the driving forces to + /// compute the internal + _pfm_history, ///< synchronization of the damage history to + /// compute the internal + _pfm_energy, ///< synchronization of the damage energy + /// density to compute the internal + + // --- CouplerSolidPhaseField tags --- + _csp_damage, ///< synchronization of the damage from phase + /// model to solid model + _csp_strain, ///< synchronization of the strain from solid + /// model to phase model + // --- LevelSet tags --- _htm_phi, ///< synchronization of the nodal level set value phi _htm_gradient_phi, ///< synchronization of the element gradient phi //--- Material non local --- _mnl_for_average, ///< synchronization of data to average in non local /// material _mnl_weight, ///< synchronization of data for the weight computations // --- NeighborhoodSynchronization tags --- _nh_criterion, // --- General tags --- _test, ///< Test tag _user_1, ///< tag for user simulations _user_2, ///< tag for user simulations _material_id, ///< synchronization of the material ids _for_dump, ///< everything that needs to be synch before dump // --- Contact & Friction --- _cf_nodal, ///< synchronization of disp, velo, and current position _cf_incr, ///< synchronization of increment // --- Solver tags --- _solver_solution ///< synchronization of the solution obained with the /// PETSc solver }; #endif /// @enum GhostType type of ghost enum GhostType { _not_ghost = 0, _ghost = 1, _casper // not used but a real cute ghost }; /// Define the flag that can be set to a node enum class NodeFlag : std::uint8_t { _normal = 0x00, _distributed = 0x01, _master = 0x03, _slave = 0x05, _pure_ghost = 0x09, _shared_mask = 0x0F, _periodic = 0x10, _periodic_master = 0x30, _periodic_slave = 0x50, _periodic_mask = 0xF0, _local_master_mask = 0xCC, // ~(_master & _periodic_mask) }; inline NodeFlag operator&(const NodeFlag & a, const NodeFlag & b) { using under = std::underlying_type_t; return NodeFlag(under(a) & under(b)); } inline NodeFlag operator|(const NodeFlag & a, const NodeFlag & b) { using under = std::underlying_type_t; return NodeFlag(under(a) | under(b)); } inline NodeFlag & operator|=(NodeFlag & a, const NodeFlag & b) { a = a | b; return a; } inline NodeFlag & operator&=(NodeFlag & a, const NodeFlag & b) { a = a & b; return a; } inline NodeFlag operator~(const NodeFlag & a) { using under = std::underlying_type_t; return NodeFlag(~under(a)); } std::ostream & operator<<(std::ostream & stream, NodeFlag flag); } // namespace akantu AKANTU_ENUM_HASH(GhostType) namespace akantu { /* -------------------------------------------------------------------------- */ struct GhostType_def { using type = GhostType; static const type _begin_ = _not_ghost; static const type _end_ = _casper; }; using ghost_type_t = safe_enum; extern ghost_type_t ghost_types; /// standard output stream operator for GhostType inline std::ostream & operator<<(std::ostream & stream, GhostType type); /* -------------------------------------------------------------------------- */ /* Global defines */ /* -------------------------------------------------------------------------- */ #define AKANTU_MIN_ALLOCATION 2000 #define AKANTU_INDENT ' ' #define AKANTU_INCLUDE_INLINE_IMPL /* -------------------------------------------------------------------------- */ #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_DEREF_PTR(name, ptr) \ inline decltype(auto) get##name() const { \ if (not ptr) { \ AKANTU_EXCEPTION("The member " << #ptr << " is not initialized"); \ } \ return (*ptr); \ } #define AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, support, con) \ inline con Array & get##name( \ const support & el_type, const GhostType & ghost_type = _not_ghost) \ con { \ return variable(el_type, ghost_type); \ } #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, ) #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, const) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, ) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, const) /* -------------------------------------------------------------------------- */ /// initialize the static part of akantu void initialize(int & argc, char **& argv); /// initialize the static part of akantu and read the global input_file void initialize(const std::string & input_file, int & argc, char **& argv); /* -------------------------------------------------------------------------- */ /// finilize correctly akantu and clean the memory void finalize(); /* -------------------------------------------------------------------------- */ /// Read an new input file void readInputFile(const std::string & input_file); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* string manipulation */ /* -------------------------------------------------------------------------- */ inline std::string to_lower(const std::string & str); /* -------------------------------------------------------------------------- */ inline std::string trim(const std::string & to_trim); inline std::string trim(const std::string & to_trim, char c); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /// give a string representation of the a human readable size in bit template std::string printMemorySize(UInt size); /* -------------------------------------------------------------------------- */ struct TensorTrait {}; struct TensorProxyTrait {}; } // namespace akantu /* -------------------------------------------------------------------------- */ /* Type traits */ /* -------------------------------------------------------------------------- */ namespace aka { /* ------------------------------------------------------------------------ */ template using is_tensor = std::is_base_of; template using is_tensor_proxy = std::is_base_of; /* ------------------------------------------------------------------------ */ template using is_scalar = std::is_arithmetic; /* ------------------------------------------------------------------------ */ template ::value> * = nullptr> bool is_of_type(T && t) { return ( dynamic_cast>::value, std::add_const_t, R>>>(&t) != nullptr); } /* -------------------------------------------------------------------------- */ template bool is_of_type(std::unique_ptr & t) { return ( dynamic_cast::value, std::add_const_t, R>>>( t.get()) != nullptr); } /* ------------------------------------------------------------------------ */ template ::value> * = nullptr> decltype(auto) as_type(T && t) { static_assert( disjunction< std::is_base_of, std::decay_t>, // down-cast std::is_base_of, std::decay_t> // up-cast >::value, "Type T and R are not valid for a as_type conversion"); return dynamic_cast>::value, std::add_const_t, R>>>(t); } /* -------------------------------------------------------------------------- */ template ::value> * = nullptr> decltype(auto) as_type(T && t) { return &as_type(*t); } /* -------------------------------------------------------------------------- */ template decltype(auto) as_type(const std::shared_ptr & t) { return std::dynamic_pointer_cast(t); } } // namespace aka #include "aka_fwd.hh" namespace 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(); #define AKANTU_CURRENT_FUNCTION \ (std::string(__func__) + "():" + std::to_string(__LINE__)) } // namespace akantu #include "aka_common_inline_impl.cc" /* -------------------------------------------------------------------------- */ #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 namespace std { /** * 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 struct hash> { hash() = default; size_t operator()(const std::pair & p) const { size_t seed = ah(p.first); return bh(p.second) + AKANTU_HASH_COMBINE_MAGIC_NUMBER + (seed << 6) + (seed >> 2); } private: const hash ah{}; const hash bh{}; }; } // namespace std #endif /* __AKANTU_COMMON_HH__ */ diff --git a/src/model/model_couplers/coupler_solid_phasefield.cc b/src/model/model_couplers/coupler_solid_phasefield.cc index dd67b7036..83e420964 100644 --- a/src/model/model_couplers/coupler_solid_phasefield.cc +++ b/src/model/model_couplers/coupler_solid_phasefield.cc @@ -1,512 +1,590 @@ /** * @file coupler_solid_phasefield.cc * * @author Mohit Pundir * * @date creation: Fri Sep 28 2018 * @date last modification: Thu Jun 20 2019 * * @brief class for coupling of solid mechancis and phasefield model * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "coupler_solid_phasefield.hh" #include "dumpable_inline_impl.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" +#include "element_synchronizer.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { CouplerSolidPhaseField::CouplerSolidPhaseField(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id, const ModelType model_type) : Model(mesh, model_type, dim, id, memory_id) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("CouplerSolidPhaseField", mesh, Model::spatial_dimension); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("coupler_solid_phasefield", id, true); this->mesh.addDumpMeshToDumper("coupler_solid_phasefield", mesh, Model::spatial_dimension, _not_ghost, _ek_regular); #endif this->registerDataAccessor(*this); solid = new SolidMechanicsModel(mesh, Model::spatial_dimension, "solid_mechanics_model", 0); phasefield = new PhaseFieldModel(mesh, Model::spatial_dimension, "phase_field_model", 0); + if (this->mesh.isDistributed()) { + auto & synchronizer = this->mesh.getElementSynchronizer(); + this->registerSynchronizer(synchronizer, SynchronizationTag::_csp_damage); + this->registerSynchronizer(synchronizer, SynchronizationTag::_csp_strain); + + } + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ CouplerSolidPhaseField::~CouplerSolidPhaseField() {} /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); this->initBC(*this, *displacement, *displacement_increment, *external_force); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::initModel() { getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ FEEngine & CouplerSolidPhaseField::getFEEngineBoundary(const ID & name) { return dynamic_cast( getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::initSolver(TimeStepSolverType, NonLinearSolverType) {} /* -------------------------------------------------------------------------- */ std::tuple CouplerSolidPhaseField::getDefaultSolverID(const AnalysisMethod & method) { return std::make_tuple("unkown", TimeStepSolverType::_not_defined); } /* -------------------------------------------------------------------------- */ TimeStepSolverType CouplerSolidPhaseField::getDefaultSolverType() const { return TimeStepSolverType::_dynamic_lumped; } /* -------------------------------------------------------------------------- */ ModelSolverOptions CouplerSolidPhaseField::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_dynamic_lumped: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_dynamic: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_linear; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); break; } return options; } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleResidual() { // computes the internal forces this->assembleInternalForces(); auto & solid_internal_force = solid->getInternalForce(); auto & solid_external_force = solid->getExternalForce(); auto & phasefield_internal_force = phasefield->getInternalForce(); auto & phasefield_external_force = phasefield->getExternalForce(); /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", solid_external_force, 1); this->getDOFManager().assembleToResidual("displacement", solid_internal_force, 1); this->getDOFManager().assembleToResidual("damage", phasefield_external_force, 1); this->getDOFManager().assembleToResidual("damage", phasefield_internal_force, 1); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleResidual(const ID & residual_part) { AKANTU_DEBUG_IN(); auto & solid_internal_force = solid->getInternalForce(); auto & solid_external_force = solid->getExternalForce(); auto & phasefield_internal_force = phasefield->getInternalForce(); auto & phasefield_external_force = phasefield->getExternalForce(); if ("external" == residual_part) { this->getDOFManager().assembleToResidual("displacement", solid_external_force, 1); this->getDOFManager().assembleToResidual("displacement", solid_internal_force, 1); AKANTU_DEBUG_OUT(); return; } if ("internal" == residual_part) { this->getDOFManager().assembleToResidual("damage", phasefield_external_force, 1); this->getDOFManager().assembleToResidual("damage", phasefield_internal_force, 1); AKANTU_DEBUG_OUT(); return; } AKANTU_CUSTOM_EXCEPTION( debug::SolverCallbackResidualPartUnknown(residual_part)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MatrixType CouplerSolidPhaseField::getMatrixType(const ID & matrix_id) { if (matrix_id == "K") return _symmetric; if (matrix_id == "M") { return _symmetric; } return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } else if (matrix_id == "M") { solid->assembleMass(); } } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleLumpedMatrix(const ID & matrix_id) { if (matrix_id == "M") { solid->assembleMassLumped(); } } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); solid->assembleInternalForces(); phasefield->assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); solid->assembleStiffnessMatrix(); phasefield->assembleStiffnessMatrix(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleMassLumped() { solid->assembleMassLumped(); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleMass() { solid->assembleMass(); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleMassLumped(GhostType ghost_type) { //solid->assembleMassLumped(ghost_type); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleMass(GhostType ghost_type) { //solid->assembleMass(ghost_type); } /* ------------------------------------------------------------------------- */ void CouplerSolidPhaseField::computeDamageOnQuadPoints( const GhostType & ghost_type) { AKANTU_DEBUG_IN(); auto & fem = phasefield->getFEEngine(); auto & mesh = phasefield->getMesh(); switch (spatial_dimension) { case 1: { auto & mat = static_cast &>(solid->getMaterial(0)); auto & damage = mat.getDamage(); for (auto & type : mesh.elementTypes(Model::spatial_dimension, ghost_type)) { auto & damage_on_qpoints_vect = damage(type, ghost_type); fem.interpolateOnIntegrationPoints( phasefield->getDamage(), damage_on_qpoints_vect, 1, type, ghost_type); } break; } case 2: { auto & mat = static_cast &>(solid->getMaterial(0)); auto & damage = mat.getDamage(); for (auto & type : mesh.elementTypes(Model::spatial_dimension, ghost_type)) { auto & damage_on_qpoints_vect = damage(type, ghost_type); fem.interpolateOnIntegrationPoints( phasefield->getDamage(), damage_on_qpoints_vect, 1, type, ghost_type); } break; } default: auto & mat = static_cast &>(solid->getMaterial(0)); break; } AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------- */ void CouplerSolidPhaseField::computeStrainOnQuadPoints( const GhostType & ghost_type) { AKANTU_DEBUG_IN(); auto & strain_on_qpoints = phasefield->getStrain(); auto & gradu_on_qpoints = solid->getMaterial(0).getGradU(); for (auto & type : mesh.elementTypes(Model::spatial_dimension, ghost_type)) { auto & strain_on_qpoints_vect = strain_on_qpoints(type, ghost_type); auto & gradu_on_qpoints_vect = gradu_on_qpoints(type, ghost_type); for (auto && values : zip(make_view(strain_on_qpoints_vect, Model::spatial_dimension, Model::spatial_dimension), make_view(gradu_on_qpoints_vect, Model::spatial_dimension, Model::spatial_dimension))) { auto & strain = std::get<0>(values); auto & grad_u = std::get<1>(values); this->gradUToEpsilon(grad_u, strain); } } AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------- */ void CouplerSolidPhaseField::solve() { solid->solveStep(); this->computeStrainOnQuadPoints(_not_ghost); phasefield->solveStep(); this->computeDamageOnQuadPoints(_not_ghost); solid->assembleInternalForces(); } /* ------------------------------------------------------------------------- */ void CouplerSolidPhaseField::gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon) { for (UInt i = 0; i < Model::spatial_dimension; ++i) { for (UInt j = 0; j < Model::spatial_dimension; ++j) epsilon(i, j) = 0.5 * (grad_u(i, j) + grad_u(j, i)); } } /* ------------------------------------------------------------------------- */ bool CouplerSolidPhaseField::checkConvergence(Array & u_new, Array & u_old, Array & d_new, Array & d_old) { const Array & blocked_dofs = solid->getBlockedDOFs(); UInt nb_degree_of_freedom = u_new.size(); auto u_n_it = u_new.begin(); auto u_o_it = u_old.begin(); auto bld_it = blocked_dofs.begin(); Real norm = 0; for (UInt n = 0; n < nb_degree_of_freedom; ++n, ++u_n_it, ++u_o_it, ++bld_it) { if ((!*bld_it)) { norm += (*u_n_it - *u_o_it) * (*u_n_it - *u_o_it); } } norm = std::sqrt(norm); auto d_n_it = d_new.begin(); auto d_o_it = d_old.begin(); nb_degree_of_freedom = d_new.size(); Real norm2 = 0; for (UInt i = 0; i < nb_degree_of_freedom; ++i) { norm2 += (*d_n_it - *d_o_it); } norm2 = std::sqrt(norm2); Real error = std::max(norm, norm2); Real tolerance = 1e-8; if (error < tolerance) { return true; } return false; } + +/* -------------------------------------------------------------------------- */ +UInt CouplerSolidPhaseField::getNbData(const Array & elements, + const SynchronizationTag & tag) const { + + AKANTU_DEBUG_IN(); + + UInt size = 0; + + switch (tag) { + case SynchronizationTag::_csp_damage: { + size += getNbIntegrationPoints(elements) * sizeof(Real); + break; + } + case SynchronizationTag::_csp_strain: { + size += getNbIntegrationPoints(elements) * spatial_dimension * spatial_dimension * sizeof(Real); + break; + } + default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } + } + + AKANTU_DEBUG_OUT(); + return size; +} + + +/* -------------------------------------------------------------------------- */ +void CouplerSolidPhaseField::packData(CommunicationBuffer & buffer, + const Array & elements, + const SynchronizationTag & tag) const { + + auto & mat = static_cast &>(solid->getMaterial(0)); + + switch (tag) { + case SynchronizationTag::_csp_damage: { + packElementalDataHelper(mat.getDamage(), buffer, elements, + true, getFEEngine()); + break; + } + case SynchronizationTag::_csp_strain: { + packElementalDataHelper(phasefield->getStrain(), buffer, elements, + true, getFEEngine()); + break; + } + default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } + } +} + +/* -------------------------------------------------------------------------- */ +void CouplerSolidPhaseField::unpackData(CommunicationBuffer & buffer, + const Array & elements, + const SynchronizationTag & tag) { + + auto & mat = static_cast &>(solid->getMaterial(0)); + + switch (tag) { + case SynchronizationTag::_csp_damage: { + unpackElementalDataHelper(mat.getDamage(), buffer, elements, true, getFEEngine()); + break; + } + case SynchronizationTag::_csp_strain: { + unpackElementalDataHelper(phasefield->getStrain(), buffer, elements, + true, getFEEngine()); + break; + } + default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } + } +} + + /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidPhaseField::createElementalField( const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) { return solid->createElementalField(field_name, group_name, padding_flag, spatial_dimension, kind); std::shared_ptr field; return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidPhaseField::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { return solid->createNodalFieldReal(field_name, group_name, padding_flag); std::shared_ptr field; return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidPhaseField::createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) { return solid->createNodalFieldBool(field_name, group_name, padding_flag); std::shared_ptr field; return field; } #else /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidPhaseField::createElementalField( const std::string &, const std::string &, bool, const UInt &, const ElementKind &) { return nullptr; } /* ----------------------------------------------------------------------- */ std::shared_ptr CouplerSolidPhaseField::createNodalFieldReal(const std::string &, const std::string &, bool) { return nullptr; } /*-------------------------------------------------------------------*/ std::shared_ptr CouplerSolidPhaseField::createNodalFieldBool(const std::string &, const std::string &, bool) { return nullptr; } #endif /* -----------------------------------------------------------------------*/ void CouplerSolidPhaseField::dump(const std::string & dumper_name) { solid->onDump(); mesh.dump(dumper_name); } /* ------------------------------------------------------------------------*/ void CouplerSolidPhaseField::dump(const std::string & dumper_name, UInt step) { solid->onDump(); mesh.dump(dumper_name, step); } /* ----------------------------------------------------------------------- */ void CouplerSolidPhaseField::dump(const std::string & dumper_name, Real time, UInt step) { solid->onDump(); mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::dump() { solid->onDump(); mesh.dump(); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::dump(UInt step) { solid->onDump(); mesh.dump(step); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::dump(Real time, UInt step) { solid->onDump(); mesh.dump(time, step); } } // namespace akantu diff --git a/src/model/model_couplers/coupler_solid_phasefield.hh b/src/model/model_couplers/coupler_solid_phasefield.hh index a6f2a90e8..4deaba1cb 100644 --- a/src/model/model_couplers/coupler_solid_phasefield.hh +++ b/src/model/model_couplers/coupler_solid_phasefield.hh @@ -1,275 +1,273 @@ /** * @file solid_phase_coupler.hh * * @author Mohit Pundir * * @date creation: Fri Sep 28 2018 * @date last modification: Fri Sep 28 2018 * * @brief class for coupling of solid mechancis and phasefield model * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "boundary_condition.hh" #include "data_accessor.hh" #include "fe_engine.hh" #include "material.hh" #include "material_phasefield.hh" #include "model.hh" #include "phase_field_model.hh" #include "solid_mechanics_model.hh" #include "sparse_matrix.hh" #include "time_step_solver.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COUPLER_SOLID_PHASEFIELD_HH__ #define __AKANTU_COUPLER_SOLID_PHASEFIELD_HH__ /* ------------------------------------------------------------------------ */ /* Coupling : Solid Mechanics / PhaseField */ /* ------------------------------------------------------------------------ */ namespace akantu { template class IntegratorGauss; template class ShapeLagrange; class DOFManager; } // namespace akantu namespace akantu { class CouplerSolidPhaseField : public Model, public DataAccessor, public DataAccessor, public BoundaryCondition { /* ------------------------------------------------------------------------ */ /* Constructor/Destructors */ /* ------------------------------------------------------------------------ */ using MyFEEngineType = FEEngineTemplate; public: CouplerSolidPhaseField( Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "coupler_solid_phasefield", const MemoryID & memory_id = 0, const ModelType model_type = ModelType::_coupler_solid_phasefield); ~CouplerSolidPhaseField() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize the complete model void initFullImpl(const ModelOptions & options) override; /// initialize the modelType void initModel() override; /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; /* ------------------------------------------------------------------------ */ /* Solver Interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the contact stiffness matrix virtual void assembleStiffnessMatrix(); /// assembles the contant internal forces virtual void assembleInternalForces(); public: /// computes damage on quad points for solid mechanics model from /// damage array from phasefield model void computeDamageOnQuadPoints(const GhostType &); /// computes strain on quadrature points for phasefield model from /// displacement gradient from solid mechanics model void computeStrainOnQuadPoints(const GhostType & ghost_type); /// solve the coupled model void solve(); private: /// computes small strain from displacement gradient void gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon); /// test the convergence criteria bool checkConvergence(Array &, Array &, Array &, Array &); protected: /// callback for the solver, this adds f_{ext} - f_{int} to the residual void assembleResidual() override; /// callback for the solver, this adds f_{ext} or f_{int} to the residual void assembleResidual(const ID & residual_part) override; bool canSplitResidual() override { return true; } /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) override; /// callback for the solver, this assembles different matrices void assembleMatrix(const ID & matrix_id) override; /// callback for the solver, this assembles the stiffness matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// callback for the model to instantiate the matricess when needed void initSolver(TimeStepSolverType, NonLinearSolverType) override; /// solve the coupled model //void solveStep(const ID & solver_id = "") override; /// solve a step using a given pre instantiated time step solver and /// non linear solver with a user defined callback instead of the /// model itself /!\ This can mess up everything //void solveStep(SolverCallback & callback, const ID & solver_id = "") override; /* ------------------------------------------------------------------------ */ /* Mass matrix for solid mechanics model */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); protected: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMass(GhostType ghost_type); protected: /* -------------------------------------------------------------------------- */ TimeStepSolverType getDefaultSolverType() const override; /* -------------------------------------------------------------------------- */ ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const; public: bool isDefaultSolverExplicit() { return method == _explicit_lumped_mass; } /* ------------------------------------------------------------------------ */ public: // DataAccessor UInt getNbData(const Array &, - const SynchronizationTag &) const override { - return 0; - } + const SynchronizationTag &) const override; void packData(CommunicationBuffer &, const Array &, - const SynchronizationTag &) const override {} + const SynchronizationTag &) const override; void unpackData(CommunicationBuffer &, const Array &, - const SynchronizationTag &) override {} - - // DataAccessor nodes - UInt getNbData(const Array &, - const SynchronizationTag &) const override { - return 0; - } - void packData(CommunicationBuffer &, const Array &, - const SynchronizationTag &) const override {} - void unpackData(CommunicationBuffer &, const Array &, - const SynchronizationTag &) override {} + const SynchronizationTag &) override; + + UInt getNbData(const Array & indexes, + const SynchronizationTag & tag) const override {} + + void packData(CommunicationBuffer & buffer, const Array & dofs, + const SynchronizationTag & tag) const override{} + + void unpackData(CommunicationBuffer & buffer, const Array & dofs, + const SynchronizationTag & tag) override {} + /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: FEEngine & getFEEngineBoundary(const ID & name = "") override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt); /// get the solid mechanics model AKANTU_GET_MACRO(SolidMechanicsModel, *solid, SolidMechanicsModel &); /// get the contact mechanics model AKANTU_GET_MACRO(PhaseFieldModel, *phasefield, PhaseFieldModel &); /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: std::shared_ptr createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) override; virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); void dump() override; virtual void dump(UInt step); virtual void dump(Real time, UInt step); /* ------------------------------------------------------------------------ */ /* Members */ /* ------------------------------------------------------------------------ */ private: /// solid mechanics model SolidMechanicsModel * solid{nullptr}; /// phasefield model PhaseFieldModel * phasefield{nullptr}; Array * displacement{nullptr}; /// Array * displacement_increment{nullptr}; /// external forces array Array * external_force{nullptr}; }; } // namespace akantu #endif /* __AKANTU_COUPLER_SOLID_PHASEFIELD_HH__ */ diff --git a/src/model/phase_field/phase_field_model.cc b/src/model/phase_field/phase_field_model.cc index 40f8fa4cb..54e9de237 100644 --- a/src/model/phase_field/phase_field_model.cc +++ b/src/model/phase_field/phase_field_model.cc @@ -1,855 +1,920 @@ /** * @file phase_field_model.cc * * @author Mohit Pundir * * @date creation: Wed Aug 01 2018 * @date last modification: Wed Aug 01 2018 * * @brief Implementation of PhaseFieldModel class * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "phase_field_model.hh" #include "dumpable_inline_impl.hh" #include "element_synchronizer.hh" #include "fe_engine_template.hh" #include "generalized_trapezoidal.hh" #include "group_manager_inline_impl.cc" #include "integrator_gauss.hh" #include "mesh.hh" #include "parser.hh" #include "shape_lagrange.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_element_partition.hh" #include "dumper_elemental_field.hh" #include "dumper_internal_material_field.hh" #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { namespace phasefield { namespace details { class ComputeDamageFunctor { public: ComputeDamageFunctor(const PhaseFieldModel & model) : model(model){}; // material functor can be created void operator()(Matrix & rho, const Element &) const { rho.set(1.); } private: const PhaseFieldModel & model; }; } // namespace details } // namespace phasefield /* -------------------------------------------------------------------------- */ PhaseFieldModel::PhaseFieldModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id, const ModelType model_type) : Model(mesh, model_type, dim, id, memory_id), BoundaryCondition(), damage_on_qpoints("damage_on_qpoints", id), damage_energy_on_qpoints("damage_energy_on_qpoints", id), damage_energy_density_on_qpoints("damage_energy_density_on_qpoints", id), damage_gradient("damage_gradient", id), strain_on_qpoints("strain_on_qpoints", id), driving_force_on_qpoints("driving_force_on_qpoints", id), phi_history_on_qpoints("phi_history_on_qpoints", id) { AKANTU_DEBUG_IN(); - - + this->initDOFManager(); this->registerDataAccessor(*this); if (this->mesh.isDistributed()) { auto & synchronizer = this->mesh.getElementSynchronizer(); this->registerSynchronizer(synchronizer, SynchronizationTag::_pfm_damage); - this->registerSynchronizer(synchronizer, - SynchronizationTag::_pfm_gradient_damage); + this->registerSynchronizer(synchronizer, SynchronizationTag::_pfm_driving); + this->registerSynchronizer(synchronizer, SynchronizationTag::_pfm_history); + this->registerSynchronizer(synchronizer, SynchronizationTag::_pfm_energy); } this->registerFEEngineObject("PhaseFieldFEEngine", mesh, Model::spatial_dimension); #ifdef AKANTU_USE_IOHELPER this->mesh.registerDumper("phase_field", id, true); this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular); #endif // AKANTU_USE_IOHELPER this->registerParam("l0", l_0, 0., _pat_parsmod, "length scale"); this->registerParam("gc", g_c, _pat_parsmod, "critical local fracture energy density"); this->registerParam("E", E, _pat_parsmod, "Young's modulus"); this->registerParam("nu", nu, _pat_parsmod, "Poisson ratio"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ PhaseFieldModel::~PhaseFieldModel() = default; /* -------------------------------------------------------------------------- */ MatrixType PhaseFieldModel::getMatrixType(const ID & matrix_id) { if (matrix_id == "K" or matrix_id == "M") { return _symmetric; } return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::initModel() { auto & fem = this->getFEEngine(); fem.initShapeFunctions(_not_ghost); fem.initShapeFunctions(_ghost); damage_on_qpoints.initialize(fem, _nb_component = 1); damage_energy_on_qpoints.initialize(fem, _nb_component = spatial_dimension * spatial_dimension); damage_energy_density_on_qpoints.initialize(fem, _nb_component = 1); damage_gradient.initialize(fem, _nb_component = spatial_dimension); strain_on_qpoints.initialize(fem, _nb_component = spatial_dimension * spatial_dimension); driving_force_on_qpoints.initialize(fem, _nb_component = 1); phi_history_on_qpoints.initialize(fem, _nb_component = 1); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); readMaterials(); this->initBC(*this, *damage, *external_force); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::readMaterials() { auto sect = this->getParserSection(); if (not std::get<1>(sect)) { const auto & section = std::get<0>(sect); this->parseSection(section); } Matrix d(spatial_dimension, spatial_dimension); d.eye(g_c * l_0); damage_energy_on_qpoints.set(d); this->updateInternalParameters(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::updateInternalParameters() { this->lambda = this->nu * this->E / ((1 + this->nu) * (1 - 2 * this->nu)); this->mu = this->E / (2 * (1 + this->nu)); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } else { AKANTU_ERROR("Unknown Matrix ID for PhaseFieldModel : " << matrix_id); } } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::predictor() { // AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::corrector() { // AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType) { DOFManager & dof_manager = this->getDOFManager(); this->allocNodalField(this->damage, 1, "damage"); this->allocNodalField(this->external_force, 1, "external_force"); this->allocNodalField(this->internal_force, 1, "internal_force"); this->allocNodalField(this->blocked_dofs, 1, "blocked_dofs"); this->allocNodalField(this->previous_damage, 1, "previous_damage"); this->allocNodalField(this->damage_increment, 1, "damage_increment"); if (!dof_manager.hasDOFs("damage")) { dof_manager.registerDOFs("damage", *this->damage, _dst_nodal); dof_manager.registerBlockedDOFs("damage", *this->blocked_dofs); dof_manager.registerDOFsIncrement("damage", *this->damage_increment); dof_manager.registerDOFsPrevious("damage", *this->previous_damage); } if (time_step_solver_type == TimeStepSolverType::_dynamic) { AKANTU_TO_IMPLEMENT(); } } /* -------------------------------------------------------------------------- */ FEEngine & PhaseFieldModel::getFEEngineBoundary(const ID & name) { return dynamic_cast(getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ std::tuple PhaseFieldModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_lumped_mass: { return std::make_tuple("explicit_lumped", TimeStepSolverType::_dynamic_lumped); } case _static: { return std::make_tuple("static", TimeStepSolverType::_static); } case _implicit_dynamic: { return std::make_tuple("implicit", TimeStepSolverType::_dynamic); } default: return std::make_tuple("unknown", TimeStepSolverType::_not_defined); } } /* -------------------------------------------------------------------------- */ ModelSolverOptions PhaseFieldModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_linear; options.integration_scheme_type["damage"] = IntegrationSchemeType::_pseudo_time; options.solution_type["damage"] = IntegrationScheme::_not_defined; break; } case TimeStepSolverType::_dynamic: { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["damage"] = IntegrationSchemeType::_backward_euler; options.solution_type["damage"] = IntegrationScheme::_damage; break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } return options; } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::beforeSolveStep() { + // compute the history of local elements + AKANTU_DEBUG_INFO("Compute phi history"); this->computePhiHistoryOnQuadPoints(_not_ghost); + + // communicate the history + AKANTU_DEBUG_INFO("Send data for synchronization"); + this->asynchronousSynchronize(SynchronizationTag::_pfm_history); + + // finalize communications + AKANTU_DEBUG_INFO("Wait distant history"); + this->waitEndSynchronize(SynchronizationTag::_pfm_history); + this->computeDamageEnergyDensityOnQuadPoints(_not_ghost); + + // communicate the energy density + AKANTU_DEBUG_INFO("Send data for synchronization"); + this->asynchronousSynchronize(SynchronizationTag::_pfm_energy); + + // finalize communications + AKANTU_DEBUG_INFO("Wait distant energy density"); + this->waitEndSynchronize(SynchronizationTag::_pfm_energy); + } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::afterSolveStep(bool converged) { if (not converged) return ; for (auto && values : zip(*damage, *previous_damage)) { auto & dam = std::get<0>(values); auto & prev_dam = std::get<1>(values); dam -= prev_dam; dam = std::min(1., 2 * dam - dam * dam); prev_dam = dam; } } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); if (!this->getDOFManager().hasMatrix("K")) { this->getDOFManager().getNewMatrix("K", getMatrixType("K")); } this->getDOFManager().clearMatrix("K"); this->assembleDamageMatrix(); this->assembleDamageGradMatrix(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::assembleDamageMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new damage matrix"); auto & fem = this->getFEEngine(); for (auto && type : mesh.elementTypes(spatial_dimension)) { auto nb_element = mesh.getNbElement(type); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = fem.getNbIntegrationPoints(type); auto nt_b_n = std::make_unique>( nb_element * nb_quadrature_points, nb_nodes_per_element * nb_nodes_per_element, "N^t*b*N"); auto & damage_energy_density_on_qpoints_vect = damage_energy_density_on_qpoints(type); // damage_energy_on_qpoints = gc/l0 + phi = scalar fem.computeNtbN(damage_energy_density_on_qpoints_vect, *nt_b_n, 2, type); /// compute @f$ K_{\grad d} = \int_e \mathbf{N}^t * \mathbf{w} * /// \mathbf{N}@f$ auto K_n = std::make_unique>( nb_element, nb_nodes_per_element * nb_nodes_per_element, "K_n"); fem.integrate(*nt_b_n, *K_n, nb_nodes_per_element * nb_nodes_per_element, type); this->getDOFManager().assembleElementalMatricesToMatrix( "K", "damage", *K_n, type, _not_ghost, _symmetric); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::assembleDamageGradMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new damage gradient matrix"); auto & fem = this->getFEEngine(); for (auto && type : mesh.elementTypes(spatial_dimension)) { auto nb_element = mesh.getNbElement(type); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = fem.getNbIntegrationPoints(type); auto bt_d_b = std::make_unique>( nb_element * nb_quadrature_points, nb_nodes_per_element * nb_nodes_per_element, "B^t*D*B"); auto & damage_energy_on_qpoints_vect = damage_energy_on_qpoints(type); // damage_energy_on_qpoints = gc*l0 = scalar fem.computeBtDB(damage_energy_on_qpoints_vect, *bt_d_b, 2, type); /// compute @f$ K_{\grad d} = \int_e \mathbf{B}^t * \mathbf{W} * /// \mathbf{B}@f$ auto K_b = std::make_unique>( nb_element, nb_nodes_per_element * nb_nodes_per_element, "K_b"); fem.integrate(*bt_d_b, *K_b, nb_nodes_per_element * nb_nodes_per_element, type); this->getDOFManager().assembleElementalMatricesToMatrix( "K", "damage", *K_b, type, _not_ghost, _symmetric); } AKANTU_DEBUG_OUT(); } +/* -------------------------------------------------------------------------- */ +void PhaseFieldModel::assembleResidual() { + + AKANTU_DEBUG_IN(); + + this->assembleInternalForces(); + + this->getDOFManager().assembleToResidual("damage", *this->external_force, 1); + this->getDOFManager().assembleToResidual("damage", *this->internal_force, 1); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void PhaseFieldModel::assembleInternalForces() { + AKANTU_DEBUG_IN(); + + AKANTU_DEBUG_INFO("Assemble the internal forces"); + + this->internal_force->clear(); + + // compute the driving force of local elements + AKANTU_DEBUG_INFO("Compute local driving forces"); + this->computeDrivingForce(_not_ghost); + + // communicate the driving forces + AKANTU_DEBUG_INFO("Send data for residual assembly"); + this->asynchronousSynchronize(SynchronizationTag::_pfm_driving); + + // assemble the forces due to local driving forces + AKANTU_DEBUG_INFO("Assemble residual for local elements"); + this->assembleInternalForces(_not_ghost); + + // finalize communications + AKANTU_DEBUG_INFO("Wait distant driving forces"); + this->waitEndSynchronize(SynchronizationTag::_pfm_driving); + + // assemble the residual due to ghost elements + AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); + this->assembleInternalForces(_ghost); + + AKANTU_DEBUG_OUT(); +} + + +/* -------------------------------------------------------------------------- */ +void PhaseFieldModel::assembleInternalForces(const GhostType & ghost_type) { + AKANTU_DEBUG_IN(); + + this->synchronize(SynchronizationTag::_pfm_damage); + auto & fem = this->getFEEngine(); + + for (auto type : mesh.elementTypes(spatial_dimension, ghost_type)) { + UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); + + auto & driving_force_on_qpoints_vect = + driving_force_on_qpoints(type, ghost_type); + + UInt nb_quad_points = driving_force_on_qpoints_vect.size(); + Array nt_driving_force(nb_quad_points, nb_nodes_per_element); + fem.computeNtb(driving_force_on_qpoints_vect, nt_driving_force, type, + ghost_type); + + UInt nb_elements = mesh.getNbElement(type, ghost_type); + Array int_nt_driving_force(nb_elements, nb_nodes_per_element); + + fem.integrate(nt_driving_force, int_nt_driving_force, + nb_nodes_per_element, type, ghost_type); + + this->getDOFManager().assembleElementalArrayLocalArray( + int_nt_driving_force, *this->internal_force, type, ghost_type, 1); + } + AKANTU_DEBUG_OUT(); +} + + + +/* -------------------------------------------------------------------------- */ + void PhaseFieldModel::assembleLumpedMatrix(const ID & /*matrix_id*/) {} + +/* -------------------------------------------------------------------------- */ +void PhaseFieldModel::setTimeStep(Real time_step, const ID & solver_id) { + Model::setTimeStep(time_step, solver_id); + +#if defined(AKANTU_USE_IOHELPER) + this->mesh.getDumper("phase_field").setTimeStep(time_step); +#endif +} +/* -------------------------------------------------------------------------- */ +void PhaseFieldModel::computeDrivingForce(const GhostType & ghost_type) { + + for (auto & type : mesh.elementTypes(spatial_dimension, ghost_type)) { + + for (auto && values : + zip(make_view(phi_history_on_qpoints(type, ghost_type)), + make_view(driving_force_on_qpoints(type, ghost_type)))) { + auto & phi_history = std::get<0>(values); + auto & driving_force = std::get<1>(values); + + driving_force = 2.0 * phi_history; + } + } + + AKANTU_DEBUG_OUT(); +} + + /* -------------------------------------------------------------------------- */ void PhaseFieldModel::computeDamageEnergyDensityOnQuadPoints( const GhostType & ghost_type) { AKANTU_DEBUG_IN(); for (auto & type : mesh.elementTypes(spatial_dimension, ghost_type)) { for (auto && values : zip(make_view(damage_energy_density_on_qpoints(type, ghost_type)), make_view(phi_history_on_qpoints(type, ghost_type)))) { auto & dam_energy_density = std::get<0>(values); auto & phi_history = std::get<1>(values); dam_energy_density = g_c / l_0 + 2.0 * phi_history; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::computePhiHistoryOnQuadPoints( const GhostType & ghost_type) { Matrix strain_plus(spatial_dimension, spatial_dimension); Matrix strain_minus(spatial_dimension, spatial_dimension); Matrix strain_dir(spatial_dimension, spatial_dimension); Matrix strain_diag_plus(spatial_dimension, spatial_dimension); Matrix strain_diag_minus(spatial_dimension, spatial_dimension); Vector strain_values(spatial_dimension); Real trace_plus, trace_minus, phi_plus; for (auto & type : mesh.elementTypes(spatial_dimension, ghost_type)) { for (auto && values : zip(make_view(strain_on_qpoints(type, ghost_type), spatial_dimension, spatial_dimension), make_view(phi_history_on_qpoints(type, ghost_type)))) { auto & strain = std::get<0>(values); auto & phi_history = std::get<1>(values); strain_plus.clear(); strain_minus.clear(); strain_dir.clear(); strain_values.clear(); strain_diag_plus.clear(); strain_diag_minus.clear(); strain.eig(strain_values, strain_dir); for (UInt i = 0; i < spatial_dimension; i++) { strain_diag_plus(i, i) = std::max(Real(0.), strain_values(i)); strain_diag_minus(i, i) = std::min(Real(0.), strain_values(i)); } Matrix mat_tmp(spatial_dimension, spatial_dimension); Matrix sigma_plus(spatial_dimension, spatial_dimension); Matrix sigma_minus(spatial_dimension, spatial_dimension); mat_tmp.mul(strain_diag_plus, strain_dir); strain_plus.mul(strain_dir, mat_tmp); mat_tmp.mul(strain_diag_minus, strain_dir); strain_minus.mul(strain_dir, mat_tmp); trace_plus = std::max(Real(0.), strain.trace()); trace_minus = std::min(Real(0.), strain.trace()); for (UInt i = 0; i < spatial_dimension; i++) { for (UInt j = 0; j < spatial_dimension; j++) { sigma_plus(i, j) = (i == j) * lambda * trace_plus + 2 * mu * strain_plus(i, j); sigma_minus(i, j) = (i == j) * lambda * trace_minus + 2 * mu * strain_minus(i, j); } } phi_plus = 1. / 2 * sigma_plus.doubleDot(strain); if (phi_plus > phi_history) { phi_history = phi_plus; } } } } - -/* -------------------------------------------------------------------------- */ -void PhaseFieldModel::assembleResidual() { - - AKANTU_DEBUG_IN(); - - this->assembleInternalForces(); - - this->getDOFManager().assembleToResidual("damage", *this->external_force, 1); - this->getDOFManager().assembleToResidual("damage", *this->internal_force, 1); - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -void PhaseFieldModel::assembleInternalForces() { - AKANTU_DEBUG_IN(); - - this->internal_force->clear(); - - this->synchronize(SynchronizationTag::_pfm_damage); - auto & fem = this->getFEEngine(); - - for (auto ghost_type : ghost_types) { - - computeDrivingForce(ghost_type); - - for (auto type : mesh.elementTypes(spatial_dimension, ghost_type)) { - UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); - - auto & driving_force_on_qpoints_vect = - driving_force_on_qpoints(type, ghost_type); - - UInt nb_quad_points = driving_force_on_qpoints_vect.size(); - Array nt_driving_force(nb_quad_points, nb_nodes_per_element); - fem.computeNtb(driving_force_on_qpoints_vect, nt_driving_force, type, - ghost_type); - - UInt nb_elements = mesh.getNbElement(type, ghost_type); - Array int_nt_driving_force(nb_elements, nb_nodes_per_element); - - fem.integrate(nt_driving_force, int_nt_driving_force, - nb_nodes_per_element, type, ghost_type); - - this->getDOFManager().assembleElementalArrayLocalArray( - int_nt_driving_force, *this->internal_force, type, ghost_type, 1); - } - } - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ - void PhaseFieldModel::assembleLumpedMatrix(const ID & /*matrix_id*/) {} - -/* -------------------------------------------------------------------------- */ -void PhaseFieldModel::setTimeStep(Real time_step, const ID & solver_id) { - Model::setTimeStep(time_step, solver_id); - -#if defined(AKANTU_USE_IOHELPER) - this->mesh.getDumper("phase_field").setTimeStep(time_step); -#endif -} -/* -------------------------------------------------------------------------- */ -void PhaseFieldModel::computeDrivingForce(const GhostType & ghost_type) { - - for (auto & type : mesh.elementTypes(spatial_dimension, ghost_type)) { - - for (auto && values : - zip(make_view(phi_history_on_qpoints(type, ghost_type)), - make_view(driving_force_on_qpoints(type, ghost_type)))) { - auto & phi_history = std::get<0>(values); - auto & driving_force = std::get<1>(values); - - driving_force = 2.0 * phi_history; - } - } - - AKANTU_DEBUG_OUT(); -} - + /* -------------------------------------------------------------------------- */ void PhaseFieldModel::computeDamageOnQuadPoints(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); auto & fem = this->getFEEngine(); for (auto & type : mesh.elementTypes(spatial_dimension, ghost_type)) { auto & damage_on_qpoints_vect = damage_on_qpoints(type, ghost_type); fem.interpolateOnIntegrationPoints(*damage, damage_on_qpoints_vect, 1, type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt PhaseFieldModel::getNbData(const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes_per_element = 0; for (const Element & el : elements) { nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); } switch (tag) { case SynchronizationTag::_pfm_damage: { size += nb_nodes_per_element * sizeof(Real); // damage break; } - case SynchronizationTag::_pfm_gradient_damage: { - size += getNbIntegrationPoints(elements) * spatial_dimension * sizeof(Real); - size += nb_nodes_per_element * sizeof(Real); + case SynchronizationTag::_pfm_driving: { + size += getNbIntegrationPoints(elements) * sizeof(Real); + break; + } + case SynchronizationTag::_pfm_history: { + size += getNbIntegrationPoints(elements) * sizeof(Real); break; } + case SynchronizationTag::_pfm_energy: { + size += getNbIntegrationPoints(elements) * sizeof(Real); + break; + } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { switch (tag) { case SynchronizationTag::_pfm_damage: { packNodalDataHelper(*damage, buffer, elements, mesh); break; } - case SynchronizationTag::_pfm_gradient_damage: { - packElementalDataHelper(damage_gradient, buffer, elements, true, - getFEEngine()); - packNodalDataHelper(*damage, buffer, elements, mesh); + case SynchronizationTag::_pfm_driving: { + packElementalDataHelper(driving_force_on_qpoints, buffer, elements, true, getFEEngine()); + break; + } + case SynchronizationTag::_pfm_history: { + packElementalDataHelper(phi_history_on_qpoints, buffer, elements, true, getFEEngine()); + break; + } + case SynchronizationTag::_pfm_energy: { + packElementalDataHelper(damage_energy_density_on_qpoints, buffer, elements, true, getFEEngine()); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { switch (tag) { case SynchronizationTag::_pfm_damage: { unpackNodalDataHelper(*damage, buffer, elements, mesh); break; } - case SynchronizationTag::_pfm_gradient_damage: { - unpackElementalDataHelper(damage_gradient, buffer, elements, true, - getFEEngine()); - unpackNodalDataHelper(*damage, buffer, elements, mesh); + case SynchronizationTag::_pfm_driving: { + unpackElementalDataHelper(driving_force_on_qpoints, buffer, elements, true, getFEEngine()); + break; + } + case SynchronizationTag::_pfm_history: { + unpackElementalDataHelper(phi_history_on_qpoints, buffer, elements, true, getFEEngine()); + break; + } + case SynchronizationTag::_pfm_energy: { + unpackElementalDataHelper(damage_energy_density_on_qpoints, buffer, elements, true, getFEEngine()); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } /* -------------------------------------------------------------------------- */ UInt PhaseFieldModel::getNbData(const Array & indexes, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes = indexes.size(); switch (tag) { case SynchronizationTag::_pfm_damage: { size += nb_nodes * sizeof(Real); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::packData(CommunicationBuffer & buffer, const Array & indexes, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); for (auto index : indexes) { switch (tag) { case SynchronizationTag::_pfm_damage: { buffer << (*damage)(index); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::unpackData(CommunicationBuffer & buffer, const Array & indexes, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); for (auto index : indexes) { switch (tag) { case SynchronizationTag::_pfm_damage: { buffer >> (*damage)(index); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } AKANTU_DEBUG_OUT(); } -/* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER std::shared_ptr PhaseFieldModel::createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool /*padding_flag*/) { std::map *> uint_nodal_fields; uint_nodal_fields["blocked_dofs"] = blocked_dofs; return mesh.createNodalField(uint_nodal_fields[field_name], group_name); std::shared_ptr field; return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr PhaseFieldModel::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool /*padding_flag*/) { if (field_name == "capacity_lumped") { AKANTU_EXCEPTION( "Capacity lumped is a nodal field now stored in the DOF manager." "Therefore it cannot be used by a dumper anymore"); } std::map *> real_nodal_fields; real_nodal_fields["damage"] = damage; real_nodal_fields["external_force"] = external_force; real_nodal_fields["internal_force"] = internal_force; return mesh.createNodalField(real_nodal_fields[field_name], group_name); std::shared_ptr field; return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr PhaseFieldModel::createElementalField( const std::string & field_name, const std::string & group_name, bool /*padding_flag*/, const UInt & /*spatial_dimension*/, const ElementKind & element_kind) { if (field_name == "partitions") { return mesh.createElementalField( mesh.getConnectivities(), group_name, this->spatial_dimension, element_kind); } else if (field_name == "damage_gradient") { ElementTypeMap nb_data_per_elem = this->mesh.getNbDataPerElem(damage_gradient); return mesh.createElementalField( damage_gradient, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); } else if (field_name == "damage_energy") { ElementTypeMap nb_data_per_elem = this->mesh.getNbDataPerElem(damage_energy_on_qpoints); return mesh.createElementalField( damage_energy_on_qpoints, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); } std::shared_ptr field; return field; } /* -------------------------------------------------------------------------- */ #else /* -------------------------------------------------------------------------- */ std::shared_ptr PhaseFieldModel::createElementalField( const std::string & /*field_name*/, const std::string & /*group_name*/, bool /*padding_flag*/, const ElementKind & /*element_kind*/) { return nullptr; } /* -------------------------------------------------------------------------- */ std::shared_ptr PhaseFieldModel::createNodalFieldBool(const std::string & /*field_name*/, const std::string & /*group_name*/, bool /*padding_flag*/) { return nullptr; } /* -------------------------------------------------------------------------- */ std::shared_ptr PhaseFieldModel::createNodalFieldReal(const std::string & /*field_name*/, const std::string & /*group_name*/, bool /*padding_flag*/) { return nullptr; } #endif /* -------------------------------------------------------------------------- */ void PhaseFieldModel::dump(const std::string & dumper_name) { mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::dump(const std::string & dumper_name, UInt step) { mesh.dump(dumper_name, step); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::dump(const std::string & dumper_name, Real time, UInt step) { mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::dump() { mesh.dump(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::dump(UInt step) { mesh.dump(step); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::dump(Real time, UInt step) { mesh.dump(time, step); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "Phase Field Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << Model::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; damage->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; stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + materials [" << std::endl; stream << space << " length scale parameter : " << l_0 << std::endl; stream << space << " critical energy release rate : " << g_c << std::endl; stream << space << " Young's modulus : " << E << std::endl; stream << space << " Poisson's ratio : " << nu << std::endl; stream << space << " Lame's first parameter : " << lambda << std::endl; stream << space << " Lame's second parameter : " << mu << std::endl; stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } } // namespace akantu diff --git a/src/model/phase_field/phase_field_model.hh b/src/model/phase_field/phase_field_model.hh index e94567d49..e3b509b2a 100644 --- a/src/model/phase_field/phase_field_model.hh +++ b/src/model/phase_field/phase_field_model.hh @@ -1,325 +1,328 @@ /** * @file phase_field_model.hh * * @author Mohit Pundir * * @date creation: Sun Jul 30 2018 * @date last modification: Mon Feb 05 2018 * * @brief Model of Phase Field * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "boundary_condition.hh" #include "data_accessor.hh" #include "fe_engine.hh" #include "model.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PHASE_FIELD_MODEL_HH__ #define __AKANTU_PHASE_FIELD_MODEL_HH__ namespace akantu { template class IntegratorGauss; template class ShapeLagrange; } // namespace akantu namespace akantu { class PhaseFieldModel : public Model, public DataAccessor, public DataAccessor, public BoundaryCondition { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: using FEEngineType = FEEngineTemplate; PhaseFieldModel(Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "phase_field_model", const MemoryID & memory_id = 0, const ModelType model_type = ModelType::_phase_field_model); ~PhaseFieldModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// generic function to initialize everything ready for explicit dynamics void initFullImpl(const ModelOptions & options) override; /// read one material file to instantiate all the materials void readMaterials(); /// allocate all vectors void initSolver(TimeStepSolverType, NonLinearSolverType) override; /// initialize the model void initModel() override; /// predictor void predictor() override; /// corrector void corrector() override; /// compute the heat flux void assembleResidual() override; /// get the type of matrix needed MatrixType getMatrixType(const ID &) override; /// callback to assemble a Matrix void assembleMatrix(const ID &) override; /// callback to assemble a lumped Matrix void assembleLumpedMatrix(const ID &) override; std::tuple getDefaultSolverID(const AnalysisMethod & method) override; ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const; /// void updateInternalParameters(); /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /* ------------------------------------------------------------------------ */ /* Methods for static */ /* ------------------------------------------------------------------------ */ public: /// assembles the phasefield stiffness matrix virtual void assembleStiffnessMatrix(); /// compute the internal forces virtual void assembleInternalForces(); + + // compute the internal forces + void assembleInternalForces(const GhostType & ghost_type); /// assemble damage matrix void assembleDamageMatrix(); /// assemble damage gradient matrix void assembleDamageGradMatrix(); /// coupling parameters damage and strains from solid mechanics model void setCouplingParameters(ElementTypeMapArray & strain_on_qpoints, Array & damage); private: /// compute vector strain history field for each quadrature point void computePhiHistoryOnQuadPoints(const GhostType & ghost_type); /// compute vector strain history field for each quadrature point void computeDamageEnergyDensityOnQuadPoints(const GhostType & ghost_type); /// compute driving force for each quadrature point void computeDrivingForce(const GhostType & ghost_type); /// compute the damage on quadrature points void computeDamageOnQuadPoints(const GhostType & ghost_type); /// compute the fracture energy Real computeFractureEnergyByNode(); /* ------------------------------------------------------------------------ */ /* Methods for static */ /* ------------------------------------------------------------------------ */ public: /// set the stable timestep void setTimeStep(Real time_step, const ID & solver_id = "") override; protected: /// callback for the solver, this is called at beginning of solve void beforeSolveStep() override; /// callback for the solver, this is called at end of solve void afterSolveStep(bool converged = true) override; /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; UInt getNbData(const Array & indexes, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt); /// return the damage array AKANTU_GET_MACRO(Damage, *damage, Array &); /// get the ContactMechanics::internal_force vector (internal forces) AKANTU_GET_MACRO(InternalForce, *internal_force, Array &); /// get the ContactMechanicsModel::external_force vector (external forces) AKANTU_GET_MACRO(ExternalForce, *external_force, Array &); /// get the ContactMechanicsModel::force vector (external forces) Array & getForce() { AKANTU_DEBUG_WARNING("getForce was maintained for backward compatibility, " "use getExternalForce instead"); return *external_force; } /// AKANTU_GET_MACRO_NOT_CONST(Strain, strain_on_qpoints, ElementTypeMapArray &); /// get the PhaseFieldModel::blocked_dofs vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); /* ------------------------------------------------------------------------ */ /* Dumpable Interface */ /* ------------------------------------------------------------------------ */ public: std::shared_ptr createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) override; virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); void dump() override; virtual void dump(UInt step); virtual void dump(Real time, UInt step); protected: /* ------------------------------------------------------------------------ */ FEEngine & getFEEngineBoundary(const ID & name = "") override; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// number of iterations UInt n_iter; /// damage array Array * damage{nullptr}; /// damage array at the previous time step (used in finite deformation) Array * previous_damage{nullptr}; /// increment of damage Array * damage_increment{nullptr}; /// damage field on quadrature points ElementTypeMapArray damage_on_qpoints; /// critical local damage energy on quadrature points for \mathbf{B}^t * /// \mathbf{W} * \mathbf{B}@f$ ElementTypeMapArray damage_energy_on_qpoints; /// critical local damage energy density on quadrature points for /// \mathbf{N}^t * \mathbf{w} * \mathbf{N}@f$ ElementTypeMapArray damage_energy_density_on_qpoints; /// the speed of change in damage ElementTypeMapArray damage_gradient; /// strain on quadrature points ElementTypeMapArray strain_on_qpoints; /// driving force on quadrature points for internal forces ElementTypeMapArray driving_force_on_qpoints; /// vector \phi plus on quadrature points ElementTypeMapArray phi_history_on_qpoints; /// boundary vector Array * blocked_dofs{nullptr}; /// external force vector Array * external_force{nullptr}; /// residuals array Array * internal_force{nullptr}; /// lengthscale parameter Real l_0; /// critical energy release rate Real g_c; /// Young's modulus Real E; /// Poisson ratio Real nu; /// Lame's first parameter Real lambda; /// Lame's second paramter Real mu; }; } // namespace akantu /* ------------------------------------------------------------------------ */ /* inline functions */ /* ------------------------------------------------------------------------ */ -#include "phase_field_model_inline_impl.cc" +//#include "phase_field_model_inline_impl.cc" #endif