Finite Elemente Methode



  • Hallo, ich versuche gerade mir selbst die FEM näher zu bringen.
    Aktuell versuche ich mich an einem Referenzdreieck T, das aus den Punkten (0,0),(1,0) und (0,1) besteht.

    Die nodalen Basisfunktionen sind dann:

    \phi_1 = 1-x-y \\ \phi_2 = x \\ \phi_3 = y \\

    Nun sollte über dieses Dreieck eine Matrix A(i,j) bestimmen können:
    Aij=_Tϕ_iϕjdTA_{ij} = \int\_T \nabla \phi\_i \nabla \phi_j dT

    als Ergebis erhalte ich hier:

    12(211110101)\dfrac{1}{2} \begin{pmatrix} 2 & -1 & -1 \\ -1 & 1 & 0 \\ -1 & 0 & 1 \end{pmatrix}

    Ich weiß nicht, was ich genau erwartet habe, aber irgendwas in Richtung 2d5p-Stencil

    ( 1 141 1 )\begin{pmatrix} ~ & -1 & ~ \\ -1 & 4 & -1 \\ ~ & -1 & ~ \end{pmatrix}

    oder einem anderen mir bekannten Differenzenschema wäre schön gewesen.

    Habe ich schon irgendwo einen Fehler drin, oder ist ein Dreieck einfach zu wenig?



  • Das Ergebnis kann schon stimmen.
    Wenn du den 5-Punkte-Stern wieder herausbekommen willst, schau dir ein regelmäßiges Schachbrettmuster aus diesem und dem gespiegelten Dreieck an
    und summiere einfach deine Ergebnisse für alle beteiligten Dreiecke auf.
    Für A_ii hast du z.B. zwei Dreiecke mit rechten Winkel am Knoten und 4 mit 45 Grad. Macht laut deinen Berechnungen 2*1 + 4*0.5 = 4.
    Sieht schonmal vielversprechend aus.



  • Ok. Mein Ergebnis scheint zu stimmen. Habe es auchin einem Buch so gefunden. Das Beispiel mit der regelmäßigen Zerlegung war auch anschaulich.

    Leider habe ich das Prinzip noch nicht ganz verstanden.
    Ich wollte das ganz mal auf eine Zerlegung von Matlab anwenden, wie si zB hier zu finden ist:
    https://de.mathworks.com/help/pde/examples/poisson-s-equation-on-a-unit-disk.html

    Nun sind in der Zerlegung keine Dreiecke mehr mit rechten Winkeln.
    Wie ordne ich die Knoten einer solch komplexen Zerlegung den Knoten meines Referenz-Dreiecks zu? Das muss ja irgendwie einheitlich geschehen, damit ich nicht einfach alle 1 Elemente benachbarter Dreiecke zusammenlege.



  • Du transformierst dein Dreieck affin auf das Referenzdreieck.
    Siehe z.B. hier: https://www.uni-muenster.de/AMM/num/Vorlesungen/WissenschaftlichesRechnen_WS1213/Dateien/Fem-intro.pdf
    Sollte eigentlich auch in jedem Buch beschrieben werden, z.B. im Braess: http://www.springer.com/de/book/9783642347962
    Ansonsten such mal nach geTeXten FEM-Skripten (um den meisten Ingenieurs-Schrott gleich rauszufiltern ;))



  • ALso ich versuche nun ein Quadrat zu berechnen, das in 8 reguläre Dreiecke unterteilt ist:

    ___
    |\|\|
    |\|\|
    

    Ich nummeriere meine 9 Knoten beginnend von unten links:

    6 7 8
    3 4 5
    0 1 2

    Und bilde damit die folgenden Dreiecke, die ich gleichzeitig auch als Matrix L interpretiere:

    0,1,3
    1,2,4
    1,4,3
    2,5,4
    3,4,6
    4,5,7
    4,7,6
    5,8,7

    Für jedes der Dreiecke habe ich die Standard-Matrix A aus meinem ersten Beitrag.

    Nun berechne ich die globale Matrix folgendermaßen:

    Für jedes Dreieck e addiere ich den Wert der lokalen Matrix KijeK^e_{ij} auf das gloable Element KLe,i,Le,jK_{L_{e,i}, L_{e,j}} .

    In Python mal:

    A = [[1, -0.5, -0.5], [-0.5, 0.5, 0], [-0.5, 0, 0.5]]
    
    L = [[0, 1, 3],
     [1, 2, 4],
     [1, 4, 3],
     [2, 5, 4],
     [3, 4, 6],
     [4, 5, 7],
     [4, 7, 6],
     [5, 8, 7]]
    
    K = [[0 for col in range(9)] for row in range(9)]
    
    for e in range(len(L)):
           for i in range(3):
               for j in range(3):
                    ii = L[e][i]
                    jj = L[e][j]
                    K[ii][jj] += A[i][j]
    

    Leider erhalte ich als Ergebnis:

    [[  1,  -0.5,    0, -0.5,    0,    0,    0,    0,    0],
     [-0.5,  2.5, -0.5, -0.5, -1.0,    0,    0,    0,    0],
     [   0, -0.5,  1.5,    0, -0.5, -0.5,    0,    0,    0],
     [-0.5, -0.5,    0,  2.0, -0.5,    0, -0.5,    0,    0],
     [   0, -1.0, -0.5, -0.5,  4.0, -0.5, -0.5, -1.0,    0],
     [   0,    0, -0.5,    0, -0.5,  2.0,    0, -0.5, -0.5],
     [   0,    0,    0, -0.5, -0.5,    0,  1.0,    0,    0],
     [   0   , 0,    0,    0, -1.0, -0.5,    0,  1.5,    0],
     [   0,    0,    0,    0,    0, -0.5,    0,    0,  0.5]]
    

    was in keinster Weise meinen Erwartungen entspricht.

    Wo liegt hier mein Fehler?



  • Mit deinen Dreiecken stimmt was nicht, der rechte Winkel muss beim ersten Knoten liegen, also z.B. 5,4,2 statt 2,5,4.
    Wenn du das korrigiert hast, solltest du in den Einträgen K_i,4 den 5-Punkte-Stern wiederfinden.



  • Wieso "muss" der rechte Winkel beim ersten Knoten liegen?

    Oder allgemeiner: Was muss ich beim Mapping beachten? Anscheinend macht die Zuordnung lokale Nummerieung zu globaler ja einen entscheidenden Unterschied.
    Wenn ich dann eine allgemeine Zerlegung eines Gebietes vorliegen habe, dann habe ich ja meist auch keine rechten Winkel mehr. Irgendwie muss ich dann ja trotzdem wissen, wie ich nummerieren darf, oder?



  • Weil bei dem Integral i.A. was anderes herauskommt, wenn du über ein anderes Dreieck integrierst.
    Es kommt nur das gleiche heraus, wenn dein anderes Dreieck eine Drehung + Verschiebung des Referenzdreiecks ist, und bei dem ist der rechte Winkel nunmal beim ersten Knoten.
    Wie man das Integral für allgemeine Dreiecke ausrechnen kann, findest du in passenden FEM Skripten.
    Siehe meinen Link 1.34 und http://numerik.mathematik.uni-mainz.de/~juengel/scripts/femscript.pdf Seite 46
    (da ist auch im Folgenden genau dieses Beispiel beschrieben)



  • Hallo nochmal.

    Ich glaube, dass ich nun eine Matrix erfolgreich aufstellen konnte (sie sieht nach 5-PunktStern aus), habe aber nun Probleme mit der rechten Seite.

    Lösen möchte ich das Problem

    Δu=f-\Delta u = f auf dem Rechteck von [-1,-1] bis [1,1].

    Auf dem Rand soll u = 0 gelten.

    f=x2+y22f = x^2 +y^2 -2

    Das Gebiet ist in Courant-Dreecke (gleichschenklig mit rechtem Winkel) zerlegt.

    Wenn ich die FEM richtig verstehe, muss ich für die rechte Seite meines Gleichungssystems den Vektor b berechnen mit:

    b_i=_Ωfϕidxdyb\_i = \displaystyle\int\_\Omega f \, \phi_i \;dx \, dy

    Ich weiß nun leider nicht, wie ich das Produkt von f und den Testfunktionen integrieren soll.

    Ich würde mich über eine Art "Musterlösung" für die rechte Seite freuen und hänge dafür mal die Zerlegung meines Gebietes im vtk-Format an:

    # vtk DataFile Version 3.0
    
    ASCII
    DATASET POLYDATA
    POINTS 145 double
    -1 -1 0
    1 1 0
    -1 1 0
    1 -1 0
    0.75 1 0
    0.5 1 0
    0.25 1 0
    2.7528e-12 1 0
    -0.25 1 0
    -0.5 1 0
    -0.75 1 0
    -1 0.75 0
    -1 0.5 0
    -1 0.25 0
    -1 2.7528e-12 0
    -1 -0.25 0
    -1 -0.5 0
    -1 -0.75 0
    0.75 -1 0
    0.5 -1 0
    0.25 -1 0
    2.7528e-12 -1 0
    -0.25 -1 0
    -0.5 -1 0
    -0.75 -1 0
    1 0.75 0
    1 0.5 0
    1 0.25 0
    1 2.7528e-12 0
    1 -0.25 0
    1 -0.5 0
    1 -0.75 0
    1.03229e-12 1.37639e-12 0
    0.5 -0.5 0
    -0.5 0.5 -0
    0.5 0.5 0
    -0.5 -0.5 0
    1.89254e-12 0.5 0
    -0.25 0.25 -0
    -0.25 0.75 -0
    0.25 0.25 0
    0.25 0.75 0
    -0.5 2.0646e-12 -0
    -0.75 0.25 -0
    0.75 -0.25 0
    0.25 -0.25 0
    0.5 2.0646e-12 0
    0.75 0.25 0
    -0.75 -0.25 0
    -0.25 -0.25 0
    -0.75 0.75 -0
    0.75 -0.75 0
    -0.75 -0.75 0
    -0.25 -0.75 0
    1.89254e-12 -0.5 0
    0.25 -0.75 0
    0.75 0.75 0
    2.32267e-12 0.75 0
    -0.125 0.625 -0
    -0.125 0.875 -0
    -0.125 0.375 -0
    -0.25 0.5 -0
    1.46241e-12 0.25 0
    -0.125 0.125 -0
    -0.375 0.375 -0
    -0.375 0.625 -0
    0.375 0.375 0
    0.25 0.5 0
    0.375 0.625 0
    0.125 0.375 0
    0.125 0.625 0
    0.125 0.125 0
    0.125 0.875 0
    -0.25 1.72049e-12 -0
    -0.375 0.125 -0
    -0.625 0.125 -0
    -0.5 0.25 -0
    -0.75 2.4087e-12 -0
    -0.875 0.125 -0
    -0.625 0.375 -0
    0.875 -0.125 0
    0.625 -0.125 0
    0.75 2.4087e-12 0
    0.5 -0.25 0
    0.375 -0.125 0
    0.625 -0.375 0
    0.375 -0.375 0
    0.125 -0.125 0
    0.25 1.72049e-12 0
    0.625 0.375 0
    0.5 0.25 0
    0.625 0.125 0
    0.375 0.125 0
    0.875 0.125 0
    -0.625 -0.375 0
    -0.5 -0.25 0
    -0.375 -0.375 0
    -0.625 -0.125 0
    -0.375 -0.125 0
    -0.875 -0.125 0
    -0.125 -0.125 0
    -0.625 0.875 -0
    -0.875 0.875 -0
    -0.375 0.875 -0
    -0.5 0.75 -0
    -0.625 0.625 -0
    0.875 -0.875 0
    0.875 -0.625 0
    0.75 -0.5 0
    0.875 -0.375 0
    0.625 -0.625 0
    -0.875 -0.625 0
    -0.875 -0.875 0
    -0.875 -0.375 0
    -0.75 -0.5 0
    -0.625 -0.625 0
    -0.625 -0.875 0
    -0.5 -0.75 0
    -0.375 -0.875 0
    -0.375 -0.625 0
    -0.125 -0.875 0
    2.32267e-12 -0.75 0
    0.125 -0.625 0
    0.125 -0.875 0
    0.125 -0.375 0
    0.25 -0.5 0
    1.46241e-12 -0.25 0
    0.375 -0.625 0
    -0.125 -0.625 0
    -0.25 -0.5 0
    -0.125 -0.375 0
    -0.875 0.625 -0
    -0.75 0.5 -0
    -0.875 0.375 -0
    0.625 -0.875 0
    0.375 -0.875 0
    0.5 -0.75 0
    0.875 0.875 0
    0.625 0.875 0
    0.5 0.75 0
    0.375 0.875 0
    0.625 0.625 0
    0.875 0.625 0
    0.875 0.375 0
    0.75 0.5 0
    
    POLYGONS 256 1024
    3 7 57 59
    3 57 58 59
    3 57 37 58
    3 59 58 39
    3 37 60 58
    3 60 61 58
    3 60 38 61
    3 58 61 39
    3 37 62 60
    3 62 63 60
    3 62 32 63
    3 60 63 38
    3 39 61 65
    3 61 64 65
    3 61 38 64
    3 65 64 34
    3 35 66 68
    3 66 67 68
    3 66 40 67
    3 68 67 41
    3 40 69 67
    3 69 70 67
    3 69 37 70
    3 67 70 41
    3 40 71 69
    3 71 62 69
    3 71 32 62
    3 69 62 37
    3 41 70 72
    3 70 57 72
    3 70 37 57
    3 72 57 7
    3 32 73 63
    3 73 74 63
    3 73 42 74
    3 63 74 38
    3 42 75 74
    3 75 76 74
    3 75 43 76
    3 74 76 38
    3 42 77 75
    3 77 78 75
    3 77 14 78
    3 75 78 43
    3 38 76 64
    3 76 79 64
    3 76 43 79
    3 64 79 34
    3 28 80 82
    3 80 81 82
    3 80 44 81
    3 82 81 46
    3 44 83 81
    3 83 84 81
    3 83 45 84
    3 81 84 46
    3 44 85 83
    3 85 86 83
    3 85 33 86
    3 83 86 45
    3 46 84 88
    3 84 87 88
    3 84 45 87
    3 88 87 32
    3 35 89 66
    3 89 90 66
    3 89 47 90
    3 66 90 40
    3 47 91 90
    3 91 92 90
    3 91 46 92
    3 90 92 40
    3 47 93 91
    3 93 82 91
    3 93 28 82
    3 91 82 46
    3 40 92 71
    3 92 88 71
    3 92 46 88
    3 71 88 32
    3 36 94 96
    3 94 95 96
    3 94 48 95
    3 96 95 49
    3 48 97 95
    3 97 98 95
    3 97 42 98
    3 95 98 49
    3 48 99 97
    3 99 77 97
    3 99 14 77
    3 97 77 42
    3 49 98 100
    3 98 73 100
    3 98 42 73
    3 100 73 32
    3 2 10 102
    3 10 101 102
    3 10 9 101
    3 102 101 50
    3 9 103 101
    3 103 104 101
    3 103 39 104
    3 101 104 50
    3 9 8 103
    3 8 59 103
    3 8 7 59
    3 103 59 39
    3 50 104 105
    3 104 65 105
    3 104 39 65
    3 105 65 34
    3 3 106 31
    3 106 107 31
    3 106 51 107
    3 31 107 30
    3 51 108 107
    3 108 109 107
    3 108 44 109
    3 107 109 30
    3 51 110 108
    3 110 85 108
    3 110 33 85
    3 108 85 44
    3 30 109 29
    3 109 80 29
    3 109 44 80
    3 29 80 28
    3 0 17 112
    3 17 111 112
    3 17 16 111
    3 112 111 52
    3 16 113 111
    3 113 114 111
    3 113 48 114
    3 111 114 52
    3 16 15 113
    3 15 99 113
    3 15 14 99
    3 113 99 48
    3 52 114 115
    3 114 94 115
    3 114 48 94
    3 115 94 36
    3 0 112 24
    3 112 116 24
    3 112 52 116
    3 24 116 23
    3 52 117 116
    3 117 118 116
    3 117 53 118
    3 116 118 23
    3 52 115 117
    3 115 119 117
    3 115 36 119
    3 117 119 53
    3 23 118 22
    3 118 120 22
    3 118 53 120
    3 22 120 21
    3 21 121 123
    3 121 122 123
    3 121 54 122
    3 123 122 55
    3 54 124 122
    3 124 125 122
    3 124 45 125
    3 122 125 55
    3 54 126 124
    3 126 87 124
    3 126 32 87
    3 124 87 45
    3 55 125 127
    3 125 86 127
    3 125 45 86
    3 127 86 33
    3 21 120 121
    3 120 128 121
    3 120 53 128
    3 121 128 54
    3 53 129 128
    3 129 130 128
    3 129 49 130
    3 128 130 54
    3 53 119 129
    3 119 96 129
    3 119 36 96
    3 129 96 49
    3 54 130 126
    3 130 100 126
    3 130 49 100
    3 126 100 32
    3 2 102 11
    3 102 131 11
    3 102 50 131
    3 11 131 12
    3 50 132 131
    3 132 133 131
    3 132 43 133
    3 131 133 12
    3 50 105 132
    3 105 79 132
    3 105 34 79
    3 132 79 43
    3 12 133 13
    3 133 78 13
    3 133 43 78
    3 13 78 14
    3 3 18 106
    3 18 134 106
    3 18 19 134
    3 106 134 51
    3 19 135 134
    3 135 136 134
    3 135 55 136
    3 134 136 51
    3 19 20 135
    3 20 123 135
    3 20 21 123
    3 135 123 55
    3 51 136 110
    3 136 127 110
    3 136 55 127
    3 110 127 33
    3 1 137 4
    3 137 138 4
    3 137 56 138
    3 4 138 5
    3 56 139 138
    3 139 140 138
    3 139 41 140
    3 138 140 5
    3 56 141 139
    3 141 68 139
    3 141 35 68
    3 139 68 41
    3 5 140 6
    3 140 72 6
    3 140 41 72
    3 6 72 7
    3 1 25 137
    3 25 142 137
    3 25 26 142
    3 137 142 56
    3 26 143 142
    3 143 144 142
    3 143 47 144
    3 142 144 56
    3 26 27 143
    3 27 93 143
    3 27 28 93
    3 143 93 47
    3 56 144 141
    3 144 89 141
    3 144 47 89
    3 141 89 35
    

Anmelden zum Antworten