Mercurial > repos > public > sbplib
changeset 349:cd6a29ab3746 feature/hypsyst
A 3D is added and an attempt to imlement 3D transfinit interpolation has been initialized
| author | Ylva Rydin <ylva.rydin@telia.com> |
|---|---|
| date | Thu, 13 Oct 2016 09:34:30 +0200 |
| parents | d9860ebc3148 |
| children | 5d5652fe826a |
| files | +grid/Ti3D.m +scheme/Hypsyst3d.m |
| diffstat | 2 files changed, 570 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/Ti3D.m Thu Oct 13 09:34:30 2016 +0200 @@ -0,0 +1,252 @@ +classdef Ti3D + properties + gs % {6}Surfaces + V % FunctionHandle(XI,ETA,ZETA) + end + + methods + % TODO function to label boundary names. + % function to find largest and smallest delta h in the grid. Maybe shouldnt live here + function obj = Ti3D(CW,CE,CS,CN,CB,CT) + obj.gs = {CE,CW,CS,CN,CB,CT}; + + gw = CW.g; + ge = CE.g; + gs = CS.g; + gn = CN.g; + gb = CB.g; + gt = CT.g; + + function o = V_fun(XI,ETA,ZETA) + XI=XI'; + ETA=ETA'; + ZETA=ZETA'; + + one=0*ETA+1; + zero=0*ETA; + + Sw = gw((1-ETA),(1-ZETA)); + Se = ge(ETA,ZETA); + Ss = gs(XI,(1-ZETA)); + Sn = gn((1-XI),ZETA); + Sb = gb(XI,ETA); + St = gt((1-XI),(1-ETA)); + + Ewt = gw(1-ETA,zero); + Ewb = gw(1-ETA,one); + Ews = gw(one,1-ZETA); + Ewn = gw(zero,1-ZETA); + Eet = ge(ETA,one); + Eeb = ge(ETA,zero); + Ees = ge(0*one,ZETA); + Een = ge(one,ZETA); + Enb = gn(1-XI,zero); + Ent = gn(1-XI,one); + Est = gs(XI,zero); + Esb = gs(XI,one); + + Cwbs = gw(one,one); + Cwbn = gw(zero,one); + Cwts = gw(one,zero); + Cwtn = gw(zero,zero); + Cebs = ge(zero,zero); + Cebn = ge(one,zero); + Cets = ge(zero,one); + Cetn = ge(one,one); + + + X1 = (1-XI).*Sw(1,:,:) + XI.*Se(1,:,:); + X2 = (1-ETA).*Ss(1,:,:) + ETA.*Sn(1,:,:); + X3 = (1-ZETA).*Sb(1,:,:) + ZETA.*St(1,:,:); + + X12 = (1-XI).*(1-ETA).*Ews(1,:,:) + (1-XI).*ETA.*Ewn(1,:,:) + XI.*(1-ETA).*Ees(1,:,:) + XI.*ETA.*Een(1,:,:); + X13 = (1-XI).*(1-ZETA).*Ewb(1,:,:) + (1-XI).*ZETA.*Ewt(1,:,:) + XI.*(1-ZETA).*Eeb(1,:,:) + XI.*ZETA.*Eet(1,:,:); + X23 = (1-ETA).*(1-ZETA).*Esb(1,:,:) + (1-ETA).*ZETA.*Est(1,:,:) + ETA.*(1-ZETA).*Enb(1,:,:) + ETA.*ZETA.*Ent(1,:,:); + + X123 = (1-XI).*(1-ETA).*(1-ZETA).*Cwbs(1,:,:) + (1-XI).*(1-ETA).*ZETA.*Cwts(1,:,:) + (1-XI).*ETA.*(1-ZETA).*Cwbn(1,:,:) + ... + (1-XI).*ETA.*ZETA.*Cwtn(1,:,:) + XI.*(1-ETA).*(1-ZETA).*Cebs(1,:,:) + XI.*(1-ETA).*ZETA.*Cets(1,:,:) + ... + XI.*ETA.*(1-ZETA).*Cebn(1,:,:) + XI.*ETA.*ZETA.*Cetn(1,:,:); + + X = X1 + X2 + X3 - X12 - X13 - X23 + X123; + + + Y1 = (1-XI).*Sw(2,:,:) + XI.*Se(2,:,:); + Y2 = (1-ETA).*Ss(2,:,:) + ETA.*Sn(2,:,:); + Y3 = (1-ZETA).*Sb(2,:,:) + ZETA.*St(2,:,:); + + Y12 = (1-XI).*(1-ETA).*Ews(2,:,:) + (1-XI).*ETA.*Ewn(2,:,:) + XI.*(1-ETA).*Ees(2,:,:) + XI.*ETA.*Een(2,:,:); + Y13 = (1-XI).*(1-ZETA).*Ewb(2,:,:) + (1-XI).*ZETA.*Ewt(2,:,:) + XI.*(1-ZETA).*Eeb(2,:,:) + XI.*ZETA.*Eet(2,:,:); + Y23 = (1-ETA).*(1-ZETA).*Esb(2,:,:) + (1-ETA).*ZETA.*Est(2,:,:) + ETA.*(1-ZETA).*Enb(2,:,:) + ETA.*ZETA.*Ent(2,:,:); + + Y123 = (1-XI).*(1-ETA).*(1-ZETA).*Cwbs(2,:,:) + (1-XI).*(1-ETA).*ZETA.*Cwts(2,:,:) + (1-XI).*ETA.*(1-ZETA).*Cwbn(2,:,:) + ... + (1-XI).*ETA.*ZETA.*Cwtn(2,:,:) + XI.*(1-ETA).*(1-ZETA).*Cebs(2,:,:) + XI.*(1-ETA).*ZETA.*Cets(2,:,:) + ... + XI.*ETA.*(1-ZETA).*Cebn(2,:,:) + XI.*ETA.*ZETA.*Cetn(2,:,:); + + Y = Y1 + Y2 + Y3 - Y12 - Y13 - Y23 + Y123; + + + Z1 = (1-XI).*Sw(3,:,:) + XI.*Se(3,:,:); + Z2 = (1-ETA).*Ss(3,:,:) + ETA.*Sn(3,:,:); + Z3 = (1-ZETA).*Sb(3,:,:) + ZETA.*St(3,:,:); + + Z12 = (1-XI).*(1-ETA).*Ews(3,:,:) + (1-XI).*ETA.*Ewn(3,:,:) + XI.*(1-ETA).*Ees(3,:,:) + XI.*ETA.*Een(3,:,:); + Z13 = (1-XI).*(1-ZETA).*Ewb(3,:,:) + (1-XI).*ZETA.*Ewt(3,:,:) + XI.*(1-ZETA).*Eeb(3,:,:) + XI.*ZETA.*Eet(3,:,:); + Z23 = (1-ETA).*(1-ZETA).*Esb(3,:,:) + (1-ETA).*ZETA.*Est(3,:,:) + ETA.*(1-ZETA).*Enb(3,:,:) + ETA.*ZETA.*Ent(3,:,:); + + Z123 = (1-XI).*(1-ETA).*(1-ZETA).*Cwbs(3,:,:) + (1-XI).*(1-ETA).*ZETA.*Cwts(3,:,:) + (1-XI).*ETA.*(1-ZETA).*Cwbn(3,:,:) + ... + (1-XI).*ETA.*ZETA.*Cwtn(3,:,:) + XI.*(1-ETA).*(1-ZETA).*Cebs(3,:,:) + XI.*(1-ETA).*ZETA.*Cets(3,:,:) + ... + XI.*ETA.*(1-ZETA).*Cebn(3,:,:) + XI.*ETA.*ZETA.*Cetn(3,:,:); + + Z = Z1 + Z2 + Z3 - Z12 - Z13 - Z23 + Z123; + o = [X;Y;Z]; + end + + obj.V = @V_fun; + end + + + function [X,Y,Z] = map(obj,XI,ETA,ZETA) + + V = obj.V; + + p = V(XI,ETA,ZETA); + X = p(1,:)'; + Y = p(2,:)'; + Z = p(3,:)'; + + end + + % function h = plot(obj,nu,nv) + % S = obj.S; + % + % default_arg('nv',nu) + % + % u = linspace(0,1,nu); + % v = linspace(0,1,nv); + % + % m = 100; + % + % X = zeros(nu+nv,m); + % Y = zeros(nu+nv,m); + % + % + % t = linspace(0,1,m); + % for i = 1:nu + % p = S(u(i),t); + % X(i,:) = p(1,:); + % Y(i,:) = p(2,:); + % end + % + % for i = 1:nv + % p = S(t,v(i)); + % X(i+nu,:) = p(1,:); + % Y(i+nu,:) = p(2,:); + % end + % + % h = line(X',Y'); + % end + % + % + % function h = show(obj,nu,nv) + % default_arg('nv',nu) + % S = obj.S; + % + % if(nu>2 || nv>2) + % h_grid = obj.plot(nu,nv); + % set(h_grid,'Color',[0 0.4470 0.7410]); + % end + % + % h_bord = obj.plot(2,2); + % set(h_bord,'Color',[0.8500 0.3250 0.0980]); + % set(h_bord,'LineWidth',2); + % end + % + % + % % TRANSFORMATIONS + % function ti = translate(obj,a) + % gs = obj.gs; + % + % for i = 1:length(gs) + % new_gs{i} = gs{i}.translate(a); + % end + % + % ti = grid.Ti(new_gs{:}); + % end + % + % % Mirrors the Ti so that the resulting Ti is still left handed. + % % (Corrected by reversing curves and switching e and w) + % function ti = mirror(obj, a, b) + % gs = obj.gs; + % + % new_gs = cell(1,4); + % + % new_gs{1} = gs{1}.mirror(a,b).reverse(); + % new_gs{3} = gs{3}.mirror(a,b).reverse(); + % new_gs{2} = gs{4}.mirror(a,b).reverse(); + % new_gs{4} = gs{2}.mirror(a,b).reverse(); + % + % ti = grid.Ti(new_gs{:}); + % end + % + % function ti = rotate(obj,a,rad) + % gs = obj.gs; + % + % for i = 1:length(gs) + % new_gs{i} = gs{i}.rotate(a,rad); + % end + % + % ti = grid.Ti(new_gs{:}); + % end + % + % function ti = rotate_edges(obj,n); + % new_gs = cell(1,4); + % for i = 0:3 + % new_i = mod(i - n,4); + % new_gs{new_i+1} = obj.gs{i+1}; + % end + % ti = grid.Ti(new_gs{:}); + % end + % end + % + % methods(Static) + % function obj = points(p1, p2, p3, p4) + % g1 = grid.Curve.line(p1,p2); + % g2 = grid.Curve.line(p2,p3); + % g3 = grid.Curve.line(p3,p4); + % g4 = grid.Curve.line(p4,p1); + % + % obj = grid.Ti(g1,g2,g3,g4); + % end + % + % function label(varargin) + % if nargin == 2 && ischar(varargin{2}) + % label_impl(varargin{:}); + % else + % for i = 1:length(varargin) + % label_impl(varargin{i},inputname(i)); + % end + % end + % + % + % function label_impl(ti,str) + % S = ti.S; + % + % pc = S(0.5,0.5); + % + % margin = 0.1; + % pw = S( margin, 0.5); + % pe = S(1-margin, 0.5); + % ps = S( 0.5, margin); + % pn = S( 0.5, 1-margin); + % + % + % ti.show(2,2); + % grid.place_label(pc,str); + % grid.place_label(pw,'w'); + % grid.place_label(pe,'e'); + % grid.place_label(ps,'s'); + % grid.place_label(pn,'n'); + % end + % end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Hypsyst3d.m Thu Oct 13 09:34:30 2016 +0200 @@ -0,0 +1,318 @@ +classdef Hypsyst3d < scheme.Scheme + properties + m % Number of points in each direction, possibly a vector + n %size of system + h % Grid spacing + x, y, z % Grid + X, Y, Z% Values of x and y for each grid point + Yx, Zx, Xy, Zy, Xz, Yz %Grid values for boundary surfaces + order % Order accuracy for the approximation + + D % non-stabalized scheme operator + A, B, C, E + + H % Discrete norm + % Norms in the x, y and z directions + Hxi,Hyi, Hzi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. + I_x,I_y, I_z, I_N + e_w, e_e, e_s, e_n, e_b, e_t + params %parameters for the coeficient matrice + end + + + methods + function obj = Hypsyst3d(m, lim, order, A, B,C, E, params) + default_arg('E', []) + xlim = lim{1}; + ylim = lim{2}; + zlim = lim{3}; + + if length(m) == 1 + m = [m m m]; + end + + obj.A=A; + obj.B=B; + obj.C=C; + obj.E=E; + m_x = m(1); + m_y = m(2); + m_z=m(3); + obj.params = params; + + ops_x = sbp.D2Standard(m_x,xlim,order); + ops_y = sbp.D2Standard(m_y,ylim,order); + ops_z = sbp.D2Standard(m_z,zlim,order); + + obj.x = ops_x.x; + obj.y = ops_y.x; + obj.z = ops_z.x; + + obj.X = kr(obj.x,ones(m_y,1),ones(m_z,1));%% Que pasa? + obj.Y = kr(ones(m_x,1),obj.y,ones(m_z,1)); + obj.Z = kr(ones(m_x,1),ones(m_y,1),obj.z); + + obj.Yx=kr(obj.y,ones(m_z,1)); + obj.Zx=kr(ones(m_y,1),obj.z); + + obj.Xy=kr(obj.x,ones(m_z,1)); + obj.Zy=kr(ones(m_x,1),obj.z); + + obj.Xz=kr(obj.x,ones(m_y,1)); + obj.Yz=kr(ones(m_z,1),obj.y); + + Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y,obj.Z); + Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y,obj.Z); + Cevaluated = obj.evaluateCoefficientMatrix(C, obj.X, obj.Y,obj.Z); + Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y,obj.Z); + + obj.n = length(A(obj.params,0,0,0)); + + I_n = eye(obj.n); + I_x = speye(m_x); + obj.I_x = I_x; + I_y = speye(m_y); + obj.I_y = I_y; + I_z = speye(m_z); + obj.I_z = I_z; + + + D1_x = kr(I_n, ops_x.D1, I_y,I_z); + obj.Hxi = kr(I_n, ops_x.HI, I_y,I_z); + D1_y = kr(I_n, I_x, ops_y.D1,I_z); + obj.Hyi = kr(I_n, I_x, ops_y.HI,I_z); + D1_z = kr(I_n, I_x, I_y,ops_z.D1); + obj.Hzi = kr(I_n, I_x,I_y, ops_y.HI); + + obj.e_w = kr(I_n, ops_x.e_l, I_y,I_z); + obj.e_e = kr(I_n, ops_x.e_r, I_y,I_z); + obj.e_s = kr(I_n, I_x, ops_y.e_l,I_z); + obj.e_n = kr(I_n, I_x, ops_y.e_r,I_z); + obj.e_b = kr(I_n, I_x, I_y, ops_z.e_l); + obj.e_t = kr(I_n, I_x, I_y, ops_z.e_r); + + obj.m=m; + obj.h=[ops_x.h ops_y.h ops_x.h]; + obj.order=order; + + obj.D=-Aevaluated*D1_x-Bevaluated*D1_y-Cevaluated*D1_z-Eevaluated; + end + + % Closure functions return the opertors applied to the own doamin to close the boundary + % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. + % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. + % type is a string specifying the type of boundary condition if there are several. + % data is a function returning the data that should be applied at the boundary. + function [closure, penalty] = boundary_condition(obj,boundary,type,L) + default_arg('type','char'); + BM=boundary_matrices(obj,boundary); + + switch type + case{'c','char'} + [closure,penalty]=boundary_condition_char(obj,BM); + case{'general'} + [closure,penalty]=boundary_condition_general(obj,BM,boundary,L); + otherwise + error('No such boundary condition') + end + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + error('An interface function does not exist yet'); + end + + function N = size(obj) + N = obj.m; + end + + function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y, Z) + params=obj.params; + side=max(length(X),length(Y)); + if isa(mat,'function_handle') + [rows,cols]=size(mat(params,0,0,0)); + matVec=mat(params,X',Y',Z'); + matVec=sparse(matVec); + else + matVec=mat; + [rows,cols]=size(matVec); + side=max(length(X),length(Y)); + cols=cols/side; + end + ret=kron(ones(rows,cols),speye(side)); + + for ii=1:rows + for jj=1:cols + ret((ii-1)*side+1:ii*side,(jj-1)*side+1:jj*side)=diag(matVec(ii,(jj-1)*side+1:jj*side)); + end + end + end + + + function [BM]=boundary_matrices(obj,boundary) + params=obj.params; + + switch boundary + case {'w','W','west'} + BM.e_=obj.e_w; + mat=obj.A; + BM.boundpos='l'; + BM.Hi=obj.Hxi; + [BM.V,BM.Vi,BM.D,signVec]=obj.matrixDiag(mat,obj.X(1),obj.Yx,obj.Zx); + BM.side=length(obj.Yx); + case {'e','E','east'} + BM.e_=obj.e_e; + mat=obj.A; + BM.boundpos='r'; + BM.Hi=obj.Hxi; + [BM.V,BM.Vi,BM.D,signVec]=obj.matrixDiag(mat,obj.X(end),obj.Yx,obj.Zx); + BM.side=length(obj.Yx); + case {'s','S','south'} + BM.e_=obj.e_s; + mat=obj.B; + BM.boundpos='l'; + BM.Hi=obj.Hyi; + [BM.V,BM.Vi,BM.D,signVec]=obj.matrixDiag(mat,obj.Xy,obj.Y(1),obj.Zy); + BM.side=length(obj.Xy); + case {'n','N','north'} + BM.e_=obj.e_n; + mat=obj.B; + BM.boundpos='r'; + BM.Hi=obj.Hyi; + [BM.V,BM.Vi,BM.D,signVec]=obj.matrixDiag(mat,obj.Xy,obj.Y(end),obj.Zy); + BM.side=length(obj.Xy); + case{'b','B','Bottom'} + BM.e_=obj.e_b; + mat=obj.C; + BM.boundpos='l'; + BM.Hi=obj.Hzi; + [BM.V,BM.Vi,BM.D,signVec]=obj.matrixDiag(mat,obj.Xz,obj.Yz,obj.Z(1)); + BM.side=length(obj.Xz); + case{'t','T','Top'} + BM.e_=obj.e_t; + mat=obj.C; + BM.boundpos='r'; + BM.Hi=obj.Hzi; + [BM.V,BM.Vi,BM.D,signVec]=obj.matrixDiag(mat,obj.Xz,obj.Yz,obj.Z(end)); + BM.side=length(obj.Xz); + end + + BM.pos=signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3); + end + + + function [closure, penalty]=boundary_condition_char(obj,BM) + side = BM.side; + pos = BM.pos; + neg = BM.neg; + zeroval=BM.zeroval; + V = BM.V; + Vi = BM.Vi; + Hi=BM.Hi; + D=BM.D; + e_=BM.e_; + + switch BM.boundpos + case {'l'} + tau=sparse(obj.n*side,pos); + Vi_plus=Vi(1:pos,:); + tau(1:pos,:)=-abs(D(1:pos,1:pos)); + closure=Hi*e_*V*tau*Vi_plus*e_'; + penalty=-Hi*e_*V*tau*Vi_plus; + case {'r'} + tau=sparse(obj.n*side,neg); + tau((pos+zeroval)+1:obj.n*side,:)=-abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); + Vi_minus=Vi((pos+zeroval)+1:obj.n*side,:); + closure=Hi*e_*V*tau*Vi_minus*e_'; + penalty=-Hi*e_*V*tau*Vi_minus; + end + end + + + function [closure,penalty]=boundary_condition_general(obj,BM,boundary,L) + side = BM.side; + pos = BM.pos; + neg = BM.neg; + zeroval=BM.zeroval; + V = BM.V; + Vi = BM.Vi; + Hi=BM.Hi; + D=BM.D; + e_=BM.e_; + switch boundary + case {'w','W','west'} + L=obj.evaluateCoefficientMatrix(L,obj.x(1),obj.Yx,obj.Zx); + case {'e','E','east'} + L=obj.evaluateCoefficientMatrix(L,obj.x(end),obj.Yx,obj.Zx); + case {'s','S','south'} + L=obj.evaluateCoefficientMatrix(L,obj.Xy,obj.y(1),obj.Zy); + case {'n','N','north'} + L=obj.evaluateCoefficientMatrix(L,obj.Xy,obj.y(end),obj.Zy); + case {'b','B','bottom'} + L=obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(1)); + case {'t','T','top'} + L=obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(end)); + end + + switch BM.boundpos + case {'l'} + tau=sparse(obj.n*side,pos); + Vi_plus=Vi(1:pos,:); + Vi_minus=Vi(pos+zeroval+1:obj.n*side,:); + V_plus=V(:,1:pos); + V_minus=V(:,(pos+zeroval)+1:obj.n*side); + + tau(1:pos,:)=-abs(D(1:pos,1:pos)); + R=-inv(L*V_plus)*(L*V_minus); + closure=Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; + penalty=-Hi*e_*V*tau*inv(L*V_plus)*L; + case {'r'} + tau=sparse(obj.n*side,neg); + tau((pos+zeroval)+1:obj.n*side,:)=-abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); + Vi_plus=Vi(1:pos,:); + Vi_minus=Vi((pos+zeroval)+1:obj.n*side,:); + + V_plus=V(:,1:pos); + V_minus=V(:,(pos+zeroval)+1:obj.n*side); + R=-inv(L*V_minus)*(L*V_plus); + closure=Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_'; + penalty=-Hi*e_*V*tau*inv(L*V_minus)*L; + end + end + + + function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y,z) + params=obj.params; + syms xs ys zs + [V, D]=eig(mat(params,xs,ys,zs)); + xs=x; + ys=y; + zs=z; + + + side=max(length(x),length(y)); + Dret=zeros(obj.n,side*obj.n); + Vret=zeros(obj.n,side*obj.n); + for ii=1:obj.n + for jj=1:obj.n + Dret(jj,(ii-1)*side+1:side*ii)=eval(D(jj,ii)); + Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii)); + end + end + + D=sparse(Dret); + V=sparse(Vret); + V=obj.evaluateCoefficientMatrix(V,x,y,z); + D=obj.evaluateCoefficientMatrix(D,x,y,z); + DD=diag(D); + + poseig=(DD>0); + zeroeig=(DD==0); + negeig=(DD<0); + + D=diag([DD(poseig); DD(zeroeig); DD(negeig)]); + V=[V(:,poseig) V(:,zeroeig) V(:,negeig)]; + Vi=inv(V); + signVec=[sum(poseig),sum(zeroeig),sum(negeig)]; + end + end +end \ No newline at end of file