Pokročilé programovací techniky v Matlabu T....

Post on 24-Feb-2021

5 views 0 download

transcript

Pokro čilé programovací techniky v Matlabu

T. Kozubek

MI21 14.11. 2010

Katedra aplikované matematikyFakulta elektrotechniky a informatikyVŠB-Technická univerzita Ostrava

Obsah1. Matlab

1. Funkce, řízení toku2. Homogenní pole (číselné matice, řetězce)3. Struktury4. Pole buněk5. OOP

MI21 24.11. 2010

6. GUI7. Import, export, MEX

2. Fourierova řada (spektrální analýza)

3. Diskrétní Fourierova transformace (spektrální analýza)

4. Konvoluce (vyhlazování signálu)

5. Laplaceova transformace (řešení DR)

6. Úlohy mechaniky – MatSol, paralelní programování

Matlab – funkcefunction [mean_x,std_x,var_x]=statistics(x)% STATISTICS Spocte zakladni statistiky vstupniho v ektoru x% UZITI: [mean,std,var]=statistics(x)

% Test na pocet vstupnich argumentuif nargin~=1; error( 'Spatny pocet vstupnich argumentu!' ); end;

% Test na korektnost vstupnich argumentu (x je cise lny vektor)if ~isvector(x); error( 'Vstupni argument x neni vektor!' ); end;

MI21 34.11. 2010

if ~isvector(x); error( 'Vstupni argument x neni vektor!' ); end;if ~isnumeric(x); error( 'Vstupni argument x neni ciselny vektor!' ); end;

% Spocti zakladni statistikymean_x=mean(x); std_x=std(x); var_x=var(x);

end

% Podfunkce MEAN (upraveny vypocet stredni hodnoty)function f=mean(x)% Vypocet stredni hodnotyf=sum(x)/(length(x)+1);end

Matlab – homogenní poleVytvo ření matice

A=[1 2 32 3 45 6 7];

B=[0,0,0; 2,3,4; 5,6,7];

Přístup k prvk ům, rušení

A(1:2,3)=[1;3]

A(A<7)=0

A(2,:)=[]

MI21 44.11. 2010

B=[0,0,0; 2,3,4; 5,6,7];

C=[A,B; B,A]; S=[‘abc’;’def’]

A(2,:)=[]

Řídké maticee=ones(5,1); A = spdiags([-e 2*e -e], -1:1, 5, 5)B=sprand(5,4,0.1)I=[1 1 2 3]; J=[1 3 2 4]; V=[1 1.5 2 3.7]; S=sparse(I,J,V,5,6)

spy(S,'o',8); imagesc(S); colorbar;

Matlab - struktury• Struktura je heterogenní pole, jehož prvky jsou libovolné objekty –

instance tříd (řetězce, číselné pole, buněčné pole, struktury, atd). • Přístup k položkám je prostřednictvím názvů položek užitých v

