pytanie zwi?zane z MES (FEM)

Jakbyś przytoczył to zdanie to byłoby to pomocne.

Nie bardzo rozumiem problem. Jeśli kłopot jest na poziomie elementu to dlaczego szukasz przykładu dla całego układu?

Z materiałów w sieci można polecić np.

dokąd nie opiszesz o co ci dokładnie chodzi to trudno zgadnąć na ile będą pomocne.

Bardzo ułatwiłoby odpowiedź gdybyś sprecyzował w czym jest problem. Kiedyś implementowałem element izoparametryczny, więc może będę w stanie Ci pomóc (choćby w znalezieniu literatury).

max

dziękuję za szybką "reakcję" - zaskoczenie :) Może powiem co mnie nie interesuje:

-fizyka problemu

-etap rozwiązywania układu równań

-etap global matrix assembly

-etap ustalania war. brzegowych

-postprocessing

-całkowanie numeryczne (teoria)

Interesuje mnie natomiast to co dzieje się pomiędzy otrzymaniem współrzędnych i parametrów fizycznych elementu, a uzyskaniem na wyjściu macierzy elementu i loading forces ( czyli w zasadzie całkowanie num. mnie bardzo interesuje, ale wydaje mi się, że posiadam już wystarczającą wiedzę i nie chcę poruszać tego tematu). Aktualnie piszę program (w zasadzie bibliotekę) służącą do powyższego: dostaje element po elemencie,a wyrzucam macierze lokalne i siły. Implementuje to w środowisku równoległym, nie na "zwykłym" procesorze, ale to akurat nie ma tu znaczenia.

To czego nie bardzo rozumiem: różnice pomiędzy iso-sub-superparametric

- o co chodzi?; czy mogę dla mapowania curvilinear wykorzystywać jako fcje kształtu fcje hierarchiczne (hierarchical shape functions) mimo, że mapuje przy wykorzystaniu Lagrangowskich wyższego rzędu?; sytuacje kiedy jakobian mapowania nie ma stałego znaku i jakie ograniczenia wprowadzić by tego uniknąć?; mapując fcjami Lagranga wprowadzam dodatkowe node'y kontrolujące geometrię - czy są to node'y,dla których będzie parametr w układzie równań - jeśli tak to co skoro użyję fcji kształtu wyższego rzędu Lagrange'a do mapowania, a niższego do interpolacji, albo wogóle użyję hierarchicznych gdzie node nie ma tak istotnego znaczenia?; czy mogę do mapowania użyć fcji hierarchicznych?; czy niezależnie od nawet najwymyślniejszego mapowania zawsze mogę całkować po prostym elemencie (prostokącie/ trójkącie prostokątnym i odpowiednikach w 3D)?;i, szczególnie dla mnie ważne, jaka będzie ogólna postać wyrażenia podcałkowego niezależnie od fizyki problemu - czy taka jak np. dla problemów elastic: B(transp.)DB?

heh, to w gruncie tyle :) nadal proszę o pomoc - może gdzieś, ktoś dobrze, krok-po-kroku, jak dziecku, to wytłumaczył? Od razu mówię, że w Zienkiewiczu, Reddym, Flahertym, Huttonie nie bardzo :/ (tzn. teoria niby jest, ale za mało szczegółów implementacyjnych jak dla mnie). Pozdrawiam i jeszcze raz dzięki!

[..]

Tak na szybko, z pamięci, bo strasznie mnie różne rzeczy gonią....

BTW. Piszesz to tylko dla siebie, czy upubliczniasz?

isop. Tak samo mapowana jest geometria jak rozwiązanie, subp. geometria jest mniej ogólna niż rozwiązanie, superp. geometria jest bardziej ogólna niż funkcja rozwiązania

Na przykładzie, wyobraź sobie trójkąt z sześcioma węzłami, naroża plus środki boków. Załóżmy, że poszukujemy rozkładu funkcji temperatury (bo najprostsza) T(x,y) wewnątrz elementu.

