Changeset 17472

Timestamp:

03/18/14 22:44:50 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: working on Tetra elements

Location:

issm/trunk-jpl/src/c

Files:

• Makefile.am (modified) (1 diff)
• analyses/StressbalanceAnalysis.cpp (modified) (26 diffs)
• classes/Elements/Element.cpp (modified) (3 diffs)
• classes/Elements/Element.h (modified) (6 diffs)
• classes/Elements/Penta.cpp (modified) (17 diffs)
• classes/Elements/Penta.h (modified) (5 diffs)
• classes/Elements/Seg.cpp (modified) (2 diffs)
• classes/Elements/Seg.h (modified) (5 diffs)
• classes/Elements/Tetra.cpp (modified) (6 diffs)
• classes/Elements/Tetra.h (modified) (7 diffs)
• classes/Elements/TetraRef.cpp (modified) (1 diff)
• classes/Elements/TetraRef.h (modified) (1 diff)
• classes/Elements/Tria.cpp (modified) (9 diffs)
• classes/Elements/Tria.h (modified) (5 diffs)
• classes/Inputs/TetraInput.cpp (added)
• classes/Inputs/TetraInput.h (added)
• classes/classes.h (modified) (1 diff)
• classes/gauss/GaussTetra.cpp (modified) (4 diffs)
• classes/gauss/GaussTetra.h (modified) (1 diff)
• modules/MeshPartitionx/MeshPartitionx.h (modified) (1 diff)
• modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (modified) (2 diffs)
• modules/ModelProcessorx/CreateNumberNodeToElementConnectivity.cpp (modified) (1 diff)
• modules/ModelProcessorx/ElementsAndVerticesPartitioning.cpp (modified) (1 diff)
• shared/Enum/EnumDefinitions.h (modified) (1 diff)
• shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• shared/Enum/StringToEnumx.cpp (modified) (5 diffs)

issm/trunk-jpl/src/c/Makefile.am

@@ -446 +446 @@
                                         ./classes/Inputs/PentaInput.h\
                                         ./classes/Inputs/PentaInput.cpp\
+                                        ./classes/Inputs/TetraInput.h\
+                                        ./classes/Inputs/TetraInput.cpp\
                                         #}}}
 #DAKOTA sources  {{{

issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

@@ -20 +20 @@
                          switch(meshtype){
                                  case Mesh3DEnum:           numdofs=2; break;
+                                 case Mesh3DtetrasEnum:     numdofs=2; break;
                                  case Mesh2DhorizontalEnum: numdofs=2; break;
                                  case Mesh2DverticalEnum:   numdofs=1; break;
@@ -29 +30 @@
                          switch(meshtype){
                                  case Mesh3DEnum:         numdofs=2; break;
+                                 case Mesh3DtetrasEnum:   numdofs=2; break;
                                  case Mesh2DverticalEnum: numdofs=1; break;
                                  default: _error_("mesh type not supported yet");
@@ -37 +39 @@
                          switch(meshtype){
                                  case Mesh3DEnum:         numdofs=3; break;
+                                 case Mesh3DtetrasEnum:   numdofs=3; break;
                                  case Mesh2DverticalEnum: numdofs=2; break;
                                  default: _error_("mesh type not supported yet");
@@ -45 +48 @@
                          switch(meshtype){
                                  case Mesh3DEnum:         numdofs=4; break;
+                                 case Mesh3DtetrasEnum:   numdofs=4; break;
                                  case Mesh2DverticalEnum: numdofs=3; break;
                                  default: _error_("mesh type not supported yet");
@@ -215 +219 @@
                 if(dakota_analysis)elements->InputDuplicate(VzEnum,QmuVzEnum);
         }
+        if(iomodel->meshtype==Mesh3DtetrasEnum){
+                iomodel->FetchDataToInput(elements,MeshVertexonbedEnum);
+                iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum);
+                iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
+                iomodel->FetchDataToInput(elements,LoadingforceZEnum);
+                iomodel->FetchDataToInput(elements,VzEnum,0.);
+                if(dakota_analysis)elements->InputDuplicate(VzEnum,QmuVzEnum);
+        }
         if(iomodel->meshtype==Mesh2DverticalEnum){
-              iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum);
+                iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum);
         }
         if(isFS){
@@ -963 +975 @@
                 case Mesh2DverticalEnum:   dim = 2; dofpernode = 1; break;
                 case Mesh3DEnum:           dim = 3; dofpernode = 2; break;
+                case Mesh3DtetrasEnum:     dim = 3; dofpernode = 2; break;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
@@ -1180 +1193 @@
                 case Mesh2DhorizontalEnum: dim = 2;break;
                 case Mesh3DEnum:           dim = 2;break;
+                case Mesh3DtetrasEnum:     dim = 2;break;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
@@ -1265 +1279 @@
                 case Mesh2DhorizontalEnum: dim = 2; bsize = 3; break;
                 case Mesh3DEnum:           dim = 2; bsize = 3; break;
+                case Mesh3DtetrasEnum:     dim = 2; bsize = 3; break;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
@@ -1377 +1392 @@
                 case Mesh2DhorizontalEnum: dim = 2;break;
                 case Mesh3DEnum:           dim = 2;break;
+                case Mesh3DtetrasEnum:     dim = 2;break;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
@@ -1442 +1458 @@
                 case Mesh2DhorizontalEnum: dim = 2;break;
                 case Mesh3DEnum:           dim = 2;break;
+                case Mesh3DtetrasEnum:     dim = 2;break;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
@@ -2140 +2157 @@
                 case Mesh2DverticalEnum: dim = 2; bsize = 2; break;
                 case Mesh3DEnum:         dim = 3; bsize = 5; break;
+                case Mesh3DtetrasEnum:   dim = 3; bsize = 5; break;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
@@ -2219 +2237 @@
                 case Mesh2DverticalEnum: dim = 2; break;
                 case Mesh3DEnum:         dim = 3; break;
+                case Mesh3DtetrasEnum:   dim = 3; break;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
@@ -2304 +2323 @@
                 case Mesh2DverticalEnum: dim = 2; break;
                 case Mesh3DEnum:         dim = 3; break;
+                case Mesh3DtetrasEnum:   dim = 3; break;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
@@ -2380 +2400 @@
                 case Mesh2DverticalEnum: dim = 2; break;
                 case Mesh3DEnum:         dim = 3; break;
+                case Mesh3DtetrasEnum:   dim = 3; break;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
@@ -2441 +2462 @@
                 case Mesh2DverticalEnum: dim = 2; break;
                 case Mesh3DEnum:         dim = 3; break;
+                case Mesh3DtetrasEnum:   dim = 3; break;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
@@ -2639 +2661 @@
                 case Mesh2DverticalEnum: dim = 2; break;
                 case Mesh3DEnum:         dim = 3; break;
+                case Mesh3DtetrasEnum:   dim = 3; break;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
@@ -2726 +2749 @@
                 case Mesh2DverticalEnum: dim = 2; break;
                 case Mesh3DEnum:         dim = 3; break;
+                case Mesh3DtetrasEnum:   dim = 3; break;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
@@ -2857 +2881 @@
                 case Mesh2DverticalEnum: dim = 2; epssize = 3; break;
                 case Mesh3DEnum:         dim = 3; epssize = 6; break;
+                case Mesh3DtetrasEnum:   dim = 3; epssize = 6; break;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
@@ -2940 +2965 @@
                 case Mesh2DverticalEnum: dim = 2;break;
                 case Mesh3DEnum:         dim = 3;break;
+                case Mesh3DtetrasEnum:   dim = 3;break;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
@@ -3034 +3060 @@
                 case Mesh2DverticalEnum: dim = 2;break;
                 case Mesh3DEnum:         dim = 3;break;
+                case Mesh3DtetrasEnum:   dim = 3;break;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
@@ -3092 +3119 @@
         element->FindParam(&meshtype,MeshTypeEnum);
         switch(meshtype){
+                case Mesh2DverticalEnum: dim = 2; break;
                 case Mesh3DEnum:         dim = 3; break;
-                case Mesh2DverticalEnum: dim = 2; break;
+                case Mesh3DtetrasEnum:   dim = 3; break;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
@@ -3176 +3204 @@
                 case Mesh2DverticalEnum: dim = 2; break;
                 case Mesh3DEnum:         dim = 3; break;
+                case Mesh3DtetrasEnum:   dim = 3; break;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
@@ -3254 +3283 @@
                 case Mesh2DverticalEnum: dim = 2; break;
                 case Mesh3DEnum:         dim = 3; break;
+                case Mesh3DtetrasEnum:   dim = 3; break;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
@@ -3330 +3360 @@
                 case Mesh2DverticalEnum: dim = 2; break;
                 case Mesh3DEnum:         dim = 3; break;
+                case Mesh3DtetrasEnum:   dim = 3; break;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
@@ -3682 +3713 @@
                 case Mesh2DverticalEnum: dim = 2; break;
                 case Mesh3DEnum:         dim = 3; break;
+                case Mesh3DtetrasEnum:   dim = 3; break;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
@@ -3749 +3781 @@
                 case Mesh2DverticalEnum: dim = 2; break;
                 case Mesh3DEnum:         dim = 3; break;
+                case Mesh3DtetrasEnum:   dim = 3; break;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }

issm/trunk-jpl/src/c/classes/Elements/Element.cpp

@@ -429 +429 @@
 
 }/*}}}*/
+void Element::GetNodesSidList(int* sidlist){/*{{{*/
+
+        _assert_(sidlist);
+        _assert_(nodes);
+        int numnodes = this->GetNumberOfNodes();
+        for(int i=0;i<numnodes;i++){
+                sidlist[i]=nodes[i]->Sid();
+        }
+}
+/*}}}*/
+void Element::GetNodesLidList(int* lidlist){/*{{{*/
+
+        _assert_(lidlist);
+        _assert_(nodes);
+        int numnodes = this->GetNumberOfNodes();
+        for(int i=0;i<numnodes;i++){
+                lidlist[i]=nodes[i]->Lid();
+        }
+}
+/*}}}*/
+void Element::GetVerticesCoordinates(IssmDouble** pxyz_list){/*{{{*/
+
+        int         numvertices = this->GetNumberOfVertices();
+        IssmDouble* xyz_list    = xNew<IssmDouble>(numvertices*3);
+        ::GetVerticesCoordinates(xyz_list,this->vertices,numvertices);
+
+        *pxyz_list = xyz_list;
+
+}/*}}}*/
 void Element::GetVerticesSidList(int* sidlist){/*{{{*/
 
@@ -439 +468 @@
         int numvertices = this->GetNumberOfVertices();
         for(int i=0;i<numvertices;i++) connectivity[i]=this->vertices[i]->Connectivity();
+}
+/*}}}*/
+void Element::InputChangeName(int original_enum,int new_enum){/*{{{*/
+        this->inputs->ChangeEnum(original_enum,new_enum);
+}
+/*}}}*/
+void Element::InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){/*{{{*/
+
+        /*Intermediaries*/
+        int        i,t,row;
+        IssmDouble time;
+
+        /*Branch on type of vector: nodal or elementary: */
+        if(vector_type==1){ //nodal vector
+
+                int         numvertices = this->GetNumberOfVertices();
+                int        *vertexids   = xNew<int>(numvertices);
+                IssmDouble *values      = xNew<IssmDouble>(numvertices);
+
+                /*Recover vertices ids needed to initialize inputs*/
+                _assert_(iomodel->elements);
+                for(i=0;i<numvertices;i++){ 
+                        vertexids[i]=reCast<int>(iomodel->elements[numvertices*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
+                }
+
+                /*Are we in transient or static? */
+                if(M==iomodel->numberofvertices){
+                        for(i=0;i<numvertices;i++) values[i]=vector[vertexids[i]-1];
+                        this->AddInput(vector_enum,values,P1Enum);
+                }
+                else if(M==iomodel->numberofvertices+1){
+                        /*create transient input: */
+                        IssmDouble* times = xNew<IssmDouble>(N);
+                        for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
+                        TransientInput* transientinput=new TransientInput(vector_enum,times,N);
+                        for(t=0;t<N;t++){
+                                for(i=0;i<numvertices;i++) values[i]=vector[N*(vertexids[i]-1)+t];
+                                switch(this->ObjectEnum()){
+                                        case TriaEnum:  transientinput->AddTimeInput(new TriaInput( vector_enum,values,P1Enum)); break;
+                                        case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,values,P1Enum)); break;
+                                        case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,values,P1Enum)); break;
+                                        default: _error_("Not implemented yet");
+                                }
+                        }
+                        this->inputs->AddInput(transientinput);
+                        xDelete<IssmDouble>(times);
+                }
+                else _error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
+
+                xDelete<IssmDouble>(values);
+                xDelete<int>(vertexids);
+        }
+        else if(vector_type==2){ //element vector
+                /*Are we in transient or static? */
+                if(M==iomodel->numberofelements){
+                        if (code==5){ //boolean
+                                this->inputs->AddInput(new BoolInput(vector_enum,reCast<bool>(vector[this->Sid()])));
+                        }
+                        else if (code==6){ //integer
+                                this->inputs->AddInput(new IntInput(vector_enum,reCast<int>(vector[this->Sid()])));
+                        }
+                        else if (code==7){ //IssmDouble
+                                this->inputs->AddInput(new DoubleInput(vector_enum,vector[this->Sid()]));
+                        }
+                        else _error_("could not recognize nature of vector from code " << code);
+                }
+                else {
+                        _error_("transient element inputs not supported yet!");
+                }
+        }
+        else{
+                _error_("Cannot add input for vector type " << vector_type << " (not supported)");
+        }
+}/*}}}*/
+void Element::InputUpdateFromConstant(int constant, int name){/*{{{*/
+
+        /*Check that name is an element input*/
+        if(!IsInput(name)) return;
+
+        /*update input*/
+        this->inputs->AddInput(new IntInput(name,constant));
+}
+/*}}}*/
+void Element::InputUpdateFromConstant(IssmDouble constant, int name){/*{{{*/
+
+        /*Check that name is an element input*/
+        if(!IsInput(name)) return;
+
+        /*update input*/
+        this->inputs->AddInput(new DoubleInput(name,constant));
+}
+/*}}}*/
+void Element::InputUpdateFromConstant(bool constant, int name){/*{{{*/
+
+        /*Check that name is an element input*/
+        if(!IsInput(name)) return;
+
+        /*update input*/
+        this->inputs->AddInput(new BoolInput(name,constant));
 }
 /*}}}*/
