source: issm/oecreview/Archive/18296-19100/ISSM-18343-18344.diff

../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp

@@ -70 +70 @@
         if(stabilization==3){
                 iomodel->FetchDataToInput(elements,MasstransportSpcthicknessEnum); //for DG, we need the spc in the element
         }
+        if(stabilization==4){
+                iomodel->FetchDataToInput(elements,MasstransportSpcthicknessEnum); //for FCT, we need the spc in the element (penlaties)
+        }
 
         if(iomodel->domaintype!=Domain2DhorizontalEnum){
                 iomodel->FetchDataToInput(elements,MeshVertexonbaseEnum);
@@ -147 +150 @@
         int stabilization;
         iomodel->Constant(&stabilization,MasstransportStabilizationEnum);
 
-        /*Do not add constraints in DG, they are weakly imposed*/
+        /*Do not add constraints in DG,  they are weakly imposed*/
         if(stabilization!=3){
                 IoModelToConstraintsx(constraints,iomodel,MasstransportSpcthicknessEnum,MasstransportAnalysisEnum,P1Enum);
         }
+        /*Do not add constraints in FCT, they are imposed using penalties*/
 }/*}}}*/
 void MasstransportAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
 
@@ -674 +678 @@
 
 }/*}}}*/
 void MasstransportAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
-           _error_("not implemented yet");
+        element->GetSolutionFromInputsOneDof(solution,ThicknessEnum);
 }/*}}}*/
 void MasstransportAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
         _error_("Not implemented yet");
@@ -772 +776 @@
         }
         return;
 }/*}}}*/
