function mnt = plate_mount(P, dom)
% PLATE_MOUNT  How the plate is held, as data on the grid.
%
%   mnt = plate_mount(P, dom), with P from plate_setup and dom from
%   plate_domain, turns params.mount into what the solvers need: the
%   direction of gravity in the plate's coordinates, the nodes the mount holds
%   still, the springs, dashpots and one-sided contacts it puts under the
%   plate, the in-plane constraints for the self-weight stress solve
%   (plate_prestress), and the rigid-body basis the results are decomposed in
%   (plate_rigid_basis).  Every mount is something that can be built on a
%   bench, and each carries an idealization that is stated in mnt.notes.
%
%   Gravity.  The plate lies in the (x, y) plane and w is along its normal z.
%   mnt.gx, mnt.gy, mnt.gz is the unit vector along which gravity pulls, in
%   these coordinates, and P.g its magnitude.  The in-plane part (gx, gy) puts
%   the plate under its own weight in tension or compression (plate_prestress);
%   the normal part gz loads it transversely (plate_static).
%
%   Mounts.
%
%   'thread'    One thread through a hole at (x, y), of length L, from a fixed
%               hook.  A plate hung from one point turns in its own plane until
%               its centre of mass is straight below the hole, which fixes
%               (gx, gy).  For small out-of-plane motion the thread is a
%               transverse spring of stiffness T / L at the hole, T = m g its
%               tension (the linearized pendulum), and a torsional spring
%               G J / L about the vertical axis with J = pi d^4 / 32 from the
%               thread's diameter d and shear modulus G (mount.thread).
%   'threads'   Two threads of the same length L from hooks at the same
%               height, through holes 1 and 2.  The line of the holes is then
%               horizontal, and the tensions come from statics: the plate is
%               refused when its centre of mass is not between the holes,
%               because one thread would go slack.
%   'pads'      A horizontal plate resting on pads (foam, felt or a knife
%               edge): each pad is a rectangle of half-widths (rx, ry) about
%               (x, y), with a stiffness k [N/m^3] and damping c [N s/m^3] per
%               unit area.  A pad pushes only while compressed, so the plate
%               can lift off and land again.
%   'cymbal'    A horizontal plate on a stand: a rod of radius r_rod through
%               the centre (no transverse effect) between two felt washers of
%               outer radius r_felt, each of stiffness k and damping c per
%               unit area and both compressed by the wing nut's preload
%               [m].  Each felt is a one-sided contact.
%   'vise'      One edge gripped in a vise to a depth [m]: the nodes within
%               the jaws are held at w = 0 and the clamped condition follows
%               from the bending energy (plate_bending).  The angle tilt [deg]
%               sets the orientation: 0 hangs straight down from the jaw, 90 is
%               a horizontal cantilever, 180 stands up on the jaw with the
%               self-weight in compression.
%   'pinched'   Held between finger and thumb over a patch of radius r about
%               (x, y): a two-sided spring k [N/m] and dashpot c [N s/m] spread
%               over the patch, the patch also fixed in the plane.  tilt [deg]
%               runs from 0 (hanging below the hand) to 90 (held flat).
%
%   Fields.
%     type                      the mount
%     gx, gy, gz                unit direction of gravity in plate coordinates
%     held                      ny x nx logical: nodes held at w = 0 (vise)
%     ks, cs                    ny x nx: linear transverse spring [N/m] and
%                               dashpot [N s/m] at each node (threads, pinch)
%     kb, cb, ob                bottom contacts: stiffness [N/m], dashpot
%                               [N s/m] and surface offset [m] per node; the
%                               contact pushes up while ob - w > 0
%     kt, ct, ot                top contacts: push down while ot + w > 0
%     ktor                      torsional stiffness about the vertical axis
%                               through the pivot [N m/rad] (one thread)
%     fix                       in-plane constraints for plate_prestress: rows
%                               [node, dir_x, dir_y] meaning the displacement
%                               of node along (dir_x, dir_y) is zero
%     free_rotation             the in-plane rotation about the first fixed
%                               node is free (one thread) and must be
%                               regularized, which needs the load to have no
%                               moment about it
%     pivot                     [x y] the rigid-body basis is taken about
%     ed, eh                    in-plane unit vectors: down and horizontal
%     nr                        number of rigid-body coordinates (0 or 3)
%     rigid_names               {1 x nr} names of the coordinates
%     supports                  struct array: name, unit, nodes (index list),
%                               kind ('spring', 'bottom' or 'top': which force
%                               the series reports at those nodes)
%     mass, xc, yc              plate mass [kg] and centre of mass, from the
%                               node weights
%     notes                     remarks about the idealization

