Matroids Matheplanet Forum Index
Moderiert von matroid
Mathematik » Numerik & Optimierung » Steifigkeitsmatrix FEM
Druckversion
Druckversion
Antworten
Antworten
Autor
Universität/Hochschule Steifigkeitsmatrix FEM
Drgglbchr
Aktiv Letzter Besuch: in der letzten Woche
Dabei seit: 15.11.2019
Mitteilungen: 154
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Themenstart: 2020-05-14


hallo!
ich bräuchte bitte eine Erklärung, wie man bei der FEM für ein L-förmiges gebiet $\Omega = (-1,1)^2/[0,1]^2$ die Matrix A berechnen kann, wobei auf dem "herausgeschnittenen" teil Dirichletdaten und sonst Neumanndaten zur Anwendung kommen.
also
$-\triangle u = 0 \qquad in\ \Omega$
$u=0\qquad auf\ \Gamma_D$
$\partial / {\partial n}\ u=g \qquad auf\ \Gamma_N$

im internet konnte ich leider nichts finden, was mir weiterhilft :/



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
Carmageddon
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 22.12.2009
Mitteilungen: 650
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.1, eingetragen 2020-05-14


Hallo Drgglbchr,


wenn man gemischte Randbedingungen hat (keine Robin Boundary Conditions) ist das Vorgehen im Prinzip genauso, nur muss man die Testräume verändern:

Sei <math>H_{D_0} := \{ u \in H^1 (\Omega): \; u = 0 \; auf \; \Gamma_D \}</math>

So ist das Variationsproblem gegeben als

Finde <math>u \in H_{D_0}</math>, so dass <math>\int\limits_\Omega \nabla u \cdot \nabla v = \int\limits_{\Gamma_N} v g</math> für alle <math>v \in H_{D_0}</math> gilt.

Damit kannst du dir nun deine Matrizen berechnen - beachte allerdings den geänderten Raum der Testfunktionen im Vergleich zu durchgehenden Dirichlet Daten - das musst du natürlich in der Numerik übernehmen!

Grob gesagt liegt das Problem "zwischen" einem reinem Dirichlet und einem reinem Neumann-Problem. Genauso verhält es sich auch für den Raum der Testfunktionen.

Quelle z.b.:
http://people.inf.ethz.ch/arbenz/FEM17/pdfs/0-19-852868-X.pdf



-----------------
Zitat: "Es gibt einen Beweis aus der Physik: Er ist kurz, er ist elegant... und falsch"



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
Drgglbchr
Aktiv Letzter Besuch: in der letzten Woche
Dabei seit: 15.11.2019
Mitteilungen: 154
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.2, vom Themenstarter, eingetragen 2020-05-16


vielen dank für deine Antwort :)

in der Theorie ist mir das alles eigentlich klar. ich habe vor allem noch Schwierigkeiten bei der numerischen Umsetzung.
aber mir ist jetzt klar, dass ich die matrix A als summe der lokalen steifigkeitsmatrizen schreiben kann, wobei ich die Dreiecke auf ein referenzdreieck überführen kann mithilfe der affinen Abbildung $\Phi (s,t) = x_1 + s(x_2-x_1)+t(x_3-x_1)$ - das müsste die Idee für die Implementierung sein, oder? :)

zur rechten Seite (also dem Volumen term) hätte ich auch noch eine frage:

in meiner vorläge für die Implementierung werden nämlich 2 Funktionen für die randterme verwendet.
also b = rhs_Volume_2d_tri(...)
b = b + rhs_Neumann....
wie könnte ich rhs_volume implementieren? bzw wäre das in meinem Beispiel nicht eh Null?🤯
Matlab
addpath('input_func')
addpath('error')
addpath('matrix_rhs')
addpath('mesh')
 
 
clear;
close all;
clc
 
 
case_flag=9; % 0 - 4 = poly Bsp, 5 = sin, 6 free, 7 = exp, 8 = singular bsp
% homo Dirichlet: case 4,7,8
 
