Weak Dirichlet Boundary Conditions for Nedelec spaces

[hughcars]

Hi,

Thank you in advance for any help, I'm new to MFEM and looking forward to learning more about it.

I'm looking to apply a weak Dirichlet boundary condition on a Nedelec space, Nitsche style. In particular, this is a homogeneous boundary condition, so i just need to add in a bilinear term. The reason being, without weak boundary conditions, adjoint consistency is lost and the convergence of output functionals drops from 2p to p+1, which is not ideal for higher order.

To do this, I'm building off of the Symmetric Interior Penalty method, similar in spirit to Buffa. The basic form ends up looking like:

- < n x v, Q curl(u) > - < Q^T curl(v), n x u> + alpha < n x v, n x u>

where alpha scales with p^2 / h. I'm constructing the length scale 1/h with |S| / |K|, the volume of the surface divided by the volume of the attached element. To this end, I have made another Integrator based off of the DGDiffusionIntegrator, which achieves a similar goal for L2/H1 spaces. My issue is, the resulting boundary condition doesn't seem to be being applied. I have attached my new integrator.

I believe there must be an alternative to achieve this, the third term is very close to a VectorFEMassIntegrator with an element dependent length scale, the first seems a MixedVectorGradientIntegrator for the flux term and the second with a Transpose operator too.

My questions then boil down to :

1. Why is my integrator below not working?
2. How could I get a mesh dependent coefficient into a VectorFEMassIntegrator?

