Changeset 483

Timestamp:

05/18/09 16:55:54 (16 years ago)

Author:

Eric.Larour

Message:

Finished thermal and melting elemnets, as well as thermal constraints

Location:

issm/trunk/src/c

Files:

• DataSet/DataSet.cpp (modified) (3 diffs)
• Makefile.am (modified) (2 diffs)
• ModelProcessorx/CreateDataSets.cpp (modified) (1 diff)
• ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (modified) (2 diffs)
• ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp (modified) (2 diffs)
• ModelProcessorx/Melting (added)
• ModelProcessorx/Melting/CreateConstraintsMelting.cpp (added)
• ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp (added)
• ModelProcessorx/Melting/CreateLoadsMelting.cpp (added)
• ModelProcessorx/Melting/CreateParametersMelting.cpp (added)
• ModelProcessorx/Model.h (modified) (1 diff)
• ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp (modified) (2 diffs)
• ModelProcessorx/Thermal (added)
• ModelProcessorx/Thermal/CreateConstraintsThermal.cpp (added)
• ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp (added)
• ModelProcessorx/Thermal/CreateLoadsThermal.cpp (added)
• ModelProcessorx/Thermal/CreateParametersThermal.cpp (added)
• objects/Friction.cpp (modified) (2 diffs)
• objects/Matpar.cpp (modified) (1 diff)
• objects/Matpar.h (modified) (1 diff)
• objects/Pengrid.cpp (modified) (15 diffs)
• objects/Pengrid.h (modified) (3 diffs)
• objects/Penpair.cpp (modified) (1 diff)
• objects/Penpair.h (modified) (1 diff)
• objects/Penta.cpp (modified) (7 diffs)
• objects/Penta.h (modified) (1 diff)
• objects/Tria.cpp (modified) (59 diffs)
• objects/Tria.h (modified) (3 diffs)

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

@@ -1022 +1022 @@
 int   DataSet::RiftIsPresent(){
 
-        return 0; //for now
-
+        int i;
+        
+        vector<Object*>::iterator object;
+        Penpair* penpair=NULL;
+        int found=0;
+        int mpi_found=0;
+
+        /*go though loads, and figure out if one of the loads is a PenPair with numdof=2: */
+        for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+                if((*object)->Enum()==PenpairEnum()){
+
+                        penpair=(Penpair*)(*object);
+                        if (penpair->GetNumDofs()==2){
+                                found=1;
+                                break;
+                        }
+                }
+        }
+
+        #ifdef _PARALLEL_
+        MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
+        MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD);                
+        found=mpi_found;
+        #endif
+
+        return found;
 }
 
@@ -1038 +1063 @@
 int   DataSet::MeltingIsPresent(){
 
-
-        return 0; //for now
-        
+        int found=0;
+        int mpi_found=0;
+
+        vector<Object*>::iterator object;
+
+        for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+                if((*object)->Enum()==PengridEnum()){
+                        found=1;
+                        break;
+                }
+        }       
+        
+        #ifdef _PARALLEL_
+        MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
+        MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD);                
+        found=mpi_found;
+        #endif
+
+        return found;
 }
 
@@ -1047 +1089 @@
 void  DataSet::MeltingConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type){
 
-        throw ErrorException(__FUNCT__," not implemented yet!");
-
+        /* generic object pointer: */
+        vector<Object*>::iterator object;
+        Pengrid* pengrid=NULL;
+
+        int unstable;
+        int num_unstable_constraints=0;
+        int sum_num_unstable_constraints=0;
+
+        num_unstable_constraints=0;     
+
+        /*Enforce constraints: */
+        for ( object=objects.begin() ; object < objects.end(); object++ ){
+                
+                if((*object)->Enum()==PengridEnum()){
+
+                        pengrid=(Pengrid*)(*object);
+                        pengrid->PenaltyConstrain(&unstable,inputs,analysis_type,sub_analysis_type);
+
+                        num_unstable_constraints+=unstable;
+                }
+        }
+
+        #ifdef _PARALLEL_
+        MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
+        MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD);                
+        num_unstable_constraints=sum_num_unstable_constraints;
+        #endif
+
+        /*Assign output pointers: */
+        *pnum_unstable_constraints=num_unstable_constraints;
 }
 

issm/trunk/src/c/Makefile.am

@@ -188 +188 @@
                                         ./ModelProcessorx/SlopeCompute/CreateLoadsSlopeCompute.cpp\
                                         ./ModelProcessorx/Control/CreateParametersControl.cpp\
+                                        ./ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp\
+                                        ./ModelProcessorx/Thermal/CreateConstraintsThermal.cpp\
+                                        ./ModelProcessorx/Thermal/CreateLoadsThermal.cpp\
+                                        ./ModelProcessorx/Thermal/CreateParametersThermal.cpp\
+                                        ./ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp\
+                                        ./ModelProcessorx/Melting/CreateConstraintsMelting.cpp\
+                                        ./ModelProcessorx/Melting/CreateLoadsMelting.cpp\
+                                        ./ModelProcessorx/Melting/CreateParametersMelting.cpp\
                                         ./Dofx/Dofx.h\
                                         ./Dofx/Dofx.cpp\
@@ -434 +442 @@
                                         ./ModelProcessorx/SlopeCompute/CreateLoadsSlopeCompute.cpp\
                                         ./ModelProcessorx/Control/CreateParametersControl.cpp\
+                                        ./ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp\
+                                        ./ModelProcessorx/Thermal/CreateConstraintsThermal.cpp\
+                                        ./ModelProcessorx/Thermal/CreateLoadsThermal.cpp\
+                                        ./ModelProcessorx/Thermal/CreateParametersThermal.cpp\
+                                        ./ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp\
+                                        ./ModelProcessorx/Melting/CreateConstraintsMelting.cpp\
+                                        ./ModelProcessorx/Melting/CreateLoadsMelting.cpp\
+                                        ./ModelProcessorx/Melting/CreateParametersMelting.cpp\
                                         ./Dofx/Dofx.h\
                                         ./Dofx/Dofx.cpp\

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

@@ -73 +73 @@
         else if (strcmp(model->analysis_type,"thermal")==0){
 
-                if ((strcmp(model->sub_analysis_type,"steady")==0) || (strcmp(model->sub_analysis_type,"transient")==0)){
-                
-                        //CreateElementsNodesAndMaterialsThermal(pelements,pnodes,pmaterials, model,model_handle);
-                        //CreateConstraintsThermal(pconstraints,model,model_handle);
-                        //CreateLoadsThermal(ploads,model,model_handle);
+                CreateElementsNodesAndMaterialsThermal(pelements,pnodes,pmaterials, model,model_handle);
+                CreateConstraintsThermal(pconstraints,model,model_handle);
+                CreateLoadsThermal(ploads,model,model_handle);
                                         
-                }
-                else if (strcmp(model->sub_analysis_type,"melting")==0){
+        }
+        else if (strcmp(model->analysis_type,"melting")==0){
                         
-                        //CreateElementsNodesAndMaterialsMelting(pelements,pnodes,pmaterials, model,model_handle);
-                        //CreateConstraintsMelting(pconstraints,model,model_handle);
-                        //CreateLoadsMelting(ploads,model,model_handle);
-
-                }
+                CreateElementsNodesAndMaterialsMelting(pelements,pnodes,pmaterials, model,model_handle);
+                CreateConstraintsMelting(pconstraints,model,model_handle);
+                CreateLoadsMelting(ploads,model,model_handle);
         }
         else if (strcmp(model->analysis_type,"prognostic")==0){

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

@@ -62 +62 @@
         double tria_melting[3];
         double tria_accumulation[3];
+        double tria_geothermalflux[3];
         int    tria_friction_type;
         double tria_p;
@@ -293 +294 @@
 
                         /*Create tria element using its constructor:*/
-                        tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot);
+                        tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_geothermalflux,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot);
 
                         /*Add tria element to elements dataset: */

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

@@ -43 +43 @@
         /*pengrid intermediary data: */
         int pengrid_id;
+        int pengrid_mparid;
         int pengrid_node_id;
         int pengrid_dof;
@@ -141 +142 @@
                         pengrid_dof=1;
                         pengrid_penalty_offset=model->penalty_offset;
+                        pengrid_mparid=model->numberofelements+1;//refers to the corresponding parmat property card
 
-                        pengrid= new Pengrid(pengrid_id, pengrid_node_id,pengrid_dof, pengrid_active, pengrid_penalty_offset,pengrid_thermal_steadystate);
+                        pengrid= new Pengrid(pengrid_id, pengrid_node_id,pengrid_mparid,pengrid_dof, pengrid_active, pengrid_penalty_offset,pengrid_thermal_steadystate);
                         
                         loads->AddObject(pengrid);

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

@@ -210 +210 @@
         void    CreateParametersControl(DataSet** pparameters,Model* model,ConstDataHandle model_handle);
 
+        /*thermal:*/
+        void    CreateElementsNodesAndMaterialsThermal(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle);
+        void    CreateConstraintsThermal(DataSet** pconstraints,Model* model,ConstDataHandle model_handle);
+        void    CreateLoadsThermal(DataSet** ploads, Model* model, ConstDataHandle model_handle);
+        void    CreateParametersThermal(DataSet** pparameters,Model* model,ConstDataHandle model_handle);
+
+        /*melting:*/
+        void    CreateElementsNodesAndMaterialsMelting(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle);
+        void    CreateConstraintsMelting(DataSet** pconstraints,Model* model,ConstDataHandle model_handle);
+        void    CreateLoadsMelting(DataSet** ploads, Model* model, ConstDataHandle model_handle);
+        void    CreateParametersMelting(DataSet** pparameters,Model* model,ConstDataHandle model_handle);
+
 
 #endif  /* _MODEL_H */

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

@@ -58 +58 @@
         double tria_melting[3];
         double tria_accumulation[3];
+        double tria_geothermalflux[3];
         int    tria_friction_type;
         double tria_p;
@@ -226 +227 @@
 
                         /*Create tria element using its constructor:*/
-                        tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot);
+                        tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_geothermalflux,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot);
 
                         /*Add tria element to elements dataset: */

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

@@ -72 +72 @@
         double  velocity_x[numgrids];
         double  velocity_y[numgrids];
+        double  velocity_z[numgrids];
         double  velocity_mag[numgrids];
 
@@ -90 +91 @@
 
                 //We need the velocity magnitude to evaluate the basal stress:
-                velocity_x[i]=*(friction->velocities+2*i+0); //velocities of size numgridsx2
-                velocity_y[i]=*(friction->velocities+2*i+1);
-                velocity_mag[i]=sqrt(pow(velocity_x[i],2)+pow(velocity_y[i],2));
+                if(strcmp(friction->element_type,"2d")){
+                        velocity_x[i]=*(friction->velocities+2*i+0); //velocities of size numgridsx2
+                        velocity_y[i]=*(friction->velocities+2*i+1);
+                        velocity_mag[i]=sqrt(pow(velocity_x[i],2)+pow(velocity_y[i],2));
+                }
+                else{
+                        velocity_x[i]=*(friction->velocities+3*i+0); //velocities of size numgridsx3
+                        velocity_y[i]=*(friction->velocities+3*i+1);
+                        velocity_z[i]=*(friction->velocities+3*i+2);
+                        velocity_mag[i]=sqrt(pow(velocity_x[i],2)+pow(velocity_y[i],2)+pow(velocity_z[i],2));
+                }
         
                 alpha2[i]=pow(friction->K[i],2)*pow(Neff[i],r)*pow(velocity_mag[i],(s-1));

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