basicName=['data/lshape_case_',num2str(case_flag),'_adaptive_h_version'];
scaling=1; % skalierung des gebiets = scaling*[-1,1]^2
 
draw_mesh_flag = 0; % draw mesh
draw_sol_flag  = 0; % draw solution
draw_eta_flag  = 0; % draw error est distribution
 
n     = 2; % n= nr. of Elements per edge => n^2 Elemente,
p     = 1; % p=polynomgrad
gqn_f = 5; % zusaetzliche Quadraturknoten
 
k     = 15;  % anzahl der Gitterverfeinerungen
theta = 0.5; % groesser gleich mehr markiert
 
 
[element,vertex]      = MeshGenerator2d_tri_simple_Lshape(n,scaling);
[element,face,vertex] = generate_faces2(element,vertex);
element               = check_local_faceNr(element,face);
marking               = ones(length(element),3);
 
for i=1:k
    tic
 
    disp('----------------------------')
    disp(['      Durchlauf Nr. ', num2str(i)] )
    disp('----------------------------')
 
    name=[basicName,'_p',num2str(p),'_meshNr',num2str(i)];
 
    % generate mesh
    marking               = mark_conform(marking,element,face,vertex);
    [element,vertex]      = refine_simple(marking,element,face,vertex);
    [element,face,vertex] = generate_faces2(element,vertex);
    element               = check_local_faceNr(element,face);
    disp('mesh done')
 
    % generate connectivity data for u = hat functions
    [element,face,vertex,innerDofs,totaldof]=generate_connectivity_data(element,face,vertex,p);
    disp('dofs done')
 
    if draw_mesh_flag==1
        fig_handle=drawMesh(element,face,vertex,0);
%         fig_handle=drawMesh_lshape_zoom(element,face,vertex)
 
        axis equal
        axis off
        saveas(fig_handle,[name,'_mesh.fig'],'fig')
        saveas(fig_handle,[name,'_mesh.jpg'],'jpg')
        close(fig_handle);
        disp('draw mesh done')
    end
 
    % Interpolation der Dirichlet Daten u
    G_u=dirichletdaten_interpolation_u(element,face,vertex,@(x,y)func_u(x,y,case_flag),innerDofs, totaldof,p);
 
    % volumen terms, rhs
    b=rhsVolume_2d_tri(element,vertex,@(x,y)volumeFunction(x,y,case_flag),gqn_f,innerDofs,p); 
    b=b+rhsNeumann_2d_tri(element,face,vertex,@(x,y,normal)func_NormalGrad_u(x,y,normal,case_flag),gqn_f,innerDofs,p);
    disp('rhs done')
 
    % matrizen for u and sigma
    [A,b]=stiffnessMatrix_2d_tri_opt(element,vertex,innerDofs, b,G_u,p);    
    disp('matrix done')
 
 
    % solve linear system
    % symmetric positive definite
    u=A\b;
 
    % include Dirichlet Data
    u_all=G_u;
    u_all(1:innerDofs)=u;
 
    dof(i,1)   = length(u);
    disp('solver done')
 
    if draw_sol_flag==1
        % u
        fig_handle=drawSolution(element,vertex,u_all,3,p);
        xlabel('x')
        ylabel('y')
        axis auto
        view([-16 36])
        saveas(fig_handle,[name,'_sol_u.fig'],'fig')
        saveas(fig_handle,[name,'_sol_u.jpg'],'jpg')
        close(fig_handle);
 
        % grad u
        fig_handle=drawSolution_gradient(element,face,vertex,u_all,0,scaling,p);
        saveas(fig_handle,[name,'_sol_grad_u.fig'],'fig')
        saveas(fig_handle,[name,'_sol_grad_u.jpg'],'jpg')
        close(fig_handle);        
    end
 
    % Fehlerschaetzer ohne Dirichlet-Daten-Fehler
    [eta_res_loc,eta_res_edge_loc,eta_Diri_approx]=compute_residual_error_estimate(element,face,vertex, @(x,y)volumeFunction(x,y,case_flag), @(x,y,normal)func_NormalGrad_u(x,y,normal,case_flag), u_all, gqn_f,p);
    eta_res(i,1)       = sqrt(sum(eta_res_loc));
    eta_res_edge(i,1)  = sqrt(sum(eta_res_edge_loc));    
    eta_error_est(i,1) = sqrt(sum(eta_res_loc + eta_res_edge_loc) )
 
 
    disp('error est done')
 
 
    if draw_eta_flag==1
        fig_handle=drawEta(element,vertex,eta_res_loc + eta_res_edge_loc);
        xlabel('x')
        ylabel('y')
        axis auto
        colorbar
        saveas(fig_handle,[name,'_etaDistr.fig'],'fig')
        saveas(fig_handle,[name,'_etaDistr.jpg'],'jpg')
        close(fig_handle);
    end  
 
    % Fehler Berechnung
    [L2_error_u(i,1),H1_error_u(i,1)]=error_u_func(element,vertex,u_all,@(x,y)func_u(x,y,case_flag),@(x,y)func_grad_u(x,y,case_flag),gqn_f,p);
 
    % Doerfler-Marking
    marking=markElement(eta_res_loc + eta_res_edge_loc,theta);
 
    save(name)
    time(i)=toc
