diff --git a/bin/args.py b/bin/args.py
index 3a911204878ccb77a6ae782c435159241296beb2..aa1d113423719584d0dcbfa3946e9555f48b7fc3 100644
--- a/bin/args.py
+++ b/bin/args.py
@@ -50,15 +50,6 @@ def parse():
help="set constant number NAME=VALUE")
parser.add_argument("-quiet", action='store_true',
help="should I be quiet?")
- parser.add_argument("-non_holomorphic", action='store_true',
- help='start a fixed-point loop for non-holomorphic '+
- 'functions')
- parser.add_argument("-nhMaxIt", type=int, default=5,
- help='maximum iteration of non-holomorphic solver')
- parser.add_argument("-nhMinRratio", type=float, default=10,
- help='minimum ratio between non-holomorphic and ' +
- 'holomorphic residuals for ending fixed-point ' +
- 'iterations')
# Done
return parser.parse_args()
diff --git a/bin/cim.py b/bin/cim.py
index 4d9a1c61d81b6f26ca8b1499d522e7fac109560d..f9d52d11534e8f7b6eb6e1da5256b8e3eeb6599c 100755
--- a/bin/cim.py
+++ b/bin/cim.py
@@ -13,12 +13,11 @@ try:
except ImportError:
from mpiseq import MPI # If not available use a dummy mpiseq
-from solver import GetDPWave
-from beyn import Beyn
-from fixed_point import FixedPoint
-import args as args
-import numpy as np
-import time as timer
+from solver import GetDPWave
+from beyn import Beyn
+import args as args
+import numpy as np
+import time as timer
# Start timer
tStart = timer.time()
@@ -46,18 +45,8 @@ beyn = Beyn(operator, arg.origin, arg.radius,
not arg.quiet)
# Compute
-if(not arg.non_holomorphic):
- if(not arg.quiet): beyn.display()
- l, V, r = beyn.simple()
-
-else:
- fixed = FixedPoint(beyn, arg.nhMaxIt, arg.nhMinRratio, not arg.quiet)
- # This variable ('fixed') is meant to be called
- if(verbose): # elsewhere (i.e. Python[] call in GetDP) thanks
- fixed.display() # to the python interpreter.
- beyn.display() # Please don't change its name!
-
- l, V, r = fixed.solve()
+if(not arg.quiet): beyn.display()
+l, V, r = beyn.simple()
# Display the computed eigenvalues
if(verbose):
diff --git a/bin/fixed_point.py b/bin/fixed_point.py
deleted file mode 100644
index 284b743ab89e6405508861dda01b43afabeac1d5..0000000000000000000000000000000000000000
--- a/bin/fixed_point.py
+++ /dev/null
@@ -1,196 +0,0 @@
-"""
-cim.py, a non-linear eigenvalue solver.
-Copyright (C) 2017 N. Marsic, F. Wolf, S. Schoeps and H. De Gersem,
-Institut fuer Theorie Elektromagnetischer Felder (TEMF),
-Technische Universitaet Darmstadt.
-
-See the LICENSE.txt and README.md for more license and copyright information.
-"""
-
-"""
-TODO: cluster identification / can a spline be holomorphic?
-"""
-
-try:
- from mpi4py import MPI # Try to import mpi4pi
-except ImportError:
- from mpiseq import MPI # If not available use a dummy mpiseq
-
-import sys as sys
-import numpy as np
-
-class FixedPoint:
- """A fixed-point solver for handling non-holomorphic operator
-
- This solver allows to approximate a non-holomorphic operator by an
- holomorphic one. This approximation is refined by a fixed-point iteration,
- until the residual of the original (non-holomorphic) problem reaches the
- residual of the corrected (holomorphic) problem. The holomorphic correction
- is preformed by calling a python handler from GetDP (using the Python[]
- function), which is responsible for implementing the holomorphic correction.
- The handler can communicate with this fixed-point solver. By convention,
- this solver is refred to as 'fixed', and is accessible from the handler.
- An example is provided in: ./example/maxwell_sibc_lagrange/
- """
-
- def __init__(self, solver, maxIt=5, minNHHRatio=10, verbose=True):
- """Instantiates a new FixedPoint non-holomorphic solver
-
- Keyword arguments:
- solver -- the solver used to compute holomorphic solutions
- maxIt -- the maximum number of fixed-point iteration (optional)
- minNHHRatio -- the solver stopping criterion (see below, optional)
- verbose -- should I be verbose? (optional)
-
- What is minNHHRatio?
- The FixedPoint solver computes two residuals for each solution it finds:
- 1. rH, computed with respect to the corrected, holomorphic, operator;
- 2. rNH, computed with respect to the original, non-holomorphic, operator
- Ultimately, the FixedPoint solver tries to impose:
- rNH = rH,
- for all eigenpairs. Thus, we use max(rNH/rH) as an indicator of the
- overall accuracy. Therefore, the FixedPoint solver will stop as soon as:
- max(rNH/rH) > minNHHRatio,
- or if the number of iteration exceeds maxIt. In other words, the closer
- is minNHHRatio to 1, the sharper the solution. Please note that it does
- not make sense to have minNHHRatio < 1. That is, the corrected
- holomorphic solution cannot be more precise that the original
- non-holomorphic one.
- """
- # Save
- self.__solver = solver
- self.__maxIt = maxIt
- self.__minNHHRatio = minNHHRatio
- self.__myProc = MPI.COMM_WORLD.Get_rank()
- self.__verbose = verbose and (self.__myProc == 0)
-
- # State
- self.__cluster = []
- self.__holomorphic = True
-
- def solve(self):
- """Solves the eigenvalue problem by a fixed-point iteration
-
- It returns the computed eigenvalues, eigenvectors and the norm of the
- residual associated to the original, non-holomorphic, operator
- """
- # Pretty copy
- verbose = self.__verbose
-
- # (Re)initialize
- self.__cluster = []
- self.__holomorphic = True
-
- # First run with corrector initial guess
- if(verbose): print ">> Non-holomorphic fixed-point solver: initial run"
-
- l, V, rH = self.__solver.simple() # Holomorphic solution
- rNH = self.__nhNorm(l, V) # Non-holomorphic residual
- maxRNHH = self.__maxRNHH(rNH, rH) # Maximim rNH/rH
-
- if(verbose): print rH
- if(verbose): print rNH
-
- if(verbose): print ">> Initial max(rNH/rH): " + format(maxRNHH, ".2e")
-
- # Fixed-point iteration
- it = 0
- if(verbose): print ">> Non-holomorphic fixed-point solver: fixed-point"
- while(maxRNHH > self.__minNHHRatio and it != self.__maxIt):
- self.__cluster = self.__findCluster(l) # Clusters for correction
- l, V, rH = self.__solver.simple() # Solve corrected system
- rNH = self.__nhNorm(l, V) # Norm of NH residual
- maxRNHH = self.__maxRNHH(rNH, rH) # Maximim rNH/rH
-
- if(verbose): print rH
- if(verbose): print rNH
-
- if(verbose): print ">> Fixed-point iteration " + str(it+1),
- if(verbose): print "done: " + format(maxRNHH, ".2e")
-
- it = it + 1 # One more iteration :-)
-
- # Warn if too many fixed-point iterations
- if(it == self.__maxIt):
- if(verbose): ">> Maximum number of fixed-point iterations reached"
-
- # Return last solution with the norm of the non-holomorphic residual
- return l, V, rNH
-
- def cluster(self):
- """Returns the eigenvalues (by cluster) found by the last iteration
-
- This method returns the eigenvalues "by cluster". That is, it some
- eigenvalues are too close from each other, only one will be considered
- """
- return self.__cluster
-
- def needNonHolomorphic(self):
- """This method returns True if this FixedPoint needs the true,
- non-holomorphic, operator and not its holomorphic approximation"""
- return not self.__holomorphic
-
- def display(self):
- """Prints this FixedPoint solver parameters"""
- print "Fixed-point non-holomorphic correction "
- print "---------------------------------------"
- print ">> Maximum number of interation: "+ " " + str(self.__maxIt)
- print ">> Minimum NH-H residual ratio: ",
- print format(self.__minNHHRatio, ".2g")
- print "---------------------------------------"
-
- def __nhNorm(self, l, V):
- """Returns the norm of the non-holomorphic residual
- for all eigenpairs (l, V)"""
- operator = self.__solver.operator() # Get operator
-
- nFound = l.shape[0] # Number of found eigenpair
- nhNorm = np.empty(nFound) # Non-holomorphic residual
-
- self.__holomorphic = False # Non-holomorphic (NH) operator
- for n in range(nFound): # Iterate on solutions
- operator.apply(V[:, n], l[n]) # -> Apply solution on NH operator
- nhRes = operator.solution(0) # -> residual
- nhNorm[n] = np.linalg.norm(nhRes) # -> norm
- self.__holomorphic = True # Back to holomorphic operator
-
- return nhNorm # Done
-
- @staticmethod
- def __findCluster(data):
- """Retruns the cluster contained in the given *sorted* list"""
- rank = MPI.COMM_WORLD.Get_rank()
- # Get clusters
- cluster = []
- delta = 0.01
- i = 1
- while i < len(data):
- j = i
- sub = [data[i-1]]
- while j < len(data):
- d = (data[i-1] - data[j]) / data[j]
- if abs(d.real) < delta and abs(d.imag) < delta:
- sub.append(data[j])
- j = j + 1
- else:
- break
-
- cluster.append(sub)
- i = j + 1
-
- # Mean of each cluster
- out = []
- for i in range(len(cluster)):
- out.append(np.mean(cluster.pop()))
- out.reverse()
-
- if(rank == 0): print "Cluster:",
- if(rank == 0): print out
-
- return out
-# return [data[0], data[3], data[6], data[11]]
-
- @staticmethod
- def __maxRNHH(rNH, rH):
- """ Returns the maximal rNH/rH"""
- return max([a/b for a, b in zip(rNH, rH)])
diff --git a/bin/interpolator.py b/bin/interpolator.py
deleted file mode 100644
index 61d63aa97e11e668f3acbd1979babefe8345d10a..0000000000000000000000000000000000000000
--- a/bin/interpolator.py
+++ /dev/null
@@ -1,33 +0,0 @@
-"""
-cim.py, a non-linear eigenvalue solver.
-Copyright (C) 2017 N. Marsic, F. Wolf, S. Schoeps and H. De Gersem,
-Institut fuer Theorie Elektromagnetischer Felder (TEMF),
-Technische Universitaet Darmstadt.
-
-See the LICENSE.txt and README.md for more license and copyright information.
-"""
-
-from numpy import poly1d
-from scipy import interpolate
-from abc import ABCMeta, abstractmethod
-
-class Interpolator:
- """A generic class for interpolation"""
- @abstractmethod
- def __init__(self):
- pass
-
- def at(self, x):
- """Returns the value of this Interpolation at a given point"""
- return self._p(x) # We use poly1d from scipy
-
- __metaclass__ = ABCMeta
-
-
-class Lagrange(Interpolator):
- """A class for Lagrange interpolation"""
- def __init__(self, xdata, ydata):
- """Instantiates a new Interpolator, such that:
- p(xdata) = ydata for all xdata, where p is build with a Lagrange basis
- """
- self._p = interpolate.lagrange(xdata, ydata)
diff --git a/example/maxwell_sibc_lagrange/Makefile b/example/maxwell_sibc_lagrange/Makefile
deleted file mode 100644
index b54d01b5f6c88d726eebce840a80185e2eb26a41..0000000000000000000000000000000000000000
--- a/example/maxwell_sibc_lagrange/Makefile
+++ /dev/null
@@ -1,21 +0,0 @@
-BIN=../../bin
-
-sphere: init cim
-
-init:
- gmsh sphere.geo -3
-
-cim:
-# python ${BIN}/cim.py sphere.pro sphere.msh Solve 1e10+1e9j 4e9 -nodes 50 -lStart 20 -lStep 30 -setnumber isSC 0 -non_holomorphic
-
- python ${BIN}/cim.py sphere.pro sphere.msh Solve 7.4e9+7.3e8j 1e9 -nodes 10 -lStart 4 -setnumber isSC 0 -non_holomorphic
-#7.474355560164e+09 + 7.742536715699e+08j
-
-# python ${BIN}/cim.py sphere.pro sphere.msh Solve 1.0e10+1.1e9j 1e9 -nodes 10 -lStart 6 -setnumber isSC 0 -non_holomorphic
-#1.050336615299e+10 + 1.186115651023e+09j
-
-clean:
- rm -f *.pos
- rm -f *.msh
- rm -f *.pre
- rm -f *.res
diff --git a/example/maxwell_sibc_lagrange/cimParameters.pro b/example/maxwell_sibc_lagrange/cimParameters.pro
deleted file mode 100644
index 80b91c50abbb3e88012b971fd11ed04ad33c876a..0000000000000000000000000000000000000000
--- a/example/maxwell_sibc_lagrange/cimParameters.pro
+++ /dev/null
@@ -1,26 +0,0 @@
-/////////////////////////////////////////////////////////////////
-// Parameters for contour integral method. //
-// The master file should inclue "cimResolution.pro" later on. //
-/////////////////////////////////////////////////////////////////
-Function{
- // Physical data //
- DefineConstant[angularFreqRe = 1, // [rad/s]
- angularFreqIm = 0]; // [rad/s]
-
- Puls[] = Complex[angularFreqRe,
- angularFreqIm]; // [rad/m]
-
- // Algebraic data //
- DefineConstant[nRHS = 1]; // Number of RHS for this run
- For i In {0:nRHS-1}
- DefineConstant[x~{i}() = {}, // Solution
- b~{i}() = {}]; // Right hand side
- EndFor
-
- // Control data //
- DefineConstant[doInit = 0, // Should I initialize some stuff?
- doSolve = 0, // Should I solve Ax = b?
- doPostpro = 0, // Should I only create a view for x()?
- doApply = 0, // Should I only apply x(): x <- Ax?
- fileName = "eig.pos"]; // Postpro file name
-}
diff --git a/example/maxwell_sibc_lagrange/cimResolution.pro b/example/maxwell_sibc_lagrange/cimResolution.pro
deleted file mode 100644
index 186472cc6b6257ff9ebc6a4807a18802ad25bcaa..0000000000000000000000000000000000000000
--- a/example/maxwell_sibc_lagrange/cimResolution.pro
+++ /dev/null
@@ -1,47 +0,0 @@
-/////////////////////////////////////////////////////////////////////
-// Resolution for contour integral method. //
-// Warning, the master file should include "cimParameters.pro", //
-// and should also define "FormulationCIM" the formulation to use. //
-/////////////////////////////////////////////////////////////////////
-Resolution{
- { Name Solve;
- System{ { Name A; NameOfFormulation FormulationCIM; Type ComplexValue; } }
- Operation{
- MPI_SetCommSelf;
- If(!doInit)
- Evaluate[Python[]{"real_corrector_init.py"}];
- EndIf
-
- Generate[A];
-
- If(doInit)
- CopySolution[A, x~{0}()];
- EndIf
-
- If(doSolve)
- // Full solve for first RHS
- CopyRightHandSide[b~{0}(), A];
- Solve[A];
- CopySolution[A, x~{0}()];
-
- // SolveAgain for other RHSs
- For i In {1:nRHS-1}
- CopyRightHandSide[b~{i}(), A];
- SolveAgain[A];
- CopySolution[A, x~{i}()];
- EndFor
- EndIf
-
- If(doApply)
- CopySolution[x~{0}(), A];
- Apply[A];
- CopySolution[A, x~{0}()];
- EndIf
-
- If(doPostpro)
- CopySolution[x~{0}(), A];
- PostOperation[Post];
- EndIf
- }
- }
-}
diff --git a/example/maxwell_sibc_lagrange/real_corrector_handl.py b/example/maxwell_sibc_lagrange/real_corrector_handl.py
deleted file mode 100644
index 269f94761999ff12809e131d88dfac5ab598f977..0000000000000000000000000000000000000000
--- a/example/maxwell_sibc_lagrange/real_corrector_handl.py
+++ /dev/null
@@ -1,22 +0,0 @@
-"""Real-part operator corrector
-
-Approximates the real-part operator by an holomorphic function, constructed with
-a Lagrange polynomial defined at some correction points. Please make sure that
-GetDP calls real_corrector_init.py before calling this file.
-
-The naming convention below must be followed:
--- a FixedPoint object named 'fixed' is available and prepared by cim.py
--- a Lagrange Interpolator object named 'lagrange' coming from initialization
--- the variable 'output' is read back by GetDP as the Pyhton[] call output
-"""
-
-# Input angular frequency
-omega = input[0]
-
-# Compute holomorphic approximation of the real part operator
-if(fixed.needNonHolomorphic()):
- output = np.real(omega) # If need non-holomorphic, do it
-elif(not fixed.cluster()):
- output = omega # If no cluster: Re(omega) ~ omega
-else:
- output = (omega + lagrange.at(omega)) / 2 # Else: correct
diff --git a/example/maxwell_sibc_lagrange/real_corrector_init.py b/example/maxwell_sibc_lagrange/real_corrector_init.py
deleted file mode 100644
index 942aa3b9bd16b11c4fc6a22250d18ece71eb4cb1..0000000000000000000000000000000000000000
--- a/example/maxwell_sibc_lagrange/real_corrector_init.py
+++ /dev/null
@@ -1,13 +0,0 @@
-"""Real-part operator corrector: initialization
-
-See: real_corrector_handl.py for more information
-"""
-
-from interpolator import Lagrange
-
-# Previous eigenvalue clusters
-cluster = fixed.cluster()
-
-# Initialize interpolator
-if(cluster):
- lagrange = Lagrange(cluster, np.conj(cluster).tolist())
diff --git a/example/maxwell_sibc_lagrange/sphere.dat b/example/maxwell_sibc_lagrange/sphere.dat
deleted file mode 100644
index db37ad8585b2209ae7880972c6fc1857206727ea..0000000000000000000000000000000000000000
--- a/example/maxwell_sibc_lagrange/sphere.dat
+++ /dev/null
@@ -1,22 +0,0 @@
-// Geometry //
-DefineConstant[R = { 100e-3, Name "Input/00Geometry/00Radius [m]" }];
-
-// Mesh //
-DefineConstant[md = { 10, Name "Input/01Mesh/00Density [Radius^-1]" }];
-cl = R/md;
-
-// Material //
-DefineConstant[isSC = { 0,
- Name "Input/02Material/00Is superconductor?",
- GmshOption "Reset",
- Choices {0 = "No", 1 = "Yes"} }];
-If(isSC == 0)
-DefineConstant[Sigma = { 1e00,
- Name "Input/02Material/01Conductivity [S m^-1]" }];
-EndIf
-If(isSC == 1)
-DefineConstant[Sigma = { 1e9,
- Name "Input/02Material/01Conductivity [S m^-1]" }];
-DefineConstant[Lambda = { 1e-8,
- Name "Input/02Material/02London depth [m]" }];
-EndIf
diff --git a/example/maxwell_sibc_lagrange/sphere.geo b/example/maxwell_sibc_lagrange/sphere.geo
deleted file mode 100644
index 5ad10a3b551aa48ff4a5da29ae48d2c965ec4cbc..0000000000000000000000000000000000000000
--- a/example/maxwell_sibc_lagrange/sphere.geo
+++ /dev/null
@@ -1,50 +0,0 @@
-Include "sphere.dat";
-
-Point(0) = { 0, 0, 0, cl};
-
-Point(1) = {+R, 0, 0, cl};
-Point(2) = { 0, +R, 0, cl};
-Point(3) = {-R, 0, 0, cl};
-Point(4) = { 0, -R, 0, cl};
-Point(5) = { 0, 0, +R, cl};
-Point(6) = { 0, 0, -R, cl};
-
-Circle(01) = {1, 0, 2};
-Circle(02) = {2, 0, 3};
-Circle(03) = {3, 0, 4};
-Circle(04) = {4, 0, 1};
-
-Circle(05) = {5, 0, 1};
-Circle(06) = {1, 0, 6};
-Circle(07) = {6, 0, 3};
-Circle(08) = {3, 0, 5};
-
-Circle(09) = {2, 0, 6};
-Circle(10) = {6, 0, 4};
-Circle(11) = {4, 0, 5};
-Circle(12) = {5, 0, 2};
-
-Line Loop(1) = {1, -12, 5};
-Line Loop(2) = {1, 9, -6};
-Line Loop(3) = {9, 7, -2};
-Line Loop(4) = {2, 8, 12};
-Line Loop(5) = {11, 5, -4};
-Line Loop(6) = {4, 6, 10};
-Line Loop(7) = {10, -3, -7};
-Line Loop(8) = {3, 11, -8};
-
-Ruled Surface(1) = {1};
-Ruled Surface(2) = {2};
-Ruled Surface(3) = {3};
-Ruled Surface(4) = {4};
-Ruled Surface(5) = {5};
-Ruled Surface(6) = {6};
-Ruled Surface(7) = {7};
-Ruled Surface(8) = {8};
-
-Surface Loop(1) = {3, 2, 1, 4, 8, 7, 6, 5};
-Volume(1) = {1};
-
-Physical Volume(1) = {1};
-Physical Surface(2) = {1, 2, 3, 4, 5, 6, 7, 8};
-Physical Point(3) = {0};
diff --git a/example/maxwell_sibc_lagrange/sphere.pro b/example/maxwell_sibc_lagrange/sphere.pro
deleted file mode 100644
index a7c28b0d7c3e16a58be5d9aff7929aec2e1c29ed..0000000000000000000000000000000000000000
--- a/example/maxwell_sibc_lagrange/sphere.pro
+++ /dev/null
@@ -1,107 +0,0 @@
-Include "sphere.dat";
-
-Group{
- Omega = Region[1];
- Bndry = Region[2];
-}
-
-Include "./cimParameters.pro"; // Functions for cim.py
-
-Function{
- // Imaginary Unit //
- J[] = Complex[0, 1];
-
- // Physical Constants //
- C0 = 299792458; // [m/s]
- Mu0 = 4*Pi*1e-7; // [H/m]
- Eps0 = 1/(Mu0 * C0^2); // [F/m]
-
- // Correction for real-part operator //
- If(doInit)
- PyPulsRe[] = 1;
- EndIf
- If(!doInit)
- PyPulsRe[] = Python[Puls[]]{"real_corrector_handl.py"};
- EndIf
-
- // Conductivity (coming from sphere.dat) //
- If(isSC == 1) // Superconductor
- Sigma[] = Sigma - J[] / (Mu0 * Lambda^2 * PyPulsRe[]]); // [S/m]
- EndIf
- If(isSC == 0) // Normal conductor
- Sigma[] = Sigma; // [S/m]
- EndIf
-
- // Leontovich SIBC (H-Formulation) //
- SL[] = Sqrt[(J[] * PyPulsRe[] * Mu0) / Sigma[]]; // H-Formulation
- SIbc[] = -1.0 / SL[]; // E-Formulation
-}
-
-Integration{
- { Name I1;
- Case{
- { Type Gauss;
- Case{
- { GeoElement Tetrahedron; NumberOfPoints 4; } // O1 * O1
- { GeoElement Triangle; NumberOfPoints 3; } // O1 * O1
- { GeoElement Line; NumberOfPoints 2; } // O1 * O1
- }
- }
- }
- }
-}
-
-Jacobian{
- { Name Jac;
- Case{
- { Region Omega; Jacobian Vol; }
- { Region Bndry; Jacobian Sur; }
- }
- }
-}
-
-FunctionSpace{
- { Name HCurl; Type Form1;
- BasisFunction{
- { Name sx; NameOfCoef ux; Function BF_Edge;
- Support Region[{Omega, Bndry}]; Entity EdgesOf[All]; }
- }
- }
-}
-
-Formulation{
- { Name FormulationCIM; Type FemEquation;
- Quantity{
- { Name e; Type Local; NameOfSpace HCurl; }
- }
- Equation{
- // Maxwell
- Galerkin{ [ 1/Mu0 * Dof{d e}, {d e}];
- In Omega; Jacobian Jac; Integration I1; }
- Galerkin{ [ -Puls[]^2 * Eps0 * Dof{ e}, { e}];
- In Omega; Jacobian Jac; Integration I1; }
-
- // IBC
- Galerkin{ [-J[] * Puls[] * SIbc[] * Dof{ e}, { e}];
- In Bndry; Jacobian Jac; Integration I1; }
- }
- }
-}
-
-Include "./cimResolution.pro"; // Resolution for cim.py
-
-PostProcessing{
- { Name Post; NameOfFormulation FormulationCIM;
- Quantity{
- { Name E; Value{ Term{ [{e}]; In Omega; Jacobian Jac; } } }
- }
- }
-}
-
-PostOperation{
- { Name Post; NameOfPostProcessing Post;
- Operation{
- Print[E, OnElementsOf Omega, File fileName];
- }
- }
-}