@@ -535 +663 @@
 }
 /*}}}*/
+ElementVector* Element::NewElementVector(int approximation_enum){/*{{{*/
+        return new ElementVector(nodes,this->GetNumberOfNodes(),this->parameters,approximation_enum);
+}
+/*}}}*/
+ElementMatrix* Element::NewElementMatrix(int approximation_enum){/*{{{*/
+        return new ElementMatrix(nodes,this->GetNumberOfNodes(),this->parameters,approximation_enum);
+}
+/*}}}*/
+ElementMatrix* Element::NewElementMatrixCoupling(int number_nodes,int approximation_enum){/*{{{*/
+        return new ElementMatrix(nodes,number_nodes,this->parameters,approximation_enum);
+}
+/*}}}*/
 void Element::ResultInterpolation(int* pinterpolation,int* pnodesperelement,int output_enum){/*{{{*/
 

issm/trunk-jpl/src/c/classes/Elements/Element.h

@@ -75 +75 @@
                 void       GetInputValue(IssmDouble* pvalue,Gauss* gauss,int enum_type);
                 IssmDouble GetMaterialParameter(int enum_in);
+                void       GetNodesSidList(int* sidlist);
+                void       GetNodesLidList(int* lidlist);
                 void       GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity);
+                void       GetVerticesCoordinates(IssmDouble** xyz_list);
                 void       GetVerticesSidList(int* sidlist);
                 void       GetVerticesConnectivityList(int* connectivitylist);
+                void       InputChangeName(int enum_type,int enum_type_old);
+                void       InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code);
+                void       InputUpdateFromConstant(IssmDouble constant, int name);
+                void       InputUpdateFromConstant(int constant, int name);
+                void       InputUpdateFromConstant(bool constant, int name);
                 bool         IsInput(int name);
                 bool       IsFloating(); 
+                ElementVector*  NewElementVector(int approximation_enum=NoneApproximationEnum);
+                ElementMatrix*  NewElementMatrix(int approximation_enum=NoneApproximationEnum);
+                ElementMatrix*  NewElementMatrixCoupling(int number_nodes,int approximation_enum=NoneApproximationEnum);
                 void       ResultInterpolation(int* pinterpolation,int*nodesperelement,int output_enum);
                 void       ResultToVector(Vector<IssmDouble>* vector,int output_enum);
@@ -161 +172 @@
                 virtual int    GetNumberOfNodesPressure(void)=0;
                 virtual int    GetNumberOfVertices(void)=0;
-                virtual void   GetNodesSidList(int* sidlist)=0;
-                virtual void   GetNodesLidList(int* sidlist)=0;
 
                 virtual int    Id()=0;
@@ -174 +183 @@
                 virtual void   GetInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0;
                 virtual Node*  GetNode(int node_number)=0;
-                virtual void   GetVerticesCoordinates(IssmDouble** xyz_list)=0;
                 virtual void   GetVerticesCoordinatesBase(IssmDouble** xyz_list)=0;
                 virtual void   GetVerticesCoordinatesTop(IssmDouble** xyz_list)=0;
@@ -184 +192 @@
                 virtual IssmDouble SurfaceArea(void)=0;
                 virtual void   InputDepthAverageAtBase(int enum_type,int average_enum_type)=0;
-                virtual void   InputChangeName(int enum_type,int enum_type_old)=0;
                 virtual void   ComputeBasalStress(Vector<IssmDouble>* sigma_b)=0;
                 virtual void   ComputeStrainRate(Vector<IssmDouble>* eps)=0;
@@ -191 +198 @@
                 virtual void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finite_element)=0;
                 virtual void   InputDuplicate(int original_enum,int new_enum)=0;
-                virtual void   InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code)=0;
                 virtual void   InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int inputenum)=0;
                 virtual void   InputUpdateFromSolutionOneDof(IssmDouble* solution,int inputenum)=0;
@@ -206 +212 @@
                 virtual Gauss* NewGaussLine(int vertex1,int vertex2,int order)=0;
                 virtual Gauss* NewGaussTop(int order)=0;
-                virtual ElementVector*  NewElementVector(int approximation_enum=NoneApproximationEnum)=0;
-                virtual ElementMatrix*  NewElementMatrix(int approximation_enum=NoneApproximationEnum)=0;
-                virtual ElementMatrix*  NewElementMatrixCoupling(int number_nodes,int approximation_enum=NoneApproximationEnum)=0;
+
                 virtual void   InputScale(int enum_type,IssmDouble scale_factor)=0;
                 virtual void   GetVectorFromInputs(Vector<IssmDouble>* vector, int name_enum)=0;

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

@@ -28 +28 @@
 Penta::Penta(int penta_id, int penta_sid, int index, IoModel* iomodel,int nummodels)
         :PentaRef(nummodels)
-        ,ElementHook(nummodels,index+1,6,iomodel){
+        ,ElementHook(nummodels,index+1,NUMVERTICES,iomodel){
 
         int penta_elements_ids[2];
@@ -843 +843 @@
 }
 /*}}}*/
-/*FUNCTION Penta::GetNodesSidList{{{*/
-void Penta::GetNodesSidList(int* sidlist){
-
-        _assert_(sidlist);
-        _assert_(nodes);
-        int numnodes = this->NumberofNodes();
-
-        for(int i=0;i<numnodes;i++){
-                sidlist[i]=nodes[i]->Sid();
-        }
-}
-/*}}}*/
-/*FUNCTION Penta::GetNodesLidList{{{*/
-void Penta::GetNodesLidList(int* lidlist){
-
-        _assert_(lidlist);
-        _assert_(nodes);
-        int numnodes = this->NumberofNodes();
-
-        for(int i=0;i<numnodes;i++){
-                lidlist[i]=nodes[i]->Lid();
-        }
-}
-/*}}}*/
 /*FUNCTION Penta::GetNumberOfNodes;{{{*/
 int Penta::GetNumberOfNodes(void){
@@ -907 +883 @@
 }
 /*}}}*/
-/*FUNCTION Penta::GetVerticesCoordinates(IssmDouble** pxyz_list){{{*/
-void Penta::GetVerticesCoordinates(IssmDouble** pxyz_list){
-
-        IssmDouble* xyz_list = xNew<IssmDouble>(NUMVERTICES*3);
-        ::GetVerticesCoordinates(xyz_list,this->vertices,NUMVERTICES);
-
-        /*Assign output pointer*/
-        *pxyz_list = xyz_list;
-
-}/*}}}*/
 /*FUNCTION Penta::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){{{*/
 void Penta::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){
@@ -953 +919 @@
 
         cross(normal,AB,AC);
-        norm=sqrt(pow(normal[0],2.0)+pow(normal[1],2.0)+pow(normal[2],2.0));
+        norm=sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
 
         for(int i=0;i<3;i++) normal[i]=normal[i]/norm;
@@ -1206 +1172 @@
 int    Penta::Id(void){
         return id; 
-}
-/*}}}*/
-/*FUNCTION Penta::InputChangeName{{{*/
-void  Penta::InputChangeName(int original_enum,int new_enum){
-
-        /*Call inputs method*/
-        this->inputs->ChangeEnum(original_enum,new_enum);
-}
-/*}}}*/
-/*FUNCTION Penta::InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){{{*/
-void Penta::InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){
-
-        /*Intermediaries*/
-        int             i,t;
-        int             penta_vertex_ids[6];
-        int             row;
-        IssmDouble      nodeinputs[6];
-        IssmDouble      time;
-        TransientInput *transientinput      = NULL;
-
-        /*Branch on type of vector: nodal or elementary: */
-        if(vector_type==1){ //nodal vector
-
-                /*Recover vertices ids needed to initialize inputs*/
-                for(i=0;i<6;i++){ 
-                        _assert_(iomodel->elements);
-                        penta_vertex_ids[i]=iomodel->elements[6*this->sid+i]; //ids for vertices are in the elements array from Matlab
-                }
-
-                /*Are we in transient or static? */
-                if(M==iomodel->numberofvertices){
-
-                        /*create input values: */
-                        for(i=0;i<6;i++)nodeinputs[i]=(IssmDouble)vector[penta_vertex_ids[i]-1];
-
-                        /*create static input: */
-                        this->inputs->AddInput(new PentaInput(vector_enum,nodeinputs,P1Enum));
-                }
-                else if(M==iomodel->numberofvertices+1){
-                        /*create transient input: */
-                        IssmDouble* times = xNew<IssmDouble>(N);
-                        for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
-                        transientinput=new TransientInput(vector_enum,times,N);
-                        for(t=0;t<N;t++){
-                                for(i=0;i<NUMVERTICES;i++) nodeinputs[i]=vector[N*(penta_vertex_ids[i]-1)+t];
-                                transientinput->AddTimeInput(new PentaInput(vector_enum,nodeinputs,P1Enum));
-                        }
-                        this->inputs->AddInput(transientinput);
-                        xDelete<IssmDouble>(times);
-                }
-                else _error_("nodal vector is either numberofvertices (" << iomodel->numberofvertices << "), or numberofvertices+1 long. Field provided is " << M << " long. Enum " << EnumToStringx(vector_enum));
-        }
-        else if(vector_type==2){ //element vector
-                /*Are we in transient or static? */
-                if(M==iomodel->numberofelements){
-
-                        /*static mode: create an input out of the element value: */
-
-                        if (code==5){ //boolean
-                                this->inputs->AddInput(new BoolInput(vector_enum,reCast<bool,IssmDouble>(vector[this->sid])));
-                        }
-                        else if (code==6){ //integer
-                                this->inputs->AddInput(new IntInput(vector_enum,reCast<int,IssmDouble>(vector[this->sid])));
-                        }
-                        else if (code==7){ //IssmDouble
-                                this->inputs->AddInput(new DoubleInput(vector_enum,vector[this->sid]));
-                        }
-                        else _error_("could not recognize nature of vector from code " << code);
-                }
-                else {
-                        _error_("transient elementary inputs not supported yet!");
-                }
-        }
-        else{
-                _error_("Cannot add input for vector type " << vector_type << " (not supported)");
-        }
-
 }
 /*}}}*/
@@ -1452 +1341 @@
 }
 /*}}}*/
