Changeset 3417

Timestamp:

04/07/10 12:23:51 (15 years ago)

Author:

Eric.Larour

Message:

Massive temporary commit: redesign to include vertices (for galerkin discontinous), to simplify ModelProcessor.

Location:

issm/trunk/src

Files:

• c/DataSet/DataSet.h (modified) (1 diff)
• c/Dofx/Dofx.cpp (modified) (4 diffs)
• c/Dofx/Dofx.h (modified) (1 diff)
• c/EnumDefinitions/EnumDefinitions.cpp (modified) (2 diffs)
• c/EnumDefinitions/EnumDefinitions.h (modified) (2 diffs)
• c/Makefile.am (modified) (4 diffs)
• c/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp (modified) (8 diffs)
• c/ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp (modified) (10 diffs)
• c/ModelProcessorx/CreateDataSets.cpp (modified) (12 diffs)
• c/ModelProcessorx/CreateParameters.cpp (modified) (1 diff)
• c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (modified) (6 diffs)
• c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp (modified) (9 diffs)
• c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp (modified) (9 diffs)
• c/ModelProcessorx/IoModel.cpp (modified) (6 diffs)
• c/ModelProcessorx/IoModel.h (modified) (13 diffs)
• c/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp (modified) (9 diffs)
• c/ModelProcessorx/Partitioning.cpp (added)
• c/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp (modified) (9 diffs)
• c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp (modified) (9 diffs)
• c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp (modified) (9 diffs)
• c/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp (modified) (9 diffs)
• c/objects/DofIndexing.cpp (added)
• c/objects/DofIndexing.h (added)
• c/objects/ElementProperties.cpp (modified) (1 diff)
• c/objects/ElementProperties.h (modified) (1 diff)
• c/objects/FemModel.cpp (modified) (9 diffs)
• c/objects/FemModel.h (modified) (6 diffs)
• c/objects/Hook.cpp (modified) (4 diffs)
• c/objects/Hook.h (modified) (2 diffs)
• c/objects/Matice.cpp (modified) (1 diff)
• c/objects/Matice.h (modified) (1 diff)
• c/objects/Matpar.cpp (modified) (1 diff)
• c/objects/Matpar.h (modified) (1 diff)
• c/objects/Model.cpp (modified) (8 diffs)
• c/objects/Node.cpp (modified) (32 diffs)
• c/objects/Node.h (modified) (3 diffs)
• c/objects/NodeProperties.cpp (added)
• c/objects/NodeProperties.h (added)
• c/objects/Penta.cpp (modified) (1 diff)
• c/objects/Penta.h (modified) (2 diffs)
• c/objects/Tria.cpp (modified) (1 diff)
• c/objects/Tria.h (modified) (2 diffs)
• c/objects/Vertex.cpp (modified) (9 diffs)
• c/objects/Vertex.h (modified) (3 diffs)
• m/enum/BeamEnum.m (modified) (1 diff)
• m/enum/DofIndexingEnum.m (added)
• m/enum/NodeEnum.m (modified) (1 diff)
• m/enum/NodePropertiesEnum.m (added)
• m/enum/PentaEnum.m (modified) (1 diff)
• m/enum/SingEnum.m (modified) (1 diff)
• m/enum/VertexEnum.m (added)
• m/enum/VerticesEnum.m (added)
• m/solutions/jpl/CreateFemModel.m (modified) (2 diffs)
• mex/Dof/Dof.cpp (modified) (4 diffs)
• mex/Dof/Dof.h (modified) (1 diff)
• mex/ModelProcessor/ModelProcessor.cpp (modified) (3 diffs)
• mex/ModelProcessor/ModelProcessor.h (modified) (1 diff)

issm/trunk/src/c/DataSet/DataSet.h

@@ -7 +7 @@
 
 
-#ifndef DATASET_H_
-#define DATASET_H_
+#ifndef _DATASET_H_
+#define _DATASET_H_
 
 #include <vector>

issm/trunk/src/c/Dofx/Dofx.cpp

@@ -10 +10 @@
 #include "../EnumDefinitions/EnumDefinitions.h"
 
-int Dofx( DofVec** ppartition, DofVec** ptpartition,DataSet* elements,DataSet* nodes, DataSet* params) {
+int Dofx( DofVec** ppartition, DofVec** ptpartition, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* params) {
 
         int noerr=1;
@@ -21 +21 @@
         /*intermediary: */
         int  numberofnodes;
+        int  numberofvertices;
         int  numberofdofspernode;
 
@@ -31 +32 @@
         tpartition=new DofVec("tpartition");
 
-        /*First, recover number of grids from parameters: */
+        /*First, recover number of vertices and nodes from parameters: */
+        found=params->FindParam(&numberofvertices,"numberofvertices");
+        if(!found)ISSMERROR(" could not find numberofvertices in parameters");
+        
         found=params->FindParam(&numberofnodes,"numberofnodes");
         if(!found)ISSMERROR(" could not find numberofnodes in parameters");
@@ -41 +45 @@
         /*Ensure that only for each cpu, the partition border nodes only will be taken into account once 
          * across the cluster. To do so, we flag all the clone nodes: */
+        vertices->FlagClones(numberofvertices);
         nodes->FlagClones(numberofnodes);
 

issm/trunk/src/c/Dofx/Dofx.h

@@ -9 +9 @@
 
 /* local prototypes: */
-int             Dofx( DofVec** partition, DofVec** ptpartition,DataSet* elements,DataSet* nodesin,DataSet* params);
+int             Dofx( DofVec** partition, DofVec** ptpartition,DataSet* elements,DataSet* nodesin, DataSet* verticesin, DataSet* params);
 
 #endif  /* _DOFX_H */

issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp

@@ -15 +15 @@
 int ParametersEnum(void){                   return          106; }
 int ResultsEnum(void){                      return          107; }
+int VerticesEnum(void){                     return          108; }
 
 /*ANALYSIS TYPES: */
@@ -71 +72 @@
 int ElementEnum(void){                      return          410; }
 int TriaEnum(void){                         return          411; }
-int ElementPropertiesEnum(void){               return          412; }
-int PentaEnum(void){                        return          413; }
-int SingEnum(void){                         return          414; }
-int BeamEnum(void){                         return          415; }
+int ElementPropertiesEnum(void){            return          412; }
+int NodePropertiesEnum(void){               return          413; }
+int PentaEnum(void){                        return          414; }
+int SingEnum(void){                         return          415; }
+int BeamEnum(void){                         return          416; }
+int DofIndexingEnum(void){                  return          417; }
+
 /*Grids: */
-int NodeEnum(void){                         return          420; }
+int VertexEnum(void){                         return        420; }
+int NodeEnum(void){                         return          421; }
 /*Loads: */
 int LoadEnum(void){                         return          430; }

issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h

@@ -16 +16 @@
 int ParametersEnum(void);
 int ResultsEnum(void);
+int VerticesEnum(void);
 
 /*ANALYSIS TYPES: */
@@ -73 +74 @@
 int TriaEnum(void);
 int ElementPropertiesEnum(void);
+int NodePropertiesEnum(void);
 int PentaEnum(void);
 int SingEnum(void);
 int BeamEnum(void);
+int DofIndexingEnum(void);
 /*Grids: */
 int NodeEnum(void);
+int VertexEnum(void);
 /*Loads: */
 int LoadEnum(void);

issm/trunk/src/c/Makefile.am

@@ -52 +52 @@
                                         ./objects/DakotaPlugin.h\
                                         ./objects/DakotaPlugin.cpp\
-                                        ./objects/Node.h\
-                                        ./objects/Node.cpp\
                                         ./objects/Vertex.h\
                                         ./objects/Vertex.cpp\
                                         ./objects/Hook.h\
                                         ./objects/Hook.cpp\
+                                        ./objects/DofIndexing.h\
+                                        ./objects/DofIndexing.cpp\
+                                        ./objects/NodeProperties.h\
+                                        ./objects/NodeProperties.cpp\
+                                        ./objects/Node.h\
+                                        ./objects/Node.cpp\
                                         ./objects/Result.h\
                                         ./objects/Result.cpp\
@@ -201 +205 @@
                                         ./ModelProcessorx/IoModel.h\
                                         ./ModelProcessorx/IoModel.cpp\
+                                        ./ModelProcessorx/Partitioning.cpp\
                                         ./ModelProcessorx/CreateDataSets.cpp\
                                         ./ModelProcessorx/CreateParameters.cpp\
@@ -446 +451 @@
                                         ./objects/DakotaPlugin.h\
                                         ./objects/DakotaPlugin.cpp\
+                                        ./objects/DofIndexing.h\
+                                        ./objects/DofIndexing.cpp\
+                                        ./objects/NodeProperties.h\
+                                        ./objects/NodeProperties.cpp\
                                         ./objects/Node.h\
                                         ./objects/Node.cpp\
@@ -592 +601 @@
                                         ./ModelProcessorx/IoModel.h\
                                         ./ModelProcessorx/IoModel.cpp\
+                                        ./ModelProcessorx/Partitioning.cpp\
                                         ./ModelProcessorx/CreateDataSets.cpp\
                                         ./ModelProcessorx/CreateParameters.cpp\

issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp

@@ -13 +13 @@
 
 
-void    CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+void    CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
 
 
@@ -26 +26 @@
         DataSet*    elements  = NULL;
         DataSet*    nodes = NULL;
+        DataSet*    vertices = NULL;
         DataSet*    materials = NULL;
         
         /*Objects: */
         Node*       node   = NULL;
+        NodeProperties node_properties;
+        Vertex*     vertex = NULL;
         Tria*       tria = NULL;
         Penta*      penta = NULL;
@@ -96 +99 @@
         /* node constructor input: */
         int node_id;
-        int node_partitionborder=0;
-        double node_x[3];
-        double node_sigma;
+        int vertex_id;
+        int partitionborder=0;
         int node_onbed;
         int node_onsurface;
@@ -121 +123 @@
         elements  = new DataSet(ElementsEnum());
         nodes     = new DataSet(NodesEnum());
+        vertices  = new DataSet(VerticesEnum());
         materials = new DataSet(MaterialsEnum());
 
@@ -436 +439 @@
         #endif
 
+                #ifdef _PARALLEL_
+                if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border
+                        partitionborder=1;
+                }
+                else{
+                        partitionborder=0;
+                }
+                #else
+                        partitionborder=0;
+                #endif
+
+                /*create vertex: */
+                vertex_id=i+1;
+                vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]),partitionborder);
+                vertices->AddObject(vertex);
+
+                /*create node: */
                 node_id=i+1; //matlab indexing
                         
                 
-                
-                #ifdef _PARALLEL_
-                if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border
-                        node_partitionborder=1;
-                }
-                else{
-                        node_partitionborder=0;
-                }
-                #else
-                        node_partitionborder=0;
-                #endif
-
-                node_x[0]=iomodel->x[i];
-                node_x[1]=iomodel->y[i];
-                node_x[2]=iomodel->z[i];
-                node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]);
-                
-                node_onbed=(int)iomodel->gridonbed[i];
-                node_onsurface=(int)iomodel->gridonsurface[i];  
-                node_onshelf=(int)iomodel->gridoniceshelf[i];   
-                node_onsheet=(int)iomodel->gridonicesheet[i];   
+                node_properties.SetProperties(
+                                (int)iomodel->gridonbed[i],
+                                (int)iomodel->gridonsurface[i],
+                                (int)iomodel->gridoniceshelf[i],
+                                (int)iomodel->gridonicesheet[i]);
 
                 if (strcmp(iomodel->meshtype,"3d")==0){
@@ -475 +479 @@
 
                 /*Create node using its constructor: */
-                node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
+                node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties);   
+
 
                 /*set single point constraints.: */
@@ -497 +502 @@
          * datasets, it will not be redone: */
         elements->Presort();
+        vertices->Presort();
         nodes->Presort();
         materials->Presort();
@@ -531 +537 @@
         *pelements=elements;
         *pnodes=nodes;
+        *pvertices=vertices;
         *pmaterials=materials;
 

issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp

@@ -13 +13 @@
 
 
-void    CreateElementsNodesAndMaterialsBalancedvelocities(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+void    CreateElementsNodesAndMaterialsBalancedvelocities(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
 
 
@@ -26 +26 @@
         DataSet*    elements  = NULL;
         DataSet*    nodes = NULL;
+        DataSet*    vertices = NULL;
         DataSet*    materials = NULL;
         ElementProperties* tria_properties=NULL;
@@ -32 +33 @@
         /*Objects: */
         Node*       node   = NULL;
+        Vertex*     vertex = NULL;
         Tria*       tria = NULL;
         Penta*      penta = NULL;
@@ -96 +98 @@
         /* node constructor input: */
         int node_id;
-        int node_partitionborder=0;
-        double node_x[3];
-        double node_sigma;
+        int vertex_id;
+        int partitionborder=0;
         int node_onbed;
         int node_onsurface;
@@ -121 +122 @@
         elements  = new DataSet(ElementsEnum());
         nodes     = new DataSet(NodesEnum());
+        vertices  = new DataSet(VerticesEnum());
         materials = new DataSet(MaterialsEnum());
 
@@ -435 +437 @@
         #endif
 
+                /*create vertex: */
+                vertex_id=i+1;
+                vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]));
+
+                vertices->AddObject(vertex);
+
                 node_id=i+1; //matlab indexing
                         
-                
-                
                 #ifdef _PARALLEL_
                 if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border
@@ -450 +456 @@
                 #endif
 
-                node_x[0]=iomodel->x[i];
-                node_x[1]=iomodel->y[i];
-                node_x[2]=iomodel->z[i];
-                node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]);
-                
-                node_onbed=(int)iomodel->gridonbed[i];
-                node_onsurface=(int)iomodel->gridonsurface[i];  
-                node_onshelf=(int)iomodel->gridoniceshelf[i];   
-                node_onsheet=(int)iomodel->gridonicesheet[i];   
+                node_properties.SetProperties(
+                        (int)iomodel->gridonbed[i],
+                        (int)iomodel->gridonsurface[i],
+                        (int)iomodel->gridoniceshelf[i],
+                        (int)iomodel->gridonicesheet[i]);
+
 
                 if (strcmp(iomodel->meshtype,"3d")==0){
@@ -474 +477 @@
 
                 /*Create node using its constructor: */
-                node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
+                node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties);   
 
                 /*set single point constraints.: */
@@ -497 +500 @@
         elements->Presort();
         nodes->Presort();
+        vertices->Presort();
         materials->Presort();
 
@@ -530 +534 @@
         *pelements=elements;
         *pnodes=nodes;
+        *pvertices=vertices;
         *pmaterials=materials;
 

issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp

@@ -15 +15 @@
 
 
-void CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle){
+void CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle){
 
         /*create parameters common to all solutions: */
@@ -26 +26 @@
                 if (iomodel->sub_analysis_type==HorizAnalysisEnum()){
 
-                        CreateElementsNodesAndMaterialsDiagnosticHoriz(pelements,pnodes,pmaterials, iomodel,iomodel_handle);
+                        CreateElementsNodesAndMaterialsDiagnosticHoriz(pelements,pnodes, pvertices, pmaterials, iomodel,iomodel_handle);
                         CreateConstraintsDiagnosticHoriz(pconstraints,iomodel,iomodel_handle);
                         CreateLoadsDiagnosticHoriz(ploads,iomodel,iomodel_handle);
@@ -34 +34 @@
                 else if (iomodel->sub_analysis_type==VertAnalysisEnum()){
                 
-                        CreateElementsNodesAndMaterialsDiagnosticVert(pelements,pnodes,pmaterials, iomodel,iomodel_handle);
+                        CreateElementsNodesAndMaterialsDiagnosticVert(pelements,pnodes, pvertices, pmaterials, iomodel,iomodel_handle);
                         CreateConstraintsDiagnosticVert(pconstraints,iomodel,iomodel_handle);
                         CreateLoadsDiagnosticVert(ploads,iomodel,iomodel_handle);
@@ -41 +41 @@
                 else if (iomodel->sub_analysis_type==StokesAnalysisEnum()){
 
-                        CreateElementsNodesAndMaterialsDiagnosticStokes(pelements,pnodes,pmaterials, iomodel,iomodel_handle);
+                        CreateElementsNodesAndMaterialsDiagnosticStokes(pelements,pnodes, pvertices, pmaterials, iomodel,iomodel_handle);
                         CreateConstraintsDiagnosticStokes(pconstraints,iomodel,iomodel_handle);
                         CreateLoadsDiagnosticStokes(ploads,iomodel,iomodel_handle);
@@ -48 +48 @@
                 else if (iomodel->sub_analysis_type==HutterAnalysisEnum()){
 
-                        CreateElementsNodesAndMaterialsDiagnosticHutter(pelements,pnodes,pmaterials, iomodel,iomodel_handle);
+                        CreateElementsNodesAndMaterialsDiagnosticHutter(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle);
                         CreateConstraintsDiagnosticHutter(pconstraints,iomodel,iomodel_handle);
                         CreateLoadsDiagnosticHutter(ploads,iomodel,iomodel_handle);
@@ -56 +56 @@
         else if (iomodel->analysis_type==SlopecomputeAnalysisEnum()){
 
-                CreateElementsNodesAndMaterialsSlopeCompute(pelements,pnodes,pmaterials, iomodel,iomodel_handle);
+                CreateElementsNodesAndMaterialsSlopeCompute(pelements,pnodes, pvertices,pmaterials, iomodel,iomodel_handle);
                 CreateConstraintsSlopeCompute(pconstraints,iomodel,iomodel_handle);
                 CreateLoadsSlopeCompute(ploads,iomodel,iomodel_handle);
@@ -63 +63 @@
         else if (iomodel->analysis_type==ThermalAnalysisEnum()){
 
-                CreateElementsNodesAndMaterialsThermal(pelements,pnodes,pmaterials, iomodel,iomodel_handle);
+                CreateElementsNodesAndMaterialsThermal(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle);
                 CreateConstraintsThermal(pconstraints,iomodel,iomodel_handle);
                 CreateLoadsThermal(ploads,iomodel,iomodel_handle);
@@ -71 +71 @@
         else if (iomodel->analysis_type==MeltingAnalysisEnum()){
                         
-                CreateElementsNodesAndMaterialsMelting(pelements,pnodes,pmaterials, iomodel,iomodel_handle);
+                CreateElementsNodesAndMaterialsMelting(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle);
                 CreateConstraintsMelting(pconstraints,iomodel,iomodel_handle);
                 CreateLoadsMelting(ploads,iomodel,iomodel_handle);
@@ -78 +78 @@
         else if (iomodel->analysis_type==PrognosticAnalysisEnum()){
 
-                CreateElementsNodesAndMaterialsPrognostic(pelements,pnodes,pmaterials, iomodel,iomodel_handle);
+                CreateElementsNodesAndMaterialsPrognostic(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle);
                 CreateConstraintsPrognostic(pconstraints,iomodel,iomodel_handle);
                 CreateLoadsPrognostic(ploads,iomodel,iomodel_handle);
@@ -86 +86 @@
         else if (iomodel->analysis_type==Prognostic2AnalysisEnum()){
 
-                CreateElementsNodesAndMaterialsPrognostic2(pelements,pnodes,pmaterials, iomodel,iomodel_handle);
+                CreateElementsNodesAndMaterialsPrognostic2(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle);
                 CreateConstraintsPrognostic2(pconstraints,iomodel,iomodel_handle);
                 CreateLoadsPrognostic2(ploads,iomodel,iomodel_handle);
@@ -94 +94 @@
         else if (iomodel->analysis_type==BalancedthicknessAnalysisEnum()){
 
-                CreateElementsNodesAndMaterialsBalancedthickness(pelements,pnodes,pmaterials, iomodel,iomodel_handle);
+                CreateElementsNodesAndMaterialsBalancedthickness(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle);
                 CreateConstraintsBalancedthickness(pconstraints,iomodel,iomodel_handle);
                 CreateLoadsBalancedthickness(ploads,iomodel,iomodel_handle);
@@ -102 +102 @@
         else if (iomodel->analysis_type==BalancedvelocitiesAnalysisEnum()){
 
-                CreateElementsNodesAndMaterialsBalancedvelocities(pelements,pnodes,pmaterials, iomodel,iomodel_handle);
+                CreateElementsNodesAndMaterialsBalancedvelocities(pelements,pnodes,pvertices,pmaterials, iomodel,iomodel_handle);
                 CreateConstraintsBalancedvelocities(pconstraints,iomodel,iomodel_handle);
                 CreateLoadsBalancedvelocities(ploads,iomodel,iomodel_handle);

issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp

@@ -238 +238 @@
         parameters->AddObject(param);
 
-        /*numberofgrids: */
+        /*numberofvertices: */
+        count++;
+        param= new Param(count,"numberofvertices",DOUBLE);
+        param->SetDouble(iomodel->numberofvertices);
+        parameters->AddObject(param);
+
+        /*numberofnodes: */
         count++;
         param= new Param(count,"numberofnodes",DOUBLE);
         if (iomodel->analysis_type==Prognostic2AnalysisEnum()) param->SetDouble(3*iomodel->numberofelements);
-        else param->SetDouble(iomodel->numberofnodes);
+        else param->SetDouble(iomodel->numberofvertices);
         parameters->AddObject(param);
 

issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp

@@ -13 +13 @@
 
 
-void    CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
-
-
-        /*output: int* epart, int* my_grids, double* my_bordergrids*/
+void    CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
 
 
         int i,j,k,n;
-        extern int my_rank;
-        extern int num_procs;
 
         /*DataSets: */
         DataSet*    elements  = NULL;
         DataSet*    nodes = NULL;
+        DataSet*    vertices = NULL;
         DataSet*    materials = NULL;
         
         /*Objects: */
         Node*       node   = NULL;
+        Vertex*     vertex = NULL;
         Tria*       tria = NULL;
         Penta*      penta = NULL;
         Matice*     matice  = NULL;
         Matpar*     matpar  = NULL;
-        ElementProperties* tria_properties=NULL;
-        ElementProperties* penta_properties=NULL;
-
-        /*output: */
-        int* epart=NULL; //element partitioning.
-        int* npart=NULL; //node partitioning.
-        int* my_grids=NULL;
-        double* my_bordergrids=NULL;
-
-
-        /*intermediary: */
-        int elements_width; //size of elements
-        double B_avg;
-                        
-        /*tria constructor input: */
-        int tria_id;
-        int tria_matice_id;
-        int tria_matpar_id;
-        int tria_numpar_id;
-        int tria_node_ids[3];
-        double tria_h[3];
-        double tria_s[3];
-        double tria_b[3];
-        double tria_k[3];
-        double tria_melting[3];
-        double tria_accumulation[3];
-        int    tria_friction_type;
-        double tria_p;
-        double tria_q;
-        int    tria_shelf;
-        bool   tria_onwater; 
-        
-        /*matice constructor input: */
-        int    matice_mid;
-        double matice_B;
-        double matice_n;
-        
-        /*penta constructor input: */
-        int penta_id;
-        int penta_matice_id;
-        int penta_matpar_id;
-        int penta_numpar_id;
-        int penta_node_ids[6];
-        double penta_h[6];
-        double penta_s[6];
-        double penta_b[6];
-        double penta_k[6];
-        int penta_friction_type;
-        double penta_p;
-        double penta_q;
-        int penta_shelf;
-        int penta_onbed;
-        int penta_onsurface;
-        int penta_collapse;
-        double penta_melting[6];
-        double penta_accumulation[6];
-        bool   penta_onwater;
-
-        /*matpar constructor input: */
-        int     matpar_mid;
-        double matpar_rho_ice; 
-        double matpar_rho_water;
-        double matpar_heatcapacity;
-        double matpar_thermalconductivity;
-        double matpar_latentheat;
-        double matpar_beta;
-        double matpar_meltingpoint;
-        double matpar_mixed_layer_capacity;
-        double matpar_thermal_exchange_velocity;
-        double matpar_g;
-
-        /* node constructor input: */
-        int node_id;
-        int node_partitionborder=0;
-        double node_x[3];
-        double node_sigma;
-        int node_onbed;
-        int node_onsurface;
-        int node_onsheet;
-        int node_onshelf;
-        int node_upper_node_id;
-        int node_numdofs;
-
-                
-        #ifdef _PARALLEL_
-        /*Metis partitioning: */
-        int  range;
-        Vec  gridborder=NULL;
-        int  my_numgrids;
-        int* all_numgrids=NULL;
-        int  gridcount;
-        int  count;
-        #endif
-        int  first_grid_index;
-
-        /*Rifts:*/
-        int el1,el2;
-         
+
+        /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */
+        if (!iomodel->ismacayealpattyn)goto cleanup_and_return;
+
         /*First create the elements, nodes and material properties: */
         elements  = new DataSet(ElementsEnum());
         nodes     = new DataSet(NodesEnum());
+        vertices  = new DataSet(VerticesEnum());
         materials = new DataSet(MaterialsEnum());
 
-
-        /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */
-        if (!iomodel->ismacayealpattyn)goto cleanup_and_return;
-
-        /*Width of elements: */
-        if(strcmp(iomodel->meshtype,"2d")==0){
-                elements_width=3; //tria elements
-        }
-        else{
-                elements_width=6; //penta elements
-        }
-
-
-
-        #ifdef _PARALLEL_
-        /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/
-        if(strcmp(iomodel->meshtype,"2d")==0){
-                /*load elements: */
-                IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
-        }
-        else{
-                /*load elements2d: */
-                IoModelFetchData(&iomodel->elements2d,NULL,NULL,iomodel_handle,"elements2d");
-        }
-
-        MeshPartitionx(&epart, &npart,iomodel->numberofelements,iomodel->numberofnodes,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofnodes2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs);
-
-        /*Free elements and elements2d: */
-        xfree((void**)&iomodel->elements);
-        xfree((void**)&iomodel->elements2d);
-
-
-        /*Deal with rifts, they have to be included into one partition only, not several: */
-        IoModelFetchData(&iomodel->riftinfo,&iomodel->numrifts,NULL,iomodel_handle,"riftinfo");
-
-        for(i=0;i<iomodel->numrifts;i++){
-                el1=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+2)-1; //matlab indexing to c indexing
-                el2=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+3)-1; //matlab indexing to c indexing
-                epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids;
-        }
-        /*Free rifts: */
-        xfree((void**)&iomodel->riftinfo); 
-
-        /*Used later on: */
-        my_grids=(int*)xcalloc(iomodel->numberofnodes,sizeof(int));
-        #endif
-
-
-
-        /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */
+        /*Partition elements and vertices and nodes: */
+        Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle);
 
         /*2d mesh: */
@@ -209 +65 @@
                 for (i=0;i<iomodel->numberofelements;i++){
 
-                #ifdef _PARALLEL_
-                /*!All elements have been partitioned above, only create elements for this CPU: */
-                if(my_rank==epart[i]){ 
-                #endif
-                        
-                        if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum()){ //elements of type 1 are Hutter type Tria. Don't create this elements.
+                        if(iomodel->my_elements[i]){
                                 
-                                /*ids: */
-                                tria_id=i+1; //matlab indexing.
-                                tria_matice_id=i+1; //refers to the corresponding material property card
-                                tria_matpar_id=iomodel->numberofelements+1;//refers to the corresponding parmat property card
-                                tria_numpar_id=1;
-
-                                /*vertices ids: */
-                                tria_node_ids[0]=(int)*(iomodel->elements+elements_width*i+0);
-                                tria_node_ids[1]=(int)*(iomodel->elements+elements_width*i+1);
-                                tria_node_ids[2]=(int)*(iomodel->elements+elements_width*i+2);
-
-                                /*thickness,surface and bed:*/
-                                tria_h[0]= *(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+0)-1)); //remember, elements is an index of vertices offsets, in matlab indexing.
-                                tria_h[1]=*(iomodel->thickness+  ((int)*(iomodel->elements+elements_width*i+1)-1)); 
-                                tria_h[2]=*(iomodel->thickness+  ((int)*(iomodel->elements+elements_width*i+2)-1)) ;
-
-                                tria_s[0]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+0)-1)); 
-                                tria_s[1]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+1)-1)); 
-                                tria_s[2]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+2)-1)); 
-
-                                tria_b[0]=*(iomodel->bed+        ((int)*(iomodel->elements+elements_width*i+0)-1)); 
-                                tria_b[1]=*(iomodel->bed+        ((int)*(iomodel->elements+elements_width*i+1)-1)); 
-                                tria_b[2]=*(iomodel->bed+        ((int)*(iomodel->elements+elements_width*i+2)-1)); 
-
-                                /*basal drag:*/
-                                tria_friction_type=(int)iomodel->drag_type;
-
-                                tria_k[0]=*(iomodel->drag+        ((int)*(iomodel->elements+elements_width*i+0)-1)); 
-                                tria_k[1]=*(iomodel->drag+        ((int)*(iomodel->elements+elements_width*i+1)-1)); 
-                                tria_k[2]=*(iomodel->drag+        ((int)*(iomodel->elements+elements_width*i+2)-1)); 
-                                
-                                tria_p=iomodel->p[i];
-                                tria_q=iomodel->q[i];
-
-                                /*meling and accumulation*/
-                                tria_melting[0]=*(iomodel->melting+        ((int)*(iomodel->elements+elements_width*i+0)-1));
-                                tria_melting[1]=*(iomodel->melting+        ((int)*(iomodel->elements+elements_width*i+1)-1));
-                                tria_melting[2]=*(iomodel->melting+        ((int)*(iomodel->elements+elements_width*i+2)-1));
-
-                                tria_accumulation[0]=*(iomodel->accumulation+        ((int)*(iomodel->elements+elements_width*i+0)-1));
-                                tria_accumulation[1]=*(iomodel->accumulation+        ((int)*(iomodel->elements+elements_width*i+1)-1));
-                                tria_accumulation[2]=*(iomodel->accumulation+        ((int)*(iomodel->elements+elements_width*i+2)-1));
-
-                                /*element on iceshelf, water?:*/
-                                tria_shelf=(int)*(iomodel->elementoniceshelf+i);
-                                tria_onwater=(bool)*(iomodel->elementonwater+i);
-
-                                /*Create properties: */
-                                tria_properties=new ElementProperties(3,tria_h, tria_s, tria_b, tria_k, tria_melting, tria_accumulation, NULL,
-                                                tria_friction_type, tria_p, tria_q, tria_shelf, UNDEF,tria_onwater, UNDEF,UNDEF,UNDEF);
-
-                                /*Create tria element using its constructor:*/
-                                tria=new Tria(tria_id, tria_node_ids, tria_matice_id, tria_matpar_id, tria_numpar_id, tria_properties);
-
-                                /*delete properties: */
-                                delete tria_properties;
-
-                                /*Add tria element to elements dataset: */
-                                elements->AddObject(tria);
-
-                                /*Deal with material property card: */
-                                matice_mid=i+1; //same as the material id from the geom2 elements.
-                                
-                                /*Average B over 3 grid elements: */
-                                B_avg=0;
-                                for(j=0;j<3;j++){
-                                        B_avg+=*(iomodel->B+((int)*(iomodel->elements+elements_width*i+j)-1));
+                                if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum()){ //elements of type 1 are Hutter type Tria. Don't create this elements.
+
+                                        /*Create and add tria element to elements dataset: */
+                                        elements->AddObject(new Tria(i,iomodel));
+                                        
+                                        /*Create and add material property to materials dataset: */
+                                        materials->AddObject(new Matice(i,iomodel,3));
                                 }
-                                B_avg=B_avg/3;
-                                matice_B=B_avg;
-                                matice_n=(double)*(iomodel->n+i);
-                                
-                                /*Create matice using its constructor:*/
-                                matice= new Matice(matice_mid,matice_B,matice_n);
-                
-                                /*Add matice element to materials dataset: */
-                                materials->AddObject(matice);
-                
-                        } //if(!MacAyealFormulationEnum)
-
-                        #ifdef _PARALLEL_
-                        /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use 
-                         *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) 
-                         into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids 
-                         will hold which grids belong to this partition*/
-                        my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1;
-                        my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1;
-                        my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1;
-                        #endif
-
-                #ifdef _PARALLEL_
-                }//if(my_rank==epart[i])
-                #endif
-
+                        }
                 }//for (i=0;i<numberofelements;i++)
-
         
                 /*Free data : */
@@ -349 +114 @@
                 
                 for (i=0;i<iomodel->numberofelements;i++){
-                #ifdef _PARALLEL_
-                /*We are using our element partition to decide which elements will be created on this node: */
-                if(my_rank==epart[i]){
-                #endif
-
-                        if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum() | *(iomodel->elements_type+2*i+0)==PattynFormulationEnum()){ //elements of type 1 are Hutter type Tria. Don't create this elements.
-                        
-                                /*name and id: */
-                                penta_id=i+1; //matlab indexing.
-                                penta_matice_id=i+1; //refers to the corresponding material property card
-                                penta_matpar_id=iomodel->numberofelements+1;//refers to the corresponding parmat property card
-                                penta_numpar_id=1;
-
-                                /*vertices,thickness,surface,bed and drag: */
-                                for(j=0;j<6;j++){
-                                        penta_node_ids[j]=(int)*(iomodel->elements+elements_width*i+j);
-                                        penta_h[j]=*(iomodel->thickness+    ((int)*(iomodel->elements+elements_width*i+j)-1)); 
-                                        penta_s[j]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+j)-1)); 
-                                        penta_b[j]=*(iomodel->bed+    ((int)*(iomodel->elements+elements_width*i+j)-1)); 
-                                        penta_k[j]=*(iomodel->drag+        ((int)*(iomodel->elements+elements_width*i+j)-1)); 
-                                        penta_melting[j]=*(iomodel->melting+        ((int)*(iomodel->elements+elements_width*i+j)-1));
-                                        penta_accumulation[j]=*(iomodel->accumulation+        ((int)*(iomodel->elements+elements_width*i+j)-1));
+                        if(my_elements[i]){
+                                if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum() | *(iomodel->elements_type+2*i+0)==PattynFormulationEnum()){ //elements of type 1 are Hutter type Tria. Don't create this elements.
+                                        /*Create and add penta element to elements dataset: */
+                                        elements->AddObject(new Penta(i,iomodel));
+                                        
+                                        /*Create and add material property to materials dataset: */
+                                        materials->AddObject(new Matice(i,iomodel,6));
+                                        
                                 }
-
-                                /*basal drag:*/
-                                penta_friction_type=(int)iomodel->drag_type;
-                
-                                penta_p=iomodel->p[i];
-                                penta_q=iomodel->q[i];
-
-                                /*diverse: */
-                                penta_shelf=(int)*(iomodel->elementoniceshelf+i);
-                                penta_onbed=(int)*(iomodel->elementonbed+i);
-                                penta_onsurface=(int)*(iomodel->elementonsurface+i);
-                                penta_onwater=(bool)*(iomodel->elementonwater+i);
-                                
-                                if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum()){ //elements of type 3 are MacAyeal type Penta. We collapse the formulation on their base.
-                                        penta_collapse=1;
-                                }
-                                else{
-                                        penta_collapse=0;
-                                }
-
-                                /*Create element properties: */
-                                penta_properties=new ElementProperties(6,penta_h, penta_s, penta_b, penta_k, penta_melting, penta_accumulation, NULL, penta_friction_type, penta_p, penta_q, penta_shelf, penta_onbed, penta_onwater, penta_onsurface, penta_collapse, UNDEF);
-
-                                /*Create Penta using its constructor:*/
-                                penta= new Penta(penta_id,penta_node_ids, penta_matice_id, penta_matpar_id, penta_numpar_id, penta_properties);
-
-                                /*delete properties: */
-                                delete penta_properties;
-
-                                /*Add penta element to elements dataset: */
-                                elements->AddObject(penta);
-                
-
-                                /*Deal with material:*/
-                                matice_mid=i+1; //same as the material id from the geom2 elements.
-                                /*Average B over 6 element grids: */
-                                B_avg=0;
-                                for(j=0;j<6;j++){
-                                        B_avg+=*(iomodel->B+((int)*(iomodel->elements+elements_width*i+j)-1));
-                                }
-                                B_avg=B_avg/6;
-                                matice_B= B_avg;
-                                matice_n=(double)*(iomodel->n+i);
-                
-                                /*Create matice using its constructor:*/
-                                matice= new Matice(matice_mid,matice_B,matice_n);
-                
-                                /*Add matice element to materials dataset: */
-                                materials->AddObject(matice);
-                                
-                        }
-
-                        #ifdef _PARALLEL_
-                        /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use 
-                         *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) 
-                         into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids 
-                         will hold which grids belong to this partition*/
-                        my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1;
-                        my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1;
-                        my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1;
-                        my_grids[(int)*(iomodel->elements+elements_width*i+3)-1]=1;
-                        my_grids[(int)*(iomodel->elements+elements_width*i+4)-1]=1;
-                        my_grids[(int)*(iomodel->elements+elements_width*i+5)-1]=1;
-                        #endif
-
-                #ifdef _PARALLEL_
-                }//if(my_rank==epart[i])
-                #endif
+                        }//if(my_elements[i])
 
                 }//for (i=0;i<numberofelements;i++)
@@ -463 +147 @@
         } //if (strcmp(meshtype,"2d")==0)
 
