diff --git a/CMakeLists.txt b/CMakeLists.txt index 2e16fd30b..9c58be97a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,280 +1,280 @@ #=============================================================================== # @file CMakeLists.txt # @author Nicolas Richart # @date Fri Jun 11 09:46:59 2010 # # @section LICENSE # # # # @section DESCRIPTION # #=============================================================================== #=============================================================================== # _ ,-, _ # ,--, /: :\/': :`\/: :\ # |`; ' `,' `.; `: | _ _ # | | | ' | |. | | | | # | : | F E | A S | T++ || __ _| | ____ _ _ __ | |_ _ _ # | :. | : | : | : | \ / _` | |/ / _` | '_ \| __| | | | # \__/: :.. : :.. | :.. | ) | (_| | < (_| | | | | |_| |_| | # `---',\___/,\___/ /' \__,_|_|\_\__,_|_| |_|\__|\__,_| # `==._ .. . /' # `-::-' #=============================================================================== #=============================================================================== # CMake Project #=============================================================================== cmake_minimum_required(VERSION 2.6) project(AKANTU) enable_language(CXX) #=============================================================================== # Misc. #=============================================================================== set(AKANTU_CMAKE_DIR "${AKANTU_SOURCE_DIR}/cmake") set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake") set(BUILD_SHARED_LIBS ON CACHE BOOL "Build shared libraries.") #=============================================================================== # Version Number #=============================================================================== # AKANTU version number. An even minor number corresponds to releases. set(AKANTU_MAJOR_VERSION 0) set(AKANTU_MINOR_VERSION 1) set(AKANTU_BUILD_VERSION 0) set(AKANTU_VERSION "${AKANTU_MAJOR_VERSION}.${AKANTU_MINOR_VERSION}.${AKANTU_BUILD_VERSION}" ) # Append the library version information to the library target properties if(NOT AKANTU_NO_LIBRARY_VERSION) set(AKANTU_LIBRARY_PROPERTIES ${AKANTU_LIBRARY_PROPERTIES} VERSION "${AKANTU_VERSION}" SOVERSION "${AKANTU_MAJOR_VERSION}.${AKANTU_MINOR_VERSION}" ) endif(NOT AKANTU_NO_LIBRARY_VERSION) #=============================================================================== # Options #=============================================================================== option(AKANTU_DEBUG "Compiles akantu with the debug messages" ON) option(AKANTU_DOCUMENTATION "Build source documentation using Doxygen." OFF) # compilation switch option(AKANTU_BUILD_CONTACT "Build the contact algorithm" OFF) # features option(AKANTU_USE_IOHELPER "Add IOHelper support in akantu" OFF) option(AKANTU_USE_MPI "Add MPI support in akantu" OFF) option(AKANTU_USE_SCOTCH "Add Scotch support in akantu" OFF) option(AKANTU_USE_BLAS "Use BLAS for arithmetic operations" OFF) #=============================================================================== # Options handling #=============================================================================== # Debug if(NOT AKANTU_DEBUG) add_definitions(-DAKANTU_NDEBUG) endif(NOT AKANTU_DEBUG) # IOHelper if(AKANTU_USE_IOHELPER) find_package(IOHelper REQUIRED) if(IOHELPER_FOUND) add_definitions(-DAKANTU_USE_IOHELPER) set(AKANTU_EXTERNAL_LIB_INCLUDE_PATH ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH} ${IOHELPER_INCLUDE_PATH} ) set(AKANTU_EXTERNAL_LIBRARIES ${AKANTU_EXTERNAL_LIBRARIES} ${IOHELPER_LIBRARIES} ) endif(IOHELPER_FOUND) endif(AKANTU_USE_IOHELPER) # MPI set(AKANTU_MPI_ON FALSE) if(AKANTU_USE_MPI) find_package(MPI REQUIRED) if(MPI_FOUND) add_definitions(-DAKANTU_USE_MPI) set(AKANTU_EXTERNAL_LIB_INCLUDE_PATH ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH} ${MPI_INCLUDE_PATH} ) set(AKANTU_EXTERNAL_LIBRARIES ${AKANTU_EXTERNAL_LIBRARIES} ${MPI_LIBRARY} ${MPI_EXTRA_LIBRARY} ) set(AKANTU_MPI_ON TRUE) endif(MPI_FOUND) endif(AKANTU_USE_MPI) # Scotch set(AKANTU_SCOTCH_ON FALSE) if(AKANTU_USE_SCOTCH) find_package(Scotch REQUIRED) if(SCOTCH_FOUND) add_definitions(-DAKANTU_USE_SCOTCH) set(AKANTU_EXTERNAL_LIB_INCLUDE_PATH ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH} ${SCOTCH_INCLUDE_PATH} ) set(AKANTU_EXTERNAL_LIBRARIES ${AKANTU_EXTERNAL_LIBRARIES} ${SCOTCH_LIBRARIES} ) set(AKANTU_SCOTCH_ON TRUE) endif(SCOTCH_FOUND) endif(AKANTU_USE_SCOTCH) # BLAS set(AKANTU_BLAS_ON FALSE) if(AKANTU_USE_BLAS) enable_language(Fortran) find_package(BLAS) if(BLAS_FOUND) add_definitions(-DAKANTU_USE_CBLAS) set(AKANTU_EXTERNAL_LIBRARIES ${AKANTU_EXTERNAL_LIBRARIES} ${BLAS_LIBRARIES} ) set(AKANTU_BLAS_ON TRUE) endif(BLAS_FOUND) endif(AKANTU_USE_BLAS) #=============================================================================== # Library #=============================================================================== set(AKANTU_COMMON_SRC common/aka_common.cc common/aka_error.cc common/aka_extern.cc common/aka_static_memory.cc common/aka_memory.cc common/aka_vector.cc common/aka_math.cc fem/mesh.cc fem/fem.cc fem/element_class.cc # model/integration_scheme/central_difference.cc model/solid_mechanics/solid_mechanics_model.cc model/solid_mechanics/solid_mechanics_model_mass.cc model/solid_mechanics/solid_mechanics_model_boundary.cc model/solid_mechanics/solid_mechanics_model_material.cc model/solid_mechanics/material.cc model/solid_mechanics/materials/material_elastic.cc mesh_utils/mesh_io.cc mesh_utils/mesh_io/mesh_io_msh.cc mesh_utils/mesh_partition.cc mesh_utils/mesh_utils.cc - solver/sparse_matrix.cc +# solver/sparse_matrix.cc synchronizer/ghost_synchronizer.cc synchronizer/synchronizer.cc synchronizer/communicator.cc synchronizer/static_communicator.cc ) set(AKANTU_CONTACT_SRC model/solid_mechanics/contact.cc model/solid_mechanics/contact_search.cc model/solid_mechanics/contact_neighbor_structure.cc model/solid_mechanics/contact/contact_2d_explicit.cc model/solid_mechanics/contact/contact_search_2d_explicit.cc model/solid_mechanics/contact/regular_grid_neighbor_structure.cc model/solid_mechanics/contact/contact_search_3d_explicit.cc model/solid_mechanics/contact/contact_3d_explicit.cc ) if(AKANTU_BUILD_CONTACT) set(AKANTU_COMMON_SRC ${AKANTU_COMMON_SRC} ${AKANTU_CONTACT_SRC}) endif(AKANTU_BUILD_CONTACT) if(AKANTU_USE_MPI AND MPI_FOUND) set(AKANTU_COMMON_SRC ${AKANTU_COMMON_SRC} synchronizer/static_communicator_mpi.cc ) endif(AKANTU_USE_MPI AND MPI_FOUND) if(AKANTU_SCOTCH_ON) set(AKANTU_COMMON_SRC ${AKANTU_COMMON_SRC} mesh_utils/mesh_partition/mesh_partition_scotch.cc ) endif(AKANTU_SCOTCH_ON) set(AKANTU_INCLUDE_DIRS common fem/ mesh_utils/ mesh_utils/mesh_io/ mesh_utils/mesh_partition/ model/ model/integration_scheme model/solid_mechanics model/solid_mechanics/materials model/solid_mechanics/contact synchronizer/ solver/ ) include_directories( ${AKANTU_INCLUDE_DIRS} ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH} ) add_library(akantu ${AKANTU_COMMON_SRC}) set_target_properties(akantu PROPERTIES ${AKANTU_LIBRARY_PROPERTIES}) #=========================================================================== # Documentation #=========================================================================== if(AKANTU_DOCUMENTATION) find_package(Doxygen REQUIRED) if(DOXYGEN_FOUND) set(DOXYGEN_WARNINGS NO) set(DOXYGEN_QUIET YES) if(CMAKE_VERBOSE_MAKEFILE) set(DOXYGEN_WARNINGS YES) set(DOXYGEN_QUIET NO) endif(CMAKE_VERBOSE_MAKEFILE) add_subdirectory(doc/doxygen) endif(DOXYGEN_FOUND) endif(AKANTU_DOCUMENTATION) #=============================================================================== # Tests #=============================================================================== option(AKANTU_TESTS "Activate tests" OFF) if(AKANTU_TESTS) find_package(GMSH REQUIRED) add_subdirectory(test) endif(AKANTU_TESTS) #=============================================================================== # Examples #=============================================================================== option(AKANTU_EXAMPLES "Activate examples" OFF) if(AKANTU_EXAMPLES) find_package(GMSH REQUIRED) add_subdirectory(examples) endif(AKANTU_EXAMPLES) diff --git a/common/aka_common.hh b/common/aka_common.hh index 3f51aa212..0e5a6d2cb 100644 --- a/common/aka_common.hh +++ b/common/aka_common.hh @@ -1,267 +1,274 @@ /** * @file aka_common.hh * @author Nicolas Richart * @date Fri Jun 11 09:48:06 2010 * * @namespace akantu * * @section LICENSE * * \ * * @section DESCRIPTION * * All common things to be included in the projects files * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COMMON_HH__ #define __AKANTU_COMMON_HH__ /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include #include #include #include #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ #define __BEGIN_AKANTU__ namespace akantu { #define __END_AKANTU__ } /* -------------------------------------------------------------------------- */ #include "aka_error.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* Common types */ /* -------------------------------------------------------------------------- */ typedef double Real; typedef unsigned int UInt; typedef int Int; typedef std::string ID; #ifdef AKANTU_NDEBUG static const Real REAL_INIT_VALUE = 0; #else static const Real REAL_INIT_VALUE = std::numeric_limits::quiet_NaN(); #endif /* -------------------------------------------------------------------------- */ /* Memory types */ /* -------------------------------------------------------------------------- */ typedef UInt MemoryID; typedef ID VectorID; /* -------------------------------------------------------------------------- */ /* Mesh/FEM/Model types */ /* -------------------------------------------------------------------------- */ typedef ID MeshID; typedef ID FEMID; typedef ID ModelID; typedef ID MaterialID; typedef ID SparseMatrixID; typedef UInt Surface; /// @enum ElementType type of element potentially contained in a Mesh enum ElementType { _not_defined = 0, _segment_2 = 1, /// first order segment _segment_3 = 2, /// second order segment _triangle_3 = 3, /// first order triangle _triangle_6 = 4, /// second order triangle _tetrahedron_4 = 5, /// first order tetrahedron _tetrahedron_10 = 6, /// second order tetrahedron @remark not implemented yet _max_element_type, _point /// point only for some algorithm to be generic like mesh partitioning }; /// @enum MaterialType different materials implemented enum MaterialType { _elastic = 0, _max_material_type }; + +/// @enum BoundaryFunctionType type of function passed for boundary conditions +enum BoundaryFunctionType { + _bft_forces, + _bft_stress +}; + /// @enum SparseMatrixType type of sparse matrix used enum SparseMatrixType { _unsymmetric, _symmetric }; /* -------------------------------------------------------------------------- */ /* Contact */ /* -------------------------------------------------------------------------- */ typedef ID ContactID; typedef ID ContactSearchID; typedef ID ContactNeighborStructureID; enum ContactType { _ct_not_defined = 0, _ct_2d_expli = 1, _ct_3d_expli = 2 }; enum ContactSearchType { _cst_not_defined = 0, _cst_2d_expli = 1, _cst_3d_expli = 2 }; enum ContactNeighborStructureType { _cnst_not_defined = 0, _cnst_regular_grid = 1 }; /* -------------------------------------------------------------------------- */ /* Ghosts handling */ /* -------------------------------------------------------------------------- */ typedef ID SynchronizerID; /// @enum CommunicatorType type of communication method to use enum CommunicatorType { _communicator_mpi, _communicator_dummy }; /// @enum GhostSynchronizationTag type of synchronizations enum GhostSynchronizationTag { /// SolidMechanicsModel tags _gst_smm_mass, /// synchronization of the SolidMechanicsModel.mass _gst_smm_residual, /// synchronization of the SolidMechanicsModel.current_position _gst_smm_boundary, /// synchronization of the boundary, forces, velocities and displacement /// Test tag _gst_test }; /// @enum GhostType type of ghost enum GhostType { _not_ghost, _ghost, _casper // not used but a real cute ghost }; /// @enum SynchronizerOperation reduce operation that the synchronizer can perform enum SynchronizerOperation { _so_sum, _so_min, _so_max, _so_null }; /* -------------------------------------------------------------------------- */ /* Global defines */ /* -------------------------------------------------------------------------- */ #define AKANTU_MIN_ALLOCATION 2000 #define AKANTU_INDENT " " /* -------------------------------------------------------------------------- */ #define AKANTU_SET_MACRO(name, variable, type) \ inline void set##name (type variable) { \ this->variable = variable; \ } #define AKANTU_GET_MACRO(name, variable, type) \ inline type get##name () const { \ return variable; \ } #define AKANTU_GET_MACRO_NOT_CONST(name, variable, type) \ inline type get##name () { \ return variable; \ } #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE(name, variable, type) \ inline type get##name (const ::akantu::ElementType & el_type) const { \ AKANTU_DEBUG_IN(); \ AKANTU_DEBUG_ASSERT(variable[el_type] != NULL, \ "No " << #variable << " of type " \ << el_type << " in " << id); \ AKANTU_DEBUG_OUT(); \ return *variable[el_type]; \ } /* -------------------------------------------------------------------------- */ //! standard output stream operator for ElementType inline std::ostream & operator <<(std::ostream & stream, ElementType type) { switch(type) { case _segment_2 : stream << "segment_2" ; break; case _segment_3 : stream << "segment_3" ; break; case _triangle_3 : stream << "triangle_3" ; break; case _triangle_6 : stream << "triangle_6" ; break; case _tetrahedron_4 : stream << "tetrahedron_4"; break; case _tetrahedron_10 : stream << "tetrahedron_10"; break; case _not_defined : stream << "undefined" ; break; case _max_element_type : stream << "ElementType(" << (int) type << ")"; break; case _point : stream << "point"; break; } return stream; } /* -------------------------------------------------------------------------- */ void initialize(int * argc, char *** argv); /* -------------------------------------------------------------------------- */ void finalize (); /* -------------------------------------------------------------------------- */ /* * For intel compiler annoying remark */ #if defined(__INTEL_COMPILER) /// remark #981: operands are evaluated in unspecified order #pragma warning ( disable : 981 ) /// remark #383: value copied to temporary, reference to temporary used #pragma warning ( disable : 383 ) #endif //defined(__INTEL_COMPILER) /* -------------------------------------------------------------------------- */ /* string manipulation */ /* -------------------------------------------------------------------------- */ inline void to_lower(std::string & str) { std::transform(str.begin(), str.end(), str.begin(), (int(*)(int))std::tolower); } /* -------------------------------------------------------------------------- */ inline void trim(std::string & to_trim) { size_t first = to_trim.find_first_not_of(" \t"); if (first != std::string::npos) { size_t last = to_trim.find_last_not_of(" \t"); to_trim = to_trim.substr(first, last - first + 1); } else to_trim = ""; } __END_AKANTU__ #endif /* __AKANTU_COMMON_HH__ */ diff --git a/examples/3d_cube_compression/3d_cube_compression.cc b/examples/3d_cube_compression/3d_cube_compression.cc index f986362a4..a31d9ec4c 100644 --- a/examples/3d_cube_compression/3d_cube_compression.cc +++ b/examples/3d_cube_compression/3d_cube_compression.cc @@ -1,146 +1,147 @@ /** * @file 3d_cube_compression.cc * @author Nicolas Richart * @date Wed Nov 24 23:15:47 2010 * * @brief 3d cube compression * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" +#include "mesh_io_msh.hh" #include "solid_mechanics_model.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER int main(int argc, char *argv[]) { akantu::ElementType type = akantu::_tetrahedron_4; #ifdef AKANTU_USE_IOHELPER akantu::UInt paraview_type = TETRA1; #endif //AKANTU_USE_IOHELPER akantu::UInt spatial_dimension = 3; akantu::UInt max_steps = 5000; akantu::Real time_step = 1e-6; // akantu::Real epot, ekin; akantu::Mesh mesh(spatial_dimension); akantu::MeshIOMSH mesh_io; mesh_io.read("cube.msh", mesh); akantu::SolidMechanicsModel * model = new akantu::SolidMechanicsModel(mesh); akantu::UInt nb_nodes = model->getFEM().getMesh().getNbNodes(); akantu::UInt nb_element = model->getFEM().getMesh().getNbElement(type); /// model initialization model->initVectors(); /// set vectors to 0 memset(model->getForce().values, 0, spatial_dimension*nb_nodes*sizeof(akantu::Real)); memset(model->getVelocity().values, 0, spatial_dimension*nb_nodes*sizeof(akantu::Real)); memset(model->getAcceleration().values, 0, spatial_dimension*nb_nodes*sizeof(akantu::Real)); memset(model->getDisplacement().values, 0, spatial_dimension*nb_nodes*sizeof(akantu::Real)); model->readMaterials("material.dat"); model->initMaterials(); model->initModel(); std::cout << model->getMaterial(0) << std::endl; model->assembleMassLumped(); #ifdef AKANTU_USE_IOHELPER /// set to 0 only for the first paraview dump memset(model->getResidual().values, 0, spatial_dimension*nb_nodes*sizeof(akantu::Real)); memset(model->getMaterial(0).getStrain(type).values, 0, spatial_dimension*spatial_dimension*nb_element*sizeof(akantu::Real)); memset(model->getMaterial(0).getStress(type).values, 0, spatial_dimension*spatial_dimension*nb_element*sizeof(akantu::Real)); #endif //AKANTU_USE_IOHELPER /// boundary conditions akantu::Real eps = 1e-16; for (akantu::UInt i = 0; i < nb_nodes; ++i) { // if(model->getFEM().getMesh().getNodes().values[spatial_dimension*i] >= 9) // model->getDisplacement().values[spatial_dimension*i] = (model->getFEM().getMesh().getNodes().values[spatial_dimension*i] - 9) / 100. ; if(fabs(model->getFEM().getMesh().getNodes().values[spatial_dimension*i + 2] - 1) <= eps) model->getForce().values[spatial_dimension*i + 2] = -250; if(model->getFEM().getMesh().getNodes().values[spatial_dimension*i + 2] <= eps) model->getBoundary().values[spatial_dimension*i + 2] = true; if(model->getFEM().getMesh().getNodes().values[spatial_dimension*i + 0] <= eps) { model->getBoundary().values[spatial_dimension*i + 0] = true; } } // akantu::Real time_step = model->getStableTimeStep() * time_factor; std::cout << "Time Step = " << time_step << "s" << std::endl; model->setTimeStep(time_step); // model->setTimeStep(3.54379e-07); #ifdef AKANTU_USE_IOHELPER DumperParaview dumper; dumper.SetMode(TEXT); dumper.SetPoints(model->getFEM().getMesh().getNodes().values, spatial_dimension, nb_nodes, "coordinates"); dumper.SetConnectivity((int *)model->getFEM().getMesh().getConnectivity(type).values, paraview_type, nb_element, C_MODE); dumper.AddNodeDataField(model->getDisplacement().values, spatial_dimension, "displacements"); dumper.AddNodeDataField(model->getVelocity().values, spatial_dimension, "velocity"); dumper.AddNodeDataField(model->getResidual().values, spatial_dimension, "force"); dumper.AddElemDataField(model->getMaterial(0).getStrain(type).values, spatial_dimension*spatial_dimension, "strain"); dumper.AddElemDataField(model->getMaterial(0).getStress(type).values, spatial_dimension*spatial_dimension, "stress"); dumper.SetEmbeddedValue("displacements", 1); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); #endif //AKANTU_USE_IOHELPER for(akantu::UInt s = 1; s <= max_steps; ++s) { model->explicitPred(); model->updateResidual(); model->updateAcceleration(); model->explicitCorr(); #ifdef AKANTU_USE_IOHELPER if(s % 1 == 0) dumper.Dump(); #endif //AKANTU_USE_IOHELPER if(s % 10 == 0) std::cout << "passing step " << s << "/" << max_steps << std::endl; } return EXIT_SUCCESS; } diff --git a/model/solid_mechanics/solid_mechanics_model.hh b/model/solid_mechanics/solid_mechanics_model.hh index 96256ea11..ea28c4efd 100644 --- a/model/solid_mechanics/solid_mechanics_model.hh +++ b/model/solid_mechanics/solid_mechanics_model.hh @@ -1,256 +1,256 @@ /** * @file solid_mechanics_model.hh * @author Nicolas Richart * @date[creation] Thu Jul 22 11:51:06 2010 * @date[last modification] Thu Oct 14 14:00:06 2010 * * @brief Model of Solid Mechanics * * @section LICENSE * * \ * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_HH__ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "model.hh" #include "material.hh" #include "material_parser.hh" /* -------------------------------------------------------------------------- */ namespace akantu { // class Material; class IntegrationScheme2ndOrder; class Contact; } __BEGIN_AKANTU__ class SolidMechanicsModel : public Model { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: SolidMechanicsModel(UInt spatial_dimension, const ModelID & id = "solid_mechanics_model", const MemoryID & memory_id = 0); SolidMechanicsModel(Mesh & mesh, UInt spatial_dimension = 0, const ModelID & id = "solid_mechanics_model", const MemoryID & memory_id = 0); virtual ~SolidMechanicsModel(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// allocate all vectors void initVectors(); /// read the material files to instantiate all the materials void readMaterials(const std::string & filename); /// read a custom material with a keyword and class as template template UInt readCustomMaterial(const std::string & filename, const std::string & keyword); /// read properties part of a material file and create the material template Material * readMaterialProperties(std::ifstream & infile, MaterialID mat_id, UInt ¤t_line); /// initialize all internal arrays for materials void initMaterials(); /// initialize the model void initModel(); /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the residual for the explicit scheme void updateResidual(); /// compute the acceleration from the residual void updateAcceleration(); /// explicit integration predictor void explicitPred(); /// explicit integration corrector void explicitCorr(); /// synchronize the ghost element boundaries values void synchronizeBoundaries(); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /// integrate a force on the boundary by providing a stress tensor void computeForcesByStressTensor(const Vector & stresses, const ElementType & type); /// integrate a force on the boundary by providing a traction vector void computeForcesByTractionVector(const Vector & tractions, const ElementType & type); /// compute force vector from a function(x,y,z) that describe stresses - void computeForcesFromFunction(void (*myf)(double *,double *), UInt function_type); + void computeForcesFromFunction(void (*myf)(double *,double *), BoundaryFunctionType function_type); private: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumpedRowSum(GhostType ghost_type, const ElementType type); /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumpedDiagonalScaling(GhostType ghost_type, const ElementType type); /* ------------------------------------------------------------------------ */ /* Ghost Synchronizer inherited members */ /* ------------------------------------------------------------------------ */ public: inline virtual UInt getNbDataToPack(const Element & element, GhostSynchronizationTag tag) const; inline virtual UInt getNbDataToUnpack(const Element & element, GhostSynchronizationTag tag) const; inline virtual void packData(Real ** buffer, const Element & element, GhostSynchronizationTag tag) const; inline virtual void unpackData(Real ** buffer, const Element & element, GhostSynchronizationTag tag) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(TimeStep, time_step, Real); AKANTU_SET_MACRO(TimeStep, time_step, Real); AKANTU_GET_MACRO(F_M2A, f_m2a, Real); AKANTU_SET_MACRO(F_M2A, f_m2a, Real); AKANTU_GET_MACRO(Displacement, *displacement, Vector &); AKANTU_GET_MACRO(CurrentPosition, *current_position, const Vector &); AKANTU_GET_MACRO(Increment, *increment, const Vector &); AKANTU_GET_MACRO(Mass, *mass, const Vector &); AKANTU_GET_MACRO(Velocity, *velocity, Vector &); AKANTU_GET_MACRO(Acceleration, *acceleration, Vector &); AKANTU_GET_MACRO(Force, *force, Vector &); AKANTU_GET_MACRO(Residual, *residual, const Vector &); AKANTU_GET_MACRO(Boundary, *boundary, Vector &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementMaterial, element_material, Vector &); inline Material & getMaterial(UInt mat_index); /// compute the stable time step Real getStableTimeStep(); // void setPotentialEnergyFlagOn(); // void setPotentialEnergyFlagOff(); Real getPotentialEnergy(); Real getKineticEnergy(); AKANTU_SET_MACRO(Contact, contact, Contact *); void setIncrementFlagOn(); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// time step Real time_step; /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array Vector * displacement; /// lumped mass array Vector * mass; /// velocities array Vector * velocity; /// accelerations array Vector * acceleration; /// forces array Vector * force; /// residuals array Vector * residual; /// boundaries array Vector * boundary; /// array of current position used during update residual Vector * current_position; /// materials of all elements ByElementTypeUInt element_material; /// materials of all ghost elements ByElementTypeUInt ghost_element_material; /// list of used materials std::vector materials; /// integration scheme of second order used IntegrationScheme2ndOrder * integrator; /// increment of displacement Vector * increment; /// flag defining if the increment must be computed or not bool increment_flag; /// object to resolve the contact Contact * contact; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_inline_impl.cc" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const SolidMechanicsModel & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */ diff --git a/model/solid_mechanics/solid_mechanics_model_boundary.cc b/model/solid_mechanics/solid_mechanics_model_boundary.cc index 6601baa24..b5b9fad8b 100644 --- a/model/solid_mechanics/solid_mechanics_model_boundary.cc +++ b/model/solid_mechanics/solid_mechanics_model_boundary.cc @@ -1,166 +1,181 @@ /** * @file solid_mechanics_model_boundary.cc * @author Guillaume ANCIAUX * @date Fri Nov 19 10:23:03 2010 * * @brief * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ +/** + * @param myf pointer to a function that fills a vector/tensor with respect to passed coordinates + * @param function_type flag to specify the take of function: 0 is tensor like and 1 is traction like. + */ void SolidMechanicsModel::computeForcesFromFunction(void (*myf)(double *,double *), - UInt function_type = 0){ + BoundaryFunctionType function_type){ /** function type is ** 0 : stress function is given ** 1 : traction function is given */ std::stringstream name; name << id << ":solidmechanics:imposed_traction"; Vector traction_funct(0, spatial_dimension,name.str()); name.clear(); name << id << ":solidmechanics:imposed_stresses"; Vector stress_funct(0, spatial_dimension*spatial_dimension,name.str()); name.clear(); name << id << ":solidmechanics:quad_coords"; Vector quad_coords(0,spatial_dimension,"quad_coords"); + UInt offset = 0; + switch(function_type) { + case _bft_stress: + offset = spatial_dimension*spatial_dimension; break; + case _bft_forces: + offset = spatial_dimension; break; + } + + //prepare the loop over element types const Mesh::ConnectivityTypeList & type_list = fem_boundary->getMesh().getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != fem_boundary->getElementDimension()) continue; UInt nb_quad = FEM::getNbQuadraturePoints(*it); UInt nb_element = fem_boundary->getMesh().getNbElement(*it); fem_boundary->interpolateOnQuadraturePoints(fem_boundary->getMesh().getNodes(), quad_coords, spatial_dimension, (*it)); Real * imposed_val = NULL; - if (function_type == 0){ + switch(function_type) { + case _bft_stress: stress_funct.resize(nb_element*nb_quad); imposed_val = stress_funct.values; - } - else { + break; + case _bft_forces: traction_funct.resize(nb_element*nb_quad); imposed_val = traction_funct.values; + break; } /// sigma/tractions on each quadrature points for (UInt el = 0; el < nb_element; ++el) { Real * qcoord = quad_coords.values+el*nb_quad*spatial_dimension; for (UInt q = 0; q < nb_quad; ++q) { myf(qcoord+q,imposed_val); - if (function_type == 0) - imposed_val += spatial_dimension*spatial_dimension; - else - imposed_val += spatial_dimension; + imposed_val += offset; } } - if (function_type == 0) - computeForcesByStressTensor(stress_funct,(*it)); - else - computeForcesByTractionVector(traction_funct,(*it)); + + switch(function_type) { + case _bft_stress: + computeForcesByStressTensor(stress_funct,(*it)); break; + case _bft_forces: + computeForcesByTractionVector(traction_funct,(*it)); break; + } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeForcesByStressTensor(const Vector & stresses, const ElementType & type){ AKANTU_DEBUG_IN(); UInt nb_element = fem_boundary->getMesh().getNbElement(type); UInt nb_quad = FEM::getNbQuadraturePoints(type); // check dimension match AKANTU_DEBUG_ASSERT(Mesh::getSpatialDimension(type) == fem_boundary->getElementDimension(), "element type dimension does not match the dimension of boundaries : " << fem_boundary->getElementDimension() << " != " << Mesh::getSpatialDimension(type)); // check size of the vector AKANTU_DEBUG_ASSERT(stresses.getSize() == nb_quad*nb_element, "the size of the vector should be the total number of quadrature points"); // check number of components AKANTU_DEBUG_ASSERT(stresses.getNbComponent() == spatial_dimension*spatial_dimension, "the number of components should be the dimension of 2-tensors"); std::stringstream name; name << id << ":solidmechanics:" << type << ":traction_boundary"; Vector funct(nb_element*nb_quad, spatial_dimension,name.str()); const Vector & normals_on_quad = fem_boundary->getNormalsOnQuadPoints(type); Math::matrix_vector(spatial_dimension,spatial_dimension,stresses,normals_on_quad,funct); computeForcesByTractionVector(funct,type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeForcesByTractionVector(const Vector & tractions, const ElementType & type){ AKANTU_DEBUG_IN(); UInt nb_element = fem_boundary->getMesh().getNbElement(type); UInt nb_nodes_per_element = fem_boundary->getMesh().getNbNodesPerElement(type); UInt nb_quad = FEM::getNbQuadraturePoints(type); // check dimension match AKANTU_DEBUG_ASSERT(Mesh::getSpatialDimension(type) == fem_boundary->getElementDimension(), "element type dimension does not match the dimension of boundaries : " << fem_boundary->getElementDimension() << " != " << Mesh::getSpatialDimension(type)); // check size of the vector AKANTU_DEBUG_ASSERT(tractions.getSize() == nb_quad*nb_element, "the size of the vector should be the total number of quadrature points"); // check number of components AKANTU_DEBUG_ASSERT(tractions.getNbComponent() == spatial_dimension, "the number of components should be the spatial dimension of the problem"); // do a complete copy of the vector Vector funct(tractions,true); // extend the vector to multiply by the shapes (prepare to assembly) funct.extendComponentsInterlaced(nb_nodes_per_element,spatial_dimension); // multiply by the shapes Real * funct_val = funct.values; Real * shapes_val = (fem_boundary->getShapes(type)).values; for (UInt el = 0; el < nb_element; ++el) { for (UInt q = 0; q < nb_quad; ++q) { for (UInt n = 0; n < nb_nodes_per_element; ++n,++shapes_val) { for (UInt i = 0; i < spatial_dimension; ++i) { *funct_val++ *= *shapes_val; } } } } // allocate the vector that will contain the integrated values std::stringstream name; name << id << ":solidmechanics:" << type << ":integral_boundary"; Vector int_funct(nb_element, spatial_dimension*nb_nodes_per_element,name.str()); //do the integration fem_boundary->integrate(funct, int_funct, spatial_dimension*nb_nodes_per_element, type); // assemble the result into force vector fem_boundary->assembleVector(int_funct,const_cast &>(getForce()), spatial_dimension, type); AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_local_material.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_local_material.cc index 4a7b60171..3187efb07 100644 --- a/test/test_model/test_solid_mechanics_model/test_materials/test_local_material.cc +++ b/test/test_model/test_solid_mechanics_model/test_materials/test_local_material.cc @@ -1,124 +1,124 @@ /** * @file test_solid_mechanics_model.cc * @author Nicolas Richart * @date Tue Jul 27 14:34:13 2010 * * @brief test of the class SolidMechanicsModel * * @section LICENSE * * \ * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "fem.hh" #include "material_damage.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; static void trac(__attribute__ ((unused)) double * position,double * stress){ memset(stress,0,sizeof(Real)*4); stress[0] = 1000; stress[3] = 1000; } int main(int argc, char *argv[]) { UInt max_steps = 10000; Real epot, ekin; Mesh mesh(2); MeshIOMSH mesh_io; mesh_io.read("circle2.msh", mesh); SolidMechanicsModel * model = new SolidMechanicsModel(mesh); /// model initialization model->initVectors(); UInt nb_nodes = model->getFEM().getMesh().getNbNodes(); memset(model->getForce().values, 0, 2*nb_nodes*sizeof(Real)); memset(model->getVelocity().values, 0, 2*nb_nodes*sizeof(Real)); memset(model->getAcceleration().values, 0, 2*nb_nodes*sizeof(Real)); memset(model->getDisplacement().values, 0, 2*nb_nodes*sizeof(Real)); memset(model->getResidual().values, 0, 2*nb_nodes*sizeof(Real)); memset(model->getMass().values, 1, nb_nodes*sizeof(Real)); model->readCustomMaterial("material.dat","DAMAGE"); model->initMaterials(); model->initModel(); Real time_step = model->getStableTimeStep(); model->setTimeStep(time_step/10.); model->assembleMassLumped(); std::cout << *model << std::endl; FEM & fem_boundary = model->getFEMBoundary(); fem_boundary.initShapeFunctions(); fem_boundary.computeNormalsOnQuadPoints(); - model->computeForcesFromFunction(trac,0); + model->computeForcesFromFunction(trac, akantu::_bft_forces); #ifdef AKANTU_USE_IOHELPER model->updateResidual(); DumperParaview dumper; dumper.SetMode(BASE64); dumper.SetPoints(model->getFEM().getMesh().getNodes().values, 2, nb_nodes, "coordinates"); dumper.SetConnectivity((int *)model->getFEM().getMesh().getConnectivity(_triangle_6).values, TRIANGLE2, model->getFEM().getMesh().getNbElement(_triangle_6), C_MODE); dumper.AddNodeDataField(model->getDisplacement().values, 2, "displacements"); dumper.AddNodeDataField(model->getVelocity().values, 2, "velocity"); dumper.AddNodeDataField(model->getForce().values, 2, "force"); dumper.AddNodeDataField(model->getMass().values, 1, "Mass"); dumper.AddNodeDataField(model->getResidual().values, 2, "residual"); dumper.AddElemDataField(model->getMaterial(0).getStrain(_triangle_6).values, 4, "strain"); dumper.AddElemDataField(model->getMaterial(0).getStress(_triangle_6).values, 4, "stress"); MaterialDamage & mat = dynamic_cast(model->getMaterial(0)); Real * dam = mat.getDamage(_triangle_6).values; dumper.AddElemDataField(dam, 1, "damage"); dumper.SetEmbeddedValue("displacements", 1); dumper.SetEmbeddedValue("force", 1); dumper.SetEmbeddedValue("residual", 1); dumper.SetEmbeddedValue("velocity", 1); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); #endif //AKANTU_USE_IOHELPER for(UInt s = 0; s < max_steps; ++s) { model->explicitPred(); model->updateResidual(); model->updateAcceleration(); model->explicitCorr(); epot = model->getPotentialEnergy(); ekin = model->getKineticEnergy(); std::cout << s << " " << epot << " " << ekin << " " << epot + ekin << std::endl; #ifdef AKANTU_USE_IOHELPER if(s % 100 == 0) dumper.Dump(); #endif //AKANTU_USE_IOHELPER } return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_bar_traction2d.cc b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_bar_traction2d.cc index 932a73d88..fea9edf9b 100644 --- a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_bar_traction2d.cc +++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_bar_traction2d.cc @@ -1,224 +1,224 @@ /** * @file test_solid_mechanics_model.cc * @author Nicolas Richart * @date Tue Jul 27 14:34:13 2010 * * @brief test of the class SolidMechanicsModel * * @section LICENSE * * \ * */ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "solid_mechanics_model.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER #define CHECK_STRESS static void trac(double * position,double * traction){ memset(traction,0,sizeof(akantu::Real)*2); if (fabs(position[0] - 10) < 1e-2){ traction[0] = 1000; traction[1] = 1000; } } int main(int argc, char *argv[]) { akantu::ElementType type = akantu::_triangle_3; #ifdef AKANTU_USE_IOHELPER akantu::UInt paraview_type = TRIANGLE1; #endif //AKANTU_USE_IOHELPER akantu::UInt spatial_dimension = 2; akantu::UInt max_steps = 10000; akantu::Real time_factor = 0.2; // akantu::Real epot, ekin; akantu::Mesh mesh(spatial_dimension); akantu::MeshIOMSH mesh_io; mesh_io.read("bar1.msh", mesh); akantu::SolidMechanicsModel * model = new akantu::SolidMechanicsModel(mesh); akantu::UInt nb_nodes = model->getFEM().getMesh().getNbNodes(); akantu::UInt nb_element = model->getFEM().getMesh().getNbElement(type); /// model initialization model->initVectors(); /// set vectors to 0 memset(model->getForce().values, 0, spatial_dimension*nb_nodes*sizeof(akantu::Real)); memset(model->getVelocity().values, 0, spatial_dimension*nb_nodes*sizeof(akantu::Real)); memset(model->getAcceleration().values, 0, spatial_dimension*nb_nodes*sizeof(akantu::Real)); memset(model->getDisplacement().values, 0, spatial_dimension*nb_nodes*sizeof(akantu::Real)); model->readMaterials("material.dat"); model->initMaterials(); model->initModel(); std::cout << model->getMaterial(0) << std::endl; model->assembleMassLumped(); #ifdef AKANTU_USE_IOHELPER /// set to 0 only for the first paraview dump memset(model->getResidual().values, 0, spatial_dimension*nb_nodes*sizeof(akantu::Real)); memset(model->getMaterial(0).getStrain(type).values, 0, spatial_dimension*spatial_dimension*nb_element*sizeof(akantu::Real)); memset(model->getMaterial(0).getStress(type).values, 0, spatial_dimension*spatial_dimension*nb_element*sizeof(akantu::Real)); #endif //AKANTU_USE_IOHELPER /// boundary conditions akantu::Real eps = 1e-16; for (akantu::UInt i = 0; i < nb_nodes; ++i) { // if(model->getFEM().getMesh().getNodes().values[spatial_dimension*i] >= 9) // model->getDisplacement().values[spatial_dimension*i] = (model->getFEM().getMesh().getNodes().values[spatial_dimension*i] - 9) / 100. ; if(model->getFEM().getMesh().getNodes().values[spatial_dimension*i] <= eps) model->getBoundary().values[spatial_dimension*i] = true; if(model->getFEM().getMesh().getNodes().values[spatial_dimension*i + 1] <= eps || model->getFEM().getMesh().getNodes().values[spatial_dimension*i + 1] >= 1 - eps ) { model->getBoundary().values[spatial_dimension*i + 1] = true; } } akantu::FEM & fem_boundary = model->getFEMBoundary(); fem_boundary.initShapeFunctions(); fem_boundary.computeNormalsOnQuadPoints(); - model->computeForcesFromFunction(trac,0); + model->computeForcesFromFunction(trac, akantu::_bft_forces); akantu::Real time_step = model->getStableTimeStep() * time_factor; std::cout << "Time Step = " << time_step << "s" << std::endl; model->setTimeStep(time_step); // model->setTimeStep(3.54379e-07); #ifdef AKANTU_USE_IOHELPER DumperParaview dumper; dumper.SetMode(TEXT); dumper.SetPoints(model->getFEM().getMesh().getNodes().values, spatial_dimension, nb_nodes, "coordinates"); dumper.SetConnectivity((int *)model->getFEM().getMesh().getConnectivity(type).values, paraview_type, nb_element, C_MODE); dumper.AddNodeDataField(model->getDisplacement().values, spatial_dimension, "displacements"); dumper.AddNodeDataField(model->getVelocity().values, spatial_dimension, "velocity"); dumper.AddNodeDataField(model->getResidual().values, spatial_dimension, "force"); dumper.AddNodeDataField(model->getForce().values, spatial_dimension, "applied_force"); dumper.AddElemDataField(model->getMaterial(0).getStrain(type).values, spatial_dimension*spatial_dimension, "strain"); dumper.AddElemDataField(model->getMaterial(0).getStress(type).values, spatial_dimension*spatial_dimension, "stress"); dumper.SetEmbeddedValue("displacements", 1); dumper.SetEmbeddedValue("applied_force", 1); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); #endif //AKANTU_USE_IOHELPER #ifdef CHECK_STRESS std::ofstream outfile; outfile.open("stress"); #endif // CHECK_STRESS for(akantu::UInt s = 1; s <= max_steps; ++s) { model->explicitPred(); model->updateResidual(); model->updateAcceleration(); model->explicitCorr(); // epot = model->getPotentialEnergy(); // ekin = model->getKineticEnergy(); // std::cout << s << " " << epot << " " << ekin << " " << epot + ekin // << std::endl; #ifdef CHECK_STRESS akantu::Real max_stress = std::numeric_limits::min(); //akantu::UInt max_el = 0; akantu::Real * stress = model->getMaterial(0).getStress(type).values; for (akantu::UInt i = 0; i < nb_element; ++i) { if(max_stress < stress[i*spatial_dimension*spatial_dimension]) { max_stress = stress[i*spatial_dimension*spatial_dimension]; // max_el = i; } } akantu::Real * coord = model->getFEM().getMesh().getNodes().values; akantu::Real * disp_val = model->getDisplacement().values; akantu::UInt * conn = model->getFEM().getMesh().getConnectivity(type).values; akantu::UInt nb_nodes_per_element = model->getFEM().getMesh().getNbNodesPerElement(type); akantu::Real * coords = new akantu::Real[spatial_dimension]; akantu::Real min_x = std::numeric_limits::max(); akantu::Real max_x = std::numeric_limits::min(); akantu::Real stress_range = 5e7; for (akantu::UInt el = 0; el < nb_element; ++el) { if(stress[el*spatial_dimension*spatial_dimension] > max_stress - stress_range) { akantu::UInt el_offset = el * nb_nodes_per_element; memset(coords, 0, spatial_dimension*sizeof(akantu::Real)); for (akantu::UInt n = 0; n < nb_nodes_per_element; ++n) { for (akantu::UInt i = 0; i < spatial_dimension; ++i) { akantu::UInt node = conn[el_offset + n] * spatial_dimension; coords[i] += (coord[node + i] + disp_val[node + i]) / ((akantu::Real) nb_nodes_per_element); } } min_x = min_x < coords[0] ? min_x : coords[0]; max_x = max_x > coords[0] ? max_x : coords[0]; } } outfile << s << " " << .5 * (min_x + max_x) << " " << min_x << " " << max_x << " " << max_x - min_x << " " << max_stress << std::endl; delete [] coords; #endif // CHECK_STRESS #ifdef AKANTU_USE_IOHELPER if(s % 100 == 0) dumper.Dump(); #endif //AKANTU_USE_IOHELPER if(s % 10 == 0) std::cout << "passing step " << s << "/" << max_steps << std::endl; } #ifdef CHECK_STRESS outfile.close(); #endif // CHECK_STRESS return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_circle_2.cc b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_circle_2.cc index 6edfefeea..92132abc2 100644 --- a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_circle_2.cc +++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_circle_2.cc @@ -1,117 +1,117 @@ /** * @file test_solid_mechanics_model.cc * @author Nicolas Richart * @date Tue Jul 27 14:34:13 2010 * * @brief test of the class SolidMechanicsModel * * @section LICENSE * * \ * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "fem.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; static void trac(__attribute__ ((unused)) double * position,double * stress){ memset(stress,0,sizeof(Real)*4); stress[0] = 1000; stress[3] = 1000; } int main(int argc, char *argv[]) { UInt max_steps = 10000; Real epot, ekin; Mesh mesh(2); MeshIOMSH mesh_io; mesh_io.read("circle2.msh", mesh); SolidMechanicsModel * model = new SolidMechanicsModel(mesh); /// model initialization model->initVectors(); UInt nb_nodes = model->getFEM().getMesh().getNbNodes(); memset(model->getForce().values, 0, 2*nb_nodes*sizeof(Real)); memset(model->getVelocity().values, 0, 2*nb_nodes*sizeof(Real)); memset(model->getAcceleration().values, 0, 2*nb_nodes*sizeof(Real)); memset(model->getDisplacement().values, 0, 2*nb_nodes*sizeof(Real)); memset(model->getResidual().values, 0, 2*nb_nodes*sizeof(Real)); memset(model->getMass().values, 1, nb_nodes*sizeof(Real)); model->readMaterials("material.dat"); model->initMaterials(); model->initModel(); Real time_step = model->getStableTimeStep(); model->setTimeStep(time_step/10.); model->assembleMassLumped(); std::cout << *model << std::endl; FEM & fem_boundary = model->getFEMBoundary(); fem_boundary.initShapeFunctions(); fem_boundary.computeNormalsOnQuadPoints(); - model->computeForcesFromFunction(trac,0); + model->computeForcesFromFunction(trac, akantu::_bft_forces); #ifdef AKANTU_USE_IOHELPER DumperParaview dumper; dumper.SetMode(BASE64); dumper.SetPoints(model->getFEM().getMesh().getNodes().values, 2, nb_nodes, "coordinates"); dumper.SetConnectivity((int *)model->getFEM().getMesh().getConnectivity(_triangle_6).values, TRIANGLE2, model->getFEM().getMesh().getNbElement(_triangle_6), C_MODE); dumper.AddNodeDataField(model->getDisplacement().values, 2, "displacements"); dumper.AddNodeDataField(model->getVelocity().values, 2, "velocity"); dumper.AddNodeDataField(model->getForce().values, 2, "force"); dumper.AddNodeDataField(model->getMass().values, 1, "Mass"); dumper.AddNodeDataField(model->getResidual().values, 2, "residual"); dumper.AddElemDataField(model->getMaterial(0).getStrain(_triangle_6).values, 4, "strain"); dumper.AddElemDataField(model->getMaterial(0).getStress(_triangle_6).values, 4, "stress"); dumper.SetEmbeddedValue("displacements", 1); dumper.SetEmbeddedValue("force", 1); dumper.SetEmbeddedValue("residual", 1); dumper.SetEmbeddedValue("velocity", 1); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); #endif //AKANTU_USE_IOHELPER for(UInt s = 0; s < max_steps; ++s) { model->explicitPred(); model->updateResidual(); model->updateAcceleration(); model->explicitCorr(); epot = model->getPotentialEnergy(); ekin = model->getKineticEnergy(); std::cout << s << " " << epot << " " << ekin << " " << epot + ekin << std::endl; #ifdef AKANTU_USE_IOHELPER if(s % 100 == 0) dumper.Dump(); #endif //AKANTU_USE_IOHELPER } return EXIT_SUCCESS; }