message: revert the element restriction for nodal finite elements to the original method - this is needed for coarsening mesh nodes in the H1 space - still need to figure out how to handle L2 spaces with nodal elements
diff: @@ -473,7 +473,6 @@ void ScalarFiniteElement::ScalarLocalRestriction(
  const IntegrationRule &ir = IntRules.Get(coarse_fe.GetGeomType(), ir_order);
  // integrate coarse_mass in the coarse space
  -
  for (int i = 0; i < ir.GetNPoints(); i++)
  {
  const IntegrationPoint &c_ip = ir.IntPoint(i);
  @@ -482,7 +481,6 @@ void ScalarFiniteElement::ScalarLocalRestriction(
  }
  // integrate coarse_fine_mass in the fine space
  -
  Trans.SetIntPoint(&Geometries.GetCenter(geom_type));
  for (int i = 0; i < ir.GetNPoints(); i++)
  {
  @@ -638,54 +636,30 @@ void InvertLinearTrans(ElementTransformation &trans,
  void NodalFiniteElement::GetLocalRestriction(ElementTransformation &Trans,
  DenseMatrix &R) const
  {
  - // todo: figure out what to do with continuous spaces, and remove
  - // duplication with ScalarFiniteElement::GetLocalRestriction.
  -
  - // General "restriction", defined by L2 projection
  - double v[Geometry::MaxDim];
  - Vector vv(v, dim);
  - IntegrationPoint c_ip;
  + IntegrationPoint ipt;
  + Vector pt(&ipt.x, dim);
  - const FiniteElement& coarse_fe{*this};
  - const int cs = coarse_fe.GetDof(), fs = this->GetDof();
  - R.SetSize(cs, fs);
  - Vector fine_shape(fs), coarse_shape(cs);
  - DenseMatrix coarse_mass(cs), coarse_fine_mass(cs, fs); // initialized with 0
  - const int ir_order = GetOrder() + coarse_fe.GetOrder();
  - const IntegrationRule &ir = IntRules.Get(coarse_fe.GetGeomType(), ir_order);
  +#ifdef MFEM_THREAD_SAFE
  + Vector c_shape(dof);
  +#endif
  - // integrate coarse_mass in the coarse space
  + Trans.SetIntPoint(&Nodes[0]);
  - for (int i = 0; i < ir.GetNPoints(); i++)
  + for (int j = 0; j < dof; j++)
  {
  - const IntegrationPoint &c_ip = ir.IntPoint(i);
  - coarse_fe.CalcShape(c_ip, coarse_shape);
  - AddMult_a_VVt(c_ip.weight, coarse_shape, coarse_mass);
  - }
  -
  - // integrate coarse_fine_mass in the fine space
  -
  - Trans.SetIntPoint(&Geometries.GetCenter(geom_type));
  - for (int i = 0; i < ir.GetNPoints(); i++)
  + InvertLinearTrans(Trans, Nodes[j], pt);
  + if (Geometries.CheckPoint(geom_type, ipt)) // do we need an epsilon here?
  {
  - const IntegrationPoint &f_ip = ir.IntPoint(i);
  - this->CalcShape(f_ip, fine_shape);
  - Trans.Transform(f_ip, vv);
  - c_ip.Set(v, dim);
  - coarse_fe.CalcShape(c_ip, coarse_shape);
  - AddMult_a_VWt(f_ip.weight*Trans.Weight(), coarse_shape, fine_shape, coarse_fine_mass);
  + CalcShape(ipt, c_shape);
  + R.SetRow(j, c_shape);
  }
  -
  - DenseMatrixInverse coarse_mass_inv(coarse_mass);
  - coarse_mass_inv.Mult(coarse_fine_mass, R);
  -
  - if (map_type == INTEGRAL)
  + else
  {
  - // todo: is this necessary?
  - // assuming Trans is linear; this should be ok for all refinement types
  - Trans.SetIntPoint(&Geometries.GetCenter(geom_type));
  - R *= 1.0 / Trans.Weight();
  + // Set the whole row to avoid valgrind warnings in R.Threshold().
  + R.SetRow(j, infinity());
  + }
  }
  + R.Threshold(1e-12);
  }
  void NodalFiniteElement::Project (
  