[issm-svn] r14067 - issm/trunk
morlighe at issm.ess.uci.edu
morlighe at issm.ess.uci.edu
Fri Nov 30 10:09:49 PST 2012
Author: morlighe
Date: 2012-11-30 10:09:49 -0800 (Fri, 30 Nov 2012)
New Revision: 14067
Added:
issm/trunk/externalpackages/gsl/Makefile.am.patch
issm/trunk/src/android/ISSM/jni/Application.mk
issm/trunk/src/android/ISSM/res/layout/activity_mapselection.xml
issm/trunk/src/android/ISSM/src/com/example/issm/MapSelection.java
issm/trunk/test/NightlyRun/test1501.py
issm/trunk/test/NightlyRun/test1502.py
Removed:
issm/trunk/src/android/helloworld/
issm/trunk/src/android/two-libs/
issm/trunk/test/NightlyRun/InNeedOfDebugging/test1501.py
issm/trunk/test/NightlyRun/InNeedOfDebugging/test1502.py
issm/trunk/test/NightlyRun/InNeedOfDebugging/test1601.py
issm/trunk/test/NightlyRun/InNeedOfDebugging/test1602.py
issm/trunk/test/NightlyRun/InNeedOfDebugging/test3001.py
issm/trunk/test/NightlyRun/InNeedOfDebugging/test3002.py
issm/trunk/test/NightlyRun/InNeedOfDebugging/test3003.py
issm/trunk/test/NightlyRun/InNeedOfDebugging/test3004.py
issm/trunk/test/NightlyRun/InNeedOfDebugging/test3005.py
issm/trunk/test/NightlyRun/InNeedOfDebugging/test3006.py
issm/trunk/test/NightlyRun/InNeedOfDebugging/test3007.py
issm/trunk/test/NightlyRun/InNeedOfDebugging/test3008.py
issm/trunk/test/NightlyRun/InNeedOfDebugging/test3009.py
issm/trunk/test/NightlyRun/InNeedOfDebugging/test3010.py
Modified:
issm/trunk/
issm/trunk/aux-config/
issm/trunk/configs/config-greenplanet.sh
issm/trunk/configure.ac
issm/trunk/etc/environment.sh
issm/trunk/externalpackages/chaco/
issm/trunk/externalpackages/doxygen/install.sh
issm/trunk/externalpackages/gsl/install-android.sh
issm/trunk/externalpackages/triangle/configs/android/configure.make
issm/trunk/externalpackages/triangle/install-android.sh
issm/trunk/src/
issm/trunk/src/android/ISSM/
issm/trunk/src/android/ISSM/.classpath
issm/trunk/src/android/ISSM/AndroidManifest.xml
issm/trunk/src/android/ISSM/bin/
issm/trunk/src/android/ISSM/bin/AndroidManifest.xml
issm/trunk/src/android/ISSM/gen/com/example/issm/R.java
issm/trunk/src/android/ISSM/jni/Android.mk
issm/trunk/src/android/ISSM/jni/Main.cpp
issm/trunk/src/android/ISSM/jni/issmlib/Android.mk
issm/trunk/src/android/ISSM/src/com/example/issm/ISSM.java
issm/trunk/src/android/ISSM/src/com/example/issm/IssmJni.java
issm/trunk/src/android/ISSM_Visual/bin/
issm/trunk/src/android/ISSM_Visual/bin/resources.ap_
issm/trunk/src/android/ISSM_Visual/src/com/example/issm_visual/
issm/trunk/src/c/Container/Observations.cpp
issm/trunk/src/c/Container/Observations.h
issm/trunk/src/c/Container/Options.h
issm/trunk/src/c/android/fac.cpp
issm/trunk/src/c/classes/FemModel.cpp
issm/trunk/src/c/classes/FemModel.h
issm/trunk/src/c/classes/objects/Elements/Penta.cpp
issm/trunk/src/c/modules/ElementConnectivityx/ElementConnectivityx.cpp
issm/trunk/src/c/modules/ElementConnectivityx/ElementConnectivityx.h
issm/trunk/src/c/modules/Krigingx/Krigingx.cpp
issm/trunk/src/c/modules/Krigingx/Krigingx.h
issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.cpp
issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.h
issm/trunk/src/c/shared/Threads/LaunchThread.cpp
issm/trunk/src/dox/issm.dox
issm/trunk/src/m/boundaryconditions/SetIceShelfBC.py
issm/trunk/src/m/boundaryconditions/SetMarineIceSheetBC.py
issm/trunk/src/m/classes/clusters/greenplanet.m
issm/trunk/src/m/classes/flowequation.py
issm/trunk/src/m/classes/initialization.py
issm/trunk/src/m/classes/mask.py
issm/trunk/src/m/classes/mesh.py
issm/trunk/src/m/classes/model/model.py
issm/trunk/src/m/classes/transient.m
issm/trunk/src/m/classes/transient.py
issm/trunk/src/m/consistency/checkfield.py
issm/trunk/src/m/exp/expcoarsen.m
issm/trunk/src/m/extrusion/project3d.py
issm/trunk/src/m/geometry/FlagElements.py
issm/trunk/src/m/geometry/GetAreas.py
issm/trunk/src/m/mesh/ComputeHessian.py
issm/trunk/src/m/mesh/GetNodalFunctionsCoeff.py
issm/trunk/src/m/mesh/bamg.py
issm/trunk/src/m/mesh/meshconvert.py
issm/trunk/src/m/mesh/rifts/meshprocessoutsiderifts.py
issm/trunk/src/m/mesh/rifts/meshprocessrifts.py
issm/trunk/src/m/mesh/triangle.py
issm/trunk/src/m/parameterization/contourenvelope.py
issm/trunk/src/m/parameterization/setflowequation.py
issm/trunk/src/m/parameterization/setmask.py
issm/trunk/src/m/plot/applyoptions.m
issm/trunk/src/m/plot/colormaps/getcolormap.m
issm/trunk/src/m/plot/radarpower.m
issm/trunk/src/m/solve/WriteData.py
issm/trunk/src/wrappers/ElementConnectivity/ElementConnectivity.cpp
issm/trunk/src/wrappers/Kriging/Kriging.cpp
issm/trunk/src/wrappers/NodeConnectivity/NodeConnectivity.cpp
issm/trunk/src/wrappers/matlab/io/FetchMatlabData.cpp
issm/trunk/src/wrappers/python/io/FetchPythonData.cpp
issm/trunk/src/wrappers/python/io/WritePythonData.cpp
issm/trunk/src/wrappers/python/io/pythonio.h
issm/trunk/test/
issm/trunk/test/Archives/Archive1205.nc
issm/trunk/test/Archives/Archive1206.nc
issm/trunk/test/Archives/Archive1207.nc
issm/trunk/test/NightlyRun/runme.m
issm/trunk/test/NightlyRun/runme.py
issm/trunk/test/NightlyRun/test109.py
issm/trunk/test/NightlyRun/test121.py
issm/trunk/test/NightlyRun/test1501.m
issm/trunk/test/NightlyRun/test1502.m
issm/trunk/test/NightlyRun/test207.py
issm/trunk/test/NightlyRun/test228.py
issm/trunk/test/NightlyRun/test229.py
issm/trunk/test/NightlyRun/test230.py
issm/trunk/test/NightlyRun/test231.py
issm/trunk/test/NightlyRun/test232.py
issm/trunk/test/NightlyRun/test3009.py
issm/trunk/test/NightlyRun/test313.py
issm/trunk/test/NightlyRun/test326.py
issm/trunk/test/NightlyRun/test407.py
issm/trunk/test/NightlyRun/test423.py
issm/trunk/test/NightlyRun/test424.py
issm/trunk/test/NightlyRun/test425.py
issm/trunk/test/NightlyRun/test426.py
issm/trunk/test/NightlyRun/test427.py
issm/trunk/test/NightlyRun/test515.py
Log:
merged trunk-jpl and trunk for revision 14066
Property changes on: issm/trunk
___________________________________________________________________
Modified: svn:mergeinfo
- /issm/trunk-jpl:10936-11994,11996-12326,12333-12336,12338-12703,12705-12707,12710-13393,13396-13974
+ /issm/trunk-jpl:10936-11994,11996-12326,12333-12336,12338-12703,12705-12707,12710-13393,13396-13974,13976-14066
Property changes on: issm/trunk/aux-config
___________________________________________________________________
Modified: svn:ignore
-
+ ar-lib
depcomp
compile
missing
config.guess
config.sub
ltmain.sh
install-sh
Modified: issm/trunk/configs/config-greenplanet.sh
===================================================================
--- issm/trunk/configs/config-greenplanet.sh 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/configs/config-greenplanet.sh 2012-11-30 18:09:49 UTC (rev 14067)
@@ -21,5 +21,5 @@
--with-scalapack-dir=$ISSM_DIR/externalpackages/petsc/install/ \
--with-blacs-dir=$ISSM_DIR/externalpackages/petsc/install/ \
--with-graphics-lib=/usr/lib64/libX11.so \
- --with-cxxoptflags="-O3 -xS" \
+ --with-cxxoptflags="-O3" \
--with-vendor=intel-linux
Modified: issm/trunk/configure.ac
===================================================================
--- issm/trunk/configure.ac 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/configure.ac 2012-11-30 18:09:49 UTC (rev 14067)
@@ -1,7 +1,7 @@
# Process this file with autoconf to produce a configure script.
#AUTOCONF
-AC_INIT([ISSM],[4.2.3],[issm at jpl.nasa.gov],[issm],[http://issm.jpl.nasa.gov]) #Initializing configure
+AC_INIT([ISSM],[4.2.4],[issm at jpl.nasa.gov],[issm],[http://issm.jpl.nasa.gov]) #Initializing configure
AC_CONFIG_AUX_DIR([./aux-config]) #Put config files in aux-config
AC_CONFIG_MACRO_DIR([m4]) #m4 macros are located in m4
m4_include([m4/issm_options.m4])
Modified: issm/trunk/etc/environment.sh
===================================================================
--- issm/trunk/etc/environment.sh 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/etc/environment.sh 2012-11-30 18:09:49 UTC (rev 14067)
@@ -85,7 +85,7 @@
pathappend "$DAKOTA_DIR/bin"
DOXYGEN_DIR="$ISSM_DIR/externalpackages/doxygen/install"
-pathappend "$DOXYGEN_DIR/bin"
+pathprepend "$DOXYGEN_DIR/bin"
AUTOTOOLS_DIR="$ISSM_DIR/externalpackages/autotools/install"
pathprepend "$AUTOTOOLS_DIR/bin"
Property changes on: issm/trunk/externalpackages/chaco
___________________________________________________________________
Modified: svn:ignore
- install
src
.ignore.txt
old
+ *.pdf
install
src
.ignore.txt
old
*.tar.gz
Modified: issm/trunk/externalpackages/doxygen/install.sh
===================================================================
--- issm/trunk/externalpackages/doxygen/install.sh 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/externalpackages/doxygen/install.sh 2012-11-30 18:09:49 UTC (rev 14067)
@@ -2,10 +2,10 @@
set -eu
#Some cleanup
-rm -rf install
+rm -rf install src
#Download latest version
-svn co https://doxygen.svn.sourceforge.net/svnroot/doxygen/trunk install
+svn co https://doxygen.svn.sourceforge.net/svnroot/doxygen/trunk src
#Configure doxygen
cd install && ./configure --prefix "$ISSM_DIR/externalpackages/doxygen/install"
Modified: issm/trunk/externalpackages/gsl/install-android.sh
===================================================================
--- issm/trunk/externalpackages/gsl/install-android.sh 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/externalpackages/gsl/install-android.sh 2012-11-30 18:09:49 UTC (rev 14067)
@@ -24,8 +24,9 @@
if [[ $step == "2" || $step == "0" ]]; then
cd src
+ mv ./../Makefile.am.patch ./
patch Makefile.am < Makefile.am.patch
-
+
autoreconf -iv --force -I $ISSM_DIR/externalpackages/autotools/install/share/aclocal
./configure \
@@ -37,9 +38,10 @@
#Compile gsl
if [[ $step == "3" || $step == "0" ]]; then
- cd src
+ cd $ISSM_DIR/externalpackages/gsl/src
+
if [ $# -eq 0 ]; then
- make $j
+ make
else
make -j $j
fi
Modified: issm/trunk/externalpackages/triangle/configs/android/configure.make
===================================================================
--- issm/trunk/externalpackages/triangle/configs/android/configure.make 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/externalpackages/triangle/configs/android/configure.make 2012-11-30 18:09:49 UTC (rev 14067)
@@ -9,11 +9,9 @@
#
# http://www.codesourcery.com/gnu_toolchains/arm/arm_gnu_linux_abi.pdf
-source $ANDROID_DIR/android_aux.sh
-
-CC=${toolchain_path}-gcc
-AR=${toolchain_path}-ar
-RANLIB=${toolchain_path}-ranlib
+CC=${host_triplet}-gcc
+AR=${host_triplet}-ar
+RANLIB=${host_triplet}-ranlib
CSWITCHES = $(CFLAGS)
TRILIBDEFS = -DTRILIBRARY
OBJ_EXT=o
Modified: issm/trunk/externalpackages/triangle/install-android.sh
===================================================================
--- issm/trunk/externalpackages/triangle/install-android.sh 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/externalpackages/triangle/install-android.sh 2012-11-30 18:09:49 UTC (rev 14067)
@@ -1,6 +1,8 @@
#!/bin/bash
set -eu
+source $ANDROID_DIR/android_aux.sh
+
#Some cleanup
rm -rf install triangle
mkdir install
Property changes on: issm/trunk/src
___________________________________________________________________
Modified: svn:mergeinfo
- /issm/branches/trunk-jpl-damage/src:11427-13118
/issm/trunk-jpl/src:10936-11994,11996-12326,12333-12336,12338-12703,12705-12707,12710-13393,13396-13974
+ /issm/branches/trunk-jpl-damage/src:11427-13118
/issm/trunk-jpl/src:10936-11994,11996-12326,12333-12336,12338-12703,12705-12707,12710-13393,13396-13974,13976-14066
Property changes on: issm/trunk/src/android/ISSM
___________________________________________________________________
Modified: svn:ignore
- Makefile.in
Makefile
libs
obj
+ gen
.settings
libs
obj
Makefile
Makefile.in
Modified: issm/trunk/src/android/ISSM/.classpath
===================================================================
--- issm/trunk/src/android/ISSM/.classpath 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/android/ISSM/.classpath 2012-11-30 18:09:49 UTC (rev 14067)
@@ -1,8 +1,8 @@
<?xml version="1.0" encoding="UTF-8"?>
<classpath>
+ <classpathentry kind="con" path="com.android.ide.eclipse.adt.ANDROID_FRAMEWORK"/>
+ <classpathentry kind="con" path="com.android.ide.eclipse.adt.LIBRARIES"/>
<classpathentry kind="src" path="src"/>
<classpathentry kind="src" path="gen"/>
- <classpathentry kind="con" path="com.android.ide.eclipse.adt.ANDROID_FRAMEWORK"/>
- <classpathentry kind="con" path="com.android.ide.eclipse.adt.LIBRARIES"/>
<classpathentry kind="output" path="bin/classes"/>
</classpath>
Modified: issm/trunk/src/android/ISSM/AndroidManifest.xml
===================================================================
--- issm/trunk/src/android/ISSM/AndroidManifest.xml 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/android/ISSM/AndroidManifest.xml 2012-11-30 18:09:49 UTC (rev 14067)
@@ -4,21 +4,23 @@
android:versionName="1.0" >
<uses-sdk
- android:minSdkVersion="15"
+ android:minSdkVersion="10"
android:targetSdkVersion="15" />
-
+ <uses-permission android:name="android.permission.WRITE_EXTERNAL_STORAGE" />
<application
android:icon="@drawable/ic_launcher"
android:label="@string/app_name"
android:theme="@style/AppTheme" >
- <activity
- android:name=".ISSM"
- android:label="@string/title_activity_issm" >
+ <activity android:name=".MapSelection"
+ android:label="@string/title_activity_issm" >
<intent-filter>
<action android:name="android.intent.action.MAIN" />
-
- <category android:name="android.intent.category.LAUNCHER" />
+ <category android:name="android.intent.category.LAUNCHER" />
</intent-filter>
+ </activity>
+ <activity
+ android:name=".ISSM"
+ android:label="@string/title_activity_issm" >
</activity>
</application>
Property changes on: issm/trunk/src/android/ISSM/bin
___________________________________________________________________
Modified: svn:ignore
- ISSM.apk
resources.ap_
res
jarlist.cache
classes.dex
classes
+ *ISSM.apk
resources.ap_
res
jarlist.cache
classes.dex
classes
Modified: issm/trunk/src/android/ISSM/bin/AndroidManifest.xml
===================================================================
--- issm/trunk/src/android/ISSM/bin/AndroidManifest.xml 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/android/ISSM/bin/AndroidManifest.xml 2012-11-30 18:09:49 UTC (rev 14067)
@@ -4,21 +4,23 @@
android:versionName="1.0" >
<uses-sdk
- android:minSdkVersion="15"
+ android:minSdkVersion="10"
android:targetSdkVersion="15" />
-
+ <uses-permission android:name="android.permission.WRITE_EXTERNAL_STORAGE" />
<application
android:icon="@drawable/ic_launcher"
android:label="@string/app_name"
android:theme="@style/AppTheme" >
- <activity
- android:name=".ISSM"
- android:label="@string/title_activity_issm" >
+ <activity android:name=".MapSelection"
+ android:label="@string/title_activity_issm" >
<intent-filter>
<action android:name="android.intent.action.MAIN" />
-
- <category android:name="android.intent.category.LAUNCHER" />
+ <category android:name="android.intent.category.LAUNCHER" />
</intent-filter>
+ </activity>
+ <activity
+ android:name=".ISSM"
+ android:label="@string/title_activity_issm" >
</activity>
</application>
Modified: issm/trunk/src/android/ISSM/gen/com/example/issm/R.java
===================================================================
--- issm/trunk/src/android/ISSM/gen/com/example/issm/R.java 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/android/ISSM/gen/com/example/issm/R.java 2012-11-30 18:09:49 UTC (rev 14067)
@@ -21,12 +21,14 @@
}
public static final class id {
public static final int button1=0x7f080001;
+ public static final int button2=0x7f080003;
public static final int input=0x7f080000;
- public static final int menu_settings=0x7f080003;
+ public static final int menu_settings=0x7f080004;
public static final int output=0x7f080002;
}
public static final class layout {
public static final int activity_issm=0x7f030000;
+ public static final int activity_mapselection=0x7f030001;
}
public static final class menu {
public static final int activity_issm=0x7f070000;
Modified: issm/trunk/src/android/ISSM/jni/Android.mk
===================================================================
--- issm/trunk/src/android/ISSM/jni/Android.mk 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/android/ISSM/jni/Android.mk 2012-11-30 18:09:49 UTC (rev 14067)
@@ -10,7 +10,7 @@
include $(CLEAR_VARS)
LOCAL_PATH := $(my_LOCAL_PATH)
LOCAL_ALLOW_UNDEFINED_SYMBOLS := true
-LOCAL_CFLAGS := -Wno-psabi
+LOCAL_CFLAGS := -Wno-psabi -DHAVE_CONFIG_H
LOCAL_LDLIBS := -llog -ldl
LOCAL_MODULE := IssmJni
LOCAL_SRC_FILES := Main.cpp
Modified: issm/trunk/src/android/ISSM/jni/Main.cpp
===================================================================
--- issm/trunk/src/android/ISSM/jni/Main.cpp 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/android/ISSM/jni/Main.cpp 2012-11-30 18:09:49 UTC (rev 14067)
@@ -1,20 +1,56 @@
#include <jni.h>
#include "../../../c/android/fac.h"
+#include "../../../c/issm.h"
+#include <cstddef>
+#include <stdio.h>
+///////////////////////////////////////////////////////////////////////////////////////////
namespace com_example_issm
{
- static jlong factorial(JNIEnv *env, jclass clazz, jlong n)
+ fac* f;
+ FemModel *fm;
+//------------------------------------------------------------------------------------
+ jint initilize(JNIEnv *env, jclass clazz, jstring file)
{
- fac *f = new fac();
- jlong result = (jlong) (f->factorial(n));
- delete(f);
- return (jlong) result;
+ char *nativefile = (char*)env->GetStringUTFChars(file,0);
+
+ f = new fac();
+
+ //call Model constructor passing in infile as File Descriptor parameter.
+ fm = new FemModel(nativefile);
+
+ env->ReleaseStringUTFChars(file, nativefile); //must realease the char*
+ //jint size = (jint) fm->bufferSize();
+ //return size;
+
+ return 0; //return 0 for now.
}
//------------------------------------------------------------------------------------
- static JNINativeMethod method_table[] = {
- {"fac" ,"(J)J" , (void *) factorial},
+ //fill out the first two slots, extract the same way in java using get(int position)
+ void solve(JNIEnv *env, jclass clazz , jint alpha, jobject buf)
+ {
+ jdouble *dBuf = (jdouble *)env->GetDirectBufferAddress(buf);
+
+ //pass bBuff to fem model to allocate data
+ // fm -> solve(dBuf, alpha);
+
+ }
+//------------------------------------------------------------------------------------
+ jlong factorial(JNIEnv *env, jclass clazz, jlong n)
+ {
+ if( f != NULL)
+ return (jlong) (f->factorial(n));
+ return 0;
+ }
+
+//------------------------------------------------------------------------------------
+ static JNINativeMethod method_table[] =
+ {
+ {"fac" , "(J)J" , (void *) factorial},
+ {"createISSMModel" ,"(Ljava/lang/String;)I" , (void *) initilize},
+ {"processBuffer", "(ILjava/nio/DoubleBuffer;)V", (void *) solve}
};
}
-
+//------------------------------------------------------------------------------------
using namespace com_example_issm;
extern "C" jint JNI_OnLoad(JavaVM* vm, void* reserved)
@@ -28,13 +64,11 @@
jclass clazz = env->FindClass("com/example/issm/IssmJni");
if(clazz)
{
- env->RegisterNatives(clazz, method_table, 1);
+ env->RegisterNatives(clazz, method_table, 3);
env->DeleteLocalRef(clazz);
return JNI_VERSION_1_6;
}
else return -1;
}
-
- // Get jclass with env->FindClass.
- // Register methods with env->RegisterNatives.
-}
\ No newline at end of file
+}
+///////////////////////////////////////////////////////////////////////////////////////////
Modified: issm/trunk/src/android/ISSM/jni/issmlib/Android.mk
===================================================================
--- issm/trunk/src/android/ISSM/jni/issmlib/Android.mk 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/android/ISSM/jni/issmlib/Android.mk 2012-11-30 18:09:49 UTC (rev 14067)
@@ -1,6 +1,6 @@
LOCAL_PATH:= $(call my-dir)
include $(CLEAR_VARS)
LOCAL_MODULE := libISSMCore
-LOCAL_SRC_FILES := ../../../../c/libISSMCore.a
-LOCAL_EXPORT_C_INCLUDES := ../../../../c/android/
+LOCAL_SRC_FILES := libISSMCore.a
+LOCAL_EXPORT_C_INCLUDES := $(ISSM_DIR)
include $(PREBUILT_STATIC_LIBRARY)
Modified: issm/trunk/src/android/ISSM/src/com/example/issm/ISSM.java
===================================================================
--- issm/trunk/src/android/ISSM/src/com/example/issm/ISSM.java 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/android/ISSM/src/com/example/issm/ISSM.java 2012-11-30 18:09:49 UTC (rev 14067)
@@ -1,7 +1,14 @@
package com.example.issm;
+import java.io.IOException;
+import java.io.InputStream;
+import java.nio.ByteBuffer;
+import java.nio.ByteOrder;
+import java.nio.DoubleBuffer;
+
import android.os.Bundle;
import android.app.Activity;
+import android.content.res.AssetManager;
import android.text.TextUtils;
import android.view.Menu;
import android.view.View;
@@ -11,29 +18,72 @@
import android.widget.TextView;
-public class ISSM extends Activity implements OnClickListener {
+public class ISSM extends Activity implements OnClickListener
+{
private EditText input;
private TextView output;
+ private DoubleBuffer buff;
+ private IssmJni issmNative;
+ private String mapName;
+ private String issmFolder;
+ private int size;
@Override
+ //------------------------------------------------------------------------------------------------
public void onCreate(Bundle savedInstanceState) {
super.onCreate(savedInstanceState);
+ Bundle map = getIntent().getExtras();
+ {
+ if(map!= null)
+ {
+ mapName = map.getString("map");
+ issmFolder = map.getString("pathToFile");
+ }
+ }
setContentView(R.layout.activity_issm);
this.input = (EditText) super.findViewById(R.id.input);
this.output = (TextView) super.findViewById(R.id.output);
Button button = (Button) super.findViewById(R.id.button1);
button.setOnClickListener(this);
+ //load up the ISSM library and create double buffer in java
+ //which later on will be pass for native allocation.
+ issmNative = new IssmJni();
+ buff = ByteBuffer.allocateDirect(15*8).order(ByteOrder.nativeOrder()).asDoubleBuffer();
+ this.createModel();
}
-
- public void onClick(View view)
+//------------------------------------------------------------------------------------------------
+ public void createModel()
+ {
+ String file = "";
+ if( mapName.equals("greenland"))
+ {
+ file = "test102.petsc";
+ }
+ else file = "test102.petsc";
+ size = issmNative.createISSMModel(issmFolder + file);
+ }
+//------------------------------------------------------------------------------------------------
+ public void fillBuffer()
+ {
+ issmNative.processBuffer(5,buff);
+ }
+ //------------------------------------------------------------------------------------------------
+ public void onClick(View view)
{
- // TODO Auto-generated method stub
+ //factorial method
String input = this.input.getText().toString();
if(TextUtils.isEmpty(input))
{
return;
- }
- long result = IssmJni.facIterative(Long.parseLong(input));
- this.output.setText("Result = " + result + "\n");
+ }
+ long resultfromFac = issmNative.fac(Long.parseLong(input));
+ //example of how to fill buffer Native
+ this.fillBuffer();
+
+ //print result from fac and the first two slot of filled buffer.
+ this.output.setText("Result = " + resultfromFac + "\n"
+ + "First slot from buffer:\n"
+ + buff.get(0) + " size " + size);
+
}
}
Modified: issm/trunk/src/android/ISSM/src/com/example/issm/IssmJni.java
===================================================================
--- issm/trunk/src/android/ISSM/src/com/example/issm/IssmJni.java 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/android/ISSM/src/com/example/issm/IssmJni.java 2012-11-30 18:09:49 UTC (rev 14067)
@@ -1,13 +1,13 @@
package com.example.issm;
+import java.nio.DoubleBuffer;
-public class IssmJni
+class IssmJni
{
- public static native long fac(long n);
- static {
+ public native long fac(long n);
+ public native void processBuffer(int alpha, DoubleBuffer buff);
+ public native int createISSMModel(String file);
+ static
+ {
System.loadLibrary("IssmJni");
}
- public static long facIterative(long n)
- {
- return fac(n);
- }
}
Property changes on: issm/trunk/src/android/ISSM_Visual/bin
___________________________________________________________________
Added: svn:ignore
+ *.apk
Modified: issm/trunk/src/android/ISSM_Visual/bin/resources.ap_
===================================================================
(Binary files differ)
Property changes on: issm/trunk/src/android/ISSM_Visual/src/com/example/issm_visual
___________________________________________________________________
Added: svn:ignore
+ *.java
*.class
Modified: issm/trunk/src/c/Container/Observations.cpp
===================================================================
--- issm/trunk/src/c/Container/Observations.cpp 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/c/Container/Observations.cpp 2012-11-30 18:09:49 UTC (rev 14067)
@@ -423,6 +423,90 @@
/*Assign output pointer*/
*pprediction = obs;
}/*}}}*/
+/*FUNCTION Observations::InterpolationV4{{{*/
+void Observations::InterpolationV4(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata){
+ /* Reference: David T. Sandwell, Biharmonic spline interpolation of GEOS-3
+ * and SEASAT altimeter data, Geophysical Research Letters, 2, 139-142,
+ * 1987. Describes interpolation using value or gradient of value in any
+ * dimension.*/
+
+ /*Intermediaries*/
+ int i,j,n_obs;
+ IssmPDouble prediction,h;
+ IssmPDouble *x = NULL;
+ IssmPDouble *y = NULL;
+ IssmPDouble *obs = NULL;
+ IssmPDouble *Green = NULL;
+ IssmPDouble *weights = NULL;
+ IssmPDouble *g = NULL;
+
+ /*Some checks*/
+ _assert_(maxdata>0);
+ _assert_(pprediction);
+
+ /*If radius is not provided or is 0, return all observations*/
+ if(radius==0) radius=this->quadtree->root->length;
+
+ /*Get list of observations for current point*/
+ this->ObservationList(&x,&y,&obs,&n_obs,x_interp,y_interp,radius,maxdata);
+
+ /*If we have less observations than mindata, return UNDEF*/
+ if(n_obs<mindata || n_obs<2){
+ prediction = UNDEF;
+ }
+ else{
+
+ /*Allocate intermediary matrix and vectors*/
+ Green = xNew<IssmPDouble>(n_obs*n_obs);
+ g = xNew<IssmPDouble>(n_obs);
+
+ /*First: distance vector*/
+ for(i=0;i<n_obs;i++){
+ h = sqrt( (x[i]-x_interp)*(x[i]-x_interp) + (y[i]-y_interp)*(y[i]-y_interp) );
+ if(h>0){
+ g[i] = h*h*(log(h)-1.);
+ }
+ else{
+ g[i] = 0.;
+ }
+ }
+
+ /*Build Green function matrix*/
+ for(i=0;i<n_obs;i++){
+ for(j=0;j<=i;j++){
+ h = sqrt( (x[i]-x[j])*(x[i]-x[j]) + (y[i]-y[j])*(y[i]-y[j]) );
+ if(h>0){
+ Green[j*n_obs+i] = h*h*(log(h)-1.);
+ }
+ else{
+ Green[j*n_obs+i] = 0.;
+ }
+ Green[i*n_obs+j] = Green[j*n_obs+i];
+ }
+ }
+
+ /*Compute weights*/
+#if _HAVE_GSL_
+ SolverxSeq(&weights,Green,obs,n_obs); // Green^-1 obs
+#else
+ _error_("GSL is required");
+#endif
+
+ /*Interpolate*/
+ prediction = 0;
+ for(i=0;i<n_obs;i++) prediction += weights[i]*g[i];
+
+ }
+
+ /*clean-up*/
+ *pprediction = prediction;
+ xDelete<IssmPDouble>(x);
+ xDelete<IssmPDouble>(y);
+ xDelete<IssmPDouble>(obs);
+ xDelete<IssmPDouble>(Green);
+ xDelete<IssmPDouble>(g);
+ xDelete<IssmPDouble>(weights);
+}/*}}}*/
/*FUNCTION Observations::QuadtreeColoring{{{*/
void Observations::QuadtreeColoring(IssmPDouble* A,IssmPDouble *x,IssmPDouble *y,int n){
Modified: issm/trunk/src/c/Container/Observations.h
===================================================================
--- issm/trunk/src/c/Container/Observations.h 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/c/Container/Observations.h 2012-11-30 18:09:49 UTC (rev 14067)
@@ -27,6 +27,7 @@
/*Methods*/
void ClosestObservation(IssmDouble *px,IssmDouble *py,IssmDouble *pobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);
void InterpolationIDW(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata,IssmDouble power);
+ void InterpolationV4(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata);
void InterpolationKriging(IssmDouble *pprediction,IssmDouble *perror,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata,Variogram* variogram);
void InterpolationNearestNeighbor(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);
void ObservationList(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs);
Modified: issm/trunk/src/c/Container/Options.h
===================================================================
--- issm/trunk/src/c/Container/Options.h 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/c/Container/Options.h 2012-11-30 18:09:49 UTC (rev 14067)
@@ -63,6 +63,7 @@
genericoption->Get(pvalue);
}
else{
+ if(GetOption(name)) _printLine_("WARNING: option "<<name<<" found but fetched format not consistent, defaulting...");
*pvalue=default_value;
}
}
Modified: issm/trunk/src/c/android/fac.cpp
===================================================================
--- issm/trunk/src/c/android/fac.cpp 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/c/android/fac.cpp 2012-11-30 18:09:49 UTC (rev 14067)
@@ -12,14 +12,16 @@
#include "./fac.h"
/*}}}*/
-
+#include "stdio.h"
long fac::factorial(long n) {
long f = 1;
long i;
+ printf("ok1\n");
for(i = 1; i <= n; i++){
f *= i;
}
+ printf("ok2\n");
return f;
}
Modified: issm/trunk/src/c/classes/FemModel.cpp
===================================================================
--- issm/trunk/src/c/classes/FemModel.cpp 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/c/classes/FemModel.cpp 2012-11-30 18:09:49 UTC (rev 14067)
@@ -19,6 +19,11 @@
#include "../modules/modules.h"
/*Object constructors and destructor*/
+/*FUNCTION FemModel::FemModel(char* filename) {{{*/
+FemModel::FemModel(char* filename){
+
+}
+/*}}}*/
/*FUNCTION FemModel::FemModel(int argc,char** argv){{{*/
FemModel::FemModel(int argc,char** argv,COMM incomm){
@@ -501,7 +506,9 @@
for(j=0;j<numnodes;j++) flags[j]=false;
/*Loop over elements that hold node number i*/
- _assert_(head_e[node->Sid()]!=-1 || head_l[node->Sid()]!=-1);
+ //if(head_e[node->Sid()]==-1 && head_l[node->Sid()]==-1){
+ // printf("[%i] vertex %i\n",IssmComm::GetRank(),node->Sid()+1);
+ //}
for(j=head_e[node->Sid()];j!=-1;j=next_e[j]){
offset=count2offset_e[j];
element=dynamic_cast<Element*>(elements->GetObjectByOffset(offset));
Modified: issm/trunk/src/c/classes/FemModel.h
===================================================================
--- issm/trunk/src/c/classes/FemModel.h 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/c/classes/FemModel.h 2012-11-30 18:09:49 UTC (rev 14067)
@@ -46,6 +46,7 @@
COMM comm; //communicator for this particular model
/*constructors, destructors: */
+ FemModel(char* filename);
FemModel(int argc,char** argv,COMM comm_init);
FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, const int solution_type,const int* analyses,const int nummodels);
~FemModel();
Modified: issm/trunk/src/c/classes/objects/Elements/Penta.cpp
===================================================================
--- issm/trunk/src/c/classes/objects/Elements/Penta.cpp 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/c/classes/objects/Elements/Penta.cpp 2012-11-30 18:09:49 UTC (rev 14067)
@@ -65,9 +65,9 @@
this->sid=penta_sid;
/*Build neighbors list*/
- if (xIsNan<IssmDouble>(iomodel->Data(MeshUpperelementsEnum)[index])) penta_elements_ids[1]=this->id; //upper penta is the same penta
+ if (xIsNan<IssmDouble>(iomodel->Data(MeshUpperelementsEnum)[index]) || iomodel->Data(MeshUpperelementsEnum)[index]==-1.) penta_elements_ids[1]=this->id; //upper penta is the same penta
else penta_elements_ids[1]=reCast<int,IssmDouble>((iomodel->Data(MeshUpperelementsEnum)[index]));
- if (xIsNan<IssmDouble>(iomodel->Data(MeshLowerelementsEnum)[index])) penta_elements_ids[0]=this->id; //lower penta is the same penta
+ if (xIsNan<IssmDouble>(iomodel->Data(MeshLowerelementsEnum)[index]) || iomodel->Data(MeshLowerelementsEnum)[index]==-1.) penta_elements_ids[0]=this->id; //lower penta is the same penta
else penta_elements_ids[0]=reCast<int,IssmDouble>((iomodel->Data(MeshLowerelementsEnum)[index]));
this->InitHookNeighbors(penta_elements_ids);
Modified: issm/trunk/src/c/modules/ElementConnectivityx/ElementConnectivityx.cpp
===================================================================
--- issm/trunk/src/c/modules/ElementConnectivityx/ElementConnectivityx.cpp 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/c/modules/ElementConnectivityx/ElementConnectivityx.cpp 2012-11-30 18:09:49 UTC (rev 14067)
@@ -14,61 +14,57 @@
#include "../../toolkits/toolkits.h"
#include "../../EnumDefinitions/EnumDefinitions.h"
-int hascommondedge(double* element1,double* element2);
+int hascommondedge(int* element1,int* element2);
-void ElementConnectivityx( double** pelementconnectivity, double* elements, int nel, double* nodeconnectivity, int nods, int width){
+void ElementConnectivityx(int** pelementconnectivity,int* elements, int nels,int* nodeconnectivity, int nods, int width){
int i,j,k,n;
/*intermediary: */
int maxels;
- double element;
- double connectedelement;
+ int connectedelement;
int connectedelementindex;
int node;
int index;
int num_elements;
- /*output: */
- double* elementconnectivity=NULL;
-
/*maxels: */
maxels=width-1;
+
/*Allocate connectivity: */
- elementconnectivity=xNewZeroInit<double>(nel*3);
+ int* elementconnectivity=xNewZeroInit<int>(nels*3);
/*Go through all elements, and for each element, go through its nodes, to get the neighbouring elements.
* Once we get the neighbouring elements, figure out if they share a segment with the current element. If so,
* plug them in the connectivity, unless they are already there.: */
+ for(n=0;n<nels;n++){
- for(n=0;n<nel;n++){
+ //element=n+1; //matlab indexing
- element=(double)(n+1); //matlab indexing
-
for(i=0;i<3;i++){
- node=(int)*(elements+n*3+i); //already matlab indexed, elements comes directly from the workspace.
+ node=elements[n*3+i]; //already matlab indexed, elements comes directly from the workspace.
index=node-1;
- num_elements=(int)*(nodeconnectivity+width*index+maxels); //retrieve number of elements already plugged into the connectivity of this node.
+ num_elements=nodeconnectivity[width*index+maxels]; //retrieve number of elements already plugged into the connectivity of this node.
for(j=0;j<num_elements;j++){
/*for each element connected to node, figure out if it has a commond edge with element: */
- connectedelement=*(nodeconnectivity+width*index+j);
- connectedelementindex=(int)(connectedelement-1); //go from matlab indexing to c indexing.
+ connectedelement=nodeconnectivity[width*index+j];
+ connectedelementindex=connectedelement-1; //go from matlab indexing to c indexing.
- if(hascommondedge(elements+n*3+0,elements+connectedelementindex*3+0)){
+ if(hascommondedge(&elements[n*3+0],&elements[connectedelementindex*3+0])){
/*Ok, this connected element has a commond edge with element, plug it into elementconnectivity, unless
*it is already there: */
for(k=0;k<3;k++){
- if (*(elementconnectivity+3*n+k)==0){
- *(elementconnectivity+3*n+k)=connectedelement;
+ if(elementconnectivity[3*n+k]==0){
+ elementconnectivity[3*n+k]=connectedelement;
break;
}
else{
- if(connectedelement==*(elementconnectivity+3*n+k))break;
+ if(connectedelement==elementconnectivity[3*n+k]) break;
}
}
}
@@ -80,17 +76,16 @@
*pelementconnectivity=elementconnectivity;
}
-int hascommondedge(double* element1,double* element2){
+int hascommondedge(int* el1,int* el2){
- int i,j;
- int count;
-
- count=0;
- for(i=0;i<3;i++){
- for(j=0;j<3;j++){
- if (*(element1+i)==*(element2+j))count++;
+ int count=0;
+ for(int i=0;i<3;i++){
+ for(int j=0;j<3;j++){
+ if(el1[i]==el2[j]) count++;
}
}
- if (count==2)return 1;
- else return 0;
+ if(count==2)
+ return 1;
+ else
+ return 0;
}
Modified: issm/trunk/src/c/modules/ElementConnectivityx/ElementConnectivityx.h
===================================================================
--- issm/trunk/src/c/modules/ElementConnectivityx/ElementConnectivityx.h 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/c/modules/ElementConnectivityx/ElementConnectivityx.h 2012-11-30 18:09:49 UTC (rev 14067)
@@ -6,6 +6,6 @@
#define _ELEMENTCONNECTIVITYX_H
/* local prototypes: */
-void ElementConnectivityx( double** pelementconnectivity, double* elements, int nel, double* nodeconnectivity, int nods, int width);
+void ElementConnectivityx(int** pelementconnectivity,int* elements,int nels,int* nodeconnectivity, int nods, int width);
#endif /* _ELEMENTCONNECTIVITYX_H */
Modified: issm/trunk/src/c/modules/Krigingx/Krigingx.cpp
===================================================================
--- issm/trunk/src/c/modules/Krigingx/Krigingx.cpp 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/c/modules/Krigingx/Krigingx.cpp 2012-11-30 18:09:49 UTC (rev 14067)
@@ -21,6 +21,7 @@
/*Intermediaries*/
int mindata,maxdata;
+ double dmindata,dmaxdata,dnumthreads; //FIXME (Options come as double but we want to retrive integers)
double radius;
char *output = NULL;
Variogram *variogram = NULL;
@@ -36,8 +37,9 @@
/*Get Variogram from Options*/
ProcessVariogram(&variogram,options);
options->Get(&radius,"searchradius",0.);
- options->Get(&mindata,"mindata",1);
- options->Get(&maxdata,"maxdata",50);
+ options->Get(&dmindata,"mindata",1.); mindata=int(dmindata);//FIXME (Options come as double but we want to retrive integers)
+ options->Get(&dmaxdata,"maxdata",50.); maxdata=int(dmaxdata);//FIXME (Options come as double but we want to retrive integers)
+ options->Get(&dnumthreads,"numthreads",double(num)); num=int(dnumthreads);//FIXME (Options come as double but we want to retrive integers)
/*Process observation dataset*/
observations=new Observations(obs_list,obs_x,obs_y,obs_length,options);
@@ -87,12 +89,12 @@
gate.observations = observations;
gate.predictions = predictions;
gate.error = error;
- gate.percent = xNewZeroInit<double>(num);
+ gate.numdone = xNewZeroInit<int>(num);
/*launch the thread manager with Krigingxt as a core: */
LaunchThread(NearestNeighbort,(void*)&gate,num);
_printLine_("\r interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%");
- xDelete<double>(gate.percent);
+ xDelete<int>(gate.numdone);
}
else if(strcmp(output,"idw")==0){ //Inverse distance weighting
double power;
@@ -108,15 +110,40 @@
gate.observations = observations;
gate.predictions = predictions;
gate.error = error;
- gate.percent = xNewZeroInit<double>(num);
+ gate.numdone = xNewZeroInit<int>(num);
gate.power = power;
/*launch the thread manager with Krigingxt as a core: */
LaunchThread(idwt,(void*)&gate,num);
_printLine_("\r interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%");
- xDelete<double>(gate.percent);
+ xDelete<int>(gate.numdone);
}
+ else if(strcmp(output,"v4")==0){ //Inverse distance weighting
+#if !defined(_HAVE_GSL_)
+ _error_("GSL is required for v4 interpolation");
+#endif
+ /*initialize thread parameters: */
+ gate.n_interp = n_interp;
+ gate.x_interp = x_interp;
+ gate.y_interp = y_interp;
+ gate.radius = radius;
+ gate.mindata = mindata;
+ gate.maxdata = maxdata;
+ gate.variogram = variogram;
+ gate.observations = observations;
+ gate.predictions = predictions;
+ gate.error = error;
+ gate.numdone = xNewZeroInit<int>(num);
+
+ /*launch the thread manager with Krigingxt as a core: */
+ LaunchThread(v4t,(void*)&gate,num);
+ _printLine_("\r interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%");
+ xDelete<int>(gate.numdone);
+ }
else if(strcmp(output,"prediction")==0){
+#if !defined(_HAVE_GSL_)
+ _error_("GSL is required for v4 interpolation");
+#endif
/*initialize thread parameters: */
gate.n_interp = n_interp;
@@ -129,12 +156,12 @@
gate.observations = observations;
gate.predictions = predictions;
gate.error = error;
- gate.percent = xNewZeroInit<double>(num);
+ gate.numdone = xNewZeroInit<int>(num);
/*launch the thread manager with Krigingxt as a core: */
LaunchThread(Krigingxt,(void*)&gate,num);
_printLine_("\r interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%");
- xDelete<double>(gate.percent);
+ xDelete<int>(gate.numdone);
}
else{
_error_("output '" << output << "' not supported yet");
@@ -175,20 +202,19 @@
Observations *observations = gate->observations;
double *predictions = gate->predictions;
double *error = gate->error;
- double *percent = gate->percent;
+ int *numdone = gate->numdone;
- /*Intermediaries*/
- double localpercent;
-
/*partition loop across threads: */
PartitionRange(&i0,&i1,n_interp,num_threads,my_thread);
for(int idx=i0;idx<i1;idx++){
/*Print info*/
- percent[my_thread]=double(idx-i0)/double(i1-i0)*100.;
- localpercent=percent[0];
- for(int i=1;i<num_threads;i++) localpercent=min(localpercent,percent[i]);
- if(my_thread==0) _printString_("\r interpolation progress: "<<setw(6)<<setprecision(2)<<localpercent<<"%");
+ numdone[my_thread]=idx-i0;
+ if(my_thread==0){
+ int alldone=numdone[0];
+ for(int i=1;i<num_threads;i++) alldone+=numdone[i];
+ _printString_("\r interpolation progress: "<<setw(6)<<setprecision(2)<<double(alldone)/double(n_interp)*100.<<"%");
+ }
/*Kriging interpolation*/
observations->InterpolationKriging(&predictions[idx],&error[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata,variogram);
@@ -223,21 +249,19 @@
Observations *observations = gate->observations;
double *predictions = gate->predictions;
double *error = gate->error;
- double *percent = gate->percent;
+ int *numdone = gate->numdone;
- /*Intermediaries*/
- int i;
- double localpercent;
-
/*partition loop across threads: */
PartitionRange(&i0,&i1,n_interp,num_threads,my_thread);
for(int idx=i0;idx<i1;idx++){
/*Print info*/
- percent[my_thread]=double(idx-i0)/double(i1-i0)*100.;
- localpercent=percent[0];
- for(i=1;i<num_threads;i++) localpercent=min(localpercent,percent[i]);
- if(my_thread==0) _printString_("\r interpolation progress: "<<setw(6)<<setprecision(2)<<localpercent<<"%");
+ numdone[my_thread]=idx-i0;
+ if(my_thread==0){
+ int alldone=numdone[0];
+ for(int i=1;i<num_threads;i++) alldone+=numdone[i];
+ _printString_("\r interpolation progress: "<<setw(6)<<setprecision(2)<<double(alldone)/double(n_interp)*100.<<"%");
+ }
observations->InterpolationNearestNeighbor(&predictions[idx],x_interp[idx],y_interp[idx],radius);
}
@@ -271,27 +295,71 @@
Observations *observations = gate->observations;
double *predictions = gate->predictions;
double *error = gate->error;
- double *percent = gate->percent;
+ int *numdone = gate->numdone;
double power = gate->power;
- /*Intermediaries*/
- double localpercent;
-
/*partition loop across threads: */
PartitionRange(&i0,&i1,n_interp,num_threads,my_thread);
for(int idx=i0;idx<i1;idx++){
/*Print info*/
- percent[my_thread]=double(idx-i0)/double(i1-i0)*100.;
- localpercent=percent[0];
- for(int i=1;i<num_threads;i++) localpercent=min(localpercent,percent[i]);
- if(my_thread==0) _printString_("\r interpolation progress: "<<setw(6)<<setprecision(2)<<localpercent<<"%");
+ numdone[my_thread]=idx-i0;
+ if(my_thread==0){
+ int alldone=numdone[0];
+ for(int i=1;i<num_threads;i++) alldone+=numdone[i];
+ _printString_("\r interpolation progress: "<<setw(6)<<setprecision(2)<<double(alldone)/double(n_interp)*100.<<"%");
+ }
observations->InterpolationIDW(&predictions[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata,power);
}
return NULL;
}/*}}}*/
+/*FUNCTION v4t{{{*/
+void* v4t(void* vpthread_handle){
+ /*gate variables :*/
+ KrigingxThreadStruct *gate = NULL;
+ pthread_handle *handle = NULL;
+ int my_thread;
+ int num_threads;
+ int i0,i1;
+
+ /*recover handle and gate: */
+ handle = (pthread_handle*)vpthread_handle;
+ gate = (KrigingxThreadStruct*)handle->gate;
+ my_thread = handle->id;
+ num_threads = handle->num;
+
+ /*recover parameters :*/
+ int n_interp = gate->n_interp;
+ double *x_interp = gate->x_interp;
+ double *y_interp = gate->y_interp;
+ double radius = gate->radius;
+ int mindata = gate->mindata;
+ int maxdata = gate->maxdata;
+ Variogram *variogram = gate->variogram;
+ Observations *observations = gate->observations;
+ double *predictions = gate->predictions;
+ double *error = gate->error;
+ int *numdone = gate->numdone;
+
+ /*partition loop across threads: */
+ PartitionRange(&i0,&i1,n_interp,num_threads,my_thread);
+ for(int idx=i0;idx<i1;idx++){
+
+ /*Print info*/
+ numdone[my_thread]=idx-i0;
+ if(my_thread==0){
+ int alldone=numdone[0];
+ for(int i=1;i<num_threads;i++) alldone+=numdone[i];
+ _printString_("\r interpolation progress: "<<setw(6)<<setprecision(2)<<double(alldone)/double(n_interp)*100.<<"%");
+ }
+
+ observations->InterpolationV4(&predictions[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata);
+ }
+ return NULL;
+}/*}}}*/
+
void ProcessVariogram(Variogram **pvariogram,Options* options){/*{{{*/
/*Intermediaries*/
Modified: issm/trunk/src/c/modules/Krigingx/Krigingx.h
===================================================================
--- issm/trunk/src/c/modules/Krigingx/Krigingx.h 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/c/modules/Krigingx/Krigingx.h 2012-11-30 18:09:49 UTC (rev 14067)
@@ -27,11 +27,12 @@
Observations *observations;
double *predictions;
double *error;
- double *percent;
+ int *numdone;
double power;//for idw
}KrigingxThreadStruct;
void* Krigingxt(void*);
void* NearestNeighbort(void*);
void* idwt(void*);
+void* v4t(void*);
#endif /* _KRIGINGX_H */
Modified: issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.cpp
===================================================================
--- issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.cpp 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.cpp 2012-11-30 18:09:49 UTC (rev 14067)
@@ -17,7 +17,7 @@
#include "../../toolkits/toolkits.h"
#include "../../EnumDefinitions/EnumDefinitions.h"
-void NodeConnectivityx( double** pconnectivity, int* pwidth, double* elements, int nel, int nods){
+void NodeConnectivityx(int** pconnectivity,int* pwidth,int* elements, int nels, int nods){
int i,j,n;
const int maxels=25;
@@ -28,26 +28,23 @@
int index;
int num_elements;
int already_plugged=0;
- double element;
+ int element;
- /*output: */
- double* connectivity=NULL;
-
/*Allocate connectivity: */
- connectivity=xNewZeroInit<double>(nods*width);
+ int* connectivity=xNewZeroInit<int>(nods*width);
/*Go through all elements, and for each elements, plug into the connectivity, all the nodes.
* If nodes are already plugged into the connectivity, skip them.: */
- for(n=0;n<nel;n++){
+ for(n=0;n<nels;n++){
- element=(double)(n+1); //matlab indexing
+ element=n+1; //matlab indexing
for(i=0;i<3;i++){
- node=(int)*(elements+n*3+i); //already matlab indexed, elements comes directly from the workspace.
+ node=elements[n*3+i]; //already matlab indexed, elements comes directly from the workspace.
index=node-1;
- num_elements=(int)*(connectivity+width*index+maxels); //retrieve number of elements already plugged into the connectivity of this node.
+ num_elements=connectivity[width*index+maxels]; //retrieve number of elements already plugged into the connectivity of this node.
already_plugged=0;
for(j=0;j<num_elements;j++){
@@ -59,8 +56,8 @@
if(already_plugged)break;
/*this elements is not yet plugged into the connectivity for this node, do it, and increase counter: */
- *(connectivity+width*index+num_elements)=element;
- *(connectivity+width*index+maxels)=(double)(num_elements+1);
+ connectivity[width*index+num_elements]=element;
+ connectivity[width*index+maxels]=num_elements+1;
}
}
@@ -68,7 +65,8 @@
/*Last check: is the number of elements on last column of the connectivity superior to maxels? If so, then error out and
* warn the user to increase the connectivity width: */
for(i=0;i<nods;i++){
- if (*(connectivity+width*i+maxels)>maxels)_error_("max connectivity width reached (" << *(connectivity+width*i+maxels) << ")! increase width of connectivity table");
+ if (*(connectivity+width*i+maxels)>maxels)
+ _error_("max connectivity width reached (" << *(connectivity+width*i+maxels) << ")! increase width of connectivity table");
}
/*Assign output pointers: */
Modified: issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.h
===================================================================
--- issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.h 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.h 2012-11-30 18:09:49 UTC (rev 14067)
@@ -6,6 +6,6 @@
#define _NODECONNECTIVITYX_H
/* local prototypes: */
-void NodeConnectivityx( double** pconnectivity, int* pwidth,double* elements, int nel, int nods);
+void NodeConnectivityx(int** pconnectivity,int* pwidth,int* elements,int nels, int nods);
#endif /* _NODECONNECTIVITYX_H */
Modified: issm/trunk/src/c/shared/Threads/LaunchThread.cpp
===================================================================
--- issm/trunk/src/c/shared/Threads/LaunchThread.cpp 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/c/shared/Threads/LaunchThread.cpp 2012-11-30 18:09:49 UTC (rev 14067)
@@ -31,6 +31,10 @@
pthread_t *threads = NULL;
pthread_handle *handles = NULL;
+ /*check number of threads*/
+ if(num_threads<1) _error_("number of threads must be at least 1");
+ if(num_threads>2000) _error_("number of threads seems to be too high ("<<num_threads<<">2000)");
+
/*dynamically allocate: */
threads=xNew<pthread_t>(num_threads);
handles=xNew<pthread_handle>(num_threads);
Modified: issm/trunk/src/dox/issm.dox
===================================================================
--- issm/trunk/src/dox/issm.dox 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/dox/issm.dox 2012-11-30 18:09:49 UTC (rev 14067)
@@ -45,25 +45,25 @@
<th bgcolor=#7AA9DD style="text-align:left;">Language</th><th bgcolor=#7AA9DD style="text-align:right;">files</th><th bgcolor=#7AA9DD style="text-align:right;">blank</th><th bgcolor=#7AA9DD style="text-align:right;">comment</th><th bgcolor=#7AA9DD style="text-align:right;">code</th><th bgcolor=#7AA9DD style="text-align:right;">Total</th>
</th>
<tr>
-<th bgcolor=#FFFFFF style="text-align:left;"> C++ </th><td bgcolor=#FFFFFF style="text-align:right;">505</td><td bgcolor=#FFFFFF style="text-align:right;">14316</td><td bgcolor=#FFFFFF style="text-align:right;">16564</td><td bgcolor=#FFFFFF style="text-align:right;">56522</td><td bgcolor=#FFFFFF style="text-align:right;">87402</td>
+<th bgcolor=#FFFFFF style="text-align:left;"> C++ </th><td bgcolor=#FFFFFF style="text-align:right;">469</td><td bgcolor=#FFFFFF style="text-align:right;">14087</td><td bgcolor=#FFFFFF style="text-align:right;">16669</td><td bgcolor=#FFFFFF style="text-align:right;">56706</td><td bgcolor=#FFFFFF style="text-align:right;">87462</td>
</tr>
<tr>
-<th bgcolor=#C6E2FF style="text-align:left;"> MATLAB </th><td bgcolor=#C6E2FF style="text-align:right;">913</td><td bgcolor=#C6E2FF style="text-align:right;">6625</td><td bgcolor=#C6E2FF style="text-align:right;">13054</td><td bgcolor=#C6E2FF style="text-align:right;">30160</td><td bgcolor=#C6E2FF style="text-align:right;">49839</td>
+<th bgcolor=#C6E2FF style="text-align:left;"> MATLAB </th><td bgcolor=#C6E2FF style="text-align:right;">945</td><td bgcolor=#C6E2FF style="text-align:right;">6583</td><td bgcolor=#C6E2FF style="text-align:right;">13346</td><td bgcolor=#C6E2FF style="text-align:right;">30647</td><td bgcolor=#C6E2FF style="text-align:right;">50576</td>
</tr>
<tr>
-<th bgcolor=#FFFFFF style="text-align:left;"> C/C++ Header </th><td bgcolor=#FFFFFF style="text-align:right;">388</td><td bgcolor=#FFFFFF style="text-align:right;">3069</td><td bgcolor=#FFFFFF style="text-align:right;">2983</td><td bgcolor=#FFFFFF style="text-align:right;">11855</td><td bgcolor=#FFFFFF style="text-align:right;">17907</td>
+<th bgcolor=#FFFFFF style="text-align:left;"> C/C++ Header </th><td bgcolor=#FFFFFF style="text-align:right;">362</td><td bgcolor=#FFFFFF style="text-align:right;">2772</td><td bgcolor=#FFFFFF style="text-align:right;">2894</td><td bgcolor=#FFFFFF style="text-align:right;">11914</td><td bgcolor=#FFFFFF style="text-align:right;">17580</td>
</tr>
<tr>
-<th bgcolor=#C6E2FF style="text-align:left;"> m4 </th><td bgcolor=#C6E2FF style="text-align:right;">6</td><td bgcolor=#C6E2FF style="text-align:right;">1052</td><td bgcolor=#C6E2FF style="text-align:right;">79</td><td bgcolor=#C6E2FF style="text-align:right;">8993</td><td bgcolor=#C6E2FF style="text-align:right;">10124</td>
+<th bgcolor=#C6E2FF style="text-align:left;"> Python </th><td bgcolor=#C6E2FF style="text-align:right;">109</td><td bgcolor=#C6E2FF style="text-align:right;">3361</td><td bgcolor=#C6E2FF style="text-align:right;">4763</td><td bgcolor=#C6E2FF style="text-align:right;">7050</td><td bgcolor=#C6E2FF style="text-align:right;">15174</td>
</tr>
<tr>
-<th bgcolor=#FFFFFF style="text-align:left;"> Python </th><td bgcolor=#FFFFFF style="text-align:right;">81</td><td bgcolor=#FFFFFF style="text-align:right;">2677</td><td bgcolor=#FFFFFF style="text-align:right;">3936</td><td bgcolor=#FFFFFF style="text-align:right;">5169</td><td bgcolor=#FFFFFF style="text-align:right;">11782</td>
+<th bgcolor=#FFFFFF style="text-align:left;"> m4 </th><td bgcolor=#FFFFFF style="text-align:right;">1</td><td bgcolor=#FFFFFF style="text-align:right;">187</td><td bgcolor=#FFFFFF style="text-align:right;">5</td><td bgcolor=#FFFFFF style="text-align:right;">1247</td><td bgcolor=#FFFFFF style="text-align:right;">1439</td>
</tr>
<tr>
-<th bgcolor=#C6E2FF style="text-align:left;"> Objective C </th><td bgcolor=#C6E2FF style="text-align:right;">9</td><td bgcolor=#C6E2FF style="text-align:right;">98</td><td bgcolor=#C6E2FF style="text-align:right;">0</td><td bgcolor=#C6E2FF style="text-align:right;">370</td><td bgcolor=#C6E2FF style="text-align:right;">468</td>
+<th bgcolor=#C6E2FF style="text-align:left;"> Objective C </th><td bgcolor=#C6E2FF style="text-align:right;">9</td><td bgcolor=#C6E2FF style="text-align:right;">96</td><td bgcolor=#C6E2FF style="text-align:right;">0</td><td bgcolor=#C6E2FF style="text-align:right;">370</td><td bgcolor=#C6E2FF style="text-align:right;">466</td>
</tr>
<tr>
-<th bgcolor=#FFFFFF style="text-align:left;"> Bourne Shell </th><td bgcolor=#FFFFFF style="text-align:right;">5</td><td bgcolor=#FFFFFF style="text-align:right;">58</td><td bgcolor=#FFFFFF style="text-align:right;">81</td><td bgcolor=#FFFFFF style="text-align:right;">268</td><td bgcolor=#FFFFFF style="text-align:right;">407</td>
+<th bgcolor=#FFFFFF style="text-align:left;"> Bourne Shell </th><td bgcolor=#FFFFFF style="text-align:right;">5</td><td bgcolor=#FFFFFF style="text-align:right;">57</td><td bgcolor=#FFFFFF style="text-align:right;">79</td><td bgcolor=#FFFFFF style="text-align:right;">256</td><td bgcolor=#FFFFFF style="text-align:right;">392</td>
</tr>
<tr>
<th bgcolor=#C6E2FF style="text-align:left;"> Perl </th><td bgcolor=#C6E2FF style="text-align:right;">3</td><td bgcolor=#C6E2FF style="text-align:right;">21</td><td bgcolor=#C6E2FF style="text-align:right;">23</td><td bgcolor=#C6E2FF style="text-align:right;">240</td><td bgcolor=#C6E2FF style="text-align:right;">284</td>
@@ -75,7 +75,7 @@
<th bgcolor=#C6E2FF style="text-align:left;"> C </th><td bgcolor=#C6E2FF style="text-align:right;">1</td><td bgcolor=#C6E2FF style="text-align:right;">2</td><td bgcolor=#C6E2FF style="text-align:right;">3</td><td bgcolor=#C6E2FF style="text-align:right;">6</td><td bgcolor=#C6E2FF style="text-align:right;">11</td>
</tr>
<tr>
-<th bgcolor=#FFFFFF style="text-align:left;"> SUM: </th><td bgcolor=#FFFFFF style="text-align:right;">1912</td><td bgcolor=#FFFFFF style="text-align:right;">27936</td><td bgcolor=#FFFFFF style="text-align:right;">36727</td><td bgcolor=#FFFFFF style="text-align:right;">113707</td><td bgcolor=#FFFFFF style="text-align:right;">178370</td>
+<th bgcolor=#FFFFFF style="text-align:left;"> SUM: </th><td bgcolor=#FFFFFF style="text-align:right;">1905</td><td bgcolor=#FFFFFF style="text-align:right;">27184</td><td bgcolor=#FFFFFF style="text-align:right;">37786</td><td bgcolor=#FFFFFF style="text-align:right;">108560</td><td bgcolor=#FFFFFF style="text-align:right;">173530</td>
</tr>
</table>
Modified: issm/trunk/src/m/boundaryconditions/SetIceShelfBC.py
===================================================================
--- issm/trunk/src/m/boundaryconditions/SetIceShelfBC.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/boundaryconditions/SetIceShelfBC.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -25,9 +25,9 @@
if not os.path.exists(icefrontfile):
raise IOError("SetIceShelfBC error message: ice front file '%s' not found." % icefrontfile)
[nodeinsideicefront,dum]=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),icefrontfile,'node',2)
- nodeonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1)).astype(float)
+ nodeonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1))
else:
- nodeonicefront=numpy.zeros((md.mesh.numberofvertices))
+ nodeonicefront=numpy.zeros((md.mesh.numberofvertices),bool)
# pos=find(md.mesh.vertexonboundary & ~nodeonicefront);
pos=numpy.nonzero(numpy.logical_and(md.mesh.vertexonboundary,numpy.logical_not(nodeonicefront)))[0]
@@ -55,20 +55,20 @@
#segment on Ice Front
#segment on Neumann (Ice Front)
# pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2)));
- pos=numpy.nonzero(numpy.logical_or(nodeonicefront[md.mesh.segments[:,0].astype(int)-1],nodeonicefront[md.mesh.segments[:,1].astype(int)-1]))[0]
+ pos=numpy.nonzero(numpy.logical_or(nodeonicefront[md.mesh.segments[:,0]-1],nodeonicefront[md.mesh.segments[:,1]-1]))[0]
if md.mesh.dimension==2:
pressureload=md.mesh.segments[pos,:]
elif md.mesh.dimension==3:
# pressureload_layer1=[md.mesh.segments(pos,1:2) md.mesh.segments(pos,2)+md.mesh.numberofvertices2d md.mesh.segments(pos,1)+md.mesh.numberofvertices2d md.mesh.segments(pos,3)];
pressureload_layer1=numpy.hstack((md.mesh.segments[pos,0:2],md.mesh.segments[pos,1]+md.mesh.numberofvertices2d,md.mesh.segments[pos,0]+md.mesh.numberofvertices2d,md.mesh.segments[pos,2]))
- pressureload=numpy.zeros((0,5))
+ pressureload=numpy.zeros((0,5),int)
for i in xrange(1,md.mesh.numberoflayers):
# pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ];
pressureload=numpy.vstack((pressureload,numpy.hstack((pressureload_layer1[:,0:4]+(i-1)*md.mesh.numberofvertices2d,pressureload_layer1[:,4]+(i-1)*md.mesh.numberofelements2d))))
#Add water or air enum depending on the element
# pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end))];
- pressureload=numpy.hstack((pressureload,1.*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape(-1,1)))
+ pressureload=numpy.hstack((pressureload,1*md.mask.elementonfloatingice[pressureload[:,-1]-1].reshape(-1,1)))
#plug onto model
md.diagnostic.icefront=pressureload
Modified: issm/trunk/src/m/boundaryconditions/SetMarineIceSheetBC.py
===================================================================
--- issm/trunk/src/m/boundaryconditions/SetMarineIceSheetBC.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/boundaryconditions/SetMarineIceSheetBC.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -27,12 +27,12 @@
if not os.path.exists(icefrontfile):
raise IOError("SetMarineIceSheetBC error message: ice front file '%s' not found." % icefrontfile)
[nodeinsideicefront,dum]=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),icefrontfile,'node',2)
- vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1)).astype(float)
+ vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1))
else:
#Guess where the ice front is
- vertexonfloatingice=numpy.zeros((md.mesh.numberofvertices))
- vertexonfloatingice[md.mesh.elements[numpy.nonzero(md.mask.elementonfloatingice),:].astype(int)-1]=1
- vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,vertexonfloatingice).astype(float)
+ vertexonfloatingice=numpy.zeros((md.mesh.numberofvertices),bool)
+ vertexonfloatingice[md.mesh.elements[numpy.nonzero(md.mask.elementonfloatingice),:]-1]=True
+ vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,vertexonfloatingice)
# pos=find(md.mesh.vertexonboundary & ~vertexonicefront);
pos=numpy.nonzero(numpy.logical_and(md.mesh.vertexonboundary,numpy.logical_not(vertexonicefront)))[0]
@@ -61,20 +61,20 @@
#segment on Neumann (Ice Front)
# pos=find(vertexonicefront(md.mesh.segments(:,1)) | vertexonicefront(md.mesh.segments(:,2)));
- pos=numpy.nonzero(numpy.logical_or(vertexonicefront[md.mesh.segments[:,0].astype(int)-1],vertexonicefront[md.mesh.segments[:,1].astype(int)-1]))[0]
+ pos=numpy.nonzero(numpy.logical_or(vertexonicefront[md.mesh.segments[:,0]-1],vertexonicefront[md.mesh.segments[:,1]-1]))[0]
if md.mesh.dimension==2:
pressureload=md.mesh.segments[pos,:]
elif md.mesh.dimension==3:
# pressureload_layer1=[md.mesh.segments(pos,1:2) md.mesh.segments(pos,2)+md.mesh.numberofvertices2d md.mesh.segments(pos,1)+md.mesh.numberofvertices2d md.mesh.segments(pos,3)];
pressureload_layer1=numpy.hstack((md.mesh.segments[pos,0:2],md.mesh.segments[pos,1]+md.mesh.numberofvertices2d,md.mesh.segments[pos,0]+md.mesh.numberofvertices2d,md.mesh.segments[pos,2]))
- pressureload=numpy.zeros((0,5))
+ pressureload=numpy.zeros((0,5),int)
for i in xrange(1,md.mesh.numberoflayers):
# pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ];
pressureload=numpy.vstack((pressureload,numpy.hstack((pressureload_layer1[:,0:4]+(i-1)*md.mesh.numberofvertices2d,pressureload_layer1[:,4]+(i-1)*md.mesh.numberofelements2d))))
#Add water or air enum depending on the element
# pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end))+ 0*md.mask.elementongroundedice(pressureload(:,end))];
- pressureload=numpy.hstack((pressureload,1.*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape(-1,1)+0.*md.mask.elementongroundedice[pressureload[:,-1].astype('int')-1].reshape(-1,1)))
+ pressureload=numpy.hstack((pressureload,1*md.mask.elementonfloatingice[pressureload[:,-1]-1].reshape(-1,1)+0*md.mask.elementongroundedice[pressureload[:,-1]-1].reshape(-1,1)))
#plug onto model
md.diagnostic.icefront=pressureload
Modified: issm/trunk/src/m/classes/clusters/greenplanet.m
===================================================================
--- issm/trunk/src/m/classes/clusters/greenplanet.m 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/classes/clusters/greenplanet.m 2012-11-30 18:09:49 UTC (rev 14067)
@@ -13,7 +13,7 @@
numnodes=20;
cpuspernode=8;
port=8000;
- queue='rignot';
+ queue='c6145';
codepath='';
executionpath='';
interactive=0;
@@ -49,7 +49,7 @@
%}}}
function md = checkconsistency(cluster,md,solution,analyses) % {{{
- available_queues={'rignot','default'};
+ available_queues={'c6145','default'};
queue_requirements_time=[Inf Inf];
queue_requirements_np=[80 80];
Modified: issm/trunk/src/m/classes/flowequation.py
===================================================================
--- issm/trunk/src/m/classes/flowequation.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/classes/flowequation.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -18,10 +18,10 @@
def __init__(self):
# {{{ Properties
- self.ismacayealpattyn = 0;
- self.ishutter = 0;
- self.isl1l2 = 0;
- self.isstokes = 0;
+ self.ismacayealpattyn = 0
+ self.ishutter = 0
+ self.isl1l2 = 0
+ self.isstokes = 0
self.vertex_equation = float('NaN')
self.element_equation = float('NaN')
self.bordermacayeal = float('NaN')
@@ -36,15 +36,15 @@
# {{{ Display
string=' flow equation parameters:'
- string="%s\n\n%s"%(string,fielddisplay(self,'ismacayealpattyn','is the macayeal or pattyn approximation used ?'))
- string="%s\n%s"%(string,fielddisplay(self,'ishutter','is the shallow ice approximation used ?'))
- string="%s\n%s"%(string,fielddisplay(self,'isl1l2','are l1l2 equations used ?'))
- string="%s\n%s"%(string,fielddisplay(self,'isstokes','are the Full-Stokes equations used ?'))
- string="%s\n%s"%(string,fielddisplay(self,'vertex_equation','flow equation for each vertex'))
- string="%s\n%s"%(string,fielddisplay(self,'element_equation','flow equation for each element'))
- string="%s\n%s"%(string,fielddisplay(self,'bordermacayeal','vertices on MacAyeal''s border (for tiling)'))
- string="%s\n%s"%(string,fielddisplay(self,'borderpattyn','vertices on Pattyn''s border (for tiling)'))
- string="%s\n%s"%(string,fielddisplay(self,'borderstokes','vertices on Stokes'' border (for tiling)'))
+ string="%s\n\n%s"%(string,fielddisplay(self,'ismacayealpattyn',"is the macayeal or pattyn approximation used ?"))
+ string="%s\n%s"%(string,fielddisplay(self,'ishutter',"is the shallow ice approximation used ?"))
+ string="%s\n%s"%(string,fielddisplay(self,'isl1l2',"are l1l2 equations used ?"))
+ string="%s\n%s"%(string,fielddisplay(self,'isstokes',"are the Full-Stokes equations used ?"))
+ string="%s\n%s"%(string,fielddisplay(self,'vertex_equation',"flow equation for each vertex"))
+ string="%s\n%s"%(string,fielddisplay(self,'element_equation',"flow equation for each element"))
+ string="%s\n%s"%(string,fielddisplay(self,'bordermacayeal',"vertices on MacAyeal's border (for tiling)"))
+ string="%s\n%s"%(string,fielddisplay(self,'borderpattyn',"vertices on Pattyn's border (for tiling)"))
+ string="%s\n%s"%(string,fielddisplay(self,'borderstokes',"vertices on Stokes' border (for tiling)"))
return string
#}}}
Modified: issm/trunk/src/m/classes/initialization.py
===================================================================
--- issm/trunk/src/m/classes/initialization.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/classes/initialization.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -69,8 +69,8 @@
md = checkfield(md,'initialization.vx','NaN',1,'size',[md.mesh.numberofvertices])
md = checkfield(md,'initialization.vy','NaN',1,'size',[md.mesh.numberofvertices])
#Triangle with zero velocity
- if numpy.any(numpy.logical_and(numpy.sum(numpy.abs(md.initialization.vx[md.mesh.elements.astype(int)-1]),axis=1)==0,\
- numpy.sum(numpy.abs(md.initialization.vy[md.mesh.elements.astype(int)-1]),axis=1)==0)):
+ if numpy.any(numpy.logical_and(numpy.sum(numpy.abs(md.initialization.vx[md.mesh.elements-1]),axis=1)==0,\
+ numpy.sum(numpy.abs(md.initialization.vy[md.mesh.elements-1]),axis=1)==0)):
md.checkmessage("at least one triangle has all its vertices with a zero velocity")
if ThermalAnalysisEnum() in analyses:
md = checkfield(md,'initialization.vx','NaN',1,'size',[md.mesh.numberofvertices])
Modified: issm/trunk/src/m/classes/mask.py
===================================================================
--- issm/trunk/src/m/classes/mask.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/classes/mask.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -32,7 +32,7 @@
string="";
string="%s\n%s"%(string,fielddisplay(self,"elementonfloatingice","element on floating ice flags list"))
string="%s\n%s"%(string,fielddisplay(self,"vertexonfloatingice","vertex on floating ice flags list"))
- string="%s\n%s"%(string,fielddisplay(self,"elementongroundedice","element on grounded ice list"))
+ string="%s\n%s"%(string,fielddisplay(self,"elementongroundedice","element on grounded ice list"))
string="%s\n%s"%(string,fielddisplay(self,"vertexongroundedice","vertex on grounded ice flags list"))
string="%s\n%s"%(string,fielddisplay(self,"elementonwater","element on water flags list"))
string="%s\n%s"%(string,fielddisplay(self,"vertexonwater","vertex on water flags list"))
Modified: issm/trunk/src/m/classes/mesh.py
===================================================================
--- issm/trunk/src/m/classes/mesh.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/classes/mesh.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -91,10 +91,10 @@
string="%s\n%s"%(string,fielddisplay(self,"elementonbed","lower elements flags list"))
string="%s\n%s"%(string,fielddisplay(self,"vertexonsurface","upper vertices flags list"))
string="%s\n%s"%(string,fielddisplay(self,"elementonsurface","upper elements flags list"))
- string="%s\n%s"%(string,fielddisplay(self,"uppervertex","upper vertex list (NaN for vertex on the upper surface)"))
- string="%s\n%s"%(string,fielddisplay(self,"upperelements","upper element list (NaN for element on the upper layer)"))
- string="%s\n%s"%(string,fielddisplay(self,"lowervertex","lower vertex list (NaN for vertex on the lower surface)"))
- string="%s\n%s"%(string,fielddisplay(self,"lowerelements","lower element list (NaN for element on the lower layer"))
+ string="%s\n%s"%(string,fielddisplay(self,"uppervertex","upper vertex list (-1 for vertex on the upper surface)"))
+ string="%s\n%s"%(string,fielddisplay(self,"upperelements","upper element list (-1 for element on the upper layer)"))
+ string="%s\n%s"%(string,fielddisplay(self,"lowervertex","lower vertex list (-1 for vertex on the lower surface)"))
+ string="%s\n%s"%(string,fielddisplay(self,"lowerelements","lower element list (-1 for element on the lower layer)"))
string="%s\n%s"%(string,fielddisplay(self,"vertexonboundary","vertices on the boundary of the domain flag list"))
string="%s\n%s"%(string,fielddisplay(self,"segments","edges on domain boundary (vertex1 vertex2 element)"))
string="%s\n%s"%(string,fielddisplay(self,"segmentmarkers","number associated to each segment"))
Modified: issm/trunk/src/m/classes/model/model.py
===================================================================
--- issm/trunk/src/m/classes/model/model.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/classes/model/model.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -1,6 +1,7 @@
#module imports {{{
import numpy
import copy
+import sys
from mesh import mesh
from mask import mask
from geometry import geometry
@@ -206,15 +207,15 @@
#kick out all elements with 3 dirichlets
spc_elem=numpy.nonzero(numpy.logical_not(flag_elem))[0]
- spc_node=numpy.unique(md1.mesh.elements[spc_elem,:]).astype(int)-1
+ spc_node=numpy.unique(md1.mesh.elements[spc_elem,:])-1
flag=numpy.ones(md1.mesh.numberofvertices)
flag[spc_node]=0
- pos=numpy.nonzero(numpy.logical_not(numpy.sum(flag[md1.mesh.elements.astype(int)-1],axis=1)))[0]
+ pos=numpy.nonzero(numpy.logical_not(numpy.sum(flag[md1.mesh.elements-1],axis=1)))[0]
flag_elem[pos]=0
#extracted elements and nodes lists
pos_elem=numpy.nonzero(flag_elem)[0]
- pos_node=numpy.unique(md1.mesh.elements[pos_elem,:]).astype(int)-1
+ pos_node=numpy.unique(md1.mesh.elements[pos_elem,:])-1
#keep track of some fields
numberofvertices1=md1.mesh.numberofvertices
@@ -225,21 +226,21 @@
flag_node[pos_node]=1
#Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements)
- Pelem=numpy.zeros(numberofelements1)
+ Pelem=numpy.zeros(numberofelements1,int)
Pelem[pos_elem]=numpy.arange(1,numberofelements2+1)
- Pnode=numpy.zeros(numberofvertices1)
+ Pnode=numpy.zeros(numberofvertices1,int)
Pnode[pos_node]=numpy.arange(1,numberofvertices2+1)
#renumber the elements (some node won't exist anymore)
elements_1=copy.deepcopy(md1.mesh.elements)
elements_2=elements_1[pos_elem,:]
- elements_2[:,0]=Pnode[elements_2[:,0].astype(int)-1]
- elements_2[:,1]=Pnode[elements_2[:,1].astype(int)-1]
- elements_2[:,2]=Pnode[elements_2[:,2].astype(int)-1]
+ elements_2[:,0]=Pnode[elements_2[:,0]-1]
+ elements_2[:,1]=Pnode[elements_2[:,1]-1]
+ elements_2[:,2]=Pnode[elements_2[:,2]-1]
if md1.mesh.dimension==3:
- elements_2[:,3]=Pnode[elements_2[:,3].astype(int)-1]
- elements_2[:,4]=Pnode[elements_2[:,4].astype(int)-1]
- elements_2[:,5]=Pnode[elements_2[:,5].astype(int)-1]
+ elements_2[:,3]=Pnode[elements_2[:,3]-1]
+ elements_2[:,4]=Pnode[elements_2[:,4]-1]
+ elements_2[:,5]=Pnode[elements_2[:,5]-1]
#OK, now create the new model!
@@ -290,20 +291,20 @@
#mesh.uppervertex mesh.lowervertex
if md1.mesh.dimension==3:
md2.mesh.uppervertex=md1.mesh.uppervertex[pos_node]
- pos=numpy.nonzero(numpy.logical_not(numpy.isnan(md2.mesh.uppervertex)))[0]
- md2.mesh.uppervertex[pos]=Pnode[md2.mesh.uppervertex[pos].astype(int)-1]
+ pos=numpy.nonzero(numpy.logical_not(md2.mesh.uppervertex==-1))[0]
+ md2.mesh.uppervertex[pos]=Pnode[md2.mesh.uppervertex[pos]-1]
md2.mesh.lowervertex=md1.mesh.lowervertex[pos_node]
- pos=numpy.nonzero(numpy.logical_not(numpy.isnan(md2.mesh.lowervertex)))[0]
- md2.mesh.lowervertex[pos]=Pnode[md2.mesh.lowervertex[pos].astype(int)-1]
+ pos=numpy.nonzero(numpy.logical_not(md2.mesh.lowervertex==-1))[0]
+ md2.mesh.lowervertex[pos]=Pnode[md2.mesh.lowervertex[pos]-1]
md2.mesh.upperelements=md1.mesh.upperelements[pos_elem]
- pos=numpy.nonzero(numpy.logical_not(numpy.isnan(md2.mesh.upperelements)))[0]
- md2.mesh.upperelements[pos]=Pelem[md2.mesh.upperelements[pos].astype(int)-1]
+ pos=numpy.nonzero(numpy.logical_not(md2.mesh.upperelements==-1))[0]
+ md2.mesh.upperelements[pos]=Pelem[md2.mesh.upperelements[pos]-1]
md2.mesh.lowerelements=md1.mesh.lowerelements[pos_elem]
- pos=numpy.nonzero(numpy.logical_not(numpy.isnan(md2.mesh.lowerelements)))[0]
- md2.mesh.lowerelements[pos]=Pelem[md2.mesh.lowerelements[pos].astype(int)-1]
+ pos=numpy.nonzero(numpy.logical_not(md2.mesh.lowerelements==-1))[0]
+ md2.mesh.lowerelements[pos]=Pelem[md2.mesh.lowerelements[pos]-1]
#Initial 2d mesh
if md1.mesh.dimension==3:
@@ -315,21 +316,21 @@
md2.mesh.numberofelements2d=numpy.size(pos_elem_2d)
md2.mesh.numberofvertices2d=numpy.size(pos_node_2d)
md2.mesh.elements2d=md1.mesh.elements2d[pos_elem_2d,:]
- md2.mesh.elements2d[:,0]=Pnode[md2.mesh.elements2d[:,0].astype(int)-1]
- md2.mesh.elements2d[:,1]=Pnode[md2.mesh.elements2d[:,1].astype(int)-1]
- md2.mesh.elements2d[:,2]=Pnode[md2.mesh.elements2d[:,2].astype(int)-1]
+ md2.mesh.elements2d[:,0]=Pnode[md2.mesh.elements2d[:,0]-1]
+ md2.mesh.elements2d[:,1]=Pnode[md2.mesh.elements2d[:,1]-1]
+ md2.mesh.elements2d[:,2]=Pnode[md2.mesh.elements2d[:,2]-1]
md2.mesh.x2d=md1.mesh.x[pos_node_2d]
md2.mesh.y2d=md1.mesh.y[pos_node_2d]
#Edges
- if len(numpy.shape(md2.mesh.edges))>1 and numpy.size(md2.mesh.edges,axis=1)>1: #do not use ~isnan because there are some NaNs...
+ if numpy.ndim(md2.mesh.edges)>1 and numpy.size(md2.mesh.edges,axis=1)>1: #do not use ~isnan because there are some NaNs...
#renumber first two columns
pos=numpy.nonzero(md2.mesh.edges[:,3]!=-1)[0]
- md2.mesh.edges[: ,0]=Pnode[md2.mesh.edges[:,0].astype(int)-1]
- md2.mesh.edges[: ,1]=Pnode[md2.mesh.edges[:,1].astype(int)-1]
- md2.mesh.edges[: ,2]=Pelem[md2.mesh.edges[:,2].astype(int)-1]
- md2.mesh.edges[pos,3]=Pelem[md2.mesh.edges[pos,3].astype(int)-1]
+ md2.mesh.edges[: ,0]=Pnode[md2.mesh.edges[:,0]-1]
+ md2.mesh.edges[: ,1]=Pnode[md2.mesh.edges[:,1]-1]
+ md2.mesh.edges[: ,2]=Pelem[md2.mesh.edges[:,2]-1]
+ md2.mesh.edges[pos,3]=Pelem[md2.mesh.edges[pos,3]-1]
#remove edges when the 2 vertices are not in the domain.
md2.mesh.edges=md2.mesh.edges[numpy.nonzero(numpy.logical_and(md2.mesh.edges[:,0],md2.mesh.edges[:,1]))[0],:]
#Replace all zeros by -1 in the last two columns
@@ -363,15 +364,15 @@
[md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices)
[md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity)
md2.mesh.segments=contourenvelope(md2)
- md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2)
- md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2].astype(int)-1]=1
+ md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2,bool)
+ md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2]-1]=True
else:
#First do the connectivity for the contourenvelope in 2d
[md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements2d,md2.mesh.numberofvertices2d)
[md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity)
md2.mesh.segments=contourenvelope(md2)
- md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2/md2.mesh.numberoflayers)
- md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2].astype(int)-1]=1
+ md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2/md2.mesh.numberoflayers,bool)
+ md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2]-1]=True
md2.mesh.vertexonboundary=numpy.tile(md2.mesh.vertexonboundary,md2.mesh.numberoflayers)
#Then do it for 3d as usual
[md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices)
@@ -380,7 +381,7 @@
#Boundary conditions: Dirichlets on new boundary
#Catch the elements that have not been extracted
orphans_elem=numpy.nonzero(numpy.logical_not(flag_elem))[0]
- orphans_node=numpy.unique(md1.mesh.elements[orphans_elem,:]).astype(int)-1
+ orphans_node=numpy.unique(md1.mesh.elements[orphans_elem,:])-1
#Figure out which node are on the boundary between md2 and md1
nodestoflag1=numpy.intersect1d(orphans_node,pos_node)
nodestoflag2=Pnode[nodestoflag1].astype(int)-1
@@ -442,8 +443,8 @@
setattr(fieldr,solutionsubfield,subfield)
#Keep track of pos_node and pos_elem
- md2.mesh.extractedvertices=pos_node.astype(float)+1
- md2.mesh.extractedelements=pos_elem.astype(float)+1
+ md2.mesh.extractedvertices=pos_node+1
+ md2.mesh.extractedelements=pos_elem+1
return md2
# }}}
@@ -525,22 +526,22 @@
number_nodes3d=numpy.size(x3d) #number of 3d nodes for the non extruded part of the mesh
#Extrude elements
- elements3d=numpy.empty((0,6))
+ elements3d=numpy.empty((0,6),int)
for i in xrange(numlayers-1):
elements3d=numpy.vstack((elements3d,numpy.hstack((md.mesh.elements+i*md.mesh.numberofvertices,md.mesh.elements+(i+1)*md.mesh.numberofvertices)))) #Create the elements of the 3d mesh for the non extruded part
number_el3d=numpy.size(elements3d,axis=0) #number of 3d nodes for the non extruded part of the mesh
#Keep a trace of lower and upper nodes
- mesh.lowervertex=float('NaN')*numpy.ones(number_nodes3d)
- mesh.uppervertex=float('NaN')*numpy.ones(number_nodes3d)
+ mesh.lowervertex=-1*numpy.ones(number_nodes3d,int)
+ mesh.uppervertex=-1*numpy.ones(number_nodes3d,int)
mesh.lowervertex[md.mesh.numberofvertices:]=numpy.arange(1,(numlayers-1)*md.mesh.numberofvertices+1)
mesh.uppervertex[:(numlayers-1)*md.mesh.numberofvertices]=numpy.arange(md.mesh.numberofvertices+1,number_nodes3d+1)
md.mesh.lowervertex=mesh.lowervertex
md.mesh.uppervertex=mesh.uppervertex
#same for lower and upper elements
- mesh.lowerelements=float('NaN')*numpy.ones(number_el3d)
- mesh.upperelements=float('NaN')*numpy.ones(number_el3d)
+ mesh.lowerelements=-1*numpy.ones(number_el3d,int)
+ mesh.upperelements=-1*numpy.ones(number_el3d,int)
mesh.lowerelements[md.mesh.numberofelements:]=numpy.arange(1,(numlayers-2)*md.mesh.numberofelements+1)
mesh.upperelements[:(numlayers-2)*md.mesh.numberofelements]=numpy.arange(md.mesh.numberofelements+1,(numlayers-1)*md.mesh.numberofelements+1)
md.mesh.lowerelements=mesh.lowerelements
@@ -602,21 +603,21 @@
md.initialization.waterfraction=project3d(md,'vector',md.initialization.waterfraction,'type','node')
#bedinfo and surface info
- md.mesh.elementonbed=project3d(md,'vector',numpy.ones(md.mesh.numberofelements2d),'type','element','layer',1)
- md.mesh.elementonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofelements2d),'type','element','layer',md.mesh.numberoflayers-1)
- md.mesh.vertexonbed=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d),'type','node','layer',1)
- md.mesh.vertexonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d),'type','node','layer',md.mesh.numberoflayers)
+ md.mesh.elementonbed=project3d(md,'vector',numpy.ones(md.mesh.numberofelements2d,bool),'type','element','layer',1)
+ md.mesh.elementonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofelements2d,bool),'type','element','layer',md.mesh.numberoflayers-1)
+ md.mesh.vertexonbed=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',1)
+ md.mesh.vertexonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',md.mesh.numberoflayers)
#elementstype
if not numpy.any(numpy.isnan(md.flowequation.element_equation)):
oldelements_type=md.flowequation.element_equation
- md.flowequation.element_equation=numpy.zeros(number_el3d)
+ md.flowequation.element_equation=numpy.zeros(number_el3d,int)
md.flowequation.element_equation=project3d(md,'vector',oldelements_type,'type','element')
#verticestype
if not numpy.any(numpy.isnan(md.flowequation.vertex_equation)):
oldvertices_type=md.flowequation.vertex_equation
- md.flowequation.vertex_equation=numpy.zeros(number_nodes3d)
+ md.flowequation.vertex_equation=numpy.zeros(number_nodes3d,int)
md.flowequation.vertex_equation=project3d(md,'vector',oldvertices_type,'type','node')
md.flowequation.bordermacayeal=project3d(md,'vector',md.flowequation.bordermacayeal,'type','node')
@@ -634,18 +635,18 @@
#in 3d, pressureload: [node1 node2 node3 node4 element]
pressureload_layer1=numpy.hstack((md.diagnostic.icefront[:,0:2],md.diagnostic.icefront[:,1:2]+md.mesh.numberofvertices2d,md.diagnostic.icefront[:,0:1]+md.mesh.numberofvertices2d,md.diagnostic.icefront[:,2:4])) #Add two columns on the first layer
- pressureload=numpy.empty((0,6))
+ pressureload=numpy.empty((0,6),int)
for i in xrange(numlayers-1):
pressureload=numpy.vstack((pressureload,numpy.hstack((pressureload_layer1[:,0:4]+i*md.mesh.numberofvertices2d,pressureload_layer1[:,4:5]+i*md.mesh.numberofelements2d,pressureload_layer1[:,5:6]))))
md.diagnostic.icefront=pressureload
#connectivity
md.mesh.elementconnectivity=numpy.tile(md.mesh.elementconnectivity,(numlayers-1,1))
- md.mesh.elementconnectivity[numpy.nonzero(md.mesh.elementconnectivity==0)]=float('NaN')
+ md.mesh.elementconnectivity[numpy.nonzero(md.mesh.elementconnectivity==0)]=-sys.maxint-1
for i in xrange(1,numlayers-1):
md.mesh.elementconnectivity[i*md.mesh.numberofelements2d:(i+1)*md.mesh.numberofelements2d,:] \
=md.mesh.elementconnectivity[i*md.mesh.numberofelements2d:(i+1)*md.mesh.numberofelements2d,:]+md.mesh.numberofelements2d
- md.mesh.elementconnectivity[numpy.nonzero(numpy.isnan(md.mesh.elementconnectivity))]=0
+ md.mesh.elementconnectivity[numpy.nonzero(md.mesh.elementconnectivity<0)]=0
#materials
md.materials.rheology_B=project3d(md,'vector',md.materials.rheology_B,'type','node')
Modified: issm/trunk/src/m/classes/transient.m
===================================================================
--- issm/trunk/src/m/classes/transient.m 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/classes/transient.m 2012-11-30 18:09:49 UTC (rev 14067)
@@ -45,8 +45,8 @@
disp(sprintf(' transient solution parameters:'));
fielddisplay(obj,'isprognostic','indicates if a prognostic solution is used in the transient');
+ fielddisplay(obj,'isdiagnostic','indicates if a diagnostic solution is used in the transient');
fielddisplay(obj,'isthermal','indicates if a thermal solution is used in the transient');
- fielddisplay(obj,'isdiagnostic','indicates if a diagnostic solution is used in the transient');
fielddisplay(obj,'isgroundingline','indicates if a groundingline migration is used in the transient');
fielddisplay(obj,'requested_outputs','list of additional outputs requested');
Modified: issm/trunk/src/m/classes/transient.py
===================================================================
--- issm/trunk/src/m/classes/transient.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/classes/transient.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -15,10 +15,10 @@
#properties
def __init__(self):
# {{{ Properties
- self.isprognostic = 0
- self.isdiagnostic = 0
- self.isthermal = 0
- self.isgroundingline = 0
+ self.isprognostic = False
+ self.isdiagnostic = False
+ self.isthermal = False
+ self.isgroundingline = False
self.requested_outputs = float('NaN')
#set defaults
@@ -29,8 +29,8 @@
# {{{ Display
string=' transient solution parameters:'
string="%s\n%s"%(string,fielddisplay(self,'isprognostic','indicates if a prognostic solution is used in the transient'))
+ string="%s\n%s"%(string,fielddisplay(self,'isdiagnostic','indicates if a diagnostic solution is used in the transient'))
string="%s\n%s"%(string,fielddisplay(self,'isthermal','indicates if a thermal solution is used in the transient'))
- string="%s\n%s"%(string,fielddisplay(self,'isdiagnostic','indicates if a diagnostic solution is used in the transient'))
string="%s\n%s"%(string,fielddisplay(self,'isgroundingline','indicates if a groundingline migration is used in the transient'))
string="%s\n%s"%(string,fielddisplay(self,'requested_outputs','list of additional outputs requested'))
return string
@@ -40,10 +40,10 @@
# {{{setdefaultparameters
#full analysis: Diagnostic, Prognostic and Thermal but no groundingline migration for now
- self.isprognostic=1
- self.isdiagnostic=1
- self.isthermal=1
- self.isgroundingline=0
+ self.isprognostic=True
+ self.isdiagnostic=True
+ self.isthermal=True
+ self.isgroundingline=False
return self
#}}}
Modified: issm/trunk/src/m/consistency/checkfield.py
===================================================================
--- issm/trunk/src/m/consistency/checkfield.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/consistency/checkfield.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -142,7 +142,7 @@
#Check forcings (size and times)
if options.getfieldvalue('forcing',0):
if numpy.size(field,0)==md.mesh.numberofvertices:
- if len(numpy.shape(field))>1 and not numpy.size(field,1)==1:
+ if numpy.ndim(field)>1 and not numpy.size(field,1)==1:
md = md.checkmessage(options.getfieldvalue('message',\
"field '%s' should have only one column as there are md.mesh.numberofvertices lines" % fieldname))
elif numpy.size(field,0)==md.mesh.numberofvertices+1:
Modified: issm/trunk/src/m/exp/expcoarsen.m
===================================================================
--- issm/trunk/src/m/exp/expcoarsen.m 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/exp/expcoarsen.m 2012-11-30 18:09:49 UTC (rev 14067)
@@ -24,7 +24,7 @@
end
%Get exp oldfile
-[path root ext ver]=fileparts(oldfile);
+[path root ext]=fileparts(oldfile);
A=expread(oldfile);
numprofiles=size(A,2);
Modified: issm/trunk/src/m/extrusion/project3d.py
===================================================================
--- issm/trunk/src/m/extrusion/project3d.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/extrusion/project3d.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -37,7 +37,7 @@
paddingvalue = options.getfieldvalue('padding',0) #0 by default
vector1d=False
- if isinstance(vector2d,numpy.ndarray) and len(vector2d.shape)==1:
+ if isinstance(vector2d,numpy.ndarray) and numpy.ndim(vector2d)==1:
vector1d=True
vector2d=vector2d.reshape(numpy.size(vector2d),1)
@@ -48,9 +48,9 @@
#Initialize 3d vector
if numpy.size(vector2d,axis=0)==md.mesh.numberofvertices2d:
- projected_vector=paddingvalue*numpy.ones((md.mesh.numberofvertices, numpy.size(vector2d,axis=1)))
+ projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofvertices, numpy.size(vector2d,axis=1)))).astype(vector2d.dtype)
elif numpy.size(vector2d,axis=0)==md.mesh.numberofvertices2d+1:
- projected_vector=paddingvalue*numpy.ones((md.mesh.numberofvertices+1,numpy.size(vector2d,axis=1)))
+ projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofvertices+1,numpy.size(vector2d,axis=1)))).astype(vector2d.dtype)
projected_vector[-1,:]=vector2d[-1,:]
vector2d=vector2d[:-1,:]
else:
@@ -67,9 +67,9 @@
#Initialize 3d vector
if numpy.size(vector2d,axis=0)==md.mesh.numberofelements2d:
- projected_vector=paddingvalue*numpy.ones((md.mesh.numberofelements, numpy.size(vector2d,axis=1)))
+ projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofelements, numpy.size(vector2d,axis=1)))).astype(vector2d.dtype)
elif numpy.size(vector2d,axis=0)==md.mesh.numberofelements2d+1:
- projected_vector=paddingvalue*numpy.ones((md.mesh.numberofelements+1,numpy.size(vector2d,axis=1)))
+ projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofelements+1,numpy.size(vector2d,axis=1)))).astype(vector2d.dtype)
projected_vector[-1,:]=vector2d[-1,:]
vector2d=vector2d[:-1,:]
else:
Modified: issm/trunk/src/m/geometry/FlagElements.py
===================================================================
--- issm/trunk/src/m/geometry/FlagElements.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/geometry/FlagElements.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -23,10 +23,10 @@
if isinstance(region,(str,unicode)):
if not region:
- flag=numpy.zeros(md.mesh.numberofelements,'bool')
+ flag=numpy.zeros(md.mesh.numberofelements,bool)
invert=0
elif strcmpi(region,'all'):
- flag=numpy.ones(md.mesh.numberofelements,'bool')
+ flag=numpy.ones(md.mesh.numberofelements,bool)
invert=0
else:
#make sure that we actually don't want the elements outside the domain outline!
@@ -42,11 +42,12 @@
raise IOError("Error: File 'region' not found!" % region)
raise RuntimeError("FlagElements.py calling basinzoom.py is not complete.")
xlim,ylim=basinzoom('basin',region)
- flag_nodes=numpy.logical_and(numpy.logical_and(md.mesh.x<xlim[1],md.mesh.x>xlim[0]),numpy.logical_and(md.mesh.y<ylim[1],md.mesh.y>ylim[0])).astype(float)
- flag=numpy.prod(flag_nodes[md.mesh.elements],axis=1)
+ flag_nodes=numpy.logical_and(numpy.logical_and(md.mesh.x<xlim[1],md.mesh.x>xlim[0]),numpy.logical_and(md.mesh.y<ylim[1],md.mesh.y>ylim[0]))
+ flag=numpy.prod(flag_nodes[md.mesh.elements],axis=1).astype(bool)
else:
#ok, flag elements
[flag,dum]=ContourToMesh(md.mesh.elements[:,0:3].copy(),md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),region,'element',1)
+ flag=flag.astype(bool)
if invert:
flag=numpy.logical_not(flag)
Modified: issm/trunk/src/m/geometry/GetAreas.py
===================================================================
--- issm/trunk/src/m/geometry/GetAreas.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/geometry/GetAreas.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -32,12 +32,12 @@
#initialization
areas=numpy.zeros(nels)
- x1=x[index[:,0].astype(int)-1]
- x2=x[index[:,1].astype(int)-1]
- x3=x[index[:,2].astype(int)-1]
- y1=y[index[:,0].astype(int)-1]
- y2=y[index[:,1].astype(int)-1]
- y3=y[index[:,2].astype(int)-1]
+ x1=x[index[:,0]-1]
+ x2=x[index[:,1]-1]
+ x3=x[index[:,2]-1]
+ y1=y[index[:,0]-1]
+ y2=y[index[:,1]-1]
+ y3=y[index[:,2]-1]
#compute the volume of each element
if not z:
@@ -45,7 +45,7 @@
areas=(0.5*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)))
else:
#V=area(triangle)*1/3(z1+z2+z3)
- thickness=numpy.mean(z[index[:,3:6].astype(int)-1])-numpy.mean(z[index[:,0:3].astype(int)-1])
+ thickness=numpy.mean(z[index[:,3:6]-1])-numpy.mean(z[index[:,0:3]-1])
areas=(0.5*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)))*thickness
return areas
Modified: issm/trunk/src/m/mesh/ComputeHessian.py
===================================================================
--- issm/trunk/src/m/mesh/ComputeHessian.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/mesh/ComputeHessian.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -44,8 +44,8 @@
field=sparse(line,numpy.ones((linesize,1)),numpy.tile(areas*field,(3,1)),numberofnodes,1)/weights
#Compute gradient for each element
- grad_elx=numpy.sum(field[index.astype(int)-1,0]*alpha,axis=1)
- grad_ely=numpy.sum(field[index.astype(int)-1,0]*beta,axis=1)
+ grad_elx=numpy.sum(field[index-1,0]*alpha,axis=1)
+ grad_ely=numpy.sum(field[index-1,0]*beta,axis=1)
#Compute gradient for each node (average of the elements around)
gradx=sparse(line,numpy.ones((linesize,1)),numpy.tile((areas*grad_elx).reshape(-1,1),(3,1)),numberofnodes,1)
@@ -54,7 +54,7 @@
grady=grady/weights
#Compute hessian for each element
- hessian=numpy.hstack((numpy.sum(gradx[index.astype(int)-1,0]*alpha,axis=1).reshape(-1,1),numpy.sum(grady[index.astype(int)-1,0]*alpha,axis=1).reshape(-1,1),numpy.sum(grady[index.astype(int)-1,0]*beta,axis=1).reshape(-1,1)))
+ hessian=numpy.hstack((numpy.sum(gradx[index-1,0]*alpha,axis=1).reshape(-1,1),numpy.sum(grady[index-1,0]*alpha,axis=1).reshape(-1,1),numpy.sum(grady[index-1,0]*beta,axis=1).reshape(-1,1)))
if strcmpi(type,'node'):
#Compute Hessian on the nodes (average of the elements around)
Modified: issm/trunk/src/m/mesh/GetNodalFunctionsCoeff.py
===================================================================
--- issm/trunk/src/m/mesh/GetNodalFunctionsCoeff.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/mesh/GetNodalFunctionsCoeff.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -38,12 +38,12 @@
beta=numpy.zeros((nels,3))
#compute nodal functions coefficients N(x,y)=alpha x + beta y +gamma
- x1=x[index[:,0].astype(int)-1]
- x2=x[index[:,1].astype(int)-1]
- x3=x[index[:,2].astype(int)-1]
- y1=y[index[:,0].astype(int)-1]
- y2=y[index[:,1].astype(int)-1]
- y3=y[index[:,2].astype(int)-1]
+ x1=x[index[:,0]-1]
+ x2=x[index[:,1]-1]
+ x3=x[index[:,2]-1]
+ y1=y[index[:,0]-1]
+ y2=y[index[:,1]-1]
+ y3=y[index[:,2]-1]
invdet=1./(x1*(y2-y3)-x2*(y1-y3)+x3*(y1-y2))
#get alpha and beta
Modified: issm/trunk/src/m/mesh/bamg.py
===================================================================
--- issm/trunk/src/m/mesh/bamg.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/mesh/bamg.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -319,10 +319,10 @@
md.private.bamg['geometry']=bamggeom(bamggeom_out)
md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
md.mesh.y=bamgmesh_out['Vertices'][:,1].copy()
- md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].copy()
- md.mesh.edges=bamgmesh_out['IssmEdges'].copy()
- md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].copy()
- md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].copy()
+ md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int)
+ md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int)
+ md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int)
+ md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int)
#Fill in rest of fields:
md.mesh.dimension=2
@@ -330,15 +330,16 @@
md.mesh.numberofvertices=numpy.size(md.mesh.x)
md.mesh.numberofedges=numpy.size(md.mesh.edges,axis=0)
md.mesh.z=numpy.zeros(md.mesh.numberofvertices)
- md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices)
- md.mask.vertexonwater=numpy.zeros(md.mesh.numberofvertices)
- md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices)
- md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements)
- md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements)
- md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices)
- md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1]=1
+ md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices,bool)
+ md.mask.vertexonwater=numpy.zeros(md.mesh.numberofvertices,bool)
+ md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices,bool)
+ md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements,bool)
+ md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements,bool)
+ md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices,bool)
+ md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
md.mesh.elementconnectivity=md.private.bamg['mesh'].ElementConnectivity
md.mesh.elementconnectivity[numpy.nonzero(numpy.isnan(md.mesh.elementconnectivity))]=0
+ md.mesh.elementconnectivity=md.mesh.elementconnectivity.astype(int)
#Check for orphan
if numpy.any(numpy.logical_not(numpy.in1d(numpy.arange(1,md.mesh.numberofvertices+1),md.mesh.elements.flat))):
Modified: issm/trunk/src/m/mesh/meshconvert.py
===================================================================
--- issm/trunk/src/m/mesh/meshconvert.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/mesh/meshconvert.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -34,10 +34,10 @@
md.private.bamg['geometry']=bamggeom(bamggeom_out)
md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
md.mesh.y=bamgmesh_out['Vertices'][:,1].copy()
- md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].copy()
- md.mesh.edges=bamgmesh_out['IssmEdges'].copy()
- md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].copy()
- md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].copy()
+ md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int)
+ md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int)
+ md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int)
+ md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int)
#Fill in rest of fields:
md.mesh.dimension=2
@@ -45,13 +45,13 @@
md.mesh.numberofvertices=numpy.size(md.mesh.x)
md.mesh.numberofedges=numpy.size(md.mesh.edges,axis=0)
md.mesh.z=numpy.zeros(md.mesh.numberofvertices)
- md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices)
- md.mask.vertexonwater=numpy.zeros(md.mesh.numberofvertices)
- md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices)
- md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements)
- md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements)
- md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices)
- md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1]=1
+ md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices,bool)
+ md.mask.vertexonwater=numpy.zeros(md.mesh.numberofvertices,bool)
+ md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices,bool)
+ md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements,bool)
+ md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements,bool)
+ md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices,bool)
+ md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
return md
Modified: issm/trunk/src/m/mesh/rifts/meshprocessoutsiderifts.py
===================================================================
--- issm/trunk/src/m/mesh/rifts/meshprocessoutsiderifts.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/mesh/rifts/meshprocessoutsiderifts.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -40,7 +40,7 @@
A=tip
B=node_connected_to_tip
- elements=numpy.empty(0)
+ elements=numpy.empty(0,int)
while flags(B): #as long as B does not belong to the domain outline, keep looking.
#detect elements on edge A,B:
@@ -61,10 +61,10 @@
md.mesh.numberofvertices=num
#replace tip in elements
- newelements=md.mesh.elements[elements.astype(int)-1,:]
+ newelements=md.mesh.elements[elements-1,:]
pos=numpy.nonzero(newelements==tip)
newelements[pos]=num
- md.mesh.elements[elements.astype(int)-1,:]=newelements
+ md.mesh.elements[elements-1,:]=newelements
rift.tips=numpy.concatenate((rift.tips,num))
#deal with segments
@@ -80,14 +80,14 @@
md.mesh.numberofelements=numpy.size(md.mesh.elements,axis=0)
md.mesh.numberofvertices=numpy.size(md.mesh.x)
md.mesh.z=numpy.zeros(md.mesh.numberofvertices)
- md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x))
- md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1]=1
+ md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x),bool)
+ md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
md.rifts.numrifts=length(md.rifts.riftstruct)
- md.flowequation.element_equation=3.*numpy.ones(md.mesh.numberofelements)
- md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices)
- md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices)
- md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements)
- md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements)
+ md.flowequation.element_equation=3*numpy.ones(md.mesh.numberofelements,int)
+ md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices,bool)
+ md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices,bool)
+ md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements,bool)
+ md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements,bool)
return md
Modified: issm/trunk/src/m/mesh/rifts/meshprocessrifts.py
===================================================================
--- issm/trunk/src/m/mesh/rifts/meshprocessrifts.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/mesh/rifts/meshprocessrifts.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -22,8 +22,11 @@
#Call MEX file
[md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers,md.rifts.riftstruct]=TriMeshProcessRifts(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers)
+ md.mesh.elements=md.mesh.elements.astype(int)
md.mesh.x=md.mesh.x.reshape(-1)
md.mesh.y=md.mesh.y.reshape(-1)
+ md.mesh.segments=md.mesh.segments.astype(int)
+ md.mesh.segmentmarkers=md.mesh.segmentmarkers.astype(int)
if not isinstance(md.rifts.riftstruct,list) or not md.rifts.riftstruct:
raise RuntimeError("TriMeshProcessRifts did not find any rift")
@@ -32,12 +35,12 @@
md.mesh.numberofelements=numpy.size(md.mesh.elements,axis=0)
md.mesh.numberofvertices=numpy.size(md.mesh.x)
md.mesh.z=numpy.zeros(md.mesh.numberofvertices)
- md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x))
- md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1]=1
- md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices)
- md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices)
- md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements)
- md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements)
+ md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x),bool)
+ md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
+ md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices,bool)
+ md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices,bool)
+ md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements,bool)
+ md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements,bool)
#get coordinates of rift tips
for rift in md.rifts.riftstruct:
Modified: issm/trunk/src/m/mesh/triangle.py
===================================================================
--- issm/trunk/src/m/mesh/triangle.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/mesh/triangle.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -43,20 +43,23 @@
#Mesh using TriMesh
[md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers]=TriMesh(domainname,riftname,area)
+ md.mesh.elements=md.mesh.elements.astype(int)
+ md.mesh.segments=md.mesh.segments.astype(int)
+ md.mesh.segmentmarkers=md.mesh.segmentmarkers.astype(int)
#Fill in rest of fields:
md.mesh.numberofelements = numpy.size(md.mesh.elements,axis=0)
md.mesh.numberofvertices = numpy.size(md.mesh.x)
md.mesh.z = numpy.zeros(md.mesh.numberofvertices)
- md.mesh.vertexonboundary = numpy.zeros(md.mesh.numberofvertices)
- md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1] = 1.
- md.mesh.vertexonbed = numpy.ones(md.mesh.numberofvertices)
- md.mesh.vertexonsurface = numpy.ones(md.mesh.numberofvertices)
- md.mesh.elementonbed = numpy.ones(md.mesh.numberofelements)
- md.mesh.elementonsurface = numpy.ones(md.mesh.numberofelements)
+ md.mesh.vertexonboundary = numpy.zeros(md.mesh.numberofvertices,bool)
+ md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1] = True
+ md.mesh.vertexonbed = numpy.ones(md.mesh.numberofvertices,bool)
+ md.mesh.vertexonsurface = numpy.ones(md.mesh.numberofvertices,bool)
+ md.mesh.elementonbed = numpy.ones(md.mesh.numberofelements,bool)
+ md.mesh.elementonsurface = numpy.ones(md.mesh.numberofelements,bool)
#Now, build the connectivity tables for this mesh.
- [md.mesh.vertexconnectivity]= NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)
+ [md.mesh.vertexconnectivity] = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)
[md.mesh.elementconnectivity] = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)
#type of model
Modified: issm/trunk/src/m/parameterization/contourenvelope.py
===================================================================
--- issm/trunk/src/m/parameterization/contourenvelope.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/parameterization/contourenvelope.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -75,7 +75,7 @@
pos=numpy.nonzero(flags)
elemin[pos]=1
- nodein[mesh.elements[pos,:].astype(int)-1]=1
+ nodein[mesh.elements[pos,:]-1]=1
#modify element connectivity
elemout=numpy.nonzero(numpy.logical_not(elemin))[0]
@@ -87,16 +87,16 @@
flag=copy.deepcopy(mesh.elementconnectivity)
if len(args)==1:
flag[numpy.nonzero(flag)]=elemin[flag[numpy.nonzero(flag)]]
- elementonboundary=numpy.logical_and(numpy.prod(flag,axis=1)==0,numpy.sum(flag,axis=1)>0).astype(float)
+ elementonboundary=numpy.logical_and(numpy.prod(flag,axis=1)==0,numpy.sum(flag,axis=1)>0)
#Find segments on boundary
pos=numpy.nonzero(elementonboundary)[0]
num_segments=numpy.size(pos)
- segments=numpy.zeros((num_segments*3,3))
+ segments=numpy.zeros((num_segments*3,3),int)
count=0
for el1 in pos:
- els2=mesh.elementconnectivity[el1,numpy.nonzero(mesh.elementconnectivity[el1,:])[0]].astype(int)-1
+ els2=mesh.elementconnectivity[el1,numpy.nonzero(mesh.elementconnectivity[el1,:])[0]]-1
if numpy.size(els2)>1:
flag=numpy.intersect1d(mesh.elements[els2[0],:],mesh.elements[els2[1],:])
nods1=mesh.elements[el1,:]
Modified: issm/trunk/src/m/parameterization/setflowequation.py
===================================================================
--- issm/trunk/src/m/parameterization/setflowequation.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/parameterization/setflowequation.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -49,11 +49,11 @@
#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
+ hutterflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(macayealflag,pattynflag)))]=True
elif strcmpi(filltype,'macayeal'):
- macayealflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(pattynflag,stokesflag))))]=1
+ macayealflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(pattynflag,stokesflag))))]=True
elif strcmpi(filltype,'pattyn'):
- pattynflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(macayealflag,stokesflag))))]=1
+ pattynflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(macayealflag,stokesflag))))]=True
#check that each element has at least one flag
if not any(hutterflag+macayealflag+l1l2flag+pattynflag+stokesflag):
@@ -62,9 +62,9 @@
#check that each element has only one flag
if any(hutterflag+macayealflag+l1l2flag+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
+ hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,macayealflag))]=False
+ hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,pattynflag))]=False
+ macayealflag[numpy.nonzero(numpy.logical_and(macayealflag,pattynflag))]=False
#Check that no pattyn or stokes for 2d mesh
if md.mesh.dimension==2:
@@ -76,16 +76,16 @@
raise TypeError("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
- nodeonl1l2=numpy.zeros(md.mesh.numberofvertices)
- nodeonl1l2[md.mesh.elements[numpy.nonzero(l1l2flag),:].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)
+ nodeonhutter=numpy.zeros(md.mesh.numberofvertices,bool)
+ nodeonhutter[md.mesh.elements[numpy.nonzero(hutterflag),:]-1]=True
+ nodeonmacayeal=numpy.zeros(md.mesh.numberofvertices,bool)
+ nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True
+ nodeonl1l2=numpy.zeros(md.mesh.numberofvertices,bool)
+ nodeonl1l2[md.mesh.elements[numpy.nonzero(l1l2flag),:]-1]=True
+ nodeonpattyn=numpy.zeros(md.mesh.numberofvertices,bool)
+ nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True
+ nodeonstokes=numpy.zeros(md.mesh.numberofvertices,bool)
+ noneflag=numpy.zeros(md.mesh.numberofelements,bool)
#First modify stokesflag to get rid of elements contrained everywhere (spc + border with pattyn or macayeal)
if any(stokesflag):
@@ -95,29 +95,29 @@
numpy.logical_not(numpy.isnan(md.diagnostic.spcvz)).astype(int)==3, \
numpy.logical_and(nodeonpattyn,nodeonstokes).reshape(-1,1)).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.reshape(-1))]=0
- nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1
+ fullspcelems=(numpy.sum(fullspcnodes[md.mesh.elements-1],axis=1)==6).astype(int) #find all the nodes on the boundary of the domain without icefront
+ stokesflag[numpy.nonzero(fullspcelems.reshape(-1))]=False
+ nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True
#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
+ pattynflag[numpy.logical_not(stokesflag)]=True
+ nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True
elif any(macayealflag): #fill with macayeal
- macayealflag[numpy.logical_not(stokesflag)]=1
- nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1
+ macayealflag[numpy.logical_not(stokesflag)]=True
+ nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True
else: #fill with none
- noneflag[numpy.nonzero(numpy.logical_not(stokesflag))]=1
+ noneflag[numpy.nonzero(numpy.logical_not(stokesflag))]=True
#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)
+ nodeonmacayealpattyn=numpy.zeros(md.mesh.numberofvertices,bool)
+ nodeonpattynstokes=numpy.zeros(md.mesh.numberofvertices,bool)
+ nodeonmacayealstokes=numpy.zeros(md.mesh.numberofvertices,bool)
+ macayealpattynflag=numpy.zeros(md.mesh.numberofelements,bool)
+ macayealstokesflag=numpy.zeros(md.mesh.numberofelements,bool)
+ pattynstokesflag=numpy.zeros(md.mesh.numberofelements,bool)
if strcmpi(coupling_method,'penalties'):
#Create the border nodes between Pattyn and MacAyeal and extrude them
numnodes2d=md.mesh.numberofvertices2d
@@ -134,105 +134,105 @@
elif strcmpi(coupling_method,'tiling'):
if any(macayealflag) and any(pattynflag): #coupling macayeal pattyn
#Find node at the border
- nodeonmacayealpattyn[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonpattyn))]=1
+ nodeonmacayealpattyn[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonpattyn))]=True
#Macayeal elements in contact with this layer become MacAyealPattyn elements
- matrixelements=ismember(md.mesh.elements.astype(int)-1,numpy.nonzero(nodeonmacayealpattyn)[0])
+ matrixelements=ismember(md.mesh.elements-1,numpy.nonzero(nodeonmacayealpattyn)[0])
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
+ commonelements[numpy.nonzero(pattynflag)]=False #only one layer: the elements previously in macayeal
+ macayealflag[numpy.nonzero(commonelements)]=False #these elements are now macayealpattynelements
+ macayealpattynflag[numpy.nonzero(commonelements)]=True
+ nodeonmacayeal[:]=False
+ nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True
#rule out elements that don't touch the 2 boundaries
pos=numpy.nonzero(macayealpattynflag)[0]
elist=numpy.zeros(numpy.size(pos),dtype=int)
- elist = elist + numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:].astype(int)-1],axis=1).astype(bool)
- elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:].astype(int)-1] ,axis=1).astype(bool)
+ elist = elist + numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:]-1],axis=1).astype(bool)
+ elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:]-1] ,axis=1).astype(bool)
pos1=numpy.nonzero(elist==1)[0]
- macayealflag[pos[pos1]]=1
- macayealpattynflag[pos[pos1]]=0
+ macayealflag[pos[pos1]]=True
+ macayealpattynflag[pos[pos1]]=False
pos2=numpy.nonzero(elist==-1)[0]
- pattynflag[pos[pos2]]=1
- macayealpattynflag[pos[pos2]]=0
+ pattynflag[pos[pos2]]=True
+ macayealpattynflag[pos[pos2]]=False
#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
+ nodeonmacayeal[:]=False
+ nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True
+ nodeonpattyn[:]=False
+ nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True
+ nodeonmacayealpattyn[:]=False
+ nodeonmacayealpattyn[md.mesh.elements[numpy.nonzero(macayealpattynflag),:]-1]=True
elif any(pattynflag) and any(stokesflag): #coupling pattyn stokes
#Find node at the border
- nodeonpattynstokes[numpy.nonzero(numpy.logical_and(nodeonpattyn,nodeonstokes))]=1
+ nodeonpattynstokes[numpy.nonzero(numpy.logical_and(nodeonpattyn,nodeonstokes))]=True
#Stokes elements in contact with this layer become PattynStokes elements
- matrixelements=ismember(md.mesh.elements.astype(int)-1,numpy.nonzero(nodeonpattynstokes)[0])
+ matrixelements=ismember(md.mesh.elements-1,numpy.nonzero(nodeonpattynstokes)[0])
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
+ commonelements[numpy.nonzero(pattynflag)]=False #only one layer: the elements previously in macayeal
+ stokesflag[numpy.nonzero(commonelements)]=False #these elements are now macayealpattynelements
+ pattynstokesflag[numpy.nonzero(commonelements)]=True
+ nodeonstokes=numpy.zeros(md.mesh.numberofvertices,bool)
+ nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True
#rule out elements that don't touch the 2 boundaries
pos=numpy.nonzero(pattynstokesflag)[0]
elist=numpy.zeros(numpy.size(pos),dtype=int)
- elist = elist + numpy.sum(nodeonstokes[md.mesh.elements[pos,:].astype(int)-1],axis=1).astype(bool)
- elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:].astype(int)-1],axis=1).astype(bool)
+ elist = elist + numpy.sum(nodeonstokes[md.mesh.elements[pos,:]-1],axis=1).astype(bool)
+ elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:]-1],axis=1).astype(bool)
pos1=numpy.nonzero(elist==1)[0]
- stokesflag[pos[pos1]]=1
- pattynstokesflag[pos[pos1]]=0
+ stokesflag[pos[pos1]]=True
+ pattynstokesflag[pos[pos1]]=False
pos2=numpy.nonzero(elist==-1)[0]
- pattynflag[pos[pos2]]=1
- pattynstokesflag[pos[pos2]]=0
+ pattynflag[pos[pos2]]=True
+ pattynstokesflag[pos[pos2]]=False
#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
+ nodeonstokes[:]=False
+ nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True
+ nodeonpattyn[:]=False
+ nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True
+ nodeonpattynstokes[:]=False
+ nodeonpattynstokes[md.mesh.elements[numpy.nonzero(pattynstokesflag),:]-1]=True
elif any(stokesflag) and any(macayealflag):
#Find node at the border
- nodeonmacayealstokes[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonstokes))]=1
+ nodeonmacayealstokes[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonstokes))]=True
#Stokes elements in contact with this layer become MacAyealStokes elements
- matrixelements=ismember(md.mesh.elements.astype(int)-1,numpy.nonzero(nodeonmacayealstokes)[0])
+ matrixelements=ismember(md.mesh.elements-1,numpy.nonzero(nodeonmacayealstokes)[0])
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
+ commonelements[numpy.nonzero(macayealflag)]=False #only one layer: the elements previously in macayeal
+ stokesflag[numpy.nonzero(commonelements)]=False #these elements are now macayealmacayealelements
+ macayealstokesflag[numpy.nonzero(commonelements)]=True
+ nodeonstokes=numpy.zeros(md.mesh.numberofvertices,bool)
+ nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True
#rule out elements that don't touch the 2 boundaries
pos=numpy.nonzero(macayealstokesflag)[0]
elist=numpy.zeros(numpy.size(pos),dtype=int)
- elist = elist + numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:].astype(int)-1],axis=1).astype(bool)
- elist = elist - numpy.sum(nodeonstokes[md.mesh.elements[pos,:].astype(int)-1] ,axis=1).astype(bool)
+ elist = elist + numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:]-1],axis=1).astype(bool)
+ elist = elist - numpy.sum(nodeonstokes[md.mesh.elements[pos,:]-1] ,axis=1).astype(bool)
pos1=numpy.nonzero(elist==1)[0]
- macayealflag[pos[pos1]]=1
- macayealstokesflag[pos[pos1]]=0
+ macayealflag[pos[pos1]]=True
+ macayealstokesflag[pos[pos1]]=False
pos2=numpy.nonzero(elist==-1)[0]
- stokesflag[pos[pos2]]=1
- macayealstokesflag[pos[pos2]]=0
+ stokesflag[pos[pos2]]=True
+ macayealstokesflag[pos[pos2]]=False
#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
+ nodeonmacayeal[:]=False
+ nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True
+ nodeonstokes[:]=False
+ nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True
+ nodeonmacayealstokes[:]=False
+ nodeonmacayealstokes[md.mesh.elements[numpy.nonzero(macayealstokesflag),:]-1]=True
elif any(stokesflag) and any(hutterflag):
raise TypeError("type of coupling not supported yet")
- #Create MacaAyealPattynApproximation where needed
- md.flowequation.element_equation=numpy.zeros(md.mesh.numberofelements)
+ #Create MacAyealPattynApproximation where needed
+ md.flowequation.element_equation=numpy.zeros(md.mesh.numberofelements,int)
md.flowequation.element_equation[numpy.nonzero(noneflag)]=0
md.flowequation.element_equation[numpy.nonzero(hutterflag)]=1
md.flowequation.element_equation[numpy.nonzero(macayealflag)]=2
@@ -249,7 +249,7 @@
md.flowequation.borderstokes=nodeonstokes
#Create vertices_type
- md.flowequation.vertex_equation=numpy.zeros(md.mesh.numberofvertices)
+ md.flowequation.vertex_equation=numpy.zeros(md.mesh.numberofvertices,int)
pos=numpy.nonzero(nodeonmacayeal)
md.flowequation.vertex_equation[pos]=2
pos=numpy.nonzero(nodeonl1l2)
@@ -272,10 +272,10 @@
md.flowequation.vertex_equation[pos]=6
#figure out solution types
- md.flowequation.ishutter=float(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.isl1l2=float(any(md.flowequation.element_equation==8))
- md.flowequation.isstokes=float(any(md.flowequation.element_equation==4))
+ md.flowequation.ishutter=any(md.flowequation.element_equation==1)
+ md.flowequation.ismacayealpattyn=bool(numpy.any(numpy.logical_or(md.flowequation.element_equation==2,md.flowequation.element_equation==3)))
+ md.flowequation.isl1l2=any(md.flowequation.element_equation==8)
+ md.flowequation.isstokes=any(md.flowequation.element_equation==4)
return md
Modified: issm/trunk/src/m/parameterization/setmask.py
===================================================================
--- issm/trunk/src/m/parameterization/setmask.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/parameterization/setmask.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -36,16 +36,16 @@
#the order here is important. we choose vertexongroundedice as default on the grounding line.
vertexonfloatingice = numpy.zeros(md.mesh.numberofvertices,'bool')
vertexongroundedice = numpy.zeros(md.mesh.numberofvertices,'bool')
- vertexongroundedice[md.mesh.elements[numpy.nonzero(elementongroundedice),:].astype(int)-1]=True
+ vertexongroundedice[md.mesh.elements[numpy.nonzero(elementongroundedice),:]-1]=True
vertexonfloatingice[numpy.nonzero(numpy.logical_not(vertexongroundedice))]=True
#}}}
#Return:
- md.mask.elementonfloatingice = elementonfloatingice.astype(float)
- md.mask.vertexonfloatingice = vertexonfloatingice.astype(float)
- md.mask.elementongroundedice = elementongroundedice.astype(float)
- md.mask.vertexongroundedice = vertexongroundedice.astype(float)
- md.mask.vertexonwater = numpy.zeros(md.mesh.numberofvertices)
- md.mask.elementonwater = numpy.zeros(md.mesh.numberofelements)
+ md.mask.elementonfloatingice = elementonfloatingice
+ md.mask.vertexonfloatingice = vertexonfloatingice
+ md.mask.elementongroundedice = elementongroundedice
+ md.mask.vertexongroundedice = vertexongroundedice
+ md.mask.vertexonwater = numpy.zeros(md.mesh.numberofvertices,bool)
+ md.mask.elementonwater = numpy.zeros(md.mesh.numberofelements,bool)
return md
Modified: issm/trunk/src/m/plot/applyoptions.m
===================================================================
--- issm/trunk/src/m/plot/applyoptions.m 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/plot/applyoptions.m 2012-11-30 18:09:49 UTC (rev 14067)
@@ -204,14 +204,21 @@
tick_vals=get(c,'YTick')';
if exist(options,'log')
logval= getfieldvalue(options,'log');
- for i= 1:size(tick_vals,1)
+ for i= 1:numel(tick_vals)
tick_vals(i)= logval^(tick_vals(i));
end
- elseif size(tick_vals,1) == 3
+ elseif numel(tick_vals) == 3
tick_vals=[tick_vals(1); mean(tick_vals(1:2)); tick_vals(2); ...
mean(tick_vals(2:3)); tick_vals(3)];
set(c,'YTick',tick_vals);
end
+ else
+ if exist(options,'log')
+ logvalue=getfieldvalue(options,'log');
+ set(c,'YTick',log(tick_vals)./log(logvalue));
+ else
+ set(c,'YTick',tick_vals);
+ end
end
labels = cell(1,numel(tick_vals));
for i = 1:numel(tick_vals)
Modified: issm/trunk/src/m/plot/colormaps/getcolormap.m
===================================================================
--- issm/trunk/src/m/plot/colormaps/getcolormap.m 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/plot/colormaps/getcolormap.m 2012-11-30 18:09:49 UTC (rev 14067)
@@ -23,7 +23,7 @@
map = jet(256);
map = map(128:end,:);
elseif strcmpi(map,'redblue'),
- map = hsv;
+ map = hsv(256);
map = rgb2hsv(map);
map(:,2) = max(min( abs(map(:,1)-0.5)/0.5 ,1),0);
map(1:128,1) = 0.7;
Modified: issm/trunk/src/m/plot/radarpower.m
===================================================================
--- issm/trunk/src/m/plot/radarpower.m 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/plot/radarpower.m 2012-11-30 18:09:49 UTC (rev 14067)
@@ -19,6 +19,7 @@
end
end
+geotiff_name = getfieldvalue(options,'geotiff_name','');
highres = getfieldvalue(options,'highres',0);
xlim = getfieldvalue(options,'xlim',[min(md.mesh.x) max(md.mesh.x)]);
ylim = getfieldvalue(options,'ylim',[min(md.mesh.y) max(md.mesh.y)]);
@@ -63,16 +64,18 @@
%md.radaroverlay.pwr=double(flipud(im(1:pixelskip:end,1:pixelskip:end)));
%md.radaroverlay.x=(xmin:(xmax-xmin)/(size(md.radaroverlay.pwr,2)-1):xmax);
%md.radaroverlay.y=(ymin:(ymax-ymin)/(size(md.radaroverlay.pwr,1)-1):ymax);
- if highres,
- if ~exist(['/u/astrid-r1b/ModelData/MOG/mog100_r2_hp1.tif']),
- error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MOG/mog100_r2_hp1.tif not found.']);
+ if ~exist(options,'geotiff_name'),
+ if highres,
+ if ~exist(['/u/astrid-r1b/ModelData/MOG/mog100_r2_hp1.tif']),
+ error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MOG/mog100_r2_hp1.tif not found.']);
+ end
+ geotiff_name=['/u/astrid-r1b/ModelData/MOG/mog100_r2_hp1.tif'];
+ else
+ if ~exist(['/u/astrid-r1b/ModelData/MOG/mog500_r2_hp1.tif']),
+ error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MOG/mog500_r2_hp1.tif not found.']);
+ end
+ geotiff_name=['/u/astrid-r1b/ModelData/MOG/mog500_r2_hp1.tif'];
end
- geotiff_name=['/u/astrid-r1b/ModelData/MOG/mog100_r2_hp1.tif'];
- else
- if ~exist(['/u/astrid-r1b/ModelData/MOG/mog500_r2_hp1.tif']),
- error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MOG/mog500_r2_hp1.tif not found.']);
- end
- geotiff_name=['/u/astrid-r1b/ModelData/MOG/mog500_r2_hp1.tif'];
end
%Name of image
@@ -90,16 +93,18 @@
system('rm -rf ./temp.tif');
elseif strcmpi(md.mesh.hemisphere,'s'),
- if highres,
- if ~exist(['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif']),
- error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif not found.']);
+ if ~exist(options,'geotiff_name'),
+ if highres,
+ if ~exist(['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif']),
+ error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif not found.']);
+ end
+ geotiff_name=['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif'];
+ else
+ if ~exist(['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif']),
+ error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif not found.']);
+ end
+ geotiff_name=['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif'];
end
- geotiff_name=['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif'];
- else
- if ~exist(['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif']),
- error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif not found.']);
- end
- geotiff_name=['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif'];
end
%Name of image
Modified: issm/trunk/src/m/solve/WriteData.py
===================================================================
--- issm/trunk/src/m/solve/WriteData.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/m/solve/WriteData.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -104,7 +104,7 @@
data=numpy.array([data])
elif isinstance(data,(list,tuple)):
data=numpy.array(data).reshape(-1,1)
- if len(data.shape) == 1:
+ if numpy.ndim(data) == 1:
if numpy.size(data):
data=data.reshape(numpy.size(data),1)
else:
@@ -137,7 +137,7 @@
data=numpy.array([data])
elif isinstance(data,(list,tuple)):
data=numpy.array(data).reshape(-1,1)
- if len(data.shape) == 1:
+ if numpy.ndim(data) == 1:
if numpy.size(data):
data=data.reshape(numpy.size(data),1)
else:
@@ -170,7 +170,7 @@
data=numpy.array([data])
elif isinstance(data,(list,tuple)):
data=numpy.array(data).reshape(-1,1)
- if len(data.shape) == 1:
+ if numpy.ndim(data) == 1:
if numpy.size(data):
data=data.reshape(numpy.size(data),1)
else:
@@ -206,7 +206,7 @@
matrix=numpy.array([matrix])
elif isinstance(matrix,(list,tuple)):
matrix=numpy.array(matrix).reshape(-1,1)
- if len(matrix.shape) == 1:
+ if numpy.ndim(matrix) == 1:
if numpy.size(matrix):
matrix=matrix.reshape(numpy.size(matrix),1)
else:
@@ -230,7 +230,7 @@
matrix=numpy.array([matrix])
elif isinstance(matrix,(list,tuple)):
matrix=numpy.array(matrix).reshape(-1,1)
- if len(matrix.shape) == 1:
+ if numpy.ndim(matrix) == 1:
matrix=matrix.reshape(numpy.size(matrix),1)
s=matrix.shape
Modified: issm/trunk/src/wrappers/ElementConnectivity/ElementConnectivity.cpp
===================================================================
--- issm/trunk/src/wrappers/ElementConnectivity/ElementConnectivity.cpp 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/wrappers/ElementConnectivity/ElementConnectivity.cpp 2012-11-30 18:09:49 UTC (rev 14067)
@@ -12,13 +12,13 @@
WRAPPER(ElementConnectivity){
/*inputs: */
- double* elements=NULL;
- double* nodeconnectivity=NULL;
- int nel,nods;
- int width;
+ int* elements=NULL;
+ int* nodeconnectivity=NULL;
+ int nels,nods;
+ int width;
/*outputs: */
- double* elementconnectivity=NULL;
+ int* elementconnectivity=NULL;
/*Boot module: */
MODULEBOOT();
@@ -27,14 +27,14 @@
CHECKARGUMENTS(NLHS,NRHS,&ElementConnectivityUsage);
/*Input datasets: */
- FetchData(&elements,&nel,NULL,ELEMENTS);
+ FetchData(&elements,&nels,NULL,ELEMENTS);
FetchData(&nodeconnectivity,&nods,&width,NODECONNECTIVITY);
/*!Generate internal degree of freedom numbers: */
- ElementConnectivityx(&elementconnectivity, elements,nel, nodeconnectivity, nods, width);
+ ElementConnectivityx(&elementconnectivity,elements,nels,nodeconnectivity,nods,width);
/*write output datasets: */
- WriteData(ELEMENTCONNECTIVITY,elementconnectivity,nel,3);
+ WriteData(ELEMENTCONNECTIVITY,elementconnectivity,nels,3);
/*end module: */
MODULEEND();
Modified: issm/trunk/src/wrappers/Kriging/Kriging.cpp
===================================================================
--- issm/trunk/src/wrappers/Kriging/Kriging.cpp 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/wrappers/Kriging/Kriging.cpp 2012-11-30 18:09:49 UTC (rev 14067)
@@ -4,6 +4,10 @@
#include "./Kriging.h"
void KrigingUsage(void){/*{{{*/
+ int num=1;
+#ifdef _MULTITHREADING_
+ num=_NUMTHREADS_;
+#endif
_pprintLine_("");
_pprintLine_(" usage: predictions=" << __FUNCT__ << "(x,y,observations,x_interp,y_interp,'options');");
_pprintLine_(" available options:");
@@ -20,6 +24,7 @@
_pprintLine_(" -'maxtrimming': maximum trimming value (default is -1.e+21)");
_pprintLine_(" -'mintrimming': minimum trimming value (default is +1.e+21)");
_pprintLine_(" -'minspacing': minimum distance between observation (default is 0.01)");
+ _pprintLine_(" -'numthreads': number of threads, default is "<<num);
_pprintLine_("");
}/*}}}*/
WRAPPER(Kriging){
Modified: issm/trunk/src/wrappers/NodeConnectivity/NodeConnectivity.cpp
===================================================================
--- issm/trunk/src/wrappers/NodeConnectivity/NodeConnectivity.cpp 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/wrappers/NodeConnectivity/NodeConnectivity.cpp 2012-11-30 18:09:49 UTC (rev 14067)
@@ -12,13 +12,13 @@
WRAPPER(NodeConnectivity){
/*inputs: */
- double* elements=NULL;
- int nel;
- int nods;
+ int* elements=NULL;
+ int nels;
+ int nods;
/*outputs: */
- double* connectivity=NULL;
- int width;
+ int* connectivity=NULL;
+ int width;
/*Boot module: */
MODULEBOOT();
@@ -27,11 +27,11 @@
CHECKARGUMENTS(NLHS,NRHS,&NodeConnectivityUsage);
/*Input datasets: */
- FetchData(&elements,&nel,NULL,ELEMENTS);
+ FetchData(&elements,&nels,NULL,ELEMENTS);
FetchData(&nods,NUMNODES);
/*!Generate internal degree of freedom numbers: */
- NodeConnectivityx(&connectivity, &width,elements,nel, nods);
+ NodeConnectivityx(&connectivity,&width,elements,nels,nods);
/*write output datasets: */
WriteData(CONNECTIVITY,connectivity,nods,width);
Modified: issm/trunk/src/wrappers/matlab/io/FetchMatlabData.cpp
===================================================================
--- issm/trunk/src/wrappers/matlab/io/FetchMatlabData.cpp 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/wrappers/matlab/io/FetchMatlabData.cpp 2012-11-30 18:09:49 UTC (rev 14067)
@@ -263,27 +263,19 @@
void FetchData(double** pvector,int* pM,const mxArray* dataref){
double* outvector=NULL;
- int outvector_rows;
+ int M,N;
- if(mxIsEmpty(dataref)){
- /*Nothing to pick up. Just initialize matrix pointer to NULL: */
- outvector_rows=0;
- outvector=NULL;
- }
- else if (mxIsClass(dataref,"double") ){
+ /*Use Fetch matrix*/
+ FetchData(&outvector,&M,&N,dataref) ;
- /*Convert matlab vector to double* vector: */
- MatlabVectorToDoubleVector(&outvector,&outvector_rows,dataref);
-
+ /*Check that it is a vector*/
+ if(M>0){
+ if(N!=1) _error_("input vector of size " << M << "x" << N << " should have only one column");
}
- else{
- /*This is an error: we don't have the correct input!: */
- _error_("Input parameter of class " << mxGetClassName(dataref) << " not supported yet");
- }
/*Assign output pointers:*/
*pvector=outvector;
- if (pM)*pM=outvector_rows;
+ if(pM)*pM=M;
}
/*}}}*/
/*FUNCTION FetchData(int** pvector,int* pM,const mxArray* dataref){{{*/
Modified: issm/trunk/src/wrappers/python/io/FetchPythonData.cpp
===================================================================
--- issm/trunk/src/wrappers/python/io/FetchPythonData.cpp 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/wrappers/python/io/FetchPythonData.cpp 2012-11-30 18:09:49 UTC (rev 14067)
@@ -66,19 +66,47 @@
int ndim;
npy_intp* dims=NULL;
- /*retrive dimensions: */
+ /*intermediary:*/
+ long* lmatrix=NULL;
+ bool* bmatrix=NULL;
+ int i;
+
+ /*retrieve dimensions: */
ndim=PyArray_NDIM((const PyArrayObject*)py_matrix);
if(ndim!=2)_error_("expecting an MxN matrix in input!");
dims=PyArray_DIMS((PyArrayObject*)py_matrix);
M=dims[0]; N=dims[1];
if (M && N) {
- /*retrieve internal value: */
- dmatrix=(double*)PyArray_DATA((PyArrayObject*)py_matrix);
+ if (PyArray_TYPE((PyArrayObject*)py_matrix) == NPY_DOUBLE) {
+ /*retrieve internal value: */
+ dmatrix=(double*)PyArray_DATA((PyArrayObject*)py_matrix);
- /*copy matrix: */
- matrix=xNew<double>(M*N);
- memcpy(matrix,dmatrix,(M*N)*sizeof(double));
+ /*copy matrix: */
+ matrix=xNew<double>(M*N);
+ memcpy(matrix,dmatrix,(M*N)*sizeof(double));
+ }
+
+ else if (PyArray_TYPE((PyArrayObject*)py_matrix) == NPY_INT64) {
+ /*retrieve internal value: */
+ lmatrix=(long*)PyArray_DATA((PyArrayObject*)py_matrix);
+
+ /*transform into double matrix: */
+ matrix=xNew<double>(M*N);
+ for(i=0;i<M*N;i++)matrix[i]=(double)lmatrix[i];
+ }
+
+ else if (PyArray_TYPE((PyArrayObject*)py_matrix) == NPY_BOOL) {
+ /*retrieve internal value: */
+ bmatrix=(bool*)PyArray_DATA((PyArrayObject*)py_matrix);
+
+ /*transform into double matrix: */
+ matrix=xNew<double>(M*N);
+ for(i=0;i<M*N;i++)matrix[i]=(double)bmatrix[i];
+ }
+
+ else
+ _error_("unrecognized matrix type in input!");
}
else
matrix=NULL;
@@ -93,28 +121,53 @@
void FetchData(int** pmatrix,int* pM,int *pN,PyObject* py_matrix){
/*output: */
- double* dmatrix=NULL;
int* matrix=NULL;
int M,N;
+ int ndim;
+ npy_intp* dims=NULL;
/*intermediary:*/
+ double* dmatrix=NULL;
+ long* lmatrix=NULL;
+ bool* bmatrix=NULL;
int i;
- int ndim;
- npy_intp* dims=NULL;
- /*retrive dimensions: */
+ /*retrieve dimensions: */
ndim=PyArray_NDIM((const PyArrayObject*)py_matrix);
if(ndim!=2)_error_("expecting an MxN matrix in input!");
dims=PyArray_DIMS((PyArrayObject*)py_matrix);
M=dims[0]; N=dims[1];
if (M && N) {
- /*retrieve internal value: */
- dmatrix=(double*)PyArray_DATA((PyArrayObject*)py_matrix);
+ if (PyArray_TYPE((PyArrayObject*)py_matrix) == NPY_DOUBLE) {
+ /*retrieve internal value: */
+ dmatrix=(double*)PyArray_DATA((PyArrayObject*)py_matrix);
- /*transform into integer matrix: */
- matrix=xNew<int>(M*N);
- for(i=0;i<M*N;i++)matrix[i]=(int)dmatrix[i];
+ /*transform into integer matrix: */
+ matrix=xNew<int>(M*N);
+ for(i=0;i<M*N;i++)matrix[i]=(int)dmatrix[i];
+ }
+
+ else if (PyArray_TYPE((PyArrayObject*)py_matrix) == NPY_INT64) {
+ /*retrieve internal value: */
+ lmatrix=(long*)PyArray_DATA((PyArrayObject*)py_matrix);
+
+ /*transform into integer matrix: */
+ matrix=xNew<int>(M*N);
+ for(i=0;i<M*N;i++)matrix[i]=(int)lmatrix[i];
+ }
+
+ else if (PyArray_TYPE((PyArrayObject*)py_matrix) == NPY_BOOL) {
+ /*retrieve internal value: */
+ bmatrix=(bool*)PyArray_DATA((PyArrayObject*)py_matrix);
+
+ /*transform into integer matrix: */
+ matrix=xNew<int>(M*N);
+ for(i=0;i<M*N;i++)matrix[i]=(int)bmatrix[i];
+ }
+
+ else
+ _error_("unrecognized matrix type in input!");
}
else
matrix=NULL;
@@ -135,19 +188,47 @@
int ndim;
npy_intp* dims=NULL;
- /*retrive dimensions: */
+ /*intermediary:*/
+ long* lvector=NULL;
+ bool* bvector=NULL;
+ int i;
+
+ /*retrieve dimensions: */
ndim=PyArray_NDIM((const PyArrayObject*)py_vector);
if(ndim!=1)_error_("expecting an Mx1 vector in input!");
dims=PyArray_DIMS((PyArrayObject*)py_vector);
M=dims[0];
if (M) {
- /*retrieve internal value: */
- dvector=(double*)PyArray_DATA((PyArrayObject*)py_vector);
+ if (PyArray_TYPE((PyArrayObject*)py_vector) == NPY_DOUBLE) {
+ /*retrieve internal value: */
+ dvector=(double*)PyArray_DATA((PyArrayObject*)py_vector);
- /*copy vector: */
- vector=xNew<double>(M);
- memcpy(vector,dvector,(M)*sizeof(double));
+ /*copy vector: */
+ vector=xNew<double>(M);
+ memcpy(vector,dvector,(M)*sizeof(double));
+ }
+
+ else if (PyArray_TYPE((PyArrayObject*)py_vector) == NPY_INT64) {
+ /*retrieve internal value: */
+ lvector=(long*)PyArray_DATA((PyArrayObject*)py_vector);
+
+ /*transform into double vector: */
+ vector=xNew<double>(M);
+ for(i=0;i<M;i++)vector[i]=(double)lvector[i];
+ }
+
+ else if (PyArray_TYPE((PyArrayObject*)py_vector) == NPY_BOOL) {
+ /*retrieve internal value: */
+ bvector=(bool*)PyArray_DATA((PyArrayObject*)py_vector);
+
+ /*transform into double vector: */
+ vector=xNew<double>(M);
+ for(i=0;i<M;i++)vector[i]=(double)bvector[i];
+ }
+
+ else
+ _error_("unrecognized vector type in input!");
}
else
vector=NULL;
Modified: issm/trunk/src/wrappers/python/io/WritePythonData.cpp
===================================================================
--- issm/trunk/src/wrappers/python/io/WritePythonData.cpp 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/wrappers/python/io/WritePythonData.cpp 2012-11-30 18:09:49 UTC (rev 14067)
@@ -42,6 +42,38 @@
PyTuple_SetItem(tuple, index, array);
}/*}}}*/
+/*FUNCTION WriteData(PyObject* py_tuple,int index, int* matrix, int M, int N){{{*/
+void WriteData(PyObject* tuple, int index, int* matrix, int M,int N){
+
+ npy_intp dims[2]={0,0};
+ PyObject* array=NULL;
+
+ /*intermediary:*/
+ long* lmatrix=NULL;
+ int i;
+
+ /*transform into long matrix: */
+ lmatrix=xNew<long>(M*N);
+ for(i=0;i<M*N;i++)lmatrix[i]=(long)matrix[i];
+
+ dims[0]=(npy_intp)M;
+ dims[1]=(npy_intp)N;
+ array=PyArray_SimpleNewFromData(2,dims,NPY_INT64,lmatrix);
+
+ PyTuple_SetItem(tuple, index, array);
+}/*}}}*/
+/*FUNCTION WriteData(PyObject* py_tuple,int index, bool* matrix, int M, int N){{{*/
+void WriteData(PyObject* tuple, int index, bool* matrix, int M,int N){
+
+ npy_intp dims[2]={0,0};
+ PyObject* array=NULL;
+
+ dims[0]=(npy_intp)M;
+ dims[1]=(npy_intp)N;
+ array=PyArray_SimpleNewFromData(2,dims,NPY_BOOL,matrix);
+
+ PyTuple_SetItem(tuple, index, array);
+}/*}}}*/
/*FUNCTION WriteData(PyObject* py_tuple,int index){{{*/
void WriteData(PyObject* py_tuple, int index){
Modified: issm/trunk/src/wrappers/python/io/pythonio.h
===================================================================
--- issm/trunk/src/wrappers/python/io/pythonio.h 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/src/wrappers/python/io/pythonio.h 2012-11-30 18:09:49 UTC (rev 14067)
@@ -17,6 +17,8 @@
#include "../../c/include/include.h"
void WriteData(PyObject* py_tuple,int index, double* matrix, int M,int N);
+void WriteData(PyObject* py_tuple,int index, int* matrix, int M,int N);
+void WriteData(PyObject* py_tuple,int index, bool* matrix, int M,int N);
void WriteData(PyObject* py_tuple,int index, int integer);
void WriteData(PyObject* py_tuple,int index, char* string);
void WriteData(PyObject* py_tuple,int index);
Property changes on: issm/trunk/test
___________________________________________________________________
Modified: svn:mergeinfo
- /issm/branches/trunk-jpl-damage/test:11427-13118
/issm/trunk-jpl/test:10936-11994,11996-12326,12333-12336,12338-12703,12705-12707,12710-13393,13396-13974
+ /issm/branches/trunk-jpl-damage/test:11427-13118
/issm/trunk-jpl/test:10936-11994,11996-12326,12333-12336,12338-12703,12705-12707,12710-13393,13396-13974,13976-14066
Modified: issm/trunk/test/Archives/Archive1205.nc
===================================================================
(Binary files differ)
Modified: issm/trunk/test/Archives/Archive1206.nc
===================================================================
(Binary files differ)
Modified: issm/trunk/test/Archives/Archive1207.nc
===================================================================
(Binary files differ)
Modified: issm/trunk/test/NightlyRun/runme.m
===================================================================
--- issm/trunk/test/NightlyRun/runme.m 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/test/NightlyRun/runme.m 2012-11-30 18:09:49 UTC (rev 14067)
@@ -96,6 +96,8 @@
%Process Ids according to benchmarks{{{
if strcmpi(benchmark,'nightly'),
test_ids=intersect(test_ids,[1:999]);
+elseif strcmpi(benchmark,'validation'),
+ test_ids=intersect(test_ids,[1001:1999]);
elseif strcmpi(benchmark,'ismip'),
test_ids=intersect(test_ids,[1101:1199]);
elseif strcmpi(benchmark,'eismint'),
@@ -104,12 +106,12 @@
test_ids=intersect(test_ids,[1301:1399]);
elseif strcmpi(benchmark,'mesh'),
test_ids=intersect(test_ids,[1401:1499]);
+elseif strcmpi(benchmark,'tranforcing'),
+ test_ids=intersect(test_ids,[1501:1502]);
+elseif strcmpi(benchmark,'referential'),
+ test_ids=intersect(test_ids,[1601:1602]);
elseif strcmpi(benchmark,'adolc'),
test_ids=intersect(test_ids,[3001:3020]);
-elseif strcmpi(benchmark,'validation'),
- test_ids=intersect(test_ids,[1001:1999]);
-elseif strcmpi(benchmark,'tranforcing'),
- test_ids=intersect(test_ids,[1501:1502]);
end
% }}}
Modified: issm/trunk/test/NightlyRun/runme.py
===================================================================
--- issm/trunk/test/NightlyRun/runme.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/test/NightlyRun/runme.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -98,6 +98,8 @@
#Process Ids according to benchmarks {{{
if strcmpi(benchmark,'nightly'):
test_ids=test_ids.intersection(set(range(1,1000)))
+ elif strcmpi(benchmark,'validation'):
+ test_ids=test_ids.intersection(set(range(1001,2000)))
elif strcmpi(benchmark,'ismip'):
test_ids=test_ids.intersection(set(range(1101,1200)))
elif strcmpi(benchmark,'eismint'):
@@ -106,12 +108,12 @@
test_ids=test_ids.intersection(set(range(1301,1400)))
elif strcmpi(benchmark,'mesh'):
test_ids=test_ids.intersection(set(range(1401,1500)))
+ elif strcmpi(benchmark,'tranforcing'):
+ test_ids=test_ids.intersection(set(range(1501,1503)))
+ elif strcmpi(benchmark,'referential'):
+ test_ids=test_ids.intersection(set(range(1601,1603)))
elif strcmpi(benchmark,'adolc'):
test_ids=test_ids.intersection(set(range(3001,3020)))
- elif strcmpi(benchmark,'validation'):
- test_ids=test_ids.intersection(set(range(1001,2000)))
- elif strcmpi(benchmark,'tranforcing'):
- test_ids=test_ids.intersection(set(range(1501,1503)))
#print 'test_ids after benchmark =',test_ids
test_ids=list(test_ids)
test_ids.sort()
Modified: issm/trunk/test/NightlyRun/test109.py
===================================================================
--- issm/trunk/test/NightlyRun/test109.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/test/NightlyRun/test109.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -13,10 +13,10 @@
md.extrude(3,1.)
md=setflowequation(md,'macayeal','all')
md.cluster=generic('name',oshostname(),'np',3)
-md.transient.isdiagnostic=0
-md.transient.isprognostic=0
-md.transient.isthermal=1
-md.transient.isgroundingline=0
+md.transient.isdiagnostic=False
+md.transient.isprognostic=False
+md.transient.isthermal=True
+md.transient.isgroundingline=False
md=solve(md,TransientSolutionEnum())
#Fields and tolerances to track changes
Modified: issm/trunk/test/NightlyRun/test121.py
===================================================================
--- issm/trunk/test/NightlyRun/test121.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/test/NightlyRun/test121.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -14,10 +14,10 @@
md=setflowequation(md,'macayeal','all')
md.cluster=generic('name',oshostname(),'np',3);
md.initialization.waterfraction=numpy.zeros((md.mesh.numberofvertices,1))
-md.transient.isdiagnostic=0
-md.transient.isprognostic=0
-md.transient.isthermal=1
-md.transient.isgroundingline=0
+md.transient.isdiagnostic=False
+md.transient.isprognostic=False
+md.transient.isthermal=True
+md.transient.isgroundingline=False
md.thermal.isenthalpy=1
md=solve(md,TransientSolutionEnum())
Modified: issm/trunk/test/NightlyRun/test1501.m
===================================================================
--- issm/trunk/test/NightlyRun/test1501.m 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/test/NightlyRun/test1501.m 2012-11-30 18:09:49 UTC (rev 14067)
@@ -12,7 +12,7 @@
md.timestepping.final_time=2000.;
%Solve for thinning rate -> -1 * surface mass balance
-smb= 2.*ones(md.mesh.numberofvertices,1);
+smb= 2.*ones(md.mesh.numberofvertices,1);
md.surfaceforcings.mass_balance= smb;
md.basalforcings.melting_rate= smb;
Modified: issm/trunk/test/NightlyRun/test1502.m
===================================================================
--- issm/trunk/test/NightlyRun/test1502.m 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/test/NightlyRun/test1502.m 2012-11-30 18:09:49 UTC (rev 14067)
@@ -13,7 +13,7 @@
md.timestepping.final_time=2000.;
%Solve for thinning rate -> -1 * surface mass balance
-smb= 2.*ones(md.mesh.numberofvertices,1);
+smb= 2.*ones(md.mesh.numberofvertices,1);
md.surfaceforcings.mass_balance= smb;
md.basalforcings.melting_rate= smb;
Modified: issm/trunk/test/NightlyRun/test207.py
===================================================================
--- issm/trunk/test/NightlyRun/test207.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/test/NightlyRun/test207.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -14,10 +14,10 @@
md.extrude(3,1.)
md=setflowequation(md,'macayeal','all')
md.cluster=generic('name',oshostname(),'np',3)
-md.transient.isdiagnostic=0
-md.transient.isprognostic=0
-md.transient.isthermal=1
-md.transient.isgroundingline=0
+md.transient.isdiagnostic=False
+md.transient.isprognostic=False
+md.transient.isthermal=True
+md.transient.isgroundingline=False
md=solve(md,TransientSolutionEnum())
# Fields and tolerances to track changes
Modified: issm/trunk/test/NightlyRun/test228.py
===================================================================
--- issm/trunk/test/NightlyRun/test228.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/test/NightlyRun/test228.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -23,7 +23,7 @@
smb=numpy.hstack((smb,smb*-1.))
md.surfaceforcings.mass_balance=numpy.vstack((smb,[1.5,3.]))
-md.transient.isthermal=0
+md.transient.isthermal=False
md=solve(md,TransientSolutionEnum())
Modified: issm/trunk/test/NightlyRun/test229.py
===================================================================
--- issm/trunk/test/NightlyRun/test229.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/test/NightlyRun/test229.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -23,7 +23,7 @@
smb=numpy.hstack((smb,smb*2.))
md.surfaceforcings.mass_balance=numpy.vstack((smb,[1.5,3.]))
-md.transient.isthermal=0
+md.transient.isthermal=False
md=solve(md,TransientSolutionEnum())
Modified: issm/trunk/test/NightlyRun/test230.py
===================================================================
--- issm/trunk/test/NightlyRun/test230.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/test/NightlyRun/test230.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -24,7 +24,7 @@
smb=numpy.hstack((smb,smb*-1.))
md.surfaceforcings.mass_balance=numpy.vstack((smb,[1.5,3.]))
-md.transient.isthermal=0
+md.transient.isthermal=False
md=solve(md,TransientSolutionEnum())
Modified: issm/trunk/test/NightlyRun/test231.py
===================================================================
--- issm/trunk/test/NightlyRun/test231.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/test/NightlyRun/test231.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -24,7 +24,7 @@
smb=numpy.hstack((smb,smb*2.))
md.surfaceforcings.mass_balance=numpy.vstack((smb,[1.5,3.]))
-md.transient.isthermal=0
+md.transient.isthermal=False
md=solve(md,TransientSolutionEnum())
Modified: issm/trunk/test/NightlyRun/test232.py
===================================================================
--- issm/trunk/test/NightlyRun/test232.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/test/NightlyRun/test232.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -17,10 +17,10 @@
md.thermal.spctemperature=numpy.vstack((numpy.hstack((md.thermal.spctemperature, md.thermal.spctemperature+5., md.thermal.spctemperature+10., md.thermal.spctemperature+15.)), [1.5,2.5,3.5,4.]))
md.timestepping.time_step=1.
md.timestepping.final_time=4.
-md.transient.isdiagnostic=0
-md.transient.isprognostic=0
-md.transient.isthermal=1
-md.transient.isgroundingline=0
+md.transient.isdiagnostic=False
+md.transient.isprognostic=False
+md.transient.isthermal=True
+md.transient.isgroundingline=False
md=solve(md,TransientSolutionEnum())
#Fields and tolerances to track changes
Modified: issm/trunk/test/NightlyRun/test3009.py
===================================================================
--- issm/trunk/test/NightlyRun/test3009.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/test/NightlyRun/test3009.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -13,10 +13,10 @@
md.extrude(3,1.)
md=setflowequation(md,'macayeal','all')
md.cluster=generic('name',oshostname(),'np',3)
-md.transient.isdiagnostic=0
-md.transient.isprognostic=0
-md.transient.isthermal=1
-md.transient.isgroundingline=0
+md.transient.isdiagnostic=False
+md.transient.isprognostic=False
+md.transient.isthermal=True
+md.transient.isgroundingline=False
md.autodiff.isautodiff=True
md=solve(md,TransientSolutionEnum())
Modified: issm/trunk/test/NightlyRun/test313.py
===================================================================
--- issm/trunk/test/NightlyRun/test313.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/test/NightlyRun/test313.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -14,10 +14,10 @@
md=setflowequation(md,'macayeal','all')
md.cluster=generic('name',oshostname(),'np',3)
md.verbose=verbose('convergence',True,'solution',True)
-md.transient.isdiagnostic=0
-md.transient.isprognostic=0
-md.transient.isthermal=1
-md.transient.isgroundingline=0
+md.transient.isdiagnostic=False
+md.transient.isprognostic=False
+md.transient.isthermal=True
+md.transient.isgroundingline=False
md=solve(md,TransientSolutionEnum())
#Fields and tolerances to track changes
Modified: issm/trunk/test/NightlyRun/test326.py
===================================================================
--- issm/trunk/test/NightlyRun/test326.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/test/NightlyRun/test326.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -15,10 +15,10 @@
md=setflowequation(md,'macayeal','all')
md.cluster=generic('name',oshostname(),'np',3)
md.initialization.waterfraction=numpy.zeros((md.mesh.numberofvertices,1))
-md.transient.isdiagnostic=0
-md.transient.isprognostic=0
-md.transient.isthermal=1
-md.transient.isgroundingline=0
+md.transient.isdiagnostic=False
+md.transient.isprognostic=False
+md.transient.isthermal=True
+md.transient.isgroundingline=False
md.thermal.isenthalpy=1
md=solve(md,TransientSolutionEnum())
Modified: issm/trunk/test/NightlyRun/test407.py
===================================================================
--- issm/trunk/test/NightlyRun/test407.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/test/NightlyRun/test407.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -14,10 +14,10 @@
md.extrude(4,1.)
md=setflowequation(md,'pattyn','all')
md.cluster=generic('name',oshostname(),'np',3)
-md.transient.isdiagnostic=0
-md.transient.isprognostic=0
-md.transient.isthermal=1
-md.transient.isgroundingline=0
+md.transient.isdiagnostic=False
+md.transient.isprognostic=False
+md.transient.isthermal=True
+md.transient.isgroundingline=False
md=solve(md,TransientSolutionEnum())
#Fields and tolerances to track changes
Modified: issm/trunk/test/NightlyRun/test423.py
===================================================================
--- issm/trunk/test/NightlyRun/test423.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/test/NightlyRun/test423.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -28,10 +28,10 @@
md=setflowequation(md,'macayeal','all')
md.cluster=generic('name',oshostname(),'np',3)
-md.transient.isthermal=0
-md.transient.isprognostic=0
-md.transient.isdiagnostic=0
-md.transient.isgroundingline=1
+md.transient.isthermal=False
+md.transient.isprognostic=False
+md.transient.isdiagnostic=False
+md.transient.isgroundingline=True
#test different grounding line dynamics.
md.groundingline.migration='AgressiveMigration'
Modified: issm/trunk/test/NightlyRun/test424.py
===================================================================
--- issm/trunk/test/NightlyRun/test424.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/test/NightlyRun/test424.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -19,8 +19,8 @@
md.geometry.thickness[:]=1000.
md.geometry.surface=md.geometry.bed+md.geometry.thickness
md.surfaceforcings.mass_balance[:]=100.
-md.transient.isdiagnostic=0
-md.transient.isgroundingline=1
+md.transient.isdiagnostic=False
+md.transient.isgroundingline=True
md.groundingline.migration='AgressiveMigration'
md.cluster=generic('name',oshostname(),'np',3)
Modified: issm/trunk/test/NightlyRun/test425.py
===================================================================
--- issm/trunk/test/NightlyRun/test425.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/test/NightlyRun/test425.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -19,8 +19,8 @@
md.geometry.thickness[:]=1300.
md.geometry.surface=md.geometry.bed+md.geometry.thickness
md.surfaceforcings.mass_balance[:]=-150.
-md.transient.isdiagnostic=0
-md.transient.isgroundingline=1
+md.transient.isdiagnostic=False
+md.transient.isgroundingline=True
md.groundingline.migration='SoftMigration'
md.cluster=generic('name',oshostname(),'np',3)
Modified: issm/trunk/test/NightlyRun/test426.py
===================================================================
--- issm/trunk/test/NightlyRun/test426.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/test/NightlyRun/test426.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -20,8 +20,8 @@
md.surfaceforcings.mass_balance[:]=100.
md.extrude(3,1.);
md=setflowequation(md,'macayeal','all');
-md.transient.isdiagnostic=0
-md.transient.isgroundingline=1
+md.transient.isdiagnostic=False
+md.transient.isgroundingline=True
md.groundingline.migration='AgressiveMigration'
md.cluster=generic('name',oshostname(),'np',3)
Modified: issm/trunk/test/NightlyRun/test427.py
===================================================================
--- issm/trunk/test/NightlyRun/test427.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/test/NightlyRun/test427.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -21,8 +21,8 @@
md.extrude(3,1.)
md.surfaceforcings.mass_balance[:]=-150
-md.transient.isdiagnostic=0
-md.transient.isgroundingline=1
+md.transient.isdiagnostic=False
+md.transient.isgroundingline=True
md.groundingline.migration='SoftMigration'
md.cluster=generic('name',oshostname(),'np',3)
md=solve(md,TransientSolutionEnum())
Modified: issm/trunk/test/NightlyRun/test515.py
===================================================================
--- issm/trunk/test/NightlyRun/test515.py 2012-11-30 18:06:27 UTC (rev 14066)
+++ issm/trunk/test/NightlyRun/test515.py 2012-11-30 18:09:49 UTC (rev 14067)
@@ -14,10 +14,10 @@
md=setflowequation(md,'pattyn','all')
md.thermal.stabilization=2
md.cluster=generic('name',oshostname(),'np',3)
-md.transient.isdiagnostic=0
-md.transient.isprognostic=0
-md.transient.isthermal=1
-md.transient.isgroundingline=0
+md.transient.isdiagnostic=False
+md.transient.isprognostic=False
+md.transient.isthermal=True
+md.transient.isgroundingline=False
md=solve(md,TransientSolutionEnum())
# Fields and tolerances to track changes
More information about the issm-svn
mailing list