Changeset 848

Timestamp:

06/08/09 15:52:02 (15 years ago)

Author:

Eric.Larour

Message:

Added FieldExtrudex, FieldDepthAveragex, UpdateNodePositionsx modules

Location:

issm/trunk/src/c

Files:

• DataSet/DataSet.cpp (modified) (4 diffs)
• DataSet/DataSet.h (modified) (1 diff)
• FieldDepthAveragex (added)
• FieldDepthAveragex/FieldDepthAveragex.cpp (added)
• FieldDepthAveragex/FieldDepthAveragex.h (added)
• FieldExtrudex (added)
• FieldExtrudex/FieldExtrudex.cpp (added)
• FieldExtrudex/FieldExtrudex.h (added)
• Makefile.am (modified) (2 diffs)
• SlopeExtrudex (deleted)
• ThicknessExtrudex (deleted)
• UpdateNodePositionsx (added)
• UpdateNodePositionsx/UpdateNodePositionsx.cpp (added)
• UpdateNodePositionsx/UpdateNodePositionsx.h (added)
• VelocityDepthAveragex (deleted)
• VelocityExtrudex (deleted)
• issm.h (modified) (2 diffs)
• objects/Node.cpp (modified) (3 diffs)
• objects/Node.h (modified) (1 diff)
• objects/Penta.cpp (modified) (1 diff)
• objects/Penta.h (modified) (1 diff)
• parallel/diagnostic_core.cpp (modified) (3 diffs)
• parallel/prognostic_core.cpp (modified) (2 diffs)

issm/trunk/src/c/DataSet/DataSet.cpp

@@ -1299 +1299 @@
 
 }
-                
-void  DataSet::VelocityExtrude(Vec ug,double* ug_serial){
+
+void  DataSet::FieldExtrude(Vec field,double* field_serial,char* field_name, int collapse){
 
         vector<Object*>::iterator object;
         Penta* penta=NULL;
+        Node*  node=NULL;
 
         for ( object=objects.begin() ; object < objects.end(); object++ ){
@@ -1310 +1311 @@
 
                         penta=(Penta*)(*object);
-                        penta->VelocityExtrude(ug,ug_serial);
-
-                }
-        }
-
-}
-void  DataSet::ThicknessExtrude(Vec tg,double* tg_serial){
-
-        vector<Object*>::iterator object;
-        Penta* penta=NULL;
-
-        for ( object=objects.begin() ; object < objects.end(); object++ ){
-                
-                if((*object)->Enum()==PentaEnum()){
-
-                        penta=(Penta*)(*object);
-                        penta->ThicknessExtrude(tg,tg_serial);
-
-                }
-        }
-
-}
-void  DataSet::VelocityExtrudeAllElements(Vec ug,double* ug_serial){
-
-        vector<Object*>::iterator object;
-        Penta* penta=NULL;
-
-        for ( object=objects.begin() ; object < objects.end(); object++ ){
-                
-                if((*object)->Enum()==PentaEnum()){
-
-                        penta=(Penta*)(*object);
-                        penta->VelocityExtrudeAllElements(ug,ug_serial);
-
-                }
-        }
-
-}
-void  DataSet::VelocityDepthAverageAtBase(Vec ug,double* ug_serial){
+                        penta->FieldExtrude(field,field_serial,field_name,collapse);
+
+                }
+                if((*object)->Enum()==NodeEnum()){
+
+                        node=(Node*)(*object);
+                        node->FieldExtrude(field,field_serial,field_name);
+
+                }
+        }
+
+}
+
+void  DataSet::FieldDepthAverageAtBase(Vec field,double* field_serial,char* fieldname){
 
         vector<Object*>::iterator object;
@@ -1358 +1334 @@
 
                         node=(Node*)(*object);
-                        node->VelocityDepthAverageAtBase(ug,ug_serial);
-                }
-        }
-
-}
-void  DataSet::SlopeExtrude(Vec sg,double* sg_serial){
-
-        vector<Object*>::iterator object;
-        Penta* penta=NULL;
-
-        for ( object=objects.begin() ; object < objects.end(); object++ ){
-                
-                if((*object)->Enum()==PentaEnum()){
-
-                        penta=(Penta*)(*object);
-                        penta->SlopeExtrude(sg,sg_serial);
-
+                        node->FieldDepthAverageAtBase(field,field_serial,fieldname);
                 }
         }
@@ -1397 +1357 @@
 
 
+void  DataSet::UpdateNodePositions(double* thickness,double* bed){
+
+        vector<Object*>::iterator object;
+        Node* node=NULL;
+
+        for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+                if((*object)->Enum()==NodeEnum()){
+
+                        node=(Node*)(*object);
+                        node->UpdateNodePosition(thickness,bed);
+                }
+        }
+
+
+
+}

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