end
 
figure
loglog(dof,H1_error_u,'-+')



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
Monom
Aktiv Letzter Besuch: in der letzten Woche
Dabei seit: 29.05.2015
Mitteilungen: 99
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.3, eingetragen 2020-05-24


Hallo Drgglbchr,

die rechte Seite der PDE ist 0. Daher gibt es keinen Volumenterm, wie Du festgestellt hast.

Viele Grüße
Monom

PS: Aus Interesse: Was für ein Matlab FEM Code ist das?



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
Drgglbchr
Aktiv Letzter Besuch: in der letzten Woche
Dabei seit: 15.11.2019
Mitteilungen: 154
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.4, vom Themenstarter, eingetragen 2020-05-25


Hallo Monom!
Den Code habe ich von meinem Prof erhalten -  ich soll ihn für mein spezifisches Problem fertigstellen.

Zu der Steifigkeitsmatrix hätte ich noch eine Frage:
Ich habe versucht die Funktion aelem zu komplettieren, aber die Dimensionen stimmen nicht... habe schon paar Versionen durchprobiert... :/
Und vl kannst du mir auch noch verraten, ob meine lokale Steifigkeitsmatrix so korrekt wäre? 😉
Matlab
function [A,b]=stiffnessMatrix_2d_tri_opt(element,vertex,innerDofs, b,G, p)
 
d=innerDofs;
 
[gqk, gqw]=dunavant_triangle_quad(2);
basisPrime=evallocalbasis_prime(98,gqk,p);
basisPrimeX=basisPrime{1};
basisPrimeY=basisPrime{2};
 
cnt=0;
for k=1:length(element) %over all possible elements
    espalte=element(k).connectivity;
    cnt=cnt+sum(espalte<=innerDofs)^2;
end
 
e1=zeros(1,cnt); e2=zeros(1,cnt); e3=zeros(1,cnt);
 
 
cnt=1;
for k=1:length(element) %over all possible elements
    aelem=compLocalStiffnessMatrix2d(element,vertex,k,p, gqk,gqw,basisPrimeX,basisPrimeY);
 
    espalte = element(k).connectivity;
    ezeile  = 1:length(espalte);
    eval    = ones(size(espalte));                                                                                                                           
 
    for i=1:length(ezeile) % test function
        for j=1:length(ezeile) % ansatz function
            if espalte(i)<=innerDofs && espalte(j)<=innerDofs
                e1(cnt)=espalte(i);
                e2(cnt)=espalte(j);
                e3(cnt)=eval(i)*eval(j)*aelem(ezeile(i),ezeile(j));
                cnt=cnt+1;
            end
 
            if espalte(i)>innerDofs && espalte(j)<=innerDofs
                b(espalte(j))=b(espalte(j))-G(espalte(i))*eval(i)*eval(j)*aelem(ezeile(i),ezeile(j));
            end
 
        end
    end
