Grid generation fails: matlab interp1

I am creating Navier-Stokes simulation of flow around particle aggregates imitating marine snow. I need the simulations to be programmatically produced to look at stochastic effects. When creating the grid, I get the error:

"Error using matlab.internal.math.interp1
Sample points must be unique and sorted in ascending order.

Error in featool"

Interestingly, the code sometimes works but other times not - I assume it has something to do with the stochastic production of sub-particles.

Below is the code that I use to generate the geometry, and how I wish to create the mesh. I am not sure if the error can be solved by "sorting" some aspect of the fea.geom.object properties or if the error is located in the gridgen function. Any help would be appreciated.

Best
Ben




clear all
close all
clc

% Geometry and grid parameters.
h         = 0.02;         % Height of rectangular domain, m.
l         = 0.02;         % Length of rectangular domain, m.
xc        = h/2;          % x-coordinate of particle center.
yc        = l/2;          % y-coordinate of particle center.
diam      = 0.003;        % Particle, m.

% Grid generation.
fea.sdim = { 'x' 'y' };

% Rectangular background object
gobj1 = gobj_rectangle( 0, h, 0, l, 'R1' );
fea = geom_add_gobj( fea, gobj1 );

% Particles
npart=10; % Number of sub particles

for kk=1:npart
    particleName=[ 'CE' num2str(kk)];
    orient=rand*2*pi;
    dist=rand*diam/2;
   
    partX=cos(orient)*dist+xc;
    partY=sin(orient)*dist+yc;
    particleRadius=rand*200E-6+50E-6;
   
    gobj = gobj_circle( [partX partY], particleRadius,particleName);
    fea = geom_add_gobj( fea, gobj );
end

% Subtract particles from background geometry
str = strjoin(fea.geom.tags,' - ');
[fea.geom, stat] = geom_apply_formula( fea.geom, str,1E-12);

% Grid generation
fea.grid = gridgen( fea, 'gridgen', 'gridgen2d', 'hmax', 10e-5, 'nsm', 25, 'q', 0.65 );

Dear Ben,

Thank you for the issue report. As the geometry is randomly generated I have not been able to produce a configuration that gives this error. Could you run your script again, and save (and attach) the resulting data and workspace produced after the error with a command such as:

  save my_workspace.mat

That is very interesting indeed. I attached the mat file with a FEA struct that does not work on my machine.

Brief info: Running Matlab 2020b with FEATool 1.13.0

( Hope this way of uploading works)

 FEA_Output_NoMesh.mat

Thank you, yes with this data and I believe I have located the issue (precision issue due to too small cut geometry edge). I will push an update soon:ish.

In the meantime perhaps you can try to skip and re-generate configurations that fail, you could use try-catch blocks to avoid triggering errors, for example something like:

  try
     fea.grid = gridgen( ...
  catch
     disp('grid generation failed, retrying...')
  end

Perfect thanks a lot!

Just to update, this issue should have been corrected with the new build, available from the GitHub repository (v1.13.1 Build 21.02.048):

https://github.com/precise-simulation/featool-multiphysics/releases/latest

That is awesome thanks a lot! This bug does not appear any more.

The geometry generation now works great, however there seems to be an issue with the I/O interface to the Fortran Gridgen2D framework. I get the error:

startio: illegal unit number
apparent state: unit -413691250

Seems like there are negative numbers in the data being sent to Gridgen?

Attached is a geometry that causes this error.
FEATool_FortranIOBug.mat

Interesting, it looks like a memory overflow type error of some sort. However, I could not reproduce this particular error on my system with a gridgen call such as

  fea.grid = gridgen( fea, 'gridgen', 'gridgen2d', 'hmax', 50e-3, 'nsm', 25, 'q', 0.65 );

but any lower 'hmax' value seems to take forever. Please share the exact gridgen function call signature for the attached mat file if you still observe this issue. You can also try using the 'triangle' and 'gmsh' grid/mesh generators to see if any of them work better in these cases.

bborer wrote

The geometry generation now works great, however there seems to be an issue with the I/O interface to the Fortran Gridgen2D framework. I get the error:

startio: illegal unit number
apparent state: unit -413691250

Seems like there are negative numbers in the data being sent to Gridgen?

Attached is a geometry that causes this error.
FEATool_FortranIOBug.mat

Interesting - my hmax is not that small however I try to get a very good resolution close to the main boundary of the particles. The gridgen call is:

        % Grid generation
        n_bdr = length(fea.geom.objects{1,1}.boundaries);   % Number of boundaries.
        gridsize=100E-6;
        bndSize=20E-6;
        bndVec=ones(1,n_bdr)*bndSize;
        bndVec(1:4)=gridsize;
        fea.grid = gridgen( fea, 'gridgen', 'gridgen2d', 'hmax', gridsize ,'hmaxb',bndVec, 'nsm', 25, 'q', 0.65 );


Thus an average grid of 100 microns and boundary grid of 20 microns (which is important due to the small features inside the domain).

I will check how I could make the overall grid requirement smaller and/or change the gridgen algorithm to see if that solves the issue.

Thank you for including the grid generation command, I still don't get any crash but the grid generator never seems too finish. It could potentially be a system/os dependent error in this case (if memory runs out or so).

What I do notice that the initial boundary mesh distribution already contains a million or so grid points which will lead to an extremely dense mesh. Are you sure you have the geometry dimensions in consistent units? As the geometry in the FEATool_FortranIOBug.mat file is of dimension ~5, and with a target boundary mesh size is 20-100e-6 would eventually lead to a mesh with the order of (5/50e-6)^2 = 1e10 grid cells which would be hard to fit in to memory (and take very long to solve).

If you really do need this mesh size and the other mesh generators (triangle or gmesh) doesn't work another work around might be to first create a coarser grid (but refined in regions of interest), and then use the "gridrefine" command to make a finer grid. For example

  n_bdr = length(fea.geom.objects{1}.boundaries);
  gridsize = 1;
  bndSize = 0.05;
  bndVec = ones(1,n_bdr)*bndSize;
  bndVec(1:4) = gridsize;
  fea.grid = gridgen( fea, 'gridgen', 'gridgen2d', 'hmax', gridsize ,'hmaxb', bndVec, 'q', 0.9 );
  
  for i=1:2
    fea = gridrefine(fea);
    fea = gridsmooth(fea);
  end
  
  plotgrid(fea)