@@ -76 +76 @@
                 void  Gradj(Vec grad_g,double*  u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type);
                 void  Misfit(double* pJ, double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type);
-                void  VelocityExtrude(Vec ug,double* ug_serial);
-                void  ThicknessExtrude(Vec ug,double* ug_serial);
-                void  VelocityExtrudeAllElements(Vec ug,double* ug_serial);
-                void  VelocityDepthAverageAtBase(Vec ug,double* ug_serial);
-                void  SlopeExtrude(Vec sg,double* sg_serial);
+                void  FieldDepthAverageAtBase(Vec field,double* field_serial,char* fieldname);
                 int   DeleteObject(Object* object);
                 void  ComputePressure(Vec p_g);
                 int   FindResult(void* pvalue, char* name);
+                void  FieldExtrude(Vec field,double* field_serial,char* field_name, int collapse);
+                void  UpdateNodePositions(double* thickness,double* bed);
 
 };

issm/trunk/src/c/Makefile.am

@@ -264 +264 @@
                                         ./ProcessParamsx/ProcessParamsx.cpp\
                                         ./ProcessParamsx/ProcessParamsx.h\
-                                        ./ThicknessExtrudex/ThicknessExtrudex.cpp\
-                                        ./ThicknessExtrudex/ThicknessExtrudex.h\
-                                        ./VelocityExtrudex/VelocityExtrudex.cpp\
-                                        ./VelocityExtrudex/VelocityExtrudex.h\
-                                        ./VelocityDepthAveragex/VelocityDepthAveragex.cpp\
-                                        ./VelocityDepthAveragex/VelocityDepthAveragex.h\
-                                        ./SlopeExtrudex/SlopeExtrudex.cpp\
-                                        ./SlopeExtrudex/SlopeExtrudex.h
+                                        ./FieldDepthAveragex/FieldDepthAveragex.cpp\
+                                        ./FieldDepthAveragex/FieldDepthAveragex.h\
+                                        ./FieldExtrudex/FieldExtrudex.cpp\
+                                        ./FieldExtrudex/FieldExtrudex.h
 
 
@@ -532 +528 @@
                                         ./Solverx/Solverx.cpp\
                                         ./Solverx/Solverx.h\
-                                        ./ThicknessExtrudex/ThicknessExtrudex.cpp\
-                                        ./ThicknessExtrudex/ThicknessExtrudex.h\
                                         ./Mergesolutionfromftogx/Mergesolutionfromftogx.cpp\
                                         ./Mergesolutionfromftogx/Mergesolutionfromftogx.h\
                                         ./ProcessParamsx/ProcessParamsx.cpp\
                                         ./ProcessParamsx/ProcessParamsx.h\
-                                        ./VelocityExtrudex/VelocityExtrudex.cpp\
-                                        ./VelocityExtrudex/VelocityExtrudex.h\
-                                        ./VelocityDepthAveragex/VelocityDepthAveragex.cpp\
-                                        ./VelocityDepthAveragex/VelocityDepthAveragex.h\
-                                        ./SlopeExtrudex/SlopeExtrudex.cpp\
-                                        ./SlopeExtrudex/SlopeExtrudex.h\
+                                        ./FieldDepthAveragex/FieldDepthAveragex.cpp\
+                                        ./FieldDepthAveragex/FieldDepthAveragex.h\
+                                        ./FieldExtrudex/FieldExtrudex.cpp\
+                                        ./FieldExtrudex/FieldExtrudex.h\
                                         ./parallel/diagnostic_core.cpp\
                                         ./parallel/diagnostic_core_linear.cpp\

