Mercurial > forge
changeset 3636:90d9933cbf2c octave-forge
Rename run.m example files to run_test.m so that indexing will ignore them
| author | adb014 |
|---|---|
| date | Wed, 11 Jul 2007 20:27:45 +0000 |
| parents | 3aea065ad216 |
| children | 7a498e61660c |
| files | extra/bim/examples/COMSON/run.m extra/bim/examples/COMSON/run_test.m extra/bim/examples/FIUME/run.m extra/bim/examples/FIUME/run_test.m extra/bim/examples/SQUARE/run.m extra/bim/examples/SQUARE/run_test.m extra/bim/examples/SSHAPE/run.m extra/bim/examples/SSHAPE/test.m |
| diffstat | 8 files changed, 353 insertions(+), 353 deletions(-) [+] |
line wrap: on
line diff
--- a/extra/bim/examples/COMSON/run.m Wed Jul 11 20:21:26 2007 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,89 +0,0 @@ -## Mesh generation and properties -[mesh] = MSH2Mgmsh('comson',1,20); -[mesh] = BIM2Cmeshproperties(mesh); - -## Construct initial guess -[xu, yu] = BIM2Cunknowncoord(mesh); -nelem = size(mesh.t,2); -nnodi = size(mesh.p,2); -uin = 0*xu; - -## set the coefficients for the problem: -## -div ( \alpha \gamma ( \eta \nabla u - \beta u ) )+ \delta \zeta u = f g -## FIX ME: CAMBIARE I VALORI CON QUALCOSA CHE DIA DEI DISEGNINI PIU' BELLI! -epsilon = .1; -alfa = ones(nelem,1); -gamma = ones(nnodi,1); -eta = epsilon*ones(nnodi,1); -beta = xu+yu; -delta = .1*ones(nelem,1); -zeta = .1*ones(nnodi,1); -f = ones(nnodi,1); -g = ones(nelem,1); - -## Construction of the discretized operators -AdvDiff = BIM2Aadvdiff(mesh,alfa,gamma,eta,beta); -Mass = BIM2Areaction(mesh,delta,zeta); -b = BIM2Arhs(mesh,f,g); - - -A = AdvDiff + Mass; - - -## BOUNDARY CONDITIONS -Dlist = BIM2Cunknownsonside(mesh, [55 38 56 57 35 54,... - 29 30 1:17 24 25 31 32,... - 33 34 18:23 26:28]); -## DIRICHLET LIST -Nlist = [] ## NEUMANN LIST -Nlist = setdiff(Nlist,Dlist); -Fn = zeros(length(Nlist),1); ## PRESCRIBED NEUMANN FLUXES -Ilist = setdiff(1:length(uin),union(Dlist,Nlist)); ## INTERNAL LIST - -## Partition the LHS and RHS - -Add = A(Dlist,Dlist); -Adn = A(Dlist,Nlist); ## shoud be all zeros hopefully!! -Adi = A(Dlist,Ilist); - -And = A(Nlist,Dlist); ## shoud be all zeros hopefully!! -Ann = A(Nlist,Nlist); -Ani = A(Nlist,Ilist); - -Aid = A(Ilist,Dlist); -Ain = A(Ilist,Nlist); -Aii = A(Ilist,Ilist); - -bd = b(Dlist); -bn = b(Nlist); -bi = b(Ilist); - -ud = uin(Dlist); -un = uin(Nlist); -ui = uin(Ilist); - -## SOLUTION -temp = [Ann Ani ; Ain Aii ] \ [ Fn+bn-And*ud ; bi-Aid*ud]; -un = temp(1:length(un)); -ui = temp(length(un)+1:end); - -## Compute fluxes through Dirichlet sides - -Fd = Add * ud + Adi * ui + Adn*un - bd; - -u = zeros(nnodi,1); -u(Dlist) = ud; u(Nlist) = un; u(Ilist) = ui; - -## Reconstruction of the currents - -## GRADIENT OF THE SOLUTION -[gx, gy] = BIM2Cpdegrad(mesh,u); - - -## FLUXES -[jxglob,jyglob] = BIM2Cglobalflux(mesh,u,alfa,gamma,eta,beta); - -## VISUALIZATION -FPL2pdesurf(mesh,u); -FPL2quiver(mesh,jxglob,jyglob,"sample_density","10000"); -FPL2quiver(mesh,gx,gy,"sample_density","10000"); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/bim/examples/COMSON/run_test.m Wed Jul 11 20:27:45 2007 +0000 @@ -0,0 +1,89 @@ +## Mesh generation and properties +[mesh] = MSH2Mgmsh('comson',1,20); +[mesh] = BIM2Cmeshproperties(mesh); + +## Construct initial guess +[xu, yu] = BIM2Cunknowncoord(mesh); +nelem = size(mesh.t,2); +nnodi = size(mesh.p,2); +uin = 0*xu; + +## set the coefficients for the problem: +## -div ( \alpha \gamma ( \eta \nabla u - \beta u ) )+ \delta \zeta u = f g +## FIX ME: CAMBIARE I VALORI CON QUALCOSA CHE DIA DEI DISEGNINI PIU' BELLI! +epsilon = .1; +alfa = ones(nelem,1); +gamma = ones(nnodi,1); +eta = epsilon*ones(nnodi,1); +beta = xu+yu; +delta = .1*ones(nelem,1); +zeta = .1*ones(nnodi,1); +f = ones(nnodi,1); +g = ones(nelem,1); + +## Construction of the discretized operators +AdvDiff = BIM2Aadvdiff(mesh,alfa,gamma,eta,beta); +Mass = BIM2Areaction(mesh,delta,zeta); +b = BIM2Arhs(mesh,f,g); + + +A = AdvDiff + Mass; + + +## BOUNDARY CONDITIONS +Dlist = BIM2Cunknownsonside(mesh, [55 38 56 57 35 54,... + 29 30 1:17 24 25 31 32,... + 33 34 18:23 26:28]); +## DIRICHLET LIST +Nlist = [] ## NEUMANN LIST +Nlist = setdiff(Nlist,Dlist); +Fn = zeros(length(Nlist),1); ## PRESCRIBED NEUMANN FLUXES +Ilist = setdiff(1:length(uin),union(Dlist,Nlist)); ## INTERNAL LIST + +## Partition the LHS and RHS + +Add = A(Dlist,Dlist); +Adn = A(Dlist,Nlist); ## shoud be all zeros hopefully!! +Adi = A(Dlist,Ilist); + +And = A(Nlist,Dlist); ## shoud be all zeros hopefully!! +Ann = A(Nlist,Nlist); +Ani = A(Nlist,Ilist); + +Aid = A(Ilist,Dlist); +Ain = A(Ilist,Nlist); +Aii = A(Ilist,Ilist); + +bd = b(Dlist); +bn = b(Nlist); +bi = b(Ilist); + +ud = uin(Dlist); +un = uin(Nlist); +ui = uin(Ilist); + +## SOLUTION +temp = [Ann Ani ; Ain Aii ] \ [ Fn+bn-And*ud ; bi-Aid*ud]; +un = temp(1:length(un)); +ui = temp(length(un)+1:end); + +## Compute fluxes through Dirichlet sides + +Fd = Add * ud + Adi * ui + Adn*un - bd; + +u = zeros(nnodi,1); +u(Dlist) = ud; u(Nlist) = un; u(Ilist) = ui; + +## Reconstruction of the currents + +## GRADIENT OF THE SOLUTION +[gx, gy] = BIM2Cpdegrad(mesh,u); + + +## FLUXES +[jxglob,jyglob] = BIM2Cglobalflux(mesh,u,alfa,gamma,eta,beta); + +## VISUALIZATION +FPL2pdesurf(mesh,u); +FPL2quiver(mesh,jxglob,jyglob,"sample_density","10000"); +FPL2quiver(mesh,gx,gy,"sample_density","10000"); \ No newline at end of file
--- a/extra/bim/examples/FIUME/run.m Wed Jul 11 20:21:26 2007 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,90 +0,0 @@ -pkg load msh -pkg load fpl -pkg load bim - -## Mesh generation and properties -[mesh] = MSH2Mgmsh('fiume',1,.2); -[mesh] = BIM2Cmeshproperties(mesh); - -## Construct initial guess -[xu, yu] = BIM2Cunknowncoord(mesh); -nelem = size(mesh.t,2); -nnodi = size(mesh.p,2); -uin = 3*xu; - -## set the coefficients for the problem: -## -div ( \alpha \gamma ( \eta \nabla u - \beta u ) )+ \delta \zeta u = f g -## FIX ME: CAMBIARE I VALORI CON QUALCOSA CHE DIA DEI DISEGNINI PIU' BELLI! -epsilon = .1; -alfa = ones(nelem,1); -gamma = ones(nnodi,1); -eta = epsilon*ones(nnodi,1); -beta = xu+yu; -delta = ones(nelem,1); -zeta = ones(nnodi,1); -f = ones(nnodi,1); -g = ones(nelem,1); - -## Construction of the discretized operators -AdvDiff = BIM2Aadvdiff(mesh,alfa,gamma,eta,beta); -Mass = BIM2Areaction(mesh,delta,zeta); -b = BIM2Arhs(mesh,f,g); - - -A = AdvDiff + Mass; - - -## BOUNDARY CONDITIONS -Dlist = BIM2Cunknownsonside(mesh, [8 18]); ## DIRICHLET LIST -Nlist = BIM2Cunknownsonside(mesh, [23 24]); ## NEUMANN LIST -Nlist = setdiff(Nlist,Dlist); -Fn = zeros(length(Nlist),1); ## PRESCRIBED NEUMANN FLUXES -Ilist = setdiff(1:length(uin),union(Dlist,Nlist)); ## INTERNAL LIST - -## Partition the LHS and RHS - -Add = A(Dlist,Dlist); -Adn = A(Dlist,Nlist); ## shoud be all zeros hopefully!! -Adi = A(Dlist,Ilist); - -And = A(Nlist,Dlist); ## shoud be all zeros hopefully!! -Ann = A(Nlist,Nlist); -Ani = A(Nlist,Ilist); - -Aid = A(Ilist,Dlist); -Ain = A(Ilist,Nlist); -Aii = A(Ilist,Ilist); - -bd = b(Dlist); -bn = b(Nlist); -bi = b(Ilist); - -ud = uin(Dlist); -un = uin(Nlist); -ui = uin(Ilist); - -## SOLUTION -temp = [Ann Ani ; Ain Aii ] \ [ Fn+bn-And*ud ; bi-Aid*ud]; -un = temp(1:length(un)); -ui = temp(length(un)+1:end); - -## Compute fluxes through Dirichlet sides - -Fd = Add * ud + Adi * ui + Adn*un - bd; - -u = zeros(nnodi,1); -u(Dlist) = ud; u(Nlist) = un; u(Ilist) = ui; - -## Reconstruction of the currents - -## GRADIENT OF THE SOLUTION -[gx, gy] = BIM2Cpdegrad(mesh,u); - - -## FLUXES -[jxglob,jyglob] = BIM2Cglobalflux(mesh,u,alfa,gamma,eta,beta); - -## VISUALIZATION -FPL2pdesurf(mesh,u); -FPL2quiver(mesh,jxglob,jyglob,"sample_density","1000"); -FPL2quiver(mesh,gx,gy,"sample_density","1000");
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/bim/examples/FIUME/run_test.m Wed Jul 11 20:27:45 2007 +0000 @@ -0,0 +1,90 @@ +pkg load msh +pkg load fpl +pkg load bim + +## Mesh generation and properties +[mesh] = MSH2Mgmsh('fiume',1,.2); +[mesh] = BIM2Cmeshproperties(mesh); + +## Construct initial guess +[xu, yu] = BIM2Cunknowncoord(mesh); +nelem = size(mesh.t,2); +nnodi = size(mesh.p,2); +uin = 3*xu; + +## set the coefficients for the problem: +## -div ( \alpha \gamma ( \eta \nabla u - \beta u ) )+ \delta \zeta u = f g +## FIX ME: CAMBIARE I VALORI CON QUALCOSA CHE DIA DEI DISEGNINI PIU' BELLI! +epsilon = .1; +alfa = ones(nelem,1); +gamma = ones(nnodi,1); +eta = epsilon*ones(nnodi,1); +beta = xu+yu; +delta = ones(nelem,1); +zeta = ones(nnodi,1); +f = ones(nnodi,1); +g = ones(nelem,1); + +## Construction of the discretized operators +AdvDiff = BIM2Aadvdiff(mesh,alfa,gamma,eta,beta); +Mass = BIM2Areaction(mesh,delta,zeta); +b = BIM2Arhs(mesh,f,g); + + +A = AdvDiff + Mass; + + +## BOUNDARY CONDITIONS +Dlist = BIM2Cunknownsonside(mesh, [8 18]); ## DIRICHLET LIST +Nlist = BIM2Cunknownsonside(mesh, [23 24]); ## NEUMANN LIST +Nlist = setdiff(Nlist,Dlist); +Fn = zeros(length(Nlist),1); ## PRESCRIBED NEUMANN FLUXES +Ilist = setdiff(1:length(uin),union(Dlist,Nlist)); ## INTERNAL LIST + +## Partition the LHS and RHS + +Add = A(Dlist,Dlist); +Adn = A(Dlist,Nlist); ## shoud be all zeros hopefully!! +Adi = A(Dlist,Ilist); + +And = A(Nlist,Dlist); ## shoud be all zeros hopefully!! +Ann = A(Nlist,Nlist); +Ani = A(Nlist,Ilist); + +Aid = A(Ilist,Dlist); +Ain = A(Ilist,Nlist); +Aii = A(Ilist,Ilist); + +bd = b(Dlist); +bn = b(Nlist); +bi = b(Ilist); + +ud = uin(Dlist); +un = uin(Nlist); +ui = uin(Ilist); + +## SOLUTION +temp = [Ann Ani ; Ain Aii ] \ [ Fn+bn-And*ud ; bi-Aid*ud]; +un = temp(1:length(un)); +ui = temp(length(un)+1:end); + +## Compute fluxes through Dirichlet sides + +Fd = Add * ud + Adi * ui + Adn*un - bd; + +u = zeros(nnodi,1); +u(Dlist) = ud; u(Nlist) = un; u(Ilist) = ui; + +## Reconstruction of the currents + +## GRADIENT OF THE SOLUTION +[gx, gy] = BIM2Cpdegrad(mesh,u); + + +## FLUXES +[jxglob,jyglob] = BIM2Cglobalflux(mesh,u,alfa,gamma,eta,beta); + +## VISUALIZATION +FPL2pdesurf(mesh,u); +FPL2quiver(mesh,jxglob,jyglob,"sample_density","1000"); +FPL2quiver(mesh,gx,gy,"sample_density","1000");
--- a/extra/bim/examples/SQUARE/run.m Wed Jul 11 20:21:26 2007 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,89 +0,0 @@ -pkg load msh -pkg load fpl -pkg load bim - -## MESH GENERATION AND PROPERTIES -[mesh] = MSH2Mgmsh('quadrato',1,1.2); -[mesh] = BIM2Cmeshproperties(mesh); - -## INITIAL GUESS (IF NEEDED) -[xu, yu] = BIM2Cunknowncoord(mesh); -[nelem] = columns(mesh.t); -[nnodes] = columns(mesh.p); -[uin] = zeros(nnodes,1); - -## COEFFICIENTS FOR THE PROBLEM: -## -div ( \alpha \gamma ( \eta \nabla u - \beta u ) )+ \delta \zeta u = f g -epsilon = .01; -alfa = ones(nelem,1); -gamma = ones(nnodes,1); -eta = epsilon*ones(nnodes,1); -potbeta = xu/epsilon; -delta = zeros(nelem,1); -zeta = zeros(nnodes,1); -f = ones(nnodes,1); -g = ones(nelem,1); - -## ADVECTION AND REACTION MATRIX, RHS -AdvDiff = BIM2Aadvdiff(mesh,alfa,gamma,eta,potbeta); -Mass = BIM2Areaction(mesh,delta,zeta); -b = BIM2Arhs(mesh,f,g); - -A = AdvDiff + Mass; - - -## BOUNDARY CONDITIONS -Dlist = BIM2Cunknownsonside(mesh, [2 4]); ## LIST OF DIRICHLET NODES -Nlist = BIM2Cunknownsonside(mesh, [1 3]); ## LIST OF NEUMANN NODES -Nlist = setdiff(Nlist,Dlist); -Fn = zeros(length(Nlist),1); ## PRESCRIBED NEUMANN FLUXES -Ilist = setdiff(1:length(uin),union(Dlist,Nlist)); ## LIST OF INTERNAL NODES - -## LHS AND RHS PARTITION - -## LHS -## DIRICHLET ROW -Add = A(Dlist,Dlist); -Adn = A(Dlist,Nlist); -Adi = A(Dlist,Ilist); -## NEUMANN ROW -And = A(Nlist,Dlist); -Ann = A(Nlist,Nlist); -Ani = A(Nlist,Ilist); -## INTERNAL ROW -Aid = A(Ilist,Dlist); -Ain = A(Ilist,Nlist); -Aii = A(Ilist,Ilist); - -## RHS -bd = b(Dlist); -bn = b(Nlist); -bi = b(Ilist); - -## UNKNOWNS -ud = uin(Dlist); -un = uin(Nlist); -ui = uin(Ilist); - -## DISPLACEMENT SOLUTION -temp = [Ann Ani ; Ain Aii ] \ [ Fn+bn-And*ud ; bi-Aid*ud]; -un = temp(1:length(un)); -ui = temp(length(un)+1:end); - -## FLUXES THROUGH DIRICHLET SIDES -Fd = Add * ud + Adi * ui + Adn*un - bd; -u = zeros(nnodes,1); -u(Dlist) = ud; u(Nlist) = un; u(Ilist) = ui; - -## GRADIENT OF THE SOLUTION -[gx, gy] = BIM2Cpdegrad(mesh,u); - - -## FLUXES RECONSTRUCTION -[jxglob,jyglob] = BIM2Cglobalflux(mesh,u,alfa,gamma,eta,potbeta); - -## RESULTS VISUALIZATION -FPL2pdesurf(mesh,u); -FPL2pdesurf(mesh,u,"plot_field","gradient"); -FPL2quiver(mesh,jxglob,jyglob,"sample_density","1000"); -FPL2quiver(mesh,gx,gy,"sample_density","1000");
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/bim/examples/SQUARE/run_test.m Wed Jul 11 20:27:45 2007 +0000 @@ -0,0 +1,89 @@ +pkg load msh +pkg load fpl +pkg load bim + +## MESH GENERATION AND PROPERTIES +[mesh] = MSH2Mgmsh('quadrato',1,1.2); +[mesh] = BIM2Cmeshproperties(mesh); + +## INITIAL GUESS (IF NEEDED) +[xu, yu] = BIM2Cunknowncoord(mesh); +[nelem] = columns(mesh.t); +[nnodes] = columns(mesh.p); +[uin] = zeros(nnodes,1); + +## COEFFICIENTS FOR THE PROBLEM: +## -div ( \alpha \gamma ( \eta \nabla u - \beta u ) )+ \delta \zeta u = f g +epsilon = .01; +alfa = ones(nelem,1); +gamma = ones(nnodes,1); +eta = epsilon*ones(nnodes,1); +potbeta = xu/epsilon; +delta = zeros(nelem,1); +zeta = zeros(nnodes,1); +f = ones(nnodes,1); +g = ones(nelem,1); + +## ADVECTION AND REACTION MATRIX, RHS +AdvDiff = BIM2Aadvdiff(mesh,alfa,gamma,eta,potbeta); +Mass = BIM2Areaction(mesh,delta,zeta); +b = BIM2Arhs(mesh,f,g); + +A = AdvDiff + Mass; + + +## BOUNDARY CONDITIONS +Dlist = BIM2Cunknownsonside(mesh, [2 4]); ## LIST OF DIRICHLET NODES +Nlist = BIM2Cunknownsonside(mesh, [1 3]); ## LIST OF NEUMANN NODES +Nlist = setdiff(Nlist,Dlist); +Fn = zeros(length(Nlist),1); ## PRESCRIBED NEUMANN FLUXES +Ilist = setdiff(1:length(uin),union(Dlist,Nlist)); ## LIST OF INTERNAL NODES + +## LHS AND RHS PARTITION + +## LHS +## DIRICHLET ROW +Add = A(Dlist,Dlist); +Adn = A(Dlist,Nlist); +Adi = A(Dlist,Ilist); +## NEUMANN ROW +And = A(Nlist,Dlist); +Ann = A(Nlist,Nlist); +Ani = A(Nlist,Ilist); +## INTERNAL ROW +Aid = A(Ilist,Dlist); +Ain = A(Ilist,Nlist); +Aii = A(Ilist,Ilist); + +## RHS +bd = b(Dlist); +bn = b(Nlist); +bi = b(Ilist); + +## UNKNOWNS +ud = uin(Dlist); +un = uin(Nlist); +ui = uin(Ilist); + +## DISPLACEMENT SOLUTION +temp = [Ann Ani ; Ain Aii ] \ [ Fn+bn-And*ud ; bi-Aid*ud]; +un = temp(1:length(un)); +ui = temp(length(un)+1:end); + +## FLUXES THROUGH DIRICHLET SIDES +Fd = Add * ud + Adi * ui + Adn*un - bd; +u = zeros(nnodes,1); +u(Dlist) = ud; u(Nlist) = un; u(Ilist) = ui; + +## GRADIENT OF THE SOLUTION +[gx, gy] = BIM2Cpdegrad(mesh,u); + + +## FLUXES RECONSTRUCTION +[jxglob,jyglob] = BIM2Cglobalflux(mesh,u,alfa,gamma,eta,potbeta); + +## RESULTS VISUALIZATION +FPL2pdesurf(mesh,u); +FPL2pdesurf(mesh,u,"plot_field","gradient"); +FPL2quiver(mesh,jxglob,jyglob,"sample_density","1000"); +FPL2quiver(mesh,gx,gy,"sample_density","1000");
--- a/extra/bim/examples/SSHAPE/run.m Wed Jul 11 20:21:26 2007 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,85 +0,0 @@ -## MESH GENERATION AND PROPERTIES -[mesh] = MSH2Mgmsh('sshape',1,1.2); -[mesh] = BIM2Cmeshproperties(mesh); - -## INITIAL GUESS (IF NEEDED) -[xu, yu] = BIM2Cunknowncoord(mesh); -[nelem] = columns(mesh.t); -[nnodes] = columns(mesh.p); -[uin] = yu; - -## COEFFICIENTS FOR THE PROBLEM: -## -div ( \alpha \gamma ( \eta \nabla u - \beta u ) )+ \delta \zeta u = f g -epsilon = .01; -alfa = ones(nelem,1); -gamma = ones(nnodes,1); -eta = epsilon*ones(nnodes,1); -potbeta = 0; -delta = zeros(nelem,1); -zeta = zeros(nnodes,1); -f = zeros(nnodes,1); -g = zeros(nelem,1); - -## ADVECTION AND REACTION MATRIX, RHS -AdvDiff = BIM2Aadvdiff(mesh,alfa,gamma,eta,potbeta); -Mass = BIM2Areaction(mesh,delta,zeta); -b = BIM2Arhs(mesh,f,g); - -A = AdvDiff + Mass; - - -## BOUNDARY CONDITIONS -Dlist = BIM2Cunknownsonside(mesh, [5 28]);## LIST OF DIRICHLET NODES -Nlist = BIM2Cunknownsonside(mesh, setdiff([1:28],Dlist));## LIST OF NEUMANN NODES -Nlist = setdiff(Nlist,Dlist); -Fn = zeros(length(Nlist),1); ## PRESCRIBED NEUMANN FLUXES -Ilist = setdiff(1:length(uin),union(Dlist,Nlist)); ## LIST OF INTERNAL NODES - -## LHS AND RHS PARTITION - -## LHS -## DIRICHLET ROW -Add = A(Dlist,Dlist); -Adn = A(Dlist,Nlist); -Adi = A(Dlist,Ilist); -## NEUMANN ROW -And = A(Nlist,Dlist); -Ann = A(Nlist,Nlist); -Ani = A(Nlist,Ilist); -## INTERNAL ROW -Aid = A(Ilist,Dlist); -Ain = A(Ilist,Nlist); -Aii = A(Ilist,Ilist); - -## RHS -bd = b(Dlist); -bn = b(Nlist); -bi = b(Ilist); - -## UNKNOWNS -ud = uin(Dlist); -un = uin(Nlist); -ui = uin(Ilist); - -## DISPLACEMENT SOLUTION -temp = [Ann Ani ; Ain Aii ] \ [ Fn+bn-And*ud ; bi-Aid*ud]; -un = temp(1:length(un)); -ui = temp(length(un)+1:end); - -## FLUXES THROUGH DIRICHLET SIDES -Fd = Add * ud + Adi * ui + Adn*un - bd; -u = zeros(nnodes,1); -u(Dlist) = ud; u(Nlist) = un; u(Ilist) = ui; - -## GRADIENT OF THE SOLUTION -[gx, gy] = BIM2Cpdegrad(mesh,u); - - -## FLUXES RECONSTRUCTION -[jxglob,jyglob] = BIM2Cglobalflux(mesh,u,alfa,gamma,eta,potbeta); - -## RESULTS VISUALIZATION -FPL2pdesurf(mesh,u); -FPL2pdesurf(mesh,u,"plot_field","gradient"); -FPL2quiver(mesh,jxglob,jyglob,"sample_density","100"); -FPL2quiver(mesh,gx,gy,"sample_density","100");
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/bim/examples/SSHAPE/test.m Wed Jul 11 20:27:45 2007 +0000 @@ -0,0 +1,85 @@ +## MESH GENERATION AND PROPERTIES +[mesh] = MSH2Mgmsh('sshape',1,1.2); +[mesh] = BIM2Cmeshproperties(mesh); + +## INITIAL GUESS (IF NEEDED) +[xu, yu] = BIM2Cunknowncoord(mesh); +[nelem] = columns(mesh.t); +[nnodes] = columns(mesh.p); +[uin] = yu; + +## COEFFICIENTS FOR THE PROBLEM: +## -div ( \alpha \gamma ( \eta \nabla u - \beta u ) )+ \delta \zeta u = f g +epsilon = .01; +alfa = ones(nelem,1); +gamma = ones(nnodes,1); +eta = epsilon*ones(nnodes,1); +potbeta = 0; +delta = zeros(nelem,1); +zeta = zeros(nnodes,1); +f = zeros(nnodes,1); +g = zeros(nelem,1); + +## ADVECTION AND REACTION MATRIX, RHS +AdvDiff = BIM2Aadvdiff(mesh,alfa,gamma,eta,potbeta); +Mass = BIM2Areaction(mesh,delta,zeta); +b = BIM2Arhs(mesh,f,g); + +A = AdvDiff + Mass; + + +## BOUNDARY CONDITIONS +Dlist = BIM2Cunknownsonside(mesh, [5 28]);## LIST OF DIRICHLET NODES +Nlist = BIM2Cunknownsonside(mesh, setdiff([1:28],Dlist));## LIST OF NEUMANN NODES +Nlist = setdiff(Nlist,Dlist); +Fn = zeros(length(Nlist),1); ## PRESCRIBED NEUMANN FLUXES +Ilist = setdiff(1:length(uin),union(Dlist,Nlist)); ## LIST OF INTERNAL NODES + +## LHS AND RHS PARTITION + +## LHS +## DIRICHLET ROW +Add = A(Dlist,Dlist); +Adn = A(Dlist,Nlist); +Adi = A(Dlist,Ilist); +## NEUMANN ROW +And = A(Nlist,Dlist); +Ann = A(Nlist,Nlist); +Ani = A(Nlist,Ilist); +## INTERNAL ROW +Aid = A(Ilist,Dlist); +Ain = A(Ilist,Nlist); +Aii = A(Ilist,Ilist); + +## RHS +bd = b(Dlist); +bn = b(Nlist); +bi = b(Ilist); + +## UNKNOWNS +ud = uin(Dlist); +un = uin(Nlist); +ui = uin(Ilist); + +## DISPLACEMENT SOLUTION +temp = [Ann Ani ; Ain Aii ] \ [ Fn+bn-And*ud ; bi-Aid*ud]; +un = temp(1:length(un)); +ui = temp(length(un)+1:end); + +## FLUXES THROUGH DIRICHLET SIDES +Fd = Add * ud + Adi * ui + Adn*un - bd; +u = zeros(nnodes,1); +u(Dlist) = ud; u(Nlist) = un; u(Ilist) = ui; + +## GRADIENT OF THE SOLUTION +[gx, gy] = BIM2Cpdegrad(mesh,u); + + +## FLUXES RECONSTRUCTION +[jxglob,jyglob] = BIM2Cglobalflux(mesh,u,alfa,gamma,eta,potbeta); + +## RESULTS VISUALIZATION +FPL2pdesurf(mesh,u); +FPL2pdesurf(mesh,u,"plot_field","gradient"); +FPL2quiver(mesh,jxglob,jyglob,"sample_density","100"); +FPL2quiver(mesh,gx,gy,"sample_density","100");