ahbarnett / mpspack Goto Github PK

What steps will reproduce the problem?
Set the value of 'ob' to 't' in testqpscatt_fmm.m and run it.


What is the expected output? What do you see instead?
Error using iterapp (line 59)
user supplied function ==> @(co)pr.applybcmatrixFMM(co,o) failed with the 
following error:

 Matrix dimensions must agree.

Error in gmres (line 213)
r = b - iterapp('mtimes',afun,atype,afcnstr,x,varargin{:});

Error in bvp/linsolve (line 99)
        [co] = gmres(A, pr.rhs, o.restart, o.eps, maxit);

Error in bvp/solvecoeffs (line 138)
      co = pr.linsolve(o); % do the linear solve (o controls which type)

Error in qpscatt/solvecoeffs (line 463)
      co = solvecoeffs@bvp(pr,o);      % pass to usual solver

Error in testqpscatt_fmm (line 48)
tic; p.solvecoeffs(o); toc,      % (fill and) least-squares soln

What version of the product are you using? On what operating system?
Matlab2014b on MacOS 10.9.5 with installed and working fmmlib2d-1.2 and lp2d 
tools.

Please provide any additional information below.

In the attached test code recursivebessel.m returns a matrix containing NaN
entries on a 64 Bit Linux machine under Matlab 2009a

Since 'fd' is the premier method for low eigenvalues, should identify 
degeneracies closer than sqrt(eps) (?) and handle separately using local search 
('fd' or 'ms'?)
to try and regain accuracy close to eps.

Currently accuracy limited to sqrt(eps), ie 1e-8, by the Boyd analytic 
rootfinder.

Body and appendix of manual don't say much about bvp, scattering, and qpscatt 
classes and how to use them.

Since layerpot.eval with Jfilter scheme seems to demand high relative
accuracy on J Bessels near the origin, which recursivebessel doesn't yet
provide, we should extend regfb.usegsl to allow 3 ways to eval J Bessels:
Matlab native (slow), recursivebessel (fast, but not good rel error near
origin), GSL (fast, good rel error, but may not be installed).
If GSL is selected, but not installed, it should downgrade to Matlab native.

Code below Basis Functions:
Definition of Zp is wrong.

Running examples/twodieldiscspointsrcscat.m but with the point source lying 
inside a dielectric fails because u_inc is only nonzero in the airdom(s). Need 
to check for this and change u_inc to exist in non-air domains. A crude fix is 
to always make u_inc exist in all domains.

that's it. Have them check the fast property of the basis first.

The Laplace cases need to be added here, so Alpert can work.

What steps will reproduce the problem?
Run showsolution with a qpscatt problem, with >1e4 points.
Whenever problem.pointsolution decides to use blocking (Timo's code),
something interacts badly with qpscatt, and it crashes. 

What is the expected output? What do you see instead?
Crash. Expect plot.

Currently S,D,T kernel evaluations in @layerpot use the old fundsol and
fundsol_deriv codes (when opts.fast=0). They should be phased out in favor
of @utils/fundsol. The complication will be that the old fundsol_deriv
included computing (or having passed in) the cosine-phi factors, and this
will have to be explicitly coded since the new fundsol only computes radial
derivative.
This is largely moot for opts.fast=1,2 since then Rokhlin routines are
called directly (however fundsol_deriv is still used to compute cosine-phi,
clumsily).

MEX files eg for GSL bessel, hank103, etc, end with .mexglx for 32bit
intel, but .mexa64 for 64bit AMD, etc. Currently only a test is done for the
.mexglx suffix. Need a better way to do this ... in layerpot.m and mfsbasis.m

see summary.
Currently crashes if ask for [A Ax Ay] = b.eval(p) where p is pointset without
normals. Need not happen.

take the single airdom and use its scattered field rep to compute far field at 
a given single angle, or list of angles.
Make a user-friendly scattering.showfarfield which automates the call and plots 
the result, on polar diagram.
Make a scattering.totalcrosssection which refines angle sampling until converged
to a given tol.

Test it on triangle eigenvalues, and square scattering.

Ultimately will be replaced by a compressed corner rep from fast direct solver?

What steps will reproduce the problem?
1. run examples/tut_square.m

What is the expected output? What do you see instead?

Warning: Direct access of structure fields returned by a function call (e.g.,
 call to tut_square) is not allowed. See MATLAB 7.10 Release Notes, "Subscripting Into Function Return Values" for details. 
??? Attempt to reference field of non-structure array.

What version of the product are you using? On what operating system?
Matlab 2010a    linux.

Please provide any additional information below.

Currently segment.addinoutlayerpots only works for TM matching.
Different weights of layerpots in and out would be needed for TE or
other matching conditions - need to work this out. This may not be possible in 
the current segment-domain-basis framework - a new basis property such as 
domsweights may have to be added (analogous to doms)
Maybe have addinoutlayerpots check if matching conds exist and flag an error if 
not.

regfbbasis outputs NaN at the origin.

There are also some other minor things on in/out format to discuss.
I'd advocate for -N...N ordering since I use that in fast FB basis.
Also I prefer [A Ax Ay] as opposed to [A An Ax Ay], since if you request
Ax and Ay you don't also need An, in my experience. This saves memory.
Happy to discuss...

I never coded the cot hypersingular term in the Maue-Kress spectral
quadrature for the T boundary integral operator on closed analytic curves.
This wasn't needed for dielectric interfaces since in the Rokhlin scheme
the hypersingular terms cancel anyway. However, if we are to use D +
i.eta.S as an exterior representation coupled in value and normal
derivative to the subdomains, this will be needed.

recursivebessel should some day have better error checking,
and report if it returns inaccurate values, either absolute or
relative error. It is not designed to have good relative error
when values are very small anyway.

need to fix up special case of evaluation at zero radius from the origin,
as was done with regfbbasis. Important since quadrature pts may fall
at origin for quadr.clencurt

s = segment.smoothstar(900, 0.3, 5);
d = domain(s, 1); s.setbc(-1, 'D'); p = evp(d);
o.eps = 0.1; o.modes = 1; o.khat = 'r'; o.fhat = 's';
p.solvespectrum([1000 1000.1], 'ntd', o);
o = []; o.inds = 1; o.dx = 0.002; o.fmm = 1; o.col = 'bw';showmodes(p,o);

This is the code from Alex Barnett and Andrew Hassel's article "Fast 
Computation of the Dirichlet eigenmodes via the spectral flow of the interior 
Neumann-to-Dirichlet map".

I get the following error in Matlab


What version of the product are you using? On what operating system?


Please provide any additional information below.

??? Method 'hfmm2dparttarg' is not defined for class 'utils' or is removed from
MATLAB's search path.

Error in ==> layerpot.layerpot>layerpot.evalFMM at 338
        U = utils.hfmm2dparttarg(iprec,k,N,node,ifslp,b.a(1)*charge,...

Error in ==> problem.problem>problem.pointsolution at 334
              u(ii) = u(ii) + b.evalFMM(co, p.x(ii), opts); % opts needed, dom!

Error in ==> problem.problem>problem.gridsolution at 376
      [u di] = pr.pointsolution(pointset(zz(:)), o); % make zz a col vec

Error in ==> evp.evp>evp.showmodes at 571
        [u gx gy di] = pe.gridsolution(o); % either conventional or GRF evalu

Error in ==> trifoil_excited_state_mspack at 5
o = []; o.inds = 1; o.dx = 0.002; o.fmm = 1; o.col = 'bw';showmodes(p,o);

We need to be able to make copies of domains. Also I do this internally in
the domain constructor, when checking for intersections of various pieces.
So currently in my version it's a value class. (Segment is still handle
of course).

What do you think?

No way to supress the  output figure of [uj,x,y]=p.showmode()
I want to take only the data. In addition to that the matrix uj is not a square 
matrix. How can I specify my own x,y grid values ?

Also, investigate why help doesn't find class methods.

I cannot run test/testbasis.m without errors, because the file 
greengardrokhlinhank106 is missing.

Possible solutions:

a) Replace in testbasis.m in line 50 and 51 'o' with 'opt'.  <- actually a bug 
i guess
b) Commit the file hank106.f, which is needed to compile 
greengardrokhlinhank106.