-/*FUNCTION Penta::InputUpdateFromConstant(bool value, int name);{{{*/
-void  Penta::InputUpdateFromConstant(bool constant, int name){
-
-        /*Check that name is an element input*/
-        if (!IsInput(name)) return;
-
-        /*update input*/
-        this->inputs->AddInput(new BoolInput(name,constant));
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromConstant(IssmDouble value, int name);{{{*/
-void  Penta::InputUpdateFromConstant(IssmDouble constant, int name){
-        /*Check that name is an element input*/
-        if (!IsInput(name)) return;
-
-        /*update input*/
-        this->inputs->AddInput(new DoubleInput(name,constant));
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromConstant(int value, int name);{{{*/
-void  Penta::InputUpdateFromConstant(int constant, int name){
-        /*Check that name is an element input*/
-        if (!IsInput(name)) return;
-
-        /*update input*/
-        this->inputs->AddInput(new IntInput(name,constant));
-}
-/*}}}*/
 /*FUNCTION Penta::InputUpdateFromIoModel {{{*/
 void Penta::InputUpdateFromIoModel(int index,IoModel* iomodel){ 
 
         /*Intermediaries*/
-        IssmInt i,j;
-        int     penta_vertex_ids[6];
-        IssmDouble  nodeinputs[6];
-        IssmDouble  cmmininputs[6];
-        IssmDouble  cmmaxinputs[6];
+        int         i,j;
+        int         penta_vertex_ids[NUMVERTICES];
+        IssmDouble  nodeinputs[NUMVERTICES];
+        IssmDouble  cmmininputs[NUMVERTICES];
+        IssmDouble  cmmaxinputs[NUMVERTICES];
 
         IssmDouble  yts;
@@ -1500 +1361 @@
         if(control_analysis) iomodel->Constant(&num_responses,InversionNumCostFunctionsEnum);
 
-        /*Checks if debuging*/
-        /*{{{*/
+        /*Recover vertices ids needed to initialize inputs*/
         _assert_(iomodel->elements);
-        /*}}}*/
-
-        /*Recover vertices ids needed to initialize inputs*/
-        for(i=0;i<6;i++){ 
-                penta_vertex_ids[i]=iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab
+        for(i=0;i<NUMVERTICES;i++){ 
+                penta_vertex_ids[i]=iomodel->elements[NUMVERTICES*index+i]; //ids for vertices are in the elements array from Matlab
         }
 
@@ -1516 +1373 @@
                                 case BalancethicknessThickeningRateEnum:
                                         if (iomodel->Data(BalancethicknessThickeningRateEnum)){
-                                                for(j=0;j<6;j++)nodeinputs[j]=iomodel->Data(BalancethicknessThickeningRateEnum)[penta_vertex_ids[j]-1]/yts;
-                                                for(j=0;j<6;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
-                                                for(j=0;j<6;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
+                                                for(j=0;j<NUMVERTICES;j++)nodeinputs[j]=iomodel->Data(BalancethicknessThickeningRateEnum)[penta_vertex_ids[j]-1]/yts;
+                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
+                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
                                                 this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
                                         }
@@ -1524 +1381 @@
                                 case VxEnum:
                                         if (iomodel->Data(VxEnum)){
-                                                for(j=0;j<6;j++)nodeinputs[j]=iomodel->Data(VxEnum)[penta_vertex_ids[j]-1]/yts;
-                                                for(j=0;j<6;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
-                                                for(j=0;j<6;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
+                                                for(j=0;j<NUMVERTICES;j++)nodeinputs[j]=iomodel->Data(VxEnum)[penta_vertex_ids[j]-1]/yts;
+                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
+                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
                                                 this->inputs->AddInput(new ControlInput(VxEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
                                         }
@@ -1532 +1389 @@
                                 case VyEnum:
                                         if (iomodel->Data(VyEnum)){
-                                                for(j=0;j<6;j++)nodeinputs[j]=iomodel->Data(VyEnum)[penta_vertex_ids[j]-1]/yts;
-                                                for(j=0;j<6;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
-                                                for(j=0;j<6;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
+                                                for(j=0;j<NUMVERTICES;j++)nodeinputs[j]=iomodel->Data(VyEnum)[penta_vertex_ids[j]-1]/yts;
+                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
+                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
                                                 this->inputs->AddInput(new ControlInput(VyEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
                                         }
@@ -1540 +1397 @@
                                 case FrictionCoefficientEnum:
                                         if (iomodel->Data(FrictionCoefficientEnum)){
-                                                for(j=0;j<6;j++)nodeinputs[j]=iomodel->Data(FrictionCoefficientEnum)[penta_vertex_ids[j]-1];
-                                                for(j=0;j<6;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i];
-                                                for(j=0;j<6;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i];
+                                                for(j=0;j<NUMVERTICES;j++)nodeinputs[j]=iomodel->Data(FrictionCoefficientEnum)[penta_vertex_ids[j]-1];
+                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i];
+                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i];
                                                 this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
                                         }
@@ -1548 +1405 @@
                                 case MaterialsRheologyBbarEnum:
                                         if(iomodel->Data(MaterialsRheologyBEnum)){
-                                                for(j=0;j<6;j++) nodeinputs[j]=iomodel->Data(MaterialsRheologyBEnum)[penta_vertex_ids[j]-1];
-                                                for(j=0;j<6;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i];
-                                                for(j=0;j<6;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i];
+                                                for(j=0;j<NUMVERTICES;j++) nodeinputs[j]=iomodel->Data(MaterialsRheologyBEnum)[penta_vertex_ids[j]-1];
+                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i];
+                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i];
                                                 this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
                                         }
@@ -1556 +1413 @@
                                 case DamageDbarEnum:
                                         if(iomodel->Data(DamageDEnum)){
-                                                for(j=0;j<6;j++) nodeinputs[j]=iomodel->Data(DamageDEnum)[penta_vertex_ids[j]-1];
-                                                for(j=0;j<6;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i];
-                                                for(j=0;j<6;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i];
+                                                for(j=0;j<NUMVERTICES;j++) nodeinputs[j]=iomodel->Data(DamageDEnum)[penta_vertex_ids[j]-1];
+                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i];
+                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i];
                                                 this->inputs->AddInput(new ControlInput(DamageDEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
                                         }
@@ -1579 +1436 @@
                 DatasetInput* datasetinput=new DatasetInput(InversionCostFunctionsCoefficientsEnum);
                 for(i=0;i<num_responses;i++){
-                        for(j=0;j<6;j++)nodeinputs[j]=iomodel->Data(InversionCostFunctionsCoefficientsEnum)[(penta_vertex_ids[j]-1)*num_responses+i];
+                        for(j=0;j<NUMVERTICES;j++)nodeinputs[j]=iomodel->Data(InversionCostFunctionsCoefficientsEnum)[(penta_vertex_ids[j]-1)*num_responses+i];
                         datasetinput->AddInput(new PentaInput(InversionCostFunctionsCoefficientsEnum,nodeinputs,P1Enum),reCast<int>(iomodel->Data(InversionCostFunctionsEnum)[i]));
                 }
@@ -1887 +1744 @@
 Gauss* Penta::NewGaussTop(int order){
         return new GaussPenta(3,4,5,order);
-}
-/*}}}*/
-/*FUNCTION Penta::NewElementVector{{{*/
-ElementVector* Penta::NewElementVector(int approximation_enum){
-        return new ElementVector(nodes,this->NumberofNodes(),this->parameters,approximation_enum);
-}
-/*}}}*/
-/*FUNCTION Penta::NewElementMatrix{{{*/
-ElementMatrix* Penta::NewElementMatrix(int approximation_enum){
-        return new ElementMatrix(nodes,this->NumberofNodes(),this->parameters,approximation_enum);
-}
-/*}}}*/
-/*FUNCTION Penta::NewElementMatrixCoupling{{{*/
-ElementMatrix* Penta::NewElementMatrixCoupling(int number_nodes,int approximation_enum){
-        return new ElementMatrix(nodes,number_nodes,this->parameters,approximation_enum);
 }
 /*}}}*/
@@ -2212 +2054 @@
 void  Penta::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, Nodes* nodesin, Materials* materialsin, Parameters* parametersin){
 
+        /*go into parameters and get the analysis_counter: */
         int analysis_counter;
-
-        /*go into parameters and get the analysis_counter: */
         parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);
 
@@ -2220 +2061 @@
         this->element_type=this->element_type_list[analysis_counter];
 
-        /*Pick up nodes */
-        if (this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
+        /*Pick up nodes*/
+        if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
         else this->nodes=NULL;
+
 }
 /*}}}*/

issm/trunk-jpl/src/c/classes/Elements/Penta.h

@@ -50 +50 @@
                 /*}}}*/
                 /*Update virtual functions definitions: {{{*/
-                void  InputUpdateFromConstant(bool constant, int name);
-                void  InputUpdateFromConstant(IssmDouble constant, int name);
-                void  InputUpdateFromConstant(int constant, int name);
                 void  InputUpdateFromVector(IssmDouble* vector, int name, int type);
                 #ifdef _HAVE_DAKOTA_
@@ -85 +82 @@
                 IssmDouble GetGroundedPortion(IssmDouble* xyz_list);
                 int    GetNodeIndex(Node* node);
-                void   GetNodesSidList(int* sidlist);
-                void   GetNodesLidList(int* lidlist);
                 int    GetNumberOfNodes(void);
                 int    GetNumberOfNodesPressure(void);
@@ -96 +91 @@
                 IssmDouble GetZcoord(Gauss* gauss);
                 void   GetVectorFromInputs(Vector<IssmDouble>* vector,int name_enum);
-                void   GetVerticesCoordinates(IssmDouble** pxyz_list);
                 void   GetVerticesCoordinatesBase(IssmDouble** pxyz_list);
                 void   GetVerticesCoordinatesTop(IssmDouble** pxyz_list);
 
                 int    Sid();
-                void   InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code);
                 void   InputDepthAverageAtBase(int enum_type,int average_enum_type);
                 void   InputDuplicate(int original_enum,int new_enum);
@@ -204 +197 @@
                 void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
                 Node*          GetNode(int node_number);
-                void           InputChangeName(int input_enum, int enum_type_old);
                 void             InputExtrude(int enum_type);
                 void           InputUpdateFromSolutionOneDof(IssmDouble* solutiong,int enum_type);
@@ -225 +217 @@
                 Gauss*         NewGaussLine(int vertex1,int vertex2,int order);
                 Gauss*         NewGaussTop(int order);
-                ElementVector* NewElementVector(int approximation_enum);
-                ElementMatrix* NewElementMatrix(int approximation_enum);
-                ElementMatrix* NewElementMatrixCoupling(int number_nodes,int approximation_enum);
                 void           NodalFunctions(IssmDouble* basis,Gauss* gauss);
                 void           NodalFunctionsP1(IssmDouble* basis,Gauss* gauss);

issm/trunk-jpl/src/c/classes/Elements/Seg.cpp

@@ -21 +21 @@
 /*FUNCTION Seg::Seg(int id, int sid,int index, IoModel* iomodel,int nummodels){{{*/
 Seg::Seg(int seg_id, int seg_sid, int index, IoModel* iomodel,int nummodels)
-                :SegRef(nummodels),ElementHook(nummodels,index+1,2,iomodel){
+                :SegRef(nummodels),ElementHook(nummodels,index+1,NUMVERTICES,iomodel){
 
                         /*id: */
@@ -263 +263 @@
 }
 /*}}}*/
-/*FUNCTION Seg::NewElementVector{{{*/
-ElementVector* Seg::NewElementVector(int approximation_enum){
-        return new ElementVector(nodes,this->NumberofNodes(),this->parameters,approximation_enum);
-}
-/*}}}*/
-/*FUNCTION Seg::NewElementMatrix{{{*/
-ElementMatrix* Seg::NewElementMatrix(int approximation_enum){
-        return new ElementMatrix(nodes,this->NumberofNodes(),this->parameters,approximation_enum);
-}
-/*}}}*/
 /*FUNCTION Seg::NodalFunctions{{{*/
 void Seg::NodalFunctions(IssmDouble* basis, Gauss* gauss){

issm/trunk-jpl/src/c/classes/Elements/Seg.h

@@ -53 +53 @@
                 void  InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type){_error_("not implemented yet");};
 #endif
-                void  InputUpdateFromConstant(IssmDouble constant, int name){_error_("not implemented yet");};
-                void  InputUpdateFromConstant(int constant, int name){_error_("not implemented yet");};
-                void  InputUpdateFromConstant(bool constant, int name){_error_("not implemented yet");};
                 void  InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
                 /*}}}*/
@@ -80 +77 @@
                 Element*    GetBasalElement(void){_error_("not implemented yet");};
                 int         GetNodeIndex(Node* node){_error_("not implemented yet");};
-                void        GetNodesSidList(int* sidlist){_error_("not implemented yet");};
-                void        GetNodesLidList(int* lidlist){_error_("not implemented yet");};
                 int         GetNumberOfNodes(void);
                 int         GetNumberOfNodesVelocity(void){_error_("not implemented yet");};
@@ -90 +85 @@
                 void        GetVerticesCoordinatesTop(IssmDouble** pxyz_list){_error_("not implemented yet");};
                 int         Sid(){_error_("not implemented yet");};
-                void        InputChangeName(int input_enum, int enum_type_old){_error_("not implemented yet");};
                 bool        IsOnBed(){_error_("not implemented yet");};
                 bool        IsOnSurface(){_error_("not implemented yet");};
@@ -138 +132 @@
                 Gauss*      NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");};
                 Gauss*      NewGaussTop(int order){_error_("not implemented yet");};
-                ElementVector* NewElementVector(int approximation_enum);
-                ElementMatrix* NewElementMatrix(int approximation_enum);
-                ElementMatrix* NewElementMatrixCoupling(int number_nodes,int approximation_enum){_error_("not implemented yet");};
                 int         VertexConnectivity(int vertexindex){_error_("not implemented yet");};
                 void        VerticalSegmentIndices(int** pindices,int* pnumseg){_error_("not implemented yet");};
@@ -152 +143 @@
                 void        GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type){_error_("not implemented yet");};
                 void        GetVectorFromInputs(Vector<IssmDouble>* vector, int name_enum){_error_("not implemented yet");};
-                void        InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){_error_("not implemented yet");};
                 void        InputDepthAverageAtBase(int enum_type,int average_enum_type){_error_("not implemented yet");};
                 void        InputDuplicate(int original_enum,int new_enum){_error_("not implemented yet");};

issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp

@@ -22 +22 @@
 /*FUNCTION Tetra::Tetra(int id, int sid,int index, IoModel* iomodel,int nummodels){{{*/
 Tetra::Tetra(int seg_id, int seg_sid, int index, IoModel* iomodel,int nummodels)
-                :TetraRef(nummodels),ElementHook(nummodels,index+1,2,iomodel){
+                :TetraRef(nummodels),ElementHook(nummodels,index+1,NUMVERTICES,iomodel){
 
                         /*id: */
@@ -127 +127 @@
 /*}}}*/
 
+/*FUNCTION Tetra::AddInput{{{*/
+void  Tetra::AddInput(int input_enum,IssmDouble* values, int interpolation_enum){
+
+        /*Call inputs method*/
+        _assert_(this->inputs);
+        this->inputs->AddInput(new TetraInput(input_enum,values,interpolation_enum));
+}
+/*}}}*/
+/*FUNCTION Tetra::Configure {{{*/
+void  Tetra::Configure(Elements* elementsin, Loads* loadsin, Nodes* nodesin,Vertices* verticesin, Materials* materialsin, Parameters* parametersin){
+
+        int analysis_counter;
+
+        /*go into parameters and get the analysis_counter: */
+        parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);
+
+        /*Get Element type*/
+        this->element_type=this->element_type_list[analysis_counter];
+
+        /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective 
+         * datasets, using internal ids and offsets hidden in hooks: */
+        if (this->hnodes[analysis_counter]) this->hnodes[analysis_counter]->configure(nodesin);
+        this->hvertices->configure(verticesin);
+        this->hmaterial->configure(materialsin);
+        this->hmatpar->configure(materialsin);
+
+        /*Now, go pick up the objects inside the hooks: */
+        if (this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
+        else this->nodes=NULL;
+        this->vertices          = (Vertex**)this->hvertices->deliverp();
+        this->material          = (Material*)this->hmaterial->delivers();
+        this->matpar            = (Matpar*)this->hmatpar->delivers();
+
+        /*point parameters to real dataset: */
+        this->parameters=parametersin;
+
+        /*get inputs configured too: */
+        this->inputs->Configure(parameters);
+}
+/*}}}*/
+/*FUNCTION Tetra::FaceOnBedIndices{{{*/
+void Tetra::FaceOnBedIndices(int* pindex1,int* pindex2,int* pindex3){
+
+        IssmDouble values[NUMVERTICES];
+        int        indices[4][3] = {{0,1,2},{0,1,3},{1,2,3},{2,0,3}};
+
+        /*Retrieve all inputs and parameters*/
+        GetInputListOnVertices(&values[0],MeshVertexonbedEnum);
+
+        for(int i=0;i<4;i++){
+                if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1. && values[indices[i][2]] == 1.){
+                        *pindex1 = indices[i][0];
+                        *pindex2 = indices[i][1];
+                        *pindex3 = indices[i][2];
+                        return;
+                }
+        }
+
+        _error_("Could not find 3 vertices on bed");
+}
+/*}}}*/
+/*FUNCTION Tetra::FaceOnSurfaceIndices{{{*/
+void Tetra::FaceOnSurfaceIndices(int* pindex1,int* pindex2,int* pindex3){
+
+        IssmDouble values[NUMVERTICES];
+        int        indices[4][3] = {{0,1,2},{0,1,3},{1,2,3},{2,0,3}};
+
+        /*Retrieve all inputs and parameters*/
+        GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
+
+        for(int i=0;i<4;i++){
+                if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1. && values[indices[i][2]] == 1.){
+                        *pindex1 = indices[i][0];
+                        *pindex2 = indices[i][1];
+                        *pindex3 = indices[i][2];
+                        return;
+                }
+        }
+
+        _error_("Could not find 3 vertices on bed");
+}
+/*}}}*/
+/*FUNCTION Tetra::FaceOnFranceIndices{{{*/
+void Tetra::FaceOnFrontIndices(int* pindex1,int* pindex2,int* pindex3){
+
+        IssmDouble values[NUMVERTICES];
+        int        indices[4][3] = {{0,1,2},{0,1,3},{1,2,3},{2,0,3}};
+
+        /*Retrieve all inputs and parameters*/
+        GetInputListOnVertices(&values[0],MaskIceLevelsetEnum);
+
+        for(int i=0;i<4;i++){
+                if(values[indices[i][0]] == 0. && values[indices[i][1]] == 0. && values[indices[i][2]] == 0.){
+                        *pindex1 = indices[i][0];
+                        *pindex2 = indices[i][1];
+                        *pindex3 = indices[i][2];
+                        return;
+                }
+        }
+
+        _error_("Could not find 3 vertices on bed");
+}
+/*}}}*/
 /*FUNCTION Tetra::GetNumberOfNodes;{{{*/
 int Tetra::GetNumberOfNodes(void){
@@ -137 +240 @@
 }
 /*}}}*/
-/*FUNCTION Tetra::GetVerticesCoordinates(IssmDouble** pxyz_list)   THIS ONE{{{*/
-void Tetra::GetVerticesCoordinates(IssmDouble** pxyz_list){
-
-        IssmDouble* xyz_list = xNew<IssmDouble>(NUMVERTICES*3);
-        ::GetVerticesCoordinates(xyz_list,this->vertices,NUMVERTICES);
+/*FUNCTION Tetra::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){{{*/
+void Tetra::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){
+
+        int        indices[3];
+        IssmDouble xyz_list[NUMVERTICES][3];
+
+        /*Element XYZ list*/
+        ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
+
+        /*Allocate Output*/
+        IssmDouble* xyz_list_edge = xNew<IssmDouble>(3*3);
+        this->FaceOnBedIndices(&indices[0],&indices[1],&indices[2]);
+        for(int i=0;i<3;i++) for(int j=0;j<3;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j];
 
         /*Assign output pointer*/
-        *pxyz_list = xyz_list;
+        *pxyz_list = xyz_list_edge;
 
 }/*}}}*/
-/*FUNCTION Tetra::IsIceInElement {{{*/
+/*FUNCTION Tetra::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){{{*/
+void Tetra::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){
+
+        int        indices[3];
+        IssmDouble xyz_list[NUMVERTICES][3];
+
+        /*Element XYZ list*/
+        ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
+
+        /*Allocate Output*/
+        IssmDouble* xyz_list_edge = xNew<IssmDouble>(3*3);
+        this->FaceOnSurfaceIndices(&indices[0],&indices[1],&indices[2]);
+        for(int i=0;i<3;i++) for(int j=0;j<3;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j];
+
+        /*Assign output pointer*/
+        *pxyz_list = xyz_list_edge;
+
+}/*}}}*/
+/*FUNCTION Tetra::GetZcoord {{{*/
+IssmDouble Tetra::GetZcoord(Gauss* gauss){
+
+        int    i;
+        IssmDouble z;
+        IssmDouble xyz_list[NUMVERTICES][3];
+        IssmDouble z_list[NUMVERTICES];
+
+        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+        for(i=0;i<NUMVERTICES;i++) z_list[i]=xyz_list[i][2];
+        TetraRef::GetInputValue(&z,&z_list[0],(GaussTetra*)gauss,P1Enum);
+
+        return z;
+}
+/*}}}*/
+/*FUNCTION Tetra::HasFaceOnBed{{{*/
+bool Tetra::HasFaceOnBed(){
+
+        IssmDouble values[NUMVERTICES];
+        IssmDouble sum;
+
+        /*Retrieve all inputs and parameters*/
+        GetInputListOnVertices(&values[0],MeshVertexonbedEnum);
+        sum = values[0]+values[1]+values[2]+values[3];
+
+        _assert_(sum==0. || sum==1. || sum==2. || sum==3.);
+
+        if(sum==3){
+                return true;
+        }
+        else{
+                return false;
+        }
+}
+/*}}}*/
+/*FUNCTION Tetra::HasFaceOnSurface{{{*/
+bool Tetra::HasFaceOnSurface(){
+
+        IssmDouble values[NUMVERTICES];
+        IssmDouble sum;
+
+        /*Retrieve all inputs and parameters*/
+        GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
+        sum = values[0]+values[1]+values[2]+values[3];
+
+        _assert_(sum==0. || sum==1. || sum==2. || sum==3.);
+
+        if(sum==3){
+                return true;
+        }
+        else{
+                return false;
+        }
+}
+/*}}}*/
+/*FUNCTION Tetra::InputUpdateFromIoModel {{{*/
+void Tetra::InputUpdateFromIoModel(int index,IoModel* iomodel){ 
+
+        /*Intermediaries*/
+        int         i,j;
+        int         tetra_vertex_ids[NUMVERTICES];
+        IssmDouble  nodeinputs[NUMVERTICES];
+        IssmDouble  cmmininputs[NUMVERTICES];
+        IssmDouble  cmmaxinputs[NUMVERTICES];
+
+        IssmDouble  yts;
+        bool    control_analysis;
+        int     num_control_type,num_responses;
+
+        /*Fetch parameters: */
+        iomodel->Constant(&yts,ConstantsYtsEnum);
+        iomodel->Constant(&control_analysis,InversionIscontrolEnum);
+        if(control_analysis) iomodel->Constant(&num_control_type,InversionNumControlParametersEnum);
+        if(control_analysis) iomodel->Constant(&num_responses,InversionNumCostFunctionsEnum);
+
+        /*Recover vertices ids needed to initialize inputs*/
+        _assert_(iomodel->elements);
+        for(i=0;i<NUMVERTICES;i++){ 
+                tetra_vertex_ids[i]=iomodel->elements[NUMVERTICES*index+i]; //ids for vertices are in the elements array from Matlab
+        }
+
+        /*Control Inputs*/
+        if (control_analysis && iomodel->Data(InversionControlParametersEnum)){
+                for(i=0;i<num_control_type;i++){
+                        switch(reCast<int,IssmDouble>(iomodel->Data(InversionControlParametersEnum)[i])){
+                                case BalancethicknessThickeningRateEnum:
+                                        if (iomodel->Data(BalancethicknessThickeningRateEnum)){
+                                                for(j=0;j<NUMVERTICES;j++)nodeinputs[j]=iomodel->Data(BalancethicknessThickeningRateEnum)[tetra_vertex_ids[j]-1]/yts;
+                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tetra_vertex_ids[j]-1)*num_control_type+i]/yts;
+                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tetra_vertex_ids[j]-1)*num_control_type+i]/yts;
+                                                this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
+                                        }
+                                        break;
+                                case VxEnum:
+                                        if (iomodel->Data(VxEnum)){
+                                                for(j=0;j<NUMVERTICES;j++)nodeinputs[j]=iomodel->Data(VxEnum)[tetra_vertex_ids[j]-1]/yts;
+                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tetra_vertex_ids[j]-1)*num_control_type+i]/yts;
+                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tetra_vertex_ids[j]-1)*num_control_type+i]/yts;
+                                                this->inputs->AddInput(new ControlInput(VxEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
+                                        }
+                                        break;
+                                case VyEnum:
+                                        if (iomodel->Data(VyEnum)){
+                                                for(j=0;j<NUMVERTICES;j++)nodeinputs[j]=iomodel->Data(VyEnum)[tetra_vertex_ids[j]-1]/yts;
+                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tetra_vertex_ids[j]-1)*num_control_type+i]/yts;
+                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tetra_vertex_ids[j]-1)*num_control_type+i]/yts;
+                                                this->inputs->AddInput(new ControlInput(VyEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
+                                        }
+                                        break;
+                                case FrictionCoefficientEnum:
+                                        if (iomodel->Data(FrictionCoefficientEnum)){
+                                                for(j=0;j<NUMVERTICES;j++)nodeinputs[j]=iomodel->Data(FrictionCoefficientEnum)[tetra_vertex_ids[j]-1];
+                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tetra_vertex_ids[j]-1)*num_control_type+i];
+                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tetra_vertex_ids[j]-1)*num_control_type+i];
+                                                this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
+                                        }
+                                        break;
+                                case MaterialsRheologyBbarEnum:
+                                        if(iomodel->Data(MaterialsRheologyBEnum)){
+                                                for(j=0;j<NUMVERTICES;j++) nodeinputs[j]=iomodel->Data(MaterialsRheologyBEnum)[tetra_vertex_ids[j]-1];
+                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tetra_vertex_ids[j]-1)*num_control_type+i];
+                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tetra_vertex_ids[j]-1)*num_control_type+i];
+                                                this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
+                                        }
+                                        break;
+                                case DamageDbarEnum:
+                                        if(iomodel->Data(DamageDEnum)){
+                                                for(j=0;j<NUMVERTICES;j++) nodeinputs[j]=iomodel->Data(DamageDEnum)[tetra_vertex_ids[j]-1];
+                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tetra_vertex_ids[j]-1)*num_control_type+i];
+                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tetra_vertex_ids[j]-1)*num_control_type+i];
+                                                this->inputs->AddInput(new ControlInput(DamageDEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
+                                        }
+                                        break;
+                                default:
+                                        _error_("Control " << EnumToStringx(reCast<int,IssmDouble>(iomodel->Data(InversionControlParametersEnum)[i])) << " not implemented yet");
+                        }
+                }
+        }
+
+        /*Need to know the type of approximation for this element*/
+        if(iomodel->Data(FlowequationElementEquationEnum)){
+                this->inputs->AddInput(new IntInput(ApproximationEnum,reCast<int>(iomodel->Data(FlowequationElementEquationEnum)[index])));
+        }
+
+        /*DatasetInputs*/
+        if (control_analysis && iomodel->Data(InversionCostFunctionsCoefficientsEnum)) {
+
+                /*Create inputs and add to DataSetInput*/
+                DatasetInput* datasetinput=new DatasetInput(InversionCostFunctionsCoefficientsEnum);
+                for(i=0;i<num_responses;i++){
+                        for(j=0;j<NUMVERTICES;j++)nodeinputs[j]=iomodel->Data(InversionCostFunctionsCoefficientsEnum)[(tetra_vertex_ids[j]-1)*num_responses+i];
+                        datasetinput->AddInput(new TetraInput(InversionCostFunctionsCoefficientsEnum,nodeinputs,P1Enum),reCast<int>(iomodel->Data(InversionCostFunctionsEnum)[i]));
+                }
+
+                /*Add datasetinput to element inputs*/
+                this->inputs->AddInput(datasetinput);
+        }
+}
+/*}}}*/
+/*FUNCTION Tetra::IsIceInElement    THIS ONE{{{*/
 bool   Tetra::IsIceInElement(){
 
-        _error_("not implemented yet");
+        /*Get levelset*/
+        IssmDouble ls[NUMVERTICES];
+        GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
+
+        /*If the level set on one of the nodes is <0, ice is present in this element*/
+        for(int i=0;i<NUMVERTICES;i++){
+                if(ls[i]<0.) return true;
+        }
+
+        return false;
+}
+/*}}}*/
+/*FUNCTION Tetra::IsOnBed {{{*/
+bool Tetra::IsOnBed(){
+        return HasFaceOnBed();
+}
+/*}}}*/
+/*FUNCTION Tetra::IsOnSurface {{{*/
+bool Tetra::IsOnSurface(){
+        return HasFaceOnSurface();
 }
 /*}}}*/
 bool Tetra::IsIcefront(void){/*{{{*/
 
-        _error_("not implemented yet");
+        /*Retrieve all inputs and parameters*/
+        IssmDouble ls[NUMVERTICES];
+        GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
+
+        /* If only one vertex has ice, there is an ice front here */
+        if(IsIceInElement()){
+                int nrice=0;       
+                for(int i=0;i<NUMVERTICES;i++) if(ls[i]<0.) nrice++;
+                if(nrice==1) return true;
+        }
+        return false;
 }/*}}}*/
 /*FUNCTION Tetra::JacobianDeterminant{{{*/
