diff --git a/src/geometry/mesh_sphere_intersector.hh b/src/geometry/mesh_sphere_intersector.hh index f14be7f34..31d0ba62d 100644 --- a/src/geometry/mesh_sphere_intersector.hh +++ b/src/geometry/mesh_sphere_intersector.hh @@ -1,104 +1,101 @@ /** * @file mesh_sphere_intersector.hh * * @author Clement Roux-Langlois * * @date creation: Wed Jun 10 2015 * * @brief Computation of mesh intersection with sphere(s) * * @section LICENSE * * Copyright (©) 2010-2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_SPHERE_INTERSECTOR_HH__ #define __AKANTU_MESH_SPHERE_INTERSECTOR_HH__ #include "aka_common.hh" #include "mesh_geom_intersector.hh" #include "mesh_geom_common.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /// Here, we know what kernel we have to use typedef Spherical SK; template class MeshSphereIntersector : public MeshGeomIntersector, SK::Sphere_3, SK> { /// Parent class type typedef MeshGeomIntersector, SK::Sphere_3, SK> parent_type; /// Result of intersection function type typedef typename IntersectionTypeHelper, K>, K::Segment_3>::intersection_type result_type; /// Pair of intersection points and element id typedef std::pair pair_type; public: /// Construct from mesh explicit MeshSphereIntersector(Mesh & mesh); /// Destructor virtual ~MeshSphereIntersector(); public: /// Construct the primitive tree object virtual void constructData(GhostType ghost_type = _not_ghost); /** * @brief Computes the intersection of the mesh with a sphere * * @param query (sphere) to compute the intersections with the mesh */ virtual void computeIntersectionQuery(const SK::Sphere_3 & query){ AKANTU_DEBUG_ERROR("This function is not implemented for spheres (It was to generic and has been replaced by computeMeshQueryIntersectionPoint"); } /// Compute intersection points between the mesh and a query virtual void computeMeshQueryIntersectionPoint(const SK::Sphere_3 & query, UInt nb_old_nodes); /// Build the IGFEM mesh virtual void buildResultFromQueryList(const std::list & query){ AKANTU_DEBUG_ERROR("This function is no longer implemented to split geometrical operations and dedicated result construction"); } /// Set the tolerance void setToleranceIntersectionOnNode(UInt tol) { this->tol_intersection_on_node = tol; } protected: /// tolerance for which the intersection is considered on the mesh node (relative to the segment lenght) Real tol_intersection_on_node; - /// number of primitive in an element of the template type - UInt nb_prim_by_el; - }; __END_AKANTU__ #include "mesh_sphere_intersector_tmpl.hh" #endif // __AKANTU_MESH_SPHERE_INTERSECTOR_HH__ diff --git a/src/geometry/mesh_sphere_intersector_tmpl.hh b/src/geometry/mesh_sphere_intersector_tmpl.hh index 2a075806b..4edb4ee5b 100644 --- a/src/geometry/mesh_sphere_intersector_tmpl.hh +++ b/src/geometry/mesh_sphere_intersector_tmpl.hh @@ -1,192 +1,190 @@ /** * @file mesh_sphere_intersector_tmpl.hh * * @author Clément Roux-Langlois * * @date creation: Wed june 10 2015 * @date last modification: Wed June 17 2015 * * @brief Computation of mesh intersection with spheres * * @section LICENSE * * Copyright (©) 2010-2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_SPHERE_INTERSECTOR_TMPL_HH__ #define __AKANTU_MESH_SPHERE_INTERSECTOR_TMPL_HH__ #include "aka_common.hh" #include "mesh_geom_common.hh" #include "tree_type_helper.hh" #include "mesh_sphere_intersector.hh" #include "static_communicator.hh" __BEGIN_AKANTU__ template MeshSphereIntersector::MeshSphereIntersector(Mesh & mesh): parent_type(mesh), - tol_intersection_on_node(1e-10), - nb_prim_by_el(0) + tol_intersection_on_node(1e-10) { this->intersection_points = new Array(0,dim); #if defined(AKANTU_IGFEM) if( (type == _triangle_3) || (type == _igfem_triangle_4) || (type == _igfem_triangle_5) ){ - this->nb_prim_by_el = 3; const_cast(this->nb_seg_by_el) = 3; } else { AKANTU_DEBUG_ERROR("Not ready for mesh type " << type); } #else if( (type != _triangle_3) ) AKANTU_DEBUG_ERROR("Not ready for mesh type " << *it); - this->nb_prim_by_el = 3; #endif - this->new_node_per_elem = new Array(0, 1 + nb_prim_by_el * 2); + this->new_node_per_elem = new Array(0, 1 + 4 * (dim-1)); } template MeshSphereIntersector::~MeshSphereIntersector() { delete this->new_node_per_elem; } template void MeshSphereIntersector::constructData(GhostType ghost_type) { this->new_node_per_elem->resize(this->mesh.getNbElement(type, ghost_type)); this->new_node_per_elem->clear(); MeshGeomIntersector, SK::Sphere_3, SK>::constructData(ghost_type); } template void MeshSphereIntersector:: computeMeshQueryIntersectionPoint(const SK::Sphere_3 & query, UInt nb_old_nodes) { /// function to replace computeIntersectionQuery in a more generic geometry module version // The newNodeEvent is not send from this method who only compute the intersection points AKANTU_DEBUG_IN(); Array & nodes = this->mesh.getNodes(); UInt nb_node = nodes.getSize() + this->intersection_points->getSize(); // Tolerance for proximity checks should be defined by user Real global_tolerance = Math::getTolerance(); Math::setTolerance(tol_intersection_on_node); typedef boost::variant sk_inter_res; TreeTypeHelper, Spherical>::const_iterator it = this->factory.getPrimitiveList().begin(), end= this->factory.getPrimitiveList().end(); for (; it != end ; ++it) { std::list s_results; CGAL::intersection(*it, query, std::back_inserter(s_results)); if (s_results.size() == 1) { // just one point if (pair_type * pair = boost::get(&s_results.front())) { if (pair->second == 1) { // not a point tangent to the sphere // Addition of the new node Vector new_node(dim, 0.0); Cartesian::Point_3 point(CGAL::to_double(pair->first.x()), CGAL::to_double(pair->first.y()), CGAL::to_double(pair->first.z())); for (UInt i = 0 ; i < dim ; i++) { new_node(i) = point[i]; } bool is_on_mesh = false, is_new = true; // check if we already compute this intersection for a neighboor element UInt n = nb_old_nodes; Array::vector_iterator existing_node = nodes.begin(dim); for (; n < nodes.getSize() ; ++n) { if (Math::are_vector_equal(dim, new_node.storage(), existing_node[n].storage())) { is_new = false; break; } } if(is_new){ Array::vector_iterator intersection_points_it = this->intersection_points->begin(dim); Array::vector_iterator intersection_points_end = this->intersection_points->end(dim); for (; intersection_points_it != intersection_points_end ; ++intersection_points_it, ++n) { if (Math::are_vector_equal(dim, new_node.storage(), intersection_points_it->storage())) { is_new = false; break; } } } Cartesian::Point_3 source_cgal(CGAL::to_double(it->source().x()), CGAL::to_double(it->source().y()), CGAL::to_double(it->source().z())); Cartesian::Point_3 target_cgal(CGAL::to_double(it->target().x()), CGAL::to_double(it->target().y()), CGAL::to_double(it->target().z())); Vector source(dim), target(dim); for (UInt i = 0 ; i < dim ; i++) { source(i) = source_cgal[i]; target(i) = target_cgal[i]; } // Check if we are close from a node of the segment if (Math::are_vector_equal(dim, source.storage(), new_node.storage()) || Math::are_vector_equal(dim, target.storage(), new_node.storage())) { is_on_mesh = true; is_new = false; } if (is_new) { this->intersection_points->push_back(new_node); nb_node++; } // deduce ghost type and element id UInt element_id = it->id(); if (!is_on_mesh) { - (*this->new_node_per_elem)(element_id, 0) += 1; - (*this->new_node_per_elem)(element_id, (2 * (*this->new_node_per_elem)(element_id, 0)) - 1) = n; - (*this->new_node_per_elem)(element_id, 2 * (*this->new_node_per_elem)(element_id, 0)) = it->segId(); + UInt & nb_new_nodes_per_el = (*this->new_node_per_elem)(element_id, 0); + nb_new_nodes_per_el += 1; + (*this->new_node_per_elem)(element_id, (2 * nb_new_nodes_per_el) - 1) = n; + (*this->new_node_per_elem)(element_id, 2 * nb_new_nodes_per_el) = it->segId(); } else { // if intersection is at a node, write node number (in el) in pennultimate position if (Math::are_vector_equal(dim, source.storage(), new_node.storage())) { (*this->new_node_per_elem)(element_id, ((*this->new_node_per_elem).getNbComponent() - 2)) = it->segId(); } else { (*this->new_node_per_elem)(element_id, ((*this->new_node_per_elem).getNbComponent() - 2)) = - (it->segId()+1) % this->nb_prim_by_el; + (it->segId()+1) % this->nb_seg_by_el; } } } } } } Math::setTolerance(global_tolerance); AKANTU_DEBUG_OUT(); } __END_AKANTU__ #endif // __AKANTU_MESH_SPHERE_INTERSECTOR_TMPL_HH__ diff --git a/src/mesh_utils/cohesive_element_inserter.cc b/src/mesh_utils/cohesive_element_inserter.cc index 44f7e7413..6fa11dc03 100644 --- a/src/mesh_utils/cohesive_element_inserter.cc +++ b/src/mesh_utils/cohesive_element_inserter.cc @@ -1,429 +1,426 @@ /** * @file cohesive_element_inserter.cc * * @author Marco Vocialta * * @date creation: Wed Dec 04 2013 * @date last modification: Tue Jul 29 2014 * * @brief Cohesive element inserter functions * * @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 #include #include "cohesive_element_inserter.hh" #include "element_group.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ CohesiveElementInserter::CohesiveElementInserter(Mesh & mesh, bool is_extrinsic, DistributedSynchronizer * synchronizer, const ID & id) : id(id), mesh(mesh), mesh_facets(mesh.initMeshFacets()), insertion_facets("insertion_facets", id), insertion_limits(mesh.getSpatialDimension(), 2), check_facets("check_facets", id) { MeshUtils::buildAllFacets(mesh, mesh_facets, 0, synchronizer); init(is_extrinsic); } /* -------------------------------------------------------------------------- */ CohesiveElementInserter::~CohesiveElementInserter() { #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) delete global_ids_updater; #endif } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::init(bool is_extrinsic) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); MeshUtils::resetFacetToDouble(mesh_facets); /// initialize facet insertion array mesh_facets.initElementTypeMapArray(insertion_facets, 1, spatial_dimension - 1, false, _ek_regular, true); /// init insertion limits for (UInt dim = 0; dim < spatial_dimension; ++dim) { insertion_limits(dim, 0) = std::numeric_limits::max() * (-1.); insertion_limits(dim, 1) = std::numeric_limits::max(); } if (is_extrinsic) { mesh_facets.initElementTypeMapArray(check_facets, 1, spatial_dimension - 1); initFacetsCheck(); } #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) facet_synchronizer = NULL; global_ids_updater = NULL; #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::initFacetsCheck() { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType facet_gt = *gt; Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, facet_gt); Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1, facet_gt); for (; it != last; ++it) { ElementType facet_type = *it; Array & f_check = check_facets(facet_type, facet_gt); const Array< std::vector > & element_to_facet = mesh_facets.getElementToSubelement(facet_type, facet_gt); UInt nb_facet = element_to_facet.getSize(); f_check.resize(nb_facet); for (UInt f = 0; f < nb_facet; ++f) { if (element_to_facet(f)[1] == ElementNull || (element_to_facet(f)[0].ghost_type == _ghost && element_to_facet(f)[1].ghost_type == _ghost)) { f_check(f) = false; } else f_check(f) = true; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::limitCheckFacets() { AKANTU_DEBUG_IN(); 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 = check_facets(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); UInt coord_in_limit = 0; while (coord_in_limit < spatial_dimension && bary_facet(coord_in_limit) > insertion_limits(coord_in_limit, 0) && bary_facet(coord_in_limit) < insertion_limits(coord_in_limit, 1)) ++coord_in_limit; if (coord_in_limit != spatial_dimension) f_check(f) = false; } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::setLimit(SpacialDirection axis, Real first_limit, Real second_limit) { AKANTU_DEBUG_ASSERT(axis < mesh.getSpatialDimension(), "You are trying to limit insertion in a direction that doesn't exist"); insertion_limits(axis, 0) = std::min(first_limit, second_limit); insertion_limits(axis, 1) = std::max(first_limit, second_limit); } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::insertIntrinsicElements() { AKANTU_DEBUG_IN(); 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) { const ElementType type_facet = *it; Array & f_insertion = insertion_facets(type_facet, ghost_type); Array > & element_to_facet = mesh_facets.getElementToSubelement(type_facet, ghost_type); UInt nb_facet = mesh_facets.getNbElement(type_facet, ghost_type); for (UInt f = 0; f < nb_facet; ++f) { if (element_to_facet(f)[1] == ElementNull) continue; mesh_facets.getBarycenter(f, type_facet, bary_facet.storage(), ghost_type); UInt coord_in_limit = 0; while (coord_in_limit < spatial_dimension && bary_facet(coord_in_limit) > insertion_limits(coord_in_limit, 0) && bary_facet(coord_in_limit) < insertion_limits(coord_in_limit, 1)) ++coord_in_limit; if (coord_in_limit == spatial_dimension) f_insertion(f) = true; } } } insertElements(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::insertIntrinsicElements(std::string physname, UInt material_index) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); ElementTypeMapArray * phys_data; try { phys_data = &(mesh_facets.getData("physical_names")); } catch(...){ phys_data = &(mesh_facets.registerData("physical_names")); mesh_facets.initElementTypeMapArray(*phys_data, 1, spatial_dimension-1, false, _ek_regular, true); } Vector bary_facet(spatial_dimension); mesh_facets.createElementGroup(physname); GhostType ghost_type = _not_ghost; 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) { const ElementType type_facet = *it; Array & f_insertion = insertion_facets(type_facet, ghost_type); Array > & element_to_facet = mesh_facets.getElementToSubelement(type_facet, ghost_type); UInt nb_facet = mesh_facets.getNbElement(type_facet, ghost_type); UInt coord_in_limit = 0; ElementGroup & group = mesh.getElementGroup(physname); ElementGroup & group_facet = mesh_facets.getElementGroup(physname); Vector bary_physgroup(spatial_dimension); Real norm_bary; for(ElementGroup::const_element_iterator el_it(group.element_begin (type_facet, ghost_type)); el_it!= group.element_end(type_facet, ghost_type); ++el_it) { UInt e = *el_it; mesh.getBarycenter(e, type_facet, bary_physgroup.storage(), ghost_type); bool find_a_partner = false; norm_bary = bary_physgroup.norm(); Array & material_id = (*phys_data)(type_facet, ghost_type); for (UInt f = 0; f < nb_facet; ++f) { if (element_to_facet(f)[1] == ElementNull) continue; mesh_facets.getBarycenter(f, type_facet, bary_facet.storage(), ghost_type); coord_in_limit = 0; while (coord_in_limit < spatial_dimension && (std::abs(bary_facet(coord_in_limit) - bary_physgroup(coord_in_limit))/norm_bary < Math::getTolerance())) ++coord_in_limit; if (coord_in_limit == spatial_dimension) { f_insertion(f) = true; find_a_partner = true; group_facet.add(type_facet, f, ghost_type,false); material_id(f) = material_index; break; } } AKANTU_DEBUG_ASSERT(find_a_partner, "The element nO " << e << " of physical group " << physname << " did not find its associated facet!" << " Try to decrease math tolerance. " << std::endl); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt CohesiveElementInserter::insertElements(bool only_double_facets) { NewNodesEvent node_event; node_event.getList().extendComponentsInterlaced(2, 1); NewElementsEvent element_event; UInt nb_new_elements = MeshUtils::insertCohesiveElements(mesh, mesh_facets, insertion_facets, node_event.getList(), element_event.getList(), only_double_facets); UInt nb_new_nodes = node_event.getList().getSize(); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (mesh.getNodesType().getSize()) { /// update nodes type updateNodesType(mesh, node_event); updateNodesType(mesh_facets, node_event); /// update global ids nb_new_nodes = updateGlobalIDs(node_event); /// compute total number of new elements StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); comm.allReduce(&nb_new_elements, 1, _so_sum); } #endif - if (nb_new_nodes > 0) { - mesh.nb_global_nodes += nb_new_nodes; - mesh_facets.nb_global_nodes += nb_new_nodes; + if (nb_new_nodes > 0) mesh.sendEvent(node_event); - } if (nb_new_elements > 0) { updateInsertionFacets(); mesh.updateTypesOffsets(_not_ghost); mesh.sendEvent(element_event); MeshUtils::resetFacetToDouble(mesh_facets); } return nb_new_elements; } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::updateInsertionFacets() { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType facet_gt = *gt; Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, facet_gt); Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1, facet_gt); for (; it != last; ++it) { ElementType facet_type = *it; Array & ins_facets = insertion_facets(facet_type, facet_gt); // this is the extrinsic case if (check_facets.exists(facet_type, facet_gt)) { Array & f_check = check_facets(facet_type, facet_gt); UInt nb_facets = f_check.getSize(); for (UInt f = 0; f < ins_facets.getSize(); ++f) { if (ins_facets(f)) { ++nb_facets; ins_facets(f) = false; f_check(f) = false; } } f_check.resize(nb_facets); } // and this the intrinsic one else { ins_facets.resize(mesh_facets.getNbElement(facet_type, facet_gt)); ins_facets.set(false); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "CohesiveElementInserter [" << std::endl; stream << space << " + mesh [" << std::endl; mesh.printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + mesh_facets [" << std::endl; mesh_facets.printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } __END_AKANTU__ diff --git a/src/mesh_utils/global_ids_updater_inline_impl.cc b/src/mesh_utils/global_ids_updater_inline_impl.cc index ab5f5a5a4..05e7909a0 100644 --- a/src/mesh_utils/global_ids_updater_inline_impl.cc +++ b/src/mesh_utils/global_ids_updater_inline_impl.cc @@ -1,120 +1,120 @@ /** * @file global_ids_updater_inline_impl.cc * @author Marco Vocialta * @date Fri Oct 2 13:33:44 2015 * * @brief Implementation of the inline functions of GlobalIdsUpdater * * @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 . * */ /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ inline UInt GlobalIdsUpdater::getNbDataForElements(const Array & elements, SynchronizationTag tag) const { UInt size = 0; if (elements.getSize() == 0) return size; if (tag == _gst_giu_global_conn) size += Mesh::getNbNodesPerElementList(elements) * sizeof(UInt); return size; } /* -------------------------------------------------------------------------- */ inline void GlobalIdsUpdater::packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const { if (tag == _gst_giu_global_conn) packUnpackGlobalConnectivity(buffer, elements); } /* -------------------------------------------------------------------------- */ inline void GlobalIdsUpdater::unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) { if (tag == _gst_giu_global_conn) packUnpackGlobalConnectivity(buffer, elements); } /* -------------------------------------------------------------------------- */ template inline void GlobalIdsUpdater::packUnpackGlobalConnectivity(CommunicationBuffer & buffer, const Array & elements) const { AKANTU_DEBUG_IN(); ElementType current_element_type = _not_defined; GhostType current_ghost_type = _casper; - Array::vector_iterator conn_begin; + Array::const_vector_iterator conn_begin; UInt nb_nodes_per_elem = 0; UInt index; Array & global_nodes_ids = mesh.getGlobalNodesIds(); Array::const_scalar_iterator it = elements.begin(); Array::const_scalar_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; - Array & connectivity = mesh.getConnectivity(current_element_type, - current_ghost_type); + const Array & connectivity = mesh.getConnectivity(current_element_type, + current_ghost_type); nb_nodes_per_elem = connectivity.getNbComponent(); conn_begin = connectivity.begin(nb_nodes_per_elem); } /// get element connectivity - Vector current_conn = conn_begin[el.element]; + const Vector current_conn = conn_begin[el.element]; /// loop on all connectivity nodes for (UInt n = 0; n < nb_nodes_per_elem; ++n) { UInt node = current_conn(n); if (pack_mode) { /// if node is local or master pack its global id, otherwise /// dummy data if (mesh.isLocalOrMasterNode(node)) - index = mesh.getNodeGlobalId(node); + index = global_nodes_ids(node); else index = UInt(-1); buffer << index; } else { buffer >> index; /// update slave nodes' index if (index != UInt(-1) && mesh.isSlaveNode(node)) global_nodes_ids(node) = index; } } } AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/src/mesh_utils/mesh_utils.cc b/src/mesh_utils/mesh_utils.cc index b83a478c3..3900669ec 100644 --- a/src/mesh_utils/mesh_utils.cc +++ b/src/mesh_utils/mesh_utils.cc @@ -1,2198 +1,2196 @@ /** * @file mesh_utils.cc * * @author Guillaume Anciaux * @author Leonardo Snozzi * @author Marco Vocialta * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Fri Aug 20 2010 * @date last modification: Mon Jun 09 2014 * * @brief All mesh utils necessary for various tasks * * @section LICENSE * * Copyright (©) 2010-2012, 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 "mesh_utils.hh" #include "aka_safe_enum.hh" #include "fe_engine.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ void MeshUtils::buildNode2Elements(const Mesh & mesh, CSR & node_to_elem, UInt spatial_dimension) { AKANTU_DEBUG_IN(); if (spatial_dimension == _all_dimensions) spatial_dimension = mesh.getSpatialDimension(); /// count number of occurrence of each node UInt nb_nodes = mesh.getNbNodes(); /// array for the node-element list node_to_elem.resizeRows(nb_nodes); node_to_elem.clearRows(); AKANTU_DEBUG_ASSERT(mesh.firstType(spatial_dimension) != mesh.lastType(spatial_dimension), "Some elements must be found in right dimension to compute facets!"); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { Mesh::type_iterator first = mesh.firstType(spatial_dimension, *gt, _ek_not_defined); Mesh::type_iterator last = mesh.lastType(spatial_dimension, *gt, _ek_not_defined); for (; first != last; ++first) { ElementType type = *first; UInt nb_element = mesh.getNbElement(type, *gt); Array::const_iterator< Vector > conn_it = mesh.getConnectivity(type, *gt).begin(Mesh::getNbNodesPerElement(type)); for (UInt el = 0; el < nb_element; ++el, ++conn_it) for (UInt n = 0; n < conn_it->size(); ++n) ++node_to_elem.rowOffset((*conn_it)(n)); } } node_to_elem.countToCSR(); node_to_elem.resizeCols(); /// rearrange element to get the node-element list Element e; node_to_elem.beginInsertions(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { Mesh::type_iterator first = mesh.firstType(spatial_dimension, *gt, _ek_not_defined); Mesh::type_iterator last = mesh.lastType(spatial_dimension, *gt, _ek_not_defined); e.ghost_type = *gt; for (; first != last; ++first) { ElementType type = *first; e.type = type; e.kind = Mesh::getKind(type); UInt nb_element = mesh.getNbElement(type, *gt); Array::const_iterator< Vector > conn_it = mesh.getConnectivity(type, *gt).begin(Mesh::getNbNodesPerElement(type)); for (UInt el = 0; el < nb_element; ++el, ++conn_it) { e.element = el; for (UInt n = 0; n < conn_it->size(); ++n) node_to_elem.insertInRow((*conn_it)(n), e); } } } node_to_elem.endInsertions(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * This function should disappear in the future (used in mesh partitioning) */ void MeshUtils::buildNode2Elements(const Mesh & mesh, CSR & node_to_elem, UInt spatial_dimension) { AKANTU_DEBUG_IN(); if (spatial_dimension == _all_dimensions) spatial_dimension = mesh.getSpatialDimension(); UInt nb_nodes = mesh.getNbNodes(); const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; UInt nb_types = type_list.size(); UInt nb_good_types = 0; UInt nb_nodes_per_element[nb_types]; UInt * conn_val[nb_types]; UInt nb_element[nb_types]; for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; if(Mesh::getSpatialDimension(type) != spatial_dimension) continue; nb_nodes_per_element[nb_good_types] = Mesh::getNbNodesPerElement(type); conn_val[nb_good_types] = mesh.getConnectivity(type, _not_ghost).storage(); nb_element[nb_good_types] = mesh.getConnectivity(type, _not_ghost).getSize(); nb_good_types++; } AKANTU_DEBUG_ASSERT(nb_good_types != 0, "Some elements must be found in right dimension to compute facets!"); /// array for the node-element list node_to_elem.resizeRows(nb_nodes); node_to_elem.clearRows(); /// count number of occurrence of each node for (UInt t = 0; t < nb_good_types; ++t) { for (UInt el = 0; el < nb_element[t]; ++el) { UInt el_offset = el*nb_nodes_per_element[t]; for (UInt n = 0; n < nb_nodes_per_element[t]; ++n) { ++node_to_elem.rowOffset(conn_val[t][el_offset + n]); } } } node_to_elem.countToCSR(); node_to_elem.resizeCols(); node_to_elem.beginInsertions(); /// rearrange element to get the node-element list for (UInt t = 0, linearized_el = 0; t < nb_good_types; ++t) for (UInt el = 0; el < nb_element[t]; ++el, ++linearized_el) { UInt el_offset = el*nb_nodes_per_element[t]; for (UInt n = 0; n < nb_nodes_per_element[t]; ++n) node_to_elem.insertInRow(conn_val[t][el_offset + n], linearized_el); } node_to_elem.endInsertions(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildNode2ElementsElementTypeMap(const Mesh & mesh, CSR & node_to_elem, const ElementType & type, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_elements = mesh.getConnectivity(type, ghost_type).getSize(); UInt * conn_val = mesh.getConnectivity(type, ghost_type).storage(); /// array for the node-element list node_to_elem.resizeRows(nb_nodes); node_to_elem.clearRows(); /// count number of occurrence of each node for (UInt el = 0; el < nb_elements; ++el) { UInt el_offset = el*nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) ++node_to_elem.rowOffset(conn_val[el_offset + n]); } /// convert the occurrence array in a csr one node_to_elem.countToCSR(); node_to_elem.resizeCols(); node_to_elem.beginInsertions(); /// save the element index in the node-element list for (UInt el = 0; el < nb_elements; ++el) { UInt el_offset = el*nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) { node_to_elem.insertInRow(conn_val[el_offset + n], el); } } node_to_elem.endInsertions(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildFacets(Mesh & mesh){ AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType gt_facet = *gt; Mesh::type_iterator it = mesh.firstType(spatial_dimension - 1, gt_facet); Mesh::type_iterator end = mesh.lastType(spatial_dimension - 1, gt_facet); for(; it != end; ++it) { mesh.getConnectivity(*it, *gt).resize(0); // \todo inform the mesh event handler } } buildFacetsDimension(mesh, mesh, true, spatial_dimension); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildAllFacets(const Mesh & mesh, Mesh & mesh_facets, UInt to_dimension, DistributedSynchronizer * synchronizer) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); buildAllFacets(mesh, mesh_facets, spatial_dimension, to_dimension, synchronizer); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildAllFacets(const Mesh & mesh, Mesh & mesh_facets, UInt from_dimension, UInt to_dimension, DistributedSynchronizer * synchronizer) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(mesh_facets.isMeshFacets(), "The mesh_facets should be initialized with initMeshFacets"); const ElementTypeMapArray * prank_to_element = NULL; if (synchronizer) { synchronizer->buildPrankToElement(); prank_to_element = &synchronizer->getPrankToElement(); } /// generate facets buildFacetsDimension(mesh, mesh_facets, false, from_dimension, prank_to_element); /// copy nodes type mesh_facets.nodes_type.resize(mesh.nodes_type.getSize()); mesh_facets.nodes_type.copy(mesh.nodes_type); /// sort facets and generate subfacets for (UInt i = from_dimension - 1; i > to_dimension; --i) { buildFacetsDimension(mesh_facets, mesh_facets, false, i, prank_to_element); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildFacetsDimension(const Mesh & mesh, Mesh & mesh_facets, bool boundary_only, UInt dimension, const ElementTypeMapArray * prank_to_element){ AKANTU_DEBUG_IN(); // save the current parent of mesh_facets and set it temporarly to mesh since // mesh is the one containing the elements for which mesh_facets has the subelements // example: if the function is called with mesh = mesh_facets const Mesh & mesh_facets_parent = mesh_facets.getMeshParent(); mesh_facets.defineMeshParent(mesh); UInt spatial_dimension = mesh.getSpatialDimension(); const Array & mesh_facets_nodes = mesh_facets.getNodes(); const Array::const_vector_iterator mesh_facets_nodes_it = mesh_facets_nodes.begin(spatial_dimension); CSR node_to_elem; buildNode2Elements(mesh, node_to_elem, dimension); Array counter; std::vector connected_elements; // init the SubelementToElement data to improve performance for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; Mesh::type_iterator first = mesh.firstType(dimension, ghost_type); Mesh::type_iterator last = mesh.lastType(dimension, ghost_type); for(; first != last; ++first) { ElementType type = *first; mesh_facets.getSubelementToElementPointer(type, ghost_type); Vector facet_types = mesh.getAllFacetTypes(type); for (UInt ft = 0; ft < facet_types.size(); ++ft) { ElementType facet_type = facet_types(ft); mesh_facets.getElementToSubelementPointer(facet_type, ghost_type); mesh_facets.getConnectivityPointer(facet_type, ghost_type); } } } Element current_element; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; GhostType facet_ghost_type = ghost_type; current_element.ghost_type = ghost_type; Mesh::type_iterator first = mesh.firstType(dimension, ghost_type); Mesh::type_iterator last = mesh.lastType(dimension, ghost_type); for(; first != last; ++first) { ElementType type = *first; Vector facet_types = mesh.getAllFacetTypes(type); current_element.type = type; for (UInt ft = 0; ft < facet_types.size(); ++ft) { ElementType facet_type = facet_types(ft); UInt nb_element = mesh.getNbElement(type, ghost_type); Array< std::vector > * element_to_subelement = &mesh_facets.getElementToSubelement(facet_type, ghost_type); Array * connectivity_facets = &mesh_facets.getConnectivity(facet_type, ghost_type); UInt nb_facet_per_element = mesh.getNbFacetsPerElement(type, ft); const Array & element_connectivity = mesh.getConnectivity(type, ghost_type); const Matrix facet_local_connectivity = mesh.getFacetLocalConnectivity(type, ft); UInt nb_nodes_per_facet = connectivity_facets->getNbComponent(); Vector facet(nb_nodes_per_facet); for (UInt el = 0; el < nb_element; ++el) { current_element.element = el; for (UInt f = 0; f < nb_facet_per_element; ++f) { for (UInt n = 0; n < nb_nodes_per_facet; ++n) facet(n) = element_connectivity(el, facet_local_connectivity(f, n)); UInt first_node_nb_elements = node_to_elem.getNbCols(facet(0)); counter.resize(first_node_nb_elements); counter.clear(); //loop over the other nodes to search intersecting elements, //which are the elements that share another node with the //starting element after first_node CSR::iterator first_node_elements = node_to_elem.begin(facet(0)); CSR::iterator first_node_elements_end = node_to_elem.end(facet(0)); UInt local_el = 0; for (; first_node_elements != first_node_elements_end; ++first_node_elements, ++local_el) { for (UInt n = 1; n < nb_nodes_per_facet; ++n) { CSR::iterator node_elements_begin = node_to_elem.begin(facet(n)); CSR::iterator node_elements_end = node_to_elem.end (facet(n)); counter(local_el) += std::count(node_elements_begin, node_elements_end, *first_node_elements); } } // counting the number of elements connected to the facets and // taking the minimum element number, because the facet should // be inserted just once UInt nb_element_connected_to_facet = 0; Element minimum_el = ElementNull; connected_elements.clear(); for (UInt el_f = 0; el_f < first_node_nb_elements; el_f++) { Element real_el = node_to_elem(facet(0), el_f); if (counter(el_f) == nb_nodes_per_facet - 1) { ++nb_element_connected_to_facet; minimum_el = std::min(minimum_el, real_el); connected_elements.push_back(real_el); } } if (minimum_el == current_element) { bool full_ghost_facet = false; UInt n = 0; while (n < nb_nodes_per_facet && mesh.isPureGhostNode(facet(n))) ++n; if (n == nb_nodes_per_facet) full_ghost_facet = true; if (!full_ghost_facet) { if (!boundary_only || (boundary_only && nb_element_connected_to_facet == 1)) { std::vector elements; // build elements_on_facets: linearized_el must come first // in order to store the facet in the correct direction // and avoid to invert the sign in the normal computation elements.push_back(current_element); /// boundary facet if (nb_element_connected_to_facet == 1) elements.push_back(ElementNull); /// internal facet else if (nb_element_connected_to_facet == 2) { elements.push_back(connected_elements[1]); /// check if facet is in between ghost and normal /// elements: if it's the case, the facet is either /// ghost or not ghost. The criterion to decide this /// is arbitrary. It was chosen to check the processor /// id (prank) of the two neighboring elements. If /// prank of the ghost element is lower than prank of /// the normal one, the facet is not ghost, otherwise /// it's ghost GhostType gt[2] = { _not_ghost, _not_ghost }; for (UInt el = 0; el < connected_elements.size(); ++el) gt[el] = connected_elements[el].ghost_type; if (gt[0] + gt[1] == 1) { if (prank_to_element) { UInt prank[2]; for (UInt el = 0; el < 2; ++el) { UInt current_el = connected_elements[el].element; ElementType current_type = connected_elements[el].type; GhostType current_gt = connected_elements[el].ghost_type; const Array & prank_to_el = (*prank_to_element)(current_type, current_gt); prank[el] = prank_to_el(current_el); } bool ghost_one = (gt[0] != _ghost); if (prank[ghost_one] > prank[!ghost_one]) facet_ghost_type = _not_ghost; else facet_ghost_type = _ghost; connectivity_facets = &mesh_facets.getConnectivity(facet_type, facet_ghost_type); element_to_subelement = &mesh_facets.getElementToSubelement(facet_type, facet_ghost_type); } } } /// facet of facet else { for (UInt i = 1; i < nb_element_connected_to_facet; ++i) { elements.push_back(connected_elements[i]); } } element_to_subelement->push_back(elements); connectivity_facets->push_back(facet); /// current facet index UInt current_facet = connectivity_facets->getSize() - 1; /// loop on every element connected to current facet and /// insert current facet in the first free spot of the /// subelement_to_element vector for (UInt elem = 0; elem < elements.size(); ++elem) { Element loc_el = elements[elem]; if (loc_el.type != _not_defined) { Array & subelement_to_element = mesh_facets.getSubelementToElement(loc_el.type, loc_el.ghost_type); UInt nb_facet_per_loc_element = subelement_to_element.getNbComponent(); for (UInt f_in = 0; f_in < nb_facet_per_loc_element; ++f_in) { if (subelement_to_element(loc_el.element, f_in).type == _not_defined) { subelement_to_element(loc_el.element, f_in).type = facet_type; subelement_to_element(loc_el.element, f_in).element = current_facet; subelement_to_element(loc_el.element, f_in).ghost_type = facet_ghost_type; break; } } } } /// reset connectivity in case a facet was found in /// between ghost and normal elements if (facet_ghost_type != ghost_type) { facet_ghost_type = ghost_type; connectivity_facets = mesh_facets.getConnectivityPointer(facet_type, facet_ghost_type); element_to_subelement = mesh_facets.getElementToSubelementPointer(facet_type, facet_ghost_type); } } } } } } } } } // restore the parent of mesh_facet mesh_facets.defineMeshParent(mesh_facets_parent); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::renumberMeshNodes(Mesh & mesh, UInt * local_connectivities, UInt nb_local_element, UInt nb_ghost_element, ElementType type, Array & old_nodes_numbers) { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); std::map renumbering_map; for (UInt i = 0; i < old_nodes_numbers.getSize(); ++i) { renumbering_map[old_nodes_numbers(i)] = i; } /// renumber the nodes renumberNodesInConnectivity(local_connectivities, (nb_local_element + nb_ghost_element)*nb_nodes_per_element, renumbering_map); std::map::iterator it = renumbering_map.begin(); std::map::iterator end = renumbering_map.end(); old_nodes_numbers.resize(renumbering_map.size()); for (;it != end; ++it) { old_nodes_numbers(it->second) = it->first; } renumbering_map.clear(); /// copy the renumbered connectivity to the right place Array * local_conn = mesh.getConnectivityPointer(type); local_conn->resize(nb_local_element); memcpy(local_conn->storage(), local_connectivities, nb_local_element * nb_nodes_per_element * sizeof(UInt)); Array * ghost_conn = mesh.getConnectivityPointer(type, _ghost); ghost_conn->resize(nb_ghost_element); memcpy(ghost_conn->storage(), local_connectivities + nb_local_element * nb_nodes_per_element, nb_ghost_element * nb_nodes_per_element * sizeof(UInt)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::renumberNodesInConnectivity(UInt * list_nodes, UInt nb_nodes, std::map & renumbering_map) { AKANTU_DEBUG_IN(); UInt * connectivity = list_nodes; UInt new_node_num = renumbering_map.size(); for (UInt n = 0; n < nb_nodes; ++n, ++connectivity) { UInt & node = *connectivity; std::map::iterator it = renumbering_map.find(node); if(it == renumbering_map.end()) { UInt old_node = node; renumbering_map[old_node] = new_node_num; node = new_node_num; ++new_node_num; } else { node = it->second; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::purifyMesh(Mesh & mesh) { AKANTU_DEBUG_IN(); std::map renumbering_map; RemovedNodesEvent remove_nodes(mesh); Array & nodes_removed = remove_nodes.getList(); for (UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type = (GhostType) gt; Mesh::type_iterator it = mesh.firstType(_all_dimensions, ghost_type, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(_all_dimensions, ghost_type, _ek_not_defined); for(; it != end; ++it) { ElementType type(*it); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const Array & connectivity_vect = mesh.getConnectivity(type, ghost_type); UInt nb_element(connectivity_vect.getSize()); UInt * connectivity = connectivity_vect.storage(); renumberNodesInConnectivity (connectivity, nb_element*nb_nodes_per_element, renumbering_map); } } Array & new_numbering = remove_nodes.getNewNumbering(); std::fill(new_numbering.begin(), new_numbering.end(), UInt(-1)); std::map::iterator it = renumbering_map.begin(); std::map::iterator end = renumbering_map.end(); for (; it != end; ++it) { new_numbering(it->first) = it->second; } for (UInt i = 0; i < new_numbering.getSize(); ++i) { if(new_numbering(i) == UInt(-1)) nodes_removed.push_back(i); } mesh.sendEvent(remove_nodes); AKANTU_DEBUG_OUT(); } #if defined(AKANTU_COHESIVE_ELEMENT) /* -------------------------------------------------------------------------- */ UInt MeshUtils::insertCohesiveElements(Mesh & mesh, Mesh & mesh_facets, const ElementTypeMapArray & facet_insertion, Array & doubled_nodes, Array & new_elements, bool only_double_facets) { UInt spatial_dimension = mesh.getSpatialDimension(); UInt elements_to_insert = updateFacetToDouble(mesh_facets, facet_insertion); if (elements_to_insert > 0) { if (spatial_dimension == 1) { doublePointFacet(mesh, mesh_facets, doubled_nodes); } else { doubleFacet(mesh, mesh_facets, spatial_dimension - 1, doubled_nodes, true); findSubfacetToDouble(mesh, mesh_facets); if (spatial_dimension == 2) { doubleSubfacet<2>(mesh, mesh_facets, doubled_nodes); } else if (spatial_dimension == 3) { doubleFacet(mesh, mesh_facets, 1, doubled_nodes, false); findSubfacetToDouble(mesh, mesh_facets); doubleSubfacet<3>(mesh, mesh_facets, doubled_nodes); } } if (!only_double_facets) updateCohesiveData(mesh, mesh_facets, new_elements); } return elements_to_insert; } #endif /* -------------------------------------------------------------------------- */ void MeshUtils::doubleNodes(Mesh & mesh, const std::vector & old_nodes, Array & doubled_nodes) { AKANTU_DEBUG_IN(); Array & position = mesh.getNodes(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt old_nb_nodes = position.getSize(); UInt new_nb_nodes = old_nb_nodes + old_nodes.size(); UInt old_nb_doubled_nodes = doubled_nodes.getSize(); UInt new_nb_doubled_nodes = old_nb_doubled_nodes + old_nodes.size(); position.resize(new_nb_nodes); doubled_nodes.resize(new_nb_doubled_nodes); Array::iterator > position_begin = position.begin(spatial_dimension); for (UInt n = 0; n < old_nodes.size(); ++n) { UInt new_node = old_nb_nodes + n; /// store doubled nodes doubled_nodes(old_nb_doubled_nodes + n, 0) = old_nodes[n]; doubled_nodes(old_nb_doubled_nodes + n, 1) = new_node; /// update position std::copy(position_begin + old_nodes[n], position_begin + old_nodes[n] + 1, position_begin + new_node); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::doubleFacet(Mesh & mesh, Mesh & mesh_facets, UInt facet_dimension, Array & doubled_nodes, bool facet_mode) { AKANTU_DEBUG_IN(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType gt_facet = *gt; Mesh::type_iterator it = mesh_facets.firstType(facet_dimension, gt_facet); Mesh::type_iterator end = mesh_facets.lastType(facet_dimension, gt_facet); for(; it != end; ++it) { ElementType type_facet = *it; Array & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); UInt nb_facet_to_double = f_to_double.getSize(); if (nb_facet_to_double == 0) continue; ElementType type_subfacet = Mesh::getFacetType(type_facet); const UInt nb_subfacet_per_facet = Mesh::getNbFacetsPerElement(type_facet); GhostType gt_subfacet = _casper; Array > * f_to_subfacet = NULL; Array & subfacet_to_facet = mesh_facets.getSubelementToElement(type_facet, gt_facet); Array & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); UInt nb_nodes_per_facet = conn_facet.getNbComponent(); UInt old_nb_facet = conn_facet.getSize(); UInt new_nb_facet = old_nb_facet + nb_facet_to_double; conn_facet.resize(new_nb_facet); subfacet_to_facet.resize(new_nb_facet); UInt new_facet = old_nb_facet - 1; Element new_facet_el(type_facet, 0, gt_facet); Array::iterator > subfacet_to_facet_begin = subfacet_to_facet.begin(nb_subfacet_per_facet); Array::iterator > conn_facet_begin = conn_facet.begin(nb_nodes_per_facet); for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { UInt old_facet = f_to_double(facet); ++new_facet; /// adding a new facet by copying original one /// copy connectivity in new facet std::copy(conn_facet_begin + old_facet, conn_facet_begin + old_facet + 1, conn_facet_begin + new_facet); /// update subfacet_to_facet std::copy(subfacet_to_facet_begin + old_facet, subfacet_to_facet_begin + old_facet + 1, subfacet_to_facet_begin + new_facet); new_facet_el.element = new_facet; /// loop on every subfacet for (UInt sf = 0; sf < nb_subfacet_per_facet; ++sf) { Element & subfacet = subfacet_to_facet(old_facet, sf); if (subfacet == ElementNull) continue; if (gt_subfacet != subfacet.ghost_type) { gt_subfacet = subfacet.ghost_type; f_to_subfacet = & mesh_facets.getElementToSubelement(type_subfacet, subfacet.ghost_type); } /// update facet_to_subfacet array (*f_to_subfacet)(subfacet.element).push_back(new_facet_el); } } /// update facet_to_subfacet and _segment_3 facets if any if (!facet_mode) { updateSubfacetToFacet(mesh_facets, type_facet, gt_facet, true); updateFacetToSubfacet(mesh_facets, type_facet, gt_facet, true); updateQuadraticSegments(mesh, mesh_facets, type_facet, gt_facet, doubled_nodes); } else updateQuadraticSegments(mesh, mesh_facets, type_facet, gt_facet, doubled_nodes); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt MeshUtils::updateFacetToDouble(Mesh & mesh_facets, const ElementTypeMapArray & facet_insertion) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh_facets.getSpatialDimension(); UInt nb_facets_to_double = 0.; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType gt_facet = *gt; Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, gt_facet); Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1, gt_facet); for(; it != end; ++it) { ElementType type_facet = *it; const Array & f_insertion = facet_insertion(type_facet, gt_facet); Array & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); Array< std::vector > & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); ElementType el_type = _not_defined; GhostType el_gt = _casper; UInt nb_facet_per_element = 0; Element old_facet_el(type_facet, 0, gt_facet); Array * facet_to_element = NULL; for (UInt f = 0; f < f_insertion.getSize(); ++f) { if (f_insertion(f) == false) continue; ++nb_facets_to_double; if (element_to_facet(f)[1].type == _not_defined #if defined(AKANTU_COHESIVE_ELEMENT) || element_to_facet(f)[1].kind == _ek_cohesive #endif ) { AKANTU_DEBUG_WARNING("attempt to double a facet on the boundary"); continue; } f_to_double.push_back(f); UInt new_facet = mesh_facets.getNbElement(type_facet, gt_facet) + f_to_double.getSize() - 1; old_facet_el.element = f; /// update facet_to_element vector Element & elem_to_update = element_to_facet(f)[1]; UInt el = elem_to_update.element; if (elem_to_update.ghost_type != el_gt || elem_to_update.type != el_type) { el_type = elem_to_update.type; el_gt = elem_to_update.ghost_type; facet_to_element = & mesh_facets.getSubelementToElement(el_type, el_gt); nb_facet_per_element = facet_to_element->getNbComponent(); } Element * f_update = std::find(facet_to_element->storage() + el * nb_facet_per_element, facet_to_element->storage() + el * nb_facet_per_element + nb_facet_per_element, old_facet_el); AKANTU_DEBUG_ASSERT(facet_to_element->storage() + el * nb_facet_per_element != facet_to_element->storage() + el * nb_facet_per_element + nb_facet_per_element, "Facet not found"); f_update->element = new_facet; /// update elements connected to facet std::vector first_facet_list = element_to_facet(f); element_to_facet.push_back(first_facet_list); /// set new and original facets as boundary facets element_to_facet(new_facet)[0] = element_to_facet(new_facet)[1]; element_to_facet(f)[1] = ElementNull; element_to_facet(new_facet)[1] = ElementNull; } } } AKANTU_DEBUG_OUT(); return nb_facets_to_double; } /* -------------------------------------------------------------------------- */ void MeshUtils::resetFacetToDouble(Mesh & mesh_facets) { AKANTU_DEBUG_IN(); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = mesh_facets.firstType(_all_dimensions, gt); Mesh::type_iterator end = mesh_facets.lastType(_all_dimensions, gt); for(; it != end; ++it) { ElementType type = *it; mesh_facets.getDataPointer("facet_to_double", type, gt, 1, false); mesh_facets.getDataPointer > ("facets_to_subfacet_double", type, gt, 1, false); mesh_facets.getDataPointer > ("elements_to_subfacet_double", type, gt, 1, false); mesh_facets.getDataPointer > ("subfacets_to_subsubfacet_double", type, gt, 1, false); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MeshUtils::findSubfacetToDouble(Mesh & mesh, Mesh & mesh_facets) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh_facets.getSpatialDimension(); if (spatial_dimension == 1) return; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType gt_facet = *gt; Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, gt_facet); Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1, gt_facet); for(; it != end; ++it) { ElementType type_facet = *it; Array & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); UInt nb_facet_to_double = f_to_double.getSize(); if (nb_facet_to_double == 0) continue; ElementType type_subfacet = Mesh::getFacetType(type_facet); GhostType gt_subfacet = _casper; ElementType type_subsubfacet = Mesh::getFacetType(type_subfacet); GhostType gt_subsubfacet = _casper; Array * conn_subfacet = NULL; Array * sf_to_double = NULL; Array > * sf_to_subfacet_double = NULL; Array > * f_to_subfacet_double = NULL; Array > * el_to_subfacet_double = NULL; UInt nb_subfacet = Mesh::getNbFacetsPerElement(type_facet); UInt nb_subsubfacet; UInt nb_nodes_per_sf_el; if (subsubfacet_mode) { nb_nodes_per_sf_el = Mesh::getNbNodesPerElement(type_subsubfacet); nb_subsubfacet = Mesh::getNbFacetsPerElement(type_subfacet); } else nb_nodes_per_sf_el = Mesh::getNbNodesPerElement(type_subfacet); Array & subfacet_to_facet = mesh_facets.getSubelementToElement(type_facet, gt_facet); Array< std::vector > & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); Array * subsubfacet_to_subfacet = NULL; UInt old_nb_facet = subfacet_to_facet.getSize() - nb_facet_to_double; Element current_facet(type_facet, 0, gt_facet); std::vector element_list; std::vector facet_list; std::vector * subfacet_list; if (subsubfacet_mode) subfacet_list = new std::vector; /// map to filter subfacets Array< std::vector > * facet_to_subfacet = NULL; /// this is used to make sure that both new and old facets are /// checked UInt facets[2]; /// loop on every facet for (UInt f_index = 0; f_index < 2; ++f_index) { for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { facets[bool(f_index)] = f_to_double(facet); facets[!bool(f_index)] = old_nb_facet + facet; UInt old_facet = facets[0]; UInt new_facet = facets[1]; Element & starting_element = element_to_facet(new_facet)[0]; current_facet.element = old_facet; /// loop on every subfacet for (UInt sf = 0; sf < nb_subfacet; ++sf) { Element & subfacet = subfacet_to_facet(old_facet, sf); if (subfacet == ElementNull) continue; if (gt_subfacet != subfacet.ghost_type) { gt_subfacet = subfacet.ghost_type; if (subsubfacet_mode) { subsubfacet_to_subfacet = & mesh_facets.getSubelementToElement(type_subfacet, gt_subfacet); } else { conn_subfacet = & mesh_facets.getConnectivity(type_subfacet, gt_subfacet); sf_to_double = & mesh_facets.getData("facet_to_double", type_subfacet, gt_subfacet); f_to_subfacet_double = & mesh_facets.getData >("facets_to_subfacet_double", type_subfacet, gt_subfacet); el_to_subfacet_double = & mesh_facets.getData >("elements_to_subfacet_double", type_subfacet, gt_subfacet); facet_to_subfacet = & mesh_facets.getElementToSubelement(type_subfacet, gt_subfacet); } } if (subsubfacet_mode) { /// loop on every subsubfacet for (UInt ssf = 0; ssf < nb_subsubfacet; ++ssf) { Element & subsubfacet = (*subsubfacet_to_subfacet)(subfacet.element, ssf); if (subsubfacet == ElementNull) continue; if (gt_subsubfacet != subsubfacet.ghost_type) { gt_subsubfacet = subsubfacet.ghost_type; conn_subfacet = & mesh_facets.getConnectivity(type_subsubfacet, gt_subsubfacet); sf_to_double = & mesh_facets.getData("facet_to_double", type_subsubfacet, gt_subsubfacet); sf_to_subfacet_double = & mesh_facets.getData >("subfacets_to_subsubfacet_double", type_subsubfacet, gt_subsubfacet); f_to_subfacet_double = & mesh_facets.getData >("facets_to_subfacet_double", type_subsubfacet, gt_subsubfacet); el_to_subfacet_double = & mesh_facets.getData >("elements_to_subfacet_double", type_subsubfacet, gt_subsubfacet); facet_to_subfacet = & mesh_facets.getElementToSubelement(type_subsubfacet, gt_subsubfacet); } UInt global_ssf = subsubfacet.element; Vector subsubfacet_connectivity(conn_subfacet->storage() + global_ssf * nb_nodes_per_sf_el, nb_nodes_per_sf_el); /// check if subsubfacet is to be doubled if (findElementsAroundSubfacet(mesh, mesh_facets, starting_element, current_facet, subsubfacet_connectivity, element_list, facet_list, subfacet_list) == false && removeElementsInVector(*subfacet_list, (*facet_to_subfacet)(global_ssf)) == false) { sf_to_double->push_back(global_ssf); sf_to_subfacet_double->push_back(*subfacet_list); f_to_subfacet_double->push_back(facet_list); el_to_subfacet_double->push_back(element_list); } } } else { const UInt global_sf = subfacet.element; Vector subfacet_connectivity(conn_subfacet->storage() + global_sf * nb_nodes_per_sf_el, nb_nodes_per_sf_el); /// check if subfacet is to be doubled if (findElementsAroundSubfacet(mesh, mesh_facets, starting_element, current_facet, subfacet_connectivity, element_list, facet_list) == false && removeElementsInVector(facet_list, (*facet_to_subfacet)(global_sf)) == false) { sf_to_double->push_back(global_sf); f_to_subfacet_double->push_back(facet_list); el_to_subfacet_double->push_back(element_list); } } } } } if (subsubfacet_mode) delete subfacet_list; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ #if defined(AKANTU_COHESIVE_ELEMENT) void MeshUtils::updateCohesiveData(Mesh & mesh, Mesh & mesh_facets, Array & new_elements) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); bool third_dimension = spatial_dimension == 3; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType gt_facet = *gt; Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, gt_facet); Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1, gt_facet); for(; it != end; ++it) { ElementType type_facet = *it; Array & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); UInt nb_facet_to_double = f_to_double.getSize(); if (nb_facet_to_double == 0) continue; ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); Array * facet_to_coh_element = mesh_facets.getSubelementToElementPointer(type_cohesive, gt_facet); Array & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); Array * conn_cohesive = mesh.getConnectivityPointer(type_cohesive, gt_facet); UInt nb_nodes_per_facet = Mesh::getNbNodesPerElement(type_facet); Array< std::vector > & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); UInt old_nb_cohesive_elements = conn_cohesive->getSize(); UInt new_nb_cohesive_elements = conn_cohesive->getSize() + nb_facet_to_double; UInt old_nb_facet = element_to_facet.getSize() - nb_facet_to_double; facet_to_coh_element->resize(new_nb_cohesive_elements); conn_cohesive->resize(new_nb_cohesive_elements); UInt new_elements_old_size = new_elements.getSize(); new_elements.resize(new_elements_old_size + nb_facet_to_double); Element c_element(type_cohesive, 0, gt_facet, _ek_cohesive); Element f_element(type_facet, 0, gt_facet); UInt facets[2]; for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { /// (in 3D cohesive elements connectivity is inverted) facets[third_dimension] = f_to_double(facet); facets[!third_dimension] = old_nb_facet + facet; UInt cohesive_element = old_nb_cohesive_elements + facet; /// store doubled facets f_element.element = facets[0]; (*facet_to_coh_element)(cohesive_element, 0) = f_element; f_element.element = facets[1]; (*facet_to_coh_element)(cohesive_element, 1) = f_element; /// modify cohesive elements' connectivity for (UInt n = 0; n < nb_nodes_per_facet; ++n) { (*conn_cohesive)(cohesive_element, n) = conn_facet(facets[0], n); (*conn_cohesive)(cohesive_element, n + nb_nodes_per_facet) = conn_facet(facets[1], n); } /// update element_to_facet vectors c_element.element = cohesive_element; element_to_facet(facets[0])[1] = c_element; element_to_facet(facets[1])[1] = c_element; /// add cohesive element to the element event list new_elements(new_elements_old_size + facet) = c_element; } } } AKANTU_DEBUG_OUT(); } #endif /* -------------------------------------------------------------------------- */ void MeshUtils::doublePointFacet(Mesh & mesh, Mesh & mesh_facets, Array & doubled_nodes) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); if (spatial_dimension != 1) return; Array & position = mesh.getNodes(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType gt_facet = *gt; Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, gt_facet); Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1, gt_facet); for(; it != end; ++it) { ElementType type_facet = *it; Array & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); Array< std::vector > & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); const Array & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); UInt nb_facet_to_double = f_to_double.getSize(); UInt old_nb_facet = element_to_facet.getSize() - nb_facet_to_double; UInt new_nb_facet = element_to_facet.getSize(); UInt old_nb_nodes = position.getSize(); UInt new_nb_nodes = old_nb_nodes + nb_facet_to_double; position.resize(new_nb_nodes); conn_facet.resize(new_nb_facet); UInt old_nb_doubled_nodes = doubled_nodes.getSize(); doubled_nodes.resize(old_nb_doubled_nodes + nb_facet_to_double); for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { UInt old_facet = f_to_double(facet); UInt new_facet = old_nb_facet + facet; ElementType type = element_to_facet(new_facet)[0].type; UInt el = element_to_facet(new_facet)[0].element; GhostType gt = element_to_facet(new_facet)[0].ghost_type; UInt old_node = conn_facet(old_facet); UInt new_node = old_nb_nodes + facet; /// update position position(new_node) = position(old_node); conn_facet(new_facet) = new_node; Array & conn_segment = mesh.getConnectivity(type, gt); UInt nb_nodes_per_segment = conn_segment.getNbComponent(); /// update facet connectivity UInt i; for (i = 0; conn_segment(el, i) != old_node && i <= nb_nodes_per_segment; ++i); conn_segment(el, i) = new_node; doubled_nodes(old_nb_doubled_nodes + facet, 0) = old_node; doubled_nodes(old_nb_doubled_nodes + facet, 1) = new_node; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MeshUtils::updateQuadraticSegments(Mesh & mesh, Mesh & mesh_facets, ElementType type_facet, GhostType gt_facet, Array & doubled_nodes) { AKANTU_DEBUG_IN(); if (type_facet != _segment_3) return; Array & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); UInt nb_facet_to_double = f_to_double.getSize(); UInt old_nb_facet = mesh_facets.getNbElement(type_facet, gt_facet) - nb_facet_to_double; Array & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); Array< std::vector > & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); /// this ones matter only for segments in 3D Array< std::vector > * el_to_subfacet_double = NULL; Array< std::vector > * f_to_subfacet_double = NULL; if (third_dim_segments) { el_to_subfacet_double = & mesh_facets.getData >("elements_to_subfacet_double", type_facet, gt_facet); f_to_subfacet_double = & mesh_facets.getData >("facets_to_subfacet_double", type_facet, gt_facet); } std::vector middle_nodes; for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { UInt old_facet = f_to_double(facet); UInt node = conn_facet(old_facet, 2); if (!mesh.isPureGhostNode(node)) middle_nodes.push_back(node); } UInt n = doubled_nodes.getSize(); doubleNodes(mesh, middle_nodes, doubled_nodes); for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { UInt old_facet = f_to_double(facet); UInt old_node = conn_facet(old_facet, 2); if (mesh.isPureGhostNode(old_node)) continue; UInt new_node = doubled_nodes(n, 1); UInt new_facet = old_nb_facet + facet; conn_facet(new_facet, 2) = new_node; if (third_dim_segments) { updateElementalConnectivity(mesh_facets, old_node, new_node, element_to_facet(new_facet)); updateElementalConnectivity(mesh, old_node, new_node, (*el_to_subfacet_double)(facet), &(*f_to_subfacet_double)(facet)); } else { updateElementalConnectivity(mesh, old_node, new_node, element_to_facet(new_facet)); } ++n; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::updateSubfacetToFacet(Mesh & mesh_facets, ElementType type_subfacet, GhostType gt_subfacet, bool facet_mode) { AKANTU_DEBUG_IN(); Array & sf_to_double = mesh_facets.getData("facet_to_double", type_subfacet, gt_subfacet); UInt nb_subfacet_to_double = sf_to_double.getSize(); /// update subfacet_to_facet vector ElementType type_facet = _not_defined; GhostType gt_facet = _casper; Array * subfacet_to_facet = NULL; UInt nb_subfacet_per_facet = 0; UInt old_nb_subfacet = mesh_facets.getNbElement(type_subfacet, gt_subfacet) - nb_subfacet_to_double; Array > * facet_list = NULL; if (facet_mode) facet_list = & mesh_facets.getData >("facets_to_subfacet_double", type_subfacet, gt_subfacet); else facet_list = & mesh_facets.getData >("subfacets_to_subsubfacet_double", type_subfacet, gt_subfacet); Element old_subfacet_el(type_subfacet, 0, gt_subfacet); Element new_subfacet_el(type_subfacet, 0, gt_subfacet); for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) { old_subfacet_el.element = sf_to_double(sf); new_subfacet_el.element = old_nb_subfacet + sf; for (UInt f = 0; f < (*facet_list)(sf).size(); ++f) { Element & facet = (*facet_list)(sf)[f]; if (facet.type != type_facet || facet.ghost_type != gt_facet) { type_facet = facet.type; gt_facet = facet.ghost_type; subfacet_to_facet = & mesh_facets.getSubelementToElement(type_facet, gt_facet); nb_subfacet_per_facet = subfacet_to_facet->getNbComponent(); } Element * sf_update = std::find(subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet, subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet + nb_subfacet_per_facet, old_subfacet_el); AKANTU_DEBUG_ASSERT(subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet != subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet + nb_subfacet_per_facet, "Subfacet not found"); *sf_update = new_subfacet_el; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::updateFacetToSubfacet(Mesh & mesh_facets, ElementType type_subfacet, GhostType gt_subfacet, bool facet_mode) { AKANTU_DEBUG_IN(); Array & sf_to_double = mesh_facets.getData("facet_to_double", type_subfacet, gt_subfacet); UInt nb_subfacet_to_double = sf_to_double.getSize(); Array< std::vector > & facet_to_subfacet = mesh_facets.getElementToSubelement(type_subfacet, gt_subfacet); Array< std::vector > * facet_to_subfacet_double = NULL; if (facet_mode) { facet_to_subfacet_double = & mesh_facets.getData >("facets_to_subfacet_double", type_subfacet, gt_subfacet); } else { facet_to_subfacet_double = & mesh_facets.getData >("subfacets_to_subsubfacet_double", type_subfacet, gt_subfacet); } UInt old_nb_subfacet = facet_to_subfacet.getSize(); facet_to_subfacet.resize(old_nb_subfacet + nb_subfacet_to_double); for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) facet_to_subfacet(old_nb_subfacet + sf) = (*facet_to_subfacet_double)(sf); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MeshUtils::doubleSubfacet(Mesh & mesh, Mesh & mesh_facets, Array & doubled_nodes) { AKANTU_DEBUG_IN(); if (spatial_dimension == 1) return; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType gt_subfacet = *gt; Mesh::type_iterator it = mesh_facets.firstType(0, gt_subfacet); Mesh::type_iterator end = mesh_facets.lastType(0, gt_subfacet); for(; it != end; ++it) { ElementType type_subfacet = *it; Array & sf_to_double = mesh_facets.getData("facet_to_double", type_subfacet, gt_subfacet); UInt nb_subfacet_to_double = sf_to_double.getSize(); if (nb_subfacet_to_double == 0) continue; AKANTU_DEBUG_ASSERT(type_subfacet == _point_1, "Only _point_1 subfacet doubling is supported at the moment"); Array & conn_subfacet = mesh_facets.getConnectivity(type_subfacet, gt_subfacet); UInt old_nb_subfacet = conn_subfacet.getSize(); UInt new_nb_subfacet = old_nb_subfacet + nb_subfacet_to_double; conn_subfacet.resize(new_nb_subfacet); std::vector nodes_to_double; UInt old_nb_doubled_nodes = doubled_nodes.getSize(); /// double nodes for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) { UInt old_subfacet = sf_to_double(sf); nodes_to_double.push_back(conn_subfacet(old_subfacet)); } doubleNodes(mesh, nodes_to_double, doubled_nodes); /// add new nodes in connectivity for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) { UInt new_subfacet = old_nb_subfacet + sf; UInt new_node = doubled_nodes(old_nb_doubled_nodes + sf, 1); conn_subfacet(new_subfacet) = new_node; } /// update facet and element connectivity Array > & f_to_subfacet_double = mesh_facets.getData >("facets_to_subfacet_double", type_subfacet, gt_subfacet); Array > & el_to_subfacet_double = mesh_facets.getData >("elements_to_subfacet_double", type_subfacet, gt_subfacet); Array > * sf_to_subfacet_double = NULL; if (spatial_dimension == 3) sf_to_subfacet_double = & mesh_facets.getData >("subfacets_to_subsubfacet_double", type_subfacet, gt_subfacet); for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) { UInt old_node = doubled_nodes(old_nb_doubled_nodes + sf, 0); UInt new_node = doubled_nodes(old_nb_doubled_nodes + sf, 1); updateElementalConnectivity(mesh, old_node, new_node, el_to_subfacet_double(sf), &f_to_subfacet_double(sf)); updateElementalConnectivity(mesh_facets, old_node, new_node, f_to_subfacet_double(sf)); if (spatial_dimension == 3) updateElementalConnectivity(mesh_facets, old_node, new_node, (*sf_to_subfacet_double)(sf)); } if (spatial_dimension == 2) { updateSubfacetToFacet(mesh_facets, type_subfacet, gt_subfacet, true); updateFacetToSubfacet(mesh_facets, type_subfacet, gt_subfacet, true); } else if (spatial_dimension == 3) { updateSubfacetToFacet(mesh_facets, type_subfacet, gt_subfacet, false); updateFacetToSubfacet(mesh_facets, type_subfacet, gt_subfacet, false); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::flipFacets(Mesh & mesh_facets, const ElementTypeMapArray & global_connectivity, GhostType gt_facet) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh_facets.getSpatialDimension(); /// get global connectivity for local mesh ElementTypeMapArray global_connectivity_tmp; mesh_facets.initElementTypeMapArray(global_connectivity_tmp, 1, spatial_dimension - 1, gt_facet, true, _ek_regular, true); mesh_facets.getGlobalConnectivity(global_connectivity_tmp, spatial_dimension - 1, gt_facet); Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, gt_facet); Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1, gt_facet); /// loop on every facet for(; it != end; ++it) { ElementType type_facet = *it; Array & connectivity = mesh_facets.getConnectivity(type_facet, gt_facet); const Array & g_connectivity = global_connectivity(type_facet, gt_facet); Array > & el_to_f = mesh_facets.getElementToSubelement(type_facet, gt_facet); Array & subfacet_to_facet = mesh_facets.getSubelementToElement(type_facet, gt_facet); UInt nb_subfacet_per_facet = subfacet_to_facet.getNbComponent(); UInt nb_nodes_per_facet = connectivity.getNbComponent(); UInt nb_facet = connectivity.getSize(); UInt nb_nodes_per_P1_facet = Mesh::getNbNodesPerElement(Mesh::getP1ElementType(type_facet)); Array & global_conn_tmp = global_connectivity_tmp(type_facet, gt_facet); Array::iterator > conn_it = connectivity.begin(nb_nodes_per_facet); Array::iterator > gconn_tmp_it = global_conn_tmp.begin(nb_nodes_per_facet); Array::const_iterator > conn_glob_it = g_connectivity.begin(nb_nodes_per_facet); Array::iterator > subf_to_f = subfacet_to_facet.begin(nb_subfacet_per_facet); UInt * conn_tmp_pointer = new UInt[nb_nodes_per_facet]; Vector conn_tmp(conn_tmp_pointer, nb_nodes_per_facet); Element el_tmp; Element * subf_tmp_pointer = new Element[nb_subfacet_per_facet]; Vector subf_tmp(subf_tmp_pointer, nb_subfacet_per_facet); for (UInt f = 0; f < nb_facet; ++f, ++conn_it, ++gconn_tmp_it, ++subf_to_f, ++conn_glob_it) { Vector & gconn_tmp = *gconn_tmp_it; const Vector & conn_glob = *conn_glob_it; Vector & conn_local = *conn_it; /// skip facet if connectivities are the same if (gconn_tmp == conn_glob) continue; /// re-arrange connectivity conn_tmp = conn_local; UInt * begin = conn_glob.storage(); UInt * end = conn_glob.storage() + nb_nodes_per_facet; for (UInt n = 0; n < nb_nodes_per_facet; ++n) { UInt * new_node = std::find(begin, end, gconn_tmp(n)); AKANTU_DEBUG_ASSERT(new_node != end, "Node not found"); UInt new_position = new_node - begin; conn_local(new_position) = conn_tmp(n); } /// if 3D, check if facets are just rotated if (spatial_dimension == 3) { /// find first node UInt * new_node = std::find(begin, end, gconn_tmp(0)); AKANTU_DEBUG_ASSERT(new_node != end, "Node not found"); UInt new_position = new_node - begin; UInt n = 1; /// count how many nodes in the received connectivity follow /// the same order of those in the local connectivity for (; n < nb_nodes_per_P1_facet && gconn_tmp(n) == conn_glob((new_position + n) % nb_nodes_per_P1_facet); ++n); /// skip the facet inversion if facet is just rotated if (n == nb_nodes_per_P1_facet) continue; } /// update data to invert facet el_tmp = el_to_f(f)[0]; el_to_f(f)[0] = el_to_f(f)[1]; el_to_f(f)[1] = el_tmp; subf_tmp = (*subf_to_f); (*subf_to_f)(0) = subf_tmp(1); (*subf_to_f)(1) = subf_tmp(0); } delete [] subf_tmp_pointer; delete [] conn_tmp_pointer; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::fillElementToSubElementsData(Mesh & mesh) { AKANTU_DEBUG_IN(); if(mesh.getNbElement(mesh.getSpatialDimension() - 1) == 0) { AKANTU_DEBUG_INFO("There are not facets, add them in the mesh file or call the buildFacet method."); return; } UInt spatial_dimension = mesh.getSpatialDimension(); ElementTypeMapArray barycenters("barycenter_tmp", mesh.getID()); mesh.initElementTypeMapArray(barycenters, spatial_dimension, _all_dimensions); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { Mesh::type_iterator it = mesh.firstType(_all_dimensions, *gt); Mesh::type_iterator end = mesh.lastType(_all_dimensions, *gt); for(; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it, *gt); Array & barycenters_arr = barycenters(*it, *gt); barycenters_arr.resize(nb_element); Array::vector_iterator bary = barycenters_arr.begin(spatial_dimension); Array::vector_iterator bary_end = barycenters_arr.end(spatial_dimension); for (UInt el = 0; bary != bary_end; ++bary, ++el) { mesh.getBarycenter(el, *it, bary->storage(), *gt); } } } for(Int sp(spatial_dimension); sp >= 1; --sp) { if(mesh.getNbElement(sp) == 0) continue; for (ghost_type_t::iterator git = ghost_type_t::begin(); git != ghost_type_t::end(); ++git) { Mesh::type_iterator tit = mesh.firstType(sp, *git); Mesh::type_iterator tend = mesh.lastType(sp, *git); for (;tit != tend; ++tit) { mesh.getSubelementToElementPointer(*tit, *git)->resize(mesh.getNbElement(*tit, *git)); mesh.getSubelementToElement(*tit, *git).clear(); } tit = mesh.firstType(sp - 1, *git); tend = mesh.lastType(sp - 1, *git); for (;tit != tend; ++tit) { mesh.getElementToSubelementPointer(*tit, *git)->resize(mesh.getNbElement(*tit, *git)); mesh.getElementToSubelement(*tit, *git).clear(); } } CSR nodes_to_elements; buildNode2Elements(mesh, nodes_to_elements, sp); Element facet_element; for (ghost_type_t::iterator git = ghost_type_t::begin(); git != ghost_type_t::end(); ++git) { Mesh::type_iterator tit = mesh.firstType(sp - 1, *git); Mesh::type_iterator tend = mesh.lastType(sp - 1, *git); facet_element.ghost_type = *git; for (;tit != tend; ++tit) { facet_element.type = *tit; Array< std::vector > & element_to_subelement = mesh.getElementToSubelement(*tit, *git); const Array & connectivity = mesh.getConnectivity(*tit, *git); Array::const_iterator< Vector > fit = connectivity.begin(mesh.getNbNodesPerElement(*tit)); Array::const_iterator< Vector > fend = connectivity.end(mesh.getNbNodesPerElement(*tit)); UInt fid = 0; for (;fit != fend; ++fit, ++fid) { const Vector & facet = *fit; facet_element.element = fid; std::map element_seen_counter; UInt nb_nodes_per_facet = mesh.getNbNodesPerElement(Mesh::getP1ElementType(*tit)); for (UInt n(0); n < nb_nodes_per_facet; ++n) { CSR::iterator eit = nodes_to_elements.begin(facet(n)); CSR::iterator eend = nodes_to_elements.end(facet(n)); for(;eit != eend; ++eit) { Element & elem = *eit; std::map::iterator cit = element_seen_counter.find(elem); if(cit != element_seen_counter.end()) { cit->second++; } else { element_seen_counter[elem] = 1; } } } std::vector connected_elements; std::map::iterator cit = element_seen_counter.begin(); std::map::iterator cend = element_seen_counter.end(); for(;cit != cend; ++cit) { if(cit->second == nb_nodes_per_facet) connected_elements.push_back(cit->first); } std::vector::iterator ceit = connected_elements.begin(); std::vector::iterator ceend = connected_elements.end(); for(;ceit != ceend; ++ceit) element_to_subelement(fid).push_back(*ceit); for (UInt ce = 0; ce < connected_elements.size(); ++ce) { Element & elem = connected_elements[ce]; Array & subelement_to_element = *(mesh.getSubelementToElementPointer(elem.type, elem.ghost_type)); UInt f(0); for(; f < mesh.getNbFacetsPerElement(elem.type) && subelement_to_element(elem.element, f) != ElementNull; ++f); AKANTU_DEBUG_ASSERT(f < mesh.getNbFacetsPerElement(elem.type), "The element " << elem << " seems to have too many facets!! (" << f << " < " << mesh.getNbFacetsPerElement(elem.type) << ")"); subelement_to_element(elem.element, f) = facet_element; } } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template bool MeshUtils::findElementsAroundSubfacet(const Mesh & mesh, const Mesh & mesh_facets, const Element & starting_element, const Element & end_facet, const Vector & subfacet_connectivity, std::vector & elem_list, std::vector & facet_list, std::vector * subfacet_list) { AKANTU_DEBUG_IN(); /// preallocated stuff before starting bool facet_matched = false; elem_list.clear(); facet_list.clear(); if (third_dim_points) subfacet_list->clear(); elem_list.push_back(starting_element); const Array * facet_connectivity = NULL; const Array * sf_connectivity = NULL; const Array * facet_to_element = NULL; const Array * subfacet_to_facet = NULL; ElementType current_type = _not_defined; GhostType current_ghost_type = _casper; ElementType current_facet_type = _not_defined; GhostType current_facet_ghost_type = _casper; ElementType current_subfacet_type = _not_defined; GhostType current_subfacet_ghost_type = _casper; const Array< std::vector > * element_to_facet = NULL; const Element * opposing_el = NULL; std::queue elements_to_check; elements_to_check.push(starting_element); /// keep going until there are elements to check while (!elements_to_check.empty()) { /// check current element Element & current_el = elements_to_check.front(); if (current_el.type != current_type || current_el.ghost_type != current_ghost_type) { current_type = current_el.type; current_ghost_type = current_el.ghost_type; facet_to_element = & mesh_facets.getSubelementToElement(current_type, current_ghost_type); } /// loop over each facet of the element for (UInt f = 0; f < facet_to_element->getNbComponent(); ++f) { const Element & current_facet = (*facet_to_element)(current_el.element, f); if (current_facet == ElementNull) continue; if (current_facet_type != current_facet.type || current_facet_ghost_type != current_facet.ghost_type) { current_facet_type = current_facet.type; current_facet_ghost_type = current_facet.ghost_type; element_to_facet = & mesh_facets.getElementToSubelement(current_facet_type, current_facet_ghost_type); facet_connectivity = & mesh_facets.getConnectivity(current_facet_type, current_facet_ghost_type); if (third_dim_points) subfacet_to_facet = & mesh_facets.getSubelementToElement(current_facet_type, current_facet_ghost_type); } /// check if end facet is reached if (current_facet == end_facet) facet_matched = true; /// add this facet if not already passed if (std::find(facet_list.begin(), facet_list.end(), current_facet) == facet_list.end() && hasElement(*facet_connectivity, current_facet, subfacet_connectivity)) { facet_list.push_back(current_facet); if (third_dim_points) { /// check subfacets for (UInt sf = 0; sf < subfacet_to_facet->getNbComponent(); ++sf) { const Element & current_subfacet = (*subfacet_to_facet)(current_facet.element, sf); if (current_subfacet == ElementNull) continue; if (current_subfacet_type != current_subfacet.type || current_subfacet_ghost_type != current_subfacet.ghost_type) { current_subfacet_type = current_subfacet.type; current_subfacet_ghost_type = current_subfacet.ghost_type; sf_connectivity = & mesh_facets.getConnectivity(current_subfacet_type, current_subfacet_ghost_type); } if (std::find(subfacet_list->begin(), subfacet_list->end(), current_subfacet) == subfacet_list->end() && hasElement(*sf_connectivity, current_subfacet, subfacet_connectivity)) subfacet_list->push_back(current_subfacet); } } } else continue; /// consider opposing element if ( (*element_to_facet)(current_facet.element)[0] == current_el) opposing_el = & (*element_to_facet)(current_facet.element)[1]; else opposing_el = & (*element_to_facet)(current_facet.element)[0]; /// skip null elements since they are on a boundary if (*opposing_el == ElementNull) continue; /// skip this element if already added if ( std::find(elem_list.begin(), elem_list.end(), *opposing_el) != elem_list.end() ) continue; /// only regular elements have to be checked if (opposing_el->kind == _ek_regular) elements_to_check.push(*opposing_el); elem_list.push_back(*opposing_el); #ifndef AKANTU_NDEBUG const Array & conn_elem = mesh.getConnectivity(opposing_el->type, opposing_el->ghost_type); AKANTU_DEBUG_ASSERT(hasElement(conn_elem, *opposing_el, subfacet_connectivity), "Subfacet doesn't belong to this element"); #endif } /// erased checked element from the list elements_to_check.pop(); } AKANTU_DEBUG_OUT(); return facet_matched; } /* -------------------------------------------------------------------------- */ void MeshUtils::buildSegmentToNodeType(const Mesh & mesh, Mesh & mesh_facets, DistributedSynchronizer * synchronizer) { buildAllFacets(mesh, mesh_facets, 1, synchronizer); UInt spatial_dimension = mesh.getSpatialDimension(); const ElementTypeMapArray & element_to_rank = synchronizer->getPrankToElement(); Int local_rank = StaticCommunicator::getStaticCommunicator().whoAmI(); 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(1, ghost_type); Mesh::type_iterator end = mesh_facets.lastType(1, ghost_type); for(; it != end; ++it) { ElementType type = *it; UInt nb_segments = mesh_facets.getNbElement(type, ghost_type); // allocate the data Array & segment_to_nodetype = *(mesh_facets.getDataPointer("segment_to_nodetype", type, ghost_type)); std::set connected_elements; const Array< std::vector > & segment_to_2Delement = mesh_facets.getElementToSubelement(type, ghost_type); // loop over segments for (UInt s = 0; s < nb_segments; ++s) { // determine the elements connected to the segment connected_elements.clear(); const std::vector & twoD_elements = segment_to_2Delement(s); if (spatial_dimension == 2) { // if 2D just take the elements connected to the segments connected_elements.insert(twoD_elements.begin(), twoD_elements.end()); } else if (spatial_dimension == 3) { // if 3D a second loop is needed to get to the 3D elements std::vector::const_iterator facet = twoD_elements.begin(); for (; facet != twoD_elements.end(); ++facet) { const std::vector & threeD_elements = mesh_facets.getElementToSubelement(facet->type, facet->ghost_type)(facet->element); connected_elements.insert(threeD_elements.begin(), threeD_elements.end()); } } // get the minimum processor rank associated to the connected // elements and verify if ghost and not ghost elements are // found Int minimum_rank = std::numeric_limits::max(); // two booleans saying if not ghost and ghost elements are found in the loop bool ghost_found[2]; ghost_found[0] = false; ghost_found[1] = false; std::set::iterator connected_elements_it = connected_elements.begin(); for (; connected_elements_it != connected_elements.end(); ++connected_elements_it) { if (*connected_elements_it == ElementNull) continue; ghost_found[connected_elements_it->ghost_type] = true; const Array & el_to_rank_array = element_to_rank(connected_elements_it->type, connected_elements_it->ghost_type); minimum_rank = std::min(minimum_rank, Int(el_to_rank_array(connected_elements_it->element))); } // if no ghost elements are found the segment is local if (!ghost_found[1]) segment_to_nodetype(s) = -1; // if no not ghost elements are found the segment is pure ghost else if (!ghost_found[0]) segment_to_nodetype(s) = -3; // if the minimum rank is equal to the local rank, the segment is master else if (local_rank == minimum_rank) segment_to_nodetype(s) = -2; // if the minimum rank is less than the local rank, the segment is slave else if (local_rank > minimum_rank) segment_to_nodetype(s) = minimum_rank; else AKANTU_DEBUG_ERROR("The local rank cannot be smaller than the minimum rank if both ghost and not ghost elements are found"); } } } } /* -------------------------------------------------------------------------- */ -UInt MeshUtils::updateLocalMasterGlobalConnectivity(Mesh & mesh, UInt old_nb_nodes) { - UInt local_nb_new_nodes = mesh.getNbNodes() - old_nb_nodes; - +UInt MeshUtils::updateLocalMasterGlobalConnectivity(Mesh & mesh, UInt local_nb_new_nodes) { StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Int rank = comm.whoAmI(); Int nb_proc = comm.getNbProc(); if (nb_proc == 1) return local_nb_new_nodes; /// resize global ids array Array & nodes_global_ids = mesh.getGlobalNodesIds(); - UInt nb_old_nodes = nodes_global_ids.getSize(); + UInt old_nb_nodes = mesh.getNbNodes() - local_nb_new_nodes; - nodes_global_ids.resize(nb_old_nodes + local_nb_new_nodes); + nodes_global_ids.resize(mesh.getNbNodes()); /// compute amount of local or master doubled nodes Vector local_master_nodes(nb_proc); for (UInt n = old_nb_nodes; n < mesh.getNbNodes(); ++n) if (mesh.isLocalOrMasterNode(n)) ++local_master_nodes(rank); comm.allGather(local_master_nodes.storage(), 1); /// update global number of nodes UInt total_nb_new_nodes = std::accumulate(local_master_nodes.storage(), local_master_nodes.storage() + nb_proc, 0); if (total_nb_new_nodes == 0) return 0; /// set global ids of local and master nodes UInt starting_index = std::accumulate(local_master_nodes.storage(), local_master_nodes.storage() + rank, mesh.getNbGlobalNodes()); for (UInt n = old_nb_nodes; n < mesh.getNbNodes(); ++n) { if (mesh.isLocalOrMasterNode(n)) { nodes_global_ids(n) = starting_index; ++starting_index; } } mesh.nb_global_nodes += total_nb_new_nodes; return total_nb_new_nodes; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ // LocalWords: ElementType diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index 0baeed877..3185a012c 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1824 +1,1833 @@ /** * @file solid_mechanics_model.cc * * @author Guillaume Anciaux * @author David Simon Kammer * @author Aurelia Isabel Cuba Ramos * @author Daniel Pino Muñoz * @author Nicolas Richart * * @date creation: Tue Jul 27 2010 * @date last modification: Fri Sep 19 2014 * * @brief Implementation of the SolidMechanicsModel class * * @section LICENSE * * Copyright (©) 2010-2012, 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_math.hh" #include "aka_common.hh" #include "solid_mechanics_model.hh" #include "group_manager_inline_impl.cc" #include "dumpable_inline_impl.hh" #include "integration_scheme_2nd_order.hh" #include "element_group.hh" #include "static_communicator.hh" #include "dof_synchronizer.hh" #include "element_group.hh" #include #ifdef AKANTU_USE_MUMPS #include "solver_mumps.hh" #endif #ifdef AKANTU_USE_PETSC #include "solver_petsc.hh" #include "petsc_matrix.hh" #endif #ifdef AKANTU_USE_IOHELPER # include "dumper_field.hh" # include "dumper_paraview.hh" # include "dumper_homogenizing_field.hh" # include "dumper_internal_material_field.hh" # include "dumper_elemental_field.hh" # include "dumper_material_padders.hh" # include "dumper_element_partition.hh" # include "dumper_iohelper.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ const SolidMechanicsModelOptions default_solid_mechanics_model_options(_explicit_lumped_mass, false); /* -------------------------------------------------------------------------- */ /** * A solid mechanics model need a mesh and a dimension to be created. the model * by it self can not do a lot, the good init functions should be called in * order to configure the model depending on what we want to do. * * @param mesh mesh representing the model we want to simulate * @param dim spatial dimension of the problem, if dim = 0 (default value) the * dimension of the problem is assumed to be the on of the mesh * @param id an id to identify the model */ SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : Model(mesh, dim, id, memory_id), BoundaryCondition(), time_step(NAN), f_m2a(1.0), mass_matrix(NULL), velocity_damping_matrix(NULL), stiffness_matrix(NULL), jacobian_matrix(NULL), material_index("material index", id), material_local_numbering("material local numbering", id), material_selector(new DefaultMaterialSelector(material_index)), is_default_material_selector(true), integrator(NULL), increment_flag(false), solver(NULL), synch_parallel(NULL), are_materials_instantiated(false) { AKANTU_DEBUG_IN(); createSynchronizerRegistry(this); registerFEEngineObject("SolidMechanicsFEEngine", mesh, spatial_dimension); this->displacement = NULL; this->mass = NULL; this->velocity = NULL; this->acceleration = NULL; this->force = NULL; this->residual = NULL; this->blocked_dofs = NULL; this->increment = NULL; this->increment_acceleration = NULL; this->dof_synchronizer = NULL; this->previous_displacement = NULL; materials.clear(); mesh.registerEventHandler(*this); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("paraview_all", id, true); this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular); #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModel::~SolidMechanicsModel() { AKANTU_DEBUG_IN(); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { delete *mat_it; } materials.clear(); delete integrator; delete solver; delete mass_matrix; delete velocity_damping_matrix; if(stiffness_matrix && stiffness_matrix != jacobian_matrix) delete stiffness_matrix; delete jacobian_matrix; delete synch_parallel; if(is_default_material_selector) { delete material_selector; material_selector = NULL; } AKANTU_DEBUG_OUT(); } void SolidMechanicsModel::setTimeStep(Real time_step) { this->time_step = time_step; #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper().setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ /** * This function groups many of the initialization in on function. For most of * basics case the function should be enough. The functions initialize the * model, the internal vectors, set them to 0, and depending on the parameters * it also initialize the explicit or implicit solver. * * @param material_file the file containing the materials to use * @param method the analysis method wanted. See the akantu::AnalysisMethod for * the different possibilities */ void SolidMechanicsModel::initFull(const ModelOptions & options) { Model::initFull(options); const SolidMechanicsModelOptions & smm_options = dynamic_cast(options); method = smm_options.analysis_method; // initialize the vectors initArrays(); // set the initial condition to 0 force->clear(); velocity->clear(); acceleration->clear(); displacement->clear(); // initialize pbc if(pbc_pair.size()!=0) initPBC(); // initialize the time integration schemes switch(method) { case _explicit_lumped_mass: initExplicit(); break; case _explicit_consistent_mass: initSolver(); initExplicit(); break; case _implicit_dynamic: initImplicit(true); break; case _static: initImplicit(false); break; default: AKANTU_EXCEPTION("analysis method not recognised by SolidMechanicsModel"); break; } // initialize the materials if(this->parser->getLastParsedFile() != "") { instantiateMaterials(); } if(!smm_options.no_init_materials) { initMaterials(); } if(increment_flag) initBC(*this, *displacement, *increment, *force); else initBC(*this, *displacement, *force); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initParallel(MeshPartition * partition, DataAccessor * data_accessor) { AKANTU_DEBUG_IN(); if (data_accessor == NULL) data_accessor = this; synch_parallel = &createParallelSynch(partition,data_accessor); synch_registry->registerSynchronizer(*synch_parallel, _gst_material_id); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_mass); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_stress); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_boundary); synch_registry->registerSynchronizer(*synch_parallel, _gst_for_dump); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initFEEngineBoundary() { FEEngine & fem_boundary = getFEEngineBoundary(); fem_boundary.initShapeFunctions(_not_ghost); fem_boundary.initShapeFunctions(_ghost); fem_boundary.computeNormalsOnControlPoints(_not_ghost); fem_boundary.computeNormalsOnControlPoints(_ghost); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initExplicit(AnalysisMethod analysis_method) { AKANTU_DEBUG_IN(); //in case of switch from implicit to explicit if(!this->isExplicit()) method = analysis_method; if (integrator) delete integrator; integrator = new CentralDifference(); UInt nb_nodes = acceleration->getSize(); UInt nb_degree_of_freedom = acceleration->getNbComponent(); std::stringstream sstr; sstr << id << ":increment_acceleration"; increment_acceleration = &(alloc(sstr.str(), nb_nodes, nb_degree_of_freedom, Real())); AKANTU_DEBUG_OUT(); } void SolidMechanicsModel::initArraysPreviousDisplacment() { AKANTU_DEBUG_IN(); SolidMechanicsModel::setIncrementFlagOn(); UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_disp_t; sstr_disp_t << id << ":previous_displacement"; previous_displacement = &(alloc (sstr_disp_t.str(), nb_nodes, spatial_dimension, 0.)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Allocate all the needed vectors. By default their are not necessarily set to * 0 * */ void SolidMechanicsModel::initArrays() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_disp; sstr_disp << id << ":displacement"; // std::stringstream sstr_mass; sstr_mass << id << ":mass"; std::stringstream sstr_velo; sstr_velo << id << ":velocity"; std::stringstream sstr_acce; sstr_acce << id << ":acceleration"; std::stringstream sstr_forc; sstr_forc << id << ":force"; std::stringstream sstr_resi; sstr_resi << id << ":residual"; std::stringstream sstr_boun; sstr_boun << id << ":blocked_dofs"; displacement = &(alloc(sstr_disp.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); // mass = &(alloc(sstr_mass.str(), nb_nodes, spatial_dimension, 0)); velocity = &(alloc(sstr_velo.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); acceleration = &(alloc(sstr_acce.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); force = &(alloc(sstr_forc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); residual = &(alloc(sstr_resi.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); blocked_dofs = &(alloc(sstr_boun.str(), nb_nodes, spatial_dimension, false)); std::stringstream sstr_curp; sstr_curp << id << ":current_position"; current_position = &(alloc(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE)); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_not_defined); for(; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it, gt); material_index.alloc(nb_element, 1, *it, gt); material_local_numbering.alloc(nb_element, 1, *it, gt); } } dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian * */ void SolidMechanicsModel::initModel() { /// \todo add the current position as a parameter to initShapeFunctions for /// large deformation getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initPBC() { Model::initPBC(); registerPBCSynchronizer(); // as long as there are ones on the diagonal of the matrix, we can put boudandary true for slaves std::map::iterator it = pbc_pair.begin(); std::map::iterator end = pbc_pair.end(); UInt dim = mesh.getSpatialDimension(); while(it != end) { for (UInt i=0; iregisterSynchronizer(*synch, _gst_smm_uv); synch_registry->registerSynchronizer(*synch, _gst_smm_mass); synch_registry->registerSynchronizer(*synch, _gst_smm_res); synch_registry->registerSynchronizer(*synch, _gst_for_dump); changeLocalEquationNumberForPBC(pbc_pair, mesh.getSpatialDimension()); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateCurrentPosition() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); current_position->resize(nb_nodes); Real * current_position_val = current_position->storage(); Real * position_val = mesh.getNodes().storage(); Real * displacement_val = displacement->storage(); /// compute current_position = initial_position + displacement memcpy(current_position_val, position_val, nb_nodes*spatial_dimension*sizeof(Real)); for (UInt n = 0; n < nb_nodes*spatial_dimension; ++n) { *current_position_val++ += *displacement_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initializeUpdateResidualData() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); residual->resize(nb_nodes); /// copy the forces in residual for boundary conditions memcpy(residual->storage(), force->storage(), nb_nodes*spatial_dimension*sizeof(Real)); // start synchronization synch_registry->asynchronousSynchronize(_gst_smm_uv); synch_registry->waitEndSynchronize(_gst_smm_uv); updateCurrentPosition(); AKANTU_DEBUG_OUT(); } /*----------------------------------------------------------------------------*/ void SolidMechanicsModel::reInitialize() { } /* -------------------------------------------------------------------------- */ /* Explicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /** * This function compute the second member of the motion equation. That is to * say the sum of forces @f$ r = F_{ext} - F_{int} @f$. @f$ F_{ext} @f$ is given * by the user in the force vector, and @f$ F_{int} @f$ is computed as @f$ * F_{int} = \int_{\Omega} N \sigma d\Omega@f$ * */ void SolidMechanicsModel::updateResidual(bool need_initialize) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); // f = f_ext - f_int // f = f_ext if(need_initialize) initializeUpdateResidualData(); AKANTU_DEBUG_INFO("Compute local stresses"); std::vector::iterator mat_it; for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllStresses(_not_ghost); } #ifdef AKANTU_DAMAGE_NON_LOCAL /* ------------------------------------------------------------------------ */ /* Computation of the non local part */ synch_registry->asynchronousSynchronize(_gst_mnl_for_average); AKANTU_DEBUG_INFO("Compute non local stresses for local elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_not_ghost); } AKANTU_DEBUG_INFO("Wait distant non local stresses"); synch_registry->waitEndSynchronize(_gst_mnl_for_average); AKANTU_DEBUG_INFO("Compute non local stresses for ghosts elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_ghost); } #endif /* ------------------------------------------------------------------------ */ /* assembling the forces internal */ // communicate the stress AKANTU_DEBUG_INFO("Send data for residual assembly"); synch_registry->asynchronousSynchronize(_gst_smm_stress); AKANTU_DEBUG_INFO("Assemble residual for local elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.assembleResidual(_not_ghost); } AKANTU_DEBUG_INFO("Wait distant stresses"); // finalize communications synch_registry->waitEndSynchronize(_gst_smm_stress); AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.assembleResidual(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeStresses() { if (isExplicit()) { // start synchronization synch_registry->asynchronousSynchronize(_gst_smm_uv); synch_registry->waitEndSynchronize(_gst_smm_uv); // compute stresses on all local elements for each materials std::vector::iterator mat_it; for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllStresses(_not_ghost); } /* ------------------------------------------------------------------------ */ #ifdef AKANTU_DAMAGE_NON_LOCAL /* Computation of the non local part */ synch_registry->asynchronousSynchronize(_gst_mnl_for_average); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_not_ghost); } synch_registry->waitEndSynchronize(_gst_mnl_for_average); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_ghost); } #endif } else { std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllStressesFromTangentModuli(_not_ghost); } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateResidualInternal() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Update the residual"); // f = f_ext - f_int - Ma - Cv = r - Ma - Cv; if(method != _static) { // f -= Ma if(mass_matrix) { // if full mass_matrix Array * Ma = new Array(*acceleration, true, "Ma"); *Ma *= *mass_matrix; /// \todo check unit conversion for implicit dynamics // *Ma /= f_m2a *residual -= *Ma; delete Ma; } else if (mass) { // else lumped mass UInt nb_nodes = acceleration->getSize(); UInt nb_degree_of_freedom = acceleration->getNbComponent(); Real * mass_val = mass->storage(); Real * accel_val = acceleration->storage(); Real * res_val = residual->storage(); bool * blocked_dofs_val = blocked_dofs->storage(); for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) { if(!(*blocked_dofs_val)) { *res_val -= *accel_val * *mass_val /f_m2a; } blocked_dofs_val++; res_val++; mass_val++; accel_val++; } } else { AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix."); } // f -= Cv if(velocity_damping_matrix) { Array * Cv = new Array(*velocity); *Cv *= *velocity_damping_matrix; *residual -= *Cv; delete Cv; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateAcceleration() { AKANTU_DEBUG_IN(); updateResidualInternal(); if(method == _explicit_lumped_mass) { /* residual = residual_{n+1} - M * acceleration_n therefore solution = increment acceleration not acceleration */ solveLumped(*increment_acceleration, *mass, *residual, *blocked_dofs, f_m2a); } else if (method == _explicit_consistent_mass) { solve(*increment_acceleration); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveLumped(Array & x, const Array & A, const Array & b, const Array & blocked_dofs, Real alpha) { Real * A_val = A.storage(); Real * b_val = b.storage(); Real * x_val = x.storage(); bool * blocked_dofs_val = blocked_dofs.storage(); UInt nb_degrees_of_freedom = x.getSize() * x.getNbComponent(); for (UInt n = 0; n < nb_degrees_of_freedom; ++n) { if(!(*blocked_dofs_val)) { *x_val = alpha * (*b_val / *A_val); } x_val++; A_val++; b_val++; blocked_dofs_val++; } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitPred() { AKANTU_DEBUG_IN(); if(increment_flag) { if(previous_displacement) increment->copy(*previous_displacement); else increment->copy(*displacement); } AKANTU_DEBUG_ASSERT(integrator,"itegrator should have been allocated: " << "have called initExplicit ? " << "or initImplicit ?"); integrator->integrationSchemePred(time_step, *displacement, *velocity, *acceleration, *blocked_dofs); if(increment_flag) { Real * inc_val = increment->storage(); Real * dis_val = displacement->storage(); UInt nb_degree_of_freedom = displacement->getSize() * displacement->getNbComponent(); for (UInt n = 0; n < nb_degree_of_freedom; ++n) { *inc_val = *dis_val - *inc_val; inc_val++; dis_val++; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitCorr() { AKANTU_DEBUG_IN(); integrator->integrationSchemeCorrAccel(time_step, *displacement, *velocity, *acceleration, *blocked_dofs, *increment_acceleration); if(previous_displacement) previous_displacement->copy(*displacement); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveStep() { AKANTU_DEBUG_IN(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeSolveStepEvent(method)); this->explicitPred(); this->updateResidual(); this->updateAcceleration(); this->explicitCorr(); EventManager::sendEvent(SolidMechanicsModelEvent::AfterSolveStepEvent(method)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Implicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /** * Initialize the solver and create the sparse matrices needed. * */ void SolidMechanicsModel::initSolver(__attribute__((unused)) SolverOptions & options) { #if !defined(AKANTU_USE_MUMPS) && !defined(AKANTU_USE_PETSC)// or other solver in the future \todo add AKANTU_HAS_SOLVER in CMake AKANTU_DEBUG_ERROR("You should at least activate one solver."); #else UInt nb_global_nodes = mesh.getNbGlobalNodes(); delete jacobian_matrix; std::stringstream sstr; sstr << id << ":jacobian_matrix"; #ifdef AKANTU_USE_PETSC jacobian_matrix = new PETScMatrix(nb_global_nodes * spatial_dimension, _symmetric, sstr.str(), memory_id); #else jacobian_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, _symmetric, sstr.str(), memory_id); #endif //AKANTU_USE PETSC jacobian_matrix->buildProfile(mesh, *dof_synchronizer, spatial_dimension); if (!isExplicit()) { delete stiffness_matrix; std::stringstream sstr_sti; sstr_sti << id << ":stiffness_matrix"; #ifdef AKANTU_USE_PETSC stiffness_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, _symmetric, sstr.str(), memory_id); stiffness_matrix->buildProfile(mesh, *dof_synchronizer, spatial_dimension); #else stiffness_matrix = new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id); #endif //AKANTU_USE_PETSC } delete solver; std::stringstream sstr_solv; sstr_solv << id << ":solver"; #ifdef AKANTU_USE_PETSC solver = new SolverPETSc(*jacobian_matrix, sstr_solv.str()); #elif defined(AKANTU_USE_MUMPS) solver = new SolverMumps(*jacobian_matrix, sstr_solv.str()); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS if(solver) solver->initialize(options); #endif //AKANTU_HAS_SOLVER } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initJacobianMatrix() { #if defined(AKANTU_USE_MUMPS) && !defined(AKANTU_USE_PETSC) // @todo make it more flexible: this is an ugly patch to treat the case of non // fix profile of the K matrix delete jacobian_matrix; std::stringstream sstr_sti; sstr_sti << id << ":jacobian_matrix"; jacobian_matrix = new SparseMatrix(*stiffness_matrix, sstr_sti.str(), memory_id); std::stringstream sstr_solv; sstr_solv << id << ":solver"; delete solver; solver = new SolverMumps(*jacobian_matrix, sstr_solv.str()); if(solver) solver->initialize(_solver_no_options); #else AKANTU_DEBUG_ERROR("You need to activate the solver mumps."); #endif } /* -------------------------------------------------------------------------- */ /** * Initialize the implicit solver, either for dynamic or static cases, * * @param dynamic */ void SolidMechanicsModel::initImplicit(bool dynamic, SolverOptions & solver_options) { AKANTU_DEBUG_IN(); method = dynamic ? _implicit_dynamic : _static; if (!increment) setIncrementFlagOn(); initSolver(solver_options); if(method == _implicit_dynamic) { if(integrator) delete integrator; integrator = new TrapezoidalRule2(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initialAcceleration() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving Ma = f"); Solver * acc_solver = NULL; std::stringstream sstr; sstr << id << ":tmp_mass_matrix"; SparseMatrix * tmp_mass = new SparseMatrix(*mass_matrix, sstr.str(), memory_id); #ifdef AKANTU_USE_MUMPS std::stringstream sstr_solver; sstr << id << ":solver_mass_matrix"; acc_solver = new SolverMumps(*mass_matrix, sstr_solver.str()); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS acc_solver->initialize(); tmp_mass->applyBoundary(*blocked_dofs); acc_solver->setRHS(*residual); acc_solver->solve(*acceleration); delete acc_solver; delete tmp_mass; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix."); stiffness_matrix->clear(); // call compute stiffness matrix on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->assembleStiffnessMatrix(_not_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SparseMatrix & SolidMechanicsModel::initVelocityDampingMatrix() { if(!velocity_damping_matrix) velocity_damping_matrix = new SparseMatrix(*jacobian_matrix, id + ":velocity_damping_matrix", memory_id); return *velocity_damping_matrix; } /* -------------------------------------------------------------------------- */ template<> bool SolidMechanicsModel::testConvergence<_scc_increment>(Real tolerance, Real & error){ AKANTU_DEBUG_IN(); UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent(); error = 0; Real norm[2] = {0., 0.}; Real * increment_val = increment->storage(); bool * blocked_dofs_val = blocked_dofs->storage(); Real * displacement_val = displacement->storage(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); for (UInt d = 0; d < nb_degree_of_freedom; ++d) { if(!(*blocked_dofs_val) && is_local_node) { norm[0] += *increment_val * *increment_val; norm[1] += *displacement_val * *displacement_val; } blocked_dofs_val++; increment_val++; displacement_val++; } } StaticCommunicator::getStaticCommunicator().allReduce(norm, 2, _so_sum); norm[0] = sqrt(norm[0]); norm[1] = sqrt(norm[1]); AKANTU_DEBUG_ASSERT(!Math::isnan(norm[0]), "Something goes wrong in the solve phase"); if (norm[1] < Math::getTolerance()) { error = norm[0]; AKANTU_DEBUG_OUT(); // cout<<"Error 1: "< Math::getTolerance()) error = norm[0] / norm[1]; else error = norm[0]; //In case the total displacement is zero! // cout<<"Error 2: "< bool SolidMechanicsModel::testConvergence<_scc_residual>(Real tolerance, Real & norm) { AKANTU_DEBUG_IN(); UInt nb_nodes = residual->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent(); norm = 0; Real * residual_val = residual->storage(); bool * blocked_dofs_val = blocked_dofs->storage(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); if(is_local_node) { for (UInt d = 0; d < nb_degree_of_freedom; ++d) { if(!(*blocked_dofs_val)) { norm += *residual_val * *residual_val; } blocked_dofs_val++; residual_val++; } } else { blocked_dofs_val += spatial_dimension; residual_val += spatial_dimension; } } StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum); norm = sqrt(norm); AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); return (norm < tolerance); } /* -------------------------------------------------------------------------- */ template<> bool SolidMechanicsModel::testConvergence<_scc_residual_mass_wgh>(Real tolerance, Real & norm) { AKANTU_DEBUG_IN(); UInt nb_nodes = residual->getSize(); norm = 0; Real * residual_val = residual->storage(); Real * mass_val = this->mass->storage(); bool * blocked_dofs_val = blocked_dofs->storage(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); if(is_local_node) { for (UInt d = 0; d < spatial_dimension; ++d) { if(!(*blocked_dofs_val)) { norm += *residual_val * *residual_val/(*mass_val * *mass_val); } blocked_dofs_val++; residual_val++; mass_val++; } } else { blocked_dofs_val += spatial_dimension; residual_val += spatial_dimension; mass_val += spatial_dimension; } } StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum); norm = sqrt(norm); AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); return (norm < tolerance); } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceResidual(Real tolerance){ AKANTU_DEBUG_IN(); Real error=0; bool res = this->testConvergence<_scc_residual>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceResidual(Real tolerance, Real & error){ AKANTU_DEBUG_IN(); bool res = this->testConvergence<_scc_residual>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance){ AKANTU_DEBUG_IN(); Real error=0; bool res = this->testConvergence<_scc_increment>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance, Real & error){ AKANTU_DEBUG_IN(); bool res = this->testConvergence<_scc_increment>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::implicitPred() { AKANTU_DEBUG_IN(); if(previous_displacement) previous_displacement->copy(*displacement); if(method == _implicit_dynamic) integrator->integrationSchemePred(time_step, *displacement, *velocity, *acceleration, *blocked_dofs); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::implicitCorr() { AKANTU_DEBUG_IN(); if(method == _implicit_dynamic) { integrator->integrationSchemeCorrDispl(time_step, *displacement, *velocity, *acceleration, *blocked_dofs, *increment); } else { UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes; Real * incr_val = increment->storage(); Real * disp_val = displacement->storage(); bool * boun_val = blocked_dofs->storage(); for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++disp_val, ++incr_val, ++boun_val){ *incr_val *= (1. - *boun_val); *disp_val += *incr_val; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateIncrement() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(previous_displacement,"The previous displacement has to be initialized." << " Are you working with Finite or Ineslactic deformations?"); UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes; Real * incr_val = increment->storage(); Real * disp_val = displacement->storage(); Real * prev_disp_val = previous_displacement->storage(); for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++disp_val, ++incr_val, ++prev_disp_val) *incr_val = (*disp_val - *prev_disp_val); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updatePreviousDisplacement() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(previous_displacement,"The previous displacement has to be initialized." << " Are you working with Finite or Ineslactic deformations?"); previous_displacement->copy(*displacement); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Information */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::synchronizeBoundaries() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(synch_registry,"Synchronizer registry was not initialized." << " Did you call initParallel?"); synch_registry->synchronize(_gst_smm_boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::synchronizeResidual() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(synch_registry,"Synchronizer registry was not initialized." << " Did you call initPBC?"); synch_registry->synchronize(_gst_smm_res); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setIncrementFlagOn() { AKANTU_DEBUG_IN(); if(!increment) { UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_inc; sstr_inc << id << ":increment"; increment = &(alloc(sstr_inc.str(), nb_nodes, spatial_dimension, 0.)); } increment_flag = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Real min_dt = getStableTimeStep(_not_ghost); /// reduction min over all processors StaticCommunicator::getStaticCommunicator().allReduce(&min_dt, 1, _so_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); Material ** mat_val = &(materials.at(0)); Real min_dt = std::numeric_limits::max(); updateCurrentPosition(); Element elem; elem.ghost_type = ghost_type; elem.kind = _ek_regular; Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator end = mesh.lastType(spatial_dimension, ghost_type); for(; it != end; ++it) { elem.type = *it; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*it); UInt nb_element = mesh.getNbElement(*it); Array::const_scalar_iterator mat_indexes = material_index(*it, ghost_type).begin(); Array::const_scalar_iterator mat_loc_num = material_local_numbering(*it, ghost_type).begin(); Array X(0, nb_nodes_per_element*spatial_dimension); FEEngine::extractNodalToElementField(mesh, *current_position, X, *it, _not_ghost); Array::matrix_iterator X_el = X.begin(spatial_dimension, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++X_el, ++mat_indexes, ++mat_loc_num) { elem.element = *mat_loc_num; Real el_h = getFEEngine().getElementInradius(*X_el, *it); Real el_c = mat_val[*mat_indexes]->getCelerity(elem); Real el_dt = el_h / el_c; min_dt = std::min(min_dt, el_dt); } } AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy() { AKANTU_DEBUG_IN(); if (!mass && !mass_matrix) AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix."); Real ekin = 0.; UInt nb_nodes = mesh.getNbNodes(); Real * vel_val = velocity->storage(); Real * mass_val = mass->storage(); for (UInt n = 0; n < nb_nodes; ++n) { Real mv2 = 0; bool is_local_node = mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !isPBCSlaveNode(n); bool count_node = is_local_node && is_not_pbc_slave_node; for (UInt i = 0; i < spatial_dimension; ++i) { if (count_node) mv2 += *vel_val * *vel_val * *mass_val; vel_val++; mass_val++; } ekin += mv2; } StaticCommunicator::getStaticCommunicator().allReduce(&ekin, 1, _so_sum); AKANTU_DEBUG_OUT(); return ekin * .5; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy(const ElementType & type, UInt index) { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = getFEEngine().getNbQuadraturePoints(type); Array vel_on_quad(nb_quadrature_points, spatial_dimension); Array filter_element(1, 1, index); getFEEngine().interpolateOnQuadraturePoints(*velocity, vel_on_quad, spatial_dimension, type, _not_ghost, filter_element); Array::vector_iterator vit = vel_on_quad.begin(spatial_dimension); Array::vector_iterator vend = vel_on_quad.end(spatial_dimension); Vector rho_v2(nb_quadrature_points); Real rho = materials[material_index(type)(index)]->getRho(); for (UInt q = 0; vit != vend; ++vit, ++q) { rho_v2(q) = rho * vit->dot(*vit); } AKANTU_DEBUG_OUT(); return .5*getFEEngine().integrate(rho_v2, type, index); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getExternalWork() { AKANTU_DEBUG_IN(); Real * velo = velocity->storage(); Real * forc = force->storage(); Real * resi = residual->storage(); bool * boun = blocked_dofs->storage(); Real work = 0.; UInt nb_nodes = mesh.getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !isPBCSlaveNode(n); bool count_node = is_local_node && is_not_pbc_slave_node; for (UInt i = 0; i < spatial_dimension; ++i) { if (count_node) { if(*boun) work -= *resi * *velo * time_step; else work += *forc * *velo * time_step; } ++velo; ++forc; ++resi; ++boun; } } StaticCommunicator::getStaticCommunicator().allReduce(&work, 1, _so_sum); AKANTU_DEBUG_OUT(); return work; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id) { AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(); } else if (energy_id == "external work"){ return getExternalWork(); } Real energy = 0.; std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { energy += (*mat_it)->getEnergy(energy_id); } /// reduction sum over all processors StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id, const ElementType & type, UInt index){ AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(type, index); } std::vector::iterator mat_it; UInt mat_index = this->material_index(type, _not_ghost)(index); UInt mat_loc_num = this->material_local_numbering(type, _not_ghost)(index); Real energy = this->materials[mat_index]->getEnergy(energy_id, type, mat_loc_num); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsAdded(const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); this->getFEEngine().initShapeFunctions(_not_ghost); this->getFEEngine().initShapeFunctions(_ghost); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = this->mesh.firstType(spatial_dimension, gt, _ek_not_defined); Mesh::type_iterator end = this->mesh.lastType(spatial_dimension, gt, _ek_not_defined); for(; it != end; ++it) { UInt nb_element = this->mesh.getNbElement(*it, gt); if(!material_index.exists(*it, gt)) { this->material_index .alloc(nb_element, 1, *it, gt); this->material_local_numbering.alloc(nb_element, 1, *it, gt); } else { this->material_index (*it, gt).resize(nb_element); this->material_local_numbering(*it, gt).resize(nb_element); } } } Array::const_iterator it = element_list.begin(); Array::const_iterator end = element_list.end(); ElementTypeMapArray filter("new_element_filter", this->getID()); for (UInt el = 0; it != end; ++it, ++el) { const Element & elem = *it; if(!filter.exists(elem.type, elem.ghost_type)) filter.alloc(0, 1, elem.type, elem.ghost_type); filter(elem.type, elem.ghost_type).push_back(elem.element); } this->assignMaterialToElements(&filter); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onElementsAdded(element_list, event); } if(method == _explicit_lumped_mass) this->assembleMassLumped(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsRemoved(__attribute__((unused)) const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) { this->getFEEngine().initShapeFunctions(_not_ghost); this->getFEEngine().initShapeFunctions(_ghost); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onElementsRemoved(element_list, new_numbering, event); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesAdded(const Array & nodes_list, __attribute__((unused)) const NewNodesEvent & event) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); if(displacement) displacement->resize(nb_nodes); if(mass ) mass ->resize(nb_nodes); if(velocity ) velocity ->resize(nb_nodes); if(acceleration) acceleration->resize(nb_nodes); if(force ) force ->resize(nb_nodes); if(residual ) residual ->resize(nb_nodes); if(blocked_dofs) blocked_dofs->resize(nb_nodes); if(previous_displacement) previous_displacement->resize(nb_nodes); if(increment_acceleration) increment_acceleration->resize(nb_nodes); if(increment) increment->resize(nb_nodes); if(current_position) current_position->resize(nb_nodes); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onNodesAdded(nodes_list, event); } + delete dof_synchronizer; + dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); + dof_synchronizer->initLocalDOFEquationNumbers(); + dof_synchronizer->initGlobalDOFEquationNumbers(); + + if (method != _explicit_lumped_mass) { + this->initSolver(); + } + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesRemoved(__attribute__((unused)) const Array & element_list, const Array & new_numbering, __attribute__((unused)) const RemovedNodesEvent & event) { if(displacement) mesh.removeNodesFromArray(*displacement, new_numbering); if(mass ) mesh.removeNodesFromArray(*mass , new_numbering); if(velocity ) mesh.removeNodesFromArray(*velocity , new_numbering); if(acceleration) mesh.removeNodesFromArray(*acceleration, new_numbering); if(force ) mesh.removeNodesFromArray(*force , new_numbering); if(residual ) mesh.removeNodesFromArray(*residual , new_numbering); if(blocked_dofs) mesh.removeNodesFromArray(*blocked_dofs, new_numbering); if(increment_acceleration) mesh.removeNodesFromArray(*increment_acceleration, new_numbering); if(increment) mesh.removeNodesFromArray(*increment , new_numbering); delete dof_synchronizer; dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); if (method != _explicit_lumped_mass) { this->initSolver(); } } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::isInternal(const std::string & field_name, const ElementKind & element_kind){ bool is_internal = false; /// check if at least one material contains field_id as an internal for (UInt m = 0; m < materials.size() && !is_internal; ++m) { is_internal |= materials[m]->isInternal(field_name, element_kind); } return is_internal; } /* -------------------------------------------------------------------------- */ ElementTypeMap SolidMechanicsModel::getInternalDataPerElem(const std::string & field_name, const ElementKind & element_kind){ if (!(this->isInternal(field_name,element_kind))) AKANTU_EXCEPTION("unknown internal " << field_name); for (UInt m = 0; m < materials.size() ; ++m) { if (materials[m]->isInternal(field_name, element_kind)) return materials[m]->getInternalDataPerElem(field_name, element_kind); } return ElementTypeMap(); } /* -------------------------------------------------------------------------- */ ElementTypeMapArray & SolidMechanicsModel::flattenInternal(const std::string & field_name, const ElementKind & kind, const GhostType ghost_type){ std::pair key(field_name,kind); if (this->registered_internals.count(key) == 0){ this->registered_internals[key] = new ElementTypeMapArray(field_name, this->id); } ElementTypeMapArray * internal_flat = this->registered_internals[key]; typedef ElementTypeMapArray::type_iterator iterator; iterator tit = internal_flat->firstType(spatial_dimension, ghost_type, kind); iterator end = internal_flat->lastType(spatial_dimension, ghost_type, kind); for (; tit != end; ++tit) { ElementType type = *tit; (*internal_flat)(type,ghost_type).clear(); } for (UInt m = 0; m < materials.size(); ++m) { if (materials[m]->isInternal(field_name, kind)) materials[m]->flattenInternal(field_name, *internal_flat, ghost_type, kind); } return *internal_flat; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::flattenAllRegisteredInternals(const ElementKind & kind){ std::map, ElementTypeMapArray *>::iterator it = this->registered_internals.begin(); std::map, ElementTypeMapArray *>::iterator end = this->registered_internals.end(); while (it != end){ if (kind == it->first.second) this->flattenInternal(it->first.first, kind); ++it; } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onDump(){ this->flattenAllRegisteredInternals(_ek_regular); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER dumper::Field * SolidMechanicsModel ::createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) { dumper::Field * field = NULL; if(field_name == "partitions") field = mesh.createElementalField(mesh.getConnectivities(),group_name,spatial_dimension,kind); else if(field_name == "material_index") field = mesh.createElementalField(material_index,group_name,spatial_dimension,kind); else { // this copy of field_name is used to compute derivated data such as // strain and von mises stress that are based on grad_u and stress std::string field_name_copy(field_name); if (field_name == "strain" || field_name == "Green strain" || field_name == "principal strain" || field_name == "principal Green strain") field_name_copy = "grad_u"; else if (field_name == "Von Mises stress") field_name_copy = "stress"; bool is_internal = this->isInternal(field_name_copy,kind); if (is_internal) { ElementTypeMap nb_data_per_elem = this->getInternalDataPerElem(field_name_copy, kind); ElementTypeMapArray & internal_flat = this->flattenInternal(field_name_copy,kind); field = mesh.createElementalField(internal_flat, group_name, spatial_dimension,kind,nb_data_per_elem); if (field_name == "strain"){ dumper::ComputeStrain * foo = new dumper::ComputeStrain(*this); field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); } else if (field_name == "Von Mises stress") { dumper::ComputeVonMisesStress * foo = new dumper::ComputeVonMisesStress(*this); field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); } else if (field_name == "Green strain") { dumper::ComputeStrain * foo = new dumper::ComputeStrain(*this); field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); } else if (field_name == "principal strain") { dumper::ComputePrincipalStrain * foo = new dumper::ComputePrincipalStrain(*this); field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); } else if (field_name == "principal Green strain") { dumper::ComputePrincipalStrain * foo = new dumper::ComputePrincipalStrain(*this); field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); } //treat the paddings if (padding_flag){ if (field_name == "stress"){ if (spatial_dimension == 2) { dumper::StressPadder<2> * foo = new dumper::StressPadder<2>(*this); field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); } } else if (field_name == "strain" || field_name == "Green strain"){ if (spatial_dimension == 2) { dumper::StrainPadder<2> * foo = new dumper::StrainPadder<2>(*this); field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); } } } // homogenize the field dumper::ComputeFunctorInterface * foo = dumper::HomogenizerProxy::createHomogenizer(*field); field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); } } return field; } /* -------------------------------------------------------------------------- */ dumper::Field * SolidMechanicsModel::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { std::map* > real_nodal_fields; real_nodal_fields["displacement" ] = displacement; real_nodal_fields["mass" ] = mass; real_nodal_fields["velocity" ] = velocity; real_nodal_fields["acceleration" ] = acceleration; real_nodal_fields["force" ] = force; real_nodal_fields["residual" ] = residual; real_nodal_fields["increment" ] = increment; dumper::Field * field = NULL; if (padding_flag) field = mesh.createNodalField(real_nodal_fields[field_name],group_name, 3); else field = mesh.createNodalField(real_nodal_fields[field_name],group_name); return field; } /* -------------------------------------------------------------------------- */ dumper::Field * SolidMechanicsModel::createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) { std::map* > uint_nodal_fields; uint_nodal_fields["blocked_dofs" ] = blocked_dofs; dumper::Field * field = NULL; field = mesh.createNodalField(uint_nodal_fields[field_name],group_name); return field; } /* -------------------------------------------------------------------------- */ #else /* -------------------------------------------------------------------------- */ dumper::Field * SolidMechanicsModel ::createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind){ return NULL; } /* -------------------------------------------------------------------------- */ dumper::Field * SolidMechanicsModel::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { return NULL; } /* -------------------------------------------------------------------------- */ dumper::Field * SolidMechanicsModel::createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) { return NULL; } #endif /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(const std::string & dumper_name) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(const std::string & dumper_name, UInt step) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(dumper_name, step); } /* ------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(const std::string & dumper_name, Real time, UInt step) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump() { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(UInt step) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(step); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(Real time, UInt step) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(time, step); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeCauchyStresses() { AKANTU_DEBUG_IN(); // call compute stiffness matrix on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; if(mat.isFiniteDeformation()) mat.computeAllCauchyStresses(_not_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::saveStressAndStrainBeforeDamage() { EventManager::sendEvent(SolidMechanicsModelEvent::BeginningOfDamageIterationEvent()); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateEnergiesAfterDamage() { EventManager::sendEvent(SolidMechanicsModelEvent::AfterDamageEvent()); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Solid Mechanics Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEEngine().printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + nodals information [" << std::endl; displacement->printself(stream, indent + 2); mass ->printself(stream, indent + 2); velocity ->printself(stream, indent + 2); acceleration->printself(stream, indent + 2); force ->printself(stream, indent + 2); residual ->printself(stream, indent + 2); blocked_dofs->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + material information [" << std::endl; material_index.printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + materials [" << std::endl; std::vector::const_iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { const Material & mat = *(*mat_it); mat.printself(stream, indent + 1); } stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc index aed7bf3bc..c55a288c8 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc @@ -1,755 +1,746 @@ /** * @file solid_mechanics_model_cohesive.cc * * @author Marco Vocialta * @author Nicolas Richart * * @date creation: Tue May 08 2012 * @date last modification: Fri Sep 05 2014 * * @brief Solid mechanics model for 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 . * */ /* -------------------------------------------------------------------------- */ #include #include "shape_cohesive.hh" #include "solid_mechanics_model_cohesive.hh" #include "dumpable_inline_impl.hh" #include "material_cohesive.hh" #ifdef AKANTU_USE_IOHELPER # include "dumper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ const SolidMechanicsModelCohesiveOptions default_solid_mechanics_model_cohesive_options(_explicit_lumped_mass, false, false); /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::SolidMechanicsModelCohesive(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : SolidMechanicsModel(mesh, dim, id, memory_id), tangents("tangents", id), facet_stress("facet_stress", id), facet_material("facet_material", id) { AKANTU_DEBUG_IN(); inserter = NULL; #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) facet_synchronizer = NULL; facet_stress_synchronizer = NULL; cohesive_distributed_synchronizer = NULL; global_connectivity = NULL; #endif delete material_selector; material_selector = new DefaultMaterialCohesiveSelector(*this); this->registerEventHandler(*this); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("cohesive elements", id); this->mesh.addDumpMeshToDumper("cohesive elements", mesh, spatial_dimension, _not_ghost, _ek_cohesive); #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::~SolidMechanicsModelCohesive() { AKANTU_DEBUG_IN(); delete inserter; #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) delete cohesive_distributed_synchronizer; delete facet_synchronizer; delete facet_stress_synchronizer; #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::setTimeStep(Real time_step) { SolidMechanicsModel::setTimeStep(time_step); #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper("cohesive elements").setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initFull(const ModelOptions & options) { AKANTU_DEBUG_IN(); const SolidMechanicsModelCohesiveOptions & smmc_options = dynamic_cast(options); this->is_extrinsic = smmc_options.extrinsic; if (!inserter) inserter = new CohesiveElementInserter(mesh, is_extrinsic, synch_parallel, id+":cohesive_element_inserter"); SolidMechanicsModel::initFull(options); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initMaterials() { AKANTU_DEBUG_IN(); // make sure the material are instantiated if(!are_materials_instantiated) instantiateMaterials(); /// find the first cohesive material UInt cohesive_index = 0; while ((dynamic_cast(materials[cohesive_index]) == NULL) && cohesive_index <= materials.size()) ++cohesive_index; AKANTU_DEBUG_ASSERT(cohesive_index != materials.size(), "No cohesive materials in the material input file"); material_selector->setFallback(cohesive_index); // set the facet information in the material in case of dynamic insertion if (is_extrinsic) { const Mesh & mesh_facets = inserter->getMeshFacets(); mesh_facets.initElementTypeMapArray(facet_material, 1, spatial_dimension - 1); Element element; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { element.ghost_type = *gt; Mesh::type_iterator first = mesh_facets.firstType(spatial_dimension - 1, *gt); Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1, *gt); for(;first != last; ++first) { element.type = *first; Array & f_material = facet_material(*first, *gt); UInt nb_element = mesh_facets.getNbElement(*first, *gt); f_material.resize(nb_element); f_material.set(cohesive_index); for (UInt el = 0; el < nb_element; ++el) { element.element = el; UInt mat_index = (*material_selector)(element); f_material(el) = mat_index; MaterialCohesive & mat = dynamic_cast(*materials[mat_index]); mat.addFacet(element); } } } SolidMechanicsModel::initMaterials(); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (facet_synchronizer != NULL) inserter->initParallel(facet_synchronizer, synch_parallel); #endif initAutomaticInsertion(); } else { try { /// In case of insertion of intrinsic elements along mesh physical surfaces const ParserSection & mesh_section = *(this->parser->getSubSections(_st_mesh).first); std::string cohesive_surfaces = mesh_section.getParameter("cohesive_surfaces"); initIntrinsicCohesiveMaterials(cohesive_surfaces); } catch(...) { /// Default intrinsic insertion initIntrinsicCohesiveMaterials(cohesive_index); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initIntrinsicCohesiveMaterials(std::string cohesive_surfaces) { AKANTU_DEBUG_IN(); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (facet_synchronizer != NULL) inserter->initParallel(facet_synchronizer, synch_parallel); #endif std::istringstream split(cohesive_surfaces); std::string physname; while(std::getline(split,physname,',')){ AKANTU_DEBUG_INFO("Pre-inserting cohesive elements along facets from physical surface: " << physname); insertElementsFromMeshData(physname); } #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (facet_synchronizer != NULL){ facet_synchronizer-> asynchronousSynchronize(*inserter, _gst_ce_groups); facet_synchronizer-> waitEndSynchronize(*inserter, _gst_ce_groups); } #endif SolidMechanicsModel::initMaterials(); if(is_default_material_selector) delete material_selector; material_selector = new MeshDataMaterialCohesiveSelector(*this); inserter->insertElements(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initIntrinsicCohesiveMaterials(UInt cohesive_index) { AKANTU_DEBUG_IN(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { Mesh::type_iterator first = mesh.firstType(spatial_dimension, *gt, _ek_cohesive); Mesh::type_iterator last = mesh.lastType(spatial_dimension, *gt, _ek_cohesive); for(;first != last; ++first) { Array & mat_indexes = this->material_index(*first, *gt); Array & mat_loc_num = this->material_local_numbering(*first, *gt); mat_indexes.set(cohesive_index); mat_loc_num.clear(); } } #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (facet_synchronizer != NULL) inserter->initParallel(facet_synchronizer, synch_parallel); #endif SolidMechanicsModel::initMaterials(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian * */ void SolidMechanicsModelCohesive::initModel() { AKANTU_DEBUG_IN(); SolidMechanicsModel::initModel(); registerFEEngineObject("CohesiveFEEngine", mesh, spatial_dimension); /// add cohesive type connectivity ElementType type = _not_defined; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType type_ghost = *gt; Mesh::type_iterator it = mesh.firstType(spatial_dimension, type_ghost); Mesh::type_iterator last = mesh.lastType(spatial_dimension, type_ghost); for (; it != last; ++it) { const Array & connectivity = mesh.getConnectivity(*it, type_ghost); if (connectivity.getSize() != 0) { type = *it; ElementType type_facet = Mesh::getFacetType(type); ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); mesh.addConnectivityType(type_cohesive, type_ghost); } } } AKANTU_DEBUG_ASSERT(type != _not_defined, "No elements in the mesh"); getFEEngine("CohesiveFEEngine").initShapeFunctions(_not_ghost); getFEEngine("CohesiveFEEngine").initShapeFunctions(_ghost); registerFEEngineObject("FacetsFEEngine", mesh.getMeshFacets(), spatial_dimension - 1); if (is_extrinsic) { getFEEngine("FacetsFEEngine").initShapeFunctions(_not_ghost); getFEEngine("FacetsFEEngine").initShapeFunctions(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::limitInsertion(BC::Axis axis, Real first_limit, Real second_limit) { AKANTU_DEBUG_IN(); inserter->setLimit(axis, first_limit, second_limit); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::insertIntrinsicElements() { AKANTU_DEBUG_IN(); inserter->insertIntrinsicElements(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::insertElementsFromMeshData(std::string physname) { AKANTU_DEBUG_IN(); UInt material_index = SolidMechanicsModel::getMaterialIndex(physname); inserter->insertIntrinsicElements(physname, material_index); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initAutomaticInsertion() { AKANTU_DEBUG_IN(); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (facet_stress_synchronizer != NULL) { DataAccessor * data_accessor = this; const ElementTypeMapArray & rank_to_element = synch_parallel->getPrankToElement(); facet_stress_synchronizer->updateFacetStressSynchronizer(*inserter, rank_to_element, *data_accessor); } #endif inserter->getMeshFacets().initElementTypeMapArray(facet_stress, 2 * spatial_dimension * spatial_dimension, spatial_dimension - 1); resizeFacetStress(); /// compute normals on facets computeNormals(); initStressInterpolation(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::updateAutomaticInsertion() { AKANTU_DEBUG_IN(); inserter->limitCheckFacets(); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (facet_stress_synchronizer != NULL) { DataAccessor * data_accessor = this; const ElementTypeMapArray & rank_to_element = synch_parallel->getPrankToElement(); facet_stress_synchronizer->updateFacetStressSynchronizer(*inserter, rank_to_element, *data_accessor); } #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initStressInterpolation() { Mesh & mesh_facets = inserter->getMeshFacets(); /// compute quadrature points coordinates on facets Array & position = mesh.getNodes(); ElementTypeMapArray quad_facets("quad_facets", id); mesh_facets.initElementTypeMapArray(quad_facets, spatial_dimension, spatial_dimension - 1); getFEEngine("FacetsFEEngine").interpolateOnQuadraturePoints(position, quad_facets); /// compute elements quadrature point positions and build /// element-facet quadrature points data structure ElementTypeMapArray elements_quad_facets("elements_quad_facets", id); mesh.initElementTypeMapArray(elements_quad_facets, spatial_dimension, spatial_dimension); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType elem_gt = *gt; Mesh::type_iterator it = mesh.firstType(spatial_dimension, elem_gt); Mesh::type_iterator last = mesh.lastType(spatial_dimension, elem_gt); for (; it != last; ++it) { ElementType type = *it; UInt nb_element = mesh.getNbElement(type, elem_gt); if (nb_element == 0) continue; /// compute elements' quadrature points and list of facet /// quadrature points positions by element Array & facet_to_element = mesh_facets.getSubelementToElement(type, elem_gt); UInt nb_facet_per_elem = facet_to_element.getNbComponent(); Array & el_q_facet = elements_quad_facets(type, elem_gt); ElementType facet_type = Mesh::getFacetType(type); UInt nb_quad_per_facet = getFEEngine("FacetsFEEngine").getNbQuadraturePoints(facet_type); el_q_facet.resize(nb_element * nb_facet_per_elem * nb_quad_per_facet); for (UInt el = 0; el < nb_element; ++el) { for (UInt f = 0; f < nb_facet_per_elem; ++f) { Element global_facet_elem = facet_to_element(el, f); UInt global_facet = global_facet_elem.element; GhostType facet_gt = global_facet_elem.ghost_type; const Array & quad_f = quad_facets(facet_type, facet_gt); for (UInt q = 0; q < nb_quad_per_facet; ++q) { for (UInt s = 0; s < spatial_dimension; ++s) { el_q_facet(el * nb_facet_per_elem * nb_quad_per_facet + f * nb_quad_per_facet + q, s) = quad_f(global_facet * nb_quad_per_facet + q, s); } } } } } } /// loop over non cohesive materials for (UInt m = 0; m < materials.size(); ++m) { try { MaterialCohesive & mat __attribute__((unused)) = dynamic_cast(*materials[m]); } catch(std::bad_cast&) { /// initialize the interpolation function materials[m]->initElementalFieldInterpolation(elements_quad_facets); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::updateResidual(bool need_initialize) { AKANTU_DEBUG_IN(); if (need_initialize) initializeUpdateResidualData(); // f -= fint std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { try { MaterialCohesive & mat = dynamic_cast(**mat_it); mat.computeTraction(_not_ghost); } catch (std::bad_cast & bce) { } } SolidMechanicsModel::updateResidual(false); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { try { MaterialCohesive & mat = dynamic_cast(**mat_it); mat.computeEnergies(); } catch (std::bad_cast & bce) { } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::computeNormals() { AKANTU_DEBUG_IN(); Mesh & mesh_facets = this->inserter->getMeshFacets(); this->getFEEngine("FacetsFEEngine").computeNormalsOnControlPoints(_not_ghost); /** * @todo store tangents while computing normals instead of * recomputing them as follows: */ /* ------------------------------------------------------------------------ */ UInt tangent_components = spatial_dimension * (spatial_dimension - 1); mesh_facets.initElementTypeMapArray(tangents, tangent_components, spatial_dimension - 1); Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1); Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1); for (; it != last; ++it) { ElementType facet_type = *it; const Array & normals = this->getFEEngine("FacetsFEEngine").getNormalsOnQuadPoints(facet_type); Array & tangents = this->tangents(facet_type); Math::compute_tangents(normals, tangents); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::interpolateStress() { ElementTypeMapArray by_elem_result("temporary_stress_by_facets", id); for (UInt m = 0; m < materials.size(); ++m) { try { MaterialCohesive & mat __attribute__((unused)) = dynamic_cast(*materials[m]); } catch(std::bad_cast&) { /// interpolate stress on facet quadrature points positions materials[m]->interpolateStressOnFacets(facet_stress, by_elem_result); } } #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData(debug::_dm_model_cohesive, "Interpolated stresses before", facet_stress); #endif synch_registry->synchronize(_gst_smmc_facets_stress); #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData(debug::_dm_model_cohesive, "Interpolated stresses", facet_stress); #endif } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModelCohesive::checkCohesiveStress() { interpolateStress(); for (UInt m = 0; m < materials.size(); ++m) { try { MaterialCohesive & mat_cohesive = dynamic_cast(*materials[m]); /// check which not ghost cohesive elements are to be created mat_cohesive.checkInsertion(); } catch(std::bad_cast&) { } } /* if(static and extrinsic) { check max mean stresses and change inserter.getInsertionFacets(type_facet); } */ /// communicate data among processors synch_registry->synchronize(_gst_smmc_facets); /// insert cohesive elements return inserter->insertElements(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onElementsAdded(const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) updateCohesiveSynchronizers(); #endif SolidMechanicsModel::onElementsAdded(element_list, event); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (cohesive_distributed_synchronizer != NULL) cohesive_distributed_synchronizer->computeAllBufferSizes(*this); #endif /// update shape functions getFEEngine("CohesiveFEEngine").initShapeFunctions(_not_ghost); getFEEngine("CohesiveFEEngine").initShapeFunctions(_ghost); if (is_extrinsic) resizeFacetStress(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onNodesAdded(const Array & doubled_nodes, __attribute__((unused)) const NewNodesEvent & event) { AKANTU_DEBUG_IN(); UInt nb_new_nodes = doubled_nodes.getSize(); Array nodes_list(nb_new_nodes); for (UInt n = 0; n < nb_new_nodes; ++n) nodes_list(n) = doubled_nodes(n, 1); SolidMechanicsModel::onNodesAdded(nodes_list, event); for (UInt n = 0; n < nb_new_nodes; ++n) { UInt old_node = doubled_nodes(n, 0); UInt new_node = doubled_nodes(n, 1); for (UInt dim = 0; dim < spatial_dimension; ++dim) { (*displacement)(new_node, dim) = (*displacement)(old_node, dim); (*velocity) (new_node, dim) = (*velocity) (old_node, dim); (*acceleration)(new_node, dim) = (*acceleration)(old_node, dim); (*blocked_dofs)(new_node, dim) = (*blocked_dofs)(old_node, dim); if (current_position) (*current_position)(new_node, dim) = (*current_position)(old_node, dim); if (increment_acceleration) (*increment_acceleration)(new_node, dim) = (*increment_acceleration)(old_node, dim); if (increment) (*increment)(new_node, dim) = (*increment)(old_node, dim); if (previous_displacement) (*previous_displacement)(new_node, dim) = (*previous_displacement)(old_node, dim); } } - delete dof_synchronizer; - dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); - dof_synchronizer->initLocalDOFEquationNumbers(); - dof_synchronizer->initGlobalDOFEquationNumbers(); - - if (method != _explicit_lumped_mass) { - this->initSolver(); - } - AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onEndSolveStep(const AnalysisMethod & method) { AKANTU_DEBUG_IN(); /****************************************************************************** This is required because the Cauchy stress is the stress measure that is used to check the insertion of cohesive elements ******************************************************************************/ std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; if(mat.isFiniteDeformation()) mat.computeAllCauchyStresses(_not_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "SolidMechanicsModelCohesive [" << std::endl; SolidMechanicsModel::printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::resizeFacetStress() { AKANTU_DEBUG_IN(); Mesh & mesh_facets = inserter->getMeshFacets(); 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; UInt nb_facet = mesh_facets.getNbElement(type, ghost_type); UInt nb_quadrature_points = getFEEngine("FacetsFEEngine").getNbQuadraturePoints(type, ghost_type); UInt new_size = nb_facet * nb_quadrature_points; facet_stress(type, ghost_type).resize(new_size); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, const ElementKind & element_kind, bool padding_flag) { AKANTU_DEBUG_IN(); UInt spatial_dimension = this->spatial_dimension; ElementKind _element_kind = element_kind; if (dumper_name == "cohesive elements") { _element_kind = _ek_cohesive; } else if (dumper_name == "facets") { spatial_dimension = this->spatial_dimension - 1; } SolidMechanicsModel::addDumpGroupFieldToDumper(dumper_name, field_id, group_name, spatial_dimension, _element_kind, padding_flag); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onDump(){ this->flattenAllRegisteredInternals(_ek_cohesive); SolidMechanicsModel::onDump(); } /* -------------------------------------------------------------------------- */ __END_AKANTU__