definici struktury.>>student=struct('jmeno','Jan Madlo', 'predmet',‘LAM','znamky',[1 2 1]);

MI21 54.11. 2010

• Paměťové požadavky : není třeba souvislý blok paměti pro celou strukturu, ale pouze pro položky struktury.

• Využití : snížení počtu vstupních a výstupních argumentů funkce, zvýšení čitelnosti kódu.

• Přístup k položkám>>student.jmeno>>student.jmeno=‘Miro Vilik'

Matlab – pole buněk• Pole bun ěk je heterogenní pole, jehož prvky jsou libovolné objekty

– instance tříd (řetězce, číselné pole, buněčné pole, struktury, atd). • Přístup k položkám je prostřednictvím indexů buněk.>>c={‘Hello’,rand(3),student};

• Paměťové požadavky : není třeba souvislý blok paměti pro celé

MI21 64.11. 2010

• Paměťové požadavky : není třeba souvislý blok paměti pro celé pole, ale pouze pro jeho položky.

• Využití : uložení řetězců různé velikosti, uložení matic různého řádu u metody rozložení oblastí, nahrazení vstupních a výstupních argumentů funkce, zvýšení čitelnosti kódu.

• Přístup k položkám>>c{1}, c{1}(2:end), c{3}.jmeno>>c{1}=‘Ahoj’, c{3}.jmeno=‘Ivan Leden’

Matlab – OOP

MATLAB má zabudovánu kompletní podporu OO přístupu

• Jak systém roste, je čím dál složitější sledovat předávání dat mezi funkcemi.

MI21 74.11. 2010

• Opakování kódu copy&paste –redundance – nekonzistence.

• Deklarace třídních proměnných – předem se ví, jaké má proměnné.

• Srozumitelnost a přirozenost návrhu.

Matlab – OOPclassdef MyClass < MySuperClass

propertiesProperty1Property2 = sin(pi/12); % default value

endproperties (SetAccess = private, GetAccess = private)

StressStrain

end

MI21 84.11. 2010

endmethods

function value = get.Property1(obj)%optionally implement get value

endfunction obj = set.Property1(obj,value)

%optionally implement set valueend

endend

Dědičnost, přetěžování operátorů, konstruktory, proměnné s omezeným přístupem.

Matlab – GUI

1. Pomocí LAYOUT editoru (.fig, .m)

2. Pomocí Matlabovských příkazů

MI21 94.11. 2010

(figure, uicontrol, uibuttongroup, guidata, guihandles, uimenu, …)

Matlab – import, export, MEX

1. Podpora práce s textovými i binárními soubory2. Načtení obrázků, zvuků, videa, MS Excell, …

Načtení a zobrazení obrázkurgb = imread('hokej.jpg');

Načtení a p řehrání zvukového souboru[y,fs,nbits]=wavread('babycry.wav');

MI21 104.11. 2010

rgb = imread('hokej.jpg');image(rgb);

[y,fs,nbits]=wavread('babycry.wav');sound(y,fs); plot(y);

MEX – přilinkování k Matlabu samostatných funkcí i celých knihoven v C, C++, Fortranu, Javě

( ) CaCzzazf ∈∈= ,,

Lineární zobrazení

MI21 114.11. 2010

( ) ))arg()(arg()arg()arg( zaiziai ezaezeazazf +===

Fourierova řada (spektrální analýza){ }

( ) RtectfTLtf

TT

,TLTe

ntin

Znnti

∈=∈

=

∑∞+

∞−

,),,0()(

perioda...rychlost,úhlová...2

)0(fcísystémortonorm..../

2

2

ω

ω

πω

( ) ,,1

ZndtetfcT

nti ∈∫= − ω

MI21 124.11. 2010

( ) ,,1

0

ZndtetfT

c ntin ∈∫= − ω

• Dvoustranné spektrum : {cn,ϕ

n}

Fourierova řada

( ) ,,1

0

ZndtetfT

cT

ntin ∈∫= − ω( ) ,, Rtectf nti

n ∈= ∑+∞

∞−

ω

( ) 1,...,0,11 1

0

21

0

2

−==∆∆

= ∑∑−

=

−−

=

∆∆

−Nnef

Ntetf

tNc

N

k

nkN

i

k

N

k

tnktN

i

kn

ππ

Diskrétní Fourierova transformace (spektr. analýza)

MI21 134.11. 2010

00∆ == NtN kk

1,...,0,)(1

0

21

0

2

−=== ∑∑−

=

=

∆∆ Nkecectf

N

n

nkN

i

n

N

n

tnktN

i

nk

ππ

1,...,0,1

,1,...,0,1

0

21

0

2

−==−== ∑∑−

=

−−

=

NnefN

cNkecfN

k

nkN

i

kn

N

n

nkN

i

nk

ππ

DFT

Konvoluce posloupností (vyhlazování signálu)

{ } { }∞=

∞= 00 , llkk ba

{ } { } { } { }lknMN

nn bacc *,1

0 =−+=

{ } { } { } { } ∑=

−∞

= =∗=n

kkknnnnlkn baccbac

00,,

{ } { } 1

0

1

0 , −=

−=

N

llM

kk ba

MI21 144.11. 2010

==

−−

1

2

2

1

0

1

21

01

01

0

2

1

0

000

00

00

000

0000

00000

0

0

0

N

N

M

MM

M

b

b

b

b

b

a

aa

aa

aa

a

a

a

a

Abc

⋯⋯

⋯⋯

⋮⋱⋮⋱⋮⋮

⋯⋯

( ) ( ) ( ) ,tfyayayaya nn

nn =+′+++ −

− 011

1 ⋯

( ) ( ) ( ) ( )10 0 0 1 0 1

nny t c , y t c , ,y t c .−−′= = =⋯

( ) ( )∞

Laplaceova transformace ( řešení DR)

MI21 154.11. 2010

L(DR)DR Y(s) y(t)

( ) ( )∫∞

−==0

,.}{ dtetfsFfLT st

ANSA

modelování geometrie

ImporterSTL, IGIS,

Specifikace problému a generování sít ě.

(typ úlohy, materiálové vlastnosti, počáteční a okrajové podmínky, …)

Exporter

MI21 164.11. 2010

CAD System Integration

•Autodesk Inventor

•CATIA

•Pro/ENGINEER

•Solid Edge

•SolidWorks

•Unigraphics

STL, IGIS, VRML, … CDB format

T. Kozubek, T. Brzobohatý, A. MarkopoulosZ. Dostál, V. Vondrák, R. Ku čera, M. Sadowská

knihovna se škálovatelnými algoritmy založena na metodách rozložení oblasti FETI a BETI

Fproblem

FETI

1.podoblasti jsou volné

nebo uchycené

(předem není znám

defekt matic K)

F F

Metody rozložení oblasti

MI21 174.11. 2010

FETI-

DP

TFETI

2.

3.

defekt matic K)

FETI-DP (částečné

rozdělení, regulární

matice tuhosti)

všechny podoblasti

defekt matic K je znám

F

F

F

F

Fproblem

FETI

1.subdomains are fixed or free but

with different defectsF F

BETIIntroduction to Total FETI/BETI

MI21 184.11. 2010

FETI-DP

TBETI

2.

3.

FETI-DP (partial splitting,

nonsingular)

všechny podoblasti

jsou volné se shodnou

dimenzí jádra

FF

F

fuKuu TT −21min I I

G G

B B

B u c

B u c

B u c

≤==

u pole posuvůK matice tuhostiB matice omezeníc vektor omezení

funkcionál energie

Total FETI – primární proměnné

MI21 194.11. 2010

FF

obstacle

c vektor omezení

FETITotalDirichletovy okrajové podmínky jsou také

vynuceny pomoci LM, B

Bu= c

B

Vše a mnohem více m ůžete najít Table of Contents - Preface.

Part I. Background1. Linear Algebra.- 2. Optimization.

Part II. Algorithms3. CG for Unconstrained Minimization4. Equality Constrained Minimization 5. Bound Constrained Minimization 6. Bound

MI21 204.11. 2010

Bound Constrained Minimization 6. Boundand Equality Constrained Minimization

Part III. Applications to VariationalInequalities7. Solution of a Coercive VariationalInequality by FETI-DP method8. Solution to a Semicoercive VariationalInequality by TFETI Method.- References.-Index.

Seminář výpočetní mechaniky, 1-3. prosince 2010, VŠB-TU & Ústav geoniky AV ČR, v.v.i.

Škálovatelnost: kostka nad tuhou překážkou

MI21 214.11. 2010

MI21 224.11. 2010

23

OOK ⋮

1

blokově diagonální struktura

matice tuhosti =>

vhodná pro paralelní

implementaci

1. ingredience: paralelní implementace

MI21 234.11. 2010

=

nKOO

OKO

OOK

K

⋯⋱⋯⋯

2

Mohou být řešeny koercivní

a semikoercivní úlohy.

• Paralelně smyčky for– mnoho cyklů

Paralelní programování

Typické případyCluster HP model BLc7000 (c-class).

Konfigurace výpo četních uzl ů:2x dual core CPU AMD Opteron 2210 HE 8GB ECC DDR2 667MHz RAM

MI21 244.11. 2010

– dlouhé iterační cykly

• Offloading Work• Objemné datové sady

IT4Innovations, Centrum excelence, Ostrava, Czech Republic, cluster with more than 30,000 processors.http://www.it4i.eu

Matlab Distributed Computing Server

ClientParallel

Scheduler

WorkerDistributed

Computing Server

Worker

MI21 254.11. 2010

ParallelComputing

Toolbox

Scheduleror

Job manager

WorkerDistributed

Computing Server

WorkerDistributed

Computing Server

• Local job manager

• MathWorks job manager

• Third-party schedulers

(Windows CCS, PBS Pro, LSF)

Job

Job

Job

Job

Job

Job

Scheduler

Pending

Queued RunningWorker

Worker

Life cycle of a job

MI21 264.11. 2010

Client

Job

Job

Job

Job

Job

Job

Job

Job

Job

JobgetAllOutputArguments

submitFinished

Worker

Worker

Worker

createJob

Distributed computations1. Find scheduler or

job manager2. Create a job3. Create tasks and

associate them withthe job

>>sched=findResource(‘scheduler‘,

’type’,’local’);

>>job=createJob(sched);

>>createTask(job,@rand,1,{3,3});

>>createTask (job,@eye,1,{4});

MI21 274.11. 2010

the job

4. Send job to the front

5. Wait until job finishes

6. Gather results7. Destroy job

>>createTask (job,@eye,1,{4});

>>createTask (job,@ones,1,{{4},{3}});

>>submit(job);

>>waitForState(job);

>>results=getAllOutputArguments(job )

>>destroy(job);

Parallel computations>>sched=findResource(‘scheduler‘,…’type’,’local’);

>>pjob=createParallelJob(sched);set(pjob,‘MaximumNumberOfWorkers,4);set(pjob,‘MinimumNumberOfWorkers,4);

1. Find scheduler or job manager

2. Create a parallel job and sets the numer of tasks

3. Create parallel tasks

MI21 284.11. 2010

>>createTask(pjob,@MyParTask,1,{});

>>submit(pjob);

>>waitForState(pjob);

>>results=getAllOutputArguments(pjob)

>>destroy(pjob);

MPT functions: numlabs, labindex, labBarrierlabBroadcast, labProbe, labReceive, labSend,

labSendReceive

tasks4. Send job to the

front5. Wait until job

finishes6. Gather results7. Destroy job

Parallel computation scenario

MasterSlave 1

(subdom. 1)Slave 2

(subdom. 2)Slave 3

(subdom. 3)

Send FEM

Assemble K,K+,R,f

Build

Send FEM models to workers and parallelize assembling of K, K+, R, fv, stresses, searching

contact pairs, multiplication procedures.

MI21 294.11. 2010

Build

contact

conditions Gather f

Multiply by K+

Gather resultIterate

Compute stress

Gather stress

X

Důlní pr ůmysl: sv ěrný spoj d ůlní výztuže

MI21 304.11. 2010

tlak vyvolaný okolním materiálem

ocelová výztužpřekryv

R1

R2

t2

t1

Coulombovské t ření

f = 0.1

E = 2.1e5 MPa,

µ = 0.3

107m

m

BB

C

F

AA

B

C

A

Důlní pr ůmysl: sv ěrný spoj d ůlní výztuže

MI21 314.11. 2010

rovina symetrie

tA = tD = 10 mm

tB = tC = 40 mm

107m

m

F

DD

počet:

binárních proměnných 65562

duálních proměnných 3112

D

Důlní pr ůmysl: sv ěrný spoj d ůlní výztuže

MI21 324.11. 2010

8 těles7 volných

redukované napětí HMHPosunutí

2D

Důlní pr ůmysl: sv ěrný spoj d ůlní výztuže

MI21 334.11. 2010

3D

1NODAL SOLUTION

bez korekcí kontaktních směru 3 korekce

σN = 985.4 MPa σN = 1020.6 MPa

MI21 344.11. 2010

MN

MX

PR03.db FKOP= -0.04 f=0.0 PLANE182

.012349126.416

252.821379.225

505.629632.033

758.437948.043

DEC 9 200811:06:26

PLOT NO. 8

STEP=2SUB =34TIME=.5CONTSTOT (AVG)DMX =.100117SMN =.012349SMX =948.043

ANSYS (paralel)dynamický výpočet: 27 hodin

ANSYS (paralel) statický výpočet: 1 hodina

MATSOL statická úloha, 3 korekce:

sekvenční verze 15 min.

paralelní 7 min.

TBETI TFETI1.548

1.611

0.774

0.387Celková posunutí [mm]

Srovnání TBETI a TFETI

MI21 354.11. 2010

TBETI TFETI

0.387

0.017

1.548

1.611

0.774

0.387

0.017

Celková posunutí [mm]

Čas výpočtu 2.41 h.

Celkový čas 2.55 h.

násobení matice

Čas výpočtu 2.54 h.

Celkový čas 3.33 h.

násobení matice

250 podoblastí f = 0.1

Srovnání TBETI a TFETI

MI21 364.11. 2010

násobení matice vektor

2 126

primární proměnné

1 592 853

duální proměnné 261 553

max(abs(u)) 1.303 mm

max(totalU) 1.548 mm

násobení matice vektor

1 882

primární proměnné

713 751

duální proměnné 261 553

max(abs(u)) 1.378 mm

max(totalU) 1.616 mm

FETI BETI

vnější kroužek

klec

vnitřní kroužek

hřídel

zatížení

Automobilový pr ůmysl: kuli čkové ložisko

MI21 374.11. 2010

uchycení

3

Válečkové ložisko v ětrné elektrárny

MI21 384.11. 2010

Statisticspočet těles 73

počet podoblastí 700

primární prom ěnné 2,73 M

duální prom ěnné 459,8 k

počet iterací 4270

Automobilový pr ůmysl: pneumatika 3.9 mil stupňů volnosti,4000 iterací

MI21 394.11. 2010

celková posunutí redukované napětí HMH

Světlomet – Visteon (Autopal)

12 těles

MI21 404.11. 2010 VŠB – TU Ostrava

Světlomet – Visteon

64 Oblastí

MI21 414.11. 2010 VŠB – TU Ostrava

Světlomet – Visteon (Test solution-Total Displacement)

MI21 424.11. 2010 VŠB – TU Ostrava

Sestava složená z 29 těles

Dveře auta

MI21 434.11. 2010 VŠB – TU Ostrava

Dveře auta

MI21 444.11. 2010 VŠB – TU Ostrava

Dveře auta

128 Oblastí

MI21 454.11. 2010 VŠB – TU Ostrava

Dveře auta (Test solution-Total Displacement)

MI21 464.11. 2010 VŠB – TU Ostrava

Auto – Sestava obsahuje okolo 700 těles

MI21 474.11. 2010 VŠB – TU Ostrava

3000 subdomains

MI21 484.11. 2010 VŠB – TU Ostrava

700 bodies

MI21 494.11. 2010 VŠB – TU Ostrava

3000 subdomains

Rozšíření na d ynamic ké úlohy

MI21 504.11. 2010

Děkuji za pozornost

MI21 514.11. 2010