source: trunk/sample/miscellaneous/ocfem/MODEL_FE.mso @ 976

#*---------------------------------------------------------------------
* EMSO Model Library (EML) Copyright (C) 2004 - 2016 ALSOC.
*
* This LIBRARY is free software; you can distribute it and/or modify
* it under the therms of the ALSOC FREE LICENSE as available at
* http://www.enq.ufrgs.br/alsoc.
*
* EMSO Copyright (C) 2004 - 2016 ALSOC, original code
* from http://www.rps.eng.br Copyright (C) 2002-2004.
* All rights reserved.
*
* EMSO is distributed under the therms of the ALSOC LICENSE as
* available at http://www.enq.ufrgs.br/alsoc.
*
*----------------------------------------------------------------------
* Author: Argimiro R. Secchi
* COPPE/UFRJ - Group of Modeling, Simulation, Control,
*                               and Optimization of Processes
* $Id$
*--------------------------------------------------------------------*#

Model MCO_EF
                
        PARAMETERS
        MCO as Plugin(Type="OCFEM",Boundary="BOTH", 
                                  InternalPoints=4, 
                                  alfa=1, beta=1);
        ne as Integer(Brief="number of elements");
        np as Integer;
        N as Integer;
        A(N) as Real(Brief="column-wise stacked matrix A");
        B(N) as Real(Brief="column-wise stacked matrix B");
        r(np+1) as Real(Brief="roots of Jacobi");
                
        SET
        
        ne = 4; #Number of finite elements
        
        np = MCO.NodalPoints;
        N = np * np;
        A = MCO.matrixA;        
        B = MCO.matrixB;
        r(1:np) = MCO.roots;

        VARIABLES
        h(ne+1) as Real(Brief="elements boundaries");
        
        rr(np,ne) as Real; # roots in termos of the whole domain
        s(ne*np-(ne-1)+1) as Real(Brief="Indep. var. for Lagrange interp.");

        y(np,ne) as Real(Brief="dependent vasriable");
        resp(ne*np-(ne-1)) as Real(Brief="final answer");
        
        mA(np,np) as Real(Brief="Matrix A");
        mB(np,np) as Real(Brief="Matrix B");
        
        dif1x(np,ne) as Real(Brief="First derivative of y");
        dif2x(np,ne) as Real(Brief="Second derivative of y");   
        
        EQUATIONS

        #mapping vectors into matrices
                for j in [1:np] do
                        for i in [1:np] do
                                mA(i,j) = A(i+(j-1)*np);
                                mB(i,j) = B(i+(j-1)*np);
                        end
                end
        
        #Derivatives:
        for k in [1:ne] do
                for i in [1:np] do
                        dif1x(i,k)=sum(mA(i,:)*y(:,k)) * 1/(h(k+1)-h(k));
                        dif2x(i,k)=sum(mB(i,:)*y(:,k)) * 1/((h(k+1)-h(k))^2);
                end
        end
        
        #Continuity
        for k in [1:ne-1] do
                y(np,k) = y(1,k+1);
                dif1x(np,k) = dif1x(1,k+1);
        end
        
#-------------------------------------------------------------------------------------- 

        #roots in terms of the whole domain
        for k in [1:ne-1] do
                for i in [1:np] do
                        rr(i,k) = r(i) * (h(k+1)-h(k)) + h(k);#rr results in rows
                end
        end
        for i in [1:np] do
                rr(i,ne) = r(i) * (h(ne+1)-h(ne)) + h(ne);
        end
        
        #Independent variable in sequence
        for i in [1:np] do
                s(i) = rr(i,1);
        end
        for j in [2:ne] do
                for i in [1:np-1] do
                                s((j-1)*np+i-(j-2)) = rr(i+1,j);
                end
        end
        s(ne*np-(ne-1)+1) = h(ne+1);
        
        #Full answer in sequence
        for i in [1:np] do
                resp(i) = y(i,1);
        end
        for j in [2:ne] do
                for i in [1:np-1] do
                                resp((j-1)*np+i-(j-2)) = y(i+1,j);
                end
        end
        
end