From df929b41e77a8428764aaecea7ea78d11f338982 Mon Sep 17 00:00:00 2001
From: Nicolas Marsic <marsic@temf.tu-darmstadt.de>
Date: Wed, 21 Sep 2016 17:00:05 +0000
Subject: [PATCH] getdp: computing eigenvectors and dumping them as GetDP
post-processing views. seems ok.
---
beyn.py | 7 ++++---
main.py | 24 +++++++++++++++++++-----
maxwell.pro | 15 +++++++++++----
solver.py | 25 +++++++++++++++++++++++--
4 files changed, 57 insertions(+), 14 deletions(-)
diff --git a/beyn.py b/beyn.py
index 2126782..8d9a197 100644
--- a/beyn.py
+++ b/beyn.py
@@ -17,7 +17,7 @@ def simple(operator, origin, radius,
rankTol -- the tolerance for the rank test (optional)
verbose -- should I be verbose? (optional)
- Returns the computed eigenvalues
+ Returns the computed eigenvalues and eigenvectors
"""
# Display the parameter used
@@ -71,13 +71,14 @@ def simple(operator, origin, radius,
A1 = integrate(operator, myPath, 1, vHat)
B = V0.H * A1 * W0 * S0Inv
- # Eigenvalues of B
+ # Eigenvalues & eigenvectors
if(verbose): print "Solving linear EVP..."
myLambda, QHat = numpy.linalg.eig(B)
+ Q = V0 * QHat;
# Done
if(verbose): print "Done!"
- return myLambda
+ return myLambda, Q
## Import only simple (other functions are just helpers)
diff --git a/main.py b/main.py
index cbb1e25..de97335 100644
--- a/main.py
+++ b/main.py
@@ -6,19 +6,33 @@ import numpy as np
pro = "maxwell.pro"
mesh = "square.msh"
resolution = "Maxwell"
+post = "Maxwell"
# Initialise GetDP
operator = Solver.GetDPWave(pro, mesh, resolution)
# Compute
-v = Beyn.simple(operator, complex(9e8, 0), 1e8)
+l, V = Beyn.simple(operator, complex(9e8, 0), 1e8)
+# Display the computed eigenvalues
print
print "Eigenvalues:"
print "----------- "
-for i in range(v.shape[0]):
- print "(" + format(np.real(v[i]).tolist(), '+e') + ")",
+for i in range(l.shape[0]):
+ print " # " + str(i) + ": ",
+ print "(" + format(np.real(l[i]).tolist(), '+e') + ")",
print "+",
- print "(" + format(np.imag(v[i]).tolist(), '+e') + ")*j"
+ print "(" + format(np.imag(l[i]).tolist(), '+e') + ")*j"
print "----------- "
-print "Found: " + str(v.shape[0])
+print "Found: " + str(l.shape[0])
+
+# Generate GetDP views for eigenvectors
+print
+print "Generating eigenvector views..."
+for i in range(l.shape[0]):
+ operator.view(V[:, i], post, "beyn" + str(i) + ".pos")
+ print " # View: " + str(i) + " done!"
+
+# Done
+print
+print "return 0;"
diff --git a/maxwell.pro b/maxwell.pro
index 9802e26..eca4aa3 100644
--- a/maxwell.pro
+++ b/maxwell.pro
@@ -21,7 +21,9 @@ Function{
d() = {}]; // DoFs
// Control data //
- DefineConstant[imposeRHS = 0]; // Should I use an imposed RHS?
+ DefineConstant[imposeRHS = 0, // Should I use an imposed RHS?
+ viewOnly = 0, // Should I use X for postpro only?
+ fileName = "eig.pos"]; // Postpro file name
}
Jacobian{
@@ -93,8 +95,13 @@ Resolution{
CopyRightHandSide[b(), A];
EndIf
- Solve[A];
- CopySolution[A, x()];
+ If(!viewOnly)
+ Solve[A];
+ CopySolution[A, x()];
+ EndIf
+ If(viewOnly)
+ CopySolution[x(), A];
+ EndIf
}
}
}
@@ -110,7 +117,7 @@ PostProcessing{
PostOperation{
{ Name Maxwell; NameOfPostProcessing Maxwell;
Operation{
- Print[e, OnElementsOf Omega, File "maxwell.pos"];
+ Print[e, OnElementsOf Omega, File fileName];
}
}
}
diff --git a/solver.py b/solver.py
index 6fe6946..eb01d1e 100644
--- a/solver.py
+++ b/solver.py
@@ -12,6 +12,8 @@ class GetDPWave:
b -- the right hand side (RHS) to use
x -- the solution vector computed by GetDP
imposeRHS -- should the imposed RHS be taken into account?
+ viewOnly -- should x be used for post-processing only
+ fileName -- post-processing file name
"""
def __init__(self, pro, mesh, resolution):
@@ -49,8 +51,27 @@ class GetDPWave:
GetDP.GetDPSetNumber("imposeRHS", 1)
GetDP.GetDP(["getdp", self.__pro,
"-msh", self.__mesh,
- "-cal", self.__resolution,
- "-v", "2"])
+ "-cal"])
+ GetDP.GetDPSetNumber("imposeRHS", 0)
+
+ def view(self, x, posName, fileName):
+ """Generates a post-processing view
+
+ Keyword arguments:
+ x -- the solution vector to use
+ posName -- the post-operation to use (from the .pro file)
+ fileName -- the post-precessing file name
+
+ This method generates a linear system
+ """
+ GetDP.GetDPSetNumber("x", self.__toGetDP(x))
+ GetDP.GetDPSetNumber("viewOnly", 1)
+ GetDP.GetDPSetString("fileName", fileName)
+ GetDP.GetDP(["getdp", self.__pro,
+ "-msh", self.__mesh,
+ "-cal",
+ "-pos", posName])
+ GetDP.GetDPSetNumber("viewOnly", 0)
@staticmethod
def __toNumpy(vGetDP):
--
GitLab