-        /*Add one constant material property to materials: */
-        matpar_mid=iomodel->numberofelements+1; //put it at the end of the materials
-        matpar_g=iomodel->g; 
-        matpar_rho_ice=iomodel->rho_ice; 
-        matpar_rho_water=iomodel->rho_water; 
-        matpar_thermalconductivity=iomodel->thermalconductivity; 
-        matpar_heatcapacity=iomodel->heatcapacity; 
-        matpar_latentheat=iomodel->latentheat; 
-        matpar_beta=iomodel->beta; 
-        matpar_meltingpoint=iomodel->meltingpoint; 
-        matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity; 
-        matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity; 
-
-        /*Create matpar object using its constructor: */
-        matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity,
-                        matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity,
-                        matpar_thermal_exchange_velocity,matpar_g);
-                
-        /*Add to materials datset: */
-        materials->AddObject(matpar);
-        
-        #ifdef _PARALLEL_
-                /*From the element partitioning, we can determine which grids are on the inside of this cpu's 
-                 *element partition, and which are on its border with other nodes:*/
-                gridborder=NewVec(iomodel->numberofnodes);
-
-                for (i=0;i<iomodel->numberofnodes;i++){
-                        if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES);
-                }
-                VecAssemblyBegin(gridborder);
-                VecAssemblyEnd(gridborder);
-
-                #ifdef _ISSM_DEBUG_
-                VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD);
-                #endif
-                
-                VecToMPISerial(&my_bordergrids,gridborder);
-
-                #ifdef _ISSM_DEBUG_
-                if(my_rank==0){
-                        for (i=0;i<iomodel->numberofnodes;i++){
-                                printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]);
-                        }
-                }
-                #endif
-        #endif
-
-        /*Partition penalties in 3d: */
-        if(strcmp(iomodel->meshtype,"3d")==0){
-        
-                /*Get penalties: */
-                IoModelFetchData(&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties");
-
-                if(iomodel->numpenalties){
-
-                        iomodel->penaltypartitioning=(int*)xmalloc(iomodel->numpenalties*sizeof(int));
-                        #ifdef _SERIAL_
-                        for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=1;
-                        #else
-                        for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=-1;
-
-                        for(i=0;i<iomodel->numpenalties;i++){
-                                first_grid_index=(int)(*(iomodel->penalties+i*iomodel->numlayers+0)-1);
-                                if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition  grids
-                                        /*All grids that are being penalised belong to this node's internal grid partition.:*/
-                                        iomodel->penaltypartitioning[i]=1;
-                                }
-                                if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border
-                                        iomodel->penaltypartitioning[i]=0;
-                                }
-                        }
-                        #endif
-                }
-
-                /*Free penalties: */
-                xfree((void**)&iomodel->penalties);
-        }
-
-        /*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids. 
-         We can therefore determine  which grids are internal to this node's partition 
-         and which ones are shared with other nodes because they are on the border of this node's partition. Knowing 
-         that, go and create the grids*/
-
-        /*Create nodes from x,y,z, as well as the spc values on those grids: */
-                
-        /*First fetch data: */
-        if (strcmp(iomodel->meshtype,"3d")==0){
-                IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids");
-                IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids");
-        }
+        
+        /*Add new constrant material property tgo materials, at the end: */
+        materials->AddObject(new Matpar(iomodel));
+        
+        /*Create nodes and vertices: */
         IoModelFetchData(&iomodel->x,NULL,NULL,iomodel_handle,"x");
         IoModelFetchData(&iomodel->y,NULL,NULL,iomodel_handle,"y");
@@ -563 +162 @@
         IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet");
         IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf");
-
-        /*Get number of dofs per node: */
-        DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type);
-
-        for (i=0;i<iomodel->numberofnodes;i++){
-        #ifdef _PARALLEL_
-
-        /*keep only this partition's nodes:*/
-        if((my_grids[i]==1)){
-        #endif
-
-                node_id=i+1; //matlab indexing
+        if (strcmp(iomodel->meshtype,"3d")==0){
+                IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids");
+                IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids");
+        }
+
+        
+        for (i=0;i<iomodel->numberovertices;i++){
+
+                /*vertices and nodes (same number, as we are running continuous galerkin formulation: */
+                if(iomodel->my_vertices[i]){
                         
-                
-                #ifdef _PARALLEL_
-                if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border
-                        node_partitionborder=1;
+                        /*Add vertex to vertices dataset: */
+                        vertices->AddObject(new Vertex(i,iomodel));
+
+                        /*Add node to nodes dataset: */
+                        nodes->AddObject(new Node(i,iomodel));
+
                 }
-                else{
-                        node_partitionborder=0;
-                }
-                #else
-                        node_partitionborder=0;
-                #endif
-
-                node_x[0]=iomodel->x[i];
-                node_x[1]=iomodel->y[i];
-                node_x[2]=iomodel->z[i];
-                node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]);
-                
-                node_onbed=(int)iomodel->gridonbed[i];
-                node_onsurface=(int)iomodel->gridonsurface[i];  
-                node_onshelf=(int)iomodel->gridoniceshelf[i];   
-                node_onsheet=(int)iomodel->gridonicesheet[i];   
-
-                if (strcmp(iomodel->meshtype,"3d")==0){
-                        if (isnan(iomodel->uppernodes[i])){
-                                node_upper_node_id=node_id;  //nodes on surface do not have upper nodes, only themselves.
-                        }
-                        else{
-                                node_upper_node_id=(int)iomodel->uppernodes[i];
-                        }
-                }
-                else{
-                        /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/
-                        node_upper_node_id=node_id;
-                }
-
-                /*Create node using its constructor: */
-                node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
-
-                /*set single point constraints.: */
-                if (strcmp(iomodel->meshtype,"3d")==0){
-                        /*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */
-                        if (iomodel->deadgrids[i]){
-                                for(k=1;k<=node_numdofs;k++){
-                                        node->FreezeDof(k);
-                                }
-                        }
-                }
-                if (iomodel->gridonhutter[i]){
-                        for(k=1;k<=node_numdofs;k++){
-                                node->FreezeDof(k);
-                        }
-                }
-                /*Add node to nodes dataset: */
-                nodes->AddObject(node);
-
-        #ifdef _PARALLEL_
-        } //if((my_grids[i]==1))
-        #endif
         }
-
-        /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 
-         * datasets, it will not be redone: */
-        elements->Presort();
-        nodes->Presort();
-        materials->Presort();
 
         /*Clean fetched data: */
@@ -655 +195 @@
         xfree((void**)&iomodel->gridonicesheet);
         xfree((void**)&iomodel->gridoniceshelf);
-        
-
-        /*Keep partitioning information into iomodel*/
-        iomodel->epart=epart;
-        iomodel->my_grids=my_grids;
-        iomodel->my_bordergrids=my_bordergrids;
-
-        /*Free ressources:*/
-        #ifdef _PARALLEL_
-        xfree((void**)&all_numgrids);
-        xfree((void**)&npart);
-        VecFree(&gridborder);
-        #endif
-
-        cleanup_and_return:
-
+
+        /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 
+         * datasets, it will not be redone: */
+        elements->Presort();
+        nodes->Presort();
+        vertices->Presort();
+        materials->Presort();
+        
         /*Assign output pointer: */
         *pelements=elements;
         *pnodes=nodes;
+        *pvertices=vertices;
         *pmaterials=materials;
 

issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp

@@ -13 +13 @@
 #include "../IoModel.h"
 
-void    CreateElementsNodesAndMaterialsDiagnosticStokes(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+void    CreateElementsNodesAndMaterialsDiagnosticStokes(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
 
         
@@ -26 +26 @@
         DataSet*    elements  = NULL;
         DataSet*    nodes = NULL;
+        DataSet*    vertices = NULL;
         DataSet*    materials = NULL;
         
         /*Objects: */
         Node*       node   = NULL;
+        Vertex*     vertex = NULL;
         Penta*      penta = NULL;
         Matice*     matice  = NULL;
@@ -88 +90 @@
         /* node constructor input: */
         int node_id;
+        int vertex_id;
         int node_partitionborder=0;
-        double node_x[3];
-        double node_sigma;
         int node_onbed;
         int node_onsurface;
@@ -113 +114 @@
         elements  = new DataSet(ElementsEnum());
         nodes     = new DataSet(NodesEnum());
+        vertices  = new DataSet(VerticesEnum());
         materials = new DataSet(MaterialsEnum());
 
@@ -400 +402 @@
         #endif
 
+                /*create vertex: */
+                vertex_id=i+1;
+                vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]));
+
+                vertices->AddObject(vertex);
+
                 node_id=i+1; //matlab indexing
-                        
-                
                 
                 #ifdef _PARALLEL_
@@ -415 +421 @@
                 #endif
 
-                node_x[0]=iomodel->x[i];
-                node_x[1]=iomodel->y[i];
-                node_x[2]=iomodel->z[i];
-                node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]);
-                
-                node_onbed=(int)iomodel->gridonbed[i];
-                node_onsurface=(int)iomodel->gridonsurface[i];  
-                node_onshelf=(int)iomodel->gridoniceshelf[i];   
-                node_onsheet=(int)iomodel->gridonicesheet[i];   
-
+                node_properties.SetProperties(
+                        (int)iomodel->gridonbed[i],
+                        (int)iomodel->gridonsurface[i],
+                        (int)iomodel->gridoniceshelf[i],
+                        (int)iomodel->gridonicesheet[i]);
+
+        
                 if (strcmp(iomodel->meshtype,"3d")==0){
                         if (isnan(iomodel->uppernodes[i])){
@@ -439 +442 @@
 
                 /*Create node using its constructor: */
-                node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
+                node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties);   
 
                 /*set single point constraints.: */
@@ -468 +471 @@
         elements->Presort();
         nodes->Presort();
+        vertices->Presort();
         materials->Presort();
 
@@ -503 +507 @@
         *pelements=elements;
         *pnodes=nodes;
+        *pvertices=vertices;
         *pmaterials=materials;
 

issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp

@@ -14 +14 @@
 
 
-void    CreateElementsNodesAndMaterialsDiagnosticVert(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+void    CreateElementsNodesAndMaterialsDiagnosticVert(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
 
 
@@ -27 +27 @@
         DataSet*    elements  = NULL;
         DataSet*    nodes = NULL;
+        DataSet*    vertices = NULL;
         DataSet*    materials = NULL;
         
         /*Objects: */
         Node*       node   = NULL;
+        Vertex*     vertex = NULL;
         Penta*      penta = NULL;
         Matice*     matice  = NULL;
@@ -88 +90 @@
         /* node constructor input: */
         int node_id;
+        int vertex_id;
         int node_partitionborder=0;
-        double node_x[3];
-        double node_sigma;
         int node_onbed;
         int node_onsurface;
@@ -123 +124 @@
         elements  = new DataSet(ElementsEnum());
         nodes     = new DataSet(NodesEnum());
+        vertices  = new DataSet(VerticesEnum());
         materials = new DataSet(MaterialsEnum());
 
@@ -329 +331 @@
         #endif
 
+                /*create vertex: */
+                vertex_id=i+1;
+                vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]));
+
+                vertices->AddObject(vertex);
+
+
                 node_id=i+1; //matlab indexing
-                        
-                
                 
                 #ifdef _PARALLEL_
@@ -344 +351 @@
                 #endif
 
-                node_x[0]=iomodel->x[i];
-                node_x[1]=iomodel->y[i];
-                node_x[2]=iomodel->z[i];
-                node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]);
-                
-                node_onbed=(int)iomodel->gridonbed[i];
-                node_onsurface=(int)iomodel->gridonsurface[i];
-                node_onshelf=(int)iomodel->gridoniceshelf[i];   
-                node_onsheet=(int)iomodel->gridonicesheet[i];   
-        
+                node_properties.SetProperties(
+                        (int)iomodel->gridonbed[i],
+                        (int)iomodel->gridonsurface[i],
+                        (int)iomodel->gridoniceshelf[i],
+                        (int)iomodel->gridonicesheet[i]);
+
                 if (isnan(iomodel->uppernodes[i])){
                         node_upper_node_id=node_id;  //nodes on surface do not have upper nodes, only themselves.
@@ -362 +365 @@
 
                 /*Create node using its constructor: */
-                node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
+                node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties);   
 
                 /*Add node to nodes dataset: */
@@ -376 +379 @@
         elements->Presort();
         nodes->Presort();
+        vertices->Presort();
         materials->Presort();
 
@@ -409 +413 @@
         *pelements=elements;
         *pnodes=nodes;
+        *pvertices=vertices;
         *pmaterials=materials;
 

issm/trunk/src/c/ModelProcessorx/IoModel.cpp

@@ -46 +46 @@
         iomodel->qmu_npart=0; 
         iomodel->numberofelements=0;
-        iomodel->numberofnodes=0;
+        iomodel->numberofvertices=0;
         iomodel->x=NULL; 
         iomodel->y=NULL;
@@ -52 +52 @@
         iomodel->elements=NULL;
         iomodel->elements_type=NULL;
-        iomodel->numberofnodes2d=0;
+        iomodel->numberofvertices2d=0;
         iomodel->elements2d=NULL;
         iomodel->deadgrids=NULL;
@@ -180 +180 @@
         iomodel->isstokes=0;
 
-
-        iomodel->epart=NULL;
-        iomodel->npart=NULL;
-        iomodel->my_grids=NULL;
-        iomodel->my_bordergrids=NULL;
+        /*exterior data: */
+        iomodel->my_elements=NULL;
+        iomodel->my_vertices=NULL;
+        iomodel->my_nodes=NULL;
+        iomodel->my_bordervertices=NULL;
+        iomodel->penaltypartitioning=NULL;
 
         return iomodel;
@@ -276 +277 @@
         xfree((void**)&iomodel->control_type);
         
-        xfree((void**)&iomodel->epart);
-        xfree((void**)&iomodel->npart);
-        xfree((void**)&iomodel->my_grids);
-        xfree((void**)&iomodel->my_bordergrids);
-        
+        /*exterior data: */
+        xfree((void**)&iomodel->my_elements);
+        xfree((void**)&iomodel->my_vertices);
+        xfree((void**)&iomodel->my_nodes);
+        xfree((void**)&iomodel->my_bordervertices);
+        xfree((void**)&iomodel->penaltypartitioning);
+
         /*!Delete entire structure: */
         xfree((void**)piomodel);
@@ -309 +312 @@
         IoModelFetchData(&iomodel->control_analysis,iomodel_handle,"control_analysis"); 
         IoModelFetchData(&iomodel->meshtype,iomodel_handle,"type");
-        /*!Get numberofelements and numberofnodes: */
-        IoModelFetchData(&iomodel->numberofnodes,iomodel_handle,"numberofgrids");
+        /*!Get numberofelements and numberofvertices: */
+        IoModelFetchData(&iomodel->numberofvertices,iomodel_handle,"numberofgrids");
         IoModelFetchData(&iomodel->numberofelements,iomodel_handle,"numberofelements");
         /*!In case we are running 3d, we are going to need the collapsed and non-collapsed 2d meshes, from which the 3d mesh was extruded: */
@@ -317 +320 @@
                 /*!Deal with 2d mesh: */
                 IoModelFetchData(&iomodel->numberofelements2d,iomodel_handle,"numberofelements2d");
-                IoModelFetchData(&iomodel->numberofnodes2d,iomodel_handle,"numberofgrids2d");
+                IoModelFetchData(&iomodel->numberofvertices2d,iomodel_handle,"numberofgrids2d");
                 IoModelFetchData(&iomodel->numlayers,iomodel_handle,"numlayers");
         }

issm/trunk/src/c/ModelProcessorx/IoModel.h

@@ -28 +28 @@
         /*2d mesh: */
         int     numberofelements;
-        int     numberofnodes;
+        int     numberofvertices;
         double* x;
         double* y;
@@ -36 +36 @@
 
         /*3d mesh: */
-        int     numberofnodes2d;
+        int     numberofvertices2d;
         int     numberofelements2d;
         double* elements2d;
@@ -176 +176 @@
         int      numoutput;
 
-        /*exterior data: */
-        int* epart; /*!element partitioning.*/
-        int* npart; /*!node partitioning.*/
-        int* my_grids; /*! grids that belong to this cpu*/
-        double* my_bordergrids; /*! grids that belong to this cpu, and some other cpu also*/
+        /*exterior partitioning data, to be carried around: */
+        bool*   my_elements;
+        bool*   my_vertices;
+        bool*   my_nodes;
+        double* my_bordervertices;
         int*    penaltypartitioning;
 
@@ -197 +197 @@
 
         /*Creation of fem datasets: general drivers*/