issm/trunk/src/c/issm.h

@@ -34 +34 @@
 #include "./UpdateFromInputsx/UpdateFromInputsx.h"
 #include "./UpdateGeometryx/UpdateGeometryx.h"
+#include "./UpdateNodePositionsx/UpdateNodePositionsx.h"
 #include "./PenaltySystemMatricesx/PenaltySystemMatricesx.h"
 #include "./Reducematrixfromgtofx/Reducematrixfromgtofx.h"
@@ -45 +46 @@
 #include "./Misfitx/Misfitx.h"
 #include "./ControlConstrainx/ControlConstrainx.h"
-#include "./VelocityExtrudex/VelocityExtrudex.h"
-#include "./VelocityDepthAveragex/VelocityDepthAveragex.h"
+#include "./FieldDepthAveragex/FieldDepthAveragex.h"
 #include "./GriddataMeshToGridx/GriddataMeshToGridx.h"
 #include "./ComputePressurex/ComputePressurex.h"
-#include "./SlopeExtrudex/SlopeExtrudex.h"
-#include "./ThicknessExtrudex/ThicknessExtrudex.h"
+#include "./FieldExtrudex/FieldExtrudex.h"
 
 

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

@@ -550 +550 @@
                 
 #undef __FUNCT__ 
-#define __FUNCT__ "Node::VelocityDepthAverageAtBase"
-void  Node::VelocityDepthAverageAtBase(Vec ug,double* ug_serial){
+#define __FUNCT__ "Node::FieldDepthAverageAtBase"
+void  Node::FieldDepthAverageAtBase(Vec field,double* field_serial,char* fieldname){
 
         /* node data: */
@@ -560 +560 @@
         Node* node=NULL;
         Node* upper_node=NULL;
-        double velocity2[2];
-        double velocity1[2];
-        double velocity_average[2];
-        double sum[2];
         double z1,z2,dz;
         double thickness;
 
+        /*Are we on the base, not on the surface, and not on a clone node?:*/
         
         if(onbed==1 & clone==0 &onsurface==0){
@@ -572 +569 @@
                 doflist1=GetDofList1();
 
-                /*this node is on the bed. We are going to:
-                 * follow the upper nodes until we reach the surface. At each upper node, we'll grab the 
-                 * velocity for this node, and add it to ug. We'll finish by 
-                 * we grab the velocity. Once we know the velocity, we follow the upper nodes, 
-                 * inserting the same velocity value into ug, until we reach the surface: */
-                sum[0]=0;
-                sum[1]=0;
-                thickness=0;
-
-                /*get dofs for this base node velocity: we know there are two dofs in ug_serial */
-                dofx=2*doflist1;
-                dofy=2*doflist1+1;
-
-                node=this;
-                for(;;){
-
-                        if (node->IsOnSurface())break;
-
-                        doflist1=node->GetDofList1();
+                /*this node is on the bed. We are going to, follow the upper nodes until we reach the surface. At each upper node, 
+                 * we'll grab the * field for this node, and add it to overall field: */
+
+                if(strcmp(fieldname,"velocity")==0){
+
+                        /*field is a velocity, 2 dofs per node: */
+                        double velocity2[2];
+                        double velocity1[2];
+                        double velocity_average[2];
+                        double sum[2];
+
+                        sum[0]=0;
+                        sum[1]=0;
+                        thickness=0;
+
+                        /*get dofs for this base node velocity: we know there are two dofs in field_serial */
+                        dofx=2*doflist1;
+                        dofy=2*doflist1+1;
+
+                        node=this;
+                        for(;;){
+
+                                if (node->IsOnSurface())break;
+
+                                doflist1=node->GetDofList1();
+                                
+                                velocity1[0]=field_serial[2*doflist1];
+                                velocity1[1]=field_serial[2*doflist1+1];
+                                z1=node->GetZ();
+
+                                upper_node=node->GetUpperNode();
+                                doflist1=upper_node->GetDofList1();
                         
-                        velocity1[0]=ug_serial[2*doflist1];
-                        velocity1[1]=ug_serial[2*doflist1+1];
-                        z1=node->GetZ();
-
-                        upper_node=node->GetUpperNode();
-                        doflist1=upper_node->GetDofList1();
-                
-                        velocity2[0]=ug_serial[2*doflist1];
-                        velocity2[1]=ug_serial[2*doflist1+1];
-                        z2=upper_node->GetZ();
-
-                        dz=(z2-z1);
-                        thickness+=dz;
-                        velocity_average[0]=(velocity1[0]+velocity2[0])/2.0;
-                        velocity_average[1]=(velocity1[1]+velocity2[1])/2.0;
-
-                        sum[0]+=velocity_average[0]*dz;
-                        sum[1]+=velocity_average[1]*dz;
-
-                        /* get next node: */
-                        node=node->GetUpperNode();
+                                velocity2[0]=field_serial[2*doflist1];
+                                velocity2[1]=field_serial[2*doflist1+1];
+                                z2=upper_node->GetZ();
+
+                                dz=(z2-z1);
+                                thickness+=dz;
+                                velocity_average[0]=(velocity1[0]+velocity2[0])/2.0;
+                                velocity_average[1]=(velocity1[1]+velocity2[1])/2.0;
+
+                                sum[0]+=velocity_average[0]*dz;
+                                sum[1]+=velocity_average[1]*dz;
+
+                                /* get next node: */
+                                node=node->GetUpperNode();
+                        }
+
+                        sum[0]=sum[0]/thickness;
+                        sum[1]=sum[1]/thickness;
+
+                        /* Plfield velocity_average*deltaH/H into base of field: */
+                        VecSetValues(field,1,&dofx,&sum[0],INSERT_VALUES);
+                        VecSetValues(field,1,&dofy,&sum[1],INSERT_VALUES);
                 }
-
-                sum[0]=sum[0]/thickness;
-                sum[1]=sum[1]/thickness;
-
-                /* Plug velocity_average*deltaH/H into base of ug: */
-                VecSetValues(ug,1,&dofx,&sum[0],INSERT_VALUES);
-                VecSetValues(ug,1,&dofy,&sum[1],INSERT_VALUES);
-        }
-}
+                else{
+                        /*field is regular, 1 dof per node: */
+                        double field2;
+                        double field1;
+                        double field_average;
+                        double sum;
+
+                        sum=0;
+                        thickness=0;
+
+                        /*get dofs for this base node velocity: we know there are two dofs in field_serial */
+                        dofx=doflist1;
+
+                        node=this;
+                        for(;;){
+
+                                if (node->IsOnSurface())break;
+
+                                doflist1=node->GetDofList1();
+                                
+                                field1=field_serial[doflist1];
+                                z1=node->GetZ();
+
+                                upper_node=node->GetUpperNode();
+                                doflist1=upper_node->GetDofList1();
+                        
+                                field2=field_serial[2*doflist1];
+                                z2=upper_node->GetZ();
+
+                                dz=(z2-z1);
+                                thickness+=dz;
+                                field_average=(field1+field2)/2.0;
+
+                                sum+=field_average*dz;
+
+                                /* get next node: */
+                                node=node->GetUpperNode();
+                        }
+
+                        sum=sum/thickness;
+
+                        /* Plug field_average*deltH/H into base of field: */
+                        VecSetValues(field,1,&dofx,&sum,INSERT_VALUES);
+                }
+        }
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Node::FieldExtrude"
+void  Node::FieldExtrude(Vec field,double* field_serial,char* field_name){
+                
+        /* node data: */
+        int   numberofdofspernode;
+        Node* node=NULL;
+        int   i;
+
+        /*Is this node on bed? :*/
+        if (onbed){
+
+                if (strcmp(field_name,"velocity")==0){
+
+                        /* node data: */
+                        const int    numdof=2;
+                        int          doflist[numdof];
+                        int          nodedofs[2];
+                        double       fieldnode[2];
+
+
+                        GetDofList(&doflist[0],&numberofdofspernode);
+
+                        //initilaize node
+                        node=this;
+                
+                        /*get field for this base node: */
+                        fieldnode[0]=field_serial[doflist[0]];
+                        fieldnode[1]=field_serial[doflist[1]];
+
+                        //go through all nodes which sit on top of this node, until we reach the surface, 
+                        //and plfield  field in field
+                        for(;;){
+
+                                node->GetDofList(&nodedofs[0],&numberofdofspernode);
+                                VecSetValues(field,1,&nodedofs[0],&fieldnode[0],INSERT_VALUES);
+                                VecSetValues(field,1,&nodedofs[1],&fieldnode[1],INSERT_VALUES);
+
+                                if (node->IsOnSurface())break;
+                                /*get next node: */
+                                node=node->GetUpperNode();
+                        }
+                } //if (strcmp(field_name,"velocity")==0)
+                else if ( 
+                                (strcmp(field_name,"thickness")==0) ||
+                                (strcmp(field_name,"surface")==0)  ||
+                                (strcmp(field_name,"bed")==0)  || 
+                                (strcmp(field_name,"slopex")==0)  ||
+                                (strcmp(field_name,"slopey")==0) 
+                                ){
+
+                        /* node data: */
+                        int          doflist;
+                        int          nodedofs;
+                        double       fieldnode;
+                                
+                        GetDofList(&doflist,&numberofdofspernode);
+
+                        //initilaize node
+                        node=this;
+
+                        /*get field for this bed node: */
+                        fieldnode=field_serial[doflist];
+
+                        //go through all nodes which sit on top of this node, until we reach the surface, 
+                        //and pltg  fieldnode in field:
+                        for(;;){
+
+                                node->GetDofList(&nodedofs,&numberofdofspernode);
+                                VecSetValues(field,1,&nodedofs,&fieldnode,INSERT_VALUES);
+
+                                if (node->IsOnSurface())break;
+                                /*get next node: */
+                                node=node->GetUpperNode();
+                        }
+                }
+                else throw ErrorException(__FUNCT__,exprintf("%s%s%s"," field ",field_name," not supported yet!"));
+
+        } //if (extrude)
+}
+                
+void  Node::UpdateNodePosition(double* thickness,double* bed){
+
+        int dof;
+
+        dof=this->GetDofList1();
+
+        if(onbed){
+
+                /*this node is on the bed, set its z position to the new bed: */
+                this->x[2]=bed[dof];
+        }
+        else if (onsurface){
+
+                /*this node is on the surface, set its z position to the new bed+ new thickness: */
+                this->x[2]=bed[dof]+thickness[dof];
+        }
+        else{
+                /*do nothing, we only update grids belonging to the exterior layers: */
+        }
+}

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

@@ -81 +81 @@
                 int   IsOnSurface();
                 void  FreezeDof(int dof);
-                void  VelocityDepthAverageAtBase(Vec ug,double* ug_serial);
+                void  FieldDepthAverageAtBase(Vec field,double* field_serial,char* fieldname);
                 int   IsOnShelf();
                 int   IsOnSheet();
+                void  FieldExtrude(Vec field,double* field_serial,char* field_name);
+                void  UpdateNodePosition(double* thickness,double* bed);
 };
 

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