1. Izop. Wszystkie sześć węzłów definiują geometrię, czyli brzeg trójkąta może być parabolą. W każdym z sześciu węzłów jest niewiadoma temperatura węzłowa, czyli masz sześć niewiadomych a tak się składa, że tyleż współczynników ma wielomian kwadratowy ze względu na x i y, czyli T(x,y)=a+bx+cy+dx^2+exy+fy^2
2. Superp. Wszystkie sześć węzłów definiuje geometrię ale temperaturę węzłową masz tylko w wierzchołkach trójkąta, czyli funkcja temperatury jest T(x,y)=a+bx+cy. Węzły na środkach boków służą jedynie do definicji geometrii i nie posiadają stopni swobody.
3. Subparametryczne, odwrotnie. Tylko węzły w wierzchołkach definiują geometrię, czyli kształt elementu jest zawsze _trójkątny_ pomimo, że ma sześć węzłów, ale węzły na środkach boków muszą leżeć w połowie odległości między wierzchołkami. Natomiast T(x,y) jest opisana przez pełny wielomian kwadratowy, jak przy izoparametrycznym.

Teraz kluczowe: *po co tak?*.

Izop. nie wymaga wyjaśnień.

Subparam. pojawia się w sposób naturalny przy elementach klasy C1 (typowo płyty zginane), kiedy mamy wysokie wymagania odnośnie ciągłości rozwiązania ale geometrię już chcemy mieć klasy C0.

Superp. pojawia się w sposób naturalny w zagadnieniach sprzężonych, kiedy w elemencie szukamy więcej niż jednego pola i z matematyki wynika (warunek LBB), że jedno musi być funkcją niższego rzędu niż drugie (w Zienkiewiczu to pojawia się bodaj przy okazji nieściśliwości), wtedy z uwagi na jedną funkcję element jest izop. a właśnie z uwagi na tą drugą staje się superp.

niewiele pamiętam o hierarchicznych f.k., musiałbym zajrzeć tu i tam.

Jakobian odwzorowania nie może zmieniać znaku, to jest stosunek powierzchni w jednym układzie do powierzchni w drugim, zmiana znaku oznacza ,,zawinięcie się'' się powierzchni. Jeden wyjątek od tej reguły kiedy jakobian staje się zerem w punkcie, gdy z elementu czworokątnego poprzez złączenie węzłów robisz trójkąt i wtedy w tym podwójnym węźle jakobian ma prawo być zero.

Nieznane są jawne wyrażenia ograniczające aby zachować poprawną geometrię (element się nie zawinął), to się sprawdza przy okazji liczenia jakobianu, jeśli zmieni znak to exit(1) =:-). W praktyce problem jest daleko mniej dolegliwy niż by się mogło wydawać, fizyka narzuca o wiele ostrzejsze ograniczenia niż odwzorowanie. Np. mówi się, że element prostokątny nie powinien mieć stosunku boków większego niż 2:1. Nie oznacza to, że potem nie będzie liczył ale błędy będą rosły.

hierarchiczne j.w.

Tak, na tym właśnie polega urok xxx_parametrycznych sformułowań.

Ogólna postać nie istnieje, bardzo zależy od fizyki i zagadnienia, które rozwiązujemy. Z Twojego punktu widzenia, ważne jest jedynie, że pod całką stoi macierz funkcyjna, która może zależeć od dziwnych parametrów (np. miary spoistości gliny piaszczystej albo parametru Lodego) *oraz* od x i y a dokładnie od izopararametrycznych (r,s) czyli (ksi,eta). Ważne jest, że musisz wiedzieć jedynie jakiego rzędu wielomian daje macierz podcałkowa (np. ta BtDB) aby właściwie dobrać rząd całkowania. Jeśli piszesz bibliotekę, to miej świadomość, że czasami nie całkujemy zgodnie z rzędem ale inaczej (całkowanie zredukowane lub selektywne).

Spróbuję jeszcze coś poszukać o hierarchicznych, jeśli coś napisałem za skrótowo to daj znać.

