From 054d58edb633256bf0768a6b7f80fcfc1d5927cb Mon Sep 17 00:00:00 2001
From: Maximilian Bergbauer <maximilian.bergbauer@tum.de>
Date: Thu, 24 Oct 2024 18:11:24 +0200
Subject: [PATCH] Implement curl in FEPointEvaluation

---
 .../deal.II/matrix_free/fe_point_evaluation.h | 38 +++++++++++++++++++
 1 file changed, 38 insertions(+)

diff --git a/include/deal.II/matrix_free/fe_point_evaluation.h b/include/deal.II/matrix_free/fe_point_evaluation.h
index 3878f076e9..5fd0de9d71 100644
--- a/include/deal.II/matrix_free/fe_point_evaluation.h
+++ b/include/deal.II/matrix_free/fe_point_evaluation.h
@@ -773,6 +773,15 @@ public:
   void
   submit_divergence(const Number &value, const unsigned int point_index);
 
+  /**
+   * Return the curl in real coordinates at the point with index
+   * `point_index` after a call to FEPointEvaluation::evaluate() with
+   * EvaluationFlags::gradients set. This functions only makes sense for a
+   * vector field with dim components and dim > 1.
+   */
+  Tensor<1, (dim == 2 ? 1 : dim), Number>
+  get_curl(const unsigned int point_index) const;
+
   /**
    * Return the Jacobian of the transformation on the current cell with the
    * given point index. Prerequisite: This class needs to be constructed with
@@ -2233,6 +2242,35 @@ FEPointEvaluationBase<n_components_, dim, spacedim, Number>::submit_divergence(
 
 
 
+template <int n_components_, int dim, int spacedim, typename Number>
+Tensor<1, (dim == 2 ? 1 : dim), Number>
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::get_curl(
+  const unsigned int point_index) const
+{
+  static_assert(
+    dim > 1 && n_components == dim,
+    "Only makes sense for a vector field with dim components and dim > 1");
+
+  const Tensor<2, dim, Number>            grad = get_gradient(point_index);
+  Tensor<1, (dim == 2 ? 1 : dim), Number> curl;
+  switch (dim)
+    {
+      case 2:
+        curl[0] = grad[1][0] - grad[0][1];
+        break;
+      case 3:
+        curl[0] = grad[2][1] - grad[1][2];
+        curl[1] = grad[0][2] - grad[2][0];
+        curl[2] = grad[1][0] - grad[0][1];
+        break;
+      default:
+        DEAL_II_NOT_IMPLEMENTED();
+    }
+  return curl;
+}
+
+
+
 template <int n_components_, int dim, int spacedim, typename Number>
 inline DerivativeForm<1, dim, spacedim, Number>
 FEPointEvaluationBase<n_components_, dim, spacedim, Number>::jacobian(
-- 
2.39.5