Consider
Th....triangulation
Xh....finite element space on Xh
lab....label for the boundary part on Th
flab....any function that runs on the boundary, one-to-one and takes non cero values.
zoi... a integer vector that will save the indices of the nods on the boundary with label lab
ndofonlab .... a integer that will save the number of dofs on the boundary with label lab
then, in this macro ....
macro boundarydofs(lab,flab,Th,Xh,zoi,ndofonlab)
{int[int] ib(0:Xh.ndof-1);
varf vbord(u1,v1) = on(lab,u1=flab);
real[int] xb=vbord(0,Xh,tgv=1);
sort(xb,ib);
zoi=ib;
xb = xb ? 1 : 0;
ndofonlab = xb.sum;
}//
in order of appearance we have ....
ib.....is just a vector of size Xh.ndof with values 0,1,2,...,Xh.ndof-1
xb....is a finite element function 0 everywhere except on the piece of boundary lab where it takes the values given by the function flab. Here it is important that flab is nonzero on boundary lab. Then, all ceros in xb are out of the boundary lab.
sort(xb,ib);....sorts the vector xb from max value to min value and the vector ib is also sorted accordingly. Here it is important that flab is one-to-one on boundary lab. Then, all values in xb on boundary lab are different.
zoi....now contains the set of indices 0,1,2,...,Xh.ndof-1 reordered such that it points from the maximun to the minimum of the values of xb
xb = xb ? 1 : 0;....just tells how many non zero elements in xb are. These are suppose to be the values on the boundary lab
ndofonlab.....integer that keeps the number of nods on the boundary lab
Finally, only the first ndofonlab values of zoi will be pointing to the dofs on the boundary lab
For instance, a call for this macro when the boundary is the vertical segment \(\{x=0\} \times \{0\le y\le 1\}\) labelled by g1 is ...
int mm1=u1[].n; int[int] zoi1(mm1); int ndofon1=0; boundarydofs(g1,y-10.,th1,Vh1,zoi1,ndofon1); zoi1.resize(ndofon1);