8 #include <unordered_set>
9 #include <boost/functional/hash.hpp>
27 template<
unsigned int dimA,
unsigned int dimB>
44 template<
unsigned int dim>
50 template<
unsigned int dim>
56 template<
unsigned int dim>
60 last_slave_for_3D_elements.assign(mesh->n_elements(), undefined_elm_idx_);
61 closed_elements.assign(mesh->n_elements(),
false);
67 template<
unsigned int dim>
71 if(elements_bb.size() == 0){
72 elements_bb.resize(mesh->n_elements());
73 bool first_3d_element =
true;
74 for (
auto elm : mesh->elements_range()) {
76 elements_bb[elm.idx()] = elm.bounding_box();
80 first_3d_element =
false;
81 mesh_3D_bb = elements_bb[elm.idx()];
83 mesh_3D_bb.expand(elements_bb[elm.idx()].min());
84 mesh_3D_bb.expand(elements_bb[elm.idx()].max());
92 template<
unsigned int dim>
96 unsigned int component_ele_idx = comp_ele.
idx(),
97 bulk_ele_idx = bulk_ele.
idx();
106 last_slave_for_3D_elements[bulk_ele_idx] = component_ele_idx;
109 intersection_list_[component_ele_idx].push_back(is);
119 template<
unsigned int dim>
127 for(; i < intersection_list_[component_ele_idx].size();i++){
129 if(intersection_list_[component_ele_idx][i].bulk_ele_idx() == bulk_ele_idx){
141 template<
unsigned int dim>
150 for (
auto elm : mesh->elements_range()) {
151 unsigned int component_ele_idx = elm.idx();
153 if (elm->dim() == dim &&
154 !closed_elements[component_ele_idx] &&
168 unsigned int bulk_ele_idx = *
it;
175 if (ele_3D->
dim() == 3 &&
176 (last_slave_for_3D_elements[bulk_ele_idx] != component_ele_idx &&
177 !intersection_exists(component_ele_idx,bulk_ele_idx) )
200 bool found = compute_initial_CI(elm, ele_3D);
203 unsigned int current_component_element_idx = component_ele_idx;
207 prolongation_decide(elm, ele_3D, intersection_list_[component_ele_idx].back());
212 bool element_covered =
true;
214 while(!bulk_queue_.empty()){
215 Prolongation pr = bulk_queue_.front();
218 if( pr.elm_3D_idx == undefined_elm_idx_)
221 element_covered =
false;
228 if(! closed_elements[current_component_element_idx])
229 closed_elements[current_component_element_idx] = element_covered;
232 if(!component_queue_.empty()){
233 Prolongation pr = component_queue_.front();
236 current_component_element_idx = pr.component_elm_idx;
240 component_queue_.pop();
243 while( !(component_queue_.empty() && bulk_queue_.empty()) );
247 if(closed_elements[component_ele_idx])
252 END_TIMER(
"Bounding box element iteration");
258 MessageOut().fmt(
"{}D-3D: number of intersections = {}\n", dim, n_intersections_);
265 template<
unsigned int dim>
268 DebugOut() <<
"######### ALGORITHM: compute_intersections_BIHtree #########\n";
274 for (
auto elm : mesh->elements_range()) {
275 unsigned int component_ele_idx = elm.idx();
277 if (elm.dim() == dim &&
291 unsigned int bulk_ele_idx = *
it;
294 if (ele_3D.
dim() == 3
304 if(is.
points().size() > 0) {
306 intersection_list_[component_ele_idx].push_back(is);
309 closed_elements[component_ele_idx] =
true;
313 END_TIMER(
"Bounding box element iteration");
320 template<
unsigned int dim>
325 compute_bounding_boxes();
330 for (
auto elm : mesh->elements_range()) {
331 unsigned int component_ele_idx = elm.idx();
333 if (elm.dim() == dim &&
334 !closed_elements[component_ele_idx] &&
335 elements_bb[component_ele_idx].intersect(mesh_3D_bb))
341 for (
auto ele_3D : mesh->elements_range()) {
342 unsigned int bulk_ele_idx = ele_3D.
idx();
349 if (ele_3D.
dim() == 3 &&
350 (last_slave_for_3D_elements[bulk_ele_idx] != component_ele_idx &&
351 elements_bb[component_ele_idx].intersect(elements_bb[bulk_ele_idx]) &&
352 !intersection_exists(component_ele_idx,bulk_ele_idx) )
376 bool found = compute_initial_CI(elm, ele_3D);
379 unsigned int current_component_element_idx = component_ele_idx;
384 prolongation_decide(elm, ele_3D, intersection_list_[component_ele_idx].back());
389 bool element_covered =
true;
391 while(!bulk_queue_.empty()){
398 element_covered =
false;
405 if(! closed_elements[current_component_element_idx])
406 closed_elements[current_component_element_idx] = element_covered;
409 if(!component_queue_.empty()){
417 component_queue_.pop();
420 while( !(component_queue_.empty() && bulk_queue_.empty()) );
424 if(closed_elements[component_ele_idx])
430 END_TIMER(
"Bounding box element iteration");
444 template<
unsigned int dim>
445 template<
unsigned int ele_dim>
448 unsigned int ip_obj_idx)
451 edges.reserve(ele_dim - ip_dim);
457 case 0:
if(ele_dim == 1) {
458 edges.push_back(mesh->edge(ele->
edge_idx(ip_obj_idx)));
462 for(
unsigned int j=0; j < RefElement<ele_dim>::n_sides_per_node; j++){
464 edges.push_back(mesh->edge(ele->
edge_idx(local_edge)));
470 case 1:
if(ele_dim == 2) {
471 edges.push_back(mesh->edge(ele->
edge_idx(ip_obj_idx)));
476 for(
unsigned int j=0; j < RefElement<ele_dim>::n_sides_per_line; j++){
478 edges.push_back(mesh->edge(ele->
edge_idx(local_edge)));
484 case 2:
ASSERT(ele_dim == 3);
485 edges.push_back(mesh->edge(ele->
edge_idx(ip_obj_idx)));
493 elements_idx.reserve(2*edges.size());
494 for(
Edge edg : edges)
495 for(
uint j=0; j < edg.n_sides();j++) {
496 if ( edg.side(j)->element().idx() != ele.
idx() )
497 elements_idx.push_back(edg.side(j)->element().idx());
504 template<
unsigned int dim>
506 unsigned int component_ele_idx,
507 std::queue< Prolongation >& queue)
512 last_slave_for_3D_elements[bulk_ele_idx] = component_ele_idx;
518 intersection_list_[component_ele_idx].push_back(il_other);
520 Prolongation pr = {component_ele_idx, bulk_ele_idx, (
unsigned int)intersection_list_[component_ele_idx].size() - 1};
530 template<
unsigned int dim>
542 unsigned int n_ip_vertices = 0;
577 if(IP.dim_A() < dim) {
578 if(IP.dim_A() == 0) n_ip_vertices++;
589 unsigned int bulk_current = bulk_ele.
idx();
592 for(
unsigned int& comp_neighbor_idx : comp_neighbors) {
593 if(!intersection_exists(comp_neighbor_idx,bulk_current))
594 create_prolongation(bulk_current, comp_neighbor_idx, component_queue_);
610 unsigned int comp_current = comp_ele.
idx();
611 unsigned int n_prolongations = 0;
613 for(
unsigned int& bulk_neighbor_idx : bulk_neighbors)
615 if(last_slave_for_3D_elements[bulk_neighbor_idx] == undefined_elm_idx_ ||
616 (last_slave_for_3D_elements[bulk_neighbor_idx] != comp_current &&
617 !intersection_exists(comp_current,bulk_neighbor_idx)))
618 n_prolongations += create_prolongation(bulk_neighbor_idx,
625 if(n_prolongations == 0)
627 Prolongation pr = {comp_ele.
idx(), undefined_elm_idx_, undefined_elm_idx_};
628 bulk_queue_.push(pr);
634 if(n_ip_vertices == is.
size()) closed_elements[comp_ele.
idx()] =
true;
647 ASSERT(0).add_value(bulk_ele_idx,
"bulk_ele_idx").error(
"Want to add the same intersection!");
653 template<
unsigned int dim>
685 prolongation_decide(elm, ele_3D,is);
715 unsigned int ele_idx, eleA_idx, eleB_idx,
716 componentA_idx, componentB_idx,
719 typedef std::pair<unsigned int, unsigned int> ipair;
720 std::unordered_set<ipair, boost::hash<ipair>> computed_pairs;
727 if(intersection_map[ele_idx].size() < 2)
continue;
732 for(
unsigned int i=0; i < local_map.size(); i++)
737 eleA_idx = local_map[i].first;
739 if(eleA->
dim() !=2 )
continue;
745 for(
unsigned int j=i+1; j < local_map.size(); j++)
747 eleB_idx = local_map[j].first;
750 if(componentA_idx == componentB_idx)
continue;
754 if(eleB->
dim() !=2 )
continue;
758 temp_eleA_idx = eleA_idx;
760 if (componentA_idx < componentB_idx){
766 ipair ip = std::make_pair(temp_eleA_idx, eleB_idx);
767 if(computed_pairs.count(ip) == 1){
773 computed_pairs.emplace(ip);
800 MessageOut() <<
"2D-2D: number of intersections = " << storage.size() <<
"\n";
816 unsigned int n_local_intersection = CI.
compute(is);
819 if(n_local_intersection > 1){
831 std::queue<unsigned int> queue;
834 if (ele->dim() == 2 &&
838 queue.push(ele.idx());
840 while(!queue.empty()){
841 unsigned int ele_idx = queue.front();
844 for(
unsigned int sid=0; sid < elm->
n_sides(); sid++) {
851 queue.push(neigh_idx);
889 unsigned int ele_idx = ele.idx();
891 if(intersection_map[ele_idx].size() < 2)
continue;
896 for(
unsigned int i=0; i < local_map.size(); i++)
901 unsigned int eleA_idx = local_map[i].first;
904 if(eleA.
dim() !=1 )
continue;
906 for(
unsigned int j=0; j < local_map.size(); j++)
908 unsigned int eleB_idx = local_map[j].first;
910 if(eleB.
dim() !=2 )
continue;
914 for(
unsigned int i=0; i<intersection_map[eleA_idx].size(); i++)
916 if(intersection_map[eleA_idx][i].first == eleB_idx) {
933 if(n_local_intersection > 0)
936 intersection_map[eleA_idx].push_back(std::make_pair(
940 intersection_map[eleB_idx].push_back(std::make_pair(
952 MessageOut() <<
"1D-2D [3]: number of intersections = " << storage.size() <<
"\n";
996 unsigned int component_ele_idx = elm.idx();
1012 unsigned int bulk_ele_idx = *
it;
1015 if (ele_2D.
dim() == 2) {
1023 if(is.
points().size() > 0) {
1028 END_TIMER(
"Bounding box element iteration");
1042 unsigned int component_ele_idx = elm.idx();
1044 if (elm->dim() == 1)
1053 for(
unsigned int bulk_ele_idx : candidate_list) {
1056 if (ele_2D->
dim() == 2) {
1064 if(is.
points().size() > 0) {
1069 END_TIMER(
"Bounding box element iteration");
#define ASSERT_PERMANENT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
Class for O(log N) lookup for intersections with a set of bounding boxes.
void find_bounding_box(const BoundingBox &boundingBox, std::vector< unsigned int > &result_list, bool full_list=false) const
const BoundingBox & ele_bounding_box(unsigned int ele_idx) const
Gets bounding box of element of given index ele_index.
const BoundingBox & tree_box() const
bool intersect(const BoundingBox &b2) const
Class for 1D-2D intersections.
unsigned int compute_final_in_plane(std::vector< IPAux12 > &IP12s)
unsigned int compute_final(std::vector< IPAux12 > &IP12s)
Class for 2D-2D intersections.
void init()
Initializes lower dimensional objects. Sets correctly the pointers to Plucker coordinates and product...
unsigned int compute(IntersectionAux< 2, 2 > &intersection)
Computes final 2D-2D intersection. Computes CIs (side vs triangle) in following order: [A0_B,...
unsigned int n_sides() const
Returns number of sides aligned with the edge.
SideIter side(const unsigned int i) const
Gets side iterator of the i -th side.
SideIter side(const unsigned int loc_index)
unsigned int idx() const
We need this method after replacing Region by RegionIdx, and movinf RegionDB instance into particular...
unsigned int edge_idx(unsigned int edg_idx) const
Return edge_idx of given index.
unsigned int n_neighs_vb() const
Return number of neighbours.
unsigned int n_sides() const
void compute_intersections_3(std::vector< std::vector< ILpair >> &intersection_map, std::vector< IntersectionLocal< 1, 2 >> &storage)
Runs the algorithm (3): compute 1D-2D intersection inside 3D mesh. It directly fills the intersection...
void compute_intersections_1(const BIHTree &bih)
Runs the algorithm (1): compute 1D-2D intersection in 2D plane. BIH is used to find intersection cand...
std::vector< IntersectionAux< 1, 2 > > intersectionaux_storage12_
Stores temporarily 1D-2D intersections.
InspectElementsAlgorithm12(Mesh *input_mesh)
void compute_intersections_2(const BIHTree &bih)
Runs the algorithm (2): compute 1D-2D intersection in 3D ambient space BIH is used to find intersecti...
std::vector< unsigned int > component_idx_
InspectElementsAlgorithm22(Mesh *input_mesh)
void compute_single_intersection(const ElementAccessor< 3 > &eleA, const ElementAccessor< 3 > &eleB, std::vector< IntersectionLocal< 2, 2 >> &storage)
Computes fundamental intersection of two 2D elements.
void create_component_numbering()
Creates numbering of the 2D components and fills component_idx_ vector.
unsigned int component_counter_
const unsigned int unset_comp
void compute_intersections(std::vector< std::vector< ILpair >> &intersection_map, std::vector< IntersectionLocal< 2, 2 >> &storage)
Runs the core algorithm for computing 2D-2D intersection in 3D.
void compute_intersections_BB()
bool intersection_exists(unsigned int component_ele_idx, unsigned int bulk_ele_idx)
A hard way to find whether the intersection of two elements has already been computed,...
unsigned int create_prolongation(unsigned int bulk_ele_idx, unsigned int component_ele_idx, std::queue< Prolongation > &queue)
~InspectElementsAlgorithm()
InspectElementsAlgorithm(Mesh *_mesh)
void assert_same_intersection(unsigned int comp_ele_idx, unsigned int bulk_ele_idx)
void prolongate(const Prolongation &pr)
Computes the intersection for a candidate in a queue and calls prolongation_decide again.
std::vector< unsigned int > get_element_neighbors(const ElementAccessor< 3 > &ele, unsigned int ip_dim, unsigned int ip_obj_idx)
void compute_intersections(const BIHTree &bih)
Uses BIHtree to find the initial candidate of a component and then prolongates the component interset...
void prolongation_decide(const ElementAccessor< 3 > &comp_ele, const ElementAccessor< 3 > &bulk_ele, IntersectionAux< dim, 3 > is)
void compute_bounding_boxes()
Computes bounding boxes of all elements. Fills elements_bb and mesh_3D_bb.
bool compute_initial_CI(const ElementAccessor< 3 > &comp_ele, const ElementAccessor< 3 > &bulk_ele)
Computes the first intersection, from which we then prolongate.
std::vector< std::vector< IntersectionAux< dim, 3 > > > intersection_list_
Resulting vector of intersections.
void compute_intersections_BIHtree(const BIHTree &bih)
Uses only BIHtree to find intersection candidates. (No prolongation).
Mesh * mesh
Auxiliary function that translates ElementAccessor<3> to Simplex<simplex_dim>.
IntersectionAlgorithmBase(Mesh *mesh)
Internal auxiliary class representing intersection object of simplex<dimA> and simplex<dimB>.
unsigned int ips_in_face() const
Returns idx of face when all IPs lie on it; -1 otherwise.
unsigned int size() const
Returns number of intersection points.
std::vector< IntersectionPointAux< dimA, dimB > > & points()
Returns intersection points by a reference.
void set_duplicities(unsigned int n_duplicities)
Class represents intersection of two elements.
Internal auxiliary class represents an intersection point of simplex<N> and simplex<M>.
unsigned int n_elements() const
Range< ElementAccessor< 3 > > elements_range() const
Returns range of mesh elements.
ElementAccessor< 3 > element_accessor(unsigned int idx) const
Create and return ElementAccessor to element of given idx.
static const IdxVector<(InDim >OutDim ? InDim+1 :dim-InDim) > interact(TInteraction< OutDim, InDim > interaction, bool inv=false)
ElementAccessor< 3 > element() const
Returns iterator to the element of the side.
Edge edge() const
Returns pointer to the edge connected to the side.
Fundamental simplicial intersections.
Internal class representing intersection object.
Classes with algorithms for computation of intersections of meshes.
Internal class representing intersection point.
#define MessageOut()
Macro defining 'message' record of log.
#define DebugOut()
Macro defining 'debug' record of log.
void swap(nlohmann::json &j1, nlohmann::json &j2) noexcept(is_nothrow_move_constructible< nlohmann::json >::value and is_nothrow_move_assignable< nlohmann::json >::value)
exchanges the values of two JSON objects
Implementation of range helper class.
Class RefElement defines numbering of vertices, sides, calculation of normal vectors etc.
Auxiliary structure for prolongation process.
unsigned int dictionary_idx
unsigned int component_elm_idx
#define END_TIMER(tag)
Ends a timer with specified tag.
#define ADD_CALLS(n_calls)
Increase number of calls in actual timer.
#define START_TIMER(tag)
Starts a timer with specified tag.