diff --git a/common/aka_common.hh b/common/aka_common.hh index c0b65bb63..dda80f0c8 100644 --- a/common/aka_common.hh +++ b/common/aka_common.hh @@ -1,351 +1,353 @@ /** * @file aka_common.hh * @author Nicolas Richart * @date Fri Jun 11 09:48:06 2010 * * @namespace akantu * * @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 * * 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 unsigned long long UInt64; typedef signed int Int; typedef std::string ID; static const Real UINT_INIT_VALUE = 0; #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 ID SolverID; typedef ID ShapeID; typedef ID IntegratorID; typedef UInt Surface; /* -------------------------------------------------------------------------- */ // BOOST PART: TOUCH ONLY IF YOU KNOW WHAT YOU ARE DOING #include #define AKANTU_BOOST_CASE_MACRO(r,macro,type) \ case type : { macro(type); break;} #define AKANTU_BOOST_ELEMENT_SWITCH(macro); \ do { \ switch(type) { \ BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_CASE_MACRO,macro,AKANTU_ELEMENT_TYPE) \ case _not_defined: \ case _max_element_type: { \ AKANTU_DEBUG_ERROR("Wrong type : " << type); \ break; \ } \ } \ } while(0) #define AKANTU_BOOST_LIST_MACRO(r,macro,type) \ macro(type) #define AKANTU_BOOST_ELEMENT_LIST(macro) \ BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_LIST_MACRO,macro,AKANTU_ELEMENT_TYPE) /* -------------------------------------------------------------------------- */ /// @boost sequence of element to loop on in global tasks #define AKANTU_ELEMENT_TYPE \ (_segment_2) \ (_segment_3) \ (_triangle_3) \ (_triangle_6) \ (_tetrahedron_4) \ (_tetrahedron_10) \ (_quadrangle_4) \ (_quadrangle_8) \ (_hexahedron_8) \ (_point) /// @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 _quadrangle_4, /// first order quadrangle _quadrangle_8, /// second order quadrangle _hexahedron_8, /// first order hexahedron _point, /// point only for some algorithm to be generic like mesh partitioning _max_element_type }; /// @enum MaterialType different materials implemented enum MaterialType { _elastic = 0, _max_material_type }; typedef void (*BoundaryFunction)(double *,double *); /// @enum BoundaryFunctionType type of function passed for boundary conditions enum BoundaryFunctionType { _bft_stress, _bft_forces }; /// @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, _ct_rigid = 3 }; enum ContactSearchType { _cst_not_defined = 0, _cst_2d_expli = 1, _cst_expli = 2 }; enum ContactNeighborStructureType { _cnst_not_defined = 0, _cnst_regular_grid = 1, _cnst_2d_grid = 2 }; /* -------------------------------------------------------------------------- */ /* Ghosts handling */ /* -------------------------------------------------------------------------- */ typedef ID SynchronizerID; /// @enum CommunicatorType type of communication method to use enum CommunicatorType { _communicator_mpi, _communicator_dummy }; /// @enum SynchronizationTag type of synchronizations enum SynchronizationTag { /// SolidMechanicsModel tags _gst_smm_mass, /// synchronization of the SolidMechanicsModel.mass _gst_smm_for_strain, /// synchronization of the SolidMechanicsModel.current_position _gst_smm_boundary, /// synchronization of the boundary, forces, velocities and displacement _gst_smm_uv, /// synchronization of the nodal velocities and displacement + _gst_htm_capacity, /// synchronization of the nodal heat capacity + _gst_htm_temperature, /// synchronization of the nodal temperature /// 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 _quadrangle_4 : stream << "quadrangle_4" ; break; case _quadrangle_8 : stream << "quadrangle_8" ; break; case _hexahedron_8 : stream << "hexahedron_8" ; 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__ //#include "aka_types.hh" #endif /* __AKANTU_COMMON_HH__ */ diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 1d0161e80..1c59eb158 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,46 +1,47 @@ #=============================================================================== # @file CMakeLists.txt # @author Nicolas Richart # @author Guillaume Anciaux # @date Fri Jun 11 09:46:59 2010 # # @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 # #=============================================================================== #=============================================================================== # List of examples #=============================================================================== add_example(2d_compression "Compression example") add_example(regular_meshes "Generate a regular mesh") +add_example(heat_propagation_explicit "Heat propagation through a bottleneck") if(AKANTU_PARTITIONER_ON AND AKANTU_MPI_ON) add_example(3d_cube_compression "Cube compression example") add_example(scalability_test "Scalability test") endif() if(AKANTU_EPSN_ON) add_example(epsn "Example of steering akantu with EPSN") endif(AKANTU_EPSN_ON) if(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/local_examples.cmake") include("local_examples.cmake") endif(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/local_examples.cmake") diff --git a/examples/heat_propagation_explicit/heat_propa_pbc.cc b/examples/heat_propagation_explicit/heat_propa_pbc.cc index 35ca2ebef..23ffb7103 100644 --- a/examples/heat_propagation_explicit/heat_propa_pbc.cc +++ b/examples/heat_propagation_explicit/heat_propa_pbc.cc @@ -1,199 +1,165 @@ /** * @file test_heat_transfer_model_cube3d.cc * @author Rui WANG * @date Tue May 17 11:31:22 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; /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER #include "io_helper.h" #endif //AKANTU_USE_IOHELPER void paraviewInit(akantu::HeatTransferModel * model,Dumper & dumper); void paraviewDump(Dumper & dumper); akantu::UInt spatial_dimension = 3; akantu:: ElementType type = akantu::_tetrahedron_4; akantu::UInt paraview_type = TETRA1; /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { akantu::initialize(&argc,&argv); akantu::Mesh mesh(spatial_dimension); akantu::MeshIOMSH mesh_io; mesh_io.read("double_cube_tet4.msh", mesh); akantu::HeatTransferModel * model; akantu::UInt nb_nodes; akantu::UInt nb_element; - mesh.computeBoundingBox(); - std::map pbc_pair; - akantu::MeshUtils::computePBCMap(mesh,0,pbc_pair); - akantu::MeshUtils::computePBCMap(mesh,1,pbc_pair); - - { - std::map::iterator it = pbc_pair.begin(); - std::map::iterator end = pbc_pair.end(); - - akantu::Real * coords = mesh.getNodes().values; - akantu::UInt dim = mesh.getSpatialDimension(); - while(it != end){ - akantu::UInt i1 = (*it).first; - akantu::UInt i2 = (*it).second; - - AKANTU_DEBUG_INFO("pairing " << i1 << "(" - << coords[dim*i1] << "," << coords[dim*i1+1] << "," - << coords[dim*i1+2] - << ") with" - << i2 << "(" - << coords[dim*i2] << "," << coords[dim*i2+1] << "," - << coords[dim*i2+2] - << ")"); - ++it; - } - } - - akantu::PBCSynchronizer synch(pbc_pair); model = new akantu::HeatTransferModel(mesh); - model->createSynchronizerRegistry(model); - model->getSynchronizerRegistry().registerSynchronizer(synch,akantu::_gst_htm_capacity); - model->getSynchronizerRegistry().registerSynchronizer(synch,akantu::_gst_htm_temperature); - + model->initPBC(1,1,0); /* -------------------------------------------------------------------------- */ model->readMaterials("material.dat"); model->initModel(); model->initVectors(); model->getHeatFlux().clear(); model->getCapacityLumped().clear(); model->getTemperatureGradient(type).clear(); /* -------------------------------------------------------------------------- */ - model->changeLocalEquationNumberforPBC(pbc_pair,1); model->assembleCapacityLumped(type); - model->getSynchronizerRegistry().synchronize(akantu::_gst_htm_capacity); /* -------------------------------------------------------------------------- */ nb_nodes = model->getFEM().getMesh().getNbNodes(); nb_element = model->getFEM().getMesh().getNbElement(type); nb_nodes = model->getFEM().getMesh().getNbNodes(); /* ------------------------------------------------------------------------ */ //get stable time step akantu::Real time_step = model->getStableTimeStep()*0.8; cout<<"time step is:"<setTimeStep(time_step); /* -------------------------------------------------------------------------- */ /// boundary conditions const akantu::Vector & nodes = model->getFEM().getMesh().getNodes(); akantu::Vector & boundary = model->getBoundary(); akantu::Vector & temperature = model->getTemperature(); akantu::Vector & heat_flux = model->getHeatFlux(); akantu::Real eps = 1e-15; double t1, t2, length; t1 = 300.; t2 = 100.; length = 1.; for (akantu::UInt i = 0; i < nb_nodes; ++i) { temperature(i) = 100.; akantu::Real dz = nodes(i,2) - mesh.getZMin(); if(fabs(dz) < 0.1){ boundary(i) = true; temperature(i) = 150.; } } /* -------------------------------------------------------------------------- */ DumperParaview dumper; paraviewInit(model,dumper); /* ------------------------------------------------------------------------ */ // //for testing int max_steps = 1000000; /* ------------------------------------------------------------------------ */ for(int i=0; iupdateHeatFlux(); model->updateTemperature(); - model->getSynchronizerRegistry().synchronize(akantu::_gst_htm_temperature); if(i % 1000 == 0){ paraviewDump(dumper); std::cout << "Step " << i << "/" << max_steps << std::endl; } } cout<< "\n\n Stable Time Step is : " << time_step << "\n \n" <getFEM().getMesh().getNbNodes(); akantu::UInt nb_element = model->getFEM().getMesh().getNbElement(type); dumper.SetMode(TEXT); dumper.SetPoints(model->getFEM().getMesh().getNodes().values, spatial_dimension, nb_nodes, "coordinates2"); dumper.SetConnectivity((int *)model->getFEM().getMesh().getConnectivity(type).values, paraview_type, nb_element, C_MODE); dumper.AddNodeDataField(model->getTemperature().values, 1, "temperature"); dumper.AddNodeDataField(model->getHeatFlux().values, 1, "heat_flux"); dumper.AddNodeDataField(model->getCapacityLumped().values, 1, "capacity_lumped"); // dumper.AddElemDataField(model->getTemperatureGradient(type).values, // spatial_dimension, "temperature_gradient"); // dumper.AddElemDataField(model->getbtkgt().values, // 4, "btkgt"); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); } /* -------------------------------------------------------------------------- */ void paraviewDump(Dumper & dumper) { dumper.Dump(); } /* -------------------------------------------------------------------------- */ diff --git a/fem/fem.cc b/fem/fem.cc index 9b77fc989..996159ba6 100644 --- a/fem/fem.cc +++ b/fem/fem.cc @@ -1,312 +1,245 @@ /** * @file fem.cc * @author Nicolas Richart * @author Guillaume Anciaux * @date Fri Jul 16 11:03:02 2010 * * @brief Implementation of the FEM class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "fem.hh" #include "mesh.hh" #include "element_class.hh" #include "static_communicator.hh" #include "aka_math.hh" #include "dof_synchronizer.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ FEM::FEM(Mesh & mesh, UInt element_dimension, FEMID id, MemoryID memory_id) : Memory(memory_id), id(id) { AKANTU_DEBUG_IN(); this->element_dimension = (element_dimension != 0) ? element_dimension : mesh.getSpatialDimension(); init(); this->mesh = &mesh; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FEM::init() { for(UInt t = _not_defined; t < _max_element_type; ++t) { this->normals_on_quad_points[t] = NULL; } } /* -------------------------------------------------------------------------- */ FEM::~FEM() { AKANTU_DEBUG_IN(); mesh = NULL; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FEM::assembleVector(const Vector & elementary_vect, Vector & nodal_values, const Vector & equation_number, UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type, const Vector * filter_elements, Real scale_factor) const { AKANTU_DEBUG_IN(); UInt nb_element; UInt * conn_val; if(ghost_type == _not_ghost) { nb_element = mesh->getNbElement(type); conn_val = mesh->getConnectivity(type).values; } else { nb_element = mesh->getNbGhostElement(type); conn_val = mesh->getGhostConnectivity(type).values; } UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_nodes = mesh->getNbNodes(); UInt * filter_elem_val = NULL; if(filter_elements != NULL) { nb_element = filter_elements->getSize(); filter_elem_val = filter_elements->values; } AKANTU_DEBUG_ASSERT(elementary_vect.getSize() == nb_element, "The vector elementary_vect(" << elementary_vect.getID() << ") has not the good size."); AKANTU_DEBUG_ASSERT(elementary_vect.getNbComponent() == nb_degre_of_freedom*nb_nodes_per_element, "The vector elementary_vect(" << elementary_vect.getID() << ") has not the good number of component." << "(" << elementary_vect.getNbComponent() << " != " << nb_degre_of_freedom*nb_nodes_per_element << ")"); AKANTU_DEBUG_ASSERT(nodal_values.getNbComponent() == nb_degre_of_freedom, "The vector nodal_values(" << nodal_values.getID() << ") has not the good number of component." << "(" << nodal_values.getNbComponent() << " != " << nb_degre_of_freedom << ")"); nodal_values.resize(nb_nodes); Real * elementary_vect_val = elementary_vect.values; Real * nodal_values_val = nodal_values.values; for (UInt el = 0; el < nb_element; ++el) { UInt el_offset = el * nb_nodes_per_element; if(filter_elements != NULL) { el_offset = filter_elem_val[el] * nb_nodes_per_element; } for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = conn_val[el_offset + n]; UInt offset_node = node * nb_degre_of_freedom; for (UInt d = 0; d < nb_degre_of_freedom; ++d) { nodal_values_val[equation_number.values[offset_node + d]] += scale_factor * elementary_vect_val[d]; } elementary_vect_val += nb_degre_of_freedom; } } AKANTU_DEBUG_OUT(); } -/* -------------------------------------------------------------------------- */ -void FEM::assembleVector(const Vector & elementary_vect, - Vector & nodal_values, - UInt nb_degre_of_freedom, - const ElementType & type, - GhostType ghost_type, - const Vector * filter_elements, - Real scale_factor) const { - AKANTU_DEBUG_IN(); - - UInt nb_element; - UInt * conn_val; - - if(ghost_type == _not_ghost) { - nb_element = mesh->getNbElement(type); - conn_val = mesh->getConnectivity(type).values; - } else { - nb_element = mesh->getNbGhostElement(type); - conn_val = mesh->getGhostConnectivity(type).values; - } - - UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); - UInt nb_nodes = mesh->getNbNodes(); - - UInt * filter_elem_val = NULL; - if(filter_elements != NULL) { - nb_element = filter_elements->getSize(); - filter_elem_val = filter_elements->values; - } - - AKANTU_DEBUG_ASSERT(elementary_vect.getSize() == nb_element, - "The vector elementary_vect(" << elementary_vect.getID() - << ") has not the good size."); - - AKANTU_DEBUG_ASSERT(elementary_vect.getNbComponent() - == nb_degre_of_freedom*nb_nodes_per_element, - "The vector elementary_vect(" << elementary_vect.getID() - << ") has not the good number of component."); - - AKANTU_DEBUG_ASSERT(nodal_values.getNbComponent() == nb_degre_of_freedom, - "The vector nodal_values(" << nodal_values.getID() - << ") has not the good number of component."); - - nodal_values.resize(nb_nodes); - - Real * elementary_vect_val = elementary_vect.values; - Real * nodal_values_val = nodal_values.values; - - for (UInt el = 0; el < nb_element; ++el) { - UInt el_offset = el * nb_nodes_per_element; - if(filter_elements != NULL) { - el_offset = filter_elem_val[el] * nb_nodes_per_element; - } - for (UInt n = 0; n < nb_nodes_per_element; ++n) { - UInt node = conn_val[el_offset + n]; - UInt offset_node = node * nb_degre_of_freedom; - for (UInt d = 0; d < nb_degre_of_freedom; ++d) { - nodal_values_val[offset_node + d] += scale_factor * elementary_vect_val[d]; - } - elementary_vect_val += nb_degre_of_freedom; - } - } - - AKANTU_DEBUG_OUT(); -} - - /* -------------------------------------------------------------------------- */ void FEM::assembleMatrix(const Vector & elementary_mat, SparseMatrix & matrix, UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type, const Vector * filter_elements) const { AKANTU_DEBUG_IN(); UInt nb_element; if(ghost_type == _not_ghost) { nb_element = mesh->getNbElement(type); } else { AKANTU_DEBUG_TO_IMPLEMENT(); } UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt * filter_elem_val = NULL; if(filter_elements != NULL) { nb_element = filter_elements->getSize(); filter_elem_val = filter_elements->values; } AKANTU_DEBUG_ASSERT(elementary_mat.getSize() == nb_element, "The vector elementary_mat(" << elementary_mat.getID() << ") has not the good size."); AKANTU_DEBUG_ASSERT(elementary_mat.getNbComponent() == nb_degre_of_freedom * nb_nodes_per_element * nb_degre_of_freedom * nb_nodes_per_element, "The vector elementary_mat(" << elementary_mat.getID() << ") has not the good number of component."); Real * elementary_mat_val = elementary_mat.values; UInt offset_elementary_mat = elementary_mat.getNbComponent(); UInt * connectivity_val = mesh->getConnectivity(type).values; UInt size_mat = nb_nodes_per_element * nb_degre_of_freedom; UInt size = mesh->getNbGlobalNodes() * nb_degre_of_freedom; Int * eq_nb_val = matrix.getDOFSynchronizer().getGlobalDOFEquationNumbers().values; Int * local_eq_nb_val = new Int[nb_degre_of_freedom * nb_nodes_per_element]; for (UInt e = 0; e < nb_element; ++e) { UInt el = e; if(filter_elements != NULL) el = filter_elem_val[e]; Int * tmp_local_eq_nb_val = local_eq_nb_val; UInt * conn_val = connectivity_val + el * nb_nodes_per_element; for (UInt i = 0; i < nb_nodes_per_element; ++i) { UInt n = conn_val[i]; for (UInt d = 0; d < nb_degre_of_freedom; ++d) { *tmp_local_eq_nb_val++ = eq_nb_val[n * nb_degre_of_freedom + d]; } // memcpy(tmp_local_eq_nb_val, eq_nb_val + n * nb_degre_of_freedom, nb_degre_of_freedom * sizeof(Int)); // tmp_local_eq_nb_val += nb_degre_of_freedom; } for (UInt i = 0; i < size_mat; ++i) { UInt c_irn = local_eq_nb_val[i]; if(c_irn < size) { UInt j_start = (matrix.getSparseMatrixType() == _symmetric) ? i : 0; for (UInt j = j_start; j < size_mat; ++j) { UInt c_jcn = local_eq_nb_val[j]; if(c_jcn < size) { matrix(c_irn, c_jcn) += elementary_mat_val[i * size_mat + j]; } } } } elementary_mat_val += offset_elementary_mat; } delete [] local_eq_nb_val; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FEM::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "FEM [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + element dimension : " << element_dimension << std::endl; stream << space << " + mesh [" << std::endl; mesh->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + mesh [" << std::endl; mesh->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } __END_AKANTU__ diff --git a/fem/fem.hh b/fem/fem.hh index afce28ec4..c73bf4284 100644 --- a/fem/fem.hh +++ b/fem/fem.hh @@ -1,244 +1,235 @@ /** * @file fem.hh * @author Nicolas Richart * @date Fri Jul 16 10:24:24 2010 * * @brief FEM 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_FEM_HH__ #define __AKANTU_FEM_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "mesh.hh" #include "element_class.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class FEM : public Memory { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: FEM(Mesh & mesh, UInt spatial_dimension = 0, FEMID id = "fem", MemoryID memory_id = 0); virtual ~FEM(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// build the profile of the sparse matrix corresponding to the mesh void initSparseMatrixProfile(SparseMatrixType sparse_matrix_type = _unsymmetric); /// pre-compute all the shape functions, their derivatives and the jacobians virtual void initShapeFunctions(GhostType ghost_type = _not_ghost)=0; /* ------------------------------------------------------------------------ */ /* Integration method bridges */ /* ------------------------------------------------------------------------ */ /// integrate f for all elements of type "type" virtual void integrate(const Vector & f, Vector &intf, UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL) const = 0; /// integrate a scalar value on all elements of type "type" virtual Real integrate(const Vector & f, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL) const = 0; /* ------------------------------------------------------------------------ */ /* compatibility with old FEM fashion */ /* ------------------------------------------------------------------------ */ /// get the number of quadrature points virtual UInt getNbQuadraturePoints(const ElementType & type)=0; /// get the precomputed shapes const virtual Vector & getShapes(const ElementType & type)=0; /// get the precomputed ghost shapes const virtual Vector & getGhostShapes(const ElementType & type)=0; /// get the derivatives of shapes const virtual Vector & getShapesDerivatives(const ElementType & type)=0; /// get the ghost derivatives of shapes const virtual Vector & getGhostShapesDerivatives(const ElementType & type)=0; /// get quadrature points const virtual Vector & getQuadraturePoints(const ElementType & type)=0; /* ------------------------------------------------------------------------ */ /* Shape method bridges */ /* ------------------------------------------------------------------------ */ virtual void gradientOnQuadraturePoints(const Vector &u, Vector &nablauq, const UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL)=0; virtual void interpolateOnQuadraturePoints(const Vector &u, Vector &uq, UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL) const =0; /* ------------------------------------------------------------------------ */ /* Other methods */ /* ------------------------------------------------------------------------ */ /// pre-compute normals on control points virtual void computeNormalsOnControlPoints(GhostType ghost_type = _not_ghost)=0; /// assemble vectors void assembleVector(const Vector & elementary_vect, Vector & nodal_values, const Vector & equation_number, UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL, Real scale_factor = 1) const; - /// assemble vectors - void assembleVector(const Vector & elementary_vect, - Vector & nodal_values, - UInt nb_degre_of_freedom, - const ElementType & type, - GhostType ghost_type = _not_ghost, - const Vector * filter_elements = NULL, - Real scale_factor = 1) const; - /// assemble matrix in the complete sparse matrix void assembleMatrix(const Vector & elementary_mat, SparseMatrix & matrix, UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL) const; /// assemble a field as a lumped matrix (ex. rho in lumped mass) virtual void assembleFieldLumped(__attribute__ ((unused)) const Vector & field_1, __attribute__ ((unused)) UInt nb_degree_of_freedom, __attribute__ ((unused)) Vector & lumped, __attribute__ ((unused)) const Vector & equation_number, __attribute__ ((unused)) ElementType type, __attribute__ ((unused)) __attribute__ ((unused)) GhostType ghost_type) { AKANTU_DEBUG_TO_IMPLEMENT(); }; /// assemble a field as a matrix (ex. rho to mass matrix) virtual void assembleFieldMatrix(__attribute__ ((unused)) const Vector & field_1, __attribute__ ((unused)) UInt nb_degree_of_freedom, __attribute__ ((unused)) const Vector & equation_number, __attribute__ ((unused)) SparseMatrix & matrix, __attribute__ ((unused)) ElementType type, __attribute__ ((unused)) GhostType ghost_type) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; private: /// initialise the class void init(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(ElementDimension, element_dimension, UInt); /// get the mesh contained in the fem object inline Mesh & getMesh() const; /// get the in-radius of an element static inline Real getElementInradius(Real * coord, const ElementType & type); /// get the normals on quadrature points AKANTU_GET_MACRO_BY_ELEMENT_TYPE(NormalsOnQuadPoints, normals_on_quad_points, const Vector &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id of the fem object FEMID id; /// spatial dimension of the problem UInt element_dimension; /// the mesh on which all computation are made Mesh * mesh; /// normals at quadrature points ByElementTypeReal normals_on_quad_points; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "fem_inline_impl.cc" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const FEM & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #include "fem_template.hh" #endif /* __AKANTU_FEM_HH__ */ diff --git a/mesh_utils/mesh_pbc.cc b/mesh_utils/mesh_pbc.cc index ff4a39f0d..bd260b92e 100644 --- a/mesh_utils/mesh_pbc.cc +++ b/mesh_utils/mesh_pbc.cc @@ -1,239 +1,254 @@ /** * @file mesh_pbc.cc * @author Guillaume Anciaux * @date Tue Feb 8 10:48:01 2011 * * @brief periodic boundary condition connectivity tweak * * @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 "mesh_utils.hh" #include __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /// class that sorts a set of nodes of same coordinates in 'dir' direction class CoordinatesComparison { public: CoordinatesComparison (const UInt dimension, const UInt dirx, const UInt diry, Real * coords): dim(dimension),dir_x(dirx),dir_y(diry),coordinates(coords){} bool operator() (UInt n1, UInt n2){ Real p1_x = coordinates[dim*n1+dir_x]; Real p2_x = coordinates[dim*n2+dir_x]; - Real p1_y = coordinates[dim*n1+dir_y]; - Real p2_y = coordinates[dim*n2+dir_y]; Real diff_x = p1_x - p2_x; - if (fabs(diff_x) > Math::tolerance) + if (dim == 2 || fabs(diff_x) > Math::tolerance) return diff_x > 0.0 ? false : true; - else { + else if (dim > 2){ + Real p1_y = coordinates[dim*n1+dir_y]; + Real p2_y = coordinates[dim*n2+dir_y]; Real diff_y = p1_y - p2_y; return diff_y > 0 ? false : true; } } private: UInt dim; UInt dir_x; UInt dir_y; Real * coordinates; }; /* -------------------------------------------------------------------------- */ // void MeshUtils::tweakConnectivityForPBC(Mesh & mesh, // bool flag_x, // bool flag_y, // bool flag_z){ // std::map pbc_pair; // mesh.computeBoundingBox(); // mesh.pbc_directions[0] = flag_x; // mesh.pbc_directions[1] = flag_y; // mesh.pbc_directions[2] = flag_z; // if (flag_x) computePBCMap(mesh,0,pbc_pair); // if (flag_y) computePBCMap(mesh,1,pbc_pair); // if (flag_z) computePBCMap(mesh,2,pbc_pair); // { // std::map::iterator it = pbc_pair.begin(); // std::map::iterator end = pbc_pair.end(); // Real * coords = mesh.nodes->values; // UInt dim = mesh.getSpatialDimension(); // while(it != end){ // UInt i1 = (*it).first; // UInt i2 = (*it).second; // AKANTU_DEBUG_INFO("pairing " << i1 << "(" // << coords[dim*i1] << "," << coords[dim*i1+1] << "," // << coords[dim*i1+2] // << ") with" // << i2 << "(" // << coords[dim*i2] << "," << coords[dim*i2+1] << "," // << coords[dim*i2+2] // << ")"); // ++it; // } // } // //allocate and initialize list of reversed elements // mesh.initByElementTypeUIntVector(mesh.reversed_elements_pbc,1,0,mesh.id,"reversed"); // // now loop over the elements to change the connectivity of some elements // const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); // Mesh::ConnectivityTypeList::const_iterator it; // for(it = type_list.begin(); it != type_list.end(); ++it) { // ElementType type = *it; // UInt nb_elem = mesh.getNbElement(type); // UInt nb_nodes_per_elem = mesh.getNbNodesPerElement(type); // UInt * conn = mesh.getConnectivityPointer(type)->values; // UInt index = 0; // Vector & list = *(mesh.reversed_elements_pbc[type]); // for (UInt el = 0; el < nb_elem; el++) { // UInt flag_should_register_elem = false; // for (UInt k = 0; k < nb_nodes_per_elem; ++k,++index){ // if (pbc_pair.count(conn[index])){ // flag_should_register_elem = true; // AKANTU_DEBUG_INFO("elem list size " << list.getSize()); // AKANTU_DEBUG_INFO("node " << conn[index] +1 // << " switch to " // << pbc_pair[conn[index]]+1); // // for (UInt toto = 0; toto < 3; ++toto) { // // AKANTU_DEBUG_INFO("dir " << toto << " coords " // // << mesh.nodes->values[conn[index]*3+toto] // // << " switch to " // // << mesh.nodes->values[pbc_pair[conn[index]]*3+toto]); // // } // std::stringstream str_temp; // str_temp << "initial elem(" << el << ") is "; // for (UInt l = 0 ; l < nb_nodes_per_elem ; ++ l){ // str_temp << conn[el*nb_nodes_per_elem+l]+1 << " "; // } // AKANTU_DEBUG_INFO(str_temp.str()); // conn[index] = pbc_pair[conn[index]]; // } // } // if (flag_should_register_elem) list.push_back(el); // } // } // } /* -------------------------------------------------------------------------- */ void MeshUtils::computePBCMap(const Mesh & mymesh,const UInt dir, std::map & pbc_pair){ std::vector selected_left; std::vector selected_right; Real * coords = mymesh.nodes->values; const UInt nb_nodes = mymesh.nodes->getSize(); const UInt dim = mymesh.getSpatialDimension(); + + if (dim <= dir) return; + AKANTU_DEBUG_INFO("min " << mymesh.xmin[dir]); AKANTU_DEBUG_INFO("max " << mymesh.xmax[dir]); for (UInt i = 0; i < nb_nodes; ++i) { AKANTU_DEBUG_TRACE("treating " << coords[dim*i+dir]); if(Math::are_float_equal(coords[dim*i+dir],mymesh.xmin[dir])){ AKANTU_DEBUG_TRACE("pushing node " << i << " on the left side"); selected_left.push_back(i); } else if(Math::are_float_equal(coords[dim*i+dir],mymesh.xmax[dir])){ selected_right.push_back(i); AKANTU_DEBUG_TRACE("pushing node " << i << " on the right side"); } } AKANTU_DEBUG_INFO("found " << selected_left.size() << " and " << selected_right.size() << " nodes at each boundary for direction " << dir); UInt dir_x,dir_y; - if (dir == 0){ - dir_x = 1;dir_y = 2; - } - else if (dir == 1){ - dir_x = 0;dir_y = 2; + if (dim == 3){ + if (dir == 0){ + dir_x = 1;dir_y = 2; + } + else if (dir == 1){ + dir_x = 0;dir_y = 2; + } + else if (dir == 2){ + dir_x = 0;dir_y = 1; + } } - else if (dir == 2){ - dir_x = 0;dir_y = 1; + else if (dim == 2){ + if (dir == 0){ + dir_x = 1; + } + else if (dir == 1){ + dir_x = 0; + } } - + CoordinatesComparison compare_nodes(dim,dir_x,dir_y,coords); std::sort(selected_left.begin(),selected_left.end(),compare_nodes); std::sort(selected_right.begin(),selected_right.end(),compare_nodes); std::vector::iterator it_left = selected_left.begin(); std::vector::iterator end_left = selected_left.end(); std::vector::iterator it_right = selected_right.begin(); std::vector::iterator end_right = selected_right.end(); while ((it_left != end_left) && (it_right != end_right) ){ UInt i1 = *it_left; UInt i2 = *it_right; AKANTU_DEBUG_TRACE("do I pair? " << i1 << "(" << coords[dim*i1] << "," << coords[dim*i1+1] << "," << coords[dim*i1+2] << ") with" << i2 << "(" << coords[dim*i2] << "," << coords[dim*i2+1] << "," << coords[dim*i2+2] << ") in direction " << dir); - Real dx = coords[dim*i1+dir_x] - coords[dim*i2+dir_x]; - Real dy = coords[dim*i1+dir_y] - coords[dim*i2+dir_y]; - + Real dx = 0.0; + Real dy = 0.0; + if (dim == 2) dx = coords[dim*i1+dir_x] - coords[dim*i2+dir_x]; + if (dim == 3) dy = coords[dim*i1+dir_y] - coords[dim*i2+dir_y]; + if (fabs(dx*dx+dy*dy) < Math::tolerance) { //then i match these pairs if (pbc_pair.count(i2)){ i2 = pbc_pair[i2]; } pbc_pair[i1] = i2; AKANTU_DEBUG_TRACE("pairing " << i1 << "(" << coords[dim*i1] << "," << coords[dim*i1+1] << "," << coords[dim*i1+2] << ") with" << i2 << "(" << coords[dim*i2] << "," << coords[dim*i2+1] << "," << coords[dim*i2+2] << ") in direction " << dir); ++it_left; ++it_right; } else if (fabs(dy) < Math::tolerance && dx > 0) ++it_right; else if (fabs(dy) < Math::tolerance && dx < 0) ++it_left; else if (dy > 0) ++it_right; else if (dy < 0) ++it_left; else { AKANTU_DEBUG_ERROR("this should not append"); } } AKANTU_DEBUG_INFO("found " << pbc_pair.size() << " pairs for direction " << dir); } __END_AKANTU__ diff --git a/model/heat_transfer/heat_transfer_model.cc b/model/heat_transfer/heat_transfer_model.cc index dd27cc533..3deb951bc 100644 --- a/model/heat_transfer/heat_transfer_model.cc +++ b/model/heat_transfer/heat_transfer_model.cc @@ -1,420 +1,433 @@ /** * @file heat_transfer_model.cc * @author Rui WANG * @date Fri May 4 13:46:43 2011 * * @brief Implementation of 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 . *el_size */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "heat_transfer_model.hh" #include "aka_math.hh" //#include "integration_scheme_2nd_order.hh" #include "aka_common.hh" #include #include #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "static_communicator.hh" #include "sparse_matrix.hh" #include "solver.hh" //#include "solid_mechanics_model.hh" #include"sparse_matrix.hh" // #include "fem_template.cc" #include "fem_template.hh" using namespace std; #ifdef AKANTU_USE_MUMPS #include "solver_mumps.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ HeatTransferModel::HeatTransferModel(Mesh & mesh, UInt dim, const ModelID & id, const MemoryID & memory_id): Model(id, memory_id), spatial_dimension(dim) { AKANTU_DEBUG_IN(); + createSynchronizerRegistry(this); + if (spatial_dimension == 0) spatial_dimension = mesh.getSpatialDimension(); registerFEMObject("HeatTransferFEM",mesh,spatial_dimension); this->temperature= NULL; for (UInt t = _not_defined; t < _max_element_type; ++t) { this->temperature_gradient[t] = NULL; } this->heat_flux = NULL; this->boundary = NULL; + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void HeatTransferModel::initModel() { getFEM().initShapeFunctions(_not_ghost); getFEM().initShapeFunctions(_ghost); } - +/* -------------------------------------------------------------------------- */ +void HeatTransferModel::initPBC(UInt x, UInt y, UInt z){ + Model::initPBC(x,y,z); + PBCSynchronizer * synch = new PBCSynchronizer(pbc_pair); + synch_registry->registerSynchronizer(*synch,akantu::_gst_htm_capacity); + synch_registry->registerSynchronizer(*synch,akantu::_gst_htm_temperature); + changeLocalEquationNumberforPBC(pbc_pair,1); +} +/* -------------------------------------------------------------------------- */ void HeatTransferModel::initVectors() { AKANTU_DEBUG_IN(); UInt nb_nodes = getFEM().getMesh().getNbNodes(); std::stringstream sstr_temp; sstr_temp << id << ":temperature"; std::stringstream sstr_heat_flux; sstr_heat_flux << id << ":heat flux"; std::stringstream sstr_lump; sstr_lump<(sstr_temp.str(), nb_nodes, 1, REAL_INIT_VALUE)); heat_flux= &(alloc(sstr_heat_flux.str(), nb_nodes, 1, REAL_INIT_VALUE)); capacity_lumped= &(alloc(sstr_lump.str(), nb_nodes, 1, REAL_INIT_VALUE)); boundary= &(alloc(sstr_boun.str(), nb_nodes, 1, false)); const Mesh::ConnectivityTypeList & type_list = this->getFEM().getMesh().getConnectivityTypeList(_not_ghost); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; std::stringstream sstr_temp_grad; sstr_temp_grad << id << ":temperature gradient:" << *it; UInt nb_element = getFEM().getMesh().getNbElement(*it); UInt nb_quad_points = this->getFEM().getNbQuadraturePoints(*it) * nb_element; std::cout << "Allocating " << sstr_temp_grad.str() << std::endl; temperature_gradient[*it] = &(alloc(sstr_temp_grad.str(), nb_quad_points, spatial_dimension, REAL_INIT_VALUE)); } dof_synchronizer = new DOFSynchronizer(getFEM().getMesh(),1); dof_synchronizer->initLocalDOFEquationNumbers(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ HeatTransferModel::~HeatTransferModel() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ //set boundary condition void setBoundaryCondition() { AKANTU_DEBUG_IN(); // Mesh::getcomputeBoundingBox(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleCapacityLumped(const ElementType &el_type) { AKANTU_DEBUG_IN(); FEM & fem = getFEM(); UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(el_type); UInt nb_element = getFEM().getMesh().getNbElement(el_type); //??? Vector rho_1 (nb_element * nb_quadrature_points,1, capacity * density); const Mesh::ConnectivityTypeList & type_list = fem.getMesh().getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; fem.assembleFieldLumped(rho_1,1,*capacity_lumped, dof_synchronizer->getLocalDOFEquationNumbers(), *it, _not_ghost); } + getSynchronizerRegistry().synchronize(akantu::_gst_htm_capacity); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ // compute the heat flux - void HeatTransferModel:: updateHeatFlux() + void HeatTransferModel::updateHeatFlux() { AKANTU_DEBUG_IN(); heat_flux->clear(); const Mesh::ConnectivityTypeList & type_list = this->getFEM().getMesh().getConnectivityTypeList(_not_ghost); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; const Vector * shapes_derivatives; shapes_derivatives = &(this->getFEM().getShapesDerivatives(*it)); UInt nb_element = getFEM().getMesh().getNbElement(*it); UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(*it); UInt size_of_shapes_derivatives = shapes_derivatives->getNbComponent(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); /*problem is included in computeTemperatureGradient*/ //const Vector *gT = computeTemperatureGradient(*it); computeTemperatureGradient(*it); Real * gT_val = temperature_gradient[*it]->values; Real * k_gT_val = new Real[spatial_dimension]; Real * shapes_derivatives_val = shapes_derivatives->values; UInt bt_k_gT_size = nb_nodes_per_element; - Vector *bt_k_gT = new Vector (nb_quadrature_points*nb_element, bt_k_gT_size); + Vector * bt_k_gT = new Vector (nb_quadrature_points*nb_element, bt_k_gT_size); Real * bt_k_gT_val = bt_k_gT->values; for (UInt el = 0; el < nb_element; ++el) { for (UInt i = 0; i < nb_quadrature_points; ++i) { Math::matrixt_vector(spatial_dimension, spatial_dimension, conductivity, gT_val, k_gT_val); gT_val += spatial_dimension; Math::matrix_vector(nb_nodes_per_element, spatial_dimension, shapes_derivatives_val, k_gT_val, bt_k_gT_val); shapes_derivatives_val += nb_nodes_per_element * spatial_dimension; bt_k_gT_val += bt_k_gT_size; } } Vector * q_e = new Vector(0,bt_k_gT_size,"q_e"); this->getFEM().integrate(*bt_k_gT, *q_e, bt_k_gT_size, *it); delete bt_k_gT; // SparseMatrix & K = new SparseMatrix(this->getFEM().getMesh().getNbNodes()); // const Vector & equation_number = this->getEquationNumber(); // this->getFEM().assembleMatrix(*q_e, K, equation_number, spatial_dimension, *it); + //the previous version of the assemble Vector, it is right, just ignore the external heat flux - this->getFEM().assembleVector(*q_e, *heat_flux, 1, *it); - //this->getFEM().assembleVector(*q_e, *heat_flux, 1, *it,_not_ghost,NULL); + this->getFEM().assembleVector(*q_e, *heat_flux, + dof_synchronizer->getLocalDOFEquationNumbers(), + 1, *it,_not_ghost,NULL,-1); delete q_e; AKANTU_DEBUG_OUT(); } } /* -------------------------------------------------------------------------- */ //compute temperature gradient void HeatTransferModel:: computeTemperatureGradient(const ElementType &el_type ) { AKANTU_DEBUG_IN(); this->getFEM().gradientOnQuadraturePoints(*temperature, *temperature_gradient[el_type], 1 ,el_type); AKANTU_DEBUG_OUT(); } /*except updateflux, everything is fine*/ //compute the temperature void HeatTransferModel::updateTemperature() { AKANTU_DEBUG_IN(); - - - // updateHeatFlux(); UInt nb_nodes=getFEM().getMesh().getNbNodes(); for(UInt i=0; i < nb_nodes; i++) { if(!((*boundary)(i))) { - //if((*temperature)(i)<0) (*temperature)(i)=0.0; - (*temperature)(i,0) = (*temperature)(i,0) - (*heat_flux)(i,0) * time_step / (*capacity_lumped)(i,0); - (*temperature)(i,0) = std::max((*temperature)(i,0), 5.0); - //for testing - //(*temperature)(i) = std::min((*temperature)(i), 300.0); + (*temperature)(i,0) += (*heat_flux)(i,0) * time_step / (*capacity_lumped)(i,0); + (*temperature)(i,0) = std::max((*temperature)(i,0), 0.0); } - - - //(*shapes_conductivity_shapes)(i)=(* shapes_conductivity_shapes)(i)*temperature(i); - //(*temperature)(i) -= (*shapes_conductivity_shapes)(i) * time_step / (*lumped)(i); - } - + + synch_registry->synchronize(akantu::_gst_htm_temperature); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ //choose scheme for temperature iteration /*a=0.0, the forward Euler; a=0.5, the Crank-Nicolson;a=2.0/3, Galerkin;a=1.0,backward Euler*/ void HeatTransferModel::integrationScheme1stOrder(Real thelta, UInt N, Vector * temperature) { AKANTU_DEBUG_IN(); //remember to use time_step if(thelta==0.0) { // temperature->values[i] += temperature->values[i] + heat_flux->values[i]/lumped->values[i]; } else if(thelta==0.5) { } else if(thelta==1.0) { } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Real conductivitymax = -std::numeric_limits::max(); Real el_size; Real min_el_size = std::numeric_limits::max(); Real * coord = getFEM().getMesh().getNodes().values; const Mesh::ConnectivityTypeList & type_list = getFEM().getMesh().getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(getFEM().getMesh().getSpatialDimension(*it) != spatial_dimension) continue; UInt nb_nodes_per_element = getFEM().getMesh().getNbNodesPerElement(*it); UInt nb_element = getFEM().getMesh().getNbElement(*it); UInt * conn = getFEM().getMesh().getConnectivity(*it).values; Real * u = new Real[nb_nodes_per_element*spatial_dimension]; for (UInt el = 0; el < nb_element; ++el) { UInt el_offset = el * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_conn = conn[el_offset + n] * spatial_dimension; memcpy(u + n * spatial_dimension, coord + offset_conn, spatial_dimension * sizeof(Real)); } el_size = getFEM().getElementInradius(u, *it); //get the biggest parameter from k11 until k33// for(UInt i = 0; i < spatial_dimension * spatial_dimension; i++) { if(conductivity[i] > conductivitymax) conductivitymax=conductivity[i]; } min_el_size = std::min(min_el_size, el_size); } cout<(*this); } else AKANTU_DEBUG_ERROR("did not find any section with material info for heat conduction"); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::setParam(const std::string & key, const std::string & value) { std::stringstream str(value); if (key == "conductivity"){ conductivity = new Real[spatial_dimension * spatial_dimension]; - for(int i=0;i> conductivity[i*spatial_dimension+j]; + if (i< spatial_dimension && j < spatial_dimension){ + str >> conductivity[i*spatial_dimension+j]; + AKANTU_DEBUG_INFO("conductivity(" << i << "," << j << ") = " + << conductivity[i*spatial_dimension+j]); + } + else { + Real tmp; + str >> tmp; + } } } else if (key == "capacity"){ str >> capacity; AKANTU_DEBUG_INFO("The capacity of the material is:" << capacity); } else if (key == "density"){ str >> density; AKANTU_DEBUG_INFO("The density of the material is:" << density); } } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/model/heat_transfer/heat_transfer_model.hh b/model/heat_transfer/heat_transfer_model.hh index a17d7c46a..f98fa2c60 100644 --- a/model/heat_transfer/heat_transfer_model.hh +++ b/model/heat_transfer/heat_transfer_model.hh @@ -1,260 +1,257 @@ /** * @file heat_transfer_model.hh * @author Rui WANG * @date Fri May 4 13:35:55 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 /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "model.hh" #include "material.hh" #include "parser.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #include "fem.hh" #include "mesh.hh" #include "aka_memory.hh" #include "element_class.hh" #include "sparse_matrix.hh" // namespace akantu { // class Solver; // class SparseMatrix; // } __BEGIN_AKANTU__ -class HeatTransferModel : public Model { +class HeatTransferModel : public Model, public DataAccessor { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef FEMTemplate MyFEMType; HeatTransferModel(UInt spatial_dimension, const ModelID & id = "heat_transfer_model", const MemoryID & memory_id = 0) ; HeatTransferModel(Mesh & mesh, UInt spatial_dimension = 0, const ModelID & id = "heat_transfer_model", const MemoryID & memory_id = 0); virtual ~HeatTransferModel() ; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// set the parameters void 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(); /// initialize the model void initModel(); + + /// init PBC synchronizer + void initPBC(UInt x,UInt y, UInt z); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const {}; /* ------------------------------------------------------------------------ */ - /* Ghost Synchronizer inherited members */ + /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: inline virtual UInt getNbDataToPack(const Element & element, SynchronizationTag tag) const { AKANTU_DEBUG_TO_IMPLEMENT(); }; inline virtual UInt getNbDataToUnpack(const Element & element, SynchronizationTag tag) const { AKANTU_DEBUG_TO_IMPLEMENT(); }; inline virtual void packData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag) const { AKANTU_DEBUG_TO_IMPLEMENT(); }; inline virtual void unpackData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag) const { AKANTU_DEBUG_TO_IMPLEMENT(); }; - - + inline virtual UInt getNbDataToPack(SynchronizationTag tag) const; + inline virtual UInt getNbDataToUnpack(SynchronizationTag tag) const; + inline virtual void packData(CommunicationBuffer & buffer, + const UInt index, + SynchronizationTag tag) const; + inline virtual void unpackData(CommunicationBuffer & buffer, + const UInt index, + SynchronizationTag tag) const; + /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); - - //* all the implementation of time step ----------------------------------- */ /// 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); - /// set the value of the time step + /// get the assembled heat flux AKANTU_GET_MACRO(HeatFlux,* heat_flux, Vector&); -/// set the value of the time step + /// get the lumped capacity AKANTU_GET_MACRO(CapacityLumped, * capacity_lumped, Vector&); + /// get the boundary vector AKANTU_GET_MACRO(Boundary, * boundary, Vector&); - /// get the SolidMechanicsModel::velocity vector + /// get the temperature gradient AKANTU_GET_MACRO_BY_ELEMENT_TYPE(TemperatureGradient, temperature_gradient, Vector &); - /// get the SolidMechanicsModel::velocity vector + /// get the temperature AKANTU_GET_MACRO(Temperature, *temperature, Vector &); - - /// get the equation number Vector + /// get the equation number Vector AKANTU_GET_MACRO(EquationNumber, *equation_number, const Vector &); - - // AKANTU_GET_MACRO(HeatFlux, *heat_flux, Vector &); -/// compute the stable time step - Real getStableTimeStep(); - - + /* ------------------------------------------------------------------------ */ + /* Methods */ + /* ------------------------------------------------------------------------ */ +public: - //initialize the heat flux + /// compute and get the stable time step + Real getStableTimeStep(); + //initialize the heat flux void initializeHeatFlux(Vector &temp); //initialize temperature void initializeTemperature(Vector &temp); //set boundary condition void setBoundaryCondition(); //compute temperature gradient void computeTemperatureGradient(const ElementType &el_type); - - //compute the heat flux + //compute the heat flux void updateHeatFlux(); - //compute the temperature void updateTemperature(); - //put the scheme into iteration void integrationScheme1stOrder(Real thelta, UInt N, Vector * temperature); - //calculate the capacity matrix of heat transfer problem void assembleCapacityLumped(const ElementType &el_type); - - - - - - /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// time step Real time_step; /// temperatures array Vector * temperature ; /// the spatial dimension UInt spatial_dimension; /// the density Real density; /// the speed of the changing temperature ByElementTypeReal temperature_gradient; /// K*T internal flux vector Vector * heat_flux; //external flux vector Vector * externalFlux; /// 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; //the biggest parameter of conductivity matrix Real conductivitymax; - }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "heat_transfer_model_inline_impl.cc" /// 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/model/heat_transfer/heat_transfer_model_inline_impl.cc b/model/heat_transfer/heat_transfer_model_inline_impl.cc index 1214c4071..a4eab06ed 100644 --- a/model/heat_transfer/heat_transfer_model_inline_impl.cc +++ b/model/heat_transfer/heat_transfer_model_inline_impl.cc @@ -1,28 +1,115 @@ /** * @file heat_transfer_model_inline_impl.cc * @author Srinivasa Babu Ramisetti * @date Fri Mar 4 17:04:25 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 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) 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(); +}; +/* -------------------------------------------------------------------------- */ + diff --git a/model/model.hh b/model/model.hh index ccb68bc4a..c15aa6876 100644 --- a/model/model.hh +++ b/model/model.hh @@ -1,158 +1,166 @@ /** * @file model.hh * @author Nicolas Richart * @date Thu Jul 22 11:02:42 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: inline Model(const ModelID & id = "model", const MemoryID & memory_id = 0); inline 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 + virtual void initPBC(UInt x, UInt y, UInt z); + /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// 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(){return *synch_registry;}; + SynchronizerRegistry & getSynchronizerRegistry(); protected: /// 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 = ""); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id ModelID 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; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "model_inline_impl.cc" /// 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/model/model_inline_impl.cc b/model/model_inline_impl.cc index 299157163..e1a7e1057 100644 --- a/model/model_inline_impl.cc +++ b/model/model_inline_impl.cc @@ -1,194 +1,243 @@ /** * @file model_inline_impl.cc * @author Guillaume ANCIAUX * @date Fri Aug 20 17:18:08 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 Model::Model(const ModelID & id, const MemoryID & memory_id) : Memory(memory_id), id(id),synch_registry(NULL) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } - +/* -------------------------------------------------------------------------- */ +inline SynchronizerRegistry & Model::getSynchronizerRegistry(){ + AKANTU_DEBUG_ASSERT(synch_registry,"synchronizer registry not initialized:" + << " did you call createSynchronizerRegistry ?"); + return *synch_registry; +} /* -------------------------------------------------------------------------- */ inline void Model::createSynchronizerRegistry(DataAccessor * data_accessor){ synch_registry = new SynchronizerRegistry(*data_accessor); } +/* -------------------------------------------------------------------------- */ +inline void Model::initPBC(UInt x, UInt y, UInt z){ + Mesh & mesh = getFEM().getMesh(); + mesh.computeBoundingBox(); + if (x) MeshUtils::computePBCMap(mesh,0,pbc_pair); + if (y) MeshUtils::computePBCMap(mesh,1,pbc_pair); + if (z) MeshUtils::computePBCMap(mesh,2,pbc_pair); + + std::map::iterator it = pbc_pair.begin(); + std::map::iterator end = pbc_pair.end(); + + Real * coords = mesh.getNodes().values; + UInt dim = mesh.getSpatialDimension(); + while(it != end){ + UInt i1 = (*it).first; + UInt i2 = (*it).second; + + AKANTU_DEBUG_INFO("pairing " << i1 << "(" + << coords[dim*i1] << "," << coords[dim*i1+1] << "," + << coords[dim*i1+2] + << ") with" + << i2 << "(" + << coords[dim*i2] << "," << coords[dim*i2+1] << "," + << coords[dim*i2+2] + << ")"); + ++it; + } + +} /* -------------------------------------------------------------------------- */ inline Synchronizer & Model::createParallelSynch(MeshPartition * partition, DataAccessor * data_accessor){ AKANTU_DEBUG_IN(); /* ------------------------------------------------------------------------ */ /* Parallel initialization */ /* ------------------------------------------------------------------------ */ StaticCommunicator * comm = StaticCommunicator::getStaticCommunicator(); Int prank = comm->whoAmI(); DistributedSynchronizer * synch = NULL; if(prank == 0) synch = DistributedSynchronizer::createDistributedSynchronizerMesh(getFEM().getMesh(), partition); else synch = DistributedSynchronizer::createDistributedSynchronizerMesh(getFEM().getMesh(), NULL); AKANTU_DEBUG_OUT(); return *synch; } /* -------------------------------------------------------------------------- */ inline Model::~Model() { AKANTU_DEBUG_IN(); FEMMap::iterator it; for (it = fems.begin(); it != fems.end(); ++it) { if(it->second) delete it->second; } for (it = fems_boundary.begin(); it != fems_boundary.end(); ++it) { if(it->second) delete it->second; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ 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; MeshUtils::buildFacets(it->second->getMesh()); 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; FEMMap::const_iterator it = fems.find(name); AKANTU_DEBUG_ASSERT(it == fems.end(), "FEM object with name " << name << " was already created"); 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); } +/* -------------------------------------------------------------------------- */ +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; + } + } +} +/* -------------------------------------------------------------------------- */ + diff --git a/model/solid_mechanics/solid_mechanics_model.cc b/model/solid_mechanics/solid_mechanics_model.cc index 1edac6ec3..c28c41628 100644 --- a/model/solid_mechanics/solid_mechanics_model.cc +++ b/model/solid_mechanics/solid_mechanics_model.cc @@ -1,851 +1,844 @@ /** * @file solid_mechanics_model.cc * @author Nicolas Richart * @date Thu Jul 22 14:35:38 2010 * * @brief Implementation of the SolidMechanicsModel class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" #include "aka_math.hh" #include "integration_scheme_2nd_order.hh" #include "static_communicator.hh" #include "sparse_matrix.hh" #include "solver.hh" #include "dof_synchronizer.hh" #ifdef AKANTU_USE_MUMPS #include "solver_mumps.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh, UInt dim, const ModelID & id, const MemoryID & memory_id) : Model(id, memory_id), time_step(NAN), f_m2a(1.0), mass_matrix(NULL), velocity_damping_matrix(NULL), stiffness_matrix(NULL),jacobian_matrix(NULL), integrator(NULL), increment_flag(false), solver(NULL), spatial_dimension(dim), mesh(mesh) { AKANTU_DEBUG_IN(); createSynchronizerRegistry(this); if (spatial_dimension == 0) spatial_dimension = mesh.getSpatialDimension(); registerFEMObject("SolidMechanicsFEM", mesh, spatial_dimension); this->displacement = NULL; this->mass = NULL; this->velocity = NULL; this->acceleration = NULL; this->force = NULL; this->residual = NULL; this->boundary = NULL; this->increment = NULL; this->increment_acceleration = NULL; for(UInt t = _not_defined; t < _max_element_type; ++t) { this->element_material[t] = NULL; this->ghost_element_material[t] = NULL; } this->dof_synchronizer = NULL; materials.clear(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModel::~SolidMechanicsModel() { AKANTU_DEBUG_IN(); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { delete *mat_it; } materials.clear(); delete integrator; if(solver) delete solver; if(mass_matrix) delete mass_matrix; if(velocity_damping_matrix) delete velocity_damping_matrix; if(stiffness_matrix) delete stiffness_matrix; if(jacobian_matrix) delete jacobian_matrix; if(dof_synchronizer) delete dof_synchronizer; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initParallel(MeshPartition * partition, DataAccessor * data_accessor) { AKANTU_DEBUG_IN(); if (data_accessor == NULL) data_accessor = this; Synchronizer & synch_parallel = createParallelSynch(partition,data_accessor); synch_registry->registerSynchronizer(synch_parallel,_gst_smm_mass); synch_registry->registerSynchronizer(synch_parallel,_gst_smm_for_strain); synch_registry->registerSynchronizer(synch_parallel,_gst_smm_boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initExplicit() { AKANTU_DEBUG_IN(); if (integrator) delete integrator; integrator = new CentralDifference(); // UInt nb_nodes = mesh.getNbNodes(); // std::stringstream sstr_eq; sstr_eq << id << ":equation_number"; // equation_number = &(alloc(sstr_eq.str(), nb_nodes, spatial_dimension, 0)); // Int * equation_number_val = equation_number->values; // for (UInt n = 0; n < nb_nodes; ++n) { // for (UInt d = 0; d < spatial_dimension; ++d) { // UInt eq_num = n * spatial_dimension + d; // *(equation_number_val++) = eq_num; // } // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initVectors() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_disp; sstr_disp << id << ":displacement"; std::stringstream sstr_mass; sstr_mass << id << ":mass"; std::stringstream sstr_velo; sstr_velo << id << ":velocity"; std::stringstream sstr_acce; sstr_acce << id << ":acceleration"; std::stringstream sstr_forc; sstr_forc << id << ":force"; std::stringstream sstr_resi; sstr_resi << id << ":residual"; std::stringstream sstr_boun; sstr_boun << id << ":boundary"; displacement = &(alloc(sstr_disp.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); mass = &(alloc(sstr_mass.str(), nb_nodes, spatial_dimension)); // \todo see if it must not be spatial_dimension velocity = &(alloc(sstr_velo.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); acceleration = &(alloc(sstr_acce.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); force = &(alloc(sstr_forc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); residual = &(alloc(sstr_resi.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); boundary = &(alloc(sstr_boun.str(), nb_nodes, spatial_dimension, false)); std::stringstream sstr_curp; sstr_curp << id << ":current_position"; current_position = &(alloc(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE)); const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { UInt nb_element = mesh.getNbElement(*it); if(!element_material[*it]) { std::stringstream sstr_elma; sstr_elma << id << ":element_material:" << *it; element_material[*it] = &(alloc(sstr_elma.str(), nb_element, 1, 0)); } } const Mesh::ConnectivityTypeList & ghost_type_list = mesh.getConnectivityTypeList(_ghost); for(it = ghost_type_list.begin(); it != ghost_type_list.end(); ++it) { UInt nb_element = mesh.getNbGhostElement(*it); std::stringstream sstr_elma; sstr_elma << id << ":ghost_element_material:" << *it; ghost_element_material[*it] = &(alloc(sstr_elma.str(), nb_element, 1, 0)); } dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initModel() { /// \todo add the current position as a parameter to initShapeFunctions for /// large deformation getFEM().initShapeFunctions(_not_ghost); getFEM().initShapeFunctions(_ghost); } - - - +/* -------------------------------------------------------------------------- */ +void SolidMechanicsModel::initPBC(UInt x, UInt y, UInt z){ + Model::initPBC(x,y,z); + PBCSynchronizer * synch = new PBCSynchronizer(pbc_pair); + synch_registry->registerSynchronizer(*synch,_gst_smm_uv); + synch_registry->registerSynchronizer(*synch,_gst_smm_mass); + changeLocalEquationNumberforPBC(pbc_pair,mesh.getSpatialDimension()); +} /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateCurrentPosition() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); current_position->resize(nb_nodes); //Vector * current_position = new Vector(nb_nodes, spatial_dimension, NAN, "position"); Real * current_position_val = current_position->values; Real * position_val = mesh.getNodes().values; Real * displacement_val = displacement->values; /// compute current_position = initial_position + displacement memcpy(current_position_val, position_val, nb_nodes*spatial_dimension*sizeof(Real)); for (UInt n = 0; n < nb_nodes*spatial_dimension; ++n) { *current_position_val++ += *displacement_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initializeUpdateResidualData() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); residual->resize(nb_nodes); /// copy the forces in residual for boundary conditions memcpy(residual->values, force->values, nb_nodes*spatial_dimension*sizeof(Real)); updateCurrentPosition(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Explicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateResidual(bool need_initialize) { AKANTU_DEBUG_IN(); if (need_initialize) initializeUpdateResidualData(); /// start synchronization synch_registry->asynchronousSynchronize(_gst_smm_for_strain); /// call update residual on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { //(*mat_it)->updateResidual(*current_position, _not_ghost); (*mat_it)->updateResidual(*displacement, _not_ghost); } /// finalize communications synch_registry->waitEndSynchronize(_gst_smm_for_strain); /// call update residual on each ghost elements for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { //(*mat_it)->updateResidual(*current_position, _ghost); (*mat_it)->updateResidual(*displacement, _ghost); } // current_position; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateAcceleration() { AKANTU_DEBUG_IN(); UInt nb_nodes = acceleration->getSize(); UInt nb_degre_of_freedom = acceleration->getNbComponent(); if(!increment_acceleration) increment_acceleration = new Vector(nb_nodes, nb_degre_of_freedom); increment_acceleration->resize(nb_nodes); increment_acceleration->clear(); Real * mass_val = mass->values; Real * residual_val = residual->values; Real * accel_val = acceleration->values; bool * boundary_val = boundary->values; Real * inc = increment_acceleration->values; for (UInt n = 0; n < nb_nodes; ++n) { for (UInt d = 0; d < nb_degre_of_freedom; d++) { if(!(*boundary_val)) { *inc = f_m2a * (*residual_val / *mass_val) - *accel_val; } residual_val++; accel_val++; boundary_val++; inc++; mass_val++; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitPred() { AKANTU_DEBUG_IN(); if(increment_flag) { memcpy(increment->values, displacement->values, displacement->getSize()*displacement->getNbComponent()*sizeof(Real)); } AKANTU_DEBUG_ASSERT(integrator,"itegrator should have been allocated: " << "have called initExplicit ? " << "or initImplicit ?"); integrator->integrationSchemePred(time_step, *displacement, *velocity, *acceleration, *boundary); if(increment_flag) { Real * inc_val = increment->values; Real * dis_val = displacement->values; UInt nb_nodes = displacement->getSize(); for (UInt n = 0; n < nb_nodes; ++n) { *inc_val = *dis_val - *inc_val; inc_val++; dis_val++; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitCorr() { AKANTU_DEBUG_IN(); Vector tmp(acceleration->getSize(), acceleration->getNbComponent()); tmp.clear(); integrator->integrationSchemeCorrAccel(time_step, *displacement, *velocity, *acceleration, *boundary, *increment_acceleration); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Implicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initImplicit(bool dynamic) { AKANTU_DEBUG_IN(); UInt nb_global_node = mesh.getNbGlobalNodes(); std::stringstream sstr; sstr << id << ":stiffness_matrix"; stiffness_matrix = new SparseMatrix(nb_global_node * spatial_dimension, _symmetric, spatial_dimension, sstr.str(), memory_id); // if(!dof_synchronizer) dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initGlobalDOFEquationNumbers(); stiffness_matrix->buildProfile(mesh, *dof_synchronizer); SparseMatrix * matrix; if(dynamic) { if(integrator) delete integrator; integrator = new TrapezoidalRule(); std::stringstream sstr_jac; sstr_jac << id << ":jacobian_matrix"; jacobian_matrix = new SparseMatrix(*stiffness_matrix, sstr_jac.str(), memory_id); // jacobian_matrix->buildProfile(mesh, *dof_synchronizer); matrix = jacobian_matrix; } else { matrix = stiffness_matrix; } #ifdef AKANTU_USE_MUMPS std::stringstream sstr_solv; sstr_solv << id << ":solver_stiffness_matrix"; solver = new SolverMumps(*matrix, sstr_solv.str()); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS solver->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initialAcceleration() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving Ma = f"); Solver * acc_solver; #ifdef AKANTU_USE_MUMPS std::stringstream sstr; sstr << id << ":solver_mass_matrix"; acc_solver = new SolverMumps(*mass_matrix, sstr.str()); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS acc_solver->initialize(); mass_matrix->applyBoundary(*boundary); acc_solver->setRHS(*residual); acc_solver->solve(*acceleration); delete acc_solver; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); stiffness_matrix->clear(); /// call compute stiffness matrix on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->assembleStiffnessMatrix(*displacement, _not_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveDynamic() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving Ma + Cv + Ku = f"); AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL, "You should first initialize the implicit solver and assemble the stiffness matrix"); AKANTU_DEBUG_ASSERT(mass_matrix != NULL, "You should first initialize the implicit solver and assemble the mass matrix"); NewmarkBeta * nmb_int = dynamic_cast(integrator); Real c = nmb_int->getAccelerationCoefficient(time_step); Real d = nmb_int->getVelocityCoefficient(time_step); Real e = nmb_int->getDisplacementCoefficient(time_step); // A = c M + d C + e K jacobian_matrix->clear(); jacobian_matrix->add(*stiffness_matrix, e); jacobian_matrix->add(*mass_matrix, c); if(velocity_damping_matrix) jacobian_matrix->add(*velocity_damping_matrix, d); jacobian_matrix->applyBoundary(*boundary); #ifndef AKANTU_NDEBUG if(AKANTU_DEBUG_TEST(dblDump)) jacobian_matrix->saveMatrix("J.mtx"); #endif // f = f_ext - f_int - Ma - Cv Vector * Ma = new Vector(*acceleration, true, "Ma"); *Ma *= *mass_matrix; *residual -= *Ma; delete Ma; if(velocity_damping_matrix) { Vector * Cv = new Vector(*velocity); *Cv *= *velocity_damping_matrix; *residual -= *Cv; delete Cv; } solver->setRHS(*residual); if(!increment) setIncrementFlagOn(); // solve A w = f solver->solve(*increment); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveStatic() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving Ku = f"); AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL, "You should first initialize the implicit solver and assemble the stiffness matrix"); UInt nb_nodes = displacement->getSize(); UInt nb_degre_of_freedom = displacement->getNbComponent(); stiffness_matrix->applyBoundary(*boundary); solver->setRHS(*residual); if(!increment) setIncrementFlagOn(); solver->solve(*increment); Real * increment_val = increment->values; Real * displacement_val = displacement->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes * nb_degre_of_freedom; ++n) { if(!(*boundary_val)) { *displacement_val += *increment_val; } displacement_val++; boundary_val++; increment_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance) { Real error; bool tmp = testConvergenceIncrement(tolerance, error); AKANTU_DEBUG_INFO("Norm of increment : " << error); return tmp; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance, Real & error) { AKANTU_DEBUG_IN(); UInt nb_nodes = displacement->getSize(); UInt nb_degre_of_freedom = displacement->getNbComponent(); Real norm = 0; Real * increment_val = increment->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); for (UInt d = 0; d < nb_degre_of_freedom; ++d) { if(!(*boundary_val) && is_local_node) { norm += *increment_val * *increment_val; } boundary_val++; increment_val++; } } synch_registry->allReduce(&norm, _so_sum); error = sqrt(norm); AKANTU_DEBUG_ASSERT(!isnan(norm), "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); return (error < tolerance); } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceResidual(Real tolerance) { Real error; bool tmp = testConvergenceResidual(tolerance, error); AKANTU_DEBUG_INFO("Norm of residual : " << error); return tmp; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceResidual(Real tolerance, Real & error) { AKANTU_DEBUG_IN(); UInt nb_nodes = residual->getSize(); Real norm = 0; Real * residual_val = residual->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); if(is_local_node) { for (UInt d = 0; d < spatial_dimension; ++d) { if(!(*boundary_val)) { norm += *residual_val * *residual_val; } boundary_val++; residual_val++; } } else { boundary_val += spatial_dimension; residual_val += spatial_dimension; } } synch_registry->allReduce(&norm, _so_sum); error = sqrt(norm); AKANTU_DEBUG_ASSERT(!isnan(norm), "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); return (error < tolerance); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::implicitPred() { AKANTU_DEBUG_IN(); integrator->integrationSchemePred(time_step, *displacement, *velocity, *acceleration, *boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::implicitCorr() { AKANTU_DEBUG_IN(); integrator->integrationSchemeCorrDispl(time_step, *displacement, *velocity, *acceleration, *boundary, *increment); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Information */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::synchronizeBoundaries() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(synch_registry,"Synchronizer registry was not initialized." << " Did you call initParallel ?"); synch_registry->synchronize(_gst_smm_boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setIncrementFlagOn() { AKANTU_DEBUG_IN(); if(!increment) { UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_inc; sstr_inc << id << ":increment"; increment = &(alloc(sstr_inc.str(), nb_nodes, spatial_dimension, 0)); } increment_flag = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Material ** mat_val = &(materials.at(0)); Real min_dt = std::numeric_limits::max(); Real * coord = mesh.getNodes().values; Real * disp_val = displacement->values; const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; Element elem; for(it = type_list.begin(); it != type_list.end(); ++it) { if(mesh.getSpatialDimension(*it) != spatial_dimension) continue; elem.type = *it; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*it); UInt nb_element = mesh.getNbElement(*it); UInt * conn = mesh.getConnectivity(*it).values; UInt * elem_mat_val = element_material[*it]->values; Real * u = new Real[nb_nodes_per_element*spatial_dimension]; for (UInt el = 0; el < nb_element; ++el) { UInt el_offset = el * nb_nodes_per_element; elem.element = el; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_conn = conn[el_offset + n] * spatial_dimension; memcpy(u + n * spatial_dimension, coord + offset_conn, spatial_dimension * sizeof(Real)); for (UInt i = 0; i < spatial_dimension; ++i) { u[n * spatial_dimension + i] += disp_val[offset_conn + i]; } } Real el_size = getFEM().getElementInradius(u, *it); Real el_dt = mat_val[elem_mat_val[el]]->getStableTimeStep(el_size, elem); min_dt = min_dt > el_dt ? el_dt : min_dt; } delete [] u; } /// reduction min over all processors synch_registry->allReduce(&min_dt, _so_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real epot = 0.; /// call update residual on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { epot += (*mat_it)->getPotentialEnergy(); } /// reduction sum over all processors synch_registry->allReduce(&epot, _so_sum); AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy() { AKANTU_DEBUG_IN(); Real ekin = 0.; UInt nb_nodes = mesh.getNbNodes(); Real * vel_val = velocity->values; Real * mass_val = mass->values; for (UInt n = 0; n < nb_nodes; ++n) { Real v2 = 0; bool is_local_node = mesh.isLocalOrMasterNode(n); for (UInt i = 0; i < spatial_dimension; ++i) { if(is_local_node) { v2 += *vel_val * *vel_val; } vel_val++; } ekin += *mass_val * v2; mass_val++; } synch_registry->allReduce(&ekin, _so_sum); AKANTU_DEBUG_OUT(); return ekin * .5; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Solid Mechanics Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEM().printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + nodals information [" << std::endl; displacement->printself(stream, indent + 2); mass ->printself(stream, indent + 2); velocity ->printself(stream, indent + 2); acceleration->printself(stream, indent + 2); force ->printself(stream, indent + 2); residual ->printself(stream, indent + 2); boundary ->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + connectivity type information [" << std::endl; const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { stream << space << AKANTU_INDENT << AKANTU_INDENT << " + " << *it <<" [" << std::endl; element_material[*it]->printself(stream, indent + 3); stream << space << AKANTU_INDENT << AKANTU_INDENT << "]" << std::endl; } stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } -/* -------------------------------------------------------------------------- */ -void SolidMechanicsModel::changeEquationNumberforPBC(std::map & pbc_pair){ - 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 < spatial_dimension; ++i) { - (*dof_synchronizer->getLocalDOFEquationNumbersPointer())(node_slave*spatial_dimension+i) = spatial_dimension*node_master+i; - } - } -} - /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/model/solid_mechanics/solid_mechanics_model.hh b/model/solid_mechanics/solid_mechanics_model.hh index f5afe3b94..d747a8c0b 100644 --- a/model/solid_mechanics/solid_mechanics_model.hh +++ b/model/solid_mechanics/solid_mechanics_model.hh @@ -1,447 +1,450 @@ /** * @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 * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_HH__ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "model.hh" #include "material.hh" #include "../parser.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #include "aka_types.hh" /* -------------------------------------------------------------------------- */ namespace akantu { // class Material; class IntegrationScheme2ndOrder; class Contact; class Solver; class SparseMatrix; } __BEGIN_AKANTU__ class SolidMechanicsModel : public Model, public DataAccessor { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef FEMTemplate MyFEMType; // 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); SolidMechanicsModel(Mesh & mesh,MeshPartition * partition, UInt spatial_dimension = 0, const ModelID & id = "solid_mechanics_model", const MemoryID & memory_id = 0); virtual ~SolidMechanicsModel(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// register the tags associated with the parallel synchronizer void initParallel(MeshPartition * partition, DataAccessor * data_accessor=NULL); /// allocate all vectors void initVectors(); /// initialize all internal arrays for materials void initMaterials(); /// initialize the model void initModel(); + /// init PBC synchronizer + void initPBC(UInt x,UInt y, UInt z); + /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* PBC */ /* ------------------------------------------------------------------------ */ public: /// change the equation number for proper assembly when using PBC void changeEquationNumberforPBC(std::map & pbc_pair); /* ------------------------------------------------------------------------ */ /* Explicit */ /* ------------------------------------------------------------------------ */ public: /// initialize the array needed by updateResidual (residual, current_position) void initializeUpdateResidualData(); /// update the current position vector void updateCurrentPosition(); /// assemble the residual for the explicit scheme void updateResidual(bool need_initialize = true); /// compute the acceleration from the residual void updateAcceleration(); /// explicit integration predictor void explicitPred(); /// explicit integration corrector void explicitCorr(); /// initialize the stuff for the explicit scheme void initExplicit(); /* ------------------------------------------------------------------------ */ /* Implicit */ /* ------------------------------------------------------------------------ */ public: /// initialize the stuff for the implicit solver void initImplicit(bool dynamic = false); /// solve Ma = f to get the initial acceleration void initialAcceleration(); /// assemble the stiffness matrix void assembleStiffnessMatrix(); /// solve @f[ A\delta u = f_ext - f_int @f] void solveDynamic(); /// solve Ku = f void solveStatic(); /// test the convergence (norm of increment) bool testConvergenceIncrement(Real tolerance); bool testConvergenceIncrement(Real tolerance, Real & error); /// test the convergence (norm of residual) bool testConvergenceResidual(Real tolerance); bool testConvergenceResidual(Real tolerance, Real & error); /// implicit time integration predictor void implicitPred(); /// implicit time integration corrector void implicitCorr(); /* ------------------------------------------------------------------------ */ /* Boundaries (solid_mechanics_model_boundary.cc) */ /* ------------------------------------------------------------------------ */ public: /// 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(BoundaryFunction in_function, BoundaryFunctionType function_type); /// synchronize the ghost element boundaries values void synchronizeBoundaries(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: /// 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); private: /// read properties part of a material file and create the material template Material * readMaterialProperties(std::ifstream & infile, MaterialID mat_id, UInt ¤t_line); /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix void assembleMass(); private: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); void assembleMass(GhostType ghost_type); /// fill a vector of rho void computeRho(Vector & rho, ElementType type, GhostType ghost_type); /* ------------------------------------------------------------------------ */ - /* Ghost Synchronizer inherited members */ + /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: inline virtual UInt getNbDataToPack(const Element & element, SynchronizationTag tag) const; inline virtual UInt getNbDataToUnpack(const Element & element, SynchronizationTag tag) const; inline virtual void packData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag) const; inline virtual void unpackData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag) const; inline virtual UInt getNbDataToPack(SynchronizationTag tag) const; inline virtual UInt getNbDataToUnpack(SynchronizationTag tag) const; inline virtual void packData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag) const; inline virtual void unpackData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the current value of the time step AKANTU_GET_MACRO(TimeStep, time_step, Real); /// set the value of the time step AKANTU_SET_MACRO(TimeStep, time_step, Real); /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement, Vector &); /** * @brief get the SolidMechanicsModel::current_position vector \warn only * consistent after a call to SolidMechanicsModel::updateCurrentPosition */ AKANTU_GET_MACRO(CurrentPosition, *current_position, const Vector &); /** * @brief get the SolidMechanicsModel::increment vector \warn only consistent * if SolidMechanicsModel::setIncrementFlagOn has been called before */ AKANTU_GET_MACRO(Increment, *increment, const Vector &); /// get the lumped SolidMechanicsModel::mass vector AKANTU_GET_MACRO(Mass, *mass, const Vector &); /// get the SolidMechanicsModel::velocity vector AKANTU_GET_MACRO(Velocity, *velocity, Vector &); /** * @brief get the SolidMechanicsModel::acceleration vector, updated by * SolidMechanicsModel::updateAcceleration */ AKANTU_GET_MACRO(Acceleration, *acceleration, Vector &); /// get the SolidMechanicsModel::force vector (boundary forces) AKANTU_GET_MACRO(Force, *force, Vector &); /// get the SolidMechanicsModel::residual vector, computed by SolidMechanicsModel::updateResidual AKANTU_GET_MACRO(Residual, *residual, const Vector &); /// get the SolidMechanicsModel::boundary vector AKANTU_GET_MACRO(Boundary, *boundary, Vector &); /// get the value of the SolidMechanicsModel::increment_flag AKANTU_GET_MACRO(IncrementFlag, increment_flag, bool); /// get the SolidMechanicsModel::element_material vector corresponding to a given akantu::ElementType AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementMaterial, element_material, Vector &); /// get a particular material inline Material & getMaterial(UInt mat_index); inline const Material & getMaterial(UInt mat_index) const; /// give the number of materials inline UInt getNbMaterials() { return materials.size(); }; /// compute the stable time step Real getStableTimeStep(); /// compute the potential energy Real getPotentialEnergy(); /// compute the kinetic energy Real getKineticEnergy(); /// set the Contact object AKANTU_SET_MACRO(Contact, contact, Contact *); /** * @brief set the SolidMechanicsModel::increment_flag to on, the activate the * update of the SolidMechanicsModel::increment vector */ void setIncrementFlagOn(); /// get the stiffness matrix AKANTU_GET_MACRO(StiffnessMatrix, *stiffness_matrix, SparseMatrix &); /// get the mass matrix AKANTU_GET_MACRO(MassMatrix, *mass_matrix, SparseMatrix &); inline FEM & getFEMBoundary(std::string name = ""); /* ------------------------------------------------------------------------ */ /* 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; /// accelerations array Vector * increment_acceleration; /// forces array Vector * force; /// residuals array Vector * residual; /// boundaries array Vector * boundary; /// array of current position used during update residual Vector * current_position; /// mass matrix SparseMatrix * mass_matrix; /// velocity damping matrix SparseMatrix * velocity_damping_matrix; /// stiffness matrix SparseMatrix * stiffness_matrix; /// jacobian matrix @f[A = cM + dD + K@f] with @f[c = \frac{1}{\beta \Dela t^2}, d = \frac{\gamma}{\beta \Delta t} @f] SparseMatrix * jacobian_matrix; /// 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; /// solver for implicit Solver * solver; /// object to resolve the contact Contact * contact; /// the spatial dimension UInt spatial_dimension; /// Mesh Mesh & mesh; }; /* -------------------------------------------------------------------------- */ /* 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_inline_impl.cc b/model/solid_mechanics/solid_mechanics_model_inline_impl.cc index f873f5d37..66b8a34e6 100644 --- a/model/solid_mechanics/solid_mechanics_model_inline_impl.cc +++ b/model/solid_mechanics/solid_mechanics_model_inline_impl.cc @@ -1,389 +1,389 @@ /** * @file solid_mechanics_model_inline_impl.cc * @author Nicolas Richart * @date Thu Jul 29 12:07:04 2010 * * @brief Implementation of the inline functions of the SolidMechanicsModel class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ inline Material & SolidMechanicsModel::getMaterial(UInt mat_index) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(mat_index < materials.size(), "The model " << id << " has no material no "<< mat_index); AKANTU_DEBUG_OUT(); return *materials[mat_index]; } /* -------------------------------------------------------------------------- */ inline const Material & SolidMechanicsModel::getMaterial(UInt mat_index) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(mat_index < materials.size(), "The model " << id << " has no material no "<< mat_index); AKANTU_DEBUG_OUT(); return *materials[mat_index]; } /* -------------------------------------------------------------------------- */ template UInt SolidMechanicsModel::readCustomMaterial(const std::string & filename, const std::string & keyword) { Parser parser; parser.open(filename); std::string key = keyword; to_lower(key); std::string mat_name = parser.getNextSection("material"); while (mat_name != ""){ if (mat_name == key) break; mat_name = parser.getNextSection("material"); } if (mat_name != key) AKANTU_DEBUG_ERROR("material " << key << " not found in file " << filename); std::stringstream sstr_mat; sstr_mat << id << ":" << materials.size() << ":" << key; MaterialID mat_id = sstr_mat.str(); Material * mat = parser.readSection(*this,mat_id); materials.push_back(mat); return materials.size();; } /* -------------------------------------------------------------------------- */ inline FEM & SolidMechanicsModel::getFEMBoundary(std::string name) { return dynamic_cast(getFEMClassBoundary(name)); } /* -------------------------------------------------------------------------- */ inline UInt SolidMechanicsModel::getNbDataToPack(const Element & element, SynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); #ifdef AKANTU_DEBUG size += spatial_dimension * sizeof(Real); /// position of the barycenter of the element (only for check) #endif switch(tag) { case _gst_smm_mass: { size += nb_nodes_per_element * sizeof(Real); // mass vector break; } case _gst_smm_for_strain: { size += nb_nodes_per_element * spatial_dimension * sizeof(Real); // displacement UInt mat = element_material[element.type]->values[element.element]; size += materials[mat]->getNbDataToPack(element, tag); break; } case _gst_smm_boundary: { // force, displacement, boundary size += nb_nodes_per_element * spatial_dimension * (2 * sizeof(Real) + sizeof(bool)); break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ inline UInt SolidMechanicsModel::getNbDataToUnpack(const Element & element, SynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); #ifdef AKANTU_DEBUG size += spatial_dimension * sizeof(Real); /// position of the barycenter of the element (only for check) #endif switch(tag) { case _gst_smm_mass: { size += nb_nodes_per_element * sizeof(Real); // mass vector break; } case _gst_smm_for_strain: { size += nb_nodes_per_element * spatial_dimension * sizeof(Real); // displacement UInt mat = ghost_element_material[element.type]->values[element.element]; size += materials[mat]->getNbDataToPack(element, tag); break; } case _gst_smm_boundary: { // force, displacement, boundary size += nb_nodes_per_element * spatial_dimension * (2 * sizeof(Real) + sizeof(bool)); break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::packData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); UInt el_offset = element.element * nb_nodes_per_element; UInt * conn = mesh.getConnectivity(element.type).values; #ifdef AKANTU_DEBUG types::RVector barycenter(spatial_dimension); mesh.getBarycenter(element.element, element.type, barycenter.storage()); buffer << barycenter; #endif switch(tag) { case _gst_smm_mass: { for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_conn = conn[el_offset + n]; buffer << (*mass)(offset_conn); } break; } case _gst_smm_for_strain: { Vector::iterator it_disp = displacement->begin(spatial_dimension); for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_conn = conn[el_offset + n]; buffer << it_disp[offset_conn]; } UInt mat = element_material[element.type]->values[element.element]; materials[mat]->packData(buffer, element, tag); break; } case _gst_smm_boundary: { Vector::iterator it_force = force->begin(spatial_dimension); Vector::iterator it_velocity = velocity->begin(spatial_dimension); Vector::iterator > it_boundary = boundary->begin(spatial_dimension); for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_conn = conn[el_offset + n]; buffer << it_force [offset_conn]; buffer << it_velocity[offset_conn]; buffer << it_boundary[offset_conn]; } break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); UInt el_offset = element.element * nb_nodes_per_element; UInt * conn = mesh.getGhostConnectivity(element.type).values; #ifdef AKANTU_DEBUG types::RVector barycenter_loc(spatial_dimension); mesh.getBarycenter(element.element, element.type, barycenter_loc.storage(), _ghost); 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) << ")"); } #endif switch(tag) { case _gst_smm_mass: { for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_conn = conn[el_offset + n]; buffer >> (*mass)(offset_conn); } break; } case _gst_smm_for_strain: { Vector::iterator it_disp = displacement->begin(spatial_dimension); for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_conn = conn[el_offset + n]; buffer >> it_disp[offset_conn]; } UInt mat = ghost_element_material[element.type]->values[element.element]; materials[mat]->unpackData(buffer, element, tag); break; } case _gst_smm_boundary: { Vector::iterator it_force = force->begin(spatial_dimension); Vector::iterator it_velocity = velocity->begin(spatial_dimension); Vector::iterator > it_boundary = boundary->begin(spatial_dimension); for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_conn = conn[el_offset + n]; buffer >> it_force [offset_conn]; buffer >> it_velocity[offset_conn]; buffer >> it_boundary[offset_conn]; } break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline UInt SolidMechanicsModel::getNbDataToPack(SynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes = mesh.getNbNodes(); switch(tag) { case _gst_smm_uv: { size += nb_nodes * sizeof(Real) * spatial_dimension * 2; break; } case _gst_smm_mass: { size += nb_nodes * sizeof(Real); break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ inline UInt SolidMechanicsModel::getNbDataToUnpack(SynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes = mesh.getNbNodes(); switch(tag) { case _gst_smm_uv: { size += nb_nodes * sizeof(Real) * spatial_dimension * 2; break; } case _gst_smm_mass: { size += nb_nodes * sizeof(Real); break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::packData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag) const { AKANTU_DEBUG_IN(); switch(tag) { case _gst_smm_uv: { for (UInt d = 0; d < spatial_dimension; ++d) { buffer << (*displacement)(index,d); buffer << (*velocity)(index,d); } break; } case _gst_smm_mass: { - AKANTU_DEBUG_INFO("pack mass of node " << index << " which is " << (*mass)(index)); + AKANTU_DEBUG_INFO("pack mass of node " << index << " which is " << (*mass)(index,0)); buffer << (*mass)(index,0); buffer << (*mass)(index,1); buffer << (*mass)(index,2); break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag) const { AKANTU_DEBUG_IN(); switch(tag) { case _gst_smm_uv: { for (UInt d = 0; d < spatial_dimension; ++d) { buffer >> (*displacement)(index,d); buffer >> (*velocity)(index,d); } break; } case _gst_smm_mass: { - AKANTU_DEBUG_INFO("mass of node " << index << " was " << (*mass)(index)); + AKANTU_DEBUG_INFO("mass of node " << index << " was " << (*mass)(index,0)); buffer >> (*mass)(index,0); buffer >> (*mass)(index,1); buffer >> (*mass)(index,2); - AKANTU_DEBUG_INFO("mass of node " << index << " is now " << (*mass)(index)); + AKANTU_DEBUG_INFO("mass of node " << index << " is now " << (*mass)(index,0)); break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ diff --git a/model/solid_mechanics/solid_mechanics_model_mass.cc b/model/solid_mechanics/solid_mechanics_model_mass.cc index 7953ff184..3fb29bc60 100644 --- a/model/solid_mechanics/solid_mechanics_model_mass.cc +++ b/model/solid_mechanics/solid_mechanics_model_mass.cc @@ -1,166 +1,165 @@ /** * @file solid_mechanics_model_mass.cc * @author Nicolas Richart * @date Mon Oct 4 15:21:23 2010 * * @brief function handling mass computation * * @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 "solid_mechanics_model.hh" #include "material.hh" - -//#include "static_communicator.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMassLumped() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); mass->clear(); assembleMassLumped(_not_ghost); assembleMassLumped(_ghost); /// for not connected nodes put mass to one in order to avoid /// wrong range in paraview Real * mass_values = mass->values; for (UInt i = 0; i < nb_nodes; ++i) { if (fabs(mass_values[i]) < std::numeric_limits::epsilon() || Math::isnan(mass_values[i])) mass_values[i] = 1.; } + synch_registry->synchronize(_gst_smm_mass); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMassLumped(GhostType ghost_type) { AKANTU_DEBUG_IN(); FEM & fem = getFEM(); Vector rho_1(0,1); const Mesh::ConnectivityTypeList & type_list = fem.getMesh().getConnectivityTypeList(ghost_type); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; ElementType type = *it; computeRho(rho_1, type, ghost_type); AKANTU_DEBUG_ASSERT(dof_synchronizer, "DOFSynchronizer number must not be initialized"); fem.assembleFieldLumped(rho_1, spatial_dimension,*mass, dof_synchronizer->getLocalDOFEquationNumbers(), type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMass() { AKANTU_DEBUG_IN(); if(!dof_synchronizer) dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); UInt nb_global_node = mesh.getNbGlobalNodes(); std::stringstream sstr; sstr << id << ":mass_matrix"; mass_matrix = new SparseMatrix(nb_global_node * spatial_dimension, _symmetric, spatial_dimension, sstr.str(), memory_id); mass_matrix->buildProfile(mesh, *dof_synchronizer); assembleMass(_not_ghost); // assembleMass(_ghost); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMass(GhostType ghost_type) { AKANTU_DEBUG_IN(); MyFEMType & fem = getFEMClass(); Vector rho_1(0,1); //UInt nb_element; mass_matrix->clear(); const Mesh::ConnectivityTypeList & type_list = fem.getMesh().getConnectivityTypeList(ghost_type); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; ElementType type = *it; computeRho(rho_1, type, ghost_type); fem.assembleFieldMatrix(rho_1, spatial_dimension, *mass_matrix, type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeRho(Vector & rho, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Material ** mat_val = &(materials.at(0)); UInt * elem_mat_val; UInt nb_element; FEM & fem = getFEM(); if(ghost_type == _not_ghost) { nb_element = fem.getMesh().getNbElement(type); elem_mat_val = element_material[type]->values; } else { nb_element = fem.getMesh().getNbGhostElement(type); elem_mat_val = ghost_element_material[type]->values; } UInt nb_quadrature_points = fem.getNbQuadraturePoints(type); rho.resize(nb_element * nb_quadrature_points); Real * rho_1_val = rho.values; /// compute @f$ rho @f$ for each nodes of each element for (UInt el = 0; el < nb_element; ++el) { Real mat_rho = mat_val[elem_mat_val[el]]->getRho(); /// here rho is constant in an element for (UInt n = 0; n < nb_quadrature_points; ++n) { *rho_1_val++ = mat_rho; } } AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/test/test_mesh_utils/test_pbc_tweak/CMakeLists.txt b/test/test_mesh_utils/test_pbc_tweak/CMakeLists.txt index 5f2974e39..4afb186ff 100644 --- a/test/test_mesh_utils/test_pbc_tweak/CMakeLists.txt +++ b/test/test_mesh_utils/test_pbc_tweak/CMakeLists.txt @@ -1,34 +1,37 @@ #=============================================================================== # @file CMakeLists.txt # @author Guillaume Anciaux # @date Fri Jun 11 09:46:59 2010 # # @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 # #=============================================================================== #=============================================================================== add_mesh(test_pbc_cube_mesh cube.geo 3 1) register_test(test_pbc_tweak test_pbc_tweak.cc) add_dependencies(test_pbc_tweak test_pbc_cube_mesh) +#=============================================================================== +file(COPY material.dat DESTINATION .) +file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/paraview) \ No newline at end of file diff --git a/test/test_mesh_utils/test_pbc_tweak/cube.geo b/test/test_mesh_utils/test_pbc_tweak/cube.geo index 0fca6b7a9..82c9c1e59 100644 --- a/test/test_mesh_utils/test_pbc_tweak/cube.geo +++ b/test/test_mesh_utils/test_pbc_tweak/cube.geo @@ -1,37 +1,43 @@ -h = 1; +h = 0.1; Point(1) = {0.0, 0.0, 0.0, h}; Point(2) = {1.0, 0.0, 0.0, h}; Point(3) = {0.0, 1.0, 0.0, h}; Point(4) = {1.0, 1.0, 0.0, h}; Point(5) = {0.0, 0.0, 1.0, h}; Point(6) = {1.0, 0.0, 1.0, h}; Point(7) = {0.0, 1.0, 1.0, h}; Point(8) = {1.0, 1.0, 1.0, h}; Line(1) = {7, 8}; Line(2) = {8, 6}; Line(3) = {6, 5}; Line(4) = {5, 7}; Line(5) = {3, 7}; Line(6) = {5, 1}; Line(7) = {1, 3}; Line(8) = {2, 1}; Line(9) = {3, 4}; Line(10) = {2, 4}; Line(11) = {2, 6}; Line(12) = {8, 4}; Line Loop(13) = {3, 6, -8, 11}; Plane Surface(14) = {13}; Line Loop(15) = {4, -5, -7, -6}; Plane Surface(16) = {15}; Line Loop(17) = {1, 12, -9, 5}; Plane Surface(18) = {17}; Line Loop(19) = {1, 2, 3, 4}; Plane Surface(20) = {19}; Line Loop(21) = {12, -10, 11, -2}; Plane Surface(22) = {21}; Line Loop(23) = {9, -10, 8, 7}; Plane Surface(24) = {23}; Surface Loop(25) = {24, 18, 20, 22, 14, 16}; Volume(26) = {25}; + +Transfinite Surface "*"; +Recombine Surface "*"; +Transfinite Volume "*"; + +Physical Volume (100) = {26}; \ No newline at end of file diff --git a/test/test_mesh_utils/test_pbc_tweak/material.dat b/test/test_mesh_utils/test_pbc_tweak/material.dat new file mode 100644 index 000000000..a5a80c228 --- /dev/null +++ b/test/test_mesh_utils/test_pbc_tweak/material.dat @@ -0,0 +1,13 @@ +material elastic [ + name = steel + rho = 7800 # density + E = 2.1e11 # young's modulus + nu = 0.3 # poisson's ratio +] + +material elastic [ + name = copper + rho = 8940 # density + E = 125e9 # young's modulus + nu = 0.34 # poisson's ratio +] diff --git a/test/test_mesh_utils/test_pbc_tweak/test_pbc_tweak.cc b/test/test_mesh_utils/test_pbc_tweak/test_pbc_tweak.cc index 88305e243..63297a91a 100644 --- a/test/test_mesh_utils/test_pbc_tweak/test_pbc_tweak.cc +++ b/test/test_mesh_utils/test_pbc_tweak/test_pbc_tweak.cc @@ -1,124 +1,92 @@ /** * @file test_facet_extraction_tetra1.cc * @author Guillaume ANCIAUX * @date Thu Aug 19 13:05:27 2010 * * @brief test of internal facet extraction * * @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 "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "pbc_synchronizer.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; int main(int argc, char *argv[]) { int dim = 3; akantu::initialize(&argc, &argv); akantu::debug::setDebugLevel(akantu::dblInfo); Mesh mesh(dim); MeshIOMSH mesh_io; mesh_io.read("cube.msh", mesh); - mesh.computeBoundingBox(); - std::map pbc_pair; - akantu::Math::tolerance = 1e-3; - MeshUtils::computePBCMap(mesh,0,pbc_pair); - MeshUtils::computePBCMap(mesh,1,pbc_pair); - MeshUtils::computePBCMap(mesh,2,pbc_pair); - { - std::map::iterator it = pbc_pair.begin(); - std::map::iterator end = pbc_pair.end(); - - Real * coords = mesh.getNodes().values; - UInt dim = mesh.getSpatialDimension(); - while(it != end){ - UInt i1 = (*it).first; - UInt i2 = (*it).second; - - AKANTU_DEBUG_INFO("pairing " << i1 << "(" - << coords[dim*i1] << "," << coords[dim*i1+1] << "," - << coords[dim*i1+2] - << ") with" - << i2 << "(" - << coords[dim*i2] << "," << coords[dim*i2+1] << "," - << coords[dim*i2+2] - << ")"); - ++it; - } - } - - PBCSynchronizer synch(pbc_pair); - SolidMechanicsModel * model = new SolidMechanicsModel(mesh); - model->getSynchronizerRegistry().registerSynchronizer(synch,_gst_smm_uv); - model->getSynchronizerRegistry().registerSynchronizer(synch,_gst_smm_mass); - + /* -------------------------------------------------------------------------- */ model->initVectors(); - UInt nb_nodes = model->getFEM().getMesh().getNbNodes(); - model->getForce().clear(); model->getVelocity().clear(); model->getAcceleration().clear(); model->getDisplacement().clear(); - + /* ------------------------------------------------------------------------ */ model->initExplicit(); model->initModel(); model->readMaterials("material.dat"); model->initMaterials(); - - model->changeLocalEquationNumberforPBC(pbc_pair,mesh.getSpatialDimension()); + /* -------------------------------------------------------------------------- */ + model->initPBC(1,1,1); model->assembleMassLumped(); - - #ifdef AKANTU_USE_IOHELPER + /* -------------------------------------------------------------------------- */ + +#ifdef AKANTU_USE_IOHELPER DumperParaview dumper; dumper.SetMode(TEXT); + UInt nb_nodes = model->getFEM().getMesh().getNbNodes(); dumper.SetPoints(mesh.getNodes().values, dim, nb_nodes, "test-pbc-tweak"); dumper.SetConnectivity((int*)mesh.getConnectivity(_hexahedron_8).values, HEX1, mesh.getNbElement(_hexahedron_8), C_MODE); dumper.AddNodeDataField(model->getMass().values, 3, "mass"); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); #endif //AKANTU_USE_IOHELPER return EXIT_SUCCESS; } diff --git a/test/test_model/test_heat_transfer_model/CMakeLists.txt b/test/test_model/test_heat_transfer_model/CMakeLists.txt index 2693f0a2f..c2dbee262 100644 --- a/test/test_model/test_heat_transfer_model/CMakeLists.txt +++ b/test/test_model/test_heat_transfer_model/CMakeLists.txt @@ -1,65 +1,65 @@ #=============================================================================== # @file CMakeLists.txt # @author Guillaume Anciaux # @date Fri Jun 11 09:46:59 2010 # # @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 # #=============================================================================== #=============================================================================== add_mesh(test_heat_cube3d_mesh1 cube.geo 3 1 OUTPUT cube1.msh) add_mesh(test_heat_square_mesh1 square.geo 2 1 OUTPUT square1.msh) #add_mesh(test_heat_bar1_mesh1 bar_structured.geo 2 1 OUTPUT bar1.msh) add_mesh(test_heat_line_mesh line.geo 1 1 OUTPUT line.msh) register_test(test_heat_transfer_model_cube3d test_heat_transfer_model_cube3d.cc) add_dependencies(test_heat_transfer_model_cube3d test_heat_cube3d_mesh1 test_heat_square_mesh1 test_heat_line_mesh) add_mesh(test_heat_cube_tet4 cube_tet.geo 3 1 OUTPUT cube_tet4.msh) register_test(test_heat_transfer_model_cube3d_pbc test_heat_transfer_model_cube3d_pbc.cc) add_dependencies(test_heat_transfer_model_cube3d_pbc test_heat_cube_tet4) add_mesh(test_heat_square_tri3 square_tri.geo 2 1 OUTPUT square_tri3.msh) register_test(test_heat_transfer_model_square2d_pbc test_heat_transfer_model_square2d_pbc.cc) -add_dependencies(test_heat_transfer_model_cube3d_pbc test_heat_square_tri4) +add_dependencies(test_heat_transfer_model_square2d_pbc test_heat_square_tri3) register_test(test_heat_transfer_model_cube3d_istropic_conductivity test_heat_transfer_model_cube3d_istropic_conductivity.cc) add_dependencies(test_heat_transfer_model_cube3d_istropic_conductivity test_heat_transfer_model_cube3d_istropic_conductivity_mesh1) #register_test(test_heat_transfer_model_cube3d_istropic_conductivity # test_heat_transfer_model_cube3d_anisotropy_conductivity.cc) #add_dependencies(test_heat_transfer_model_cube3d_anisotropy_conductivity test_heat_transfer_model_cube3d_anisotropy_conductivity_mesh1) # #register_test(test_heat_transfer_model_cube3d_tetra10 # test_heat_transfer_model_cube3d_tetra10.cc) #add_dependencies(test_heat_transfer_model_cube3d_tetra10 test_heat_cube3d_mesh2) #=============================================================================== file(COPY material.dat DESTINATION .) file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/paraview) diff --git a/test/test_model/test_heat_transfer_model/test_heat_transfer_model_cube3d_pbc.cc b/test/test_model/test_heat_transfer_model/test_heat_transfer_model_cube3d_pbc.cc index 1f7780bd9..d8d82649c 100644 --- a/test/test_model/test_heat_transfer_model/test_heat_transfer_model_cube3d_pbc.cc +++ b/test/test_model/test_heat_transfer_model/test_heat_transfer_model_cube3d_pbc.cc @@ -1,205 +1,169 @@ /** * @file test_heat_transfer_model_cube3d.cc * @author Rui WANG * @date Tue May 17 11:31:22 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; /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER #include "io_helper.h" #endif //AKANTU_USE_IOHELPER void paraviewInit(akantu::HeatTransferModel * model,Dumper & dumper); void paraviewDump(Dumper & dumper); akantu::UInt spatial_dimension = 3; akantu:: ElementType type = akantu::_tetrahedron_4; akantu::UInt paraview_type = TETRA1; /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { akantu::initialize(&argc,&argv); akantu::Mesh mesh(spatial_dimension); akantu::MeshIOMSH mesh_io; mesh_io.read("cube_tet4.msh", mesh); akantu::HeatTransferModel * model; akantu::UInt nb_nodes; akantu::UInt nb_element; - mesh.computeBoundingBox(); - std::map pbc_pair; - akantu::MeshUtils::computePBCMap(mesh,0,pbc_pair); - akantu::MeshUtils::computePBCMap(mesh,1,pbc_pair); - akantu::MeshUtils::computePBCMap(mesh,2,pbc_pair); - - { - std::map::iterator it = pbc_pair.begin(); - std::map::iterator end = pbc_pair.end(); - - akantu::Real * coords = mesh.getNodes().values; - akantu::UInt dim = mesh.getSpatialDimension(); - while(it != end){ - akantu::UInt i1 = (*it).first; - akantu::UInt i2 = (*it).second; - - AKANTU_DEBUG_INFO("pairing " << i1 << "(" - << coords[dim*i1] << "," << coords[dim*i1+1] << "," - << coords[dim*i1+2] - << ") with" - << i2 << "(" - << coords[dim*i2] << "," << coords[dim*i2+1] << "," - << coords[dim*i2+2] - << ")"); - ++it; - } - } - - akantu::PBCSynchronizer synch(pbc_pair); - model = new akantu::HeatTransferModel(mesh); - model->createSynchronizerRegistry(model); - model->getSynchronizerRegistry().registerSynchronizer(synch,akantu::_gst_htm_capacity); - model->getSynchronizerRegistry().registerSynchronizer(synch,akantu::_gst_htm_temperature); - /* -------------------------------------------------------------------------- */ model->readMaterials("material.dat"); model->initModel(); model->initVectors(); model->getHeatFlux().clear(); model->getCapacityLumped().clear(); model->getTemperatureGradient(type).clear(); /* -------------------------------------------------------------------------- */ - model->changeLocalEquationNumberforPBC(pbc_pair,1); + model->initPBC(1,1,1); model->assembleCapacityLumped(type); - model->getSynchronizerRegistry().synchronize(akantu::_gst_htm_capacity); /* -------------------------------------------------------------------------- */ nb_nodes = model->getFEM().getMesh().getNbNodes(); nb_element = model->getFEM().getMesh().getNbElement(type); nb_nodes = model->getFEM().getMesh().getNbNodes(); /* ------------------------------------------------------------------------ */ //get stable time step akantu::Real time_step = model->getStableTimeStep()*0.8; cout<<"time step is:"<setTimeStep(time_step); /* -------------------------------------------------------------------------- */ /// boundary conditions const akantu::Vector & nodes = model->getFEM().getMesh().getNodes(); akantu::Vector & boundary = model->getBoundary(); akantu::Vector & temperature = model->getTemperature(); akantu::Vector & heat_flux = model->getHeatFlux(); akantu::Real eps = 1e-15; double t1, t2, length; t1 = 300.; t2 = 100.; length = 1.; for (akantu::UInt i = 0; i < nb_nodes; ++i) { temperature(i) = 100.; akantu::Real dx = nodes(i,0) - length/4.; akantu::Real dy = nodes(i,1) - length/4.; akantu::Real dz = nodes(i,2) - length/4.; akantu::Real d = sqrt(dx*dx + dy*dy + dz*dz); if(d < 0.1){ boundary(i) = true; temperature(i) = 300.; } } /* -------------------------------------------------------------------------- */ DumperParaview dumper; paraviewInit(model,dumper); /* ------------------------------------------------------------------------ */ // //for testing int max_steps = 3000; /* ------------------------------------------------------------------------ */ for(int i=0; iupdateHeatFlux(); model->updateTemperature(); - model->getSynchronizerRegistry().synchronize(akantu::_gst_htm_temperature); if(i % 100 == 0) paraviewDump(dumper); if(i % 10 == 0) std::cout << "Step " << i << "/" << max_steps << std::endl; } cout<< "\n\n Stable Time Step is : " << time_step << "\n \n" <getFEM().getMesh().getNbNodes(); akantu::UInt nb_element = model->getFEM().getMesh().getNbElement(type); dumper.SetMode(TEXT); dumper.SetPoints(model->getFEM().getMesh().getNodes().values, spatial_dimension, nb_nodes, "coordinates2"); dumper.SetConnectivity((int *)model->getFEM().getMesh().getConnectivity(type).values, paraview_type, nb_element, C_MODE); dumper.AddNodeDataField(model->getTemperature().values, 1, "temperature"); dumper.AddNodeDataField(model->getHeatFlux().values, 1, "heat_flux"); dumper.AddNodeDataField(model->getCapacityLumped().values, 1, "capacity_lumped"); // dumper.AddElemDataField(model->getTemperatureGradient(type).values, // spatial_dimension, "temperature_gradient"); // dumper.AddElemDataField(model->getbtkgt().values, // 4, "btkgt"); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); } /* -------------------------------------------------------------------------- */ void paraviewDump(Dumper & dumper) { dumper.Dump(); } /* -------------------------------------------------------------------------- */ diff --git a/test/test_model/test_heat_transfer_model/test_heat_transfer_model_square2d_pbc.cc b/test/test_model/test_heat_transfer_model/test_heat_transfer_model_square2d_pbc.cc index 74d28b194..df461b8db 100644 --- a/test/test_model/test_heat_transfer_model/test_heat_transfer_model_square2d_pbc.cc +++ b/test/test_model/test_heat_transfer_model/test_heat_transfer_model_square2d_pbc.cc @@ -1,212 +1,169 @@ /** * @file test_heat_transfer_model_cube3d.cc * @author Rui WANG * @date Tue May 17 11:31:22 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; /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER #include "io_helper.h" #endif //AKANTU_USE_IOHELPER void paraviewInit(akantu::HeatTransferModel * model,Dumper & dumper); void paraviewDump(Dumper & dumper); akantu::UInt spatial_dimension = 2; akantu:: ElementType type = akantu::_triangle_3; akantu::UInt paraview_type = TRIANGLE1; /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { akantu::initialize(&argc,&argv); akantu::Mesh mesh(spatial_dimension); akantu::MeshIOMSH mesh_io; mesh_io.read("square_tri3.msh", mesh); akantu::HeatTransferModel * model; akantu::UInt nb_nodes; akantu::UInt nb_element; - mesh.computeBoundingBox(); - std::map pbc_pair; - akantu::MeshUtils::computePBCMap(mesh,0,pbc_pair); - akantu::MeshUtils::computePBCMap(mesh,1,pbc_pair); - akantu::MeshUtils::computePBCMap(mesh,2,pbc_pair); - - { - std::map::iterator it = pbc_pair.begin(); - std::map::iterator end = pbc_pair.end(); - - akantu::Real * coords = mesh.getNodes().values; - akantu::UInt dim = mesh.getSpatialDimension(); - while(it != end){ - akantu::UInt i1 = (*it).first; - akantu::UInt i2 = (*it).second; - - if (dim == 3) - AKANTU_DEBUG_INFO("pairing " << i1 << "(" - << coords[dim*i1] << "," << coords[dim*i1+1] << "," - << coords[dim*i1+2] - << ") with" - << i2 << "(" - << coords[dim*i2] << "," << coords[dim*i2+1] << "," - << coords[dim*i2+2] - << ")"); - else if (dim == 2) - AKANTU_DEBUG_INFO("pairing " << i1 << "(" - << coords[dim*i1] << "," << coords[dim*i1+1] - << ") with" - << i2 << "(" - << coords[dim*i2] << "," << coords[dim*i2+1] - << ")"); - ++it; - } - } - - akantu::PBCSynchronizer synch(pbc_pair); - model = new akantu::HeatTransferModel(mesh); - model->createSynchronizerRegistry(model); - model->getSynchronizerRegistry().registerSynchronizer(synch,akantu::_gst_htm_capacity); - model->getSynchronizerRegistry().registerSynchronizer(synch,akantu::_gst_htm_temperature); /* -------------------------------------------------------------------------- */ model->readMaterials("material.dat"); model->initModel(); model->initVectors(); model->getHeatFlux().clear(); model->getCapacityLumped().clear(); model->getTemperatureGradient(type).clear(); /* -------------------------------------------------------------------------- */ - model->changeLocalEquationNumberforPBC(pbc_pair,1); + model->initPBC(1,1,1); model->assembleCapacityLumped(type); - model->getSynchronizerRegistry().synchronize(akantu::_gst_htm_capacity); /* -------------------------------------------------------------------------- */ nb_nodes = model->getFEM().getMesh().getNbNodes(); nb_element = model->getFEM().getMesh().getNbElement(type); nb_nodes = model->getFEM().getMesh().getNbNodes(); /* ------------------------------------------------------------------------ */ //get stable time step akantu::Real time_step = model->getStableTimeStep()*0.8; cout<<"time step is:"<setTimeStep(time_step); /* -------------------------------------------------------------------------- */ /// boundary conditions const akantu::Vector & nodes = model->getFEM().getMesh().getNodes(); akantu::Vector & boundary = model->getBoundary(); akantu::Vector & temperature = model->getTemperature(); akantu::Vector & heat_flux = model->getHeatFlux(); akantu::Real eps = 1e-15; double t1, t2, length; t1 = 300.; t2 = 100.; length = 1.; 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.; } } /* -------------------------------------------------------------------------- */ DumperParaview dumper; paraviewInit(model,dumper); /* ------------------------------------------------------------------------ */ // //for testing int max_steps = 100000; /* ------------------------------------------------------------------------ */ for(int i=0; iupdateHeatFlux(); model->updateTemperature(); - model->getSynchronizerRegistry().synchronize(akantu::_gst_htm_temperature); if(i % 100 == 0) paraviewDump(dumper); if(i % 10 == 0) std::cout << "Step " << i << "/" << max_steps << std::endl; } cout<< "\n\n Stable Time Step is : " << time_step << "\n \n" <getFEM().getMesh().getNbNodes(); akantu::UInt nb_element = model->getFEM().getMesh().getNbElement(type); dumper.SetMode(TEXT); dumper.SetPoints(model->getFEM().getMesh().getNodes().values, spatial_dimension, nb_nodes, "coordinates2"); dumper.SetConnectivity((int *)model->getFEM().getMesh().getConnectivity(type).values, paraview_type, nb_element, C_MODE); dumper.AddNodeDataField(model->getTemperature().values, 1, "temperature"); dumper.AddNodeDataField(model->getHeatFlux().values, 1, "heat_flux"); dumper.AddNodeDataField(model->getCapacityLumped().values, 1, "capacity_lumped"); dumper.AddElemDataField(model->getTemperatureGradient(type).values, spatial_dimension, "temperature_gradient"); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); } /* -------------------------------------------------------------------------- */ void paraviewDump(Dumper & dumper) { dumper.Dump(); } /* -------------------------------------------------------------------------- */