@@ -168 +485 @@
 void Tetra::JacobianDeterminantSurface(IssmDouble* pJdet,IssmDouble* xyz_list,Gauss* gauss){
 
-        _error_("not implemented yet");
+        _assert_(gauss->Enum()==GaussTetraEnum);
+        this->GetJacobianDeterminantFace(pJdet,xyz_list,(GaussTetra*)gauss);
+
+}
+/*}}}*/
+/*FUNCTION Tetra::JacobianDeterminantBase{{{*/
+void Tetra::JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){
+
+        _assert_(gauss->Enum()==GaussTetraEnum);
+        this->GetJacobianDeterminantFace(pJdet,xyz_list_base,(GaussTetra*)gauss);
+
+}
+/*}}}*/
+/*FUNCTION Tetra::JacobianDeterminantTop{{{*/
+void Tetra::JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){
+
+        _assert_(gauss->Enum()==GaussTetraEnum);
+        this->GetJacobianDeterminantFace(pJdet,xyz_list_base,(GaussTetra*)gauss);
 
 }
@@ -182 +516 @@
 }
 /*}}}*/
-/*FUNCTION Tetra::NewElementVector THIS ONE{{{*/
-ElementVector* Tetra::NewElementVector(int approximation_enum){
-        return new ElementVector(nodes,this->NumberofNodes(),this->parameters,approximation_enum);
-}
-/*}}}*/
-/*FUNCTION Tetra::NewElementMatrix  THIS ONE{{{*/
-ElementMatrix* Tetra::NewElementMatrix(int approximation_enum){
-        return new ElementMatrix(nodes,this->NumberofNodes(),this->parameters,approximation_enum);
+/*FUNCTION Tetra::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){{{*/
+Gauss* Tetra::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){
+        /*FIXME: this is messed up, should provide indices and not xyz_list!*/
+        int indices[3];
+        this->FaceOnFrontIndices(&indices[0],&indices[1],&indices[2]);
+        return new GaussTetra(indices[0],indices[1],indices[2],max(order_horiz,order_vert));
+}
+/*}}}*/
+/*FUNCTION Tetra::NewGaussBase(int order){{{*/
+Gauss* Tetra::NewGaussBase(int order){
+
+        int indices[3];
+        this->FaceOnBedIndices(&indices[0],&indices[1],&indices[2]);
+        return new GaussTetra(indices[0],indices[1],indices[2],order);
+}
+/*}}}*/
+/*FUNCTION Tetra::NewGaussTop(int order){{{*/
+Gauss* Tetra::NewGaussTop(int order){
+
+        int indices[3];
+        this->FaceOnSurfaceIndices(&indices[0],&indices[1],&indices[2]);
+        return new GaussTetra(indices[0],indices[1],indices[2],order);
 }
 /*}}}*/
