function mesh_3tri(varargin)
% mesh_3tri computes a mesh with 3 nodes per element
%
global numlps
global meltyp
global ndiv
global join
global rn
global xcor
global ycor
global xlc ylc nplc
global lnr
global lnp
global inp
global xord 
global yord
global npn
global numele 
global numnp
global numel
global np
global numsnp
global edgearray

if nargin>0
    infile=varargin{1};
    outfile=varargin{2};
else
    infile = [];
    outfile=[];
end

read_in_data(infile);
%calc_numnp;
 
numele=0;
inp=0;          

for i=1:numlps 
	
	i1=ndiv(i,1);   
	i4=ndiv(i,2); 
	
	switch meltyp
	case 3
	i2=i1+1;
	i5=i4+1;
	i7=i2*i5;
    case 4
	i2=i1+1;
	i5=i4+1;
	i7=i2*i5;
    case 6
	i2=2*i1+1;
	i5=2*i4+1;
	i7=i2*i5;
    case 9
    i2=2*i1+1;
	i5=2*i4+1;
	i7=i2*i5;
	end
	

	side_array_calc(i,i1,i2,i5,i7); 

	npn_array_calc(i,i7);
 
	lnp_array_calc(i,i2,i5);

	coord_calc(i,i1,i2,i4,i5);

	np_array_calc(i,i1,i2,i4);

end

numnp=inp;

edge_array_calc;

[idiff,jmem,memjt] = setup(numnp,numel,np);
[ib,jnt] = optnum(numnp,numel,idiff,jmem,memjt);
update_nodes(jnt);


out_data(outfile,ib,jnt);




