Function GetIBernsteinBasis

comp/webgui.cpp:47–79

}

Matrix<double> GetIBernsteinBasis(ngfem::ELEMENT_TYPE etype, int order) {
  std::map<std::pair<ngfem::ELEMENT_TYPE, int>, Matrix<double>> cache;
  if(cache.count({etype, order})) return cache[{etype, order}];
  if (etype == ET_SEGM) {
    Matrix<double> ret(order + 1, order + 1);
    ret = 0.;
    for (auto i : Range(order + 1))
      for (auto j : Range(order + 1))
        ret(i, j) = Bernstein(1. * i / order, j, order);
    CalcInverse(ret);
    cache.insert({{etype, order}, ret});
    return ret;
  }

  if (etype == ET_TRIG) {
    int n = (order + 1) * (order + 2) / 2;
    Matrix<double> ret(n, n);
    int ii = 0;
    for (auto ix : Range(order + 1))
      for (auto iy : Range(order + 1 - ix)) {
        int jj = 0;
        for (auto jx : Range(order + 1))
          for (auto jy : Range(order + 1 - jx))
            ret(ii, jj++) =
                BernsteinTrig(1. * ix / order, 1. * iy / order, jx, jy, order);
        ii++;
      }
    CalcInverse(ret);
    cache.insert({{etype, order}, ret});
    return ret;
  }
  throw Exception("Element type not supported");
}

const IntegrationRule & GetElementPoints(ngfem::ELEMENT_TYPE etype, int order) {
  static std::map<std::pair<ngfem::ELEMENT_TYPE, int>, IntegrationRule> cache;