@@ -209 +557 @@
 /*}}}*/
 /*FUNCTION Tetra::NormalSection{{{*/
-void Tetra::NormalSection(IssmDouble* normal,IssmDouble* xyz_list_front){
-
-        _error_("Not implemented yet");
-}
-/*}}}*/
+void Tetra::NormalSection(IssmDouble* normal,IssmDouble* xyz_list){
+
+        /*Build unit outward pointing vector*/
+        IssmDouble AB[3];
+        IssmDouble AC[3];
+        IssmDouble norm;
+
+        AB[0]=xyz_list[1*3+0] - xyz_list[0*3+0];
+        AB[1]=xyz_list[1*3+1] - xyz_list[0*3+1];
+        AB[2]=xyz_list[1*3+2] - xyz_list[0*3+2];
+        AC[0]=xyz_list[2*3+0] - xyz_list[0*3+0];
+        AC[1]=xyz_list[2*3+1] - xyz_list[0*3+1];
+        AC[2]=xyz_list[2*3+2] - xyz_list[0*3+2];
+
+        cross(normal,AB,AC);
+        norm=sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
+
+        for(int i=0;i<3;i++) normal[i]=normal[i]/norm;
+}
+/*}}}*/
+/*FUNCTION Tetra::ReduceMatrices{{{*/
+void Tetra::ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){
+
+        int analysis_type;
+        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+
+        if(pe){
+                if(analysis_type==StressbalanceAnalysisEnum){
+                        if(this->element_type==MINIcondensedEnum){
+                                int approximation;
+                                inputs->GetInputValue(&approximation,ApproximationEnum);
+                                if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
+                                        //Do nothing, condensation already done in PVectorCoupling
+                                }
+                                else{
+                                        _error_("Not implemented");
+                                }
+                        }
+                        else if(this->element_type==P1bubblecondensedEnum){
+                                _error_("Not implemented");
+                        }
+                }
+        }
+
+        if(Ke){
+                if(analysis_type==StressbalanceAnalysisEnum){
+                        int approximation;
+                        inputs->GetInputValue(&approximation,ApproximationEnum);
+                        if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
+                                //Do nothing condensatino already done for Stokes part
+                        }
+                        else{
+                                if(this->element_type==MINIcondensedEnum){
+                                        _error_("Not implemented");
+                                }
+                                else if(this->element_type==P1bubblecondensedEnum){
+                                        _error_("Not implemented");
+                                }
+                        }
+                }
+        }
+}
+/*}}}*/
+/*FUNCTION Tetra::SetCurrentConfiguration {{{*/
+void  Tetra::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, Nodes* nodesin, Materials* materialsin, Parameters* parametersin){
+
+        /*go into parameters and get the analysis_counter: */
+        int analysis_counter;
+        parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);
+
+        /*Get Element type*/
+        this->element_type=this->element_type_list[analysis_counter];
+
+        /*Pick up nodes*/
+        if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
+        else this->nodes=NULL;
+
+}
+/*}}}*/
+/*FUNCTION Tetra::SetwiseNodeConnectivity   THIS ONE (take from Penta.cpp){{{*/
+void Tetra::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){
+
+        /*Intermediaries*/
+        const int numnodes = this->NumberofNodes();
+
+        /*Output */
+        int d_nz = 0;
+        int o_nz = 0;
+
+        /*Loop over all nodes*/
+        for(int i=0;i<numnodes;i++){
+
+                if(!flags[this->nodes[i]->Lid()]){
+
+                        /*flag current node so that no other element processes it*/
+                        flags[this->nodes[i]->Lid()]=true;
+
+                        int counter=0;
+                        while(flagsindices[counter]>=0) counter++;
+                        flagsindices[counter]=this->nodes[i]->Lid();
+
+                        /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
+                        switch(set2_enum){
+                                case FsetEnum:
+                                        if(nodes[i]->indexing.fsize){
+                                                if(this->nodes[i]->IsClone())
+                                                 o_nz += 1;
+                                                else
+                                                 d_nz += 1;
+                                        }
+                                        break;
+                                case GsetEnum:
+                                        if(nodes[i]->indexing.gsize){
+                                                if(this->nodes[i]->IsClone())
+                                                 o_nz += 1;
+                                                else
+                                                 d_nz += 1;
+                                        }
+                                        break;
+                                case SsetEnum:
+                                        if(nodes[i]->indexing.ssize){
+                                                if(this->nodes[i]->IsClone())
+                                                 o_nz += 1;
+                                                else
+                                                 d_nz += 1;
+                                        }
+                                        break;
+                                default: _error_("not supported");
+                        }
+                }
+        }
+
+        /*Assign output pointers: */
+        *pd_nz=d_nz;
+        *po_nz=o_nz;
+}
+/*}}}*/
+/*FUNCTION Tetra::Sid  THIS ONE{{{*/
+int    Tetra::Sid(){
+
+        return sid;
+
+}
+/*}}}*/
+/*FUNCTION Tetra::Update {{{*/
+void Tetra::Update(int index,IoModel* iomodel,int analysis_counter,int analysis_type,int finiteelement_type){ 
+
+        /*Intermediaries*/
+        int        i;
+        int        tetra_vertex_ids[6];
+        IssmDouble nodeinputs[6];
+        IssmDouble yts;
+        bool       dakota_analysis;
+        bool       isFS;
+        int        numnodes;
+        int*       tetra_node_ids = NULL;
+
+        /*Fetch parameters: */
+        iomodel->Constant(&yts,ConstantsYtsEnum);
+        iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
+        iomodel->Constant(&isFS,FlowequationIsFSEnum);
+
+        /*Checks if debuging*/
+        _assert_(iomodel->elements);
+
+        /*Recover element type*/
+        this->SetElementType(finiteelement_type,analysis_counter);
+
+        /*Recover vertices ids needed to initialize inputs*/
+        for(i=0;i<6;i++) tetra_vertex_ids[i]=iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab
+
+        /*Recover nodes ids needed to initialize the node hook.*/
+        switch(finiteelement_type){
+                case P1Enum:
+                        numnodes         = 4;
+                        tetra_node_ids   = xNew<int>(numnodes);
+                        tetra_node_ids[0]=iomodel->nodecounter+iomodel->elements[4*index+0];
+                        tetra_node_ids[1]=iomodel->nodecounter+iomodel->elements[4*index+1];
+                        tetra_node_ids[2]=iomodel->nodecounter+iomodel->elements[4*index+2];
+                        tetra_node_ids[3]=iomodel->nodecounter+iomodel->elements[4*index+3];
+                        break;
+                default:
+                        _error_("Finite element "<<EnumToStringx(finiteelement_type)<<" not supported yet");
+        }
+
+        /*hooks: */
+        this->SetHookNodes(tetra_node_ids,numnodes,analysis_counter); this->nodes=NULL;
+        xDelete<int>(tetra_node_ids);
+
+        /*Fill with IoModel*/
+        this->InputUpdateFromIoModel(index,iomodel);
+}
+/*}}}*/
+/*FUNCTION Tetra::ZeroLevelsetCoordinates{{{*/
+void Tetra::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){
+        /*Compute portion of the element that is grounded*/ 
+
+        IssmDouble  levelset[NUMVERTICES];
+        IssmDouble* xyz_zero = xNew<IssmDouble>(3*3);
+
+        /*Recover parameters and values*/
+        GetInputListOnVertices(&levelset[0],levelsetenum);
+
+        if(levelset[0]==0. && levelset[1]==0. && levelset[2]==0.){ 
+                xyz_zero[3*0+0]=xyz_list[0*3+0];
+                xyz_zero[3*0+1]=xyz_list[0*3+1];
+                xyz_zero[3*0+2]=xyz_list[0*3+2];
+
+                /*New point 2*/
+                xyz_zero[3*1+0]=xyz_list[1*3+0];
+                xyz_zero[3*1+1]=xyz_list[1*3+1];
+                xyz_zero[3*1+2]=xyz_list[1*3+2];
+
+                /*New point 3*/
+                xyz_zero[3*2+0]=xyz_list[2*3+0];
+                xyz_zero[3*2+1]=xyz_list[2*3+1];
+                xyz_zero[3*2+2]=xyz_list[2*3+2];
+        }
+        else if(levelset[0]==0. && levelset[1]==0. && levelset[3]==0.){ 
+                xyz_zero[3*0+0]=xyz_list[0*3+0];
+                xyz_zero[3*0+1]=xyz_list[0*3+1];
+                xyz_zero[3*0+2]=xyz_list[0*3+2];
+
+                /*New point 2*/
+                xyz_zero[3*1+0]=xyz_list[1*3+0];
+                xyz_zero[3*1+1]=xyz_list[1*3+1];
+                xyz_zero[3*1+2]=xyz_list[1*3+2];
+
+                /*New point 3*/
+                xyz_zero[3*2+0]=xyz_list[3*3+0];
+                xyz_zero[3*2+1]=xyz_list[3*3+1];
+                xyz_zero[3*2+2]=xyz_list[3*3+2];
+        }
+        else if(levelset[1]==0. && levelset[2]==0. && levelset[3]==0.){ 
+                xyz_zero[3*0+0]=xyz_list[1*3+0];
+                xyz_zero[3*0+1]=xyz_list[1*3+1];
+                xyz_zero[3*0+2]=xyz_list[1*3+2];
+
+                /*New point 2*/
+                xyz_zero[3*1+0]=xyz_list[2*3+0];
+                xyz_zero[3*1+1]=xyz_list[2*3+1];
+                xyz_zero[3*1+2]=xyz_list[2*3+2];
+
+                /*New point 3*/
+                xyz_zero[3*2+0]=xyz_list[3*3+0];
+                xyz_zero[3*2+1]=xyz_list[3*3+1];
+                xyz_zero[3*2+2]=xyz_list[3*3+2];
+        }
+        else if(levelset[2]==0. && levelset[0]==0. && levelset[3]==0.){ 
+                xyz_zero[3*0+0]=xyz_list[2*3+0];
+                xyz_zero[3*0+1]=xyz_list[2*3+1];
+                xyz_zero[3*0+2]=xyz_list[2*3+2];
+
+                /*New point 2*/
+                xyz_zero[3*1+0]=xyz_list[0*3+0];
+                xyz_zero[3*1+1]=xyz_list[0*3+1];
+                xyz_zero[3*1+2]=xyz_list[0*3+2];
+
+                /*New point 3*/
+                xyz_zero[3*2+0]=xyz_list[3*3+0];
+                xyz_zero[3*2+1]=xyz_list[3*3+1];
+                xyz_zero[3*2+2]=xyz_list[3*3+2];
+        }
+        else _error_("Case not covered");
+
+        /*Assign output pointer*/
+        *pxyz_zero= xyz_zero;
+}
+/*}}}*/

