Index: ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 18343) +++ ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 18344) @@ -70,6 +70,9 @@ 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,10 +150,11 @@ 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,7 +678,7 @@ }/*}}}*/ void MasstransportAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ - _error_("not implemented yet"); + element->GetSolutionFromInputsOneDof(solution,ThicknessEnum); }/*}}}*/ void MasstransportAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ _error_("Not implemented yet"); @@ -772,3 +776,148 @@ } return; }/*}}}*/ + +/*Flux Correction Transport*/ +void MasstransportAnalysis::LumpedMassMatrix(Vector** 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* Mlff = new Vector(flocalsize,fsize); + + /*Create and assemble matrix*/ + for(int i=0;ielements->Size();i++){ + Element* element = dynamic_cast(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(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();igend();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(xyz_list); + xDelete(basis); + delete gauss; + return Me; +}/*}}}*/ +void MasstransportAnalysis::FctKMatrix(Matrix** pKff,Matrix** pKfs,FemModel* femmodel){/*{{{*/ + + /*Output*/ + Matrix* Kff = NULL; + Matrix* Kfs = NULL; + + /*Initialize Jacobian Matrix*/ + AllocateSystemMatricesx(&Kff,&Kfs,NULL,NULL,femmodel); + + /*Create and assemble matrix*/ + for(int i=0;ielements->Size();i++){ + Element* element = dynamic_cast(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(dim*numnodes); + IssmDouble* Bprime = xNew(dim*numnodes); + IssmDouble* D = xNewZeroInit(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();igend();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(xyz_list); + xDelete(B); + xDelete(Bprime); + xDelete(D); + delete gauss; + return Ke; +}/*}}}*/ Index: ../trunk-jpl/src/c/analyses/MasstransportAnalysis.h =================================================================== --- ../trunk-jpl/src/c/analyses/MasstransportAnalysis.h (revision 18343) +++ ../trunk-jpl/src/c/analyses/MasstransportAnalysis.h (revision 18344) @@ -35,5 +35,11 @@ void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); void InputUpdateFromSolution(IssmDouble* solution,Element* element); void UpdateConstraints(FemModel* femmodel); + + /*FCT*/ + void LumpedMassMatrix(Vector** pMLff,FemModel* femmodel); + void FctKMatrix(Matrix** pKff,Matrix** pKfs,FemModel* femmodel); + ElementMatrix* CreateMassMatrix(Element* element); + ElementMatrix* CreateFctKMatrix(Element* element); }; #endif Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp =================================================================== --- ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp (revision 0) +++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp (revision 18344) @@ -0,0 +1,161 @@ +/*!\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* Mlf = NULL; + Matrix* Kff = NULL; + Matrix* Kfs = NULL; + Vector* ug = NULL; + Vector* uf = NULL; + Vector* 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; rowdmax) 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;iconstraints->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;iconstraints->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(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; + +} Index: ../trunk-jpl/src/c/solutionsequences/solutionsequences.h =================================================================== --- ../trunk-jpl/src/c/solutionsequences/solutionsequences.h (revision 18343) +++ ../trunk-jpl/src/c/solutionsequences/solutionsequences.h (revision 18344) @@ -15,6 +15,7 @@ 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); Index: ../trunk-jpl/src/c/Makefile.am =================================================================== --- ../trunk-jpl/src/c/Makefile.am (revision 18343) +++ ../trunk-jpl/src/c/Makefile.am (revision 18344) @@ -353,6 +353,7 @@ ./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\ Index: ../trunk-jpl/src/c/cores/masstransport_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/masstransport_core.cpp (revision 18343) +++ ../trunk-jpl/src/c/cores/masstransport_core.cpp (revision 18344) @@ -15,7 +15,7 @@ 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,6 +28,7 @@ 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,7 +49,12 @@ } 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); Index: ../trunk-jpl/src/c/classes/matrix/ElementMatrix.cpp =================================================================== --- ../trunk-jpl/src/c/classes/matrix/ElementMatrix.cpp (revision 18343) +++ ../trunk-jpl/src/c/classes/matrix/ElementMatrix.cpp (revision 18344) @@ -328,6 +328,39 @@ } /*}}}*/ +void ElementMatrix::AddDiagonalToGlobal(Vector* 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(this->row_fsize); + for(int i=0;irow_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(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_ Index: ../trunk-jpl/src/c/classes/matrix/ElementMatrix.h =================================================================== --- ../trunk-jpl/src/c/classes/matrix/ElementMatrix.h (revision 18343) +++ ../trunk-jpl/src/c/classes/matrix/ElementMatrix.h (revision 18344) @@ -59,6 +59,7 @@ /*ElementMatrix specific routines*/ void AddToGlobal(Matrix* Kff, Matrix* Kfs); void AddToGlobal(Matrix* Jff); + void AddDiagonalToGlobal(Vector* pf); void Echo(void); void CheckConsistency(void); void StaticCondensation(int numindices,int* indices); Index: ../trunk-jpl/src/c/classes/Constraints/SpcDynamic.h =================================================================== --- ../trunk-jpl/src/c/classes/Constraints/SpcDynamic.h (revision 18343) +++ ../trunk-jpl/src/c/classes/Constraints/SpcDynamic.h (revision 18344) @@ -37,6 +37,7 @@ /*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(); Index: ../trunk-jpl/src/c/classes/Constraints/SpcStatic.h =================================================================== --- ../trunk-jpl/src/c/classes/Constraints/SpcStatic.h (revision 18343) +++ ../trunk-jpl/src/c/classes/Constraints/SpcStatic.h (revision 18344) @@ -36,6 +36,7 @@ /*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(); Index: ../trunk-jpl/src/c/classes/Constraints/SpcTransient.h =================================================================== --- ../trunk-jpl/src/c/classes/Constraints/SpcTransient.h (revision 18343) +++ ../trunk-jpl/src/c/classes/Constraints/SpcTransient.h (revision 18344) @@ -13,13 +13,14 @@ 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,6 +39,7 @@ /*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(); Index: ../trunk-jpl/src/c/classes/Constraints/Constraint.h =================================================================== --- ../trunk-jpl/src/c/classes/Constraints/Constraint.h (revision 18343) +++ ../trunk-jpl/src/c/classes/Constraints/Constraint.h (revision 18344) @@ -20,6 +20,7 @@ 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; }; Index: ../trunk-jpl/src/c/classes/Constraints/SpcTransient.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Constraints/SpcTransient.cpp (revision 18343) +++ ../trunk-jpl/src/c/classes/Constraints/SpcTransient.cpp (revision 18344) @@ -14,25 +14,27 @@ /*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(spc_nsteps); - times=xNew(spc_nsteps); + values = xNew(spc_nsteps); + times = xNew(spc_nsteps); xMemCpy(values,spc_values,nsteps); xMemCpy(times,spc_times,nsteps); } @@ -68,7 +70,9 @@ this->Echo(); } /*}}}*/ -int SpcTransient::Id(void){ return sid; }/*{{{*/ +int SpcTransient::Id(void){/*{{{*/ + return sid; +} /*}}}*/ int SpcTransient::ObjectEnum(void){/*{{{*/ @@ -100,7 +104,7 @@ /*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,7 +140,63 @@ } } /*}}}*/ +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;iGetDof(dof,GsetEnum); + if(xIsNan(value)){ + gdof = -1; + } + } + else{ + value = NAN; + gdof = -1; + } + + /*Assign output pointers*/ + *pdof = gdof; + *pvalue = value; +} +/*}}}*/ + /*SpcTransient functions*/ int SpcTransient::GetDof(){/*{{{*/ return dof; Index: ../trunk-jpl/src/c =================================================================== --- ../trunk-jpl/src/c (revision 18343) +++ ../trunk-jpl/src/c (revision 18344) Property changes on: ../trunk-jpl/src/c ___________________________________________________________________ Modified: svn:ignore ## -4,6 +4,8 ## *.o *.obj *.exe +issm +kriging appscan.* ouncemake* list