Seite 1 von 1

FENICS: Zeile einer Matrix auf 0 setzen

Verfasst: Montag 21. November 2016, 11:23
von MatheMagie
Hallo zusammen,

Mit Hilfe von FENICS möchte ich eine Zeile einer Matrix auf 0 setzen.
Hierzu soll die zero() Funktion verwendet werden:
https://fenicsproject.org/documentation ... atrix.html

/edit: Alternativ kann das Ganze natürlich auch mit setrow() implementiert werden.

Um meine Idee zu veranschaulichen, habe ich versucht das ganze zu programmieren:

Code: Alles auswählen

from dolfin import *
import numpy

mesh = UnitSquareMesh(50,50)
F = FunctionSpace(mesh, "Lagrange", 1)
x = TrialFunction(F)
y = TestFunction(F)

m = inner(grad(x),grad(y))*dx
M = assemble(m)

n = F.dim()
d = mesh.geometry().dim()
dof = F.tabulate_dof_coordinates().reshape(n,d)

o = numpy.array([], dtype=np.intc)
for i in xrange(0, len(dof)):
   if (i%100==0):
      u = numpy.array([i], dtype=np.intc)
      o = numpy.append(o, u)

n = len(o)
M.zero(n, o)
Allerdings erhalte ich jedemal den Fehler "NotImplementedError: Wrong number or type of arguments for overloaded function 'Matrix_zero'."
und habe keine Idee diesen zu beheben. Kann mir jemand weiterhelfen?

Re: FENICS: Zeile einer Matrix auf 0 setzen

Verfasst: Montag 21. November 2016, 11:49
von BlackJack
@MatheMagie: Rein vom Bauchgefühl würde ich sagen lass das erste Argument weg weil das keinen Sinn macht. In der C++-API ist das die Anzahl weil als zweites Argument ein Zeiger auf den ersten Index übergeben wird. In Python weiss ein Numpy-Array seine eigene Länge, da muss man die nicht noch mal zusätzlich als Argument übergeben.

Das Array erstellst Du übrigens auf extrem umständliche Weise. Die Zeilen 16 bis 20 kann man auch so schreiben: ``o = numpy.arange(0, len(dof), 100, numpy.intc)``. Ohne das dauernd temporäre Arrays, teilweise mit nur einem Element, erzeugt und im Speicher herum kopiert werden müssen.

Re: FENICS: Zeile einer Matrix auf 0 setzen

Verfasst: Montag 21. November 2016, 11:58
von MatheMagie

Code: Alles auswählen

from dolfin import *
import numpy
 
mesh = UnitSquareMesh(50,50)
F = FunctionSpace(mesh, "Lagrange", 1)
x = TrialFunction(F)
y = TestFunction(F)
 
m = inner(grad(x),grad(y))*dx
M = assemble(m)
 
n = F.dim()
d = mesh.geometry().dim()
dof = F.tabulate_dof_coordinates().reshape(n,d)
 
o = numpy.array([], dtype=numpy.intc)
for i in xrange(0, len(dof)):
   if (i%100==0):
      u = numpy.array([i], dtype=numpy.intc)
      o = numpy.append(o, u)

M.zero(o)
Vielen Dank erstmal,
Ich habe beim kopieren des Codes noch 2 kleine Fehler gefunden.
Nun habe ich das erste Argument weggelassen.
Bekomme keine Fehlermeldung, allerdings weiß ich nicht so recht, wie ich das Ganze nun überprüfen kann.
Da bei Fenics eine Ausgabe der Matrix nicht so einfach möglich ist .


Das mit dem Erstellen des Array habe ich nur so gemacht, weil die IF-Bedingung späte rnoch um einiges komplizierter werden soll.

Re: FENICS: Zeile einer Matrix auf 0 setzen

Verfasst: Montag 21. November 2016, 12:11
von MatheMagie
Ich glaube, dass es so wie du sagtest funktioniert!
Gibt es eine Möglichkeit dieses nun auch auf einen Vektor anzuwenden?
In dem Array stecken alle Indizes die 0 gesetzt werden sollen.
Der Befehl v.zero(o) geht leider nicht.

Re: FENICS: Zeile einer Matrix auf 0 setzen

Verfasst: Montag 21. November 2016, 12:57
von BlackJack
@MatheMagie: Auch wenn die ``if``-Bedingung komplizierter wird ist es keine gute Idee das so mit den Numpy-Arrays zu machen das Du das Array jedes mal um einen Wert erweiterst. Den Numpy-Arrays kann man nicht erweitern. `numpy.append()` erstellt da jedes mal ein neues Array was einen Speicherplatz länger ist, kopiert das vorherige Array in das neue Array und den neuen Wert an den letzten Index. Das ist superineffizient.

Um die ``if``-Bedingung gehören übrigens keine Klammern.

Was ist denn ein „Vektor“? Welche Methoden auf einem Objekt existieren hängt davon ab welchen Datentyp das Objekt hat.

Re: FENICS: Zeile einer Matrix auf 0 setzen

Verfasst: Montag 21. November 2016, 13:08
von MatheMagie
MatheMagie hat geschrieben:

Code: Alles auswählen

from dolfin import *
import numpy
 
mesh = UnitSquareMesh(50,50)
F = FunctionSpace(mesh, "Lagrange", 1)
x = TrialFunction(F)
y = TestFunction(F)
 
m = inner(grad(x),grad(y))*dx
M = assemble(m)
 
n = F.dim()
d = mesh.geometry().dim()
dof = F.tabulate_dof_coordinates().reshape(n,d)
 
o = numpy.array([], dtype=numpy.intc)
for i in xrange(0, len(dof)):
   if (i%100==0):
      u = numpy.array([i], dtype=numpy.intc)
      o = numpy.append(o, u)

M.zero(o)

f = Expression("5*sinx[0]")
n = f*y*dx
N = assemble(n)

a) Ich habe den Vektor in den obigen Code mit eingebaut.
Mittels assemble(n) wird der Vektor N vom Typ 'dolfin.cpp.la.Vector' erstellt.
Hier sollen nun die selben Zeilen auf 0 gesetzt werden wie in der obigen Matrix.

b) Wie würde ich das mit dem array denn effizienter machen?
Auch bei komplexerer Bedingung