Writing classes with parallel properties (R2016b)

I am trying to write a class to hold finite element model data. I have decided to write a 'Grid' class to represent the grid points in the FEA mesh, with 'GridID' and 'Position' as its properties. I have also made it a handle-type class, as I want to do cross-referencing with a 'TriangularElement' class. This means that if I change the ID of one of my grid points, the corresponding grid ID in all elements will be updated automatically.

My initial approach was to have one Grid object for each grid point in my model, and I've quickly realized that this does not scale very well to the models that I have. It took 15 mins to create an object array for 100k grid points, and when saving, the .mat file came up to around 10Mb. This is an issue for me as I have models with 1m+ grid points.

The most sensible solution I've found so far is from Andrew in this thread Is MATLAB OOP slow or am I doing something wrong?, where he talks about writing classes such that an instance wraps an array of objects. He also brings up plane organization, as seen in Ways to Organize Data in Structure Arrays. The key is making the properties of my class 'parallel' in his terms. I was not able to come up with a straightforward solution to implement this however. How should I define a requirement to ensure that my GridID and Position arrays always have, say, the same number of rows?

[Guillaume]

It certainly isn't particularly fast, but it takes about 21s on my computer (R2019a) to create a 1 million elements object array and the mat files is only around 38 kB. Nowhere near 15 minutes.

You're still better off having a scalar class though. I'm a bit unclear on where you're stuck though. Having some example class code may help clarify.

For what is worth, this is the class I tested with (not production worthy code at all):

classdef gridelem < handle
    properties (SetAccess = private)
        ID(:, 1) uint64;
        Location(:, 3) double;
    end
    methods
        function this = gridelem(locations)
            if nargin~= 0    %nargin == 0 happens when the array is created on the next line.
                this(size(locations, 1), 1) = gridelem;  %allocate object array, before filling up location
                locations  = num2cell(locations, 2);
                [this.Location] = locations{:};
            end
            IDs = num2cell(1:numel(this));  %fill up IDs even if no location was provided.
            [this.ID] = IDs{:};
        end
    end
end

%testing class
g = gridelem(rand(1e6, 3))