issm/trunk-jpl/src/c/classes/Elements/Tetra.h

@@ -53 +53 @@
                 void  InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type){_error_("not implemented yet");};
 #endif
-                void  InputUpdateFromConstant(IssmDouble constant, int name){_error_("not implemented yet");};
-                void  InputUpdateFromConstant(int constant, int name){_error_("not implemented yet");};
-                void  InputUpdateFromConstant(bool constant, int name){_error_("not implemented yet");};
-                void  InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
+                void  InputUpdateFromIoModel(int index, IoModel* iomodel);
                 /*}}}*/
                 /*Element virtual functions definitions: {{{*/
                 void        AddBasalInput(int input_enum, IssmDouble* values, int interpolation_enum){_error_("not implemented yet");};
-                void        AddInput(int input_enum, IssmDouble* values, int interpolation_enum){_error_("not implemented yet");};
+                void        AddInput(int input_enum, IssmDouble* values, int interpolation_enum);
                 IssmDouble  CharacteristicLength(void){_error_("not implemented yet");};
                 void        ComputeBasalStress(Vector<IssmDouble>* sigma_b){_error_("not implemented yet");};
                 void        ComputeStrainRate(Vector<IssmDouble>* eps){_error_("not implemented yet");};
                 void        ComputeStressTensor(){_error_("not implemented yet");};
-                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters){_error_("not implemented yet");};
-                void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters){_error_("not implemented yet");};
-                void        SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){_error_("not implemented yet");};
+                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
+                void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
+                void        SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
                 void        Delta18oParameterization(void){_error_("not implemented yet");};
                 void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");};
@@ -74 +71 @@
                 IssmDouble  EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");};
                 IssmDouble  EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");};
+                void        FaceOnFrontIndices(int* pindex1,int* pindex2,int* pindex3);
+                void        FaceOnBedIndices(int* pindex1,int* pindex2,int* pindex3);
+                void        FaceOnSurfaceIndices(int* pindex1,int* pindex2,int* pindex3);
                 int         FiniteElement(void);
                 Element*    GetUpperElement(void){_error_("not implemented yet");};
@@ -80 +80 @@
                 Element*    GetBasalElement(void){_error_("not implemented yet");};
                 int         GetNodeIndex(Node* node){_error_("not implemented yet");};
-                void        GetNodesSidList(int* sidlist){_error_("not implemented yet");};
-                void        GetNodesLidList(int* lidlist){_error_("not implemented yet");};
                 int         GetNumberOfNodes(void);
                 int         GetNumberOfNodesVelocity(void){_error_("not implemented yet");};
                 int         GetNumberOfNodesPressure(void){_error_("not implemented yet");};
                 int         GetNumberOfVertices(void);
-                void        GetVerticesCoordinates(IssmDouble** pxyz_list);
-                void        GetVerticesCoordinatesBase(IssmDouble** pxyz_list){_error_("not implemented yet");};
-                void        GetVerticesCoordinatesTop(IssmDouble** pxyz_list){_error_("not implemented yet");};
-                int         Sid(){_error_("not implemented yet");};
-                void        InputChangeName(int input_enum, int enum_type_old){_error_("not implemented yet");};
-                bool        IsOnBed(){_error_("not implemented yet");};
-                bool        IsOnSurface(){_error_("not implemented yet");};
+                void        GetVerticesCoordinatesBase(IssmDouble** pxyz_list);
+                void        GetVerticesCoordinatesTop(IssmDouble** pxyz_list);
+                bool        HasFaceOnBed();
+                bool        HasFaceOnSurface();
+                int         Sid();
+                bool        IsOnBed();
+                bool        IsOnSurface();
                 bool        IsNodeOnShelfFromFlags(IssmDouble* flags){_error_("not implemented yet");};
                 void        JacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
                 void        JacobianDeterminantLine(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
                 void        JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
-                void        JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");};
-                void        JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");};
+                void        JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
+                void        JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
                 IssmDouble  MinEdgeLength(IssmDouble* xyz_list){_error_("not implemented yet");};
                 void        NodalFunctions(IssmDouble* basis,Gauss* gauss);
@@ -128 +126 @@
                 IssmDouble  GetXcoord(Gauss* gauss){_error_("Not implemented");};
                 IssmDouble  GetYcoord(Gauss* gauss){_error_("Not implemented");};
-                IssmDouble  GetZcoord(Gauss* gauss){_error_("not implemented yet");};
+                IssmDouble  GetZcoord(Gauss* gauss);
                 int         GetElementType(void){_error_("not implemented yet");};
                 Gauss*      NewGauss(void);
                 Gauss*      NewGauss(int order);
       Gauss*      NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");};
-      Gauss*      NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){_error_("not implemented yet");};
+      Gauss*      NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert);
       Gauss*      NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order){_error_("not implemented yet");};
-                Gauss*      NewGaussBase(int order){_error_("not implemented yet");};
+                Gauss*      NewGaussBase(int order);
                 Gauss*      NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");};
-                Gauss*      NewGaussTop(int order){_error_("not implemented yet");};
-                ElementVector* NewElementVector(int approximation_enum);
-                ElementMatrix* NewElementMatrix(int approximation_enum);
-                ElementMatrix* NewElementMatrixCoupling(int number_nodes,int approximation_enum){_error_("not implemented yet");};
+                Gauss*      NewGaussTop(int order);
                 int         VertexConnectivity(int vertexindex){_error_("not implemented yet");};
                 void        VerticalSegmentIndices(int** pindices,int* pnumseg){_error_("not implemented yet");};
@@ -146 +141 @@
                 bool        IsZeroLevelset(int levelset_enum){_error_("not implemented");};
                 bool               IsIcefront(void);
-                void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented");};
+                void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
                 void               GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");};
                 void        GetNormalFromLSF(IssmDouble *pnormal){_error_("not implemented yet");};
@@ -152 +147 @@
                 void        GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type){_error_("not implemented yet");};
                 void        GetVectorFromInputs(Vector<IssmDouble>* vector, int name_enum){_error_("not implemented yet");};
