diff --git a/src/geometry/mesh_abstract_intersector.hh b/src/geometry/mesh_abstract_intersector.hh index 3ddb57f20..f4d3c846c 100644 --- a/src/geometry/mesh_abstract_intersector.hh +++ b/src/geometry/mesh_abstract_intersector.hh @@ -1,108 +1,108 @@ /** * @file mesh_abstract_intersector.hh * * @author Lucas Frerot * * @date creation: Mon Jul 13 2015 * @date last modification: Mon Jul 13 2015 * * @brief Abstract class for intersection computations * * @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_ABSTRACT_INTERSECTOR_HH__ #define __AKANTU_MESH_ABSTRACT_INTERSECTOR_HH__ #include "aka_common.hh" #include "mesh_geom_abstract.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /** * @brief Class used to perform intersections on a mesh and construct output data */ template class MeshAbstractIntersector : public MeshGeomAbstract { public: /// Construct from mesh explicit MeshAbstractIntersector(Mesh & mesh, const ID & id = "mesh_abstract_intersector", const MemoryID & memory_id = 0); /// Destructor virtual ~MeshAbstractIntersector(); public: /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the new_node_per_elem array AKANTU_GET_MACRO(NewNodePerElem, new_node_per_elem, const ElementTypeMapUInt &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(NewNodePerElem, new_node_per_elem, UInt); - /// get the new_nodes array - AKANTU_GET_MACRO(NewNodes, new_nodes, const Array *); + /// get the intersection_points array + AKANTU_GET_MACRO(IntersectionPoints, intersection_points, const Array *); /// get the nb_seg_by_el UInt AKANTU_GET_MACRO(NbSegByEl, nb_seg_by_el, const UInt); /** * @brief Compute the intersection with a query object * * This function needs to be implemented for every subclass. It computes the intersections * with the tree of primitives and creates the data for the user. * * @param query the CGAL primitive of the query object */ virtual void computeIntersectionQuery(const Query & query) = 0; /// Compute intersection points between the mesh primitives (segments) and a query (surface in 3D or a curve in 2D), double intersection points for the same primitives are not considered virtual void computeMeshQueryIntersectionPoint(const Query & query) = 0; /// Compute intersection between the mesh and a list of queries virtual void computeIntersectionQueryList(const std::list & query_list); /// Compute intersection points between the mesh and a list of queries virtual void computeMeshQueryListIntersectionPoint(const std::list & query_list); /// Compute whatever result is needed from the user virtual void buildResultFromQueryList(const std::list & query_list) = 0; protected: /// new node per element (column 0: number of new nodes, then odd is the intersection node number and even the ID of the intersected segment) ElementTypeMapUInt new_node_per_elem; /// intersection output: new intersection points (computeMeshQueryListIntersectionPoint) - Array * new_nodes; + Array * intersection_points; /// number of segment in a considered element of the templated type of element specialized intersector const UInt nb_seg_by_el; }; __END_AKANTU__ #include "mesh_abstract_intersector_tmpl.hh" #endif // __AKANTU_MESH_ABSTRACT_INTERSECTOR_HH__ diff --git a/src/geometry/mesh_abstract_intersector_tmpl.hh b/src/geometry/mesh_abstract_intersector_tmpl.hh index 412cbad70..df037f012 100644 --- a/src/geometry/mesh_abstract_intersector_tmpl.hh +++ b/src/geometry/mesh_abstract_intersector_tmpl.hh @@ -1,89 +1,89 @@ /** * @file mesh_abstract_intersector_tmpl.hh * * @author Lucas Frerot * * @date creation: Mon Jul 13 2015 * @date last modification: Mon Jul 13 2015 * * @brief General class for intersection computations * * @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_ABSTRACT_INTERSECTOR_TMPL_HH__ #define __AKANTU_MESH_ABSTRACT_INTERSECTOR_TMPL_HH__ #include "aka_common.hh" #include "mesh_abstract_intersector.hh" __BEGIN_AKANTU__ template MeshAbstractIntersector::MeshAbstractIntersector(Mesh & mesh, const ID & id, const MemoryID & memory_id): MeshGeomAbstract(mesh, id, memory_id), new_node_per_elem("new_node_per_elem", id), - new_nodes(NULL), + intersection_points(NULL), nb_seg_by_el(0) {} template MeshAbstractIntersector::~MeshAbstractIntersector() {} template void MeshAbstractIntersector::computeIntersectionQueryList( const std::list & query_list) { AKANTU_DEBUG_IN(); typename std::list::const_iterator query_it = query_list.begin(), query_end = query_list.end(); for (; query_it != query_end ; ++query_it) { computeIntersectionQuery(*query_it); } AKANTU_DEBUG_OUT(); } template void MeshAbstractIntersector::computeMeshQueryListIntersectionPoint( const std::list & query_list) { AKANTU_DEBUG_IN(); typename std::list::const_iterator query_it = query_list.begin(), query_end = query_list.end(); for (; query_it != query_end ; ++query_it) { computeMeshQueryIntersectionPoint(*query_it); } AKANTU_DEBUG_OUT(); } __END_AKANTU__ #endif // __AKANTU_MESH_ABSTRACT_INTERSECTOR_TMPL_HH__ diff --git a/src/geometry/mesh_geom_intersector_tmpl.hh b/src/geometry/mesh_geom_intersector_tmpl.hh index b4ebdbb00..f378627c8 100644 --- a/src/geometry/mesh_geom_intersector_tmpl.hh +++ b/src/geometry/mesh_geom_intersector_tmpl.hh @@ -1,64 +1,64 @@ /** * @file mesh_geom_intersector_tmpl.hh * * @author Lucas Frerot * * @date creation: Wed Apr 29 2015 * @date last modification: Wed Apr 29 2015 * * @brief General class for intersection computations * * @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_GEOM_INTERSECTOR_TMPL_HH__ #define __AKANTU_MESH_GEOM_INTERSECTOR_TMPL_HH__ #include "aka_common.hh" #include "mesh_geom_intersector.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ template MeshGeomIntersector::MeshGeomIntersector(Mesh & mesh, const ID & id, const MemoryID & memory_id) : MeshAbstractIntersector(mesh, id, memory_id), factory(mesh) {} template MeshGeomIntersector::~MeshGeomIntersector() {} template void MeshGeomIntersector::constructData() { - this->new_nodes->resize(0); + this->intersection_points->resize(0); factory.constructData(); } __END_AKANTU__ #endif // __AKANTU_MESH_GEOM_INTERSECTOR_TMPL_HH__ diff --git a/src/geometry/mesh_segment_intersector_tmpl.hh b/src/geometry/mesh_segment_intersector_tmpl.hh index 64cc4c124..0e96b40a2 100644 --- a/src/geometry/mesh_segment_intersector_tmpl.hh +++ b/src/geometry/mesh_segment_intersector_tmpl.hh @@ -1,267 +1,267 @@ /** * @file mesh_segment_intersector_tmpl.hh * * @author Lucas Frerot * * @date creation: Wed Apr 29 2015 * @date last modification: Wed Apr 29 2015 * * @brief Computation of mesh intersection with segments * * @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_SEGMENT_INTERSECTOR_TMPL_HH__ #define __AKANTU_MESH_SEGMENT_INTERSECTOR_TMPL_HH__ #include "aka_common.hh" #include "mesh_geom_common.hh" #include "tree_type_helper.hh" __BEGIN_AKANTU__ template MeshSegmentIntersector::MeshSegmentIntersector(Mesh & mesh, Mesh & result_mesh): parent_type(mesh), result_mesh(result_mesh), current_physical_name() { - this->new_nodes = new Array(0,dim); + this->intersection_points = new Array(0,dim); this->constructData(); } template MeshSegmentIntersector::~MeshSegmentIntersector() {} template void MeshSegmentIntersector::computeIntersectionQuery(const K::Segment_3 & query) { AKANTU_DEBUG_IN(); result_mesh.addConnectivityType(_segment_2, _not_ghost); result_mesh.addConnectivityType(_segment_2, _ghost); std::list result_list; std::set, segmentPairsLess> segment_set; this->factory.getTree().all_intersections(query, std::back_inserter(result_list)); this->computeSegments(result_list, segment_set, query); // Arrays for storing nodes and connectivity Array & nodes = result_mesh.getNodes(); Array & connectivity = result_mesh.getConnectivity(_segment_2); // Arrays for storing associated element and physical name bool valid_elemental_data = true; Array * associated_element = NULL; Array * associated_physical_name = NULL; try { associated_element = &result_mesh.getData("associated_element", _segment_2); associated_physical_name = &result_mesh.getData("physical_names", _segment_2); } catch (debug::Exception & e) { valid_elemental_data = false; } std::set::iterator it = segment_set.begin(), end = segment_set.end(); // Loop over the segment pairs for (; it != end ; ++it) { if (!it->first.is_degenerate()) { Vector segment_connectivity(2); segment_connectivity(0) = result_mesh.getNbNodes(); segment_connectivity(1) = result_mesh.getNbNodes() + 1; connectivity.push_back(segment_connectivity); // Copy nodes Vector source(dim), target(dim); for (UInt j = 0 ; j < dim ; j++) { source(j) = it->first.source()[j]; target(j) = it->first.target()[j]; } nodes.push_back(source); nodes.push_back(target); // Copy associated element info if (valid_elemental_data) { associated_element->push_back(Element(type, it->second)); associated_physical_name->push_back(current_physical_name); } } } AKANTU_DEBUG_OUT(); } template void MeshSegmentIntersector:: computeMeshQueryIntersectionPoint(const K::Segment_3 & query) { AKANTU_DEBUG_ERROR("The method: computeMeshQueryIntersectionPoint has not been implemented in class MeshSegmentIntersector!"); } template void MeshSegmentIntersector::buildResultFromQueryList(const std::list & query_list) { AKANTU_DEBUG_IN(); this->computeIntersectionQueryList(query_list); AKANTU_DEBUG_OUT(); } template void MeshSegmentIntersector::computeSegments(const std::list & intersections, std::set & segments, const K::Segment_3 & query) { AKANTU_DEBUG_IN(); /* * Number of intersections = 0 means * * - query is completely outside mesh * - query is completely inside primitive * * We try to determine the case and still construct the segment list */ if (intersections.size() == 0) { // We look at all the primitives intersected by two rays // If there is one primitive in common, then query is inside // that primitive K::Ray_3 ray1(query.source(), query.target()); K::Ray_3 ray2(query.target(), query.source()); std::set ray1_results, ray2_results; this->factory.getTree().all_intersected_primitives(ray1, std::inserter(ray1_results, ray1_results.begin())); this->factory.getTree().all_intersected_primitives(ray2, std::inserter(ray2_results, ray2_results.begin())); bool inside_primitive = false; UInt primitive_id = 0; std::set::iterator ray2_it = ray2_results.begin(), ray2_end = ray2_results.end(); // Test if first list contains an element of second list for (; ray2_it != ray2_end && !inside_primitive ; ++ray2_it) { if (ray1_results.find(*ray2_it) != ray1_results.end()) { inside_primitive = true; primitive_id = *ray2_it; } } if (inside_primitive) { segments.insert(std::make_pair(query, primitive_id)); } } else { typename std::list::const_iterator it = intersections.begin(), end = intersections.end(); for(; it != end ; ++it) { UInt el = (*it)->second; // Result of intersection is a segment if (const K::Segment_3 * segment = boost::get(&((*it)->first))) { // Check if the segment was alread created segments.insert(std::make_pair(*segment, el)); } // Result of intersection is a point else if (const K::Point_3 * point = boost::get(&((*it)->first))) { // We only want to treat points differently if we're in 3D with Tetra4 elements // This should be optimized by compilator if (dim == 3 && type == _tetrahedron_4) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); TreeTypeHelper, K>::container_type facets; const Array & nodes = this->mesh.getNodes(); Array::const_vector_iterator connectivity_vec = this->mesh.getConnectivity(type).begin(nb_nodes_per_element); const Vector & el_connectivity = connectivity_vec[el]; Matrix node_coordinates(dim, nb_nodes_per_element); for (UInt i = 0 ; i < nb_nodes_per_element ; i++) for (UInt j = 0 ; j < dim ; j++) node_coordinates(j, i) = nodes(el_connectivity(i), j); this->factory.addPrimitive(node_coordinates, el, facets); // Local tree TreeTypeHelper, K>::tree * local_tree = new TreeTypeHelper, K>::tree(facets.begin(), facets.end()); // Compute local intersections (with current element) std::list local_intersections; local_tree->all_intersections(query, std::back_inserter(local_intersections)); bool out_point_found = false; typename std::list::const_iterator local_it = local_intersections.begin(), local_end = local_intersections.end(); for (; local_it != local_end ; ++local_it) { if (const K::Point_3 * local_point = boost::get(&((*local_it)->first))) { if (!comparePoints(*point, *local_point)) { K::Segment_3 seg(*point, *local_point); segments.insert(std::make_pair(seg, el)); out_point_found = true; } } } if (!out_point_found) { TreeTypeHelper, K>::point_type a(node_coordinates(0, 0), node_coordinates(1, 0), node_coordinates(2, 0)), b(node_coordinates(0, 1), node_coordinates(1, 1), node_coordinates(2, 1)), c(node_coordinates(0, 2), node_coordinates(1, 2), node_coordinates(2, 2)), d(node_coordinates(0, 3), node_coordinates(1, 3), node_coordinates(2, 3)); K::Tetrahedron_3 tetra(a, b, c, d); const K::Point_3 * inside_point = NULL; if (tetra.has_on_bounded_side(query.source()) && !tetra.has_on_boundary(query.source())) inside_point = &query.source(); else if (tetra.has_on_bounded_side(query.target()) && !tetra.has_on_boundary(query.target())) inside_point = &query.target(); if (inside_point) { K::Segment_3 seg(*inside_point, *point); segments.insert(std::make_pair(seg, el)); } } delete local_tree; } } } } AKANTU_DEBUG_OUT(); } __END_AKANTU__ #endif // __AKANTU_MESH_SEGMENT_INTERSECTOR_TMPL_HH__ diff --git a/src/geometry/mesh_sphere_intersector_tmpl.hh b/src/geometry/mesh_sphere_intersector_tmpl.hh index 310caad48..ed60e7fc6 100644 --- a/src/geometry/mesh_sphere_intersector_tmpl.hh +++ b/src/geometry/mesh_sphere_intersector_tmpl.hh @@ -1,573 +1,573 @@ /** * @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, const ID & id, const MemoryID & memory_id): parent_type(mesh, id, memory_id), tol_intersection_on_node(1e-10), nb_nodes_fem(mesh.getNodes().getSize()), nb_prim_by_el(0) { - this->new_nodes = new Array(0,dim); + this->intersection_points = new Array(0,dim); for(Mesh::type_iterator it = mesh.firstType(dim); it != mesh.lastType(dim); ++it){ #if defined(AKANTU_IGFEM) if(*it == _triangle_3){ // Addition of the igfem types in the mesh this->mesh.addConnectivityType(_igfem_triangle_4, _not_ghost); this->mesh.addConnectivityType(_igfem_triangle_4, _ghost); this->mesh.addConnectivityType(_igfem_triangle_5, _not_ghost); this->mesh.addConnectivityType(_igfem_triangle_5, _ghost); } else if( (*it != _triangle_3) && (*it != _igfem_triangle_4) && (*it != _igfem_triangle_5) ) { AKANTU_DEBUG_ERROR("Not ready for mesh type " << *it); } 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( (*it != _triangle_3) ) AKANTU_DEBUG_ERROR("Not ready for mesh type " << *it); this->nb_prim_by_el = 3; #endif } for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; UInt nb_comp = 1 + nb_prim_by_el * 2; this->new_node_per_elem.alloc(0, nb_comp, type, ghost_type, 0); // std::cout << "new_node_per_elem size = " << this->new_node_per_elem(type,ghost_type).getSize() << " (type = " << type << ", ghost_type = " << ghost_type << ")" << std::endl; } this->constructData(); } template MeshSphereIntersector::~MeshSphereIntersector() {} template void MeshSphereIntersector::constructData() { for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; this->new_node_per_elem(type, ghost_type).resize(this->mesh.getNbElement(type, ghost_type)); this->new_node_per_elem(type, ghost_type).clear(); } MeshGeomIntersector, SK::Sphere_3, SK>::constructData(); } template void MeshSphereIntersector::computeIntersectionQuery(const SK::Sphere_3 & query) { AKANTU_DEBUG_IN(); Array & nodes = this->mesh.getNodes(); UInt nb_node = nodes.getSize(); UInt nb_not_ghost_elements = this->mesh.getNbElement(type); // Tolerance for proximity checks should be defined by user 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(); #if defined(AKANTU_IGFEM) NewIGFEMNodesEvent new_nodes; #else NewNodesEvent new_nodes; #endif 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_nodes_fem-1; for (; n < nb_node ; ++n) { Array::vector_iterator existing_node = nodes.begin(dim) + n; if (Math::are_vector_equal(dim, new_node.storage(), existing_node->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) { nodes.push_back(new_node); new_nodes.getList().push_back(nb_node); nb_node++; } // deduce ghost type and element id UInt element_id = it->id(); GhostType ghost_type = _not_ghost; if (element_id >= nb_not_ghost_elements) { element_id -= nb_not_ghost_elements; ghost_type = _ghost; } Array & new_node_per_elem_array = this->new_node_per_elem(type, ghost_type); if (!is_on_mesh) { new_node_per_elem_array(element_id, 0) += 1; new_node_per_elem_array(element_id, (2 * new_node_per_elem_array(element_id, 0)) - 1) = n; new_node_per_elem_array(element_id, 2 * new_node_per_elem_array(element_id, 0)) = 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())) { new_node_per_elem_array(element_id, (new_node_per_elem_array.getNbComponent() - 2)) = it->segId() - 1; } else { new_node_per_elem_array(element_id, (new_node_per_elem_array.getNbComponent() - 2)) = it->segId() % this->nb_prim_by_el; } } } } } } if(new_nodes.getList().getSize()) { #if defined(AKANTU_IGFEM) new_nodes.setNewNodePerElem(this->new_node_per_elem); new_nodes.setType(type); #endif this->mesh.sendEvent(new_nodes); } AKANTU_DEBUG_OUT(); } template void MeshSphereIntersector:: computeMeshQueryIntersectionPoint(const SK::Sphere_3 & query) { /// 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->new_nodes->getSize(); + UInt nb_node = nodes.getSize() + this->intersection_points->getSize(); UInt nb_not_ghost_elements = this->mesh.getNbElement(type); // Tolerance for proximity checks should be defined by user 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_nodes_fem; 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 new_nodes_it = this->new_nodes->begin(dim); - Array::vector_iterator new_nodes_end = this->new_nodes->end(dim); - for (; new_nodes_it != new_nodes_end ; ++new_nodes_it, ++n) { - if (Math::are_vector_equal(dim, new_node.storage(), new_nodes_it->storage())) { + 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->new_nodes->push_back(new_node); + this->intersection_points->push_back(new_node); nb_node++; } // deduce ghost type and element id UInt element_id = it->id(); GhostType ghost_type = _not_ghost; if (element_id >= nb_not_ghost_elements) { element_id -= nb_not_ghost_elements; ghost_type = _ghost; } Array & new_node_per_elem_array = this->new_node_per_elem(type, ghost_type); if (!is_on_mesh) { new_node_per_elem_array(element_id, 0) += 1; new_node_per_elem_array(element_id, (2 * new_node_per_elem_array(element_id, 0)) - 1) = n; new_node_per_elem_array(element_id, 2 * new_node_per_elem_array(element_id, 0)) = 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())) { new_node_per_elem_array(element_id, (new_node_per_elem_array.getNbComponent() - 2)) = it->segId() - 1; } else { new_node_per_elem_array(element_id, (new_node_per_elem_array.getNbComponent() - 2)) = it->segId() % this->nb_prim_by_el; } } } } } } AKANTU_DEBUG_OUT(); } template void MeshSphereIntersector::removeAdditionnalNodes() { AKANTU_DEBUG_IN(); RemovedNodesEvent remove_nodes(this->mesh); Array & nodes_removed = remove_nodes.getList(); Array & new_numbering = remove_nodes.getNewNumbering(); for(UInt i = 0 ; i < nb_nodes_fem ; ++i){ new_numbering(i) = i; } for(UInt i = nb_nodes_fem ; i < this->mesh.getNodes().getSize(); ++i){ new_numbering(i) = UInt(-1) ; nodes_removed.push_back(i); } this->mesh.sendEvent(remove_nodes); AKANTU_DEBUG_OUT(); } #if defined(AKANTU_IGFEM) template void MeshSphereIntersector::buildResultFromQueryList(const std::list & query_list) { AKANTU_DEBUG_IN(); this->computeIntersectionQueryList(query_list); NewIGFEMElementsEvent new_elements; UInt total_new_elements = 0; Array & new_elements_list = new_elements.getList(); Array & old_elements_list = new_elements.getOldElementsList(); RemovedElementsEvent remove_elem(this->mesh); Int nb_proc = StaticCommunicator::getStaticCommunicator().getNbProc(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; if (nb_proc == 1 && ghost_type == _ghost) continue; if (!this->mesh.getConnectivities().exists(type, ghost_type)) continue; UInt n_new_el = 0; Array & connec_type_tmpl = this->mesh.getConnectivity(type, ghost_type); Array & connec_igfem_tri4 = this->mesh.getConnectivity(_igfem_triangle_4, ghost_type), & connec_igfem_tri5 = this->mesh.getConnectivity(_igfem_triangle_5, ghost_type), & connec_tri3 = this->mesh.getConnectivity(_triangle_3, ghost_type); Element element_tri3(_triangle_3, 0, ghost_type, _ek_regular); Element element_tri4(_igfem_triangle_4, 0, ghost_type, _ek_igfem); Element element_tri5(_igfem_triangle_5, 0, ghost_type, _ek_igfem); remove_elem.getNewNumbering().alloc(connec_type_tmpl.getSize(), 1, type, ghost_type); Array & new_numbering = remove_elem.getNewNumbering(type, ghost_type); UInt new_nb_type_tmpl = 0; // Meter of the number of element (type) will we loop on elements Array & new_node_per_elem_array = this->new_node_per_elem(type, ghost_type); for (UInt nel = 0 ; nel != new_node_per_elem_array.getSize(); ++nel) { if( (type != _triangle_3) && ( (new_node_per_elem_array(nel,0)==0) || ( (new_node_per_elem_array(nel,0) == 1) && ( ( (new_node_per_elem_array(nel,2)+1) % this->nb_prim_by_el ) != new_node_per_elem_array(nel, new_node_per_elem_array.getNbComponent() - 2) ) ) ) ){ Element element_type_tmpl(type, 0, ghost_type, Mesh::getKind(type)); new_elements_list.resize(n_new_el+1); old_elements_list.resize(n_new_el+1); Vector ctri3(3); ctri3(0) = connec_type_tmpl(nel,0); ctri3(1) = connec_type_tmpl(nel,1); ctri3(2) = connec_type_tmpl(nel,2); /// add the new element in the mesh UInt nb_tri3 = connec_tri3.getSize(); connec_tri3.push_back(ctri3); element_tri3.element = nb_tri3; new_elements_list(n_new_el) = element_tri3; /// the number of the old element in the mesh element_type_tmpl.element = nel; old_elements_list(n_new_el) = element_type_tmpl; new_numbering(nel) = UInt(-1); ++n_new_el; remove_elem.getList().push_back(element_type_tmpl); } else if( (new_node_per_elem_array(nel,0)!=0) && !( (new_node_per_elem_array(nel,0) == 1) && ( ( (new_node_per_elem_array(nel,2)+1) % this->nb_prim_by_el ) != new_node_per_elem_array(nel, new_node_per_elem_array.getNbComponent() - 2) ) ) ){ Element element_type_tmpl(type, 0, ghost_type); element_type_tmpl.kind = Mesh::getKind(type); new_elements_list.resize(n_new_el+1); old_elements_list.resize(n_new_el+1); switch(new_node_per_elem_array(nel,0)){ case 1 :{ Vector ctri4(4); switch(new_node_per_elem_array(nel,2)){ case 1 : ctri4(0) = connec_type_tmpl(nel,2); ctri4(1) = connec_type_tmpl(nel,0); ctri4(2) = connec_type_tmpl(nel,1); break; case 2 : ctri4(0) = connec_type_tmpl(nel,0); ctri4(1) = connec_type_tmpl(nel,1); ctri4(2) = connec_type_tmpl(nel,2); break; case 3 : ctri4(0) = connec_type_tmpl(nel,1); ctri4(1) = connec_type_tmpl(nel,2); ctri4(2) = connec_type_tmpl(nel,0); break; default : AKANTU_DEBUG_ERROR("A triangle have only 3 segments not "<< new_node_per_elem_array(nel,2)); break; } ctri4(3) = new_node_per_elem_array(nel,1); UInt nb_tri4 = connec_igfem_tri4.getSize(); connec_igfem_tri4.push_back(ctri4); element_tri4.element = nb_tri4; //new_elements.getList().push_back(element_tri4); new_elements_list(n_new_el) = element_tri4; if(type == _igfem_triangle_4) new_numbering.push_back(connec_igfem_tri4.getSize()-2); break; } case 2 :{ Vector ctri5(5); if( (new_node_per_elem_array(nel,2)==1) && (new_node_per_elem_array(nel,4)==2) ){ ctri5(0) = connec_type_tmpl(nel,1); ctri5(1) = connec_type_tmpl(nel,2); ctri5(2) = connec_type_tmpl(nel,0); ctri5(3) = new_node_per_elem_array(nel,3); ctri5(4) = new_node_per_elem_array(nel,1); } else if((new_node_per_elem_array(nel,2)==1) && (new_node_per_elem_array(nel,4)==3)){ ctri5(0) = connec_type_tmpl(nel,0); ctri5(1) = connec_type_tmpl(nel,1); ctri5(2) = connec_type_tmpl(nel,2); ctri5(3) = new_node_per_elem_array(nel,1); ctri5(4) = new_node_per_elem_array(nel,3); } else if((new_node_per_elem_array(nel,2)==2) && (new_node_per_elem_array(nel,4)==3)){ ctri5(0) = connec_type_tmpl(nel,2); ctri5(1) = connec_type_tmpl(nel,0); ctri5(2) = connec_type_tmpl(nel,1); ctri5(3) = new_node_per_elem_array(nel,3); ctri5(4) = new_node_per_elem_array(nel,1); } else if((new_node_per_elem_array(nel,2)==2) && (new_node_per_elem_array(nel,4)==1)){ ctri5(0) = connec_type_tmpl(nel,1); ctri5(1) = connec_type_tmpl(nel,2); ctri5(2) = connec_type_tmpl(nel,0); ctri5(3) = new_node_per_elem_array(nel,1); ctri5(4) = new_node_per_elem_array(nel,3); } else if((new_node_per_elem_array(nel,2)==3) && (new_node_per_elem_array(nel,4)==1)){ ctri5(0) = connec_type_tmpl(nel,0); ctri5(1) = connec_type_tmpl(nel,1); ctri5(2) = connec_type_tmpl(nel,2); ctri5(3) = new_node_per_elem_array(nel,3); ctri5(4) = new_node_per_elem_array(nel,1); } else if((new_node_per_elem_array(nel,2)==3) && (new_node_per_elem_array(nel,4)==2)){ ctri5(0) = connec_type_tmpl(nel,2); ctri5(1) = connec_type_tmpl(nel,0); ctri5(2) = connec_type_tmpl(nel,1); ctri5(3) = new_node_per_elem_array(nel,1); ctri5(4) = new_node_per_elem_array(nel,3); } else{ AKANTU_DEBUG_ERROR("A triangle have only 3 segments not "<< new_node_per_elem_array(nel,2) << "and" << new_node_per_elem_array(nel,4)); } UInt nb_tri5 = connec_igfem_tri5.getSize(); connec_igfem_tri5.push_back(ctri5); element_tri5.element = nb_tri5; //new_elements.getList().push_back(element_tri5); new_elements_list(n_new_el) = element_tri5; if(type == _igfem_triangle_5){ new_numbering.push_back(new_nb_type_tmpl); ++new_nb_type_tmpl; } break; } default: AKANTU_DEBUG_ERROR("Igfem cannot add "<< new_node_per_elem_array(nel,0) << " nodes to triangles"); break; } element_type_tmpl.element = nel; old_elements_list(n_new_el) = element_type_tmpl; remove_elem.getList().push_back(element_type_tmpl); new_numbering(nel) = UInt(-1); ++n_new_el; } else{ new_numbering(nel) = new_nb_type_tmpl; ++new_nb_type_tmpl; } } // for(UInt nel=new_node_per_elem_array.getSize(); nel < this->mesh.getNbElement(type); ++nel) { // new_numbering(nel) = new_nb_type_tmpl; // ++new_nb_type_tmpl; // } UInt el_index = 0; for (UInt e = 0; e < this->mesh.getNbElement(type, ghost_type); ++e) { if (new_numbering(e) != UInt(-1)) { new_numbering(e) = el_index; ++el_index; } } total_new_elements += n_new_el; } if(total_new_elements > 0){ this->mesh.sendEvent(new_elements); this->mesh.sendEvent(remove_elem); } AKANTU_DEBUG_OUT(); } #else template void MeshSphereIntersector::buildResultFromQueryList(const std::list & query_list) { AKANTU_DEBUG_TO_IMPLEMENT(); } #endif __END_AKANTU__ #endif // __AKANTU_MESH_SPHERE_INTERSECTOR_TMPL_HH__