diff --git a/extra_packages/parallel-cohesive-element/package.cmake b/extra_packages/parallel-cohesive-element/package.cmake new file mode 100644 index 000000000..dd9d24a1e --- /dev/null +++ b/extra_packages/parallel-cohesive-element/package.cmake @@ -0,0 +1,34 @@ +#=============================================================================== +# @file package.cmake +# +# @author Nicolas Richart +# +# +# @brief package description for parallel cohesive elements +# +# @section LICENSE +# +# Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) +# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) +# +#=============================================================================== + +package_declare(parallel_cohesive_element + DESCRIPTION "Use parallel cohesive_element package of Akantu" + DEPENDS cohesive_element parallel) + +package_declare_sources(parallel_cohesive_element + cohesive_element_inserter_parallel.cc + cohesive_element_inserter_inline_impl.cc + solid_mechanics_model_cohesive_parallel.hh + solid_mechanics_model_cohesive_parallel.cc + solid_mechanics_model_cohesive_parallel_inline_impl.cc + facet_synchronizer.cc + facet_synchronizer.hh + facet_synchronizer_inline_impl.cc + facet_stress_synchronizer.cc + facet_stress_synchronizer.hh + ) + +package_declare_documentation(parallel_cohesive_element + "This option activates the parallel cohesive elements' features of AKANTU.") diff --git a/extra_packages/parallel-cohesive-element/src/cohesive_element_inserter_inline_impl.cc b/extra_packages/parallel-cohesive-element/src/cohesive_element_inserter_inline_impl.cc new file mode 100644 index 000000000..4869c7a58 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/src/cohesive_element_inserter_inline_impl.cc @@ -0,0 +1,92 @@ +/** + * @file cohesive_element_inserter_inline_impl.cc + * + * @author Marco Vocialta + * + * + * @brief Cohesive element inserter inline functions + * + * @section LICENSE + * + * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + */ + +/* -------------------------------------------------------------------------- */ +inline UInt CohesiveElementInserter::getNbDataForElements(const Array & elements, + SynchronizationTag tag) const { + AKANTU_DEBUG_IN(); + + UInt size = 0; + + if (tag == _gst_ce_groups) { + size = elements.getSize() * (sizeof(bool) + sizeof(unsigned int)); + } + + AKANTU_DEBUG_OUT(); + return size; +} + +/* -------------------------------------------------------------------------- */ +inline void CohesiveElementInserter::packElementData(CommunicationBuffer & buffer, + const Array & elements, + SynchronizationTag tag) const { + AKANTU_DEBUG_IN(); + if (tag == _gst_ce_groups) + packUnpackGroupedInsertionData(buffer,elements); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +inline void CohesiveElementInserter::unpackElementData(CommunicationBuffer & buffer, + const Array & elements, + SynchronizationTag tag) { + AKANTU_DEBUG_IN(); + if (tag == _gst_ce_groups) + packUnpackGroupedInsertionData(buffer,elements); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +template +inline void CohesiveElementInserter::packUnpackGroupedInsertionData(CommunicationBuffer & buffer, + const Array & elements) const { + + AKANTU_DEBUG_IN(); + + ElementType current_element_type = _not_defined; + GhostType current_ghost_type = _casper; + ElementTypeMapArray & physical_names = mesh_facets.registerData("physical_names"); + + Array *vect = NULL; + Array *vect2 = NULL; + + Array::const_iterator it = elements.begin(); + Array::const_iterator end = elements.end(); + for (; it != end; ++it) { + const Element & el = *it; + if(el.type != current_element_type || el.ghost_type != current_ghost_type) { + current_element_type = el.type; + current_ghost_type = el.ghost_type; + vect = &const_cast &>(insertion_facets(el.type, el.ghost_type)); + vect2 = &(physical_names(el.type, el.ghost_type)); + } + + Vector data(vect->storage() + el.element, 1); + Vector data2(vect2->storage() + el.element, 1); + + if(pack_mode) { + buffer << data; + buffer << data2; + } + else { + buffer >> data; + buffer >> data2; + } + } + + AKANTU_DEBUG_OUT(); +} diff --git a/extra_packages/parallel-cohesive-element/src/cohesive_element_inserter_parallel.cc b/extra_packages/parallel-cohesive-element/src/cohesive_element_inserter_parallel.cc new file mode 100644 index 000000000..8175da288 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/src/cohesive_element_inserter_parallel.cc @@ -0,0 +1,77 @@ +/** + * @file cohesive_element_inserter_parallel.cc + * + * @author Marco Vocialta + * + * + * @brief Parallel functions for the cohesive element inserter + * + * @section LICENSE + * + * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + */ + +/* -------------------------------------------------------------------------- */ +#include "cohesive_element_inserter.hh" + +__BEGIN_AKANTU__ + +/* -------------------------------------------------------------------------- */ +void CohesiveElementInserter::initParallel(FacetSynchronizer * facet_synchronizer, + DistributedSynchronizer * distributed_synchronizer) { + AKANTU_DEBUG_IN(); + + this->facet_synchronizer = facet_synchronizer; + + global_ids_updater = new GlobalIdsUpdater(mesh, distributed_synchronizer); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void CohesiveElementInserter::updateNodesType(Mesh & mesh, + NewNodesEvent & node_event) { + AKANTU_DEBUG_IN(); + + Array & doubled_nodes = node_event.getList(); + UInt local_nb_new_nodes = doubled_nodes.getSize(); + + Array & nodes_type = mesh.getNodesType(); + UInt nb_old_nodes = nodes_type.getSize(); + nodes_type.resize(nb_old_nodes + local_nb_new_nodes); + + for (UInt n = 0; n < local_nb_new_nodes; ++n) { + UInt old_node = doubled_nodes(n, 0); + UInt new_node = doubled_nodes(n, 1); + nodes_type(new_node) = nodes_type(old_node); + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +UInt CohesiveElementInserter::updateGlobalIDs(NewNodesEvent & node_event) { + AKANTU_DEBUG_IN(); + + Array & doubled_nodes = node_event.getList(); + + UInt total_nb_new_nodes + = global_ids_updater->updateGlobalIDsLocally(doubled_nodes.getSize()); + + AKANTU_DEBUG_OUT(); + return total_nb_new_nodes; +} + +void CohesiveElementInserter::synchronizeGlobalIDs(NewNodesEvent & node_event) { + AKANTU_DEBUG_IN(); + + global_ids_updater->synchronizeGlobalIDs(); + + AKANTU_DEBUG_OUT(); +} + + + +__END_AKANTU__ diff --git a/extra_packages/parallel-cohesive-element/src/facet_stress_synchronizer.cc b/extra_packages/parallel-cohesive-element/src/facet_stress_synchronizer.cc new file mode 100644 index 000000000..2943fbb1f --- /dev/null +++ b/extra_packages/parallel-cohesive-element/src/facet_stress_synchronizer.cc @@ -0,0 +1,127 @@ +/** + * @file facet_stress_synchronizer.cc + * + * @author Marco Vocialta + * + * + * @brief Stress check on facets synchronizer + * + * @section LICENSE + * + * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + */ + +/* -------------------------------------------------------------------------- */ + +#include "facet_stress_synchronizer.hh" + +/* -------------------------------------------------------------------------- */ + +__BEGIN_AKANTU__ + +/* -------------------------------------------------------------------------- */ +FacetStressSynchronizer::FacetStressSynchronizer(Mesh & mesh_facets, + SynchronizerID id, + MemoryID memory_id) : + DistributedSynchronizer(mesh_facets, id, memory_id) { + AKANTU_DEBUG_IN(); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +FacetStressSynchronizer * FacetStressSynchronizer:: +createFacetStressSynchronizer(FacetSynchronizer & facet_synchronizer, + Mesh & mesh_facets, + SynchronizerID id, + MemoryID memory_id) { + AKANTU_DEBUG_IN(); + + FacetStressSynchronizer & fs_synchronizer = + *(new FacetStressSynchronizer(mesh_facets, id, memory_id)); + + for (UInt p = 0; p < fs_synchronizer.nb_proc; ++p) { + fs_synchronizer.send_element[p].copy(facet_synchronizer.recv_element[p]); + fs_synchronizer.recv_element[p].copy(facet_synchronizer.send_element[p]); + } + + AKANTU_DEBUG_OUT(); + return &fs_synchronizer; +} + +/* -------------------------------------------------------------------------- */ +void FacetStressSynchronizer::updateFacetStressSynchronizer(const CohesiveElementInserter & inserter, + const ElementTypeMapArray & rank_to_element, + DataAccessor & data_accessor) { + AKANTU_DEBUG_IN(); + + updateElementList(send_element, inserter, rank_to_element); + updateElementList(recv_element, inserter, rank_to_element); + + std::map::iterator it = communications.begin(); + std::map::iterator end = communications.end(); + + for (; it != end; ++it) { + SynchronizationTag tag = it->first; + computeBufferSize(data_accessor, tag); + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void FacetStressSynchronizer::updateElementList(Array * elements, + const CohesiveElementInserter & inserter, + const ElementTypeMapArray & rank_to_element) { + AKANTU_DEBUG_IN(); + + for (UInt p = 0; p < nb_proc; ++p) { + + ElementType current_element_type = _not_defined; + GhostType current_ghost_type = _casper; + const Array * f_check = NULL; + const Array< std::vector > * element_to_facet = NULL; + + UInt nb_element = 0; + + Array::iterator it_new = elements[p].begin(); + Array::iterator it = elements[p].begin(); + Array::iterator end = elements[p].end(); + for (; it != end; ++it) { + const Element & el = *it; + + if(el.type != current_element_type || el.ghost_type != current_ghost_type) { + current_element_type = el.type; + current_ghost_type = el.ghost_type; + + element_to_facet = & mesh.getElementToSubelement(el.type, el.ghost_type); + + f_check = & inserter.getCheckFacets(el.type, el.ghost_type); + } + + if ( (*f_check)(el.element) ) { + + const Element & el_0 = (*element_to_facet)(el.element)[0]; + const Array & rank_el_0 = rank_to_element(el_0.type, el_0.ghost_type); + + const Element & el_1 = (*element_to_facet)(el.element)[1]; + const Array & rank_el_1 = rank_to_element(el_1.type, el_1.ghost_type); + + if ( (rank_el_0(el_0.element) == rank && rank_el_1(el_1.element) == p) || + (rank_el_0(el_0.element) == p && rank_el_1(el_1.element) == rank) ) { + *it_new = el; + ++it_new; + ++nb_element; + } + } + } + elements[p].resize(nb_element); + } + + AKANTU_DEBUG_OUT(); +} + + +__END_AKANTU__ diff --git a/extra_packages/parallel-cohesive-element/src/facet_stress_synchronizer.hh b/extra_packages/parallel-cohesive-element/src/facet_stress_synchronizer.hh new file mode 100644 index 000000000..1dcf5d3d7 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/src/facet_stress_synchronizer.hh @@ -0,0 +1,83 @@ +/** + * @file facet_stress_synchronizer.hh + * + * @author Marco Vocialta + * + * + * @brief Stress check on facets synchronizer + * + * @section LICENSE + * + * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + */ + +/* -------------------------------------------------------------------------- */ + +#ifndef __AKANTU_FACET_STRESS_SYNCHRONIZER_HH__ +#define __AKANTU_FACET_STRESS_SYNCHRONIZER_HH__ + +#include "distributed_synchronizer.hh" +#include "facet_synchronizer.hh" +#include "cohesive_element_inserter.hh" + +/* -------------------------------------------------------------------------- */ + +__BEGIN_AKANTU__ + +class FacetStressSynchronizer : public DistributedSynchronizer { + /* ------------------------------------------------------------------------ */ + /* Constructors/Destructors */ + /* ------------------------------------------------------------------------ */ +protected: + + FacetStressSynchronizer(Mesh & mesh_facets, + SynchronizerID id = "facet_stress_synchronizer", + MemoryID memory_id = 0); + + // virtual ~FacetStressSynchronizer(); + + /* ------------------------------------------------------------------------ */ + /* Methods */ + /* ------------------------------------------------------------------------ */ +public: + + static FacetStressSynchronizer * + createFacetStressSynchronizer(FacetSynchronizer & facet_synchronizer, + Mesh & mesh_facets, + SynchronizerID id = "facet_stress_synchronizer", + MemoryID memory_id = 0); + + void updateFacetStressSynchronizer(const CohesiveElementInserter & inserter, + const ElementTypeMapArray & rank_to_element, + DataAccessor & data_accessor); + +protected: + + void updateElementList(Array * elements, + const CohesiveElementInserter & inserter, + const ElementTypeMapArray & rank_to_element); + + /* ------------------------------------------------------------------------ */ + /* Accessors */ + /* ------------------------------------------------------------------------ */ +public: + + /* ------------------------------------------------------------------------ */ + /* Class Members */ + /* ------------------------------------------------------------------------ */ +private: + +}; + + +/* -------------------------------------------------------------------------- */ +/* inline functions */ +/* -------------------------------------------------------------------------- */ + +// #include "facet_stress_synchronizer_inline_impl.cc" + +__END_AKANTU__ + +#endif /* __AKANTU_FACET_STRESS_SYNCHRONIZER_HH__ */ diff --git a/extra_packages/parallel-cohesive-element/src/facet_synchronizer.cc b/extra_packages/parallel-cohesive-element/src/facet_synchronizer.cc new file mode 100644 index 000000000..867dbb5f1 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/src/facet_synchronizer.cc @@ -0,0 +1,425 @@ +/** + * @file facet_synchronizer.cc + * + * @author Marco Vocialta + * + * + * @brief Facet synchronizer for parallel simulations with cohesive elments + * + * @section LICENSE + * + * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + */ + +/* -------------------------------------------------------------------------- */ + +#include "facet_synchronizer.hh" + +/* -------------------------------------------------------------------------- */ + +__BEGIN_AKANTU__ + +/* -------------------------------------------------------------------------- */ +FacetSynchronizer::FacetSynchronizer(Mesh & mesh, + SynchronizerID id, + MemoryID memory_id) : + DistributedSynchronizer(mesh, id, memory_id) { + AKANTU_DEBUG_IN(); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +FacetSynchronizer * FacetSynchronizer:: +createFacetSynchronizer(DistributedSynchronizer & distributed_synchronizer, + Mesh & mesh, + SynchronizerID id, + MemoryID memory_id) { + AKANTU_DEBUG_IN(); + + FacetSynchronizer & f_synchronizer = *(new FacetSynchronizer(mesh, id, memory_id)); + + f_synchronizer.setupFacetSynchronization(distributed_synchronizer); + + AKANTU_DEBUG_OUT(); + return &f_synchronizer; +} + +/* -------------------------------------------------------------------------- */ +void FacetSynchronizer::updateDistributedSynchronizer(DistributedSynchronizer & distributed_synchronizer, + DataAccessor & data_accessor, + const Mesh & mesh_cohesive) { + AKANTU_DEBUG_IN(); + + Array * distrib_send_element = distributed_synchronizer.send_element; + Array * distrib_recv_element = distributed_synchronizer.recv_element; + + updateElementList(distrib_send_element, send_element, mesh_cohesive); + updateElementList(distrib_recv_element, recv_element, mesh_cohesive); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void FacetSynchronizer::updateElementList(Array * elements, + const Array * facets, + const Mesh & mesh_cohesive) { + AKANTU_DEBUG_IN(); + + ElementType current_element_type = _not_defined; + GhostType current_ghost_type = _casper; + + UInt old_nb_cohesive_elements = 0; + const Array< std::vector > * element_to_facet = NULL; + + for (UInt p = 0; p < nb_proc; ++p) { + const Array & fa = facets[p]; + Array & el = elements[p]; + + Array::const_iterator it = fa.begin(); + Array::const_iterator end = fa.end(); + + for (; it != end; ++it) { + const Element & facet = *it; + + if(facet.type != current_element_type || facet.ghost_type != current_ghost_type) { + current_element_type = facet.type; + current_ghost_type = facet.ghost_type; + + ElementType current_coh_element_type + = FEEngine::getCohesiveElementType(current_element_type); + + /// compute old number of cohesive elements + old_nb_cohesive_elements = mesh_cohesive.getNbElement(current_coh_element_type, + current_ghost_type); + old_nb_cohesive_elements -= mesh.getData("facet_to_double", + current_element_type, + current_ghost_type).getSize(); + + element_to_facet = & mesh.getData >("element_to_subelement", + current_element_type, + current_ghost_type); + } + + const Element & cohesive_element = (*element_to_facet)(facet.element)[1]; + + if (cohesive_element.kind == _ek_cohesive && + cohesive_element.element >= old_nb_cohesive_elements) + el.push_back(cohesive_element); + } + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void FacetSynchronizer::setupFacetSynchronization(DistributedSynchronizer & distributed_synchronizer) { + AKANTU_DEBUG_IN(); + + Array * distrib_send_element = distributed_synchronizer.send_element; + Array * distrib_recv_element = distributed_synchronizer.recv_element; + + /// build rank to facet correspondance + ElementTypeMapArray rank_to_facet("rank_to_facet", id); + initRankToFacet(rank_to_facet); + buildRankToFacet(rank_to_facet, distrib_recv_element); + + /// generate temp_send/recv element arrays with their connectivity + Array *> temp_send_element(nb_proc); + Array *> temp_recv_element(nb_proc); + Array *> send_connectivity(nb_proc); + Array *> recv_connectivity(nb_proc); + + UInt spatial_dimension = mesh.getSpatialDimension(); + + for (UInt p = 0; p < nb_proc; ++p) { + if (p == rank) continue; + std::stringstream sstr; sstr << p; + + temp_send_element(p) = + new ElementTypeMapArray("temp_send_element_proc_"+sstr.str(), id); + mesh.initElementTypeMapArray(*temp_send_element(p), 1, spatial_dimension - 1); + + temp_recv_element(p) = + new ElementTypeMapArray("temp_recv_element_proc_"+sstr.str(), id); + mesh.initElementTypeMapArray(*temp_recv_element(p), 1, spatial_dimension - 1); + + send_connectivity(p) = + new ElementTypeMapArray("send_connectivity_proc_"+sstr.str(), id); + mesh.initElementTypeMapArray(*send_connectivity(p), 1, spatial_dimension - 1, true); + + recv_connectivity(p) = + new ElementTypeMapArray("recv_connectivity_proc_"+sstr.str(), id); + mesh.initElementTypeMapArray(*recv_connectivity(p), 1, spatial_dimension - 1, true); + } + + /// build global connectivity arrays + getFacetGlobalConnectivity<_not_ghost>(distributed_synchronizer, + rank_to_facet, + distrib_send_element, + send_connectivity, + temp_send_element); + + getFacetGlobalConnectivity<_ghost>(distributed_synchronizer, + rank_to_facet, + distrib_recv_element, + recv_connectivity, + temp_recv_element); + + /// build send/recv facet arrays + buildSendElementList(send_connectivity, recv_connectivity, temp_send_element); + buildRecvElementList(temp_recv_element); + +#ifndef AKANTU_NDEBUG + /// count recv facets for each processor + Array nb_facets_recv(nb_proc); + nb_facets_recv.clear(); + + Mesh::type_iterator first = mesh.firstType(spatial_dimension - 1, _ghost); + Mesh::type_iterator last = mesh.lastType(spatial_dimension - 1, _ghost); + + for (; first != last; ++first) { + const Array & r_to_f = rank_to_facet(*first, _ghost); + UInt nb_facet = r_to_f.getSize(); + + for (UInt f = 0; f < nb_facet; ++f) { + UInt proc = r_to_f(f); + if (proc != rank) + ++nb_facets_recv(proc); + } + } + + for (UInt p = 0; p < nb_proc; ++p) { + AKANTU_DEBUG_ASSERT(nb_facets_recv(p) == recv_element[p].getSize(), + "Wrong number of recv facets"); + } + +#endif + + for (UInt p = 0; p < nb_proc; ++p) { + delete temp_send_element(p); + delete temp_recv_element(p); + delete send_connectivity(p); + delete recv_connectivity(p); + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void FacetSynchronizer::buildSendElementList(const Array *> & send_connectivity, + const Array *> & recv_connectivity, + const Array *> & temp_send_element) { + AKANTU_DEBUG_IN(); + + StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + + UInt spatial_dimension = mesh.getSpatialDimension(); + + GhostType ghost_type = _ghost; + + Mesh::type_iterator first = mesh.firstType(spatial_dimension - 1, ghost_type); + Mesh::type_iterator last = mesh.lastType(spatial_dimension - 1, ghost_type); + + /// do every communication by element type + for (; first != last; ++first) { + + ElementType facet_type = *first; + + std::vector send_requests; + UInt * send_size = new UInt[nb_proc]; + + /// send asynchronous data + for (UInt p = 0; p < nb_proc; ++p) { + if (p == rank) continue; + + const Array & recv_conn = (*recv_connectivity(p))(facet_type, _ghost); + send_size[p] = recv_conn.getSize(); + + /// send connectivity size + send_requests.push_back(comm.asyncSend(send_size + p, + 1, + p, + Tag::genTag(rank, p, 0))); + + /// send connectivity data + send_requests.push_back(comm.asyncSend(recv_conn.storage(), + recv_conn.getSize() * + recv_conn.getNbComponent(), + p, + Tag::genTag(rank, p, 1))); + } + + UInt * recv_size = new UInt[nb_proc]; + UInt nb_nodes_per_facet = Mesh::getNbNodesPerElement(facet_type); + + /// receive data + for (UInt p = 0; p < nb_proc; ++p) { + if (p == rank) continue; + + /// receive connectivity size + comm.receive(recv_size + p, 1, p, Tag::genTag(p, rank, 0)); + + Array conn_to_match(recv_size[p], nb_nodes_per_facet); + + /// receive connectivity + comm.receive(conn_to_match.storage(), + conn_to_match.getSize() * conn_to_match.getNbComponent(), + p, + Tag::genTag(p, rank, 1)); + + const Array & send_conn = + (*send_connectivity(p))(facet_type, _not_ghost); + const Array & list = (*temp_send_element(p))(facet_type, _not_ghost); + UInt nb_local_facets = send_conn.getSize(); + + AKANTU_DEBUG_ASSERT(nb_local_facets == list.getSize(), + "connectivity and facet list have different sizes"); + + Array checked(nb_local_facets); + checked.clear(); + + Element facet(facet_type, 0, _not_ghost, _ek_regular); + + Array::iterator > c_to_match_it = + conn_to_match.begin(nb_nodes_per_facet); + Array::iterator > c_to_match_end = + conn_to_match.end(nb_nodes_per_facet); + + /// for every sent facet of other processors, find the + /// corresponding one in the local send connectivity data in + /// order to build the send_element arrays + for (; c_to_match_it != c_to_match_end; ++c_to_match_it) { + + Array::const_iterator > c_local_it = + send_conn.begin(nb_nodes_per_facet); + Array::const_iterator > c_local_end = + send_conn.end(nb_nodes_per_facet); + + for (UInt f = 0; f < nb_local_facets; ++f, ++c_local_it) { + if (checked(f)) continue; + + if ( (*c_to_match_it) == (*c_local_it) ) { + checked(f) = true; + facet.element = list(f); + send_element[p].push_back(facet); + break; + } + } + AKANTU_DEBUG_ASSERT(c_local_it != c_local_end, "facet not found"); + } + } + + /// wait for all communications to be done and free the + /// communication request array + comm.waitAll(send_requests); + comm.freeCommunicationRequest(send_requests); + + delete [] send_size; + delete [] recv_size; + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void FacetSynchronizer::buildRecvElementList(const Array *> & temp_recv_element) { + AKANTU_DEBUG_IN(); + + UInt spatial_dimension = mesh.getSpatialDimension(); + + for (UInt p = 0; p < nb_proc; ++p) { + if (p == rank) continue; + + GhostType ghost_type = _ghost; + + Mesh::type_iterator first = mesh.firstType(spatial_dimension - 1, ghost_type); + Mesh::type_iterator last = mesh.lastType(spatial_dimension - 1, ghost_type); + + for (; first != last; ++first) { + ElementType facet_type = *first; + + const Array & list = (*temp_recv_element(p))(facet_type, ghost_type); + UInt nb_local_facets = list.getSize(); + + Element facet(facet_type, 0, ghost_type, _ek_regular); + + for (UInt f = 0; f < nb_local_facets; ++f) { + facet.element = list(f); + recv_element[p].push_back(facet); + } + } + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void FacetSynchronizer::initRankToFacet(ElementTypeMapArray & rank_to_facet) { + AKANTU_DEBUG_IN(); + + UInt spatial_dimension = mesh.getSpatialDimension(); + + mesh.initElementTypeMapArray(rank_to_facet, 1, spatial_dimension - 1); + + GhostType ghost_type = _ghost; + + Mesh::type_iterator first = mesh.firstType(spatial_dimension - 1, ghost_type); + Mesh::type_iterator last = mesh.lastType(spatial_dimension - 1, ghost_type); + + for (; first != last; ++first) { + ElementType type = *first; + UInt nb_facet = mesh.getNbElement(type, ghost_type); + + Array & rank_to_f = rank_to_facet(type, ghost_type); + rank_to_f.resize(nb_facet); + + for (UInt f = 0; f < nb_facet; ++f) + rank_to_f(f) = rank; + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void FacetSynchronizer::buildRankToFacet(ElementTypeMapArray & rank_to_facet, + const Array * elements) { + AKANTU_DEBUG_IN(); + + for (UInt p = 0; p < nb_proc; ++p) { + if (p == rank) continue; + const Array & elem = elements[p]; + UInt nb_element = elem.getSize(); + + for (UInt el = 0; el < nb_element; ++el) { + ElementType type = elem(el).type; + GhostType gt = elem(el).ghost_type; + UInt el_index = elem(el).element; + + const Array & facet_to_element + = mesh.getSubelementToElement(type, gt); + UInt nb_facets_per_element = Mesh::getNbFacetsPerElement(type); + ElementType facet_type = Mesh::getFacetType(type); + + for (UInt f = 0; f < nb_facets_per_element; ++f) { + const Element & facet = facet_to_element(el_index, f); + if (facet == ElementNull) continue; + UInt facet_index = facet.element; + GhostType facet_gt = facet.ghost_type; + + if (facet_gt == _not_ghost) continue; + + Array & t_to_f = rank_to_facet(facet_type, facet_gt); + if ((p < t_to_f(facet_index)) || (t_to_f(facet_index) == rank)) + t_to_f(facet_index) = p; + } + } + } + + AKANTU_DEBUG_OUT(); +} + + +__END_AKANTU__ diff --git a/extra_packages/parallel-cohesive-element/src/facet_synchronizer.hh b/extra_packages/parallel-cohesive-element/src/facet_synchronizer.hh new file mode 100644 index 000000000..02b4a7878 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/src/facet_synchronizer.hh @@ -0,0 +1,117 @@ + /** + * @file facet_synchronizer.hh + * + * @author Marco Vocialta + * + * + * @brief Facet synchronizer for parallel simulations with cohesive elments + * + * @section LICENSE + * + * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + */ + +/* -------------------------------------------------------------------------- */ + +#ifndef __AKANTU_FACET_SYNCHRONIZER_HH__ +#define __AKANTU_FACET_SYNCHRONIZER_HH__ + +#include "distributed_synchronizer.hh" +#include "fe_engine.hh" + +/* -------------------------------------------------------------------------- */ + +__BEGIN_AKANTU__ + +class FacetSynchronizer : public DistributedSynchronizer { + /* ------------------------------------------------------------------------ */ + /* Constructors/Destructors */ + /* ------------------------------------------------------------------------ */ +protected: + + FacetSynchronizer(Mesh & mesh, + SynchronizerID id = "facet_synchronizer", + MemoryID memory_id = 0); + +// public: + +// virtual ~FacetSynchronizer() {}; + + /* ------------------------------------------------------------------------ */ + /* Methods */ + /* ------------------------------------------------------------------------ */ +public: + + /// get a distributed synchronizer and a facet mesh, and create the + /// associated FacetSynchronizer + static FacetSynchronizer * + createFacetSynchronizer(DistributedSynchronizer & distributed_synchronizer, + Mesh & mesh, + SynchronizerID id = "facet_synchronizer", + MemoryID memory_id = 0); + + /// update distributed synchronizer after elements' insertion + void updateDistributedSynchronizer(DistributedSynchronizer & distributed_synchronizer, + DataAccessor & data_accessor, + const Mesh & mesh_cohesive); + +protected: + + /// update elements list based on facets list + void updateElementList(Array * elements, + const Array * facets, + const Mesh & mesh_cohesive); + + /// setup facet synchronization + void setupFacetSynchronization(DistributedSynchronizer & distributed_synchronizer); + + /// build send facet arrays + void buildSendElementList(const Array *> & send_connectivity, + const Array *> & recv_connectivity, + const Array *> & temp_send_element); + + /// build recv facet arrays + void buildRecvElementList(const Array *> & temp_recv_element); + + /// get facets' global connectivity for a list of elements + template + inline void getFacetGlobalConnectivity(const DistributedSynchronizer & distributed_synchronizer, + const ElementTypeMapArray & rank_to_facet, + const Array * elements, + Array *> & connectivity, + Array *> & facets); + + /// initialize ElementTypeMap containing correspondance between + /// facets and processors + void initRankToFacet(ElementTypeMapArray & rank_to_facet); + + /// find which processor a facet is assigned to + void buildRankToFacet(ElementTypeMapArray & rank_to_facet, + const Array * elements); + + /* ------------------------------------------------------------------------ */ + /* Accessors */ + /* ------------------------------------------------------------------------ */ +public: + + /* ------------------------------------------------------------------------ */ + /* Class Members */ + /* ------------------------------------------------------------------------ */ +private: + + friend class FacetStressSynchronizer; + +}; + + +/* -------------------------------------------------------------------------- */ +/* inline functions */ +/* -------------------------------------------------------------------------- */ + +#include "facet_synchronizer_inline_impl.cc" + +__END_AKANTU__ + +#endif /* __AKANTU_FACET_SYNCHRONIZER_HH__ */ diff --git a/extra_packages/parallel-cohesive-element/src/facet_synchronizer_inline_impl.cc b/extra_packages/parallel-cohesive-element/src/facet_synchronizer_inline_impl.cc new file mode 100644 index 000000000..337cb4ed1 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/src/facet_synchronizer_inline_impl.cc @@ -0,0 +1,117 @@ +/** + * @file facet_synchronizer_inline_impl.cc + * + * @author Marco Vocialta + * + * + * @brief facet synchronizer inline implementation + * + * @section LICENSE + * + * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + */ + +/* -------------------------------------------------------------------------- */ +template +inline void FacetSynchronizer::getFacetGlobalConnectivity(const DistributedSynchronizer & distributed_synchronizer, + const ElementTypeMapArray & rank_to_facet, + const Array * elements, + Array *> & connectivity, + Array *> & facets) { + AKANTU_DEBUG_IN(); + + UInt spatial_dimension = mesh.getSpatialDimension(); + + /// init facet check tracking + ElementTypeMapArray facet_checked("facet_checked", id); + + mesh.initElementTypeMapArray(facet_checked, 1, spatial_dimension - 1); + + Mesh::type_iterator first = mesh.firstType(spatial_dimension - 1, ghost_facets); + Mesh::type_iterator last = mesh.lastType(spatial_dimension - 1, ghost_facets); + + for (; first != last; ++first) { + ElementType type = *first; + Array & f_checked = facet_checked(type, ghost_facets); + UInt nb_facet = mesh.getNbElement(type, ghost_facets); + f_checked.resize(nb_facet); + } + + const Array & nodes_global_ids = + distributed_synchronizer.mesh.getGlobalNodesIds(); + + /// loop on every processor + for (UInt p = 0; p < nb_proc; ++p) { + if (p == rank) continue; + + /// reset facet check + Mesh::type_iterator first = mesh.firstType(spatial_dimension - 1, ghost_facets); + Mesh::type_iterator last = mesh.lastType(spatial_dimension - 1, ghost_facets); + + for (; first != last; ++first) { + ElementType type = *first; + Array & f_checked = facet_checked(type, ghost_facets); + f_checked.clear(); + } + + ElementTypeMapArray & global_conn = (*connectivity(p)); + const Array & elem = elements[p]; + ElementTypeMapArray & facet_list = (*facets(p)); + + UInt nb_element = elem.getSize(); + + /// loop on every send/recv element + for (UInt el = 0; el < nb_element; ++el) { + ElementType type = elem(el).type; + GhostType gt = elem(el).ghost_type; + UInt el_index = elem(el).element; + + const Array & facet_to_element = + mesh.getSubelementToElement(type, gt); + UInt nb_facets_per_element = Mesh::getNbFacetsPerElement(type); + ElementType facet_type = Mesh::getFacetType(type); + UInt nb_nodes_per_facet = Mesh::getNbNodesPerElement(facet_type); + Vector conn_tmp(nb_nodes_per_facet); + + /// loop on every facet of the element + for (UInt f = 0; f < nb_facets_per_element; ++f) { + + const Element & facet = facet_to_element(el_index, f); + if (facet == ElementNull) continue; + UInt facet_index = facet.element; + GhostType facet_gt = facet.ghost_type; + + const Array & t_to_f = rank_to_facet(facet_type, facet_gt); + + /// exclude not ghost facets, facets assigned to other + /// processors + if (facet_gt != ghost_facets) continue; + if ((facet_gt == _ghost) && (t_to_f(facet_index) != p)) continue; + + /// exclude facets that have already been added + Array & f_checked = facet_checked(facet_type, facet_gt); + if (f_checked(facet_index)) continue; + else f_checked(facet_index) = true; + + /// add facet index + Array & f_list = facet_list(facet_type, facet_gt); + f_list.push_back(facet_index); + + /// add sorted facet global connectivity + const Array & conn = mesh.getConnectivity(facet_type, facet_gt); + Array & g_conn = global_conn(facet_type, facet_gt); + + for (UInt n = 0; n < nb_nodes_per_facet; ++n) + conn_tmp(n) = nodes_global_ids(conn(facet_index, n)); + + std::sort(conn_tmp.storage(), conn_tmp.storage() + nb_nodes_per_facet); + + g_conn.push_back(conn_tmp); + } + } + } + + AKANTU_DEBUG_OUT(); +} diff --git a/extra_packages/parallel-cohesive-element/src/solid_mechanics_model_cohesive_parallel.cc b/extra_packages/parallel-cohesive-element/src/solid_mechanics_model_cohesive_parallel.cc new file mode 100644 index 000000000..81ac1792e --- /dev/null +++ b/extra_packages/parallel-cohesive-element/src/solid_mechanics_model_cohesive_parallel.cc @@ -0,0 +1,127 @@ +/** + * @file solid_mechanics_model_cohesive_parallel.cc + * + * @author Marco Vocialta + * + * + * @brief Functions for parallel cohesive elements + * + * @section LICENSE + * + * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + */ + +/* -------------------------------------------------------------------------- */ +#include "solid_mechanics_model_cohesive.hh" + +__BEGIN_AKANTU__ + +/* -------------------------------------------------------------------------- */ +void SolidMechanicsModelCohesive::initParallel(MeshPartition * partition, + DataAccessor * data_accessor, + bool extrinsic) { + AKANTU_DEBUG_IN(); + + SolidMechanicsModel::initParallel(partition, data_accessor); + + /// create the distributed synchronizer for cohesive elements + cohesive_distributed_synchronizer = + new DistributedSynchronizer(mesh, "cohesive_distributed_synchronizer"); + + synch_registry->registerSynchronizer(*cohesive_distributed_synchronizer, + _gst_material_id); + + synch_registry->registerSynchronizer(*cohesive_distributed_synchronizer, + _gst_smm_stress); + + synch_registry->registerSynchronizer(*cohesive_distributed_synchronizer, + _gst_smm_boundary); + + synch_parallel->filterElementsByKind(cohesive_distributed_synchronizer, + _ek_cohesive); + + inserter = new CohesiveElementInserter(mesh, extrinsic, synch_parallel, + id+":cohesive_element_inserter"); + Mesh & mesh_facets = inserter->getMeshFacets(); + + facet_synchronizer = + FacetSynchronizer::createFacetSynchronizer(*synch_parallel, + mesh_facets); + + synch_registry->registerSynchronizer(*facet_synchronizer, _gst_smmc_facets); + synch_registry->registerSynchronizer(*facet_synchronizer, _gst_smmc_facets_conn); + + synchronizeGhostFacetsConnectivity(); + + /// create the facet synchronizer for extrinsic simulations + if (extrinsic) { + facet_stress_synchronizer = + FacetStressSynchronizer::createFacetStressSynchronizer(*facet_synchronizer, + mesh_facets); + + synch_registry->registerSynchronizer(*facet_stress_synchronizer, + _gst_smmc_facets_stress); + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void SolidMechanicsModelCohesive::synchronizeGhostFacetsConnectivity() { + AKANTU_DEBUG_IN(); + + StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + Int psize = comm.getNbProc(); + + if (psize > 1) { + + /// get global connectivity for not ghost facets + global_connectivity = new ElementTypeMapArray("global_connectivity", id); + + Mesh & mesh_facets = inserter->getMeshFacets(); + + mesh_facets.initElementTypeMapArray(*global_connectivity, 1, + spatial_dimension - 1, true, + _ek_regular, true); + + mesh_facets.getGlobalConnectivity(*global_connectivity, + spatial_dimension - 1, _not_ghost); + + /// communicate + synch_registry->synchronize(_gst_smmc_facets_conn); + + /// flip facets + MeshUtils::flipFacets(mesh_facets, *global_connectivity, _ghost); + + delete global_connectivity; + } + + AKANTU_DEBUG_OUT(); +} + + +/* -------------------------------------------------------------------------- */ +void SolidMechanicsModelCohesive::updateCohesiveSynchronizers() { + /// update synchronizers if needed + + if (facet_synchronizer != NULL) { + facet_synchronizer->updateDistributedSynchronizer(*cohesive_distributed_synchronizer, + *this, + mesh); + + cohesive_distributed_synchronizer->computeBufferSize(*this, _gst_material_id); + } + + if (facet_stress_synchronizer != NULL) { + const ElementTypeMapArray & prank_to_element + = synch_parallel->getPrankToElement(); + + facet_stress_synchronizer->updateFacetStressSynchronizer(*inserter, + prank_to_element, + *this); + } +} + +__END_AKANTU__ diff --git a/extra_packages/parallel-cohesive-element/src/solid_mechanics_model_cohesive_parallel.hh b/extra_packages/parallel-cohesive-element/src/solid_mechanics_model_cohesive_parallel.hh new file mode 100644 index 000000000..0b9f3e3bf --- /dev/null +++ b/extra_packages/parallel-cohesive-element/src/solid_mechanics_model_cohesive_parallel.hh @@ -0,0 +1,91 @@ +/** + * @file solid_mechanics_model_cohesive_parallel.hh + * + * @author Marco Vocialta + * @author Nicolas Richart + * + * + * @brief Include of the class members for parallelism + * + * @section LICENSE + * + * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + */ + +/* -------------------------------------------------------------------------- */ + +/* ------------------------------------------------------------------------ */ +/* Methods */ +/* ------------------------------------------------------------------------ */ +public: +/// register the tags associated with the parallel synchronizer for +/// cohesive elements +void initParallel(MeshPartition * partition, + DataAccessor * data_accessor = NULL, + bool extrinsic = false); + +protected: +void synchronizeGhostFacetsConnectivity(); + +void updateCohesiveSynchronizers(); + +/* ------------------------------------------------------------------------ */ +/* Accessors */ +/* ------------------------------------------------------------------------ */ +public: + +/// get cohesive elements synchronizer +AKANTU_GET_MACRO(CohesiveSynchronizer, + cohesive_distributed_synchronizer, + const DistributedSynchronizer *); + +/* ------------------------------------------------------------------------ */ +/* Data Accessor inherited members */ +/* ------------------------------------------------------------------------ */ +public: + +inline UInt getNbQuadsForFacetCheck(const Array & elements) const; + +inline virtual UInt getNbDataForElements(const Array & elements, + SynchronizationTag tag) const; + +inline virtual void packElementData(CommunicationBuffer & buffer, + const Array & elements, + SynchronizationTag tag) const; + +inline virtual void unpackElementData(CommunicationBuffer & buffer, + const Array & elements, + SynchronizationTag tag); + +template +inline void packFacetStressDataHelper(const ElementTypeMapArray & data_to_pack, + CommunicationBuffer & buffer, + const Array & elements) const; + +template +inline void unpackFacetStressDataHelper(ElementTypeMapArray & data_to_unpack, + CommunicationBuffer & buffer, + const Array & elements) const; + +template +inline void packUnpackFacetStressDataHelper(ElementTypeMapArray & data_to_pack, + CommunicationBuffer & buffer, + const Array & element) const; + +/* ------------------------------------------------------------------------ */ +/* Class Members */ +/* ------------------------------------------------------------------------ */ +private: +/// facet synchronizer +FacetSynchronizer * facet_synchronizer; + +/// facet stress synchronizer +FacetStressSynchronizer * facet_stress_synchronizer; + +/// cohesive elements synchronizer +DistributedSynchronizer * cohesive_distributed_synchronizer; + +/// global connectivity +ElementTypeMapArray * global_connectivity; diff --git a/extra_packages/parallel-cohesive-element/src/solid_mechanics_model_cohesive_parallel_inline_impl.cc b/extra_packages/parallel-cohesive-element/src/solid_mechanics_model_cohesive_parallel_inline_impl.cc new file mode 100644 index 000000000..cf0632825 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/src/solid_mechanics_model_cohesive_parallel_inline_impl.cc @@ -0,0 +1,313 @@ +/** + * @file solid_mechanics_model_cohesive_parallel_inline_impl.cc + * + * @author Marco Vocialta + * + * + * @brief Solid mechanics model inline implementation + * + * @section LICENSE + * + * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + */ + +/* -------------------------------------------------------------------------- */ + +__BEGIN_AKANTU__ + +/* -------------------------------------------------------------------------- */ +inline UInt SolidMechanicsModelCohesive::getNbQuadsForFacetCheck(const Array & elements) const { + UInt nb_quads = 0; + UInt nb_quad_per_facet = 0; + + ElementType current_element_type = _not_defined; + GhostType current_ghost_type = _casper; + + Array::const_iterator it = elements.begin(); + Array::const_iterator end = elements.end(); + for (; it != end; ++it) { + const Element & el = *it; + + if(el.type != current_element_type || el.ghost_type != current_ghost_type) { + current_element_type = el.type; + current_ghost_type = el.ghost_type; + + nb_quad_per_facet = this->getFEEngine("FacetsFEEngine").getNbIntegrationPoints(el.type, + el.ghost_type); + } + + nb_quads += nb_quad_per_facet; + } + + return nb_quads; +} + +/* -------------------------------------------------------------------------- */ +inline UInt SolidMechanicsModelCohesive::getNbDataForElements(const Array & elements, + SynchronizationTag tag) const { + AKANTU_DEBUG_IN(); + + UInt size = 0; + + if (elements.getSize() == 0) return size; + + /// regular element case + if (elements(0).kind == _ek_regular) { + + switch(tag) { + case _gst_smmc_facets: { + size += elements.getSize() * sizeof(bool); + break; + } + case _gst_smmc_facets_conn: { + UInt nb_nodes = Mesh::getNbNodesPerElementList(elements); + size += nb_nodes * sizeof(UInt); + break; + } + case _gst_smmc_facets_stress: { + UInt nb_quads = getNbQuadsForFacetCheck(elements); + size += nb_quads * spatial_dimension * spatial_dimension * sizeof(Real); + break; + } + default: { + size += SolidMechanicsModel::getNbDataForElements(elements, tag); + } + } + } + /// cohesive element case + else if (elements(0).kind == _ek_cohesive) { + + switch(tag) { + case _gst_material_id: { + size += elements.getSize() * sizeof(UInt); + break; + } + case _gst_smm_boundary: { + UInt nb_nodes_per_element = 0; + + Array::const_iterator it = elements.begin(); + Array::const_iterator end = elements.end(); + for (; it != end; ++it) { + const Element & el = *it; + nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); + } + + // force, displacement, boundary + size += nb_nodes_per_element * spatial_dimension * (2 * sizeof(Real) + sizeof(bool)); + break; + } + default: { } + } + + if(tag != _gst_material_id && tag != _gst_smmc_facets) { + Array * elements_per_mat = new Array[materials.size()]; + this->splitElementByMaterial(elements, elements_per_mat); + + for (UInt i = 0; i < materials.size(); ++i) { + size += materials[i]->getNbDataForElements(elements_per_mat[i], tag); + } + delete [] elements_per_mat; + } + } + + AKANTU_DEBUG_OUT(); + return size; +} + +/* -------------------------------------------------------------------------- */ +inline void SolidMechanicsModelCohesive::packElementData(CommunicationBuffer & buffer, + const Array & elements, + SynchronizationTag tag) const { + AKANTU_DEBUG_IN(); + + if (elements.getSize() == 0) return; + + if (elements(0).kind == _ek_regular) { + + switch(tag) { + + case _gst_smmc_facets: { + packElementalDataHelper(inserter->getInsertionFacetsByElement(), + buffer, elements, false, getFEEngine()); + break; + } + case _gst_smmc_facets_conn: { + packElementalDataHelper(*global_connectivity, buffer, elements, + false, getFEEngine()); + break; + } + case _gst_smmc_facets_stress: { + packFacetStressDataHelper(facet_stress, buffer, elements); + break; + } + default: { + SolidMechanicsModel::packElementData(buffer, elements, tag); + } + } + } + else if (elements(0).kind == _ek_cohesive) { + + switch(tag) { + + case _gst_material_id: { + packElementalDataHelper(material_index, buffer, + elements, false, getFEEngine("CohesiveFEEngine")); + break; + } + case _gst_smm_boundary: { + packNodalDataHelper(*force, buffer, elements, mesh); + packNodalDataHelper(*velocity, buffer, elements, mesh); + packNodalDataHelper(*blocked_dofs, buffer, elements, mesh); + break; + } + default: { } + } + + if(tag != _gst_material_id && tag != _gst_smmc_facets) { + Array * elements_per_mat = new Array[materials.size()]; + splitElementByMaterial(elements, elements_per_mat); + + for (UInt i = 0; i < materials.size(); ++i) { + materials[i]->packElementData(buffer, elements_per_mat[i], tag); + } + + delete [] elements_per_mat; + } + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +inline void SolidMechanicsModelCohesive::unpackElementData(CommunicationBuffer & buffer, + const Array & elements, + SynchronizationTag tag) { + AKANTU_DEBUG_IN(); + + if (elements.getSize() == 0) return; + + if (elements(0).kind == _ek_regular) { + + switch(tag) { + case _gst_smmc_facets: { + unpackElementalDataHelper(inserter->getInsertionFacetsByElement(), + buffer, elements, false, getFEEngine()); + break; + } + case _gst_smmc_facets_conn: { + unpackElementalDataHelper(*global_connectivity, buffer, elements, + false, getFEEngine()); + break; + } + case _gst_smmc_facets_stress: { + unpackFacetStressDataHelper(facet_stress, buffer, elements); + break; + } + default: { + SolidMechanicsModel::unpackElementData(buffer, elements, tag); + } + } + } + else if (elements(0).kind == _ek_cohesive) { + + switch(tag) { + case _gst_material_id: { + unpackElementalDataHelper(material_index, buffer, + elements, false, getFEEngine("CohesiveFEEngine")); + break; + } + case _gst_smm_boundary: { + unpackNodalDataHelper(*force, buffer, elements, mesh); + unpackNodalDataHelper(*velocity, buffer, elements, mesh); + unpackNodalDataHelper(*blocked_dofs, buffer, elements, mesh); + break; + } + default: { } + } + + if(tag != _gst_material_id && tag != _gst_smmc_facets) { + Array * elements_per_mat = new Array[materials.size()]; + splitElementByMaterial(elements, elements_per_mat); + + for (UInt i = 0; i < materials.size(); ++i) { + materials[i]->unpackElementData(buffer, elements_per_mat[i], tag); + } + + delete [] elements_per_mat; + } + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +template +inline void SolidMechanicsModelCohesive::packFacetStressDataHelper(const ElementTypeMapArray & data_to_pack, + CommunicationBuffer & buffer, + const Array & elements) const { + packUnpackFacetStressDataHelper(const_cast &>(data_to_pack), + buffer, elements); +} + +/* -------------------------------------------------------------------------- */ +template +inline void SolidMechanicsModelCohesive::unpackFacetStressDataHelper(ElementTypeMapArray & data_to_unpack, + CommunicationBuffer & buffer, + const Array & elements) const { + packUnpackFacetStressDataHelper(data_to_unpack, buffer, elements); +} + +/* -------------------------------------------------------------------------- */ +template +inline void SolidMechanicsModelCohesive::packUnpackFacetStressDataHelper(ElementTypeMapArray & data_to_pack, + CommunicationBuffer & buffer, + const Array & element) const { + ElementType current_element_type = _not_defined; + GhostType current_ghost_type = _casper; + UInt nb_quad_per_elem = 0; + UInt sp2 = spatial_dimension * spatial_dimension; + UInt nb_component = sp2 * 2; + bool element_rank = 0; + Mesh & mesh_facets = inserter->getMeshFacets(); + + Array * vect = NULL; + Array > * element_to_facet = NULL; + + Array::const_iterator it = element.begin(); + Array::const_iterator end = element.end(); + for (; it != end; ++it) { + const Element & el = *it; + if(el.type != current_element_type || el.ghost_type != current_ghost_type) { + current_element_type = el.type; + current_ghost_type = el.ghost_type; + vect = &data_to_pack(el.type, el.ghost_type); + + element_to_facet = + &( mesh_facets.getElementToSubelement(el.type, el.ghost_type) ); + + nb_quad_per_elem = this->getFEEngine("FacetsFEEngine").getNbIntegrationPoints(el.type, + el.ghost_type); + } + + if (pack_helper) + element_rank = (*element_to_facet)(el.element)[0].ghost_type != _not_ghost; + else + element_rank = (*element_to_facet)(el.element)[0].ghost_type == _not_ghost; + + for (UInt q = 0; q < nb_quad_per_elem; ++q) { + Vector data(vect->storage() + + (el.element * nb_quad_per_elem + q) * nb_component + + element_rank * sp2, + sp2); + + if(pack_helper) + buffer << data; + else + buffer >> data; + } + } +} + +__END_AKANTU__ diff --git a/extra_packages/parallel-cohesive-element/test/CMakeLists.txt b/extra_packages/parallel-cohesive-element/test/CMakeLists.txt new file mode 100644 index 000000000..6390e374d --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/CMakeLists.txt @@ -0,0 +1,27 @@ +#=============================================================================== +# @file CMakeLists.txt +# +# @author Nicolas Richart +# +# +# @brief configuration for synchronizer tests +# +# @section LICENSE +# +# Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) +# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) +# +# @section DESCRIPTION +# +#=============================================================================== + +add_akantu_test(test_parallel_cohesive "parallel_cohesive_test") + +add_mesh(test_facet_synchronizer_mesh + facet.geo 3 2) + +register_test(test_facet_synchronizer + SOURCES test_facet_synchronizer.cc test_data_accessor.hh + DEPENDS test_facet_synchronizer_mesh + PACKAGE parallel_cohesive_element + ) diff --git a/extra_packages/parallel-cohesive-element/test/facet.geo b/extra_packages/parallel-cohesive-element/test/facet.geo new file mode 100644 index 000000000..f92619ff5 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/facet.geo @@ -0,0 +1,41 @@ +h = 1.; + +Point(1) = { 1, 1, 1, h}; +Point(2) = {-1, 1, 1, h}; +Point(3) = {-1,-1, 1, h}; +Point(4) = { 1,-1, 1, h}; +Point(5) = { 1, 1,-1, h}; +Point(6) = {-1, 1,-1, h}; +Point(7) = {-1,-1,-1, h}; +Point(8) = { 1,-1,-1, h}; +Line(1) = {1, 2}; +Line(2) = {2, 3}; +Line(3) = {3, 4}; +Line(4) = {4, 1}; +Line(5) = {5, 6}; +Line(6) = {6, 7}; +Line(7) = {7, 8}; +Line(8) = {8, 5}; +Line(9) = {5, 1}; +Line(10) = {8, 4}; +Line(11) = {7, 3}; +Line(12) = {6, 2}; +Line Loop(13) = {10, 4, -9, -8}; +Plane Surface(14) = {13}; +Line Loop(15) = {7, 8, 5, 6}; +Plane Surface(16) = {15}; +Line Loop(17) = {11, -2, -12, 6}; +Plane Surface(18) = {17}; +Line Loop(19) = {1, -12, -5, 9}; +Plane Surface(20) = {19}; +Line Loop(21) = {3, -10, -7, 11}; +Plane Surface(22) = {21}; +Line Loop(23) = {3, 4, 1, 2}; +Plane Surface(24) = {23}; +Surface Loop(25) = {14, 22, 24, 20, 18, 16}; +Volume(26) = {25}; + +// Transfinite Volume "*"; +// Transfinite Line {1, 3} = 3; +// Transfinite Line {2, 4} = 3; +// Transfinite Surface "*"; diff --git a/extra_packages/parallel-cohesive-element/test/test_data_accessor.hh b/extra_packages/parallel-cohesive-element/test/test_data_accessor.hh new file mode 100644 index 000000000..8bdef93b4 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_data_accessor.hh @@ -0,0 +1,134 @@ +/** + * @file test_data_accessor.hh + * + * @author Marco Vocialta + * @author Nicolas Richart + * + * @date creation: Thu Apr 11 2013 + * @date last modification: Thu Jun 05 2014 + * + * @brief Data Accessor class for testing + * + * @section LICENSE + * + * Copyright (©) 2014 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 "data_accessor.hh" + + +/* -------------------------------------------------------------------------- */ + +__BEGIN_AKANTU__ + +class TestAccessor : public DataAccessor { + /* ------------------------------------------------------------------------ */ + /* Constructors/Destructors */ + /* ------------------------------------------------------------------------ */ +public: + inline TestAccessor(const Mesh & mesh, const ElementTypeMapArray & barycenters); + + AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Barycenter, barycenters, Real); + + /* ------------------------------------------------------------------------ */ + /* Ghost Synchronizer inherited members */ + /* ------------------------------------------------------------------------ */ +protected: + inline UInt getNbDataForElements(const Array & elements, + SynchronizationTag tag) const; + inline void packElementData(CommunicationBuffer & buffer, + const Array & elements, + SynchronizationTag tag) const; + inline void unpackElementData(CommunicationBuffer & buffer, + const Array & elements, + SynchronizationTag tag); + + /* ------------------------------------------------------------------------ */ + /* Class Members */ + /* ------------------------------------------------------------------------ */ +protected: + const ElementTypeMapArray & barycenters; + const Mesh & mesh; +}; + + +/* -------------------------------------------------------------------------- */ +/* TestSynchronizer implementation */ +/* -------------------------------------------------------------------------- */ +inline TestAccessor::TestAccessor(const Mesh & mesh, + const ElementTypeMapArray & barycenters) + : barycenters(barycenters), mesh(mesh) { } + +inline UInt TestAccessor::getNbDataForElements(const Array & elements, + __attribute__ ((unused)) SynchronizationTag tag) const { + if(elements.getSize()) + // return Mesh::getSpatialDimension(elements(0).type) * sizeof(Real) * elements.getSize(); + return mesh.getSpatialDimension() * sizeof(Real) * elements.getSize(); + else + return 0; +} + +inline void TestAccessor::packElementData(CommunicationBuffer & buffer, + const Array & elements, + __attribute__ ((unused)) SynchronizationTag tag) const { + UInt spatial_dimension = mesh.getSpatialDimension(); + Array::const_iterator bit = elements.begin(); + Array::const_iterator bend = elements.end(); + for (; bit != bend; ++bit) { + const Element & element = *bit; + + Vector bary(this->barycenters(element.type, element.ghost_type).storage() + + element.element * spatial_dimension, + spatial_dimension); + buffer << bary; + } +} + +inline void TestAccessor::unpackElementData(CommunicationBuffer & buffer, + const Array & elements, + __attribute__ ((unused)) SynchronizationTag tag) { + UInt spatial_dimension = mesh.getSpatialDimension(); + Array::const_iterator bit = elements.begin(); + Array::const_iterator bend = elements.end(); + for (; bit != bend; ++bit) { + const Element & element = *bit; + + Vector barycenter_loc(this->barycenters(element.type, element.ghost_type).storage() + + element.element * spatial_dimension, + spatial_dimension); + + Vector bary(spatial_dimension); + buffer >> bary; + std::cout << element << barycenter_loc << std::endl; + Real tolerance = 1e-15; + for (UInt i = 0; i < spatial_dimension; ++i) { + if(!(std::abs(bary(i) - barycenter_loc(i)) <= tolerance)) + AKANTU_DEBUG_ERROR("Unpacking an unknown value for the element: " + << element + << "(barycenter = " << barycenter_loc + << " and buffer = " << bary << ") direction (" << i << ")- tag: " << tag); + } + } +} + + +__END_AKANTU__ diff --git a/extra_packages/parallel-cohesive-element/test/test_facet_synchronizer.cc b/extra_packages/parallel-cohesive-element/test/test_facet_synchronizer.cc new file mode 100644 index 000000000..d3b2c8d36 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_facet_synchronizer.cc @@ -0,0 +1,100 @@ +/** + * @file test_facet_synchronizer.cc + * + * @author Marco Vocialta + * + * + * @brief Facet synchronizer test + * + * @section LICENSE + * + * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + */ + +/* -------------------------------------------------------------------------- */ + +#include "test_data_accessor.hh" +#include "facet_synchronizer.hh" +#include "mesh_utils.hh" +#include "synchronizer_registry.hh" +#include + +/* -------------------------------------------------------------------------- */ + +using namespace akantu; + +int main(int argc, char *argv[]) { + akantu::initialize(argc, argv); + + UInt spatial_dimension = 3; + Mesh mesh(spatial_dimension); + + StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + Int psize = comm.getNbProc(); + Int prank = comm.whoAmI(); + DistributedSynchronizer * dist = NULL; + + /// partition the mesh + if(prank == 0) { + mesh.read("facet.msh"); + MeshPartition * partition = new MeshPartitionScotch(mesh, spatial_dimension); + partition->partitionate(psize); + dist = DistributedSynchronizer::createDistributedSynchronizerMesh(mesh, partition); + delete partition; + } else { + dist = DistributedSynchronizer::createDistributedSynchronizerMesh(mesh, NULL); + } + + /// create facets + Mesh mesh_facets(mesh.initMeshFacets("mesh_facets")); + MeshUtils::buildAllFacets(mesh, mesh_facets, 0, dist); + + /// compute barycenter for each facet + ElementTypeMapArray barycenters("barycenters", "", 0); + mesh_facets.initElementTypeMapArray(barycenters, + spatial_dimension, + spatial_dimension - 1); + + for (ghost_type_t::iterator gt = ghost_type_t::begin(); + gt != ghost_type_t::end(); ++gt) { + GhostType ghost_type = *gt; + Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, + ghost_type); + Mesh::type_iterator last_type = mesh_facets.lastType(spatial_dimension - 1, + ghost_type); + + for(; it != last_type; ++it) { + UInt nb_element = mesh_facets.getNbElement(*it, ghost_type); + Array & barycenter = barycenters(*it, ghost_type); + barycenter.resize(nb_element); + + Array::iterator< Vector > bary_it + = barycenter.begin(spatial_dimension); + + for (UInt elem = 0; elem < nb_element; ++elem, ++bary_it) + mesh_facets.getBarycenter(elem, *it, bary_it->storage(), ghost_type); + } + } + + /// setup facet communications + FacetSynchronizer * facet_synchronizer = + FacetSynchronizer::createFacetSynchronizer(*dist, mesh_facets); + + AKANTU_DEBUG_INFO("Creating TestAccessor"); + TestAccessor test_accessor(mesh, barycenters); + SynchronizerRegistry synch_registry(test_accessor); + + synch_registry.registerSynchronizer(*facet_synchronizer, _gst_test); + + /// synchronize facets and check results + AKANTU_DEBUG_INFO("Synchronizing tag"); + synch_registry.synchronize(_gst_test); + + delete facet_synchronizer; + delete dist; + akantu::finalize(); + + return EXIT_SUCCESS; +} diff --git a/extra_packages/parallel-cohesive-element/test/test_facet_synchronizer.sh b/extra_packages/parallel-cohesive-element/test/test_facet_synchronizer.sh new file mode 100755 index 000000000..0cad18aa8 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_facet_synchronizer.sh @@ -0,0 +1,3 @@ +#!/bin/bash + +mpirun -np 3 ./test_facet_synchronizer diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/CMakeLists.txt b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/CMakeLists.txt new file mode 100644 index 000000000..0d1d2034f --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/CMakeLists.txt @@ -0,0 +1,23 @@ +#=============================================================================== +# @file CMakeLists.txt +# +# @author Marco Vocialta +# +# +# @brief configuration for cohesive elements tests +# +# @section LICENSE +# +# Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) +# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) +# +# @section DESCRIPTION +# +#=============================================================================== + +add_akantu_test(test_cohesive_parallel_intrinsic "test_cohesive_parallel_intrinsic") +add_akantu_test(test_cohesive_facet_stress_synchronizer "test_cohesive_facet_stress_synchronizer") +add_akantu_test(test_cohesive_parallel_extrinsic "test_cohesive_parallel_extrinsic") +add_akantu_test(test_cohesive_parallel_extrinsic_IG_TG "test_cohesive_parallel_extrinsic_IG_TG") +add_akantu_test(test_cohesive_parallel_buildfragments "test_cohesive_parallel_buildfragments") +add_akantu_test(test_cohesive_parallel_insertion "test_cohesive_parallel_insertion") diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_facet_stress_synchronizer/CMakeLists.txt b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_facet_stress_synchronizer/CMakeLists.txt new file mode 100644 index 000000000..d89b05e49 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_facet_stress_synchronizer/CMakeLists.txt @@ -0,0 +1,24 @@ +#=============================================================================== +# @file CMakeLists.txt +# +# @author Marco Vocialta +# +# +# @brief configuration for parallel test of facet stress synchronizer +# +# @section LICENSE +# +# Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) +# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) +# +# @section DESCRIPTION +# +#=============================================================================== + +add_mesh(test_cohesive_facet_stress_synchronizer_mesh tetrahedron.geo 3 2) + +register_test(test_cohesive_facet_stress_synchronizer + SOURCES test_cohesive_facet_stress_synchronizer.cc + DEPENDS test_cohesive_facet_stress_synchronizer_mesh + PACKAGE parallel_cohesive_element + FILES_TO_COPY material.dat) diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_facet_stress_synchronizer/material.dat b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_facet_stress_synchronizer/material.dat new file mode 100644 index 000000000..1014533aa --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_facet_stress_synchronizer/material.dat @@ -0,0 +1,21 @@ +# material plastic [ +# name = steel +# rho = 1e3 # density +# E = 1e3 # young's modulus +# nu = 0.001 # poisson's ratio +# finite_deformation = 1 +# ] + +material elastic [ + name = steel + rho = 1e3 # density + E = 1e3 # young's modulus + nu = 0.001 # poisson's ratio +] + +material cohesive_linear [ + name = cohesive + beta = 1 + G_c = 100 + sigma_c = 100 +] diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_facet_stress_synchronizer/test_cohesive_facet_stress_synchronizer.cc b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_facet_stress_synchronizer/test_cohesive_facet_stress_synchronizer.cc new file mode 100644 index 000000000..1a32681bf --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_facet_stress_synchronizer/test_cohesive_facet_stress_synchronizer.cc @@ -0,0 +1,205 @@ +/** + * @file test_cohesive_facet_stress_synchronizer.cc + * + * @author Marco Vocialta + * + * + * @brief Test for facet stress synchronizer + * + * @section LICENSE + * + * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + */ + +/* -------------------------------------------------------------------------- */ +#include +#include +#include + + +/* -------------------------------------------------------------------------- */ +#include "solid_mechanics_model_cohesive.hh" +/* -------------------------------------------------------------------------- */ + +using namespace akantu; + +Real function(Real constant, Real x, Real y, Real z) { + return constant + 2. * x + 3. * y + 4 * z; +} + +int main(int argc, char *argv[]) { + initialize("material.dat", argc, argv); + + const UInt spatial_dimension = 3; + + ElementType type = _tetrahedron_10; + ElementType type_facet = Mesh::getFacetType(type); + Math::setTolerance(1.e-11); + + Mesh mesh(spatial_dimension); + + StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + Int psize = comm.getNbProc(); + Int prank = comm.whoAmI(); + + akantu::MeshPartition * partition = NULL; + if(prank == 0) { + // Read the mesh + mesh.read("tetrahedron.msh"); + + /// partition the mesh + partition = new MeshPartitionScotch(mesh, spatial_dimension); + partition->partitionate(psize); + } + + SolidMechanicsModelCohesive model(mesh); + model.initParallel(partition, NULL, true); + model.initFull(SolidMechanicsModelCohesiveOptions(_explicit_lumped_mass, true)); + + Array & position = mesh.getNodes(); + + /* ------------------------------------------------------------------------ */ + /* Facet part */ + /* ------------------------------------------------------------------------ */ + + /// compute quadrature points positions on facets + const Mesh & mesh_facets = model.getMeshFacets(); + UInt nb_facet = mesh_facets.getNbElement(type_facet); + UInt nb_quad_per_facet = model.getFEEngine("FacetsFEEngine").getNbIntegrationPoints(type_facet); + UInt nb_tot_quad = nb_quad_per_facet * nb_facet; + + Array quad_facets(nb_tot_quad, spatial_dimension); + + model.getFEEngine("FacetsFEEngine").interpolateOnIntegrationPoints(position, + quad_facets, + spatial_dimension, + type_facet); + + /* ------------------------------------------------------------------------ */ + /* End of facet part */ + /* ------------------------------------------------------------------------ */ + + + /// compute quadrature points position of the elements + UInt nb_quad_per_element = model.getFEEngine().getNbIntegrationPoints(type); + UInt nb_element = mesh.getNbElement(type); + UInt nb_tot_quad_el = nb_quad_per_element * nb_element; + + Array quad_elements(nb_tot_quad_el, spatial_dimension); + + + model.getFEEngine().interpolateOnIntegrationPoints(position, + quad_elements, + spatial_dimension, + type); + + /// assign some values to stresses + Array & stress + = const_cast&>(model.getMaterial(0).getStress(type)); + + Array::iterator > stress_it + = stress.begin(spatial_dimension, spatial_dimension); + + for (UInt q = 0; q < nb_tot_quad_el; ++q, ++stress_it) { + + /// compute values + for (UInt i = 0; i < spatial_dimension; ++i) { + for (UInt j = i; j < spatial_dimension; ++j) { + UInt index = i * spatial_dimension + j; + (*stress_it)(i, j) = function(index, + quad_elements(q, 0), + quad_elements(q, 1), + quad_elements(q, 2)); + } + } + + /// fill symmetrical part + for (UInt i = 0; i < spatial_dimension; ++i) { + for (UInt j = 0; j < i; ++j) { + (*stress_it)(i, j) = (*stress_it)(j, i); + } + } + + // stress_it->clear(); + // for (UInt i = 0; i < spatial_dimension; ++i) + // (*stress_it)(i, i) = sigma_c * 5; + } + + + /// compute and communicate stress on facets + model.checkCohesiveStress(); + + /* ------------------------------------------------------------------------ */ + /* Check facet stress */ + /* ------------------------------------------------------------------------ */ + + const Array & facet_stress = model.getStressOnFacets(type_facet); + const Array & facet_check + = model.getElementInserter().getCheckFacets(type_facet); + const Array > & elements_to_facet + = model.getMeshFacets().getElementToSubelement(type_facet); + + Array::iterator > quad_facet_it + = quad_facets.begin(spatial_dimension); + Array::const_iterator > facet_stress_it + = facet_stress.begin(spatial_dimension, spatial_dimension * 2); + + Matrix current_stress(spatial_dimension, spatial_dimension); + + for (UInt f = 0; f < nb_facet; ++f) { + + if (!facet_check(f) || + (elements_to_facet(f)[0].ghost_type == _not_ghost && + elements_to_facet(f)[1].ghost_type == _not_ghost)) { + quad_facet_it += nb_quad_per_facet; + facet_stress_it += nb_quad_per_facet; + continue; + } + + for (UInt q = 0; q < nb_quad_per_facet; ++q, ++quad_facet_it, ++facet_stress_it) { + /// compute current_stress + for (UInt i = 0; i < spatial_dimension; ++i) { + for (UInt j = i; j < spatial_dimension; ++j) { + UInt index = i * spatial_dimension + j; + current_stress(i, j) = function(index, + (*quad_facet_it)(0), + (*quad_facet_it)(1), + (*quad_facet_it)(2)); + } + } + + /// fill symmetrical part + for (UInt i = 0; i < spatial_dimension; ++i) { + for (UInt j = 0; j < i; ++j) { + current_stress(i, j) = current_stress(j, i); + } + } + + /// compare it to interpolated one + for (UInt s = 0; s < 2; ++s) { + Matrix stress_to_check(facet_stress_it->storage() + + s * spatial_dimension * spatial_dimension, + spatial_dimension, + spatial_dimension); + + + for (UInt i = 0; i < spatial_dimension; ++i) { + for (UInt j = 0; j < spatial_dimension; ++j) { + if (!Math::are_float_equal(stress_to_check(i, j), current_stress(i, j))) { + std::cout << "Stress doesn't match" << std::endl; + finalize(); + return EXIT_FAILURE; + } + } + } + } + } + } + + finalize(); + if (prank == 0) + std::cout << "OK: test_cohesive_facet_stress_synchronizer passed!" << std::endl; + return EXIT_SUCCESS; +} diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_facet_stress_synchronizer/test_cohesive_facet_stress_synchronizer.sh b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_facet_stress_synchronizer/test_cohesive_facet_stress_synchronizer.sh new file mode 100755 index 000000000..2d231d523 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_facet_stress_synchronizer/test_cohesive_facet_stress_synchronizer.sh @@ -0,0 +1,3 @@ +#! /bin/bash + +mpirun -np 8 ./test_cohesive_facet_stress_synchronizer diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_facet_stress_synchronizer/tetrahedron.geo b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_facet_stress_synchronizer/tetrahedron.geo new file mode 100644 index 000000000..73ecd4687 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_facet_stress_synchronizer/tetrahedron.geo @@ -0,0 +1,63 @@ +lc = 1; +h = 1; +Point(1) = {lc, lc, lc, h}; +Point(2) = {-lc, lc, lc, h}; +Point(3) = {-lc, -lc, lc, h}; +Point(4) = {lc, -lc, lc, h}; +Point(5) = {lc, lc, -lc, h}; +Point(6) = {-lc, lc, -lc, h}; +Point(7) = {-lc, -lc, -lc, h}; +Point(8) = {lc, -lc, -lc, h}; +Point(9) = {0, lc, lc, h}; +Point(10) = {0, -lc, lc, h}; +Point(11) = {0, -lc, -lc, h}; +Point(12) = {0, lc, -lc, h}; +Line(1) = {2, 3}; +Line(2) = {3, 10}; +Line(3) = {10, 9}; +Line(4) = {9, 2}; +Line(5) = {10, 4}; +Line(6) = {4, 1}; +Line(7) = {1, 9}; +Line(8) = {10, 11}; +Line(9) = {4, 8}; +Line(10) = {3, 7}; +Line(11) = {7, 11}; +Line(12) = {11, 8}; +Line(13) = {8, 5}; +Line(14) = {5, 1}; +Line(15) = {12, 9}; +Line(16) = {6, 2}; +Line(17) = {7, 6}; +Line(18) = {6, 12}; +Line(19) = {12, 5}; +Line(20) = {11, 12}; +Line Loop(21) = {8, 20, 15, -3}; +Plane Surface(22) = {21}; +Line Loop(23) = {9, 13, 14, -6}; +Plane Surface(24) = {23}; +Line Loop(25) = {10, 17, 16, 1}; +Plane Surface(26) = {25}; +Line Loop(27) = {8, -11, -10, 2}; +Plane Surface(28) = {27}; +Line Loop(29) = {9, -12, -8, 5}; +Plane Surface(30) = {29}; +Line Loop(31) = {7, -15, 19, 14}; +Plane Surface(32) = {31}; +Line Loop(33) = {4, -16, 18, 15}; +Plane Surface(34) = {33}; +Line Loop(35) = {5, 6, 7, -3}; +Plane Surface(36) = {35}; +Line Loop(37) = {2, 3, 4, 1}; +Plane Surface(38) = {37}; +Line Loop(39) = {18, -20, -11, 17}; +Plane Surface(40) = {39}; +Line Loop(41) = {19, -13, -12, 20}; +Plane Surface(42) = {41}; +Surface Loop(43) = {34, 38, 28, 40, 26, 22}; +Volume(44) = {43}; +Surface Loop(45) = {42, 32, 36, 30, 24, 22}; +Volume(46) = {45}; + +Physical Volume(47) = {44}; +Physical Volume(48) = {46}; diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_buildfragments/CMakeLists.txt b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_buildfragments/CMakeLists.txt new file mode 100644 index 000000000..d475f758b --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_buildfragments/CMakeLists.txt @@ -0,0 +1,25 @@ +#=============================================================================== +# @file CMakeLists.txt +# +# @author Marco Vocialta +# +# +# @brief configuration for build fragments tests +# +# @section LICENSE +# +# Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) +# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) +# +# @section DESCRIPTION +# +#=============================================================================== + +add_mesh(test_cohesive_parallel_buildfragments_mesh mesh.geo 3 2) + +register_test(test_cohesive_parallel_buildfragments + SOURCES test_cohesive_parallel_buildfragments.cc + DEPENDS test_cohesive_parallel_buildfragments_mesh + FILES_TO_COPY material.dat + PACKAGE parallel_cohesive_element + DIRECTORIES_TO_CREATE paraview) diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_buildfragments/material.dat b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_buildfragments/material.dat new file mode 100644 index 000000000..10adbb716 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_buildfragments/material.dat @@ -0,0 +1,15 @@ +seed = 1 + +material elastic [ + name = bulk + rho = 2500 + E = 42.e9 + nu = 0.2 +] + +material cohesive_linear [ + name = cohesive + beta = 1 + G_c = 1e-6 + sigma_c = 1 +] diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_buildfragments/mesh.geo b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_buildfragments/mesh.geo new file mode 100644 index 000000000..11c5ecfc1 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_buildfragments/mesh.geo @@ -0,0 +1,14 @@ +L = 0.03; +h = L/5; +s = L/20; + +Point(1) = {-L/2., 0, 0, s}; + +line[] = Extrude{L,0,0}{ Point{1}; }; +surface[] = Extrude{0,h,0}{ Line{line[1]}; }; +volume[] = Extrude{0,0,h}{ Surface{surface[1]}; }; + +Physical Volume(1) = {volume[1]}; + +Transfinite Surface "*"; +Transfinite Volume "*"; diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_buildfragments/test_cohesive_parallel_buildfragments.cc b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_buildfragments/test_cohesive_parallel_buildfragments.cc new file mode 100644 index 000000000..0484a2351 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_buildfragments/test_cohesive_parallel_buildfragments.cc @@ -0,0 +1,432 @@ +/** + * @file test_cohesive_parallel_buildfragments.cc + * + * @author Marco Vocialta + * + * + * @brief Test to build fragments in parallel + * + * @section LICENSE + * + * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + */ + +/* -------------------------------------------------------------------------- */ +#include +#include +#include +#include +#include +/* -------------------------------------------------------------------------- */ +#include "solid_mechanics_model_cohesive.hh" +#include "material_cohesive.hh" +#include "fragment_manager.hh" +/* -------------------------------------------------------------------------- */ + +using namespace akantu; + +void verticalInsertionLimit(SolidMechanicsModelCohesive &); +void displaceElements(SolidMechanicsModelCohesive &, const Real, const Real); +bool isInertiaEqual(const Vector &, const Vector &); +void rotateArray(Array & array, Real angle); +UInt getNbBigFragments(FragmentManager &, UInt); + +const UInt spatial_dimension = 3; +const UInt total_nb_fragment = 4; +const Real rotation_angle = M_PI / 4.; +const Real global_tolerance = 1.e-9; + +int main(int argc, char *argv[]) { + initialize("material.dat", argc, argv); + + Math::setTolerance(global_tolerance); + + Mesh mesh(spatial_dimension); + + StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + Int psize = comm.getNbProc(); + Int prank = comm.whoAmI(); + + akantu::MeshPartition * partition = NULL; + if(prank == 0) { + // Read the mesh + mesh.read("mesh.msh"); + + /// partition the mesh + MeshUtils::purifyMesh(mesh); + partition = new MeshPartitionScotch(mesh, spatial_dimension); + partition->partitionate(psize); + } + + SolidMechanicsModelCohesive model(mesh); + model.initParallel(partition, NULL, true); + + delete partition; + + /// model initialization + model.initFull(SolidMechanicsModelCohesiveOptions(_explicit_lumped_mass, true)); + + mesh.computeBoundingBox(); + Real L = mesh.getUpperBounds()(0) - mesh.getLowerBounds()(0); + Real h = mesh.getUpperBounds()(1) - mesh.getLowerBounds()(1); + Real rho = model.getMaterial("bulk").getParam("rho"); + + Real theoretical_mass = L * h * h * rho; + Real frag_theo_mass = theoretical_mass / total_nb_fragment; + + UInt nb_element = mesh.getNbElement(spatial_dimension, _not_ghost, _ek_regular); + comm.allReduce(&nb_element, 1, _so_sum); + UInt nb_element_per_fragment = nb_element / total_nb_fragment; + + FragmentManager fragment_manager(model); + + fragment_manager.computeAllData(); + getNbBigFragments(fragment_manager, nb_element_per_fragment + 1); + + model.setBaseName("extrinsic"); + model.addDumpFieldVector("displacement"); + model.addDumpField("velocity"); + model.addDumpField("stress"); + model.addDumpField("partitions"); + model.addDumpField("fragments"); + model.addDumpField("fragments mass"); + model.addDumpField("moments of inertia"); + model.addDumpField("elements per fragment"); + model.dump(); + + model.setBaseNameToDumper("cohesive elements", "cohesive_elements_test"); + model.addDumpFieldVectorToDumper("cohesive elements", "displacement"); + model.addDumpFieldToDumper("cohesive elements", "damage"); + model.dump("cohesive elements"); + + /// set check facets + verticalInsertionLimit(model); + + model.assembleMassLumped(); + model.synchronizeBoundaries(); + + /// impose initial displacement + Array & displacement = model.getDisplacement(); + Array & velocity = model.getVelocity(); + const Array & position = mesh.getNodes(); + UInt nb_nodes = mesh.getNbNodes(); + + for (UInt n = 0; n < nb_nodes; ++n) { + displacement(n, 0) = position(n, 0) * 0.1; + velocity(n, 0) = position(n, 0); + } + + rotateArray(mesh.getNodes(), rotation_angle); + // rotateArray(displacement, rotation_angle); + // rotateArray(velocity, rotation_angle); + + model.updateResidual(); + model.checkCohesiveStress(); + model.dump(); + model.dump("cohesive elements"); + + const Array & fragment_mass = fragment_manager.getMass(); + const Array & fragment_center = fragment_manager.getCenterOfMass(); + + Real el_size = L / total_nb_fragment; + Real lim = -L/2 + el_size * 0.99; + + /// define theoretical inertia moments + Vector small_frag_inertia(spatial_dimension); + small_frag_inertia(0) = frag_theo_mass * (h*h + h*h) / 12.; + small_frag_inertia(1) = frag_theo_mass * (el_size*el_size + h*h) / 12.; + small_frag_inertia(2) = frag_theo_mass * (el_size*el_size + h*h) / 12.; + + std::sort(small_frag_inertia.storage(), + small_frag_inertia.storage() + spatial_dimension, + std::greater()); + + const Array & inertia_moments = fragment_manager.getMomentsOfInertia(); + Array::const_iterator< Vector > inertia_moments_begin + = inertia_moments.begin(spatial_dimension); + + /// displace one fragment each time + for (UInt frag = 1; frag <= total_nb_fragment; ++frag) { + if (prank == 0) + std::cout << "Generating fragment: " << frag << std::endl; + + fragment_manager.computeAllData(); + + /// check number of big fragments + UInt nb_big_fragment = getNbBigFragments(fragment_manager, + nb_element_per_fragment + 1); + + model.dump(); + model.dump("cohesive elements"); + + + if (frag < total_nb_fragment) { + if (nb_big_fragment != 1) { + AKANTU_DEBUG_ERROR("The number of big fragments is wrong: " << nb_big_fragment); + return EXIT_FAILURE; + } + } + else { + if (nb_big_fragment != 0) { + AKANTU_DEBUG_ERROR("The number of big fragments is wrong: " << nb_big_fragment); + return EXIT_FAILURE; + } + } + + /// check number of fragments + UInt nb_fragment_num = fragment_manager.getNbFragment(); + + if (nb_fragment_num != frag) { + AKANTU_DEBUG_ERROR("The number of fragments is wrong! Numerical: " << nb_fragment_num << " Theoretical: " << frag); + return EXIT_FAILURE; + } + + /// check mass computation + if (frag < total_nb_fragment) { + Real total_mass = 0.; + UInt small_fragments = 0; + + for (UInt f = 0; f < nb_fragment_num; ++f) { + const Vector & current_inertia = inertia_moments_begin[f]; + + if (Math::are_float_equal(fragment_mass(f, 0), frag_theo_mass)) { + + /// check center of mass + if (fragment_center(f, 0) > (L * frag / total_nb_fragment - L / 2)) { + AKANTU_DEBUG_ERROR("Fragment center is wrong!"); + return EXIT_FAILURE; + } + + /// check moment of inertia + if (!isInertiaEqual(current_inertia, small_frag_inertia)) { + AKANTU_DEBUG_ERROR("Inertia moments are wrong"); + return EXIT_FAILURE; + } + + ++small_fragments; + total_mass += frag_theo_mass; + } + else { + /// check the moment of inertia for the biggest fragment + Real big_frag_mass = frag_theo_mass * (total_nb_fragment - frag + 1); + Real big_frag_size = el_size * (total_nb_fragment - frag + 1); + + Vector big_frag_inertia(spatial_dimension); + big_frag_inertia(0) = big_frag_mass * (h*h + h*h) / 12.; + big_frag_inertia(1) + = big_frag_mass * (big_frag_size*big_frag_size + h*h) / 12.; + big_frag_inertia(2) + = big_frag_mass * (big_frag_size*big_frag_size + h*h) / 12.; + + std::sort(big_frag_inertia.storage(), + big_frag_inertia.storage() + spatial_dimension, + std::greater()); + + if (!isInertiaEqual(current_inertia, big_frag_inertia)) { + AKANTU_DEBUG_ERROR("Inertia moments are wrong"); + return EXIT_FAILURE; + } + + } + } + + if (small_fragments != nb_fragment_num - 1) { + AKANTU_DEBUG_ERROR("The number of small fragments is wrong!"); + return EXIT_FAILURE; + } + + if (!Math::are_float_equal(total_mass, + small_fragments * frag_theo_mass)) { + AKANTU_DEBUG_ERROR("The mass of small fragments is wrong!"); + return EXIT_FAILURE; + } + } + + + /// displace fragments + rotateArray(mesh.getNodes(), -rotation_angle); + // rotateArray(displacement, -rotation_angle); + displaceElements(model, lim, el_size * 2); + rotateArray(mesh.getNodes(), rotation_angle); + // rotateArray(displacement, rotation_angle); + + model.updateResidual(); + + lim += el_size; + } + + model.dump(); + model.dump("cohesive elements"); + + /// check centers + const Array & fragment_velocity = fragment_manager.getVelocity(); + + Real initial_position = -L / 2. + el_size / 2.; + + for (UInt frag = 0; frag < total_nb_fragment; ++frag) { + Real theoretical_center = initial_position + el_size * frag; + + UInt f_index = 0; + while (f_index < total_nb_fragment && + !Math::are_float_equal(fragment_center(f_index, 0), theoretical_center)) + ++f_index; + + if (f_index == total_nb_fragment) { + AKANTU_DEBUG_ERROR("The fragments' center is wrong!"); + return EXIT_FAILURE; + } + + f_index = 0; + while (f_index < total_nb_fragment && + !Math::are_float_equal(fragment_velocity(f_index, 0), + theoretical_center)) + ++f_index; + + if (f_index == total_nb_fragment) { + AKANTU_DEBUG_ERROR("The fragments' velocity is wrong!"); + return EXIT_FAILURE; + } + } + + finalize(); + + if (prank == 0) + std::cout << "OK: test_cohesive_buildfragments was passed!" << std::endl; + return EXIT_SUCCESS; +} + +void verticalInsertionLimit(SolidMechanicsModelCohesive & model) { + UInt spatial_dimension = model.getSpatialDimension(); + const Mesh & mesh_facets = model.getMeshFacets(); + const Array & position = mesh_facets.getNodes(); + + for (ghost_type_t::iterator gt = ghost_type_t::begin(); + gt != ghost_type_t::end(); + ++gt) { + GhostType ghost_type = *gt; + + Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, ghost_type); + Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1, ghost_type); + for(; it != end; ++it) { + ElementType type = *it; + + Array & check_facets + = model.getElementInserter().getCheckFacets(type, ghost_type); + + const Array & connectivity = mesh_facets.getConnectivity(type, ghost_type); + UInt nb_nodes_per_facet = connectivity.getNbComponent(); + + for (UInt f = 0; f < check_facets.getSize(); ++f) { + if (!check_facets(f)) continue; + + UInt nb_aligned_nodes = 1; + Real first_node_pos = position(connectivity(f, 0), 0); + + for (; nb_aligned_nodes < nb_nodes_per_facet; ++nb_aligned_nodes) { + Real other_node_pos = position(connectivity(f, nb_aligned_nodes), 0); + + if (! Math::are_float_equal(first_node_pos, other_node_pos)) + break; + } + + if (nb_aligned_nodes != nb_nodes_per_facet) { + check_facets(f) = false; + } + } + } + } +} + +void displaceElements(SolidMechanicsModelCohesive & model, + const Real lim, + const Real amount) { + UInt spatial_dimension = model.getSpatialDimension(); + Array & displacement = model.getDisplacement(); + Mesh & mesh = model.getMesh(); + UInt nb_nodes = mesh.getNbNodes(); + Array displaced(nb_nodes); + displaced.clear(); + Vector barycenter(spatial_dimension); + + for (ghost_type_t::iterator gt = ghost_type_t::begin(); + gt != ghost_type_t::end(); + ++gt) { + GhostType ghost_type = *gt; + + Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); + Mesh::type_iterator end = mesh.lastType(spatial_dimension, ghost_type); + for(; it != end; ++it) { + ElementType type = *it; + const Array & connectivity = mesh.getConnectivity(type, ghost_type); + UInt nb_element = connectivity.getSize(); + UInt nb_nodes_per_element = connectivity.getNbComponent(); + + Array::const_vector_iterator conn_el + = connectivity.begin(nb_nodes_per_element); + + for (UInt el = 0; el < nb_element; ++el) { + mesh.getBarycenter(el, type, barycenter.storage(), ghost_type); + + if (barycenter(0) < lim) { + const Vector & conn = conn_el[el]; + + for (UInt n = 0; n < nb_nodes_per_element; ++n) { + UInt node = conn(n); + if (!displaced(node)) { + displacement(node, 0) -= amount; + displaced(node) = true; + } + } + } + } + } + } +} + +bool isInertiaEqual(const Vector & a, const Vector & b) { + UInt nb_terms = a.size(); + UInt equal_terms = 0; + + while (equal_terms < nb_terms && + std::abs(a(equal_terms) - b(equal_terms)) / a(equal_terms) < Math::getTolerance()) + ++equal_terms; + + return equal_terms == nb_terms; +} + +void rotateArray(Array & array, Real angle) { + UInt spatial_dimension = array.getNbComponent(); + + Real rotation_values[] = {std::cos(angle), std::sin(angle), 0, + -std::sin(angle), std::cos(angle), 0, + 0, 0, 1}; + Matrix rotation(rotation_values, spatial_dimension, spatial_dimension); + + RVector displaced_node(spatial_dimension); + Array::vector_iterator node_it = array.begin(spatial_dimension); + Array::vector_iterator node_end = array.end(spatial_dimension); + + for (; node_it != node_end; ++node_it) { + displaced_node.mul(rotation, *node_it); + *node_it = displaced_node; + } +} + +UInt getNbBigFragments(FragmentManager & fragment_manager, + UInt minimum_nb_elements) { + fragment_manager.computeNbElementsPerFragment(); + const Array & nb_elements_per_fragment + = fragment_manager.getNbElementsPerFragment(); + UInt nb_fragment = fragment_manager.getNbFragment(); + UInt nb_big_fragment = 0; + + for (UInt frag = 0; frag < nb_fragment; ++frag) { + if (nb_elements_per_fragment(frag) >= minimum_nb_elements) { + ++nb_big_fragment; + } + } + + return nb_big_fragment; +} diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_buildfragments/test_cohesive_parallel_buildfragments.sh b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_buildfragments/test_cohesive_parallel_buildfragments.sh new file mode 100755 index 000000000..a07b277fe --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_buildfragments/test_cohesive_parallel_buildfragments.sh @@ -0,0 +1,6 @@ +#!/bin/bash + +rm -r paraview +mkdir paraview + +mpirun -np 8 ./test_cohesive_parallel_buildfragments diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/CMakeLists.txt b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/CMakeLists.txt new file mode 100644 index 000000000..59d9da668 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/CMakeLists.txt @@ -0,0 +1,40 @@ +#=============================================================================== +# @file CMakeLists.txt +# +# @author Marco Vocialta +# +# +# @brief configuration for parallel test for extrinsic cohesive elements +# +# @section LICENSE +# +# Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) +# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) +# +# @section DESCRIPTION +# +#=============================================================================== + +add_mesh(test_cohesive_parallel_extrinsic_mesh mesh.geo 2 2) +add_mesh(test_cohesive_parallel_extrinsic_tetrahedron_mesh tetrahedron.geo 3 2) + +register_test(test_cohesive_parallel_extrinsic + SOURCES test_cohesive_parallel_extrinsic.cc + DEPENDS test_cohesive_parallel_extrinsic_mesh + PACKAGE parallel_cohesive_element + FILES_TO_COPY material.dat + DIRECTORIES_TO_CREATE paraview) + +register_test(test_cohesive_parallel_extrinsic_tetrahedron + SOURCES test_cohesive_parallel_extrinsic_tetrahedron.cc + DEPENDS test_cohesive_parallel_extrinsic_tetrahedron_mesh + PACKAGE parallel_cohesive_element + FILES_TO_COPY material.dat + DIRECTORIES_TO_CREATE paraview) + +register_test(test_cohesive_parallel_extrinsic_tetrahedron_displacement + SOURCES test_cohesive_parallel_extrinsic_tetrahedron_displacement.cc + DEPENDS test_cohesive_parallel_extrinsic_tetrahedron_mesh + PACKAGE parallel_cohesive_element + FILES_TO_COPY material.dat + DIRECTORIES_TO_CREATE paraview displacement) diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/material.dat b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/material.dat new file mode 100644 index 000000000..1014533aa --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/material.dat @@ -0,0 +1,21 @@ +# material plastic [ +# name = steel +# rho = 1e3 # density +# E = 1e3 # young's modulus +# nu = 0.001 # poisson's ratio +# finite_deformation = 1 +# ] + +material elastic [ + name = steel + rho = 1e3 # density + E = 1e3 # young's modulus + nu = 0.001 # poisson's ratio +] + +material cohesive_linear [ + name = cohesive + beta = 1 + G_c = 100 + sigma_c = 100 +] diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/mesh.geo b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/mesh.geo new file mode 100644 index 000000000..70b751674 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/mesh.geo @@ -0,0 +1,17 @@ +h = 1.; + +Point(1) = { 1, 1, 0, h}; +Point(2) = {-1, 1, 0, h}; +Point(3) = {-1,-1, 0, h}; +Point(4) = { 1,-1, 0, h}; +Line(1) = {1, 2}; +Line(2) = {2, 3}; +Line(3) = {3, 4}; +Line(4) = {4, 1}; +Line Loop(5) = {1, 2, 3, 4}; +Plane Surface(6) = {5}; +Physical Surface(7) = {6}; + +Transfinite Line {1, 3} = 3; +Transfinite Line {2, 4} = 3; +// Transfinite Surface "*"; diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/test_cohesive_parallel_extrinsic.cc b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/test_cohesive_parallel_extrinsic.cc new file mode 100644 index 000000000..f5b83daac --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/test_cohesive_parallel_extrinsic.cc @@ -0,0 +1,172 @@ +/** + * @file test_cohesive_parallel_extrinsic.cc + * + * @author Marco Vocialta + * + * + * @brief parallel test for cohesive elements + * + * @section LICENSE + * + * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + */ + +/* -------------------------------------------------------------------------- */ +#include "solid_mechanics_model_cohesive.hh" +#include "dumper_paraview.hh" +/* -------------------------------------------------------------------------- */ +using namespace akantu; + +int main(int argc, char *argv[]) { + initialize("material.dat", argc, argv); + + const UInt max_steps = 500; + + UInt spatial_dimension = 2; + Mesh mesh(spatial_dimension); + + StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + Int psize = comm.getNbProc(); + Int prank = comm.whoAmI(); + + akantu::MeshPartition * partition = NULL; + if(prank == 0) { + // Read the mesh + mesh.read("mesh.msh"); + + /// partition the mesh + partition = new MeshPartitionScotch(mesh, spatial_dimension); + // debug::setDebugLevel(dblDump); + partition->partitionate(psize); + // debug::setDebugLevel(dblWarning); + } + + SolidMechanicsModelCohesive model(mesh); + + model.initParallel(partition, NULL, true); + + // debug::setDebugLevel(dblDump); + // std::cout << mesh << std::endl; + // debug::setDebugLevel(dblWarning); + + + model.initFull(SolidMechanicsModelCohesiveOptions(_explicit_lumped_mass, true)); + model.limitInsertion(_y, -0.30, -0.20); + model.updateAutomaticInsertion(); + + // debug::setDebugLevel(dblDump); + // std::cout << mesh_facets << std::endl; + // debug::setDebugLevel(dblWarning); + + Real time_step = model.getStableTimeStep()*0.1; + model.setTimeStep(time_step); + std::cout << "Time step: " << time_step << std::endl; + + model.assembleMassLumped(); + + + Array & position = mesh.getNodes(); + Array & velocity = model.getVelocity(); + Array & boundary = model.getBlockedDOFs(); + Array & displacement = model.getDisplacement(); + // const Array & residual = model.getResidual(); + + UInt nb_nodes = mesh.getNbNodes(); + + /// boundary conditions + for (UInt n = 0; n < nb_nodes; ++n) { + if (position(n, 1) > 0.99 || position(n, 1) < -0.99) + boundary(n, 1) = true; + + if (position(n, 0) > 0.99 || position(n, 0) < -0.99) + boundary(n, 0) = true; + } + + /// initial conditions + Real loading_rate = 0.5; + Real disp_update = loading_rate * time_step; + for (UInt n = 0; n < nb_nodes; ++n) { + velocity(n, 1) = loading_rate * position(n, 1); + } + + model.synchronizeBoundaries(); + model.updateResidual(); + + model.setBaseName("extrinsic_parallel"); + model.addDumpFieldVector("displacement"); + model.addDumpField("velocity" ); + model.addDumpField("acceleration"); + model.addDumpField("residual" ); + model.addDumpField("stress"); + model.addDumpField("grad_u"); + model.addDumpField("partitions"); + // model.getDumper().getDumper().setMode(iohelper::BASE64); + model.dump(); + + model.setBaseNameToDumper("cohesive elements", + "extrinsic_parallel_cohesive_elements"); + model.addDumpFieldVectorToDumper("cohesive elements", "displacement"); + model.addDumpFieldToDumper("cohesive elements", "damage"); + model.dump("cohesive elements"); + + /// Main loop + for (UInt s = 1; s <= max_steps; ++s) { + + /// update displacement on extreme nodes + for (UInt n = 0; n < nb_nodes; ++n) { + if (position(n, 1) > 0.99 || position(n, 1) < -0.99) + displacement(n, 1) += disp_update * position(n, 1); + } + + model.checkCohesiveStress(); + model.solveStep(); + + // model.dump(); + if(s % 10 == 0) { + if(prank == 0) std::cout << "passing step " << s << "/" << max_steps << std::endl; + } + + // // update displacement + // for (UInt n = 0; n < nb_nodes; ++n) { + // if (position(n, 1) + displacement(n, 1) > 0) { + // displacement(n, 0) -= 0.01; + // } + // } + + // Real Ed = dynamic_cast (model.getMaterial(1)).getDissipatedEnergy(); + // Real Er = dynamic_cast (model.getMaterial(1)).getReversibleEnergy(); + + // edis << s << " " + // << Ed << std::endl; + + // erev << s << " " + // << Er << std::endl; + + } + + model.dump(); + model.dump("cohesive elements"); + + // edis.close(); + // erev.close(); + + Real Ed = model.getEnergy("dissipated"); + + Real Edt = 200 * sqrt(2); + + + if(prank == 0) { + std::cout << Ed << " " << Edt << std::endl; + + if (Ed < Edt * 0.999 || Ed > Edt * 1.001 || std::isnan(Ed)) { + std::cout << "The dissipated energy is incorrect" << std::endl; + finalize(); + return EXIT_FAILURE; + } + } + + finalize(); + return EXIT_SUCCESS; +} diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/test_cohesive_parallel_extrinsic.sh b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/test_cohesive_parallel_extrinsic.sh new file mode 100755 index 000000000..6b5f5e70d --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/test_cohesive_parallel_extrinsic.sh @@ -0,0 +1,6 @@ +#!/bin/bash + +rm -r paraview +mkdir paraview + +mpirun -np 8 ./test_cohesive_parallel_extrinsic diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/test_cohesive_parallel_extrinsic_tetrahedron.cc b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/test_cohesive_parallel_extrinsic_tetrahedron.cc new file mode 100644 index 000000000..05e10df9c --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/test_cohesive_parallel_extrinsic_tetrahedron.cc @@ -0,0 +1,230 @@ +/** + * @file test_cohesive_parallel_extrinsic_tetrahedron.cc + * + * @author Marco Vocialta + * + * + * @brief 3D extrinsic cohesive elements test + * + * @section LICENSE + * + * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + */ + +/* -------------------------------------------------------------------------- */ +#include "solid_mechanics_model_cohesive.hh" +#include "material_cohesive_linear.hh" + +/* -------------------------------------------------------------------------- */ +using namespace akantu; + +Real function(Real constant, Real x, Real y, Real z) { + return constant + 2. * x + 3. * y + 4 * z; +} + +int main(int argc, char *argv[]) { + initialize("material.dat", argc, argv); + + debug::setDebugLevel(dblWarning); + + // const UInt max_steps = 1000; + // Real increment = 0.005; + const UInt spatial_dimension = 3; + Math::setTolerance(1.e-12); + + ElementType type = _tetrahedron_10; + ElementType type_facet = Mesh::getFacetType(type); + ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); + + Mesh mesh(spatial_dimension); + + StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + Int psize = comm.getNbProc(); + Int prank = comm.whoAmI(); + + akantu::MeshPartition * partition = NULL; + if(prank == 0) { + // Read the mesh + mesh.read("tetrahedron.msh"); + + /// partition the mesh + partition = new MeshPartitionScotch(mesh, spatial_dimension); + partition->partitionate(psize); + } + + SolidMechanicsModelCohesive model(mesh); + + /// model initialization + model.initParallel(partition, NULL, true); + model.initFull(SolidMechanicsModelCohesiveOptions(_explicit_lumped_mass, true)); + + const MaterialCohesiveLinear<3> & mat_cohesive + = dynamic_cast < const MaterialCohesiveLinear<3> & > (model.getMaterial(1)); + + const Real sigma_c = mat_cohesive.getParam< RandomInternalField >("sigma_c"); + const Real beta = mat_cohesive.getParam("beta"); + // const Real G_cI = mat_cohesive.getParam("G_cI"); + + Array & position = mesh.getNodes(); + + /* ------------------------------------------------------------------------ */ + /* Facet part */ + /* ------------------------------------------------------------------------ */ + + /// compute quadrature points positions on facets + const Mesh & mesh_facets = model.getMeshFacets(); + UInt nb_facet = mesh_facets.getNbElement(type_facet); + UInt nb_quad_per_facet = model.getFEEngine("FacetsFEEngine").getNbIntegrationPoints(type_facet); + UInt nb_tot_quad = nb_quad_per_facet * nb_facet; + + Array quad_facets(nb_tot_quad, spatial_dimension); + + model.getFEEngine("FacetsFEEngine").interpolateOnIntegrationPoints(position, + quad_facets, + spatial_dimension, + type_facet); + + /* ------------------------------------------------------------------------ */ + /* End of facet part */ + /* ------------------------------------------------------------------------ */ + + + /// compute quadrature points position of the elements + UInt nb_quad_per_element = model.getFEEngine().getNbIntegrationPoints(type); + UInt nb_element = mesh.getNbElement(type); + UInt nb_tot_quad_el = nb_quad_per_element * nb_element; + + Array quad_elements(nb_tot_quad_el, spatial_dimension); + + + model.getFEEngine().interpolateOnIntegrationPoints(position, + quad_elements, + spatial_dimension, + type); + + /// assign some values to stresses + Array & stress + = const_cast&>(model.getMaterial(0).getStress(type)); + + Array::iterator > stress_it + = stress.begin(spatial_dimension, spatial_dimension); + + for (UInt q = 0; q < nb_tot_quad_el; ++q, ++stress_it) { + + /// compute values + for (UInt i = 0; i < spatial_dimension; ++i) { + for (UInt j = i; j < spatial_dimension; ++j) { + UInt index = i * spatial_dimension + j; + (*stress_it)(i, j) = index * function(sigma_c * 5, + quad_elements(q, 0), + quad_elements(q, 1), + quad_elements(q, 2)); + } + } + + /// fill symmetrical part + for (UInt i = 0; i < spatial_dimension; ++i) { + for (UInt j = 0; j < i; ++j) { + (*stress_it)(i, j) = (*stress_it)(j, i); + } + } + } + + + /// compute stress on facet quads + Array stress_facets(nb_tot_quad, spatial_dimension * spatial_dimension); + + Array::iterator > stress_facets_it + = stress_facets.begin(spatial_dimension, spatial_dimension); + + for (UInt q = 0; q < nb_tot_quad; ++q, ++stress_facets_it) { + /// compute values + for (UInt i = 0; i < spatial_dimension; ++i) { + for (UInt j = i; j < spatial_dimension; ++j) { + UInt index = i * spatial_dimension + j; + (*stress_facets_it)(i, j) = index * function(sigma_c * 5, + quad_facets(q, 0), + quad_facets(q, 1), + quad_facets(q, 2)); + } + } + + /// fill symmetrical part + for (UInt i = 0; i < spatial_dimension; ++i) { + for (UInt j = 0; j < i; ++j) { + (*stress_facets_it)(i, j) = (*stress_facets_it)(j, i); + } + } + } + + + /// insert cohesive elements + model.checkCohesiveStress(); + + + /// check insertion stress + const Array & normals = + model.getFEEngine("FacetsFEEngine").getNormalsOnIntegrationPoints(type_facet); + const Array & tangents = model.getTangents(type_facet); + const Array & sigma_c_eff = mat_cohesive.getInsertionTraction(type_cohesive); + + Vector normal_stress(spatial_dimension); + + const Array > & coh_element_to_facet + = mesh_facets.getElementToSubelement(type_facet); + + Array::iterator > quad_facet_stress + = stress_facets.begin(spatial_dimension, spatial_dimension); + + Array::const_iterator > quad_normal + = normals.begin(spatial_dimension); + + Array::const_iterator > quad_tangents + = tangents.begin(tangents.getNbComponent()); + + for (UInt f = 0; f < nb_facet; ++f) { + const Element & cohesive_element = coh_element_to_facet(f)[1]; + + for (UInt q = 0; q < nb_quad_per_facet; ++q, ++quad_facet_stress, + ++quad_normal, ++quad_tangents) { + if (cohesive_element == ElementNull) continue; + + normal_stress.mul(*quad_facet_stress, *quad_normal); + + Real normal_contrib = normal_stress.dot(*quad_normal); + + Real first_tangent_contrib = 0; + + for (UInt dim = 0; dim < spatial_dimension; ++dim) + first_tangent_contrib += normal_stress(dim) * (*quad_tangents)(dim); + + Real second_tangent_contrib = 0; + + for (UInt dim = 0; dim < spatial_dimension; ++dim) + second_tangent_contrib + += normal_stress(dim) * (*quad_tangents)(dim + spatial_dimension); + + Real tangent_contrib = std::sqrt(first_tangent_contrib * first_tangent_contrib + + second_tangent_contrib * second_tangent_contrib); + + normal_contrib = std::max(0., normal_contrib); + + Real effective_norm = std::sqrt(normal_contrib * normal_contrib + + tangent_contrib * tangent_contrib / beta / beta); + + if (effective_norm < sigma_c) continue; + + if (!Math::are_float_equal(effective_norm, + sigma_c_eff(cohesive_element.element + * nb_quad_per_facet + q))) { + std::cout << "Insertion tractions do not match" << std::endl; + finalize(); + return EXIT_FAILURE; + } + } + } + finalize(); + return EXIT_SUCCESS; +} diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/test_cohesive_parallel_extrinsic_tetrahedron.sh b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/test_cohesive_parallel_extrinsic_tetrahedron.sh new file mode 100755 index 000000000..02f0ea54b --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/test_cohesive_parallel_extrinsic_tetrahedron.sh @@ -0,0 +1,6 @@ +#!/bin/bash + +# rm -r paraview +# mkdir paraview + +mpirun -np 8 ./test_cohesive_parallel_extrinsic_tetrahedron diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/test_cohesive_parallel_extrinsic_tetrahedron_displacement.cc b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/test_cohesive_parallel_extrinsic_tetrahedron_displacement.cc new file mode 100644 index 000000000..2abbd0b12 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/test_cohesive_parallel_extrinsic_tetrahedron_displacement.cc @@ -0,0 +1,395 @@ +/** + * @file test_cohesive_parallel_extrinsic_tetrahedron_displacement.cc + * + * @author Marco Vocialta + * + * + * @brief Displacement test for 3D cohesive elements + * + * @section LICENSE + * + * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + */ + +/* -------------------------------------------------------------------------- */ + +/* -------------------------------------------------------------------------- */ +#include "solid_mechanics_model_cohesive.hh" +#include "dumper_paraview.hh" +#include "material_cohesive_linear.hh" + +#ifdef AKANTU_USE_IOHELPER +# include "dumper_paraview.hh" +#endif + +/* -------------------------------------------------------------------------- */ +using namespace akantu; + +bool checkDisplacement(SolidMechanicsModelCohesive & model, + ElementType type, + std::ofstream & error_output, + UInt step, + bool barycenters); + +int main(int argc, char *argv[]) { + initialize("material.dat", argc, argv); + + debug::setDebugLevel(dblWarning); + + const UInt max_steps = 500; + Math::setTolerance(1.e-12); + + UInt spatial_dimension = 3; + ElementType type = _tetrahedron_10; + + Mesh mesh(spatial_dimension); + + StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + Int psize = comm.getNbProc(); + Int prank = comm.whoAmI(); + + akantu::MeshPartition * partition = NULL; + if(prank == 0) { + // Read the mesh + mesh.read("tetrahedron.msh"); + + /// partition the mesh + partition = new MeshPartitionScotch(mesh, spatial_dimension); + // debug::setDebugLevel(dblDump); + partition->partitionate(psize); + // debug::setDebugLevel(dblWarning); + } + + SolidMechanicsModelCohesive model(mesh); + + model.initParallel(partition, NULL, true); + + // debug::setDebugLevel(dblDump); + // std::cout << mesh << std::endl; + // debug::setDebugLevel(dblWarning); + + model.initFull(SolidMechanicsModelCohesiveOptions(_explicit_lumped_mass, true)); + + /* ------------------------------------------------------------------------ */ + /* Facet part */ + /* ------------------------------------------------------------------------ */ + + // Array limits(spatial_dimension, 2); + // limits(0, 0) = -0.01; + // limits(0, 1) = 0.01; + // limits(1, 0) = -100; + // limits(1, 1) = 100; + // limits(2, 0) = -100; + // limits(2, 1) = 100; + + // model.enableFacetsCheckOnArea(limits); + + /* ------------------------------------------------------------------------ */ + /* End of facet part */ + /* ------------------------------------------------------------------------ */ + + // debug::setDebugLevel(dblDump); + // std::cout << mesh_facets << std::endl; + // debug::setDebugLevel(dblWarning); + + Real time_step = model.getStableTimeStep()*0.1; + model.setTimeStep(time_step); + std::cout << "Time step: " << time_step << std::endl; + + model.assembleMassLumped(); + + + Array & position = mesh.getNodes(); + Array & velocity = model.getVelocity(); + Array & boundary = model.getBlockedDOFs(); + Array & displacement = model.getDisplacement(); + // const Array & residual = model.getResidual(); + + UInt nb_nodes = mesh.getNbNodes(); + + /// boundary conditions + for (UInt n = 0; n < nb_nodes; ++n) { + if (position(n, 0) > 0.99 || position(n, 0) < -0.99) { + for (UInt dim = 0; dim < spatial_dimension; ++dim) { + boundary(n, dim) = true; + } + } + + if (position(n, 0) > 0.99 || position(n, 0) < -0.99) { + for (UInt dim = 0; dim < spatial_dimension; ++dim) { + boundary(n, dim) = true; + } + } + } + +// #if defined (AKANTU_DEBUG_TOOLS) +// Vector facet_center(spatial_dimension); +// facet_center(0) = 0; +// facet_center(1) = -0.16666667; +// facet_center(2) = 0.5; + +// debug::element_manager.setMesh(mesh); +// debug::element_manager.addModule(debug::_dm_material_cohesive); +// debug::element_manager.addModule(debug::_dm_debug_tools); +// //debug::element_manager.addModule(debug::_dm_integrator); +// #endif + + /// initial conditions + Real loading_rate = 1; + Real disp_update = loading_rate * time_step; + for (UInt n = 0; n < nb_nodes; ++n) { + velocity(n, 0) = loading_rate * position(n, 0); + velocity(n, 1) = loading_rate * position(n, 0); + } + + model.synchronizeBoundaries(); + model.updateResidual(); + + std::stringstream paraview_output; + paraview_output << "extrinsic_parallel_tetrahedron_" << psize; + + model.setBaseName(paraview_output.str()); + model.addDumpFieldVector("displacement"); + model.addDumpFieldVector("velocity" ); + model.addDumpFieldVector("acceleration"); + model.addDumpFieldVector("residual" ); + model.addDumpFieldTensor("stress"); + model.addDumpFieldTensor("grad_u"); + model.addDumpField("partitions"); + // model.getDumper().getDumper().setMode(iohelper::BASE64); + model.dump(); + + + model.setBaseNameToDumper("cohesive elements", + paraview_output.str()+"_cohesive_elements"); + model.addDumpFieldVectorToDumper("cohesive elements", "displacement"); + model.addDumpFieldToDumper("cohesive elements", "damage"); + model.dump("cohesive elements"); + + + std::stringstream error_stream; + error_stream << "error" << ".csv"; + std::ofstream error_output; + error_output.open(error_stream.str().c_str()); + error_output << "# Step, Average, Max, Min" << std::endl; + + if (checkDisplacement(model, type, error_output, 0, true)) {} + + /// Main loop + for (UInt s = 1; s <= max_steps; ++s) { + + /// update displacement on extreme nodes + for (UInt n = 0; n < mesh.getNbNodes(); ++n) { + if (position(n, 0) > 0.99 || position(n, 0) < -0.99) { + displacement(n, 0) += disp_update * position(n, 0); + displacement(n, 1) += disp_update * position(n, 0); + } + } + + model.checkCohesiveStress(); + + model.solveStep(); + + if(s % 100 == 0) { + if(prank == 0) std::cout << "passing step " << s << "/" << max_steps << std::endl; + } + } + + model.dump(); + model.dump("cohesive elements"); + + if (!checkDisplacement(model, type, error_output, max_steps, false)) { + finalize(); + return EXIT_FAILURE; + } + + finalize(); + return EXIT_SUCCESS; +} + +bool checkDisplacement(SolidMechanicsModelCohesive & model, + ElementType type, + std::ofstream & error_output, + UInt step, + bool barycenters) { + + Mesh & mesh = model.getMesh(); + UInt spatial_dimension = mesh.getSpatialDimension(); + const Array & connectivity = mesh.getConnectivity(type); + const Array & displacement = model.getDisplacement(); + UInt nb_element = mesh.getNbElement(type); + UInt nb_nodes_per_elem = Mesh::getNbNodesPerElement(type); + + StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + Int psize = comm.getNbProc(); + Int prank = comm.whoAmI(); + + if (psize == 1) { + std::stringstream displacement_file; + displacement_file << "displacement/displacement_" + << std::setfill('0') << std::setw(6) + << step; + std::ofstream displacement_output; + displacement_output.open(displacement_file.str().c_str()); + + for (UInt el = 0; el < nb_element; ++el) { + for (UInt n = 0; n < nb_nodes_per_elem; ++n) { + UInt node = connectivity(el, n); + + for (UInt dim = 0; dim < spatial_dimension; ++dim) { + displacement_output << std::setprecision(15) + << displacement(node, dim) << " "; + } + displacement_output << std::endl; + } + } + + displacement_output.close(); + + if (barycenters) { + std::stringstream barycenter_file; + barycenter_file << "displacement/barycenters"; + std::ofstream barycenter_output; + barycenter_output.open(barycenter_file.str().c_str()); + + Element element(type, 0); + Vector bary(spatial_dimension); + + for (UInt el = 0; el < nb_element; ++el) { + element.element = el; + mesh.getBarycenter(element, bary); + + for (UInt dim = 0; dim < spatial_dimension; ++dim) { + barycenter_output << std::setprecision(15) + << bary(dim) << " "; + } + barycenter_output << std::endl; + } + + barycenter_output.close(); + } + } + else { + + if (barycenters) return true; + + /// read data + std::stringstream displacement_file; + displacement_file << "displacement/displacement_" + << std::setfill('0') << std::setw(6) + << step; + std::ifstream displacement_input; + displacement_input.open(displacement_file.str().c_str()); + + Array displacement_serial(0, spatial_dimension); + Vector disp_tmp(spatial_dimension); + + while (displacement_input.good()) { + for (UInt i = 0; i < spatial_dimension; ++i) + displacement_input >> disp_tmp(i); + + displacement_serial.push_back(disp_tmp); + } + + std::stringstream barycenter_file; + barycenter_file << "displacement/barycenters"; + std::ifstream barycenter_input; + barycenter_input.open(barycenter_file.str().c_str()); + + Array barycenter_serial(0, spatial_dimension); + + while (barycenter_input.good()) { + for (UInt dim = 0; dim < spatial_dimension; ++dim) + barycenter_input >> disp_tmp(dim); + + barycenter_serial.push_back(disp_tmp); + } + + Element element(type, 0); + Vector bary(spatial_dimension); + + Array::iterator > it; + Array::iterator > begin + = barycenter_serial.begin(spatial_dimension); + Array::iterator > end + = barycenter_serial.end(spatial_dimension); + + Array::const_iterator > disp_it; + Array::iterator > disp_serial_it; + + Vector difference(spatial_dimension); + Array error; + + /// compute error + for (UInt el = 0; el < nb_element; ++el) { + element.element = el; + mesh.getBarycenter(element, bary); + + /// find element + for (it = begin; it != end; ++it) { + UInt matched_dim = 0; + + while (matched_dim < spatial_dimension && + Math::are_float_equal(bary(matched_dim), (*it)(matched_dim))) + ++matched_dim; + + if (matched_dim == spatial_dimension) break; + } + + if (it == end) { + std::cout << "Element barycenter not found!" << std::endl; + return false; + } + + UInt matched_el = it - begin; + + disp_serial_it = displacement_serial.begin(spatial_dimension) + + matched_el * nb_nodes_per_elem; + + for (UInt n = 0; n < nb_nodes_per_elem; ++n, ++disp_serial_it) { + UInt node = connectivity(el, n); + if (!mesh.isLocalOrMasterNode(node)) continue; + + disp_it = displacement.begin(spatial_dimension) + node; + + difference = *disp_it; + difference -= *disp_serial_it; + + error.push_back(difference.norm()); + } + } + + /// compute average error + Real average_error = std::accumulate(error.begin(), error.end(), 0.); + comm.allReduce(&average_error, 1, _so_sum); + + UInt error_size = error.getSize(); + comm.allReduce(&error_size, 1, _so_sum); + + average_error /= error_size; + + /// compute maximum and minimum + Real max_error = *std::max_element(error.begin(), error.end()); + comm.allReduce(&max_error, 1, _so_max); + + Real min_error = *std::min_element(error.begin(), error.end()); + comm.allReduce(&min_error, 1, _so_min); + + /// output data + if (prank == 0) { + error_output << step << ", " + << average_error << ", " + << max_error << ", " + << min_error << std::endl; + } + + if (max_error > 1.e-9) { + std::cout << "Displacement error is too big!" << std::endl; + return false; + } + } + + return true; +} diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/test_cohesive_parallel_extrinsic_tetrahedron_displacement.sh b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/test_cohesive_parallel_extrinsic_tetrahedron_displacement.sh new file mode 100755 index 000000000..e4319e555 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/test_cohesive_parallel_extrinsic_tetrahedron_displacement.sh @@ -0,0 +1,10 @@ +#! /bin/bash + +rm -fr paraview +mkdir paraview + +rm -fr displacement +mkdir displacement + +./test_cohesive_parallel_extrinsic_tetrahedron_displacement +mpirun -np 8 ./test_cohesive_parallel_extrinsic_tetrahedron_displacement diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/tetrahedron.geo b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/tetrahedron.geo new file mode 100644 index 000000000..3e0a4c4fa --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic/tetrahedron.geo @@ -0,0 +1,13 @@ +lc = 1; +h = 1; + +Point(1) = {lc, lc, lc, h}; + +line[] = Extrude{-2*lc,0,0}{ Point{1}; }; +surface[] = Extrude{0,-2*lc,0}{ Line{line[1]}; }; +volume[] = Extrude{0,0,-2*lc}{ Surface{surface[1]}; }; + +Physical Volume(1) = {volume[1]}; + +Transfinite Surface "*"; +Transfinite Volume "*"; diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic_IG_TG/CMakeLists.txt b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic_IG_TG/CMakeLists.txt new file mode 100644 index 000000000..abb81ae57 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic_IG_TG/CMakeLists.txt @@ -0,0 +1,23 @@ +#=============================================================================== +# @file CMakeLists.txt +# +# @author Seyedeh Mohadeseh Taheri Mousavi +# @author Marco Vocialta +# +# +# @brief Configuration for test_cohesive_parallel_extrinsic_IG_TG +# +# @section LICENSE +# +# Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) +# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) +# +# @section DESCRIPTION +# +#=============================================================================== + +register_test(test_cohesive_parallel_extrinsic_IG_TG + SOURCES test_cohesive_parallel_extrinsic_IG_TG.cc + PACKAGE parallel_cohesive_element + FILES_TO_COPY material.dat square.msh + DIRECTORIES_TO_CREATE paraview) diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic_IG_TG/material.dat b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic_IG_TG/material.dat new file mode 100644 index 000000000..9a274f61a --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic_IG_TG/material.dat @@ -0,0 +1,20 @@ +material elastic [ + name = steel + rho = 1e3 # density + E = 1e3 # young's modulus + nu = 0.0 # poisson's ratio +] + +material cohesive_linear [ + name = TG_cohesive + beta = 0 + G_c = 10 + sigma_c = 100 +] + +material cohesive_linear [ + name = IG_cohesive + beta = 0 + G_c = 10 + sigma_c = 20 +] diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic_IG_TG/square.msh b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic_IG_TG/square.msh new file mode 100644 index 000000000..e34221c8e --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic_IG_TG/square.msh @@ -0,0 +1,66 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +41 +1 1 1 0 +2 -1 1 0 +3 -1 -1 0 +4 1 -1 0 +5 2.081668171172169e-12 1 0 +6 0.5000000000009475 1 0 +7 -0.4999999999989593 1 0 +8 -1 2.081668171172169e-12 0 +9 -1 0.5000000000009475 0 +10 -1 -0.4999999999989593 0 +11 -2.081668171172169e-12 -1 0 +12 -0.5000000000009475 -1 0 +13 0.4999999999989593 -1 0 +14 1 -2.081668171172169e-12 0 +15 1 -0.5000000000009475 0 +16 1 0.4999999999989593 0 +17 -2.888892828660043e-14 2.890280607440824e-14 0 +18 0.4999999999993061 -0.500000000000694 0 +19 0.500000000000694 0.4999999999993061 0 +20 -0.4999999999993061 0.500000000000694 0 +21 -0.5000000000006987 -0.4999999999993013 0 +22 0.7499999999996531 -0.7500000000003471 0 +23 0.2499999999986122 -0.7500000000003471 0 +24 0.7500000000003471 0.7499999999996531 0 +25 0.7500000000003471 0.2499999999986122 0 +26 -0.7499999999996531 0.7500000000003471 0 +27 -0.2499999999986122 0.7500000000003471 0 +28 -0.7500000000003493 -0.7499999999996506 0 +29 -0.7500000000003493 -0.2499999999986098 0 +30 1.026389621442784e-12 0.5000000000000144 0 +31 0.2500000000003326 0.2499999999996675 0 +32 0.2500000000013878 0.7499999999996531 0 +33 0.4999999999999856 -1.02638268254888e-12 0 +34 0.2499999999996386 -0.2500000000003325 0 +35 0.7499999999996531 -0.2500000000013878 0 +36 -0.2499999999996675 0.2500000000003614 0 +37 -0.5000000000000144 1.055285488623288e-12 0 +38 -0.7499999999996531 0.2500000000013878 0 +39 -1.055278549729385e-12 -0.4999999999999855 0 +40 -0.2500000000003638 -0.2499999999996362 0 +41 -0.2500000000013902 -0.7499999999996507 0 +$EndNodes +$Elements +16 +1 9 2 1 1 4 18 11 22 23 13 +2 9 2 2 2 1 19 14 24 25 16 +3 9 2 2 2 2 20 5 26 27 7 +4 9 2 1 1 3 21 8 28 29 10 +5 9 2 2 2 5 17 19 30 31 32 +6 9 2 1 1 14 17 18 33 34 35 +7 9 2 2 2 14 19 17 25 31 33 +8 9 2 2 2 5 20 17 27 36 30 +9 9 2 2 2 8 17 20 37 36 38 +10 9 2 1 1 11 18 17 23 34 39 +11 9 2 1 1 8 21 17 29 40 37 +12 9 2 1 1 11 17 21 39 40 41 +13 9 2 1 1 4 14 18 15 35 22 +14 9 2 2 2 2 8 20 9 38 26 +15 9 2 1 1 3 11 21 12 41 28 +16 9 2 2 2 1 5 19 6 32 24 +$EndElements diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic_IG_TG/test_cohesive_parallel_extrinsic_IG_TG.cc b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic_IG_TG/test_cohesive_parallel_extrinsic_IG_TG.cc new file mode 100644 index 000000000..84d56d478 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic_IG_TG/test_cohesive_parallel_extrinsic_IG_TG.cc @@ -0,0 +1,295 @@ +/** + * @file test_cohesive_parallel_extrinsic_IG_TG.cc + * + * @author Seyedeh Mohadeseh Taheri Mousavi + * @author Marco Vocialta + * + * + * @brief Test for considering different cohesive properties for + * intergranular (IG) and transgranular (TG) fractures in extrinsic + * cohesive elements + * + * @section LICENSE + * + * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + */ + +/* -------------------------------------------------------------------------- */ + +#include +#include +#include + + +/* -------------------------------------------------------------------------- */ +#include "solid_mechanics_model_cohesive.hh" +#include "material_cohesive_linear.hh" +/* -------------------------------------------------------------------------- */ + +using namespace akantu; + +class MultiGrainMaterialSelector : public DefaultMaterialCohesiveSelector { +public: + MultiGrainMaterialSelector(const SolidMechanicsModelCohesive & model, const ID & transgranular_id, const ID & intergranular_id) : + DefaultMaterialCohesiveSelector(model), + transgranular_id(transgranular_id), + intergranular_id(intergranular_id), + model(model), + mesh(model.getMesh()), + mesh_facets(model.getMeshFacets()), + spatial_dimension(model.getSpatialDimension()), + nb_IG(0), nb_TG(0) { + } + + UInt operator()(const Element & element) { + if(mesh_facets.getSpatialDimension(element.type) == (spatial_dimension - 1)) { + const std::vector & element_to_subelement = mesh_facets.getElementToSubelement(element.type, element.ghost_type)(element.element); + + const Element & el1 = element_to_subelement[0]; + const Element & el2 = element_to_subelement[1]; + + UInt grain_id1 = mesh.getData("tag_0", el1.type, el1.ghost_type)(el1.element); + if(el2 != ElementNull) { + UInt grain_id2 = mesh.getData("tag_0", el2.type, el2.ghost_type)(el2.element); + if (grain_id1 == grain_id2){ + //transgranular = 0 indicator + nb_TG++; + return model.getMaterialIndex(transgranular_id); + } else { + //intergranular = 1 indicator + nb_IG++; + return model.getMaterialIndex(intergranular_id); + } + } else { + //transgranular = 0 indicator + nb_TG++; + return model.getMaterialIndex(transgranular_id); + } + } else { + return DefaultMaterialCohesiveSelector::operator()(element); + } + } + +private: + ID transgranular_id, intergranular_id; + const SolidMechanicsModelCohesive & model; + const Mesh & mesh; + const Mesh & mesh_facets; + UInt spatial_dimension; + + UInt nb_IG; + UInt nb_TG; +}; + +/* -------------------------------------------------------------------------- */ +void limitInsertion(SolidMechanicsModelCohesive & model) { + + Real tolerance = 0.1; + + const Mesh & mesh = model.getMesh(); + const Mesh & mesh_facets = model.getMeshFacets(); + CohesiveElementInserter & inserter = model.getElementInserter(); + + UInt spatial_dimension = mesh.getSpatialDimension(); + Vector bary_facet(spatial_dimension); + + for (ghost_type_t::iterator gt = ghost_type_t::begin(); + gt != ghost_type_t::end(); + ++gt) { + GhostType ghost_type = *gt; + + Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, ghost_type); + Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1, ghost_type); + for(; it != end; ++it) { + ElementType type = *it; + Array & f_check = inserter.getCheckFacets(type, ghost_type); + UInt nb_facet = mesh_facets.getNbElement(type, ghost_type); + + for (UInt f = 0; f < nb_facet; ++f) { + if (f_check(f)) { + + mesh_facets.getBarycenter(f, type, bary_facet.storage(), ghost_type); + + if ( !(bary_facet(0) > -tolerance && bary_facet(0) < tolerance) && + !(bary_facet(1) > -tolerance && bary_facet(1) < tolerance) ) + f_check(f) = false; + } + } + } + } + + model.updateAutomaticInsertion(); +} + + +int main(int argc, char *argv[]) { + initialize("material.dat", argc, argv); + + debug::setDebugLevel(dblWarning); + + const UInt spatial_dimension = 2; + const UInt max_steps = 600; + + Mesh mesh(spatial_dimension); + + StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + Int psize = comm.getNbProc(); + Int prank = comm.whoAmI(); + + akantu::MeshPartition * partition = NULL; + if(prank == 0) { + mesh.read("square.msh"); + partition = new MeshPartitionScotch(mesh, spatial_dimension); + partition->partitionate(psize); + } + + SolidMechanicsModelCohesive model(mesh); + + /// model initialization + model.initParallel(partition, NULL, true); + + delete partition; + + MultiGrainMaterialSelector material_selector(model, "TG_cohesive", "IG_cohesive"); + model.setMaterialSelector(material_selector); + model.initFull(SolidMechanicsModelCohesiveOptions(_explicit_lumped_mass, true, false)); + + Real time_step = model.getStableTimeStep()*0.1; + model.setTimeStep(time_step); + // std::cout << "Time step: " << time_step << std::endl; + + limitInsertion(model); + + // std::cout << mesh << std::endl; + + Array & position = mesh.getNodes(); + Array & velocity = model.getVelocity(); + Array & boundary = model.getBlockedDOFs(); + Array & displacement = model.getDisplacement(); + // const Array & residual = model.getResidual(); + + UInt nb_nodes = mesh.getNbNodes(); + + /// boundary conditions + for (UInt n = 0; n < nb_nodes; ++n) { + if (position(n, 1) > 0.99|| position(n, 1) < -0.99) + boundary(n, 1) = true; + + if (position(n, 0) > 0.99 || position(n, 0) < -0.99) + boundary(n, 0) = true; + } + + model.synchronizeBoundaries(); + model.updateResidual(); + + model.setBaseName("extrinsic"); + model.addDumpFieldVector("displacement"); + model.addDumpField("velocity" ); + model.addDumpField("acceleration"); + model.addDumpField("residual" ); + model.addDumpField("stress"); + model.addDumpField("strain"); + model.addDumpField("partitions"); + + model.setBaseNameToDumper("cohesive elements", "extrinsic_cohesive"); + model.addDumpFieldVectorToDumper("cohesive elements", "displacement"); + model.addDumpFieldToDumper("cohesive elements", "damage"); + + model.dump(); + model.dump("cohesive elements"); + + /// initial conditions + Real loading_rate = 0.1; + // bar_height = 2 + Real VI = loading_rate * 2* 0.5; + for (UInt n = 0; n < nb_nodes; ++n) { + velocity(n, 1) = loading_rate * position(n, 1); + velocity(n, 0) = loading_rate * position(n, 0); + } + + // std::ofstream edis("edis.txt"); + // std::ofstream erev("erev.txt"); + + // Array & residual = model.getResidual(); + // model.dump(); + // const Array & stress = model.getMaterial(0).getStress(type); + Real dispy = 0; + // UInt nb_coh_elem = 0; + + /// Main loop + for (UInt s = 1; s <= max_steps; ++s) { + dispy += VI * time_step; + /// update displacement on extreme nodes + for (UInt n = 0; n < mesh.getNbNodes(); ++n) { + if (position(n, 1) > 0.99){ + displacement(n, 1) = dispy; + velocity(n,1) = VI;} + if (position(n, 1) < -0.99){ + displacement(n, 1) = -dispy; + velocity(n,1) = -VI;} + if (position(n, 0) > 0.99){ + displacement(n, 0) = dispy; + velocity(n,0) = VI;} + if (position(n, 0) < -0.99){ + displacement(n, 0) = -dispy; + velocity(n,0) = -VI;} + } + + model.checkCohesiveStress(); + + model.solveStep(); + + if(s % 10 == 0) { + if(prank == 0) + std::cout << "passing step " << s << "/" << max_steps << std::endl; + + // model.dump(); + // model.dump("cohesive elements"); + } + + // Real Ed = model.getEnergy("dissipated"); + + // edis << s << " " + // << Ed << std::endl; + + // erev << s << " " + // << Er << std::endl; + } + + model.dump(); + model.dump("cohesive elements"); + + // edis.close(); + // erev.close(); + + // mesh.write("mesh_final.msh"); + + Real Ed = model.getEnergy("dissipated"); + + Real Edt = 40; + + if(prank == 0) + std::cout << Ed << " " << Edt << std::endl; + + if (Ed < Edt * 0.99 || Ed > Edt * 1.01 || std::isnan(Ed)) { + if(prank == 0) + std::cout << "The dissipated energy is incorrect" << std::endl; + finalize(); + return EXIT_FAILURE; + } + + // for (UInt n = 0; n < position.getSize(); ++n) { + // for (UInt s = 0; s < spatial_dimension; ++s) { + // position(n, s) += displacement(n, s); + // } + // } + + + finalize(); + + if(prank == 0) + std::cout << "OK: test_cohesive_extrinsic_IG_TG was passed!" << std::endl; + return EXIT_SUCCESS; +} diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic_IG_TG/test_cohesive_parallel_extrinsic_IG_TG.sh b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic_IG_TG/test_cohesive_parallel_extrinsic_IG_TG.sh new file mode 100755 index 000000000..3fa1ff251 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_extrinsic_IG_TG/test_cohesive_parallel_extrinsic_IG_TG.sh @@ -0,0 +1,6 @@ +#!/bin/bash + +rm -r paraview +mkdir paraview + +mpirun -np 8 ./test_cohesive_parallel_extrinsic_IG_TG diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/2d_basic_interface.msh b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/2d_basic_interface.msh new file mode 100644 index 000000000..1c64afb98 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/2d_basic_interface.msh @@ -0,0 +1,46 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$PhysicalNames +3 +1 2 "coh_2" +1 3 "coh_1" +2 1 "bulk" +$EndPhysicalNames +$Nodes +13 +1 0.5 0 0 +2 0 0.5 0 +3 0 -0.5 0 +4 1 0 0 +5 1 0.5 0 +6 1 -0.5 0 +7 0 0 0 +8 0.499999999998694 0.5 0 +9 0.499999999998694 -0.5 0 +10 0.7499999999996735 0.25 0 +11 0.249999999999347 0.25 0 +12 0.7499999999996734 -0.25 0 +13 0.249999999999347 -0.25 0 +$EndNodes +$Elements +18 +1 1 2 2 4 1 4 +2 1 2 3 8 7 1 +3 2 2 1 9 5 8 10 +4 2 2 1 9 1 8 11 +5 2 2 1 9 2 7 11 +6 2 2 1 9 1 10 8 +7 2 2 1 9 1 11 7 +8 2 2 1 9 1 4 10 +9 2 2 1 9 4 5 10 +10 2 2 1 9 2 11 8 +11 2 2 1 11 6 12 9 +12 2 2 1 11 1 13 9 +13 2 2 1 11 3 13 7 +14 2 2 1 11 1 9 12 +15 2 2 1 11 1 7 13 +16 2 2 1 11 1 12 4 +17 2 2 1 11 4 12 6 +18 2 2 1 11 3 9 13 +$EndElements diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/3d_spherical_inclusion.geo b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/3d_spherical_inclusion.geo new file mode 100644 index 000000000..ad6b42ee9 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/3d_spherical_inclusion.geo @@ -0,0 +1,100 @@ +dx = 1; +w = 2; +radius = 0.2; + +Point(1) = {0,0,0,dx}; +Point(2) = {0,.3,0,dx}; +Point(3) = {0,-.3,0,dx}; +Point(4) = {w,0,0,dx}; +Point(5) = {w,.3,0,dx}; +Point(6) = {w,-.3,0,dx}; +Point(7) = {0,0,w,dx}; +Point(8) = {0,.3,w,dx}; +Point(9) = {0,-.3,w,dx}; +Point(10) = {w,0,w,dx}; +Point(11) = {w,.3,w,dx}; +Point(12) = {w,-.3,w,dx}; +Point(13) = {0.5*w,0,0.5*w,dx}; +Point(14) = {0.5*w+radius,0,0.5*w+radius,dx}; +Point(15) = {0.5*w-radius,0,0.5*w+radius,dx}; +Point(16) = {0.5*w+radius,0,0.5*w-radius,dx}; +Point(17) = {0.5*w-radius,0,0.5*w-radius,dx}; +Line(1) = {1, 2}; +Line(2) = {2, 5}; +Line(3) = {5, 4}; +Line(4) = {1, 4}; +Line(5) = {1, 3}; +Line(6) = {6, 4}; +Line(7) = {3, 6}; +Line(8) = {8, 11}; +Line(9) = {7, 10}; +Line(10) = {9, 12}; +Line(11) = {8, 7}; +Line(12) = {11, 10}; +Line(13) = {10, 12}; +Line(14) = {7, 9}; +Line(15) = {2, 8}; +Line(16) = {1, 7}; +Line(17) = {3, 9}; +Line(18) = {5, 11}; +Line(19) = {4, 10}; +Line(20) = {6, 12}; +Line Loop(21) = {18, 12, -19, -3}; +Plane Surface(22) = {21}; +Line Loop(23) = {19, 13, -20, 6}; +Plane Surface(24) = {23}; +Line Loop(25) = {11, -16, 1, 15}; +Plane Surface(26) = {25}; +Line Loop(27) = {14, -17, -5, 16}; +Plane Surface(28) = {27}; +Line Loop(29) = {8, 12, -9, -11}; +Plane Surface(30) = {29}; +Line Loop(31) = {9, 13, -10, -14}; +Plane Surface(32) = {31}; +Line Loop(33) = {2, 3, -4, 1}; +Plane Surface(34) = {33}; +Line Loop(35) = {7, 6, -4, 5}; +Plane Surface(36) = {35}; +Line Loop(37) = {18, -8, -15, 2}; +Plane Surface(38) = {37}; +Line Loop(39) = {7, 20, -10, -17}; +Plane Surface(40) = {39}; +Circle(41) = {16, 13, 14}; +Circle(42) = {14, 13, 15}; +Circle(43) = {15, 13, 17}; +Circle(44) = {17, 13, 16}; +Translate {0.5, 0, 0.5} { + Duplicata { Line{44, 43, 42, 41}; } +} +Translate {-0.5, 0, 0.5} { + Duplicata { Line{44, 43, 42, 41}; } +} +Translate {0.5, 0, -0.5} { + Duplicata { Line{44, 43, 42, 41}; } +} +Translate {-0.5, 0, -0.5} { + Duplicata { Line{44, 43, 42, 41}; } +} +Line Loop(61) = {9, -19, -4, 16}; +Line Loop(62) = {52, 51, 50, 49}; +Line Loop(63) = {57, 60, 59, 58}; +Line Loop(64) = {44, 41, 42, 43}; +Line Loop(65) = {55, 54, 53, 56}; +Line Loop(66) = {48, 47, 46, 45}; +Plane Surface(67) = {61, 62, 63, 64, 65, 66}; +Plane Surface(68) = {62}; +Plane Surface(69) = {63}; +Plane Surface(70) = {64}; +Plane Surface(71) = {66}; +Plane Surface(72) = {65}; +Surface Loop(73) = {26, 30, 38, 22, 34, 67, 68, 69, 70, 72, 71}; +Volume(74) = {73}; +Surface Loop(75) = {32, 24, 40, 36, 28, 67, 68, 69, 70, 72, 71}; +Volume(76) = {75}; +Physical Surface("interface") = {67}; +Physical Surface("coh1") = {68}; +Physical Surface("coh2") = {69}; +Physical Surface("coh3") = {70}; +Physical Surface("coh4") = {71}; +Physical Surface("coh5") = {72}; +Physical Volume("bulk") = {74, 76}; diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/CMakeLists.txt b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/CMakeLists.txt new file mode 100644 index 000000000..bb65b90ba --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/CMakeLists.txt @@ -0,0 +1,52 @@ +#=============================================================================== +# @file CMakeLists.txt +# +# @author Fabian Barras +# +# @date creation: Fri Aug 7 09:07:44 2015 +# +# @brief Tests parallel insertion of cohesive elements +# +# @section LICENSE +# +# Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) +# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) +# +# Akantu is free software: you can redistribute it and/or modify it under the +# terms of the GNU Lesser General Public License as published by the Free +# Software Foundation, either version 3 of the License, or (at your option) any +# later version. +# +# Akantu is distributed in the hope that it will be useful, but WITHOUT ANY +# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +# A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +# details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Akantu. If not, see . +# +# @section DESCRIPTION +# +#=============================================================================== + +add_mesh(3d_spherical_inclusion_parallel 3d_spherical_inclusion.geo 3 2) + +register_test(test_cohesive_parallel_insertion_along_physical_surfaces + SOURCES test_cohesive_parallel_insertion_along_physical_surfaces.cc + DEPENDS 3d_spherical_inclusion_parallel + FILES_TO_COPY input_file.dat + PACKAGE parallel_cohesive_element + ) + +register_test(test_cohesive_parallel_intrinsic_implicit_insertion + PARALLEL + POSTPROCESS verify_insertion.sh + PARALLEL_LEVEL "2 4" + SOURCES test_cohesive_parallel_intrinsic_implicit_insertion.cc + FILES_TO_COPY input_file_iii.dat + 2d_basic_interface.msh + output_dir_verified + PACKAGE parallel_cohesive_element + DIRECTORIES_TO_CREATE output_dir + ) + diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/input_file.dat b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/input_file.dat new file mode 100644 index 000000000..bbd74ab72 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/input_file.dat @@ -0,0 +1,67 @@ +material elastic [ + name = bulk + rho = 2500 + nu = 0.29 + E = 70e9 +# finite_deformation = 1 +] + +material cohesive_exponential [ + name = coh1 + sigma_c = 1.5e6 + beta = 1 + delta_c = 1e-4 + exponential_penalty = true + contact_tangent = 1.0 +] + +material cohesive_exponential [ + name = coh2 + sigma_c = 1.5e6 + beta = 1 + delta_c = 1e-4 + exponential_penalty = true + contact_tangent = 1.0 +] + +material cohesive_exponential [ + name = coh3 + sigma_c = 1.5e6 + beta = 1 + delta_c = 1e-4 + exponential_penalty = true + contact_tangent = 1.0 +] + +material cohesive_exponential [ + name = coh4 + sigma_c = 1.5e6 + beta = 1 + delta_c = 1e-4 + exponential_penalty = true + contact_tangent = 1.0 +] + +material cohesive_exponential [ + name = coh5 + sigma_c = 1.5e6 + beta = 1 + delta_c = 1e-4 + exponential_penalty = true + contact_tangent = 1.0 +] + +material cohesive_exponential [ + name = interface + sigma_c = 1.5e6 + beta = 1 + delta_c = 1e-4 + exponential_penalty = true + contact_tangent = 1.0 +] + +mesh parameters [ + + cohesive_surfaces = coh1,coh2,coh3,coh4,coh5,interface + +] \ No newline at end of file diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/input_file_iii.dat b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/input_file_iii.dat new file mode 100644 index 000000000..d02bdb7ce --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/input_file_iii.dat @@ -0,0 +1,32 @@ +seed = 5.0 +material elastic [ + name = head + rho = 2500 + nu = 0.29 + E = 70e9 +# finite_deformation = 1 +] + +material cohesive_exponential [ + name = coh_1 + sigma_c = 30e6 + beta = 1 + delta_c = 4e-7 + exponential_penalty = true + contact_tangent = 1.0 +] + +material cohesive_exponential [ + name = coh_2 + sigma_c = 60e6 + beta = 1 + delta_c = 8e-7 + exponential_penalty = true + contact_tangent = 1.0 +] + +mesh parameters [ + + cohesive_surfaces = coh_1,coh_2 + +] \ No newline at end of file diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/output_dir_verified/output_from_proc_0_out_of_2.out b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/output_dir_verified/output_from_proc_0_out_of_2.out new file mode 100644 index 000000000..bf4ad3d41 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/output_dir_verified/output_from_proc_0_out_of_2.out @@ -0,0 +1,289 @@ +The root processor read the mesh. +Only GMSH physical objects are created in the mesh. +Element type: _segment_2 ghost type: not_ghost +2 to partitionate between 2 processsors +0 1 +1 0 +Element type: _triangle_3 ghost type: not_ghost +16 to partitionate between 2 processsors +0 1 +1 0 +2 0 +3 1 +4 0 +5 1 +6 1 +7 0 +8 1 +9 0 +10 0 +11 1 +12 0 +13 1 +14 1 +15 0 +Nodes are also read and set with type -1 (normal node) +Number of nodes: 13 +Node # 0, x-coord: 0.5, y-coord: 0, of type: -1 +Node # 1, x-coord: 0, y-coord: 0.5, of type: -1 +Node # 2, x-coord: 0, y-coord: -0.5, of type: -1 +Node # 3, x-coord: 1, y-coord: 0, of type: -1 +Node # 4, x-coord: 1, y-coord: 0.5, of type: -1 +Node # 5, x-coord: 1, y-coord: -0.5, of type: -1 +Node # 6, x-coord: 0, y-coord: 0, of type: -1 +Node # 7, x-coord: 0.5, y-coord: 0.5, of type: -1 +Node # 8, x-coord: 0.5, y-coord: -0.5, of type: -1 +Node # 9, x-coord: 0.75, y-coord: 0.25, of type: -1 +Node # 10, x-coord: 0.25, y-coord: 0.25, of type: -1 +Node # 11, x-coord: 0.75, y-coord: -0.25, of type: -1 +Node # 12, x-coord: 0.25, y-coord: -0.25, of type: -1 + +Before initParallel(), non-root processors have empty Mesh object + +Element type: _segment_2, not_ghost: 2 in the mesh of processor 0 +Element no 0 0 3 +Element no 1 6 0 + +Element type: _triangle_3, not_ghost: 16 in the mesh of processor 0 +Element no 0 4 7 9 +Element no 1 0 7 10 +Element no 2 1 6 10 +Element no 3 0 9 7 +Element no 4 0 10 6 +Element no 5 0 3 9 +Element no 6 3 4 9 +Element no 7 1 10 7 +Element no 8 5 11 8 +Element no 9 0 12 8 +Element no 10 2 12 6 +Element no 11 0 8 11 +Element no 12 0 6 12 +Element no 13 0 11 3 +Element no 14 3 11 5 +Element no 15 2 8 12 +After initParallel(), Mesh object on each processor is a local partionated mesh containing ghost elements + +Element type: _segment_2, not_ghost: 1 in the mesh of processor 0 +Element no 0 0 1 + +Element type: _triangle_3, not_ghost: 8 in the mesh of processor 0 +Element no 0 1 3 4 +Element no 1 5 0 4 +Element no 2 1 4 0 +Element no 3 5 4 3 +Element no 4 1 6 7 +Element no 5 8 6 0 +Element no 6 1 0 6 +Element no 7 8 7 6 + +Element type: _segment_2, ghost: 1 in the mesh of processor 0 +Element no 0 1 2 + +Element type: _triangle_3, ghost: 6 in the mesh of processor 0 +Element no 0 9 3 10 +Element no 1 1 10 3 +Element no 2 1 2 10 +Element no 3 11 12 7 +Element no 4 1 7 12 +Element no 5 1 12 2 +Nodes are also partionated and new node types are defined: +Number of nodes: 13 +Node # 0, x-coord: 0, y-coord: 0, of type: -1 +Node # 1, x-coord: 0.5, y-coord: 0, of type: -2 +Node # 2, x-coord: 1, y-coord: 0, of type: -3 +Node # 3, x-coord: 0.5, y-coord: 0.5, of type: -2 +Node # 4, x-coord: 0.25, y-coord: 0.25, of type: -1 +Node # 5, x-coord: 0, y-coord: 0.5, of type: -1 +Node # 6, x-coord: 0.25, y-coord: -0.25, of type: -1 +Node # 7, x-coord: 0.5, y-coord: -0.5, of type: -2 +Node # 8, x-coord: 0, y-coord: -0.5, of type: -1 +Node # 9, x-coord: 1, y-coord: 0.5, of type: -3 +Node # 10, x-coord: 0.75, y-coord: 0.25, of type: -3 +Node # 11, x-coord: 1, y-coord: -0.5, of type: -3 +Node # 12, x-coord: 0.75, y-coord: -0.25, of type: -3 + +-3: pure ghost node -> not a local node +-2: master node -> node shared with other processor(s) -> local and global node +>0: slave node -> -> node shared with other processor(s) -> only local node (its id is the rank of the processor owning the master node) +Each local node has a corresponding global id used during assembly: +Global nodes ID: +0 6 +1 0 +2 3 +3 7 +4 10 +5 1 +6 12 +7 8 +8 2 +9 4 +10 9 +11 5 +12 11 + +Within cohesive element model, initParallel() creates a second Mesh object usually called mesh_facet +This Mesh object contains all sub-dimensional elements where potential cohesive element can be inserted + +Element type: _point_1, not_ghost: 8 in the mesh of processor 0 +Element no 0 1 +Element no 1 3 +Element no 2 4 +Element no 3 5 +Element no 4 0 +Element no 5 6 +Element no 6 7 +Element no 7 8 + +Element type: _segment_2, not_ghost: 15 in the mesh of processor 0 +Element no 0 1 3 +Element no 1 3 4 +Element no 2 4 1 +Element no 3 5 0 +Element no 4 0 4 +Element no 5 4 5 +Element no 6 0 1 +Element no 7 3 5 +Element no 8 1 6 +Element no 9 6 7 +Element no 10 7 1 +Element no 11 8 6 +Element no 12 6 0 +Element no 13 0 8 +Element no 14 8 7 + +Element type: _point_1, ghost: 0 in the mesh of processor 0 + +Element type: _segment_2, ghost: 7 in the mesh of processor 0 +Element no 0 9 3 +Element no 1 3 10 +Element no 2 1 10 +Element no 3 1 2 +Element no 4 12 7 +Element no 5 7 11 +Element no 6 12 1 +The distributed synchronizer of solid mechanics model is used to synchronize fields with ghost element: +From processor 0 to processor 1 + Sending element(s): +Element [_triangle_3, 0, not_ghost] +Element [_triangle_3, 2, not_ghost] +Element [_triangle_3, 3, not_ghost] +Element [_triangle_3, 4, not_ghost] +Element [_triangle_3, 6, not_ghost] +Element [_triangle_3, 7, not_ghost] + Receiving element(s): +Element [_triangle_3, 0, ghost] +Element [_triangle_3, 1, ghost] +Element [_triangle_3, 2, ghost] +Element [_triangle_3, 3, ghost] +Element [_triangle_3, 4, ghost] +Element [_triangle_3, 5, ghost] + +In case of insertion along physical objects, cohesive elements are created during initFull() +Elements list after insertion + +Element type: _cohesive_2d_4, not_ghost: 1 in the mesh of processor 0 +Element no 0 0 1 14 13 + +Element type: _segment_2, not_ghost: 1 in the mesh of processor 0 +Element no 0 0 1 + +Element type: _triangle_3, not_ghost: 8 in the mesh of processor 0 +Element no 0 1 3 4 +Element no 1 5 0 4 +Element no 2 1 4 0 +Element no 3 5 4 3 +Element no 4 13 6 7 +Element no 5 8 6 14 +Element no 6 13 14 6 +Element no 7 8 7 6 + +Element type: _cohesive_2d_4, ghost: 1 in the mesh of processor 0 +Element no 0 1 2 13 2 + +Element type: _segment_2, ghost: 1 in the mesh of processor 0 +Element no 0 1 2 + +Element type: _triangle_3, ghost: 6 in the mesh of processor 0 +Element no 0 9 3 10 +Element no 1 1 10 3 +Element no 2 1 2 10 +Element no 3 11 12 7 +Element no 4 13 7 12 +Element no 5 13 12 2 +Node list after insertion: (Total number of nodes 15) +Number of nodes: 15 +Node # 0, x-coord: 0, y-coord: 0, of type: -1 +Node # 1, x-coord: 0.5, y-coord: 0, of type: -2 +Node # 2, x-coord: 1, y-coord: 0, of type: -3 +Node # 3, x-coord: 0.5, y-coord: 0.5, of type: -2 +Node # 4, x-coord: 0.25, y-coord: 0.25, of type: -1 +Node # 5, x-coord: 0, y-coord: 0.5, of type: -1 +Node # 6, x-coord: 0.25, y-coord: -0.25, of type: -1 +Node # 7, x-coord: 0.5, y-coord: -0.5, of type: -2 +Node # 8, x-coord: 0, y-coord: -0.5, of type: -1 +Node # 9, x-coord: 1, y-coord: 0.5, of type: -3 +Node # 10, x-coord: 0.75, y-coord: 0.25, of type: -3 +Node # 11, x-coord: 1, y-coord: -0.5, of type: -3 +Node # 12, x-coord: 0.75, y-coord: -0.25, of type: -3 +Node # 13, x-coord: 0.5, y-coord: 0, of type: -2 +Node # 14, x-coord: 0, y-coord: 0, of type: -1 + +Node global ids after insertion: (Total number of nodes 16) +Global nodes ID: +0 6 +1 0 +2 3 +3 7 +4 10 +5 1 +6 12 +7 8 +8 2 +9 4 +10 9 +11 5 +12 11 +13 13 +14 14 + +Solid mechanics model cohesive has its own distributed synchronizer to handle ghost cohesive element: +From processor 0 to processor 1 + Sending element(s): +Element [_cohesive_2d_4, 0, not_ghost] + Receiving element(s): +Element [_cohesive_2d_4, 0, ghost] + +A synchronizer dedicated to degrees of freedom (DOFs) is used by the solver to build matrices in parallel: +This DOFSynchronizer is built based on nodes global id +Number of global dofs 32 for processor 0 +Local dof 0, global id: 12 +Local dof 1, global id: 13 +Local dof 2, global id: 0 +Local dof 3, global id: 1 +Local dof 4, global id: 6 +Local dof 5, global id: 7 +Local dof 6, global id: 14 +Local dof 7, global id: 15 +Local dof 8, global id: 20 +Local dof 9, global id: 21 +Local dof 10, global id: 2 +Local dof 11, global id: 3 +Local dof 12, global id: 24 +Local dof 13, global id: 25 +Local dof 14, global id: 16 +Local dof 15, global id: 17 +Local dof 16, global id: 4 +Local dof 17, global id: 5 +Local dof 18, global id: 8 +Local dof 19, global id: 9 +Local dof 20, global id: 18 +Local dof 21, global id: 19 +Local dof 22, global id: 10 +Local dof 23, global id: 11 +Local dof 24, global id: 22 +Local dof 25, global id: 23 +Local dof 26, global id: 26 +Local dof 27, global id: 27 +Local dof 28, global id: 28 +Local dof 29, global id: 29 + diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/output_dir_verified/output_from_proc_0_out_of_4.out b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/output_dir_verified/output_from_proc_0_out_of_4.out new file mode 100644 index 000000000..5f704014a --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/output_dir_verified/output_from_proc_0_out_of_4.out @@ -0,0 +1,289 @@ +The root processor read the mesh. +Only GMSH physical objects are created in the mesh. +Element type: _segment_2 ghost type: not_ghost +2 to partitionate between 4 processsors +0 2 +1 0 +Element type: _triangle_3 ghost type: not_ghost +16 to partitionate between 4 processsors +0 2 +1 0 +2 0 +3 2 +4 0 +5 2 +6 2 +7 0 +8 3 +9 1 +10 1 +11 3 +12 1 +13 3 +14 3 +15 1 +Nodes are also read and set with type -1 (normal node) +Number of nodes: 13 +Node # 0, x-coord: 0.5, y-coord: 0, of type: -1 +Node # 1, x-coord: 0, y-coord: 0.5, of type: -1 +Node # 2, x-coord: 0, y-coord: -0.5, of type: -1 +Node # 3, x-coord: 1, y-coord: 0, of type: -1 +Node # 4, x-coord: 1, y-coord: 0.5, of type: -1 +Node # 5, x-coord: 1, y-coord: -0.5, of type: -1 +Node # 6, x-coord: 0, y-coord: 0, of type: -1 +Node # 7, x-coord: 0.5, y-coord: 0.5, of type: -1 +Node # 8, x-coord: 0.5, y-coord: -0.5, of type: -1 +Node # 9, x-coord: 0.75, y-coord: 0.25, of type: -1 +Node # 10, x-coord: 0.25, y-coord: 0.25, of type: -1 +Node # 11, x-coord: 0.75, y-coord: -0.25, of type: -1 +Node # 12, x-coord: 0.25, y-coord: -0.25, of type: -1 + +Before initParallel(), non-root processors have empty Mesh object + +Element type: _segment_2, not_ghost: 2 in the mesh of processor 0 +Element no 0 0 3 +Element no 1 6 0 + +Element type: _triangle_3, not_ghost: 16 in the mesh of processor 0 +Element no 0 4 7 9 +Element no 1 0 7 10 +Element no 2 1 6 10 +Element no 3 0 9 7 +Element no 4 0 10 6 +Element no 5 0 3 9 +Element no 6 3 4 9 +Element no 7 1 10 7 +Element no 8 5 11 8 +Element no 9 0 12 8 +Element no 10 2 12 6 +Element no 11 0 8 11 +Element no 12 0 6 12 +Element no 13 0 11 3 +Element no 14 3 11 5 +Element no 15 2 8 12 +After initParallel(), Mesh object on each processor is a local partionated mesh containing ghost elements + +Element type: _segment_2, not_ghost: 1 in the mesh of processor 0 +Element no 0 0 1 + +Element type: _triangle_3, not_ghost: 4 in the mesh of processor 0 +Element no 0 1 3 4 +Element no 1 5 0 4 +Element no 2 1 4 0 +Element no 3 5 4 3 + +Element type: _segment_2, ghost: 1 in the mesh of processor 0 +Element no 0 1 2 + +Element type: _triangle_3, ghost: 8 in the mesh of processor 0 +Element no 0 6 3 7 +Element no 1 1 7 3 +Element no 2 1 2 7 +Element no 3 1 8 9 +Element no 4 10 8 0 +Element no 5 1 9 11 +Element no 6 1 0 8 +Element no 7 1 11 2 +Nodes are also partionated and new node types are defined: +Number of nodes: 12 +Node # 0, x-coord: 0, y-coord: 0, of type: -2 +Node # 1, x-coord: 0.5, y-coord: 0, of type: -2 +Node # 2, x-coord: 1, y-coord: 0, of type: -3 +Node # 3, x-coord: 0.5, y-coord: 0.5, of type: -2 +Node # 4, x-coord: 0.25, y-coord: 0.25, of type: -1 +Node # 5, x-coord: 0, y-coord: 0.5, of type: -1 +Node # 6, x-coord: 1, y-coord: 0.5, of type: -3 +Node # 7, x-coord: 0.75, y-coord: 0.25, of type: -3 +Node # 8, x-coord: 0.25, y-coord: -0.25, of type: -3 +Node # 9, x-coord: 0.5, y-coord: -0.5, of type: -3 +Node # 10, x-coord: 0, y-coord: -0.5, of type: -3 +Node # 11, x-coord: 0.75, y-coord: -0.25, of type: -3 + +-3: pure ghost node -> not a local node +-2: master node -> node shared with other processor(s) -> local and global node +>0: slave node -> -> node shared with other processor(s) -> only local node (its id is the rank of the processor owning the master node) +Each local node has a corresponding global id used during assembly: +Global nodes ID: +0 6 +1 0 +2 3 +3 7 +4 10 +5 1 +6 4 +7 9 +8 12 +9 8 +10 2 +11 11 + +Within cohesive element model, initParallel() creates a second Mesh object usually called mesh_facet +This Mesh object contains all sub-dimensional elements where potential cohesive element can be inserted + +Element type: _point_1, not_ghost: 5 in the mesh of processor 0 +Element no 0 1 +Element no 1 3 +Element no 2 4 +Element no 3 5 +Element no 4 0 + +Element type: _segment_2, not_ghost: 8 in the mesh of processor 0 +Element no 0 1 3 +Element no 1 3 4 +Element no 2 4 1 +Element no 3 5 0 +Element no 4 0 4 +Element no 5 4 5 +Element no 6 0 1 +Element no 7 3 5 + +Element type: _point_1, ghost: 0 in the mesh of processor 0 + +Element type: _segment_2, ghost: 9 in the mesh of processor 0 +Element no 0 6 3 +Element no 1 3 7 +Element no 2 1 7 +Element no 3 1 2 +Element no 4 1 8 +Element no 5 9 1 +Element no 6 8 0 +Element no 7 0 10 +Element no 8 11 1 +The distributed synchronizer of solid mechanics model is used to synchronize fields with ghost element: +From processor 0 to processor 1 + Sending element(s): +Element [_triangle_3, 0, not_ghost] +Element [_triangle_3, 1, not_ghost] +Element [_triangle_3, 2, not_ghost] + Receiving element(s): +Element [_triangle_3, 3, ghost] +Element [_triangle_3, 4, ghost] +Element [_triangle_3, 6, ghost] +From processor 0 to processor 2 + Sending element(s): +Element [_triangle_3, 0, not_ghost] +Element [_triangle_3, 2, not_ghost] +Element [_triangle_3, 3, not_ghost] + Receiving element(s): +Element [_triangle_3, 0, ghost] +Element [_triangle_3, 1, ghost] +Element [_triangle_3, 2, ghost] +From processor 0 to processor 3 + Sending element(s): +Element [_triangle_3, 0, not_ghost] +Element [_triangle_3, 2, not_ghost] + Receiving element(s): +Element [_triangle_3, 5, ghost] +Element [_triangle_3, 7, ghost] + +In case of insertion along physical objects, cohesive elements are created during initFull() +Elements list after insertion + +Element type: _cohesive_2d_4, not_ghost: 1 in the mesh of processor 0 +Element no 0 0 1 13 12 + +Element type: _segment_2, not_ghost: 1 in the mesh of processor 0 +Element no 0 0 1 + +Element type: _triangle_3, not_ghost: 4 in the mesh of processor 0 +Element no 0 1 3 4 +Element no 1 5 0 4 +Element no 2 1 4 0 +Element no 3 5 4 3 + +Element type: _cohesive_2d_4, ghost: 1 in the mesh of processor 0 +Element no 0 1 2 12 2 + +Element type: _segment_2, ghost: 1 in the mesh of processor 0 +Element no 0 1 2 + +Element type: _triangle_3, ghost: 8 in the mesh of processor 0 +Element no 0 6 3 7 +Element no 1 1 7 3 +Element no 2 1 2 7 +Element no 3 12 8 9 +Element no 4 10 8 13 +Element no 5 12 9 11 +Element no 6 12 13 8 +Element no 7 12 11 2 +Node list after insertion: (Total number of nodes 14) +Number of nodes: 14 +Node # 0, x-coord: 0, y-coord: 0, of type: -2 +Node # 1, x-coord: 0.5, y-coord: 0, of type: -2 +Node # 2, x-coord: 1, y-coord: 0, of type: -3 +Node # 3, x-coord: 0.5, y-coord: 0.5, of type: -2 +Node # 4, x-coord: 0.25, y-coord: 0.25, of type: -1 +Node # 5, x-coord: 0, y-coord: 0.5, of type: -1 +Node # 6, x-coord: 1, y-coord: 0.5, of type: -3 +Node # 7, x-coord: 0.75, y-coord: 0.25, of type: -3 +Node # 8, x-coord: 0.25, y-coord: -0.25, of type: -3 +Node # 9, x-coord: 0.5, y-coord: -0.5, of type: -3 +Node # 10, x-coord: 0, y-coord: -0.5, of type: -3 +Node # 11, x-coord: 0.75, y-coord: -0.25, of type: -3 +Node # 12, x-coord: 0.5, y-coord: 0, of type: -2 +Node # 13, x-coord: 0, y-coord: 0, of type: -2 + +Node global ids after insertion: (Total number of nodes 16) +Global nodes ID: +0 6 +1 0 +2 3 +3 7 +4 10 +5 1 +6 4 +7 9 +8 12 +9 8 +10 2 +11 11 +12 13 +13 14 + +Solid mechanics model cohesive has its own distributed synchronizer to handle ghost cohesive element: +From processor 0 to processor 1 + Sending element(s): +Element [_cohesive_2d_4, 0, not_ghost] + Receiving element(s): +From processor 0 to processor 2 + Sending element(s): +Element [_cohesive_2d_4, 0, not_ghost] + Receiving element(s): +Element [_cohesive_2d_4, 0, ghost] +From processor 0 to processor 3 + Sending element(s): +Element [_cohesive_2d_4, 0, not_ghost] + Receiving element(s): + +A synchronizer dedicated to degrees of freedom (DOFs) is used by the solver to build matrices in parallel: +This DOFSynchronizer is built based on nodes global id +Number of global dofs 32 for processor 0 +Local dof 0, global id: 12 +Local dof 1, global id: 13 +Local dof 2, global id: 0 +Local dof 3, global id: 1 +Local dof 4, global id: 6 +Local dof 5, global id: 7 +Local dof 6, global id: 14 +Local dof 7, global id: 15 +Local dof 8, global id: 20 +Local dof 9, global id: 21 +Local dof 10, global id: 2 +Local dof 11, global id: 3 +Local dof 12, global id: 8 +Local dof 13, global id: 9 +Local dof 14, global id: 18 +Local dof 15, global id: 19 +Local dof 16, global id: 24 +Local dof 17, global id: 25 +Local dof 18, global id: 16 +Local dof 19, global id: 17 +Local dof 20, global id: 4 +Local dof 21, global id: 5 +Local dof 22, global id: 22 +Local dof 23, global id: 23 +Local dof 24, global id: 26 +Local dof 25, global id: 27 +Local dof 26, global id: 28 +Local dof 27, global id: 29 + diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/output_dir_verified/output_from_proc_1_out_of_2.out b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/output_dir_verified/output_from_proc_1_out_of_2.out new file mode 100644 index 000000000..974c7e0b9 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/output_dir_verified/output_from_proc_1_out_of_2.out @@ -0,0 +1,227 @@ +Before initParallel(), non-root processors have empty Mesh object +After initParallel(), Mesh object on each processor is a local partionated mesh containing ghost elements + +Element type: _segment_2, not_ghost: 1 in the mesh of processor 1 +Element no 0 0 1 + +Element type: _triangle_3, not_ghost: 8 in the mesh of processor 1 +Element no 0 3 4 5 +Element no 1 0 5 4 +Element no 2 0 1 5 +Element no 3 1 3 5 +Element no 4 6 7 8 +Element no 5 0 8 7 +Element no 6 0 7 1 +Element no 7 1 7 6 + +Element type: _segment_2, ghost: 1 in the mesh of processor 1 +Element no 0 2 0 + +Element type: _triangle_3, ghost: 6 in the mesh of processor 1 +Element no 0 0 4 9 +Element no 1 0 9 2 +Element no 2 10 9 4 +Element no 3 0 11 8 +Element no 4 0 2 11 +Element no 5 12 8 11 +Nodes are also partionated and new node types are defined: +Number of nodes: 13 +Node # 0, x-coord: 0.5, y-coord: 0, of type: 0 +Node # 1, x-coord: 1, y-coord: 0, of type: -1 +Node # 2, x-coord: 0, y-coord: 0, of type: -3 +Node # 3, x-coord: 1, y-coord: 0.5, of type: -1 +Node # 4, x-coord: 0.5, y-coord: 0.5, of type: 0 +Node # 5, x-coord: 0.75, y-coord: 0.25, of type: -1 +Node # 6, x-coord: 1, y-coord: -0.5, of type: -1 +Node # 7, x-coord: 0.75, y-coord: -0.25, of type: -1 +Node # 8, x-coord: 0.5, y-coord: -0.5, of type: 0 +Node # 9, x-coord: 0.25, y-coord: 0.25, of type: -3 +Node # 10, x-coord: 0, y-coord: 0.5, of type: -3 +Node # 11, x-coord: 0.25, y-coord: -0.25, of type: -3 +Node # 12, x-coord: 0, y-coord: -0.5, of type: -3 + +-3: pure ghost node -> not a local node +-2: master node -> node shared with other processor(s) -> local and global node +>0: slave node -> -> node shared with other processor(s) -> only local node (its id is the rank of the processor owning the master node) +Each local node has a corresponding global id used during assembly: +Global nodes ID: +0 0 +1 3 +2 6 +3 4 +4 7 +5 9 +6 5 +7 11 +8 8 +9 10 +10 1 +11 12 +12 2 + +Within cohesive element model, initParallel() creates a second Mesh object usually called mesh_facet +This Mesh object contains all sub-dimensional elements where potential cohesive element can be inserted + +Element type: _point_1, not_ghost: 8 in the mesh of processor 1 +Element no 0 3 +Element no 1 4 +Element no 2 5 +Element no 3 0 +Element no 4 1 +Element no 5 6 +Element no 6 7 +Element no 7 8 + +Element type: _segment_2, not_ghost: 13 in the mesh of processor 1 +Element no 0 3 4 +Element no 1 4 5 +Element no 2 5 3 +Element no 3 0 5 +Element no 4 0 1 +Element no 5 1 5 +Element no 6 1 3 +Element no 7 6 7 +Element no 8 7 8 +Element no 9 8 6 +Element no 10 7 0 +Element no 11 7 1 +Element no 12 6 1 + +Element type: _point_1, ghost: 0 in the mesh of processor 1 + +Element type: _segment_2, ghost: 9 in the mesh of processor 1 +Element no 0 0 4 +Element no 1 8 0 +Element no 2 4 9 +Element no 3 9 0 +Element no 4 2 0 +Element no 5 4 10 +Element no 6 0 11 +Element no 7 11 8 +Element no 8 12 8 +The distributed synchronizer of solid mechanics model is used to synchronize fields with ghost element: +From processor 1 to processor 0 + Sending element(s): +Element [_triangle_3, 0, not_ghost] +Element [_triangle_3, 1, not_ghost] +Element [_triangle_3, 2, not_ghost] +Element [_triangle_3, 4, not_ghost] +Element [_triangle_3, 5, not_ghost] +Element [_triangle_3, 6, not_ghost] + Receiving element(s): +Element [_triangle_3, 0, ghost] +Element [_triangle_3, 1, ghost] +Element [_triangle_3, 2, ghost] +Element [_triangle_3, 3, ghost] +Element [_triangle_3, 4, ghost] +Element [_triangle_3, 5, ghost] + +In case of insertion along physical objects, cohesive elements are created during initFull() +Elements list after insertion + +Element type: _cohesive_2d_4, not_ghost: 1 in the mesh of processor 1 +Element no 0 0 1 13 14 + +Element type: _segment_2, not_ghost: 1 in the mesh of processor 1 +Element no 0 0 1 + +Element type: _triangle_3, not_ghost: 8 in the mesh of processor 1 +Element no 0 3 4 5 +Element no 1 0 5 4 +Element no 2 0 1 5 +Element no 3 1 3 5 +Element no 4 6 7 8 +Element no 5 13 8 7 +Element no 6 13 7 14 +Element no 7 14 7 6 + +Element type: _cohesive_2d_4, ghost: 1 in the mesh of processor 1 +Element no 0 2 0 2 13 + +Element type: _segment_2, ghost: 1 in the mesh of processor 1 +Element no 0 2 0 + +Element type: _triangle_3, ghost: 6 in the mesh of processor 1 +Element no 0 0 4 9 +Element no 1 0 9 2 +Element no 2 10 9 4 +Element no 3 13 11 8 +Element no 4 13 2 11 +Element no 5 12 8 11 +Node list after insertion: (Total number of nodes 15) +Number of nodes: 15 +Node # 0, x-coord: 0.5, y-coord: 0, of type: 0 +Node # 1, x-coord: 1, y-coord: 0, of type: -1 +Node # 2, x-coord: 0, y-coord: 0, of type: -3 +Node # 3, x-coord: 1, y-coord: 0.5, of type: -1 +Node # 4, x-coord: 0.5, y-coord: 0.5, of type: 0 +Node # 5, x-coord: 0.75, y-coord: 0.25, of type: -1 +Node # 6, x-coord: 1, y-coord: -0.5, of type: -1 +Node # 7, x-coord: 0.75, y-coord: -0.25, of type: -1 +Node # 8, x-coord: 0.5, y-coord: -0.5, of type: 0 +Node # 9, x-coord: 0.25, y-coord: 0.25, of type: -3 +Node # 10, x-coord: 0, y-coord: 0.5, of type: -3 +Node # 11, x-coord: 0.25, y-coord: -0.25, of type: -3 +Node # 12, x-coord: 0, y-coord: -0.5, of type: -3 +Node # 13, x-coord: 0.5, y-coord: 0, of type: 0 +Node # 14, x-coord: 1, y-coord: 0, of type: -1 + +Node global ids after insertion: (Total number of nodes 16) +Global nodes ID: +0 0 +1 3 +2 6 +3 4 +4 7 +5 9 +6 5 +7 11 +8 8 +9 10 +10 1 +11 12 +12 2 +13 13 +14 15 + +Solid mechanics model cohesive has its own distributed synchronizer to handle ghost cohesive element: +From processor 1 to processor 0 + Sending element(s): +Element [_cohesive_2d_4, 0, not_ghost] + Receiving element(s): +Element [_cohesive_2d_4, 0, ghost] + +A synchronizer dedicated to degrees of freedom (DOFs) is used by the solver to build matrices in parallel: +This DOFSynchronizer is built based on nodes global id +Number of global dofs 32 for processor 1 +Local dof 0, global id: 0 +Local dof 1, global id: 1 +Local dof 2, global id: 6 +Local dof 3, global id: 7 +Local dof 4, global id: 12 +Local dof 5, global id: 13 +Local dof 6, global id: 8 +Local dof 7, global id: 9 +Local dof 8, global id: 14 +Local dof 9, global id: 15 +Local dof 10, global id: 18 +Local dof 11, global id: 19 +Local dof 12, global id: 10 +Local dof 13, global id: 11 +Local dof 14, global id: 22 +Local dof 15, global id: 23 +Local dof 16, global id: 16 +Local dof 17, global id: 17 +Local dof 18, global id: 20 +Local dof 19, global id: 21 +Local dof 20, global id: 2 +Local dof 21, global id: 3 +Local dof 22, global id: 24 +Local dof 23, global id: 25 +Local dof 24, global id: 4 +Local dof 25, global id: 5 +Local dof 26, global id: 26 +Local dof 27, global id: 27 +Local dof 28, global id: 30 +Local dof 29, global id: 31 + diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/output_dir_verified/output_from_proc_1_out_of_4.out b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/output_dir_verified/output_from_proc_1_out_of_4.out new file mode 100644 index 000000000..96c8e4c20 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/output_dir_verified/output_from_proc_1_out_of_4.out @@ -0,0 +1,225 @@ +Before initParallel(), non-root processors have empty Mesh object +After initParallel(), Mesh object on each processor is a local partionated mesh containing ghost elements + +Element type: _segment_2, not_ghost: 0 in the mesh of processor 1 + +Element type: _triangle_3, not_ghost: 4 in the mesh of processor 1 +Element no 0 0 3 4 +Element no 1 5 3 2 +Element no 2 0 2 3 +Element no 3 5 4 3 + +Element type: _segment_2, ghost: 2 in the mesh of processor 1 +Element no 0 0 1 +Element no 1 2 0 + +Element type: _triangle_3, ghost: 8 in the mesh of processor 1 +Element no 0 0 6 7 +Element no 1 8 2 7 +Element no 2 0 9 6 +Element no 3 0 7 2 +Element no 4 0 1 9 +Element no 5 10 11 4 +Element no 6 0 4 11 +Element no 7 0 11 1 +Nodes are also partionated and new node types are defined: +Number of nodes: 12 +Node # 0, x-coord: 0.5, y-coord: 0, of type: 0 +Node # 1, x-coord: 1, y-coord: 0, of type: -3 +Node # 2, x-coord: 0, y-coord: 0, of type: 0 +Node # 3, x-coord: 0.25, y-coord: -0.25, of type: -1 +Node # 4, x-coord: 0.5, y-coord: -0.5, of type: -2 +Node # 5, x-coord: 0, y-coord: -0.5, of type: -1 +Node # 6, x-coord: 0.5, y-coord: 0.5, of type: -3 +Node # 7, x-coord: 0.25, y-coord: 0.25, of type: -3 +Node # 8, x-coord: 0, y-coord: 0.5, of type: -3 +Node # 9, x-coord: 0.75, y-coord: 0.25, of type: -3 +Node # 10, x-coord: 1, y-coord: -0.5, of type: -3 +Node # 11, x-coord: 0.75, y-coord: -0.25, of type: -3 + +-3: pure ghost node -> not a local node +-2: master node -> node shared with other processor(s) -> local and global node +>0: slave node -> -> node shared with other processor(s) -> only local node (its id is the rank of the processor owning the master node) +Each local node has a corresponding global id used during assembly: +Global nodes ID: +0 0 +1 3 +2 6 +3 12 +4 8 +5 2 +6 7 +7 10 +8 1 +9 9 +10 5 +11 11 + +Within cohesive element model, initParallel() creates a second Mesh object usually called mesh_facet +This Mesh object contains all sub-dimensional elements where potential cohesive element can be inserted + +Element type: _point_1, not_ghost: 5 in the mesh of processor 1 +Element no 0 0 +Element no 1 3 +Element no 2 4 +Element no 3 5 +Element no 4 2 + +Element type: _segment_2, not_ghost: 7 in the mesh of processor 1 +Element no 0 0 3 +Element no 1 3 4 +Element no 2 4 0 +Element no 3 5 3 +Element no 4 3 2 +Element no 5 2 5 +Element no 6 5 4 + +Element type: _point_1, ghost: 0 in the mesh of processor 1 + +Element type: _segment_2, ghost: 10 in the mesh of processor 1 +Element no 0 2 0 +Element no 1 0 6 +Element no 2 7 0 +Element no 3 8 2 +Element no 4 2 7 +Element no 5 0 9 +Element no 6 0 1 +Element no 7 11 4 +Element no 8 4 10 +Element no 9 11 0 +The distributed synchronizer of solid mechanics model is used to synchronize fields with ghost element: +From processor 1 to processor 0 + Sending element(s): +Element [_triangle_3, 0, not_ghost] +Element [_triangle_3, 1, not_ghost] +Element [_triangle_3, 2, not_ghost] + Receiving element(s): +Element [_triangle_3, 0, ghost] +Element [_triangle_3, 1, ghost] +Element [_triangle_3, 3, ghost] +From processor 1 to processor 2 + Sending element(s): +Element [_triangle_3, 0, not_ghost] +Element [_triangle_3, 2, not_ghost] + Receiving element(s): +Element [_triangle_3, 2, ghost] +Element [_triangle_3, 4, ghost] +From processor 1 to processor 3 + Sending element(s): +Element [_triangle_3, 0, not_ghost] +Element [_triangle_3, 2, not_ghost] +Element [_triangle_3, 3, not_ghost] + Receiving element(s): +Element [_triangle_3, 5, ghost] +Element [_triangle_3, 6, ghost] +Element [_triangle_3, 7, ghost] + +In case of insertion along physical objects, cohesive elements are created during initFull() +Elements list after insertion + +Element type: _cohesive_2d_4, not_ghost: 0 in the mesh of processor 1 + +Element type: _segment_2, not_ghost: 0 in the mesh of processor 1 + +Element type: _triangle_3, not_ghost: 4 in the mesh of processor 1 +Element no 0 13 3 4 +Element no 1 5 3 12 +Element no 2 13 12 3 +Element no 3 5 4 3 + +Element type: _cohesive_2d_4, ghost: 2 in the mesh of processor 1 +Element no 0 2 0 12 13 +Element no 1 0 1 13 1 + +Element type: _segment_2, ghost: 2 in the mesh of processor 1 +Element no 0 0 1 +Element no 1 2 0 + +Element type: _triangle_3, ghost: 8 in the mesh of processor 1 +Element no 0 0 6 7 +Element no 1 8 2 7 +Element no 2 0 9 6 +Element no 3 0 7 2 +Element no 4 0 1 9 +Element no 5 10 11 4 +Element no 6 13 4 11 +Element no 7 13 11 1 +Node list after insertion: (Total number of nodes 14) +Number of nodes: 14 +Node # 0, x-coord: 0.5, y-coord: 0, of type: 0 +Node # 1, x-coord: 1, y-coord: 0, of type: -3 +Node # 2, x-coord: 0, y-coord: 0, of type: 0 +Node # 3, x-coord: 0.25, y-coord: -0.25, of type: -1 +Node # 4, x-coord: 0.5, y-coord: -0.5, of type: -2 +Node # 5, x-coord: 0, y-coord: -0.5, of type: -1 +Node # 6, x-coord: 0.5, y-coord: 0.5, of type: -3 +Node # 7, x-coord: 0.25, y-coord: 0.25, of type: -3 +Node # 8, x-coord: 0, y-coord: 0.5, of type: -3 +Node # 9, x-coord: 0.75, y-coord: 0.25, of type: -3 +Node # 10, x-coord: 1, y-coord: -0.5, of type: -3 +Node # 11, x-coord: 0.75, y-coord: -0.25, of type: -3 +Node # 12, x-coord: 0, y-coord: 0, of type: 0 +Node # 13, x-coord: 0.5, y-coord: 0, of type: 0 + +Node global ids after insertion: (Total number of nodes 16) +Global nodes ID: +0 0 +1 3 +2 6 +3 12 +4 8 +5 2 +6 7 +7 10 +8 1 +9 9 +10 5 +11 11 +12 14 +13 13 + +Solid mechanics model cohesive has its own distributed synchronizer to handle ghost cohesive element: +From processor 1 to processor 0 + Sending element(s): + Receiving element(s): +Element [_cohesive_2d_4, 0, ghost] +From processor 1 to processor 2 + Sending element(s): + Receiving element(s): +Element [_cohesive_2d_4, 1, ghost] +From processor 1 to processor 3 + Sending element(s): + Receiving element(s): + +A synchronizer dedicated to degrees of freedom (DOFs) is used by the solver to build matrices in parallel: +This DOFSynchronizer is built based on nodes global id +Number of global dofs 32 for processor 1 +Local dof 0, global id: 0 +Local dof 1, global id: 1 +Local dof 2, global id: 6 +Local dof 3, global id: 7 +Local dof 4, global id: 12 +Local dof 5, global id: 13 +Local dof 6, global id: 24 +Local dof 7, global id: 25 +Local dof 8, global id: 16 +Local dof 9, global id: 17 +Local dof 10, global id: 4 +Local dof 11, global id: 5 +Local dof 12, global id: 14 +Local dof 13, global id: 15 +Local dof 14, global id: 20 +Local dof 15, global id: 21 +Local dof 16, global id: 2 +Local dof 17, global id: 3 +Local dof 18, global id: 18 +Local dof 19, global id: 19 +Local dof 20, global id: 10 +Local dof 21, global id: 11 +Local dof 22, global id: 22 +Local dof 23, global id: 23 +Local dof 24, global id: 28 +Local dof 25, global id: 29 +Local dof 26, global id: 26 +Local dof 27, global id: 27 + diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/output_dir_verified/output_from_proc_2_out_of_4.out b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/output_dir_verified/output_from_proc_2_out_of_4.out new file mode 100644 index 000000000..9f25b14a2 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/output_dir_verified/output_from_proc_2_out_of_4.out @@ -0,0 +1,227 @@ +Before initParallel(), non-root processors have empty Mesh object +After initParallel(), Mesh object on each processor is a local partionated mesh containing ghost elements + +Element type: _segment_2, not_ghost: 1 in the mesh of processor 2 +Element no 0 0 1 + +Element type: _triangle_3, not_ghost: 4 in the mesh of processor 2 +Element no 0 3 4 5 +Element no 1 0 5 4 +Element no 2 0 1 5 +Element no 3 1 3 5 + +Element type: _segment_2, ghost: 1 in the mesh of processor 2 +Element no 0 2 0 + +Element type: _triangle_3, ghost: 8 in the mesh of processor 2 +Element no 0 0 4 6 +Element no 1 0 6 2 +Element no 2 7 6 4 +Element no 3 0 8 9 +Element no 4 0 9 10 +Element no 5 0 2 8 +Element no 6 0 10 1 +Element no 7 1 10 11 +Nodes are also partionated and new node types are defined: +Number of nodes: 12 +Node # 0, x-coord: 0.5, y-coord: 0, of type: 0 +Node # 1, x-coord: 1, y-coord: 0, of type: -2 +Node # 2, x-coord: 0, y-coord: 0, of type: -3 +Node # 3, x-coord: 1, y-coord: 0.5, of type: -1 +Node # 4, x-coord: 0.5, y-coord: 0.5, of type: 0 +Node # 5, x-coord: 0.75, y-coord: 0.25, of type: -1 +Node # 6, x-coord: 0.25, y-coord: 0.25, of type: -3 +Node # 7, x-coord: 0, y-coord: 0.5, of type: -3 +Node # 8, x-coord: 0.25, y-coord: -0.25, of type: -3 +Node # 9, x-coord: 0.5, y-coord: -0.5, of type: -3 +Node # 10, x-coord: 0.75, y-coord: -0.25, of type: -3 +Node # 11, x-coord: 1, y-coord: -0.5, of type: -3 + +-3: pure ghost node -> not a local node +-2: master node -> node shared with other processor(s) -> local and global node +>0: slave node -> -> node shared with other processor(s) -> only local node (its id is the rank of the processor owning the master node) +Each local node has a corresponding global id used during assembly: +Global nodes ID: +0 0 +1 3 +2 6 +3 4 +4 7 +5 9 +6 10 +7 1 +8 12 +9 8 +10 11 +11 5 + +Within cohesive element model, initParallel() creates a second Mesh object usually called mesh_facet +This Mesh object contains all sub-dimensional elements where potential cohesive element can be inserted + +Element type: _point_1, not_ghost: 5 in the mesh of processor 2 +Element no 0 3 +Element no 1 4 +Element no 2 5 +Element no 3 0 +Element no 4 1 + +Element type: _segment_2, not_ghost: 7 in the mesh of processor 2 +Element no 0 3 4 +Element no 1 4 5 +Element no 2 5 3 +Element no 3 0 5 +Element no 4 0 1 +Element no 5 1 5 +Element no 6 1 3 + +Element type: _point_1, ghost: 0 in the mesh of processor 2 + +Element type: _segment_2, ghost: 10 in the mesh of processor 2 +Element no 0 0 4 +Element no 1 4 6 +Element no 2 6 0 +Element no 3 2 0 +Element no 4 4 7 +Element no 5 0 8 +Element no 6 9 0 +Element no 7 10 0 +Element no 8 10 1 +Element no 9 11 1 +The distributed synchronizer of solid mechanics model is used to synchronize fields with ghost element: +From processor 2 to processor 0 + Sending element(s): +Element [_triangle_3, 0, not_ghost] +Element [_triangle_3, 1, not_ghost] +Element [_triangle_3, 2, not_ghost] + Receiving element(s): +Element [_triangle_3, 0, ghost] +Element [_triangle_3, 1, ghost] +Element [_triangle_3, 2, ghost] +From processor 2 to processor 1 + Sending element(s): +Element [_triangle_3, 1, not_ghost] +Element [_triangle_3, 2, not_ghost] + Receiving element(s): +Element [_triangle_3, 3, ghost] +Element [_triangle_3, 5, ghost] +From processor 2 to processor 3 + Sending element(s): +Element [_triangle_3, 1, not_ghost] +Element [_triangle_3, 2, not_ghost] +Element [_triangle_3, 3, not_ghost] + Receiving element(s): +Element [_triangle_3, 4, ghost] +Element [_triangle_3, 6, ghost] +Element [_triangle_3, 7, ghost] + +In case of insertion along physical objects, cohesive elements are created during initFull() +Elements list after insertion + +Element type: _cohesive_2d_4, not_ghost: 1 in the mesh of processor 2 +Element no 0 0 1 12 13 + +Element type: _segment_2, not_ghost: 1 in the mesh of processor 2 +Element no 0 0 1 + +Element type: _triangle_3, not_ghost: 4 in the mesh of processor 2 +Element no 0 3 4 5 +Element no 1 0 5 4 +Element no 2 0 1 5 +Element no 3 1 3 5 + +Element type: _cohesive_2d_4, ghost: 1 in the mesh of processor 2 +Element no 0 2 0 2 12 + +Element type: _segment_2, ghost: 1 in the mesh of processor 2 +Element no 0 2 0 + +Element type: _triangle_3, ghost: 8 in the mesh of processor 2 +Element no 0 0 4 6 +Element no 1 0 6 2 +Element no 2 7 6 4 +Element no 3 12 8 9 +Element no 4 12 9 10 +Element no 5 12 2 8 +Element no 6 12 10 13 +Element no 7 13 10 11 +Node list after insertion: (Total number of nodes 14) +Number of nodes: 14 +Node # 0, x-coord: 0.5, y-coord: 0, of type: 0 +Node # 1, x-coord: 1, y-coord: 0, of type: -2 +Node # 2, x-coord: 0, y-coord: 0, of type: -3 +Node # 3, x-coord: 1, y-coord: 0.5, of type: -1 +Node # 4, x-coord: 0.5, y-coord: 0.5, of type: 0 +Node # 5, x-coord: 0.75, y-coord: 0.25, of type: -1 +Node # 6, x-coord: 0.25, y-coord: 0.25, of type: -3 +Node # 7, x-coord: 0, y-coord: 0.5, of type: -3 +Node # 8, x-coord: 0.25, y-coord: -0.25, of type: -3 +Node # 9, x-coord: 0.5, y-coord: -0.5, of type: -3 +Node # 10, x-coord: 0.75, y-coord: -0.25, of type: -3 +Node # 11, x-coord: 1, y-coord: -0.5, of type: -3 +Node # 12, x-coord: 0.5, y-coord: 0, of type: 0 +Node # 13, x-coord: 1, y-coord: 0, of type: -2 + +Node global ids after insertion: (Total number of nodes 16) +Global nodes ID: +0 0 +1 3 +2 6 +3 4 +4 7 +5 9 +6 10 +7 1 +8 12 +9 8 +10 11 +11 5 +12 13 +13 15 + +Solid mechanics model cohesive has its own distributed synchronizer to handle ghost cohesive element: +From processor 2 to processor 0 + Sending element(s): +Element [_cohesive_2d_4, 0, not_ghost] + Receiving element(s): +Element [_cohesive_2d_4, 0, ghost] +From processor 2 to processor 1 + Sending element(s): +Element [_cohesive_2d_4, 0, not_ghost] + Receiving element(s): +From processor 2 to processor 3 + Sending element(s): +Element [_cohesive_2d_4, 0, not_ghost] + Receiving element(s): + +A synchronizer dedicated to degrees of freedom (DOFs) is used by the solver to build matrices in parallel: +This DOFSynchronizer is built based on nodes global id +Number of global dofs 32 for processor 2 +Local dof 0, global id: 0 +Local dof 1, global id: 1 +Local dof 2, global id: 6 +Local dof 3, global id: 7 +Local dof 4, global id: 12 +Local dof 5, global id: 13 +Local dof 6, global id: 8 +Local dof 7, global id: 9 +Local dof 8, global id: 14 +Local dof 9, global id: 15 +Local dof 10, global id: 18 +Local dof 11, global id: 19 +Local dof 12, global id: 20 +Local dof 13, global id: 21 +Local dof 14, global id: 2 +Local dof 15, global id: 3 +Local dof 16, global id: 24 +Local dof 17, global id: 25 +Local dof 18, global id: 16 +Local dof 19, global id: 17 +Local dof 20, global id: 22 +Local dof 21, global id: 23 +Local dof 22, global id: 10 +Local dof 23, global id: 11 +Local dof 24, global id: 26 +Local dof 25, global id: 27 +Local dof 26, global id: 30 +Local dof 27, global id: 31 + diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/output_dir_verified/output_from_proc_3_out_of_4.out b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/output_dir_verified/output_from_proc_3_out_of_4.out new file mode 100644 index 000000000..c92fd5bb1 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/output_dir_verified/output_from_proc_3_out_of_4.out @@ -0,0 +1,225 @@ +Before initParallel(), non-root processors have empty Mesh object +After initParallel(), Mesh object on each processor is a local partionated mesh containing ghost elements + +Element type: _segment_2, not_ghost: 0 in the mesh of processor 3 + +Element type: _triangle_3, not_ghost: 4 in the mesh of processor 3 +Element no 0 3 4 5 +Element no 1 0 5 4 +Element no 2 0 4 1 +Element no 3 1 4 3 + +Element type: _segment_2, ghost: 2 in the mesh of processor 3 +Element no 0 0 1 +Element no 1 2 0 + +Element type: _triangle_3, ghost: 8 in the mesh of processor 3 +Element no 0 0 6 7 +Element no 1 0 8 6 +Element no 2 0 7 2 +Element no 3 0 1 8 +Element no 4 1 9 8 +Element no 5 0 10 5 +Element no 6 0 2 10 +Element no 7 11 5 10 +Nodes are also partionated and new node types are defined: +Number of nodes: 12 +Node # 0, x-coord: 0.5, y-coord: 0, of type: 0 +Node # 1, x-coord: 1, y-coord: 0, of type: 2 +Node # 2, x-coord: 0, y-coord: 0, of type: -3 +Node # 3, x-coord: 1, y-coord: -0.5, of type: -1 +Node # 4, x-coord: 0.75, y-coord: -0.25, of type: -1 +Node # 5, x-coord: 0.5, y-coord: -0.5, of type: 1 +Node # 6, x-coord: 0.5, y-coord: 0.5, of type: -3 +Node # 7, x-coord: 0.25, y-coord: 0.25, of type: -3 +Node # 8, x-coord: 0.75, y-coord: 0.25, of type: -3 +Node # 9, x-coord: 1, y-coord: 0.5, of type: -3 +Node # 10, x-coord: 0.25, y-coord: -0.25, of type: -3 +Node # 11, x-coord: 0, y-coord: -0.5, of type: -3 + +-3: pure ghost node -> not a local node +-2: master node -> node shared with other processor(s) -> local and global node +>0: slave node -> -> node shared with other processor(s) -> only local node (its id is the rank of the processor owning the master node) +Each local node has a corresponding global id used during assembly: +Global nodes ID: +0 0 +1 3 +2 6 +3 5 +4 11 +5 8 +6 7 +7 10 +8 9 +9 4 +10 12 +11 2 + +Within cohesive element model, initParallel() creates a second Mesh object usually called mesh_facet +This Mesh object contains all sub-dimensional elements where potential cohesive element can be inserted + +Element type: _point_1, not_ghost: 5 in the mesh of processor 3 +Element no 0 3 +Element no 1 4 +Element no 2 5 +Element no 3 0 +Element no 4 1 + +Element type: _segment_2, not_ghost: 6 in the mesh of processor 3 +Element no 0 3 4 +Element no 1 4 5 +Element no 2 5 3 +Element no 3 4 0 +Element no 4 4 1 +Element no 5 3 1 + +Element type: _point_1, ghost: 0 in the mesh of processor 3 + +Element type: _segment_2, ghost: 11 in the mesh of processor 3 +Element no 0 5 0 +Element no 1 0 1 +Element no 2 0 6 +Element no 3 7 0 +Element no 4 0 8 +Element no 5 2 0 +Element no 6 1 8 +Element no 7 1 9 +Element no 8 0 10 +Element no 9 10 5 +Element no 10 11 5 +The distributed synchronizer of solid mechanics model is used to synchronize fields with ghost element: +From processor 3 to processor 0 + Sending element(s): +Element [_triangle_3, 1, not_ghost] +Element [_triangle_3, 2, not_ghost] + Receiving element(s): +Element [_triangle_3, 0, ghost] +Element [_triangle_3, 2, ghost] +From processor 3 to processor 1 + Sending element(s): +Element [_triangle_3, 0, not_ghost] +Element [_triangle_3, 1, not_ghost] +Element [_triangle_3, 2, not_ghost] + Receiving element(s): +Element [_triangle_3, 5, ghost] +Element [_triangle_3, 6, ghost] +Element [_triangle_3, 7, ghost] +From processor 3 to processor 2 + Sending element(s): +Element [_triangle_3, 1, not_ghost] +Element [_triangle_3, 2, not_ghost] +Element [_triangle_3, 3, not_ghost] + Receiving element(s): +Element [_triangle_3, 1, ghost] +Element [_triangle_3, 3, ghost] +Element [_triangle_3, 4, ghost] + +In case of insertion along physical objects, cohesive elements are created during initFull() +Elements list after insertion + +Element type: _cohesive_2d_4, not_ghost: 0 in the mesh of processor 3 + +Element type: _segment_2, not_ghost: 0 in the mesh of processor 3 + +Element type: _triangle_3, not_ghost: 4 in the mesh of processor 3 +Element no 0 3 4 5 +Element no 1 13 5 4 +Element no 2 13 4 12 +Element no 3 12 4 3 + +Element type: _cohesive_2d_4, ghost: 2 in the mesh of processor 3 +Element no 0 0 1 13 12 +Element no 1 2 0 2 13 + +Element type: _segment_2, ghost: 2 in the mesh of processor 3 +Element no 0 0 1 +Element no 1 2 0 + +Element type: _triangle_3, ghost: 8 in the mesh of processor 3 +Element no 0 0 6 7 +Element no 1 0 8 6 +Element no 2 0 7 2 +Element no 3 0 1 8 +Element no 4 1 9 8 +Element no 5 13 10 5 +Element no 6 13 2 10 +Element no 7 11 5 10 +Node list after insertion: (Total number of nodes 14) +Number of nodes: 14 +Node # 0, x-coord: 0.5, y-coord: 0, of type: 0 +Node # 1, x-coord: 1, y-coord: 0, of type: 2 +Node # 2, x-coord: 0, y-coord: 0, of type: -3 +Node # 3, x-coord: 1, y-coord: -0.5, of type: -1 +Node # 4, x-coord: 0.75, y-coord: -0.25, of type: -1 +Node # 5, x-coord: 0.5, y-coord: -0.5, of type: 1 +Node # 6, x-coord: 0.5, y-coord: 0.5, of type: -3 +Node # 7, x-coord: 0.25, y-coord: 0.25, of type: -3 +Node # 8, x-coord: 0.75, y-coord: 0.25, of type: -3 +Node # 9, x-coord: 1, y-coord: 0.5, of type: -3 +Node # 10, x-coord: 0.25, y-coord: -0.25, of type: -3 +Node # 11, x-coord: 0, y-coord: -0.5, of type: -3 +Node # 12, x-coord: 1, y-coord: 0, of type: 2 +Node # 13, x-coord: 0.5, y-coord: 0, of type: 0 + +Node global ids after insertion: (Total number of nodes 16) +Global nodes ID: +0 0 +1 3 +2 6 +3 5 +4 11 +5 8 +6 7 +7 10 +8 9 +9 4 +10 12 +11 2 +12 15 +13 13 + +Solid mechanics model cohesive has its own distributed synchronizer to handle ghost cohesive element: +From processor 3 to processor 0 + Sending element(s): + Receiving element(s): +Element [_cohesive_2d_4, 1, ghost] +From processor 3 to processor 1 + Sending element(s): + Receiving element(s): +From processor 3 to processor 2 + Sending element(s): + Receiving element(s): +Element [_cohesive_2d_4, 0, ghost] + +A synchronizer dedicated to degrees of freedom (DOFs) is used by the solver to build matrices in parallel: +This DOFSynchronizer is built based on nodes global id +Number of global dofs 32 for processor 3 +Local dof 0, global id: 0 +Local dof 1, global id: 1 +Local dof 2, global id: 6 +Local dof 3, global id: 7 +Local dof 4, global id: 12 +Local dof 5, global id: 13 +Local dof 6, global id: 10 +Local dof 7, global id: 11 +Local dof 8, global id: 22 +Local dof 9, global id: 23 +Local dof 10, global id: 16 +Local dof 11, global id: 17 +Local dof 12, global id: 14 +Local dof 13, global id: 15 +Local dof 14, global id: 20 +Local dof 15, global id: 21 +Local dof 16, global id: 18 +Local dof 17, global id: 19 +Local dof 18, global id: 8 +Local dof 19, global id: 9 +Local dof 20, global id: 24 +Local dof 21, global id: 25 +Local dof 22, global id: 4 +Local dof 23, global id: 5 +Local dof 24, global id: 30 +Local dof 25, global id: 31 +Local dof 26, global id: 26 +Local dof 27, global id: 27 + diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/test_cohesive_parallel_insertion_along_physical_surfaces.cc b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/test_cohesive_parallel_insertion_along_physical_surfaces.cc new file mode 100644 index 000000000..74d52adbf --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/test_cohesive_parallel_insertion_along_physical_surfaces.cc @@ -0,0 +1,118 @@ +/** + * @file test_cohesive_insertion_along_physical_surfaces.cc + * @author Fabian Barras + * @date Fri Aug 7 09:07:44 2015 + * + * @brief Test parallel intrinsic insertion of cohesive elements along physical surfaces + * + * @section LICENSE + * + * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) any + * later version. + * + * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR + * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with Akantu. If not, see . + * + */ + +/* -------------------------------------------------------------------------- */ +#include +#include +#include +/* -------------------------------------------------------------------------- */ +#include "aka_common.hh" +#include "mesh.hh" +#include "mesh_io.hh" +#include "mesh_io_msh.hh" +#include "mesh_utils.hh" +#include "solid_mechanics_model_cohesive.hh" +#include "material.hh" +#include "material_cohesive.hh" + +/* -------------------------------------------------------------------------- */ +using namespace akantu; + +int main(int argc, char *argv[]) { + + initialize("input_file.dat", argc, argv); + + Math::setTolerance(1e-15); + + const UInt spatial_dimension = 3; + + Mesh mesh(spatial_dimension); + + StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + Int psize = comm.getNbProc(); + Int prank = comm.whoAmI(); + akantu::MeshPartition * partition = NULL; + + if(prank==0){ + + mesh.read("3d_spherical_inclusion.msh"); + + partition = new MeshPartitionScotch(mesh, spatial_dimension); + partition->partitionate(psize); + + } + SolidMechanicsModelCohesive model(mesh); + model.initParallel(partition); + mesh.createGroupsFromMeshData("physical_names"); + model.initFull(SolidMechanicsModelCohesiveOptions(_static)); + + std::vector surfaces_name = {"interface", "coh1", "coh2", "coh3", "coh4", "coh5"}; + UInt nb_surf = surfaces_name.size(); + + for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { + + std::string ghost_str; + + if(*gt == 1) ghost_str = "ghost"; + else ghost_str = "not ghost"; + + Mesh::type_iterator it = mesh.firstType(spatial_dimension, *gt, _ek_cohesive); + Mesh::type_iterator end = mesh.lastType(spatial_dimension, *gt, _ek_cohesive); + + for(; it != end; ++it) { + + Array & material_id = mesh.getMeshFacets().getData("physical_names")(mesh.getFacetType(*it), *gt); + + for (UInt i = 0; i < nb_surf; ++i) { + + UInt expected_insertion = 0; + + for(UInt m = 0; m " << inserted_elements << " inserted elements of type " << ghost_str + << " out of " + << expected_insertion << std::endl; + return EXIT_FAILURE; + } + } + } + } + + model.assembleStiffnessMatrix(); + + finalize(); + + return EXIT_SUCCESS; +} diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/test_cohesive_parallel_insertion_along_physical_surfaces.sh b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/test_cohesive_parallel_insertion_along_physical_surfaces.sh new file mode 100755 index 000000000..47c6b9eb4 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/test_cohesive_parallel_insertion_along_physical_surfaces.sh @@ -0,0 +1,3 @@ +#!/bin/bash + +mpirun -np 8 test_cohesive_parallel_insertion_along_physical_surfaces diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/test_cohesive_parallel_intrinsic_implicit_insertion.cc b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/test_cohesive_parallel_intrinsic_implicit_insertion.cc new file mode 100644 index 000000000..d6643a5a1 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/test_cohesive_parallel_intrinsic_implicit_insertion.cc @@ -0,0 +1,270 @@ +/** + * @file test_cohesive_parallel_intrinsic_implicit_insertion.cc + * @author Fabian Barras + * @date Fri Aug 5 17:05:59 2016 + * + * @brief Verifying the proper insertion and synchronization of intrinsic cohesive elements + * + * @section LICENSE + * + * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) any + * later version. + * + * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR + * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with Akantu. If not, see . + * + */ + +/* -------------------------------------------------------------------------- */ +#include +#include +#include +/* -------------------------------------------------------------------------- */ +#include "aka_common.hh" +#include "mesh.hh" +#include "mesh_io.hh" +#include "mesh_io_msh.hh" +#include "mesh_utils.hh" +#include "solid_mechanics_model_cohesive.hh" +#include "material.hh" +#include "material_cohesive.hh" +/* -------------------------------------------------------------------------- */ +using namespace akantu; +std::ofstream output; + +/* -------------------------------------------------------------------------- */ +void printMeshContent(Mesh & mesh) { + + StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + Int prank = comm.whoAmI(); + comm.barrier(); + for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { + + Mesh::type_iterator first = mesh.firstType(_all_dimensions, *gt, _ek_not_defined); + Mesh::type_iterator last = mesh.lastType(_all_dimensions, *gt, _ek_not_defined); + + for(; first != last; ++first) { + UInt nb_element = mesh.getNbElement(*first,*gt); + output << std::endl + << "Element type: " << *first << ", " << *gt << ": " + << nb_element << " in the mesh of processor " << prank << std::endl; + Array & conn = mesh.getConnectivity(*first,*gt); + for (UInt i = 0; i < conn.getSize(); ++i) { + output << "Element no " << i << " "; + for (UInt j = 0; j < conn.getNbComponent(); ++j) { + output << conn(i,j) << " "; + } + output << std::endl; + } + } + } +} + +/* -------------------------------------------------------------------------- */ +void printNodeList(Mesh & mesh) { + + Array & nodes = mesh.getNodes(); + + output << "Number of nodes: " << mesh.getNbNodes()<< std::endl; + for (UInt i = 0; i < mesh.getNbNodes(); ++i) { + output<<"Node # " << i + << ", x-coord: " + << nodes(i,0) + << ", y-coord: " + << nodes(i,1) + << ", of type: " << mesh.getNodeType(i) << std::endl; + } + output << std::endl; +} + +/* -------------------------------------------------------------------------- */ +void getGlobalIDs(Mesh & mesh) { + + const Array & glob_id = mesh.getGlobalNodesIds(); + if(&glob_id) { + output << "Global nodes ID: " << std::endl; + for (UInt i = 0; i < glob_id.getSize(); ++i) { + output << i << " " << glob_id(i) << std::endl; + } + } + output << std::endl; +} + +/* -------------------------------------------------------------------------- */ +void printSynchroinfo(Mesh & mesh, const DistributedSynchronizer & synch) { + + StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + Int prank = comm.whoAmI(); + Int psize = comm.getNbProc(); + + if(comm.getNbProc()==1) + return; + + for (Int p = 0; p < psize; ++p) { + + if (p==prank) + continue; + output << "From processor " << prank << " to processor " << p << std::endl; + const Array & sele = *(synch.getSendElement()+p); + output << " Sending element(s): " << std::endl; + for (UInt i = 0; i < sele.getSize(); ++i) { + output << sele(i) << std::endl; + } + + const Array & rele = *(synch.getReceiveElement()+p); + output <<" Receiving element(s): " << std::endl; + for (UInt i = 0; i < rele.getSize(); ++i) { + output << rele(i) << std::endl; + } + } + output << std::endl; +} + +/* -------------------------------------------------------------------------- */ +void printDOF(SolidMechanicsModelCohesive & model) { + + StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + if(comm.getNbProc()==1) + return; + Int prank = comm.whoAmI(); + + const DOFSynchronizer & dof = model.getDOFSynchronizer(); + + output << "Number of global dofs " << dof.getNbGlobalDOFs() << " for processor "<< prank << std::endl; + + const Array & dof_global_ids = dof.getDOFGlobalIDs(); + + for (UInt i = 0; i < dof_global_ids.getSize(); ++i) { + + output << "Local dof " << i << ", global id: " << dof_global_ids(i) << std::endl; + } + output << std::endl; +} + +/* -------------------------------------------------------------------------- */ +int main(int argc, char *argv[]) { + + std::string input_file = "input_file_iii.dat"; + std::string mesh_file = "2d_basic_interface.msh"; + std::string dir = "output_dir/"; + + initialize(input_file, argc, argv); + + debug::setDebugLevel(dbl0); + + const UInt spatial_dimension = 2; + + Mesh mesh(spatial_dimension); + + StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + Int psize = comm.getNbProc(); + Int prank = comm.whoAmI(); + akantu::MeshPartition * partition = NULL; + + std::stringstream filename; + filename << dir.c_str() << "output_from_proc_" << prank << "_out_of_" << psize << ".out"; + output.open(filename.str()); + + if(prank==0){ + + mesh.read(mesh_file); + partition = new MeshPartitionScotch(mesh, spatial_dimension); + partition->partitionate(psize); + + const ElementTypeMapArray & partitions = partition->getPartitions(); + + output << "The root processor read the mesh." << std::endl + << "Only GMSH physical objects are created in the mesh." < part = partitions(*first, *gt); + + for (UInt i = 0; i < part.getSize(); ++i) { + output << i << " " << part(i) << std::endl; + } + } + } + + output << "Nodes are also read and set with type -1 (normal node)" << std::endl; + printNodeList(mesh); + } + + SolidMechanicsModelCohesive model(mesh); + + output << "Before initParallel(), non-root processors have empty Mesh object" << std::endl; + printMeshContent(mesh); + + model.initParallel(partition); + + output << "After initParallel(), Mesh object on each processor is a local partionated mesh containing ghost elements" << std::endl; + printMeshContent(mesh); + + output << "Nodes are also partionated and new node types are defined:" << std::endl; + printNodeList(mesh); + output << "-3: pure ghost node -> not a local node" << std::endl + << "-2: master node -> node shared with other processor(s) -> local and global node" <0: slave node -> -> node shared with other processor(s) -> only local node (its id is the rank of the processor owning the master node)" <("physical_names"); + model.initFull(SolidMechanicsModelCohesiveOptions(_static)); + + output << "In case of insertion along physical objects, cohesive elements are created during initFull()" << std::endl; + output << "Elements list after insertion" << std::endl; + printMeshContent(mesh); + + output << "Node list after insertion: (Total number of nodes " << mesh.getNbNodes()<< ")"<< std::endl; + printNodeList(mesh); + + output << "Node global ids after insertion: (Total number of nodes " << mesh.getNbGlobalNodes()<< ")"<< std::endl; + getGlobalIDs(mesh); + + const DistributedSynchronizer & coh_synch_model = *(model.getCohesiveSynchronizer()); + output << "Solid mechanics model cohesive has its own distributed synchronizer to handle ghost cohesive element:" << std::endl; + printSynchroinfo(mesh, coh_synch_model); + + output << "A synchronizer dedicated to degrees of freedom (DOFs) is used by the solver to build matrices in parallel:" + << std::endl + << "This DOFSynchronizer is built based on nodes global id " << std::endl; + printDOF(model); + + output.close(); + + finalize(); + + return EXIT_SUCCESS; +} diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/verify_insertion.sh b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/verify_insertion.sh new file mode 100644 index 000000000..7968e0897 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_insertion/verify_insertion.sh @@ -0,0 +1,18 @@ +#!/bin/bash + +status=0 + +PARALLEL_LEVEL=( 2 4 ) +pwd +for (( p=0; p<$PARALLEL_LEVEL; ++p)); do + PRANK=${PARALLEL_LEVEL[${p}]} + + for (( i=0; i<$PRANK; ++i)); do + echo "Checking output_from_proc_"${i}"_out_of_"${PRANK}".out" + FILE_NAME="output_from_proc_"${i}"_out_of_"${PRANK}".out" + diff "output_dir/"${FILE_NAME} "output_dir_verified/"${FILE_NAME} + status=$(($? + ${status})) + done +done + +exit ${status} diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/CMakeLists.txt b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/CMakeLists.txt new file mode 100644 index 000000000..a7dbaf598 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/CMakeLists.txt @@ -0,0 +1,33 @@ +#=============================================================================== +# @file CMakeLists.txt +# +# @author Marco Vocialta +# +# +# @brief configuration for parallel test for intrinsic cohesive elements +# +# @section LICENSE +# +# Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) +# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) +# +# @section DESCRIPTION +# +#=============================================================================== + +add_mesh(test_cohesive_parallel_intrinsic_mesh mesh.geo 2 2) +add_mesh(test_cohesive_parallel_intrinsic_tetrahedron_mesh tetrahedron.geo 3 2) + +register_test(test_cohesive_parallel_intrinsic + SOURCES test_cohesive_parallel_intrinsic.cc + DEPENDS test_cohesive_parallel_intrinsic_mesh + PACKAGE parallel_cohesive_element + FILES_TO_COPY material.dat + DIRECTORIES_TO_CREATE paraview) + +register_test(test_cohesive_parallel_intrinsic_tetrahedron + SOURCES test_cohesive_parallel_intrinsic_tetrahedron.cc + DEPENDS test_cohesive_parallel_intrinsic_tetrahedron_mesh + PACKAGE parallel_cohesive_element + FILES_TO_COPY material_tetrahedron.dat + DIRECTORIES_TO_CREATE paraview) diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/material.dat b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/material.dat new file mode 100644 index 000000000..eecf76c42 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/material.dat @@ -0,0 +1,14 @@ +material elastic [ + name = bulk + rho = 1e3 # density + E = 1e3 # young's modulus + nu = 0.1 # poisson's ratio +] + +material cohesive_bilinear [ + name = cohesive + sigma_c = 1 + beta = 1.5 + G_c = 1 + delta_0 = 0.1 +] diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/material_tetrahedron.dat b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/material_tetrahedron.dat new file mode 100644 index 000000000..ba242d682 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/material_tetrahedron.dat @@ -0,0 +1,16 @@ +material elastic [ + name = steel + rho = 10 # density + E = 10 # young's modulus + nu = 0.1 # poisson's ratio +] + +material cohesive_bilinear [ + name = cohesive + sigma_c = 1 + beta = 1.5 + G_c = 1 + kappa = 1.2 + delta_0 = 0.1 +# penalty = 1e10 +] diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/mesh.geo b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/mesh.geo new file mode 100644 index 000000000..cbe3dd7a7 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/mesh.geo @@ -0,0 +1,15 @@ +h = 1.; + +Point(1) = { 1, 1, 0, h}; +Point(2) = {-1, 1, 0, h}; +Point(3) = {-1,-1, 0, h}; +Point(4) = { 1,-1, 0, h}; +Line(1) = {1, 2}; +Line(2) = {2, 3}; +Line(3) = {3, 4}; +Line(4) = {4, 1}; +Line Loop(5) = {1, 2, 3, 4}; +Plane Surface(6) = {5}; +Physical Surface(7) = {6}; + +Transfinite Line {1, 2, 3, 4} = 3; \ No newline at end of file diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/test_cohesive_parallel_intrinsic.cc b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/test_cohesive_parallel_intrinsic.cc new file mode 100644 index 000000000..9bd19c8b7 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/test_cohesive_parallel_intrinsic.cc @@ -0,0 +1,160 @@ +/** + * @file test_cohesive_parallel_intrinsic.cc + * + * @author Marco Vocialta + * + * + * @brief parallel test for intrinsic cohesive elements + * + * @section LICENSE + * + * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + */ + +/* -------------------------------------------------------------------------- */ +#include "solid_mechanics_model_cohesive.hh" + +/* -------------------------------------------------------------------------- */ +using namespace akantu; + +int main(int argc, char *argv[]) { + initialize("material.dat", argc, argv); + + const UInt max_steps = 350; + + UInt spatial_dimension = 2; + Mesh mesh(spatial_dimension); + + StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + Int psize = comm.getNbProc(); + Int prank = comm.whoAmI(); + + akantu::MeshPartition * partition = NULL; + if(prank == 0) { + // Read the mesh + mesh.read("mesh.msh"); + + // /// insert cohesive elements + // CohesiveElementInserter inserter(mesh); + // inserter.setLimit('x', -0.26, -0.24); + // inserter.insertIntrinsicElements(); + + /// partition the mesh + partition = new MeshPartitionScotch(mesh, spatial_dimension); + // debug::setDebugLevel(dblDump); + partition->partitionate(psize); + // debug::setDebugLevel(dblWarning); + } + + SolidMechanicsModelCohesive model(mesh); + + model.initParallel(partition); + + model.initFull(); + + model.limitInsertion(_x, -0.26, -0.24); + model.insertIntrinsicElements(); + + debug::setDebugLevel(dblDump); + std::cout << mesh << std::endl; + debug::setDebugLevel(dblWarning); + + Real time_step = model.getStableTimeStep()*0.8; + model.setTimeStep(time_step); + // std::cout << "Time step: " << time_step << std::endl; + + model.assembleMassLumped(); + + + Array & position = mesh.getNodes(); + Array & velocity = model.getVelocity(); + Array & boundary = model.getBlockedDOFs(); + // Array & displacement = model.getDisplacement(); + // const Array & residual = model.getResidual(); + + UInt nb_nodes = mesh.getNbNodes(); + Real epsilon = std::numeric_limits::epsilon(); + + for (UInt n = 0; n < nb_nodes; ++n) { + if (std::abs(position(n, 0) - 1.) < epsilon) + boundary(n, 0) = true; + } + + model.synchronizeBoundaries(); + model.updateResidual(); + + model.setBaseName("intrinsic_parallel"); + model.addDumpFieldVector("displacement"); + model.addDumpField("velocity" ); + model.addDumpField("acceleration"); + model.addDumpField("residual" ); + model.addDumpField("stress"); + model.addDumpField("strain"); + model.addDumpField("partitions"); + model.addDumpField("force"); + model.dump(); + + + model.setBaseNameToDumper("cohesive elements", + "cohesive_elements_parallel_intrinsic"); + model.addDumpFieldVectorToDumper("cohesive elements", "displacement"); + model.dump("cohesive elements"); + + /// initial conditions + Real loading_rate = .2; + for (UInt n = 0; n < nb_nodes; ++n) { + velocity(n, 0) = loading_rate * position(n, 0); + } + + /// Main loop + for (UInt s = 1; s <= max_steps; ++s) { + + model.solveStep(); + + if(s % 20 == 0) { + model.dump(); + model.dump("cohesive elements"); + if(prank == 0) std::cout << "passing step " << s << "/" << max_steps << std::endl; + } + + // // update displacement + // for (UInt n = 0; n < nb_nodes; ++n) { + // if (position(n, 1) + displacement(n, 1) > 0) { + // displacement(n, 0) -= 0.01; + // } + // } + + // Real Ed = dynamic_cast (model.getMaterial(1)).getDissipatedEnergy(); + // Real Er = dynamic_cast (model.getMaterial(1)).getReversibleEnergy(); + + // edis << s << " " + // << Ed << std::endl; + + // erev << s << " " + // << Er << std::endl; + + } + + // edis.close(); + // erev.close(); + + Real Ed = model.getEnergy("dissipated"); + + Real Edt = 2 * sqrt(2); + + + if(prank == 0) { + std::cout << Ed << " " << Edt << std::endl; + + if (std::abs((Ed - Edt) / Edt) > 0.01 || std::isnan(Ed)) { + std::cout << "The dissipated energy is incorrect" << std::endl; + return EXIT_FAILURE; + } + } + + finalize(); + if(prank == 0) std::cout << "OK: Test passed!" << std::endl; + return EXIT_SUCCESS; +} diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/test_cohesive_parallel_intrinsic.sh b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/test_cohesive_parallel_intrinsic.sh new file mode 100755 index 000000000..d24f982d0 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/test_cohesive_parallel_intrinsic.sh @@ -0,0 +1,3 @@ +#!/bin/bash + +mpirun -np 3 ./test_cohesive_parallel_intrinsic diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/test_cohesive_parallel_intrinsic_tetrahedron.cc b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/test_cohesive_parallel_intrinsic_tetrahedron.cc new file mode 100644 index 000000000..f6607255d --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/test_cohesive_parallel_intrinsic_tetrahedron.cc @@ -0,0 +1,704 @@ +/** + * @file test_cohesive_parallel_intrinsic_tetrahedron.cc + * + * @author Marco Vocialta + * + * + * @brief Test for 3D intrinsic cohesive elements simulation in parallel + * + * @section LICENSE + * + * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + */ + +/* -------------------------------------------------------------------------- */ +#include "solid_mechanics_model_cohesive.hh" +#include "dumper_paraview.hh" +#include "material_cohesive.hh" + +/* -------------------------------------------------------------------------- */ +using namespace akantu; + +void updateDisplacement(SolidMechanicsModelCohesive & model, + const ElementTypeMapArray & elements, + Vector & increment); + +bool checkTractions(SolidMechanicsModelCohesive & model, + Vector & opening, + Vector & theoretical_traction, + Matrix & rotation); + +void findNodesToCheck(const Mesh & mesh, + const ElementTypeMapArray & elements, + Array & nodes_to_check, + Int psize); + +bool checkEquilibrium(const Mesh & mesh, + const Array & residual); + +bool checkResidual(const Array & residual, + const Vector & traction, + const Array & nodes_to_check, + const Matrix & rotation); + +void findElementsToDisplace(const Mesh & mesh, + ElementTypeMapArray & elements); + +int main(int argc, char *argv[]) { + initialize("material_tetrahedron.dat", argc, argv); + + const UInt spatial_dimension = 3; + const UInt max_steps = 60; + const Real increment_constant = 0.01; + ElementType type = _tetrahedron_10; + Math::setTolerance(1.e-10); + + Mesh mesh(spatial_dimension); + + StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + Int psize = comm.getNbProc(); + Int prank = comm.whoAmI(); + + UInt nb_nodes_to_check_serial = 0; + UInt total_nb_nodes = 0; + UInt nb_elements_check_serial = 0; + + akantu::MeshPartition * partition = NULL; + if(prank == 0) { + // Read the mesh + mesh.read("tetrahedron.msh"); + + /// count nodes with zero position + const Array & position = mesh.getNodes(); + for (UInt n = 0; n < position.getSize(); ++n) { + if (std::abs(position(n, 0) - 0.) < 1e-6) + ++nb_nodes_to_check_serial; + } + + // /// insert cohesive elements + // CohesiveElementInserter inserter(mesh); + // inserter.setLimit(0, -0.01, 0.01); + // inserter.insertIntrinsicElements(); + + /// find nodes to check in serial + ElementTypeMapArray elements_serial("elements_serial", ""); + findElementsToDisplace(mesh, elements_serial); + nb_elements_check_serial = elements_serial(type).getSize(); + + total_nb_nodes = mesh.getNbNodes() + nb_nodes_to_check_serial; + + /// partition the mesh + partition = new MeshPartitionScotch(mesh, spatial_dimension); + debug::setDebugLevel(dblDump); + partition->partitionate(psize); + debug::setDebugLevel(dblInfo); + } + + comm.broadcast(&nb_nodes_to_check_serial, 1, 0); + comm.broadcast(&nb_elements_check_serial, 1, 0); + + SolidMechanicsModelCohesive model(mesh); + + model.initParallel(partition); + model.initFull(); + + model.limitInsertion(_x, -0.01, 0.01); + model.insertIntrinsicElements(); + + { + comm.broadcast(&total_nb_nodes, 1, 0); + + Array nb_local_nodes(psize); + nb_local_nodes.clear(); + + for (UInt n = 0; n < mesh.getNbNodes(); ++n) { + if (mesh.isLocalOrMasterNode(n)) ++nb_local_nodes(prank); + } + + comm.allGather(nb_local_nodes.storage(), 1); + + UInt total_nb_nodes_parallel = std::accumulate(nb_local_nodes.begin(), + nb_local_nodes.end(), 0); + + Array global_nodes_list(total_nb_nodes_parallel); + + UInt first_global_node = std::accumulate(nb_local_nodes.begin(), + nb_local_nodes.begin() + prank, 0); + + for (UInt n = 0; n < mesh.getNbNodes(); ++n) { + if (mesh.isLocalOrMasterNode(n)) { + global_nodes_list(first_global_node) = mesh.getNodeGlobalId(n); + ++first_global_node; + } + } + + comm.allGatherV(global_nodes_list.storage(), nb_local_nodes.storage()); + + if (prank == 0) + std::cout << "Maximum node index: " + << *(std::max_element(global_nodes_list.begin(), + global_nodes_list.end())) << std::endl; + + Array repeated_nodes; + repeated_nodes.resize(0); + + for (UInt n = 0; n < total_nb_nodes_parallel; ++n) { + UInt appearances = std::count(global_nodes_list.begin() + n, + global_nodes_list.end(), + global_nodes_list(n)); + + if (appearances > 1) { + std::cout << "Node " << global_nodes_list(n) + << " appears " << appearances << " times" << std::endl; + + std::cout << " in position: " << n; + + repeated_nodes.push_back(global_nodes_list(n)); + + UInt * node_position = global_nodes_list.storage() + n; + + for (UInt i = 1; i < appearances; ++i) { + node_position = std::find(node_position + 1, + global_nodes_list.storage() + + total_nb_nodes_parallel, + global_nodes_list(n)); + + UInt current_index = node_position - global_nodes_list.storage(); + + std::cout << ", " << current_index; + } + std::cout << std::endl << std::endl; + } + } + + + for (UInt n = 0; n < mesh.getNbNodes(); ++n) { + UInt global_node = mesh.getNodeGlobalId(n); + + if (std::find(repeated_nodes.begin(), repeated_nodes.end(), global_node) + != repeated_nodes.end()) { + std::cout << "Repeated global node " << global_node + << " corresponds to local node " << n << std::endl; + } + } + + + if (total_nb_nodes != total_nb_nodes_parallel) { + if (prank == 0) { + std::cout << "Error: total number of nodes is wrong in parallel" << std::endl; + std::cout << "Serial: " << total_nb_nodes + << " Parallel: " << total_nb_nodes_parallel << std::endl; + } + finalize(); + return EXIT_FAILURE; + } + } + + model.updateResidual(); + + model.setBaseName("intrinsic_parallel_tetrahedron"); + model.addDumpFieldVector("displacement"); + model.addDumpField("residual"); + model.addDumpField("partitions"); + model.dump(); + + model.setBaseNameToDumper("cohesive elements", + "cohesive_elements_parallel_tetrahedron"); + model.addDumpFieldVectorToDumper("cohesive elements", "displacement"); + model.dump("cohesive elements"); + + /// find elements to displace + ElementTypeMapArray elements("elements", ""); + findElementsToDisplace(mesh, elements); + + UInt nb_elements_check = elements(type).getSize(); + comm.allReduce(&nb_elements_check, 1, _so_sum); + + if (nb_elements_check != nb_elements_check_serial) { + if (prank == 0) { + std::cout << "Error: number of elements to check is wrong" << std::endl; + std::cout << "Serial: " << nb_elements_check_serial + << " Parallel: " << nb_elements_check << std::endl; + } + finalize(); + return EXIT_FAILURE; + } + + /// find nodes to check + Array nodes_to_check; + findNodesToCheck(mesh, elements, nodes_to_check, psize); + + Vector nodes_to_check_size(psize); + nodes_to_check_size(prank) = nodes_to_check.getSize(); + comm.allGather(nodes_to_check_size.storage(), 1); + + UInt nodes_to_check_global_size = std::accumulate(nodes_to_check_size.storage(), + nodes_to_check_size.storage() + + psize, 0); + + if (nodes_to_check_global_size != nb_nodes_to_check_serial) { + if (prank == 0) { + std::cout << "Error: number of nodes to check is wrong in parallel" << std::endl; + std::cout << "Serial: " << nb_nodes_to_check_serial + << " Parallel: " << nodes_to_check_global_size << std::endl; + } + finalize(); + return EXIT_FAILURE; + } + + /// rotate mesh + Real angle = 1.; + + Matrix rotation(spatial_dimension, spatial_dimension); + rotation.clear(); + rotation(0, 0) = std::cos(angle); + rotation(0, 1) = std::sin(angle) * -1.; + rotation(1, 0) = std::sin(angle); + rotation(1, 1) = std::cos(angle); + rotation(2, 2) = 1.; + + + Vector increment_tmp(spatial_dimension); + for (UInt dim = 0; dim < spatial_dimension; ++dim) { + increment_tmp(dim) = (dim + 1) * increment_constant; + } + + Vector increment(spatial_dimension); + increment.mul(rotation, increment_tmp); + + Array & position = mesh.getNodes(); + Array position_tmp(position); + + Array::iterator > position_it = position.begin(spatial_dimension); + Array::iterator > position_end = position.end(spatial_dimension); + Array::iterator > position_tmp_it + = position_tmp.begin(spatial_dimension); + + for (; position_it != position_end; ++position_it, ++position_tmp_it) + position_it->mul(rotation, *position_tmp_it); + + model.dump(); + model.dump("cohesive elements"); + + updateDisplacement(model, elements, increment); + + + Real theoretical_Ed = 0; + + Vector opening(spatial_dimension); + Vector traction(spatial_dimension); + Vector opening_old(spatial_dimension); + Vector traction_old(spatial_dimension); + + opening.clear(); + traction.clear(); + opening_old.clear(); + traction_old.clear(); + + Vector Dt(spatial_dimension); + Vector Do(spatial_dimension); + + const Array & residual = model.getResidual(); + + /// Main loop + for (UInt s = 1; s <= max_steps; ++s) { + + model.updateResidual(); + + opening += increment_tmp; + if (checkTractions(model, opening, traction, rotation) || + checkEquilibrium(mesh, residual) || + checkResidual(residual, traction, nodes_to_check, rotation)) { + finalize(); + return EXIT_FAILURE; + } + + /// compute energy + Do = opening; + Do -= opening_old; + + Dt = traction_old; + Dt += traction; + + theoretical_Ed += .5 * Do.dot(Dt); + + opening_old = opening; + traction_old = traction; + + + updateDisplacement(model, elements, increment); + + if(s % 10 == 0) { + if (prank == 0) + std::cout << "passing step " << s << "/" << max_steps << std::endl; + model.dump(); + model.dump("cohesive elements"); + } + } + + model.dump(); + model.dump("cohesive elements"); + + Real Ed = model.getEnergy("dissipated"); + + theoretical_Ed *= 4.; + + if (prank == 0) + std::cout << "Dissipated energy: " << Ed + << ", theoretical value: " << theoretical_Ed << std::endl; + + if (!Math::are_float_equal(Ed, theoretical_Ed) || std::isnan(Ed)) { + + if (prank == 0) + std::cout << "Error: the dissipated energy is incorrect" << std::endl; + + finalize(); + return EXIT_FAILURE; + } + + finalize(); + if(prank == 0) std::cout << "OK: Test passed!" << std::endl; + return EXIT_SUCCESS; +} + +/* -------------------------------------------------------------------------- */ + +void updateDisplacement(SolidMechanicsModelCohesive & model, + const ElementTypeMapArray & elements, + Vector & increment) { + + UInt spatial_dimension = model.getSpatialDimension(); + Mesh & mesh = model.getFEEngine().getMesh(); + UInt nb_nodes = mesh.getNbNodes(); + + Array & displacement = model.getDisplacement(); + Array update(nb_nodes); + update.clear(); + + for (ghost_type_t::iterator gt = ghost_type_t::begin(); + gt != ghost_type_t::end(); + ++gt) { + + GhostType ghost_type = *gt; + Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); + Mesh::type_iterator last = mesh.lastType(spatial_dimension, ghost_type); + + for (; it != last; ++it) { + ElementType type = *it; + + const Array & elem = elements(type, ghost_type); + const Array & connectivity = mesh.getConnectivity(type, ghost_type); + UInt nb_nodes_per_element = connectivity.getNbComponent(); + + for (UInt el = 0; el < elem.getSize(); ++el) { + for (UInt n = 0; n < nb_nodes_per_element; ++n) { + UInt node = connectivity(elem(el), n); + if (!update(node)) { + Vector node_disp(displacement.storage() + node * spatial_dimension, + spatial_dimension); + node_disp += increment; + update(node) = true; + } + } + } + } + } +} + +/* -------------------------------------------------------------------------- */ + +bool checkTractions(SolidMechanicsModelCohesive & model, + Vector & opening, + Vector & theoretical_traction, + Matrix & rotation) { + UInt spatial_dimension = model.getSpatialDimension(); + const Mesh & mesh = model.getMesh(); + + const MaterialCohesive & mat_cohesive + = dynamic_cast < const MaterialCohesive & > (model.getMaterial(1)); + + Real sigma_c = mat_cohesive.getParam< RandomInternalField >("sigma_c"); + const Real beta = mat_cohesive.getParam("beta"); + const Real G_cI = mat_cohesive.getParam("G_c"); + // Real G_cII = mat_cohesive.getParam("G_cII"); + const Real delta_0 = mat_cohesive.getParam("delta_0"); + const Real kappa = mat_cohesive.getParam("kappa"); + Real delta_c = 2 * G_cI / sigma_c; + sigma_c *= delta_c / (delta_c - delta_0); + + Vector normal_opening(spatial_dimension); + normal_opening.clear(); + normal_opening(0) = opening(0); + Real normal_opening_norm = normal_opening.norm(); + + Vector tangential_opening(spatial_dimension); + tangential_opening.clear(); + for (UInt dim = 1; dim < spatial_dimension; ++dim) + tangential_opening(dim) = opening(dim); + + Real tangential_opening_norm = tangential_opening.norm(); + + Real beta2_kappa2 = beta * beta/kappa/kappa; + Real beta2_kappa = beta * beta/kappa; + + Real delta = std::sqrt(tangential_opening_norm * tangential_opening_norm + * beta2_kappa2 + + normal_opening_norm * normal_opening_norm); + + delta = std::max(delta, delta_0); + + Real theoretical_damage = std::min(delta / delta_c, 1.); + + if (Math::are_float_equal(theoretical_damage, 1.)) + theoretical_traction.clear(); + else { + theoretical_traction = tangential_opening; + theoretical_traction *= beta2_kappa; + theoretical_traction += normal_opening; + theoretical_traction *= sigma_c / delta * (1. - theoretical_damage); + } + + Vector theoretical_traction_rotated(spatial_dimension); + theoretical_traction_rotated.mul(rotation, theoretical_traction); + + // adjust damage + theoretical_damage = std::max((delta - delta_0) / (delta_c - delta_0), 0.); + theoretical_damage = std::min(theoretical_damage, 1.); + + for (ghost_type_t::iterator gt = ghost_type_t::begin(); + gt != ghost_type_t::end(); + ++gt) { + + GhostType ghost_type = *gt; + Mesh::type_iterator it + = mesh.firstType(spatial_dimension, ghost_type, _ek_cohesive); + Mesh::type_iterator last + = mesh.lastType(spatial_dimension, ghost_type, _ek_cohesive); + + for (; it != last; ++it) { + ElementType type = *it; + + const Array & traction = mat_cohesive.getTraction(type, ghost_type); + const Array & damage = mat_cohesive.getDamage(type, ghost_type); + + UInt nb_quad_per_el + = model.getFEEngine("CohesiveFEEngine").getNbIntegrationPoints(type); + UInt nb_element = model.getMesh().getNbElement(type, ghost_type); + UInt tot_nb_quad = nb_element * nb_quad_per_el; + + for (UInt q = 0; q < tot_nb_quad; ++q) { + for (UInt dim = 0; dim < spatial_dimension; ++dim) { + if (!Math::are_float_equal(std::abs(theoretical_traction_rotated(dim)), + std::abs(traction(q, dim)))) { + std::cout << "Error: tractions are incorrect" << std::endl; + return 1; + } + } + + if (ghost_type == _not_ghost) + if (!Math::are_float_equal(theoretical_damage, damage(q))) { + std::cout << "Error: damage is incorrect" << std::endl; + return 1; + } + } + } + } + + return 0; +} + +/* -------------------------------------------------------------------------- */ + +void findNodesToCheck(const Mesh & mesh, + const ElementTypeMapArray & elements, + Array & nodes_to_check, + Int psize) { + + StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + Int prank = comm.whoAmI(); + + nodes_to_check.resize(0); + + Array global_nodes_to_check; + + UInt spatial_dimension = mesh.getSpatialDimension(); + const Array & position = mesh.getNodes(); + + UInt nb_nodes = position.getSize(); + + Array checked_nodes(nb_nodes); + checked_nodes.clear(); + + Mesh::type_iterator it = mesh.firstType(spatial_dimension); + Mesh::type_iterator last = mesh.lastType(spatial_dimension); + + for (; it != last; ++it) { + ElementType type = *it; + + const Array & elem = elements(type); + const Array & connectivity = mesh.getConnectivity(type); + UInt nb_nodes_per_elem = connectivity.getNbComponent(); + + for (UInt el = 0; el < elem.getSize(); ++el) { + + UInt element = elem(el); + Vector conn_el(connectivity.storage() + nb_nodes_per_elem * element, + nb_nodes_per_elem); + + for (UInt n = 0; n < nb_nodes_per_elem; ++n) { + UInt node = conn_el(n); + if (std::abs(position(node, 0) - 0.) < 1.e-6 + && !checked_nodes(node)) { + checked_nodes(node) = true; + nodes_to_check.push_back(node); + global_nodes_to_check.push_back(mesh.getNodeGlobalId(node)); + } + } + } + } + + std::vector requests; + + for (Int p = prank + 1; p < psize; ++p) { + requests.push_back(comm.asyncSend(global_nodes_to_check.storage(), + global_nodes_to_check.getSize(), + p, prank)); + } + + Array recv_nodes; + + for (Int p = 0; p < prank; ++p) { + CommunicationStatus status; + comm.probe(p, p, status); + + UInt recv_nodes_size = recv_nodes.getSize(); + recv_nodes.resize(recv_nodes_size + status.getSize()); + + comm.receive(recv_nodes.storage() + recv_nodes_size, status.getSize(), p, p); + } + + comm.waitAll(requests); + comm.freeCommunicationRequest(requests); + + for (UInt i = 0; i < recv_nodes.getSize(); ++i) { + Array::iterator node_position + = std::find(global_nodes_to_check.begin(), + global_nodes_to_check.end(), + recv_nodes(i)); + + if (node_position != global_nodes_to_check.end()) { + UInt index = node_position - global_nodes_to_check.begin(); + nodes_to_check.erase(index); + global_nodes_to_check.erase(index); + } + } +} + +/* -------------------------------------------------------------------------- */ + +bool checkEquilibrium(const Mesh & mesh, + const Array & residual) { + + UInt spatial_dimension = residual.getNbComponent(); + + Vector residual_sum(spatial_dimension); + residual_sum.clear(); + + Array::const_iterator > res_it + = residual.begin(spatial_dimension); + + for (UInt n = 0; n < residual.getSize(); ++n, ++res_it) { + if (mesh.isLocalOrMasterNode(n)) + residual_sum += *res_it; + } + + StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + comm.allReduce(residual_sum.storage(), spatial_dimension, _so_sum); + + for (UInt s = 0; s < spatial_dimension; ++s) { + if (!Math::are_float_equal(residual_sum(s), 0.)) { + + if (comm.whoAmI() == 0) + std::cout << "Error: system is not in equilibrium!" << std::endl; + + return 1; + } + } + + return 0; +} + +/* -------------------------------------------------------------------------- */ + +bool checkResidual(const Array & residual, + const Vector & traction, + const Array & nodes_to_check, + const Matrix & rotation) { + + UInt spatial_dimension = residual.getNbComponent(); + + Vector total_force(spatial_dimension); + total_force.clear(); + + for (UInt n = 0; n < nodes_to_check.getSize(); ++n) { + UInt node = nodes_to_check(n); + + Vector res(residual.storage() + node * spatial_dimension, + spatial_dimension); + + total_force += res; + } + + StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + comm.allReduce(total_force.storage(), spatial_dimension, _so_sum); + + Vector theoretical_total_force(spatial_dimension); + theoretical_total_force.mul(rotation, traction); + theoretical_total_force *= -1 * 2 * 2; + + for (UInt s = 0; s < spatial_dimension; ++s) { + if (!Math::are_float_equal(total_force(s), theoretical_total_force(s))) { + + if (comm.whoAmI() == 0) + std::cout << "Error: total force isn't correct!" << std::endl; + + return 1; + } + } + + return 0; +} + +/* -------------------------------------------------------------------------- */ + +void findElementsToDisplace(const Mesh & mesh, + ElementTypeMapArray & elements) { + UInt spatial_dimension = mesh.getSpatialDimension(); + + mesh.initElementTypeMapArray(elements, 1, spatial_dimension); + + Vector bary(spatial_dimension); + + for (ghost_type_t::iterator gt = ghost_type_t::begin(); + gt != ghost_type_t::end(); + ++gt) { + + GhostType ghost_type = *gt; + Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); + Mesh::type_iterator last = mesh.lastType(spatial_dimension, ghost_type); + + for (; it != last; ++it) { + ElementType type = *it; + + Array & elem = elements(type, ghost_type); + UInt nb_element = mesh.getNbElement(type, ghost_type); + + for (UInt el = 0; el < nb_element; ++el) { + mesh.getBarycenter(el, type, bary.storage(), ghost_type); + if (bary(0) > 0.0001) elem.push_back(el); + } + } + } +} diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/test_cohesive_parallel_intrinsic_tetrahedron.sh b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/test_cohesive_parallel_intrinsic_tetrahedron.sh new file mode 100755 index 000000000..538722641 --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/test_cohesive_parallel_intrinsic_tetrahedron.sh @@ -0,0 +1,3 @@ +#!/bin/bash + +mpirun -np 8 ./test_cohesive_parallel_intrinsic_tetrahedron diff --git a/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/tetrahedron.geo b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/tetrahedron.geo new file mode 100644 index 000000000..3e0a4c4fa --- /dev/null +++ b/extra_packages/parallel-cohesive-element/test/test_parallel_cohesive/test_cohesive_parallel_intrinsic/tetrahedron.geo @@ -0,0 +1,13 @@ +lc = 1; +h = 1; + +Point(1) = {lc, lc, lc, h}; + +line[] = Extrude{-2*lc,0,0}{ Point{1}; }; +surface[] = Extrude{0,-2*lc,0}{ Line{line[1]}; }; +volume[] = Extrude{0,0,-2*lc}{ Surface{surface[1]}; }; + +Physical Volume(1) = {volume[1]}; + +Transfinite Surface "*"; +Transfinite Volume "*";