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.

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.

[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]

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 to importdata(). With that correction both now return the same data, just retained the code for comparison of methods.

I also augmented the Mat_Coord_Elemt array to deal with the 3-node elements. The array must retain the 13 column width already defined so I kept the same structure with the divider "0" in column 5; just NaN fill elements for the fourth node and associated last two coordinate values. Added an if ... else if ... end clause to deal with the undefined indices.

The refined code then looks like

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
    % 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);                             % ensure in bounds of array
  Mat_Coord_Elemt = zeros(num_elements, 13);                  % Preallocate matrix
  barycentre = zeros(num_elements, 2);                        % preallocate for COG, x, y
  for i = 1:num_elements
    nodes=Mat_top(i,1:4);
    % make up somethng to handle triangular nodes
    nnodes=nnz(isfinite(nodes));                              % count how many are defined (4, 3 only handled now)
    if nnodes==4
      Mat_Coord_Elemt(i,:)=[nodes 0 reshape(coord_nodes(nodes,2:3).',1,[])];
      barycentre(i,:)=mean(coord_nodes(Mat_Coord_Elemt(i,1:nnodes),2:3));
    elseif nnodes==3
      Mat_Coord_Elemt(i,:)=[nodes 0 reshape(coord_nodes(nodes(1:3),2:3).',1,[]) nan nan]; % fill missing locaions w/NaN
      barycentre(i,:)=mean(coord_nodes(Mat_Coord_Elemt(i,1:nnodes),2:3));
    else
      msg=sprintf('Bad nodes, element %d: %d %d %d %d',i,nodes);
      warning(msg)
    end
 end
 classement=[1:num_elements].';   % Assign a unique number to the center of gravity -- sequential numbering scheme
 clear Mat_top s                                             % Clear some memory
% Display results
disp('Barycentre Coordinates:');
%disp(barycentre);
disp('Node Numbering:');
%disp(classement);
  
% Create a figure
hF=figure; hAx=axes(hF);
hold(hAx,'on')

% Plot each element
for i = 1:num_elements
  x=Mat_Coord_Elemt(i,6:2:end);
  y=Mat_Coord_Elemt(i,7:2:end);
  fill(x,y, 'b', 'FaceAlpha', 0.5);
  % Calculate the center of gravity for labeling
  center_x = mean(x);
  center_y = mean(y);
  % 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 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');                        % converts row
    try
      top_data(i,1:numel(vals))=vals;
    catch
      disp(top_text(i))
      vals
    end
  end
end
end

This plots the figure of the mesh; the attempt to put the element numbers in and the axis equal turns the figure into just a black streak so I commented those out leaving the attached figure image; play with as wish.

Not sure about what the three node elements really are; it appears they are all just three points on a line so they don't have any area; you'll have to look at what they mean as elements in the actual code and modify as/if necessry; this just handles the indexing to keep from crashing with undefined indices in the data arrays.

It turns out in the end while there are 2425 elements there are actually some 7000+ nodes/points so once got rid of the undefined indices, that problem did go away other than figuring out what the real interpretation of those is in the actual model.

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