@@ -172 +172 @@
         return new Matpar(*this); 
 }
+
+double Matpar::GetMixedLayerCapacity(){
+        return mixed_layer_capacity;
+}
+
+double Matpar::GetThermalExchangeVelocity(){
+        return thermal_exchange_velocity;
+}
+
+double Matpar::GetHeatCapacity(){
+        return heatcapacity;
+}
+                
+double Matpar::GetThermalConductivity(){
+        return thermalconductivity;
+}
+                
+double Matpar::GetLatentHeat(){
+        return latentheat;
+}
+double Matpar::GetBeta(){
+        return beta;
+}
+double Matpar::GetMeltingPoint(){
+        return meltingpoint;
+}
+

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

@@ -44 +44 @@
                 double GetRhoIce();
                 double GetRhoWater();
+                double GetMixedLayerCapacity();
+                double GetThermalExchangeVelocity();
+                double GetHeatCapacity();
+                double GetThermalConductivity();
+                double GetLatentHeat();
+                double GetBeta();
+                double GetMeltingPoint();
                 Object* copy();
 

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

@@ -22 +22 @@
 }
 
-Pengrid::Pengrid(int    pengrid_id, int pengrid_node_id,int pengrid_dof, int pengrid_active, double pengrid_penalty_offset,int pengrid_thermal_steadystate){
+Pengrid::Pengrid(int    pengrid_id, int pengrid_mparid,int pengrid_node_id,int pengrid_dof, int pengrid_active, double pengrid_penalty_offset,int pengrid_thermal_steadystate){
         
         id=pengrid_id;
+        mparid=pengrid_mparid;
         dof=pengrid_dof;
         active=pengrid_active;
@@ -33 +34 @@
         node_offset=UNDEF;
         node=NULL;
+        matpar=NULL;
+        matpar_offset=UNDEF;
 
         return;
@@ -45 +48 @@
         printf("Pengrid:\n");
         printf("   id: %i\n",id);
+        printf("   mparid: %i\n",mparid);
         printf("   dof: %i\n",dof);
         printf("   active: %i\n",active);
@@ -51 +55 @@
         printf("   node_id: [%i]\n",node_id);
         printf("   node_offset: [%i]\n",node_offset);
+        printf("   matpar_offset=%i\n",matpar_offset);
         
         if(node)node->Echo();
+        if(matpar)matpar->Echo();
         return;
 }
@@ -72 +78 @@
         /*marshall Pengrid data: */
         memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
+        memcpy(marshalled_dataset,&mparid,sizeof(mparid));marshalled_dataset+=sizeof(mparid);
         memcpy(marshalled_dataset,&dof,sizeof(dof));marshalled_dataset+=sizeof(dof);
         memcpy(marshalled_dataset,&active,sizeof(active));marshalled_dataset+=sizeof(active);
@@ -78 +85 @@
         memcpy(marshalled_dataset,&node_id,sizeof(node_id));marshalled_dataset+=sizeof(node_id);
         memcpy(marshalled_dataset,&node_offset,sizeof(node_offset));marshalled_dataset+=sizeof(node_offset);
+        memcpy(marshalled_dataset,&matpar,sizeof(matpar));marshalled_dataset+=sizeof(matpar);
+        memcpy(marshalled_dataset,&matpar_offset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset);
 
         *pmarshalled_dataset=marshalled_dataset;
@@ -86 +95 @@
 
         return sizeof(id)+
+                sizeof(mparid)+
                 sizeof(dof)+
                 sizeof(active)+
@@ -92 +102 @@
                 sizeof(node_id)+
                 sizeof(node_offset)+
+                +sizeof(matpar)
+                +sizeof(matpar_offset)+
                 sizeof(int); //sizeof(int) for enum type
 }
@@ -117 +129 @@
         memcpy(&node_id,marshalled_dataset,sizeof(node_id));marshalled_dataset+=sizeof(node_id);
         memcpy(&node_offset,marshalled_dataset,sizeof(node_offset));marshalled_dataset+=sizeof(node_offset);
+        memcpy(&matpar,marshalled_dataset,sizeof(matpar));marshalled_dataset+=sizeof(matpar);
+        memcpy(&matpar_offset,marshalled_dataset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset);
+
 
         node=NULL;
+        matpar=NULL;
 
         /*return: */
@@ -143 +159 @@
 
         DataSet* nodesin=NULL;
+        DataSet* materialsin=NULL;
 
         /*Recover pointers :*/
         nodesin=(DataSet*)pnodesin;
+        materialsin=(DataSet*)pmaterialsin;
 
         /*Link this load with its nodes: */
         ResolvePointers((Object**)&node,&node_id,&node_offset,1,nodesin);
-
+        ResolvePointers((Object**)&matpar,&mparid,&matpar_offset,1,materialsin);
 }
 
@@ -184 +202 @@
 
                 PenaltyCreateKMatrixDiagnosticStokes( Kgg,inputs,kmax,analysis_type,sub_analysis_type);
+        }
+        else if (analysis_type==ThermalAnalysisEnum()){
+                
+                PenaltyCreateKMatrixThermal( Kgg,inputs,kmax,analysis_type,sub_analysis_type);
+                
+        }
+        else if (analysis_type==MeltingAnalysisEnum()){
+                        
+                PenaltyCreateKMatrixMelting( Kgg,inputs,kmax,analysis_type,sub_analysis_type);
 
         }
@@ -202 +229 @@
         int       numberofdofspernode;
 
-        int dofs1=0;
-        int dofs2=1;
+        int dofs1[1]={0};
+        int dofs2[1]={1};
         double slope[2];
         int found=0;
@@ -217 +244 @@
         
         /*recover slope: */
-        found=inputs->Recover("bedslopex",&slope[0],1,&dofs1,numgrids,(void**)&node);
+        found=inputs->Recover("bedslopex",&slope[0],1,dofs1,numgrids,(void**)&node);
         if(!found)throw ErrorException(__FUNCT__," bedslopex needed in inputs!");
-        found=inputs->Recover("bedslopey",&slope[1],1,&dofs2,numgrids,(void**)&node);
+        found=inputs->Recover("bedslopey",&slope[1],1,dofs2,numgrids,(void**)&node);
         if(!found)throw ErrorException(__FUNCT__," bedslopey needed in inputs!");
 
@@ -232 +259 @@
 
 #undef __FUNCT__ 
+#define __FUNCT__ "Pengrid::PenaltyCreateKMatrixThermal"
+void  Pengrid::PenaltyCreateKMatrixThermal(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
+
+        int found=0;
+        
+        const int numgrids=1;
+        const int NDOF1=1;
+        const int numdof=numgrids*NDOF1;
+        double Ke[numdof][numdof];
+        int       doflist[numdof];
+        int       numberofdofspernode;
+
+        ParameterInputs* inputs=NULL;
+
+        /*recover pointers: */
+        inputs=(ParameterInputs*)vinputs;
+
+
+        if(!active)return;
+
+        /*Get dof list: */
+        GetDofList(&doflist[0],&numberofdofspernode);
+
+        Ke[0][0]=kmax*pow(10,penalty_offset);
+        
+        /*Add Ke to global matrix Kgg: */
+        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Pengrid::PenaltyCreateKMatrixMelting"
+void  Pengrid::PenaltyCreateKMatrixMelting(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
+
+        int found=0;
+        const int numgrids=1;
+        const int NDOF1=1;
+        const int numdof=numgrids*NDOF1;
+        double Ke[numdof][numdof]={0.0};
+        int     dofs1[1]={0};
+        int       doflist[numdof];
+        int      numberofdofspernode;
+        double  meltingpoint;
+
+        double pressure;
+        double temperature;
+        double beta,t_pmp;
+
+        ParameterInputs* inputs=NULL;
+
+        /*recover pointers: */
+        inputs=(ParameterInputs*)vinputs;
+
+        found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
+        if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
+
+        found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node);
+        if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!");
+
+        /*Get dof list: */
+        GetDofList(&doflist[0],&numberofdofspernode);
+        
+        //Compute pressure melting point
+        meltingpoint=matpar->GetMeltingPoint();
+        beta=matpar->GetBeta();
+        t_pmp=meltingpoint-beta*pressure;
+
+        //Add penalty load
+        if (temperature<t_pmp){ //If T<Tpmp, there must be no melting. Therefore, melting should be  constrained to 0 when T<Tpmp, instead of using spcs, use penalties
+                Ke[0][0]=kmax*pow(10,penalty_offset);
+        }
+        
+        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
+}
+
+#undef __FUNCT__ 
 #define __FUNCT__ "Pengrid::PenaltyCreatePVector"
 void  Pengrid::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
 
-        /*No penalty applied, do nothing: */
-        return;
+        if (analysis_type==ThermalAnalysisEnum()){
+                
+                PenaltyCreatePVectorThermal( pg,inputs,kmax,analysis_type,sub_analysis_type);
+                
+        }
+        else if (analysis_type==MeltingAnalysisEnum()){
+                        
+                PenaltyCreatePVectorMelting( pg,inputs,kmax,analysis_type,sub_analysis_type);
+
+        }
+        else{
+                throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet"));
+        }
 
 }
@@ -259 +372 @@
         *pnumberofdofspernode=numberofdofspernode;
 }
