[976] | 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 | |
---|
| 22 | Model 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 | |
---|
| 117 | end |
---|