Method DiffJacobi

fem/coefficient.cpp:3734–3787

  shared_ptr<CoefficientFunction> DiffJacobi (const CoefficientFunction * var, typename BASE::T_DJC & cache) const override
  {
    auto thisptr = const_pointer_cast<CoefficientFunction>(this->shared_from_this());
    if (cache.find(thisptr) != cache.end())
      return cache[thisptr];
    
    int h = c1->Dimension();
    int dimvar = var->Dimension();
    
    if (this == var)
      return IdentityCF(h);

    Array<int> dimres{h};
    dimres += var->Dimensions();

    shared_ptr<CoefficientFunction> res;
    if (false) //first version, a little bit slower
      {
        auto diffc1 = c1 -> DiffJacobi (var, cache);
        auto diffc2 = c2 -> DiffJacobi (var, cache);
        
        Array<shared_ptr<CoefficientFunction>> cross1(dimvar);
        Array<shared_ptr<CoefficientFunction>> cross2(dimvar);
        for (auto i : Range(dimvar))
          {
            cross1[i] = CrossProduct(MakeSubTensorCoefficientFunction(diffc1,i,Array<int>({3}),Array<int>({dimvar})),c2);
            cross2[i] = CrossProduct(c1,MakeSubTensorCoefficientFunction(diffc2,i,Array<int>({3}),Array<int>({dimvar})));
          }
        
        auto dcross1 = MakeVectorialCoefficientFunction (std::move(cross1))->Reshape(dimvar,h) -> Transpose() -> Reshape(dimres);
        auto dcross2 = MakeVectorialCoefficientFunction (std::move(cross2))->Reshape(dimvar,h) -> Transpose() -> Reshape(dimres);
    
        res = dcross1 + dcross2;
      }
    else //second version, a little bit faster
      {
        auto c11 = MakeComponentCoefficientFunction(c1,0);
        auto c12 = MakeComponentCoefficientFunction(c1,1);
        auto c13 = MakeComponentCoefficientFunction(c1,2);
        auto c21 = MakeComponentCoefficientFunction(c2,0);
        auto c22 = MakeComponentCoefficientFunction(c2,1);
        auto c23 = MakeComponentCoefficientFunction(c2,2);

        Array<shared_ptr<CoefficientFunction>> cross(3);
        cross[0] = c12*c23 - c13*c22;
        cross[1] = c13*c21 - c11*c23;
        cross[2] = c11*c22 - c12*c21;

        res = MakeVectorialCoefficientFunction(std::move(cross)) -> DiffJacobi (var, cache);
        
      }
    cache[thisptr] = res;
    return res;
  }
  
};