max

Tymbardziej bardzo dziękuję Ci za odpowiedź!

To jest moja praca magisterska - mam nadzieję, że uczelnia nie zabroni mi jej upublicznić, bo oczywiście mam taki zamiar.

Świetne wyjaśniłeś te różnice - bardzo Ci dziękuję, teraz powoli to rozumiem. Małe pytanie - piszesz, że w przypadku subparam. węzły na bokach muszą leżeć w środku boków - dlaczego? I jeszcze, czy używa się krzywoliniowości wyższych rzędów?-spotkałeś się z czymś takim, czy paraboliczne boki na ogół wystarczają?

Akurat hierarchiczne, z tego co czytam, sprawują się o niebo lepiej niż Lagrange'owskie - ponoć macierz jest dużo lepiej "conditioned". Generalnie sprawa polega na tym, że węzłów masz tylko tyle ile wymaga tego opis geometrii niezależnie od rzędu fcji kształtu (nie tak jak przy Lagr. gdzie przy wyższych rzędach masz element upchany węzłami). Przy hierarchicznych np. rzędu 3, trójkąt nadal ma tylko 3 węzły a owe hierarchiczne fcje kształtu wyglądają następująco:

- 3 liniowe (takie same jak Lagrange'a dla wierzchołków trójkąta - zerują się na pozostałych wierzchołkach)

- 3 kwadratowe dla boków - po jednej na każdy bok (zerują się na pozostałych 2 bokach)

- 3 cubic dla boków - po jednej na każdy bok (zerują się na pozostałych 2 bokach)

- 1 cubic dla elementu (zeruje się na wszystkich bokach) Czyli dla trójkąta przy żadanej aprox. 3 stopnia mamy 10 fcji kształtu, a cały czas tylko 3 węzły. Oczywiście wartość szukanej nadal zdefiniowana jest jako: F(x,y) = SUMA(OD i=1 DO ILOSC_FCJI_KSZTALTU) Fcja_i(x,y)*param_i z tymże tylko 3 parametry (te dla fcji liniowych) związane są z "realnymi" węzłami (wierzchołkami trójkąta) pozostałe parametry nie ma ją już tak jasnego znaczenia.

No i tu trochę sprawa mi się komplikuje w przypadku mapowania geometrii wyższego rzędu, a że tych fcji hierarchicznych jak widzę używa się obecnie znacznie częściej niż Lagr. wobec tego próbuję to rozkminić...

ok, czyli sprawdzam znak podczas liczenia jakobianu + liczę na to, że siatka, którą dostaje na wejściu nie jest zbyt "pogięta"

wiesz, ja wyobrażam sobie to tak, że:

1. Na wejściu dostaję strukturę, element po elemecie, z informacjami: współrzędne el.+rząd mapowania+numeracja węzłów;rząd fcji kształtu;rodzaj fcji kształtu (Lagr/Hier);oraz pewna ustandaryzowana informacja o "fizyce" problemu (np. macierz fizyki problemu - coś w tym stylu, jeszcze tego dokładnie nie czuję, w każdym razie chodzi o to, żebym uniezależnił się od konkretnego problemu, dostał na wejściu po prostu dane opakowane w odpowiedniej formie, a czym one są to już mnie nie obchodzi)
2. Następnie w miarę nierozgałęzionym algorytmem obliczam dla każdego elementu local matrix i loading forces i wyrzucam je na wyjściu. (ewentualnie buduję od razu macierz globalną)

Jeszcze raz dzięki za pomoc. Niespodziewałem się, że ktoś zdecyduje mi się to tłumaczyć dlatego prosiłem w zasadzie tylko o linki - tymbardziej wielkie wielkie dzięki! Gdybyś znalazł coś o kwestii tych fcji hierarchicznych w kontekście elementów curvilinear to będę bardzo wdzięczny (sam oczywiście cały czas szukam) - rzecz jasna nie poganiam.

Pozdrawiam serdecznie! ps. gdybyś kiedyś pisał książkę o MES to daj znać - kupuję :)

