path: root/components.c

/*
 * 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;
}