Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 4014) +++ /issm/trunk/src/c/Makefile.am (revision 4015) @@ -963,5 +963,6 @@ bin_PROGRAMS = else -bin_PROGRAMS = DiagnosticAnalysis.exe ThermalAnalysis.exe PrognosticAnalysis.exe Prognostic2Analysis.exe BalancedthicknessAnalysis.exe Balancedthickness2Analysis.exe BalancedvelocitiesAnalysis.exe TransientAnalysis.exe SteadystateAnalysis.exe SlopecomputeAnalysis.exe +bin_PROGRAMS = DiagnosticAnalysis.exe +dnl bin_PROGRAMS = DiagnosticAnalysis.exe ThermalAnalysis.exe PrognosticAnalysis.exe Prognostic2Analysis.exe BalancedthicknessAnalysis.exe Balancedthickness2Analysis.exe BalancedvelocitiesAnalysis.exe TransientAnalysis.exe SteadystateAnalysis.exe SlopecomputeAnalysis.exe endif Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 4014) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 4015) @@ -574,4 +574,8 @@ UpdateInputsFromSolutionDiagnosticHoriz( solution,analysis_type,sub_analysis_type); } + else if (sub_analysis_type==StokesAnalysisEnum){ + + UpdateInputsFromSolutionDiagnosticStokes( solution,analysis_type,sub_analysis_type); + } else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"); @@ -621,9 +625,17 @@ double vx[numvertices]; double vy[numvertices]; - int dummy; + double pressure[numvertices]; + double surface[numvertices]; + double rho_ice,g; + double xyz_list[numvertices][3]; + double gauss[numvertices][numvertices]={{1,0,0,0},{0,1,0,0},{0,0,1,0},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; + /*Get dof list: */ GetDofList(&doflist[0],&dummy); + + /*Get node data: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices); /*Use the dof list to index into the solution vector: */ @@ -638,12 +650,80 @@ } + /*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D, + *so the pressure is just the pressure at the z elevation: */ + rho_ice=matpar->GetRhoIce(); + g=matpar->GetG(); + inputs->GetParameterValues(&surface[0],&gauss[0][0],6,SurfaceEnum); + + for(i=0;iinputs->ChangeEnum(VxEnum,VxOldEnum); this->inputs->ChangeEnum(VyEnum,VyOldEnum); + this->inputs->ChangeEnum(PressureEnum,PressureOldEnum); /*Add vx and vy as inputs to the tria element: */ this->inputs->AddInput(new PentaVertexInput(VxEnum,vx)); this->inputs->AddInput(new PentaVertexInput(VyEnum,vy)); + this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure)); +} + +/*}}}*/ +/*FUNCTION Penta::UpdateInputsFromSolutionDiagnosticStokes {{{1*/ +void Penta::UpdateInputsFromSolutionDiagnosticStokes(double* solution, int analysis_type, int sub_analysis_type){ + + + int i; + + const int numvertices=6; + const int numdofpervertex=4; + const int numdof=numdofpervertex*numvertices; + + int doflist[numdof]; + double values[numdof]; + double vx[numvertices]; + double vy[numvertices]; + double vz[numvertices]; + double pressure[numvertices]; + int dummy; + double stokesreconditioning; + + /*Get dof list: */ + GetDofList(&doflist[0],&dummy); + + /*Use the dof list to index into the solution vector: */ + for(i=0;iparameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); + for(i=0;iinputs->ChangeEnum(VxEnum,VxOldEnum); + this->inputs->ChangeEnum(VyEnum,VyOldEnum); + this->inputs->ChangeEnum(VzEnum,VzOldEnum); + this->inputs->ChangeEnum(PressureEnum,PressureOldEnum); + + /*Add vx and vy as inputs to the tria element: */ + this->inputs->AddInput(new PentaVertexInput(VxEnum,vx)); + this->inputs->AddInput(new PentaVertexInput(VyEnum,vy)); + this->inputs->AddInput(new PentaVertexInput(VzEnum,vz)); + this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure)); } Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 4014) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 4015) @@ -155,4 +155,5 @@ void UpdateInputsFromSolution(double* solution, int analysis_type, int sub_analysis_type); void UpdateInputsFromSolutionDiagnosticHoriz( double* solution,int analysis_type,int sub_analysis_type); + void UpdateInputsFromSolutionDiagnosticStokes( double* solution,int analysis_type,int sub_analysis_type); void UpdateInputsFromSolutionSlopeCompute( double* solution,int analysis_type,int sub_analysis_type); void UpdateInputsFromSolutionPrognostic( double* solution,int analysis_type,int sub_analysis_type); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 4014) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 4015) @@ -537,4 +537,9 @@ double vx[numvertices]; double vy[numvertices]; + double pressure[numvertices]; + double thickness[numvertices]; + double rho_ice,g; + double gauss[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}}; + int dummy; @@ -554,12 +559,26 @@ } + /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D, + *so the pressure is just the pressure at the bedrock: */ + rho_ice=matpar->GetRhoIce(); + g=matpar->GetG(); + inputs->GetParameterValues(&thickness[0],&gauss[0][0],3,ThicknessEnum); + + for(i=0;iinputs->ChangeEnum(VxEnum,VxOldEnum); this->inputs->ChangeEnum(VyEnum,VyOldEnum); + this->inputs->ChangeEnum(PressureEnum,PressureOldEnum); /*Add vx and vy as inputs to the tria element: */ this->inputs->AddInput(new TriaVertexInput(VxEnum,vx)); this->inputs->AddInput(new TriaVertexInput(VyEnum,vy)); + this->inputs->AddInput(new TriaVertexInput(PressureEnum,pressure)); + + } Index: /issm/trunk/src/c/objects/Inputs/TriaVertexInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/TriaVertexInput.h (revision 4014) +++ /issm/trunk/src/c/objects/Inputs/TriaVertexInput.h (revision 4015) @@ -44,5 +44,4 @@ void UpdateInputsFromConstant(int constant, int name){ISSMERROR("Not implemented yet!");} void UpdateInputsFromConstant(bool constant, int name){ISSMERROR("Not implemented yet!");} - void UpdateInputsFromSolution(double* solution, int analysis_type, int sub_analysis_type){ISSMERROR("Not implemented yet!");} Index: /issm/trunk/src/c/solutions/ResetBoundaryConditions.cpp =================================================================== --- /issm/trunk/src/c/solutions/ResetBoundaryConditions.cpp (revision 4014) +++ /issm/trunk/src/c/solutions/ResetBoundaryConditions.cpp (revision 4015) @@ -7,14 +7,18 @@ #include "../EnumDefinitions/EnumDefinitions.h" -void ResetBoundaryConditions(FemModel* femmodel, int analysis_type, Vec ug){ +void ResetBoundaryConditions(FemModel* femmodel, int analysis_type, int sub_analysis_type){ int verbose=0; + Vec ug=NULL; femmodel->parameters->FindParam(&verbose,VerboseEnum); if(verbose)_printf_("%s\n"," updating boundary conditions..."); + + GetSolutionFromInputsx( &ug, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters, analysis_type,sub_analysis_type); + /*set current analysis: */ femmodel->SetCurrentAnalysis(analysis_type); - /*First, for this analysis_type, free existing boundary condition vectors: */ + /*For this analysis_type, free existing boundary condition vectors: */ //global dof set @@ -29,3 +33,6 @@ Reducevectorgtosx(&femmodel->ys[analysis_counter],femmodel->yg[analysis_counter]->vector,femmodel->nodesets[analysis_counter]); + /*Free ressources:*/ + VecFree(&ug); + } Index: /issm/trunk/src/c/solutions/diagnostic.cpp =================================================================== --- /issm/trunk/src/c/solutions/diagnostic.cpp (revision 4014) +++ /issm/trunk/src/c/solutions/diagnostic.cpp (revision 4015) @@ -66,5 +66,5 @@ _printf_("create finite element model, using analyses types statically defined above:\n"); - femmodel=new FemModel(iomodel,analyses,5); + femmodel=new FemModel(fid,analyses,5); /*get parameters: */ Index: /issm/trunk/src/c/solutions/diagnostic_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/diagnostic_core.cpp (revision 4014) +++ /issm/trunk/src/c/solutions/diagnostic_core.cpp (revision 4015) @@ -12,7 +12,4 @@ void* diagnostic_core(FemModel* femmodel){ - - /*solutions: */ - Vec ug=NULL; /*parameters: */ @@ -39,4 +36,5 @@ ReinitializeInputx(femmodel,VyEnum,QmuVyEnum); ReinitializeInputx(femmodel,VzEnum,QmuVzEnum); + ReinitializeInputx(femmodel,PressureEnum,QmuPressureEnum); } @@ -51,12 +49,10 @@ } - if(ishutter){ if(verbose)_printf_("%s\n"," computing hutter velocities..."); - solve_linear(&ug,femmodel,DiagnosticAnalysisEnum,HutterAnalysisEnum); + solve_linear(NULL,femmodel,DiagnosticAnalysisEnum,HutterAnalysisEnum); - if (ismacayealpattyn) ResetBoundaryConditions(femmodel,DiagnosticHorizAnalysisEnum,ug); - VecFree(&ug); + if (ismacayealpattyn) ResetBoundaryConditions(femmodel,DiagnosticAnalysisEnum,HorizAnalysisEnum); } @@ -80,6 +76,5 @@ if(verbose)_printf_("%s\n"," update boundary conditions for stokes using velocities previously computed..."); GetSolutionFromInputsx( &ug, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters, DiagnosticAnalysisEnum,StokesAnalysisEnum); - ResetBoundaryConditions(femmodel,DiagnosticStokesAnalysisEnum,ug); - VecFree(&ug); + ResetBoundaryConditions(femmodel,DiagnosticAnalysisEnum,StokesAnalysisEnum); if(verbose)_printf_("%s\n"," computing stokes velocities and pressure ..."); @@ -87,6 +82,3 @@ } } - - /*Free ressources: */ - VecFree(&ug); } Index: /issm/trunk/src/c/solutions/slope_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/slope_core.cpp (revision 4014) +++ /issm/trunk/src/c/solutions/slope_core.cpp (revision 4015) @@ -16,7 +16,4 @@ bool isstokes; bool ishutter; - - /*output: */ - Vec slope=NULL; /*Recover some parameters: */ Index: /issm/trunk/src/c/solutions/solutions.h =================================================================== --- /issm/trunk/src/c/solutions/solutions.h (revision 4014) +++ /issm/trunk/src/c/solutions/solutions.h (revision 4015) @@ -33,5 +33,5 @@ void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* fem,bool conserve_loads,int analysis_type,int sub_analysis_type); void diagnostic_core_linear(Vec* ppg,FemModel* fem,int analysis_type,int sub_analysis_type); -void solver_linear(Vec* ppg,FemModel* fem,int analysis_type,int sub_analysis_type); +void solver_linear(Vec* pug, FemModel* femmodel,int analysis_type,int sub_analysis_type); @@ -56,5 +56,5 @@ void CreateFemModel(FemModel* femmodel,ConstDataHandle MODEL,int analysis_type,int sub_analysis_type); //int BatchDebug(Mat* Kgg,Vec* pg,FemModel* femmodel,char* filename); -void ResetBoundaryConditions(FemModel* femmodel, int analysis_type, Vec ug); +void ResetBoundaryConditions(FemModel* femmodel, int analysis_type, int sub_analysis_type); #endif Index: /issm/trunk/src/c/solutions/solver_linear.cpp =================================================================== --- /issm/trunk/src/c/solutions/solver_linear.cpp (revision 4014) +++ /issm/trunk/src/c/solutions/solver_linear.cpp (revision 4015) @@ -8,5 +8,5 @@ #include "../modules/modules.h" -void solver_linear(Vec* pug,FemModel* fem,int analysis_type,int sub_analysis_type){ +void solver_linear(Vec* pug, FemModel* fem,int analysis_type,int sub_analysis_type){ /*parameters:*/ @@ -61,5 +61,5 @@ xfree((void**)&solver_string); - /*Assign output pointers if requested:*/ - if (pug) *pug=ug; + /*Assign output pointers:*/ + if(pug)*pug=ug; }