diff --git a/src/fem/element_group.cc b/src/fem/element_group.cc index 40b0861f1..93a83c3dc 100644 --- a/src/fem/element_group.cc +++ b/src/fem/element_group.cc @@ -1,93 +1,95 @@ /** * @file element_group.cc * * @author Dana Christen * * @date Wed Mar 06 09:30:00 2013 * * @brief Stores information relevent to the notion of domain boundary and surfaces. * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include #include #include "mesh.hh" #include "aka_csr.hh" #include "mesh_utils.hh" #include "element_group.hh" #if defined(AKANTU_USE_IOHELPER) # include "dumper_paraview.hh" #endif __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ ElementGroup::ElementGroup(const std::string & boundary_name, const Mesh & mesh, NodeGroup & node_group, UInt dimension, const std::string & id, const MemoryID & mem_id) : Memory(mem_id), mesh(mesh), id(id), name(boundary_name), elements("elements", id, mem_id), node_group(node_group), dimension(dimension) { AKANTU_DEBUG_IN(); +#if defined(AKANTU_USE_IOHELPER) this->registerDumper("paraview_" + boundary_name, boundary_name, true); this->addDumpFilteredMesh(mesh, elements, node_group.getNodes(), dimension); +#endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ElementGroup::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "ElementGroup [" << std::endl; stream << space << " + name: " << name << std::endl; elements.printself(stream, indent + 1); node_group.printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void ElementGroup::cleanUpNodeList() { // UInt size_before = nodes.getSize(); // UInt * newEnd; // UInt * begin = nodes.storage(); // UInt * end = nodes.storage()+nodes.getSize(); // std::sort(begin, end); // newEnd = std::unique(begin, end); // UInt crop_size = end-newEnd; // UInt size_after = size_before - crop_size; // nodes.resize(size_after); } __END_AKANTU__ diff --git a/src/io/parser/parsable.hh b/src/io/parser/parsable.hh index f808acd9a..b011888fa 100644 --- a/src/io/parser/parsable.hh +++ b/src/io/parser/parsable.hh @@ -1,175 +1,173 @@ /** * @file parsable.hh * * @author Nicolas Richart * * @date Fri Aug 2 13:13:31 2013 * * @brief Interface of parsable objects * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "parser.hh" /* -------------------------------------------------------------------------- */ - #ifndef __AKANTU_PARSABLE_HH__ #define __AKANTU_PARSABLE_HH__ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ enum ParamAccessType { _pat_internal = 0x0001, _pat_writable = 0x0010, _pat_readable = 0x0100, _pat_modifiable = 0x0110, //_pat_readable | _pat_writable, _pat_parsable = 0x1000, _pat_parsmod = 0x1110 }; inline ParamAccessType operator|(const ParamAccessType & a, const ParamAccessType & b) { ParamAccessType tmp = ParamAccessType(UInt(a) | UInt(b)); return tmp; } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template class ParsableParamTyped; /* -------------------------------------------------------------------------- */ /** * Interface for the ParsableParamTyped */ class ParsableParam { public: ParsableParam(); ParsableParam(std::string name, std::string description, ParamAccessType param_type); virtual ~ParsableParam() {}; /* ------------------------------------------------------------------------ */ bool isInternal() const; bool isWritable() const; bool isReadable() const; bool isParsable() const; /* ------------------------------------------------------------------------ */ template void set(T & value); template const T & get() const; virtual void parseParam(const ParserParameter & param); /* ------------------------------------------------------------------------ */ virtual void printself(std::ostream & stream) const; protected: template const ParsableParamTyped & getParsableParamTyped() const; template ParsableParamTyped & getParsableParamTyped(); private: std::string name; std::string description; ParamAccessType param_type; }; /* -------------------------------------------------------------------------- */ /* Typed Parameter */ /* -------------------------------------------------------------------------- */ /** * Type parameter transfering a ParserParameter (string: string) to a typed parameter in the memory of the p */ template class ParsableParamTyped : public ParsableParam { public: ParsableParamTyped(std::string name, std::string description, ParamAccessType param_type, T & param); /* ------------------------------------------------------------------------ */ void setTyped(T & value); const T & getTyped() const; void parseParam(const ParserParameter & param); virtual void printself(std::ostream & stream) const; private: T & param; }; - /* -------------------------------------------------------------------------- */ /* Parsable Interface */ /* -------------------------------------------------------------------------- */ class Parsable { public: Parsable(const SectionType & section_type, const ID & id = std::string()) : section_type(section_type), id(id) {}; virtual ~Parsable(); /* ------------------------------------------------------------------------ */ template void registerParam(std::string name, T & variable, ParamAccessType type, const std::string description = ""); template void registerParam(std::string name, T & variable, T default_value, ParamAccessType type, const std::string description = ""); void registerSubSection(const SectionType & type, const std::string & name, Parsable & sub_section); /* ------------------------------------------------------------------------ */ template void setParam(std::string name, T value); template const T & getParam(std::string name) const; /* ------------------------------------------------------------------------ */ virtual void parseSection(const ParserSection & section); virtual void parseSubSection(const ParserSection & section); virtual void parseParam(const ParserParameter & parameter); /* ------------------------------------------------------------------------ */ virtual void printself(std::ostream & stream, int indent) const; private: SectionType section_type; ID id; std::map params; typedef std::pair SubSectionKey; typedef std::map SubSections; SubSections sub_sections; }; __END_AKANTU__ #include "parsable_tmpl.hh" #endif /* __AKANTU_PARSABLE_HH__ */ diff --git a/src/io/parser/parser.hh b/src/io/parser/parser.hh index d4f728da9..f2b3abd01 100644 --- a/src/io/parser/parser.hh +++ b/src/io/parser/parser.hh @@ -1,369 +1,368 @@ - /** * @file boost_spirit_parser.cc * * @author Nicolas Richart * * @date Wed Jul 24 11:34:23 2013 * * @brief * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_random_generator.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PARSER_HH__ #define __AKANTU_PARSER_HH__ __BEGIN_AKANTU__ #define AKANTU_SECTION_TYPES \ (global) \ (material) \ (model) \ (heat) \ (rules) \ (non_local) \ (not_defined) #define AKANTU_SECTION_TYPES_PREFIX(elem) BOOST_PP_CAT(_st_, elem) #define AKANTU_SECT_PREFIX(s, data, elem) AKANTU_SECTION_TYPES_PREFIX(elem) enum SectionType { BOOST_PP_SEQ_ENUM(BOOST_PP_SEQ_TRANSFORM(AKANTU_SECT_PREFIX, _, AKANTU_SECTION_TYPES)) }; #undef AKANTU_SECT_PREFIX #define AKANTU_SECTION_TYPE_PRINT_CASE(r, data, elem) \ case AKANTU_SECTION_TYPES_PREFIX(elem): { \ stream << BOOST_PP_STRINGIZE(AKANTU_SECTION_TYPES_PREFIX(elem)); \ break; \ } inline std::ostream & operator <<(std::ostream & stream, SectionType type) { switch(type) { BOOST_PP_SEQ_FOR_EACH(AKANTU_SECTION_TYPE_PRINT_CASE, _, AKANTU_SECTION_TYPES) default: stream << "not a SectionType"; break; } return stream; } #undef AKANTU_SECTION_TYPE_PRINT_CASE enum ParserParameterSearchCxt { _ppsc_current_scope = 0x1, _ppsc_parent_scope = 0x2, _ppsc_current_and_parent_scope = 0x3 }; /* ------------------------------------------------------------------------ */ /* Parameters Class */ /* ------------------------------------------------------------------------ */ class ParserSection; class ParserParameter { public: ParserParameter() : parent_section(NULL), name(std::string()), value(std::string()), dbg_filename(std::string()) {} ParserParameter(const std::string & name, const std::string & value, const ParserSection & parent_section) : parent_section(&parent_section), name(name), value(value), dbg_filename(std::string()) {} ParserParameter(const ParserParameter & param) : parent_section(param.parent_section), name(param.name), value(param.value), dbg_filename(param.dbg_filename), dbg_line(param.dbg_line), dbg_column(param.dbg_column) {} virtual ~ParserParameter() {} const std::string & getName() const { return name; } const std::string & getValue() const { return value; } void setDebugInfo(const std::string & filename, UInt line, UInt column) { dbg_filename = filename; dbg_line = line; dbg_column = column; } template inline operator T() const; // template inline operator Vector() const; // template inline operator Matrix() const; void printself(std::ostream & stream, unsigned int indent = 0) const { stream << name << ": " << value << " (" << dbg_filename << ":" << dbg_line << ":" << dbg_column << ")"; } private: void setParent(const ParserSection & sect) { parent_section = § } friend ParserSection; private: const ParserSection * parent_section; std::string name; std::string value; std::string dbg_filename; UInt dbg_line, dbg_column; }; /* ------------------------------------------------------------------------ */ /* Sections Class */ /* ------------------------------------------------------------------------ */ class ParserSection { public: typedef std::multimap SubSections; typedef std::map Parameters; private: typedef SubSections::const_iterator const_section_iterator_; public: /* ------------------------------------------------------------------------ */ /* SubSection iterator */ /* ------------------------------------------------------------------------ */ class const_section_iterator { public: const_section_iterator(const const_section_iterator & other) : it(other.it) { } const_section_iterator(const const_section_iterator_ & it) : it(it) { } const_section_iterator & operator=(const const_section_iterator & other) { if(this != &other) { it = other.it; } return *this; } const ParserSection & operator* () const { return it->second; } const ParserSection * operator-> () const { return &(it->second); } bool operator==(const const_section_iterator & other) const { return it == other.it; } bool operator!=(const const_section_iterator & other) const { return it != other.it; } const_section_iterator & operator++() { ++it; return *this; } const_section_iterator operator++(int) { const_section_iterator tmp = *this; operator++(); return tmp; } private: const_section_iterator_ it; }; /* ------------------------------------------------------------------------ */ /* Parameters iterator */ /* ------------------------------------------------------------------------ */ class const_parameter_iterator { public: const_parameter_iterator(const const_parameter_iterator & other) : it(other.it) { } const_parameter_iterator(const Parameters::const_iterator & it) : it(it) { } const_parameter_iterator & operator=(const const_parameter_iterator & other) { if(this != &other) { it = other.it; } return *this; } const ParserParameter & operator* () const { return it->second; } const ParserParameter * operator->() { return &(it->second); }; bool operator==(const const_parameter_iterator & other) const { return it == other.it; } bool operator!=(const const_parameter_iterator & other) const { return it != other.it; } const_parameter_iterator & operator++() { ++it; return *this; } const_parameter_iterator operator++(int) { const_parameter_iterator tmp = *this; operator++(); return tmp; } private: Parameters::const_iterator it; }; /* ---------------------------------------------------------------------- */ ParserSection() : parent_section(NULL), name(std::string()), type(_st_not_defined){} ParserSection(const std::string & name, SectionType type) : parent_section(NULL), name(name), type(type) {} ParserSection(const std::string & name, SectionType type, const std::string & option, const ParserSection & parent_section) : parent_section(&parent_section), name(name), type(type), option(option) {} ParserSection(const ParserSection & section) : parent_section(section.parent_section), name(section.name), type(section.type), option(section.option), parameters(section.parameters), sub_sections_by_type(section.sub_sections_by_type) { setChldrenPointers(); } ParserSection & operator=(const ParserSection & other) { if(&other != this) { parent_section = other.parent_section; name = other.name; type = other.type; option = other.option; parameters = other.parameters; sub_sections_by_type = sub_sections_by_type; setChldrenPointers(); } return *this; } virtual ~ParserSection(); virtual void printself(std::ostream & stream, unsigned int indent = 0) const; /* ---------------------------------------------------------------------- */ /* Creation functions */ /* ---------------------------------------------------------------------- */ public: ParserParameter & addParameter(const ParserParameter & param); ParserSection & addSubSection(const ParserSection & section); private: void setChldrenPointers() { Parameters::iterator pit = this->parameters.begin(); for(;pit != this->parameters.end(); ++pit) pit->second.setParent(*this); SubSections::iterator sit = this->sub_sections_by_type.begin(); for (;sit != this->sub_sections_by_type.end(); ++sit) sit->second.setParent(*this); } /* ---------------------------------------------------------------------- */ /* Accessors */ /* ---------------------------------------------------------------------- */ public: std::pair getSubSections(SectionType type = _st_not_defined) const { if(type != _st_not_defined) { std::pair range = sub_sections_by_type.equal_range(type); return std::pair(range.first, range.second); } else { return std::pair(sub_sections_by_type.begin(), sub_sections_by_type.end()); } } std::pair getParameters() const { return std::pair(parameters.begin(), parameters.end()); } /* ---------------------------------------------------------------------- */ const ParserParameter & getParameter(const std::string & name, ParserParameterSearchCxt search_ctx = _ppsc_current_scope) const { Parameters::const_iterator it; if(search_ctx & _ppsc_current_scope) it = parameters.find(name); if(it == parameters.end()) { if((search_ctx & _ppsc_parent_scope) && parent_section) return parent_section->getParameter(name, search_ctx); else { AKANTU_SILENT_EXCEPTION("The parameter " << name << " has not been found in the specified context"); } } return it->second; } const std::string & getName() const { return name; } const SectionType & getType() const { return type; } const std::string & getOption() const { return option; } protected: void setParent(const ParserSection & sect) { parent_section = § } /* ---------------------------------------------------------------------- */ /* Members */ /* ---------------------------------------------------------------------- */ private: const ParserSection * parent_section; std::string name; SectionType type; std::string option; Parameters parameters; SubSections sub_sections_by_type; }; /* ------------------------------------------------------------------------ */ /* Parser Class */ /* ------------------------------------------------------------------------ */ class Parser : public ParserSection { public: Parser() : ParserSection("global", _st_global) {} void parse(const std::string & filename); public: static Real parseReal (const std::string & value, const ParserSection & section); static Vector parseVector(const std::string & value, const ParserSection & section); static Matrix parseMatrix(const std::string & value, const ParserSection & section); static RandomParameter parseRandomParameter(const std::string & value, const ParserSection & section); protected: template static T parseType(const std::string & value, Grammar & grammar); public: static bool parser_permissive; }; inline std::ostream & operator <<(std::ostream & stream, const ParserParameter &_this) { _this.printself(stream); return stream; } inline std::ostream & operator <<(std::ostream & stream, const ParserSection & section) { section.printself(stream); return stream; } __END_AKANTU__ #include "parser_tmpl.hh" #endif /* __AKANTU_PARSER_HH__ */ diff --git a/src/model/solid_mechanics/material.hh b/src/model/solid_mechanics/material.hh index bae8e3054..16b4a17d6 100644 --- a/src/model/solid_mechanics/material.hh +++ b/src/model/solid_mechanics/material.hh @@ -1,596 +1,597 @@ /** * @file material.hh * * @author Marco Vocialta * @author Nicolas Richart * * @date Tue Jul 27 18:15:37 2010 * * @brief Mother class for all materials * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "parser.hh" #include "parsable.hh" #include "data_accessor.hh" #include "internal_field.hh" #include "random_internal_field.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_HH__ #define __AKANTU_MATERIAL_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class Model; class SolidMechanicsModel; } __BEGIN_AKANTU__ template class VoigtHelper { public: const static UInt size; // matrix of vector index I as function of tensor indices i,j const static UInt mat[dim][dim]; // array of matrix indices ij as function of vector index I const static UInt vec[dim*dim][2]; // factors to multiply the strain by for voigt notation const static Real factors[dim*(dim-(dim-1)/2)]; }; template const UInt VoigtHelper::size = dim*(dim-(dim-1)/2); /** * Interface of all materials * Prerequisites for a new material * - inherit from this class * - implement the following methods: * \code * virtual Real getStableTimeStep(Real h, const Element & element = ElementNull); * * virtual void computeStress(ElementType el_type, * GhostType ghost_type = _not_ghost); * * virtual void computeTangentStiffness(const ElementType & el_type, * Array & tangent_matrix, * GhostType ghost_type = _not_ghost); * \endcode * */ + class Material : public Memory, public DataAccessor, public Parsable, public MeshEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Material(SolidMechanicsModel & model, const ID & id = ""); virtual ~Material(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: template void registerInternal(__attribute__((unused)) InternalField & vect) { AKANTU_DEBUG_TO_IMPLEMENT(); } template void unregisterInternal(__attribute__((unused)) InternalField & vect) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// function called to update the internal parameters when the modifiable /// parameters are modified virtual void updateInternalParameters() {} /// initialize the material computed parameter virtual void initMaterial(); /// compute the residual for this material virtual void updateResidual(GhostType ghost_type = _not_ghost); /// assemble the residual for this material virtual void assembleResidual(GhostType ghost_type); /// save the stress in the precious_stress if needed virtual void savePreviousState(GhostType ghost_type); /// compute the stresses for this material virtual void computeAllStresses(GhostType ghost_type = _not_ghost); virtual void computeAllNonLocalStresses(__attribute__((unused)) GhostType ghost_type = _not_ghost) {}; virtual void computeAllStressesFromTangentModuli(GhostType ghost_type = _not_ghost); virtual void computeAllCauchyStresses(GhostType ghost_type = _not_ghost); /// set material to steady state void setToSteadyState(GhostType ghost_type = _not_ghost); /// compute the stiffness matrix virtual void assembleStiffnessMatrix(GhostType ghost_type); /// compute the stable time step for an element of size h virtual Real getStableTimeStep(__attribute__((unused)) Real h, __attribute__((unused)) const Element & element = ElementNull) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the p-wave speed in the material virtual Real getPushWaveSpeed() const { AKANTU_DEBUG_TO_IMPLEMENT(); }; /// compute the s-wave speed in the material virtual Real getShearWaveSpeed() const { AKANTU_DEBUG_TO_IMPLEMENT(); }; /// add an element to the local mesh filter inline UInt addElement(const ElementType & type, UInt element, const GhostType & ghost_type); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /** * interpolate stress on given positions for each element by means * of a geometrical interpolation on quadrature points */ virtual void interpolateStress(ByElementTypeReal & result, const GhostType ghost_type = _not_ghost); /** * function to initialize the elemental field interpolation * function by inverting the quadrature points' coordinates */ virtual void initElementalFieldInterpolation(ByElementTypeReal & interpolation_points_coordinates); protected: /// constitutive law virtual void computeStress(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// template void assembleResidual(GhostType ghost_type); /// Computation of Cauchy stress tensor in the case of finite deformation template void computeCauchyStress(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost); template inline void computeCauchyStressOnQuad(const Matrix & F, const Matrix & S, Matrix & cauchy); template void computeAllStressesFromTangentModuli(const ElementType & type, GhostType ghost_type); /// set the material to steady state (to be implemented for materials that need it) virtual void setToSteadyState(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) {} /// compute the tangent stiffness matrix virtual void computeTangentModuli(__attribute__((unused)) const ElementType & el_type, __attribute__((unused)) Array & tangent_matrix, __attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the potential energy virtual void computePotentialEnergy(ElementType el_type, GhostType ghost_type = _not_ghost); template void assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type); /// assembling in finite deformation template void assembleStiffnessMatrixNL(const ElementType & type, GhostType ghost_type); template void assembleStiffnessMatrixL2(const ElementType & type, GhostType ghost_type); /// transfer the B matrix to a Voigt notation B matrix template inline void transferBMatrixToSymVoigtBMatrix(const Matrix & B, Matrix & Bvoigt, UInt nb_nodes_per_element) const; /// transfer the BNL and BL2 matrix to a Voigt notation B matrix (See Bathe et al. IJNME vol 9, 1975) template inline void transferBMatrixToBNL(const Matrix & B, Matrix & Bvoigt, UInt nb_nodes_per_element) const; template inline void transferBMatrixToBL2(const Matrix & B, const Matrix & grad_u, Matrix & Bvoigt, UInt nb_nodes_per_element) const; /// write the stress tensor in the Voigt notation. template inline void SetCauchyStressArray(const Matrix & S_t, Matrix & Stress_vect); inline UInt getTangentStiffnessVoigtSize(UInt spatial_dimension) const; /// Size of the Stress matrix for the case of finite deformation see: Bathe et al, IJNME, Vol 9, 353-386, 1975 inline UInt getCauchyStressMatrixSize(UInt spatial_dimension) const; /// Sets the stress matrix according to Bathe et al, IJNME, Vol 9, 353-386, 1975 template inline void SetCauchyStressMatrix(const Matrix & S_t, Matrix & Stress_matrix); /// compute the potential energy by element void computePotentialEnergyByElement(); public: /// compute the coordinates of the quadrature points void computeQuadraturePointsCoordinates(ByElementTypeReal & quadrature_points_coordinates, const GhostType & ghost_type) const; protected: /// interpolate an elemental field on given points for each element template void interpolateElementalField(const Array & field, Array & result, const GhostType ghost_type); /// template function to initialize the elemental field interpolation template void initElementalFieldInterpolation(const Array & quad_coordinates, const Array & interpolation_points_coordinates, const GhostType ghost_type); /// build the coordinate matrix for the interpolation on elemental field template inline void buildElementalFieldInterpolationCoodinates(const Matrix & coordinates, Matrix & coordMatrix); /// build interpolation coordinates for basic linear elements inline void buildElementalFieldInterpolationCoodinatesLinear(const Matrix & coordinates, Matrix & coordMatrix); /// build interpolation coordinates for basic quadratic elements inline void buildElementalFieldInterpolationCoodinatesQuadratic(const Matrix & coordinates, Matrix & coordMatrix); /// get the size of the coordiante matrix used in the interpolation template inline UInt getSizeElementalFieldInterpolationCoodinates(GhostType ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Function for all materials */ /* ------------------------------------------------------------------------ */ protected: /// compute the potential energy for a quadrature point virtual inline void computePotentialEnergyOnQuad(Matrix & grad_u, Matrix & sigma, Real & epot); /// compute the potential energy for an element virtual void computePotentialEnergyByElement(ElementType type, UInt index, Vector & epot_on_quad_points); public: /* ------------------------------------------------------------------------ */ template inline void gradUToF(const Matrix & grad_u, Matrix & F); inline void rightCauchy(const Matrix & F, Matrix & C); inline void leftCauchy (const Matrix & F, Matrix & B); template inline void deformationJacobian(const Matrix & F, Real & J); inline void gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon); /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: virtual inline UInt getNbDataForElements(const Array & elements, SynchronizationTag tag) const; virtual inline void packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const; virtual inline void unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag); template inline void packElementDataHelper(const ByElementTypeArray & data_to_pack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id = ID()) const; template inline void unpackElementDataHelper(ByElementTypeArray & data_to_unpack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id = ID()); public: /* ------------------------------------------------------------------------ */ virtual inline void onElementsAdded(const Array & element_list, const NewElementsEvent & event); virtual inline void onElementsRemoved(const Array & element_list, const ByElementTypeUInt & new_numbering, const RemovedElementsEvent & event); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Model, *model, const SolidMechanicsModel &) AKANTU_GET_MACRO(ID, id, const ID &); AKANTU_GET_MACRO(Rho, rho, Real); AKANTU_SET_MACRO(Rho, rho, Real); /// return the potential energy for the subset of elements contained by the material Real getPotentialEnergy(); /// return the potential energy for the provided element Real getPotentialEnergy(ElementType & type, UInt index); /// return the energy (identified by id) for the subset of elements contained by the material virtual Real getEnergy(std::string energy_id); /// return the energy (identified by id) for the provided element virtual Real getEnergy(std::string energy_id, ElementType type, UInt index); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementFilter, element_filter, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Strain, strain, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PotentialEnergy, potential_energy, Real); AKANTU_GET_MACRO(Strain, strain, const ByElementTypeReal &); AKANTU_GET_MACRO(Stress, stress, const ByElementTypeReal &); AKANTU_GET_MACRO(ElementFilter, element_filter, const ByElementTypeUInt &); bool isNonLocal() const { return is_non_local; } const Array & getArray(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost) const; Array & getArray(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost); const ByElementTypeArray & getInternal(const ID & id) const; ByElementTypeArray & getInternal(const ID & id); bool isFiniteDeformation() const { return finite_deformation; } bool isInelasticDeformation() const { return inelastic_deformation; } protected: bool isInit() const { return is_init; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// boolean to know if the material has been initialized bool is_init; std::map *> internal_vectors_real; std::map *> internal_vectors_uint; protected: /// id of the material ID id; /// Finite deformation bool finite_deformation; /// Finite deformation bool inelastic_deformation; /// material name std::string name; /// The model to witch the material belong SolidMechanicsModel * model; /// density : rho Real rho; /// spatial dimension UInt spatial_dimension; /// list of element handled by the material ByElementTypeArray element_filter; /// stresses arrays ordered by element types InternalField stress; /// strains arrays ordered by element types InternalField strain; /// stress increment arrays ordered by element types (Finite deformation) InternalField delta_stress; /// strain increment arrays ordered by element types (Finite deformation) InternalField delta_strain; /// strain increment arrays ordered by element types (inelastic deformation) InternalField inelas_strain; /// Second Piola-Kirchhoff stress tensor arrays ordered by element types (Finite deformation) InternalField piola_kirchhoff_stress; /// potential energy by element InternalField potential_energy; /// tell if using in non local mode or not bool is_non_local; /// tell if the material need the previous stress state bool use_previous_stress; /// previous stresses InternalField previous_stress; /// tell if the material need the previous strain state bool use_previous_strain; /// previous strain InternalField previous_strain; /// elemental field interpolation coordinates InternalField interpolation_inverse_coordinates; /// elemental field interpolation points InternalField interpolation_points_matrices; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_inline_impl.cc" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Material & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #include "internal_field_tmpl.hh" #include "random_internal_field_tmpl.hh" /* -------------------------------------------------------------------------- */ /* Auto loop */ /* -------------------------------------------------------------------------- */ #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type) \ Array::iterator< Matrix > strain_it = \ this->strain(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ Array::iterator< Matrix > strain_end = \ this->strain(el_type, ghost_type).end(spatial_dimension, \ spatial_dimension); \ \ this->stress(el_type, \ ghost_type).resize(this->strain(el_type, \ ghost_type).getSize()); \ \ Array::iterator< Matrix > stress_it = \ this->stress(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ \ if(this->isFiniteDeformation()) \ Array::iterator< Matrix > stress_it = \ this->piola_kirchhoff_stress(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ for(;strain_it != strain_end; ++strain_it, ++stress_it) { \ Matrix & __attribute__((unused)) grad_u = *strain_it; \ Matrix & __attribute__((unused)) sigma = *stress_it #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END \ } \ #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_mat) \ Array::iterator< Matrix > strain_it = \ this->strain(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ Array::iterator< Matrix > strain_end = \ this->strain(el_type, ghost_type).end(spatial_dimension, \ spatial_dimension); \ Array::iterator< Matrix > sigma_it = \ this->stress(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ \ tangent_mat.resize(this->strain(el_type, ghost_type).getSize()); \ \ UInt tangent_size = \ this->getTangentStiffnessVoigtSize(spatial_dimension); \ Array::iterator< Matrix > tangent_it = \ tangent_mat.begin(tangent_size, \ tangent_size); \ \ for(;strain_it != strain_end; ++strain_it, ++sigma_it, ++tangent_it) { \ Matrix & __attribute__((unused)) grad_u = *strain_it; \ Matrix & __attribute__((unused)) sigma_tensor = *sigma_it; \ Matrix & tangent = *tangent_it #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END \ } \ #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_FINITE_DEFORMATION_BEGIN(tangent_mat) \ Array::iterator< Matrix > strain_it = \ this->strain(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ Array::iterator< Matrix > strain_end = \ this->strain(el_type, ghost_type).end(spatial_dimension, \ spatial_dimension); \ Array::iterator< Matrix > green_it = \ this->delta_strain(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ Array::iterator< Matrix > sigma_it = \ this->stress(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ \ tangent_mat.resize(this->strain(el_type, ghost_type).getSize()); \ \ UInt tangent_size = \ this->getTangentStiffnessVoigtSize(spatial_dimension); \ Array::iterator< Matrix > tangent_it = \ tangent_mat.begin(tangent_size, \ tangent_size); \ \ for(;strain_it != strain_end; ++strain_it, ++green_it, ++sigma_it, ++tangent_it) { \ Matrix & __attribute__((unused)) grad_u = *strain_it; \ Matrix & __attribute__((unused)) grad_delta_u = *green_it; \ Matrix & __attribute__((unused)) sigma_tensor = *sigma_it; \ Matrix & tangent = *tangent_it #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_FINITE_DEFORMATION_END \ } \ /* -------------------------------------------------------------------------- */ #define INSTANSIATE_MATERIAL(mat_name) \ template class mat_name<1>; \ template class mat_name<2>; \ template class mat_name<3> #endif /* __AKANTU_MATERIAL_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index 11a90bbfc..20eb238d9 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1853 +1,1853 @@ /** * @file solid_mechanics_model.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date Tue Jul 27 18:15:37 2010 * * @brief Implementation of the SolidMechanicsModel class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_math.hh" #include "aka_common.hh" #include "solid_mechanics_model.hh" #include "integration_scheme_2nd_order.hh" -#include "dumper_paraview.hh" #include "element_group.hh" #include "static_communicator.hh" #include "dof_synchronizer.hh" #include #ifdef AKANTU_USE_MUMPS #include "solver_mumps.hh" #endif #ifdef AKANTU_USE_IOHELPER # include "dumper_paraview.hh" # include "dumper_iohelper_tmpl.hh" # include "dumper_iohelper_tmpl_homogenizing_field.hh" # include "dumper_iohelper_tmpl_material_internal_field.hh" #endif /* -------------------------------------------------------------------------- */ __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), Dumpable(), BoundaryCondition(), time_step(NAN), f_m2a(1.0), mass_matrix(NULL), velocity_damping_matrix(NULL), stiffness_matrix(NULL), jacobian_matrix(NULL), element_index_by_material("element index by material", id), material_selector(new DefaultMaterialSelector(element_index_by_material)), is_default_material_selector(true), integrator(NULL), increment_flag(false), solver(NULL), synch_parallel(NULL), are_materials_instantiated(false) { AKANTU_DEBUG_IN(); createSynchronizerRegistry(this); registerFEMObject("SolidMechanicsFEM", mesh, spatial_dimension); this->displacement = NULL; this->mass = NULL; this->velocity = NULL; this->acceleration = NULL; this->force = NULL; this->residual = NULL; this->boundary = NULL; this->increment = NULL; this->increment_acceleration = NULL; this->dof_synchronizer = NULL; this->previous_displacement = NULL; materials.clear(); mesh.registerEventHandler(*this); +#if defined(AKANTU_USE_IOHELPER) this->registerDumper("paraview_all", id, true); this->addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular); - +#endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModel::~SolidMechanicsModel() { AKANTU_DEBUG_IN(); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { delete *mat_it; } materials.clear(); delete integrator; if(solver) delete solver; if(mass_matrix) delete mass_matrix; if(velocity_damping_matrix) delete velocity_damping_matrix; if(stiffness_matrix && stiffness_matrix != jacobian_matrix) delete stiffness_matrix; if(jacobian_matrix) delete jacobian_matrix; if(dof_synchronizer) delete dof_synchronizer; if(synch_parallel) delete synch_parallel; if(is_default_material_selector) delete material_selector; AKANTU_DEBUG_OUT(); } void SolidMechanicsModel::setTimeStep(Real time_step) { this->time_step = time_step; #if defined(AKANTU_USE_IOHELPER) getDumper().setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ /** * 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(std::string input_file, const ModelOptions & options) { Model::initFull(input_file, options); const SolidMechanicsModelOptions & smm_options = dynamic_cast(options); method = smm_options.analysis_method; // initialize the vectors initArrays(); // set the initial condition to 0 force->clear(); velocity->clear(); acceleration->clear(); displacement->clear(); // initialize pcb if(pbc_pair.size()!=0) initPBC(); // initialize the time integration schemes switch(method) { case _explicit_lumped_mass: initExplicit(); break; case _explicit_consistent_mass: initSolver(); initExplicit(); break; case _implicit_dynamic: initImplicit(true); break; case _static: initImplicit(false); break; } // initialize the materials if(input_file != "") { instantiateMaterials(); } if(!smm_options.no_init_materials) { initMaterials(); } if(increment_flag) initBC(*this, *displacement, *increment, *force); else initBC(*this, *displacement, *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); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initFEMBoundary() { FEM & fem_boundary = getFEMBoundary(); fem_boundary.initShapeFunctions(_not_ghost); fem_boundary.initShapeFunctions(_ghost); fem_boundary.computeNormalsOnControlPoints(_not_ghost); fem_boundary.computeNormalsOnControlPoints(_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(sstr.str(), nb_nodes, nb_degree_of_freedom, Real())); AKANTU_DEBUG_OUT(); } void SolidMechanicsModel::initArraysPreviousDisplacment() { AKANTU_DEBUG_IN(); SolidMechanicsModel::setIncrementFlagOn(); UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_disp_t; sstr_disp_t << id << ":previous_displacement"; previous_displacement = &(alloc (sstr_disp_t.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); 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_mass; sstr_mass << id << ":mass"; std::stringstream sstr_velo; sstr_velo << id << ":velocity"; std::stringstream sstr_acce; sstr_acce << id << ":acceleration"; std::stringstream sstr_forc; sstr_forc << id << ":force"; std::stringstream sstr_resi; sstr_resi << id << ":residual"; std::stringstream sstr_boun; sstr_boun << id << ":boundary"; displacement = &(alloc(sstr_disp.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); // mass = &(alloc(sstr_mass.str(), nb_nodes, spatial_dimension, 0)); velocity = &(alloc(sstr_velo.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); acceleration = &(alloc(sstr_acce.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); force = &(alloc(sstr_forc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); residual = &(alloc(sstr_resi.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); boundary = &(alloc(sstr_boun.str(), nb_nodes, spatial_dimension, false)); std::stringstream sstr_curp; sstr_curp << id << ":current_position"; current_position = &(alloc(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE)); 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); element_index_by_material.alloc(nb_element, 2, *it, gt); } } dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); 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 getFEM().initShapeFunctions(_not_ghost); getFEM().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::iterator it = pbc_pair.begin(); std::map::iterator end = pbc_pair.end(); UInt dim = mesh.getSpatialDimension(); while(it != end) { for (UInt i=0; iregisterSynchronizer(*synch, _gst_smm_uv); synch_registry->registerSynchronizer(*synch, _gst_smm_mass); synch_registry->registerSynchronizer(*synch, _gst_smm_res); changeLocalEquationNumberForPBC(pbc_pair, mesh.getSpatialDimension()); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateCurrentPosition() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); current_position->resize(nb_nodes); Real * current_position_val = current_position->values; Real * position_val = mesh.getNodes().values; Real * displacement_val = displacement->values; /// compute current_position = initial_position + displacement memcpy(current_position_val, position_val, nb_nodes*spatial_dimension*sizeof(Real)); for (UInt n = 0; n < nb_nodes*spatial_dimension; ++n) { *current_position_val++ += *displacement_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initializeUpdateResidualData() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); residual->resize(nb_nodes); /// copy the forces in residual for boundary conditions memcpy(residual->values, force->values, nb_nodes*spatial_dimension*sizeof(Real)); // start synchronization synch_registry->asynchronousSynchronize(_gst_smm_uv); synch_registry->waitEndSynchronize(_gst_smm_uv); updateCurrentPosition(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Explicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /** * This function compute the second member of the motion equation. That is to * say the sum of forces @f$ r = F_{ext} - F_{int} @f$. @f$ F_{ext} @f$ is given * by the user in the force vector, and @f$ F_{int} @f$ is computed as @f$ * F_{int} = \int_{\Omega} N \sigma d\Omega@f$ * */ void SolidMechanicsModel::updateResidual(bool need_initialize) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); // f = f_ext - f_int // f = f_ext if(need_initialize) initializeUpdateResidualData(); AKANTU_DEBUG_INFO("Compute local stresses"); std::vector::iterator mat_it; for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.savePreviousState(_not_ghost); mat.savePreviousState(_ghost); mat.computeAllStresses(_not_ghost); if(mat.isFiniteDeformation()) mat.computeAllCauchyStresses(_not_ghost); } #ifdef AKANTU_DAMAGE_NON_LOCAL /* ------------------------------------------------------------------------ */ /* Computation of the non local part */ synch_registry->asynchronousSynchronize(_gst_mnl_for_average); AKANTU_DEBUG_INFO("Compute non local stresses for local elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_not_ghost); } AKANTU_DEBUG_INFO("Wait distant non local stresses"); synch_registry->waitEndSynchronize(_gst_mnl_for_average); AKANTU_DEBUG_INFO("Compute non local stresses for ghosts elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_ghost); } #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.assembleResidual(_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.assembleResidual(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ 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::iterator mat_it; for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.savePreviousState(_not_ghost); mat.savePreviousState(_ghost); mat.computeAllStresses(_not_ghost); if(mat.isFiniteDeformation()) mat.computeAllCauchyStresses(_not_ghost); } /* ------------------------------------------------------------------------ */ #ifdef AKANTU_DAMAGE_NON_LOCAL /* Computation of the non local part */ synch_registry->asynchronousSynchronize(_gst_mnl_for_average); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_not_ghost); } synch_registry->waitEndSynchronize(_gst_mnl_for_average); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_ghost); } #endif } else { std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllStressesFromTangentModuli(_not_ghost); } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateResidualInternal() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Update the residual"); // f = f_ext - f_int - Ma - Cv = r - Ma - Cv; if(method != _static) { // f -= Ma if(mass_matrix) { // if full mass_matrix Array * Ma = new Array(*acceleration, true, "Ma"); *Ma *= *mass_matrix; /// \todo check unit conversion for implicit dynamics // *Ma /= f_m2a *residual -= *Ma; delete Ma; } else if (mass) { // else lumped mass UInt nb_nodes = acceleration->getSize(); UInt nb_degree_of_freedom = acceleration->getNbComponent(); Real * mass_val = mass->values; Real * accel_val = acceleration->values; Real * res_val = residual->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) { if(!(*boundary_val)) { *res_val -= *accel_val * *mass_val /f_m2a; } boundary_val++; res_val++; mass_val++; accel_val++; } } else { AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix."); } // f -= Cv if(velocity_damping_matrix) { Array * Cv = new Array(*velocity); *Cv *= *velocity_damping_matrix; *residual -= *Cv; delete Cv; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::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, *boundary, f_m2a); } else if (method == _explicit_consistent_mass) { solve(*increment_acceleration); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveLumped(Array & x, const Array & A, const Array & b, const Array & boundary, Real alpha) { Real * A_val = A.storage(); Real * b_val = b.storage(); Real * x_val = x.storage(); bool * boundary_val = boundary.storage(); UInt nb_degrees_of_freedom = x.getSize() * x.getNbComponent(); for (UInt n = 0; n < nb_degrees_of_freedom; ++n) { if(!(*boundary_val)) { *x_val = alpha * (*b_val / *A_val); } x_val++; A_val++; b_val++; boundary_val++; } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitPred() { AKANTU_DEBUG_IN(); if(increment_flag) increment->copy(*displacement); if(previous_displacement) previous_displacement->copy(*displacement); AKANTU_DEBUG_ASSERT(integrator,"itegrator should have been allocated: " << "have called initExplicit ? " << "or initImplicit ?"); integrator->integrationSchemePred(time_step, *displacement, *velocity, *acceleration, *boundary); if(increment_flag) { Real * inc_val = increment->values; Real * dis_val = displacement->values; UInt nb_degree_of_freedom = displacement->getSize() * displacement->getNbComponent(); for (UInt n = 0; n < nb_degree_of_freedom; ++n) { *inc_val = *dis_val - *inc_val; inc_val++; dis_val++; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitCorr() { AKANTU_DEBUG_IN(); integrator->integrationSchemeCorrAccel(time_step, *displacement, *velocity, *acceleration, *boundary, *increment_acceleration); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveStep() { AKANTU_DEBUG_IN(); this->explicitPred(); this->updateResidual(); this->updateAcceleration(); this->explicitCorr(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Implicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /** * Initialize the solver and create the sparse matrices needed. * */ void SolidMechanicsModel::initSolver(__attribute__((unused)) SolverOptions & options) { #if !defined(AKANTU_USE_MUMPS) // or other solver in the future \todo add AKANTU_HAS_SOLVER in CMake AKANTU_DEBUG_ERROR("You should at least activate one solver."); #else UInt nb_global_nodes = mesh.getNbGlobalNodes(); delete jacobian_matrix; std::stringstream sstr; sstr << id << ":jacobian_matrix"; jacobian_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, _symmetric, spatial_dimension, sstr.str(), memory_id); jacobian_matrix->buildProfile(mesh, *dof_synchronizer); if (!isExplicit()) { delete stiffness_matrix; std::stringstream sstr_sti; sstr_sti << id << ":stiffness_matrix"; stiffness_matrix = new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id); } #ifdef AKANTU_USE_MUMPS std::stringstream sstr_solv; sstr_solv << id << ":solver"; solver = new SolverMumps(*jacobian_matrix, sstr_solv.str()); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS if(solver) solver->initialize(options); #endif //AKANTU_HAS_SOLVER } /* -------------------------------------------------------------------------- */ /** * Initialize the implicit solver, either for dynamic or static cases, * * @param dynamic */ void SolidMechanicsModel::initImplicit(bool dynamic, SolverOptions & solver_options) { AKANTU_DEBUG_IN(); method = dynamic ? _implicit_dynamic : _static; if (!increment) setIncrementFlagOn(); initSolver(solver_options); if(method == _implicit_dynamic) { if(integrator) delete integrator; integrator = new TrapezoidalRule2(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initialAcceleration() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving Ma = f"); Solver * acc_solver = NULL; std::stringstream sstr; sstr << id << ":tmp_mass_matrix"; SparseMatrix * tmp_mass = new SparseMatrix(*mass_matrix, sstr.str(), memory_id); #ifdef AKANTU_USE_MUMPS std::stringstream sstr_solver; sstr << id << ":solver_mass_matrix"; acc_solver = new SolverMumps(*mass_matrix, sstr_solver.str()); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS acc_solver->initialize(); tmp_mass->applyBoundary(*boundary); acc_solver->setRHS(*residual); acc_solver->solve(*acceleration); delete acc_solver; delete tmp_mass; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix."); stiffness_matrix->clear(); // call compute stiffness matrix on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->assembleStiffnessMatrix(_not_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveDynamic() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL, "You should first initialize the implicit solver and assemble the stiffness matrix"); AKANTU_DEBUG_ASSERT(mass_matrix != NULL, "You should first initialize the implicit solver and assemble the mass matrix"); // updateResidualInternal(); solve(*increment); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveStatic() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving Ku = f"); AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL, "You should first initialize the implicit solver and assemble the stiffness matrix"); UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes; jacobian_matrix->copyContent(*stiffness_matrix); jacobian_matrix->applyBoundary(*boundary); increment->clear(); solver->setRHS(*residual); solver->solve(*increment); Real * increment_val = increment->storage(); Real * displacement_val = displacement->storage(); bool * boundary_val = boundary->storage(); for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++displacement_val, ++increment_val, ++boundary_val) { if (!(*boundary_val)) { *displacement_val += *increment_val; } else { *increment_val = 0.0; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveStatic(Array & boundary_normal, Array & EulerAngles) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving Ku = f"); AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL, "You should first initialize the implicit solver and assemble the stiffness matrix"); UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent(); // if(method != _static) jacobian_matrix->copyContent(*stiffness_matrix); Array * residual_rotated = new Array (nb_nodes, nb_degree_of_freedom, "residual_rotated"); //stiffness_matrix->saveMatrix("stiffness_original.out"); jacobian_matrix->applyBoundaryNormal(boundary_normal, EulerAngles, *residual, (*stiffness_matrix).getA(), *residual_rotated); //jacobian_matrix->saveMatrix("stiffness_rotated_dir.out"); jacobian_matrix->applyBoundary(*boundary); solver->setRHS(*residual_rotated); delete residual_rotated; if (!increment) setIncrementFlagOn(); solver->solve(*increment); Matrix T(nb_degree_of_freedom, nb_degree_of_freedom); Matrix small_rhs(nb_degree_of_freedom, nb_degree_of_freedom); Matrix T_small_rhs(nb_degree_of_freedom, nb_degree_of_freedom); Real * increment_val = increment->values; Real * displacement_val = displacement->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes; ++n) { bool constrain_ij = false; for (UInt j = 0; j < nb_degree_of_freedom; j++) { if (boundary_normal(n, j)) { constrain_ij = true; break; } } if (constrain_ij) { if (nb_degree_of_freedom == 2) { Real Theta = EulerAngles(n, 0); T(0, 0) = cos(Theta); T(0, 1) = -sin(Theta); T(1, 1) = cos(Theta); T(1, 0) = sin(Theta); } else if (nb_degree_of_freedom == 3) { Real Theta_x = EulerAngles(n, 0); Real Theta_y = EulerAngles(n, 1); Real Theta_z = EulerAngles(n, 2); T(0, 0) = cos(Theta_y) * cos(Theta_z); T(0, 1) = -cos(Theta_y) * sin(Theta_z); T(0, 2) = sin(Theta_y); T(1, 0) = cos(Theta_x) * sin(Theta_z) + cos(Theta_z) * sin(Theta_x) * sin(Theta_y); T(1, 1) = cos(Theta_x) * cos(Theta_z) - sin(Theta_x) * sin(Theta_y) * sin(Theta_z); T(1, 2) = -cos(Theta_y) * sin(Theta_x); T(2, 0) = sin(Theta_x) * sin(Theta_z) - cos(Theta_x) * cos(Theta_z) * sin(Theta_y); T(2, 1) = cos(Theta_z) * sin(Theta_x) + cos(Theta_x) * sin(Theta_y) * sin(Theta_z); T(2, 2) = cos(Theta_x) * cos(Theta_y); } small_rhs.clear(); T_small_rhs.clear(); for (UInt j = 0; j < nb_degree_of_freedom; j++) if(!(boundary_normal(n,j)) ) small_rhs(j,j)=increment_val[j]; T_small_rhs.mul(T,small_rhs); for (UInt j = 0; j < nb_degree_of_freedom; j++){ if(!(boundary_normal(n,j))){ for (UInt k = 0; k < nb_degree_of_freedom; k++) displacement_val[k]+=T_small_rhs(k,j); } } displacement_val += nb_degree_of_freedom; boundary_val += nb_degree_of_freedom; increment_val += nb_degree_of_freedom; } else { for (UInt j = 0; j < nb_degree_of_freedom; j++, ++displacement_val, ++increment_val, ++boundary_val) { if (!(*boundary_val)) { *displacement_val += *increment_val; } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SparseMatrix & SolidMechanicsModel::initVelocityDampingMatrix() { if(!velocity_damping_matrix) velocity_damping_matrix = new SparseMatrix(*jacobian_matrix, id + ":velocity_damping_matrix", memory_id); return *velocity_damping_matrix; } /* -------------------------------------------------------------------------- */ template<> bool SolidMechanicsModel::testConvergence<_scc_increment>(Real tolerance, Real & error){ AKANTU_DEBUG_IN(); UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent(); error = 0; Real norm[2] = {0., 0.}; Real * increment_val = increment->storage(); bool * boundary_val = boundary->storage(); Real * displacement_val = displacement->storage(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); for (UInt d = 0; d < nb_degree_of_freedom; ++d) { if(!(*boundary_val) && is_local_node) { norm[0] += *increment_val * *increment_val; norm[1] += *displacement_val * *displacement_val; } boundary_val++; increment_val++; displacement_val++; } } StaticCommunicator::getStaticCommunicator().allReduce(norm, 2, _so_sum); norm[0] = sqrt(norm[0]); norm[1] = sqrt(norm[1]); AKANTU_DEBUG_ASSERT(!Math::isnan(norm[0]), "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); error = norm[0] / norm[1]; return (error < tolerance); } /* -------------------------------------------------------------------------- */ template<> bool SolidMechanicsModel::testConvergence<_scc_residual>(Real tolerance, Real & norm) { AKANTU_DEBUG_IN(); UInt nb_nodes = residual->getSize(); norm = 0; Real * residual_val = residual->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); if(is_local_node) { for (UInt d = 0; d < spatial_dimension; ++d) { if(!(*boundary_val)) { norm += *residual_val * *residual_val; } boundary_val++; residual_val++; } } else { boundary_val += spatial_dimension; residual_val += spatial_dimension; } } StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum); norm = sqrt(norm); AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); return (norm < tolerance); } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceResidual(Real tolerance){ AKANTU_DEBUG_IN(); Real error=0; bool res = this->testConvergence<_scc_residual>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceResidual(Real tolerance, Real & error){ AKANTU_DEBUG_IN(); bool res = this->testConvergence<_scc_residual>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance){ AKANTU_DEBUG_IN(); Real error=0; bool res = this->testConvergence<_scc_increment>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance, Real & error){ AKANTU_DEBUG_IN(); bool res = this->testConvergence<_scc_increment>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::implicitPred() { AKANTU_DEBUG_IN(); if(previous_displacement) previous_displacement->copy(*displacement); if(method == _implicit_dynamic) integrator->integrationSchemePred(time_step, *displacement, *velocity, *acceleration, *boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::implicitCorr() { AKANTU_DEBUG_IN(); if(method == _implicit_dynamic) { integrator->integrationSchemeCorrDispl(time_step, *displacement, *velocity, *acceleration, *boundary, *increment); } else { UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes; Real * incr_val = increment->values; Real * disp_val = displacement->values; bool * boun_val = boundary->values; for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++disp_val, ++incr_val, ++boun_val){ *disp_val += (1. - *boun_val) * *incr_val; *incr_val *= (1. - *boun_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->values; Real * disp_val = displacement->values; Real * prev_disp_val = previous_displacement->values; 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(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::max(); updateCurrentPosition(); Element elem; elem.ghost_type = ghost_type; elem.kind = _ek_regular; Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator end = mesh.lastType(spatial_dimension, ghost_type); for(; it != end; ++it) { elem.type = *it; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*it); UInt nb_element = mesh.getNbElement(*it); Array::iterator< Vector > eibm = element_index_by_material(*it, ghost_type).begin(2); Array X(0, nb_nodes_per_element*spatial_dimension); FEM::extractNodalToElementField(mesh, *current_position, X, *it, _not_ghost); Array::iterator< Matrix > X_el = X.begin(spatial_dimension, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++X_el, ++eibm) { elem.element = el; Real el_size = getFEM().getElementInradius(*X_el, *it); Real el_dt = mat_val[(*eibm)(0)]->getStableTimeStep(el_size, elem); min_dt = std::min(min_dt, el_dt); } } AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real energy = 0.; /// call update residual on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { energy += (*mat_it)->getPotentialEnergy(); } /// reduction sum over all processors StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ 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->values; Real * mass_val = mass->values; for (UInt n = 0; n < nb_nodes; ++n) { Real mv2 = 0; bool is_local_node = mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !getIsPBCSlaveNode(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 = getFEM().getNbQuadraturePoints(type); Array vel_on_quad(nb_quadrature_points, spatial_dimension); Array filter_element(1, 1, index); getFEM().interpolateOnQuadraturePoints(*velocity, vel_on_quad, spatial_dimension, type, _not_ghost, filter_element); Array::iterator< Vector > vit = vel_on_quad.begin(spatial_dimension); Array::iterator< Vector > vend = vel_on_quad.end(spatial_dimension); Vector rho_v2(nb_quadrature_points); Real rho = materials[element_index_by_material(type)(index, 0)]->getRho(); for (UInt q = 0; vit != vend; ++vit, ++q) { rho_v2(q) = rho * vit->dot(*vit); } AKANTU_DEBUG_OUT(); return .5*getFEM().integrate(rho_v2, type, index); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getExternalWork() { AKANTU_DEBUG_IN(); Real * velo = velocity->storage(); Real * forc = force->storage(); Real * resi = residual->storage(); bool * boun = boundary->storage(); Real work = 0.; UInt nb_nodes = mesh.getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !getIsPBCSlaveNode(n); bool count_node = is_local_node && is_not_pbc_slave_node; for (UInt i = 0; i < spatial_dimension; ++i) { if (count_node) { if(*boun) work -= *resi * *velo * time_step; else work += *forc * *velo * time_step; } ++velo; ++forc; ++resi; ++boun; } } StaticCommunicator::getStaticCommunicator().allReduce(&work, 1, _so_sum); 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::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::iterator mat_it; Vector mat = element_index_by_material(type, _not_ghost).begin(2)[index]; Real energy = materials[mat(0)]->getEnergy(energy_id, type, mat(1)); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesAdded(const Array & 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(force ) force ->resize(nb_nodes); if(residual ) residual ->resize(nb_nodes); if(boundary ) boundary ->resize(nb_nodes); if(increment_acceleration) increment_acceleration->resize(nb_nodes); if(increment) increment->resize(nb_nodes); if(current_position) current_position->resize(nb_nodes); delete dof_synchronizer; dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onNodesAdded(nodes_list, event); } if (method != _explicit_lumped_mass) { delete stiffness_matrix; delete jacobian_matrix; delete solver; SolverOptions solver_options; initImplicit((method == _implicit_dynamic), solver_options); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsAdded(const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); getFEM().initShapeFunctions(_not_ghost); getFEM().initShapeFunctions(_ghost); Array::const_iterator it = element_list.begin(); Array::const_iterator end = element_list.end(); /// \todo have rules to choose the correct material UInt mat_id = 0; UInt * mat_id_vect = NULL; try { const NewMaterialElementsEvent & event_mat = dynamic_cast(event); mat_id_vect = event_mat.getMaterialList().storage(); } catch(...) { } for (UInt el = 0; it != end; ++it, ++el) { const Element & elem = *it; if(mat_id_vect) mat_id = mat_id_vect[el]; else mat_id = (*material_selector)(elem); Material & mat = *materials[mat_id]; UInt mat_index = mat.addElement(elem.type, elem.element, elem.ghost_type); Vector id(2); id[0] = mat_id; id[1] = mat_index; if(!element_index_by_material.exists(elem.type, elem.ghost_type)) element_index_by_material.alloc(0, 2, elem.type, elem.ghost_type); element_index_by_material(elem.type, elem.ghost_type).push_back(id); } std::vector::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) AKANTU_DEBUG_TO_IMPLEMENT(); assembleMassLumped(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsRemoved(__attribute__((unused)) const Array & element_list, const ByElementTypeUInt & new_numbering, const RemovedElementsEvent & event) { // MeshUtils::purifyMesh(mesh); getFEM().initShapeFunctions(_not_ghost); getFEM().initShapeFunctions(_ghost); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onElementsRemoved(element_list, new_numbering, event); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesRemoved(__attribute__((unused)) const Array & element_list, const Array & 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(force ) mesh.removeNodesFromArray(*force , new_numbering); if(residual ) mesh.removeNodesFromArray(*residual , new_numbering); if(boundary ) mesh.removeNodesFromArray(*boundary , new_numbering); if(increment_acceleration) mesh.removeNodesFromArray(*increment_acceleration, new_numbering); if(increment) mesh.removeNodesFromArray(*increment , new_numbering); delete dof_synchronizer; dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER #define ADD_FIELD(field, type) \ internalAddDumpFieldToDumper(dumper_name, \ BOOST_PP_STRINGIZE(field), \ new DumperIOHelper::NodalField(*field)) if(field_id == "displacement") { ADD_FIELD(displacement, Real); } else if(field_id == "mass" ) { ADD_FIELD(mass , Real); } else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); } else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); } else if(field_id == "force" ) { ADD_FIELD(force , Real); } else if(field_id == "residual" ) { ADD_FIELD(residual , Real); } else if(field_id == "boundary" ) { ADD_FIELD(boundary , bool); } else if(field_id == "increment" ) { ADD_FIELD(increment , Real); } else if(field_id == "partitions" ) { internalAddDumpFieldToDumper(dumper_name, field_id, new DumperIOHelper::ElementPartitionField<>(mesh, spatial_dimension, _not_ghost, _ek_regular)); } else if(field_id == "element_index_by_material") { internalAddDumpFieldToDumper(dumper_name, field_id, new DumperIOHelper::ElementalField(element_index_by_material, spatial_dimension, _not_ghost, _ek_regular)); } else { internalAddDumpFieldToDumper(dumper_name, field_id, new DumperIOHelper::HomogenizedField(*this, field_id, spatial_dimension, _not_ghost, _ek_regular)); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpGroupField(const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); this->addDumpGroupFieldToDumper(group.getDefaultDumperName(), field_id, group_name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name) { #ifdef AKANTU_USE_IOHELPER ElementGroup & group = mesh.getElementGroup(group_name); UInt dim = group.getDimension(); const Array * nodal_filter = &(group.getNodes()); const ByElementTypeArray * elemental_filter = &(group.getElements()); #define ADD_FIELD(field, type) \ group.addDumpFieldExternalToDumper(dumper_name, \ BOOST_PP_STRINGIZE(field), \ new DumperIOHelper::NodalField(*field, 0, 0, nodal_filter)) if(field_id == "displacement") { ADD_FIELD(displacement, Real); } else if(field_id == "mass" ) { ADD_FIELD(mass , Real); } else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); } else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); } else if(field_id == "force" ) { ADD_FIELD(force , Real); } else if(field_id == "residual" ) { ADD_FIELD(residual , Real); } else if(field_id == "boundary" ) { ADD_FIELD(boundary , bool); } else if(field_id == "increment" ) { ADD_FIELD(increment , Real); } else if(field_id == "partitions" ) { group.addDumpFieldExternalToDumper(dumper_name, field_id, new DumperIOHelper::ElementPartitionField(mesh, dim, _not_ghost, _ek_regular, elemental_filter)); } else if(field_id == "element_index_by_material") { group.addDumpFieldExternalToDumper(dumper_name, field_id, new DumperIOHelper::ElementalField(element_index_by_material, dim, _not_ghost, _ek_regular, elemental_filter)); } else { typedef DumperIOHelper::HomogenizedField Field; group.addDumpFieldExternalToDumper(dumper_name, field_id, new Field(*this, field_id, dim, _not_ghost, _ek_regular, elemental_filter)); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::removeDumpGroupField(const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); this->removeDumpGroupFieldFromDumper(group.getDefaultDumperName(), field_id, group_name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::removeDumpGroupFieldFromDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.removeDumpFieldFromDumper(dumper_name, field_id); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER #define ADD_FIELD(field, type) \ DumperIOHelper::Field * f = \ new DumperIOHelper::NodalField(*field); \ f->setPadding(3); \ internalAddDumpFieldToDumper(dumper_name, BOOST_PP_STRINGIZE(field), f) if(field_id == "displacement") { ADD_FIELD(displacement, Real); } else if(field_id == "mass" ) { ADD_FIELD(mass , Real); } else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); } else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); } else if(field_id == "force" ) { ADD_FIELD(force , Real); } else if(field_id == "residual" ) { ADD_FIELD(residual , Real); } else { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular); field->setPadding(3); internalAddDumpFieldToDumper(dumper_name, field_id, field); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpGroupFieldVector(const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); this->addDumpGroupFieldVectorToDumper(group.getDefaultDumperName(), field_id, group_name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpGroupFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name) { #ifdef AKANTU_USE_IOHELPER ElementGroup & group = mesh.getElementGroup(group_name); UInt dim = group.getDimension(); const Array * nodal_filter = &(group.getNodes()); const ByElementTypeArray * elemental_filter = &(group.getElements()); #define ADD_FIELD(field, type) \ DumperIOHelper::Field * f = \ new DumperIOHelper::NodalField(*field, 0, 0, nodal_filter); \ f->setPadding(3); \ group.addDumpFieldExternalToDumper(dumper_name, BOOST_PP_STRINGIZE(field), f) if(field_id == "displacement") { ADD_FIELD(displacement, Real); } else if(field_id == "mass" ) { ADD_FIELD(mass , Real); } else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); } else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); } else if(field_id == "force" ) { ADD_FIELD(force , Real); } else if(field_id == "residual" ) { ADD_FIELD(residual , Real); } else { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, dim, _not_ghost, _ek_regular, elemental_filter); field->setPadding(3); group.addDumpFieldExternalToDumper(dumper_name, field_id, field); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpFieldTensorToDumper(const std::string & dumper_name, const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER if(field_id == "stress") { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular); field->setPadding(3, 3); internalAddDumpFieldToDumper(dumper_name, field_id, field); } else if(field_id == "strain") { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular); field->setPadding(3, 3); internalAddDumpFieldToDumper(dumper_name, field_id, field); } else { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular); field->setPadding(3, 3); internalAddDumpFieldToDumper(dumper_name, field_id, field); } #endif } /* -------------------------------------------------------------------------- */ 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; getFEM().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); force ->printself(stream, indent + 2); residual ->printself(stream, indent + 2); boundary ->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + connectivity type information [" << std::endl; element_index_by_material.printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index ba209e97a..d58704a4c 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,673 +1,675 @@ /** * @file solid_mechanics_model.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date Tue Jul 27 18:15:37 2010 * * @brief Model of Solid Mechanics * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_HH__ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_types.hh" #include "model.hh" #include "data_accessor.hh" #include "mesh.hh" #include "dumpable.hh" #include "boundary_condition.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #include "integration_scheme_2nd_order.hh" #include "solver.hh" #include "material_selector.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class Material; class IntegrationScheme2ndOrder; class Contact; class SparseMatrix; } __BEGIN_AKANTU__ struct SolidMechanicsModelOptions : public ModelOptions { SolidMechanicsModelOptions(AnalysisMethod analysis_method = _explicit_lumped_mass, bool no_init_materials = false) : analysis_method(analysis_method), no_init_materials(no_init_materials) {} AnalysisMethod analysis_method; bool no_init_materials; }; extern const SolidMechanicsModelOptions default_solid_mechanics_model_options; class SolidMechanicsModel : public Model, public DataAccessor, public MeshEventHandler, public Dumpable, public BoundaryCondition { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: class NewMaterialElementsEvent : public NewElementsEvent { public: AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Array &); AKANTU_GET_MACRO(MaterialList, material, const Array &); protected: Array material; }; typedef FEMTemplate MyFEMType; SolidMechanicsModel(Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "solid_mechanics_model", const MemoryID & memory_id = 0); virtual ~SolidMechanicsModel(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize completely the model virtual void initFull(std::string material_file, const ModelOptions & options = default_solid_mechanics_model_options); /// initialize the fem object needed for boundary conditions void initFEMBoundary(); /// register the tags associated with the parallel synchronizer void initParallel(MeshPartition * partition, DataAccessor * data_accessor=NULL); /// allocate all vectors void initArrays(); /// allocate all vectors void initArraysPreviousDisplacment(); /// initialize all internal arrays for materials virtual void initMaterials(); /// initialize the model virtual void initModel(); /// init PBC synchronizer void initPBC(); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* PBC */ /* ------------------------------------------------------------------------ */ public: /// change the equation number for proper assembly when using PBC void changeEquationNumberforPBC(std::map & pbc_pair); /// synchronize Residual for output void synchronizeResidual(); protected: /// register PBC synchronizer void registerPBCSynchronizer(); /* ------------------------------------------------------------------------ */ /* Explicit */ /* ------------------------------------------------------------------------ */ public: /// initialize the stuff for the explicit scheme void initExplicit(AnalysisMethod analysis_method = _explicit_lumped_mass); bool isExplicit() { return method == _explicit_lumped_mass || method == _explicit_consistent_mass; } /// initialize the array needed by updateResidual (residual, current_position) void initializeUpdateResidualData(); /// update the current position vector void updateCurrentPosition(); /// assemble the residual for the explicit scheme virtual void updateResidual(bool need_initialize = true); /** * \brief compute the acceleration from the residual * this function is the explicit equivalent to solveDynamic in implicit * In the case of lumped mass just divide the residual by the mass * In the case of not lumped mass call solveDynamic<_acceleration_corrector> */ void updateAcceleration(); void updateIncrement(); void updatePreviousDisplacement(); /// Solve the system @f[ A x = \alpha b @f] with A a lumped matrix void solveLumped(Array & x, const Array & A, const Array & b, const Array & boundary, Real alpha); /// explicit integration predictor void explicitPred(); /// explicit integration corrector void explicitCorr(); public: void solveStep(); /* ------------------------------------------------------------------------ */ /* Implicit */ /* ------------------------------------------------------------------------ */ public: /// initialize the solver and the jacobian_matrix (called by initImplicit) void initSolver(SolverOptions & options = _solver_no_options); /// initialize the stuff for the implicit solver void initImplicit(bool dynamic = false, SolverOptions & solver_options = _solver_no_options); /// solve Ma = f to get the initial acceleration void initialAcceleration(); /// assemble the stiffness matrix void assembleStiffnessMatrix(); public: template bool solveStep(Real tolerance, UInt max_iteration = 100); public: /// solve @f[ A\delta u = f_ext - f_int @f] in displacement void solveDynamic(); /// solve Ku = f void solveStatic(); /// solve Ku = f void solveStatic(Array & boundary_normal, Array & EulerAngles); /// test if the system is converged template bool testConvergence(Real tolerance, Real & error); /// test the convergence (norm of increment) bool testConvergenceIncrement(Real tolerance) __attribute__((deprecated)); bool testConvergenceIncrement(Real tolerance, Real & error) __attribute__((deprecated)); /// test the convergence (norm of residual) bool testConvergenceResidual(Real tolerance) __attribute__((deprecated)); bool testConvergenceResidual(Real tolerance, Real & error) __attribute__((deprecated)); /// create and return the velocity damping matrix SparseMatrix & initVelocityDampingMatrix(); /// implicit time integration predictor void implicitPred(); /// implicit time integration corrector void implicitCorr(); protected: /// finish the computation of residual to solve in increment void updateResidualInternal(); /// compute the support reaction and store it in force void updateSupportReaction(); public: //protected: Daniel changed it just for a test /// compute A and solve @f[ A\delta u = f_ext - f_int @f] template void solve(Array & increment); /* ------------------------------------------------------------------------ */ /* Explicit/Implicit */ /* ------------------------------------------------------------------------ */ public: /// Update the stresses for the computation of the residual of the Stiffness /// matrix in the case of finite deformation void computeStresses(); /// synchronize the ghost element boundaries values void synchronizeBoundaries(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: - /// registers all the custom materials of a given types present in the input file + /// registers all the custom materials of a given type present in the input file template void registerNewCustomMaterials(const ID & mat_type); /// register an empty material of a given type template - Material & registerNewEmptyMaterial(const ID & mat_type, const std::string & name); + Material & registerNewEmptyMaterial(const std::string & name); // /// Use a UIntData in the mesh to specify the material to use per element // void setMaterialIDsFromIntData(const std::string & data_name); protected: /// register a material in the dynamic database template Material & registerNewMaterial(const ParserSection & mat_section); /// read the material files to instantiate all the materials void instantiateMaterials(); /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ 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); /// fill a vector of rho void computeRho(Array & rho, ElementType type, GhostType ghost_type); /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: inline virtual UInt getNbDataForElements(const Array & elements, SynchronizationTag tag) const; inline virtual void packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const; inline virtual void unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag); inline virtual UInt getNbDataToPack(SynchronizationTag tag) const; inline virtual UInt getNbDataToUnpack(SynchronizationTag tag) const; inline virtual void packData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag) const; inline virtual void unpackData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag); protected: inline void splitElementByMaterial(const Array & elements, Array * elements_per_mat) const; inline virtual void packBarycenter(const Mesh & mesh, CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const; inline virtual void unpackBarycenter(const Mesh & mesh, CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: virtual void onNodesAdded (const Array & nodes_list, const NewNodesEvent & event); virtual void onNodesRemoved(const Array & element_list, const Array & new_numbering, const RemovedNodesEvent & event); virtual void onElementsAdded (const Array & nodes_list, const NewElementsEvent & event); virtual void onElementsRemoved(const Array & element_list, const ByElementTypeUInt & new_numbering, const RemovedElementsEvent & event); /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); virtual void addDumpGroupField(const std::string & field_id, const std::string & group_name); virtual void addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name); virtual void removeDumpGroupField(const std::string & field_id, const std::string & group_name); virtual void removeDumpGroupFieldFromDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name); virtual void addDumpFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id); virtual void addDumpGroupFieldVector(const std::string & field_id, const std::string & group_name); virtual void addDumpGroupFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name); virtual void addDumpFieldTensorToDumper(const std::string & dumper_name, const std::string & field_id); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the current value of the time step AKANTU_GET_MACRO(TimeStep, time_step, Real); /// set the value of the time step void setTimeStep(Real time_step); /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement, Array &); /// get the SolidMechanicsModel::previous_displacement vector AKANTU_GET_MACRO(PreviousDisplacement, *previous_displacement, Array &); /// get the SolidMechanicsModel::current_position vector \warn only consistent /// after a call to SolidMechanicsModel::updateCurrentPosition AKANTU_GET_MACRO(CurrentPosition, *current_position, const Array &); /// get the SolidMechanicsModel::increment vector \warn only consistent if /// SolidMechanicsModel::setIncrementFlagOn has been called before AKANTU_GET_MACRO(Increment, *increment, Array &); /// get the lumped SolidMechanicsModel::mass vector AKANTU_GET_MACRO(Mass, *mass, Array &); /// get the SolidMechanicsModel::velocity vector AKANTU_GET_MACRO(Velocity, *velocity, Array &); /// get the SolidMechanicsModel::acceleration vector, updated by /// SolidMechanicsModel::updateAcceleration AKANTU_GET_MACRO(Acceleration, *acceleration, Array &); /// get the SolidMechanicsModel::force vector (boundary forces) AKANTU_GET_MACRO(Force, *force, Array &); /// get the SolidMechanicsModel::residual vector, computed by /// SolidMechanicsModel::updateResidual AKANTU_GET_MACRO(Residual, *residual, Array &); /// get the SolidMechanicsModel::boundary vector AKANTU_GET_MACRO(Boundary, *boundary, Array &); /// get the SolidMechnicsModel::incrementAcceleration vector AKANTU_GET_MACRO(IncrementAcceleration, *increment_acceleration, Array &); /// get the value of the SolidMechanicsModel::increment_flag AKANTU_GET_MACRO(IncrementFlag, increment_flag, bool); /// get a particular material (by material index) inline Material & getMaterial(UInt mat_index); /// get a particular material (by material index) inline const Material & getMaterial(UInt mat_index) const; /// get a particular material (by material name) inline Material & getMaterial(const std::string & name); /// get a particular material (by material name) inline const Material & getMaterial(const std::string & name) const; /// get a particular material id from is name inline UInt getMaterialIndex(const std::string & name) const; /// give the number of materials inline UInt getNbMaterials() const { return materials.size(); }; inline void setMaterialSelector(MaterialSelector & selector); /// give the material internal index from its id Int getInternalIndexFromID(const ID & id) const; /// compute the stable time step Real getStableTimeStep(); /// compute the potential energy Real getPotentialEnergy(); /// compute the kinetic energy Real getKineticEnergy(); Real getKineticEnergy(const ElementType & type, UInt index); /// compute the external work (for impose displacement, the velocity should be given too) Real getExternalWork(); /// get the energies Real getEnergy(const std::string & energy_id); /// compute the energy for energy Real getEnergy(const std::string & energy_id, const ElementType & type, UInt index); /// set the Contact object AKANTU_SET_MACRO(Contact, contact, Contact *); /** * @brief set the SolidMechanicsModel::increment_flag to on, the activate the * update of the SolidMechanicsModel::increment vector */ void setIncrementFlagOn(); /// get the stiffness matrix AKANTU_GET_MACRO(StiffnessMatrix, *stiffness_matrix, SparseMatrix &); /// get the mass matrix AKANTU_GET_MACRO(MassMatrix, *mass_matrix, SparseMatrix &); /// get the velocity damping matrix AKANTU_GET_MACRO(VelocityDampingMatrix, *velocity_damping_matrix, SparseMatrix &); /// get the FEM object to integrate or interpolate on the boundary inline FEM & getFEMBoundary(const ID & name = ""); /// get integrator AKANTU_GET_MACRO(Integrator, *integrator, const IntegrationScheme2ndOrder &); /// get access to the internal solver AKANTU_GET_MACRO(Solver, *solver, Solver &); /// get synchronizer AKANTU_GET_MACRO(Synchronizer, *synch_parallel, const DistributedSynchronizer &); AKANTU_GET_MACRO(ElementIndexByMaterial, element_index_by_material, const ByElementTypeArray &); -protected: + /// vectors containing local material element index for each global element index AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementIndexByMaterial, element_index_by_material, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementIndexByMaterial, element_index_by_material, UInt); +protected: + friend Material; template class WeightFunction> friend class MaterialNonLocal; protected: /// compute the stable time step Real getStableTimeStep(const GhostType & ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// time step Real time_step; /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array Array * displacement; /// displacements array at the previous time step (used in finite deformation) Array * previous_displacement; /// lumped mass array Array * mass; /// velocities array Array * velocity; /// accelerations array Array * acceleration; /// accelerations array Array * increment_acceleration; /// forces array Array * force; /// residuals array Array * residual; /// boundaries array Array * boundary; /// array of current position used during update residual Array * current_position; /// mass matrix SparseMatrix * mass_matrix; /// velocity damping matrix SparseMatrix * velocity_damping_matrix; /// stiffness matrix SparseMatrix * stiffness_matrix; /// jacobian matrix @f[A = cM + dD + K@f] with @f[c = \frac{1}{\beta \Delta /// t^2}, d = \frac{\gamma}{\beta \Delta t} @f] SparseMatrix * jacobian_matrix; /// vectors containing local material element index for each global element index ByElementTypeUInt element_index_by_material; /// list of used materials std::vector materials; /// mapping between material name and material internal id std::map materials_names_to_id; /// class defining of to choose a material MaterialSelector * material_selector; /// define if it is the default selector or not bool is_default_material_selector; /// integration scheme of second order used IntegrationScheme2ndOrder * integrator; /// increment of displacement Array * increment; /// flag defining if the increment must be computed or not bool increment_flag; /// solver for implicit Solver * solver; /// object to resolve the contact Contact * contact; /// analysis method check the list in akantu::AnalysisMethod AnalysisMethod method; /// internal synchronizer for parallel computations DistributedSynchronizer * synch_parallel; /// tells if the material are instantiated bool are_materials_instantiated; }; __END_AKANTU__ /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "parser.hh" #include "material.hh" __BEGIN_AKANTU__ #include "solid_mechanics_model_tmpl.hh" #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "solid_mechanics_model_inline_impl.cc" #endif /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const SolidMechanicsModel & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #include "material_selector_tmpl.hh" #endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_tmpl.hh b/src/model/solid_mechanics/solid_mechanics_model_tmpl.hh index 0ff2456b0..baf20941a 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_tmpl.hh +++ b/src/model/solid_mechanics/solid_mechanics_model_tmpl.hh @@ -1,104 +1,103 @@ /** * @file solid_mechanics_model_tmpl.hh * * @author Guillaume Anciaux * @author Nicolas Richart * @author Dana Christen * * @date Thu Nov 24 09:36:33 2011 * * @brief template part of solid mechanics model * * @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 . * */ /* -------------------------------------------------------------------------- */ template Material & SolidMechanicsModel::registerNewMaterial(const ParserSection & section) { std::string mat_name; std::string mat_type = section.getName(); try { std::string tmp = section.getParameter("name"); mat_name = tmp; /** this can seam weird, but there is an ambiguous operator * overload that i couldn't solve. @todo remove the * weirdness of this code */ } catch(debug::Exception & ){ AKANTU_DEBUG_ERROR("A material of type \'" << mat_type << "\' in the input file has been defined without a name!"); } AKANTU_DEBUG_ASSERT(materials_names_to_id.find(mat_name) == materials_names_to_id.end(), "A material with this name '" << mat_name << "' has already been registered. " << "Please use unique names for materials"); UInt mat_count = materials.size(); materials_names_to_id[mat_name] = mat_count; std::stringstream sstr_mat; sstr_mat << id << ":" << mat_count << ":" << mat_type; ID mat_id = sstr_mat.str(); Material * material; material = new M(*this, mat_id); materials.push_back(material); material->parseSection(section); return *material; } /* -------------------------------------------------------------------------- */ template -Material & SolidMechanicsModel::registerNewEmptyMaterial(const ID & mat_type, - const std::string & mat_name) { +Material & SolidMechanicsModel::registerNewEmptyMaterial(const std::string & mat_name) { AKANTU_DEBUG_ASSERT(materials_names_to_id.find(mat_name) == materials_names_to_id.end(), "A material with this name '" << mat_name << "' has already been registered. " << "Please use unique names for materials"); UInt mat_count = materials.size(); materials_names_to_id[mat_name] = mat_count; - std::stringstream sstr_mat; sstr_mat << id << ":" << mat_count << ":" << mat_type; + std::stringstream sstr_mat; sstr_mat << id << ":" << mat_count << ":" << mat_name; ID mat_id = sstr_mat.str(); Material * material; material = new M(*this, mat_id); materials.push_back(material); return *material; } /* -------------------------------------------------------------------------- */ template void SolidMechanicsModel::registerNewCustomMaterials(const ID & mat_type) { std::pair sub_sect = parser.getSubSections(_st_material); Parser::const_section_iterator it = sub_sect.first; for (; it != sub_sect.second; ++it) { if(it->getName() == mat_type) { registerNewMaterial(*it); } } }