Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F90034112
FreeFemHelmholtzEngine.py
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Mon, Oct 28, 16:43
Size
30 KB
Mime Type
text/x-python
Expires
Wed, Oct 30, 16:43 (2 d)
Engine
blob
Format
Raw Data
Handle
21997019
Attached To
R6746 RationalROMPy
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
{}
\n
See 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
{}
\n
See 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
{}
\n
See 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
{}
\n
See 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
{}
\n
See 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
:]])
Event Timeline
Log In to Comment