%----------------------
%   PROGRAM Poisson.m
%----------------------

      clear
%------------------------------
%     LOAD DATA
%------------------------------
      load MESHo   -ASCII
      load NODES   -ASCII
      load NP      -ASCII

      NUMNP = MESHo(1);
      NUMEL = MESHo(2);
      NNPE  = MESHo(3);
      if NNPE ~= 3
	error('NNPE in MESHo  must equal 3')
      end
      for I=1:NUMNP;
          XORD(I)  =NODES(I,1);
          YORD(I)  =NODES(I,2);
          NPcode(I)=NODES(I,3);
      end
%     ----------------------- 
%     Declare sparse matrices
%     ----------------------- 
      SK  = sparse(NUMNP,NUMNP);

%     -----------------------
%     General Initialization
%     -----------------------
      Q(1:NUMNP,1)  =0.0;
      PHI(1:NUMNP)  =0.0;
      NPBC(1:NUMNP) =0;

%     -----------------------
%     User's Initialization
%     -----------------------
      INITIAL

%     ---------------
%     Place Q in RHS
%     ---------------
      RHS=Q;

%     ------------------------------------
%     Formation of Finite Element Matrices
%     Element by Element
%     ------------------------------------
      for I=1:NUMEL;
         rc=NP(I,:);  % row and column numbers
         XJ=XORD(NP(I,2))-XORD(NP(I,1));
         YJ=YORD(NP(I,2))-YORD(NP(I,1));
         XK=XORD(NP(I,3))-XORD(NP(I,1));
         YK=YORD(NP(I,3))-YORD(NP(I,1));
         XC=XORD(NP(I,1))+(XJ+XK)/3.0;
         YC=YORD(NP(I,1))+(YJ+YK)/3.0;

         VOL=(XJ*YK-XK*YJ)/2.0 ;
         if VOL < 0.0
	    error(' Element VOL is less than zero ')
	 end

%        --------------------------
%        INCLUDE user's COEF.m code
%        --------------------------
         COEF

         COMM=1.0/(2.0*VOL);

         DNDX(1)=-(YK-YJ)*COMM; 
         DNDX(2)=+(YK   )*COMM;
         DNDX(3)=-(YJ   )*COMM;
         DNDY(1)=+(XK-XJ)*COMM;
         DNDY(2)=-(XK   )*COMM;
         DNDY(3)=+(XJ   )*COMM;

         S = ((DNDX')*RXI*(DNDX)+(DNDY')*RYI*(DNDY))*VOL;
	 Qe=(VOL*QVI/3.0)*[1;1;1];

%        ---------------------------------
%        Place in Global SK and Q matrices
%        ---------------------------------
         SK(rc,rc)=SK(rc,rc)+S;
         RHS(rc)=RHS(rc)+Qe;

      end

%     ---------------------------
%     Specify boundary condtions
%     ---------------------------
      for I=1:NUMNP;
        if NPBC(I) == 1
           SK(I,I)=SK(I,I)*1.0E+10;  
           RHS(I)=PHI(I)*SK(I,I);
        end
      end

%     ------------------------------------
%     Solution of Finite Element Equations
%     ------------------------------------
      PHI = SK\RHS;

%     -------------------
%     Save results
%     -------------------
      save PHI PHI -ASCII

