Commit 2d74d37e authored by Adam Powell's avatar Adam Powell

Imported Upstream version 5.5.0.svn.4844.dfsg

parent e4b79721
......@@ -21,6 +21,9 @@ INSTALLS += target
edf.path = $${ELMERGUI_HOME}/edf
edf.files = edf/*
INSTALLS += edf
edf-extra.path = $${ELMERGUI_HOME}/edf-extra
edf-extra.files = edf-extra/*
INSTALLS += edf-extra
#------------------------------------------------------------------------------
# Compiler flags:
......
......@@ -96,6 +96,12 @@
<StatusTip> Set value to maximum solubility. </StatusTip>
<Whatis> Give value to maximum solubility of solute. </Whatis>
</Parameter>
<Parameter Widget="Edit" >
<Name> Ion Charge </Name>
<SifName> Concentration Ion Charge </SifName>
<StatusTip> Set value to charge of solute ions, e.g. 1, -2 etc. </StatusTip>
<Whatis> Give value of ion charge. </Whatis>
</Parameter>
</Material>
<InitialCondition>
......
......@@ -53,6 +53,12 @@ using namespace std;
#define MY_PI 3.14159265
#define ZSHIFT -5.0
// Get qreal regardless of whether it's float or double
static inline void glGetQrealv(GLenum e, GLfloat* data) { glGetFloatv(e,data); }
static inline void glGetQrealv(GLenum e, GLdouble* data) { glGetDoublev(e,data); }
static inline void glMultMatrixq( const GLdouble *m ) { glMultMatrixd(m); }
static inline void glMultMatrixq( const GLfloat *m ) { glMultMatrixf(m); }
list_t::list_t()
{
nature = PDE_UNKNOWN;
......@@ -730,7 +736,7 @@ void GLWidget::mouseMoveEvent(QMouseEvent *event)
double az = 0.0;
glLoadIdentity();
glTranslated(ax, ay, az);
glMultMatrixd(matrix);
glMultMatrixq(matrix);
updateGL();
}
......@@ -1011,7 +1017,7 @@ void GLWidget::mouseDoubleClickEvent(QMouseEvent *event)
//-----------------------------------------------------------------------------
void GLWidget::getMatrix()
{
glGetDoublev(GL_MODELVIEW_MATRIX, matrix);
glGetQrealv(GL_MODELVIEW_MATRIX, matrix);
helpers->invertMatrix(matrix, invmatrix);
}
......
......@@ -183,8 +183,8 @@ private:
GLuint makeLists();
GLdouble matrix[16];
GLdouble invmatrix[16];
qreal matrix[16];
qreal invmatrix[16];
void getMatrix();
QPoint lastPos;
......
......@@ -57,7 +57,7 @@ Helpers :: ~Helpers()
// Normalize
//====================================================================
void Helpers::normalize(double *a)
void Helpers::normalize(qreal *a)
{
double b;
......@@ -71,7 +71,7 @@ void Helpers::normalize(double *a)
// Length
//====================================================================
double Helpers::vlen(double *a)
qreal Helpers::vlen(qreal *a)
{
return sqrt(a[0]*a[0] + a[1]*a[1] + a[2]*a[2]);
}
......@@ -80,7 +80,7 @@ double Helpers::vlen(double *a)
// Cross product
//====================================================================
void Helpers::crossProduct(double *a, double *b, double *c)
void Helpers::crossProduct(qreal *a, qreal *b, qreal *c)
{
c[0] = a[1]*b[2] - a[2]*b[1];
c[1] = a[2]*b[0] - a[0]*b[2];
......@@ -90,7 +90,7 @@ void Helpers::crossProduct(double *a, double *b, double *c)
//====================================================================
// Invert 4x4 matrix (for visualiztion only)
//====================================================================
void Helpers::invertMatrix(const double *a, double *inva)
void Helpers::invertMatrix(const qreal *a, qreal *inva)
{
QMatrix4x4 matrix(a);
......
......@@ -42,6 +42,7 @@
#define HELPERS_H
#include "meshtype.h"
#include <QMatrix4x4>
class Helpers
{
......@@ -49,10 +50,10 @@ class Helpers
Helpers();
~Helpers();
void invertMatrix(const double *a, double *b);
void crossProduct(double *a, double *b, double *c);
double vlen(double *a);
void normalize(double *a);
void invertMatrix(const qreal *a, qreal *b);
void crossProduct(qreal *a, qreal *b, qreal *c);
qreal vlen(qreal *a);
void normalize(qreal *a);
private:
......
......@@ -1062,7 +1062,7 @@ void Meshutils::findSurfaceElementEdges(mesh_t *mesh)
//-----------------------------------------------------------------------------
void Meshutils::findSharpPoints(mesh_t *mesh, double limit)
{
double t0[3], t1[3];
qreal t0[3], t1[3];
cout << "Limit: " << limit << " degrees" << endl;
cout.flush();
......@@ -1569,7 +1569,7 @@ int Meshutils::divideSurfaceBySharpEdges(mesh_t *mesh)
//-----------------------------------------------------------------------------
void Meshutils::findSurfaceElementNormals(mesh_t *mesh)
{
static double a[3], b[3], c[3];
static qreal a[3], b[3], c[3];
double center_surface[3], center_element[3], center_difference[3];
Helpers *helpers = new Helpers;
int u, v, w, e0, e1, i0, i1, bigger;
......
......@@ -377,10 +377,10 @@ QPointF RenderArea::mapToRenderport(QPointF point) const
void RenderArea::fitSlot()
{
double xmin = 9e9;
double xmax = -9e9;
double ymin = 9e9;
double ymax = -9e9;
qreal xmin = 9e9;
qreal xmax = -9e9;
qreal ymin = 9e9;
qreal ymax = -9e9;
for(int i = 0; i < points.keys().size(); i++) {
int idx = points.keys().at(i);
......
......@@ -77,9 +77,21 @@ macx {
unix {
VTK_INCLUDEPATH = /usr/include/vtk-5.2
VTK_LIBPATH = /usr/lib
VTK_LIBS = -lvtkHybrid \
-lvtkWidgets \
-lQVTK
VTK_LIBS = -lQVTK \
-lvtkCommon \
-lvtkDICOMParser \
-lvtkFiltering \
-lvtkGenericFiltering \
-lvtkGraphics \
-lvtkHybrid \
-lvtkIO \
-lvtkImaging \
-lvtkInfovis \
-lvtkNetCDF \
-lvtkRendering \
-lvtkViews \
-lvtkVolumeRendering \
-lvtkWidgets
}
win32 {
......@@ -129,7 +141,23 @@ macx {
unix {
OCC_INCLUDEPATH = /usr/include/opencascade
OCC_LIBPATH = /usr/lib
OCC_LIBS = -lTKBRep -lTKSTL -lTKSTEP -lTKIGES
OCC_LIBS = -lTKSTL \
-lTKBRep \
-lTKernel \
-lTKG2d \
-lTKG3d \
-lTKGeomAlgo \
-lTKGeomBase \
-lTKMath \
-lTKMesh \
-lTKShHealing \
-lTKSTEP \
-lTKSTEP209 \
-lTKSTEPAttr \
-lTKSTEPBase \
-lTKIGES \
-lTKTopAlgo \
-lTKXSBase
}
win32 {
......
......@@ -233,7 +233,7 @@ PythonQtImporter_load_module(PyObject *obj, PyObject *args)
Py_DECREF(code);
if (Py_VerboseFlag)
PySys_WriteStderr("import %s # loaded from %s\n",
fullname, modpath);
fullname, (char*)modpath.toLatin1().data());
return mod;
error:
Py_DECREF(code);
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -706,6 +706,94 @@ else
fi
])dnl ACX_HYPRE
dnl
dnl @synopsis ACX_MUMPS([ACTION-IF-FOUND[, ACTION-IF-NOT-FOUND]])
dnl
dnl Look for MUMPS library
dnl
AC_DEFUN([ACX_MUMPS], [
AC_PREREQ(2.50)
AC_REQUIRE([ACX_MPI])
AC_REQUIRE([AC_FC_LIBRARY_LDFLAGS])
AC_REQUIRE([ACX_LANG_COMPILER_MS])
acx_mumps_ok=no
AC_ARG_WITH(mumps,
[AC_HELP_STRING([--with-mumps=<lib>], [Specify location of MUMPS])])
case $with_mumps in
yes | "") ;;
no) acx_mumps_ok=disable ;;
-* | */* | *.a | *.so | *.so.* | *.o) MUMPS_LIBS="$with_mumps" ;;
*) MUMPS_LIBS="-l$with_mumps" ;;
esac
# Get fortran linker names of MUMPS functions to check for.
#AC_FC_FUNC(mumps_function)
mumps_function="dmumps_"
acx_mumps_save_LIBS="$LIBS"
LIBS="$MPI_LIBS $LAPACK_LIBS $BLAS_LIBS $LIBS $FCLIBS $FLIBS"
# First, check MUMPS_LIBS environment variable
if test $acx_mumps_ok = no; then
if test "x$MUMPS_LIBS" != x; then
save_LIBS="$LIBS"; LIBS="$MUMPS_LIBS $LIBS"
AC_MSG_CHECKING([for $mumps_function in $MUMPS_LIBS])
if test "$acx_cv_c_compiler_ms" = "yes"; then
# windose shite
save_CFLAGS="$CFLAGS"
CFLAGS="$CFLAGS -Gz"
AC_LINK_IFELSE(
[int main ()
{
$mumps_function(1);
return 0;
}
],
[
acx_mumps_ok=yes
],
[
MUMPS_LIBS=""
])
CFLAGS="$save_CFLAGS"
else
AC_TRY_LINK_FUNC($mumps_function, [acx_mumps_ok=yes], [MUMPS_LIBS=""])
fi
AC_MSG_RESULT($acx_mumps_ok)
LIBS="$save_LIBS"
fi
fi
# General MUMPS test
if test $acx_mumps_ok = no; then
AC_CHECK_LIB(dmumps, $mumps_function, [acx_mumps_ok=yes; MUMPS_LIBS="-ldmumps $MPI_LIBS"],,[-lm])
if test $acx_mumps_ok = no; then
save_LIBS="$LIBS"; LIBS="-lmumps_common -lpord -lscalapack -lblacs $LIBS"
AC_CHECK_LIB(dmumps, $mumps_function, [acx_mumps_ok=yes; MUMPS_LIBS="-ldmumps -lmumps_common -lpord -lscalapack -lblacs $MPI_LIBS"],,[-lm])
LIBS="$save_LIBS"
fi
fi
AC_SUBST(MUMPS_LIBS)
LIBS="$acx_mumps_save_LIBS"
# Finally, execute ACTION-IF-FOUND/ACTION-IF-NOT-FOUND:
if test x"$acx_mumps_ok" = xyes; then
ifelse([$1],,AC_DEFINE(HAVE_MUMPS,1,[Define if you have a MUMPS library.]),[$1])
FCPPFLAGS="$FCPPFLAGS -DHAVE_MUMPS"
FCFLAGS="$FCFLAGS -I/usr/include"
:
else
acx_mumps_ok=no
$2
fi
])dnl ACX_MUMPS
dnl
dnl look for the std c libraries:
dnl
......
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.
......@@ -102,6 +102,7 @@ ACX_EIOF([], [AC_MSG_ERROR([libeiof $ELMER_LIBERRORMSG])])
ACX_ARPACK([], [AC_MSG_ERROR([We need arpack!])])
ACX_PARPACK([], [AC_MSG_WARN([No parallel arpack found.])])
ACX_HYPRE([], [AC_MSG_WARN([HYPRE not found, some functionaly will be disabled.])])
ACX_MUMPS([], [AC_MSG_WARN([MUMPS not found, some functionaly will be disabled.])])
ACX_UMFPACK([], [AC_MSG_WARN([UMFPACK not found, some functionaly will be disabled.])])
ACX_MATC([], [AC_MSG_ERROR([libmatc $ELMER_LIBERRORMSG])])
AC_CHECK_LIB(m,main,[UNIX_MATH_LIB="-lm"; LIBS="$LIBS -lm"],[UNIX_MATH_LIB=""])
......@@ -174,7 +175,7 @@ else
fi
fi
SOLVER_LIBS="$LIBS $ARPACK_LIBS $PARPACK_LIBS $HYPRE_LIBS $UMFPACK_LIBS $HUTI_LIBS $MATC_LIBS $EIOF_LIBS $LIBS $LAPACK_LIBS $BLAS_LIBS $FCLIBS $FLIBS $STDCXX_LIBS"
SOLVER_LIBS="$LIBS $ARPACK_LIBS $PARPACK_LIBS $HYPRE_LIBS $MUMPS_LIBS $UMFPACK_LIBS $HUTI_LIBS $MATC_LIBS $EIOF_LIBS $LIBS $LAPACK_LIBS $BLAS_LIBS $FCLIBS $FLIBS $STDCXX_LIBS"
dnl Host specific hacks.
dnl
......@@ -337,6 +338,7 @@ Libs:
BLAS $BLAS_LIBS
LAPACK $LAPACK_LIBS
HYPRE $HYPRE_LIBS
MUMPS $MUMPS_LIBS
UMFPACK $UMFPACK_LIBS
ARPACK $ARPACK_LIBS
PARPACK $PARPACK_LIBS
......
This diff is collapsed.
......@@ -5037,12 +5037,8 @@ END IF
Surface(2) = BoundaryNodes % y(1)
Surface(3) = 0.0_dp
t1(1) = BoundaryNodes % x(2) - Surface(1)
t1(2) = BoundaryNodes % y(2) - Surface(2)
t1(3) = 0.0_dp
Normal(1) = -t1(2)
Normal(2) = t1(1)
Normal(1) = Surface(2) - BoundaryNodes % y(2)
Normal(2) = BoundaryNodes % x(2) - Surface(1)
Normal(3) = 0.0_dp
END IF
......@@ -5059,7 +5055,41 @@ END IF
END FUNCTION LineFaceIntersection
!---------------------------------------------------------------------------
FUNCTION PointFaceDistance(BoundaryElement,BoundaryNodes,&
Coord,Normal,u0,v0) RESULT ( Dist )
!---------------------------------------------------------------------------
! This subroutine computes the signed distance of a point from a surface.
!---------------------------------------------------------------------------
TYPE(Nodes_t) :: BoundaryNodes
TYPE(Element_t), POINTER :: BoundaryElement
REAL(KIND=dp) :: Coord(3),Normal(3)
REAL(KIND=dp),OPTIONAL :: u0,v0
REAL(KIND=dp) :: Dist
REAL (KIND=dp) :: Surface(3),t1(3),t2(3),u,v
! For higher order elements this may be a necessity
IF( PRESENT( u0 ) .AND. PRESENT(v0) ) THEN
u = u0
v = v0
Surface = SurfaceVector( BoundaryElement, BoundaryNodes, u, v )
ELSE
u = 0.0_dp
v = 0.0_dp
! Any point known to be at the surface, even corner node
Surface(1) = BoundaryNodes % x(1)
Surface(2) = BoundaryNodes % y(1)
Surface(3) = BoundaryNodes % z(1)
END IF
Normal = NormalVector( BoundaryElement, BoundaryNodes, u, v, .TRUE. )
! Project of the line to the face normal
Dist = SUM( (Surface - Coord ) * Normal )
END FUNCTION PointFaceDistance
......
......@@ -549,6 +549,9 @@
CASE( 'user defined' )
CompressibilityModel = UserDefined1
CASE( 'pressure dependent' )
CompressibilityModel = UserDefined2
CASE( 'artificial compressible' )
CompressibilityModel = Incompressible
PseudoCompressible = .TRUE.
......@@ -730,6 +733,10 @@
END IF
END IF
CASE( UserDefined2 )
Density(1:n) = GetReal( Material,'Density' )
Pressure(1:nd) = FlowSolution(NSDOFs*FlowPerm(Indexes(1:nd)))
CASE( Thermal )
Pressure(1:n) = FlowSolution(NSDOFs*FlowPerm(Indexes(1:nd)))
......@@ -882,23 +889,17 @@
!------------------------------------------------------------------------------
SELECT CASE( CompressibilityModel )
!------------------------------------------------------------------------------
CASE( Incompressible,PerfectGas1,UserDefined1,Thermal)
CASE( Incompressible,PerfectGas1,UserDefined1,UserDefined2,Thermal)
!------------------------------------------------------------------------------
! Density needed for steady-state, also pressure for transient
!------------------------------------------------------------------------------
CALL NavierStokesCompose( &
MASS,STIFF,FORCE, LoadVector, &
Viscosity,Density,U,V,W,MU,MV,MW, &
ReferencePressure+Pressure(1:n), &
LocalTemperature, Convect, StabilizeFlag, &
CompressibilityModel == PerfectGas1, &
CompressibilityModel == Thermal .OR. &
CompressibilityModel == UserDefined1, &
PseudoCompressible, PseudoCompressibility, GasConstant, Porous, Drag, &
PotentialForce, PotentialField, PotentialCoefficient, &
MagneticForce, Rotating, AngularVelocity, &
DivDiscretization, GradPDiscretization, NewtonLinearization, &
Element,n,ElementNodes)
CALL NavierStokesCompose( MASS,STIFF,FORCE, LoadVector, &
Viscosity,Density,U,V,W,MU,MV,MW,ReferencePressure+Pressure(1:n), &
LocalTemperature, Convect, StabilizeFlag, CompressibilityModel, &
PseudoCompressible, PseudoCompressibility, GasConstant, Porous, &
Drag, PotentialForce, PotentialField, PotentialCoefficient, &
MagneticForce, Rotating, AngularVelocity, DivDiscretization, &
GradPDiscretization, NewtonLinearization, Element,n,ElementNodes)
!------------------------------------------------------------------------------
END SELECT
!------------------------------------------------------------------------------
......@@ -976,7 +977,7 @@
!------------------------------------------------------------------------------
! Neumann & Newton boundary conditions
!------------------------------------------------------------------------------
DO t = 1,Solver % Mesh % NumberOfBoundaryElements
DO t = 1,GetNOFBoundaryElements()
Element => GetBoundaryElement(t)
IF ( .NOT. ActiveBoundaryElement() ) CYCLE
......@@ -1090,8 +1091,8 @@
! * 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' )
Density(1:n) = GetParentMatProp( 'Density' )
Viscosity(1:n) = GetParentMatProp( 'Viscosity' )
LayerThickness(1:n) = GetReal( BC, 'Boundary Layer Thickness' )
SurfaceRoughness(1:n) = GetReal( BC, 'Surface Roughness',GotIt )
......@@ -1111,6 +1112,28 @@
CALL NavierStokesWallLaw( STIFF,FORCE, &
LayerThickness,SurfaceRoughness,Viscosity,Density,U,V,W, &
Element,n, ElementNodes )
ELSE IF ( GetLogical( BC, 'VMS Wall', GotIt ) ) THEN
Density(1:n) = GetParentMatProp( 'Density' )
Viscosity(1:n) = GetParentMatProp( 'Viscosity' )
LayerThickness(1:n) = GetReal( BC, 'Boundary Layer Thickness', GotIt )
SurfaceRoughness(1:n) = GetReal( BC, 'Surface Roughness',GotIt )
SELECT CASE( NSDOFs )
CASE(3)
U(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-2 )
V(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-1 )
W(1:n) = 0.0d0
CASE(4)
U(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-3 )
V(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-2 )
W(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-1 )
END SELECT
CALL VMSWalls( STIFF,FORCE, &
LayerThickness,SurfaceRoughness,Viscosity,Density,U,V,W, &
Element,n, ElementNodes )
END IF
!------------------------------------------------------------------------------
......
This diff is collapsed.
......@@ -872,7 +872,18 @@ CONTAINS
!------------------------------------------------------------------------------
CHARACTER(LEN=MAX_NAME_LEN) :: str
!------------------------------------------------------------------------------
str = ComponentName(Var % Name, Component)
IF ( Var % Name == 'flow solution' ) THEN
str='flow solution'
IF ( .NOT. PRESENT(Component) ) RETURN
IF ( Component == Var % DOFs ) THEN
str = 'pressure'
RETURN
ELSE
str = 'velocity ' // TRIM(IntegerToString(Component))
END IF
ELSE
str = ComponentName(Var % Name, Component)
END IF
!------------------------------------------------------------------------------
END FUNCTION ComponentNameVar
!------------------------------------------------------------------------------
......@@ -894,8 +905,9 @@ END FUNCTION ComponentNameVar
IF ( ind<=0 ) THEN
str = BaseName
IF ( Component > 0 ) &
IF ( Component > 0 ) THEN
str = TRIM(str) // ' ' // TRIM( IntegerToString(Component) )
END IF
ELSE
DOFsTot = 0
DO WHILE( .TRUE. )
......@@ -1625,6 +1637,14 @@ END FUNCTION ComponentNameVar
END FUNCTION NormalRandom
!---------------------------------------------------------
! Returns values from a even distribution [0,1]
!---------------------------------------------------------
FUNCTION EvenRandom() RESULT ( rand )
REAL(KIND=dp) :: rand
CALL RANDOM_NUMBER(rand)
END FUNCTION EvenRandom
END MODULE GeneralUtils
This diff is collapsed.
......@@ -274,7 +274,7 @@ EXTRA_DIST = \
NonphysicalMeshSolve.src \
NormalSolver.src \
ParallelEigenSolve.src \
ParallelUtils.src \
ParticleUtils.src \
PElementBase.src \
PElementMaps.src \
PoissonBEM.src \
......
......@@ -150,6 +150,7 @@ MKOCTFILE_DL_LDFLAGS = @MKOCTFILE_DL_LDFLAGS@
MPI_F90 = @MPI_F90@
MPI_INCLUDE_DIR = @MPI_INCLUDE_DIR@
MPI_LIBS = @MPI_LIBS@
MUMPS_LIBS = @MUMPS_LIBS@
MV = @MV@
NO_OCT_FILE_STRIP = @NO_OCT_FILE_STRIP@
OBJEXT = @OBJEXT@
......@@ -406,6 +407,7 @@ SOLVEROBJS = \
iso_varying_string$(OBJ_EXT) \
umf4_f77wrapper$(OBJ_EXT) \
VankaCreate$(OBJ_EXT) \
ParticleUtils$(OBJ_EXT) \
ElmerSolver$(OBJ_EXT)
NORMAL_TARGETS = mpif libelmersolver$(SHL_EXT) ElmerSolver$(EXE_EXT) $(BINARIES) ViewFactors$(EXE_EXT) GebhardtFactors$(EXE_EXT)
......@@ -500,7 +502,7 @@ EXTRA_DIST = \
NonphysicalMeshSolve.src \
NormalSolver.src \
ParallelEigenSolve.src \
ParallelUtils.src \
ParticleUtils.src \
PElementBase.src \
PElementMaps.src \
PoissonBEM.src \
......@@ -996,7 +998,7 @@ Komega.$(OBJEXT): Komega.f90 DefUtils.$(OBJEXT) MaterialModels.$(OBJEXT)
LevelSet.$(OBJEXT): LevelSet.f90 MaterialModels.$(OBJEXT) DefUtils.$(OBJEXT) Integration.$(OBJEXT) Lists.$(OBJEXT) ElementDescription.$(OBJEXT) Types.$(OBJEXT) SolverUtils.$(OBJEXT)
LinearAlgebra.$(OBJEXT): LinearAlgebra.f90 Types.$(OBJEXT)
ListMatrix.$(OBJEXT): ListMatrix.f90 CRSMatrix.$(OBJEXT) GeneralUtils.$(OBJEXT)
Lists.$(OBJEXT): Lists.f90 GeneralUtils.$(OBJEXT) Messages.$(OBJEXT) PElementMaps.$(OBJEXT) CoordinateSystems.$(OBJEXT) CRSMatrix.$(OBJEXT) Types.$(OBJEXT) Interpolation.$(OBJEXT)
Lists.$(OBJEXT): Lists.f90 GeneralUtils.$(OBJEXT) Messages.$(OBJEXT) PElementMaps.$(OBJEXT) CoordinateSystems.$(OBJEXT) CRSMatrix.$(OBJEXT) Types.$(OBJEXT) Interpolation.$(OBJEXT) SParIterComm.$(OBJEXT)
LUDecomposition.$(OBJEXT): LUDecomposition.f90 Types.$(OBJEXT)
MagneticSolve.$(OBJEXT): MagneticSolve.f90 MaxwellGeneral.$(OBJEXT) Maxwell.$(OBJEXT) Differentials.$(OBJEXT) DefUtils.$(OBJEXT) MaxwellAxiS.$(OBJEXT)
MagneticW1Solve.$(OBJEXT): MagneticW1Solve.f90 Differentials.$(OBJEXT) DirectSolve.$(OBJEXT) ElementUtils.$(OBJEXT) Integration.$(OBJEXT) Lists.$(OBJEXT) BandwidthOptimize.$(OBJEXT) IterSolve.$(OBJEXT) FreeSurface.$(OBJEXT) ElementDescription.$(OBJEXT) CoordinateSystems.$(OBJEXT) TimeIntegrate.$(OBJEXT) Types.$(OBJEXT) SolverUtils.$(OBJEXT)
......@@ -1020,6 +1022,7 @@ NonphysicalMeshSolve.$(OBJEXT): NonphysicalMeshSolve.f90 DefUtils.$(OBJEXT)
NormalSolver.$(OBJEXT): NormalSolver.f90 DefUtils.$(OBJEXT) CoordinateSystems.$(OBJEXT)
ParallelEigenSolve.$(OBJEXT): ParallelEigenSolve.f90 CRSMatrix.$(OBJEXT) IterSolve.$(OBJEXT) Multigrid.$(OBJEXT) ParallelUtils.$(OBJEXT)
ParallelUtils.$(OBJEXT): ParallelUtils.f90 SParIterSolver.$(OBJEXT)
ParticleUtils.$(OBJEXT): ParticleUtils.src DefUtils.$(OBJEXT)
PElementBase.$(OBJEXT): PElementBase.f90 PElementMaps.$(OBJEXT) Messages.$(OBJEXT) Types.$(OBJEXT)
PElementMaps.$(OBJEXT): PElementMaps.f90 Types.$(OBJEXT)
PhaseChangeSolve.$(OBJEXT): PhaseChangeSolve.f90 DefUtils.$(OBJEXT) Lists.$(OBJEXT) ElementDescription.$(OBJEXT) Types.$(OBJEXT)
......
......@@ -231,6 +231,7 @@ this ise not in USE
END INTERFACE
!------------------------------------------------------------------------------
mu = Viscosity
IF ( PRESENT(muder) ) muder=0
k = ListGetInteger( CurrentModel % Bodies(Element % BodyId) % Values, 'Material', &
minv=1, maxv=CurrentModel % NumberOFMaterials )
......@@ -265,7 +266,6 @@ this ise not in USE
ss=SecondInvariant(Velo,dVelodx,Metric,Symb)/2
IF ( PRESENT(muder) ) muder = 0._dp
SELECT CASE( ViscosityFlag )
CASE('power law')
......@@ -288,7 +288,9 @@ this ise not in USE
CASE('power law too')
c2n = ListGetReal( Material, 'Viscosity Exponent', n, Element % NodeIndexes )
c2 = SUM( Basis(1:n) * c2n(1:n) )
mu = Viscosity **(-1.0d0/c2)* ss**(-(c2-1)/(2*c2)) / 2
mu = Viscosity **(-1/c2)* ss**(-(c2-1)/(2*c2)) / 2
IF ( PRESENT(muder) ) muder = &
Viscosity**(-1/c2)*(-(c2-1)/(2*c2))*ss*(-(c2-1)/(2*c2)-1) / 2
CASE ('carreau')
c1n = ListGetReal( Material, 'Viscosity Difference',n,Element % NodeIndexes )
......@@ -300,9 +302,13 @@ this ise not in USE
c4 = ListGetConstReal( Material, 'Yasuda Coefficient',gotIt)
IF(gotIt) THEN
s = SQRT(ss)
mu = Viscosity + c1 * (1 + (c3*s)**c4)**((c2-1)/c4)
mu = Viscosity + c1 * (1 + c3**c4*ss**(c4/2))**((c2-1)/c4)
IF ( PRESENT(muder ) ) muder = &
c1*(1+c3**c4*ss**(c4/2))**((c2-1)/c4-1)*(c2-1)/2*c3**c4*ss**(c4/2-1)
ELSE
mu = Viscosity + c1 * (1 + (c3*c3*ss))**((c2-1)/2.0)
mu = Viscosity + c1 * (1 + c3*c3*ss)**((c2-1)/2)
IF ( PRESENT(muder) ) muder = &
c1*(c2-1)/2*c3**2*(1+c3**2*ss)**((c2-1)/2-1)
END IF
CASE ('cross')
......@@ -313,6 +319,8 @@ this ise not in USE
c3n = ListGetReal( Material, 'Viscosity Transition',n,Element % NodeIndexes )
c3 = SUM( Basis(1:n) * c3n(1:n) )
mu = Viscosity + c1 / (1 + c3*ss**(c2/2))
IF ( PRESENT(muder) ) muder = &
-c1*c3*ss**(c2/2)*c2 / (2*(1+c3*ss**(c2/2))**2*ss)
CASE ('powell eyring')
c1n = ListGetReal( Material, 'Viscosity Difference',n,Element % NodeIndexes)
......@@ -323,6 +331,9 @@ this ise not in USE
mu = Viscosity + c1
ELSE
mu = Viscosity + c1 * LOG(c2*s+SQRT(c2*c2*ss+1))/(c2*s)
IF ( PRESENT(muder) ) muder = &
c1*(c2/(2*s)+c2**2/(2*SQRT(c2**2*ss+1)))/((c2*s+SQRT(c2*ss+1))*c2*s) - &
c1*LOG(c2*s+SQRT(c2**2*ss+1))/(c2*s**3)/2
END IF
CASE ('thermal carreau')
......@@ -352,8 +363,12 @@ this ise not in USE
IF(gotIt) THEN
s = SQRT(ss)
mu = Viscosity + c1 * (1.0d0 + (c3*s)**c4)**((c2-1)/c4)
IF ( PRESENT(muder) ) muder = &
c1*(1+(c3*s)**c4)**((c2-1)/c4-1)*(c2-1)*(c3*s)**c4/(2*ss)
ELSE
mu = Viscosity + c1 * (1.0d0 + (c3*c3*ss))**((c2-1)/2.0)
IF ( PRESENT(muder) ) muder = &
c1*(1+c3**2*ss)**((c2-1)/2-1)*(c2-1)/2*c3**2
END IF
CASE( 'smagorinsky' )
......@@ -361,7 +376,9 @@ this ise not in USE
n, Element % NodeIndexes,gotit )
c2 = SUM( Basis(1:n) * c2n(1:n) )
h = ElementDiameter( Element, Nodes )
Viscosity = Viscosity + Density * c2 * h**2 * SQRT( 2*ss ) / 2
Viscosity = Viscosity + Density * c2 * h**2 * SQRT(2*ss) / 2
IF ( PRESENT(muder) ) muder = &
Density*c2*h**2*SQRT(2._dp)/(4*SQRT(ss))
CASE( 'ke','k-epsilon' )
IF (ListGetString(Material,'KE Model',gotIt)/='v2-f' ) THEN
......