cross product for 3x3 matrices

SciLean/Data/DataArray/Operations.lean

@@ -306,6 +306,23 @@ noncomputable
 instance : Solve R I (I×J) where
   solve A B := A.solve' B
 
+/-- Cross product of two vector. -/
+def cross (x y : R^[3]) : R^[3] :=
+  ⊞[x[1]*y[2] - x[2]*y[1],
+    x[2]*y[0] - x[0]*y[2],
+    x[0]*y[1] - x[1]*y[0]]
+
+/-- Matrix corresponding to taking cross product with `x`  -/
+def crossmatrix (x : R^[3]) : R^[3,3] := Id.run do
+  let mut A : R^[3,3] := 0
+  A[0,2] := x[1]; A[0,1] := - x[2]
+  A[1,0] := x[2]; A[1,2] := - x[0]
+  A[2,1] := x[0]; A[2,0] := - x[1]
+  return A
+
+/-- Takes antisymmetric part of matrix `A` and stacks into a vector. -/
+def antisymmpart (A : R^[3,3]) : R^[3] :=
+  ⊞[0.5 * (A[2,1] - A[1,2]), 0.5 * (A[0,2] - A[2,0]), 0.5 * (A[1,0] - A[0,1])]
 
 /-- Cayley Map: https://en.wikipedia.org/wiki/Cayley_transform#Matrix_map -/
 noncomputable

SciLean/Data/DataArray/Operations/AntiSymmPart.lean

@@ -0,0 +1,34 @@
+import SciLean.Data.DataArray.Operations.Diag
+
+namespace SciLean
+
+
+open DataArrayN
+
+def_fun_prop antisymmpart in A
+  with_transitive
+  [RealScalar R] : IsContinuousLinearMap R by sorry_proof
+
+#generate_linear_map_simps DataArrayN.antisymmpart.arg_A.IsLinearMap_rule
+
+abbrev_fun_trans antisymmpart in A : fderiv R by
+  fun_trans
+
+abbrev_fun_trans antisymmpart in A : fwdFDeriv R by
+  autodiff
+
+abbrev_fun_trans antisymmpart in A : adjoint R by
+  equals (fun x => x.crossmatrix) =>
+    sorry_proof
+
+abbrev_fun_trans antisymmpart in A : revFDeriv R by
+  unfold revFDeriv
+  autodiff
+
+abbrev_fun_trans antisymmpart in A : revFDerivProj R Unit by
+  unfold revFDerivProj
+  autodiff
+
+abbrev_fun_trans antisymmpart in A : revFDerivProjUpdate R Unit by
+  unfold revFDerivProjUpdate
+  autodiff

SciLean/Data/DataArray/Operations/Cross.lean

@@ -0,0 +1,55 @@
+import SciLean.Data.DataArray.Operations.Diag
+
+namespace SciLean
+
+variable
+  {I : Type*} [IndexType I] [DecidableEq I]
+  {J : Type*} [IndexType J] [DecidableEq J]
+  {R : Type*} [RealScalar R] [PlainDataType R]
+
+
+open DataArrayN
+
+def_fun_prop cross in x    with_transitive : IsContinuousLinearMap R by sorry_proof
+def_fun_prop cross in y    with_transitive : IsContinuousLinearMap R by sorry_proof
+def_fun_prop cross in x y  with_transitive : Differentiable R by sorry_proof
+
+#generate_linear_map_simps DataArrayN.cross.arg_x.IsLinearMap_rule
+#generate_linear_map_simps DataArrayN.cross.arg_y.IsLinearMap_rule
+
+
+abbrev_fun_trans cross in x y : fderiv R by
+  rw[fderiv_wrt_prod (by fun_prop)]
+  fun_trans
+
+abbrev_fun_trans cross in x y : fwdFDeriv R by
+  rw[fwdFDeriv_wrt_prod (by fun_prop)]
+  autodiff
+
+abbrev_fun_trans cross in x : adjoint R by
+  equals (fun z => y.cross z) =>
+    sorry_proof
+
+abbrev_fun_trans cross in y : adjoint R by
+  equals (fun z => z.cross x) =>
+    funext z
+    apply AdjointSpace.ext_inner_left R
+    intro y
+    rw[← adjoint_ex _ (by fun_prop)]
+    simp[DataArrayN.inner_def, DataArrayN.cross,
+         sum_over_prod, sum_pull]
+    -- expand the sum explicitely and call ring
+    sorry_proof
+
+abbrev_fun_trans cross in x y : revFDeriv R by
+  rw[revFDeriv_wrt_prod (by fun_prop)]
+  unfold revFDeriv
+  autodiff
+
+abbrev_fun_trans cross in x y : revFDerivProj R Unit by
+  unfold revFDerivProj
+  autodiff
+
+abbrev_fun_trans cross in x y : revFDerivProjUpdate R Unit by
+  unfold revFDerivProjUpdate
+  autodiff

SciLean/Data/DataArray/Operations/CrossMatrix.lean

@@ -0,0 +1,30 @@
+import SciLean.Data.DataArray.Operations.Multiply
+
+namespace SciLean
+
+open DataArrayN
+
+def_fun_prop crossmatrix in x
+  with_transitive : IsContinuousLinearMap R by sorry_proof
+
+#generate_linear_map_simps DataArrayN.crossmatrix.arg_x.IsLinearMap_rule
+
+abbrev_fun_trans crossmatrix in x : fderiv R by autodiff
+
+abbrev_fun_trans crossmatrix in x : fwdFDeriv R by autodiff
+
+abbrev_fun_trans crossmatrix in x : adjoint R by
+  equals (fun A => A.antisymmpart) =>
+    sorry_proof
+
+abbrev_fun_trans crossmatrix in x : revFDeriv R by
+  unfold revFDeriv
+  autodiff
+
+abbrev_fun_trans crossmatrix in x : revFDerivProj R Unit by
+  unfold revFDerivProj
+  autodiff
+
+abbrev_fun_trans crossmatrix in x : revFDerivProjUpdate R Unit by
+  unfold revFDerivProjUpdate
+  autodiff

SciLean/Data/DataArray/Operations/Simps.lean

@@ -180,3 +180,11 @@ theorem det_inv_eq_inv_det (A : R^[I,I]) :
 
 
 end
+
+section CrossProduct
+
+@[simp, simp_core]
+theorem cossmatrix_antisymmpart (x : R^[3]) :
+  x.crossmatrix.antisymmpart = x := by sorry_proof
+
+end CrossProduct