end
 
A=sparse(e1(1:cnt-1),e2(1:cnt-1),e3(1:cnt-1),d,d);
end
 
 
function aelem=compLocalStiffnessMatrix2d(element,vertex,elementNr,p, gqk,gqw,basisPrimeX,basisPrimeY)
 
p1=vertex(element(elementNr).vertex(1)).coor;
p2=vertex(element(elementNr).vertex(2)).coor;
p3=vertex(element(elementNr).vertex(3)).coor;
 
% a  =p1; 
a1 =p2-p1;
a2 =p3-p1;
%-------------
a = p1;
kreuz = 2*abs(a1(1)*a2(2)-a1(2)*a2(1));
 
aelem = [(a2*a2.'+a1*a1.'-2*a1*a2.')/kreuz (-(a2*a2.')+a1*a2.')/kreuz (-(a1*a1.')+a1*a2.')/kreuz,0 (a2*a2.')/kreuz -a1*a2.'/kreuz, 0 0 (a1*a1.')/kreuz]
 
 
 
 
 
 
end



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
Folgende Antworten hat der Fragesteller vermutlich noch nicht gesehen.
Carmageddon
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 22.12.2009
Mitteilungen: 650
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.5, eingetragen 2020-05-25


2020-05-25 09:59 - Drgglbchr in Beitrag No. 4 schreibt:
Hallo Monom!
Den Code habe ich von meinem Prof erhalten -  ich soll ihn für mein spezifisches Problem fertigstellen.

Zu der Steifigkeitsmatrix hätte ich noch eine Frage:
Ich habe versucht die Funktion aelem zu komplettieren, aber die Dimensionen stimmen nicht... habe schon paar Versionen durchprobiert... :/
Und vl kannst du mir auch noch verraten, ob meine lokale Steifigkeitsmatrix so korrekt wäre? 😉

Hallo,

Wo genau fliegt der Code denn raus? Es ist nicht so einfach, deinen Code zu debuggen, wenn man die ganzen zusätzlichen Infos (wie Grid-Infos, etc) nicht hat.

Ansonsten ist es eine gute Idee, sich die Matrix mal ausgeben zu lassen - wenn man eine lexikographische Sortierung verwendet, sollte man eine gute Struktur erkennen können und dort eventuell Fehler finden.
Der entsprechende Matlab-Befehl lautet "spy".

Viele Grüße


-----------------
Zitat: "Es gibt einen Beweis aus der Physik: Er ist kurz, er ist elegant... und falsch"



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
Drgglbchr wird per Mail über neue Antworten informiert.
Neues Thema [Neues Thema] Antworten [Antworten]    Druckversion [Druckversion]

 


Wechsel in ein anderes Forum:
 Suchen    
 
All logos and trademarks in this site are property of their respective owner. The comments are property of their posters, all the rest © 2001-2020 by Matroids Matheplanet
This web site was originally made with PHP-Nuke, a former web portal system written in PHP that seems no longer to be maintained nor supported. PHP-Nuke is Free Software released under the GNU/GPL license.
Ich distanziere mich von rechtswidrigen oder anstößigen Inhalten, die sich trotz aufmerksamer Prüfung hinter hier verwendeten Links verbergen mögen.
Lesen Sie die Nutzungsbedingungen, die Distanzierung, die Datenschutzerklärung und das Impressum.
[Seitenanfang]