Ja bym radził skupić się na ele. isop. lub subp. Elementy z hierarchicznymi funk. kształtu umożliwiają adaptacje typu p bez dodawania węzłów. Jak ja bym coś takiego robił to uogólnił bym implementacje do wykorzystania faktu że mes jest szczególnym przypadkiem metody podziału jedności (PUM). A wtedy bez problemu wyprowadzi się hierarchiczne funkcje kształtu i wiele innych. Tu musisz poprosić swojego promotora o literaturę i pomoc.

Poza tym, nie wiem czy jest sens stosować hierarchiczne funkcje kształtu do opisu geometrii. Proszę pamiętać że es musi być wypukły.

Max mnie poprawi ale z tego co pamiętam, to w przypadku ele. superp. czasami bywają problemy ze zbieżnością, lub jest ona marna, dlatego skreślił bym ten typ elementów na wstępie.

Bo wykorzystując idę podziału jedności można wzbogacić aproksymacje o wielomiany czybyszewa które są względem siebie ortogonalne.

Swoją drogą warto by sprawdzić jak taka aproksymacja zachowuje się dla materiałów prawie-nieściśliwych. Co z blokadą objętościową. Blokada ścinania pewnie nie ma.

Te elementy z hier. funkcj. kaszt. wyglądają fajnie, ale nie nadają się do problemów z nieliniowościami fizycznymi (np. plastyczność) bo nie zapewniają odpowiedniego rzędu zbieżności. W tym wypadku najlepsze są ,,proste'' elementy i adaptacja typu h.

Albo dla elementów czterobocznych ośmiowęzłowych (serendipity element) stosuje się kwadraturę 2x2 mimo że z rzędu wielomiany wynika że powinna być 3x3.

Ja bym radził zapomnieć albo o nieliniowościach fizycznych albo/i o przetwarzaniu równoległym.

Coś o elementach hierarchicznych znajdziesz na tej stronie,

pewno coś więcej będzie w czasopismach naukowych, poproś promotora o pomoc.

Ja popieram ten pomysł :)

likask

Przepraszam za zwłokę ale akurat mam tzw. ,,urwanie głowy''.

Tekstu pracy nie może zabraniać upublicznienia, z kodem bywa różnie.

[..]

Formalnie to tylko boki muszą być proste. No ze względu na jednoznaczność węzłowych (fizycznych) stopni swobody węzły pośrednie nie mogą pokrywać się z narożnymi. Czyli, znowu czysto formalnie, węzły pośrednie mogłyby leżeć gdziekolwiek na długości boku. Z uwagi na uwarunkowanie i z czysto praktycznych względów najlepiej je umieścić w środku. Ale gdyby były jakieś ważne przyczyny mogłyby leżec np. w 1/3.

Są takie elementy, przynajmniej opisane. Nie widziałem aby ktoś je stosował. W praktyce paraboliczne wystarczają, co można uzasadnić następująco. Jeśli brzeg jest krzywoliniowy (pomyślmy o funkcji sin(x)) a bok elementu prostoliniowy to trzeba bardzo wielu elementów aby łamana dobrze przybliżała brzeg. Jeśli ten sam brzeg przybliżymy parabolami, to trzeba stosunkowo niewiele elementów aby przybliżenie ,,optycznie'' było nieodróżnialne od rzeczywistego kształtu brzegu.

znalazłem jakieś materiały o hier. ale nie mam czasu aby poczytać,

Coś tu nie gra, gdzie st. sw.?

Dokładnie, oprócz tego możesz liczyć na to że na etapie prezentacji graficznej użytkownik zobaczy ,,zawiniętą'' siatkę albo że generator siaki nie zezwoli na coś takiego.

