Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 4012) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 4013) @@ -46,5 +46,5 @@ //slope SurfaceAnalysisEnum, - SlopecomputeAnalysisEnum, + SlopeComputeAnalysisEnum, BedXAnalysisEnum, BedYAnalysisEnum, @@ -154,6 +154,6 @@ ArtDiffEnum, BedEnum, - BedSlopexEnum, - BedSlopeyEnum, + BedSlopeXEnum, + BedSlopeYEnum, BoundaryEnum, CmMaxDmpSlopeEnum, @@ -204,6 +204,6 @@ SurfaceAreaEnum, SurfaceEnum, - SurfaceSlopexEnum, - SurfaceSlopeyEnum, + SurfaceSlopeXEnum, + SurfaceSlopeYEnum, TemperatureEnum, TemperatureAverageEnum, Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 4012) +++ /issm/trunk/src/c/Makefile.am (revision 4013) @@ -951,4 +951,6 @@ ./solutions/transient_core_3d.cpp\ ./solutions/steadystate_core.cpp\ + ./solutions/ResetBoundaryConditions.cpp\ + ./solutions/solve_linear.cpp\ ./modules/Bamgx/Bamgx.cpp\ ./modules/Bamgx/Bamgx.h\ Index: /issm/trunk/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.h =================================================================== --- /issm/trunk/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.h (revision 4012) +++ /issm/trunk/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.h (revision 4013) @@ -10,5 +10,5 @@ /* local prototypes: */ -void GetSolutionFromInputsx( Vec* psolution, DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials, Parameters* parameters, int analysis_type, int sub_analysis_type); +void GetSolutionFromInputsx( Vec* psolution, DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials, Parameters* parameters, int analysis_type, int sub_analysis_type); #endif /* _GETSOLUTIONFROMINPUTSXX_H */ Index: /issm/trunk/src/c/modules/OutputResultsx/OutputResultsx.cpp =================================================================== --- /issm/trunk/src/c/modules/OutputResultsx/OutputResultsx.cpp (revision 4012) +++ /issm/trunk/src/c/modules/OutputResultsx/OutputResultsx.cpp (revision 4013) @@ -35,3 +35,12 @@ pfclose(fid,filename); } + + + /*output if we have rifts: */ + if(numrifts){ + OutputRiftsx( &riftproperties,fem_dh->loads,numrifts); + results->AddObject(new Result(results->Size()+1,0,1,"riftproperties",riftproperties)); + } + + } Index: /issm/trunk/src/c/objects/Elements/Beam.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Beam.cpp (revision 4012) +++ /issm/trunk/src/c/objects/Elements/Beam.cpp (revision 4013) @@ -119,6 +119,6 @@ /*FUNCTION Beam::IsInput{{{1*/ bool Beam::IsInput(int name){ - if (name==SurfaceSlopexEnum || - name==SurfaceSlopeyEnum){ + if (name==SurfaceSlopeXEnum || + name==SurfaceSlopeYEnum){ return true; } @@ -412,6 +412,6 @@ //recover extra inputs - inputs->GetParameterValue(&slope[0],&gauss1[0],SurfaceSlopexEnum); - inputs->GetParameterValue(&slope[1],&gauss1[0],SurfaceSlopeyEnum); + inputs->GetParameterValue(&slope[0],&gauss1[0],SurfaceSlopeXEnum); + inputs->GetParameterValue(&slope[1],&gauss1[0],SurfaceSlopeYEnum); inputs->GetParameterValue(&surface,&gauss1[0],SurfaceEnum); inputs->GetParameterValue(&thickness,&gauss1[0],ThicknessEnum); Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 4012) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 4013) @@ -577,5 +577,5 @@ } - else if (analysis_type==SlopecomputeAnalysisEnum){ + else if (analysis_type==SlopeComputeAnalysisEnum){ UpdateInputsFromSolutionSlopeCompute( solution,analysis_type,sub_analysis_type); @@ -813,6 +813,6 @@ name==SurfaceEnum || name==BedEnum || - name==SurfaceSlopexEnum || - name==SurfaceSlopeyEnum || + name==SurfaceSlopeXEnum || + name==SurfaceSlopeYEnum || name==MeltingRateEnum || name==AccumulationRateEnum || Index: /issm/trunk/src/c/objects/Elements/Sing.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Sing.cpp (revision 4012) +++ /issm/trunk/src/c/objects/Elements/Sing.cpp (revision 4013) @@ -117,6 +117,6 @@ /*FUNCTION Sing::IsInput{{{1*/ bool Sing::IsInput(int name){ - if (name==SurfaceSlopexEnum || - name==SurfaceSlopeyEnum){ + if (name==SurfaceSlopeXEnum || + name==SurfaceSlopeYEnum){ return true; } @@ -324,6 +324,6 @@ matice=(Matice*)hmatice.delivers(); - inputs->GetParameterValue(&slope[0],SurfaceSlopexEnum); - inputs->GetParameterValue(&slope[1],SurfaceSlopeyEnum); + inputs->GetParameterValue(&slope[0],SurfaceSlopeXEnum); + inputs->GetParameterValue(&slope[1],SurfaceSlopeYEnum); GetDofList(&doflist[0],&numberofdofspernode); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 4012) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 4013) @@ -501,5 +501,5 @@ else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"); } - else if (analysis_type==SlopecomputeAnalysisEnum){ + else if (analysis_type==SlopeComputeAnalysisEnum){ UpdateInputsFromSolutionSlopeCompute( solution,analysis_type,sub_analysis_type); } @@ -693,6 +693,6 @@ /*FUNCTION Tria::IsInput{{{1*/ bool Tria::IsInput(int name){ - if (name==SurfaceSlopexEnum || - name==SurfaceSlopeyEnum){ + if (name==SurfaceSlopeXEnum || + name==SurfaceSlopeYEnum){ return true; } @@ -878,5 +878,5 @@ else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"); } - else if (analysis_type==SlopecomputeAnalysisEnum){ + else if (analysis_type==SlopeComputeAnalysisEnum){ CreateKMatrixSlopeCompute( Kgg,analysis_type,sub_analysis_type); } Index: /issm/trunk/src/c/solutions/ResetBoundaryConditions.cpp =================================================================== --- /issm/trunk/src/c/solutions/ResetBoundaryConditions.cpp (revision 4013) +++ /issm/trunk/src/c/solutions/ResetBoundaryConditions.cpp (revision 4013) @@ -0,0 +1,31 @@ +/*!\file: ResetBoundaryConditions.cpp + * \brief: change boundary conditions of a model, using a solution vector from another analysis + */ + +#include "../objects/objects.h" +#include "../modules/modules.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +void ResetBoundaryConditions(FemModel* femmodel, int analysis_type, Vec ug){ + + int verbose=0; + + femmodel->parameters->FindParam(&verbose,VerboseEnum); + if(verbose)_printf_("%s\n"," updating boundary conditions..."); + + femmodel->SetCurrentAnalysis(analysis_type); + + /*First, for this analysis_type, free existing boundary condition vectors: */ + + //global dof set + VecFree(&femmodel->yg[femmodel->analysis_counter]->vector); + //in the s-set + VecFree(&femmodel->ys[analysis_counter]); + + //Now, duplicate ug (the solution vector) into the boundary conditions vector on the g-set + VecDuplicatePatch(&femmodel->yg[analysis_counter]->vector,ug); + + //Reduce from g to s set + Reducevectorgtosx(&femmodel->ys[analysis_counter],femmodel->yg[analysis_counter]->vector,femmodel->nodesets[analysis_counter]); + +} Index: /issm/trunk/src/c/solutions/diagnostic.cpp =================================================================== --- /issm/trunk/src/c/solutions/diagnostic.cpp (revision 4012) +++ /issm/trunk/src/c/solutions/diagnostic.cpp (revision 4013) @@ -31,6 +31,4 @@ /*Results: */ - Results* results=NULL; - bool waitonlock=false; @@ -83,5 +81,5 @@ _printf_("call computational core:\n"); MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); - results=diagnostic_core(femmodel); + iagnostic_core(femmodel); MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); @@ -91,5 +89,5 @@ _printf_("call computational core:\n"); MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); - results=control_core(femmodel); + ontrol_core(femmodel); MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); @@ -97,5 +95,5 @@ _printf_("write results to disk:\n"); - OutputResults(results,outputfilename); + OutputResults(outputfilename); } else{ @@ -119,5 +117,4 @@ /*Free ressources */ delete femmodel; - delete results; /*Get finish time and close*/ Index: /issm/trunk/src/c/solutions/diagnostic_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/diagnostic_core.cpp (revision 4012) +++ /issm/trunk/src/c/solutions/diagnostic_core.cpp (revision 4013) @@ -11,55 +11,17 @@ #include "../include/include.h" -Results* diagnostic_core(FemModel* femmodel){ - - extern int my_rank; - int dummy; - - /*output: */ - Results* results=NULL; - Result* result=NULL; +void* diagnostic_core(FemModel* femmodel){ /*solutions: */ Vec ug=NULL; - Vec ug_horiz=NULL; - Vec ug_vert=NULL; - Vec ug_stokes=NULL; - Vec pg=NULL; - Vec riftproperties=NULL; - double* u_g_initial=NULL; - double* vx=NULL; - double* vy=NULL; - double* vz=NULL; - /*flags: */ - int verbose=0; + /*parameters: */ + int verbose=0; bool qmu_analysis=false; - int dim=-1; + int dim=-1; bool ishutter=false; bool ismacayealpattyn=false; bool isstokes=false; - int numberofdofspernode_sl; - int numberofdofspernode_dh; - int numberofdofspernode_ds; - int numberofnodes; - int numrifts=0; - - /*slopes: */ - Vec surfaceslopex=NULL; - Vec surfaceslopey=NULL; - Vec bedslopex=NULL; - Vec bedslopey=NULL; - double stokesreconditioning; - - /*dof recovery: */ - int dof01[2]={0,1}; - int dof2[1]={2}; - int dof012[3]={0,1,2}; - int dof3[1]={3}; - double* dofset=NULL; - - //initialize results - results=new Results(); //first recover parameters needed to drive the solution @@ -68,14 +30,7 @@ femmodel->parameters->FindParam(&ishutter,IsHutterEnum); femmodel->parameters->FindParam(&ismacayealpattyn,IsMacAyealPattynEnum); - femmodel->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); femmodel->parameters->FindParam(&isstokes,IsStokesEnum); femmodel->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); - femmodel->parameters->FindParam(&numrifts,NumRiftsEnum); femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum); - - //specific parameters for specific models - femmodel->parameters->FindParam(&numberofdofspernode_dh,NumberOfDofsPerNodeEnum); - femmodel->parameters->FindParam(&numberofdofspernode_sl,NumberOfDofsPerNodeEnum); - femmodel->parameters->FindParam(&numberofdofspernode_ds,NumberOfDofsPerNodeEnum); /*for qmu analysis, reinitialize velocity so that fake sensitivities do not show up as a result of a different restart of the convergence at each trial.*/ @@ -87,8 +42,12 @@ /*Compute slopes: */ - slope_core(femmodel,SurfaceXAnalysisEnum); - slope_core(femmodel,SurfaceYAnalysisEnum); - slope_core(femmodel,BedSlopeXAnalysisEnum); - slope_core(femmodel,BedSlopeYAnalysisEnum); + if(ishutter){ + slope_core(femmodel,SurfaceXAnalysisEnum); + slope_core(femmodel,SurfaceYAnalysisEnum); + } + if(isstokes){ + slope_core(femmodel,BedSlopeXAnalysisEnum); + slope_core(femmodel,BedSlopeYAnalysisEnum); + } @@ -96,16 +55,8 @@ if(verbose)_printf_("%s\n"," computing hutter velocities..."); - diagnostic_core_linear(&ug,fem_dhu,DiagnosticAnalysisEnum,HutterAnalysisEnum); - - if(verbose)_printf_("%s\n"," computing pressure according to MacAyeal..."); - ComputePressurex(&pg,fem_dhu->elements,fem_dhu->nodes, fem_dhu->vertices,fem_dhu->loads,fem_dhu->materials,fem_dhu->parameters,DiagnosticAnalysisEnum,HutterAnalysisEnum); - - if(verbose)_printf_("%s\n"," update boundary conditions for macyeal pattyn using hutter results..."); - if (ismacayealpattyn){ - VecFree(&fem_dh->yg->vector); VecFree(&fem_dh->ys); - VecDuplicatePatch(&fem_dh->yg->vector,ug); - Reducevectorgtosx(&fem_dh->ys,fem_dh->yg->vector,fem_dh->nodesets); - } - + solve_linear(&ug,femmodel,DiagnosticAnalysisEnum,HutterAnalysisEnum); + + if (ismacayealpattyn) ResetBoundaryConditions(femmodel,DiagnosticHorizAnalysisEnum,ug); + VecFree(&ug); } @@ -113,11 +64,5 @@ if(verbose)_printf_("%s\n"," computing horizontal velocities..."); - diagnostic_core_nonlinear(&ug,NULL,NULL,fem_dh->loads,fem_dh,DiagnosticAnalysisEnum,HorizAnalysisEnum); - - if(dim==2){ - if(verbose)_printf_("%s\n"," computing pressure according to MacAyeal..."); - ComputePressurex(&pg,fem_dh->elements,fem_dh->nodes, fem_dh->vertices,fem_dh->loads,fem_dh->materials,fem_dh->parameters,DiagnosticAnalysisEnum,HorizAnalysisEnum); - } - + diagnostic_core_nonlinear(NULL,NULL,NULL,femmodel,false,DiagnosticAnalysisEnum,HorizAnalysisEnum); //false means we expose loads to changes inside the solution } @@ -125,97 +70,23 @@ if (dim==3){ - if(verbose)_printf_("%s\n"," extruding horizontal velocities..."); - FieldExtrudex(ug,fem_dh->elements,fem_dh->nodes, fem_dh->vertices,fem_dh->loads,fem_dh->materials,fem_dh->parameters,"velocity",1); - SplitSolutionVectorx(ug,numberofnodes,numberofdofspernode_dh,&vx,&vy); - model->UpdateInputsFromVector(vx,VxEnum,VertexEnum); - model->UpdateInputsFromVector(vy,VyEnum,VertexEnum); - if(verbose)_printf_("%s\n"," computing vertical velocities..."); - diagnostic_core_linear(&ug_vert,fem_dv,DiagnosticAnalysisEnum,VertAnalysisEnum); - model->UpdateInputsFromVector(ug_vert,VzEnum,VertexEnum); + solve_linear(NULL,femmodel,DiagnosticAnalysisEnum,VertAnalysisEnum); - if(verbose)_printf_("%s\n"," computing pressure according to Pattyn..."); - ComputePressurex(&pg,fem_dh->elements, fem_dh->nodes, fem_dh->vertices,fem_dh->loads, fem_dh->materials,fem_dh->parameters,DiagnosticAnalysisEnum,HorizAnalysisEnum); - model->UpdateInputsFromVector(pg,PressureEnum,VertexEnum); - if (isstokes){ - //"recondition" pressure - VecScale(pg,1.0/stokesreconditioning); + //"recondition" pressure computed previously: + DuplicateInputx(femmodel,PressureEnum,PressureStokesEnum); ScaleInputx(femmmodel,PressureStokesEnum,1.0/stokesreconditioning); if(verbose)_printf_("%s\n"," update boundary conditions for stokes using velocities previously computed..."); - xfree((void**)&dofset);dofset=dofsetgen(2,dof01,4,numberofnodes*4); VecMerge(fem_ds->yg->vector,ug,dofset,2*numberofnodes); - xfree((void**)&dofset);dofset=dofsetgen(1,dof2,4,numberofnodes*4); VecMerge(fem_ds->yg->vector,ug_vert,dofset,1*numberofnodes); - VecFree(&fem_ds->ys); - Reducevectorgtosx(&fem_ds->ys,fem_ds->yg->vector,fem_ds->nodesets); + GetSolutionFromInputsx( &ug, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters, DiagnosticAnalysisEnum,StokesAnalysisEnum); + ResetBoundaryConditions(femmodel,DiagnosticStokesAnalysisEnum,ug); + VecFree(&ug); if(verbose)_printf_("%s\n"," computing stokes velocities and pressure ..."); - VecFree(&ug); - diagnostic_core_nonlinear(&ug,NULL,NULL,NULL,fem_ds,DiagnosticAnalysisEnum,StokesAnalysisEnum); - - //decondition" pressure - VecFree(&pg); - xfree((void**)&dofset);dofset=dofsetgen(1,dof3,4,numberofnodes*4); VecPartition(&pg, ug, dofset, numberofnodes*1); - VecScale(pg,stokesreconditioning); + diagnostic_core_nonlinear(NULL,NULL,NULL,NULL,fem_ds,DiagnosticAnalysisEnum,StokesAnalysisEnum); } } - /*Plug results into output dataset: */ - if(dim==2){ - if(ismacayealpattyn){ - InputToResultx(&result,fem_dh->elements,fem_dh->nodes,fem_dh->vertices, fem_dh->loads, fem_dh->materials,fem_dh->parameters,VxEnum,results->Size()+1,0,1); results->AddObject(result); - InputToResultx(&result,fem_dh->elements,fem_dh->nodes,fem_dh->vertices, fem_dh->loads, fem_dh->materials,fem_dh->parameters,VyEnum,results->Size()+1,0,1); results->AddObject(result); - InputToResultx(&result,fem_dh->elements,fem_dh->nodes,fem_dh->vertices, fem_dh->loads, fem_dh->materials,fem_dh->parameters,PressureEnum,results->Size()+1,0,1); results->AddObject(result); - } - else{ - InputToResultx(&result,fem_dhu->elements,fem_dhu->nodes,fem_dhu->vertices, fem_dhu->loads, fem_dhu->materials,fem_dhu->parameters,VxEnum,results->Size()+1,0,1); results->AddObject(result); - InputToResultx(&result,fem_dhu->elements,fem_dhu->nodes,fem_dhu->vertices, fem_dhu->loads, fem_dhu->materials,fem_dhu->parameters,VyEnum,results->Size()+1,0,1); results->AddObject(result); - InputToResultx(&result,fem_dhu->elements,fem_dhu->nodes,fem_dhu->vertices, fem_dhu->loads, fem_dhu->materials,fem_dhu->parameters,PressureEnum,results->Size()+1,0,1); results->AddObject(result); - } - } - else{ - if(isstokes){ - InputToResultx(&result,fem_ds->elements,fem_ds->nodes,fem_ds->vertices, fem_ds->loads, fem_ds->materials,fem_ds->parameters,VxEnum,results->Size()+1,0,1); results->AddObject(result); - InputToResultx(&result,fem_ds->elements,fem_ds->nodes,fem_ds->vertices, fem_ds->loads, fem_ds->materials,fem_ds->parameters,VyEnum,results->Size()+1,0,1); results->AddObject(result); - InputToResultx(&result,fem_ds->elements,fem_ds->nodes,fem_ds->vertices, fem_ds->loads, fem_ds->materials,fem_ds->parameters,VzEnum,results->Size()+1,0,1); results->AddObject(result); - InputToResultx(&result,fem_ds->elements,fem_ds->nodes,fem_ds->vertices, fem_ds->loads, fem_ds->materials,fem_ds->parameters,PressureEnum,results->Size()+1,0,1); results->AddObject(result); - } - else{ - if(ismacayealpattyn){ - InputToResultx(&result,fem_dh->elements,fem_dh->nodes,fem_dh->vertices, fem_dh->loads, fem_dh->materials,fem_dh->parameters,VxEnum,results->Size()+1,0,1); results->AddObject(result); - InputToResultx(&result,fem_dh->elements,fem_dh->nodes,fem_dh->vertices, fem_dh->loads, fem_dh->materials,fem_dh->parameters,VyEnum,results->Size()+1,0,1); results->AddObject(result); - InputToResultx(&result,fem_dh->elements,fem_dh->nodes,fem_dh->vertices, fem_dh->loads, fem_dh->materials,fem_dh->parameters,VzEnum,results->Size()+1,0,1); results->AddObject(result); - InputToResultx(&result,fem_dh->elements,fem_dh->nodes,fem_dh->vertices, fem_dh->loads, fem_dh->materials,fem_dh->parameters,PressureEnum,results->Size()+1,0,1); results->AddObject(result); - } - else{ - InputToResultx(&result,fem_dhu->elements,fem_dhu->nodes,fem_dhu->vertices, fem_dhu->loads, fem_dhu->materials,fem_dhu->parameters,VxEnum,results->Size()+1,0,1); results->AddObject(result); - InputToResultx(&result,fem_dhu->elements,fem_dhu->nodes,fem_dhu->vertices, fem_dhu->loads, fem_dhu->materials,fem_dhu->parameters,VyEnum,results->Size()+1,0,1); results->AddObject(result); - InputToResultx(&result,fem_dhu->elements,fem_dhu->nodes,fem_dhu->vertices, fem_dhu->loads, fem_dhu->materials,fem_dhu->parameters,VzEnum,results->Size()+1,0,1); results->AddObject(result); - InputToResultx(&result,fem_dhu->elements,fem_dhu->nodes,fem_dhu->vertices, fem_dhu->loads, fem_dhu->materials,fem_dhu->parameters,PressureEnum,results->Size()+1,0,1); results->AddObject(result); - } - - } - } - - results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(DiagnosticAnalysisEnum))); - - /*output if we have rifts: */ - if(numrifts){ - OutputRiftsx( &riftproperties,fem_dh->loads,numrifts); - results->AddObject(new Result(results->Size()+1,0,1,"riftproperties",riftproperties)); - } - /*Free ressources: */ VecFree(&ug); - VecFree(&pg); - xfree((void**)&dofset); - xfree((void**)&u_g_initial); - - VecFree(&surfaceslopex); - VecFree(&surfaceslopey); - VecFree(&bedslopex); - VecFree(&bedslopey); - - /*return: */ - return results; } Index: /issm/trunk/src/c/solutions/diagnostic_core_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/solutions/diagnostic_core_nonlinear.cpp (revision 4012) +++ /issm/trunk/src/c/solutions/diagnostic_core_nonlinear.cpp (revision 4013) @@ -9,5 +9,5 @@ #include "./solutions.h" -void diagnostic_core_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, DataSet* input_loads,FemModel* fem,int analysis_type,int sub_analysis_type){ +void diagnostic_core_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, FemModel* fem,bool conserve_loads, int analysis_type,int sub_analysis_type){ @@ -38,4 +38,7 @@ int verbose=0; + //Set active analysis_type: + fem->SetActiveAnalysis(analysis_type); + /*Recover parameters: */ kflag=1; pflag=1; @@ -47,9 +50,9 @@ /*Were loads requested as output? : */ - if(!input_loads){ - loads=fem->loads->Copy(); //we don't want to clobber loads, as they are not needed in output. + if(conserve_loads){ + loads=fem->loads->Copy(); //protect loads from being modified by the solution } else{ - loads=input_loads; //we are going to modify the loads! + loads=fem->loads; //modify loads in this solution } @@ -130,4 +133,11 @@ } } + + /*extrude if we are in 3D: */ + if (dim==3){ + if(verbose)_printf_("%s\n","extruding velocity and pressure in 3d..."); + InputsExtrudex( fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,VxEnum,0); + InputsExtrudex( fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,VyEnum,0); + } //more output might be needed, when running in control @@ -143,5 +153,5 @@ /*Delete loads only if no ouput was requested: */ - if(!input_loads)delete loads; + if(!conserve_loads)delete loads; /*clean up*/ @@ -152,5 +162,7 @@ /*Assign output pointers: */ - *pug=ug; + if(pug)*pug=ug; + else VecFree(&ug); + if(pKff0)*pKff0=Kff; if(pKfs0)*pKfs0=Kfs; Index: /issm/trunk/src/c/solutions/slope_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/slope_core.cpp (revision 4012) +++ /issm/trunk/src/c/solutions/slope_core.cpp (revision 4013) @@ -28,8 +28,4 @@ if(verbose)_printf_("%s\n","computing slope..."); - /*Early return possible? */ - if(!ishutter && AnalysisEnum==SurfaceAnalysisEnum) return; - if(!isstokes && AnalysisEnum==BedAnalysisEnum) return; - /*Call on core computations: */ solver_linear(&slope,femmodel,SlopecomputeAnalysisEnum,sub_analysis_type); Index: /issm/trunk/src/c/solutions/solutions.h =================================================================== --- /issm/trunk/src/c/solutions/solutions.h (revision 4012) +++ /issm/trunk/src/c/solutions/solutions.h (revision 4013) @@ -14,5 +14,5 @@ /*cores: */ Results* gradjcompute_core(Model* model); -Results* diagnostic_core(Model* model); +void* diagnostic_core(Model* model); Results* prognostic_core(Model* model); Results* prognostic2_core(Model* model); @@ -27,11 +27,11 @@ Results* transient_core_3d(Model* model); Results* thermal_core(Model* model); +void slope_core(FemModel* fem,int sub_analysis_type); /*computational cores: */ void thermal_core_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem,int analysis_type,int sub_analysis_type); -void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, DataSet* loads, FemModel* fem,int analysis_type,int sub_analysis_type); +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 slope_core(FemModel* fem,int sub_analysis_type); @@ -56,4 +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); #endif Index: /issm/trunk/src/c/solutions/solver_linear.cpp =================================================================== --- /issm/trunk/src/c/solutions/solver_linear.cpp (revision 4012) +++ /issm/trunk/src/c/solutions/solver_linear.cpp (revision 4013) @@ -61,5 +61,5 @@ xfree((void**)&solver_string); - /*Assign output pointers:*/ - *pug=ug; + /*Assign output pointers if requested:*/ + if (pug) *pug=ug; }