Page MenuHomec4science
Files F85444462

FreeFemHSEngine.py
No OneTemporary
Actions

File Metadata

Created
Sun, Sep 29, 06:52

FreeFemHSEngine.py
View Options

#!/usr/bin/python
import os
import subprocess
import numpy as np
import warnings
import utilities
import FreeFemConversionTools as FFCT
from RROMPyTypes import Np1D, strLst, N2FSExpr
class FreeFemHSEngine:
"""
FreeFem++-based Hilbert space engine.
"""
def __init__(self, V:str, compl : bool = True, dim : int = 2,
meshname : str = "Th", Vname : str = "V"):
if compl:
self.compl = "<complex>"
else:
self.compl = ""
self.V = V
self.dim = dim
self.meshname = meshname
self.Vname = Vname
def name(self) -> str:
"""Class label."""
return self.__class__.__name__
def norm(self, u:Np1D, normType : N2FSExpr = "H10") -> float:
"""
Compute general norm of complex-valued function with given dofs.
Args:
u: numpy complex array with function dofs.
normType(optional): Target norm identifier. If matrix, target norm
is one induced by normType. If number, target norm is weighted
H^1 norm with given weight. If string, must be among 'H1',
'L2', and 'H10'. Defaults to 'H10'.
Returns:
Norm of the function (non-negative).
"""
if type(normType).__name__[-6:] == "matrix":
return np.abs(u.dot(normType.dot(u).conj())) ** .5
filename = utilities.getNewFilename("temp", "edp")
ufile = FFCT.npvec2ffvec(u)
if normType == 'L2':
coeff1 = 0.
coeff2 = 1.
else:
coeff1 = 1.
if isinstance(normType, (int, float)):
coeff2 = normType ** 2.
elif normType == 'H10':
coeff2 = 0.
elif normType == 'H1':
coeff2 = 1.
else:
raise Exception("Norm type not recognized.")
with open(filename, 'w') as f:
f.write(self.V + "\n")
f.write("ifstream fin(\"{}\");\n".format(ufile))
f.write("{}{} u;\n".format(self.Vname, self.compl))
f.write("fin >> u[];\n")
f.write("cout.precision(15);")
f.write("cout.scientific << sqrt(")
if np.isclose(coeff1, 0.):
f.write("0")
else:
f.write("{} * int{}d({})({})".format(coeff1, self.dim,
self.meshname,
FFCT.diff("u", self.dim)))
if np.isclose(coeff2, 0.):
f.write(" + 0);\n")
else:
f.write(" + {} * int{}d({})(abs(u)^2));\n".format(
coeff2,
self.dim,
self.meshname))
bashCommand = "FreeFem++ -v 0 {}".format(filename)
process = subprocess.Popen(bashCommand.split(), stdout=subprocess.PIPE)
output, error = process.communicate()
try:
nrm = np.float(output)
except:
raise Exception(("Error was encountered in the execution of "
"FreeFem++ code:\n{}\nSee files {} and {}.")\
.format(output.decode('utf-8'), filename, ufile))
os.remove(filename)
os.remove(ufile)
return nrm
def plot(self, u:Np1D, name : str = "u", save : bool = False,
what : strLst = "ALL", figspecs : str = "fill = 1, value = 1"):
"""
Do some nice plots of the complex-valued function with given dofs.
Args:
u: numpy complex array with function dofs.
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
save(optional): Whether to save plot(s). Defaults to False.
figspecs(optional): Optional arguments for FreeFem++ figure
creation. Defaults to "fill = 1, value = 1".
"""
if isinstance(what, (str,)):
what = what.upper()
if what != 'ALL':
what = [what]
elif save:
warnings.warn(("Saving and 'ALL' are not compatible. "
"Overriding 'ALL' to ['ABS', 'PHASE', 'REAL', "
"'IMAG']."))
what = ['ABS', 'PHASE', 'REAL', 'IMAG']
if what != 'ALL':
what = utilities.purgeList(what, ['ABS', 'PHASE', 'REAL', 'IMAG'],
listname = self.name() + ".what")
if len(what) == 0: return
filename = utilities.getNewFilename("temp", "edp")
ufile = FFCT.npvec2ffvec(u)
figspecs = figspecs.strip()
if len(figspecs) > 0 and figspecs[0] != ",":
figspecs = ", " + figspecs
with open(filename, 'w') as f:
f.write(self.V + "\n")
f.write("ifstream fin(\"{}\");\n".format(ufile))
f.write("{}{} u;\n".format(self.Vname, self.compl))
f.write("fin >> u[];\n")
if isinstance(what, (str,)):
f.write("plot(u, wait = 1, cmm=\"{}\"{});\n".format(
name, figspecs))
else:
for i, w in enumerate(what):
if i == 0:
VnameEff = self.Vname
else:
VnameEff = ""
if w == 'ABS':
outspecs = ""
if save:
outspecs = (", ps = \""
+ utilities.getNewFilename("fig", "abs.eps")
+ "\"")
f.write("{} v = abs(u);\n".format(VnameEff))
f.write("plot(v,wait = 1,cmm=\"|{}|\"{}{});\n".format(
name, figspecs, outspecs))
if w == 'PHASE':
outspecs = ""
if save:
outspecs = (", ps = \""
+ utilities.getNewFilename("fig", "ph.eps")
+ "\"")
f.write("{} v = arg(u);\n".format(VnameEff))
f.write("plot(v,wait = 1,cmm=\"ph {}\"{}{});\n".format(
name, figspecs, outspecs))
if w == 'REAL':
outspecs = ""
if save:
outspecs = (", ps = \""
+ utilities.getNewFilename("fig", "re.eps")
+ "\"")
f.write("{} v = real(u);\n".format(VnameEff))
f.write("plot(v,wait = 1,cmm=\"Re {}\"{}{});\n".format(
name, figspecs, outspecs))
if w == 'IMAG':
outspecs = ""
if save:
outspecs = (", ps = \""
+ utilities.getNewFilename("fig", "im.eps")
+ "\"")
f.write("{} v = imag(u);\n".format(VnameEff))
f.write("plot(v,wait = 1,cmm=\"Im {}\"{}{});\n".format(
name, figspecs, outspecs))
bashCommand = "FreeFem++ -v 0 {}".format(filename)
process = subprocess.Popen(bashCommand.split(), stdout=subprocess.PIPE)
output, error = process.communicate()
if len(output) > 0:
raise Exception(("Error was encountered in the execution of "
"FreeFem++ code:\n{}\nSee files {} and {}.")\
.format(output.decode('utf-8'), filename, ufile))
os.remove(filename)
os.remove(ufile)
def plotmesh(self, name : str = "Mesh", save : bool = False,
figspecs : str = ""):
"""
Do a nice plot of the mesh.
Args:
name(optional): Name to be shown as title of the plots. Defaults to
'Mesh'.
save(optional): Whether to save plot(s). Defaults to False.
figspecs(optional): Optional arguments for FreeFem++ figure
creation.
"""
filename = utilities.getNewFilename("temp", "edp")
figspecs = figspecs.strip()
if len(figspecs) > 0 and figspecs[0] != ",":
figspecs = ", " + figspecs
outspecs = ""
if save:
outspecs = (", ps = \"" + utilities.getNewFilename("msh", "eps")
+ "\"")
with open(filename, 'w') as f:
f.write(self.V + "\n")
f.write("plot({}, wait = 1, cmm=\"{}\"{}{});\n".format(
self.meshname, name, figspecs, outspecs))
bashCommand = "FreeFem++ -v 0 {}".format(filename)
process = subprocess.Popen(bashCommand.split(), stdout=subprocess.PIPE)
output, error = process.communicate()
if len(output) > 0:
raise Exception(("Error was encountered in the execution of "
"FreeFem++ code:\n{}\nSee file {}.")\
.format(output.decode('utf-8'), filename))
os.remove(filename)
def convExpression(self, expr:str) -> Np1D:
"""
Evaluate dofs of general expression.
Args:
expr: expression.
Returns:
Numpy 1D array with dofs.
"""
filename = utilities.getNewFilename("temp", "edp")
efile = utilities.getNewFilename("e", "tmp")
with open(filename, 'w') as f:
f.write(self.V + "\n")
f.write("{}{} expr = {};\n".format(self.Vname, self.compl,
expr.replace("j", "i")))
f.write("ofstream eout(\"{}\");\n".format(efile))
f.write("eout << expr[];\n")
bashCommand = "FreeFem++ -v 0 {}".format(filename)
process = subprocess.Popen(bashCommand.split(), stdout=subprocess.PIPE)
output, error = process.communicate()
if len(output) > 0:
raise Exception(("Error was encountered in the execution of "
"FreeFem++ code:\n{}\nSee files {} and {}.")\
.format(output.decode('utf-8'), filename, efile))
out = FFCT.ffvec2npvec(efile)
os.remove(filename)
os.remove(efile)
return out
class FreeFemHSAugmentedEngine(FreeFemHSEngine):
"""
FreeFem-based Hilbert space engine for augmented spaces.
"""
def __init__(self, V:str, d : int = 2, compl : bool = True, dim : int = 2,
meshname : str = "Th", Vname : str = "V"):
self.d = d
FreeFemHSEngine.__init__(self, V, compl, dim, meshname, Vname)
def norm(self, u:Np1D, normType : N2FSExpr = "H10") -> Np1D:
"""
Compute general norm of complex-valued function with given dofs.
Args:
u: numpy complex array with function dofs.
normType(optional): Target norm identifier. If matrix, target norm
is one induced by normType. If number, target norm is weighted
H^1 norm with given weight. If string, must be among 'H1',
'L2', and 'H10'. Defaults to 'H10'.
Returns:
Norms of the function (non-negative).
"""
if isinstance(normType, (int, float)):
u = utilities.listify(u, self.d)
norms = [None] * self.d
for j in range(self.d):
norms[j] = FreeFemHSEngine.norm(self, u[j], normType)
return norms
elif type(normType).__name__[-6:] == "matrix":
u = utilities.vectorify(u, self.d)
if normType.shape[0] % u[0] == 0:
N = int(normType.shape[0] / u[0])
else:
raise Exception("Energy matrix dimension mismatch.")
return FreeFemHSEngine.norm(self, np.concatenate(u[: N]), normType)
raise Exception("Norm type not recognized.")
def plot(self, u:Np1D, split : bool = True, name : str = "u",
save : bool = False, what : strLst = "ALL",
figspecs : str = "fill = 1, value = 1"):
"""
Do some nice plots of the complex-valued function with given dofs.
Args:
u: numpy complex array with function dofs.
split: whether to split the given array.
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
save(optional): Whether to save plot(s). Defaults to False.
what(optional): Which plots to do. If list, can contain 'ABS',
'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'.
Defaults to 'ALL'.
figspecs(optional): Optional arguments for FreeFem++ figure
creation. Defaults to "fill = 1, value = 1".
"""
if split:
u = utilities.listify(u, self.d)
for j in range(self.d):
FreeFemHSEngine.plot(self, u[j], what = what, save = save,
name = "{0}, comp. {1}".format(name, j),
figspecs = figspecs)
else:
FreeFemHSEngine.plot(self, u, what = what, save = save,
name = name, figspecs = figspecs)

Event Timeline

c4science · Help