Index: /issm/trunk/src/c/DataSet/DataSet.cpp =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.cpp (revision 4054) +++ /issm/trunk/src/c/DataSet/DataSet.cpp (revision 4055) @@ -684,29 +684,4 @@ xfree((void**)&truedofs); xfree((void**)&alltruedofs); - -} -/*}}}*/ -/*FUNCTION DataSet::FieldExtrude{{{1*/ -void DataSet::FieldExtrude(Vec field,double* field_serial,char* field_name, int collapse){ - - vector::iterator object; - Penta* penta=NULL; - Node* node=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - if((*object)->Enum()==PentaEnum){ - - penta=(Penta*)(*object); - penta->FieldExtrude(field,field_serial,field_name,collapse); - - } - if((*object)->Enum()==NodeEnum){ - - node=(Node*)(*object); - node->FieldExtrude(field,field_serial,field_name); - - } - } } Index: /issm/trunk/src/c/DataSet/DataSet.h =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.h (revision 4054) +++ /issm/trunk/src/c/DataSet/DataSet.h (revision 4055) @@ -80,5 +80,4 @@ void InputExtrude(int enum_type); int DeleteObject(Object* object); - void FieldExtrude(Vec field,double* field_serial,char* field_name, int collapse); void OutputRifts(Vec riftproperties); /*}}}*/ Index: /issm/trunk/src/c/EnumDefinitions/EnumAsString.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumAsString.cpp (revision 4054) +++ /issm/trunk/src/c/EnumDefinitions/EnumAsString.cpp (revision 4055) @@ -26,7 +26,4 @@ case VerticesEnum : return "Vertices"; case SolutionTypeEnum : return "SolutionType"; - case DiagnosticSolutionEnum : return "DiagnosticSolution"; - case SurfaceSlopeSolutionEnum : return "SurfaceSlopeSolution"; - case BedSlopeSolutionEnum : return "BedSlopeSolution"; case AnalysisTypeEnum : return "AnalysisType"; case SubAnalysisTypeEnum : return "SubAnalysisType"; @@ -46,7 +43,8 @@ case InverseAnalysisEnum : return "InverseAnalysis"; case ThermalAnalysisEnum : return "ThermalAnalysis"; - case TransientAnalysisEnum : return "TransientAnalysis"; + case Transient2DAnalysisEnum : return "Transient2DAnalysis"; + case Transient3DAnalysisEnum : return "Transient3DAnalysis"; case SteadyAnalysisEnum : return "SteadyAnalysis"; - case SlopeComputeAnalysisEnum : return "SlopeComputeAnalysis"; + case SlopeAnalysisEnum : return "SlopeAnalysis"; case BedSlopeXAnalysisEnum : return "BedSlopeXAnalysis"; case BedSlopeYAnalysisEnum : return "BedSlopeYAnalysis"; @@ -220,4 +218,7 @@ case StringExternalResultEnum : return "StringExternalResult"; case JEnum : return "J"; + case RelativeEnum : return "Relative"; + case ResidualEnum : return "Residual"; + case AbsoluteEnum : return "Absolute"; case BetaEnum : return "Beta"; case CmGradientEnum : return "CmGradient"; @@ -225,5 +226,7 @@ case CmMaxEnum : return "CmMax"; case CmMinEnum : return "CmMin"; + case AdjointEnum : return "Adjoint"; case GradientEnum : return "Gradient"; + case OldGradientEnum : return "OldGradient"; case ConnectivityEnum : return "Connectivity"; case ControlParameterEnum : return "ControlParameter"; Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 4054) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 4055) @@ -22,7 +22,4 @@ /*Solution types {{{1 */ SolutionTypeEnum, - DiagnosticSolutionEnum, - SurfaceSlopeSolutionEnum, - BedSlopeSolutionEnum, /*}}}*/ /*Analysis types {{{1 */ @@ -48,8 +45,9 @@ ThermalAnalysisEnum, //transient - TransientAnalysisEnum, + Transient2DAnalysisEnum, + Transient3DAnalysisEnum, SteadyAnalysisEnum, //slope - SlopeComputeAnalysisEnum, + SlopeAnalysisEnum, BedSlopeXAnalysisEnum, BedSlopeYAnalysisEnum, @@ -254,4 +252,9 @@ JEnum, /*}}}*/ + /*Convergence{{{1*/ + RelativeEnum, + ResidualEnum, + AbsoluteEnum, + /*}}}*/ /*Parameters{{{1*/ BetaEnum, @@ -260,5 +263,7 @@ CmMaxEnum, CmMinEnum, + AdjointEnum, GradientEnum, + OldGradientEnum, ConnectivityEnum, ControlParameterEnum, Index: /issm/trunk/src/c/EnumDefinitions/StringAsEnum.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/StringAsEnum.cpp (revision 4054) +++ /issm/trunk/src/c/EnumDefinitions/StringAsEnum.cpp (revision 4055) @@ -24,7 +24,4 @@ else if (strcmp(name,"Vertices")==0) return VerticesEnum; else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum; - else if (strcmp(name,"DiagnosticSolution")==0) return DiagnosticSolutionEnum; - else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum; - else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum; else if (strcmp(name,"AnalysisType")==0) return AnalysisTypeEnum; else if (strcmp(name,"SubAnalysisType")==0) return SubAnalysisTypeEnum; @@ -44,7 +41,8 @@ else if (strcmp(name,"InverseAnalysis")==0) return InverseAnalysisEnum; else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum; - else if (strcmp(name,"TransientAnalysis")==0) return TransientAnalysisEnum; + else if (strcmp(name,"Transient2DAnalysis")==0) return Transient2DAnalysisEnum; + else if (strcmp(name,"Transient3DAnalysis")==0) return Transient3DAnalysisEnum; else if (strcmp(name,"SteadyAnalysis")==0) return SteadyAnalysisEnum; - else if (strcmp(name,"SlopeComputeAnalysis")==0) return SlopeComputeAnalysisEnum; + else if (strcmp(name,"SlopeAnalysis")==0) return SlopeAnalysisEnum; else if (strcmp(name,"BedSlopeXAnalysis")==0) return BedSlopeXAnalysisEnum; else if (strcmp(name,"BedSlopeYAnalysis")==0) return BedSlopeYAnalysisEnum; @@ -218,4 +216,7 @@ else if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum; else if (strcmp(name,"J")==0) return JEnum; + else if (strcmp(name,"Relative")==0) return RelativeEnum; + else if (strcmp(name,"Residual")==0) return ResidualEnum; + else if (strcmp(name,"Absolute")==0) return AbsoluteEnum; else if (strcmp(name,"Beta")==0) return BetaEnum; else if (strcmp(name,"CmGradient")==0) return CmGradientEnum; @@ -223,5 +224,7 @@ else if (strcmp(name,"CmMax")==0) return CmMaxEnum; else if (strcmp(name,"CmMin")==0) return CmMinEnum; + else if (strcmp(name,"Adjoint")==0) return AdjointEnum; else if (strcmp(name,"Gradient")==0) return GradientEnum; + else if (strcmp(name,"OldGradient")==0) return OldGradientEnum; else if (strcmp(name,"Connectivity")==0) return ConnectivityEnum; else if (strcmp(name,"ControlParameter")==0) return ControlParameterEnum; Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 4054) +++ /issm/trunk/src/c/Makefile.am (revision 4055) @@ -510,8 +510,10 @@ ./modules/DepthAverageInputx/DepthAverageInputx.cpp\ ./modules/DepthAverageInputx/DepthAverageInputx.h\ - ./modules/FieldExtrudex/FieldExtrudex.cpp\ - ./modules/FieldExtrudex/FieldExtrudex.h\ + ./modules/VecExtrudex/VecExtrudex.cpp\ + ./modules/VecExtrudex/VecExtrudex.h\ ./modules/InputToResultx/InputToResultx.cpp\ ./modules/InputToResultx/InputToResultx.h\ + ./modules/InputConvergencex/InputConvergencex.cpp\ + ./modules/InputConvergencex/InputConvergencex.h\ ./modules/ExtrudeInputx/ExtrudeInputx.cpp\ ./modules/ExtrudeInputx/ExtrudeInputx.h\ @@ -1014,8 +1016,10 @@ ./modules/DepthAverageInputx/DepthAverageInputx.cpp\ ./modules/DepthAverageInputx/DepthAverageInputx.h\ - ./modules/FieldExtrudex/FieldExtrudex.cpp\ - ./modules/FieldExtrudex/FieldExtrudex.h\ + ./modules/VecExtrudex/VecExtrudex.cpp\ + ./modules/VecExtrudex/VecExtrudex.h\ ./modules/InputToResultx/InputToResultx.cpp\ ./modules/InputToResultx/InputToResultx.h\ + ./modules/InputConvergencex/InputConvergencex.cpp\ + ./modules/InputConvergencex/InputConvergencex.h\ ./modules/ExtrudeInputx/ExtrudeInputx.cpp\ ./modules/ExtrudeInputx/ExtrudeInputx.h\ @@ -1030,5 +1034,6 @@ ./solutions/controlrestart.cpp\ ./solutions/objectivefunctionC.cpp\ - ./solutions/gradjcompute_core.cpp\ + ./solutions/gradient_core.cpp\ + ./solutions/adjoint_core.cpp\ ./solutions/prognostic_core.cpp\ ./solutions/prognostic2_core.cpp\ @@ -1040,7 +1045,6 @@ ./solutions/bedslope.cpp\ ./solutions/bedslope_core.cpp\ - ./solutions/transient_core.cpp\ - ./solutions/transient_core_2d.cpp\ - ./solutions/transient_core_3d.cpp\ + ./solutions/transient2d_core.cpp\ + ./solutions/transient3d_core.cpp\ ./solutions/steadystate_core.cpp\ ./solutions/ResetBoundaryConditions.cpp\ @@ -1058,5 +1062,5 @@ bin_PROGRAMS = else -bin_PROGRAMS = DiagnosticSolution.exe ThermalSolution.exe PrognosticSolution.exe Prognostic2Solution.exe BalancedthicknessSolution.exe Balancedthickness2Solution.exe BalancedvelocitiesSolution.exe TransientSolution.exe SteadystateSolution.exe SurfaceSlopeSolution.exe BedSlopeSolution.exe +bin_PROGRAMS = DiagnosticSolution.exe ThermalSolution.exe PrognosticSolution.exe Prognostic2Solution.exe BalancedthicknessSolution.exe Balancedthickness2Solution.exe BalancedvelocitiesSolution.exe Transient2DSolution.exe Transient23Solution.exe SteadystateSolution.exe SurfaceSlopeSolution.exe BedSlopeSolution.exe endif @@ -1093,4 +1097,7 @@ BedSlopeSolution_exe_CXXFLAGS= -fPIC -D_PARALLEL_ -TransientSolution_exe_SOURCES = solutions/transient.cpp -TransientSolution_exe_CXXFLAGS= -fPIC -D_PARALLEL_ +Transient2DSolution_exe_SOURCES = solutions/transient2d.cpp +Transient2DSolution_exe_CXXFLAGS= -fPIC -D_PARALLEL_ + +Transient3DSolution_exe_SOURCES = solutions/transient3d.cpp +Transient3DSolution_exe_CXXFLAGS= -fPIC -D_PARALLEL_ Index: /issm/trunk/src/c/modules/CostFunctionx/CostFunctionx.h =================================================================== --- /issm/trunk/src/c/modules/CostFunctionx/CostFunctionx.h (revision 4054) +++ /issm/trunk/src/c/modules/CostFunctionx/CostFunctionx.h (revision 4055) @@ -10,6 +10,5 @@ /* local prototypes: */ -void CostFunctionx( double* pJ, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters, - int analysis_type,int sub_analysis_type); +void CostFunctionx( double* pJ, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters); #endif /* _MISFITX_H */ Index: /issm/trunk/src/c/modules/DepthAverageInputx/DepthAverageInputx.cpp =================================================================== --- /issm/trunk/src/c/modules/DepthAverageInputx/DepthAverageInputx.cpp (revision 4054) +++ /issm/trunk/src/c/modules/DepthAverageInputx/DepthAverageInputx.cpp (revision 4055) @@ -10,5 +10,5 @@ #include "../../EnumDefinitions/EnumDefinitions.h" -void DepthAverageInputx( DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int enum_type){ +void DepthAverageInputx( DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int enum_type,int average_enum_type){ /*Intermediary*/ @@ -22,9 +22,8 @@ for (i=0;iSize();i++){ element=(Element*)elements->GetObjectByOffset(i); - element->DepthAverageInputAtBase(enum_type); + element->DepthAverageInputAtBase(enum_type,average_enum_type); } /*Then extrude vertically the new inputs*/ - elements->InputExtrude(enum_type); - + elements->InputExtrude(average_enum_type); } Index: /issm/trunk/src/c/modules/DepthAverageInputx/DepthAverageInputx.h =================================================================== --- /issm/trunk/src/c/modules/DepthAverageInputx/DepthAverageInputx.h (revision 4054) +++ /issm/trunk/src/c/modules/DepthAverageInputx/DepthAverageInputx.h (revision 4055) @@ -9,5 +9,5 @@ /* local prototypes: */ -void DepthAverageInputx( DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int enum_type); +void DepthAverageInputx( DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int enum_type,int average_enum_type); #endif /* _DEPTHAVERAGEINPUTX_H */ Index: /issm/trunk/src/c/modules/Dux/Dux.h =================================================================== --- /issm/trunk/src/c/modules/Dux/Dux.h (revision 4054) +++ /issm/trunk/src/c/modules/Dux/Dux.h (revision 4055) @@ -10,7 +10,5 @@ /* local prototypes: */ -void Dux( Vec* pdu_g, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters, - int analysis_type,int sub_analysis_type); - +void Dux( Vec* pdu_g, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters); #endif /* _DUX_H */ Index: /issm/trunk/src/c/modules/FieldDepthAveragex/FieldDepthAveragex.cpp =================================================================== --- /issm/trunk/src/c/modules/FieldDepthAveragex/FieldDepthAveragex.cpp (revision 4054) +++ /issm/trunk/src/c/modules/FieldDepthAveragex/FieldDepthAveragex.cpp (revision 4055) @@ -43,5 +43,8 @@ /*Extrude field to all nodes: */ - nodes->FieldExtrude(field,field_serial,fieldname,0); + for(i=0;iSize();i++){ + node=(Node*)nodes->GetObjectByOffset(i); + node->VecExtrude(field,field_serial); + } /*Assemble vector: */ Index: /issm/trunk/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.h =================================================================== --- /issm/trunk/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.h (revision 4054) +++ /issm/trunk/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.h (revision 4055) @@ -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); #endif /* _GETSOLUTIONFROMINPUTSXX_H */ Index: /issm/trunk/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp =================================================================== --- /issm/trunk/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp (revision 4054) +++ /issm/trunk/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp (revision 4055) @@ -38,2 +38,20 @@ } + +void GetVectorFromInputsx( double* pvector, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters, int NameEnum, int TypeEnum){ + + /*output: */ + double* vector=NULL; + + /*intermediary: */ + Vec vec_vector=NULL; + + GetVectorFromInputsx( &vec_vector, elements,nodes, vertices, loads, materials, parameters, NameEnum, TypeEnum); + VecToMPISerial(&vector,vec_vector); + + /*Free ressources:*/ + VecFree(&vec_vector); + + /*Assign output pointers:*/ + *pvector=vector; +} Index: /issm/trunk/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.h =================================================================== --- /issm/trunk/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.h (revision 4054) +++ /issm/trunk/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.h (revision 4055) @@ -10,4 +10,5 @@ /* local prototypes: */ void GetVectorFromInputsx( Vec* pvector, DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials, Parameters* parameters,int NameEnum,int TypeEnum); +void GetVectorFromInputsx( double* pvector, DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials, Parameters* parameters,int NameEnum,int TypeEnum); #endif /* _GETVECTORFROMINPUTSXX_H */ Index: /issm/trunk/src/c/modules/Gradjx/Gradjx.cpp =================================================================== --- /issm/trunk/src/c/modules/Gradjx/Gradjx.cpp (revision 4054) +++ /issm/trunk/src/c/modules/Gradjx/Gradjx.cpp (revision 4055) @@ -10,10 +10,11 @@ #include "../../EnumDefinitions/EnumDefinitions.h" -void Gradjx( DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int control_type){ +void Gradjx( Vec* pgradient, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int control_type){ int i; int dim; - - + int numberofvertices; + Vec gradient=NULL; + /*First, get elements and loads configured: */ elements-> Configure(elements,loads, nodes,vertices, materials,parameters); @@ -23,13 +24,20 @@ /*retrieve some parameters: */ parameters->FindParam(&dim,DimEnum); + parameters->FindParam(&numberofvertices,NumberOfVerticesEnum); /*Compute gradients: */ for (i=0;iSize();i++){ Element* element=(Element*)elements->GetObjectByOffset(i); - element->Gradj(control_type); + element->Gradj(gradient,control_type); } + /*Assemble vector: */ + VecAssemblyBegin(gradient); + VecAssemblyEnd(gradient); + /*Extrude if needed: */ - if(dim==3) ExtrudeInputx( elements,nodes, vertices, loads, materials, parameters,GradientEnum); + if(dim==3) VecExtrudex(gradient, elements,nodes, vertices, loads, materials, parameters,0); + /*Assign output pointers: */ + *pgradient=gradient; } Index: /issm/trunk/src/c/modules/Gradjx/Gradjx.h =================================================================== --- /issm/trunk/src/c/modules/Gradjx/Gradjx.h (revision 4054) +++ /issm/trunk/src/c/modules/Gradjx/Gradjx.h (revision 4055) @@ -10,5 +10,5 @@ /* local prototypes: */ -void Gradjx( DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters, int control_type); +void Gradjx(Vec* pgrad_g, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters, int control_type); #endif /* _GRADJX_H */ Index: /issm/trunk/src/c/modules/InputConvergencex/InputConvergencex.cpp =================================================================== --- /issm/trunk/src/c/modules/InputConvergencex/InputConvergencex.cpp (revision 4055) +++ /issm/trunk/src/c/modules/InputConvergencex/InputConvergencex.cpp (revision 4055) @@ -0,0 +1,58 @@ +/*!\file InputConvergencex + */ + +#include "./InputConvergencex.h" +#include "../../shared/shared.h" +#include "../../include/include.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" + +int InputConvergencex(DataSet* elements,DataSet* nodes,DataSet* vertices,DataSet* loads,DataSet* materials,Parameters* parameters,int enums, int num_enums, int criterionenums, double* criterionvalues,int num_criterionenums){ + + /*intermediary:*/ + int i; + int converged = 1; + int node_converged=1; + double *eps = NULL; + extern int my_rank; + int rank_allowed; + + /*First, configure our datasets: */ + elements->Configure(elements,loads, nodes,vertices, materials,parameters); + nodes-> Configure(elements,loads, nodes,vertices, materials,parameters); + loads-> Configure(elements, loads, nodes,vertices, materials,parameters); + parameters->Configure(elements,loads, nodes,vertices, materials,parameters); + + /*allocate dynamic memory: */ + eps=(double*)xmalloc(num_criterionenums*sizeof(double)); + + /*Go through elements, and ask them to do the job: */ + for(i=0;iSize();i++){ + Element* element=(Element*)elements->GetObjectByOffset(i); + element->InputConvergence(&converged,eps,enums,num_enums,criterionenums,criterionvalues,num_criterionenums); + if(!converged)break; + } + + if(!converged){ + /*get convergence information to user: */ + if(verbose>2)rank_allowed=my_rank; + else rank_allowed=0; + if(my_rank==rank_allowed){ + printf("%s%i%s"," non convergence on element ",element->Id()," for inputs ("); for(i=0;idim==2) parameters->AddObject(new IntParam(DimEnum,2)); else parameters->AddObject(new IntParam(DimEnum,3)); - parameters->AddObject(new StringParam(OutputFileNameEnum,iomodel->outputfilename)); parameters->AddObject(new BoolParam(IsHutterEnum,iomodel->ishutter)); parameters->AddObject(new BoolParam(IsMacAyealPattynEnum,iomodel->ismacayealpattyn)); Index: /issm/trunk/src/c/modules/ModelProcessorx/SlopeCompute/CreateNodesSlopeCompute.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/SlopeCompute/CreateNodesSlopeCompute.cpp (revision 4054) +++ /issm/trunk/src/c/modules/ModelProcessorx/SlopeCompute/CreateNodesSlopeCompute.cpp (revision 4055) @@ -46,5 +46,5 @@ /*Add node to nodes dataset: */ - nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,iomodel,SlopeComputeAnalysisEnum)); + nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,iomodel,SlopeAnalysisEnum)); } Index: /issm/trunk/src/c/modules/OutputResultsx/OutputResultsx.cpp =================================================================== --- /issm/trunk/src/c/modules/OutputResultsx/OutputResultsx.cpp (revision 4054) +++ /issm/trunk/src/c/modules/OutputResultsx/OutputResultsx.cpp (revision 4055) @@ -15,9 +15,10 @@ #include "../../objects/objects.h" -void OutputResults(DataSet* elements, DataSet* loads, DataSet* nodes, DataSet* vertices, DataSet* materials, Parameters* parameters, char* filename){ +void OutputResults(DataSet* elements, DataSet* loads, DataSet* nodes, DataSet* vertices, DataSet* materials, Parameters* parameters){ int i; Patch* patch=NULL; + char* filename=NULL; int solutiontype; @@ -37,4 +38,5 @@ //Recover solutiontype which will be output to disk: parameters->FindParam(&solutiontype,SolutionTypeEnum); + parameters->FindParam(&filename,OutputFileNameEnum); //Process results to be output in the correct units @@ -104,4 +106,5 @@ /*Free ressources:*/ delete patch; + xfree((void**)&filename); } Index: /issm/trunk/src/c/modules/OutputResultsx/OutputResultsx.h =================================================================== --- /issm/trunk/src/c/modules/OutputResultsx/OutputResultsx.h (revision 4054) +++ /issm/trunk/src/c/modules/OutputResultsx/OutputResultsx.h (revision 4055) @@ -9,5 +9,5 @@ /* local prototypes: */ -void OutputResults(DataSet* elements, DataSet* loads, DataSet* nodes, DataSet* vertices, DataSet* materials, Parameters* parameters, char* filename); +void OutputResults(DataSet* elements, DataSet* loads, DataSet* nodes, DataSet* vertices, DataSet* materials, Parameters* parameters); #endif /* _OUTPUTRESULTS_H */ Index: /issm/trunk/src/c/modules/UpdateGeometryx/UpdateGeometryx.cpp =================================================================== --- /issm/trunk/src/c/modules/UpdateGeometryx/UpdateGeometryx.cpp (revision 4054) +++ /issm/trunk/src/c/modules/UpdateGeometryx/UpdateGeometryx.cpp (revision 4055) @@ -11,12 +11,13 @@ #include "../../DataSet/DataSet.h" -void UpdateGeometryx(Vec* poutthickness,Vec* poutbed,Vec* poutsurface, - DataSet* elements, DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials, Parameters* parameters, - Vec vec_newthickness,Vec vec_bed,Vec vec_surface){ +void UpdateGeometryx( DataSet* elements, DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials, Parameters* parameters){ int i; int dof; - /*serialized input: */ + /*vectorized inputs: */ + Vec vec_newthickness=NULL; + Vec vec_bed=NULL; + Vec vec_surface=NULL: double* newthickness=NULL; double* bed=NULL; @@ -24,5 +25,5 @@ /*objects: */ - Node* node=NULL; + Vertex* vertex=NULL; Object* object=NULL; Matpar* matpar=NULL; @@ -31,14 +32,4 @@ double h,b,s; double rho_ice,rho_water; - - /*output: */ - Vec outthickness=NULL; - Vec outbed=NULL; - Vec outsurface=NULL; - - /*Duplicate inputs to outputs, do not copy values!: */ - VecDuplicate(vec_newthickness,&outthickness); - VecDuplicate(vec_bed,&outbed); - VecDuplicate(vec_surface,&outsurface); /*First, get elements and loads configured: */ @@ -49,5 +40,9 @@ parameters->Configure(elements,loads, nodes,vertices, materials,parameters); - /*Serialize inputs :*/ + /*retrieve inputs: */ + GetVectorFromInputsx( &vec_newthickness, elements, nodes, vertices, loads, materials, parameters, ThicknessEnum,VertexEnum); + GetVectorFromInputsx( &vec_bed, elements, nodes, vertices, loads, materials, parameters, ThicknessEnum,VertexEnum); + GetVectorFromInputsx( &vec_surface, elements, nodes, vertices, loads, materials, parameters, ThicknessEnum,VertexEnum); + VecToMPISerial(&newthickness,vec_newthickness); VecToMPISerial(&bed,vec_bed); @@ -67,14 +62,14 @@ rho_water=matpar->GetRhoWater(); - /*Go through nodes and for each node, update the thickness, bed and surface, using the + /*Go through vertices and for each vertex, update the thickness, bed and surface, using the * new thickness computed in the prognostic_core part of the transient solution sequence: */ - for(i=0;iSize();i++){ + for(i=0;iSize();i++){ /*get i'th node: */ - node=(Node*)nodes->GetObjectByOffset(i); + vertex=(Vertex*)vertices->GetObjectByOffset(i); - /*first, recover thickness, surface and bed for this node: */ - dof=node->GetDofList1(); + /*first, recover thickness, surface and bed for this vertex: */ + dof=vertex->GetDofList1(); h=newthickness[dof]; s=surface[dof]; @@ -85,10 +80,10 @@ //For grids on ice sheet, the new bed remains the same, the new surface becomes bed+new_thickness. - if (node->IsOnSheet()){ + if (vertex->IsOnSheet()){ s=b+h; } //For grids on ice shelt, we have hydrostatic equilibrium (for now) - if (node->IsOnShelf()){ + if (vertex->IsOnShelf()){ b=-rho_ice/rho_water*h; s=(1-rho_ice/rho_water)*h; @@ -96,27 +91,31 @@ /*Ok, plug values of thickness, surafce and bed into our outputs: */ - VecSetValues(outthickness,1,&dof,&h,INSERT_VALUES); - VecSetValues(outbed,1,&dof,&b,INSERT_VALUES); - VecSetValues(outsurface,1,&dof,&s,INSERT_VALUES); + VecSetValues(vec_newthickness,1,&dof,&h,INSERT_VALUES); + VecSetValues(vec_bed,1,&dof,&b,INSERT_VALUES); + VecSetValues(vec_surface,1,&dof,&s,INSERT_VALUES); } /*Assemble vectors: */ - VecAssemblyBegin(outthickness); - VecAssemblyEnd(outthickness); + VecAssemblyBegin(vec_newthickness); + VecAssemblyEnd(vec_newthickness); - VecAssemblyBegin(outbed); - VecAssemblyEnd(outbed); + VecAssemblyBegin(vec_bed); + VecAssemblyEnd(vec_bed); - VecAssemblyBegin(outsurface); - VecAssemblyEnd(outsurface); + VecAssemblyBegin(vec_surface); + VecAssemblyEnd(vec_surface); + + /*Put back into inputs:*/ + UpdateInputsFromVectorx( elements,nodes, vertices,loads, materials, parameters,vec_newthickness,ThicknessEnum,VertexEnum); //clobber existing thickness input with new one + UpdateInputsFromVectorx( elements,nodes, vertices,loads, materials, parameters,vec_bed,BedEnum,VertexEnum); + UpdateInputsFromVectorx( elements,nodes, vertices,loads, materials, parameters,vec_surface,SurfaceEnum,VertexEnum); /*Free ressources:*/ + VecFree(&vec_newthickness); + VecFree(&vec_bed); + VecFree(&vec_surface); xfree((void**)&newthickness); xfree((void**)&bed); xfree((void**)&surface); - /*Assign output pointers: */ - *poutthickness=outthickness; - *poutbed=outbed; - *poutsurface=outsurface; } Index: /issm/trunk/src/c/modules/UpdateGeometryx/UpdateGeometryx.h =================================================================== --- /issm/trunk/src/c/modules/UpdateGeometryx/UpdateGeometryx.h (revision 4054) +++ /issm/trunk/src/c/modules/UpdateGeometryx/UpdateGeometryx.h (revision 4055) @@ -10,7 +10,5 @@ /* local prototypes: */ -void UpdateGeometryx(Vec* poutthickness,Vec* poutbed,Vec* poutsurface, - DataSet* elements, DataSet* nodes,DataSet* vertices,DataSet* loads, DataSet* materials, Parameters* parameters, - Vec newthickness,Vec bed,Vec surface); +void UpdateGeometryx(DataSet* elements, DataSet* nodes,DataSet* vertices,DataSet* loads, DataSet* materials, Parameters* parameters); #endif /* _UPDATEGEOMETRYX_H */ Index: /issm/trunk/src/c/modules/VecExtrudex/VecExtrudex.cpp =================================================================== --- /issm/trunk/src/c/modules/VecExtrudex/VecExtrudex.cpp (revision 4055) +++ /issm/trunk/src/c/modules/VecExtrudex/VecExtrudex.cpp (revision 4055) @@ -0,0 +1,56 @@ +/*!\file VecExtrudex + * \brief: vertical velocity extrusion + */ + +#include "./VecExtrudex.h" + +#include "../../shared/shared.h" +#include "../../include/include.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" + +void VecExtrudex( Vec vector, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int collapse){ + + double* vector_serial=NULL; + int i; + Element* element=NULL; + Penta* penta=NULL; + Node* node=NULL; + + /*First, get elements and nodes configured: */ + elements->Configure(elements,loads, nodes,vertices, materials,parameters); + nodes->Configure(elements,loads, nodes,vertices, materials,parameters); + parameters->Configure(elements,loads, nodes,vertices, materials,parameters); + + /*Serialize vector: */ + VecToMPISerial(&vector_serial,vector); + + /*ok, is the extrusion dependent upon elements or just nodes? If collapse is 1, it means we only want + * to extrude the vector on elemenets that are collapsed (element operation). If collapse is 0, it means + * we want to extrude vector on every node (node operation): */ + + if(collapse){ + /*Extrude vector vertically: */ + for(i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + if(element->Enum()==PentaEnum){ + penta=(Penta*)element; + penta->VecExtrude(vector,vector_serial,collapse); + } + } + } + else{ + for(i=0;iSize();i++){ + node=(Node*)nodes->GetObjectByOffset(i); + node->VecExtrude(vector,vector_serial); + } + } + + /*Assemble vector: */ + VecAssemblyBegin(vector); + VecAssemblyEnd(vector); + + /*Free ressources:*/ + xfree((void**)&vector_serial); + +} Index: /issm/trunk/src/c/modules/VecExtrudex/VecExtrudex.h =================================================================== --- /issm/trunk/src/c/modules/VecExtrudex/VecExtrudex.h (revision 4055) +++ /issm/trunk/src/c/modules/VecExtrudex/VecExtrudex.h (revision 4055) @@ -0,0 +1,14 @@ +/*!\file: VecExtrudex.h + * \brief header file for field extrusion + */ + +#ifndef _VECEXTRUDEX_H +#define _VECEXTRUDEX_H + +#include "../../DataSet/DataSet.h" + +/* local prototypes: */ +void VecExtrudex( Vec vector, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int collapse); + +#endif /* _VECEXTRUDEX_H */ + Index: /issm/trunk/src/c/modules/modules.h =================================================================== --- /issm/trunk/src/c/modules/modules.h (revision 4054) +++ /issm/trunk/src/c/modules/modules.h (revision 4055) @@ -53,5 +53,5 @@ #include "./ComputePressurex/ComputePressurex.h" #include "./ComputeStrainRatex/ComputeStrainRatex.h" -#include "./FieldExtrudex/FieldExtrudex.h" +#include "./VecExtrudex/VecExtrudex.h" #include "./Qmux/Qmux.h" #include "./NodeConnectivityx/NodeConnectivityx.h" @@ -81,4 +81,6 @@ #include "./ScaleInputx/ScaleInputx.h" #include "./AXPYInputx/AXPYInputx.h" +#include "./GetVectorFromInputsx/GetVectorFromInputsx.h" +#include "./InputConvergencex/InputConvergencex.h" #endif Index: /issm/trunk/src/c/objects/Elements/Beam.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Beam.cpp (revision 4054) +++ /issm/trunk/src/c/objects/Elements/Beam.cpp (revision 4055) @@ -562,15 +562,15 @@ /*}}}*/ /*FUNCTION Beam::Gradj{{{1*/ -void Beam::Gradj(int control_type){ +void Beam::Gradj(Vec gradient,int control_type){ ISSMERROR(" not supported yet!"); } /*}}}*/ /*FUNCTION Beam::GradjB{{{1*/ -void Beam::GradjB(void){ +void Beam::GradjB(Vec gradient){ ISSMERROR(" not supported yet!"); } /*}}}*/ /*FUNCTION Beam::GradjDrag{{{1*/ -void Beam::GradjDrag(void){ +void Beam::GradjDrag(Vec gradient){ ISSMERROR(" not supported yet!"); } @@ -1055,2 +1055,30 @@ } /*}}}*/ +/*FUNCTION Beam::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/ +void Beam::InputConvergence(int* pconverged,double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){a + + Input** new_inputs=NULL; + Input** old_inputs=NULL; + int converged=1; + + new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs + old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs + + for(i=0;iinputs->GetInput(new_enums[2*i+0]); + old_inputs[i]=(Input*)this->inputs->GetInput(old_enums[2*i+1]); + if(!new_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(new_enums[2*i+0])); + if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(old_enums[2*i+0])); + } + + /*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/ + for(i=0;icriterionvalues[i]) converged=0; + } + + /*Assign output pointers:*/ + *pconverged=converged; + +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Beam.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Beam.h (revision 4054) +++ /issm/trunk/src/c/objects/Elements/Beam.h (revision 4055) @@ -62,5 +62,5 @@ void UpdateInputsFromConstant(bool constant, int name){ISSMERROR("Not implemented yet!");} void UpdateInputsFromSolution(double* solution); - void DepthAverageInputAtBase(int enum_type){ISSMERROR("not implemented yet");}; + void DepthAverageInputAtBase(int enum_type,int average_enum_type){ISSMERROR("not implemented yet");}; void InputToResult(int enum_type,int step,double time); void ProcessResultsUnits(void); @@ -98,4 +98,5 @@ void AXPYInput(int YEnum, double scalar, int XEnum); void ControlConstrainInput(int control_type,double cm_min, double cm_max); + void InputConvergence(int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums); /*}}}*/ @@ -106,7 +107,7 @@ void GetThicknessList(double* thickness_list); void Du(Vec); - void Gradj(int control_type); - void GradjDrag(void); - void GradjB(void); + void Gradj(Vec gradient,int control_type); + void GradjDrag(Vec gradient); + void GradjB(Vec gradient); double Misfit(void); double SurfaceArea(void); Index: /issm/trunk/src/c/objects/Elements/Element.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Element.h (revision 4054) +++ /issm/trunk/src/c/objects/Elements/Element.h (revision 4055) @@ -37,11 +37,11 @@ virtual void GetBedList(double* bed_list)=0; virtual void Du(Vec du_g)=0; - virtual void Gradj(int control_type)=0; - virtual void GradjDrag(void)=0; - virtual void GradjB(void)=0; + virtual void Gradj(Vec gradient,int control_type)=0; + virtual void GradjDrag(Vec gradient)=0; + virtual void GradjB(Vec gradient)=0; virtual double Misfit(void)=0; virtual double CostFunction(void)=0; virtual double SurfaceArea(void)=0; - virtual void DepthAverageInputAtBase(int enum_type)=0; + virtual void DepthAverageInputAtBase(int enum_type,int average_enum_type)=0; virtual void ComputeBasalStress(Vec sigma_b)=0; virtual void ComputePressure(Vec p_g)=0; @@ -69,4 +69,6 @@ virtual void AXPYInput(int YEnum, double scalar, int XEnum)=0; virtual void ControlConstrainInput(int control_type,double cm_min, double cm_max)=0; + virtual void InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double criterionvalues,double* criterionvalues,int num_criterionenums)=0; + /*Implementation: */ Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 4054) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 4055) @@ -658,5 +658,5 @@ } - else if (analysis_type==SlopeComputeAnalysisEnum){ + else if (analysis_type==SlopeAnalysisEnum){ UpdateInputsFromSolutionSlopeCompute( solution); @@ -1270,5 +1270,5 @@ else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"); } - else if (analysis_type==SlopeComputeAnalysisEnum){ + else if (analysis_type==SlopeAnalysisEnum){ CreateKMatrixSlopeCompute( Kgg); } @@ -2357,5 +2357,5 @@ else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"); } - else if (analysis_type==SlopeComputeAnalysisEnum){ + else if (analysis_type==SlopeAnalysisEnum){ CreatePVectorSlopeCompute( pg); } @@ -3263,12 +3263,10 @@ } /*}}}*/ -/*FUNCTION Penta::FieldExtrude {{{1*/ -void Penta::FieldExtrude(Vec field,double* field_serial,char* field_name, int iscollapsed){ +/*FUNCTION Penta::VecExtrude {{{1*/ +void Penta::VecExtrude(Vec vector,double* vector_serial,int iscollapsed){ /* node data: */ - const int numgrids=6; - int numberofdofspernode; + int i; Node* node=NULL; - int i; int extrude=0; @@ -3296,116 +3294,35 @@ if(onbed==0)extrude=0; - /*Go on and extrude field: */ + /*Go on and extrude vector: */ if (extrude){ - if (strcmp(field_name,"velocity")==0){ - - /* node data: */ - const int numdof=2*numgrids; - int doflist[numdof]; - int nodedofs[2]; - double fieldel[2]; - - - GetDofList(&doflist[0],&numberofdofspernode); - - /*this penta is a collapsed macayeal. For each node on the base of this penta, - * we grab the field. Once we know the field, we follow the upper nodes, - * inserting the same field value into field, until we reach the surface: */ - for(i=0;i<3;i++){ - - node=nodes[i]; //base nodes - - /*get field for this base node: */ - fieldel[0]=field_serial[doflist[numberofdofspernode*i+0]]; - fieldel[1]=field_serial[doflist[numberofdofspernode*i+1]]; - - //go throfieldn all nodes which sit on top of this node, until we reach the surface, - //and plfield field in field - for(;;){ - - node->GetDofList(&nodedofs[0],&numberofdofspernode); - VecSetValues(field,1,&nodedofs[0],&fieldel[0],INSERT_VALUES); - VecSetValues(field,1,&nodedofs[1],&fieldel[1],INSERT_VALUES); - - if (node->IsOnSurface())break; - /*get next node: */ - node=node->GetUpperNode(); - } - } - } //if (strcmp(field_name,"velocity")==0) - else if (strcmp(field_name,"gradj")==0){ - - /* node data: */ - int dof1; - double fieldel; - - /*this penta is a collapsed macayeal. For each node on the base of this penta, - * we grab the field. Once we know the field, we follow the upper nodes, - * inserting the same field value into field, until we reach the surface: */ - for(i=0;i<3;i++){ - - node=nodes[i]; //base nodes + /* node data: */ + int dof1; + double vectorel; + + /*this penta is a collapsed macayeal. For each node on the base of this penta, + * we grab the vector. Once we know the vector, we follow the upper nodes, + * inserting the same vector value into vector, until we reach the surface: */ + for(i=0;i<3;i++){ + + node=nodes[i]; //base nodes + dof1=node->GetDofList1(); + + /*get vector for this base node: */ + vectorel=vector_serial[dof1]; + + //go throvectorn all nodes which sit on top of this node, until we reach the surface, + //and plvector vector in vector + for(;;){ + dof1=node->GetDofList1(); - - /*get field for this base node: */ - fieldel=field_serial[dof1]; - - //go throfieldn all nodes which sit on top of this node, until we reach the surface, - //and plfield field in field - for(;;){ - - dof1=node->GetDofList1(); - VecSetValues(field,1,&dof1,&fieldel,INSERT_VALUES); - - if (node->IsOnSurface())break; - /*get next node: */ - node=node->GetUpperNode(); - } + VecSetValues(vector,1,&dof1,&vectorel,INSERT_VALUES); + + if (node->IsOnSurface())break; + /*get next node: */ + node=node->GetUpperNode(); } } - else if ( - (strcmp(field_name,"thickness")==0) || - (strcmp(field_name,"surface")==0) || - (strcmp(field_name,"bed")==0) || - (strcmp(field_name,"slopex")==0) || - (strcmp(field_name,"slopey")==0) - ){ - - /* node data: */ - const int numdof=1*numgrids; - int doflist[numdof]; - int nodedofs; - double fieldel; - - GetDofList(&doflist[0],&numberofdofspernode); - - /*this penta is on the bed. For each node on the base of this penta, - * we grab the thickness. Once we know the thickness, we follow the upper nodes, - * inserting the same thickness value into tg, until we reach the surface: */ - for(i=0;i<3;i++){ - - node=nodes[i]; //base nodes - - /*get velocity for this base node: */ - fieldel=field_serial[doflist[numberofdofspernode*i+0]]; - - //go through all nodes which sit on top of this node, until we reach the surface, - //and pltg fieldel in field: - for(;;){ - - node->GetDofList(&nodedofs,&numberofdofspernode); - VecSetValues(field,1,&nodedofs,&fieldel,INSERT_VALUES); - - if (node->IsOnSurface())break; - /*get next node: */ - node=node->GetUpperNode(); - } - } - - } - else ISSMERROR("%s%s%s"," field ",field_name," not supported yet!"); - - } //if (extrude) + } } /*}}}*/ @@ -3451,5 +3368,5 @@ /*}}}*/ /*FUNCTION Penta::DepthAverageInputAtBase{{{1*/ -void Penta::DepthAverageInputAtBase(int enum_type){ +void Penta::DepthAverageInputAtBase(int enum_type,int average_enum_type){ ISSMERROR("Not implemented yet (see Penta::InputExtrude and Node::FieldDepthAverageAtBase)"); } @@ -4555,5 +4472,5 @@ /*}}}*/ /*FUNCTION Penta::Gradj {{{1*/ -void Penta::Gradj(int control_type){ +void Penta::Gradj(Vec gradient,int control_type){ /*inputs: */ @@ -4567,8 +4484,8 @@ if (control_type==DragCoefficientEnum){ - GradjDrag(); + GradjDrag(gradient); } else if (control_type=RheologyBEnum){ - GradjB(); + GradjB(gradient); } else ISSMERROR("%s%i","control type not supported yet: ",control_type); @@ -4576,10 +4493,10 @@ /*}}}*/ /*FUNCTION Penta::GradjDrag {{{1*/ -void Penta::GradjDrag(void){ +void Penta::GradjDrag(Vec gradient){ int i; Tria* tria=NULL; TriaVertexInput* triavertexinput=NULL; - double gradient[6]={0,0,0,0,0,0}; + double temp_gradient[6]={0,0,0,0,0,0}; /*inputs: */ @@ -4611,5 +4528,6 @@ /*MacAyeal or Pattyn*/ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->GradjDrag(); + tria->GradjDrag(gradient); + delete tria; } else if (sub_analysis_type==StokesAnalysisEnum){ @@ -4617,24 +4535,18 @@ /*Stokes*/ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->GradjDragStokes(); + tria->GradjDragStokes(gradient); + delete tria; } else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"); - /*Now, the tria has a GradientEnum input. Take it, and make a PentaVertexInput out of it, then delete the tria: */ - triavertexinput=(TriaVertexInput*)tria->inputs->GetInput(GradientEnum); - for(i=0;i<3;i++)gradient[i]=triavertexinput->values[i]; - this->inputs->AddInput(new PentaVertexInput(GradientEnum,&gradient[0])); - - delete tria; } /*}}}*/ /*FUNCTION Penta::GradjB {{{1*/ -void Penta::GradjB(void){ +void Penta::GradjB(Vec gradient){ int i; Tria* tria=NULL; TriaVertexInput* triavertexinput=NULL; - double gradient[6]={0,0,0,0,0,0}; /*inputs: */ @@ -4658,18 +4570,14 @@ * and compute gardj*/ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - tria->GradjB(); + tria->GradjB(gradient); + delete tria; } else{ /*B is a 2d field, use MacAyeal(2d) gradient even if it is Stokes or Pattyn*/ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - tria->GradjB(); + tria->GradjB(gradient); + delete tria; } - /*Now, the tria has a GradientEnum input. Take it, and make a PentaVertexInput out of it, then delete the tria: */ - triavertexinput=(TriaVertexInput*)tria->inputs->GetInput(GradientEnum); - for(i=0;i<3;i++)gradient[i]=triavertexinput->values[i]; - this->inputs->AddInput(new PentaVertexInput(GradientEnum,&gradient[0])); - - delete tria; } @@ -5407,2 +5315,30 @@ } /*}}}*/ +/*FUNCTION Penta::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/ +void Penta::InputConvergence(int* pconverged,double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){a + + Input** new_inputs=NULL; + Input** old_inputs=NULL; + int converged=1; + + new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs + old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs + + for(i=0;iinputs->GetInput(new_enums[2*i+0]); + old_inputs[i]=(Input*)this->inputs->GetInput(old_enums[2*i+1]); + if(!new_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(new_enums[2*i+0])); + if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(old_enums[2*i+0])); + } + + /*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/ + for(i=0;icriterionvalues[i]) converged=0; + } + + /*Assign output pointers:*/ + *pconverged=converged; + +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 4054) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 4055) @@ -86,7 +86,7 @@ bool GetOnBed(); void Du(Vec du_gg); - void Gradj(int control_type); - void GradjDrag(void); - void GradjB(void); + void Gradj(Vec gradient,int control_type); + void GradjDrag(Vec gradient); + void GradjB(Vec gradient); double Misfit(void); double SurfaceArea(void); @@ -110,7 +110,7 @@ void GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_coord); void GetNodalFunctions(double* l1l6, double* gauss_coord); - void FieldExtrude(Vec field,double* field_serial,char* field_name, int iscollapsed); + void VecExtrude(Vec vector,double* vector_serial,int iscollapsed); void InputExtrude(int enum_type); - void DepthAverageInputAtBase(int enum_type); + void DepthAverageInputAtBase(int enum_type,int average_enum_type); void ComputeBasalStress(Vec sigma_bg); void ComputePressure(Vec p_gg); @@ -165,4 +165,5 @@ void AXPYInput(int YEnum, double scalar, int XEnum); void ControlConstrainInput(int control_type,double cm_min, double cm_max); + void InputConvergence(int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums); void GetVectorFromInputs(Vec vector,int NameEnum); Index: /issm/trunk/src/c/objects/Elements/Sing.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Sing.cpp (revision 4054) +++ /issm/trunk/src/c/objects/Elements/Sing.cpp (revision 4055) @@ -383,15 +383,15 @@ /*}}}*/ /*FUNCTION Sing::Gradj {{{1*/ -void Sing::Gradj(int control_type){ +void Sing::Gradj(Vec gradient,int control_type){ ISSMERROR(" not supported yet!"); } /*}}}*/ /*FUNCTION Sing::GradB {{{1*/ -void Sing::GradjB(void){ +void Sing::GradjB(Vec gradient){ ISSMERROR(" not supported yet!"); } /*}}}*/ /*FUNCTION Sing::GradjDrag {{{1*/ -void Sing::GradjDrag(void){ +void Sing::GradjDrag(Vec gradient){ ISSMERROR(" not supported yet!"); } @@ -761,2 +761,30 @@ } /*}}}*/ +/*FUNCTION Sing::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/ +void Sing::InputConvergence(int* pconverged,double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){a + + Input** new_inputs=NULL; + Input** old_inputs=NULL; + int converged=1; + + new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs + old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs + + for(i=0;iinputs->GetInput(new_enums[2*i+0]); + old_inputs[i]=(Input*)this->inputs->GetInput(old_enums[2*i+1]); + if(!new_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(new_enums[2*i+0])); + if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(old_enums[2*i+0])); + } + + /*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/ + for(i=0;icriterionvalues[i]) converged=0; + } + + /*Assign output pointers:*/ + *pconverged=converged; + +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Sing.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Sing.h (revision 4054) +++ /issm/trunk/src/c/objects/Elements/Sing.h (revision 4055) @@ -62,5 +62,5 @@ void UpdateInputsFromConstant(bool constant, int name){ISSMERROR("Not implemented yet!");} void UpdateInputsFromSolution(double* solutiong); - void DepthAverageInputAtBase(int enum_type){ISSMERROR("not implemented yet");}; + void DepthAverageInputAtBase(int enum_type,int average_enum_type){ISSMERROR("not implemented yet");}; void InputToResult(int enum_type,int step,double time); void ProcessResultsUnits(void); @@ -96,4 +96,5 @@ void ScaleInput(int enum_type,double scale_factor); void AXPYInput(int YEnum, double scalar, int XEnum); + void InputConvergence(int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums); void ControlConstrainInput(int control_type,double cm_min, double cm_max); @@ -106,7 +107,7 @@ void GetThicknessList(double* thickness_list); void Du(Vec); - void Gradj(int control_type); - void GradjDrag(void); - void GradjB(void); + void Gradj(Vec gradient,int control_type); + void GradjDrag(Vec gradient); + void GradjB(Vec gradient); double Misfit(void); double SurfaceArea(void); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 4054) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 4055) @@ -569,5 +569,5 @@ else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"); } - else if (analysis_type==SlopeComputeAnalysisEnum){ + else if (analysis_type==SlopeAnalysisEnum){ UpdateInputsFromSolutionSlopeCompute( solution); } @@ -979,5 +979,5 @@ else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"); } - else if (analysis_type==SlopeComputeAnalysisEnum){ + else if (analysis_type==SlopeAnalysisEnum){ CreateKMatrixSlopeCompute( Kgg); } @@ -2394,5 +2394,5 @@ else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"); } - else if (analysis_type==SlopeComputeAnalysisEnum){ + else if (analysis_type==SlopeAnalysisEnum){ CreatePVectorSlopeCompute( pg); } @@ -4125,5 +4125,5 @@ /*}}}*/ /*FUNCTION Tria::Gradj {{{1*/ -void Tria::Gradj(int control_type){ +void Tria::Gradj(Vec gradient,int control_type){ /*inputs: */ @@ -4137,8 +4137,8 @@ if (control_type==DragCoefficientEnum){ - GradjDrag(); + GradjDrag(gradient); } else if (control_type==RheologyBEnum){ - GradjB(); + GradjB(gradient); } else ISSMERROR("%s%i","control type not supported yet: ",control_type); @@ -4146,5 +4146,5 @@ /*}}}*/ /*FUNCTION Tria::GradjDrag {{{1*/ -void Tria::GradjDrag(void){ +void Tria::GradjDrag(Vec gradient){ @@ -4295,6 +4295,6 @@ - /*Add grade_g to the inputs of this element: */ - this->inputs->AddInput(new TriaVertexInput(GradientEnum,&grade_g[0])); + /*Add grade_g to global vector gradient: */ + VecSetValues(gradient,numgrids,doflist1,(const double*)grade_g,ADD_VALUES); cleanup_and_return: @@ -4308,5 +4308,5 @@ /*}}}*/ /*FUNCTION Tria::GradjDragStokes {{{1*/ -void Tria::GradjDragStokes(void){ +void Tria::GradjDragStokes(Vec gradient){ int i; @@ -4465,4 +4465,7 @@ } + /*Add grade_g to global vector gradient: */ + VecSetValues(gradient,numgrids,doflist1,(const double*)grade_g,ADD_VALUES); + /*Add grade_g to the inputs of this element: */ this->inputs->AddInput(new TriaVertexInput(GradientEnum,&grade_g[0])); @@ -4478,5 +4481,5 @@ /*}}}*/ /*FUNCTION Tria::GradjB{{{1*/ -void Tria::GradjB(void){ +void Tria::GradjB(Vec gradient){ int i; @@ -4613,6 +4616,6 @@ } - /*Add grade_g to the inputs of this element: */ - this->inputs->AddInput(new TriaVertexInput(GradientEnum,&grade_g[0])); + /*Add grade_g to global vector gradient: */ + VecSetValues(gradient,numgrids,doflist1,(const double*)grade_g,ADD_VALUES); cleanup_and_return: @@ -5014,5 +5017,5 @@ /*}}}*/ /*FUNCTION Tria::DepthAverageInputAtBase {{{1*/ -void Tria::DepthAverageInputAtBase(int enum_type){ +void Tria::DepthAverageInputAtBase(int enum_type,int average_enum_type){ /*New input*/ @@ -5026,17 +5029,5 @@ /*Assign new name (average)*/ - switch(enum_type){ - - case VxEnum: - newinput->ChangeEnum(VxAverageEnum); - break; - - case VyEnum: - newinput->ChangeEnum(VyAverageEnum); - break; - - default: - ISSMERROR("enum_type: %i (%s) not suppoerted yet",enum_type,EnumAsString(enum_type)); - } + newinput->ChangeEnum(average_enum_type); /*Add new input to current element*/ @@ -5503,2 +5494,30 @@ } /*}}}*/ +/*FUNCTION Tria::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/ +void Tria::InputConvergence(int* pconverged,double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){a + + Input** new_inputs=NULL; + Input** old_inputs=NULL; + int converged=1; + + new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs + old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs + + for(i=0;iinputs->GetInput(new_enums[2*i+0]); + old_inputs[i]=(Input*)this->inputs->GetInput(old_enums[2*i+1]); + if(!new_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(new_enums[2*i+0])); + if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(old_enums[2*i+0])); + } + + /*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/ + for(i=0;icriterionvalues[i]) converged=0; + } + + /*Assign output pointers:*/ + *pconverged=converged; + +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 4054) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 4055) @@ -55,5 +55,5 @@ int MyRank(void); void SetClone(int* minranks); - void DepthAverageInputAtBase(int enum_type); + void DepthAverageInputAtBase(int enum_type,int average_enum_type); void* SpawnSing(int g0); void* SpawnBeam(int g0, int g1); @@ -88,8 +88,8 @@ void GetJacobian(double* J, double* xyz_list,double* gauss_l1l2l3); void Du(Vec du_g); - void Gradj(int control_type); - void GradjDrag(void); - void GradjDragStokes(void); - void GradjB(void); + void Gradj(Vec gradient,int control_type); + void GradjDrag(Vec gradient); + void GradjDragStokes(Vec gradient); + void GradjB(Vec gradient); void SurfaceNormal(double* surface_normal, double xyz_list[3][3]); double Misfit(void); @@ -142,4 +142,5 @@ void ScaleInput(int enum_type,double scale_factor); void AXPYInput(int YEnum, double scalar, int XEnum); + void InputConvergence(int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums); void ControlConstrainInput(int control_type,double cm_min, double cm_max); void GetVectorFromInputs(Vec vector,int NameEnum); Index: /issm/trunk/src/c/objects/FemModel.cpp =================================================================== --- /issm/trunk/src/c/objects/FemModel.cpp (revision 4054) +++ /issm/trunk/src/c/objects/FemModel.cpp (revision 4055) @@ -35,17 +35,17 @@ /*Dynamically allocate whatever is a list of length nummodels: */ analysis_type_list=(int*)xmalloc(nummodels*sizeof(int)); - Rmg=(Mat*)xmalloc(nummodels*sizeof(Mat)); - Gmn=(Mat*)xmalloc(nummodels*sizeof(Mat)); - nodesets=(NodeSets**)xmalloc(nummodels*sizeof(NodeSets*)); - yg=(Vec*)xmalloc(nummodels*sizeof(Vec)); - ys=(Vec*)xmalloc(nummodels*sizeof(Vec)); + m_Rmg=(Mat*)xmalloc(nummodels*sizeof(Mat)); + m_Gmn=(Mat*)xmalloc(nummodels*sizeof(Mat)); + m_nodesets=(NodeSets**)xmalloc(nummodels*sizeof(NodeSets*)); + m_yg=(Vec*)xmalloc(nummodels*sizeof(Vec)); + m_ys=(Vec*)xmalloc(nummodels*sizeof(Vec)); /*Initialize: */ for(i=0;ivalues; + *pnum_values=2; + +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Inputs/BeamVertexInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/BeamVertexInput.h (revision 4054) +++ /issm/trunk/src/c/objects/Inputs/BeamVertexInput.h (revision 4055) @@ -78,4 +78,5 @@ void Constrain(double cm_min, double cm_max); void GetVectorFromInputs(Vec vector,int* doflist); + void GetValuesPtr(double* pvalues,int* pnum_values); /*}}}*/ Index: /issm/trunk/src/c/objects/Inputs/BoolInput.cpp =================================================================== --- /issm/trunk/src/c/objects/Inputs/BoolInput.cpp (revision 4054) +++ /issm/trunk/src/c/objects/Inputs/BoolInput.cpp (revision 4055) @@ -256,2 +256,9 @@ } /*}}}*/ +/*FUNCTION BoolInput::GetValuesPtr(double* pvalues,int* pnum_values);{{{1*/ +void BoolInput::GetValuesPtr(double* pvalues,int* pnum_values){ + + ISSMERROR(" not supported yet!"); + +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Inputs/BoolInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/BoolInput.h (revision 4054) +++ /issm/trunk/src/c/objects/Inputs/BoolInput.h (revision 4055) @@ -78,4 +78,5 @@ void Constrain(double cm_min, double cm_max); void GetVectorFromInputs(Vec vector,int* doflist); + void GetValuesPtr(double* pvalues,int* pnum_values); /*}}}*/ Index: /issm/trunk/src/c/objects/Inputs/DoubleInput.cpp =================================================================== --- /issm/trunk/src/c/objects/Inputs/DoubleInput.cpp (revision 4054) +++ /issm/trunk/src/c/objects/Inputs/DoubleInput.cpp (revision 4055) @@ -267,2 +267,9 @@ } /*}}}*/ +/*FUNCTION DoubleInput::GetValuesPtr(double* pvalues,int* pnum_values);{{{1*/ +void DoubleInput::GetValuesPtr(double* pvalues,int* pnum_values){ + + ISSMERROR(" not supported yet!"); + +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Inputs/DoubleInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/DoubleInput.h (revision 4054) +++ /issm/trunk/src/c/objects/Inputs/DoubleInput.h (revision 4055) @@ -79,4 +79,5 @@ void Constrain(double cm_min, double cm_max); void GetVectorFromInputs(Vec vector,int* doflist); + void GetValuesPtr(double* pvalues,int* pnum_values); /*}}}*/ Index: /issm/trunk/src/c/objects/Inputs/Input.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/Input.h (revision 4054) +++ /issm/trunk/src/c/objects/Inputs/Input.h (revision 4055) @@ -55,4 +55,5 @@ virtual void Constrain(double cm_min, double cm_max)=0; virtual void GetVectorFromInputs(Vec vector,int* doflist)=0; + virtual void GetValuesPtr(double* pvalues,int* pnum_values)=0; /*}}}*/ Index: /issm/trunk/src/c/objects/Inputs/IntInput.cpp =================================================================== --- /issm/trunk/src/c/objects/Inputs/IntInput.cpp (revision 4054) +++ /issm/trunk/src/c/objects/Inputs/IntInput.cpp (revision 4055) @@ -257,2 +257,9 @@ } /*}}}*/ +/*FUNCTION IntInput::GetValuesPtr(double* pvalues,int* pnum_values);{{{1*/ +void IntInput::GetValuesPtr(double* pvalues,int* pnum_values){ + + ISSMERROR(" not supported yet!"); + +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Inputs/IntInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/IntInput.h (revision 4054) +++ /issm/trunk/src/c/objects/Inputs/IntInput.h (revision 4055) @@ -78,4 +78,5 @@ void Constrain(double cm_min, double cm_max); void GetVectorFromInputs(Vec vector,int* doflist); + void GetValuesPtr(double* pvalues,int* pnum_values); /*}}}*/ Index: /issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp =================================================================== --- /issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp (revision 4054) +++ /issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp (revision 4055) @@ -943,2 +943,10 @@ } /*}}}*/ +/*FUNCTION PentaVertexInput::GetValuesPtr(double* pvalues,int* pnum_values);{{{1*/ +void PentaVertexInput::GetValuesPtr(double* pvalues,int* pnum_values){ + + *pvalues=this->values; + *pnum_values=6; + +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Inputs/PentaVertexInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/PentaVertexInput.h (revision 4054) +++ /issm/trunk/src/c/objects/Inputs/PentaVertexInput.h (revision 4055) @@ -87,4 +87,5 @@ void Constrain(double cm_min, double cm_max); void GetVectorFromInputs(Vec vector,int* doflist); + void GetValuesPtr(double* pvalues,int* pnum_values); /*}}}*/ Index: /issm/trunk/src/c/objects/Inputs/SingVertexInput.cpp =================================================================== --- /issm/trunk/src/c/objects/Inputs/SingVertexInput.cpp (revision 4054) +++ /issm/trunk/src/c/objects/Inputs/SingVertexInput.cpp (revision 4055) @@ -258,3 +258,11 @@ } -/*}}}*/ +/2*}}}*/ +/*FUNCTION SingVertexInput::GetValuesPtr(double* pvalues,int* pnum_values);{{{1*/ +void SingVertexInput::GetValuesPtr(double* pvalues,int* pnum_values){ + + *pvalues=&this->value; + *pnum_values=1; + +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Inputs/SingVertexInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/SingVertexInput.h (revision 4054) +++ /issm/trunk/src/c/objects/Inputs/SingVertexInput.h (revision 4055) @@ -77,4 +77,5 @@ void Constrain(double cm_min, double cm_max); void GetVectorFromInputs(Vec vector,int* doflist); + void GetValuesPtr(double* pvalues,int* pnum_values); /*}}}*/ Index: /issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp =================================================================== --- /issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp (revision 4054) +++ /issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp (revision 4055) @@ -517,2 +517,10 @@ } /*}}}*/ +/*FUNCTION TriaVertexInput::GetValuesPtr(double* pvalues,int* pnum_values);{{{1*/ +void TriaVertexInput::GetValuesPtr(double* pvalues,int* pnum_values){ + + *pvalues=this->values; + *pnum_values=3; + +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Inputs/TriaVertexInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/TriaVertexInput.h (revision 4054) +++ /issm/trunk/src/c/objects/Inputs/TriaVertexInput.h (revision 4055) @@ -84,4 +84,5 @@ void Constrain(double cm_min, double cm_max); void GetVectorFromInputs(Vec vector,int* doflist); + void GetValuesPtr(double* pvalues,int* pnum_values); /*}}}*/ Index: /issm/trunk/src/c/objects/Node.cpp =================================================================== --- /issm/trunk/src/c/objects/Node.cpp (revision 4054) +++ /issm/trunk/src/c/objects/Node.cpp (revision 4055) @@ -144,5 +144,5 @@ analysis_type==PrognosticAnalysisEnum || analysis_type==MeltingAnalysisEnum || - analysis_type==SlopeComputeAnalysisEnum || + analysis_type==SlopeAnalysisEnum || analysis_type==BalancedvelocitiesAnalysisEnum || analysis_type==BalancedthicknessAnalysisEnum @@ -668,6 +668,6 @@ } /*}}}*/ -/*FUNCTION Node::FieldExtrude {{{2*/ -void Node::FieldExtrude(Vec field,double* field_serial,char* field_name){ +/*FUNCTION Node::VecExtrude {{{2*/ +void Node::VecExtrude(Vec vector,double* vector_serial){ /* node data: */ @@ -684,103 +684,27 @@ if (onbed){ - if (strcmp(field_name,"velocity")==0){ - - /* node data: */ - int vertexdof; - int doflist[2]; - double fieldnode[2]; - - vertexdof=GetVertexDof(); - - /*get dofs for this base node velocity: we know there are two dofs in field_serial */ - doflist[0]=2*vertexdof; - doflist[1]=2*vertexdof+1; - - //initilaize node - node=this; - - /*get field for this base node: */ - fieldnode[0]=field_serial[doflist[0]]; - fieldnode[1]=field_serial[doflist[1]]; - - //go through all nodes which sit on top of this node, until we reach the surface, - //and plfield field in field - for(;;){ - - vertexdof=node->GetVertexDof(); - doflist[0]=2*vertexdof; - doflist[1]=2*vertexdof+1; - VecSetValues(field,1,&doflist[0],&fieldnode[0],INSERT_VALUES); - VecSetValues(field,1,&doflist[1],&fieldnode[1],INSERT_VALUES); - - if (node->IsOnSurface())break; - /*get next node: */ - node=node->GetUpperNode(); - } - } //if (strcmp(field_name,"velocity")==0) - else if (strcmp(field_name,"gradj")==0){ - - /* node data: */ - int dof1; - double fieldel; - - //initilaize node and get dof1 - node=this; + /* node data: */ + int dof1; + double vectorel; + + //initilaize node and get dof1 + node=this; + dof1=node->GetVertexDof(); + + /*get vector for this base node: */ + vectorel=vector_serial[dof1]; + + //go throvectorn all nodes which sit on top of this node, until we reach the surface, + //and plvector vector in vector + for(;;){ + dof1=node->GetVertexDof(); - - /*get field for this base node: */ - fieldel=field_serial[dof1]; - - //go throfieldn all nodes which sit on top of this node, until we reach the surface, - //and plfield field in field - for(;;){ - - dof1=node->GetVertexDof(); - VecSetValues(field,1,&dof1,&fieldel,INSERT_VALUES); - - if (node->IsOnSurface())break; - /*get next node: */ - node=node->GetUpperNode(); - } - } - else if ( - (strcmp(field_name,"vx")==0) || - (strcmp(field_name,"vy")==0) || - (strcmp(field_name,"thickness")==0) || - (strcmp(field_name,"temperature")==0) || - (strcmp(field_name,"surface")==0) || - (strcmp(field_name,"bed")==0) || - (strcmp(field_name,"slopex")==0) || - (strcmp(field_name,"slopey")==0) - ){ - - /* node data: */ - int doflist; - int nodedofs; - double fieldnode; - - GetDofList(&doflist,&numberofdofspernode); - - //initialize node - node=this; - - /*get field for this bed node: */ - fieldnode=field_serial[doflist]; - - //go through all nodes which sit on top of this node, until we reach the surface, - //and plug fieldnode in field: - for(;;){ - - node->GetDofList(&nodedofs,&numberofdofspernode); - VecSetValues(field,1,&nodedofs,&fieldnode,INSERT_VALUES); - - if (node->IsOnSurface())break; - /*get next node: */ - node=node->GetUpperNode(); - } - } - else ISSMERROR("%s%s%s"," field ",field_name," not supported yet!"); - - } //if (extrude) + VecSetValues(vector,1,&dof1,&vectorel,INSERT_VALUES); + + if (node->IsOnSurface())break; + /*get next node: */ + node=node->GetUpperNode(); + } + } } /*}}}*/ Index: /issm/trunk/src/c/objects/Node.h =================================================================== --- /issm/trunk/src/c/objects/Node.h (revision 4054) +++ /issm/trunk/src/c/objects/Node.h (revision 4055) @@ -89,5 +89,5 @@ int IsOnShelf(); int IsOnSheet(); - void FieldExtrude(Vec field,double* field_serial,char* field_name); + void VecExtrude(Vec vector,double* vector_serial); /*}}}*/ /*FUNCTION DofObject routines {{{1*/ Index: /issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp =================================================================== --- /issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp (revision 4054) +++ /issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp (revision 4055) @@ -25,5 +25,5 @@ numdofs=2; } - else if (analysis_type==SlopeComputeAnalysisEnum){ + else if (analysis_type==SlopeAnalysisEnum){ numdofs=1; } Index: /issm/trunk/src/c/shared/Numerics/IsInputConverged.cpp =================================================================== --- /issm/trunk/src/c/shared/Numerics/IsInputConverged.cpp (revision 4055) +++ /issm/trunk/src/c/shared/Numerics/IsInputConverged.cpp (revision 4055) @@ -0,0 +1,54 @@ +/*!\file: IsInputConverged.cpp + * \brief: figure out convergence at the input level. + * We are given a list of inputs, new and old, and a criterion, and using the inputs, we + * return the value of the criterion test, which will be used at the solution level + * to determine convergence. + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../../objects/objects.h" + +void IsInputConverged(double* peps, Input** new_inputs,Input** old_inputs,int num_inputs,int criterion_enum){ + + int i; + + /*output: */ + double eps; + + /*intermediary: */ + double *newvalues = NULL; + double *oldvalues = NULL; + int num_values; + double ndu = 0; + double nu = 0; + + if(criterion_enum==RelativeEnum){ + + /*conpute ndu/du (where u could be velocity, pressure, temperature, etc ...): */ + for(i=0;iGetValuesPtr(&newvalues,&num_values); + old_inputs[i]->GetValuesPtr(&oldvalues,&num_values); + for(j=0;jSetAnalysis(DiagnosticAnalysisEnum,DiagnosticStokesAnalysisEnum); - else femmodel->SetAnalysis(DiagnosticAnalysisEnum,DiagnosticHorizAnalysisEnum); + if(isstokes)femmodel->SetCurrentAnalysis(DiagnosticStokesAnalysisEnum); + else femmodel->SetCurrentAnalysis(DiagnosticHorizAnalysisEnum); _printf_("%s\n"," recover solution for this stiffness and right hand side:"); - solver_diagnostic_nonlinear(NULL,&K_ff0,&K_fs0,NULL, femmodel,DiagnosticAnalysisEnum,sub_analysis_type); + solver_diagnostic_nonlinear(NULL,&K_ff0,&K_fs0, femmodel,conserve_loads); _printf_("%s\n"," buid Du, difference between observed velocity and model velocity:"); - Dux(du_g, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + Dux(&du_g, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); _printf_("%s\n"," reduce adjoint load from g-set to f-set:"); @@ -66,14 +67,14 @@ /*Update inputs using adjoint solution, and same type of setup as diagnostic solution: */ - if(isstokes)femmodel->SetAnalysisAlias(DiagnosticAnalysisEnum,DiagnosticStokesAnalysisEnum,AdjointAnalysisEnum,NoneAnalysisEnum); - else femmodel->SetAnalysisAlias(DiagnosticAnalysisEnum,DiagnosticHorizAnalysisEnum,AdjointAnalysisEnum,NoneAnalysisEnum); + if(isstokes)femmodel->SetCurrentAnalysisAlias(DiagnosticStokesAnalysisEnum,AdjointAnalysisEnum); + else femmodel->SetCurrentAnalysisAlias(DiagnosticHorizAnalysisEnum,AdjointAnalysisEnum); - UpdateInputsFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,adjoint_g,AdjointAnalysisEnum, NoneAnalysisEnum); + UpdateInputsFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,adjoint_g); /*Free ressources:*/ - char* solverstring=NULL; - VecFree(&ug); + xfree((void**)&solverstring); + VecFree(&u_g); MatFree(&K_ff0); - MatFree(&Mat K_fs0); + MatFree(&K_fs0); VecFree(&du_g); VecFree(&du_f); Index: /issm/trunk/src/c/solutions/balancedthickness.cpp =================================================================== --- /issm/trunk/src/c/solutions/balancedthickness.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/balancedthickness.cpp (revision 4055) @@ -22,17 +22,10 @@ char* inputfilename=NULL; char* outputfilename=NULL; + bool qmu_analysis=false; char* lockname=NULL; - int numberofnodes; bool waitonlock=false; - Model* model=NULL; - - double dt; - double yts; - bool qmu_analysis; - - /*Results: */ - Results* results=NULL; - Param* param=NULL; + /*FemModel: */ + FemModel* femmodel=NULL; /*time*/ @@ -40,4 +33,7 @@ double start_core, finish_core; double start_init, finish_init; + + int analyses[1]={BalancedthicknessAnalysisEnum}; + int solution_type=BalancedthicknessAnalysisEnum; MODULEBOOT(); @@ -66,13 +62,14 @@ fid=pfopen(inputfilename,"rb"); - /*Initialize model structure: */ - model=new Model(); + _printf_("create finite element model:\n"); + femmodel=new FemModel(fid,solution_type,analyses,1); - _printf_("read and create finite element model:\n"); - model->AddFormulation(fid,BalancedthicknessAnalysisEnum); - - /*recover parameters: */ - model->FindParam(&waitonlock,WaitOnLockEnum); - model->FindParam(&qmu_analysis,QmuAnalysisEnum); + /*add outputfilename in parameters: */ + femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename); + + /*get parameters: */ + femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum); + femmodel->parameters->FindParam(&control_analysis,ControlAnalysisEnum); + femmodel->parameters->FindParam(&waitonlock,WaitOnLockEnum); MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime(); @@ -84,6 +81,10 @@ _printf_("call computational core:\n"); MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); - results=balancedthickness_core(model); + balancedthickness_core(model); MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + + _printf_("write results to disk:\n"); + OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters, outputfilename); + } } @@ -95,5 +96,5 @@ #ifdef _HAVE_DAKOTA_ MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); - Qmux(model,BalancedthicknessAnalysisEnum,NoneAnalysisEnum); + Qmux(femmodel,BalancedthicknessAnalysisEnum,NoneAnalysisEnum); MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); #else @@ -101,7 +102,4 @@ #endif } - - _printf_("write results to disk:\n"); - OutputResults(results,outputfilename); if (waitonlock>0){ @@ -111,6 +109,5 @@ /*Free ressources:*/ - delete results; - delete model; + delete femmodel; /*Get finish time and close*/ Index: /issm/trunk/src/c/solutions/balancedthickness2.cpp =================================================================== --- /issm/trunk/src/c/solutions/balancedthickness2.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/balancedthickness2.cpp (revision 4055) @@ -23,18 +23,10 @@ char* inputfilename=NULL; char* outputfilename=NULL; + bool qmu_analysis=false; char* lockname=NULL; bool waitonlock=false; - Model* model=NULL; - - double dt; - double yts; - bool qmu_analysis; - - /*Results: */ - Results* results=NULL; - Result* result=NULL; - - Param* param=NULL; + /*FemModel: */ + FemModel* femmodel=NULL; /*time*/ @@ -42,4 +34,8 @@ double start_core, finish_core; double start_init, finish_init; + + int analyses[1]={Balancedthickness2AnalysisEnum}; + int solution_type=Balancedthickness2AnalysisEnum; + MODULEBOOT(); @@ -68,13 +64,9 @@ fid=pfopen(inputfilename,"rb"); - /*Initialize model structure: */ - model=new Model(); + _printf_("create finite element model:\n"); + femmodel=new FemModel(fid,solution_type,analyses,1); - _printf_("read and create finite element model:\n"); - model->AddFormulation(fid,Balancedthickness2AnalysisEnum); - - /*recover parameters: */ - model->FindParam(&waitonlock,WaitOnLockEnum); - model->FindParam(&qmu_analysis,QmuAnalysisEnum); + /*add outputfilename in parameters: */ + femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename); MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime(); @@ -86,6 +78,9 @@ _printf_("call computational core:\n"); MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); - results=balancedthickness2_core(model); + balancedthickness2_core(femmodel); MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + + _printf_("write results to disk:\n"); + OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters, outputfilename); } @@ -97,5 +92,5 @@ #ifdef _HAVE_DAKOTA_ MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); - Qmux(model,Balancedthickness2AnalysisEnum,NoneAnalysisEnum); + Qmux(femmodel,Balancedthickness2AnalysisEnum,NoneAnalysisEnum); MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); #else @@ -103,7 +98,4 @@ #endif } - - _printf_("write results to disk:\n"); - OutputResults(results,outputfilename); if (waitonlock>0){ @@ -113,6 +105,5 @@ /*Free ressources:*/ - delete results; - delete model; + delete femmodel; /*Get finish time and close*/ Index: /issm/trunk/src/c/solutions/balancedthickness2_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/balancedthickness2_core.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/balancedthickness2_core.cpp (revision 4055) @@ -12,71 +12,35 @@ #include "../solvers/solvers.h" -Results* balancedthickness2_core(Model* model){ - - extern int my_rank; - - /*output: */ - Results* results=NULL; - Result* result=NULL; - - /*intermediary: */ - Vec vx_g=NULL; - Vec vy_g=NULL; - - /*solutions: */ - Vec h_g=NULL; +void balancedthickness2_core(FemModel* femmodel){ /*flags: */ int verbose=0; - int numberofdofspernode; - int numberofnodes; - int numberofvertices; - int dofs[1]={1}; + int dim; - /*fem balancedthickness2 model: */ - FemModel* fem_p=NULL; - - //initialize results: - results=new Results(); - - fem_p=model->GetFormulation(Balancedthickness2AnalysisEnum); - - //first recover parameters common to all solutions - model->FindParam(&verbose,VerboseEnum); - model->FindParam(&numberofnodes,NumberOfNodesEnum); - model->FindParam(&numberofvertices,NumberOfVerticesEnum); - model->FindParam(&numberofdofspernode,NumberOfDofsPerNodeEnum); + /*activate formulation: */ + femmodel->SetCurrentAnalysis(BalancedthicknessAnalysisEnum); + + /*recover parameters: */ + femmodel->parameters->FindParam(&verbose,VerboseEnum); + femmodel->parameters->FindParam(&dim,DimEnum); _printf_("depth averaging velocity...\n"); - //vx_g=inputs->Get("vx",&dofs[0],1); - //vy_g=inputs->Get("vy",&dofs[0],1); - /* NOT WORKING YET.... - FieldDepthAveragex(vx_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"vx"); - FieldDepthAveragex(vy_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"vy"); - */ - //inputs->Add("vx_average",vx_g,1,numberofvertices); - //inputs->Add("vy_average",vy_g,1,numberofvertices); - ISSMERROR("not supported yet!"); - + DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VxEnum); + DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VyEnum); + if(dim==3) DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VzEnum); + + _printf_("call computational core:\n"); - diagnostic_core_linear(&h_g,fem_p,Balancedthickness2AnalysisEnum,NoneAnalysisEnum); + solver_linear(NULL,femmodel); _printf_("Averaging over vertices:\n"); - FieldAverageOntoVerticesx(&h_g,fem_p->elements,fem_p->nodes,fem_p->vertices,fem_p->loads,fem_p->materials,fem_p->parameters); + ISSMERROR(" not supported yet!"); +//FieldAverageOntoVerticesx(&h_g,fem_p->elements,fem_p->nodes,fem_p->vertices,fem_p->loads,fem_p->materials,fem_p->parameters); // _printf_("extrude computed thickness on all layers:\n"); // FieldExtrudex( h_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0); - /*Plug results into output dataset: */ - InputToResultx(&result,fem_p->elements,fem_p->nodes,fem_p->vertices, fem_p->loads, fem_p->materials,fem_p->parameters,ThicknessEnum,results->Size()+1,0,1); results->AddObject(result); - results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(Balancedthickness2AnalysisEnum))); - - /*Free ressources:*/ - VecFree(&vx_g); - VecFree(&vy_g); - VecFree(&h_g); - - /*return: */ - return results; + if(verbose)_printf_("saving results:\n"); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); } Index: /issm/trunk/src/c/solutions/balancedthickness_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/balancedthickness_core.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/balancedthickness_core.cpp (revision 4055) @@ -12,60 +12,30 @@ #include "../solvers/solvers.h" -Results* balancedthickness_core(Model* model){ +void balancedthickness_core(FemModel* femmodel){ - extern int my_rank; + /*parameters: */ + int verbose=0; + int dim; - /*output: */ - Results* results=NULL; - Result* result=NULL; - - /*intermediary: */ - Vec u_g=NULL; - - /*solutions: */ - Vec h_g=NULL; - - /*flags: */ - int verbose=0; - int numberofdofspernode; - int numberofnodes; - int dofs[2]={1,1}; - - /*fem balancedthickness model: */ - FemModel* fem_p=NULL; - - //initialize results: - results=new Results(); - - /*recover fem model: */ - fem_p=model->GetFormulation(BalancedthicknessAnalysisEnum); - - //first recover parameters common to all solutions - model->FindParam(&verbose,VerboseEnum); - model->FindParam(&numberofnodes,NumberOfNodesEnum); - model->FindParam(&numberofdofspernode,NumberOfDofsPerNodeEnum); + /*activate formulation: */ + femmodel->SetCurrentAnalysis(BalancedthicknessAnalysisEnum); + + /*recover parameters: */ + femmodel->parameters->FindParam(&verbose,VerboseEnum); + femmodel->parameters->FindParam(&dim,DimEnum); _printf_("depth averaging velocity...\n"); - //u_g=inputs->Get("velocity",&dofs[0],2); //take (vx,vy) from inputs velocity - //FieldDepthAveragex( u_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"velocity"); - //inputs->Add("velocity_average",u_g,2,numberofnodes); - ISSMERROR("not supported yet!"); - + DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VxEnum); + DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VyEnum); + if(dim==3) DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VzEnum); + _printf_("call computational core:\n"); - diagnostic_core_linear(&h_g,fem_p,BalancedthicknessAnalysisEnum,NoneAnalysisEnum); + solver_linear(NULL,femmodel); _printf_("extrude computed thickness on all layers:\n"); - FieldExtrudex( h_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0); + ExtrudeInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,ThicknessEnum); - /*Plug results into output dataset: */ - InputToResultx(&result,fem_p->elements,fem_p->nodes,fem_p->vertices, fem_p->loads, fem_p->materials,fem_p->parameters,ThicknessEnum,results->Size()+1,0,1); results->AddObject(result); - results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(BalancedthicknessAnalysisEnum))); - - /*Free ressources:*/ - VecFree(&u_g); - VecFree(&h_g); - - /*return: */ - return results; + if(verbose)_printf_("saving results:\n"); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); } Index: /issm/trunk/src/c/solutions/balancedvelocities.cpp =================================================================== --- /issm/trunk/src/c/solutions/balancedvelocities.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/balancedvelocities.cpp (revision 4055) @@ -23,16 +23,10 @@ char* inputfilename=NULL; char* outputfilename=NULL; + bool qmu_analysis; char* lockname=NULL; - int numberofnodes; bool waitonlock=false; - Model* model=NULL; - - bool qmu_analysis; - - /*Results: */ - Results* results=NULL; - Result* result=NULL; - Param* param=NULL; + /*FemModel: */ + FemModel* femmodel=NULL; /*time*/ @@ -40,4 +34,7 @@ double start_core, finish_core; double start_init, finish_init; + + int analyses[1]={BalancedvelocitiesAnalysisEnum}; + int solution_type=BalancedvelocitiesAnalysisEnum; MODULEBOOT(); @@ -66,13 +63,14 @@ fid=pfopen(inputfilename,"rb"); - /*Initialize model structure: */ - model=new Model(); + _printf_("create finite element model:\n"); + femmodel=new FemModel(fid,solution_type,analyses,1); - _printf_("read and create finite element model:\n"); - model->AddFormulation(fid,BalancedvelocitiesAnalysisEnum); - - /*recover parameters: */ - model->FindParam(&waitonlock,WaitOnLockEnum); - model->FindParam(&qmu_analysis,QmuAnalysisEnum); + /*add outputfilename in parameters: */ + femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename); + + /*get parameters: */ + femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum); + femmodel->parameters->FindParam(&control_analysis,ControlAnalysisEnum); + femmodel->parameters->FindParam(&waitonlock,WaitOnLockEnum); MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime(); @@ -84,6 +82,9 @@ _printf_("call computational core:\n"); MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); - results=balancedvelocities_core(model); + balancedvelocities_core(femmodel); MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + + _printf_("write results to disk:\n"); + OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters, outputfilename); } @@ -95,5 +96,5 @@ #ifdef _HAVE_DAKOTA_ MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); - Qmux(model,BalancedvelocitiesAnalysisEnum,NoneAnalysisEnum); + Qmux(femmodel,BalancedvelocitiesAnalysisEnum,NoneAnalysisEnum); MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); #else @@ -101,7 +102,4 @@ #endif } - - _printf_("write results to disk:\n"); - OutputResults(results,outputfilename); if (waitonlock>0){ @@ -111,6 +109,5 @@ /*Free ressources:*/ - delete results; - delete model; + delete femmodel; /*Get finish time and close*/ Index: /issm/trunk/src/c/solutions/balancedvelocities_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/balancedvelocities_core.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/balancedvelocities_core.cpp (revision 4055) @@ -11,60 +11,32 @@ #include "../solvers/solvers.h" -Results* balancedvelocities_core(Model* model){ - - extern int my_rank; - - /*output: */ - Results* results=NULL; - Result* result=NULL; - - /*intermediary: */ - Vec u_g=NULL; - - /*solutions: */ - Vec v_g=NULL; +void balancedvelocities_core(FemModel* femmodel){ /*flags: */ int verbose=0; - int numberofdofspernode; - int numberofnodes; - int dofs[2]={1,1}; + int dim; - /*fem balancedvelocities model: */ - FemModel* fem_p=NULL; - - //initialize results: - results=new Results(); - - /*recover fem model: */ - fem_p=model->GetFormulation(BalancedvelocitiesAnalysisEnum); - - //first recover parameters common to all solutions - model->FindParam(&verbose,VerboseEnum); - model->FindParam(&numberofnodes,NumberOfNodesEnum); - model->FindParam(&numberofdofspernode,NumberOfDofsPerNodeEnum); + /*activate formulation: */ + femmodel->SetCurrentAnalysis(BalancedvelocitiesAnalysisEnum); + + /*recover parameters: */ + femmodel->parameters->FindParam(&verbose,VerboseEnum); + femmodel->parameters->FindParam(&dim,DimEnum); _printf_("depth averaging velocity...\n"); - //u_g=inputs->Get("velocity",&dofs[0],2); //take (vx,vy) from inputs velocity - //FieldDepthAveragex( u_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"velocity"); - //inputs->Add("velocity_average",u_g,2,numberofnodes); - ISSMERROR(" not supported yet!"); - + DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VxEnum); + DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VyEnum); + if(dim==3) DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VzEnum); + _printf_("call computational core:\n"); - diagnostic_core_linear(&v_g,fem_p,BalancedvelocitiesAnalysisEnum,NoneAnalysisEnum); + solver_linear(NULL,femmodel); - _printf_("extrude computed thickness on all layers:\n"); - FieldExtrudex( v_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0); + _printf_("extrude computed velocity on all layers:\n"); + ExtrudeInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VxEnum); + ExtrudeInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VyEnum); - /*Plug results into output dataset: */ - InputToResultx(&result,fem_p->elements,fem_p->nodes,fem_p->vertices, fem_p->loads, fem_p->materials,fem_p->parameters,VxEnum,results->Size()+1,0,1); results->AddObject(result); - InputToResultx(&result,fem_p->elements,fem_p->nodes,fem_p->vertices, fem_p->loads, fem_p->materials,fem_p->parameters,VyEnum,results->Size()+1,0,1); results->AddObject(result); - results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(BalancedvelocitiesAnalysisEnum))); + if(verbose)_printf_("saving results:\n"); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum); - /*Free ressources:*/ - VecFree(&u_g); - VecFree(&v_g); - - /*return: */ - return results; } Index: /issm/trunk/src/c/solutions/bedslope.cpp =================================================================== --- /issm/trunk/src/c/solutions/bedslope.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/bedslope.cpp (revision 4055) @@ -37,6 +37,6 @@ double start_init, finish_init; - int analyses[1]={BedSlopeComputeAnalysisEnum}; - int solution_type=BedSlopeComputeSolutionEnum; + int analyses[1]={SlopeAnalysisEnum}; + int solution_type=SlopeAnalysisEnum; MODULEBOOT(); @@ -67,4 +67,7 @@ _printf_("create finite element model, using analyses types statically defined above:\n"); femmodel=new FemModel(fid,solution_type,analyses,1); + + /*add outputfilename in parameters: */ + femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename)); /*get parameters: */ @@ -83,5 +86,5 @@ _printf_("write results to disk:\n"); - OutputResults(femmodel,outputfilename,DiagnosticAnalysisEnum); + OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters); } else{ @@ -92,5 +95,5 @@ #ifdef _HAVE_DAKOTA_ MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); - Qmux(femmodel,SlopeComputeAnalysisEnum,NoneAnalysisEnum); + Qmux(femmodel,BedSlopeAnalysisEnum,NoneAnalysisEnum); MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); #else Index: /issm/trunk/src/c/solutions/bedslope_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/bedslope_core.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/bedslope_core.cpp (revision 4055) @@ -25,12 +25,19 @@ /*Call on core computations: */ - solver_linear(NULL,femmodel,SlopeComputeAnalysisEnum,BedSlopeXAnalysisEnum); - solver_linear(NULL,femmodel,SlopeComputeAnalysisEnum,BedSlopeYAnalysisEnum); + femmodel->SetCurrentAnalysisAlias(SlopeAnalysisEnum,BedSlopeXAnalysisEnum); + solver_linear(NULL,femmodel); + femmodel->SetCurrentAnalysisAlias(SlopeAnalysisEnum,BedSlopeYAnalysisEnum); + solver_linear(NULL,femmodel); /*extrude inputs 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); + if(verbose)_printf_("%s\n","extruding bed in 3d..."); + ExtrudeInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,BedSlopeXEnum); + ExtrudeInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,BedSlopeYEnum); } + + if(verbose)_printf_("saving results:\n"); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedSlopeXEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedSlopeYEnum); + } Index: /issm/trunk/src/c/solutions/control_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/control_core.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/control_core.cpp (revision 4055) @@ -104,5 +104,5 @@ if (((n+1)%5)==0){ _printf_("%s"," saving temporary results..."); - controlrestart(femmodel); + controlrestart(femmodel,J); } } Index: /issm/trunk/src/c/solutions/controlrestart.cpp =================================================================== --- /issm/trunk/src/c/solutions/controlrestart.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/controlrestart.cpp (revision 4055) @@ -7,20 +7,20 @@ #include "../EnumDefinitions/EnumDefinitions.h" -void controlrestart(FemModel* femmodel){ +void controlrestart(FemModel* femmodel,double* J){ - char* outputfilename=NULL; - + int control_type; + int nsteps; + /*retrieve output file name: */ - model->FindParam(&outputfilename,OutputFileNameEnum); + femmodel->parameters->FindParam(&control_type,ControlTypeEnum); + femmodel->parameters->FindParam(&nsteps,NStepsEnum); /*we essentially want J and the parameter: */ InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type); //the parameter itself! - femmodel->otherresults->AddObject(new DoubleResult(femmodel->otherresults->Size()+1,0,1,"J",J,nsteps)); - femmodel->otherresults->AddObject(new StringResult(results->Size()+1,ControlTypeEnum,0,1,EnumAsString(control_type))); + femmodel->results->AddObject(new DoubleVecExternalResult(femmodel->results->Size()+1,JEnum,J,nsteps,1,0)); + femmodel->results->AddObject(new StringExternalResult(femmodel->results->Size()+1,ControlTypeEnum,EnumAsString(control_type),1,0)); /*write to disk: */ - OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters, outputfilename); - - /*Free ressources:*/ - xfree((void**)&outputfilename); + OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters); + } Index: /issm/trunk/src/c/solutions/diagnostic.cpp =================================================================== --- /issm/trunk/src/c/solutions/diagnostic.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/diagnostic.cpp (revision 4055) @@ -36,5 +36,5 @@ double start_init, finish_init; - int analyses[5]={DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticStokesAnalysisEnum,DiagnosticHutterAnalysisEnum,SlopeComputeAnalysisEnum}; + int analyses[5]={DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticStokesAnalysisEnum,DiagnosticHutterAnalysisEnum,SlopeAnalysisEnum}; int solution_type=DiagnosticAnalysisEnum; @@ -58,5 +58,4 @@ lockname=argv[4]; - /*Initialize femmodel structure: */ MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime(); @@ -66,4 +65,7 @@ _printf_("create finite element model:\n"); femmodel=new FemModel(fid,solution_type,analyses,5); + + /*add outputfilename in parameters: */ + femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename)); /*get parameters: */ @@ -94,5 +96,5 @@ _printf_("write results to disk:\n"); - OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters, outputfilename); + OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters); } else{ Index: /issm/trunk/src/c/solutions/diagnostic_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/diagnostic_core.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/diagnostic_core.cpp (revision 4055) @@ -22,6 +22,8 @@ bool isstokes=false; double stokesreconditioning; + bool conserve_loads=true; + bool modify_loads=true; - /* recover parameters: {{{1*/ + /* recover parameters:*/ femmodel->parameters->FindParam(&verbose,VerboseEnum); femmodel->parameters->FindParam(&dim,DimEnum); @@ -31,5 +33,4 @@ femmodel->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum); - /*}}}*/ /*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.*/ @@ -48,5 +49,6 @@ if(verbose)_printf_("%s\n"," computing hutter velocities..."); - solver_linear(NULL,femmodel,DiagnosticAnalysisEnum,HutterAnalysisEnum); + femmodel->SetCurrentAnalysis(DiagnosticHutterAnalysisEnum); + solver_linear(NULL,femmodel); if (ismacayealpattyn) ResetBoundaryConditions(femmodel,DiagnosticAnalysisEnum,HorizAnalysisEnum); @@ -56,5 +58,6 @@ if(verbose)_printf_("%s\n"," computing horizontal velocities..."); - solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,false,DiagnosticAnalysisEnum,HorizAnalysisEnum); //false means we expose loads to changes inside the solution + femmodel->SetCurrentAnalysis(DiagnosticHorizAnalysisEnum); + solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,modify_loads); } @@ -63,5 +66,6 @@ if(verbose)_printf_("%s\n"," computing vertical velocities..."); - solver_linear(NULL,femmodel,DiagnosticAnalysisEnum,VertAnalysisEnum); + femmodel->SetCurrentAnalysis(DiagnosticVertAnalysisEnum); + solver_linear(NULL,femmodel); if (isstokes){ @@ -75,5 +79,6 @@ if(verbose)_printf_("%s\n"," computing stokes velocities and pressure ..."); - solver_diagnostic_nonlinear(NULL,NULL,NULL,NULL,femmodel,DiagnosticAnalysisEnum,StokesAnalysisEnum); + femmodel->SetCurrentAnalysis(DiagnosticStokesAnalysisEnum); + solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,conserve_loads); } } Index: /issm/trunk/src/c/solutions/gradient_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/gradient_core.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/gradient_core.cpp (revision 4055) @@ -32,10 +32,10 @@ if(verbose)_printf_("%s\n"," compute gradient:"); - Gradjx(&gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, control_type,AdjointEnum); + Gradjx(&gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, control_type); if(verbose)_printf_("%s\n"," retrieve old gradient:"); GetVectorFromInputsx(&old_gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, OldGradientEnum,VertexEnum); - if(control_steady)diagnostic_core(model); + if(control_steady)diagnostic_core(femmodel); if (step>0 && search_scalar==0){ @@ -55,6 +55,6 @@ /*plug back into inputs: */ - UpdateInputsFromVectorx( femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials, femmoel->parameters,gradient,GradientEnum,VertexEnum); - UpdateInputsFromVectorx( femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials, femmoel->parameters,old_gradient,OldGradientEnum,VertexEnum); + UpdateInputsFromVectorx( femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials, femmodel->parameters,gradient,GradientEnum,VertexEnum); + UpdateInputsFromVectorx( femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials, femmodel->parameters,old_gradient,OldGradientEnum,VertexEnum); /*Free ressources and return:*/ Index: /issm/trunk/src/c/solutions/objectivefunctionC.cpp =================================================================== --- /issm/trunk/src/c/solutions/objectivefunctionC.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/objectivefunctionC.cpp (revision 4055) @@ -3,8 +3,5 @@ */ -#include "../modules/modules.h" -#include "./solutions.h" -#include "../EnumDefinitions/EnumDefinitions.h" - +/*include files: {{{1*/ #ifdef HAVE_CONFIG_H #include "config.h" @@ -12,4 +9,14 @@ #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" #endif + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../solvers/solvers.h" +#include "./solutions.h" +#include "../modules/modules.h" +#include "../include/include.h" +/*}}}*/ double objectivefunctionC(double search_scalar,OptArgs* optargs){ @@ -22,5 +29,5 @@ /*parameters: */ FemModel* femmodel=NULL; - int n; + int n; double* optscal=NULL; double* fit=NULL; @@ -29,7 +36,7 @@ int control_type; int analysis_type; - int sub_analysis_type; int control_steady; - + bool conserve_loads=true; + /*Recover finite element model: */ femmodel=optargs->femmodel; @@ -45,5 +52,4 @@ femmodel->parameters->FindParam(&control_steady,ControlSteadyEnum); femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum); - femmodel->parameters->FindParam(&sub_analysis_type,SubAnalysisTypeEnum); /*Use ControlParameterEnum input to reinitialize our input parameter: */ @@ -57,10 +63,10 @@ /*Run diagnostic with updated inputs: */ - if(!control_steady) solver_diagnostic_nonlinear(NULL,NULL,NULL,true,femmodel,DiagnosticAnalysisEnum,sub_analysis_type); //true means we conserve loads at each diagnostic run + if(!control_steady) solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,conserve_loads); //true means we conserve loads at each diagnostic run else diagnostic_core(femmodel); //We need a 3D velocity!! (vz is required for the next thermal run) /*Compute misfit for this velocity field.*/ UpdateInputsFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,fit[n],FitEnum); - CostFunctionx( &J, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters,analysis_type,sub_analysis_type); + CostFunctionx( &J, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); /*Free ressources:*/ Index: /issm/trunk/src/c/solutions/prognostic.cpp =================================================================== --- /issm/trunk/src/c/solutions/prognostic.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/prognostic.cpp (revision 4055) @@ -24,15 +24,9 @@ char* outputfilename=NULL; char* lockname=NULL; - int numberofnodes; + bool qmu_analysis=false; bool waitonlock=false; - Model* model=NULL; - - bool qmu_analysis; - - /*Results: */ - Results* results=NULL; - - Param* param=NULL; + /*FemModel: */ + FemModel* femmodel=NULL; /*time*/ @@ -40,4 +34,7 @@ double start_core, finish_core; double start_init, finish_init; + + int analyses[1]={PrognosticAnalysisEnum}; + int solution_type=PrognosticAnalysisEnum; MODULEBOOT(); @@ -62,15 +59,18 @@ /*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"); - _printf_("read and create finite element model:\n"); - model->AddFormulation(fid,PrognosticAnalysisEnum); + _printf_("create finite element model:\n"); + femmodel=new FemModel(fid,solution_type,analyses,1); - /*recover parameters: */ - model->FindParam(&waitonlock,WaitOnLockEnum); - model->FindParam(&qmu_analysis,QmuAnalysisEnum); + /*add outputfilename in parameters: */ + femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename)); + + /*get parameters: */ + femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum); + femmodel->parameters->FindParam(&waitonlock,WaitOnLockEnum); + MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime(); @@ -82,6 +82,9 @@ _printf_("call computational core:\n"); MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); - results=prognostic_core(model); + prognostic_core(femmodel); MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + + _printf_("write results to disk:\n"); + OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters); } @@ -93,5 +96,5 @@ #ifdef _HAVE_DAKOTA_ MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); - Qmux(model,PrognosticAnalysisEnum,NoneAnalysisEnum); + Qmux(femmodel,PrognosticAnalysisEnum,NoneAnalysisEnum); MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); #else @@ -101,7 +104,4 @@ - _printf_("write results to disk:\n"); - OutputResults(results,outputfilename); - if (waitonlock>0){ _printf_("write lock file:\n"); @@ -110,6 +110,5 @@ /*Free ressources:*/ - delete results; - delete model; + delete femmodel; /*Get finish time and close*/ Index: /issm/trunk/src/c/solutions/prognostic2.cpp =================================================================== --- /issm/trunk/src/c/solutions/prognostic2.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/prognostic2.cpp (revision 4055) @@ -24,15 +24,9 @@ char* outputfilename=NULL; char* lockname=NULL; - int numberofnodes; + bool qmu_analysis=false; bool waitonlock=false; - Model* model=NULL; - - bool qmu_analysis; - - /*Results: */ - Results* results=NULL; - - Param* param=NULL; + /*FemModel: */ + FemModel* femmodel=NULL; /*time*/ @@ -40,4 +34,7 @@ double start_core, finish_core; double start_init, finish_init; + + int analyses[1]={Prognostic2AnalysisEnum}; + int solution_type=Prognosti2cAnalysisEnum; MODULEBOOT(); @@ -62,16 +59,18 @@ /*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"); - _printf_("read and create finite element model:\n"); - model->AddFormulation(fid,Prognostic2AnalysisEnum); + _printf_("create finite element model:\n"); + femmodel=new FemModel(fid,solution_type,analyses,1); + + /*add outputfilename in parameters: */ + femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename)); /*recover parameters: */ model->FindParam(&waitonlock,WaitOnLockEnum); model->FindParam(&qmu_analysis,QmuAnalysisEnum); - model->FindParam(&numberofnodes,NumberOfNodesEnum); + MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime(); @@ -83,6 +82,9 @@ _printf_("call computational core:\n"); MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); - results=prognostic2_core(model); + prognostic2_core(femmodel); MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + + _printf_("write results to disk:\n"); + OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters); } @@ -94,5 +96,5 @@ #ifdef _HAVE_DAKOTA_ MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); - Qmux(model,Prognostic2AnalysisEnum,NoneAnalysisEnum); + Qmux(femmodel,Prognostic2AnalysisEnum,NoneAnalysisEnum); MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); #else @@ -100,7 +102,4 @@ #endif } - - _printf_("write results to disk:\n"); - OutputResults(results,outputfilename); if (waitonlock>0){ @@ -110,6 +109,5 @@ /*Free ressources:*/ - delete results; - delete model; + delete femmodel; /*Get finish time and close*/ Index: /issm/trunk/src/c/solutions/prognostic2_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/prognostic2_core.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/prognostic2_core.cpp (revision 4055) @@ -11,64 +11,29 @@ #include "../solvers/solvers.h" -Results* prognostic2_core(Model* model){ - - extern int my_rank; - - /*output: */ - Results* results=NULL; - Result* result=NULL; - - /*intermediary: */ - Vec vx_g=NULL; - Vec vy_g=NULL; - - /*solutions: */ - Vec h_g=NULL; +void prognostic2_core(FemModel* femmodel){ /*flags: */ int verbose=0; - int numberofdofspernode; - int numberofnodes; - int numberofvertices; - int dofs[1]={1}; - /*fem prognostic model: */ - FemModel* fem_p=NULL; + /*activate formulation: */ + femmodel->SetCurrentAnalysis(Prognostic2AnalysisEnum); - //initialize results: - results=new Results(); - - /*recover fem model: */ - fem_p=model->GetFormulation(Prognostic2AnalysisEnum); - - //first recover parameters common to all solutions - model->FindParam(&verbose,VerboseEnum); - model->FindParam(&numberofnodes,NumberOfNodesEnum); - model->FindParam(&numberofvertices,NumberOfVerticesEnum); - model->FindParam(&numberofdofspernode,NumberOfNodesEnum); + /*recover parameters: */ + femmodel->parameters->FindParam(&verbose,VerboseEnum); _printf_("depth averaging velocity...\n"); /*Where is it done?*/ + + _printf_("call computational core:\n"); + solver_linear(NULL,femmodel); - _printf_("call computational core:\n"); - diagnostic_core_linear(&h_g,fem_p,Prognostic2AnalysisEnum,NoneAnalysisEnum); - _printf_("Averaging over vertices:\n"); - FieldAverageOntoVerticesx(&h_g,fem_p->elements,fem_p->nodes,fem_p->vertices,fem_p->loads,fem_p->materials,fem_p->parameters); + //FieldAverageOntoVerticesx(&h_g,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); //_printf_("extrude computed thickness on all layers:\n"); //FieldExtrudex(h_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0); - /*Plug results into output dataset: */ - InputToResultx(&result,fem_p->elements,fem_p->nodes,fem_p->vertices, fem_p->loads, fem_p->materials,fem_p->parameters,ThicknessEnum,results->Size()+1,0,1); results->AddObject(result); - results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(Prognostic2AnalysisEnum))); - - /*Free ressources:*/ - VecFree(&vx_g); - VecFree(&vy_g); - VecFree(&h_g); - - /*return: */ - return results; + if(verbose)_printf_("saving results:\n"); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); } Index: /issm/trunk/src/c/solutions/prognostic_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/prognostic_core.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/prognostic_core.cpp (revision 4055) @@ -11,52 +11,28 @@ #include "../solvers/solvers.h" -Results* prognostic_core(Model* model){ +void prognostic_core(FemModel* femmodel){ - extern int my_rank; - - /*output: */ - Results* results=NULL; - Result* result=NULL; - - /*solutions: */ - Vec h_g=NULL; - - /*flags: */ + /*parameters: */ int verbose=0; - /*fem prognostic model: */ - FemModel* fem_p=NULL; - - //initialize results - results=new Results(); - - /*recover fem model: */ - fem_p=model->GetFormulation(PrognosticAnalysisEnum); ISSMASSERT(fem_p); - - //first recover parameters common to all solutions - model->FindParam(&verbose,VerboseEnum); + /*activate formulation: */ + femmodel->SetCurrentAnalysis(PrognosticAnalysisEnum); + + /*recover parameters: */ + femmodel->parameters->FindParam(&verbose,VerboseEnum); _printf_("depth averaging velocity...\n"); - DepthAverageInputx(fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,VxEnum); - DepthAverageInputx(fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,VyEnum); + DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VxEnum,VxAverageEnum); + DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VyEnum,VyAverageEnum); _printf_("call computational core:\n"); - diagnostic_core_linear(&h_g,fem_p,PrognosticAnalysisEnum,NoneAnalysisEnum); + solver_linear(NULL,femmodel); - _printf_("update inputs:\n"); - UpdateInputsFromSolutionx( fem_p->elements,fem_p->nodes, fem_p->vertices, fem_p->loads, fem_p->materials, fem_p->parameters,h_g,PrognosticAnalysisEnum,NoneAnalysisEnum); + _printf_("extrude computed thickness on all layers:\n"); + ExtrudeInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,ThicknessEnum); - _printf_("extrude computed thickness on all layers:\n"); - ExtrudeInputx(fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,ThicknessEnum); + if(verbose)_printf_("saving results:\n"); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); - _printf_("extract result from extruded inputs: \n"); - InputToResultx(&result,fem_p->elements,fem_p->nodes,fem_p->vertices, fem_p->loads, fem_p->materials,fem_p->parameters,ThicknessEnum,results->Size()+1,0,1); results->AddObject(result); - results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(PrognosticAnalysisEnum))); +} - /*Free ressources:*/ - VecFree(&h_g); - - //return: - return results; - -} Index: /issm/trunk/src/c/solutions/solutions.h =================================================================== --- /issm/trunk/src/c/solutions/solutions.h (revision 4054) +++ /issm/trunk/src/c/solutions/solutions.h (revision 4055) @@ -13,44 +13,37 @@ /*cores: */ -Results* prognostic_core(FemModel* femmodel); -Results* prognostic2_core(FemModel* femmodel); -Results* balancedthickness_core(FemModel* femmodel); -Results* balancedthickness2_core(FemModel* femmodel); -Results* balancedvelocities_core(FemModel* femmodel); -Results* slopecompute_core(FemModel* femmodel); -Results* steadystate_core(FemModel* femmodel); -Results* transient_core(FemModel* femmodel); -Results* transient_core_2d(FemModel* femmodel); -Results* transient_core_3d(FemModel* femmodel); void adjoint_core(FemModel* femmodel); void gradient_core(FemModel* femmodel,int step=0, double search_scalar=0); void diagnostic_core(FemModel* femmodel); void thermal_core(FemModel* femmodel); -void surfaceslope_core(FemModel* femfemmodel); -void bedslope_core(FemModel* femfemmodel); +void thermal_core_step(FemModel* femmodel,int step, double time); +void surfaceslope_core(FemModel* femmodel); +void bedslope_core(FemModel* femmodel); void control_core(FemModel* femmodel); +void prognostic_core(FemModel* femmodel); +void prognostic2_core(FemModel* femmodel); +void balancedthickness_core(FemModel* femmodel); +void balancedthickness2_core(FemModel* femmodel); +void balancedvelocities_core(FemModel* femmodel); +void slopecompute_core(FemModel* femmodel); +void steadystate_core(FemModel* femmodel); +void transient2d_core(FemModel* femmodel); +void transient3d_core(FemModel* femmodel); +double objectivefunctionC(double search_scalar,OptArgs* optargs); -//int GradJOrth(WorkspaceParams* workspaceparams); +//convergence: void convergence(int* pconverged, Mat K_ff,Vec p_f,Vec u_f,Vec u_f_old,Parameters* parameters); int controlconvergence(double* J, double* fit, double eps_cm, int n); -int GoldenSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*),FemModel* femfemmodel); +//optimization +int GoldenSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*),FemModel* femmodel); +int BrentSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*),FemModel* femmodel); +int GradJSearch(double* search_vector,FemModel* femmodel,int step); -int BrentSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*),FemModel* femfemmodel); - -double objectivefunctionC(double search_scalar,OptArgs* optargs); - -int GradJSearch(double* search_vector,FemModel* femfemmodel,int step); -//int GradJCheck(WorkspaceParams* workspaceparams,int step,int status); - -//int ParameterUpdate(double* search_vector,int step, WorkspaceParams* workspaceparams,BatchParams* batchparams); +//diverse void WriteLockFile(char* filename); - -void stokescontrolinit(FemModel* femfemmodel); -void controlrestart(FemModel* femfemmodel); - -void CreateFemModel(FemModel* femfemmodel,ConstDataHandle MODEL,int analysis_type,int sub_analysis_type); -//int BatchDebug(Mat* Kgg,Vec* pg,FemModel* femfemmodel,char* filename); -void ResetBoundaryConditions(FemModel* femfemmodel, int analysis_type, int sub_analysis_type); +void stokescontrolinit(FemModel* femmodel); +void controlrestart(FemModel* femmodel,double* J); +void ResetBoundaryConditions(FemModel* femmodel, int analysis_type, int sub_analysis_type); #endif Index: /issm/trunk/src/c/solutions/steadystate.cpp =================================================================== --- /issm/trunk/src/c/solutions/steadystate.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/steadystate.cpp (revision 4055) @@ -24,25 +24,10 @@ char* outputfilename=NULL; char* lockname=NULL; - int numberofnodes; bool qmu_analysis=false; bool control_analysis=false; - char* control_type=NULL; + bool waitonlock=false; - /*Model: */ - Model* model=NULL; - FemModel* fem_dh=NULL; - FemModel* fem_ds=NULL; - - /*Results: */ - Results* results=NULL; - - bool waitonlock=false; - - double* u_g_initial=NULL; - double* p_g_initial=NULL; - double* u_g_obs=NULL; - double* weights=NULL; - double dt; - BoolParam* param=NULL; + /*FemModel: */ + FemModel* femmodel=NULL; /*time*/ @@ -50,4 +35,10 @@ double start_core, finish_core; double start_init, finish_init; + + /*intermediary: */ + BoolParam* param=NULL; + + int analyses[7]={DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticStokesAnalysisEnum,DiagnosticHutterAnalysisEnum,SlopeAnalysisEnum,ThermalAnalysisEnum,MeltingAnalysisEnum}; + int solution_type=SteadystateAnalysisEnum; MODULEBOOT(); @@ -70,32 +61,15 @@ lockname=argv[4]; + /*Initialize model structure: */ + MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime(); + /*Open handle to data on disk: */ fid=pfopen(inputfilename,"rb"); - /*Initialize model structure: */ - MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime(); - model=new Model(); + _printf_("create finite element model:\n"); + femmodel=new FemModel(fid,solution_type,analyses,5); - _printf_("read and create finite element model:\n"); - - _printf_("\n reading diagnostic horiz model data:\n"); - model->AddFormulation(fid,DiagnosticAnalysisEnum,HorizAnalysisEnum); - - _printf_("\n reading diagnostic vert model data:\n"); - model->AddFormulation(fid,DiagnosticAnalysisEnum,VertAnalysisEnum); - - _printf_("\n reading diagnostic stokes model data:\n"); - model->AddFormulation(fid,DiagnosticAnalysisEnum,StokesAnalysisEnum); - - _printf_("\n reading diagnostic hutter model data:\n"); - model->AddFormulation(fid,DiagnosticAnalysisEnum,HutterAnalysisEnum); - - _printf_("\n reading surface and bed slope computation model data:\n"); - model->AddFormulation(fid,SlopecomputeAnalysisEnum); - - _printf_("\n read and create thermal finite element model:\n"); - model->AddFormulation(fid,ThermalAnalysisEnum); - _printf_("\n read and create melting finite element model:\n"); - model->AddFormulation(fid,MeltingAnalysisEnum); + /*add outputfilename in parameters: */ + femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename)); /*recover parameters: */ @@ -113,5 +87,5 @@ _printf_("call computational core:\n"); MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); - results=steadystate_core(model); + steadystate_core(femmodel); MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); @@ -119,20 +93,16 @@ else{ /*change control_steady to 1 to know we are doing steadystate*/ - fem_dh=model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum); - fem_ds=model->GetFormulation(DiagnosticAnalysisEnum,StokesAnalysisEnum); - - param=(BoolParam*)fem_dh->parameters->FindParamObject(ControlSteadyEnum); param->value=true; - param=(BoolParam*)fem_ds->parameters->FindParamObject(ControlSteadyEnum); param->value=true; + param=(BoolParam*)femmodel->parameters->FindParamObject(ControlSteadyEnum); param->value=true; /*run control analysis: */ _printf_("call computational core:\n"); MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); - results=control_core(model); + control_core(femmodel); MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); - } _printf_("write results to disk:\n"); - OutputResults(results,outputfilename); + OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters); + } else{ @@ -142,5 +112,5 @@ #ifdef _HAVE_DAKOTA_ MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); - Qmux(model,SteadystateAnalysisEnum,NoneAnalysisEnum); + Qmux(femmodel,SteadystateAnalysisEnum,NoneAnalysisEnum); MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); #else @@ -155,9 +125,5 @@ /*Free ressources */ - xfree((void**)&u_g_initial); - xfree((void**)&u_g_obs); - xfree((void**)&p_g_initial); - delete model; - delete results; + delete femmodel; /*Get finish time and close*/ Index: /issm/trunk/src/c/solutions/steadystate_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/steadystate_core.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/steadystate_core.cpp (revision 4055) @@ -9,136 +9,54 @@ #include "./solutions.h" #include "../modules/modules.h" +#include "../include/include.h" #include "../solvers/solvers.h" -Results* steadystate_core(Model* model){ +void steadystate_core(FemModel* femmodel){ - extern int my_rank; + /*intermediary: */ + int step; - /*fem models: */ - FemModel* fem_dh=NULL; - FemModel* fem_dv=NULL; - FemModel* fem_dhu=NULL; - FemModel* fem_ds=NULL; - FemModel* fem_sl=NULL; - FemModel* fem_t=NULL; - FemModel* fem_m=NULL; + /*parameters: */ + int verbose; + int dim; + + /* recover parameters:*/ + femmodel->parameters->FindParam(&verbose,VerboseEnum); + femmodel->parameters->FindParam(&dim,DimEnum); - /*output: */ - Results* results=NULL; - Results* results_thermal=NULL; - Results* results_diagnostic=NULL; - - /*solutions: */ - Vec u_g=NULL; - Vec old_u_g=NULL; - Vec t_g=NULL; - Vec t_g_average=NULL; - Vec old_t_g=NULL; - Vec p_g=NULL; - Vec m_g=NULL; - Vec du_g=NULL; - Vec dt_g=NULL; - double ndu,nu; - double normdt,normt; - double eps_rel; - - /*flags: */ - int verbose=0; - int isstokes=0; - int numberofnodes; - int ndof; - int converged; - int step; - - /*recover fem models: */ - fem_dh=model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum); - fem_dv=model->GetFormulation(DiagnosticAnalysisEnum,VertAnalysisEnum); - fem_ds=model->GetFormulation(DiagnosticAnalysisEnum,StokesAnalysisEnum); - fem_dhu=model->GetFormulation(DiagnosticAnalysisEnum,HutterAnalysisEnum); - fem_sl=model->GetFormulation(SlopecomputeAnalysisEnum); - fem_t=model->GetFormulation(ThermalAnalysisEnum); - fem_m=model->GetFormulation(MeltingAnalysisEnum); - - - //first recover parameters common to all solutions - model->FindParam(&verbose,VerboseEnum); - model->FindParam(&numberofnodes,NumberOfNodesEnum); - model->FindParam(&eps_rel,EpsResEnum); - model->FindParam(&isstokes,IsStokesEnum); - - //initialize: - converged=0; + /*intialize counters: */ step=1; - - if (isstokes)ndof=4; - else ndof=3; for(;;){ if(verbose)_printf_("%s%i\n"," computing temperature and velocity for step: ",step); + thermal_core(femmodel); - //first compute temperature at steady state. - results_thermal=thermal_core(model); + if(verbose)_printf_("%s\n","computing depth average temperature"); + DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,TemperatureEnum,TemperatureAverageEnum); - //get t_g and m_g; - VecFree(&t_g);results_thermal->FindResult(&t_g,"t_g"); - VecFree(&m_g);results_thermal->FindResult(&m_g,"m_g"); - delete results_thermal; + if(verbose)_printf_("%s\n","computing new velocity"); + diagnostic_core(femmodel); - //Add temperature to inputs. - //compute depth averaged temperature and add to inputs - VecDuplicatePatch(&t_g_average,t_g); - FieldDepthAveragex( t_g_average, fem_t->elements,fem_t->nodes, fem_t->vertices,fem_t->loads, fem_t->materials,fem_t->parameters,"temperature"); - model->UpdateInputsFromVector(t_g_average,TemperatureAverageEnum,VertexEnum); - model->UpdateInputsFromVector(t_g,TemperatureEnum,VertexEnum); - VecFree(&t_g_average); //not needed anymore - - //now compute diagnostic velocity using the steady state temperature. - results_diagnostic=diagnostic_core(model); - - //get p_g and u_g - VecFree(&u_g);results_diagnostic->FindResult(&u_g,"u_g"); - VecFree(&p_g);results_diagnostic->FindResult(&p_g,"p_g"); - delete results_diagnostic; - - //convergence? - if(step>1){ - VecDuplicatePatch(&du_g,old_u_g);VecAYPX(du_g,-1.0,u_g); - VecNorm(du_g,NORM_2,&ndu); VecNorm(old_u_g,NORM_2,&nu); VecFree(&du_g); - - VecDuplicatePatch(&dt_g,old_t_g); VecAYPX(dt_g,-1.0,t_g); - VecNorm(dt_g,NORM_2,&normdt); VecNorm(old_t_g,NORM_2,&normt);VecFree(&dt_g); - - if (verbose) _printf_("%-60s%g\n %s%g\n %s%g%s\n", - " relative convergence criterion: velocity -> norm(du)/norm(u)= ",ndu/nu*100," temperature -> norm(dt)/norm(t)=",normdt/normt*100," eps_rel: ",eps_rel*100," %"); + if(verbose)_printf_("%s\n","checking velocity, temperature and pressure convergence"); + if (step>1) if(steadystateconvergence(femmodel))goto cleanup_and_return; - if ((ndu/nu<=eps_rel) && (normdt/normt<=eps_rel)) converged=1; - else converged=0; - } - else{ - converged=0; - } - - VecFree(&old_u_g);VecDuplicatePatch(&old_u_g,u_g); - VecFree(&old_t_g);VecDuplicatePatch(&old_t_g,t_g); - + if(verbose)_printf_("%s\n","saving velocity, temperature and pressure to check for convergence at next step"); + DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,OldVxEnum); + DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,OldVyEnum); + if(dim==3)DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum,OldVzEnum); + DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,OldPressureEnum); + DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,OldTemperatureEnum); + + //increase counter step++; - if (converged)break; } - /*Plug results into output dataset: */ - results->AddObject(new Result(results->Size()+1,0,1,"u_g",u_g)); - results->AddObject(new Result(results->Size()+1,0,1,"p_g",p_g)); - results->AddObject(new Result(results->Size()+1,0,1,"t_g",t_g)); - results->AddObject(new Result(results->Size()+1,0,1,"m_g",m_g)); - results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(SteadyAnalysisEnum))); - - - /*Free ressource*/ - VecFree(&old_u_g); - VecFree(&old_t_g); - VecFree(&u_g); - VecFree(&p_g); - VecFree(&t_g); - VecFree(&m_g); + cleanup_and_return: + + if(verbose)_printf_("saving results:\n"); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum); + if(dim==3) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); } Index: /issm/trunk/src/c/solutions/steadystateconvergence.cpp =================================================================== --- /issm/trunk/src/c/solutions/steadystateconvergence.cpp (revision 4055) +++ /issm/trunk/src/c/solutions/steadystateconvergence.cpp (revision 4055) @@ -0,0 +1,36 @@ +/*!\file: steadystateconvergence.cpp + * \brief: determine convergence of steady state solution + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +int steadystateconvergence(FemModel* femmodel){ + + /*output: */ + int converged=0; + int velocity_converged=0; + int temperature_converged=0; + + /*intermediary: */ + int velocityenums[8]={VxEnum,OldVxEnum,VyEnum,OldVyEnum,VzEnum,OldVzEnum,PressureEnum,OldPressureEnum}; //pairs of enums (new and old) on which to carry out the converence tests + int temperatureenums[2]={TemperatureEnum,OldTemperatureEnum}; + int convergencecriterion[1]={RelativeEnum}; //criterions for convergence, RelativeEnum or AbsoluteEnum + int convergencecriterionvalue[1]; //value of criterion to be respected + + /*retrieve parameters: */ + femmodel->parameters->FindParam(&convergencecriterionvalue[0],EpsRelEnum); + + /*figure out convergence at the input level, because we don't have the solution vectors!: */ + velocity_converged=InputConvergencex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,velocityenums,8,convergencecriterion,convergencecriterionvalue,1); + temperature_converged=InputConvergencex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,temperatureenums,2,convergencecriterion,convergencecriterionvalue,1); + + if(velocity_converged && temperature_converged)converged=1; + + /*return: */ + return converged; + +} Index: /issm/trunk/src/c/solutions/stokescontrolinit.cpp =================================================================== --- /issm/trunk/src/c/solutions/stokescontrolinit.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/stokescontrolinit.cpp (revision 4055) @@ -18,15 +18,16 @@ int isstokes=0; double stokesreconditioning; - + bool conserve_loads=true; + /*first recover parameters common to all solutions:*/ - model->FindParam(&verbose,VerboseEnum); - model->FindParam(&isstokes,IsStokesEnum); - model->FindParam(&stokesreconditioning,StokesReconditioningEnum); + femmodel->parameters->FindParam(&verbose,VerboseEnum); + femmodel->parameters->FindParam(&isstokes,IsStokesEnum); + femmodel->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); /*if no Stokes analysis carried out, assign output and return*/ - if (!isstokes)femmodel->SetAnalysisType(DiagnosticAnalysisEnum,DiagnosticHorizAnalysisEnum); + if (!isstokes)femmodel->SetCurrentAnalysis(DiagnosticHorizAnalysisEnum); /* For Stokes inverse control method, we are going to carry out the inversion only on the Stokes part. So we need - * to solve the Hutter or MacAyeal/Pattyn model here, and constrain the Stokes model using the Hutter + * to solve the Hutter or MacAyeal/Pattyn femmodel here, and constrain the Stokes femmodel using the Hutter * or MacAyeal/Pattyn at the boundary. We don't want to have to do that at every inversion step, as * it needs be done only once: */ @@ -36,8 +37,9 @@ /*Run a complete diagnostic to update the Stokes spcs: */ - solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,false,DiagnosticAnalysisEnum,HorizAnalysisEnum); //false means we expose loads to changes inside the solution + femmodel->SetCurrentAnalysis(DiagnosticHorizAnalysisEnum); + solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,conserve_loads); //vertical velocity - solver_linear(NULL,femmodel,DiagnosticAnalysisEnum,VertAnalysisEnum); + solver_linear(NULL,femmodel); //recondition" pressure computed previously: @@ -49,7 +51,5 @@ if(verbose)_printf_("%s\n"," computing stokes velocities and pressure ..."); - solver_diagnostic_nonlinear(NULL,NULL,NULL,NULL,femmodel,DiagnosticAnalysisEnum,StokesAnalysisEnum); - - /*Assign output*/ - model->SetActiveFormulation(fem_ds); + femmodel->SetCurrentAnalysis(DiagnosticStokesAnalysisEnum); + solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,conserve_loads); } Index: /issm/trunk/src/c/solutions/surfaceslope.cpp =================================================================== --- /issm/trunk/src/c/solutions/surfaceslope.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/surfaceslope.cpp (revision 4055) @@ -37,6 +37,6 @@ double start_init, finish_init; - int analyses[1]={SurfaceSlopeComputeAnalysisEnum}; - int solution_type=SurfaceSlopeComputeSolutionEnum; + int analyses[1]={SlopeAnalysisEnum}; + int solution_type=SlopeAnalysisEnum; MODULEBOOT(); @@ -67,4 +67,7 @@ _printf_("create finite element model, using analyses types statically defined above:\n"); femmodel=new FemModel(fid,solution_type,analyses,1); + + /*add outputfilename in parameters: */ + femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename)); /*get parameters: */ @@ -83,5 +86,5 @@ _printf_("write results to disk:\n"); - OutputResults(femmodel,outputfilename,DiagnosticAnalysisEnum); + OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters); } else{ @@ -92,5 +95,5 @@ #ifdef _HAVE_DAKOTA_ MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); - Qmux(femmodel,SlopeComputeAnalysisEnum,NoneAnalysisEnum); + Qmux(femmodel,SurfaceSlopeAnalysisEnum,NoneAnalysisEnum); MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); #else Index: /issm/trunk/src/c/solutions/surfaceslope_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/surfaceslope_core.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/surfaceslope_core.cpp (revision 4055) @@ -25,12 +25,19 @@ /*Call on core computations: */ - solver_linear(NULL,femmodel,SlopeComputeAnalysisEnum,SurfaceSlopeXAnalysisEnum); - solver_linear(NULL,femmodel,SlopeComputeAnalysisEnum,SurfaceSlopeYAnalysisEnum); + femmodel->SetCurrentAnalysisAlias(SlopeAnalysisEnum,SurfaceSlopeXAnalysisEnum); + solver_linear(NULL,femmodel); + femmodel->SetCurrentAnalysisAlias(SlopeAnalysisEnum,SurfaceSlopeYAnalysisEnum); + solver_linear(NULL,femmodel); /*extrude inputs 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); + ExtrudeInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,SurfaceSlopeXEnum); + ExtrudeInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,SurfaceSlopeYEnum); } + + if(verbose)_printf_("saving results:\n"); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceSlopeXEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceSlopeYEnum); + } Index: /issm/trunk/src/c/solutions/thermal.cpp =================================================================== --- /issm/trunk/src/c/solutions/thermal.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/thermal.cpp (revision 4055) @@ -68,4 +68,7 @@ femmodel=new FemModel(fid,solution_type,analyses,2); + /*add outputfilename in parameters: */ + femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename)); + /*get parameters: */ femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum); @@ -85,5 +88,5 @@ _printf_("write results to disk:\n"); - OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters, outputfilename); + OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters); } Index: /issm/trunk/src/c/solutions/thermal_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/thermal_core.cpp (revision 4054) +++ /issm/trunk/src/c/solutions/thermal_core.cpp (revision 4055) @@ -42,9 +42,6 @@ time=(i+1)*dt; - if(verbose)_printf_("computing temperatures:\n"); - solver_thermal_nonlinear(NULL,NULL,femmodel,ThermalAnalysisEnum,NoneAnalysisEnum); - - if(verbose)_printf_("computing melting:\n"); - solver_linear(NULL,femmodel,MeltingAnalysisEnum,NoneAnalysisEnum); + /*call thermal_core_step: */ + thermal_core_step(femmodel,i,time); if(verbose)_printf_("saving results:\n"); Index: /issm/trunk/src/c/solutions/thermal_core_step.cpp =================================================================== --- /issm/trunk/src/c/solutions/thermal_core_step.cpp (revision 4055) +++ /issm/trunk/src/c/solutions/thermal_core_step.cpp (revision 4055) @@ -0,0 +1,26 @@ +/*!\file: thermal_core_step.cpp + * \brief: core of the thermal solution + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./solutions.h" +#include "../modules/modules.h" +#include "../include/include.h" +#include "../solvers/solvers.h" + +void thermal_core_step(FemModel* femmodel,int step, double time){ + + //first recover parameters common to all solutions + femmodel->parameters->FindParam(&verbose,VerboseEnum); + + if(verbose)_printf_("computing temperatures:\n"); + femmodel->SetCurrentAnalysis(ThermalAnalysisEnum); + solver_thermal_nonlinear(NULL,NULL,femmodel); + + if(verbose)_printf_("computing melting:\n"); + femmodel->SetCurrentAnalysis(MeltingAnalysisEnum); + solver_linear(NULL,femmodel); +} Index: sm/trunk/src/c/solutions/transient.cpp =================================================================== --- /issm/trunk/src/c/solutions/transient.cpp (revision 4054) +++ (revision ) @@ -1,150 +1,0 @@ -/*!\file: transient.cpp - * \brief: transient 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 qmu_analysis=false; - bool waitonlock=false; - - /*Model: */ - Model* model=NULL; - int dim=-1; - - /*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]; - - /*Open handle to data on disk: */ - fid=pfopen(inputfilename,"rb"); - - /*Initialize model structure: */ - MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime(); - model=new Model(); - - _printf_("read and create finite element model:\n"); - _printf_("\n reading diagnostic horiz model data:\n"); - model->AddFormulation(fid,DiagnosticAnalysisEnum,HorizAnalysisEnum); - - _printf_("\n reading diagnostic vert model data:\n"); - model->AddFormulation(fid,DiagnosticAnalysisEnum,VertAnalysisEnum); - - _printf_("\n reading diagnostic stokes model data:\n"); - model->AddFormulation(fid,DiagnosticAnalysisEnum,StokesAnalysisEnum); - - _printf_("\n reading diagnostic hutter model data:\n"); - model->AddFormulation(fid,DiagnosticAnalysisEnum,HutterAnalysisEnum); - - _printf_("\n reading surface and bed slope computation model data:\n"); - model->AddFormulation(fid,SlopecomputeAnalysisEnum); - - _printf_("\n reading prognositc model data:\n"); - model->AddFormulation(fid,PrognosticAnalysisEnum); - - /*Do we run in 3d?, in which case we need thermal and melting also:*/ - model->FindParam(&dim,DimEnum); - if(dim==3){ - _printf_("read and create thermal finite element model:\n"); - model->AddFormulation(fid,ThermalAnalysisEnum,TransientAnalysisEnum); - _printf_("read and create melting finite element model:\n"); - model->AddFormulation(fid,MeltingAnalysisEnum,TransientAnalysisEnum); - } - - /*recover parameters: */ - model->FindParam(&waitonlock,WaitOnLockEnum); - model->FindParam(&qmu_analysis,QmuAnalysisEnum); - - MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime(); - - /*are we running the solution sequence, or a qmu wrapper around it? : */ - if(!qmu_analysis){ - - /*run diagnostic analysis: */ - _printf_("call computational core:\n"); - MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); - results=transient_core(model); - MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); - - _printf_("write results to disk:\n"); - OutputResults(results,outputfilename); - } - else{ - /*run qmu analysis: */ - _printf_("calling qmu analysis on transient core:\n"); - - #ifdef _HAVE_DAKOTA_ - MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); - Qmux(model,TransientAnalysisEnum,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 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: /issm/trunk/src/c/solutions/transient2d.cpp =================================================================== --- /issm/trunk/src/c/solutions/transient2d.cpp (revision 4055) +++ /issm/trunk/src/c/solutions/transient2d.cpp (revision 4055) @@ -0,0 +1,120 @@ +/*!\file: transient2d.cpp + * \brief: transient2d 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; + bool waitonlock=false; + + /*FemModel: */ + FemModel* femmodel=NULL; + + /*time*/ + double start, finish; + double start_core, finish_core; + double start_init, finish_init; + + int analyses[2]={DiagnosticHorizAnalysisEnum,PrognosticAnalysisEnum}; + int solution_type=TransientAnalysisEnum; + + 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]; + + /*Open handle to data on disk: */ + fid=pfopen(inputfilename,"rb"); + + _printf_("create finite element model:\n"); + femmodel=new FemModel(fid,solution_type,analyses,7); + + /*add outputfilename in parameters: */ + femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename)); + + /*get parameters: */ + femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum); + femmodel->parameters->FindParam(&control_analysis,ControlAnalysisEnum); + 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){ + + /*run diagnostic analysis: */ + _printf_("call computational core:\n"); + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + transient2d_core(model); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + + _printf_("write results to disk:\n"); + OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters); + } + else{ + /*run qmu analysis: */ + _printf_("calling qmu analysis on transient2d core:\n"); + + #ifdef _HAVE_DAKOTA_ + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + Qmux(femmodel,Transient2DAnalysisEnum,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","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: /issm/trunk/src/c/solutions/transient2d_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/transient2d_core.cpp (revision 4055) +++ /issm/trunk/src/c/solutions/transient2d_core.cpp (revision 4055) @@ -0,0 +1,59 @@ +/*!\file: transient_2d_core.cpp + * \brief: core of the transient_2d solution + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./solutions.h" +#include "../modules/modules.h" + +void transient2d_core(FemModel* femmodel){ + + int i; + + /*parameters: */ + int verbose=0; + double finaltime; + double dt; + + /*intermediary: */ + int step; + double time; + + /* recover parameters: */ + femmodel->parameters->FindParam(&verbose,VerboseEnum); + femmodel->parameters->FindParam(&finaltime,NdtEnum); + femmodel->parameters->FindParam(&dt,DtEnum); + + /*initialize: */ + step=0; + time=0; + + while(timeelements, femmodel->nodes,femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); + + if(verbose)_printf_("%s\n","saving results:\n"); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time); + } + +} Index: /issm/trunk/src/c/solutions/transient3d.cpp =================================================================== --- /issm/trunk/src/c/solutions/transient3d.cpp (revision 4055) +++ /issm/trunk/src/c/solutions/transient3d.cpp (revision 4055) @@ -0,0 +1,120 @@ +/*!\file: transient.cpp + * \brief: transient 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; + bool waitonlock=false; + + /*FemModel: */ + FemModel* femmodel=NULL; + + /*time*/ + double start, finish; + double start_core, finish_core; + double start_init, finish_init; + + int analyses[8]={DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticStokesAnalysisEnum,DiagnosticHutterAnalysisEnum,SlopeAnalysisEnum,PrognosticAnalysisEnum,ThermalAnalysisEnum,MeltingAnalysisEnum}; + int solution_type=TransientAnalysisEnum; + + 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]; + + /*Open handle to data on disk: */ + fid=pfopen(inputfilename,"rb"); + + _printf_("create finite element model:\n"); + femmodel=new FemModel(fid,solution_type,analyses,8); + + /*add outputfilename in parameters: */ + femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename)); + + /*get parameters: */ + femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum); + femmodel->parameters->FindParam(&control_analysis,ControlAnalysisEnum); + 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){ + + /*run diagnostic analysis: */ + _printf_("call computational core:\n"); + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + transient3d_core(model); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + + _printf_("write results to disk:\n"); + OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters); + } + else{ + /*run qmu analysis: */ + _printf_("calling qmu analysis on transient core:\n"); + + #ifdef _HAVE_DAKOTA_ + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + Qmux(femmodel,Transient3DAnalysisEnum,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","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: /issm/trunk/src/c/solutions/transient3d_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/transient3d_core.cpp (revision 4055) +++ /issm/trunk/src/c/solutions/transient3d_core.cpp (revision 4055) @@ -0,0 +1,75 @@ +/*!\file: transient_3d_core.cpp + * \brief: core of the transient_3d solution + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./solutions.h" +#include "../modules/modules.h" +#include "../include/include.h" +#include "../solvers/solvers.h" + +void transient3d_core(FemModel* femmodel){ + + /*parameters: */ + int verbose=0; + double finaltime; + double dt; + + /*intermediary: */ + int step; + double time; + + //first recover parameters common to all solutions + femmodel->parameters->FindParam(&verbose,VerboseEnum); + femmodel->parameters->FindParam(&finaltime,NdtEnum); + femmodel->parameters->FindParam(&dt,DtEnum); + + /*initialize: */ + step=0; + time=0; + + while(timeelements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,TemperatureEnum,TemperatureAverageEnum); + + if(verbose)_printf_("%s\n","computing new velocity"); + diagnostic_core(femmodel); + + if(verbose)_printf_("%s\n","computing new thickness"); + prognostic_core(femmodel); + + if(verbose)_printf_(" updating geometry\n"); + UpdateGeometryx(femmodel->elements, femmodel->nodes,femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); + + if(verbose)_printf_("%s\n","updating vertices positions"); + UpdateVertexPositionsx(femmodel->elements, femmodel->nodes,femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); + + if(verbose)_printf_("%s\n","saving results:\n"); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MeltingRateEnum,step,time); + + if (step%5==0){ + if(verbose)_printf_("%s"," saving temporary results..."); + OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters); + } + } + +} Index: sm/trunk/src/c/solutions/transient_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/transient_core.cpp (revision 4054) +++ (revision ) @@ -1,30 +1,0 @@ -/*!\file: transient_core.cpp - * \brief: core of the transient solution - */ - -#include "../toolkits/toolkits.h" -#include "../objects/objects.h" -#include "../shared/shared.h" -#include "../EnumDefinitions/EnumDefinitions.h" -#include "./solutions.h" -#include "../modules/modules.h" -#include "../solvers/solvers.h" - -Results* transient_core(Model* model){ - - int dim=-1; - Results* results=NULL; - - //first recover parameters common to all solutions - model->FindParam(&dim,DimEnum); - - //branch out - if(dim==2){ - results=transient_core_2d(model); - } - else if(dim==3){ - results=transient_core_3d(model); - } - else ISSMERROR("%s%i%s"," dimension ",dim," not supported yet!"); - -} Index: sm/trunk/src/c/solutions/transient_core_2d.cpp =================================================================== --- /issm/trunk/src/c/solutions/transient_core_2d.cpp (revision 4054) +++ (revision ) @@ -1,155 +1,0 @@ -/*!\file: transient_2d_core.cpp - * \brief: core of the transient_2d 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* transient_core_2d(Model* model){ - - extern int my_rank; - - /*fem models: */ - FemModel* fem_dh=NULL; - FemModel* fem_dv=NULL; - FemModel* fem_dhu=NULL; - FemModel* fem_ds=NULL; - FemModel* fem_sl=NULL; - FemModel* fem_p=NULL; - - /*output: */ - Results* results=NULL; - Result * result =NULL; - - /*solutions: */ - int step,i; - double time; - - Vec u_g=NULL; - Vec p_g=NULL; - Vec h_g=NULL; - Vec new_h_g=NULL; - Vec h_g_intermediary=NULL; - Vec s_g=NULL; - Vec new_s_g=NULL; - Vec b_g=NULL; - Vec new_b_g=NULL; - double* vx=NULL; - double* vy=NULL; - - /*intermediary: */ - double finaltime; - double dt; - Results* diagnostic_results=NULL; - Results* prognostic_results=NULL; - - - /*flags: */ - int verbose=0; - int numberofnodes; - - /*dof recovery: */ - int dof01[2]={0,1}; - int dof2[1]={2}; - int dof012[3]={0,1,2}; - int dof3[1]={3}; - - //initialize results - results=new Results(); - - /*recover fem models: */ - fem_dh=model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum); - fem_dv=model->GetFormulation(DiagnosticAnalysisEnum,VertAnalysisEnum); - fem_ds=model->GetFormulation(DiagnosticAnalysisEnum,StokesAnalysisEnum); - fem_dhu=model->GetFormulation(DiagnosticAnalysisEnum,HutterAnalysisEnum); - fem_sl=model->GetFormulation(SlopecomputeAnalysisEnum); - fem_p=model->GetFormulation(PrognosticAnalysisEnum); - - - //first recover parameters common to all solutions - model->FindParam(&verbose,VerboseEnum); - model->FindParam(&finaltime,NdtEnum); - model->FindParam(&dt,DtEnum); - model->FindParam(&numberofnodes,NumberOfNodesEnum); - - /*initialize: */ - step=1; - time=0; - - while(timeFindResult(&u_g,"u_g"); - VecFree(&p_g); diagnostic_results->FindResult(&p_g,"p_g"); - delete diagnostic_results; - - //compute new thickness - _printf_("%s\n","computing new thickness"); - - SplitSolutionVectorx(u_g,numberofnodes,2,&vx,&vy); - model->UpdateInputsFromVector(vx,VxEnum,VertexEnum); - model->UpdateInputsFromVector(vy,VyEnum,VertexEnum); - xfree((void**)&vx); xfree((void**)&vy); - - prognostic_results= prognostic_core(model); - - //extract h_g prognostic_results, and erase prognostic_results; - prognostic_results->FindResult(&h_g_intermediary,"h_g"); - delete prognostic_results; - - //update surface and bed using the new thickness - _printf_(" updating geometry\n"); - UpdateGeometryx(&new_h_g,&new_b_g,&new_s_g, fem_p->elements, fem_p->nodes,fem_p->vertices,fem_p->loads, fem_p->materials, fem_p->parameters,h_g_intermediary,b_g,s_g); - - VecFree(&h_g);h_g=new_h_g; - VecFree(&s_g);s_g=new_s_g; - VecFree(&b_g);b_g=new_b_g; - VecFree(&h_g_intermediary); - - //plug into results. - results->AddObject(new Result(results->Size()+1,time,step,"u_g",u_g)); - results->AddObject(new Result(results->Size()+1,time,step,"p_g",p_g)); - results->AddObject(new Result(results->Size()+1,time,step,"h_g",h_g)); - results->AddObject(new Result(results->Size()+1,time,step,"s_g",s_g)); - results->AddObject(new Result(results->Size()+1,time,step,"b_g",b_g)); - - //update inputs - model->UpdateInputsFromVector(h_g,ThicknessEnum,VertexEnum); - model->UpdateInputsFromVector(s_g,SurfaceEnum,VertexEnum); - model->UpdateInputsFromVector(b_g,BedEnum,VertexEnum); - - SplitSolutionVectorx(u_g,numberofnodes,2,&vx,&vy); - model->UpdateInputsFromVector(vx,VxEnum,VertexEnum); - model->UpdateInputsFromVector(vy,VyEnum,VertexEnum); - - } - - /*Add analysis_type to results: */ - results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(TransientAnalysisEnum))); - - /*Free ressources:*/ - VecFree(&u_g); - VecFree(&p_g); - VecFree(&h_g); - VecFree(&s_g); - VecFree(&b_g); - xfree((void**)&vx); - xfree((void**)&vy); - - //return: - return results; - -} Index: sm/trunk/src/c/solutions/transient_core_3d.cpp =================================================================== --- /issm/trunk/src/c/solutions/transient_core_3d.cpp (revision 4054) +++ (revision ) @@ -1,206 +1,0 @@ -/*!\file: transient_3d_core.cpp - * \brief: core of the transient_3d 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* transient_core_3d(Model* model){ - - extern int my_rank; - - /*fem models: */ - FemModel* fem_dh=NULL; - FemModel* fem_dv=NULL; - FemModel* fem_dhu=NULL; - FemModel* fem_ds=NULL; - FemModel* fem_sl=NULL; - FemModel* fem_p=NULL; - FemModel* fem_t=NULL; - FemModel* fem_m=NULL; - - /*output: */ - Results* results=NULL; - - /*solutions: */ - int step; - double time; - - double* vx=NULL; - double* vy=NULL; - double* vz=NULL; - - Vec u_g=NULL; - Vec p_g=NULL; - Vec h_g=NULL; - Vec h_g_intermediary=NULL; - Vec s_g=NULL; - Vec b_g=NULL; - Vec t_g=NULL; - Vec t_g_average=NULL; - double melting_offset; - Vec m_g=NULL; - - /*intermediary: */ - double finaltime; - double dt; - Results* diagnostic_results=NULL; - Results* prognostic_results=NULL; - - - /*flags: */ - int verbose=0; - int numberofnodes; - - /*dof recovery: */ - int dof01[2]={0,1}; - int dof2[1]={2}; - int dof012[3]={0,1,2}; - int dof3[1]={3}; - - //initialize results - results=new Results(); - - /*recover fem models: */ - fem_dh=model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum); - fem_dv=model->GetFormulation(DiagnosticAnalysisEnum,VertAnalysisEnum); - fem_ds=model->GetFormulation(DiagnosticAnalysisEnum,StokesAnalysisEnum); - fem_dhu=model->GetFormulation(DiagnosticAnalysisEnum,HutterAnalysisEnum); - fem_sl=model->GetFormulation(SlopecomputeAnalysisEnum); - fem_p=model->GetFormulation(PrognosticAnalysisEnum); - fem_t=model->GetFormulation(ThermalAnalysisEnum); - fem_m=model->GetFormulation(MeltingAnalysisEnum); - - - //first recover parameters common to all solutions - model->FindParam(&verbose,VerboseEnum); - model->FindParam(&finaltime,NdtEnum); - model->FindParam(&dt,DtEnum); - model->FindParam(&numberofnodes,NumberOfDofsPerNodeEnum); - - /*initialize: */ - step=1; - time=0; - - while(timeUpdateInputsFromVector(t_g,TemperatureEnum,VertexEnum); - model->UpdateInputsFromConstant(melting_offset,MeltingOffsetEnum); - diagnostic_core_linear(&m_g,fem_m,MeltingAnalysisEnum,TransientAnalysisEnum); - - //Compute depth averaged temperature and add to inputs - if(verbose)_printf_("%s\n","computing depth average temperature"); - VecDuplicatePatch(&t_g_average,t_g); - FieldDepthAveragex( t_g_average, fem_t->elements,fem_t->nodes, fem_t->vertices,fem_t->loads, fem_t->materials,fem_t->parameters,"temperature"); - model->UpdateInputsFromVector(t_g_average,TemperatureAverageEnum,VertexEnum); - VecFree(&t_g_average); //not needed anymore - - //Deal with velocities. - diagnostic_results=diagnostic_core(model); - - //extract u_g and p_g from diagnostic_results, and erase diagnostic_results; - diagnostic_results->FindResult(&u_g,"u_g"); - diagnostic_results->FindResult(&p_g,"p_g"); - delete diagnostic_results; - - //compute new thickness - if(verbose)_printf_("%s\n","computing new thickness"); - - SplitSolutionVectorx(u_g,numberofnodes,3,&vx,&vy,&vz); - model->UpdateInputsFromVector(vx,VxEnum,VertexEnum); - model->UpdateInputsFromVector(vy,VyEnum,VertexEnum); - model->UpdateInputsFromVector(vz,VzEnum,VertexEnum); - xfree((void**)&vx); xfree((void**)&vy); xfree((void**)&vz); - - prognostic_results=prognostic_core(model); - - //extract h_g prognostic_results, and erase prognostic_results; - prognostic_results->FindResult(&h_g_intermediary,"h_g"); - delete prognostic_results; - - //update surface and bed using the new thickness - if(verbose)_printf_(" updating geometry\n"); - UpdateGeometryx(&h_g,&b_g,&s_g, fem_p->elements, fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials, fem_p->parameters,h_g_intermediary,b_g,s_g); - VecFree(&h_g_intermediary); - - if(verbose)_printf_("%s\n","updating node positions"); - UpdateVertexPositionsx( fem_dh ->vertices,h_g,b_g); - UpdateVertexPositionsx( fem_dv ->vertices,h_g,b_g); - UpdateVertexPositionsx( fem_dhu->vertices,h_g,b_g); - UpdateVertexPositionsx( fem_ds ->vertices,h_g,b_g); - UpdateVertexPositionsx( fem_sl ->vertices,h_g,b_g); - UpdateVertexPositionsx( fem_p ->vertices,h_g,b_g); - UpdateVertexPositionsx( fem_t ->vertices,h_g,b_g); - UpdateVertexPositionsx( fem_m ->vertices,h_g,b_g); - - //plug into results. - results->AddObject(new Result(results->Size()+1,time,step,"u_g",u_g)); - results->AddObject(new Result(results->Size()+1,time,step,"p_g",p_g)); - results->AddObject(new Result(results->Size()+1,time,step,"h_g",h_g)); - results->AddObject(new Result(results->Size()+1,time,step,"s_g",s_g)); - results->AddObject(new Result(results->Size()+1,time,step,"b_g",b_g)); - results->AddObject(new Result(results->Size()+1,time,step,"t_g",t_g)); - results->AddObject(new Result(results->Size()+1,time,step,"m_g",m_g)); - - //update inputs - model->UpdateInputsFromVector(h_g,ThicknessEnum,VertexEnum); - model->UpdateInputsFromVector(s_g,SurfaceEnum,VertexEnum); - model->UpdateInputsFromVector(b_g,BedEnum,VertexEnum); - - SplitSolutionVectorx(u_g,numberofnodes,3,&vx,&vy,&vz); - model->UpdateInputsFromVector(vx,VxEnum,VertexEnum); - model->UpdateInputsFromVector(vy,VyEnum,VertexEnum); - model->UpdateInputsFromVector(vz,VzEnum,VertexEnum); - xfree((void**)&vx); xfree((void**)&vy); xfree((void**)&vz); - - model->UpdateInputsFromVector(p_g,PressureEnum,VertexEnum); - model->UpdateInputsFromVector(t_g,TemperatureEnum,VertexEnum); - - //Save temporary results every 5 steps - if (step%5==0){ - _printf_("%s"," saving temporary results..."); - - char* outputfilename=NULL; - - model->FindParam(&outputfilename,OutputFileNameEnum); - OutputResults(results,outputfilename); - - xfree((void**)&outputfilename); - - _printf_("%s\n"," done."); - } - } - - /*Add analysis_type to results: */ - results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(TransientAnalysisEnum))); - - - /*Free ressources:*/ - VecFree(&u_g); - VecFree(&p_g); - VecFree(&h_g); - VecFree(&s_g); - VecFree(&b_g); - VecFree(&t_g); - VecFree(&m_g); - xfree((void**)&vx); - xfree((void**)&vy); - xfree((void**)&vz); - - //return: - return results; - -} Index: /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp (revision 4054) +++ /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp (revision 4055) @@ -9,5 +9,5 @@ #include "./solutions.h" -void solver_diagnostic_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, FemModel* fem,bool conserve_loads, int analysis_type,int sub_analysis_type){ +void solver_diagnostic_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, FemModel* fem,bool conserve_loads){ @@ -61,5 +61,5 @@ /*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); + GetSolutionFromInputsx(&ug, fem->elements, fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters); Reducevectorgtofx(&uf, ug, fem->nodesets); @@ -72,9 +72,9 @@ 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); + SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag); 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); + PenaltySystemMatricesx(Kgg, pg,NULL,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag); if (verbose) _printf_(" reducing matrix from g to f set\n"); @@ -102,9 +102,9 @@ //Update inputs using new solution: - UpdateInputsFromSolutionx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,ug,analysis_type, sub_analysis_type); + UpdateInputsFromSolutionx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,ug); //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); + PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters); if(verbose)_printf_(" number of unstable constraints: %i\n",num_unstable_constraints); @@ -146,5 +146,5 @@ 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); + SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag); Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets); MatFree(&Kgg);VecFree(&pg); Index: /issm/trunk/src/c/solvers/solver_linear.cpp =================================================================== --- /issm/trunk/src/c/solvers/solver_linear.cpp (revision 4054) +++ /issm/trunk/src/c/solvers/solver_linear.cpp (revision 4055) @@ -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){ /*parameters:*/ @@ -33,9 +33,9 @@ //*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); + SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag); 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); + PenaltySystemMatricesx(Kgg, pg,NULL,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag); /*!Reduce matrix from g to f size:*/ @@ -56,5 +56,5 @@ //Update inputs using new solution: - UpdateInputsFromSolutionx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,ug,analysis_type, sub_analysis_type); + UpdateInputsFromSolutionx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,ug); /*free ressources: */ Index: /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp (revision 4054) +++ /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp (revision 4055) @@ -8,5 +8,5 @@ #include "../modules/modules.h" -void solver_thermal_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem,int analysis_type,int sub_analysis_type){ +void solver_thermal_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem){ /*solution : */ @@ -63,5 +63,5 @@ /*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); + SystemMatricesx(&Kgg_nopenalty, &pg_nopenalty,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag); } @@ -71,10 +71,10 @@ //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); + PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag); } else{ - SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); + SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag); //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); + PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag); } @@ -108,9 +108,9 @@ //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); + PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters); //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); + UpdateInputsFromSolutionx(fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,tg); if (!converged){ Index: /issm/trunk/src/c/solvers/solvers.h =================================================================== --- /issm/trunk/src/c/solvers/solvers.h (revision 4054) +++ /issm/trunk/src/c/solvers/solvers.h (revision 4055) @@ -12,8 +12,8 @@ class FemModel; -void solver_thermal_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem,int analysis_type,int sub_analysis_type); -void solver_diagnostic_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* fem,bool conserve_loads,int analysis_type,int sub_analysis_type); -void solver_linear(Vec* pug, FemModel* femmodel,int analysis_type,int sub_analysis_type); -void solver_adjoint(Vec* pug, FemModel* femmodel,int analysis_type,int sub_analysis_type); +void solver_thermal_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* femmodel); +void solver_diagnostic_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* femmodel,bool conserve_loads); +void solver_linear(Vec* pug, FemModel* femmodel); +void solver_adjoint(Vec* pug, FemModel* femmodel);