FreeFemHelmholtzEngine.py
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Mon, Oct 28, 16:43

FreeFemHelmholtzEngine.py
View Options

	#!/usr/bin/python

	import warnings
	import os
	import subprocess
	import numpy as np
	import scipy.sparse as scsp
	import utilities
	import FreeFemConversionTools as FFCT
	from RROMPyTypes import Np1D, Np2D, DictAny, Tuple, List

	PI = np.pi

	class DomainError(Exception):
	"""
	Domain error exception.

	Args:
	value: Human readable string describing the exception.

	Attributes:
	value: Human readable string describing the exception.
	"""
	def __init__(self, value:str):
	self.value = value
	def __str__(self):
	return self.value

	class FreeFemHelmholtzEngine:
	"""
	FreeFem++-based solver for Helmholtz problems.
	- \nabla \cdot (a \nabla u) - k^2 * n*2 u = f in \Omega
	u = u0 on \Gamma_D
	\partial_nu = g1 on \Gamma_N
	\partial_nu + h u = g2 on \Gamma_R

	Args:

	Attributes:
	"""

	def __init__(self, V:str,
	wavenumber:complex,
	compl : bool = True,
	dim : int = 2,
	meshname : str = "Th",
	Vname : str = "V",
	diffusivity : str = "1",
	refractionIndex : str = "1",
	forcingTerm : str = "0",
	DirichletBoundary : str = "",
	NeumannBoundary : str = "",
	RobinBoundary : str = "",
	DirichletDatum : str = "0",
	NeumannDatum : str = "0",
	RobinDatum : str = "0",
	RobinDatum2 : str = "0"):
	# process boundaries
	boundariesList = [DirichletBoundary, NeumannBoundary, RobinBoundary]
	for i in range(3):
	boundariesList[i] = boundariesList[i].strip()
	DirichletBoundary, NeumannBoundary, RobinBoundary = boundariesList
	if boundariesList.count("") == 3:
	raise DomainError("At least one boundary must be prescribed.")
	self.DirichletBoundary = DirichletBoundary
	self.NeumannBoundary = NeumannBoundary
	self.RobinBoundary = RobinBoundary

	if compl:
	self.compl = "<complex>"
	else:
	self.compl = ""
	self.V = V
	self.dim = dim
	self.meshname = meshname
	self.Vname = Vname
	self.wavenumber2 = wavenumber ** 2.
	self.diffusivity = diffusivity
	self.refractionIndex = refractionIndex
	self.forcingTerm = forcingTerm
	self.DirichletDatum = DirichletDatum
	self.NeumannDatum = NeumannDatum
	self.RobinDatum = RobinDatum
	self.RobinDatum2 = RobinDatum2

	functional = lambda self, u: 0.

	def Vdim(self) -> int:
	"""Dimension of finite element space."""
	if not hasattr(self, "liftDData"): self.liftDirichletData()
	return len(self.liftDData)

	def energyNormMatrix(self, w:float) -> Np2D:
	"""
	Sparse matrix (in CSR format) assoociated to w-weighted H10 norm.

	Args:
	w: Weight.

	Returns:
	Sparse matrix in CSR format.
	"""
	filename = utilities.getNewFilename("temp", "edp")
	Mfile = utilities.getNewFilename("m", "tmp")
	with open(filename, 'w') as f:
	f.write(self.V + "\n")
	f.write("{}{} u, v;\n".format(self.Vname, self.compl))
	f.write("varf energy(u, v) = int{}d({})({})".format(
	self.dim,
	self.meshname,
	FFCT.diff2(dim = self.dim,
	compl = not self.compl == "")))
	f.write(" + int{}d({})({} * {});\n".format(
	self.dim,
	self.meshname,
	w ** 2,
	FFCT.prod2(compl = not self.compl == "")))
	f.write("matrix{0} M = energy({1}, {1});\n".format(self.compl,
	self.Vname))
	f.write("ofstream Mout(\"{}\");\n".format(Mfile))
	f.write("Mout << M;\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 file {}.")\
	.format(output.decode('utf-8'), filename))
	M = FFCT.ffmat2spmat(Mfile)
	os.remove(filename)
	os.remove(Mfile)
	return M.tocsr()

	def problemData(self) -> DictAny:
	"""List of HF problem data."""
	dataDict = {}
	dataDict["forcingTerm"] = self.forcingTerm
	dataDict["DirichletDatum"] = self.DirichletDatum
	dataDict["NeumannDatum"] = self.NeumannDatum
	dataDict["RobinDatum2"] = self.RobinDatum2
	return dataDict

	def setProblemData(self, dataDict:DictAny, k2:complex):
	"""
	Set problem data.

	Args:
	dataDict: Dictionary of problem data.
	k2: Parameter value.
	"""
	self.wavenumber2 = k2
	self.forcingTerm = dataDict["forcingTerm"]
	self.DirichletDatum = dataDict["DirichletDatum"]
	self.NeumannDatum = dataDict["NeumannDatum"]
	self.RobinDatum2 = dataDict["RobinDatum2"]
	if hasattr(self, "A0"): del self.A0
	if hasattr(self, "A1"): del self.A1
	if hasattr(self, "b0"): del self.b0

	def liftDirichletData(self):
	"""Lift Dirichlet datum."""
	if not hasattr(self, "liftDData"):
	filename = utilities.getNewFilename("temp", "edp")
	u0file = utilities.getNewFilename("u0", "tmp")
	with open(filename, 'w') as f:
	f.write(self.V + "\n")
	f.write("{}{} u0 = {};\n".format(self.Vname, self.compl,
	self.DirichletDatum.replace('j', 'i')))
	f.write("ofstream u0out(\"{}\");\n".format(u0file))
	f.write("u0out << u0[];\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 file {}.")\
	.format(output.decode('utf-8'), filename))
	self.liftDData = FFCT.ffvec2npvec(u0file)
	os.remove(filename)
	os.remove(u0file)
	return self.liftDData

	@property
	def diffusivity(self):
	"""Value of a."""
	return self._diffusivity
	@diffusivity.setter
	def diffusivity(self, diffusivity):
	if hasattr(self, "A0"): del self.A0
	self._diffusivity = diffusivity

	@property
	def refractionIndex(self):
	"""Value of n."""
	return self._refractionIndex
	@refractionIndex.setter
	def refractionIndex(self, refractionIndex):
	if hasattr(self, "A1"): del self.A1
	self._refractionIndex = refractionIndex

	@property
	def forcingTerm(self):
	"""Value of f."""
	return self._forcingTerm
	@forcingTerm.setter
	def forcingTerm(self, forcingTerm):
	if hasattr(self, "b0"): del self.b0
	self._forcingTerm = forcingTerm

	@property
	def DirichletDatum(self):
	"""Value of u0."""
	return self._DirichletDatum
	@DirichletDatum.setter
	def DirichletDatum(self, DirichletDatum):
	if hasattr(self, "b0"): del self.b0
	if hasattr(self, "liftDData"): del self.liftDData
	self._DirichletDatum = DirichletDatum

	@property
	def NeumannDatum(self):
	"""Value of g1."""
	return self._NeumannDatum
	@NeumannDatum.setter
	def NeumannDatum(self, NeumannDatum):
	if hasattr(self, "b0"): del self.b0
	self._NeumannDatum = NeumannDatum

	@property
	def RobinDatum(self):
	"""Value of h."""
	return self._RobinDatum
	@RobinDatum.setter
	def RobinDatum(self, RobinDatum):
	if hasattr(self, "A0"): del self.A0
	self._RobinDatum = RobinDatum

	@property
	def RobinDatum2(self):
	"""Value of g2."""
	return self._RobinDatum2
	@RobinDatum2.setter
	def RobinDatum2(self, RobinDatum2):
	if hasattr(self, "b0"): del self.b0
	self._RobinDatum2 = RobinDatum2

	def assembleA(self):
	"""Assemble matrix blocks of linear system."""
	if not (hasattr(self, "A0") and hasattr(self, "A1")):
	filename = utilities.getNewFilename("temp", "edp")
	with open(filename, 'w') as f:
	f.write(self.V + "\n")
	f.write("{}{} u, v;\n".format(self.Vname, self.compl))
	if not hasattr(self, "A0"):
	A0file = utilities.getNewFilename("A0", "tmp")
	f.write("varf a0(u, v) = int{}d({})({} * {})".format(
	self.dim,
	self.meshname,
	self.diffusivity.replace('j', 'i'),
	FFCT.diff2(dim = self.dim,
	compl = not self.compl == "")))
	if len(self.RobinBoundary) > 0:
	f.write(" + int{}d({}, {})({} * {})".format(
	self.dim - 1,
	self.meshname,
	self.RobinBoundary,
	self.RobinDatum.replace('j', 'i'),
	FFCT.prod2(compl = not self.compl == "")))
	if len(self.DirichletBoundary) > 0:
	f.write(" + on({}, u = 0.);\n".format(
	self.DirichletBoundary))
	else:
	f.write(";\n")
	f.write("matrix{0} A0 = a0({1},{1});\n".format(self.compl,
	self.Vname))
	f.write("ofstream A0out(\"{}\");\n".format(A0file))
	f.write("A0out << A0;\n")
	if not hasattr(self, "A1"):
	A1file = utilities.getNewFilename("A1", "tmp")
	f.write("varf a1(u, v) = - int{}d({})(({})^2 * {})".format(
	self.dim,
	self.meshname,
	self.refractionIndex.replace('j', 'i'),
	FFCT.prod2(compl = not self.compl == "")))
	if len(self.DirichletBoundary) > 0:
	f.write(" + on({}, u = 0.);\n".format(
	self.DirichletBoundary))
	else:
	f.write(";\n")
	f.write("matrix{0} A1 = a1({1},{1}, tgv = 0.);\n".format(
	self.compl, self.Vname))
	f.write("ofstream A1out(\"{}\");\n".format(A1file))
	f.write("A1out << A1;\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 file {}.")\
	.format(output.decode('utf-8'), filename))
	if not hasattr(self, "A0"):
	self.A0 = FFCT.ffmat2spmat(A0file)
	os.remove(A0file)
	if not hasattr(self, "A1"):
	self.A1 = FFCT.ffmat2spmat(A1file)
	os.remove(A1file)
	os.remove(filename)

	def assembleb(self):
	"""Assemble RHS blocks of linear system."""
	if not hasattr(self, "b0"):
	filename = utilities.getNewFilename("temp", "edp")
	b0file = utilities.getNewFilename("b0", "tmp")
	if self.compl == "":
	complEff = "real"
	else:
	complEff = "complex"
	with open(filename, 'w') as f:
	f.write(self.V + "\n")
	f.write("{}{} u, v;\n".format(self.Vname, self.compl))
	f.write("varf b0(u, v) = int{}d({})({})".format(
	self.dim,
	self.meshname,
	FFCT.prod2(self.forcingTerm.replace('j', 'i'),
	compl = not self.compl == "")))
	if len(self.NeumannBoundary) > 0:
	f.write(" + int{}d({}, {})({})".format(
	self.dim - 1,
	self.meshname,
	self.NeumannBoundary,
	FFCT.prod2(self.NeumannDatum.replace('j', 'i'),
	compl = not self.compl == "")))
	if len(self.RobinBoundary) > 0:
	f.write(" + int{}d({}, {})({})".format(
	self.dim - 1,
	self.meshname,
	self.RobinBoundary,
	FFCT.prod2(self.RobinDatum2.replace('j', 'i'),
	compl = not self.compl == "")))
	if len(self.DirichletBoundary) > 0:
	f.write(" + on({}, u = {});\n".format(
	self.DirichletBoundary,
	self.DirichletDatum.replace('j', 'i')))
	else:
	f.write(";\n")
	f.write("{}[int] B0 = b0(0, {});\n".format(complEff,
	self.Vname))
	f.write("ofstream b0out(\"{}\");\n".format(b0file))
	f.write("b0out << B0;\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 file {}.")\
	.format(output.decode('utf-8'), filename))
	self.b0 = FFCT.ffvec2npvec(b0file)
	os.remove(b0file)
	os.remove(filename)

	parDegree = 1

	def As(self) -> List[Np2D]:
	"""Matrix blocks of linear system."""
	self.assembleA()
	return [self.A0, self.A1]

	def bs(self) -> List[Np1D]:
	"""RHS blocks of linear system."""
	self.assembleb()
	return [self.b0]

	def A(self) -> Np2D:
	"""Assemble matrix of linear system."""
	self.assembleA()
	return self.A0 + self.wavenumber2 * self.A1

	def b(self) -> Np1D:
	"""Assemble RHS of linear system."""
	self.assembleb()
	return self.b0

	def solve(self) -> Np1D:
	"""
	Find solution of linear system.

	Returns:
	Solution vector.
	"""
	return scsp.linalg.spsolve(self.A(), self.b())

	def getLSBlocks(self) -> Tuple[List[Np2D], List[Np1D]]:
	"""
	Get blocks of linear system.

	Returns:
	Blocks of LS (\sum_{j=0}^J k^j A_j = \sum_{j=0}^K k^j b_j).
	"""
	return self.As(), self.bs()


	class FreeFemHelmholtzScatteringEngine(FreeFemHelmholtzEngine):
	"""
	FreeFem++-based solver for Helmholtz scattering problems with Robin ABC.
	- \nabla \cdot (a \nabla u) - k^2 * n*2 u = f in \Omega
	u = u0 on \Gamma_D
	\partial_nu = g1 on \Gamma_N
	\partial_nu - i k u = g2 on \Gamma_R

	Args:

	Attributes:
	"""

	def __init__(self, V:str,
	wavenumber:complex,
	compl : bool = True,
	dim : int = 2,
	meshname : str = "Th",
	Vname : str = "V",
	diffusivity : str = "1",
	refractionIndex : str = "1",
	forcingTerm : str = "0",
	DirichletBoundary : str = "",
	NeumannBoundary : str = "",
	RobinBoundary : str = "",
	DirichletDatum : str = "0",
	NeumannDatum : str = "0",
	RobinDatum2 : str = "0"):
	self.wavenumber = wavenumber
	FreeFemHelmholtzEngine.__init__(self, V = V,
	wavenumber = wavenumber,
	compl = compl,
	dim = dim,
	meshname = meshname,
	Vname = Vname,
	diffusivity = diffusivity,
	refractionIndex = refractionIndex,
	forcingTerm = forcingTerm,
	DirichletBoundary = DirichletBoundary,
	NeumannBoundary = NeumannBoundary,
	RobinBoundary = RobinBoundary,
	DirichletDatum = DirichletDatum,
	NeumannDatum = NeumannDatum,
	RobinDatum = -1.j * wavenumber,
	RobinDatum2 = RobinDatum2)
	if len(self.RobinBoundary) == 0:
	raise Exception(("Empty Robin boundary not allowed in scattering "
	"engine. Use standard engine instead."))

	def setProblemData(self, dataDict:DictAny, k:complex):
	"""
	Set problem data.

	Args:
	dataDict: Dictionary of problem data.
	k: Parameter value.
	"""
	self.wavenumber = k
	self.forcingTerm = dataDict["forcingTerm"]
	self.DirichletDatum = dataDict["DirichletDatum"]
	self.NeumannDatum = dataDict["NeumannDatum"]
	self.RobinDatum2 = dataDict["RobinDatum2"]
	if hasattr(self, "A0"): del self.A0
	if hasattr(self, "A1"): del self.A1
	if hasattr(self, "A2"): del self.A2
	if hasattr(self, "b0"): del self.b0

	@property
	def wavenumber2(self):
	"""Value of k^2."""
	return self._wavenumber2
	@wavenumber2.setter
	def wavenumber2(self, wavenumber2):
	self._wavenumber2 = wavenumber2
	self._wavenumber = self.wavenumber2 ** .5
	self.RobinDatum = "{}".format(-1.j*self.wavenumber)

	@property
	def wavenumber(self):
	"""Value of k."""
	return self._wavenumber
	@wavenumber.setter
	def wavenumber(self, wavenumber):
	self.wavenumber2 = wavenumber ** 2.

	@property
	def RobinDatum(self):
	"""Value of h."""
	return super().RobinDatum
	@RobinDatum.setter
	def RobinDatum(self, RobinDatum):
	try:
	self._RobinDatum = "{}".format(complex(RobinDatum))
	except:
	raise Exception(("Value of RobinDatum not recognized. Must be "
	"complex number of str convertible to complex."))
	if not (hasattr(self, 'wavenumber')
	and np.isclose(self.wavenumber, 1.j * complex(self.RobinDatum))):
	self.wavenumber = 1.j * complex(self.RobinDatum)

	def assembleA(self):
	"""Assemble matrix blocks of linear system."""
	if not (hasattr(self, "A0") and hasattr(self, "A1")
	and hasattr(self, "A2")):
	filename = utilities.getNewFilename("temp", "edp")
	with open(filename, 'w') as f:
	f.write(self.V + "\n")
	f.write("{}{} u, v;\n".format(self.Vname, self.compl))
	if not hasattr(self, "A0"):
	A0file = utilities.getNewFilename("A0", "tmp")
	f.write("varf a0(u, v) = int{}d({})({} * {})".format(
	self.dim,
	self.meshname,
	self.diffusivity.replace('j', 'i'),
	FFCT.diff2(dim = self.dim,
	compl = not self.compl == "")))
	if len(self.DirichletBoundary) > 0:
	f.write(" + on({}, u = 0.);\n".format(
	self.DirichletBoundary))
	else:
	f.write(";\n")
	f.write("matrix{0} A0 = a0({1},{1});\n".format(self.compl,
	self.Vname))
	f.write("ofstream A0out(\"{}\");\n".format(A0file))
	f.write("A0out << A0;\n")
	if not hasattr(self, "A1"):
	A1file = utilities.getNewFilename("A1", "tmp")
	f.write("varf a1(u, v) = - int{}d({}, {})(1.i * {})".format(
	self.dim - 1,
	self.meshname,
	self.RobinBoundary,
	FFCT.prod2(compl = not self.compl == "")))
	if len(self.DirichletBoundary) > 0:
	f.write(" + on({}, u = 0.);\n".format(
	self.DirichletBoundary))
	else:
	f.write(";\n")
	f.write("matrix{0} A1 = a1({1},{1}, tgv = 0.);\n".format(
	self.compl, self.Vname))
	f.write("ofstream A1out(\"{}\");\n".format(A1file))
	f.write("A1out << A1;\n")
	if not hasattr(self, "A2"):
	A2file = utilities.getNewFilename("A2", "tmp")
	f.write("varf a2(u, v) = - int{}d({})({} * {})".format(
	self.dim,
	self.meshname,
	self.refractionIndex.replace('j', 'i'),
	FFCT.prod2(compl = not self.compl == "")))
	if len(self.DirichletBoundary) > 0:
	f.write(" + on({}, u = 0.);\n".format(
	self.DirichletBoundary))
	else:
	f.write(";\n")
	f.write("matrix{0} A2 = a2({1},{1}, tgv = 0.);\n".format(
	self.compl, self.Vname))
	f.write("ofstream A2out(\"{}\");\n".format(A2file))
	f.write("A2out << A2;\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 file {}.")\
	.format(output.decode('utf-8'), filename))
	if not hasattr(self, "A0"):
	self.A0 = FFCT.ffmat2spmat(A0file)
	os.remove(A0file)
	if not hasattr(self, "A1"):
	self.A1 = FFCT.ffmat2spmat(A1file)
	os.remove(A1file)
	if not hasattr(self, "A2"):
	self.A2 = FFCT.ffmat2spmat(A2file)
	os.remove(A2file)
	os.remove(filename)

	def As(self) -> List[Np2D]:
	"""Matrix blocks of linear system."""
	self.assembleA()
	return [self.A0, self.A1, self.A2]

	def A(self) -> Np2D:
	"""Assemble matrix of linear system."""
	self.assembleA()
	return self.A0 + self.wavenumber * self.A1 + self.wavenumber2 * self.A2


	class FreeFemHelmholtzScatteringAugmentedEngine(
	FreeFemHelmholtzScatteringEngine):
	"""
	FreeFem++-based solver for Helmholtz scattering problems with Robin ABC.
	- \nabla \cdot (a \nabla u) - k^2 * n*2 u = f in \Omega
	u = u0 on \Gamma_D
	\partial_nu = g1 on \Gamma_N
	\partial_nu - i k u = g2 on \Gamma_R
	Linear dependence on k is achieved by introducing an additional variable.

	Args:

	Attributes:
	"""

	def __init__(self, V:str,
	wavenumber:complex,
	compl : bool = True,
	dim : int = 2,
	meshname : str = "Th",
	Vname : str = "V",
	diffusivity : str = "1",
	refractionIndex : str = "1",
	forcingTerm : str = "0",
	DirichletBoundary : str = "",
	NeumannBoundary : str = "",
	RobinBoundary : str = "",
	DirichletDatum : str = "0",
	NeumannDatum : str = "0",
	RobinDatum2 : str = "0",
	constraintType : str = "IDENTITY"):
	self.constraintType = constraintType
	FreeFemHelmholtzScatteringEngine.__init__(self, V = V,
	wavenumber = wavenumber,
	compl = compl,
	dim = dim,
	meshname = meshname,
	Vname = Vname,
	diffusivity = diffusivity,
	refractionIndex = refractionIndex,
	forcingTerm = forcingTerm,
	DirichletBoundary = DirichletBoundary,
	NeumannBoundary = NeumannBoundary,
	RobinBoundary = RobinBoundary,
	DirichletDatum = DirichletDatum,
	NeumannDatum = NeumannDatum,
	RobinDatum2 = RobinDatum2)

	def energyNormMatrix(self, w:float) -> Np2D:
	"""
	Sparse matrix (in CSR format) assoociated to w-weighted H10 norm.

	Args:
	w: Weight.

	Returns:
	Sparse matrix in CSR format.
	"""
	normMatFen = FreeFemHelmholtzEngine.energyNormMatrix(self, w)
	return scsp.block_diag((normMatFen, 1/w * normMatFen))

	def liftDirichletData(self):
	"""Lift Dirichlet datum."""
	lift1 = FreeFemHelmholtzScatteringEngine.liftDirichletData(self)
	return np.pad(lift1, (0, len(lift1)), 'constant')

	def setProblemData(self, dataDict:DictAny, k:complex):
	"""
	Set problem data.

	Args:
	dataDict: Dictionary of problem data.
	k: Parameter value.
	"""
	self.wavenumber = k
	self.forcingTerm = dataDict["forcingTerm"]
	self.DirichletDatum = dataDict["DirichletDatum"]
	self.NeumannDatum = dataDict["NeumannDatum"]
	self.RobinDatum2 = dataDict["RobinDatum2"]
	if hasattr(self, "A0"): del self.A0
	if hasattr(self, "A1"): del self.A1
	if hasattr(self, "b0"): del self.b0

	@property
	def constraintType(self):
	"""Value of constraintType."""
	return self._constraintType
	@constraintType.setter
	def constraintType(self, constraintType):
	try:
	constraintType = constraintType.upper().strip().replace(" ","")
	if constraintType not in ["IDENTITY", "MASS"]: raise
	self._constraintType = constraintType
	except:
	warnings.warn(("Prescribed constraintType not recognized. "
	"Overriding to 'IDENTITY'."))
	self.constraintType = "IDENTITY"

	def assembleA(self):
	"""Assemble matrix blocks of linear system."""
	builtA0 = hasattr(self, "A0")
	builtA1 = hasattr(self, "A1")
	if not (builtA0 and builtA1):
	if builtA1 and self.constraintType == "IDENTITY":
	self.A2 = None
	FreeFemHelmholtzScatteringEngine.assembleA(self)
	if self.constraintType == "IDENTITY":
	if not builtA0:
	block = scsp.eye(self.A0.shape[0])
	else:
	block = scsp.eye(self.A1.shape[0])
	else: #MASS
	block = - self.A2
	I = list(self.DirichletBC0.get_boundary_values().keys())
	warnings.simplefilter("ignore", scsp.SparseEfficiencyWarning)
	block[I, I] = self.A0[I[0], I[0]]
	warnings.simplefilter("default", scsp.SparseEfficiencyWarning)
	if not builtA0:
	self.A0 = scsp.block_diag((self.A0, block))
	if not builtA1:
	self.A1 = scsp.bmat([[self.A1, self.A2], [- block, None]])
	if hasattr(self, 'A2'): del self.A2

	def assembleb(self):
	"""Assemble RHS of linear system."""
	if not hasattr(self, "b0"):
	FreeFemHelmholtzScatteringEngine.assembleb(self)
	self.b0 = np.pad(self.b0, (0, self.b0.shape[0]), 'constant')

	parDegree = 1

	def As(self) -> List[Np2D]:
	"""Matrix blocks of linear system."""
	self.assembleA()
	return [self.A0, self.A1]

	def A(self) -> Np2D:
	"""Assemble matrix of linear system."""
	self.assembleA()
	return self.A0 + self.wavenumber * self.A1

	def solve(self) -> Tuple[Np1D, Np1D]:
	"""
	Find solution of linear system.

	Returns:
	Solution vector.
	"""
	u = FreeFemHelmholtzScatteringEngine.solve(self)
	lu2 = int(len(u) / 2)
	return np.array([u[: lu2], u[lu2 :]])

FreeFemHelmholtzEngine.pyNo OneTemporaryActions

File Metadata

FreeFemHelmholtzEngine.pyView Options

Event Timeline

FreeFemHelmholtzEngine.py
No OneTemporary
Actions

FreeFemHelmholtzEngine.py
View Options