mount = P.mount;
nx = dom.nx; ny = dom.ny; dx = dom.dx; dy = dom.dy;
X = dom.X; Y = dom.Y;
inside = dom.inside;
a = dom.m * dx * dy;                       % node areas [m^2]
Mi = P.mu * a;                             % node masses [kg]
mass = sum(Mi(:));
xc = sum(Mi(:) .* X(:)) / mass;
yc = sum(Mi(:) .* Y(:)) / mass;
zero = zeros(ny, nx);
mnt.type = mount.type;
mnt.held = false(ny, nx);
mnt.ks = zero; mnt.cs = zero;
mnt.kb = zero; mnt.cb = zero; mnt.ob = zero;
mnt.kt = zero; mnt.ct = zero; mnt.ot = zero;
mnt.ktor = 0;
mnt.fix = zeros(0, 3);
mnt.free_rotation = false;
mnt.mass = mass; mnt.xc = xc; mnt.yc = yc;
supports = struct('name', {}, 'unit', {}, 'nodes', {}, 'kind', {});
notes = {};

switch mount.type
    case 'thread'
        kA = nearest_node([mount.x mount.y], dom);
        A = [X(kA) Y(kA)];                                 % the plate hangs from the node the hole snaps to
        d = [xc yc] - A;
        if norm(d) < 1e-3 * max(P.plate.Lx, P.plate.Ly)
            error('plate_mount: the hole is at the centre of mass; a plate hung there has no preferred orientation');
        end
        ed = d / norm(d);
        eh = [-ed(2) ed(1)];
        mnt.gx = ed(1); mnt.gy = ed(2); mnt.gz = 0;
        T = mass * P.g;
        mnt.ks(kA) = T / mount.length;
        mnt.fix = [kA 1 0; kA 0 1];
        mnt.free_rotation = true;
        mnt.ktor = torsion(mount);
        mnt.pivot = [X(kA) Y(kA)];
        mnt.ed = ed; mnt.eh = eh;
        mnt.nr = 3;
        mnt.rigid_names = {'sway of the hole', 'swing about the hole', 'spin about the thread'};
        supports(end + 1) = struct('name', 'thread', 'unit', 'N', 'nodes', kA, 'kind', 'spring');
        ang = atan2(ed(1), -ed(2)) * 180 / pi;
        notes{end + 1} = sprintf('hung from the node (%.4g, %.4g) m; the plate turns %.3g degrees in its plane so the centre of mass hangs below the hole; tension %.3g N gives a transverse spring of %.3g N/m at the hole', X(kA), Y(kA), ang, T, T / mount.length);
        if mnt.ktor > 0
            notes{end + 1} = sprintf('spin about the thread is restrained only by its torsion, %.3g N m/rad', mnt.ktor);
        else
            notes{end + 1} = 'spin about the thread is free: nothing but the sigma0 loss slows it';
        end
    case 'threads'
        A1 = [mount.x1 mount.y1]; A2 = [mount.x2 mount.y2];
        k1 = nearest_node(A1, dom); k2 = nearest_node(A2, dom);
        if k1 == k2
            error('plate_mount: the two holes fall on the same grid node; move them apart or refine the grid');
        end
        A1 = [X(k1) Y(k1)]; A2 = [X(k2) Y(k2)];
        e = (A2 - A1) / norm(A2 - A1);
        ed = [-e(2) e(1)];                                 % perpendicular to the line of the holes
        if dot([xc yc] - A1, ed) < 0
            ed = -ed;                                      % the plate hangs below its threads
        end
        eh = [-ed(2) ed(1)];
        s2 = norm(A2 - A1);
        sc = dot([xc yc] - A1, eh) / dot(A2 - A1, eh);      % centre of mass along the line, 0 at hole 1, 1 at hole 2
        if ~(sc > 0 && sc < 1)
            error('plate_mount: the centre of mass is not between the two holes, so one thread would go slack');
        end
        T1 = mass * P.g * (1 - sc);
        T2 = mass * P.g * sc;
        mnt.gx = ed(1); mnt.gy = ed(2); mnt.gz = 0;
        mnt.ks(k1) = T1 / mount.length;
        mnt.ks(k2) = T2 / mount.length;
        mnt.fix = [k1 ed; k2 ed; k1 eh];                   % threads inextensible, sway fixed at one hole
        mnt.pivot = (A1 + A2) / 2;
        mnt.ed = ed; mnt.eh = eh;
        mnt.nr = 3;
        mnt.rigid_names = {'sway of the threads', 'swing about the holes', 'yaw about the vertical'};
        supports(end + 1) = struct('name', 'thread 1', 'unit', 'N', 'nodes', k1, 'kind', 'spring');
        supports(end + 1) = struct('name', 'thread 2', 'unit', 'N', 'nodes', k2, 'kind', 'spring');
        notes{end + 1} = sprintf('hung from the nodes (%.4g, %.4g) and (%.4g, %.4g) m on threads of %.3g m; tensions %.3g and %.3g N; the line of the holes is horizontal', A1(1), A1(2), A2(1), A2(2), mount.length, T1, T2);
        notes{end + 1} = sprintf('span %.3g m between the holes', s2);
    case 'pads'
        mnt.gx = 0; mnt.gy = 0; mnt.gz = -1;
        if ~isfield(mount, 'pads') || isempty(mount.pads)
            error('plate_mount: the pads mount needs at least one pad');
        end
        npads = numel(mount.pads);
        for k = 1:npads
            pd = mount.pads(k);
            in = inside & abs(X - pd.x) <= pd.rx + 1e-12 & abs(Y - pd.y) <= pd.ry + 1e-12;
            if ~any(in(:))
                kn = nearest_node([pd.x pd.y], dom);
                in = false(ny, nx); in(kn) = true;
                notes{end + 1} = sprintf('pad %d is smaller than a grid cell, so it holds the one node nearest its centre', k);
            end
            mnt.kb(in) = mnt.kb(in) + mount.k * a(in);
            mnt.cb(in) = mnt.cb(in) + mount.c * a(in);
            supports(end + 1) = struct('name', sprintf('pad %d', k), 'unit', 'N', 'nodes', find(in), 'kind', 'bottom');
        end
        mnt.pivot = [xc yc];
        mnt.ed = [1 0]; mnt.eh = [0 1];
        mnt.nr = 3;
        mnt.rigid_names = {'heave', 'tilt along x', 'tilt along y'};
        notes{end + 1} = sprintf('resting on %d pad(s); each pushes only while compressed, so the plate can lift off', npads);
    case 'cymbal'
        mnt.gx = 0; mnt.gy = 0; mnt.gz = -1;
        r = sqrt((X - P.plate.Lx / 2).^2 + (Y - P.plate.Ly / 2).^2);
        in = inside & r <= mount.r_felt + 1e-12 & r > mount.r_rod;
        if ~any(in(:))
            error('plate_mount: the felt washer covers no node; enlarge r_felt or refine the grid');
        end
        mnt.kb(in) = mount.k * a(in); mnt.cb(in) = mount.c * a(in); mnt.ob(in) = mount.preload;
        mnt.kt(in) = mount.k * a(in); mnt.ct(in) = mount.c * a(in); mnt.ot(in) = mount.preload;
        supports(end + 1) = struct('name', 'lower felt', 'unit', 'N', 'nodes', find(in), 'kind', 'bottom');
        supports(end + 1) = struct('name', 'upper felt', 'unit', 'N', 'nodes', find(in), 'kind', 'top');
        mnt.pivot = [P.plate.Lx / 2, P.plate.Ly / 2];
        mnt.ed = [1 0]; mnt.eh = [0 1];
        mnt.nr = 3;
        mnt.rigid_names = {'heave', 'tilt along x', 'tilt along y'};
        notes{end + 1} = sprintf('on a stand: felts over %d node(s) between radii %.3g and %.3g m, each preloaded by %.3g mm; the rod and the hole it passes through are not modelled', nnz(in), mount.r_rod, mount.r_felt, mount.preload * 1e3);
    case 'vise'
        depth = mount.depth;
        switch mount.edge
            case 'x0'
                coord = X; edir = [1 0];
            case 'x1'
                coord = P.plate.Lx - X; edir = [-1 0];
            case 'y0'
                coord = Y; edir = [0 1];
            case 'y1'
                coord = P.plate.Ly - Y; edir = [0 -1];
            otherwise
                error('plate_mount: unknown vise edge "%s" (expected x0, x1, y0 or y1)', mount.edge);
        end
        held = inside & coord <= depth + 1e-12;
        % the clamped condition needs a held node behind every held edge node
        h = min(dx, dy);
        if nnz(held) < 2 || max(coord(held)) < h - 1e-12
            held = inside & coord <= h + 1e-12;
            notes{end + 1} = sprintf('jaw depth widened to %.3g mm, two rows of nodes, which the clamped condition needs', h * 1e3);
        end
        if ~any(held(:))
            error('plate_mount: the vise grips no node of the plate');
        end
        if nnz(inside & ~held) < 9
            error('plate_mount: the vise leaves almost nothing of the plate free');
        end
        mnt.held = held;
        tilt = mount.tilt * pi / 180;
        mnt.gx = cos(tilt) * edir(1); mnt.gy = cos(tilt) * edir(2); mnt.gz = -sin(tilt);
        kh = find(held);
        mnt.fix = [kh ones(numel(kh), 1) zeros(numel(kh), 1); kh zeros(numel(kh), 1) ones(numel(kh), 1)];
        mnt.pivot = [xc yc];
        mnt.ed = edir; mnt.eh = [-edir(2) edir(1)];
        mnt.nr = 0;
        mnt.rigid_names = {};
        notes{end + 1} = sprintf('%d nodes within %.3g mm of the %s edge are held in the vise; tilt %g degrees', nnz(held), depth * 1e3, mount.edge, mount.tilt);
        if mount.tilt > 90
            notes{end + 1} = 'standing on the jaw: the self-weight is compressive, which softens the lowest mode and can buckle a long thin plate';
        end
    case 'pinched'
        A = [mount.x mount.y];
        r = sqrt((X - A(1)).^2 + (Y - A(2)).^2);
        in = inside & r <= mount.radius + 1e-12;
        if ~any(in(:))
            kA = nearest_node(A, dom);
            in = false(ny, nx); in(kA) = true;
            notes{end + 1} = 'the pinch patch is smaller than a grid cell, so it holds one node';
        end
        atot = sum(a(in));
        mnt.ks(in) = mount.k * a(in) / atot;
        mnt.cs(in) = mount.c * a(in) / atot;
        kin = find(in);
        mnt.fix = [kin ones(numel(kin), 1) zeros(numel(kin), 1); kin zeros(numel(kin), 1) ones(numel(kin), 1)];
        d = [xc yc] - A;
        if norm(d) < 1e-9
            d = [0 -1];
        end
        ed = d / norm(d);
        tilt = mount.tilt * pi / 180;
        mnt.gx = cos(tilt) * ed(1); mnt.gy = cos(tilt) * ed(2); mnt.gz = -sin(tilt);
        mnt.pivot = A;
        mnt.ed = ed; mnt.eh = [-ed(2) ed(1)];
        mnt.nr = 3;
        mnt.rigid_names = {'heave at the hand', 'swing about the hand', 'turn about the hand'};
        supports(end + 1) = struct('name', 'fingers', 'unit', 'N', 'nodes', kin, 'kind', 'spring');
        notes{end + 1} = sprintf('pinched over %d node(s) about (%.4g, %.4g) m: a two-sided spring of %.3g N/m with %.3g N s/m of damping, tilt %g degrees', nnz(in), A(1), A(2), mount.k, mount.c, mount.tilt);
end
mnt.supports = supports;
mnt.notes = notes;
end

function k = nearest_node(A, dom)
d2 = (dom.X - A(1)).^2 + (dom.Y - A(2)).^2;
d2(~dom.inside) = Inf;
[dmin, k] = min(d2(:));
if ~isfinite(dmin)
    error('plate_mount: no node of the plate near (%.4g, %.4g) m', A(1), A(2));
end
end

function kt = torsion(mount)
% Torsional stiffness G J / L of a round thread; zero without thread data.
kt = 0;
if isfield(mount, 'thread') && isfield(mount.thread, 'diameter') && isfield(mount.thread, 'G')
    J = pi * mount.thread.diameter^4 / 32;
    kt = mount.thread.G * J / mount.length;
end
end