Index: /issm/trunk/src/c/solutions/bedslope.cpp =================================================================== --- /issm/trunk/src/c/solutions/bedslope.cpp (revision 4030) +++ /issm/trunk/src/c/solutions/bedslope.cpp (revision 4030) @@ -0,0 +1,121 @@ +/*!\file: bedslope.cpp + * \brief: bedslope computation solution + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../DataSet/DataSet.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../include/include.h" +#include "../modules/modules.h" +#include "./solutions.h" + +int main(int argc,char* *argv){ + + /*I/O: */ + FILE* fid=NULL; + char* inputfilename=NULL; + char* outputfilename=NULL; + char* lockname=NULL; + bool qmu_analysis=false; + + /*FemModel: */ + FemModel* femmodel=NULL; + + /*Results: */ + bool waitonlock=false; + + /*time*/ + double start, finish; + double start_core, finish_core; + double start_init, finish_init; + + int analyses[1]={BedSlopeComputeAnalysisEnum}; + int solution_type=BedSlopeComputeSolutionEnum; + + MODULEBOOT(); + + #if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_)) + ISSMERROR(" parallel executable was compiled without support of parallel libraries!"); + #endif + + /*Initialize Petsc and get start time*/ + PetscInitialize(&argc,&argv,(char *)0,""); + MPI_Barrier(MPI_COMM_WORLD); start=MPI_Wtime(); + + /*Size and rank: */ + MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); + MPI_Comm_size(MPI_COMM_WORLD,&num_procs); + + _printf_("recover , input file name and output file name:\n"); + inputfilename=argv[2]; + outputfilename=argv[3]; + lockname=argv[4]; + + /*Initialize femmodel structure: */ + MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime(); + + /*Open handle to data on disk: */ + fid=pfopen(inputfilename,"rb"); + + _printf_("create finite element model, using analyses types statically defined above:\n"); + femmodel=new FemModel(fid,solution_type,analyses,1); + + /*get parameters: */ + femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum); + femmodel->parameters->FindParam(&waitonlock,WaitOnLockEnum); + + MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime(); + + /*are we running the solution sequence, or a qmu wrapper around it? : */ + if(!qmu_analysis){ + + _printf_("call computational core:\n"); + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + bedslope_core(femmodel); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + + _printf_("write results to disk:\n"); + OutputResults(femmodel,outputfilename,DiagnosticAnalysisEnum); + } + else{ + + /*run qmu analysis: */ + _printf_("calling qmu analysis on bed slope core:\n"); + + #ifdef _HAVE_DAKOTA_ + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + Qmux(femmodel,SlopeComputeAnalysisEnum,NoneAnalysisEnum); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + #else + ISSMERROR(" Dakota not present, cannot do qmu!"); + #endif + } + + if (waitonlock>0){ + _printf_("write lock file:\n"); + WriteLockFile(lockname); + } + + /*Free ressources */ + delete femmodel; + + /*Get finish time and close*/ + MPI_Barrier(MPI_COMM_WORLD); finish = MPI_Wtime( ); + _printf_("\n %-34s %f seconds \n","FemModel initialization elapsed time:",finish_init-start_init); + _printf_(" %-34s %f seconds \n","Core solution elapsed time:",finish_core-start_core); + _printf_("\n %s %i hrs %i min %i sec\n\n","Total elapsed time:",int((finish-start)/3600),int(int(finish-start)%3600/60),int(finish-start)%60); + _printf_("closing MPI and Petsc\n"); + PetscFinalize(); + + /*end module: */ + MODULEEND(); + + return 0; //unix success return; +} Index: /issm/trunk/src/c/solutions/bedslope_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/bedslope_core.cpp (revision 4030) +++ /issm/trunk/src/c/solutions/bedslope_core.cpp (revision 4030) @@ -0,0 +1,46 @@ +/*!\file: bedslope_core.cpp + * \brief: core of the bed slope solution + */ + +#include "./solutions.h" +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../modules/modules.h" + +void bedslope_core(FemModel* femmodel){ + + /*parameters: */ + int verbose; + int dim; + bool isstokes; + bool ishutter; + + Vec slopex=NULL; + Vec slopey=NULL; + + /*Recover some parameters: */ + femmodel->parameters->FindParam(&verbose,VerboseEnum); + femmodel->parameters->FindParam(&dim,DimEnum); + + if(verbose)_printf_("%s\n","computing slope..."); + + /*Call on core computations: */ + solver_linear(&slopex,femmodel,SlopeComputeAnalysisEnum,BedSlopeXAnalysisEnum); + solver_linear(&slopey,femmodel,SlopeComputeAnalysisEnum,BedSlopeYAnalysisEnum); + + /*Plug solution into inputs: */ + UpdateInputsFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,slopex,SlopeComputeAnalysisEnum,BedSlopeXAnalysisEnum); + UpdateInputsFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,slopey,SlopeComputeAnalysisEnum,BedSlopeYAnalysisEnum); + + /*extrude if we are in 3D: */ + if (dim==3){ + if(verbose)_printf_("%s\n","extruding slope in 3d..."); + InputsExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedSlopeXEnum,0); + InputsExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedSlopeYEnum,0); + } + + /*Free ressources:*/ + VecFree(&slopex); + VecFree(&slopey); +} Index: sm/trunk/src/c/solutions/diagnostic_core_linear.cpp =================================================================== --- /issm/trunk/src/c/solutions/diagnostic_core_linear.cpp (revision 4029) +++ (revision ) @@ -1,65 +1,0 @@ -/*!\file: diagnostic_core_nonlinear.cpp - * \brief: core of the diagnostic solution for non linear materials - */ - -#include "../toolkits/toolkits.h" -#include "../objects/objects.h" -#include "../EnumDefinitions/EnumDefinitions.h" -#include "../modules/modules.h" - -void diagnostic_core_linear(Vec* pug,FemModel* fem,int analysis_type,int sub_analysis_type){ - - /*parameters:*/ - int kflag,pflag; - int verbose=0; - char* solver_string=NULL; - - /*output: */ - Vec ug=NULL; - Vec uf=NULL; - - /*intermediary: */ - Mat Kgg=NULL; - Mat Kff=NULL; - Mat Kfs=NULL; - Vec pg=NULL; - Vec pf=NULL; - - /*Recover parameters: */ - kflag=1; pflag=1; - fem->parameters->FindParam(&verbose,VerboseEnum); - fem->parameters->FindParam(&solver_string,SolverStringEnum); - - //*Generate system matrices - if (verbose) _printf_(" Generating matrices\n"); - SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); - - if (verbose) _printf_(" Generating penalty matrices\n"); - //*Generate penalty system matrices - PenaltySystemMatricesx(Kgg, pg,NULL,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); - - /*!Reduce matrix from g to f size:*/ - if (verbose) _printf_(" reducing matrix from g to f set\n"); - Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets); MatFree(&Kgg); - - /*!Reduce load from g to f size: */ - if (verbose) _printf_(" reducing load from g to f set\n"); - Reduceloadfromgtofx(&pf, pg, fem->Gmn, Kfs, fem->ys, fem->nodesets);VecFree(&pg); MatFree(&Kfs); - - /*Solve: */ - if (verbose) _printf_(" solving\n"); - Solverx(&uf, Kff, pf, NULL, solver_string); MatFree(&Kff); VecFree(&pf); - - //Merge back to g set - if (verbose) _printf_(" merging solution from f to g set\n"); - Mergesolutionfromftogx(&ug, uf,fem->Gmn,fem->ys,fem->nodesets);VecFree(&uf); - - //Update inputs using new solution: - UpdateInputsFromSolutionx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,ug,analysis_type, sub_analysis_type); - - /*free ressources: */ - xfree((void**)&solver_string); - - /*Assign output pointers:*/ - *pug=ug; -} Index: sm/trunk/src/c/solutions/diagnostic_core_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/solutions/diagnostic_core_nonlinear.cpp (revision 4029) +++ (revision ) @@ -1,170 +1,0 @@ -/*!\file: diagnostic_core_nonlinear.cpp - * \brief: core of the diagnostic solution for non linear materials - */ - -#include "../toolkits/toolkits.h" -#include "../objects/objects.h" -#include "../EnumDefinitions/EnumDefinitions.h" -#include "../modules/modules.h" -#include "./solutions.h" - -void diagnostic_core_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, FemModel* fem,bool conserve_loads, int analysis_type,int sub_analysis_type){ - - - /*solution : */ - Vec ug=NULL; - Vec uf=NULL; - Vec old_ug=NULL; - Vec old_uf=NULL; - DataSet* loads=NULL; - - /*intermediary: */ - Mat Kgg=NULL; - Mat Kff=NULL; - Mat Kfs=NULL; - Vec pg=NULL; - Vec pf=NULL; - int converged; - int constraints_converged; - int num_unstable_constraints; - int count; - int numberofnodes; - int min_mechanical_constraints; - int max_nonlinear_iterations; - - /*parameters:*/ - int kflag,pflag; - char* solver_string=NULL; - int verbose=0; - - //Set active analysis_type: - fem->SetActiveAnalysis(analysis_type); - - /*Recover parameters: */ - kflag=1; pflag=1; - fem->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); - fem->parameters->FindParam(&solver_string,SolverStringEnum); - fem->parameters->FindParam(&verbose,VerboseEnum); - fem->parameters->FindParam(&min_mechanical_constraints,MinMechanicalConstraintsEnum); - fem->parameters->FindParam(&max_nonlinear_iterations,MaxNonlinearIterationsEnum); - - /*Were loads requested as output? : */ - if(conserve_loads){ - loads=fem->loads->Copy(); //protect loads from being modified by the solution - } - else{ - loads=fem->loads; //modify loads in this solution - } - - count=1; - converged=0; - - /*Start non-linear iteration using input velocity: */ - GetSolutionFromInputsx(&ug, fem->elements, fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters, analysis_type, sub_analysis_type); - Reducevectorgtofx(&uf, ug, fem->nodesets); - - for(;;){ - - //save pointer to old velocity - VecFree(&old_ug);old_ug=ug; - VecFree(&old_uf);old_uf=uf; - - if (verbose) _printf_(" Generating matrices\n"); - //*Generate system matrices - SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); - - if (verbose) _printf_(" Generating penalty matrices\n"); - //*Generate penalty system matrices - PenaltySystemMatricesx(Kgg, pg,NULL,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); - - if (verbose) _printf_(" reducing matrix from g to f set\n"); - /*!Reduce matrix from g to f size:*/ - Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets); - - /*Free ressources: */ - MatFree(&Kgg); - - if (verbose) _printf_(" reducing load from g to f set\n"); - /*!Reduce load from g to f size: */ - Reduceloadfromgtofx(&pf, pg, fem->Gmn, Kfs, fem->ys, fem->nodesets); - - //no need for pg and Kfs anymore - VecFree(&pg); - MatFree(&Kfs); - - /*Solve: */ - if (verbose) _printf_(" solving\n"); - Solverx(&uf, Kff, pf, old_uf, solver_string); - - //Merge back to g set - if (verbose) _printf_(" merging solution from f to g set\n"); - Mergesolutionfromftogx(&ug, uf,fem->Gmn,fem->ys,fem->nodesets); - - //Update inputs using new solution: - UpdateInputsFromSolutionx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,ug,analysis_type, sub_analysis_type); - - //Deal with penalty loads - if (verbose) _printf_(" penalty constraints\n"); - PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,analysis_type,sub_analysis_type); - - if(verbose)_printf_(" number of unstable constraints: %i\n",num_unstable_constraints); - - /*Figure out if convergence is reached.*/ - convergence(&converged,Kff,pf,uf,old_uf,fem->parameters); - MatFree(&Kff);VecFree(&pf); - - /*add converged to inputs: */ - UpdateInputsFromConstantx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,converged,ConvergedEnum); - - //rift convergence - if (!constraints_converged) { - if (converged){ - if (num_unstable_constraints <= min_mechanical_constraints) converged=1; - else converged=0; - } - } - - /*Increase count: */ - count++; - if(converged==1)break; - if(count>=max_nonlinear_iterations){ - _printf_(" maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations); - break; - } - } - - /*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 - if(pKff0){ - - kflag=1; pflag=0; //stiffness generation only - - SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); - Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets); - MatFree(&Kgg);VecFree(&pg); - - } - - /*Delete loads only if no ouput was requested: */ - if(!conserve_loads)delete loads; - - /*clean up*/ - VecFree(&uf); - VecFree(&old_uf); - VecFree(&old_ug); - xfree((void**)&solver_string); - - /*Assign output pointers: */ - if(pug)*pug=ug; - else VecFree(&ug); - - if(pKff0)*pKff0=Kff; - if(pKfs0)*pKfs0=Kfs; - -} Index: sm/trunk/src/c/solutions/slope_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/slope_core.cpp (revision 4029) +++ (revision ) @@ -1,51 +1,0 @@ -/*!\file: slope_core.cpp - * \brief: core of the slope solution - */ - -#include "./solutions.h" -#include "../toolkits/toolkits.h" -#include "../objects/objects.h" -#include "../EnumDefinitions/EnumDefinitions.h" -#include "../modules/modules.h" - -void slope_core(FemModel* femmodel, int sub_analysis_type){ - - /*parameters: */ - int verbose; - int dim; - bool isstokes; - bool ishutter; - - /*Recover some parameters: */ - femmodel->parameters->FindParam(&verbose,VerboseEnum); - femmodel->parameters->FindParam(&dim,DimEnum); - femmodel->parameters->FindParam(&isstokes,IsStokesEnum); - femmodel->parameters->FindParam(&ishutter,IsHutterEnum); - - if(verbose)_printf_("%s\n","computing slope..."); - - /*Call on core computations: */ - solver_linear(&slope,femmodel,SlopecomputeAnalysisEnum,sub_analysis_type); - - /*Plug solution into inputs: */ - UpdateInputsFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,slope,analysis_type, sub_analysis_type); - - /*extrude if we are in 3D: */ - if (dim==3){ - if(verbose)_printf_("%s\n","extruding slope in 3d..."); - switch(sub_analysis_type){ - case SurfaceXAnalysisEnum: - InputsExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceSlopeXEnum,0); - case SurfaceYAnalysisEnum: - InputsExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceSlopeYEnum,0); - case BedXAnalysisEnum: - InputsExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedSlopeXEnum,0); - case BedYAnalysisEnum: - InputsExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedSlopeYEnum,0); - default: ISSMERROR("%s%s%s"," sub_analysis_type ",EnumAsString(sub_analysis_type)," not supported yet!"); - } - } - - /*Free ressources:*/ - VecFree(&slope); -} Index: sm/trunk/src/c/solutions/slopecompute.cpp =================================================================== --- /issm/trunk/src/c/solutions/slopecompute.cpp (revision 4029) +++ (revision ) @@ -1,109 +1,0 @@ -/*!\file: slopecompute.cpp - * \brief: slopecompute solution - */ - -#ifdef HAVE_CONFIG_H - #include "config.h" -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif - -#include "../objects/objects.h" -#include "../shared/shared.h" -#include "../DataSet/DataSet.h" -#include "../EnumDefinitions/EnumDefinitions.h" -#include "../include/include.h" -#include "../modules/modules.h" -#include "./solutions.h" - - -int main(int argc,char* *argv){ - - /*I/O: */ - FILE* fid=NULL; - char* inputfilename=NULL; - char* outputfilename=NULL; - char* lockname=NULL; - int numberofnodes; - bool waitonlock=false; - - Model* model=NULL; - - /*Results: */ - Results* results=NULL; - - Param* param=NULL; - - /*time*/ - double start, finish; - double start_core, finish_core; - double start_init, finish_init; - - MODULEBOOT(); - - #if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_)) - ISSMERROR(" parallel executable was compiled without support of parallel libraries!"); - #endif - - /*Initialize Petsc and get start time*/ - PetscInitialize(&argc,&argv,(char *)0,""); - MPI_Barrier(MPI_COMM_WORLD); start=MPI_Wtime(); - - /*Size and rank: */ - MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); - MPI_Comm_size(MPI_COMM_WORLD,&num_procs); - - _printf_("recover , input file name and output file name:\n"); - inputfilename=argv[2]; - outputfilename=argv[3]; - lockname=argv[4]; - - /*Initialize model structure: */ - MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime(); - model=new Model(); - - /*Open handle to data on disk: */ - fid=pfopen(inputfilename,"rb"); - - /*Initialize model structure: */ - model=new Model(); - - _printf_("read and create finite element model:\n"); - model->AddFormulation(fid,SlopecomputeAnalysisEnum); - - /*recover parameters: */ - model->FindParam(&waitonlock,WaitOnLockEnum); - - MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime(); - - /*run slopecompute analysis: */ - _printf_("call computational core:\n"); - MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); - results=slopecompute_core(model); - MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); - - _printf_("write results to disk:\n"); - OutputResults(results,outputfilename); - - if (waitonlock>0){ - _printf_("write lock file:\n"); - WriteLockFile(lockname); - } - - /*Free ressources:*/ - delete results; - delete model; - - /*Get finish time and close*/ - MPI_Barrier(MPI_COMM_WORLD); finish = MPI_Wtime( ); - _printf_("\n %-34s %f seconds \n","Model initialization elapsed time:",finish_init-start_init); - _printf_(" %-34s %f seconds \n","Core solution elapsed time:",finish_core-start_core); - _printf_("\n %s %i hrs %i min %i sec\n\n","Total elapsed time:",int((finish-start)/3600),int(int(finish-start)%3600/60),int(finish-start)%60); - _printf_("closing MPI and Petsc\n"); - PetscFinalize(); - - /*end module: */ - MODULEEND(); - - return 0; //unix success return; -} Index: sm/trunk/src/c/solutions/slopecompute_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/slopecompute_core.cpp (revision 4029) +++ (revision ) @@ -1,63 +1,0 @@ -/*!\file: slopecompute_core.cpp - * \brief: core of the slopecompute solution - */ - -#include "../toolkits/toolkits.h" -#include "../objects/objects.h" -#include "../shared/shared.h" -#include "../EnumDefinitions/EnumDefinitions.h" -#include "./solutions.h" -#include "../modules/modules.h" - -Results* slopecompute_core(Model* model){ - - extern int my_rank; - - /*output: */ - Results* results=NULL; - Result * result =NULL; - - /*solutions: */ - Vec sx_g=NULL; - Vec sy_g=NULL; - - /*flags: */ - int verbose=0; - int numberofdofspernode; - int numberofnodes; - int dofs[2]={1,1}; - - /*fem slopecompute model: */ - FemModel* fem_sl=NULL; - - //initialize results - results=new Results(); - - /*recover fem model: */ - fem_sl=model->GetFormulation(SlopecomputeAnalysisEnum); - - //first recover parameters common to all solutions - model->FindParam(&verbose,VerboseEnum); - model->FindParam(&numberofnodes,NumberOfNodesEnum); - model->FindParam(&numberofdofspernode,NumberOfDofsPerNodeEnum); - - _printf_("call computational core:\n"); - diagnostic_core_linear(&sx_g,fem_sl,SlopecomputeAnalysisEnum,SurfaceXAnalysisEnum); - diagnostic_core_linear(&sy_g,fem_sl,SlopecomputeAnalysisEnum,SurfaceYAnalysisEnum); - - _printf_("extrude computed slope on all layers:\n"); - FieldExtrudex( sx_g, fem_sl->elements,fem_sl->nodes, fem_sl->vertices,fem_sl->loads, fem_sl->materials,fem_sl->parameters,"slopex",0); - FieldExtrudex( sy_g, fem_sl->elements,fem_sl->nodes, fem_sl->vertices,fem_sl->loads, fem_sl->materials,fem_sl->parameters,"slopey",0); - - /*Plug results into output dataset: */ - InputToResultx(&result,fem_sl->elements,fem_sl->nodes,fem_sl->vertices, fem_sl->loads, fem_sl->materials,fem_sl->parameters,SurfaceXAnalysisEnum,results->Size()+1,0,1); results->AddObject(result); - InputToResultx(&result,fem_sl->elements,fem_sl->nodes,fem_sl->vertices, fem_sl->loads, fem_sl->materials,fem_sl->parameters,SurfaceYAnalysisEnum,results->Size()+1,0,1); results->AddObject(result); - results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(SlopecomputeAnalysisEnum))); - - /*Free ressources:*/ - VecFree(&sx_g); - VecFree(&sy_g); - - //return: - return results; -} Index: sm/trunk/src/c/solutions/solver_linear.cpp =================================================================== --- /issm/trunk/src/c/solutions/solver_linear.cpp (revision 4029) +++ (revision ) @@ -1,65 +1,0 @@ -/*!\file: solver_linear.cpp - * \brief: numerical core of linear solutions - */ - -#include "../toolkits/toolkits.h" -#include "../objects/objects.h" -#include "../EnumDefinitions/EnumDefinitions.h" -#include "../modules/modules.h" - -void solver_linear(Vec* pug, FemModel* fem,int analysis_type,int sub_analysis_type){ - - /*parameters:*/ - int kflag,pflag; - int verbose=0; - char* solver_string=NULL; - - /*output: */ - Vec ug=NULL; - Vec uf=NULL; - - /*intermediary: */ - Mat Kgg=NULL; - Mat Kff=NULL; - Mat Kfs=NULL; - Vec pg=NULL; - Vec pf=NULL; - - /*Recover parameters: */ - kflag=1; pflag=1; - fem->parameters->FindParam(&verbose,VerboseEnum); - fem->parameters->FindParam(&solver_string,SolverStringEnum); - - //*Generate system matrices - if (verbose) _printf_(" Generating matrices\n"); - SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); - - if (verbose) _printf_(" Generating penalty matrices\n"); - //*Generate penalty system matrices - PenaltySystemMatricesx(Kgg, pg,NULL,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); - - /*!Reduce matrix from g to f size:*/ - if (verbose) _printf_(" reducing matrix from g to f set\n"); - Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets); MatFree(&Kgg); - - /*!Reduce load from g to f size: */ - if (verbose) _printf_(" reducing load from g to f set\n"); - Reduceloadfromgtofx(&pf, pg, fem->Gmn, Kfs, fem->ys, fem->nodesets);VecFree(&pg); MatFree(&Kfs); - - /*Solve: */ - if (verbose) _printf_(" solving\n"); - Solverx(&uf, Kff, pf, NULL, solver_string); MatFree(&Kff); VecFree(&pf); - - //Merge back to g set - if (verbose) _printf_(" merging solution from f to g set\n"); - Mergesolutionfromftogx(&ug, uf,fem->Gmn,fem->ys,fem->nodesets);VecFree(&uf); - - //Update inputs using new solution: - UpdateInputsFromSolutionx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,ug,analysis_type, sub_analysis_type); - - /*free ressources: */ - xfree((void**)&solver_string); - - /*Assign output pointers:*/ - if(pug)*pug=ug; -} Index: /issm/trunk/src/c/solutions/surfaceslope.cpp =================================================================== --- /issm/trunk/src/c/solutions/surfaceslope.cpp (revision 4030) +++ /issm/trunk/src/c/solutions/surfaceslope.cpp (revision 4030) @@ -0,0 +1,121 @@ +/*!\file: surfaceslope.cpp + * \brief: surfaceslope computation solution + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../DataSet/DataSet.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../include/include.h" +#include "../modules/modules.h" +#include "./solutions.h" + +int main(int argc,char* *argv){ + + /*I/O: */ + FILE* fid=NULL; + char* inputfilename=NULL; + char* outputfilename=NULL; + char* lockname=NULL; + bool qmu_analysis=false; + + /*FemModel: */ + FemModel* femmodel=NULL; + + /*Results: */ + bool waitonlock=false; + + /*time*/ + double start, finish; + double start_core, finish_core; + double start_init, finish_init; + + int analyses[1]={SurfaceSlopeComputeAnalysisEnum}; + int solution_type=SurfaceSlopeComputeSolutionEnum; + + MODULEBOOT(); + + #if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_)) + ISSMERROR(" parallel executable was compiled without support of parallel libraries!"); + #endif + + /*Initialize Petsc and get start time*/ + PetscInitialize(&argc,&argv,(char *)0,""); + MPI_Barrier(MPI_COMM_WORLD); start=MPI_Wtime(); + + /*Size and rank: */ + MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); + MPI_Comm_size(MPI_COMM_WORLD,&num_procs); + + _printf_("recover , input file name and output file name:\n"); + inputfilename=argv[2]; + outputfilename=argv[3]; + lockname=argv[4]; + + /*Initialize femmodel structure: */ + MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime(); + + /*Open handle to data on disk: */ + fid=pfopen(inputfilename,"rb"); + + _printf_("create finite element model, using analyses types statically defined above:\n"); + femmodel=new FemModel(fid,solution_type,analyses,1); + + /*get parameters: */ + femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum); + femmodel->parameters->FindParam(&waitonlock,WaitOnLockEnum); + + MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime(); + + /*are we running the solution sequence, or a qmu wrapper around it? : */ + if(!qmu_analysis){ + + _printf_("call computational core:\n"); + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + surfaceslope_core(femmodel); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + + _printf_("write results to disk:\n"); + OutputResults(femmodel,outputfilename,DiagnosticAnalysisEnum); + } + else{ + + /*run qmu analysis: */ + _printf_("calling qmu analysis on surface slope core:\n"); + + #ifdef _HAVE_DAKOTA_ + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + Qmux(femmodel,SlopeComputeAnalysisEnum,NoneAnalysisEnum); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + #else + ISSMERROR(" Dakota not present, cannot do qmu!"); + #endif + } + + if (waitonlock>0){ + _printf_("write lock file:\n"); + WriteLockFile(lockname); + } + + /*Free ressources */ + delete femmodel; + + /*Get finish time and close*/ + MPI_Barrier(MPI_COMM_WORLD); finish = MPI_Wtime( ); + _printf_("\n %-34s %f seconds \n","FemModel initialization elapsed time:",finish_init-start_init); + _printf_(" %-34s %f seconds \n","Core solution elapsed time:",finish_core-start_core); + _printf_("\n %s %i hrs %i min %i sec\n\n","Total elapsed time:",int((finish-start)/3600),int(int(finish-start)%3600/60),int(finish-start)%60); + _printf_("closing MPI and Petsc\n"); + PetscFinalize(); + + /*end module: */ + MODULEEND(); + + return 0; //unix success return; +} Index: /issm/trunk/src/c/solutions/surfaceslope_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/surfaceslope_core.cpp (revision 4030) +++ /issm/trunk/src/c/solutions/surfaceslope_core.cpp (revision 4030) @@ -0,0 +1,46 @@ +/*!\file: surfaceslope_core.cpp + * \brief: core of the slope solution + */ + +#include "./solutions.h" +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../modules/modules.h" + +void surfaceslope_core(FemModel* femmodel){ + + /*parameters: */ + int verbose; + int dim; + bool isstokes; + bool ishutter; + + Vec slopex=NULL; + Vec slopey=NULL; + + /*Recover some parameters: */ + femmodel->parameters->FindParam(&verbose,VerboseEnum); + femmodel->parameters->FindParam(&dim,DimEnum); + + if(verbose)_printf_("%s\n","computing slope..."); + + /*Call on core computations: */ + solver_linear(&slopex,femmodel,SlopeComputeAnalysisEnum,SurfaceSlopeXAnalysisEnum); + solver_linear(&slopey,femmodel,SlopeComputeAnalysisEnum,SurfaceSlopeYAnalysisEnum); + + /*Plug solution into inputs: */ + UpdateInputsFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,slopex,SlopeComputeAnalysisEnum,SurfaceSlopeXAnalysisEnum); + UpdateInputsFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,slopey,SlopeComputeAnalysisEnum,SurfaceSlopeYAnalysisEnum); + + /*extrude if we are in 3D: */ + if (dim==3){ + if(verbose)_printf_("%s\n","extruding slope in 3d..."); + InputsExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceSlopeXEnum,0); + InputsExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceSlopeYEnum,0); + } + + /*Free ressources:*/ + VecFree(&slopex); + VecFree(&slopey); +} Index: sm/trunk/src/c/solutions/thermal_core_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/solutions/thermal_core_nonlinear.cpp (revision 4029) +++ (revision ) @@ -1,135 +1,0 @@ -/*!\file: thermal_core_nonlinear.cpp - * \brief: core of the thermal solution - */ - -#include "../toolkits/toolkits.h" -#include "../objects/objects.h" -#include "../EnumDefinitions/EnumDefinitions.h" -#include "../modules/modules.h" - -void thermal_core_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem,int analysis_type,int sub_analysis_type){ - - /*solution : */ - Vec tg=NULL; - Vec tf=NULL; - Vec tf_old=NULL; - double melting_offset; - - /*intermediary: */ - Mat Kgg=NULL; - Mat Kgg_nopenalty=NULL; - Mat Kff=NULL; - Mat Kfs=NULL; - Vec pg=NULL; - Vec pg_nopenalty=NULL; - Vec pf=NULL; - - int converged; - int constraints_converged; - int num_unstable_constraints; - int count; - int numberofnodes; - int min_thermal_constraints; - bool reset_penalties; - - /*parameters:*/ - int kflag,pflag; - char* solver_string=NULL; - int verbose=0; - bool lowmem=0; - - /*Recover parameters: */ - kflag=1; pflag=1; - - fem->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); - fem->parameters->FindParam(&solver_string,SolverStringEnum); - fem->parameters->FindParam(&verbose,VerboseEnum); - fem->parameters->FindParam(&lowmem,LowmemEnum); - fem->parameters->FindParam(&min_thermal_constraints,MinThermalConstraintsEnum); - - count=1; - converged=0; - - for(;;){ - - if(verbose)_printf_("%s\n","starting direct shooting method"); - - if(count==1) reset_penalties=1; else reset_penalties=0; - UpdateInputsFromConstantx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,reset_penalties,ResetPenaltiesEnum); - - //*Generate system matrices - if (!lowmem){ - - /*Compute Kgg_nopenalty and pg_nopenalty once for all: */ - if (count==1){ - SystemMatricesx(&Kgg_nopenalty, &pg_nopenalty,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); - } - - /*Copy K_gg_nopenalty into Kgg, same for pg: */ - MatDuplicate(Kgg_nopenalty,MAT_COPY_VALUES,&Kgg); - VecDuplicatePatch(&pg,pg_nopenalty); - - //apply penalties each time - PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); - } - else{ - SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); - //apply penalties - PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); - } - - /*!Reduce matrix from g to f size:*/ - Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets); - - /*Free ressources: */ - MatFree(&Kgg); - - if (verbose) _printf_(" reducing load from g to f set\n"); - /*!Reduce load from g to f size: */ - Reduceloadfromgtofx(&pf, pg, fem->Gmn, Kfs, fem->ys, fem->nodesets); - - //no need for pg and Kfs anymore - VecFree(&pg); - MatFree(&Kfs); - - /*Solve: */ - if(verbose)_printf_("%s\n","solving"); - VecFree(&tf); - Solverx(&tf, Kff, pf,tf_old, solver_string); - VecFree(&tf_old); VecDuplicatePatch(&tf_old,tf); - - //no need for Kff and pf anymore - MatFree(&Kff);VecFree(&pf);VecFree(&tg); - - if (verbose) _printf_(" merging solution from f to g set\n"); - //Merge back to g set - Mergesolutionfromftogx(&tg, tf,fem->Gmn,fem->ys,fem->nodesets); - - //Deal with penalty loads - if (verbose) _printf_(" penalty constraints\n"); - PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,analysis_type,sub_analysis_type); - - //Update inputs using new solution: - UpdateInputsFromVectorx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,tg,TemperatureEnum,VertexEnum); - UpdateInputsFromSolutionx(fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,tg,analysis_type, sub_analysis_type); - - if (!converged){ - if(verbose)_printf_("%s%i\n"," #unstable constraints = ",num_unstable_constraints); - if (num_unstable_constraints <= min_thermal_constraints)converged=1; - } - count++; - - if(converged==1)break; - } - - /*Free ressources: */ - MatFree(&Kgg_nopenalty); - VecFree(&pg_nopenalty); - VecFree(&tf); - VecFree(&tf_old); - delete solver_string; - - /*Assign output pointers: */ - *ptg=tg; - *pmelting_offset=melting_offset; -}