1"""Minimal ctypes interface to the fmm3dbie shared library.
3Only the routines needed for the interior Helmholtz impedance problem are
4wrapped. All Fortran arguments are passed by reference; integers are
5integer*8, and arrays are Fortran (column-major) ordered.
6"""
8import contextlib
9import ctypes
10import os
11import sys
13# OpenBLAS threads compete with the library's OpenMP threads
14os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
15# write Fortran console output immediately so it can be redirected to a log
16os.environ.setdefault("GFORTRAN_UNBUFFERED_PRECONNECTED", "y")
18import numpy as np
20LIBPATH = os.environ.get(
21 "FMM3DBIE_LIB", os.path.expanduser("~/src/fmm3dbie/lib/libfmm3dbie.so")
22)
23_lib = ctypes.CDLL(LIBPATH)
25_i8 = ctypes.c_int64
26_f8 = ctypes.c_double
30def fortran_output_to(path):
31 """Redirect the library's console output (Fortran unit 6, i.e. fd 1) to a file."""
32 sys.stdout.flush()
33 saved = os.dup(1)
34 with open(path, "a") as f:
35 os.dup2(f.fileno(), 1)
36 try:
37 yield
38 finally:
39 _libc.fflush(None)
40 os.dup2(saved, 1)
41 os.close(saved)
44_libc = ctypes.CDLL(None)
47def _ref_i8(v):
48 return ctypes.byref(_i8(int(v)))
51def _ref_f8(v):
52 return ctypes.byref(_f8(float(v)))
55def _ptr(a):
56 return a.ctypes.data_as(ctypes.c_void_p)
59def _farr(a, dtype):
60 return np.asfortranarray(a, dtype=dtype)
63class Surface:
64 """Discretized surface in the fmm3dbie format (all patches share one order and type)."""
66 def __init__(self, norder, iptype, srcvals):
67 self.norder = norder
68 self.iptype = iptype
69 self.srcvals = _farr(srcvals, np.float64) # (12, npts)
70 self.npts = self.srcvals.shape[1]
71 self.npols = npols_for(norder, iptype)
72 self.npatches = self.npts // self.npols
73 self.norders = np.full(self.npatches, norder, dtype=np.int64)
74 self.iptypes = np.full(self.npatches, iptype, dtype=np.int64)
75 self.ixyzs = (np.arange(self.npatches + 1, dtype=np.int64) * self.npols + 1)
76 self.srccoefs = surf_vals_to_coefs(self, self.srcvals[0:9, :])
77 self.wts = get_qwts(self)
79 def args(self):
80 return (_ref_i8(self.npatches), _ptr(self.norders), _ptr(self.ixyzs),
81 _ptr(self.iptypes), _ref_i8(self.npts), _ptr(self.srccoefs),
82 _ptr(self.srcvals))
85def npols_for(norder, iptype):
86 if iptype == 1:
87 return (norder + 1) * (norder + 2) // 2
88 return (norder + 1) ** 2
91def get_disc_nodes(norder, iptype):
92 npols = npols_for(norder, iptype)
93 uvs = np.zeros((2, npols), order="F")
94 _lib.get_disc_nodes_(_ref_i8(norder), _ref_i8(npols), _ref_i8(iptype), _ptr(uvs))
95 return uvs
98def surf_vals_to_coefs(surf, vals):
99 vals = _farr(vals, np.float64)
100 nd = vals.shape[0]
101 coefs = np.zeros_like(vals, order="F")
102 _lib.surf_vals_to_coefs_(_ref_i8(nd), _ref_i8(surf.npatches), _ptr(surf.norders),
103 _ptr(surf.ixyzs), _ptr(surf.iptypes), _ref_i8(surf.npts),
104 _ptr(vals), _ptr(coefs))
105 return coefs
108def get_qwts(surf):
109 wts = np.zeros(surf.npts)
110 _lib.get_qwts_(_ref_i8(surf.npatches), _ptr(surf.norders), _ptr(surf.ixyzs),
111 _ptr(surf.iptypes), _ref_i8(surf.npts), _ptr(surf.srcvals), _ptr(wts))
112 return wts
115def helm_rpcomb_imp_solver(surf, zk, alpha, zlams, rhs, ifinout, eps=1e-6,
116 eps_gmres=1e-6, numit=400):
117 """Solve n.grad(u) + i k lambda u = rhs with u = S_k[s] + i alpha D_k S_{i|k|}[s].
119 ifinout = 0 for the interior problem, 1 for the exterior problem.
120 Returns (sigma, siksigma, niter, errs, rres).
121 """
122 zpars = np.array([zk, alpha], dtype=np.complex128)
123 zlams = _farr(zlams, np.complex128)
124 rhs = _farr(rhs, np.complex128)
125 sigma = np.zeros(surf.npts, dtype=np.complex128)
126 siksigma = np.zeros(surf.npts, dtype=np.complex128)
127 errs = np.zeros(numit + 1)
128 niter = _i8(0)
129 rres = _f8(0.0)
130 _lib.helm_rpcomb_imp_solver_(*surf.args(), _ref_f8(eps), _ptr(zpars), _ptr(zlams),
131 _ref_i8(numit), _ref_i8(ifinout), _ptr(rhs),
132 _ref_f8(eps_gmres), ctypes.byref(niter), _ptr(errs),
133 ctypes.byref(rres), _ptr(sigma), _ptr(siksigma))
134 return sigma, siksigma, niter.value, errs[: niter.value], rres.value
137def helm_rpcomb_eval(surf, zk, alpha, sigma, siksigma, targs, eps=1e-6):
138 """Evaluate u = S_k[s] + i alpha D_k S_{i|k|}[s] at off-surface targets (3, ntarg)."""
139 targs = _farr(targs, np.float64)
140 ntarg = targs.shape[1]
141 zpars = np.array([zk, alpha], dtype=np.complex128)
142 ipatch_id = np.full(ntarg, -1, dtype=np.int64)
143 uvs_targ = np.zeros((2, ntarg), order="F")
144 pot = np.zeros(ntarg, dtype=np.complex128)
145 _lib.helm_rpcomb_eval_(*surf.args(), _ref_i8(3), _ref_i8(ntarg), _ptr(targs),
146 _ptr(ipatch_id), _ptr(uvs_targ), _ref_f8(eps), _ptr(zpars),
147 _ptr(_farr(sigma, np.complex128)),
148 _ptr(_farr(siksigma, np.complex128)), _ptr(pot))
149 return pot
152def helm_rpcomb_eval_boundary(surf, zk, alpha, sigma, siksigma, eps=1e-6):
153 """Interior boundary limit of u = S_k[s] + i alpha D_k S_{i|k|}[s] at the nodes.
155 The on-surface evaluation gives the principal value of D_k; with outward
156 normals the interior limit adds -1/2 times the double-layer density.
157 """
158 uvs = get_disc_nodes(surf.norder, surf.iptype)
159 ipatch_id = np.repeat(np.arange(1, surf.npatches + 1), surf.npols).astype(np.int64)
160 uvs_targ = np.asfortranarray(np.tile(uvs, (1, surf.npatches)))
161 zpars = np.array([zk, alpha], dtype=np.complex128)
162 sigma = _farr(sigma, np.complex128)
163 siksigma = _farr(siksigma, np.complex128)
164 pot = np.zeros(surf.npts, dtype=np.complex128)
165 _lib.helm_rpcomb_eval_(*surf.args(), _ref_i8(12), _ref_i8(surf.npts), _ptr(surf.srcvals),
166 _ptr(ipatch_id), _ptr(uvs_targ), _ref_f8(eps), _ptr(zpars),
167 _ptr(sigma), _ptr(siksigma), _ptr(pot))
168 return pot - 0.5j * alpha * siksigma