Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 433) +++ /issm/trunk/src/c/Makefile.am (revision 434) @@ -77,4 +77,6 @@ ./shared/Matrix/matrix.h\ ./shared/Matrix/MatrixUtils.cpp\ + ./shared/Dofs/dofs.h\ + ./shared/Dofs/dofsetgen.cpp\ ./shared/Numerics/numerics.h\ ./shared/Numerics/GaussPoints.h\ @@ -125,4 +127,5 @@ ./toolkits/petsc/patches/NewMat.cpp\ ./toolkits/petsc/patches/VecFree.cpp\ + ./toolkits/petsc/patches/VecDuplicatePatch.cpp\ ./toolkits/petsc/patches/KSPFree.cpp\ ./toolkits/petsc/patches/ISFree.cpp\ @@ -320,4 +323,6 @@ ./shared/Matrix/matrix.h\ ./shared/Matrix/MatrixUtils.cpp\ + ./shared/Dofs/dofs.h\ + ./shared/Dofs/dofsetgen.cpp\ ./shared/Numerics/numerics.h\ ./shared/Numerics/GaussPoints.h\ @@ -368,4 +373,5 @@ ./toolkits/petsc/patches/NewMat.cpp\ ./toolkits/petsc/patches/VecFree.cpp\ + ./toolkits/petsc/patches/VecDuplicatePatch.cpp\ ./toolkits/petsc/patches/KSPFree.cpp\ ./toolkits/petsc/patches/ISFree.cpp\ @@ -488,4 +494,6 @@ ./SlopeExtrudex/SlopeExtrudex.cpp\ ./SlopeExtrudex/SlopeExtrudex.h\ + ./parallel/diagnostic_core.cpp\ + ./parallel/diagnostic_core_linear.cpp\ ./parallel/diagnostic_core_nonlinear.cpp\ ./parallel/CreateFemModel.cpp\ Index: /issm/trunk/src/c/parallel/GradJCompute.cpp =================================================================== --- /issm/trunk/src/c/parallel/GradJCompute.cpp (revision 433) +++ /issm/trunk/src/c/parallel/GradJCompute.cpp (revision 434) @@ -46,5 +46,5 @@ _printf_("%s\n"," recover solution for this stiffness and right hand side:"); - diagnostic_core_nonlinear(&u_g,&K_ff0,&K_fs0,inputs,femmodel); + diagnostic_core_nonlinear(&u_g,&K_ff0,&K_fs0,femmodel,inputs,analysis_type); VecToMPISerial(&u_g_double,u_g); VecFree(&u_g); Index: /issm/trunk/src/c/parallel/OutputDiagnostic.cpp =================================================================== --- /issm/trunk/src/c/parallel/OutputDiagnostic.cpp (revision 433) +++ /issm/trunk/src/c/parallel/OutputDiagnostic.cpp (revision 434) @@ -11,5 +11,5 @@ #include "../objects/objects.h" -void OutputDiagnostic(Vec u_g,Vec partition,char* filename,NodeSets* nodesets,char* analysis_type){ +void OutputDiagnostic(Vec u_g,Vec p_g, FemModel* femmodels,char* filename){ int i; @@ -19,8 +19,13 @@ FILE* fid=NULL; + /*intermediary: */ + FemModel femmodel; + NodeSets* nodesets=NULL; + Vec partition=NULL; + char* analysis_type="diagnostic"; + /* standard output: */ Vec partition_shifted=NULL; double* serial_partition=NULL; - int analysis_size; double* serial_ug=NULL; @@ -28,4 +33,9 @@ int gsize; int nods; + + /*Recover diagnostic horiz femmodel: */ + femmodel=femmodels[0]; + partition=femmodel.partition; + nodesets=femmodel.nodesets; /*serialize outputs: */ Index: /issm/trunk/src/c/parallel/control.cpp =================================================================== --- /issm/trunk/src/c/parallel/control.cpp (revision 433) +++ /issm/trunk/src/c/parallel/control.cpp (revision 434) @@ -25,5 +25,5 @@ char* outputfilename=NULL; char* lockname=NULL; - char* analysis_type="control"; + int analysis_type; /*Intermediary: */ @@ -78,5 +78,5 @@ _printf_("read and create finite element model:\n"); - CreateFemModel(&femmodel,fid,analysis_type); + CreateFemModel(&femmodel,fid,"control"); /*Recover parameters used throughout the solution:*/ @@ -92,4 +92,5 @@ femmodel.parameters->FindParam((void*)&p_g,"p_g"); femmodel.parameters->FindParam((void*)&u_g_obs,"u_g_obs"); + femmodel.parameters->FindParam((void*)&analysis_type,"analysis_type"); gsize=femmodel.nodes->NumberOfDofs(); @@ -155,5 +156,5 @@ inputs->Add(control_type,p_g,2,numberofnodes); inputs->Add("fit",fit[n]); - diagnostic_core_nonlinear(&u_g,NULL,NULL,inputs,&femmodel); + diagnostic_core_nonlinear(&u_g,NULL,NULL,&femmodel,inputs,analysis_type); OutputControl(u_g,p_g,J,nsteps,femmodel.partition,outputfilename,femmodel.nodesets); _printf_("%s\n"," done."); @@ -174,5 +175,5 @@ UpdateFromInputsx(femmodel.elements,femmodel.nodes,femmodel.loads, femmodel.materials,inputs); - diagnostic_core_nonlinear(&u_g,NULL,NULL,inputs,&femmodel); + diagnostic_core_nonlinear(&u_g,NULL,NULL,&femmodel,inputs,analysis_type); _printf_("%s\n"," saving final results..."); Index: /issm/trunk/src/c/parallel/diagnostic.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic.cpp (revision 433) +++ /issm/trunk/src/c/parallel/diagnostic.cpp (revision 434) @@ -23,10 +23,11 @@ char* outputfilename=NULL; char* lockname=NULL; - char* analysis_type="diagnostic_horiz"; int numberofnodes; - /*Intermediary: */ - FemModel femmodel; + /*Fem models : */ + FemModel femmodels[5]; + Vec u_g=NULL; + Vec p_g=NULL; ParameterInputs* inputs=NULL; int waitonlock=0; @@ -57,24 +58,29 @@ _printf_("read and create finite element model:\n"); - CreateFemModel(&femmodel,fid,analysis_type); + CreateFemModel(&femmodels[0],fid,"diagnostic_horiz"); + CreateFemModel(&femmodels[1],fid,"diagnostic_vert"); + CreateFemModel(&femmodels[2],fid,"diagnostic_stokes"); + CreateFemModel(&femmodels[3],fid,"diagnostic_hutter"); + CreateFemModel(&femmodels[4],fid,"slope_compute"); _printf_("initialize inputs:\n"); - femmodel.parameters->FindParam((void*)&u_g_initial,"u_g"); - femmodel.parameters->FindParam((void*)&numberofnodes,"numberofnodes"); + femmodels[0].parameters->FindParam((void*)&u_g_initial,"u_g"); + femmodels[0].parameters->FindParam((void*)&numberofnodes,"numberofnodes"); inputs=new ParameterInputs; inputs->Add("velocity",u_g_initial,3,numberofnodes); - param=(Param*)femmodel.parameters->FindParamObject("u_g"); - femmodel.parameters->DeleteObject((Object*)param); + /*erase velocities: */ + param=(Param*)femmodels[0].parameters->FindParamObject("u_g"); + femmodels[0].parameters->DeleteObject((Object*)param); _printf_("call computational core:\n"); - diagnostic_core_nonlinear(&u_g,NULL,NULL,inputs,&femmodel); + diagnostic_core(&u_g,&p_g,femmodels,inputs); _printf_("write results to disk:\n"); - OutputDiagnostic(u_g,femmodel.partition,outputfilename,femmodel.nodesets,analysis_type); + OutputDiagnostic(u_g,p_g,&femmodels[0],outputfilename); _printf_("write lock file:\n"); - femmodel.parameters->FindParam((void*)&waitonlock,"waitonlock"); + femmodels[0].parameters->FindParam((void*)&waitonlock,"waitonlock"); if (waitonlock){ WriteLockFile(lockname); Index: /issm/trunk/src/c/parallel/diagnostic_core.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic_core.cpp (revision 434) +++ /issm/trunk/src/c/parallel/diagnostic_core.cpp (revision 434) @@ -0,0 +1,174 @@ +/*!\file: diagnostic_core.cpp + * \brief: core of the diagnostic solution + */ + +#undef __FUNCT__ +#define __FUNCT__ "cielodiagnostic_core" + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./parallel.h" +#include "../issm.h" + +void diagnostic_core(Vec* pug, Vec* ppg,FemModel* fems, ParameterInputs* inputs){ + + /*fem models: */ + FemModel* fem_dh=NULL; + FemModel* fem_dv=NULL; + FemModel* fem_dhu=NULL; + FemModel* fem_ds=NULL; + FemModel* fem_sl=NULL; + + /*solutions: */ + Vec ug=NULL; + Vec ug_horiz=NULL; + Vec ug_vert=NULL; + Vec ug_stokes=NULL; + Vec pg=NULL; + Vec slopex=NULL; + Vec slopey=NULL; + + /*flags: */ + int debug=0; + int dim=-1; + int ishutter=0; + int ismacayealpattyn=0; + int isstokes=0; + int numberofdofspernode_sl; + int numberofdofspernode_dh; + int numberofdofspernode_ds; + int numberofnodes; + + double stokesreconditioning; + + /*dof recovery: */ + int dof01[2]={0,1}; + int dof2[1]={2}; + int dof012[3]={0,1,2}; + int dof3[1]={3}; + + /*recover fem models: */ + fem_dh=fems+0; + fem_dv=fems+1; + fem_ds=fems+2; + fem_dhu=fems+3; + fem_sl=fems+4; + + //first recover parameters common to all solutions + fem_dh->parameters->FindParam((void*)&debug,"debug"); + fem_dh->parameters->FindParam((void*)&dim,"dim"); + fem_dhu->parameters->FindParam((void*)&ishutter,"ishutter"); + fem_dh->parameters->FindParam((void*)&ismacayealpattyn,"ismacayealpattyn"); + fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes"); + fem_ds->parameters->FindParam((void*)&isstokes,"isstokes"); + fem_ds->parameters->FindParam((void*)&stokesreconditioning,"stokesreconditioning"); + + //specific parameters for specific models + fem_dh->parameters->FindParam((void*)&numberofdofspernode_dh,"numberofdofspernode"); + fem_sl->parameters->FindParam((void*)&numberofdofspernode_sl,"numberofdofspernode"); + fem_ds->parameters->FindParam((void*)&numberofdofspernode_ds,"numberofdofspernode"); + + if(ishutter){ + + if(debug)_printf_("%s\n","computing surface slope (x and y derivatives)..."); + diagnostic_core_linear(&slopex,fem_sl,inputs,SurfaceSlopeComputeXAnalysisEnum()); + diagnostic_core_linear(&slopey,fem_sl,inputs,SurfaceSlopeComputeYAnalysisEnum()); + + if (dim==3){ + + if(debug)_printf_("%s\n","extruding slopes in 3d..."); + SlopeExtrudex( slopex, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials); + SlopeExtrudex( slopey, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials); + } + + if(debug)_printf_("%s\n"," adding slopes in inputs..."); + inputs->Add("surfaceslopex",slopex,numberofdofspernode_sl,numberofnodes); + inputs->Add("surfaceslopey",slopey,numberofdofspernode_sl,numberofnodes); + + if(debug)_printf_("%s\n"," computing hutter velocities..."); + diagnostic_core_linear(&ug,fem_dhu,inputs,DiagnosticHutterAnalysisEnum()); + + if(debug)_printf_("%s\n"," computing pressure according to MacAyeal..."); + ComputePressurex( &pg,fem_dhu->elements,fem_dhu->nodes,fem_dhu->loads,fem_dhu->materials, numberofnodes); + + if(debug)_printf_("%s\n"," update boundary conditions for macyeal pattyn using hutter results..."); + if (ismacayealpattyn){ + VecFree(&fem_dh->yg); VecFree(&fem_dh->ys);VecFree(&fem_dh->ys0); + VecDuplicatePatch(&fem_dh->yg,ug); + Reducevectorgtosx(&fem_dh->ys,&fem_dh->ys0, fem_dh->yg,fem_dh->nodesets); + } + + } + + if (ismacayealpattyn){ + + if(debug)_printf_("%s\n"," computing horizontal velocities..."); + diagnostic_core_nonlinear(&ug,NULL,NULL,fem_dh,inputs,DiagnosticHorizAnalysisEnum()); + + if(debug)_printf_("%s\n"," computing pressure according to MacAyeal..."); + ComputePressurex( &pg,fem_dh->elements, fem_dh->nodes, fem_dh->loads, fem_dh->materials, numberofnodes); + + } + + + if (dim==3){ + + if(debug)_printf_("%s\n"," extruding horizontal velocities..."); + VecDuplicatePatch(&ug_horiz,ug); VelocityExtrudex( ug_horiz,fem_dh->elements,fem_dh->nodes, fem_dh->loads,fem_dh-> materials); + + if(debug)_printf_("%s\n"," computing vertical velocities..."); + inputs->Add("velocity",ug_horiz,numberofdofspernode_dh,numberofnodes); + diagnostic_core_linear(&ug_vert,fem_dv,inputs,DiagnosticVertAnalysisEnum()); + + if(debug)_printf_("%s\n"," combining horizontal and vertical velocities..."); + VecFree(&ug); ug=NewVec(numberofnodes*3); + + VecMerge(ug,ug_horiz,dofsetgen(2,&dof01[0],3,numberofnodes*3),numberofnodes*2); + VecMerge(ug,ug_vert,dofsetgen(1,&dof2[0],3,numberofnodes*3),numberofnodes*1); + + if(debug)_printf_("%s\n"," computing pressure according to Pattyn..."); + ComputePressurex( &pg,fem_dh->elements, fem_dh->nodes, fem_dh->loads, fem_dh->materials, numberofnodes); + + if (isstokes){ + + //"recondition" pressure + VecScale(pg,1.0/stokesreconditioning); + + if(debug)_printf_("%s\n","computing bed slope (x and y derivatives)..."); + diagnostic_core_linear(&slopex,fem_sl,inputs,BedSlopeComputeXAnalysisEnum()); + diagnostic_core_linear(&slopey,fem_sl,inputs,BedSlopeComputeYAnalysisEnum()); + SlopeExtrudex( slopex, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials); + SlopeExtrudex( slopey, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials); + + inputs->Add("bedslopex",slopex,numberofdofspernode_sl,numberofnodes); + inputs->Add("bedslopey",slopey,numberofdofspernode_sl,numberofnodes); + + //recombine ug and pg: + ug_stokes=NewVec(fem_ds->nodesets->GetGSize()); + VecMerge(ug_stokes,ug,dofsetgen(3,dof012,4,numberofnodes*4),numberofnodes*3); + VecMerge(ug_stokes,pg,dofsetgen(1,dof3,4,numberofnodes*4),numberofnodes); + + inputs->Add("velocity",ug_stokes,numberofdofspernode_ds,numberofnodes); + + if(debug)_printf_("%s\n"," update boundary conditions for stokes using velocities previously computed..."); + VecMerge(fem_ds->yg,ug,dofsetgen(3,dof012,4,numberofnodes*4),3*numberofnodes); + VecFree(&fem_ds->ys); VecFree(&fem_ds->ys0); + Reducevectorgtosx(&fem_ds->ys,&fem_ds->ys0, fem_ds->yg,fem_ds->nodesets); + + if(debug)_printf_("%s\n"," computing stokes velocities and pressure ..."); + VecFree(&ug); + diagnostic_core_nonlinear(&ug,NULL,NULL,fem_ds,inputs,DiagnosticStokesAnalysisEnum()); + + //decondition" pressure + VecFree(&pg); + VecPartition(&pg, ug, dofsetgen(1,dof3,4,numberofnodes*4), numberofnodes*1); + VecScale(pg,stokesreconditioning); + } + } + + /*Assign output pointers: */ + *pug=ug; + *ppg=pg; +} Index: /issm/trunk/src/c/parallel/diagnostic_core_linear.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic_core_linear.cpp (revision 434) +++ /issm/trunk/src/c/parallel/diagnostic_core_linear.cpp (revision 434) @@ -0,0 +1,63 @@ +/*!\file: diagnostic_core_nonlinear.cpp + * \brief: core of the diagnostic solution for non linear materials + */ + +#undef __FUNCT__ +#define __FUNCT__ "diagnostic_core_inear" + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../issm.h" + +void diagnostic_core_linear(Vec* pug,FemModel* fem,ParameterInputs* inputs,int analysis_type){ + + /*parameters:*/ + int kflag,pflag,connectivity,numberofdofspernode; + int debug=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((void*)&connectivity,"connectivity"); + fem->parameters->FindParam((void*)&numberofdofspernode,"numberofdofspernode"); + fem->parameters->FindParam((void*)&debug,"debug"); + fem->parameters->FindParam((void*)&solver_string,"solverstring"); + + /*Update parameters: */ + UpdateFromInputsx(fem->elements,fem->nodes,fem->loads, fem->materials,inputs); + + //*Generate system matrices + if (debug) _printf_(" Generating matrices\n"); + SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type); + + /*!Reduce matrix from g to f size:*/ + if (debug) _printf_(" reducing matrix from g to f set\n"); + Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets); + + /*!Reduce load from g to f size: */ + if (debug) _printf_(" reducing load from g to f set\n"); + Reduceloadfromgtofx(&pf, pg, fem->Gmn, Kfs, fem->ys, fem->nodesets); + + /*Solve: */ + if (debug) _printf_(" solving\n"); + Solverx(&uf, Kff, pf, NULL, solver_string); + + //Merge back to g set + if (debug) _printf_(" merging solution from f to g set\n"); + Mergesolutionfromftogx(&ug, uf,fem->Gmn,fem->ys,fem->nodesets); + + /*Assign output pointers:*/ + *pug=ug; +} Index: /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp (revision 433) +++ /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp (revision 434) @@ -11,5 +11,5 @@ #include "../issm.h" -void diagnostic_core_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, ParameterInputs* inputs,FemModel* fem){ +void diagnostic_core_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, FemModel* fem,ParameterInputs* inputs,int analysis_type){ @@ -38,5 +38,4 @@ /*parameters:*/ int kflag,pflag,connectivity,numberofdofspernode; - int analysis_type; char* solver_string=NULL; int debug=0; @@ -55,5 +54,4 @@ fem->parameters->FindParam((void*)&yts,"yts"); fem->parameters->FindParam((void*)&debug,"debug"); - fem->parameters->FindParam((void*)&analysis_type,"analysis_type"); /*Copy loads for backup: */ Index: /issm/trunk/src/c/parallel/objectivefunctionC.cpp =================================================================== --- /issm/trunk/src/c/parallel/objectivefunctionC.cpp (revision 433) +++ /issm/trunk/src/c/parallel/objectivefunctionC.cpp (revision 434) @@ -75,5 +75,5 @@ //Run diagnostic with updated parameters. - diagnostic_core_nonlinear(&u_g,NULL,NULL,inputs,femmodel); + diagnostic_core_nonlinear(&u_g,NULL,NULL,femmodel,inputs,analysis_type); VecToMPISerial(&u_g_double,u_g); VecFree(&u_g); Index: /issm/trunk/src/c/parallel/parallel.h =================================================================== --- /issm/trunk/src/c/parallel/parallel.h (revision 433) +++ /issm/trunk/src/c/parallel/parallel.h (revision 434) @@ -7,8 +7,11 @@ #include "../objects/objects.h" +#include "../io/io.h" Vec GradJCompute(ParameterInputs* inputs,FemModel* femmodel,double* u_g_obs); -void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, ParameterInputs* inputs,FemModel* fem); +void diagnostic_core(Vec* pug, Vec* ppg,FemModel* fems, ParameterInputs* inputs); +void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* fem,ParameterInputs* inputs,int analysis_type); +void diagnostic_core_linear(Vec* ppg,FemModel* fem,ParameterInputs* inputs,int analysis_type); int cielothermal_core(Vec** pt_g,ParameterInputs* inputs,FemModel* femmodel); @@ -25,5 +28,5 @@ //int ParameterUpdate(double* search_vector,int step, WorkspaceParams* workspaceparams,BatchParams* batchparams); -void OutputDiagnostic(Vec u_g,Vec tpartition,char* filename,NodeSets* nodesets,char* analysis_type); +void OutputDiagnostic(Vec u_g,Vec p_g, FemModel* femmodels,char* filename); int OutputThermal(Vec* t_g,Vec* tpartition,char* filename,char* analysis_type); void OutputControl(Vec u_g,double* p_g, double* J, int nsteps, Vec partition,char* filename,NodeSets* nodesets); Index: /issm/trunk/src/c/shared/Dofs/dofs.h =================================================================== --- /issm/trunk/src/c/shared/Dofs/dofs.h (revision 434) +++ /issm/trunk/src/c/shared/Dofs/dofs.h (revision 434) @@ -0,0 +1,12 @@ +/*!\file: dofs.h + * \brief prototypes for dofs.h + */ + +#ifndef _DOFS_H_ +#define _DOFS_H_ + +double* dofsetgen(int numdofs,int* doflist,int dofspernode,int totaldofs); + + +#endif //ifndef _DOFS_H_ + Index: /issm/trunk/src/c/shared/Dofs/dofsetgen.cpp =================================================================== --- /issm/trunk/src/c/shared/Dofs/dofsetgen.cpp (revision 434) +++ /issm/trunk/src/c/shared/Dofs/dofsetgen.cpp (revision 434) @@ -0,0 +1,47 @@ +/*!\file: dofsetgen.cpp + * \brief create list of dofs. + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../Alloc/alloc.h" +#include "../../include/macros.h" +#include "stdio.h" + +#undef __FUNCT__ +#define __FUNCT__ "dofsetgen" + +double* dofsetgen(int numdofs,int* doflist,int dofspernode,int totaldofs){ + + int i,j; + int count=0; + + /*output: */ + double* outdoflist=NULL; + + /*number of nodes :*/ + int numberofnodes; + + /*allocate: */ + numberofnodes=(int)totaldofs/dofspernode; + outdoflist=(double*)xmalloc(numberofnodes*numdofs*sizeof(double)); + + count=0; + for(i=0;i