Method AddTransSIMDIR

comp/l2hofespace.cpp:3580–3639

    using DiffOp<DiffOpGradVectorL2Piola>::AddTransSIMDIR;
    static void AddTransSIMDIR (const FiniteElement & bfel, const SIMD_BaseMappedIntegrationRule & bmir,
                                BareSliceMatrix<SIMD<double>> y, BareSliceVector<double> x)
    {
      auto & mir = static_cast<const SIMD_MappedIntegrationRule<DIM_SPACE,DIM_SPACE>&> (bmir);
      auto & fel = static_cast<const VectorFiniteElement&> (bfel);
      auto & feli = static_cast<const BaseScalarFiniteElement&> (fel[0]);
      size_t ndofi = feli.GetNDof();

      STACK_ARRAY(SIMD<double>, mem, DIM_SPC*DIM_SPC*mir.Size());
      FlatMatrix<SIMD<double>> grad(DIM_SPC*DIM_SPC, mir.Size(), &mem[0]);
      grad = SIMD<double>(0.0);

      for (size_t i = 0; i < mir.Size(); i++)
        {
          auto trans = 1/(mir[i].GetJacobiDet())*mir[i].GetJacobian();
          for (int j = 0; j < DIM_SPC; j++)
            for (int k = 0; k < DIM_SPC; k++)
              for (int l = 0; l < DIM_SPC; l++)
                grad(j*DIM_SPC+k, i) += trans(l,j)*y(l*DIM_SPC+k, i);
        }

      for (size_t k = 0; k < DIM_SPC; k++)
        feli.AddGradTrans (mir, grad.Rows(k*DIM_SPC, (k+1)*DIM_SPC), x.Range(k*ndofi, (k+1)*ndofi));

      if (!mir.GetTransformation().IsCurvedElement())
        return;

      STACK_ARRAY(SIMD<double>, mem2, DIM_SPC*mir.Size());
      FlatMatrix<SIMD<double>> val(DIM_SPC, mir.Size(), &mem2[0]);
      val = SIMD<double>(0.0);

      for (size_t ip = 0; ip < mir.Size(); ip++)
        {
          auto jac = mir[ip].GetJacobian();
          auto inv = mir[ip].GetJacobianInverse();
          auto invJ = 1/mir[ip].GetJacobiDet();
          Vec<DIM_SPC, Mat<DIM_SPC,DIM_SPC,SIMD<double>>> hesse;
          mir[ip].CalcHesse(hesse);

          Vec<DIM_SPC, Mat<DIM_SPC,DIM_SPC,SIMD<double>>> invjac_hesse;
          for (int i = 0; i < DIM_SPC; i++)
            invjac_hesse(i) = Trans(inv) * hesse(i);

          Vec<DIM_SPC,SIMD<double>> inv_hesse = SIMD<double>(0.0);
          for (int i = 0; i < DIM_SPC; i++)
            for (int j = 0; j < DIM_SPC; j++)
              inv_hesse(i) += invjac_hesse(j)(j,i);
          inv_hesse = Trans(inv) * inv_hesse;

          for (int i = 0; i < DIM_SPC; i++)
            for (int j = 0; j < DIM_SPC; j++)
              for (int k = 0; k < DIM_SPC; k++)
                val(k,ip) +=
                  invJ * (invjac_hesse(i)(j,k)-inv_hesse(j)*jac(i,k)) *
                  y(i*DIM_SPC+j,ip);
        }

      for (size_t k = 0; k < DIM_SPC; k++)