-        void  CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
+        void  CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
         void  CreateParameters(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
 
@@ -204 +204 @@
         
         /*diagnostic horizontal*/
-        void    CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+        void    CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
         void    CreateConstraintsDiagnosticHoriz(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
         void  CreateLoadsDiagnosticHoriz(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
@@ -210 +210 @@
 
         /*diagnostic vertical*/
-        void    CreateElementsNodesAndMaterialsDiagnosticVert(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+        void    CreateElementsNodesAndMaterialsDiagnosticVert(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
         void    CreateConstraintsDiagnosticVert(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
         void  CreateLoadsDiagnosticVert(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
 
         /*diagnostic hutter*/
-        void    CreateElementsNodesAndMaterialsDiagnosticHutter(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+        void    CreateElementsNodesAndMaterialsDiagnosticHutter(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
         void    CreateConstraintsDiagnosticHutter(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
         void  CreateLoadsDiagnosticHutter(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
 
         /*diagnostic stokes*/
-        void    CreateElementsNodesAndMaterialsDiagnosticStokes(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+        void    CreateElementsNodesAndMaterialsDiagnosticStokes(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
         void    CreateConstraintsDiagnosticStokes(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
         void  CreateLoadsDiagnosticStokes(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
 
         /*slope compute*/
-        void    CreateElementsNodesAndMaterialsSlopeCompute(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+        void    CreateElementsNodesAndMaterialsSlopeCompute(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
         void    CreateConstraintsSlopeCompute(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
         void  CreateLoadsSlopeCompute(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
@@ -233 +233 @@
 
         /*thermal:*/
-        void    CreateElementsNodesAndMaterialsThermal(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+        void    CreateElementsNodesAndMaterialsThermal(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
         void    CreateConstraintsThermal(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
         void  CreateLoadsThermal(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
@@ -239 +239 @@
 
         /*melting:*/
-        void    CreateElementsNodesAndMaterialsMelting(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+        void    CreateElementsNodesAndMaterialsMelting(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
         void    CreateConstraintsMelting(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
         void  CreateLoadsMelting(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
@@ -245 +245 @@
 
         /*prognostic:*/
-        void    CreateElementsNodesAndMaterialsPrognostic(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+        void    CreateElementsNodesAndMaterialsPrognostic(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
         void    CreateConstraintsPrognostic(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
         void  CreateLoadsPrognostic(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
@@ -251 +251 @@
 
         /*prognostic2:*/
-        void    CreateElementsNodesAndMaterialsPrognostic2(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+        void    CreateElementsNodesAndMaterialsPrognostic2(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
         void    CreateConstraintsPrognostic2(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
         void  CreateLoadsPrognostic2(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
@@ -257 +257 @@
 
         /*balancedthickness:*/
-        void    CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+        void    CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
         void    CreateConstraintsBalancedthickness(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
         void  CreateLoadsBalancedthickness(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
@@ -263 +263 @@
 
         /*balancedvelocities:*/
-        void    CreateElementsNodesAndMaterialsBalancedvelocities(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+        void    CreateElementsNodesAndMaterialsBalancedvelocities(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
         void    CreateConstraintsBalancedvelocities(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
         void  CreateLoadsBalancedvelocities(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
@@ -271 +271 @@
         void CreateParametersQmu(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
 
+        
+        /*partitioning: */
+        void  Partitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_nodes, double** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle);
+
 #endif  /* _IOMODEL_H */

issm/trunk/src/c/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp

@@ -13 +13 @@
 
 
-void    CreateElementsNodesAndMaterialsMelting(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+void    CreateElementsNodesAndMaterialsMelting(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
 
 
@@ -26 +26 @@
         DataSet*    elements  = NULL;
         DataSet*    nodes = NULL;
+        DataSet*    vertices = NULL;
         DataSet*    materials = NULL;
         
         /*Objects: */
         Node*       node   = NULL;
+        Vertex*     vertex = NULL;
         Penta*      penta = NULL;
         Matice*     matice  = NULL;
@@ -88 +90 @@
         /* node constructor input: */
         int node_id;
+        int vertex_id;
         int node_partitionborder=0;
-        double node_x[3];
-        double node_sigma;
         int node_onbed;
         int node_onsurface;
@@ -113 +114 @@
         elements  = new DataSet(ElementsEnum());
         nodes     = new DataSet(NodesEnum());
+        vertices  = new DataSet(VerticesEnum());
         materials = new DataSet(MaterialsEnum());
 
@@ -371 +373 @@
         #endif
 
+                /*create vertex: */
+                vertex_id=i+1;
+                vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]));
+                vertices->AddObject(vertex);
+
+
                 node_id=i+1; //matlab indexing
-                        
-                
                 
                 #ifdef _PARALLEL_
@@ -386 +392 @@
                 #endif
 
-                node_x[0]=iomodel->x[i];
-                node_x[1]=iomodel->y[i];
-                node_x[2]=iomodel->z[i];
-                node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]);
-                
-                node_onbed=(int)iomodel->gridonbed[i];
-                node_onsurface=(int)iomodel->gridonsurface[i];  
-                node_onshelf=(int)iomodel->gridoniceshelf[i];   
-                node_onsheet=(int)iomodel->gridonicesheet[i];   
+                node_properties.SetProperties(
+                        (int)iomodel->gridonbed[i],
+                        (int)iomodel->gridonsurface[i],
+                        (int)iomodel->gridoniceshelf[i],
+                        (int)iomodel->gridonicesheet[i]);
 
                 if (strcmp(iomodel->meshtype,"3d")==0){
@@ -410 +412 @@
 
                 /*Create node using its constructor: */
-                node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
+                node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties);   
 
                 /*set single point constraints.: */
@@ -430 +432 @@
         elements->Presort();
         nodes->Presort();
+        vertices->Presort();
         materials->Presort();
 
@@ -462 +465 @@
         *pelements=elements;
         *pnodes=nodes;
+        *pvertices=vertices;
         *pmaterials=materials;
 }

issm/trunk/src/c/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp

@@ -14 +14 @@
 
 
-void    CreateElementsNodesAndMaterialsPrognostic(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+void    CreateElementsNodesAndMaterialsPrognostic(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
 
 
@@ -27 +27 @@
         DataSet*    elements  = NULL;
         DataSet*    nodes = NULL;
+        DataSet*    vertices = NULL;
         DataSet*    materials = NULL;
         
         /*Objects: */
         Node*       node   = NULL;
+        Vertex*     vertex = NULL;
         Tria*       tria = NULL;
         Penta*      penta = NULL;
@@ -97 +99 @@
         /* node constructor input: */
         int node_id;
+        int vertex_id;
         int node_partitionborder=0;
-        double node_x[3];
-        double node_sigma;
         int node_onbed;
         int node_onsurface;
@@ -122 +123 @@
         elements  = new DataSet(ElementsEnum());
         nodes     = new DataSet(NodesEnum());
+        vertices  = new DataSet(VerticesEnum());
         materials = new DataSet(MaterialsEnum());
 
@@ -436 +438 @@
         #endif
 
+                /*create vertex: */
+                vertex_id=i+1;
+                vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]));
+                vertices->AddObject(vertex);
+
+
                 node_id=i+1; //matlab indexing
-                        
-                
                 
                 #ifdef _PARALLEL_
@@ -451 +457 @@
                 #endif
 
-                node_x[0]=iomodel->x[i];
-                node_x[1]=iomodel->y[i];
-                node_x[2]=iomodel->z[i];
-                node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]);
-                
-                node_onbed=(int)iomodel->gridonbed[i];
-                node_onsurface=(int)iomodel->gridonsurface[i];  
-                node_onshelf=(int)iomodel->gridoniceshelf[i];   
-                node_onsheet=(int)iomodel->gridonicesheet[i];   
-
+                node_properties.SetProperties(
+                        (int)iomodel->gridonbed[i],
+                        (int)iomodel->gridonsurface[i],
+                        (int)iomodel->gridoniceshelf[i],
+                        (int)iomodel->gridonicesheet[i]);
+
+        
                 if (strcmp(iomodel->meshtype,"3d")==0){
                         if (isnan(iomodel->uppernodes[i])){
@@ -475 +478 @@
 
                 /*Create node using its constructor: */
-                node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
+                node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties);   
 
                 /*set single point constraints.: */
@@ -498 +501 @@
         elements->Presort();
         nodes->Presort();
+        vertices->Presort();
         materials->Presort();
 
@@ -531 +535 @@
         *pelements=elements;
         *pnodes=nodes;
+        *pvertices=vertices;
         *pmaterials=materials;
 

issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp

@@ -13 +13 @@
 #include "../IoModel.h"
 
-void    CreateElementsNodesAndMaterialsPrognostic2(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+void    CreateElementsNodesAndMaterialsPrognostic2(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
 
         int i,j,k,n;
@@ -22 +22 @@
         DataSet* elements  = NULL;
         DataSet* nodes = NULL;
+        DataSet*    vertices = NULL;
         DataSet* materials = NULL;
         
         /*Objects: */
         Node*   node   = NULL;
+        Vertex*     vertex = NULL;
         Tria*   tria = NULL;
         Penta*  penta = NULL;
@@ -96 +98 @@
         /* node constructor input: */
         int node_id;
+        int vertex_id;
+        int vertex_partitionborder=0;
         int node_partitionborder=0;
-        double node_x[3];
-        double node_sigma;
         int node_onbed;
         int node_onsurface;
@@ -125 +127 @@
         elements  = new DataSet(ElementsEnum());
         nodes     = new DataSet(NodesEnum());
+        vertices  = new DataSet(VerticesEnum());
         materials = new DataSet(MaterialsEnum());
 
@@ -349 +352 @@
         DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type);
 
+        /*Create vertices: */
+        for (i=0;i<iomodel->numberofnodes;i++){
+        #ifdef _PARALLEL_
+        /*keep only this partition's nodes:*/
+        if((my_grids[i]==1)){
+        #endif
+
+        /*create vertex: */
+        vertex_id=i+1;
+        vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]));
+
+        vertices->AddObject(vertex);
+
+        #ifdef _PARALLEL_
+        } //if((my_grids[i]==1))
+        #endif
+        }
+
         /*Build Nodes dataset -> 3 for each element*/
         for (i=0;i<iomodel->numberofelements;i++){
@@ -374 +395 @@
                                         node_partitionborder=0;
                                 #endif
-                                node_x[0]=iomodel->x[k];
-                                node_x[1]=iomodel->y[k];
-                                node_x[2]=iomodel->z[k];
-                                node_sigma=(iomodel->z[k]-iomodel->bed[k])/(iomodel->thickness[k]);
-                                node_onbed=(int)iomodel->gridonbed[k];
-                                node_onsurface=(int)iomodel->gridonsurface[k];  
-                                node_onshelf=(int)iomodel->gridoniceshelf[k];   
-                                node_onsheet=(int)iomodel->gridonicesheet[k];   
+
+                                node_properties.SetProperties(
+                                        (int)iomodel->gridonbed[k],
+                                        (int)iomodel->gridonsurface[k],
+                                        (int)iomodel->gridoniceshelf[k],
+                                        (int)iomodel->gridonicesheet[k]);
+
                                 if (strcmp(iomodel->meshtype,"3d")==0){
                                         if (isnan(iomodel->uppernodes[k])){
@@ -396 +416 @@
 
                                 /*Create node using its constructor: */
-                                node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
+                                node=new Node(node_id, k, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties);   
 
                                 /*Add node to nodes dataset: */
@@ -410 +430 @@
         elements->Presort();
         nodes->Presort();
+        vertices->Presort();
         materials->Presort();
 
@@ -443 +464 @@
         *pelements=elements;
         *pnodes=nodes;
+        *pvertices=vertices;
         *pmaterials=materials;
 }

issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp

@@ -13 +13 @@
 #include "../IoModel.h"
 
-void    CreateElementsNodesAndMaterialsSlopeCompute(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+void    CreateElementsNodesAndMaterialsSlopeCompute(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
 
 
@@ -26 +26 @@
         DataSet*    elements  = NULL;
         DataSet*    nodes = NULL;
+        DataSet*    vertices = NULL;
         DataSet*    materials = NULL;
         
         /*Objects: */
         Node*       node   = NULL;
+        Vertex*     vertex = NULL;
         Tria*       tria = NULL;
         Penta*      penta = NULL;
@@ -68 +70 @@
         /* node constructor input: */
         int node_id;
+        int vertex_id;
         int node_partitionborder=0;
-        double node_x[3];
-        double node_sigma;
         int node_onbed;
         int node_onsurface;
@@ -93 +94 @@
         elements  = new DataSet(ElementsEnum());
         nodes     = new DataSet(NodesEnum());
+        vertices  = new DataSet(VerticesEnum());
         materials = new DataSet(MaterialsEnum());
 
@@ -372 +374 @@
         #endif
 
+                /*create vertex: */
+                vertex_id=i+1;
+                vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]));
+                vertices->AddObject(vertex);
+
+
                 node_id=i+1; //matlab indexing
-                        
-                
                 
                 #ifdef _PARALLEL_
@@ -387 +393 @@
                 #endif
 
-                node_x[0]=iomodel->x[i];
-                node_x[1]=iomodel->y[i];
-                node_x[2]=iomodel->z[i];
-                node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]);
-                
-                node_onbed=(int)iomodel->gridonbed[i];
-                node_onsurface=(int)iomodel->gridonsurface[i];  
-                node_onshelf=(int)iomodel->gridoniceshelf[i];   
-                node_onsheet=(int)iomodel->gridonicesheet[i];   
+                node_properties.SetProperties(
+                        (int)iomodel->gridonbed[i],
+                        (int)iomodel->gridonsurface[i],
+                        (int)iomodel->gridoniceshelf[i],
+                        (int)iomodel->gridonicesheet[i]);
 
                 if (strcmp(iomodel->meshtype,"3d")==0){
@@ -411 +413 @@
 
                 /*Create node using its constructor: */
-                node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
+                node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties);   
 
                 /*set single point constraints.: */
@@ -434 +436 @@
         elements->Presort();
         nodes->Presort();
+        vertices->Presort();
         materials->Presort();
 
@@ -467 +470 @@
         *pelements=elements;
         *pnodes=nodes;
+        *pvertices=vertices;
         *pmaterials=materials;
 

issm/trunk/src/c/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp

@@ -12 +12 @@
 #include "../IoModel.h"
 
-void    CreateElementsNodesAndMaterialsThermal(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+void    CreateElementsNodesAndMaterialsThermal(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
 
 
@@ -25 +25 @@
         DataSet*    elements  = NULL;
         DataSet*    nodes = NULL;
+        DataSet*    vertices = NULL;
         DataSet*    materials = NULL;
         
         /*Objects: */
         Node*       node   = NULL;
+        Vertex*     vertex = NULL;
         Penta*      penta = NULL;
         Matice*     matice  = NULL;
@@ -86 +88 @@
         /* node constructor input: */
         int node_id;
+        int vertex_id;
         int node_partitionborder=0;
-        double node_x[3];
-        double node_sigma;
         int node_onbed;
         int node_onsurface;
@@ -111 +112 @@
         elements  = new DataSet(ElementsEnum());
         nodes     = new DataSet(NodesEnum());
+        vertices  = new DataSet(VerticesEnum());
         materials = new DataSet(MaterialsEnum());
 
@@ -369 +371 @@
         #endif
 
+                /*create vertex: */
+                vertex_id=i+1;
+                vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]));
+                vertices->AddObject(vertex);
+
+
                 node_id=i+1; //matlab indexing
                 
@@ -382 +390 @@
                 #endif
 
-                node_x[0]=iomodel->x[i];
-                node_x[1]=iomodel->y[i];
-                node_x[2]=iomodel->z[i];
-                node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]);
-                
-                node_onbed=(int)iomodel->gridonbed[i];
-                node_onsurface=(int)iomodel->gridonsurface[i];  
-                node_onshelf=(int)iomodel->gridoniceshelf[i];   
-                node_onsheet=(int)iomodel->gridonicesheet[i];   
+                node_properties.SetProperties(
+                        (int)iomodel->gridonbed[i],
+                        (int)iomodel->gridonsurface[i],
+                        (int)iomodel->gridoniceshelf[i],
+                        (int)iomodel->gridonicesheet[i]);
+
 
                 if (strcmp(iomodel->meshtype,"3d")==0){
@@ -406 +411 @@
 
                 /*Create node using its constructor: */
-                node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
+                node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties);   
 
                 /*Add node to nodes dataset: */
@@ -420 +425 @@
         elements->Presort();
         nodes->Presort();
+        vertices->Presort();
         materials->Presort();
 
@@ -440 +446 @@
         *pelements=elements;
         *pnodes=nodes;
+        *pvertices=vertices;
         *pmaterials=materials;
 

issm/trunk/src/c/objects/ElementProperties.cpp

@@ -96 +96 @@
 }
 /*}}}*/
-/*FUNCTION ElementProperties copy constructor{{{1*/
-ElementProperties::ElementProperties(ElementProperties* prop){ //copy constructor
-        
+/*FUNCTION ElementProperties constructor {{{1*/
+
+ElementProperties::ElementProperties(int elementproperties_numnodes, double* elementproperties_h,double* elementproperties_s,double* elementproperties_b,
+                double* elementproperties_k,double* elementproperties_melting,double* elementproperties_accumulation,
+                double* elementproperties_geothermalflux, int elementproperties_friction_type,double elementproperties_p,
+                double elementproperties_q, int elementproperties_shelf, int elementproperties_onbed, bool elementproperties_onwater, 
+                int elementproperties_onsurface, int elementproperties_collapse, int elementproperties_thermal_steadystate){
+
+        this->Init(elementproperties_numnodes, elementproperties_h,elementproperties_s,elementproperties_b,
+                elementproperties_k,elementproperties_melting,elementproperties_accumulation,
+                elementproperties_geothermalflux, elementproperties_friction_type,elementproperties_p,
+                elementproperties_q, elementproperties_shelf, elementproperties_onbed, elementproperties_onwater,
+                elementproperties_onsurface, elementproperties_collapse, elementproperties_thermal_steadystate);
+}
+/*}}}*/
+/*FUNCTION ElementProperties Initialize propreties, used by constructor{{{1*/
+ElementProperties::Init(int elementproperties_numnodes, double* elementproperties_h,double* elementproperties_s,double* elementproperties_b,
+                double* elementproperties_k,double* elementproperties_melting,double* elementproperties_accumulation,
+                double* elementproperties_geothermalflux, int elementproperties_friction_type,double elementproperties_p,
+                double elementproperties_q, int elementproperties_shelf, int elementproperties_onbed, bool elementproperties_onwater, 
+                int elementproperties_onsurface, int elementproperties_collapse, int elementproperties_thermal_steadystate){
+
         int i;
 

issm/trunk/src/c/objects/ElementProperties.h

@@ -33 +33 @@
                 ElementProperties();
                 ElementProperties(int numnodes, double* h, double* s, double* b, double* k, double* melting, double* accumulation, double* geothermalflux, int friction_type, double p, double q, int shelf, int onbed, bool onwater, int onsurface, int collapse, int thermal_steadystate);
+                Init(int numnodes, double* h, double* s, double* b, double* k, double* melting, double* accumulation, double* geothermalflux, int friction_type, double p, double q, int shelf, int onbed, bool onwater, int onsurface, int collapse, int thermal_steadystate);
                 ElementProperties(ElementProperties* properties);
                 ~ElementProperties();

issm/trunk/src/c/objects/FemModel.cpp

@@ -9 +9 @@
 #endif
 
-#include "./FemModel.h"
 #include "stdio.h"
 #include "../shared/shared.h"
 #include "../include/macros.h"
-
+#include "./FemModel.h"
+
+/*Object constructors, destructors: {{{1:*/
 FemModel::FemModel(){
 
         elements=NULL;
         nodes=NULL;
+        vertices=NULL;
         constraints=NULL;
         loads=NULL;
@@ -34 +36 @@
 }
 
-FemModel::FemModel(DataSet* femmodel_elements,DataSet* femmodel_nodes,DataSet* femmodel_constraints,DataSet* femmodel_loads,
+FemModel::FemModel(DataSet* femmodel_elements,DataSet* femmodel_nodes,DataSet* femmodel_vertices, DataSet* femmodel_constraints,DataSet* femmodel_loads,
                 DataSet* femmodel_materials,DataSet* femmodel_parameters, DofVec* femmodel_partition,DofVec* femmodel_tpartition,DofVec* femmodel_yg,
                 Mat femmodel_Rmg,Mat femmodel_Gmn,NodeSets* femmodel_nodesets,Vec femmodel_ys,Vec femmodel_ys0){
@@ -41 +43 @@
         elements=femmodel_elements;
         nodes=femmodel_nodes;
+        vertices=femmodel_vertices;
         constraints=femmodel_constraints;
         loads=femmodel_loads;
@@ -61 +64 @@
         delete elements;
         delete nodes;
+        delete vertices;
         delete loads;
         delete constraints;
@@ -76 +80 @@
 
 }
-
-
+/*}}}*/
+/* Object management: {{{1*/
 void FemModel::Echo(void){
 
@@ -83 +87 @@
         printf("   elements: %p\n",elements);
         printf("   nodes: %p\n",nodes);
+        printf("   vertices: %p\n",vertices);
         printf("   loads: %p\n",loads);
         printf("   materials: %p\n",materials);
@@ -105 +110 @@
         printf("   nodes: \n");
         nodes->Echo();
+        printf("   vertices: \n");
+        vertices->Echo();
         printf("   loads: \n");
         nodes->Echo();
@@ -204 +211 @@
 DataSet*            FemModel::get_elements(void){return elements;}
 DataSet*            FemModel::get_nodes(void){return nodes ;}
+DataSet*            FemModel::get_vertices(void){return vertices ;}
 DataSet*            FemModel::get_constraints(void){return constraints ;}
 DataSet*            FemModel::get_loads(void){return loads;}
@@ -216 +224 @@
 Vec                 FemModel::get_ys0(void){return ys0;}
 Mat                 FemModel::get_Gmn(void){return Gmn;}
+/*}}}*/

issm/trunk/src/c/objects/FemModel.h

@@ -6 +6 @@
 #define FEMMODEL_H_
 
+#include "./Object.h"
+#include "../DataSet/DataSet.h"
+#include "./DofVec.h"
 #include "../toolkits/toolkits.h"
-#include "../DataSet/DataSet.h"
 #include "../parallel/parallel.h"
-#include "./DofVec.h"
-
-class DataSet;
-struct OptArgs;
 
 class FemModel: public Object{
@@ -22 +20 @@
                 DataSet*            elements;
                 DataSet*            nodes;
+                DataSet*            vertices;
                 DataSet*            constraints;
                 DataSet*            loads;
@@ -27 +26 @@
                 DataSet*            parameters;
 
-                DofVec*                 partition;
-                DofVec*                 tpartition;
-                DofVec*                 yg;
+                DofVec*             partition;
+                DofVec*             tpartition;
+                DofVec*             yg;
+
                 Mat                 Rmg;
                 NodeSets*           nodesets;
@@ -38 +38 @@
                 FemModel();
                 ~FemModel();
-                FemModel(DataSet* elements,DataSet* nodes,DataSet* constraints,DataSet* loads,DataSet* materials,DataSet* parameters,
+                FemModel(DataSet* elements,DataSet* nodes,DataSet* vertices, DataSet* constraints,DataSet* loads,DataSet* materials,DataSet* parameters,
                                       DofVec* partition,DofVec* tpartition,DofVec* yg,Mat Rmg,Mat Gmn,NodeSets* nodesets,Vec ys,Vec ys0);
      
@@ -62 +62 @@
                 DataSet* get_elements(void);
                 DataSet* get_nodes(void);
+                DataSet* get_vertices(void){return vertices ;}
                 DataSet* get_constraints(void);
                 DataSet* get_loads(void);
@@ -75 +76 @@
                 Mat      get_Gmn(void);
 
-                
-
-
 };
 

issm/trunk/src/c/objects/Hook.cpp

@@ -12 +12 @@
 #include "stdio.h"
 #include <string.h>
+#include "./Object.h"
+#include "../DataSet/DataSet.h"
+#include "./Hook.h"
 #include "../EnumDefinitions/EnumDefinitions.h"
 #include "./ParameterInputs.h"
@@ -17 +20 @@
 #include "../include/typedefs.h"
 #include "../include/macros.h"
-#include "./Hook.h"
 
 
@@ -31 +33 @@
 /*}}}*/
 /*FUNCTION Hook::Hook(int* ids, int num){{{1*/
-Hook::Hook(int* ids, int num){
-
-        int i;
-        this->num=num;
+Hook::Hook(int* in_ids, int in_num){
+        this->Init(in_ids,in_num);
+}
+/*}}}*/
+/*FUNCTION Hook::Init(int* ids, int num){{{1*/
+Hook::Init(int* in_ids, int in_num){
+
+        int i;
+        this->num=in_num;
         
         /*Allocate: */
@@ -43 +50 @@
         /*Copy ids: */
         for (i=0;i<this->num;i++){
-                this->ids[i]=ids[i];
+                this->ids[i]=in_ids[i];
                 this->objects[i]=NULL;
                 this->offsets[i]=0;

issm/trunk/src/c/objects/Hook.h

@@ -10 +10 @@
 
 #include "./Object.h"
-#include "./../DataSet/DataSet.h"
+#include "../DataSet/DataSet.h"
 #include "../toolkits/toolkits.h"
 
 class DataSet;
-class Object;
-
 class Hook{
 
@@ -29 +27 @@
                 Hook();
                 Hook(int* ids, int num);
+                Init(int* ids, int num);
                 Hook(Object** objects, int* ids, int* offsets,int num);
                 Hook(Hook* hook);

issm/trunk/src/c/objects/Matice.cpp

@@ -17 +17 @@
                 
 /*Object constructors and destructor*/
+/*FUNCTION Matice::default constructor {{{1*/
+Matice::Matice(){
+        return;
+}
+/*}}}*/
 /*FUNCTION Matice::constructor {{{1*/
-Matice::Matice(){
-        return;
-}
-/*}}}*/
-/*FUNCTION Matice::creation {{{1*/
-Matice::Matice(int matice_mid,double matice_B,double matice_n){
-        mid=matice_mid;
-        B=matice_B;
-        n=matice_n;
-        return;
+Matice::Matice(int in_mid,double in_B,double in_n){
+        this->Init(in_mid,in_B,in_n);
+}
+/*}}}*/
+/*FUNCTION Matice::init {{{1*/
+Matice::Init(int in_mid,double in_B,double in_n){
+        this->mid=in_mid;
+        this->B=in_B;
+        this->n=in_n;
+}
+/*}}}*/
+/*FUNCTION Matice::constructor from iomodel {{{1*/
+Matice::Matice(int i, IoModel* iomodel, int num_vertices){
+
+        int j;
+        
+        /*needed for Init routine:*/
+        int matice_mid;
+        double matice_B;
+        double matice_n;
+
+        /*intermediary: */
+        double B_avg;
+        
+        /*id: */
+        matice_mid=i+1;  //same as element it is created for
+ 
+        /*B and n: */
+        B_avg=0;
+        for(j=0;j<num_vertices;j++){
+                B_avg+=*(iomodel->B+((int)*(iomodel->elements+num_vertices*i+j)-1));
+        }
+        B_avg=B_avg/num_vertices;
+        
+        matice_B=B_avg;
+        matice_n=(double)*(iomodel->n+i);
+
+        this->Init(matice_mid,matice_B,matice_n);
 }
 /*}}}*/

issm/trunk/src/c/objects/Matice.h

@@ -19 +19 @@
                 Matice();
                 Matice(int mid,double B,double n);
+                Matice(int i, IoModel* iomodel, int num_vertices);
+                Init(int mid,double B,double n);
                 ~Matice();
 

issm/trunk/src/c/objects/Matpar.cpp

@@ -18 +18 @@
                 
 /*Object constructors and destructor*/
-/*FUNCTION Matpar::constructor {{{1*/
+/*FUNCTION Matpar::default constructor {{{1*/
 Matpar::Matpar(){
         return;
 }
 /*}}}1*/
-/*FUNCTION Matpar::creation {{{1*/
+/*FUNCTION Matpar::constructorr {{{1*/
 Matpar::Matpar(int      matpar_mid, double      matpar_rho_ice, double  matpar_rho_water, double  matpar_heatcapacity, double  matpar_thermalconductivity, double  matpar_latentheat, double  matpar_beta, double  matpar_meltingpoint, double  matpar_mixed_layer_capacity, double  matpar_thermal_exchange_velocity, double  matpar_g){
 
-
-        mid=matpar_mid; 
-        rho_ice=matpar_rho_ice; 
-        rho_water=matpar_rho_water; 
-        heatcapacity=matpar_heatcapacity; 
-        thermalconductivity=matpar_thermalconductivity; 
-        latentheat=matpar_latentheat; 
-        beta=matpar_beta; 
-        meltingpoint=matpar_meltingpoint; 
-        mixed_layer_capacity=matpar_mixed_layer_capacity; 
-        thermal_exchange_velocity=matpar_thermal_exchange_velocity; 
-        g=matpar_g; 
+        this->Init(matpar_mid, matpar_rho_ice, matpar_rho_water, matpar_heatcapacity, matpar_thermalconductivity, matpar_latentheat, matpar_beta, matpar_meltingpoint, matpar_mixed_layer_capacity, matpar_thermal_exchange_velocity, matpar_g);
+
+}
+/*}}}1*/
+/*FUNCTION Matpar::constructor  from iomodel{{{1*/
+Matpar::Matpar(IoModel* iomodel){
+
+        int     matpar_mid;
+        double  matpar_rho_ice;
+        double  matpar_rho_water;
+        double  matpar_heatcapacity;
+        double  matpar_thermalconductivity;
+        double  matpar_latentheat;
+        double  matpar_beta;
+        double  matpar_meltingpoint;
+        double  matpar_mixed_layer_capacity;
+        double  matpar_thermal_exchange_velocity;
+        double  matpar_g;
+
+        matpar_mid=iomodel->numberofelements+1; //put it at the end of the materials
+        matpar_g=iomodel->g; 
+        matpar_rho_ice=iomodel->rho_ice; 
+        matpar_rho_water=iomodel->rho_water; 
+        matpar_thermalconductivity=iomodel->thermalconductivity; 
+        matpar_heatcapacity=iomodel->heatcapacity; 
+        matpar_latentheat=iomodel->latentheat; 
+        matpar_beta=iomodel->beta; 
+        matpar_meltingpoint=iomodel->meltingpoint; 
+        matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity; 
+        matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity; 
+
+        this->Init(matpar_mid, matpar_rho_ice, matpar_rho_water, matpar_heatcapacity, matpar_thermalconductivity, matpar_latentheat, matpar_beta, matpar_meltingpoint, matpar_mixed_layer_capacity, matpar_thermal_exchange_velocity, matpar_g);
+
+}
+/*}}}1*/
+/*FUNCTION Matpar::Init {{{1*/
+Matpar::Init(int        matpar_mid, double      matpar_rho_ice, double  matpar_rho_water, double  matpar_heatcapacity, double  matpar_thermalconductivity, double  matpar_latentheat, double  matpar_beta, double  matpar_meltingpoint, double  matpar_mixed_layer_capacity, double  matpar_thermal_exchange_velocity, double  matpar_g){
+
+        this->mid=matpar_mid; 
+        this->rho_ice=matpar_rho_ice; 
+        this->rho_water=matpar_rho_water; 
+        this->heatcapacity=matpar_heatcapacity; 
+        this->thermalconductivity=matpar_thermalconductivity; 
+        this->latentheat=matpar_latentheat; 
+        this->beta=matpar_beta; 
+        this->meltingpoint=matpar_meltingpoint; 
+        this->mixed_layer_capacity=matpar_mixed_layer_capacity; 
+        this->thermal_exchange_velocity=matpar_thermal_exchange_velocity; 
+        this->g=matpar_g; 
 
         return;

issm/trunk/src/c/objects/Matpar.h

@@ -28 +28 @@
         
                 Matpar(int      mid, double     rho_ice, double rho_water, double  heatcapacity, double  thermalconductivity, double  latentheat, double  beta, double  meltingpoint, double  mixed_layer_capacity, double  thermal_exchange_velocity, double  g);
+                Init(int        mid, double     rho_ice, double rho_water, double  heatcapacity, double  thermalconductivity, double  latentheat, double  beta, double  meltingpoint, double  mixed_layer_capacity, double  thermal_exchange_velocity, double  g);
+                Matpar(IoModel* iomodel);
                 
                 ~Matpar();

issm/trunk/src/c/objects/Model.cpp

@@ -43 +43 @@
         DataSet*            elements=NULL;
         DataSet*            nodes=NULL;
+        DataSet*            vertices=NULL;
         DataSet*            constraints=NULL;
         DataSet*            loads=NULL;
@@ -71 +72 @@
 
         _printf_("   create datasets:\n");
-        CreateDataSets(&elements,&nodes,&materials,&constraints,&loads,&parameters,iomodel,MODEL);
+        CreateDataSets(&elements,&nodes,&vertices,&materials,&constraints,&loads,&parameters,iomodel,MODEL);
 
         _printf_("   create degrees of freedom: \n");
-        Dofx( &partition,&tpartition,elements,nodes, parameters);
+        Dofx( &partition,&tpartition,elements,nodes, vertices,parameters);
         
         _printf_("   create single point constraints: \n");
@@ -92 +93 @@
 
         _printf_("   configuring element and loads:\n");
-        ConfigureObjectsx(elements, loads, nodes, materials,parameters);
+        ConfigureObjectsx(elements, loads, nodes, vertices, materials,parameters);
 
         _printf_("   process parameters:\n");
@@ -101 +102 @@
 
         /*Use all the data created to create an femmodel: */
-        femmodel=new FemModel(elements,nodes,constraints,loads,materials,parameters,
+        femmodel=new FemModel(elements, nodes, vertices, constraints,loads,materials,parameters,
                                       partition,tpartition,yg,Rmg,Gmn,nodesets,ys,ys0);
 
@@ -116 +117 @@
         DataSet*            elements=NULL;
         DataSet*            nodes=NULL;
+        DataSet*            vertices=NULL;
         DataSet*            constraints=NULL;
         DataSet*            loads=NULL;
@@ -148 +150 @@
 
         _printf_("   create datasets:\n");
-        CreateDataSets(&elements,&nodes,&materials,&constraints,&loads,&parameters,iomodel,MODEL);
+        CreateDataSets(&elements,&nodes,&vertices,&materials,&constraints,&loads,&parameters,iomodel,MODEL);
 
         _printf_("   create degrees of freedom: \n");
-        Dofx( &partition,&tpartition,elements,nodes, parameters);
+        Dofx( &partition,&tpartition,elements,nodes, vertices,parameters);
         
         _printf_("   create single point constraints: \n");
@@ -169 +171 @@
 
         _printf_("   configuring element and loads:\n");
-        ConfigureObjectsx(elements, loads, nodes, materials,parameters);
+        ConfigureObjectsx(elements, loads, nodes, vertices, materials,parameters);
 
         _printf_("   process parameters:\n");
@@ -178 +180 @@
 
         /*Use all the data created to create an femmodel: */
-        femmodel=new FemModel(elements,nodes,constraints,loads,materials,parameters,
+        femmodel=new FemModel(elements,nodes,vertices, constraints,loads,materials,parameters,
                                       partition,tpartition,yg,Rmg,Gmn,nodesets,ys,ys0);
 

issm/trunk/src/c/objects/Node.cpp

@@ -3 +3 @@
  */
 
-
+/*Include files: {{{1*/
 #ifdef HAVE_CONFIG_H
         #include "config.h"
@@ -11 +11 @@
 
 #include "stdio.h"
+#include "./Vertex.h"
 #include "./Node.h"
+#include "./Hook.h"
+#include "./DofIndexing.h"
 #include <string.h>
 #include "../EnumDefinitions/EnumDefinitions.h"
@@ -18 +21 @@
 #include "../include/typedefs.h"
 #include "../include/macros.h"
-
-/*Object constructors and destructor*/
-/*FUNCTION Node constructor {{{1*/
+/*}}}*/
+/*Object constructors and destructors: {{{1*/
+/*FUNCTION Node default constructor {{{2*/
 Node::Node(){
         return;
 }
 /*}}}*/
-/*FUNCTION Node constructor {{{1*/
-Node::Node(int node_id,int node_partitionborder,int node_numdofs, double node_x[3],double node_sigma,int node_onbed,int node_onsurface,int node_upper_node_id,int node_onshelf,int node_onsheet){
-
-        int i;
-        
-        id=node_id;
-        partitionborder=node_partitionborder;
-        numberofdofs=node_numdofs;
-        x[0]=node_x[0];
-        x[1]=node_x[1];
-        x[2]=node_x[2];
-        sigma=node_sigma;
-        onbed=node_onbed;
-        onsurface=node_onsurface;
-        onshelf=node_onshelf;
-        onsheet=node_onsheet;
-
-        /*Initialize sets: */
-        for(i=0;i<numberofdofs;i++){
-                mset[i]=0;
-                nset[i]=1; 
-                fset[i]=1; //we assume new nodes are not constrained in rigid body mode, or single point constraint mode.
-                sset[i]=0;
-        }
-
-        /*Initialize upper node:*/
-        upper_node_id=node_upper_node_id;
-        upper_node=NULL;
-        upper_node_offset=UNDEF;
+/*FUNCTION Node constructor {{{2*/
+Node::Node(int node_id,int node_vertex_id, int node_upper_node_id, int node_partitionborder,int node_numdofs, NodeProperties* node_properties):
+        indexing(node_numdofs,node_partitionborder),
+        properties(node_properties),
+    hvertex(&node_vertex_id,1),
+    hupper_node(&node_upper_node_id,1){
+
+        this->id=node_id;
 
         return;
 }
 /*}}}*/
-/*FUNCTION Node destructor{{{1*/
+/*FUNCTION Node other constructor {{{2*/
+Node::Node(int node_id,DofIndexing* node_indexing, NodeProperties* node_properties, Hook* node_vertex, Hook* node_upper_node):
+            indexing(node_indexing),
+                properties(node_properties),
+                hvertex(node_vertex),
+                hupper_node(node_upper_node) {
+
+            /*all the initialization has been done by the initializer, just fill in the id: */
+            this->id=node_id;
+
+}
+/*}}}*/
+/*FUNCTION Node constructor  from iomodel{{{2*/
+Node::Node(int i, IoModel* iomodel){ //i is the node index
+
+        int k;
+        
+        int numdofs;
+        int partitionborder;
+        int vertex_id;
+
+
+        /*id: */
+        this->id=i+1; //matlab indexing
+
+        /*indexing:*/
+        DistributeNumDofs(&numdofs,iomodel->analysis_type,iomodel->sub_analysis_type); //number of dofs per node
+        if(my_bordervertices[i])partitionborder=1; else partitionborder=0;//is this node on a partition border?
+
+        this->indexing.Init(numdofs,partitionborder);
+
+        /*properties: */
+        this->properties.Init(
+                        (int)iomodel->gridonbed[i],
+                        (int)iomodel->gridonsurface[i],
+                        (int)iomodel->gridoniceshelf[i],
+                        (int)iomodel->gridonicesheet[i]);
+
+        /*hooks: */
+        vertex_id=this->id; //node and vertex have the same id, as we are running galerkin continuous, with same number of nodes and vertices.
+
+        if (strcmp(iomodel->meshtype,"3d")==0){
+                if (isnan(iomodel->uppernodes[i])){
+                        upper_node_id=this->id; //nodes on surface do not have upper nodes, only themselves.
+                }
+                else{
+                        upper_node_id=(int)iomodel->uppernodes[i];
+                }
+        }
+        else{
+                /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/
+                upper_node_id=this->id;
+        }
+
+        this->hvertex->Init(vertex_id,1); //node id is the same as the vertex id, continuous galerkin!
+        this->hupper_node->Init(upper_node_id,1);
+        
+
+        /*set single point constraints.: */
+        if (strcmp(iomodel->meshtype,"3d")==0){
+                /*We have a  3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */
+                if (iomodel->deadgrids[i]){
+                        for(k=1;k<=numdofs;k++){
+                                node->FreezeDof(k);
+                        }
+                }
+        }
+        if (iomodel->gridonhutter[i]){
+                for(k=1;k<=numdofs;k++){
+                        node->FreezeDof(k);
+                }
+        }
+}
+/*}}}*/
+/*FUNCTION Node destructor{{{2*/
 Node::~Node(){
         return;
 }
 /*}}}*/
-
-/*Object management*/
-/*FUNCTION Node Marshall{{{1*/
-void  Node::Marshall(char** pmarshalled_dataset){
+/*}}}*/
+/*Object management: {{{1*/
+/*FUNCTION Node Configure {{{2*/
+void  Node::Configure(void* pnodes, void* pvertices){
+
+        int i;
+
+        DataSet* nodesin=NULL;
+        DataSet* verticesin=NULL;
+
+        /*Recover pointers :*/
+        nodesin=(DataSet*)pnodes;
+        verticesin=(DataSet*)pvertices;
+
+        /*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 ofthis->indexing.f_sets hidden in hooks: */
+        hvertex.configure(verticesin);
+        hupper_node.configure(nodesin);
+
+}
+/*}}}*/
+/*FUNCTION Node copy {{{2*/
+Object* Node::copy() {
+                
+        return new Node(this->id,&this->indexing, &this->properties, &this->hvertex,&this->hupper_node);
+
+}
+
+/*}}}*/
+/*FUNCTION Node DeepEcho{{{2*/
+void Node::DeepEcho(void){
+
+        printf("Node:\n");
+        printf("   id: %i\n",id);
+        indexing.DeepEcho();
+        properties.DeepEcho();
+        hvertex.DeepEcho();
+        hupper_node.DeepEcho();
+
+}
+/*}}}*/
+/*FUNCTION Node Demarshall{{{2*/
+void  Node::Demarshall(char** pmarshalled_dataset){
 
         char* marshalled_dataset=NULL;
-        int   enum_type=0;
 
         /*recover marshalled_dataset: */
         marshalled_dataset=*pmarshalled_dataset;
 
-        /*get enum type of Node: */
-        enum_type=NodeEnum();
-        
-        /*marshall enum: */
-        memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
-        
-        /*marshall Node data: */
-        memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
-        memcpy(marshalled_dataset,&partitionborder,sizeof(partitionborder));marshalled_dataset+=sizeof(partitionborder);
-        memcpy(marshalled_dataset,&clone,sizeof(clone));marshalled_dataset+=sizeof(clone);
-        memcpy(marshalled_dataset,&numberofdofs,sizeof(numberofdofs));marshalled_dataset+=sizeof(numberofdofs);
-        memcpy(marshalled_dataset,&x,sizeof(x));marshalled_dataset+=sizeof(x);
-        memcpy(marshalled_dataset,&sigma,sizeof(sigma));marshalled_dataset+=sizeof(sigma);
-        memcpy(marshalled_dataset,&onbed,sizeof(onbed));marshalled_dataset+=sizeof(onbed);
-        memcpy(marshalled_dataset,&onsurface,sizeof(onsurface));marshalled_dataset+=sizeof(onsurface);
-        memcpy(marshalled_dataset,&onshelf,sizeof(onshelf));marshalled_dataset+=sizeof(onshelf);
-        memcpy(marshalled_dataset,&onsheet,sizeof(onsheet));marshalled_dataset+=sizeof(onsheet);
-        memcpy(marshalled_dataset,&doflist,sizeof(doflist));marshalled_dataset+=sizeof(doflist);
-        memcpy(marshalled_dataset,&doflist1,sizeof(doflist1));marshalled_dataset+=sizeof(doflist1);
-        memcpy(marshalled_dataset,&mset,sizeof(mset));marshalled_dataset+=sizeof(mset);
-        memcpy(marshalled_dataset,&nset,sizeof(nset));marshalled_dataset+=sizeof(nset);
-        memcpy(marshalled_dataset,&fset,sizeof(fset));marshalled_dataset+=sizeof(fset);
-        memcpy(marshalled_dataset,&sset,sizeof(sset));marshalled_dataset+=sizeof(sset);
-        memcpy(marshalled_dataset,&upper_node_id,sizeof(upper_node_id));marshalled_dataset+=sizeof(upper_node_id);
-        memcpy(marshalled_dataset,&upper_node,sizeof(upper_node));marshalled_dataset+=sizeof(upper_node);
-        memcpy(marshalled_dataset,&upper_node_offset,sizeof(upper_node_offset));marshalled_dataset+=sizeof(upper_node_offset);
-
-        *pmarshalled_dataset=marshalled_dataset;
-        return;
-}
-/*}}}*/
-/*FUNCTION Node MarshallSize{{{1*/
-int   Node::MarshallSize(){
-
-        return sizeof(id)+
-                sizeof(partitionborder)+
-                sizeof(clone)+
-                sizeof(numberofdofs)+
-                sizeof(x)+
-                sizeof(sigma)+
-                sizeof(onbed)+
-                sizeof(onsurface)+
-                sizeof(onshelf)+
-                sizeof(onsheet)+
-                sizeof(doflist)+
-                sizeof(doflist1)+
-                sizeof(mset)+
-                sizeof(nset)+
-                sizeof(fset)+
-                sizeof(sset)+
-                sizeof(upper_node_id)+
-                sizeof(upper_node)+
-                sizeof(upper_node_offset)+
-                sizeof(int); //sizeof(int) for enum type
-}
-/*}}}*/
-/*FUNCTION Node Demarshall{{{1*/
-void  Node::Demarshall(char** pmarshalled_dataset){
-
-        char* marshalled_dataset=NULL;
-
-        /*recover marshalled_dataset: */
-        marshalled_dataset=*pmarshalled_dataset;
-
         /*this time, no need to get enum type, the pointer directly points to the beginning of the 
          *object data (thanks to DataSet::Demarshall):*/
 
         memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id);
-        memcpy(&partitionborder,marshalled_dataset,sizeof(partitionborder));marshalled_dataset+=sizeof(partitionborder);
-        memcpy(&clone,marshalled_dataset,sizeof(clone));marshalled_dataset+=sizeof(clone);
-        memcpy(&numberofdofs,marshalled_dataset,sizeof(numberofdofs));marshalled_dataset+=sizeof(numberofdofs);
-        memcpy(&x,marshalled_dataset,sizeof(x));marshalled_dataset+=sizeof(x);
-        memcpy(&sigma,marshalled_dataset,sizeof(sigma));marshalled_dataset+=sizeof(sigma);
-        memcpy(&onbed,marshalled_dataset,sizeof(onbed));marshalled_dataset+=sizeof(onbed);
-        memcpy(&onsurface,marshalled_dataset,sizeof(onsurface));marshalled_dataset+=sizeof(onsurface);
-        memcpy(&onshelf,marshalled_dataset,sizeof(onshelf));marshalled_dataset+=sizeof(onshelf);
-        memcpy(&onsheet,marshalled_dataset,sizeof(onsheet));marshalled_dataset+=sizeof(onsheet);
-        memcpy(&doflist,marshalled_dataset,sizeof(doflist));marshalled_dataset+=sizeof(doflist);
-        memcpy(&doflist1,marshalled_dataset,sizeof(doflist1));marshalled_dataset+=sizeof(doflist1);
-        memcpy(&mset,marshalled_dataset,sizeof(mset));marshalled_dataset+=sizeof(mset);
-        memcpy(&nset,marshalled_dataset,sizeof(nset));marshalled_dataset+=sizeof(nset);
-        memcpy(&fset,marshalled_dataset,sizeof(fset));marshalled_dataset+=sizeof(fset);
-        memcpy(&sset,marshalled_dataset,sizeof(sset));marshalled_dataset+=sizeof(sset);
-        memcpy(&upper_node_id,marshalled_dataset,sizeof(upper_node_id));marshalled_dataset+=sizeof(upper_node_id);
-        memcpy(&upper_node,marshalled_dataset,sizeof(upper_node));marshalled_dataset+=sizeof(upper_node);
-        memcpy(&upper_node_offset,marshalled_dataset,sizeof(upper_node_offset));marshalled_dataset+=sizeof(upper_node_offset);
-        
-        /*upper node is not pointing to correct object anymore: */
-        upper_node=NULL;
-
+        
+        /*demarshall objects: */
+        indexing.Demarshall(&marshalled_dataset);
+        properties.Demarshall(&marshalled_dataset);
+        hvertex.Demarshall(&marshalled_dataset);
+        hupper_node.Demarshall(&marshalled_dataset);
+        
         /*return: */
         *pmarshalled_dataset=marshalled_dataset;
@@ -169 +184 @@
 }
 /*}}}*/
-/*FUNCTION Node GetId{{{1*/
+/*FUNCTION Node Echo{{{2*/
+void Node::Echo(void){
+
+        printf("Node:\n");
+        printf("   id: %i\n",id);
+        indexing.Echo();
+        properties.Echo();
+        hvertex.Echo();
+        hupper_node.Echo();
+
+}
+/*}}}*/
+/*FUNCTION Node Enum{{{2*/
+int Node::Enum(void){
+
+        return NodeEnum();
+
+}
+/*}}}*/
+/*FUNCTION Node GetId{{{2*/
 int    Node::GetId(void){ return id; }
 /*}}}*/
-/*FUNCTION Node GetName{{{1*/
+/*FUNCTION Node GetName{{{2*/
 char* Node::GetName(void){
         return "node";
 }
 /*}}}*/
-                
-/*Object functions*/
-/*FUNCTION Node ApplyConstraints{{{1*/
+/*FUNCTION Node GetVertexDof {{{2*/
+int   Node::GetVertexDof(void){
+
+        Vertex*  vertex=NULL;
+
+        vertex=(Vertex*)hvertex.delivers();
+        return vertex->dof;
+}
+/*}}}*/
+/*FUNCTION Node Marshall{{{2*/
+void  Node::Marshall(char** pmarshalled_dataset){
+
+        char* marshalled_dataset=NULL;
+        int   enum_type=0;
+
+        /*recover marshalled_dataset: */
+        marshalled_dataset=*pmarshalled_dataset;
+
+        /*get enum type of Node: */
+        enum_type=NodeEnum();
+        
+        /*marshall enum: */
+        memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
+        
+        /*marshall Node data: */
+        memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
+        
+        /*marshall objects: */
+        indexing.Marshall(&marshalled_dataset);
+        properties.Marshall(&marshalled_dataset);
+        hvertex.Marshall(&marshalled_dataset);
+        hupper_node.Marshall(&marshalled_dataset);
+
+        *pmarshalled_dataset=marshalled_dataset;
+        return;
+}
+/*}}}*/
+/*FUNCTION Node MarshallSize{{{2*/
+int   Node::MarshallSize(){
+
+        return sizeof(id)+
+                +indexing.MarshallSize()+
+                +properties.MarshallSize()+
+                +hvertex.MarshallSize()+
+                +hupper_node.MarshallSize()+
+                sizeof(int); //sizeof(int) for enum type
+}
+/*}}}*/
+/*FUNCTION Node SetVertexDof {{{2*/
+void   Node::SetVertexDof(int in_dof){
+
+        Vertex*  vertex=NULL;
+
+        vertex=(Vertex*)hvertex.delivers();
+        vertex->dof=in_dof;
+
+}
+/*}}}*/
+/*}}}*/
+/*Object numerics: {{{1*/
+/*FUNCTION Node ApplyConstraints{{{2*/
 void  Node::ApplyConstraint(Vec yg,int dof,double value){
 
@@ -192 +284 @@
          *  we are a clone!*/
 
-        if(!clone){
-
-                index=doflist[dof-1]; //matlab indexing
+        if(!indexing.clone){
+
+                index=indexing.doflist[dof-1]; //matlab indexing
 
                 VecSetValues(yg,1,&index,&value,INSERT_VALUES);
@@ -202 +294 @@
 }
 /*}}}*/
-/*FUNCTION Node Configure {{{1*/
-void  Node::Configure(void* pnodes){
-
-        DataSet* nodes=NULL;
-
-        /*Recover pointers :*/
-        nodes=(DataSet*)pnodes;
-
-        /*Link this node with its upper node: */
-        ResolvePointers((Object**)&upper_node,&upper_node_id,&upper_node_offset,1,nodes);
-        
-}
-/*}}}*/
-/*FUNCTION Node copy {{{1*/
-Object* Node::copy() {
-        return new Node(*this); 
-}
-/*}}}*/
-/*FUNCTION Node CreatePartition{{{1*/
+/*FUNCTION Node CreatePartition{{{2*/
 void  Node::CreatePartition(Vec partition){ 
 
@@ -227 +301 @@
 
         idxm=(id-1);
-        value=(double)doflist1[0];
+        value=(double)this->GetVertexDof();
         ISSMASSERT(value>=0);
 
@@ -235 +309 @@
 }
 /*}}}*/
-/*FUNCTION Node CreateVecSets {{{1*/
+/*FUNCTION Node CreateVecSets {{{2*/
 void  Node::CreateVecSets(Vec pv_g,Vec pv_m,Vec pv_n,Vec pv_f,Vec pv_s){
 
@@ -243 +317 @@
         int i;
 
-        for(i=0;i<numberofdofs;i++){
+        for(i=0;i<this->indexing.numberofdofs;i++){
 
                 /*g set: */
-                VecSetValues(pv_g,1,&doflist[i],&gvalue,INSERT_VALUES);
+                VecSetValues(pv_g,1,&indexing.doflist[i],&gvalue,INSERT_VALUES);
                 
                 /*m set: */
-                value=(double)mset[i];
-                VecSetValues(pv_m,1,&doflist[i],&value,INSERT_VALUES);
+                value=(double)this->indexing.m_set[i];
+                VecSetValues(pv_m,1,&indexing.doflist[i],&value,INSERT_VALUES);
 
                 /*n set: */
-                value=(double)nset[i];
-                VecSetValues(pv_n,1,&doflist[i],&value,INSERT_VALUES);
+                value=(double)this->indexing.n_set[i];
+                VecSetValues(pv_n,1,&indexing.doflist[i],&value,INSERT_VALUES);
 
                 /*f set: */
-                value=(double)fset[i];
-                VecSetValues(pv_f,1,&doflist[i],&value,INSERT_VALUES);
+                value=(double)this->indexing.f_set[i];
+                VecSetValues(pv_f,1,&indexing.doflist[i],&value,INSERT_VALUES);
 
                 /*s set: */
-                value=(double)sset[i];
-                VecSetValues(pv_s,1,&doflist[i],&value,INSERT_VALUES);
-
-        }
-
-
-}
-/*}}}*/
-/*FUNCTION Node DeepEcho{{{1*/
-void Node::DeepEcho(void){
-
-        int i;
-
-        printf("Node:\n");
-        printf("   id: %i\n",id);
-        printf("   partitionborder: %i\n",partitionborder);
-        printf("   clone: %i\n",clone);
-        printf("   numberofdofs: %i\n",numberofdofs);
-        printf("   x=[%g,%g,%g]\n",x[0],x[1],x[2]);
-        printf("   sigma=%g\n",sigma);
-        printf("   onbed: %i\n",onbed);
-        printf("   onsurface: %i\n",onsurface);
-        printf("   onshelf: %i\n",onshelf);
-        printf("   onsheet: %i\n",onsheet);
-        printf("   upper_node_id=%i\n",upper_node_id);
-        printf("   upper_node_offset=%i\n",upper_node_offset);
-        printf("   doflist:|");
-        for(i=0;i<numberofdofs;i++){
-                if(i>MAXDOFSPERNODE)break;
-                printf("%i|",doflist[i]);
-        }
-        printf("   doflist1:|");
-        printf("%i|",doflist1[0]);
-        printf("\n");
-
-        printf("   set membership: m,n,f,s sets \n");
-        for(i=0;i<numberofdofs;i++){
-                if(i>MAXDOFSPERNODE)break;
-                printf("      dof %i: %i %i %i %i\n",i,mset[i],nset[i],fset[i],sset[i]);
-        }
-        printf("\n");
-        if(upper_node)printf("   upper_node pointer: %p\n",upper_node);
-
-        return;
-}               
-/*}}}*/
-/*FUNCTION Node DistributeDofs{{{1*/
+                value=(double)this->indexing.s_set[i];
+                VecSetValues(pv_s,1,&indexing.doflist[i],&value,INSERT_VALUES);
+
+        }
+
+
+}
+/*}}}*/
+/*FUNCTION Node DistributeDofs{{{2*/
 void  Node::DistributeDofs(int* pdofcount,int* pdofcount1){
 
@@ -318 +354 @@
         dofcount1=*pdofcount1;
         
-        if(clone){
+        if(indexing.clone){
                 /*This node is a clone! Don't distribute dofs, it will get them from another cpu!*/
                 return;
@@ -324 +360 @@
 
         /*This node should distribute dofs, go ahead: */
-        for(i=0;i<numberofdofs;i++){
-                doflist[i]=dofcount+i;
-        }
-        dofcount+=numberofdofs;
-
-        doflist1[0]=dofcount1;
+        for(i=0;i<this->indexing.numberofdofs;i++){
+                indexing.doflist[i]=dofcount+i;
+        }
+        dofcount+=this->indexing.numberofdofs;
+
+        SetVertexDof(dofcount1);
         dofcount1+=1;
 
@@ -339 +375 @@
 }
 /*}}}*/
-/*FUNCTION Node DofInMSet{{{1*/
+/*FUNCTION Node DofInMSet{{{2*/
 void  Node::DofInMSet(int dof){
 
         /*Put dof for this node into the m set (m set is for rigid body modes)*/
 
-        mset[dof]=1; //m and n are mutually exclusive (m for rigid body modes)
-        nset[dof]=0;
-        fset[dof]=0; //n splits into f (for which we solve) and s (single point constraints)
-        sset[dof]=0;
-}
-/*}}}*/
-/*FUNCTION Node DofInSSet {{{1*/
+        this->indexing.m_set[dof]=1; //m and n are mutually exclusive (m for rigid body modes)
+        this->indexing.n_set[dof]=0;
+        this->indexing.f_set[dof]=0; //n splits into f (for which we solve) and s (single point constraints)
+        this->indexing.s_set[dof]=0;
+}
+/*}}}*/
+/*FUNCTION Node DofInSSet {{{2*/
 void  Node::DofInSSet(int dof){
 
@@ -356 +392 @@
          * to a fixed value during computations. */
 
-        mset[dof]=0; //m and n are mutually exclusive (m for rigid body modes)
-        nset[dof]=1;
-        fset[dof]=0; //n splits into f (for which we solve) and s (single point constraints)
-        sset[dof]=1;
-}
-/*}}}*/
-/*FUNCTION Node DofIsInMSet{{{1*/
+        this->indexing.m_set[dof]=0; //m and n are mutually exclusive (m for rigid body modes)
+        this->indexing.n_set[dof]=1;
+        this->indexing.f_set[dof]=0; //n splits into f (for which we solve) and s (single point constraints)
+        this->indexing.s_set[dof]=1;
+}
+/*}}}*/
+/*FUNCTION Node DofIsInMSet{{{2*/
 int  Node::DofIsInMSet(int dof){
 
-        if (mset[dof])return 1;
+        if (this->indexing.m_set[dof])return 1;
         else return 0;
 
 }
 /*}}}*/
-/*FUNCTION Node Echo{{{1*/
-void Node::Echo(void){
-
-        int i;
-
-        printf("Node:\n");
-        printf("   id: %i\n",id);
-        printf("   partitionborder: %i\n",partitionborder);
-        printf("   clone: %i\n",clone);
-        printf("   numberofdofs: %i\n",numberofdofs);
-        printf("   x=[%g,%g,%g]\n",x[0],x[1],x[2]);
-        printf("   sigma=%g\n",sigma);
-        printf("   onbed: %i\n",onbed);
-        printf("   onsurface: %i\n",onsurface);
-        printf("   onshelf: %i\n",onshelf);
-        printf("   onsheet: %i\n",onsheet);
-        printf("   upper_node_id=%i\n",upper_node_id);
-        printf("   upper_node_offset=%i\n",upper_node_offset);
-        printf("   doflist:|");
-        for(i=0;i<numberofdofs;i++){
-                if(i>MAXDOFSPERNODE)break;
-                printf("%i|",doflist[i]);
-        }
-        printf("   doflist1:|");
-        printf("%i|",doflist1[0]);
-        printf("\n");
-
-        printf("   set membership: m,n,f,s sets \n");
-        for(i=0;i<numberofdofs;i++){
-                if(i>MAXDOFSPERNODE)break;
-                printf("      dof %i: %i %i %i %i\n",i,mset[i],nset[i],fset[i],sset[i]);
-        }
-        printf("\n");
-        if(upper_node)printf("   upper_node pointer: %p\n",upper_node);
-
-        return;
-}
-/*}}}*/
-/*FUNCTION Node Enum{{{1*/
-int Node::Enum(void){
-
-        return NodeEnum();
-
-}
-/*}}}*/
-/*FUNCTION Node FieldAverageOntoVertices{{{1*/
+/*FUNCTION Node FieldAverageOntoVertices{{{2*/
 void  Node::FieldAverageOntoVertices(Vec field,double* field_serial,char* fieldname){
 
@@ -421 +412 @@
 }
 /*}}}*/
-/*FUNCTION Node FieldDepthAverageAtBase{{{1*/
+/*FUNCTION Node FieldDepthAverageAtBase{{{2*/
 void  Node::FieldDepthAverageAtBase(Vec field,double* field_serial,char* fieldname){
 
         /* node data: */
-        int          doflist1;
+        int          vertexdof;
         int          dofx,dofy;
         int          isnodeonsurface;
@@ -436 +427 @@
         /*Are we on the base, not on the surface, and not on a clone node?:*/
         
-        if(onbed==1 & clone==0 &onsurface==0){
+        if(properties.onbed==1 & indexing.clone==0 &properties.onsurface==0){
                         
-                doflist1=GetDofList1();
+                vertexdof=this->GetVertexDof();
 
                 /*this node is on the bed. We are going to, follow the upper nodes until we reach the surface. At each upper node, 
@@ -456 +447 @@
 
                         /*get dofs for this base node velocity: we know there are two dofs in field_serial */
-                        dofx=2*doflist1;
-                        dofy=2*doflist1+1;
+                        dofx=2*vertexdof;
+                        dofy=2*vertexdof+1;
 
                         node=this;
@@ -464 +455 @@
                                 if (node->IsOnSurface())break;
 
-                                doflist1=node->GetDofList1();
+                                vertexdof=node->GetVertexDof();
                                 
-                                velocity1[0]=field_serial[2*doflist1];
-                                velocity1[1]=field_serial[2*doflist1+1];
+                                velocity1[0]=field_serial[2*vertexdof];
+                                velocity1[1]=field_serial[2*vertexdof+1];
                                 z1=node->GetZ();
 
                                 upper_node=node->GetUpperNode();
-                                doflist1=upper_node->GetDofList1();
+                                vertexdof=upper_node->GetVertexDof();
                         
-                                velocity2[0]=field_serial[2*doflist1];
-                                velocity2[1]=field_serial[2*doflist1+1];
+                                velocity2[0]=field_serial[2*vertexdof];
+                                velocity2[1]=field_serial[2*vertexdof+1];
                                 z2=upper_node->GetZ();
 
@@ -507 +498 @@
 
                         /*get dofs for this base node velocity: we know there are two dofs in field_serial */
-                        dofx=doflist1;
+                        dofx=vertexdof;
 
                         node=this;
@@ -514 +505 @@
                                 if (node->IsOnSurface())break;
 
-                                doflist1=node->GetDofList1();
+                                vertexdof=node->GetVertexDof();
                                 
-                                field1=field_serial[doflist1];
+                                field1=field_serial[vertexdof];
                                 z1=node->GetZ();
 
                                 upper_node=node->GetUpperNode();
-                                doflist1=upper_node->GetDofList1();
+                                vertexdof=upper_node->GetVertexDof();
                         
-                                field2=field_serial[doflist1];
+                                field2=field_serial[vertexdof];
                                 z2=upper_node->GetZ();
 
@@ -543 +534 @@
 }
 /*}}}*/
-/*FUNCTION Node FieldExtrude {{{1*/
+/*FUNCTION Node FieldExtrude {{{2*/
 void  Node::FieldExtrude(Vec field,double* field_serial,char* field_name){
                 
@@ -552 +543 @@
 
         /*Is this node on bed? :*/
-        if (onbed){
+        if (properties.onbed){
 
                 if (strcmp(field_name,"velocity")==0){
 
                         /* node data: */
-                        int          doflist1;
+                        int          vertexdof;
                         int          doflist[2];
                         double       fieldnode[2];
 
-                        doflist1=GetDofList1();
+                        vertexdof=GetVertexDof();
 
                         /*get dofs for this base node velocity: we know there are two dofs in field_serial */
-                        doflist[0]=2*doflist1;
-                        doflist[1]=2*doflist1+1;
+                        doflist[0]=2*vertexdof;
+                        doflist[1]=2*vertexdof+1;
 
                         //initilaize node
@@ -578 +569 @@
                         for(;;){
 
-                                doflist1=node->GetDofList1();
-                                doflist[0]=2*doflist1;
-                                doflist[1]=2*doflist1+1;
+                                vertexdof=node->GetVertexDof();
+                                doflist[0]=2*vertexdof;
+                                doflist[1]=2*vertexdof+1;
                                 VecSetValues(field,1,&doflist[0],&fieldnode[0],INSERT_VALUES);
                                 VecSetValues(field,1,&doflist[1],&fieldnode[1],INSERT_VALUES);
@@ -597 +588 @@
                         //initilaize node and get dof1
                         node=this;
-                        dof1=node->GetDofList1();
+                        dof1=node->GetVertexDof();
 
                         /*get field for this base node: */
@@ -606 +597 @@
                         for(;;){
 
-                                dof1=node->GetDofList1();
+                                dof1=node->GetVertexDof();
                                 VecSetValues(field,1,&dof1,&fieldel,INSERT_VALUES);
 
@@ -655 +646 @@
 }
 /*}}}*/
-/*FUNCTION Node FreezeDof{{{1*/
+/*FUNCTION Node FreezeDof{{{2*/
 void  Node::FreezeDof(int dof){
         
@@ -662 +653 @@
 }
 /*}}}*/
-/*FUNCTION Node GetDof {{{1*/
+/*FUNCTION Node GetDof {{{2*/
 int   Node::GetDof(int dofindex){
 
-        return doflist[dofindex];
-
-}
-/*}}}*/
-/*FUNCTION Node GetDofList{{{1*/
+        return indexing.doflist[dofindex];
+
+}
+/*}}}*/
+/*FUNCTION Node GetDofList{{{2*/
 void  Node::GetDofList(int* outdoflist,int* pnumberofdofspernode){
 
         int i;
-        for(i=0;i<numberofdofs;i++){
-                outdoflist[i]=doflist[i];
+        for(i=0;i<this->indexing.numberofdofs;i++){
+                outdoflist[i]=indexing.doflist[i];
         }
         /*Assign output pointers:*/
-        *pnumberofdofspernode=numberofdofs;
-}
-/*}}}*/
-/*FUNCTION Node GetDOfList1{{{1*/
-int   Node::GetDofList1(void){
-        return doflist1[0];
-}
-/*}}}*/
-/*FUNCTION Node GetNumberOfDofs{{{1*/
+        *pnumberofdofspernode=this->indexing.numberofdofs;
+}
+/*}}}*/
+/*FUNCTION Node GetNumberOfDofs{{{2*/
 int   Node::GetNumberOfDofs(){
         
-        return numberofdofs;
-
-}
-/*}}}*/
-/*FUNCTION Node GetSigma {{{1*/
-double Node::GetSigma(){return sigma;}
-/*}}}*/
-/*FUNCTION Node GetUpperNode {{{1*/
+        return this->indexing.numberofdofs;
+
+}
+/*}}}*/
+/*FUNCTION Node GetSigma {{{2*/
+double Node::GetSigma(){
+        Vertex* vertex=NULL;
+
+        vertex=(Vertex*)hvertex.delivers();
+        return vertex->sigma;
+}
+/*}}}*/
+/*FUNCTION Node GetUpperNode {{{2*/
 Node* Node::GetUpperNode(){
+        Node* upper_node=NULL;
+        upper_node=(Node*)hupper_node.delivers();
         return upper_node;
 }
 /*}}}*/
-/*FUNCTION Node GetX {{{1*/
-double Node::GetX(){return x[0];}
-/*}}}*/
-/*FUNCTION Node GetY {{{1*/
-double Node::GetY(){return x[1];}
-/*}}}*/
-/*FUNCTION Node GetZ {{{1*/
-double Node::GetZ(){return x[2];}
-/*}}}*/
-/*FUNCTION Node IsClone {{{1*/
+/*FUNCTION Node GetX {{{2*/
+double Node::GetX(){
+        Vertex* vertex=NULL;
+
+        vertex=(Vertex*)hvertex.delivers();
+        return vertex->x;
+}
+/*}}}*/
+/*FUNCTION Node GetY {{{2*/
+double Node::GetY(){
+        Vertex* vertex=NULL;
+
+        vertex=(Vertex*)hvertex.delivers();
+        return vertex->y;
+}
+/*}}}*/
+/*FUNCTION Node GetZ {{{2*/
+double Node::GetZ(){
+        Vertex* vertex=NULL;
+
+        vertex=(Vertex*)hvertex.delivers();
+        return vertex->z;
+}
+/*}}}*/
+/*FUNCTION Node IsClone {{{2*/
 int   Node::IsClone(){
         
-        return clone;
-
-}
-/*}}}*/
-/*FUNCTION Node IsOnBed {{{1*/
+        return indexing.clone;
+
+}
+/*}}}*/
+/*FUNCTION Node IsOnBed {{{2*/
 int   Node::IsOnBed(){
-        return onbed;
-}
-/*}}}*/
-/*FUNCTION Node IsOnSheet {{{1*/
+        return properties.onbed;
+}
+/*}}}*/
+/*FUNCTION Node IsOnSheet {{{2*/
 int   Node::IsOnSheet(){
-        return onsheet;
+        return properties.onsheet;
 }               
 /*}}}*/
-/*FUNCTION Node IsOnShelf {{{1*/
+/*FUNCTION Node IsOnShelf {{{2*/
 int   Node::IsOnShelf(){
-        return onshelf;
-}
-/*}}}*/
-/*FUNCTION Node IsOnSurface {{{1*/
+        return properties.onshelf;
+}
+/*}}}*/
+/*FUNCTION Node IsOnSurface {{{2*/
 int   Node::IsOnSurface(){
-        return onsurface;
-}
-/*}}}*/
-/*FUNCTION Node MyRank{{{1*/
+        return properties.onsurface;
+}
+/*}}}*/
+/*FUNCTION Node MyRank{{{2*/
 int    Node::MyRank(void){ 
         extern int my_rank;
@@ -743 +751 @@
 }
 /*}}}*/
-/*FUNCTION Node SetClone {{{1*/
+/*FUNCTION Node SetClone {{{2*/
 void  Node::SetClone(int* minranks){
 
@@ -749 +757 @@
 
         if (minranks[id-1]==my_rank){
-                clone=0;
+                indexing.clone=0;
         }
         else{
                 /*!there is a cpu with lower rank that has the same node, 
                 therefore, I am a clone*/
-                clone=1;        
-        }
-
-}
-/*}}}*/
-/*FUNCTION Node ShowBorderDofs{{{1*/
+                indexing.clone=1;       
+        }
+
+}
+/*}}}*/
+/*FUNCTION Node ShowBorderDofs{{{2*/
 void  Node::ShowBorderDofs(int* borderdofs,int* borderdofs1){
 
@@ -766 +774 @@
         
         /*Is this node on the partition border? */
-        if(!partitionborder)return;
+        if(!indexing.partitionborder)return;
 
         /*Are we the cpu that created this node's dof list? */
-        if(clone)return;
+        if(indexing.clone)return;
 
         /*Ok, we are on the partition border, and we did create the 
          * dofs for this node, plug the doflist into borderdofs: */
-        for(j=0;j<numberofdofs;j++){
-                *(borderdofs+numberofdofs*(id-1)+j)=doflist[j];
-        }
-        *(borderdofs1+(id-1)+0)=doflist1[0];
+        for(j=0;j<this->indexing.numberofdofs;j++){
+                *(borderdofs+this->indexing.numberofdofs*(id-1)+j)=indexing.doflist[j];
+        }
+        *(borderdofs1+(id-1)+0)=this->GetVertexDof();
 
         return;
 }
 /*}}}*/
-/*FUNCTION Node UpdateBorderDofs{{{1*/
+/*FUNCTION Node UpdateBorderDofs{{{2*/
 void  Node::UpdateBorderDofs(int* allborderdofs,int* allborderdofs1){ 
 
@@ -788 +796 @@
         
         /*Is this node on the partition border? */
-        if(!partitionborder)return;
+        if(!indexing.partitionborder)return;
         
         /*Are we the cpu that created this node's dof list? */
-        if(clone==0)return;
+        if(indexing.clone==0)return;
 
         /*Ok, we are on the partition border, but we did not create 
          * the dofs for this node. Therefore, our doflist is garbage right 
          * now. Go pick it up in the allborderdofs: */
-        for(j=0;j<numberofdofs;j++){
-                doflist[j]=*(allborderdofs+numberofdofs*(id-1)+j);
-        }
-        doflist1[0]=*(allborderdofs1+(id-1)+0);
+        for(j=0;j<this->indexing.numberofdofs;j++){
+                indexing.doflist[j]=*(allborderdofs+this->indexing.numberofdofs*(id-1)+j);
+        }
+        this->SetVertexDof(*(allborderdofs1+(id-1)+0));
         return; 
 }
 /*}}}*/
-/*FUNCTION Node UpdateDofs{{{1*/
+/*FUNCTION Node UpdateDofs{{{2*/
 void  Node::UpdateDofs(int dofcount,int dofcount1){
         
@@ -809 +817 @@
         extern int my_rank;
         
-        if(clone){
+        if(indexing.clone){
                 /*This node is a clone, don't update the dofs!: */
                 return;
@@ -815 +823 @@
 
         /*This node should update the dofs, go ahead: */
-        for(i=0;i<numberofdofs;i++){
-                doflist[i]+=dofcount;
-        }
-        doflist1[0]+=dofcount1;
+        for(i=0;i<this->indexing.numberofdofs;i++){
+                indexing.doflist[i]+=dofcount;
+        }
+        this->SetVertexDof(this->GetVertexDof()+dofcount1);
 
         return; 
 }
 /*}}}*/
-/*FUNCTION Node UpdateFromInputs {{{1*/
+/*FUNCTION Node UpdateFromInputs {{{2*/
 void  Node::UpdateFromInputs(void* vinputs){
-        
-        ParameterInputs* inputs=NULL;
-        Node* node=this;
-        int dof[1]={0};
-
-        /*Recover parameter inputs: */
-        inputs=(ParameterInputs*)vinputs;
-
-        /*Update internal data if inputs holds new values: */
-        inputs->Recover("x",&x[0],1,dof,1,(void**)&node);
-        inputs->Recover("y",&x[1],1,dof,1,(void**)&node);
-        inputs->Recover("z",&x[2],1,dof,1,(void**)&node);
-
-        
-}
-/*}}}*/
-/*FUNCTION Node UpdateNodePosition {{{1*/
-void  Node::UpdateNodePosition(double* thickness,double* bed){
-
-        int dof;
-
-        dof=this->GetDofList1();
-
-        /*sigma remains constant. z=bed+sigma*thickness*/
-        this->x[2]=bed[dof]+sigma*thickness[dof];
-
-}
-/*}}}*/
+
+        ISSMERROR("not used yet!");
+        
+}
+/*}}}*/
+/*}}}*/

issm/trunk/src/c/objects/Node.h

@@ -6 +6 @@
 #define _NODE_H_
 
+/*indefinitions: */
+class Object;
+class Vertex;
+class Hook;
+
 #include "./Object.h"
+#include "../DataSet/DataSet.h"
+#include "./Hook.h"
+#include "./Vertex.h"
+#include "./DofIndexing.h"
+#include "./NodeProperties.h"
 #include "../toolkits/toolkits.h"
-
-#define MAXDOFSPERNODE 4
-#define PSSTRINGLENGTH 20
 
 class Node: public Object{
 
         private: 
-                /*raw data coming from the model processor: */
-                int         id; /*! id, to track it*/
-                int     partitionborder; /*! during parallel partitioning, does this grid belong to a partition border?*/
-                int     numberofdofs; //real number of dofs, between 0 and MAXDOFSPERNODE
-                int     clone;  /*!this nodes is one the partition border, and is cloned*/
-                double  x[3]; /*! coordinates*/
-                double  sigma; /*! sigma = (z-bed)/thickness*/
+
+                int         id; 
+                                
+                DofIndexing    indexing;
+                NodeProperties properties;
+                Hook           hvertex;
+                Hook           hupper_node;
+
                 
-                int         onbed; /*! for 3d, on bedrock*/
-                int         onsurface; /*! for 3d, on surface*/
-                int         onshelf; 
-                int         onsheet;
-
-                /*for dof constraining: */
-                int     mset[MAXDOFSPERNODE];
-                int     nset[MAXDOFSPERNODE];
-                int     fset[MAXDOFSPERNODE];
-                int     sset[MAXDOFSPERNODE];
-
-                int   upper_node_id;
-                Node* upper_node;
-                int   upper_node_offset;
-
-                /*data that is post processed : */
-                int doflist[MAXDOFSPERNODE]; //dof list on which we solve
-                int doflist1[1]; //1d dof list to recover input parameters
-
         public:
 
+                /*FUNCTION constructors, destructors {{{1*/
                 Node();
-                Node(int node_id,int node_partitionborder,int node_numdofs, double node_x[3],double sigma,int node_onbed,int node_onsurface,int upper_node_id,int onshelf,int onsheet);
+                Node(int id,int vertex_id, int upper_node_id, int partitionborder,int numberofdofs, NodeProperties* node_properties);
+                Node(int id,DofIndexing* indexing, NodeProperties* properties, Hook* vertex, Hook* upper_node);
+                Node(int i, IoModel* iomodel);
                 ~Node();
-
+                /*}}}*/
+                /*FUNCTION object management {{{1*/
+                void  Configure(void* pnodes, void* pvertices);
+                void  DeepEcho();
+                void  Demarshall(char** pmarshalled_dataset);
                 void  Echo();
-                void  DeepEcho();
+                int   Enum();
+                int   GetId(void); 
+                char* GetName();
+                int   GetVertexDof(void);
                 void  Marshall(char** pmarshalled_dataset);
                 int   MarshallSize();
-                char* GetName();
-                void  Demarshall(char** pmarshalled_dataset);
-                int   Enum();
-                int   GetId(void); 
                 int   MyRank(void);
+                void  SetVertexDof(int in_dof);
+                /*}}}*/
+                /*FUNCTION numerical routines {{{1*/
                 void  DistributeDofs(int* pdofcount,int* pdofcount1);
                 void  UpdateDofs(int dofcount,int dofcount1);
@@ -79 +77 @@
                 Object* copy();
                 void  UpdateFromInputs(void* inputs);
-                void  Configure(void* pnodes);
                 Node* GetUpperNode();
                 int   IsOnBed();
@@ -90 +87 @@
                 void  FieldExtrude(Vec field,double* field_serial,char* field_name);
                 void  UpdateNodePosition(double* thickness,double* bed);
+                /*}}}*/
 };
 
 #endif  /* _NODE_H_ */
-

issm/trunk/src/c/objects/Penta.cpp

@@ -51 +51 @@
 
         return;
+}
+/*}}}*/
+/*FUNCTION Penta other constructor {{{1*/
+Penta::Penta(int i, IoModel* iomodel){ //i is the element index
+
+        int j;
+        int offset;
+        int penta_node_ids[6];
+        double penta_h[6];
+        double penta_s[6];
+        double penta_b[6];
+        double penta_k[6];
+        double penta_melting[6];
+        double penta_accumulation[6];
+
+        /*id: */
+        this->id=i+1;
+
+        /*hooks: */
+        for(j=0;j<6;j++){ //go recover node ids, needed to initialize the node hook.
+                penta_node_ids[j]=(int)*(iomodel->elements+6*i+j); //ids for vertices are in the elements array from Matlab
+        }
+        penta_matice_id=i+1; //refers to the corresponding ice material object
+        penta_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object
+        penta_numpar_id=1; //refers to numerical parameters object
+
+        this->hnodes->Init(penta_node_ids,6);
+        this->hmatice->Init(penta_matice_id,1);
+        this->hmatpar->Init(penta_matpar_id,1);
+        this->hnumpar->Init(penta_numpar_id,1);
+
+        /*properties: */
+        for(j=0;j<6;j++){
+
+                offset=(int)*(iomodel->elements+6*i+j)-1; //get index of this vertex, in C numbering. 
+                penta_h[j]= *(iomodel->thickness+offset); 
+                penta_s[j]= *(iomodel->surface+offset); 
+                penta_b[j]=*(iomodel->bed+offset);
+                penta_k[j]=*(iomodel->drag+offset);
+                penta_melting[j]=*(iomodel->melting+offset);
+                penta_accumulation[j]=*(iomodel->accumulation+offset);
+        }
+        penta_friction_type=(int)iomodel->drag_type;
+        penta_p=iomodel->p[i];
+        penta_q=iomodel->q[i];
+        penta_shelf=(int)*(iomodel->elementoniceshelf+i);
+        penta_onwater=(bool)*(iomodel->elementonwater+i);
+        penta_onsurface=(int)*(iomodel->elementonsurface+i);
+        penta_onwater=(bool)*(iomodel->elementonwater+i);
+
+        if (*(iomodel->elements_type+2*i+0)==macayealformulationenum()){ //elements of type 3 are macayeal type penta. we collapse the formulation on their base.
+                penta_collapse=1;
+        }
+        else{
+                penta_collapse=0;
+        }
+
+        this->properties->Init(6,penta_h, penta_s, penta_b, penta_k, penta_melting, penta_accumulation, NULL, penta_friction_type, penta_p, penta_q, penta_shelf, penta_onbed, penta_onwater, penta_onsurface, penta_collapse, UNDEF);
+
+
 }
 /*}}}*/

issm/trunk/src/c/objects/Penta.h

@@ -11 +11 @@
 #include "./Numpar.h"
 #include "./Matice.h"
-#include "./Node.h"
 #include "./Tria.h"
 #include "./Hook.h"
 #include "./ElementProperties.h"
 #include "./ParameterInputs.h"
+#include "./Node.h"
 
-class Numpar;
+class Object;
 class Node;
-class Matice;
-class Matpar;
-class  ElementProperties;
+class Hook;
+class ElementProperties;
+class DataSet;
 
 class Penta: public Element{
@@ -41 +41 @@
                 Penta(int penta_id,int* penta_node_ids, int penta_matice_id, int penta_matpar_id, int penta_numpar_id, ElementProperties* properties);
                 Penta(int penta_id,Hook* penta_hnodes, Hook* penta_hmatice, Hook* penta_hmatpar, Hook* penta_hnumpar, ElementProperties* penta_properties);
+                Penta(int i, IoModel* iomodel);
                 ~Penta();
                 /*}}}*/

issm/trunk/src/c/objects/Tria.cpp

@@ -59 +59 @@
 
         return;
+}
+/*}}}*/
+/*FUNCTION Tria constructor  from iomodel{{{1*/
+Tria::Tria(int i, IoModel* iomodel){ //i is the element index
+
+        int j;
+        int offset;
+        int tria_node_ids[3];
+        double tria_h[3];
+        double tria_s[3];
+        double tria_b[3];
+        double tria_k[3];
+        double tria_melting[3];
+        double tria_accumulation[3];
+
+        /*id: */
+        this->id=i+1;
+
+        /*hooks: */
+        for(j=0;j<3;j++){ //go recover node ids, needed to initialize the node hook.
+                tria_node_ids[j]=(int)*(iomodel->elements+3*i+j); //ids for vertices are in the elements array from Matlab
+        }
+        tria_matice_id=i+1; //refers to the corresponding ice material object
+        tria_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object
+        tria_numpar_id=1; //refers to numerical parameters object
+
+        this->hnodes->Init(tria_node_ids,3);
+        this->hmatice->Init(tria_matice_id,1);
+        this->hmatpar->Init(tria_matpar_id,1);
+        this->hnumpar->Init(tria_numpar_id,1);
+
+        /*properties: */
+        for(j=0;j<3;j++){
+
+                offset=(int)*(iomodel->elements+3*i+j)-1; //get index of this vertex, in C numbering. 
+                tria_h[j]= *(iomodel->thickness+offset); 
+                tria_s[j]= *(iomodel->surface+offset); 
+                tria_b[j]=*(iomodel->bed+offset);
+                tria_k[j]=*(iomodel->drag+offset);
+                tria_melting[j]=*(iomodel->melting+offset);
+                tria_accumulation[j]=*(iomodel->accumulation+offset);
+        }
+        tria_friction_type=(int)iomodel->drag_type;
+        tria_p=iomodel->p[i];
+        tria_q=iomodel->q[i];
+        tria_shelf=(int)*(iomodel->elementoniceshelf+i);
+        tria_onwater=(bool)*(iomodel->elementonwater+i);
+        
+        this->properties->Init(3,tria_h, tria_s, tria_b, tria_k, tria_melting, tria_accumulation, NULL,
+                        tria_friction_type, tria_p, tria_q, tria_shelf, UNDEF,tria_onwater, UNDEF,UNDEF,UNDEF);
+
+
 }
 /*}}}*/

issm/trunk/src/c/objects/Tria.h

@@ -7 +7 @@
 
 
+#include "./Object.h"
 #include "./Element.h"
+#include "./Hook.h"
+#include "./Node.h"
+#include "./ElementProperties.h"
 #include "../DataSet/DataSet.h"
-#include "./Object.h"
-#include "./Node.h"
-#include "./Matice.h"
-#include "./Matpar.h"
-#include "./Numpar.h"
-#include "./ParameterInputs.h"
-#include "./ElementProperties.h"
-#include "./Hook.h"
 
 class Object;
+class Element;
+class Hook;
+class ElementProperties;
+class DataSet;
+class Node;
 
 class Tria: public Element{
@@ -39 +40 @@
                 Tria(int tria_id,int* tria_node_ids, int tria_matice_id, int tria_matpar_id, int tria_numpar_id, ElementProperties* tria_properties);
                 Tria(int tria_id,Hook* tria_hnodes, Hook* tria_hmatice, Hook* tria_hmatpar, Hook* tria_hnumpar, ElementProperties* tria_properties);
+                Tria(int i, IoModel* iomodel);
                 ~Tria();
                 /*}}}*/

issm/trunk/src/c/objects/Vertex.cpp

@@ -3 +3 @@
  */
 
+/*Include files: {{{1*/
 #ifdef HAVE_CONFIG_H
         #include "config.h"
@@ -16 +17 @@
 #include "../include/typedefs.h"
 #include "../include/macros.h"
+/*}}}*/
 
 /*Object constructors and destructor*/
@@ -24 +26 @@
 /*}}}*/
 /*FUNCTION Vertex constructor {{{1*/
-Vertex::Vertex(int tria_id, double tria_x, double tria_y, double tria_z){
+Vertex::Vertex(int tria_id, double tria_x, double tria_y, double tria_z, double tria_sigma, int partitionborder){
+        this->Init(tria_id, tria_x, tria_y, tria_z, tria_sigma, partitionborder);
+}
+/*}}}*/
+/*FUNCTION Vertex init, used by constructor {{{1*/
+Vertex::Init(int tria_id, double tria_x, double tria_y, double tria_z, double tria_sigma, int partitionborder){
 
         /*all the initialization has been done by the initializer, just fill in the id: */
@@ -31 +38 @@
         this->y=tria_y;
         this->z=tria_z;
-
-        return;
+        this->sigma=tria_sigma;
+        this->dof=UNDEF;
+
+        return;
+}
+/*}}}*/
+/*FUNCTION Vertex constructor  from iomodel{{{1*/
+Vertex::Vertex(int i, IoModel* iomodel){
+
+        int partitionborder;
+
+        /*is this vertex on a partition border? */
+        if(my_bordervertices[i])partitionborder=1;
+        else partitionborder=0;
+
+        this->Init(i+1, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]),partitionborder);
+
 }
 /*}}}*/
@@ -71 +93 @@
         memcpy(&y,marshalled_dataset,sizeof(y));marshalled_dataset+=sizeof(y);
         memcpy(&z,marshalled_dataset,sizeof(z));marshalled_dataset+=sizeof(z);
+        memcpy(&sigma,marshalled_dataset,sizeof(sigma));marshalled_dataset+=sizeof(sigma);
+        memcpy(&dof,marshalled_dataset,sizeof(dof));marshalled_dataset+=sizeof(dof);
+        memcpy(&partitionborder,marshalled_dataset,sizeof(partitionborder));marshalled_dataset+=sizeof(partitionborder);
+        memcpy(&clone,marshalled_dataset,sizeof(clone));marshalled_dataset+=sizeof(clone);
 
         /*return: */
@@ -86 +112 @@
         printf("   y: %g\n",y);
         printf("   z: %g\n",z);
+        printf("   sigma: %g\n",sigma);
+        printf("   dof: %g\n",dof);
+        printf("   partitionborder: %g\n",partitionborder);
+        printf("   clone: %g\n",clone);
 
         return;
@@ -125 +155 @@
         memcpy(marshalled_dataset,&y,sizeof(y));marshalled_dataset+=sizeof(y);
         memcpy(marshalled_dataset,&z,sizeof(z));marshalled_dataset+=sizeof(z);
+        memcpy(marshalled_dataset,&sigma,sizeof(sigma));marshalled_dataset+=sizeof(sigma);
+        memcpy(marshalled_dataset,&dof,sizeof(dof));marshalled_dataset+=sizeof(dof);
+        memcpy(marshalled_dataset,&partitionborder,sizeof(partitionborder));marshalled_dataset+=sizeof(partitionborder);
+        memcpy(marshalled_dataset,&clone,sizeof(clone));marshalled_dataset+=sizeof(clone);
 
         *pmarshalled_dataset=marshalled_dataset;
@@ -137 +171 @@
                 sizeof(y)+
                 sizeof(z)+
+                sizeof(sigma)+
+                sizeof(dof)+
+                sizeof(partitionborder)+
+                sizeof(clone)+
                 +sizeof(int); //sizeof(int) for enum type
 }
@@ -155 +193 @@
 /*FUNCTION UpdateFromInputs {{{1*/
 void  Vertex::UpdateFromInputs(void* vinputs){
-
+        
+        ParameterInputs* inputs=NULL;
+        Vertex* vertex=NULL;
+        int dof[1]={0};
+        double coord[3];
+
+        vertex=this;
+
+        coord[0]=this->x; coord[1]=this->y; coord[2]=this->z;
+
+        /*Recover parameter inputs: */
+        inputs=(ParameterInputs*)vinputs;
+
+        /*Update internal data if inputs holds new values: */
+        inputs->Recover("x",&coord[0],1,dof,1,(void**)&vertex);
+        inputs->Recover("y",&coord[1],1,dof,1,(void**)&vertex);
+        inputs->Recover("z",&coord[2],1,dof,1,(void**)&vertex);
+        
         ISSMERROR("not supported yet!");
         
 }
 /*}}}*/
+/*FUNCTION UpdateVertexPosition {{{1*/
+void  Vertex::UpdatePosition(double* thickness,double* bed){
+
+        /*sigma remains constant. z=bed+sigma*thickness*/
+        this->z=bed[dof]+sigma*thickness[this->dof];
+
+}
+/*}}}*/

issm/trunk/src/c/objects/Vertex.h

@@ -10 +10 @@
 class Vertex: public Object{
 
-        private: 
+        public: 
 
                 int  id;
@@ -16 +16 @@
                 double y;
                 double z;
+                double sigma; //sigma coordinate: (z-bed)/thickness
 
-        public:
+                /*dof management: */
+                int    partitionborder;
+                int    clone;
+                int    dof; //dof to recover values in a vertex indexed vector
 
                 /*FUNCTION constructors, destructors {{{1*/
                 Vertex();
-                Vertex(int id, double x, double y, double z);
+                Vertex(int id, double x, double y, double z, double sigma,int partitionborder); 
+                Init(int id, double x, double y, double z, double sigma,int partitionborder);
+                Vertex(int i, IoModel* iomodel);
                 ~Vertex();
                 /*}}}*/
@@ -37 +43 @@
                 void  UpdateFromDakota(void* vinputs);
                 void  UpdateFromInputs(void* vinputs);
+                void  UpdatePosition(double* thickness,double* bed);
+
+
                 /*}}}*/
 

issm/trunk/src/m/enum/BeamEnum.m

@@ -7 +7 @@
 %      macro=BeamEnum()
 
-macro=415;
+macro=416;

issm/trunk/src/m/enum/NodeEnum.m

@@ -7 +7 @@
 %      macro=NodeEnum()
 
-macro=420;
+macro=421;

issm/trunk/src/m/enum/PentaEnum.m

@@ -7 +7 @@
 %      macro=PentaEnum()
 
-macro=413;
+macro=414;

issm/trunk/src/m/enum/SingEnum.m

@@ -7 +7 @@
 %      macro=SingEnum()
 
-macro=414;
+macro=415;

issm/trunk/src/m/solutions/jpl/CreateFemModel.m

@@ -9 +9 @@
         
         displaystring(md.verbose,'\n   reading data from model %s...',md.name);
-        [m.elements,m.nodes,m.constraints,m.loads,m.materials,m.parameters]=ModelProcessor(md);
+        [m.elements,m.nodes,m.vertices,m.constraints,m.loads,m.materials,m.parameters]=ModelProcessor(md);
 
         displaystring(md.verbose,'%s','   generating degrees of freedom...');
-        [m.nodes,m.part,m.tpart]=Dof(m.elements,m.nodes,m.parameters);
+        [m.nodes,m.vertices,m.part,m.tpart]=Dof(m.elements,m.nodes,m.vertices, m.parameters);
 
         displaystring(md.verbose,'%s','   generating single point constraints...');
@@ -30 +30 @@
 
         displaystring(md.verbose,'%s','   configuring element and loads...');
-        [m.elements,m.loads,m.nodes,m.parameters] = ConfigureObjects( m.elements, m.loads, m.nodes, m.materials,m.parameters);
+        [m.elements,m.loads,m.nodes,m.vertices,m.parameters] = ConfigureObjects( m.elements, m.loads, m.nodes, m.materials,m.parameters);
 
         displaystring(md.verbose,'%s','   processing parameters...');

issm/trunk/src/mex/Dof/Dof.cpp

@@ -12 +12 @@
         /*input datasets: */
         DataSet* nodes=NULL;
+        DataSet* vertices=NULL;
         DataSet* elements=NULL;
         DataSet* params=NULL;
@@ -28 +29 @@
         FetchData(&elements,ELEMENTS);
         FetchData(&nodes,NODESIN);
+        FetchData(&vertices,VERTICESIN);
         FetchData(&params,PARAMS);
 
         /*!Generate internal degree of freedom numbers: */
-        Dofx(&partition, &tpartition,elements,nodes,params); 
+        Dofx(&partition, &tpartition,elements,nodes, vertices, params); 
 
         /*partition and tpartition should be incremented by 1: */
@@ -39 +41 @@
         /*write output datasets: */
         WriteData(NODES,nodes);
+        WriteData(VERTICES,vertices);
         WriteData(PARTITION,partition);
         WriteData(TPARTITION,tpartition);
@@ -44 +47 @@
         /*Free ressources: */
         delete nodes;
+        delete vertices;
         delete elements;
         delete params;

issm/trunk/src/mex/Dof/Dof.h

@@ -19 +19 @@
 #define ELEMENTS (mxArray*)prhs[0]
 #define NODESIN (mxArray*)prhs[1]
-#define PARAMS (mxArray*)prhs[2]
+#define VERTICESIN (mxArray*)prhs[2]
+#define PARAMS (mxArray*)prhs[3]
 
 /* serial output macros: */
 #define NODES (mxArray**)&plhs[0]
-#define PARTITION (mxArray**)&plhs[1]
-#define TPARTITION (mxArray**)&plhs[2]
+#define VERTICES (mxArray**)&plhs[1]
+#define PARTITION (mxArray**)&plhs[2]
+#define TPARTITION (mxArray**)&plhs[3]
 
 /* serial arg counts: */
 #undef NLHS
-#define NLHS  3
+#define NLHS  4
 #undef NRHS
-#define NRHS  3
+#define NRHS  4
 
 

issm/trunk/src/mex/ModelProcessor/ModelProcessor.cpp

@@ -14 +14 @@
         DataSet* elements=NULL;
         DataSet* nodes=NULL;
+        DataSet* vertices=NULL;
         DataSet* constraints=NULL;
         DataSet* loads=NULL;
@@ -32 +33 @@
 
         /*Create elements, nodes and materials: */
-        CreateDataSets(&elements,&nodes,&materials,&constraints, &loads, &parameters, iomodel,MODEL);
+        CreateDataSets(&elements,&nodes,&vertices,&materials,&constraints, &loads, &parameters, iomodel,MODEL);
         
         /*Write output data: */
         WriteData(ELEMENTS,elements);
         WriteData(NODES,nodes);
+        WriteData(VERTICES,vertices);
         WriteData(CONSTRAINTS,constraints);
         WriteData(LOADS,loads);
@@ -47 +49 @@
         delete elements;
         delete nodes;
+        delete vertices;
         delete constraints;
         delete loads;

issm/trunk/src/mex/ModelProcessor/ModelProcessor.h

@@ -26 +26 @@
 #define ELEMENTS (mxArray**)&plhs[0]
 #define NODES (mxArray**)&plhs[1]
-#define CONSTRAINTS (mxArray**)&plhs[2]
-#define LOADS (mxArray**)&plhs[3]
-#define MATERIALS (mxArray**)&plhs[4]
-#define PARAMETERS (mxArray**)&plhs[5]
+#define VERTICES (mxArray**)&plhs[2]
+#define CONSTRAINTS (mxArray**)&plhs[3]
+#define LOADS (mxArray**)&plhs[4]
+#define MATERIALS (mxArray**)&plhs[5]
+#define PARAMETERS (mxArray**)&plhs[6]
                 
 /* serial arg counts: */
 #undef NLHS
-#define NLHS  6
+#define NLHS  7
 #undef NRHS
 #define NRHS  1