Implement own bilinear form

[AIBSCT80]

Implement own bilinear form

[atallah727]

Hello @YHDAUT, you can check the shifted miniapp as an example. There we define our own integrators over the boundary. But the main idea is to create a class that inherits from BilinearFormIntegrator (or LinearFormIntegrator ) and then you would overload the assemble function. For integrals on the domain, as in your case, you would need to overload the AssembleElementMatrix function.

[AIBSCT80]

Hello @atallah727 . Many thanks. May I ask where I can see this shifted miniapp example. Are there any example defining own integrators over the domain? Your comments help me a lot. Thanks..

[atallah727]

You can find the shifted miniapp in the miniapps directory. For more details on writing domain integrators you can refer to this documentation page: https://mfem.org/integration/. Look for the section titled: Writing Custom Integrators .

[AIBSCT80]

Hi @atallah727 . Many thanks for your comments. This is really a great help to me. :)

[atallah727]

Great. I will go ahead and close this issue. But please feel free to reopen it.

[AIBSCT80]

Hi @atallah727 it seems to me not trivial to implement the mixed bilinear form :

\partial_x u \partial_x v + \partial_y u \partial_y v - \partial_z u \partial_z v

I am trying to make this new bilinear form from the standard
diffusion integrator in MFEM.

However, is there any way/example to show how to
select derivative in each direction x, y, z

\partial_x u \partial_x v
\partial_y u \partial_y v
\partial_z u \partial_z v

Then perform summation or subtraction of these terms.

Thanks

void MyDiffusionIntegrator::AssembleElementMatrix2(
const FiniteElement &trial_fe, const FiniteElement &test_fe,
ElementTransformation &Trans, DenseMatrix &elmat)
{
int tr_nd = trial_fe.GetDof();
int te_nd = test_fe.GetDof();
dim = trial_fe.GetDim();
int spaceDim = Trans.GetSpaceDim();
bool square = (dim == spaceDim);
double w;

if (VQ)
{
MFEM_VERIFY(VQ->GetVDim() == spaceDim,
"Unexpected dimension for VectorCoefficient");
}
if (MQ)
{
MFEM_VERIFY(MQ->GetWidth() == spaceDim,
"Unexpected width for MatrixCoefficient");
MFEM_VERIFY(MQ->GetHeight() == spaceDim,
"Unexpected height for MatrixCoefficient");
}

#ifdef MFEM_THREAD_SAFE
DenseMatrix dshape(tr_nd, dim), dshapedxt(tr_nd, spaceDim);
DenseMatrix te_dshape(te_nd, dim), te_dshapedxt(te_nd, spaceDim);
DenseMatrix invdfdx(dim, spaceDim);
DenseMatrix dshapedxt_m(te_nd, MQ ? spaceDim : 0);
DenseMatrix M(MQ ? spaceDim : 0);
Vector D(VQ ? VQ->GetVDim() : 0);
#else
dshape.SetSize(tr_nd, dim);
dshapedxt.SetSize(tr_nd, spaceDim);
te_dshape.SetSize(te_nd, dim);
te_dshapedxt.SetSize(te_nd, spaceDim);
invdfdx.SetSize(dim, spaceDim);
dshapedxt_m.SetSize(te_nd, MQ ? spaceDim : 0);
M.SetSize(MQ ? spaceDim : 0);
D.SetSize(VQ ? VQ->GetVDim() : 0);
#endif
elmat.SetSize(te_nd, tr_nd);

const IntegrationRule *ir = IntRule ? IntRule : &GetRule(trial_fe, test_fe);

elmat = 0.0;
for (int i = 0; i < ir->GetNPoints(); i++)
{
const IntegrationPoint &ip = ir->IntPoint(i);
trial_fe.CalcDShape(ip, dshape);
test_fe.CalcDShape(ip, te_dshape);

  Trans.SetIntPoint(&ip);
  CalcAdjugate(Trans.Jacobian(), invdfdx);
  w = Trans.Weight();
  w = ip.weight / (square ? w : w*w*w);
  Mult(dshape, invdfdx, dshapedxt);
  Mult(te_dshape, invdfdx, te_dshapedxt);
  // invdfdx, dshape, and te_dshape no longer needed
  if (MQ)
  {
     MQ->Eval(M, Trans, ip);
     M *= w;
     Mult(te_dshapedxt, M, dshapedxt_m);
     AddMultABt(dshapedxt_m, dshapedxt, elmat);
  }
  else if (VQ)
  {
     VQ->Eval(D, Trans, ip);
     D *= w;
     AddMultADAt(dshapedxt, D, elmat);
  }
  else
  {
     if (Q)
     {
        w *= Q->Eval(Trans, ip);
     }
     dshapedxt *= w;
     AddMultABt(te_dshapedxt, dshapedxt, elmat);
  }

}
}