@@ -1922 +1922 @@
 
 #undef __FUNCT__ 
-#define __FUNCT__ "Penta::VelocityExtrude"
-void  Penta::VelocityExtrude(Vec ug,double* ug_serial){
-
+#define __FUNCT__ "Penta::FieldExtrude"
+void  Penta::FieldExtrude(Vec field,double* field_serial,char* field_name, int iscollapsed){
+                
         /* node data: */
         const int    numgrids=6;
-        const int    numdof=2*numgrids;
-        int          doflist[numdof];
-        int          nodedofs[2];
         int          numberofdofspernode;
-        
         Node* node=NULL;
         int   i;
-        double velocity[2];
-                
-
-        if((collapse==1) && (onbed==1)){
-                        
-                GetDofList(&doflist[0],&numberofdofspernode);
-
-                /*this penta is a collapsed macayeal. For each node on the base of this penta, 
-                 * we grab the velocity. Once we know the velocity, we follow the upper nodes, 
-                 * inserting the same velocity value into ug, until we reach the surface: */
-                for(i=0;i<3;i++){
-        
-                        node=nodes[i]; //base nodes
-        
-                        /*get velocity for this base node: */
-                        velocity[0]=ug_serial[doflist[numberofdofspernode*i+0]];
-                        velocity[1]=ug_serial[doflist[numberofdofspernode*i+1]];
-
-                        //go througn all nodes which sit on top of this node, until we reach the surface, 
-                        //and plug  velocity in ug
-                        for(;;){
-
-                                node->GetDofList(&nodedofs[0],&numberofdofspernode);
-                                VecSetValues(ug,1,&nodedofs[0],&velocity[0],INSERT_VALUES);
-                                VecSetValues(ug,1,&nodedofs[1],&velocity[1],INSERT_VALUES);
-
-                                if (node->IsOnSurface())break;
-                                /*get next node: */
-                                node=node->GetUpperNode();
+        int   extrude=0;
+
+        /*Figure out if we should extrude for this element: */
+        if (iscollapsed){
+                /*From higher level, we are told to extrude only elements that have the collapse flag on: */
+                if (collapse)extrude=1;
+                else extrude=0;
+        }
+        else{
+                /*From higher level, we are told to extrude all elements: */
+                extrude=1;
+        }
+
+        /*Now, extrusion starts from the bed on, so double check this element is on 
+         * the bedrock: */
+        if(onbed==0)extrude=0;
+
+        /*Go on and extrude field: */
+        if (extrude){
+
+                if (strcmp(field_name,"velocity")==0){
+
+                        /* node data: */
+                        const int    numdof=2*numgrids;
+                        int          doflist[numdof];
+                        int          nodedofs[2];
+                        double       fieldel[2];
+
+
+                        GetDofList(&doflist[0],&numberofdofspernode);
+
+                        /*this penta is a collapsed macayeal. For each node on the base of this penta, 
+                         * we grab the field. Once we know the field, we follow the upper nodes, 
+                         * inserting the same field value into field, until we reach the surface: */
+                        for(i=0;i<3;i++){
+                
+                                node=nodes[i]; //base nodes
+                
+                                /*get field for this base node: */
+                                fieldel[0]=field_serial[doflist[numberofdofspernode*i+0]];
+                                fieldel[1]=field_serial[doflist[numberofdofspernode*i+1]];
+
+                                //go throfieldn all nodes which sit on top of this node, until we reach the surface, 
+                                //and plfield  field in field
+                                for(;;){
+
+                                        node->GetDofList(&nodedofs[0],&numberofdofspernode);
+                                        VecSetValues(field,1,&nodedofs[0],&fieldel[0],INSERT_VALUES);
+                                        VecSetValues(field,1,&nodedofs[1],&fieldel[1],INSERT_VALUES);
+
+                                        if (node->IsOnSurface())break;
+                                        /*get next node: */
+                                        node=node->GetUpperNode();
+                                }
                         }
+                } //if (strcmp(field_name,"velocity")==0)
+                else if ( 
+                                (strcmp(field_name,"thickness")==0) ||
+                                (strcmp(field_name,"surface")==0)  ||
+                                (strcmp(field_name,"bed")==0)  ||
+                                (strcmp(field_name,"slopex")==0)  ||
+                                (strcmp(field_name,"slopey")==0) 
+                                ){
+
+                        /* node data: */
+                        const int    numdof=1*numgrids;
+                        int          doflist[numdof];
+                        int          nodedofs;
+                        double       fieldel;
+                                
+                        GetDofList(&doflist[0],&numberofdofspernode);
+
+                        /*this penta is on the bed. For each node on the base of this penta, 
+                         * we grab the thickness. Once we know the thickness, we follow the upper nodes, 
+                         * inserting the same thickness value into tg, until we reach the surface: */
+                        for(i=0;i<3;i++){
+                
+                                node=nodes[i]; //base nodes
+                
+                                /*get velocity for this base node: */
+                                fieldel=field_serial[doflist[numberofdofspernode*i+0]];
+
+                                //go through all nodes which sit on top of this node, until we reach the surface, 
+                                //and pltg  fieldel in field:
+                                for(;;){
+
+                                        node->GetDofList(&nodedofs,&numberofdofspernode);
+                                        VecSetValues(field,1,&nodedofs,&fieldel,INSERT_VALUES);
+
+                                        if (node->IsOnSurface())break;
+                                        /*get next node: */
+                                        node=node->GetUpperNode();
+                                }
+                        }
+
                 }
-
-        }
-
-}
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Penta::VelocityExtrudeAllElements"
-void  Penta::VelocityExtrudeAllElements(Vec ug,double* ug_serial){
-
-        /* node data: */
-        const int    numgrids=6;
-        int          doflist[numgrids];
-        int          doflist1;
-        int          nodedofs[2];
-        
-        Node* node=NULL;
-        int   i;
-        double velocity[2];
-                
-
-        if(onbed==1){
-                        
-                GetDofList1(&doflist[0]);
-
-                /*this penta is on the bed. For each node on the base of this penta, 
-                 * we grab the velocity. Once we know the velocity, we follow the upper nodes, 
-                 * inserting the same velocity value into ug, until we reach the surface: 
-                 * the velocity has two ddls vx and vy*/
-                for(i=0;i<3;i++){
-        
-                        node=nodes[i]; //base nodes
-        
-                        /*get velocity for this base node: */
-                        velocity[0]=ug_serial[2*doflist[i]+0];
-                        velocity[1]=ug_serial[2*doflist[i]+1];
-
-                        //go througn all nodes which sit on top of this node, until we reach the surface, 
-                        //and plug  velocity in ug
-                        for(;;){
-
-                                doflist1=node->GetDofList1();
-                                nodedofs[0]=2*doflist1;
-                                nodedofs[1]=2*doflist1+1;
-                                VecSetValues(ug,1,&nodedofs[0],&velocity[0],INSERT_VALUES);
-                                VecSetValues(ug,1,&nodedofs[1],&velocity[1],INSERT_VALUES);
-
-                                if (node->IsOnSurface())break;
-                                /*get next node: */
-                                node=node->GetUpperNode();
-                        }
-                }
-
-        }
-
-}
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Penta::ThicknessExtrude"
-void  Penta::ThicknessExtrude(Vec tg,double* tg_serial){
-
-        /* node data: */
-        const int    numgrids=6;
-        const int    numdof=1*numgrids;
-        int          doflist[numdof];
-        int          nodedofs;
-        int          numberofdofspernode;
-        
-        Node* node=NULL;
-        int   i;
-        double thickness;
-                
-
-        if(onbed==1){
-                        
-                GetDofList(&doflist[0],&numberofdofspernode);
-
-                /*this penta is on the bed. For each node on the base of this penta, 
-                 * we grab the thickness. Once we know the thickness, we follow the upper nodes, 
-                 * inserting the same thickness value into tg, until we reach the surface: */
-                for(i=0;i<3;i++){
-        
-                        node=nodes[i]; //base nodes
-        
-                        /*get velocity for this base node: */
-                        thickness=tg_serial[doflist[numberofdofspernode*i+0]];
-
-                        //go through all nodes which sit on top of this node, until we reach the surface, 
-                        //and pltg  velocity in tg
-                        for(;;){
-
-                                node->GetDofList(&nodedofs,&numberofdofspernode);
-                                VecSetValues(tg,1,&nodedofs,&thickness,INSERT_VALUES);
-
-                                if (node->IsOnSurface())break;
-                                /*get next node: */
-                                node=node->GetUpperNode();
-                        }
-                }
-
-        }
-
-}
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Penta::SlopeExtrude"
-void  Penta::SlopeExtrude(Vec sg,double* sg_serial){
-
-        /* node data: */
-        const int    numgrids=6;
-        const int    NDOF1=1;
-        const int    numdof=NDOF1*numgrids;
-        int          doflist[numdof];
-        int          nodedof;
-        int          numberofdofspernode;
-        
-        Node* node=NULL;
-        int   i;
-        double slope;
-                
-
-        if(onbed==1){
-                        
-                GetDofList(&doflist[0],&numberofdofspernode);
-
-                if(numberofdofspernode!=1)throw ErrorException(__FUNCT__," slope can only be extruded on 1 dof per node");
-
-                /*For each node on the base of this penta,  we grab the slope. Once we know the slope, we follow the upper nodes, 
-                 * inserting the same slope value into sg, until we reach the surface: */
-                
-                for(i=0;i<3;i++){
-        
-                        node=nodes[i]; //base nodes
-        
-                        /*get velocity for this base node: */
-                        slope=sg_serial[doflist[i]];
-
-                        //go throsgn all nodes which sit on top of this node, until we reach the surface, 
-                        //and plsg  slope in sg
-                        for(;;){
-
-                                node->GetDofList(&nodedof,&numberofdofspernode);
-                                VecSetValues(sg,1,&nodedof,&slope,INSERT_VALUES);
-
-                                if (node->IsOnSurface())break;
-                                /*get next node: */
-                                node=node->GetUpperNode();
-                        }
-                }
-
-        }
-
-}
-
+                else throw ErrorException(__FUNCT__,exprintf("%s%s%s"," field ",field_name," not supported yet!"));
+
+        } //if (extrude)
+}
 
 #undef __FUNCT__ 

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

@@ -111 +111 @@
                 void  GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_l1l2l3l4);
                 void  GetNodalFunctions(double* l1l2l3l4l5l6, double* gauss_l1l2l3l4);
-                void  VelocityExtrudeAllElements(Vec ug,double* ug_serial);
-                void  VelocityExtrude(Vec ug,double* ug_serial);
-                void  ThicknessExtrude(Vec ug,double* ug_serial);
-                void  SlopeExtrude(Vec sg,double* sg_serial);
+                void  FieldExtrude(Vec field,double* field_serial,char* field_name, int iscollapsed);
                 void  ComputePressure(Vec p_g);
                 void  CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);

