21 #ifndef PATCH_POINT_VALUES_HH_
22 #define PATCH_POINT_VALUES_HH_
24 #include <Eigen/Dense>
35 template<
unsigned int spacedim>
class ElOp;
121 template<
unsigned int spacedim = 3>
191 elOp.allocate_result(sizes[elOp.size_type()], *
patch_arena_);
216 for (
uint i_col=0; i_col<coords.n_cols; ++i_col)
217 for (
uint i_row=0; i_row<coords.n_rows; ++i_row) {
219 coords_mat(i_row, i_col)(i_elem) = coords(i_row, i_col);
239 for (
uint i_col=0; i_col<elm_coords.n_cols; ++i_col)
240 for (
uint i_row=0; i_row<elm_coords.n_rows; ++i_row) {
242 coords_mat(i_row, i_col)(i_elem) = elm_coords(i_row, i_col);
252 for (
uint i_col=0; i_col<side_coords.n_cols; ++i_col)
253 for (
uint i_row=0; i_row<side_coords.n_rows; ++i_row) {
255 coords_mat(i_row, i_col)(i_elem) = side_coords(i_row, i_col);
272 uint point_pos = i_point_on_elem *
n_elems_ + elem_table_row;
273 int_vals_(0)(point_pos) = value_patch_idx;
274 int_vals_(1)(point_pos) = elem_table_row;
392 for (
uint i=0; i<3; ++i)
398 const auto &op_matrix =
operations_[op_idx].result_matrix();
400 for (
uint i=0; i<3; ++i)
401 val(i) = op_matrix(i)(op_matrix_idx);
413 for (
uint i=0; i<3; ++i)
414 for (
uint j=0; j<3; ++j)
420 const auto &op_matrix =
operations_[op_idx].result_matrix();
422 for (
uint i=0; i<3; ++i)
423 for (
uint j=0; j<3; ++j)
424 val(i,j) = op_matrix(i,j)(op_matrix_idx);
439 const auto &
mat = op.result_matrix();
440 for (
uint i_mat=0; i_mat<
mat.rows()*
mat.cols(); ++i_mat) {
441 if (
mat(i_mat).data_size()==0) stream <<
"<empty>";
443 const double *vals =
mat(i_mat).data_ptr();
444 for (
size_t i_val=0; i_val<
mat(i_mat).data_size(); ++i_val)
445 stream << vals[i_val] <<
" ";
455 for (
uint i_col=0; i_col<3; ++i_col)
456 stream <<
int_vals_(i_col)(i_row) <<
" ";
472 {
"el_coords",
"jacobian",
"inv_jac",
"jac_det",
"pt_coords",
"weights",
"JxW",
"",
"",
"",
"",
"" },
473 {
"el_coords",
"el_jac",
"el_inv_jac",
"side_coords",
"side_jac",
"side_jac_det",
"exp_el_coords",
"exp_el_jac",
"exp_el_inv_jac",
474 "exp_side_coords",
"exp_side_jac",
"exp_side_jac_det",
"pt_coords",
"weights",
"JxW",
"normal_vec",
"",
"",
"",
"",
"" }
476 stream << std::setfill(
' ') <<
" Operation" << setw(12) <<
"" <<
"Shape" << setw(2) <<
""
477 <<
"Result row" << setw(2) <<
"" <<
"n DOFs" << setw(2) <<
"" <<
"Input operations" << endl;
479 stream <<
" " << std::left << setw(20) << op_names[bulk_side][i] <<
"" <<
" " << setw(6) <<
operations_[i].format_shape() <<
"" <<
" "
480 << setw(11) <<
operations_[i].result_row() <<
"" <<
" " << setw(7) <<
operations_[i].n_dofs() <<
"" <<
" ";
482 for (
auto i_o : input_ops) stream << op_names[bulk_side][i_o] <<
" ";
523 friend class ElOp<spacedim>;
524 template<
unsigned int dim>
526 template<
unsigned int dim>
528 template<
unsigned int dim>
535 template<
unsigned int spacedim = 3>
550 if (shape_vec.size() == 1) shape_vec.push_back(1);
611 result_ = Eigen::Matrix<ArenaVec<double>, Eigen::Dynamic, Eigen::Dynamic>(
shape_[0],
shape_[1]);
617 Eigen::Matrix<ArenaVec<double>, Eigen::Dynamic, Eigen::Dynamic> &
result_matrix() {
625 const Eigen::Matrix<ArenaVec<double>, Eigen::Dynamic, Eigen::Dynamic> &
result_matrix()
const {
633 template<
unsigned int dim1,
unsigned int dim2>
634 Eigen::Map<Eigen::Matrix<ArrayDbl, dim1, dim2>>
value(
TableDbl &op_results,
uint i_dof = 0)
const {
635 return Eigen::Map<Eigen::Matrix<ArrayDbl, dim1, dim2>>(op_results.data() +
result_row_ + i_dof *
n_comp(), dim1, dim2);
640 return Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>>(op_results(
result_row_).data(), dim1, dim2);
645 return Eigen::Map<Eigen::Vector<double, Eigen::Dynamic>>(op_results(
result_row_).data(), op_results(
result_row_).rows());
653 Eigen::Matrix<ArenaVec<double>, Eigen::Dynamic, Eigen::Dynamic>
result_;
674 uint size = op_results(row_begin).rows();
675 for (
int i_pt=size-1; i_pt>=0; --i_pt) {
676 uint el_table_idx = el_table(1)(i_pt);
677 for (
uint i_q=row_begin; i_q<row_end; ++i_q)
678 op_results(i_q)(i_pt) = op_results(i_q)(el_table_idx);
686 template<
unsigned int dim>
690 auto &jac_value = op.result_matrix();
691 const auto &coords_value = operations[ op.input_ops()[0] ].result_matrix();
692 for (
unsigned int i=0; i<3; i++)
693 for (
unsigned int j=0; j<dim; j++)
694 jac_value(i,j) = coords_value(i,j+1) - coords_value(i,0);
696 template<
unsigned int dim>
700 auto &inv_jac_value = op.result_matrix();
701 const auto &jac_value = operations[ op.input_ops()[0] ].result_matrix();
702 inv_jac_value = eigen_arena_tools::inverse<3, dim>(jac_value);
704 template<
unsigned int dim>
708 auto &jac_det_value = op.result_matrix();
709 const auto &jac_value = operations[ op.input_ops()[0] ].result_matrix();
710 jac_det_value(0,0) = eigen_arena_tools::determinant<3, dim>(jac_value).abs();
719 op.allocate_result(point_weights.size(), *arena);
720 auto &weights_value = op.result_matrix();
721 for (
uint i=0; i<point_weights.size(); ++i)
722 weights_value(0,0)(i) = point_weights[i];
726 auto &weights_value = operations[ op.input_ops()[0] ].result_matrix();
727 auto &jac_det_value = operations[ op.input_ops()[1] ].result_matrix();
728 ArenaOVec<double> weights_ovec( weights_value(0,0), weights_value(0,0).data_size() );
729 ArenaOVec<double> jac_det_ovec( jac_det_value(0,0), jac_det_value(0,0).data_size() );
731 auto &jxw_value = op.result_matrix();
732 jxw_value(0,0) = jxw_ovec.
get_vec();
736 auto &op = operations[scalar_shape_op_idx];
737 uint n_dofs = shape_values.size();
738 uint n_points = shape_values[0].size();
739 uint n_elem = op.result_matrix()(0).data_size() / n_points;
741 auto &shape_matrix = op.result_matrix();
744 for (
uint i=0; i<n_elem; ++i) {
749 for (
uint i_dof=0; i_dof<n_dofs; ++i_dof)
750 for (
uint i_p=0; i_p<n_points; ++i_p)
751 ref_vec(i_dof * n_points + i_p) = shape_values[i_dof][i_p];
753 shape_matrix(0) = shape_ovec.
get_vec();
755 template<
unsigned int dim>
758 auto &op = operations[scalar_shape_grads_op_idx];
759 auto &inv_jac_vec = operations[ op.input_ops()[0] ].result_matrix();
760 uint n_points = ref_shape_grads.size();
761 uint n_dofs = ref_shape_grads[0].size();
763 Eigen::Matrix<ArenaVec<double>, dim, 1> ref_grads_vec;
764 Eigen::Matrix<ArenaOVec<double>, dim, 1> ref_grads_ovec;
765 for (
uint i=0; i<ref_grads_vec.rows(); ++i) {
766 ref_grads_vec(i) =
ArenaVec<double>(n_points * n_dofs, op.result_matrix()(0).arena());
767 for (
uint i_dof=0; i_dof<n_dofs; ++i_dof)
768 for (
uint i_p=0; i_p<n_points; ++i_p)
769 ref_grads_vec(i)(i_dof * n_points + i_p) = ref_shape_grads[i_dof][i_p](i);
770 ref_grads_ovec(i) =
ArenaOVec(ref_grads_vec(i), ref_grads_vec(i).data_size());
773 Eigen::Matrix<ArenaOVec<double>, dim, 3> inv_jac_ovec;
774 for (
uint i=0; i<dim*3; ++i) {
775 inv_jac_ovec(i) =
ArenaOVec(inv_jac_vec(i), inv_jac_vec(i).data_size());
778 auto &result_vec = op.result_matrix();
779 Eigen::Matrix<ArenaOVec<double>, 3, 1> result_ovec = inv_jac_ovec.transpose() * ref_grads_ovec;
780 for (
uint i=0; i<3; ++i) {
781 result_vec(i) = result_ovec(i).get_vec();
791 template<
unsigned int dim>
795 auto &jac_value = op.result_matrix();
796 const auto &coords_value = operations[ op.input_ops()[0] ].result_matrix();
797 for (
unsigned int i=0; i<3; i++)
798 for (
unsigned int j=0; j<dim; j++)
799 jac_value(i,j) = coords_value(i,j+1) - coords_value(i,0);
801 template<
unsigned int dim>
804 auto &inv_jac_value = op.result_matrix();
805 const auto &jac_value = operations[ op.input_ops()[0] ].result_matrix();
806 inv_jac_value = eigen_arena_tools::inverse<3, dim>(jac_value);
808 template<
unsigned int dim>
812 auto &jac_value = op.result_matrix();
813 const auto &coords_value = operations[ op.input_ops()[0] ].result_matrix();
814 for (
unsigned int i=0; i<3; i++)
815 for (
unsigned int j=0; j<dim-1; j++)
816 jac_value(i,j) = coords_value(i,j+1) - coords_value(i,0);
818 template<
unsigned int dim>
822 auto &jac_det_value = op.result_matrix();
823 const auto &jac_value = operations[ op.input_ops()[0] ].result_matrix();
861 op.allocate_result(point_weights.size(), *arena);
862 auto &weights_value = op.result_matrix();
863 for (
uint i=0; i<point_weights.size(); ++i)
864 weights_value(0,0)(i) = point_weights[i];
868 auto &weights_value = operations[ op.input_ops()[0] ].result_matrix();
869 auto &jac_det_value = operations[ op.input_ops()[1] ].result_matrix();
870 ArenaOVec<double> weights_ovec( weights_value(0,0), weights_value(0,0).data_size() );
871 ArenaOVec<double> jac_det_ovec( jac_det_value(0,0), jac_det_value(0,0).data_size() );
873 auto &jxw_value = op.result_matrix();
874 jxw_value(0,0) = jxw_ovec.
get_vec();
876 template<
unsigned int dim>
879 auto normal_value = op.value<3, 1>(op_results);
880 auto inv_jac_mat_value = operations[ op.input_ops()[0] ].value<dim, 3>(op_results);
884 norm_vec.resize(normal_value(0).
rows());
885 Eigen::VectorXd A(3);
887 for (
uint i=0; i<normal_value(0).rows(); ++i) {
888 A(0) = normal_value(0)(i);
889 A(1) = normal_value(1)(i);
890 A(2) = normal_value(2)(i);
891 norm_vec(i) = A.norm();
893 for (
uint i=0; i<3; ++i) {
894 normal_value(i) /= norm_vec;
899 auto &op = operations[scalar_shape_op_idx];
900 uint n_points = shape_values[0].size();
902 for (
uint i_row=0; i_row<shape_values[0][0].size(); ++i_row) {
903 ArrayDbl &result_row = op_results( op.result_row()+i_row );
904 for (
uint i_pt=0; i_pt<result_row.rows(); ++i_pt)
905 result_row(i_pt) = shape_values[el_table(3)(i_pt)][i_pt % n_points][i_row];
908 template<
unsigned int dim>
911 auto &op = operations[scalar_shape_grads_op_idx];
912 uint n_points = ref_shape_grads[0].size();
913 uint n_dofs = ref_shape_grads[0][0].size();
915 Eigen::Vector<ArrayDbl, Eigen::Dynamic> ref_shape_grads_expd;
916 ref_shape_grads_expd.resize(dim*n_dofs);
917 for (
uint i=0; i<ref_shape_grads_expd.rows(); ++i)
918 ref_shape_grads_expd(i).resize(op_results(0).
rows());
920 for (
uint i_dof=0; i_dof<n_dofs; ++i_dof)
921 for (
uint i_c=0; i_c<dim; ++i_c) {
922 ArrayDbl &shape_grad_row = ref_shape_grads_expd(i_dof*dim+i_c);
923 for (
uint i_pt=0; i_pt<shape_grad_row.rows(); ++i_pt)
924 shape_grad_row(i_pt) = ref_shape_grads[el_table(3)(i_pt)][i_pt % n_points][i_dof][i_c];
927 auto inv_jac_value = operations[ op.input_ops()[0] ].value<dim, 3>(op_results);
928 for (
uint i_dof=0; i_dof<n_dofs; ++i_dof) {
929 auto shape_grad_value = op.value<3, 1>(op_results, i_dof);
930 Eigen::Map<Eigen::Matrix<ArrayDbl, dim, 1>> ref_shape_grads_dof_value(ref_shape_grads_expd.data() + dim*i_dof, dim, 1);
931 shape_grad_value = inv_jac_value.transpose() * ref_shape_grads_dof_value;
945 auto &result_vec = op.result_matrix();
946 for (
uint i=0;i<result_vec(0,0).data_size(); ++i) {
947 result_vec(0,0)(i) = 1.0;
960 template<
unsigned int spacedim = 3>
982 template<
unsigned int dim>
1014 template<
unsigned int spacedim = 3>
1036 template<
unsigned int dim>
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Outer product - only proposal of multi operator.
ArenaVec< T > get_vec() const
Convert ArenaOVec to ArenaVec and its.
size_t data_size() const
Getter for data_size_.
Class represents element or FE operations.
Eigen::Matrix< ArenaVec< double >, Eigen::Dynamic, Eigen::Dynamic > & result_matrix()
Return map referenced result as Eigen::Matrix.
void allocate_result(size_t data_size, PatchArena &arena)
std::string format_shape() const
std::vector< uint > input_ops_
Indices of operations in PatchPointValues::operations_ vector on which ElOp is depended.
std::vector< uint > set_shape_vec(std::initializer_list< uint > shape) const
Aligns shape_vec to 2 items (equal to matrix number of dimensions)
Eigen::Matrix< ArenaVec< double >, Eigen::Dynamic, Eigen::Dynamic > result_
Result matrix of operation.
std::vector< uint > shape_
Shape of stored data (size of vector or number of rows and cols of matrix)
uint result_row() const
Getter for result_row_.
const std::vector< uint > & shape() const
Getter for shape_.
const std::vector< uint > & input_ops() const
Getter for input_ops_.
OpSizeType size_type_
Type of operation by size of vector (element, point or fixed size)
uint result_row_
First row to scalar, vector or matrix result TODO replace.
const Eigen::Matrix< ArenaVec< double >, Eigen::Dynamic, Eigen::Dynamic > & result_matrix() const
Same as previous but return const reference.
uint n_dofs() const
Getter for n_dofs_.
ReinitFunction reinit_func
Pointer to patch reinit function of element data table specialized by operation.
uint n_dofs_
Number of DOFs of FE operations (or 1 in case of element operations)
uint dim() const
Getter for dimension.
Eigen::Map< Eigen::Matrix< ArrayDbl, dim1, dim2 > > value(TableDbl &op_results, uint i_dof=0) const
Return map referenced Eigen::Matrix of given dimension.
Eigen::Map< Eigen::Vector< double, Eigen::Dynamic > > vector_value(TableDbl &op_results) const
Return map referenced Eigen::Matrix of given dimensions.
ElOp(uint dim, std::initializer_list< uint > shape, uint result_row, ReinitFunction reinit_f, OpSizeType size_type, std::vector< uint > input_ops={}, uint n_dofs=1)
Eigen::Map< Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > > matrix_value(TableDbl &op_results, uint dim1, uint dim2) const
Return map referenced Eigen::Matrix of given dimensions.
OpSizeType size_type() const
Getter for size_type_.
void reinit_function(std::vector< ElOp< spacedim >> &operations, TableDbl &data_table, TableInt &int_table)
Call reinit function on element table if function is defined.
Bulk data specialization, order of item in operations_ vector corresponds to the BulkOps enum.
PatchPointValues(uint dim, uint quad_order, AssemblyArena &asm_arena)
Constructor.
void init()
Initialize operations vector.
Bulk Side specialization, order of item in operations_ vector corresponds to the SideOps enum.
PatchPointValues(uint dim, uint quad_order, AssemblyArena &asm_arena)
Constructor.
void init()
Initialize operations vector.
void print_data_tables(ostream &stream, bool points, bool ints) const
PatchPointValues(uint dim, AssemblyArena &asm_arena)
ElOp< spacedim > & make_expansion(ElOp< spacedim > &el_op, std::initializer_list< uint > shape, ReinitFunction reinit_f)
AssemblyArena & asm_arena_
Reference to global assembly arena of PatchFeValues.
std::vector< OpSizeType > row_sizes_
hold sizes of rows by type of operation
void initialize(uint int_cols)
uint n_rows() const
Getter for n_rows_.
std::vector< uint > points_map_
Map of point patch indices to point_vals_ and int_vals_ tables.
uint register_side_point(uint elem_table_row, uint value_patch_idx, uint elem_idx, uint side_idx)
uint n_rows_
Number of columns of point_vals table.
uint register_bulk_point(uint elem_table_row, uint value_patch_idx, uint elem_idx, uint i_point_on_elem)
uint register_element(arma::mat coords, uint element_patch_idx)
void resize_tables(uint n_elems, uint n_points)
Resize data tables. Method is called before reinit of patch.
Tensor tensor_val(uint result_row, uint point_idx) const
uint register_side(arma::mat elm_coords, arma::mat side_coords)
ElOp< spacedim > & make_fe_op(std::initializer_list< uint > shape, ReinitFunction reinit_f, std::vector< uint > input_ops_vec, uint n_dofs, OpSizeType size_type=pointOp)
std::vector< ElOp< spacedim > > operations_
Vector of all defined operations.
Quadrature * get_quadrature() const
Getter for quadrature.
uint n_elems() const
Getter for n_elems_.
uint i_elem_
Index of registered element in table, helper value used during patch creating.
void reset()
Reset number of columns (points and elements)
uint n_points() const
Getter for n_points_.
ElOp< spacedim > & make_fixed_op(std::initializer_list< uint > shape, ReinitFunction reinit_f)
std::vector< uint > elements_map_
Map of element patch indices to el_vals_ table.
uint dim() const
Getter for dim_.
Tensor tensor_value(uint op_idx, uint point_idx, uint i_dof=0) const
Scalar scalar_value(uint op_idx, uint point_idx, uint i_dof=0) const
void init_finalize(PatchArena *patch_arena)
Vector vector_val(uint result_row, uint point_idx) const
uint n_points_
Number of points in patch.
void print_operations(ostream &stream, uint bulk_side) const
Quadrature * quad_
Quadrature of given dimension and order passed in constructor.
Scalar scalar_val(uint result_row, uint point_idx) const
ElOp< spacedim > & make_new_op(std::initializer_list< uint > shape, ReinitFunction reinit_f, std::vector< uint > input_ops_vec, OpSizeType size_type=pointOp)
Vector vector_value(uint op_idx, uint point_idx, uint i_dof=0) const
PatchArena * patch_arena_
Pointer to global patch arena of PatchFeValues.
uint n_elems_
Number of elements in patch.
Symmetric Gauss-Legendre quadrature formulae on simplices.
Base class for quadrature rules on simplices in arbitrary dimensions.
const std::vector< double > & get_weights() const
Return a reference to the whole array of weights.
static Eigen::Vector< ArrayDbl, dim > normal_vector_array(Eigen::Array< uint, Eigen::Dynamic, 1 > loc_side_idx_array)
Class ElementValues calculates data related to transformation of reference cell to actual cell (Jacob...
ArmaMat< double, N, M > mat
@ opCoords
operations evaluated on quadrature points
@ opInvJac
inverse Jacobian
@ opWeights
weight of quadrature point
@ opJac
Jacobian of element.
@ opElCoords
operations evaluated on elements
@ opJxW
JxW value of quadrature point.
@ opNormalVec
normal vector of quadrature point
@ opSideCoords
operations evaluated on sides
@ opJxW
JxW value of quadrature point.
@ opCoords
operation executed expansion to quadrature point (value of element / side > values on quadrature poin...
@ opElCoords
operations evaluated on elements
@ opElJac
Jacobian of element.
@ opWeights
weight of quadrature point
@ opSideJac
Jacobian of element.
std::function< void(std::vector< ElOp< 3 > > &, TableDbl &, TableInt &)> ReinitFunction
Type for conciseness.
OpSizeType
Distinguishes operations by type and size of output rows.
@ pointOp
operation is evaluated on quadrature points
@ elemOp
operation is evaluated on elements or sides
@ fixedSizeOp
operation has fixed size and it is filled during initialization
Definitions of particular quadrature rules on simplices.
Defines reinit operations on bulk points.
static void elop_jac(std::vector< ElOp< 3 >> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table)
static void ptop_scalar_shape_grads(std::vector< ElOp< 3 >> &operations, std::vector< std::vector< arma::mat > > ref_shape_grads, uint scalar_shape_grads_op_idx)
static void ptop_JxW(std::vector< ElOp< 3 >> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table)
static void ptop_coords(FMT_UNUSED std::vector< ElOp< 3 >> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table)
static void elop_jac_det(std::vector< ElOp< 3 >> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table)
static void ptop_weights(std::vector< ElOp< 3 >> &operations, PatchArena *arena, const std::vector< double > &point_weights)
static void elop_inv_jac(std::vector< ElOp< 3 >> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table)
static void ptop_scalar_shape(std::vector< ElOp< 3 >> &operations, std::vector< std::vector< double > > shape_values, uint scalar_shape_op_idx)
Defines common functionality of reinit operations.
static void expand_data(ElOp< 3 > &op, TableDbl &op_results, TableInt &el_table)
static void op_base(FMT_UNUSED std::vector< ElOp< 3 >> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table)
Defines reinit operations on side points.
static void elop_sd_jac_det(std::vector< ElOp< 3 >> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table)
static void ptop_normal_vec(std::vector< ElOp< 3 >> &operations, TableDbl &op_results, TableInt &el_table)
static void ptop_weights(std::vector< ElOp< 3 >> &operations, PatchArena *arena, const std::vector< double > &point_weights)
static void ptop_scalar_shape(std::vector< ElOp< 3 >> &operations, TableDbl &op_results, TableInt &el_table, std::vector< std::vector< std::vector< double > > > shape_values, uint scalar_shape_op_idx)
static void ptop_JxW(std::vector< ElOp< 3 >> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table)
static void ptop_coords(FMT_UNUSED std::vector< ElOp< 3 >> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table)
static void elop_sd_jac(std::vector< ElOp< 3 >> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table)
static void ptop_scalar_shape_grads(std::vector< ElOp< 3 >> &operations, TableDbl &op_results, TableInt &el_table, std::vector< std::vector< std::vector< arma::mat > > > ref_shape_grads, uint scalar_shape_grads_op_idx)
static void elop_el_jac(std::vector< ElOp< 3 >> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table)
static void elop_el_inv_jac(std::vector< ElOp< 3 >> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table)