How to import mesh data *.PF3 To Matlab

Hi,

I imported mesh data from a commercial finite element method software into MATLAB. I am now attempting to construct the topological matrix.

Is it possible for you to help me program these instructions in MATLAB? Please and thank you.

[david]

mesh data

[dpb]

"... for you to help me program ..."

The key word there is "help". Folks here will be glad to help when you get stuck, but the forum isn't a free consulting service. What have you done so far and where, specifically, did you get stymied?

These are very explicit instructions that appear to have been prepared for some existing package -- from whence did you obtain the pasted information; that would seem to be the lodestone. If you're trying to replicate a commercial product to avoid purchase, that's a whole different matter.

[dpb]

fn="data mesh.zip";
unzip(fn)
dir

.              ..             data mesh      data mesh.zip  

fn1="data mesh";  fn2="datamesh";
movefile(fn1,fn2)
d=dir(fn2);
d=d(~[d.isdir].');
cd(d(1).folder)
!ls -ls

total 860
848 -rw-rw-rw- 1 matlab matlabgroup 865088 Nov 12 09:28 Distrib.PF3
 12 -rw-rw-rw- 1 matlab matlabgroup   9561 Nov 25 08:25 untitled4.m

type 'untitled4'

clc
clear all

% Define the file name
filename = 'Distrib.PF3'; % Replace with your actual file name

% Open the file and read its content
fid = fopen(filename, 'r');
s = fread(fid, '*char')'; % Read the entire file as a character array
fclose(fid);

% Find the index for "DESCRIPTEUR DE TOPOLOGIE DES ELEMENTS"
index_topology = strfind(s, 'DESCRIPTEUR DE TOPOLOGIE DES ELEMENTS');
topology_data = importdata(filename, ' ', index_topology);

% Extract node numbers from topology.data (even lines only)
Mat_top = topology_data.data(2:2:end, :); % Get even lines

% Find the index for "COORDONNEES DES NOEUDS"
index_coords = strfind(s, 'COORDONNEES DES NOEUDS');

coord_data = importdata(filename, ' ', index_coords);

% Extract node coordinates
class(coord_data)
%coord_nodes = coord_data.data;
%if isempty(idx_coord)
        % Créer des coordonnées fictives pour la démonstration
       % num_nodes = max(Mat_top(:));
      %  coord_nodes = [(1:num_nodes)', rand(num_nodes,1)*100, ...
                     %  rand(num_nodes,1)*100, zeros(num_nodes,1)];
    %else
coord_section = s(index_coords:min(index_coords+5000, length(s)));
        coord_lines = strsplit(coord_section, '\n');
        coord_nodes = [];
        for i = 1:length(coord_lines)
            nums = str2num(coord_lines{i});
            if length(nums) == 4
                coord_nodes = [coord_nodes; nums];
            end
        end
  

% Define a cell array
%coord_data.data = {'one', 'two', 'three'; 100, 200, rand(3,3)};
% Accessing all elements in the first row
%all_elements = [coord_data.data{:}]; % Concatenate all elements into a single array
%disp(all_elements); % Displays all elements
% Accessing all elements in the cell array
%disp(class(coord_data)); % Check the type of coord_data

% Shift coordinates by X and Y
X = 0; % Example shift in x direction
Y = 0;  % Example shift in y direction
o=size(coord_nodes);

coord_nodes(:, 2) = coord_nodes(:, 2) + X; % Shift x coordinates
coord_nodes(:, 3) = coord_nodes(:, 3) + Y; % Shift y coordinates

% Create Mat_Coord_Elemt matrix
num_elements = size(Mat_top, 1);
Mat_Coord_Elemt = zeros(num_elements, 12); % Preallocate matrix

for i = 1:Mat_Coord_Elemt
    node1 = Mat_top(i, 1);
    node2 = Mat_top(i, 2);
    node3 = Mat_top(i, 3);
    node4 = Mat_top(i, 4);

    Mat_Coord_Elemt(i, :) = [
        node1, node2, node3, node4, 0, ...
        coord_nodes(node1, 2), coord_nodes(node1, 3), ...
        coord_nodes(node2, 2), coord_nodes(node2, 3), ...
        coord_nodes(node3, 2), coord_nodes(node3, 3), ...
        coord_nodes(node4, 2), coord_nodes(node4, 3)
    ];
end

% Calculate center of gravity for each block
barycentre = zeros(num_elements, 2);
for i = 1:num_elements
    barycentre(i, :) = [
        mean([coord_nodes(Mat_Coord_Elemt(i, 1), 2), ...
 coord_nodes(Mat_Coord_Elemt(i, 2), 2), ...
              coord_nodes(Mat_Coord_Elemt(i, 3), 2), ...
              coord_nodes(Mat_Coord_Elemt(i, 4), 2)]), ...
        mean([coord_nodes(Mat_Coord_Elemt(i, 1), 3), ...
              coord_nodes(Mat_Coord_Elemt(i, 2), 3), ...
              coord_nodes(Mat_Coord_Elemt(i, 3), 3), ...
              coord_nodes(Mat_Coord_Elemt(i, 4), 3)])
    ];
end

% Determine number of elements in x and y directions
x_coords = barycentre(:, 1);
y_coords = barycentre(:, 2);
xref = min(x_coords);
yref = min(y_coords);

Nbr_Elmt_x_Rt = sum(x_coords == xref);
Nbr_Elmt_y_Rt = sum(y_coords == yref);

% Clear Mat_top to free memory
clear Mat_top;
% Assuming Mat_Coord_Elemt is already defined and contains the node data
num_elements = size(Mat_Coord_Elemt, 1);
classement = zeros(num_elements, 1); % Initialize classement vector
barycentre = zeros(num_elements, 2); % Initialize barycentre matrix

for i = 1:num_elements
    % Extract node coordinates for the current element
    nodes = Mat_Coord_Elemt(i, 1:4); % Assuming first four columns are node indices
    x_coords = Mat_Coord_Elemt(nodes, end-1); % Get x coordinates
    y_coords = Mat_Coord_Elemt(nodes, end );  % Get y coordinates

    % Calculate center of gravity
    barycentre(i, 1) = mean(x_coords); % Average x
    barycentre(i, 2) = mean(y_coords); % Average y

    % Store the center gravity in Mat_Coord_Elemt
    Mat_Coord_Elemt(i, end-1) = barycentre(i, 1);
    Mat_Coord_Elemt(i, end) = barycentre(i, 2);

    % Assign a unique number to the center of gravity
    classement(i) = i; % Simple numbering scheme
end

% Display results
disp('Barycentre Coordinates:');
disp(barycentre);
disp('Node Numbering:');
disp(classement);
% Example data for Mat_Coord_Elemt and classement
Mat_Coord_Elemt = [1, 2, 3, 4;  % Example node indices for elements
                    5, 6, 7, 8]; % More elements
classement = [1; 2]; % Node numbering

% Coordinates for visualization (example)
x_coords = [0, 1, 1, 0;  % Rectangle 1
            1, 2, 2, 1]; % Rectangle 2
y_coords = [0, 0, 1, 1;  % Rectangle 1
            0, 0, 1, 1]; % Rectangle 2

% Create a figure
figure;
hold on;

% Plot each element
for i = 1:size(x_coords, 1)
    fill(x_coords(i, :), y_coords(i, :), 'b', 'FaceAlpha', 0.5); % Draw rectangles
    % Calculate the center of gravity for labeling
    center_x = mean(x_coords(i, :));
    center_y = mean(y_coords(i, :));
    % Annotate with the node number
    text(center_x, center_y, num2str(classement(i)), 'HorizontalAlignment', 'center', 'Color', 'k');
end

% Set axis properties
axis equal;
xlabel('X Coordinate');
ylabel('Y Coordinate');
title('Visualization of Reluctance Network');
grid on;
hold off;
function [barycentre_Reg, Regions, coord_nodes_Reg] = calculateMeshCenters(meshFiles, rep, delimiterIn, headerlinesIn)
    % Initialize output variables
    barycentre_Reg = [];
    Regions = {};
    coord_nodes_Reg = [];

    % Loop through each mesh file
    for i = 1:length(meshFiles)
        % Construct the full file path
        filePath = fullfile(rep, meshFiles{i});

        % Read the mesh data, skipping header lines
        data = readtable(filePath, 'Delimiter', delimiterIn, 'HeaderLines', headerlinesIn);

        % Extract region names and coordinates
        regionNames = data.RegionName; % Assuming 'RegionName' is a column in your data
        vertices = data{:, {'X', 'Y', 'Z'}}; % Assuming 'X', 'Y', 'Z' are coordinate columns

        % Calculate the center of gravity for the region
        barycentre = mean(vertices, 1); % Average of the vertices

        % Store results
        barycentre_Reg = [barycentre_Reg; barycentre]; % Append barycenter
        Regions = [Regions; regionNames]; % Append region names
        coord_nodes_Reg = [coord_nodes_Reg; vertices]; % Append coordinates
    end

    % Adjust coordinates for model shift (if needed)
    % Example: shift upward by a certain value (e.g., 1 unit)
    shiftValue = 1; % Define your shift value
    coord_nodes_Reg = coord_nodes_Reg + shiftValue; % Modify coordinates

    % Convert Regions to a cell array if necessary
    Regions = unique(Regions); % Ensure unique region names
end
function connectionMatrix = generateConnectionMatrix(numNodes)
    % Initialize the connection matrix
    connectionMatrix = zeros(numNodes, 4); % Each row corresponds to [lower, left, right, upper]

    for node = 1:numNodes
        % Determine the type of node (off-border, left border, right border, etc.)
        if isOffBorderNode(node, numNodes)
            connectionMatrix(node, :) = getOffBorderConnections(node);
        elseif isBorderNode(node, numNodes)
            connectionMatrix(node, :) = getBorderConnections(node);
        else
            connectionMatrix(node, :) = getInterpolationNodeConnections(node);
        end
    end
end

function isOffBorder = isOffBorderNode(node, numNodes)
    % Logic to determine if the node is an off-border node
    isOffBorder = (node > 1 && node < numNodes); % Example condition
end

function isBorder = isBorderNode(node, numNodes)
    % Logic to determine if the node is a border node
    isBorder = (node == 1 || node == numNodes); % Example condition
end

function connections = getOffBorderConnections(node)
    % Define connections for an off-border node
    connections = [lowerNode(node), leftNode(node), rightNode(node), upperNode(node)];
end

function connections = getBorderConnections(node)
    % Define connections for a border node
    connections = [lowerNode(node), rightNode(node), leftNode(node), upperNode(node)];
end

function connections = getInterpolationNodeConnections(node)
    % Special case for interpolation nodes
    connections = [centralNode(node), 0, 0, 0];
end

function ln = lowerNode(node)
    % Logic to get the lower node connection
    ln = node - 1; % Example logic
end

function ln = leftNode(node)
    % Logic to get the left node connection
    ln = node - 1; % Example logic
end

function ln = rightNode(node)
    % Logic to get the right node connection
    ln = node + 1; % Example logic
end

function ln = upperNode(node)
    % Logic to get the upper node connection
    ln = node + 1; % Example logic
end

function cn = centralNode(node)
    % Logic to get the central node for interpolation
    cn = node; % Example logic
end
disp(class(coord_data));
disp(fieldnames(coord_data));
coord_data.data = someData; % where someData is your actual data
coord_nodes = coord_data{1}; % if coord_data is a cell array

f=readlines('Distrib.PF3');
f(1:20)

ans = 20×1 string array
    " Fichier cree par F3DXPE 1.0. Date: 11/12/25 10:28:07                            "
    "       3           NOMBRE DE DIMENSIONS DU DECOUPAGE                           "
    "    2425           NOMBRE  D'ELEMENTS                                          "
    "       0           NOMBRE  D'ELEMENTS VOLUMIQUES                               "
    "    2275           NOMBRE  D'ELEMENTS SURFACIQUES                              "
    "     150           NOMBRE  D'ELEMENTS LINEIQUES                                "
    "       0           NOMBRE  D'ELEMENTS PONCTUELS                                "
    "       0           NOMBRE DE MACRO-ELEMENTS                                    "
    "    7026           NOMBRE DE POINTS                                            "
    "      13           NOMBRE DE REGIONS                                           "
    "       0           NOMBRE DE REGIONS VOLUMIQUES                                "
    "       9           NOMBRE DE REGIONS SURFACIQUES                               "
    "       4           NOMBRE DE REGIONS LINEIQUES                                 "
    "       0           NOMBRE DE REGIONS PONCTUELLES                               "
    "       0           NOMBRE DE REGIONS MACRO-ELEMENTAIRES                        "
    "      20           NOMBRE DE NOEUDS DANS 1 ELEMENT (MAX)                       "
    "      20           NOMBRE DE POINTS D'INTEGRATION / ELEMENT (MAX)              "
    " NOMS DES REGIONS                                            "
    "     REGIONS SURFACIQUES                                         "
    " GAP                                                                             "

% well, let's try it...
untitled4

ans = 'cell'

Error using :
Colon operands must be real scalars.

Error in untitled4 (line 64)
for i = 1:Mat_Coord_Elemt
^^^^^^^^^^^^^^^^^^^^^^^^^

Very difficult to debug in this environment, but it appears you have the basic code already -- set a breakpoint and see what looks to be going wrong -- it appears from the above error message that the Mat_top array isn't correct since that's where the size elements are coming from that error.

Have you modified this code from what it was from wherever it came?

[dpb]

This code is very inefficient given that it reads the file at least three times very early on if not more times later than have looked through quickly. importdata is also somewhat flaky in what and how it returns results. It might be easier debugging to rewrite it from its basic structure and instead of reading the file in initially as one long char() string variable, s as is now done, read the file as a string array with readlines and then one could find the specific lines and parse them directly.

Just a thought, but you should be able to see where/why the present fails -- one would guess that importdata didn't return what was expected at some point...either the file doesn't quite follow the expected pattern or perhaps there's a behavior difference in importdata since the version the code was written/debugged with.

importdata doesn't refer to any specific changes in its functionality, but some of the helper functions may have that aren't reflected. While not (at least yet) disparaged, Mathworks has certainly shifted focus away from using it to the newer toolsets like readcell and friends so it wouldn't be terribly surprising if there were some subtle changes in behavior.

Answer 1

[dpb]

OK, I had a little time once discovered there was code already there and not just starting from the description.

I recast initially into a function; using scripts is very dangerous as well as the use of the clear words wreaks havoc on existing workspace if doing anything at all else and so should be avoided when possible.

The reworked function through the nodes now works -- obviously renaming it to something more apropos to its purpose would be a starting detail <vbg>.

The first point where your code breaks is where the length of the string to be parsed was artifically limited to be 5000 characters; in the following line...

  coord_section = s(index_coords:min(index_coords+5000, length(s)));

Instead, I added

  ixEnd=index_coords+strfind(s(index_coords:end),' ==== DECOUPAGE  TERMINE')-1; % find end of section -- probably better way but works
  coord_section = s(index_coords:ixEnd);                                        % only get to where want

that limits the search to the section itself and backs up to not pickup the trailing text.

Then, having the raw data of the file at hand, parse the numerical data directly...MATLAB works in column order (down looking a a 2D array), so to read the data as is written we must transpose the lines to be columns and then convert. Once converted, then transpose the output array back to having four columns.

  % let's work on the raw char data directly instead...
  ixStart=strfind(coord_section,newline)+1;                                     % skip the header line
  coord_nodes=reshape(sscanf(coord_section(ixStart:end).','%f'),4,[]).';        % convert to double, reshape

After doing the above, there's no need for converting the raw data into cellstr arrays and using the very slow str2num to parse.

That doesn't answer the problem above trying to run the initial code; that comes later in the following section --

  % Create Mat_Coord_Elemt matrix
  num_elements = size(Mat_top, 1);
  Mat_Coord_Elemt = zeros(num_elements, 12); % Preallocate matrix
  for i = 1:Mat_Coord_Elemt

That is clearly wrong; Mat_Coord_Elemt is the array; it's height() is num_elements from above so that loop range should be

   for i = 1:num_elements           % use height, not array name to iterate into
      node1 = Mat_top(i, 1);
      node2 = Mat_top(i, 2);
      node3 = Mat_top(i, 3);
      node4 = Mat_top(i, 4);
  
      Mat_Coord_Elemt(i, :) = [
          node1, node2, node3, node4, 0, ...
          coord_nodes(node1, 2), coord_nodes(node1, 3), ...
          coord_nodes(node2, 2), coord_nodes(node2, 3), ...
          coord_nodes(node3, 2), coord_nodes(node3, 3), ...
          coord_nodes(node4, 2), coord_nodes(node4, 3)
      ];

The above then still fails on the assignment line, however, with

Unable to perform assignment because the size of the left side is 1-by-12 and 
                                     the size of the right side is 1-by-13.
Error in readSomeUnknownPF3 (line 78)
      Mat_Coord_Elemt(i, :) = [ 

which shows that there are 13 elements in the RH array but only room for 12. Reading the specification, the "0" is a separator of some function so the issue is the array isn't defined large enough to hold the needed 13 elements.

Mat_Coord_Elemt = zeros(num_elements, 13); % Preallocate matrix leave room for separator "0"

But, with the presumption it is the extraneous 0 that is the problem, the above can also be written more succinctly as

for i = 1:num_elements           % use height, not array name to iterate into
  nodes=Mat_top(i,1:4);
  Mat_Coord_Elemt(i,:)=[nodes 0 reshape(coord_nodes(nodes,2:3).',1,[])];

I made those changes in the attached file; you can revert back to the prior if the presumption is incorrect and figure out what should be if not.

I didn't try any further, but that should get you started.

As noted before, if this was supposed to have been working code, is there a real original copy somewhere? The "untitled4" makes one wonder...

Good luck..

function res=readSomeUnknownPF3(filename)
  % Define the file name
  if nargin==0
    filename = 'Distrib.PF3'; % Replace or call with your actual file name
  end
  % Open the file and read its content
  fid = fopen(filename, 'r');
  s = fread(fid, '*char')'; % Read the entire file as a character array
  fclose(fid);
  
  % Find the index for "DESCRIPTEUR DE TOPOLOGIE DES ELEMENTS"
  index_topology = strfind(s, 'DESCRIPTEUR DE TOPOLOGIE DES ELEMENTS');
  topology_data = importdata(filename, ' ', index_topology);
  
  % Extract node numbers from topology.data (even lines only)
  Mat_top = topology_data.data(2:2:end, :); % Get even lines
  
  % Find the index for "COORDONNEES DES NOEUDS"
  index_coords = strfind(s, 'COORDONNEES DES NOEUDS');
  
  coord_data = importdata(filename, ' ', index_coords);
  
  % Extract node coordinates
  %class(coord_data)
  %coord_nodes = coord_data.data;
  %if isempty(idx_coord)
          % Créer des coordonnées fictives pour la démonstration
         % num_nodes = max(Mat_top(:));
        %  coord_nodes = [(1:num_nodes)', rand(num_nodes,1)*100, ...
                       %  rand(num_nodes,1)*100, zeros(num_nodes,1)];
      %else
  %coord_section = s(index_coords:min(index_coords+5000, length(s)));           % cuts off before end
  ixEnd=index_coords+strfind(s(index_coords:end),' ==== DECOUPAGE  TERMINE')-1; % find end of section -- probably better way but works
  coord_section = s(index_coords:ixEnd);                                        % only get to where want
  % let's work on the raw char data directly instead...
  ixStart=strfind(coord_section,newline)+1;                                     % skip the header line
  coord_nodes=reshape(sscanf(coord_section(ixStart:end).','%f'),4,[]).';        % convert to double, reshape
  % then don't need to do all this converting back and forth...
  % more efficient yet would be to find the location in the file and read in directly, but...
  %coord_lines = strip(split(coord_section,newline));                            % and strip blanks
  %coord_lines=coord_lines(~cellfun(@isempty,coord_lines));                      % remove empty lines
  %coord_nodes = [];
  %for i = 1:length(coord_lines)
  %    nums = str2num(coord_lines{i});
  %    if length(nums) == 4
  %        coord_nodes = [coord_nodes; nums];
  %    end
  %end
    
  
  % Define a cell array
  %coord_data.data = {'one', 'two', 'three'; 100, 200, rand(3,3)};
  % Accessing all elements in the first row
  %all_elements = [coord_data.data{:}]; % Concatenate all elements into a single array
  %disp(all_elements); % Displays all elements
  % Accessing all elements in the cell array
  %disp(class(coord_data)); % Check the type of coord_data
  
  % Shift coordinates by X and Y
  X = 0; % Example shift in x direction
  Y = 0;  % Example shift in y direction
  o=size(coord_nodes);
  
  coord_nodes(:, 2) = coord_nodes(:, 2) + X; % Shift x coordinates
  coord_nodes(:, 3) = coord_nodes(:, 3) + Y; % Shift y coordinates
  
  % Create Mat_Coord_Elemt matrix
  num_elements = size(Mat_top, 1);
  Mat_Coord_Elemt = zeros(num_elements, 13); % Preallocate matrix leave room for separator "0"
  for i = 1:num_elements        % height(), not array itself Mat_Coord_Elemt
      % node1 = Mat_top(i, 1);
      % node2 = Mat_top(i, 2);
      % node3 = Mat_top(i, 3);
      % node4 = Mat_top(i, 4);
      % 
      % Mat_Coord_Elemt(i, :) = [
      %     node1, node2, node3, node4, 0, ...
      %     coord_nodes(node1, 2), coord_nodes(node1, 3), ...
      %     coord_nodes(node2, 2), coord_nodes(node2, 3), ...
      %     coord_nodes(node3, 2), coord_nodes(node3, 3), ...
      %     coord_nodes(node4, 2), coord_nodes(node4, 3)
      % ];
    nodes=Mat_top(i,1:4);
    Mat_Coord_Elemt(i,:)=[nodes reshape(coord_nodes(nodes,2:3).',1,[])];
  end
  
  % Calculate center of gravity for each block
  barycentre = zeros(num_elements, 2);
  for i = 1:num_elements
      barycentre(i, :) = [
          mean([coord_nodes(Mat_Coord_Elemt(i, 1), 2), ...
   coord_nodes(Mat_Coord_Elemt(i, 2), 2), ...
                coord_nodes(Mat_Coord_Elemt(i, 3), 2), ...
                coord_nodes(Mat_Coord_Elemt(i, 4), 2)]), ...
          mean([coord_nodes(Mat_Coord_Elemt(i, 1), 3), ...
                coord_nodes(Mat_Coord_Elemt(i, 2), 3), ...
                coord_nodes(Mat_Coord_Elemt(i, 3), 3), ...
                coord_nodes(Mat_Coord_Elemt(i, 4), 3)])
      ];
  end
  
  % Determine number of elements in x and y directions
  x_coords = barycentre(:, 1);
  y_coords = barycentre(:, 2);
  xref = min(x_coords);
  yref = min(y_coords);
  
  Nbr_Elmt_x_Rt = sum(x_coords == xref);
  Nbr_Elmt_y_Rt = sum(y_coords == yref);
  
  % Clear Mat_top to free memory
  clear Mat_top;
  % Assuming Mat_Coord_Elemt is already defined and contains the node data
  num_elements = size(Mat_Coord_Elemt, 1);
  classement = zeros(num_elements, 1); % Initialize classement vector
  barycentre = zeros(num_elements, 2); % Initialize barycentre matrix
  
  for i = 1:num_elements
      % Extract node coordinates for the current element
      nodes = Mat_Coord_Elemt(i, 1:4); % Assuming first four columns are node indices
      x_coords = Mat_Coord_Elemt(nodes, end-1); % Get x coordinates
      y_coords = Mat_Coord_Elemt(nodes, end );  % Get y coordinates
  
      % Calculate center of gravity
      barycentre(i, 1) = mean(x_coords); % Average x
      barycentre(i, 2) = mean(y_coords); % Average y
  
      % Store the center gravity in Mat_Coord_Elemt
      Mat_Coord_Elemt(i, end-1) = barycentre(i, 1);
      Mat_Coord_Elemt(i, end) = barycentre(i, 2);
  
      % Assign a unique number to the center of gravity
      classement(i) = i; % Simple numbering scheme
  end
  
  % Display results
  disp('Barycentre Coordinates:');
  disp(barycentre);
  disp('Node Numbering:');
  disp(classement);
  % Example data for Mat_Coord_Elemt and classement
  Mat_Coord_Elemt = [1, 2, 3, 4;  % Example node indices for elements
                      5, 6, 7, 8]; % More elements
  classement = [1; 2]; % Node numbering
  
  % Coordinates for visualization (example)
  x_coords = [0, 1, 1, 0;  % Rectangle 1
              1, 2, 2, 1]; % Rectangle 2
  y_coords = [0, 0, 1, 1;  % Rectangle 1
              0, 0, 1, 1]; % Rectangle 2
  
  % Create a figure
  figure;
  hold on;
  
  % Plot each element
  for i = 1:size(x_coords, 1)
      fill(x_coords(i, :), y_coords(i, :), 'b', 'FaceAlpha', 0.5); % Draw rectangles
      % Calculate the center of gravity for labeling
      center_x = mean(x_coords(i, :));
      center_y = mean(y_coords(i, :));
      % Annotate with the node number
      text(center_x, center_y, num2str(classement(i)), 'HorizontalAlignment', 'center', 'Color', 'k');
  end
  
  % Set axis properties
  axis equal;
  xlabel('X Coordinate');
  ylabel('Y Coordinate');
  title('Visualization of Reluctance Network');
  grid on;
  hold off;
  
    function [barycentre_Reg, Regions, coord_nodes_Reg] = calculateMeshCenters(meshFiles, rep, delimiterIn, headerlinesIn)
      % Initialize output variables
      barycentre_Reg = [];
      Regions = {};
      coord_nodes_Reg = [];
  
      % Loop through each mesh file
      for i = 1:length(meshFiles)
          % Construct the full file path
          filePath = fullfile(rep, meshFiles{i});
  
          % Read the mesh data, skipping header lines
          data = readtable(filePath, 'Delimiter', delimiterIn, 'HeaderLines', headerlinesIn);
  
          % Extract region names and coordinates
          regionNames = data.RegionName; % Assuming 'RegionName' is a column in your data
          vertices = data{:, {'X', 'Y', 'Z'}}; % Assuming 'X', 'Y', 'Z' are coordinate columns
  
          % Calculate the center of gravity for the region
          barycentre = mean(vertices, 1); % Average of the vertices
  
          % Store results
          barycentre_Reg = [barycentre_Reg; barycentre]; % Append barycenter
          Regions = [Regions; regionNames]; % Append region names
          coord_nodes_Reg = [coord_nodes_Reg; vertices]; % Append coordinates
      end
  
      % Adjust coordinates for model shift (if needed)
      % Example: shift upward by a certain value (e.g., 1 unit)
      shiftValue = 1; % Define your shift value
      coord_nodes_Reg = coord_nodes_Reg + shiftValue; % Modify coordinates
  
      % Convert Regions to a cell array if necessary
      Regions = unique(Regions); % Ensure unique region names
    end
    function connectionMatrix = generateConnectionMatrix(numNodes)
      % Initialize the connection matrix
      connectionMatrix = zeros(numNodes, 4); % Each row corresponds to [lower, left, right, upper]
  
      for node = 1:numNodes
          % Determine the type of node (off-border, left border, right border, etc.)
          if isOffBorderNode(node, numNodes)
              connectionMatrix(node, :) = getOffBorderConnections(node);
          elseif isBorderNode(node, numNodes)
              connectionMatrix(node, :) = getBorderConnections(node);
          else
              connectionMatrix(node, :) = getInterpolationNodeConnections(node);
          end
      end
    end
  
    function isOffBorder = isOffBorderNode(node, numNodes)
      % Logic to determine if the node is an off-border node
      isOffBorder = (node > 1 && node < numNodes); % Example condition
    end
  
    function isBorder = isBorderNode(node, numNodes)
      % Logic to determine if the node is a border node
      isBorder = (node == 1 || node == numNodes); % Example condition
    end
  
    function connections = getOffBorderConnections(node)
      % Define connections for an off-border node
      connections = [lowerNode(node), leftNode(node), rightNode(node), upperNode(node)];
    end
  
    function connections = getBorderConnections(node)
      % Define connections for a border node
      connections = [lowerNode(node), rightNode(node), leftNode(node), upperNode(node)];
    end
  
    function connections = getInterpolationNodeConnections(node)
      % Special case for interpolation nodes
      connections = [centralNode(node), 0, 0, 0];
    end
  
    function ln = lowerNode(node)
      % Logic to get the lower node connection
      ln = node - 1; % Example logic
    end
  
    function ln = leftNode(node)
      % Logic to get the left node connection
      ln = node - 1; % Example logic
    end
  
    function ln = rightNode(node)
      % Logic to get the right node connection
      ln = node + 1; % Example logic
    end
  
    function ln = upperNode(node)
      % Logic to get the upper node connection
      ln = node + 1; % Example logic
    end
  
    function cn = centralNode(node)
      % Logic to get the central node for interpolation
      cn = node; % Example logic
    end
  disp(class(coord_data));
  disp(fieldnames(coord_data));
  coord_data.data = someData; % where someData is your actual data
  coord_nodes = coord_data{1}; % if coord_data is a cell array
end

[dpb]

ADDENDUM:

There's a problem that the assignment of

node1 = Mat_top(i, 1);
node2 = Mat_top(i, 2);
node3 = Mat_top(i, 3);
node4 = Mat_top(i, 4);

Mat_Coord_Elemt(i, :) = [
   node1, node2, node3, node4, 0,                ...
   coord_nodes(node1, 2), coord_nodes(node1, 3), ...
   coord_nodes(node2, 2), coord_nodes(node2, 3), ...
   coord_nodes(node3, 2), coord_nodes(node3, 3), ...
   coord_nodes(node4, 2), coord_nodes(node4, 3)];

fails for quite a number owing to one or more of the nodes being 0 in Mat_top(i,1:4) and so are invalid array indices.

I don't know if this indicates a problem in Mat_top itself or whether the zero elements should be special-cased; you'll have to use the base code specification document to figure that out.

ADDENDUM

It is, in fact a problem in reading the data file -- as I suspected using the result of strfind as an index into importdata as a number of lines to skip is wrong -- strfind() is counting characters, not lines.

In this case, the result for the numeric topology data is

K>> topology_data.data(1:20,:)
ans =
        1399         575         576        1400        4265        4214        4266        4259
        1443           4         303           4           3           0           8           8
           0           0           0           0         NaN         NaN         NaN         NaN
        1398         574         575        1399        4264        4191        4265        4257
        1444           4         303           4           3           0           8           8
           0           0           0           0         NaN         NaN         NaN         NaN
        1397         573         574        1398        4263        4168        4264        4255
        1445           4         303           4           3           0           8           8
           0           0           0           0         NaN         NaN         NaN         NaN
         113           3         573        1397        2568        4145        4263        4253
        1446           4         303           4           3           0           8           8
           0           0           0           0         NaN         NaN         NaN         NaN
        1396        1400         455         456        4260        4261        4262        4250
        1447           4         303           4           3           0           8           8
           0           0           0           0         NaN         NaN         NaN         NaN
        1395        1399        1400        1396        4258        4259        4260        4247
        1448           4         303           4           3           0           8           8
           0           0           0           0         NaN         NaN         NaN         NaN
        1394        1398        1399        1395        4256        4257        4258        4244
        1449           4         303           4           3           0           8           8
K>>

whereas the actual data file content at that point is

K>> s(index_topology:index_topology+20*80)
ans =
    'DESCRIPTEUR DE TOPOLOGIE DES ELEMENTS                       
                   1        4      303        1        3        0        8        8        0        0        0        0
                2376      761       79      805     7026     6175     3452     7016
                   2        4      303        1        3        0        8        8        0        0        0        0
                2375      762      761     2376     7025     6177     7026     7014
                   3        4      303        1        3        0        8        8        0        0        0        0
                2374      763      762     2375     7024     6179     7025     7012
                   4        4      303        1        3        0        8        8        0        0        0        0
                2373      764      763     2374     7023     6181     7024     7010
                   5        4      303        1        3        0        8        8        0        0        0        0
                2372      765      764     2373     7022     6183     7023     7008
                   6        4      303        1        3        0        8        8        0        0        0        0
                2371      766      765     2372     7021     6185     7022     7006
                   7        4      303        1        3        0        8        8        0        0        0        0
                2370      767      766     2371     7020     6187     7021     7004
                   8        4      303        1        3        0        8        8        0        0        0        0
                2369      768      767     2370     7019'
K>> 

One needs to count lines, not characters. As noted before, this could be solved by reading the whole file as either a string array by line or a cell array by individual cells or by parsing the file content and reading directly.

Error using feof
Invalid file identifier. Use fopen to generate a valid file identifier.

Answer 2

[dpb]

Other comments getting rather long so will add a new finding as another Answer...the above comments/findings are all still true; this is one more needed patch-up. I attached my local copy with this update as well...

Try replacing the topology reading section

  % Find the index for "DESCRIPTEUR DE TOPOLOGIE DES ELEMENTS"
  %index_topology = strfind(s, 'DESCRIPTEUR DE TOPOLOGIE DES ELEMENTS');
  %topology_data = importdata(filename, ' ', index_topology);
  
  % Extract node numbers from topology.data (even lines only)
  %Mat_top = topology_data.data(2:2:end, :); % Get even lines

with

Mat_top=get_top(s);   % do it more better

where get_top() is

function top_data=get_top(s)
  s=split(s,newline);                                     % turns into cellstr array by line
  ix=find(contains(s,"NOMBRE  D'ELEMENTS"));              % in case number header lines can change
  Nelements=sscanf(s{ix(1)},'%d');                        % the number elements; avoid other element types
  ix=find(contains(s,'DESCRIPTEUR DE TOPOLOGIE DES ELEMENTS'))  ; % we found the heading
  i1=ix+2; i2=i1+Nelements;
  top_data=s(i1:2:i2);
  % now read the topology data
  top_data=reshape(sscanf(char(top_data).','%f'),8,[]).'; % convert by column, then reshape back
end

This now runs past the prior failure point of having zero-valued nodes down to

>> readSomeUnknownPF3
Index in position 1 exceeds array bounds. Index must not exceed 1213.
Error in readSomeUnknownPF3 (line 129)
      x_coords = Mat_Coord_Elemt(nodes, end-1); % Get x coordinates 
>> 

I don't have any more time to try to dig further but I suspect the issue is probably that the Mat_Coord_Elemt array isn't correct.owing to similar issues in indexing into the file as was done by character position instead of by line number.

If you open the output file in a text editor and look, you'll probably be able to spot the problem in what is being interpreted as those values vis a vis what should be. Again guessing, probably replacing that code with something more on the lines of what I show above will be the cure.

NOTA BENE: in the above code the use of the input problem size data of the number of elements to compute where in the file the end of the desired section is. This, undoubtedly, could be used effectively for the other data sections if one will look at the model and the output and observe how they are connected and make use of that information.

[dpb]

ERRATUM:

I forgot to account for there being two lines per element in computing the second index range above, it should be

i1=ix+2; i2=i1+2*(Nelements-1);

which returns

K>> [i1 i2]
ans =
49        4897
K>> 

which agrees with looking at the output file in editor.

However, this then fails because at line 4599 the number of nodes in the odd-numbered lines switches from 8 to only 3 per record so the end number of nodes is no longer divisible by eight.

I don't know what this means nor how to work around it; only knowing more about the actual model and digging into the description will be able to fix up this, methinks.

The original code runs to get a regular array, but we've already seen it fails later with zero node numbers.

[dpb]

The original code using importdata NaN-fills the short rows; if this is the intended behavior, then use

function top_data=get_top(s)
  s=split(s,newline);                                     % turns into cellstr array by line
  ix=find(contains(s,"NOMBRE  D'ELEMENTS"));              % in case number header lines can change
  Nelements=sscanf(s{ix(1)},'%d');                        % the number elements; avoid other element types
  ix=find(contains(s,'DESCRIPTEUR DE TOPOLOGIE DES ELEMENTS'))  ; % we found the heading
  i1=ix+2; i2=i1+2*(Nelements-1);
  top_text=s(i1:2:i2);
  % following can't be used if thre are non-qauadrilateral elements with fewer than eight values
  %top_data=reshape(sscanf(char(top_text).','%f'),8,[]).'; % convert by column, then reshape back
  top_data=nan(height(top_text),8);                       % preallocate nan-filled for short rows
  for i=1:height(top_text)
    vals=sscanf(top_text{i},'%f');
    top_data(i,1:numel(vals))=vals;
  end
end

With this, then back to the prior result that the code doesn't know how to handle 0 or NaN node numbers.

Answer 3

[dpb]

OK, one more iteration -- I cleaned up the input section significantly and added the ability to choose between the original idea of using importdata() mulitple times versus parsing the in-memory data directly. I corrected the erroneous logic of the original that used the character position in the file instead of the line number of the target record as the headerlines argument. Both return the same data, just retained for comparison.

function res=readSomeUnknownPF3(filename,USEIMPORTDATA,X,Y)
% read unidentified FEM input file
% USAGE:
%   readSomeUnknownPF3(file, useimport, Xoffset, YOffset)
%       file      - input file, default Distr.PF3
%       useimport - reread file multiple times with importdata(), default False
%       X, YOffset- lateral shift of X, Y coordinates; default 0

% modified from source from unknown author -- dpb 11Nov25
% original reread data from file with importdata() but used character offsets in file instead of line numbers
% fixed that to do either, but default is to parse the file internally since already been read into memory

  arguments                                                     % validate, provide default arguments
    filename (1,:) char = 'Distrib.PF3'
    USEIMPORTDATA (1,1) logical = false
    X (1,1) double = 0
    Y (1,1) double = 0
  end
  % Open the file and read its content
  fid = fopen(filename, 'r');
  s = fread(fid, '*char')';                                   % Read the entire file as a character array
  if USEIMPORTDATA
    % Find the ROW indices for topology and nodes sections, NOT character position for importdata() use
    frewind(fid)
    index_topology=0;                                         % topology header line number
    index_coords=0;                                           % coordinates header line number
    index=1;                                                  % file line counter
    while index_topology==0 || index_coords==0                % until find both
      l=fgetl(fid);
      if contains(l, 'DESCRIPTEUR DE TOPOLOGIE DES ELEMENTS'), index_topology=index; end
      if contains(l, 'COORDONNEES DES NOEUDS'),                index_coords=index; end
      index=index+1;
    end
    fclose(fid);                                              % close file either way
    topology_data = importdata(filename, ' ', index_topology);
    % Extract node numbers from topology.data (even lines only)
    Mat_top = topology_data.data(2:2:end, :);                 % Get even lines
    % now get coordinates data
    coord_data = importdata(filename, ' ', index_coords);     % reread the file from header line
    coord_nodes = coord_data.data;
  else                                                        % parse the char string containing full file read
    fclose(fid);                                              % close file either way
    Mat_top=get_top(s);                                       % parse the string -- must by row since not all 8 elements
    % Later code doesn't handle triangular elements, remove the Mat_top elements with NaN members
    %Mat_top_all=Mat_top;                                     % hang  onto all for later maybe...
    %Mat_top=Mat_top(all(isfinite(Mat_top),2),:);             % keep only those without NaN
    % now the coordinates from internal file
    index_coords = strfind(s, 'COORDONNEES DES NOEUDS');      % in terms of char() string, not file rows
    ixEnd=index_coords+strfind(s(index_coords:end),' ==== DECOUPAGE  TERMINE')-1; % find end of section -- probably better way but works
    coord_section = s(index_coords:ixEnd);                                        % only get to where want
    % let's work on the raw char data directly instead...
    ixStart=strfind(coord_section,newline)+1;                                     % skip the header line
    coord_nodes=reshape(sscanf(coord_section(ixStart(1):end).','%f'),4,[]).';     % convert to double, reshape
  end
  
  % Shift coordinates by X and Y
  coord_nodes(:, 2) = coord_nodes(:, 2) + X;                  % Shift x coordinates
  coord_nodes(:, 3) = coord_nodes(:, 3) + Y;                  % Shift y coordinates
  
  % Create Mat_Coord_Elemt matrix
  num_elements = height(Mat_top);
  Mat_Coord_Elemt = zeros(num_elements, 13);                  % Preallocate matrix
  for i = 1:num_elements
    nodes=Mat_top(i,1:4);
    try
      Mat_Coord_Elemt(i,:)=[nodes 0 reshape(coord_nodes(nodes,2:3).',1,[])];
    catch
      msg=sprintf('Bad nodes, element %d: %d %d %d %d',i,nodes);
      warning(msg)
    end
  end
  % Calculate center of gravity for each block
  barycentre = zeros(num_elements, 2);
  for i = 1:num_elements
    try
      barycentre(i, :) = mean(coord_nodes(Mat_Coord_Elemt(i,1:4),2:3));
    catch
      msg=sprintf('Bad nodes COG, element %d: ',i);
      warning(msg)
    end
  end
  ...

The above is through the calculation of the COG positions; you will notice I also simplified that to use MATLAB vector syntax instead of the initial enumeration of all the elements.

This still has the issue that the original code and the documentation have no provisions to handle anything but quadrilateral elements with the 8 nodes/element; the elements with only three nodes at the end of the model thus create invalid addressing errors with only three points. Deciding how to and implementing code to handle those is, to quote a famous phrase, left as "exercise for Student". <g>

I also again attached the whole code as is locally...