Dokładnie tak, tylko jest jedna subtelność. właśnie owo coś o fizyce problemu. Tak naprawdę musisz na tym etapie przekazać porcję informacji (nie wiesz jaką!) i wskaźnik do funkcji, która to będzie przetwarzać. Trick polega na tym, że to jest tzw. poziom punktu czyli w tym miejscu musisz wygenerować macierz D (dla liniowej statyki) lub D_styczne plus siły wewnętrzne (dla nieliniowej statyki). Te wielkości mogą zależeć od nieznanych (dla ciebie w chwili tworzenia programu) parametrów i funkcji, którą dopiero napisze przyszły użytkownik.

Odradzam to ewentualnie. Macierz elementu jest mała (powiedzmy do

100x100) i reprezentuje się ja jako pełną tablicę. Macierz układu przy 100 tys. stopni swobody zajęłaby 8*10^5*10^5 B = 80*10^9 B ~ 80 GB RAMu. Na poziomie układu _musimy_ używać rzadkich reprezentacji macierzy i dlatego nie należy mieszać tych dwu rzeczy.

max

likask wrote: [..]

Ale, o ile dobrze rozumiem, to ma być biblioteka, więc nie może czegoś nie być bo nie. Ktoś może potrzebować tego akurat w jakimś specyficznym końtekście. To, że ogólnie albo w innym szczególnym przypadku coś może się nie zachowywać dobrze niewiele znaczy.

Widzę, że fascynujesz się Partition of Unity, ale to jest praca magisterska a nie doktorska. Na to jest pół roku, nie da się zrobić jeszcze PUM w tym czasie =:-)

Obiecuję, że za dwa tygodnie doczytam hierarchiczne i będę w stanie dyskutować =:-)

Superp. nie stosuje się samodzielnie ale dla pól sprzężonych więc nie bardzo mamy wyjście.

Nie to jest cecha elementu *i* materiału. Tu mówimy o szkielecie do implementacji elementów. Więc takie badanie to inny etap.

j.w., to ma wiedzieć ten który się będzie bawił plastycznością. IMHO autor wątku pisze getFEMa/Dolphina/.../

To jest słuszne! Do tego, zapewne twoja biblioteka ma dostęp elektroniczny do czasopism, nic tylko szukać i czytać. Szczególnie nowe rzeczy łatwo znaleźć.

Słuchajcie, nie podpuszczajcie mnie bo jeszcze napiszę i będziecie musieli to czytać =:-)

max

Nie warto skupiać się na mało praktycznych metodach. Biblioteka nie jest dobra gdy zawiera wszystkie możliwe metody (to jest niemożliwe). Jest dobra, gdy łatwo się jej używa, ma dobrze przemyślaną strukturę, jest szybka itd. I najważniejsze można ją rozszerzyć o nowe podejścia.

PUM w FEM ( nazywanym często XFEM ) już się nie zajmuje. Mimo, że bardzo modne teraz. Fajne to jest, ale ja wole coś innego.

Znalazłem coś znaczenie fajniejszego, odkryłem że Trefftz bije na głowę XFEM. Po prostu dzieli się element, żadnych trików z całkowaniem po obszarze es, bo go nie ma :). Brak obaw o jakość siatki. A jak wiesz tam nie ma PUM wcale. Dodatkowo, udało się sformułowanie HTS-Discontinus rozszerzyć do nieliniowości geometrycznych. Artykuł już napisany, ostatnie poprawki w toku.

Tam nie ma filozofii. A es musi być wypukły, jak dobrze pamiętam, wystarczy sprawdzić co się dzieje z jakobianem jak nie jest.

Tu nie rozumiem, co ma to wspólnego z geometrią. Subparam. tak ale nie superparamet.

Racja. Ja sobie pomyślałem że to dobry temat na prace magisterską.

Piesze to bo nie wiem do czego ma być ta biblioteka. Warto być świadomym ograniczeń.

Tutaj mam wątpliwości czy temat jest dobrze przemyślany.

Witam, Jeżeli chodzi o książke to podpisuję się dwiema rękami i dwoma nogami. Też z chęcią kupię.

Pozdrawiam Szymon Hernik