From f2b5425a879cd879fc1299783811d78770a1e781 Mon Sep 17 00:00:00 2001
From: Douglas Guptill <douglas.guptill@dal.ca>
Date: Mon, 13 Jul 2009 14:04:00 +0000
Subject: [PATCH] incorporate updates (2009-07-10) from Jean Braun
---
DOUAR.f90 | 207 +++++++++++++++++++--------------------
Makefile | 3 +-
NN/volume.c | 63 ++++++------
OCTREE/OctreeBitPlus.f90 | 39 ++++++++
create_surfaces.f90 | 8 +-
improve_osolve.f90 | 39 ++++----
wsmp_setup.f90 | 2 +-
7 files changed, 200 insertions(+), 161 deletions(-)
diff --git a/DOUAR.f90 b/DOUAR.f90
index f28dcbf8..4d18bca2 100644
--- a/DOUAR.f90
+++ b/DOUAR.f90
@@ -308,22 +308,22 @@ do while (istep.le.params%nstep)
!---------------------------------------------------------------------------|
do is=1,params%ns
surface0(is)%nsurface=surface(is)%nsurface
- allocate (surface0(is)%r(surface0(is)%nsurface),stat=threadinfo%err)
- call heap (threadinfo,'surface0(is)%r','main',size(surface0(is)%r),'dp',+1)
- allocate (surface0(is)%s(surface0(is)%nsurface),stat=threadinfo%err)
- call heap (threadinfo,'surface0(is)%s','main',size(surface0(is)%s),'dp',+1)
- allocate (surface0(is)%x(surface0(is)%nsurface),stat=threadinfo%err)
- call heap (threadinfo,'surface0(is)%x','main',size(surface0(is)%x),'dp',+1)
- allocate (surface0(is)%y(surface0(is)%nsurface),stat=threadinfo%err)
- call heap (threadinfo,'surface0(is)%y','main',size(surface0(is)%y),'dp',+1)
- allocate (surface0(is)%z(surface0(is)%nsurface),stat=threadinfo%err)
- call heap (threadinfo,'surface0(is)%z','main',size(surface0(is)%z),'dp',+1)
- allocate (surface0(is)%xn(surface0(is)%nsurface),stat=threadinfo%err)
- call heap (threadinfo,'surface0(is)%xn','main',size(surface0(is)%xn),'dp',+1)
- allocate (surface0(is)%yn(surface0(is)%nsurface),stat=threadinfo%err)
- call heap (threadinfo,'surface0(is)%yn','main',size(surface0(is)%yn),'dp',+1)
- allocate (surface0(is)%zn(surface0(is)%nsurface),stat=threadinfo%err)
- call heap (threadinfo,'surface0(is)%zn','main',size(surface0(is)%zn),'dp',+1)
+ allocate (surface0(is)%r(surface0(is)%nsurface),stat=threadinfo%err)
+ call heap (threadinfo,'surface0(is)%r','main',size(surface0(is)%r),'dp',+1)
+ allocate (surface0(is)%s(surface0(is)%nsurface),stat=threadinfo%err)
+ call heap (threadinfo,'surface0(is)%s','main',size(surface0(is)%s),'dp',+1)
+ allocate (surface0(is)%x(surface0(is)%nsurface),stat=threadinfo%err)
+ call heap (threadinfo,'surface0(is)%x','main',size(surface0(is)%x),'dp',+1)
+ allocate (surface0(is)%y(surface0(is)%nsurface),stat=threadinfo%err)
+ call heap (threadinfo,'surface0(is)%y','main',size(surface0(is)%y),'dp',+1)
+ allocate (surface0(is)%z(surface0(is)%nsurface),stat=threadinfo%err)
+ call heap (threadinfo,'surface0(is)%z','main',size(surface0(is)%z),'dp',+1)
+ allocate (surface0(is)%xn(surface0(is)%nsurface),stat=threadinfo%err)
+ call heap (threadinfo,'surface0(is)%xn','main',size(surface0(is)%xn),'dp',+1)
+ allocate (surface0(is)%yn(surface0(is)%nsurface),stat=threadinfo%err)
+ call heap (threadinfo,'surface0(is)%yn','main',size(surface0(is)%yn),'dp',+1)
+ allocate (surface0(is)%zn(surface0(is)%nsurface),stat=threadinfo%err)
+ call heap (threadinfo,'surface0(is)%zn','main',size(surface0(is)%zn),'dp',+1)
surface0(is)%r=surface(is)%r
surface0(is)%s=surface(is)%s
surface0(is)%x=surface(is)%x
@@ -360,8 +360,7 @@ do while (istep.le.params%nstep)
write(*,'(a,L1)') shift//'(1) params%griditer < 0 ',(params%griditer.lt.0)
write(*,'(a,L1)') shift//'(2) current_level==params%levelmax_oct',(current_level==params%levelmax_oct)
write(*,'(a,L1)') shift//'(3) increase_current_level ',increase_current_level
- write(*,'(a,L1)') shift//'converge = (1) & (2) & (3) -> ',&
- (params%griditer .lt. 0 .and. current_level==params%levelmax_oct .and. increase_current_level)
+ write(*,'(a,L1)') shift//'converge = (1) & (2) & (3) -> ',(params%griditer .lt. 0 .and. current_level==params%levelmax_oct .and. increase_current_level)
end if
@@ -383,8 +382,8 @@ do while (istep.le.params%nstep)
osolve%noctree=params%noctreemax
- allocate (osolve%octree(osolve%noctree),stat=threadinfo%err)
- call heap (threadinfo,'osolve%octree','main',size(osolve%octree),'int',+1)
+ allocate (osolve%octree(osolve%noctree),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%octree','main',size(osolve%octree),'int',+1)
call octree_init (osolve%octree,osolve%noctree)
call octree_create_uniform (osolve%octree,osolve%noctree,params%leveluniform_oct)
@@ -412,16 +411,16 @@ do while (istep.le.params%nstep)
if (iproc.eq.0) write (8,*) osolve%nleaves,' leaves in solve octree '
- allocate (osolve%icon(8,osolve%nleaves),stat=threadinfo%err)
- call heap (threadinfo,'osolve%icon','main',size(osolve%icon),'int',+1)
- allocate (osolve%x(osolve%nnode),stat=threadinfo%err)
- call heap (threadinfo,'osolve%x', 'main',size(osolve%x),'dp',+1)
- allocate (osolve%y(osolve%nnode),stat=threadinfo%err)
- call heap (threadinfo,'osolve%y', 'main',size(osolve%y),'dp',+1)
- allocate (osolve%z(osolve%nnode),stat=threadinfo%err)
- call heap (threadinfo,'osolve%z', 'main',size(osolve%z),'dp',+1)
- allocate (osolve%lsf(osolve%nnode,osolve%nlsf),stat=threadinfo%err)
- call heap (threadinfo,'osolve%lsf', 'main',size(osolve%lsf),'dp',+1)
+ allocate (osolve%icon(8,osolve%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%icon','main',size(osolve%icon),'int',+1)
+ allocate (osolve%x(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%x', 'main',size(osolve%x),'dp',+1)
+ allocate (osolve%y(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%y', 'main',size(osolve%y),'dp',+1)
+ allocate (osolve%z(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%z', 'main',size(osolve%z),'dp',+1)
+ allocate (osolve%lsf(osolve%nnode,osolve%nlsf),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%lsf', 'main',size(osolve%lsf),'dp',+1)
osolve%lsf=0.d0
call octree_find_node_connectivity (osolve%octree,osolve%noctree, &
@@ -475,20 +474,20 @@ do while (istep.le.params%nstep)
! allocate fields of osolve
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
- allocate (osolve%u(osolve%nnode),stat=threadinfo%err)
- call heap (threadinfo,'osolve%u','main',size(osolve%u),'dp',+1)
- allocate (osolve%v(osolve%nnode),stat=threadinfo%err)
- call heap (threadinfo,'osolve%v','main',size(osolve%v),'dp',+1)
- allocate (osolve%w(osolve%nnode),stat=threadinfo%err)
- call heap (threadinfo,'osolve%w','main',size(osolve%w),'dp',+1)
- allocate (osolve%temp(osolve%nnode),stat=threadinfo%err)
- call heap (threadinfo,'osolve%temp','main',size(osolve%temp),'dp',+1)
- allocate (osolve%pressure(osolve%nleaves),stat=threadinfo%err)
- call heap (threadinfo,'osolve%pressure','main',size(osolve%pressure),'dp',+1)
- allocate (osolve%spressure(osolve%nleaves),stat=threadinfo%err)
- call heap (threadinfo,'osolve%spressure','main',size(osolve%spressure),'dp',+1)
- allocate (osolve%strain(osolve%nnode),stat=threadinfo%err)
- call heap (threadinfo,'osolve%strain','main',size(osolve%strain),'dp',+1)
+ allocate (osolve%u(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%u','main',size(osolve%u),'dp',+1)
+ allocate (osolve%v(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%v','main',size(osolve%v),'dp',+1)
+ allocate (osolve%w(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%w','main',size(osolve%w),'dp',+1)
+ allocate (osolve%temp(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%temp','main',size(osolve%temp),'dp',+1)
+ allocate (osolve%pressure(osolve%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%pressure','main',size(osolve%pressure),'dp',+1)
+ allocate (osolve%spressure(osolve%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%spressure','main',size(osolve%spressure),'dp',+1)
+ allocate (osolve%strain(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%strain','main',size(osolve%strain),'dp',+1)
osolve%u = 0.d0
osolve%v = 0.d0
osolve%w = 0.d0
@@ -527,8 +526,8 @@ do while (istep.le.params%nstep)
call show_time (total,step,inc,1,'Find bad faces$')
osolve%nface=osolve%nleaves
- allocate (osolve%iface(9,osolve%nface),stat=threadinfo%err)
- call heap (threadinfo,'osolve%iface','main',size(osolve%iface),'int',+1)
+ allocate (osolve%iface(9,osolve%nface),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%iface','main',size(osolve%iface),'int',+1)
call octree_find_bad_faces (osolve%octree,osolve%noctree, &
osolve%iface,osolve%nface, &
@@ -546,18 +545,18 @@ do while (istep.le.params%nstep)
!------------------------------------------------------------------------|
call show_time (total,step,inc,1,'Find void$')
- allocate (vo%node(osolve%nnode),stat=threadinfo%err)
- call heap (threadinfo,'vo%node','main',size(vo%node),'int',+1)
- allocate (vo%leaf(osolve%nnode),stat=threadinfo%err)
- call heap (threadinfo,'vo%leaf','main',size(vo%leaf),'int',+1)
- allocate (vo%face(osolve%nface),stat=threadinfo%err)
- call heap (threadinfo,'vo%face','main',size(vo%face),'int',+1)
- allocate (vo%rtf(osolve%nnode),stat=threadinfo%err)
- call heap (threadinfo,'vo%rtf','main',size(vo%rtf),'int',+1)
- allocate (vo%ftr(osolve%nnode),stat=threadinfo%err)
- call heap (threadinfo,'vo%ftr','main',size(vo%ftr),'int',+1)
- allocate (vo%influid(osolve%nnode),stat=threadinfo%err)
- call heap (threadinfo,'vo%influid','main',size(vo%influid),'bool',+1)
+ allocate (vo%node(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'vo%node','main',size(vo%node),'int',+1)
+ allocate (vo%leaf(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'vo%leaf','main',size(vo%leaf),'int',+1)
+ allocate (vo%face(osolve%nface),stat=threadinfo%err)
+ call heap (threadinfo,'vo%face','main',size(vo%face),'int',+1)
+ allocate (vo%rtf(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'vo%rtf','main',size(vo%rtf),'int',+1)
+ allocate (vo%ftr(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'vo%ftr','main',size(vo%ftr),'int',+1)
+ allocate (vo%influid(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'vo%influid','main',size(vo%influid),'bool',+1)
call find_void_nodes (params,osolve,vo,istep,ov,threadinfo)
@@ -567,10 +566,10 @@ do while (istep.le.params%nstep)
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
call show_time (total,step,inc,1,'Define boundary conditions$')
- allocate (osolve%kfix(osolve%nnode*3),stat=threadinfo%err)
- call heap (threadinfo,'osolve%kfix','main',size(osolve%kfix),'int',+1)
- allocate (osolve%kfixt(osolve%nnode),stat=threadinfo%err)
- call heap (threadinfo,'osolve%kfixt','main',size(osolve%kfixt),'int',+1)
+ allocate (osolve%kfix(osolve%nnode*3),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%kfix','main',size(osolve%kfix),'int',+1)
+ allocate (osolve%kfixt(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%kfixt','main',size(osolve%kfixt),'int',+1)
call define_bc (params,osolve,vo)
@@ -619,22 +618,22 @@ do while (istep.le.params%nstep)
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
- allocate (osolve%eviscosity(osolve%nleaves),stat=threadinfo%err)
- call heap (threadinfo,'osolve%eviscosity','main',size(osolve%eviscosity),'dp',+1)
- allocate (osolve%q(osolve%nleaves),stat=threadinfo%err)
- call heap (threadinfo,'osolve%q','main',size(osolve%q),'dp',+1)
- allocate (osolve%crit(osolve%nleaves),stat=threadinfo%err)
- call heap (threadinfo,'osolve%crit','main',size(osolve%crit),'dp',+1)
- allocate (osolve%e2d(osolve%nleaves),stat=threadinfo%err)
- call heap (threadinfo,'osolve%e2d','main',size(osolve%e2d),'dp',+1)
- allocate (osolve%e3d(osolve%nleaves),stat=threadinfo%err)
- call heap (threadinfo,'osolve%e3d','main',size(osolve%e3d),'dp',+1)
- allocate (osolve%lode(osolve%nleaves),stat=threadinfo%err)
- call heap (threadinfo,'osolve%lode','main',size(osolve%lode),'dp',+1)
- allocate (osolve%is_plastic(osolve%nleaves),stat=threadinfo%err)
- call heap (threadinfo,'osolve%is_plastic','main',size(osolve%is_plastic),'bool',+1)
- allocate (weightel(osolve%nleaves))
- call heap (threadinfo,'weightel','main',size(weightel),'bool',+1)
+ allocate (osolve%eviscosity(osolve%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%eviscosity','main',size(osolve%eviscosity),'dp',+1)
+ allocate (osolve%q(osolve%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%q','main',size(osolve%q),'dp',+1)
+ allocate (osolve%crit(osolve%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%crit','main',size(osolve%crit),'dp',+1)
+ allocate (osolve%e2d(osolve%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%e2d','main',size(osolve%e2d),'dp',+1)
+ allocate (osolve%e3d(osolve%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%e3d','main',size(osolve%e3d),'dp',+1)
+ allocate (osolve%lode(osolve%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%lode','main',size(osolve%lode),'dp',+1)
+ allocate (osolve%is_plastic(osolve%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%is_plastic','main',size(osolve%is_plastic),'bool',+1)
+ allocate (weightel(osolve%nleaves))
+ call heap (threadinfo,'weightel','main',size(weightel),'bool',+1)
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
@@ -770,8 +769,8 @@ do while (istep.le.params%nstep)
call show_time (total,step,inc,1,'slicing the cube$')
! call slices(params,osolve,ov,vo,sections,istep,iter,iter_nl)
- call heap (threadinfo,'osolve%eviscosity','main',size(osolve%eviscosity),'dp',-1) ; deallocate (osolve%eviscosity)
- call heap (threadinfo,'osolve%q','main',size(osolve%q),'dp',-1) ; deallocate (osolve%q)
+ call heap (threadinfo,'osolve%eviscosity','main',size(osolve%eviscosity),'dp',-1) ; deallocate (osolve%eviscosity)
+ call heap (threadinfo,'osolve%q','main',size(osolve%q),'dp',-1) ; deallocate (osolve%q)
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! compute isostasy
@@ -823,7 +822,7 @@ do while (istep.le.params%nstep)
if (.not.converge) then
call heap (threadinfo,'osolve%octree','main',size(osolve%octree),'int',-1) ; deallocate (osolve%octree)
- call heap (threadinfo,'osolve%icon','main',size(osolve%icon),'int',-1) ; deallocate (osolve%icon)
+ call heap (threadinfo,'osolve%icon','main',size(osolve%icon),'int',-1) ; deallocate (osolve%icon)
call heap (threadinfo,'osolve%x','main',size(osolve%x),'dp',-1) ; deallocate (osolve%x)
call heap (threadinfo,'osolve%y','main',size(osolve%y),'dp',-1) ; deallocate (osolve%y)
call heap (threadinfo,'osolve%z','main',size(osolve%z),'dp',-1) ; deallocate (osolve%z)
@@ -897,18 +896,18 @@ do while (istep.le.params%nstep)
olsf%nleaves=ov%nleaves
olsf%noctree=ov%noctree
- allocate (olsf%octree(olsf%noctree),stat=threadinfo%err)
- call heap (threadinfo,'olsf%octree','main',size(olsf%octree),'int',+1)
- allocate (olsf%x(olsf%nnode),stat=threadinfo%err)
- call heap (threadinfo,'olsf%x','main',size(olsf%x),'dp',+1)
- allocate (olsf%y(olsf%nnode),stat=threadinfo%err)
- call heap (threadinfo,'olsf%y','main',size(olsf%y),'dp',+1)
- allocate (olsf%z(olsf%nnode),stat=threadinfo%err)
- call heap (threadinfo,'olsf%z','main',size(olsf%z),'dp',+1)
- allocate (olsf%lsf(olsf%nnode),stat=threadinfo%err)
- call heap (threadinfo,'olsf%lsf','main',size(olsf%lsf),'dp',+1)
- allocate (olsf%icon(params%mpe,olsf%nleaves),stat=threadinfo%err)
- call heap (threadinfo,'olsf%icon','main',size(olsf%icon),'int',+1)
+ allocate (olsf%octree(olsf%noctree),stat=threadinfo%err)
+ call heap (threadinfo,'olsf%octree','main',size(olsf%octree),'int',+1)
+ allocate (olsf%x(olsf%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'olsf%x','main',size(olsf%x),'dp',+1)
+ allocate (olsf%y(olsf%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'olsf%y','main',size(olsf%y),'dp',+1)
+ allocate (olsf%z(olsf%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'olsf%z','main',size(olsf%z),'dp',+1)
+ allocate (olsf%lsf(olsf%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'olsf%lsf','main',size(olsf%lsf),'dp',+1)
+ allocate (olsf%icon(params%mpe,olsf%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'olsf%icon','main',size(olsf%icon),'int',+1)
olsf%octree(1:olsf%noctree)=ov%octree(1:ov%noctree)
olsf%x(1:olsf%nnode)=ov%x(1:ov%nnode)
@@ -1124,13 +1123,13 @@ do while (istep.le.params%nstep)
olsf%nleaves=ov%nleaves
olsf%noctree=ov%noctree
- allocate (olsf%octree(olsf%noctree),stat=threadinfo%err) ; call heap(threadinfo,'olsf%octree','main',size(olsf%octree),'int',+1)
- allocate (olsf%x(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%x','main',size(olsf%x),'dp',+1)
- allocate (olsf%y(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%y','main',size(olsf%y),'dp',+1)
- allocate (olsf%z(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%z','main',size(olsf%z),'dp',+1)
- allocate (olsf%lsf(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%lsf','main',size(olsf%lsf),'dp',+1)
- allocate (olsf%icon(params%mpe,olsf%nleaves),stat=threadinfo%err)
- call heap (threadinfo,'olsf%icon','main',size(olsf%icon),'int',+1)
+ allocate (olsf%octree(olsf%noctree),stat=threadinfo%err) ; call heap (threadinfo,'olsf%octree','main',size(olsf%octree),'int',+1)
+ allocate (olsf%x(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%x','main',size(olsf%x),'dp',+1)
+ allocate (olsf%y(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%y','main',size(olsf%y),'dp',+1)
+ allocate (olsf%z(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%z','main',size(olsf%z),'dp',+1)
+ allocate (olsf%lsf(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%lsf','main',size(olsf%lsf),'dp',+1)
+ allocate (olsf%icon(params%mpe,olsf%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'olsf%icon','main',size(olsf%icon),'int',+1)
olsf%octree(1:olsf%noctree)=ov%octree(1:ov%noctree)
olsf%x(1:olsf%nnode)=ov%x(1:ov%nnode)
@@ -1304,10 +1303,10 @@ call heap (threadinfo,'ov%vnode','main',size(ov%vnode),'dp',-1) ; deallocat
call heap (threadinfo,'ov%wnode','main',size(ov%wnode),'dp',-1) ; deallocate (ov%wnode)
call heap (threadinfo,'ov%icon','main',size(ov%icon),'int',-1) ; deallocate (ov%icon)
call heap (threadinfo,'ov%temp','main',size(ov%temp),'dp',-1) ; deallocate (ov%temp)
-call heap (threadinfo,'ov%temporary_nodal...','main',size(ov%temporary_nodal_pressure),'dp',-1)
-deallocate (ov%temporary_nodal_pressure)
-call heap (threadinfo,'ov%whole_leaf...','main',size(ov%whole_leaf_in_fluid),'bool',-1)
-deallocate (ov%whole_leaf_in_fluid)
+call heap (threadinfo,'ov%temporary_nodal...','main',size(ov%temporary_nodal_pressure),'dp',-1)
+ deallocate (ov%temporary_nodal_pressure)
+call heap (threadinfo,'ov%whole_leaf...','main',size(ov%whole_leaf_in_fluid),'bool',-1)
+ deallocate (ov%whole_leaf_in_fluid)
do is=1,params%ns
diff --git a/Makefile b/Makefile
index b7030b5b..1f90b23b 100644
--- a/Makefile
+++ b/Makefile
@@ -63,11 +63,10 @@ refine_surface.o \
solve_with_pwssmp.o \
solve_with_pwgsmp.o \
strain_history.o \
-scanfile.o slices.o smooth_pressures.o \
+scanfile.o smooth_pressures.o \
toolbox.o \
update_cloud_structure.o \
update_cloud_fields.o \
-visualise_matrix.o \
remove_point.o \
wsmp_setup.o \
DOUAR.o \
diff --git a/NN/volume.c b/NN/volume.c
index e827973f..ba0ac82f 100644
--- a/NN/volume.c
+++ b/NN/volume.c
@@ -115,7 +115,7 @@ if (nn == 1)
}
else if ( *(b+l) < 0.) /* Constraints are inconsistent */
{
- printf("Inconsistent constraints found after reduction to n = 1 \n");
+/* printf("Inconsistent constraints found after reduction to n = 1 \n");*/
return(0.);
}
}
@@ -167,8 +167,8 @@ for (i=0;i<mm;i++)
/* allocate memory */
- ap = (float *) malloc((sizeof (float))*nm1*mm1);
- bp = (float *) malloc((sizeof (float))*mm1);
+ ap = (float *) malloc((long)nm1*(long)mm1*sizeof *ap);
+ bp = (float *) malloc((long)mm1*sizeof *bp);
/* reduce a and b into ap and bp eliminating variable t and constraint i */
@@ -221,8 +221,8 @@ void volume_();
void volume();
void volume_(a,b,m,n,mmax,nmax,result)
- int *n, *m, *mmax, *nmax;
- float *a, *b;
+int *n, *m, *mmax, *nmax;
+float *a, *b;
float *result;
{
volume (a,b,m,n,mmax,nmax,result);
@@ -317,7 +317,7 @@ if (nn == 1)
/* Constraints are inconsistent.
Set volume and derivative to zero. */
- printf("Inconsistent constraints found after reduction to n = 1 \n");
+/* printf("Inconsistent constraints found after reduction to n = 1 \n");*/
*dvdb = 0.;
return(0.);
}
@@ -377,8 +377,8 @@ for (i=0;i<mm;i++)
/* allocate memory */
- ap = (float *) malloc((sizeof (float))*nm1*mm1);
- bp = (float *) malloc((sizeof (float))*mm1);
+ ap = (float *) malloc((long)nm1*(long)mm1*sizeof *ap);
+ bp = (float *) malloc((long)mm1*sizeof *bp);
/* reduce a and b into ap and bp eliminating variable t and constraint i */
@@ -559,7 +559,7 @@ if (nn == 1)
/* Constraints are inconsistent.
Set volume and derivative to zero. */
- printf("Inconsistent constraints found after reduction to n = 1 \n");
+/* printf("Inconsistent constraints found after reduction to n = 1 \n");*/
*dvdb = 0.;
return(0.);
}
@@ -623,8 +623,8 @@ v=0.;
/* allocate memory */
- ap = (float *) malloc((sizeof (float))*nm1*mm1);
- bp = (float *) malloc((sizeof (float))*mm1);
+ ap = (float *) malloc((long)nm1*(long)mm1*sizeof *ap);
+ bp = (float *) malloc((long)mm1*sizeof *bp);
/* reduce a and b into ap and bp eliminating variable t and constraint i */
@@ -796,7 +796,7 @@ if (nn == 1)
/* Constraint is inconsistent.
Set derivative to zero. */
- printf("Inconsistent constraints found after reduction to n = 1 \n");
+/* printf("Inconsistent constraints found after reduction to n = 1 \n");*/
*code = 1;
return(0.);
}
@@ -900,7 +900,7 @@ for (i=0;i<mm;i++)
if (special == 0)
{
- ttemp = (float *) malloc((sizeof (float))*mm1);
+ ttemp = (float *) malloc((long)mm1*sizeof *ttemp);
for (j=0;j<mm1;j++) *(ttemp+j) = - *(a+j+1+tmm)/pivot;
kval = 1;
}
@@ -916,7 +916,7 @@ for (i=0;i<mm;i++)
if(special == 0)
{
- ttemp = (float *) malloc((sizeof (float))*mm1);
+ ttemp = (float *) malloc((long)mm1*sizeof *ttemp);
for (j=0;j<i;j++) *(ttemp+j) = *(temp+j)
-(*(temp+i) * *(a+j+tmm)/pivot);
for (j=i;j<mm1;j++) *(ttemp+j) = *(temp+j+1)
@@ -932,8 +932,8 @@ for (i=0;i<mm;i++)
/* allocate memory */
- ap = (float *) malloc((sizeof (float))*nm1*mm1);
- bp = (float *) malloc((sizeof (float))*mm1);
+ ap = (float *) malloc((long)nm1*(long)mm1*sizeof *ap);
+ bp = (float *) malloc((long)mm1*sizeof *bp);
/* reduce a and b into ap and bp eliminating variable t and constraint i */
@@ -1064,8 +1064,8 @@ void dvda_();
void dvda();
void dvda_ (a,b,m,n,mmax,nmax,idim,result,code)
- int *n, *m, *mmax, *nmax, *idim, *code;
- float *a, *b;
+int *n, *m, *mmax, *nmax, *idim, *code;
+float *a, *b;
float *result;
{
dvda(a,b,m,n,mmax,nmax,idim,result,code);
@@ -1091,7 +1091,6 @@ tdim = *idim - 1;
free(temp);
}
-
/*--------------------------------------------------------------
ROUTINE: cvolumef
@@ -1167,7 +1166,7 @@ if (nn == 1)
}
else if ( *(b+l) < 0.) /* Constraints are inconsistent */
{
- printf("Inconsistent constraints found after reduction to n = 1 \n");
+/* printf("Inconsistent constraints found after reduction to n = 1 \n");*/
return(0.);
}
}
@@ -1219,8 +1218,8 @@ for (i=0;i<mm;i++)
/* allocate memory */
- ap = (float *) malloc((sizeof (float))*nm1*mm1);
- bp = (float *) malloc((sizeof (float))*mm1);
+ ap = (float *) malloc((long)nm1*(long)mm1*sizeof *ap);
+ bp = (float *) malloc((long)mm1*sizeof *bp);
/* reduce a and b into ap and bp eliminating variable t and constraint i */
@@ -1384,7 +1383,7 @@ if (nn == 1)
/* Constraint is inconsistent.
Set derivative to zero. */
- printf("Inconsistent constraints found after reduction to n = 1 \n");
+ /*printf("Inconsistent constraints found after reduction to n = 1 \n");*/
*code = 1;
return(0.);
}
@@ -1488,7 +1487,7 @@ for (i=0;i<mm;i++)
if (special == 0)
{
- ttemp = (float *) malloc((sizeof (float))*mm1);
+ ttemp = (float *) malloc((long)mm1*sizeof *ttemp);
for (j=0;j<mm1;j++) *(ttemp+j) = - *(a+j+1+tmm)/pivot;
kval = 1;
}
@@ -1504,7 +1503,7 @@ for (i=0;i<mm;i++)
if(special == 0)
{
- ttemp = (float *) malloc((sizeof (float))*mm1);
+ ttemp = (float *) malloc((long)mm1*sizeof *ttemp);
for (j=0;j<i;j++) *(ttemp+j) = *(temp+j)
-(*(temp+i) * *(a+j+tmm)/pivot);
for (j=i;j<mm1;j++) *(ttemp+j) = *(temp+j+1)
@@ -1520,8 +1519,8 @@ for (i=0;i<mm;i++)
/* allocate memory */
- ap = (float *) malloc((sizeof (float))*nm1*mm1);
- bp = (float *) malloc((sizeof (float))*mm1);
+ ap = (float *) malloc((long)nm1*(long)mm1*sizeof *ap);
+ bp = (float *) malloc((long)mm1*sizeof *bp);
/* reduce a and b into ap and bp eliminating variable t and constraint i */
@@ -1787,7 +1786,7 @@ if (nn == 1)
/* Constraint is inconsistent.
Set derivative to zero. */
- printf("Inconsistent constraints found after reduction to n = 1 \n");
+ /*printf("Inconsistent constraints found after reduction to n = 1 \n");*/
*code = 1;
return(0.);
}
@@ -1907,7 +1906,7 @@ for (i=0;i<mm;i++)
if (special == 0)
{
- ttemp = (float *) malloc((sizeof (float))*mm);
+ ttemp = (float *) malloc((long)mm*sizeof *ttemp);
for (j=0;j<mm1;j++) *(ttemp+j) = - *(a+j+1+tmm)/pivot;
kval = 1;
printf("\n Generating temp n = %d m = %d i = %d \n",nn,mm,i);
@@ -1927,7 +1926,7 @@ for (i=0;i<mm;i++)
if(special == 0)
{
- ttemp = (float *) malloc((sizeof (float))*mm1);
+ ttemp = (float *) malloc((long)mm1*sizeof *ttemp);
for (j=0;j<i;j++) *(ttemp+j) = *(temp+j)
-(*(temp+i) * *(a+j+tmm)/pivot);
for (j=i;j<mm1;j++) *(ttemp+j) = *(temp+j+1)
@@ -1946,8 +1945,8 @@ for (i=0;i<mm;i++)
/* allocate memory */
- ap = (float *) malloc((sizeof (float))*nm1*mm1);
- bp = (float *) malloc((sizeof (float))*mm1);
+ ap = (float *) malloc((long)nm1*(long)mm1*sizeof *ap);
+ bp = (float *) malloc((long)mm1*sizeof *bp);
/* reduce a and b into ap and bp eliminating variable t and constraint i */
diff --git a/OCTREE/OctreeBitPlus.f90 b/OCTREE/OctreeBitPlus.f90
index 8a11dba2..c37a95e3 100644
--- a/OCTREE/OctreeBitPlus.f90
+++ b/OCTREE/OctreeBitPlus.f90
@@ -112,6 +112,45 @@ z=dfloat(iz)/powermax
return
end
+!========================================!
+!=====[OCTREE_CREATE_FROM_PARTICLE]=====!
+!========================================!
+
+subroutine octree_create_from_particle (octree,noctree,x,y,z,np,level)
+
+! updates the octree by creating a leaf at point (x,y,z)
+! of level level
+! if the leaf (or a cube of smaller level) exists, the routine has no effect on the octree
+! note that x,y,z must belong to [0,1[
+
+implicit none
+
+integer noctree,octree(noctree),np
+double precision x,y,z
+integer level
+double precision xp,yp,zp
+
+integer levelin,ix,iy,iz,levelmax,loc,ip
+
+levelmax=octree(1)
+
+xp=x
+yp=y
+zp=z
+levelin=0
+loc=4
+ if (xp.eq.1.d0) xp=1.d0-1.d-20
+ if (yp.eq.1.d0) yp=1.d0-1.d-20
+ if (zp.eq.1.d0) zp=1.d0-1.d-20
+ if (xp.eq.0.d0) xp=1.d-20
+ if (yp.eq.0.d0) yp=1.d-20
+ if (zp.eq.0.d0) zp=1.d-20
+ if (xp*(xp-1.d0).lt.0.d0 .and. yp*(yp-1.d0).lt.0.d0 .and. zp*(zp-1.d0).lt.0.d0) &
+ call update (octree,noctree,xp,yp,zp,level,levelin,loc)
+
+return
+end
+
!========================================!
!=====[OCTREE_CREATE_FROM_PARTICLES]=====!
!========================================!
diff --git a/create_surfaces.f90 b/create_surfaces.f90
index 47c05c22..7be73cb4 100644
--- a/create_surfaces.f90
+++ b/create_surfaces.f90
@@ -46,7 +46,7 @@ integer debug
integer :: levelt, surface_type, nproc,iproc,ierr,i,ij,j,ie,i1,i2,i3,k
integer :: seed, indix,err
integer :: ntmax,nhmax,npmax,nmax,nmode,nvmax,nnpnmax
-integer :: nh,nohalt_hull,loc,nnn,nnlist,ntrilist
+integer :: nh,nohalt_hull,loc
integer :: dmode,nt,nbmax,it,nhull,icircles,irect
integer,dimension(:), allocatable :: hulltriangles
integer,dimension(:), allocatable :: vis_tlist,vis_elist,add_tlist,nb
@@ -64,6 +64,7 @@ logical clockwise
character(len=4) :: cs
character ch
character*72 :: shift
+integer nnn(1),nnlist(1),ntrilist(1)
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
@@ -126,7 +127,6 @@ else
!----------------
seed=1234567
-
sp01=surface%sp01
sp02=surface%sp02
sp03=surface%sp03
@@ -401,8 +401,10 @@ else
points(2,:)=surface%y
dmode=-2
nmode=0
+ loc=1
+ nohalt_hull=0
clockwise=.true.
- eps=0.d0
+ eps=1.d-8
call nn2d_setup (surface%nsurface,ntmax,nhmax,npmax,nnpnmax,nmax, &
points,dmode,nmode,clockwise,field,nt,vertices, &
centres,neighbour,nh,hulltriangles,nohalt_hull, &
diff --git a/improve_osolve.f90 b/improve_osolve.f90
index ca598e68..20f83c6f 100644
--- a/improve_osolve.f90
+++ b/improve_osolve.f90
@@ -107,16 +107,16 @@ shift=' '
call show_time (total,step,inc,1,'compute ref. criterion on ov$')
-allocate (x(mpe),y(mpe),z(mpe),stat=threadinfo%err)
-call heap (threadinfo,'x/y/z','improve_osolve',size(x)+size(y)+size(z),'dp',+1)
-allocate (vx(mpe),vy(mpe),vz(mpe),stat=threadinfo%err)
-call heap (threadinfo,'vx/y/z','improve_osolve',size(vx)+size(vy)+size(vz),'dp',+1)
-allocate (dhdx(mpe),dhdy(mpe),dhdz(mpe),stat=threadinfo%err)
-call heap (threadinfo,'dhdx/y/z','improve_osolve',size(dhdx)+size(dhdy)+size(dhdz),'dp',+1)
-allocate (crit(ov%nleaves),stat=threadinfo%err)
-call heap (threadinfo,'crit','improve_osolve',size(crit),'dp',+1)
-allocate (crittemp(ov%nleaves),stat=threadinfo%err)
-call heap (threadinfo,'crittemp','improve_osolve',size(crittemp),'dp',+1)
+allocate (x(mpe),y(mpe),z(mpe),stat=threadinfo%err)
+ call heap (threadinfo,'x/y/z','improve_osolve',size(x)+size(y)+size(z),'dp',+1)
+allocate (vx(mpe),vy(mpe),vz(mpe),stat=threadinfo%err)
+ call heap (threadinfo,'vx/y/z','improve_osolve',size(vx)+size(vy)+size(vz),'dp',+1)
+allocate (dhdx(mpe),dhdy(mpe),dhdz(mpe),stat=threadinfo%err)
+ call heap (threadinfo,'dhdx/y/z','improve_osolve',size(dhdx)+size(dhdy)+size(dhdz),'dp',+1)
+allocate (crit(ov%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'crit','improve_osolve',size(crit),'dp',+1)
+allocate (crittemp(ov%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'crittemp','improve_osolve',size(crittemp),'dp',+1)
crit=0.d0
crittemp=0.d0
@@ -208,7 +208,7 @@ if (critmax.gt.tiny(critmax)) then
! yyy=0.5d0*(ov%y(ov%icon(1,i))+ov%y(ov%icon(8,i)))
! zzz=0.5d0*(ov%z(ov%icon(1,i))+ov%z(ov%icon(8,i)))
! if (ov%whole_leaf_in_fluid(i)) &
-! call octree_create_from_particles(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,refine_level)
+! call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,refine_level)
! end if
!else
if (crit(i).gt. critmax*params%refine_ratio ) then
@@ -216,7 +216,7 @@ if (critmax.gt.tiny(critmax)) then
yyy=0.5d0*(ov%y(ov%icon(1,i))+ov%y(ov%icon(8,i)))
zzz=0.5d0*(ov%z(ov%icon(1,i))+ov%z(ov%icon(8,i)))
if (ov%whole_leaf_in_fluid(i)) &
- call octree_create_from_particles(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,refine_level)
+ call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,refine_level)
end if
!end if
@@ -258,7 +258,8 @@ if (params%nboxes.gt.0) then
if (.not.((xxx>=x0).and.(xxx<=x1) .and. &
(yyy>=y0).and.(yyy<=y1) .and. &
(zzz>=z0).and.(zzz<=z1))) call stop_run ('pb with box cloud$')
- call octree_create_from_particles(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level_box)
+if (iproc.eq.0) print*,i,j,k
+ call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level_box)
enddo
enddo
enddo
@@ -300,7 +301,7 @@ if (params%ref_on_faces) then
do j=1,nv
yyy=l+(real(i)-.5d0)*2.d0**(-(level+1))
zzz=b+(real(j)-.5d0)*2.d0**(-(level+1))
- call octree_create_from_particles(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level)
+ call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level)
enddo
enddo
@@ -310,7 +311,7 @@ if (params%ref_on_faces) then
do j=1,nv
yyy=l+(real(i)-.5d0)*2.d0**(-(level+1))
zzz=b+(real(j)-.5d0)*2.d0**(-(level+1))
- call octree_create_from_particles(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level)
+ call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level)
enddo
enddo
@@ -320,7 +321,7 @@ if (params%ref_on_faces) then
do j=1,nv
xxx=l+(real(i)-.5d0)*2.d0**(-(level+1))
zzz=b+(real(j)-.5d0)*2.d0**(-(level+1))
- call octree_create_from_particles(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level)
+ call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level)
enddo
enddo
@@ -330,7 +331,7 @@ if (params%ref_on_faces) then
do j=1,nv
xxx=l+(real(i)-.5d0)*2.d0**(-(level+1))
zzz=b+(real(j)-.5d0)*2.d0**(-(level+1))
- call octree_create_from_particles(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level)
+ call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level)
enddo
enddo
@@ -340,7 +341,7 @@ if (params%ref_on_faces) then
do j=1,nv
xxx=l+(real(i)-.5d0)*2.d0**(-(level+1))
yyy=b+(real(j)-.5d0)*2.d0**(-(level+1))
- call octree_create_from_particles(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level)
+ call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level)
enddo
enddo
case (6) ! z=1 face
@@ -349,7 +350,7 @@ if (params%ref_on_faces) then
do j=1,nv
xxx=l+(real(i)-.5d0)*2.d0**(-(level+1))
yyy=b+(real(j)-.5d0)*2.d0**(-(level+1))
- call octree_create_from_particles(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level)
+ call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level)
enddo
enddo
end select
diff --git a/wsmp_setup.f90 b/wsmp_setup.f90
index fca4e9af..65f3b95a 100644
--- a/wsmp_setup.f90
+++ b/wsmp_setup.f90
@@ -210,7 +210,7 @@ if (params%visualise_matrix) then
jcn(kk)=i
enddo
enddo
- call visualise_matrix (nz,irn,jcn,n,istep,iter,ndof)
+! call visualise_matrix (nz,irn,jcn,n,istep,iter,ndof)
call heap (threadinfo,'irn','wsmp_setup',size(irn),'int',-1) ; deallocate (irn)
call heap (threadinfo,'jcn','wsmp_setup',size(jcn),'int',-1) ; deallocate (jcn)
end if
--
GitLab