From 70228a3cc1cfde5a357aa1fddc7792cc72a8616e Mon Sep 17 00:00:00 2001 From: Jan Mandel Date: Wed, 21 Feb 2024 12:46:19 -0700 Subject: [PATCH] separated constant part in module_w_assembly.f90 --- femwind/fortran/module_w_assembly.f90 | 37 ++++++++++++++++++++--------------- 1 file changed, 21 insertions(+), 16 deletions(-) diff --git a/femwind/fortran/module_w_assembly.f90 b/femwind/fortran/module_w_assembly.f90 index d9e8559..de47fa5 100644 --- a/femwind/fortran/module_w_assembly.f90 +++ b/femwind/fortran/module_w_assembly.f90 @@ -22,7 +22,6 @@ subroutine w_assembly( & implicit none - !*** arguments integer, intent(in):: & @@ -31,28 +30,27 @@ integer, intent(in):: & ifps, ifpe, kfps, kfpe, jfps, jfpe, & ! fire patch bounds ifts, ifte, kfts, kfte, jfts,jfte - - - - real, intent(in), dimension(ifms:ifme,kfms:kfme,jfms:jfme):: lambda -!Input for hexa - +!Input for hexa real, intent(in) :: A(3,3) real, intent(in), dimension(ifms:ifme, kfms:kfme, jfms:jfme)::X,Y,Z,u0, v0, w0 +!Output real, intent(out), dimension(ifms:ifme, kfms:kfme, jfms:jfme)::U,V,W + !*** local -integer:: ie1, ie2, ie3, ic1, ic2, ic3, iloc, i, & - k1, k2, k3 +integer:: ie1, ie2, ie3, ic1, ic2, ic3, iloc, i, k1, k2, k3 real :: Kloc(8,8), Floc(8), Jg(8,3) real :: Xloc(3,8), u0loc(3) real :: grad(3) real :: lambda_loc(8) real :: A_inv(3,3) real :: vol + +!*** executable + lambda_loc = 0. Xloc = 99999. Jg = 0. @@ -60,8 +58,7 @@ Kloc = 0. Floc = 0. grad = 0. u0loc =0. - - + !*** u0loc is an input for module_hexa, but is not used to construct K. Do I need to define this? !** executable @@ -69,17 +66,13 @@ u0loc =0. do ie2=jfts,jfte do ie3=kfts, kfte do ie1=ifts, ifte + ! constant part do ic2=0,1 do ic3=0,1 do ic1=0,1 - iloc=1+ic1+2*(ic2+2*ic3); !local index of the node in the element Xloc(1,iloc)=X(ie1 + ic1, ie3 + ic3, ie2 + ic2) Xloc(2,iloc)=Y(ie1 + ic1, ie3 + ic3, ie2 + ic2) Xloc(3,iloc)=Z(ie1 + ic1, ie3 + ic3, ie2 + ic2) - k1 = ie1 + ic1 - k2 = ie2 + ic2 - k3 = ie3 + ic3 - lambda_loc(iloc) = lambda(k1,k3,k2) enddo enddo enddo @@ -92,6 +85,18 @@ do ie2=jfts,jfte call hexa(A,Xloc,u0loc,Kloc,Floc,Jg,vol,3) !print*, 'Jg is', Jg !print*, shape(jg) + !*** end of constant part + do ic2=0,1 + do ic3=0,1 + do ic1=0,1 + iloc=1+ic1+2*(ic2+2*ic3); !local index of the node in the element + k1 = ie1 + ic1 + k2 = ie2 + ic2 + k3 = ie3 + ic3 + lambda_loc(iloc) = lambda(k1,k3,k2) + enddo + enddo + enddo grad = matmul(transpose(Jg),lambda_loc) !not ok print *,'Grad before multiplication by A_inv is', grad -- 2.11.4.GIT