Function MakeTensorTraceCoefficientFunction

fem/coefficient.cpp:4927–4951

shared_ptr<CoefficientFunction>
MakeTensorTraceCoefficientFunction (shared_ptr<CoefficientFunction> c1, int i1, int i2)
{
  if (i1 < 0 || i1 >= c1->Dimensions().Size())
    throw Exception ("TensorTrace, i1 out of range");
  if (i2 < 0 || i2 >= c1->Dimensions().Size())
      throw Exception ("TensorTrace, i2 out of range");
  
  // a simple workaround for the beginning:
  if ( (i1 == 0 && i2 == 1) || (i1 == 1 && i2 == 0) )
    {
      Array<int> dims1 { c1->Dimensions() };
      auto dimsres = dims1.Range(2, END);

      // auto shapematT = TransposeCF (ReshapeCF (std::move(c1), sqr(dims1[0]), -1) );
      // auto prod = shapematT * ReshapeCF ( IdentityCF(dims1[0]), sqr (dims1[0]) );

      auto shapematT = c1 -> Reshape (sqr(dims1[0]), -1) -> Transpose(); 
      auto prod = shapematT * IdentityCF(dims1[0]) -> Reshape (sqr (dims1[0]));

      return prod -> Reshape (Array<int> (dimsres));
    }

  throw Exception("MakeTensorTraceCF not implemented for general case");
}