Skip to content
Snippets Groups Projects
Commit df929b41 authored by Nicolas Marsic's avatar Nicolas Marsic
Browse files

getdp: computing eigenvectors and dumping them as GetDP post-processing views. seems ok.

parent 2a06ad21
No related branches found
No related tags found
No related merge requests found
...@@ -17,7 +17,7 @@ def simple(operator, origin, radius, ...@@ -17,7 +17,7 @@ def simple(operator, origin, radius,
rankTol -- the tolerance for the rank test (optional) rankTol -- the tolerance for the rank test (optional)
verbose -- should I be verbose? (optional) verbose -- should I be verbose? (optional)
Returns the computed eigenvalues Returns the computed eigenvalues and eigenvectors
""" """
# Display the parameter used # Display the parameter used
...@@ -71,13 +71,14 @@ def simple(operator, origin, radius, ...@@ -71,13 +71,14 @@ def simple(operator, origin, radius,
A1 = integrate(operator, myPath, 1, vHat) A1 = integrate(operator, myPath, 1, vHat)
B = V0.H * A1 * W0 * S0Inv B = V0.H * A1 * W0 * S0Inv
# Eigenvalues of B # Eigenvalues & eigenvectors
if(verbose): print "Solving linear EVP..." if(verbose): print "Solving linear EVP..."
myLambda, QHat = numpy.linalg.eig(B) myLambda, QHat = numpy.linalg.eig(B)
Q = V0 * QHat;
# Done # Done
if(verbose): print "Done!" if(verbose): print "Done!"
return myLambda return myLambda, Q
## Import only simple (other functions are just helpers) ## Import only simple (other functions are just helpers)
......
...@@ -6,19 +6,33 @@ import numpy as np ...@@ -6,19 +6,33 @@ import numpy as np
pro = "maxwell.pro" pro = "maxwell.pro"
mesh = "square.msh" mesh = "square.msh"
resolution = "Maxwell" resolution = "Maxwell"
post = "Maxwell"
# Initialise GetDP # Initialise GetDP
operator = Solver.GetDPWave(pro, mesh, resolution) operator = Solver.GetDPWave(pro, mesh, resolution)
# Compute # Compute
v = Beyn.simple(operator, complex(9e8, 0), 1e8) l, V = Beyn.simple(operator, complex(9e8, 0), 1e8)
# Display the computed eigenvalues
print print
print "Eigenvalues:" print "Eigenvalues:"
print "----------- " print "----------- "
for i in range(v.shape[0]): for i in range(l.shape[0]):
print "(" + format(np.real(v[i]).tolist(), '+e') + ")", print " # " + str(i) + ": ",
print "(" + format(np.real(l[i]).tolist(), '+e') + ")",
print "+", print "+",
print "(" + format(np.imag(v[i]).tolist(), '+e') + ")*j" print "(" + format(np.imag(l[i]).tolist(), '+e') + ")*j"
print "----------- " 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;"
...@@ -21,7 +21,9 @@ Function{ ...@@ -21,7 +21,9 @@ Function{
d() = {}]; // DoFs d() = {}]; // DoFs
// Control data // // 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{ Jacobian{
...@@ -93,8 +95,13 @@ Resolution{ ...@@ -93,8 +95,13 @@ Resolution{
CopyRightHandSide[b(), A]; CopyRightHandSide[b(), A];
EndIf EndIf
Solve[A]; If(!viewOnly)
CopySolution[A, x()]; Solve[A];
CopySolution[A, x()];
EndIf
If(viewOnly)
CopySolution[x(), A];
EndIf
} }
} }
} }
...@@ -110,7 +117,7 @@ PostProcessing{ ...@@ -110,7 +117,7 @@ PostProcessing{
PostOperation{ PostOperation{
{ Name Maxwell; NameOfPostProcessing Maxwell; { Name Maxwell; NameOfPostProcessing Maxwell;
Operation{ Operation{
Print[e, OnElementsOf Omega, File "maxwell.pos"]; Print[e, OnElementsOf Omega, File fileName];
} }
} }
} }
...@@ -12,6 +12,8 @@ class GetDPWave: ...@@ -12,6 +12,8 @@ class GetDPWave:
b -- the right hand side (RHS) to use b -- the right hand side (RHS) to use
x -- the solution vector computed by GetDP x -- the solution vector computed by GetDP
imposeRHS -- should the imposed RHS be taken into account? 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): def __init__(self, pro, mesh, resolution):
...@@ -49,8 +51,27 @@ class GetDPWave: ...@@ -49,8 +51,27 @@ class GetDPWave:
GetDP.GetDPSetNumber("imposeRHS", 1) GetDP.GetDPSetNumber("imposeRHS", 1)
GetDP.GetDP(["getdp", self.__pro, GetDP.GetDP(["getdp", self.__pro,
"-msh", self.__mesh, "-msh", self.__mesh,
"-cal", self.__resolution, "-cal"])
"-v", "2"]) 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 @staticmethod
def __toNumpy(vGetDP): def __toNumpy(vGetDP):
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment