star-enclosures-2d

Extract enclosures from 2D geometry
git clone git://git.meso-star.fr/star-enclosures-2d.git
Log | Files | Refs | README | LICENSE

commit 0ed0c75ff3ff34f7002e87b5485e5c5f65b92614
parent 3802be212827192b310bb14dc5d2a3a7242ff318
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date:   Fri,  6 Apr 2018 14:12:51 +0200

Improve perf. by approx. 20% by easing omp thread syncs.

Diffstat:
Msrc/senc2d_scene_analyze.c | 1303++++++++++++++++++++++++++++++++++++++++++-------------------------------------
Asrc/senc2d_scene_analyze.c.save | 1254+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Msrc/senc2d_scene_analyze_c.h | 23+++++++++++++++--------
3 files changed, 1970 insertions(+), 610 deletions(-)

diff --git a/src/senc2d_scene_analyze.c b/src/senc2d_scene_analyze.c @@ -188,7 +188,7 @@ self_hit_filter hit_side = (hit->normal[1] > 0) ? SIDE_FRONT : SIDE_BACK; hit_component = hit_seg_comp->component[hit_side]; - /* Not self hit or distance should be small */ + /* If self hit, distance should be small */ ASSERT(hit_component != *origin_component || hit->distance < 1e-6); return (hit_component == *origin_component); } @@ -202,297 +202,300 @@ extract_connex_components struct darray_segment_comp* segments_comp, struct s2d_scene_view** s2d_view) { + /* This function is called from an omp parallel block and executed + * concurrently. */ res_T res = RES_OK; const struct senc2d_scene* scn; struct mem_allocator* alloc; - ATOMIC component_count = 0; - volatile int exit_for = 0; int64_t mm; + struct darray_side_id stack; #ifndef NDEBUG seg_id_t s_; component_id_t c; #endif + /* shared between threads */ + static ATOMIC component_count = 0; + static volatile int exit_for = 0; ASSERT(segsides && desc && connex_components && segments_tmp_array); - ASSERT(darray_ptr_component_descriptor_size_get(connex_components) == 0); alloc = descriptor_get_allocator(desc); scn = desc->scene; - /* Just a hint; to avoid contention on first loop */ - OK2(darray_ptr_component_descriptor_reserve(connex_components, 2 * scn->nmeds), - error_); /* Cannot goto into openmp block */ - #ifndef NDEBUG - FOR_EACH(s_, 0, scn->nusegs) { - const struct segment_in* seg_in = - darray_segment_in_cdata_get(&scn->segments_in) + s_; - const struct side_range* media_use - = darray_side_range_cdata_get(&scn->media_use); - FOR_EACH(mm, 0, 2) { - const side_id_t side = SEGIDxSIDE_2_SEGSIDE(s_, mm); - const medium_id_t medium = seg_in->medium[mm]; - ASSERT(media_use[medium].first <= side && side <= media_use[medium].last); + #pragma omp single + { + ASSERT(darray_ptr_component_descriptor_size_get(connex_components) == 0); + FOR_EACH(s_, 0, scn->nusegs) { + const struct segment_in* seg_in = + darray_segment_in_cdata_get(&scn->segments_in) + s_; + const struct side_range* media_use + = darray_side_range_cdata_get(&scn->media_use); + FOR_EACH(mm, 0, 2) { + const side_id_t side = SEGIDxSIDE_2_SEGSIDE(s_, mm); + const medium_id_t medium = seg_in->medium[mm]; + ASSERT(media_use[medium].first <= side && side <= media_use[medium].last); + } } - } + } /* Implicit barrier here */ #endif - #pragma omp parallel - { - struct darray_side_id stack; - darray_side_id_init(alloc, &stack); + darray_side_id_init(alloc, &stack); - #pragma omp for schedule(dynamic) nowait - for(mm = 0; mm < (int64_t)scn->nmeds; mm++) { /* Process all media */ - const medium_id_t m = (medium_id_t)mm; - struct cc_descriptor* cc; - /* Any not-already-used side is used as a starting point */ - side_id_t first_side_not_in_component; + #pragma omp for schedule(dynamic) nowait + for(mm = 0; mm < (int64_t)scn->nmeds; mm++) { /* Process all media */ + const medium_id_t m = (medium_id_t)mm; + struct cc_descriptor* cc; + /* Any not-already-used side is used as a starting point */ + side_id_t first_side_not_in_component; - if(exit_for) continue; - first_side_not_in_component - = darray_side_range_cdata_get(&scn->media_use)[m].first; - if(first_side_not_in_component == SIDE_NULL__) - continue; /* Unused medium */ - ASSERT(first_side_not_in_component < 2 * scn->nusegs); - ASSERT(darray_side_id_size_get(&stack) == 0); - for(;;) { /* Process all components for this medium */ - const side_id_t start_side_id = get_side_not_in_connex_component(scn, - segsides, &first_side_not_in_component, m); - side_id_t crt_side_id = start_side_id; - side_id_t last_side_id = start_side_id; - ASSERT(start_side_id == SIDE_NULL__ || start_side_id < 2 * scn->nusegs); - if(start_side_id == SIDE_NULL__) - break; /* start_side_id=SIDE_NULL__ => done! */ - ASSERT(segsides[start_side_id].list_id == FLAG_LIST_SIDE_LIST); + if(exit_for) continue; + first_side_not_in_component + = darray_side_range_cdata_get(&scn->media_use)[m].first; + if(first_side_not_in_component == SIDE_NULL__) + continue; /* Unused medium */ + ASSERT(first_side_not_in_component < 2 * scn->nusegs); + ASSERT(darray_side_id_size_get(&stack) == 0); + for(;;) { /* Process all components for this medium */ + const side_id_t start_side_id = get_side_not_in_connex_component(scn, + segsides, &first_side_not_in_component, m); + side_id_t crt_side_id = start_side_id; + side_id_t last_side_id = start_side_id; + ASSERT(start_side_id == SIDE_NULL__ || start_side_id < 2 * scn->nusegs); + if(start_side_id == SIDE_NULL__) + break; /* start_side_id=SIDE_NULL__ => done! */ + ASSERT(segsides[start_side_id].list_id == FLAG_LIST_SIDE_LIST); #ifndef NDEBUG - { - seg_id_t tid = SEGSIDE_2_SEG(start_side_id); - enum side_flag s = SEGSIDE_2_SIDE(start_side_id); - medium_id_t side_med - = darray_segment_in_data_get(&desc->scene->segments_in)[tid].medium[s]; - ASSERT(side_med == m); - } + { + seg_id_t tid = SEGSIDE_2_SEG(start_side_id); + enum side_flag s = SEGSIDE_2_SIDE(start_side_id); + medium_id_t side_med + = darray_segment_in_data_get(&desc->scene->segments_in)[tid].medium[s]; + ASSERT(side_med == m); + } #endif - /* Create and init a new component */ - cc = MEM_ALLOC(alloc, sizeof(struct cc_descriptor)); - if(!cc) { - res = RES_MEM_ERR; - goto error1; - } - cc_descriptor_init(alloc, cc); - ASSERT(m == segsides[start_side_id].medium); - cc->cc_id = (component_id_t)(ATOMIC_INCR(&component_count) - 1); - cc->medium = m; - cc->side_range.first = start_side_id; - - for(;;) { /* Process all sides of this component */ - int i; - enum side_flag crt_side_flag = SEGSIDE_2_SIDE(crt_side_id); - struct segside* crt_side = segsides + crt_side_id; - const seg_id_t crt_seg_id = SEGSIDE_2_SEG(crt_side_id); - const struct segment_in* seg_in = - darray_segment_in_cdata_get(&scn->segments_in) + crt_seg_id; - struct segment_comp* seg_comp = - darray_segment_comp_data_get(segments_comp) + crt_seg_id; - const struct segment_tmp* const seg_tmp = - darray_segment_tmp_cdata_get(segments_tmp_array) + crt_seg_id; - const union double2* vertices = - darray_position_cdata_get(&scn->vertices); - ASSERT(crt_seg_id < scn->nusegs); - ASSERT(segsides[crt_side_id].medium == m); - - /* Record Ymax information - * Keep track of the appropriate vertex/side of the connex component - * in order to cast a ray at the component grouping step of the - * algorithm. - * The most appropriate vertex is (the) one with the greater Y - * coordinate. - * If more than one vertex/side has the same Y, we want the side that - * most faces Y (that is the one with the greater ny). - * This is mandatory to select the correct side when both sides of a - * segment are candidate. */ - if(cc->max_vrtx[1] <= seg_tmp->max_y) { - /* Can either improve y or ny */ - - if(cc->max_y_side_id == SEGSIDE_OPPOSITE(crt_side_id)) { - /* Both sides are in cc and the opposite side is currently selected - * Just keep the side with ny>0 */ - if(cc->max_y_ny < 0) { - /* Change side! */ - cc->max_y_ny = -cc->max_y_ny; - cc->max_y_side_id = crt_side_id; - } + /* Create and init a new component */ + cc = MEM_ALLOC(alloc, sizeof(struct cc_descriptor)); + if(!cc) { + res = RES_MEM_ERR; + goto error1; + } + cc_descriptor_init(alloc, cc); + ASSERT(m == segsides[start_side_id].medium); + cc->cc_id = (component_id_t)(ATOMIC_INCR(&component_count) - 1); + cc->medium = m; + cc->side_range.first = start_side_id; + + for(;;) { /* Process all sides of this component */ + int i; + enum side_flag crt_side_flag = SEGSIDE_2_SIDE(crt_side_id); + struct segside* crt_side = segsides + crt_side_id; + const seg_id_t crt_seg_id = SEGSIDE_2_SEG(crt_side_id); + const struct segment_in* seg_in = + darray_segment_in_cdata_get(&scn->segments_in) + crt_seg_id; + struct segment_comp* seg_comp = + darray_segment_comp_data_get(segments_comp) + crt_seg_id; + const struct segment_tmp* const seg_tmp = + darray_segment_tmp_cdata_get(segments_tmp_array) + crt_seg_id; + const union double2* vertices = + darray_position_cdata_get(&scn->vertices); + ASSERT(crt_seg_id < scn->nusegs); + ASSERT(segsides[crt_side_id].medium == m); + + /* Record Ymax information + * Keep track of the appropriate vertex/side of the connex component + * in order to cast a ray at the component grouping step of the + * algorithm. + * The most appropriate vertex is (the) one with the greater Y + * coordinate. + * If more than one vertex/side has the same Y, we want the side that + * most faces Y (that is the one with the greater ny). + * This is mandatory to select the correct side when both sides of a + * segment are candidate. */ + if(cc->max_vrtx[1] <= seg_tmp->max_y) { + /* Can either improve y or ny */ + + if(cc->max_y_side_id == SEGSIDE_OPPOSITE(crt_side_id)) { + /* Both sides are in cc and the opposite side is currently selected + * Just keep the side with ny>0 */ + if(cc->max_y_ny < 0) { + /* Change side! */ + cc->max_y_ny = -cc->max_y_ny; + cc->max_y_side_id = crt_side_id; + } + } else { + double edge[2], normal[2], norm, side_ny; + int process = 0; + + d2_sub(edge, vertices[seg_in->vertice_id[1]].vec, + vertices[seg_in->vertice_id[0]].vec); + d2(normal, edge[1], -edge[0]); + norm = d2_normalize(normal, normal); + ASSERT(norm); (void) norm; + + /* Geometrical normal points toward the right */ + if(SEGSIDE_IS_FRONT(crt_side_id)) { + side_ny = -normal[1]; + process = 1; } else { - double edge[2], normal[2], norm, side_ny; - int process = 0; - - d2_sub(edge, vertices[seg_in->vertice_id[1]].vec, - vertices[seg_in->vertice_id[0]].vec); - d2(normal, edge[1], -edge[0]); - norm = d2_normalize(normal, normal); - ASSERT(norm); (void) norm; - - /* Geometrical normal points toward the right */ - if(SEGSIDE_IS_FRONT(crt_side_id)) { - side_ny = -normal[1]; - process = 1; - } else { - side_ny = normal[1]; - process = 1; + side_ny = normal[1]; + process = 1; + } + + if(process) { + int change = 0; + if(cc->max_vrtx[1] < seg_tmp->max_y) { + change = 1; /* Try first to improve y */ + } else if(cc->max_y_ny <= 0 && fabs(cc->max_y_ny) < fabs(side_ny)) { + change = 1; /* If ny <= 0, the more negative the better */ + } else if(cc->max_y_ny > 0 && cc->max_y_ny < side_ny) { + change = 1; /* If ny > 0, the more positive the better */ } - if(process) { - int change = 0; - if(cc->max_vrtx[1] < seg_tmp->max_y) { - change = 1; /* Try first to improve y */ - } else if(cc->max_y_ny <= 0 && fabs(cc->max_y_ny) < fabs(side_ny)) { - change = 1; /* If ny <= 0, the more negative the better */ - } else if(cc->max_y_ny > 0 && cc->max_y_ny < side_ny) { - change = 1; /* If ny > 0, the more positive the better */ - } - - if(change) { - cc->max_y_ny = side_ny; - cc->max_y_side_id = crt_side_id; - ASSERT(seg_tmp->max_y_vrtx_rank < 2); - ASSERT(seg_in->vertice_id[seg_tmp->max_y_vrtx_rank] < scn->nverts); - cc->max_y_vrtx_id = seg_in->vertice_id[seg_tmp->max_y_vrtx_rank]; - ASSERT(seg_tmp->max_y == vertices[cc->max_y_vrtx_id].pos.y); - d2_set(cc->max_vrtx, vertices[cc->max_y_vrtx_id].vec); - } + if(change) { + cc->max_y_ny = side_ny; + cc->max_y_side_id = crt_side_id; + ASSERT(seg_tmp->max_y_vrtx_rank < 2); + ASSERT(seg_in->vertice_id[seg_tmp->max_y_vrtx_rank] < scn->nverts); + cc->max_y_vrtx_id = seg_in->vertice_id[seg_tmp->max_y_vrtx_rank]; + ASSERT(seg_tmp->max_y == vertices[cc->max_y_vrtx_id].pos.y); + d2_set(cc->max_vrtx, vertices[cc->max_y_vrtx_id].vec); } } } + } - /* Record crt_side both as component and segment level */ - cc->side_count++; - segsides[crt_side_id].list_id = FLAG_LIST_COMPONENT; - ASSERT(seg_comp->component[crt_side_flag] == COMPONENT_NULL__); - seg_comp->component[crt_side_flag] = cc->cc_id; + /* Record crt_side both as component and segment level */ + cc->side_count++; + segsides[crt_side_id].list_id = FLAG_LIST_COMPONENT; + ASSERT(seg_comp->component[crt_side_flag] == COMPONENT_NULL__); + seg_comp->component[crt_side_flag] = cc->cc_id; #ifndef NDEBUG - crt_side->member_of_cc = cc->cc_id; + crt_side->member_of_cc = cc->cc_id; #endif - /* Store neighbour sides in a waiting stack */ - FOR_EACH(i, 0, 2) { - side_id_t neighbour_id = crt_side->facing_side_id[i]; - seg_id_t nbour_seg_id = SEGSIDE_2_SEG(neighbour_id); - const struct segside* neighbour = segsides + neighbour_id; - ASSERT(m == crt_side->medium); - if(neighbour->medium != crt_side->medium) { - /* Found medium discontinuity! Model topology is broken. */ - const struct segment_in* segments_in - = darray_segment_in_cdata_get(&scn->segments_in); - const union double2* positions - = darray_position_cdata_get(&scn->vertices); - log_err(scn->dev, - "Medium mismatch found between neighbour segments %lu %s" - " side and %u %s side.\n", - (unsigned long) segments_in[crt_seg_id].global_id, - SEGSIDE_IS_FRONT(crt_side_id) ? "front" : "back", - segments_in[nbour_seg_id].global_id, - SEGSIDE_IS_FRONT(neighbour_id) ? "front" : "back"); - log_err(scn->dev, - "Segment %lu:\n (%g %g) (%g %g))\n", - (unsigned long) segments_in[crt_seg_id].global_id, - SPLIT2(positions[segments_in[crt_seg_id].vertice_id[0]].vec), - SPLIT2(positions[segments_in[crt_seg_id].vertice_id[1]].vec)); - log_err(scn->dev, - "Segment %lu:\n (%g %g) (%g %g)\n", - (unsigned long) segments_in[nbour_seg_id].global_id, - SPLIT2(positions[segments_in[nbour_seg_id].vertice_id[0]].vec), - SPLIT2(positions[segments_in[nbour_seg_id].vertice_id[1]].vec)); - log_err(desc->scene->dev, "Media: %lu VS %lu\n", - (unsigned long)neighbour->medium, (unsigned long)crt_side->medium); - res = RES_BAD_ARG; - goto error1; - } - if(neighbour->list_id == FLAG_LIST_COMPONENT) { - /* Already processed */ + /* Store neighbour sides in a waiting stack */ + FOR_EACH(i, 0, 2) { + side_id_t neighbour_id = crt_side->facing_side_id[i]; + seg_id_t nbour_seg_id = SEGSIDE_2_SEG(neighbour_id); + const struct segside* neighbour = segsides + neighbour_id; + ASSERT(m == crt_side->medium); + if(neighbour->medium != crt_side->medium) { + /* Found medium discontinuity! Model topology is broken. */ + const struct segment_in* segments_in + = darray_segment_in_cdata_get(&scn->segments_in); + const union double2* positions + = darray_position_cdata_get(&scn->vertices); + log_err(scn->dev, + "Medium mismatch found between neighbour segments %lu %s" + " side and %u %s side.\n", + (unsigned long)segments_in[crt_seg_id].global_id, + SEGSIDE_IS_FRONT(crt_side_id) ? "front" : "back", + segments_in[nbour_seg_id].global_id, + SEGSIDE_IS_FRONT(neighbour_id) ? "front" : "back"); + log_err(scn->dev, + "Segment %lu:\n (%g %g) (%g %g))\n", + (unsigned long)segments_in[crt_seg_id].global_id, + SPLIT2(positions[segments_in[crt_seg_id].vertice_id[0]].vec), + SPLIT2(positions[segments_in[crt_seg_id].vertice_id[1]].vec)); + log_err(scn->dev, + "Segment %lu:\n (%g %g) (%g %g)\n", + (unsigned long)segments_in[nbour_seg_id].global_id, + SPLIT2(positions[segments_in[nbour_seg_id].vertice_id[0]].vec), + SPLIT2(positions[segments_in[nbour_seg_id].vertice_id[1]].vec)); + log_err(desc->scene->dev, "Media: %lu VS %lu\n", + (unsigned long)neighbour->medium, (unsigned long)crt_side->medium); + res = RES_BAD_ARG; + goto error1; + } + if(neighbour->list_id == FLAG_LIST_COMPONENT) { + /* Already processed */ #ifndef NDEBUG - ASSERT(neighbour->member_of_cc == cc->cc_id); + ASSERT(neighbour->member_of_cc == cc->cc_id); #endif - continue; - } - if(neighbour->list_id == FLAG_WAITING_STACK) { - continue; /* Already processed */ - } - add_side_to_stack(scn, &stack, segsides, neighbour_id); - } - if(darray_side_id_size_get(&stack) == 0) - break; /* Empty stack => connex component is done! */ - crt_side_id = get_side_from_stack(segsides, &stack); - last_side_id = MMAX(last_side_id, crt_side_id); - } - /* Keep track of this new connex component */ - cc->side_range.last = last_side_id; - /* Need to synchronize connex_components growth as this global structure - * is accessed by multipe threads */ - #pragma omp critical - { - struct cc_descriptor** components; - size_t sz = darray_ptr_component_descriptor_size_get(connex_components); - if(sz <= cc->cc_id) { - res_T tmp_res = darray_ptr_component_descriptor_resize(connex_components, - 1 + cc->cc_id); - if(tmp_res != RES_OK) res = tmp_res; + continue; } - if(res == RES_OK) { - /* Don't set the pointer before resize as this can lead to move data */ - components = - darray_ptr_component_descriptor_data_get(connex_components); - ASSERT(components[cc->cc_id] == NULL); - components[cc->cc_id] = cc; + if(neighbour->list_id == FLAG_WAITING_STACK) { + continue; /* Already processed */ } + add_side_to_stack(scn, &stack, segsides, neighbour_id); } - OK2(res, error1); + if(darray_side_id_size_get(&stack) == 0) + break; /* Empty stack => connex component is done! */ + crt_side_id = get_side_from_stack(segsides, &stack); + last_side_id = MMAX(last_side_id, crt_side_id); } - continue; - error1: - /* Cannot goto out of openmp block */ - exit_for = 1; - continue; - } - /* No barrier here (nowait clause). - * The first thread executes the single block */ - darray_side_id_release(&stack); - #pragma omp single nowait - if(res == RES_OK) { - struct s2d_device* s2d = NULL; - struct s2d_scene* s2d_scn = NULL; - struct s2d_shape* s2d_shp = NULL; - struct s2d_vertex_data attribs; - - attribs.type = S2D_FLOAT2; - attribs.usage = S2D_POSITION; - attribs.get = get_scn_position; - - /* Put geometry in a 2D view */ - OK(s2d_device_create(desc->scene->dev->logger, alloc, 0, &s2d)); - OK(s2d_scene_create(s2d, &s2d_scn)); - OK(s2d_shape_create_line_segments(s2d, &s2d_shp)); - - /* Back to API type for ntris and nverts */ - ASSERT(desc->scene->nusegs < UINT_MAX); - ASSERT(desc->scene->nuverts < UINT_MAX); - OK(s2d_line_segments_setup_indexed_vertices(s2d_shp, (unsigned)desc->scene->nusegs, - get_scn_indices, (unsigned)desc->scene->nuverts, &attribs, 1, desc->scene)); - s2d_line_segments_set_hit_filter_function(s2d_shp, self_hit_filter, segments_comp); - OK(s2d_scene_attach_shape(s2d_scn, s2d_shp)); - OK(s2d_scene_view_create(s2d_scn, S2D_TRACE, s2d_view)); - error: - if(s2d) S2D(device_ref_put(s2d)); - if(s2d_scn) S2D(scene_ref_put(s2d_scn)); - if(s2d_shp) S2D(shape_ref_put(s2d_shp)); + /* Keep track of this new connex component */ + cc->side_range.last = last_side_id; + /* Need to synchronize connex_components growth as this global structure + * is accessed by multipe threads */ + #pragma omp critical + { + struct cc_descriptor** components; + size_t sz = darray_ptr_component_descriptor_size_get(connex_components); + if(sz <= cc->cc_id) { + res_T tmp_res = darray_ptr_component_descriptor_resize(connex_components, + 1 + cc->cc_id); + if(tmp_res != RES_OK) res = tmp_res; + } + if(res == RES_OK) { + /* Don't set the pointer before resize as this can lead to move data */ + components = + darray_ptr_component_descriptor_data_get(connex_components); + ASSERT(components[cc->cc_id] == NULL); + components[cc->cc_id] = cc; + } + } + OK2(res, error1); } + continue; + error1: + /* Cannot goto out of openmp block */ + exit_for = 1; + continue; + } + /* No barrier here (nowait clause). + * The first thread here creates the s2d view */ + darray_side_id_release(&stack); + #pragma omp single nowait + if(res == RES_OK) { + struct s2d_device* s2d = NULL; + struct s2d_scene* s2d_scn = NULL; + struct s2d_shape* s2d_shp = NULL; + struct s2d_vertex_data attribs; + + attribs.type = S2D_FLOAT2; + attribs.usage = S2D_POSITION; + attribs.get = get_scn_position; + + /* Put geometry in a 2D view */ + OK(s2d_device_create(desc->scene->dev->logger, alloc, 0, &s2d)); + OK(s2d_scene_create(s2d, &s2d_scn)); + OK(s2d_shape_create_line_segments(s2d, &s2d_shp)); + + /* Back to API type for ntris and nverts */ + ASSERT(desc->scene->nusegs < UINT_MAX); + ASSERT(desc->scene->nuverts < UINT_MAX); + OK(s2d_line_segments_setup_indexed_vertices(s2d_shp, + (unsigned)desc->scene->nusegs, get_scn_indices, + (unsigned)desc->scene->nuverts, &attribs, 1, desc->scene)); + s2d_line_segments_set_hit_filter_function(s2d_shp, self_hit_filter, + segments_comp); + OK(s2d_scene_attach_shape(s2d_scn, s2d_shp)); + OK(s2d_scene_view_create(s2d_scn, S2D_TRACE, s2d_view)); + error: + if(s2d) S2D(device_ref_put(s2d)); + if(s2d_scn) S2D(scene_ref_put(s2d_scn)); + if(s2d_shp) S2D(shape_ref_put(s2d_shp)); } OK2(res, error_); +#ifndef NDEBUG + /* Need to wait for all threads done to be able to check stuff */ + #pragma omp barrier ASSERT(component_count == (int)darray_ptr_component_descriptor_size_get(connex_components)); -#ifndef NDEBUG FOR_EACH(s_, 0, scn->nusegs) { struct segment_comp* seg_comp = darray_segment_comp_data_get(segments_comp) + s_; @@ -505,14 +508,14 @@ extract_connex_components ASSERT(components[c] != NULL && components[c]->cc_id == c); } + ASSERT(desc->segment_count + == darray_segment_comp_size_get(segments_comp)); #endif exit: - ASSERT(desc->segment_count - == darray_segment_comp_size_get(segments_comp)); - /* segments_enc is still unused: no size to assert */ return res; error_: + exit_for = 1; goto exit; } @@ -524,16 +527,21 @@ group_connex_components struct darray_ptr_component_descriptor* connex_components, struct s2d_scene_view* s2d_view) { + /* This function is called from an omp parallel block and executed + * concurrently. */ res_T res = RES_OK; struct cc_descriptor** descriptors; size_t tmp; component_id_t cc_count; int64_t ccc; - volatile int exit_for = 0; - ATOMIC next_enclosure_id = desc->enclosures_count; - struct cc_descriptor* infinity_first_cc = NULL; + /* shared between threads */ + static struct cc_descriptor* infinity_first_cc = NULL; + static volatile int exit_for = 0; + static ATOMIC next_enclosure_id = 1; + (void)segsides; ASSERT(desc && segsides && segments_comp && connex_components); + ASSERT(desc->enclosures_count == 1); descriptors = darray_ptr_component_descriptor_data_get(connex_components); tmp = darray_ptr_component_descriptor_size_get(connex_components); @@ -541,7 +549,7 @@ group_connex_components cc_count = (component_id_t)tmp; /* Cast rays to find links between connex components */ - #pragma omp parallel for + #pragma omp for for(ccc = 0; ccc < (int64_t)cc_count; ccc++) { component_id_t c = (component_id_t)ccc; struct s2d_hit hit = S2D_HIT_NULL; @@ -599,18 +607,18 @@ group_connex_components log_err(desc->scene->dev, "Medium mismatch found between segment %lu %s side and segment" " %lu %s side, both facing infinity.\n", - (unsigned long) segments_in[t1].global_id, + (unsigned long)segments_in[t1].global_id, SEGSIDE_IS_FRONT(infinity_first_side) ? "front" : "back", - (unsigned long) segments_in[t2].global_id, + (unsigned long)segments_in[t2].global_id, SEGSIDE_IS_FRONT(cc->max_y_side_id) ? "front" : "back"); log_err(desc->scene->dev, "Segment %lu:\n (%g %g) (%g %g)\n", - (unsigned long) segments_in[t1].global_id, + (unsigned long)segments_in[t1].global_id, SPLIT2(positions[segments_in[t1].vertice_id[0]].vec), SPLIT2(positions[segments_in[t1].vertice_id[1]].vec)); log_err(desc->scene->dev, "Segment %lu:\n (%g %g) (%g %g)\n", - (unsigned long) segments_in[t2].global_id, + (unsigned long)segments_in[t2].global_id, SPLIT2(positions[segments_in[t2].vertice_id[0]].vec), SPLIT2(positions[segments_in[t2].vertice_id[1]].vec)); log_err(desc->scene->dev, "Media: %lu VS %lu\n", @@ -655,18 +663,18 @@ group_connex_components log_err(desc->scene->dev, "Medium mismatch found between segment %lu %s side and segment" " %lu %s side facing each other.\n", - (unsigned long) segments_in[t1].global_id, + (unsigned long)segments_in[t1].global_id, SEGSIDE_IS_FRONT(hit_side) ? "front" : "back", - (unsigned long) segments_in[t2].global_id, + (unsigned long)segments_in[t2].global_id, SEGSIDE_IS_FRONT(cc->max_y_side_id) ? "front" : "back"); log_err(desc->scene->dev, "Segment %lu:\n (%g %g) (%g %g)\n", - (unsigned long) segments_in[t1].global_id, + (unsigned long)segments_in[t1].global_id, SPLIT2(positions[segments_in[t1].vertice_id[0]].vec), SPLIT2(positions[segments_in[t1].vertice_id[1]].vec)); log_err(desc->scene->dev, "Segment %lu:\n (%g %g) (%g %g)\n", - (unsigned long) segments_in[t2].global_id, + (unsigned long)segments_in[t2].global_id, SPLIT2(positions[segments_in[t2].vertice_id[0]].vec), SPLIT2(positions[segments_in[t2].vertice_id[1]].vec)); log_err(desc->scene->dev, "Media: %lu VS %lu\n", @@ -683,44 +691,53 @@ group_connex_components exit_for = 1; continue; } + /* Implicit barrier here */ ASSERT(next_enclosure_id < ENCLOSURE_MAX__); - desc->enclosures_count = (enclosure_id_t)next_enclosure_id; OK(res); - /* Post-process links to group connex components */ - OK(darray_enclosure_resize(&desc->enclosures, desc->enclosures_count)); - FOR_EACH(ccc, 0, cc_count) { - component_id_t c = (component_id_t)ccc; - struct cc_descriptor* const cc = descriptors[c]; - const struct cc_descriptor* other_desc = cc; - struct enclosure_data* enclosures - = darray_enclosure_data_get(&desc->enclosures); - component_id_t fst; - - while(other_desc->enclosure_id == CC_GROUP_ID_NONE) { - other_desc = *(darray_ptr_component_descriptor_cdata_get(connex_components) - + other_desc->cc_group_root); + /* One thread post-processes links to group connex components */ + #pragma omp single + { + desc->enclosures_count = (enclosure_id_t)next_enclosure_id; + res = darray_enclosure_resize(&desc->enclosures, desc->enclosures_count); + if(res == RES_OK) { + FOR_EACH(ccc, 0, cc_count) { + component_id_t c = (component_id_t)ccc; + struct cc_descriptor* const cc = descriptors[c]; + const struct cc_descriptor* other_desc = cc; + struct enclosure_data* enclosures + = darray_enclosure_data_get(&desc->enclosures); + component_id_t fst; + + while(other_desc->enclosure_id == CC_GROUP_ID_NONE) { + other_desc = *(darray_ptr_component_descriptor_cdata_get(connex_components) + + other_desc->cc_group_root); + } + ASSERT(other_desc->cc_group_root != CC_GROUP_ROOT_NONE); + ASSERT(other_desc->enclosure_id != CC_GROUP_ID_NONE); + ASSERT(cc->medium == other_desc->medium); + cc->cc_group_root = other_desc->cc_group_root; + cc->enclosure_id = other_desc->enclosure_id; + ++enclosures[cc->enclosure_id].cc_count; + /* Linked list of componnents */ + fst = enclosures[cc->enclosure_id].first_component; + cc->enclosure_next_component = fst; + enclosures[cc->enclosure_id].first_component = cc->cc_id; + enclosures[cc->enclosure_id].side_range.first + = MMIN(enclosures[cc->enclosure_id].side_range.first, cc->side_range.first); + enclosures[cc->enclosure_id].side_range.last + = MMAX(enclosures[cc->enclosure_id].side_range.last, cc->side_range.last); + enclosures[cc->enclosure_id].side_count += cc->side_count; + } } - ASSERT(other_desc->cc_group_root != CC_GROUP_ROOT_NONE); - ASSERT(other_desc->enclosure_id != CC_GROUP_ID_NONE); - ASSERT(cc->medium == other_desc->medium); - cc->cc_group_root = other_desc->cc_group_root; - cc->enclosure_id = other_desc->enclosure_id; - ++enclosures[cc->enclosure_id].cc_count; - /* Linked list of componnents */ - fst = enclosures[cc->enclosure_id].first_component; - cc->enclosure_next_component = fst; - enclosures[cc->enclosure_id].first_component = cc->cc_id; - enclosures[cc->enclosure_id].side_range.first - = MMIN(enclosures[cc->enclosure_id].side_range.first, cc->side_range.first); - enclosures[cc->enclosure_id].side_range.last - = MMAX(enclosures[cc->enclosure_id].side_range.last, cc->side_range.last); - enclosures[cc->enclosure_id].side_count += cc->side_count; } + /* Implicit barrier here */ + OK(res); exit: return res; error: + exit_for = 1; goto exit; } @@ -728,19 +745,23 @@ static res_T collect_and_link_neighbours (struct senc2d_scene* scn, struct segside* segsides, - struct darray_segment_tmp* segments_tmp_array) + struct darray_segment_tmp* segments_tmp_array, + struct darray_neighbourhood* neighbourhood_by_vertex) { + /* This function is called from an omp parallel block and executed + * concurrently. + * Resize / Push operations on neighbourhood_by_vertex are valid + * because each neighbourhood is processes by an unique thread */ res_T res = RES_OK; const struct segment_in *segments_in; struct segment_tmp *segments_tmp; const union double2* vertices; - /* Array to keep neighbourhood of vertices. - * We create a neighbourhood for every single unique vertex, - * regardless the fact it is used by a segment - * Resize/Push operations on neighbourhood_by_vertex are valid in the - * openmp block because each neighbourhood is processes by an unique thread */ - struct darray_neighbourhood neighbourhood_by_vertex; - volatile int exit_for = 0; + const int thread_count = omp_get_num_threads(); + const int rank = omp_get_thread_num(); + vrtx_id_t v; + seg_id_t s; + /* shared between threads */ + static volatile int exit_for = 0; ASSERT(scn && segsides && segments_tmp_array); ASSERT((size_t)scn->nuverts + (size_t)scn->nusegs + 2 <= EDGE_MAX__); @@ -751,148 +772,131 @@ collect_and_link_neighbours ASSERT(scn->nusegs == darray_segment_tmp_size_get(segments_tmp_array)); - darray_neighbourhood_init(scn->dev->allocator, &neighbourhood_by_vertex); - OK2(darray_neighbourhood_resize(&neighbourhood_by_vertex, scn->nuverts), - error_); - - printf("\n"); - #pragma omp parallel - { - const int thread_count = omp_get_num_threads(); - const int rank = omp_get_thread_num(); - /* Array to keep neighbourhood of vertices - * Resize/Push operations on neighbourhood_by_vertex are valid in the - * openmp block because it is thread local data */ - vrtx_id_t v; - seg_id_t s; - /* Loop on segments to collect edges' neighbours. - * All threads considering all the vertices and processing some */ - FOR_EACH(s, 0, scn->nusegs) { - unsigned char vv; - FOR_EACH(vv, 0, 2) { - struct darray_neighbour* neighbourhood; - struct neighbour_info* info; - const vrtx_id_t vertex = segments_in[s].vertice_id[vv]; - size_t sz; - ASSERT(vertex < scn->nuverts); - if(exit_for) continue; - /* Process only "my" vertices! */ - if((int64_t)vertex % thread_count != rank) continue; - /* Find neighbourhood */ - ASSERT(vertex < darray_neighbourhood_size_get(&neighbourhood_by_vertex)); - neighbourhood = - darray_neighbourhood_data_get(&neighbourhood_by_vertex) + vertex; - sz = darray_neighbour_size_get(neighbourhood); - /* Make room for this neighbour */ - if(darray_neighbour_capacity(neighbourhood) == sz) { - /* 2 seems to be a good guess for initial capacity */ - size_t new_sz = sz ? sz + 1 : 2; - OK(darray_neighbour_reserve(neighbourhood, new_sz)); - } - OK(darray_neighbour_resize(neighbourhood, 1 + sz)); - /* Add neighbour info to vertex's neighbour list */ - info = darray_neighbour_data_get(neighbourhood) + sz; - info->seg_id = s; - info->common_vertex_rank = vv; - } - } - /* No implicit barrier here. */ - - /* For each of "my" vertices sort segments sides by rotation angle - * and connect neighbours. - * All threads considering all the vertices and processing some */ - FOR_EACH(v, 0, scn->nuverts) { - const vrtx_id_t common_vrtx = v; - vrtx_id_t other_vrtx; + /* Loop on segments to collect edges' neighbours. + * All threads considering all the vertices and processing some */ + FOR_EACH(s, 0, scn->nusegs) { + unsigned char vv; + FOR_EACH(vv, 0, 2) { struct darray_neighbour* neighbourhood; - side_id_t i, neighbour_count; + struct neighbour_info* info; + const vrtx_id_t vertex = segments_in[s].vertice_id[vv]; size_t sz; + ASSERT(vertex < scn->nuverts); + if(exit_for) continue; /* Process only "my" vertices! */ - if((int64_t)v % thread_count != rank) continue; - neighbourhood - = darray_neighbourhood_data_get(&neighbourhood_by_vertex) + v; + if((int64_t)vertex % thread_count != rank) continue; + /* Find neighbourhood */ + ASSERT(vertex < darray_neighbourhood_size_get(neighbourhood_by_vertex)); + neighbourhood = + darray_neighbourhood_data_get(neighbourhood_by_vertex) + vertex; sz = darray_neighbour_size_get(neighbourhood); - /* sz can be 0 as a vertex can be unused */ - if(!sz) continue; - ASSERT(sz <= SIDE_MAX__); - neighbour_count = (side_id_t)sz; - FOR_EACH(i, 0, neighbour_count) { - double max_y, disp[2]; - unsigned char max_y_vrank; - struct neighbour_info* neighbour_info - = darray_neighbour_data_get(neighbourhood) + i; - const struct segment_in* seg_in = segments_in + neighbour_info->seg_id; - struct segment_tmp* neighbour = segments_tmp + neighbour_info->seg_id; - other_vrtx = - seg_in->vertice_id[(neighbour_info->common_vertex_rank + 1) % 2]; - if(vertices[other_vrtx].pos.y > vertices[common_vrtx].pos.y) { - max_y = vertices[other_vrtx].pos.y; - max_y_vrank = 1 - neighbour_info->common_vertex_rank; - } else { - max_y = vertices[common_vrtx].pos.y; - max_y_vrank = neighbour_info->common_vertex_rank; - } - ASSERT(neighbour->max_y <= max_y); - neighbour->max_y = max_y; - neighbour->max_y_vrtx_rank = max_y_vrank; - /* Compute rotation angle around common vertex (in world system) */ - d2_sub(disp, vertices[other_vrtx].vec, vertices[common_vrtx].vec); - ASSERT(disp[0] || disp[1]); - neighbour_info->angle = atan2(disp[1], disp[0]); - if(neighbour_info->angle < 0) neighbour_info->angle += 2 * PI; - /* Due to catastrophic cancelation, -eps+2pi translates to 2pi */ - ASSERT(0 <= neighbour_info->angle && neighbour_info->angle <= 2 * PI); + /* Make room for this neighbour */ + if(darray_neighbour_capacity(neighbourhood) == sz) { + /* 2 seems to be a good guess for initial capacity */ + size_t new_sz = sz ? sz + 1 : 2; + OK(darray_neighbour_reserve(neighbourhood, new_sz)); } - /* Sort segments by rotation angle */ - qsort(darray_neighbour_data_get(neighbourhood), neighbour_count, - sizeof(struct neighbour_info), neighbour_cmp); - /* Link sides. - * Create cycles of sides by neighbourhood around common vertex. */ - FOR_EACH(i, 0, neighbour_count) { - /* Neighbourhood info for current pair of segments */ - const struct neighbour_info* current - = darray_neighbour_cdata_get(neighbourhood) + i; - const struct neighbour_info* ccw_neighbour - = darray_neighbour_cdata_get(neighbourhood) + (i+1) % neighbour_count; - /* Rank of the end of interest in segments */ - const unsigned char crt_end = current->common_vertex_rank; - const unsigned char ccw_end = ccw_neighbour->common_vertex_rank; - /* User id of current segments */ - const seg_id_t crt_id = current->seg_id; - const seg_id_t ccw_id = ccw_neighbour->seg_id; - /* Facing sides of segments */ - const enum side_id crt_side = crt_end ? SIDE_BACK : SIDE_FRONT; - const enum side_id ccw_side = ccw_end ? SIDE_FRONT : SIDE_BACK; - /* Index of sides in segsides */ - const side_id_t crt_side_idx = SEGIDxSIDE_2_SEGSIDE(crt_id, crt_side); - const side_id_t ccw_side_idx = SEGIDxSIDE_2_SEGSIDE(ccw_id, ccw_side); - /* Side ptrs */ - struct segside* const p_crt_side = segsides + crt_side_idx; - struct segside* const p_ccw_side = segsides + ccw_side_idx; - /* Link sides */ - p_crt_side->facing_side_id[crt_end] = ccw_side_idx; - p_ccw_side->facing_side_id[ccw_end] = crt_side_idx; - /* Record media */ - p_crt_side->medium = segments_in[crt_id].medium[crt_side]; - p_ccw_side->medium = segments_in[ccw_id].medium[ccw_side]; - ASSERT(p_crt_side->medium < scn->nmeds); - ASSERT(p_ccw_side->medium < scn->nmeds); - p_crt_side->list_id = FLAG_LIST_SIDE_LIST; - p_ccw_side->list_id = FLAG_LIST_SIDE_LIST; + OK(darray_neighbour_resize(neighbourhood, 1 + sz)); + /* Add neighbour info to vertex's neighbour list */ + info = darray_neighbour_data_get(neighbourhood) + sz; + info->seg_id = s; + info->common_vertex_rank = vv; + } + } + /* When a thread has build his share of neighbourhoods + * it can process them whithout waiting for other threads + * (no barrier needed here). */ + + /* For each of "my" vertices sort segments sides by rotation angle + * and connect neighbours. + * All threads considering all the vertices and processing some */ + FOR_EACH(v, 0, scn->nuverts) { + const vrtx_id_t common_vrtx = v; + vrtx_id_t other_vrtx; + struct darray_neighbour* neighbourhood; + side_id_t i, neighbour_count; + size_t sz; + /* Process only "my" neighbourhoods! */ + if((int64_t)v % thread_count != rank) continue; + neighbourhood + = darray_neighbourhood_data_get(neighbourhood_by_vertex) + v; + sz = darray_neighbour_size_get(neighbourhood); + /* sz can be 0 as a vertex can be unused */ + if(!sz) continue; + ASSERT(sz <= SIDE_MAX__); + neighbour_count = (side_id_t)sz; + FOR_EACH(i, 0, neighbour_count) { + double max_y, disp[2]; + unsigned char max_y_vrank; + struct neighbour_info* neighbour_info + = darray_neighbour_data_get(neighbourhood) + i; + const struct segment_in* seg_in = segments_in + neighbour_info->seg_id; + struct segment_tmp* neighbour = segments_tmp + neighbour_info->seg_id; + other_vrtx = + seg_in->vertice_id[(neighbour_info->common_vertex_rank + 1) % 2]; + if(vertices[other_vrtx].pos.y > vertices[common_vrtx].pos.y) { + max_y = vertices[other_vrtx].pos.y; + max_y_vrank = 1 - neighbour_info->common_vertex_rank; + } else { + max_y = vertices[common_vrtx].pos.y; + max_y_vrank = neighbour_info->common_vertex_rank; } + ASSERT(neighbour->max_y <= max_y); + neighbour->max_y = max_y; + neighbour->max_y_vrtx_rank = max_y_vrank; + /* Compute rotation angle around common vertex (in world system) */ + d2_sub(disp, vertices[other_vrtx].vec, vertices[common_vrtx].vec); + ASSERT(disp[0] || disp[1]); + neighbour_info->angle = atan2(disp[1], disp[0]); + if(neighbour_info->angle < 0) neighbour_info->angle += 2 * PI; + /* Due to catastrophic cancelation, -eps+2pi translates to 2pi */ + ASSERT(0 <= neighbour_info->angle && neighbour_info->angle <= 2 * PI); + } + /* Sort segments by rotation angle */ + qsort(darray_neighbour_data_get(neighbourhood), neighbour_count, + sizeof(struct neighbour_info), neighbour_cmp); + /* Link sides. + * Create cycles of sides by neighbourhood around common vertex. */ + FOR_EACH(i, 0, neighbour_count) { + /* Neighbourhood info for current pair of segments */ + const struct neighbour_info* current + = darray_neighbour_cdata_get(neighbourhood) + i; + const struct neighbour_info* ccw_neighbour + = darray_neighbour_cdata_get(neighbourhood) + (i + 1) % neighbour_count; + /* Rank of the end of interest in segments */ + const unsigned char crt_end = current->common_vertex_rank; + const unsigned char ccw_end = ccw_neighbour->common_vertex_rank; + /* User id of current segments */ + const seg_id_t crt_id = current->seg_id; + const seg_id_t ccw_id = ccw_neighbour->seg_id; + /* Facing sides of segments */ + const enum side_id crt_side = crt_end ? SIDE_BACK : SIDE_FRONT; + const enum side_id ccw_side = ccw_end ? SIDE_FRONT : SIDE_BACK; + /* Index of sides in segsides */ + const side_id_t crt_side_idx = SEGIDxSIDE_2_SEGSIDE(crt_id, crt_side); + const side_id_t ccw_side_idx = SEGIDxSIDE_2_SEGSIDE(ccw_id, ccw_side); + /* Side ptrs */ + struct segside* const p_crt_side = segsides + crt_side_idx; + struct segside* const p_ccw_side = segsides + ccw_side_idx; + /* Link sides */ + p_crt_side->facing_side_id[crt_end] = ccw_side_idx; + p_ccw_side->facing_side_id[ccw_end] = crt_side_idx; + /* Record media */ + p_crt_side->medium = segments_in[crt_id].medium[crt_side]; + p_ccw_side->medium = segments_in[ccw_id].medium[ccw_side]; + ASSERT(p_crt_side->medium < scn->nmeds); + ASSERT(p_ccw_side->medium < scn->nmeds); + p_crt_side->list_id = FLAG_LIST_SIDE_LIST; + p_ccw_side->list_id = FLAG_LIST_SIDE_LIST; } - /* jump error block */ - goto after_error; -error: - /* Cannot goto out of openmp block */ - exit_for = 1; -after_error: - ; } -error_: - darray_neighbourhood_release(&neighbourhood_by_vertex); - +end: + /* Threads are allowed to return whitout sync. */ return res; +error: + /* Cannot goto out of openmp block */ + exit_for = 1; + goto end; } static res_T @@ -901,6 +905,8 @@ build_result const struct darray_ptr_component_descriptor* connex_components, const struct darray_segment_comp* segments_comp_array) { + /* This function is called from an omp parallel block and executed + * concurrently. */ res_T res = RES_OK; struct mem_allocator* alloc; struct cc_descriptor* const* cc_descriptors; @@ -908,7 +914,11 @@ build_result const struct segment_in* segments_in; struct segment_enc* segments_enc; const struct segment_comp* segments_comp; - volatile int exit_for = 0; + struct htable_vrtx_id vtable; + int64_t sg; + int64_t ee; + /* shared between threads */ + static volatile int exit_for = 0; ASSERT(desc && connex_components && segments_comp_array); @@ -918,126 +928,123 @@ build_result enclosures = darray_enclosure_data_get(&desc->enclosures); segments_in = darray_segment_in_cdata_get(&desc->scene->segments_in); segments_comp = darray_segment_comp_cdata_get(segments_comp_array); - OK2(darray_segment_enc_resize(&desc->segments_enc, desc->scene->nusegs), - error_); - segments_enc = darray_segment_enc_data_get(&desc->segments_enc); - - #pragma omp parallel + #pragma omp single { - struct htable_vrtx_id vtable; - int64_t sg; - int64_t ee; - - /* Build global enclosure information */ - #pragma omp for - for(sg = 0; sg < (int64_t) desc->scene->nusegs; sg++) { - seg_id_t s = (seg_id_t) sg; - const component_id_t cf_id = segments_comp[s].component[SIDE_FRONT]; - const component_id_t cb_id = segments_comp[s].component[SIDE_BACK]; - const struct cc_descriptor* cf = cc_descriptors[cf_id]; - const struct cc_descriptor* cb = cc_descriptors[cb_id]; - const enclosure_id_t ef_id = cf->enclosure_id; - const enclosure_id_t eb_id = cb->enclosure_id; - ASSERT(segments_enc[s].enclosure[SIDE_FRONT] == ENCLOSURE_NULL__); - segments_enc[s].enclosure[SIDE_FRONT] = ef_id; - ASSERT(segments_enc[s].enclosure[SIDE_BACK] == ENCLOSURE_NULL__); - segments_enc[s].enclosure[SIDE_BACK] = eb_id; - } - /* Implicit barrier here */ - - /* Resize/push operations on enclosure's fields are valid in the - * openmp block as a given enclosure is processed by a single thread */ - htable_vrtx_id_init(alloc, &vtable); - - ASSERT(desc->enclosures_count <= ENCLOSURE_MAX__); - #pragma omp for schedule(dynamic) nowait - for(ee = 0; ee < (int64_t)desc->enclosures_count; ee++) { - const enclosure_id_t e = (enclosure_id_t)ee; - struct enclosure_data* enc = enclosures + e; - const struct cc_descriptor* current = cc_descriptors[enc->first_component]; - seg_id_t fst_idx = 0; - seg_id_t sgd_idx = enc->side_count; - seg_id_t s; - ASSERT(enc->first_component - < darray_ptr_component_descriptor_size_get(connex_components)); - ASSERT(current->cc_id == enc->first_component); + res = darray_segment_enc_resize(&desc->segments_enc, desc->scene->nusegs); + }/* Implicit barrier here. */ + OK2(res, error_); + segments_enc = darray_segment_enc_data_get(&desc->segments_enc); - if(exit_for) continue; - ASSERT(e <= UINT_MAX); - enc->header.enclosure_id = (unsigned)e; /* Back to API type */ - ASSERT(current->enclosure_id == enc->header.enclosure_id); - enc->header.is_infinite = (e == 0); - enc->header.enclosed_medium - = (unsigned)current->medium; /* Back to API type */ - ASSERT(enc->header.enclosed_medium < desc->scene->nmeds); - - /* Build side and vertex lists. */ - OK(darray_segment_in_resize(&enc->sides, enc->side_count)); - /* Size is just a int */ - OK(darray_vrtx_id_reserve(&enc->vertices, enc->side_count + 1)); - /* New vertex numbering scheme local to the enclosure */ - htable_vrtx_id_clear(&vtable); - ASSERT(desc->scene->nusegs - == darray_segment_in_size_get(&desc->scene->segments_in)); - /* Put at the end the back-faces of segments that also have their - * front-face in the list. */ - for(s = SEGSIDE_2_SEG(enc->side_range.first); - s <= SEGSIDE_2_SEG(enc->side_range.last); - s++) - { - const struct segment_in* seg_in = segments_in + s; - struct segment_in* seg; - unsigned vertice_id[2]; - int i; - if(segments_enc[s].enclosure[SIDE_FRONT] != e - && segments_enc[s].enclosure[SIDE_BACK] != e) - continue; - ++enc->header.unique_segment_count; + /* Build global enclosure information */ + #pragma omp for + for(sg = 0; sg < (int64_t)desc->scene->nusegs; sg++) { + seg_id_t s = (seg_id_t)sg; + const component_id_t cf_id = segments_comp[s].component[SIDE_FRONT]; + const component_id_t cb_id = segments_comp[s].component[SIDE_BACK]; + const struct cc_descriptor* cf = cc_descriptors[cf_id]; + const struct cc_descriptor* cb = cc_descriptors[cb_id]; + const enclosure_id_t ef_id = cf->enclosure_id; + const enclosure_id_t eb_id = cb->enclosure_id; + ASSERT(segments_enc[s].enclosure[SIDE_FRONT] == ENCLOSURE_NULL__); + segments_enc[s].enclosure[SIDE_FRONT] = ef_id; + ASSERT(segments_enc[s].enclosure[SIDE_BACK] == ENCLOSURE_NULL__); + segments_enc[s].enclosure[SIDE_BACK] = eb_id; + } + /* Implicit barrier here */ + + /* Resize/push operations on enclosure's fields are valid in the + * openmp block as a given enclosure is processed by a single thread */ + htable_vrtx_id_init(alloc, &vtable); + + ASSERT(desc->enclosures_count <= ENCLOSURE_MAX__); + #pragma omp for schedule(dynamic) nowait + for(ee = 0; ee < (int64_t)desc->enclosures_count; ee++) { + const enclosure_id_t e = (enclosure_id_t)ee; + struct enclosure_data* enc = enclosures + e; + const struct cc_descriptor* current = cc_descriptors[enc->first_component]; + seg_id_t fst_idx = 0; + seg_id_t sgd_idx = enc->side_count; + seg_id_t s; + ASSERT(enc->first_component + < darray_ptr_component_descriptor_size_get(connex_components)); + ASSERT(current->cc_id == enc->first_component); - FOR_EACH(i, 0, 2) { - vrtx_id_t* id = htable_vrtx_id_find(&vtable, seg_in->vertice_id + i); - if(id) { - vertice_id[i] = *id; /* Known vertex */ - } else { - /* Create new association */ - size_t tmp = htable_vrtx_id_size_get(&vtable); - ASSERT(tmp == darray_vrtx_id_size_get(&enc->vertices)); - ASSERT(tmp < VRTX_MAX__); - vertice_id[i] = (vrtx_id_t)tmp; - OK(htable_vrtx_id_set(&vtable, seg_in->vertice_id + i, - vertice_id + i)); - OK(darray_vrtx_id_push_back(&enc->vertices, seg_in->vertice_id + i)); - ++enc->header.vertices_count; - } - } - ASSERT(segments_enc[s].enclosure[SIDE_FRONT] == e - || segments_enc[s].enclosure[SIDE_BACK] == e); - if(segments_enc[s].enclosure[SIDE_FRONT] == e) { - ++enc->header.segment_count; - seg = darray_segment_in_data_get(&enc->sides) + fst_idx++; - - FOR_EACH(i, 0, 2) seg->medium[i] = seg_in->medium[i]; - seg->global_id = seg_in->global_id; - FOR_EACH(i, 0, 2) seg->vertice_id[i] = vertice_id[i]; - } - if(segments_enc[s].enclosure[SIDE_BACK] == e) { - ++enc->header.segment_count; - seg = darray_segment_in_data_get(&enc->sides) + - ((segments_enc[s].enclosure[SIDE_FRONT] == e) ? --sgd_idx : fst_idx++); - - FOR_EACH(i, 0, 2) seg->medium[i] = seg_in->medium[1 - i]; - seg->global_id = seg_in->global_id; - FOR_EACH(i, 0, 2) seg->vertice_id[i] = vertice_id[1 - i]; + if(exit_for) continue; + ASSERT(e <= UINT_MAX); + enc->header.enclosure_id = (unsigned)e; /* Back to API type */ + ASSERT(current->enclosure_id == enc->header.enclosure_id); + enc->header.is_infinite = (e == 0); + enc->header.enclosed_medium + = (unsigned)current->medium; /* Back to API type */ + ASSERT(enc->header.enclosed_medium < desc->scene->nmeds); + + /* Build side and vertex lists. */ + OK(darray_segment_in_resize(&enc->sides, enc->side_count)); + /* Size is just a int */ + OK(darray_vrtx_id_reserve(&enc->vertices, enc->side_count + 1)); + /* New vertex numbering scheme local to the enclosure */ + htable_vrtx_id_clear(&vtable); + ASSERT(desc->scene->nusegs + == darray_segment_in_size_get(&desc->scene->segments_in)); + /* Put at the end the back-faces of segments that also have their + * front-face in the list. */ + for(s = SEGSIDE_2_SEG(enc->side_range.first); + s <= SEGSIDE_2_SEG(enc->side_range.last); + s++) + { + const struct segment_in* seg_in = segments_in + s; + struct segment_in* seg; + unsigned vertice_id[2]; + int i; + if(segments_enc[s].enclosure[SIDE_FRONT] != e + && segments_enc[s].enclosure[SIDE_BACK] != e) + continue; + ++enc->header.unique_segment_count; + + FOR_EACH(i, 0, 2) { + vrtx_id_t* id = htable_vrtx_id_find(&vtable, seg_in->vertice_id + i); + if(id) { + vertice_id[i] = *id; /* Known vertex */ + } else { + /* Create new association */ + size_t tmp = htable_vrtx_id_size_get(&vtable); + ASSERT(tmp == darray_vrtx_id_size_get(&enc->vertices)); + ASSERT(tmp < VRTX_MAX__); + vertice_id[i] = (vrtx_id_t)tmp; + OK(htable_vrtx_id_set(&vtable, seg_in->vertice_id + i, + vertice_id + i)); + OK(darray_vrtx_id_push_back(&enc->vertices, seg_in->vertice_id + i)); + ++enc->header.vertices_count; } - if(fst_idx == sgd_idx) break; } - continue; - error: - /* Cannot goto out of openmp block */ - exit_for = 1; + ASSERT(segments_enc[s].enclosure[SIDE_FRONT] == e + || segments_enc[s].enclosure[SIDE_BACK] == e); + if(segments_enc[s].enclosure[SIDE_FRONT] == e) { + ++enc->header.segment_count; + seg = darray_segment_in_data_get(&enc->sides) + fst_idx++; + + FOR_EACH(i, 0, 2) seg->medium[i] = seg_in->medium[i]; + seg->global_id = seg_in->global_id; + FOR_EACH(i, 0, 2) seg->vertice_id[i] = vertice_id[i]; + } + if(segments_enc[s].enclosure[SIDE_BACK] == e) { + ++enc->header.segment_count; + seg = darray_segment_in_data_get(&enc->sides) + + ((segments_enc[s].enclosure[SIDE_FRONT] == e) ? --sgd_idx : fst_idx++); + + FOR_EACH(i, 0, 2) seg->medium[i] = seg_in->medium[1 - i]; + seg->global_id = seg_in->global_id; + FOR_EACH(i, 0, 2) seg->vertice_id[i] = vertice_id[1 - i]; + } + if(fst_idx == sgd_idx) break; } - htable_vrtx_id_release(&vtable); - } + continue; + error: + /* Cannot goto out of openmp block */ + exit_for = 1; + } /* No barrier here */ + htable_vrtx_id_release(&vtable); + OK2(res, error_); exit: @@ -1061,18 +1068,19 @@ senc2d_scene_analyze(struct senc2d_scene* scn, struct senc2d_descriptor** out_de * They are refered to by arrays of ids. */ struct darray_ptr_component_descriptor connex_components; char connex_components_initialized = 0; - /* Segment neighbourhood by edge. */ - struct darray_neighbourhood neighbourhood_by_edge; - char neighbourhood_by_edge_initialized = 0; /* Store by-segment components */ struct darray_segment_comp segments_comp; char segments_comp_initialized = 0; /* Array of vertices (segment ends). */ struct segside* segsides = NULL; struct s2d_scene_view* s2d_view = NULL; + /* Array to keep neighbourhood of vertices */ + struct darray_neighbourhood neighbourhood_by_vertex; + char neighbourhood_by_vertex_initialized = 0; if(!scn || !out_desc) return RES_BAD_ARG; + /* The first part of the analyze is single threaded */ desc = descriptor_create(scn); if(!desc) { res = RES_MEM_ERR; @@ -1092,66 +1100,157 @@ senc2d_scene_analyze(struct senc2d_scene* scn, struct senc2d_descriptor** out_de goto error; } - /* Step 1: build neighbourhoods */ - res = collect_and_link_neighbours(scn, segsides, &segments_tmp); + /* We create a neighbourhood for every single unique vertex, + * regardless the fact it is used by a segment. + * This allows threads to use these neighbourhoods whithout syn. */ + darray_neighbourhood_init(scn->dev->allocator, &neighbourhood_by_vertex); + neighbourhood_by_vertex_initialized = 1; + OK(darray_neighbourhood_resize(&neighbourhood_by_vertex, scn->nuverts)); - if(res != RES_OK) { - log_err(scn->dev, - "%s: could not build neighbourhoods from scene.\n", FUNC_NAME); - goto error; - } + /* The end of the analyze is multithreaded */ + #pragma omp parallel + { + res_T thread_res = RES_OK; + + /* Step 1: build neighbourhoods */ + thread_res = collect_and_link_neighbours(scn, segsides, &segments_tmp, + &neighbourhood_by_vertex); + /* No barrier at the end of step 1: data used in step 1 cannot be + * released / data produced by step 1 cannot be used + * until next sync point */ + + if(thread_res != RES_OK) { + res = thread_res; + #pragma omp single nowait + { + log_err(scn->dev, + "%s: could not build neighbourhoods from scene.\n", FUNC_NAME); + } /* No barrier here */ + } + OK2(res, error_); + + /* The first thread here allocates some data */ + #pragma omp single + { + darray_ptr_component_descriptor_init(scn->dev->allocator, + &connex_components); + connex_components_initialized = 1; + /* Just a hint; to limit contention */ + thread_res = darray_ptr_component_descriptor_reserve(&connex_components, + 2 * scn->nmeds); + if(thread_res != RES_OK) res = thread_res; + darray_segment_comp_init(scn->dev->allocator, &segments_comp); + segments_comp_initialized = 1; + thread_res = darray_segment_comp_resize(&segments_comp, scn->nusegs); + if(thread_res != RES_OK) res = thread_res; + } + /* Implicit barrier here: constraints on step 1 data are now met */ + OK2(res, error_); - darray_ptr_component_descriptor_init(scn->dev->allocator, &connex_components); - connex_components_initialized = 1; - darray_segment_comp_init(scn->dev->allocator, &segments_comp); - segments_comp_initialized = 1; - OK(darray_segment_comp_resize(&segments_comp, scn->nusegs)); - - /* Step 2: extract segment connex components */ - res = extract_connex_components(desc, segsides, &connex_components, - &segments_tmp, &segments_comp, &s2d_view); - if(res != RES_OK) { - log_err(scn->dev, - "%s: could not extract connex components from scene.\n", FUNC_NAME); - goto error; - } + /* One thread releases some data before going to step 2, + * the others go to step 2 without sync */ + #pragma omp single nowait + { + darray_neighbourhood_release(&neighbourhood_by_vertex); + neighbourhood_by_vertex_initialized = 0; + } /* Barrier here */ + + /* Step 2: extract segment connex components */ + thread_res = extract_connex_components(desc, segsides, &connex_components, + &segments_tmp, &segments_comp, &s2d_view); + /* No barrier at the end of step 2: data used in step 2 cannot be + * released / data produced by step 2 cannot be used + * until next sync point */ + + if(thread_res != RES_OK) { + res = thread_res; + #pragma omp single nowait + { + log_err(scn->dev, + "%s: could not extract connex components from scene.\n", FUNC_NAME); + } /* No barrier here */ + } + OK2(res, error_); - darray_segment_tmp_release(&segments_tmp); - segments_tmp_initialized = 0; + #pragma omp barrier + /* Constraints on step 2 data are now met */ - /* Step 3: group components */ - res = group_connex_components(desc, segsides, &segments_comp, - &connex_components, s2d_view); - if (s2d_view) S2D(scene_view_ref_put(s2d_view)); - if(res != RES_OK) { - log_err(scn->dev, - "%s: could not group connex components from scene.\n", FUNC_NAME); - goto error; - } + /* One thread releases some data before going to step 3, + * the others go to step 3 without sync */ + #pragma omp single nowait + { + darray_segment_tmp_release(&segments_tmp); + segments_tmp_initialized = 0; + } /* No barrier here */ + + /* Step 3: group components */ + thread_res = group_connex_components(desc, segsides, &segments_comp, + &connex_components, s2d_view); + /* Barrier at the end of step 3: data used in step 3 can be released / + * data produced by step 2 can be used */ + + /* One thread releases some data before going to step 4, + * the others go to step 4 without sync */ + #pragma omp single nowait + { + if(s2d_view) S2D(scene_view_ref_put(s2d_view)); + s2d_view = NULL; + } /* No barrier here */ + + if(res != RES_OK) { + res = thread_res; + #pragma omp single nowait + { + log_err(scn->dev, + "%s: could not group connex components from scene.\n", FUNC_NAME); + } /* No barrier here */ + } + OK2(res, error_); - /* Build result. */ - res = build_result(desc, &connex_components, &segments_comp); - if(res != RES_OK) { - log_err(scn->dev, "%s: could not build result.\n", FUNC_NAME); - goto error; - } + /* Step 4: Build result */ + thread_res = build_result(desc, &connex_components, &segments_comp); + /* No barrier at the end of step 4: data used in step 4 cannot be + * released / data produced by step 4 cannot be used + * until next sync point */ - darray_segment_comp_release(&segments_comp); - segments_comp_initialized = 0; + if(thread_res != RES_OK) { + res = thread_res; + #pragma omp single nowait + { + log_err(scn->dev, "%s: could not build result.\n", FUNC_NAME); + } /* No barrier here */ + } + OK2(res, error_); + + #pragma omp barrier + /* Constraints on step 4 data are now met */ + /* Some threads release data */ + #pragma omp sections nowait + { + #pragma omp section + { + ASSERT(connex_components_initialized); + custom_darray_ptr_component_descriptor_release(&connex_components); + connex_components_initialized = 0; + } + #pragma omp section + { + darray_segment_comp_release(&segments_comp); + segments_comp_initialized = 0; + } + } /* No barrier here */ +error_: + ; + } /* Implicit barrier here */ + +if(res != RES_OK) goto error; exit: - if(connex_components_initialized) { - size_t c, cc_count = - darray_ptr_component_descriptor_size_get(&connex_components); - struct cc_descriptor** components = - darray_ptr_component_descriptor_data_get(&connex_components); - FOR_EACH(c, 0, cc_count) { - ptr_component_descriptor_release(scn->dev->allocator, components + c); - } - darray_ptr_component_descriptor_release(&connex_components); - } - if(neighbourhood_by_edge_initialized) - darray_neighbourhood_release(&neighbourhood_by_edge); + if(connex_components_initialized) + custom_darray_ptr_component_descriptor_release(&connex_components); + if(s2d_view) S2D(scene_view_ref_put(s2d_view)); + if(neighbourhood_by_vertex_initialized) + darray_neighbourhood_release(&neighbourhood_by_vertex); if(segments_tmp_initialized) darray_segment_tmp_release(&segments_tmp); if(segments_comp_initialized) darray_segment_comp_release(&segments_comp); if(segsides) MEM_RM(scn->dev->allocator, segsides); diff --git a/src/senc2d_scene_analyze.c.save b/src/senc2d_scene_analyze.c.save @@ -0,0 +1,1254 @@ +/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program 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 General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. */ + +#include "senc2d.h" +#include "senc2d_descriptor_c.h" +#include "senc2d_device_c.h" +#include "senc2d_scene_c.h" +#include "senc2d_scene_analyze_c.h" +#include "senc2d_internal_types.h" + +#include <rsys/rsys.h> +#include <rsys/float2.h> +#include <rsys/mem_allocator.h> +#include <rsys/hash_table.h> +#include <rsys/dynamic_array.h> + +#if defined(COMPILER_GCC) +#define ATOMIC_CAS_PTR(Atom, NewVal, Comparand) /* Return the initial value */ \ + ATOMIC_CAS((Atom), (NewVal), (Comparand)) +#elif defined(COMPILER_CL) +#define ATOMIC_CAS_PTR(Atom, NewVal, Comparand) /* Return the initial value */ \ + (InterlockedCompareExchangePointer(Atom, NewVal, Comparand)) +#else +#error "Undefined atomic operations" +#endif + +#include <star/s2d.h> + +#include <omp.h> +#include <limits.h> +#include <stdlib.h> + +#define CC_DESCRIPTOR_NULL__ {\ + {0,-DBL_MAX}, -1, SIDE_NULL__, VRTX_NULL__, 0, MEDIUM_NULL__,\ + CC_ID_NONE, CC_GROUP_ROOT_NONE, CC_GROUP_ID_NONE, CC_ID_NONE,\ + { SEG_NULL__, 0}\ +} +const struct cc_descriptor CC_DESCRIPTOR_NULL = CC_DESCRIPTOR_NULL__; + +#define DARRAY_NAME component_id +#define DARRAY_DATA component_id_t +#include <rsys/dynamic_array.h> + +#define DARRAY_NAME side_id +#define DARRAY_DATA side_id_t +#include <rsys/dynamic_array.h> + +/******************************************************************************* + * Helper function + ******************************************************************************/ +static INLINE int +find_side_in_list + (const struct segside* segsides, + const struct darray_side_id* side_ids, + const side_id_t side_id, + const enum list_id list_id) +{ + side_id_t i; + size_t tmp; + (void)list_id; + (void)segsides; + ASSERT(segsides && side_ids); + tmp = darray_side_id_size_get(side_ids); + ASSERT(tmp <= SIDE_MAX__); + FOR_EACH(i, 0, (side_id_t)tmp) { + const side_id_t id = darray_side_id_cdata_get(side_ids)[i]; + ASSERT(segsides[id].list_id & list_id); + if(id == side_id) return 1; + } + return 0; +} + +static INLINE int +neighbour_cmp(const void* w1, const void* w2) +{ + const double a1 = ((struct neighbour_info*)w1)->angle; + const double a2 = ((struct neighbour_info*)w2)->angle; + return (a1 > a2) - (a1 < a2); +} + +static FINLINE void +add_side_to_stack + (const struct senc2d_scene* scn, + struct darray_side_id* stack, + struct segside* segsides, + const side_id_t side_id) +{ + (void)scn; + ASSERT(scn && segsides && stack + && side_id < SIDE_MAX__ && side_id < 2 * scn->nusegs); + ASSERT((darray_side_id_size_get(stack) > 128) + || !find_side_in_list(segsides, stack, side_id, FLAG_WAITING_STACK)); + darray_side_id_push_back(stack, &side_id); + segsides[side_id].list_id = FLAG_WAITING_STACK; +} + +static side_id_t +get_side_not_in_connex_component + (const struct senc2d_scene* scn, + const struct segside* segsides, + side_id_t* first_side_not_in_component, + const medium_id_t medium) +{ + side_id_t i; + const seg_id_t last_side + = darray_side_range_cdata_get(&scn->media_use)[medium].last; + ASSERT(scn && segsides && first_side_not_in_component); + i = *first_side_not_in_component; + while(i <= last_side + && (segsides[i].medium != medium + || segsides[i].list_id != FLAG_LIST_SIDE_LIST)) { + ++i; + } + + *first_side_not_in_component = i+1; + if(i > last_side) return SIDE_NULL__; + return i; +} + +static FINLINE side_id_t +get_side_from_stack + (const struct segside* segsides, + struct darray_side_id* stack) +{ + side_id_t id; + size_t sz; + (void)segsides; + ASSERT(segsides && stack); + sz = darray_side_id_size_get(stack); + ASSERT(sz); + id = darray_side_id_cdata_get(stack)[sz - 1]; + ASSERT(segsides[id].list_id == FLAG_WAITING_STACK); + darray_side_id_pop_back(stack); + return id; +} + +static void +get_scn_indices(const unsigned iseg, unsigned ids[2], void* ctx) { + int i; + const struct senc2d_scene* scene = ctx; + const struct segment_in* seg = + darray_segment_in_cdata_get(&scene->segments_in) + iseg; + FOR_EACH(i, 0, 2) { + ASSERT(seg->vertice_id[i] < scene->nverts); + ids[i] = (unsigned)seg->vertice_id[i]; /* Back to API type */ + } +} + +static void +get_scn_position(const unsigned ivert, float pos[2], void* ctx) { + const struct senc2d_scene* scene = ctx; + const union double2* pt = + darray_position_cdata_get(&scene->vertices) + ivert; + f2_set_d2(pos, pt->vec); +} + +static int +self_hit_filter + (const struct s2d_hit* hit, + const float ray_org[2], + const float ray_dir[2], + void* ray_data, + void* filter_data) +{ + const struct darray_segment_comp* segments_comp = filter_data; + const component_id_t* origin_component = ray_data; + const struct segment_comp* hit_seg_comp; + enum side_id hit_side; + component_id_t hit_component; + + (void) ray_org; (void) ray_dir; + ASSERT(hit && segments_comp && origin_component); + ASSERT(hit->prim.prim_id < darray_segment_comp_size_get(segments_comp)); + hit_seg_comp = darray_segment_comp_cdata_get(segments_comp) + + hit->prim.prim_id; + hit_side = (hit->normal[1] > 0) ? SIDE_FRONT : SIDE_BACK; + hit_component = hit_seg_comp->component[hit_side]; + + /* Not self hit or distance should be small */ + ASSERT(hit_component != *origin_component || hit->distance < 1e-6); + return (hit_component == *origin_component); +} + +static res_T +extract_connex_components + (struct senc2d_descriptor* desc, + struct segside* segsides, + struct darray_ptr_component_descriptor* connex_components, + const struct darray_segment_tmp* segments_tmp_array, + struct darray_segment_comp* segments_comp, + struct s2d_scene_view** s2d_view) +{ + res_T res = RES_OK; + const struct senc2d_scene* scn; + struct mem_allocator* alloc; + static ATOMIC component_count = 0; + static volatile int exit_for = 0; + int64_t mm; +#ifndef NDEBUG + seg_id_t s_; + component_id_t c; +#endif + + ASSERT(segsides && desc && connex_components && segments_tmp_array); + ASSERT(darray_ptr_component_descriptor_size_get(connex_components) == 0); + alloc = descriptor_get_allocator(desc); + scn = desc->scene; + +#ifndef NDEBUG + #pragma omp single + { + FOR_EACH(s_, 0, scn->nusegs) { + const struct segment_in* seg_in = + darray_segment_in_cdata_get(&scn->segments_in) + s_; + const struct side_range* media_use + = darray_side_range_cdata_get(&scn->media_use); + FOR_EACH(mm, 0, 2) { + const side_id_t side = SEGIDxSIDE_2_SEGSIDE(s_, mm); + const medium_id_t medium = seg_in->medium[mm]; + ASSERT(media_use[medium].first <= side && side <= media_use[medium].last); + } + } + } /* Implicit barrier here */ +#endif + + /* #pragma omp parallel */ + { + struct darray_side_id stack; + darray_side_id_init(alloc, &stack); + + #pragma omp for schedule(dynamic) nowait + for(mm = 0; mm < (int64_t)scn->nmeds; mm++) { /* Process all media */ + const medium_id_t m = (medium_id_t)mm; + struct cc_descriptor* cc; + /* Any not-already-used side is used as a starting point */ + side_id_t first_side_not_in_component; + + if(exit_for) continue; + first_side_not_in_component + = darray_side_range_cdata_get(&scn->media_use)[m].first; + if(first_side_not_in_component == SIDE_NULL__) + continue; /* Unused medium */ + ASSERT(first_side_not_in_component < 2 * scn->nusegs); + ASSERT(darray_side_id_size_get(&stack) == 0); + for(;;) { /* Process all components for this medium */ + const side_id_t start_side_id = get_side_not_in_connex_component(scn, + segsides, &first_side_not_in_component, m); + side_id_t crt_side_id = start_side_id; + side_id_t last_side_id = start_side_id; + ASSERT(start_side_id == SIDE_NULL__ || start_side_id < 2 * scn->nusegs); + if(start_side_id == SIDE_NULL__) + break; /* start_side_id=SIDE_NULL__ => done! */ + ASSERT(segsides[start_side_id].list_id == FLAG_LIST_SIDE_LIST); + +#ifndef NDEBUG + { + seg_id_t tid = SEGSIDE_2_SEG(start_side_id); + enum side_flag s = SEGSIDE_2_SIDE(start_side_id); + medium_id_t side_med + = darray_segment_in_data_get(&desc->scene->segments_in)[tid].medium[s]; + ASSERT(side_med == m); + } +#endif + + /* Create and init a new component */ + cc = MEM_ALLOC(alloc, sizeof(struct cc_descriptor)); + if(!cc) { + res = RES_MEM_ERR; + goto error1; + } + cc_descriptor_init(alloc, cc); + ASSERT(m == segsides[start_side_id].medium); + cc->cc_id = (component_id_t)(ATOMIC_INCR(&component_count) - 1); + cc->medium = m; + cc->side_range.first = start_side_id; + + for(;;) { /* Process all sides of this component */ + int i; + enum side_flag crt_side_flag = SEGSIDE_2_SIDE(crt_side_id); + struct segside* crt_side = segsides + crt_side_id; + const seg_id_t crt_seg_id = SEGSIDE_2_SEG(crt_side_id); + const struct segment_in* seg_in = + darray_segment_in_cdata_get(&scn->segments_in) + crt_seg_id; + struct segment_comp* seg_comp = + darray_segment_comp_data_get(segments_comp) + crt_seg_id; + const struct segment_tmp* const seg_tmp = + darray_segment_tmp_cdata_get(segments_tmp_array) + crt_seg_id; + const union double2* vertices = + darray_position_cdata_get(&scn->vertices); + ASSERT(crt_seg_id < scn->nusegs); + ASSERT(segsides[crt_side_id].medium == m); + + /* Record Ymax information + * Keep track of the appropriate vertex/side of the connex component + * in order to cast a ray at the component grouping step of the + * algorithm. + * The most appropriate vertex is (the) one with the greater Y + * coordinate. + * If more than one vertex/side has the same Y, we want the side that + * most faces Y (that is the one with the greater ny). + * This is mandatory to select the correct side when both sides of a + * segment are candidate. */ + if(cc->max_vrtx[1] <= seg_tmp->max_y) { + /* Can either improve y or ny */ + + if(cc->max_y_side_id == SEGSIDE_OPPOSITE(crt_side_id)) { + /* Both sides are in cc and the opposite side is currently selected + * Just keep the side with ny>0 */ + if(cc->max_y_ny < 0) { + /* Change side! */ + cc->max_y_ny = -cc->max_y_ny; + cc->max_y_side_id = crt_side_id; + } + } else { + double edge[2], normal[2], norm, side_ny; + int process = 0; + + d2_sub(edge, vertices[seg_in->vertice_id[1]].vec, + vertices[seg_in->vertice_id[0]].vec); + d2(normal, edge[1], -edge[0]); + norm = d2_normalize(normal, normal); + ASSERT(norm); (void) norm; + + /* Geometrical normal points toward the right */ + if(SEGSIDE_IS_FRONT(crt_side_id)) { + side_ny = -normal[1]; + process = 1; + } else { + side_ny = normal[1]; + process = 1; + } + + if(process) { + int change = 0; + if(cc->max_vrtx[1] < seg_tmp->max_y) { + change = 1; /* Try first to improve y */ + } else if(cc->max_y_ny <= 0 && fabs(cc->max_y_ny) < fabs(side_ny)) { + change = 1; /* If ny <= 0, the more negative the better */ + } else if(cc->max_y_ny > 0 && cc->max_y_ny < side_ny) { + change = 1; /* If ny > 0, the more positive the better */ + } + + if(change) { + cc->max_y_ny = side_ny; + cc->max_y_side_id = crt_side_id; + ASSERT(seg_tmp->max_y_vrtx_rank < 2); + ASSERT(seg_in->vertice_id[seg_tmp->max_y_vrtx_rank] < scn->nverts); + cc->max_y_vrtx_id = seg_in->vertice_id[seg_tmp->max_y_vrtx_rank]; + ASSERT(seg_tmp->max_y == vertices[cc->max_y_vrtx_id].pos.y); + d2_set(cc->max_vrtx, vertices[cc->max_y_vrtx_id].vec); + } + } + } + } + + /* Record crt_side both as component and segment level */ + cc->side_count++; + segsides[crt_side_id].list_id = FLAG_LIST_COMPONENT; + ASSERT(seg_comp->component[crt_side_flag] == COMPONENT_NULL__); + seg_comp->component[crt_side_flag] = cc->cc_id; +#ifndef NDEBUG + crt_side->member_of_cc = cc->cc_id; +#endif + + /* Store neighbour sides in a waiting stack */ + FOR_EACH(i, 0, 2) { + side_id_t neighbour_id = crt_side->facing_side_id[i]; + seg_id_t nbour_seg_id = SEGSIDE_2_SEG(neighbour_id); + const struct segside* neighbour = segsides + neighbour_id; + ASSERT(m == crt_side->medium); + if(neighbour->medium != crt_side->medium) { + /* Found medium discontinuity! Model topology is broken. */ + const struct segment_in* segments_in + = darray_segment_in_cdata_get(&scn->segments_in); + const union double2* positions + = darray_position_cdata_get(&scn->vertices); + log_err(scn->dev, + "Medium mismatch found between neighbour segments %lu %s" + " side and %u %s side.\n", + (unsigned long) segments_in[crt_seg_id].global_id, + SEGSIDE_IS_FRONT(crt_side_id) ? "front" : "back", + segments_in[nbour_seg_id].global_id, + SEGSIDE_IS_FRONT(neighbour_id) ? "front" : "back"); + log_err(scn->dev, + "Segment %lu:\n (%g %g) (%g %g))\n", + (unsigned long) segments_in[crt_seg_id].global_id, + SPLIT2(positions[segments_in[crt_seg_id].vertice_id[0]].vec), + SPLIT2(positions[segments_in[crt_seg_id].vertice_id[1]].vec)); + log_err(scn->dev, + "Segment %lu:\n (%g %g) (%g %g)\n", + (unsigned long) segments_in[nbour_seg_id].global_id, + SPLIT2(positions[segments_in[nbour_seg_id].vertice_id[0]].vec), + SPLIT2(positions[segments_in[nbour_seg_id].vertice_id[1]].vec)); + log_err(desc->scene->dev, "Media: %lu VS %lu\n", + (unsigned long)neighbour->medium, (unsigned long)crt_side->medium); + res = RES_BAD_ARG; + goto error1; + } + if(neighbour->list_id == FLAG_LIST_COMPONENT) { + /* Already processed */ +#ifndef NDEBUG + ASSERT(neighbour->member_of_cc == cc->cc_id); +#endif + continue; + } + if(neighbour->list_id == FLAG_WAITING_STACK) { + continue; /* Already processed */ + } + add_side_to_stack(scn, &stack, segsides, neighbour_id); + } + if(darray_side_id_size_get(&stack) == 0) + break; /* Empty stack => connex component is done! */ + crt_side_id = get_side_from_stack(segsides, &stack); + last_side_id = MMAX(last_side_id, crt_side_id); + } + /* Keep track of this new connex component */ + cc->side_range.last = last_side_id; + /* Need to synchronize connex_components growth as this global structure + * is accessed by multipe threads */ + #pragma omp critical + { + struct cc_descriptor** components; + size_t sz = darray_ptr_component_descriptor_size_get(connex_components); + if(sz <= cc->cc_id) { + res_T tmp_res = darray_ptr_component_descriptor_resize(connex_components, + 1 + cc->cc_id); + if(tmp_res != RES_OK) res = tmp_res; + } + if(res == RES_OK) { + /* Don't set the pointer before resize as this can lead to move data */ + components = + darray_ptr_component_descriptor_data_get(connex_components); + ASSERT(components[cc->cc_id] == NULL); + components[cc->cc_id] = cc; + } + } + OK2(res, error1); + } + continue; + error1: + /* Cannot goto out of openmp block */ + exit_for = 1; + continue; + } + /* No barrier here (nowait clause). + * The first thread executes the single block */ + darray_side_id_release(&stack); + #pragma omp single nowait + if(res == RES_OK) { + struct s2d_device* s2d = NULL; + struct s2d_scene* s2d_scn = NULL; + struct s2d_shape* s2d_shp = NULL; + struct s2d_vertex_data attribs; + + attribs.type = S2D_FLOAT2; + attribs.usage = S2D_POSITION; + attribs.get = get_scn_position; + + /* Put geometry in a 2D view */ + OK(s2d_device_create(desc->scene->dev->logger, alloc, 0, &s2d)); + OK(s2d_scene_create(s2d, &s2d_scn)); + OK(s2d_shape_create_line_segments(s2d, &s2d_shp)); + + /* Back to API type for ntris and nverts */ + ASSERT(desc->scene->nusegs < UINT_MAX); + ASSERT(desc->scene->nuverts < UINT_MAX); + OK(s2d_line_segments_setup_indexed_vertices(s2d_shp, (unsigned)desc->scene->nusegs, + get_scn_indices, (unsigned)desc->scene->nuverts, &attribs, 1, desc->scene)); + s2d_line_segments_set_hit_filter_function(s2d_shp, self_hit_filter, segments_comp); + OK(s2d_scene_attach_shape(s2d_scn, s2d_shp)); + OK(s2d_scene_view_create(s2d_scn, S2D_TRACE, s2d_view)); + error: + if(s2d) S2D(device_ref_put(s2d)); + if(s2d_scn) S2D(scene_ref_put(s2d_scn)); + if(s2d_shp) S2D(shape_ref_put(s2d_shp)); + } + }/* No barrier here */ + OK2(res, error_); + +#ifndef NDEBUG + #pragma omp barrier + ASSERT(component_count == + (int)darray_ptr_component_descriptor_size_get(connex_components)); + FOR_EACH(s_, 0, scn->nusegs) { + struct segment_comp* seg_comp = + darray_segment_comp_data_get(segments_comp) + s_; + ASSERT(seg_comp->component[SIDE_FRONT] != COMPONENT_NULL__); + ASSERT(seg_comp->component[SIDE_BACK] != COMPONENT_NULL__); + } + FOR_EACH(c, 0, component_count) { + struct cc_descriptor** components = + darray_ptr_component_descriptor_data_get(connex_components); + ASSERT(components[c] != NULL && + components[c]->cc_id == c); + } + ASSERT(desc->segment_count + == darray_segment_comp_size_get(segments_comp)); +#endif + +exit: + /* segments_enc is still unused: no size to assert */ + return res; +error_: + goto exit; +} + +static res_T +group_connex_components + (struct senc2d_descriptor* desc, + struct segside* segsides, + struct darray_segment_comp* segments_comp, + struct darray_ptr_component_descriptor* connex_components, + struct s2d_scene_view* s2d_view) +{ + res_T res = RES_OK; + struct cc_descriptor** descriptors; + size_t tmp; + component_id_t cc_count; + int64_t ccc; + static volatile int exit_for = 0; + static ATOMIC next_enclosure_id = 1; + struct cc_descriptor* infinity_first_cc = NULL; + (void)segsides; + ASSERT(desc && segsides && segments_comp && connex_components); + ASSERT(desc->enclosures_count == 1); + + descriptors = darray_ptr_component_descriptor_data_get(connex_components); + tmp = darray_ptr_component_descriptor_size_get(connex_components); + ASSERT(tmp <= COMPONENT_MAX__); + cc_count = (component_id_t)tmp; + + /* Cast rays to find links between connex components */ + #pragma omp /* parallel */ for + for(ccc = 0; ccc < (int64_t)cc_count; ccc++) { + component_id_t c = (component_id_t)ccc; + struct s2d_hit hit = S2D_HIT_NULL; + float origin[2]; + const float dir[2] = { 0, 1 }; + const float range[2] = { 0, FLT_MAX }; + struct cc_descriptor* const cc = descriptors[c]; + const struct segment_comp* origin_seg = + darray_segment_comp_cdata_get(segments_comp) + cc->max_y_vrtx_id; + component_id_t self_hit_component + = origin_seg->component[1 - SEGSIDE_2_SIDE(cc->max_y_side_id)]; + + if(exit_for) continue; + ASSERT(cc->cc_id == c); + ASSERT(cc->cc_group_root == CC_GROUP_ID_NONE); + + if(cc->max_y_ny < 0) { + int64_t id; + /* Don't need to cast a ray */ + cc->cc_group_root = cc->cc_id; /* New group with self as root */ + id = ATOMIC_INCR(&next_enclosure_id) - 1; + ASSERT(id < ENCLOSURE_MAX__); + cc->enclosure_id = (enclosure_id_t)id; + continue; + } + + ASSERT(cc->max_y_ny != 0 + /* The only situation with ny==0 we can think of is this one: */ + || (segsides[cc->max_y_side_id].medium + == segsides[SEGSIDE_OPPOSITE(cc->max_y_side_id)].medium + && (segsides[cc->max_y_side_id].facing_side_id[0] + == SEGSIDE_OPPOSITE(cc->max_y_side_id) + || segsides[cc->max_y_side_id].facing_side_id[1] + == SEGSIDE_OPPOSITE(cc->max_y_side_id)))); + f2_set_d2(origin, cc->max_vrtx); + /* Self-hit data: self hit if hit this component "on the other side" */ + OK2(s2d_scene_view_trace_ray(s2d_view, origin, dir, range, + &self_hit_component, &hit), error_); + /* If no hit, the component is facing an infinite medium */ + if(S2D_HIT_NONE(&hit)) { + cc->cc_group_root = CC_GROUP_ROOT_INFINITE; + cc->enclosure_id = 0; + /* Keep track of the first component facing infinity */ + ATOMIC_CAS_PTR(&infinity_first_cc, cc, NULL); + if(infinity_first_cc->medium != cc->medium) { + const side_id_t infinity_first_side = infinity_first_cc->max_y_side_id; + const medium_id_t infinity_medium = infinity_first_cc->medium; + /* Medium mismatch! Model topology is broken. */ + const seg_id_t t1 = SEGSIDE_2_SEG(infinity_first_side); + const seg_id_t t2 = SEGSIDE_2_SEG(cc->max_y_side_id); + const struct segment_in* segments_in + = darray_segment_in_cdata_get(&desc->scene->segments_in); + const union double2* positions + = darray_position_cdata_get(&desc->scene->vertices); + log_err(desc->scene->dev, + "Medium mismatch found between segment %lu %s side and segment" + " %lu %s side, both facing infinity.\n", + (unsigned long)segments_in[t1].global_id, + SEGSIDE_IS_FRONT(infinity_first_side) ? "front" : "back", + (unsigned long) segments_in[t2].global_id, + SEGSIDE_IS_FRONT(cc->max_y_side_id) ? "front" : "back"); + log_err(desc->scene->dev, + "Segment %lu:\n (%g %g) (%g %g)\n", + (unsigned long)segments_in[t1].global_id, + SPLIT2(positions[segments_in[t1].vertice_id[0]].vec), + SPLIT2(positions[segments_in[t1].vertice_id[1]].vec)); + log_err(desc->scene->dev, + "Segment %lu:\n (%g %g) (%g %g)\n", + (unsigned long)segments_in[t2].global_id, + SPLIT2(positions[segments_in[t2].vertice_id[0]].vec), + SPLIT2(positions[segments_in[t2].vertice_id[1]].vec)); + log_err(desc->scene->dev, "Media: %lu VS %lu\n", + (unsigned long)infinity_medium, (unsigned long)cc->medium); + res = RES_BAD_ARG; + goto error_; + } + /* Same medium as previous members of the group: OK */ + continue; + } else { + /* If hit, group this component */ + const seg_id_t hit_seg_id = (seg_id_t)hit.prim.prim_id; + const struct segment_in* hit_seg_in = + darray_segment_in_cdata_get(&desc->scene->segments_in) + hit_seg_id; + const struct segment_comp* hit_seg_comp = + darray_segment_comp_cdata_get(segments_comp) + hit_seg_id; + enum side_id hit_side = (hit.normal[1] > 0) ? SIDE_FRONT : SIDE_BACK; + const side_id_t hit_side_id = SEGIDxSIDE_2_SEGSIDE(hit_seg_id, hit_side); + + ASSERT(hit_seg_id < desc->scene->nusegs); + + /* Not really the root until following links */ + cc->cc_group_root = hit_seg_comp->component[hit_side]; +#ifndef NDEBUG + { + const struct cc_descriptor* hit_desc; + ASSERT(cc->cc_group_root + < darray_ptr_component_descriptor_size_get(connex_components)); + hit_desc = *(darray_ptr_component_descriptor_cdata_get(connex_components) + + cc->cc_group_root); + ASSERT(hit_desc->medium == hit_seg_in->medium[hit_side]); + } +#endif + if(hit_seg_in->medium[hit_side] != cc->medium) { + /* Medium mismatch! Model topology is broken. */ + const seg_id_t t1 = SEGSIDE_2_SEG(hit_side_id); + const seg_id_t t2 = SEGSIDE_2_SEG(cc->max_y_side_id); + const struct segment_in* segments_in + = darray_segment_in_cdata_get(&desc->scene->segments_in); + const union double2* positions + = darray_position_cdata_get(&desc->scene->vertices); + log_err(desc->scene->dev, + "Medium mismatch found between segment %lu %s side and segment" + " %lu %s side facing each other.\n", + (unsigned long)segments_in[t1].global_id, + SEGSIDE_IS_FRONT(hit_side) ? "front" : "back", + (unsigned long)segments_in[t2].global_id, + SEGSIDE_IS_FRONT(cc->max_y_side_id) ? "front" : "back"); + log_err(desc->scene->dev, + "Segment %lu:\n (%g %g) (%g %g)\n", + (unsigned long)segments_in[t1].global_id, + SPLIT2(positions[segments_in[t1].vertice_id[0]].vec), + SPLIT2(positions[segments_in[t1].vertice_id[1]].vec)); + log_err(desc->scene->dev, + "Segment %lu:\n (%g %g) (%g %g)\n", + (unsigned long)segments_in[t2].global_id, + SPLIT2(positions[segments_in[t2].vertice_id[0]].vec), + SPLIT2(positions[segments_in[t2].vertice_id[1]].vec)); + log_err(desc->scene->dev, "Media: %lu VS %lu\n", + (unsigned long)hit_seg_in->medium[hit_side], + (unsigned long)cc->medium); + + res = RES_BAD_ARG; + goto error_; + } + } + continue; + error_: + /* Cannot goto out of openmp block */ + exit_for = 1; + continue; + } /* Implicit barrier here */ + ASSERT(next_enclosure_id < ENCLOSURE_MAX__); + OK(res); + + /* Post-process links to group connex components */ + #pragma omp single + { + desc->enclosures_count = (enclosure_id_t)next_enclosure_id; + res = darray_enclosure_resize(&desc->enclosures, desc->enclosures_count); + if(res == RES_OK) { + FOR_EACH(ccc, 0, cc_count) { + component_id_t c = (component_id_t) ccc; + struct cc_descriptor* const cc = descriptors[c]; + const struct cc_descriptor* other_desc = cc; + struct enclosure_data* enclosures + = darray_enclosure_data_get(&desc->enclosures); + component_id_t fst; + + while(other_desc->enclosure_id == CC_GROUP_ID_NONE) { + other_desc = *(darray_ptr_component_descriptor_cdata_get(connex_components) + + other_desc->cc_group_root); + } + ASSERT(other_desc->cc_group_root != CC_GROUP_ROOT_NONE); + ASSERT(other_desc->enclosure_id != CC_GROUP_ID_NONE); + ASSERT(cc->medium == other_desc->medium); + cc->cc_group_root = other_desc->cc_group_root; + cc->enclosure_id = other_desc->enclosure_id; + ++enclosures[cc->enclosure_id].cc_count; + /* Linked list of componnents */ + fst = enclosures[cc->enclosure_id].first_component; + cc->enclosure_next_component = fst; + enclosures[cc->enclosure_id].first_component = cc->cc_id; + enclosures[cc->enclosure_id].side_range.first + = MMIN(enclosures[cc->enclosure_id].side_range.first, cc->side_range.first); + enclosures[cc->enclosure_id].side_range.last + = MMAX(enclosures[cc->enclosure_id].side_range.last, cc->side_range.last); + enclosures[cc->enclosure_id].side_count += cc->side_count; + } + } + } /* Implicit barrier here */ + OK(res); + +exit: + return res; +error: + goto exit; +} + +static res_T +collect_and_link_neighbours + (struct senc2d_scene* scn, + struct segside* segsides, + struct darray_segment_tmp* segments_tmp_array, + struct darray_neighbourhood* neighbourhood_by_vertex) +{ + res_T res = RES_OK; + const struct segment_in *segments_in; + struct segment_tmp *segments_tmp; + const union double2* vertices; + static volatile int exit_for = 0; + + ASSERT(scn && segsides && segments_tmp_array); + ASSERT((size_t)scn->nuverts + (size_t)scn->nusegs + 2 <= EDGE_MAX__); + + segments_in = darray_segment_in_cdata_get(&scn->segments_in); + segments_tmp = darray_segment_tmp_data_get(segments_tmp_array); + vertices = darray_position_cdata_get(&scn->vertices); + + ASSERT(scn->nusegs == darray_segment_tmp_size_get(segments_tmp_array)); + + /* Resize / Push operations on neighbourhood_by_vertex are valid in the + * openmp block because each neighbourhood is processes by an unique thread */ + /* #pragma omp parallel */ + { + const int thread_count = omp_get_num_threads(); + const int rank = omp_get_thread_num(); + /* Array to keep neighbourhood of vertices + * Resize/Push operations on neighbourhood_by_vertex are valid in the + * openmp block because it is thread local data */ + vrtx_id_t v; + seg_id_t s; + /* Loop on segments to collect edges' neighbours. + * All threads considering all the vertices and processing some */ + FOR_EACH(s, 0, scn->nusegs) { + unsigned char vv; + FOR_EACH(vv, 0, 2) { + struct darray_neighbour* neighbourhood; + struct neighbour_info* info; + const vrtx_id_t vertex = segments_in[s].vertice_id[vv]; + size_t sz; + ASSERT(vertex < scn->nuverts); + if(exit_for) continue; + /* Process only "my" vertices! */ + if((int64_t)vertex % thread_count != rank) continue; + /* Find neighbourhood */ + ASSERT(vertex < darray_neighbourhood_size_get(neighbourhood_by_vertex)); + neighbourhood = + darray_neighbourhood_data_get(neighbourhood_by_vertex) + vertex; + sz = darray_neighbour_size_get(neighbourhood); + /* Make room for this neighbour */ + if(darray_neighbour_capacity(neighbourhood) == sz) { + /* 2 seems to be a good guess for initial capacity */ + size_t new_sz = sz ? sz + 1 : 2; + OK(darray_neighbour_reserve(neighbourhood, new_sz)); + } + OK(darray_neighbour_resize(neighbourhood, 1 + sz)); + /* Add neighbour info to vertex's neighbour list */ + info = darray_neighbour_data_get(neighbourhood) + sz; + info->seg_id = s; + info->common_vertex_rank = vv; + } + } /* No barrier here. */ + + /* For each of "my" vertices sort segments sides by rotation angle + * and connect neighbours. + * All threads considering all the vertices and processing some */ + FOR_EACH(v, 0, scn->nuverts) { + const vrtx_id_t common_vrtx = v; + vrtx_id_t other_vrtx; + struct darray_neighbour* neighbourhood; + side_id_t i, neighbour_count; + size_t sz; + /* Process only "my" vertices! */ + if((int64_t)v % thread_count != rank) continue; + neighbourhood + = darray_neighbourhood_data_get(neighbourhood_by_vertex) + v; + sz = darray_neighbour_size_get(neighbourhood); + /* sz can be 0 as a vertex can be unused */ + if(!sz) continue; + ASSERT(sz <= SIDE_MAX__); + neighbour_count = (side_id_t)sz; + FOR_EACH(i, 0, neighbour_count) { + double max_y, disp[2]; + unsigned char max_y_vrank; + struct neighbour_info* neighbour_info + = darray_neighbour_data_get(neighbourhood) + i; + const struct segment_in* seg_in = segments_in + neighbour_info->seg_id; + struct segment_tmp* neighbour = segments_tmp + neighbour_info->seg_id; + other_vrtx = + seg_in->vertice_id[(neighbour_info->common_vertex_rank + 1) % 2]; + if(vertices[other_vrtx].pos.y > vertices[common_vrtx].pos.y) { + max_y = vertices[other_vrtx].pos.y; + max_y_vrank = 1 - neighbour_info->common_vertex_rank; + } else { + max_y = vertices[common_vrtx].pos.y; + max_y_vrank = neighbour_info->common_vertex_rank; + } + ASSERT(neighbour->max_y <= max_y); + neighbour->max_y = max_y; + neighbour->max_y_vrtx_rank = max_y_vrank; + /* Compute rotation angle around common vertex (in world system) */ + d2_sub(disp, vertices[other_vrtx].vec, vertices[common_vrtx].vec); + ASSERT(disp[0] || disp[1]); + neighbour_info->angle = atan2(disp[1], disp[0]); + if(neighbour_info->angle < 0) neighbour_info->angle += 2 * PI; + /* Due to catastrophic cancelation, -eps+2pi translates to 2pi */ + ASSERT(0 <= neighbour_info->angle && neighbour_info->angle <= 2 * PI); + } + /* Sort segments by rotation angle */ + qsort(darray_neighbour_data_get(neighbourhood), neighbour_count, + sizeof(struct neighbour_info), neighbour_cmp); + /* Link sides. + * Create cycles of sides by neighbourhood around common vertex. */ + FOR_EACH(i, 0, neighbour_count) { + /* Neighbourhood info for current pair of segments */ + const struct neighbour_info* current + = darray_neighbour_cdata_get(neighbourhood) + i; + const struct neighbour_info* ccw_neighbour + = darray_neighbour_cdata_get(neighbourhood) + (i+1) % neighbour_count; + /* Rank of the end of interest in segments */ + const unsigned char crt_end = current->common_vertex_rank; + const unsigned char ccw_end = ccw_neighbour->common_vertex_rank; + /* User id of current segments */ + const seg_id_t crt_id = current->seg_id; + const seg_id_t ccw_id = ccw_neighbour->seg_id; + /* Facing sides of segments */ + const enum side_id crt_side = crt_end ? SIDE_BACK : SIDE_FRONT; + const enum side_id ccw_side = ccw_end ? SIDE_FRONT : SIDE_BACK; + /* Index of sides in segsides */ + const side_id_t crt_side_idx = SEGIDxSIDE_2_SEGSIDE(crt_id, crt_side); + const side_id_t ccw_side_idx = SEGIDxSIDE_2_SEGSIDE(ccw_id, ccw_side); + /* Side ptrs */ + struct segside* const p_crt_side = segsides + crt_side_idx; + struct segside* const p_ccw_side = segsides + ccw_side_idx; + /* Link sides */ + p_crt_side->facing_side_id[crt_end] = ccw_side_idx; + p_ccw_side->facing_side_id[ccw_end] = crt_side_idx; + /* Record media */ + p_crt_side->medium = segments_in[crt_id].medium[crt_side]; + p_ccw_side->medium = segments_in[ccw_id].medium[ccw_side]; + ASSERT(p_crt_side->medium < scn->nmeds); + ASSERT(p_ccw_side->medium < scn->nmeds); + p_crt_side->list_id = FLAG_LIST_SIDE_LIST; + p_ccw_side->list_id = FLAG_LIST_SIDE_LIST; + } + } /* No barrier here. */ + /* jump error block */ + goto after_error; +error: + /* Cannot goto out of openmp block */ + exit_for = 1; +after_error: + ; + } + return res; +} + +static res_T +build_result + (struct senc2d_descriptor* desc, + const struct darray_ptr_component_descriptor* connex_components, + const struct darray_segment_comp* segments_comp_array) +{ + res_T res = RES_OK; + struct mem_allocator* alloc; + struct cc_descriptor* const* cc_descriptors; + struct enclosure_data* enclosures; + const struct segment_in* segments_in; + struct segment_enc* segments_enc; + const struct segment_comp* segments_comp; + static volatile int exit_for = 0; + + ASSERT(desc && connex_components && segments_comp_array); + + alloc = descriptor_get_allocator(desc); + ASSERT(darray_ptr_component_descriptor_size_get(connex_components) < COMPONENT_MAX__); + cc_descriptors = darray_ptr_component_descriptor_cdata_get(connex_components); + enclosures = darray_enclosure_data_get(&desc->enclosures); + segments_in = darray_segment_in_cdata_get(&desc->scene->segments_in); + segments_comp = darray_segment_comp_cdata_get(segments_comp_array); + #pragma omp single + { + res = darray_segment_enc_resize(&desc->segments_enc, desc->scene->nusegs); + }/* Implicit barrier here. */ + OK2(res, error_); + segments_enc = darray_segment_enc_data_get(&desc->segments_enc); + + /* #pragma omp parallel */ + { + struct htable_vrtx_id vtable; + int64_t sg; + int64_t ee; + + /* Build global enclosure information */ + #pragma omp for + for(sg = 0; sg < (int64_t) desc->scene->nusegs; sg++) { + seg_id_t s = (seg_id_t) sg; + const component_id_t cf_id = segments_comp[s].component[SIDE_FRONT]; + const component_id_t cb_id = segments_comp[s].component[SIDE_BACK]; + const struct cc_descriptor* cf = cc_descriptors[cf_id]; + const struct cc_descriptor* cb = cc_descriptors[cb_id]; + const enclosure_id_t ef_id = cf->enclosure_id; + const enclosure_id_t eb_id = cb->enclosure_id; + ASSERT(segments_enc[s].enclosure[SIDE_FRONT] == ENCLOSURE_NULL__); + segments_enc[s].enclosure[SIDE_FRONT] = ef_id; + ASSERT(segments_enc[s].enclosure[SIDE_BACK] == ENCLOSURE_NULL__); + segments_enc[s].enclosure[SIDE_BACK] = eb_id; + } + /* Implicit barrier here */ + + /* Resize/push operations on enclosure's fields are valid in the + * openmp block as a given enclosure is processed by a single thread */ + htable_vrtx_id_init(alloc, &vtable); + + ASSERT(desc->enclosures_count <= ENCLOSURE_MAX__); + #pragma omp for schedule(dynamic) nowait + for(ee = 0; ee < (int64_t)desc->enclosures_count; ee++) { + const enclosure_id_t e = (enclosure_id_t)ee; + struct enclosure_data* enc = enclosures + e; + const struct cc_descriptor* current = cc_descriptors[enc->first_component]; + seg_id_t fst_idx = 0; + seg_id_t sgd_idx = enc->side_count; + seg_id_t s; + ASSERT(enc->first_component + < darray_ptr_component_descriptor_size_get(connex_components)); + ASSERT(current->cc_id == enc->first_component); + + if(exit_for) continue; + ASSERT(e <= UINT_MAX); + enc->header.enclosure_id = (unsigned)e; /* Back to API type */ + ASSERT(current->enclosure_id == enc->header.enclosure_id); + enc->header.is_infinite = (e == 0); + enc->header.enclosed_medium + = (unsigned)current->medium; /* Back to API type */ + ASSERT(enc->header.enclosed_medium < desc->scene->nmeds); + + /* Build side and vertex lists. */ + OK(darray_segment_in_resize(&enc->sides, enc->side_count)); + /* Size is just a int */ + OK(darray_vrtx_id_reserve(&enc->vertices, enc->side_count + 1)); + /* New vertex numbering scheme local to the enclosure */ + htable_vrtx_id_clear(&vtable); + ASSERT(desc->scene->nusegs + == darray_segment_in_size_get(&desc->scene->segments_in)); + /* Put at the end the back-faces of segments that also have their + * front-face in the list. */ + for(s = SEGSIDE_2_SEG(enc->side_range.first); + s <= SEGSIDE_2_SEG(enc->side_range.last); + s++) + { + const struct segment_in* seg_in = segments_in + s; + struct segment_in* seg; + unsigned vertice_id[2]; + int i; + if(segments_enc[s].enclosure[SIDE_FRONT] != e + && segments_enc[s].enclosure[SIDE_BACK] != e) + continue; + ++enc->header.unique_segment_count; + + FOR_EACH(i, 0, 2) { + vrtx_id_t* id = htable_vrtx_id_find(&vtable, seg_in->vertice_id + i); + if(id) { + vertice_id[i] = *id; /* Known vertex */ + } else { + /* Create new association */ + size_t tmp = htable_vrtx_id_size_get(&vtable); + ASSERT(tmp == darray_vrtx_id_size_get(&enc->vertices)); + ASSERT(tmp < VRTX_MAX__); + vertice_id[i] = (vrtx_id_t)tmp; + OK(htable_vrtx_id_set(&vtable, seg_in->vertice_id + i, + vertice_id + i)); + OK(darray_vrtx_id_push_back(&enc->vertices, seg_in->vertice_id + i)); + ++enc->header.vertices_count; + } + } + ASSERT(segments_enc[s].enclosure[SIDE_FRONT] == e + || segments_enc[s].enclosure[SIDE_BACK] == e); + if(segments_enc[s].enclosure[SIDE_FRONT] == e) { + ++enc->header.segment_count; + seg = darray_segment_in_data_get(&enc->sides) + fst_idx++; + + FOR_EACH(i, 0, 2) seg->medium[i] = seg_in->medium[i]; + seg->global_id = seg_in->global_id; + FOR_EACH(i, 0, 2) seg->vertice_id[i] = vertice_id[i]; + } + if(segments_enc[s].enclosure[SIDE_BACK] == e) { + ++enc->header.segment_count; + seg = darray_segment_in_data_get(&enc->sides) + + ((segments_enc[s].enclosure[SIDE_FRONT] == e) ? --sgd_idx : fst_idx++); + + FOR_EACH(i, 0, 2) seg->medium[i] = seg_in->medium[1 - i]; + seg->global_id = seg_in->global_id; + FOR_EACH(i, 0, 2) seg->vertice_id[i] = vertice_id[1 - i]; + } + if(fst_idx == sgd_idx) break; + } + continue; + error: + /* Cannot goto out of openmp block */ + exit_for = 1; + } /* No barrier here */ + htable_vrtx_id_release(&vtable); + } + OK2(res, error_); + +exit: + return res; +error_: + goto exit; +} + +/******************************************************************************* + * Exported functions + ******************************************************************************/ +res_T +senc2d_scene_analyze(struct senc2d_scene* scn, struct senc2d_descriptor** out_desc) +{ + res_T res = RES_OK; + struct senc2d_descriptor* desc = NULL; + /* By segment tmp data */ + struct darray_segment_tmp segments_tmp; + char segments_tmp_initialized = 0; + /* Array of connex components. + * They are refered to by arrays of ids. */ + struct darray_ptr_component_descriptor connex_components; + char connex_components_initialized = 0; + /* Store by-segment components */ + struct darray_segment_comp segments_comp; + char segments_comp_initialized = 0; + /* Array of vertices (segment ends). */ + struct segside* segsides = NULL; + struct s2d_scene_view* s2d_view = NULL; + /* Array to keep neighbourhood of vertices */ + struct darray_neighbourhood neighbourhood_by_vertex; + char neighbourhood_by_vertex_initialized = 0; + + if(!scn || !out_desc) return RES_BAD_ARG; + + desc = descriptor_create(scn); + if(!desc) { + res = RES_MEM_ERR; + goto error; + } + + if(!scn->nusegs) goto exit; + + darray_segment_tmp_init(scn->dev->allocator, &segments_tmp); + segments_tmp_initialized = 1; + + OK(darray_segment_tmp_resize(&segments_tmp, scn->nusegs)); + segsides + = MEM_CALLOC(scn->dev->allocator, 2 * scn->nusegs, sizeof(struct segside)); + if(!segsides) { + res = RES_MEM_ERR; + goto error; + } + + /* We create a neighbourhood for every single unique vertex, + * regardless the fact it is used by a segment */ + darray_neighbourhood_init(scn->dev->allocator, &neighbourhood_by_vertex); + neighbourhood_by_vertex_initialized = 1; + OK(darray_neighbourhood_resize(&neighbourhood_by_vertex, scn->nuverts)); + + #pragma omp parallel + { + res_T tmpres = RES_OK; + + /* Step 1: build neighbourhoods */ + tmpres = collect_and_link_neighbours(scn, segsides, &segments_tmp, + &neighbourhood_by_vertex); + /* No barrier at the end of step 1 */ + + if(tmpres != RES_OK) { + res = tmpres; + #pragma omp single nowait + { + log_err(scn->dev, + "%s: could not build neighbourhoods from scene.\n", FUNC_NAME); + } /* No barrier here */ + } + OK2(res, error_); + + /* The first thread here allocates some data */ + #pragma omp single + { + darray_ptr_component_descriptor_init(scn->dev->allocator, + &connex_components); + connex_components_initialized = 1; + /* Just a hint; to limit contention */ + tmpres = darray_ptr_component_descriptor_reserve(&connex_components, + 2 * scn->nmeds); + if(tmpres != RES_OK) res = tmpres; + darray_segment_comp_init(scn->dev->allocator, &segments_comp); + segments_comp_initialized = 1; + tmpres = darray_segment_comp_resize(&segments_comp, scn->nusegs); + if(tmpres != RES_OK) res = tmpres; + } /* Implicit barrier here */ + OK2(res, error_); + + /* One thread releases some data, the others go to step 2 */ + #pragma omp single nowait + { + darray_neighbourhood_release(&neighbourhood_by_vertex); + neighbourhood_by_vertex_initialized = 0; + } /* No implicit barrier here */ + + + /* Step 2: extract segment connex components */ + tmpres = extract_connex_components(desc, segsides, &connex_components, + &segments_tmp, &segments_comp, &s2d_view); + /* No barrier at the end of step 2 */ + + if(tmpres != RES_OK) { + res = tmpres; + #pragma omp single nowait + { + log_err(scn->dev, + "%s: could not extract connex components from scene.\n", FUNC_NAME); + } /* No implicit barrier here */ + } + OK2(res, error_); + + #pragma omp barrier + /* One thread releases some data, the others go to step 3 */ + #pragma omp single nowait + { + darray_segment_tmp_release(&segments_tmp); + segments_tmp_initialized = 0; + } /* No implicit barrier here */ + + /* Step 3: group components */ + tmpres = group_connex_components(desc, segsides, &segments_comp, + &connex_components, s2d_view); + /* Barrier at the end of step 3 */ + + /* The first thread here releases some data */ + #pragma omp single nowait + { + if(s2d_view) S2D(scene_view_ref_put(s2d_view)); + s2d_view = NULL; + } /* No implicit barrier here */ + + if(res != RES_OK) { + res = tmpres; + #pragma omp single nowait + { + log_err(scn->dev, + "%s: could not group connex components from scene.\n", FUNC_NAME); + } /* No implicit barrier here */ + } + OK2(res, error_); + + /* Build result. */ + tmpres = build_result(desc, &connex_components, &segments_comp); + /* No barrier at the end of build step */ + + if(tmpres != RES_OK) { + res = tmpres; + #pragma omp single nowait + { + log_err(scn->dev, "%s: could not build result.\n", FUNC_NAME); + } /* No implicit barrier here */ + } + OK2(res, error_); + + #pragma omp barrier + /* Some threads release data */ + #pragma omp sections nowait + { + #pragma omp section + { + darray_segment_comp_release(&segments_comp); + segments_comp_initialized = 0; + } + #pragma omp section + { + if(connex_components_initialized) { + size_t c, cc_count = + darray_ptr_component_descriptor_size_get(&connex_components); + struct cc_descriptor** components = + darray_ptr_component_descriptor_data_get(&connex_components); + FOR_EACH(c, 0, cc_count) { + ptr_component_descriptor_release(scn->dev->allocator, components + c); + } + darray_ptr_component_descriptor_release(&connex_components); + connex_components_initialized = 0; + } + } + } /* No implicit barrier here */ +error_: + ; + } /* Implicit barrier here */ + +if(res != RES_OK) goto error; +exit: + if(connex_components_initialized) { + size_t c, cc_count = + darray_ptr_component_descriptor_size_get(&connex_components); + struct cc_descriptor** components = + darray_ptr_component_descriptor_data_get(&connex_components); + FOR_EACH(c, 0, cc_count) { + ptr_component_descriptor_release(scn->dev->allocator, components + c); + } + darray_ptr_component_descriptor_release(&connex_components); + } + if(s2d_view) S2D(scene_view_ref_put(s2d_view)); + if(neighbourhood_by_vertex_initialized) + darray_neighbourhood_release(&neighbourhood_by_vertex); + if(segments_tmp_initialized) darray_segment_tmp_release(&segments_tmp); + if(segments_comp_initialized) darray_segment_comp_release(&segments_comp); + if(segsides) MEM_RM(scn->dev->allocator, segsides); + if(desc) *out_desc = desc; + + return res; +error: + if(desc) SENC2D(descriptor_ref_put(desc)); + desc = NULL; + goto exit; +} diff --git a/src/senc2d_scene_analyze_c.h b/src/senc2d_scene_analyze_c.h @@ -105,20 +105,27 @@ ptr_component_descriptor_init ASSERT(data); *data = NULL; } -static FINLINE void -ptr_component_descriptor_release - (struct mem_allocator* allocator, - struct cc_descriptor** data) -{ - ASSERT(allocator && data); - MEM_RM(allocator, *data); -} #define DARRAY_NAME ptr_component_descriptor #define DARRAY_DATA struct cc_descriptor* #define DARRAY_FUNCTOR_INIT ptr_component_descriptor_init #include <rsys/dynamic_array.h> +/* Need allocator to free array elts: cannot rely on standard + * darray release stuff */ +static FINLINE void +custom_darray_ptr_component_descriptor_release + (struct darray_ptr_component_descriptor* array) +{ + size_t c, cc_count; + struct cc_descriptor** components; + if(!array) return; + cc_count = darray_ptr_component_descriptor_size_get(array); + components = darray_ptr_component_descriptor_data_get(array); + FOR_EACH(c, 0, cc_count) MEM_RM(array->allocator, components[c]); + darray_ptr_component_descriptor_release(array); +} + /* Segment information. * Depending on lifespan, information is kept in different places: * - segment_in for user provided information (kept in scene)