Function to gather mesh information from gmesh files

Description

function [O1,O2,O3,O4,O5,O6,O7,O8,O9,O10]=getMesh(input)
%Function to gather mesh information from gmesh files
%---------------------------
% DEFINE SHORT-NAME OUTPUTS:
%---------------------------
% O1=msh; Mesh file from gmesh
% O2=Mdpts; Midpoint Database
% O3=cons; Various important constants
% O4=solidpts; Nodes (vertices and midpoints) on solid domain
% O5=fluidpts; Nodes (vertices and midpoints) on fluid domain
% O6=solperimeter; Nodes (vertices and midpoints) on solid/fluid bdry
% O7=outerbdrynodes; Vertices on the outer boundary
% O8=ELsolid; Global element numbers of vertices on solid domain
% O9=ELfluid; Global element numbers of vertices on fluid domain
% O10=edgetri; List of global element #'s in solid touching bdry
%--------------------------------------------------------------------------
global searchtriangles
global solidbdrynodes
disp('Gathering mesh information...')
%Choose particular mesh to use:
if input==1
usemesh='solidfluid.msh';
elseif input==2

usemesh='solidfluid2.msh';
elseif input==3
usemesh='solidfluid3.msh';
elseif input==4
usemesh='solidfluid4.msh';
elseif input==5
usemesh='solidfluid5.msh';
elseif input==6
usemesh='solidfluid6.msh';
elseif input==7
usemesh='solidfluid7.msh';
elseif input==-1 % <----debugging with mesh by hand
usemesh='fixedhandmesh1.msh';
elseif input==-2 % <----debugging with mesh by hand
usemesh='fixedhandmesh2.msh';
else
usemesh='solidfluid.msh';
end
disp([' - Using Mesh: ',usemesh])
msh = load_gmsh2(usemesh,[1 2]); % [1,2] requests variables for lines and
triangles only
M=msh.nbNod; % Number of nodes in mesh
NODECO=(1/3)*msh.POS(1:M,1:2); % NODE COordinates

msh.POS(1:M,1:2)=NODECO;
L=msh.nbTriangles; % Number of elements (triangular)
ELNODE=msh.TRIANGLES(1:L,1:4); % Global node numbers for all elements; element domain
chtype=ELNODE(1,end); % get tag for domain type; assume first type listed is fluid
k=0;
j=0;
for i=1:L % get a list of elements in each domain
if ELNODE(i,4)~=chtype
k=k+1;
ELsolid(k,1:3)=ELNODE(i,1:3); %coordinates of elements in solid domain
else
j=j+1;
ELfluid(j,1:3)=ELNODE(i,1:3); %coordinates of elements in fluid domain
end
end
nELsolid=k; % number of elements in the solid domain
nELfluid=j; % number of elements in the fluid domain
% The following information can be used to subdivide the boundary into segments
%global L1
%L1=msh.nbLines; % Number of element edges on the natural boundary

%global BDRYEDGES
BDRYEDGES=msh.LINES; % EDGES (node, node, ...)
[nBdryNodes,dum]=size(BDRYEDGES); % total number of nodes on a boundary
k=1;
k2=1;
chtype=BDRYEDGES(1,end);
for i=1:nBdryNodes
if BDRYEDGES(i,end)==chtype
outerbdrynodes(k,1)=BDRYEDGES(i,1);
k=k+1;
else
solidbdrynodes(k2,1)=BDRYEDGES(i,1);
k2=k2+1;
end
end
%----------------------
%Determine fluid/solid nodes
%----------------------
nfluidpts=0;
nsolidpts=0;
temp=ELfluid(:,1:3);
temp2=ELsolid(:,1:3);
for k=1:M
check=isempty(find(k==temp));

if check==0
nfluidpts=nfluidpts+1;
fluidpts(nfluidpts)=k;
end
check=isempty(find(k==temp2));
if check==0
nsolidpts=nsolidpts+1;
solidpts(nsolidpts)=k;
end
end
otherDim=nfluidpts;
nsolidnodes=nsolidpts;
%-------------------
% Midpoint Database
%-------------------
function [out]=mdptcalc(v1,v2)
% Gives midpoint (x,y) coordinates between global nodes v1 and v2
% also places a flag for the domain
% out(3)=0 means pt in interior of solid
% out(3)=1 means pt on solid bdry
% out(3)=2 means pt outside solid domain not on solid bdry
v1coord=NODECO(v1,1:2);
v2coord=NODECO(v2,1:2);
out(1)=.5*(v1coord(1)+v2coord(1));
out(2)=.5*(v1coord(2)+v2coord(2));
out(3)=0;

