FENICS: Zeile einer Matrix auf 0 setzen

mit matplotlib, NumPy, pandas, SciPy, SymPy und weiteren mathematischen Programmbibliotheken.
Antworten
MatheMagie
User
Beiträge: 7
Registriert: Montag 21. November 2016, 11:16

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?
Zuletzt geändert von Anonymous am Montag 21. November 2016, 11:34, insgesamt 1-mal geändert.
Grund: Quelltext in Python-Codebox-Tags gesetzt.
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.
MatheMagie
User
Beiträge: 7
Registriert: Montag 21. November 2016, 11:16

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.
MatheMagie
User
Beiträge: 7
Registriert: Montag 21. November 2016, 11:16

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.
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.
MatheMagie
User
Beiträge: 7
Registriert: Montag 21. November 2016, 11:16

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
Antworten