+
+void  Pengrid::PenaltyCreatePVectorThermal(Vec pg, void* vinputs, double kmax,int analysis_type,int sub_analysis_type){
+
+        const int numgrids=1;
+        const int NDOF1=1;
+        const int numdof=numgrids*NDOF1;
+        int       doflist[numdof];
+        double  P_terms[numdof]={0.0};
+        int    numberofdofspernode;
+        int    found=0;
+        double pressure;
+        int dofs1[1]={0};
+        double meltingpoint;
+        double beta;
+        double t_pmp;
+
+        ParameterInputs* inputs=NULL;
+
+        /*recover pointers: */
+        inputs=(ParameterInputs*)vinputs;
+
+        if(!active)return;
+
+        /*Get dof list: */
+        GetDofList(&doflist[0],&numberofdofspernode);
+        
+                //First recover pressure
+                found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
+                if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
+
+                //Compute pressure melting point
+                meltingpoint=matpar->GetMeltingPoint();
+                beta=matpar->GetBeta();
+                t_pmp=meltingpoint-beta*pressure;
+
+                //Add penalty load
+                P_terms[0]=kmax*pow(10,penalty_offset)*t_pmp;
+
+                /*Add P_terms to global vector pg: */
+                VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
+}
+
+void  Pengrid::PenaltyCreatePVectorMelting(Vec pg, void* vinputs, double kmax,int analysis_type,int sub_analysis_type){
+        
+        const int numgrids=1;
+        const int NDOF1=1;
+        const int numdof=numgrids*NDOF1;
+        int       doflist[numdof];
+        double  P_terms[numdof]={0.0};
+        int numberofdofspernode;
+        int    found=0;
+        int dofs1[1]={0};
+        double pressure;
+        double temperature;
+        double melting_offset;
+        double meltingpoint;
+        double beta, heatcapacity;
+        double latentheat;
+        double t_pmp;
+        double dt;
+
+        ParameterInputs* inputs=NULL;
+
+        /*recover pointers: */
+        inputs=(ParameterInputs*)vinputs;
+
+        //First recover pressure,melting offset and temperature vectors
+        found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
+        if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
+
+        found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node);
+        if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!");
+
+        found=inputs->Recover("melting_offset",&melting_offset);
+        if(!found)throw ErrorException(__FUNCT__," could not find melting_offset in inputs!");
+
+        found=inputs->Recover("dt",&dt);
+        if((!found) && (sub_analysis_type==TransientAnalysisEnum()))throw ErrorException(__FUNCT__," could not find dt in inputs!");
+
+
+        meltingpoint=matpar->GetMeltingPoint();
+        beta=matpar->GetBeta();
+        heatcapacity=matpar->GetHeatCapacity();
+        latentheat=matpar->GetLatentHeat();
+
+        //Compute pressure melting point
+        t_pmp=meltingpoint-beta*pressure;
+
+        //Add penalty load
+        //This time, the penalty must have the same value as the one used for the thermal computation
+        //so that the corresponding melting can be computed correctly
+        //In the thermal computation, we used kmax=melting_offset, and the same penalty_offset
+        if (temperature<t_pmp){ //%no melting
+                P_terms[0]=0;
+        }
+        else{
+                if (sub_analysis_type==SteadyAnalysisEnum()){
+                        P_terms[0]=melting_offset*pow(10,penalty_offset)*(temperature-t_pmp);
+                }
+                else{
+                        P_terms[0]=melting_offset*pow(10,penalty_offset)*(temperature-t_pmp)/dt;
+                }
+        }
+        /*Add P_terms to global vector pg: */
+        VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
+}
+                
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Pengrid::PenaltyConstrain"
+void  Pengrid::PenaltyConstrain(int* punstable,void* vinputs,int analysis_type,int sub_analysis_type){
+
+        //   The penalty is stable if it doesn't change during to successive iterations.   
+
+        int found=0;
+        const int numgrids=1;
+
+
+        double pressure;
+        double temperature;
+        double beta,t_pmp;
+        double meltingpoint;
+        int    new_active;
+        int* dofs1={0};
+        int  unstable=0;
+        
+        ParameterInputs* inputs=NULL;
+
+        /*recover pointers: */
+        inputs=(ParameterInputs*)vinputs;
+
+
+        //First recover beta, pressure and temperature vectors;
+        found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
+        if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
+
+        found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node);
+        if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!");
+
+
+        //Compute pressure melting point
+        meltingpoint=matpar->GetMeltingPoint();
+        beta=matpar->GetBeta();
+
+        t_pmp=meltingpoint-beta*pressure;
+
+        //Figure out if temperature is over melting_point, in which case, this penalty needs to be activated.
+
+        if (temperature>t_pmp){
+                new_active=1;
+        }
+        else{
+                new_active=0;
+        }
+
+
+        //Figure out stability of this penalty
+        if (active==new_active){
+                unstable=0;
+        }
+        else{
+                unstable=1;
+        }
+        
+        //Set penalty flag
+        active=new_active;
+
+        //*Assign output pointers:*/
+        *punstable=unstable;
+}

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

@@ -8 +8 @@
 #include "./Node.h"
 #include "./Element.h"
+#include "./Matpar.h"
 
 class Element;
@@ -25 +26 @@
                 int     node_offset;
 
+                int mparid;
+                Matpar* matpar; 
+                int   matpar_offset;
+
         public:
 
                 Pengrid();
-                Pengrid(int     id, int node_id,int dof, int active, double penalty_offset,int thermal_steadystate);
+                Pengrid(int     id, int node_id,int mparid,int dof, int active, double penalty_offset,int thermal_steadystate);
                 ~Pengrid();
 
@@ -49 +54 @@
                 void  GetDofList(int* doflist,int* pnumberofdofspernode);
                 Object* copy();
-
+                void  PenaltyCreateKMatrixThermal(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type);
+                void  PenaltyCreateKMatrixMelting(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type);
+                void  MatparConfiguration(Matpar* matpar,int matpar_offset);
+                void  PenaltyCreatePVectorThermal(Vec pg, void* inputs, double kmax,int analysis_type,int sub_analysis_type);
+                void  PenaltyCreatePVectorMelting(Vec pg, void* inputs, double kmax,int analysis_type,int sub_analysis_type);
+                void  PenaltyConstrain(int* punstable,void* inputs,int analysis_type,int sub_analysis_type);
 };
 

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

@@ -260 +260 @@
         return new Penpair(*this); 
 }
-
+                
+int   Penpair::GetNumDofs(){
+        return numdofs;
+}
+

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

@@ -58 +58 @@
                 void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
                 void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
+                int   GetNumDofs();
                 Object* copy();
 };

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

@@ -313 +313 @@
                 CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type,sub_analysis_type);
         }
+        else if (analysis_type==ThermalAnalysisEnum()){
+                
+                CreateKMatrixThermal( Kgg,inputs,analysis_type,sub_analysis_type);
+        }
+        else if (analysis_type==MeltingAnalysisEnum()){
+                
+                CreateKMatrixMelting( Kgg,inputs,analysis_type,sub_analysis_type);
+        }
         else{
                 throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
@@ -834 +842 @@
 
                         /*Compute thickness: */
-                        GetParameterValue(&thickness, &h[0],gauss_coord);
 
                         /*Compute strain rate: */
@@ -1014 +1021 @@
                 CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type);
         }