issm/trunk/src/c/parallel/diagnostic_core.cpp

@@ -85 +85 @@
                 
                         if(debug)_printf_("%s\n","extruding slopes in 3d...");
-                        SlopeExtrudex( slopex, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials);
-                        SlopeExtrudex( slopey, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials);
+                        FieldExtrudex( slopex, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials,"slopex",0);
+                        FieldExtrudex( slopey, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials,"slopex",0);
                 }
 
@@ -122 +122 @@
 
                 if(debug)_printf_("%s\n"," extruding horizontal velocities...");
-                VecDuplicatePatch(&ug_horiz,ug); VelocityExtrudex( ug_horiz,fem_dh->elements,fem_dh->nodes, fem_dh->loads,fem_dh-> materials);
+                VecDuplicatePatch(&ug_horiz,ug); FieldExtrudex( ug_horiz,fem_dh->elements,fem_dh->nodes, fem_dh->loads,fem_dh-> materials,"velocity",1);
 
                 if(debug)_printf_("%s\n"," computing vertical velocities...");
@@ -145 +145 @@
                         diagnostic_core_linear(&slopex,fem_sl,inputs,SlopeComputeAnalysisEnum(),BedXAnalysisEnum());
                         diagnostic_core_linear(&slopey,fem_sl,inputs,SlopeComputeAnalysisEnum(),BedYAnalysisEnum());
-                        SlopeExtrudex( slopex, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials);
-                        SlopeExtrudex( slopey, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials);
+                        FieldExtrudex( slopex, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials,"slopex",0);
+                        FieldExtrudex( slopey, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials,"slopey",0);
 
                         inputs->Add("bedslopex",slopex,numberofdofspernode_sl,numberofnodes);

issm/trunk/src/c/parallel/prognostic_core.cpp

@@ -39 +39 @@
         _printf_("depth averaging velocity...");
         u_g=inputs->Get("velocity",&dofs[0],2); //take (vx,vy) from inputs velocity
-        VelocityDepthAveragex( u_g, fem->elements,fem->nodes, fem->loads, fem->materials);
+        FieldDepthAveragex( u_g, fem->elements,fem->nodes, fem->loads, fem->materials,"velocity");
         inputs->Add("velocity_average",u_g,2,numberofnodes);
         
@@ -46 +46 @@
 
         _printf_("extrude computed thickness on all layers:\n");
-        ThicknessExtrudex( h_g, fem->elements,fem->nodes, fem->loads, fem->materials);
+        FieldExtrudex( h_g, fem->elements,fem->nodes, fem->loads, fem->materials,"thickness",0);
 
         /*Plug results into output dataset: */