Implementing smoothness indicators with ProjectGrad

[FH-9]

Implementing smoothness indicators with ProjectGrad

[FH-9]

I realized the integral computation was way off. I think after projecting the gradient degrees of freedom back onto the finite dimensional space of the solution, the integral, $\int \left( \frac{d}{dx} \tilde{p} \right)^2 dx \approx \int \left(\sum_k \hat{dp}_k \phi_k \right)^2 dx$, can be cast into the the quadratic form $\hat{x}^T M \hat{x}$, where $\hat{x}$ is the vector of gradient degrees of freedom, and $M$ is the extrapolated mass integrator (since the integral is over the target element only). I have updated the code but it is still blowing up. Is my reasoning here correct? And do you notice anything wrong with the code?
Here is the updated code:

vfes->GetElementVDofs(indx, vdof_indices);
sol.GetSubVector(vdof_indices, el_vdofs);
fe->ProjectGrad(*fe, *Tr, grad_matrix);
DenseMatrix vdof_mat(el_vdofs.GetData(), ndofs, num_equations);
// beta is num_equations x num_stencil array
// columns of which are the beta values for the given element in the stencil
beta.GetColumnReference(j, beta_col); 
 .
 .
 .
Vector grad_vec(ndofs), Mdx(ndofs);
DenseMatrix mass_matrix(ndofs);
if (dim == 1)
{
  grad_vdof_mat = vdof_mat;
  for (int grad_level = 1; grad_level < order + 1; grad_level++)
  {
    int_order = 2 * fe->GetOrder() + Tr->OrderW();
    const IntegrationRule* ir = &IntRules.Get(fe->GetGeomType(), int_order);
    mfem::Mult(grad_matrix, grad_vdof_mat, tmp_mat);
    grad_vdof_mat = tmp_mat;
    integrand = 0.0;
    mass_matrix = 0.0;
// Calculating the (extrapolated) mass integrator
    for (int k = 0; k < ir->GetNPoints(); k++)
    {
      const IntegrationPoint &ip = ir->IntPoint(k);
      // j = 0 is the index of the target/troubled element; this is where the integral is evaluated
      if (j != 0)
      {
        Tr_target->Transform(ip, phys_ip);
        Tr->TransformBack(phys_ip, extrap_ip);
        Tr->SetIntPoint(&extrap_ip);
        fe->CalcPhysShape(*Tr, shape);
      }
      else
      {
        Tr->SetIntPoint(&ip);
        fe->CalcPhysShape(*Tr, shape);
      }
      w = Tr->Weight() * ip.weight;
      AddMult_a_VVt(w, shape, mass_matrix);
    }
// x^T M x quadratic form calculation where x are the gradient degrees of freedom
    for (int eq = 0 ; eq < num_equations; eq++)
    {
      grad_vdof_mat.GetColumn(eq, grad_vect);
      mass_matrix.Mult(grad_vect, Mdx);
      integrand(eq) = grad_vect * Mdx;
    }
    integrand /= pow(factorial(grad_level), 2); 
    integrand *= pow(target_el_size, 2 * grad_level - 1);
    beta_col += integrand;
  }
}

Thank you very much.

[FH-9]

I would appreciate if anyone could offer some insight. Thank you.