Remove plist (nodal parameter lists) from Mesh class, closes toastpp#7 (

toastpp#8)

* Removed Mesh::plist and all references to mesh-based parameter sets

* Fixed remaining plist references that were missed when configured with --enable-threads --enable-fdot

* fixed plist reference in libfdot/util.cc; fixed makefile which skipped most targets in libfdot

* updated .gitignore


Former-commit-id: 4d93966

.gitignore

@@ -1,4 +1,6 @@
 
+*~
+
 *.dylib
 
 *.mexmaci64
@@ -25,12 +27,16 @@ lib
 
 linux
 
+linux64
+
 *.xcbkptlist
 
 Makefile.incl
 
 numerics/cuda/cusp
 
+script/matlab/demos/*.mat
+
 toast.env
 
 toastenv.csh

script/matlab/demos/toast_demo1.m

@@ -89,7 +89,7 @@ function init(handles)
 prm.basis.hMesh.ReadQM('circle25_1x32.qm');
 prm.basis.hBasis = toastBasis(prm.basis.hMesh,[prm.bx prm.by],'Linear');
 prm.qvec = prm.basis.hMesh.Qvec('Neumann','Gaussian',2);
-prm.mvec = prm.basis.hMesh.Mvec('Gaussian',2);
+prm.mvec = prm.basis.hMesh.Mvec('Gaussian',2,1.4);
 n = prm.basis.hMesh.NodeCount();
 prm.mua = ones(n,1)*0.025;
 prm.mus = ones(n,1)*2;
@@ -189,7 +189,7 @@ function setsource(handles,phi)
 Q = [rad*cos(phi) rad*sin(phi)];
 prm.basis.hMesh.SetQM(Q,M);
 prm.qvec = prm.basis.hMesh.Qvec('Neumann','Gaussian',2);
-prm.mvec = prm.basis.hMesh.Mvec('Gaussian',2);
+prm.mvec = prm.basis.hMesh.Mvec('Gaussian',2,1.4);
 
 axes(handles.axes7);
 n = prm.basis.hMesh.NodeCount();

script/matlab/demos/toast_demo1a.m

@@ -84,7 +84,7 @@ function init(handles)
 prm.basis.hMesh.ReadQM('circle25_1x32.qm');
 prm.basis.hBasis = toastBasis(prm.basis.hMesh,prm.basis.bdim,'Linear');
 prm.qvec = real(prm.basis.hMesh.Qvec('Isotropic','Gaussian',2));
-prm.mvec = real(prm.basis.hMesh.Mvec('Gaussian',2));
+prm.mvec = real(prm.basis.hMesh.Mvec('Gaussian',2,1.4));
 prm.mvec = prm.mvec(:,[4,8,12,16]);   % only run with 4 detectors
 %prm.smask = toastSolutionMask(prm.basis.hBasis);
 n = prm.basis.hMesh.NodeCount();

script/matlab/demos/toast_demo2.m

@@ -93,7 +93,7 @@ function init(handles)
 prm.basis.hMesh.SetQM(Q,M);
 
 prm.qvec = prm.basis.hMesh.Qvec('Neumann','Gaussian',2);
-prm.mvec = prm.basis.hMesh.Mvec('Gaussian',2);
+prm.mvec = prm.basis.hMesh.Mvec('Gaussian',2,1.4);
 prm.basis.hBasis = toastBasis(prm.basis.hMesh,[prm.bx prm.by],'Linear');
 n = prm.basis.hMesh.NodeCount();
 prm.mua = ones(n,1)*0.025;
@@ -176,7 +176,7 @@ function set_qm(handles,phiq,phim)
 M = [radm*cos(phim) radm*sin(phim)];
 prm.basis.hMesh.SetQM(Q,M);
 prm.qvec = prm.basis.hMesh.Qvec('Neumann','Gaussian',2);
-prm.mvec = prm.basis.hMesh.Mvec('Gaussian',2);
+prm.mvec = prm.basis.hMesh.Mvec('Gaussian',2,1.4);
 prm.phiq = phiq;
 prm.phim = phim;
 

script/matlab/toast2/toastRecon.m

@@ -88,7 +88,7 @@ function toastRecon(prm)
 
 % ----------------------------------------------------------------------
 % Generate measurement vectors
-LPRM.mvec = LPRM.hMesh.Mvec (prm.meas.det.prof, prm.meas.det.width);
+LPRM.mvec = LPRM.hMesh.Mvec (prm.meas.det.prof, prm.meas.det.width, ref);
 fprintf ('Detector vector setup: %d detectors\n', size(LPRM.mvec,2));
 
 % ----------------------------------------------------------------------

script/matlab/toast2/toastSolver.m

@@ -275,7 +275,7 @@ function echo_iteration_state (obj, itr, x, err)
             fprintf ('Source vector set up: %d sources\n', size(localprm.qvec,2));
 
             % Generate measurement vectors
-            localprm.mvec = localprm.hMesh.Mvec (prm.meas.det.prof, prm.meas.det.width);
+            localprm.mvec = localprm.hMesh.Mvec (prm.meas.det.prof, prm.meas.det.width, localprm.ref);
             fprintf ('Detector vector setup: %d detectors\n', size(localprm.mvec,2));
 
             % Initial data set f[x0]
@@ -415,4 +415,4 @@ function solve_iter (this, prm, itr, x, err)
         x0 = [];
         c0 = 0.3;
     end
-end
+end

src/bintools/Makefile.src

@@ -8,7 +8,6 @@ BINS =  getnim \
 	lin2quadmesh \
 	makeqm \
 	makeqm3d \
-	mesh2nim \
 	meshinfo \
 	mkmesh_brick \
 	mkmesh_cyl \

src/bintools/lin2cubicmesh.cc

@@ -52,7 +52,6 @@ void Convert2DMesh (Mesh &mesh)
     nlen = mesh.nlen();
 
     NodeList nnlist (7*elen); // max number of new nodes
-    ParameterList nplist (7*elen);
     int nnlen = 0;
 
     for (el = 0; el < elen; el++) {
@@ -93,15 +92,6 @@ void Convert2DMesh (Mesh &mesh)
 		        nnlist[nnlen].SetBndTp (mesh.nlist[n0].BndTp());
 		    else
 		        nnlist[nnlen].SetBndTp (BND_NONE);
-
-		    nplist[nnlen].SetMua ((2.0-s)/3.0*mesh.plist[n0].Mua() +
-					  (1.0+s)/3.0*mesh.plist[n1].Mua());
-		    nplist[nnlen].SetKappa ((2.0-s)/3.0*mesh.plist[n0].Kappa()+
-					   (1.0+s)/3.0*mesh.plist[n1].Kappa());
-		    nplist[nnlen].SetN ((2.0-s)/3.0*mesh.plist[n0].N() +
-					(1.0+s)/3.0*mesh.plist[n1].N());
-		    nplist[nnlen].SetA ((2.0-s)/3.0*mesh.plist[n0].A() +
-					(1.0+s)/3.0*mesh.plist[n1].A());
 		    nnlen++;
 		} else {
 		  nd[inew] += nlen;
@@ -117,18 +107,6 @@ void Convert2DMesh (Mesh &mesh)
 					  mesh.nlist[pel->Node[1]][j] +
 					  mesh.nlist[pel->Node[2]][j]);
 	}
-	nplist[nnlen].SetMua (1.0/3.0 * (mesh.plist[pel->Node[0]].Mua() +
-					 mesh.plist[pel->Node[1]].Mua() +
-					 mesh.plist[pel->Node[2]].Mua()));
-	nplist[nnlen].SetKappa (1.0/3.0 * (mesh.plist[pel->Node[0]].Kappa() +
-					 mesh.plist[pel->Node[1]].Kappa() +
-					 mesh.plist[pel->Node[2]].Kappa()));
-	nplist[nnlen].SetN (1.0/3.0 * (mesh.plist[pel->Node[0]].N() +
-					 mesh.plist[pel->Node[1]].N() +
-					 mesh.plist[pel->Node[2]].N()));
-	nplist[nnlen].SetA (1.0/3.0 * (mesh.plist[pel->Node[0]].A() +
-					 mesh.plist[pel->Node[1]].A() +
-					 mesh.plist[pel->Node[2]].A()));
 	nd[inew++] = nlen + nnlen;
 	nnlen++;
 
@@ -147,10 +125,8 @@ void Convert2DMesh (Mesh &mesh)
     }
     if (nnlen) {
         mesh.nlist.Append (nnlen);
-	mesh.plist.Append (nnlen);
 	for (i = 0; i < nnlen; i++) {
 	    mesh.nlist[nlen+i].Copy (nnlist[i]);
-	    mesh.plist[nlen+i] = nplist[i];
 	}
     }
 }
@@ -167,7 +143,6 @@ void Convert3DMesh (Mesh &mesh)
     nlen = mesh.nlen();
 
     NodeList nnlist (6*elen); // max number of new nodes
-    ParameterList nplist (6*elen);
     int nnlen = 0;
 
     for (el = 0; el < elen; el++) {
@@ -191,15 +166,6 @@ void Convert3DMesh (Mesh &mesh)
 		    nnlist[nnlen].SetBndTp (mesh.nlist[n0].BndTp());
 		else
 		    nnlist[nnlen].SetBndTp (BND_NONE);
-
-		nplist[nnlen].SetMua ((mesh.plist[n0].Mua() +
-				       mesh.plist[n1].Mua())*0.5);
-		nplist[nnlen].SetKappa ((mesh.plist[n0].Kappa() +
-					 mesh.plist[n1].Kappa())*0.5);
-		nplist[nnlen].SetN ((mesh.plist[n0].N() +
-				     mesh.plist[n1].N())*0.5);
-		nplist[nnlen].SetA ((mesh.plist[n0].A() +
-				     mesh.plist[n1].A())*0.5);
 		nnlen++;
 	    } else {
 	        nd[edge] += nlen;
@@ -214,10 +180,8 @@ void Convert3DMesh (Mesh &mesh)
     }
     if (nnlen) {
         mesh.nlist.Append (nnlen);
-	mesh.plist.Append (nnlen);
 	for (i = 0; i < nnlen; i++) {
 	    mesh.nlist[nlen+i].Copy (nnlist[i]);
-	    mesh.plist[nlen+i] = nplist[i];
 	}
     }
 #endif

src/bintools/lin2quadmesh.cc

@@ -61,7 +61,6 @@ void Lin2Quad (Mesh &mesh)
 
     // allocate lists for new nodes
     NodeList nnlist (maxnew);
-    ParameterList nplist (maxnew);
 
     // find neighbour nodes for each node in the mesh
     cerr << "Setting up node neighbour list\n";
@@ -130,23 +129,16 @@ void Lin2Quad (Mesh &mesh)
         nn = nd_nnbhr[i];
 	nb = nd_nbhr[i];
 	Node &ni = mesh.nlist[i];
-	Parameter &pi = mesh.plist[i];
 	for (j = 0; j < nn; j++) {
 	    if ((nbj = nb[j]) <= i) continue; // to avoid counting twice
 	    Node &nj = mesh.nlist[nbj];
-	    Parameter &pj = mesh.plist[nbj];
 	    // create a new node between i and its jth neighbour
 	    nnlist[nnlen].New(dim);
 	    for (k = 0; k < dim; k++)
 	        nnlist[nnlen][k] = 0.5*(ni[k]+nj[k]);
 	    if (ni.BndTp() == nj.BndTp()) nnlist[nnlen].SetBndTp (ni.BndTp());
 	    else                          nnlist[nnlen].SetBndTp (BND_NONE);
 
-	    nplist[nnlen].SetMua ((pi.Mua()+pj.Mua())*0.5);
-	    nplist[nnlen].SetKappa ((pi.Kappa()+pj.Kappa())*0.5);
-	    nplist[nnlen].SetN ((pi.N()+pj.N())*0.5);
-	    nplist[nnlen].SetA ((pi.A()+pj.A())*0.5);
-
 	    // now find all elements sharing the new node
 	    for (k = 0; k < nelref[i]; k++) {
 	        el = elref[i][k];
@@ -207,10 +199,8 @@ void Lin2Quad (Mesh &mesh)
     cerr << "Finalising mesh\n";
     if (nnlen) {
         mesh.nlist.Append (nnlen);
-	mesh.plist.Append (nnlen);
 	for (i = 0; i < nnlen; i++) {
 	    mesh.nlist[nlen+i].Copy (nnlist[i]);
-	    mesh.plist[nlen+i] = nplist[i];
 	}
     }
 

src/bintools/mkmesh_brick.cc

@@ -121,21 +121,14 @@ int main (void)
     }
     cout << "Generated element list (" << nels << " elements)" << endl;
 
-    // create parameter list
-    mesh.plist.New (nnode);
-    mesh.plist.SetMua (hmua);
-    mesh.plist.SetMus (hmus);
-    mesh.plist.SetN (href);
     for (i = 0; i < mesh.nlist.Len(); i++)
         mesh.nlist[i].SetRegion (hreg);
-    cout << "Generated parameter list" << endl;
 
     for (;;) {
         int nd, cmd;
 	int xmin, xmax, ymin, ymax, zmin, zmax;
-	double mua, mus, refind;
 
-        cout << "Add parameter perturbation (1/0) ? ";
+        cout << "Set region index on node subset (1/0) ? ";
         cin  >> cmd;
 	if (cmd == 0) break;
 	cout << "x-extension of perturbation:" << endl;
@@ -147,8 +140,8 @@ int main (void)
 	cout << "z-extension of perturbation:" << endl;
 	cout << "zmin[0-" << dimz << "] zmax[0-" << dimz << "]: ";
 	cin  >> zmin >> zmax;
-	cout << "mua mus refind region: ";
-	cin  >> mua >> mus >> refind >> reg;
+	cout << "region index: ";
+	cin  >> reg;
 
 	for (k = zmin; k <= zmax; k++) {
 	    if (k < 0 || k > dimz) continue;
@@ -157,9 +150,6 @@ int main (void)
 		for (i = xmin; i <= xmax; i++) {
 		    if (i < 0 || i > dimx) continue;
 		    nd = k * (nnode_slice) + j * (dimy+1) + i;
-		    mesh.plist[nd].SetMua (mua);
-		    mesh.plist[nd].SetMus (mus);
-		    mesh.plist[nd].SetN (refind);
 		    mesh.nlist[nd].SetRegion (reg);
 		}
 	    }