+
+/*Flux Correction Transport*/
+void           MasstransportAnalysis::LumpedMassMatrix(Vector<IssmDouble>** pMlff,FemModel* femmodel){/*{{{*/
+
+        /*Intermediaries*/
+        int  configuration_type;
+
+        /*Initialize Lumped mass matrix (actually we just save its diagonal)*/
+        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+        int fsize      = femmodel->nodes->NumberOfDofs(configuration_type,FsetEnum);
+        int flocalsize = femmodel->nodes->NumberOfDofsLocal(configuration_type,FsetEnum);
+        Vector<IssmDouble>* Mlff = new Vector<IssmDouble>(flocalsize,fsize);
+
+        /*Create and assemble matrix*/
+        for(int i=0;i<femmodel->elements->Size();i++){
+                Element*       element = dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                ElementMatrix* MLe     = this->CreateMassMatrix(element);
+                if(MLe){
+                        MLe->Lump();
+                        MLe->AddDiagonalToGlobal(Mlff);
+                }
+                delete MLe;
+        }
+        Mlff->Assemble();
+
+        /*Assign output pointer*/
+        *pMlff=Mlff;
+}/*}}}*/
+ElementMatrix* MasstransportAnalysis::CreateMassMatrix(Element* element){/*{{{*/
+
+        /* Check if ice in element */
+        if(!element->IsIceInElement()) return NULL;
+
+        /*Intermediaries*/
+        IssmDouble  D,Jdet;
+        IssmDouble* xyz_list = NULL;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Initialize Element vector and other vectors*/
+        ElementMatrix* Me     = element->NewElementMatrix();
+        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinates(&xyz_list);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGauss(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                element->NodalFunctions(basis,gauss);
+
+                D=gauss->weight*Jdet;
+                TripleMultiply(basis,1,numnodes,1,
+                                        &D,1,1,0,
+                                        basis,1,numnodes,0,
+                                        &Me->values[0],1);
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(basis);
+        delete gauss;
+        return Me;
+}/*}}}*/
+void           MasstransportAnalysis::FctKMatrix(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,FemModel* femmodel){/*{{{*/
+
+        /*Output*/
+        Matrix<IssmDouble>* Kff = NULL;
+        Matrix<IssmDouble>* Kfs = NULL;
+
+        /*Initialize Jacobian Matrix*/
+        AllocateSystemMatricesx(&Kff,&Kfs,NULL,NULL,femmodel);
+
+        /*Create and assemble matrix*/
+        for(int i=0;i<femmodel->elements->Size();i++){
+                Element*       element = dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                ElementMatrix* Ke     = this->CreateFctKMatrix(element);
+                if(Ke) Ke->AddToGlobal(Kff,Kfs);
+                delete Ke;
+        }
+        Kff->Assemble();
+        Kfs->Assemble();
+
+        /*Assign output pointer*/
+        *pKff=Kff;
+        *pKfs=Kfs;
+}/*}}}*/
+ElementMatrix* MasstransportAnalysis::CreateFctKMatrix(Element* element){/*{{{*/
+
+        /* Check if ice in element */
+        if(!element->IsIceInElement()) return NULL;
+
+        /*Intermediaries */
+        IssmDouble Jdet;
+        IssmDouble vx,vy;
+        IssmDouble* xyz_list = NULL;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+        int dim      = 2;
+
+        /*Initialize Element vector and other vectors*/
+        ElementMatrix* Ke     = element->NewElementMatrix();
+        IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
+        IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
+        IssmDouble*    D      = xNewZeroInit<IssmDouble>(dim*dim);
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinates(&xyz_list);
+        Input* vxaverage_input=element->GetInput(VxEnum); _assert_(vxaverage_input);
+        Input* vyaverage_input=element->GetInput(VyEnum); _assert_(vyaverage_input);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGauss(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                GetB(B,element,dim,xyz_list,gauss);
+                GetBprime(Bprime,element,dim,xyz_list,gauss);
+                vxaverage_input->GetInputValue(&vx,gauss);
+                vyaverage_input->GetInputValue(&vy,gauss);
+
+                D[0*dim+0] = -gauss->weight*vx*Jdet;
+                D[1*dim+1] = -gauss->weight*vy*Jdet;
+
+                TripleMultiply(B,dim,numnodes,1,
+                                        D,dim,dim,0,
+                                        Bprime,dim,numnodes,0,
+                                        &Ke->values[0],1);
+
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(B);
+        xDelete<IssmDouble>(Bprime);
+        xDelete<IssmDouble>(D);
+        delete gauss;
+        return Ke;
+}/*}}}*/

../trunk-jpl/src/c/analyses/MasstransportAnalysis.h

@@ -35 +35 @@
                 void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
                 void InputUpdateFromSolution(IssmDouble* solution,Element* element);
                 void UpdateConstraints(FemModel* femmodel);
+
+                /*FCT*/
+                void           LumpedMassMatrix(Vector<IssmDouble>** pMLff,FemModel* femmodel);
+                void           FctKMatrix(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,FemModel* femmodel);
+                ElementMatrix* CreateMassMatrix(Element* element);
+                ElementMatrix* CreateFctKMatrix(Element* element);
 };
 #endif

../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp

@@ +1 @@
+/*!\file: solutionsequence_linear.cpp
+ * \brief: numerical core of linear solutions
+ */ 
+
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+
+void solutionsequence_fct(FemModel* femmodel){
+
+        /*intermediary: */
+        Vector<IssmDouble>*  Mlf = NULL;
+        Matrix<IssmDouble>*  Kff = NULL;
+        Matrix<IssmDouble>*  Kfs = NULL;
+        Vector<IssmDouble>*  ug  = NULL;
+        Vector<IssmDouble>*  uf  = NULL;
+        Vector<IssmDouble>*  ys  = NULL;
+
+        IssmDouble theta,deltat;
+        int        dof,ncols,ncols2,ncols3,rstart,rend;
+        int        configuration_type,analysis_type;
+        double     d,diagD,mi;
+        double     dmax = 0.;
+        int*       cols  = NULL;
+        int*       cols2 = NULL;
+        int*       cols3 = NULL;
+        double*    vals  = NULL;
+        double*    vals2 = NULL;
+        double*    vals3 = NULL;
+
+        /*Create analysis*/
+        MasstransportAnalysis* analysis = new MasstransportAnalysis();
+
+        /*Recover parameters: */
+        femmodel->parameters->FindParam(&deltat,TimesteppingTimeStepEnum);
+        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+        femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+        femmodel->UpdateConstraintsx();
+        theta = 0.5;
+
+        /*Create lumped mass matrix*/
+        analysis->LumpedMassMatrix(&Mlf,femmodel);
+        analysis->FctKMatrix(&Kff,&Kfs,femmodel);
+
+        /*Create Dff Matrix*/
+        #ifdef _HAVE_PETSC_
+        Mat Kff_transp = NULL;
+        Mat Dff_petsc  = NULL;
+        Mat LHS        = NULL;
+        Mat Kff_petsc  = Kff->pmatrix->matrix;
+        Vec Mlf_petsc  = Mlf->pvector->vector;
+        MatTranspose(Kff_petsc,MAT_INITIAL_MATRIX,&Kff_transp);
+        MatDuplicate(Kff_petsc,MAT_SHARE_NONZERO_PATTERN,&Dff_petsc);
+        MatGetOwnershipRange(Kff_transp,&rstart,&rend);
+        for(int row=rstart; row<rend; row++){
+                diagD = 0.;
+                MatGetRow(Kff_petsc ,row,&ncols, (const int**)&cols, (const double**)&vals);
+                MatGetRow(Kff_transp,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
+                _assert_(ncols==ncols2);
+                for(int j=0; j<ncols; j++) {
+                        _assert_(cols[j]==cols2[j]);
+                        d = max(max(-vals[j],-vals2[j]),0.);
+                        MatSetValues(Dff_petsc,1,&row,1,&cols[j],(const double*)&d,INSERT_VALUES);
+                        if(cols[j]!=row) diagD -= d;
+                }
+                MatSetValues(Dff_petsc,1,&row,1,&row,(const double*)&diagD,INSERT_VALUES);
+                MatRestoreRow(Kff_petsc, row,&ncols, (const int**)&cols, (const double**)&vals);
+                MatRestoreRow(Kff_transp,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
+        }
+        MatAssemblyBegin(Dff_petsc,MAT_FINAL_ASSEMBLY);
+        MatAssemblyEnd(  Dff_petsc,MAT_FINAL_ASSEMBLY);
+        MatFree(&Kff_transp);
+
+        /*Create LHS: [ML − theta*detlat *(K+D)^n+1]*/
+        MatDuplicate(Kff_petsc,MAT_SHARE_NONZERO_PATTERN,&LHS);
+        for(int row=rstart; row<rend; row++){
+                MatGetRow(Kff_petsc,row,&ncols, (const int**)&cols, (const double**)&vals);
+                MatGetRow(Dff_petsc,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
+                _assert_(ncols==ncols2);
+                for(int j=0; j<ncols; j++) {
+                        _assert_(cols[j]==cols2[j]);
+                        d = -theta*deltat*(vals[j] + vals2[j]);
+                        if(cols[j]==row){
+                                VecGetValues(Mlf_petsc,1,(const int*)&cols[j],&mi);
+                                d += mi;
+                        }
+                        if(fabs(d)>dmax) dmax = fabs(d);
+                        MatSetValues(LHS,1,&row,1,&cols[j],(const double*)&d,INSERT_VALUES);
+                }
+                MatRestoreRow(Kff_petsc,row,&ncols, (const int**)&cols, (const double**)&vals);
+                MatRestoreRow(Dff_petsc,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
+        }
+
+        /*Penalize Dirichlet boundary*/
+        dmax = dmax * 1.e+3;
+        for(int i=0;i<femmodel->constraints->Size();i++){
+                Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
+                if(constraint->InAnalysis(analysis_type)){
+                        constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
+                        if(dof!=-1){
+                                MatSetValues(LHS,1,&dof,1,&dof,(const double*)&dmax,INSERT_VALUES);
+                        }
+                }
+        }
+        MatAssemblyBegin(LHS,MAT_FINAL_ASSEMBLY);
+        MatAssemblyEnd(  LHS,MAT_FINAL_ASSEMBLY);
+
+        /*Create RHS: [ML + (1 − theta) deltaT L^n] u^n */
+        GetSolutionFromInputsx(&ug,femmodel);
+        Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
+        delete ug;
+        Vec u  = uf->pvector->vector;
+        Vec Ku = NULL;
+        Vec Du = NULL;
+        Vec RHS = NULL;
+        VecDuplicate(u,&Ku);
+        VecDuplicate(u,&Du);
+        VecDuplicate(u,&RHS);
+        MatMult(Kff_petsc,u,Ku);
+        MatMult(Dff_petsc,u,Du);
+        VecPointwiseMult(RHS,Mlf_petsc,u);
+        VecAXPBYPCZ(RHS,(1-theta)*deltat,(1-theta)*deltat,1,Ku,Du);// RHS = M*u + (1-theta)*deltat*K*u + (1-theta)*deltat*D*u
+        VecFree(&Ku);
+        VecFree(&Du);
+        MatFree(&Dff_petsc);
+        delete uf;
+
+        /*Penalize Dirichlet boundary*/
+        for(int i=0;i<femmodel->constraints->Size();i++){
+                Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
+                if(constraint->InAnalysis(analysis_type)){
+                        constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
+                        d = d*dmax;
+                        if(dof!=-1){
+                                VecSetValues(RHS,1,&dof,(const double*)&d,INSERT_VALUES);
+                        }
+                }
+        }
+
+        /*Go solve!*/
+        SolverxPetsc(&u,LHS,RHS,NULL,NULL, femmodel->parameters); 
+        MatFree(&LHS);
+        VecFree(&RHS);
+        uf =new Vector<IssmDouble>(u);
+        VecFree(&u);
+
+        Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys;
+        InputUpdateFromSolutionx(femmodel,ug); 
+        delete ug;  
+
+        #else
+        _error_("PETSc needs to be installed");
+        #endif
+
+        delete Mlf;
+        delete Kff;
+        delete Kfs;
+        delete analysis;
+
+}

../trunk-jpl/src/c/solutionsequences/solutionsequences.h

@@ -15 +15 @@
 void solutionsequence_hydro_nonlinear(FemModel* femmodel);
 void solutionsequence_nonlinear(FemModel* femmodel,bool conserve_loads);
 void solutionsequence_newton(FemModel* femmodel);
+void solutionsequence_fct(FemModel* femmodel);
 void solutionsequence_FScoupling_nonlinear(FemModel* femmodel,bool conserve_loads);
 void solutionsequence_linear(FemModel* femmodel);
 void solutionsequence_la(FemModel* femmodel);

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

@@ -353 +353 @@
                                         ./solutionsequences/solutionsequence_linear.cpp\
                                         ./solutionsequences/solutionsequence_nonlinear.cpp\
                                         ./solutionsequences/solutionsequence_newton.cpp\
+                                        ./solutionsequences/solutionsequence_fct.cpp\
                                         ./solutionsequences/convergence.cpp\
                                         ./classes/Options/Options.h\
                                         ./classes/Options/Options.cpp\

../trunk-jpl/src/c/cores/masstransport_core.cpp

@@ -15 +15 @@
         int    numoutputs,domaintype;
         bool   save_results;
         bool   isFS,isfreesurface,dakota_analysis;
-        int    solution_type;
+        int    solution_type,stabilization;
         char** requested_outputs = NULL;
 
         /*activate configuration*/
@@ -28 +28 @@
         femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
         femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
         femmodel->parameters->FindParam(&numoutputs,MasstransportNumRequestedOutputsEnum);
+        femmodel->parameters->FindParam(&stabilization,MasstransportStabilizationEnum);
         if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,MasstransportRequestedOutputsEnum);
 
         /*Calculate new Surface Mass Balance (SMB)*/
@@ -48 +49 @@
         }
         else{
                 if(VerboseSolution()) _printf0_("   call computational core\n");
-                solutionsequence_linear(femmodel);
+                if(stabilization==4){
+                        solutionsequence_fct(femmodel);
+                }
+                else{
+                        solutionsequence_linear(femmodel);
+                }
                 femmodel->parameters->SetParam(ThicknessEnum,InputToExtrudeEnum);
                 extrudefrombase_core(femmodel);
                 femmodel->parameters->SetParam(BaseEnum,InputToExtrudeEnum);

../trunk-jpl/src/c/classes/matrix/ElementMatrix.cpp

@@ -328 +328 @@
 
 }
 /*}}}*/
+void ElementMatrix::AddDiagonalToGlobal(Vector<IssmDouble>* pf){/*{{{*/
+
+        IssmDouble* localvalues=NULL;
+
+        /*Check that pf is not NULL*/
+        _assert_(pf); 
+
+        /*In debugging mode, check consistency (no NaN, and values not too big)*/
+        this->CheckConsistency();
+
+        if(this->dofsymmetrical){
+                /*only use row dofs to add values into global matrices: */
+
+                if(this->row_fsize){
+                        /*first, retrieve values that are in the f-set from the g-set values matrix: */
+                        localvalues=xNew<IssmDouble>(this->row_fsize);
+                        for(int i=0;i<this->row_fsize;i++){
+                                localvalues[i] = this->values[this->ncols*this->row_flocaldoflist[i]+ this->row_flocaldoflist[i]];
+                        }
+
+                        /*add local values into global  matrix, using the fglobaldoflist: */
+                        pf->SetValues(this->row_fsize,this->row_fglobaldoflist,localvalues,ADD_VAL);
+
+                        /*Free ressources:*/
+                        xDelete<IssmDouble>(localvalues);
+                }
+        }
+        else{
+                _error_("non dofsymmetrical matrix AddToGlobal routine not support yet!");
+        }
+
+}
+/*}}}*/
 void ElementMatrix::CheckConsistency(void){/*{{{*/
         /*Check element matrix values, only in debugging mode*/
         #ifdef _ISSM_DEBUG_ 

../trunk-jpl/src/c/classes/matrix/ElementMatrix.h

@@ -59 +59 @@
                 /*ElementMatrix specific routines*/
                 void AddToGlobal(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs);
                 void AddToGlobal(Matrix<IssmDouble>* Jff);
+                void AddDiagonalToGlobal(Vector<IssmDouble>* pf);
                 void Echo(void);
                 void CheckConsistency(void);
                 void StaticCondensation(int numindices,int* indices);

../trunk-jpl/src/c/classes/Constraints/SpcDynamic.h

@@ -37 +37 @@
                 /*Constraint virtual functions definitions: {{{*/
                 void   ConstrainNode(Nodes* nodes,Parameters* parameters);
                 bool   InAnalysis(int analysis_type);
+                void   PenaltyDofAndValue(int* dof,IssmDouble* value,Nodes* nodes,Parameters* parameters){_error_("not implemented yet");};
                 /*}}}*/
                 /*SpcDynamic management:{{{ */
                 int    GetNodeId();

../trunk-jpl/src/c/classes/Constraints/SpcStatic.h

@@ -36 +36 @@
                 /*Constraint virtual functions definitions: {{{*/
                 void   ConstrainNode(Nodes* nodes,Parameters* parameters);
                 bool   InAnalysis(int analysis_type);
+                void   PenaltyDofAndValue(int* dof,IssmDouble* value,Nodes* nodes,Parameters* parameters){_error_("not implemented yet");};
                 /*}}}*/
                 /*SpcStatic management:{{{ */
                 int    GetNodeId();

../trunk-jpl/src/c/classes/Constraints/SpcTransient.h

@@ -13 +13 @@
 class SpcTransient: public Constraint{
 
         private: 
-                int     sid; /*! id, to track it*/
-                int     nodeid; /*!node id*/
-                int dof; /*!component*/
-                IssmDouble* values; /*different values in time*/
-                IssmDouble* times; /*different time steps*/
-                int nsteps; /*number of time steps*/
-                int analysis_type;
+                bool        penalty;         /*Is this a penalty constraint */
+                int         sid;             /* id, to track it             */
+                int         nodeid;          /*node id                      */
+                int         dof;             /*component                    */
+                IssmDouble *values;          /*different values in time     */
+                IssmDouble *times;           /*different time steps         */
+                int         nsteps;          /*number of time steps         */
+                int         analysis_type;
 
         public:
 
@@ -38 +39 @@
                 /*Constraint virtual functions definitions: {{{*/
                 void   ConstrainNode(Nodes* nodes,Parameters* parameters);
                 bool   InAnalysis(int analysis_type);
+                void   PenaltyDofAndValue(int* dof,IssmDouble* value,Nodes* nodes,Parameters* parameters);
                 /*}}}*/
                 /*SpcTransient management:{{{ */
                 int    GetNodeId();

../trunk-jpl/src/c/classes/Constraints/Constraint.h

@@ -20 +20 @@
 
                 virtual      ~Constraint(){};
                 virtual void ConstrainNode(Nodes* nodes,Parameters* parameters)=0;
+                virtual void PenaltyDofAndValue(int* dof,IssmDouble* value,Nodes* nodes,Parameters* parameters)=0;
                 virtual bool InAnalysis(int analysis_type)=0;
 
 };

../trunk-jpl/src/c/classes/Constraints/SpcTransient.cpp

@@ -14 +14 @@
 
 /*SpcTransient constructors and destructor*/
 SpcTransient::SpcTransient(){/*{{{*/
-        sid=-1;
-        nodeid=-1;
-        dof=-1;
-        values=NULL;
-        times=NULL;
-        nsteps=-1;
-        analysis_type=-1;
+        penalty       = false;
+        sid           = -1;
+        nodeid        = -1;
+        dof           = -1;
+        values        = NULL;
+        times         = NULL;
+        nsteps        = -1;
+        analysis_type = -1;
         return;
 }
 /*}}}*/
 SpcTransient::SpcTransient(int spc_sid,int spc_nodeid, int spc_dof,int spc_nsteps, IssmDouble* spc_times, IssmDouble* spc_values,int spc_analysis_type){/*{{{*/
 
-        sid=spc_sid;
-        nodeid=spc_nodeid;
-        dof=spc_dof;
-        nsteps=spc_nsteps;
+        penalty= false;
+        sid    = spc_sid;
+        nodeid = spc_nodeid;
+        dof    = spc_dof;
+        nsteps = spc_nsteps;
         if(spc_nsteps){
-                values=xNew<IssmDouble>(spc_nsteps);
-                times=xNew<IssmDouble>(spc_nsteps);
+                values = xNew<IssmDouble>(spc_nsteps);
+                times  = xNew<IssmDouble>(spc_nsteps);
                 xMemCpy<IssmDouble>(values,spc_values,nsteps);
                 xMemCpy<IssmDouble>(times,spc_times,nsteps);
         }
@@ -68 +70 @@
         this->Echo();
 }               
 /*}}}*/
-int    SpcTransient::Id(void){ return sid; }/*{{{*/
+int  SpcTransient::Id(void){/*{{{*/
+        return sid;
+}
 /*}}}*/
 int SpcTransient::ObjectEnum(void){/*{{{*/
 
@@ -100 +104 @@
         /*Chase through nodes and find the node to which this SpcTransient applys: */
         node=(Node*)nodes->GetObjectById(NULL,nodeid);
 
-        if(node){ //in case the spc is dealing with a node on another cpu
+        if(!this->penalty && node){ //in case the spc is dealing with a node on another cpu
 
                 /*Retrieve time in parameters: */
                 parameters->FindParam(&time,TimeEnum);
@@ -136 +140 @@
         }
 }
 /*}}}*/
+void SpcTransient::PenaltyDofAndValue(int* pdof,IssmDouble* pvalue,Nodes* nodes,Parameters* parameters){/*{{{*/
 
+        if(!this->penalty) _error_("cannot return dof and value for non penalty constraint");
+
+        Node       *node  = NULL;
+        IssmDouble  time  = 0.;
+        int         i,gdof;
+        IssmDouble  alpha = -1.;
+        IssmDouble  value;
+        bool        found = false;
+
+        /*Chase through nodes and find the node to which this SpcTransient applys: */
+        node=(Node*)nodes->GetObjectById(NULL,nodeid);
+
+        if(node){ //in case the spc is dealing with a node on another cpu
+
+                /*Retrieve time in parameters: */
+                parameters->FindParam(&time,TimeEnum);
+
+                /*Now, go fetch value for this time: */
+                if (time<=times[0]){
+                        value=values[0];
+                        found=true;
+                }
+                else if (time>=times[nsteps-1]){
+                        value=values[nsteps-1];
+                        found=true;
+                }
+                else{
+                        for(i=0;i<nsteps-1;i++){
+                                if (times[i]<=time && time<times[i+1]){
+                                        alpha=(time-times[i])/(times[i+1]-times[i]);
+                                        value=(1-alpha)*values[i]+alpha*values[i+1];
+                                        found=true;
+                                        break;
+                                }
+                        }
+                }
+                if(!found)_error_("could not find time segment for constraint");
+
+                /*Get gdof */
+                gdof = node->GetDof(dof,GsetEnum);
+                if(xIsNan<IssmDouble>(value)){
+                        gdof = -1;
+                }
+        }
+        else{
+                value = NAN;
+                gdof = -1;
+        }
+
+        /*Assign output pointers*/
+        *pdof   = gdof;
+        *pvalue = value;
+}
+/*}}}*/
+
 /*SpcTransient functions*/
 int SpcTransient::GetDof(){/*{{{*/
         return dof;