Index: /issm/branches/trunk-jpl-damage/configs/config-arm-linux.sh =================================================================== --- /issm/branches/trunk-jpl-damage/configs/config-arm-linux.sh (revision 12945) +++ /issm/branches/trunk-jpl-damage/configs/config-arm-linux.sh (revision 12946) @@ -9,3 +9,4 @@ --without-modules \ --with-gsl-dir=$ISSM_DIR/externalpackages/gsl/install \ - --with-triangle-dir=$ISSM_DIR/externalpackages/triangle/install + --with-triangle-dir=$ISSM_DIR/externalpackages/triangle/install \ + --disable-static Index: /issm/branches/trunk-jpl-damage/configure.ac =================================================================== --- /issm/branches/trunk-jpl-damage/configure.ac (revision 12945) +++ /issm/branches/trunk-jpl-damage/configure.ac (revision 12946) @@ -21,18 +21,14 @@ AM_SILENT_RULES([yes]) #Do not show compilation command by default AM_PROG_CC_C_O - AM_PROG_AR -AC_PROG_RANLIB +AC_PROG_RANLIB # To produce shared libraries uncomment LT_INIT. In theory, Libtool can be initialized # even if all linking is static, since, if not shared libraries are to be produce, libtool # will act simply as a wrapper around AR and Ranlib. -# #LT_INIT #Run issm_options.m4 ISSM_OPTIONS - -#AM_PROG_AR must be placed after ISSM_OPTIONS #List all Makefiles Index: /issm/branches/trunk-jpl-damage/etc/environment.sh =================================================================== --- /issm/branches/trunk-jpl-damage/etc/environment.sh (revision 12945) +++ /issm/branches/trunk-jpl-damage/etc/environment.sh (revision 12946) @@ -278,4 +278,7 @@ export LD_RUN_PATH="$LD_RUN_PATH:$TCLX_DIR" +#TCL/TK +export LD_LIBRARY_PATH="$LD_LIBRARY_PATH:/usr/local/lib" +export LD_RUN_PATH="$LD_RUN_PATH:/usr/local/lib" #ASPELL Index: /issm/branches/trunk-jpl-damage/externalpackages/gsl/install-android.sh =================================================================== --- /issm/branches/trunk-jpl-damage/externalpackages/gsl/install-android.sh (revision 12945) +++ /issm/branches/trunk-jpl-damage/externalpackages/gsl/install-android.sh (revision 12946) @@ -10,5 +10,5 @@ #Download from ISSM server - $ISSM_DIR/scripts/DownloadExternalPackage.py 'http://issm.jpl.nasa.gov/files/externalpackages/gsl-1.15.tar.gz' 'gsl-1.15.tar.gz' + $ISSM_DIR/scripts/DownloadExternalPackage.py 'http://issm.jpl.nasa.gov/files/externalpackages/gsl-1.' 'gsl-1.15.tar.gz' #Untar @@ -28,6 +28,7 @@ ./configure \ --build="i386-apple-darwin10.8.0" \ - --host=$host_triplet\ - --prefix="$ISSM_DIR/externalpackages/gsl/install/" + --host=$host_triplet \ + --prefix="$ISSM_DIR/externalpackages/gsl/install/" \ + --disable-static fi Index: /issm/branches/trunk-jpl-damage/externalpackages/libtool/install.sh =================================================================== --- /issm/branches/trunk-jpl-damage/externalpackages/libtool/install.sh (revision 12945) +++ /issm/branches/trunk-jpl-damage/externalpackages/libtool/install.sh (revision 12946) @@ -6,13 +6,12 @@ #Download from ISSM server -$ISSM_DIR/scripts/DownloadExternalPackage.py 'http://issm.jpl.nasa.gov/files/externalpackages/libtool' 'libtool' +$ISSM_DIR/scripts/DownloadExternalPackage.py 'http://issm.jpl.nasa.gov/files/externalpackages/libtool-2.4.2.tar.gz' 'libtool-2.4.2.tar.gz' #Untar -tar -zxvf libtool[\w.-]*.tar.gz -rm libtool[\w.-]*.tar.gz - +tar -zxvf libtool-2.4.2.tar.gz +rm libtool-2.4.2.tar.gz #Move libtool into src directory -mv libtool[\w.-]* src +mv libtool-2.4.2 src #Compile libtool Index: /issm/branches/trunk-jpl-damage/externalpackages/matplotlib/install-osx.sh =================================================================== --- /issm/branches/trunk-jpl-damage/externalpackages/matplotlib/install-osx.sh (revision 12945) +++ /issm/branches/trunk-jpl-damage/externalpackages/matplotlib/install-osx.sh (revision 12946) @@ -2,4 +2,12 @@ pythonversion=2.7 PYTHON=python${pythonversion} + +export GIT_SSL_NO_VERIFY=true +export CC="gcc -fPIC" +export CXX="g++ -fPIC" +export F77="gfortran -fPIC" +export FC="gfortran -fPIC" +export FFLAGS=-ff2c + git clone https://github.com/matplotlib/matplotlib @@ -7,4 +15,5 @@ mkdir deps cd src +exit make -f make.osx PREFIX=$ISSM_DIR/externalpackages/matplotlib/deps PYVERSION=$pythonversion fetch deps mpl_install_std ${PYTHON} -c "import matplotlib; print 'Installed matplotlib', matplotlib.__version__, matplotlib.__file__" Index: /issm/branches/trunk-jpl-damage/externalpackages/matplotlib/install.sh =================================================================== --- /issm/branches/trunk-jpl-damage/externalpackages/matplotlib/install.sh (revision 12945) +++ /issm/branches/trunk-jpl-damage/externalpackages/matplotlib/install.sh (revision 12946) @@ -1,7 +1,17 @@ #/bin/bash +pythonversion=2.7 +PYTHON=python${pythonversion} + +export GIT_SSL_NO_VERIFY=true +export CC="gcc -fPIC " +export CXX="g++ -fPIC -L$ISSM_DIR/externalpackages/tcl/install/lib" +export F77="gfortran -fPIC" +export FC="gfortran -fPIC" +export FFLAGS=-ff2c + git clone https://github.com/matplotlib/matplotlib mv matplotlib src cd src -python setup.py build -python setup.py install +python setup.py build +python setup.py install Index: /issm/branches/trunk-jpl-damage/externalpackages/mpich2/install-1.3.1-linux64-berg.sh =================================================================== --- /issm/branches/trunk-jpl-damage/externalpackages/mpich2/install-1.3.1-linux64-berg.sh (revision 12946) +++ /issm/branches/trunk-jpl-damage/externalpackages/mpich2/install-1.3.1-linux64-berg.sh (revision 12946) @@ -0,0 +1,51 @@ +#!/bin/bash + +#Get number of cpus +NUMCPUS=$1; + +#version of mpich2 +#version=1.0.2p1 +version=1.3.1 +#version=1.4 + +#Some cleanup +rm -rf src install mpich2-$version + +#Create src and install directories +mkdir src install + +#Untar +tar -zxvf mpich2-$version.tar.gz + +#Move mpich2 into src directory +mv mpich2-$version/* src +rm -rf mpich2-$version + +#Configure mpich2 +cd src +if [[ "$version" == "1.0.2p1" ]]; +then + make distclean +fi + +source ../configure.sh + +if [[ "$version" == "1.0.2p1" ]]; +then + #Apply Patch to examples + patch ./examples/Makefile ../Examples.Makefile.patch + patch -R ./src/binding/cxx/mpicxx.h ../mpicxx.h.patch + patch -R ./src/include/mpicxx.h ../mpicxx.h.patch +fi + +#Compile mpich2 +if [ -z $NUMCPUS ]; +then + make +else + make -j $NUMCPUS +fi +make install + +cd ../install/lib +rm -rf *.so Index: /issm/branches/trunk-jpl-damage/externalpackages/mpich2/install-1.4-linux64.sh =================================================================== --- /issm/branches/trunk-jpl-damage/externalpackages/mpich2/install-1.4-linux64.sh (revision 12945) +++ /issm/branches/trunk-jpl-damage/externalpackages/mpich2/install-1.4-linux64.sh (revision 12946) @@ -32,6 +32,2 @@ fi make install - -#remove so files -cd ../install/lib -rm -rf *.so Index: /issm/branches/trunk-jpl-damage/externalpackages/mpich2/install-1.4-macosx32.sh =================================================================== --- /issm/branches/trunk-jpl-damage/externalpackages/mpich2/install-1.4-macosx32.sh (revision 12945) +++ /issm/branches/trunk-jpl-damage/externalpackages/mpich2/install-1.4-macosx32.sh (revision 12946) @@ -35,6 +35,2 @@ fi make install - -#remove so files -cd ../install/lib -rm -rf *.so Index: /issm/branches/trunk-jpl-damage/externalpackages/mpich2/install-1.4-macosx64.sh =================================================================== --- /issm/branches/trunk-jpl-damage/externalpackages/mpich2/install-1.4-macosx64.sh (revision 12945) +++ /issm/branches/trunk-jpl-damage/externalpackages/mpich2/install-1.4-macosx64.sh (revision 12946) @@ -35,6 +35,2 @@ fi make install - -#remove so files -cd ../install/lib -rm -rf *.so Index: /issm/branches/trunk-jpl-damage/externalpackages/tcl/install.sh =================================================================== --- /issm/branches/trunk-jpl-damage/externalpackages/tcl/install.sh (revision 12945) +++ /issm/branches/trunk-jpl-damage/externalpackages/tcl/install.sh (revision 12946) @@ -1,3 +1,5 @@ #!/bin/bash + +sudochoice=0; #Some cleanup @@ -20,4 +22,27 @@ ver="8.4.12" -cd src -make -C ./macosx install INSTALL_ROOT="$ISSM_DIR/externalpackages/tcl/install" +cd src/unix + +#User mode: +if [[ $sudochoice == "0" ]]; +then + ./configure --prefix=$ISSM_DIR/externalpackages/tcl/install + if [ -z $1 ]; then + make + else + make -j $1 + fi + make install +fi + +#sudo mode: +if [[ $sudochoice == "1" ]]; +then + sudo ./configure + if [ -z $1 ]; then + sudo make + else + sudo make -j $1 + fi + sudo make install +fi Index: /issm/branches/trunk-jpl-damage/externalpackages/tk/install.sh =================================================================== --- /issm/branches/trunk-jpl-damage/externalpackages/tk/install.sh (revision 12946) +++ /issm/branches/trunk-jpl-damage/externalpackages/tk/install.sh (revision 12946) @@ -0,0 +1,45 @@ +#!/bin/bash + +sudochoice=0; + +#Some cleanup +rm -rf src +rm -rf install +rm -rf tk8.5.12 +mkdir src install + +#Download from ISSM server +$ISSM_DIR/scripts/DownloadExternalPackage.py 'http://issm.jpl.nasa.gov/files/externalpackages/tk8.5.12-src.tar.gz' 'tk8.5.12.tar.gz' + +#Untar +tar -zxvf tk8.5.12.tar.gz + +#Move tk into src directory +mv tk8.5.12/* src +rm -rf tk8.5.12 + +cd src/unix + +#User mode: +if [[ $sudochoice == "0" ]]; +then + ./configure --prefix=$ISSM_DIR/externalpackages/tk/install + if [ -z $1 ]; then + make + else + make -j $1 + fi + make install +fi + +#sudo mode: +if [[ $sudochoice == "1" ]]; +then + sudo ./configure + if [ -z $1 ]; then + sudo make + else + sudo make -j $1 + fi + sudo make install +fi Index: /issm/branches/trunk-jpl-damage/externalpackages/triangle/install-android.sh =================================================================== --- /issm/branches/trunk-jpl-damage/externalpackages/triangle/install-android.sh (revision 12945) +++ /issm/branches/trunk-jpl-damage/externalpackages/triangle/install-android.sh (revision 12946) @@ -1,34 +1,26 @@ #!/bin/bash -source $ANDROID_DIR/android_aux.sh -# use matlab? -matlab=0 - -# Some cleanup +#Some cleanup rm -rf install triangle mkdir install -# Download from ISSM server +#Download from ISSM server $ISSM_DIR/scripts/DownloadExternalPackage.py 'http://issm.jpl.nasa.gov/files/externalpackages/triangle.zip' 'triangle.zip' -# Untar +#Untar cd install cp ../triangle.zip ./ unzip triangle.zip -# copy new makefile -cp ../configs//android/configure.make ./ +#copy new makefile +cp ../configs/android/configure.make ./ cp ../makefile ./ -# Patch triangle.c -if [[ $matlab == "1" ]];then - patch triangle.c ../triangle.c.patch.matlab -else - patch triangle.c ../triangle.c.patch.python -fi +#Patch triangle.c +patch triangle.c ../triangle.c.patch -# Compile triangle -make -j $j +#Compile triangle +make -# Patch triangle.h +#Patch triangle.h patch triangle.h ../triangle.h.patch Index: /issm/branches/trunk-jpl-damage/externalpackages/vim/addons/vim/syntax/Synchronize.sh =================================================================== --- /issm/branches/trunk-jpl-damage/externalpackages/vim/addons/vim/syntax/Synchronize.sh (revision 12945) +++ /issm/branches/trunk-jpl-damage/externalpackages/vim/addons/vim/syntax/Synchronize.sh (revision 12946) @@ -15,6 +15,5 @@ cat c.vim | sed "/ISSM's objects begin/,/vim: ts=8/d" > temp echo "\"ISSM's objects begin" >> temp -cat ../../../../../src/c/objects/objects.h | grep "\.h" | sed -e "s/\// /g" -e "s/\.h\"//" | awk '{print $(NF)}' | awk '{ printf "syn keyword cType " $1 "\n"}' >> temp -echo "\"ISSM's objects end" >> temp +find ../../../../../src/c/classes -name "*.h" | sed -e "s/\// /g" -e "s/\.h//"| awk '{print $(NF)}' | awk '{ printf "syn keyword cType " $1 "\n"}'>> temp cat c.vim | sed "1,/ISSM's objects end/d" >> temp Index: /issm/branches/trunk-jpl-damage/externalpackages/vim/addons/vim/syntax/c.vim =================================================================== --- /issm/branches/trunk-jpl-damage/externalpackages/vim/addons/vim/syntax/c.vim (revision 12945) +++ /issm/branches/trunk-jpl-damage/externalpackages/vim/addons/vim/syntax/c.vim (revision 12946) @@ -411,27 +411,66 @@ "ISSM's objects begin -syn keyword cType objects.h -syn keyword cType Object +syn keyword cType AdjacentTriangle +syn keyword cType BamgGeom +syn keyword cType BamgMesh +syn keyword cType bamgobjects +syn keyword cType BamgOpts +syn keyword cType BamgQuadtree +syn keyword cType BamgVertex +syn keyword cType CrackedEdge +syn keyword cType Curve +syn keyword cType Direction +syn keyword cType DoubleAndInt +syn keyword cType Edge +syn keyword cType GeomEdge +syn keyword cType Geometry +syn keyword cType GeomSubDomain +syn keyword cType GeomVertex +syn keyword cType include +syn keyword cType ListofIntersectionTriangles +syn keyword cType macros +syn keyword cType Mesh +syn keyword cType Metric +syn keyword cType R2 +syn keyword cType SetOfE4 +syn keyword cType SubDomain +syn keyword cType Triangle +syn keyword cType typedefs +syn keyword cType VertexOnEdge +syn keyword cType VertexOnGeom +syn keyword cType VertexOnVertex +syn keyword cType classes +syn keyword cType DakotaPlugin +syn keyword cType DofIndexing +syn keyword cType FemModel +syn keyword cType gaussobjects +syn keyword cType GaussPenta +syn keyword cType GaussTria syn keyword cType Hook -syn keyword cType DofIndexing -syn keyword cType Vertex -syn keyword cType Node -syn keyword cType Segment syn keyword cType IoModel -syn keyword cType Patch -syn keyword cType Update +syn keyword cType ExponentialVariogram +syn keyword cType GaussianVariogram +syn keyword cType krigingobjects +syn keyword cType Observation +syn keyword cType PowerVariogram +syn keyword cType Quadtree +syn keyword cType SphericalVariogram +syn keyword cType Variogram +syn keyword cType ElementMatrix +syn keyword cType ElementVector +syn keyword cType Matrix +syn keyword cType matrixobjects +syn keyword cType Vector syn keyword cType Constraint +syn keyword cType SpcDynamic syn keyword cType SpcStatic syn keyword cType SpcTransient -syn keyword cType SpcDynamic -syn keyword cType GaussTria -syn keyword cType GaussPenta -syn keyword cType Load -syn keyword cType Friction -syn keyword cType Icefront -syn keyword cType Numericalflux -syn keyword cType Riftfront -syn keyword cType Penpair -syn keyword cType Pengrid +syn keyword cType Contour +syn keyword cType BoolElementResult +syn keyword cType DoubleElementResult +syn keyword cType ElementResult +syn keyword cType ElementResultLocal +syn keyword cType PentaP1ElementResult +syn keyword cType TriaP1ElementResult syn keyword cType Element syn keyword cType Penta @@ -441,7 +480,24 @@ syn keyword cType TriaHook syn keyword cType TriaRef +syn keyword cType BoolExternalResult +syn keyword cType DoubleExternalResult +syn keyword cType DoubleMatExternalResult +syn keyword cType DoubleVecExternalResult +syn keyword cType ExternalResult +syn keyword cType IntExternalResult +syn keyword cType PetscVecExternalResult +syn keyword cType StringExternalResult +syn keyword cType BoolInput +syn keyword cType ControlInput +syn keyword cType DatasetInput +syn keyword cType DoubleInput +syn keyword cType Input +syn keyword cType IntInput +syn keyword cType PentaP1Input +syn keyword cType TransientInput +syn keyword cType TriaP1Input syn keyword cType KML_Attribute +syn keyword cType KML_ColorStyle syn keyword cType KML_Comment -syn keyword cType KML_ColorStyle syn keyword cType KML_Container syn keyword cType KML_Document @@ -459,6 +515,6 @@ syn keyword cType KML_Object syn keyword cType KML_Overlay +syn keyword cType KML_Placemark syn keyword cType KML_Point -syn keyword cType KML_Placemark syn keyword cType KML_Polygon syn keyword cType KML_PolyStyle @@ -468,57 +524,43 @@ syn keyword cType KML_Unknown syn keyword cType KMLFileReadUtils -syn keyword cType Option -syn keyword cType OptionDouble -syn keyword cType OptionLogical -syn keyword cType OptionChar -syn keyword cType OptionStruct -syn keyword cType OptionCell -syn keyword cType OptionUtilities -syn keyword cType Input -syn keyword cType BoolInput -syn keyword cType DoubleInput -syn keyword cType IntInput -syn keyword cType PentaP1Input -syn keyword cType TriaP1Input -syn keyword cType ControlInput -syn keyword cType DatasetInput -syn keyword cType TransientInput -syn keyword cType ElementResult -syn keyword cType DoubleElementResult -syn keyword cType TriaP1ElementResult -syn keyword cType PentaP1ElementResult -syn keyword cType BoolElementResult -syn keyword cType ExternalResult -syn keyword cType BoolExternalResult -syn keyword cType DoubleExternalResult -syn keyword cType DoubleVecExternalResult -syn keyword cType DoubleMatExternalResult -syn keyword cType IntExternalResult -syn keyword cType PetscVecExternalResult -syn keyword cType StringExternalResult +syn keyword cType Friction +syn keyword cType Icefront +syn keyword cType Load +syn keyword cType Numericalflux +syn keyword cType Pengrid +syn keyword cType Penpair +syn keyword cType Riftfront syn keyword cType Material syn keyword cType Matice syn keyword cType Matpar -syn keyword cType ElementMatrix -syn keyword cType ElementVector -syn keyword cType Vector -syn keyword cType Matrix +syn keyword cType Node +syn keyword cType Object +syn keyword cType objects +syn keyword cType Option +syn keyword cType OptionCell +syn keyword cType OptionChar +syn keyword cType OptionDouble +syn keyword cType OptionLogical +syn keyword cType OptionStruct +syn keyword cType OptionUtilities syn keyword cType BoolParam +syn keyword cType DoubleMatArrayParam syn keyword cType DoubleMatParam +syn keyword cType DoubleParam syn keyword cType DoubleTransientMatParam -syn keyword cType DoubleMatArrayParam -syn keyword cType DoubleParam syn keyword cType DoubleVecParam +syn keyword cType FileParam +syn keyword cType IntMatParam syn keyword cType IntParam syn keyword cType IntVecParam -syn keyword cType IntMatParam -syn keyword cType FileParam +syn keyword cType MatrixParam syn keyword cType Param -syn keyword cType MatrixParam -syn keyword cType VectorParam syn keyword cType StringArrayParam syn keyword cType StringParam -syn keyword cType Contour -syn keyword cType FemModel +syn keyword cType TransientParam +syn keyword cType VectorParam +syn keyword cType Segment +syn keyword cType Vertex +syn keyword cType objects syn keyword cType OptArgs syn keyword cType OptPars Index: /issm/branches/trunk-jpl-damage/m4/issm_options.m4 =================================================================== --- /issm/branches/trunk-jpl-damage/m4/issm_options.m4 (revision 12945) +++ /issm/branches/trunk-jpl-damage/m4/issm_options.m4 (revision 12946) @@ -173,42 +173,46 @@ MATLABINCL=-I"$MATLAB_ROOT/extern/include"; - dnl OS-dependent variables and checks - case "${host_os}" in - *linux*) - if test "${host_cpu}" = "x86_64"; - then - MEXLIB=-L"$MATLAB_ROOT/bin/glnxa64/ -lmex" - MEXLINK="-pthread -shared -W2,--version-script,${MATLAB_ROOT}/extern/lib/glnxa64/mexFunction.map"; - else - MEXLIB=-L"$MATLAB_ROOT/bin/glnx86/ -lmex" - MEXLINK="-pthread -shared -W2,--version-script,${MATLAB_ROOT}/extern/lib/glnx86/mexFunction.map"; - fi - MEXEXT=`$MATLAB_ROOT/bin/mexext` - MEXEXT=".$MEXEXT" - ;; - *darwin*) - dnl mex -v gives all the flags for compilation of mex files - dnl if matlab version is 7.9 or more, we must use mexmaci64 (64 bits) - MEXLINK="-O -Wl,-flat_namespace -undefined suppress -arch i386 -bundle -Wl,-exported_symbols_list,$MATLAB_ROOT/extern/lib/maci/mexFunction.map" - MEXLIB=" -L$MATLAB_ROOT/bin/maci/ -lmx -lmex -lmat -lstdc++ -largeArrayDims" - if test $MATLAB_MAJOR -ge 7; then - if test $MATLAB_MINOR -ge 9; then - MEXLINK="-O -Wl,-flat_namespace -undefined suppress -bundle -Wl,-exported_symbols_list,$MATLAB_ROOT/extern/lib/maci64/mexFunction.map" - MEXLIB=" -L$MATLAB_ROOT/bin/maci64/ -lmx -lmex -lmat -lstdc++" - fi - fi - MEXEXT=`$MATLAB_ROOT/bin/mexext` - MEXEXT=".$MEXEXT" - ;; - *cygwin*) - if test $VENDOR = intel-win7-32; then - MEXLIB="-dll -export:mexFunction -LIBPATH:\"$MATLAB_ROOT\extern\lib\win32\microsoft\" libmx.lib libmex.lib libmat.lib kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib -NOLOGO -INCREMENTAL:NO -manifest" - elif test $VENDOR = intel-win7-64; then - MEXLIB="-dll -export:mexFunction -LIBPATH:\"$MATLAB_ROOT\extern\lib\win64\microsoft\" libmx.lib libmex.lib libmat.lib kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib -NOLOGO -INCREMENTAL:NO -manifest" - fi - MEXEXT=`$MATLAB_ROOT/bin/mexext.bat` - MEXEXT=".$MEXEXT" - ;; - esac + dnl 4. get MEXLIB MEXLINK and MEXEXT (experimental) + MEXLINK=$(mex -v 2>&1 < /dev/null | grep LDFLAGS | sed -e "s/ LDFLAGS = //g") + MEXLIB=$( mex -v 2>&1 < /dev/null | grep CXXLIBS | sed -e "s/ CXXLIBS = //g") + MEXEXT=$( mex -v 2>&1 < /dev/null | grep LDEXTENSION | sed -e "s/ LDEXTENSION = //g") + dnl dnl OS-dependent variables and checks (old stuff) + dnl case "${host_os}" in + dnl *linux*) + dnl if test "${host_cpu}" = "x86_64"; + dnl then + dnl MEXLIB="-Wl,-rpath-link,$MATLAB_ROOT/bin/glnxa64 -L$MATLAB_ROOT/bin/glnxa64/ -lmx -lmex -lmat -lm" + dnl MEXLINK="-pthread -shared -W2,--version-script,${MATLAB_ROOT}/extern/lib/glnxa64/mexFunction.map"; + dnl else + dnl MEXLIB=-L"$MATLAB_ROOT/bin/glnx86/ -lmex" + dnl MEXLINK="-pthread -shared -W2,--version-script,${MATLAB_ROOT}/extern/lib/glnx86/mexFunction.map"; + dnl fi + dnl MEXEXT=`$MATLAB_ROOT/bin/mexext` + dnl MEXEXT=".$MEXEXT" + dnl ;; + dnl *darwin*) + dnl dnl mex -v gives all the flags for compilation of mex files + dnl dnl if matlab version is 7.9 or more, we must use mexmaci64 (64 bits) + dnl MEXLINK="-O -Wl,-flat_namespace -undefined suppress -arch i386 -bundle -Wl,-exported_symbols_list,$MATLAB_ROOT/extern/lib/maci/mexFunction.map" + dnl MEXLIB=" -L$MATLAB_ROOT/bin/maci/ -lmx -lmex -lmat -lstdc++ -largeArrayDims" + dnl if test $MATLAB_MAJOR -ge 7; then + dnl if test $MATLAB_MINOR -ge 9; then + dnl MEXLINK="-O -Wl,-flat_namespace -undefined suppress -bundle -Wl,-exported_symbols_list,$MATLAB_ROOT/extern/lib/maci64/mexFunction.map" + dnl MEXLIB=" -L$MATLAB_ROOT/bin/maci64/ -lmx -lmex -lmat -lstdc++" + dnl fi + dnl fi + dnl MEXEXT=`$MATLAB_ROOT/bin/mexext` + dnl MEXEXT=".$MEXEXT" + dnl ;; + dnl *cygwin*) + dnl if test $VENDOR = intel-win7-32; then + dnl MEXLIB="-dll -export:mexFunction -LIBPATH:\"$MATLAB_ROOT\extern\lib\win32\microsoft\" libmx.lib libmex.lib libmat.lib kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib -NOLOGO -INCREMENTAL:NO -manifest" + dnl elif test $VENDOR = intel-win7-64; then + dnl MEXLIB="-dll -export:mexFunction -LIBPATH:\"$MATLAB_ROOT\extern\lib\win64\microsoft\" libmx.lib libmex.lib libmat.lib kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib -NOLOGO -INCREMENTAL:NO -manifest" + dnl fi + dnl MEXEXT=`$MATLAB_ROOT/bin/mexext.bat` + dnl MEXEXT=".$MEXEXT" + dnl ;; + dnl esac AC_MSG_RESULT(done) @@ -501,5 +505,5 @@ HAVE_GSL=yes GSLINCL="-I$GSL_ROOT/include" - GSLLIB="-dy -L$GSL_ROOT/lib/ -lgsl -lgslcblas -lm" + GSLLIB="-dy -L$GSL_ROOT/lib -lgsl -lgslcblas -lm" AC_DEFINE([_HAVE_GSL_],[1],[with gsl in ISSM src]) @@ -671,5 +675,5 @@ PETSCLIB="-L$PETSC_ROOT/$PETSC_ARCH/lib -L$PETSC_ROOT/lib -L$PETSC_ROOT/lib/$PETSC_ARCH -lpetscksp -lpetscdm -lpetscmat -lpetscvec -lpetsc -lpetscsnes -lpetscts" else - PETSCLIB="-L$PETSC_ROOT/$PETSC_ARCH/lib -L$PETSC_ROOT/lib -L$PETSC_ROOT/lib/$PETSC_ARCH -lpetsc" + PETSCLIB="-L$PETSC_ROOT/$PETSC_ARCH/lib -L$PETSC_ROOT/lib -L$PETSC_ROOT/lib/$PETSC_ARCH -lpetsc -ldl" fi ;; Index: /issm/branches/trunk-jpl-damage/scripts/DownloadExternalPackage.py =================================================================== --- /issm/branches/trunk-jpl-damage/scripts/DownloadExternalPackage.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/scripts/DownloadExternalPackage.py (revision 12946) @@ -1,4 +1,8 @@ #!/usr/bin/env python # -*- coding: ISO-8859-1 -*- +# +# TODO: Improve upon extension recognition by checking for mismatches in found targets +# and specified local file. +# import os,sys,re @@ -21,18 +25,78 @@ # End class myHTMLParser -pivot = sys.argv[1].rfind("/") -url = (sys.argv[1])[:pivot] -localFile=sys.argv[2] +# Separates the URL into a directory and the file or pattern based on the +# last appearance of '/'. +if len(sys.argv) > 1: + pivot = sys.argv[1].rfind("/") + url = (sys.argv[1])[:pivot] + pivot += 1 + find = (sys.argv[1])[pivot:] +else: + print "******************************************************************************************************************************" + print "* Invalid input! *" + print "* *" + print "* Try: 'DownloadExternalPackage.py url [localFile]' *" + print "* *" + print "* Where 'URL' is the URL with an explicit package name or the URL followed by the truncated package name. And 'localFile' is *" + print "* the file name (including extension) that you would like to save as. *" + print "* *" + print "* Examples: *" + print "* *" + print "* DownloadExternalPackage.py 'http://issm.jpl.nasa.gov/files/externalpackages/petsc-2.3.2-p3.tar.gz' 'petsc-2.3.2-p3.tar.gz' *" + print "* *" + print "* This is the old style and the safest way to download a package. *" + print "* *" + print "* DownloadExternalPackage.py 'http://issm.jpl.nasa.gov/files/externalpackages/libtool' 'libtool.tar.gz' *" + print "* *" + print "* This is the new style. For packages like 'Libtool', which we never expect to be using multiple versions, this will *" + print "* download the most recent version and save it as the generic 'libtool.tar.gz'. *" + print "* *" + print "* DownloadExternalPackage.py 'http://issm.jpl.nasa.gov/files/externalpackages/gsl-1.' 'gsl-1.15.tar.gz' *" + print "* *" + print "* This is the new style. This is a demonstration of how this script can be used to disambiguate a package name if there *" + print "* are more than once package matching 'gsl-'. *" + print "* *" + print "* DownloadExternalPackage.py 'http://issm.jpl.nasa.gov/files/externalpackages/libtool' *" + print "* *" + print "* This is the new style. This will download a package with 'libtool' as a prefix and save it as its canonical name. *" + print "* *" + print "* *" + print "******************************************************************************************************************************" -pivot += 1; -find = (sys.argv[1])[pivot:]; +if len(sys.argv) > 2: + localFile=sys.argv[2] + print "Downloaded file will saved as: " + localFile +else: + localFile = None + print "Downloaded file will saved with the same file name." + + print "Looking for " + find + +# As an extra precaution, if no extension is given for a particular package +# such as '.../libtool', then ensure that files found are of appropriate +# file extensions. +# +# WARNING: The external packages directory includes executable binaries with +# '.exe' extensions. As such, '.exe' is an acceptable suffix, but this is +# inherently dangerous since this script can be used to download from any +# valid website. Furthermore, if an individual attempts a "man-in-the-middle" +# attack, then the user would be capable of downloading executables from +# an untrusted source. pattern = find + "[\w.-]*(\.tar\.gz|tar\.gz2|tgz|zip|exe)?" parser = MyHTMLParser(pattern) +# Creates a 'FancyURL' which allows the script to fail gracefully by catching +# HTTP error codes 30X and several 40X(where 'X' is a natural number). urlObject = FancyURLopener() obj = urlObject.open(url) parser.feed(obj.read()) +# If a file pattern was used to describe the file that should be downloaded, +# then there is the potential for multiple file matches. Currently, the script +# will detect this ambiguity and print out all the matches, while informing +# the user that he must refine his search. +# +# TODO: Prompt the user to select from a list his/her preferred target. if len(parser.targets) > 1: print "Could not resolve your download due to the number of hits." @@ -42,17 +106,27 @@ elif len(parser.targets) == 1: + print "Found: " + parser.targets[0] url += "/" + parser.targets[0] - if os.path.exists(localFile): - print "File "+ localFile +" already exists and will not be downloaded..." - elif parser.targets[0] == localFile: - urllib.urlretrieve(url, localFile) - print "Found: " + parser.targets[0] - elif parser.matcher.match(localFile) != "None": - urllib.urlretrieve(url,parser.targets[0]); - print "Found: " + parser.targets[0] + + if localFile is None: + if os.path.exists(parser.targets[0]): + print "File " + parser.targets[0] + " already exists and will not be downloaded..." + else: + urllib.urlretrieve(url, parser.targets[0]) + print "File saved as: " + parser.targets[0] else: - urllib.urlretrieve(url, parser.targets[0]); - print "WARNING: the file found \'" + parser.targets[0] + "\' does not match \'" + localFile + "\'" - print "Ensure the downloaded version is suitable." + if os.path.exists(localFile): + print "File "+ localFile +" already exists and will not be downloaded..." + else: + if parser.targets[0] == localFile: + print "File found and destination match." + elif parser.matcher.match(localFile) != "None": + print "File found matches destination pattern." + else: + print "WARNING: the file found \'" + parser.targets[0] + "\' does not match \'" + localFile + "\'" + print "Ensure the downloaded version is suitable." + + urllib.urlretrieve(url, localFile) + print "File saved as: " + localFile else: Index: /issm/branches/trunk-jpl-damage/src/c/Container/Observations.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/Container/Observations.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/Container/Observations.cpp (revision 12946) @@ -379,7 +379,7 @@ /*Solve the three linear systems*/ #if _HAVE_GSL_ - SolverxGsl(&GinvG0,Gamma,gamma0,n_obs); // Gamma^-1 gamma0 - SolverxGsl(&Ginv1, Gamma,ones,n_obs); // Gamma^-1 ones - SolverxGsl(&GinvZ, Gamma,obs,n_obs); // Gamma^-1 Z + SolverxSeq(&GinvG0,Gamma,gamma0,n_obs); // Gamma^-1 gamma0 + SolverxSeq(&Ginv1, Gamma,ones,n_obs); // Gamma^-1 ones + SolverxSeq(&GinvZ, Gamma,obs,n_obs); // Gamma^-1 Z #else _error2_("GSL is required"); Index: /issm/branches/trunk-jpl-damage/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/EnumDefinitions/EnumDefinitions.h (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/EnumDefinitions/EnumDefinitions.h (revision 12946) @@ -405,6 +405,8 @@ StressTensoryzEnum, StressTensorzzEnum, - IceVolumeEnum, - TotalSmbEnum, + IceVolumeEnum, //FIXME reposition + TotalSmbEnum, //FIXME reposition + ThicknessAlongGradientEnum, + ThicknessAcrossGradientEnum, /*}}}*/ /*Element Interpolations{{{1*/ Index: /issm/branches/trunk-jpl-damage/src/c/Makefile.am =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/Makefile.am (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/Makefile.am (revision 12946) @@ -7,5 +7,5 @@ #if SHARED -#lib_LTLIBRARIES = libISSM.la libISSMCore.la libISSMOverload.la +#lib_LTLIBRARIES = libISSMCore.la libISSMOverload.la libISSM.la #endif if PYTHON @@ -461,4 +461,8 @@ ./modules/ThicknessAbsGradientx/ThicknessAbsGradientx.cpp\ ./modules/ThicknessAbsGradientx/ThicknessAbsGradientx.h\ + ./modules/ThicknessAlongGradientx/ThicknessAlongGradientx.cpp\ + ./modules/ThicknessAlongGradientx/ThicknessAlongGradientx.h\ + ./modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.cpp\ + ./modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.h\ ./modules/RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.cpp\ ./modules/RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.h\ @@ -1004,8 +1008,9 @@ #if SHARED #libISSM_la_SOURCES = solutions/issm.cpp -#libISSM_la_LIBADD = @LTLIBOBJS@ -#libISSM_la_LDFLAGS = -prefer-pic - +#libISSM_la_LIBADD = libISSMCore.la libISSMOverload.la +# #libISSMCore_la_SOURCES = $(issm_sources) +#libISSMCore_la_LIBADD = $(PETSCLIB) $(TAOLIB) $(PLAPACKLIB) $(MUMPSLIB) $(SUPERLULIB) $(SPOOLESLIB) $(SCALAPACKLIB) $(BLACSLIB) $(HYPRELIB) $(SPAILIB) $(PROMETHEUSLIB) $(PASTIXLIB) $(MLLIB) $(DAKOTALIB) $(METISLIB) $(CHACOLIB) $(SCOTCHLIB) $(BLASLAPACKLIB) $(MKLLIB) $(MPILIB) $(MATHLIB) $(FORTRANLIB) $(GRAPHICSLIB) $(MULTITHREADINGLIB) $(OSLIBS) $(GSLLIB) $(ADOLCLIB) +# #libISSMOverload_la_SOURCES = ./shared/String/stricmp.c #endif Index: /issm/branches/trunk-jpl-damage/src/c/classes/bamg/Mesh.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/bamg/Mesh.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/bamg/Mesh.cpp (revision 12946) @@ -5644,153 +5644,4 @@ } /*}}}*/ - /*FUNCTION CloseBoundaryEdgeV2{{{*/ - AdjacentTriangle CloseBoundaryEdgeV2(I2 C,Triangle *t, double &a,double &b) { - /*Original code from Frederic Hecht (BAMG v1.01, Mesh2.cpp/CloseBoundaryEdgeV2)*/ - // walk around the vertex - // version 2 for remove the probleme if we fill the hole - //int bug=1; - // Triangle *torigine = t; - // restart: - // int dir=0; - if (t->link != 0){ - _error2_("t->link != 0"); - } - // to have a starting edges - // try the 3 edge bourna-- in case of internal hole - // and choice the best - // - // the probleme is in case of the fine and long internal hole - // for exemple neart the training edge of a wing - BamgVertex * s=0,*s1=0, *s0=0; - Icoor2 imax = MaxICoor22; - Icoor2 l0 = imax,l1 = imax; - double dd2 = imax;// infinity - AdjacentTriangle er; - int cas=-2; - for (int j=0;j<3;j++) - { - AdjacentTriangle ta=t->FindBoundaryEdge(j); - if (! (Triangle *) ta) continue; - s0 = ta.EdgeVertex(0); - s1 = ta.EdgeVertex(1); - I2 A = * s0; - I2 B = *ta.EdgeVertex(1); - I2 AB = B-A,AC=C-A,BC=B-C; - Icoor2 ACAC = (AC,AC), BCBC = (BC,BC); - Icoor2 AB2 = Norme2_2(AB); // ||AB||^2 - Icoor2 ABAC = (AB,AC); // AB.AC| - - double d2; - if ( ABAC < 0 ) // DIST A - { - if ( (d2=(double) ACAC) < dd2) - { - er = ta; - l0 = ACAC; - l1 = BCBC; - cas = 0; - s = s0; - } - } - else if (ABAC > AB2) // DIST B - { - if ( (d2=(double) BCBC) < dd2) - { - dd2 = d2; - er = Adj(ta); // other direction - l0 = BCBC; - l1 = ACAC; - cas = 1; - s = s1; - } - } - else // DIST AB - { - - double det_2 = (double) Det(AB,AC); - det_2 *= det_2; // square of area*2 of triangle ABC - d2 = det_2/ (double) AB2; // hauteur^2 in C of of triangle ABC - - if (d2 < dd2) - { - dd2 = d2; - er = ta; - l0 = (AC,AC); - l1 = (BC,BC); - s = 0; - cas = -1; - b = ((double) ABAC/(double) AB2); - a = 1 - b; - } - } - } - if (cas ==-2){ - _error2_("cas==-2"); - } - // l1 = ||C s1|| , l0 = ||C s0|| - // where s0,s1 are the vertex of the edge er - - if ( s) - { - t=er; - AdjacentTriangle edge(er); - - int kkk=0; - int linkp = t->link == 0; - - Triangle * tt=t=edge=Adj(Previous(edge)); - do { // loop over vertex s - kkk++; - if (edge.EdgeVertex(0)!=s && kkk>=10000){ - _error2_("edge.EdgeVertex(0)!=s && kkk>=10000"); - } - - int link = tt->link == 0; - if ((link + linkp) == 1) - { // a boundary edge - BamgVertex * st = edge.EdgeVertex(1); - I2 I=*st; - Icoor2 ll = Norme2_2 (C-I); - if (ll < l1) { // the other vertex is neart - s1=st; - l1=ll; - er = edge; - if(llpvector=NULL; + #endif this->svector=NULL; @@ -38,5 +40,7 @@ Vector::Vector(int M,bool fromlocalsize,int in_type){ + #ifdef _HAVE_PETSC_ pvector=NULL; + #endif svector=NULL; type=in_type; @@ -57,8 +61,22 @@ } /*}}}*/ +#ifdef _HAVE_PETSC_ +/*FUNCTION Vector::Vector(Vec petsc_vector){{{*/ +Vector::Vector(Vec petsc_vector){ + + this->type=PetscVecType; + this->svector=NULL; + this->pvector=new PetscVec(petsc_vector); + +} +/*}}}*/ +#endif /*FUNCTION Vector::Vector(IssmDouble* serial_vec,int M,int type){{{*/ Vector::Vector(IssmDouble* serial_vec,int M,int in_type){ + #ifdef _HAVE_PETSC_ pvector=NULL; + #endif + svector=NULL; type=in_type; Index: /issm/branches/trunk-jpl-damage/src/c/classes/matrix/Vector.h =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/matrix/Vector.h (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/matrix/Vector.h (revision 12946) @@ -36,4 +36,7 @@ Vector(int M,bool fromlocalsize=false,int type=PetscVecType); Vector(IssmDouble* serial_vec,int pM,int type=PetscVecType); + #ifdef _HAVE_PETSC_ + Vector(Vec petsc_vector); + #endif ~Vector(); /*}}}*/ Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/Elements/Element.h =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/Elements/Element.h (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/Elements/Element.h (revision 12946) @@ -105,4 +105,6 @@ virtual IssmDouble SurfaceAverageVelMisfit(bool process_units,int weight_index)=0; virtual IssmDouble ThicknessAbsGradient(bool process_units,int weight_index)=0; + virtual IssmDouble ThicknessAlongGradient(bool process_units,int weight_index)=0; + virtual IssmDouble ThicknessAcrossGradient(bool process_units,int weight_index)=0; virtual IssmDouble RheologyBbarAbsGradient(bool process_units,int weight_index)=0; virtual IssmDouble DragCoefficientAbsGradient(bool process_units,int weight_index)=0; Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/Elements/Penta.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/Elements/Penta.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/Elements/Penta.cpp (revision 12946) @@ -687,8 +687,12 @@ /*FUNCTION Penta::Delta18oParameterization{{{*/ void Penta::Delta18oParameterization(void){ - + /*Are we on the base? If not, return*/ + if(!IsOnBed()) return; + + int i; IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12]; IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12]; IssmDouble PrecipitationsPresentday[NUMVERTICES][12]; + IssmDouble tmp[NUMVERTICES]; IssmDouble Delta18oPresent,Delta18oLgm,Delta18oTime; IssmDouble Delta18oSurfacePresent,Delta18oSurfaceLgm,Delta18oSurfaceTime; @@ -699,6 +703,6 @@ /*Recover present day temperature and precipitation*/ - Input* input=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum); _assert_(input); - Input* input2=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum); _assert_(input2); + Input* input=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum); _assert_(input); + Input* input2=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum); _assert_(input2); Input* input3=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input3); GaussPenta* gauss=new GaussPenta(); @@ -709,15 +713,15 @@ input2->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts); input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts); - monthlyprec[iv][month]=monthlyprec[iv][month]*yts; // convertion to m/yr + PrecipitationsPresentday[iv][month]=PrecipitationsPresentday[iv][month]/yts; // converion in m/sec } } /*Recover delta18o and Delta18oSurface at present day, lgm and at time t*/ - this->parameters->FindParam(&Delta18oPresent,SurfaceforcingsDelta18oEnum,finaltime*yts); - this->parameters->FindParam(&Delta18oLgm,SurfaceforcingsDelta18oEnum,(finaltime-21000)*yts); - this->parameters->FindParam(&Delta18oTime,SurfaceforcingsDelta18oEnum,time*yts); - this->parameters->FindParam(&Delta18oSurfacePresent,SurfaceforcingsDelta18oSurfaceEnum,finaltime*yts); - this->parameters->FindParam(&Delta18oSurfaceLgm,SurfaceforcingsDelta18oSurfaceEnum,(finaltime-21000)*yts); - this->parameters->FindParam(&Delta18oSurfaceTime,SurfaceforcingsDelta18oSurfaceEnum,time*yts); + this->parameters->FindParam(&Delta18oPresent,SurfaceforcingsDelta18oEnum,finaltime); + this->parameters->FindParam(&Delta18oLgm,SurfaceforcingsDelta18oEnum,finaltime-(21000*yts)); + this->parameters->FindParam(&Delta18oTime,SurfaceforcingsDelta18oEnum,time); + this->parameters->FindParam(&Delta18oSurfacePresent,SurfaceforcingsDelta18oSurfaceEnum,finaltime); + this->parameters->FindParam(&Delta18oSurfaceLgm,SurfaceforcingsDelta18oSurfaceEnum,finaltime-(21000*yts)); + this->parameters->FindParam(&Delta18oSurfaceTime,SurfaceforcingsDelta18oSurfaceEnum,time); /*Compute the temperature and precipitation*/ @@ -734,13 +738,18 @@ TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum); for (int imonth=0;imonth<12;imonth++) { - for(int iv=0;ivAddTimeInput(newmonthinput1,imonth/12.*yts); - PentaP1Input* newmonthinput2 = new PentaP1Input(SurfaceforcingsPrecipitationEnum,&monthlyprec[iv][imonth]); - NewPrecipitationInput->AddTimeInput(newmonthinput2,imonth/12.*yts); - } - } + for(i=0;iAddTimeInput(newmonthinput1,time+imonth/12.*yts); + + for(i=0;iAddTimeInput(newmonthinput2,time+imonth/12.*yts); + } + this->inputs->AddInput(NewTemperatureInput); this->inputs->AddInput(NewPrecipitationInput); + + this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,ElementEnum); + this->InputExtrude(SurfaceforcingsPrecipitationEnum,ElementEnum); /*clean-up*/ @@ -1559,5 +1568,4 @@ } #endif - } /*}}}*/ @@ -2009,4 +2017,5 @@ name==QmuVxEnum || name==QmuVyEnum || + name==QmuVzEnum || name==QmuPressureEnum || name==QmuBedEnum || Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/Elements/Penta.h =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/Elements/Penta.h (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/Elements/Penta.h (revision 12946) @@ -162,4 +162,6 @@ IssmDouble SurfaceAverageVelMisfit(bool process_units,int weight_index); IssmDouble ThicknessAbsGradient(bool process_units,int weight_index); + IssmDouble ThicknessAlongGradient( bool process_units,int weight_index){_error2_("not supported");}; + IssmDouble ThicknessAcrossGradient(bool process_units,int weight_index){_error2_("not supported");}; void InputControlUpdate(IssmDouble scalar,bool save_parameter); #endif Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/Elements/Tria.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/Elements/Tria.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/Elements/Tria.cpp (revision 12946) @@ -915,5 +915,5 @@ IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12]; IssmDouble PrecipitationsPresentday[NUMVERTICES][12]; - IssmDouble temp[NUMVERTICES]; + IssmDouble tmp[NUMVERTICES]; IssmDouble Delta18oPresent,Delta18oLgm,Delta18oTime; IssmDouble Delta18oSurfacePresent,Delta18oSurfaceLgm,Delta18oSurfaceTime; @@ -961,12 +961,13 @@ TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum); for (int imonth=0;imonth<12;imonth++) { - for(i=0;iAddTimeInput(newmonthinput1,time+imonth/12.*yts); - for(i=0;iAddTimeInput(newmonthinput2,time+imonth/12.*yts); } + this->inputs->AddInput(NewTemperatureInput); this->inputs->AddInput(NewPrecipitationInput); @@ -2192,7 +2193,4 @@ rho_water=matpar->GetRhoFreshwater(); - //if(id==1) printf(" monthlytemperatures %f\n",monthlytemperatures[1][0]); - //if(id==1) printf("monthlyprec %f\n",monthlyprec[1][0]); - /*measure the surface mass balance*/ for (int iv = 0; ivGetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* thickness_input= inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + + /*Get all parameters at gaussian point*/ + weights_input->GetInputValue(&weight,gauss,weight_index); + thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + + /*J = 1/2 ( vx*dH/dx + vy*dH/dy )^2 */ + Jelem+=weight*1/2*(vx*dp[0] + vy*dp[1])*(vx*dp[0] + vy*dp[1])*Jdet*gauss->weight; + } + + /*Clean up and return*/ + delete gauss; + return Jelem; +} +/*}}}*/ +/*FUNCTION Tria::ThicknessAcrossGradient{{{*/ +IssmDouble Tria::ThicknessAcrossGradient(bool process_units,int weight_index){ + + /* Intermediaries */ + int ig; + IssmDouble Jelem = 0; + IssmDouble weight; + IssmDouble Jdet; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble dp[NDOF2]; + IssmDouble vx,vy; + GaussTria *gauss = NULL; + + /*retrieve parameters and inputs*/ + + /*If on water, return 0: */ + if(IsOnWater()) return 0; + + /*Retrieve all inputs we will be needing: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + Input* weights_input = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* thickness_input= inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + + /*Get all parameters at gaussian point*/ + weights_input->GetInputValue(&weight,gauss,weight_index); + thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + + /*J = 1/2 ( -vy*dH/dx + vx*dH/dy )^2 */ + Jelem+=weight*1/2*(-vy*dp[0] + vx*dp[1])*(-vy*dp[0] + vx*dp[1])*Jdet*gauss->weight; + } + + /*Clean up and return*/ + delete gauss; + return Jelem; +} +/*}}}*/ /*FUNCTION Tria::ThicknessAbsMisfit {{{*/ IssmDouble Tria::ThicknessAbsMisfit(bool process_units,int weight_index){ @@ -4453,13 +4551,14 @@ /*Intermediaries */ int i,ig,resp; - IssmDouble Jdet; - IssmDouble thickness,thicknessobs,weight; - int *responses = NULL; + IssmDouble Jdet; + IssmDouble thickness,thicknessobs,weight; int num_responses; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble basis[3]; - IssmDouble dbasis[NDOF2][NUMVERTICES]; - IssmDouble dH[2]; - GaussTria* gauss=NULL; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble basis[3]; + IssmDouble dbasis[NDOF2][NUMVERTICES]; + IssmDouble dH[2]; + IssmDouble v[2]; + GaussTria *gauss = NULL; + int *responses = NULL; /*Initialize Element vector*/ @@ -4470,7 +4569,9 @@ this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum); - Input* thickness_input = inputs->GetInput(ThicknessEnum); _assert_(thickness_input); - Input* thicknessobs_input = inputs->GetInput(InversionThicknessObsEnum);_assert_(thicknessobs_input); - Input* weights_input = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* thickness_input = inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* thicknessobs_input = inputs->GetInput(InversionThicknessObsEnum); _assert_(thicknessobs_input); + Input* weights_input = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input); /* Start looping on the number of gaussian points: */ @@ -4499,4 +4600,16 @@ for(i=0;ivalues[i]+= - weight*dH[0]*dbasis[0][i]*Jdet*gauss->weight; for(i=0;ivalues[i]+= - weight*dH[1]*dbasis[1][i]*Jdet*gauss->weight; + break; + case ThicknessAlongGradientEnum: + weights_input->GetInputValue(&weight, gauss,resp); + vx_input->GetInputValue(&v[0],gauss); + vy_input->GetInputValue(&v[1],gauss); + for(i=0;ivalues[i]+= - weight*(dH[0]*v[0]+dH[1]*v[1])*(dbasis[0][i]*v[0]+dbasis[1][i]*v[1])*Jdet*gauss->weight; + break; + case ThicknessAcrossGradientEnum: + weights_input->GetInputValue(&weight, gauss,resp); + vx_input->GetInputValue(&v[0],gauss); + vy_input->GetInputValue(&v[1],gauss); + for(i=0;ivalues[i]+= - weight*(dH[0]*(-v[1])+dH[1]*v[0])*(dbasis[0][i]*(-v[1])+dbasis[1][i]*v[0])*Jdet*gauss->weight; break; default: @@ -4675,4 +4788,10 @@ /*Nothing in P vector*/ break; + case ThicknessAlongGradientEnum: + /*Nothing in P vector*/ + break; + case ThicknessAcrossGradientEnum: + /*Nothing in P vector*/ + break; case RheologyBbarAbsGradientEnum: /*Nothing in P vector*/ @@ -4852,4 +4971,10 @@ /*Nothing in P vector*/ break; + case ThicknessAcrossGradientEnum: + /*Nothing in P vector*/ + break; + case ThicknessAlongGradientEnum: + /*Nothing in P vector*/ + break; case RheologyBbarAbsGradientEnum: /*Nothing in P vector*/ Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/Elements/Tria.h =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/Elements/Tria.h (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/Elements/Tria.h (revision 12946) @@ -160,5 +160,7 @@ IssmDouble ThicknessAbsMisfit( bool process_units,int weight_index); IssmDouble SurfaceAbsVelMisfit( bool process_units,int weight_index); - IssmDouble ThicknessAbsGradient(bool process_units,int weight_index); + IssmDouble ThicknessAbsGradient( bool process_units,int weight_index); + IssmDouble ThicknessAlongGradient( bool process_units,int weight_index); + IssmDouble ThicknessAcrossGradient(bool process_units,int weight_index); IssmDouble SurfaceRelVelMisfit( bool process_units,int weight_index); IssmDouble SurfaceLogVelMisfit( bool process_units,int weight_index); Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/Inputs/TransientInput.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/Inputs/TransientInput.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/Inputs/TransientInput.cpp (revision 12946) @@ -282,4 +282,12 @@ inputs->AddObject(input); +} +/*}}}*/ +/*FUNCTION TransientInput::Extrude{{{*/ +void TransientInput::Extrude(void){ + + for(int i=0;inumtimesteps;i++){ + ((Input*)this->inputs->GetObjectByOffset(i))->Extrude(); + } } /*}}}*/ Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/Inputs/TransientInput.h =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/Inputs/TransientInput.h (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/Inputs/TransientInput.h (revision 12946) @@ -79,5 +79,5 @@ IssmDouble Min(void); IssmDouble MinAbs(void); - void Extrude(void){_error2_("not supported yet");} + void Extrude(void); void VerticallyIntegrate(Input* thickness_forcing){_error2_("not supported yet");}; void GetVectorFromInputs(Vector* vector,int* doflist); Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KMLFileReadUtils.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KMLFileReadUtils.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KMLFileReadUtils.cpp (revision 12946) @@ -29,7 +29,7 @@ int c; int ibuf=0,buflen=1024,bufblk=1024; - char *buffer=NULL,*bufferc=NULL; - - buffer=(char *) xmalloc(buflen*sizeof(char)); + char *buffer=NULL,*bufferc=NULL,**pcom2=NULL; + + buffer=xNew(buflen); buffer[0]='\0'; @@ -57,9 +57,14 @@ if (pncom && ppcom) { (*pncom)++; - *ppcom=(char **) xrealloc(*ppcom,*pncom*sizeof(char*)); + pcom2=xNew(*pncom); + memcpy(pcom2,*ppcom,(*pncom-1)*sizeof(char*)); + xDelete(*ppcom); + *ppcom=pcom2; + pcom2=NULL; +// *ppcom=(char **) xrealloc(*ppcom,*pncom*sizeof(char*)); (*ppcom)[*pncom-1]=bufferc; } else - xfree((void**)&bufferc); + xDelete(bufferc); inew=1; continue; @@ -98,9 +103,14 @@ if (pncom && ppcom) { (*pncom)++; - *ppcom=(char **) xrealloc(*ppcom,*pncom*sizeof(char*)); + pcom2=xNew(*pncom); + memcpy(pcom2,*ppcom,(*pncom-1)*sizeof(char*)); + xDelete(*ppcom); + *ppcom=pcom2; + pcom2=NULL; +// *ppcom=(char **) xrealloc(*ppcom,*pncom*sizeof(char*)); (*ppcom)[*pncom-1]=bufferc; } else - xfree((void**)&bufferc); + xDelete(bufferc); } else { @@ -133,5 +143,5 @@ if (!ibuf) - xfree((void**)&buffer); + xDelete(buffer); return(buffer); @@ -149,5 +159,5 @@ char* buffer=NULL; - buffer=(char *) xmalloc(buflen*sizeof(char)); + buffer=xNew(buflen); buffer[0]='\0'; @@ -172,5 +182,5 @@ for (i=ibuf-1; i>=0; i--) ungetc(buffer[i],fid); - xfree((void**)&buffer); + xDelete(buffer); return(buffer); } @@ -197,5 +207,5 @@ if (!ibuf) - xfree((void**)&buffer); + xDelete(buffer); return(buffer); @@ -209,5 +219,5 @@ /* add the specified character to the token buffer */ - char* buffer=NULL; + char* buffer2=NULL; /* check buffer length and realloc if necessary */ @@ -215,5 +225,10 @@ if (*pibuf+2 > *pbuflen) { *pbuflen+=bufblk; - *pbuffer=(char *) xrealloc(*pbuffer,*pbuflen*sizeof(char)); + buffer2=xNew(*pbuflen); + memcpy(buffer2,*pbuffer,(*pbuflen-bufblk)*sizeof(char)); + xDelete(*pbuffer); + *pbuffer=buffer2; + buffer2=NULL; +// *pbuffer=(char *) xrealloc(*pbuffer,*pbuflen*sizeof(char)); } @@ -248,5 +263,5 @@ /* strtok modifies ktag, so work on copy */ - ktagi=(char *) xmalloc((strlen(ktag)+1)*sizeof(char)); + ktagi=xNew(strlen(ktag)+1); memcpy(ktagi,ktag,(strlen(ktag)+1)*sizeof(char)); @@ -258,7 +273,7 @@ if (!pname) { if (maxlen) - pname=(char *) xmalloc((maxlen +1)*sizeof(char)); - else - pname=(char *) xmalloc((strlen(ktokn)+1)*sizeof(char)); + pname=xNew(maxlen +1); + else + pname=xNew(strlen(ktokn)+1); } @@ -271,5 +286,5 @@ memcpy(pname,ktokn,(strlen(ktokn)+1)*sizeof(char)); - xfree((void**)&ktagi); + xDelete(ktagi); if (m) @@ -293,5 +308,5 @@ /* strtok modifies ktag, so work on copy */ - ktagi=(char *) xmalloc((strlen(ktag)+1)*sizeof(char)); + ktagi=xNew(strlen(ktag)+1); memcpy(ktagi,ktag,(strlen(ktag)+1)*sizeof(char)); @@ -315,5 +330,5 @@ } - xfree((void**)&ktagi); + xDelete(ktagi); /* check for xml declaration, dtd declaration, or solo tag */ @@ -343,5 +358,5 @@ sscanf(kstr,"%d",pival); - xfree((void**)&kstr); + xDelete(kstr); /* get additional token and compare to closing tag */ @@ -355,5 +370,5 @@ {_error2_("KMLFileTokenParse -- Missing closing tag for " << ktag << ".\n");} else - xfree((void**)&kstr); + xDelete(kstr); // _pprintLine_("KMLFileTokenParse -- " << ktag << "=" << *pival << "."); @@ -377,5 +392,5 @@ sscanf(kstr,"%d",&ival); *pbval=(bool)ival; - xfree((void**)&kstr); + xDelete(kstr); /* get additional token and compare to closing tag */ @@ -389,5 +404,5 @@ {_error2_("KMLFileTokenParse -- Missing closing tag for " << ktag << ".\n");} else - xfree((void**)&kstr); + xDelete(kstr); // _pprintLine_("KMLFileTokenParse -- " << ktag << "=" << (*pbval ? "true" : "false") << "."); @@ -422,7 +437,7 @@ if (!pstr) { if (maxlen) - pstr=(char *) xmalloc((maxlen +1)*sizeof(char)); - else - pstr=(char *) xmalloc((strlen(kstr)+1)*sizeof(char)); + pstr=xNew(maxlen +1); + else + pstr=xNew(strlen(kstr)+1); } @@ -435,5 +450,5 @@ memcpy(pstr,kstr,(strlen(kstr)+1)*sizeof(char)); - xfree((void**)&kstr); + xDelete(kstr); if (m) @@ -450,5 +465,5 @@ {_error2_("KMLFileTokenParse -- Missing closing tag for " << ktag << ".\n");} else - xfree((void**)&kstr); + xDelete(kstr); // _pprintLine_("KMLFileTokenParse -- " << ktag << "=\"" << pstr << "\"."); @@ -472,5 +487,5 @@ sscanf(kstr,"%g",pfval); - xfree((void**)&kstr); + xDelete(kstr); /* get additional token and compare to closing tag */ @@ -484,5 +499,5 @@ {_error2_("KMLFileTokenParse -- Missing closing tag for " << ktag << ".\n");} else - xfree((void**)&kstr); + xDelete(kstr); // _pprintLine_("KMLFileTokenParse -- " << ktag << "=" << *pfval << "."); @@ -506,5 +521,5 @@ sscanf(kstr,"%lg",pdval); - xfree((void**)&kstr); + xDelete(kstr); /* get additional token and compare to closing tag */ @@ -518,5 +533,5 @@ {_error2_("KMLFileTokenParse -- Missing closing tag for " << ktag << ".\n");} else - xfree((void**)&kstr); + xDelete(kstr); // _pprintLine_("KMLFileTokenParse -- " << ktag << "=" << *pdval << "."); @@ -533,4 +548,5 @@ char* kstr; char* ktok; + double* dval2=NULL; char delim[]={' ',',','\f','\n','\r','\t','\v','\0'}; @@ -544,7 +560,7 @@ if (!*pdval) if (maxlen) - *pdval=(double *) xmalloc(maxlen *sizeof(double)); - else - *pdval=(double *) xmalloc(((strlen(kstr)+1)/2)*sizeof(double)); + *pdval=xNew(maxlen ); + else + *pdval=xNew((strlen(kstr)+1)/2); /* loop through string to get all values */ @@ -558,8 +574,13 @@ ktok=strtok(NULL,delim); } - xfree((void**)&kstr); + xDelete(kstr); if (!maxlen) - *pdval=(double *) xrealloc(*pdval,(i+1)*sizeof(double)); + dval2=xNew(i+1); + memcpy(dval2,*pdval,(i+1)*sizeof(double)); + xDelete(*pdval); + *pdval=dval2; + dval2=NULL; +// *pdval=(double *) xrealloc(*pdval,(i+1)*sizeof(double)); if (m) @@ -576,5 +597,5 @@ {_error2_("KMLFileTokenParse -- Missing closing tag for " << ktag << ".\n");} else - xfree((void**)&kstr); + xDelete(kstr); // _pprintLine_("KMLFileTokenParse -- " << ktag << "=..."); @@ -586,11 +607,12 @@ /*}}}*/ /*FUNCTION KMLFileTokenParse {{{*/ -int KMLFileTokenParse(double (**pdval3)[3],int* m,int maxlen, +int KMLFileTokenParse(double **pdval,int* m,int n,int maxlen, char* ktag, FILE* fid){ - int i=0,j=-1; + int i=-1,j=-1; char* kstr; char* ktok; + double* dval2=NULL; char delim[]={' ',',','\f','\n','\r','\t','\v','\0'}; @@ -600,11 +622,11 @@ NULL,NULL)) || (kstr[0] == '<')) - _error2_("KMLFileTokenParse -- Missing double [m x 3] field for " << ktag << ".\n"); - - if (!*pdval3) + _error2_("KMLFileTokenParse -- Missing double [m x n] field for " << ktag << ".\n"); + + if (!*pdval) if (maxlen) - *pdval3=(double (*)[3]) xmalloc((maxlen*3) *sizeof(double)); - else - *pdval3=(double (*)[3]) xmalloc(((strlen(kstr)+1)/2)*sizeof(double)); + *pdval=xNew(maxlen*n ); + else + *pdval=xNew((strlen(kstr)+1)/2); /* loop through string to get all values */ @@ -612,24 +634,26 @@ ktok=strtok(kstr,delim); while (ktok) { - j++; - if (j == 3) { - i++; - j=0; - if (maxlen && (maxlen < i+1)) - _error2_("KMLFileTokenParse -- Double [m x 3] field too short for " << ktag << ".\n"); - } - sscanf(ktok,"%lg",&((*pdval3)[i][j])); + i++; + if (maxlen && (maxlen*n < i+1)) + _error2_("KMLFileTokenParse -- Double [m x n] field too short for " << ktag << ".\n"); + j=(j+1) % n; + sscanf(ktok,"%lg",&((*pdval)[i])); ktok=strtok(NULL,delim); } - xfree((void**)&kstr); + xDelete(kstr); if (!maxlen) - *pdval3=(double (*)[3]) xrealloc(*pdval3,((i+1)*3)*sizeof(double)); + dval2=xNew((i+1)*n); + memcpy(dval2,*pdval,((i+1)*n)*sizeof(double)); + xDelete(*pdval); + *pdval=dval2; + dval2=NULL; +// *pdval=(double *) xrealloc(*pdval,((i+1)*n)*sizeof(double)); if (m) - *m=i+1; - - if (j != 2) - _pprintLine_("KMLFileTokenParse -- Double [m x 3] field for " << ktag << " does not have multiple of 3 values."); + *m=((i+1)+(n-1))/n; + + if ((i+1) % n) + _pprintLine_("KMLFileTokenParse -- Double [m x n] field for " << ktag << " does not have multiple of n values."); /* get additional token and compare to closing tag */ @@ -643,9 +667,9 @@ {_error2_("KMLFileTokenParse -- Missing closing tag for " << ktag << ".\n");} else - xfree((void**)&kstr); + xDelete(kstr); // _pprintLine_("KMLFileTokenParse -- " << ktag << "=..."); // for (j=0; j<=i; j++) -// _pprintLine_(" [" << j << "][0-2]: " << (*pdval3)[j][0] << "g," << (*pdval3)[j][1] << "g," << (*pdval3)[j][2] << "g"); +// _pprintLine_(" [" << j << "]: " << (*pdval)[j] << "g"); return(0); @@ -671,5 +695,5 @@ (!strncmp(&(kstr[2]),&(ktag[1]),(strcspn(ktag," >")-1)/sizeof(char)))) { _pprintLine_("KMLFileTagSkip -- closing tag " << kstr << "."); - xfree((void**)&kstr); + xDelete(kstr); return(0); } @@ -691,5 +715,5 @@ } - xfree((void**)&kstr); + xDelete(kstr); } Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KMLFileReadUtils.h =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KMLFileReadUtils.h (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KMLFileReadUtils.h (revision 12946) @@ -48,5 +48,5 @@ char* ktag, FILE* fid); -int KMLFileTokenParse(double (**pdval3)[3],int* m,int maxlen, +int KMLFileTokenParse(double **pdval,int* m,int n,int maxlen, char* ktag, FILE* fid); Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Attribute.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Attribute.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Attribute.cpp (revision 12946) @@ -32,6 +32,6 @@ KML_Attribute::~KML_Attribute(){ - if (name ) xfree((void**)&name); - if (value ) xfree((void**)&value); + if (name ) xDelete(name); + if (value ) xDelete(value); } @@ -100,8 +100,8 @@ void KML_Attribute::Alloc(const char* namei,const char* valuei){ - name =(char *) xmalloc((strlen(namei )+1)*sizeof(char)); + name =xNew(strlen(namei )+1); memcpy(name,namei,(strlen(namei)+1)*sizeof(char)); - value=(char *) xmalloc((strlen(valuei)+1)*sizeof(char)); + value=xNew(strlen(valuei)+1); memcpy(value,valuei,(strlen(valuei)+1)*sizeof(char)); @@ -121,9 +121,9 @@ if (!value || !strlen(value)) { - *pvalueo=(char *) xmalloc((strlen(deflt)+1)*sizeof(char)); + *pvalueo=xNew(strlen(deflt)+1); memcpy(*pvalueo,deflt,(strlen(deflt)+1)*sizeof(char)); } else { - *pvalueo=(char *) xmalloc((strlen(value)+1)*sizeof(char)); + *pvalueo=xNew(strlen(value)+1); memcpy(*pvalueo,value,(strlen(value)+1)*sizeof(char)); } Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Comment.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Comment.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Comment.cpp (revision 12946) @@ -31,5 +31,5 @@ KML_Comment::~KML_Comment(){ - if (value ) xfree((void**)&value); + if (value ) xDelete(value); } @@ -94,5 +94,5 @@ void KML_Comment::Alloc(const char* valuei){ - value=(char *) xmalloc((strlen(valuei)+1)*sizeof(char)); + value=xNew(strlen(valuei)+1); memcpy(value,valuei,(strlen(valuei)+1)*sizeof(char)); @@ -111,5 +111,5 @@ void KML_Comment::Get(char** pvalueo){ - *pvalueo=(char *) xmalloc((strlen(value)+1)*sizeof(char)); + *pvalueo=xNew(strlen(value)+1); memcpy(*pvalueo,value,(strlen(value)+1)*sizeof(char)); Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Container.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Container.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Container.cpp (revision 12946) @@ -116,5 +116,5 @@ if (!strncmp(kstr,"(kstr); return; } Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Document.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Document.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Document.cpp (revision 12946) @@ -102,5 +102,5 @@ &ncom,&pcom)) { if (!strncmp(kstri,"(kstri); break; } @@ -113,5 +113,5 @@ KML_Container::Read(fid,kstri); - xfree((void**)&kstri); + xDelete(kstri); } @@ -119,6 +119,6 @@ for (ncom; ncom>0; ncom--) - xfree((void**)&(pcom[ncom-1])); - xfree((void**)&pcom); + xDelete(pcom[ncom-1]); + xDelete(pcom); return; Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_File.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_File.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_File.cpp (revision 12946) @@ -103,5 +103,5 @@ &ncom,&pcom)) { if (!strncmp(kstri,"(kstri); break; } @@ -114,5 +114,5 @@ KML_Object::Read(fid,kstri); - xfree((void**)&kstri); + xDelete(kstri); } @@ -120,6 +120,6 @@ for (ncom; ncom>0; ncom--) - xfree((void**)&(pcom[ncom-1])); - xfree((void**)&pcom); + xDelete(pcom[ncom-1]); + xDelete(pcom); return; Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Folder.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Folder.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Folder.cpp (revision 12946) @@ -102,5 +102,5 @@ &ncom,&pcom)) { if (!strncmp(kstri,"(kstri); break; } @@ -113,5 +113,5 @@ KML_Container::Read(fid,kstri); - xfree((void**)&kstri); + xDelete(kstri); } @@ -119,6 +119,6 @@ for (ncom; ncom>0; ncom--) - xfree((void**)&(pcom[ncom-1])); - xfree((void**)&pcom); + xDelete(pcom[ncom-1]); + xDelete(pcom); return; Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_GroundOverlay.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_GroundOverlay.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_GroundOverlay.cpp (revision 12946) @@ -127,5 +127,5 @@ &ncom,&pcom)) { if (!strncmp(kstri,"(kstri); break; } @@ -151,5 +151,5 @@ KML_Overlay::Read(fid,kstri); - xfree((void**)&kstri); + xDelete(kstri); } @@ -157,6 +157,6 @@ for (ncom; ncom>0; ncom--) - xfree((void**)&(pcom[ncom-1])); - xfree((void**)&pcom); + xDelete(pcom[ncom-1]); + xDelete(pcom); return; Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Icon.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Icon.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Icon.cpp (revision 12946) @@ -141,5 +141,5 @@ &ncom,&pcom)) { if (!strncmp(kstri,"(kstri); break; } @@ -169,5 +169,5 @@ KML_Object::Read(fid,kstri); - xfree((void**)&kstri); + xDelete(kstri); } @@ -175,6 +175,6 @@ for (ncom; ncom>0; ncom--) - xfree((void**)&(pcom[ncom-1])); - xfree((void**)&pcom); + xDelete(pcom[ncom-1]); + xDelete(pcom); return; Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_LatLonBox.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_LatLonBox.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_LatLonBox.cpp (revision 12946) @@ -117,5 +117,5 @@ &ncom,&pcom)) { if (!strncmp(kstri,"(kstri); break; } @@ -149,5 +149,5 @@ KML_Object::Read(fid,kstri); - xfree((void**)&kstri); + xDelete(kstri); } @@ -155,6 +155,6 @@ for (ncom; ncom>0; ncom--) - xfree((void**)&(pcom[ncom-1])); - xfree((void**)&pcom); + xDelete(pcom[ncom-1]); + xDelete(pcom); return; Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_LineString.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_LineString.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_LineString.cpp (revision 12946) @@ -137,5 +137,5 @@ &ncom,&pcom)) { if (!strncmp(kstri,"(kstri); break; } @@ -165,5 +165,5 @@ KML_Geometry::Read(fid,kstri); - xfree((void**)&kstri); + xDelete(kstri); } @@ -171,6 +171,6 @@ for (ncom; ncom>0; ncom--) - xfree((void**)&(pcom[ncom-1])); - xfree((void**)&pcom); + xDelete(pcom[ncom-1]); + xDelete(pcom); return; @@ -186,6 +186,6 @@ /* extract latitude and longitude into vectors */ - lat=(double *) xmalloc(ncoord*sizeof(double)); - lon=(double *) xmalloc(ncoord*sizeof(double)); + lat=xNew(ncoord); + lon=xNew(ncoord); for (i=0; i(ncoord); + y =xNew(ncoord); Ll2xyx(x,y,lat,lon,ncoord,sgn,cm,sp); @@ -221,10 +221,10 @@ fprintf(fid,"\n"); - xfree((void**)&y); - xfree((void**)&x); - xfree((void**)&lon); - xfree((void**)&lat); - - return; -} -/*}}}*/ + xDelete(y); + xDelete(x); + xDelete(lon); + xDelete(lat); + + return; +} +/*}}}*/ Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_LineStyle.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_LineStyle.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_LineStyle.cpp (revision 12946) @@ -109,5 +109,5 @@ &ncom,&pcom)) { if (!strncmp(kstri,"(kstri); break; } @@ -125,5 +125,5 @@ KML_ColorStyle::Read(fid,kstri); - xfree((void**)&kstri); + xDelete(kstri); } @@ -131,6 +131,6 @@ for (ncom; ncom>0; ncom--) - xfree((void**)&(pcom[ncom-1])); - xfree((void**)&pcom); + xDelete(pcom[ncom-1]); + xDelete(pcom); return; Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_LinearRing.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_LinearRing.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_LinearRing.cpp (revision 12946) @@ -150,5 +150,5 @@ KMLFileTokenParse(altmode,NULL,KML_LINEARRING_ALTMODE_LENGTH,kstri,fid); else if (!strcmp(kstri,"")) - KMLFileTokenParse(&coords,&ncoord,0,kstri,fid); + KMLFileTokenParse(&coords,&ncoord,3,0,kstri,fid); else if (!strncmp(kstri,"<",1)) KML_Geometry::Read(fid,kstri); @@ -159,5 +159,6 @@ this->AddCommnt(ncom,pcom); - for(ncom;ncom>0;ncom--) xDelete((pcom[ncom-1])); + for(ncom; ncom>0; ncom--) + xDelete(pcom[ncom-1]); xDelete(pcom); @@ -174,6 +175,6 @@ /* extract latitude and longitude into vectors */ - lat=xNew(ncoord); - lon=xNew(ncoord); + lat=xNew(ncoord); + lon=xNew(ncoord); for (i=0; i(ncoord); - y =xNew(ncoord); + x =xNew(ncoord); + y =xNew(ncoord); Ll2xyx(x,y,lat,lon,ncoord,sgn,cm,sp); Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_MultiGeometry.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_MultiGeometry.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_MultiGeometry.cpp (revision 12946) @@ -136,5 +136,5 @@ &ncom,&pcom)) { if (!strncmp(kstri,"(kstri); break; } @@ -177,5 +177,5 @@ KML_Geometry::Read(fid,kstri); - xfree((void**)&kstri); + xDelete(kstri); } @@ -183,6 +183,6 @@ for (ncom; ncom>0; ncom--) - xfree((void**)&(pcom[ncom-1])); - xfree((void**)&pcom); + xDelete(pcom[ncom-1]); + xDelete(pcom); return; Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Object.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Object.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Object.cpp (revision 12946) @@ -304,5 +304,5 @@ katt->Get(pvalue,deflt); else { - *pvalue=(char *) xmalloc((strlen(deflt)+1)*sizeof(char)); + *pvalue=xNew(strlen(deflt)+1); memcpy(*pvalue,deflt,(strlen(deflt)+1)*sizeof(char)); } Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Overlay.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Overlay.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Overlay.cpp (revision 12946) @@ -106,5 +106,5 @@ if (!strncmp(kstr,"(kstr); return; } Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Placemark.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Placemark.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Placemark.cpp (revision 12946) @@ -136,5 +136,5 @@ &ncom,&pcom)) { if (!strncmp(kstri,"(kstri); break; } @@ -177,5 +177,5 @@ KML_Feature::Read(fid,kstri); - xfree((void**)&kstri); + xDelete(kstri); } @@ -183,6 +183,6 @@ for (ncom; ncom>0; ncom--) - xfree((void**)&(pcom[ncom-1])); - xfree((void**)&pcom); + xDelete(pcom[ncom-1]); + xDelete(pcom); return; Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Point.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Point.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Point.cpp (revision 12946) @@ -122,5 +122,5 @@ &ncom,&pcom)) { if (!strncmp(kstri,"(kstri); break; } @@ -140,5 +140,5 @@ KML_Geometry::Read(fid,kstri); - xfree((void**)&kstri); + xDelete(kstri); } @@ -146,6 +146,6 @@ for (ncom; ncom>0; ncom--) - xfree((void**)&(pcom[ncom-1])); - xfree((void**)&pcom); + xDelete(pcom[ncom-1]); + xDelete(pcom); return; Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_PolyStyle.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_PolyStyle.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_PolyStyle.cpp (revision 12946) @@ -113,5 +113,5 @@ &ncom,&pcom)) { if (!strncmp(kstri,"(kstri); break; } @@ -133,5 +133,5 @@ KML_ColorStyle::Read(fid,kstri); - xfree((void**)&kstri); + xDelete(kstri); } @@ -139,6 +139,6 @@ for (ncom; ncom>0; ncom--) - xfree((void**)&(pcom[ncom-1])); - xfree((void**)&pcom); + xDelete(pcom[ncom-1]); + xDelete(pcom); return; Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Polygon.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Polygon.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Polygon.cpp (revision 12946) @@ -175,5 +175,5 @@ &ncom,&pcom)) { if (!strncmp(kstri,"(kstri); break; } @@ -203,5 +203,5 @@ &ncom,&pcom)) { if (!strncmp(kstrj,"(kstrj); break; } @@ -220,5 +220,5 @@ KML_Geometry::Read(fid,kstrj); - xfree((void**)&kstrj); + xDelete(kstrj); } @@ -230,5 +230,5 @@ &ncom,&pcom)) { if (!strncmp(kstrj,"(kstrj); break; } @@ -247,5 +247,5 @@ KML_Geometry::Read(fid,kstrj); - xfree((void**)&kstrj); + xDelete(kstrj); } @@ -254,5 +254,5 @@ KML_Geometry::Read(fid,kstri); - xfree((void**)&kstri); + xDelete(kstri); } @@ -260,6 +260,6 @@ for (ncom; ncom>0; ncom--) - xfree((void**)&(pcom[ncom-1])); - xfree((void**)&pcom); + xDelete(pcom[ncom-1]); + xDelete(pcom); return; Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Style.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Style.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Style.cpp (revision 12946) @@ -38,10 +38,8 @@ if (list) { // delete list; - if (list ) xfree((void**)&list); list =NULL; } if (balloon) { // delete balloon; - if (balloon ) xfree((void**)&balloon); balloon =NULL; } @@ -56,10 +54,8 @@ if (label) { // delete label; - if (label ) xfree((void**)&label); label =NULL; } if (icon) { // delete icon; - if (icon ) xfree((void**)&icon); icon =NULL; } @@ -189,5 +185,5 @@ &ncom,&pcom)) { if (!strncmp(kstri,"(kstri); break; } @@ -230,5 +226,5 @@ KML_StyleSelector::Read(fid,kstri); - xfree((void**)&kstri); + xDelete(kstri); } @@ -236,8 +232,8 @@ for (ncom; ncom>0; ncom--) - xfree((void**)&(pcom[ncom-1])); - xfree((void**)&pcom); - - return; -} -/*}}}*/ + xDelete(pcom[ncom-1]); + xDelete(pcom); + + return; +} +/*}}}*/ Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Unknown.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Unknown.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/KML/KML_Unknown.cpp (revision 12946) @@ -32,6 +32,6 @@ KML_Unknown::~KML_Unknown(){ - if (name ) xfree((void**)&name); - if (value ) xfree((void**)&value); + if (name ) xDelete(name); + if (value ) xDelete(value); } @@ -77,5 +77,5 @@ if (value ) { - valuei=(char *) xmalloc((strlen(value)+1)*sizeof(char)); + valuei=xNew(strlen(value)+1); memcpy(valuei,value,(strlen(value)+1)*sizeof(char)); @@ -87,5 +87,5 @@ if(flag) _pprintLine_("\""); - xfree((void**)&valuei); + xDelete(valuei); } else @@ -108,5 +108,5 @@ if (value ) { - valuei=(char *) xmalloc((strlen(value)+1)*sizeof(char)); + valuei=xNew(strlen(value)+1); memcpy(valuei,value,(strlen(value)+1)*sizeof(char)); @@ -117,5 +117,5 @@ fprintf(filout,"%s %s\n",indent,vtoken); - xfree((void**)&valuei); + xDelete(valuei); } @@ -131,4 +131,5 @@ char* kstri; + char* value2=NULL; int ncom=0; char** pcom=NULL; @@ -155,5 +156,5 @@ !strncmp(&kstri[2],name,strlen(name))) { // _pprintLine_("KML_Unknown::Read -- closing name=" << name << "."); - xfree((void**)&kstri); + xDelete(kstri); break; } @@ -163,10 +164,15 @@ else if (strncmp(kstri,"<",1)) { if (value) { - value=(char *) xrealloc(value,(strlen(value)+1+strlen(kstri)+1)*sizeof(char)); + value2=xNew(strlen(value)+1+strlen(kstri)+1); + memcpy(value2,value,(strlen(value)+1)*sizeof(char)); + xDelete(value); + value=value2; + value2=NULL; +// value=(char *) xrealloc(value,(strlen(value)+1+strlen(kstri)+1)*sizeof(char)); strcat(value,nl); strcat(value,kstri); } else { - value=(char *) xmalloc((strlen(kstri)+1)*sizeof(char)); + value=xNew(strlen(kstri)+1); memcpy(value,kstri,(strlen(kstri)+1)*sizeof(char)); } @@ -176,5 +182,5 @@ KML_Object::Read(fid,kstri); - xfree((void**)&kstri); + xDelete(kstri); } @@ -182,6 +188,6 @@ for (ncom; ncom>0; ncom--) - xfree((void**)&(pcom[ncom-1])); - xfree((void**)&pcom); + xDelete(pcom[ncom-1]); + xDelete(pcom); return; Index: /issm/branches/trunk-jpl-damage/src/c/modules/EnumToStringx/EnumToStringx.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 12946) @@ -398,4 +398,6 @@ case IceVolumeEnum : return "IceVolume"; case TotalSmbEnum : return "TotalSmb"; + case ThicknessAlongGradientEnum : return "ThicknessAlongGradient"; + case ThicknessAcrossGradientEnum : return "ThicknessAcrossGradient"; case P0Enum : return "P0"; case P1Enum : return "P1"; Index: /issm/branches/trunk-jpl-damage/src/c/modules/Responsex/Responsex.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/modules/Responsex/Responsex.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/modules/Responsex/Responsex.cpp (revision 12946) @@ -41,4 +41,6 @@ case ThicknessAbsMisfitEnum: ThicknessAbsMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break; case ThicknessAbsGradientEnum: ThicknessAbsGradientx( responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break; + case ThicknessAlongGradientEnum: ThicknessAlongGradientx( responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break; + case ThicknessAcrossGradientEnum: ThicknessAcrossGradientx( responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break; case TotalSmbEnum: TotalSmbx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; case RheologyBbarAbsGradientEnum:RheologyBbarAbsGradientx( responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break; Index: /issm/branches/trunk-jpl-damage/src/c/modules/Solverx/Solverx.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/modules/Solverx/Solverx.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/modules/Solverx/Solverx.cpp (revision 12946) @@ -30,16 +30,20 @@ /*According to matrix type, use specific solvers: */ - if(Kff->type==PetscMatType){ - PetscVec* uf0_vector = NULL; - PetscVec* df_vector = NULL; - if(uf0) uf0_vector = uf0->pvector; - if(df) df_vector = df->pvector; - - SolverxPetsc(&uf->pvector,Kff->pmatrix,pf->pvector,uf0_vector,df_vector,parameters); + switch(Kff->type){ + #ifdef _HAVE_PETSC_ + case PetscMatType:{ + PetscVec* uf0_vector = NULL; + PetscVec* df_vector = NULL; + if(uf0) uf0_vector = uf0->pvector; + if(df) df_vector = df->pvector; + SolverxPetsc(&uf->pvector,Kff->pmatrix,pf->pvector,uf0_vector,df_vector,parameters); + break;} + #endif + case SeqMatType:{ + SolverxSeq(&uf->svector,Kff->smatrix,pf->svector); + break;} + default: + _error2_("Matrix type: " << Kff->type << " not supported yet!"); } - else if(Kff->type==SeqMatType){ - SolverxSeq(&uf->svector,Kff->smatrix,pf->svector); - } - else _error2_("Matrix type: " << Kff->type << " not supported yet!"); /*Assign output pointers:*/ Index: /issm/branches/trunk-jpl-damage/src/c/modules/StringToEnumx/StringToEnumx.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 12946) @@ -408,4 +408,6 @@ else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum; else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum; + else if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum; + else if (strcmp(name,"ThicknessAcrossGradient")==0) return ThicknessAcrossGradientEnum; else if (strcmp(name,"P0")==0) return P0Enum; else if (strcmp(name,"P1")==0) return P1Enum; Index: /issm/branches/trunk-jpl-damage/src/c/modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.cpp (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/c/modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.cpp (revision 12946) @@ -0,0 +1,37 @@ +/*!\file ThicknessAcrossGradientx + * \brief: compute misfit between observations and model + */ + +#include "./ThicknessAcrossGradientx.h" + +#include "../../shared/shared.h" +#include "../../include/include.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" + +void ThicknessAcrossGradientx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units, int weight_index){ + + /*Intermediary*/ + int i; + Element* element=NULL; + + /*output: */ + double J=0; + double J_sum; + + /*Compute Misfit: */ + for (i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + J+=element->ThicknessAcrossGradient(process_units,weight_index); + } + + /*Sum all J from all cpus of the cluster:*/ + #ifdef _HAVE_MPI_ + MPI_Reduce (&J,&J_sum,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&J_sum,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + J=J_sum; + #endif + + /*Assign output pointers: */ + *pJ=J; +} Index: /issm/branches/trunk-jpl-damage/src/c/modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.h =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.h (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/c/modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.h (revision 12946) @@ -0,0 +1,14 @@ +/*!\file: ThicknessAcrossGradientx.h + * \brief header file for inverse methods misfit computation + */ + +#ifndef _THICKNESSACROSSGRADIENT_H +#define _THICKNESSACROSSGRADIENT_H + +#include "../../Container/Container.h" +#include "../../classes/objects/objects.h" + +/* local prototypes: */ +void ThicknessAcrossGradientx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units,int weight_index); + +#endif Index: /issm/branches/trunk-jpl-damage/src/c/modules/ThicknessAlongGradientx/ThicknessAlongGradientx.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/modules/ThicknessAlongGradientx/ThicknessAlongGradientx.cpp (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/c/modules/ThicknessAlongGradientx/ThicknessAlongGradientx.cpp (revision 12946) @@ -0,0 +1,37 @@ +/*!\file ThicknessAlongGradientx + * \brief: compute misfit between observations and model + */ + +#include "./ThicknessAlongGradientx.h" + +#include "../../shared/shared.h" +#include "../../include/include.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" + +void ThicknessAlongGradientx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units, int weight_index){ + + /*Intermediary*/ + int i; + Element* element=NULL; + + /*output: */ + double J=0; + double J_sum; + + /*Compute Misfit: */ + for (i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + J+=element->ThicknessAlongGradient(process_units,weight_index); + } + + /*Sum all J from all cpus of the cluster:*/ + #ifdef _HAVE_MPI_ + MPI_Reduce (&J,&J_sum,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&J_sum,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + J=J_sum; + #endif + + /*Assign output pointers: */ + *pJ=J; +} Index: /issm/branches/trunk-jpl-damage/src/c/modules/ThicknessAlongGradientx/ThicknessAlongGradientx.h =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/modules/ThicknessAlongGradientx/ThicknessAlongGradientx.h (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/c/modules/ThicknessAlongGradientx/ThicknessAlongGradientx.h (revision 12946) @@ -0,0 +1,14 @@ +/*!\file: ThicknessAlongGradientx.h + * \brief header file for inverse methods misfit computation + */ + +#ifndef _THICKNESSALONGGRADIENT_H +#define _THICKNESSALONGGRADIENT_H + +#include "../../Container/Container.h" +#include "../../classes/objects/objects.h" + +/* local prototypes: */ +void ThicknessAlongGradientx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units,int weight_index); + +#endif Index: /issm/branches/trunk-jpl-damage/src/c/modules/modules.h =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/modules/modules.h (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/modules/modules.h (revision 12946) @@ -119,4 +119,6 @@ #include "./ThicknessAbsMisfitx/ThicknessAbsMisfitx.h" #include "./ThicknessAbsGradientx/ThicknessAbsGradientx.h" +#include "./ThicknessAlongGradientx/ThicknessAlongGradientx.h" +#include "./ThicknessAcrossGradientx/ThicknessAcrossGradientx.h" #include "./UpdateVertexPositionsx/UpdateVertexPositionsx.h" #include "./UpdateConstraintsx/UpdateConstraintsx.h" Index: /issm/branches/trunk-jpl-damage/src/c/shared/Elements/ComputeDelta18oTemperaturePrecipitation.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/shared/Elements/ComputeDelta18oTemperaturePrecipitation.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/shared/Elements/ComputeDelta18oTemperaturePrecipitation.cpp (revision 12946) @@ -16,16 +16,10 @@ IssmDouble glacialindex; // used to vary present day temperature int imonth=0; - // printf("Delta18oTime %f\n",Delta18oTime); - // printf("Delta18oPresent %f\n",Delta18oPresent); - // printf("Delta18oLgm %f\n",Delta18oLgm); - glacialindex = (Delta18oTime-Delta18oPresent-delta18oLapseRate*(Delta18oSurfaceTime-Delta18oSurfacePresent)) - /(Delta18oLgm-Delta18oPresent-delta18oLapseRate*(Delta18oSurfaceLgm-Delta18oSurfacePresent)); - // printf(" monthlytemperatureseeeeee %f\n",TemperaturesPresentday[1]); - // printf(" monthlypreceeeeeeee %f\n",TemperaturesLgm[1]); - // printf(" glacialindex %f\n",glacialindex); + glacialindex = 0;//(Delta18oTime-Delta18oPresent-delta18oLapseRate*(Delta18oSurfaceTime-Delta18oSurfacePresent)) + // /(Delta18oLgm-Delta18oPresent-delta18oLapseRate*(Delta18oSurfaceLgm-Delta18oSurfacePresent)); for (int imonth = 0; imonth<12; imonth++){ - monthlytemperaturestmp[imonth] = TemperaturesPresentday[imonth];//glacialindex*TemperaturesLgm[imonth] + (1-glacialindex)*TemperaturesPresentday[imonth]; + monthlytemperaturestmp[imonth] = glacialindex*TemperaturesLgm[imonth] + (1-glacialindex)*TemperaturesPresentday[imonth]; //monthlyprectmp[imonth] = 1.5*pow(2,((monthlytemperaturestmp[imonth]-273.15-0)/10)); //equation from rob's paper monthlyprectmp[imonth] = PrecipitationsPresentday[imonth]; Index: /issm/branches/trunk-jpl-damage/src/c/shared/Numerics/UnitConversion.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/shared/Numerics/UnitConversion.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/shared/Numerics/UnitConversion.cpp (revision 12946) @@ -51,25 +51,25 @@ case TimesteppingTimeStepEnum: scale=1.0/yts;break; //yr case TimeEnum: scale=1.0/yts;break; //yr - case HydrologyWaterVxEnum: scale=yts;break; //m/yr - case HydrologyWaterVyEnum: scale=yts;break; //m/yr - case VxEnum: scale=yts;break; //m/yr - case InversionVxObsEnum: scale=yts;break; //m/yr - case VyEnum: scale=yts;break; //m/yr - case InversionVyObsEnum: scale=yts;break; //m/yr - case VzEnum: scale=yts;break; //m/yr - case InversionVzObsEnum: scale=yts;break; //m/yr - case VelEnum: scale=yts;break; //m/yr - case InversionVelObsEnum: scale=yts;break; //m/yr - case BalancethicknessThickeningRateEnum: scale=yts;break; //m/yr - case BasalforcingsMeltingRateEnum: scale=yts;break; //m/yr - case BasalforcingsMeltingRateCorrectionEnum: scale=yts;break; //m/yr - case SurfaceforcingsPrecipitationEnum: scale=yts;break; //m/yr - case SurfaceforcingsMassBalanceEnum: scale=yts;break; //m/yr - case SurfaceforcingsSmbPosMaxEnum: scale=yts;break; //m/yr - case SurfaceforcingsSmbPosMinEnum: scale=yts;break; //m/yr - case SurfaceforcingsAPosEnum: scale=yts;break; //m/yr - case SurfaceforcingsBPosEnum: scale=yts;break; //m/yr - case SurfaceforcingsANegEnum: scale=yts;break; //m/yr - case SurfaceforcingsBNegEnum: scale=yts;break; //m/yr + case HydrologyWaterVxEnum: scale=yts;break; //m/yr + case HydrologyWaterVyEnum: scale=yts;break; //m/yr + case VxEnum: scale=yts;break; //m/yr + case InversionVxObsEnum: scale=yts;break; //m/yr + case VyEnum: scale=yts;break; //m/yr + case InversionVyObsEnum: scale=yts;break; //m/yr + case VzEnum: scale=yts;break; //m/yr + case InversionVzObsEnum: scale=yts;break; //m/yr + case VelEnum: scale=yts;break; //m/yr + case InversionVelObsEnum: scale=yts;break; //m/yr + case BalancethicknessThickeningRateEnum: scale=yts;break; //m/yr + case BasalforcingsMeltingRateEnum: scale=yts;break; //m/yr + case BasalforcingsMeltingRateCorrectionEnum: scale=yts;break; //m/yr + case SurfaceforcingsPrecipitationEnum: scale=yts;break; //m/yr + case SurfaceforcingsMassBalanceEnum: scale=yts;break; //m/yr + case SurfaceforcingsSmbPosMaxEnum: scale=yts;break; //m/yr + case SurfaceforcingsSmbPosMinEnum: scale=yts;break; //m/yr + case SurfaceforcingsAPosEnum: scale=yts;break; //m/yr + case SurfaceforcingsBPosEnum: scale=yts;break; //m/yr + case SurfaceforcingsANegEnum: scale=yts;break; //m/yr + case SurfaceforcingsBNegEnum: scale=yts;break; //m/yr case MisfitEnum: scale=pow(yts,2);break; //(m/yr)^2 case MassFluxEnum: scale=pow((IssmDouble)10,-12)*yts;break; // (GigaTon/year) Index: /issm/branches/trunk-jpl-damage/src/c/solutions/controltao_core.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/solutions/controltao_core.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/solutions/controltao_core.cpp (revision 12946) @@ -21,4 +21,5 @@ typedef struct { FemModel* femmodel; + double* J; } AppCtx; @@ -69,9 +70,10 @@ GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound"); GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound"); - TaoSetInitialVector(tao,X->vector); - TaoSetVariableBounds(tao,XL->vector,XU->vector); + TaoSetInitialVector(tao,X->pvector->vector); + TaoSetVariableBounds(tao,XL->pvector->vector,XU->pvector->vector); xdelete(&XL); xdelete(&XU); + user.J=(double*)xcalloc((maxiter+5),sizeof(double)); user.femmodel=femmodel; TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user); @@ -81,9 +83,10 @@ TaoSolve(tao); TaoView(tao,PETSC_VIEWER_STDOUT_WORLD); - TaoGetSolutionVector(tao,&X->vector); + TaoGetSolutionVector(tao,&X->pvector->vector); SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X); for(int i=0;ielements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_list[i]); } + femmodel->results->AddObject(new DoubleVecExternalResult(femmodel->results->Size()+1,JEnum,user.J,maxiter+3,1,0)); /*Finalize*/ @@ -142,5 +145,5 @@ /*Compute gradient*/ Gradjx(&gradient,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); - VecCopy(gradient->vector,G); xdelete(&gradient); + VecCopy(gradient->pvector->vector,G); xdelete(&gradient); VecScale(G,-1.); @@ -166,4 +169,5 @@ if(its==0) _pprintLine_("-----------------------------------+-----------------------"); _pprintString_(setw(4)<J[its]=f; /*Retrieve objective functions independently*/ Index: /issm/branches/trunk-jpl-damage/src/c/solutions/transient_core.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/solutions/transient_core.cpp (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/c/solutions/transient_core.cpp (revision 12946) @@ -24,5 +24,5 @@ /*parameters: */ IssmDouble starttime,finaltime,dt,yts; - bool isdiagnostic,isprognostic,isthermal,isgroundingline,isenthalpy; + bool isdiagnostic,isprognostic,isthermal,isgroundingline,isenthalpy,isdelta18o; bool save_results,dakota_analysis; bool time_adapt=false; @@ -53,4 +53,5 @@ femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum); if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,TransientRequestedOutputsEnum); + femmodel->parameters->FindParam(&isdelta18o,SurfaceforcingsIsdelta18oEnum); /*initialize: */ @@ -141,4 +142,9 @@ RequestedOutputsx(femmodel->results,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,requested_outputs,numoutputs); + if(isdelta18o){ + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceforcingsMonthlytemperaturesEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceforcingsPrecipitationEnum); + } + if(VerboseSolution()) _pprintLine_(" saving temporary results"); OutputResultsx(femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,femmodel->results); Index: /issm/branches/trunk-jpl-damage/src/m/Makefile.am =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/Makefile.am (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/Makefile.am (revision 12946) @@ -11,5 +11,4 @@ ./model/partition/*.m \ ./model/solvers/*.m \ - ./planet/*.m \ ./qmu/*.m \ ./utils/*.m \ Index: /issm/branches/trunk-jpl-damage/src/m/classes/clusters/acenet.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/classes/clusters/acenet.m (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/classes/clusters/acenet.m (revision 12946) @@ -18,5 +18,5 @@ time=10; % codepath='/usr/local/issm-r11321/bin'; % this one is for issm on acenet global - codepath='/home/klemorza/issm/bin'; % this one is for issm on my acenet directory + codepath='PATH'; % this one is for issm on my acenet directory executionpath='/home/klemorza/issm/execution'; %}}} @@ -86,7 +86,4 @@ compressstring = [compressstring ' ' filelist{i}]; end - if cluster.interactive, - compressstring = [compressstring ' ' modelname '.errlog ' modelname '.outlog ']; - end system(compressstring); Index: /issm/branches/trunk-jpl-damage/src/m/classes/clusters/generic.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/classes/clusters/generic.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/classes/clusters/generic.py (revision 12946) @@ -14,7 +14,7 @@ GENERIC cluster class definition - Usage: - cluster=generic('name','astrid','np',3); - cluster=generic('name',oshostname(),'np',3,'login','username'); + Usage: + cluster=generic('name','astrid','np',3); + cluster=generic('name',oshostname(),'np',3,'login','username'); """ Index: /issm/branches/trunk-jpl-damage/src/m/classes/diagnostic.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/classes/diagnostic.m (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/classes/diagnostic.m (revision 12946) @@ -73,5 +73,5 @@ md = checkfield(md,'diagnostic.spcvy','forcing',1); if md.mesh.dimension==3, md = checkfield(md,'diagnostic.spcvz','forcing',1); end - md = checkfield(md,'diagnostic.restol','size',[1 1],'>',0); + md = checkfield(md,'diagnostic.restol','size',[1 1],'>',0,'NaN',1); md = checkfield(md,'diagnostic.reltol','size',[1 1]); md = checkfield(md,'diagnostic.abstol','size',[1 1]); Index: /issm/branches/trunk-jpl-damage/src/m/classes/mesh.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/classes/mesh.m (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/classes/mesh.m (revision 12946) @@ -57,5 +57,5 @@ function obj = setdefaultparameters(obj) % {{{ - %the connectivity is the avergaded number of nodes linked to a + %the connectivity is the averaged number of nodes linked to a %given node through an edge. This connectivity is used to initially %allocate memory to the stiffness matrix. A value of 16 seems to Index: /issm/branches/trunk-jpl-damage/src/m/classes/mesh.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/classes/mesh.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/classes/mesh.py (revision 12946) @@ -1,4 +1,8 @@ #module imports +import numpy from fielddisplay import fielddisplay +from checkfield import * +from EnumDefinitions import * +from MatlabFuncs import * class mesh: @@ -65,10 +69,5 @@ else: string="\n%s"%(" Elements and vertices:") - string="%s\n%s"%(string,fielddisplay(obj,"numberofelements","number of elements")) - - - - string="%s\n%s"%(string,fielddisplay(obj,"numberofvertices","number of vertices")) string="%s\n%s"%(string,fielddisplay(obj,"elements","index into (x,y,z), coordinates of the vertices")) @@ -80,5 +79,4 @@ string="%s%s"%(string,"\n Properties:") - string="%s\n%s"%(string,fielddisplay(obj,"dimension","mesh dimension (2d or 3d)")) string="%s\n%s"%(string,fielddisplay(obj,"numberoflayers","number of extrusion layers")) @@ -92,5 +90,4 @@ string="%s\n%s"%(string,fielddisplay(obj,"lowerelements","lower element list (NaN for element on the lower layer")) string="%s\n%s"%(string,fielddisplay(obj,"vertexonboundary","vertices on the boundary of the domain flag list")) - string="%s\n%s"%(string,fielddisplay(obj,"segments","edges on domain boundary (vertex1 vertex2 element)")) string="%s\n%s"%(string,fielddisplay(obj,"segmentmarkers","number associated to each segment")) @@ -100,5 +97,4 @@ string="%s%s"%(string,"\n Extracted model:") - string="%s\n%s"%(string,fielddisplay(obj,"extractedvertices","vertices extracted from the model")) string="%s\n%s"%(string,fielddisplay(obj,"extractedelements","elements extracted from the model")) @@ -114,5 +110,5 @@ # {{{setdefaultparameters - #the connectivity is the avergaded number of nodes linked to a + #the connectivity is the averaged number of nodes linked to a #given node through an edge. This connectivity is used to initially #allocate memory to the stiffness matrix. A value of 16 seems to @@ -124,2 +120,75 @@ #}}} + def checkconsistency(self,md,solution,analyses): # {{{ + + md = checkfield(md,'mesh.x','NaN',1,'size',[md.mesh.numberofvertices]) + md = checkfield(md,'mesh.y','NaN',1,'size',[md.mesh.numberofvertices]) + md = checkfield(md,'mesh.z','NaN',1,'size',[md.mesh.numberofvertices]) + md = checkfield(md,'mesh.elements','NaN',1,'>',0,'values',range(1,md.mesh.numberofvertices+1)) + if md.mesh.dimension==2: + md = checkfield(md,'mesh.elements','size',[md.mesh.numberofelements,3]) + else: + md = checkfield(md,'mesh.elements','size',[md.mesh.numberofelements,6]) + if any(numpy.logical_not(ismember(range(1,md.mesh.numberofvertices+1),md.mesh.elements))): + md = checkmessage(md,"orphan nodes have been found. Check the mesh outline") + md = checkfield(md,'mesh.dimension','values',[2,3]) + md = checkfield(md,'mesh.numberoflayers','>=',0) + md = checkfield(md,'mesh.numberofelements','>',0) + md = checkfield(md,'mesh.numberofvertices','>',0) + #no checks for numberofedges lat long and hemisphere + md = checkfield(md,'mesh.elementonbed','size',[md.mesh.numberofelements],'values',[0,1]) + md = checkfield(md,'mesh.elementonsurface','size',[md.mesh.numberofelements],'values',[0,1]) + md = checkfield(md,'mesh.vertexonbed','size',[md.mesh.numberofvertices],'values',[0,1]) + md = checkfield(md,'mesh.vertexonsurface','size',[md.mesh.numberofvertices],'values',[0,1]) + if md.mesh.dimension==2: + md = checkfield(md,'mesh.average_vertex_connectivity','>=',9,'message',"'mesh.average_vertex_connectivity' should be at least 9 in 2d") + else: + md = checkfield(md,'mesh.average_vertex_connectivity','>=',24,'message',"'mesh.average_vertex_connectivity' should be at least 24 in 3d") + md = checkfield(md,'mesh.elementconnectivity','size',[md.mesh.numberofelements,3],'NaN',1) + + #Solution specific checks + if solution==PrognosticSolutionEnum: + if md.prognostic.stabilization==3: + md = checkfield(md,'mesh.dimension','values',2,'message',"Discontinuous Galerkin only supported for 2d meshes") + md = checkfield(md,'mesh.edges','size',[float('NaN'),4]) + md = checkfield(md,'mesh.edges[:,1:3]','>',0) + elif solution==BalancethicknessSolutionEnum: + if md.balancethickness.stabilization==3: + md = checkfield(md,'mesh.dimension','values',2,'message',"Discontinuous Galerkin only supported for 2d meshes") + md = checkfield(md,'mesh.edges','size',[float('NaN'),4]) + md = checkfield(md,'mesh.edges[:,1:3]','>',0) + elif solution==TransientSolutionEnum: + if md.transient.isprognostic and md.prognostic.stabilization==3: + md = checkfield(md,'mesh.dimension','values',2,'message',"Discontinuous Galerkin only supported for 2d meshes") + md = checkfield(md,'mesh.edges','size',[float('NaN'),4]) + md = checkfield(md,'mesh.edges[:,1:3]','>',0) + elif solution==ThermalSolutionEnum: + md = checkfield(md,'mesh.dimension','values',3,'message','thermal solution only supported on 3d meshes') + + return md + # }}} + + def marshall(self,fid): # {{{ + WriteData(fid,'object',self,'fieldname','x','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'fieldname','y','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'fieldname','z','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'fieldname','elements','format','DoubleMat','mattype',2) + WriteData(fid,'object',self,'fieldname','dimension','format','Integer') + WriteData(fid,'object',self,'fieldname','numberoflayers','format','Integer') + WriteData(fid,'object',self,'fieldname','numberofelements','format','Integer') + WriteData(fid,'object',self,'fieldname','numberofvertices','format','Integer') + WriteData(fid,'object',self,'fieldname','numberofedges','format','Integer') + WriteData(fid,'object',self,'fieldname','elementonbed','format','BooleanMat','mattype',2) + WriteData(fid,'object',self,'fieldname','elementonsurface','format','BooleanMat','mattype',2) + WriteData(fid,'object',self,'fieldname','vertexonbed','format','BooleanMat','mattype',1) + WriteData(fid,'object',self,'fieldname','vertexonsurface','format','BooleanMat','mattype',1) + WriteData(fid,'object',self,'fieldname','lowerelements','format','DoubleMat','mattype',2) + WriteData(fid,'object',self,'fieldname','upperelements','format','DoubleMat','mattype',2) + WriteData(fid,'object',self,'fieldname','edges','format','DoubleMat','mattype',3) + WriteData(fid,'object',self,'fieldname','elementconnectivity','format','DoubleMat','mattype',3) + WriteData(fid,'object',self,'fieldname','average_vertex_connectivity','format','Integer') + WriteData(fid,'object',self,'fieldname','elements2d','format','DoubleMat','mattype',3) + WriteData(fid,'object',self,'fieldname','numberofvertices2d','format','Integer') + WriteData(fid,'object',self,'fieldname','numberofelements2d','format','Integer') + # }}} + Index: /issm/branches/trunk-jpl-damage/src/m/classes/model/model.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/classes/model/model.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/classes/model/model.py (revision 12946) @@ -32,5 +32,7 @@ from miscellaneous import miscellaneous from private import private +from EnumDefinitions import * #}}} + class model(object): #properties @@ -51,7 +53,12 @@ self.debug = debug() - self.verbose = verbose() + self.verbose = verbose('solution',True,'qmu',True,'control',True) self.settings = settings() self.solver = solver() +# need EnumToString for these to work +# if ismumps: +# md.solver.addoptions(DiagnosticVertAnalysisEnum(),mumpsoptions) +# else: +# md.solver.addoptions(DiagnosticVertAnalysisEnum(),iluasmoptions) self.cluster = generic() @@ -70,9 +77,10 @@ self.qmu = qmu() - self.results = []; + self.results = {} self.radaroverlay = radaroverlay() self.miscellaneous = miscellaneous() self.private = private() #}}} + def __repr__(obj): # {{{ Display @@ -111,4 +119,10 @@ string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("radaroverlay","[%s,%s]" % ("1x1",obj.radaroverlay.__class__.__name__),"radar image for plot overlay")) string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("miscellaneous","[%s,%s]" % ("1x1",obj.miscellaneous.__class__.__name__),"miscellaneous fields")) - return string; - #}}} + return string + # }}} + + def checkmessage(self,string): # {{{ + print ("model not consistent: %s" % string) + self.private.isconsistent=False + # }}} + Index: /issm/branches/trunk-jpl-damage/src/m/classes/pairoptions.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/classes/pairoptions.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/classes/pairoptions.py (revision 12946) @@ -5,7 +5,7 @@ PAIROPTIONS class definition - Usage: - pairoptions=pairoptions(); - pairoptions=pairoptions('module',true,'solver',false); + Usage: + pairoptions=pairoptions(); + pairoptions=pairoptions('module',true,'solver',false); """ Index: /issm/branches/trunk-jpl-damage/src/m/classes/planetmesh.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/classes/planetmesh.m (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/m/classes/planetmesh.m (revision 12946) @@ -0,0 +1,115 @@ +%MESH class definition +% +% Usage: +% planetmesh=planetmesh(); + +classdef planetmesh + properties (SetAccess=public) + r = NaN; + theta = NaN; + phi = NaN + elements = NaN + dimension = 0; + numberoflayers = 0; + numberofelements = 0; + numberofvertices = 0; + + lat = NaN + long = NaN + + vertexconnectivity = NaN + elementconnectivity = NaN + average_vertex_connectivity = 0; + end + methods + function obj = planetmesh(varargin) % {{{ + switch nargin + case 0 + obj=setdefaultparameters(obj); + otherwise + error('constructor not supported'); + end + end % }}} + function obj = setdefaultparameters(obj) % {{{ + + %the connectivity is the avergaded number of nodes linked to a + %given node through an edge. This connectivity is used to initially + %allocate memory to the stiffness matrix. A value of 16 seems to + %give a good memory/time ration. This value can be checked in + %trunk/test/Miscellaneous/runme.m + obj.average_vertex_connectivity=25; + end % }}} + function md = checkconsistency(obj,md,solution,analyses) % {{{ + + md = checkfield(md,'planetmesh.r','NaN',1,'size',[md.planetmesh.numberofvertices 1]); + md = checkfield(md,'planetmesh.theta','NaN',1,'size',[md.planetmesh.numberofvertices 1]); + md = checkfield(md,'planetmesh.phi','NaN',1,'size',[md.planetmesh.numberofvertices 1]); + md = checkfield(md,'planetmesh.elements','NaN',1,'>',0,'values',1:md.planetmesh.numberofvertices); + if(md.planetmesh.dimension==2), + md = checkfield(md,'planetmesh.elements','size',[md.planetmesh.numberofelements 3]); + else + md = checkfield(md,'planetmesh.elements','size',[md.planetmesh.numberofelements 6]); + end + if any(~ismember(1:md.planetmesh.numberofvertices,sort(unique(md.planetmesh.elements(:))))); + md = checkmessage(md,'orphan nodes have been found. Check the planetmesh outline'); + end + md = checkfield(md,'planetmesh.dimension','values',[2 3]); + md = checkfield(md,'planetmesh.numberoflayers','>=',0); + md = checkfield(md,'planetmesh.numberofelements','>',0); + md = checkfield(md,'planetmesh.numberofvertices','>',0); + %no checks for numberofedges lat long and hemisphere + if (md.planetmesh.dimension==2), + md = checkfield(md,'planetmesh.average_vertex_connectivity','>=',9,'message','''planetmesh.average_vertex_connectivity'' should be at least 9 in 2d'); + else + md = checkfield(md,'planetmesh.average_vertex_connectivity','>=',24,'message','''planetmesh.average_vertex_connectivity'' should be at least 24 in 3d'); + end + md = checkfield(md,'planetmesh.elementconnectivity','size',[md.planetmesh.numberofelements 3],'NaN',1); + + %Solution specific checks + switch(solution), + case PrognosticSolutionEnum, + if md.prognostic.stabilization==3, + md = checkfield(md,'planetmesh.dimension','values',2,'message','Discontinuous Galerkin only supported for 2d planetmeshes'); + end + case TransientSolutionEnum, + if md.transient.isprognostic & md.prognostic.stabilization==3, + md = checkfield(md,'planetmesh.dimension','values',2,'message','Discontinuous Galerkin only supported for 2d planetmeshes'); + end + case ThermalSolutionEnum, + md = checkfield(md,'planetmesh.dimension','values',3,'message','thermal solution only supported on 3d planetmeshes'); + end + end % }}} + function disp(obj) % {{{ + disp(sprintf(' Mesh:')); + + disp(sprintf('\n Elements and vertices:')); + fielddisplay(obj,'numberofelements','number of elements'); + fielddisplay(obj,'numberofvertices','number of vertices'); + fielddisplay(obj,'elements','index into (x,y,z), coordinates of the vertices'); + fielddisplay(obj,'r','vertices r coordinate'); + fielddisplay(obj,'theta','vertices theta coordinate'); + fielddisplay(obj,'phi','vertices phi coordinate'); + + disp(sprintf('\n Properties:')); + fielddisplay(obj,'dimension','planetmesh dimension (2d or 3d)'); + fielddisplay(obj,'numberoflayers','number of extrusion layers'); + + fielddisplay(obj,'vertexconnectivity','list of vertices connected to vertex_i'); + fielddisplay(obj,'elementconnectivity','list of vertices connected to element_i'); + fielddisplay(obj,'average_vertex_connectivity','average number of vertices connected to one vertex'); + + end % }}} + function marshall(obj,fid) % {{{ + WriteData(fid,'object',obj,'fieldname','r','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'fieldname','theta','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'fieldname','phi','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'fieldname','elements','format','DoubleMat','mattype',2); + WriteData(fid,'object',obj,'fieldname','dimension','format','Integer'); + WriteData(fid,'object',obj,'fieldname','numberoflayers','format','Integer'); + WriteData(fid,'object',obj,'fieldname','numberofelements','format','Integer'); + WriteData(fid,'object',obj,'fieldname','numberofvertices','format','Integer'); + WriteData(fid,'object',obj,'fieldname','elementconnectivity','format','DoubleMat','mattype',3); + WriteData(fid,'object',obj,'fieldname','average_vertex_connectivity','format','Integer'); + end % }}} + end +end Index: /issm/branches/trunk-jpl-damage/src/m/classes/private.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/classes/private.m (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/classes/private.m (revision 12946) @@ -31,5 +31,5 @@ fielddisplay(obj,'isconsistent','is model self consistent'); fielddisplay(obj,'runtimename','name of the run launched'); - fielddisplay(obj,'bamg','structure with mesh properties construced if bamg is used to mesh the domain'); + fielddisplay(obj,'bamg','structure with mesh properties constructed if bamg is used to mesh the domain'); fielddisplay(obj,'solution','type of solution launched'); Index: /issm/branches/trunk-jpl-damage/src/m/classes/private.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/classes/private.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/classes/private.py (revision 12946) @@ -6,7 +6,8 @@ def __init__(self): # {{{ Properties - self.runtimename = '' - self.bamg = {} - self.solution = ''; + self.isconsistent = True; + self.runtimename = '' + self.bamg = {} + self.solution = ''; #set defaults @@ -17,6 +18,8 @@ # {{{ Display string=' private parameters: do not change' + + string="%s\n%s"%(string,fielddisplay(obj,'isconsistent','is model self consistent')) string="%s\n%s"%(string,fielddisplay(obj,'runtimename','name of the run launched')) - string="%s\n%s"%(string,fielddisplay(obj,'bamg','structure with mesh properties construced if bamg is used to mesh the domain')) + string="%s\n%s"%(string,fielddisplay(obj,'bamg','structure with mesh properties constructed if bamg is used to mesh the domain')) string="%s\n%s"%(string,fielddisplay(obj,'solution','type of solution launched')) return string Index: /issm/branches/trunk-jpl-damage/src/m/classes/verbose.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/classes/verbose.m (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/classes/verbose.m (revision 12946) @@ -42,5 +42,5 @@ case 1, binary=varargin{1}; - if ischar(binary), + if ischar(binary), if strcmpi(binary,'all'), binary=2^11-1; %all ones @@ -51,4 +51,6 @@ verbose=BinaryToVerbose(verbose,binary); end + elseif isnumeric(binary), + verbose=BinaryToVerbose(verbose,binary); end otherwise, @@ -88,11 +90,11 @@ %BEGINBIN2VERB - if bitand(binary,1), verbose.mprocessor=true; else verbose.mprocessor=false; end - if bitand(binary,2), verbose.module=true; else verbose.module=false; end - if bitand(binary,4), verbose.solution=true; else verbose.solution=false; end - if bitand(binary,8), verbose.solver=true; else verbose.solver=false; end - if bitand(binary,16), verbose.convergence=true; else verbose.convergence=false; end - if bitand(binary,32), verbose.control=true; else verbose.control=false; end - if bitand(binary,64), verbose.qmu=true; else verbose.qmu=false; end + verbose.mprocessor =logical(bitand(binary, 1)); + verbose.module =logical(bitand(binary, 2)); + verbose.solution =logical(bitand(binary, 4)); + verbose.solver =logical(bitand(binary, 8)); + verbose.convergence=logical(bitand(binary,16)); + verbose.control =logical(bitand(binary,32)); + verbose.qmu =logical(bitand(binary,64)); %ENDBIN2VERB Index: /issm/branches/trunk-jpl-damage/src/m/classes/verbose.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/classes/verbose.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/classes/verbose.py (revision 12946) @@ -1,9 +1,32 @@ -#module imports -from fielddisplay import fielddisplay +from pairoptions import * +from MatlabFuncs import * +from WriteData import * +from EnumDefinitions import * -class verbose: - #properties - def __init__(self): - # {{{ Properties +class verbose(object): + """ + VERBOSE class definition + + Available verbosity levels: + mprocessor : model processing + module : modules + solution : solution sequence + solver : solver info (extensive) + convergence : convergence criteria + control : control method + qmu : sensitivity analysis + + Usage: + verbose=verbose(); + verbose=verbose(3); + verbose=verbose('001100'); + verbose=verbose('module',True,'solver',False); + + WARNING: some parts of this file are Synchronized with src/c/shared/Numerics/Verbosity.h + Do not modify these sections. See src/c/shared/Numerics/README for more info + """ + + def __init__(self,*args): # {{{ + #BEGINFIELDS self.mprocessor = False self.module = False @@ -13,15 +36,94 @@ self.control = False self.qmu = False - #}}} - def __repr__(obj): - # {{{ Display - string="%s%s%s\n\n"%("class '",obj.__class__.__name__,"'=") - string="%s%s\n"%(string," %15s : %s"%("mprocessor",obj.mprocessor)) - string="%s%s\n"%(string," %15s : %s"%("module",obj.module)) - string="%s%s\n"%(string," %15s : %s"%("solution",obj.solution)) - string="%s%s\n"%(string," %15s : %s"%("solver",obj.solver)) - string="%s%s\n"%(string," %15s : %s"%("convergence",obj.convergence)) - string="%s%s\n"%(string," %15s : %s"%("control",obj.control)) - string="%s%s\n"%(string," %15s : %s"%("qmu",obj.qmu)) - return string - #}}} + #ENDFIELDS + + if not len(args): + #Don't do anything + pass + + elif len(args) == 1: + binary=args[0] + if isinstance(binary,str): + if strcmpi(binary,'all'): + binary=2**11-1 #all ones + self.BinaryToVerbose(binary) + self.solver=False #Do not use by default + else: + binary=int(binary,2) + self.BinaryToVerbose(binary) + elif isinstance(binary,(int,float)): + self.BinaryToVerbose(int(binary)) + + else: + #Use options to initialize object + self=pairoptions(*args).AssignObjectFields(self) + + #Cast to logicals + listproperties=vars(self) + for [fieldname,fieldvalue] in listproperties.iteritems(): + if isinstance(fieldvalue,bool) or isinstance(fieldvalue,(int,float)): + setattr(self,fieldname,bool(fieldvalue)) + else: + raise TypeError("verbose supported field values are logicals only (True or False)") + # }}} + + def VerboseToBinary(self): # {{{ + + #BEGINVERB2BIN + binary=0 + if self.mprocessor: + binary=binary | 1 + if self.module: + binary=binary | 2 + if self.solution: + binary=binary | 4 + if self.solver: + binary=binary | 8 + if self.convergence: + binary=binary | 16 + if self.control: + binary=binary | 32 + if self.qmu: + binary=binary | 64 + #ENDVERB2BIN + + return binary + # }}} + + def BinaryToVerbose(self,binary): # {{{ + + #BEGINBIN2VERB + self.mprocessor =bool(binary & 1) + self.module =bool(binary & 2) + self.solution =bool(binary & 4) + self.solver =bool(binary & 8) + self.convergence=bool(binary & 16) + self.control =bool(binary & 32) + self.qmu =bool(binary & 64) + #ENDBIN2VERB + # }}} + + def checkconsistency(self,md,solution,analyses): # {{{ + return md + # }}} + + def __repr__(self): # {{{ + + #BEGINDISP + s ="class '%s' = \n" % type(self) + s+=" %15s : %s\n" % ('mprocessor',self.mprocessor) + s+=" %15s : %s\n" % ('module',self.module) + s+=" %15s : %s\n" % ('solution',self.solution) + s+=" %15s : %s\n" % ('solver',self.solver) + s+=" %15s : %s\n" % ('convergence',self.convergence) + s+=" %15s : %s\n" % ('control',self.control) + s+=" %15s : %s\n" % ('qmu',self.qmu) + #ENDDISP + + return s + # }}} + + def marshall(self,fid): # {{{ + WriteData(fid,'data',VerboseToBinary(self),'enum',VerboseEnum,'format','Integer') + # }}} + Index: /issm/branches/trunk-jpl-damage/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/enum/EnumDefinitions.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/enum/EnumDefinitions.py (revision 12946) @@ -3817,4 +3817,24 @@ return StringToEnum('TotalSmb') +def ThicknessAlongGradientEnum(): + """ + THICKNESSALONGGRADIENTENUM - Enum of ThicknessAlongGradient + + Usage: + macro=ThicknessAlongGradientEnum() + """ + + return StringToEnum('ThicknessAlongGradient') + +def ThicknessAcrossGradientEnum(): + """ + THICKNESSACROSSGRADIENTENUM - Enum of ThicknessAcrossGradient + + Usage: + macro=ThicknessAcrossGradientEnum() + """ + + return StringToEnum('ThicknessAcrossGradient') + def P0Enum(): """ Index: /issm/branches/trunk-jpl-damage/src/m/enum/MaximumNumberOfEnums.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/enum/MaximumNumberOfEnums.m (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/enum/MaximumNumberOfEnums.m (revision 12946) @@ -9,3 +9,3 @@ % macro=MaximumNumberOfEnums() -macro=460; +macro=458; Index: /issm/branches/trunk-jpl-damage/src/m/enum/ThicknessAcrossGradientEnum.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/enum/ThicknessAcrossGradientEnum.m (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/m/enum/ThicknessAcrossGradientEnum.m (revision 12946) @@ -0,0 +1,11 @@ +function macro=ThicknessAcrossGradientEnum() +%THICKNESSACROSSGRADIENTENUM - Enum of ThicknessAcrossGradient +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=ThicknessAcrossGradientEnum() + +macro=StringToEnum('ThicknessAcrossGradient'); Index: /issm/branches/trunk-jpl-damage/src/m/enum/ThicknessAlongGradientEnum.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/enum/ThicknessAlongGradientEnum.m (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/m/enum/ThicknessAlongGradientEnum.m (revision 12946) @@ -0,0 +1,11 @@ +function macro=ThicknessAlongGradientEnum() +%THICKNESSALONGGRADIENTENUM - Enum of ThicknessAlongGradient +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=ThicknessAlongGradientEnum() + +macro=StringToEnum('ThicknessAlongGradient'); Index: /issm/branches/trunk-jpl-damage/src/m/model/AnalysisConfiguration.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/AnalysisConfiguration.py (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/m/model/AnalysisConfiguration.py (revision 12946) @@ -0,0 +1,59 @@ +from EnumDefinitions import * + +def AnalysisConfiguration(solutiontype): + """ + ANALYSISCONFIGURATION - return type of analyses, number of analyses + + Usage: + [analyses, numanalyses]=AnalysisConfiguration(solutiontype); + """ + + if solutiontype == DiagnosticSolutionEnum: + numanalyses=5 + analyses=[DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticHutterAnalysisEnum,SurfaceSlopeAnalysisEnum,BedSlopeAnalysisEnum] + + elif solutiontype == SteadystateSolutionEnum: + numanalyses=7 + analyses=[DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticHutterAnalysisEnum,SurfaceSlopeAnalysisEnum,BedSlopeAnalysisEnum,ThermalAnalysisEnum,MeltingAnalysisEnum] + + elif solutiontype == ThermalSolutionEnum: + numanalyses=2 + analyses=[ThermalAnalysisEnum,MeltingAnalysisEnum] + + elif solutiontype == EnthalpySolutionEnum: + numanalyses=1 + analyses=[EnthalpyAnalysisEnum] + + elif solutiontype == PrognosticSolutionEnum: + numanalyses=1 + analyses=[PrognosticAnalysisEnum] + + elif solutiontype == BalancethicknessSolutionEnum: + numanalyses=1 + analyses=[BalancethicknessAnalysisEnum] + + elif solutiontype == SurfaceSlopeSolutionEnum: + numanalyses=1 + analyses=[SurfaceSlopeAnalysisEnum] + + elif solutiontype == BedSlopeSolutionEnum: + numanalyses=1 + analyses=[BedSlopeAnalysisEnum] + + elif solutiontype == TransientSolutionEnum: + numanalyses=9 + analyses=[DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticHutterAnalysisEnum,SurfaceSlopeAnalysisEnum,BedSlopeAnalysisEnum,ThermalAnalysisEnum,MeltingAnalysisEnum,EnthalpyAnalysisEnum,PrognosticAnalysisEnum] + + elif solutiontype == FlaimSolutionEnum: + numanalyses=1 + analyses=[FlaimAnalysisEnum] + + elif solutiontype == HydrologySolutionEnum: + numanalyses=3 + analyses=[BedSlopeAnalysisEnum,SurfaceSlopeAnalysisEnum,HydrologyAnalysisEnum] + + else: + raise TypeError("solution type: '%s' not supported yet!" % EnumToString(solutiontype)) + + return analyses,numanalyses + Index: /issm/branches/trunk-jpl-damage/src/m/model/MatlabFuncs.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/MatlabFuncs.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/model/MatlabFuncs.py (revision 12946) @@ -1,2 +1,7 @@ +def oshostname(): + import socket + + return socket.gethostname().lower().split('.')[0] + def strcmp(s1,s2): Index: /issm/branches/trunk-jpl-damage/src/m/model/MatlabProcessPatch.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/MatlabProcessPatch.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/model/MatlabProcessPatch.py (revision 12946) @@ -3,6 +3,6 @@ PROCESSPATCH - create a structure from a patch - Usage: - Result=ProcessPatch(Result); + Usage: + Result=ProcessPatch(Result); """ Index: /issm/branches/trunk-jpl-damage/src/m/model/WriteData.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/WriteData.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/model/WriteData.py (revision 12946) @@ -8,10 +8,10 @@ WRITEDATA - write model field in binary file - Usage: - WriteData(fid,*args) + Usage: + WriteData(fid,varargin) """ #process options - options=pairoptions(args) + options=pairoptions(*args) #Get data properties Index: /issm/branches/trunk-jpl-damage/src/m/model/addnote.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/addnote.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/model/addnote.py (revision 12946) @@ -2,10 +2,10 @@ """ ADDNOTE - add a note to the existing model notes field - - Usage: - md=addnote(md,string); - - Example: - md=addnote(md,'Pine Island, Geometry of 2007'); + + Usage: + md=addnote(md,string); + + Example: + md=addnote(md,'Pine Island, Geometry of 2007'); """ Index: /issm/branches/trunk-jpl-damage/src/m/model/ismodelselfconsistent.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/ismodelselfconsistent.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/model/ismodelselfconsistent.py (revision 12946) @@ -1,12 +1,14 @@ +from AnalysisConfiguration import * + def ismodelselfconsistent(md): """ ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem. - Usage: - ismodelselfconsistent(md), + Usage: + ismodelselfconsistent(md), """ #initialize consistency as true - md.private.isconsistent=true + md.private.isconsistent=True #Get solution and associated analyses @@ -15,6 +17,6 @@ #Go through a model fields, check that it is a class, and call checkconsistency - fields=vars('model') - for field in fields: + fields=vars(md) + for field in fields.iterkeys(): #Some properties do not need to be checked @@ -24,5 +26,5 @@ #Check that current field is an object if not hasattr(getattr(md,field),'checkconsistency'): - md=checkmessage(md,"field '%s' is not an object." % field) + md.checkmessage("field '%s' is not an object." % field) #Check consistency of the object Index: /issm/branches/trunk-jpl-damage/src/m/model/loadresultsfromcluster.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/loadresultsfromcluster.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/model/loadresultsfromcluster.py (revision 12946) @@ -8,6 +8,6 @@ LOADRESULTSFROMCLUSTER - load results of solution sequence from cluster - Usage: - md=loadresultsfromcluster(md,runtimename); + Usage: + md=loadresultsfromcluster(md,runtimename); """ Index: /issm/branches/trunk-jpl-damage/src/m/model/loadresultsfromdisk.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/loadresultsfromdisk.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/model/loadresultsfromdisk.py (revision 12946) @@ -5,6 +5,6 @@ LOADRESULTSFROMDISK - load results of solution sequence from disk file "filename" - Usage: - md=loadresultsfromdisk(md=False,filename=False); + Usage: + md=loadresultsfromdisk(md=False,filename=False); """ Index: /issm/branches/trunk-jpl-damage/src/m/model/marshall.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/marshall.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/model/marshall.py (revision 12946) @@ -5,9 +5,9 @@ MARSHALL - outputs a compatible binary file from @model md, for certain solution type. - The routine creates a compatible binary file from @model md - This binary file will be used for parallel runs in JPL-package + The routine creates a compatible binary file from @model md + This binary file will be used for parallel runs in JPL-package - Usage: - marshall(md) + Usage: + marshall(md) """ @@ -21,10 +21,10 @@ #First, write MaximumNumberOfEnum to make sure that the Enums are synchronized - WriteData(fid,'enum',MaximumNumberOfEnums(),'data',true,'format','Boolean') + WriteData(fid,'enum',MaximumNumberOfEnums(),'data',True,'format','Boolean') #Go through all model fields: check that it is a class and call checkconsistency fields=vars(md) - for field in fields.interkeys(): + for field in fields.iterkeys(): #Some properties do not need to be marshalled Index: /issm/branches/trunk-jpl-damage/src/m/model/mesh/planet/mesh_refine_tri4.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/mesh/planet/mesh_refine_tri4.m (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/m/model/mesh/planet/mesh_refine_tri4.m (revision 12946) @@ -0,0 +1,146 @@ +function [ FV ] = mesh_refine_tri4(FV) + +% mesh_refine_tri4 - creates 4 triangle from each triangle of a mesh +% +% [ FV ] = mesh_refine_tri4( FV ) +% +% FV.vertices - mesh vertices (Nx3 matrix) +% FV.faces - faces with indices into 3 rows +% of FV.vertices (Mx3 matrix) +% +% For each face, 3 new vertices are created at the +% triangle edge midpoints. Each face is divided into 4 +% faces and returned in FV. +% +% B +% /\ +% / \ +% a/____\b Construct new triangles +% /\ /\ [A,a,c] +% / \ / \ [a,B,b] +% /____\/____\ [c,b,C] +% A c C [a,b,c] +% +% It is assumed that the vertices are listed in clockwise order in +% FV.faces (A,B,C above), as viewed from the outside in a RHS coordinate +% system. +% +% See also: mesh_refine, sphere_tri, sphere_project +% + + +% ---this method is not implemented, but the idea here remains... +% This can be done until some minimal distance (D) of the mean +% distance between vertices of all triangles is achieved. If +% no D argument is given, the function refines the mesh once. +% Alternatively, it could be done until some minimum mean +% area of faces is achieved. As is, it just refines once. + + +% $Revision: 1.1 $ $Date: 2004/11/12 01:32:35 $ + +% Licence: GNU GPL, no implied or express warranties +% History: 05/2002, Darren.Weber_at_radiology.ucsf.edu, created +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +tic; +fprintf('...refining mesh (tri4)...') + +% NOTE +% The centroid is located one third of the way from each vertex to +% the midpoint of the opposite side. Each median divides the triangle +% into two equal areas; all the medians together divide it into six +% equal parts, and the lines from the median point to the vertices +% divide the whole into three equivalent triangles. + +% Each input triangle with vertices labelled [A,B,C] as shown +% below will be turned into four new triangles: +% +% Make new midpoints +% a = (A+B)/2 +% b = (B+C)/2 +% c = (C+A)/2 +% +% B +% /\ +% / \ +% a/____\b Construct new triangles +% /\ /\ [A,a,c] +% / \ / \ [a,B,b] +% /____\/____\ [c,b,C] +% A c C [a,b,c] +% + +% Initialise a new vertices and faces matrix +Nvert = size(FV.vertices,1); +Nface = size(FV.faces,1); +V2 = zeros(Nface*3,3); +F2 = zeros(Nface*4,3); + +for f = 1:Nface, + + % Get the triangle vertex indices + NA = FV.faces(f,1); + NB = FV.faces(f,2); + NC = FV.faces(f,3); + + % Get the triangle vertex coordinates + A = FV.vertices(NA,:); + B = FV.vertices(NB,:); + C = FV.vertices(NC,:); + + % Now find the midpoints between vertices + a = (A + B) ./ 2; + b = (B + C) ./ 2; + c = (C + A) ./ 2; + + % Find the length of each median + %A2blen = sqrt ( sum( (A - b).^2, 2 ) ); + %B2clen = sqrt ( sum( (B - c).^2, 2 ) ); + %C2alen = sqrt ( sum( (C - a).^2, 2 ) ); + + % Store the midpoint vertices, while + % checking if midpoint vertex already exists + [FV, Na] = mesh_find_vertex(FV,a); + [FV, Nb] = mesh_find_vertex(FV,b); + [FV, Nc] = mesh_find_vertex(FV,c); + + % Create new faces with orig vertices plus midpoints + F2(f*4-3,:) = [ NA, Na, Nc ]; + F2(f*4-2,:) = [ Na, NB, Nb ]; + F2(f*4-1,:) = [ Nc, Nb, NC ]; + F2(f*4-0,:) = [ Na, Nb, Nc ]; + +end + +% Replace the faces matrix +FV.faces = F2; + +t=toc; fprintf('done (%5.2f sec)\n',t); + +return + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function [FV, N] = mesh_find_vertex(FV,vertex) + + Vn = size(FV.vertices,1); + Va = repmat(vertex,Vn,1); + Vexist = find( FV.vertices(:,1) == Va(:,1) & ... + FV.vertices(:,2) == Va(:,2) & ... + FV.vertices(:,3) == Va(:,3) ); + if Vexist, + if size(Vexist) == [1,1], + N = Vexist; + else, + msg = sprintf('replicated vertices'); + error(msg); + end + else + FV.vertices(end+1,:) = vertex; + N = size(FV.vertices,1); + end + +return Index: /issm/branches/trunk-jpl-damage/src/m/model/mesh/planet/planetmesher.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/mesh/planet/planetmesher.m (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/m/model/mesh/planet/planetmesher.m (revision 12946) @@ -0,0 +1,35 @@ +function md=planetmesher(md,varargin) +%PLANETMESHER - create planet mesh using several packages +% +% This routine creates a planet mesh using several packages: a custom made +% ISSM mesher (called planetmixedmesh.m, method 'mixed') and the planettrimesh +% (method 'tria') +% where md is a @planet object, varargin is a list of options +% +% Usage: md=planetmesh(md,'method','mixed','radius',6378000,'angleresol',1); +% md=planetmesh(md,'method','tria','shape','iso','radius',6378000,'refinement',5); +% +% + + %recover options + options=pairoptions(varargin{:}); + + method=getfieldvalue(options,'method','mixed'); + + if strcmpi(method,'mixed'), + %recover radius and angleresol: + radius=getfieldvalue(options,'radius',6378000); %earth radius as default + angleresol=getfieldvalue(options,'angleresol',10); %10 degree resolution + + %call mixed mesh + md=planetmixedmesh(md,radius,angleresol); + else + %recover radius, shape and level of refinmenet + radius=getfieldvalue(options,'radius',6378000); %earth radius as default + refinement=getfieldvalue(options,'refinement',5); %refine 5 times + shape=getfieldvalue(options,'shape','ico'); + + %call triangular mesh + md=planettrimesh(md,shape,radius,refinement); + end +end Index: /issm/branches/trunk-jpl-damage/src/m/model/mesh/planet/planetmixedmesh.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/mesh/planet/planetmixedmesh.m (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/m/model/mesh/planet/planetmixedmesh.m (revision 12946) @@ -0,0 +1,99 @@ +function md=planetmixedmesh(md,radius,angleresol) +%PLANETMIXEDMESH: build 2d shell mesh +% +% Usage: md=planetmixedshell(md,radius,angleresol) +% + +conv=pi/180; + +r=radius; +theta=(0:angleresol:360)'; +theta=theta*conv; +phi=(0:angleresol:180)'; +phi=phi(2:end-1)*conv; + +nphi=length(phi); +nthetha=length(theta); + +x=zeros(nphi*nthetha,1); +y=zeros(nphi*nthetha,1); +z=zeros(nphi*nthetha,1); + +for i=1:nphi, + phii=phi(i); + x((i-1)*nthetha+1:i*nthetha)=r.*cos(theta).*sin(phii); + y((i-1)*nthetha+1:i*nthetha)=r.*sin(theta).*sin(phii); + z((i-1)*nthetha+1:i*nthetha)=r.*cos(phii); +end + +quads=zeros((nthetha-1)*(nphi-1),4); + +%build elements +for i=1:nphi-1, + for j=1:nthetha-1, + count=(i-1)*(nthetha-1)+j; + + A=(i-1)*nthetha+j; + B=(i-1)*nthetha+j+1; + C=(i)*nthetha+j+1; + D=(i)*nthetha+j; + + quads(count,:)=[A B C D]; + end +end + +%now, add polar elements: +%first north: phi = 0; +x=[r.*cos(0).*sin(0);x]; +y=[r.*sin(0).*sin(0);y]; +z=[r.*cos(0);z]; + +%add 1 to quads: +quads=quads+1; + +%add tria elements: +trias=zeros(nthetha-1,4); + +for i=1:nthetha-1, + A=1; + B=i+1; + C=i+2; + trias(i,:)=[A B C NaN]; +end + +quads=[trias;quads]; + + +% now add south pole: +phii=180*conv; +x=[x;r.*cos(phii).*sin(phii)]; +y=[y;r.*sin(phii).*sin(phii)]; +z=[z;r.*cos(phii)]; +nods=length(x); + +%add tria elements: +trias=zeros(nthetha-1,4); + +start=nods-nthetha; +for i=1:nthetha-1, + A=start+i-1; + B=start+i; + C=nods; + trias(i,:)=[A B C NaN]; +end + +quads=[quads;trias]; + + +md.mesh.elements=quads; +md.x=x; +md.y=y; +md.z=z; +md.r=sqrt(x.^2+y.^2+z.^2); +md.theta=acos(z./r); +md.phi=atan2(y,x); + +md.mesh.numberofvertices=length(md.x); +md.mesh.numberofelements=size(md.mesh.elements,1); + +md.mesh.dimension=3; Index: /issm/branches/trunk-jpl-damage/src/m/model/mesh/planet/planettrimesh.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/mesh/planet/planettrimesh.m (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/m/model/mesh/planet/planettrimesh.m (revision 12946) @@ -0,0 +1,20 @@ +function md=planettrimesh(md,shape,radius,refinement) +%PLANETTRIMESH: build 2d shell mesh +% +% Usage: md=planettrimesh(md,shape,radius,refinement) +% + +results = sphere_tri(shape,refinement,radius); +md.x=results.vertices(:,1); +md.y=results.vertices(:,2); +md.z=results.vertices(:,3); +md.mesh.elements=results.faces; + +md.r=sqrt(md.x.^2+md.y.^2+md.z.^2); +md.theta=acos(md.z./md.r); +md.phi=atan2(md.y,md.x); + +md.mesh.numberofvertices=length(md.x); +md.mesh.numberofelements=size(md.mesh.elements,1); + +md.mesh.dimension=3; Index: /issm/branches/trunk-jpl-damage/src/m/model/mesh/planet/runme.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/mesh/planet/runme.m (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/m/model/mesh/planet/runme.m (revision 12946) @@ -0,0 +1,4 @@ +% 5 -> level of refinment +% 1000 -> radius +FV = sphere_tri('ico',5,1000); +patch('vertices',FV.vertices,'faces',FV.faces,'facecolor',[1 0 0],'edgecolor',[.2 .2 .6]); Index: /issm/branches/trunk-jpl-damage/src/m/model/mesh/planet/sphere_project.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/mesh/planet/sphere_project.m (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/m/model/mesh/planet/sphere_project.m (revision 12946) @@ -0,0 +1,66 @@ +function V = sphere_project(v,r,c) + +% sphere_project - project point X,Y,Z to the surface of sphere radius r +% +% V = sphere_project(v,r,c) +% +% Cartesian inputs: +% v is the vertex matrix, Nx3 (XYZ) +% r is the sphere radius, 1x1 (default 1) +% c is the sphere centroid, 1x3 (default 0,0,0) +% +% XYZ are converted to spherical coordinates and their radius is +% adjusted according to r, from c toward XYZ (defined with theta,phi) +% +% V is returned as Cartesian 3D coordinates +% + +% $Revision: 1.1 $ $Date: 2004/11/12 01:32:36 $ + +% Licence: GNU GPL, no implied or express warranties +% History: 06/2002, Darren.Weber_at_radiology.ucsf.edu, created +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +if ~exist('v','var'), + msg = sprintf('SPHERE_PROJECT: No input vertices (X,Y,Z)\n'); + error(msg); +end + +X = v(:,1); +Y = v(:,2); +Z = v(:,3); + +if ~exist('c','var'), + xo = 0; + yo = 0; + zo = 0; +else + xo = c(1); + yo = c(2); + zo = c(3); +end + +if ~exist('r','var'), r = 1; end + +% alternate method is to use unit vector of V +% [ n = 'magnitude(V)'; unitV = V ./ n; ] +% to change the radius, multiply the unitV +% by the radius required. This avoids the +% use of arctan functions, which have branches. + + +% Convert Cartesian X,Y,Z to spherical (radians) +theta = atan2( (Y-yo), (X-xo) ); +phi = atan2( sqrt( (X-xo).^2 + (Y-yo).^2 ), (Z-zo) ); +% do not calc: r = sqrt( (X-xo).^2 + (Y-yo).^2 + (Z-zo).^2); + +% Recalculate X,Y,Z for constant r, given theta & phi. +R = ones(size(phi)) * r; +x = R .* sin(phi) .* cos(theta); +y = R .* sin(phi) .* sin(theta); +z = R .* cos(phi); + +V = [x y z]; + +return Index: /issm/branches/trunk-jpl-damage/src/m/model/mesh/planet/sphere_tri.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/mesh/planet/sphere_tri.m (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/m/model/mesh/planet/sphere_tri.m (revision 12946) @@ -0,0 +1,204 @@ +function [FV] = sphere_tri(shape,maxlevel,r,winding) + +% sphere_tri - generate a triangle mesh approximating a sphere +% +% Usage: FV = sphere_tri(shape,Nrecurse,r,winding) +% +% shape is a string, either of the following: +% 'ico' starts with icosahedron (most even, default) +% 'oct' starts with octahedron +% 'tetra' starts with tetrahedron (least even) +% +% Nrecurse is int >= 0, setting the recursions (default 0) +% +% r is the radius of the sphere (default 1) +% +% winding is 0 for clockwise, 1 for counterclockwise (default 0). The +% matlab patch command gives outward surface normals for clockwise +% order of vertices in the faces (viewed from outside the surface). +% +% FV has fields FV.vertices and FV.faces. The vertices +% are listed in clockwise order in FV.faces, as viewed +% from the outside in a RHS coordinate system. +% +% The function uses recursive subdivision. The first +% approximation is an platonic solid, either an icosahedron, +% octahedron or a tetrahedron. Each level of refinement +% subdivides each triangle face by a factor of 4 (see also +% mesh_refine). At each refinement, the vertices are +% projected to the sphere surface (see sphere_project). +% +% A recursion level of 3 or 4 is a good sphere surface, if +% gouraud shading is used for rendering. +% +% The returned struct can be used in the patch command, eg: +% +% % create and plot, vertices: [2562x3] and faces: [5120x3] +% FV = sphere_tri('ico',4,1); +% lighting phong; shading interp; figure; +% patch('vertices',FV.vertices,'faces',FV.faces,... +% 'facecolor',[1 0 0],'edgecolor',[.2 .2 .6]); +% axis off; camlight infinite; camproj('perspective'); +% +% See also: mesh_refine, sphere_project +% + + + +% $Revision: 1.2 $ $Date: 2005/07/20 23:07:03 $ + +% Licence: GNU GPL, no implied or express warranties +% Jon Leech (leech @ cs.unc.edu) 3/24/89 +% icosahedral code added by Jim Buddenhagen (jb1556@daditz.sbc.com) 5/93 +% 06/2002, adapted from c to matlab by Darren.Weber_at_radiology.ucsf.edu +% 05/2004, reorder of the faces for the 'ico' surface so they are indeed +% clockwise! Now the surface normals are directed outward. Also reset the +% default recursions to zero, so we can get out just the platonic solids. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +eegversion = '$Revision: 1.2 $'; +fprintf('SPHERE_TRI [v %s]\n',eegversion(11:15)); tic + +if ~exist('shape','var') || isempty(shape), + shape = 'ico'; +end +fprintf('...creating sphere tesselation based on %s\n',shape); + +% default maximum subdivision level +if ~exist('maxlevel','var') || isempty(maxlevel) || maxlevel < 0, + maxlevel = 0; +end + +% default radius +if ~exist('r','var') || isempty(r), + r = 1; +end + +if ~exist('winding','var') || isempty(winding), + winding = 0; +end + + +% ----------------- +% define the starting shapes + +shape = lower(shape); + +switch shape, +case 'tetra', + + % Vertices of a tetrahedron + sqrt_3 = 0.5773502692; + + tetra.v = [ sqrt_3, sqrt_3, sqrt_3 ; % +X, +Y, +Z - PPP + -sqrt_3, -sqrt_3, sqrt_3 ; % -X, -Y, +Z - MMP + -sqrt_3, sqrt_3, -sqrt_3 ; % -X, +Y, -Z - MPM + sqrt_3, -sqrt_3, -sqrt_3 ]; % +X, -Y, -Z - PMM + + % Structure describing a tetrahedron + tetra.f = [ 1, 2, 3; + 1, 4, 2; + 3, 2, 4; + 4, 1, 3 ]; + + FV.vertices = tetra.v; + FV.faces = tetra.f; + +case 'oct', + + % Six equidistant points lying on the unit sphere + oct.v = [ 1, 0, 0 ; % X + -1, 0, 0 ; % -X + 0, 1, 0 ; % Y + 0, -1, 0 ; % -Y + 0, 0, 1 ; % Z + 0, 0, -1 ]; % -Z + + % Join vertices to create a unit octahedron + oct.f = [ 1 5 3 ; % X Z Y - First the top half + 3 5 2 ; % Y Z -X + 2 5 4 ; % -X Z -Y + 4 5 1 ; % -Y Z X + 1 3 6 ; % X Y -Z - Now the bottom half + 3 2 6 ; % Y Z -Z + 2 4 6 ; % -X Z -Z + 4 1 6 ]; % -Y Z -Z + + FV.vertices = oct.v; + FV.faces = oct.f; + +case 'ico', + + % Twelve vertices of icosahedron on unit sphere + tau = 0.8506508084; % t=(1+sqrt(5))/2, tau=t/sqrt(1+t^2) + one = 0.5257311121; % one=1/sqrt(1+t^2) , unit sphere + + ico.v( 1,:) = [ tau, one, 0 ]; % ZA + ico.v( 2,:) = [ -tau, one, 0 ]; % ZB + ico.v( 3,:) = [ -tau, -one, 0 ]; % ZC + ico.v( 4,:) = [ tau, -one, 0 ]; % ZD + ico.v( 5,:) = [ one, 0 , tau ]; % YA + ico.v( 6,:) = [ one, 0 , -tau ]; % YB + ico.v( 7,:) = [ -one, 0 , -tau ]; % YC + ico.v( 8,:) = [ -one, 0 , tau ]; % YD + ico.v( 9,:) = [ 0 , tau, one ]; % XA + ico.v(10,:) = [ 0 , -tau, one ]; % XB + ico.v(11,:) = [ 0 , -tau, -one ]; % XC + ico.v(12,:) = [ 0 , tau, -one ]; % XD + + % Structure for unit icosahedron + ico.f = [ 5, 8, 9 ; + 5, 10, 8 ; + 6, 12, 7 ; + 6, 7, 11 ; + 1, 4, 5 ; + 1, 6, 4 ; + 3, 2, 8 ; + 3, 7, 2 ; + 9, 12, 1 ; + 9, 2, 12 ; + 10, 4, 11 ; + 10, 11, 3 ; + 9, 1, 5 ; + 12, 6, 1 ; + 5, 4, 10 ; + 6, 11, 4 ; + 8, 2, 9 ; + 7, 12, 2 ; + 8, 10, 3 ; + 7, 3, 11 ]; + + FV.vertices = ico.v; + FV.faces = ico.f; +end + + +% ----------------- +% refine the starting shapes with subdivisions +if maxlevel, + + % Subdivide each starting triangle (maxlevel) times + for level = 1:maxlevel, + + % Subdivide each triangle and normalize the new points thus + % generated to lie on the surface of a sphere radius r. + FV = mesh_refine_tri4(FV); + FV.vertices = sphere_project(FV.vertices,r); + + % An alternative might be to define a min distance + % between vertices and recurse or use fminsearch + + end +end + +if winding, + fprintf('...returning counterclockwise vertex order (viewed from outside)\n'); + FV.faces = FV.faces(:,[1 3 2]); +else + fprintf('...returning clockwise vertex order (viewed from outside)\n'); +end + +t=toc; fprintf('...done (%6.2f sec)\n\n',t); + +return Index: /issm/branches/trunk-jpl-damage/src/m/model/parameterization/parameterize.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/parameterization/parameterize.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/model/parameterization/parameterize.py (revision 12946) @@ -6,14 +6,14 @@ """ PARAMETERIZE - parameterize a model - - from a parameter python file, start filling in all the model fields that were not - filled in by the mesh.py and mask.py model methods. - Warning: the parameter file must be able to be run in Python - - Usage: - md=parameterize(md,parametername) - - Example: - md=parameterize(md,'Square.par'); + + from a parameter python file, start filling in all the model fields that were not + filled in by the mesh.py and mask.py model methods. + Warning: the parameter file must be able to be run in Python + + Usage: + md=parameterize(md,parametername) + + Example: + md=parameterize(md,'Square.par'); """ Index: /issm/branches/trunk-jpl-damage/src/m/model/parseresultsfromdisk.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/parseresultsfromdisk.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/model/parseresultsfromdisk.py (revision 12946) @@ -7,7 +7,7 @@ """ PARSERESULTSFROMDISK - ... - - Usage: - results=parseresultsfromdisk(filename,iosplit) + + Usage: + results=parseresultsfromdisk(filename,iosplit) """ Index: /issm/branches/trunk-jpl-damage/src/m/model/petscversion.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/petscversion.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/model/petscversion.py (revision 12946) @@ -7,6 +7,6 @@ PETSCVERSION - recover petsc version number, inside config.h file - Usage: - PETSC_VERSION=petscversion(); + Usage: + PETSC_VERSION=petscversion(); """ Index: /issm/branches/trunk-jpl-damage/src/m/model/process_solve_options.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/process_solve_options.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/model/process_solve_options.py (revision 12946) @@ -1,3 +1,4 @@ import os +from EnumDefinitions import * def process_solve_options(options): @@ -5,8 +6,8 @@ DEFAULT_SOLVE_OPTIONS - set up default options for solve phase - Usage: - options=process_solve_options(options) + Usage: + options=process_solve_options(options) - See also: SOLVE + See also: SOLVE """ @@ -15,5 +16,5 @@ #solution_type: check on this option, error out otherwise solution_type=options.getfieldvalue('solution_type') - if solution_type in (DiagnosticSolutionEnum,PrognosticSolutionEnum,ThermalSolutionEnum,\ + if not solution_type in (DiagnosticSolutionEnum,PrognosticSolutionEnum,ThermalSolutionEnum,\ SteadystateSolutionEnum,TransientSolutionEnum,EnthalpySolutionEnum,\ BalancethicknessSolutionEnum,BedSlopeSolutionEnum,SurfaceSlopeSolutionEnum,HydrologySolutionEnum,FlaimSolutionEnum): @@ -42,5 +43,5 @@ outoptions['keep']=options.getfieldvalue('keep','y') outoptions['latsgn']=options.getfieldvalue('latsgn',0) - outoptions['cmap']=options.getfieldvalue('cmap',None) + outoptions['cmap']=options.getfieldvalue('cmap',[]) return outoptions Index: /issm/branches/trunk-jpl-damage/src/m/model/recover_areas.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/recover_areas.m (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/model/recover_areas.m (revision 12946) @@ -16,5 +16,4 @@ pattynflag =FlagElements(md,getfieldvalue(options,'pattyn','')); stokesflag =FlagElements(md,getfieldvalue(options,'stokes','')); - hutterflag =FlagElements(md,getfieldvalue(options,'hutter','')); filltype =getfieldvalue(options,'fill','none'); Index: /issm/branches/trunk-jpl-damage/src/m/model/recover_areas.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/recover_areas.py (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/m/model/recover_areas.py (revision 12946) @@ -0,0 +1,26 @@ +from pairoptions import * +from FlagElements import * + +def recover_areas(md,*args): + """ + RECOVER_AREAS - flag the element depending on the physical model that is assigned to them + + This routine is called by setelementstype, do not use + + Usage: + [hutterflag macayealflag pattynflag stokesflag filltype]=recover_areas(md,varargin); + """ + + #go through varargin, extract options and plug them into subtype options, by order of appearance + options=pairoptions(*args) +# options=deleteduplicates(options,1); + + #recover elements distribution + hutterflag =FlagElements(md,options.getfieldvalue('hutter','')) + macayealflag=FlagElements(md,options.getfieldvalue('macayeal','')) + pattynflag =FlagElements(md,options.getfieldvalue('pattyn','')) + stokesflag =FlagElements(md,options.getfieldvalue('stokes','')) + filltype =options.getfieldvalue('fill','none') + + return hutterflag,macayealflag,pattynflag,stokesflag,filltype + Index: /issm/branches/trunk-jpl-damage/src/m/model/setflowequation.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/setflowequation.m (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/model/setflowequation.m (revision 12946) @@ -33,5 +33,5 @@ [hutterflag macayealflag pattynflag stokesflag filltype]=recover_areas(md,varargin{:}); -%Flag the elements that has not been flagged as filltype +%Flag the elements that have not been flagged as filltype if strcmpi(filltype,'hutter'), hutterflag(find(~macayealflag & ~pattynflag))=1; @@ -58,5 +58,5 @@ if (md.mesh.dimension==2), if any(stokesflag | pattynflag) - error('setflowequation error message: stokes and pattyn elements no allowed in 2d mesh, extrude it first') + error('setflowequation error message: stokes and pattyn elements not allowed in 2d mesh, extrude it first') end end Index: /issm/branches/trunk-jpl-damage/src/m/model/setflowequation.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/setflowequation.py (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/m/model/setflowequation.py (revision 12946) @@ -0,0 +1,278 @@ +import numpy +from model import * +from pairoptions import * +from recover_areas import * +from MatlabFuncs import * + +def setflowequation(md,*args): + """ + SETELEMENTSTYPE - associate a solution type to each element + + This routine works like plotmodel: it works with an even number of inputs + 'hutter','macayeal','pattyn','stokes' and 'fill' are the possible options + that must be followed by the corresponding exp file or flags list + It can either be a domain file (argus type, .exp extension), or an array of element flags. + If user wants every element outside the domain to be + setflowequationd, add '~' to the name of the domain file (ex: '~Pattyn.exp'); + an empty string '' will be considered as an empty domain + a string 'all' will be considered as the entire domain + You can specify the type of coupling, 'penalties' or 'tiling', to use with the input 'coupling' + + Usage: + md=setflowequation(md,varargin) + + Example: + md=setflowequation(md,'pattyn','Pattyn.exp','macayeal',md.mask.elementonfloatingice,'fill','hutter'); + md=setflowequation(md,'pattyn','Pattyn.exp',fill','hutter','coupling','tiling'); + """ + + #some checks on list of arguments + if not isinstance(md,model) or not len(args): + raise TypeError("setflowequation error message") + + #Find_out what kind of coupling to use + options=pairoptions(*args) + coupling_method=options.getfieldvalue('coupling','tiling') + if not strcmpi(coupling_method,'tiling') and not strcmpi(coupling_method,'penalties'): + raise TypeError("coupling type can only be: tiling or penalties") + + hutterflag,macayealflag,pattynflag,stokesflag,filltype=recover_areas(md,*args) + + #Flag the elements that have not been flagged as filltype + if strcmpi(filltype,'hutter'): + hutterflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(macayealflag,pattynflag)))]=1 + elif strcmpi(filltype,'macayeal'): + macayealflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(pattynflag,stokesflag))))]=1 + elif strcmpi(filltype,'pattyn'): + pattynflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(macayealflag,stokesflag))))]=1 + + #check that each element has at least one flag + if not numpy.any(hutterflag+macayealflag+pattynflag+stokesflag): + raise TypeError("setflowequation error message: elements type not assigned, must be specified") + + #check that each element has only one flag + if numpy.any(hutterflag+macayealflag+pattynflag+stokesflag>1): + print "setflowequation warning message: some elements have several types, higher order type is used for them" + hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,macayealflag))]=0 + hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,pattynflag))]=0 + macayealflag[numpy.nonzero(numpy.logical_and(macayealflag,pattynflag))]=0 + + #Check that no pattyn or stokes for 2d mesh + if md.mesh.dimension==2: + if numpy.any(numpy.logical_or(stokesflag,pattynflag)): + raise TypeError("setflowequation error message: stokes and pattyn elements not allowed in 2d mesh, extrude it first") + + #Stokes can only be used alone for now: + if numpy.any(stokesflag) and numpy.any(hutterflag): + raise TypeError("setflowequation error message: stokes cannot be used with any other model for now, put stokes everywhere") + + #Initialize node fields + nodeonhutter=numpy.zeros(md.mesh.numberofvertices) + nodeonhutter[md.mesh.elements[numpy.nonzero(hutterflag),:].astype(int)-1]=1 + nodeonmacayeal=numpy.zeros(md.mesh.numberofvertices) + nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1 + nodeonpattyn=numpy.zeros(md.mesh.numberofvertices) + nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1 + nodeonstokes=numpy.zeros(md.mesh.numberofvertices) + noneflag=numpy.zeros(md.mesh.numberofelements) + + #First modify stokesflag to get rid of elements contrained everywhere (spc + border with pattyn or macayeal) + if any(stokesflag): +# fullspcnodes=double((~isnan(md.diagnostic.spcvx)+~isnan(md.diagnostic.spcvy)+~isnan(md.diagnostic.spcvz))==3 | (nodeonpattyn & nodeonstokes)); %find all the nodes on the boundary of the domain without icefront + fullspcnodes=numpy.logical_or(numpy.logical_not(numpy.isnan(md.diagnostic.spcvx))+ \ + numpy.logical_not(numpy.isnan(md.diagnostic.spcvy))+ \ + numpy.logical_not(numpy.isnan(md.diagnostic.spcvz))==3, \ + numpy.logical_and(nodeonpattyn,nodeonstokes)).astype(int) #find all the nodes on the boundary of the domain without icefront +# fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6); %find all the nodes on the boundary of the domain without icefront + fullspcelems=(numpy.sum(fullspcnodes[md.mesh.elements.astype(int)-1],axis=1)==6).astype(int) #find all the nodes on the boundary of the domain without icefront + stokesflag[numpy.nonzero(fullspcelems)]=0 + nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1 + + #Then complete with NoneApproximation or the other model used if there is no stokes + if any(stokesflag): + if any(pattynflag): #fill with pattyn + pattynflag[numpy.logical_not(stokesflag)]=1 + nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1 + elif any(macayealflag): #fill with macayeal + macayealflag[numpy.logical_not(stokesflag)]=1 + nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1 + else: #fill with none + noneflag[numpy.nonzero(numpy.logical_not(stokesflag))]=1 + + #Now take care of the coupling between MacAyeal and Pattyn + md.diagnostic.vertex_pairing=numpy.array([]) + nodeonmacayealpattyn=numpy.zeros(md.mesh.numberofvertices) + nodeonpattynstokes=numpy.zeros(md.mesh.numberofvertices) + nodeonmacayealstokes=numpy.zeros(md.mesh.numberofvertices) + macayealpattynflag=numpy.zeros(md.mesh.numberofelements) + macayealstokesflag=numpy.zeros(md.mesh.numberofelements) + pattynstokesflag=numpy.zeros(md.mesh.numberofelements) + if strcmpi(coupling_method,'penalties'): + #Create the border nodes between Pattyn and MacAyeal and extrude them + numnodes2d=md.mesh.numberofvertices2d + numlayers=md.mesh.numberoflayers + bordernodes2d=numpy.nonzero(numpy.logical_and(nodeonpattyn[1:numnodes2d],nodeonmacayeal[1:numnodes2d])) #Nodes connected to two different types of elements + + #initialize and fill in penalties structure + if numpy.all(numpy.logical_not(numpy.isnan(bordernodes2d))): + penalties=numpy.zeros((0,2)) + for i in xrange(1,numlayers): + penalties=numpy.concatenate((penalties,numpy.concatenate((bordernodes2d,bordernodes2d+md.mesh.numberofvertices2d*(i)),axis=1)),axis=0) + md.diagnostic.vertex_pairing=penalties + + elif strcmpi(coupling_method,'tiling'): + if numpy.any(macayealflag) and numpy.any(pattynflag): #coupling macayeal pattyn + #Find node at the border + nodeonmacayealpattyn[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonpattyn))]=1 + #Macayeal elements in contact with this layer become MacAyealPattyn elements + matrixelements=ismember(md.mesh.elements,numpy.nonzero(nodeonmacayealpattyn)) + commonelements=numpy.sum(matrixelements,axis=1)!=0 + commonelements[numpy.nonzero(pattynflag)]=0 #only one layer: the elements previously in macayeal + macayealflag[numpy.nonzero(commonelements)]=0 #these elements are now macayealpattynelements + macayealpattynflag[numpy.nonzero(commonelements)]=1 + nodeonmacayeal[:]=0 + nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1 + + #rule out elements that don't touch the 2 boundaries + pos=numpy.nonzero(macayealpattynflag) + elist=numpy.zeros(len(pos)) + elist = elist + numpy.any(numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:].astype(int)-1],axis=1),axis=1) + elist = elist - numpy.any(numpy.sum(nodeonpattyn[md.mesh.elements[pos,:].astype(int)-1] ,axis=1),axis=1) + pos1=[i for i,item in enumerate(elist) if item==1] + macayealflag[pos[pos1]]=1 + macayealpattynflag[pos[pos1]]=0 + pos2=[i for i,item in enumerate(elist) if item==-1] + pattynflag[pos[pos2]]=1 + macayealpattynflag[pos[pos2]]=0 + + #Recompute nodes associated to these elements + nodeonmacayeal[:]=0 + nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1 + nodeonpattyn[:]=0 + nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1 + nodeonmacayealpattyn[:]=0 + nodeonmacayealpattyn[md.mesh.elements[numpy.nonzero(macayealpattynflag),:].astype(int)-1]=1 + + elif numpy.any(pattynflag) and numpy.any(stokesflag): #coupling pattyn stokes + #Find node at the border + nodeonpattynstokes[numpy.nonzero(numpy.logical_and(nodeonpattyn,nodeonstokes))]=1 + #Stokes elements in contact with this layer become PattynStokes elements + matrixelements=ismember(md.mesh.elements,numpy.nonzero(nodeonpattynstokes)) + commonelements=numpy.sum(matrixelements,axis=1)!=0 + commonelements[numpy.nonzero(pattynflag)]=0 #only one layer: the elements previously in macayeal + stokesflag[numpy.nonzero(commonelements)]=0 #these elements are now macayealpattynelements + pattynstokesflag[numpy.nonzero(commonelements)]=1 + nodeonstokes=numpy.zeros(md.mesh.numberofvertices) + nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1 + + #rule out elements that don't touch the 2 boundaries + pos=numpy.nonzero(pattynstokesflag) + elist=numpy.zeros(len(pos)) + elist = elist + numpy.any(numpy.sum(nodeonstokes[md.mesh.elements[pos,:].astype(int)-1],axis=1),axis=1) + elist = elist - numpy.any(numpy.sum(nodeonpattyn[md.mesh.elements[pos,:].astype(int)-1],axis=1),axis=1) + pos1=[i for i,item in enumerate(elist) if item==1] + stokesflag[pos[pos1]]=1 + pattynstokesflag[pos[pos1]]=0 + pos2=[i for i,item in enumerate(elist) if item==-1] + pattynflag[pos[pos2]]=1 + pattynstokesflag[pos[pos2]]=0 + + #Recompute nodes associated to these elements + nodeonstokes[:]=0 + nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1 + nodeonpattyn[:]=0 + nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1 + nodeonpattynstokes[:]=0 + nodeonpattynstokes[md.mesh.elements[numpy.nonzero(pattynstokesflag),:].astype(int)-1]=1 + + elif numpy.any(stokesflag) and numpy.any(macayealflag): + #Find node at the border + nodeonmacayealstokes[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonstokes))]=1 + #Stokes elements in contact with this layer become MacAyealStokes elements + matrixelements=ismember(md.mesh.elements,numpy.nonzero(nodeonmacayealstokes)) + commonelements=numpy.sum(matrixelements,axis=1)!=0 + commonelements[numpy.nonzero(macayealflag)]=0 #only one layer: the elements previously in macayeal + stokesflag[numpy.nonzero(commonelements)]=0 #these elements are now macayealmacayealelements + macayealstokesflag[numpy.nonzero(commonelements)]=1 + nodeonstokes=numpy.zeros(md.mesh.numberofvertices) + nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1 + + #rule out elements that don't touch the 2 boundaries + pos=numpy.nonzero(macayealstokesflag) + elist=numpy.zeros(len(pos)) + elist = elist + numpy.any(numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:].astype(int)-1],axis=1),axis=1) + elist = elist - numpy.any(numpy.sum(nodeonstokes[md.mesh.elements[pos,:].astype(int)-1] ,axis=1),axis=1) + pos1=[i for i,item in enumerate(elist) if item==1] + macayealflag[pos[pos1]]=1 + macayealstokesflag[pos[pos1]]=0 + pos2=[i for i,item in enumerate(elist) if item==-1] + stokesflag[pos[pos2]]=1 + macayealstokesflag[pos[pos2]]=0 + + #Recompute nodes associated to these elements + nodeonmacayeal[:]=0 + nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1 + nodeonstokes[:]=0 + nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1 + nodeonmacayealstokes[:]=0 + nodeonmacayealstokes[md.mesh.elements[numpy.nonzero(macayealstokesflag),:].astype(int)-1]=1 + + elif numpy.any(stokesflag) and numpy.any(hutterflag): + raise TypeError("type of coupling not supported yet") + + #Create MacaAyealPattynApproximation where needed + md.flowequation.element_equation=numpy.zeros(md.mesh.numberofelements) + md.flowequation.element_equation[numpy.nonzero(noneflag)]=0 + md.flowequation.element_equation[numpy.nonzero(hutterflag)]=1 + md.flowequation.element_equation[numpy.nonzero(macayealflag)]=2 + md.flowequation.element_equation[numpy.nonzero(pattynflag)]=3 + md.flowequation.element_equation[numpy.nonzero(stokesflag)]=4 + md.flowequation.element_equation[numpy.nonzero(macayealpattynflag)]=5 + md.flowequation.element_equation[numpy.nonzero(macayealstokesflag)]=6 + md.flowequation.element_equation[numpy.nonzero(pattynstokesflag)]=7 + + #border + md.flowequation.borderpattyn=nodeonpattyn + md.flowequation.bordermacayeal=nodeonmacayeal + md.flowequation.borderstokes=nodeonstokes + + #Create vertices_type + md.flowequation.vertex_equation=numpy.zeros(md.mesh.numberofvertices) + pos=numpy.nonzero(nodeonhutter) + md.flowequation.vertex_equation[pos]=1 + pos=numpy.nonzero(nodeonmacayeal) + md.flowequation.vertex_equation[pos]=2 + pos=numpy.nonzero(nodeonpattyn) + md.flowequation.vertex_equation[pos]=3 + pos=numpy.nonzero(nodeonhutter) + md.flowequation.vertex_equation[pos]=1 + pos=numpy.nonzero(nodeonmacayealpattyn) + md.flowequation.vertex_equation[pos]=5 + pos=numpy.nonzero(nodeonstokes) + md.flowequation.vertex_equation[pos]=4 + if numpy.any(stokesflag): + pos=numpy.nonzero(numpy.logical_not(nodeonstokes)) + if not (numpy.any(pattynflag) or numpy.any(macayealflag)): + md.flowequation.vertex_equation[pos]=0 + pos=numpy.nonzero(nodeonpattynstokes) + md.flowequation.vertex_equation[pos]=7 + pos=numpy.nonzero(nodeonmacayealstokes) + md.flowequation.vertex_equation[pos]=6 + + #figure out solution types + md.flowequation.ishutter=float(numpy.any(md.flowequation.element_equation==1)) + md.flowequation.ismacayealpattyn=float(numpy.any(numpy.logical_or(md.flowequation.element_equation==2,md.flowequation.element_equation==3))) + md.flowequation.isstokes=float(numpy.any(md.flowequation.element_equation==4)) + + return md + + #Check that tiling can work: + if numpy.any(md.flowequation.bordermacayeal) and numpy.any(md.flowequation.borderpattyn) and numpy.any(md.flowequation.borderpattyn + md.flowequation.bordermacayeal !=1): + raise TypeError("error coupling domain too irregular") + if numpy.any(md.flowequation.bordermacayeal) and numpy.any(md.flowequation.borderstokes) and numpy.any(md.flowequation.borderstokes + md.flowequation.bordermacayeal !=1): + raise TypeError("error coupling domain too irregular") + if numpy.any(md.flowequation.borderstokes) and numpy.any(md.flowequation.borderpattyn) and numpy.any(md.flowequation.borderpattyn + md.flowequation.borderstokes !=1): + raise TypeError("error coupling domain too irregular") + + return md + Index: /issm/branches/trunk-jpl-damage/src/m/model/solve.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/solve.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/model/solve.py (revision 12946) @@ -2,4 +2,8 @@ import os import shutil +from pairoptions import * +from process_solve_options import * +from EnumDefinitions import * +from ismodelselfconsistent import * def solve(md,solutionenum,*args): @@ -7,25 +11,25 @@ SOLVE - apply solution sequence for this model - Usage: - md=solve(md,solutionenum,varargin) - where varargin is a list of paired arguments of string OR enums + Usage: + md=solve(md,solutionenum,varargin) + where varargin is a list of paired arguments of string OR enums - solution types available comprise: - - DiagnosticSolutionEnum - - PrognosticSolutionEnum - - ThermalSolutionEnum - - SteadystateSolutionEnum - - TransientSolutionEnum... - - BalancethicknessSolutionEnum - - BedSlopeSolutionEnum - - SurfaceSlopeSolutionEnum - - HydrologySolutionEnum - - FlaimSolutionEnum + solution types available comprise: + - DiagnosticSolutionEnum + - PrognosticSolutionEnum + - ThermalSolutionEnum + - SteadystateSolutionEnum + - TransientSolutionEnum... + - BalancethicknessSolutionEnum + - BedSlopeSolutionEnum + - SurfaceSlopeSolutionEnum + - HydrologySolutionEnum + - FlaimSolutionEnum extra options: - - loadonly : does not solve. only load results + - loadonly : does not solve. only load results - Examples: - md=solve(md,DiagnosticSolutionEnum); + Examples: + md=solve(md,DiagnosticSolutionEnum); """ Index: /issm/branches/trunk-jpl-damage/src/m/model/solvers/asmoptions.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/solvers/asmoptions.py (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/m/model/solvers/asmoptions.py (revision 12946) @@ -0,0 +1,31 @@ +#module imports {{{ +from pairoptions import * +from petscversion import * +#}}} +def asmoptions(*args): + #ASMOPTIONS - return ASM petsc options + # + # Usage: + # options=asmoptions; + + #retrieve options provided in varargin + arguments=pairoptions(*args) + + options=[['mat_type','aij'],['ksp_type','gmres'],['pc_type','asm'],['sub_pc_type','lu'],['pc_asm_overlap',3],['ksp_max_it',100],['ksp_rtol',1e-30]]; + + #now, go through our arguments, and write over default options. + for i in range(len(arguments.list)): + arg1=arguments.list[i][0] + arg2=arguments.list[i][1] + found=0; + for j in range(len(options)): + joption=options[j][0] + if joption==arg1: + joption[1]=arg2; + options[j]=joption; + found=1; + break + if not found: + #this option did not exist, add it: + options.append([arg1,arg2]) + return options Index: /issm/branches/trunk-jpl-damage/src/m/model/solvers/iluasmoptions.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/solvers/iluasmoptions.m (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/model/solvers/iluasmoptions.m (revision 12946) @@ -10,8 +10,8 @@ %default iluasm options -iluasm.mat_type=getfieldvalue(options,'mat_type''','aij'); -iluasm.ksp_type=getfieldvalue(options,'ksp_type''','gmres'); -iluasm.pc_type=getfieldvalue(options,'pc_type''','asm'); -iluasm.sub_pc_type=getfieldvalue(options,'sub_pc_type''','ilu'); +iluasm.mat_type=getfieldvalue(options,'mat_type','aij'); +iluasm.ksp_type=getfieldvalue(options,'ksp_type','gmres'); +iluasm.pc_type=getfieldvalue(options,'pc_type','asm'); +iluasm.sub_pc_type=getfieldvalue(options,'sub_pc_type','ilu'); iluasm.pc_asm_overlap=getfieldvalue(options,'pc_asm_overlap',5); iluasm.ksp_max_it=getfieldvalue(options,'ksp_max_it',100); Index: /issm/branches/trunk-jpl-damage/src/m/model/solvers/iluasmoptions.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/solvers/iluasmoptions.py (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/m/model/solvers/iluasmoptions.py (revision 12946) @@ -0,0 +1,25 @@ +from pairoptions import * + +def iluasmoptions(*args): + """ + ILUASMOPTIONS - + + Usage: + options=iluasmoptions; + """ + + #retrieve options provided in varargin + options=pairoptions(*args) + iluasm={} + + #default iluasm options + iluasm['mat_type']=options.getfieldvalue('mat_type','aij') + iluasm['ksp_type']=options.getfieldvalue('ksp_type','gmres') + iluasm['pc_type']=options.getfieldvalue('pc_type','asm') + iluasm['sub_pc_type']=options.getfieldvalue('sub_pc_type','ilu') + iluasm['pc_asm_overlap']=options.getfieldvalue('pc_asm_overlap',5) + iluasm['ksp_max_it']=options.getfieldvalue('ksp_max_it',100) + iluasm['ksp_rtol']=options.getfieldvalue('ksp_rtol',1e-15) + + return iluasm + Index: /issm/branches/trunk-jpl-damage/src/m/model/solvers/jacobiasmoptions.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/solvers/jacobiasmoptions.py (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/m/model/solvers/jacobiasmoptions.py (revision 12946) @@ -0,0 +1,31 @@ +#module imports {{{ +from pairoptions import * +from petscversion import * +#}}} +def jacobiasmoptions(*args): + #ASMOPTIONS - return Additive Shwartz Method with Jacobi preconditioner petsc options + # + # Usage: + # options=jacobiasmoptions; + + #retrieve options provided in varargin + arguments=pairoptions(*args) + + options=[['mat_type','aij'],['ksp_type','gmres'],['pc_type','asm'],['sub_pc_type','jacobi'],['pc_asm_overlap',3],['ksp_max_it',100],['ksp_rtol',1e-15]]; + + #now, go through our arguments, and write over default options. + for i in range(len(arguments.list)): + arg1=arguments.list[i][0] + arg2=arguments.list[i][1] + found=0; + for j in range(len(options)): + joption=options[j][0] + if joption==arg1: + joption[1]=arg2; + options[j]=joption; + found=1; + break + if not found: + #this option did not exist, add it: + options.append([arg1,arg2]) + return options Index: /issm/branches/trunk-jpl-damage/src/m/model/solvers/jacobicgoptions.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/solvers/jacobicgoptions.py (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/m/model/solvers/jacobicgoptions.py (revision 12946) @@ -0,0 +1,31 @@ +#module imports {{{ +from pairoptions import * +from petscversion import * +#}}} +def jacobicgoptions(*args): + #ASMOPTIONS - return Additive Shwartz Method with Jacobi preconditioner petsc options + # + # Usage: + # options=jacobicgoptions; + + #retrieve options provided in varargin + arguments=pairoptions(*args) + + options=[['mat_type','aij'],['ksp_type','cg'],['ksp_max_it',100],['ksp_rtol',1e-15]]; + + #now, go through our arguments, and write over default options. + for i in range(len(arguments.list)): + arg1=arguments.list[i][0] + arg2=arguments.list[i][1] + found=0; + for j in range(len(options)): + joption=options[j][0] + if joption==arg1: + joption[1]=arg2; + options[j]=joption; + found=1; + break + if not found: + #this option did not exist, add it: + options.append([arg1,arg2]) + return options Index: /issm/branches/trunk-jpl-damage/src/m/model/solvers/matlaboptions.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/solvers/matlaboptions.py (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/m/model/solvers/matlaboptions.py (revision 12946) @@ -0,0 +1,31 @@ +#module imports {{{ +from pairoptions import * +from petscversion import * +#}}} +def matlaboptions(*args): + #MATLABOPTIONS - return Matlab petsc options + # + # Usage: + # options=matlaboptions; + + #retrieve options provided in varargin + arguments=pairoptions(*args) + + options=[['ksp_type','matlab']]; + + #now, go through our arguments, and write over default options. + for i in range(len(arguments.list)): + arg1=arguments.list[i][0] + arg2=arguments.list[i][1] + found=0; + for j in range(len(options)): + joption=options[j][0] + if joption==arg1: + joption[1]=arg2; + options[j]=joption; + found=1; + break + if not found: + #this option did not exist, add it: + options.append([arg1,arg2]) + return options Index: /issm/branches/trunk-jpl-damage/src/m/model/solvers/mumpsoptions.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/solvers/mumpsoptions.py (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/m/model/solvers/mumpsoptions.py (revision 12946) @@ -0,0 +1,33 @@ +from pairoptions import * +from petscversion import * + +def mumpsoptions(*args): + """ + MUMPSOPTIONS - return MUMPS direct solver petsc options + + Usage: + options=mumpsoptions; + """ + + #retrieve options provided in varargin + options=pairoptions(*args) + mumps={} + + #default mumps options + PETSC_VERSION=petscversion() + if PETSC_VERSION==2: + mumps['mat_type']=options.getfieldvalue('mat_type','aijmumps') + mumps['ksp_type']=options.getfieldvalue('ksp_type','preonly') + mumps['pc_type']=options.getfieldvalue('pc_type','lu') + mumps['mat_mumps_icntl_14']=options.getfieldvalue('mat_mumps_icntl_14',120) + mumps['pc_factor_shift_positive_definite']=options.getfieldvalue('pc_factor_shift_positive_definite','true') + if PETSC_VERSION==3: + mumps['mat_type']=options.getfieldvalue('mat_type','mpiaij') + mumps['ksp_type']=options.getfieldvalue('ksp_type','preonly') + mumps['pc_type']=options.getfieldvalue('pc_type','lu') + mumps['pc_factor_mat_solver_package']=options.getfieldvalue('pc_factor_mat_solver_package','mumps') + mumps['mat_mumps_icntl_14']=options.getfieldvalue('mat_mumps_icntl_14',120) + mumps['pc_factor_shift_positive_definite']=options.getfieldvalue('pc_factor_shift_positive_definite','true') + + return mumps + Index: /issm/branches/trunk-jpl-damage/src/m/model/solvers/soroptions.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/solvers/soroptions.py (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/m/model/solvers/soroptions.py (revision 12946) @@ -0,0 +1,31 @@ +#module imports {{{ +from pairoptions import * +from petscversion import * +#}}} +def soroptions(*args): + #SOROPTIONS - return Relaxation Solver petsc options + # + # Usage: + # options=soroptions; + + #retrieve options provided in varargin + arguments=pairoptions(*args) + + options=[['mat_type','aij'],['ksp_type','cg'],['pc_type','sor'],['pc_sor_omega',1.1],['pc_sor_its',2]]; + + #now, go through our arguments, and write over default options. + for i in range(len(arguments.list)): + arg1=arguments.list[i][0] + arg2=arguments.list[i][1] + found=0; + for j in range(len(options)): + joption=options[j][0] + if joption==arg1: + joption[1]=arg2; + options[j]=joption; + found=1; + break + if not found: + #this option did not exist, add it: + options.append([arg1,arg2]) + return options Index: /issm/branches/trunk-jpl-damage/src/m/model/solvers/stokesoptions.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/solvers/stokesoptions.py (revision 12946) +++ /issm/branches/trunk-jpl-damage/src/m/model/solvers/stokesoptions.py (revision 12946) @@ -0,0 +1,40 @@ +#module imports {{{ +from pairoptions import * +from petscversion import * +#}}} +def stokesoptions(*args): + #STOKESOPTIONS - return STOKES multi-physics solver petsc options + # + # Usage: + # options=stokesoptions; + + #retrieve options provided in varargin + arguments=pairoptions(*args) + + + #default stokes options + PETSC_VERSION=petscversion() + + if PETSC_VERSION==2: + raise RuntimeError('stokesoptions error message: multi-physics options not supported in Petsc 2') + if PETSC_VERSION==3: + options=[['mat_type','mpiaij'],['ksp_max_it',1000],['ksp_type','gmres'],['pc_type','fieldsplit'],['pc_field_split_type','schur'],\ + ['fieldsplit_0_pc_type','hypre'],['fieldsplit_0_ksp_type','gmres'],['fieldsplit_0_pc_hypre_type','boomerang'],\ + ['fieldsplit_1_pc_type','jacobi'],['fieldsplit_1_ksp_type','preonly'],['issm_option_solver','stokes']] + + #now, go through our arguments, and write over default options. + for i in range(len(arguments.list)): + arg1=arguments.list[i][0] + arg2=arguments.list[i][1] + found=0; + for j in range(len(options)): + joption=options[j][0] + if joption==arg1: + joption[1]=arg2; + options[j]=joption; + found=1; + break + if not found: + #this option did not exist, add it: + options.append([arg1,arg2]) + return options Index: /issm/branches/trunk-jpl-damage/src/m/model/waitonlock.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/waitonlock.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/model/waitonlock.py (revision 12946) @@ -8,9 +8,9 @@ WAITONLOCK - wait for a file - This routine will return when a file named 'filename' is written to disk. - If the time limit given in input is exceeded, return 0 + This routine will return when a file named 'filename' is written to disk. + If the time limit given in input is exceeded, return 0 - Usage: - flag=waitonlock(md,executionpath) + Usage: + flag=waitonlock(md,executionpath) """ Index: /issm/branches/trunk-jpl-damage/src/m/utils/BC/SetIceSheetBC.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/utils/BC/SetIceSheetBC.m (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/utils/BC/SetIceSheetBC.m (revision 12946) @@ -34,9 +34,9 @@ %Create zeros basal melting rate and surface mass balance if not specified -if isnan(md.surfaceforcings.precipitation), +if (isnan(md.surfaceforcings.precipitation) & (md.surfaceforcings.ispdd==1) & (md.surfaceforcings.ispdd==0)), md.surfaceforcings.precipitation=zeros(md.mesh.numberofvertices,1); disp(' no surfaceforcings.precipitation specified: values set as zero'); end -if isnan(md.surfaceforcings.mass_balance), +if isnan(md.surfaceforcings.mass_balance) & (md.surfaceforcings.ispdd==0), md.surfaceforcings.mass_balance=zeros(md.mesh.numberofvertices,1); disp(' no surfaceforcings.mass_balance specified: values set as zero'); Index: /issm/branches/trunk-jpl-damage/src/m/utils/BC/SetIceShelfBC.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/utils/BC/SetIceShelfBC.m (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/utils/BC/SetIceShelfBC.m (revision 12946) @@ -65,11 +65,11 @@ %Create zeros basalforcings and surfaceforcings -if isnan(md.surfaceforcings.precipitation), +if (isnan(md.surfaceforcings.precipitation) & (md.surfaceforcings.ispdd==1) & (md.surfaceforcings.ispdd==0)), md.surfaceforcings.precipitation=zeros(md.mesh.numberofvertices,1); disp(' no surfaceforcings.precipitation specified: values set as zero'); end -if isnan(md.surfaceforcings.mass_balance), - md.surfaceforcings.mass_balance=zeros(md.mesh.numberofvertices,1); - disp(' no surfaceforcings.mass_balance specified: values set as zero'); +if isnan(md.surfaceforcings.mass_balance) & (md.surfaceforcings.ispdd==0), + md.surfaceforcings.mass_balance=zeros(md.mesh.numberofvertices,1); + disp(' no surfaceforcings.mass_balance specified: values set as zero'); end if isnan(md.basalforcings.melting_rate), Index: /issm/branches/trunk-jpl-damage/src/m/utils/BC/SetIceShelfBC.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/utils/BC/SetIceShelfBC.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/utils/BC/SetIceShelfBC.py (revision 12946) @@ -5,17 +5,17 @@ """ SETICESHELFBC - Create the boundary conditions for diagnostic and thermal models for a Ice Shelf with Ice Front - - Neumann BC are used on the ice front (an ANRGUS contour around the ice front - must be given in input) - Dirichlet BC are used elsewhere for diagnostic - - Usage: - md=SetIceShelfBC(md,varargin) - - Example: - md=SetIceShelfBC(md); - md=SetIceShelfBC(md,'Front.exp'); - - See also: SETICESHEETBC, SETMARINEICESHEETBC + + Neumann BC are used on the ice front (an ANRGUS contour around the ice front + must be given in input) + Dirichlet BC are used elsewhere for diagnostic + + Usage: + md=SetIceShelfBC(md,varargin) + + Example: + md=SetIceShelfBC(md); + md=SetIceShelfBC(md,'Front.exp'); + + See also: SETICESHEETBC, SETMARINEICESHEETBC """ Index: /issm/branches/trunk-jpl-damage/src/m/utils/BC/SetMarineIceSheetBC.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/utils/BC/SetMarineIceSheetBC.m (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/utils/BC/SetMarineIceSheetBC.m (revision 12946) @@ -76,9 +76,9 @@ %Create zeros basalforcings and surfaceforcings -if isnan(md.surfaceforcings.precipitation), +if (isnan(md.surfaceforcings.precipitation)& (md.surfaceforcings.ispdd==1) & (md.surfaceforcings.ispdd==0)), md.surfaceforcings.precipitation=zeros(md.mesh.numberofvertices,1); disp(' no surfaceforcings.precipitation specified: values set as zero'); end -if isnan(md.surfaceforcings.mass_balance), +if isnan(md.surfaceforcings.mass_balance) & (md.surfaceforcings.ispdd==0), md.surfaceforcings.mass_balance=zeros(md.mesh.numberofvertices,1); disp(' no surfaceforcings.mass_balance specified: values set as zero'); Index: /issm/branches/trunk-jpl-damage/src/m/utils/Cluster/parallelrange.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/utils/Cluster/parallelrange.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/utils/Cluster/parallelrange.py (revision 12946) @@ -4,6 +4,6 @@ PARALLELRANGE - from a rank, and a number of processors, figure out a range, for parallel tasks. - Usage: - i1,i2=parallelrange(rank,numprocs,globalsize) + Usage: + i1,i2=parallelrange(rank,numprocs,globalsize) """ Index: /issm/branches/trunk-jpl-damage/src/m/utils/Geometry/FlagElements.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/utils/Geometry/FlagElements.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/utils/Geometry/FlagElements.py (revision 12946) @@ -1,54 +1,62 @@ -from numpy import * +import numpy +import os +#from basinzoom import * +#from ContourToMesh import * +from MatlabFuncs import * + def FlagElements(md,region): -#FLAGELEMENTS - flag the elements in an region -# -# The region can be given with an exp file, a list of elements. -# -# Usage: -# flag=FlagElements(md,region); -# -# Example: -# flag=FlagElements(md,'all'); -# flag=FlagElements(md,''); -# flag=FlagElements(md,'Domain.exp'); -# flag=FlagElements(md,'~Domain.exp'); -# flag=FlagElements(md,md.mask.elementongroundedice); + """ + FLAGELEMENTS - flag the elements in an region - if isinstance(region,basestring): - if not(region): - flag=zeros(md.mesh.numberofelements,'bool') - invert=0; - elif region=='all': - flag=ones(md.mesh.numberofelements,'bool') - invert=0; + The region can be given with an exp file, a list of elements. + + Usage: + flag=FlagElements(md,region); + + Example: + flag=FlagElements(md,'all'); + flag=FlagElements(md,''); + flag=FlagElements(md,'Domain.exp'); + flag=FlagElements(md,'~Domain.exp'); + flag=FlagElements(md,md.mask.elementongroundedice); + """ + + if isinstance(region,str): + if not region: + flag=numpy.zeros(md.mesh.numberofelements,'bool') + invert=0 + elif strcmpi(region,'all'): + flag=numpy.ones(md.mesh.numberofelements,'bool') + invert=0 else: #make sure that we actually don't want the elements outside the domain outline! - if region[0]=='~': + if strcmpi(region[0],'~'): region=region[1:] - invert=1; + invert=1 else: - invert=0; - + invert=0 + #does the region domain outline exist or do we have to look for xlim,ylim in basinzoom? - if not os.path.isfile(region): - [xlim,ylim]=basinzoom('basin',region); - flag_nodes=double(md.mesh.xxlim(1) & md.mesh.yylim(1)); - flag=prod(flag_nodes(md.mesh.elements),2); + if not os.path.exists(region): + if len(region)>3 and not strcmp(region[-4:],'.exp'): + raise IOError("Error: File 'region' not found!" % region) + xlim,ylim=basinzoom('basin',region) + flag_nodes=numpy.logical_and(numpy.logical_and(md.mesh.xxlim[0]),numpy.logical_and(md.mesh.yylim[0])).astype(float) + flag=numpy.prod(flag_nodes[md.mesh.elements],axis=1) else: #ok, flag elements - flag=ContourToMesh(md.mesh.elements[:,0:3],md.mesh.x,md.mesh.y,region,'element',1); - + flag=ContourToMesh(md.mesh.elements[:,0:3],md.mesh.x,md.mesh.y,region,'element',1) + if invert: - flag=~flag; - - elif isinstance(region,nparray): - if len(region)!=md.mesh.numberofelements: - print FlagElements.__doc__ - print 'Flaglist for region must be of same size as number of elements in model' - return [] - flag=region; + flag=numpy.logical_not(flag) + + elif isinstance(region,numpy.nparray) or isinstance(region,bool): + if not numpy.size(region,0)==md.mesh.numberofelements: + raise TypeError("Flaglist for region must be of same size as number of elements in model.") + flag=region + else: - print 'Invalid region option' - return [] + raise TypeError("Invalid region option") - return flag; + return flag + Index: /issm/branches/trunk-jpl-damage/src/m/utils/OS/ismumps.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/utils/OS/ismumps.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/utils/OS/ismumps.py (revision 12946) @@ -7,6 +7,6 @@ ISMUMPS - figure out if MUMPS package was compiled with ISSM - Usage: - flag=ismumps(); + Usage: + flag=ismumps(); """ Index: /issm/branches/trunk-jpl-damage/src/m/utils/OS/ispetsc.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/utils/OS/ispetsc.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/utils/OS/ispetsc.py (revision 12946) @@ -7,6 +7,6 @@ ISPETSC - figure out if PETSC package was compiled with ISSM - Usage: - flag=ispetsc(); + Usage: + flag=ispetsc(); """ Index: /issm/branches/trunk-jpl-damage/src/m/utils/OS/issmscpin.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/utils/OS/issmscpin.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/utils/OS/issmscpin.py (revision 12946) @@ -10,7 +10,5 @@ SCPIN get packages from host, using scp on unix, and pscp on windows - usage: scpin(host,packages,path) - - + usage: scpin(host,packages,path) """ Index: /issm/branches/trunk-jpl-damage/src/m/utils/OS/issmscpout.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/utils/OS/issmscpout.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/utils/OS/issmscpout.py (revision 12946) @@ -10,7 +10,5 @@ SCPOUT send packages to a host, using scp on unix, and pscp on windows - usage: scpout(host,path,packages) - - + usage: scpout(host,path,packages) """ Index: /issm/branches/trunk-jpl-damage/src/m/utils/OS/issmssh.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/utils/OS/issmssh.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/utils/OS/issmssh.py (revision 12946) @@ -9,6 +9,6 @@ ISSMSSH - wrapper for OS independent ssh command. - usage: - issmssh(host,command) + usage: + issmssh(host,command) """ Index: /issm/branches/trunk-jpl-damage/src/m/utils/Shell/issmdir.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/utils/Shell/issmdir.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/utils/Shell/issmdir.py (revision 12946) @@ -7,6 +7,6 @@ ISSMDIR - Get ISSM_DIR environment variable - Usage: - ISSM_DIR=issmdir() + Usage: + ISSM_DIR=issmdir() """ Index: /issm/branches/trunk-jpl-damage/src/m/utils/consistency/checkfield.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/utils/consistency/checkfield.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/src/m/utils/consistency/checkfield.py (revision 12946) @@ -8,27 +8,27 @@ CHECKFIELD - check field consistency - Used to check model consistency. - Available options: - - NaN: 1 if check that there is no NaN - - size: [lines cols], NaN for non checked dimensions - - >: greater than provided value - - >=: greater or equal to provided value - - <: smallerthan provided value - - <=: smaller or equal to provided value - - < vec: smallerthan provided values on each vertex - - forcing: 1 if check forcing consistency (size and time) - - values: cell of strings or vector of acceptable values - - numel: list of acceptable number of elements - - cell: 1 if check that is cell - - empty: 1 if check that non empty - - message: overloaded error message + Used to check model consistency. + Available options: + - NaN: 1 if check that there is no NaN + - size: [lines cols], NaN for non checked dimensions + - >: greater than provided value + - >=: greater or equal to provided value + - <: smallerthan provided value + - <=: smaller or equal to provided value + - < vec: smallerthan provided values on each vertex + - forcing: 1 if check forcing consistency (size and time) + - values: cell of strings or vector of acceptable values + - numel: list of acceptable number of elements + - cell: 1 if check that is cell + - empty: 1 if check that non empty + - message: overloaded error message - Usage: - md = checkfield(md,fieldname,options); + Usage: + md = checkfield(md,fieldname,options); - Example: - md = checkfield(md,'mesh.elementonbed','size',[md.mesh.numberofelements 1],'values',[0 1]); - md = checkfield(md,'diagnostic.icefront','size',[NaN 4],'NaN',1); - md = checkfield(md,'diagnostic.icefront(:,end)','values',[0 1 2]); + Example: + md = checkfield(md,'mesh.elementonbed','size',[md.mesh.numberofelements 1],'values',[0 1]); + md = checkfield(md,'diagnostic.icefront','size',[NaN 4],'NaN',1); + md = checkfield(md,'diagnostic.icefront(:,end)','values',[0 1 2]); """ @@ -37,8 +37,9 @@ #get field from model - field=getattr(md,fieldname) +# field=getattr(md,fieldname) + exec("field=md.%s" % fieldname) #check empty - if 'empty' in options: + if options.exist('empty'): if not field: md = md.checkmessage(options.getfieldvalue('message',\ @@ -46,21 +47,26 @@ #Check size - if 'size' in options: + if options.exist('size'): fieldsize=options.getfieldvalue('size') - if numpy.isnan(fieldsize[0]): - if not numpy.size(field,1)==fieldsize[1]: + if len(fieldsize) == 1: + if (not numpy.size(field,0)==fieldsize[0]): md = md.checkmessage(options.getfieldvalue('message',\ - "field '%s' should have %d columns" % (fieldname,fieldsize[1]))) - elif numpy.isnan(fieldsize[1]): - if not numpy.size(field,0)==fieldsize[0]: - md = md.checkmessage(options.getfieldvalue('message',\ - "field '%s' should have %d lines" % (fieldname,fieldsize[0]))) - else: - if (not numpy.size(field,0)==fieldsize[0]) or (not numpy.size(field,1)==fieldsize[1]): - md = md.checkmessage(options.getfieldvalue('message',\ - "field '%s' size should be %d x %d" % (fieldname,fieldsize[0],fieldsize[1]))) + "field '%s' size should be %d" % (fieldname,fieldsize[0]))) + elif len(fieldsize) == 2: + if numpy.isnan(fieldsize[0]): + if not numpy.size(field,1)==fieldsize[1]: + md = md.checkmessage(options.getfieldvalue('message',\ + "field '%s' should have %d columns" % (fieldname,fieldsize[1]))) + elif numpy.isnan(fieldsize[1]): + if not numpy.size(field,0)==fieldsize[0]: + md = md.checkmessage(options.getfieldvalue('message',\ + "field '%s' should have %d lines" % (fieldname,fieldsize[0]))) + else: + if (not numpy.size(field,0)==fieldsize[0]) or (not numpy.size(field,1)==fieldsize[1]): + md = md.checkmessage(options.getfieldvalue('message',\ + "field '%s' size should be %d x %d" % (fieldname,fieldsize[0],fieldsize[1]))) #Check numel - if 'numel' in options: + if options.exist('numel'): fieldnumel=options.getfieldvalue('numel') if not numpy.size(field) in fieldnumel: @@ -88,5 +94,5 @@ #check values - if 'values' in options: + if options.exist('values'): fieldvalues=options.getfieldvalue('values') if False in ismember(field,fieldvalues): @@ -102,10 +108,10 @@ #check greater - if '>=' in options: + if options.exist('>='): lowerbound=options.getfieldvalue('>=') if numpy.any(field' in options: + if options.exist('>'): lowerbound=options.getfieldvalue('>') if numpy.any(field<=lowerbound): @@ -114,10 +120,10 @@ #check smaller - if '<=' in options: + if options.exist('<='): upperbound=options.getfieldvalue('<=') if numpy.any(field>upperbound): md = md.checkmessage(options.getfieldvalue('message',\ "field '%s' should have values below %d" % (fieldname,upperbound))) - if '<' in options: + if options.exist('<'): upperbound=options.getfieldvalue('<') if numpy.any(field>=upperbound): Index: sm/branches/trunk-jpl-damage/src/py/model/solvers/asmoptions.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/py/model/solvers/asmoptions.py (revision 12945) +++ (revision ) @@ -1,31 +1,0 @@ -#module imports {{{ -from pairoptions import * -from petscversion import * -#}}} -def asmoptions(*args): - #ASMOPTIONS - return ASM petsc options - # - # Usage: - # options=asmoptions; - - #retrieve options provided in varargin - arguments=pairoptions(*args) - - options=[['mat_type','aij'],['ksp_type','gmres'],['pc_type','asm'],['sub_pc_type','lu'],['pc_asm_overlap',3],['ksp_max_it',100],['ksp_rtol',1e-30]]; - - #now, go through our arguments, and write over default options. - for i in range(len(arguments.list)): - arg1=arguments.list[i][0] - arg2=arguments.list[i][1] - found=0; - for j in range(len(options)): - joption=options[j][0] - if joption==arg1: - joption[1]=arg2; - options[j]=joption; - found=1; - break - if not found: - #this option did not exist, add it: - options.append([arg1,arg2]) - return options Index: sm/branches/trunk-jpl-damage/src/py/model/solvers/iluasmoptions.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/py/model/solvers/iluasmoptions.py (revision 12945) +++ (revision ) @@ -1,31 +1,0 @@ -#module imports {{{ -from pairoptions import * -from petscversion import * -#}}} -def iluasmoptions(*args): - #ILUASMOPTIONS - return MUMPS direct solver petsc options - # - # Usage: - # options=iluasmoptions; - - #retrieve options provided in varargin - arguments=pairoptions(*args) - - options=[['mat_type','aij'],['ksp_type','gmres'],['pc_type','asm'],['sub_pc_type','ilu'],['pc_asm_overlap',5],['ksp_max_it',100],['ksp_rtol',1e-15]]; - - #now, go through our arguments, and write over default options. - for i in range(len(arguments.list)): - arg1=arguments.list[i][0] - arg2=arguments.list[i][1] - found=0; - for j in range(len(options)): - joption=options[j][0] - if joption==arg1: - joption[1]=arg2; - options[j]=joption; - found=1; - break - if not found: - #this option did not exist, add it: - options.append([arg1,arg2]) - return options Index: sm/branches/trunk-jpl-damage/src/py/model/solvers/jacobiasmoptions.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/py/model/solvers/jacobiasmoptions.py (revision 12945) +++ (revision ) @@ -1,31 +1,0 @@ -#module imports {{{ -from pairoptions import * -from petscversion import * -#}}} -def jacobiasmoptions(*args): - #ASMOPTIONS - return Additive Shwartz Method with Jacobi preconditioner petsc options - # - # Usage: - # options=jacobiasmoptions; - - #retrieve options provided in varargin - arguments=pairoptions(*args) - - options=[['mat_type','aij'],['ksp_type','gmres'],['pc_type','asm'],['sub_pc_type','jacobi'],['pc_asm_overlap',3],['ksp_max_it',100],['ksp_rtol',1e-15]]; - - #now, go through our arguments, and write over default options. - for i in range(len(arguments.list)): - arg1=arguments.list[i][0] - arg2=arguments.list[i][1] - found=0; - for j in range(len(options)): - joption=options[j][0] - if joption==arg1: - joption[1]=arg2; - options[j]=joption; - found=1; - break - if not found: - #this option did not exist, add it: - options.append([arg1,arg2]) - return options Index: sm/branches/trunk-jpl-damage/src/py/model/solvers/jacobicgoptions.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/py/model/solvers/jacobicgoptions.py (revision 12945) +++ (revision ) @@ -1,31 +1,0 @@ -#module imports {{{ -from pairoptions import * -from petscversion import * -#}}} -def jacobicgoptions(*args): - #ASMOPTIONS - return Additive Shwartz Method with Jacobi preconditioner petsc options - # - # Usage: - # options=jacobicgoptions; - - #retrieve options provided in varargin - arguments=pairoptions(*args) - - options=[['mat_type','aij'],['ksp_type','cg'],['ksp_max_it',100],['ksp_rtol',1e-15]]; - - #now, go through our arguments, and write over default options. - for i in range(len(arguments.list)): - arg1=arguments.list[i][0] - arg2=arguments.list[i][1] - found=0; - for j in range(len(options)): - joption=options[j][0] - if joption==arg1: - joption[1]=arg2; - options[j]=joption; - found=1; - break - if not found: - #this option did not exist, add it: - options.append([arg1,arg2]) - return options Index: sm/branches/trunk-jpl-damage/src/py/model/solvers/matlaboptions.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/py/model/solvers/matlaboptions.py (revision 12945) +++ (revision ) @@ -1,31 +1,0 @@ -#module imports {{{ -from pairoptions import * -from petscversion import * -#}}} -def matlaboptions(*args): - #MATLABOPTIONS - return Matlab petsc options - # - # Usage: - # options=matlaboptions; - - #retrieve options provided in varargin - arguments=pairoptions(*args) - - options=[['ksp_type','matlab']]; - - #now, go through our arguments, and write over default options. - for i in range(len(arguments.list)): - arg1=arguments.list[i][0] - arg2=arguments.list[i][1] - found=0; - for j in range(len(options)): - joption=options[j][0] - if joption==arg1: - joption[1]=arg2; - options[j]=joption; - found=1; - break - if not found: - #this option did not exist, add it: - options.append([arg1,arg2]) - return options Index: sm/branches/trunk-jpl-damage/src/py/model/solvers/mumpsoptions.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/py/model/solvers/mumpsoptions.py (revision 12945) +++ (revision ) @@ -1,38 +1,0 @@ -#module imports {{{ -from pairoptions import * -from petscversion import * -#}}} -def mumpsoptions(*args): - #MUMPSOPTIONS - return MUMPS direct solver petsc options - # - # Usage: - # options=mumpsoptions; - - #retrieve options provided in varargin - arguments=pairoptions(*args) - - - #default mumps options - PETSC_VERSION=petscversion() - - if PETSC_VERSION==2: - options=[['mat_type','aijmumps'],['ksp_type','preonly'],['pc_type','lu'],['mat_mumps_icntl_14',120],['pc_factor_shift_positive_definite','true']] - if PETSC_VERSION==3: - options=[['mat_type','mpiaij'],['ksp_type','preonly'],['pc_type','lu'],['pc_factor_mat_solver_package','mumps'],['mat_mumps_icntl_14',120],['pc_factor_shift_positive_definite','true']] - - #now, go through our arguments, and write over default options. - for i in range(len(arguments.list)): - arg1=arguments.list[i][0] - arg2=arguments.list[i][1] - found=0; - for j in range(len(options)): - joption=options[j][0] - if joption==arg1: - joption[1]=arg2; - options[j]=joption; - found=1; - break - if not found: - #this option did not exist, add it: - options.append([arg1,arg2]) - return options Index: sm/branches/trunk-jpl-damage/src/py/model/solvers/soroptions.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/py/model/solvers/soroptions.py (revision 12945) +++ (revision ) @@ -1,31 +1,0 @@ -#module imports {{{ -from pairoptions import * -from petscversion import * -#}}} -def soroptions(*args): - #SOROPTIONS - return Relaxation Solver petsc options - # - # Usage: - # options=soroptions; - - #retrieve options provided in varargin - arguments=pairoptions(*args) - - options=[['mat_type','aij'],['ksp_type','cg'],['pc_type','sor'],['pc_sor_omega',1.1],['pc_sor_its',2]]; - - #now, go through our arguments, and write over default options. - for i in range(len(arguments.list)): - arg1=arguments.list[i][0] - arg2=arguments.list[i][1] - found=0; - for j in range(len(options)): - joption=options[j][0] - if joption==arg1: - joption[1]=arg2; - options[j]=joption; - found=1; - break - if not found: - #this option did not exist, add it: - options.append([arg1,arg2]) - return options Index: sm/branches/trunk-jpl-damage/src/py/model/solvers/stokesoptions.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/py/model/solvers/stokesoptions.py (revision 12945) +++ (revision ) @@ -1,40 +1,0 @@ -#module imports {{{ -from pairoptions import * -from petscversion import * -#}}} -def stokesoptions(*args): - #STOKESOPTIONS - return STOKES multi-physics solver petsc options - # - # Usage: - # options=stokesoptions; - - #retrieve options provided in varargin - arguments=pairoptions(*args) - - - #default stokes options - PETSC_VERSION=petscversion() - - if PETSC_VERSION==2: - raise RuntimeError('stokesoptions error message: multi-physics options not supported in Petsc 2') - if PETSC_VERSION==3: - options=[['mat_type','mpiaij'],['ksp_max_it',1000],['ksp_type','gmres'],['pc_type','fieldsplit'],['pc_field_split_type','schur'],\ - ['fieldsplit_0_pc_type','hypre'],['fieldsplit_0_ksp_type','gmres'],['fieldsplit_0_pc_hypre_type','boomerang'],\ - ['fieldsplit_1_pc_type','jacobi'],['fieldsplit_1_ksp_type','preonly'],['issm_option_solver','stokes']] - - #now, go through our arguments, and write over default options. - for i in range(len(arguments.list)): - arg1=arguments.list[i][0] - arg2=arguments.list[i][1] - found=0; - for j in range(len(options)): - joption=options[j][0] - if joption==arg1: - joption[1]=arg2; - options[j]=joption; - found=1; - break - if not found: - #this option did not exist, add it: - options.append([arg1,arg2]) - return options Index: /issm/branches/trunk-jpl-damage/test/Data/delta18o.data =================================================================== --- /issm/branches/trunk-jpl-damage/test/Data/delta18o.data (revision 12946) +++ /issm/branches/trunk-jpl-damage/test/Data/delta18o.data (revision 12946) @@ -0,0 +1,2 @@ + -3.9010000e-02 -3.9180000e-02 -3.9280000e-02 -3.9320000e-02 -3.9050000e-02 -3.8880000e-02 -3.8680000e-02 -3.8560000e-02 -3.8500000e-02 -3.8520000e-02 -3.8560000e-02 + 0.0000000e+00 2.0000000e+01 4.0000000e+01 6.0000000e+01 8.0000000e+01 1.0000000e+02 1.2000000e+02 1.4000000e+02 1.6000000e+02 1.8000000e+02 2.0000000e+02 Index: /issm/branches/trunk-jpl-damage/test/NightlyRun/IdToName.py =================================================================== --- /issm/branches/trunk-jpl-damage/test/NightlyRun/IdToName.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/test/NightlyRun/IdToName.py (revision 12946) @@ -66,5 +66,7 @@ 233 : 'SquareShelfTranP3dForcTemp', 234 : 'SquareShelfTranForceNeg2dDakotaSamp', - 235 : 'SquareShelfTranForceNeg2dDakotaLocal', + 235 : 'SquareShelfTranForceNeg2dDakotaLocal', + 236 : 'SquareShelfTranIspddIsdeltaM2d'; + 237 : 'SquareShelfTranIspddIsdeltaM3d'; 301 : 'SquareSheetConstrainedDiagM2d', 302 : 'SquareSheetConstrainedDiagH2d', Index: /issm/branches/trunk-jpl-damage/test/NightlyRun/test101.py =================================================================== --- /issm/branches/trunk-jpl-damage/test/NightlyRun/test101.py (revision 12945) +++ /issm/branches/trunk-jpl-damage/test/NightlyRun/test101.py (revision 12946) @@ -3,4 +3,8 @@ from setmask import * from parameterize import * +from setflowequation import * +from EnumDefinitions import * +from solve import * +from MatlabFuncs import * md=model() @@ -8,4 +12,8 @@ md=setmask(md,'all','') md=parameterize(md,'../Par/SquareShelfConstrained.py') +md=setflowequation(md,'macayeal','all') +md.cluster=generic('name',oshostname(),'np',3) +md.diagnostic.requested_outputs=StressTensorEnum +md=solve(md,DiagnosticSolutionEnum) #Fields and tolerances to track changes Index: /issm/branches/trunk-jpl-damage/test/NightlyRun/test236.m =================================================================== --- /issm/branches/trunk-jpl-damage/test/NightlyRun/test236.m (revision 12946) +++ /issm/branches/trunk-jpl-damage/test/NightlyRun/test236.m (revision 12946) @@ -0,0 +1,87 @@ +md=triangle(model,'../Exp/Square.exp',150000); +md=setmask(md,'all',''); + +% Use of ispdd and isdelta18o methods +md.surfaceforcings.ispdd=1; +md.surfaceforcings.isdelta18o=1; + +md=parameterize(md,'../Par/SquareShelf.par'); + +% Add temperature, precipitation and delta18o needed to measure the surface mass balance +% creating delta18o +load '../Data/delta18o.data' +md.surfaceforcings.delta18o=delta18o; +% creating delta18oSurface +md.surfaceforcings.delta18o_surface(1,1:(length(delta18o))) = 0; +md.surfaceforcings.delta18o_surface(2,:) = delta18o(2,:); + +% creating Present day and lgm temperatures +% Same temperature over the all region : +imonth=0:11; +tmonth(1:12)=238.15+20; +for imonth=0:11 + md.surfaceforcings.temperatures_presentday(1:md.mesh.numberofvertices,imonth+1)=tmonth(imonth+1); + md.surfaceforcings.temperatures_lgm(1:md.mesh.numberofvertices,imonth+1)=tmonth(imonth+1)-20; + % Time for the last line: + md.surfaceforcings.temperatures_presentday(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12); + md.surfaceforcings.temperatures_lgm(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12); +end + +% creating initialization and spc temperatures initialization and spc +md.thermal.spctemperature=mean(md.surfaceforcings.temperatures_lgm(1:md.mesh.numberofvertices,1:12),2); %-10*ones(md.mesh.numberofvertices,1); +md.thermal.spctemperature=repmat(md.thermal.spctemperature,1,md.timestepping.final_time/md.timestepping.time_step); +itemp = 0:md.timestepping.time_step:md.timestepping.final_time-md.timestepping.time_step; +md.thermal.spctemperature(md.mesh.numberofvertices+1,:)=itemp; + +md.initialization.temperature=md.surfaceforcings.temperatures_lgm(1:md.mesh.numberofvertices,1); %*ones(md.mesh.numberofvertices,1); + +% creating precipitation +for imonth=0:11 + md.surfaceforcings.precipitations_presentday(1:md.mesh.numberofvertices,imonth+1)=-0.4*10^(-6)*md.mesh.y+0.5; + md.surfaceforcings.precipitations_presentday(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12); +end + +% time steps and resolution +md.timestepping.time_step=20; +md.timestepping.final_time= 60; + +% +md=setflowequation(md,'macayeal','all'); +md.cluster=generic('name',oshostname(),'np',3); +md=solve(md,TransientSolutionEnum); + + +%Fields and tolerances to track changes +field_names ={'Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1','Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2','Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3'}; +field_tolerances={1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,... + 1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,... + 1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13}; +field_values={... + (md.results.TransientSolution(1).Vx),... + (md.results.TransientSolution(1).Vy),... + (md.results.TransientSolution(1).Vel),... + (md.results.TransientSolution(1).Pressure),... + (md.results.TransientSolution(1).Bed),... + (md.results.TransientSolution(1).Surface),... + (md.results.TransientSolution(1).Thickness),... + (md.results.TransientSolution(1).SurfaceforcingsMonthlytemperatures),... + (md.results.TransientSolution(1).SurfaceforcingsMassBalance),... + (md.results.TransientSolution(2).Vx),... + (md.results.TransientSolution(2).Vy),... + (md.results.TransientSolution(2).Vel),... + (md.results.TransientSolution(2).Pressure),... + (md.results.TransientSolution(2).Bed),... + (md.results.TransientSolution(2).Surface),... + (md.results.TransientSolution(2).Thickness),... + (md.results.TransientSolution(2).SurfaceforcingsMonthlytemperatures),... + (md.results.TransientSolution(2).SurfaceforcingsMassBalance),... + (md.results.TransientSolution(3).Vx),... + (md.results.TransientSolution(3).Vy),... + (md.results.TransientSolution(3).Vel),... + (md.results.TransientSolution(3).Pressure),... + (md.results.TransientSolution(3).Bed),... + (md.results.TransientSolution(3).Surface),... + (md.results.TransientSolution(3).Thickness),... + (md.results.TransientSolution(3).SurfaceforcingsMonthlytemperatures),... + (md.results.TransientSolution(3).SurfaceforcingsMassBalance),... + }; Index: /issm/branches/trunk-jpl-damage/test/NightlyRun/test237.m =================================================================== --- /issm/branches/trunk-jpl-damage/test/NightlyRun/test237.m (revision 12946) +++ /issm/branches/trunk-jpl-damage/test/NightlyRun/test237.m (revision 12946) @@ -0,0 +1,101 @@ +md=triangle(model,'../Exp/Square.exp',600000);%180000 +md=setmask(md,'all',''); + +% Use of ispdd and isdelta18o methods +md.surfaceforcings.ispdd=1; +md.surfaceforcings.isdelta18o=1; + +md=parameterize(md,'../Par/SquareShelf.par'); + +% Add temperature, precipitation and delta18o needed to measure the surface mass balance +% creating delta18o +load '../Data/delta18o.data' +md.surfaceforcings.delta18o=delta18o; +% creating delta18oSurface +md.surfaceforcings.delta18o_surface(1,1:(length(delta18o))) = 0; +md.surfaceforcings.delta18o_surface(2,:) = delta18o(2,:); + +% creating Present day and lgm temperatures +% Same temperature over the all region : +imonth=0:11; +tmonth(1:12)=238.15+20; +for imonth=0:11 + md.surfaceforcings.temperatures_presentday(1:md.mesh.numberofvertices,imonth+1)=tmonth(imonth+1); + md.surfaceforcings.temperatures_lgm(1:md.mesh.numberofvertices,imonth+1)=tmonth(imonth+1)-20; + % Time for the last line: + md.surfaceforcings.temperatures_presentday(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12); + md.surfaceforcings.temperatures_lgm(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12); +end + +% creating initialization and spc temperatures initialization and spc +md.thermal.spctemperature=mean(md.surfaceforcings.temperatures_lgm(1:md.mesh.numberofvertices,1:12),2); %-10*ones(md.mesh.numberofvertices,1); +md.thermal.spctemperature=repmat(md.thermal.spctemperature,1,md.timestepping.final_time/md.timestepping.time_step); +itemp = 0:md.timestepping.time_step:md.timestepping.final_time-md.timestepping.time_step; +md.thermal.spctemperature(md.mesh.numberofvertices+1,:)=itemp; + +md.initialization.temperature=md.surfaceforcings.temperatures_lgm(1:md.mesh.numberofvertices,1); %*ones(md.mesh.numberofvertices,1); + +% creating precipitation +for imonth=0:11 + md.surfaceforcings.precipitations_presentday(1:md.mesh.numberofvertices,imonth+1)=-0.4*10^(-6)*md.mesh.y+0.5; + md.surfaceforcings.precipitations_presentday(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12); +end + +% time steps and resolution +md.timestepping.time_step=20; +md.settings.output_frequency=1; +md.timestepping.final_time= 60; + +% +md=extrude(md,3,1); +md=setflowequation(md,'macayeal','all'); +md.cluster=generic('name',oshostname(),'np',3); +md=solve(md,TransientSolutionEnum); + +%Fields and tolerances to track changes +field_names ={'Vx1','Vy1','Vz1','Vel1','Pressure1','Bed1','Surface1','Thickness1','Temperature1','BasalforcingsMeltingRate1', ... + 'Vx2','Vy2','Vz2','Vel2','Pressure2','Bed2','Surface2','Thickness2','Temperature2','BasalforcingsMeltingRate2', ... + 'Vx3','Vy3','Vz3','Vel3','Pressure3','Bed3','Surface3','Thickness3','Temperature3','BasalforcingsMeltingRate3'}; +field_tolerances={... + 1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-8,... + 1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-8,... + 1e-13,1e-13,1e-08,1e-13,1e-13,1e-10,1e-10,1e-10,1e-13,1e-8}; + +field_values={... + (md.results.TransientSolution(1).Vx),... + (md.results.TransientSolution(1).Vy),... + (md.results.TransientSolution(1).Vz),... + (md.results.TransientSolution(1).Vel),... + (md.results.TransientSolution(1).Pressure),... + (md.results.TransientSolution(1).Bed),... + (md.results.TransientSolution(1).Surface),... + (md.results.TransientSolution(1).Thickness),... + (md.results.TransientSolution(1).Temperature),... + (md.results.TransientSolution(1).BasalforcingsMeltingRate),... + (md.results.TransientSolution(1).SurfaceforcingsMonthlytemperatures),... + (md.results.TransientSolution(1).SurfaceforcingsMassBalance),... + (md.results.TransientSolution(2).Vx),... + (md.results.TransientSolution(2).Vy),... + (md.results.TransientSolution(2).Vz),... + (md.results.TransientSolution(2).Vel),... + (md.results.TransientSolution(2).Pressure),... + (md.results.TransientSolution(2).Bed),... + (md.results.TransientSolution(2).Surface),... + (md.results.TransientSolution(2).Thickness),... + (md.results.TransientSolution(2).Temperature),... + (md.results.TransientSolution(2).BasalforcingsMeltingRate),... + (md.results.TransientSolution(2).SurfaceforcingsMonthlytemperatures),... + (md.results.TransientSolution(2).SurfaceforcingsMassBalance),... + (md.results.TransientSolution(3).Vx),... + (md.results.TransientSolution(3).Vy),... + (md.results.TransientSolution(3).Vz),... + (md.results.TransientSolution(3).Vel),... + (md.results.TransientSolution(3).Pressure),... + (md.results.TransientSolution(3).Bed),... + (md.results.TransientSolution(3).Surface),... + (md.results.TransientSolution(3).Thickness),... + (md.results.TransientSolution(3).Temperature),... + (md.results.TransientSolution(3).BasalforcingsMeltingRate),... + (md.results.TransientSolution(3).SurfaceforcingsMonthlytemperatures),... + (md.results.TransientSolution(3).SurfaceforcingsMassBalance),... + };