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

Last change on this file was 976, checked in by Argimiro Resende Secchi, 6 years ago

Examples for polynomial approximation.

File size: 3.0 KB
Line 
1#*---------------------------------------------------------------------
2* EMSO Model Library (EML) Copyright (C) 2004 - 2016 ALSOC.
3*
4* This LIBRARY is free software; you can distribute it and/or modify
5* it under the therms of the ALSOC FREE LICENSE as available at
6* http://www.enq.ufrgs.br/alsoc.
7*
8* EMSO Copyright (C) 2004 - 2016 ALSOC, original code
9* from http://www.rps.eng.br Copyright (C) 2002-2004.
10* All rights reserved.
11*
12* EMSO is distributed under the therms of the ALSOC LICENSE as
13* available at http://www.enq.ufrgs.br/alsoc.
14*
15*----------------------------------------------------------------------
16* Author: Argimiro R. Secchi
17* COPPE/UFRJ - Group of Modeling, Simulation, Control,
18*                               and Optimization of Processes
19* $Id$
20*--------------------------------------------------------------------*#
21
22Model MCO_EF
23               
24        PARAMETERS
25        MCO as Plugin(Type="OCFEM",Boundary="BOTH",
26                                  InternalPoints=4,
27                                  alfa=1, beta=1);
28        ne as Integer(Brief="number of elements");
29        np as Integer;
30        N as Integer;
31        A(N) as Real(Brief="column-wise stacked matrix A");
32        B(N) as Real(Brief="column-wise stacked matrix B");
33        r(np+1) as Real(Brief="roots of Jacobi");
34               
35        SET
36       
37        ne = 4; #Number of finite elements
38       
39        np = MCO.NodalPoints;
40        N = np * np;
41        A = MCO.matrixA;       
42        B = MCO.matrixB;
43        r(1:np) = MCO.roots;
44
45        VARIABLES
46        h(ne+1) as Real(Brief="elements boundaries");
47       
48        rr(np,ne) as Real; # roots in termos of the whole domain
49        s(ne*np-(ne-1)+1) as Real(Brief="Indep. var. for Lagrange interp.");
50
51        y(np,ne) as Real(Brief="dependent vasriable");
52        resp(ne*np-(ne-1)) as Real(Brief="final answer");
53       
54        mA(np,np) as Real(Brief="Matrix A");
55        mB(np,np) as Real(Brief="Matrix B");
56       
57        dif1x(np,ne) as Real(Brief="First derivative of y");
58        dif2x(np,ne) as Real(Brief="Second derivative of y");   
59       
60        EQUATIONS
61
62        #mapping vectors into matrices
63                for j in [1:np] do
64                        for i in [1:np] do
65                                mA(i,j) = A(i+(j-1)*np);
66                                mB(i,j) = B(i+(j-1)*np);
67                        end
68                end
69       
70        #Derivatives:
71        for k in [1:ne] do
72                for i in [1:np] do
73                        dif1x(i,k)=sum(mA(i,:)*y(:,k)) * 1/(h(k+1)-h(k));
74                        dif2x(i,k)=sum(mB(i,:)*y(:,k)) * 1/((h(k+1)-h(k))^2);
75                end
76        end
77       
78        #Continuity
79        for k in [1:ne-1] do
80                y(np,k) = y(1,k+1);
81                dif1x(np,k) = dif1x(1,k+1);
82        end
83       
84#--------------------------------------------------------------------------------------
85
86        #roots in terms of the whole domain
87        for k in [1:ne-1] do
88                for i in [1:np] do
89                        rr(i,k) = r(i) * (h(k+1)-h(k)) + h(k);#rr results in rows
90                end
91        end
92        for i in [1:np] do
93                rr(i,ne) = r(i) * (h(ne+1)-h(ne)) + h(ne);
94        end
95       
96        #Independent variable in sequence
97        for i in [1:np] do
98                s(i) = rr(i,1);
99        end
100        for j in [2:ne] do
101                for i in [1:np-1] do
102                                s((j-1)*np+i-(j-2)) = rr(i+1,j);
103                end
104        end
105        s(ne*np-(ne-1)+1) = h(ne+1);
106       
107        #Full answer in sequence
108        for i in [1:np] do
109                resp(i) = y(i,1);
110        end
111        for j in [2:ne] do
112                for i in [1:np-1] do
113                                resp((j-1)*np+i-(j-2)) = y(i+1,j);
114                end
115        end
116       
117end
Note: See TracBrowser for help on using the repository browser.