#include "types.h"
#include "entity_data.h"

/* Field derivatives kernel, case for elements with planar faces */
void compute_field_derivatives(const entity_data_cpu<1>& ed,
    const vecxd& f, vecxd& df_dx, vecxd& df_dy, vecxd& df_dz)
{
    /* The elements must be ordered by orientation -> entity::sort_by_orientation() */
    size_t orient_elem_base = 0;
    for (size_t iO = 0; iO < ed.num_orientations; iO++)
    {
        size_t num_elems = ed.num_elems[iO];
        size_t orient_base = orient_elem_base * ed.num_bf;

        for (size_t iT = 0; iT < num_elems; iT++)
        {
            size_t ent_iT = orient_elem_base + iT;
            size_t dofofs = ed.dofs_base + orient_base + iT*ed.num_bf;

            for (size_t odof = 0; odof < ed.num_bf; odof++)
            {
                double acc_df_dx = 0.0;
                double acc_df_dy = 0.0;
                double acc_df_dz = 0.0;

                for (size_t refsys_dir = 0; refsys_dir < 3; refsys_dir++)
                {
                    for (size_t idof = 0; idof < ed.num_bf; idof++)
                    {
                        auto dm_row = iO*ed.num_bf + odof;
                        auto dm_col = refsys_dir*ed.num_bf + idof;
                        double D = ed.differentiation_matrix(dm_row, dm_col);
                        acc_df_dx += ed.jacobians(3*ent_iT + 0, refsys_dir) * D * f(dofofs+idof);
                        acc_df_dy += ed.jacobians(3*ent_iT + 1, refsys_dir) * D * f(dofofs+idof);
                        acc_df_dz += ed.jacobians(3*ent_iT + 2, refsys_dir) * D * f(dofofs+idof);
                    }
                }

                df_dx(dofofs+odof) = acc_df_dx;
                df_dy(dofofs+odof) = acc_df_dy;
                df_dz(dofofs+odof) = acc_df_dz;
            }
        }

        orient_elem_base += num_elems;
    }
}