function [cmpx, newset] = mfs_newset(cmp, settype, fct, ops, setnam)

# usage: [cmpx, newset] = mfs_newset(cmp, settype, fct, ops, setnam)
#
# Input  cmp      Structure with component data
#        settype  Type of set (nset or eset)
#        fct      Function to be performed
#        ops      Operands
#        setnam   Name of new set
#
# Output cmpx     Structure with component data
#        newset   Array with identifiers of objects in new set
#                 (optional)
#
# The function creates a new set by applying the requested operation
# fct. The new set is stored in the component. If two output arguments
# are specified, an array with identifiers of # objects in the new set
# is also returned.
#
# The following functions are supported:
#
#   "near"        node or element closest to point
#                   ops(:)  Coordinates of point
#   "box"         nodes or elements inside a box
#                   ops(:, :)  xmin, xmax; ymin, ymax; zmin, zmax
#   "union"       union of the sets:
#                   ops{:}  Cell array of names of sets
#   "intersect"   intersection of the sets
#                   ops{:}  Cell array of names of sets
#
# --------------------------------------------------------------------

# Copyright (c) 2023 by Johannes Wandinger

  newset = [];
  cmpx   = cmp;
  loop   = 0;

# Check arguments

  if (nargin != 5 || nargout < 1 || nargout > 2)
     print_usage();
  end
  if (! isstruct(cmp))
     error("mfs_newset: first argument must be a structure\n");
  end
  if (! ischar(settype))
     error("mfs_newset: settype must be a string\n");
  end
  if (! ischar(fct))
     error("mfs_newset: set function must be a string\n");
  end
  if (! ischar(setnam))
     error("mfs_newset: setname must be a string\n");
  end

# Get set type

  nset = strcmp(settype, "nset");
  if (! nset && ! strcmp(settype, "eset"))
     error("mfs_newset: unknown settype %s\n", settype);
  end

# Check if new set already exists

  setold = cmp.(settype);
  if (isfield(setold, setnam))
     error("mfs_newset: set %s already exists\n", setnam);
  end

# Get set function

  switch fct

  case "near"

     if (! isnumeric(ops))
        error("mfs_newset: ops must be a vector\n");
     end
     refpnt = ops(:)';
     if (columns(refpnt) != cmp.nodes.ncoor)
        error("mfs_newset: reference point has wrong number of coordinates");
     end
     if (nset)
        d = cmp.nodes.coor - refpnt;
        r = dot(d, d, 2);
        [rmin, setn] = min(r);
     else
        d = eltcoor(cmp.elements.elem, cmp.nodes.ncoor) - refpnt;
        r = dot(d, d, 2);
        [rmin, ixe] = min(r);
        setn = lookup(cmp.elements.index(:, 2), ixe);
     end

  case "box"

     if (! isnumeric(ops))
        error("mfs_newset: ops must be an array\n");
     end
     [mr, nc] = size(ops); 
     ncoor    = cmp.nodes.ncoor;
     if (mr != ncoor || nc != nc)
        error("mfs_newset: ops must be a %1.0d by 2 array\n", ncoor);
     end
     if (nset)
        setn = fillbox(ops, cmp.nodes.coor);
     else
        coor = eltcoor(cmp.elements.elem, cmp.nodes.ncoor);
        ixe  = fillbox(ops, coor);
        setn = lookup(cmp.elements.index(:, 2), ixe);
     end

  case "union"

     if (! iscell(ops))
        error("mfs_newset: ops must be a cell array\n");
     end
     setf = @union;
     loop = 1;

  case "intersect"

     if (! iscell(ops))
        error("mfs_newset: ops must be a cell array\n");
     end
     setf = @intersect;
     loop = 1;

  otherwise

     error("mfs_newset: function %s not supported\n", fct);

  end

# Loop over the sets

  if (loop)

     if (isfield(setold, ops{1}))
        setn = getfield(setold, ops{1});
     else
        error("mfs_newset: set %s does not exist\n", ops{1});
     end

     for n = 2 : length(ops)
         if (isfield(setold, ops{n}))
            setn = setf(setn, getfield(setold, ops{n}), "stable");
         else
            error("mfs_newset: set %s does not exist\n", ops{n});
         end
     end

  end

# Store the new set in the component

  cmpx.(settype) = setfield(setold, setnam, setn);

# Get external identifiers of objects in new set

  if (nargout == 2)
     if (nset)
        newset = cmp.nodes.ids(setn);
     else
        newset = cmp.elements.index(setn, 1);
     end
  end

end

function mcoor = eltcoor(elem, nc)

# Input  elem(:)      Structure array with element data
#        nc           Number of coordinates
#
# Output mcoor(nc, :) Coordinates of element mean points
#
# -------------------------------------------------------------------

  nofelt = length(elem);
  mcoor  = zeros(nofelt, nc);

  for n = 1 : nofelt
      mcoor(n, :) = mean(elem(n).coor);
  end

end

function ixbox = fillbox(cbox, coor)

# Input  cbox(ncoor, 2    Definition of box:
#                           xmin, xmax; ymin, ymax; zmin, zmax
#        coor(:, ncoor)   Coordinates
#
# Output ixbox(:)         Indices of points in box within coor
#
# -------------------------------------------------------------------
#
# x-coordinate

  ixbox = find(coor(:, 1) >= cbox(1, 1));
  ix    = find(coor(ixbox, 1) <= cbox(1, 2));
  ixbox = ixbox(ix);

# Remaining coordinates

  for k = 2 : rows(cbox)
     ix    = find(coor(ixbox, k) >= cbox(k, 1));
     ixbox = ixbox(ix);
     ix    = find(coor(ixbox, k) <= cbox(k, 2));
     ixbox = ixbox(ix);
  end

end
