Changeset 14067 for issm/trunk


Ignore:
Timestamp:
11/30/12 10:09:49 (12 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 14066

Location:
issm/trunk
Files:
16 deleted
102 edited
6 copied

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/aux-config

    • Property svn:ignore

      •  

        old new  
        1 
         1ar-lib
         2depcomp
         3compile
         4missing
         5config.guess
         6config.sub
         7ltmain.sh
         8install-sh

  • issm/trunk/configs/config-greenplanet.sh

    r13975 r14067  
    2222 --with-blacs-dir=$ISSM_DIR/externalpackages/petsc/install/ \
    2323 --with-graphics-lib=/usr/lib64/libX11.so \
    24  --with-cxxoptflags="-O3 -xS" \
     24 --with-cxxoptflags="-O3" \
    2525 --with-vendor=intel-linux

  • issm/trunk/configure.ac

    r13975 r14067  
    22
    33#AUTOCONF
    4 AC_INIT([ISSM],[4.2.3],[issm@jpl.nasa.gov],[issm],[http://issm.jpl.nasa.gov]) #Initializing configure
     4AC_INIT([ISSM],[4.2.4],[issm@jpl.nasa.gov],[issm],[http://issm.jpl.nasa.gov]) #Initializing configure
    55AC_CONFIG_AUX_DIR([./aux-config])         #Put config files in aux-config
    66AC_CONFIG_MACRO_DIR([m4])                 #m4 macros are located in m4

  • issm/trunk/etc/environment.sh

    r13975 r14067  
    8686
    8787DOXYGEN_DIR="$ISSM_DIR/externalpackages/doxygen/install"
    88 pathappend "$DOXYGEN_DIR/bin"
     88pathprepend "$DOXYGEN_DIR/bin"
    8989
    9090AUTOTOOLS_DIR="$ISSM_DIR/externalpackages/autotools/install"

  • issm/trunk/externalpackages/chaco

    • Property svn:ignore

      •  

        old new  
         1*.pdf
        12install
        23src
        34.ignore.txt
        45old
         6*.tar.gz

  • issm/trunk/externalpackages/doxygen/install.sh

    r13395 r14067  
    33
    44#Some cleanup
    5 rm -rf install
     5rm -rf install src
    66
    77#Download latest version
    8 svn co https://doxygen.svn.sourceforge.net/svnroot/doxygen/trunk install
     8svn co https://doxygen.svn.sourceforge.net/svnroot/doxygen/trunk src
    99
    1010#Configure doxygen

  • issm/trunk/externalpackages/gsl/install-android.sh

    r13975 r14067  
    2525    cd src
    2626
     27    mv ./../Makefile.am.patch ./
    2728    patch Makefile.am < Makefile.am.patch
    28    
     29
    2930    autoreconf -iv --force -I $ISSM_DIR/externalpackages/autotools/install/share/aclocal
    3031
    … …  
    3839#Compile gsl
    3940if [[ $step == "3" || $step == "0" ]]; then
    40         cd src
     41        cd $ISSM_DIR/externalpackages/gsl/src
     42
    4143    if [ $# -eq 0 ]; then
    42             make $j
     44            make
    4345    else
    4446            make -j $j

  • issm/trunk/externalpackages/triangle/configs/android/configure.make

    r13395 r14067  
    1010# http://www.codesourcery.com/gnu_toolchains/arm/arm_gnu_linux_abi.pdf
    1111
    12 source $ANDROID_DIR/android_aux.sh
    13 
    14 CC=${toolchain_path}-gcc
    15 AR=${toolchain_path}-ar
    16 RANLIB=${toolchain_path}-ranlib
     12CC=${host_triplet}-gcc
     13AR=${host_triplet}-ar
     14RANLIB=${host_triplet}-ranlib
    1715CSWITCHES = $(CFLAGS)
    1816TRILIBDEFS = -DTRILIBRARY

  • issm/trunk/externalpackages/triangle/install-android.sh

    r13395 r14067  
    11#!/bin/bash
    22set -eu
     3
     4source $ANDROID_DIR/android_aux.sh
    35
    46#Some cleanup

  • issm/trunk/src

  • issm/trunk/src/android/ISSM

    • Property svn:ignore

      •  

        old new  
        1 Makefile.in
        2 Makefile
         1gen
         2.settings
        33libs
        44obj
         5Makefile
         6Makefile.in

  • issm/trunk/src/android/ISSM/.classpath

    r13890 r14067  
    11<?xml version="1.0" encoding="UTF-8"?>
    22<classpath>
     3        <classpathentry kind="con" path="com.android.ide.eclipse.adt.ANDROID_FRAMEWORK"/>
     4        <classpathentry kind="con" path="com.android.ide.eclipse.adt.LIBRARIES"/>
    35        <classpathentry kind="src" path="src"/>
    46        <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"/>
    77        <classpathentry kind="output" path="bin/classes"/>
    88</classpath>

  • issm/trunk/src/android/ISSM/AndroidManifest.xml

    r13890 r14067  
    55
    66    <uses-sdk
    7         android:minSdkVersion="15"
     7        android:minSdkVersion="10"
    88        android:targetSdkVersion="15" />
    9 
     9        <uses-permission android:name="android.permission.WRITE_EXTERNAL_STORAGE" />
    1010    <application
    1111        android:icon="@drawable/ic_launcher"
    1212        android:label="@string/app_name"
    1313        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
    1522            android:name=".ISSM"
    1623            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>
    2224        </activity>
    2325    </application>

  • issm/trunk/src/android/ISSM/bin

    • Property svn:ignore

      •  

        old new  
        1 ISSM.apk
         1*ISSM.apk
        22resources.ap_
        33res

  • issm/trunk/src/android/ISSM/bin/AndroidManifest.xml

    r13891 r14067  
    55
    66    <uses-sdk
    7         android:minSdkVersion="15"
     7        android:minSdkVersion="10"
    88        android:targetSdkVersion="15" />
    9 
     9        <uses-permission android:name="android.permission.WRITE_EXTERNAL_STORAGE" />
    1010    <application
    1111        android:icon="@drawable/ic_launcher"
    1212        android:label="@string/app_name"
    1313        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
    1522            android:name=".ISSM"
    1623            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>
    2224        </activity>
    2325    </application>

  • issm/trunk/src/android/ISSM/gen/com/example/issm/R.java

    r13890 r14067  
    2222    public static final class id {
    2323        public static final int button1=0x7f080001;
     24        public static final int button2=0x7f080003;
    2425        public static final int input=0x7f080000;
    25         public static final int menu_settings=0x7f080003;
     26        public static final int menu_settings=0x7f080004;
    2627        public static final int output=0x7f080002;
    2728    }
    2829    public static final class layout {
    2930        public static final int activity_issm=0x7f030000;
     31        public static final int activity_mapselection=0x7f030001;
    3032    }
    3133    public static final class menu {

  • issm/trunk/src/android/ISSM/jni/Android.mk

    r13931 r14067  
    1111LOCAL_PATH := $(my_LOCAL_PATH)
    1212LOCAL_ALLOW_UNDEFINED_SYMBOLS := true
    13 LOCAL_CFLAGS := -Wno-psabi
     13LOCAL_CFLAGS := -Wno-psabi -DHAVE_CONFIG_H
    1414LOCAL_LDLIBS := -llog -ldl
    1515LOCAL_MODULE := IssmJni

  • issm/trunk/src/android/ISSM/jni/Main.cpp

    r13959 r14067  
    11#include <jni.h>
    22#include "../../../c/android/fac.h"
     3#include "../../../c/issm.h"
     4#include <cstddef>
     5#include <stdio.h>
     6///////////////////////////////////////////////////////////////////////////////////////////
    37namespace com_example_issm
    48{
    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)
    613        {
    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.
    1126        }
    1227//------------------------------------------------------------------------------------
    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}
    1551        };
    1652}
    17 
     53//------------------------------------------------------------------------------------
    1854using namespace com_example_issm;
    1955
    … …  
    2965        if(clazz)
    3066        {
    31                 env->RegisterNatives(clazz, method_table, 1);
     67                env->RegisterNatives(clazz, method_table, 3);
    3268                env->DeleteLocalRef(clazz);
    3369                return JNI_VERSION_1_6;
    … …  
    3571        else return -1;
    3672    }
    37 
    38     // Get jclass with env->FindClass.
    39     // Register methods with env->RegisterNatives.
    4073}
     74///////////////////////////////////////////////////////////////////////////////////////////

  • issm/trunk/src/android/ISSM/jni/issmlib/Android.mk

    r13931 r14067  
    22include $(CLEAR_VARS)
    33LOCAL_MODULE    := libISSMCore
    4 LOCAL_SRC_FILES := ../../../../c/libISSMCore.a
    5 LOCAL_EXPORT_C_INCLUDES := ../../../../c/android/
     4LOCAL_SRC_FILES := libISSMCore.a
     5LOCAL_EXPORT_C_INCLUDES := $(ISSM_DIR)
    66include $(PREBUILT_STATIC_LIBRARY)

  • issm/trunk/src/android/ISSM/src/com/example/issm/ISSM.java

    r13931 r14067  
    11package com.example.issm;
     2
     3import java.io.IOException;
     4import java.io.InputStream;
     5import java.nio.ByteBuffer;
     6import java.nio.ByteOrder;
     7import java.nio.DoubleBuffer;
    28
    39import android.os.Bundle;
    410import android.app.Activity;
     11import android.content.res.AssetManager;
    512import android.text.TextUtils;
    613import android.view.Menu;
    … …  
    1219
    1320
    14 public class ISSM extends Activity implements OnClickListener {
     21public class ISSM extends Activity implements OnClickListener
     22{
    1523        private EditText input;
    1624        private TextView output;
     25        private DoubleBuffer buff;
     26        private IssmJni issmNative;
     27        private String mapName;
     28        private String issmFolder;
     29        private int size;
    1730    @Override
     31  //------------------------------------------------------------------------------------------------   
    1832    public void onCreate(Bundle savedInstanceState) {
    1933        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        }
    2042        setContentView(R.layout.activity_issm);
    2143        this.input  = (EditText) super.findViewById(R.id.input);
    … …  
    2446        button.setOnClickListener(this);
    2547       
     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();
    2653    }
    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)
    2972        {
    30                 // TODO Auto-generated method stub
     73                //factorial method
    3174                String input = this.input.getText().toString();
    3275                if(TextUtils.isEmpty(input))
    3376                {
    3477                        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
    3888        }
    3989}

  • issm/trunk/src/android/ISSM/src/com/example/issm/IssmJni.java

    r13931 r14067  
    11package com.example.issm;
     2import java.nio.DoubleBuffer;
    23
    3 public class IssmJni
     4class IssmJni
    45{
    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        {
    711        System.loadLibrary("IssmJni");
    812    }
    9         public static long facIterative(long n)
    10         {
    11                 return fac(n);
    12         }
    1313}

  • issm/trunk/src/android/ISSM_Visual/bin

    • Property svn:ignore set to
      *.apk

  • issm/trunk/src/android/ISSM_Visual/src/com/example/issm_visual

    • Property svn:ignore set to
      *.java
      *.class

  • issm/trunk/src/c/Container/Observations.cpp

    r13975 r14067  
    424424        *pprediction = obs;
    425425}/*}}}*/
     426/*FUNCTION Observations::InterpolationV4{{{*/
     427void 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}/*}}}*/
    426510/*FUNCTION Observations::QuadtreeColoring{{{*/
    427511void Observations::QuadtreeColoring(IssmPDouble* A,IssmPDouble *x,IssmPDouble *y,int n){

  • issm/trunk/src/c/Container/Observations.h

    r13975 r14067  
    2828                void ClosestObservation(IssmDouble *px,IssmDouble *py,IssmDouble *pobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);
    2929                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);
    3031                void InterpolationKriging(IssmDouble *pprediction,IssmDouble *perror,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata,Variogram* variogram);
    3132                void InterpolationNearestNeighbor(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);

  • issm/trunk/src/c/Container/Options.h

    r13975 r14067  
    6464                        }
    6565                        else{
     66                                if(GetOption(name)) _printLine_("WARNING: option "<<name<<" found but fetched format not consistent, defaulting...");
    6667                                *pvalue=default_value;
    6768                        }

  • issm/trunk/src/c/android/fac.cpp

    r13709 r14067  
    1313
    1414/*}}}*/
    15 
     15#include "stdio.h"
    1616long fac::factorial(long n) {
    1717        long f = 1;
    1818        long i;
     19        printf("ok1\n");
    1920       
    2021        for(i = 1; i <= n; i++){
    2122                f *= i;
    2223        }
     24        printf("ok2\n");
    2325       
    2426        return f;

  • issm/trunk/src/c/classes/FemModel.cpp

    r13975 r14067  
    2020
    2121/*Object constructors and destructor*/
     22/*FUNCTION FemModel::FemModel(char* filename) {{{*/
     23FemModel::FemModel(char* filename){
     24
     25}
     26/*}}}*/
    2227/*FUNCTION FemModel::FemModel(int argc,char** argv){{{*/
    2328FemModel::FemModel(int argc,char** argv,COMM incomm){
    … …  
    502507
    503508                        /*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                        //}
    505512                        for(j=head_e[node->Sid()];j!=-1;j=next_e[j]){
    506513                                offset=count2offset_e[j];

  • issm/trunk/src/c/classes/FemModel.h

    r13975 r14067  
    4747
    4848                /*constructors, destructors: */
     49                FemModel(char* filename);
    4950                FemModel(int argc,char** argv,COMM comm_init);
    5051                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  
    6666
    6767        /*Build neighbors list*/
    68         if (xIsNan<IssmDouble>(iomodel->Data(MeshUpperelementsEnum)[index])) penta_elements_ids[1]=this->id; //upper penta is the same penta
     68        if (xIsNan<IssmDouble>(iomodel->Data(MeshUpperelementsEnum)[index]) || iomodel->Data(MeshUpperelementsEnum)[index]==-1.) penta_elements_ids[1]=this->id; //upper penta is the same penta
    6969        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 penta
     70        if (xIsNan<IssmDouble>(iomodel->Data(MeshLowerelementsEnum)[index]) || iomodel->Data(MeshLowerelementsEnum)[index]==-1.) penta_elements_ids[0]=this->id; //lower penta is the same penta
    7171        else                                    penta_elements_ids[0]=reCast<int,IssmDouble>((iomodel->Data(MeshLowerelementsEnum)[index]));
    7272        this->InitHookNeighbors(penta_elements_ids);

  • issm/trunk/src/c/modules/ElementConnectivityx/ElementConnectivityx.cpp

    r13975 r14067  
    1515#include "../../EnumDefinitions/EnumDefinitions.h"
    1616
    17 int hascommondedge(double* element1,double* element2);
     17int hascommondedge(int* element1,int* element2);
    1818
    19 void    ElementConnectivityx( double** pelementconnectivity, double* elements, int nel, double* nodeconnectivity, int nods, int width){
     19void ElementConnectivityx(int** pelementconnectivity,int* elements, int nels,int* nodeconnectivity, int nods, int width){
    2020
    2121        int i,j,k,n;
    … …  
    2323        /*intermediary: */
    2424        int    maxels;
    25         double element;
    26         double connectedelement;
     25        int  connectedelement;
    2726        int    connectedelementindex;
    2827        int    node;
    … …  
    3029        int    num_elements;
    3130
    32         /*output: */
    33         double* elementconnectivity=NULL;
    34 
    3531        /*maxels: */
    3632        maxels=width-1;
     33
    3734        /*Allocate connectivity: */
    38         elementconnectivity=xNewZeroInit<double>(nel*3);
     35        int* elementconnectivity=xNewZeroInit<int>(nels*3);
    3936
    4037        /*Go through all elements, and for each element, go through its nodes, to get the neighbouring elements.
    4138         * Once we get the neighbouring elements, figure out if they share a segment with the current element. If so,
    4239         * plug them in the connectivity, unless they are already there.: */
     40        for(n=0;n<nels;n++){
    4341
    44         for(n=0;n<nel;n++){
    45 
    46                 element=(double)(n+1); //matlab indexing
     42                //element=n+1; //matlab indexing
    4743
    4844                for(i=0;i<3;i++){
    4945
    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.
    5147                        index=node-1;
    5248
    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.
    5450
    5551                        for(j=0;j<num_elements;j++){
    5652
    5753                                /*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.
    6056
    61                                 if(hascommondedge(elements+n*3+0,elements+connectedelementindex*3+0)){
     57                                if(hascommondedge(&elements[n*3+0],&elements[connectedelementindex*3+0])){
    6258                                        /*Ok, this connected element has a commond edge  with element, plug it into elementconnectivity, unless
    6359                                         *it is already there: */
    6460
    6561                                        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;
    6864                                                        break;
    6965                                                }
    7066                                                else{
    71                                                         if(connectedelement==*(elementconnectivity+3*n+k))break;
     67                                                        if(connectedelement==elementconnectivity[3*n+k]) break;
    7268                                                }
    7369                                        }
    … …  
    8177}
    8278
    83 int hascommondedge(double* element1,double* element2){
     79int hascommondedge(int* el1,int* el2){
    8480
    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++;
    9285                }
    9386        }
    94         if (count==2)return 1;
    95         else return 0;
     87        if(count==2)
     88         return 1;
     89        else
     90         return 0;
    9691}

  • issm/trunk/src/c/modules/ElementConnectivityx/ElementConnectivityx.h

    r13975 r14067  
    77
    88/* local prototypes: */
    9 void    ElementConnectivityx( double** pelementconnectivity, double* elements, int nel, double* nodeconnectivity, int nods, int width);
     9void    ElementConnectivityx(int** pelementconnectivity,int* elements,int nels,int* nodeconnectivity, int nods, int width);
    1010
    1111#endif  /* _ELEMENTCONNECTIVITYX_H */

  • issm/trunk/src/c/modules/Krigingx/Krigingx.cpp

    r13395 r14067  
    2222        /*Intermediaries*/
    2323        int           mindata,maxdata;
     24        double        dmindata,dmaxdata,dnumthreads; //FIXME (Options come as double but we want to retrive integers)
    2425        double        radius;
    2526        char         *output       = NULL;
    … …  
    3738        ProcessVariogram(&variogram,options);
    3839        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)
    4143
    4244        /*Process observation dataset*/
    … …  
    8890                gate.predictions  = predictions;
    8991                gate.error        = error;
    90                 gate.percent      = xNewZeroInit<double>(num);
     92                gate.numdone      = xNewZeroInit<int>(num);
    9193
    9294                /*launch the thread manager with Krigingxt as a core: */
    9395                LaunchThread(NearestNeighbort,(void*)&gate,num);
    9496                _printLine_("\r      interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%");
    95                 xDelete<double>(gate.percent);
     97                xDelete<int>(gate.numdone);
    9698        }
    9799        else if(strcmp(output,"idw")==0){ //Inverse distance weighting
    … …  
    109111                gate.predictions  = predictions;
    110112                gate.error        = error;
    111                 gate.percent      = xNewZeroInit<double>(num);
     113                gate.numdone      = xNewZeroInit<int>(num);
    112114                gate.power        = power;
    113115
    … …  
    115117                LaunchThread(idwt,(void*)&gate,num);
    116118                _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
    121125                /*initialize thread parameters: */
    122126                gate.n_interp     = n_interp;
    … …  
    130134                gate.predictions  = predictions;
    131135                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);
    133160
    134161                /*launch the thread manager with Krigingxt as a core: */
    135162                LaunchThread(Krigingxt,(void*)&gate,num);
    136163                _printLine_("\r      interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%");
    137                 xDelete<double>(gate.percent);
     164                xDelete<int>(gate.numdone);
    138165        }
    139166        else{
    … …  
    176203        double       *predictions  = gate->predictions;
    177204        double       *error        = gate->error;
    178         double       *percent      = gate->percent;
    179 
    180         /*Intermediaries*/
    181         double        localpercent;
     205        int          *numdone      = gate->numdone;
    182206
    183207        /*partition loop across threads: */
    … …  
    186210
    187211                /*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                }
    192218
    193219                /*Kriging interpolation*/
    … …  
    224250        double       *predictions  = gate->predictions;
    225251        double       *error        = gate->error;
    226         double       *percent      = gate->percent;
    227 
    228         /*Intermediaries*/
    229         int           i;
    230         double        localpercent;
     252        int          *numdone      = gate->numdone;
    231253
    232254        /*partition loop across threads: */
    … …  
    235257
    236258                /*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                }
    241265
    242266                observations->InterpolationNearestNeighbor(&predictions[idx],x_interp[idx],y_interp[idx],radius);
    … …  
    272296        double       *predictions  = gate->predictions;
    273297        double       *error        = gate->error;
    274         double       *percent      = gate->percent;
     298        int          *numdone      = gate->numdone;
    275299        double        power        = gate->power;
    276 
    277         /*Intermediaries*/
    278         double localpercent;
    279300
    280301        /*partition loop across threads: */
    … …  
    283304
    284305                /*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                }
    289312
    290313                observations->InterpolationIDW(&predictions[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata,power);
     314        }
     315        return NULL;
     316}/*}}}*/
     317/*FUNCTION v4t{{{*/
     318void* 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);
    291359        }
    292360        return NULL;

  • issm/trunk/src/c/modules/Krigingx/Krigingx.h

    r13395 r14067  
    2828        double       *predictions;
    2929        double       *error;
    30         double       *percent;
     30        int          *numdone;
    3131        double        power;//for idw
    3232}KrigingxThreadStruct;
    … …  
    3535void* NearestNeighbort(void*);
    3636void* idwt(void*);
     37void* v4t(void*);
    3738#endif /* _KRIGINGX_H */

  • issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.cpp

    r13975 r14067  
    1818#include "../../EnumDefinitions/EnumDefinitions.h"
    1919
    20 void    NodeConnectivityx( double** pconnectivity, int* pwidth, double* elements, int nel, int nods){
     20void    NodeConnectivityx(int** pconnectivity,int* pwidth,int* elements, int nels, int nods){
    2121
    2222        int i,j,n;
    … …  
    2929        int     num_elements;
    3030        int     already_plugged=0;
    31         double  element;
    32 
    33         /*output: */
    34         double* connectivity=NULL;
     31        int element;
    3532
    3633        /*Allocate connectivity: */
    37         connectivity=xNewZeroInit<double>(nods*width);
     34        int* connectivity=xNewZeroInit<int>(nods*width);
    3835
    3936        /*Go through all elements, and for each elements, plug into the connectivity, all the nodes.
    4037         * If nodes are already plugged into the connectivity, skip them.: */
    41         for(n=0;n<nel;n++){
     38        for(n=0;n<nels;n++){
    4239
    43                 element=(double)(n+1); //matlab indexing
     40                element=n+1; //matlab indexing
    4441
    4542                for(i=0;i<3;i++){
    4643
    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.
    4845                        index=node-1;
    4946
    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.
    5148
    5249                        already_plugged=0;
    … …  
    6057
    6158                        /*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;
    6461
    6562                }
    … …  
    6966         * warn the user to increase the connectivity width: */
    7067        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");
    7270        }
    7371

  • issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.h

    r13975 r14067  
    77
    88/* local prototypes: */
    9 void    NodeConnectivityx( double** pconnectivity, int* pwidth,double* elements, int nel, int nods);
     9void    NodeConnectivityx(int** pconnectivity,int* pwidth,int* elements,int nels, int nods);
    1010
    1111#endif  /* _NODECONNECTIVITYX_H */

  • issm/trunk/src/c/shared/Threads/LaunchThread.cpp

    r13975 r14067  
    3131        pthread_t      *threads = NULL;
    3232        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)");
    3337
    3438        /*dynamically allocate: */

  • issm/trunk/src/dox/issm.dox

    r13975 r14067  
    4646</th>
    4747<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>
    4949</tr>
    5050<tr>
    51 <th  bgcolor=#C6E2FF style="text-align:left;"> MATLAB </th><td  bgcolor=#C6E2FF style="text-align:right;">913</td><td  bgcolor=#C6E2FF style="text-align:right;">6625</td><td  bgcolor=#C6E2FF style="text-align:right;">13054</td><td  bgcolor=#C6E2FF style="text-align:right;">30160</td><td  bgcolor=#C6E2FF style="text-align:right;">49839</td>
     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>
    5252</tr>
    5353<tr>
    54 <th  bgcolor=#FFFFFF style="text-align:left;"> C/C++  Header </th><td  bgcolor=#FFFFFF style="text-align:right;">388</td><td  bgcolor=#FFFFFF style="text-align:right;">3069</td><td  bgcolor=#FFFFFF style="text-align:right;">2983</td><td  bgcolor=#FFFFFF style="text-align:right;">11855</td><td  bgcolor=#FFFFFF style="text-align:right;">17907</td>
     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>
    5555</tr>
    5656<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>
    5858</tr>
    5959<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>
    6161</tr>
    6262<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;">98</td><td  bgcolor=#C6E2FF style="text-align:right;">0</td><td  bgcolor=#C6E2FF style="text-align:right;">370</td><td  bgcolor=#C6E2FF style="text-align:right;">468</td>
     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>
    6464</tr>
    6565<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;">58</td><td  bgcolor=#FFFFFF style="text-align:right;">81</td><td  bgcolor=#FFFFFF style="text-align:right;">268</td><td  bgcolor=#FFFFFF style="text-align:right;">407</td>
     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>
    6767</tr>
    6868<tr>
    … …  
    7676</tr>
    7777<tr>
    78 <th  bgcolor=#FFFFFF style="text-align:left;"> SUM: </th><td  bgcolor=#FFFFFF style="text-align:right;">1912</td><td  bgcolor=#FFFFFF style="text-align:right;">27936</td><td  bgcolor=#FFFFFF style="text-align:right;">36727</td><td  bgcolor=#FFFFFF style="text-align:right;">113707</td><td  bgcolor=#FFFFFF style="text-align:right;">178370</td>
     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>
    7979</tr>
    8080</table>

  • issm/trunk/src/m/boundaryconditions/SetIceShelfBC.py

    r13975 r14067  
    2626                        raise IOError("SetIceShelfBC error message: ice front file '%s' not found." % icefrontfile)
    2727                [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))
    2929        else:
    30                 nodeonicefront=numpy.zeros((md.mesh.numberofvertices))
     30                nodeonicefront=numpy.zeros((md.mesh.numberofvertices),bool)
    3131
    3232#       pos=find(md.mesh.vertexonboundary & ~nodeonicefront);
    … …  
    5656        #segment on Neumann (Ice Front)
    5757#       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]
    5959        if   md.mesh.dimension==2:
    6060                pressureload=md.mesh.segments[pos,:]
    … …  
    6262#               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)];
    6363                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)
    6565                for i in xrange(1,md.mesh.numberoflayers):
    6666#                       pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ];
    … …  
    6969        #Add water or air enum depending on the element
    7070#       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)))
    7272
    7373        #plug onto model

  • issm/trunk/src/m/boundaryconditions/SetMarineIceSheetBC.py

    r13975 r14067  
    2828                        raise IOError("SetMarineIceSheetBC error message: ice front file '%s' not found." % icefrontfile)
    2929                [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))
    3131        else:
    3232                #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]=1
    35                 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)
    3636
    3737#       pos=find(md.mesh.vertexonboundary & ~vertexonicefront);
    … …  
    6262        #segment on Neumann (Ice Front)
    6363#       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]
    6565        if   md.mesh.dimension==2:
    6666                pressureload=md.mesh.segments[pos,:]
    … …  
    6868#               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)];
    6969                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)
    7171                for i in xrange(1,md.mesh.numberoflayers):
    7272#                       pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ];
    … …  
    7575        #Add water or air enum depending on the element
    7676#       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)))
    7878
    7979        #plug onto model

  • issm/trunk/src/m/classes/clusters/greenplanet.m

    r13395 r14067  
    1414                 cpuspernode=8;
    1515                 port=8000;
    16                  queue='rignot';
     16                 queue='c6145';
    1717                 codepath='';
    1818                 executionpath='';
    … …  
    5050                 function md = checkconsistency(cluster,md,solution,analyses) % {{{
    5151
    52                          available_queues={'rignot','default'};
     52                         available_queues={'c6145','default'};
    5353                         queue_requirements_time=[Inf Inf];
    5454                         queue_requirements_np=[80 80];

  • issm/trunk/src/m/classes/flowequation.py

    r13975 r14067  
    1919                # {{{ Properties
    2020               
    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
    2525                self.vertex_equation      = float('NaN')
    2626                self.element_equation     = float('NaN')
    … …  
    3737                string='   flow equation parameters:'
    3838
    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)"))
    4848                return string
    4949                #}}}

  • issm/trunk/src/m/classes/initialization.py

    r13975 r14067  
    7070                        md = checkfield(md,'initialization.vy','NaN',1,'size',[md.mesh.numberofvertices])
    7171                        #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)):
    7474                                md.checkmessage("at least one triangle has all its vertices with a zero velocity")
    7575                if ThermalAnalysisEnum() in analyses:

  • issm/trunk/src/m/classes/mask.py

    r13395 r14067  
    3333                string="%s\n%s"%(string,fielddisplay(self,"elementonfloatingice","element on floating ice flags list"))
    3434                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  list"))
     35                string="%s\n%s"%(string,fielddisplay(self,"elementongroundedice","element on grounded ice list"))
    3636                string="%s\n%s"%(string,fielddisplay(self,"vertexongroundedice","vertex on grounded ice flags list"))
    3737                string="%s\n%s"%(string,fielddisplay(self,"elementonwater","element on water flags list"))

  • issm/trunk/src/m/classes/mesh.py

    r13975 r14067  
    9292                string="%s\n%s"%(string,fielddisplay(self,"vertexonsurface","upper vertices flags list"))
    9393                string="%s\n%s"%(string,fielddisplay(self,"elementonsurface","upper elements flags list"))
    94                 string="%s\n%s"%(string,fielddisplay(self,"uppervertex","upper vertex list (NaN for vertex on the upper surface)"))
    95                 string="%s\n%s"%(string,fielddisplay(self,"upperelements","upper element list (NaN for element on the upper layer)"))
    96                 string="%s\n%s"%(string,fielddisplay(self,"lowervertex","lower vertex list (NaN for 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)"))
    9898                string="%s\n%s"%(string,fielddisplay(self,"vertexonboundary","vertices on the boundary of the domain flag list"))
    9999                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  
    22import numpy
    33import copy
     4import sys
    45from mesh import mesh
    56from mask import mask
    … …  
    207208                #kick out all elements with 3 dirichlets
    208209                spc_elem=numpy.nonzero(numpy.logical_not(flag_elem))[0]
    209                 spc_node=numpy.unique(md1.mesh.elements[spc_elem,:]).astype(int)-1
     210                spc_node=numpy.unique(md1.mesh.elements[spc_elem,:])-1
    210211                flag=numpy.ones(md1.mesh.numberofvertices)
    211212                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]
    213214                flag_elem[pos]=0
    214215
    215216                #extracted elements and nodes lists
    216217                pos_elem=numpy.nonzero(flag_elem)[0]
    217                 pos_node=numpy.unique(md1.mesh.elements[pos_elem,:]).astype(int)-1
     218                pos_node=numpy.unique(md1.mesh.elements[pos_elem,:])-1
    218219
    219220                #keep track of some fields
    … …  
    226227
    227228                #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)
    229230                Pelem[pos_elem]=numpy.arange(1,numberofelements2+1)
    230                 Pnode=numpy.zeros(numberofvertices1)
     231                Pnode=numpy.zeros(numberofvertices1,int)
    231232                Pnode[pos_node]=numpy.arange(1,numberofvertices2+1)
    232233
    … …  
    234235                elements_1=copy.deepcopy(md1.mesh.elements)
    235236                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]
    239240                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]
    243244
    244245                #OK, now create the new model!
    … …  
    291292                if md1.mesh.dimension==3:
    292293                        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]
    295296
    296297                        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]
    299300
    300301                        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]
    303304
    304305                        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]
    307308
    308309                #Initial 2d mesh
    … …  
    316317                        md2.mesh.numberofvertices2d=numpy.size(pos_node_2d)
    317318                        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]
    321322
    322323                        md2.mesh.x2d=md1.mesh.x[pos_node_2d]
    … …  
    324325
    325326                #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...
    327328                        #renumber first two columns
    328329                        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]
    333334                        #remove edges when the 2 vertices are not in the domain.
    334335                        md2.mesh.edges=md2.mesh.edges[numpy.nonzero(numpy.logical_and(md2.mesh.edges[:,0],md2.mesh.edges[:,1]))[0],:]
    … …  
    364365                        [md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity)
    365366                        md2.mesh.segments=contourenvelope(md2)
    366                         md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2)
    367                         md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2].astype(int)-1]=1
     367                        md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2,bool)
     368                        md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2]-1]=True
    368369                else:
    369370                        #First do the connectivity for the contourenvelope in 2d
    … …  
    371372                        [md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity)
    372373                        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]=1
     374                        md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2/md2.mesh.numberoflayers,bool)
     375                        md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2]-1]=True
    375376                        md2.mesh.vertexonboundary=numpy.tile(md2.mesh.vertexonboundary,md2.mesh.numberoflayers)
    376377                        #Then do it for 3d as usual
    … …  
    381382                #Catch the elements that have not been extracted
    382383                orphans_elem=numpy.nonzero(numpy.logical_not(flag_elem))[0]
    383                 orphans_node=numpy.unique(md1.mesh.elements[orphans_elem,:]).astype(int)-1
     384                orphans_node=numpy.unique(md1.mesh.elements[orphans_elem,:])-1
    384385                #Figure out which node are on the boundary between md2 and md1
    385386                nodestoflag1=numpy.intersect1d(orphans_node,pos_node)
    … …  
    443444
    444445                #Keep track of pos_node and pos_elem
    445                 md2.mesh.extractedvertices=pos_node.astype(float)+1
    446                 md2.mesh.extractedelements=pos_elem.astype(float)+1
     446                md2.mesh.extractedvertices=pos_node+1
     447                md2.mesh.extractedelements=pos_elem+1
    447448
    448449                return md2
    … …  
    526527
    527528                #Extrude elements
    528                 elements3d=numpy.empty((0,6))
     529                elements3d=numpy.empty((0,6),int)
    529530                for i in xrange(numlayers-1):
    530531                        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
    … …  
    532533
    533534                #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)
    536537                mesh.lowervertex[md.mesh.numberofvertices:]=numpy.arange(1,(numlayers-1)*md.mesh.numberofvertices+1)
    537538                mesh.uppervertex[:(numlayers-1)*md.mesh.numberofvertices]=numpy.arange(md.mesh.numberofvertices+1,number_nodes3d+1)
    … …  
    540541
    541542                #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)
    544545                mesh.lowerelements[md.mesh.numberofelements:]=numpy.arange(1,(numlayers-2)*md.mesh.numberofelements+1)
    545546                mesh.upperelements[:(numlayers-2)*md.mesh.numberofelements]=numpy.arange(md.mesh.numberofelements+1,(numlayers-1)*md.mesh.numberofelements+1)
    … …  
    603604
    604605                #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)
    609610
    610611                #elementstype
    611612                if not numpy.any(numpy.isnan(md.flowequation.element_equation)):
    612613                        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)
    614615                        md.flowequation.element_equation=project3d(md,'vector',oldelements_type,'type','element')
    615616
    … …  
    617618                if not numpy.any(numpy.isnan(md.flowequation.vertex_equation)):
    618619                        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)
    620621                        md.flowequation.vertex_equation=project3d(md,'vector',oldvertices_type,'type','node')
    621622
    … …  
    635636                #in 3d, pressureload: [node1 node2 node3 node4 element]
    636637                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)
    638639                for i in xrange(numlayers-1):
    639640                        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]))))
    … …  
    642643                #connectivity
    643644                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
    645646                for i in xrange(1,numlayers-1):
    646647                        md.mesh.elementconnectivity[i*md.mesh.numberofelements2d:(i+1)*md.mesh.numberofelements2d,:] \
    647648                                =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))]=0
     649                md.mesh.elementconnectivity[numpy.nonzero(md.mesh.elementconnectivity<0)]=0
    649650
    650651                #materials

  • issm/trunk/src/m/classes/transient.m

    r13395 r14067  
    4646
    4747                        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');
    4849                        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');
    5050                        fielddisplay(obj,'isgroundingline','indicates if a groundingline migration is used in the transient');
    5151                        fielddisplay(obj,'requested_outputs','list of additional outputs requested');

  • issm/trunk/src/m/classes/transient.py

    r13395 r14067  
    1616        def __init__(self):
    1717                # {{{ Properties
    18                 self.isprognostic      = 0
    19                 self.isdiagnostic      = 0
    20                 self.isthermal         = 0
    21                 self.isgroundingline   = 0
     18                self.isprognostic      = False
     19                self.isdiagnostic      = False
     20                self.isthermal         = False
     21                self.isgroundingline   = False
    2222                self.requested_outputs = float('NaN')
    2323
    … …  
    3030                string='   transient solution parameters:'
    3131                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'))
    3233                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'))
    3434                string="%s\n%s"%(string,fielddisplay(self,'isgroundingline','indicates if a groundingline migration is used in the transient'))
    3535                string="%s\n%s"%(string,fielddisplay(self,'requested_outputs','list of additional outputs requested'))
    … …  
    4141               
    4242                #full analysis: Diagnostic, Prognostic and Thermal but no groundingline migration for now
    43                 self.isprognostic=1
    44                 self.isdiagnostic=1
    45                 self.isthermal=1
    46                 self.isgroundingline=0
     43                self.isprognostic=True
     44                self.isdiagnostic=True
     45                self.isthermal=True
     46                self.isgroundingline=False
    4747
    4848                return self

  • issm/trunk/src/m/consistency/checkfield.py

    r13395 r14067  
    143143        if options.getfieldvalue('forcing',0):
    144144                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:
    146146                                md = md.checkmessage(options.getfieldvalue('message',\
    147147                                        "field '%s' should have only one column as there are md.mesh.numberofvertices lines" % fieldname))

  • issm/trunk/src/m/exp/expcoarsen.m

    r13975 r14067  
    2525
    2626%Get exp oldfile
    27 [path root ext ver]=fileparts(oldfile);
     27[path root ext]=fileparts(oldfile);
    2828A=expread(oldfile);
    2929numprofiles=size(A,2);

  • issm/trunk/src/m/extrusion/project3d.py

    r13395 r14067  
    3838
    3939        vector1d=False
    40         if isinstance(vector2d,numpy.ndarray) and len(vector2d.shape)==1:
     40        if isinstance(vector2d,numpy.ndarray) and numpy.ndim(vector2d)==1:
    4141                vector1d=True
    4242                vector2d=vector2d.reshape(numpy.size(vector2d),1)
    … …  
    4949                #Initialize 3d vector
    5050                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)
    5252                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)
    5454                        projected_vector[-1,:]=vector2d[-1,:]
    5555                        vector2d=vector2d[:-1,:]
    … …  
    6868                #Initialize 3d vector
    6969                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)
    7171                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)
    7373                        projected_vector[-1,:]=vector2d[-1,:]
    7474                        vector2d=vector2d[:-1,:]

  • issm/trunk/src/m/geometry/FlagElements.py

    r13975 r14067  
    2424        if   isinstance(region,(str,unicode)):
    2525                if   not region:
    26                         flag=numpy.zeros(md.mesh.numberofelements,'bool')
     26                        flag=numpy.zeros(md.mesh.numberofelements,bool)
    2727                        invert=0
    2828                elif strcmpi(region,'all'):
    29                         flag=numpy.ones(md.mesh.numberofelements,'bool')
     29                        flag=numpy.ones(md.mesh.numberofelements,bool)
    3030                        invert=0
    3131                else:
    … …  
    4343                                raise RuntimeError("FlagElements.py calling basinzoom.py is not complete.")
    4444                                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)
    4747                        else:
    4848                                #ok, flag elements
    4949                                [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)
    5051
    5152                if invert:

  • issm/trunk/src/m/geometry/GetAreas.py

    r13975 r14067  
    3333        #initialization
    3434        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]
    4141
    4242        #compute the volume of each element
    … …  
    4646        else:
    4747                #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])
    4949                areas=(0.5*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)))*thickness
    5050

  • issm/trunk/src/m/mesh/ComputeHessian.py

    r13975 r14067  
    4545
    4646        #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)
    4949
    5050        #Compute gradient for each node (average of the elements around)
    … …  
    5555
    5656        #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)))
    5858
    5959        if strcmpi(type,'node'):

  • issm/trunk/src/m/mesh/GetNodalFunctionsCoeff.py

    r13975 r14067  
    3939
    4040        #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]
    4747        invdet=1./(x1*(y2-y3)-x2*(y1-y3)+x3*(y1-y2))
    4848

  • issm/trunk/src/m/mesh/bamg.py

    r13975 r14067  
    320320        md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
    321321        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)
    326326
    327327        #Fill in rest of fields:
    … …  
    331331        md.mesh.numberofedges=numpy.size(md.mesh.edges,axis=0)
    332332        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]=1
     333        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
    340340        md.mesh.elementconnectivity=md.private.bamg['mesh'].ElementConnectivity
    341341        md.mesh.elementconnectivity[numpy.nonzero(numpy.isnan(md.mesh.elementconnectivity))]=0
     342        md.mesh.elementconnectivity=md.mesh.elementconnectivity.astype(int)
    342343
    343344        #Check for orphan

  • issm/trunk/src/m/mesh/meshconvert.py

    r13975 r14067  
    3535        md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
    3636        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)
    4141
    4242        #Fill in rest of fields:
    … …  
    4646        md.mesh.numberofedges=numpy.size(md.mesh.edges,axis=0)
    4747        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]=1
     48        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
    5555
    5656        return md

  • issm/trunk/src/m/mesh/rifts/meshprocessoutsiderifts.py

    r13975 r14067  
    4141                        B=node_connected_to_tip
    4242
    43                         elements=numpy.empty(0)
     43                        elements=numpy.empty(0,int)
    4444
    4545                        while flags(B):    #as long as B does not belong to the domain outline, keep looking.
    … …  
    6262               
    6363                        #replace tip in elements
    64                         newelements=md.mesh.elements[elements.astype(int)-1,:]
     64                        newelements=md.mesh.elements[elements-1,:]
    6565                        pos=numpy.nonzero(newelements==tip)
    6666                        newelements[pos]=num
    67                         md.mesh.elements[elements.astype(int)-1,:]=newelements
     67                        md.mesh.elements[elements-1,:]=newelements
    6868                        rift.tips=numpy.concatenate((rift.tips,num))
    6969
    … …  
    8181        md.mesh.numberofvertices=numpy.size(md.mesh.x)
    8282        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]=1
     83        md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x),bool)
     84        md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
    8585        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)
    9191
    9292        return md

  • issm/trunk/src/m/mesh/rifts/meshprocessrifts.py

    r13975 r14067  
    2323        #Call MEX file
    2424        [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)
    2526        md.mesh.x=md.mesh.x.reshape(-1)
    2627        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)
    2730        if not isinstance(md.rifts.riftstruct,list) or not md.rifts.riftstruct:
    2831                raise RuntimeError("TriMeshProcessRifts did not find any rift")
    … …  
    3336        md.mesh.numberofvertices=numpy.size(md.mesh.x)
    3437        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]=1
    37         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)
    4144
    4245        #get coordinates of rift tips

  • issm/trunk/src/m/mesh/triangle.py

    r13975 r14067  
    4444        #Mesh using TriMesh
    4545        [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)
    4649
    4750        #Fill in rest of fields:
    … …  
    4952        md.mesh.numberofvertices = numpy.size(md.mesh.x)
    5053        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)
    5760
    5861        #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)
    6063        [md.mesh.elementconnectivity] = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)
    6164

  • issm/trunk/src/m/parameterization/contourenvelope.py

    r13975 r14067  
    7676                        pos=numpy.nonzero(flags)
    7777                        elemin[pos]=1
    78                         nodein[mesh.elements[pos,:].astype(int)-1]=1
     78                        nodein[mesh.elements[pos,:]-1]=1
    7979
    8080                        #modify element connectivity
    … …  
    8888        if len(args)==1:
    8989                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)
    9191
    9292        #Find segments on boundary
    9393        pos=numpy.nonzero(elementonboundary)[0]
    9494        num_segments=numpy.size(pos)
    95         segments=numpy.zeros((num_segments*3,3))
     95        segments=numpy.zeros((num_segments*3,3),int)
    9696        count=0
    9797
    9898        for el1 in pos:
    99                 els2=mesh.elementconnectivity[el1,numpy.nonzero(mesh.elementconnectivity[el1,:])[0]].astype(int)-1
     99                els2=mesh.elementconnectivity[el1,numpy.nonzero(mesh.elementconnectivity[el1,:])[0]]-1
    100100                if numpy.size(els2)>1:
    101101                        flag=numpy.intersect1d(mesh.elements[els2[0],:],mesh.elements[els2[1],:])

  • issm/trunk/src/m/parameterization/setflowequation.py

    r13975 r14067  
    5050        #Flag the elements that have not been flagged as filltype
    5151        if   strcmpi(filltype,'hutter'):
    52                 hutterflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(macayealflag,pattynflag)))]=1
     52                hutterflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(macayealflag,pattynflag)))]=True
    5353        elif strcmpi(filltype,'macayeal'):
    54                 macayealflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(pattynflag,stokesflag))))]=1
     54                macayealflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(pattynflag,stokesflag))))]=True
    5555        elif strcmpi(filltype,'pattyn'):
    56                 pattynflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(macayealflag,stokesflag))))]=1
     56                pattynflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(macayealflag,stokesflag))))]=True
    5757
    5858        #check that each element has at least one flag
    … …  
    6363        if any(hutterflag+macayealflag+l1l2flag+pattynflag+stokesflag>1):
    6464                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))]=0
    66                 hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,pattynflag))]=0
    67                 macayealflag[numpy.nonzero(numpy.logical_and(macayealflag,pattynflag))]=0
     65                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
    6868
    6969        #Check that no pattyn or stokes for 2d mesh
    … …  
    7777
    7878        #Initialize node fields
    79         nodeonhutter=numpy.zeros(md.mesh.numberofvertices)
    80         nodeonhutter[md.mesh.elements[numpy.nonzero(hutterflag),:].astype(int)-1]=1
    81         nodeonmacayeal=numpy.zeros(md.mesh.numberofvertices)
    82         nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1
    83         nodeonl1l2=numpy.zeros(md.mesh.numberofvertices)
    84         nodeonl1l2[md.mesh.elements[numpy.nonzero(l1l2flag),:].astype(int)-1]=1
    85         nodeonpattyn=numpy.zeros(md.mesh.numberofvertices)
    86         nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1
    87         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)
    8989
    9090        #First modify stokesflag to get rid of elements contrained everywhere (spc + border with pattyn or macayeal)
    … …  
    9696                                              numpy.logical_and(nodeonpattyn,nodeonstokes).reshape(-1,1)).astype(int)    #find all the nodes on the boundary of the domain without icefront
    9797#               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 icefront
    99                 stokesflag[numpy.nonzero(fullspcelems.reshape(-1))]=0
    100                 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1
     98                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
    101101
    102102        #Then complete with NoneApproximation or the other model used if there is no stokes
    103103        if any(stokesflag):
    104104                if   any(pattynflag):    #fill with pattyn
    105                         pattynflag[numpy.logical_not(stokesflag)]=1
    106                         nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1
     105                        pattynflag[numpy.logical_not(stokesflag)]=True
     106                        nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True
    107107                elif any(macayealflag):    #fill with macayeal
    108                         macayealflag[numpy.logical_not(stokesflag)]=1
    109                         nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1
     108                        macayealflag[numpy.logical_not(stokesflag)]=True
     109                        nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True
    110110                else:    #fill with none
    111                         noneflag[numpy.nonzero(numpy.logical_not(stokesflag))]=1
     111                        noneflag[numpy.nonzero(numpy.logical_not(stokesflag))]=True
    112112
    113113        #Now take care of the coupling between MacAyeal and Pattyn
    114114        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)
    121121        if   strcmpi(coupling_method,'penalties'):
    122122                #Create the border nodes between Pattyn and MacAyeal and extrude them
    … …  
    135135                if   any(macayealflag) and any(pattynflag):    #coupling macayeal pattyn
    136136                        #Find node at the border
    137                         nodeonmacayealpattyn[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonpattyn))]=1
     137                        nodeonmacayealpattyn[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonpattyn))]=True
    138138                        #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])
    140140                        commonelements=numpy.sum(matrixelements,axis=1)!=0
    141                         commonelements[numpy.nonzero(pattynflag)]=0    #only one layer: the elements previously in macayeal
    142                         macayealflag[numpy.nonzero(commonelements)]=0    #these elements are now macayealpattynelements
    143                         macayealpattynflag[numpy.nonzero(commonelements)]=1
    144                         nodeonmacayeal[:]=0
    145                         nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1
     141                        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
    146146
    147147                        #rule out elements that don't touch the 2 boundaries
    148148                        pos=numpy.nonzero(macayealpattynflag)[0]
    149149                        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)
    152152                        pos1=numpy.nonzero(elist==1)[0]
    153                         macayealflag[pos[pos1]]=1
    154                         macayealpattynflag[pos[pos1]]=0
     153                        macayealflag[pos[pos1]]=True
     154                        macayealpattynflag[pos[pos1]]=False
    155155                        pos2=numpy.nonzero(elist==-1)[0]
    156                         pattynflag[pos[pos2]]=1
    157                         macayealpattynflag[pos[pos2]]=0
     156                        pattynflag[pos[pos2]]=True
     157                        macayealpattynflag[pos[pos2]]=False
    158158
    159159                        #Recompute nodes associated to these elements
    160                         nodeonmacayeal[:]=0
    161                         nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1
    162                         nodeonpattyn[:]=0
    163                         nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1
    164                         nodeonmacayealpattyn[:]=0
    165                         nodeonmacayealpattyn[md.mesh.elements[numpy.nonzero(macayealpattynflag),:].astype(int)-1]=1
     160                        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
    166166
    167167                elif any(pattynflag) and any(stokesflag):    #coupling pattyn stokes
    168168                        #Find node at the border
    169                         nodeonpattynstokes[numpy.nonzero(numpy.logical_and(nodeonpattyn,nodeonstokes))]=1
     169                        nodeonpattynstokes[numpy.nonzero(numpy.logical_and(nodeonpattyn,nodeonstokes))]=True
    170170                        #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])
    172172                        commonelements=numpy.sum(matrixelements,axis=1)!=0
    173                         commonelements[numpy.nonzero(pattynflag)]=0    #only one layer: the elements previously in macayeal
    174                         stokesflag[numpy.nonzero(commonelements)]=0    #these elements are now macayealpattynelements
    175                         pattynstokesflag[numpy.nonzero(commonelements)]=1
    176                         nodeonstokes=numpy.zeros(md.mesh.numberofvertices)
    177                         nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1
     173                        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
    178178
    179179                        #rule out elements that don't touch the 2 boundaries
    180180                        pos=numpy.nonzero(pattynstokesflag)[0]
    181181                        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)
    184184                        pos1=numpy.nonzero(elist==1)[0]
    185                         stokesflag[pos[pos1]]=1
    186                         pattynstokesflag[pos[pos1]]=0
     185                        stokesflag[pos[pos1]]=True
     186                        pattynstokesflag[pos[pos1]]=False
    187187                        pos2=numpy.nonzero(elist==-1)[0]
    188                         pattynflag[pos[pos2]]=1
    189                         pattynstokesflag[pos[pos2]]=0
     188                        pattynflag[pos[pos2]]=True
     189                        pattynstokesflag[pos[pos2]]=False
    190190
    191191                        #Recompute nodes associated to these elements
    192                         nodeonstokes[:]=0
    193                         nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1
    194                         nodeonpattyn[:]=0
    195                         nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1
    196                         nodeonpattynstokes[:]=0
    197                         nodeonpattynstokes[md.mesh.elements[numpy.nonzero(pattynstokesflag),:].astype(int)-1]=1
     192                        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
    198198
    199199                elif any(stokesflag) and any(macayealflag):
    200200                        #Find node at the border
    201                         nodeonmacayealstokes[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonstokes))]=1
     201                        nodeonmacayealstokes[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonstokes))]=True
    202202                        #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])
    204204                        commonelements=numpy.sum(matrixelements,axis=1)!=0
    205                         commonelements[numpy.nonzero(macayealflag)]=0    #only one layer: the elements previously in macayeal
    206                         stokesflag[numpy.nonzero(commonelements)]=0    #these elements are now macayealmacayealelements
    207                         macayealstokesflag[numpy.nonzero(commonelements)]=1
    208                         nodeonstokes=numpy.zeros(md.mesh.numberofvertices)
    209                         nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1
     205                        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
    210210
    211211                        #rule out elements that don't touch the 2 boundaries
    212212                        pos=numpy.nonzero(macayealstokesflag)[0]
    213213                        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)
    216216                        pos1=numpy.nonzero(elist==1)[0]
    217                         macayealflag[pos[pos1]]=1
    218                         macayealstokesflag[pos[pos1]]=0
     217                        macayealflag[pos[pos1]]=True
     218                        macayealstokesflag[pos[pos1]]=False
    219219                        pos2=numpy.nonzero(elist==-1)[0]
    220                         stokesflag[pos[pos2]]=1
    221                         macayealstokesflag[pos[pos2]]=0
     220                        stokesflag[pos[pos2]]=True
     221                        macayealstokesflag[pos[pos2]]=False
    222222
    223223                        #Recompute nodes associated to these elements
    224                         nodeonmacayeal[:]=0
    225                         nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1
    226                         nodeonstokes[:]=0
    227                         nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1
    228                         nodeonmacayealstokes[:]=0
    229                         nodeonmacayealstokes[md.mesh.elements[numpy.nonzero(macayealstokesflag),:].astype(int)-1]=1
     224                        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
    230230
    231231                elif any(stokesflag) and any(hutterflag):
    232232                        raise TypeError("type of coupling not supported yet")
    233233
    234         #Create MacaAyealPattynApproximation where needed
    235         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)
    236236        md.flowequation.element_equation[numpy.nonzero(noneflag)]=0
    237237        md.flowequation.element_equation[numpy.nonzero(hutterflag)]=1
    … …  
    250250
    251251        #Create vertices_type
    252         md.flowequation.vertex_equation=numpy.zeros(md.mesh.numberofvertices)
     252        md.flowequation.vertex_equation=numpy.zeros(md.mesh.numberofvertices,int)
    253253        pos=numpy.nonzero(nodeonmacayeal)
    254254        md.flowequation.vertex_equation[pos]=2
    … …  
    273273
    274274        #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)
    279279
    280280        return md

  • issm/trunk/src/m/parameterization/setmask.py

    r13975 r14067  
    3737        vertexonfloatingice = numpy.zeros(md.mesh.numberofvertices,'bool')
    3838        vertexongroundedice = numpy.zeros(md.mesh.numberofvertices,'bool')
    39         vertexongroundedice[md.mesh.elements[numpy.nonzero(elementongroundedice),:].astype(int)-1]=True
     39        vertexongroundedice[md.mesh.elements[numpy.nonzero(elementongroundedice),:]-1]=True
    4040        vertexonfloatingice[numpy.nonzero(numpy.logical_not(vertexongroundedice))]=True
    4141        #}}}
    4242
    4343        #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)
    5050
    5151        return md

  • issm/trunk/src/m/plot/applyoptions.m

    r13975 r14067  
    205205                        if exist(options,'log')
    206206                                logval= getfieldvalue(options,'log');
    207                                 for i= 1:size(tick_vals,1)
     207                                for i= 1:numel(tick_vals)
    208208                                        tick_vals(i)= logval^(tick_vals(i));
    209209                                end
    210                         elseif size(tick_vals,1) == 3
     210                        elseif numel(tick_vals) == 3
    211211                                tick_vals=[tick_vals(1); mean(tick_vals(1:2)); tick_vals(2); ...
    212212                                        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
    213220                                set(c,'YTick',tick_vals);
    214221                        end

  • issm/trunk/src/m/plot/colormaps/getcolormap.m

    r13975 r14067  
    2424        map = map(128:end,:);
    2525elseif strcmpi(map,'redblue'),
    26         map = hsv;
     26        map = hsv(256);
    2727        map = rgb2hsv(map);
    2828        map(:,2)       = max(min( abs(map(:,1)-0.5)/0.5 ,1),0);

  • issm/trunk/src/m/plot/radarpower.m

    r13975 r14067  
    2020end
    2121
     22geotiff_name = getfieldvalue(options,'geotiff_name','');
    2223highres = getfieldvalue(options,'highres',0);
    2324xlim    = getfieldvalue(options,'xlim',[min(md.mesh.x) max(md.mesh.x)]);
    … …  
    6465                %md.radaroverlay.x=(xmin:(xmax-xmin)/(size(md.radaroverlay.pwr,2)-1):xmax);
    6566                %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'];
    6978                        end
    70                         geotiff_name=['/u/astrid-r1b/ModelData/MOG/mog100_r2_hp1.tif'];
    71                 else
    72                         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                         end
    75                         geotiff_name=['/u/astrid-r1b/ModelData/MOG/mog500_r2_hp1.tif'];
    7679                end
    7780
    … …  
    9194
    9295        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'];
    96107                        end
    97                         geotiff_name=['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif'];
    98                 else
    99                         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                         end
    102                         geotiff_name=['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif'];
    103108                end
    104109

  • issm/trunk/src/m/solve/WriteData.py

    r13975 r14067  
    105105                elif isinstance(data,(list,tuple)):
    106106                        data=numpy.array(data).reshape(-1,1)
    107                 if len(data.shape) == 1:
     107                if numpy.ndim(data) == 1:
    108108                        if numpy.size(data):
    109109                                data=data.reshape(numpy.size(data),1)
    … …  
    138138                elif isinstance(data,(list,tuple)):
    139139                        data=numpy.array(data).reshape(-1,1)
    140                 if len(data.shape) == 1:
     140                if numpy.ndim(data) == 1:
    141141                        if numpy.size(data):
    142142                                data=data.reshape(numpy.size(data),1)
    … …  
    171171                elif isinstance(data,(list,tuple)):
    172172                        data=numpy.array(data).reshape(-1,1)
    173                 if len(data.shape) == 1:
     173                if numpy.ndim(data) == 1:
    174174                        if numpy.size(data):
    175175                                data=data.reshape(numpy.size(data),1)
    … …  
    207207                        elif isinstance(matrix,(list,tuple)):
    208208                                matrix=numpy.array(matrix).reshape(-1,1)
    209                         if len(matrix.shape) == 1:
     209                        if numpy.ndim(matrix) == 1:
    210210                                if numpy.size(matrix):
    211211                                        matrix=matrix.reshape(numpy.size(matrix),1)
    … …  
    231231                        elif isinstance(matrix,(list,tuple)):
    232232                                matrix=numpy.array(matrix).reshape(-1,1)
    233                         if len(matrix.shape) == 1:
     233                        if numpy.ndim(matrix) == 1:
    234234                                matrix=matrix.reshape(numpy.size(matrix),1)
    235235

  • issm/trunk/src/wrappers/ElementConnectivity/ElementConnectivity.cpp

    r13236 r14067  
    1313
    1414        /*inputs: */
    15         double* elements=NULL;
    16         double* nodeconnectivity=NULL;
    17         int     nel,nods;
    18         int     width;
     15        int* elements=NULL;
     16        int* nodeconnectivity=NULL;
     17        int  nels,nods;
     18        int  width;
    1919
    2020        /*outputs: */
    21         double* elementconnectivity=NULL;
     21        int* elementconnectivity=NULL;
    2222
    2323        /*Boot module: */
    … …  
    2828       
    2929        /*Input datasets: */
    30         FetchData(&elements,&nel,NULL,ELEMENTS);
     30        FetchData(&elements,&nels,NULL,ELEMENTS);
    3131        FetchData(&nodeconnectivity,&nods,&width,NODECONNECTIVITY);
    3232
    3333        /*!Generate internal degree of freedom numbers: */
    34         ElementConnectivityx(&elementconnectivity, elements,nel, nodeconnectivity, nods, width);
     34        ElementConnectivityx(&elementconnectivity,elements,nels,nodeconnectivity,nods,width);
    3535
    3636        /*write output datasets: */
    37         WriteData(ELEMENTCONNECTIVITY,elementconnectivity,nel,3);
     37        WriteData(ELEMENTCONNECTIVITY,elementconnectivity,nels,3);
    3838
    3939        /*end module: */

  • issm/trunk/src/wrappers/Kriging/Kriging.cpp

    r13250 r14067  
    55
    66void KrigingUsage(void){/*{{{*/
     7        int num=1;
     8#ifdef _MULTITHREADING_
     9        num=_NUMTHREADS_;
     10#endif
    711        _pprintLine_("");
    812        _pprintLine_("   usage: predictions=" << __FUNCT__ << "(x,y,observations,x_interp,y_interp,'options');");
    … …  
    2125        _pprintLine_("      -'mintrimming':  minimum trimming value (default is +1.e+21)");
    2226        _pprintLine_("      -'minspacing':   minimum distance between observation (default is 0.01)");
     27        _pprintLine_("      -'numthreads':   number of threads, default is "<<num);
    2328        _pprintLine_("");
    2429}/*}}}*/

  • issm/trunk/src/wrappers/NodeConnectivity/NodeConnectivity.cpp

    r13236 r14067  
    1313
    1414        /*inputs: */
    15         double* elements=NULL;
    16         int     nel;
    17         int     nods;
     15        int* elements=NULL;
     16        int  nels;
     17        int  nods;
    1818
    1919        /*outputs: */
    20         double* connectivity=NULL;
    21         int     width;
     20        int* connectivity=NULL;
     21        int  width;
    2222
    2323        /*Boot module: */
    … …  
    2828       
    2929        /*Input datasets: */
    30         FetchData(&elements,&nel,NULL,ELEMENTS);
     30        FetchData(&elements,&nels,NULL,ELEMENTS);
    3131        FetchData(&nods,NUMNODES);
    3232
    3333        /*!Generate internal degree of freedom numbers: */
    34         NodeConnectivityx(&connectivity, &width,elements,nel, nods);
     34        NodeConnectivityx(&connectivity,&width,elements,nels,nods);
    3535
    3636        /*write output datasets: */

  • issm/trunk/src/wrappers/matlab/io/FetchMatlabData.cpp

    r13749 r14067  
    264264
    265265        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");
    282274        }
    283275
    284276        /*Assign output pointers:*/
    285277        *pvector=outvector;
    286         if (pM)*pM=outvector_rows;
     278        if(pM)*pM=M;
    287279}
    288280/*}}}*/

  • issm/trunk/src/wrappers/python/io/FetchPythonData.cpp

    r13749 r14067  
    6767        npy_intp*  dims=NULL;
    6868
    69         /*retrive dimensions: */
     69        /*intermediary:*/
     70        long* lmatrix=NULL;
     71        bool* bmatrix=NULL;
     72        int i;
     73
     74        /*retrieve dimensions: */
    7075        ndim=PyArray_NDIM((const PyArrayObject*)py_matrix);
    7176        if(ndim!=2)_error_("expecting an MxN matrix in input!");
    … …  
    7479
    7580        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!");
    82110        }
    83111        else
    … …  
    94122
    95123        /*output: */
    96         double* dmatrix=NULL;
    97124        int* matrix=NULL;
    98125        int M,N;
    99 
    100         /*intermediary:*/
    101         int i;
    102126        int ndim;
    103127        npy_intp*  dims=NULL;
    104128
    105         /*retrive dimensions: */
     129        /*intermediary:*/
     130        double* dmatrix=NULL;
     131        long* lmatrix=NULL;
     132        bool* bmatrix=NULL;
     133        int i;
     134
     135        /*retrieve dimensions: */
    106136        ndim=PyArray_NDIM((const PyArrayObject*)py_matrix);
    107137        if(ndim!=2)_error_("expecting an MxN matrix in input!");
    … …  
    110140
    111141        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!");
    118171        }
    119172        else
    … …  
    136189        npy_intp*  dims=NULL;
    137190
    138         /*retrive dimensions: */
     191        /*intermediary:*/
     192        long* lvector=NULL;
     193        bool* bvector=NULL;
     194        int i;
     195
     196        /*retrieve dimensions: */
    139197        ndim=PyArray_NDIM((const PyArrayObject*)py_vector);
    140198        if(ndim!=1)_error_("expecting an Mx1 vector in input!");
    … …  
    143201
    144202        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!");
    151232        }
    152233        else

  • issm/trunk/src/wrappers/python/io/WritePythonData.cpp

    r13749 r14067  
    4040        dims[1]=(npy_intp)N;
    4141        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){{{*/
     46void 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){{{*/
     66void 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);
    4274
    4375        PyTuple_SetItem(tuple, index, array);

  • issm/trunk/src/wrappers/python/io/pythonio.h

    r13749 r14067  
    1818
    1919void WriteData(PyObject* py_tuple,int index, double* matrix, int M,int N);
     20void WriteData(PyObject* py_tuple,int index, int* matrix, int M,int N);
     21void WriteData(PyObject* py_tuple,int index, bool* matrix, int M,int N);
    2022void WriteData(PyObject* py_tuple,int index, int integer);
    2123void WriteData(PyObject* py_tuple,int index, char* string);

  • issm/trunk/test

  • issm/trunk/test/NightlyRun/runme.m

    r13975 r14067  
    9797if strcmpi(benchmark,'nightly'),
    9898        test_ids=intersect(test_ids,[1:999]);
     99elseif strcmpi(benchmark,'validation'),
     100        test_ids=intersect(test_ids,[1001:1999]);
    99101elseif strcmpi(benchmark,'ismip'),
    100102        test_ids=intersect(test_ids,[1101:1199]);
    … …  
    105107elseif strcmpi(benchmark,'mesh'),
    106108        test_ids=intersect(test_ids,[1401:1499]);
     109elseif strcmpi(benchmark,'tranforcing'),
     110        test_ids=intersect(test_ids,[1501:1502]);
     111elseif strcmpi(benchmark,'referential'),
     112        test_ids=intersect(test_ids,[1601:1602]);
    107113elseif strcmpi(benchmark,'adolc'),
    108114        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]);
    113115end
    114116% }}}

  • issm/trunk/test/NightlyRun/runme.py

    r13975 r14067  
    9999        if   strcmpi(benchmark,'nightly'):
    100100                test_ids=test_ids.intersection(set(range(1,1000)))
     101        elif strcmpi(benchmark,'validation'):
     102                test_ids=test_ids.intersection(set(range(1001,2000)))
    101103        elif strcmpi(benchmark,'ismip'):
    102104                test_ids=test_ids.intersection(set(range(1101,1200)))
    … …  
    107109        elif strcmpi(benchmark,'mesh'):
    108110                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)))
    109115        elif strcmpi(benchmark,'adolc'):
    110116                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)))
    115117        #print 'test_ids after benchmark =',test_ids
    116118        test_ids=list(test_ids)

  • issm/trunk/test/NightlyRun/test109.py

    r13975 r14067  
    1414md=setflowequation(md,'macayeal','all')
    1515md.cluster=generic('name',oshostname(),'np',3)
    16 md.transient.isdiagnostic=0
    17 md.transient.isprognostic=0
    18 md.transient.isthermal=1
    19 md.transient.isgroundingline=0
     16md.transient.isdiagnostic=False
     17md.transient.isprognostic=False
     18md.transient.isthermal=True
     19md.transient.isgroundingline=False
    2020md=solve(md,TransientSolutionEnum())
    2121

  • issm/trunk/test/NightlyRun/test121.py

    r13975 r14067  
    1515md.cluster=generic('name',oshostname(),'np',3);
    1616md.initialization.waterfraction=numpy.zeros((md.mesh.numberofvertices,1))
    17 md.transient.isdiagnostic=0
    18 md.transient.isprognostic=0
    19 md.transient.isthermal=1
    20 md.transient.isgroundingline=0
     17md.transient.isdiagnostic=False
     18md.transient.isprognostic=False
     19md.transient.isthermal=True
     20md.transient.isgroundingline=False
    2121md.thermal.isenthalpy=1
    2222md=solve(md,TransientSolutionEnum())

  • issm/trunk/test/NightlyRun/test1501.m

    r13975 r14067  
    1313
    1414%Solve for thinning rate -> -1 * surface mass balance
    15 smb= 2.*ones(md.mesh.numberofvertices,1);   
     15smb= 2.*ones(md.mesh.numberofvertices,1);
    1616md.surfaceforcings.mass_balance= smb;
    1717md.basalforcings.melting_rate= smb;

  • issm/trunk/test/NightlyRun/test1502.m

    r13975 r14067  
    1414
    1515%Solve for thinning rate -> -1 * surface mass balance
    16 smb= 2.*ones(md.mesh.numberofvertices,1);   
     16smb= 2.*ones(md.mesh.numberofvertices,1);
    1717md.surfaceforcings.mass_balance= smb;
    1818md.basalforcings.melting_rate= smb;

  • issm/trunk/test/NightlyRun/test207.py

    r13975 r14067  
    1515md=setflowequation(md,'macayeal','all')
    1616md.cluster=generic('name',oshostname(),'np',3)
    17 md.transient.isdiagnostic=0
    18 md.transient.isprognostic=0
    19 md.transient.isthermal=1
    20 md.transient.isgroundingline=0
     17md.transient.isdiagnostic=False
     18md.transient.isprognostic=False
     19md.transient.isthermal=True
     20md.transient.isgroundingline=False
    2121md=solve(md,TransientSolutionEnum())
    2222

  • issm/trunk/test/NightlyRun/test228.py

    r13975 r14067  
    2424
    2525md.surfaceforcings.mass_balance=numpy.vstack((smb,[1.5,3.]))
    26 md.transient.isthermal=0
     26md.transient.isthermal=False
    2727
    2828md=solve(md,TransientSolutionEnum())

  • issm/trunk/test/NightlyRun/test229.py

    r13975 r14067  
    2424
    2525md.surfaceforcings.mass_balance=numpy.vstack((smb,[1.5,3.]))
    26 md.transient.isthermal=0
     26md.transient.isthermal=False
    2727
    2828md=solve(md,TransientSolutionEnum())

  • issm/trunk/test/NightlyRun/test230.py

    r13975 r14067  
    2525
    2626md.surfaceforcings.mass_balance=numpy.vstack((smb,[1.5,3.]))
    27 md.transient.isthermal=0
     27md.transient.isthermal=False
    2828
    2929md=solve(md,TransientSolutionEnum())

  • issm/trunk/test/NightlyRun/test231.py

    r13975 r14067  
    2525
    2626md.surfaceforcings.mass_balance=numpy.vstack((smb,[1.5,3.]))
    27 md.transient.isthermal=0
     27md.transient.isthermal=False
    2828
    2929md=solve(md,TransientSolutionEnum())

  • issm/trunk/test/NightlyRun/test232.py

    r13975 r14067  
    1818md.timestepping.time_step=1.
    1919md.timestepping.final_time=4.
    20 md.transient.isdiagnostic=0
    21 md.transient.isprognostic=0
    22 md.transient.isthermal=1
    23 md.transient.isgroundingline=0
     20md.transient.isdiagnostic=False
     21md.transient.isprognostic=False
     22md.transient.isthermal=True
     23md.transient.isgroundingline=False
    2424md=solve(md,TransientSolutionEnum())
    2525

  • issm/trunk/test/NightlyRun/test3009.py

    r13975 r14067  
    1414md=setflowequation(md,'macayeal','all')
    1515md.cluster=generic('name',oshostname(),'np',3)
    16 md.transient.isdiagnostic=0
    17 md.transient.isprognostic=0
    18 md.transient.isthermal=1
    19 md.transient.isgroundingline=0
     16md.transient.isdiagnostic=False
     17md.transient.isprognostic=False
     18md.transient.isthermal=True
     19md.transient.isgroundingline=False
    2020md.autodiff.isautodiff=True
    2121md=solve(md,TransientSolutionEnum())

  • issm/trunk/test/NightlyRun/test313.py

    r13975 r14067  
    1515md.cluster=generic('name',oshostname(),'np',3)
    1616md.verbose=verbose('convergence',True,'solution',True)
    17 md.transient.isdiagnostic=0
    18 md.transient.isprognostic=0
    19 md.transient.isthermal=1
    20 md.transient.isgroundingline=0
     17md.transient.isdiagnostic=False
     18md.transient.isprognostic=False
     19md.transient.isthermal=True
     20md.transient.isgroundingline=False
    2121md=solve(md,TransientSolutionEnum())
    2222

  • issm/trunk/test/NightlyRun/test326.py

    r13975 r14067  
    1616md.cluster=generic('name',oshostname(),'np',3)
    1717md.initialization.waterfraction=numpy.zeros((md.mesh.numberofvertices,1))
    18 md.transient.isdiagnostic=0
    19 md.transient.isprognostic=0
    20 md.transient.isthermal=1
    21 md.transient.isgroundingline=0
     18md.transient.isdiagnostic=False
     19md.transient.isprognostic=False
     20md.transient.isthermal=True
     21md.transient.isgroundingline=False
    2222md.thermal.isenthalpy=1
    2323md=solve(md,TransientSolutionEnum())

  • issm/trunk/test/NightlyRun/test407.py

    r13975 r14067  
    1515md=setflowequation(md,'pattyn','all')
    1616md.cluster=generic('name',oshostname(),'np',3)
    17 md.transient.isdiagnostic=0
    18 md.transient.isprognostic=0
    19 md.transient.isthermal=1
    20 md.transient.isgroundingline=0
     17md.transient.isdiagnostic=False
     18md.transient.isprognostic=False
     19md.transient.isthermal=True
     20md.transient.isgroundingline=False
    2121md=solve(md,TransientSolutionEnum())
    2222

  • issm/trunk/test/NightlyRun/test423.py

    r13975 r14067  
    2929md.cluster=generic('name',oshostname(),'np',3)
    3030
    31 md.transient.isthermal=0
    32 md.transient.isprognostic=0
    33 md.transient.isdiagnostic=0
    34 md.transient.isgroundingline=1
     31md.transient.isthermal=False
     32md.transient.isprognostic=False
     33md.transient.isdiagnostic=False
     34md.transient.isgroundingline=True
    3535
    3636#test different grounding line dynamics.

  • issm/trunk/test/NightlyRun/test424.py

    r13975 r14067  
    2020md.geometry.surface=md.geometry.bed+md.geometry.thickness
    2121md.surfaceforcings.mass_balance[:]=100.
    22 md.transient.isdiagnostic=0
    23 md.transient.isgroundingline=1
     22md.transient.isdiagnostic=False
     23md.transient.isgroundingline=True
    2424md.groundingline.migration='AgressiveMigration'
    2525

  • issm/trunk/test/NightlyRun/test425.py

    r13975 r14067  
    2020md.geometry.surface=md.geometry.bed+md.geometry.thickness
    2121md.surfaceforcings.mass_balance[:]=-150.
    22 md.transient.isdiagnostic=0
    23 md.transient.isgroundingline=1
     22md.transient.isdiagnostic=False
     23md.transient.isgroundingline=True
    2424md.groundingline.migration='SoftMigration'
    2525

  • issm/trunk/test/NightlyRun/test426.py

    r13975 r14067  
    2121md.extrude(3,1.);
    2222md=setflowequation(md,'macayeal','all');
    23 md.transient.isdiagnostic=0
    24 md.transient.isgroundingline=1
     23md.transient.isdiagnostic=False
     24md.transient.isgroundingline=True
    2525md.groundingline.migration='AgressiveMigration'
    2626md.cluster=generic('name',oshostname(),'np',3)

  • issm/trunk/test/NightlyRun/test427.py

    r13975 r14067  
    2222
    2323md.surfaceforcings.mass_balance[:]=-150
    24 md.transient.isdiagnostic=0
    25 md.transient.isgroundingline=1
     24md.transient.isdiagnostic=False
     25md.transient.isgroundingline=True
    2626md.groundingline.migration='SoftMigration'
    2727md.cluster=generic('name',oshostname(),'np',3)

  • issm/trunk/test/NightlyRun/test515.py

    r13975 r14067  
    1515md.thermal.stabilization=2
    1616md.cluster=generic('name',oshostname(),'np',3)
    17 md.transient.isdiagnostic=0
    18 md.transient.isprognostic=0
    19 md.transient.isthermal=1
    20 md.transient.isgroundingline=0
     17md.transient.isdiagnostic=False
     18md.transient.isprognostic=False
     19md.transient.isthermal=True
     20md.transient.isgroundingline=False
    2121md=solve(md,TransientSolutionEnum())
    2222
Note: See TracChangeset for help on using the changeset viewer.