function model = createAnalysisModel(K, M, mesh)
%CREATEANALYSISMODEL Copy assembled model data for side-effect-free solvers.
% The returned value is a plain struct. Numerical solvers receive this value
% instead of the mutable StructFEProblem facade or its FEMesh handle.

if ~ismatrix(K) || size(K, 1) ~= size(K, 2)
    error('MKEF:InvalidModel', 'The stiffness matrix must be square.');
end
if ~isequal(size(M), size(K))
    error('MKEF:InvalidModel', ...
        'The mass and stiffness matrices must have the same dimensions.');
end

numberOfDOFs = mesh.numberOfDOFs;
if size(K, 1) ~= numberOfDOFs
    error('MKEF:InvalidModel', ...
        'Matrix dimensions do not match the mesh degrees of freedom.');
end

model = struct();
model.stiffness = K;
model.mass = M;
model.dofMap = mesh.iMnod;
model.dofRegistry = mesh.dofRegistry;
model.releases = mesh.releases;
model.warnings = mesh.warnings;
model.fixedBoundaryConditions = normalizeSupports(mesh.allFixBCs, mesh.dofPerNode);
model.forceBoundaryConditions = mesh.allForceBCs;
model.elementLoads = mesh.elementLoads;
model.multiPointConstraints = mesh.multiPointConstraints;
model.nodeCoordinates = mesh.allNodes(:, 1:2);
model.elementData = mesh.allMeshElems;
model.sections = mesh.sections;
model.numberOfNodes = mesh.numberOfNodes;
model.dofPerNode = mesh.dofPerNode;
model.numberOfDOFs = numberOfDOFs;
end