check1=isempty(find(v1==fluidpts)); %is v1 in fluid domain, 0=yes,1=no
check2=isempty(find(v2==fluidpts)); %is v2 in fluid domain, 0=yes,1=no
if check1+check2==0 %if both vertices are in fluid domain
index1=find(v1==solidbdrynodes); % where is v1 on solid bdry if at
all?
index2=find(v2==solidbdrynodes);
check3=isempty(index1);
check4=isempty(index2);
if check3+check4==0 %if both vertices are on the solid bdry
if abs(index1-index2)==1 || abs(index1-index2)==length(
solidbdrynodes)-1
out(3)=1;
end
else
out(3)=2;
end
end
end
disp('Generating midpoint database...')
disp(' ')
%Mdpts= [x-coord, y-coord, bdry?, local node #1, global elem #1, local node
#2,
% global element #2, fluid/solid]
% loctype is from 1-6 depending on what node is across
% midpoint could be used twice, so have 2 loctype and node2

% If a point is in both domains, notated as 1, fluid=2, solid=0
tol=10^(-8);
Mdpts=zeros(3*M,8); %pre-allocate space since no more than triple the nodes
for mdpts
Mdpts(:,3)=ones(3*M,1); %assume midpoint is on boundary until proven not
nMdpts=1; % will count the number of unique midpoints in the mesh
for i=1:L % cycle through each triangle (can we skip solid domain midpts?)
v1=ELNODE(i,1); v2=ELNODE(i,2); v3=ELNODE(i,3);
mid=zeros(3,4);
mid(1,1:3)=mdptcalc(v2,v3); %organize midpoints in same way as local basis
phi's
mid(2,1:3)=mdptcalc(v3,v1);
mid(3,1:3)=mdptcalc(v1,v2);
mid(:,4)=[4;5;6];
%Check to see if any of the above are repeats
for j=1:3
vec=ones(nMdpts,1);
diff_x=abs(mid(j,1)*vec-Mdpts(1:nMdpts,1));
diff_y=abs(mid(j,2)*vec-Mdpts(1:nMdpts,2));
check=min(diff_x+diff_y);
if check>tol
Mdpts(nMdpts,1:2)=mid(j,1:2);
Mdpts(nMdpts,4)=mid(j,4);
Mdpts(nMdpts,5)=i;
Mdpts(nMdpts,8)=mid(j,3);
nMdpts=nMdpts+1;

else
[dum,place]=min(diff_x+diff_y);
Mdpts(place,6:7)=[j+3,i];
Mdpts(place,3)=0; %mark that your midpoint isn't on bdry
end
end
end
nMdpts=nMdpts-1;
Mdpts=Mdpts(1:nMdpts,1:8);
%----------------------------------------
%Determine the fluid/solid midpoints
%----------------------------------------
temp=Mdpts(:,8);
count=0;
clear solperimeter
for k=1:nMdpts
if temp(k)~=0
nfluidpts=nfluidpts+1;
fluidpts(nfluidpts)=k+M;
end
if temp(k)~=2
nsolidpts=nsolidpts+1;
solidpts(nsolidpts)=k+M;
if temp(k)==1
count=count+1;

solperimeter(count,1)=k+M;
end
end
end
solperimeter=[solidbdrynodes; solperimeter];
nsolperimeter=length(solperimeter);
%Build a list of global node numbers that intersect the solid boundary,
%edgetri
edgetri=zeros(0,2);
count=0;
chtype=ELNODE(1,end); % get tag for domain type; assume first type
listed is fluid
for q=1:L
if ELNODE(q,end)~=chtype % if we are looking at an element in solid
domain
count=count+1;
v1=ELNODE(q,1);
test1=isempty(find(v1==solidbdrynodes)); %see if any of the vertices
are on solid bdry
v2=ELNODE(q,2);
test2=isempty(find(v2==solidbdrynodes)); %test=0 means yes on the bdry
v3=ELNODE(q,3);
test3=isempty(find(v3==solidbdrynodes));
if test1+test2+test3==1 % if two are on bdry



edgetri=[edgetri;[q,count]]; %record the global element number and
solid global element number
end
end
end
nedgetri=length(edgetri); % number of solid elements on bdry
cons=[nELsolid,nELfluid,otherDim,nsolidpts,nfluidpts,nMdpts,nsolperimeter,
nedgetri];
searchtriangles=[Mdpts(:,5),Mdpts(:,7)]; %Global element numbers for each
midpoint
% DEFINE SHORT-NAME OUTPUTS:
O1=msh;
O2=Mdpts;
O3=cons;
O4=solidpts;
O5=fluidpts;
O6=solperimeter;
O7=outerbdrynodes;
O8=ELsolid;
O9=ELfluid;
O10=edgetri;
https://matlab1.com/jupyter-using-it-and-its-great-functions/
http://gmsh.info/doc/texinfo/gmsh.html