-                void        InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){_error_("not implemented yet");};
                 void        InputDepthAverageAtBase(int enum_type,int average_enum_type){_error_("not implemented yet");};
                 void        InputDuplicate(int original_enum,int new_enum){_error_("not implemented yet");};
@@ -160 +154 @@
                 void        PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){_error_("not implemented yet");};
                 void        ResetFSBasalBoundaryCondition(void){_error_("not implemented yet");};
-                void        ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){_error_("not implemented yet");};
+                void        ReduceMatrices(ElementMatrix* Ke,ElementVector* pe);
                 void        SetTemporaryElementType(int element_type_in){_error_("not implemented yet");};
                 void          SmbGradients(){_error_("not implemented yet");};
                 IssmDouble  SurfaceArea(void){_error_("not implemented yet");};
-                void        Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement){_error_("not implemented yet");};
+                void        Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);
                 IssmDouble  TimeAdapt(){_error_("not implemented yet");};
                 void UpdateConstraintsExtrudeFromBase(){_error_("not implemented");};

issm/trunk-jpl/src/c/classes/Elements/TetraRef.cpp

@@ -329 +329 @@
 }
 /*}}}*/
+/*FUNCTION TetraRef::GetJacobianDeterminantFace{{{*/
+void TetraRef::GetJacobianDeterminantFace(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussTetra* gauss){
+        /*The Jacobian determinant is constant over the element, discard the gaussian points. 
+         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
+
+        IssmDouble x1=xyz_list[3*0+0];
+        IssmDouble y1=xyz_list[3*0+1];
+        IssmDouble z1=xyz_list[3*0+2];
+        IssmDouble x2=xyz_list[3*1+0];
+        IssmDouble y2=xyz_list[3*1+1];
+        IssmDouble z2=xyz_list[3*1+2];
+        IssmDouble x3=xyz_list[3*2+0];
+        IssmDouble y3=xyz_list[3*2+1];
+        IssmDouble z3=xyz_list[3*2+2];
+
+        /*Jdet = norm( AB ^ AC ) / (2 * area of the reference triangle), with areaRef=sqrt(3) */
+        *Jdet=SQRT3/6.*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2),0.5);
+        if(*Jdet<0) _error_("negative jacobian determinant!");
+}
+/*}}}*/
 /*FUNCTION TetraRef::GetJacobianInvert {{{*/
 void TetraRef::GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussTetra* gauss){

issm/trunk-jpl/src/c/classes/Elements/TetraRef.h

@@ -25 +25 @@
                 void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussTetra* gauss);
                 void GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussTetra* gauss);
+                void GetJacobianDeterminantFace(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussTetra* gauss);
                 void GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussTetra* gauss);
                 void GetNodalFunctions(IssmDouble* basis,GaussTetra* gauss);

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -26 +26 @@
 /*FUNCTION Tria::Tria(int id, int sid,int index, IoModel* iomodel,int nummodels){{{*/
 Tria::Tria(int tria_id, int tria_sid, int index, IoModel* iomodel,int nummodels)
-        :TriaRef(nummodels),ElementHook(nummodels,index+1,3,iomodel){
+        :TriaRef(nummodels),ElementHook(nummodels,index+1,NUMVERTICES,iomodel){
 
                 /*id: */
@@ -687 +687 @@
 }
 /*}}}*/
-/*FUNCTION Tria::GetVerticesCoordinates(IssmDouble** pxyz_list){{{*/
-void Tria::GetVerticesCoordinates(IssmDouble** pxyz_list){
-
-        IssmDouble* xyz_list = xNew<IssmDouble>(NUMVERTICES*3);
-        ::GetVerticesCoordinates(xyz_list,this->vertices,NUMVERTICES);
-
-        /*Assign output pointer*/
-        *pxyz_list = xyz_list;
-
-}/*}}}*/
 /*FUNCTION Tria::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){{{*/
 void Tria::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){
@@ -927 +917 @@
 }
 /*}}}*/
-/*FUNCTION Tria::GetNodesSidList{{{*/
-void Tria::GetNodesSidList(int* sidlist){
-
-        _assert_(sidlist);
-        _assert_(nodes);
-        int numnodes = this->NumberofNodes();
-
-        for(int i=0;i<numnodes;i++){
-                sidlist[i]=nodes[i]->Sid();
-        }
-}
-/*}}}*/
-/*FUNCTION Tria::GetNodesLidList{{{*/
-void Tria::GetNodesLidList(int* lidlist){
-
-        _assert_(lidlist);
-        _assert_(nodes);
-        int numnodes = this->NumberofNodes();
-
-        for(int i=0;i<numnodes;i++){
-                lidlist[i]=nodes[i]->Lid();
-        }
-}
-/*}}}*/
 /*FUNCTION Tria::GetNumberOfNodes;{{{*/
 int Tria::GetNumberOfNodes(void){
@@ -1086 +1052 @@
 }
 /*}}}*/
-/*FUNCTION Tria::InputChangeName{{{*/
-void  Tria::InputChangeName(int original_enum,int new_enum){
-
-        /*Call inputs method*/
-        this->inputs->ChangeEnum(original_enum,new_enum);
-}
-/*}}}*/
 /*FUNCTION Tria::InputScale{{{*/
 void  Tria::InputScale(int enum_type,IssmDouble scale_factor){
@@ -1104 +1063 @@
         /*Scale: */
         input->Scale(scale_factor);
-}
-/*}}}*/
-/*FUNCTION Tria::InputUpdateFromConstant(int value, int name);{{{*/
-void  Tria::InputUpdateFromConstant(int constant, int name){
-
-        /*Check that name is an element input*/
-        if(!IsInput(name)) return;
-
-        /*update input*/
-        this->inputs->AddInput(new IntInput(name,constant));
-}
-/*}}}*/
-/*FUNCTION Tria::InputUpdateFromConstant(IssmDouble value, int name);{{{*/
-void  Tria::InputUpdateFromConstant(IssmDouble constant, int name){
-
-        /*Check that name is an element input*/
-        if(!IsInput(name)) return;
-
-        /*update input*/
-        this->inputs->AddInput(new DoubleInput(name,constant));
-}
-/*}}}*/
-/*FUNCTION Tria::InputUpdateFromConstant(bool value, int name);{{{*/
-void  Tria::InputUpdateFromConstant(bool constant, int name){
-
-        /*Check that name is an element input*/
-        if(!IsInput(name)) return;
-
-        /*update input*/
-        this->inputs->AddInput(new BoolInput(name,constant));
 }
 /*}}}*/
@@ -1340 +1269 @@
 }
 /*}}}*/
-/*FUNCTION Tria::InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){{{*/
-void Tria::InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){
-
-        /*Intermediaries*/
-        int    i,t;
-        int    tria_vertex_ids[3];
-        int    row;
-        IssmDouble nodeinputs[3];
-        IssmDouble time;
-        TransientInput* transientinput=NULL;
-
-        /*Branch on type of vector: nodal or elementary: */
-        if(vector_type==1){ //nodal vector
-
-                /*Recover vertices ids needed to initialize inputs*/
-                for(i=0;i<3;i++){ 
-                        _assert_(iomodel->elements);
-                        tria_vertex_ids[i]=reCast<int>(iomodel->elements[3*this->sid+i]); //ids for vertices are in the elements array from Matlab
-                }
-
-                /*Are we in transient or static? */
-                if(M==iomodel->numberofvertices){
-
-                        /*create input values: */
-                        for(i=0;i<3;i++)nodeinputs[i]=vector[tria_vertex_ids[i]-1];
-
-                        /*create static input: */
-                        this->inputs->AddInput(new TriaInput(vector_enum,nodeinputs,P1Enum));
-                }
-                else if(M==iomodel->numberofvertices+1){
-                        /*create transient input: */
-                        IssmDouble* times = xNew<IssmDouble>(N);
-                        for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
-                        transientinput=new TransientInput(vector_enum,times,N);
-                        for(t=0;t<N;t++){
-                                for(i=0;i<NUMVERTICES;i++) nodeinputs[i]=vector[N*(tria_vertex_ids[i]-1)+t];
-                                transientinput->AddTimeInput(new TriaInput(vector_enum,nodeinputs,P1Enum));
-                        }
-                        this->inputs->AddInput(transientinput);
-                        xDelete<IssmDouble>(times);
-                }
-                else _error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
-        }
-        else if(vector_type==2){ //element vector
-                /*Are we in transient or static? */
-                if(M==iomodel->numberofelements){
-
-                        /*static mode: create an input out of the element value: */
-
-                        if (code==5){ //boolean
-                                this->inputs->AddInput(new BoolInput(vector_enum,reCast<bool>(vector[this->sid])));
-                        }
-                        else if (code==6){ //integer
-                                this->inputs->AddInput(new IntInput(vector_enum,reCast<int>(vector[this->sid])));
-                        }
-                        else if (code==7){ //IssmDouble
-                                this->inputs->AddInput(new DoubleInput(vector_enum,vector[this->sid]));
-                        }
-                        else _error_("could not recognize nature of vector from code " << code);
-                }
-                else {
-                        _error_("transient element inputs not supported yet!");
-                }
-        }
-        else{
-                _error_("Cannot add input for vector type " << vector_type << " (not supported)");
-        }
-
-}
-/*}}}*/
 /*FUNCTION Tria::IsOnBed {{{*/
 bool Tria::IsOnBed(){
@@ -1681 +1540 @@
         this->EdgeOnSurfaceIndices(&indices[0],&indices[1]);
         return new GaussTria(indices[0],indices[1],order);
-}
-/*}}}*/
-/*FUNCTION Tria::NewElementVector{{{*/
-ElementVector* Tria::NewElementVector(int approximation_enum){
-        return new ElementVector(nodes,this->NumberofNodes(),this->parameters,approximation_enum);
-}
-/*}}}*/
-/*FUNCTION Tria::NewElementMatrix{{{*/
-ElementMatrix* Tria::NewElementMatrix(int approximation_enum){
-        return new ElementMatrix(nodes,this->NumberofNodes(),this->parameters,approximation_enum);
 }
 /*}}}*/
@@ -2024 +1873 @@
 }
 /*}}}*/
-/*FUNCTION Tria::SetCurrentConfiguration {{{*/
-void  Tria::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, Nodes* nodesin, Materials* materialsin, Parameters* parametersin){
-
-        /*go into parameters and get the analysis_counter: */
-        int analysis_counter;
-        parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);
-
-        /*Get Element type*/
-        this->element_type=this->element_type_list[analysis_counter];
-
-        /*Pick up nodes*/
-        if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
-        else this->nodes=NULL;
-
-}
-/*}}}*/
 /*FUNCTION Tria::SpawnSeg {{{*/
 Seg*  Tria::SpawnSeg(int index1,int index2){
@@ -2104 +1937 @@
                         _error_("not implemented yet");
         }
+}
+/*}}}*/
+/*FUNCTION Tria::SetCurrentConfiguration {{{*/
+void  Tria::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, Nodes* nodesin, Materials* materialsin, Parameters* parametersin){
+
+        /*go into parameters and get the analysis_counter: */
+        int analysis_counter;
+        parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);
+
+        /*Get Element type*/
+        this->element_type=this->element_type_list[analysis_counter];
+
+        /*Pick up nodes*/
+        if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
+        else this->nodes=NULL;
+
 }
 /*}}}*/

issm/trunk-jpl/src/c/classes/Elements/Tria.h

@@ -53 +53 @@
                 void  InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type);
                 #endif
-                void  InputUpdateFromConstant(IssmDouble constant, int name);
-                void  InputUpdateFromConstant(int constant, int name);
-                void  InputUpdateFromConstant(bool constant, int name);
                 void  InputUpdateFromIoModel(int index, IoModel* iomodel);
                 /*}}}*/
@@ -79 +76 @@
                 IssmDouble  GetGroundedPortion(IssmDouble* xyz_list);
                 int         GetNodeIndex(Node* node);
-                void        GetNodesSidList(int* sidlist);
-                void        GetNodesLidList(int* lidlist);
                 int         GetNumberOfNodes(void);
                 int         GetNumberOfNodesPressure(void);
@@ -102 +97 @@
                 void        GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type);
                 void        GetVectorFromInputs(Vector<IssmDouble>* vector, int name_enum);
-                void        GetVerticesCoordinates(IssmDouble** pxyz_list);
                 void        GetVerticesCoordinatesBase(IssmDouble** pxyz_list);
                 void        GetVerticesCoordinatesTop(IssmDouble** pxyz_list);