void TangentialDirichletBoundaryIntegrator::AssembleFaceMatrix(
   const mfem::FiniteElement &el1, const mfem::FiniteElement &el2,
   mfem::FaceElementTransformations &Trans, mfem::DenseMatrix &elmat)
{
  using namespace mfem;
  MFEM_ASSERT(el1.GetMapType() == mfem::FiniteElement::H_CURL, "Valid elements must be in H(curl)");
  MFEM_ASSERT(Trans.Elem2No < 0, "Can only be used for boundary faces");

  MFEM_ASSERT( alpha > 0, "")

  auto const dimc = el1.GetCurlDim();
  auto const ndof = el1.GetDof();
  auto const dimv = el1.GetVDim(); // the dimension the vector dof exist in.
  
  MFEM_ASSERT(MQ == nullptr || (dimc == dim && MQ->GetHeight() == dimc && MQ->GetWidth() == dimc),
    "If a matrix coefficient is provided, must be 3D and MQ must match.");

  mfem::Vector normal; ///< outward pointing normal, ndim
  mfem::DenseMatrix ncross; ///< matrix equivalent to left cross product with the normal, dimc x dimv 
  mfem::DenseMatrix shape; ///< vector shape function, ndof x dimv 
  mfem::DenseMatrix curlshape; ///< curl of vector shape function, ndof x dimc
  mfem::DenseMatrix qt_ncross_shape; ///< q^T times normal cross with vector shape function, ndof x dimc
  mfem::DenseMatrix ncross_shape; ///< temporary matrix for storing tangential flux, ndof x dimc
  mfem::DenseMatrix mq; ///< evaluation of matrix coefficient, dimv x dimv
  mfem::DenseMatrix jmat; ///< storage for the Dirichlet penalty term (a.k.a jump term), ndof x ndof.

  normal.SetSize(dimv);
  ncross.SetSize(dimc, dimv); // in 3D = 3x3, in 2D = 1x2

  if (MQ != nullptr){
    MFEM_ASSERT(dim > 2, "Using a matrix coefficent, need to be in 3D");
    mq.SetSize(dimc);
  }
    
  shape.SetSize(ndof, dimv);
  curlshape.SetSize(ndof, dimc);
  ncross_shape.SetSize(ndof, dimc);
  ncross_shape.SetSize(ndof, dimc);

  elmat.SetSize(ndof);
  elmat = 0.0;

  jmat.SetSize(ndof);
  jmat = 0.;  

  const mfem::IntegrationRule *ir = IntRule;
  if (ir == NULL)
  {
    // a simple choice for the integration order; is this OK?
    int const order = 2 * el1.GetOrder() + Trans.OrderW();
        
    ir = &IntRules.Get(Trans.GetGeometryType(), order);
  }

  // construct matrix multiplication equivalent to left cross product
  // a x b = cross_mat(a) * b, where a, b are column vectors
  // For a 3 vector, this is a dense 3x3 matrix
  // For a 2 vector, this is a 1x2 matrix. 
  // Reason being that the cross product exists purely in the out-of-plane direction,
  // and the cross is thus stored as a scalar
  auto cross_mat = [dimv](mfem::Vector const& a, mfem::DenseMatrix& m) {
    if (dimv == 3) {
        // skew-symmetric matrix 
        m(0,0) =     0; m(0,1) = -a[2]; m(0,2) =  a[1];
        m(1,0) =  a[2]; m(1,1) =     0; m(1,2) = -a[0];
        m(2,0) = -a[1]; m(2,1) =  a[0]; m(2,2) =     0;
    } else if (dimv == 2) {
        m(0,0) = -a[1]; m(0,1) = a[0]; // The (2,:) row, skipping the diagonal.
    }
  };

  /*
    assemble: < n x v, Q \curl(u) >  --> elmat
              alpha * p^2 / h < n x v, Q (n x u) >  --> jmat
    
    elmat is constructed by building the first integral, and adding the transposition.
    This makes the adjoint consistency readily apparent. 
        
    jmat can constructed by writing it as weighted inner product of the basis functions, 
    where the weight matrix is rank deficient in the normal direction. 
    < n x v, Q (n x u) > = < [n]_x v, Q [n]_x u> 
                        = v^T W u
                        = (v, u)_X
    where [n]_x is the left cross product of n, written as a matrix operator, 
    and W = [n]_x^T Q [n]_x is rank deficient. 

    Alternatively, we rewrite using the common weight 
    < Q^T (n x v), curl(u) > --> elmat 
    < Q^T (n x v), (alpha * p^2 / h) * (n x u) > --> jmat  
    
    Ultimately, elmat will be set using:
    elmat := - elmat - elmat^t + jmat
  */
  for (int p = 0; p < ir->GetNPoints(); p++)
  {
    auto const& ip = ir->IntPoint(p);

    // Set the integration point in the face and the neighboring elements
    Trans.SetAllIntPoints(&ip);

    // g index is for quadrature points. Common weight for evaluating the integral on the reference FACE element.
    double const w = ip.weight * Trans.Weight(); // w_g * det(J)

    // Access the neighboring element's integration points
    auto const& eip1 = Trans.GetElement1IntPoint();

    CalcOrtho(Trans.Jacobian(), normal); // normal is perpendicular to the face
    cross_mat(normal, ncross); // construct the matrix operator equivalent to left cross product

    // Use Phys versions because then we don't need to worry about the reference basis scaling.
    // This is computationally potentially costly, but helps reduce one unknown
    el1.CalcPhysVShape(*Trans.Elem1, shape); // \vec{phi}_i -- ndof x dimv
    el1.CalcPhysCurlShape(*Trans.Elem1, curlshape); // \curl(\vec{\phi}_i) -- ndof x dimc

    // (n x phi) = (ndof x dimc) -> the tangential component of the basis vectors
    // Note: in 3D this is a full vector (with no component in n).
    // In 2D, the result is always "out of plane", and stored as a scalar
    // NB: shape is a stack of row vectors corresponding to each basis vector, 
    // thus evaluating (n x phi)^T = phi^T * [n]_x^T
    MFEM_VERIFY( shape.Width() == ncross.Width()
              && shape.Height() == ncross_shape.Height()
              && ncross.Height() == ncross_shape.Width(), "Matrix Product invalid");
    MultABt(shape, ncross, ncross_shape); // re-use the tangential flux vector for this intermediate evaluation

    // evaluate the scaling coefficient, before applying to the ncross_shape temporary variable
    // in order to fill the weighting function qt_ncross_shape
    if (MQ == nullptr)
    {
        // Given the weight is a scalar, can fast evaluate the weighting
        qt_ncross_shape = ncross_shape;
        qt_ncross_shape *= (Q != nullptr ? Q->Eval(*Trans.Elem1, eip1) : 1.0); 
    }
    else
    {
        // If there is a matrix coefficient, evaluate it.
        MQ->Eval(mq, *Trans.Elem1, eip1);

        MFEM_VERIFY( ncross_shape.Width() == mq.Height()
                  && ncross_shape.Height() == qt_ncross_shape.Height()
                  && mq.Width() == qt_ncross_shape.Width(), "Matrix Product invalid");
        Mult(ncross_shape, mq, qt_ncross_shape);
    }

    // At this point, ncross_shape is loaded with (n x v), thus we first evaluate jmat 
    double const invh = Trans.Weight() / Trans.Elem1->Weight(); // (n^T n) / |J_1|

    double const wq = alpha * w * std::pow(el1.GetOrder(), 2.0) * invh; // alpha * w_g * p^2 / h

    // jmat
    MFEM_VERIFY( qt_ncross_shape.Width() == ncross_shape.Width()
              && qt_ncross_shape.Height() == jmat.Height()
              && ncross_shape.Height() == jmat.Width(), "Matrix Product invalid");
    AddMult_a_ABt(wq, qt_ncross_shape, ncross_shape, jmat); 

    // elmat
    MFEM_VERIFY( qt_ncross_shape.Width() == curlshape.Width()
              && qt_ncross_shape.Height() == elmat.Height()
              && curlshape.Height() == elmat.Width(), "Matrix Product invalid");
    AddMult_a_ABt(w, qt_ncross_shape, curlshape, elmat); // += w * (ndof x dimc) * (dimc x ndof)
  }

  // elmat := -elmat - elmat^t + jmat

  // efficiently
  // for (int i = 0; i < ndof; i++)
  // {
  //   for (int j = 0; j < i; j++)
  //   {
  //       double aij = elmat(i, j), aji = elmat(j, i), mij = jmat(i, j);
  //       elmat(i, j) = - aji - aij + mij;
  //       elmat(j, i) = - aij - aji + mij;
  //   }
  //   elmat(i, i) = - 2 * elmat(i, i) + jmat(i, i);
  // }

  // lazily
  elmat *= -1;

  DenseMatrix elmatt = elmat; 
  elmatt.Transpose();

  elmat += elmatt;
  elmat += jmat;
}