"Diameter" is a misleading name for what is actually the radius of a domain,
so either change the name or make it return twice the current answer.

With the intention of getting used to the package, i decided to work myself 
through the tutorial. Unfortunately, i failed right at the start, i.e. at 
defining the domain for the introductory example (laplace's equation) on page 2.

While defining the segment "s" worked quite well (looked like figure 1), i got 
an error message when trying to put d=domain(s, +1):

"??? Attempt to execute SCRIPT utils.inpolyc as a function:
C:\Users\glatz\Documents\MATLAB\mpspack-1.1beta\@utils\inpolyc.m

Error in ==> utils.inpolywrapper at 31
i = logical(utils.inpolyc(p, v));   % converts int to logical

Error in ==> domain.domain>domain.inside at 114
          i = utils.inpolywrapper(p(:), v);

Error in ==> domain.domain>domain.normalscheck at 97
        norout = ~d.inside(p);                 % true if points away from domain

Error in ==> domain.domain>domain.domain at 75
        norout = normalscheck(d);              % normals sense of whole seg list"

'ms' mode: if degeneracy detected via multiple small singvals in gridminfit, 
should also make solvemodescoeffs handle this first in iter=0 mode (easy: SVD), 
and then in iter=1 mode (need to invent block version of iterative 
utils.minsingval).

Currently epwbasis interprets non-unity direction vector as giving the 
evanescent parameter, and chooses a particular decay direction (upwards).
This is a hack to get qpscatt polefix going.

Merge epwbasis and rpwbasis to become pwbasis, which can handle propagating and 
evanescent waves in all directions.

We should make a platform-independent Makefile
that checks if GSL is installed. I don't know how to do this.
It would make our potential users happier.

I've downloaded mpspack but I can't get the examples to work:


>> addpath /home/angus/mpspack % seems to work okay
Warning: Function /home/angus/mpspack/plot.m has the same name as a MATLAB 
builtin. We suggest you rename the function to avoid a potential name conflict. 
> In path at 110
  In addpath at 87

>> s = segment.radialfunc(250, {@(q) 1 + 0.3*cos(3*q), @(q) -0.9*sin(3*q)});
Error using segment.radialfunc
Method 'radialfunc' is not defined for class 'segment' oris removed from 
MATLAB's search path.


The radialfunc.m file is there, it's just not being accessed for some reason.

Matlab Version 8.1.0.604 (R2013a) on 64bit Arch Linux.

hank106.f source file is referenced in the @utils Makefile, but it is not in 
the repository.

in evp make a method for evaluating phi_j on the bdry then using its L2 norm to 
return absolute error in k eigenfrequency. Only do for Dirichlet BC for now.

In class domain constructor,
currently the topology-checking is only that connected segment pieces don't
fall outside the primary domain. We should eventually check for
self-intersections of each analytic segment, and intersections with other
segments in the same list.

Docs need:
review of various options, eg bb, dx in all the methods calling 
problem.gridsolution

EVP examples in Tutorial.

A multipole H-expansion basis object should be included for comparison
against MFS and layer-potential methods to represent outgoing fields.
For k=0 this would become a Laurent expansion