+        else if (analysis_type==ThermalAnalysisEnum()){
+                
+                CreatePVectorThermal( pg,inputs,analysis_type,sub_analysis_type);
+        }
+        else if (analysis_type==MeltingAnalysisEnum()){
+                
+                CreatePVectorMelting( pg,inputs,analysis_type,sub_analysis_type);
+        }
         else{
                 throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
@@ -1208 +1223 @@
         double   tria_melting[3];
         double   tria_accumulation[3];
+        double   tria_geothermalflux[3];
         
         /*configuration: */
@@ -1238 +1254 @@
         tria_accumulation[2]=accumulation[g2];
 
+        tria_geothermalflux[0]=geothermalflux[g0];
+        tria_geothermalflux[1]=geothermalflux[g1];
+        tria_geothermalflux[2]=geothermalflux[g2];
+
         tria_node_ids[0]=node_ids[g0];
         tria_node_ids[1]=node_ids[g1];
@@ -1250 +1270 @@
         tria_node_offsets[2]=node_offsets[g2];
 
-        tria= new Tria(id,mid,mparid,tria_node_ids,tria_h,tria_s,tria_b,tria_k, tria_melting, tria_accumulation, friction_type,p,q,shelf,meanvel,epsvel,viscosity_overshoot);
+        tria= new Tria(id,mid,mparid,tria_node_ids,tria_h,tria_s,tria_b,tria_k, tria_melting, tria_accumulation, tria_geothermalflux,friction_type,p,q,shelf,meanvel,epsvel,viscosity_overshoot);
 
         tria->NodeConfiguration(tria_node_ids,tria_nodes,tria_node_offsets);
@@ -3128 +3148 @@
 
 }
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Penta::CreateKMatrixThermal"
+void  Penta::CreateKMatrixThermal(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
+
+        /* local declarations */
+        int i,j;
+        int found=0;
+
+        /* node data: */
+        const int    numgrids=6;
+        const int    NDOF1=1;
+        const int    numdof=NDOF1*numgrids;
+        double       xyz_list[numgrids][3];
+        int          doflist[numdof];
+        int          numberofdofspernode;
+
+        /* gaussian points: */
+        int     num_area_gauss,igarea,igvert;
+        double* first_gauss_area_coord  =  NULL;
+        double* second_gauss_area_coord =  NULL;
+        double* third_gauss_area_coord  =  NULL;
+        double* vert_gauss_coord = NULL;
+        double* area_gauss_weights  =  NULL;
+        double* vert_gauss_weights  =  NULL;
+        double  gauss_weight,area_gauss_weight,vert_gauss_weight;
+        double  gauss_coord[4];
+
+        int area_order=5;
+        int     num_vert_gauss=5;
+        
+        int            dofs[3]={0,1,2};
+        double         dt;
+        int            dt_exists;
+        
+        double         vxvyvz_list[numgrids][3];
+        double         vx_list[numgrids];
+        int            vx_list_exists;
+        double         u;
+        double         vy_list[numgrids];
+        int            vy_list_exists;
+        double         v;
+        double         vz_list[numgrids];
+        int            vz_list_exists;
+        double         w;
+
+        /* element parameters: */
+        int         onbed;
+        int         shelf;
+        int     steady_state;
+
+        /*matrices: */
+        double     K_terms[numdof][numdof]={0.0};
+        double     Ke_gaussian_conduct[numdof][numdof];
+        double     Ke_gaussian_advec[numdof][numdof];
+        double     Ke_gaussian_transient[numdof][numdof];
+        double     B[3][numdof];
+        double     Bprime[3][numdof];
+        double     B_conduct[3][numdof];
+        double     B_advec[3][numdof];
+        double     Bprime_advec[3][numdof];
+        double     L[numdof];
+        double     D_scalar;
+        double     D[3][3];
+        double     l1l2l3[3];
+        double     tl1l2l3D[3];
+        double     tBD[3][numdof];
+        double     tBD_conduct[3][numdof];
+        double     tBD_advec[3][numdof];
+        double     tLD[numdof];
+
+        double     Jdet;
+
+        /*Material properties: */
+        double         gravity,rho_ice,rho_water;
+        double         heatcapacity,thermalconductivity;
+        double         mixed_layer_capacity,thermal_exchange_velocity;
+
+        /*Collapsed formulation: */
+        Tria*  tria=NULL;
+        ParameterInputs* inputs=NULL;
+
+        /*recover pointers: */
+        inputs=(ParameterInputs*)vinputs;
+
+        /* Get node coordinates and dof list: */
+        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+        GetDofList(&doflist[0],&numberofdofspernode);
+        
+        // /*recovre material parameters: */
+        rho_water=matpar->GetRhoWater();
+        rho_ice=matpar->GetRhoIce();
+        gravity=matpar->GetG();
+        heatcapacity=matpar->GetHeatCapacity();
+        thermalconductivity=matpar->GetThermalConductivity();
+
+                
+        /*recover extra inputs from users, dt and velocity: */
+        found=inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes);
+        if(!found)throw ErrorException(__FUNCT__," could not find velocity in inputs!");
+        found=inputs->Recover("dt",&dt);
+        if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!");
+
+        for(i=0;i<numgrids;i++){
+                vx_list[i]=vxvyvz_list[i][0];
+                vy_list[i]=vxvyvz_list[i][1];
+                vz_list[i]=vxvyvz_list[i][2];
+        }
+
+        /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
+           get tria gaussian points as well as segment gaussian points. For tria gaussian 
+           points, order of integration is 2, because we need to integrate the product tB*D*B' 
+           which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian 
+           points, same deal, which yields 3 gaussian points.: */
+        
+        /*Get gaussian points: */
+        area_order=2;
+        num_vert_gauss=2;
+
+        GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss);
+
+        /* Start  looping on the number of gaussian points: */
+        for (igarea=0; igarea<num_area_gauss; igarea++){
+                for (igvert=0; igvert<num_vert_gauss; igvert++){
+                        /*Pick up the gaussian point: */
+                        area_gauss_weight=*(area_gauss_weights+igarea);
+                        vert_gauss_weight=*(vert_gauss_weights+igvert);
+                        gauss_weight=area_gauss_weight*vert_gauss_weight;
+                        gauss_coord[0]=*(first_gauss_area_coord+igarea); 
+                        gauss_coord[1]=*(second_gauss_area_coord+igarea);
+                        gauss_coord[2]=*(third_gauss_area_coord+igarea);
+                        gauss_coord[3]=*(vert_gauss_coord+igvert);
+
+                        /* Get Jacobian determinant: */
+                        GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
+
+                        /*Conduction: */
+
+                        /*Get B_conduct matrix: */
+                        GetB_conduct(&B_conduct[0][0],&xyz_list[0][0],gauss_coord); 
+
+                        /*Build D: */
+                        D_scalar=gauss_weight*Jdet*(thermalconductivity/(rho_ice*heatcapacity));
+
+                        if(sub_analysis_type!=SteadyAnalysisEnum()){
+                                D_scalar=D_scalar*dt;
+                        }
+
+                        D[0][0]=D_scalar; D[0][1]=0; D[0][2]=0;
+                        D[1][0]=0; D[1][1]=D_scalar; D[1][2]=0;
+                        D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar;
+
+
+                        /*  Do the triple product B'*D*B: */
+                        MatrixMultiply(&B_conduct[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_conduct[0][0],0);
+                        MatrixMultiply(&tBD_conduct[0][0],numdof,3,0,&B_conduct[0][0],3,numdof,0,&Ke_gaussian_conduct[0][0],0);
+
+                        /*Advection: */
+
+                        /*Get B_advec and Bprime_advec matrices: */
+                        GetB_advec(&B_advec[0][0],&xyz_list[0][0],gauss_coord); 
+                        GetBprime_advec(&Bprime_advec[0][0],&xyz_list[0][0],gauss_coord); 
+
+
+                        //Build the D matrix
+                        GetParameterValue(&u, &vx_list[0],gauss_coord);
+                        GetParameterValue(&v, &vy_list[0],gauss_coord);
+                        GetParameterValue(&w, &vz_list[0],gauss_coord);
+
+                        D_scalar=gauss_weight*Jdet;
+
+                        if(sub_analysis_type!=SteadyAnalysisEnum()){
+                                D_scalar=D_scalar*dt;
+                        }
+
+                        D[0][0]=D_scalar*u;D[0][1]=0;         D[0][2]=0;
+                        D[1][0]=0;         D[1][1]=D_scalar*v;D[1][2]=0;
+                        D[2][0]=0;         D[2][1]=0;         D[2][2]=D_scalar*w;
+
+                        /*  Do the triple product B'*D*Bprime: */
+                        MatrixMultiply(&B_advec[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_advec[0][0],0);
+                        MatrixMultiply(&tBD_advec[0][0],numdof,3,0,&Bprime_advec[0][0],3,numdof,0,&Ke_gaussian_advec[0][0],0);
+
+
+                        /*Transient: */
+                        if(sub_analysis_type!=SteadyAnalysisEnum()){
+                                GetNodalFunctions(&L[0], gauss_coord);
+                                D_scalar=gauss_weight*Jdet;
+                                D_scalar=D_scalar;
+
+                                /*  Do the triple product L'*D*L: */
+                                MatrixMultiply(&L[0],numdof,1,0,&D_scalar,1,1,0,&tLD[0],0);
+                                MatrixMultiply(&tLD[0],numdof,1,0,&L[0],1,numdof,0,&Ke_gaussian_transient[0][0],0);
+                        }
+                        else{
+                                for(i=0;i<numdof;i++){
+                                        for(j=0;j<numdof;j++){
+                                                Ke_gaussian_transient[i][j]=0;
+                                        }
+                                }
+                        }
+                                        
+                        /*Add Ke_gaussian to pKe: */
+                        for(i=0;i<numdof;i++){
+                                for(j=0;j<numdof;j++){
+                                        K_terms[i][j]+=Ke_gaussian_conduct[i][j]+Ke_gaussian_advec[i][j]+Ke_gaussian_transient[i][j];
+                                }
+                        }
+                        
+                }
+        }
+        
+        //Ice/ocean heat exchange flux on ice shelf base 
+
+        if(onbed && shelf){
+
+                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
+                tria->CreateKMatrixThermal(Kgg,inputs, analysis_type,sub_analysis_type);
+                delete tria;
+                return;
+        }
+
+        /*Add Ke_gg to global matrix Kgg: */
+        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES);
+
+        cleanup_and_return: 
+        xfree((void**)&first_gauss_area_coord);
+        xfree((void**)&second_gauss_area_coord);
+        xfree((void**)&third_gauss_area_coord);
+        xfree((void**)&area_gauss_weights);
+        xfree((void**)&vert_gauss_weights);
+        xfree((void**)&vert_gauss_coord);
+
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Penta::GetB_conduct"
+void Penta::GetB_conduct(double* B_conduct, double* xyz_list, double* gauss_coord){
+
+        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. 
+         * For grid i, Bi' can be expressed in the basic coordinate system
+         * by: 
+         *       Bi_conduct_basic=[ dh/dx ]
+         *                       [ dh/dy ]
+         *                       [ dh/dz ]
+         * where h is the interpolation function for grid i.
+         *
+         * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
+         */
+        
+        int i;
+        int calculationdof=3;
+        int numgrids=6;
+        int DOFPERGRID=1;
+
+        /*Same thing in the basic coordinate system: */
+        double dh1dh6_basic[calculationdof][numgrids];
+
+        /*Get dh1dh2dh3 in basic coordinates system : */
+        GetNodalFunctionsDerivativesBasic(&dh1dh6_basic[0][0],xyz_list,gauss_coord);
+
+        /*Build B': */
+        for (i=0;i<numgrids;i++){
+                *(B_conduct+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6_basic[0][i]; 
+                *(B_conduct+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6_basic[1][i]; 
+                *(B_conduct+DOFPERGRID*numgrids*2+DOFPERGRID*i)=dh1dh6_basic[2][i]; 
+        }
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Penta::GetB_advec"
+void Penta::GetB_advec(double* B_advec, double* xyz_list, double* gauss_coord){
+
+        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. 
+         * For grid i, Bi' can be expressed in the basic coordinate system
+         * by: 
+         *       Bi_conduct_basic=[ h ]
+         *                       [ h ]
+         *                       [ h ]
+         * where h is the interpolation function for grid i.
+         *
+         * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
+         */
+        
+        int i;
+        int calculationdof=3;
+        int numgrids=6;
+        int DOFPERGRID=1;
+
+        /*Same thing in the basic coordinate system: */
+        double l1l6[6];
+
+        /*Get dh1dh2dh3 in basic coordinates system : */
+        GetNodalFunctions(l1l6, gauss_coord);
+
+        /*Build B': */
+        for (i=0;i<numgrids;i++){
+                *(B_advec+DOFPERGRID*numgrids*0+DOFPERGRID*i)=l1l6[i]; 
+                *(B_advec+DOFPERGRID*numgrids*1+DOFPERGRID*i)=l1l6[i]; 
+                *(B_advec+DOFPERGRID*numgrids*2+DOFPERGRID*i)=l1l6[i]; 
+        }
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Penta::GetBprime_advec"
+void Penta::GetBprime_advec(double* Bprime_advec, double* xyz_list, double* gauss_coord){
+
+
+        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. 
+         * For grid i, Bi' can be expressed in the basic coordinate system
+         * by: 
+         *       Biprime_advec=[ dh/dx ]
+         *                       [ dh/dy ]
+         *                       [ dh/dz ]
+         * where h is the interpolation function for grid i.
+         *
+         * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
+         */
+        
+        int i;
+        int calculationdof=3;
+        int numgrids=6;
+        int DOFPERGRID=1;
+
+        /*Same thing in the basic coordinate system: */
+        double dh1dh6_basic[calculationdof][numgrids];
+
+        /*Get dh1dh2dh3 in basic coordinates system : */
+        GetNodalFunctionsDerivativesBasic(&dh1dh6_basic[0][0],xyz_list,gauss_coord);
+
+        /*Build B': */
+        for (i=0;i<numgrids;i++){
+                *(Bprime_advec+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6_basic[0][i]; 
+                *(Bprime_advec+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6_basic[1][i]; 
+                *(Bprime_advec+DOFPERGRID*numgrids*2+DOFPERGRID*i)=dh1dh6_basic[2][i]; 
+        }
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Penta::CreateKMatrixMelting"
+void  Penta::CreateKMatrixMelting(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
+        
+        Tria* tria=NULL;
+        if (!onbed){
+                return;
+        }
+        else{
+                
+                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
+                tria->CreateKMatrixMelting(Kgg,inputs, analysis_type,sub_analysis_type);
+                delete tria;
+                return;
+        }
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Penta::CreatePVectorThermal"
+void Penta::CreatePVectorThermal( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type){
+
+
+        /*indexing: */
+        int i,j;
+        int found=0;
+
+        const int  numgrids=6;
+        const int  NDOF1=1;
+        const int  numdof=numgrids*NDOF1;
+        int        doflist[numdof];
+        int        numberofdofspernode;
+
+        /*Grid data: */
+        double        xyz_list[numgrids][3];
+        
+        /* gaussian points: */
+        int     num_area_gauss,igarea,igvert;
+        double* first_gauss_area_coord  =  NULL;
+        double* second_gauss_area_coord =  NULL;
+        double* third_gauss_area_coord  =  NULL;
+        double* vert_gauss_coord = NULL;
+        double* area_gauss_weights  =  NULL;
+        double* vert_gauss_weights  =  NULL;
+        double  gauss_weight,area_gauss_weight,vert_gauss_weight;
+        double  gauss_coord[4];
+        int area_order=2;
+        int     num_vert_gauss=3;
+        
+        double         dt;
+        double         vx_list[numgrids];
+        double         vy_list[numgrids];
+        double         vz_list[numgrids];
+        double         vxvyvz_list[numgrids][3];
+        double         pressure_list[3];
+        double         pressure;
+        double         temperature_list[numgrids];
+        double         temperature;
+
+        /*Material properties: */
+        double         gravity,rho_ice,rho_water;
+        double         mixed_layer_capacity,heatcapacity;
+        double         beta,meltingpoint,thermal_exchange_velocity;
+
+        /* element parameters: */
+        int     friction_type;
+        
+        int            dofs[3]={0,1,2};
+        int            dofs1[1]={0};
+
+        /*matrices: */
+        double     P_terms[numdof]={0.0};
+        double     L[numdof];
+        double     l1l2l3[3];
+        double     alpha2_list[3];
+        double     basalfriction_list[3]={0.0};
+        double     basalfriction;
+        double     epsilon[6];
+        double     epsilon_sqr[3][3];
+        double     epsilon_matrix[3][3];
+
+        double     Jdet;
+        double     viscosity;
+        double     epsilon_eff;
+        double     phi;
+        double     t_pmp;
+        double     scalar;
+        double     scalar_def;
+        double     scalar_ocean;
+        double     scalar_transient;
+
+        /*Collapsed formulation: */
+        Tria*  tria=NULL;
+        ParameterInputs* inputs=NULL;
+
+        /*recover pointers: */
+        inputs=(ParameterInputs*)vinputs;
+
+        /* Get node coordinates and dof list: */
+        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+        GetDofList(&doflist[0],&numberofdofspernode);
+        
+        // /*recovre material parameters: */
+        rho_water=matpar->GetRhoWater();
+        rho_ice=matpar->GetRhoIce();
+        gravity=matpar->GetG();
+        heatcapacity=matpar->GetHeatCapacity();
+        beta=matpar->GetBeta();
+        meltingpoint=matpar->GetMeltingPoint();
+
+        /*recover extra inputs from users, dt and velocity: */
+        found=inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes);
+        if(!found)throw ErrorException(__FUNCT__," could not find velocity in inputs!");
+
+        found=inputs->Recover("dt",&dt);
+        if((!found) && (sub_analysis_type==TransientAnalysisEnum()))throw ErrorException(__FUNCT__," could not find dt in inputs!");
+        
+        found=inputs->Recover("pressure",&pressure_list[0],1,dofs1,numgrids,(void**)nodes);
+        if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
+
+        found=inputs->Recover("temperature",&temperature_list[0],1,dofs1,numgrids,(void**)nodes);
+        if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!");
+
+        for(i=0;i<numgrids;i++){
+                vx_list[i]=vxvyvz_list[i][0];
+                vy_list[i]=vxvyvz_list[i][1];
+                vz_list[i]=vxvyvz_list[i][2];
+        }       
+
+        /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
+           get tria gaussian points as well as segment gaussian points. For tria gaussian 
+           points, order of integration is 2, because we need to integrate the product tB*D*B' 
+           which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian 
+           points, same deal, which yields 3 gaussian points.: */
+        
+        /*Get gaussian points: */
+        GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss);
+
+        /* Start  looping on the number of gaussian points: */
+        for (igarea=0; igarea<num_area_gauss; igarea++){
+                for (igvert=0; igvert<num_vert_gauss; igvert++){
+                        /*Pick up the gaussian point: */
+                        area_gauss_weight=*(area_gauss_weights+igarea);
+                        vert_gauss_weight=*(vert_gauss_weights+igvert);
+                        gauss_weight=area_gauss_weight*vert_gauss_weight;
+                        gauss_coord[0]=*(first_gauss_area_coord+igarea); 
+                        gauss_coord[1]=*(second_gauss_area_coord+igarea);
+                        gauss_coord[2]=*(third_gauss_area_coord+igarea);
+                        gauss_coord[3]=*(vert_gauss_coord+igvert);
+        
+                        /*Compute strain rate and viscosity: */
+                        GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
+                        matice->GetViscosity3d(&viscosity,&epsilon[0]);
+
+                        /* Get Jacobian determinant: */
+                        GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
+
+                        /* Get nodal functions */
+                        GetNodalFunctions(&L[0], gauss_coord);
+
+                        /*Build deformational heating: */
+                        GetPhi(&phi, &epsilon[0], viscosity);
+
+                        /*Build pe_gaussian */
+                        if(sub_analysis_type==SteadyAnalysisEnum()){
+                                scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss_weight;
+                        }
+                        else{
+                                scalar_def=dt*phi/(rho_ice*heatcapacity)*Jdet*gauss_weight;
+                        }
+
+                        for(i=0;i<numgrids;i++){
+                                P_terms[i]+=scalar_def*L[i];
+                        }
+
+                        /* Build transient now */
+                        if(sub_analysis_type==TransientAnalysisEnum()){
+                                scalar_transient=temperature*Jdet*gauss_weight;
+                                for(i=0;i<numgrids;i++){
+                                        P_terms[i]+=scalar_transient*L[i];
+                                }
+                        }
+                }
+        }
+
+        /*Add pe_g to global vector pg: */
+        VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
+
+        /* Ice/ocean heat exchange flux on ice shelf base */
+        if(onbed && shelf){
+
+                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
+                tria->CreatePVectorThermalShelf(pg,inputs, analysis_type,sub_analysis_type);
+                delete tria;
+                return;
+        }
+
+        /* Geothermal flux on ice sheet base and basal friction */
+        if(onbed && !shelf){
+                
+                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
+                tria->CreatePVectorThermalSheet(pg,inputs, analysis_type,sub_analysis_type);
+                delete tria;
+                return;
+        }
+
+        cleanup_and_return:
+        xfree((void**)first_gauss_area_coord);
+        xfree((void**)second_gauss_area_coord);
+        xfree((void**)third_gauss_area_coord);
+        xfree((void**)vert_gauss_coord);
+        xfree((void**)area_gauss_weights);
+        xfree((void**)vert_gauss_weights);
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Penta::CreatePVectorMelting"
+void Penta::CreatePVectorMelting( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type){
+        return;
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Penta::GetPhi"
+void Penta::GetPhi(double* phi, double*  epsilon, double viscosity){
+        /*Compute deformational heating from epsilon and viscosity */
+
+        double epsilon_matrix[3][3];
+        double epsilon_eff;
+        double epsilon_sqr[3][3];
+
+        /* Build epsilon matrix */
+        epsilon_matrix[0][0]=*(epsilon+0);
+        epsilon_matrix[1][0]=*(epsilon+3);
+        epsilon_matrix[2][0]=*(epsilon+4);
+        epsilon_matrix[0][1]=*(epsilon+3);
+        epsilon_matrix[1][1]=*(epsilon+1);
+        epsilon_matrix[2][1]=*(epsilon+5);
+        epsilon_matrix[0][2]=*(epsilon+4);
+        epsilon_matrix[1][2]=*(epsilon+5);
+        epsilon_matrix[2][2]=*(epsilon+2);
+
+        /* Effective value of epsilon_matrix */
+        epsilon_sqr[0][0]=pow(epsilon_matrix[0][0],2);
+        epsilon_sqr[1][0]=pow(epsilon_matrix[1][0],2);
+        epsilon_sqr[2][0]=pow(epsilon_matrix[2][0],2);
+        epsilon_sqr[0][1]=pow(epsilon_matrix[0][1],2);
+        epsilon_sqr[1][1]=pow(epsilon_matrix[1][1],2);
+        epsilon_sqr[2][1]=pow(epsilon_matrix[2][1],2);
+        epsilon_sqr[0][2]=pow(epsilon_matrix[0][2],2);
+        epsilon_sqr[1][2]=pow(epsilon_matrix[1][2],2);
+        epsilon_sqr[2][2]=pow(epsilon_matrix[2][2],2);
+
+        epsilon_eff=1/pow(2,1.0/2.0)*pow((epsilon_sqr[0][0]+epsilon_sqr[0][1]+ epsilon_sqr[0][2]+ epsilon_sqr[1][0]+ epsilon_sqr[1][1]+ epsilon_sqr[1][2]+ epsilon_sqr[2][0]+ epsilon_sqr[2][1]+ epsilon_sqr[2][2]),1.0/2.0);
+        *phi=2*pow(epsilon_eff,2.0)*viscosity;
+}

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

@@ -130 +130 @@
                 void  ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp);
                 void  GetNodalFunctionsStokes(double* l1l7, double* gauss_coord);
+                void  CreateKMatrixThermal(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
+                void  GetB_conduct(double* B_conduct, double* xyz_list, double* gauss_coord);
+                void  GetB_advec(double* B_advec, double* xyz_list, double* gauss_coord);
+                void  GetBprime_advec(double* Bprime_advec, double* xyz_list, double* gauss_coord);
+                void  CreateKMatrixMelting(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
+                void  CreatePVectorThermal( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type);
+                void  CreatePVectorMelting( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type);
+                void  GetPhi(double* phi, double*  epsilon, double viscosity);
 
 

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

@@ -33 +33 @@
 
 Tria::Tria(int tria_id,int tria_mid,int tria_mparid,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,double tria_meanvel,double tria_epsvel,
+                                double tria_accumulation[3],double tria_geothermalflux[3],int tria_friction_type,double tria_p,double tria_q,int tria_shelf,double tria_meanvel,double tria_epsvel,
                                 double tria_viscosity_overshoot){
         
@@ -51 +51 @@
                 melting[i]=tria_melting[i];
                 accumulation[i]=tria_accumulation[i];
-        }
+                geothermalflux[i]=tria_geothermalflux[i]; }
         matice=NULL;
         matice_offset=UNDEF;
@@ -87 +87 @@
         printf("   melting=[%g,%g,%g]\n",melting[0],melting[1],melting[2]);
         printf("   accumulation=[%g,%g,%g]\n",accumulation[0],accumulation[1],accumulation[2]);
+        printf("   geothermalflux=[%g,%g,%g]\n",geothermalflux[0],geothermalflux[1],geothermalflux[2]);
         printf("   friction_type: %i\n",friction_type);
         printf("   p: %g\n",p);
@@ -136 +137 @@
         memcpy(marshalled_dataset,&melting,sizeof(melting));marshalled_dataset+=sizeof(melting);
         memcpy(marshalled_dataset,&accumulation,sizeof(accumulation));marshalled_dataset+=sizeof(accumulation);
+        memcpy(marshalled_dataset,&geothermalflux,sizeof(geothermalflux));marshalled_dataset+=sizeof(geothermalflux);
         memcpy(marshalled_dataset,&friction_type,sizeof(friction_type));marshalled_dataset+=sizeof(friction_type);
         memcpy(marshalled_dataset,&onbed,sizeof(onbed));marshalled_dataset+=sizeof(onbed);
@@ -166 +168 @@
                 +sizeof(melting)
                 +sizeof(accumulation)
+                +sizeof(geothermalflux)
                 +sizeof(friction_type)
                 +sizeof(onbed)
@@ -208 +211 @@
         memcpy(&melting,marshalled_dataset,sizeof(melting));marshalled_dataset+=sizeof(melting);
         memcpy(&accumulation,marshalled_dataset,sizeof(accumulation));marshalled_dataset+=sizeof(accumulation);
+        memcpy(&geothermalflux,marshalled_dataset,sizeof(geothermalflux));marshalled_dataset+=sizeof(geothermalflux);
         memcpy(&friction_type,marshalled_dataset,sizeof(friction_type));marshalled_dataset+=sizeof(friction_type);
         memcpy(&onbed,marshalled_dataset,sizeof(onbed));marshalled_dataset+=sizeof(onbed);
@@ -264 +268 @@
 
 #undef __FUNCT__ 
-#define __FUNCT__ "Tria::CreateKMatrix"
+#define __FUNCT__ "Tria::reateKMatrix"
 
 void  Tria::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
@@ -305 +309 @@
         /* node data: */
         const int    numgrids=3;
-        const int    numdofs=2*numgrids;
+        const int    numdof=2*numgrids;
         double       xyz_list[numgrids][3];
-        int          doflist[numdofs];
+        int          doflist[numdof];
         int          numberofdofspernode;
         
@@ -329 +333 @@
 
         /* matrices: */
-        double B[3][numdofs];
-        double Bprime[3][numdofs];
+        double B[3][numdof];
+        double Bprime[3][numdof];
         double D[3][3]={{ 0,0,0 },{0,0,0},{0,0,0}};              // material matrix, simple scalar matrix.
         double D_scalar;
 
         /* local element matrices: */
-        double Ke_gg[numdofs][numdofs]; //local element stiffness matrix 
-        double Ke_gg_gaussian[numdofs][numdofs]; //stiffness matrix evaluated at the gaussian point.
+        double Ke_gg[numdof][numdof]; //local element stiffness matrix 
+        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
         
         double Jdet;
@@ -360 +364 @@
 
         /* Set Ke_gg to 0: */
-        for(i=0;i<numdofs;i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]=0.0;
+        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
 
         #ifdef _DEBUGELEMENTS_
@@ -447 +451 @@
 
                 /*  Do the triple product tB*D*Bprime: */
-                TripleMultiply( &B[0][0],3,numdofs,1,
+                TripleMultiply( &B[0][0],3,numdof,1,
                                           &D[0][0],3,3,0,
-                                          &Bprime[0][0],3,numdofs,0,
+                                          &Bprime[0][0],3,numdof,0,
                                           &Ke_gg_gaussian[0][0],0);
 
                 /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
-                for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
+                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
                 
                 #ifdef _DEBUGELEMENTS_
@@ -459 +463 @@
                         printf("      B:\n");
                         for(i=0;i<3;i++){
-                                for(j=0;j<numdofs;j++){
+                                for(j=0;j<numdof;j++){
                                         printf("%g ",B[i][j]);
                                 }
@@ -466 +470 @@
                         printf("      Bprime:\n");
                         for(i=0;i<3;i++){
-                                for(j=0;j<numdofs;j++){
+                                for(j=0;j<numdof;j++){
                                         printf("%g ",Bprime[i][j]);
                                 }
@@ -476 +480 @@
 
         /*Add Ke_gg to global matrix Kgg: */
-        MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
+        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
 
 
@@ -488 +492 @@
         if(my_rank==RANK && id==ELID){ 
                 printf("      Ke_gg erms:\n");
-                for( i=0; i<numdofs; i++){
-                        for (j=0;j<numdofs;j++){
+                for( i=0; i<numdof; i++){
+                        for (j=0;j<numdof;j++){
                                 printf("%g ",Ke_gg[i][j]);
                         }
@@ -495 +499 @@
                 }
                 printf("      Ke_gg row_indices:\n");
-                for( i=0; i<numdofs; i++){
+                for( i=0; i<numdof; i++){
                         printf("%i ",doflist[i]);
                 }
@@ -520 +524 @@
         /* node data: */
         const int    numgrids=3;
-        const int    numdofs=2*numgrids;
+        const int    numdof=2*numgrids;
         double       xyz_list[numgrids][3];
-        int          doflist[numdofs];
+        int          doflist[numdof];
         int          numberofdofspernode;
         
@@ -535 +539 @@
 
         /* matrices: */
-        double L[2][numdofs];
+        double L[2][numdof];
         double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag
         double DL_scalar;
 
         /* local element matrices: */
-        double Ke_gg[numdofs][numdofs]; //local element stiffness matrix 
-        double Ke_gg_gaussian[numdofs][numdofs]; //stiffness matrix contribution from drag
+        double Ke_gg[numdof][numdof]; //local element stiffness matrix 
+        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
         
         double Jdet;
@@ -570 +574 @@
 
         /* Set Ke_gg to 0: */
-        for(i=0;i<numdofs;i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]=0.0;
+        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
 
         if (shelf){
@@ -658 +662 @@
                 
                 /*  Do the triple producte tL*D*L: */
-                TripleMultiply( &L[0][0],2,numdofs,1,
+                TripleMultiply( &L[0][0],2,numdof,1,
                                         &DL[0][0],2,2,0,
-                                        &L[0][0],2,numdofs,0,
+                                        &L[0][0],2,numdof,0,
                                         &Ke_gg_gaussian[0][0],0);
 
-                for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
+                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
 
         } // for (ig=0; ig<num_gauss; ig++)
 
         /*Add Ke_gg to global matrix Kgg: */
-        MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
+        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
 
         cleanup_and_return: 
@@ -689 +693 @@
         const int    numgrids=3;
         const int    NDOF1=1;
-        const int    numdofs=NDOF1*numgrids;
+        const int    numdof=NDOF1*numgrids;
         double       xyz_list[numgrids][3];
-        int          doflist[numdofs];
+        int          doflist[numdof];
         int          numberofdofspernode;
         
@@ -708 +712 @@
 
         /* local element matrices: */
-        double Ke_gg[numdofs][numdofs]; //local element stiffness matrix 
-        double Ke_gg_gaussian[numdofs][numdofs]; //stiffness matrix evaluated at the gaussian point.
+        double Ke_gg[numdof][numdof]; //local element stiffness matrix 
+        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
         
         double Jdet;
@@ -718 +722 @@
 
         /* Set Ke_gg to 0: */
-        for(i=0;i<numdofs;i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]=0.0;
+        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
 
         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
@@ -747 +751 @@
 
                 /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
-                for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
+                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
         } //for (ig=0; ig<num_gauss; ig++
 
         /*Add Ke_gg to global matrix Kgg: */
-        MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
+        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
                 
         cleanup_and_return: 
@@ -790 +794 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Tria::CreatePVectorDiagnosticHoriz"
-
 void Tria::CreatePVectorDiagnosticHoriz( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){
 
@@ -797 +800 @@
         /* node data: */
         const int    numgrids=3;
-        const int    numdofs=2*numgrids;
+        const int    numdof=2*numgrids;
         const int    NDOF2=2;
         double       xyz_list[numgrids][3];
-        int          doflist[numdofs];
+        int          doflist[numdof];
         int          numberofdofspernode;
         
@@ -824 +827 @@
 
         /*element vector at the gaussian points: */
-        double  pe_g[numdofs];
-        double  pe_g_gaussian[numdofs];
+        double  pe_g[numdof];
+        double  pe_g_gaussian[numdof];
 
         /*input parameters for structural analysis (diagnostic): */
@@ -840 +843 @@
 
         /* Set pe_g to 0: */
-        for(i=0;i<numdofs;i++) pe_g[i]=0.0;
+        for(i=0;i<numdof;i++) pe_g[i]=0.0;
 
 
@@ -927 +930 @@
 
                 /*Add pe_g_gaussian vector to pe_g: */
-                for( i=0; i<numdofs; i++)pe_g[i]+=pe_g_gaussian[i];
+                for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i];
 
         } //for (ig=0; ig<num_gauss; ig++)
@@ -946 +949 @@
 
         /*Add pe_g to global vector pg: */
-        VecSetValues(pg,numdofs,doflist,(const double*)pe_g,ADD_VALUES);
+        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
 
         cleanup_and_return: 
@@ -966 +969 @@
         const int    numgrids=3;
         const int    NDOF1=1;
-        const int    numdofs=NDOF1*numgrids;
+        const int    numdof=NDOF1*numgrids;
         double       xyz_list[numgrids][3];
-        int          doflist[numdofs];
+        int          doflist[numdof];
         int          numberofdofspernode;
         
@@ -987 +990 @@
 
         /*element vector at the gaussian points: */
-        double  pe_g[numdofs];
-        double  pe_g_gaussian[numdofs];
+        double  pe_g[numdof];
+        double  pe_g_gaussian[numdof];
         double  param[3];
         double  slope[2];
@@ -997 +1000 @@
 
         /* Set pe_g to 0: */
-        for(i=0;i<numdofs;i++) pe_g[i]=0.0;
+        for(i=0;i<numdof;i++) pe_g[i]=0.0;
 
         if ( (sub_analysis_type==SurfaceXAnalysisEnum()) || (sub_analysis_type==SurfaceYAnalysisEnum())){
-                for(i=0;i<numdofs;i++) param[i]=s[i];
+                for(i=0;i<numdof;i++) param[i]=s[i];
         }
         if ( (sub_analysis_type==BedXAnalysisEnum()) || (sub_analysis_type==BedYAnalysisEnum())){
-                for(i=0;i<numdofs;i++) param[i]=b[i];
+                for(i=0;i<numdof;i++) param[i]=b[i];
         }
 
@@ -1028 +1031 @@
                 /*Build pe_g_gaussian vector: */
                 if ( (sub_analysis_type==SurfaceXAnalysisEnum()) || (sub_analysis_type==BedXAnalysisEnum())){
-                        for(i=0;i<numdofs;i++) pe_g_gaussian[i]=Jdet*gauss_weight*slope[0]*l1l2l3[i];
+                        for(i=0;i<numdof;i++) pe_g_gaussian[i]=Jdet*gauss_weight*slope[0]*l1l2l3[i];
                 }
                 if ( (sub_analysis_type==SurfaceYAnalysisEnum()) || (sub_analysis_type==BedYAnalysisEnum())){
-                        for(i=0;i<numdofs;i++) pe_g_gaussian[i]=Jdet*gauss_weight*slope[1]*l1l2l3[i];
+                        for(i=0;i<numdof;i++) pe_g_gaussian[i]=Jdet*gauss_weight*slope[1]*l1l2l3[i];
                 }
 
                 /*Add pe_g_gaussian vector to pe_g: */
-                for( i=0; i<numdofs; i++)pe_g[i]+=pe_g_gaussian[i];
+                for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i];
 
         } //for (ig=0; ig<num_gauss; ig++)
 
         /*Add pe_g to global vector pg: */
-        VecSetValues(pg,numdofs,doflist,(const double*)pe_g,ADD_VALUES);
+        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
 
         cleanup_and_return: 
@@ -1071 +1074 @@
         inputs->Recover("drag",&k[0],1,dofs,3,(void**)nodes);
         inputs->Recover("melting",&melting[0],1,dofs,3,(void**)nodes);
-        inputs->Recover("accumulation_param",&accumulation[0],1,dofs,3,(void**)nodes);
+        inputs->Recover("accumulation",&accumulation[0],1,dofs,3,(void**)nodes);
+        inputs->Recover("geothermalflux",&geothermalflux[0],1,dofs,3,(void**)nodes);
         
         //Update material if necessary
@@ -1556 +1560 @@
         /* node data: */
         const int    numgrids=3;
-        const int    numdofs=2*numgrids;
+        const int    numdof=2*numgrids;
         const int    NDOF2=2;
         double       xyz_list[numgrids][3];
-        int          doflist[numdofs];
+        int          doflist[numdof];
         int          numberofdofspernode;
 
@@ -1588 +1592 @@
 
         /*element vector : */
-        double  due_g[numdofs];
-        double  due_g_gaussian[numdofs];
+        double  due_g[numdof];
+        double  due_g_gaussian[numdof];
 
         /* Jacobian: */
@@ -1613 +1617 @@
 
         /* Set due_g to 0: */
-        for(i=0;i<numdofs;i++) due_g[i]=0.0;
+        for(i=0;i<numdof;i++) due_g[i]=0.0;
 
         /* Recover input data: */
@@ -1735 +1739 @@
                 
                 /*Add due_g_gaussian vector to due_g: */
-                for( i=0; i<numdofs; i++){
+                for( i=0; i<numdof; i++){
                         due_g[i]+=due_g_gaussian[i];
                 }
@@ -1741 +1745 @@
         
         /*Add due_g to global vector du_g: */
-        VecSetValues(du_g,numdofs,doflist,(const double*)due_g,ADD_VALUES);
+        VecSetValues(du_g,numdof,doflist,(const double*)due_g,ADD_VALUES);
         
         cleanup_and_return: 
@@ -1773 +1777 @@
         const int    numgrids=3;
         const int    NDOF2=2;
-        const int    numdofs=NDOF2*numgrids;
+        const int    numdof=NDOF2*numgrids;
         double       xyz_list[numgrids][3];
-        int          doflist[numdofs];
+        int          doflist[numdof];
         int          numberofdofspernode;
 
@@ -1814 +1818 @@
 
         /*element vector at the gaussian points: */
-        double  grade_g[numdofs];
+        double  grade_g[numdof];
         double  grade_g_gaussian[numgrids];
 
@@ -1836 +1840 @@
 
         /* Set grade_g to 0: */
-        for(i=0;i<numdofs;i++) grade_g[i]=0.0;
+        for(i=0;i<numdof;i++) grade_g[i]=0.0;
 
         /* recover input parameters: */
@@ -1947 +1951 @@
 
         /*Add grade_g to global vector grad_g: */
-        VecSetValues(grad_g,numdofs,doflist,(const double*)grade_g,ADD_VALUES);
+        VecSetValues(grad_g,numdof,doflist,(const double*)grade_g,ADD_VALUES);
         
         cleanup_and_return: 
@@ -1967 +1971 @@
         const int    NDOF1=1;
         const int    NDOF2=2;
-        const int    numdofs=NDOF2*numgrids;
+        const int    numdof=NDOF2*numgrids;
         double       xyz_list[numgrids][3];
-        int          doflist[numdofs];
+        int          doflist[numdof];
         int          numberofdofspernode;
 
@@ -1989 +1993 @@
 
         /*element vector at the gaussian points: */
-        double  grade_g[numdofs];
+        double  grade_g[numdof];
         double  grade_g_gaussian[numgrids];
 
@@ -2023 +2027 @@
 
         /* Set grade_g to 0: */
-        for(i=0;i<numdofs;i++) grade_g[i]=0.0;
+        for(i=0;i<numdof;i++) grade_g[i]=0.0;
 
         /* recover input parameters: */
@@ -2088 +2092 @@
         
         /*Add grade_g to global vector grad_g: */
-        VecSetValues(grad_g,numdofs,doflist,(const double*)grade_g,ADD_VALUES);
+        VecSetValues(grad_g,numdof,doflist,(const double*)grade_g,ADD_VALUES);
         
         cleanup_and_return: 
@@ -2109 +2113 @@
         /* node data: */
         const int    numgrids=3;
-        const int    numdofs=2*numgrids;
+        const int    numdof=2*numgrids;
         const int    NDOF2=2;
         double       xyz_list[numgrids][3];
-        int          doflist[numdofs];
+        int          doflist[numdof];
         int          numberofdofspernode;
 
@@ -2291 +2295 @@
         const int    numgrids=3;
         const int    NDOF1=1;
-        const int    numdofs=NDOF1*numgrids;
+        const int    numdof=NDOF1*numgrids;
         double       xyz_list[numgrids][3];
-        int          doflist[numdofs];
+        int          doflist[numdof];
         int          numberofdofspernode;
 
@@ -2322 +2326 @@
 
         /* local element matrices: */
-        double Ke_gg[numdofs][numdofs]; //local element stiffness matrix 
-        double Ke_gg_gaussian[numdofs][numdofs]; //stiffness matrix evaluated at the gaussian point.
+        double Ke_gg[numdof][numdof]; //local element stiffness matrix 
+        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
 
         ParameterInputs* inputs=NULL;
@@ -2335 +2339 @@
 
         /* Set Ke_gg to 0: */
-        for(i=0;i<numdofs;i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]=0.0;
+        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
 
         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
@@ -2394 +2398 @@
 
                 /* Add the Ke_gg_gaussian, onto Ke_gg: */
-                for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
+                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
 
                 
@@ -2400 +2404 @@
 
         /*Add Ke_gg to global matrix Kgg: */
-        MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
+        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
         
         cleanup_and_return: 
@@ -2419 +2423 @@
         const int    numgrids=3;
         const int    NDOF1=1;
-        const int    numdofs=NDOF1*numgrids;
+        const int    numdof=NDOF1*numgrids;
         double       xyz_list[numgrids][3];
-        int          doflist[numdofs];
+        int          doflist[numdof];
         int          numberofdofspernode;
         
@@ -2440 +2444 @@
 
         /*element vector at the gaussian points: */
-        double  pe_g[numdofs];
-        double  pe_g_gaussian[numdofs];
+        double  pe_g[numdof];
+        double  pe_g_gaussian[numdof];
 
         /* matrices: */
@@ -2471 +2475 @@
 
         /* Set pe_g to 0: */
-        for(i=0;i<numdofs;i++) pe_g[i]=0.0;
+        for(i=0;i<numdof;i++) pe_g[i]=0.0;
 
         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
@@ -2511 +2515 @@
         
                 /*Add pe_g_gaussian vector to pe_g: */
-                for( i=0; i<numdofs; i++)pe_g[i]+=pe_g_gaussian[i];
+                for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i];
         
         }
 
         /*Add pe_g to global vector pg: */
-        VecSetValues(pg,numdofs,doflist,(const double*)pe_g,ADD_VALUES);
+        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
 
         cleanup_and_return: 
@@ -2550 +2554 @@
 
 }
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::CreateKMatrixThermal"
+void  Tria::CreateKMatrixThermal(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
+
+        int i,j;
+        int found=0;
+        
+        /* node data: */
+        const int    numgrids=3;
+        const int    NDOF1=1;
+        const int    numdof=NDOF1*numgrids;
+        double       xyz_list[numgrids][3];
+        int          doflist[numdof];
+        int          numberofdofspernode;
+
+        double mixed_layer_capacity;
+        double thermal_exchange_velocity;
+        double rho_water;
+        double rho_ice;
+        double heatcapacity;
+        double dt;
+
+        int     num_gauss,ig;
+        double* first_gauss_area_coord  =  NULL;
+        double* second_gauss_area_coord =  NULL;
+        double* third_gauss_area_coord  =  NULL;
+        double* gauss_weights           =  NULL;
+        double  gauss_weight;
+        double  gauss_coord[3];
+
+        /*matrices: */
+        double  Jdet;
+        double  K_terms[numdof][numdof]={0.0};
+        double  Ke_gaussian[numdof][numdof]={0.0};
+        double  l1l2l3[numgrids];
+        double     tl1l2l3D[3];
+        double  D_scalar;
+        ParameterInputs* inputs=NULL;
+
+        /*recover pointers: */
+        inputs=(ParameterInputs*)vinputs;
+
+        /*recover extra inputs from users, dt: */
+        inputs->Recover("dt",&dt);
+        if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!");
+
+        /* Get node coordinates and dof list: */
+        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+        GetDofList(&doflist[0],&numberofdofspernode);
+
+        //recover material parameters
+        mixed_layer_capacity=matpar->GetMixedLayerCapacity();
+        thermal_exchange_velocity=matpar->GetThermalExchangeVelocity();
+        rho_water=matpar->GetRhoWater();
+        rho_ice=matpar->GetRhoIce();
+        heatcapacity=matpar->GetHeatCapacity();
+
+
+        GaussTria (&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+
+        /* Start looping on the number of gauss (nodes on the bedrock) */
+        for (ig=0; ig<num_gauss; ig++){
+                gauss_weight=*(gauss_weights+ig);
+                gauss_coord[0]=*(first_gauss_area_coord+ig); 
+                gauss_coord[1]=*(second_gauss_area_coord+ig);
+                gauss_coord[2]=*(third_gauss_area_coord+ig);
+                
+                //Get the Jacobian determinant
+                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss_coord);
+                
+                /*Get nodal functions values: */
+                GetNodalFunctions(&l1l2l3[0], gauss_coord);
+                                
+                /*Calculate DL on gauss point */
+                D_scalar=gauss_weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(heatcapacity*rho_ice);
+                if(sub_analysis_type!=SteadyAnalysisEnum()){
+                        D_scalar=dt*D_scalar;
+                }
+
+                /*  Do the triple product tL*D*L: */
+                MatrixMultiply(&l1l2l3[0],numdof,1,0,&D_scalar,1,1,0,&tl1l2l3D[0],0);
+                MatrixMultiply(&tl1l2l3D[0],numdof,1,0,&l1l2l3[0],1,numdof,0,&Ke_gaussian[0][0],0);
+
+                for(i=0;i<3;i++){
+                        for(j=0;j<3;j++){
+                                K_terms[i][j]+=Ke_gaussian[i][j];
+                        }
+                }
+        }
+        
+        /*Add Ke_gg to global matrix Kgg: */
+        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES);
+
+        cleanup_and_return:
+        xfree((void**)&first_gauss_area_coord);
+        xfree((void**)&second_gauss_area_coord);
+        xfree((void**)&third_gauss_area_coord);
+        xfree((void**)&gauss_weights);
+
+}
+
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::CreateKMatrixMelting"
+void  Tria::CreateKMatrixMelting(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
+        
+        /*indexing: */
+        int i,j;
+
+        const int  numgrids=3;
+        const int  NDOF1=1;
+        const int  numdof=numgrids*NDOF1;
+        int        doflist[numdof];
+        int        numberofdofspernode;
+
+        /*Grid data: */
+        double     xyz_list[numgrids][3];
+                
+        /*Material constants */
+        double     heatcapacity,latentheat;
+
+        /* gaussian points: */
+        int     num_area_gauss,ig;
+        double* gauss_weights  =  NULL;
+        double* first_gauss_area_coord  =  NULL;
+        double* second_gauss_area_coord =  NULL;
+        double* third_gauss_area_coord  =  NULL;
+        double  gauss_weight;
+        double  gauss_coord[3];
+
+        /*matrices: */
+        double     Jdet;
+        double     D_scalar;
+        double     K_terms[numdof][numdof]={0.0};
+        double     L[3];
+        double     tLD[3];
+        double     Ke_gaussian[numdof][numdof]={0.0};
+
+        /*Recover constants of ice */
+        latentheat=matpar->GetLatentHeat();
+        heatcapacity=matpar->GetHeatCapacity();
+
+        /* Get node coordinates and dof list: */
+        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+        GetDofList(&doflist[0],&numberofdofspernode);
+
+        /* Get gaussian points and weights: */
+        GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+
+        /* Start looping on the number of gauss  (nodes on the bedrock) */
+        for (ig=0; ig<num_area_gauss; ig++){
+                gauss_weight=*(gauss_weights+ig);
+                gauss_coord[0]=*(first_gauss_area_coord+ig); 
+                gauss_coord[1]=*(second_gauss_area_coord+ig);
+                gauss_coord[2]=*(third_gauss_area_coord+ig);
+                
+                //Get the Jacobian determinant
+                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss_coord);
+                
+                /*Get L matrix : */
+                GetL(&L[0], &xyz_list[0][0], gauss_coord,NDOF1);
+                                
+                /*Calculate DL on gauss point */
+                D_scalar=latentheat/heatcapacity*gauss_weight*Jdet;
+
+                /*  Do the triple product tL*D*L: */
+                MatrixMultiply(&L[0],numdof,1,0,&D_scalar,1,1,0,&tLD[0],0);
+                MatrixMultiply(&tLD[0],numdof,1,0,&L[0],1,numdof,0,&Ke_gaussian[0][0],0);
+
+                for(i=0;i<numgrids;i++){
+                        for(j=0;j<numgrids;j++){
+                        K_terms[i][j]+=Ke_gaussian[i][j];
+                        }
+                }
+        }
+        
+        /*Add Ke_gg to global matrix Kgg: */
+        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES);
+
+        cleanup_and_return:
+        xfree((void**)&first_gauss_area_coord);
+        xfree((void**)&second_gauss_area_coord);
+        xfree((void**)&third_gauss_area_coord);
+        xfree((void**)&gauss_weights);
+
+}
+
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::CreatePVectorThermalShelf"
+void Tria::CreatePVectorThermalShelf( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){
+
+        int i,found;
+        
+        const int  numgrids=3;
+        const int  NDOF1=1;
+        const int  numdof=numgrids*NDOF1;
+        int        doflist[numdof];
+        int        numberofdofspernode;
+        double       xyz_list[numgrids][3];
+
+        double mixed_layer_capacity;
+        double thermal_exchange_velocity;
+        double rho_water;
+        double rho_ice;
+        double heatcapacity;
+        double beta;
+        double meltingpoint;
+
+        /*inputs: */
+        double dt;
+        double pressure_list[3];
+        double pressure;
+
+        /* gaussian points: */
+        int     num_area_gauss,ig;
+        double* gauss_weights  =  NULL;
+        double* first_gauss_area_coord  =  NULL;
+        double* second_gauss_area_coord =  NULL;
+        double* third_gauss_area_coord  =  NULL;
+        double  gauss_weight;
+        double  gauss_coord[3];
+        int     dofs1[1]={0};
+
+        /*matrices: */
+        double  Jdet;
+        double  P_terms[numdof];
+        double  l1l2l3[numgrids];
+
+        double  t_pmp;
+        double  scalar_ocean;
+
+        ParameterInputs* inputs=NULL;
+
+        /*recover pointers: */
+        inputs=(ParameterInputs*)vinputs;
+        
+        /* Get node coordinates and dof list: */
+        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+        GetDofList(&doflist[0],&numberofdofspernode);
+
+        //recover material parameters
+        mixed_layer_capacity=matpar->GetMixedLayerCapacity();
+        thermal_exchange_velocity=matpar->GetThermalExchangeVelocity();
+        rho_water=matpar->GetRhoWater();
+        rho_ice=matpar->GetRhoIce();
+        heatcapacity=matpar->GetHeatCapacity();
+        beta=matpar->GetBeta();
+        meltingpoint=matpar->GetMeltingPoint();
+
+
+        /*recover extra inputs from users, dt and velocity: */
+        found=inputs->Recover("dt",&dt);
+        if((!found) && (sub_analysis_type==TransientAnalysisEnum()))throw ErrorException(__FUNCT__," could not find dt in inputs!");
+        
+        found=inputs->Recover("pressure",&pressure_list[0],1,dofs1,numgrids,(void**)nodes);
+        if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
+
+        /* Ice/ocean heat exchange flux on ice shelf base */
+
+        GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+
+        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
+        for (ig=0; ig<num_area_gauss; ig++){
+                gauss_weight=*(gauss_weights+ig);
+                gauss_coord[0]=*(first_gauss_area_coord+ig); 
+                gauss_coord[1]=*(second_gauss_area_coord+ig);
+                gauss_coord[2]=*(third_gauss_area_coord+ig);
+
+                //Get the Jacobian determinant
+                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss_coord);
+
+                /*Get nodal functions values: */
+                GetNodalFunctions(&l1l2l3[0], gauss_coord);
+
+                /*Get geothermal flux and basal friction */
+                GetParameterValue(&pressure,&pressure_list[0],gauss_coord);
+                t_pmp=meltingpoint-beta*pressure;
+
+                /*Calculate scalar parameter*/
+                scalar_ocean=gauss_weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity*t_pmp/(heatcapacity*rho_ice);
+                if(sub_analysis_type==TransientAnalysisEnum()){
+                        scalar_ocean=dt*scalar_ocean;
+                }
+
+                for(i=0;i<3;i++){
+                        P_terms[i]+=scalar_ocean*l1l2l3[i];
+                }
+        }
+
+        /*Add pe_g to global vector pg: */
+        VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
+
+        cleanup_and_return: 
+        xfree((void**)&first_gauss_area_coord);
+        xfree((void**)&second_gauss_area_coord);
+        xfree((void**)&third_gauss_area_coord);
+        xfree((void**)&gauss_weights);
+
+}
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::CreatePVectorThermalSheet"
+void Tria::CreatePVectorThermalSheet( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){
+
+        int i,found;
+        
+        const int  numgrids=3;
+        const int  NDOF1=1;
+        const int  numdof=numgrids*NDOF1;
+        int        doflist[numdof];
+        int        numberofdofspernode;
+        double       xyz_list[numgrids][3];
+        double     vxvyvz_list[numgrids][3];
+        double     vx_list[numgrids];
+        double     vy_list[numgrids];
+
+        double rho_ice;
+        double heatcapacity;
+
+        /*inputs: */
+        double dt;
+        double pressure_list[3];
+        double pressure;
+        double alpha2_list[3];
+        double basalfriction_list[3];
+        double basalfriction;
+        double geothermalflux_value;
+
+        /* gaussian points: */
+        int     num_area_gauss,ig;
+        double* gauss_weights  =  NULL;
+        double* first_gauss_area_coord  =  NULL;
+        double* second_gauss_area_coord =  NULL;
+        double* third_gauss_area_coord  =  NULL;
+        double  gauss_weight;
+        double  gauss_coord[3];
+        int     dofs1[1]={0};
+
+        /*matrices: */
+        double  Jdet;
+        double  P_terms[numdof];
+        double  l1l2l3[numgrids];
+        double  scalar;
+
+        int     dofs[3]={0,1,2};
+
+        ParameterInputs* inputs=NULL;
+
+        /*recover pointers: */
+        inputs=(ParameterInputs*)vinputs;
+        
+        /* Get node coordinates and dof list: */
+        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+        GetDofList(&doflist[0],&numberofdofspernode);
+
+        //recover material parameters
+        rho_ice=matpar->GetRhoIce();
+        heatcapacity=matpar->GetHeatCapacity();
+
+
+        /*recover extra inputs from users, dt and velocity: */
+        found=inputs->Recover("dt",&dt);
+        if((!found) && (sub_analysis_type==TransientAnalysisEnum()))throw ErrorException(__FUNCT__," could not find dt in inputs!");
+        
+        found=inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes);
+        if(!found)throw ErrorException(__FUNCT__," could not find velocity in inputs!");
+
+        for(i=0;i<numgrids;i++){
+                vx_list[i]=vxvyvz_list[i][0];
+                vy_list[i]=vxvyvz_list[i][1];
+        }       
+
+        /*Build alpha2_list used by drag stiffness matrix*/
+        Friction* friction=NewFriction();
+        
+        /*Initialize all fields: */
+        if (friction_type!=2)throw ErrorException(__FUNCT__," non-viscous friction not supported yet!");
+        
+        friction->element_type=(char*)xmalloc((strlen("3d")+1)*sizeof(char));
+        strcpy(friction->element_type,"3d");
+        
+        friction->gravity=matpar->GetG();
+        friction->rho_ice=matpar->GetRhoIce();
+        friction->rho_water=matpar->GetRhoWater();
+        friction->K=&k[0];
+        friction->bed=&b[0];
+        friction->thickness=&h[0];
+        friction->velocities=&vxvyvz_list[0][0];
+        friction->p=p;
+        friction->q=q;
+
+        /*Compute alpha2_list: */
+        FrictionGetAlpha2(&alpha2_list[0],friction);
+
+        /*Erase friction object: */
+        DeleteFriction(&friction);
+
+        /* Compute basal friction */
+        for(i=0;i<numgrids;i++){
+                basalfriction_list[i]= alpha2_list[i]*(pow(vx_list[i],2.0)+pow(vy_list[i],2.0));
+        }
+        
+        /* Ice/ocean heat exchange flux on ice shelf base */
+        GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+
+        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
+        for (ig=0; ig<num_area_gauss; ig++){
+                gauss_weight=*(gauss_weights+ig);
+                gauss_coord[0]=*(first_gauss_area_coord+ig); 
+                gauss_coord[1]=*(second_gauss_area_coord+ig);
+                gauss_coord[2]=*(third_gauss_area_coord+ig);
+
+                //Get the Jacobian determinant
+                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss_coord);
+
+                /*Get nodal functions values: */
+                GetNodalFunctions(&l1l2l3[0], gauss_coord);
+
+                /*Get geothermal flux and basal friction */
+                GetParameterValue(&geothermalflux_value,&geothermalflux[0],gauss_coord);
+                GetParameterValue(&basalfriction,&basalfriction_list[0],gauss_coord);
+
+                /*Calculate scalar parameter*/
+                scalar=gauss_weight*Jdet*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice);
+                if(sub_analysis_type==TransientAnalysisEnum()){
+                        scalar=dt*scalar;
+                }
+
+                for(i=0;i<3;i++){
+                        P_terms[i]+=scalar*l1l2l3[i];
+                }
+        }
+
+        /*Add pe_g to global vector pg: */
+        VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
+
+        cleanup_and_return: 
+        xfree((void**)&first_gauss_area_coord);
+        xfree((void**)&second_gauss_area_coord);
+        xfree((void**)&third_gauss_area_coord);
+        xfree((void**)&gauss_weights);
+
+}
+
+

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

@@ -40 +40 @@
                 double melting[3];
                 double accumulation[3];
+                double geothermalflux[3];
                 int    friction_type;
                 double p;
@@ -52 +53 @@
 
                 Tria();
-                Tria(int id,int mid,int mparid,int node_ids[3],double h[3],double s[3],double b[3],double k[3],double melting[3],double accumulation[3],
+                Tria(int id,int mid,int mparid,int node_ids[3],double h[3],double s[3],double b[3],double k[3],double melting[3],double accumulation[3],double geothermalflux[3],
                                 int friction_type,double p,double q,int shelf,double meanvel,double epsvel,double viscosity_overshoot);
                 ~Tria();
@@ -108 +109 @@
                 void  MatparConfiguration(Matpar* matpar,int matpar_offset);
                 void  ComputePressure(Vec p_g);
+                void  CreateKMatrixThermal(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
+                void  CreateKMatrixMelting(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
+                void  CreatePVectorThermalShelf( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);
+                void  CreatePVectorThermalSheet( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);