diff --git a/CMakeLists.txt b/CMakeLists.txt index b3af4ebb1..63aaa4987 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,175 +1,175 @@ #=============================================================================== # @file CMakeLists.txt # # @author Nicolas Richart # # @date Mon Jun 14 19:12:20 2010 # # @brief main configuration file # # @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 . # # @section DESCRIPTION #------------------------------------------------------------------------------- # _ _ # | | | | # __ _| | ____ _ _ __ | |_ _ _ # / _` | |/ / _` | '_ \| __| | | | # | (_| | < (_| | | | | |_| |_| | # \__,_|_|\_\__,_|_| |_|\__|\__,_| # #=============================================================================== #=============================================================================== # CMake Project #=============================================================================== cmake_minimum_required(VERSION 2.8) project(akantu) enable_language(CXX) #=============================================================================== # Misc. config for cmake #=============================================================================== set(AKANTU_CMAKE_DIR "${PROJECT_SOURCE_DIR}/cmake") set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/cmake") set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/cmake/Modules") set(BUILD_SHARED_LIBS ON CACHE BOOL "Build shared libraries.") if(NOT AKANTU_TARGETS_EXPORT) set(AKANTU_TARGETS_EXPORT AkantuLibraryDepends) endif() include(CMakeVersionGenerator) include(CMakePackagesSystem) include(CMakeFlagsHandling) + include(AkantuMacros) #cmake_activate_debug_message() #=============================================================================== # Version Number #=============================================================================== # AKANTU version number. An even minor number corresponds to releases. set(AKANTU_MAJOR_VERSION 1) set(AKANTU_MINOR_VERSION 0) set(AKANTU_PATCH_VERSION 0) define_project_version() #=============================================================================== # Options #=============================================================================== # Debug option(AKANTU_DEBUG "Compiles akantu with the debug messages" ON) mark_as_advanced(AKANTU_DEBUG) add_flags(cxx "-Wall") if(NOT AKANTU_DEBUG) set(AKANTU_NDEBUG ON) else() set(AKANTU_NDEBUG OFF) endif() #Profiling set(CMAKE_CXX_FLAGS_PROFILING "-g -pg" CACHE STRING "Flags used by the compiler during profiling builds") set(CMAKE_EXE_LINKER_FLAGS_PROFILING "-pg" CACHE STRING "Flags used by the linker during profiling builds") set(CMAKE_SHARED_LINKER_FLAGS_PROFILING "-pg" CACHE STRING "Flags used by the linker during profiling builds") mark_as_advanced(CMAKE_CXX_FLAGS_PROFILING) mark_as_advanced(CMAKE_EXE_LINKER_FLAGS_PROFILING) mark_as_advanced(CMAKE_SHARED_LINKER_FLAGS_PROFILING) include(CMakeDetermineCCompiler) #=============================================================================== # Dependencies #=============================================================================== set(Boost_USE_STATIC_LIBS ON) set(Boost_USE_MULTITHREADED ON) set(Boost_USE_STATIC_RUNTIME OFF) find_package( Boost 1.36.0 COMPONENTS chrono system) if(Boost_FOUND) include_directories(${Boost_INCLUDE_DIRS}) set (AKANTU_EXTERNAL_LIBRARIES ${Boost_LIBRARIES}) endif() add_all_packages(${PROJECT_SOURCE_DIR}/packages ${PROJECT_SOURCE_DIR}/src) ## meta option \todo better way to do it when multiple package give enable the ## same feature - if(AKANTU_SCOTCH) set(AKANTU_PARTITIONER ON) else() set(AKANTU_PARTITIONER OFF) endif() if(AKANTU_MUMPS) set(AKANTU_SOLVER ON) else() set(AKANTU_SOLVER OFF) endif() #=============================================================================== # Akantu library #=============================================================================== add_subdirectory(src) #=============================================================================== # Documentation #=============================================================================== option(AKANTU_DOCUMENTATION "Build source documentation using Doxygen." OFF) if(AKANTU_DOCUMENTATION) add_subdirectory(doc) endif() #=============================================================================== # Examples and tests #=============================================================================== option(AKANTU_EXAMPLES "Activate examples" OFF) option(AKANTU_TESTS "Activate tests" OFF) if(AKANTU_EXAMPLES OR AKANTU_TESTS) include(AkantuTestAndExamples) find_package(GMSH REQUIRED) endif() if(AKANTU_EXAMPLES) add_subdirectory(examples) endif() if(AKANTU_TESTS) enable_testing() include(CTest) set(AKANTU_TESTS_EXCLUDE_FILES) add_subdirectory(test) endif() #=============================================================================== # Install and Packaging #=============================================================================== include(AkantuInstall) if(NOT AKANTU_DISABLE_CPACK) message("AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA") include(AkantuCPack) endif() set(TMP ${AKANTU_INCLUDE_DIRS}) set(AKANTU_INCLUDE_DIRS ${TMP} CACHE INTERNAL "Internal include directories to link with Akantu as a subproject") diff --git a/src/model/heat_transfer/heat_transfer_model.hh b/src/model/heat_transfer/heat_transfer_model.hh index 05e8f18da..7f74948f5 100644 --- a/src/model/heat_transfer/heat_transfer_model.hh +++ b/src/model/heat_transfer/heat_transfer_model.hh @@ -1,306 +1,306 @@ /** * @file heat_transfer_model.hh * * @author Rui Wang * @author Srinivasa Babu Ramisetti * @author Guillaume Anciaux * * @date Sun May 01 19:14:43 2011 * * @brief Model of Heat Transfer * * @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_HEAT_TRANSFER_MODEL_HH__ #define __AKANTU_HEAT_TRANSFER_MODEL_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "model.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" namespace akantu { class IntegrationScheme1stOrder; // class Solver; // class SparseMatrix; } __BEGIN_AKANTU__ class HeatTransferModel : public Model, public DataAccessor { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef FEMTemplate MyFEMType; HeatTransferModel(UInt spatial_dimension, const ID & id = "heat_transfer_model", const MemoryID & memory_id = 0) ; HeatTransferModel(Mesh & mesh, UInt spatial_dimension = 0, const ID & id = "heat_transfer_model", const MemoryID & memory_id = 0); virtual ~HeatTransferModel() ; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// generic function to initialize everything ready for explicit dynamics void initFull(const std::string & material_file); /// initialize the fem object of the boundary void initFEMBoundary(bool create_surface = true); /// set the parameters bool setParam(const std::string & key, const std::string & value); /// read one material file to instantiate all the materials void readMaterials(const std::string & filename); /// allocate all vectors void initVectors(); /// register the tags associated with the parallel synchronizer void initParallel(MeshPartition * partition, DataAccessor * data_accessor=NULL); /// initialize the model void initModel(); /// init PBC synchronizer void initPBC(); /// function to print the contain of the class virtual void printself(__attribute__ ((unused)) std::ostream & stream, __attribute__ ((unused)) int indent = 0) const {}; /* ------------------------------------------------------------------------ */ /* Methods for explicit */ /* ------------------------------------------------------------------------ */ public: /// compute and get the stable time step Real getStableTimeStep(); /// compute the heat flux void updateResidual(); /// calculate the lumped capacity vector for heat transfer problem void assembleCapacityLumped(); /// update the temperature from the temperature rate void explicitPred(); /// update the temperature rate from the increment void explicitCorr(); // /// initialize the heat flux // void initializeResidual(Vector &temp); // /// initialize temperature // void initializeTemperature(Vector &temp); private: /// solve the system in temperature rate @f$C\delta \dot T = q_{n+1} - C \dot T_{n}@f$ void solveExplicitLumped(); /// compute the heat flux on ghost types void updateResidual(const GhostType & ghost_type); /// calculate the lumped capacity vector for heat transfer problem (w ghosttype) void assembleCapacityLumped(const GhostType & ghost_type); /// compute the conductivity tensor for each quadrature point in an array void computeConductivityOnQuadPoints(const GhostType & ghost_type); /// compute vector k \grad T for each quadrature point void computeKgradT(const GhostType & ghost_type); /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: inline UInt getNbDataForElements(const Vector & elements, SynchronizationTag tag) const; - inline void packData(CommunicationBuffer & buffer, - const Vector & elements, - SynchronizationTag tag) const; - inline void unpackData(CommunicationBuffer & buffer, - const Vector & elements, - SynchronizationTag tag); + inline void packElementData(CommunicationBuffer & buffer, + const Vector & elements, + SynchronizationTag tag) const; + inline void unpackElementData(CommunicationBuffer & buffer, + const Vector & elements, + SynchronizationTag tag); inline UInt getNbDataToPack(SynchronizationTag tag) const; inline UInt getNbDataToUnpack(SynchronizationTag tag) const; inline void packData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag) const; inline void unpackData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: inline FEM & getFEMBoundary(std::string name = ""); /// get 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 AKANTU_SET_MACRO(TimeStep, time_step, Real); /// get the assembled heat flux AKANTU_GET_MACRO(Residual, *residual, Vector&); /// get the lumped capacity AKANTU_GET_MACRO(CapacityLumped, * capacity_lumped, Vector&); /// get the boundary vector AKANTU_GET_MACRO(Boundary, * boundary, Vector&); /// get the external flux vector AKANTU_GET_MACRO(ExternalFlux, * external_flux, Vector&); /// get the temperature gradient AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(TemperatureGradient, temperature_gradient, Real); /// get the conductivity on q points AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ConductivityOnQpoints, conductivity_on_qpoints, Real); /// get the conductivity on q points AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(TemperatureOnQpoints, temperature_on_qpoints, Real); /// internal variables AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(KGradtOnQpoints, k_gradt_on_qpoints, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(IntBtKgT, int_bt_k_gT, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(BtKgT, bt_k_gT, Real); /// get the temperature AKANTU_GET_MACRO(Temperature, *temperature, Vector &); /// get the temperature derivative AKANTU_GET_MACRO(TemperatureRate, *temperature_rate, Vector &); /// get the equation number Vector AKANTU_GET_MACRO(EquationNumber, *equation_number, const Vector &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: IntegrationScheme1stOrder * integrator; /// time step Real time_step; /// temperatures array Vector * temperature; /// temperatures derivatives array Vector * temperature_rate; /// increment array (@f$\delta \dot T@f$ or @f$\delta T@f$) Vector * increment; /// the spatial dimension UInt spatial_dimension; /// the density Real density; /// the speed of the changing temperature ByElementTypeReal temperature_gradient; /// temperature field on quadrature points ByElementTypeReal temperature_on_qpoints; /// conductivity tensor on quadrature points ByElementTypeReal conductivity_on_qpoints; /// vector k \grad T on quad points ByElementTypeReal k_gradt_on_qpoints; /// vector \int \grad N k \grad T ByElementTypeReal int_bt_k_gT; /// vector \grad N k \grad T ByElementTypeReal bt_k_gT; //external flux vector Vector * external_flux; /// residuals array Vector * residual; /// position of a dof in the K matrix Vector * equation_number; //lumped vector Vector * capacity_lumped; /// boundary vector Vector * boundary; //realtime Real time; ///capacity Real capacity; //conductivity matrix Real* conductivity; //linear variation of the conductivity (for temperature dependent conductivity) Real conductivity_variation; // reference temperature for the interpretation of temperature variation Real t_ref; //the biggest parameter of conductivity matrix Real conductivitymax; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "heat_transfer_model_inline_impl.cc" #endif /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const HeatTransferModel & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_HEAT_TRANSFER_MODEL_HH__ */ diff --git a/src/model/heat_transfer/heat_transfer_model_inline_impl.cc b/src/model/heat_transfer/heat_transfer_model_inline_impl.cc index 004e337b7..6b36f0da7 100644 --- a/src/model/heat_transfer/heat_transfer_model_inline_impl.cc +++ b/src/model/heat_transfer/heat_transfer_model_inline_impl.cc @@ -1,299 +1,298 @@ /** * @file heat_transfer_model_inline_impl.cc * * @author Guillaume Anciaux * @author Srinivasa Babu Ramisetti * * @date Sun May 01 19:14:43 2011 * * @brief Implementation of the inline functions of the HeatTransferModel 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 . * */ /* -------------------------------------------------------------------------- */ inline FEM & HeatTransferModel::getFEMBoundary(std::string name) { return dynamic_cast(getFEMClassBoundary(name)); } /* -------------------------------------------------------------------------- */ inline UInt HeatTransferModel::getNbDataToPack(SynchronizationTag tag) const{ AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes = getFEM().getMesh().getNbNodes(); switch(tag) { case _gst_htm_temperature: case _gst_htm_capacity: { size += nb_nodes * sizeof(Real); break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ inline UInt HeatTransferModel::getNbDataToUnpack(SynchronizationTag tag) const{ AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes = getFEM().getMesh().getNbNodes(); switch(tag) { case _gst_htm_capacity: case _gst_htm_temperature: { size += nb_nodes * sizeof(Real); break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ inline void HeatTransferModel::packData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag) const{ AKANTU_DEBUG_IN(); switch(tag) { case _gst_htm_capacity: buffer << (*capacity_lumped)(index); break; case _gst_htm_temperature: { buffer << (*temperature)(index); break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void HeatTransferModel::unpackData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag) { AKANTU_DEBUG_IN(); switch(tag) { case _gst_htm_capacity: { buffer >> (*capacity_lumped)(index); break; } case _gst_htm_temperature: { buffer >> (*temperature)(index); break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline UInt HeatTransferModel::getNbDataForElements(const Vector & elements, SynchronizationTag tag) const { AKANTU_DEBUG_IN(); - UInt size = 0; UInt nb_nodes_per_element = 0; Vector::const_iterator it = elements.begin(); Vector::const_iterator end = elements.end(); for (; it != end; ++it) { const Element & el = *it; nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); } #ifndef AKANTU_NDEBUG size += elements.getSize() * spatial_dimension * sizeof(Real); /// position of the barycenter of the element (only for check) - size += spatial_dimension * nb_nodes_per_element * sizeof(Real); /// position of the nodes of the element + // size += spatial_dimension * nb_nodes_per_element * sizeof(Real); /// position of the nodes of the element #endif switch(tag) { case _gst_htm_capacity: { size += nb_nodes_per_element * sizeof(Real); // capacity vector break; } case _gst_htm_temperature: { size += nb_nodes_per_element * sizeof(Real); // temperature break; } case _gst_htm_gradient_temperature: { - size += spatial_dimension * sizeof(Real); // temperature gradient + size += getNbQuadraturePoints(elements) * spatial_dimension * sizeof(Real); // temperature gradient size += nb_nodes_per_element * sizeof(Real); // nodal temperatures - size += spatial_dimension * nb_nodes_per_element * sizeof(Real); // shape derivatives + // size += spatial_dimension * nb_nodes_per_element * sizeof(Real); // shape derivatives break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ -inline void HeatTransferModel::packData(CommunicationBuffer & buffer, - const Vector & elements, - SynchronizationTag tag) const { +inline void HeatTransferModel::packElementData(CommunicationBuffer & buffer, + const Vector & elements, + SynchronizationTag tag) const { #ifndef AKANTU_NDEBUG Vector::const_iterator bit = elements.begin(); Vector::const_iterator bend = elements.end(); for (; bit != bend; ++bit) { const Element & element = *bit; types::RVector barycenter(spatial_dimension); mesh.getBarycenter(element.element, element.type, barycenter.storage(), element.ghost_type); buffer << barycenter; } // packNodalDataHelper(mesh.getNodes(), buffer, elements); #endif switch (tag){ case _gst_htm_capacity: { packNodalDataHelper(*capacity_lumped, buffer, elements); break; } case _gst_htm_temperature: { packNodalDataHelper(*temperature, buffer, elements); break; } case _gst_htm_gradient_temperature: { packElementalDataHelper(temperature_gradient, buffer, elements); packNodalDataHelper(*temperature, buffer, elements); // Vector::const_iterator it_shaped = // getFEM().getShapesDerivatives(element.type, ghost_type).begin(nb_nodes_per_element,spatial_dimension); // buffer << it_shaped[element.element]; break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } } /* -------------------------------------------------------------------------- */ -inline void HeatTransferModel::unpackData(CommunicationBuffer & buffer, - const Vector & elements, - SynchronizationTag tag) { +inline void HeatTransferModel::unpackElementData(CommunicationBuffer & buffer, + const Vector & elements, + SynchronizationTag tag) { #ifndef AKANTU_NDEBUG Vector::const_iterator bit = elements.begin(); Vector::const_iterator bend = elements.end(); for (; bit != bend; ++bit) { const Element & element = *bit; types::RVector barycenter_loc(spatial_dimension); mesh.getBarycenter(element.element, element.type, barycenter_loc.storage(), element.ghost_type); types::RVector barycenter(spatial_dimension); buffer >> barycenter; Real tolerance = 1e-15; for (UInt i = 0; i < spatial_dimension; ++i) { if(!(std::abs(barycenter(i) - barycenter_loc(i)) <= tolerance)) AKANTU_DEBUG_ERROR("Unpacking an unknown value for the element: " << element << "(barycenter[" << i << "] = " << barycenter_loc(i) << " and buffer[" << i << "] = " << barycenter(i) << ") - tag: " << tag); } } // types::RVector coords(spatial_dimension); // Real * nodes = getFEM().getMesh().getNodes().values; // for (UInt n = 0; n < nb_nodes_per_element; ++n) { // buffer >> coords; // UInt offset_conn = conn[el_offset + n]; // Real * coords_local = nodes+spatial_dimension*offset_conn; // for (UInt i = 0; i < spatial_dimension; ++i) { // if(!(std::abs(coords(i) - coords_local[i]) <= tolerance)) // AKANTU_EXCEPTION("Unpacking to wrong node for the element : " // << element // << "(coords[" << i << "] = " << coords_local[i] // << " and buffer[" << i << "] = " << coords(i) << ")"); // } // } #endif switch (tag){ case _gst_htm_capacity: { unpackNodalDataHelper(*capacity_lumped, buffer, elements); break; } case _gst_htm_temperature: { unpackNodalDataHelper(*temperature, buffer, elements); break; } case _gst_htm_gradient_temperature: { unpackElementalDataHelper(temperature_gradient, buffer, elements); unpackNodalDataHelper(*temperature, buffer, elements); // // Real tolerance = 1e-15; // if (!Math::are_vector_equal(spatial_dimension,gtemp.storage(),it_gtemp[element.element].storage())){ // Real dist = Math::distance_3d(gtemp.storage(), it_gtemp[element.element].storage()); // debug::debugger.getOutputStream().precision(20); // std::stringstream temperatures_str; // temperatures_str.precision(20); // temperatures_str << std::scientific << "temperatures are "; // for (UInt n = 0; n < nb_nodes_per_element; ++n) { // UInt offset_conn = conn[el_offset + n]; // temperatures_str << (*temperature)(offset_conn) << " "; // } // Vector::iterator it_shaped = // const_cast &>(getFEM().getShapesDerivatives(element.type, ghost_type)) // .begin(nb_nodes_per_element,spatial_dimension); // AKANTU_EXCEPTION("packed gradient do not match for element " << element.element << std::endl // << "buffer is " << gtemp << " local is " << it_gtemp[element.element] // << " dist is " << dist << std::endl // << temperatures_str.str() << std::endl // << std::scientific << std::setprecision(20) // << " distant temperatures " << temp_nodes // << "temperature gradient size " << temperature_gradient(element.type, ghost_type).getSize() // << " number of ghost elements " << getFEM().getMesh().getNbElement(element.type,_ghost) // << std::scientific << std::setprecision(20) // << " shaped " << shaped // << std::scientific << std::setprecision(20) // << " local shaped " << it_shaped[element.element]); // } break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } } /* -------------------------------------------------------------------------- */ diff --git a/src/model/model.hh b/src/model/model.hh index 66e93b96f..a0cdb1dc7 100644 --- a/src/model/model.hh +++ b/src/model/model.hh @@ -1,224 +1,226 @@ /** * @file model.hh * * @author Guillaume Anciaux * @author David Simon Kammer * @author Nicolas Richart * * @date Tue Jul 27 18:15:37 2010 * * @brief Interface of a 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MODEL_HH__ #define __AKANTU_MODEL_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "mesh.hh" #include "fem.hh" #include "mesh_utils.hh" #include "synchronizer_registry.hh" #include "distributed_synchronizer.hh" #include "static_communicator.hh" #include "mesh_partition.hh" #include "dof_synchronizer.hh" #include "pbc_synchronizer.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class Model : public Memory { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - + typedef Mesh mesh_type; - + Model(Mesh& mesh, const ID & id = "model", const MemoryID & memory_id = 0); - + virtual ~Model(); - + typedef std::map FEMMap; - + /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual void initModel() = 0; - + /// create the synchronizer registry object void createSynchronizerRegistry(DataAccessor * data_accessor); - + /// create a parallel synchronizer and distribute the mesh Synchronizer & createParallelSynch(MeshPartition * partition, DataAccessor * data_accessor); - + /// change local equation number so that PBC is assembled properly void changeLocalEquationNumberforPBC(std::map & pbc_pair,UInt dimension); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const = 0; - + /// initialize the model for PBC void setPBC(UInt x, UInt y, UInt z); void setPBC(SurfacePairList & surface_pairs, ElementType surface_e_type); - + virtual void initPBC(); - + /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get id of model AKANTU_GET_MACRO(ID, id, const ID) - + /// get the number of surfaces AKANTU_GET_MACRO(Mesh, mesh, Mesh&); - + /// return the object hadling equation numbers AKANTU_GET_MACRO(DOFSynchronizer, *dof_synchronizer, const DOFSynchronizer &) - + /// synchronize the boundary in case of parallel run virtual void synchronizeBoundaries() {}; - + /// return the fem object associated with a provided name inline FEM & getFEM(std::string name = "") const; - + /// return the fem boundary object associated with a provided name virtual FEM & getFEMBoundary(std::string name = ""); - + /// register a fem object associated with name template inline void registerFEMObject(const std::string & name, Mesh & mesh, UInt spatial_dimension); /// unregister a fem object associated with name inline void unRegisterFEMObject(const std::string & name); - + /// return the synchronizer registry SynchronizerRegistry & getSynchronizerRegistry(); - + public: /// return the fem object associated with a provided name template inline FEMClass & getFEMClass(std::string name = "") const; - + /// return the fem boundary object associated with a provided name template inline FEMClass & getFEMClassBoundary(std::string name = ""); - + /// get the pbc pairs std::map & getPBCPairs(){return pbc_pair;}; - + + /* ------------------------------------------------------------------------ */ + /* Pack and unpack helper functions */ + /* ------------------------------------------------------------------------ */ +public: + inline UInt getNbQuadraturePoints(const Vector & elements) const; + protected: /// returns if node is slave in pbc inline bool getIsPBCSlaveNode(const UInt node); - - /* ------------------------------------------------------------------------ */ - /* Class Members */ - /* ------------------------------------------------------------------------ */ + template inline void packElementalDataHelper(const ByElementTypeVector & data_to_pack, CommunicationBuffer & buffer, const Vector & elements) const; template inline void unpackElementalDataHelper(ByElementTypeVector & data_to_unpack, CommunicationBuffer & buffer, const Vector & elements) const; template inline void packUnpackElementalDataHelper(ByElementTypeVector & data_to_pack, CommunicationBuffer & buffer, const Vector & element) const; /* -------------------------------------------------------------------------- */ template inline void packNodalDataHelper(Vector & data_to_pack, CommunicationBuffer & buffer, const Vector & element) const; template inline void unpackNodalDataHelper(Vector & data_to_unpack, CommunicationBuffer & buffer, const Vector & element) const; template inline void packUnpackNodalDataHelper(Vector & data, CommunicationBuffer & buffer, const Vector & elements) const; - - + /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: - + /// Mesh Mesh & mesh; - + /// id ID id; - + /// the main fem object present in all models FEMMap fems; - + /// the fem object present in all models for boundaries FEMMap fems_boundary; - + /// default fem object std::string default_fem; - + /// synchronizer registry SynchronizerRegistry * synch_registry; - + /// handle the equation number things DOFSynchronizer * dof_synchronizer; - + /// pbc pairs std::map pbc_pair; - + /// flag per node to know is pbc slave Vector is_pbc_slave_node; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "model_inline_impl.cc" #endif /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Model & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_MODEL_HH__ */ diff --git a/src/model/model_inline_impl.cc b/src/model/model_inline_impl.cc index 31b389be0..c3513c4da 100644 --- a/src/model/model_inline_impl.cc +++ b/src/model/model_inline_impl.cc @@ -1,279 +1,290 @@ /** * @file model_inline_impl.cc * * @author Guillaume Anciaux * * @date Wed Aug 25 08:50:54 2010 * * @brief inline implementation of the model 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 . * */ /* -------------------------------------------------------------------------- */ inline SynchronizerRegistry & Model::getSynchronizerRegistry(){ AKANTU_DEBUG_ASSERT(synch_registry,"synchronizer registry not initialized:" << " did you call createSynchronizerRegistry ?"); return *synch_registry; } /* -------------------------------------------------------------------------- */ template inline FEMClass & Model::getFEMClassBoundary(std::string name) { AKANTU_DEBUG_IN(); if (name == "") name = default_fem; FEMMap::const_iterator it_boun = fems_boundary.find(name); FEMClass * tmp_fem_boundary; if (it_boun == fems_boundary.end()){ AKANTU_DEBUG_INFO("Creating FEM boundary " << name); FEMMap::const_iterator it = fems.find(name); AKANTU_DEBUG_ASSERT(it != fems.end(), "The FEM " << name << " is not registered"); UInt spatial_dimension = it->second->getElementDimension(); std::stringstream sstr; sstr << id << ":fem_boundary:" << name; tmp_fem_boundary = new FEMClass(it->second->getMesh(), spatial_dimension-1, sstr.str(), memory_id); fems_boundary[name] = tmp_fem_boundary; } else { tmp_fem_boundary = dynamic_cast(it_boun->second); } AKANTU_DEBUG_OUT(); return *tmp_fem_boundary; } /* -------------------------------------------------------------------------- */ template inline FEMClass & Model::getFEMClass(std::string name) const{ AKANTU_DEBUG_IN(); if (name == "") name = default_fem; FEMMap::const_iterator it = fems.find(name); AKANTU_DEBUG_ASSERT(it != fems.end(), "The FEM " << name << " is not registered"); AKANTU_DEBUG_OUT(); return dynamic_cast(*(it->second)); } /* -------------------------------------------------------------------------- */ inline void Model::unRegisterFEMObject(const std::string & name){ FEMMap::iterator it = fems.find(name); AKANTU_DEBUG_ASSERT(it != fems.end(), "FEM object with name " << name << " was not found"); delete((*it).second); fems.erase(it); if (!fems.empty()) default_fem = (*fems.begin()).first; } /* -------------------------------------------------------------------------- */ template inline void Model::registerFEMObject(const std::string & name, Mesh & mesh, UInt spatial_dimension){ if (fems.size() == 0) default_fem = name; #ifndef AKANTU_NDEBUG FEMMap::iterator it = fems.find(name); AKANTU_DEBUG_ASSERT(it == fems.end(), "FEM object with name " << name << " was already created"); #endif std::stringstream sstr; sstr << id << ":fem:" << name; fems[name] = new FEMClass(mesh, spatial_dimension, sstr.str(), memory_id); // MeshUtils::buildFacets(fems[name]->getMesh()); // std::stringstream sstr2; sstr2 << id << ":fem_boundary:" << name; // fems_boundary[name] = new FEMClass(mesh, spatial_dimension-1, sstr2.str(), memory_id); } /* -------------------------------------------------------------------------- */ inline FEM & Model::getFEM(std::string name) const{ AKANTU_DEBUG_IN(); if (name == "") name = default_fem; FEMMap::const_iterator it = fems.find(name); AKANTU_DEBUG_ASSERT(it != fems.end(),"The FEM " << name << " is not registered"); AKANTU_DEBUG_OUT(); return *(it->second); } /* -------------------------------------------------------------------------- */ inline FEM & Model::getFEMBoundary(std::string name){ AKANTU_DEBUG_IN(); if (name == "") name = default_fem; FEMMap::const_iterator it = fems_boundary.find(name); AKANTU_DEBUG_ASSERT(it != fems_boundary.end(), "The FEM boundary " << name << " is not registered"); AKANTU_DEBUG_ASSERT(it->second != NULL, "The FEM boundary " << name << " was not created"); AKANTU_DEBUG_OUT(); return *(it->second); } /* -------------------------------------------------------------------------- */ /// @todo : should merge with a single function which handles local and global inline void Model::changeLocalEquationNumberforPBC(std::map & pbc_pair, UInt dimension){ for (std::map::iterator it = pbc_pair.begin(); it != pbc_pair.end();++it) { Int node_master = (*it).second; Int node_slave = (*it).first; for (UInt i = 0; i < dimension; ++i) { (*dof_synchronizer->getLocalDOFEquationNumbersPointer()) (node_slave*dimension+i) = dimension*node_master+i; (*dof_synchronizer->getGlobalDOFEquationNumbersPointer()) (node_slave*dimension+i) = dimension*node_master+i; } } } /* -------------------------------------------------------------------------- */ inline bool Model::getIsPBCSlaveNode(const UInt node) { // if no pbc is defined, is_pbc_slave_node is of size zero if (is_pbc_slave_node.getSize() == 0) return false; else return is_pbc_slave_node(node); } +/* -------------------------------------------------------------------------- */ +inline UInt Model::getNbQuadraturePoints(const Vector & elements) const { + UInt nb_quad = 0; + Vector::const_iterator it = elements.begin(); + Vector::const_iterator end = elements.end(); + for (; it != end; ++it) { + const Element & el = *it; + nb_quad += getFEM().getNbQuadraturePoints(el.type, el.ghost_type); + } + return nb_quad; +} /* -------------------------------------------------------------------------- */ template inline void Model::packElementalDataHelper(const ByElementTypeVector & data_to_pack, CommunicationBuffer & buffer, const Vector & elements) const { packUnpackElementalDataHelper(const_cast &>(data_to_pack), buffer, elements); } /* -------------------------------------------------------------------------- */ template inline void Model::unpackElementalDataHelper(ByElementTypeVector & data_to_unpack, CommunicationBuffer & buffer, const Vector & elements) const { packUnpackElementalDataHelper(data_to_unpack, buffer, elements); } /* -------------------------------------------------------------------------- */ template inline void Model::packUnpackElementalDataHelper(ByElementTypeVector & data_to_pack, CommunicationBuffer & buffer, const Vector & element) const { ElementType current_element_type = _not_defined; GhostType current_ghost_type = _casper; UInt nb_quad_per_elem = 0; UInt nb_component = 0; Vector * vect = NULL; Vector::const_iterator it = element.begin(); Vector::const_iterator end = element.end(); for (; it != end; ++it) { const Element & el = *it; if(el.type != current_element_type || el.ghost_type != current_ghost_type) { current_element_type = el.type; current_ghost_type = el.ghost_type; vect = &data_to_pack(el.type, el.ghost_type); nb_quad_per_elem = this->getFEM().getNbQuadraturePoints(el.type, el.ghost_type); nb_component = vect->getNbComponent(); } types::Vector data(vect->storage() + el.element * nb_component * nb_quad_per_elem, nb_component * nb_quad_per_elem); if(pack_helper) buffer << data; else buffer >> data; } } /* -------------------------------------------------------------------------- */ template inline void Model::packNodalDataHelper(Vector & data_to_pack, CommunicationBuffer & buffer, const Vector & element) const { packUnpackNodalDataHelper(data_to_pack, buffer, element); } /* -------------------------------------------------------------------------- */ template inline void Model::unpackNodalDataHelper(Vector & data_to_unpack, CommunicationBuffer & buffer, const Vector & element) const { packUnpackNodalDataHelper(data_to_unpack, buffer, element); } /* -------------------------------------------------------------------------- */ template inline void Model::packUnpackNodalDataHelper(Vector & data, CommunicationBuffer & buffer, const Vector & elements) const { UInt nb_component = data.getNbComponent(); UInt nb_nodes_per_element = 0; ElementType current_element_type = _not_defined; GhostType current_ghost_type = _casper; UInt * conn = NULL; Vector::const_iterator it = elements.begin(); Vector::const_iterator end = elements.end(); for (; it != end; ++it) { const Element & el = *it; if(el.type != current_element_type || el.ghost_type != current_ghost_type) { current_element_type = el.type; current_ghost_type = el.ghost_type; conn = mesh.getConnectivity(el.type, el.ghost_type).storage(); nb_nodes_per_element = Mesh::getNbNodesPerElement(el.type); } UInt el_offset = el.element * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_conn = conn[el_offset + n]; types::Vector data_vect(data.storage() + offset_conn * nb_component, nb_component); if(pack_helper) buffer << data_vect; else buffer >> data_vect; } } } diff --git a/src/model/solid_mechanics/material.hh b/src/model/solid_mechanics/material.hh index 73faba994..d23fa302a 100644 --- a/src/model/solid_mechanics/material.hh +++ b/src/model/solid_mechanics/material.hh @@ -1,503 +1,500 @@ /** * @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 "data_accessor.hh" #include "material_parameters.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_HH__ #define __AKANTU_MATERIAL_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class Model; class SolidMechanicsModel; class CommunicationBuffer; } __BEGIN_AKANTU__ /** * 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, * Vector & tangent_matrix, * GhostType ghost_type = _not_ghost); * \endcode * */ class Material : protected Memory, public DataAccessor, public Parsable, public MeshEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Material(SolidMechanicsModel & model, const ID & id = ""); virtual ~Material(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// register a material parameter template void registerParam(std::string name, T & variable, T default_value, ParamAccessType type, std::string description = ""); template void registerParam(std::string name, T & variable, ParamAccessType type, std::string description = ""); template void registerInternal(__attribute__((unused)) ByElementTypeVector & vect) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// read parameter from file virtual bool parseParam(const std::string & key, const std::string & value, const ID & id); /// 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 void assembleResidual(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); /// 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(const ElementType type, Vector & result); /** * 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 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)) Vector & 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); /// transfer the B matrix to a Voigt notation B matrix template inline void transferBMatrixToSymVoigtBMatrix(const types::RMatrix & B, types::RMatrix & Bvoigt, UInt nb_nodes_per_element) const; inline UInt getTangentStiffnessVoigtSize(UInt spatial_dimension) const; /// compute the potential energy by element void computePotentialEnergyByElement(); /// compute the coordinates of the quadrature points void computeQuadraturePointsCoordinates(ByElementTypeReal & quadrature_points_coordinates, const GhostType & ghost_type) const; /// interpolate an elemental field on given points for each element template void interpolateElementalField(const Vector & field, Vector & result); /// template function to initialize the elemental field interpolation template void initElementalFieldInterpolation(const Vector & quad_coordinates, const Vector & interpolation_points_coordinates); /// build the coordinate matrix for the interpolation on elemental field template inline void buildElementalFieldInterpolationCoodinates(const types::RMatrix & coordinates, types::RMatrix & coordMatrix); /// get the size of the coordiante matrix used in the interpolation template inline UInt getSizeElementalFieldInterpolationCoodinates(); /* ------------------------------------------------------------------------ */ /* Function for all materials */ /* ------------------------------------------------------------------------ */ protected: /// compute the potential energy for a quadrature point inline void computePotentialEnergyOnQuad(types::RMatrix & grad_u, types::RMatrix & sigma, Real & epot); /// compute the potential energy for an element virtual void computePotentialEnergyByElement(ElementType type, UInt index, types::RVector & epot_on_quad_points); protected: /// allocate an internal vector template void initInternalVector(ByElementTypeVector & vect, UInt nb_component, bool temporary = false, ElementKind element_kind = _ek_regular); public: /// resize an internal vector template void resizeInternalVector(ByElementTypeVector & vect, ElementKind element_kind = _ek_regular) const; /* ------------------------------------------------------------------------ */ template inline void gradUToF(const types::RMatrix & grad_u, types::RMatrix & F); inline void rightCauchy(const types::RMatrix & F, types::RMatrix & C); inline void leftCauchy (const types::RMatrix & F, types::RMatrix & B); /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: virtual inline UInt getNbDataForElements(const Vector & elements, SynchronizationTag tag) const; virtual inline void packElementData(CommunicationBuffer & buffer, const Vector & elements, SynchronizationTag tag) const; virtual inline void unpackElementData(CommunicationBuffer & buffer, const Vector & elements, SynchronizationTag tag); template inline void packElementDataHelper(const ByElementTypeVector & data_to_pack, CommunicationBuffer & buffer, const Vector & elements) const; template inline void unpackElementDataHelper(ByElementTypeVector & data_to_unpack, CommunicationBuffer & buffer, const Vector & elements) const; protected: template inline void packUnpackElementDataHelper(ByElementTypeVector & data_to_pack, CommunicationBuffer & buffer, const Vector & elements) const; -public: - inline UInt getNbQuadraturePoints(const Vector & elements) const; - public: /* ------------------------------------------------------------------------ */ virtual inline void onElementsAdded(const Vector & element_list); virtual inline void onElementsRemoved(const Vector & element_list, const ByElementTypeUInt & new_numbering); protected: template void removeQuadraturePointsFromVectors(ByElementTypeVector & data, const ByElementTypeUInt & new_numbering); /* ------------------------------------------------------------------------ */ /* 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); bool isNonLocal() const { return is_non_local; } const Vector & getVector(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost) const; Vector & getVector(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost); template inline T getParam(const ID & param) const; template inline void setParam(const ID & param, T value); protected: bool isInit() const { return is_init; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id of the material ID id; /// material name std::string name; /// The model to witch the material belong SolidMechanicsModel * model; /// density : rho Real rho; /// stresses arrays ordered by element types ByElementTypeReal stress; /// strains arrays ordered by element types ByElementTypeReal strain; /// list of element handled by the material ByElementTypeUInt element_filter; /// is the vector for potential energy initialized // bool potential_energy_vector; /// potential energy by element ByElementTypeReal potential_energy; /// tell if using in non local mode or not bool is_non_local; /// spatial dimension UInt spatial_dimension; /// elemental field interpolation coordinates ByElementTypeReal interpolation_inverse_coordinates; /// elemental field interpolation points ByElementTypeReal interpolation_points_matrices; /// list of the paramters MaterialParameters params; private: /// boolean to know if the material has been initialized bool is_init; std::map internal_vectors_real; std::map internal_vectors_uint; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "material_inline_impl.cc" #endif /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Material & _this) { _this.printself(stream); return stream; } __END_AKANTU__ /* -------------------------------------------------------------------------- */ /* Auto loop */ /* -------------------------------------------------------------------------- */ #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type) \ Vector::iterator strain_it = \ this->strain(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ Vector::iterator strain_end = \ this->strain(el_type, ghost_type).end(spatial_dimension, \ spatial_dimension); \ Vector::iterator stress_it = \ this->stress(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ \ for(;strain_it != strain_end; ++strain_it, ++stress_it) { \ types::RMatrix & __attribute__((unused)) grad_u = *strain_it; \ types::RMatrix & __attribute__((unused)) sigma = *stress_it #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END \ } \ #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_mat) \ Vector::iterator strain_it = \ this->strain(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ Vector::iterator strain_end = \ this->strain(el_type, ghost_type).end(spatial_dimension, \ spatial_dimension); \ \ UInt tangent_size = \ this->getTangentStiffnessVoigtSize(spatial_dimension); \ Vector::iterator tangent_it = \ tangent_mat.begin(tangent_size, \ tangent_size); \ \ for(;strain_it != strain_end; ++strain_it, ++tangent_it) { \ types::RMatrix & __attribute__((unused)) grad_u = *strain_it; \ types::RMatrix & tangent = *tangent_it #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END \ } /* -------------------------------------------------------------------------- */ #define INSTANSIATE_MATERIAL(mat_name) \ template class mat_name<1>; \ template class mat_name<2>; \ template class mat_name<3> /* -------------------------------------------------------------------------- */ /* Material list */ /* -------------------------------------------------------------------------- */ // elastic materials #include "material_elastic.hh" #define AKANTU_CORE_MATERIAL_LIST \ ((2, (elastic , MaterialElastic ))) #if defined(AKANTU_EXTRA_MATERIALS) # include "material_extra_includes.hh" #else # define AKANTU_EXTRA_MATERIAL_LIST #endif #if defined(AKANTU_COHESIVE_ELEMENT) # include "material_cohesive_includes.hh" #else # define AKANTU_COHESIVE_MATERIAL_LIST #endif #if defined(AKANTU_DAMAGE_NON_LOCAL) # include "material_non_local_includes.hh" #else # define AKANTU_DAMAGE_NON_LOCAL_MATERIAL_LIST #endif #define AKANTU_MATERIAL_LIST \ AKANTU_CORE_MATERIAL_LIST \ AKANTU_EXTRA_MATERIAL_LIST \ AKANTU_COHESIVE_MATERIAL_LIST \ AKANTU_DAMAGE_NON_LOCAL_MATERIAL_LIST #endif /* __AKANTU_MATERIAL_HH__ */ diff --git a/src/model/solid_mechanics/material_inline_impl.cc b/src/model/solid_mechanics/material_inline_impl.cc index de6f26d25..fb0173c1c 100644 --- a/src/model/solid_mechanics/material_inline_impl.cc +++ b/src/model/solid_mechanics/material_inline_impl.cc @@ -1,462 +1,450 @@ /** * @file material_inline_impl.cc * * @author Marco Vocialta * @author Nicolas Richart * * @date Tue Jul 27 18:15:37 2010 * * @brief Implementation of the inline functions of the class material * * @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 . * */ __END_AKANTU__ #include "solid_mechanics_model.hh" #include __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ inline UInt Material::addElement(const ElementType & type, UInt element, const GhostType & ghost_type) { element_filter(type, ghost_type).push_back(element); return element_filter(type, ghost_type).getSize()-1; } /* -------------------------------------------------------------------------- */ inline UInt Material::getTangentStiffnessVoigtSize(UInt dim) const { return (dim * (dim - 1) / 2 + dim); } /* -------------------------------------------------------------------------- */ template inline void Material::gradUToF(const types::RMatrix & grad_u, types::RMatrix & F) { UInt size_F = F.size(); AKANTU_DEBUG_ASSERT(F.size() >= grad_u.size() && grad_u.size() == dim, "The dimension of the tensor F should be greater or equal to the dimension of the tensor grad_u."); for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) F(i, j) = grad_u(i, j); for (UInt i = 0; i < size_F; ++i) F(i, i) += 1; } /* -------------------------------------------------------------------------- */ inline void Material::rightCauchy(const types::RMatrix & F, types::RMatrix & C) { C.mul(F, F); } /* -------------------------------------------------------------------------- */ inline void Material::leftCauchy(const types::RMatrix & F, types::RMatrix & B) { B.mul(F, F); } /* -------------------------------------------------------------------------- */ inline void Material::computePotentialEnergyOnQuad(types::RMatrix & grad_u, types::RMatrix & sigma, Real & epot) { epot = 0.; for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) epot += sigma(i, j) * grad_u(i, j); epot *= .5; } /* -------------------------------------------------------------------------- */ template inline void Material::transferBMatrixToSymVoigtBMatrix(const types::RMatrix & B, types::RMatrix & Bvoigt, UInt nb_nodes_per_element) const { Bvoigt.clear(); for (UInt i = 0; i < dim; ++i) for (UInt n = 0; n < nb_nodes_per_element; ++n) Bvoigt(i, i + n*dim) = B(n, i); if(dim == 2) { ///in 2D, fill the @f$ [\frac{\partial N_i}{\partial x}, \frac{\partial N_i}{\partial y}]@f$ row for (UInt n = 0; n < nb_nodes_per_element; ++n) { Bvoigt(2, 1 + n*2) = B(n, 0); Bvoigt(2, 0 + n*2) = B(n, 1); } } if(dim == 3) { for (UInt n = 0; n < nb_nodes_per_element; ++n) { Real dndx = B(n, 0); Real dndy = B(n, 1); Real dndz = B(n, 2); ///in 3D, fill the @f$ [0, \frac{\partial N_i}{\partial y}, \frac{N_i}{\partial z}]@f$ row Bvoigt(3, 1 + n*3) = dndz; Bvoigt(3, 2 + n*3) = dndy; ///in 3D, fill the @f$ [\frac{\partial N_i}{\partial x}, 0, \frac{N_i}{\partial z}]@f$ row Bvoigt(4, 0 + n*3) = dndz; Bvoigt(4, 2 + n*3) = dndx; ///in 3D, fill the @f$ [\frac{\partial N_i}{\partial x}, \frac{N_i}{\partial y}, 0]@f$ row Bvoigt(5, 0 + n*3) = dndy; Bvoigt(5, 1 + n*3) = dndx; } } } /* -------------------------------------------------------------------------- */ template inline void Material::buildElementalFieldInterpolationCoodinates(__attribute__((unused)) const types::RMatrix & coordinates, __attribute__((unused)) types::RMatrix & coordMatrix) { AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_triangle_3>(const types::RMatrix & coordinates, types::RMatrix & coordMatrix) { for (UInt i = 0; i < coordinates.rows(); ++i) coordMatrix(i, 0) = 1; } /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_triangle_6>(const types::RMatrix & coordinates, types::RMatrix & coordMatrix) { UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(_triangle_6); for (UInt i = 0; i < coordinates.rows(); ++i) { coordMatrix(i, 0) = 1; for (UInt j = 1; j < nb_quadrature_points; ++j) coordMatrix(i, j) = coordinates(i, j-1); } } /** * @todo Write a more efficient interpolation for quadrangles by * dropping unnecessary quadrature points * */ /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_quadrangle_4>(const types::RMatrix & coordinates, types::RMatrix & coordMatrix) { for (UInt i = 0; i < coordinates.rows(); ++i) { Real x = coordinates(i, 0); Real y = coordinates(i, 1); coordMatrix(i, 0) = 1; coordMatrix(i, 1) = x; coordMatrix(i, 2) = y; coordMatrix(i, 3) = x * y; } } /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_quadrangle_8>(const types::RMatrix & coordinates, types::RMatrix & coordMatrix) { for (UInt i = 0; i < coordinates.rows(); ++i) { UInt j = 0; Real x = coordinates(i, 0); Real y = coordinates(i, 1); for (UInt e = 0; e <= 2; ++e) { for (UInt n = 0; n <= 2; ++n) { coordMatrix(i, j) = std::pow(x, e) * std::pow(y, n); ++j; } } } } /* -------------------------------------------------------------------------- */ template inline UInt Material::getSizeElementalFieldInterpolationCoodinates() { return model->getFEM().getNbQuadraturePoints(type); } /* -------------------------------------------------------------------------- */ template void Material::registerParam(std::string name, T & variable, T default_value, ParamAccessType type, std::string description) { AKANTU_DEBUG_IN(); params.registerParam(name, variable, default_value, type, description); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::registerParam(std::string name, T & variable, ParamAccessType type, std::string description) { AKANTU_DEBUG_IN(); params.registerParam(name, variable, type, description); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline UInt Material::getNbDataForElements(const Vector & elements, SynchronizationTag tag) const { if(tag == _gst_smm_stress) { - return spatial_dimension * spatial_dimension * sizeof(Real) * this->getNbQuadraturePoints(elements); + return spatial_dimension * spatial_dimension * sizeof(Real) * this->getModel().getNbQuadraturePoints(elements); } return 0; } /* -------------------------------------------------------------------------- */ inline void Material::packElementData(CommunicationBuffer & buffer, const Vector & elements, SynchronizationTag tag) const { if(tag == _gst_smm_stress) { packElementDataHelper(stress, buffer, elements); } } /* -------------------------------------------------------------------------- */ inline void Material::unpackElementData(CommunicationBuffer & buffer, const Vector & elements, SynchronizationTag tag) { if(tag == _gst_smm_stress) { unpackElementDataHelper(stress, buffer, elements); } } -/* -------------------------------------------------------------------------- */ -inline UInt Material::getNbQuadraturePoints(const Vector & elements) const { - UInt nb_quad = 0; - Vector::const_iterator it = elements.begin(); - Vector::const_iterator end = elements.end(); - for (; it != end; ++it) { - const Element & el = *it; - nb_quad += this->model->getFEM().getNbQuadraturePoints(el.type, el.ghost_type); - } - return nb_quad; -} - /* -------------------------------------------------------------------------- */ template inline void Material::packElementDataHelper(const ByElementTypeVector & data_to_pack, CommunicationBuffer & buffer, const Vector & elements) const { packUnpackElementDataHelper(const_cast &>(data_to_pack), buffer, elements); } /* -------------------------------------------------------------------------- */ template inline void Material::unpackElementDataHelper(ByElementTypeVector & data_to_unpack, CommunicationBuffer & buffer, const Vector & elements) const { packUnpackElementDataHelper(data_to_unpack, buffer, elements); } /* -------------------------------------------------------------------------- */ template inline void Material::packUnpackElementDataHelper(ByElementTypeVector & data_to_pack, CommunicationBuffer & buffer, const Vector & element) const { ElementType current_element_type = _not_defined; GhostType current_ghost_type = _casper; UInt nb_quad_per_elem = 0; UInt nb_component = 0; Vector * vect = NULL; Vector * element_index_material = NULL; Vector::const_iterator it = element.begin(); Vector::const_iterator end = element.end(); for (; it != end; ++it) { const Element & el = *it; if(el.type != current_element_type || el.ghost_type != current_ghost_type) { current_element_type = el.type; current_ghost_type = el.ghost_type; vect = &data_to_pack(el.type, el.ghost_type); element_index_material = &(this->model->getElementIndexByMaterial(current_element_type, current_ghost_type)); nb_quad_per_elem = this->model->getFEM().getNbQuadraturePoints(el.type, el.ghost_type); nb_component = vect->getNbComponent(); } UInt el_id = (*element_index_material)(el.element, 0); types::Vector data(vect->storage() + el_id * nb_component * nb_quad_per_elem, nb_component * nb_quad_per_elem); if(pack_helper) buffer << data; else buffer >> data; } } /* -------------------------------------------------------------------------- */ template inline T Material::getParam(const ID & param) const { try { return params.get(param); } catch (...) { AKANTU_EXCEPTION("No parameter " << param << " in the material " << id); } } /* -------------------------------------------------------------------------- */ template inline void Material::setParam(const ID & param, T value) { try { params.set(param, value); } catch(...) { AKANTU_EXCEPTION("No parameter " << param << " in the material " << id); } updateInternalParameters(); } /* -------------------------------------------------------------------------- */ template void Material::removeQuadraturePointsFromVectors(ByElementTypeVector & data, const ByElementTypeUInt & new_numbering) { for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; ByElementTypeVector::type_iterator it = new_numbering.firstType(0, gt, _ek_not_defined); ByElementTypeVector::type_iterator end = new_numbering.lastType(0, gt, _ek_not_defined); for (; it != end; ++it) { ElementType type = *it; if(data.exists(type, gt)){ const Vector & renumbering = new_numbering(type, gt); Vector & vect = data(type, gt); UInt nb_quad_per_elem = this->model->getFEM().getNbQuadraturePoints(type, gt); UInt nb_component = vect.getNbComponent(); Vector tmp(renumbering.getSize()*nb_quad_per_elem, nb_component); UInt new_size = 0; for (UInt i = 0; i < vect.getSize(); ++i) { UInt new_i = renumbering(i); if(new_i != UInt(-1)) { memcpy(tmp.storage() + new_i * nb_component * nb_quad_per_elem, vect.storage() + i * nb_component * nb_quad_per_elem, nb_component * nb_quad_per_elem * sizeof(T)); ++new_size; } } tmp.resize(new_size * nb_quad_per_elem); vect.copy(tmp); } } } } /* -------------------------------------------------------------------------- */ inline void Material::onElementsAdded(__attribute__((unused)) const Vector & element_list) { for (std::map::iterator it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) { resizeInternalVector(*(it->second)); } for (std::map::iterator it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) { resizeInternalVector(*(it->second)); } } /* -------------------------------------------------------------------------- */ inline void Material::onElementsRemoved(const Vector & element_list, const ByElementTypeUInt & new_numbering) { UInt my_num = model->getInternalIndexFromID(id); ByElementTypeUInt material_local_new_numbering("remove mat filter elem", id); Vector::const_iterator el_begin = element_list.begin(); Vector::const_iterator el_end = element_list.end(); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; ByElementTypeVector::type_iterator it = new_numbering.firstType(0, gt, _ek_not_defined); ByElementTypeVector::type_iterator end = new_numbering.lastType(0, gt, _ek_not_defined); for (; it != end; ++it) { ElementType type = *it; if(element_filter.exists(type, gt)){ Vector & elem_filter = element_filter(type, gt); Vector & element_index_material = model->getElementIndexByMaterial(type, gt); element_index_material.resize(model->getFEM().getMesh().getNbElement(type, gt)); // all materials will resize of the same size... if(!material_local_new_numbering.exists(type, gt)) material_local_new_numbering.alloc(elem_filter.getSize(), 1, type, gt); Vector & mat_renumbering = material_local_new_numbering(type, gt); const Vector & renumbering = new_numbering(type, gt); Vector elem_filter_tmp; UInt ni = 0; Element el; el.type = type; el.ghost_type = gt; for (UInt i = 0; i < elem_filter.getSize(); ++i) { el.element = elem_filter(i); if(std::find(el_begin, el_end, el) == el_end) { UInt new_el = renumbering(el.element); AKANTU_DEBUG_ASSERT(new_el != UInt(-1), "A not removed element as been badly renumbered"); elem_filter_tmp.push_back(new_el); mat_renumbering(i) = ni; element_index_material(new_el, 0) = ni; element_index_material(new_el, 1) = my_num; ++ni; } else { mat_renumbering(i) = UInt(-1); } } elem_filter.resize(elem_filter_tmp.getSize()); elem_filter.copy(elem_filter); } } } for (std::map::iterator it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) { this->removeQuadraturePointsFromVectors(*(it->second), material_local_new_numbering); } for (std::map::iterator it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) { this->removeQuadraturePointsFromVectors(*(it->second), material_local_new_numbering); } } diff --git a/src/model/solid_mechanics/materials/material_damage/material_marigo_inline_impl.cc b/src/model/solid_mechanics/materials/material_damage/material_marigo_inline_impl.cc index 6f6940d45..3c0bb17ea 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_marigo_inline_impl.cc +++ b/src/model/solid_mechanics/materials/material_damage/material_marigo_inline_impl.cc @@ -1,131 +1,131 @@ /** * @file material_marigo_inline_impl.cc * * @author Guillaume Anciaux * @author Marion Estelle Chambart * @author Nicolas Richart * * @date Thu Feb 02 11:09:36 2012 * * @brief Implementation of the inline functions of the material marigo * * @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 inline void MaterialMarigo::computeStressOnQuad(types::RMatrix & grad_u, types::RMatrix & sigma, Real & dam, Real & Y, Real &Ydq) { MaterialElastic::computeStressOnQuad(grad_u, sigma); Y = 0; for (UInt i = 0; i < spatial_dimension; ++i) { for (UInt j = 0; j < spatial_dimension; ++j) { Y += sigma(i,j) * grad_u(i,j); } } Y *= 0.5; if(damage_in_y) Y *= (1 - dam); if(yc_limit) Y = std::min(Y, Yc); if(!this->is_non_local) { computeDamageAndStressOnQuad(sigma, dam, Y, Ydq); } } /* -------------------------------------------------------------------------- */ template inline void MaterialMarigo::computeDamageAndStressOnQuad(types::RMatrix & sigma, Real & dam, Real & Y, Real &Ydq) { Real Fd = Y - Ydq - Sd * dam; if (Fd > 0) dam = (Y - Ydq) / Sd; dam = std::min(dam,1.); sigma *= 1-dam; } /* -------------------------------------------------------------------------- */ template inline UInt MaterialMarigo::getNbDataForElements(const Vector & elements, SynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt size = 0; if(tag == _gst_smm_init_mat) { - size += sizeof(Real) * this->getNbQuadraturePoints(elements); + size += sizeof(Real) * this->getModel().getNbQuadraturePoints(elements); } size += MaterialDamage::getNbDataForElements(elements, tag); AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ template inline void MaterialMarigo::packElementData(CommunicationBuffer & buffer, const Vector & elements, SynchronizationTag tag) const { AKANTU_DEBUG_IN(); if(tag == _gst_smm_init_mat) { this->packElementDataHelper(Yd_rand, buffer, elements); // UInt nb_quadrature_points = this->model->getFEM().getNbQuadraturePoints(element.type); // Vector::const_iterator Yds = Yd_rand(element.type, _not_ghost).begin(); // Yds += element.element * nb_quadrature_points; // for (UInt q = 0; q < nb_quadrature_points; ++q, ++Yds) // buffer << *Yds; } MaterialDamage::packElementData(buffer, elements, tag); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void MaterialMarigo::unpackElementData(CommunicationBuffer & buffer, const Vector & elements, SynchronizationTag tag) { AKANTU_DEBUG_IN(); if(tag == _gst_smm_init_mat) { this->unpackElementDataHelper(Yd_rand, buffer, elements); // UInt nb_quadrature_points = this->model->getFEM().getNbQuadraturePoints(element.type); // Vector::iterator Ydr = Yd_rand(element.type, _ghost).begin(); // Ydr += element.element * nb_quadrature_points; // for (UInt q = 0; q < nb_quadrature_points; ++q, ++Ydr) // buffer << *Ydr; } MaterialDamage::unpackElementData(buffer, elements, tag); AKANTU_DEBUG_OUT(); } diff --git a/src/model/solid_mechanics/materials/material_non_local_inline_impl.cc b/src/model/solid_mechanics/materials/material_non_local_inline_impl.cc index deab4ec78..380bd45ec 100644 --- a/src/model/solid_mechanics/materials/material_non_local_inline_impl.cc +++ b/src/model/solid_mechanics/materials/material_non_local_inline_impl.cc @@ -1,851 +1,851 @@ /** * @file material_non_local_inline_impl.cc * * @author Nicolas Richart * * @date Wed Aug 31 11:09:48 2011 * * @brief Non-local inline implementation * * @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 . * */ /* -------------------------------------------------------------------------- */ __END_AKANTU__ /* -------------------------------------------------------------------------- */ #include "aka_types.hh" #include "grid_synchronizer.hh" #include "synchronizer_registry.hh" #include "integrator.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template class WeightFunction> MaterialNonLocal::MaterialNonLocal(SolidMechanicsModel & model, const ID & id) : Material(model, id), radius(100.), weight_func(NULL), cell_list(NULL), update_weights(0), compute_stress_calls(0), is_creating_grid(false), grid_synchronizer(NULL) { AKANTU_DEBUG_IN(); this->registerParam("radius" , radius , 100., _pat_readable , "Non local radius"); this->registerParam("UpdateWeights" , update_weights, 0U , _pat_modifiable, "Update weights frequency"); this->registerParam("Weight function", weight_func , _pat_internal); this->is_non_local = true; this->weight_func = new WeightFunction(*this); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> MaterialNonLocal::~MaterialNonLocal() { AKANTU_DEBUG_IN(); delete cell_list; delete weight_func; delete grid_synchronizer; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::initMaterial() { AKANTU_DEBUG_IN(); // Material::initMaterial(); Mesh & mesh = this->model->getFEM().getMesh(); ByElementTypeReal quadrature_points_coordinates("quadrature_points_coordinates_tmp_nl", id); this->initInternalVector(quadrature_points_coordinates, spatial_dimension, true); computeQuadraturePointsCoordinates(quadrature_points_coordinates, _not_ghost); ByElementType nb_ghost_protected; Mesh::type_iterator it = mesh.firstType(spatial_dimension, _ghost); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _ghost); for(; it != last_type; ++it) nb_ghost_protected(mesh.getNbElement(*it, _ghost), *it, _ghost); createCellList(quadrature_points_coordinates); updatePairList(quadrature_points_coordinates); cleanupExtraGhostElement(nb_ghost_protected); weight_func->setRadius(radius); weight_func->init(); computeWeights(quadrature_points_coordinates); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::cleanupExtraGhostElement(const ByElementType & nb_ghost_protected) { AKANTU_DEBUG_IN(); // Create list of element to keep std::set relevant_ghost_element; pair_type::const_iterator first_pair_types = existing_pairs[1].begin(); pair_type::const_iterator last_pair_types = existing_pairs[1].end(); for (; first_pair_types != last_pair_types; ++first_pair_types) { ElementType type2 = first_pair_types->second; GhostType ghost_type2 = _ghost; UInt nb_quad2 = this->model->getFEM().getNbQuadraturePoints(type2); Vector & elem_filter = element_filter(type2, ghost_type2); const Vector & pairs = pair_list(first_pair_types->first, _not_ghost)(first_pair_types->second, ghost_type2); Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); for(;first_pair != last_pair; ++first_pair) { UInt _q2 = (*first_pair)(1); QuadraturePoint q2(type2, elem_filter(_q2 / nb_quad2), _q2 % nb_quad2, ghost_type2); relevant_ghost_element.insert(q2); } } // Create list of element to remove and new numbering for element to keep Mesh & mesh = this->model->getFEM().getMesh(); std::set ghost_to_erase; Mesh::type_iterator it = mesh.firstType(spatial_dimension, _ghost); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _ghost); RemovedElementsEvent remove_elem(mesh); Element element; element.ghost_type = _ghost; for(; it != last_type; ++it) { element.type = *it; UInt nb_ghost_elem = mesh.getNbElement(*it, _ghost); UInt nb_ghost_elem_protected = 0; try { nb_ghost_elem_protected = nb_ghost_protected(*it, _ghost); } catch (...) {} if(!remove_elem.getNewNumbering().exists(*it, _ghost)) remove_elem.getNewNumbering().alloc(nb_ghost_elem, 1, *it, _ghost); Vector & elem_filter = element_filter(*it, _ghost); Vector & new_numbering = remove_elem.getNewNumbering(*it, _ghost); UInt ng = 0; for (UInt g = 0; g < nb_ghost_elem; ++g) { element.element = elem_filter(g); if(element.element >= nb_ghost_elem_protected && (std::find(relevant_ghost_element.begin(), relevant_ghost_element.end(), element) == relevant_ghost_element.end())) { ghost_to_erase.insert(element); remove_elem.getList().push_back(element); new_numbering(g) = UInt(-1); } else { new_numbering(g) = ng; ++ng; } } } mesh.sendEvent(remove_elem); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::createCellList(ByElementTypeReal & quadrature_points_coordinates) { AKANTU_DEBUG_IN(); const Real safety_factor = 1.2; // for the cell grid spacing Mesh & mesh = this->model->getFEM().getMesh(); mesh.computeBoundingBox(); Real lower_bounds[spatial_dimension]; Real upper_bounds[spatial_dimension]; mesh.getLocalLowerBounds(lower_bounds); mesh.getLocalUpperBounds(upper_bounds); Real spacing[spatial_dimension]; for (UInt i = 0; i < spatial_dimension; ++i) { spacing[i] = radius * safety_factor; } cell_list = new RegularGrid(spatial_dimension, lower_bounds, upper_bounds, spacing); fillCellList(quadrature_points_coordinates, _not_ghost); is_creating_grid = true; SynchronizerRegistry & synch_registry = this->model->getSynchronizerRegistry(); std::stringstream sstr; sstr << id << ":grid_synchronizer"; grid_synchronizer = GridSynchronizer::createGridSynchronizer(mesh, *cell_list, sstr.str()); synch_registry.registerSynchronizer(*grid_synchronizer, _gst_mnl_for_average); synch_registry.registerSynchronizer(*grid_synchronizer, _gst_mnl_weight); is_creating_grid = false; this->computeQuadraturePointsCoordinates(quadrature_points_coordinates, _ghost); fillCellList(quadrature_points_coordinates, _ghost); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::fillCellList(const ByElementTypeReal & quadrature_points_coordinates, const GhostType & ghost_type) { Mesh & mesh = this->model->getFEM().getMesh(); QuadraturePoint q; q.ghost_type = ghost_type; Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType (spatial_dimension, ghost_type); for(; it != last_type; ++it) { Vector & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_quad = this->model->getFEM().getNbQuadraturePoints(*it, ghost_type); const Vector & quads = quadrature_points_coordinates(*it, ghost_type); q.type = *it; Vector::const_iterator quad = quads.begin(spatial_dimension); UInt * elem = elem_filter.storage(); for (UInt e = 0; e < nb_element; ++e) { q.element = *elem; for (UInt nq = 0; nq < nb_quad; ++nq) { q.num_point = nq; q.setPosition(*quad); cell_list->insert(q, *quad); ++quad; } ++elem; } } } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::updatePairList(const ByElementTypeReal & quadrature_points_coordinates) { AKANTU_DEBUG_IN(); Mesh & mesh = this->model->getFEM().getMesh(); GhostType ghost_type = _not_ghost; // generate the pair of neighbor depending of the cell_list Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { // Preparing datas const Vector & quads = quadrature_points_coordinates(*it, ghost_type); Vector::const_iterator first_quad = quads.begin(spatial_dimension); Vector::const_iterator last_quad = quads.end(spatial_dimension); ByElementTypeUInt & pairs = pair_list(ByElementTypeUInt("pairs", id, memory_id), *it, ghost_type); ElementType current_element_type = _not_defined; GhostType current_ghost_type = _casper; UInt existing_pairs_num = 0; Vector * neighbors = NULL; Vector * element_index_material2 = NULL; UInt my_num_quad = 0; // loop over quad points for(;first_quad != last_quad; ++first_quad, ++my_num_quad) { RegularGrid::Cell cell = cell_list->getCell(*first_quad); RegularGrid::neighbor_cells_iterator first_neigh_cell = cell_list->beginNeighborCells(cell); RegularGrid::neighbor_cells_iterator last_neigh_cell = cell_list->endNeighborCells(cell); // loop over neighbors cells of the one containing the current quadrature // point for (; first_neigh_cell != last_neigh_cell; ++first_neigh_cell) { RegularGrid::iterator first_neigh_quad = cell_list->beginCell(*first_neigh_cell); RegularGrid::iterator last_neigh_quad = cell_list->endCell(*first_neigh_cell); // loop over the quadrature point in the current cell of the cell list for (;first_neigh_quad != last_neigh_quad; ++first_neigh_quad){ QuadraturePoint & quad = *first_neigh_quad; UInt nb_quad_per_elem = this->model->getFEM().getNbQuadraturePoints(quad.type, quad.ghost_type); // little optimization to not search in the map at each quad points if(quad.type != current_element_type || quad.ghost_type != current_ghost_type) { current_element_type = quad.type; current_ghost_type = quad.ghost_type; existing_pairs_num = quad.ghost_type == _not_ghost ? 0 : 1; if(!pairs.exists(current_element_type, current_ghost_type)) { neighbors = &(pairs.alloc(0, 2, current_element_type, current_ghost_type)); } else { neighbors = &(pairs(current_element_type, current_ghost_type)); } existing_pairs[existing_pairs_num].insert(std::pair(*it, current_element_type)); element_index_material2 = &(this->model->getElementIndexByMaterial(current_element_type, current_ghost_type)); } UInt neigh_num_quad = (*element_index_material2)(quad.element) * nb_quad_per_elem + quad.num_point; const types::RVector & neigh_quad = quad.getPosition(); Real distance = first_quad->distance(neigh_quad); if(distance <= radius && (current_ghost_type == _ghost || (current_ghost_type == _not_ghost && my_num_quad <= neigh_num_quad))) { // sotring only half lists UInt pair[2]; pair[0] = my_num_quad; pair[1] = neigh_num_quad; neighbors->push_back(pair); } } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::computeWeights(const ByElementTypeReal & quadrature_points_coordinates) { AKANTU_DEBUG_IN(); GhostType ghost_type1; ghost_type1 = _not_ghost; ByElementTypeReal quadrature_points_volumes("quadrature_points_volumes", id, memory_id); this->initInternalVector(quadrature_points_volumes, 1, true); resizeInternalVector(quadrature_points_volumes); const ByElementTypeReal & jacobians_by_type = this->model->getFEM().getIntegratorInterface().getJacobians(); weight_func->updateInternals(quadrature_points_volumes); for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; UInt existing_pairs_num = gt - _not_ghost; pair_type::iterator first_pair_types = existing_pairs[existing_pairs_num].begin(); pair_type::iterator last_pair_types = existing_pairs[existing_pairs_num].end(); // Compute the weights for (; first_pair_types != last_pair_types; ++first_pair_types) { ElementType type1 = first_pair_types->first; ElementType type2 = first_pair_types->second; const Vector & pairs = pair_list(type1, ghost_type1)(type2, ghost_type2); std::string ghost_id = ""; if (ghost_type1 == _ghost) ghost_id = ":ghost"; ByElementTypeReal & weights_type_1 = pair_weight(type1, ghost_type1); std::stringstream sstr; sstr << id << ":pair_weight:" << type1 << ghost_id; weights_type_1.setID(sstr.str()); Vector * tmp_weight = NULL; if(!weights_type_1.exists(type2, ghost_type2)) { tmp_weight = &(weights_type_1.alloc(0, 2, type2, ghost_type2)); } else { tmp_weight = &(weights_type_1(type2, ghost_type2)); } Vector & weights = *tmp_weight; weights.resize(pairs.getSize()); weights.clear(); const Vector & jacobians_1 = jacobians_by_type(type1, ghost_type1); const Vector & jacobians_2 = jacobians_by_type(type2, ghost_type2); UInt nb_quad1 = this->model->getFEM().getNbQuadraturePoints(type1); UInt nb_quad2 = this->model->getFEM().getNbQuadraturePoints(type2); Vector & quads_volumes1 = quadrature_points_volumes(type1, ghost_type1); Vector & quads_volumes2 = quadrature_points_volumes(type2, ghost_type2); Vector::const_iterator iquads1; Vector::const_iterator iquads2; iquads1 = quadrature_points_coordinates(type1, ghost_type1).begin(spatial_dimension); iquads2 = quadrature_points_coordinates(type2, ghost_type2).begin(spatial_dimension); Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector::iterator weight = weights.begin(2); this->weight_func->selectType(type1, ghost_type1, type2, ghost_type2); // Weight function for(;first_pair != last_pair; ++first_pair, ++weight) { UInt _q1 = (*first_pair)(0); UInt _q2 = (*first_pair)(1); const types::RVector & pos1 = iquads1[_q1]; const types::RVector & pos2 = iquads2[_q2]; QuadraturePoint q1(_q1 / nb_quad1, _q1 % nb_quad1, _q1, pos1, type1, ghost_type1); QuadraturePoint q2(_q2 / nb_quad2, _q2 % nb_quad2, _q2, pos2, type2, ghost_type2); Real r = pos1.distance(pos2); Real w2J2 = jacobians_2(_q2); (*weight)(0) = w2J2 * this->weight_func->operator()(r, q1, q2); if(ghost_type2 != _ghost && _q1 != _q2) { Real w1J1 = jacobians_1(_q1); (*weight)(1) = w1J1 * this->weight_func->operator()(r, q2, q1); } else (*weight)(1) = 0; quads_volumes1(_q1) += (*weight)(0); if(ghost_type2 != _ghost) quads_volumes2(_q2) += (*weight)(1); } } } //normalize the weights for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; UInt existing_pairs_num = gt - _not_ghost; pair_type::iterator first_pair_types = existing_pairs[existing_pairs_num].begin(); pair_type::iterator last_pair_types = existing_pairs[existing_pairs_num].end(); first_pair_types = existing_pairs[existing_pairs_num].begin(); for (; first_pair_types != last_pair_types; ++first_pair_types) { ElementType type1 = first_pair_types->first; ElementType type2 = first_pair_types->second; const Vector & pairs = pair_list(type1, ghost_type1)(type2, ghost_type2); Vector & weights = pair_weight(type1, ghost_type1)(type2, ghost_type2); Vector & quads_volumes1 = quadrature_points_volumes(type1, ghost_type1); Vector & quads_volumes2 = quadrature_points_volumes(type2, ghost_type2); Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector::iterator weight = weights.begin(2); for(;first_pair != last_pair; ++first_pair, ++weight) { UInt q1 = (*first_pair)(0); UInt q2 = (*first_pair)(1); (*weight)(0) *= 1. / quads_volumes1(q1); if(ghost_type2 != _ghost) (*weight)(1) *= 1. / quads_volumes2(q2); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> template void MaterialNonLocal::weightedAvergageOnNeighbours(const ByElementTypeVector & to_accumulate, ByElementTypeVector & accumulated, UInt nb_degree_of_freedom, GhostType ghost_type2) const { AKANTU_DEBUG_IN(); UInt existing_pairs_num = 0; if (ghost_type2 == _ghost) existing_pairs_num = 1; pair_type::const_iterator first_pair_types = existing_pairs[existing_pairs_num].begin(); pair_type::const_iterator last_pair_types = existing_pairs[existing_pairs_num].end(); GhostType ghost_type1 = _not_ghost; // does not make sens the ghost vs ghost so this should always by _not_ghost for (; first_pair_types != last_pair_types; ++first_pair_types) { const Vector & pairs = pair_list(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); const Vector & weights = pair_weight(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); const Vector & to_acc = to_accumulate(first_pair_types->second, ghost_type2); Vector & acc = accumulated(first_pair_types->first, ghost_type1); if(ghost_type2 == _not_ghost) acc.clear(); Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector::const_iterator< types::Vector > pair_w = weights.begin(2); for(;first_pair != last_pair; ++first_pair, ++pair_w) { UInt q1 = (*first_pair)(0); UInt q2 = (*first_pair)(1); for(UInt d = 0; d < nb_degree_of_freedom; ++d){ acc(q1, d) += (*pair_w)(0) * to_acc(q2, d); if(ghost_type2 != _ghost) acc(q2, d) += (*pair_w)(1) * to_acc(q1, d); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::updateResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); // Update the weights for the non local variable averaging if(ghost_type == _not_ghost && this->update_weights && (this->compute_stress_calls % this->update_weights == 0)) { ByElementTypeReal quadrature_points_coordinates("quadrature_points_coordinates", id); Mesh & mesh = this->model->getFEM().getMesh(); mesh.initByElementTypeVector(quadrature_points_coordinates, spatial_dimension, 0); computeQuadraturePointsCoordinates(quadrature_points_coordinates, _not_ghost); computeQuadraturePointsCoordinates(quadrature_points_coordinates, _ghost); computeWeights(quadrature_points_coordinates); } if(ghost_type == _not_ghost) ++this->compute_stress_calls; computeAllStresses(ghost_type); computeNonLocalStresses(ghost_type); assembleResidual(ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::computeAllNonLocalStresses(GhostType ghost_type) { // Update the weights for the non local variable averaging if(ghost_type == _not_ghost) { if(this->update_weights && (this->compute_stress_calls % this->update_weights == 0)) { this->model->getSynchronizerRegistry().asynchronousSynchronize(_gst_mnl_weight); ByElementTypeReal quadrature_points_coordinates("quadrature_points_coordinates", id); Mesh & mesh = this->model->getFEM().getMesh(); mesh.initByElementTypeVector(quadrature_points_coordinates, spatial_dimension, 0); computeQuadraturePointsCoordinates(quadrature_points_coordinates, _not_ghost); computeQuadraturePointsCoordinates(quadrature_points_coordinates, _ghost); this->model->getSynchronizerRegistry().waitEndSynchronize(_gst_mnl_weight); computeWeights(quadrature_points_coordinates); } typename std::map::iterator it = non_local_variables.begin(); typename std::map::iterator end = non_local_variables.end(); for(;it != end; ++it) { NonLocalVariable & non_local_variable = it->second; resizeInternalVector(*non_local_variable.non_local_variable); this->weightedAvergageOnNeighbours(*non_local_variable.local_variable, *non_local_variable.non_local_variable, non_local_variable.non_local_variable_nb_component, _not_ghost); } ++this->compute_stress_calls; } else { typename std::map::iterator it = non_local_variables.begin(); typename std::map::iterator end = non_local_variables.end(); for(;it != end; ++it) { NonLocalVariable & non_local_variable = it->second; this->weightedAvergageOnNeighbours(*non_local_variable.local_variable, *non_local_variable.non_local_variable, non_local_variable.non_local_variable_nb_component, _ghost); } computeNonLocalStresses(_not_ghost); } } /* -------------------------------------------------------------------------- */ template class WeightFunction> bool MaterialNonLocal::parseParam(const std::string & key, const std::string & value, __attribute__((unused)) const ID & id) { std::stringstream sstr(value); if(key == "radius") { sstr >> radius; } else if(key == "UpdateWeights") { sstr >> update_weights; } else if(!weight_func->parseParam(key, value)) return Material::parseParam(key, value, id); return true; } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::savePairs(const std::string & filename) const { std::ofstream pout; std::stringstream sstr; StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Int prank = comm.whoAmI(); sstr << filename << "." << prank; pout.open(sstr.str().c_str()); GhostType ghost_type1; ghost_type1 = _not_ghost; for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; UInt existing_pairs_num = gt - _not_ghost; pair_type::const_iterator first_pair_types = existing_pairs[existing_pairs_num].begin(); pair_type::const_iterator last_pair_types = existing_pairs[existing_pairs_num].end(); for (; first_pair_types != last_pair_types; ++first_pair_types) { const Vector & pairs = pair_list(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); const Vector & weights = pair_weight(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); pout << "Types : " << first_pair_types->first << " (" << ghost_type1 << ") - " << first_pair_types->second << " (" << ghost_type2 << ")" << std::endl; Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector::const_iterator pair_w = weights.begin(2); for(;first_pair != last_pair; ++first_pair, ++pair_w) { UInt q1 = (*first_pair)(0); UInt q2 = (*first_pair)(1); pout << q1 << " " << q2 << " "<< (*pair_w)(0) << " " << (*pair_w)(1) << std::endl; } } } } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::neighbourhoodStatistics(const std::string & filename) const { std::ofstream pout; pout.open(filename.c_str()); const Mesh & mesh = this->model->getFEM().getMesh(); GhostType ghost_type1; ghost_type1 = _not_ghost; StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Int prank = comm.whoAmI(); for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; UInt existing_pairs_num = gt - _not_ghost; ByElementTypeUInt nb_neighbors("nb_neighbours", id, memory_id); mesh.initByElementTypeVector(nb_neighbors, 1, spatial_dimension); resizeInternalVector(nb_neighbors); pair_type::const_iterator first_pair_types = existing_pairs[existing_pairs_num].begin(); pair_type::const_iterator last_pair_types = existing_pairs[existing_pairs_num].end(); for (; first_pair_types != last_pair_types; ++first_pair_types) { const Vector & pairs = pair_list(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); if(prank == 0) { pout << ghost_type2 << " "; pout << "Types : " << first_pair_types->first << " " << first_pair_types->second << std::endl; } Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector & nb_neigh_1 = nb_neighbors(first_pair_types->first, ghost_type1); Vector & nb_neigh_2 = nb_neighbors(first_pair_types->second, ghost_type2); for(;first_pair != last_pair; ++first_pair) { UInt q1 = (*first_pair)(0); UInt q2 = (*first_pair)(1); ++(nb_neigh_1(q1)); if(q1 != q2) ++(nb_neigh_2(q2)); } } Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type1); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type1); UInt nb_quads = 0; Real sum_nb_neig = 0; UInt max_nb_neig = 0; UInt min_nb_neig = std::numeric_limits::max(); for(; it != last_type; ++it) { Vector & nb_neighor = nb_neighbors(*it, ghost_type1); Vector::iterator nb_neigh = nb_neighor.begin(); Vector::iterator end_neigh = nb_neighor.end(); for (; nb_neigh != end_neigh; ++nb_neigh, ++nb_quads) { UInt nb = *nb_neigh; sum_nb_neig += nb; max_nb_neig = std::max(max_nb_neig, nb); min_nb_neig = std::min(min_nb_neig, nb); } } comm.allReduce(&nb_quads, 1, _so_sum); comm.allReduce(&sum_nb_neig, 1, _so_sum); comm.allReduce(&max_nb_neig, 1, _so_max); comm.allReduce(&min_nb_neig, 1, _so_min); if(prank == 0) { pout << ghost_type2 << " "; pout << "Nb quadrature points: " << nb_quads << std::endl; Real mean_nb_neig = sum_nb_neig / Real(nb_quads); pout << ghost_type2 << " "; pout << "Average nb neighbors: " << mean_nb_neig << "(" << sum_nb_neig << ")" << std::endl; pout << ghost_type2 << " "; pout << "Max nb neighbors: " << max_nb_neig << std::endl; pout << ghost_type2 << " "; pout << "Min nb neighbors: " << min_nb_neig << std::endl; } } pout.close(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> inline UInt MaterialNonLocal::getNbDataForElements(const Vector & elements, SynchronizationTag tag) const { - UInt nb_quadrature_points = this->getNbQuadraturePoints(elements); + UInt nb_quadrature_points = this->getModel().getNbQuadraturePoints(elements); UInt size = 0; if(tag == _gst_mnl_for_average) { typename std::map::const_iterator it = non_local_variables.begin(); typename std::map::const_iterator end = non_local_variables.end(); for(;it != end; ++it) { const NonLocalVariable & non_local_variable = it->second; size += non_local_variable.non_local_variable_nb_component * sizeof(Real) * nb_quadrature_points; } } size += weight_func->getNbDataForElements(elements, tag); return size; } /* -------------------------------------------------------------------------- */ template class WeightFunction> inline void MaterialNonLocal::packElementData(CommunicationBuffer & buffer, const Vector & elements, SynchronizationTag tag) const { if(tag == _gst_mnl_for_average) { typename std::map::const_iterator it = non_local_variables.begin(); typename std::map::const_iterator end = non_local_variables.end(); for(;it != end; ++it) { const NonLocalVariable & non_local_variable = it->second; this->packElementDataHelper(*non_local_variable.local_variable, buffer, elements); } } weight_func->packElementData(buffer, elements, tag); } /* -------------------------------------------------------------------------- */ template class WeightFunction> inline void MaterialNonLocal::unpackElementData(CommunicationBuffer & buffer, const Vector & elements, SynchronizationTag tag) { if(tag == _gst_mnl_for_average) { typename std::map::iterator it = non_local_variables.begin(); typename std::map::iterator end = non_local_variables.end(); for(;it != end; ++it) { NonLocalVariable & non_local_variable = it->second; this->unpackElementDataHelper(*non_local_variable.local_variable, buffer, elements); } } weight_func->unpackElementData(buffer, elements, tag); } /* -------------------------------------------------------------------------- */ template class WeightFunction> inline void MaterialNonLocal::onElementsAdded(const Vector & element_list) { AKANTU_DEBUG_IN(); Material::onElementsAdded(element_list); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> inline void MaterialNonLocal::onElementsRemoved(const Vector & element_list, const ByElementTypeUInt & new_numbering) { AKANTU_DEBUG_IN(); Material::onElementsRemoved(element_list, new_numbering); pair_type::const_iterator first_pair_types = existing_pairs[1].begin(); pair_type::const_iterator last_pair_types = existing_pairs[1].end(); // Renumber element to keep for (; first_pair_types != last_pair_types; ++first_pair_types) { ElementType type2 = first_pair_types->second; GhostType ghost_type2 = _ghost; UInt nb_quad2 = this->model->getFEM().getNbQuadraturePoints(type2); Vector & pairs = pair_list(first_pair_types->first, _not_ghost)(first_pair_types->second, ghost_type2); Vector::iterator< types::Vector > first_pair = pairs.begin(2); Vector::iterator< types::Vector > last_pair = pairs.end(2); for(;first_pair != last_pair; ++first_pair) { UInt _q2 = (*first_pair)(1); const Vector & renumbering = new_numbering(type2, ghost_type2); UInt el = _q2 / nb_quad2; UInt new_el = renumbering(el); AKANTU_DEBUG_ASSERT(new_el != UInt(-1), "A local quad as been removed instead f just renumbered"); (*first_pair)(1) = new_el * nb_quad2 + _q2 % nb_quad2; } } AKANTU_DEBUG_OUT(); } diff --git a/src/model/solid_mechanics/materials/weight_function.hh b/src/model/solid_mechanics/materials/weight_function.hh index db1d65f3c..a2beec7d0 100644 --- a/src/model/solid_mechanics/materials/weight_function.hh +++ b/src/model/solid_mechanics/materials/weight_function.hh @@ -1,391 +1,391 @@ /** * @file weight_function.hh * * @author Nicolas Richart * * @date Fri Apr 13 16:49:32 2012 * * @brief Weight functions for non local 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_types.hh" #include "solid_mechanics_model.hh" #include /* -------------------------------------------------------------------------- */ #include #ifndef __AKANTU_WEIGHT_FUNCTION_HH__ #define __AKANTU_WEIGHT_FUNCTION_HH__ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* Normal weight function */ /* -------------------------------------------------------------------------- */ template class BaseWeightFunction { public: BaseWeightFunction(Material & material) : material(material) {} virtual ~BaseWeightFunction() {} virtual void init() {}; virtual void updateInternals(__attribute__((unused)) const ByElementTypeReal & quadrature_points_coordinates) {}; /* ------------------------------------------------------------------------ */ inline void setRadius(Real radius) { R = radius; R2 = R * R; } /* ------------------------------------------------------------------------ */ inline void selectType(__attribute__((unused)) ElementType type1, __attribute__((unused)) GhostType ghost_type1, __attribute__((unused)) ElementType type2, __attribute__((unused)) GhostType ghost_type2) { } /* ------------------------------------------------------------------------ */ inline Real operator()(Real r, __attribute__((unused)) QuadraturePoint & q1, __attribute__((unused)) QuadraturePoint & q2) { Real w = 0; if(r <= R) { Real alpha = (1. - r*r / R2); w = alpha * alpha; // *weight = 1 - sqrt(r / radius); } return w; } /* ------------------------------------------------------------------------ */ bool parseParam(__attribute__((unused)) const std::string & key, __attribute__((unused)) const std::string & value) { return false; } /* ------------------------------------------------------------------------ */ virtual void printself(std::ostream & stream) const { stream << "BaseWeightFunction"; } public: virtual UInt getNbDataForElements(__attribute__((unused)) const Vector & elements, __attribute__((unused)) SynchronizationTag tag) const { return 0; } virtual inline void packElementData(__attribute__((unused)) CommunicationBuffer & buffer, __attribute__((unused)) const Vector & elements, __attribute__((unused)) SynchronizationTag tag) const {} virtual inline void unpackElementData(__attribute__((unused)) CommunicationBuffer & buffer, __attribute__((unused)) const Vector & elements, __attribute__((unused)) SynchronizationTag tag) {} protected: Material & material; Real R; Real R2; }; /* -------------------------------------------------------------------------- */ /* Damage weight function */ /* -------------------------------------------------------------------------- */ template class DamagedWeightFunction : public BaseWeightFunction { public: DamagedWeightFunction(Material & material) : BaseWeightFunction(material) {} inline void selectType(__attribute__((unused)) ElementType type1, __attribute__((unused)) GhostType ghost_type1, ElementType type2, GhostType ghost_type2) { selected_damage = &this->material.getVector("damage", type2, ghost_type2); } inline Real operator()(Real r, __attribute__((unused)) QuadraturePoint & q1, QuadraturePoint & q2) { UInt quad = q2.global_num; Real D = (*selected_damage)(quad); Real Radius = (1.-D)*(1.-D) * this->R2; if(Radius < Math::getTolerance()) { Radius = 0.01 * 0.01 * this->R2; } Real alpha = std::max(0., 1. - r*r / Radius); Real w = alpha * alpha; return w; } /* ------------------------------------------------------------------------ */ bool parseParam(__attribute__((unused)) const std::string & key, __attribute__((unused)) const std::string & value) { return false; } /* ------------------------------------------------------------------------ */ virtual void printself(std::ostream & stream) const { stream << "DamagedWeightFunction"; } private: const Vector * selected_damage; }; /* -------------------------------------------------------------------------- */ /* Remove damaged weight function */ /* -------------------------------------------------------------------------- */ template class RemoveDamagedWeightFunction : public BaseWeightFunction { public: RemoveDamagedWeightFunction(Material & material) : BaseWeightFunction(material) {} inline void selectType(__attribute__((unused)) ElementType type1, __attribute__((unused)) GhostType ghost_type1, ElementType type2, GhostType ghost_type2) { selected_damage = &this->material.getVector("damage", type2, ghost_type2); } inline Real operator()(Real r, __attribute__((unused)) QuadraturePoint & q1, QuadraturePoint & q2) { UInt quad = q2.global_num; if(q1.global_num == quad) return 1.; Real D = (*selected_damage)(quad); Real w = 0.; if(D < damage_limit) { Real alpha = std::max(0., 1. - r*r / this->R2); w = alpha * alpha; } return w; } /* ------------------------------------------------------------------------ */ bool parseParam(const std::string & key, const std::string & value) { std::stringstream sstr(value); if(key == "damage_limit") { sstr >> damage_limit; } else return BaseWeightFunction::parseParam(key, value); return true; } /* ------------------------------------------------------------------------ */ virtual void printself(std::ostream & stream) const { stream << "RemoveDamagedWeightFunction [damage_limit: " << damage_limit << "]"; } virtual UInt getNbDataForElements(const Vector & elements, SynchronizationTag tag) const { if(tag == _gst_mnl_weight) - return this->material.getNbQuadraturePoints(elements) * sizeof(Real); + return this->material.getModel().getNbQuadraturePoints(elements) * sizeof(Real); return 0; } virtual inline void packElementData(CommunicationBuffer & buffer, const Vector & elements, SynchronizationTag tag) const { if(tag == _gst_mnl_weight) this->material.packElementDataHelper(dynamic_cast &>(this->material).getDamage(), buffer, elements); } virtual inline void unpackElementData(CommunicationBuffer & buffer, const Vector & elements, SynchronizationTag tag) { if(tag == _gst_mnl_weight) this->material.unpackElementDataHelper(dynamic_cast< MaterialDamage &>(this->material).getDamage(), buffer, elements); } private: /// limit at which a point is considered as complitely broken Real damage_limit; /// internal pointer to the current damage vector const Vector * selected_damage; }; /* -------------------------------------------------------------------------- */ /* Remove damaged with damage rate weight function */ /* -------------------------------------------------------------------------- */ template class RemoveDamagedWithDamageRateWeightFunction : public BaseWeightFunction { public: RemoveDamagedWithDamageRateWeightFunction(Material & material) : BaseWeightFunction(material) {} inline void selectType(__attribute__((unused)) ElementType type1, __attribute__((unused)) GhostType ghost_type1, ElementType type2, GhostType ghost_type2) { selected_damage_with_damage_rate = &(this->material.getVector("damage",type2, ghost_type2)); selected_damage_rate_with_damage_rate = &(this->material.getVector("damage-rate",type2, ghost_type2)); } inline Real operator()(Real r, __attribute__((unused)) QuadraturePoint & q1, QuadraturePoint & q2) { UInt quad = q2.global_num; if(q1.global_num == quad) return 1.; Real D = (*selected_damage_with_damage_rate)(quad); Real w = 0.; Real alphaexp = 1.; Real betaexp = 2.; if(D < damage_limit_with_damage_rate) { Real alpha = std::max(0., 1. - pow((r*r / this->R2),alphaexp)); w = pow(alpha, betaexp); } return w; } /* ------------------------------------------------------------------------ */ bool setParam(const std::string & key, const std::string & value) { std::stringstream sstr(value); if(key == "damage_limit") { sstr >> damage_limit_with_damage_rate; } else return BaseWeightFunction::setParam(key, value); return true; } /* ------------------------------------------------------------------------ */ virtual void printself(std::ostream & stream) const { stream << "RemoveDamagedWithDamageRateWeightFunction [damage_limit: " << damage_limit_with_damage_rate << "]"; } private: /// limit at which a point is considered as complitely broken Real damage_limit_with_damage_rate; /// internal pointer to the current damage vector const Vector * selected_damage_with_damage_rate; /// internal pointer to the current damage rate vector const Vector * selected_damage_rate_with_damage_rate; }; /* -------------------------------------------------------------------------- */ /* Stress Based Weight */ /* -------------------------------------------------------------------------- */ template class StressBasedWeightFunction : public BaseWeightFunction { public: StressBasedWeightFunction(Material & material); void init(); virtual void updateInternals(__attribute__((unused)) const ByElementTypeReal & quadrature_points_coordinates) { updatePrincipalStress(_not_ghost); updatePrincipalStress(_ghost); }; void updatePrincipalStress(GhostType ghost_type); inline void updateQuadraturePointsCoordinates(ByElementTypeReal & quadrature_points_coordinates); inline void selectType(ElementType type1, GhostType ghost_type1, ElementType type2, GhostType ghost_type2); inline Real operator()(Real r, QuadraturePoint & q1, QuadraturePoint & q2); inline Real computeRhoSquare(Real r, types::RVector & eigs, types::RMatrix & eigenvects, types::RVector & x_s); /* ------------------------------------------------------------------------ */ bool parseParam(const std::string & key, const std::string & value) { std::stringstream sstr(value); if(key == "ft") { sstr >> ft; } else return BaseWeightFunction::parseParam(key, value); return true; } /* ------------------------------------------------------------------------ */ virtual void printself(std::ostream & stream) const { stream << "StressBasedWeightFunction [ft: " << ft << "]"; } private: Real ft; ByElementTypeReal stress_diag; Vector * selected_stress_diag; ByElementTypeReal stress_base; Vector * selected_stress_base; ByElementTypeReal quadrature_points_coordinates; Vector * selected_position_1; Vector * selected_position_2; ByElementTypeReal characteristic_size; Vector * selected_characteristic_size; }; template inline std::ostream & operator <<(std::ostream & stream, const BaseWeightFunction & _this) { _this.printself(stream); return stream; } template inline std::ostream & operator <<(std::ostream & stream, const BaseWeightFunction * _this) { _this->printself(stream); return stream; } template inline std::ostream & operator >>(std::ostream & stream, __attribute__((unused)) const BaseWeightFunction * _this) { AKANTU_DEBUG_TO_IMPLEMENT(); return stream; } template inline std::ostream & operator >>(std::ostream & stream, __attribute__((unused)) const RemoveDamagedWeightFunction * _this) { AKANTU_DEBUG_TO_IMPLEMENT(); return stream; } template inline std::ostream & operator >>(std::ostream & stream, __attribute__((unused)) const DamagedWeightFunction * _this) { AKANTU_DEBUG_TO_IMPLEMENT(); return stream; } template inline std::ostream & operator >>(std::ostream & stream, __attribute__((unused)) const StressBasedWeightFunction * _this) { AKANTU_DEBUG_TO_IMPLEMENT(); return stream; } #include "weight_function_tmpl.hh" __END_AKANTU__ #endif /* __AKANTU_WEIGHT_FUNCTION_HH__ */ diff --git a/test/test_model/test_heat_transfer_model/test_heat_transfer_model_square2d.cc b/test/test_model/test_heat_transfer_model/test_heat_transfer_model_square2d.cc index e2432631f..498be9103 100644 --- a/test/test_model/test_heat_transfer_model/test_heat_transfer_model_square2d.cc +++ b/test/test_model/test_heat_transfer_model/test_heat_transfer_model_square2d.cc @@ -1,168 +1,169 @@ /** * @file test_heat_transfer_model_square2d.cc * * * @date Sun May 01 19:14:43 2011 * * @brief test of the class HeatTransferModel on the 3d cube * * @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 "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "heat_transfer_model.hh" #include "pbc_synchronizer.hh" /* -------------------------------------------------------------------------- */ #include #include #include using namespace std; /* -------------------------------------------------------------------------- */ - -#include "io_helper.hh" +#ifdef AKANTU_USE_IOHELPER +# include "dumper_paraview.hh" +#endif //AKANTU_USE_IOHELPER static void paraviewInit(akantu::HeatTransferModel & model, iohelper::Dumper & dumper); static void paraviewDump(iohelper::Dumper & dumper); iohelper::ElemType paraview_type = iohelper::TRIANGLE1; /* -------------------------------------------------------------------------- */ akantu::UInt spatial_dimension = 2; akantu:: ElementType type = akantu::_triangle_3; /* -------------------------------------------------------------------------- */ std::string base_name; int main(int argc, char *argv[]) { akantu::debug::setDebugLevel(akantu::dblWarning); akantu::initialize(argc, argv); //create mesh akantu::Mesh mesh(spatial_dimension); akantu::MeshIOMSH mesh_io; mesh_io.read("square_tri3.msh", mesh); akantu::HeatTransferModel model(mesh); //initialize everything model.initFull("material.dat"); //assemble the lumped capacity model.assembleCapacityLumped(); //get stable time step akantu::Real time_step = model.getStableTimeStep()*0.8; cout<<"time step is:" << time_step << endl; model.setTimeStep(time_step); //boundary conditions const akantu::Vector & nodes = model.getFEM().getMesh().getNodes(); akantu::Vector & boundary = model.getBoundary(); akantu::Vector & temperature = model.getTemperature(); double length; length = 1.; akantu::UInt nb_nodes = model.getFEM().getMesh().getNbNodes(); for (akantu::UInt i = 0; i < nb_nodes; ++i) { temperature(i) = 100.; akantu::Real dx = nodes(i,0) - length/4.; akantu::Real dy = 0.0; akantu::Real dz = 0.0; if (spatial_dimension > 1) dy = nodes(i,1) - length/4.; if (spatial_dimension == 3) dz = nodes(i,2) - length/4.; akantu::Real d = sqrt(dx*dx + dy*dy + dz*dz); // if(dx < 0.0){ if(d < 0.1){ boundary(i) = true; temperature(i) = 300.; } } iohelper::DumperParaview dumper; paraviewInit(model,dumper); //main loop int max_steps = 1000; for(int i=0; i