template <typename T>
bool feMatrix<T>::getBoundaryNodeIndices(ot::DA* da, std::vector<unsigned int> &indices, std::vector<double>& coords, std::function<double(double,double,double, int)> f) {
  // loop over elements, writable? and maintain list of indices and coords that are close enough.
  std::map<unsigned int, Point> bdy;

  unsigned int maxD = da->getMaxDepth();
  unsigned int lev;
  double hx, hy, hz, dist, half;
  Point pt;

  double xFac = m_dLx/((double)(1<<(maxD-1)));
  double yFac = m_dLy/((double)(1<<(maxD-1)));
  double zFac = m_dLz/((double)(1<<(maxD-1)));
  double xx[8], yy[8], zz[8];
  unsigned int idx[8];

  for ( da->init<ot::DA_FLAGS::WRITABLE>(); da->curr() < da->end<ot::DA_FLAGS::WRITABLE>(); da->next<ot::DA_FLAGS::WRITABLE>() ) {
    // set the value
    lev = da->getLevel(da->curr());
    hx = xFac * (1 << (maxD - lev));
    hy = yFac * (1 << (maxD - lev));
    hz = zFac * (1 << (maxD - lev));

    half = hx;
    if (half > hy)
    half = hy;
    if (half > hz)
    half = hz;

    half = half * 0.5;

    pt = da->getCurrentOffset();

    //! get the correct coordinates of the nodes ...
    unsigned char hangingMask = da->getHangingNodeIndex(da->curr());

    xx[0] = pt.x() * xFac;
    yy[0] = pt.y() * yFac;
    zz[0] = pt.z() * zFac;
    xx[1] = pt.x() * xFac + hx;
    yy[1] = pt.y() * yFac;
    zz[1] = pt.z() * zFac;
    xx[2] = pt.x() * xFac;
    yy[2] = pt.y() * yFac + hy;
    zz[2] = pt.z() * zFac;
    xx[3] = pt.x() * xFac + hx;
    yy[3] = pt.y() * yFac + hy;
    zz[3] = pt.z() * zFac;

    xx[4] = pt.x() * xFac;
    yy[4] = pt.y() * yFac;
    zz[4] = pt.z() * zFac + hz;
    xx[5] = pt.x() * xFac + hx;
    yy[5] = pt.y() * yFac;
    zz[5] = pt.z() * zFac + hz;
    xx[6] = pt.x() * xFac;
    yy[6] = pt.y() * yFac + hy;
    zz[6] = pt.z() * zFac + hz;
    xx[7] = pt.x() * xFac + hx;
    yy[7] = pt.y() * yFac + hy;
    zz[7] = pt.z() * zFac + hz;

    da->getNodeIndices(idx);
    for (int i = 0; i < 8; ++i) {
      if (!(hangingMask & (1u << i))) {
    for (int j = 1; j <= 6; j++) {
          dist = f(xx[i], yy[i], zz[i], j);
          if (dist < 1e-6) {
            bdy[idx[i]] = Point(xx[i], yy[i], zz[i]);
        //std::cout<<"idx: "<<idx[i]<<std::endl;
      }
        }
      }
    }
  } // loop over elements

  // extract the indices and coords from map
  for(auto const& i: bdy) {
    indices.push_back(i.first);
    coords.push_back(i.second.x());
    coords.push_back(i.second.y());
    coords.push_back(i.second.z());
  }
} // getBdy