Commit 0d998314 authored by Adam Powell's avatar Adam Powell

Imported Upstream version 5.5.0.svn.4876.dfsg

parent 2d74d37e
This diff is collapsed.
#==============================================================================
#
# ElmerGUI: configuration file
#
#==============================================================================
#------------------------------------------------------------------------------
# Optional components (undefine or comment out to exclude from compilation):
#------------------------------------------------------------------------------
DEFINES += EG_QWT # Use QWT for convergence monitor?
DEFINES += EG_VTK # Use VTK for postprocessing?
DEFINES += EG_MATC # Use MATC for internal operations in postprocessing?
DEFINES += EG_OCC # Use OpenCASCADE 6.3 for importing CAD files? Needs VTK.
DEFINES -= EG_PYTHONQT # Use PythonQt for scripting in post processor?
#------------------------------------------------------------------------------
# 64 bit system?
#------------------------------------------------------------------------------
BITS = 32
#------------------------------------------------------------------------------
# Installation directory:
#------------------------------------------------------------------------------
ELMERGUI_HOME = $$(ELMERGUI_HOME)
isEmpty(ELMERGUI_HOME) {
ELMER_HOME = $$(ELMER_HOME)
isEmpty(ELMER_HOME) {
unix: ELMER_HOME = /usr/local
win32: ELMER_HOME = c:\Elmer5.4
macx: ELMER_HOME = /usr/local
}
ELMERGUI_HOME = $${ELMER_HOME}/bin
}
#------------------------------------------------------------------------------
# Python library:
#------------------------------------------------------------------------------
unix {
PY_INCLUDEPATH = /usr/include/python2.5
PY_LIBPATH = /usr/lib
PY_LIBS = -lpython2.5
}
win32 {
PY_INCLUDEPATH = c:\PYTHON\Python-2.6.1\Include
PY_LIBPATH = c:\PYTHON\Python-2.6.1\PCbuild
PY_LIBS = -lpython26
}
macx {
}
#------------------------------------------------------------------------------
# QWT library:
#------------------------------------------------------------------------------
unix {
QWT_INCLUDEPATH = /usr/include/qwt-qt4
QWT_LIBPATH = /usr/lib
QWT_LIBS = -lqwt-qt4
}
win32 {
QWT_INCLUDEPATH = c:\Source\Qwt\include
QWT_LIBPATH = c:\Source\Qwt\lib
QWT_LIBS = -lqwt5
}
macx {
QWT_INCLUDEPATH = /usr/local/qwt-5.0.2/include
QWT_LIBPATH = /usr/local/qwt-5.0.2/lib
QWT_LIBS = -lqwt5
}
#------------------------------------------------------------------------------
# VTK library:
#------------------------------------------------------------------------------
unix {
VTK_INCLUDEPATH = /usr/include/vtk-5.2
VTK_LIBPATH = /usr/lib
VTK_LIBS = -lvtkHybrid \
-lvtkWidgets \
-lQVTK
}
win32 {
VTK_INCLUDEPATH = c:\Source\VTK\include\vtk-5.4
VTK_LIBPATH = c:\Source\VTK\lib\vtk-5.4
VTK_LIBS = -lQVTK \
-lvtkCommon \
-lvtkDICOMParser \
-lvtkFiltering \
-lvtkGenericFiltering \
-lvtkGraphics \
-lvtkHybrid \
-lvtkIO \
-lvtkImaging \
-lvtkInfovis \
-lvtkNetCDF \
-lvtkRendering \
-lvtkViews \
-lvtkVolumeRendering \
-lvtkWidgets \
-lvtkexoIIc \
-lvtkexpat \
-lvtkfreetype \
-lvtkftgl \
-lvtkjpeg \
-lvtklibxml2 \
-lvtkmetaio \
-lvtkpng \
-lvtksys \
-lvtktiff \
-lvtkverdict \
-lvtkzlib \
-ladvapi32
}
macx {
VTK_INCLUDEPATH = /usr/local/include/vtk-5.0
VTK_LIBPATH = /usr/lib
VTK_LIBS = -lvtkHybrid \
-lvtkWidgets \
-lQVTK
}
#------------------------------------------------------------------------------
# OpenCASCADE library:
#------------------------------------------------------------------------------
unix {
OCC_INCLUDEPATH = /usr/include/opencascade
OCC_LIBPATH = /usr/lib
OCC_LIBS = -lTKBRep -lTKSTL -lTKSTEP -lTKIGES
}
win32 {
OCC_INCLUDEPATH = $(CASROOT)/inc
OCC_LIBPATH = $(CASROOT)/win32/lib
OCC_LIBS = $(CASROOT)/win32/lib/TKBRep.lib \
$(CASROOT)/win32/lib/TKernel.lib \
$(CASROOT)/win32/lib/TKG2d.lib \
$(CASROOT)/win32/lib/TKG3d.lib \
$(CASROOT)/win32/lib/TKGeomAlgo.lib \
$(CASROOT)/win32/lib/TKGeomBase.lib \
$(CASROOT)/win32/lib/TKMath.lib \
$(CASROOT)/win32/lib/TKMesh.lib \
$(CASROOT)/win32/lib/TKShHealing.lib \
$(CASROOT)/win32/lib/TKSTEP.lib \
$(CASROOT)/win32/lib/TKSTEP209.lib \
$(CASROOT)/win32/lib/TKSTEPAttr.lib \
$(CASROOT)/win32/lib/TKSTEPBase.lib \
$(CASROOT)/win32/lib/TKIGES.lib \
$(CASROOT)/win32/lib/TKTopAlgo.lib \
$(CASROOT)/win32/lib/TKXSBase.lib
}
macx {
OCC_INCLUDEPATH =
OCC_LIBPATH =
OCC_LIBS =
}
##### ElmerGrid input file for structured grid generation ######
Version = 210903
Coordinate System = Cartesian 2D
Subcell Divisions in 2D = 1 1
Subcell Sizes 1 = 0.06
Subcell Sizes 2 = 0.01
Material Structure in 2D
1
End
Materials Interval = 1 1
Boundary Definitions
# type out int
1 -3 1 1
2 -4 1 1
3 -1 1 1
4 -2 1 1
End
Numbering = Horizontal
Element Degree = 1
Element Innernodes = False
Surface Elements = 3000
Element Ratios 1 = 1 1
Element Ratios 2 = 1 1
Element Densities 1 = 1 1
Element Densities 2 = 1 1
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -1029,7 +1029,10 @@ END FUNCTION CRS_RowSum
TYPE(Matrix_t) :: A
!------------------------------------------------------------------------------
! INTEGER, POINTER, CONTIGUOUS :: Cols(:),Rows(:)
INTEGER, POINTER :: Cols(:),Rows(:)
! REAL(KIND=dp), POINTER, CONTIGUOUS :: Values(:)
REAL(KIND=dp), POINTER :: Values(:)
INTEGER :: i,j,n
......@@ -1041,6 +1044,11 @@ END FUNCTION CRS_RowSum
Cols => A % Cols
Values => A % Values
IF ( A % MatvecSubr /= 0 ) THEN
CALL MatVecSubr(A % MatVecSubr,A % SpMV, n,Rows,Cols,Values,u,v,0)
RETURN
END IF
!$omp parallel do private(j,rsum)
DO i=1,n
rsum = 0.0d0
......@@ -2669,7 +2677,9 @@ END FUNCTION CRS_RowSum
REAL(KIND=dp) :: u(HUTI_NDIM),v(HUTI_NDIM)
!------------------------------------------------------------------------------
! INTEGER, POINTER, CONTIGUOUS :: Cols(:),Rows(:)
INTEGER, POINTER :: Cols(:),Rows(:)
! REAL(KIND=dp), POINTER, CONTIGUOUS :: Values(:)
REAL(KIND=dp), POINTER :: Values(:)
INTEGER :: i,j,n
......@@ -2681,6 +2691,12 @@ END FUNCTION CRS_RowSum
Cols => GlobalMatrix % Cols
Values => GlobalMatrix % Values
IF ( GlobalMatrix % MatVecSubr /= 0 ) THEN
CALL MatVecSubr(GlobalMatrix % MatVecSubr, &
GlobalMatrix % SpMV, n,Rows,Cols,Values,u,v,0)
print*,globalmatrix % spMV
RETURN
END IF
!***********************************************************************
! The following #ifdefs seem really necessery, if speed is an issue:
! SGI compiler optimizer wants to know the sizes of the arrays very
......
......@@ -1089,10 +1089,7 @@ CONTAINS
GB = ListGetLogical( Solver % Values, 'Bubbles in Global System', Found )
IF (.NOT.Found) GB = .TRUE.
IF ( ASSOCIATED(Element % BoundaryInfo) .AND. .NOT. &
( ASSOCIATED(Element % EdgeIndexes) .OR. &
ASSOCIATED(Element % FaceIndexes) ) ) THEN
IF ( ASSOCIATED(Element % BoundaryInfo) ) THEN
Parent => Element % BoundaryInfo % Left
IF (.NOT.ASSOCIATED(Parent) ) &
Parent => Element % BoundaryInfo % Right
......@@ -1150,7 +1147,7 @@ CONTAINS
END IF
END SELECT
ELSE IF ( GB ) THEN
IF ( ASSOCIATED( Element % BubbleIndexes ) ) THEN
IF ( ASSOCIATED(Element % BubbleIndexes) ) THEN
DO i=1,Element % BDOFs
NB = NB + 1
Indexes(NB) = FaceDOFs*Solver % Mesh % NumberOfFaces + &
......
......@@ -268,7 +268,7 @@ SUBROUTINE DistanceSolver1( Model,Solver,dt,TransientSimulation )
!------------------------------------------------------------------------------
! Local variables
!------------------------------------------------------------------------------
LOGICAL :: AllocationsDone = .FALSE., Found
LOGICAL :: AllocationsDone = .FALSE., Found, Found1
TYPE(Element_t),POINTER :: Element
TYPE(ValueList_t), POINTER :: SolverParams, BC
......@@ -278,7 +278,7 @@ SUBROUTINE DistanceSolver1( Model,Solver,dt,TransientSimulation )
TYPE(Mesh_t), POINTER :: Mesh
TYPE(ValueList_t), POINTER :: BodyForce
REAL(KIND=dp), ALLOCATABLE :: STIFF(:,:), FORCE(:), buf(:) ,xp(:),yp(:),zp(:)
REAL(KIND=dp), ALLOCATABLE :: STIFF(:,:), FORCE(:), buf(:) ,xp(:),yp(:),zp(:),bdist(:),bd(:)
REAL(KIND=dp), POINTER :: distance(:)
INTEGER, POINTER :: gPerm(:), ibuf(:), aperm(:),bperm(:),cperm(:)
......@@ -292,7 +292,7 @@ REAL(kind=dp) :: cputime, st
Mesh => GetMesh()
n = Mesh % NumberOfNodes
ALLOCATE( aperm(n), bperm(n) ); aperm = 0; bperm = 0
ALLOCATE( aperm(n), bperm(n), bdist(n) ); aperm = 0; bperm = 0; bdist=0
nb = 0
DO t=1,Mesh % NumberOfBoundaryElements
......@@ -300,8 +300,8 @@ REAL(kind=dp) :: cputime, st
IF ( .NOT. ActiveBoundaryElement() ) CYCLE
BC => GetBC()
IF ( ListCheckPresent(BC,ComponentName(Solver % Variable)).OR. &
GetLogical( BC, 'Noslip Wall BC', Found ) ) THEN
dist = GetCReal(BC, Solver % Variable % Name, Found )
IF ( Found .OR. GetLogical( BC, 'Noslip Wall BC', Found1 ) ) THEN
nd = GetElementNOFNodes()
DO i=1,nd
j = Element % NodeIndexes(i)
......@@ -309,6 +309,7 @@ REAL(kind=dp) :: cputime, st
nb = nb + 1
aperm(nb) = j
bperm(j) = nb
bdist(j) = dist
END IF
END DO
END IF
......@@ -334,7 +335,7 @@ REAL(kind=dp) :: cputime, st
END IF
END DO
ALLOCATE( xp(nnb), yp(nnb), zp(nnb), buf(nnb) )
ALLOCATE( xp(nnb), yp(nnb), zp(nnb), bd(nnb), buf(nnb) )
buf=0._dp
buf(cperm(gperm(aperm(1:nb)))) = Mesh % Nodes % x(aperm(1:nb))
......@@ -351,14 +352,20 @@ REAL(kind=dp) :: cputime, st
CALL MPI_ALLREDUCE(buf, zp, nnb, MPI_DOUBLE_PRECISION, MPI_SUM, comm, i)
zp(1:nnb) = zp(1:nnb)/ibuf(1:nnb)
buf=0._dp
buf(cperm(gperm(aperm(1:nb)))) = bdist(aperm(1:nb))
CALL MPI_ALLREDUCE(buf, bd, nnb, MPI_DOUBLE_PRECISION, MPI_SUM, comm, i)
bd(1:nnb) = bd(1:nnb)/ibuf(1:nnb)
nb = nnb
DEALLOCATE(ibuf,buf,cperm)
ELSE
ALLOCATE( xp(nb),yp(nb),zp(nb) )
ALLOCATE( xp(nb),yp(nb),zp(nb),bd(nb) )
DO i=1,nb
xp(i) = Solver % Mesh % Nodes % x(aperm(i))
yp(i) = Solver % Mesh % Nodes % y(aperm(i))
zp(i) = Solver % Mesh % Nodes % z(aperm(i))
bd(i) = bdist(aperm(i))
END DO
END IF
......@@ -372,10 +379,7 @@ REAL(kind=dp) :: cputime, st
! if ( parenv % mype==0) print*,'algo 1:', cputime()-st
DEALLOCATE(xp,yp,zp,aperm,bperm)
Solver % Variable % Norm = SQRT(SUM(distance))
distance = SQRT(distance)
Solver % Variable % Norm = SQRT(SUM(distance**2))
CALL Info('DistanceSolver1','All done')
......@@ -407,24 +411,26 @@ CONTAINS
SUBROUTINE distcomp
REAL(KIND=dp), ALLOCATABLE :: xxp(:),yyp(:),zzp(:),d(:),dd(:)
REAL(KIND=dp) :: xl,yl,zl,dl,c(3)
INTEGER :: i,j,k,l,nnb,cnt
REAL(KIND=dp), ALLOCATABLE :: xxp(:),yyp(:),zzp(:),d(:),dd(:), bbd(:)
REAL(KIND=dp) :: xl,yl,zl,dl,c(3),ldist
INTEGER :: i,j,k,l,nnb,cnt,minl
INTEGER, ALLOCATABLE :: dperm(:), near(:)
nnb = n+nb
ALLOCATE( xxp(nnb), yyp(nnb), zzp(nnb), d(nnb), dperm(nnb), &
near(nnb), dd(nb) )
near(nnb), dd(nb), bbd(nnb) )
CALL RANDOM_NUMBER(c)
xxp(1:n) = Mesh % Nodes % x
yyp(1:n) = Mesh % Nodes % y
zzp(1:n) = Mesh % Nodes % z
bbd(1:n) = bdist(1:n)
xxp(n+1:nnb) = xp
yyp(n+1:nnb) = yp
zzp(n+1:nnb) = zp
bbd(n+1:nnb) = bd
xl = MAXVAL(xxp)-MINVAL(xxp)
yl = MAXVAL(yyp)-MINVAL(yyp)
......@@ -447,6 +453,7 @@ CONTAINS
yp(j)=yyp(dperm(i))
zp(j)=zzp(dperm(i))
dd(j)=d(i)
bd(j)=bbd(dperm(i))
END IF
END DO
......@@ -458,7 +465,7 @@ CONTAINS
IF ( k <= 0 ) CYCLE
IF ( bperm(j)/=0 ) THEN
distance(k) = 0._dp
distance(k) = bdist(j)
ELSE
dl = d(i)
xl = xxp(j)
......@@ -467,17 +474,23 @@ CONTAINS
dist = HUGE(1._dp)
DO l=near(i)+1,nb
IF( (dd(l)-dl)**2>dist ) EXIT
dist = MIN(dist,(xp(l)-xl)**2+(yp(l)-yl)**2+(zp(l)-zl)**2)
ldist = (xp(l)-xl)**2+(yp(l)-yl)**2+(zp(l)-zl)**2
IF ( dist>ldist) THEN
minl = l; dist=ldist
END IF
END DO
DO l=near(i),1,-1
IF( (dd(l)-dl)**2>dist ) EXIT
dist = MIN(dist,(xp(l)-xl)**2+(yp(l)-yl)**2+(zp(l)-zl)**2)
ldist = (xp(l)-xl)**2+(yp(l)-yl)**2+(zp(l)-zl)**2
IF ( dist>ldist) THEN
minl = l; dist=ldist
END IF
END DO
distance(k) = dist
distance(k) = SQRT(dist)+bd(minl)
END IF
END DO
DEALLOCATE( xxp, yyp, zzp, d, dperm, near, dd )
DEALLOCATE( xxp, yyp, zzp, d, dperm, near, dd, bbd )
END SUBROUTINE distcomp
!------------------------------------------------------------------------------
END SUBROUTINE DistanceSolver1
......
......@@ -909,7 +909,7 @@
INTEGER :: interval, timestep, i, j, k, n
REAL(KIND=dp) :: dt, ddt, dtfunc
INTEGER :: timeleft,cum_timestep
INTEGER, SAVE :: stepcount=0
INTEGER, SAVE :: stepcount=0, RealTimestep
LOGICAL :: ExecThis,SteadyStateReached=.FALSE.
REAL(KIND=dp) :: CumTime, MaxErr, AdaptiveLimit, &
......@@ -947,6 +947,7 @@
!------------------------------------------------------------------------------
! go trough number of timesteps within an interval
!------------------------------------------------------------------------------
RealTimestep = 1
DO timestep = 1,Timesteps(interval)
dtfunc = ListGetConstReal( CurrentModel % Simulation, &
......@@ -1079,7 +1080,7 @@
sTime(1) = s + CumTime + ddt
sSize(1) = ddt
CALL SolveEquations( CurrentModel, ddt, Transient, &
CoupledMinIter, CoupledMaxIter, SteadyStateReached )
CoupledMinIter, CoupledMaxIter, SteadyStateReached, RealTimestep )
MaxErr = ListGetConstReal( CurrentModel % Simulation, &
......@@ -1102,10 +1103,10 @@
sStep(1) = ddt / 2
sTime(1) = s + CumTime + ddt/2
CALL SolveEquations( CurrentModel, ddt/2, Transient, &
CoupledMinIter, CoupledMaxIter, SteadyStateReached )
CoupledMinIter, CoupledMaxIter, SteadyStateReached, RealTimestep )
sTime(1) = s + CumTime + ddt
CALL SolveEquations( CurrentModel, ddt/2, Transient, &
CoupledMinIter, CoupledMaxIter, SteadyStateReached )
CoupledMinIter, CoupledMaxIter, SteadyStateReached, RealTimestep )
MaxErr = ABS( MaxErr - ListGetConstReal( CurrentModel % Simulation, &
'Adaptive Error Measure', GotIt ) )
......@@ -1124,6 +1125,7 @@
IF ( MaxErr < AdaptiveLimit .OR. ddt <= AdaptiveMinTimestep ) THEN
CumTime = CumTime + ddt
RealTimestep = RealTimestep+1
IF ( SmallestCount >= AdaptiveKeepSmallest .OR. StepControl > 0 ) THEN
ddt = MIN( 2*ddt, AdaptiveMaxTimeStep )
StepControl = 1
......@@ -1154,10 +1156,11 @@
sSize(1) = dt
sTime(1) = s + dt
DEALLOCATE( xx, xxnrm, yynrm )
DEALLOCATE( xx, xxnrm, yynrm, prevxx )
ELSE ! Adaptive timestepping
CALL SolveEquations( CurrentModel, dt, Transient, &
CoupledMinIter, CoupledMaxIter, SteadyStateReached )
CoupledMinIter, CoupledMaxIter, SteadyStateReached, RealTimestep )
RealTimestep = RealTimestep+1
END IF
!------------------------------------------------------------------------------
! Save results to disk, if requested
......
......@@ -236,14 +236,11 @@ SUBROUTINE FluxSolver( Model,Solver,dt,Transient )
RETURN
END IF
PRINT *,'Dim',dim,dofs,CalculateFlux,CalculateFluxMag,CalculateGrad,CalculateGradMag
!-------------------------------------------------------------------------------
! If only one component is used use the scalar equation, otherwise use an
! auxiliary variable to store all the dimensions
!-------------------------------------------------------------------------------
VarName = GetString(SolverParams,'Flux Variable',GotIt )
IF(.NOT. GotIt) VarName = GetString(SolverParams,'Target Variable',GotIt )
IF(.NOT. gotIt) THEN
......
......@@ -1222,7 +1222,7 @@
IF(.NOT. GotIt) PowerSensitivity = 4.0_dp
PowerScaling = PowerScaling * (1 + PowerSensitivity * PowerRelax * (MeltPoint/yave - 1.0d0) )
print *,'MeltPoint',MeltPoint,yave
PRINT *,'MeltPoint',MeltPoint,yave
IF( ListGetLogical( Solver % Values,'Smart Heater Transient Speedup',GotIt ) ) THEN
Temperature = Temperature * (1 + PowerRelax * (MeltPoint/yave - 1.0d0) )
......@@ -1552,7 +1552,7 @@ CONTAINS
!------------------------------------------------------------------------------
SUBROUTINE EffectiveHeatCapacity
LOGICAL :: Found
LOGICAL :: Found, Specific
!------------------------------------------------------------------------------
! See if temperature gradient indside the element is large enough
! to use the c_p = SQRT( (dH/dx)^2 / (dT/dx)^2 ), otherwise
......@@ -1607,6 +1607,8 @@ CONTAINS
END IF
PhaseSpatial = ( PhaseChangeModel == PHASE_SPATIAL_2 )
Specific = ListCheckPresent( Material,'Specific Enthalpy')
!------------------------------------------------------------------------------
SELECT CASE( PhaseChangeModel )
!------------------------------------------------------------------------------
......@@ -1614,21 +1616,45 @@ CONTAINS
HeatCapacity(1:n) = ListGetReal( Material, &
'Effective Heat Capacity', n,Element % NodeIndexes, Found )
IF ( .NOT. Found ) THEN
HeatCapacity(1:n) = ListGetDerivValue( Material, &
'Enthalpy', n,Element % NodeIndexes )
IF( Specific ) THEN
HeatCapacity(1:n) = ListGetDerivValue( Material, &
'Specific Enthalpy', n,Element % NodeIndexes )
HeatCapacity(1:n) = Density(1:n) * HeatCapacity(1:n)
ELSE
HeatCapacity(1:n) = ListGetDerivValue( Material, &
'Enthalpy', n,Element % NodeIndexes )
END IF
END IF
!------------------------------------------------------------------------------
CASE( PHASE_SPATIAL_2 )
Enthalpy(1:n) = ListGetReal(Material,'Enthalpy',n,Element % NodeIndexes)
IF( Specific ) THEN
Enthalpy(1:n) = ListGetReal(Material,'Specific Enthalpy',n,Element % NodeIndexes)
Enthalpy(1:n) = Density(1:n) * Enthalpy(1:n)
ELSE
Enthalpy(1:n) = ListGetReal(Material,'Enthalpy',n,Element % NodeIndexes)
END IF
!------------------------------------------------------------------------------
CASE( PHASE_TEMPORAL )
TSolution1 = Temperature(1:LocalNodes)
Work(1:n) = ListGetReal( Material,'Enthalpy',n,Element % NodeIndexes )
Temperature(1:LocalNodes) = TSolution
Work(1:n) = Work(1:n) - ListGetReal( Material,'Enthalpy', &
n,Element % NodeIndexes )
IF( Specific ) THEN
Work(1:n) = ListGetReal( Material,'Specific Enthalpy',n,Element % NodeIndexes )
Temperature(1:LocalNodes) = TSolution
Work(1:n) = Work(1:n) - ListGetReal( Material,'Specific Enthalpy', &
n,Element % NodeIndexes )
ELSE
Work(1:n) = ListGetReal( Material,'Enthalpy',n,Element % NodeIndexes )
Temperature(1:LocalNodes) = TSolution
Work(1:n) = Work(1:n) - ListGetReal( Material,'Enthalpy', &
n,Element % NodeIndexes )
END IF
Temperature(1:LocalNodes) = TSolution1
HeatCapacity(1:n) = Work(1:n) / HeatCapacity(1:n)
IF( Specific ) THEN
HeatCapacity(1:n) = Density(1:n) * HeatCapacity(1:n)
END IF
!------------------------------------------------------------------------------
END SELECT
!------------------------------------------------------------------------------
......
......@@ -408,58 +408,53 @@
IF ( ASSOCIATED( BC ) ) THEN
IF ( GetLogical( BC, 'Wall Law',gotIt ) ) THEN
! /*
! * NOTE: that the following is not really valid, the pointer to
! * the Material structure is from the remains of the last of the
! * bulk elements.
! */
Density(1:n) = GetReal( Material,'Density' )
Viscosity(1:n) = GetReal( Material,'Viscosity' )
SurfaceRoughness(1:n) = GetReal( BC, 'Surface Roughness', gotIt )
LayerThickness(1:n) = GetReal( BC, 'Boundary Layer Thickness' )
DO j=1,n
k = FlowPerm(NodeIndexes(j))
IF ( k > 0 ) THEN
SELECT CASE( NSDOFs )
CASE(3)
U(j) = FlowSolution( NSDOFs*k-2 )
V(j) = FlowSolution( NSDOFs*k-1 )
W(j) = 0.0D0
CASE(4)
U(j) = FlowSolution( NSDOFs*k-3 )
V(j) = FlowSolution( NSDOFs*k-2 )
W(j) = FlowSolution( NSDOFs*k-1 )
END SELECT
ELSE
U(j) = 0.0d0
V(j) = 0.0d0
W(j) = 0.0d0
END IF
END DO
DO j=1,n
CALL KEWall( Work(1), Work(2), SQRT(U(j)**2+V(j)**2+W(j)**2), &
LayerThickness(j), SurfaceRoughness(j), Viscosity(j), &
Density(j) )
k = DOFs*(KinPerm(NodeIndexes(j))-1)
ForceVector(k+1) = Work(1)
CALL ZeroRow( StiffMatrix,k+1 )
CALL SetMatrixElement( StiffMatrix,k+1,k+1,1.0d0 )
ForceVector(k+2) = Work(2)
CALL ZeroRow( StiffMatrix,k+2 )
CALL SetMatrixElement( StiffMatrix,k+2,k+2,1.0d0 )
END DO
END IF
IF ( GetLogical( BC, 'Epsilon Wall BC', gotIt ) .OR. &
GetLogical( BC, 'Noslip Wall BC', gotIt ) ) THEN
CALL EpsilonWall( Element, n, STIFF, FORCE )
END IF
Density(1:n) = GetParentMatProp( 'Density' )
Viscosity(1:n) = GetParentMatProp( 'Viscosity' )
SurfaceRoughness(1:n) = GetReal( BC, 'Surface Roughness', gotIt )
LayerThickness(1:n) = GetReal( BC, 'Boundary Layer Thickness' )
DO j=1,n
k = FlowPerm(NodeIndexes(j))
IF ( k > 0 ) THEN
SELECT CASE( NSDOFs )
CASE(3)
U(j) = FlowSolution( NSDOFs*k-2 )
V(j) = FlowSolution( NSDOFs*k-1 )
W(j) = 0.0D0
CASE(4)
U(j) = FlowSolution( NSDOFs*k-3 )
V(j) = FlowSolution( NSDOFs*k-2 )
W(j) = FlowSolution( NSDOFs*k-1 )
END SELECT
ELSE
U(j) = 0.0d0
V(j) = 0.0d0
W(j) = 0.0d0
END IF
END DO
DO j=1,n
CALL KEWall( Work(1), Work(2), Work(3), SQRT(U(j)**2+V(j)**2+W(j)**2), &
LayerThickness(j), SurfaceRoughness(j), Viscosity(j), &
Density(j) )
k = DOFs*(KinPerm(NodeIndexes(j))-1)
ForceVector(k+1) = Work(1)
CALL ZeroRow( StiffMatrix,k+1 )
CALL SetMatrixElement( StiffMatrix,k+1,k+1,1.0d0 )
ForceVector(k+2) = Work(2)
CALL ZeroRow( StiffMatrix,k+2 )
CALL SetMatrixElement( StiffMatrix,k+2,k+2,1.0d0 )
END DO
END IF
IF ( GetLogical( BC, 'Epsilon Wall BC', gotIt ) .OR. &
GetLogical( BC, 'Noslip Wall BC', gotIt ) ) THEN
CALL EpsilonWall( Element, n, STIFF, FORCE )
END IF
END IF
END DO
!------------------------------------------------------------------------------
......
......@@ -624,4 +624,24 @@ void STDCALLBULL FC_FUNC(execlocalassembly, EXECLOCALASSEMBLY )
/*--------------------------------------------------------------------------
INTERNAL: execute complete localmatrix call
-------------------------------------------------------------------------*/
static void DoMatVecSubr( void (STDCALLBULL *matvec)(),
void **SpMV,void *n,void *rows,void *cols,void *vals,void *u, void *v, void *reinit )
{
(*matvec)( SpMV,n,rows,cols,vals,u,v,reinit);
}
/*--------------------------------------------------------------------------
This routine will call complete local matrix add-on
-------------------------------------------------------------------------*/
void STDCALLBULL FC_FUNC(matvecsubr, MMATVECSUBR)
( f_ptr matvec, void **SpMV, void *n, void *rows, void *cols, void *vals, void *u, void *v,void *reinit )
{
DoMatVecSubr( (void (STDCALLBULL *)())*matvec,SpMV,n,rows,cols,vals,u,v,reinit);
}
This diff is collapsed.
......@@ -485,7 +485,7 @@ MODULE NavierStokes
Element, Nodes, n, n, u, v, w, muder0 )
ViscNewtonLin = NewtonLinearization .AND. muder0/=0._dp
IF ( ViscNewtonLin ) Strain = Grad+TRANSPOSE(Grad)
IF ( ViscNewtonLin ) Strain = (Grad+TRANSPOSE(Grad))/2
ELSE
DO i=1,6
......@@ -733,7 +733,7 @@ MODULE NavierStokes
IF ( ViscNewtonLin ) THEN
DO i=1,dim
muder = 2*muder0*SUM(Strain(i,:)*dBasisdx(q,:))
muder = muder0*SUM(Strain(i,:)*dBasisdx(q,:))
DO j=1,dim
Jac(i,i)=Jac(i,i) + s*muder*Grad(i,j)*dBasisdx(p,j)