-                void        InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code);
                 void        InputDepthAverageAtBase(int enum_type,int average_enum_type);
                 void        InputDuplicate(int original_enum,int new_enum);
@@ -215 +208 @@
                 void           GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype);
                 Node*          GetNode(int node_number);
-                void             InputChangeName(int enum_type,int enum_type_old);
                 void             InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type);
                 void             InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int enum_type){_error_("not implemented yet");};
@@ -232 +224 @@
                 Gauss*         NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");};
                 Gauss*         NewGaussTop(int order);
-                ElementVector* NewElementVector(int approximation_enum);
-                ElementMatrix* NewElementMatrix(int approximation_enum);
-                ElementMatrix* NewElementMatrixCoupling(int number_nodes,int approximation_enum){_error_("not implemented yet");};
                 void           NodalFunctions(IssmDouble* basis,Gauss* gauss);
                 void           NodalFunctionsP1(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");};

issm/trunk-jpl/src/c/classes/classes.h

@@ -61 +61 @@
 #include "./Inputs/DoubleInput.h"
 #include "./Inputs/IntInput.h"
+#include "./Inputs/TetraInput.h"
 #include "./Inputs/PentaInput.h"
 #include "./Inputs/TriaInput.h"

issm/trunk-jpl/src/c/classes/gauss/GaussTetra.cpp

@@ -3 +3 @@
  */
 
+#include <math.h>
 #include "./GaussTetra.h"
 #include "../../shared/io/Print/Print.h"
@@ -32 +33 @@
 /*FUNCTION GaussTetra::GaussTetra(int order) {{{*/
 GaussTetra::GaussTetra(int order){
-        _error_("not implemented yet");
+        /*Get gauss points*/
+        GaussLegendreTetra(&numgauss,&coords1,&coords2,&coords3,&coords4,&weights,order);
+        for(int i=0;i<numgauss;i++) this->weights[i]=this->weights[i]*sqrt(2.)/6.; //FIXME: double check that
+
+        /*Initialize static fields as undefinite*/
+        weight=UNDEF;
+        coord1=UNDEF;
+        coord2=UNDEF;
+        coord3=UNDEF;
+}
+/*}}}*/
+/*FUNCTION GaussTetra::GaussTetra(int index1,int index2,int index3,int order) {{{*/
+GaussTetra::GaussTetra(int index1,int index2,int index3,int order){
+
+        /*Basal Tria*/
+        if(index1==0 && index2==1 && index3==2){
+                GaussLegendreTria(&numgauss,&coords1,&coords2,&coords3,&weights,order);
+                coords4=xNew<IssmDouble>(numgauss);
+                for(int i=0;i<numgauss;i++) coords4[i]=0.;
+        }
+        else if(index1==0 && index2==1 && index3==3){
+                GaussLegendreTria(&numgauss,&coords1,&coords2,&coords4,&weights,order);
+                coords3=xNew<IssmDouble>(numgauss);
+                for(int i=0;i<numgauss;i++) coords3[i]=0.;
+        }
+        else if(index1==1 && index2==2 && index3==3){
+                GaussLegendreTria(&numgauss,&coords2,&coords3,&coords4,&weights,order);
+                coords1=xNew<IssmDouble>(numgauss);
+                for(int i=0;i<numgauss;i++) coords1[i]=0.;
+        }
+        else if(index1==2 && index2==0 && index3==3){
+                GaussLegendreTria(&numgauss,&coords1,&coords3,&coords4,&weights,order);
+                coords2=xNew<IssmDouble>(numgauss);
+                for(int i=0;i<numgauss;i++) coords2[i]=0.;
+        }
+        else{
+                _error_(index1 <<" "<<index2 <<" "<<index3 <<" Not supported yet");
+        }
 }
 /*}}}*/
@@ -119 +157 @@
         /*update static arrays*/
         switch(iv){
-                default: _error_("vertex index should be in [0 5]");
+                case 0: coord1=1.; coord2=0.; coord3=0.; coord4=0.; break;
+                case 1: coord1=0.; coord2=1.; coord3=0.; coord4=0.; break;
+                case 2: coord1=0.; coord2=0.; coord3=1.; coord4=0.; break;
+                case 3: coord1=0.; coord2=0.; coord3=0.; coord4=1.; break;
+                default: _error_("vertex index should be in [0 3]");
 
         }
@@ -133 +175 @@
         /*update static arrays*/
         switch(finiteelement){
+                case P1Enum: case P1DGEnum:
+                        switch(iv){
+                                case 0: coord1=1.; coord2=0.; coord3=0.; coord4=0.; break;
+                                case 1: coord1=0.; coord2=1.; coord3=0.; coord4=0.; break;
+                                case 2: coord1=0.; coord2=0.; coord3=1.; coord4=0.; break;
+                                case 3: coord1=0.; coord2=0.; coord3=0.; coord4=1.; break;
+                                default: _error_("node index should be in [0 3]");
+                        }
+                        break;
                 default: _error_("Finite element "<<EnumToStringx(finiteelement)<<" not supported");
         }

issm/trunk-jpl/src/c/classes/gauss/GaussTetra.h

@@ -31 +31 @@
                 GaussTetra();
                 GaussTetra(int order);
+                GaussTetra(int index1,int index2,int index3,int order);
                 ~GaussTetra();
 

issm/trunk-jpl/src/c/modules/MeshPartitionx/MeshPartitionx.h

@@ -33 +33 @@
         switch(meshtype){
                 case Mesh2DhorizontalEnum:
-                        epart=xNew<int>(numberofelements);
-                        npart=xNew<int>(numberofnodes);
-                        index=xNew<int>(elements_width*numberofelements);
-                        for (i=0;i<numberofelements;i++){
-                                for (j=0;j<elements_width;j++){
-                                        *(index+elements_width*i+j)=(*(elements+elements_width*i+j))-1; //-1 for C indexing in Metis
-                                }
-                        }
-
-                        /*Partition using Metis:*/
-                        if (num_procs>1){
-#ifdef _HAVE_METIS_
-                                METIS_PartMeshNodalPatch(&numberofelements,&numberofnodes, index, &etype, &numflag, &num_procs, &edgecut, epart, npart);
-#else
-                                _error_("metis has not beed installed. Cannot run with more than 1 cpu");
-#endif
-                        }
-                        else if (num_procs==1){
-                                /*METIS does not know how to deal with one cpu only!*/
-                                for (i=0;i<numberofelements;i++) epart[i]=0;
-                                for (i=0;i<numberofnodes;i++)    npart[i]=0;
-                        }
-                        else _error_("At least one processor is required");
-                        break;
                 case Mesh2DverticalEnum:
+                case Mesh3DtetrasEnum:
                         epart=xNew<int>(numberofelements);
                         npart=xNew<int>(numberofnodes);

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp

@@ -37 +37 @@
                         }
                         break;
+                case Mesh3DtetrasEnum:
+                        for(i=0;i<iomodel->numberofelements;i++){
+                                if(iomodel->my_elements[i]) elements->AddObject(new Tetra(i+1,i,i,iomodel,nummodels));
+                        }
+                        break;
                 case Mesh3DEnum:
                         iomodel->FetchData(2,MeshUpperelementsEnum,MeshLowerelementsEnum);
@@ -60 +65 @@
                                         if(dakota_analysis) elements->InputDuplicate(MaterialsRheologyBbarEnum,QmuMaterialsRheologyBEnum);
                                         break;
-                                case Mesh3DEnum:
+                                case Mesh3DEnum: case Mesh3DtetrasEnum:
                                         if(dakota_analysis) elements->InputDuplicate(MaterialsRheologyBEnum,QmuMaterialsRheologyBEnum); 
                                         break;

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNumberNodeToElementConnectivity.cpp

@@ -39 +39 @@
                 case Mesh2DhorizontalEnum: elementswidth=3; break;
                 case Mesh2DverticalEnum:   elementswidth=3; break;
+                case Mesh3DtetrasEnum:     elementswidth=4; break;
                 case Mesh3DEnum:           elementswidth=6; break;
                 default:                   _error_("mesh not supported yet");

issm/trunk-jpl/src/c/modules/ModelProcessorx/ElementsAndVerticesPartitioning.cpp

@@ -60 +60 @@
                 case Mesh2DverticalEnum:
                         elements_width=3;
+                        numberofelements2d = 0;
+                        numberofvertices2d = 0;
+                        numlayers          = 0;
+                        break;
+                case Mesh3DtetrasEnum:
+                        elements_width=4; //penta elements
                         numberofelements2d = 0;
                         numberofvertices2d = 0;

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -215 +215 @@
         Mesh2DverticalEnum,
         Mesh3DEnum,
+        Mesh3DtetrasEnum,
         MiscellaneousNameEnum, //FIXME: only used by qmu, should not be marshalled (already in queueing script)
         MasstransportHydrostaticAdjustmentEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -223 +223 @@
                 case Mesh2DverticalEnum : return "Mesh2Dvertical";
                 case Mesh3DEnum : return "Mesh3D";
+                case Mesh3DtetrasEnum : return "Mesh3Dtetras";
                 case MiscellaneousNameEnum : return "MiscellaneousName";
                 case MasstransportHydrostaticAdjustmentEnum : return "MasstransportHydrostaticAdjustment";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -226 +226 @@
               else if (strcmp(name,"Mesh2Dvertical")==0) return Mesh2DverticalEnum;
               else if (strcmp(name,"Mesh3D")==0) return Mesh3DEnum;
+              else if (strcmp(name,"Mesh3Dtetras")==0) return Mesh3DtetrasEnum;
               else if (strcmp(name,"MiscellaneousName")==0) return MiscellaneousNameEnum;
               else if (strcmp(name,"MasstransportHydrostaticAdjustment")==0) return MasstransportHydrostaticAdjustmentEnum;
@@ -259 +260 @@
               else if (strcmp(name,"ProfilingSolutionTime")==0) return ProfilingSolutionTimeEnum;
               else if (strcmp(name,"MaxIterationConvergenceFlag")==0) return MaxIterationConvergenceFlagEnum;
-              else if (strcmp(name,"SteadystateMaxiter")==0) return SteadystateMaxiterEnum;
          else stage=3;
    }
    if(stage==3){
-              if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum;
+              if (strcmp(name,"SteadystateMaxiter")==0) return SteadystateMaxiterEnum;
+              else if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum;
               else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum;
               else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum;
@@ -382 +383 @@
               else if (strcmp(name,"Materials")==0) return MaterialsEnum;
               else if (strcmp(name,"Nodes")==0) return NodesEnum;
-              else if (strcmp(name,"Contours")==0) return ContoursEnum;
          else stage=4;
    }
    if(stage==4){
-              if (strcmp(name,"Parameters")==0) return ParametersEnum;
+              if (strcmp(name,"Contours")==0) return ContoursEnum;
+              else if (strcmp(name,"Parameters")==0) return ParametersEnum;
               else if (strcmp(name,"Vertices")==0) return VerticesEnum;
               else if (strcmp(name,"Results")==0) return ResultsEnum;
@@ -505 +506 @@
               else if (strcmp(name,"Vz")==0) return VzEnum;
               else if (strcmp(name,"VzSSA")==0) return VzSSAEnum;
-              else if (strcmp(name,"VzHO")==0) return VzHOEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"VzPicard")==0) return VzPicardEnum;
+              if (strcmp(name,"VzHO")==0) return VzHOEnum;
+              else if (strcmp(name,"VzPicard")==0) return VzPicardEnum;
               else if (strcmp(name,"VzFS")==0) return VzFSEnum;
               else if (strcmp(name,"VxMesh")==0) return VxMeshEnum;
@@ -628 +629 @@
               else if (strcmp(name,"LockFileName")==0) return LockFileNameEnum;
               else if (strcmp(name,"ToolkitsOptionsAnalyses")==0) return ToolkitsOptionsAnalysesEnum;
-              else if (strcmp(name,"ToolkitsOptionsStrings")==0) return ToolkitsOptionsStringsEnum;
          else stage=6;
    }
    if(stage==6){
-              if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum;
+              if (strcmp(name,"ToolkitsOptionsStrings")==0) return ToolkitsOptionsStringsEnum;
+              else if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum;
               else if (strcmp(name,"QmuInName")==0) return QmuInNameEnum;
               else if (strcmp(name,"QmuOutName")==0) return QmuOutNameEnum;