From 443988d2cdb0ed4c346bb6bfc9b89be0a38cc6b7 Mon Sep 17 00:00:00 2001
From: Francesc Verdugo <f.verdugo.rojano@vu.nl>
Date: Fri, 30 Aug 2024 08:09:37 +0200
Subject: [PATCH] Fix 3d case

---
 src/domain.jl          | 11 ++++++++++-
 test/assembly_tests.jl | 17 +++++++++++++++++
 2 files changed, 27 insertions(+), 1 deletion(-)

diff --git a/src/domain.jl b/src/domain.jl
index c6149ef3..e10548e7 100644
--- a/src/domain.jl
+++ b/src/domain.jl
@@ -992,6 +992,14 @@ function align_field(a::AbstractQuantity,glue::BoundaryGlue)
 end
 
 function inverse_map_impl(f,x0)
+    function pseudo_inverse_if_not_square(J)
+        m,n = size(J)
+        if m != n
+            pinv(J)
+        else
+            inv(J)
+        end
+    end
     function invf(fx)
         x = x0
         tol = 1.0e-12
@@ -999,7 +1007,8 @@ function inverse_map_impl(f,x0)
         niters = 100
         for _ in 1:niters
             J = ForwardDiff.jacobian(f,x)
-            dx = pinv(J)*(fx-f(x))
+            Jinv = pseudo_inverse_if_not_square(J)
+            dx = Jinv*(fx-f(x))
             x += dx
             if norm(dx) < tol
                 return x
diff --git a/test/assembly_tests.jl b/test/assembly_tests.jl
index b71f0bc1..7bc9573f 100644
--- a/test/assembly_tests.jl
+++ b/test/assembly_tests.jl
@@ -254,4 +254,21 @@ el2 = ∫( q->abs2(eh(q)), dΩ) |> sum |> sqrt
 eh1 = ∫( q->∇eh(q)⋅∇eh(q), dΩ) |> sum |> sqrt
 @test el2 < tol
 
+# 3d case
+
+n = 2
+domain = (0,1,0,1,0,1)
+cells = (n,n,n)
+mesh = GT.cartesian_mesh(domain,cells)
+Ω = GT.interior(mesh)
+k = 1
+V = GT.lagrange_space(Ω,k)
+dΩ = GT.measure(Ω,2*k)
+gradient(u) = x->ForwardDiff.gradient(u,x)
+∇(u,x) = GT.call(gradient,u)(x)
+a(u,v) = GT.∫( x->∇(u,x)⋅∇(v,x), dΩ)
+l(v) = 0
+x,A,b = GT.linear_problem(Float64,V,a,l)
+A |> display
+
 end # module