From f541d08b4343132f8d872b4dc00ae39258288e6f Mon Sep 17 00:00:00 2001
From: Jan Mandel <jan.mandel@gmail.com>
Date: Sat, 24 Feb 2024 18:48:54 -0700
Subject: [PATCH] taking computing of the divergence load F outside of hexa

---
 femwind/fortran/module_f_assembly.f90 |  1 +
 femwind/fortran/module_hexa.f90       | 12 ------------
 femwind/hexa.m                        | 12 ++++++++----
 femwind/hexa_test.m                   |  4 +++-
 4 files changed, 12 insertions(+), 17 deletions(-)

diff --git a/femwind/fortran/module_f_assembly.f90 b/femwind/fortran/module_f_assembly.f90
index d78824e..27402e6 100644
--- a/femwind/fortran/module_f_assembly.f90
+++ b/femwind/fortran/module_f_assembly.f90
@@ -74,6 +74,7 @@ do ie2=jfts,jfte
             u0loc(3) = Zu0(ie1,ie3,ie2)
             
             call hexa(A,Xloc,u0loc,Kloc,Floc,Jg,vol,2)
+            Floc = -vol*matmul(Jg,u0loc)
             do id2 = 0,1
                 do id3 = 0,1
                     do id1 = 0,1
diff --git a/femwind/fortran/module_hexa.f90 b/femwind/fortran/module_hexa.f90
index 952f7b4..274dde9 100644
--- a/femwind/fortran/module_hexa.f90
+++ b/femwind/fortran/module_hexa.f90
@@ -95,19 +95,7 @@ do i=1,9          ! loop over i quadrature nodes + center
         ! exact if the element is linear image of reference 
         ! i.e. all sides planar, and opposite sides are parallel 
         vol = abs(detJx)*8
-        if(p) print *, 'hexa: vol',vol
-        ! Floc =  - Jg * u0 * vol;
-        !u0_tmp = u0*vol
-        !do j = 1,8
-        !    tmp = 0
-        !    do k = 1,3
-        !        tmp = tmp+Jg(j,k)*u0_tmp(k)
-        !    end do
-        !    tmp_mat(j) = tmp        
-        !end do
-        !Floc = Floc-tmp_mat
         Floc = vol*matmul(-Jg,u0)
-        !print *,'hexa: Floc err=',maxval(abs(Floc  + tmp_mat))
     end if !end for computing Floc
 
 end do ! loop over i quadrature nodes + center
diff --git a/femwind/hexa.m b/femwind/hexa.m
index bc92409..aa220a2 100644
--- a/femwind/hexa.m
+++ b/femwind/hexa.m
@@ -1,4 +1,4 @@
-function [Kloc,Floc,Jg]=hexa(A,X,u0)
+function [Kloc,Floc,Jg,vol]=hexa(A,X,u0)
 % create local stiffness matrix for hexa 3d element
 % in:
 %   A   coefficient matrix size 3x3, symmetric positive definite
@@ -7,7 +7,8 @@ function [Kloc,Floc,Jg]=hexa(A,X,u0)
 % out:
 %   Kloc   local stiffness matrix
 %   Floc   local divergence load vector
-%   Jg     gradient at center of function with values V is V'*Jg          
+%   Jg     gradient at centoer of function with values V is V'*Jg          
+%   vol    optional, volume of the element
 
 % basis functions on reference element [-1,1]^3
 Nb = 8;  % number of basis functions
@@ -44,8 +45,11 @@ for j=1:Ng+1
         % volume = determinant at midpoint times the volume of reference element
         % exact if the element is linear image of reference 
         % or all sides planar, and sides in two directions are parallel 
-        vol = adetJx*8;   
-        Floc = Floc - Jg * u0 * vol;
+        vol_tmp = adetJx*8;   
+        Floc = Floc - Jg * u0 * vol_tmp;
+        if nargout > 3
+            vol = vol_tmp
+        end
     end
 end
 %check_symmetry(Kloc,'Kloc',eps)
diff --git a/femwind/hexa_test.m b/femwind/hexa_test.m
index be48318..b8e6334 100644
--- a/femwind/hexa_test.m
+++ b/femwind/hexa_test.m
@@ -27,7 +27,9 @@ T = rand(3);
 X = T*X0+rand(3,1)*ones(1,8);
 
 u = rand(3,1);
-[K, F, Jg] = hexa(A,X,u);
+[K, Fold, Jg, vol] = hexa(A,X,u);
+F = -vol*Jg*u;
+err_F=big(F-Fold)
 eig_K=eig(K)  % one eigenvalue zero
 
 % test F
-- 
2.11.4.GIT