diff options
Diffstat (limited to 'components.c')
-rw-r--r-- | components.c | 251 |
1 files changed, 251 insertions, 0 deletions
diff --git a/components.c b/components.c new file mode 100644 index 0000000..ef41653 --- /dev/null +++ b/components.c @@ -0,0 +1,251 @@ +/* + * Multi-component utility code + * Copyright 2019 Anton Khirnov <anton@khirnov.net> + * + * 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 <errno.h> +#include <stddef.h> +#include <stdint.h> +#include <stdlib.h> +#include <string.h> + +#include "common.h" +#include "components.h" +#include "mg2d_boundary.h" + +int mg2di_dg_copy(DomainGeometry **pdst, const DomainGeometry *src) +{ + DomainGeometry *dst; + + dst = mg2di_dg_alloc(src->nb_components); + if (!dst) + return -ENOMEM; + + memcpy(dst->components, src->components, + src->nb_components * sizeof(*src->components)); + dst->domain_size[0] = src->domain_size[0]; + dst->domain_size[1] = src->domain_size[1]; + + *pdst = dst; + + return 0; +} + +DomainGeometry *mg2di_dg_alloc(unsigned int nb_components) +{ + DomainGeometry *dg; + + dg = calloc(1, sizeof(*dg)); + if (!dg) + return NULL; + + dg->components = calloc(nb_components, sizeof(*dg->components)); + if (!dg->components) + goto fail; + dg->nb_components = nb_components; + + return dg; +fail: + mg2di_dg_free(&dg); + return NULL; +} + +void mg2di_dg_free(DomainGeometry **pdg) +{ + DomainGeometry *dg = *pdg; + + if (!dg) + return; + + free(dg->components); + + free(dg); + *pdg = NULL; +} + +int mg2di_rect_intersect(Rect *dst, const Rect *src1, const Rect *src2) +{ + ptrdiff_t intersect_start0 = MAX(src1->start[0], src2->start[0]); + ptrdiff_t intersect_start1 = MAX(src1->start[1], src2->start[1]); + ptrdiff_t intersect_end0 = MIN(src1->start[0] + src1->size[0], src2->start[0] + src2->size[0]); + ptrdiff_t intersect_end1 = MIN(src1->start[1] + src1->size[1], src2->start[1] + src2->size[1]); + + if (intersect_start0 < intersect_end0 && intersect_start1 < intersect_end1) { + dst->start[0] = intersect_start0; + dst->start[1] = intersect_start1; + dst->size[0] = intersect_end0 - intersect_start0; + dst->size[1] = intersect_end1 - intersect_start1; + + return 1; + } + + dst->size[0] = 0; + dst->size[1] = 0; + + return 0; +} + +/* merge adjacent/partially overlapping src into dst, + * if the result is rectangular */ +static int rect_merge(Rect *dst, const Rect *src) +{ + /* if dst is empty, copy src */ + if (!dst->size[0] && !dst->size[1]) { + *dst = *src; + return 0; + } + + /* if dst is fully inside src, copy src */ + if (dst->start[0] >= src->start[0] && + dst->start[0] + dst->size[0] <= src->start[0] + src->size[0] && + dst->start[1] >= src->start[1] && + dst->start[1] + dst->size[1] <= src->start[1] + src->size[1]) { + *dst = *src; + return 0; + } + + /* if src is fully inside dst, do nothing */ + if (src->start[0] >= dst->start[0] && + src->start[0] + src->size[0] <= dst->start[0] + dst->size[0] && + src->start[1] >= dst->start[1] && + src->start[1] + src->size[1] <= dst->start[1] + dst->size[1]) + return 0; + + /* if src is adjacent or partially overlaps dst, merge them */ + for (int dir = 0; dir < 2; dir++) { + if (dst->start[dir] != src->start[dir] || + dst->size[dir] != src->size[dir]) + continue; + + /* merge from above */ + if (src->start[!dir] >= dst->start[!dir] && + src->start[!dir] <= dst->start[!dir] + dst->size[!dir]) { + dst->size[!dir] = MAX(dst->size[!dir], src->start[!dir] + src->size[!dir] - dst->start[!dir]); + return 0; + } + + /* merge from below */ + if (dst->start[!dir] >= src->start[!dir] && + dst->start[!dir] <= src->start[!dir] + src->size[!dir]) { + + dst->size[!dir] = MAX(src->size[!dir], dst->start[!dir] + dst->size[!dir] - src->start[!dir]); + dst->start[!dir] = src->start[!dir]; + return 0; + } + } + + /* src and dst do not match */ + return -EINVAL; +} + +static int calc_overlaps_component(Rect *overlaps, const DomainGeometry *dg, + unsigned int comp, int edge_width) +{ + const DomainComponent *c = &dg->components[comp]; + int ret = 0; + + Rect (*overlaps_tmp)[4]; + int *nb_overlaps_tmp; + + overlaps_tmp = calloc(dg->nb_components, sizeof(*overlaps_tmp)); + nb_overlaps_tmp = calloc(dg->nb_components, sizeof(*nb_overlaps_tmp)); + if (!overlaps_tmp || !nb_overlaps_tmp) + return -ENOMEM; + + /* calculate the overlaps for each boundary layer */ + for (int bnd_idx = 0; bnd_idx < 4; bnd_idx++) { + const int ci = mg2d_bnd_coord_idx(bnd_idx); + const int upper = mg2d_bnd_is_upper(bnd_idx); + Rect bnd; + + bnd.start[ci] = c->interior.start[ci] + (upper ? c->interior.size[ci] : -edge_width); + bnd.size[ci] = edge_width; + + bnd.start[!ci] = c->interior.start[!ci] - edge_width; + bnd.size[!ci] = c->interior.size[!ci] + 2 * edge_width; + + for (unsigned int comp_idx = 0; comp_idx < dg->nb_components; comp_idx++) { + const DomainComponent *c1 = &dg->components[comp_idx]; + Rect *dst = &overlaps_tmp[comp_idx][nb_overlaps_tmp[comp_idx]]; + + if (comp_idx == comp) + continue; + + mg2di_rect_intersect(dst, &bnd, &c1->exterior); + if (dst->size[0] && dst->size[1]) + nb_overlaps_tmp[comp_idx]++; + } + } + + /* merge all the overlaps, we should get a single rectangle for each component */ + for (unsigned int comp_idx = 0; comp_idx < dg->nb_components; comp_idx++) { + while (nb_overlaps_tmp[comp_idx] > 1) { + Rect *dst = &overlaps_tmp[comp_idx][0]; + int merged = 0; + + for (int i = 1; i < nb_overlaps_tmp[comp_idx]; i++) { + Rect *src = &overlaps_tmp[comp_idx][i]; + ret = rect_merge(dst, src); + if (ret < 0) + continue; + + memmove(src, src + 1, nb_overlaps_tmp[comp_idx] - i - 1); + nb_overlaps_tmp[comp_idx]--; + merged = 1; + break; + } + + if (!merged) + goto fail; + } + overlaps[comp_idx] = overlaps_tmp[comp_idx][0]; + } + +fail: + free(overlaps_tmp); + free(nb_overlaps_tmp); + + return ret; +} + +int mg2di_dg_edge_overlaps(Rect *overlaps_recv, Rect *overlaps_send, + const DomainGeometry *dg, unsigned int comp, + unsigned int edge_width) +{ + Rect *overlaps = NULL; + int ret; + + overlaps = calloc(SQR(dg->nb_components), sizeof(*overlaps)); + if (!overlaps) + return -ENOMEM; + + for (unsigned int i = 0; i < dg->nb_components; i++) { + ret = calc_overlaps_component(overlaps + i * dg->nb_components, + dg, i, edge_width); + if (ret < 0) + goto fail; + } + + memcpy(overlaps_recv, overlaps + comp * dg->nb_components, + dg->nb_components * sizeof(*overlaps_recv)); + + for (unsigned int i = 0; i < dg->nb_components; i++) + overlaps_send[i] = overlaps[i * dg->nb_components + comp]; + +fail: + free(overlaps); + return ret; +} |