diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh index 275604dce..622a149b7 100644 --- a/src/common/aka_common.hh +++ b/src/common/aka_common.hh @@ -1,414 +1,415 @@ /** * @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 #include /* -------------------------------------------------------------------------- */ #define __BEGIN_AKANTU__ namespace akantu { #define __END_AKANTU__ } /* -------------------------------------------------------------------------- */ #include "aka_config.hh" #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; /* -------------------------------------------------------------------------- */ /* Mesh/FEM/Model types */ /* -------------------------------------------------------------------------- */ typedef UInt Surface; - +typedef std::pair SurfacePair; +typedef std::list< SurfacePair > SurfacePairList; /* -------------------------------------------------------------------------- */ /// @boost sequence of element to loop on in global tasks #define AKANTU_REGULAR_ELEMENT_TYPE \ (_segment_2) \ (_segment_3) \ (_triangle_3) \ (_triangle_6) \ (_tetrahedron_4) \ (_tetrahedron_10) \ (_quadrangle_4) \ (_quadrangle_8) \ (_hexahedron_8) \ (_point) \ (_bernoulli_beam_2) #define AKANTU_COHESIVE_ELEMENT_TYPE \ (_cohesive_2d_4) \ (_cohesive_2d_6) /// @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 _bernoulli_beam_2, /// bernoulli beam 2D _cohesive_2d_4, /// first order 2D cohesive _cohesive_2d_6, /// second order 2D cohesive _max_element_type }; /// @enum MaterialType different materials implemented enum MaterialType { _elastic = 0, _max_material_type }; /// myfunction(double * position, double * stress/force, double * normal, unsigned int material_id) typedef void (*BoundaryFunction)(double *,double *, double*, unsigned int); /// @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_smm_init_mat, /// synchronization of the data to initialize materials /// HeatTransfer tags _gst_htm_capacity, /// synchronization of the nodal heat capacity _gst_htm_temperature, /// synchronization of the nodal temperature _gst_htm_gradient_temperature, /// synchronization of the element gradient temperature /// Material non local _gst_mnl_damage, /// synchronization of data to average in non local material /// 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 " " #if !defined(__aka_inline__) # define __aka_inline__ #else # define AKANTU_INCLUDE_INLINE_IMPL #endif /* -------------------------------------------------------------------------- */ #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_CONST(name, variable, type) \ inline const Vector & \ get##name (const ::akantu::ElementType & el_type, \ const GhostType & ghost_type = _not_ghost) const { \ AKANTU_DEBUG_IN(); \ \ AKANTU_DEBUG_OUT(); \ return variable(el_type, ghost_type); \ } #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE(name, variable, type) \ inline Vector & \ get##name (const ::akantu::ElementType & el_type, \ const GhostType & ghost_type = _not_ghost) { \ AKANTU_DEBUG_IN(); \ \ AKANTU_DEBUG_OUT(); \ return variable(el_type, ghost_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 _bernoulli_beam_2 : stream << "bernoulli_beam_2"; break; case _cohesive_2d_4 : stream << "cohesive_2d_4" ; break; case _cohesive_2d_6 : stream << "cohesive_2d_6" ; break; case _not_defined : stream << "undefined" ; break; case _max_element_type : stream << "ElementType(" << (int) type << ")"; break; case _point : stream << "point"; break; } return stream; } /// standard output stream operator for GhostType inline std::ostream & operator <<(std::ostream & stream, GhostType type) { switch(type) { case _not_ghost : stream << "not_ghost"; break; case _ghost : stream << "ghost" ; break; case _casper : stream << "Casper the friendly ghost"; 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" #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "aka_common_inline_impl.cc" #endif /* -------------------------------------------------------------------------- */ // BOOST PART: TOUCH ONLY IF YOU KNOW WHAT YOU ARE DOING #include #define AKANTU_BOOST_CASE_MACRO_EXCLUDE(r,macro,type) \ case type : { \ AKANTU_DEBUG_ERROR("Type (" << type << ") not handled by this function"); \ break; \ } #define AKANTU_BOOST_CASE_MACRO(r,macro,type) \ case type : { macro(type); break; } #define AKANTU_BOOST_ELEMENT_SWITCH(macro, list, exclude_list) \ do { \ switch(type) { \ BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_CASE_MACRO,macro,list) \ BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_CASE_MACRO_EXCLUDE,macro, exclude_list) \ case _not_defined: \ case _max_element_type: { \ AKANTU_DEBUG_ERROR("Wrong type : " << type); \ break; \ } \ } \ } while(0) #define AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(macro) \ AKANTU_BOOST_ELEMENT_SWITCH(macro, \ AKANTU_REGULAR_ELEMENT_TYPE, \ AKANTU_COHESIVE_ELEMENT_TYPE) #define AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(macro) \ AKANTU_BOOST_ELEMENT_SWITCH(macro, \ AKANTU_COHESIVE_ELEMENT_TYPE, \ AKANTU_REGULAR_ELEMENT_TYPE) #define AKANTU_BOOST_LIST_MACRO(r,macro,type) \ macro(type) #define AKANTU_BOOST_ELEMENT_LIST(macro,list) \ BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_LIST_MACRO,macro,list) #define AKANTU_BOOST_REGULAR_ELEMENT_LIST(macro) \ AKANTU_BOOST_ELEMENT_LIST(macro, AKANTU_REGULAR_ELEMENT_TYPE) #endif /* __AKANTU_COMMON_HH__ */ diff --git a/src/mesh_utils/mesh_pbc.cc b/src/mesh_utils/mesh_pbc.cc index 7f67fe264..33a4910c5 100644 --- a/src/mesh_utils/mesh_pbc.cc +++ b/src/mesh_utils/mesh_pbc.cc @@ -1,368 +1,368 @@ /** * @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 /* -------------------------------------------------------------------------- */ #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ __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 diff_x = p1_x - p2_x; if (dim == 2 || fabs(diff_x) > Math::getTolerance()) return diff_x > 0.0 ? false : true; 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; } return 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.lower_bounds[dir]); AKANTU_DEBUG_INFO("max " << mymesh.upper_bounds[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.lower_bounds[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.upper_bounds[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); // match pairs MeshUtils::matchPBCPairs(mymesh, dir, selected_left, selected_right, pbc_pair); } /* -------------------------------------------------------------------------- */ void MeshUtils::computePBCMap(const Mesh & mymesh, - const std::pair & surface_pair, + const SurfacePair & surface_pair, const ElementType type, std::map & pbc_pair) { std::vector selected_first; std::vector selected_second; // get access to surface information UInt nb_surface_element = mymesh.getNbElement(type); const Vector & surface_id = mymesh.getSurfaceID(type); const Vector & connect = mymesh.getConnectivity(type); UInt nodes_per_surface_element = mymesh.getNbNodesPerElement(type); // find nodes on surfaces for(UInt i=0; i < nb_surface_element; ++i) { if (surface_id(i) == surface_pair.first) { for(UInt j=0; j::iterator it_s; std::sort(selected_first.begin(), selected_first.end()); it_s = std::unique(selected_first.begin(), selected_first.end()); selected_first.resize(it_s - selected_first.begin()); std::sort(selected_second.begin(), selected_second.end()); it_s = std::unique(selected_second.begin(), selected_second.end()); selected_second.resize(it_s - selected_second.begin()); // coordinates const Vector & coords = mymesh.getNodes(); const UInt dim = mymesh.getSpatialDimension(); // variables to find min and max of surfaces Real first_max[3], first_min[3]; Real second_max[3], second_min[3]; for (UInt i=0; i::max(); second_min[i] = std::numeric_limits::max(); first_max[i] = -std::numeric_limits::max(); second_max[i] = -std::numeric_limits::max(); } // find min and max of surface nodes for (std::vector::iterator it = selected_first.begin(); it != selected_first.end(); ++it) { for (UInt i=0; i coords(*it,i)) first_min[i] = coords(*it,i); if (first_max[i] < coords(*it,i)) first_max[i] = coords(*it,i); } } for (std::vector::iterator it = selected_second.begin(); it != selected_second.end(); ++it) { for (UInt i=0; i coords(*it,i)) second_min[i] = coords(*it,i); if (second_max[i] < coords(*it,i)) second_max[i] = coords(*it,i); } } // find direction of pbc Int first_dir=-1, second_dir=-2; for (UInt i=0; i 0) ++it_right; else if (fabs(dy) < Math::getTolerance() && 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/src/model/heat_transfer/heat_transfer_model.cc b/src/model/heat_transfer/heat_transfer_model.cc index 399d4040b..1e83f7112 100644 --- a/src/model/heat_transfer/heat_transfer_model.cc +++ b/src/model/heat_transfer/heat_transfer_model.cc @@ -1,619 +1,619 @@ /** * @file heat_transfer_model.cc * @author Rui WANG * @author Guillaume Anciaux * @author Nicolas Richart * @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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "heat_transfer_model.hh" #include "aka_math.hh" #include "aka_common.hh" #include "fem_template.hh" #include "mesh.hh" #include "static_communicator.hh" #include "parser.hh" #include "generalized_trapezoidal.hh" // #include "sparse_matrix.hh" // #include "solver.hh" // #ifdef AKANTU_USE_MUMPS // #include "solver_mumps.hh" // #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ HeatTransferModel::HeatTransferModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : Model(id, memory_id), integrator(new ForwardEuler()), spatial_dimension(dim), temperature_gradient("temperature_gradient", id) { AKANTU_DEBUG_IN(); createSynchronizerRegistry(this); if (spatial_dimension == 0) spatial_dimension = mesh.getSpatialDimension(); std::stringstream sstr; sstr << id << ":fem"; registerFEMObject(sstr.str(), mesh,spatial_dimension); this->temperature= NULL; this->residual = NULL; this->boundary = NULL; this->conductivity_variation = 0.0; this->t_ref = 0; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::initModel() { getFEM().initShapeFunctions(_not_ghost); getFEM().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::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_htm_capacity); synch_registry->registerSynchronizer(synch_parallel, _gst_htm_temperature); synch_registry->registerSynchronizer(synch_parallel, _gst_htm_gradient_temperature); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void HeatTransferModel::initPBC(UInt x, UInt y, UInt z) { +void HeatTransferModel::initPBC() { AKANTU_DEBUG_IN(); - - Model::initPBC(x,y,z); + + Model::initPBC(); PBCSynchronizer * synch = new PBCSynchronizer(pbc_pair); - + synch_registry->registerSynchronizer(*synch, _gst_htm_capacity); synch_registry->registerSynchronizer(*synch, _gst_htm_temperature); changeLocalEquationNumberforPBC(pbc_pair,1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::initVectors() { AKANTU_DEBUG_IN(); UInt nb_nodes = getFEM().getMesh().getNbNodes(); std::stringstream sstr_temp; sstr_temp << id << ":temperature"; std::stringstream sstr_temp_rate; sstr_temp << id << ":temperature_rate"; std::stringstream sstr_inc; sstr_temp_rate<< id << ":increment"; std::stringstream sstr_residual; sstr_residual << id << ":residual"; std::stringstream sstr_lump; sstr_lump << id << ":lumped"; std::stringstream sstr_boun; sstr_boun << id << ":boundary"; temperature = &(alloc(sstr_temp.str(), nb_nodes, 1, REAL_INIT_VALUE)); temperature_rate = &(alloc(sstr_temp_rate.str(), nb_nodes, 1, REAL_INIT_VALUE)); increment = &(alloc(sstr_inc.str(), nb_nodes, 1, REAL_INIT_VALUE)); residual = &(alloc(sstr_residual.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)); Mesh::ConnectivityTypeList::const_iterator it; /* -------------------------------------------------------------------------- */ // byelementtype vectors getFEM().getMesh().initByElementTypeVector(temperature_on_qpoints, 1, spatial_dimension); getFEM().getMesh().initByElementTypeVector(temperature_gradient, spatial_dimension, spatial_dimension); getFEM().getMesh().initByElementTypeVector(conductivity_on_qpoints, spatial_dimension*spatial_dimension, spatial_dimension); getFEM().getMesh().initByElementTypeVector(k_gradt_on_qpoints, spatial_dimension, spatial_dimension); getFEM().getMesh().initByElementTypeVector(bt_k_gT, 1, spatial_dimension,true); getFEM().getMesh().initByElementTypeVector(int_bt_k_gT, 1, spatial_dimension,true); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; const Mesh::ConnectivityTypeList & type_list = getFEM().getMesh().getConnectivityTypeList(gt); for(it = type_list.begin(); it != type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; UInt nb_element = getFEM().getMesh().getNbElement(*it, gt); UInt nb_quad_points = this->getFEM().getNbQuadraturePoints(*it, gt) * nb_element; temperature_on_qpoints(*it, gt).resize(nb_quad_points); temperature_on_qpoints(*it, gt).clear(); temperature_gradient(*it, gt).resize(nb_quad_points); temperature_gradient(*it, gt).clear(); conductivity_on_qpoints(*it, gt).resize(nb_quad_points); conductivity_on_qpoints(*it, gt).clear(); k_gradt_on_qpoints(*it, gt).resize(nb_quad_points); k_gradt_on_qpoints(*it, gt).clear(); bt_k_gT(*it, gt).resize(nb_quad_points); bt_k_gT(*it, gt).clear(); int_bt_k_gT(*it, gt).resize(nb_element); bt_k_gT(*it, gt).clear(); } } /* -------------------------------------------------------------------------- */ dof_synchronizer = new DOFSynchronizer(getFEM().getMesh(),1); dof_synchronizer->initLocalDOFEquationNumbers(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ HeatTransferModel::~HeatTransferModel() { AKANTU_DEBUG_IN(); delete [] conductivity; if(dof_synchronizer) delete dof_synchronizer; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleCapacityLumped(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); FEM & fem = getFEM(); 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; UInt nb_element = getFEM().getMesh().getNbElement(*it,ghost_type); UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(*it, ghost_type); Vector rho_1 (nb_element * nb_quadrature_points,1, capacity * density); fem.assembleFieldLumped(rho_1,1,*capacity_lumped, dof_synchronizer->getLocalDOFEquationNumbers(), *it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleCapacityLumped() { AKANTU_DEBUG_IN(); capacity_lumped->clear(); assembleCapacityLumped(_not_ghost); assembleCapacityLumped(_ghost); getSynchronizerRegistry().synchronize(akantu::_gst_htm_capacity); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::updateResidual() { AKANTU_DEBUG_IN(); /// @f$ r = q_{ext} - q_{int} - C \dot T @f$ // start synchronization synch_registry->asynchronousSynchronize(_gst_htm_temperature); // finalize communications synch_registry->waitEndSynchronize(_gst_htm_temperature); //clear the array /// first @f$ r = q_{ext} @f$ residual->clear(); /// then @f$ r -= q_{int} @f$ // update the not ghost ones updateResidual(_not_ghost); // update for the received ghosts updateResidual(_ghost); /// finally @f$ r -= C \dot T @f$ // lumped C UInt nb_nodes = temperature_rate->getSize(); UInt nb_degre_of_freedom = temperature_rate->getNbComponent(); Real * capacity_val = capacity_lumped->values; Real * temp_rate_val = temperature_rate->values; Real * res_val = residual->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes * nb_degre_of_freedom; ++n) { if(!(*boundary_val)) { *res_val -= *capacity_val * *temp_rate_val; } boundary_val++; res_val++; capacity_val++; temp_rate_val++; } #ifndef AKANTU_NDEBUG getSynchronizerRegistry().synchronize(akantu::_gst_htm_gradient_temperature); #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::computeConductivityOnQuadPoints(const GhostType & ghost_type) { const Mesh::ConnectivityTypeList & type_list = this->getFEM().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; //compute the temperature on quadrature points this->getFEM().interpolateOnQuadraturePoints(*temperature, temperature_on_qpoints(*it, ghost_type), 1 ,*it,ghost_type); Vector::iterator c_iterator = conductivity_on_qpoints(*it,ghost_type).begin(spatial_dimension,spatial_dimension); Real * T_val = temperature_on_qpoints(*it, ghost_type).values; UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(*it, ghost_type); UInt nb_element = getFEM().getMesh().getNbElement(*it,ghost_type); for (UInt el = 0; el < nb_element; ++el) { for (UInt q = 0; q < nb_quadrature_points; ++q) { Real * c_onquad_val = c_iterator->storage(); for (UInt i = 0; i < spatial_dimension*spatial_dimension; ++i) { // if (i == 0 && *T_val - t_ref > 0) // AKANTU_DEBUG_INFO("conductivity is " << conductivity[i] // << " variation is " << conductivity_variation // << " temp increase is " << *T_val - t_ref // << " new conductivity is " << conductivity[i] + conductivity_variation* (*T_val - t_ref)); c_onquad_val[i] = conductivity[i] + conductivity_variation* (*T_val - t_ref); } ++T_val; ++c_iterator; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::computeKgradT(const GhostType & ghost_type) { computeConductivityOnQuadPoints(ghost_type); const Mesh::ConnectivityTypeList & type_list = this->getFEM().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; UInt nb_element = getFEM().getMesh().getNbElement(*it,ghost_type); UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(*it, ghost_type); Vector & gradient = temperature_gradient(*it, ghost_type); gradient.resize(nb_element * nb_quadrature_points); this->getFEM().gradientOnQuadraturePoints(*temperature, gradient, 1 ,*it,ghost_type); Vector::iterator c_iterator = conductivity_on_qpoints(*it,ghost_type).begin(spatial_dimension,spatial_dimension); Vector::iterator gT_iterator = temperature_gradient(*it,ghost_type).begin(spatial_dimension); Vector::iterator k_gT_iterator = k_gradt_on_qpoints(*it,ghost_type).begin(spatial_dimension); for (UInt el = 0; el < nb_element; ++el) { for (UInt i = 0; i < nb_quadrature_points; ++i){ Math::matrixt_vector(spatial_dimension, spatial_dimension, c_iterator->storage(), gT_iterator->storage(), k_gT_iterator->storage()); ++c_iterator; ++gT_iterator; ++k_gT_iterator; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::updateResidual(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); const Mesh::ConnectivityTypeList & type_list = this->getFEM().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; Vector & shapes_derivatives = const_cast &>(getFEM().getShapesDerivatives(*it,ghost_type)); UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(*it, ghost_type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); UInt nb_element = getFEM().getMesh().getNbElement(*it,ghost_type); // compute k \grad T computeKgradT(ghost_type); Vector::iterator k_gT_iterator = k_gradt_on_qpoints(*it,ghost_type).begin(spatial_dimension); Vector::iterator shapesd_iterator = shapes_derivatives.begin(nb_nodes_per_element,spatial_dimension); Vector::iterator bt_k_gT_iterator = bt_k_gT(*it,ghost_type).begin(nb_nodes_per_element); // UInt bt_k_gT_size = nb_nodes_per_element; // 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::matrix_vector(nb_nodes_per_element, spatial_dimension, shapesd_iterator->storage(), k_gT_iterator->storage(), bt_k_gT_iterator->storage()); ++shapesd_iterator; ++k_gT_iterator; ++bt_k_gT_iterator; } } this->getFEM().integrate(bt_k_gT(*it,ghost_type), int_bt_k_gT(*it,ghost_type), nb_nodes_per_element, *it,ghost_type); this->getFEM().assembleVector(int_bt_k_gT(*it,ghost_type), *residual, dof_synchronizer->getLocalDOFEquationNumbers(), 1, *it,ghost_type,NULL,-1); //delete q_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::solveExplicitLumped() { AKANTU_DEBUG_IN(); UInt nb_nodes = increment->getSize(); UInt nb_degre_of_freedom = increment->getNbComponent(); Real * capa_val = capacity_lumped->values; Real * res_val = residual->values; bool * boundary_val = boundary->values; Real * inc = increment->values; for (UInt n = 0; n < nb_nodes * nb_degre_of_freedom; ++n) { if(!(*boundary_val)) { *inc = (*res_val / *capa_val); } res_val++; boundary_val++; inc++; capa_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::explicitPred() { AKANTU_DEBUG_IN(); integrator->integrationSchemePred(time_step, *temperature, *temperature_rate, *boundary); UInt nb_nodes = temperature->getSize(); UInt nb_degre_of_freedom = temperature->getNbComponent(); Real * temp = temperature->values; for (UInt n = 0; n < nb_nodes * nb_degre_of_freedom; ++n, ++temp) if(*temp < 0.) *temp = 0.; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::explicitCorr() { AKANTU_DEBUG_IN(); integrator->integrationSchemeCorrTempRate(time_step, *temperature, *temperature_rate, *boundary, *increment); 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, _not_ghost).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); } AKANTU_DEBUG_INFO("The minimum element size : " << min_el_size << " and the max conductivity is : " << conductivitymax); delete [] u; } Real min_dt = 2 * min_el_size * min_el_size * density * capacity/conductivitymax; StaticCommunicator::getStaticCommunicator()->allReduce(&min_dt, 1, _so_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::readMaterials(const std::string & filename) { Parser parser; parser.open(filename); std::string mat_type = parser.getNextSection("heat"); if (mat_type != ""){ parser.readSection(*this); } else AKANTU_DEBUG_ERROR("did not find any section with material info " <<"for heat conduction"); } /* -------------------------------------------------------------------------- */ bool 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(UInt i = 0; i < 3; i++) for(UInt j = 0; j < 3; 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 == "conductivity_variation") { str >> conductivity_variation; } else if (key == "temperature_reference") { str >> t_ref; } 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); } else { return false; } return true; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/heat_transfer/heat_transfer_model.hh b/src/model/heat_transfer/heat_transfer_model.hh index 4f7e1ff46..75a15cef3 100644 --- a/src/model/heat_transfer/heat_transfer_model.hh +++ b/src/model/heat_transfer/heat_transfer_model.hh @@ -1,289 +1,289 @@ /** * @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 "aka_common.hh" #include "aka_memory.hh" #include "model.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" namespace akantu { class IntegrationScheme1stOrder; // class Solver; // class SparseMatrix; } __BEGIN_AKANTU__ class HeatTransferModel : public Model, public DataAccessor { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef FEMTemplate MyFEMType; HeatTransferModel(UInt spatial_dimension, const ID & id = "heat_transfer_model", const MemoryID & memory_id = 0) ; HeatTransferModel(Mesh & mesh, UInt spatial_dimension = 0, const ID & id = "heat_transfer_model", const MemoryID & memory_id = 0); virtual ~HeatTransferModel() ; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// set the parameters bool setParam(const std::string & key, const std::string & value); /// read one material file to instantiate all the materials void readMaterials(const std::string & filename); /// allocate all vectors void initVectors(); /// register the tags associated with the parallel synchronizer void initParallel(MeshPartition * partition, DataAccessor * data_accessor=NULL); /// initialize the model void initModel(); /// init PBC synchronizer - void initPBC(UInt x,UInt y, UInt z); + void initPBC(); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const {}; /* ------------------------------------------------------------------------ */ /* Methods for explicit */ /* ------------------------------------------------------------------------ */ public: /// compute and get the stable time step Real getStableTimeStep(); /// compute the heat flux void updateResidual(); /// solve the system in temperature rate @f$C\delta \dot T = q_{n+1} - C \dot T_{n}@f$ void solveExplicitLumped(); /// calculate the lumped capacity vector for heat transfer problem void assembleCapacityLumped(); /// update the temperature from the temperature rate void explicitPred(); /// update the temperature rate from the increment void explicitCorr(); // /// initialize the heat flux // void initializeResidual(Vector &temp); // /// initialize temperature // void initializeTemperature(Vector &temp); private: /// compute the heat flux on ghost types void updateResidual(const GhostType & ghost_type); /// calculate the lumped capacity vector for heat transfer problem (w ghosttype) void assembleCapacityLumped(const GhostType & ghost_type); /// compute the conductivity tensor for each quadrature point in an array void computeConductivityOnQuadPoints(const GhostType & ghost_type); /// compute vector k \grad T for each quadrature point void computeKgradT(const GhostType & ghost_type); /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: __aka_inline__ UInt getNbDataToPack(const Element & element, SynchronizationTag tag) const; __aka_inline__ UInt getNbDataToUnpack(const Element & element, SynchronizationTag tag) const; __aka_inline__ void packData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag) const; __aka_inline__ void unpackData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag); __aka_inline__ UInt getNbDataToPack(SynchronizationTag tag) const; __aka_inline__ UInt getNbDataToUnpack(SynchronizationTag tag) const; __aka_inline__ void packData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag) const; __aka_inline__ void unpackData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the current value of the time step AKANTU_GET_MACRO(TimeStep, time_step, Real); /// set the value of the time step AKANTU_SET_MACRO(TimeStep, time_step, Real); /// get the assembled heat flux AKANTU_GET_MACRO(Residual, *residual, Vector&); /// get the lumped capacity AKANTU_GET_MACRO(CapacityLumped, * capacity_lumped, Vector&); /// get the boundary vector AKANTU_GET_MACRO(Boundary, * boundary, Vector&); /// get the temperature gradient AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(TemperatureGradient, temperature_gradient, Real); /// get the conductivity on q points AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ConductivityOnQpoints, conductivity_on_qpoints, Real); /// get the conductivity on q points AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(TemperatureOnQpoints, temperature_on_qpoints, Real); /// get the temperature AKANTU_GET_MACRO(Temperature, *temperature, Vector &); /// get the temperature derivative AKANTU_GET_MACRO(TemperatureRate, *temperature_rate, Vector &); /// get the equation number Vector AKANTU_GET_MACRO(EquationNumber, *equation_number, const Vector &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: IntegrationScheme1stOrder * integrator; /// time step Real time_step; /// temperatures array Vector * temperature; /// temperatures derivatives array Vector * temperature_rate; /// increment array (@f$\delta \dot T@f$ or @f$\delta T@f$) Vector * increment; /// the spatial dimension UInt spatial_dimension; /// the density Real density; /// the speed of the changing temperature ByElementTypeReal temperature_gradient; /// temperature field on quadrature points ByElementTypeReal temperature_on_qpoints; /// conductivity tensor on quadrature points ByElementTypeReal conductivity_on_qpoints; /// vector k \grad T on quad points ByElementTypeReal k_gradt_on_qpoints; /// vector \int \grad N k \grad T ByElementTypeReal int_bt_k_gT; /// vector \grad N k \grad T ByElementTypeReal bt_k_gT; //external flux vector Vector * external_flux; /// residuals array Vector * residual; /// position of a dof in the K matrix Vector * equation_number; //lumped vector Vector * capacity_lumped; /// boundary vector Vector * boundary; //realtime Real time; ///capacity Real capacity; //conductivity matrix Real* conductivity; //linear variation of the conductivity (for temperature dependent conductivity) Real conductivity_variation; // reference temperature for the interpretation of temperature variation Real t_ref; //the biggest parameter of conductivity matrix Real conductivitymax; }; /* -------------------------------------------------------------------------- */ /* __aka_inline__ functions */ /* -------------------------------------------------------------------------- */ #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "heat_transfer_model_inline_impl.cc" #endif /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const HeatTransferModel & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_HEAT_TRANSFER_MODEL_HH__ */ diff --git a/src/model/model.cc b/src/model/model.cc index d3aab5394..8af2d3478 100644 --- a/src/model/model.cc +++ b/src/model/model.cc @@ -1,145 +1,141 @@ /** * @file model.cc * @author Nicolas Richart * @date Fri Sep 16 14:46:01 2011 * * @brief * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "model.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ Model::Model(const ID & id, const MemoryID & memory_id) : Memory(memory_id), id(id),synch_registry(NULL),is_pbc_slave_node(0,1,"is_pbc_slave_node") { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Model::createSynchronizerRegistry(DataAccessor * data_accessor){ synch_registry = new SynchronizerRegistry(*data_accessor); } /* -------------------------------------------------------------------------- */ -void Model::initPBC(UInt x, UInt y, UInt z){ +void Model::setPBC(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); - - initPBC(); } /* -------------------------------------------------------------------------- */ -void Model::initPBC(std::list< std::pair > & surface_pairs, - ElementType surface_e_type){ +void Model::setPBC(SurfacePairList & surface_pairs, + ElementType surface_e_type){ Mesh & mesh = getFEM().getMesh(); - std::list< std::pair >::iterator s_it; + SurfacePairList::iterator s_it; for(s_it = surface_pairs.begin(); s_it != surface_pairs.end(); ++s_it) { MeshUtils::computePBCMap(mesh, *s_it, surface_e_type, pbc_pair); } - - initPBC(); } /* -------------------------------------------------------------------------- */ void Model::initPBC() { Mesh & mesh = getFEM().getMesh(); std::map::iterator it = pbc_pair.begin(); std::map::iterator end = pbc_pair.end(); is_pbc_slave_node.resize(mesh.getNbNodes()); Real * coords = mesh.getNodes().values; UInt dim = mesh.getSpatialDimension(); while(it != end){ UInt i1 = (*it).first; UInt i2 = (*it).second; is_pbc_slave_node(i1) = true; 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; } } /* -------------------------------------------------------------------------- */ Synchronizer & Model::createParallelSynch(MeshPartition * partition, __attribute__((unused)) 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; } /* -------------------------------------------------------------------------- */ 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(); } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/model.hh b/src/model/model.hh index bf789c690..ffbaf8f21 100644 --- a/src/model/model.hh +++ b/src/model/model.hh @@ -1,179 +1,178 @@ /** * @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: Model(const ID & id = "model", const MemoryID & memory_id = 0); virtual ~Model(); typedef std::map FEMMap; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual void initModel() = 0; /// create the synchronizer registry object void createSynchronizerRegistry(DataAccessor * data_accessor); /// create a parallel synchronizer and distribute the mesh Synchronizer & createParallelSynch(MeshPartition * partition, DataAccessor * data_accessor); /// change local equation number so that PBC is assembled properly void changeLocalEquationNumberforPBC(std::map & pbc_pair,UInt dimension); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const = 0; /// initialize the model for PBC - virtual void initPBC(UInt x, UInt y, UInt z); - virtual void initPBC(std::list< std::pair > & surface_pairs, - ElementType surface_e_type); + void setPBC(UInt x, UInt y, UInt z); + void setPBC(SurfacePairList & surface_pairs, + ElementType surface_e_type); -private: - void initPBC(); + virtual void initPBC(); /* ------------------------------------------------------------------------ */ /* 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 __aka_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 __aka_inline__ void registerFEMObject(const std::string & name, Mesh & mesh, UInt spatial_dimension); /// unregister a fem object associated with name __aka_inline__ void unRegisterFEMObject(const std::string & name); /// return the synchronizer registry SynchronizerRegistry & getSynchronizerRegistry(); protected: /// return the fem object associated with a provided name template __aka_inline__ FEMClass & getFEMClass(std::string name = "") const; /// return the fem boundary object associated with a provided name template __aka_inline__ FEMClass & getFEMClassBoundary(std::string name = ""); /// returns if node is slave in pbc __aka_inline__ bool getIsPBCSlaveNode(const UInt node); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id ID id; /// the main fem object present in all models FEMMap fems; /// the fem object present in all models for boundaries FEMMap fems_boundary; /// default fem object std::string default_fem; /// synchronizer registry SynchronizerRegistry * synch_registry; /// handle the equation number things DOFSynchronizer * dof_synchronizer; /// pbc pairs std::map pbc_pair; /// flag per node to know is pbc slave Vector is_pbc_slave_node; }; /* -------------------------------------------------------------------------- */ /* __aka_inline__ functions */ /* -------------------------------------------------------------------------- */ #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "model_inline_impl.cc" #endif /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Model & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_MODEL_HH__ */ diff --git a/src/model/model_inline_impl.cc b/src/model/model_inline_impl.cc index e65bd50d0..d9fa58211 100644 --- a/src/model/model_inline_impl.cc +++ b/src/model/model_inline_impl.cc @@ -1,167 +1,169 @@ /** * @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 SynchronizerRegistry & Model::getSynchronizerRegistry(){ AKANTU_DEBUG_ASSERT(synch_registry,"synchronizer registry not initialized:" << " did you call createSynchronizerRegistry ?"); return *synch_registry; } /* -------------------------------------------------------------------------- */ template inline FEMClass & Model::getFEMClassBoundary(std::string name){ AKANTU_DEBUG_IN(); if (name == "") name = default_fem; FEMMap::const_iterator it_boun = fems_boundary.find(name); FEMClass * tmp_fem_boundary; if (it_boun == fems_boundary.end()){ AKANTU_DEBUG_INFO("Creating FEM boundary " << name); FEMMap::const_iterator it = fems.find(name); AKANTU_DEBUG_ASSERT(it != fems.end(), "The FEM " << name << " is not registered"); UInt spatial_dimension = it->second->getElementDimension(); std::stringstream sstr; sstr << id << ":fem_boundary:" << name; tmp_fem_boundary = new FEMClass(it->second->getMesh(), spatial_dimension-1, sstr.str(), memory_id); fems_boundary[name] = tmp_fem_boundary; } else { tmp_fem_boundary = dynamic_cast(it_boun->second); } AKANTU_DEBUG_OUT(); return *tmp_fem_boundary; } /* -------------------------------------------------------------------------- */ template inline FEMClass & Model::getFEMClass(std::string name) const{ AKANTU_DEBUG_IN(); if (name == "") name = default_fem; FEMMap::const_iterator it = fems.find(name); AKANTU_DEBUG_ASSERT(it != fems.end(), "The FEM " << name << " is not registered"); AKANTU_DEBUG_OUT(); return dynamic_cast(*(it->second)); } /* -------------------------------------------------------------------------- */ inline void Model::unRegisterFEMObject(const std::string & name){ FEMMap::iterator it = fems.find(name); AKANTU_DEBUG_ASSERT(it != fems.end(), "FEM object with name " << name << " was not found"); delete((*it).second); fems.erase(it); if (!fems.empty()) default_fem = (*fems.begin()).first; } /* -------------------------------------------------------------------------- */ template inline void Model::registerFEMObject(const std::string & name, Mesh & mesh, UInt spatial_dimension){ if (fems.size() == 0) default_fem = name; 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; + (*dof_synchronizer->getGlobalDOFEquationNumbersPointer()) + (node_slave*dimension+i) = dimension*node_master+i; } } } /* -------------------------------------------------------------------------- */ inline bool Model::getIsPBCSlaveNode(const UInt node) { // if no pbc is defined, is_pbc_slave_node is of size zero if (is_pbc_slave_node.getSize() == 0) return false; else return is_pbc_slave_node(node); } diff --git a/src/model/solid_mechanics/contact/contact_search_explicit.cc b/src/model/solid_mechanics/contact/contact_search_explicit.cc index 50e1e9f8b..49688b9f3 100644 --- a/src/model/solid_mechanics/contact/contact_search_explicit.cc +++ b/src/model/solid_mechanics/contact/contact_search_explicit.cc @@ -1,635 +1,639 @@ /** * @file contact_search_explicit.cc * @author David Kammer * @date Tue Oct 26 18:49:04 2010 * * @brief Specialization of the contact search structure for 3D within an * explicit time scheme * * @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 "regular_grid_neighbor_structure.hh" #include "contact_search_explicit.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ ContactSearchExplicit::ContactSearchExplicit(Contact & contact, const ContactNeighborStructureType & neighbors_structure_type, const ContactSearchType & type, const ContactSearchID & id) : ContactSearch(contact, neighbors_structure_type, type, id), spatial_dimension(contact.getModel().getSpatialDimension()), mesh(contact.getModel().getFEM().getMesh()) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactSearchExplicit::findPenetration(const Surface & master_surface, PenetrationList & penetration_list) { AKANTU_DEBUG_IN(); /// get the NodesNeighborList for the given master surface std::map::iterator it_surface; for (it_surface = neighbors_structure.begin(); it_surface != neighbors_structure.end(); ++it_surface) { if(it_surface->first == master_surface) { break; } } AKANTU_DEBUG_ASSERT(it_surface != neighbors_structure.end(), "Master surface not found in this search object " << id); const NodesNeighborList & neighbor_list = dynamic_cast(it_surface->second->getNeighborList()); UInt nb_impactor_nodes = neighbor_list.impactor_nodes.getSize(); Vector * closest_master_nodes = new Vector(nb_impactor_nodes, 1); Vector * has_closest_master_node = new Vector(nb_impactor_nodes, 1, false); findClosestMasterNodes(master_surface, closest_master_nodes, has_closest_master_node); UInt * closest_master_nodes_val = closest_master_nodes->values; bool * has_closest_master_node_val = has_closest_master_node->values; /// get list of impactor and master nodes from neighbor list UInt * impactor_nodes_val = neighbor_list.impactor_nodes.values; //const Mesh & mesh = contact.getModel().getFEM().getMesh(); const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; /// find existing surface element types UInt nb_facet_types = 0; ElementType facet_type[_max_element_type]; for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; if(mesh.getSpatialDimension(type) == spatial_dimension) { ElementType current_facet_type = mesh.getFacetElementType(type); facet_type[nb_facet_types++] = current_facet_type; /// initialization of penetration list // std::stringstream sstr_facets_offset; // sstr_facets_offset << id << ":penetrated_facets_offset:" << current_facet_type; penetration_list.penetrated_facets_offset.alloc(0, 1, current_facet_type, _not_ghost); // std::stringstream sstr_facets; // sstr_facets << id << ":penetrated_facets:" << current_facet_type; penetration_list.penetrated_facets.alloc(0, 1, current_facet_type, _not_ghost); // std::stringstream sstr_normals; // sstr_normals << id << ":facets_normals:" << current_facet_type; penetration_list.facets_normals.alloc(0, spatial_dimension, current_facet_type, _not_ghost); // std::stringstream sstr_gaps; // sstr_gaps << id << ":gaps:" << current_facet_type; penetration_list.gaps.alloc(0, 1, current_facet_type, _not_ghost); // std::stringstream sstr_projected_positions; // sstr_projected_positions << id << ":projected_positions:" << current_facet_type; penetration_list.projected_positions.alloc(0, spatial_dimension, current_facet_type, _not_ghost); } } for(UInt in = 0; in < nb_impactor_nodes; ++in) { if (!has_closest_master_node_val[in]) continue; UInt current_impactor_node = impactor_nodes_val[in]; UInt closest_master_node = closest_master_nodes_val[in]; std::vector surface_elements; Vector * are_inside = new Vector(0, 1); Vector * are_in_projection_area = new Vector(0, 1); Element considered_element; for (UInt el_type = 0; el_type < nb_facet_types; ++el_type) { ElementType type = facet_type[el_type]; UInt * surface_id_val = mesh.getSurfaceID(type, _not_ghost).values; const Vector & node_to_elements_offset = contact.getNodeToElementsOffset(type, _not_ghost); const Vector & node_to_elements = contact.getNodeToElements(type, _not_ghost); UInt * node_to_elements_offset_val = node_to_elements_offset.values; UInt * node_to_elements_val = node_to_elements.values; UInt min_element_offset = node_to_elements_offset_val[closest_master_node]; UInt max_element_offset = node_to_elements_offset_val[closest_master_node + 1]; considered_element.type = type; for(UInt el = min_element_offset; el < max_element_offset; ++el) { UInt surface_element = node_to_elements_val[el]; if(surface_id_val[surface_element] == master_surface) { bool is_inside; bool is_in_projection_area; checkPenetrationSituation(current_impactor_node, surface_element, type, is_inside, is_in_projection_area); considered_element.element = surface_element; surface_elements.push_back(considered_element); are_inside->push_back(is_inside); are_in_projection_area->push_back(is_in_projection_area); } } } UInt nb_penetrated_elements = 0; UInt nb_elements_type[_max_element_type]; memset(nb_elements_type, 0, sizeof(UInt) * _max_element_type); bool * are_inside_val = are_inside->values; bool * are_in_projection_area_val = are_in_projection_area->values; bool is_outside = false; UInt nb_surface_elements = static_cast(surface_elements.size()); // add all elements for which impactor is inside and in projection area for(UInt el = 0; el < nb_surface_elements; ++el) { if(are_inside_val[el] && are_in_projection_area_val[el]) { ElementType current_type = surface_elements.at(el).type; UInt current_element = surface_elements.at(el).element; penetration_list.penetrated_facets(current_type, _not_ghost).push_back(current_element); Real normal[3]; Real projected_position[3]; Real gap; computeComponentsOfProjection(current_impactor_node, current_element, current_type, normal, gap, projected_position); penetration_list.facets_normals(current_type, _not_ghost).push_back(normal); penetration_list.projected_positions(current_type, _not_ghost).push_back(projected_position); penetration_list.gaps(current_type, _not_ghost).push_back(gap); nb_penetrated_elements++; nb_elements_type[current_type]++; } } if(nb_penetrated_elements > 0) { for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; if(mesh.getSpatialDimension(type) == spatial_dimension) { penetration_list.penetrating_nodes.push_back(current_impactor_node); ElementType current_facet_type = mesh.getFacetElementType(type); penetration_list.penetrated_facets_offset(current_facet_type, _not_ghost).push_back(nb_elements_type[current_facet_type]); } } } // if there was no element which is inside and in projection area // check if node is not definitly outside else { for(UInt el = 0; el < nb_surface_elements && !is_outside; ++el) { if(!are_inside_val[el] && are_in_projection_area_val[el]) { is_outside = true; } } // it is not definitly outside take all elements to which it is at least inside if(!is_outside) { bool found_inside_node = false; for(UInt el = 0; el < nb_surface_elements; ++el) { if(are_inside_val[el] && !are_in_projection_area_val[el]) { ElementType current_type = surface_elements.at(el).type; UInt current_element = surface_elements.at(el).element; penetration_list.penetrated_facets(current_type, _not_ghost).push_back(current_element); Real normal[3]; Real projected_position[3]; Real gap; computeComponentsOfProjection(current_impactor_node, current_element, current_type, normal, gap, projected_position); penetration_list.facets_normals(current_type, _not_ghost).push_back(normal); penetration_list.projected_positions(current_type, _not_ghost).push_back(projected_position); penetration_list.gaps(current_type, _not_ghost).push_back(gap); nb_penetrated_elements++; nb_elements_type[current_type]++; found_inside_node = true; } } if(found_inside_node) { for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; if(mesh.getSpatialDimension(type) == spatial_dimension) { penetration_list.penetrating_nodes.push_back(current_impactor_node); ElementType current_facet_type = mesh.getFacetElementType(type); penetration_list.penetrated_facets_offset(current_facet_type, _not_ghost).push_back(nb_elements_type[current_facet_type]); } } } } } delete are_in_projection_area; delete are_inside; } // make the offset table a real offset table for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; if(mesh.getSpatialDimension(type) == spatial_dimension) { ElementType current_facet_type = mesh.getFacetElementType(type); UInt tmp_nb_facets = penetration_list.penetrated_facets_offset(current_facet_type, _not_ghost).getSize(); penetration_list.penetrated_facets_offset(current_facet_type, _not_ghost).resize(tmp_nb_facets+1); Vector & tmp_penetrated_facets_offset = (penetration_list.penetrated_facets_offset(current_facet_type, _not_ghost)); UInt * tmp_penetrated_facets_offset_val = tmp_penetrated_facets_offset.values; for (UInt i = 1; i < tmp_nb_facets; ++i) tmp_penetrated_facets_offset_val[i] += tmp_penetrated_facets_offset_val[i - 1]; for (UInt i = tmp_nb_facets; i > 0; --i) tmp_penetrated_facets_offset_val[i] = tmp_penetrated_facets_offset_val[i - 1]; tmp_penetrated_facets_offset_val[0] = 0; } } delete closest_master_nodes; delete has_closest_master_node; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactSearchExplicit::findClosestMasterNodes(const Surface & master_surface, Vector * closest_master_nodes, Vector * has_closest_master_node) { AKANTU_DEBUG_IN(); bool * has_closest_master_node_val = has_closest_master_node->values; UInt * closest_master_nodes_val = closest_master_nodes->values; /// get the NodesNeighborList for the given master surface std::map::iterator it; for (it = neighbors_structure.begin(); it != neighbors_structure.end(); ++it) { if(it->first == master_surface) { break; } } AKANTU_DEBUG_ASSERT(it != neighbors_structure.end(), "Master surface not found in this search object " << id); const NodesNeighborList & neighbor_list = dynamic_cast(it->second->getNeighborList()); /// get list of impactor and master nodes from neighbor list UInt * impactor_nodes_val = neighbor_list.impactor_nodes.values; UInt * master_nodes_offset_val = neighbor_list.master_nodes_offset.values; UInt * master_nodes_val = neighbor_list.master_nodes.values; /// loop over all impactor nodes and find for each the closest master node UInt nb_impactor_nodes = neighbor_list.impactor_nodes.getSize(); for(UInt imp = 0; imp < nb_impactor_nodes; ++imp) { UInt current_impactor_node = impactor_nodes_val[imp]; UInt min_offset = master_nodes_offset_val[imp]; UInt max_offset = master_nodes_offset_val[imp + 1]; // if there is no master node go to next impactor node if (min_offset == max_offset) continue; Real min_square_distance = std::numeric_limits::max(); UInt closest_master_node = (UInt)-1; // for finding error for(UInt mn = min_offset; mn < max_offset; ++mn) { UInt current_master_node = master_nodes_val[mn]; Real square_distance = computeSquareDistanceBetweenNodes(current_impactor_node, current_master_node); if(min_square_distance > square_distance) { min_square_distance = square_distance; closest_master_node = current_master_node; } } AKANTU_DEBUG_ASSERT(closest_master_node != ((UInt)-1), "could not find a closest master node for impactor node: " << current_impactor_node); has_closest_master_node_val[imp] = true; closest_master_nodes_val[imp] = closest_master_node; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactSearchExplicit::computeComponentsOfProjection(const UInt impactor_node, const UInt surface_element, const ElementType type, Real * normal, Real & gap, Real * projected_position) { AKANTU_DEBUG_IN(); switch(type) { case _segment_2: { computeComponentsOfProjectionSegment2(impactor_node, surface_element, normal, gap, projected_position); break; } case _triangle_3: { computeComponentsOfProjectionTriangle3(impactor_node, surface_element, normal, gap, projected_position); break; } case _not_defined: { AKANTU_DEBUG_ERROR("Not a valid surface element type : " << type); break; } case _segment_3: case _triangle_6: case _tetrahedron_4: case _tetrahedron_10: case _quadrangle_4: case _quadrangle_8: case _hexahedron_8: case _point: case _bernoulli_beam_2: + case _cohesive_2d_4: + case _cohesive_2d_6: case _max_element_type: { AKANTU_DEBUG_ERROR("Contact search is not implemented for this surface element type : " << type); break; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactSearchExplicit::checkPenetrationSituation(const UInt impactor_node, const UInt surface_element, const ElementType type, bool & is_inside, bool & is_in_projection_area) { AKANTU_DEBUG_IN(); switch(type) { case _segment_2: { checkPenetrationSituationSegment2(impactor_node, surface_element, is_inside, is_in_projection_area); break; } case _triangle_3: { checkPenetrationSituationTriangle3(impactor_node, surface_element, is_inside, is_in_projection_area); break; } case _not_defined: { AKANTU_DEBUG_ERROR("Not a valid surface element type : " << type); break; } case _segment_3: case _triangle_6: case _tetrahedron_4: case _tetrahedron_10: case _quadrangle_4: case _quadrangle_8: case _hexahedron_8: case _point: case _bernoulli_beam_2: + case _cohesive_2d_4: + case _cohesive_2d_6: case _max_element_type: { AKANTU_DEBUG_ERROR("Contact search is not implemented for this surface element type : " << type); break; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactSearchExplicit::computeComponentsOfProjectionSegment2(const UInt impactor_node, const UInt surface_element, Real * normal, Real & gap, Real * projected_position) { AKANTU_DEBUG_IN(); const UInt dim = spatial_dimension; const ElementType type = _segment_2; const UInt nb_nodes_element = Mesh::getNbNodesPerElement(type); Real * current_position = contact.getModel().getCurrentPosition().values; UInt * connectivity = mesh.getConnectivity(type, _not_ghost).values; UInt node_1 = surface_element * nb_nodes_element; Real * position_node_1 = &(current_position[connectivity[node_1 + 0] * spatial_dimension]); Real * position_node_2 = &(current_position[connectivity[node_1 + 1] * spatial_dimension]); Real * position_impactor_node = &(current_position[impactor_node * spatial_dimension]); /// compute the normal of the face Real vector_1[2]; Math::vector_2d(position_node_1, position_node_2, vector_1); /// @todo: check if correct order of nodes !! Math::normal2(vector_1, normal); /// compute the gap between impactor and face /// gap positive if impactor outside of surface Real vector_node_1_impactor[2]; Math::vector_2d(position_node_1, position_impactor_node, vector_node_1_impactor); gap = Math::vectorDot2(vector_node_1_impactor, normal); /// compute the projected position of the impactor node onto the face for(UInt i=0; i < dim; ++i) { projected_position[i] = position_impactor_node[i] - gap * normal[i]; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactSearchExplicit::computeComponentsOfProjectionTriangle3(const UInt impactor_node, const UInt surface_element, Real * normal, Real & gap, Real * projected_position) { AKANTU_DEBUG_IN(); const UInt dim = spatial_dimension; const ElementType type = _triangle_3; const UInt nb_nodes_element = Mesh::getNbNodesPerElement(type); Real * current_position = contact.getModel().getCurrentPosition().values; UInt * connectivity = mesh.getConnectivity(type, _not_ghost).values; UInt node_1 = surface_element * nb_nodes_element; Real * position_node_1 = &(current_position[connectivity[node_1 + 0] * spatial_dimension]); Real * position_node_2 = &(current_position[connectivity[node_1 + 1] * spatial_dimension]); Real * position_node_3 = &(current_position[connectivity[node_1 + 2] * spatial_dimension]); Real * position_impactor_node = &(current_position[impactor_node * spatial_dimension]); /// compute the normal of the face Real vector_1[3]; Real vector_2[3]; Math::vector_3d(position_node_1, position_node_2, vector_1); /// @todo: check if correct order of nodes !! Math::vector_3d(position_node_1, position_node_3, vector_2); Math::normal3(vector_1, vector_2, normal); /// compute the gap between impactor and face /// gap positive if impactor outside of surface Real vector_node_1_impactor[3]; Math::vector_3d(position_node_1, position_impactor_node, vector_node_1_impactor); gap = Math::vectorDot3(vector_node_1_impactor, normal); /// compute the projected position of the impactor node onto the face for(UInt i=0; i < dim; ++i) { projected_position[i] = position_impactor_node[i] - gap * normal[i]; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactSearchExplicit::checkPenetrationSituationSegment2(const UInt impactor_node, const UInt surface_element, bool & is_inside, bool & is_in_projection_area) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(spatial_dimension == 2, "wrong spatial dimension (=" << spatial_dimension << ") for checkPenetrationSituationSegment2"); const ElementType type = _segment_2; const UInt nb_nodes_element = Mesh::getNbNodesPerElement(type); const Real tolerance = std::numeric_limits::epsilon(); Real * current_position = contact.getModel().getCurrentPosition().values; UInt * connectivity = mesh.getConnectivity(type, _not_ghost).values; Real gap; Real normal[2]; Real projected_position[2]; computeComponentsOfProjectionSegment2(impactor_node, surface_element, normal, gap, projected_position); // ------------------------------------------------------- /// Find if impactor node is inside or outside of the face // ------------------------------------------------------- if(gap < -tolerance) is_inside = true; else is_inside = false; // ---------------------------------------------------- /// Find if impactor node is in projection area of face // ---------------------------------------------------- UInt node_1 = surface_element * nb_nodes_element; Real * position_node_1 = &(current_position[connectivity[node_1 + 0] * spatial_dimension]); Real * position_node_2 = &(current_position[connectivity[node_1 + 1] * spatial_dimension]); Real tmp_vector_1_imp[2]; Real tmp_vector_1_2[2]; // find vectors from master node 1 to impactor and master node 2 Math::vector_2d(position_node_1, position_node_2, tmp_vector_1_2); Math::vector_2d(position_node_1, projected_position, tmp_vector_1_imp); Real length_difference = Math::norm2(tmp_vector_1_imp) - Math::norm2(tmp_vector_1_2); // the projection is outside if the test area is larger than the area of the triangle if(length_difference > tolerance) is_in_projection_area = false; else is_in_projection_area = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactSearchExplicit::checkPenetrationSituationTriangle3(const UInt impactor_node, const UInt surface_element, bool & is_inside, bool & is_in_projection_area) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(spatial_dimension == 3, "wrong spatial dimension (=" << spatial_dimension << ") for checkPenetrationSituationTriangle3"); const ElementType type = _triangle_3; const UInt nb_nodes_element = Mesh::getNbNodesPerElement(type); const Real tolerance = std::numeric_limits::epsilon(); Real * current_position = contact.getModel().getCurrentPosition().values; UInt * connectivity = mesh.getConnectivity(type, _not_ghost).values; Real gap; Real normal[3]; Real projected_position[3]; computeComponentsOfProjectionTriangle3(impactor_node, surface_element, normal, gap, projected_position); // ------------------------------------------------------- /// Find if impactor node is inside or outside of the face // ------------------------------------------------------- if(gap < -tolerance) is_inside = true; else is_inside = false; // ---------------------------------------------------- /// Find if impactor node is in projection area of face // ---------------------------------------------------- UInt node_1 = surface_element * nb_nodes_element; Real * position_node_1 = &(current_position[connectivity[node_1 + 0] * spatial_dimension]); Real * position_node_2 = &(current_position[connectivity[node_1 + 1] * spatial_dimension]); Real * position_node_3 = &(current_position[connectivity[node_1 + 2] * spatial_dimension]); Real triangle_area; Real test_area = 0.0; Real tmp_vector_1[3]; Real tmp_vector_2[3]; Real tmp_vector_3[3]; // find area of triangle Math::vector_3d(position_node_1, position_node_2, tmp_vector_1); Math::vector_3d(position_node_1, position_node_3, tmp_vector_2); Math::vectorProduct3(tmp_vector_1, tmp_vector_2, tmp_vector_3); triangle_area = 0.5 * Math::norm3(tmp_vector_3); // compute areas of projected position and two nodes of master triangle UInt nb_sub_areas = nb_nodes_element; UInt node_order[4]; node_order[0] = 0; node_order[1] = 1; node_order[2] = 2; node_order[3] = 0; for(UInt i=0; i < nb_sub_areas; ++i) { position_node_1 = &(current_position[connectivity[node_1 + node_order[i ]] * spatial_dimension]); position_node_2 = &(current_position[connectivity[node_1 + node_order[i+1]] * spatial_dimension]); Math::vector_3d(projected_position, position_node_1, tmp_vector_1); Math::vector_3d(projected_position, position_node_2, tmp_vector_2); Math::vectorProduct3(tmp_vector_1, tmp_vector_2, tmp_vector_3); test_area += 0.5 * Math::norm3(tmp_vector_3); } // the projection is outside if the test area is larger than the area of the triangle if((test_area - triangle_area) > tolerance) is_in_projection_area = false; else is_in_projection_area = true; AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/src/model/solid_mechanics/contact/friction_coefficient/velocity_dependent_fric_coef.cc b/src/model/solid_mechanics/contact/friction_coefficient/velocity_dependent_fric_coef.cc index b4a29ad2f..7b1ba593f 100644 --- a/src/model/solid_mechanics/contact/friction_coefficient/velocity_dependent_fric_coef.cc +++ b/src/model/solid_mechanics/contact/friction_coefficient/velocity_dependent_fric_coef.cc @@ -1,169 +1,171 @@ /** * @file velocity_dependent_fric_coef.cc * @author David Kammer * @date Mon Jun 20 15:24:21 2011 * * @brief implementation for velocity dependence of friction coefficient * * @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 "velocity_dependent_fric_coef.hh" #include "contact_rigid.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ VelocityDependentFricCoef::VelocityDependentFricCoef(ContactRigid & contact, const Surface & master_surface) : FrictionCoefficient(contact, master_surface), spatial_dimension(this->contact.getSpatialDimension()) { AKANTU_DEBUG_IN(); this->relative_sliding_velocities = new Vector(0,1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ VelocityDependentFricCoef::~VelocityDependentFricCoef() { AKANTU_DEBUG_IN(); delete this->relative_sliding_velocities; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void VelocityDependentFricCoef::initializeComputeFricCoef() { AKANTU_DEBUG_IN(); // find impactor_information for given master const ContactRigid::SurfaceToImpactInfoMap & imp_info = this->contact.getImpactorsInformation(); ContactRigid::SurfaceToImpactInfoMap::const_iterator it; it = imp_info.find(this->master_surface); AKANTU_DEBUG_ASSERT(it != imp_info.end(), "Couldn't find impactor information object for master surface " << master_surface); ContactRigid::ImpactorInformationPerMaster * impactor_info = it->second; Vector * active_nodes = impactor_info->active_impactor_nodes; UInt * active_nodes_val = active_nodes->values; Vector * master_normals = impactor_info->master_normals; Real * master_normals_val = master_normals->values; Real * velocity_val = this->contact.getModel().getVelocity().values; // resize the relative sliding velocity vector to the nb of active impactor nodes this->relative_sliding_velocities->resize(active_nodes->getSize()); Real * relative_sliding_velocities_val = this->relative_sliding_velocities->values; for(UInt n = 0; n < active_nodes->getSize(); ++n) { UInt impactor = active_nodes_val[n]; UInt master_normal_index = n * this->spatial_dimension; UInt velocity_index = impactor * this->spatial_dimension; Real projected_length; if(this->spatial_dimension == 2) projected_length = Math::vectorDot2(&master_normals_val[master_normal_index], &velocity_val[velocity_index]); else if(this->spatial_dimension == 3) projected_length = Math::vectorDot3(&master_normals_val[master_normal_index], &velocity_val[velocity_index]); Real tangential_impactor_velocity[3]; for(UInt i=0; ispatial_dimension; ++i) { tangential_impactor_velocity[i] = velocity_val[velocity_index + i] - projected_length * master_normals_val[master_normal_index + i]; } Real tangential_master_velocity[3]; //Real * tangential_master_velocity_p = &(tangential_master_velocity[0]); computeTangentialMasterVelocity(n, impactor_info, tangential_master_velocity); Real relative_sliding_speed; Real relative_sliding_velocity[3]; //Real * relative_sliding_velocity_p = &(relative_sliding_velocity_p[0]); for(UInt i =0; ispatial_dimension; ++i) relative_sliding_velocity[i] = tangential_impactor_velocity[i] - tangential_master_velocity[i]; if(this->spatial_dimension == 2) relative_sliding_speed = Math::norm2(relative_sliding_velocity); else if(this->spatial_dimension == 3) relative_sliding_speed = Math::norm3(relative_sliding_velocity); relative_sliding_velocities_val[n] = relative_sliding_speed; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void VelocityDependentFricCoef::computeTangentialMasterVelocity(UInt impactor_index, ContactRigid::ImpactorInformationPerMaster * impactor_info, Real * tangential_master_velocity) { AKANTU_DEBUG_IN(); // for ContactRigid, the master velocity is by definition zero const ContactType contact_type = this->contact.getType(); if(contact_type == _ct_rigid /* do not change for general contact */) { for (UInt i=0; ispatial_dimension; ++i) tangential_master_velocity[i] = 0.; } else { // compute the tangential master velocity ElementType type = impactor_info->master_element_type->at(impactor_index); switch(type) { /* case _segment_2: { //computeComponentsOfProjectionSegment2(impactor_node, surface_element, normal, gap, projected_position); break; } */ /* case _triangle_3: { //computeComponentsOfProjectionTriangle3(impactor_node, surface_element, normal, gap, projected_position); break; } */ case _not_defined: { AKANTU_DEBUG_ERROR("Not a valid surface element type : " << type << " for computation of tangential velocity of master element"); break; } case _segment_2: case _segment_3: case _triangle_3: case _triangle_6: case _tetrahedron_4: case _tetrahedron_10: case _quadrangle_4: case _quadrangle_8: case _hexahedron_8: case _point: case _bernoulli_beam_2: + case _cohesive_2d_4: + case _cohesive_2d_6: case _max_element_type: { AKANTU_DEBUG_ERROR("Not a valid surface element type : " << type << " for computation of tangential velocity of master element"); break; } } } AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index f56cdc5ef..45b59ea62 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,906 +1,914 @@ /** * @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 "dof_synchronizer.hh" #ifdef AKANTU_USE_MUMPS #include "solver_mumps.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh, UInt dim, const ID & 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), element_material("element_material", id), integrator(NULL), increment_flag(false), solver(NULL), spatial_dimension(dim), mesh(mesh), dynamic(false), implicit(false) { 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; 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 && stiffness_matrix != jacobian_matrix) delete stiffness_matrix; if(jacobian_matrix) delete jacobian_matrix; if(dof_synchronizer) delete dof_synchronizer; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initFull(std::string material_file, bool implicit_scheme, bool implicit_dynamic) { // initialize the model initModel(); // initialize the vectors initVectors(); // set the initial condition to 0 force->clear(); velocity->clear(); acceleration->clear(); displacement->clear(); + // initialize pcb + if(pbc_pair.size()!=0) + initPBC(); + // initialize the time integration schemes if(implicit_scheme) initImplicit(implicit_dynamic); else initExplicit(); // initialize the materials if(material_file != "") { readMaterials(material_file); initMaterials(); } } 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::initFEMBoundary(bool create_surface) { if(create_surface) MeshUtils::buildFacets(mesh, true, false); FEM & fem_boundary = getFEMBoundary(); fem_boundary.initShapeFunctions(); fem_boundary.computeNormalsOnControlPoints(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initExplicit() { AKANTU_DEBUG_IN(); if (integrator) delete integrator; integrator = new CentralDifference(); dynamic = true; implicit = false; 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, 0)); velocity = &(alloc(sstr_velo.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); acceleration = &(alloc(sstr_acce.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); force = &(alloc(sstr_forc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); residual = &(alloc(sstr_resi.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); boundary = &(alloc(sstr_boun.str(), nb_nodes, spatial_dimension, false)); std::stringstream sstr_curp; sstr_curp << id << ":current_position"; current_position = &(alloc(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE)); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(gt); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { UInt nb_element = mesh.getNbElement(*it, gt); element_material.alloc(nb_element, 1, *it, gt); } } dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); + dof_synchronizer->initGlobalDOFEquationNumbers(); 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); +void SolidMechanicsModel::initPBC() { + Model::initPBC(); registerPBCSynchronizer(); -} -/* -------------------------------------------------------------------------- */ -void SolidMechanicsModel::initPBC(std::list< std::pair > & surface_pairs, - ElementType surface_e_type){ - Model::initPBC(surface_pairs, surface_e_type); - registerPBCSynchronizer(); + // as long as there are ones on the diagonal of the matrix, we can put boudandary true for slaves + std::map::iterator it = pbc_pair.begin(); + std::map::iterator end = pbc_pair.end(); + UInt dim = mesh.getSpatialDimension(); + while(it != end) { + for (UInt i=0; iregisterSynchronizer(*synch, _gst_smm_uv); synch_registry->registerSynchronizer(*synch, _gst_smm_mass); 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(); // f = f_ext - f_int - Ma - Cv // f = f_ext if (need_initialize) initializeUpdateResidualData(); // f -= fint /// start synchronization synch_registry->asynchronousSynchronize(_gst_smm_uv); synch_registry->waitEndSynchronize(_gst_smm_uv); 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); } if(dynamic) { // f -= Ma if(mass_matrix) { // if full mass_matrix Vector * Ma = new Vector(*acceleration, true, "Ma"); *Ma *= *mass_matrix; *residual -= *Ma; delete Ma; } else { // else lumped mass UInt nb_nodes = acceleration->getSize(); UInt nb_degre_of_freedom = acceleration->getNbComponent(); Real * mass_val = mass->values; Real * accel_val = acceleration->values; Real * res_val = residual->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes * nb_degre_of_freedom; ++n) { if(!(*boundary_val)) { *res_val -= *accel_val * *mass_val; } boundary_val++; res_val++; mass_val++; accel_val++; } } // f -= Cv if(velocity_damping_matrix) { Vector * Cv = new Vector(*velocity); *Cv *= *velocity_damping_matrix; *residual -= *Cv; delete Cv; } } 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(); if(!implicit && !mass_matrix) { Real * mass_val = mass->values; Real * residual_val = residual->values; bool * boundary_val = boundary->values; Real * inc = increment_acceleration->values; Real * accel_val = 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); } residual_val++; boundary_val++; inc++; mass_val++; accel_val++; } } } else { solveDynamic(*increment_acceleration); } 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(); integrator->integrationSchemeCorrAccel(time_step, *displacement, *velocity, *acceleration, *boundary, *increment_acceleration); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Implicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initSolver() { #if !defined(AKANTU_USE_MUMPS) // or other solver in the future \todo add AKANTU_HAS_SOLVER in CMake AKANTU_DEBUG_ERROR("You should at least activate one solver."); #else UInt nb_global_node = mesh.getNbGlobalNodes(); std::stringstream sstr; sstr << id << ":jacobian_matrix"; jacobian_matrix = new SparseMatrix(nb_global_node * spatial_dimension, _symmetric, spatial_dimension, sstr.str(), memory_id); - dof_synchronizer->initGlobalDOFEquationNumbers(); + // dof_synchronizer->initGlobalDOFEquationNumbers(); jacobian_matrix->buildProfile(mesh, *dof_synchronizer); #ifdef AKANTU_USE_MUMPS std::stringstream sstr_solv; sstr_solv << id << ":solver"; solver = new SolverMumps(*jacobian_matrix, sstr_solv.str()); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS solver->initialize(); #endif //AKANTU_HAS_SOLVER } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initImplicit(bool dynamic) { AKANTU_DEBUG_IN(); this->dynamic = dynamic; implicit = true; initSolver(); if(dynamic) { if(integrator) delete integrator; integrator = new TrapezoidalRule2(); std::stringstream sstr; sstr << id << ":stiffness_matrix"; stiffness_matrix = new SparseMatrix(*jacobian_matrix, sstr.str(), memory_id); } else { stiffness_matrix = jacobian_matrix; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initialAcceleration() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving Ma = f"); Solver * acc_solver = NULL; #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_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"); if(!increment) setIncrementFlagOn(); solveDynamic(*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); if(dynamic) jacobian_matrix->copyContent(*stiffness_matrix); 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(); error = 0; Real norm[2] = {0., 0.}; Real * increment_val = increment->storage(); bool * boundary_val = boundary->storage(); Real * displacement_val = displacement->storage(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); for (UInt d = 0; d < nb_degre_of_freedom; ++d) { if(!(*boundary_val) && is_local_node) { norm[0] += *increment_val * *increment_val; norm[1] += *displacement_val * *displacement_val; } boundary_val++; increment_val++; displacement_val++; } } StaticCommunicator::getStaticCommunicator()->allReduce(norm, 2, _so_sum); norm[0] = sqrt(norm[0]); norm[1] = sqrt(norm[1]); AKANTU_DEBUG_ASSERT(!Math::isnan(norm[0]), "Something goes wrong in the solve phase"); std::cout << norm[0] << " " << norm[1] << std::endl; AKANTU_DEBUG_OUT(); error = norm[0] / norm[1]; 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 & norm) { AKANTU_DEBUG_IN(); UInt nb_nodes = residual->getSize(); norm = 0; Real * residual_val = residual->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); if(is_local_node) { for (UInt d = 0; d < spatial_dimension; ++d) { if(!(*boundary_val)) { norm += *residual_val * *residual_val; } boundary_val++; residual_val++; } } else { boundary_val += spatial_dimension; residual_val += spatial_dimension; } } StaticCommunicator::getStaticCommunicator()->allReduce(&norm, 1, _so_sum); norm = sqrt(norm); AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); return (norm < tolerance); } /* -------------------------------------------------------------------------- */ 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(); Real min_dt = getStableTimeStep(_not_ghost); /// reduction min over all processors StaticCommunicator::getStaticCommunicator()->allReduce(&min_dt, 1, _so_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); Material ** mat_val = &(materials.at(0)); Real min_dt = std::numeric_limits::max(); Real * coord = mesh.getNodes().values; Real * disp_val = displacement->values; Element elem; elem.ghost_type = ghost_type; const Mesh::ConnectivityTypeList & type_list = mesh.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; elem.type = *it; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*it); UInt nb_element = mesh.getNbElement(*it); UInt * conn = mesh.getConnectivity(*it, ghost_type).storage(); UInt * elem_mat_val = element_material(*it, ghost_type).storage(); 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; } 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 StaticCommunicator::getStaticCommunicator()->allReduce(&epot, 1, _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 mv2 = 0; bool is_local_node = mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !getIsPBCSlaveNode(n); for (UInt i = 0; i < spatial_dimension; ++i) { // if(is_local_node) { mv2 += is_local_node * is_not_pbc_slave_node * *vel_val * *vel_val * *mass_val; // } vel_val++; mass_val++; } ekin += mv2; } StaticCommunicator::getStaticCommunicator()->allReduce(&ekin, 1, _so_sum); AKANTU_DEBUG_OUT(); return ekin * .5; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Solid Mechanics Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEM().printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + nodals information [" << std::endl; displacement->printself(stream, indent + 2); mass ->printself(stream, indent + 2); velocity ->printself(stream, indent + 2); acceleration->printself(stream, indent + 2); force ->printself(stream, indent + 2); residual ->printself(stream, indent + 2); boundary ->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + connectivity type information [" << std::endl; element_material.printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index 5477992e1..e756f2cef 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,513 +1,511 @@ /** * @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 "data_accessor.hh" #include "material.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #include "aka_types.hh" #include "integration_scheme_2nd_order.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(Mesh & mesh, UInt spatial_dimension = 0, const ID & id = "solid_mechanics_model", const MemoryID & memory_id = 0); virtual ~SolidMechanicsModel(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize completely the model void initFull(std::string material_file = "", bool implicit_scheme = false, bool implicit_dynamic = false); /// initialize the fem object needed for boundary conditions void initFEMBoundary(bool create_surface = true); /// 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); - void initPBC(std::list< std::pair > & surface_pairs, - ElementType surface_e_type); + void initPBC(); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* PBC */ /* ------------------------------------------------------------------------ */ public: /// change the equation number for proper assembly when using PBC void changeEquationNumberforPBC(std::map & pbc_pair); private: /// register PBC synchronizer void registerPBCSynchronizer(); /* ------------------------------------------------------------------------ */ /* Explicit */ /* ------------------------------------------------------------------------ */ public: /// initialize the stuff for the explicit scheme void initExplicit(); /// 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); /** * \brief compute the acceleration from the residual * this function is the explicit equivalent to solveDynamic in implicit * In the case of lumped mass just divide the residual by the mass * In the case of not lumped mass call solveDynamic<_acceleration_corrector> */ void updateAcceleration(); /// explicit integration predictor void explicitPred(); /// explicit integration corrector void explicitCorr(); /* ------------------------------------------------------------------------ */ /* Implicit */ /* ------------------------------------------------------------------------ */ public: /// initialize the solver and the jacobian_matrix (called by initImplicit) void initSolver(); /// 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] in displacement 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(); private: /// compute A and solve @f[ A\delta u = f_ext - f_int @f] template void solveDynamic(Vector & increment); /* ------------------------------------------------------------------------ */ /* Boundaries (solid_mechanics_model_boundary.cc) */ /* ------------------------------------------------------------------------ */ public: class SurfaceLoadFunctor { public: virtual void operator()(const types::Vector & position, types::Vector & force, const types::Vector & normal, Surface surface_id) { AKANTU_DEBUG_TO_IMPLEMENT(); } virtual void operator()(const types::Vector & position, types::Matrix & stress, const types::Vector & normal, Surface surface_id) { AKANTU_DEBUG_TO_IMPLEMENT(); } }; /// compute force vector from a function(x,y,z) that describe stresses void computeForcesFromFunction(BoundaryFunction in_function, BoundaryFunctionType function_type) __attribute__((deprecated)); template void computeForcesFromFunction(Functor & functor, BoundaryFunctionType function_type); /// integrate a force on the boundary by providing a stress tensor void computeForcesByStressTensor(const Vector & stresses, const ElementType & type, const GhostType & ghost_type); /// integrate a force on the boundary by providing a traction vector void computeForcesByTractionVector(const Vector & tractions, const ElementType & type, const GhostType & ghost_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, ID 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); /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: __aka_inline__ virtual UInt getNbDataToPack(const Element & element, SynchronizationTag tag) const; __aka_inline__ virtual UInt getNbDataToUnpack(const Element & element, SynchronizationTag tag) const; __aka_inline__ virtual void packData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag) const; __aka_inline__ virtual void unpackData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag); __aka_inline__ virtual UInt getNbDataToPack(SynchronizationTag tag) const; __aka_inline__ virtual UInt getNbDataToUnpack(SynchronizationTag tag) const; __aka_inline__ virtual void packData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag) const; __aka_inline__ virtual void unpackData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* 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 SolidMechnicsModel::incrementAcceleration vector AKANTU_GET_MACRO(IncrementAcceleration, *increment_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, UInt); /// get a particular material __aka_inline__ Material & getMaterial(UInt mat_index); __aka_inline__ const Material & getMaterial(UInt mat_index) const; /// give the number of materials __aka_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 &); __aka_inline__ FEM & getFEMBoundary(std::string name = ""); private: /// compute the stable time step Real getStableTimeStep(const GhostType & ghost_type); /* ------------------------------------------------------------------------ */ /* 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 \Delta t^2}, d = \frac{\gamma}{\beta \Delta t} @f] SparseMatrix * jacobian_matrix; /// materials of all elements ByElementTypeUInt 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; bool dynamic; bool implicit; }; __END_AKANTU__ /* -------------------------------------------------------------------------- */ /* __aka_inline__ functions */ /* -------------------------------------------------------------------------- */ #include "parser.hh" __BEGIN_AKANTU__ #include "solid_mechanics_model_tmpl.hh" #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "solid_mechanics_model_inline_impl.cc" #endif /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const SolidMechanicsModel & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */ diff --git a/src/solver/sparse_matrix.cc b/src/solver/sparse_matrix.cc index c7675b874..026e36ce7 100644 --- a/src/solver/sparse_matrix.cc +++ b/src/solver/sparse_matrix.cc @@ -1,419 +1,430 @@ /** * @file sparse_matrix.cc * @author Nicolas Richart * @date Tue Oct 26 18:25:07 2010 * * @brief implementation of the SparseMatrix 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 /* -------------------------------------------------------------------------- */ #include "sparse_matrix.hh" #include "static_communicator.hh" #include "dof_synchronizer.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ SparseMatrix::SparseMatrix(UInt size, const SparseMatrixType & sparse_matrix_type, UInt nb_degre_of_freedom, const ID & id, const MemoryID & memory_id) : Memory(memory_id), id(id), sparse_matrix_type(sparse_matrix_type), nb_degre_of_freedom(nb_degre_of_freedom), size(size), nb_non_zero(0), irn(0,1,"irn"), jcn(0,1,"jcn"), a(0,1,"A") { AKANTU_DEBUG_IN(); StaticCommunicator * comm = StaticCommunicator::getStaticCommunicator(); nb_proc = comm->getNbProc(); dof_synchronizer = NULL; irn_save = NULL; jcn_save = NULL; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SparseMatrix::SparseMatrix(const SparseMatrix & matrix, const ID & id, const MemoryID & memory_id) : Memory(memory_id), id(id), sparse_matrix_type(matrix.sparse_matrix_type), nb_degre_of_freedom(matrix.nb_degre_of_freedom), size(matrix.size), nb_proc(matrix.nb_proc), nb_non_zero(matrix.nb_non_zero), irn(matrix.irn, true), jcn(matrix.jcn, true), a(matrix.getA(), true), irn_jcn_k(matrix.irn_jcn_k) { AKANTU_DEBUG_IN(); size_save = 0; irn_save = NULL; jcn_save = NULL; dof_synchronizer = matrix.dof_synchronizer; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SparseMatrix::~SparseMatrix() { AKANTU_DEBUG_IN(); // if (irn_jcn_to_k) delete irn_jcn_to_k; if(irn_save) delete irn_save; if(jcn_save) delete jcn_save; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrix::buildProfile(const Mesh & mesh, const DOFSynchronizer & dof_synchronizer) { AKANTU_DEBUG_IN(); // if(irn_jcn_to_k) delete irn_jcn_to_k; // irn_jcn_to_k = new std::map, UInt>; clearProfile(); this->dof_synchronizer = &const_cast(dof_synchronizer); coordinate_list_map::iterator irn_jcn_k_it; Int * eq_nb_val = dof_synchronizer.getGlobalDOFEquationNumbers().values; const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if (Mesh::getSpatialDimension(*it) != mesh.getSpatialDimension()) continue; UInt nb_element = mesh.getNbElement(*it); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); UInt size_mat = nb_nodes_per_element * nb_degre_of_freedom; UInt * conn_val = mesh.getConnectivity(*it, _not_ghost).values; Int * local_eq_nb_val = new Int[nb_degre_of_freedom * nb_nodes_per_element]; for (UInt e = 0; e < nb_element; ++e) { Int * tmp_local_eq_nb_val = local_eq_nb_val; 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 = (sparse_matrix_type == _symmetric) ? i : 0; for (UInt j = j_start; j < size_mat; ++j) { UInt c_jcn = local_eq_nb_val[j]; if(c_jcn < size) { KeyCOO irn_jcn = key(c_irn, c_jcn); irn_jcn_k_it = irn_jcn_k.find(irn_jcn); if (irn_jcn_k_it == irn_jcn_k.end()) { irn_jcn_k[irn_jcn] = nb_non_zero; irn.push_back(c_irn + 1); jcn.push_back(c_jcn + 1); nb_non_zero++; } } } } } conn_val += nb_nodes_per_element; } delete [] local_eq_nb_val; } + for (UInt i = 0; i < size; ++i) { + KeyCOO irn_jcn = key(i, i); + coordinate_list_map::const_iterator irn_jcn_k_it = irn_jcn_k.find(irn_jcn); + if(irn_jcn_k_it == irn_jcn_k.end()) { + irn_jcn_k[irn_jcn] = nb_non_zero; + irn.push_back(i + 1); + jcn.push_back(i + 1); + nb_non_zero++; + } + } + a.resize(nb_non_zero); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrix::applyBoundary(const Vector & boundary) { AKANTU_DEBUG_IN(); const DOFSynchronizer::GlobalEquationNumberMap & local_eq_num_to_global = dof_synchronizer->getGlobalEquationNumberToLocal(); Int * irn_val = irn.values; Int * jcn_val = jcn.values; Real * a_val = a.values; for (UInt i = 0; i < nb_non_zero; ++i) { UInt ni = local_eq_num_to_global.find(*irn_val - 1)->second; UInt nj = local_eq_num_to_global.find(*jcn_val - 1)->second; if(boundary.values[ni] || boundary.values[nj]) { if (*irn_val != *jcn_val) *a_val = 0; else { if(dof_synchronizer->getDOFTypes()(ni) >= 0) *a_val = 0; else *a_val = 1; } } irn_val++; jcn_val++; a_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrix::removeBoundary(const Vector & boundary) { AKANTU_DEBUG_IN(); if(irn_save) delete irn_save; if(jcn_save) delete jcn_save; irn_save = new Vector(irn, true); jcn_save = new Vector(jcn, true); UInt n = boundary.getSize()*boundary.getNbComponent(); UInt * perm = new UInt[n]; size_save = size; size = 0; for (UInt i = 0; i < n; ++i) { if(!boundary.values[i]) { perm[i] = size; // std::cout << "perm["<< i <<"] = " << size << std::endl; size++; } } for (UInt i = 0; i < nb_non_zero;) { if(boundary.values[irn.at(i) - 1] || boundary.values[jcn.at(i) - 1]) { irn.erase(i); jcn.erase(i); a.erase(i); nb_non_zero--; } else { irn.values[i] = perm[irn.values[i] - 1] + 1; jcn.values[i] = perm[jcn.values[i] - 1] + 1; i++; } } delete [] perm; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrix::restoreProfile() { AKANTU_DEBUG_IN(); irn.resize(irn_save->getSize()); jcn.resize(jcn_save->getSize()); nb_non_zero = irn.getSize(); a.resize(nb_non_zero); size = size_save; memcpy(irn.values, irn_save->values, irn.getSize()*sizeof(Int)); memcpy(jcn.values, jcn_save->values, jcn.getSize()*sizeof(Int)); delete irn_save; irn_save = NULL; delete jcn_save; jcn_save = NULL; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrix::saveProfile(const std::string & filename) const { AKANTU_DEBUG_IN(); std::ofstream outfile; outfile.open(filename.c_str()); outfile << "%%MatrixMarket matrix coordinate pattern"; if(sparse_matrix_type == _symmetric) outfile << " symmetric"; else outfile << " general"; outfile << std::endl; UInt m = size; outfile << m << " " << m << " " << nb_non_zero << std::endl; for (UInt i = 0; i < nb_non_zero; ++i) { outfile << irn.values[i] << " " << jcn.values[i] << " 1" << std::endl; } outfile.close(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrix::saveMatrix(const std::string & filename) const { AKANTU_DEBUG_IN(); std::ofstream outfile; outfile.precision(std::numeric_limits::digits10); outfile.open(filename.c_str()); outfile << "%%MatrixMarket matrix coordinate real"; if(sparse_matrix_type == _symmetric) outfile << " symmetric"; else outfile << " general"; outfile << std::endl; outfile << size << " " << size << " " << nb_non_zero << std::endl; for (UInt i = 0; i < nb_non_zero; ++i) { outfile << irn(i) << " " << jcn(i) << " " << a(i) << std::endl; } outfile.close(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Vector & operator*=(Vector & vect, const SparseMatrix & mat) { AKANTU_DEBUG_IN(); // AKANTU_DEBUG_ASSERT((vect.getSize()*vect.getNbComponent() == mat.getSize()) && // (vect.getNbComponent() == mat.getNbDegreOfFreedom()), // "The size of the matrix and the vector do not match"); const SparseMatrixType & sparse_matrix_type = mat.getSparseMatrixType(); DOFSynchronizer * dof_synchronizer = mat.getDOFSynchronizerPointer(); UInt nb_non_zero = mat.getNbNonZero(); Real * tmp = new Real [vect.getNbComponent() * vect.getSize()]; std::fill_n(tmp, vect.getNbComponent() * vect.getSize(), 0); Int * i_val = mat.getIRN().values; Int * j_val = mat.getJCN().values; Real * a_val = mat.getA().values; Real * vect_val = vect.values; for (UInt k = 0; k < nb_non_zero; ++k) { UInt i = *(i_val++); UInt j = *(j_val++); Real a = *(a_val++); UInt local_i = i - 1; UInt local_j = j - 1; if(dof_synchronizer) { local_i = dof_synchronizer->getDOFLocalID(local_i); local_j = dof_synchronizer->getDOFLocalID(local_j); } tmp[local_i] += a * vect_val[local_j]; if(sparse_matrix_type == _symmetric && (local_i != local_j)) tmp[local_j] += a * vect_val[local_i]; } memcpy(vect_val, tmp, vect.getNbComponent() * vect.getSize() * sizeof(Real)); delete [] tmp; if(dof_synchronizer) dof_synchronizer->reduceSynchronize >(vect); AKANTU_DEBUG_OUT(); return vect; } /* -------------------------------------------------------------------------- */ void SparseMatrix::copyContent(const SparseMatrix & matrix) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(nb_non_zero == matrix.getNbNonZero(), "The to matrix don't have the same profiles"); memcpy(a.values, matrix.getA().values, nb_non_zero * sizeof(Real)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrix::add(const SparseMatrix & matrix, Real alpha) { AKANTU_DEBUG_ASSERT(nb_non_zero == matrix.getNbNonZero(), "The to matrix don't have the same profiles"); Real * a_val = a.values; Real * b_val = matrix.a.values; for (UInt n = 0; n < nb_non_zero; ++n) { *a_val++ += alpha * *b_val++; } } /* -------------------------------------------------------------------------- */ void SparseMatrix::lump(Vector & lumped) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT((lumped.getNbComponent() == nb_degre_of_freedom), "The size of the matrix and the vector do not match"); UInt vect_size = size / nb_degre_of_freedom; if(dof_synchronizer) vect_size = dof_synchronizer->getNbDOFs() / nb_degre_of_freedom; lumped.resize(vect_size); lumped.clear(); Int * i_val = irn.values; Int * j_val = jcn.values; Real * a_val = a.values; Real * vect_val = lumped.values; for (UInt k = 0; k < nb_non_zero; ++k) { UInt i = *(i_val++); UInt j = *(j_val++); Real a = *(a_val++); UInt local_i = i - 1; UInt local_j = j - 1; if(dof_synchronizer) { local_i = dof_synchronizer->getDOFLocalID(local_i); local_j = dof_synchronizer->getDOFLocalID(local_j); } vect_val[local_i] += a; if(sparse_matrix_type == _symmetric && (i != j)) vect_val[local_j] += a; } if(dof_synchronizer) dof_synchronizer->reduceSynchronize >(lumped); AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/src/synchronizer/dof_synchronizer.hh b/src/synchronizer/dof_synchronizer.hh index 9046a61da..6798c34f9 100644 --- a/src/synchronizer/dof_synchronizer.hh +++ b/src/synchronizer/dof_synchronizer.hh @@ -1,201 +1,202 @@ /** * @file dof_synchronizer.hh * @author Nicolas Richart * @date Mon May 23 16:50:47 2011 * * @brief Synchronize Vector of DOFs * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_vector.hh" #include "static_communicator.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_DOF_SYNCHRONIZER_HH__ #define __AKANTU_DOF_SYNCHRONIZER_HH__ __BEGIN_AKANTU__ class Mesh; template class AddOperation { public: __aka_inline__ T operator()(T & b, T & a) { return a + b; }; }; class DOFSynchronizer { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: DOFSynchronizer(const Mesh & mesh, UInt nb_degre_of_freedom); virtual ~DOFSynchronizer(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// init the scheme for scatter and gather operation, need extra memory void initScatterGatherCommunicationScheme(); /// initialize the equation number with local ids void initLocalDOFEquationNumbers(); /// initialize the equation number with local ids void initGlobalDOFEquationNumbers(); /** * Gather the DOF value on the root proccessor * * @param to_gather data to gather * @param root processor on which data are gathered * @param gathered Vector containing the gathered data, only valid on root processor */ template void gather(const Vector & to_gather, UInt root, Vector * gathered = NULL) const; /** * Scatter a DOF Vector form root to all processors * * @param scattered data to scatter, only valid on root processor * @param root processor scattering data * @param to_scatter result of scattered data */ template void scatter(Vector & scattered, UInt root, const Vector * to_scatter = NULL) const; template void synchronize(Vector & vector) const ; template void reduceSynchronize(Vector & vector) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the equation_number Vector AKANTU_GET_MACRO(LocalDOFEquationNumbers, local_dof_equation_numbers, const Vector &); AKANTU_GET_MACRO(GlobalDOFEquationNumbers, global_dof_equation_numbers, const Vector &); Vector * getLocalDOFEquationNumbersPointer(){return &local_dof_equation_numbers;}; + Vector * getGlobalDOFEquationNumbersPointer(){return &global_dof_equation_numbers;}; typedef unordered_map::type GlobalEquationNumberMap; AKANTU_GET_MACRO(GlobalEquationNumberToLocal, global_dof_equation_number_to_local, const GlobalEquationNumberMap &) /// get the Vector of global ids of the dofs AKANTU_GET_MACRO(DOFGlobalIDs, dof_global_ids, const Vector &); /// get the global id of a dof __aka_inline__ UInt getDOFGlobalID(UInt local_id) const { return dof_global_ids(local_id); } /// get the local id of a global dof __aka_inline__ UInt getDOFLocalID(UInt global_id) const { return global_dof_to_local.find(global_id)->second; } /// get the DOF type Vector AKANTU_GET_MACRO(DOFTypes, dof_types, const Vector &); AKANTU_GET_MACRO(NbDOFs, nb_dofs, UInt); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// equation number position where a dof is synchronized in the matrix (by default = global id) Vector global_dof_equation_numbers; Vector local_dof_equation_numbers; GlobalEquationNumberMap global_dof_equation_number_to_local; /// DOF global id Vector dof_global_ids; /* * DOF type -3 pure ghost, -2 master for the dof, -1 normal dof, i in * [0-N] slave dof and master is proc i */ Vector dof_types; /// number of dofs UInt nb_dofs; UInt nb_global_dofs; unordered_map::type global_dof_to_local; UInt prank; UInt psize; StaticCommunicator * communicator; struct PerProcInformations { /// dofs to send to the proc Vector slave_dofs; /// dofs to recvs from the proc Vector master_dofs; /* ---------------------------------------------------------------------- */ /* Data for gather/scatter */ /* ---------------------------------------------------------------------- */ /// the dof that the node handle Vector dofs; /// the dof that the proc need Vector needed_dofs; }; std::vector proc_informations; /// nb dofs with type -1 or -2 UInt nb_local_dofs; /// nb dof with type >= 0 UInt nb_needed_dofs; bool gather_scatter_scheme_initialized; }; /* -------------------------------------------------------------------------- */ /* __aka_inline__ functions */ /* -------------------------------------------------------------------------- */ #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "dof_synchronizer_inline_impl.cc" #endif /// standard output stream operator // __aka_inline__ std::ostream & operator <<(std::ostream & stream, const DOFSynchronizer & _this) // { // _this.printself(stream); // return stream; // } __END_AKANTU__ #endif /* __AKANTU_DOF_SYNCHRONIZER_HH__ */ 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 5cdb31326..becb207de 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,93 +1,94 @@ /** * @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.hh" #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); SolidMechanicsModel * model = new SolidMechanicsModel(mesh); /* -------------------------------------------------------------------------- */ model->initVectors(); model->getForce().clear(); model->getVelocity().clear(); model->getAcceleration().clear(); model->getDisplacement().clear(); /* ------------------------------------------------------------------------ */ model->initExplicit(); model->initModel(); model->readMaterials("material.dat"); model->initMaterials(); /* -------------------------------------------------------------------------- */ - model->initPBC(1,1,1); + model->setPBC(1,1,1); + model->initPBC(); model->assembleMassLumped(); /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER iohelper::DumperParaview dumper; dumper.SetMode(iohelper::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, iohelper::HEX1, mesh.getNbElement(_hexahedron_8), iohelper::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_solid_mechanics_model/patch_tests/explicit/patch_test_explicit.cc b/test/test_model/test_solid_mechanics_model/patch_tests/explicit/patch_test_explicit.cc index 8f3e6a747..4de149ea8 100644 --- a/test/test_model/test_solid_mechanics_model/patch_tests/explicit/patch_test_explicit.cc +++ b/test/test_model/test_solid_mechanics_model/patch_tests/explicit/patch_test_explicit.cc @@ -1,379 +1,342 @@ /** * @file test_check_stress.cc * @author David Kammer * @date Wed Feb 16 13:56:42 2011 * * @brief patch test for elastic material in solid mechanics model * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #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 "material.hh" #include "element_class.hh" #ifdef AKANTU_USE_IOHELPER # include "io_helper.hh" #endif //AKANTU_USE_IOHELPER using namespace akantu; Real alpha [3][4] = { { 0.01, 0.02, 0.03, 0.04 }, { 0.05, 0.06, 0.07, 0.08 }, { 0.09, 0.10, 0.11, 0.12 } }; #ifdef AKANTU_USE_IOHELPER static void paraviewInit(iohelper::Dumper & dumper, const SolidMechanicsModel & model); static void paraviewDump(iohelper::Dumper & dumper); #endif /* -------------------------------------------------------------------------- */ template static types::Matrix prescribed_strain() { UInt spatial_dimension = ElementClass::getSpatialDimension(); types::Matrix strain(spatial_dimension, spatial_dimension); for (UInt i = 0; i < spatial_dimension; ++i) { for (UInt j = 0; j < spatial_dimension; ++j) { strain(i, j) = alpha[i][j + 1]; } } return strain; } template static types::Matrix prescribed_stress() { UInt spatial_dimension = ElementClass::getSpatialDimension(); types::Matrix stress(spatial_dimension, spatial_dimension); //plane strain in 2d types::Matrix strain(spatial_dimension, spatial_dimension); types::Matrix pstrain; pstrain = prescribed_strain(); Real nu = 0.3; Real E = 2.1e11; Real trace = 0; /// symetric part of the strain tensor for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) strain(i,j) = 0.5 * (pstrain(i, j) + pstrain(j, i)); for (UInt i = 0; i < spatial_dimension; ++i) trace += strain(i,i); Real Ep = E / (1 + nu); for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) { stress(i, j) = Ep * strain(i,j); if(i == j) stress(i, j) += Ep * (nu / (1 - 2*nu)) * trace; } return stress; } /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { initialize(argc, argv); debug::setDebugLevel(dblWarning); UInt dim = ElementClass::getSpatialDimension(); const ElementType element_type = TYPE; - // UInt imposing_steps = 1000; - // Real max_displacement = -0.01; - UInt damping_steps = 400000; UInt damping_interval = 50; Real damping_ratio = 0.99; UInt additional_steps = 20000; UInt max_steps = damping_steps + additional_steps; /// load mesh Mesh my_mesh(dim); MeshIOMSH mesh_io; std::stringstream filename; filename << TYPE << ".msh"; mesh_io.read(filename.str(), my_mesh); UInt nb_nodes = my_mesh.getNbNodes(); /// declaration of model SolidMechanicsModel my_model(my_mesh); /// model initialization my_model.initVectors(); // initialize the vectors my_model.getForce().clear(); my_model.getVelocity().clear(); my_model.getAcceleration().clear(); my_model.getDisplacement().clear(); my_model.initExplicit(); my_model.initModel(); my_model.readMaterials("material_check_stress.dat"); my_model.initMaterials(); Real time_step = my_model.getStableTimeStep()/5.; my_model.setTimeStep(time_step); my_model.assembleMassLumped(); std::cout << "The number of time steps is: " << max_steps << " (" << time_step << "s)" << std::endl; - // // boundary conditions - // Vector top_nodes(0, 1); - // Vector base_nodes(0, 1); - + // boundary conditions const Vector & coordinates = my_mesh.getNodes(); Vector & displacement = my_model.getDisplacement(); Vector & boundary = my_model.getBoundary(); MeshUtils::buildFacets(my_mesh); MeshUtils::buildSurfaceID(my_mesh); CSR surface_nodes; MeshUtils::buildNodesPerSurface(my_mesh, surface_nodes); - - //CSR::iterator snode = surface_nodes.begin(0); - for (UInt s = 0; s < surface_nodes.getNbRows(); ++s) { CSR::iterator snode = surface_nodes.begin(s); for(; snode != surface_nodes.end(s); ++snode) { UInt n = *snode; std::cout << "Node " << n << std::endl; for (UInt i = 0; i < dim; ++i) { displacement(n, i) = alpha[i][0]; for (UInt j = 0; j < dim; ++j) { displacement(n, i) += alpha[i][j + 1] * coordinates(n, j); } boundary(n, i) = true; } - // if (coordinates(n, 0) < Math::tolerance) { - // boundary(n, 0) = true; - // displacement(n, 0) = 0.; - // } - // if (coordinates(n, 1) < Math::tolerance) { - // boundary(n, 1) = true; - // displacement(n, 1) = 0.; - // } } } Vector & velocity = my_model.getVelocity(); #ifdef AKANTU_USE_IOHELPER my_model.updateResidual(); iohelper::DumperParaview dumper; paraviewInit(dumper, my_model); #endif std::ofstream energy; std::stringstream energy_filename; energy_filename << "energy_" << TYPE << ".csv"; energy.open(energy_filename.str().c_str()); energy << "id,time,ekin" << std::endl; Real ekin_mean = 0.; - - double * disp = displacement.values; - double * vel = velocity.values; - - /* ------------------------------------------------------------------------ */ /* Main loop */ /* ------------------------------------------------------------------------ */ for(UInt s = 1; s <= max_steps; ++s) { if(s % 10000 == 0) std::cout << "passing step " << s << "/" << max_steps << " (" << s*time_step << "s)" <(imposing_steps)) * s; - // for(UInt n = 0; n < top_nodes.getSize(); ++n) { - // UInt node = top_nodes(n); - // displacement(node, 1) = current_displacement; - // } - // } - // damp velocity in order to find equilibrium if((s < damping_steps) && (s % damping_interval == 0)) { velocity *= damping_ratio; } if(s % 1000 == 0) { ekin_mean = ekin_mean / 1000.; std::cout << "Ekin mean = " << ekin_mean << std::endl; if (ekin_mean < 1e-10) break; ekin_mean = 0.; } my_model.explicitPred(); my_model.updateResidual(); my_model.updateAcceleration(); my_model.explicitCorr(); akantu::Real ekin = my_model.getKineticEnergy();; ekin_mean += ekin; if(s % 1000 == 0) energy << s << "," << s*time_step << "," << ekin << std::endl; #ifdef AKANTU_USE_IOHELPER if(s % 10000 == 0) paraviewDump(dumper); #endif } #ifdef AKANTU_USE_IOHELPER paraviewDump(dumper); #endif energy.close(); - // UInt check_element = 0; - // UInt quadrature_point = 0; UInt nb_quadrature_points = my_model.getFEM().getNbQuadraturePoints(TYPE); Vector & stress_vect = const_cast &>(my_model.getMaterial(0).getStress(element_type)); Vector & strain_vect = const_cast &>(my_model.getMaterial(0).getStrain(element_type)); Vector::iterator stress_it = stress_vect.begin(dim, dim); Vector::iterator strain_it = strain_vect.begin(dim, dim); types::Matrix presc_stress; presc_stress = prescribed_stress(); types::Matrix presc_strain; presc_strain = prescribed_strain(); UInt nb_element = my_mesh.getNbElement(TYPE); for (UInt el = 0; el < nb_element; ++el) { for (UInt q = 0; q < nb_quadrature_points; ++q) { types::Matrix & stress = *stress_it; types::Matrix & strain = *strain_it; for (UInt i = 0; i < dim; ++i) { for (UInt j = 0; j < dim; ++j) { if(!(std::abs(strain(i, j) - presc_strain(i, j)) < 1e-9)) { std::cerr << "strain[" << i << "," << j << "] = " << strain(i, j) << " but should be = " << presc_strain(i, j) << " (-" << std::abs(strain(i, j) - presc_strain(i, j)) << ") [el : " << el<< " - q : " << q << "]" << std::endl; std::cerr << strain << presc_strain << std::endl; return EXIT_FAILURE; } if(!(std::abs(stress(i, j) - presc_stress(i, j)) < 1e2)) { std::cerr << "stress[" << i << "," << j << "] = " << stress(i, j) << " but should be = " << presc_stress(i, j) << " (-" << std::abs(stress(i, j) - presc_stress(i, j)) << ") [el : " << el<< " - q : " << q << "]" << std::endl; std::cerr << stress << presc_stress << std::endl; return EXIT_FAILURE; } } } ++stress_it; ++strain_it; } } for (UInt n = 0; n < nb_nodes; ++n) { for (UInt i = 0; i < dim; ++i) { Real disp = alpha[i][0]; for (UInt j = 0; j < dim; ++j) { disp += alpha[i][j + 1] * coordinates(n, j); } if(!(std::abs(displacement(n,i) - disp) < 1e-7)) { std::cerr << "displacement(" << n << ", " << i <<")=" << displacement(n,i) << " should be equal to " << disp << "(" << displacement(n,i) - disp << ")" << std::endl; return EXIT_FAILURE; } } } - - // std::cout << "Strain : " << strain; - // std::cout << "Stress : " << stress; - finalize(); return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ /* iohelper::Dumper vars */ /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER template static iohelper::ElemType paraviewType(); template <> iohelper::ElemType paraviewType<_segment_2>() { return iohelper::LINE1; } template <> iohelper::ElemType paraviewType<_segment_3>() { return iohelper::LINE2; } template <> iohelper::ElemType paraviewType<_triangle_3>() { return iohelper::TRIANGLE1; } template <> iohelper::ElemType paraviewType<_triangle_6>() { return iohelper::TRIANGLE2; } template <> iohelper::ElemType paraviewType<_quadrangle_4>() { return iohelper::QUAD1; } template <> iohelper::ElemType paraviewType<_quadrangle_8>() { return iohelper::QUAD2; } template <> iohelper::ElemType paraviewType<_tetrahedron_4>() { return iohelper::TETRA1; } template <> iohelper::ElemType paraviewType<_tetrahedron_10>() { return iohelper::TETRA2; } template <> iohelper::ElemType paraviewType<_hexahedron_8>() { return iohelper::HEX1; } void paraviewInit(iohelper::Dumper & dumper, const SolidMechanicsModel & model) { UInt spatial_dimension = ElementClass::getSpatialDimension(); UInt nb_nodes = model.getFEM().getMesh().getNbNodes(); UInt nb_element = model.getFEM().getMesh().getNbElement(TYPE); std::stringstream filename; filename << "out_" << TYPE; dumper.SetMode(iohelper::TEXT); dumper.SetPoints(model.getFEM().getMesh().getNodes().values, spatial_dimension, nb_nodes, filename.str().c_str()); dumper.SetConnectivity((int *)model.getFEM().getMesh().getConnectivity(TYPE).values, paraviewType(), nb_element, iohelper::C_MODE); dumper.AddNodeDataField(model.getDisplacement().values, spatial_dimension, "displacements"); dumper.AddNodeDataField(model.getVelocity().values, spatial_dimension, "velocity"); dumper.AddNodeDataField(model.getResidual().values, spatial_dimension, "force"); dumper.AddNodeDataField(model.getMass().values, spatial_dimension, "mass"); dumper.AddNodeDataField(model.getForce().values, spatial_dimension, "applied_force"); dumper.AddElemDataField(model.getMaterial(0).getStrain(TYPE).values, spatial_dimension*spatial_dimension, "strain"); dumper.AddElemDataField(model.getMaterial(0).getStrain(TYPE).values, spatial_dimension*spatial_dimension, "stress"); dumper.SetEmbeddedValue("displacements", 1); dumper.SetEmbeddedValue("applied_force", 1); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); } /* -------------------------------------------------------------------------- */ void paraviewDump(iohelper::Dumper & dumper) { dumper.Dump(); } #endif diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_bar_traction2d_structured_pbc.cc b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_bar_traction2d_structured_pbc.cc index 2cad68063..c522b34ae 100644 --- a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_bar_traction2d_structured_pbc.cc +++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_bar_traction2d_structured_pbc.cc @@ -1,194 +1,195 @@ /** * @file test_solid_mechanics_model_bar_traction2d_structured_pbc.cc * @author David Kammer * @date Thu Dec 22 13:37:06 2011 * * @brief test of pbc class for SolidMechanicsModel * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "solid_mechanics_model.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER # include "io_helper.hh" #endif //AKANTU_USE_IOHELPER #ifdef AKANTU_USE_IOHELPER akantu::ElementType type = akantu::_quadrangle_4; iohelper::ElemType paraview_type = iohelper::QUAD1; #endif //AKANTU_USE_IOHELPER #ifdef AKANTU_USE_IOHELPER static void paraviewInit(iohelper::Dumper & dumper, const akantu::SolidMechanicsModel & model); static void paraviewDump(iohelper::Dumper & dumper); #endif int main(int argc, char *argv[]) { akantu::debug::setDebugLevel(akantu::dblWarning); akantu::initialize(argc, argv); akantu::UInt spatial_dimension = 2; akantu::UInt max_steps = 1000; akantu::Real time_factor = 0.2; akantu::Real epot, ekin; akantu::Mesh mesh(spatial_dimension); akantu::MeshIOMSH mesh_io; mesh_io.read("bar_structured1.msh", mesh); akantu::SolidMechanicsModel * model = new akantu::SolidMechanicsModel(mesh); /// model initialization model->initVectors(); /// set vectors to 0 akantu::UInt nb_nodes = model->getFEM().getMesh().getNbNodes(); memset(model->getForce().values, 0, spatial_dimension*nb_nodes*sizeof(akantu::Real)); memset(model->getVelocity().values, 0, spatial_dimension*nb_nodes*sizeof(akantu::Real)); memset(model->getAcceleration().values, 0, spatial_dimension*nb_nodes*sizeof(akantu::Real)); memset(model->getDisplacement().values, 0, spatial_dimension*nb_nodes*sizeof(akantu::Real)); model->initExplicit(); model->initModel(); model->readMaterials("material.dat"); model->initMaterials(); std::cout << model->getMaterial(0) << std::endl; - model->initPBC(1,0,0); + model->setPBC(1,0,0); + model->initPBC(); model->assembleMassLumped(); /// boundary conditions mesh.computeBoundingBox(); akantu::Real eps = 1e-16; akantu::Real signal_start = 0.6*mesh.getXMax(); akantu::Real signal_end = 0.7*mesh.getXMax(); akantu::Real delta_d = signal_end - signal_start; akantu::Real signal = 1.; const akantu::Vector & coords = model->getFEM().getMesh().getNodes(); akantu::Vector & disp = model->getDisplacement(); for (akantu::UInt i = 0; i < nb_nodes; ++i) { if(coords(i,0) >= signal_start && coords(i,0) <= signal_end) { if (coords(i,0) <= 0.5 * (signal_start + signal_end)) disp(i,0) = (coords(i,0) - signal_start) * 2 * signal / delta_d; else disp(i,0) = (signal_end - coords(i,0)) * 2 * signal / delta_d; } if(coords(i,1) <= eps || coords(i,1) >= 1 - eps ) { model->getBoundary().values[spatial_dimension*i + 1] = true; } } akantu::Real time_step = model->getStableTimeStep() * time_factor; std::cout << "Time Step = " << time_step << "s" << std::endl; model->setTimeStep(time_step); std::ofstream energy; energy.open("energy_2d_pbc.csv"); energy << "id,epot,ekin,tot" << std::endl; #ifdef AKANTU_USE_IOHELPER /// initialize the paraview output model->updateResidual(); iohelper::DumperParaview dumper; paraviewInit(dumper, *model); #endif //AKANTU_USE_IOHELPER for(akantu::UInt s = 1; s <= max_steps; ++s) { model->explicitPred(); model->updateResidual(); model->updateAcceleration(); model->explicitCorr(); epot = model->getPotentialEnergy(); ekin = model->getKineticEnergy(); #ifdef AKANTU_USE_IOHELPER if(s % 20 == 0) paraviewDump(dumper); #endif //AKANTU_USE_IOHELPER std::cerr << "passing step " << s << "/" << max_steps << std::endl; energy << s << "," << epot << "," << ekin << "," << epot + ekin << std::endl; } energy.close(); return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ /* iohelper::Dumper vars */ /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER void paraviewInit(iohelper::Dumper & dumper, const akantu::SolidMechanicsModel & model) { akantu::Real spatial_dimension = model.getSpatialDimension(); akantu::UInt nb_nodes = model.getFEM().getMesh().getNbNodes(); akantu::UInt nb_element = model.getFEM().getMesh().getNbElement(type); dumper.SetPoints(model.getFEM().getMesh().getNodes().values, spatial_dimension, nb_nodes, "bar_traction_2d_structured_pbc"); dumper.SetConnectivity((int *)model.getFEM().getMesh().getConnectivity(type).values, paraview_type, nb_element, iohelper::C_MODE); dumper.AddNodeDataField(model.getDisplacement().values, spatial_dimension, "displacements"); dumper.AddNodeDataField(model.getVelocity().values, spatial_dimension, "velocity"); dumper.AddNodeDataField(model.getAcceleration().values, spatial_dimension, "acceleration"); dumper.AddNodeDataField(model.getResidual().values, spatial_dimension, "force"); dumper.AddNodeDataField(model.getMass().values, spatial_dimension, "mass"); dumper.AddNodeDataField(model.getForce().values, spatial_dimension, "applied_force"); dumper.AddElemDataField(model.getMaterial(0).getStrain(type).values, spatial_dimension*spatial_dimension, "strain"); dumper.AddElemDataField(model.getMaterial(0).getStress(type).values, spatial_dimension*spatial_dimension, "stress"); dumper.SetEmbeddedValue("displacements", 1); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); } /* -------------------------------------------------------------------------- */ void paraviewDump(iohelper::Dumper & dumper) { dumper.Dump(); } #endif diff --git a/test/test_synchronizer/test_grid_synchronizer.cc b/test/test_synchronizer/test_grid_synchronizer.cc index ae50018b0..15b61901a 100644 --- a/test/test_synchronizer/test_grid_synchronizer.cc +++ b/test/test_synchronizer/test_grid_synchronizer.cc @@ -1,331 +1,330 @@ /** * @file test_grid_synchronizer.cc * @author Nicolas Richart * @date Mon Nov 7 11:58:02 2011 * * @brief test the GridSynchronizer object * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_grid.hh" #include "mesh.hh" #include "mesh_io.hh" #include "grid_synchronizer.hh" #include "mesh_partition.hh" #include "synchronizer_registry.hh" #ifdef AKANTU_USE_IOHELPER # include "io_helper.hh" -using namespace iohelper; #endif //AKANTU_USE_IOHELPER using namespace akantu; class TestAccessor : public DataAccessor { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: TestAccessor(const Mesh & mesh); ~TestAccessor(); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Barycenter, barycenter, Real); /* ------------------------------------------------------------------------ */ /* Ghost Synchronizer inherited members */ /* ------------------------------------------------------------------------ */ protected: virtual UInt getNbDataToPack(const Element & element, SynchronizationTag tag) const; virtual UInt getNbDataToUnpack(const Element & element, SynchronizationTag tag) const; virtual void packData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag) const; virtual void unpackData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: std::string id; ByElementTypeReal barycenter; const Mesh & mesh; }; /* -------------------------------------------------------------------------- */ /* TestSynchronizer implementation */ /* -------------------------------------------------------------------------- */ TestAccessor::TestAccessor(const Mesh & mesh) : barycenter("barycenter", id), mesh(mesh) { UInt spatial_dimension = mesh.getSpatialDimension(); id = "test_synchronizer"; Mesh::ConnectivityTypeList::const_iterator it; const Mesh::ConnectivityTypeList & ghost_type_list = mesh.getConnectivityTypeList(_ghost); for(it = ghost_type_list.begin(); it != ghost_type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; UInt nb_ghost_element = mesh.getNbElement(*it,_ghost); barycenter.alloc(nb_ghost_element, spatial_dimension, *it); } } TestAccessor::~TestAccessor() { } UInt TestAccessor::getNbDataToPack(const Element & element, __attribute__ ((unused)) SynchronizationTag tag) const { return Mesh::getSpatialDimension(element.type) * sizeof(Real); } UInt TestAccessor::getNbDataToUnpack(const Element & element, __attribute__ ((unused)) SynchronizationTag tag) const { return Mesh::getSpatialDimension(element.type) * sizeof(Real); } void TestAccessor::packData(CommunicationBuffer & buffer, const Element & element, __attribute__ ((unused)) SynchronizationTag tag) const { UInt spatial_dimension = Mesh::getSpatialDimension(element.type); types::RVector bary(spatial_dimension); mesh.getBarycenter(element.element, element.type, bary.storage()); buffer << bary; } void TestAccessor::unpackData(CommunicationBuffer & buffer, const Element & element, __attribute__ ((unused)) SynchronizationTag tag) { UInt spatial_dimension = Mesh::getSpatialDimension(element.type); Vector::iterator bary = barycenter(element.type).begin(spatial_dimension); buffer >> bary[element.element]; } /* -------------------------------------------------------------------------- */ /* Main */ /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { akantu::initialize(argc, argv); UInt spatial_dimension = 2; ElementType type = _triangle_3; Mesh mesh(spatial_dimension); StaticCommunicator * comm = StaticCommunicator::getStaticCommunicator(); Int psize = comm->getNbProc(); Int prank = comm->whoAmI(); DistributedSynchronizer * communicator; if(prank == 0) { MeshIOMSH mesh_io; mesh_io.read("triangle.msh", mesh); MeshPartition * partition = new MeshPartitionScotch(mesh, spatial_dimension); partition->partitionate(psize); communicator = DistributedSynchronizer::createDistributedSynchronizerMesh(mesh, partition); delete partition; } else { communicator = DistributedSynchronizer::createDistributedSynchronizerMesh(mesh, NULL); } #ifdef AKANTU_USE_IOHELPER double * part; unsigned int nb_nodes = mesh.getNbNodes(); unsigned int nb_element = mesh.getNbElement(type); - DumperParaview dumper; - dumper.SetMode(TEXT); + iohelper::DumperParaview dumper; + dumper.SetMode(iohelper::TEXT); dumper.SetParallelContext(prank, psize); dumper.SetPoints(mesh.getNodes().values, spatial_dimension, nb_nodes, "test-grid-synchronizer"); dumper.SetConnectivity((int*) mesh.getConnectivity(type).values, - TRIANGLE1, nb_element, C_MODE); + iohelper::TRIANGLE1, nb_element, iohelper::C_MODE); part = new double[nb_element]; for (unsigned int i = 0; i < nb_element; ++i) part[i] = prank; dumper.AddElemDataField(part, 1, "partitions"); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); delete [] part; unsigned int nb_ghost_element = mesh.getNbElement(type,_ghost); - DumperParaview dumper_ghost; - dumper_ghost.SetMode(TEXT); + iohelper::DumperParaview dumper_ghost; + dumper_ghost.SetMode(iohelper::TEXT); dumper_ghost.SetParallelContext(prank, psize); dumper_ghost.SetPoints(mesh.getNodes().values, spatial_dimension, nb_nodes, "test-grid-synchronizer-ghost"); dumper_ghost.SetConnectivity((int*) mesh.getConnectivity(type,_ghost).values, - TRIANGLE1, nb_ghost_element, C_MODE); + iohelper::TRIANGLE1, nb_ghost_element, iohelper::C_MODE); part = new double[nb_ghost_element]; for (unsigned int i = 0; i < nb_ghost_element; ++i) part[i] = prank; dumper_ghost.AddElemDataField(part, 1, "partitions"); dumper_ghost.SetPrefix("paraview/"); dumper_ghost.Init(); dumper_ghost.Dump(); delete [] part; #endif comm->barrier(); mesh.computeBoundingBox(); Real lower_bounds[spatial_dimension]; Real upper_bounds[spatial_dimension]; mesh.getLowerBounds(lower_bounds); mesh.getUpperBounds(upper_bounds); Real spacing[spatial_dimension]; for (UInt i = 0; i < spatial_dimension; ++i) { spacing[i] = (upper_bounds[i] - lower_bounds[i]) / 10.; } Real local_lower_bounds[spatial_dimension]; Real local_upper_bounds[spatial_dimension]; mesh.getLocalLowerBounds(local_lower_bounds); mesh.getLocalUpperBounds(local_upper_bounds); RegularGrid grid(spatial_dimension, local_lower_bounds, local_upper_bounds, spacing); GhostType ghost_type = _not_ghost; Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type); std::cout << prank << " TTTTTTTOOOOOOOOOOOOOTTTTTTTTTTTTTTTTOOOOOOOOOOOOO" << grid < & barycenter = barycenters(*it); barycenter.resize(nb_element); Vector::iterator bary_it = barycenter.begin(spatial_dimension); for (UInt elem = 0; elem < nb_element; ++elem) { mesh.getBarycenter(elem, *it, bary_it->storage()); grid.count(*bary_it); ++bary_it; } } Element e; e.ghost_type = ghost_type; std::cout << prank << " TTTTTTTOOOOOOOOOOOOOTTTTTTTTTTTTTTTTOOOOOOOOOOOOO" << std::endl; // second insert the point in the cells grid.beginInsertions(); it = mesh.firstType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { UInt nb_element = mesh.getNbElement(*it); e.type = *it; Vector & barycenter = barycenters(*it); Vector::iterator bary_it = barycenter.begin(spatial_dimension); for (UInt elem = 0; elem < nb_element; ++elem) { e.element = elem; grid.insert(e, *bary_it); ++bary_it; } } grid.endInsertions(); MeshIOMSH mesh_io; std::stringstream sstr; sstr << "mesh_" << prank << ".msh"; mesh_io.write(sstr.str(), mesh); Mesh grid_mesh(spatial_dimension, "grid_mesh", 0); std::stringstream sstr_grid; sstr_grid << "grid_mesh_" << prank << ".msh"; grid.saveAsMesh(grid_mesh); mesh_io.write(sstr_grid.str(), grid_mesh); std::cout << "TTTTTTTOOOOOOOOOOOOOTTTTTTTTTTTTTTTTOOOOOOOOOOOOO" << std::endl; GridSynchronizer * grid_communicator = GridSynchronizer::createGridSynchronizer(mesh, grid); AKANTU_DEBUG_INFO("Creating TestAccessor"); TestAccessor test_accessor(mesh); SynchronizerRegistry synch_registry(test_accessor); synch_registry.registerSynchronizer(*grid_communicator, _gst_test); AKANTU_DEBUG_INFO("Synchronizing tag"); synch_registry.synchronize(_gst_test); #ifdef AKANTU_USE_IOHELPER nb_ghost_element = mesh.getNbElement(type,_ghost); nb_nodes = mesh.getNbNodes(); - dumper_ghost.SetMode(TEXT); + dumper_ghost.SetMode(iohelper::TEXT); dumper_ghost.SetParallelContext(prank, psize); dumper_ghost.SetPoints(mesh.getNodes().values, spatial_dimension, nb_nodes, "test-grid-synchronizer-ghost"); dumper_ghost.SetConnectivity((int*) mesh.getConnectivity(type,_ghost).values, - TRIANGLE1, nb_ghost_element, C_MODE); + iohelper::TRIANGLE1, nb_ghost_element, iohelper::C_MODE); part = new double[nb_ghost_element]; for (unsigned int i = 0; i < nb_ghost_element; ++i) part[i] = prank; dumper_ghost.AddElemDataField(part, 1, "partitions"); dumper_ghost.Dump(); delete [] part; unsigned int nb_grid_element = grid_mesh.getNbElement(_quadrangle_4); unsigned int nb_quadrature_points = 4; - DumperParaview dumper_grid; - dumper_grid.SetMode(TEXT); + iohelper::DumperParaview dumper_grid; + dumper_grid.SetMode(iohelper::TEXT); dumper_grid.SetParallelContext(prank, psize); dumper_grid.SetPoints(grid_mesh.getNodes().values, spatial_dimension, grid_mesh.getNbNodes(), "test-grid-synchronizer-grid"); dumper_grid.SetConnectivity((int*) grid_mesh.getConnectivity(_quadrangle_4).values, - QUAD1, nb_grid_element, C_MODE); + iohelper::QUAD1, nb_grid_element, iohelper::C_MODE); part = new double[nb_grid_element * nb_quadrature_points]; std::fill_n(part, nb_grid_element * nb_quadrature_points, prank); dumper_grid.AddElemDataField(part, 1, "partitions"); dumper_grid.SetPrefix("paraview/"); dumper_grid.Init(); dumper_grid.Dump(); delete [] part; #endif //AKANTU_USE_IOHELPER akantu::finalize(); return EXIT_SUCCESS; }