Matroids Matheplanet Forum Index
Moderiert von matroid
Mathematik » Numerik & Optimierung » Implementierung der rechten Seite der FEM
Druckversion
Druckversion
Antworten
Antworten
Autor
Universität/Hochschule Implementierung der rechten Seite der FEM
Drgglbchr
Aktiv Letzter Besuch: im letzten Monat
Dabei seit: 15.11.2019
Mitteilungen: 153
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Themenstart: 2020-05-22


Hallo Leute!
Ich soll den Volumen Term für den Laplaceoperator implementieren. (Mit homogenen Dirichletrandbedingungen und Neumannrandbedingungen)
Für die Neumannbedingung habe ich einen Code:
Matlab
function b=rhsNeumann_2d_tri(element,face,vertex,Neuman_func,gqn,innerDofs,p)
 
d=innerDofs; %innere regulaere dofs
b=zeros(d,1); %lege rhs-vector an
 
[gqk, gqw]=gaussvalues(p+gqn);
gqw=gqw(:);
 
for l=1:length(face)
    if face(l).bctype==2 % ist neumann kante
        k=face(l).neighbor; % hat nur einen nachbarn
        p1=vertex(face(l).vertex(1)).coor; %linker kanten knoten
        p2=vertex(face(l).vertex(2)).coor; %rechter kanten knoten
        lengthOfFace=norm(p2-p1); %laenge der Kante
 
        espalte = element(k).connectivity;
        ezeile  = 1:length(espalte);
        eval    = ones(size(espalte));
 
        belem=compute_local_neumannRHS(element,vertex,k,Neuman_func,lengthOfFace,p1,p2,gqk,gqw,face(l).normal,p);
        for i=1:length(ezeile) % test function
            if espalte(i)<=innerDofs
                b(espalte(i))=b(espalte(i))+eval(i)*belem(ezeile(i));
            end
        end
 
    end
end
 
 
end
 
 
function belem=compute_local_neumannRHS(element,vertex,elementNr,Neuman_func,lengthOfFace,p1,p2,gqk,gqw,normal,p)
 
p1E = vertex(element(elementNr).vertex(1)).coor;
p2E = vertex(element(elementNr).vertex(2)).coor;
p3E = vertex(element(elementNr).vertex(3)).coor;
aE  =p1E;
a1E =p2E-p1E;
a2E =p3E-p1E;
 
a = (p1+p2)/2;
b = (p2-p1)/2;
 
tmp=[a(1)+b(1)*gqk
     a(2)+b(2)*gqk];
 
gqkE=[a1E a2E]\(tmp-aE);
gqkE=gqkE';
basis=evallocalbasis(98,gqkE,p);
 
tmpN=Neuman_func(tmp(1,:),tmp(2,:),normal).*gqw;
belem=basis*tmpN*lengthOfFace/2;
 
end

jetzt war meine Idee, diesen einfach umzuschreiben....
also bctype == 0 -> ds wäre dann eine innere Kante
aber wie kann ich das mit den Nachbarn umsetzen?

lg Drgglbchr



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
Folgende Antworten hat der Fragesteller vermutlich noch nicht gesehen.
Monom
Aktiv Letzter Besuch: im letzten Monat
Dabei seit: 29.05.2015
Mitteilungen: 98
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.1, eingetragen 2020-05-24


Hallo Drgglbchr,

der Code ist für die Randintegrale auf dem Neumann-Rand des Gebiets.
Es werden also über bestimmte Dreieckskanten integriert (bctype == 2), die am Rand des Gebiets liegen.

Für einen Volumenterm müsste über ganze Dreiecke und nicht bloß über eine (oder zwei) ihrer Seiten integriert werden. Es macht also wenig Sinn, diesen Code als Vorlage zu nehmen.

Viele Grüße
Monom



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]