Changeset 14067 for issm/trunk
- Timestamp:
- 11/30/12 10:09:49 (12 years ago)
- Location:
- issm/trunk
- Files:
-
- 16 deleted
- 102 edited
- 6 copied
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk
- Property svn:mergeinfo changed
/issm/trunk-jpl merged: 13976-13979,13982-14022,14026,14028,14034,14046-14063,14066
- Property svn:mergeinfo changed
-
issm/trunk/aux-config
- Property svn:ignore
-
old new 1 1 ar-lib 2 depcomp 3 compile 4 missing 5 config.guess 6 config.sub 7 ltmain.sh 8 install-sh
-
- Property svn:ignore
-
issm/trunk/configs/config-greenplanet.sh
r13975 r14067 22 22 --with-blacs-dir=$ISSM_DIR/externalpackages/petsc/install/ \ 23 23 --with-graphics-lib=/usr/lib64/libX11.so \ 24 --with-cxxoptflags="-O3 -xS" \24 --with-cxxoptflags="-O3" \ 25 25 --with-vendor=intel-linux -
issm/trunk/configure.ac
r13975 r14067 2 2 3 3 #AUTOCONF 4 AC_INIT([ISSM],[4.2. 3],[issm@jpl.nasa.gov],[issm],[http://issm.jpl.nasa.gov]) #Initializing configure4 AC_INIT([ISSM],[4.2.4],[issm@jpl.nasa.gov],[issm],[http://issm.jpl.nasa.gov]) #Initializing configure 5 5 AC_CONFIG_AUX_DIR([./aux-config]) #Put config files in aux-config 6 6 AC_CONFIG_MACRO_DIR([m4]) #m4 macros are located in m4 -
issm/trunk/etc/environment.sh
r13975 r14067 86 86 87 87 DOXYGEN_DIR="$ISSM_DIR/externalpackages/doxygen/install" 88 path append "$DOXYGEN_DIR/bin"88 pathprepend "$DOXYGEN_DIR/bin" 89 89 90 90 AUTOTOOLS_DIR="$ISSM_DIR/externalpackages/autotools/install" -
issm/trunk/externalpackages/chaco
- Property svn:ignore
-
old new 1 *.pdf 1 2 install 2 3 src 3 4 .ignore.txt 4 5 old 6 *.tar.gz
-
- Property svn:ignore
-
issm/trunk/externalpackages/doxygen/install.sh
r13395 r14067 3 3 4 4 #Some cleanup 5 rm -rf install 5 rm -rf install src 6 6 7 7 #Download latest version 8 svn co https://doxygen.svn.sourceforge.net/svnroot/doxygen/trunk install8 svn co https://doxygen.svn.sourceforge.net/svnroot/doxygen/trunk src 9 9 10 10 #Configure doxygen -
issm/trunk/externalpackages/gsl/install-android.sh
r13975 r14067 25 25 cd src 26 26 27 mv ./../Makefile.am.patch ./ 27 28 patch Makefile.am < Makefile.am.patch 28 29 29 30 autoreconf -iv --force -I $ISSM_DIR/externalpackages/autotools/install/share/aclocal 30 31 … … 38 39 #Compile gsl 39 40 if [[ $step == "3" || $step == "0" ]]; then 40 cd src 41 cd $ISSM_DIR/externalpackages/gsl/src 42 41 43 if [ $# -eq 0 ]; then 42 make $j44 make 43 45 else 44 46 make -j $j -
issm/trunk/externalpackages/triangle/configs/android/configure.make
r13395 r14067 10 10 # http://www.codesourcery.com/gnu_toolchains/arm/arm_gnu_linux_abi.pdf 11 11 12 source $ANDROID_DIR/android_aux.sh 13 14 CC=${toolchain_path}-gcc 15 AR=${toolchain_path}-ar 16 RANLIB=${toolchain_path}-ranlib 12 CC=${host_triplet}-gcc 13 AR=${host_triplet}-ar 14 RANLIB=${host_triplet}-ranlib 17 15 CSWITCHES = $(CFLAGS) 18 16 TRILIBDEFS = -DTRILIBRARY -
issm/trunk/externalpackages/triangle/install-android.sh
r13395 r14067 1 1 #!/bin/bash 2 2 set -eu 3 4 source $ANDROID_DIR/android_aux.sh 3 5 4 6 #Some cleanup -
issm/trunk/src
- Property svn:mergeinfo changed
/issm/trunk-jpl/src merged: 13976-13977,13983-13984,13988-14002,14004-14020,14026,14028,14034,14046-14059,14062
- Property svn:mergeinfo changed
-
issm/trunk/src/android/ISSM
- Property svn:ignore
-
old new 1 Makefile.in2 Makefile 1 gen 2 .settings 3 3 libs 4 4 obj 5 Makefile 6 Makefile.in
-
- Property svn:ignore
-
issm/trunk/src/android/ISSM/.classpath
r13890 r14067 1 1 <?xml version="1.0" encoding="UTF-8"?> 2 2 <classpath> 3 <classpathentry kind="con" path="com.android.ide.eclipse.adt.ANDROID_FRAMEWORK"/> 4 <classpathentry kind="con" path="com.android.ide.eclipse.adt.LIBRARIES"/> 3 5 <classpathentry kind="src" path="src"/> 4 6 <classpathentry kind="src" path="gen"/> 5 <classpathentry kind="con" path="com.android.ide.eclipse.adt.ANDROID_FRAMEWORK"/>6 <classpathentry kind="con" path="com.android.ide.eclipse.adt.LIBRARIES"/>7 7 <classpathentry kind="output" path="bin/classes"/> 8 8 </classpath> -
issm/trunk/src/android/ISSM/AndroidManifest.xml
r13890 r14067 5 5 6 6 <uses-sdk 7 android:minSdkVersion="1 5"7 android:minSdkVersion="10" 8 8 android:targetSdkVersion="15" /> 9 9 <uses-permission android:name="android.permission.WRITE_EXTERNAL_STORAGE" /> 10 10 <application 11 11 android:icon="@drawable/ic_launcher" 12 12 android:label="@string/app_name" 13 13 android:theme="@style/AppTheme" > 14 <activity 14 <activity android:name=".MapSelection" 15 android:label="@string/title_activity_issm" > 16 <intent-filter> 17 <action android:name="android.intent.action.MAIN" /> 18 <category android:name="android.intent.category.LAUNCHER" /> 19 </intent-filter> 20 </activity> 21 <activity 15 22 android:name=".ISSM" 16 23 android:label="@string/title_activity_issm" > 17 <intent-filter>18 <action android:name="android.intent.action.MAIN" />19 20 <category android:name="android.intent.category.LAUNCHER" />21 </intent-filter>22 24 </activity> 23 25 </application> -
issm/trunk/src/android/ISSM/bin
- Property svn:ignore
-
old new 1 ISSM.apk1 *ISSM.apk 2 2 resources.ap_ 3 3 res
-
- Property svn:ignore
-
issm/trunk/src/android/ISSM/bin/AndroidManifest.xml
r13891 r14067 5 5 6 6 <uses-sdk 7 android:minSdkVersion="1 5"7 android:minSdkVersion="10" 8 8 android:targetSdkVersion="15" /> 9 9 <uses-permission android:name="android.permission.WRITE_EXTERNAL_STORAGE" /> 10 10 <application 11 11 android:icon="@drawable/ic_launcher" 12 12 android:label="@string/app_name" 13 13 android:theme="@style/AppTheme" > 14 <activity 14 <activity android:name=".MapSelection" 15 android:label="@string/title_activity_issm" > 16 <intent-filter> 17 <action android:name="android.intent.action.MAIN" /> 18 <category android:name="android.intent.category.LAUNCHER" /> 19 </intent-filter> 20 </activity> 21 <activity 15 22 android:name=".ISSM" 16 23 android:label="@string/title_activity_issm" > 17 <intent-filter>18 <action android:name="android.intent.action.MAIN" />19 20 <category android:name="android.intent.category.LAUNCHER" />21 </intent-filter>22 24 </activity> 23 25 </application> -
issm/trunk/src/android/ISSM/gen/com/example/issm/R.java
r13890 r14067 22 22 public static final class id { 23 23 public static final int button1=0x7f080001; 24 public static final int button2=0x7f080003; 24 25 public static final int input=0x7f080000; 25 public static final int menu_settings=0x7f08000 3;26 public static final int menu_settings=0x7f080004; 26 27 public static final int output=0x7f080002; 27 28 } 28 29 public static final class layout { 29 30 public static final int activity_issm=0x7f030000; 31 public static final int activity_mapselection=0x7f030001; 30 32 } 31 33 public static final class menu { -
issm/trunk/src/android/ISSM/jni/Android.mk
r13931 r14067 11 11 LOCAL_PATH := $(my_LOCAL_PATH) 12 12 LOCAL_ALLOW_UNDEFINED_SYMBOLS := true 13 LOCAL_CFLAGS := -Wno-psabi 13 LOCAL_CFLAGS := -Wno-psabi -DHAVE_CONFIG_H 14 14 LOCAL_LDLIBS := -llog -ldl 15 15 LOCAL_MODULE := IssmJni -
issm/trunk/src/android/ISSM/jni/Main.cpp
r13959 r14067 1 1 #include <jni.h> 2 2 #include "../../../c/android/fac.h" 3 #include "../../../c/issm.h" 4 #include <cstddef> 5 #include <stdio.h> 6 /////////////////////////////////////////////////////////////////////////////////////////// 3 7 namespace com_example_issm 4 8 { 5 static jlong factorial(JNIEnv *env, jclass clazz, jlong n) 9 fac* f; 10 FemModel *fm; 11 //------------------------------------------------------------------------------------ 12 jint initilize(JNIEnv *env, jclass clazz, jstring file) 6 13 { 7 fac *f = new fac(); 8 jlong result = (jlong) (f->factorial(n)); 9 delete(f); 10 return (jlong) result; 14 char *nativefile = (char*)env->GetStringUTFChars(file,0); 15 16 f = new fac(); 17 18 //call Model constructor passing in infile as File Descriptor parameter. 19 fm = new FemModel(nativefile); 20 21 env->ReleaseStringUTFChars(file, nativefile); //must realease the char* 22 //jint size = (jint) fm->bufferSize(); 23 //return size; 24 25 return 0; //return 0 for now. 11 26 } 12 27 //------------------------------------------------------------------------------------ 13 static JNINativeMethod method_table[] = { 14 {"fac" ,"(J)J" , (void *) factorial}, 28 //fill out the first two slots, extract the same way in java using get(int position) 29 void solve(JNIEnv *env, jclass clazz , jint alpha, jobject buf) 30 { 31 jdouble *dBuf = (jdouble *)env->GetDirectBufferAddress(buf); 32 33 //pass bBuff to fem model to allocate data 34 // fm -> solve(dBuf, alpha); 35 36 } 37 //------------------------------------------------------------------------------------ 38 jlong factorial(JNIEnv *env, jclass clazz, jlong n) 39 { 40 if( f != NULL) 41 return (jlong) (f->factorial(n)); 42 return 0; 43 } 44 45 //------------------------------------------------------------------------------------ 46 static JNINativeMethod method_table[] = 47 { 48 {"fac" , "(J)J" , (void *) factorial}, 49 {"createISSMModel" ,"(Ljava/lang/String;)I" , (void *) initilize}, 50 {"processBuffer", "(ILjava/nio/DoubleBuffer;)V", (void *) solve} 15 51 }; 16 52 } 17 53 //------------------------------------------------------------------------------------ 18 54 using namespace com_example_issm; 19 55 … … 29 65 if(clazz) 30 66 { 31 env->RegisterNatives(clazz, method_table, 1);67 env->RegisterNatives(clazz, method_table, 3); 32 68 env->DeleteLocalRef(clazz); 33 69 return JNI_VERSION_1_6; … … 35 71 else return -1; 36 72 } 37 38 // Get jclass with env->FindClass.39 // Register methods with env->RegisterNatives.40 73 } 74 /////////////////////////////////////////////////////////////////////////////////////////// -
issm/trunk/src/android/ISSM/jni/issmlib/Android.mk
r13931 r14067 2 2 include $(CLEAR_VARS) 3 3 LOCAL_MODULE := libISSMCore 4 LOCAL_SRC_FILES := ../../../../c/libISSMCore.a5 LOCAL_EXPORT_C_INCLUDES := ../../../../c/android/4 LOCAL_SRC_FILES := libISSMCore.a 5 LOCAL_EXPORT_C_INCLUDES := $(ISSM_DIR) 6 6 include $(PREBUILT_STATIC_LIBRARY) -
issm/trunk/src/android/ISSM/src/com/example/issm/ISSM.java
r13931 r14067 1 1 package com.example.issm; 2 3 import java.io.IOException; 4 import java.io.InputStream; 5 import java.nio.ByteBuffer; 6 import java.nio.ByteOrder; 7 import java.nio.DoubleBuffer; 2 8 3 9 import android.os.Bundle; 4 10 import android.app.Activity; 11 import android.content.res.AssetManager; 5 12 import android.text.TextUtils; 6 13 import android.view.Menu; … … 12 19 13 20 14 public class ISSM extends Activity implements OnClickListener { 21 public class ISSM extends Activity implements OnClickListener 22 { 15 23 private EditText input; 16 24 private TextView output; 25 private DoubleBuffer buff; 26 private IssmJni issmNative; 27 private String mapName; 28 private String issmFolder; 29 private int size; 17 30 @Override 31 //------------------------------------------------------------------------------------------------ 18 32 public void onCreate(Bundle savedInstanceState) { 19 33 super.onCreate(savedInstanceState); 34 Bundle map = getIntent().getExtras(); 35 { 36 if(map!= null) 37 { 38 mapName = map.getString("map"); 39 issmFolder = map.getString("pathToFile"); 40 } 41 } 20 42 setContentView(R.layout.activity_issm); 21 43 this.input = (EditText) super.findViewById(R.id.input); … … 24 46 button.setOnClickListener(this); 25 47 48 //load up the ISSM library and create double buffer in java 49 //which later on will be pass for native allocation. 50 issmNative = new IssmJni(); 51 buff = ByteBuffer.allocateDirect(15*8).order(ByteOrder.nativeOrder()).asDoubleBuffer(); 52 this.createModel(); 26 53 } 27 28 public void onClick(View view) 54 //------------------------------------------------------------------------------------------------ 55 public void createModel() 56 { 57 String file = ""; 58 if( mapName.equals("greenland")) 59 { 60 file = "test102.petsc"; 61 } 62 else file = "test102.petsc"; 63 size = issmNative.createISSMModel(issmFolder + file); 64 } 65 //------------------------------------------------------------------------------------------------ 66 public void fillBuffer() 67 { 68 issmNative.processBuffer(5,buff); 69 } 70 //------------------------------------------------------------------------------------------------ 71 public void onClick(View view) 29 72 { 30 // TODO Auto-generated method stub73 //factorial method 31 74 String input = this.input.getText().toString(); 32 75 if(TextUtils.isEmpty(input)) 33 76 { 34 77 return; 35 } 36 long result = IssmJni.facIterative(Long.parseLong(input)); 37 this.output.setText("Result = " + result + "\n"); 78 } 79 long resultfromFac = issmNative.fac(Long.parseLong(input)); 80 //example of how to fill buffer Native 81 this.fillBuffer(); 82 83 //print result from fac and the first two slot of filled buffer. 84 this.output.setText("Result = " + resultfromFac + "\n" 85 + "First slot from buffer:\n" 86 + buff.get(0) + " size " + size); 87 38 88 } 39 89 } -
issm/trunk/src/android/ISSM/src/com/example/issm/IssmJni.java
r13931 r14067 1 1 package com.example.issm; 2 import java.nio.DoubleBuffer; 2 3 3 publicclass IssmJni4 class IssmJni 4 5 { 5 public static native long fac(long n); 6 static { 6 public native long fac(long n); 7 public native void processBuffer(int alpha, DoubleBuffer buff); 8 public native int createISSMModel(String file); 9 static 10 { 7 11 System.loadLibrary("IssmJni"); 8 12 } 9 public static long facIterative(long n)10 {11 return fac(n);12 }13 13 } -
issm/trunk/src/android/ISSM_Visual/bin
-
Property svn:ignore
set to
*.apk
-
Property svn:ignore
set to
-
issm/trunk/src/android/ISSM_Visual/src/com/example/issm_visual
-
Property svn:ignore
set to
*.java
*.class
-
Property svn:ignore
set to
-
issm/trunk/src/c/Container/Observations.cpp
r13975 r14067 424 424 *pprediction = obs; 425 425 }/*}}}*/ 426 /*FUNCTION Observations::InterpolationV4{{{*/ 427 void Observations::InterpolationV4(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata){ 428 /* Reference: David T. Sandwell, Biharmonic spline interpolation of GEOS-3 429 * and SEASAT altimeter data, Geophysical Research Letters, 2, 139-142, 430 * 1987. Describes interpolation using value or gradient of value in any 431 * dimension.*/ 432 433 /*Intermediaries*/ 434 int i,j,n_obs; 435 IssmPDouble prediction,h; 436 IssmPDouble *x = NULL; 437 IssmPDouble *y = NULL; 438 IssmPDouble *obs = NULL; 439 IssmPDouble *Green = NULL; 440 IssmPDouble *weights = NULL; 441 IssmPDouble *g = NULL; 442 443 /*Some checks*/ 444 _assert_(maxdata>0); 445 _assert_(pprediction); 446 447 /*If radius is not provided or is 0, return all observations*/ 448 if(radius==0) radius=this->quadtree->root->length; 449 450 /*Get list of observations for current point*/ 451 this->ObservationList(&x,&y,&obs,&n_obs,x_interp,y_interp,radius,maxdata); 452 453 /*If we have less observations than mindata, return UNDEF*/ 454 if(n_obs<mindata || n_obs<2){ 455 prediction = UNDEF; 456 } 457 else{ 458 459 /*Allocate intermediary matrix and vectors*/ 460 Green = xNew<IssmPDouble>(n_obs*n_obs); 461 g = xNew<IssmPDouble>(n_obs); 462 463 /*First: distance vector*/ 464 for(i=0;i<n_obs;i++){ 465 h = sqrt( (x[i]-x_interp)*(x[i]-x_interp) + (y[i]-y_interp)*(y[i]-y_interp) ); 466 if(h>0){ 467 g[i] = h*h*(log(h)-1.); 468 } 469 else{ 470 g[i] = 0.; 471 } 472 } 473 474 /*Build Green function matrix*/ 475 for(i=0;i<n_obs;i++){ 476 for(j=0;j<=i;j++){ 477 h = sqrt( (x[i]-x[j])*(x[i]-x[j]) + (y[i]-y[j])*(y[i]-y[j]) ); 478 if(h>0){ 479 Green[j*n_obs+i] = h*h*(log(h)-1.); 480 } 481 else{ 482 Green[j*n_obs+i] = 0.; 483 } 484 Green[i*n_obs+j] = Green[j*n_obs+i]; 485 } 486 } 487 488 /*Compute weights*/ 489 #if _HAVE_GSL_ 490 SolverxSeq(&weights,Green,obs,n_obs); // Green^-1 obs 491 #else 492 _error_("GSL is required"); 493 #endif 494 495 /*Interpolate*/ 496 prediction = 0; 497 for(i=0;i<n_obs;i++) prediction += weights[i]*g[i]; 498 499 } 500 501 /*clean-up*/ 502 *pprediction = prediction; 503 xDelete<IssmPDouble>(x); 504 xDelete<IssmPDouble>(y); 505 xDelete<IssmPDouble>(obs); 506 xDelete<IssmPDouble>(Green); 507 xDelete<IssmPDouble>(g); 508 xDelete<IssmPDouble>(weights); 509 }/*}}}*/ 426 510 /*FUNCTION Observations::QuadtreeColoring{{{*/ 427 511 void Observations::QuadtreeColoring(IssmPDouble* A,IssmPDouble *x,IssmPDouble *y,int n){ -
issm/trunk/src/c/Container/Observations.h
r13975 r14067 28 28 void ClosestObservation(IssmDouble *px,IssmDouble *py,IssmDouble *pobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius); 29 29 void InterpolationIDW(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata,IssmDouble power); 30 void InterpolationV4(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata); 30 31 void InterpolationKriging(IssmDouble *pprediction,IssmDouble *perror,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata,Variogram* variogram); 31 32 void InterpolationNearestNeighbor(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius); -
issm/trunk/src/c/Container/Options.h
r13975 r14067 64 64 } 65 65 else{ 66 if(GetOption(name)) _printLine_("WARNING: option "<<name<<" found but fetched format not consistent, defaulting..."); 66 67 *pvalue=default_value; 67 68 } -
issm/trunk/src/c/android/fac.cpp
r13709 r14067 13 13 14 14 /*}}}*/ 15 15 #include "stdio.h" 16 16 long fac::factorial(long n) { 17 17 long f = 1; 18 18 long i; 19 printf("ok1\n"); 19 20 20 21 for(i = 1; i <= n; i++){ 21 22 f *= i; 22 23 } 24 printf("ok2\n"); 23 25 24 26 return f; -
issm/trunk/src/c/classes/FemModel.cpp
r13975 r14067 20 20 21 21 /*Object constructors and destructor*/ 22 /*FUNCTION FemModel::FemModel(char* filename) {{{*/ 23 FemModel::FemModel(char* filename){ 24 25 } 26 /*}}}*/ 22 27 /*FUNCTION FemModel::FemModel(int argc,char** argv){{{*/ 23 28 FemModel::FemModel(int argc,char** argv,COMM incomm){ … … 502 507 503 508 /*Loop over elements that hold node number i*/ 504 _assert_(head_e[node->Sid()]!=-1 || head_l[node->Sid()]!=-1); 509 //if(head_e[node->Sid()]==-1 && head_l[node->Sid()]==-1){ 510 // printf("[%i] vertex %i\n",IssmComm::GetRank(),node->Sid()+1); 511 //} 505 512 for(j=head_e[node->Sid()];j!=-1;j=next_e[j]){ 506 513 offset=count2offset_e[j]; -
issm/trunk/src/c/classes/FemModel.h
r13975 r14067 47 47 48 48 /*constructors, destructors: */ 49 FemModel(char* filename); 49 50 FemModel(int argc,char** argv,COMM comm_init); 50 51 FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, const int solution_type,const int* analyses,const int nummodels); -
issm/trunk/src/c/classes/objects/Elements/Penta.cpp
r13975 r14067 66 66 67 67 /*Build neighbors list*/ 68 if (xIsNan<IssmDouble>(iomodel->Data(MeshUpperelementsEnum)[index]) ) penta_elements_ids[1]=this->id; //upper penta is the same penta68 if (xIsNan<IssmDouble>(iomodel->Data(MeshUpperelementsEnum)[index]) || iomodel->Data(MeshUpperelementsEnum)[index]==-1.) penta_elements_ids[1]=this->id; //upper penta is the same penta 69 69 else penta_elements_ids[1]=reCast<int,IssmDouble>((iomodel->Data(MeshUpperelementsEnum)[index])); 70 if (xIsNan<IssmDouble>(iomodel->Data(MeshLowerelementsEnum)[index]) ) penta_elements_ids[0]=this->id; //lower penta is the same penta70 if (xIsNan<IssmDouble>(iomodel->Data(MeshLowerelementsEnum)[index]) || iomodel->Data(MeshLowerelementsEnum)[index]==-1.) penta_elements_ids[0]=this->id; //lower penta is the same penta 71 71 else penta_elements_ids[0]=reCast<int,IssmDouble>((iomodel->Data(MeshLowerelementsEnum)[index])); 72 72 this->InitHookNeighbors(penta_elements_ids); -
issm/trunk/src/c/modules/ElementConnectivityx/ElementConnectivityx.cpp
r13975 r14067 15 15 #include "../../EnumDefinitions/EnumDefinitions.h" 16 16 17 int hascommondedge( double* element1,double* element2);17 int hascommondedge(int* element1,int* element2); 18 18 19 void ElementConnectivityx( double** pelementconnectivity, double* elements, int nel, double* nodeconnectivity, int nods, int width){19 void ElementConnectivityx(int** pelementconnectivity,int* elements, int nels,int* nodeconnectivity, int nods, int width){ 20 20 21 21 int i,j,k,n; … … 23 23 /*intermediary: */ 24 24 int maxels; 25 double element; 26 double connectedelement; 25 int connectedelement; 27 26 int connectedelementindex; 28 27 int node; … … 30 29 int num_elements; 31 30 32 /*output: */33 double* elementconnectivity=NULL;34 35 31 /*maxels: */ 36 32 maxels=width-1; 33 37 34 /*Allocate connectivity: */ 38 elementconnectivity=xNewZeroInit<double>(nel*3);35 int* elementconnectivity=xNewZeroInit<int>(nels*3); 39 36 40 37 /*Go through all elements, and for each element, go through its nodes, to get the neighbouring elements. 41 38 * Once we get the neighbouring elements, figure out if they share a segment with the current element. If so, 42 39 * plug them in the connectivity, unless they are already there.: */ 40 for(n=0;n<nels;n++){ 43 41 44 for(n=0;n<nel;n++){ 45 46 element=(double)(n+1); //matlab indexing 42 //element=n+1; //matlab indexing 47 43 48 44 for(i=0;i<3;i++){ 49 45 50 node= (int)*(elements+n*3+i); //already matlab indexed, elements comes directly from the workspace.46 node=elements[n*3+i]; //already matlab indexed, elements comes directly from the workspace. 51 47 index=node-1; 52 48 53 num_elements= (int)*(nodeconnectivity+width*index+maxels); //retrieve number of elements already plugged into the connectivity of this node.49 num_elements=nodeconnectivity[width*index+maxels]; //retrieve number of elements already plugged into the connectivity of this node. 54 50 55 51 for(j=0;j<num_elements;j++){ 56 52 57 53 /*for each element connected to node, figure out if it has a commond edge with element: */ 58 connectedelement= *(nodeconnectivity+width*index+j);59 connectedelementindex= (int)(connectedelement-1); //go from matlab indexing to c indexing.54 connectedelement=nodeconnectivity[width*index+j]; 55 connectedelementindex=connectedelement-1; //go from matlab indexing to c indexing. 60 56 61 if(hascommondedge( elements+n*3+0,elements+connectedelementindex*3+0)){57 if(hascommondedge(&elements[n*3+0],&elements[connectedelementindex*3+0])){ 62 58 /*Ok, this connected element has a commond edge with element, plug it into elementconnectivity, unless 63 59 *it is already there: */ 64 60 65 61 for(k=0;k<3;k++){ 66 if (*(elementconnectivity+3*n+k)==0){67 *(elementconnectivity+3*n+k)=connectedelement;62 if(elementconnectivity[3*n+k]==0){ 63 elementconnectivity[3*n+k]=connectedelement; 68 64 break; 69 65 } 70 66 else{ 71 if(connectedelement== *(elementconnectivity+3*n+k))break;67 if(connectedelement==elementconnectivity[3*n+k]) break; 72 68 } 73 69 } … … 81 77 } 82 78 83 int hascommondedge( double* element1,double* element2){79 int hascommondedge(int* el1,int* el2){ 84 80 85 int i,j; 86 int count; 87 88 count=0; 89 for(i=0;i<3;i++){ 90 for(j=0;j<3;j++){ 91 if (*(element1+i)==*(element2+j))count++; 81 int count=0; 82 for(int i=0;i<3;i++){ 83 for(int j=0;j<3;j++){ 84 if(el1[i]==el2[j]) count++; 92 85 } 93 86 } 94 if (count==2)return 1; 95 else return 0; 87 if(count==2) 88 return 1; 89 else 90 return 0; 96 91 } -
issm/trunk/src/c/modules/ElementConnectivityx/ElementConnectivityx.h
r13975 r14067 7 7 8 8 /* local prototypes: */ 9 void ElementConnectivityx( double** pelementconnectivity, double* elements, int nel, double* nodeconnectivity, int nods, int width);9 void ElementConnectivityx(int** pelementconnectivity,int* elements,int nels,int* nodeconnectivity, int nods, int width); 10 10 11 11 #endif /* _ELEMENTCONNECTIVITYX_H */ -
issm/trunk/src/c/modules/Krigingx/Krigingx.cpp
r13395 r14067 22 22 /*Intermediaries*/ 23 23 int mindata,maxdata; 24 double dmindata,dmaxdata,dnumthreads; //FIXME (Options come as double but we want to retrive integers) 24 25 double radius; 25 26 char *output = NULL; … … 37 38 ProcessVariogram(&variogram,options); 38 39 options->Get(&radius,"searchradius",0.); 39 options->Get(&mindata,"mindata",1); 40 options->Get(&maxdata,"maxdata",50); 40 options->Get(&dmindata,"mindata",1.); mindata=int(dmindata);//FIXME (Options come as double but we want to retrive integers) 41 options->Get(&dmaxdata,"maxdata",50.); maxdata=int(dmaxdata);//FIXME (Options come as double but we want to retrive integers) 42 options->Get(&dnumthreads,"numthreads",double(num)); num=int(dnumthreads);//FIXME (Options come as double but we want to retrive integers) 41 43 42 44 /*Process observation dataset*/ … … 88 90 gate.predictions = predictions; 89 91 gate.error = error; 90 gate. percent = xNewZeroInit<double>(num);92 gate.numdone = xNewZeroInit<int>(num); 91 93 92 94 /*launch the thread manager with Krigingxt as a core: */ 93 95 LaunchThread(NearestNeighbort,(void*)&gate,num); 94 96 _printLine_("\r interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%"); 95 xDelete< double>(gate.percent);97 xDelete<int>(gate.numdone); 96 98 } 97 99 else if(strcmp(output,"idw")==0){ //Inverse distance weighting … … 109 111 gate.predictions = predictions; 110 112 gate.error = error; 111 gate. percent = xNewZeroInit<double>(num);113 gate.numdone = xNewZeroInit<int>(num); 112 114 gate.power = power; 113 115 … … 115 117 LaunchThread(idwt,(void*)&gate,num); 116 118 _printLine_("\r interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%"); 117 xDelete<double>(gate.percent); 118 } 119 else if(strcmp(output,"prediction")==0){ 120 119 xDelete<int>(gate.numdone); 120 } 121 else if(strcmp(output,"v4")==0){ //Inverse distance weighting 122 #if !defined(_HAVE_GSL_) 123 _error_("GSL is required for v4 interpolation"); 124 #endif 121 125 /*initialize thread parameters: */ 122 126 gate.n_interp = n_interp; … … 130 134 gate.predictions = predictions; 131 135 gate.error = error; 132 gate.percent = xNewZeroInit<double>(num); 136 gate.numdone = xNewZeroInit<int>(num); 137 138 /*launch the thread manager with Krigingxt as a core: */ 139 LaunchThread(v4t,(void*)&gate,num); 140 _printLine_("\r interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%"); 141 xDelete<int>(gate.numdone); 142 } 143 else if(strcmp(output,"prediction")==0){ 144 #if !defined(_HAVE_GSL_) 145 _error_("GSL is required for v4 interpolation"); 146 #endif 147 148 /*initialize thread parameters: */ 149 gate.n_interp = n_interp; 150 gate.x_interp = x_interp; 151 gate.y_interp = y_interp; 152 gate.radius = radius; 153 gate.mindata = mindata; 154 gate.maxdata = maxdata; 155 gate.variogram = variogram; 156 gate.observations = observations; 157 gate.predictions = predictions; 158 gate.error = error; 159 gate.numdone = xNewZeroInit<int>(num); 133 160 134 161 /*launch the thread manager with Krigingxt as a core: */ 135 162 LaunchThread(Krigingxt,(void*)&gate,num); 136 163 _printLine_("\r interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%"); 137 xDelete< double>(gate.percent);164 xDelete<int>(gate.numdone); 138 165 } 139 166 else{ … … 176 203 double *predictions = gate->predictions; 177 204 double *error = gate->error; 178 double *percent = gate->percent; 179 180 /*Intermediaries*/ 181 double localpercent; 205 int *numdone = gate->numdone; 182 206 183 207 /*partition loop across threads: */ … … 186 210 187 211 /*Print info*/ 188 percent[my_thread]=double(idx-i0)/double(i1-i0)*100.; 189 localpercent=percent[0]; 190 for(int i=1;i<num_threads;i++) localpercent=min(localpercent,percent[i]); 191 if(my_thread==0) _printString_("\r interpolation progress: "<<setw(6)<<setprecision(2)<<localpercent<<"%"); 212 numdone[my_thread]=idx-i0; 213 if(my_thread==0){ 214 int alldone=numdone[0]; 215 for(int i=1;i<num_threads;i++) alldone+=numdone[i]; 216 _printString_("\r interpolation progress: "<<setw(6)<<setprecision(2)<<double(alldone)/double(n_interp)*100.<<"%"); 217 } 192 218 193 219 /*Kriging interpolation*/ … … 224 250 double *predictions = gate->predictions; 225 251 double *error = gate->error; 226 double *percent = gate->percent; 227 228 /*Intermediaries*/ 229 int i; 230 double localpercent; 252 int *numdone = gate->numdone; 231 253 232 254 /*partition loop across threads: */ … … 235 257 236 258 /*Print info*/ 237 percent[my_thread]=double(idx-i0)/double(i1-i0)*100.; 238 localpercent=percent[0]; 239 for(i=1;i<num_threads;i++) localpercent=min(localpercent,percent[i]); 240 if(my_thread==0) _printString_("\r interpolation progress: "<<setw(6)<<setprecision(2)<<localpercent<<"%"); 259 numdone[my_thread]=idx-i0; 260 if(my_thread==0){ 261 int alldone=numdone[0]; 262 for(int i=1;i<num_threads;i++) alldone+=numdone[i]; 263 _printString_("\r interpolation progress: "<<setw(6)<<setprecision(2)<<double(alldone)/double(n_interp)*100.<<"%"); 264 } 241 265 242 266 observations->InterpolationNearestNeighbor(&predictions[idx],x_interp[idx],y_interp[idx],radius); … … 272 296 double *predictions = gate->predictions; 273 297 double *error = gate->error; 274 double *percent = gate->percent;298 int *numdone = gate->numdone; 275 299 double power = gate->power; 276 277 /*Intermediaries*/278 double localpercent;279 300 280 301 /*partition loop across threads: */ … … 283 304 284 305 /*Print info*/ 285 percent[my_thread]=double(idx-i0)/double(i1-i0)*100.; 286 localpercent=percent[0]; 287 for(int i=1;i<num_threads;i++) localpercent=min(localpercent,percent[i]); 288 if(my_thread==0) _printString_("\r interpolation progress: "<<setw(6)<<setprecision(2)<<localpercent<<"%"); 306 numdone[my_thread]=idx-i0; 307 if(my_thread==0){ 308 int alldone=numdone[0]; 309 for(int i=1;i<num_threads;i++) alldone+=numdone[i]; 310 _printString_("\r interpolation progress: "<<setw(6)<<setprecision(2)<<double(alldone)/double(n_interp)*100.<<"%"); 311 } 289 312 290 313 observations->InterpolationIDW(&predictions[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata,power); 314 } 315 return NULL; 316 }/*}}}*/ 317 /*FUNCTION v4t{{{*/ 318 void* v4t(void* vpthread_handle){ 319 320 /*gate variables :*/ 321 KrigingxThreadStruct *gate = NULL; 322 pthread_handle *handle = NULL; 323 int my_thread; 324 int num_threads; 325 int i0,i1; 326 327 /*recover handle and gate: */ 328 handle = (pthread_handle*)vpthread_handle; 329 gate = (KrigingxThreadStruct*)handle->gate; 330 my_thread = handle->id; 331 num_threads = handle->num; 332 333 /*recover parameters :*/ 334 int n_interp = gate->n_interp; 335 double *x_interp = gate->x_interp; 336 double *y_interp = gate->y_interp; 337 double radius = gate->radius; 338 int mindata = gate->mindata; 339 int maxdata = gate->maxdata; 340 Variogram *variogram = gate->variogram; 341 Observations *observations = gate->observations; 342 double *predictions = gate->predictions; 343 double *error = gate->error; 344 int *numdone = gate->numdone; 345 346 /*partition loop across threads: */ 347 PartitionRange(&i0,&i1,n_interp,num_threads,my_thread); 348 for(int idx=i0;idx<i1;idx++){ 349 350 /*Print info*/ 351 numdone[my_thread]=idx-i0; 352 if(my_thread==0){ 353 int alldone=numdone[0]; 354 for(int i=1;i<num_threads;i++) alldone+=numdone[i]; 355 _printString_("\r interpolation progress: "<<setw(6)<<setprecision(2)<<double(alldone)/double(n_interp)*100.<<"%"); 356 } 357 358 observations->InterpolationV4(&predictions[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata); 291 359 } 292 360 return NULL; -
issm/trunk/src/c/modules/Krigingx/Krigingx.h
r13395 r14067 28 28 double *predictions; 29 29 double *error; 30 double *percent;30 int *numdone; 31 31 double power;//for idw 32 32 }KrigingxThreadStruct; … … 35 35 void* NearestNeighbort(void*); 36 36 void* idwt(void*); 37 void* v4t(void*); 37 38 #endif /* _KRIGINGX_H */ -
issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.cpp
r13975 r14067 18 18 #include "../../EnumDefinitions/EnumDefinitions.h" 19 19 20 void NodeConnectivityx( double** pconnectivity, int* pwidth, double* elements, int nel, int nods){20 void NodeConnectivityx(int** pconnectivity,int* pwidth,int* elements, int nels, int nods){ 21 21 22 22 int i,j,n; … … 29 29 int num_elements; 30 30 int already_plugged=0; 31 double element; 32 33 /*output: */ 34 double* connectivity=NULL; 31 int element; 35 32 36 33 /*Allocate connectivity: */ 37 connectivity=xNewZeroInit<double>(nods*width);34 int* connectivity=xNewZeroInit<int>(nods*width); 38 35 39 36 /*Go through all elements, and for each elements, plug into the connectivity, all the nodes. 40 37 * If nodes are already plugged into the connectivity, skip them.: */ 41 for(n=0;n<nel ;n++){38 for(n=0;n<nels;n++){ 42 39 43 element= (double)(n+1); //matlab indexing40 element=n+1; //matlab indexing 44 41 45 42 for(i=0;i<3;i++){ 46 43 47 node= (int)*(elements+n*3+i); //already matlab indexed, elements comes directly from the workspace.44 node=elements[n*3+i]; //already matlab indexed, elements comes directly from the workspace. 48 45 index=node-1; 49 46 50 num_elements= (int)*(connectivity+width*index+maxels); //retrieve number of elements already plugged into the connectivity of this node.47 num_elements=connectivity[width*index+maxels]; //retrieve number of elements already plugged into the connectivity of this node. 51 48 52 49 already_plugged=0; … … 60 57 61 58 /*this elements is not yet plugged into the connectivity for this node, do it, and increase counter: */ 62 *(connectivity+width*index+num_elements)=element;63 *(connectivity+width*index+maxels)=(double)(num_elements+1);59 connectivity[width*index+num_elements]=element; 60 connectivity[width*index+maxels]=num_elements+1; 64 61 65 62 } … … 69 66 * warn the user to increase the connectivity width: */ 70 67 for(i=0;i<nods;i++){ 71 if (*(connectivity+width*i+maxels)>maxels)_error_("max connectivity width reached (" << *(connectivity+width*i+maxels) << ")! increase width of connectivity table"); 68 if (*(connectivity+width*i+maxels)>maxels) 69 _error_("max connectivity width reached (" << *(connectivity+width*i+maxels) << ")! increase width of connectivity table"); 72 70 } 73 71 -
issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.h
r13975 r14067 7 7 8 8 /* local prototypes: */ 9 void NodeConnectivityx( double** pconnectivity, int* pwidth,double* elements, int nel, int nods);9 void NodeConnectivityx(int** pconnectivity,int* pwidth,int* elements,int nels, int nods); 10 10 11 11 #endif /* _NODECONNECTIVITYX_H */ -
issm/trunk/src/c/shared/Threads/LaunchThread.cpp
r13975 r14067 31 31 pthread_t *threads = NULL; 32 32 pthread_handle *handles = NULL; 33 34 /*check number of threads*/ 35 if(num_threads<1) _error_("number of threads must be at least 1"); 36 if(num_threads>2000) _error_("number of threads seems to be too high ("<<num_threads<<">2000)"); 33 37 34 38 /*dynamically allocate: */ -
issm/trunk/src/dox/issm.dox
r13975 r14067 46 46 </th> 47 47 <tr> 48 <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>48 <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> 49 49 </tr> 50 50 <tr> 51 <th bgcolor=#C6E2FF style="text-align:left;"> MATLAB </th><td bgcolor=#C6E2FF style="text-align:right;">9 13</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>51 <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> 52 52 </tr> 53 53 <tr> 54 <th bgcolor=#FFFFFF style="text-align:left;"> C/C++ Header </th><td bgcolor=#FFFFFF style="text-align:right;">3 88</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>54 <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> 55 55 </tr> 56 56 <tr> 57 <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>57 <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> 58 58 </tr> 59 59 <tr> 60 <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>60 <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> 61 61 </tr> 62 62 <tr> 63 <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;">9 8</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>63 <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> 64 64 </tr> 65 65 <tr> 66 <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;">5 8</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>66 <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> 67 67 </tr> 68 68 <tr> … … 76 76 </tr> 77 77 <tr> 78 <th bgcolor=#FFFFFF style="text-align:left;"> SUM: </th><td bgcolor=#FFFFFF style="text-align:right;">19 12</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>78 <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> 79 79 </tr> 80 80 </table> -
issm/trunk/src/m/boundaryconditions/SetIceShelfBC.py
r13975 r14067 26 26 raise IOError("SetIceShelfBC error message: ice front file '%s' not found." % icefrontfile) 27 27 [nodeinsideicefront,dum]=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),icefrontfile,'node',2) 28 nodeonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1)) .astype(float)28 nodeonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1)) 29 29 else: 30 nodeonicefront=numpy.zeros((md.mesh.numberofvertices) )30 nodeonicefront=numpy.zeros((md.mesh.numberofvertices),bool) 31 31 32 32 # pos=find(md.mesh.vertexonboundary & ~nodeonicefront); … … 56 56 #segment on Neumann (Ice Front) 57 57 # pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2))); 58 pos=numpy.nonzero(numpy.logical_or(nodeonicefront[md.mesh.segments[:,0] .astype(int)-1],nodeonicefront[md.mesh.segments[:,1].astype(int)-1]))[0]58 pos=numpy.nonzero(numpy.logical_or(nodeonicefront[md.mesh.segments[:,0]-1],nodeonicefront[md.mesh.segments[:,1]-1]))[0] 59 59 if md.mesh.dimension==2: 60 60 pressureload=md.mesh.segments[pos,:] … … 62 62 # 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)]; 63 63 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])) 64 pressureload=numpy.zeros((0,5) )64 pressureload=numpy.zeros((0,5),int) 65 65 for i in xrange(1,md.mesh.numberoflayers): 66 66 # pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ]; … … 69 69 #Add water or air enum depending on the element 70 70 # pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end))]; 71 pressureload=numpy.hstack((pressureload,1 .*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape(-1,1)))71 pressureload=numpy.hstack((pressureload,1*md.mask.elementonfloatingice[pressureload[:,-1]-1].reshape(-1,1))) 72 72 73 73 #plug onto model -
issm/trunk/src/m/boundaryconditions/SetMarineIceSheetBC.py
r13975 r14067 28 28 raise IOError("SetMarineIceSheetBC error message: ice front file '%s' not found." % icefrontfile) 29 29 [nodeinsideicefront,dum]=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),icefrontfile,'node',2) 30 vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1)) .astype(float)30 vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1)) 31 31 else: 32 32 #Guess where the ice front is 33 vertexonfloatingice=numpy.zeros((md.mesh.numberofvertices) )34 vertexonfloatingice[md.mesh.elements[numpy.nonzero(md.mask.elementonfloatingice),:] .astype(int)-1]=135 vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,vertexonfloatingice) .astype(float)33 vertexonfloatingice=numpy.zeros((md.mesh.numberofvertices),bool) 34 vertexonfloatingice[md.mesh.elements[numpy.nonzero(md.mask.elementonfloatingice),:]-1]=True 35 vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,vertexonfloatingice) 36 36 37 37 # pos=find(md.mesh.vertexonboundary & ~vertexonicefront); … … 62 62 #segment on Neumann (Ice Front) 63 63 # pos=find(vertexonicefront(md.mesh.segments(:,1)) | vertexonicefront(md.mesh.segments(:,2))); 64 pos=numpy.nonzero(numpy.logical_or(vertexonicefront[md.mesh.segments[:,0] .astype(int)-1],vertexonicefront[md.mesh.segments[:,1].astype(int)-1]))[0]64 pos=numpy.nonzero(numpy.logical_or(vertexonicefront[md.mesh.segments[:,0]-1],vertexonicefront[md.mesh.segments[:,1]-1]))[0] 65 65 if md.mesh.dimension==2: 66 66 pressureload=md.mesh.segments[pos,:] … … 68 68 # 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)]; 69 69 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])) 70 pressureload=numpy.zeros((0,5) )70 pressureload=numpy.zeros((0,5),int) 71 71 for i in xrange(1,md.mesh.numberoflayers): 72 72 # pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ]; … … 75 75 #Add water or air enum depending on the element 76 76 # pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end))+ 0*md.mask.elementongroundedice(pressureload(:,end))]; 77 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)))77 pressureload=numpy.hstack((pressureload,1*md.mask.elementonfloatingice[pressureload[:,-1]-1].reshape(-1,1)+0*md.mask.elementongroundedice[pressureload[:,-1]-1].reshape(-1,1))) 78 78 79 79 #plug onto model -
issm/trunk/src/m/classes/clusters/greenplanet.m
r13395 r14067 14 14 cpuspernode=8; 15 15 port=8000; 16 queue=' rignot';16 queue='c6145'; 17 17 codepath=''; 18 18 executionpath=''; … … 50 50 function md = checkconsistency(cluster,md,solution,analyses) % {{{ 51 51 52 available_queues={' rignot','default'};52 available_queues={'c6145','default'}; 53 53 queue_requirements_time=[Inf Inf]; 54 54 queue_requirements_np=[80 80]; -
issm/trunk/src/m/classes/flowequation.py
r13975 r14067 19 19 # {{{ Properties 20 20 21 self.ismacayealpattyn = 0 ;22 self.ishutter = 0 ;23 self.isl1l2 = 0;24 self.isstokes = 0 ;21 self.ismacayealpattyn = 0 22 self.ishutter = 0 23 self.isl1l2 = 0 24 self.isstokes = 0 25 25 self.vertex_equation = float('NaN') 26 26 self.element_equation = float('NaN') … … 37 37 string=' flow equation parameters:' 38 38 39 string="%s\n\n%s"%(string,fielddisplay(self,'ismacayealpattyn', 'is the macayeal or pattyn approximation used ?'))40 string="%s\n%s"%(string,fielddisplay(self,'ishutter', 'is the shallow ice approximation used ?'))41 string="%s\n%s"%(string,fielddisplay(self,'isl1l2', 'are l1l2 equations used ?'))42 string="%s\n%s"%(string,fielddisplay(self,'isstokes', 'are the Full-Stokes equations used ?'))43 string="%s\n%s"%(string,fielddisplay(self,'vertex_equation', 'flow equation for each vertex'))44 string="%s\n%s"%(string,fielddisplay(self,'element_equation', 'flow equation for each element'))45 string="%s\n%s"%(string,fielddisplay(self,'bordermacayeal', 'vertices on MacAyeal''s border (for tiling)'))46 string="%s\n%s"%(string,fielddisplay(self,'borderpattyn', 'vertices on Pattyn''s border (for tiling)'))47 string="%s\n%s"%(string,fielddisplay(self,'borderstokes', 'vertices on Stokes'' border (for tiling)'))39 string="%s\n\n%s"%(string,fielddisplay(self,'ismacayealpattyn',"is the macayeal or pattyn approximation used ?")) 40 string="%s\n%s"%(string,fielddisplay(self,'ishutter',"is the shallow ice approximation used ?")) 41 string="%s\n%s"%(string,fielddisplay(self,'isl1l2',"are l1l2 equations used ?")) 42 string="%s\n%s"%(string,fielddisplay(self,'isstokes',"are the Full-Stokes equations used ?")) 43 string="%s\n%s"%(string,fielddisplay(self,'vertex_equation',"flow equation for each vertex")) 44 string="%s\n%s"%(string,fielddisplay(self,'element_equation',"flow equation for each element")) 45 string="%s\n%s"%(string,fielddisplay(self,'bordermacayeal',"vertices on MacAyeal's border (for tiling)")) 46 string="%s\n%s"%(string,fielddisplay(self,'borderpattyn',"vertices on Pattyn's border (for tiling)")) 47 string="%s\n%s"%(string,fielddisplay(self,'borderstokes',"vertices on Stokes' border (for tiling)")) 48 48 return string 49 49 #}}} -
issm/trunk/src/m/classes/initialization.py
r13975 r14067 70 70 md = checkfield(md,'initialization.vy','NaN',1,'size',[md.mesh.numberofvertices]) 71 71 #Triangle with zero velocity 72 if numpy.any(numpy.logical_and(numpy.sum(numpy.abs(md.initialization.vx[md.mesh.elements .astype(int)-1]),axis=1)==0,\73 numpy.sum(numpy.abs(md.initialization.vy[md.mesh.elements .astype(int)-1]),axis=1)==0)):72 if numpy.any(numpy.logical_and(numpy.sum(numpy.abs(md.initialization.vx[md.mesh.elements-1]),axis=1)==0,\ 73 numpy.sum(numpy.abs(md.initialization.vy[md.mesh.elements-1]),axis=1)==0)): 74 74 md.checkmessage("at least one triangle has all its vertices with a zero velocity") 75 75 if ThermalAnalysisEnum() in analyses: -
issm/trunk/src/m/classes/mask.py
r13395 r14067 33 33 string="%s\n%s"%(string,fielddisplay(self,"elementonfloatingice","element on floating ice flags list")) 34 34 string="%s\n%s"%(string,fielddisplay(self,"vertexonfloatingice","vertex on floating ice flags list")) 35 string="%s\n%s"%(string,fielddisplay(self,"elementongroundedice","element on grounded ice 35 string="%s\n%s"%(string,fielddisplay(self,"elementongroundedice","element on grounded ice list")) 36 36 string="%s\n%s"%(string,fielddisplay(self,"vertexongroundedice","vertex on grounded ice flags list")) 37 37 string="%s\n%s"%(string,fielddisplay(self,"elementonwater","element on water flags list")) -
issm/trunk/src/m/classes/mesh.py
r13975 r14067 92 92 string="%s\n%s"%(string,fielddisplay(self,"vertexonsurface","upper vertices flags list")) 93 93 string="%s\n%s"%(string,fielddisplay(self,"elementonsurface","upper elements flags list")) 94 string="%s\n%s"%(string,fielddisplay(self,"uppervertex","upper vertex list ( NaNfor vertex on the upper surface)"))95 string="%s\n%s"%(string,fielddisplay(self,"upperelements","upper element list ( NaNfor element on the upper layer)"))96 string="%s\n%s"%(string,fielddisplay(self,"lowervertex","lower vertex list ( NaNfor vertex on the lower surface)"))97 string="%s\n%s"%(string,fielddisplay(self,"lowerelements","lower element list ( NaN for element on the lower layer"))94 string="%s\n%s"%(string,fielddisplay(self,"uppervertex","upper vertex list (-1 for vertex on the upper surface)")) 95 string="%s\n%s"%(string,fielddisplay(self,"upperelements","upper element list (-1 for element on the upper layer)")) 96 string="%s\n%s"%(string,fielddisplay(self,"lowervertex","lower vertex list (-1 for vertex on the lower surface)")) 97 string="%s\n%s"%(string,fielddisplay(self,"lowerelements","lower element list (-1 for element on the lower layer)")) 98 98 string="%s\n%s"%(string,fielddisplay(self,"vertexonboundary","vertices on the boundary of the domain flag list")) 99 99 string="%s\n%s"%(string,fielddisplay(self,"segments","edges on domain boundary (vertex1 vertex2 element)")) -
issm/trunk/src/m/classes/model/model.py
r13975 r14067 2 2 import numpy 3 3 import copy 4 import sys 4 5 from mesh import mesh 5 6 from mask import mask … … 207 208 #kick out all elements with 3 dirichlets 208 209 spc_elem=numpy.nonzero(numpy.logical_not(flag_elem))[0] 209 spc_node=numpy.unique(md1.mesh.elements[spc_elem,:]) .astype(int)-1210 spc_node=numpy.unique(md1.mesh.elements[spc_elem,:])-1 210 211 flag=numpy.ones(md1.mesh.numberofvertices) 211 212 flag[spc_node]=0 212 pos=numpy.nonzero(numpy.logical_not(numpy.sum(flag[md1.mesh.elements .astype(int)-1],axis=1)))[0]213 pos=numpy.nonzero(numpy.logical_not(numpy.sum(flag[md1.mesh.elements-1],axis=1)))[0] 213 214 flag_elem[pos]=0 214 215 215 216 #extracted elements and nodes lists 216 217 pos_elem=numpy.nonzero(flag_elem)[0] 217 pos_node=numpy.unique(md1.mesh.elements[pos_elem,:]) .astype(int)-1218 pos_node=numpy.unique(md1.mesh.elements[pos_elem,:])-1 218 219 219 220 #keep track of some fields … … 226 227 227 228 #Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements) 228 Pelem=numpy.zeros(numberofelements1 )229 Pelem=numpy.zeros(numberofelements1,int) 229 230 Pelem[pos_elem]=numpy.arange(1,numberofelements2+1) 230 Pnode=numpy.zeros(numberofvertices1 )231 Pnode=numpy.zeros(numberofvertices1,int) 231 232 Pnode[pos_node]=numpy.arange(1,numberofvertices2+1) 232 233 … … 234 235 elements_1=copy.deepcopy(md1.mesh.elements) 235 236 elements_2=elements_1[pos_elem,:] 236 elements_2[:,0]=Pnode[elements_2[:,0] .astype(int)-1]237 elements_2[:,1]=Pnode[elements_2[:,1] .astype(int)-1]238 elements_2[:,2]=Pnode[elements_2[:,2] .astype(int)-1]237 elements_2[:,0]=Pnode[elements_2[:,0]-1] 238 elements_2[:,1]=Pnode[elements_2[:,1]-1] 239 elements_2[:,2]=Pnode[elements_2[:,2]-1] 239 240 if md1.mesh.dimension==3: 240 elements_2[:,3]=Pnode[elements_2[:,3] .astype(int)-1]241 elements_2[:,4]=Pnode[elements_2[:,4] .astype(int)-1]242 elements_2[:,5]=Pnode[elements_2[:,5] .astype(int)-1]241 elements_2[:,3]=Pnode[elements_2[:,3]-1] 242 elements_2[:,4]=Pnode[elements_2[:,4]-1] 243 elements_2[:,5]=Pnode[elements_2[:,5]-1] 243 244 244 245 #OK, now create the new model! … … 291 292 if md1.mesh.dimension==3: 292 293 md2.mesh.uppervertex=md1.mesh.uppervertex[pos_node] 293 pos=numpy.nonzero(numpy.logical_not( numpy.isnan(md2.mesh.uppervertex)))[0]294 md2.mesh.uppervertex[pos]=Pnode[md2.mesh.uppervertex[pos] .astype(int)-1]294 pos=numpy.nonzero(numpy.logical_not(md2.mesh.uppervertex==-1))[0] 295 md2.mesh.uppervertex[pos]=Pnode[md2.mesh.uppervertex[pos]-1] 295 296 296 297 md2.mesh.lowervertex=md1.mesh.lowervertex[pos_node] 297 pos=numpy.nonzero(numpy.logical_not( numpy.isnan(md2.mesh.lowervertex)))[0]298 md2.mesh.lowervertex[pos]=Pnode[md2.mesh.lowervertex[pos] .astype(int)-1]298 pos=numpy.nonzero(numpy.logical_not(md2.mesh.lowervertex==-1))[0] 299 md2.mesh.lowervertex[pos]=Pnode[md2.mesh.lowervertex[pos]-1] 299 300 300 301 md2.mesh.upperelements=md1.mesh.upperelements[pos_elem] 301 pos=numpy.nonzero(numpy.logical_not( numpy.isnan(md2.mesh.upperelements)))[0]302 md2.mesh.upperelements[pos]=Pelem[md2.mesh.upperelements[pos] .astype(int)-1]302 pos=numpy.nonzero(numpy.logical_not(md2.mesh.upperelements==-1))[0] 303 md2.mesh.upperelements[pos]=Pelem[md2.mesh.upperelements[pos]-1] 303 304 304 305 md2.mesh.lowerelements=md1.mesh.lowerelements[pos_elem] 305 pos=numpy.nonzero(numpy.logical_not( numpy.isnan(md2.mesh.lowerelements)))[0]306 md2.mesh.lowerelements[pos]=Pelem[md2.mesh.lowerelements[pos] .astype(int)-1]306 pos=numpy.nonzero(numpy.logical_not(md2.mesh.lowerelements==-1))[0] 307 md2.mesh.lowerelements[pos]=Pelem[md2.mesh.lowerelements[pos]-1] 307 308 308 309 #Initial 2d mesh … … 316 317 md2.mesh.numberofvertices2d=numpy.size(pos_node_2d) 317 318 md2.mesh.elements2d=md1.mesh.elements2d[pos_elem_2d,:] 318 md2.mesh.elements2d[:,0]=Pnode[md2.mesh.elements2d[:,0] .astype(int)-1]319 md2.mesh.elements2d[:,1]=Pnode[md2.mesh.elements2d[:,1] .astype(int)-1]320 md2.mesh.elements2d[:,2]=Pnode[md2.mesh.elements2d[:,2] .astype(int)-1]319 md2.mesh.elements2d[:,0]=Pnode[md2.mesh.elements2d[:,0]-1] 320 md2.mesh.elements2d[:,1]=Pnode[md2.mesh.elements2d[:,1]-1] 321 md2.mesh.elements2d[:,2]=Pnode[md2.mesh.elements2d[:,2]-1] 321 322 322 323 md2.mesh.x2d=md1.mesh.x[pos_node_2d] … … 324 325 325 326 #Edges 326 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...327 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... 327 328 #renumber first two columns 328 329 pos=numpy.nonzero(md2.mesh.edges[:,3]!=-1)[0] 329 md2.mesh.edges[: ,0]=Pnode[md2.mesh.edges[:,0] .astype(int)-1]330 md2.mesh.edges[: ,1]=Pnode[md2.mesh.edges[:,1] .astype(int)-1]331 md2.mesh.edges[: ,2]=Pelem[md2.mesh.edges[:,2] .astype(int)-1]332 md2.mesh.edges[pos,3]=Pelem[md2.mesh.edges[pos,3] .astype(int)-1]330 md2.mesh.edges[: ,0]=Pnode[md2.mesh.edges[:,0]-1] 331 md2.mesh.edges[: ,1]=Pnode[md2.mesh.edges[:,1]-1] 332 md2.mesh.edges[: ,2]=Pelem[md2.mesh.edges[:,2]-1] 333 md2.mesh.edges[pos,3]=Pelem[md2.mesh.edges[pos,3]-1] 333 334 #remove edges when the 2 vertices are not in the domain. 334 335 md2.mesh.edges=md2.mesh.edges[numpy.nonzero(numpy.logical_and(md2.mesh.edges[:,0],md2.mesh.edges[:,1]))[0],:] … … 364 365 [md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity) 365 366 md2.mesh.segments=contourenvelope(md2) 366 md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2 )367 md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2] .astype(int)-1]=1367 md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2,bool) 368 md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2]-1]=True 368 369 else: 369 370 #First do the connectivity for the contourenvelope in 2d … … 371 372 [md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity) 372 373 md2.mesh.segments=contourenvelope(md2) 373 md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2/md2.mesh.numberoflayers )374 md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2] .astype(int)-1]=1374 md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2/md2.mesh.numberoflayers,bool) 375 md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2]-1]=True 375 376 md2.mesh.vertexonboundary=numpy.tile(md2.mesh.vertexonboundary,md2.mesh.numberoflayers) 376 377 #Then do it for 3d as usual … … 381 382 #Catch the elements that have not been extracted 382 383 orphans_elem=numpy.nonzero(numpy.logical_not(flag_elem))[0] 383 orphans_node=numpy.unique(md1.mesh.elements[orphans_elem,:]) .astype(int)-1384 orphans_node=numpy.unique(md1.mesh.elements[orphans_elem,:])-1 384 385 #Figure out which node are on the boundary between md2 and md1 385 386 nodestoflag1=numpy.intersect1d(orphans_node,pos_node) … … 443 444 444 445 #Keep track of pos_node and pos_elem 445 md2.mesh.extractedvertices=pos_node .astype(float)+1446 md2.mesh.extractedelements=pos_elem .astype(float)+1446 md2.mesh.extractedvertices=pos_node+1 447 md2.mesh.extractedelements=pos_elem+1 447 448 448 449 return md2 … … 526 527 527 528 #Extrude elements 528 elements3d=numpy.empty((0,6) )529 elements3d=numpy.empty((0,6),int) 529 530 for i in xrange(numlayers-1): 530 531 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 … … 532 533 533 534 #Keep a trace of lower and upper nodes 534 mesh.lowervertex= float('NaN')*numpy.ones(number_nodes3d)535 mesh.uppervertex= float('NaN')*numpy.ones(number_nodes3d)535 mesh.lowervertex=-1*numpy.ones(number_nodes3d,int) 536 mesh.uppervertex=-1*numpy.ones(number_nodes3d,int) 536 537 mesh.lowervertex[md.mesh.numberofvertices:]=numpy.arange(1,(numlayers-1)*md.mesh.numberofvertices+1) 537 538 mesh.uppervertex[:(numlayers-1)*md.mesh.numberofvertices]=numpy.arange(md.mesh.numberofvertices+1,number_nodes3d+1) … … 540 541 541 542 #same for lower and upper elements 542 mesh.lowerelements= float('NaN')*numpy.ones(number_el3d)543 mesh.upperelements= float('NaN')*numpy.ones(number_el3d)543 mesh.lowerelements=-1*numpy.ones(number_el3d,int) 544 mesh.upperelements=-1*numpy.ones(number_el3d,int) 544 545 mesh.lowerelements[md.mesh.numberofelements:]=numpy.arange(1,(numlayers-2)*md.mesh.numberofelements+1) 545 546 mesh.upperelements[:(numlayers-2)*md.mesh.numberofelements]=numpy.arange(md.mesh.numberofelements+1,(numlayers-1)*md.mesh.numberofelements+1) … … 603 604 604 605 #bedinfo and surface info 605 md.mesh.elementonbed=project3d(md,'vector',numpy.ones(md.mesh.numberofelements2d ),'type','element','layer',1)606 md.mesh.elementonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofelements2d ),'type','element','layer',md.mesh.numberoflayers-1)607 md.mesh.vertexonbed=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d ),'type','node','layer',1)608 md.mesh.vertexonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d ),'type','node','layer',md.mesh.numberoflayers)606 md.mesh.elementonbed=project3d(md,'vector',numpy.ones(md.mesh.numberofelements2d,bool),'type','element','layer',1) 607 md.mesh.elementonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofelements2d,bool),'type','element','layer',md.mesh.numberoflayers-1) 608 md.mesh.vertexonbed=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',1) 609 md.mesh.vertexonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',md.mesh.numberoflayers) 609 610 610 611 #elementstype 611 612 if not numpy.any(numpy.isnan(md.flowequation.element_equation)): 612 613 oldelements_type=md.flowequation.element_equation 613 md.flowequation.element_equation=numpy.zeros(number_el3d )614 md.flowequation.element_equation=numpy.zeros(number_el3d,int) 614 615 md.flowequation.element_equation=project3d(md,'vector',oldelements_type,'type','element') 615 616 … … 617 618 if not numpy.any(numpy.isnan(md.flowequation.vertex_equation)): 618 619 oldvertices_type=md.flowequation.vertex_equation 619 md.flowequation.vertex_equation=numpy.zeros(number_nodes3d )620 md.flowequation.vertex_equation=numpy.zeros(number_nodes3d,int) 620 621 md.flowequation.vertex_equation=project3d(md,'vector',oldvertices_type,'type','node') 621 622 … … 635 636 #in 3d, pressureload: [node1 node2 node3 node4 element] 636 637 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 637 pressureload=numpy.empty((0,6) )638 pressureload=numpy.empty((0,6),int) 638 639 for i in xrange(numlayers-1): 639 640 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])))) … … 642 643 #connectivity 643 644 md.mesh.elementconnectivity=numpy.tile(md.mesh.elementconnectivity,(numlayers-1,1)) 644 md.mesh.elementconnectivity[numpy.nonzero(md.mesh.elementconnectivity==0)]= float('NaN')645 md.mesh.elementconnectivity[numpy.nonzero(md.mesh.elementconnectivity==0)]=-sys.maxint-1 645 646 for i in xrange(1,numlayers-1): 646 647 md.mesh.elementconnectivity[i*md.mesh.numberofelements2d:(i+1)*md.mesh.numberofelements2d,:] \ 647 648 =md.mesh.elementconnectivity[i*md.mesh.numberofelements2d:(i+1)*md.mesh.numberofelements2d,:]+md.mesh.numberofelements2d 648 md.mesh.elementconnectivity[numpy.nonzero( numpy.isnan(md.mesh.elementconnectivity))]=0649 md.mesh.elementconnectivity[numpy.nonzero(md.mesh.elementconnectivity<0)]=0 649 650 650 651 #materials -
issm/trunk/src/m/classes/transient.m
r13395 r14067 46 46 47 47 fielddisplay(obj,'isprognostic','indicates if a prognostic solution is used in the transient'); 48 fielddisplay(obj,'isdiagnostic','indicates if a diagnostic solution is used in the transient'); 48 49 fielddisplay(obj,'isthermal','indicates if a thermal solution is used in the transient'); 49 fielddisplay(obj,'isdiagnostic','indicates if a diagnostic solution is used in the transient');50 50 fielddisplay(obj,'isgroundingline','indicates if a groundingline migration is used in the transient'); 51 51 fielddisplay(obj,'requested_outputs','list of additional outputs requested'); -
issm/trunk/src/m/classes/transient.py
r13395 r14067 16 16 def __init__(self): 17 17 # {{{ Properties 18 self.isprognostic = 019 self.isdiagnostic = 020 self.isthermal = 021 self.isgroundingline = 018 self.isprognostic = False 19 self.isdiagnostic = False 20 self.isthermal = False 21 self.isgroundingline = False 22 22 self.requested_outputs = float('NaN') 23 23 … … 30 30 string=' transient solution parameters:' 31 31 string="%s\n%s"%(string,fielddisplay(self,'isprognostic','indicates if a prognostic solution is used in the transient')) 32 string="%s\n%s"%(string,fielddisplay(self,'isdiagnostic','indicates if a diagnostic solution is used in the transient')) 32 33 string="%s\n%s"%(string,fielddisplay(self,'isthermal','indicates if a thermal solution is used in the transient')) 33 string="%s\n%s"%(string,fielddisplay(self,'isdiagnostic','indicates if a diagnostic solution is used in the transient'))34 34 string="%s\n%s"%(string,fielddisplay(self,'isgroundingline','indicates if a groundingline migration is used in the transient')) 35 35 string="%s\n%s"%(string,fielddisplay(self,'requested_outputs','list of additional outputs requested')) … … 41 41 42 42 #full analysis: Diagnostic, Prognostic and Thermal but no groundingline migration for now 43 self.isprognostic= 144 self.isdiagnostic= 145 self.isthermal= 146 self.isgroundingline= 043 self.isprognostic=True 44 self.isdiagnostic=True 45 self.isthermal=True 46 self.isgroundingline=False 47 47 48 48 return self -
issm/trunk/src/m/consistency/checkfield.py
r13395 r14067 143 143 if options.getfieldvalue('forcing',0): 144 144 if numpy.size(field,0)==md.mesh.numberofvertices: 145 if len(numpy.shape(field))>1 and not numpy.size(field,1)==1:145 if numpy.ndim(field)>1 and not numpy.size(field,1)==1: 146 146 md = md.checkmessage(options.getfieldvalue('message',\ 147 147 "field '%s' should have only one column as there are md.mesh.numberofvertices lines" % fieldname)) -
issm/trunk/src/m/exp/expcoarsen.m
r13975 r14067 25 25 26 26 %Get exp oldfile 27 [path root ext ver]=fileparts(oldfile);27 [path root ext]=fileparts(oldfile); 28 28 A=expread(oldfile); 29 29 numprofiles=size(A,2); -
issm/trunk/src/m/extrusion/project3d.py
r13395 r14067 38 38 39 39 vector1d=False 40 if isinstance(vector2d,numpy.ndarray) and len(vector2d.shape)==1:40 if isinstance(vector2d,numpy.ndarray) and numpy.ndim(vector2d)==1: 41 41 vector1d=True 42 42 vector2d=vector2d.reshape(numpy.size(vector2d),1) … … 49 49 #Initialize 3d vector 50 50 if numpy.size(vector2d,axis=0)==md.mesh.numberofvertices2d: 51 projected_vector= paddingvalue*numpy.ones((md.mesh.numberofvertices, numpy.size(vector2d,axis=1)))51 projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofvertices, numpy.size(vector2d,axis=1)))).astype(vector2d.dtype) 52 52 elif numpy.size(vector2d,axis=0)==md.mesh.numberofvertices2d+1: 53 projected_vector= paddingvalue*numpy.ones((md.mesh.numberofvertices+1,numpy.size(vector2d,axis=1)))53 projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofvertices+1,numpy.size(vector2d,axis=1)))).astype(vector2d.dtype) 54 54 projected_vector[-1,:]=vector2d[-1,:] 55 55 vector2d=vector2d[:-1,:] … … 68 68 #Initialize 3d vector 69 69 if numpy.size(vector2d,axis=0)==md.mesh.numberofelements2d: 70 projected_vector= paddingvalue*numpy.ones((md.mesh.numberofelements, numpy.size(vector2d,axis=1)))70 projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofelements, numpy.size(vector2d,axis=1)))).astype(vector2d.dtype) 71 71 elif numpy.size(vector2d,axis=0)==md.mesh.numberofelements2d+1: 72 projected_vector= paddingvalue*numpy.ones((md.mesh.numberofelements+1,numpy.size(vector2d,axis=1)))72 projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofelements+1,numpy.size(vector2d,axis=1)))).astype(vector2d.dtype) 73 73 projected_vector[-1,:]=vector2d[-1,:] 74 74 vector2d=vector2d[:-1,:] -
issm/trunk/src/m/geometry/FlagElements.py
r13975 r14067 24 24 if isinstance(region,(str,unicode)): 25 25 if not region: 26 flag=numpy.zeros(md.mesh.numberofelements, 'bool')26 flag=numpy.zeros(md.mesh.numberofelements,bool) 27 27 invert=0 28 28 elif strcmpi(region,'all'): 29 flag=numpy.ones(md.mesh.numberofelements, 'bool')29 flag=numpy.ones(md.mesh.numberofelements,bool) 30 30 invert=0 31 31 else: … … 43 43 raise RuntimeError("FlagElements.py calling basinzoom.py is not complete.") 44 44 xlim,ylim=basinzoom('basin',region) 45 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)46 flag=numpy.prod(flag_nodes[md.mesh.elements],axis=1) 45 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])) 46 flag=numpy.prod(flag_nodes[md.mesh.elements],axis=1).astype(bool) 47 47 else: 48 48 #ok, flag elements 49 49 [flag,dum]=ContourToMesh(md.mesh.elements[:,0:3].copy(),md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),region,'element',1) 50 flag=flag.astype(bool) 50 51 51 52 if invert: -
issm/trunk/src/m/geometry/GetAreas.py
r13975 r14067 33 33 #initialization 34 34 areas=numpy.zeros(nels) 35 x1=x[index[:,0] .astype(int)-1]36 x2=x[index[:,1] .astype(int)-1]37 x3=x[index[:,2] .astype(int)-1]38 y1=y[index[:,0] .astype(int)-1]39 y2=y[index[:,1] .astype(int)-1]40 y3=y[index[:,2] .astype(int)-1]35 x1=x[index[:,0]-1] 36 x2=x[index[:,1]-1] 37 x3=x[index[:,2]-1] 38 y1=y[index[:,0]-1] 39 y2=y[index[:,1]-1] 40 y3=y[index[:,2]-1] 41 41 42 42 #compute the volume of each element … … 46 46 else: 47 47 #V=area(triangle)*1/3(z1+z2+z3) 48 thickness=numpy.mean(z[index[:,3:6] .astype(int)-1])-numpy.mean(z[index[:,0:3].astype(int)-1])48 thickness=numpy.mean(z[index[:,3:6]-1])-numpy.mean(z[index[:,0:3]-1]) 49 49 areas=(0.5*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)))*thickness 50 50 -
issm/trunk/src/m/mesh/ComputeHessian.py
r13975 r14067 45 45 46 46 #Compute gradient for each element 47 grad_elx=numpy.sum(field[index .astype(int)-1,0]*alpha,axis=1)48 grad_ely=numpy.sum(field[index .astype(int)-1,0]*beta,axis=1)47 grad_elx=numpy.sum(field[index-1,0]*alpha,axis=1) 48 grad_ely=numpy.sum(field[index-1,0]*beta,axis=1) 49 49 50 50 #Compute gradient for each node (average of the elements around) … … 55 55 56 56 #Compute hessian for each element 57 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)))57 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))) 58 58 59 59 if strcmpi(type,'node'): -
issm/trunk/src/m/mesh/GetNodalFunctionsCoeff.py
r13975 r14067 39 39 40 40 #compute nodal functions coefficients N(x,y)=alpha x + beta y +gamma 41 x1=x[index[:,0] .astype(int)-1]42 x2=x[index[:,1] .astype(int)-1]43 x3=x[index[:,2] .astype(int)-1]44 y1=y[index[:,0] .astype(int)-1]45 y2=y[index[:,1] .astype(int)-1]46 y3=y[index[:,2] .astype(int)-1]41 x1=x[index[:,0]-1] 42 x2=x[index[:,1]-1] 43 x3=x[index[:,2]-1] 44 y1=y[index[:,0]-1] 45 y2=y[index[:,1]-1] 46 y3=y[index[:,2]-1] 47 47 invdet=1./(x1*(y2-y3)-x2*(y1-y3)+x3*(y1-y2)) 48 48 -
issm/trunk/src/m/mesh/bamg.py
r13975 r14067 320 320 md.mesh.x=bamgmesh_out['Vertices'][:,0].copy() 321 321 md.mesh.y=bamgmesh_out['Vertices'][:,1].copy() 322 md.mesh.elements=bamgmesh_out['Triangles'][:,0:3]. copy()323 md.mesh.edges=bamgmesh_out['IssmEdges']. copy()324 md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3]. copy()325 md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3]. copy()322 md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int) 323 md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int) 324 md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int) 325 md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int) 326 326 327 327 #Fill in rest of fields: … … 331 331 md.mesh.numberofedges=numpy.size(md.mesh.edges,axis=0) 332 332 md.mesh.z=numpy.zeros(md.mesh.numberofvertices) 333 md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices )334 md.mask.vertexonwater=numpy.zeros(md.mesh.numberofvertices )335 md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices )336 md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements )337 md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements )338 md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices )339 md.mesh.vertexonboundary[md.mesh.segments[:,0:2] .astype(int)-1]=1333 md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices,bool) 334 md.mask.vertexonwater=numpy.zeros(md.mesh.numberofvertices,bool) 335 md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices,bool) 336 md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements,bool) 337 md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements,bool) 338 md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices,bool) 339 md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True 340 340 md.mesh.elementconnectivity=md.private.bamg['mesh'].ElementConnectivity 341 341 md.mesh.elementconnectivity[numpy.nonzero(numpy.isnan(md.mesh.elementconnectivity))]=0 342 md.mesh.elementconnectivity=md.mesh.elementconnectivity.astype(int) 342 343 343 344 #Check for orphan -
issm/trunk/src/m/mesh/meshconvert.py
r13975 r14067 35 35 md.mesh.x=bamgmesh_out['Vertices'][:,0].copy() 36 36 md.mesh.y=bamgmesh_out['Vertices'][:,1].copy() 37 md.mesh.elements=bamgmesh_out['Triangles'][:,0:3]. copy()38 md.mesh.edges=bamgmesh_out['IssmEdges']. copy()39 md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3]. copy()40 md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3]. copy()37 md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int) 38 md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int) 39 md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int) 40 md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int) 41 41 42 42 #Fill in rest of fields: … … 46 46 md.mesh.numberofedges=numpy.size(md.mesh.edges,axis=0) 47 47 md.mesh.z=numpy.zeros(md.mesh.numberofvertices) 48 md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices )49 md.mask.vertexonwater=numpy.zeros(md.mesh.numberofvertices )50 md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices )51 md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements )52 md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements )53 md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices )54 md.mesh.vertexonboundary[md.mesh.segments[:,0:2] .astype(int)-1]=148 md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices,bool) 49 md.mask.vertexonwater=numpy.zeros(md.mesh.numberofvertices,bool) 50 md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices,bool) 51 md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements,bool) 52 md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements,bool) 53 md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices,bool) 54 md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True 55 55 56 56 return md -
issm/trunk/src/m/mesh/rifts/meshprocessoutsiderifts.py
r13975 r14067 41 41 B=node_connected_to_tip 42 42 43 elements=numpy.empty(0 )43 elements=numpy.empty(0,int) 44 44 45 45 while flags(B): #as long as B does not belong to the domain outline, keep looking. … … 62 62 63 63 #replace tip in elements 64 newelements=md.mesh.elements[elements .astype(int)-1,:]64 newelements=md.mesh.elements[elements-1,:] 65 65 pos=numpy.nonzero(newelements==tip) 66 66 newelements[pos]=num 67 md.mesh.elements[elements .astype(int)-1,:]=newelements67 md.mesh.elements[elements-1,:]=newelements 68 68 rift.tips=numpy.concatenate((rift.tips,num)) 69 69 … … 81 81 md.mesh.numberofvertices=numpy.size(md.mesh.x) 82 82 md.mesh.z=numpy.zeros(md.mesh.numberofvertices) 83 md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x) )84 md.mesh.vertexonboundary[md.mesh.segments[:,0:2] .astype(int)-1]=183 md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x),bool) 84 md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True 85 85 md.rifts.numrifts=length(md.rifts.riftstruct) 86 md.flowequation.element_equation=3 .*numpy.ones(md.mesh.numberofelements)87 md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices )88 md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices )89 md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements )90 md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements )86 md.flowequation.element_equation=3*numpy.ones(md.mesh.numberofelements,int) 87 md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices,bool) 88 md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices,bool) 89 md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements,bool) 90 md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements,bool) 91 91 92 92 return md -
issm/trunk/src/m/mesh/rifts/meshprocessrifts.py
r13975 r14067 23 23 #Call MEX file 24 24 [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) 25 md.mesh.elements=md.mesh.elements.astype(int) 25 26 md.mesh.x=md.mesh.x.reshape(-1) 26 27 md.mesh.y=md.mesh.y.reshape(-1) 28 md.mesh.segments=md.mesh.segments.astype(int) 29 md.mesh.segmentmarkers=md.mesh.segmentmarkers.astype(int) 27 30 if not isinstance(md.rifts.riftstruct,list) or not md.rifts.riftstruct: 28 31 raise RuntimeError("TriMeshProcessRifts did not find any rift") … … 33 36 md.mesh.numberofvertices=numpy.size(md.mesh.x) 34 37 md.mesh.z=numpy.zeros(md.mesh.numberofvertices) 35 md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x) )36 md.mesh.vertexonboundary[md.mesh.segments[:,0:2] .astype(int)-1]=137 md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices )38 md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices )39 md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements )40 md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements )38 md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x),bool) 39 md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True 40 md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices,bool) 41 md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices,bool) 42 md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements,bool) 43 md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements,bool) 41 44 42 45 #get coordinates of rift tips -
issm/trunk/src/m/mesh/triangle.py
r13975 r14067 44 44 #Mesh using TriMesh 45 45 [md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers]=TriMesh(domainname,riftname,area) 46 md.mesh.elements=md.mesh.elements.astype(int) 47 md.mesh.segments=md.mesh.segments.astype(int) 48 md.mesh.segmentmarkers=md.mesh.segmentmarkers.astype(int) 46 49 47 50 #Fill in rest of fields: … … 49 52 md.mesh.numberofvertices = numpy.size(md.mesh.x) 50 53 md.mesh.z = numpy.zeros(md.mesh.numberofvertices) 51 md.mesh.vertexonboundary = numpy.zeros(md.mesh.numberofvertices )52 md.mesh.vertexonboundary[md.mesh.segments[:,0:2] .astype(int)-1] = 1.53 md.mesh.vertexonbed = numpy.ones(md.mesh.numberofvertices )54 md.mesh.vertexonsurface = numpy.ones(md.mesh.numberofvertices )55 md.mesh.elementonbed = numpy.ones(md.mesh.numberofelements )56 md.mesh.elementonsurface = numpy.ones(md.mesh.numberofelements )54 md.mesh.vertexonboundary = numpy.zeros(md.mesh.numberofvertices,bool) 55 md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1] = True 56 md.mesh.vertexonbed = numpy.ones(md.mesh.numberofvertices,bool) 57 md.mesh.vertexonsurface = numpy.ones(md.mesh.numberofvertices,bool) 58 md.mesh.elementonbed = numpy.ones(md.mesh.numberofelements,bool) 59 md.mesh.elementonsurface = numpy.ones(md.mesh.numberofelements,bool) 57 60 58 61 #Now, build the connectivity tables for this mesh. 59 [md.mesh.vertexconnectivity] = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)62 [md.mesh.vertexconnectivity] = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices) 60 63 [md.mesh.elementconnectivity] = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity) 61 64 -
issm/trunk/src/m/parameterization/contourenvelope.py
r13975 r14067 76 76 pos=numpy.nonzero(flags) 77 77 elemin[pos]=1 78 nodein[mesh.elements[pos,:] .astype(int)-1]=178 nodein[mesh.elements[pos,:]-1]=1 79 79 80 80 #modify element connectivity … … 88 88 if len(args)==1: 89 89 flag[numpy.nonzero(flag)]=elemin[flag[numpy.nonzero(flag)]] 90 elementonboundary=numpy.logical_and(numpy.prod(flag,axis=1)==0,numpy.sum(flag,axis=1)>0) .astype(float)90 elementonboundary=numpy.logical_and(numpy.prod(flag,axis=1)==0,numpy.sum(flag,axis=1)>0) 91 91 92 92 #Find segments on boundary 93 93 pos=numpy.nonzero(elementonboundary)[0] 94 94 num_segments=numpy.size(pos) 95 segments=numpy.zeros((num_segments*3,3) )95 segments=numpy.zeros((num_segments*3,3),int) 96 96 count=0 97 97 98 98 for el1 in pos: 99 els2=mesh.elementconnectivity[el1,numpy.nonzero(mesh.elementconnectivity[el1,:])[0]] .astype(int)-199 els2=mesh.elementconnectivity[el1,numpy.nonzero(mesh.elementconnectivity[el1,:])[0]]-1 100 100 if numpy.size(els2)>1: 101 101 flag=numpy.intersect1d(mesh.elements[els2[0],:],mesh.elements[els2[1],:]) -
issm/trunk/src/m/parameterization/setflowequation.py
r13975 r14067 50 50 #Flag the elements that have not been flagged as filltype 51 51 if strcmpi(filltype,'hutter'): 52 hutterflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(macayealflag,pattynflag)))]= 152 hutterflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(macayealflag,pattynflag)))]=True 53 53 elif strcmpi(filltype,'macayeal'): 54 macayealflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(pattynflag,stokesflag))))]= 154 macayealflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(pattynflag,stokesflag))))]=True 55 55 elif strcmpi(filltype,'pattyn'): 56 pattynflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(macayealflag,stokesflag))))]= 156 pattynflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(macayealflag,stokesflag))))]=True 57 57 58 58 #check that each element has at least one flag … … 63 63 if any(hutterflag+macayealflag+l1l2flag+pattynflag+stokesflag>1): 64 64 print "setflowequation warning message: some elements have several types, higher order type is used for them" 65 hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,macayealflag))]= 066 hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,pattynflag))]= 067 macayealflag[numpy.nonzero(numpy.logical_and(macayealflag,pattynflag))]= 065 hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,macayealflag))]=False 66 hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,pattynflag))]=False 67 macayealflag[numpy.nonzero(numpy.logical_and(macayealflag,pattynflag))]=False 68 68 69 69 #Check that no pattyn or stokes for 2d mesh … … 77 77 78 78 #Initialize node fields 79 nodeonhutter=numpy.zeros(md.mesh.numberofvertices )80 nodeonhutter[md.mesh.elements[numpy.nonzero(hutterflag),:] .astype(int)-1]=181 nodeonmacayeal=numpy.zeros(md.mesh.numberofvertices )82 nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:] .astype(int)-1]=183 nodeonl1l2=numpy.zeros(md.mesh.numberofvertices )84 nodeonl1l2[md.mesh.elements[numpy.nonzero(l1l2flag),:] .astype(int)-1]=185 nodeonpattyn=numpy.zeros(md.mesh.numberofvertices )86 nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:] .astype(int)-1]=187 nodeonstokes=numpy.zeros(md.mesh.numberofvertices )88 noneflag=numpy.zeros(md.mesh.numberofelements )79 nodeonhutter=numpy.zeros(md.mesh.numberofvertices,bool) 80 nodeonhutter[md.mesh.elements[numpy.nonzero(hutterflag),:]-1]=True 81 nodeonmacayeal=numpy.zeros(md.mesh.numberofvertices,bool) 82 nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True 83 nodeonl1l2=numpy.zeros(md.mesh.numberofvertices,bool) 84 nodeonl1l2[md.mesh.elements[numpy.nonzero(l1l2flag),:]-1]=True 85 nodeonpattyn=numpy.zeros(md.mesh.numberofvertices,bool) 86 nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True 87 nodeonstokes=numpy.zeros(md.mesh.numberofvertices,bool) 88 noneflag=numpy.zeros(md.mesh.numberofelements,bool) 89 89 90 90 #First modify stokesflag to get rid of elements contrained everywhere (spc + border with pattyn or macayeal) … … 96 96 numpy.logical_and(nodeonpattyn,nodeonstokes).reshape(-1,1)).astype(int) #find all the nodes on the boundary of the domain without icefront 97 97 # fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6); %find all the nodes on the boundary of the domain without icefront 98 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 icefront99 stokesflag[numpy.nonzero(fullspcelems.reshape(-1))]= 0100 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:] .astype(int)-1]=198 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 99 stokesflag[numpy.nonzero(fullspcelems.reshape(-1))]=False 100 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True 101 101 102 102 #Then complete with NoneApproximation or the other model used if there is no stokes 103 103 if any(stokesflag): 104 104 if any(pattynflag): #fill with pattyn 105 pattynflag[numpy.logical_not(stokesflag)]= 1106 nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:] .astype(int)-1]=1105 pattynflag[numpy.logical_not(stokesflag)]=True 106 nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True 107 107 elif any(macayealflag): #fill with macayeal 108 macayealflag[numpy.logical_not(stokesflag)]= 1109 nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:] .astype(int)-1]=1108 macayealflag[numpy.logical_not(stokesflag)]=True 109 nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True 110 110 else: #fill with none 111 noneflag[numpy.nonzero(numpy.logical_not(stokesflag))]= 1111 noneflag[numpy.nonzero(numpy.logical_not(stokesflag))]=True 112 112 113 113 #Now take care of the coupling between MacAyeal and Pattyn 114 114 md.diagnostic.vertex_pairing=numpy.array([]) 115 nodeonmacayealpattyn=numpy.zeros(md.mesh.numberofvertices )116 nodeonpattynstokes=numpy.zeros(md.mesh.numberofvertices )117 nodeonmacayealstokes=numpy.zeros(md.mesh.numberofvertices )118 macayealpattynflag=numpy.zeros(md.mesh.numberofelements )119 macayealstokesflag=numpy.zeros(md.mesh.numberofelements )120 pattynstokesflag=numpy.zeros(md.mesh.numberofelements )115 nodeonmacayealpattyn=numpy.zeros(md.mesh.numberofvertices,bool) 116 nodeonpattynstokes=numpy.zeros(md.mesh.numberofvertices,bool) 117 nodeonmacayealstokes=numpy.zeros(md.mesh.numberofvertices,bool) 118 macayealpattynflag=numpy.zeros(md.mesh.numberofelements,bool) 119 macayealstokesflag=numpy.zeros(md.mesh.numberofelements,bool) 120 pattynstokesflag=numpy.zeros(md.mesh.numberofelements,bool) 121 121 if strcmpi(coupling_method,'penalties'): 122 122 #Create the border nodes between Pattyn and MacAyeal and extrude them … … 135 135 if any(macayealflag) and any(pattynflag): #coupling macayeal pattyn 136 136 #Find node at the border 137 nodeonmacayealpattyn[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonpattyn))]= 1137 nodeonmacayealpattyn[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonpattyn))]=True 138 138 #Macayeal elements in contact with this layer become MacAyealPattyn elements 139 matrixelements=ismember(md.mesh.elements .astype(int)-1,numpy.nonzero(nodeonmacayealpattyn)[0])139 matrixelements=ismember(md.mesh.elements-1,numpy.nonzero(nodeonmacayealpattyn)[0]) 140 140 commonelements=numpy.sum(matrixelements,axis=1)!=0 141 commonelements[numpy.nonzero(pattynflag)]= 0#only one layer: the elements previously in macayeal142 macayealflag[numpy.nonzero(commonelements)]= 0#these elements are now macayealpattynelements143 macayealpattynflag[numpy.nonzero(commonelements)]= 1144 nodeonmacayeal[:]= 0145 nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:] .astype(int)-1]=1141 commonelements[numpy.nonzero(pattynflag)]=False #only one layer: the elements previously in macayeal 142 macayealflag[numpy.nonzero(commonelements)]=False #these elements are now macayealpattynelements 143 macayealpattynflag[numpy.nonzero(commonelements)]=True 144 nodeonmacayeal[:]=False 145 nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True 146 146 147 147 #rule out elements that don't touch the 2 boundaries 148 148 pos=numpy.nonzero(macayealpattynflag)[0] 149 149 elist=numpy.zeros(numpy.size(pos),dtype=int) 150 elist = elist + numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:] .astype(int)-1],axis=1).astype(bool)151 elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:] .astype(int)-1] ,axis=1).astype(bool)150 elist = elist + numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:]-1],axis=1).astype(bool) 151 elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:]-1] ,axis=1).astype(bool) 152 152 pos1=numpy.nonzero(elist==1)[0] 153 macayealflag[pos[pos1]]= 1154 macayealpattynflag[pos[pos1]]= 0153 macayealflag[pos[pos1]]=True 154 macayealpattynflag[pos[pos1]]=False 155 155 pos2=numpy.nonzero(elist==-1)[0] 156 pattynflag[pos[pos2]]= 1157 macayealpattynflag[pos[pos2]]= 0156 pattynflag[pos[pos2]]=True 157 macayealpattynflag[pos[pos2]]=False 158 158 159 159 #Recompute nodes associated to these elements 160 nodeonmacayeal[:]= 0161 nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:] .astype(int)-1]=1162 nodeonpattyn[:]= 0163 nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:] .astype(int)-1]=1164 nodeonmacayealpattyn[:]= 0165 nodeonmacayealpattyn[md.mesh.elements[numpy.nonzero(macayealpattynflag),:] .astype(int)-1]=1160 nodeonmacayeal[:]=False 161 nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True 162 nodeonpattyn[:]=False 163 nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True 164 nodeonmacayealpattyn[:]=False 165 nodeonmacayealpattyn[md.mesh.elements[numpy.nonzero(macayealpattynflag),:]-1]=True 166 166 167 167 elif any(pattynflag) and any(stokesflag): #coupling pattyn stokes 168 168 #Find node at the border 169 nodeonpattynstokes[numpy.nonzero(numpy.logical_and(nodeonpattyn,nodeonstokes))]= 1169 nodeonpattynstokes[numpy.nonzero(numpy.logical_and(nodeonpattyn,nodeonstokes))]=True 170 170 #Stokes elements in contact with this layer become PattynStokes elements 171 matrixelements=ismember(md.mesh.elements .astype(int)-1,numpy.nonzero(nodeonpattynstokes)[0])171 matrixelements=ismember(md.mesh.elements-1,numpy.nonzero(nodeonpattynstokes)[0]) 172 172 commonelements=numpy.sum(matrixelements,axis=1)!=0 173 commonelements[numpy.nonzero(pattynflag)]= 0#only one layer: the elements previously in macayeal174 stokesflag[numpy.nonzero(commonelements)]= 0#these elements are now macayealpattynelements175 pattynstokesflag[numpy.nonzero(commonelements)]= 1176 nodeonstokes=numpy.zeros(md.mesh.numberofvertices )177 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:] .astype(int)-1]=1173 commonelements[numpy.nonzero(pattynflag)]=False #only one layer: the elements previously in macayeal 174 stokesflag[numpy.nonzero(commonelements)]=False #these elements are now macayealpattynelements 175 pattynstokesflag[numpy.nonzero(commonelements)]=True 176 nodeonstokes=numpy.zeros(md.mesh.numberofvertices,bool) 177 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True 178 178 179 179 #rule out elements that don't touch the 2 boundaries 180 180 pos=numpy.nonzero(pattynstokesflag)[0] 181 181 elist=numpy.zeros(numpy.size(pos),dtype=int) 182 elist = elist + numpy.sum(nodeonstokes[md.mesh.elements[pos,:] .astype(int)-1],axis=1).astype(bool)183 elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:] .astype(int)-1],axis=1).astype(bool)182 elist = elist + numpy.sum(nodeonstokes[md.mesh.elements[pos,:]-1],axis=1).astype(bool) 183 elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:]-1],axis=1).astype(bool) 184 184 pos1=numpy.nonzero(elist==1)[0] 185 stokesflag[pos[pos1]]= 1186 pattynstokesflag[pos[pos1]]= 0185 stokesflag[pos[pos1]]=True 186 pattynstokesflag[pos[pos1]]=False 187 187 pos2=numpy.nonzero(elist==-1)[0] 188 pattynflag[pos[pos2]]= 1189 pattynstokesflag[pos[pos2]]= 0188 pattynflag[pos[pos2]]=True 189 pattynstokesflag[pos[pos2]]=False 190 190 191 191 #Recompute nodes associated to these elements 192 nodeonstokes[:]= 0193 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:] .astype(int)-1]=1194 nodeonpattyn[:]= 0195 nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:] .astype(int)-1]=1196 nodeonpattynstokes[:]= 0197 nodeonpattynstokes[md.mesh.elements[numpy.nonzero(pattynstokesflag),:] .astype(int)-1]=1192 nodeonstokes[:]=False 193 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True 194 nodeonpattyn[:]=False 195 nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True 196 nodeonpattynstokes[:]=False 197 nodeonpattynstokes[md.mesh.elements[numpy.nonzero(pattynstokesflag),:]-1]=True 198 198 199 199 elif any(stokesflag) and any(macayealflag): 200 200 #Find node at the border 201 nodeonmacayealstokes[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonstokes))]= 1201 nodeonmacayealstokes[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonstokes))]=True 202 202 #Stokes elements in contact with this layer become MacAyealStokes elements 203 matrixelements=ismember(md.mesh.elements .astype(int)-1,numpy.nonzero(nodeonmacayealstokes)[0])203 matrixelements=ismember(md.mesh.elements-1,numpy.nonzero(nodeonmacayealstokes)[0]) 204 204 commonelements=numpy.sum(matrixelements,axis=1)!=0 205 commonelements[numpy.nonzero(macayealflag)]= 0#only one layer: the elements previously in macayeal206 stokesflag[numpy.nonzero(commonelements)]= 0#these elements are now macayealmacayealelements207 macayealstokesflag[numpy.nonzero(commonelements)]= 1208 nodeonstokes=numpy.zeros(md.mesh.numberofvertices )209 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:] .astype(int)-1]=1205 commonelements[numpy.nonzero(macayealflag)]=False #only one layer: the elements previously in macayeal 206 stokesflag[numpy.nonzero(commonelements)]=False #these elements are now macayealmacayealelements 207 macayealstokesflag[numpy.nonzero(commonelements)]=True 208 nodeonstokes=numpy.zeros(md.mesh.numberofvertices,bool) 209 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True 210 210 211 211 #rule out elements that don't touch the 2 boundaries 212 212 pos=numpy.nonzero(macayealstokesflag)[0] 213 213 elist=numpy.zeros(numpy.size(pos),dtype=int) 214 elist = elist + numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:] .astype(int)-1],axis=1).astype(bool)215 elist = elist - numpy.sum(nodeonstokes[md.mesh.elements[pos,:] .astype(int)-1] ,axis=1).astype(bool)214 elist = elist + numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:]-1],axis=1).astype(bool) 215 elist = elist - numpy.sum(nodeonstokes[md.mesh.elements[pos,:]-1] ,axis=1).astype(bool) 216 216 pos1=numpy.nonzero(elist==1)[0] 217 macayealflag[pos[pos1]]= 1218 macayealstokesflag[pos[pos1]]= 0217 macayealflag[pos[pos1]]=True 218 macayealstokesflag[pos[pos1]]=False 219 219 pos2=numpy.nonzero(elist==-1)[0] 220 stokesflag[pos[pos2]]= 1221 macayealstokesflag[pos[pos2]]= 0220 stokesflag[pos[pos2]]=True 221 macayealstokesflag[pos[pos2]]=False 222 222 223 223 #Recompute nodes associated to these elements 224 nodeonmacayeal[:]= 0225 nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:] .astype(int)-1]=1226 nodeonstokes[:]= 0227 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:] .astype(int)-1]=1228 nodeonmacayealstokes[:]= 0229 nodeonmacayealstokes[md.mesh.elements[numpy.nonzero(macayealstokesflag),:] .astype(int)-1]=1224 nodeonmacayeal[:]=False 225 nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True 226 nodeonstokes[:]=False 227 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True 228 nodeonmacayealstokes[:]=False 229 nodeonmacayealstokes[md.mesh.elements[numpy.nonzero(macayealstokesflag),:]-1]=True 230 230 231 231 elif any(stokesflag) and any(hutterflag): 232 232 raise TypeError("type of coupling not supported yet") 233 233 234 #Create Mac aAyealPattynApproximation where needed235 md.flowequation.element_equation=numpy.zeros(md.mesh.numberofelements )234 #Create MacAyealPattynApproximation where needed 235 md.flowequation.element_equation=numpy.zeros(md.mesh.numberofelements,int) 236 236 md.flowequation.element_equation[numpy.nonzero(noneflag)]=0 237 237 md.flowequation.element_equation[numpy.nonzero(hutterflag)]=1 … … 250 250 251 251 #Create vertices_type 252 md.flowequation.vertex_equation=numpy.zeros(md.mesh.numberofvertices )252 md.flowequation.vertex_equation=numpy.zeros(md.mesh.numberofvertices,int) 253 253 pos=numpy.nonzero(nodeonmacayeal) 254 254 md.flowequation.vertex_equation[pos]=2 … … 273 273 274 274 #figure out solution types 275 md.flowequation.ishutter= float(any(md.flowequation.element_equation==1))276 md.flowequation.ismacayealpattyn= float(numpy.any(numpy.logical_or(md.flowequation.element_equation==2,md.flowequation.element_equation==3)))277 md.flowequation.isl1l2= float(any(md.flowequation.element_equation==8))278 md.flowequation.isstokes= float(any(md.flowequation.element_equation==4))275 md.flowequation.ishutter=any(md.flowequation.element_equation==1) 276 md.flowequation.ismacayealpattyn=bool(numpy.any(numpy.logical_or(md.flowequation.element_equation==2,md.flowequation.element_equation==3))) 277 md.flowequation.isl1l2=any(md.flowequation.element_equation==8) 278 md.flowequation.isstokes=any(md.flowequation.element_equation==4) 279 279 280 280 return md -
issm/trunk/src/m/parameterization/setmask.py
r13975 r14067 37 37 vertexonfloatingice = numpy.zeros(md.mesh.numberofvertices,'bool') 38 38 vertexongroundedice = numpy.zeros(md.mesh.numberofvertices,'bool') 39 vertexongroundedice[md.mesh.elements[numpy.nonzero(elementongroundedice),:] .astype(int)-1]=True39 vertexongroundedice[md.mesh.elements[numpy.nonzero(elementongroundedice),:]-1]=True 40 40 vertexonfloatingice[numpy.nonzero(numpy.logical_not(vertexongroundedice))]=True 41 41 #}}} 42 42 43 43 #Return: 44 md.mask.elementonfloatingice = elementonfloatingice .astype(float)45 md.mask.vertexonfloatingice = vertexonfloatingice .astype(float)46 md.mask.elementongroundedice = elementongroundedice .astype(float)47 md.mask.vertexongroundedice = vertexongroundedice .astype(float)48 md.mask.vertexonwater = numpy.zeros(md.mesh.numberofvertices )49 md.mask.elementonwater = numpy.zeros(md.mesh.numberofelements )44 md.mask.elementonfloatingice = elementonfloatingice 45 md.mask.vertexonfloatingice = vertexonfloatingice 46 md.mask.elementongroundedice = elementongroundedice 47 md.mask.vertexongroundedice = vertexongroundedice 48 md.mask.vertexonwater = numpy.zeros(md.mesh.numberofvertices,bool) 49 md.mask.elementonwater = numpy.zeros(md.mesh.numberofelements,bool) 50 50 51 51 return md -
issm/trunk/src/m/plot/applyoptions.m
r13975 r14067 205 205 if exist(options,'log') 206 206 logval= getfieldvalue(options,'log'); 207 for i= 1: size(tick_vals,1)207 for i= 1:numel(tick_vals) 208 208 tick_vals(i)= logval^(tick_vals(i)); 209 209 end 210 elseif size(tick_vals,1) == 3210 elseif numel(tick_vals) == 3 211 211 tick_vals=[tick_vals(1); mean(tick_vals(1:2)); tick_vals(2); ... 212 212 mean(tick_vals(2:3)); tick_vals(3)]; 213 set(c,'YTick',tick_vals); 214 end 215 else 216 if exist(options,'log') 217 logvalue=getfieldvalue(options,'log'); 218 set(c,'YTick',log(tick_vals)./log(logvalue)); 219 else 213 220 set(c,'YTick',tick_vals); 214 221 end -
issm/trunk/src/m/plot/colormaps/getcolormap.m
r13975 r14067 24 24 map = map(128:end,:); 25 25 elseif strcmpi(map,'redblue'), 26 map = hsv ;26 map = hsv(256); 27 27 map = rgb2hsv(map); 28 28 map(:,2) = max(min( abs(map(:,1)-0.5)/0.5 ,1),0); -
issm/trunk/src/m/plot/radarpower.m
r13975 r14067 20 20 end 21 21 22 geotiff_name = getfieldvalue(options,'geotiff_name',''); 22 23 highres = getfieldvalue(options,'highres',0); 23 24 xlim = getfieldvalue(options,'xlim',[min(md.mesh.x) max(md.mesh.x)]); … … 64 65 %md.radaroverlay.x=(xmin:(xmax-xmin)/(size(md.radaroverlay.pwr,2)-1):xmax); 65 66 %md.radaroverlay.y=(ymin:(ymax-ymin)/(size(md.radaroverlay.pwr,1)-1):ymax); 66 if highres, 67 if ~exist(['/u/astrid-r1b/ModelData/MOG/mog100_r2_hp1.tif']), 68 error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MOG/mog100_r2_hp1.tif not found.']); 67 if ~exist(options,'geotiff_name'), 68 if highres, 69 if ~exist(['/u/astrid-r1b/ModelData/MOG/mog100_r2_hp1.tif']), 70 error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MOG/mog100_r2_hp1.tif not found.']); 71 end 72 geotiff_name=['/u/astrid-r1b/ModelData/MOG/mog100_r2_hp1.tif']; 73 else 74 if ~exist(['/u/astrid-r1b/ModelData/MOG/mog500_r2_hp1.tif']), 75 error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MOG/mog500_r2_hp1.tif not found.']); 76 end 77 geotiff_name=['/u/astrid-r1b/ModelData/MOG/mog500_r2_hp1.tif']; 69 78 end 70 geotiff_name=['/u/astrid-r1b/ModelData/MOG/mog100_r2_hp1.tif'];71 else72 if ~exist(['/u/astrid-r1b/ModelData/MOG/mog500_r2_hp1.tif']),73 error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MOG/mog500_r2_hp1.tif not found.']);74 end75 geotiff_name=['/u/astrid-r1b/ModelData/MOG/mog500_r2_hp1.tif'];76 79 end 77 80 … … 91 94 92 95 elseif strcmpi(md.mesh.hemisphere,'s'), 93 if highres, 94 if ~exist(['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif']), 95 error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif not found.']); 96 if ~exist(options,'geotiff_name'), 97 if highres, 98 if ~exist(['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif']), 99 error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif not found.']); 100 end 101 geotiff_name=['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif']; 102 else 103 if ~exist(['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif']), 104 error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif not found.']); 105 end 106 geotiff_name=['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif']; 96 107 end 97 geotiff_name=['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif'];98 else99 if ~exist(['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif']),100 error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif not found.']);101 end102 geotiff_name=['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif'];103 108 end 104 109 -
issm/trunk/src/m/solve/WriteData.py
r13975 r14067 105 105 elif isinstance(data,(list,tuple)): 106 106 data=numpy.array(data).reshape(-1,1) 107 if len(data.shape) == 1:107 if numpy.ndim(data) == 1: 108 108 if numpy.size(data): 109 109 data=data.reshape(numpy.size(data),1) … … 138 138 elif isinstance(data,(list,tuple)): 139 139 data=numpy.array(data).reshape(-1,1) 140 if len(data.shape) == 1:140 if numpy.ndim(data) == 1: 141 141 if numpy.size(data): 142 142 data=data.reshape(numpy.size(data),1) … … 171 171 elif isinstance(data,(list,tuple)): 172 172 data=numpy.array(data).reshape(-1,1) 173 if len(data.shape) == 1:173 if numpy.ndim(data) == 1: 174 174 if numpy.size(data): 175 175 data=data.reshape(numpy.size(data),1) … … 207 207 elif isinstance(matrix,(list,tuple)): 208 208 matrix=numpy.array(matrix).reshape(-1,1) 209 if len(matrix.shape) == 1:209 if numpy.ndim(matrix) == 1: 210 210 if numpy.size(matrix): 211 211 matrix=matrix.reshape(numpy.size(matrix),1) … … 231 231 elif isinstance(matrix,(list,tuple)): 232 232 matrix=numpy.array(matrix).reshape(-1,1) 233 if len(matrix.shape) == 1:233 if numpy.ndim(matrix) == 1: 234 234 matrix=matrix.reshape(numpy.size(matrix),1) 235 235 -
issm/trunk/src/wrappers/ElementConnectivity/ElementConnectivity.cpp
r13236 r14067 13 13 14 14 /*inputs: */ 15 double* elements=NULL;16 double* nodeconnectivity=NULL;17 int nel,nods;18 int 15 int* elements=NULL; 16 int* nodeconnectivity=NULL; 17 int nels,nods; 18 int width; 19 19 20 20 /*outputs: */ 21 double* elementconnectivity=NULL;21 int* elementconnectivity=NULL; 22 22 23 23 /*Boot module: */ … … 28 28 29 29 /*Input datasets: */ 30 FetchData(&elements,&nel ,NULL,ELEMENTS);30 FetchData(&elements,&nels,NULL,ELEMENTS); 31 31 FetchData(&nodeconnectivity,&nods,&width,NODECONNECTIVITY); 32 32 33 33 /*!Generate internal degree of freedom numbers: */ 34 ElementConnectivityx(&elementconnectivity, elements,nel, nodeconnectivity, nods,width);34 ElementConnectivityx(&elementconnectivity,elements,nels,nodeconnectivity,nods,width); 35 35 36 36 /*write output datasets: */ 37 WriteData(ELEMENTCONNECTIVITY,elementconnectivity,nel ,3);37 WriteData(ELEMENTCONNECTIVITY,elementconnectivity,nels,3); 38 38 39 39 /*end module: */ -
issm/trunk/src/wrappers/Kriging/Kriging.cpp
r13250 r14067 5 5 6 6 void KrigingUsage(void){/*{{{*/ 7 int num=1; 8 #ifdef _MULTITHREADING_ 9 num=_NUMTHREADS_; 10 #endif 7 11 _pprintLine_(""); 8 12 _pprintLine_(" usage: predictions=" << __FUNCT__ << "(x,y,observations,x_interp,y_interp,'options');"); … … 21 25 _pprintLine_(" -'mintrimming': minimum trimming value (default is +1.e+21)"); 22 26 _pprintLine_(" -'minspacing': minimum distance between observation (default is 0.01)"); 27 _pprintLine_(" -'numthreads': number of threads, default is "<<num); 23 28 _pprintLine_(""); 24 29 }/*}}}*/ -
issm/trunk/src/wrappers/NodeConnectivity/NodeConnectivity.cpp
r13236 r14067 13 13 14 14 /*inputs: */ 15 double* elements=NULL;16 int nel;17 int 15 int* elements=NULL; 16 int nels; 17 int nods; 18 18 19 19 /*outputs: */ 20 double* connectivity=NULL;21 int 20 int* connectivity=NULL; 21 int width; 22 22 23 23 /*Boot module: */ … … 28 28 29 29 /*Input datasets: */ 30 FetchData(&elements,&nel ,NULL,ELEMENTS);30 FetchData(&elements,&nels,NULL,ELEMENTS); 31 31 FetchData(&nods,NUMNODES); 32 32 33 33 /*!Generate internal degree of freedom numbers: */ 34 NodeConnectivityx(&connectivity, &width,elements,nel,nods);34 NodeConnectivityx(&connectivity,&width,elements,nels,nods); 35 35 36 36 /*write output datasets: */ -
issm/trunk/src/wrappers/matlab/io/FetchMatlabData.cpp
r13749 r14067 264 264 265 265 double* outvector=NULL; 266 int outvector_rows; 267 268 if(mxIsEmpty(dataref)){ 269 /*Nothing to pick up. Just initialize matrix pointer to NULL: */ 270 outvector_rows=0; 271 outvector=NULL; 272 } 273 else if (mxIsClass(dataref,"double") ){ 274 275 /*Convert matlab vector to double* vector: */ 276 MatlabVectorToDoubleVector(&outvector,&outvector_rows,dataref); 277 278 } 279 else{ 280 /*This is an error: we don't have the correct input!: */ 281 _error_("Input parameter of class " << mxGetClassName(dataref) << " not supported yet"); 266 int M,N; 267 268 /*Use Fetch matrix*/ 269 FetchData(&outvector,&M,&N,dataref) ; 270 271 /*Check that it is a vector*/ 272 if(M>0){ 273 if(N!=1) _error_("input vector of size " << M << "x" << N << " should have only one column"); 282 274 } 283 275 284 276 /*Assign output pointers:*/ 285 277 *pvector=outvector; 286 if (pM)*pM=outvector_rows;278 if(pM)*pM=M; 287 279 } 288 280 /*}}}*/ -
issm/trunk/src/wrappers/python/io/FetchPythonData.cpp
r13749 r14067 67 67 npy_intp* dims=NULL; 68 68 69 /*retrive dimensions: */ 69 /*intermediary:*/ 70 long* lmatrix=NULL; 71 bool* bmatrix=NULL; 72 int i; 73 74 /*retrieve dimensions: */ 70 75 ndim=PyArray_NDIM((const PyArrayObject*)py_matrix); 71 76 if(ndim!=2)_error_("expecting an MxN matrix in input!"); … … 74 79 75 80 if (M && N) { 76 /*retrieve internal value: */ 77 dmatrix=(double*)PyArray_DATA((PyArrayObject*)py_matrix); 78 79 /*copy matrix: */ 80 matrix=xNew<double>(M*N); 81 memcpy(matrix,dmatrix,(M*N)*sizeof(double)); 81 if (PyArray_TYPE((PyArrayObject*)py_matrix) == NPY_DOUBLE) { 82 /*retrieve internal value: */ 83 dmatrix=(double*)PyArray_DATA((PyArrayObject*)py_matrix); 84 85 /*copy matrix: */ 86 matrix=xNew<double>(M*N); 87 memcpy(matrix,dmatrix,(M*N)*sizeof(double)); 88 } 89 90 else if (PyArray_TYPE((PyArrayObject*)py_matrix) == NPY_INT64) { 91 /*retrieve internal value: */ 92 lmatrix=(long*)PyArray_DATA((PyArrayObject*)py_matrix); 93 94 /*transform into double matrix: */ 95 matrix=xNew<double>(M*N); 96 for(i=0;i<M*N;i++)matrix[i]=(double)lmatrix[i]; 97 } 98 99 else if (PyArray_TYPE((PyArrayObject*)py_matrix) == NPY_BOOL) { 100 /*retrieve internal value: */ 101 bmatrix=(bool*)PyArray_DATA((PyArrayObject*)py_matrix); 102 103 /*transform into double matrix: */ 104 matrix=xNew<double>(M*N); 105 for(i=0;i<M*N;i++)matrix[i]=(double)bmatrix[i]; 106 } 107 108 else 109 _error_("unrecognized matrix type in input!"); 82 110 } 83 111 else … … 94 122 95 123 /*output: */ 96 double* dmatrix=NULL;97 124 int* matrix=NULL; 98 125 int M,N; 99 100 /*intermediary:*/101 int i;102 126 int ndim; 103 127 npy_intp* dims=NULL; 104 128 105 /*retrive dimensions: */ 129 /*intermediary:*/ 130 double* dmatrix=NULL; 131 long* lmatrix=NULL; 132 bool* bmatrix=NULL; 133 int i; 134 135 /*retrieve dimensions: */ 106 136 ndim=PyArray_NDIM((const PyArrayObject*)py_matrix); 107 137 if(ndim!=2)_error_("expecting an MxN matrix in input!"); … … 110 140 111 141 if (M && N) { 112 /*retrieve internal value: */ 113 dmatrix=(double*)PyArray_DATA((PyArrayObject*)py_matrix); 114 115 /*transform into integer matrix: */ 116 matrix=xNew<int>(M*N); 117 for(i=0;i<M*N;i++)matrix[i]=(int)dmatrix[i]; 142 if (PyArray_TYPE((PyArrayObject*)py_matrix) == NPY_DOUBLE) { 143 /*retrieve internal value: */ 144 dmatrix=(double*)PyArray_DATA((PyArrayObject*)py_matrix); 145 146 /*transform into integer matrix: */ 147 matrix=xNew<int>(M*N); 148 for(i=0;i<M*N;i++)matrix[i]=(int)dmatrix[i]; 149 } 150 151 else if (PyArray_TYPE((PyArrayObject*)py_matrix) == NPY_INT64) { 152 /*retrieve internal value: */ 153 lmatrix=(long*)PyArray_DATA((PyArrayObject*)py_matrix); 154 155 /*transform into integer matrix: */ 156 matrix=xNew<int>(M*N); 157 for(i=0;i<M*N;i++)matrix[i]=(int)lmatrix[i]; 158 } 159 160 else if (PyArray_TYPE((PyArrayObject*)py_matrix) == NPY_BOOL) { 161 /*retrieve internal value: */ 162 bmatrix=(bool*)PyArray_DATA((PyArrayObject*)py_matrix); 163 164 /*transform into integer matrix: */ 165 matrix=xNew<int>(M*N); 166 for(i=0;i<M*N;i++)matrix[i]=(int)bmatrix[i]; 167 } 168 169 else 170 _error_("unrecognized matrix type in input!"); 118 171 } 119 172 else … … 136 189 npy_intp* dims=NULL; 137 190 138 /*retrive dimensions: */ 191 /*intermediary:*/ 192 long* lvector=NULL; 193 bool* bvector=NULL; 194 int i; 195 196 /*retrieve dimensions: */ 139 197 ndim=PyArray_NDIM((const PyArrayObject*)py_vector); 140 198 if(ndim!=1)_error_("expecting an Mx1 vector in input!"); … … 143 201 144 202 if (M) { 145 /*retrieve internal value: */ 146 dvector=(double*)PyArray_DATA((PyArrayObject*)py_vector); 147 148 /*copy vector: */ 149 vector=xNew<double>(M); 150 memcpy(vector,dvector,(M)*sizeof(double)); 203 if (PyArray_TYPE((PyArrayObject*)py_vector) == NPY_DOUBLE) { 204 /*retrieve internal value: */ 205 dvector=(double*)PyArray_DATA((PyArrayObject*)py_vector); 206 207 /*copy vector: */ 208 vector=xNew<double>(M); 209 memcpy(vector,dvector,(M)*sizeof(double)); 210 } 211 212 else if (PyArray_TYPE((PyArrayObject*)py_vector) == NPY_INT64) { 213 /*retrieve internal value: */ 214 lvector=(long*)PyArray_DATA((PyArrayObject*)py_vector); 215 216 /*transform into double vector: */ 217 vector=xNew<double>(M); 218 for(i=0;i<M;i++)vector[i]=(double)lvector[i]; 219 } 220 221 else if (PyArray_TYPE((PyArrayObject*)py_vector) == NPY_BOOL) { 222 /*retrieve internal value: */ 223 bvector=(bool*)PyArray_DATA((PyArrayObject*)py_vector); 224 225 /*transform into double vector: */ 226 vector=xNew<double>(M); 227 for(i=0;i<M;i++)vector[i]=(double)bvector[i]; 228 } 229 230 else 231 _error_("unrecognized vector type in input!"); 151 232 } 152 233 else -
issm/trunk/src/wrappers/python/io/WritePythonData.cpp
r13749 r14067 40 40 dims[1]=(npy_intp)N; 41 41 array=PyArray_SimpleNewFromData(2,dims,NPY_DOUBLE,matrix); 42 43 PyTuple_SetItem(tuple, index, array); 44 }/*}}}*/ 45 /*FUNCTION WriteData(PyObject* py_tuple,int index, int* matrix, int M, int N){{{*/ 46 void WriteData(PyObject* tuple, int index, int* matrix, int M,int N){ 47 48 npy_intp dims[2]={0,0}; 49 PyObject* array=NULL; 50 51 /*intermediary:*/ 52 long* lmatrix=NULL; 53 int i; 54 55 /*transform into long matrix: */ 56 lmatrix=xNew<long>(M*N); 57 for(i=0;i<M*N;i++)lmatrix[i]=(long)matrix[i]; 58 59 dims[0]=(npy_intp)M; 60 dims[1]=(npy_intp)N; 61 array=PyArray_SimpleNewFromData(2,dims,NPY_INT64,lmatrix); 62 63 PyTuple_SetItem(tuple, index, array); 64 }/*}}}*/ 65 /*FUNCTION WriteData(PyObject* py_tuple,int index, bool* matrix, int M, int N){{{*/ 66 void WriteData(PyObject* tuple, int index, bool* matrix, int M,int N){ 67 68 npy_intp dims[2]={0,0}; 69 PyObject* array=NULL; 70 71 dims[0]=(npy_intp)M; 72 dims[1]=(npy_intp)N; 73 array=PyArray_SimpleNewFromData(2,dims,NPY_BOOL,matrix); 42 74 43 75 PyTuple_SetItem(tuple, index, array); -
issm/trunk/src/wrappers/python/io/pythonio.h
r13749 r14067 18 18 19 19 void WriteData(PyObject* py_tuple,int index, double* matrix, int M,int N); 20 void WriteData(PyObject* py_tuple,int index, int* matrix, int M,int N); 21 void WriteData(PyObject* py_tuple,int index, bool* matrix, int M,int N); 20 22 void WriteData(PyObject* py_tuple,int index, int integer); 21 23 void WriteData(PyObject* py_tuple,int index, char* string); -
issm/trunk/test
- Property svn:mergeinfo changed
/issm/trunk-jpl/test merged: 13976,13982,14003,14021-14022
- Property svn:mergeinfo changed
-
issm/trunk/test/NightlyRun/runme.m
r13975 r14067 97 97 if strcmpi(benchmark,'nightly'), 98 98 test_ids=intersect(test_ids,[1:999]); 99 elseif strcmpi(benchmark,'validation'), 100 test_ids=intersect(test_ids,[1001:1999]); 99 101 elseif strcmpi(benchmark,'ismip'), 100 102 test_ids=intersect(test_ids,[1101:1199]); … … 105 107 elseif strcmpi(benchmark,'mesh'), 106 108 test_ids=intersect(test_ids,[1401:1499]); 109 elseif strcmpi(benchmark,'tranforcing'), 110 test_ids=intersect(test_ids,[1501:1502]); 111 elseif strcmpi(benchmark,'referential'), 112 test_ids=intersect(test_ids,[1601:1602]); 107 113 elseif strcmpi(benchmark,'adolc'), 108 114 test_ids=intersect(test_ids,[3001:3020]); 109 elseif strcmpi(benchmark,'validation'),110 test_ids=intersect(test_ids,[1001:1999]);111 elseif strcmpi(benchmark,'tranforcing'),112 test_ids=intersect(test_ids,[1501:1502]);113 115 end 114 116 % }}} -
issm/trunk/test/NightlyRun/runme.py
r13975 r14067 99 99 if strcmpi(benchmark,'nightly'): 100 100 test_ids=test_ids.intersection(set(range(1,1000))) 101 elif strcmpi(benchmark,'validation'): 102 test_ids=test_ids.intersection(set(range(1001,2000))) 101 103 elif strcmpi(benchmark,'ismip'): 102 104 test_ids=test_ids.intersection(set(range(1101,1200))) … … 107 109 elif strcmpi(benchmark,'mesh'): 108 110 test_ids=test_ids.intersection(set(range(1401,1500))) 111 elif strcmpi(benchmark,'tranforcing'): 112 test_ids=test_ids.intersection(set(range(1501,1503))) 113 elif strcmpi(benchmark,'referential'): 114 test_ids=test_ids.intersection(set(range(1601,1603))) 109 115 elif strcmpi(benchmark,'adolc'): 110 116 test_ids=test_ids.intersection(set(range(3001,3020))) 111 elif strcmpi(benchmark,'validation'):112 test_ids=test_ids.intersection(set(range(1001,2000)))113 elif strcmpi(benchmark,'tranforcing'):114 test_ids=test_ids.intersection(set(range(1501,1503)))115 117 #print 'test_ids after benchmark =',test_ids 116 118 test_ids=list(test_ids) -
issm/trunk/test/NightlyRun/test109.py
r13975 r14067 14 14 md=setflowequation(md,'macayeal','all') 15 15 md.cluster=generic('name',oshostname(),'np',3) 16 md.transient.isdiagnostic= 017 md.transient.isprognostic= 018 md.transient.isthermal= 119 md.transient.isgroundingline= 016 md.transient.isdiagnostic=False 17 md.transient.isprognostic=False 18 md.transient.isthermal=True 19 md.transient.isgroundingline=False 20 20 md=solve(md,TransientSolutionEnum()) 21 21 -
issm/trunk/test/NightlyRun/test121.py
r13975 r14067 15 15 md.cluster=generic('name',oshostname(),'np',3); 16 16 md.initialization.waterfraction=numpy.zeros((md.mesh.numberofvertices,1)) 17 md.transient.isdiagnostic= 018 md.transient.isprognostic= 019 md.transient.isthermal= 120 md.transient.isgroundingline= 017 md.transient.isdiagnostic=False 18 md.transient.isprognostic=False 19 md.transient.isthermal=True 20 md.transient.isgroundingline=False 21 21 md.thermal.isenthalpy=1 22 22 md=solve(md,TransientSolutionEnum()) -
issm/trunk/test/NightlyRun/test1501.m
r13975 r14067 13 13 14 14 %Solve for thinning rate -> -1 * surface mass balance 15 smb= 2.*ones(md.mesh.numberofvertices,1); 15 smb= 2.*ones(md.mesh.numberofvertices,1); 16 16 md.surfaceforcings.mass_balance= smb; 17 17 md.basalforcings.melting_rate= smb; -
issm/trunk/test/NightlyRun/test1502.m
r13975 r14067 14 14 15 15 %Solve for thinning rate -> -1 * surface mass balance 16 smb= 2.*ones(md.mesh.numberofvertices,1); 16 smb= 2.*ones(md.mesh.numberofvertices,1); 17 17 md.surfaceforcings.mass_balance= smb; 18 18 md.basalforcings.melting_rate= smb; -
issm/trunk/test/NightlyRun/test207.py
r13975 r14067 15 15 md=setflowequation(md,'macayeal','all') 16 16 md.cluster=generic('name',oshostname(),'np',3) 17 md.transient.isdiagnostic= 018 md.transient.isprognostic= 019 md.transient.isthermal= 120 md.transient.isgroundingline= 017 md.transient.isdiagnostic=False 18 md.transient.isprognostic=False 19 md.transient.isthermal=True 20 md.transient.isgroundingline=False 21 21 md=solve(md,TransientSolutionEnum()) 22 22 -
issm/trunk/test/NightlyRun/test228.py
r13975 r14067 24 24 25 25 md.surfaceforcings.mass_balance=numpy.vstack((smb,[1.5,3.])) 26 md.transient.isthermal= 026 md.transient.isthermal=False 27 27 28 28 md=solve(md,TransientSolutionEnum()) -
issm/trunk/test/NightlyRun/test229.py
r13975 r14067 24 24 25 25 md.surfaceforcings.mass_balance=numpy.vstack((smb,[1.5,3.])) 26 md.transient.isthermal= 026 md.transient.isthermal=False 27 27 28 28 md=solve(md,TransientSolutionEnum()) -
issm/trunk/test/NightlyRun/test230.py
r13975 r14067 25 25 26 26 md.surfaceforcings.mass_balance=numpy.vstack((smb,[1.5,3.])) 27 md.transient.isthermal= 027 md.transient.isthermal=False 28 28 29 29 md=solve(md,TransientSolutionEnum()) -
issm/trunk/test/NightlyRun/test231.py
r13975 r14067 25 25 26 26 md.surfaceforcings.mass_balance=numpy.vstack((smb,[1.5,3.])) 27 md.transient.isthermal= 027 md.transient.isthermal=False 28 28 29 29 md=solve(md,TransientSolutionEnum()) -
issm/trunk/test/NightlyRun/test232.py
r13975 r14067 18 18 md.timestepping.time_step=1. 19 19 md.timestepping.final_time=4. 20 md.transient.isdiagnostic= 021 md.transient.isprognostic= 022 md.transient.isthermal= 123 md.transient.isgroundingline= 020 md.transient.isdiagnostic=False 21 md.transient.isprognostic=False 22 md.transient.isthermal=True 23 md.transient.isgroundingline=False 24 24 md=solve(md,TransientSolutionEnum()) 25 25 -
issm/trunk/test/NightlyRun/test3009.py
r13975 r14067 14 14 md=setflowequation(md,'macayeal','all') 15 15 md.cluster=generic('name',oshostname(),'np',3) 16 md.transient.isdiagnostic= 017 md.transient.isprognostic= 018 md.transient.isthermal= 119 md.transient.isgroundingline= 016 md.transient.isdiagnostic=False 17 md.transient.isprognostic=False 18 md.transient.isthermal=True 19 md.transient.isgroundingline=False 20 20 md.autodiff.isautodiff=True 21 21 md=solve(md,TransientSolutionEnum()) -
issm/trunk/test/NightlyRun/test313.py
r13975 r14067 15 15 md.cluster=generic('name',oshostname(),'np',3) 16 16 md.verbose=verbose('convergence',True,'solution',True) 17 md.transient.isdiagnostic= 018 md.transient.isprognostic= 019 md.transient.isthermal= 120 md.transient.isgroundingline= 017 md.transient.isdiagnostic=False 18 md.transient.isprognostic=False 19 md.transient.isthermal=True 20 md.transient.isgroundingline=False 21 21 md=solve(md,TransientSolutionEnum()) 22 22 -
issm/trunk/test/NightlyRun/test326.py
r13975 r14067 16 16 md.cluster=generic('name',oshostname(),'np',3) 17 17 md.initialization.waterfraction=numpy.zeros((md.mesh.numberofvertices,1)) 18 md.transient.isdiagnostic= 019 md.transient.isprognostic= 020 md.transient.isthermal= 121 md.transient.isgroundingline= 018 md.transient.isdiagnostic=False 19 md.transient.isprognostic=False 20 md.transient.isthermal=True 21 md.transient.isgroundingline=False 22 22 md.thermal.isenthalpy=1 23 23 md=solve(md,TransientSolutionEnum()) -
issm/trunk/test/NightlyRun/test407.py
r13975 r14067 15 15 md=setflowequation(md,'pattyn','all') 16 16 md.cluster=generic('name',oshostname(),'np',3) 17 md.transient.isdiagnostic= 018 md.transient.isprognostic= 019 md.transient.isthermal= 120 md.transient.isgroundingline= 017 md.transient.isdiagnostic=False 18 md.transient.isprognostic=False 19 md.transient.isthermal=True 20 md.transient.isgroundingline=False 21 21 md=solve(md,TransientSolutionEnum()) 22 22 -
issm/trunk/test/NightlyRun/test423.py
r13975 r14067 29 29 md.cluster=generic('name',oshostname(),'np',3) 30 30 31 md.transient.isthermal= 032 md.transient.isprognostic= 033 md.transient.isdiagnostic= 034 md.transient.isgroundingline= 131 md.transient.isthermal=False 32 md.transient.isprognostic=False 33 md.transient.isdiagnostic=False 34 md.transient.isgroundingline=True 35 35 36 36 #test different grounding line dynamics. -
issm/trunk/test/NightlyRun/test424.py
r13975 r14067 20 20 md.geometry.surface=md.geometry.bed+md.geometry.thickness 21 21 md.surfaceforcings.mass_balance[:]=100. 22 md.transient.isdiagnostic= 023 md.transient.isgroundingline= 122 md.transient.isdiagnostic=False 23 md.transient.isgroundingline=True 24 24 md.groundingline.migration='AgressiveMigration' 25 25 -
issm/trunk/test/NightlyRun/test425.py
r13975 r14067 20 20 md.geometry.surface=md.geometry.bed+md.geometry.thickness 21 21 md.surfaceforcings.mass_balance[:]=-150. 22 md.transient.isdiagnostic= 023 md.transient.isgroundingline= 122 md.transient.isdiagnostic=False 23 md.transient.isgroundingline=True 24 24 md.groundingline.migration='SoftMigration' 25 25 -
issm/trunk/test/NightlyRun/test426.py
r13975 r14067 21 21 md.extrude(3,1.); 22 22 md=setflowequation(md,'macayeal','all'); 23 md.transient.isdiagnostic= 024 md.transient.isgroundingline= 123 md.transient.isdiagnostic=False 24 md.transient.isgroundingline=True 25 25 md.groundingline.migration='AgressiveMigration' 26 26 md.cluster=generic('name',oshostname(),'np',3) -
issm/trunk/test/NightlyRun/test427.py
r13975 r14067 22 22 23 23 md.surfaceforcings.mass_balance[:]=-150 24 md.transient.isdiagnostic= 025 md.transient.isgroundingline= 124 md.transient.isdiagnostic=False 25 md.transient.isgroundingline=True 26 26 md.groundingline.migration='SoftMigration' 27 27 md.cluster=generic('name',oshostname(),'np',3) -
issm/trunk/test/NightlyRun/test515.py
r13975 r14067 15 15 md.thermal.stabilization=2 16 16 md.cluster=generic('name',oshostname(),'np',3) 17 md.transient.isdiagnostic= 018 md.transient.isprognostic= 019 md.transient.isthermal= 120 md.transient.isgroundingline= 017 md.transient.isdiagnostic=False 18 md.transient.isprognostic=False 19 md.transient.isthermal=True 20 md.transient.isgroundingline=False 21 21 md=solve(md,TransientSolutionEnum()) 22 22
Note:
See TracChangeset
for help on using the changeset viewer.