Index: /issm/trunk/src/c/DataSet/DataSet.cpp =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.cpp (revision 464) +++ /issm/trunk/src/c/DataSet/DataSet.cpp (revision 465) @@ -900,5 +900,5 @@ #undef __FUNCT__ #define __FUNCT__ "DataSet::CreateKMatrix" -void DataSet::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){ +void DataSet::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ vector::iterator object; @@ -911,10 +911,10 @@ element=(Element*)(*object); - element->CreateKMatrix(Kgg,inputs,analysis_type); + element->CreateKMatrix(Kgg,inputs,analysis_type,sub_analysis_type); } if(EnumIsLoad((*object)->Enum())){ load=(Load*)(*object); - load->CreateKMatrix(Kgg,inputs,analysis_type); + load->CreateKMatrix(Kgg,inputs,analysis_type,sub_analysis_type); } } @@ -924,5 +924,5 @@ #undef __FUNCT__ #define __FUNCT__ "DataSet::CreatePVector" -void DataSet::CreatePVector(Vec pg,void* inputs,int analysis_type){ +void DataSet::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){ vector::iterator object; @@ -935,10 +935,10 @@ element=(Element*)(*object); - element->CreatePVector(pg,inputs,analysis_type); + element->CreatePVector(pg,inputs,analysis_type,sub_analysis_type); } if(EnumIsLoad((*object)->Enum())){ load=(Load*)(*object); - load->CreatePVector(pg,inputs,analysis_type); + load->CreatePVector(pg,inputs,analysis_type,sub_analysis_type); } } @@ -984,5 +984,5 @@ #undef __FUNCT__ #define __FUNCT__ "DataSet::PenaltyCreateKMatrix" -void DataSet::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type){ +void DataSet::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){ vector::iterator object; @@ -994,5 +994,5 @@ load=(Load*)(*object); - load->PenaltyCreateKMatrix(Kgg,inputs,kmax,analysis_type); + load->PenaltyCreateKMatrix(Kgg,inputs,kmax,analysis_type,sub_analysis_type); } } @@ -1002,5 +1002,5 @@ #undef __FUNCT__ #define __FUNCT__ "DataSet::PenaltyCreatePVector" -void DataSet::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type){ +void DataSet::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){ vector::iterator object; @@ -1012,5 +1012,5 @@ load=(Load*)(*object); - load->PenaltyCreatePVector(pg,inputs,kmax,analysis_type); + load->PenaltyCreatePVector(pg,inputs,kmax,analysis_type,sub_analysis_type); } } @@ -1028,5 +1028,5 @@ #undef __FUNCT__ #define __FUNCT__ "RiftConstraints" -void DataSet::RiftConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type){ +void DataSet::RiftConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type){ throw ErrorException(__FUNCT__," not implemented yet!"); @@ -1045,5 +1045,5 @@ #undef __FUNCT__ #define __FUNCT__ "MeltingConstraints" -void DataSet::MeltingConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type){ +void DataSet::MeltingConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type){ throw ErrorException(__FUNCT__," not implemented yet!"); @@ -1084,5 +1084,5 @@ -void DataSet::Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type){ +void DataSet::Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type){ @@ -1095,12 +1095,12 @@ element=(Element*)(*object); - element->Du(du_g,u_g,u_g_obs,inputs,analysis_type); - } - } - - -} - -void DataSet::Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,char* control_type){ + element->Du(du_g,u_g,u_g_obs,inputs,analysis_type,sub_analysis_type); + } + } + + +} + +void DataSet::Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type){ @@ -1113,5 +1113,5 @@ element=(Element*)(*object); - element->Gradj(grad_g,u_g,lambda_g,inputs,analysis_type,control_type); + element->Gradj(grad_g,u_g,lambda_g,inputs,analysis_type,sub_analysis_type,control_type); } } @@ -1120,5 +1120,5 @@ } -void DataSet::Misfit(double* pJ, double* u_g,double* u_g_obs,void* inputs,int analysis_type){ +void DataSet::Misfit(double* pJ, double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type){ double J=0;; @@ -1132,5 +1132,5 @@ element=(Element*)(*object); - J+=element->Misfit(u_g,u_g_obs,inputs,analysis_type); + J+=element->Misfit(u_g,u_g_obs,inputs,analysis_type,sub_analysis_type); } Index: /issm/trunk/src/c/DataSet/DataSet.h =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.h (revision 464) +++ /issm/trunk/src/c/DataSet/DataSet.h (revision 465) @@ -47,5 +47,5 @@ void DistributeDofs(int numberofnodes,int numdofspernode); void CreatePartitioningVector(Vec* ppartition,int numnods,int numdofspernode); - void DistributeNumDofs(int** pnumdofspernode,int numberofnodes,int analysis_type); + void DistributeNumDofs(int** pnumdofspernode,int numberofnodes,int analysis_type,int sub_analysis_type); void FlagClones(int numberofnodes); int NumberOfDofs(); @@ -61,17 +61,17 @@ void SetSorting(int* in_sorted_ids,int* in_id_offsets); void Sort(); - void CreateKMatrix(Mat Kgg,void* inputs, int analysis_type); - void CreatePVector(Vec pg,void* inputs, int analysis_type); + void CreateKMatrix(Mat Kgg,void* inputs, int analysis_type,int sub_analysis_type); + void CreatePVector(Vec pg,void* inputs, int analysis_type,int sub_analysis_type); void UpdateFromInputs(void* inputs); - void PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type); - void PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type); + void PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type); + void PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type); int RiftIsPresent(); - void RiftConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type); + void RiftConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type); int MeltingIsPresent(); - void MeltingConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type); + void MeltingConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type); DataSet* Copy(void); - void Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type); - void Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,char* control_type); - void Misfit(double* pJ, double* u_g,double* u_g_obs,void* inputs,int analysis_type); + void Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type); + void Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type); + void Misfit(double* pJ, double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type); void VelocityExtrude(Vec ug,double* ug_serial); void SlopeExtrude(Vec sg,double* sg_serial); Index: /issm/trunk/src/c/Dux/Dux.cpp =================================================================== --- /issm/trunk/src/c/Dux/Dux.cpp (revision 464) +++ /issm/trunk/src/c/Dux/Dux.cpp (revision 465) @@ -14,5 +14,5 @@ void Dux( Vec* pdu_g, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials, - double* u_g, double* u_g_obs, ParameterInputs* inputs,int analysis_type){ + double* u_g, double* u_g_obs, ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ int i; @@ -33,5 +33,5 @@ /*Compute velocity difference: */ - elements->Du(du_g, u_g,u_g_obs,inputs,analysis_type); + elements->Du(du_g, u_g,u_g_obs,inputs,analysis_type,sub_analysis_type); /*Assemble vector: */ Index: /issm/trunk/src/c/Dux/Dux.h =================================================================== --- /issm/trunk/src/c/Dux/Dux.h (revision 464) +++ /issm/trunk/src/c/Dux/Dux.h (revision 465) @@ -10,5 +10,5 @@ /* local prototypes: */ void Dux( Vec* pdu_g, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials, - double* u_g, double* u_g_obs, ParameterInputs* inputs,int analysis_type); + double* u_g, double* u_g_obs, ParameterInputs* inputs,int analysis_type,int sub_analysis_type); #endif /* _DUX_H */ Index: /issm/trunk/src/c/EnumDefinitions/AnalysisTypeAsEnum.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/AnalysisTypeAsEnum.cpp (revision 464) +++ /issm/trunk/src/c/EnumDefinitions/AnalysisTypeAsEnum.cpp (revision 465) @@ -15,47 +15,50 @@ } - if (strcmp(analysis_type,"control")==0){ + if (strcmp(analysis_type,"diagnostic")==0){ + return DiagnosticAnalysisEnum(); + } + else if (strcmp(analysis_type,"control")==0){ return ControlAnalysisEnum(); } - else if (strcmp(analysis_type,"diagnostic_horiz")==0){ - return DiagnosticHorizAnalysisEnum(); - } - else if (strcmp(analysis_type,"diagnostic_vert")==0){ - return DiagnosticVertAnalysisEnum(); + else if (strcmp(analysis_type,"thermal")==0){ + return ThermalAnalysisEnum(); } else if (strcmp(analysis_type,"prognostic")==0){ return PrognosticAnalysisEnum(); } - else if (strcmp(analysis_type,"thermalsteady")==0){ - return ThermalSteadyAnalysisEnum(); - } - else if (strcmp(analysis_type,"thermaltransient")==0){ - return ThermalTransientAnalysisEnum(); - } else if (strcmp(analysis_type,"melting")==0){ return MeltingAnalysisEnum(); - } - else if (strcmp(analysis_type,"diagnostic_stokes")==0){ - return DiagnosticStokesAnalysisEnum(); - } - else if (strcmp(analysis_type,"diagnostic_hutter")==0){ - return DiagnosticHutterAnalysisEnum(); } else if (strcmp(analysis_type,"slope_compute")==0){ return SlopeComputeAnalysisEnum(); } - else if (strcmp(analysis_type,"surface_slope_compute_x")==0){ - return SurfaceSlopeComputeXAnalysisEnum(); + else if (strcmp(analysis_type,"stokes")==0){ + return StokesAnalysisEnum(); } - else if (strcmp(analysis_type,"surface_slope_compute_y")==0){ - return SurfaceSlopeComputeYAnalysisEnum(); + else if (strcmp(analysis_type,"hutter")==0){ + return HutterAnalysisEnum(); } - else if (strcmp(analysis_type,"bed_slope_compute_x")==0){ - return BedSlopeComputeXAnalysisEnum(); + else if (strcmp(analysis_type,"surfacex")==0){ + return SurfaceXAnalysisEnum(); } - else if (strcmp(analysis_type,"bed_slope_compute_y")==0){ - return BedSlopeComputeYAnalysisEnum(); + else if (strcmp(analysis_type,"surfacey")==0){ + return SurfaceYAnalysisEnum(); } - else throw ErrorException(__FUNCT__," could not recognized analysis type!"); + else if (strcmp(analysis_type,"bedx")==0){ + return BedXAnalysisEnum(); + } + else if (strcmp(analysis_type,"bedy")==0){ + return BedYAnalysisEnum(); + } + else if (strcmp(analysis_type,"horiz")==0){ + return HorizAnalysisEnum(); + } + else if (strcmp(analysis_type,"vert")==0){ + return VertAnalysisEnum(); + } + else if (strcmp(analysis_type,"")==0){ + return NoneAnalysisEnum(); + } + else throw ErrorException(__FUNCT__,exprintf("%s%s"," could not recognize analysis type: ",analysis_type)); } Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp (revision 464) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp (revision 465) @@ -15,20 +15,19 @@ /*analysis types: */ -int StructuralAnalysisEnum(void){ return 20; } -int ThermalAnalysisEnum(void){ return 21; } -int ControlAnalysisEnum(void){ return 22; } -int DiagnosticHorizAnalysisEnum(void){ return 23; } -int DiagnosticVertAnalysisEnum(void){ return 24; } -int PrognosticAnalysisEnum(void){ return 26; } -int ThermalSteadyAnalysisEnum(void) { return 27; } -int ThermalTransientAnalysisEnum(void){ return 28; } -int MeltingAnalysisEnum(void){ return 29; } -int DiagnosticStokesAnalysisEnum(void){ return 30; } -int DiagnosticHutterAnalysisEnum(void){ return 31; } -int SlopeComputeAnalysisEnum(void){ return 32; } -int SurfaceSlopeComputeXAnalysisEnum(void){ return 33; } -int SurfaceSlopeComputeYAnalysisEnum(void){ return 34; } -int BedSlopeComputeXAnalysisEnum(void){ return 35; } -int BedSlopeComputeYAnalysisEnum(void){ return 36; } +int DiagnosticAnalysisEnum(void){ return 20; } +int ControlAnalysisEnum(void){ return 21; } +int ThermalAnalysisEnum(void){ return 22; } +int PrognosticAnalysisEnum(void){ return 23; } +int MeltingAnalysisEnum(void){ return 24; } +int SlopeComputeAnalysisEnum(void){ return 25; } +int StokesAnalysisEnum(void){ return 26; } +int HutterAnalysisEnum(void){ return 27; } +int SurfaceXAnalysisEnum(void){ return 28; } +int SurfaceYAnalysisEnum(void){ return 29; } +int BedXAnalysisEnum(void){ return 30; } +int BedYAnalysisEnum(void){ return 31; } +int HorizAnalysisEnum(void){ return 32; } +int VertAnalysisEnum(void){ return 33; } +int NoneAnalysisEnum(void){ return 34; } /*datasets: */ Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 464) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 465) @@ -34,18 +34,19 @@ /* analysis enums: */ +int DiagnosticAnalysisEnum(void); int ControlAnalysisEnum(void); -int DiagnosticHorizAnalysisEnum(void); -int DiagnosticVertAnalysisEnum(void); +int ThermalAnalysisEnum(void); int PrognosticAnalysisEnum(void); -int ThermalSteadyAnalysisEnum(void); -int ThermalTransientAnalysisEnum(void); int MeltingAnalysisEnum(void); -int DiagnosticStokesAnalysisEnum(void); -int DiagnosticHutterAnalysisEnum(void); int SlopeComputeAnalysisEnum(void); -int SurfaceSlopeComputeXAnalysisEnum(void); -int SurfaceSlopeComputeYAnalysisEnum(void); -int BedSlopeComputeXAnalysisEnum(void); -int BedSlopeComputeYAnalysisEnum(void); +int StokesAnalysisEnum(void); +int HutterAnalysisEnum(void); +int SurfaceXAnalysisEnum(void); +int SurfaceYAnalysisEnum(void); +int BedXAnalysisEnum(void); +int BedYAnalysisEnum(void); +int HorizAnalysisEnum(void); +int VertAnalysisEnum(void); +int NoneAnalysisEnum(void); Index: /issm/trunk/src/c/Gradjx/Gradjx.cpp =================================================================== --- /issm/trunk/src/c/Gradjx/Gradjx.cpp (revision 464) +++ /issm/trunk/src/c/Gradjx/Gradjx.cpp (revision 465) @@ -14,5 +14,5 @@ void Gradjx( Vec* pgrad_g, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials, - double* u_g, double* lambda_g, ParameterInputs* inputs,int analysis_type,char* control_type){ + double* u_g, double* lambda_g, ParameterInputs* inputs,int analysis_type,int sub_analysis_type,char* control_type){ int i; @@ -33,5 +33,5 @@ /*Compute gradients: */ - elements->Gradj(grad_g, u_g,lambda_g,inputs,analysis_type,control_type); + elements->Gradj(grad_g, u_g,lambda_g,inputs,analysis_type,sub_analysis_type,control_type); /*Assemble vector: */ Index: /issm/trunk/src/c/Gradjx/Gradjx.h =================================================================== --- /issm/trunk/src/c/Gradjx/Gradjx.h (revision 464) +++ /issm/trunk/src/c/Gradjx/Gradjx.h (revision 465) @@ -10,5 +10,5 @@ /* local prototypes: */ void Gradjx( Vec* pgrad_g, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials, - double* u_g, double* lambda_g, ParameterInputs* inputs,int analysis_type,char* control_type); + double* u_g, double* lambda_g, ParameterInputs* inputs,int analysis_type,int sub_analysis_type,char* control_type); #endif /* _GRADJX_H */ Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 464) +++ /issm/trunk/src/c/Makefile.am (revision 465) @@ -79,4 +79,5 @@ ./shared/Dofs/dofs.h\ ./shared/Dofs/dofsetgen.cpp\ + ./shared/Dofs/DistributeNumDofs.cpp\ ./shared/Numerics/numerics.h\ ./shared/Numerics/GaussPoints.h\ @@ -86,251 +87,4 @@ ./shared/Numerics/BrentSearch.cpp\ ./shared/Numerics/OptFunc.cpp\ - ./shared/Numerics/extrema.cpp\ - ./shared/Numerics/DistributeNumDofs.cpp\ - ./shared/Exceptions/exceptions.h\ - ./shared/Exceptions/Exceptions.cpp\ - ./shared/Exceptions/exprintf.cpp\ - ./shared/Exp/exp.h\ - ./shared/Exp/IsInPoly.cpp\ - ./shared/Exp/IsInPolySerial.cpp\ - ./shared/Exp/DomainOutlineRead.cpp\ - ./shared/TriMesh/trimesh.h\ - ./shared/TriMesh/AssociateSegmentToElement.cpp\ - ./shared/TriMesh/GridInsideHole.cpp\ - ./shared/TriMesh/OrderSegments.cpp\ - ./shared/TriMesh/SplitMeshForRifts.cpp\ - ./shared/TriMesh/TriMeshUtils.cpp\ - ./shared/Sorting/binary_search.cpp\ - ./shared/Sorting/sorting.h\ - ./shared/Elements/elements.h\ - ./shared/Elements/ResolvePointers.cpp\ - ./shared/Elements/Paterson.cpp\ - ./shared/Elements/GetElementNodeData.cpp\ - ./toolkits/petsc\ - ./toolkits/petsc/patches\ - ./toolkits/petsc/patches/petscpatches.h\ - ./toolkits/petsc/patches/MatlabMatrixToPetscMatrix.cpp\ - ./toolkits/petsc/patches/MatlabVectorToPetscVector.cpp\ - ./toolkits/petsc/patches/PetscMatrixToMatlabMatrix.cpp\ - ./toolkits/petsc/patches/PetscVectorToMatlabVector.cpp\ - ./toolkits/petsc/patches/MatlabMatrixToDoubleMatrix.cpp\ - ./toolkits/petsc/patches/MatlabVectorToDoubleVector.cpp\ - ./toolkits/petsc/patches/PetscDetermineLocalSize.cpp\ - ./toolkits/petsc/patches/VecTranspose.cpp\ - ./toolkits/petsc/patches/VecToMPISerial.cpp\ - ./toolkits/petsc/patches/MatToSerial.cpp\ - ./toolkits/petsc/patches/VecMerge.cpp\ - ./toolkits/petsc/patches/NewVec.cpp\ - ./toolkits/petsc/patches/NewVecFromLocalSize.cpp\ - ./toolkits/petsc/patches/PetscOptionsDetermineSolverType.cpp\ - ./toolkits/petsc/patches/PetscOptionsInsertMultipleString.cpp\ - ./toolkits/petsc/patches/NewMat.cpp\ - ./toolkits/petsc/patches/VecFree.cpp\ - ./toolkits/petsc/patches/VecDuplicatePatch.cpp\ - ./toolkits/petsc/patches/KSPFree.cpp\ - ./toolkits/petsc/patches/ISFree.cpp\ - ./toolkits/petsc/patches/MatFree.cpp\ - ./toolkits/petsc/patches/GetOwnershipBoundariesFromRange.cpp\ - ./toolkits/petsc/patches/VecPartition.cpp\ - ./toolkits/petsc/patches/MatPartition.cpp\ - ./toolkits/petsc/patches/MatInvert.cpp\ - ./toolkits/petsc/patches/MatMultPatch.cpp\ - ./toolkits/petsc/petscincludes.h\ - ./toolkits/mpi/mpiincludes.h\ - ./toolkits/mpi/patches/mpipatches.h\ - ./toolkits/mpi/patches/MPI_Upperrow.cpp\ - ./toolkits/mpi/patches/MPI_Lowerrow.cpp\ - ./toolkits/mpi/patches/MPI_Boundariesfromrange.cpp\ - ./toolkits/metis/metisincludes.h\ - ./toolkits/triangle/triangleincludes.h\ - ./toolkits.h\ - ./io/io.h\ - ./io/FetchData.cpp\ - ./io/WriteData.cpp\ - ./io/WriteDataToDisk.cpp\ - ./io/SerialFetchData.cpp\ - ./io/SerialWriteData.cpp\ - ./io/ParallelFetchData.cpp\ - ./io/ParallelFetchInteger.cpp\ - ./io/ParallelFetchMat.cpp\ - ./io/ParallelFetchScalar.cpp\ - ./io/ParallelFetchString.cpp\ - ./io/ModelFetchData.cpp\ - ./io/WriteNodeSets.cpp\ - ./io/WriteParams.cpp\ - ./io/FetchNodeSets.cpp\ - ./io/FetchRifts.cpp\ - ./io/ParameterInputsInit.cpp\ - ./EnumDefinitions/EnumDefinitions.h\ - ./EnumDefinitions/EnumDefinitions.cpp\ - ./EnumDefinitions/AnalysisTypeAsEnum.cpp\ - ./ModelProcessorx/Model.h\ - ./ModelProcessorx/Model.cpp\ - ./ModelProcessorx/CreateElementsNodesAndMaterials.cpp\ - ./ModelProcessorx/CreateConstraints.cpp \ - ./ModelProcessorx/CreateLoads.cpp\ - ./ModelProcessorx/CreateParameters.cpp\ - ./ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp\ - ./ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp \ - ./ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp\ - ./ModelProcessorx/DiagnosticHoriz/CreateParametersDiagnosticHoriz.cpp\ - ./ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp\ - ./ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp \ - ./ModelProcessorx/DiagnosticVert/CreateLoadsDiagnosticVert.cpp\ - ./ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp\ - ./ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp \ - ./ModelProcessorx/DiagnosticHutter/CreateLoadsDiagnosticHutter.cpp\ - ./ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp\ - ./ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp \ - ./ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp\ - ./ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp\ - ./ModelProcessorx/SlopeCompute/CreateConstraintsSlopeCompute.cpp \ - ./ModelProcessorx/SlopeCompute/CreateLoadsSlopeCompute.cpp\ - ./ModelProcessorx/Control/CreateParametersControl.cpp\ - ./Dofx/Dofx.h\ - ./Dofx/Dofx.cpp\ - ./Dux/Dux.h\ - ./Dux/Dux.cpp\ - ./GriddataMeshToGridx/GriddataMeshToGridx.h\ - ./GriddataMeshToGridx/GriddataMeshToGridx.cpp\ - ./ControlConstrainx/ControlConstrainx.h\ - ./ControlConstrainx/ControlConstrainx.cpp\ - ./Misfitx/Misfitx.h\ - ./Misfitx/Misfitx.cpp\ - ./Orthx/Orthx.h\ - ./Orthx/Orthx.cpp\ - ./Gradjx/Gradjx.h\ - ./Gradjx/Gradjx.cpp\ - ./UpdateFromInputsx/UpdateFromInputsx.h\ - ./UpdateFromInputsx/UpdateFromInputsx.cpp\ - ./ConfigureObjectsx/ConfigureObjectsx.h\ - ./ConfigureObjectsx/ConfigureObjectsx.cpp\ - ./ComputePressurex/ComputePressurex.h\ - ./ComputePressurex/ComputePressurex.cpp\ - ./BuildNodeSetsx/BuildNodeSetsx.h\ - ./BuildNodeSetsx/BuildNodeSetsx.cpp\ - ./BuildNodeSetsx/PartitionSets.cpp\ - ./SpcNodesx/SpcNodesx.h\ - ./SpcNodesx/SpcNodesx.cpp\ - ./MpcNodesx/MpcNodesx.h\ - ./MpcNodesx/MpcNodesx.cpp\ - ./DataInterpx/DataInterpx.cpp\ - ./DataInterpx/DataInterpx.h\ - ./HoleFillerx/HoleFillerx.cpp\ - ./HoleFillerx/HoleFillerx.h\ - ./MeshPartitionx/MeshPartitionx.cpp\ - ./MeshPartitionx/MeshPartitionx.h\ - ./ContourToMeshx/ContourToMeshx.cpp\ - ./ContourToMeshx/ContourToMeshx.h\ - ./Reducevectorgtosx/Reducevectorgtosx.cpp\ - ./Reducevectorgtosx/Reducevectorgtosx.h\ - ./Reducematrixfromgtofx/Reducematrixfromgtofx.cpp\ - ./Reducematrixfromgtofx/Reducematrixfromgton.cpp\ - ./Reducematrixfromgtofx/Reducematrixfromgtofx.h\ - ./Reduceloadfromgtofx/Reduceloadfromgtofx.h\ - ./Reduceloadfromgtofx/Reduceloadfromgtofx.cpp\ - ./NormalizeConstraintsx/NormalizeConstraintsx.cpp\ - ./NormalizeConstraintsx/NormalizeConstraintsx.h\ - ./SystemMatricesx/SystemMatricesx.cpp\ - ./SystemMatricesx/SystemMatricesx.h\ - ./PenaltyConstraintsx/PenaltyConstraintsx.cpp\ - ./PenaltyConstraintsx/PenaltyConstraintsx.h\ - ./PenaltySystemMatricesx/PenaltySystemMatricesx.cpp\ - ./PenaltySystemMatricesx/PenaltySystemMatricesx.h\ - ./Solverx/Solverx.cpp\ - ./Solverx/Solverx.h\ - ./Mergesolutionfromftogx/Mergesolutionfromftogx.cpp\ - ./Mergesolutionfromftogx/Mergesolutionfromftogx.h\ - ./ProcessParamsx/ProcessParamsx.cpp\ - ./ProcessParamsx/ProcessParamsx.h\ - ./VelocityExtrudex/VelocityExtrudex.cpp\ - ./VelocityExtrudex/VelocityExtrudex.h\ - ./SlopeExtrudex/SlopeExtrudex.cpp\ - ./SlopeExtrudex/SlopeExtrudex.h - - - - -libISSM_a_CXXFLAGS = -fPIC -DMATLAB -D_SERIAL_ -ansi -D_GNU_SOURCE -fno-omit-frame-pointer -pthread -D_CPP_ -if LARGEARRAYS -libISSM_a_CXXFLAGS += -D__GCC4BUILD__ -else -libISSM_a_CXXFLAGS += -DMX_COMPAT_32 -endif - - - - -#Parallel compilation -libpISSM_a_SOURCES = ./objects/objects.h\ - ./objects/Object.h\ - ./objects/Element.h\ - ./objects/Element.cpp\ - ./objects/Material.h\ - ./objects/Material.cpp\ - ./objects/Load.h\ - ./objects/Load.cpp\ - ./objects/SolverEnum.h\ - ./objects/Contour.h\ - ./objects/Contour.cpp\ - ./objects/OptArgs.h\ - ./objects/OptPars.h\ - ./objects/Friction.h\ - ./objects/Friction.cpp\ - ./objects/Node.h\ - ./objects/Node.cpp\ - ./objects/Tria.h\ - ./objects/Tria.cpp\ - ./objects/Sing.h\ - ./objects/Sing.cpp\ - ./objects/Beam.h\ - ./objects/Beam.cpp\ - ./objects/Penta.h\ - ./objects/Penta.cpp\ - ./objects/Matice.h\ - ./objects/Matice.cpp\ - ./objects/Matpar.h\ - ./objects/Matpar.cpp\ - ./objects/Input.h\ - ./objects/Input.cpp\ - ./objects/ParameterInputs.h\ - ./objects/ParameterInputs.cpp\ - ./objects/Spc.cpp\ - ./objects/Spc.h\ - ./objects/Rgb.cpp\ - ./objects/Rgb.h\ - ./objects/Penpair.cpp\ - ./objects/Penpair.h\ - ./objects/Pengrid.cpp\ - ./objects/Pengrid.h\ - ./objects/Icefront.cpp\ - ./objects/Icefront.h\ - ./objects/Param.cpp\ - ./objects/Param.h\ - ./objects/NodeSets.cpp\ - ./objects/NodeSets.h\ - ./DataSet/DataSet.cpp\ - ./DataSet/DataSet.h\ - ./shared/shared.h\ - ./shared/Alloc/alloc.h\ - ./shared/Alloc/alloc.cpp\ - ./shared/Matlab/matlabshared.h\ - ./shared/Matlab/PrintfFunction.cpp\ - ./shared/Matlab/ModuleBoot.cpp\ - ./shared/Matlab/mxGetAssignedField.cpp\ - ./shared/Matlab/mxGetField.cpp\ - ./shared/Matlab/CheckNumMatlabArguments.cpp\ - ./shared/Matrix/matrix.h\ - ./shared/Matrix/MatrixUtils.cpp\ - ./shared/Dofs/dofs.h\ - ./shared/Dofs/dofsetgen.cpp\ - ./shared/Numerics/numerics.h\ - ./shared/Numerics/GaussPoints.h\ - ./shared/Numerics/GaussPoints.cpp\ - ./shared/Numerics/cross.cpp\ - ./shared/Numerics/norm.cpp\ - ./shared/Numerics/BrentSearch.cpp\ - ./shared/Numerics/OptFunc.cpp\ - ./shared/Numerics/DistributeNumDofs.cpp\ ./shared/Numerics/extrema.cpp\ ./shared/Exceptions/exceptions.h\ @@ -413,7 +167,249 @@ ./ModelProcessorx/Model.h\ ./ModelProcessorx/Model.cpp\ - ./ModelProcessorx/CreateElementsNodesAndMaterials.cpp\ - ./ModelProcessorx/CreateConstraints.cpp \ - ./ModelProcessorx/CreateLoads.cpp\ + ./ModelProcessorx/CreateDataSets.cpp\ + ./ModelProcessorx/CreateParameters.cpp\ + ./ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp\ + ./ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp \ + ./ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp\ + ./ModelProcessorx/DiagnosticHoriz/CreateParametersDiagnosticHoriz.cpp\ + ./ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp\ + ./ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp \ + ./ModelProcessorx/DiagnosticVert/CreateLoadsDiagnosticVert.cpp\ + ./ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp\ + ./ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp \ + ./ModelProcessorx/DiagnosticHutter/CreateLoadsDiagnosticHutter.cpp\ + ./ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp\ + ./ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp \ + ./ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp\ + ./ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp\ + ./ModelProcessorx/SlopeCompute/CreateConstraintsSlopeCompute.cpp \ + ./ModelProcessorx/SlopeCompute/CreateLoadsSlopeCompute.cpp\ + ./ModelProcessorx/Control/CreateParametersControl.cpp\ + ./Dofx/Dofx.h\ + ./Dofx/Dofx.cpp\ + ./Dux/Dux.h\ + ./Dux/Dux.cpp\ + ./GriddataMeshToGridx/GriddataMeshToGridx.h\ + ./GriddataMeshToGridx/GriddataMeshToGridx.cpp\ + ./ControlConstrainx/ControlConstrainx.h\ + ./ControlConstrainx/ControlConstrainx.cpp\ + ./Misfitx/Misfitx.h\ + ./Misfitx/Misfitx.cpp\ + ./Orthx/Orthx.h\ + ./Orthx/Orthx.cpp\ + ./Gradjx/Gradjx.h\ + ./Gradjx/Gradjx.cpp\ + ./UpdateFromInputsx/UpdateFromInputsx.h\ + ./UpdateFromInputsx/UpdateFromInputsx.cpp\ + ./ConfigureObjectsx/ConfigureObjectsx.h\ + ./ConfigureObjectsx/ConfigureObjectsx.cpp\ + ./ComputePressurex/ComputePressurex.h\ + ./ComputePressurex/ComputePressurex.cpp\ + ./BuildNodeSetsx/BuildNodeSetsx.h\ + ./BuildNodeSetsx/BuildNodeSetsx.cpp\ + ./BuildNodeSetsx/PartitionSets.cpp\ + ./SpcNodesx/SpcNodesx.h\ + ./SpcNodesx/SpcNodesx.cpp\ + ./MpcNodesx/MpcNodesx.h\ + ./MpcNodesx/MpcNodesx.cpp\ + ./DataInterpx/DataInterpx.cpp\ + ./DataInterpx/DataInterpx.h\ + ./HoleFillerx/HoleFillerx.cpp\ + ./HoleFillerx/HoleFillerx.h\ + ./MeshPartitionx/MeshPartitionx.cpp\ + ./MeshPartitionx/MeshPartitionx.h\ + ./ContourToMeshx/ContourToMeshx.cpp\ + ./ContourToMeshx/ContourToMeshx.h\ + ./Reducevectorgtosx/Reducevectorgtosx.cpp\ + ./Reducevectorgtosx/Reducevectorgtosx.h\ + ./Reducematrixfromgtofx/Reducematrixfromgtofx.cpp\ + ./Reducematrixfromgtofx/Reducematrixfromgton.cpp\ + ./Reducematrixfromgtofx/Reducematrixfromgtofx.h\ + ./Reduceloadfromgtofx/Reduceloadfromgtofx.h\ + ./Reduceloadfromgtofx/Reduceloadfromgtofx.cpp\ + ./NormalizeConstraintsx/NormalizeConstraintsx.cpp\ + ./NormalizeConstraintsx/NormalizeConstraintsx.h\ + ./SystemMatricesx/SystemMatricesx.cpp\ + ./SystemMatricesx/SystemMatricesx.h\ + ./PenaltyConstraintsx/PenaltyConstraintsx.cpp\ + ./PenaltyConstraintsx/PenaltyConstraintsx.h\ + ./PenaltySystemMatricesx/PenaltySystemMatricesx.cpp\ + ./PenaltySystemMatricesx/PenaltySystemMatricesx.h\ + ./Solverx/Solverx.cpp\ + ./Solverx/Solverx.h\ + ./Mergesolutionfromftogx/Mergesolutionfromftogx.cpp\ + ./Mergesolutionfromftogx/Mergesolutionfromftogx.h\ + ./ProcessParamsx/ProcessParamsx.cpp\ + ./ProcessParamsx/ProcessParamsx.h\ + ./VelocityExtrudex/VelocityExtrudex.cpp\ + ./VelocityExtrudex/VelocityExtrudex.h\ + ./SlopeExtrudex/SlopeExtrudex.cpp\ + ./SlopeExtrudex/SlopeExtrudex.h + + + + +libISSM_a_CXXFLAGS = -fPIC -DMATLAB -D_SERIAL_ -ansi -D_GNU_SOURCE -fno-omit-frame-pointer -pthread -D_CPP_ +if LARGEARRAYS +libISSM_a_CXXFLAGS += -D__GCC4BUILD__ +else +libISSM_a_CXXFLAGS += -DMX_COMPAT_32 +endif + + + + +#Parallel compilation +libpISSM_a_SOURCES = ./objects/objects.h\ + ./objects/Object.h\ + ./objects/Element.h\ + ./objects/Element.cpp\ + ./objects/Material.h\ + ./objects/Material.cpp\ + ./objects/Load.h\ + ./objects/Load.cpp\ + ./objects/SolverEnum.h\ + ./objects/Contour.h\ + ./objects/Contour.cpp\ + ./objects/OptArgs.h\ + ./objects/OptPars.h\ + ./objects/Friction.h\ + ./objects/Friction.cpp\ + ./objects/Node.h\ + ./objects/Node.cpp\ + ./objects/Tria.h\ + ./objects/Tria.cpp\ + ./objects/Sing.h\ + ./objects/Sing.cpp\ + ./objects/Beam.h\ + ./objects/Beam.cpp\ + ./objects/Penta.h\ + ./objects/Penta.cpp\ + ./objects/Matice.h\ + ./objects/Matice.cpp\ + ./objects/Matpar.h\ + ./objects/Matpar.cpp\ + ./objects/Input.h\ + ./objects/Input.cpp\ + ./objects/ParameterInputs.h\ + ./objects/ParameterInputs.cpp\ + ./objects/Spc.cpp\ + ./objects/Spc.h\ + ./objects/Rgb.cpp\ + ./objects/Rgb.h\ + ./objects/Penpair.cpp\ + ./objects/Penpair.h\ + ./objects/Pengrid.cpp\ + ./objects/Pengrid.h\ + ./objects/Icefront.cpp\ + ./objects/Icefront.h\ + ./objects/Param.cpp\ + ./objects/Param.h\ + ./objects/NodeSets.cpp\ + ./objects/NodeSets.h\ + ./DataSet/DataSet.cpp\ + ./DataSet/DataSet.h\ + ./shared/shared.h\ + ./shared/Alloc/alloc.h\ + ./shared/Alloc/alloc.cpp\ + ./shared/Matlab/matlabshared.h\ + ./shared/Matlab/PrintfFunction.cpp\ + ./shared/Matlab/ModuleBoot.cpp\ + ./shared/Matlab/mxGetAssignedField.cpp\ + ./shared/Matlab/mxGetField.cpp\ + ./shared/Matlab/CheckNumMatlabArguments.cpp\ + ./shared/Matrix/matrix.h\ + ./shared/Matrix/MatrixUtils.cpp\ + ./shared/Dofs/dofs.h\ + ./shared/Dofs/dofsetgen.cpp\ + ./shared/Dofs/DistributeNumDofs.cpp\ + ./shared/Numerics/numerics.h\ + ./shared/Numerics/GaussPoints.h\ + ./shared/Numerics/GaussPoints.cpp\ + ./shared/Numerics/cross.cpp\ + ./shared/Numerics/norm.cpp\ + ./shared/Numerics/BrentSearch.cpp\ + ./shared/Numerics/OptFunc.cpp\ + ./shared/Numerics/extrema.cpp\ + ./shared/Exceptions/exceptions.h\ + ./shared/Exceptions/Exceptions.cpp\ + ./shared/Exceptions/exprintf.cpp\ + ./shared/Exp/exp.h\ + ./shared/Exp/IsInPoly.cpp\ + ./shared/Exp/IsInPolySerial.cpp\ + ./shared/Exp/DomainOutlineRead.cpp\ + ./shared/TriMesh/trimesh.h\ + ./shared/TriMesh/AssociateSegmentToElement.cpp\ + ./shared/TriMesh/GridInsideHole.cpp\ + ./shared/TriMesh/OrderSegments.cpp\ + ./shared/TriMesh/SplitMeshForRifts.cpp\ + ./shared/TriMesh/TriMeshUtils.cpp\ + ./shared/Sorting/binary_search.cpp\ + ./shared/Sorting/sorting.h\ + ./shared/Elements/elements.h\ + ./shared/Elements/ResolvePointers.cpp\ + ./shared/Elements/Paterson.cpp\ + ./shared/Elements/GetElementNodeData.cpp\ + ./toolkits/petsc\ + ./toolkits/petsc/patches\ + ./toolkits/petsc/patches/petscpatches.h\ + ./toolkits/petsc/patches/MatlabMatrixToPetscMatrix.cpp\ + ./toolkits/petsc/patches/MatlabVectorToPetscVector.cpp\ + ./toolkits/petsc/patches/PetscMatrixToMatlabMatrix.cpp\ + ./toolkits/petsc/patches/PetscVectorToMatlabVector.cpp\ + ./toolkits/petsc/patches/MatlabMatrixToDoubleMatrix.cpp\ + ./toolkits/petsc/patches/MatlabVectorToDoubleVector.cpp\ + ./toolkits/petsc/patches/PetscDetermineLocalSize.cpp\ + ./toolkits/petsc/patches/VecTranspose.cpp\ + ./toolkits/petsc/patches/VecToMPISerial.cpp\ + ./toolkits/petsc/patches/MatToSerial.cpp\ + ./toolkits/petsc/patches/VecMerge.cpp\ + ./toolkits/petsc/patches/NewVec.cpp\ + ./toolkits/petsc/patches/NewVecFromLocalSize.cpp\ + ./toolkits/petsc/patches/PetscOptionsDetermineSolverType.cpp\ + ./toolkits/petsc/patches/PetscOptionsInsertMultipleString.cpp\ + ./toolkits/petsc/patches/NewMat.cpp\ + ./toolkits/petsc/patches/VecFree.cpp\ + ./toolkits/petsc/patches/VecDuplicatePatch.cpp\ + ./toolkits/petsc/patches/KSPFree.cpp\ + ./toolkits/petsc/patches/ISFree.cpp\ + ./toolkits/petsc/patches/MatFree.cpp\ + ./toolkits/petsc/patches/GetOwnershipBoundariesFromRange.cpp\ + ./toolkits/petsc/patches/VecPartition.cpp\ + ./toolkits/petsc/patches/MatPartition.cpp\ + ./toolkits/petsc/patches/MatInvert.cpp\ + ./toolkits/petsc/patches/MatMultPatch.cpp\ + ./toolkits/petsc/petscincludes.h\ + ./toolkits/mpi/mpiincludes.h\ + ./toolkits/mpi/patches/mpipatches.h\ + ./toolkits/mpi/patches/MPI_Upperrow.cpp\ + ./toolkits/mpi/patches/MPI_Lowerrow.cpp\ + ./toolkits/mpi/patches/MPI_Boundariesfromrange.cpp\ + ./toolkits/metis/metisincludes.h\ + ./toolkits/triangle/triangleincludes.h\ + ./toolkits.h\ + ./io/io.h\ + ./io/FetchData.cpp\ + ./io/WriteData.cpp\ + ./io/WriteDataToDisk.cpp\ + ./io/SerialFetchData.cpp\ + ./io/SerialWriteData.cpp\ + ./io/ParallelFetchData.cpp\ + ./io/ParallelFetchInteger.cpp\ + ./io/ParallelFetchMat.cpp\ + ./io/ParallelFetchScalar.cpp\ + ./io/ParallelFetchString.cpp\ + ./io/ModelFetchData.cpp\ + ./io/WriteNodeSets.cpp\ + ./io/WriteParams.cpp\ + ./io/FetchNodeSets.cpp\ + ./io/FetchRifts.cpp\ + ./io/ParameterInputsInit.cpp\ + ./EnumDefinitions/EnumDefinitions.h\ + ./EnumDefinitions/EnumDefinitions.cpp\ + ./EnumDefinitions/AnalysisTypeAsEnum.cpp\ + ./ModelProcessorx/Model.h\ + ./ModelProcessorx/Model.cpp\ + ./ModelProcessorx/CreateDataSets.cpp\ ./ModelProcessorx/CreateParameters.cpp\ ./ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp\ @@ -510,6 +506,5 @@ bin_PROGRAMS = else -bin_PROGRAMS = diagnostic.exe control.exe -#thermalsteady.exe +bin_PROGRAMS = diagnostic.exe control.exe thermal.exe endif @@ -522,4 +517,4 @@ control_exe_CXXFLAGS= -fPIC -D_PARALLEL_ -thermalsteady_exe_SOURCES = parallel/thermalsteady.cpp -thermalsteady_exe_CXXFLAGS= -fPIC -D_PARALLEL_ +thermal_exe_SOURCES = parallel/thermal.cpp +thermal_exe_CXXFLAGS= -fPIC -D_PARALLEL_ Index: /issm/trunk/src/c/Misfitx/Misfitx.cpp =================================================================== --- /issm/trunk/src/c/Misfitx/Misfitx.cpp (revision 464) +++ /issm/trunk/src/c/Misfitx/Misfitx.cpp (revision 465) @@ -14,5 +14,5 @@ void Misfitx( double* pJ, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials, - double* u_g, double* u_g_obs, ParameterInputs* inputs,int analysis_type){ + double* u_g, double* u_g_obs, ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ /*output: */ @@ -24,5 +24,5 @@ /*Compute gradients: */ - elements->Misfit(&J, u_g,u_g_obs,inputs,analysis_type); + elements->Misfit(&J, u_g,u_g_obs,inputs,analysis_type,sub_analysis_type); /*Sum all J from all cpus of the cluster:*/ Index: /issm/trunk/src/c/Misfitx/Misfitx.h =================================================================== --- /issm/trunk/src/c/Misfitx/Misfitx.h (revision 464) +++ /issm/trunk/src/c/Misfitx/Misfitx.h (revision 465) @@ -10,5 +10,5 @@ /* local prototypes: */ void Misfitx( double* pJ, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials, - double* u_g, double* u_g_obs, ParameterInputs* inputs,int analysis_type); + double* u_g, double* u_g_obs, ParameterInputs* inputs,int analysis_type,int sub_analysis_type); #endif /* _MISFITX_H */ Index: /issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp (revision 464) +++ /issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp (revision 465) @@ -13,8 +13,9 @@ #include "../Model.h" -void CreateParametersControl(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount){ +void CreateParametersControl(DataSet** pparameters,Model* model,ConstDataHandle model_handle){ int i; + DataSet* parameters=NULL; Param* param = NULL; int count; @@ -29,6 +30,7 @@ double* vy_obs=NULL; - /*Get counter : */ - count=*pcount; + /*Get parameters: */ + parameters=*pparameters; + count=parameters->Size(); /*control_type: */ @@ -132,5 +134,5 @@ /*Assign output pointer: */ - *pcount=count; + *pparameters=parameters; } Index: sm/trunk/src/c/ModelProcessorx/CreateConstraints.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateConstraints.cpp (revision 464) +++ (revision ) @@ -1,48 +1,0 @@ -/*!\file: CreateConstraints.cpp - * \brief general driver for creating constraints - */ - -#ifdef HAVE_CONFIG_H - #include "config.h" -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif - -#undef __FUNCT__ -#define __FUNCT__ "CreateConstraints" - -#include "./Model.h" -#include "../shared/shared.h" - -void CreateConstraints(DataSet** pconstraints, Model* model,ConstDataHandle model_handle){ - - /*This is just a high level driver: */ - if ((strcmp(model->analysis_type,"diagnostic_horiz")==0)|| (strcmp(model->analysis_type,"control")==0)){ - CreateConstraintsDiagnosticHoriz(pconstraints,model,model_handle); - } - else if (strcmp(model->analysis_type,"diagnostic_vert")==0){ - CreateConstraintsDiagnosticVert(pconstraints,model,model_handle); - } - else if (strcmp(model->analysis_type,"diagnostic_stokes")==0){ - CreateConstraintsDiagnosticStokes(pconstraints,model,model_handle); - } - else if (strcmp(model->analysis_type,"diagnostic_hutter")==0){ - CreateConstraintsDiagnosticHutter(pconstraints,model,model_handle); - } - else if (strcmp(model->analysis_type,"slope_compute")==0){ - CreateConstraintsSlopeCompute(pconstraints,model,model_handle); - } - /* - else if (strcmp(model->analysis_type,"melting")==0){ - CreateConstraintsMelting(pconstraints,model,model_handle); - } - else if (strcmp(model->analysis_type,"prognostic")==0){ - CreateConstraintsPrognostic(pconstraints,model,model_handle); - } - else if ((strcmp(model->analysis_type,"thermalsteady")==0) || (strcmp(model->analysis_type,"thermaltransient")==0)){ - CreateConstraintsThermal(pconstraints,model,model_handle); - }*/ - else{ - throw ErrorException(__FUNCT__," analysis_type not supported yet!"); - } -} Index: /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp (revision 465) +++ /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp (revision 465) @@ -0,0 +1,104 @@ +/*!\file: CreateDataSets + * \brief general driver for creating all datasets that make a finite element model + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#undef __FUNCT__ +#define __FUNCT__ "CreateElementsNodesAndMaterials" + +#include "./Model.h" +#include "../shared/shared.h" + + +void CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,Model* model,ConstDataHandle model_handle){ + + + /*create parameters common to all solutions: */ + CreateParameters(pparameters,model,model_handle); + + /*This is just a high level driver: */ + if (strcmp(model->analysis_type,"control")==0){ + + CreateElementsNodesAndMaterialsDiagnosticHoriz(pelements,pnodes,pmaterials, model,model_handle); + CreateConstraintsDiagnosticHoriz(pconstraints,model,model_handle); + CreateLoadsDiagnosticHoriz(ploads,model,model_handle); + CreateParametersControl(pparameters,model,model_handle); + + } + else if (strcmp(model->analysis_type,"diagnostic")==0){ + + if (strcmp(model->sub_analysis_type,"horiz")==0){ + + CreateElementsNodesAndMaterialsDiagnosticHoriz(pelements,pnodes,pmaterials, model,model_handle); + CreateConstraintsDiagnosticHoriz(pconstraints,model,model_handle); + CreateLoadsDiagnosticHoriz(ploads,model,model_handle); + CreateParametersDiagnosticHoriz(pparameters,model,model_handle); + + } + else if (strcmp(model->sub_analysis_type,"vert")==0){ + + CreateElementsNodesAndMaterialsDiagnosticVert(pelements,pnodes,pmaterials, model,model_handle); + CreateConstraintsDiagnosticVert(pconstraints,model,model_handle); + CreateLoadsDiagnosticVert(ploads,model,model_handle); + + } + else if (strcmp(model->sub_analysis_type,"stokes")==0){ + + CreateElementsNodesAndMaterialsDiagnosticStokes(pelements,pnodes,pmaterials, model,model_handle); + CreateConstraintsDiagnosticStokes(pconstraints,model,model_handle); + CreateLoadsDiagnosticStokes(ploads,model,model_handle); + + } + else if (strcmp(model->sub_analysis_type,"hutter")==0){ + + CreateElementsNodesAndMaterialsDiagnosticHutter(pelements,pnodes,pmaterials, model,model_handle); + CreateConstraintsDiagnosticHutter(pconstraints,model,model_handle); + CreateLoadsDiagnosticHutter(ploads,model,model_handle); + + } + } + else if (strcmp(model->analysis_type,"slope_compute")==0){ + + CreateElementsNodesAndMaterialsSlopeCompute(pelements,pnodes,pmaterials, model,model_handle); + CreateConstraintsSlopeCompute(pconstraints,model,model_handle); + CreateLoadsSlopeCompute(ploads,model,model_handle); + + + } + else if (strcmp(model->analysis_type,"thermal")==0){ + + if ((strcmp(model->sub_analysis_type,"steady")==0) || (strcmp(model->sub_analysis_type,"transient")==0)){ + + //CreateElementsNodesAndMaterialsThermal(pelements,pnodes,pmaterials, model,model_handle); + //CreateConstraintsThermal(pconstraints,model,model_handle); + //CreateLoadsThermal(ploads,model,model_handle); + + } + else if (strcmp(model->sub_analysis_type,"melting")==0){ + + //CreateElementsNodesAndMaterialsMelting(pelements,pnodes,pmaterials, model,model_handle); + //CreateConstraintsMelting(pconstraints,model,model_handle); + //CreateLoadsMelting(ploads,model,model_handle); + + } + } + else if (strcmp(model->analysis_type,"prognostic")==0){ + + //CreateElementsNodesAndMaterialsPrognostic(pelements,pnodes,pmaterials, model,model_handle); + //CreateConstraintsPrognostic(pconstraints,model,model_handle); + //CreateLoadsPrognostic(ploads,model,model_handle); + } + else{ + throw ErrorException(__FUNCT__,exprintf("%s%s%s%s"," analysis_type: ",model->analysis_type," sub_analysis_type: ",model->sub_analysis_type," not supported yet!")); + } + + + + + +} Index: sm/trunk/src/c/ModelProcessorx/CreateElementsNodesAndMaterials.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateElementsNodesAndMaterials.cpp (revision 464) +++ (revision ) @@ -1,63 +1,0 @@ -/*!\file: CreateElementsNodesAndMaterials.cpp - * \brief general driver for creating elements, nodes and materials - */ - -#ifdef HAVE_CONFIG_H - #include "config.h" -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif - -#undef __FUNCT__ -#define __FUNCT__ "CreateElementsNodesAndMaterials" - -#include "./Model.h" -#include "../shared/shared.h" - - -void CreateElementsNodesAndMaterials(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle){ - - /*This is just a high level driver: */ - if ((strcmp(model->analysis_type,"diagnostic_horiz")==0)|| (strcmp(model->analysis_type,"control")==0)){ - - CreateElementsNodesAndMaterialsDiagnosticHoriz(pelements,pnodes,pmaterials, model,model_handle); - - } - else if ((strcmp(model->analysis_type,"diagnostic_vert")==0)){ - - CreateElementsNodesAndMaterialsDiagnosticVert(pelements,pnodes,pmaterials, model,model_handle); - - } - else if ((strcmp(model->analysis_type,"diagnostic_stokes")==0)){ - - CreateElementsNodesAndMaterialsDiagnosticStokes(pelements,pnodes,pmaterials, model,model_handle); - - } - else if ((strcmp(model->analysis_type,"diagnostic_hutter")==0)){ - - CreateElementsNodesAndMaterialsDiagnosticHutter(pelements,pnodes,pmaterials, model,model_handle); - - } - else if ((strcmp(model->analysis_type,"slope_compute")==0)){ - - CreateElementsNodesAndMaterialsSlopeCompute(pelements,pnodes,pmaterials, model,model_handle); - - } - /* - else if (strcmp(model->analysis_type,"melting")==0){ - - CreateElementsNodesAndMaterialsDataSetsMelting(pelements,pnodes,pmaterials, model,model->analysis_type); - } - else if (strcmp(model->analysis_type,"prognostic")==0){ - - CreateElementsNodesAndMaterialsDataSetsPrognostic(pelements,pnodes,pmaterials, model,model->analysis_type); - } - else if ((strcmp(model->analysis_type,"thermalsteady")==0) || (strcmp(model->analysis_type,"thermaltransient")==0)){ - - CreateElementsNodesAndMaterialsDataSetsThermal(pelements,pnodes,pmaterials, model,model->analysis_type); - }*/ - else{ - throw ErrorException(__FUNCT__,exprintf("%s%s%s"," analysis_type: ",model->analysis_type," not supported yet!")); - } - -} Index: sm/trunk/src/c/ModelProcessorx/CreateLoads.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateLoads.cpp (revision 464) +++ (revision ) @@ -1,57 +1,0 @@ -/*!\file: CreateLoads.cpp - * \brief create loads datasets. General driver. - */ - -#ifdef HAVE_CONFIG_H - #include "config.h" -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif - -#undef __FUNCT__ -#define __FUNCT__ "CreateLoads" - -#include "./Model.h" -#include "../shared/shared.h" - -void CreateLoads(DataSet** ploads, Model* model,ConstDataHandle model_handle){ - - /*This is just a high level driver: */ - if ((strcmp(model->analysis_type,"diagnostic_horiz")==0)|| (strcmp(model->analysis_type,"control")==0)){ - - CreateLoadsDiagnosticHoriz(ploads,model,model_handle); - } - else if (strcmp(model->analysis_type,"diagnostic_vert")==0){ - - CreateLoadsDiagnosticVert(ploads,model,model_handle); - } - else if (strcmp(model->analysis_type,"diagnostic_stokes")==0){ - - CreateLoadsDiagnosticStokes(ploads,model,model_handle); - } - else if (strcmp(model->analysis_type,"diagnostic_hutter")==0){ - - CreateLoadsDiagnosticHutter(ploads,model,model_handle); - } - else if (strcmp(model->analysis_type,"slope_compute")==0){ - - CreateLoadsSlopeCompute(ploads,model,model_handle); - } - /* - else if (strcmp(model->analysis_type,"melting")==0){ - - CreateLoadsMelting(ploads,model,model_handle); - } - else if (strcmp(model->analysis_type,"prognostic")==0){ - - CreateLoadsPrognostic(ploads,model,model_handle); - } - else if ((strcmp(model->analysis_type,"thermalsteady")==0) || (strcmp(model->analysis_type,"thermaltransient")==0)){ - - CreateLoadsThermal(ploads,model,model_handle); - - }*/ - else{ - throw ErrorException(__FUNCT__," analysis_type not supported yet!"); - } -} Index: /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp (revision 464) +++ /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp (revision 465) @@ -21,4 +21,5 @@ int count=1; int analysis_type; + int sub_analysis_type; int numberofdofspernode; int dim; @@ -29,8 +30,14 @@ //Get analysis_type: analysis_type=AnalysisTypeAsEnum(model->analysis_type); + sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type); count++; param= new Param(count,"analysis_type",INTEGER); param->SetInteger(analysis_type); + parameters->AddObject(param); + + count++; + param= new Param(count,"sub_analysis_type",INTEGER); + param->SetInteger(sub_analysis_type); parameters->AddObject(param); @@ -187,5 +194,5 @@ /*Deal with numberofdofspernode: */ - DistributeNumDofs(&numberofdofspernode,analysis_type); + DistributeNumDofs(&numberofdofspernode,analysis_type,sub_analysis_type); count++; @@ -194,11 +201,4 @@ parameters->AddObject(param); - /*Deal with control parameters: */ - if (analysis_type==ControlAnalysisEnum()){ - CreateParametersControl(parameters,model,model_handle,&count); - } - if (analysis_type==DiagnosticHorizAnalysisEnum()){ - CreateParametersDiagnosticHoriz(parameters,model,model_handle,&count); - } Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 464) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 465) @@ -38,4 +38,5 @@ int analysis_type; + int sub_analysis_type; /*output: */ @@ -154,4 +155,5 @@ /*Get analysis_type: */ analysis_type=AnalysisTypeAsEnum(model->analysis_type); + sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type); @@ -571,5 +573,5 @@ /*Get number of dofs per node: */ - DistributeNumDofs(&node_numdofs,analysis_type); + DistributeNumDofs(&node_numdofs,analysis_type,sub_analysis_type); for (i=0;inumberofnodes;i++){ Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateParametersDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateParametersDiagnosticHoriz.cpp (revision 464) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateParametersDiagnosticHoriz.cpp (revision 465) @@ -13,7 +13,8 @@ #include "../Model.h" -void CreateParametersDiagnosticHoriz(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount){ +void CreateParametersDiagnosticHoriz(DataSet** pparameters,Model* model,ConstDataHandle model_handle){ Param* param = NULL; + DataSet* parameters=NULL; int count; int i; @@ -23,9 +24,11 @@ double* vz=NULL; + /*recover parameters : */ + parameters=*pparameters; + + count=parameters->Size(); + /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ if (!model->ismacayealpattyn)return; - - /*Get counter : */ - count=*pcount; /*Get vx and vy: */ @@ -61,5 +64,5 @@ /*Assign output pointer: */ - *pcount=count; + *pparameters=parameters; } Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp (revision 464) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp (revision 465) @@ -40,4 +40,5 @@ int analysis_type; + int sub_analysis_type; /*output: */ @@ -106,4 +107,5 @@ /*Get analysis_type: */ analysis_type=AnalysisTypeAsEnum(model->analysis_type); + sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type); /*Width of elements: */ @@ -297,5 +299,5 @@ /*Get number of dofs per node: */ - DistributeNumDofs(&node_numdofs,analysis_type); + DistributeNumDofs(&node_numdofs,analysis_type,sub_analysis_type); for (i=0;inumberofnodes;i++){ Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp (revision 464) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp (revision 465) @@ -38,4 +38,5 @@ int analysis_type; + int sub_analysis_type; /*output: */ @@ -134,4 +135,5 @@ /*Get analysis_type: */ analysis_type=AnalysisTypeAsEnum(model->analysis_type); + sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type); /*Width of elements: */ @@ -428,5 +430,5 @@ /*Get number of dofs per node: */ - DistributeNumDofs(&node_numdofs,analysis_type); + DistributeNumDofs(&node_numdofs,analysis_type,sub_analysis_type); for (i=0;inumberofnodes;i++){ Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp (revision 464) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp (revision 465) @@ -37,4 +37,5 @@ int analysis_type; + int sub_analysis_type; /*output: */ @@ -143,4 +144,5 @@ /*Get analysis_type: */ analysis_type=AnalysisTypeAsEnum(model->analysis_type); + sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type); /*Width of elements: */ @@ -354,5 +356,5 @@ /*Get number of dofs per node: */ - DistributeNumDofs(&node_numdofs,analysis_type); + DistributeNumDofs(&node_numdofs,analysis_type,sub_analysis_type); for (i=0;inumberofnodes;i++){ Index: /issm/trunk/src/c/ModelProcessorx/Model.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Model.cpp (revision 464) +++ /issm/trunk/src/c/ModelProcessorx/Model.cpp (revision 465) @@ -37,4 +37,5 @@ model->meshtype=NULL; model->analysis_type=NULL; + model->sub_analysis_type=NULL; model->numberofelements=0; model->numberofnodes=0; @@ -243,4 +244,5 @@ xfree((void**)&model->meshtype); xfree((void**)&model->analysis_type); + xfree((void**)&model->sub_analysis_type); if(model->numrifts){ @@ -289,4 +291,5 @@ /*In ModelInit, we get all the data that is not difficult to get, and that is small: */ ModelFetchData((void**)&model->analysis_type,NULL,NULL,model_handle,"analysis_type","String",NULL); + ModelFetchData((void**)&model->sub_analysis_type,NULL,NULL,model_handle,"sub_analysis_type","String",NULL); ModelFetchData((void**)&model->meshtype,NULL,NULL,model_handle,"type","String",NULL); Index: /issm/trunk/src/c/ModelProcessorx/Model.h =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Model.h (revision 464) +++ /issm/trunk/src/c/ModelProcessorx/Model.h (revision 465) @@ -16,4 +16,5 @@ char* meshtype; char* analysis_type; + char* sub_analysis_type; char* solverstring; @@ -174,7 +175,5 @@ /*Creation of fem datasets: general drivers*/ - void CreateElementsNodesAndMaterials(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle); - void CreateConstraints(DataSet** pconstraints, Model* model,ConstDataHandle model_handle); - void CreateLoads(DataSet** ploads, Model* model, ConstDataHandle model_handle); + void CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,Model* model,ConstDataHandle model_handle); void CreateParameters(DataSet** pparameters,Model* model,ConstDataHandle model_handle); @@ -186,5 +185,5 @@ void CreateConstraintsDiagnosticHoriz(DataSet** pconstraints,Model* model,ConstDataHandle model_handle); void CreateLoadsDiagnosticHoriz(DataSet** ploads, Model* model, ConstDataHandle model_handle); - void CreateParametersDiagnosticHoriz(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount); + void CreateParametersDiagnosticHoriz(DataSet** pparameters,Model* model,ConstDataHandle model_handle); /*diagnostic vertical*/ @@ -192,5 +191,4 @@ void CreateConstraintsDiagnosticVert(DataSet** pconstraints,Model* model,ConstDataHandle model_handle); void CreateLoadsDiagnosticVert(DataSet** ploads, Model* model, ConstDataHandle model_handle); - void CreateParametersDiagnosticVert(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount); /*diagnostic hutter*/ @@ -198,5 +196,4 @@ void CreateConstraintsDiagnosticHutter(DataSet** pconstraints,Model* model,ConstDataHandle model_handle); void CreateLoadsDiagnosticHutter(DataSet** ploads, Model* model, ConstDataHandle model_handle); - void CreateParametersDiagnosticHutter(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount); /*diagnostic stokes*/ @@ -204,5 +201,4 @@ void CreateConstraintsDiagnosticStokes(DataSet** pconstraints,Model* model,ConstDataHandle model_handle); void CreateLoadsDiagnosticStokes(DataSet** ploads, Model* model, ConstDataHandle model_handle); - void CreateParametersDiagnosticStokes(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount); /*slope compute*/ @@ -210,8 +206,7 @@ void CreateConstraintsSlopeCompute(DataSet** pconstraints,Model* model,ConstDataHandle model_handle); void CreateLoadsSlopeCompute(DataSet** ploads, Model* model, ConstDataHandle model_handle); - void CreateParametersSlopeCompute(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount); /*control:*/ - void CreateParametersControl(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount); + void CreateParametersControl(DataSet** pparameters,Model* model,ConstDataHandle model_handle); Index: /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp (revision 464) +++ /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp (revision 465) @@ -36,4 +36,5 @@ int analysis_type; + int sub_analysis_type; /*output: */ @@ -131,4 +132,5 @@ /*Get analysis_type: */ analysis_type=AnalysisTypeAsEnum(model->analysis_type); + sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type); /*Width of elements: */ @@ -397,5 +399,5 @@ /*Get number of dofs per node: */ - DistributeNumDofs(&node_numdofs,analysis_type); + DistributeNumDofs(&node_numdofs,analysis_type,sub_analysis_type); for (i=0;inumberofnodes;i++){ Index: /issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.cpp =================================================================== --- /issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.cpp (revision 464) +++ /issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.cpp (revision 465) @@ -14,5 +14,5 @@ void PenaltyConstraintsx(int* pconverged, int* pnum_unstable_constraints, DataSet* elements,DataSet* nodes, - DataSet* loads,DataSet* materials, ParameterInputs* inputs,int analysis_type){ + DataSet* loads,DataSet* materials, ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ int i; @@ -31,8 +31,8 @@ * management routine, otherwise, skip : */ if (loads->RiftIsPresent()){ - loads->RiftConstraints(&num_unstable_constraints,inputs,analysis_type); + loads->RiftConstraints(&num_unstable_constraints,inputs,analysis_type,sub_analysis_type); } else if(loads->MeltingIsPresent()){ - loads->MeltingConstraints(&num_unstable_constraints,inputs,analysis_type); + loads->MeltingConstraints(&num_unstable_constraints,inputs,analysis_type,sub_analysis_type); } else{ Index: /issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.h =================================================================== --- /issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.h (revision 464) +++ /issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.h (revision 465) @@ -11,5 +11,5 @@ /* local prototypes: */ void PenaltyConstraintsx(int* pconverged, int* pnum_unstable_constraints, DataSet* elements,DataSet* nodes, - DataSet* loads,DataSet* materials, ParameterInputs* inputs,int analysis_type); + DataSet* loads,DataSet* materials, ParameterInputs* inputs,int analysis_type,int sub_analysis_type); #endif /* _PENALTYCONSTRAINTSX_H */ Index: /issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp =================================================================== --- /issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp (revision 464) +++ /issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp (revision 465) @@ -14,5 +14,5 @@ void PenaltySystemMatricesx(Mat Kgg, Vec pg,DataSet* elements,DataSet* nodes,DataSet* loads,DataSet* materials, - int kflag,int pflag,ParameterInputs* inputs,int analysis_type){ + int kflag,int pflag,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ int i; @@ -33,6 +33,6 @@ /*Add penalties to stiffnesses, from loads: */ - if(kflag)loads->PenaltyCreateKMatrix(Kgg,inputs,kmax,analysis_type); - if(pflag)loads->PenaltyCreatePVector(pg,inputs,kmax,analysis_type); + if(kflag)loads->PenaltyCreateKMatrix(Kgg,inputs,kmax,analysis_type,sub_analysis_type); + if(pflag)loads->PenaltyCreatePVector(pg,inputs,kmax,analysis_type,sub_analysis_type); /*Assemble matrices: */ Index: /issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.h =================================================================== --- /issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.h (revision 464) +++ /issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.h (revision 465) @@ -11,5 +11,5 @@ /* local prototypes: */ void PenaltySystemMatricesx(Mat Kgg, Vec pg,DataSet* elements,DataSet* nodes,DataSet* loads,DataSet* materials, - int kflag,int pflag,ParameterInputs* inputs,int analysis_type); + int kflag,int pflag,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); #endif /* _PENALTYSYSTEMMATRICESX_H */ Index: /issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp =================================================================== --- /issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp (revision 464) +++ /issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp (revision 465) @@ -54,9 +54,5 @@ - if ( (analysis_type==ControlAnalysisEnum()) || - (analysis_type==DiagnosticHorizAnalysisEnum()) || - (analysis_type==DiagnosticVertAnalysisEnum()) || - (analysis_type==DiagnosticStokesAnalysisEnum()) - ){ + if ( (analysis_type==ControlAnalysisEnum()) || (analysis_type==DiagnosticAnalysisEnum())){ parameters->FindParam((void*)&vx,"vx"); Index: /issm/trunk/src/c/SystemMatricesx/SystemMatricesx.cpp =================================================================== --- /issm/trunk/src/c/SystemMatricesx/SystemMatricesx.cpp (revision 464) +++ /issm/trunk/src/c/SystemMatricesx/SystemMatricesx.cpp (revision 465) @@ -14,5 +14,5 @@ void SystemMatricesx(Mat* pKgg, Vec* ppg,DataSet* elements,DataSet* nodes,DataSet* loads,DataSet* materials, - int kflag,int pflag,int connectivity,int numberofdofspernode,ParameterInputs* inputs,int analysis_type){ + int kflag,int pflag,int connectivity,int numberofdofspernode,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ int i; @@ -40,10 +40,10 @@ /*Fill stiffness matrix and right hand side vector, from elements: */ - if(kflag)elements->CreateKMatrix(Kgg,inputs,analysis_type); - if(pflag)elements->CreatePVector(pg,inputs,analysis_type); + if(kflag)elements->CreateKMatrix(Kgg,inputs,analysis_type,sub_analysis_type); + if(pflag)elements->CreatePVector(pg,inputs,analysis_type,sub_analysis_type); /*Fill stiffness matrix and right hand side vector, from loads: */ - if(kflag)loads->CreateKMatrix(Kgg,inputs,analysis_type); - if(pflag)loads->CreatePVector(pg,inputs,analysis_type); + if(kflag)loads->CreateKMatrix(Kgg,inputs,analysis_type,sub_analysis_type); + if(pflag)loads->CreatePVector(pg,inputs,analysis_type,sub_analysis_type); /*Assemble matrices: */ Index: /issm/trunk/src/c/SystemMatricesx/SystemMatricesx.h =================================================================== --- /issm/trunk/src/c/SystemMatricesx/SystemMatricesx.h (revision 464) +++ /issm/trunk/src/c/SystemMatricesx/SystemMatricesx.h (revision 465) @@ -11,5 +11,5 @@ /* local prototypes: */ void SystemMatricesx(Mat* pKgg, Vec* ppg,DataSet* elements,DataSet* nodes,DataSet* loads,DataSet* materials, - int kflag,int pflag,int connectivity,int numberofdofspernode,ParameterInputs* inputs,int analysis_type); + int kflag,int pflag,int connectivity,int numberofdofspernode,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); #endif /* _SYSTEMMATRICESX_H */ Index: /issm/trunk/src/c/io/SerialFetchData.cpp =================================================================== --- /issm/trunk/src/c/io/SerialFetchData.cpp (revision 464) +++ /issm/trunk/src/c/io/SerialFetchData.cpp (revision 465) @@ -27,12 +27,12 @@ int outdataset_size; - double* outmatrix=NULL; - Mat outpetscmatrix=NULL; + double* outmatrix; + Mat outpetscmatrix; double* outmatrix_workspace=NULL;; int outmatrix_rows,outmatrix_cols; int petsc; - double* outvector=NULL; - Vec outpetscvector=NULL; + double* outvector; + Vec outpetscvector; double* outvector_workspace=NULL; int outvector_rows; Index: /issm/trunk/src/c/objects/Beam.cpp =================================================================== --- /issm/trunk/src/c/objects/Beam.cpp (revision 464) +++ /issm/trunk/src/c/objects/Beam.cpp (revision 465) @@ -212,11 +212,15 @@ #define __FUNCT__ "Beam::CreateKMatrix" -void Beam::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){ +void Beam::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ - if (analysis_type==DiagnosticHutterAnalysisEnum()) { - - CreateKMatrixDiagnosticHutter( Kgg,inputs,analysis_type); - + if (analysis_type==DiagnosticAnalysisEnum()) { + + if (sub_analysis_type==HutterAnalysisEnum()) { + + CreateKMatrixDiagnosticHutter( Kgg,inputs,analysis_type,sub_analysis_type); + } + else + throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis_type: ",sub_analysis_type," not supported yet")); } else{ @@ -230,5 +234,5 @@ #define __FUNCT__ "Beam::CreateKMatrixDiagnosticHutter" -void Beam::CreateKMatrixDiagnosticHutter(Mat Kgg,void* vinputs,int analysis_type){ +void Beam::CreateKMatrixDiagnosticHutter(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ @@ -265,9 +269,13 @@ #undef __FUNCT__ #define __FUNCT__ "Beam::CreatePVector" -void Beam::CreatePVector(Vec pg,void* inputs,int analysis_type){ +void Beam::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){ /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */ - if (analysis_type==DiagnosticHutterAnalysisEnum()) { - CreatePVectorDiagnosticHutter( pg,inputs,analysis_type); + if (analysis_type==DiagnosticAnalysisEnum()) { + if (sub_analysis_type==HutterAnalysisEnum()) { + CreatePVectorDiagnosticHutter( pg,inputs,analysis_type,sub_analysis_type); + } + else + throw ErrorException(__FUNCT__,exprintf("%s%i%s"," analysis ",analysis_type," not supported yet")); } else{ @@ -281,5 +289,5 @@ #define __FUNCT__ "Beam::CreatePVectorDiagnosticHutter" -void Beam::CreatePVectorDiagnosticHutter( Vec pg, void* vinputs, int analysis_type){ +void Beam::CreatePVectorDiagnosticHutter( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){ int i,j,k; @@ -534,25 +542,25 @@ #undef __FUNCT__ #define __FUNCT__ "Beam::Du" -void Beam::Du(_p_Vec*, double*, double*, void*, int){ +void Beam::Du(_p_Vec*, double*, double*, void*, int,int){ throw ErrorException(__FUNCT__," not supported yet!"); } #undef __FUNCT__ #define __FUNCT__ "Beam::Gradj" -void Beam::Gradj(_p_Vec*, double*, double*, void*, int, char*){ +void Beam::Gradj(_p_Vec*, double*, double*, void*, int, int,char*){ throw ErrorException(__FUNCT__," not supported yet!"); } #undef __FUNCT__ #define __FUNCT__ "Beam::GradjDrag" -void Beam::GradjDrag(_p_Vec*, double*, double*, void*, int){ +void Beam::GradjDrag(_p_Vec*, double*, double*, void*, int,int ){ throw ErrorException(__FUNCT__," not supported yet!"); } #undef __FUNCT__ #define __FUNCT__ "Beam::GradjB" -void Beam::GradjB(_p_Vec*, double*, double*, void*, int){ +void Beam::GradjB(_p_Vec*, double*, double*, void*, int, int){ throw ErrorException(__FUNCT__," not supported yet!"); } #undef __FUNCT__ #define __FUNCT__ "Beam::Misfit" -double Beam::Misfit(double*, double*, void*, int){ +double Beam::Misfit(double*, double*, void*, int,int ){ throw ErrorException(__FUNCT__," not supported yet!"); } Index: /issm/trunk/src/c/objects/Beam.h =================================================================== --- /issm/trunk/src/c/objects/Beam.h (revision 464) +++ /issm/trunk/src/c/objects/Beam.h (revision 465) @@ -56,12 +56,12 @@ int MyRank(); void Configure(void* loads,void* nodes,void* materials); - void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type); - void CreatePVector(Vec pg, void* inputs, int analysis_type); + void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); + void CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type); void UpdateFromInputs(void* inputs); void GetDofList(int* doflist,int* pnumberofdofs); void GetDofList1(int* doflist); - void CreateKMatrixDiagnosticHutter(Mat Kgg,void* inputs,int analysis_type); + void CreateKMatrixDiagnosticHutter(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); void GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3); - void CreatePVectorDiagnosticHutter(Vec pg,void* inputs,int analysis_type); + void CreatePVectorDiagnosticHutter(Vec pg,void* inputs,int analysis_type,int sub_analysis_type); void* GetMatPar(); @@ -78,9 +78,9 @@ void GetBedList(double*); void GetThicknessList(double* thickness_list); - void Du(_p_Vec*, double*, double*, void*, int); - void Gradj(_p_Vec*, double*, double*, void*, int, char*); - void GradjDrag(_p_Vec*, double*, double*, void*, int); - void GradjB(_p_Vec*, double*, double*, void*, int); - double Misfit(double*, double*, void*, int); + void Du(_p_Vec*, double*, double*, void*, int,int); + void Gradj(_p_Vec*, double*, double*, void*, int, int,char*); + void GradjDrag(_p_Vec*, double*, double*, void*, int,int ); + void GradjB(_p_Vec*, double*, double*, void*, int,int ); + double Misfit(double*, double*, void*, int,int ); void GetNodalFunctions(double* l1l2, double gauss_coord); void GetParameterValue(double* pvalue, double* value_list,double gauss_coord); Index: /issm/trunk/src/c/objects/Element.h =================================================================== --- /issm/trunk/src/c/objects/Element.h (revision 464) +++ /issm/trunk/src/c/objects/Element.h (revision 465) @@ -24,6 +24,6 @@ virtual void Demarshall(char** pmarshalled_dataset)=0; virtual void Configure(void* loads,void* nodes,void* materials)=0; - virtual void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type)=0; - virtual void CreatePVector(Vec pg, void* inputs, int analysis_type)=0; + virtual void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type)=0; + virtual void CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type)=0; virtual void UpdateFromInputs(void* inputs)=0; virtual void GetNodes(void** nodes)=0; @@ -33,9 +33,9 @@ virtual void GetThicknessList(double* thickness_list)=0; virtual void GetBedList(double* bed_list)=0; - virtual void Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type)=0; - virtual void Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,char* control_type)=0; - virtual void GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type)=0; - virtual void GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type)=0; - virtual double Misfit(double* u_g,double* u_g_obs,void* inputs,int analysis_type)=0; + virtual void Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type)=0; + virtual void Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type)=0; + virtual void GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type)=0; + virtual void GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type)=0; + virtual double Misfit(double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type)=0; virtual void ComputePressure(Vec p_g)=0; Index: /issm/trunk/src/c/objects/Icefront.cpp =================================================================== --- /issm/trunk/src/c/objects/Icefront.cpp (revision 464) +++ /issm/trunk/src/c/objects/Icefront.cpp (revision 465) @@ -199,5 +199,5 @@ return my_rank; } -void Icefront::DistributeNumDofs(int* numdofspernode,int analysis_type){return;} +void Icefront::DistributeNumDofs(int* numdofspernode,int analysis_type,int sub_analysis_type){return;} #undef __FUNCT__ @@ -240,5 +240,5 @@ #define __FUNCT__ "Icefront::CreateKMatrix" -void Icefront::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){ +void Icefront::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ /*No stiffness loads applied, do nothing: */ @@ -250,13 +250,26 @@ #undef __FUNCT__ #define __FUNCT__ "Icefront::CreatePVector" -void Icefront::CreatePVector(Vec pg, void* inputs, int analysis_type){ +void Icefront::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){ /*Just branch to the correct element icefront vector generator, according to the type of analysis we are carrying out: */ - if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){ - CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type); - } - else if ((analysis_type==DiagnosticStokesAnalysisEnum())){ - - CreatePVectorDiagnosticStokes( pg,inputs,analysis_type); + if (analysis_type==ControlAnalysisEnum()){ + + CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type); + + } + else if (analysis_type==DiagnosticAnalysisEnum()){ + + if (sub_analysis_type==HorizAnalysisEnum()){ + + CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type); + + } + else if (sub_analysis_type==StokesAnalysisEnum()){ + + CreatePVectorDiagnosticStokes( pg,inputs,analysis_type,sub_analysis_type); + + } + else throw ErrorException(__FUNCT__,exprintf("%s%i%s"," sub_analysis_type: ",sub_analysis_type," not supported yet")); + } else{ @@ -268,12 +281,12 @@ #undef __FUNCT__ #define __FUNCT__ "Icefront::CreatePVectorDiagnosticHoriz" -void Icefront::CreatePVectorDiagnosticHoriz( Vec pg, void* inputs, int analysis_type){ +void Icefront::CreatePVectorDiagnosticHoriz( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){ /*Branck on the type of icefront: */ if (strcmp(type,"segment")==0){ - CreatePVectorDiagnosticHorizSegment(pg,inputs,analysis_type); + CreatePVectorDiagnosticHorizSegment(pg,inputs,analysis_type,sub_analysis_type); } else{ - CreatePVectorDiagnosticHorizQuad(pg,inputs,analysis_type); + CreatePVectorDiagnosticHorizQuad(pg,inputs,analysis_type,sub_analysis_type); } } @@ -282,5 +295,5 @@ #define __FUNCT__ "Icefront::CreatePVectorDiagnosticHorizSegment" -void Icefront::CreatePVectorDiagnosticHorizSegment( Vec pg,void* vinputs, int analysis_type){ +void Icefront::CreatePVectorDiagnosticHorizSegment( Vec pg,void* vinputs, int analysis_type,int sub_analysis_type){ int i,j; @@ -395,5 +408,5 @@ #undef __FUNCT__ #define __FUNCT__ "Icefont::CreatePVectorDiagnosticHorizQuad" -void Icefront::CreatePVectorDiagnosticHorizQuad( Vec pg,void* vinputs, int analysis_type){ +void Icefront::CreatePVectorDiagnosticHorizQuad( Vec pg,void* vinputs, int analysis_type,int sub_analysis_type){ int i,j; @@ -557,5 +570,5 @@ #undef __FUNCT__ #define __FUNCT__ "Icefont::CreatePVectorDiagnosticStokes" -void Icefront::CreatePVectorDiagnosticStokes( Vec pg,void* vinputs, int analysis_type){ +void Icefront::CreatePVectorDiagnosticStokes( Vec pg,void* vinputs, int analysis_type,int sub_analysis_type){ int i,j; @@ -853,5 +866,5 @@ #undef __FUNCT__ #define __FUNCT__ "Icefront::PenaltyCreateKMatrix" -void Icefront::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type){ +void Icefront::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){ /*do nothing: */ } @@ -859,5 +872,5 @@ #undef __FUNCT__ #define __FUNCT__ "Icefront::PenaltyCreatePVector" -void Icefront::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type){ +void Icefront::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){ /*do nothing: */ } Index: /issm/trunk/src/c/objects/Icefront.h =================================================================== --- /issm/trunk/src/c/objects/Icefront.h (revision 464) +++ /issm/trunk/src/c/objects/Icefront.h (revision 465) @@ -55,14 +55,14 @@ int GetId(); int MyRank(); - void DistributeNumDofs(int* numdofspernode,int analysis_type); + void DistributeNumDofs(int* numdofspernode,int analysis_type,int sub_analysis_type); void Configure(void* elements,void* nodes,void* materials); - void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type); + void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); void UpdateFromInputs(void* inputs); - void CreatePVector(Vec pg, void* inputs, int analysis_type); - void CreatePVectorDiagnosticHoriz( Vec pg, void* inputs, int analysis_type); - void CreatePVectorDiagnosticHorizSegment( Vec pg,void* inputs, int analysis_type); - void CreatePVectorDiagnosticHorizQuad( Vec pg,void* inputs, int analysis_type); - void CreatePVectorDiagnosticStokes( Vec pg,void* inputs, int analysis_type); + void CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type); + void CreatePVectorDiagnosticHoriz( Vec pg, void* inputs, int analysis_type,int sub_analysis_type); + void CreatePVectorDiagnosticHorizSegment( Vec pg,void* inputs, int analysis_type,int sub_analysis_type); + void CreatePVectorDiagnosticHorizQuad( Vec pg,void* inputs, int analysis_type,int sub_analysis_type); + void CreatePVectorDiagnosticStokes( Vec pg,void* inputs, int analysis_type,int sub_analysis_type); void GetDofList(int* doflist,int* pnumberofdofs); void SegmentPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, double* normal,double length,int fill); @@ -71,6 +71,6 @@ void QuadPressureLoadStokes(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list, int fill); - void PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type); - void PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type); + void PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type); + void PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type); Object* copy(); }; Index: /issm/trunk/src/c/objects/Load.h =================================================================== --- /issm/trunk/src/c/objects/Load.h (revision 464) +++ /issm/trunk/src/c/objects/Load.h (revision 465) @@ -24,9 +24,9 @@ virtual void Demarshall(char** pmarshalled_dataset)=0; virtual void Configure(void* elements,void* nodes,void* materials)=0; - virtual void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type)=0; - virtual void CreatePVector(Vec pg, void* inputs, int analysis_type)=0; + virtual void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type)=0; + virtual void CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type)=0; virtual void UpdateFromInputs(void* inputs)=0; - virtual void PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type)=0; - virtual void PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type)=0; + virtual void PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type)=0; + virtual void PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type)=0; int Enum(); Index: /issm/trunk/src/c/objects/OptArgs.h =================================================================== --- /issm/trunk/src/c/objects/OptArgs.h (revision 464) +++ /issm/trunk/src/c/objects/OptArgs.h (revision 465) @@ -20,4 +20,5 @@ mxArray* n; mxArray* analysis_type; + mxArray* sub_analysis_type; }; Index: /issm/trunk/src/c/objects/Pengrid.cpp =================================================================== --- /issm/trunk/src/c/objects/Pengrid.cpp (revision 464) +++ /issm/trunk/src/c/objects/Pengrid.cpp (revision 465) @@ -135,5 +135,5 @@ return my_rank; } -void Pengrid::DistributeNumDofs(int* numdofpernode,int analysis_type){return;} +void Pengrid::DistributeNumDofs(int* numdofpernode,int analysis_type,int sub_analysis_type){return;} #undef __FUNCT__ @@ -156,5 +156,5 @@ #define __FUNCT__ "Pengrid::CreateKMatrix" -void Pengrid::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){ +void Pengrid::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ /*No loads applied, do nothing: */ @@ -165,5 +165,5 @@ #undef __FUNCT__ #define __FUNCT__ "Pengrid::CreatePVector" -void Pengrid::CreatePVector(Vec pg, void* inputs, int analysis_type){ +void Pengrid::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){ /*No loads applied, do nothing: */ @@ -179,13 +179,13 @@ #undef __FUNCT__ #define __FUNCT__ "Pengrid::PenaltyCreateKMatrix" -void Pengrid::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type){ - - if ((analysis_type==DiagnosticStokesAnalysisEnum())){ - - PenaltyCreateKMatrixDiagnosticStokes( Kgg,inputs,kmax,analysis_type); +void Pengrid::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){ + + if ((analysis_type==DiagnosticAnalysisEnum()) && ((sub_analysis_type==StokesAnalysisEnum()))){ + + PenaltyCreateKMatrixDiagnosticStokes( Kgg,inputs,kmax,analysis_type,sub_analysis_type); } else{ - throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet")); + throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet")); } @@ -194,5 +194,5 @@ #undef __FUNCT__ #define __FUNCT__ "Pengrid::PenaltyCreateKMatrixDiagnosticStokes" -void Pengrid::PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,void* vinputs,double kmax,int analysis_type){ +void Pengrid::PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){ const int numgrids=1; Index: /issm/trunk/src/c/objects/Pengrid.h =================================================================== --- /issm/trunk/src/c/objects/Pengrid.h (revision 464) +++ /issm/trunk/src/c/objects/Pengrid.h (revision 465) @@ -39,12 +39,12 @@ int GetId(); int MyRank(); - void DistributeNumDofs(int* numdofspernode,int analysis_type); + void DistributeNumDofs(int* numdofspernode,int analysis_type,int sub_analysis_type); void Configure(void* elements,void* nodes,void* materials); - void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type); - void CreatePVector(Vec pg, void* inputs, int analysis_type); + void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); + void CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type); void UpdateFromInputs(void* inputs); - void PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type); - void PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type); - void PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,void* inputs,double kmax,int analysis_type); + void PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type); + void PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type); + void PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type); void GetDofList(int* doflist,int* pnumberofdofspernode); Object* copy(); Index: /issm/trunk/src/c/objects/Penpair.cpp =================================================================== --- /issm/trunk/src/c/objects/Penpair.cpp (revision 464) +++ /issm/trunk/src/c/objects/Penpair.cpp (revision 465) @@ -182,5 +182,5 @@ return my_rank; } -void Penpair::DistributeNumDofs(int* numdofspernode,int analysis_type){return;} +void Penpair::DistributeNumDofs(int* numdofspernode,int analysis_type,int sub_analysis_type){return;} #undef __FUNCT__ @@ -208,5 +208,5 @@ #define __FUNCT__ "Penpair::CreateKMatrix" -void Penpair::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){ +void Penpair::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ /*No loads applied, do nothing: */ @@ -217,5 +217,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penpair::CreatePVector" -void Penpair::CreatePVector(Vec pg, void* inputs, int analysis_type){ +void Penpair::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){ /*No loads applied, do nothing: */ @@ -231,5 +231,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penpair::PenaltyCreateKMatrix" -void Penpair::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type){ +void Penpair::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){ if(numdofs==1){ @@ -245,5 +245,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penpair::PenaltyCreatePVector" -void Penpair::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type){ +void Penpair::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){ if(numdofs==1){ Index: /issm/trunk/src/c/objects/Penpair.h =================================================================== --- /issm/trunk/src/c/objects/Penpair.h (revision 464) +++ /issm/trunk/src/c/objects/Penpair.h (revision 465) @@ -51,11 +51,11 @@ int GetId(); int MyRank(); - void DistributeNumDofs(int* numdofspernode,int analysis_type); + void DistributeNumDofs(int* numdofspernode,int analysis_type,int sub_analysis_type); void Configure(void* elements,void* nodes,void* materials); - void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type); - void CreatePVector(Vec pg, void* inputs, int analysis_type); + void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); + void CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type); void UpdateFromInputs(void* inputs); - void PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type); - void PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type); + void PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type); + void PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type); Object* copy(); }; Index: /issm/trunk/src/c/objects/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Penta.cpp (revision 464) +++ /issm/trunk/src/c/objects/Penta.cpp (revision 465) @@ -284,30 +284,32 @@ #define __FUNCT__ "Penta::CreateKMatrix" -void Penta::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){ +void Penta::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ - if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){ - - CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type); - - } - else if ((analysis_type==DiagnosticVertAnalysisEnum())){ - - CreateKMatrixDiagnosticVert( Kgg,inputs,analysis_type); - - } - else if ((analysis_type==DiagnosticStokesAnalysisEnum())){ - - CreateKMatrixDiagnosticStokes( Kgg,inputs,analysis_type); - - } - else if ( - (analysis_type==SurfaceSlopeComputeXAnalysisEnum()) || - (analysis_type==SurfaceSlopeComputeYAnalysisEnum()) || - (analysis_type==BedSlopeComputeXAnalysisEnum()) || - (analysis_type==BedSlopeComputeYAnalysisEnum()) - ){ - - CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type); + if (analysis_type==ControlAnalysisEnum()){ + + CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type); + + } + else if (analysis_type==DiagnosticAnalysisEnum()){ + + if (sub_analysis_type==HorizAnalysisEnum()){ + + CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type); + } + else if (sub_analysis_type==VertAnalysisEnum()){ + + CreateKMatrixDiagnosticVert( Kgg,inputs,analysis_type,sub_analysis_type); + } + else if (sub_analysis_type==StokesAnalysisEnum()){ + + CreateKMatrixDiagnosticStokes( Kgg,inputs,analysis_type,sub_analysis_type); + + } + else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet")); + } + else if (analysis_type==SlopeComputeAnalysisEnum()){ + + CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type,sub_analysis_type); } else{ @@ -320,5 +322,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penta:CreateKMatrixDiagnosticHoriz" -void Penta::CreateKMatrixDiagnosticHoriz( Mat Kgg, void* vinputs, int analysis_type){ +void Penta::CreateKMatrixDiagnosticHoriz( Mat Kgg, void* vinputs, int analysis_type,int sub_analysis_type){ @@ -421,5 +423,5 @@ *and use its CreateKMatrix functionality to fill the global stiffness matrix: */ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrix(Kgg,inputs, analysis_type); + tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type); delete tria; return; @@ -552,5 +554,5 @@ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type); + tria->CreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type,sub_analysis_type); delete tria; } @@ -573,5 +575,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penta:CreateKMatrixDiagnosticVert" -void Penta::CreateKMatrixDiagnosticVert( Mat Kgg, void* vinputs, int analysis_type){ +void Penta::CreateKMatrixDiagnosticVert( Mat Kgg, void* vinputs, int analysis_type,int sub_analysis_type){ /* local declarations */ @@ -623,5 +625,5 @@ if(onsurface){ tria=(Tria*)SpawnTria(3,4,5); //nodes 3,4 and 5 are on the surface - tria->CreateKMatrixDiagnosticSurfaceVert(Kgg,inputs, analysis_type); + tria->CreateKMatrixDiagnosticSurfaceVert(Kgg,inputs, analysis_type,sub_analysis_type); delete tria; } @@ -710,5 +712,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penta:CreateKMatrixDiagnosticStokes" -void Penta::CreateKMatrixDiagnosticStokes( Mat Kgg, void* vinputs, int analysis_type){ +void Penta::CreateKMatrixDiagnosticStokes( Mat Kgg, void* vinputs, int analysis_type,int sub_analysis_type){ int i,j; @@ -985,27 +987,30 @@ #undef __FUNCT__ #define __FUNCT__ "Penta::CreatePVector" -void Penta::CreatePVector(Vec pg, void* inputs, int analysis_type){ +void Penta::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){ /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ - if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){ - - CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type); - } - else if ((analysis_type==DiagnosticVertAnalysisEnum())){ - - CreatePVectorDiagnosticVert( pg,inputs,analysis_type); - } - else if ((analysis_type==DiagnosticStokesAnalysisEnum())){ - - CreatePVectorDiagnosticStokes( pg,inputs,analysis_type); - } - else if ( - (analysis_type==SurfaceSlopeComputeXAnalysisEnum()) || - (analysis_type==SurfaceSlopeComputeYAnalysisEnum()) || - (analysis_type==BedSlopeComputeXAnalysisEnum()) || - (analysis_type==BedSlopeComputeYAnalysisEnum()) - ){ - - CreatePVectorSlopeCompute( pg,inputs,analysis_type); + if (analysis_type==ControlAnalysisEnum()){ + + CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type); + } + else if (analysis_type==DiagnosticAnalysisEnum()){ + + if (sub_analysis_type==HorizAnalysisEnum()){ + + CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type); + } + else if (sub_analysis_type==VertAnalysisEnum()){ + + CreatePVectorDiagnosticVert( pg,inputs,analysis_type,sub_analysis_type); + } + else if (sub_analysis_type==StokesAnalysisEnum()){ + + CreatePVectorDiagnosticStokes( pg,inputs,analysis_type,sub_analysis_type); + } + else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet")); + } + else if (analysis_type==SlopeComputeAnalysisEnum()){ + + CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type); } else{ @@ -1091,5 +1096,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penta::Du" -void Penta::Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type){ +void Penta::Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type){ int i; @@ -1122,11 +1127,11 @@ #undef __FUNCT__ #define __FUNCT__ "Penta::Gradj" -void Penta::Gradj(Vec grad_g,double* velocity,double* adjoint,void* inputs,int analysis_type,char* control_type){ +void Penta::Gradj(Vec grad_g,double* velocity,double* adjoint,void* inputs,int analysis_type,int sub_analysis_type,char* control_type){ if (strcmp(control_type,"drag")==0){ - GradjDrag( grad_g,velocity,adjoint,inputs,analysis_type); + GradjDrag( grad_g,velocity,adjoint,inputs,analysis_type,sub_analysis_type); } else if (strcmp(control_type,"B")==0){ - GradjB( grad_g, velocity, adjoint, inputs,analysis_type); + GradjB( grad_g, velocity, adjoint, inputs,analysis_type,sub_analysis_type); } else throw ErrorException(__FUNCT__,exprintf("%s%s","control type not supported yet: ",control_type)); @@ -1135,5 +1140,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penta::GradjDrag" -void Penta::GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type){ +void Penta::GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type){ @@ -1147,5 +1152,5 @@ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->GradjDrag( grad_g,u_g,lambda_g,inputs,analysis_type); + tria->GradjDrag( grad_g,u_g,lambda_g,inputs,analysis_type,sub_analysis_type); delete tria; return; @@ -1155,5 +1160,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penta::GradjB" -void Penta::GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type){ +void Penta::GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type){ throw ErrorException(__FUNCT__," not supported yet!"); } @@ -1161,5 +1166,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penta::Misfit" -double Penta::Misfit(double* velocity,double* obs_velocity,void* inputs,int analysis_type){ +double Penta::Misfit(double* velocity,double* obs_velocity,void* inputs,int analysis_type,int sub_analysis_type){ double J; @@ -1641,5 +1646,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penta::CreatePVectorDiagnosticHoriz" -void Penta::CreatePVectorDiagnosticHoriz( Vec pg, void* vinputs,int analysis_type){ +void Penta::CreatePVectorDiagnosticHoriz( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type){ int i,j; @@ -1707,5 +1712,5 @@ *and use its CreatePVector functionality to return an elementary load vector: */ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVector(pg,inputs, analysis_type); + tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type); delete tria; return; @@ -1997,5 +2002,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penta:CreatePVectorDiagnosticVert" -void Penta::CreatePVectorDiagnosticVert( Vec pg, void* vinputs,int analysis_type){ +void Penta::CreatePVectorDiagnosticVert( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type){ int i; @@ -2055,5 +2060,5 @@ if(onbed){ tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 are on the bedrock - tria->CreatePVectorDiagnosticBaseVert(pg,inputs, analysis_type); + tria->CreatePVectorDiagnosticBaseVert(pg,inputs, analysis_type,sub_analysis_type); delete tria; } @@ -2176,5 +2181,5 @@ #define __FUNCT__ "Penta::CreateKMatrixSlopeCompute" -void Penta::CreateKMatrixSlopeCompute(Mat Kgg,void* inputs,int analysis_type){ +void Penta::CreateKMatrixSlopeCompute(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ /*Collapsed formulation: */ @@ -2186,5 +2191,5 @@ /*Spawn Tria element from the base of the Penta: */ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrix(Kgg,inputs, analysis_type); + tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type); delete tria; return; @@ -2195,5 +2200,5 @@ #define __FUNCT__ "Penta::CreatePVectorSlopeCompute" -void Penta::CreatePVectorSlopeCompute( Vec pg, void* inputs, int analysis_type){ +void Penta::CreatePVectorSlopeCompute( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){ /*Collapsed formulation: */ @@ -2205,5 +2210,5 @@ /*Spawn Tria element from the base of the Penta: */ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVector(pg,inputs, analysis_type); + tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type); delete tria; return; @@ -2826,5 +2831,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penta::CreatePVectorDiagnosticStokes" -void Penta::CreatePVectorDiagnosticStokes( Vec pg, void* vinputs,int analysis_type){ +void Penta::CreatePVectorDiagnosticStokes( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type){ Index: /issm/trunk/src/c/objects/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Penta.h (revision 464) +++ /issm/trunk/src/c/objects/Penta.h (revision 465) @@ -73,8 +73,8 @@ int MyRank(); void Configure(void* loads,void* nodes,void* materials); - void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type); - void CreateKMatrixDiagnosticHoriz( Mat Kgg, void* inputs, int analysis_type); - void CreateKMatrixDiagnosticVert( Mat Kgg, void* inputs, int analysis_type); - void CreatePVector(Vec pg, void* inputs, int analysis_type); + void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); + void CreateKMatrixDiagnosticHoriz( Mat Kgg, void* inputs, int analysis_type,int sub_analysis_type); + void CreateKMatrixDiagnosticVert( Mat Kgg, void* inputs, int analysis_type,int sub_analysis_type); + void CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type); void UpdateFromInputs(void* inputs); void GetDofList(int* doflist,int* pnumberofdofs); @@ -84,9 +84,9 @@ void GetNodes(void** nodes); int GetOnBed(); - void Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type); - void Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,char* control_type); - void GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type); - void GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type); - double Misfit(double* u_g,double* u_g_obs,void* inputs,int analysis_type); + void Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type); + void Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type); + void GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type); + void GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type); + double Misfit(double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type); void GetThicknessList(double* thickness_list); @@ -105,6 +105,6 @@ void GetNodalFunctionsDerivativesParams(double* dl1dl2dl3dl4dl5dl6,double* gauss_l1l2l3l4); void GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss_l1l2l3l4); - void CreatePVectorDiagnosticHoriz( Vec pg, void* inputs,int analysis_type); - void CreatePVectorDiagnosticVert( Vec pg, void* inputs,int analysis_type); + void CreatePVectorDiagnosticHoriz( Vec pg, void* inputs,int analysis_type,int sub_analysis_type); + void CreatePVectorDiagnosticVert( Vec pg, void* inputs,int analysis_type,int sub_analysis_type); void GetParameterValue(double* pvalue, double* v_list,double* gauss_l1l2l3l4); void GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_l1l2l3l4); @@ -113,9 +113,9 @@ void SlopeExtrude(Vec sg,double* sg_serial); void ComputePressure(Vec p_g); - void CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type); - void CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type); + void CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type); + void CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type); - void CreateKMatrixDiagnosticStokes( Mat Kgg, void* vinputs, int analysis_type); - void CreatePVectorDiagnosticStokes( Vec pg, void* vinputs,int analysis_type); + void CreateKMatrixDiagnosticStokes( Mat Kgg, void* vinputs, int analysis_type,int sub_analysis_type); + void CreatePVectorDiagnosticStokes( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type); void ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp); void GetMatrixInvert(double* Ke_invert, double* Ke); Index: /issm/trunk/src/c/objects/Sing.cpp =================================================================== --- /issm/trunk/src/c/objects/Sing.cpp (revision 464) +++ /issm/trunk/src/c/objects/Sing.cpp (revision 465) @@ -193,14 +193,14 @@ #define __FUNCT__ "Sing::CreateKMatrix" -void Sing::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){ +void Sing::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ - if (analysis_type==DiagnosticHutterAnalysisEnum()){ - - CreateKMatrixDiagnosticHutter( Kgg,inputs,analysis_type); + if ((analysis_type==DiagnosticAnalysisEnum()) && (sub_analysis_type==HutterAnalysisEnum())){ + + CreateKMatrixDiagnosticHutter( Kgg,inputs,analysis_type,sub_analysis_type); } else{ - throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet")); + throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s\n","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet")); } @@ -211,5 +211,5 @@ #define __FUNCT__ "Sing::CreateKMatrixDiagnosticHutter" -void Sing::CreateKMatrixDiagnosticHutter(Mat Kgg,void* vinputs,int analysis_type){ +void Sing::CreateKMatrixDiagnosticHutter(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ const int numgrids=1; @@ -234,9 +234,11 @@ #undef __FUNCT__ #define __FUNCT__ "Sing::CreatePVector" -void Sing::CreatePVector(Vec pg,void* inputs,int analysis_type){ +void Sing::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){ /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */ - if (analysis_type==DiagnosticHutterAnalysisEnum()){ - CreatePVectorDiagnosticHutter( pg,inputs,analysis_type); + if ((analysis_type==DiagnosticAnalysisEnum()) && (sub_analysis_type==HutterAnalysisEnum())){ + + CreatePVectorDiagnosticHutter( pg,inputs,analysis_type,sub_analysis_type); + } else{ @@ -251,5 +253,5 @@ #define __FUNCT__ "Sing::CreatePVectorDiagnosticHutter" -void Sing::CreatePVectorDiagnosticHutter( Vec pg, void* vinputs, int analysis_type){ +void Sing::CreatePVectorDiagnosticHutter( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){ @@ -427,25 +429,25 @@ #undef __FUNCT__ #define __FUNCT__ "Sing::Du" -void Sing::Du(_p_Vec*, double*, double*, void*, int){ +void Sing::Du(_p_Vec*, double*, double*, void*, int,int){ throw ErrorException(__FUNCT__," not supported yet!"); } #undef __FUNCT__ #define __FUNCT__ "Sing::Gradj" -void Sing::Gradj(_p_Vec*, double*, double*, void*, int, char*){ +void Sing::Gradj(_p_Vec*, double*, double*, void*, int, int ,char*){ throw ErrorException(__FUNCT__," not supported yet!"); } #undef __FUNCT__ #define __FUNCT__ "Sing::GradjDrag" -void Sing::GradjDrag(_p_Vec*, double*, double*, void*, int){ +void Sing::GradjDrag(_p_Vec*, double*, double*, void*, int,int){ throw ErrorException(__FUNCT__," not supported yet!"); } #undef __FUNCT__ #define __FUNCT__ "Sing::GradjB" -void Sing::GradjB(_p_Vec*, double*, double*, void*, int){ +void Sing::GradjB(_p_Vec*, double*, double*, void*, int,int){ throw ErrorException(__FUNCT__," not supported yet!"); } #undef __FUNCT__ #define __FUNCT__ "Sing::Misfit" -double Sing::Misfit(double*, double*, void*, int){ +double Sing::Misfit(double*, double*, void*, int,int){ throw ErrorException(__FUNCT__," not supported yet!"); } Index: /issm/trunk/src/c/objects/Sing.h =================================================================== --- /issm/trunk/src/c/objects/Sing.h (revision 464) +++ /issm/trunk/src/c/objects/Sing.h (revision 465) @@ -51,12 +51,12 @@ int MyRank(); void Configure(void* loads,void* nodes,void* materials); - void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type); - void CreatePVector(Vec pg, void* inputs, int analysis_type); + void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); + void CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type); void UpdateFromInputs(void* inputs); void GetDofList(int* doflist,int* pnumberofdofs); void GetDofList1(int* doflist); - void CreateKMatrixDiagnosticHutter(Mat Kgg,void* inputs,int analysis_type); + void CreateKMatrixDiagnosticHutter(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); void GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3); - void CreatePVectorDiagnosticHutter(Vec pg,void* inputs,int analysis_type); + void CreatePVectorDiagnosticHutter(Vec pg,void* inputs,int analysis_type,int sub_analysis_type); void* GetMatPar(); @@ -73,9 +73,9 @@ void GetBedList(double*); void GetThicknessList(double* thickness_list); - void Du(_p_Vec*, double*, double*, void*, int); - void Gradj(_p_Vec*, double*, double*, void*, int, char*); - void GradjDrag(_p_Vec*, double*, double*, void*, int); - void GradjB(_p_Vec*, double*, double*, void*, int); - double Misfit(double*, double*, void*, int); + void Du(_p_Vec*, double*, double*, void*, int,int); + void Gradj(_p_Vec*, double*, double*, void*, int, int,char*); + void GradjDrag(_p_Vec*, double*, double*, void*, int,int); + void GradjB(_p_Vec*, double*, double*, void*, int,int); + double Misfit(double*, double*, void*, int,int); Index: /issm/trunk/src/c/objects/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Tria.cpp (revision 464) +++ /issm/trunk/src/c/objects/Tria.cpp (revision 465) @@ -266,20 +266,24 @@ #define __FUNCT__ "Tria::CreateKMatrix" -void Tria::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){ +void Tria::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ - if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){ - - CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type); - - } - else if ( - (analysis_type==SurfaceSlopeComputeXAnalysisEnum()) || - (analysis_type==SurfaceSlopeComputeYAnalysisEnum()) || - (analysis_type==BedSlopeComputeXAnalysisEnum()) || - (analysis_type==BedSlopeComputeYAnalysisEnum()) - ){ - - CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type); + if (analysis_type==ControlAnalysisEnum()){ + + CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type); + } + else if (analysis_type==DiagnosticAnalysisEnum()){ + + if (sub_analysis_type==HorizAnalysisEnum()){ + + CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type); + } + else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet")); + + } + else if (analysis_type==SlopeComputeAnalysisEnum()){ + + CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type,sub_analysis_type); + } else{ @@ -293,5 +297,5 @@ #define __FUNCT__ "Tria::CreateKMatrixDiagnosticHoriz" -void Tria::CreateKMatrixDiagnosticHoriz(Mat Kgg,void* vinputs,int analysis_type){ +void Tria::CreateKMatrixDiagnosticHoriz(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ @@ -477,5 +481,5 @@ /*Do not forget to include friction: */ if(!shelf){ - CreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type); + CreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type,sub_analysis_type); } @@ -508,5 +512,5 @@ #undef __FUNCT__ #define __FUNCT__ "Tria::CreateKMatrixDiagnosticHorizFriction" -void Tria::CreateKMatrixDiagnosticHorizFriction(Mat Kgg,void* vinputs,int analysis_type){ +void Tria::CreateKMatrixDiagnosticHorizFriction(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ @@ -677,5 +681,5 @@ #define __FUNCT__ "Tria::CreateKMatrixSlopeCompute" -void Tria::CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type){ +void Tria::CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ /* local declarations */ @@ -758,18 +762,23 @@ #undef __FUNCT__ #define __FUNCT__ "Tria::CreatePVector" -void Tria::CreatePVector(Vec pg,void* inputs,int analysis_type){ +void Tria::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){ /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */ - if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){ - CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type); - } - else if ( - (analysis_type==SurfaceSlopeComputeXAnalysisEnum()) || - (analysis_type==SurfaceSlopeComputeYAnalysisEnum()) || - (analysis_type==BedSlopeComputeXAnalysisEnum()) || - (analysis_type==BedSlopeComputeYAnalysisEnum()) - ){ - - CreatePVectorSlopeCompute( pg,inputs,analysis_type); + if (analysis_type==ControlAnalysisEnum()){ + + CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type); + + } + else if (analysis_type==DiagnosticAnalysisEnum()){ + if (sub_analysis_type==HorizAnalysisEnum()){ + + CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type); + + } + else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet")); + } + else if (analysis_type==SlopeComputeAnalysisEnum()){ + + CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type); } else{ @@ -779,10 +788,8 @@ } - - #undef __FUNCT__ #define __FUNCT__ "Tria::CreatePVectorDiagnosticHoriz" -void Tria::CreatePVectorDiagnosticHoriz( Vec pg, void* vinputs, int analysis_type){ +void Tria::CreatePVectorDiagnosticHoriz( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){ int i,j; @@ -952,5 +959,5 @@ #define __FUNCT__ "Tria::CreatePVectorSlopeCompute" -void Tria::CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type){ +void Tria::CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){ int i,j; @@ -992,8 +999,8 @@ for(i=0;ipenpair.numdofs==1){ - PenpairLoadPenaltyCreateKMatrixOneDof( pKe_gg,vpthis,inputs,K_flag,kmax,analysis_type); + PenpairLoadPenaltyCreateKMatrixOneDof( pKe_gg,vpthis,inputs,K_flag,kmax,analysis_type,sub_analysis_type); } else if (this->penpair.numdofs==2){ - PenpairLoadPenaltyCreateKMatrixTwoDof( pKe_gg,vpthis,inputs,K_flag,kmax,analysis_type); + PenpairLoadPenaltyCreateKMatrixTwoDof( pKe_gg,vpthis,inputs,K_flag,kmax,analysis_type,sub_analysis_type); } else{ @@ -719,5 +719,5 @@ #undef __FUNCT__ #define __FUNCT__ "PenpairLoadPenaltyCreateKMatrixOneDof" -int PenpairLoadPenaltyCreateKMatrixOneDof( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inptus, int K_flag, double kmax, int analysis_type){ +int PenpairLoadPenaltyCreateKMatrixOneDof( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inptus, int K_flag, double kmax, int analysis_type,int sub_analysis_type){ @@ -769,5 +769,5 @@ #undef __FUNCT__ #define __FUNCT__ "PenpairLoadPenaltyCreateKMatrixTwoDof" -int PenpairLoadPenaltyCreateKMatrixTwoDof( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, double kmax, int analysis_type){ +int PenpairLoadPenaltyCreateKMatrixTwoDof( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, double kmax, int analysis_type,int sub_analysis_type){ @@ -916,5 +916,5 @@ #define __FUNCT__ "PenpairLoadPenaltyCreatePVector" -int PenpairLoadPenaltyCreatePVector( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, double kmax, int analysis_type){ +int PenpairLoadPenaltyCreatePVector( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, double kmax, int analysis_type,int sub_analysis_type){ PenpairLoad* this = NULL; @@ -931,5 +931,5 @@ } else if (this->penpair.numdofs==2){ - PenpairLoadPenaltyCreatePVectorTwoDof( ppe_g,vpthis,inputs,kmax,analysis_type); + PenpairLoadPenaltyCreatePVectorTwoDof( ppe_g,vpthis,inputs,kmax,analysis_type,sub_analysis_type); } else{ @@ -947,5 +947,5 @@ #define __FUNCT__ "PenpairLoadPenaltyCreatePVectorTwoDof" -int PenpairLoadPenaltyCreatePVectorTwoDof( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, double kmax, int analysis_type){ +int PenpairLoadPenaltyCreatePVectorTwoDof( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, double kmax, int analysis_type,int sub_analysis_type){ int i,j; @@ -1092,5 +1092,4 @@ return noerr; } -int PotentialUnstableConstraints(DataSet* lst,ParameterInputs* inputs,int analysis_type_enum); /*-------------------------------------------------- @@ -1100,5 +1099,6 @@ #define __FUNCT__ "PenpairLoadPotentialUnstableConstraint" -int PenpairLoadPotentialUnstableConstraint(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type){ +int PotentialUnstableConstraints(DataSet* lst,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); +int PenpairLoadPotentialUnstableConstraint(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type,int sub_analysis_type){ /* vpthis for polymorphic function compatibility */ @@ -1186,5 +1186,5 @@ #define __FUNCT__ "PenpairLoadPenaltyPreConstrain" -int PenpairLoadPenaltyPreConstrain(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type){ +int PenpairLoadPenaltyPreConstrain(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type,int sub_analysis_type){ /* vpthis for polymorphic function compatibility */ @@ -1290,5 +1290,5 @@ #define __FUNCT__ "PenpairLoadPenaltyConstrain" -int PenpairLoadPenaltyConstrain(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type){ +int PenpairLoadPenaltyConstrain(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type,int sub_analysis_type){ /* vpthis for polymorphic function compatibility */ @@ -1417,5 +1417,5 @@ #define __FUNCT__ "PenpairLoadPenetration" -int PenpairLoadPenetration(double* ppenetration, void* vpthis,ParameterInputs* inputs, int analysis_type){ +int PenpairLoadPenetration(double* ppenetration, void* vpthis,ParameterInputs* inputs, int analysis_type,int sub_analysis_type){ /* vpthis for polymorphic function compatibility */ Index: /issm/trunk/src/c/objects/cielo/PentaElement.c =================================================================== --- /issm/trunk/src/c/objects/cielo/PentaElement.c (revision 464) +++ /issm/trunk/src/c/objects/cielo/PentaElement.c (revision 465) @@ -818,5 +818,5 @@ #undef __FUNCT__ #define __FUNCT__ "PentaElementCreateKMatrix" -int PentaElementCreateKMatrix( ElemMatrix* *pKe_gg, ElemMatrix** pKTe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, int KT_flag, int analysis_type){ +int PentaElementCreateKMatrix( ElemMatrix* *pKe_gg, ElemMatrix** pKTe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, int KT_flag, int analysis_type,int sub_analysis_type){ @@ -2564,5 +2564,5 @@ #define __FUNCT__ "PentaElementCreatePVector" -int PentaElementCreatePVector( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, int analysis_type){ +int PentaElementCreatePVector( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, int analysis_type,int sub_analysis_type){ int noerr=1; @@ -5002,5 +5002,5 @@ #undef __FUNCT__ #define __FUNCT__ "PentaElementCreateDuVector" -int PentaElementCreateDuVector( ElemVector* *pdue_g, void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type){ +int PentaElementCreateDuVector( ElemVector* *pdue_g, void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type,int sub_analysis_type){ int noerr=1; @@ -5077,5 +5077,5 @@ /*Ok, now triaelement is correctly configured, call on its method CreateDuVector: */ - noerr=TriaElementCreateDuVector( pdue_g, (void*)triaelement, velocity, obs_velocity, inputs, analysis_type); + noerr=TriaElementCreateDuVector( pdue_g, (void*)triaelement, velocity, obs_velocity, inputs, analysis_type,sub_analysis_type); /*Now delete tria and triaelement: */ @@ -5092,13 +5092,13 @@ #undef __FUNCT__ #define __FUNCT__ "PentaElementCreateGradVectors" -int PentaElementCreateGradjVectors( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,char* control_type){ +int PentaElementCreateGradjVectors( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,int sub_analysis_type,char* control_type){ int noerr=1; if strcasecmp_eq(control_type,"drag"){ - noerr=PentaElementCreateGradjVectorsDrag( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type); + noerr=PentaElementCreateGradjVectorsDrag( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type,sub_analysis_type); } else if strcasecmp_eq(control_type,"B"){ - noerr=PentaElementCreateGradjVectorsB( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type); + noerr=PentaElementCreateGradjVectorsB( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type,sub_analysis_type); } else{ @@ -5114,5 +5114,5 @@ #undef __FUNCT__ #define __FUNCT__ "PentaElementMisfit" -int PentaElementMisfit( double* pJelem,void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type){ +int PentaElementMisfit( double* pJelem,void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type,int sub_analysis_type){ int noerr=1; @@ -5189,5 +5189,5 @@ /*Ok, now triaelement is correctly configured, call on its method Misfit: */ - noerr=TriaElementMisfit(pJelem,(void*)triaelement, velocity, obs_velocity, inputs, analysis_type); + noerr=TriaElementMisfit(pJelem,(void*)triaelement, velocity, obs_velocity, inputs, analysis_type,sub_analysis_type); /*Now delete tria and triaelement: */ @@ -5204,5 +5204,5 @@ #undef __FUNCT__ #define __FUNCT__ "PentaElementCreateGradVectorsDrag" -int PentaElementCreateGradjVectorsDrag( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type){ +int PentaElementCreateGradjVectorsDrag( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ int noerr=1; @@ -5279,5 +5279,5 @@ /*Ok, now triaelement is correctly configured, call on its method CreateGradjVectorsdrag: */ - noerr=TriaElementCreateGradjVectorsDrag( pgradje_g, (void*)triaelement,velocity, adjoint, inputs,analysis_type); + noerr=TriaElementCreateGradjVectorsDrag( pgradje_g, (void*)triaelement,velocity, adjoint, inputs,analysis_type,sub_analysis_type); /*Now delete tria and triaelement: */ @@ -5295,5 +5295,5 @@ #undef __FUNCT__ #define __FUNCT__ "PentaElementCreateGradjVectorsB" -int PentaElementCreateGradjVectorsB( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type){ +int PentaElementCreateGradjVectorsB( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ printf("Not supported yet!\n"); Index: /issm/trunk/src/c/objects/cielo/TriaElement.c =================================================================== --- /issm/trunk/src/c/objects/cielo/TriaElement.c (revision 464) +++ /issm/trunk/src/c/objects/cielo/TriaElement.c (revision 465) @@ -777,5 +777,5 @@ #undef __FUNCT__ #define __FUNCT__ "TriaElementCreateKMatrix" -int TriaElementCreateKMatrix( ElemMatrix* *pKe_gg, ElemMatrix** pKTe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, int KT_flag, int analysis_type){ +int TriaElementCreateKMatrix( ElemMatrix* *pKe_gg, ElemMatrix** pKTe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, int KT_flag, int analysis_type,int sub_analysis_type){ @@ -800,5 +800,5 @@ #define __FUNCT__ "TriaElementCreatePVector" -int TriaElementCreatePVector( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, int analysis_type){ +int TriaElementCreatePVector( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, int analysis_type,int sub_analysis_type){ int noerr=1; @@ -2044,5 +2044,5 @@ #undef __FUNCT__ #define __FUNCT__ "TriaElementCreateDuVector" -int TriaElementCreateDuVector( ElemVector* *pdue_g, void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type){ +int TriaElementCreateDuVector( ElemVector* *pdue_g, void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type,int sub_analysis_type){ int i,j; @@ -2263,13 +2263,13 @@ #undef __FUNCT__ #define __FUNCT__ "TriaElementCreateGradVectors" -int TriaElementCreateGradjVectors( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,char* control_type){ +int TriaElementCreateGradjVectors( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,int sub_analysis_type,char* control_type){ int noerr=1; if strcasecmp_eq(control_type,"drag"){ - noerr=TriaElementCreateGradjVectorsDrag( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type); + noerr=TriaElementCreateGradjVectorsDrag( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type,sub_analysis_type); } else if strcasecmp_eq(control_type,"B"){ - noerr=TriaElementCreateGradjVectorsB( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type); + noerr=TriaElementCreateGradjVectorsB( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type,sub_analysis_type); } else{ @@ -2285,5 +2285,5 @@ #undef __FUNCT__ #define __FUNCT__ "TriaElementCreateGradVectorsDrag" -int TriaElementCreateGradjVectorsDrag( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type){ +int TriaElementCreateGradjVectorsDrag( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ int i,j; @@ -2566,5 +2566,5 @@ #undef __FUNCT__ #define __FUNCT__ "TriaElementCreateGradjVectorsB" -int TriaElementCreateGradjVectorsB( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type){ +int TriaElementCreateGradjVectorsB( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ int i,j; @@ -2770,5 +2770,5 @@ #undef __FUNCT__ #define __FUNCT__ "TriaElementMisfit" -int TriaElementMisfit( double* pJelem,void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type){ +int TriaElementMisfit( double* pJelem,void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type,int sub_analysis_type){ int i,j; Index: /issm/trunk/src/c/parallel/CreateFemModel.cpp =================================================================== --- /issm/trunk/src/c/parallel/CreateFemModel.cpp (revision 464) +++ /issm/trunk/src/c/parallel/CreateFemModel.cpp (revision 465) @@ -10,5 +10,5 @@ #include "../issm.h" -void CreateFemModel(FemModel* femmodel,ConstDataHandle MODEL,char* analysis_type){ +void CreateFemModel(FemModel* femmodel,ConstDataHandle MODEL,char* analysis_type,char* sub_analysis_type){ /*Model output: */ @@ -37,16 +37,8 @@ _printf_(" specifying analysis\n"); model->analysis_type=(char*)xmalloc((strlen(analysis_type)+1)*sizeof(char)); strcpy(model->analysis_type,analysis_type); + model->sub_analysis_type=(char*)xmalloc((strlen(sub_analysis_type)+1)*sizeof(char)); strcpy(model->sub_analysis_type,sub_analysis_type); - _printf_(" create elements, nodes and materials:\n"); - CreateElementsNodesAndMaterials(&elements,&nodes,&materials,model,MODEL); - - _printf_(" create constraints: \n"); - CreateConstraints(&constraints,model,MODEL); - - _printf_(" create loads: \n"); - CreateLoads(&loads,model,MODEL); - - _printf_(" create parameters: \n"); - CreateParameters(¶meters,model,MODEL); + _printf_(" create datasets:\n"); + CreateDataSets(&elements,&nodes,&materials,&constraints,&loads,¶meters,model,MODEL); _printf_(" create degrees of freedom: \n"); Index: /issm/trunk/src/c/parallel/GradJCompute.cpp =================================================================== --- /issm/trunk/src/c/parallel/GradJCompute.cpp (revision 464) +++ /issm/trunk/src/c/parallel/GradJCompute.cpp (revision 465) @@ -20,4 +20,5 @@ /*intermediary: */ int analysis_type; + int sub_analysis_type; char* solverstring=NULL; char* control_type=NULL; @@ -42,13 +43,14 @@ /*some parameters:*/ femmodel->parameters->FindParam((void*)&analysis_type,"analysis_type"); + femmodel->parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type"); femmodel->parameters->FindParam((void*)&solverstring,"solverstring"); femmodel->parameters->FindParam((void*)&control_type,"control_type"); _printf_("%s\n"," recover solution for this stiffness and right hand side:"); - diagnostic_core_nonlinear(&u_g,&K_ff0,&K_fs0,femmodel,inputs,analysis_type); + diagnostic_core_nonlinear(&u_g,&K_ff0,&K_fs0,femmodel,inputs,analysis_type,sub_analysis_type); VecToMPISerial(&u_g_double,u_g); VecFree(&u_g); _printf_("%s\n"," buid Du, difference between observed velocity and model velocity:"); - Dux( &du_g, femmodel->elements,femmodel->nodes,femmodel->loads,femmodel->materials,u_g_double,u_g_obs, inputs,analysis_type); + Dux( &du_g, femmodel->elements,femmodel->nodes,femmodel->loads,femmodel->materials,u_g_double,u_g_obs, inputs,analysis_type,sub_analysis_type); _printf_("%s\n"," reduce adjoint load from g-set to f-set:"); @@ -68,5 +70,5 @@ Gradjx( &grad_g, femmodel->elements,femmodel->nodes, femmodel->loads, femmodel->materials, - u_g_double,lambda_g_double, inputs,analysis_type,control_type); + u_g_double,lambda_g_double, inputs,analysis_type,sub_analysis_type,control_type); /*Free ressources:*/ Index: /issm/trunk/src/c/parallel/control.cpp =================================================================== --- /issm/trunk/src/c/parallel/control.cpp (revision 464) +++ /issm/trunk/src/c/parallel/control.cpp (revision 465) @@ -26,4 +26,5 @@ char* lockname=NULL; int analysis_type; + int sub_analysis_type; /*Intermediary: */ @@ -78,5 +79,5 @@ _printf_("read and create finite element model:\n"); - CreateFemModel(&femmodel,fid,"control"); + CreateFemModel(&femmodel,fid,"control",""); /*Recover parameters used throughout the solution:*/ @@ -93,4 +94,5 @@ femmodel.parameters->FindParam((void*)&u_g_obs,"u_g_obs"); femmodel.parameters->FindParam((void*)&analysis_type,"analysis_type"); + femmodel.parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type"); gsize=femmodel.nodes->NumberOfDofs(); @@ -156,5 +158,5 @@ inputs->Add(control_type,p_g,2,numberofnodes); inputs->Add("fit",fit[n]); - diagnostic_core_nonlinear(&u_g,NULL,NULL,&femmodel,inputs,analysis_type); + diagnostic_core_nonlinear(&u_g,NULL,NULL,&femmodel,inputs,analysis_type,sub_analysis_type); OutputControl(u_g,p_g,J,nsteps,femmodel.partition,outputfilename,femmodel.nodesets); _printf_("%s\n"," done."); @@ -175,5 +177,5 @@ UpdateFromInputsx(femmodel.elements,femmodel.nodes,femmodel.loads, femmodel.materials,inputs); - diagnostic_core_nonlinear(&u_g,NULL,NULL,&femmodel,inputs,analysis_type); + diagnostic_core_nonlinear(&u_g,NULL,NULL,&femmodel,inputs,analysis_type,sub_analysis_type); _printf_("%s\n"," saving final results..."); Index: /issm/trunk/src/c/parallel/diagnostic.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic.cpp (revision 464) +++ /issm/trunk/src/c/parallel/diagnostic.cpp (revision 465) @@ -59,13 +59,13 @@ _printf_("read and create finite element model:\n"); _printf_("\n reading diagnostic horiz model data:\n"); - CreateFemModel(&femmodels[0],fid,"diagnostic_horiz"); + CreateFemModel(&femmodels[0],fid,"diagnostic","horiz"); _printf_("\n reading diagnostic vert model data:\n"); - CreateFemModel(&femmodels[1],fid,"diagnostic_vert"); + CreateFemModel(&femmodels[1],fid,"diagnostic","vert"); _printf_("\n reading diagnostic stokes model data:\n"); - CreateFemModel(&femmodels[2],fid,"diagnostic_stokes"); + CreateFemModel(&femmodels[2],fid,"diagnostic","stokes"); _printf_("\n reading diagnostic hutter model data:\n"); - CreateFemModel(&femmodels[3],fid,"diagnostic_hutter"); + CreateFemModel(&femmodels[3],fid,"diagnostic","hutter"); _printf_("\n reading surface and bed slope computation model data:\n"); - CreateFemModel(&femmodels[4],fid,"slope_compute"); + CreateFemModel(&femmodels[4],fid,"slope_compute",""); _printf_("initialize inputs:\n"); Index: /issm/trunk/src/c/parallel/diagnostic_core.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic_core.cpp (revision 464) +++ /issm/trunk/src/c/parallel/diagnostic_core.cpp (revision 465) @@ -74,6 +74,6 @@ if(debug)_printf_("%s\n","computing surface slope (x and y derivatives)..."); - diagnostic_core_linear(&slopex,fem_sl,inputs,SurfaceSlopeComputeXAnalysisEnum()); - diagnostic_core_linear(&slopey,fem_sl,inputs,SurfaceSlopeComputeYAnalysisEnum()); + diagnostic_core_linear(&slopex,fem_sl,inputs,SlopeComputeAnalysisEnum(),SurfaceXAnalysisEnum()); + diagnostic_core_linear(&slopey,fem_sl,inputs,SlopeComputeAnalysisEnum(),SurfaceYAnalysisEnum()); if (dim==3){ @@ -89,5 +89,5 @@ if(debug)_printf_("%s\n"," computing hutter velocities..."); - diagnostic_core_linear(&ug,fem_dhu,inputs,DiagnosticHutterAnalysisEnum()); + diagnostic_core_linear(&ug,fem_dhu,inputs,DiagnosticAnalysisEnum(),HutterAnalysisEnum()); if(debug)_printf_("%s\n"," computing pressure according to MacAyeal..."); @@ -106,5 +106,5 @@ if(debug)_printf_("%s\n"," computing horizontal velocities..."); - diagnostic_core_nonlinear(&ug,NULL,NULL,fem_dh,inputs,DiagnosticHorizAnalysisEnum()); + diagnostic_core_nonlinear(&ug,NULL,NULL,fem_dh,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum()); if(debug)_printf_("%s\n"," computing pressure according to MacAyeal..."); @@ -121,5 +121,5 @@ if(debug)_printf_("%s\n"," computing vertical velocities..."); inputs->Add("velocity",ug_horiz,numberofdofspernode_dh,numberofnodes); - diagnostic_core_linear(&ug_vert,fem_dv,inputs,DiagnosticVertAnalysisEnum()); + diagnostic_core_linear(&ug_vert,fem_dv,inputs,DiagnosticAnalysisEnum(),VertAnalysisEnum()); if(debug)_printf_("%s\n"," combining horizontal and vertical velocities..."); @@ -138,6 +138,6 @@ if(debug)_printf_("%s\n","computing bed slope (x and y derivatives)..."); - diagnostic_core_linear(&slopex,fem_sl,inputs,BedSlopeComputeXAnalysisEnum()); - diagnostic_core_linear(&slopey,fem_sl,inputs,BedSlopeComputeYAnalysisEnum()); + diagnostic_core_linear(&slopex,fem_sl,inputs,SlopeComputeAnalysisEnum(),BedXAnalysisEnum()); + diagnostic_core_linear(&slopey,fem_sl,inputs,SlopeComputeAnalysisEnum(),BedYAnalysisEnum()); SlopeExtrudex( slopex, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials); SlopeExtrudex( slopey, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials); @@ -160,5 +160,5 @@ if(debug)_printf_("%s\n"," computing stokes velocities and pressure ..."); VecFree(&ug); - diagnostic_core_nonlinear(&ug,NULL,NULL,fem_ds,inputs,DiagnosticStokesAnalysisEnum()); + diagnostic_core_nonlinear(&ug,NULL,NULL,fem_ds,inputs,DiagnosticAnalysisEnum(),StokesAnalysisEnum()); //decondition" pressure Index: /issm/trunk/src/c/parallel/diagnostic_core_linear.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic_core_linear.cpp (revision 464) +++ /issm/trunk/src/c/parallel/diagnostic_core_linear.cpp (revision 465) @@ -11,5 +11,5 @@ #include "../issm.h" -void diagnostic_core_linear(Vec* pug,FemModel* fem,ParameterInputs* inputs,int analysis_type){ +void diagnostic_core_linear(Vec* pug,FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ /*parameters:*/ @@ -41,5 +41,5 @@ //*Generate system matrices if (debug) _printf_(" Generating matrices\n"); - SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type); + SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type,sub_analysis_type); /*!Reduce matrix from g to f size:*/ Index: /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp (revision 464) +++ /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp (revision 465) @@ -11,5 +11,5 @@ #include "../issm.h" -void diagnostic_core_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, FemModel* fem,ParameterInputs* inputs,int analysis_type){ +void diagnostic_core_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ @@ -84,9 +84,9 @@ if (debug) _printf_(" Generating matrices\n"); //*Generate system matrices - SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type); + SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type,sub_analysis_type); if (debug) _printf_(" Generating penalty matrices\n"); //*Generate penalty system matrices - PenaltySystemMatricesx(Kgg, pg,fem->elements,fem->nodes,loads,fem->materials,kflag,pflag,inputs,analysis_type); + PenaltySystemMatricesx(Kgg, pg,fem->elements,fem->nodes,loads,fem->materials,kflag,pflag,inputs,analysis_type,sub_analysis_type); if (debug) _printf_(" reducing matrix from g to f set\n"); @@ -122,5 +122,5 @@ inputs->Add("velocity",ug,numberofdofspernode,numberofnodes); - PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,loads,fem->materials,inputs,analysis_type); + PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,loads,fem->materials,inputs,analysis_type,sub_analysis_type); //Figure out if convergence is reached. @@ -167,5 +167,5 @@ kflag=1; pflag=0; //stiffness generation only - SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type); + SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type,sub_analysis_type); Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets); Index: /issm/trunk/src/c/parallel/objectivefunctionC.cpp =================================================================== --- /issm/trunk/src/c/parallel/objectivefunctionC.cpp (revision 464) +++ /issm/trunk/src/c/parallel/objectivefunctionC.cpp (revision 465) @@ -39,4 +39,5 @@ double* p_g_copy=NULL; int analysis_type; + int sub_analysis_type; Vec u_g=NULL; double* u_g_double=NULL; @@ -59,4 +60,5 @@ femmodel->parameters->FindParam((void*)&fit,"fit"); femmodel->parameters->FindParam((void*)&analysis_type,"analysis_type"); + femmodel->parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type"); femmodel->parameters->FindParam((void*)&numberofnodes,"numberofnodes"); @@ -75,5 +77,5 @@ //Run diagnostic with updated parameters. - diagnostic_core_nonlinear(&u_g,NULL,NULL,femmodel,inputs,analysis_type); + diagnostic_core_nonlinear(&u_g,NULL,NULL,femmodel,inputs,analysis_type,sub_analysis_type); VecToMPISerial(&u_g_double,u_g); VecFree(&u_g); @@ -81,5 +83,5 @@ inputs->Add("fit",fit[n]); Misfitx( &J, femmodel->elements,femmodel->nodes, femmodel->loads, femmodel->materials, - u_g_double,u_g_obs, inputs,analysis_type); + u_g_double,u_g_obs, inputs,analysis_type,sub_analysis_type); Index: /issm/trunk/src/c/parallel/parallel.h =================================================================== --- /issm/trunk/src/c/parallel/parallel.h (revision 464) +++ /issm/trunk/src/c/parallel/parallel.h (revision 465) @@ -12,6 +12,6 @@ void diagnostic_core(Vec* pug, Vec* ppg,FemModel* fems, ParameterInputs* inputs); -void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* fem,ParameterInputs* inputs,int analysis_type); -void diagnostic_core_linear(Vec* ppg,FemModel* fem,ParameterInputs* inputs,int analysis_type); +void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); +void diagnostic_core_linear(Vec* ppg,FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); int cielothermal_core(Vec** pt_g,ParameterInputs* inputs,FemModel* femmodel); @@ -29,9 +29,9 @@ //int ParameterUpdate(double* search_vector,int step, WorkspaceParams* workspaceparams,BatchParams* batchparams); void OutputDiagnostic(Vec u_g,Vec p_g, FemModel* femmodels,char* filename); -int OutputThermal(Vec* t_g,Vec* tpartition,char* filename,char* analysis_type); +int OutputThermal(Vec* t_g,Vec* tpartition,char* filename,char* analysis_type,int sub_analysis_type); void OutputControl(Vec u_g,double* p_g, double* J, int nsteps, Vec partition,char* filename,NodeSets* nodesets); void WriteLockFile(char* filename); -void CreateFemModel(FemModel* femmodel,ConstDataHandle MODEL,char* analysis_type); +void CreateFemModel(FemModel* femmodel,ConstDataHandle MODEL,char* analysis_type,char* sub_analysis_type); //int BatchDebug(Mat* Kgg,Vec* pg,FemModel* femmodel,char* filename); Index: /issm/trunk/src/c/parallel/thermal.cpp =================================================================== --- /issm/trunk/src/c/parallel/thermal.cpp (revision 465) +++ /issm/trunk/src/c/parallel/thermal.cpp (revision 465) @@ -0,0 +1,102 @@ +/*!\file: thermal.cpp + * \brief: thermal solution + */ + +#include "../issm.h" +#include "./parallel.h" + +#undef __FUNCT__ +#define __FUNCT__ "thermal" + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +int main(int argc,char* *argv){ + + /*I/O: */ + FILE* fid=NULL; + char* inputfilename=NULL; + char* outputfilename=NULL; + char* lockname=NULL; + int numberofnodes; + + /*Fem models : */ + FemModel femmodels[2]; + + /*initial velocity and pressure: */ + Vec u_g=NULL; + Vec p_g=NULL; + + + /*solution vectors: */ + Vec t_g=NULL; + Vec m_g=NULL; + + ParameterInputs* inputs=NULL; + int waitonlock=0; + + MODULEBOOT(); + + #if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_)) + throw ErrorException(__FUNCT__," parallel executable was compiled without support of parallel libraries!"); + #endif + + PetscInitialize(&argc,&argv,(char *)0,""); + + /*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=fopen(inputfilename,"rb"); + if(fid==NULL) throw ErrorException(__FUNCT__,exprintf("%s%s%s","could not open file ",inputfilename," for binary reading")); + + _printf_("read and create finite element model:\n"); + CreateFemModel(&femmodels[0],fid,"thermal",""); + CreateFemModel(&femmodels[1],fid,"melting",""); + + _printf_("initialize inputs:\n"); + femmodels[0].parameters->FindParam((void*)&u_g,"u_g"); + femmodels[0].parameters->FindParam((void*)&p_g,"p_g"); + femmodels[0].parameters->FindParam((void*)&numberofnodes,"numberofnodes"); + + inputs=new ParameterInputs; + //inputs->Add("velocity",u_g_initial,3,numberofnodes); + + //erase velocities: + //param=(Param*)femmodels[0].parameters->FindParamObject("u_g"); + //femmodels[0].parameters->DeleteObject((Object*)param); + + _printf_("call computational core:\n"); + diagnostic_core(&u_g,&p_g,femmodels,inputs); + + _printf_("write results to disk:\n"); + OutputDiagnostic(u_g,p_g,&femmodels[0],outputfilename); + + _printf_("write lock file:\n"); + femmodels[0].parameters->FindParam((void*)&waitonlock,"waitonlock"); + if (waitonlock){ + WriteLockFile(lockname); + } + + _printf_("closing MPI and Petsc\n"); + MPI_Barrier(MPI_COMM_WORLD); + + /*Close MPI libraries: */ + PetscFinalize(); + + + /*end module: */ + MODULEEND(); + + return 0; //unix success return; +} + Index: sm/trunk/src/c/parallel/thermalsteady.c =================================================================== --- /issm/trunk/src/c/parallel/thermalsteady.c (revision 464) +++ (revision ) @@ -1,111 +1,0 @@ -/* - * cielothermalsteady.c: - */ - -#include "../include/cielo.h" -#include "../modules.h" -#include "./parallel.h" - -#undef __FUNCT__ -#define __FUNCT__ "cielothermalsteady" - -int main(int argc,char* *argv){ - - /*I/O: */ - FILE* fid=NULL; - char* inputfilename=NULL; - char* outputfilename=NULL; - char* lockname=NULL; - char* analysis_type_t="thermalsteady"; - char* analysis_type_m="melting"; - - /*Intermediary: */ - - #if defined(_PARALLEL_) && defined(_HAVE_PETSC_) - FemModel femmodel_t; - FemModel femmodel_m; Vec* t_g=NULL; - #endif - ParameterInputs* inputs=NULL; - - #if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_)) - _printf_("%s%s\n",__FUNCT__," error message: parallel executable was compiled without support of parallel libraries!"); - return 1; - #else - - /*Initialize MPI environment: */ - PetscInitialize(&argc,&argv,(char *)0,""); - - /*Size and rank: */ - MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); - MPI_Comm_size(MPI_COMM_WORLD,&num_procs); - - /*Some checks on size of cluster*/ - if (num_procs<=1){ - _printf_("\nSize of MPI COMM WORLD is 1, needs to be at least 2. Include more nodes\n"); - PetscFinalize(); - return 0; - } - - /*Recover dbdir, input file name and output file name: */ - dbdir=argv[1]; - inputfilename=argv[2]; - outputfilename=argv[3]; - lockname=argv[4]; - - - - /*Read and create thermal finite element model: */ - fid=fopen(inputfilename,"rb"); - if(fid==NULL){ - _printf_("%s%s\n",__FUNCT__," error message: could not open file ",inputfilename," for binary reading"); - return 0; - } - if(!CreateFemModel(&femmodel_t,fid,analysis_type_t)){ - _printf_("%s\n",__FUNCT__," error message: could not read melting finite element model!\n"); - return 0; - } - fclose(fid); - - /*Read and create melting finite element model: */ - fid=fopen(inputfilename,"rb"); - if(!CreateFemModel(&femmodel_m,fid,analysis_type_m)){ - _printf_("%s\n",__FUNCT__," error message: could not read thermal finite element model!\n"); - return 0; - } - fclose(fid); - - /*Initialize inputs: */ - inputs=NewParameterInputs(); - - ParameterInputsAddFromMat(inputs,WorkspaceParamsGetVelocity(femmodel_t.workspaceparams),femmodel_t.workspaceparams->gsize,"velocity"); - ParameterInputsAddFromMat(inputs,WorkspaceParamsGetPressure(femmodel_t.workspaceparams),femmodel_t.workspaceparams->gsize,"pressure"); - - ParameterInputsAddFromMat(inputs,&femmodel_t.batchparams->dt,1,"dt"); - - /*Call computational core: */ - if(!cielothermal_core(&t_g,inputs,&femmodel_t)){ - _printf_("%s\n",__FUNCT__," error message: could not run computational core!\n"); - return 0; - } - - /*Write results to disk: */ - OutputThermal(t_g,femmodel_t.tpartition,outputfilename,analysis_type_t); - - /*Write lock file if requested: */ - if (femmodel_t.batchparams->waitonlock){ - WriteLockFile(lockname); - } - - cleanup_and_return: - _printf_("done.\n"); - - /*Synchronize everyone before exiting: */ - MPI_Barrier(MPI_COMM_WORLD); - - /*Close MPI libraries: */ - PetscFinalize(); - - return 0; //unix success return; - - #endif -} Index: /issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp =================================================================== --- /issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp (revision 465) +++ /issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp (revision 465) @@ -0,0 +1,47 @@ +/*!\file: DistributeNumDofs.cpp + * \brief: figure out the maximum number of dofs per grid. + */ + +#include "../../include/macros.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../shared.h" + +#undef __FUNCT__ +#define __FUNCT__ "DistributeNumDofs" + +int DistributeNumDofs(int *pnumdofs,int analysis_type,int sub_analysis_type){ + + int numdofs=2; //default numdofs + int i; + + /*ok, according to analysis type: */ + if (analysis_type==ControlAnalysisEnum()){ + numdofs=2; + } + else if (analysis_type==DiagnosticAnalysisEnum()){ + if (sub_analysis_type==HorizAnalysisEnum()){ + numdofs=2; + } + else if (sub_analysis_type==VertAnalysisEnum()){ + numdofs=1; + } + else if (sub_analysis_type==StokesAnalysisEnum()){ + numdofs=4; + } + else if (sub_analysis_type==HutterAnalysisEnum()){ + numdofs=2; + } + } + else if (analysis_type==SlopeComputeAnalysisEnum()){ + numdofs=1; + } + else if (analysis_type==ThermalAnalysisEnum()){ + numdofs=1; + } + else throw ErrorException(__FUNCT__," analysis type not implemented yet!"); + + /*Assign output pointers:*/ + *pnumdofs=numdofs;; + + return 1; +} Index: /issm/trunk/src/c/shared/Dofs/dofs.h =================================================================== --- /issm/trunk/src/c/shared/Dofs/dofs.h (revision 464) +++ /issm/trunk/src/c/shared/Dofs/dofs.h (revision 465) @@ -7,4 +7,5 @@ double* dofsetgen(int numdofs,int* doflist,int dofspernode,int totaldofs); +int DistributeNumDofs(int *pnumdofs,int analysis_type,int sub_analysis_type); Index: /issm/trunk/src/c/shared/Numerics/OptFunc.cpp =================================================================== --- /issm/trunk/src/c/shared/Numerics/OptFunc.cpp (revision 464) +++ /issm/trunk/src/c/shared/Numerics/OptFunc.cpp (revision 465) @@ -23,5 +23,5 @@ double J; - mxArray* inputs[8]; + mxArray* inputs[9]; mxArray* psearch_scalar=NULL; mxArray* mxJ=NULL; @@ -30,5 +30,4 @@ psearch_scalar=mxCreateDoubleScalar(scalar); inputs[0]=psearch_scalar; - inputs[1]=optargs->m; inputs[2]=optargs->inputs; @@ -38,6 +37,7 @@ inputs[6]=optargs->n; inputs[7]=optargs->analysis_type; + inputs[8]=optargs->sub_analysis_type; - mexCallMATLAB( 1, &mxJ, 8, (mxArray**)inputs, optargs->function_name); + mexCallMATLAB( 1, &mxJ, 9, (mxArray**)inputs, optargs->function_name); /*extract misfit from mxArray*/ Index: /issm/trunk/src/c/shared/Numerics/numerics.h =================================================================== --- /issm/trunk/src/c/shared/Numerics/numerics.h (revision 464) +++ /issm/trunk/src/c/shared/Numerics/numerics.h (revision 465) @@ -15,5 +15,4 @@ void cross(double* result,double* vector1,double* vector2); double norm(double* vector); -int DistributeNumDofs(int *pnumdofs,int analysis_type); #endif //ifndef _NUMERICS_H_ Index: /issm/trunk/src/m/classes/@model/model.m =================================================================== --- /issm/trunk/src/m/classes/@model/model.m (revision 464) +++ /issm/trunk/src/m/classes/@model/model.m (revision 465) @@ -255,4 +255,5 @@ md.solverstring=''; md.analysis_type=''; + md.sub_analysis_type=''; %management of large models Index: /issm/trunk/src/m/classes/public/BuildQueueingScript.m =================================================================== --- /issm/trunk/src/m/classes/public/BuildQueueingScript.m (revision 464) +++ /issm/trunk/src/m/classes/public/BuildQueueingScript.m (revision 465) @@ -1,7 +1,7 @@ -function BuildQueueingScript(md,solutiontype,executionpath,codepath) +function BuildQueueingScript(md,executionpath,codepath) %BUILDQUEUEINGSCRIPT - % % Usage: -% BuildQueueingScript(md,solutiontype,executionpath,codepath) +% BuildQueueingScript(md,executionpath,codepath) disp('building queueing script'); @@ -11,7 +11,7 @@ if exist(function_name,'file'), %Call this function: - eval([function_name '(md,solutiontype,executionpath,codepath);']); + eval([function_name '(md,executionpath,codepath);']); else %Call the generic BuildQueueingScript: - BuildQueueingScriptGeneric(md,solutiontype,executionpath,codepath); + BuildQueueingScriptGeneric(md,executionpath,codepath); end Index: /issm/trunk/src/m/classes/public/BuildQueueingScriptGeneric.m =================================================================== --- /issm/trunk/src/m/classes/public/BuildQueueingScriptGeneric.m (revision 464) +++ /issm/trunk/src/m/classes/public/BuildQueueingScriptGeneric.m (revision 465) @@ -1,3 +1,3 @@ -function BuildQueueingScriptGeneric(md,solutiontype,executionpath,codepath) +function BuildQueueingScriptGeneric(md,executionpath,codepath) %BUILDQUEUEINGSCRIPTGENERIC - ... % @@ -17,15 +17,14 @@ fprintf(fid,'mpirun -np %i ',md.np); -if strcmpi(solutiontype,'diagnostic_horiz') | strcmpi(solutiontype,'diagnostic'), +if strcmpi(md.analysis_type,'diagnostic'), fprintf(fid,'%s/diagnostic.exe',codepath); -elseif strcmpi(solutiontype,'control'), +elseif strcmpi(md.analysis_type,'control'), fprintf(fid,'%s/control.exe',codepath); -elseif strcmpi(solutiontype,'thermalsteady'), - fprintf(fid,'%s/thermalsteady.exe',codepath); +elseif strcmpi(md.analysis_type,'thermal'), + fprintf(fid,'%s/thermal.exe',codepath); else error('BuildQueueingScriptGeneric error message: unsupported solution type!'); end - fprintf(fid,' %s %s.bin %s.outbin %s.lock 2> %s.errlog >%s.outlog & ',executionpath,md.name,md.name,md.name,md.name,md.name); fclose(fid); Index: /issm/trunk/src/m/classes/public/BuildQueueingScriptcosmos.m =================================================================== --- /issm/trunk/src/m/classes/public/BuildQueueingScriptcosmos.m (revision 464) +++ /issm/trunk/src/m/classes/public/BuildQueueingScriptcosmos.m (revision 465) @@ -2,5 +2,5 @@ %BUILDQUEUEINGSCRIPTCOSMOS - build queueing script for batch system when running parallel % -% solutiontype is 'diagnostic','prognostic','transient','thermalsteady','thermaltransient','parameters' or 'control' +% solutiontype is 'diagnostic','prognostic','transient','thermal','parameters' or 'control' % and varargin is an optional package name ('cielo', 'ice' or 'macayeal') % @@ -38,12 +38,12 @@ fprintf(fid,'mpirun.lsf /home/larour/bin/alloc_cleanup.exe\n'); end -if strcmpi(solutiontype,'diagnostic_horiz') | strcmpi(solutiontype,'diagnostic'), +if strcmpi(solutiontype,'diagnostic'), fprintf(fid,'mpirun.lsf %s/diagnostic.exe %s %s.bin %s.outbin %s.lock',codepath,executionpath,md.name,md.name,md.name); elseif strcmpi(solutiontype,'control') , fprintf(fid,'mpirun.lsf %s/control.exe %s %s.bin %s.outbin %s.lock',codepath,executionpath,md.name,md.name,md.name); -elseif strcmpi(solutiontype,'thermalsteady') , - fprintf(fid,'mpirun.lsf %s/thermalsteady.exe %s %s.bin %s.outbin %s.lock',codepath,executionpath,md.name,md.name,md.name); +elseif strcmpi(solutiontype,'thermal') , + fprintf(fid,'mpirun.lsf %s/thermal.exe %s %s.bin %s.outbin %s.lock',codepath,executionpath,md.name,md.name,md.name); else - error('BuildQueueingScriptcosmso error message: unsupported solution type!'); + error('BuildQueueingScriptcosmos error message: unsupported solution type!'); end fclose(fid); Index: /issm/trunk/src/m/classes/public/BuildQueueingScriptgemini.m =================================================================== --- /issm/trunk/src/m/classes/public/BuildQueueingScriptgemini.m (revision 464) +++ /issm/trunk/src/m/classes/public/BuildQueueingScriptgemini.m (revision 465) @@ -22,10 +22,10 @@ fprintf(fid,'rm -rf %s.outlog %s.errlog %s.lock\n',md.name,md.name,md.name); -if strcmpi(solutiontype,'diagnostic_horiz') | strcmpi(solutiontype,'diagnostic'), +if strcmpi(solutiontype,'diagnostic') , fprintf(fid,'mpirun -np %i %s/diagnostic.exe %s %s.bin %s.outbin %s.lock',md.np,codepath,executionpath,md.name,md.name,md.name); elseif strcmpi(solutiontype,'control') , fprintf(fid,'mpirun -np %i %s/control.exe %s %s.bin %s.outbin %s.lock',md.np,codepath,executionpath,md.name,md.name,md.name); -elseif strcmpi(solutiontype,'thermalsteady') , - fprintf(fid,'mpirun -np %i %s/thermalsteady.exe %s %s.bin %s.outbin %s.lock',md.np,codepath,executionpath,md.name,md.name,md.name); +elseif strcmpi(solutiontype,'thermal') , + fprintf(fid,'mpirun -np %i %s/thermal.exe %s %s.bin %s.outbin %s.lock',md.np,codepath,executionpath,md.name,md.name,md.name); else error('BuildQueueingScriptgemini error message: unsupported solution type!'); Index: /issm/trunk/src/m/classes/public/ismodelselfconsistent.m =================================================================== --- /issm/trunk/src/m/classes/public/ismodelselfconsistent.m (revision 464) +++ /issm/trunk/src/m/classes/public/ismodelselfconsistent.m (revision 465) @@ -1,11 +1,11 @@ -function bool=ismodelselfconsistent(md,solutiontype,package), +function bool=ismodelselfconsistent(md,package), %ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem. % % Usage: -% bool=ismodelselfconsistent(md,solutiontype,package), +% bool=ismodelselfconsistent(md,package), %tolerance we use in litmus tests for the consistency of the model tolerance=10^-12; -if (nargin~=3 ) +if (nargin~=2 ) help ismodelselfconsistent error('ismodelselfconsistent error message: wrong usage'); @@ -115,5 +115,5 @@ %SIZE NUMBEROFGRIDS -fields={'x','y','z','B','drag','gridondirichlet_diag','melting','accumulation','surface','thickness','bed','gridonbed','gridonsurface','deadgrids'}; +fields={'x','y','z','B','drag','gridondirichlet_diag','melting','accumulation','surface','thickness','bed','gridonbed','gridonsurface'}; for i=1:length(fields), if find(length(eval(['md.' char(fields(i))]))~=md.numberofgrids), @@ -159,5 +159,5 @@ %DIAGNOSTIC -if strcmpi(solutiontype,'diagnostic') +if strcmpi(md.analysis_type,'diagnostic') %HUTTER ON ICESHELF WARNING @@ -189,5 +189,5 @@ %PROGNOSTIC -if strcmp(solutiontype,'prognostic'), +if strcmp(md.analysis_type,'prognostic'), %old testing % if isempty(find(md.gridondirichlet_diag)), @@ -208,23 +208,25 @@ %THERMAL TRANSIENT -if strcmp(solutiontype,'thermaltransient') - - %INITIAL TEMPERATURE, MELTING AND ACCUMULATION - if isempty(md.temperature), - disp(['An initial temperature is needed for a transient thermal computation']) - bool=0;return; - end - if isstruct(md.temperature) | isstruct(md.melting) | isstruct(md.accumulation), - disp(['The initial temperature, melting or accumulation should be a list and not a structure']) - bool=0;return; +if strcmp(md.analysis_type,'thermal') + if strcmp(md.sub_analysis_type,'transient') + + %INITIAL TEMPERATURE, MELTING AND ACCUMULATION + if isempty(md.temperature), + disp(['An initial temperature is needed for a transient thermal computation']) + bool=0;return; + end + if isstruct(md.temperature) | isstruct(md.melting) | isstruct(md.accumulation), + disp(['The initial temperature, melting or accumulation should be a list and not a structure']) + bool=0;return; + end end end %THERMAL STEADY AND THERMAL TRANSIENT -if strcmpi(solutiontype,'thermalsteady') | strcmp(solutiontype,'thermaltransient'), +if strcmpi(md.analysis_type,'thermal'), %EXTRUSION if strcmp(md.type,'2d'), - disp(['For a ' solutiontype ' computation, the model must be 3d, extrude it first!']) + disp(['For a ' md.analysis_type ' computation, the model must be 3d, extrude it first!']) bool=0;return; end @@ -242,5 +244,5 @@ %THERMAL TRANSIENT AND TRANSIENT -if strcmp(solutiontype,'thermaltransient') | strcmp(solutiontype,'transient'), +if strcmp(md.analysis_type,'thermal'), %DT and NDT @@ -261,5 +263,5 @@ %PARAMETERS -if strcmp(solutiontype,'parameters') +if strcmp(md.analysis_type,'parameters') %OUTPUT @@ -289,5 +291,5 @@ %CONTROL -if strcmpi(solutiontype,'control'), +if strcmpi(md.analysis_type,'control'), %CONTROL TYPE @@ -320,5 +322,5 @@ %QMU -if strcmpi(solutiontype,'qmu'), +if strcmpi(md.analysis_type,'qmu'), if ~strcmpi(md.cluster,'none'), if md.waitonlock==0, @@ -330,5 +332,5 @@ %MESH -if strcmpi(solutiontype,'mesh'), +if strcmpi(md.analysis_type,'mesh'), %this solution is a little special. It should come right after the md=model; operation. So a lot less checks! @@ -338,5 +340,5 @@ %MESH2GRID -if strcmpi(solutiontype,'mesh2grid'), +if strcmpi(md.analysis_type,'mesh2grid'), if ~strcmpi(md.cluster,'none'), disp(['model is not correctly configured: mesh2grid not supported in parallel yet!']); Index: /issm/trunk/src/m/classes/public/isresultconsistent.m =================================================================== --- /issm/trunk/src/m/classes/public/isresultconsistent.m (revision 464) +++ /issm/trunk/src/m/classes/public/isresultconsistent.m (revision 465) @@ -1,3 +1,3 @@ -function bool=isresultconsistent(md,solutiontype), +function bool=isresultconsistent(md) %ISRESULTCONSISTENT: check that results are consistent % @@ -6,15 +6,15 @@ % % Usage: -% bool=IsModelSelfConsistent(model,solutiontype) +% bool=IsModelSelfConsistent(model) %tolerance we use in litmus tests for the consistency of the model tolerance=10^-2; -if (nargin~=2 ) +if (nargin~=1 ) IsResultConsistentUsage(); error(' '); end %Check velocities i(no penalties for 2d meshes) -if (strcmp(solutiontype,'diagnostic') || strcmp(solutiontype,'diagnostic_horiz') ), +if strcmp(md.analysis_type,'diagnostic'), if strcmpi(md.type,'3d'), @@ -37,9 +37,9 @@ end -if strcmp(solutiontype,'thermalsteady') +if strcmp(md.analysis_type,'thermal') %check melting if any(abs(md.melting(md.numberofgrids2d+1:end))>tolerance) - disp(['''thermalsteady'' result not consistent: increase penalty_melting (negative melting)']); + disp(['''thermal'' result not consistent: increase penalty_melting (negative melting)']); bool=0; return; end @@ -47,5 +47,5 @@ end -if strcmp(solutiontype,'thermaltransient') +if strcmp(md.analysis_type,'thermaltransient') for iter=1:length(md.thermaltransient_results) @@ -61,5 +61,5 @@ -if (strcmp(solutiontype,'transient') & strcmpi(md.type,'3d')), +if (strcmp(md.analysis_type,'transient') & strcmpi(md.type,'3d')), for iter=1:length(md.transient_results) @@ -99,4 +99,4 @@ function IsResultConsistentUsage() disp(' '); -disp(' IsResultConsistent usage: flag=IsResultConsistent(model,solutiontype)'); +disp(' IsResultConsistent usage: flag=IsResultConsistent(model)'); disp(' where model is an instance of the @model class, and flag a boolean'); Index: /issm/trunk/src/m/classes/public/loadresultsfromcluster.m =================================================================== --- /issm/trunk/src/m/classes/public/loadresultsfromcluster.m (revision 464) +++ /issm/trunk/src/m/classes/public/loadresultsfromcluster.m (revision 465) @@ -1,7 +1,7 @@ -function md=loadresultsfromcluster(md,solutiontype) +function md=loadresultsfromcluster(md) %LOADRESULTSFROMCLUSTER - load results of solution sequence from cluster % % Usage: -% md=loadresultsfromcluster(md,solutiontype); +% md=loadresultsfromcluster(md); %Get cluster.rc location @@ -31,5 +31,5 @@ %If we are here, no errors in the solution sequence, call loadresultsfromdisk. -md=loadresultsfromdisk(md,[md.name '.outbin'],solutiontype); +md=loadresultsfromdisk(md,[md.name '.outbin'],md.analysis_type); %erase the log and output files Index: /issm/trunk/src/m/classes/public/loadresultsfromdisk.m =================================================================== --- /issm/trunk/src/m/classes/public/loadresultsfromdisk.m (revision 464) +++ /issm/trunk/src/m/classes/public/loadresultsfromdisk.m (revision 465) @@ -1,3 +1,3 @@ -function md=loadresultsfromdisk(md,filename,solutiontype) +function md=loadresultsfromdisk(md,filename) %LOADRESULTSFROMDISK - load results of solution sequence from disk file "filename" % @@ -103,5 +103,5 @@ disp(sprintf('%s\n','checking result consistency')); -if ~isresultconsistent(md,solutiontype), +if ~isresultconsistent(md), disp('!! results not consistent correct the model !!') %it would be very cruel to put an error, it would kill the computed results (even if not consistent...) end Index: /issm/trunk/src/m/classes/public/marshall.m =================================================================== --- /issm/trunk/src/m/classes/public/marshall.m (revision 464) +++ /issm/trunk/src/m/classes/public/marshall.m (revision 465) @@ -1,11 +1,10 @@ -function marshall(md,solutiontype,package) +function marshall(md) %MARSHALL - output Cielo compatible binary file from @model md, for certain solution type. % % The routine creates a Cielo compatible binary file from @model md -% solutiontype can be 'diagnostic','prognostic' or 'control'. % This binary file will be used for parallel runs in cielo. % % Usage: -% marshall(md,solutiontype,package) +% marshall(md) %some checks on list of arguments @@ -23,9 +22,6 @@ end -if strcmp(solutiontype,'diagnostic'), - WriteData(fid,'diagnostic_horiz','String','analysis_type'); -else - WriteData(fid,solutiontype,'String','analysis_type'); -end +WriteData(fid,md.analysis_type','String','analysis_type'); +WriteData(fid,md.sub_analysis_type','String','sub_analysis_type'); WriteData(fid,md.type,'String','type'); WriteData(fid,md.numberofgrids,'Integer','numberofgrids'); Index: /issm/trunk/src/m/classes/public/process_solve_options.m =================================================================== --- /issm/trunk/src/m/classes/public/process_solve_options.m (revision 465) +++ /issm/trunk/src/m/classes/public/process_solve_options.m (revision 465) @@ -0,0 +1,79 @@ +function outoptions=process_solve_options(options) +%DEFAULT_SOLVE_OPTIONS - set up default options for solve phase +% +% Usage: +% options=process_solve_options(options) +% +% See also: SOLVE,RECOVER_SOLVE_OPTIONS + + + + +%analysis_type: check on this option, error out otherwise +found=0; +for i=1:size(options,1), + if strcmpi(options{i,1},'analysis_type'), + analysis_type=options{i,2}; + found=1; + end +end +if ~found, + error('recover_solve_options error message: no ''analysis_type'' was provided'); +end + +%package: is there one? check to ''cielo'' +found=0; +for i=1:size(options,1), + if strcmpi(options{i,1},'package'), + package=options{i,2}; + found=1; + end +end +if ~found, + disp('recover_solve_options info message: no ''package'' was provided, defaulting to ''cielo'''); + options(end+1,:)={'package' 'cielo'}; + package='cielo'; +end + +if ~ischar(package), + error(['process_solve_options error message: package ' package ' not supported yet']); +end + +%sub_analysis_type: check on it, not mandatory +found=0; +for i=1:size(options,1), + if strcmpi(options{i,1},'sub_analysis_type'), + sub_analysis_type=options{i,2}; + found=1; + end +end +if ~found, + sub_analysis_type=''; +end + +%check solution type is supported +if ~(strcmpi(analysis_type,'control') | ... + strcmpi(analysis_type,'diagnostic') | ... + strcmpi(analysis_type,'prognostic') | ... + strcmpi(analysis_type,'thermal') | ... + strcmpi(analysis_type,'mesh') | ... + strcmpi(analysis_type,'mesh2grid') | ... + strcmpi(analysis_type,'qmu') | ... + strcmpi(analysis_type,'transient') ), + error(['process_solve_options error message: analysis_type ' analysis_type ' not supported yet!']); +end + +if ~(strcmpi(sub_analysis_type,'') | ... + strcmpi(sub_analysis_type,'steady_state') | ... + strcmpi(sub_analysis_type,'horiz') | ... + strcmpi(sub_analysis_type,'vert') | ... + strcmpi(sub_analysis_type,'transient') ), + error(['process_solve_options error message: sub_analysis_type ' sub_analysis_type ' not supported yet!']); +end + + + +%setup final options structure +outoptions.analysis_type=analysis_type; +outoptions.package=package; +outoptions.sub_analysis_type=sub_analysis_type; Index: /issm/trunk/src/m/classes/public/recover_solve_options.m =================================================================== --- /issm/trunk/src/m/classes/public/recover_solve_options.m (revision 465) +++ /issm/trunk/src/m/classes/public/recover_solve_options.m (revision 465) @@ -0,0 +1,22 @@ +function options=recover_solve_options(md,varargin) +%RECOVER_SOLVE_OPTIONS - recover solution options +% +% Usage: +% options=recover_solve_options(md,varargin); +% +% See also: SOLVE + +%initialize options. +options=cell(0,2); + +%make sure length(varargin) is even, ie options come in pairs. +if mod(length(varargin),2), + error('recover_solve_options error message: an even number of options is necessary'); +end + +%go through varargin, extract options +for i=1:length(varargin)/2, + + options(end+1,:)={varargin{2*i-1} varargin{2*i}}; + +end Index: /issm/trunk/src/m/classes/public/solve.m =================================================================== --- /issm/trunk/src/m/classes/public/solve.m (revision 464) +++ /issm/trunk/src/m/classes/public/solve.m (revision 465) @@ -1,64 +1,47 @@ -function md=solve(md,solutiontype,varargin) +function md=solve(md,varargin) %SOLVE - apply solution sequence for this model % -% solutiontype is 'diagnostic','prognostic','transient','thermalsteady','thermaltransient','parameters' or 'control' -% and varargin is an optional package name ('cielo', 'ice' or 'macayeal') -% % Usage: -% md=solve(md,solutiontype,varargin) +% md=solve(md,varargin) +% where varargin is a lit of paired arguments. +% arguments can be: 'analysis_type': 'diagnostic','thermal','prognostic','transient' +% arguments can be: 'sub_analysis_type': 'transient','' +% arguments can be: 'package': 'macayeal','ice','cielo' % % Examples: -% md=solve(md,'diagnostic','macayeal'); -% md=solve(md,'control','cielo'); +% md=solve(md,'analysis_type','diagnostic','package','cielo'); +% md=solve(md,'analysis_type','control','package','ice'); +% md=solve(md,'analysis_type','thermal','sub_analysis_type','transient','package','ice'); +% md=solve(md,'analysis_type','thermal','sub_analysis_type','','package','cielo'); %some checks on list of arguments - global ISSM_DIR -solutions={'mesh2grid','mesh','thermaltransient','thermalsteady','qmu','diagnostic','diagnostic_horiz','prognostic','transient','parameters','control'}; -found=0; -for i=1:length(solutions), - if strcmpi(solutiontype,solutions{i}), - found=1; - break; - end -end +%recover options +options=recover_solve_options(md,varargin{:}); -if found==0, - error(['solve error message: solution type ' solutiontype ' not supported yet!']); -end +%add default options +options=process_solve_options(options); -%we are good with this solutiontype, put in the analysis_type field of md: -md.analysis_type=solutiontype; - -%Recover type of package being used: -if nargin==2, - package='Ice'; -else - package=varargin{1}; -end - -if ~ischar(package), - error('Package specified in varargin can only be ''ice'', or ''cielo'''); -end - -if ~(strcmpi(package,'ice') || strcmpi(package,'cielo') || strcmpi(package,'macayeal')) - error('Package specified in varargin can only be ''ice'', ''macayeal'', or ''cielo'''); -end +%recover some fields +md.analysis_type=options.analysis_type; +md.sub_analysis_type=options.sub_analysis_type; +package=options.package; %Use package to set solution namespace usenamespace(package); -if ~ismodelselfconsistent(md,solutiontype,package), - error(' '); %previous error messages should explain what is going on. +if ~ismodelselfconsistent(md,package), + error(' '); %previous error messages should explain what is going on. end displaystring(md.debug,'\n%s\n','launching solution sequence'); + %If running in parallel, we have a different way of launching the solution -%sequences. -if ~strcmpi(solutiontype,'qmu'), +%sequences. qmu is the only solution sequence that cannot run in parallel +if ~strcmpi(md.analysis_type,'qmu'), if ~strcmpi(md.cluster,'none'), - md=solveparallel(md,solutiontype,package); + md=solveparallel(md) return; end @@ -66,30 +49,27 @@ %Lauch correct solution sequence -if strcmpi(solutiontype,'diagnostic'), +if strcmpi(md.analysis_type,'diagnostic'), md=diagnostic(md); -elseif strcmpi(solutiontype,'mesh'), +elseif strcmpi(md.analysis_type,'mesh'), md=mesh(md); -elseif strcmpi(solutiontype,'transient'), +elseif strcmpi(md.analysis_type,'transient'), md=transient(md); -elseif strcmpi(solutiontype,'qmu'), +elseif strcmpi(md.analysis_type,'qmu'), md=qmu(md,package); -elseif strcmpi(solutiontype,'parameters'), - md=parameters(md);; - -elseif strcmpi(solutiontype,'mesh2grid'), +elseif strcmpi(md.analysis_type,'mesh2grid'), md=mesh2grid(md);; -elseif strcmpi(solutiontype,'prognostic'), +elseif strcmpi(md.analysis_type,'prognostic'), md=prognostic(md);; -elseif strcmpi(solutiontype,'control'), +elseif strcmpi(md.analysis_type,'control'), md=control(md); -elseif strcmpi(solutiontype,'thermalsteady') | strcmpi(solutiontype,'thermaltransient'), - md=thermal(md,solutiontype); +elseif strcmpi(md.analysis_type,'thermal'), + md=thermal(md); else error('solution type not supported for this model configuration.'); @@ -98,5 +78,5 @@ %Check result is consistent displaystring(md.debug,'%s\n','checking result consistency'); -if ~isresultconsistent(md,solutiontype), +if ~isresultconsistent(md), disp('!! results not consistent correct the model !!') %it would be very cruel to put an error, it would kill the computed results (even if not consistent...) end @@ -110,5 +90,5 @@ function solveusage(); disp(' '); -disp(' Solve usage: md=solve(md,solutiontype,package)'); -disp(' solutiontype can be ''diagnostic'',''transient'', ''thermalsteady'',''thermaltransient'',''parameters'',''mesh2grid'' or ''control'''); +disp(' Solve usage: md=solve(md,md.analysis_type,package)'); +disp(' md.analysis_type can be ''diagnostic'',''transient'', ''thermal'',''thermaltransient'',''parameters'',''mesh2grid'' or ''control'''); disp(' package is either ''cielo'' or ''ice'''); Index: /issm/trunk/src/m/classes/public/solveparallel.m =================================================================== --- /issm/trunk/src/m/classes/public/solveparallel.m (revision 464) +++ /issm/trunk/src/m/classes/public/solveparallel.m (revision 465) @@ -1,22 +1,7 @@ -function md=solveparallel(md,solutiontype,varargin) +function md=solveparallel(md) %SOLVEPARALLEL - solution sequence using a cluster in parallel mode % % Usage: -% md=solveparallel(md,solutiontype,varargin) - -%Recover type of package being used: -if nargin==2, - package='Ice'; -else - package=varargin{1}; -end - -if ~ischar(package), - error('Package specified in varargin can only be ''ice'', or ''cielo'''); -end - -if ~(strcmpi(package,'ice') || strcmpi(package,'cielo') || strcmpi(package,'macayeal')) - error('Package specified in varargin can only be ''ice'', ''macayeal'', or ''cielo'''); -end +% md=solveparallel(md); %Get cluster.rc location @@ -27,8 +12,8 @@ %Marshall model data into a binary file. -marshall(md,solutiontype,package); +marshall(md); %Now, we need to build the queuing script, used by the cluster to launch the job. -BuildQueueingScript(md,solutiontype,executionpath,codepath); +BuildQueueingScript(md,executionpath,codepath); %Now, launch the queueing script @@ -40,5 +25,5 @@ waitonlock([executionpath '/' md.name '.lock']); %load results - md=loadresultsfromcluster(md,solutiontype); + md=loadresultsfromcluster(md); else return; Index: /issm/trunk/src/m/solutions/cielo/GradJCompute.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/GradJCompute.m (revision 464) +++ /issm/trunk/src/m/solutions/cielo/GradJCompute.m (revision 465) @@ -1,3 +1,3 @@ -function [u_g grad_g]=GradJCompute(m,inputs, u_g_obs,analysis_type); +function [u_g grad_g]=GradJCompute(m,inputs, u_g_obs,analysis_type,sub_analysis_type); %Recover solution for this stiffness and right hand side: @@ -5,5 +5,5 @@ disp(' computing velocities...') end -[u_g K_ff0 K_fs0 ]=diagnostic_core_nonlinear(m,inputs,analysis_type); +[u_g K_ff0 K_fs0 ]=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type); %Buid Du, difference between observed velocity and model velocity. @@ -11,5 +11,5 @@ disp(' computing Du...') end -[Du_g]=Du(m.elements,m.nodes,m.loads,m.materials,m.parameters, u_g,u_g_obs,inputs); +[Du_g]=Du(m.elements,m.nodes,m.loads,m.materials,m.parameters, u_g,u_g_obs,inputs,analysis_type,sub_analysis_type); %Reduce adjoint load from g-set to f-set @@ -26,3 +26,3 @@ %Compute gradJ -grad_g=Gradj(m.elements,m.nodes,m.loads,m.materials,m.parameters, u_g, lambda_g,inputs); +grad_g=Gradj(m.elements,m.nodes,m.loads,m.materials,m.parameters, u_g, lambda_g,inputs,analysis_type,sub_analysis_type); Index: /issm/trunk/src/m/solutions/cielo/control.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/control.m (revision 464) +++ /issm/trunk/src/m/solutions/cielo/control.m (revision 465) @@ -6,4 +6,5 @@ %Build all models requested for control simulation md.analysis_type='control'; + md.sub_analysis_type=''; m=CreateFemModel(md); @@ -40,5 +41,5 @@ disp(' computing gradJ...'); - [u_g c(n).grad_g]=GradJCompute(m,inputs,u_g_obs,md.analysis_type); + [u_g c(n).grad_g]=GradJCompute(m,inputs,u_g_obs,md.analysis_type,md.sub_analysis_type); inputs=add(inputs,'velocity',u_g,'doublevec',2,m.parameters.numberofnodes); disp(' done.'); @@ -58,5 +59,5 @@ disp(' optimizing along gradient direction...'); - [search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,m,inputs,p_g,u_g_obs,c(n).grad_g,n,md.analysis_type); + [search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,m,inputs,p_g,u_g_obs,c(n).grad_g,n,md.analysis_type,md.sub_analysis_type); disp(' done.'); @@ -76,5 +77,5 @@ %some temporary saving if(mod(n,5)==0), - solution=controlfinalsol(c,m,p_g,inputs,md.analysis_type); + solution=controlfinalsol(c,m,p_g,inputs,md.analysis_type,md.sub_analysis_type); save temporary_control_results solution end @@ -85,5 +86,5 @@ %Create final solution disp(' preparing final velocity solution...'); - solution=controlfinalsol(c,m,p_g,inputs,md.analysis_type); + solution=controlfinalsol(c,m,p_g,inputs,md.analysis_type,md.sub_analysis_type); disp(' done.'); Index: /issm/trunk/src/m/solutions/cielo/controlfinalsol.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/controlfinalsol.m (revision 464) +++ /issm/trunk/src/m/solutions/cielo/controlfinalsol.m (revision 465) @@ -1,7 +1,7 @@ -function solution=controlfinalsol(c,m,p_g,inputs,analysis_type); +function solution=controlfinalsol(c,m,p_g,inputs,analysis_type,sub_analysis_type); %From parameters, build inputs for icediagnostic_core, using the final parameters inputs=add(inputs,m.parameters.control_type,p_g,'doublevec',2,m.parameters.numberofnodes); -u_g=diagnostic_core_nonlinear(m,inputs,analysis_type); +u_g=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type); %Build partitioning vectors to recover solution Index: /issm/trunk/src/m/solutions/cielo/diagnostic.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/diagnostic.m (revision 464) +++ /issm/trunk/src/m/solutions/cielo/diagnostic.m (revision 465) @@ -10,17 +10,17 @@ %Build all models requested for diagnostic simulation displaystring(md.debug,'%s',['reading diagnostic horiz model data']); - md.analysis_type='diagnostic_horiz'; m_dh=CreateFemModel(md); + md.analysis_type='diagnostic'; md.sub_analysis_type='horiz'; m_dh=CreateFemModel(md); displaystring(md.debug,'\n%s',['reading diagnostic vert model data']); - md.analysis_type='diagnostic_vert'; m_dv=CreateFemModel(md); + md.analysis_type='diagnostic'; md.sub_analysis_type='vert'; m_dv=CreateFemModel(md); displaystring(md.debug,'\n%s',['reading diagnostic stokes model data']); - md.analysis_type='diagnostic_stokes'; m_ds=CreateFemModel(md); + md.analysis_type='diagnostic'; md.sub_analysis_type='stokes'; m_ds=CreateFemModel(md); displaystring(md.debug,'\n%s',['reading diagnostic hutter model data']); - md.analysis_type='diagnostic_hutter'; m_dhu=CreateFemModel(md); + md.analysis_type='diagnostic'; md.sub_analysis_type='hutter'; m_dhu=CreateFemModel(md); displaystring(md.debug,'\n%s',['reading surface and bed slope computation model data']); - md.analysis_type='slope_compute'; m_sl=CreateFemModel(md); + md.analysis_type='slope_compute'; md.sub_analysis_type=''; m_sl=CreateFemModel(md); % figure out number of dof: just for information purposes. Index: /issm/trunk/src/m/solutions/cielo/diagnostic_core.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/diagnostic_core.m (revision 464) +++ /issm/trunk/src/m/solutions/cielo/diagnostic_core.m (revision 465) @@ -16,6 +16,6 @@ displaystring(debug,'\n%s',['computing surface slope (x and y derivatives)...']); - slopex=diagnostic_core_linear(m_sl,inputs,'surface_slope_compute_x'); - slopey=diagnostic_core_linear(m_sl,inputs,'surface_slope_compute_y'); + slopex=diagnostic_core_linear(m_sl,inputs,'slope_compute','surfacex'); + slopey=diagnostic_core_linear(m_sl,inputs,'slope_compute','surfacey'); if dim==3, @@ -30,5 +30,5 @@ displaystring(debug,'\n%s',['computing hutter velocities...']); - u_g=diagnostic_core_linear(m_dhu,inputs,'diagnostic_hutter'); + u_g=diagnostic_core_linear(m_dhu,inputs,'diagnostic','hutter'); displaystring(debug,'\n%s',['computing pressure according to MacAyeal...']); @@ -46,5 +46,5 @@ displaystring(debug,'\n%s',['computing horizontal velocities...']); - u_g=diagnostic_core_nonlinear(m_dh,inputs,'diagnostic_horiz'); + u_g=diagnostic_core_nonlinear(m_dh,inputs,'diagnostic','horiz'); displaystring(debug,'\n%s',['computing pressure according to MacAyeal...']); @@ -59,5 +59,5 @@ displaystring(debug,'\n%s',['computing vertical velocities...']); inputs=add(inputs,'velocity',u_g_horiz,'doublevec',m_dh.parameters.numberofdofspernode,m_dh.parameters.numberofnodes); - u_g_vert=diagnostic_core_linear(m_dv,inputs,'diagnostic_vert'); + u_g_vert=diagnostic_core_linear(m_dv,inputs,'diagnostic','vert'); displaystring(debug,'\n%s',['combining horizontal and vertical velocities...']); @@ -75,6 +75,6 @@ displaystring(debug,'\n%s',['computing bed slope (x and y derivatives)...']); - slopex=diagnostic_core_linear(m_sl,inputs,'bed_slope_compute_x'); - slopey=diagnostic_core_linear(m_sl,inputs,'bed_slope_compute_y'); + slopex=diagnostic_core_linear(m_sl,inputs,'slope_compute','bedx'); + slopey=diagnostic_core_linear(m_sl,inputs,'slope_compute','bedy'); slopex=SlopeExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,slopex); slopey=SlopeExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,slopey); @@ -95,5 +95,5 @@ displaystring(debug,'\n%s',['computing stokes velocities and pressure ...']); - u_g=diagnostic_core_nonlinear(m_ds,inputs,'diagnostic_stokes'); + u_g=diagnostic_core_nonlinear(m_ds,inputs,'diagnostic','stokes'); %"decondition" pressure Index: /issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m (revision 464) +++ /issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m (revision 465) @@ -1,3 +1,3 @@ -function [u_g varargout]=diagnostic_core_nonlinear(m,inputs,analysis_type) +function [u_g varargout]=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type) %INPUT function [ru_g varargout]=cielodiagnostic_core_nonlinear(m,inputs) @@ -30,8 +30,8 @@ %system matrices - [K_gg_nopenalty , p_g_nopenalty]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type); + [K_gg_nopenalty , p_g_nopenalty]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); %penalties - [K_gg , p_g]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type); + [K_gg , p_g]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); @@ -60,5 +60,5 @@ %penalty constraints - [loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints( m.elements,m.nodes, loads, m.materials,m.parameters,inputs); + [loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints( m.elements,m.nodes, loads, m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); % Figure out if convergence is reached. @@ -107,5 +107,5 @@ inputs=add(inputs,'velocity',soln(count).u_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes); m.parameters.kflag=1; m.parameters.pflag=0; - [K_gg, p_g]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type); + [K_gg, p_g]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); [K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets); varargout(1)={K_ff}; Index: /issm/trunk/src/m/solutions/cielo/objectivefunctionC.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/objectivefunctionC.m (revision 464) +++ /issm/trunk/src/m/solutions/cielo/objectivefunctionC.m (revision 465) @@ -1,3 +1,3 @@ -function J =objectivefunctionC(search_scalar,m,inputs,p_g,u_g_obs,grad_g,n,analysis_type); +function J =objectivefunctionC(search_scalar,m,inputs,p_g,u_g_obs,grad_g,n,analysis_type,sub_analysis_type); %recover some parameters @@ -13,6 +13,6 @@ %Run diagnostic with updated parameters. -u_g=diagnostic_core_nonlinear(m,inputs,analysis_type); +u_g=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type); %Compute misfit for this velocity field. -J=Misfit(m.elements,m.nodes,m.loads,m.materials,m.parameters, u_g, u_g_obs,inputs); +J=Misfit(m.elements,m.nodes,m.loads,m.materials,m.parameters, u_g, u_g_obs,inputs,analysis_type,sub_analysis_type); Index: /issm/trunk/src/m/solutions/cielo/thermal.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/thermal.m (revision 464) +++ /issm/trunk/src/m/solutions/cielo/thermal.m (revision 465) @@ -7,33 +7,32 @@ %timing t1=clock; - - %md.analysis_type is set to thermalsteady or thermaltransient. save it. - solutiontype=md.analysis_type; - %Build all models requested for thermal and melting simulation - md.analysis_type=solutiontype; m_t=CreateFemModel(md); - md.analysis_type='melting' m_m=CreateFEMModel(md); + %Build all models requested for diagnostic simulation + displaystring(md.debug,'%s',['reading thermal model data']); + md.analysis_type='thermal'; m_t=CreateFemModel(md); + + displaystring(md.debug,'%s',['reading melting model data']); + md.analysis_type='melting'; m_m=CreateFemModel(md); % figure out number of dof: just for information purposes. - md.dof=m_t.nodesets.fsize; %biggest dof number + md.dof=m_t.nodesets.fsize; - - if strcmpi(solutiontype,'thermalsteady'), - - %initialize velocity and pressure - u_g=m_t.parameters.u_g; - p_g=m_t.parameters.p_g; - - inputs=struct('pressure',p_g,'velocity',u_g); - - disp(sprintf('\n%s\n','computing temperature...')); - [t_g,m_t]=thermal_core(m_t,inputs); - - disp(sprintf('\n%s\n','computing melting...')); - inputs=struct('pressure',p_g,'temperature',t_g); + %initialize inputs + displaystring(debug,'\n%s',['setup inputs...']); + inputs=inputlist; + inputs=add(inputs,'velocity',m_t.parameters.u_g,'doublevec',3,m_t.parameters.numberofnodes); + inputs=add(inputs,'pressure',m_t.parameters.p_g,'doublevec',1,m_t.parameters.numberofnodes); + inputs=add(inputs,'dt',m_t.parameters.dt,'double'); + + if strcmpi(md.sub_analysis_type,'steady'), + + displaystring(debug,'\n%s',['computing temperatures...']); + [t_g loads_t melting_offset]=thermal_core(m_t,inputs,'thermal','steady'); - melting_g=melting_core(m_m,inputs); - - %load results onto model + displaystring(debug,'\n%s',['computing melting...']); + inputs=add(inputs,'melting_offset',melting_offset,'double'); + melting_g=melting_core(m_m,inputs,'melting','steady'); + + displaystring(debug,'\n%s',['load results...']); md.temperature=t_g; md.melting=melting_g*md.yts; %from m/s to m/a @@ -41,9 +40,9 @@ else - %initialize velocity, pressure and temperature - u_g=m.parameters.u_g; - p_g=m.parameters.p_g; - t_g=m.parameters.t_g; - melting_g=m.parameters.melting_g; + %initialize velocity, pressure, temperature and melting + u_g=m_t.parameters.u_g; + p_g=m_t.parameters.p_g; + t_g=m_t.parameters.t_g; + melting_g=m_t.parameters.melting_g; nsteps=m_t.parameters.ndt/m_t.parameters.dt; @@ -55,9 +54,9 @@ for n=1:nsteps, - disp(sprintf('\n%s%i/%i\n','time step: ',n,md.ndt/md.dt)); - disp(' computing temperature...'); - - inputs=struct('pressure',p_g,'temperature',soln(n).t_g,'velocity',u_g,'dt',m_t.parameters.dt); - [soln(n+1).t_g,m_t]=thermal_core(m_t,inputs); + displaystring(debug,'\n%s%i/%i\n','time step: ',n,nsteps); + + displaystring(debug,'\n%s',[' computing temperatures...']); + inputs=add(inputs,'temperature',soln(n).t_g,'doublevec',1,m_t.parameters.numberofnodes); + [soln(n+1).t_g loads_t melting_offset]=thermal_core(m_t,inputs); disp(' computing melting...'); @@ -80,3 +79,6 @@ md.melting=solution_melting; end - + + %stop timing + t2=clock; + displaystring(md.debug,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']); Index: /issm/trunk/src/mex/ControlOptimization/ControlOptimization.cpp =================================================================== --- /issm/trunk/src/mex/ControlOptimization/ControlOptimization.cpp (revision 464) +++ /issm/trunk/src/mex/ControlOptimization/ControlOptimization.cpp (revision 465) @@ -51,4 +51,5 @@ optargs.n=STEP; optargs.analysis_type=ANALYSIS; + optargs.sub_analysis_type=SUBANALYSIS; optpars.xmin=xmin; @@ -73,5 +74,5 @@ { _printf_("\n"); - _printf_(" usage: [search_scalar J] = %s(function_name,xmin,xmax,options,m,inputs,p_g,u_g_obs,grad_g,step,analysis_type)\n",__FUNCT__); + _printf_(" usage: [search_scalar J] = %s(function_name,xmin,xmax,options,m,inputs,p_g,u_g_obs,grad_g,step,analysis_type,sub_analysis_type)\n",__FUNCT__); _printf_("\n"); } Index: /issm/trunk/src/mex/ControlOptimization/ControlOptimization.h =================================================================== --- /issm/trunk/src/mex/ControlOptimization/ControlOptimization.h (revision 464) +++ /issm/trunk/src/mex/ControlOptimization/ControlOptimization.h (revision 465) @@ -28,4 +28,5 @@ #define STEP (mxArray*)prhs[9] #define ANALYSIS (mxArray*)prhs[10] +#define SUBANALYSIS (mxArray*)prhs[11] /* serial output macros: */ @@ -37,5 +38,5 @@ #define NLHS 2 #undef NRHS -#define NRHS 11 +#define NRHS 12 Index: /issm/trunk/src/mex/Du/Du.cpp =================================================================== --- /issm/trunk/src/mex/Du/Du.cpp (revision 464) +++ /issm/trunk/src/mex/Du/Du.cpp (revision 465) @@ -15,8 +15,11 @@ DataSet* loads=NULL; DataSet* materials=NULL; - int analysis_type; double* u_g=NULL; double* u_g_obs=NULL; ParameterInputs* inputs=NULL; + char* analysis_type_string=NULL; + int analysis_type; + char* sub_analysis_type_string=NULL; + int sub_analysis_type; /* output datasets: */ @@ -35,7 +38,10 @@ FetchData((void**)&loads,NULL,NULL,LOADS,"DataSet",NULL); FetchData((void**)&materials,NULL,NULL,MATERIALS,"DataSet",NULL); - FetchData((void**)&analysis_type,NULL,NULL,mxGetField(PARAMETERS,0,"analysis_type"),"Integer",NULL); FetchData((void**)&u_g,NULL,NULL,UG,"Vector","Vec"); FetchData((void**)&u_g_obs,NULL,NULL,UGOBS,"Vector","Vec"); + FetchData((void**)&analysis_type_string,NULL,NULL,ANALYSIS,"String",NULL); + analysis_type=AnalysisTypeAsEnum(analysis_type_string); + FetchData((void**)&sub_analysis_type_string,NULL,NULL,SUBANALYSIS,"String",NULL); + sub_analysis_type=AnalysisTypeAsEnum(sub_analysis_type_string); /*Fetch inputs: */ @@ -44,5 +50,5 @@ /*!Call core code: */ - Dux(&du_g, elements,nodes,loads,materials,u_g,u_g_obs,inputs,analysis_type); + Dux(&du_g, elements,nodes,loads,materials,u_g,u_g_obs,inputs,analysis_type,sub_analysis_type); /*write output : */ @@ -58,4 +64,6 @@ VecFree(&du_g); delete inputs; + xfree((void**)&analysis_type_string); + xfree((void**)&sub_analysis_type_string); /*end module: */ Index: /issm/trunk/src/mex/Du/Du.h =================================================================== --- /issm/trunk/src/mex/Du/Du.h (revision 464) +++ /issm/trunk/src/mex/Du/Du.h (revision 465) @@ -25,4 +25,6 @@ #define UGOBS (mxArray*)prhs[6] #define INPUTS (mxArray*)prhs[7] +#define ANALYSIS (mxArray*)prhs[8] +#define SUBANALYSIS (mxArray*)prhs[9] /* serial output macros: */ @@ -33,5 +35,5 @@ #define NLHS 1 #undef NRHS -#define NRHS 8 +#define NRHS 10 Index: /issm/trunk/src/mex/Gradj/Gradj.cpp =================================================================== --- /issm/trunk/src/mex/Gradj/Gradj.cpp (revision 464) +++ /issm/trunk/src/mex/Gradj/Gradj.cpp (revision 465) @@ -15,9 +15,12 @@ DataSet* loads=NULL; DataSet* materials=NULL; - int analysis_type; char* control_type=NULL; double* u_g=NULL; double* lambda_g=NULL; ParameterInputs* inputs=NULL; + char* analysis_type_string=NULL; + int analysis_type; + char* sub_analysis_type_string=NULL; + int sub_analysis_type; /* output datasets: */ @@ -36,8 +39,11 @@ FetchData((void**)&loads,NULL,NULL,LOADS,"DataSet",NULL); FetchData((void**)&materials,NULL,NULL,MATERIALS,"DataSet",NULL); - FetchData((void**)&analysis_type,NULL,NULL,mxGetField(PARAMETERS,0,"analysis_type"),"Integer",NULL); FetchData((void**)&control_type,NULL,NULL,mxGetField(PARAMETERS,0,"control_type"),"String",NULL); FetchData((void**)&u_g,NULL,NULL,UG,"Vector","Vec"); FetchData((void**)&lambda_g,NULL,NULL,LAMBDAG,"Vector","Vec"); + FetchData((void**)&analysis_type_string,NULL,NULL,ANALYSIS,"String",NULL); + analysis_type=AnalysisTypeAsEnum(analysis_type_string); + FetchData((void**)&sub_analysis_type_string,NULL,NULL,SUBANALYSIS,"String",NULL); + sub_analysis_type=AnalysisTypeAsEnum(sub_analysis_type_string); /*Fetch inputs: */ @@ -46,5 +52,5 @@ /*!Call core code: */ - Gradjx(&grad_g, elements,nodes,loads,materials,u_g,lambda_g,inputs,analysis_type,control_type); + Gradjx(&grad_g, elements,nodes,loads,materials,u_g,lambda_g,inputs,analysis_type,sub_analysis_type,control_type); /*write output : */ @@ -61,4 +67,6 @@ delete inputs; VecFree(&grad_g); + xfree((void**)&analysis_type_string); + xfree((void**)&sub_analysis_type_string); /*end module: */ Index: /issm/trunk/src/mex/Gradj/Gradj.h =================================================================== --- /issm/trunk/src/mex/Gradj/Gradj.h (revision 464) +++ /issm/trunk/src/mex/Gradj/Gradj.h (revision 465) @@ -25,4 +25,6 @@ #define LAMBDAG (mxArray*)prhs[6] #define INPUTS (mxArray*)prhs[7] +#define ANALYSIS (mxArray*)prhs[8] +#define SUBANALYSIS (mxArray*)prhs[9] /* serial output macros: */ @@ -33,5 +35,5 @@ #define NLHS 1 #undef NRHS -#define NRHS 8 +#define NRHS 10 Index: /issm/trunk/src/mex/Misfit/Misfit.cpp =================================================================== --- /issm/trunk/src/mex/Misfit/Misfit.cpp (revision 464) +++ /issm/trunk/src/mex/Misfit/Misfit.cpp (revision 465) @@ -15,8 +15,11 @@ DataSet* loads=NULL; DataSet* materials=NULL; - int analysis_type; double* u_g=NULL; double* u_g_obs=NULL; ParameterInputs* inputs=NULL; + char* analysis_type_string=NULL; + int analysis_type; + char* sub_analysis_type_string=NULL; + int sub_analysis_type; /* output datasets: */ @@ -34,7 +37,10 @@ FetchData((void**)&loads,NULL,NULL,LOADS,"DataSet",NULL); FetchData((void**)&materials,NULL,NULL,MATERIALS,"DataSet",NULL); - FetchData((void**)&analysis_type,NULL,NULL,mxGetField(PARAMETERS,0,"analysis_type"),"Integer",NULL); FetchData((void**)&u_g,NULL,NULL,UG,"Vector","Vec"); FetchData((void**)&u_g_obs,NULL,NULL,UGOBS,"Vector","Vec"); + FetchData((void**)&analysis_type_string,NULL,NULL,ANALYSIS,"String",NULL); + analysis_type=AnalysisTypeAsEnum(analysis_type_string); + FetchData((void**)&sub_analysis_type_string,NULL,NULL,SUBANALYSIS,"String",NULL); + sub_analysis_type=AnalysisTypeAsEnum(sub_analysis_type_string); /*Fetch inputs: */ @@ -43,5 +49,5 @@ /*!Call core code: */ - Misfitx(&J, elements,nodes,loads,materials,u_g,u_g_obs,inputs,analysis_type); + Misfitx(&J, elements,nodes,loads,materials,u_g,u_g_obs,inputs,analysis_type,sub_analysis_type); /*write output : */ @@ -55,5 +61,7 @@ xfree((void**)&u_g); xfree((void**)&u_g_obs); - delete inputs + delete inputs; + xfree((void**)&analysis_type_string); + xfree((void**)&sub_analysis_type_string); /*end module: */ Index: /issm/trunk/src/mex/Misfit/Misfit.h =================================================================== --- /issm/trunk/src/mex/Misfit/Misfit.h (revision 464) +++ /issm/trunk/src/mex/Misfit/Misfit.h (revision 465) @@ -25,4 +25,6 @@ #define UGOBS (mxArray*)prhs[6] #define INPUTS (mxArray*)prhs[7] +#define ANALYSIS (mxArray*)prhs[8] +#define SUBANALYSIS (mxArray*)prhs[9] /* serial output macros: */ @@ -33,5 +35,5 @@ #define NLHS 1 #undef NRHS -#define NRHS 8 +#define NRHS 10 Index: /issm/trunk/src/mex/ModelProcessor/ModelProcessor.cpp =================================================================== --- /issm/trunk/src/mex/ModelProcessor/ModelProcessor.cpp (revision 464) +++ /issm/trunk/src/mex/ModelProcessor/ModelProcessor.cpp (revision 465) @@ -32,13 +32,5 @@ /*Create elements, nodes and materials: */ - CreateElementsNodesAndMaterials(&elements,&nodes,&materials,model,MODEL); - - /*Create constraints: */ - CreateConstraints(&constraints,model,MODEL); - - /*Create loads: */ - CreateLoads(&loads,model,MODEL); - /*Create parameters: */ - CreateParameters(¶meters,model,MODEL); + CreateDataSets(&elements,&nodes,&materials,&constraints, &loads, ¶meters, model,MODEL); /*Write output data: */ Index: /issm/trunk/src/mex/PenaltyConstraints/PenaltyConstraints.cpp =================================================================== --- /issm/trunk/src/mex/PenaltyConstraints/PenaltyConstraints.cpp (revision 464) +++ /issm/trunk/src/mex/PenaltyConstraints/PenaltyConstraints.cpp (revision 465) @@ -16,5 +16,8 @@ DataSet* materials=NULL; ParameterInputs* inputs=NULL; + char* analysis_type_string=NULL; int analysis_type; + char* sub_analysis_type_string=NULL; + int sub_analysis_type; /*output: */ @@ -35,5 +38,9 @@ /*parameters: */ - FetchData((void**)&analysis_type,NULL,NULL,mxGetField(PARAMETERS,0,"analysis_type"),"Integer",NULL); + FetchData((void**)&analysis_type_string,NULL,NULL,ANALYSIS,"String",NULL); + analysis_type=AnalysisTypeAsEnum(analysis_type_string); + + FetchData((void**)&sub_analysis_type_string,NULL,NULL,SUBANALYSIS,"String",NULL); + sub_analysis_type=AnalysisTypeAsEnum(sub_analysis_type_string); /*Fetch inputs: */ @@ -42,5 +49,5 @@ /*!Generate internal degree of freedom numbers: */ - PenaltyConstraintsx(&converged, &num_unstable_constraints, elements,nodes,loads,materials,inputs,analysis_type); + PenaltyConstraintsx(&converged, &num_unstable_constraints, elements,nodes,loads,materials,inputs,analysis_type,sub_analysis_type); /*write output datasets: */ @@ -56,4 +63,6 @@ delete materials; delete inputs; + xfree((void**)&analysis_type_string); + xfree((void**)&sub_analysis_type_string); /*end module: */ Index: /issm/trunk/src/mex/PenaltyConstraints/PenaltyConstraints.h =================================================================== --- /issm/trunk/src/mex/PenaltyConstraints/PenaltyConstraints.h (revision 464) +++ /issm/trunk/src/mex/PenaltyConstraints/PenaltyConstraints.h (revision 465) @@ -23,4 +23,6 @@ #define PARAMETERS (mxArray*)prhs[4] #define INPUTS (mxArray*)prhs[5] +#define ANALYSIS (mxArray*)prhs[6] +#define SUBANALYSIS (mxArray*)prhs[7] /* serial output macros: */ @@ -33,5 +35,5 @@ #define NLHS 3 #undef NRHS -#define NRHS 6 +#define NRHS 8 #endif /* _PENALTYCONSTRAINTS_H */ Index: /issm/trunk/src/mex/PenaltySystemMatrices/PenaltySystemMatrices.cpp =================================================================== --- /issm/trunk/src/mex/PenaltySystemMatrices/PenaltySystemMatrices.cpp (revision 464) +++ /issm/trunk/src/mex/PenaltySystemMatrices/PenaltySystemMatrices.cpp (revision 465) @@ -21,4 +21,6 @@ int analysis_type; char* analysis_type_string=NULL; + int sub_analysis_type; + char* sub_analysis_type_string=NULL; /*Boot module: */ @@ -42,4 +44,7 @@ analysis_type=AnalysisTypeAsEnum(analysis_type_string); + FetchData((void**)&sub_analysis_type_string,NULL,NULL,SUBANALYSIS,"String",NULL); + sub_analysis_type=AnalysisTypeAsEnum(sub_analysis_type_string); + /*Fetch inputs: */ @@ -48,5 +53,5 @@ /*!Generate stiffnesses from penalties: */ - PenaltySystemMatricesx(Kgg, pg,elements,nodes,loads,materials,kflag,pflag,inputs,analysis_type); + PenaltySystemMatricesx(Kgg, pg,elements,nodes,loads,materials,kflag,pflag,inputs,analysis_type,sub_analysis_type); /*write output datasets: */ @@ -63,4 +68,5 @@ VecFree(&pg); xfree((void**)&analysis_type_string); + xfree((void**)&sub_analysis_type_string); /*end module: */ Index: /issm/trunk/src/mex/PenaltySystemMatrices/PenaltySystemMatrices.h =================================================================== --- /issm/trunk/src/mex/PenaltySystemMatrices/PenaltySystemMatrices.h (revision 464) +++ /issm/trunk/src/mex/PenaltySystemMatrices/PenaltySystemMatrices.h (revision 465) @@ -26,4 +26,5 @@ #define INPUTS (mxArray*)prhs[7] #define ANALYSIS (mxArray*)prhs[8] +#define SUBANALYSIS (mxArray*)prhs[9] /* serial output macros: */ @@ -35,5 +36,5 @@ #define NLHS 2 #undef NRHS -#define NRHS 9 +#define NRHS 10 Index: /issm/trunk/src/mex/SystemMatrices/SystemMatrices.cpp =================================================================== --- /issm/trunk/src/mex/SystemMatrices/SystemMatrices.cpp (revision 464) +++ /issm/trunk/src/mex/SystemMatrices/SystemMatrices.cpp (revision 465) @@ -21,4 +21,6 @@ char* analysis_type_string=NULL; int analysis_type; + char* sub_analysis_type_string=NULL; + int sub_analysis_type; /* output datasets: */ @@ -46,4 +48,7 @@ analysis_type=AnalysisTypeAsEnum(analysis_type_string); + FetchData((void**)&sub_analysis_type_string,NULL,NULL,SUBANALYSIS,"String",NULL); + sub_analysis_type=AnalysisTypeAsEnum(sub_analysis_type_string); + /*Fetch inputs: */ inputs=new ParameterInputs; @@ -51,5 +56,5 @@ /*!Generate internal degree of freedom numbers: */ - SystemMatricesx(&Kgg, &pg,elements,nodes,loads,materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type); + SystemMatricesx(&Kgg, &pg,elements,nodes,loads,materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type,sub_analysis_type); /*write output datasets: */ @@ -60,4 +65,5 @@ /*Free ressources: */ xfree((void**)&analysis_type_string); + xfree((void**)&sub_analysis_type_string); delete elements; delete nodes; Index: /issm/trunk/src/mex/SystemMatrices/SystemMatrices.h =================================================================== --- /issm/trunk/src/mex/SystemMatrices/SystemMatrices.h (revision 464) +++ /issm/trunk/src/mex/SystemMatrices/SystemMatrices.h (revision 465) @@ -24,4 +24,5 @@ #define INPUTS (mxArray*)prhs[5] #define ANALYSIS (mxArray*)prhs[6] +#define SUBANALYSIS (mxArray*)prhs[7] /* serial output macros: */ @@ -33,5 +34,5 @@ #define NLHS 2 #undef NRHS -#define NRHS 7 +#define NRHS 8