The use of high-order basis functions in the forward solver

[ndgird]

Hello Toast++ team,

I would like to use high-order basis functions in the forward solver. Unfortunately, I have some troubles with that. My code with explanations is provided below:

1) I am using Matlab 2016a PDE Toolbox to generate mesh:

g = [    4.0000    4.0000    4.0000    4.0000
  -25.0000    0.0000   25.0000    0.0000
    0.0000   25.0000    0.0000  -25.0000
   -0.0000  -25.0000         0   25.0000
  -25.0000         0   25.0000   -0.0000
    1.0000    1.0000    1.0000    1.0000
         0         0         0         0
         0         0         0         0
         0         0         0         0
   25.0000   25.0000   25.0000   25.0000
   25.0000   25.0000   25.0000   25.0000
         0         0         0         0];
model = createpde(1);
pg = geometryFromEdges(model,g);
msh = generateMesh(model,'Hmax',6,'GeometricOrder','quadratic');
figure; pdeplot(model);

2) The generated mesh is passed into Toast++ “toastMesh” function along with quadratic basis functions type:

vtx = msh.Nodes';
idx = msh.Elements';
eltp = ones(size(idx,1),1)*6;
mesh = toastMesh(vtx,idx,eltp);
figure; mesh.Display;

3) I set up the scattering and absorption coefficients as follows:

mua0 = 0.01;
mus0 = 1.0;
Musfun0 = [num2str(mus0) '*cos(0.15.*x).*cos(0.15.*y)/2 + ' num2str(mus0)];
Muafun0 = [num2str(mua0) '*sin(0.15.*x).*sin(0.15.*y)/2 + ' num2str(mua0)];
Mua0 = ['@(x,y)' Muafun0];
Mus0 = ['@(x,y)' Musfun0];
MuaFunc = str2func(Mua0);
MusFunc = str2func(Mus0);
ref_bkg = 1.0;
nnd = mesh.NodeCount;
ref = ones(nnd,1) * ref_bkg;
N = 1024;
[gx1,gx2] = ndgrid(linspace(-25,25,N),linspace(-25,25,N));
MUA= MuaFunc(gx1,gx2);
MUS = MusFunc(gx1,gx2);
grd = size(MUA);
bb = [-25 25
      -25 25];

4) I call “toastBasis” function:
basis = toastBasis(mesh,grd,grd,bb,'LINEAR');

and get warnings: “Node without pixel support found”. What does this message mean? What is wrong?

5) Finally, I set up sources and detectors

nq = 16;
rad = 25;  %mesh radius [mm]
for i=1:nq
  phi_q = 2*pi*(i-1)/nq;
  Q(i,:) = rad * [cos(phi_q) sin(phi_q)];
  Q(i,:) = [24 0];
  phi_m = 2*pi*(i-0.5)/nq;
  M(i,:) = rad * [cos(phi_m) sin(phi_m)];
end
mesh.SetQM(Q,M);
hold on
plot(Q(:,1),Q(:,2),'ro','MarkerFaceColor','r');
plot(M(:,1),M(:,2),'bs','MarkerFaceColor','b');

% Create the source and boundary projection vectors
qvec = mesh.Qvec('Isotropic', 'Gaussian', 4);
mvec = mesh.Mvec('Gaussian', 2);

% Solve the FEM linear system
K = dotSysmat (mesh,mua,mus,ref,0);

and end up with a system matrix K which contains “inf” and “-inf” values.

What is wrong with my code? Could you help me to understand how to use high-order basis functions in Toast++?

This is an info about my setup:
Mac OS X Yosemite 10.10.5
Matlab R2016a
Toast++ 2.0.0 with precompiled binaries downloaded from https://sourceforge.net/projects/toastpp/files/

Thank you very much for your help.

[Martin Schweiger]

Ok, first we should try to figure out the "node without pixel support" problem. This would be caused by the specified bounding box being smaller than the mesh extents. This might be caused by rounding errors. You could just try to create the basis mapper without explicitly specifying a bounding box. Toast then automatically calculates the bounding box from the mesh extents, I.e.

basis = toastBasis(mesh,grd);

The second problem of ending up with inf values in the system matrix may be a result of the first problem.

Can you let me know if this solves the problem? Otherwise I'll look into this in more detail.