Changeset 465


Ignore:
Timestamp:
05/18/09 09:43:19 (16 years ago)
Author:
Eric.Larour
Message:

New sub analysis type field in all solutions. New thermal solution sequence

Location:
issm/trunk/src
Files:
4 added
3 deleted
97 edited
1 moved

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/DataSet/DataSet.cpp

    r444 r465  
    900900#undef __FUNCT__
    901901#define __FUNCT__ "DataSet::CreateKMatrix"
    902 void  DataSet::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){
     902void  DataSet::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
    903903
    904904        vector<Object*>::iterator object;
    … …  
    911911
    912912                        element=(Element*)(*object);
    913                         element->CreateKMatrix(Kgg,inputs,analysis_type);
     913                        element->CreateKMatrix(Kgg,inputs,analysis_type,sub_analysis_type);
    914914                }
    915915                if(EnumIsLoad((*object)->Enum())){
    916916
    917917                        load=(Load*)(*object);
    918                         load->CreateKMatrix(Kgg,inputs,analysis_type);
     918                        load->CreateKMatrix(Kgg,inputs,analysis_type,sub_analysis_type);
    919919                }
    920920        }
    … …  
    924924#undef __FUNCT__
    925925#define __FUNCT__ "DataSet::CreatePVector"
    926 void  DataSet::CreatePVector(Vec pg,void* inputs,int analysis_type){
     926void  DataSet::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){
    927927
    928928        vector<Object*>::iterator object;
    … …  
    935935
    936936                        element=(Element*)(*object);
    937                         element->CreatePVector(pg,inputs,analysis_type);
     937                        element->CreatePVector(pg,inputs,analysis_type,sub_analysis_type);
    938938                }
    939939                if(EnumIsLoad((*object)->Enum())){
    940940
    941941                        load=(Load*)(*object);
    942                         load->CreatePVector(pg,inputs,analysis_type);
     942                        load->CreatePVector(pg,inputs,analysis_type,sub_analysis_type);
    943943                }               
    944944        }
    … …  
    984984#undef __FUNCT__
    985985#define __FUNCT__ "DataSet::PenaltyCreateKMatrix"
    986 void  DataSet::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type){
     986void  DataSet::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
    987987
    988988        vector<Object*>::iterator object;
    … …  
    994994
    995995                        load=(Load*)(*object);
    996                         load->PenaltyCreateKMatrix(Kgg,inputs,kmax,analysis_type);
     996                        load->PenaltyCreateKMatrix(Kgg,inputs,kmax,analysis_type,sub_analysis_type);
    997997                }
    998998        }
    … …  
    10021002#undef __FUNCT__
    10031003#define __FUNCT__ "DataSet::PenaltyCreatePVector"
    1004 void  DataSet::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type){
     1004void  DataSet::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
    10051005
    10061006        vector<Object*>::iterator object;
    … …  
    10121012
    10131013                        load=(Load*)(*object);
    1014                         load->PenaltyCreatePVector(pg,inputs,kmax,analysis_type);
     1014                        load->PenaltyCreatePVector(pg,inputs,kmax,analysis_type,sub_analysis_type);
    10151015                }               
    10161016        }
    … …  
    10281028#undef __FUNCT__
    10291029#define __FUNCT__ "RiftConstraints"
    1030 void  DataSet::RiftConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type){
     1030void  DataSet::RiftConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type){
    10311031       
    10321032        throw ErrorException(__FUNCT__," not implemented yet!");
    … …  
    10451045#undef __FUNCT__
    10461046#define __FUNCT__ "MeltingConstraints"
    1047 void  DataSet::MeltingConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type){
     1047void  DataSet::MeltingConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type){
    10481048
    10491049        throw ErrorException(__FUNCT__," not implemented yet!");
    … …  
    10841084
    10851085
    1086 void  DataSet::Du(Vec du_g,double*  u_g,double* u_g_obs,void* inputs,int analysis_type){
     1086void  DataSet::Du(Vec du_g,double*  u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type){
    10871087
    10881088
    … …  
    10951095
    10961096                        element=(Element*)(*object);
    1097                         element->Du(du_g,u_g,u_g_obs,inputs,analysis_type);
    1098                 }
    1099         }
    1100 
    1101 
    1102 }
    1103 
    1104 void  DataSet::Gradj(Vec grad_g,double*  u_g,double* lambda_g,void* inputs,int analysis_type,char* control_type){
     1097                        element->Du(du_g,u_g,u_g_obs,inputs,analysis_type,sub_analysis_type);
     1098                }
     1099        }
     1100
     1101
     1102}
     1103
     1104void  DataSet::Gradj(Vec grad_g,double*  u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type){
    11051105
    11061106
    … …  
    11131113
    11141114                        element=(Element*)(*object);
    1115                         element->Gradj(grad_g,u_g,lambda_g,inputs,analysis_type,control_type);
     1115                        element->Gradj(grad_g,u_g,lambda_g,inputs,analysis_type,sub_analysis_type,control_type);
    11161116                }
    11171117        }
    … …  
    11201120}               
    11211121               
    1122 void  DataSet::Misfit(double* pJ, double* u_g,double* u_g_obs,void* inputs,int analysis_type){
     1122void  DataSet::Misfit(double* pJ, double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type){
    11231123
    11241124        double J=0;;
    … …  
    11321132
    11331133                        element=(Element*)(*object);
    1134                         J+=element->Misfit(u_g,u_g_obs,inputs,analysis_type);
     1134                        J+=element->Misfit(u_g,u_g_obs,inputs,analysis_type,sub_analysis_type);
    11351135
    11361136                }

  • issm/trunk/src/c/DataSet/DataSet.h

    r358 r465  
    4747                void  DistributeDofs(int numberofnodes,int numdofspernode);
    4848                void  CreatePartitioningVector(Vec* ppartition,int numnods,int numdofspernode);
    49                 void  DistributeNumDofs(int** pnumdofspernode,int numberofnodes,int analysis_type);
     49                void  DistributeNumDofs(int** pnumdofspernode,int numberofnodes,int analysis_type,int sub_analysis_type);
    5050                void  FlagClones(int numberofnodes);
    5151                int   NumberOfDofs();
    … …  
    6161                void  SetSorting(int* in_sorted_ids,int* in_id_offsets);
    6262                void  Sort();
    63                 void  CreateKMatrix(Mat Kgg,void* inputs, int analysis_type);
    64                 void  CreatePVector(Vec pg,void* inputs, int analysis_type);
     63                void  CreateKMatrix(Mat Kgg,void* inputs, int analysis_type,int sub_analysis_type);
     64                void  CreatePVector(Vec pg,void* inputs, int analysis_type,int sub_analysis_type);
    6565                void  UpdateFromInputs(void* inputs);
    66                 void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type);
    67                 void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type);
     66                void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
     67                void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
    6868                int   RiftIsPresent();
    69                 void  RiftConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type);
     69                void  RiftConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type);
    7070                int   MeltingIsPresent();
    71                 void  MeltingConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type);
     71                void  MeltingConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type);
    7272                DataSet* Copy(void);
    73                 void  Du(Vec du_g,double*  u_g,double* u_g_obs,void* inputs,int analysis_type);
    74                 void  Gradj(Vec grad_g,double*  u_g,double* lambda_g,void* inputs,int analysis_type,char* control_type);
    75                 void  Misfit(double* pJ, double* u_g,double* u_g_obs,void* inputs,int analysis_type);
     73                void  Du(Vec du_g,double*  u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type);
     74                void  Gradj(Vec grad_g,double*  u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type);
     75                void  Misfit(double* pJ, double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type);
    7676                void  VelocityExtrude(Vec ug,double* ug_serial);
    7777                void  SlopeExtrude(Vec sg,double* sg_serial);

  • issm/trunk/src/c/Dux/Dux.cpp

    r1 r465  
    1414
    1515void Dux( Vec* pdu_g, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials,
    16                 double* u_g, double* u_g_obs, ParameterInputs* inputs,int analysis_type){
     16                double* u_g, double* u_g_obs, ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
    1717
    1818        int i;
    … …  
    3333
    3434        /*Compute velocity difference: */
    35         elements->Du(du_g, u_g,u_g_obs,inputs,analysis_type);
     35        elements->Du(du_g, u_g,u_g_obs,inputs,analysis_type,sub_analysis_type);
    3636
    3737        /*Assemble vector: */

  • issm/trunk/src/c/Dux/Dux.h

    r1 r465  
    1010/* local prototypes: */
    1111void Dux( Vec* pdu_g, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials,
    12                 double* u_g, double* u_g_obs, ParameterInputs* inputs,int analysis_type);
     12                double* u_g, double* u_g_obs, ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
    1313
    1414#endif  /* _DUX_H */

  • issm/trunk/src/c/EnumDefinitions/AnalysisTypeAsEnum.cpp

    r358 r465  
    1515        }
    1616
    17         if (strcmp(analysis_type,"control")==0){
     17        if (strcmp(analysis_type,"diagnostic")==0){
     18                return DiagnosticAnalysisEnum();
     19        }
     20        else if (strcmp(analysis_type,"control")==0){
    1821                return ControlAnalysisEnum();
    1922        }
    20         else if (strcmp(analysis_type,"diagnostic_horiz")==0){
    21                 return DiagnosticHorizAnalysisEnum();
    22         }
    23         else if (strcmp(analysis_type,"diagnostic_vert")==0){
    24                 return DiagnosticVertAnalysisEnum();
     23        else if (strcmp(analysis_type,"thermal")==0){
     24                return ThermalAnalysisEnum();
    2525        }
    2626        else if (strcmp(analysis_type,"prognostic")==0){
    2727                return PrognosticAnalysisEnum();
    2828        }
    29         else if (strcmp(analysis_type,"thermalsteady")==0){
    30                 return ThermalSteadyAnalysisEnum();
    31         }
    32         else if (strcmp(analysis_type,"thermaltransient")==0){
    33                 return ThermalTransientAnalysisEnum();
    34         }
    3529        else if (strcmp(analysis_type,"melting")==0){
    3630                return MeltingAnalysisEnum();
    37         }
    38         else if (strcmp(analysis_type,"diagnostic_stokes")==0){
    39                 return DiagnosticStokesAnalysisEnum();
    40         }
    41         else if (strcmp(analysis_type,"diagnostic_hutter")==0){
    42                 return DiagnosticHutterAnalysisEnum();
    4331        }
    4432        else if (strcmp(analysis_type,"slope_compute")==0){
    4533                return SlopeComputeAnalysisEnum();
    4634        }
    47         else if (strcmp(analysis_type,"surface_slope_compute_x")==0){
    48                 return SurfaceSlopeComputeXAnalysisEnum();
     35        else if (strcmp(analysis_type,"stokes")==0){
     36                return StokesAnalysisEnum();
    4937        }
    50         else if (strcmp(analysis_type,"surface_slope_compute_y")==0){
    51                 return SurfaceSlopeComputeYAnalysisEnum();
     38        else if (strcmp(analysis_type,"hutter")==0){
     39                return HutterAnalysisEnum();
    5240        }
    53         else if (strcmp(analysis_type,"bed_slope_compute_x")==0){
    54                 return BedSlopeComputeXAnalysisEnum();
     41        else if (strcmp(analysis_type,"surfacex")==0){
     42                return SurfaceXAnalysisEnum();
    5543        }
    56         else if (strcmp(analysis_type,"bed_slope_compute_y")==0){
    57                 return BedSlopeComputeYAnalysisEnum();
     44        else if (strcmp(analysis_type,"surfacey")==0){
     45                return SurfaceYAnalysisEnum();
    5846        }
    59         else throw ErrorException(__FUNCT__," could not recognized analysis type!");
     47        else if (strcmp(analysis_type,"bedx")==0){
     48                return BedXAnalysisEnum();
     49        }
     50        else if (strcmp(analysis_type,"bedy")==0){
     51                return BedYAnalysisEnum();
     52        }
     53        else if (strcmp(analysis_type,"horiz")==0){
     54                return HorizAnalysisEnum();
     55        }
     56        else if (strcmp(analysis_type,"vert")==0){
     57                return VertAnalysisEnum();
     58        }
     59        else if (strcmp(analysis_type,"")==0){
     60                return NoneAnalysisEnum();
     61        }
     62        else throw ErrorException(__FUNCT__,exprintf("%s%s"," could not recognize analysis type: ",analysis_type));
    6063
    6164}

  • issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp

    r387 r465  
    1515
    1616/*analysis types: */
    17 int StructuralAnalysisEnum(void){       return          20; }
    18 int ThermalAnalysisEnum(void){          return          21; }
    19 int ControlAnalysisEnum(void){          return          22; }
    20 int DiagnosticHorizAnalysisEnum(void){  return          23; }
    21 int DiagnosticVertAnalysisEnum(void){   return          24; }
    22 int PrognosticAnalysisEnum(void){       return          26; }
    23 int ThermalSteadyAnalysisEnum(void)     {   return          27; }
    24 int ThermalTransientAnalysisEnum(void){ return          28; }
    25 int MeltingAnalysisEnum(void){          return          29; }
    26 int DiagnosticStokesAnalysisEnum(void){ return          30; }
    27 int DiagnosticHutterAnalysisEnum(void){ return          31; }
    28 int SlopeComputeAnalysisEnum(void){ return              32; }
    29 int SurfaceSlopeComputeXAnalysisEnum(void){ return      33; }
    30 int SurfaceSlopeComputeYAnalysisEnum(void){ return      34; }
    31 int BedSlopeComputeXAnalysisEnum(void){ return          35; }
    32 int BedSlopeComputeYAnalysisEnum(void){ return          36; }
     17int DiagnosticAnalysisEnum(void){  return               20; }
     18int ControlAnalysisEnum(void){          return          21; }
     19int ThermalAnalysisEnum(void){          return          22; }
     20int PrognosticAnalysisEnum(void){          return       23; }
     21int MeltingAnalysisEnum(void){          return          24; }
     22int SlopeComputeAnalysisEnum(void){ return              25; }
     23int StokesAnalysisEnum(void){ return                    26; }
     24int HutterAnalysisEnum(void){ return                    27; }
     25int SurfaceXAnalysisEnum(void){ return                  28; }
     26int SurfaceYAnalysisEnum(void){ return                  29; }
     27int BedXAnalysisEnum(void){ return                      30; }
     28int BedYAnalysisEnum(void){ return                      31; }
     29int HorizAnalysisEnum(void){ return                     32; }
     30int VertAnalysisEnum(void){ return                      33; }
     31int NoneAnalysisEnum(void){ return                      34; }
    3332
    3433/*datasets: */

  • issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h

    r387 r465  
    3434
    3535/* analysis enums: */
     36int DiagnosticAnalysisEnum(void);
    3637int ControlAnalysisEnum(void);
    37 int DiagnosticHorizAnalysisEnum(void);
    38 int DiagnosticVertAnalysisEnum(void);
     38int ThermalAnalysisEnum(void);
    3939int PrognosticAnalysisEnum(void);
    40 int ThermalSteadyAnalysisEnum(void);
    41 int ThermalTransientAnalysisEnum(void);
    4240int MeltingAnalysisEnum(void);
    43 int DiagnosticStokesAnalysisEnum(void);
    44 int DiagnosticHutterAnalysisEnum(void);
    4541int SlopeComputeAnalysisEnum(void);
    46 int SurfaceSlopeComputeXAnalysisEnum(void);
    47 int SurfaceSlopeComputeYAnalysisEnum(void);
    48 int BedSlopeComputeXAnalysisEnum(void);
    49 int BedSlopeComputeYAnalysisEnum(void);
     42int StokesAnalysisEnum(void);
     43int HutterAnalysisEnum(void);
     44int SurfaceXAnalysisEnum(void);
     45int SurfaceYAnalysisEnum(void);
     46int BedXAnalysisEnum(void);
     47int BedYAnalysisEnum(void);
     48int HorizAnalysisEnum(void);
     49int VertAnalysisEnum(void);
     50int NoneAnalysisEnum(void);
    5051
    5152

  • issm/trunk/src/c/Gradjx/Gradjx.cpp

    r1 r465  
    1414
    1515void Gradjx( Vec* pgrad_g, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials,
    16                 double* u_g, double* lambda_g, ParameterInputs* inputs,int analysis_type,char* control_type){
     16                double* u_g, double* lambda_g, ParameterInputs* inputs,int analysis_type,int sub_analysis_type,char* control_type){
    1717
    1818        int i;
    … …  
    3333
    3434        /*Compute gradients: */
    35         elements->Gradj(grad_g, u_g,lambda_g,inputs,analysis_type,control_type);
     35        elements->Gradj(grad_g, u_g,lambda_g,inputs,analysis_type,sub_analysis_type,control_type);
    3636
    3737        /*Assemble vector: */

  • issm/trunk/src/c/Gradjx/Gradjx.h

    r1 r465  
    1010/* local prototypes: */
    1111void Gradjx( Vec* pgrad_g, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials,
    12                 double* u_g, double* lambda_g, ParameterInputs* inputs,int analysis_type,char* control_type);
     12                double* u_g, double* lambda_g, ParameterInputs* inputs,int analysis_type,int sub_analysis_type,char* control_type);
    1313
    1414#endif  /* _GRADJX_H */

  • issm/trunk/src/c/Makefile.am

    r434 r465  
    7979                                        ./shared/Dofs/dofs.h\
    8080                                        ./shared/Dofs/dofsetgen.cpp\
     81                                        ./shared/Dofs/DistributeNumDofs.cpp\
    8182                                        ./shared/Numerics/numerics.h\
    8283                                        ./shared/Numerics/GaussPoints.h\
    … …  
    8687                                        ./shared/Numerics/BrentSearch.cpp\
    8788                                        ./shared/Numerics/OptFunc.cpp\
    88                                         ./shared/Numerics/extrema.cpp\
    89                                         ./shared/Numerics/DistributeNumDofs.cpp\
    90                                         ./shared/Exceptions/exceptions.h\
    91                                         ./shared/Exceptions/Exceptions.cpp\
    92                                         ./shared/Exceptions/exprintf.cpp\
    93                                         ./shared/Exp/exp.h\
    94                                         ./shared/Exp/IsInPoly.cpp\
    95                                         ./shared/Exp/IsInPolySerial.cpp\
    96                                         ./shared/Exp/DomainOutlineRead.cpp\
    97                                         ./shared/TriMesh/trimesh.h\
    98                                         ./shared/TriMesh/AssociateSegmentToElement.cpp\
    99                                         ./shared/TriMesh/GridInsideHole.cpp\
    100                                         ./shared/TriMesh/OrderSegments.cpp\
    101                                         ./shared/TriMesh/SplitMeshForRifts.cpp\
    102                                         ./shared/TriMesh/TriMeshUtils.cpp\
    103                                         ./shared/Sorting/binary_search.cpp\
    104                                         ./shared/Sorting/sorting.h\
    105                                         ./shared/Elements/elements.h\
    106                                         ./shared/Elements/ResolvePointers.cpp\
    107                                         ./shared/Elements/Paterson.cpp\
    108                                         ./shared/Elements/GetElementNodeData.cpp\
    109                                         ./toolkits/petsc\
    110                                         ./toolkits/petsc/patches\
    111                                         ./toolkits/petsc/patches/petscpatches.h\
    112                                         ./toolkits/petsc/patches/MatlabMatrixToPetscMatrix.cpp\
    113                                         ./toolkits/petsc/patches/MatlabVectorToPetscVector.cpp\
    114                                         ./toolkits/petsc/patches/PetscMatrixToMatlabMatrix.cpp\
    115                                         ./toolkits/petsc/patches/PetscVectorToMatlabVector.cpp\
    116                                         ./toolkits/petsc/patches/MatlabMatrixToDoubleMatrix.cpp\
    117                                         ./toolkits/petsc/patches/MatlabVectorToDoubleVector.cpp\
    118                                         ./toolkits/petsc/patches/PetscDetermineLocalSize.cpp\
    119                                         ./toolkits/petsc/patches/VecTranspose.cpp\
    120                                         ./toolkits/petsc/patches/VecToMPISerial.cpp\
    121                                         ./toolkits/petsc/patches/MatToSerial.cpp\
    122                                         ./toolkits/petsc/patches/VecMerge.cpp\
    123                                         ./toolkits/petsc/patches/NewVec.cpp\
    124                                         ./toolkits/petsc/patches/NewVecFromLocalSize.cpp\
    125                                         ./toolkits/petsc/patches/PetscOptionsDetermineSolverType.cpp\
    126                                         ./toolkits/petsc/patches/PetscOptionsInsertMultipleString.cpp\
    127                                         ./toolkits/petsc/patches/NewMat.cpp\
    128                                         ./toolkits/petsc/patches/VecFree.cpp\
    129                                         ./toolkits/petsc/patches/VecDuplicatePatch.cpp\
    130                                         ./toolkits/petsc/patches/KSPFree.cpp\
    131                                         ./toolkits/petsc/patches/ISFree.cpp\
    132                                         ./toolkits/petsc/patches/MatFree.cpp\
    133                                         ./toolkits/petsc/patches/GetOwnershipBoundariesFromRange.cpp\
    134                                         ./toolkits/petsc/patches/VecPartition.cpp\
    135                                         ./toolkits/petsc/patches/MatPartition.cpp\
    136                                         ./toolkits/petsc/patches/MatInvert.cpp\
    137                                         ./toolkits/petsc/patches/MatMultPatch.cpp\
    138                                         ./toolkits/petsc/petscincludes.h\
    139                                         ./toolkits/mpi/mpiincludes.h\
    140                                         ./toolkits/mpi/patches/mpipatches.h\
    141                                         ./toolkits/mpi/patches/MPI_Upperrow.cpp\
    142                                         ./toolkits/mpi/patches/MPI_Lowerrow.cpp\
    143                                         ./toolkits/mpi/patches/MPI_Boundariesfromrange.cpp\
    144                                         ./toolkits/metis/metisincludes.h\
    145                                         ./toolkits/triangle/triangleincludes.h\
    146                                         ./toolkits.h\
    147                                         ./io/io.h\
    148                                         ./io/FetchData.cpp\
    149                                         ./io/WriteData.cpp\
    150                                         ./io/WriteDataToDisk.cpp\
    151                                         ./io/SerialFetchData.cpp\
    152                                         ./io/SerialWriteData.cpp\
    153                                         ./io/ParallelFetchData.cpp\
    154                                         ./io/ParallelFetchInteger.cpp\
    155                                         ./io/ParallelFetchMat.cpp\
    156                                         ./io/ParallelFetchScalar.cpp\
    157                                         ./io/ParallelFetchString.cpp\
    158                                         ./io/ModelFetchData.cpp\
    159                                         ./io/WriteNodeSets.cpp\
    160                                         ./io/WriteParams.cpp\
    161                                         ./io/FetchNodeSets.cpp\
    162                                         ./io/FetchRifts.cpp\
    163                                         ./io/ParameterInputsInit.cpp\
    164                                         ./EnumDefinitions/EnumDefinitions.h\
    165                                         ./EnumDefinitions/EnumDefinitions.cpp\
    166                                         ./EnumDefinitions/AnalysisTypeAsEnum.cpp\
    167                                         ./ModelProcessorx/Model.h\
    168                                         ./ModelProcessorx/Model.cpp\
    169                                         ./ModelProcessorx/CreateElementsNodesAndMaterials.cpp\
    170                                         ./ModelProcessorx/CreateConstraints.cpp \
    171                                         ./ModelProcessorx/CreateLoads.cpp\
    172                                         ./ModelProcessorx/CreateParameters.cpp\
    173                                         ./ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp\
    174                                         ./ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp \
    175                                         ./ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp\
    176                                         ./ModelProcessorx/DiagnosticHoriz/CreateParametersDiagnosticHoriz.cpp\
    177                                         ./ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp\
    178                                         ./ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp \
    179                                         ./ModelProcessorx/DiagnosticVert/CreateLoadsDiagnosticVert.cpp\
    180                                         ./ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp\
    181                                         ./ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp \
    182                                         ./ModelProcessorx/DiagnosticHutter/CreateLoadsDiagnosticHutter.cpp\
    183                                         ./ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp\
    184                                         ./ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp \
    185                                         ./ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp\
    186                                         ./ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp\
    187                                         ./ModelProcessorx/SlopeCompute/CreateConstraintsSlopeCompute.cpp \
    188                                         ./ModelProcessorx/SlopeCompute/CreateLoadsSlopeCompute.cpp\
    189                                         ./ModelProcessorx/Control/CreateParametersControl.cpp\
    190                                         ./Dofx/Dofx.h\
    191                                         ./Dofx/Dofx.cpp\
    192                                         ./Dux/Dux.h\
    193                                         ./Dux/Dux.cpp\
    194                                         ./GriddataMeshToGridx/GriddataMeshToGridx.h\
    195                                         ./GriddataMeshToGridx/GriddataMeshToGridx.cpp\
    196                                         ./ControlConstrainx/ControlConstrainx.h\
    197                                         ./ControlConstrainx/ControlConstrainx.cpp\
    198                                         ./Misfitx/Misfitx.h\
    199                                         ./Misfitx/Misfitx.cpp\
    200                                         ./Orthx/Orthx.h\
    201                                         ./Orthx/Orthx.cpp\
    202                                         ./Gradjx/Gradjx.h\
    203                                         ./Gradjx/Gradjx.cpp\
    204                                         ./UpdateFromInputsx/UpdateFromInputsx.h\
    205                                         ./UpdateFromInputsx/UpdateFromInputsx.cpp\
    206                                         ./ConfigureObjectsx/ConfigureObjectsx.h\
    207                                         ./ConfigureObjectsx/ConfigureObjectsx.cpp\
    208                                         ./ComputePressurex/ComputePressurex.h\
    209                                         ./ComputePressurex/ComputePressurex.cpp\
    210                                         ./BuildNodeSetsx/BuildNodeSetsx.h\
    211                                         ./BuildNodeSetsx/BuildNodeSetsx.cpp\
    212                                         ./BuildNodeSetsx/PartitionSets.cpp\
    213                                         ./SpcNodesx/SpcNodesx.h\
    214                                         ./SpcNodesx/SpcNodesx.cpp\
    215                                         ./MpcNodesx/MpcNodesx.h\
    216                                         ./MpcNodesx/MpcNodesx.cpp\
    217                                         ./DataInterpx/DataInterpx.cpp\
    218                                         ./DataInterpx/DataInterpx.h\
    219                                         ./HoleFillerx/HoleFillerx.cpp\
    220                                         ./HoleFillerx/HoleFillerx.h\
    221                                         ./MeshPartitionx/MeshPartitionx.cpp\
    222                                         ./MeshPartitionx/MeshPartitionx.h\
    223                                         ./ContourToMeshx/ContourToMeshx.cpp\
    224                                         ./ContourToMeshx/ContourToMeshx.h\
    225                                         ./Reducevectorgtosx/Reducevectorgtosx.cpp\
    226                                         ./Reducevectorgtosx/Reducevectorgtosx.h\
    227                                         ./Reducematrixfromgtofx/Reducematrixfromgtofx.cpp\
    228                                         ./Reducematrixfromgtofx/Reducematrixfromgton.cpp\
    229                                         ./Reducematrixfromgtofx/Reducematrixfromgtofx.h\
    230                                         ./Reduceloadfromgtofx/Reduceloadfromgtofx.h\
    231                                         ./Reduceloadfromgtofx/Reduceloadfromgtofx.cpp\
    232                                         ./NormalizeConstraintsx/NormalizeConstraintsx.cpp\
    233                                         ./NormalizeConstraintsx/NormalizeConstraintsx.h\
    234                                         ./SystemMatricesx/SystemMatricesx.cpp\
    235                                         ./SystemMatricesx/SystemMatricesx.h\
    236                                         ./PenaltyConstraintsx/PenaltyConstraintsx.cpp\
    237                                         ./PenaltyConstraintsx/PenaltyConstraintsx.h\
    238                                         ./PenaltySystemMatricesx/PenaltySystemMatricesx.cpp\
    239                                         ./PenaltySystemMatricesx/PenaltySystemMatricesx.h\
    240                                         ./Solverx/Solverx.cpp\
    241                                         ./Solverx/Solverx.h\
    242                                         ./Mergesolutionfromftogx/Mergesolutionfromftogx.cpp\
    243                                         ./Mergesolutionfromftogx/Mergesolutionfromftogx.h\
    244                                         ./ProcessParamsx/ProcessParamsx.cpp\
    245                                         ./ProcessParamsx/ProcessParamsx.h\
    246                                         ./VelocityExtrudex/VelocityExtrudex.cpp\
    247                                         ./VelocityExtrudex/VelocityExtrudex.h\
    248                                         ./SlopeExtrudex/SlopeExtrudex.cpp\
    249                                         ./SlopeExtrudex/SlopeExtrudex.h
    250 
    251 
    252 
    253 
    254 libISSM_a_CXXFLAGS = -fPIC -DMATLAB -D_SERIAL_ -ansi -D_GNU_SOURCE -fno-omit-frame-pointer -pthread -D_CPP_
    255 if LARGEARRAYS
    256 libISSM_a_CXXFLAGS += -D__GCC4BUILD__ 
    257 else
    258 libISSM_a_CXXFLAGS += -DMX_COMPAT_32
    259 endif
    260 
    261 
    262 
    263 
    264 #Parallel compilation
    265 libpISSM_a_SOURCES = ./objects/objects.h\
    266                                         ./objects/Object.h\
    267                                         ./objects/Element.h\
    268                                         ./objects/Element.cpp\
    269                                         ./objects/Material.h\
    270                                         ./objects/Material.cpp\
    271                                         ./objects/Load.h\
    272                                         ./objects/Load.cpp\
    273                                         ./objects/SolverEnum.h\
    274                                         ./objects/Contour.h\
    275                                         ./objects/Contour.cpp\
    276                                         ./objects/OptArgs.h\
    277                                         ./objects/OptPars.h\
    278                                         ./objects/Friction.h\
    279                                         ./objects/Friction.cpp\
    280                                         ./objects/Node.h\
    281                                         ./objects/Node.cpp\
    282                                         ./objects/Tria.h\
    283                                         ./objects/Tria.cpp\
    284                                         ./objects/Sing.h\
    285                                         ./objects/Sing.cpp\
    286                                         ./objects/Beam.h\
    287                                         ./objects/Beam.cpp\
    288                                         ./objects/Penta.h\
    289                                         ./objects/Penta.cpp\
    290                                         ./objects/Matice.h\
    291                                         ./objects/Matice.cpp\
    292                                         ./objects/Matpar.h\
    293                                         ./objects/Matpar.cpp\
    294                                         ./objects/Input.h\
    295                                         ./objects/Input.cpp\
    296                                         ./objects/ParameterInputs.h\
    297                                         ./objects/ParameterInputs.cpp\
    298                                         ./objects/Spc.cpp\
    299                                         ./objects/Spc.h\
    300                                         ./objects/Rgb.cpp\
    301                                         ./objects/Rgb.h\
    302                                         ./objects/Penpair.cpp\
    303                                         ./objects/Penpair.h\
    304                                         ./objects/Pengrid.cpp\
    305                                         ./objects/Pengrid.h\
    306                                         ./objects/Icefront.cpp\
    307                                         ./objects/Icefront.h\
    308                                         ./objects/Param.cpp\
    309                                         ./objects/Param.h\
    310                                         ./objects/NodeSets.cpp\
    311                                         ./objects/NodeSets.h\
    312                                         ./DataSet/DataSet.cpp\
    313                                         ./DataSet/DataSet.h\
    314                                         ./shared/shared.h\
    315                                         ./shared/Alloc/alloc.h\
    316                                         ./shared/Alloc/alloc.cpp\
    317                                         ./shared/Matlab/matlabshared.h\
    318                                         ./shared/Matlab/PrintfFunction.cpp\
    319                                         ./shared/Matlab/ModuleBoot.cpp\
    320                                         ./shared/Matlab/mxGetAssignedField.cpp\
    321                                         ./shared/Matlab/mxGetField.cpp\
    322                                         ./shared/Matlab/CheckNumMatlabArguments.cpp\
    323                                         ./shared/Matrix/matrix.h\
    324                                         ./shared/Matrix/MatrixUtils.cpp\
    325                                         ./shared/Dofs/dofs.h\
    326                                         ./shared/Dofs/dofsetgen.cpp\
    327                                         ./shared/Numerics/numerics.h\
    328                                         ./shared/Numerics/GaussPoints.h\
    329                                         ./shared/Numerics/GaussPoints.cpp\
    330                                         ./shared/Numerics/cross.cpp\
    331                                         ./shared/Numerics/norm.cpp\
    332                                         ./shared/Numerics/BrentSearch.cpp\
    333                                         ./shared/Numerics/OptFunc.cpp\
    334                                         ./shared/Numerics/DistributeNumDofs.cpp\
    33589                                        ./shared/Numerics/extrema.cpp\
    33690                                        ./shared/Exceptions/exceptions.h\
    … …  
    413167                                        ./ModelProcessorx/Model.h\
    414168                                        ./ModelProcessorx/Model.cpp\
    415                                         ./ModelProcessorx/CreateElementsNodesAndMaterials.cpp\
    416                                         ./ModelProcessorx/CreateConstraints.cpp \
    417                                         ./ModelProcessorx/CreateLoads.cpp\
     169                                        ./ModelProcessorx/CreateDataSets.cpp\
     170                                        ./ModelProcessorx/CreateParameters.cpp\
     171                                        ./ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp\
     172                                        ./ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp \
     173                                        ./ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp\
     174                                        ./ModelProcessorx/DiagnosticHoriz/CreateParametersDiagnosticHoriz.cpp\
     175                                        ./ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp\
     176                                        ./ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp \
     177                                        ./ModelProcessorx/DiagnosticVert/CreateLoadsDiagnosticVert.cpp\
     178                                        ./ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp\
     179                                        ./ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp \
     180                                        ./ModelProcessorx/DiagnosticHutter/CreateLoadsDiagnosticHutter.cpp\
     181                                        ./ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp\
     182                                        ./ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp \
     183                                        ./ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp\
     184                                        ./ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp\
     185                                        ./ModelProcessorx/SlopeCompute/CreateConstraintsSlopeCompute.cpp \
     186                                        ./ModelProcessorx/SlopeCompute/CreateLoadsSlopeCompute.cpp\
     187                                        ./ModelProcessorx/Control/CreateParametersControl.cpp\
     188                                        ./Dofx/Dofx.h\
     189                                        ./Dofx/Dofx.cpp\
     190                                        ./Dux/Dux.h\
     191                                        ./Dux/Dux.cpp\
     192                                        ./GriddataMeshToGridx/GriddataMeshToGridx.h\
     193                                        ./GriddataMeshToGridx/GriddataMeshToGridx.cpp\
     194                                        ./ControlConstrainx/ControlConstrainx.h\
     195                                        ./ControlConstrainx/ControlConstrainx.cpp\
     196                                        ./Misfitx/Misfitx.h\
     197                                        ./Misfitx/Misfitx.cpp\
     198                                        ./Orthx/Orthx.h\
     199                                        ./Orthx/Orthx.cpp\
     200                                        ./Gradjx/Gradjx.h\
     201                                        ./Gradjx/Gradjx.cpp\
     202                                        ./UpdateFromInputsx/UpdateFromInputsx.h\
     203                                        ./UpdateFromInputsx/UpdateFromInputsx.cpp\
     204                                        ./ConfigureObjectsx/ConfigureObjectsx.h\
     205                                        ./ConfigureObjectsx/ConfigureObjectsx.cpp\
     206                                        ./ComputePressurex/ComputePressurex.h\
     207                                        ./ComputePressurex/ComputePressurex.cpp\
     208                                        ./BuildNodeSetsx/BuildNodeSetsx.h\
     209                                        ./BuildNodeSetsx/BuildNodeSetsx.cpp\
     210                                        ./BuildNodeSetsx/PartitionSets.cpp\
     211                                        ./SpcNodesx/SpcNodesx.h\
     212                                        ./SpcNodesx/SpcNodesx.cpp\
     213                                        ./MpcNodesx/MpcNodesx.h\
     214                                        ./MpcNodesx/MpcNodesx.cpp\
     215                                        ./DataInterpx/DataInterpx.cpp\
     216                                        ./DataInterpx/DataInterpx.h\
     217                                        ./HoleFillerx/HoleFillerx.cpp\
     218                                        ./HoleFillerx/HoleFillerx.h\
     219                                        ./MeshPartitionx/MeshPartitionx.cpp\
     220                                        ./MeshPartitionx/MeshPartitionx.h\
     221                                        ./ContourToMeshx/ContourToMeshx.cpp\
     222                                        ./ContourToMeshx/ContourToMeshx.h\
     223                                        ./Reducevectorgtosx/Reducevectorgtosx.cpp\
     224                                        ./Reducevectorgtosx/Reducevectorgtosx.h\
     225                                        ./Reducematrixfromgtofx/Reducematrixfromgtofx.cpp\
     226                                        ./Reducematrixfromgtofx/Reducematrixfromgton.cpp\
     227                                        ./Reducematrixfromgtofx/Reducematrixfromgtofx.h\
     228                                        ./Reduceloadfromgtofx/Reduceloadfromgtofx.h\
     229                                        ./Reduceloadfromgtofx/Reduceloadfromgtofx.cpp\
     230                                        ./NormalizeConstraintsx/NormalizeConstraintsx.cpp\
     231                                        ./NormalizeConstraintsx/NormalizeConstraintsx.h\
     232                                        ./SystemMatricesx/SystemMatricesx.cpp\
     233                                        ./SystemMatricesx/SystemMatricesx.h\
     234                                        ./PenaltyConstraintsx/PenaltyConstraintsx.cpp\
     235                                        ./PenaltyConstraintsx/PenaltyConstraintsx.h\
     236                                        ./PenaltySystemMatricesx/PenaltySystemMatricesx.cpp\
     237                                        ./PenaltySystemMatricesx/PenaltySystemMatricesx.h\
     238                                        ./Solverx/Solverx.cpp\
     239                                        ./Solverx/Solverx.h\
     240                                        ./Mergesolutionfromftogx/Mergesolutionfromftogx.cpp\
     241                                        ./Mergesolutionfromftogx/Mergesolutionfromftogx.h\
     242                                        ./ProcessParamsx/ProcessParamsx.cpp\
     243                                        ./ProcessParamsx/ProcessParamsx.h\
     244                                        ./VelocityExtrudex/VelocityExtrudex.cpp\
     245                                        ./VelocityExtrudex/VelocityExtrudex.h\
     246                                        ./SlopeExtrudex/SlopeExtrudex.cpp\
     247                                        ./SlopeExtrudex/SlopeExtrudex.h
     248
     249
     250
     251
     252libISSM_a_CXXFLAGS = -fPIC -DMATLAB -D_SERIAL_ -ansi -D_GNU_SOURCE -fno-omit-frame-pointer -pthread -D_CPP_
     253if LARGEARRAYS
     254libISSM_a_CXXFLAGS += -D__GCC4BUILD__ 
     255else
     256libISSM_a_CXXFLAGS += -DMX_COMPAT_32
     257endif
     258
     259
     260
     261
     262#Parallel compilation
     263libpISSM_a_SOURCES = ./objects/objects.h\
     264                                        ./objects/Object.h\
     265                                        ./objects/Element.h\
     266                                        ./objects/Element.cpp\
     267                                        ./objects/Material.h\
     268                                        ./objects/Material.cpp\
     269                                        ./objects/Load.h\
     270                                        ./objects/Load.cpp\
     271                                        ./objects/SolverEnum.h\
     272                                        ./objects/Contour.h\
     273                                        ./objects/Contour.cpp\
     274                                        ./objects/OptArgs.h\
     275                                        ./objects/OptPars.h\
     276                                        ./objects/Friction.h\
     277                                        ./objects/Friction.cpp\
     278                                        ./objects/Node.h\
     279                                        ./objects/Node.cpp\
     280                                        ./objects/Tria.h\
     281                                        ./objects/Tria.cpp\
     282                                        ./objects/Sing.h\
     283                                        ./objects/Sing.cpp\
     284                                        ./objects/Beam.h\
     285                                        ./objects/Beam.cpp\
     286                                        ./objects/Penta.h\
     287                                        ./objects/Penta.cpp\
     288                                        ./objects/Matice.h\
     289                                        ./objects/Matice.cpp\
     290                                        ./objects/Matpar.h\
     291                                        ./objects/Matpar.cpp\
     292                                        ./objects/Input.h\
     293                                        ./objects/Input.cpp\
     294                                        ./objects/ParameterInputs.h\
     295                                        ./objects/ParameterInputs.cpp\
     296                                        ./objects/Spc.cpp\
     297                                        ./objects/Spc.h\
     298                                        ./objects/Rgb.cpp\
     299                                        ./objects/Rgb.h\
     300                                        ./objects/Penpair.cpp\
     301                                        ./objects/Penpair.h\
     302                                        ./objects/Pengrid.cpp\
     303                                        ./objects/Pengrid.h\
     304                                        ./objects/Icefront.cpp\
     305                                        ./objects/Icefront.h\
     306                                        ./objects/Param.cpp\
     307                                        ./objects/Param.h\
     308                                        ./objects/NodeSets.cpp\
     309                                        ./objects/NodeSets.h\
     310                                        ./DataSet/DataSet.cpp\
     311                                        ./DataSet/DataSet.h\
     312                                        ./shared/shared.h\
     313                                        ./shared/Alloc/alloc.h\
     314                                        ./shared/Alloc/alloc.cpp\
     315                                        ./shared/Matlab/matlabshared.h\
     316                                        ./shared/Matlab/PrintfFunction.cpp\
     317                                        ./shared/Matlab/ModuleBoot.cpp\
     318                                        ./shared/Matlab/mxGetAssignedField.cpp\
     319                                        ./shared/Matlab/mxGetField.cpp\
     320                                        ./shared/Matlab/CheckNumMatlabArguments.cpp\
     321                                        ./shared/Matrix/matrix.h\
     322                                        ./shared/Matrix/MatrixUtils.cpp\
     323                                        ./shared/Dofs/dofs.h\
     324                                        ./shared/Dofs/dofsetgen.cpp\
     325                                        ./shared/Dofs/DistributeNumDofs.cpp\
     326                                        ./shared/Numerics/numerics.h\
     327                                        ./shared/Numerics/GaussPoints.h\
     328                                        ./shared/Numerics/GaussPoints.cpp\
     329                                        ./shared/Numerics/cross.cpp\
     330                                        ./shared/Numerics/norm.cpp\
     331                                        ./shared/Numerics/BrentSearch.cpp\
     332                                        ./shared/Numerics/OptFunc.cpp\
     333                                        ./shared/Numerics/extrema.cpp\
     334                                        ./shared/Exceptions/exceptions.h\
     335                                        ./shared/Exceptions/Exceptions.cpp\
     336                                        ./shared/Exceptions/exprintf.cpp\
     337                                        ./shared/Exp/exp.h\
     338                                        ./shared/Exp/IsInPoly.cpp\
     339                                        ./shared/Exp/IsInPolySerial.cpp\
     340                                        ./shared/Exp/DomainOutlineRead.cpp\
     341                                        ./shared/TriMesh/trimesh.h\
     342                                        ./shared/TriMesh/AssociateSegmentToElement.cpp\
     343                                        ./shared/TriMesh/GridInsideHole.cpp\
     344                                        ./shared/TriMesh/OrderSegments.cpp\
     345                                        ./shared/TriMesh/SplitMeshForRifts.cpp\
     346                                        ./shared/TriMesh/TriMeshUtils.cpp\
     347                                        ./shared/Sorting/binary_search.cpp\
     348                                        ./shared/Sorting/sorting.h\
     349                                        ./shared/Elements/elements.h\
     350                                        ./shared/Elements/ResolvePointers.cpp\
     351                                        ./shared/Elements/Paterson.cpp\
     352                                        ./shared/Elements/GetElementNodeData.cpp\
     353                                        ./toolkits/petsc\
     354                                        ./toolkits/petsc/patches\
     355                                        ./toolkits/petsc/patches/petscpatches.h\
     356                                        ./toolkits/petsc/patches/MatlabMatrixToPetscMatrix.cpp\
     357                                        ./toolkits/petsc/patches/MatlabVectorToPetscVector.cpp\
     358                                        ./toolkits/petsc/patches/PetscMatrixToMatlabMatrix.cpp\
     359                                        ./toolkits/petsc/patches/PetscVectorToMatlabVector.cpp\
     360                                        ./toolkits/petsc/patches/MatlabMatrixToDoubleMatrix.cpp\
     361                                        ./toolkits/petsc/patches/MatlabVectorToDoubleVector.cpp\
     362                                        ./toolkits/petsc/patches/PetscDetermineLocalSize.cpp\
     363                                        ./toolkits/petsc/patches/VecTranspose.cpp\
     364                                        ./toolkits/petsc/patches/VecToMPISerial.cpp\
     365                                        ./toolkits/petsc/patches/MatToSerial.cpp\
     366                                        ./toolkits/petsc/patches/VecMerge.cpp\
     367                                        ./toolkits/petsc/patches/NewVec.cpp\
     368                                        ./toolkits/petsc/patches/NewVecFromLocalSize.cpp\
     369                                        ./toolkits/petsc/patches/PetscOptionsDetermineSolverType.cpp\
     370                                        ./toolkits/petsc/patches/PetscOptionsInsertMultipleString.cpp\
     371                                        ./toolkits/petsc/patches/NewMat.cpp\
     372                                        ./toolkits/petsc/patches/VecFree.cpp\
     373                                        ./toolkits/petsc/patches/VecDuplicatePatch.cpp\
     374                                        ./toolkits/petsc/patches/KSPFree.cpp\
     375                                        ./toolkits/petsc/patches/ISFree.cpp\
     376                                        ./toolkits/petsc/patches/MatFree.cpp\
     377                                        ./toolkits/petsc/patches/GetOwnershipBoundariesFromRange.cpp\
     378                                        ./toolkits/petsc/patches/VecPartition.cpp\
     379                                        ./toolkits/petsc/patches/MatPartition.cpp\
     380                                        ./toolkits/petsc/patches/MatInvert.cpp\
     381                                        ./toolkits/petsc/patches/MatMultPatch.cpp\
     382                                        ./toolkits/petsc/petscincludes.h\
     383                                        ./toolkits/mpi/mpiincludes.h\
     384                                        ./toolkits/mpi/patches/mpipatches.h\
     385                                        ./toolkits/mpi/patches/MPI_Upperrow.cpp\
     386                                        ./toolkits/mpi/patches/MPI_Lowerrow.cpp\
     387                                        ./toolkits/mpi/patches/MPI_Boundariesfromrange.cpp\
     388                                        ./toolkits/metis/metisincludes.h\
     389                                        ./toolkits/triangle/triangleincludes.h\
     390                                        ./toolkits.h\
     391                                        ./io/io.h\
     392                                        ./io/FetchData.cpp\
     393                                        ./io/WriteData.cpp\
     394                                        ./io/WriteDataToDisk.cpp\
     395                                        ./io/SerialFetchData.cpp\
     396                                        ./io/SerialWriteData.cpp\
     397                                        ./io/ParallelFetchData.cpp\
     398                                        ./io/ParallelFetchInteger.cpp\
     399                                        ./io/ParallelFetchMat.cpp\
     400                                        ./io/ParallelFetchScalar.cpp\
     401                                        ./io/ParallelFetchString.cpp\
     402                                        ./io/ModelFetchData.cpp\
     403                                        ./io/WriteNodeSets.cpp\
     404                                        ./io/WriteParams.cpp\
     405                                        ./io/FetchNodeSets.cpp\
     406                                        ./io/FetchRifts.cpp\
     407                                        ./io/ParameterInputsInit.cpp\
     408                                        ./EnumDefinitions/EnumDefinitions.h\
     409                                        ./EnumDefinitions/EnumDefinitions.cpp\
     410                                        ./EnumDefinitions/AnalysisTypeAsEnum.cpp\
     411                                        ./ModelProcessorx/Model.h\
     412                                        ./ModelProcessorx/Model.cpp\
     413                                        ./ModelProcessorx/CreateDataSets.cpp\
    418414                                        ./ModelProcessorx/CreateParameters.cpp\
    419415                                        ./ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp\
    … …  
    510506bin_PROGRAMS =
    511507else
    512 bin_PROGRAMS = diagnostic.exe control.exe
    513 #thermalsteady.exe
     508bin_PROGRAMS = diagnostic.exe control.exe thermal.exe
    514509endif
    515510
    … …  
    522517control_exe_CXXFLAGS= -fPIC -D_PARALLEL_
    523518
    524 thermalsteady_exe_SOURCES = parallel/thermalsteady.cpp
    525 thermalsteady_exe_CXXFLAGS= -fPIC -D_PARALLEL_
     519thermal_exe_SOURCES = parallel/thermal.cpp
     520thermal_exe_CXXFLAGS= -fPIC -D_PARALLEL_

  • issm/trunk/src/c/Misfitx/Misfitx.cpp

    r1 r465  
    1414
    1515void Misfitx( double* pJ, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials,
    16                 double* u_g, double* u_g_obs, ParameterInputs* inputs,int analysis_type){
     16                double* u_g, double* u_g_obs, ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
    1717       
    1818        /*output: */
    … …  
    2424
    2525        /*Compute gradients: */
    26         elements->Misfit(&J, u_g,u_g_obs,inputs,analysis_type);
     26        elements->Misfit(&J, u_g,u_g_obs,inputs,analysis_type,sub_analysis_type);
    2727
    2828        /*Sum all J from all cpus of the cluster:*/

  • issm/trunk/src/c/Misfitx/Misfitx.h

    r1 r465  
    1010/* local prototypes: */
    1111void Misfitx( double* pJ, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials,
    12                 double* u_g, double* u_g_obs, ParameterInputs* inputs,int analysis_type);
     12                double* u_g, double* u_g_obs, ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
    1313
    1414#endif  /* _MISFITX_H */

  • issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp

    r209 r465  
    1313#include "../Model.h"
    1414
    15 void CreateParametersControl(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount){
     15void CreateParametersControl(DataSet** pparameters,Model* model,ConstDataHandle model_handle){
    1616       
    1717        int i;
    1818       
     19        DataSet* parameters=NULL;
    1920        Param*   param = NULL;
    2021        int      count;
    … …  
    2930        double* vy_obs=NULL;
    3031
    31         /*Get counter : */
    32         count=*pcount;
     32        /*Get parameters: */
     33        parameters=*pparameters;
     34        count=parameters->Size();
    3335       
    3436        /*control_type: */
    … …  
    132134
    133135        /*Assign output pointer: */
    134         *pcount=count;
     136        *pparameters=parameters;
    135137}
    136138

  • issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp

    r387 r465  
    2121        int      count=1;
    2222        int      analysis_type;
     23        int      sub_analysis_type;
    2324        int      numberofdofspernode;
    2425        int      dim;
    … …  
    2930        //Get analysis_type:
    3031        analysis_type=AnalysisTypeAsEnum(model->analysis_type);
     32        sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type);
    3133       
    3234        count++;
    3335        param= new Param(count,"analysis_type",INTEGER);
    3436        param->SetInteger(analysis_type);
     37        parameters->AddObject(param);
     38
     39        count++;
     40        param= new Param(count,"sub_analysis_type",INTEGER);
     41        param->SetInteger(sub_analysis_type);
    3542        parameters->AddObject(param);
    3643
    … …  
    187194
    188195        /*Deal with numberofdofspernode: */
    189         DistributeNumDofs(&numberofdofspernode,analysis_type);
     196        DistributeNumDofs(&numberofdofspernode,analysis_type,sub_analysis_type);
    190197
    191198        count++;
    … …  
    194201        parameters->AddObject(param);
    195202
    196         /*Deal with control parameters: */
    197         if (analysis_type==ControlAnalysisEnum()){
    198                 CreateParametersControl(parameters,model,model_handle,&count);
    199         }
    200         if (analysis_type==DiagnosticHorizAnalysisEnum()){
    201                 CreateParametersDiagnosticHoriz(parameters,model,model_handle,&count);
    202         }
    203203
    204204       

  • issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp

    r454 r465  
    3838
    3939        int         analysis_type;
     40        int         sub_analysis_type;
    4041       
    4142        /*output: */
    … …  
    154155        /*Get analysis_type: */
    155156        analysis_type=AnalysisTypeAsEnum(model->analysis_type);
     157        sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type);
    156158
    157159       
    … …  
    571573
    572574        /*Get number of dofs per node: */
    573         DistributeNumDofs(&node_numdofs,analysis_type);
     575        DistributeNumDofs(&node_numdofs,analysis_type,sub_analysis_type);
    574576
    575577        for (i=0;i<model->numberofnodes;i++){

  • issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateParametersDiagnosticHoriz.cpp

    r300 r465  
    1313#include "../Model.h"
    1414
    15 void CreateParametersDiagnosticHoriz(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount){
     15void CreateParametersDiagnosticHoriz(DataSet** pparameters,Model* model,ConstDataHandle model_handle){
    1616       
    1717        Param*   param = NULL;
     18        DataSet* parameters=NULL;
    1819        int      count;
    1920        int i;
    … …  
    2324        double* vz=NULL;
    2425
     26        /*recover parameters : */
     27        parameters=*pparameters;
     28
     29        count=parameters->Size();
     30
    2531        /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */
    2632        if (!model->ismacayealpattyn)return;
    27 
    28         /*Get counter : */
    29         count=*pcount;
    3033
    3134        /*Get vx and vy: */
    … …  
    6164       
    6265        /*Assign output pointer: */
    63         *pcount=count;
     66        *pparameters=parameters;
    6467}
    6568

  • issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp

    r367 r465  
    4040
    4141        int         analysis_type;
     42        int         sub_analysis_type;
    4243       
    4344        /*output: */
    … …  
    106107        /*Get analysis_type: */
    107108        analysis_type=AnalysisTypeAsEnum(model->analysis_type);
     109        sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type);
    108110
    109111        /*Width of elements: */
    … …  
    297299       
    298300        /*Get number of dofs per node: */
    299         DistributeNumDofs(&node_numdofs,analysis_type);
     301        DistributeNumDofs(&node_numdofs,analysis_type,sub_analysis_type);
    300302
    301303        for (i=0;i<model->numberofnodes;i++){

  • issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp

    r454 r465  
    3838
    3939        int         analysis_type;
     40        int         sub_analysis_type;
    4041       
    4142        /*output: */
    … …  
    134135        /*Get analysis_type: */
    135136        analysis_type=AnalysisTypeAsEnum(model->analysis_type);
     137        sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type);
    136138
    137139        /*Width of elements: */
    … …  
    428430
    429431        /*Get number of dofs per node: */
    430         DistributeNumDofs(&node_numdofs,analysis_type);
     432        DistributeNumDofs(&node_numdofs,analysis_type,sub_analysis_type);
    431433
    432434        for (i=0;i<model->numberofnodes;i++){

  • issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp

    r454 r465  
    3737
    3838        int         analysis_type;
     39        int         sub_analysis_type;
    3940       
    4041        /*output: */
    … …  
    143144        /*Get analysis_type: */
    144145        analysis_type=AnalysisTypeAsEnum(model->analysis_type);
     146        sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type);
    145147
    146148        /*Width of elements: */
    … …  
    354356       
    355357        /*Get number of dofs per node: */
    356         DistributeNumDofs(&node_numdofs,analysis_type);
     358        DistributeNumDofs(&node_numdofs,analysis_type,sub_analysis_type);
    357359
    358360        for (i=0;i<model->numberofnodes;i++){

  • issm/trunk/src/c/ModelProcessorx/Model.cpp

    r394 r465  
    3737        model->meshtype=NULL;
    3838        model->analysis_type=NULL;
     39        model->sub_analysis_type=NULL;
    3940        model->numberofelements=0;
    4041        model->numberofnodes=0;
    … …  
    243244        xfree((void**)&model->meshtype);
    244245        xfree((void**)&model->analysis_type);
     246        xfree((void**)&model->sub_analysis_type);
    245247       
    246248        if(model->numrifts){
    … …  
    289291        /*In ModelInit, we get all the data that is not difficult to get, and that is small: */
    290292        ModelFetchData((void**)&model->analysis_type,NULL,NULL,model_handle,"analysis_type","String",NULL);
     293        ModelFetchData((void**)&model->sub_analysis_type,NULL,NULL,model_handle,"sub_analysis_type","String",NULL);
    291294
    292295        ModelFetchData((void**)&model->meshtype,NULL,NULL,model_handle,"type","String",NULL);

  • issm/trunk/src/c/ModelProcessorx/Model.h

    r387 r465  
    1616        char*   meshtype;
    1717        char*   analysis_type;
     18        char*   sub_analysis_type;
    1819        char*   solverstring;
    1920
    … …  
    174175
    175176        /*Creation of fem datasets: general drivers*/
    176         void    CreateElementsNodesAndMaterials(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle);
    177         void    CreateConstraints(DataSet** pconstraints, Model* model,ConstDataHandle model_handle);
    178         void    CreateLoads(DataSet** ploads, Model* model, ConstDataHandle model_handle);
     177        void    CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,Model* model,ConstDataHandle model_handle);
    179178        void    CreateParameters(DataSet** pparameters,Model* model,ConstDataHandle model_handle);
    180179
    … …  
    186185        void    CreateConstraintsDiagnosticHoriz(DataSet** pconstraints,Model* model,ConstDataHandle model_handle);
    187186        void    CreateLoadsDiagnosticHoriz(DataSet** ploads, Model* model, ConstDataHandle model_handle);
    188         void    CreateParametersDiagnosticHoriz(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount);
     187        void    CreateParametersDiagnosticHoriz(DataSet** pparameters,Model* model,ConstDataHandle model_handle);
    189188
    190189        /*diagnostic vertical*/
    … …  
    192191        void    CreateConstraintsDiagnosticVert(DataSet** pconstraints,Model* model,ConstDataHandle model_handle);
    193192        void    CreateLoadsDiagnosticVert(DataSet** ploads, Model* model, ConstDataHandle model_handle);
    194         void    CreateParametersDiagnosticVert(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount);
    195193
    196194        /*diagnostic hutter*/
    … …  
    198196        void    CreateConstraintsDiagnosticHutter(DataSet** pconstraints,Model* model,ConstDataHandle model_handle);
    199197        void    CreateLoadsDiagnosticHutter(DataSet** ploads, Model* model, ConstDataHandle model_handle);
    200         void    CreateParametersDiagnosticHutter(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount);
    201198
    202199        /*diagnostic stokes*/
    … …  
    204201        void    CreateConstraintsDiagnosticStokes(DataSet** pconstraints,Model* model,ConstDataHandle model_handle);
    205202        void    CreateLoadsDiagnosticStokes(DataSet** ploads, Model* model, ConstDataHandle model_handle);
    206         void    CreateParametersDiagnosticStokes(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount);
    207203
    208204        /*slope compute*/
    … …  
    210206        void    CreateConstraintsSlopeCompute(DataSet** pconstraints,Model* model,ConstDataHandle model_handle);
    211207        void    CreateLoadsSlopeCompute(DataSet** ploads, Model* model, ConstDataHandle model_handle);
    212         void    CreateParametersSlopeCompute(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount);
    213208
    214209        /*control:*/
    215         void    CreateParametersControl(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount);
     210        void    CreateParametersControl(DataSet** pparameters,Model* model,ConstDataHandle model_handle);
    216211
    217212

  • issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp

    r454 r465  
    3636
    3737        int         analysis_type;
     38        int         sub_analysis_type;
    3839       
    3940        /*output: */
    … …  
    131132        /*Get analysis_type: */
    132133        analysis_type=AnalysisTypeAsEnum(model->analysis_type);
     134        sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type);
    133135
    134136        /*Width of elements: */
    … …  
    397399
    398400        /*Get number of dofs per node: */
    399         DistributeNumDofs(&node_numdofs,analysis_type);
     401        DistributeNumDofs(&node_numdofs,analysis_type,sub_analysis_type);
    400402
    401403        for (i=0;i<model->numberofnodes;i++){

  • issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.cpp

    r1 r465  
    1414
    1515void PenaltyConstraintsx(int* pconverged, int* pnum_unstable_constraints, DataSet* elements,DataSet* nodes,
    16                 DataSet* loads,DataSet* materials, ParameterInputs* inputs,int analysis_type){
     16                DataSet* loads,DataSet* materials, ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
    1717
    1818        int i;
    … …  
    3131         * management routine, otherwise, skip : */
    3232        if (loads->RiftIsPresent()){
    33                 loads->RiftConstraints(&num_unstable_constraints,inputs,analysis_type);
     33                loads->RiftConstraints(&num_unstable_constraints,inputs,analysis_type,sub_analysis_type);
    3434        }
    3535        else if(loads->MeltingIsPresent()){
    36                 loads->MeltingConstraints(&num_unstable_constraints,inputs,analysis_type);
     36                loads->MeltingConstraints(&num_unstable_constraints,inputs,analysis_type,sub_analysis_type);
    3737        }
    3838        else{

  • issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.h

    r1 r465  
    1111/* local prototypes: */
    1212void PenaltyConstraintsx(int* pconverged, int* pnum_unstable_constraints, DataSet* elements,DataSet* nodes,
    13                 DataSet* loads,DataSet* materials, ParameterInputs* inputs,int analysis_type);
     13                DataSet* loads,DataSet* materials, ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
    1414
    1515#endif  /* _PENALTYCONSTRAINTSX_H */

  • issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp

    r1 r465  
    1414
    1515void PenaltySystemMatricesx(Mat Kgg, Vec pg,DataSet* elements,DataSet* nodes,DataSet* loads,DataSet* materials,
    16                 int kflag,int pflag,ParameterInputs* inputs,int analysis_type){
     16                int kflag,int pflag,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
    1717       
    1818        int i;
    … …  
    3333
    3434        /*Add penalties to stiffnesses, from loads: */
    35         if(kflag)loads->PenaltyCreateKMatrix(Kgg,inputs,kmax,analysis_type);
    36         if(pflag)loads->PenaltyCreatePVector(pg,inputs,kmax,analysis_type);
     35        if(kflag)loads->PenaltyCreateKMatrix(Kgg,inputs,kmax,analysis_type,sub_analysis_type);
     36        if(pflag)loads->PenaltyCreatePVector(pg,inputs,kmax,analysis_type,sub_analysis_type);
    3737       
    3838        /*Assemble matrices: */

  • issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.h

    r1 r465  
    1111/* local prototypes: */
    1212void PenaltySystemMatricesx(Mat Kgg, Vec pg,DataSet* elements,DataSet* nodes,DataSet* loads,DataSet* materials,
    13                 int kflag,int pflag,ParameterInputs* inputs,int analysis_type);
     13                int kflag,int pflag,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
    1414
    1515#endif  /* _PENALTYSYSTEMMATRICESX_H */

  • issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp

    r304 r465  
    5454
    5555       
    56         if (   (analysis_type==ControlAnalysisEnum()) ||
    57                (analysis_type==DiagnosticHorizAnalysisEnum()) ||
    58                (analysis_type==DiagnosticVertAnalysisEnum()) ||
    59                (analysis_type==DiagnosticStokesAnalysisEnum())
    60            ){
     56        if (   (analysis_type==ControlAnalysisEnum()) ||  (analysis_type==DiagnosticAnalysisEnum())){
    6157
    6258                parameters->FindParam((void*)&vx,"vx");

  • issm/trunk/src/c/SystemMatricesx/SystemMatricesx.cpp

    r128 r465  
    1414
    1515void SystemMatricesx(Mat* pKgg, Vec* ppg,DataSet* elements,DataSet* nodes,DataSet* loads,DataSet* materials,
    16                 int kflag,int pflag,int connectivity,int numberofdofspernode,ParameterInputs* inputs,int analysis_type){
     16                int kflag,int pflag,int connectivity,int numberofdofspernode,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
    1717       
    1818        int i;
    … …  
    4040
    4141        /*Fill stiffness matrix and right hand side vector, from elements: */
    42         if(kflag)elements->CreateKMatrix(Kgg,inputs,analysis_type);
    43         if(pflag)elements->CreatePVector(pg,inputs,analysis_type);
     42        if(kflag)elements->CreateKMatrix(Kgg,inputs,analysis_type,sub_analysis_type);
     43        if(pflag)elements->CreatePVector(pg,inputs,analysis_type,sub_analysis_type);
    4444       
    4545        /*Fill stiffness matrix and right hand side vector, from loads: */
    46         if(kflag)loads->CreateKMatrix(Kgg,inputs,analysis_type);
    47         if(pflag)loads->CreatePVector(pg,inputs,analysis_type);
     46        if(kflag)loads->CreateKMatrix(Kgg,inputs,analysis_type,sub_analysis_type);
     47        if(pflag)loads->CreatePVector(pg,inputs,analysis_type,sub_analysis_type);
    4848
    4949        /*Assemble matrices: */

  • issm/trunk/src/c/SystemMatricesx/SystemMatricesx.h

    r1 r465  
    1111/* local prototypes: */
    1212void SystemMatricesx(Mat* pKgg, Vec* ppg,DataSet* elements,DataSet* nodes,DataSet* loads,DataSet* materials,
    13                 int kflag,int pflag,int connectivity,int numberofdofspernode,ParameterInputs* inputs,int analysis_type);
     13                int kflag,int pflag,int connectivity,int numberofdofspernode,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
    1414
    1515#endif  /* _SYSTEMMATRICESX_H */

  • issm/trunk/src/c/io/SerialFetchData.cpp

    r464 r465  
    2727        int      outdataset_size;
    2828
    29         double*  outmatrix=NULL;
    30         Mat      outpetscmatrix=NULL;
     29        double*  outmatrix;
     30        Mat      outpetscmatrix;
    3131        double*  outmatrix_workspace=NULL;;
    3232        int      outmatrix_rows,outmatrix_cols;
    3333        int      petsc;
    3434
    35         double*  outvector=NULL;
    36         Vec      outpetscvector=NULL;
     35        double*  outvector;
     36        Vec      outpetscvector;
    3737        double*  outvector_workspace=NULL;
    3838        int      outvector_rows;

  • issm/trunk/src/c/objects/Beam.cpp

    r350 r465  
    212212#define __FUNCT__ "Beam::CreateKMatrix"
    213213
    214 void  Beam::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){
     214void  Beam::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
    215215
    216216        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    217         if (analysis_type==DiagnosticHutterAnalysisEnum()) {
    218 
    219                 CreateKMatrixDiagnosticHutter( Kgg,inputs,analysis_type);
    220 
     217        if (analysis_type==DiagnosticAnalysisEnum()) {
     218       
     219                if (sub_analysis_type==HutterAnalysisEnum()) {
     220
     221                        CreateKMatrixDiagnosticHutter( Kgg,inputs,analysis_type,sub_analysis_type);
     222                }
     223                else
     224                        throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis_type: ",sub_analysis_type," not supported yet"));
    221225        }
    222226        else{
    … …  
    230234#define __FUNCT__ "Beam::CreateKMatrixDiagnosticHutter"
    231235
    232 void  Beam::CreateKMatrixDiagnosticHutter(Mat Kgg,void* vinputs,int analysis_type){
     236void  Beam::CreateKMatrixDiagnosticHutter(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
    233237       
    234238       
    … …  
    265269#undef __FUNCT__
    266270#define __FUNCT__ "Beam::CreatePVector"
    267 void  Beam::CreatePVector(Vec pg,void* inputs,int analysis_type){
     271void  Beam::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){
    268272       
    269273        /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
    270         if (analysis_type==DiagnosticHutterAnalysisEnum()) {
    271                 CreatePVectorDiagnosticHutter( pg,inputs,analysis_type);
     274        if (analysis_type==DiagnosticAnalysisEnum()) {
     275                if (sub_analysis_type==HutterAnalysisEnum()) {
     276                        CreatePVectorDiagnosticHutter( pg,inputs,analysis_type,sub_analysis_type);
     277                }
     278                else
     279                        throw ErrorException(__FUNCT__,exprintf("%s%i%s"," analysis ",analysis_type," not supported yet"));
    272280        }
    273281        else{
    … …  
    281289#define __FUNCT__ "Beam::CreatePVectorDiagnosticHutter"
    282290
    283 void Beam::CreatePVectorDiagnosticHutter( Vec pg, void* vinputs, int analysis_type){
     291void Beam::CreatePVectorDiagnosticHutter( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){
    284292
    285293        int i,j,k;
    … …  
    534542#undef __FUNCT__
    535543#define __FUNCT__ "Beam::Du"
    536 void  Beam::Du(_p_Vec*, double*, double*, void*, int){
     544void  Beam::Du(_p_Vec*, double*, double*, void*, int,int){
    537545        throw ErrorException(__FUNCT__," not supported yet!");
    538546}
    539547#undef __FUNCT__
    540548#define __FUNCT__ "Beam::Gradj"
    541 void  Beam::Gradj(_p_Vec*, double*, double*, void*, int, char*){
     549void  Beam::Gradj(_p_Vec*, double*, double*, void*, int, int,char*){
    542550        throw ErrorException(__FUNCT__," not supported yet!");
    543551}
    544552#undef __FUNCT__
    545553#define __FUNCT__ "Beam::GradjDrag"
    546 void  Beam::GradjDrag(_p_Vec*, double*, double*, void*, int){
     554void  Beam::GradjDrag(_p_Vec*, double*, double*, void*, int,int ){
    547555        throw ErrorException(__FUNCT__," not supported yet!");
    548556}
    549557#undef __FUNCT__
    550558#define __FUNCT__ "Beam::GradjB"
    551 void  Beam::GradjB(_p_Vec*, double*, double*, void*, int){
     559void  Beam::GradjB(_p_Vec*, double*, double*, void*, int, int){
    552560        throw ErrorException(__FUNCT__," not supported yet!");
    553561}
    554562#undef __FUNCT__
    555563#define __FUNCT__ "Beam::Misfit"
    556 double Beam::Misfit(double*, double*, void*, int){
     564double Beam::Misfit(double*, double*, void*, int,int ){
    557565        throw ErrorException(__FUNCT__," not supported yet!");
    558566}

  • issm/trunk/src/c/objects/Beam.h

    r308 r465  
    5656                int   MyRank();
    5757                void  Configure(void* loads,void* nodes,void* materials);
    58                 void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type);
    59                 void  CreatePVector(Vec pg, void* inputs, int analysis_type);
     58                void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
     59                void  CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type);
    6060                void  UpdateFromInputs(void* inputs);
    6161                void  GetDofList(int* doflist,int* pnumberofdofs);
    6262                void  GetDofList1(int* doflist);
    63                 void  CreateKMatrixDiagnosticHutter(Mat Kgg,void* inputs,int analysis_type);
     63                void  CreateKMatrixDiagnosticHutter(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
    6464                void  GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3);
    65                 void  CreatePVectorDiagnosticHutter(Vec pg,void* inputs,int analysis_type);
     65                void  CreatePVectorDiagnosticHutter(Vec pg,void* inputs,int analysis_type,int sub_analysis_type);
    6666                void* GetMatPar();
    6767
    … …  
    7878                void  GetBedList(double*);
    7979                void  GetThicknessList(double* thickness_list);
    80                 void  Du(_p_Vec*, double*, double*, void*, int);
    81                 void  Gradj(_p_Vec*, double*, double*, void*, int, char*);
    82                 void  GradjDrag(_p_Vec*, double*, double*, void*, int);
    83                 void  GradjB(_p_Vec*, double*, double*, void*, int);
    84                 double Misfit(double*, double*, void*, int);
     80                void  Du(_p_Vec*, double*, double*, void*, int,int);
     81                void  Gradj(_p_Vec*, double*, double*, void*, int, int,char*);
     82                void  GradjDrag(_p_Vec*, double*, double*, void*, int,int );
     83                void  GradjB(_p_Vec*, double*, double*, void*, int,int );
     84                double Misfit(double*, double*, void*, int,int );
    8585                void  GetNodalFunctions(double* l1l2, double gauss_coord);
    8686                void  GetParameterValue(double* pvalue, double* value_list,double gauss_coord);

  • issm/trunk/src/c/objects/Element.h

    r301 r465  
    2424                virtual void  Demarshall(char** pmarshalled_dataset)=0;
    2525                virtual void  Configure(void* loads,void* nodes,void* materials)=0;
    26                 virtual void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type)=0;
    27                 virtual void  CreatePVector(Vec pg, void* inputs, int analysis_type)=0;
     26                virtual void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type)=0;
     27                virtual void  CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type)=0;
    2828                virtual void  UpdateFromInputs(void* inputs)=0;
    2929                virtual void  GetNodes(void** nodes)=0;
    … …  
    3333                virtual void  GetThicknessList(double* thickness_list)=0;
    3434                virtual void  GetBedList(double* bed_list)=0;
    35                 virtual void  Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type)=0;
    36                 virtual void  Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,char* control_type)=0;
    37                 virtual void  GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type)=0;
    38                 virtual void  GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type)=0;
    39         virtual double Misfit(double* u_g,double* u_g_obs,void* inputs,int analysis_type)=0;
     35                virtual void  Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type)=0;
     36                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;
     37                virtual void  GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type)=0;
     38                virtual void  GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type)=0;
     39        virtual double Misfit(double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type)=0;
    4040        virtual void  ComputePressure(Vec p_g)=0;
    4141

  • issm/trunk/src/c/objects/Icefront.cpp

    r442 r465  
    199199        return my_rank;
    200200}
    201 void  Icefront::DistributeNumDofs(int* numdofspernode,int analysis_type){return;}
     201void  Icefront::DistributeNumDofs(int* numdofspernode,int analysis_type,int sub_analysis_type){return;}
    202202               
    203203#undef __FUNCT__
    … …  
    240240#define __FUNCT__ "Icefront::CreateKMatrix"
    241241
    242 void  Icefront::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){
     242void  Icefront::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
    243243
    244244        /*No stiffness loads applied, do nothing: */
    … …  
    250250#undef __FUNCT__
    251251#define __FUNCT__ "Icefront::CreatePVector"
    252 void  Icefront::CreatePVector(Vec pg, void* inputs, int analysis_type){
     252void  Icefront::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
    253253
    254254        /*Just branch to the correct element icefront vector generator, according to the type of analysis we are carrying out: */
    255         if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){
    256                 CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type);
    257         }
    258         else if ((analysis_type==DiagnosticStokesAnalysisEnum())){
    259 
    260                 CreatePVectorDiagnosticStokes( pg,inputs,analysis_type);
     255        if (analysis_type==ControlAnalysisEnum()){
     256               
     257                CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
     258
     259        }
     260        else if (analysis_type==DiagnosticAnalysisEnum()){
     261       
     262                if (sub_analysis_type==HorizAnalysisEnum()){
     263               
     264                        CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
     265
     266                }
     267                else if (sub_analysis_type==StokesAnalysisEnum()){
     268                       
     269                        CreatePVectorDiagnosticStokes( pg,inputs,analysis_type,sub_analysis_type);
     270
     271                }
     272                else throw ErrorException(__FUNCT__,exprintf("%s%i%s"," sub_analysis_type: ",sub_analysis_type," not supported yet"));
     273
    261274        }
    262275        else{
    … …  
    268281#undef __FUNCT__
    269282#define __FUNCT__ "Icefront::CreatePVectorDiagnosticHoriz"
    270 void Icefront::CreatePVectorDiagnosticHoriz( Vec pg, void* inputs, int analysis_type){
     283void Icefront::CreatePVectorDiagnosticHoriz( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
    271284
    272285        /*Branck on the type of icefront: */
    273286        if (strcmp(type,"segment")==0){
    274                 CreatePVectorDiagnosticHorizSegment(pg,inputs,analysis_type);
     287                CreatePVectorDiagnosticHorizSegment(pg,inputs,analysis_type,sub_analysis_type);
    275288        }
    276289        else{
    277                 CreatePVectorDiagnosticHorizQuad(pg,inputs,analysis_type);
     290                CreatePVectorDiagnosticHorizQuad(pg,inputs,analysis_type,sub_analysis_type);
    278291        }
    279292}       
    … …  
    282295#define __FUNCT__ "Icefront::CreatePVectorDiagnosticHorizSegment"
    283296
    284 void Icefront::CreatePVectorDiagnosticHorizSegment( Vec pg,void* vinputs, int analysis_type){
     297void Icefront::CreatePVectorDiagnosticHorizSegment( Vec pg,void* vinputs, int analysis_type,int sub_analysis_type){
    285298
    286299        int i,j;
    … …  
    395408#undef __FUNCT__
    396409#define __FUNCT__ "Icefont::CreatePVectorDiagnosticHorizQuad"
    397 void Icefront::CreatePVectorDiagnosticHorizQuad( Vec pg,void* vinputs, int analysis_type){
     410void Icefront::CreatePVectorDiagnosticHorizQuad( Vec pg,void* vinputs, int analysis_type,int sub_analysis_type){
    398411
    399412        int i,j;
    … …  
    557570#undef __FUNCT__
    558571#define __FUNCT__ "Icefont::CreatePVectorDiagnosticStokes"
    559 void Icefront::CreatePVectorDiagnosticStokes( Vec pg,void* vinputs, int analysis_type){
     572void Icefront::CreatePVectorDiagnosticStokes( Vec pg,void* vinputs, int analysis_type,int sub_analysis_type){
    560573
    561574        int i,j;
    … …  
    853866#undef __FUNCT__
    854867#define __FUNCT__ "Icefront::PenaltyCreateKMatrix"
    855 void  Icefront::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type){
     868void  Icefront::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
    856869        /*do nothing: */
    857870}
    … …  
    859872#undef __FUNCT__
    860873#define __FUNCT__ "Icefront::PenaltyCreatePVector"
    861 void  Icefront::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type){
     874void  Icefront::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
    862875        /*do nothing: */
    863876}

  • issm/trunk/src/c/objects/Icefront.h

    r416 r465  
    5555                int   GetId();
    5656                int   MyRank();
    57                 void  DistributeNumDofs(int* numdofspernode,int analysis_type);
     57                void  DistributeNumDofs(int* numdofspernode,int analysis_type,int sub_analysis_type);
    5858                void  Configure(void* elements,void* nodes,void* materials);
    59                 void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type);
     59                void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
    6060                void  UpdateFromInputs(void* inputs);
    6161               
    62                 void  CreatePVector(Vec pg, void* inputs, int analysis_type);
    63                 void  CreatePVectorDiagnosticHoriz( Vec pg, void* inputs, int analysis_type);
    64                 void  CreatePVectorDiagnosticHorizSegment( Vec pg,void* inputs, int analysis_type);
    65                 void  CreatePVectorDiagnosticHorizQuad( Vec pg,void* inputs, int analysis_type);
    66                 void  CreatePVectorDiagnosticStokes( Vec pg,void* inputs, int analysis_type);
     62                void  CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type);
     63                void  CreatePVectorDiagnosticHoriz( Vec pg, void* inputs, int analysis_type,int sub_analysis_type);
     64                void  CreatePVectorDiagnosticHorizSegment( Vec pg,void* inputs, int analysis_type,int sub_analysis_type);
     65                void  CreatePVectorDiagnosticHorizQuad( Vec pg,void* inputs, int analysis_type,int sub_analysis_type);
     66                void  CreatePVectorDiagnosticStokes( Vec pg,void* inputs, int analysis_type,int sub_analysis_type);
    6767                void  GetDofList(int* doflist,int* pnumberofdofs);
    6868                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);
    … …  
    7171                void  QuadPressureLoadStokes(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list,
    7272                                              double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list, int fill);
    73                 void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type);
    74                 void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type);
     73                void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
     74                void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
    7575                Object* copy();
    7676};

  • issm/trunk/src/c/objects/Load.h

    r246 r465  
    2424                virtual void  Demarshall(char** pmarshalled_dataset)=0;
    2525                virtual void  Configure(void* elements,void* nodes,void* materials)=0;
    26                 virtual void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type)=0;
    27                 virtual void  CreatePVector(Vec pg, void* inputs, int analysis_type)=0;
     26                virtual void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type)=0;
     27                virtual void  CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type)=0;
    2828                virtual void  UpdateFromInputs(void* inputs)=0;
    29                 virtual void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type)=0;
    30                 virtual void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type)=0;
     29                virtual void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type)=0;
     30                virtual void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type)=0;
    3131
    3232                int           Enum();

  • issm/trunk/src/c/objects/OptArgs.h

    r377 r465  
    2020        mxArray* n;
    2121        mxArray* analysis_type;
     22        mxArray* sub_analysis_type;
    2223
    2324};

  • issm/trunk/src/c/objects/Pengrid.cpp

    r461 r465  
    135135        return my_rank;
    136136}
    137 void  Pengrid::DistributeNumDofs(int* numdofpernode,int analysis_type){return;}
     137void  Pengrid::DistributeNumDofs(int* numdofpernode,int analysis_type,int sub_analysis_type){return;}
    138138
    139139#undef __FUNCT__
    … …  
    156156#define __FUNCT__ "Pengrid::CreateKMatrix"
    157157
    158 void  Pengrid::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){
     158void  Pengrid::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
    159159
    160160        /*No loads applied, do nothing: */
    … …  
    165165#undef __FUNCT__
    166166#define __FUNCT__ "Pengrid::CreatePVector"
    167 void  Pengrid::CreatePVector(Vec pg, void* inputs, int analysis_type){
     167void  Pengrid::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
    168168
    169169        /*No loads applied, do nothing: */
    … …  
    179179#undef __FUNCT__
    180180#define __FUNCT__ "Pengrid::PenaltyCreateKMatrix"
    181 void  Pengrid::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type){
    182 
    183         if ((analysis_type==DiagnosticStokesAnalysisEnum())){
    184 
    185                 PenaltyCreateKMatrixDiagnosticStokes( Kgg,inputs,kmax,analysis_type);
     181void  Pengrid::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
     182
     183        if ((analysis_type==DiagnosticAnalysisEnum()) && ((sub_analysis_type==StokesAnalysisEnum()))){
     184
     185                PenaltyCreateKMatrixDiagnosticStokes( Kgg,inputs,kmax,analysis_type,sub_analysis_type);
    186186
    187187        }
    188188        else{
    189                 throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
     189                throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet"));
    190190        }
    191191
    … …  
    194194#undef __FUNCT__
    195195#define __FUNCT__ "Pengrid::PenaltyCreateKMatrixDiagnosticStokes"
    196 void  Pengrid::PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,void* vinputs,double kmax,int analysis_type){
     196void  Pengrid::PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
    197197       
    198198        const int numgrids=1;

  • issm/trunk/src/c/objects/Pengrid.h

    r394 r465  
    3939                int   GetId();
    4040                int   MyRank();
    41                 void  DistributeNumDofs(int* numdofspernode,int analysis_type);
     41                void  DistributeNumDofs(int* numdofspernode,int analysis_type,int sub_analysis_type);
    4242                void  Configure(void* elements,void* nodes,void* materials);
    43                 void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type);
    44                 void  CreatePVector(Vec pg, void* inputs, int analysis_type);
     43                void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
     44                void  CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type);
    4545                void  UpdateFromInputs(void* inputs);
    46                 void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type);
    47                 void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type);
    48                 void  PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,void* inputs,double kmax,int analysis_type);
     46                void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
     47                void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
     48                void  PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
    4949                void  GetDofList(int* doflist,int* pnumberofdofspernode);
    5050                Object* copy();

  • issm/trunk/src/c/objects/Penpair.cpp

    r246 r465  
    182182        return my_rank;
    183183}
    184 void  Penpair::DistributeNumDofs(int* numdofspernode,int analysis_type){return;}
     184void  Penpair::DistributeNumDofs(int* numdofspernode,int analysis_type,int sub_analysis_type){return;}
    185185
    186186#undef __FUNCT__
    … …  
    208208#define __FUNCT__ "Penpair::CreateKMatrix"
    209209
    210 void  Penpair::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){
     210void  Penpair::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
    211211
    212212        /*No loads applied, do nothing: */
    … …  
    217217#undef __FUNCT__
    218218#define __FUNCT__ "Penpair::CreatePVector"
    219 void  Penpair::CreatePVector(Vec pg, void* inputs, int analysis_type){
     219void  Penpair::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
    220220
    221221        /*No loads applied, do nothing: */
    … …  
    231231#undef __FUNCT__
    232232#define __FUNCT__ "Penpair::PenaltyCreateKMatrix"
    233 void  Penpair::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type){
     233void  Penpair::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
    234234        if(numdofs==1){
    235235               
    … …  
    245245#undef __FUNCT__
    246246#define __FUNCT__ "Penpair::PenaltyCreatePVector"
    247 void  Penpair::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type){
     247void  Penpair::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
    248248        if(numdofs==1){
    249249               

  • issm/trunk/src/c/objects/Penpair.h

    r246 r465  
    5151                int   GetId();
    5252                int   MyRank();
    53                 void  DistributeNumDofs(int* numdofspernode,int analysis_type);
     53                void  DistributeNumDofs(int* numdofspernode,int analysis_type,int sub_analysis_type);
    5454                void  Configure(void* elements,void* nodes,void* materials);
    55                 void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type);
    56                 void  CreatePVector(Vec pg, void* inputs, int analysis_type);
     55                void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
     56                void  CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type);
    5757                void  UpdateFromInputs(void* inputs);
    58                 void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type);
    59                 void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type);
     58                void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
     59                void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
    6060                Object* copy();
    6161};

  • issm/trunk/src/c/objects/Penta.cpp

    r442 r465  
    284284#define __FUNCT__ "Penta::CreateKMatrix"
    285285
    286 void  Penta::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){
     286void  Penta::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
    287287
    288288        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    289         if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){
    290 
    291                 CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type);
    292 
    293         }
    294         else if ((analysis_type==DiagnosticVertAnalysisEnum())){
    295 
    296                 CreateKMatrixDiagnosticVert( Kgg,inputs,analysis_type);
    297 
    298         }
    299         else if ((analysis_type==DiagnosticStokesAnalysisEnum())){
    300 
    301                 CreateKMatrixDiagnosticStokes( Kgg,inputs,analysis_type);
    302 
    303         }
    304         else if (
    305                         (analysis_type==SurfaceSlopeComputeXAnalysisEnum()) ||
    306                         (analysis_type==SurfaceSlopeComputeYAnalysisEnum()) ||
    307                         (analysis_type==BedSlopeComputeXAnalysisEnum()) ||
    308                         (analysis_type==BedSlopeComputeYAnalysisEnum())
    309                         ){
    310                
    311                 CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type);
     289        if (analysis_type==ControlAnalysisEnum()){
     290
     291                CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type);
     292
     293        }
     294        else if (analysis_type==DiagnosticAnalysisEnum()){
     295       
     296                if (sub_analysis_type==HorizAnalysisEnum()){
     297               
     298                        CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type);
     299                }
     300                else if (sub_analysis_type==VertAnalysisEnum()){
     301               
     302                        CreateKMatrixDiagnosticVert( Kgg,inputs,analysis_type,sub_analysis_type);
     303                }
     304                else if (sub_analysis_type==StokesAnalysisEnum()){
     305               
     306                        CreateKMatrixDiagnosticStokes( Kgg,inputs,analysis_type,sub_analysis_type);
     307               
     308                }
     309                else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"));
     310        }
     311        else if (analysis_type==SlopeComputeAnalysisEnum()){
     312               
     313                CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type,sub_analysis_type);
    312314        }
    313315        else{
    … …  
    320322#undef __FUNCT__
    321323#define __FUNCT__ "Penta:CreateKMatrixDiagnosticHoriz"
    322 void Penta::CreateKMatrixDiagnosticHoriz( Mat Kgg, void* vinputs, int analysis_type){
     324void Penta::CreateKMatrixDiagnosticHoriz( Mat Kgg, void* vinputs, int analysis_type,int sub_analysis_type){
    323325
    324326
    … …  
    421423                 *and use its CreateKMatrix functionality to fill the global stiffness matrix: */
    422424                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    423                 tria->CreateKMatrix(Kgg,inputs, analysis_type);
     425                tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type);
    424426                delete tria;
    425427                return;
    … …  
    552554
    553555                        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    554                         tria->CreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type);
     556                        tria->CreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type,sub_analysis_type);
    555557                        delete tria;
    556558                }
    … …  
    573575#undef __FUNCT__
    574576#define __FUNCT__ "Penta:CreateKMatrixDiagnosticVert"
    575 void Penta::CreateKMatrixDiagnosticVert( Mat Kgg, void* vinputs, int analysis_type){
     577void Penta::CreateKMatrixDiagnosticVert( Mat Kgg, void* vinputs, int analysis_type,int sub_analysis_type){
    576578
    577579        /* local declarations */
    … …  
    623625        if(onsurface){
    624626                tria=(Tria*)SpawnTria(3,4,5); //nodes 3,4 and 5 are on the surface
    625                 tria->CreateKMatrixDiagnosticSurfaceVert(Kgg,inputs, analysis_type);
     627                tria->CreateKMatrixDiagnosticSurfaceVert(Kgg,inputs, analysis_type,sub_analysis_type);
    626628                delete tria;
    627629        }
    … …  
    710712#undef __FUNCT__
    711713#define __FUNCT__ "Penta:CreateKMatrixDiagnosticStokes"
    712 void Penta::CreateKMatrixDiagnosticStokes( Mat Kgg, void* vinputs, int analysis_type){
     714void Penta::CreateKMatrixDiagnosticStokes( Mat Kgg, void* vinputs, int analysis_type,int sub_analysis_type){
    713715
    714716        int i,j;
    … …  
    985987#undef __FUNCT__
    986988#define __FUNCT__ "Penta::CreatePVector"
    987 void  Penta::CreatePVector(Vec pg, void* inputs, int analysis_type){
     989void  Penta::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
    988990
    989991        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    990         if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){
    991 
    992                 CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type);
    993         }
    994         else if ((analysis_type==DiagnosticVertAnalysisEnum())){
    995 
    996                 CreatePVectorDiagnosticVert( pg,inputs,analysis_type);
    997         }
    998         else if ((analysis_type==DiagnosticStokesAnalysisEnum())){
    999 
    1000                 CreatePVectorDiagnosticStokes( pg,inputs,analysis_type);
    1001         }
    1002         else if (
    1003                         (analysis_type==SurfaceSlopeComputeXAnalysisEnum()) ||
    1004                         (analysis_type==SurfaceSlopeComputeYAnalysisEnum()) ||
    1005                         (analysis_type==BedSlopeComputeXAnalysisEnum()) ||
    1006                         (analysis_type==BedSlopeComputeYAnalysisEnum())
    1007                         ){
    1008                
    1009                 CreatePVectorSlopeCompute( pg,inputs,analysis_type);
     992        if (analysis_type==ControlAnalysisEnum()){
     993               
     994                CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
     995        }
     996        else if (analysis_type==DiagnosticAnalysisEnum()){
     997       
     998                if (sub_analysis_type==HorizAnalysisEnum()){
     999
     1000                        CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
     1001                }
     1002                else if (sub_analysis_type==VertAnalysisEnum()){
     1003               
     1004                        CreatePVectorDiagnosticVert( pg,inputs,analysis_type,sub_analysis_type);
     1005                }
     1006                else if (sub_analysis_type==StokesAnalysisEnum()){
     1007               
     1008                        CreatePVectorDiagnosticStokes( pg,inputs,analysis_type,sub_analysis_type);
     1009                }
     1010                else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"));
     1011        }
     1012        else if (analysis_type==SlopeComputeAnalysisEnum()){
     1013               
     1014                CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type);
    10101015        }
    10111016        else{
    … …  
    10911096#undef __FUNCT__
    10921097#define __FUNCT__ "Penta::Du"
    1093 void  Penta::Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type){
     1098void  Penta::Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type){
    10941099
    10951100        int i;
    … …  
    11221127#undef __FUNCT__
    11231128#define __FUNCT__ "Penta::Gradj"
    1124 void  Penta::Gradj(Vec grad_g,double* velocity,double* adjoint,void* inputs,int analysis_type,char* control_type){
     1129void  Penta::Gradj(Vec grad_g,double* velocity,double* adjoint,void* inputs,int analysis_type,int sub_analysis_type,char* control_type){
    11251130
    11261131        if (strcmp(control_type,"drag")==0){
    1127                 GradjDrag( grad_g,velocity,adjoint,inputs,analysis_type);
     1132                GradjDrag( grad_g,velocity,adjoint,inputs,analysis_type,sub_analysis_type);
    11281133        }
    11291134        else if (strcmp(control_type,"B")==0){
    1130                 GradjB( grad_g, velocity, adjoint, inputs,analysis_type);
     1135                GradjB( grad_g, velocity, adjoint, inputs,analysis_type,sub_analysis_type);
    11311136        }
    11321137        else throw ErrorException(__FUNCT__,exprintf("%s%s","control type not supported yet: ",control_type));
    … …  
    11351140#undef __FUNCT__
    11361141#define __FUNCT__ "Penta::GradjDrag"
    1137 void  Penta::GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type){
     1142void  Penta::GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type){
    11381143       
    11391144       
    … …  
    11471152               
    11481153                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    1149                 tria->GradjDrag( grad_g,u_g,lambda_g,inputs,analysis_type);
     1154                tria->GradjDrag( grad_g,u_g,lambda_g,inputs,analysis_type,sub_analysis_type);
    11501155                delete tria;
    11511156                return;
    … …  
    11551160#undef __FUNCT__
    11561161#define __FUNCT__ "Penta::GradjB"
    1157 void  Penta::GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type){
     1162void  Penta::GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type){
    11581163        throw ErrorException(__FUNCT__," not supported yet!");
    11591164}
    … …  
    11611166#undef __FUNCT__
    11621167#define __FUNCT__ "Penta::Misfit"
    1163 double Penta::Misfit(double* velocity,double* obs_velocity,void* inputs,int analysis_type){
     1168double Penta::Misfit(double* velocity,double* obs_velocity,void* inputs,int analysis_type,int sub_analysis_type){
    11641169       
    11651170        double J;
    … …  
    16411646#undef __FUNCT__
    16421647#define __FUNCT__ "Penta::CreatePVectorDiagnosticHoriz"
    1643 void Penta::CreatePVectorDiagnosticHoriz( Vec pg, void* vinputs,int analysis_type){
     1648void Penta::CreatePVectorDiagnosticHoriz( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type){
    16441649
    16451650        int i,j;
    … …  
    17071712                 *and use its CreatePVector functionality to return an elementary load vector: */
    17081713                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    1709                 tria->CreatePVector(pg,inputs, analysis_type);
     1714                tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type);
    17101715                delete tria;
    17111716                return;
    … …  
    19972002#undef __FUNCT__
    19982003#define __FUNCT__ "Penta:CreatePVectorDiagnosticVert"
    1999 void  Penta::CreatePVectorDiagnosticVert( Vec pg, void* vinputs,int analysis_type){
     2004void  Penta::CreatePVectorDiagnosticVert( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type){
    20002005
    20012006        int i;
    … …  
    20552060        if(onbed){
    20562061                tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 are on the bedrock
    2057                 tria->CreatePVectorDiagnosticBaseVert(pg,inputs, analysis_type);
     2062                tria->CreatePVectorDiagnosticBaseVert(pg,inputs, analysis_type,sub_analysis_type);
    20582063                delete tria;
    20592064        }
    … …  
    21762181#define __FUNCT__ "Penta::CreateKMatrixSlopeCompute"
    21772182
    2178 void  Penta::CreateKMatrixSlopeCompute(Mat Kgg,void* inputs,int analysis_type){
     2183void  Penta::CreateKMatrixSlopeCompute(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
    21792184
    21802185        /*Collapsed formulation: */
    … …  
    21862191        /*Spawn Tria element from the base of the Penta: */
    21872192        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    2188         tria->CreateKMatrix(Kgg,inputs, analysis_type);
     2193        tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type);
    21892194        delete tria;
    21902195        return;
    … …  
    21952200#define __FUNCT__ "Penta::CreatePVectorSlopeCompute"
    21962201
    2197 void Penta::CreatePVectorSlopeCompute( Vec pg, void* inputs, int analysis_type){
     2202void Penta::CreatePVectorSlopeCompute( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
    21982203       
    21992204        /*Collapsed formulation: */
    … …  
    22052210        /*Spawn Tria element from the base of the Penta: */
    22062211        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    2207         tria->CreatePVector(pg,inputs, analysis_type);
     2212        tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type);
    22082213        delete tria;
    22092214        return;
    … …  
    28262831#undef __FUNCT__
    28272832#define __FUNCT__ "Penta::CreatePVectorDiagnosticStokes"
    2828 void Penta::CreatePVectorDiagnosticStokes( Vec pg, void* vinputs,int analysis_type){
     2833void Penta::CreatePVectorDiagnosticStokes( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type){
    28292834
    28302835

  • issm/trunk/src/c/objects/Penta.h

    r394 r465  
    7373                int   MyRank();
    7474                void  Configure(void* loads,void* nodes,void* materials);
    75                 void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type);
    76                 void  CreateKMatrixDiagnosticHoriz( Mat Kgg, void* inputs, int analysis_type);
    77                 void  CreateKMatrixDiagnosticVert( Mat Kgg, void* inputs, int analysis_type);
    78                 void  CreatePVector(Vec pg, void* inputs, int analysis_type);
     75                void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
     76                void  CreateKMatrixDiagnosticHoriz( Mat Kgg, void* inputs, int analysis_type,int sub_analysis_type);
     77                void  CreateKMatrixDiagnosticVert( Mat Kgg, void* inputs, int analysis_type,int sub_analysis_type);
     78                void  CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type);
    7979                void  UpdateFromInputs(void* inputs);
    8080                void  GetDofList(int* doflist,int* pnumberofdofs);
    … …  
    8484                void  GetNodes(void** nodes);
    8585                int   GetOnBed();
    86                 void  Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type);
    87                 void  Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,char* control_type);
    88                 void  GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type);
    89                 void  GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type);
    90         double Misfit(double* u_g,double* u_g_obs,void* inputs,int analysis_type);
     86                void  Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type);
     87                void  Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type);
     88                void  GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type);
     89                void  GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type);
     90        double Misfit(double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type);
    9191               
    9292                void          GetThicknessList(double* thickness_list);
    … …  
    105105                void  GetNodalFunctionsDerivativesParams(double* dl1dl2dl3dl4dl5dl6,double* gauss_l1l2l3l4);
    106106                void  GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss_l1l2l3l4);
    107                 void  CreatePVectorDiagnosticHoriz( Vec pg, void* inputs,int analysis_type);
    108                 void  CreatePVectorDiagnosticVert( Vec pg, void* inputs,int analysis_type);
     107                void  CreatePVectorDiagnosticHoriz( Vec pg, void* inputs,int analysis_type,int sub_analysis_type);
     108                void  CreatePVectorDiagnosticVert( Vec pg, void* inputs,int analysis_type,int sub_analysis_type);
    109109                void  GetParameterValue(double* pvalue, double* v_list,double* gauss_l1l2l3l4);
    110110                void  GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_l1l2l3l4);
    … …  
    113113                void  SlopeExtrude(Vec sg,double* sg_serial);
    114114                void  ComputePressure(Vec p_g);
    115                 void  CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type);
    116                 void  CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type);
     115                void  CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
     116                void  CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);
    117117
    118                 void  CreateKMatrixDiagnosticStokes( Mat Kgg, void* vinputs, int analysis_type);
    119                 void  CreatePVectorDiagnosticStokes( Vec pg, void* vinputs,int analysis_type);
     118                void  CreateKMatrixDiagnosticStokes( Mat Kgg, void* vinputs, int analysis_type,int sub_analysis_type);
     119                void  CreatePVectorDiagnosticStokes( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type);
    120120                void  ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp);
    121121                void  GetMatrixInvert(double*  Ke_invert, double* Ke);

  • issm/trunk/src/c/objects/Sing.cpp

    r350 r465  
    193193#define __FUNCT__ "Sing::CreateKMatrix"
    194194
    195 void  Sing::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){
     195void  Sing::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
    196196
    197197        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    198         if (analysis_type==DiagnosticHutterAnalysisEnum()){
    199 
    200                 CreateKMatrixDiagnosticHutter( Kgg,inputs,analysis_type);
     198        if ((analysis_type==DiagnosticAnalysisEnum()) && (sub_analysis_type==HutterAnalysisEnum())){
     199
     200                CreateKMatrixDiagnosticHutter( Kgg,inputs,analysis_type,sub_analysis_type);
    201201
    202202        }
    203203        else{
    204                 throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
     204                throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s\n","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet"));
    205205        }
    206206
    … …  
    211211#define __FUNCT__ "Sing::CreateKMatrixDiagnosticHutter"
    212212
    213 void  Sing::CreateKMatrixDiagnosticHutter(Mat Kgg,void* vinputs,int analysis_type){
     213void  Sing::CreateKMatrixDiagnosticHutter(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
    214214       
    215215        const int numgrids=1;
    … …  
    234234#undef __FUNCT__
    235235#define __FUNCT__ "Sing::CreatePVector"
    236 void  Sing::CreatePVector(Vec pg,void* inputs,int analysis_type){
     236void  Sing::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){
    237237       
    238238        /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
    239         if (analysis_type==DiagnosticHutterAnalysisEnum()){
    240                 CreatePVectorDiagnosticHutter( pg,inputs,analysis_type);
     239        if ((analysis_type==DiagnosticAnalysisEnum()) && (sub_analysis_type==HutterAnalysisEnum())){
     240       
     241                        CreatePVectorDiagnosticHutter( pg,inputs,analysis_type,sub_analysis_type);
     242
    241243        }
    242244        else{
    … …  
    251253#define __FUNCT__ "Sing::CreatePVectorDiagnosticHutter"
    252254
    253 void Sing::CreatePVectorDiagnosticHutter( Vec pg, void* vinputs, int analysis_type){
     255void Sing::CreatePVectorDiagnosticHutter( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){
    254256       
    255257       
    … …  
    427429#undef __FUNCT__
    428430#define __FUNCT__ "Sing::Du"
    429 void  Sing::Du(_p_Vec*, double*, double*, void*, int){
     431void  Sing::Du(_p_Vec*, double*, double*, void*, int,int){
    430432        throw ErrorException(__FUNCT__," not supported yet!");
    431433}
    432434#undef __FUNCT__
    433435#define __FUNCT__ "Sing::Gradj"
    434 void  Sing::Gradj(_p_Vec*, double*, double*, void*, int, char*){
     436void  Sing::Gradj(_p_Vec*, double*, double*, void*, int, int ,char*){
    435437        throw ErrorException(__FUNCT__," not supported yet!");
    436438}
    437439#undef __FUNCT__
    438440#define __FUNCT__ "Sing::GradjDrag"
    439 void  Sing::GradjDrag(_p_Vec*, double*, double*, void*, int){
     441void  Sing::GradjDrag(_p_Vec*, double*, double*, void*, int,int){
    440442        throw ErrorException(__FUNCT__," not supported yet!");
    441443}
    442444#undef __FUNCT__
    443445#define __FUNCT__ "Sing::GradjB"
    444 void  Sing::GradjB(_p_Vec*, double*, double*, void*, int){
     446void  Sing::GradjB(_p_Vec*, double*, double*, void*, int,int){
    445447        throw ErrorException(__FUNCT__," not supported yet!");
    446448}
    447449#undef __FUNCT__
    448450#define __FUNCT__ "Sing::Misfit"
    449 double Sing::Misfit(double*, double*, void*, int){
     451double Sing::Misfit(double*, double*, void*, int,int){
    450452        throw ErrorException(__FUNCT__," not supported yet!");
    451453}

  • issm/trunk/src/c/objects/Sing.h

    r308 r465  
    5151                int   MyRank();
    5252                void  Configure(void* loads,void* nodes,void* materials);
    53                 void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type);
    54                 void  CreatePVector(Vec pg, void* inputs, int analysis_type);
     53                void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
     54                void  CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type);
    5555                void  UpdateFromInputs(void* inputs);
    5656                void  GetDofList(int* doflist,int* pnumberofdofs);
    5757                void  GetDofList1(int* doflist);
    58                 void  CreateKMatrixDiagnosticHutter(Mat Kgg,void* inputs,int analysis_type);
     58                void  CreateKMatrixDiagnosticHutter(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
    5959                void  GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3);
    60                 void  CreatePVectorDiagnosticHutter(Vec pg,void* inputs,int analysis_type);
     60                void  CreatePVectorDiagnosticHutter(Vec pg,void* inputs,int analysis_type,int sub_analysis_type);
    6161                void* GetMatPar();
    6262
    … …  
    7373                void  GetBedList(double*);
    7474                void  GetThicknessList(double* thickness_list);
    75                 void  Du(_p_Vec*, double*, double*, void*, int);
    76                 void  Gradj(_p_Vec*, double*, double*, void*, int, char*);
    77                 void  GradjDrag(_p_Vec*, double*, double*, void*, int);
    78                 void  GradjB(_p_Vec*, double*, double*, void*, int);
    79                 double Misfit(double*, double*, void*, int);
     75                void  Du(_p_Vec*, double*, double*, void*, int,int);
     76                void  Gradj(_p_Vec*, double*, double*, void*, int, int,char*);
     77                void  GradjDrag(_p_Vec*, double*, double*, void*, int,int);
     78                void  GradjB(_p_Vec*, double*, double*, void*, int,int);
     79                double Misfit(double*, double*, void*, int,int);
    8080
    8181

  • issm/trunk/src/c/objects/Tria.cpp

    r423 r465  
    266266#define __FUNCT__ "Tria::CreateKMatrix"
    267267
    268 void  Tria::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){
     268void  Tria::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
    269269
    270270        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    271         if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){
    272 
    273                 CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type);
    274 
    275         }
    276         else if (
    277                         (analysis_type==SurfaceSlopeComputeXAnalysisEnum()) ||
    278                         (analysis_type==SurfaceSlopeComputeYAnalysisEnum()) ||
    279                         (analysis_type==BedSlopeComputeXAnalysisEnum()) ||
    280                         (analysis_type==BedSlopeComputeYAnalysisEnum())
    281                         ){
    282                
    283                 CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type);
     271        if (analysis_type==ControlAnalysisEnum()){
     272               
     273                CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type);
     274        }
     275        else if (analysis_type==DiagnosticAnalysisEnum()){
     276       
     277                if (sub_analysis_type==HorizAnalysisEnum()){
     278
     279                        CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type);
     280                }
     281                else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"));
     282
     283        }
     284        else if (analysis_type==SlopeComputeAnalysisEnum()){
     285
     286                CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type,sub_analysis_type);
     287
    284288        }
    285289        else{
    … …  
    293297#define __FUNCT__ "Tria::CreateKMatrixDiagnosticHoriz"
    294298
    295 void  Tria::CreateKMatrixDiagnosticHoriz(Mat Kgg,void* vinputs,int analysis_type){
     299void  Tria::CreateKMatrixDiagnosticHoriz(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
    296300
    297301
    … …  
    477481        /*Do not forget to include friction: */
    478482        if(!shelf){
    479                 CreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type);
     483                CreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type,sub_analysis_type);
    480484        }
    481485
    … …  
    508512#undef __FUNCT__
    509513#define __FUNCT__ "Tria::CreateKMatrixDiagnosticHorizFriction"
    510 void  Tria::CreateKMatrixDiagnosticHorizFriction(Mat Kgg,void* vinputs,int analysis_type){
     514void  Tria::CreateKMatrixDiagnosticHorizFriction(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
    511515
    512516
    … …  
    677681#define __FUNCT__ "Tria::CreateKMatrixSlopeCompute"
    678682
    679 void  Tria::CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type){
     683void  Tria::CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
    680684
    681685        /* local declarations */
    … …  
    758762#undef __FUNCT__
    759763#define __FUNCT__ "Tria::CreatePVector"
    760 void  Tria::CreatePVector(Vec pg,void* inputs,int analysis_type){
     764void  Tria::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){
    761765       
    762766        /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
    763         if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){
    764                 CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type);
    765         }
    766         else if (
    767                         (analysis_type==SurfaceSlopeComputeXAnalysisEnum()) ||
    768                         (analysis_type==SurfaceSlopeComputeYAnalysisEnum()) ||
    769                         (analysis_type==BedSlopeComputeXAnalysisEnum()) ||
    770                         (analysis_type==BedSlopeComputeYAnalysisEnum())
    771                         ){
    772                
    773                 CreatePVectorSlopeCompute( pg,inputs,analysis_type);
     767        if (analysis_type==ControlAnalysisEnum()){
     768               
     769                CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
     770       
     771        }
     772        else if (analysis_type==DiagnosticAnalysisEnum()){
     773                if (sub_analysis_type==HorizAnalysisEnum()){
     774               
     775                        CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
     776               
     777                }
     778                else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"));
     779        }
     780        else if (analysis_type==SlopeComputeAnalysisEnum()){
     781               
     782                CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type);
    774783        }
    775784        else{
    … …  
    779788}
    780789
    781 
    782 
    783790#undef __FUNCT__
    784791#define __FUNCT__ "Tria::CreatePVectorDiagnosticHoriz"
    785792
    786 void Tria::CreatePVectorDiagnosticHoriz( Vec pg, void* vinputs, int analysis_type){
     793void Tria::CreatePVectorDiagnosticHoriz( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){
    787794
    788795        int             i,j;
    … …  
    952959#define __FUNCT__ "Tria::CreatePVectorSlopeCompute"
    953960
    954 void Tria::CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type){
     961void Tria::CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){
    955962
    956963        int             i,j;
    … …  
    992999        for(i=0;i<numdofs;i++) pe_g[i]=0.0;
    9931000
    994         if ( (analysis_type==SurfaceSlopeComputeXAnalysisEnum()) || (analysis_type==SurfaceSlopeComputeYAnalysisEnum())){
     1001        if ( (sub_analysis_type==SurfaceXAnalysisEnum()) || (sub_analysis_type==SurfaceYAnalysisEnum())){
    9951002                for(i=0;i<numdofs;i++) param[i]=s[i];
    9961003        }
    997         if ( (analysis_type==BedSlopeComputeXAnalysisEnum()) || (analysis_type==BedSlopeComputeYAnalysisEnum())){
     1004        if ( (sub_analysis_type==BedXAnalysisEnum()) || (sub_analysis_type==BedYAnalysisEnum())){
    9981005                for(i=0;i<numdofs;i++) param[i]=b[i];
    9991006        }
    … …  
    10201027
    10211028                /*Build pe_g_gaussian vector: */
    1022                 if ( (analysis_type==SurfaceSlopeComputeXAnalysisEnum()) || (analysis_type==BedSlopeComputeXAnalysisEnum())){
     1029                if ( (sub_analysis_type==SurfaceXAnalysisEnum()) || (sub_analysis_type==BedXAnalysisEnum())){
    10231030                        for(i=0;i<numdofs;i++) pe_g_gaussian[i]=Jdet*gauss_weight*slope[0]*l1l2l3[i];
    10241031                }
    1025                 if ( (analysis_type==SurfaceSlopeComputeYAnalysisEnum()) || (analysis_type==BedSlopeComputeYAnalysisEnum())){
     1032                if ( (sub_analysis_type==SurfaceYAnalysisEnum()) || (sub_analysis_type==BedYAnalysisEnum())){
    10261033                        for(i=0;i<numdofs;i++) pe_g_gaussian[i]=Jdet*gauss_weight*slope[1]*l1l2l3[i];
    10271034                }
    … …  
    15431550#undef __FUNCT__
    15441551#define __FUNCT__ "Tria::Du"
    1545 void Tria::Du(Vec du_g,double* velocity,double* obs_velocity,void* vinputs,int analysis_type){
     1552void Tria::Du(Vec du_g,double* velocity,double* obs_velocity,void* vinputs,int analysis_type,int sub_analysis_type){
    15461553
    15471554        int i;
    … …  
    17451752#undef __FUNCT__
    17461753#define __FUNCT__ "Tria::Gradj"
    1747 void  Tria::Gradj(Vec grad_g,double* velocity,double* adjoint,void* inputs,int analysis_type,char* control_type){
     1754void  Tria::Gradj(Vec grad_g,double* velocity,double* adjoint,void* inputs,int analysis_type,int sub_analysis_type,char* control_type){
    17481755
    17491756        if (strcmp(control_type,"drag")==0){
    1750                 GradjDrag( grad_g,velocity,adjoint,inputs,analysis_type);
     1757                GradjDrag( grad_g,velocity,adjoint,inputs,analysis_type,sub_analysis_type);
    17511758        }
    17521759        else if (strcmp(control_type,"B")==0){
    1753                 GradjB( grad_g, velocity, adjoint, inputs,analysis_type);
     1760                GradjB( grad_g, velocity, adjoint, inputs,analysis_type,sub_analysis_type);
    17541761        }
    17551762        else throw ErrorException(__FUNCT__,exprintf("%s%s","control type not supported yet: ",control_type));
    … …  
    17581765#undef __FUNCT__
    17591766#define __FUNCT__ "Tria::GradjDrag"
    1760 void  Tria::GradjDrag(Vec grad_g,double* velocity,double* adjoint,void* vinputs,int analysis_type){
     1767void  Tria::GradjDrag(Vec grad_g,double* velocity,double* adjoint,void* vinputs,int analysis_type,int sub_analysis_type){
    17611768
    17621769
    … …  
    19521959#undef __FUNCT__
    19531960#define __FUNCT__ "Tria::GradjB"
    1954 void  Tria::GradjB(Vec grad_g,double* velocity,double* adjoint,void* vinputs,int analysis_type){
     1961void  Tria::GradjB(Vec grad_g,double* velocity,double* adjoint,void* vinputs,int analysis_type,int sub_analysis_type){
    19551962
    19561963        int i;
    … …  
    20932100#undef __FUNCT__
    20942101#define __FUNCT__ "Tria::Misfit"
    2095 double Tria::Misfit(double* velocity,double* obs_velocity,void* vinputs,int analysis_type){
     2102double Tria::Misfit(double* velocity,double* obs_velocity,void* vinputs,int analysis_type,int sub_analysis_type){
    20962103
    20972104        int i;
    … …  
    22772284#undef __FUNCT__
    22782285#define __FUNCT__ "Tria::CreateKMatrixDiagnosticSurfaceVert"
    2279 void  Tria::CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,void* vinputs,int analysis_type){
     2286void  Tria::CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
    22802287
    22812288        int i,j;
    … …  
    24052412#undef __FUNCT__
    24062413#define __FUNCT__ "Tria::CreatePVectorDiagnosticBaseVert"
    2407 void  Tria::CreatePVectorDiagnosticBaseVert(Vec pg,void* vinputs,int analysis_type){
     2414void  Tria::CreatePVectorDiagnosticBaseVert(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type){
    24082415
    24092416        int             i,j;

  • issm/trunk/src/c/objects/Tria.h

    r386 r465  
    6565                int   MyRank();
    6666                void  Configure(void* loads,void* nodes,void* materials);
    67                 void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type);
    68                 void  CreatePVector(Vec pg, void* inputs, int analysis_type);
     67                void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
     68                void  CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type);
    6969                void  UpdateFromInputs(void* inputs);
    7070                void  GetDofList(int* doflist,int* pnumberofdofs);
    7171                void  GetDofList1(int* doflist);
    72                 void  CreateKMatrixDiagnosticHoriz(Mat Kgg,void* inputs,int analysis_type);
    73                 void  CreateKMatrixDiagnosticHorizFriction(Mat Kgg,void* inputs,int analysis_type);
    74                 void  CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,void* inputs,int analysis_type);
    75                 void  CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type);
     72                void  CreateKMatrixDiagnosticHoriz(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
     73                void  CreateKMatrixDiagnosticHorizFriction(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
     74                void  CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
     75                void  CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
    7676                void  GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3);
    7777                void  GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss_l1l2l3);
    … …  
    8787                void  GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss_l1l2l3);
    8888                void  GetJacobian(double* J, double* xyz_list,double* gauss_l1l2l3);
    89                 void  Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type);
    90                 void  Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,char* control_type);
    91                 void  GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type);
    92                 void  GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type);
    93         double Misfit(double* u_g,double* u_g_obs,void* inputs,int analysis_type);
     89                void  Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type);
     90                void  Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type);
     91                void  GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type);
     92                void  GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type);
     93        double Misfit(double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type);
    9494
    95                 void  CreatePVectorDiagnosticHoriz(Vec pg,void* inputs,int analysis_type);
    96                 void  CreatePVectorDiagnosticBaseVert(Vec pg,void* inputs,int analysis_type);
    97                 void  CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type);
     95                void  CreatePVectorDiagnosticHoriz(Vec pg,void* inputs,int analysis_type,int sub_analysis_type);
     96                void  CreatePVectorDiagnosticBaseVert(Vec pg,void* inputs,int analysis_type,int sub_analysis_type);
     97                void  CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);
    9898                void* GetMatPar();
    9999                int   GetShelf();

  • issm/trunk/src/c/objects/cielo/PenpairLoad.c

    r215 r465  
    689689#undef __FUNCT__
    690690#define __FUNCT__ "PenpairLoadPenaltyCreateKMatrix"
    691 int PenpairLoadPenaltyCreateKMatrix( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, double kmax, int analysis_type){
     691int PenpairLoadPenaltyCreateKMatrix( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, double kmax, int analysis_type,int sub_analysis_type){
    692692
    693693       
    … …  
    702702
    703703        if (this->penpair.numdofs==1){
    704                 PenpairLoadPenaltyCreateKMatrixOneDof( pKe_gg,vpthis,inputs,K_flag,kmax,analysis_type);
     704                PenpairLoadPenaltyCreateKMatrixOneDof( pKe_gg,vpthis,inputs,K_flag,kmax,analysis_type,sub_analysis_type);
    705705        }
    706706        else if (this->penpair.numdofs==2){
    707                 PenpairLoadPenaltyCreateKMatrixTwoDof( pKe_gg,vpthis,inputs,K_flag,kmax,analysis_type);
     707                PenpairLoadPenaltyCreateKMatrixTwoDof( pKe_gg,vpthis,inputs,K_flag,kmax,analysis_type,sub_analysis_type);
    708708        }
    709709        else{
    … …  
    719719#undef __FUNCT__
    720720#define __FUNCT__ "PenpairLoadPenaltyCreateKMatrixOneDof"
    721 int PenpairLoadPenaltyCreateKMatrixOneDof( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inptus, int K_flag, double kmax, int analysis_type){
     721int PenpairLoadPenaltyCreateKMatrixOneDof( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inptus, int K_flag, double kmax, int analysis_type,int sub_analysis_type){
    722722
    723723
    … …  
    769769#undef __FUNCT__
    770770#define __FUNCT__ "PenpairLoadPenaltyCreateKMatrixTwoDof"
    771 int PenpairLoadPenaltyCreateKMatrixTwoDof( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, double kmax, int analysis_type){
     771int PenpairLoadPenaltyCreateKMatrixTwoDof( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, double kmax, int analysis_type,int sub_analysis_type){
    772772
    773773
    … …  
    916916#define __FUNCT__ "PenpairLoadPenaltyCreatePVector"
    917917
    918 int PenpairLoadPenaltyCreatePVector( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, double kmax, int analysis_type){
     918int PenpairLoadPenaltyCreatePVector( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, double kmax, int analysis_type,int sub_analysis_type){
    919919
    920920        PenpairLoad* this   = NULL;
    … …  
    931931        }
    932932        else if (this->penpair.numdofs==2){
    933                 PenpairLoadPenaltyCreatePVectorTwoDof( ppe_g,vpthis,inputs,kmax,analysis_type);
     933                PenpairLoadPenaltyCreatePVectorTwoDof( ppe_g,vpthis,inputs,kmax,analysis_type,sub_analysis_type);
    934934        }
    935935        else{
    … …  
    947947#define __FUNCT__ "PenpairLoadPenaltyCreatePVectorTwoDof"
    948948
    949 int PenpairLoadPenaltyCreatePVectorTwoDof( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, double kmax, int analysis_type){
     949int PenpairLoadPenaltyCreatePVectorTwoDof( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, double kmax, int analysis_type,int sub_analysis_type){
    950950
    951951        int             i,j;
    … …  
    10921092        return noerr;
    10931093}
    1094 int PotentialUnstableConstraints(DataSet* lst,ParameterInputs* inputs,int analysis_type_enum);
    10951094
    10961095/*--------------------------------------------------
    … …  
    11001099#define __FUNCT__ "PenpairLoadPotentialUnstableConstraint"
    11011100
    1102 int PenpairLoadPotentialUnstableConstraint(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type){
     1101int PotentialUnstableConstraints(DataSet* lst,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
     1102int PenpairLoadPotentialUnstableConstraint(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type,int sub_analysis_type){
    11031103
    11041104        /* vpthis for polymorphic function compatibility */
    … …  
    11861186#define __FUNCT__ "PenpairLoadPenaltyPreConstrain"
    11871187
    1188 int PenpairLoadPenaltyPreConstrain(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type){
     1188int PenpairLoadPenaltyPreConstrain(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type,int sub_analysis_type){
    11891189
    11901190        /* vpthis for polymorphic function compatibility */
    … …  
    12901290#define __FUNCT__ "PenpairLoadPenaltyConstrain"
    12911291
    1292 int PenpairLoadPenaltyConstrain(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type){
     1292int PenpairLoadPenaltyConstrain(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type,int sub_analysis_type){
    12931293
    12941294        /* vpthis for polymorphic function compatibility */
    … …  
    14171417#define __FUNCT__ "PenpairLoadPenetration"
    14181418
    1419 int PenpairLoadPenetration(double* ppenetration, void* vpthis,ParameterInputs* inputs, int analysis_type){
     1419int PenpairLoadPenetration(double* ppenetration, void* vpthis,ParameterInputs* inputs, int analysis_type,int sub_analysis_type){
    14201420
    14211421        /* vpthis for polymorphic function compatibility */

  • issm/trunk/src/c/objects/cielo/PentaElement.c

    r216 r465  
    818818#undef __FUNCT__
    819819#define __FUNCT__ "PentaElementCreateKMatrix"
    820 int PentaElementCreateKMatrix( ElemMatrix* *pKe_gg, ElemMatrix** pKTe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, int KT_flag, int analysis_type){
     820int PentaElementCreateKMatrix( ElemMatrix* *pKe_gg, ElemMatrix** pKTe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, int KT_flag, int analysis_type,int sub_analysis_type){
    821821
    822822
    … …  
    25642564#define __FUNCT__ "PentaElementCreatePVector"
    25652565
    2566 int PentaElementCreatePVector( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, int analysis_type){
     2566int PentaElementCreatePVector( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, int analysis_type,int sub_analysis_type){
    25672567
    25682568        int noerr=1;
    … …  
    50025002#undef __FUNCT__
    50035003#define __FUNCT__ "PentaElementCreateDuVector"
    5004 int PentaElementCreateDuVector( ElemVector* *pdue_g, void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type){
     5004int PentaElementCreateDuVector( ElemVector* *pdue_g, void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type,int sub_analysis_type){
    50055005
    50065006        int noerr=1;
    … …  
    50775077
    50785078                /*Ok, now triaelement is correctly configured, call on its method CreateDuVector: */
    5079                 noerr=TriaElementCreateDuVector( pdue_g, (void*)triaelement, velocity, obs_velocity, inputs, analysis_type);
     5079                noerr=TriaElementCreateDuVector( pdue_g, (void*)triaelement, velocity, obs_velocity, inputs, analysis_type,sub_analysis_type);
    50805080
    50815081                /*Now delete tria and triaelement: */
    … …  
    50925092#undef __FUNCT__
    50935093#define __FUNCT__ "PentaElementCreateGradVectors"
    5094 int PentaElementCreateGradjVectors( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,char* control_type){
     5094int PentaElementCreateGradjVectors( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,int sub_analysis_type,char* control_type){
    50955095
    50965096        int noerr=1;
    50975097       
    50985098        if strcasecmp_eq(control_type,"drag"){
    5099                 noerr=PentaElementCreateGradjVectorsDrag( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type);
     5099                noerr=PentaElementCreateGradjVectorsDrag( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type,sub_analysis_type);
    51005100        }
    51015101        else if strcasecmp_eq(control_type,"B"){
    5102                 noerr=PentaElementCreateGradjVectorsB( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type);
     5102                noerr=PentaElementCreateGradjVectorsB( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type,sub_analysis_type);
    51035103        }
    51045104        else{
    … …  
    51145114#undef __FUNCT__
    51155115#define __FUNCT__ "PentaElementMisfit"
    5116 int PentaElementMisfit( double* pJelem,void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type){
     5116int PentaElementMisfit( double* pJelem,void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type,int sub_analysis_type){
    51175117
    51185118        int noerr=1;
    … …  
    51895189
    51905190                /*Ok, now triaelement is correctly configured, call on its method Misfit: */
    5191                 noerr=TriaElementMisfit(pJelem,(void*)triaelement, velocity, obs_velocity, inputs, analysis_type);
     5191                noerr=TriaElementMisfit(pJelem,(void*)triaelement, velocity, obs_velocity, inputs, analysis_type,sub_analysis_type);
    51925192
    51935193                /*Now delete tria and triaelement: */
    … …  
    52045204#undef __FUNCT__
    52055205#define __FUNCT__ "PentaElementCreateGradVectorsDrag"
    5206 int PentaElementCreateGradjVectorsDrag( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type){
     5206int PentaElementCreateGradjVectorsDrag( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
    52075207
    52085208        int noerr=1;
    … …  
    52795279
    52805280                /*Ok, now triaelement is correctly configured, call on its method CreateGradjVectorsdrag: */
    5281                 noerr=TriaElementCreateGradjVectorsDrag( pgradje_g, (void*)triaelement,velocity, adjoint, inputs,analysis_type);
     5281                noerr=TriaElementCreateGradjVectorsDrag( pgradje_g, (void*)triaelement,velocity, adjoint, inputs,analysis_type,sub_analysis_type);
    52825282
    52835283                /*Now delete tria and triaelement: */
    … …  
    52955295#undef __FUNCT__
    52965296#define __FUNCT__ "PentaElementCreateGradjVectorsB"
    5297 int PentaElementCreateGradjVectorsB( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type){
     5297int PentaElementCreateGradjVectorsB( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
    52985298
    52995299        printf("Not supported yet!\n");

  • issm/trunk/src/c/objects/cielo/TriaElement.c

    r216 r465  
    777777#undef __FUNCT__
    778778#define __FUNCT__ "TriaElementCreateKMatrix"
    779 int TriaElementCreateKMatrix( ElemMatrix* *pKe_gg, ElemMatrix** pKTe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, int KT_flag, int analysis_type){
     779int TriaElementCreateKMatrix( ElemMatrix* *pKe_gg, ElemMatrix** pKTe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, int KT_flag, int analysis_type,int sub_analysis_type){
    780780
    781781
    … …  
    800800#define __FUNCT__ "TriaElementCreatePVector"
    801801
    802 int TriaElementCreatePVector( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, int analysis_type){
     802int TriaElementCreatePVector( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, int analysis_type,int sub_analysis_type){
    803803
    804804        int noerr=1;
    … …  
    20442044#undef __FUNCT__
    20452045#define __FUNCT__ "TriaElementCreateDuVector"
    2046 int TriaElementCreateDuVector( ElemVector* *pdue_g, void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type){
     2046int TriaElementCreateDuVector( ElemVector* *pdue_g, void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type,int sub_analysis_type){
    20472047                                       
    20482048        int             i,j;
    … …  
    22632263#undef __FUNCT__
    22642264#define __FUNCT__ "TriaElementCreateGradVectors"
    2265 int TriaElementCreateGradjVectors( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,char* control_type){
     2265int TriaElementCreateGradjVectors( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,int sub_analysis_type,char* control_type){
    22662266
    22672267        int noerr=1;
    22682268       
    22692269        if strcasecmp_eq(control_type,"drag"){
    2270                 noerr=TriaElementCreateGradjVectorsDrag( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type);
     2270                noerr=TriaElementCreateGradjVectorsDrag( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type,sub_analysis_type);
    22712271        }
    22722272        else if strcasecmp_eq(control_type,"B"){
    2273                 noerr=TriaElementCreateGradjVectorsB( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type);
     2273                noerr=TriaElementCreateGradjVectorsB( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type,sub_analysis_type);
    22742274        }
    22752275        else{
    … …  
    22852285#undef __FUNCT__
    22862286#define __FUNCT__ "TriaElementCreateGradVectorsDrag"
    2287 int TriaElementCreateGradjVectorsDrag( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type){
     2287int TriaElementCreateGradjVectorsDrag( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
    22882288
    22892289        int             i,j;
    … …  
    25662566#undef __FUNCT__
    25672567#define __FUNCT__ "TriaElementCreateGradjVectorsB"
    2568 int TriaElementCreateGradjVectorsB( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type){
     2568int TriaElementCreateGradjVectorsB( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
    25692569
    25702570        int             i,j;
    … …  
    27702770#undef __FUNCT__
    27712771#define __FUNCT__ "TriaElementMisfit"
    2772 int TriaElementMisfit( double* pJelem,void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type){
     2772int TriaElementMisfit( double* pJelem,void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type,int sub_analysis_type){
    27732773
    27742774        int             i,j;

  • issm/trunk/src/c/parallel/CreateFemModel.cpp

    r464 r465  
    1010#include "../issm.h"
    1111
    12 void CreateFemModel(FemModel* femmodel,ConstDataHandle MODEL,char* analysis_type){
     12void CreateFemModel(FemModel* femmodel,ConstDataHandle MODEL,char* analysis_type,char* sub_analysis_type){
    1313
    1414        /*Model output: */
    … …  
    3737        _printf_("   specifying analysis\n");
    3838        model->analysis_type=(char*)xmalloc((strlen(analysis_type)+1)*sizeof(char)); strcpy(model->analysis_type,analysis_type);
     39        model->sub_analysis_type=(char*)xmalloc((strlen(sub_analysis_type)+1)*sizeof(char)); strcpy(model->sub_analysis_type,sub_analysis_type);
    3940
    40         _printf_("   create elements, nodes and materials:\n");
    41         CreateElementsNodesAndMaterials(&elements,&nodes,&materials,model,MODEL);
    42 
    43         _printf_("   create constraints: \n");
    44         CreateConstraints(&constraints,model,MODEL);
    45        
    46         _printf_("   create loads: \n");
    47         CreateLoads(&loads,model,MODEL);
    48 
    49         _printf_("   create parameters: \n");
    50         CreateParameters(&parameters,model,MODEL);
     41        _printf_("   create datasets:\n");
     42        CreateDataSets(&elements,&nodes,&materials,&constraints,&loads,&parameters,model,MODEL);
    5143
    5244        _printf_("   create degrees of freedom: \n");

  • issm/trunk/src/c/parallel/GradJCompute.cpp

    r434 r465  
    2020        /*intermediary: */
    2121        int analysis_type;
     22        int sub_analysis_type;
    2223        char* solverstring=NULL;
    2324        char* control_type=NULL;
    … …  
    4243        /*some parameters:*/
    4344        femmodel->parameters->FindParam((void*)&analysis_type,"analysis_type");
     45        femmodel->parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type");
    4446        femmodel->parameters->FindParam((void*)&solverstring,"solverstring");
    4547        femmodel->parameters->FindParam((void*)&control_type,"control_type");
    4648
    4749        _printf_("%s\n","      recover solution for this stiffness and right hand side:");
    48         diagnostic_core_nonlinear(&u_g,&K_ff0,&K_fs0,femmodel,inputs,analysis_type);
     50        diagnostic_core_nonlinear(&u_g,&K_ff0,&K_fs0,femmodel,inputs,analysis_type,sub_analysis_type);
    4951        VecToMPISerial(&u_g_double,u_g); VecFree(&u_g);
    5052
    5153        _printf_("%s\n","      buid Du, difference between observed velocity and model velocity:");
    52         Dux( &du_g, femmodel->elements,femmodel->nodes,femmodel->loads,femmodel->materials,u_g_double,u_g_obs, inputs,analysis_type);
     54        Dux( &du_g, femmodel->elements,femmodel->nodes,femmodel->loads,femmodel->materials,u_g_double,u_g_obs, inputs,analysis_type,sub_analysis_type);
    5355
    5456        _printf_("%s\n","      reduce adjoint load from g-set to f-set:");
    … …  
    6870
    6971        Gradjx( &grad_g, femmodel->elements,femmodel->nodes, femmodel->loads, femmodel->materials,
    70                 u_g_double,lambda_g_double, inputs,analysis_type,control_type);
     72                u_g_double,lambda_g_double, inputs,analysis_type,sub_analysis_type,control_type);
    7173       
    7274        /*Free ressources:*/

  • issm/trunk/src/c/parallel/control.cpp

    r434 r465  
    2626        char* lockname=NULL;
    2727        int   analysis_type;
     28        int   sub_analysis_type;
    2829
    2930        /*Intermediary: */
    … …  
    7879       
    7980        _printf_("read and create finite element model:\n");
    80         CreateFemModel(&femmodel,fid,"control");
     81        CreateFemModel(&femmodel,fid,"control","");
    8182
    8283        /*Recover parameters used throughout the solution:*/
    … …  
    9394        femmodel.parameters->FindParam((void*)&u_g_obs,"u_g_obs");
    9495        femmodel.parameters->FindParam((void*)&analysis_type,"analysis_type");
     96        femmodel.parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type");
    9597        gsize=femmodel.nodes->NumberOfDofs();
    9698
    … …  
    156158                        inputs->Add(control_type,p_g,2,numberofnodes);
    157159                        inputs->Add("fit",fit[n]);
    158                         diagnostic_core_nonlinear(&u_g,NULL,NULL,&femmodel,inputs,analysis_type);
     160                        diagnostic_core_nonlinear(&u_g,NULL,NULL,&femmodel,inputs,analysis_type,sub_analysis_type);
    159161                        OutputControl(u_g,p_g,J,nsteps,femmodel.partition,outputfilename,femmodel.nodesets);
    160162                        _printf_("%s\n","      done.");
    … …  
    175177        UpdateFromInputsx(femmodel.elements,femmodel.nodes,femmodel.loads, femmodel.materials,inputs);
    176178       
    177         diagnostic_core_nonlinear(&u_g,NULL,NULL,&femmodel,inputs,analysis_type);
     179        diagnostic_core_nonlinear(&u_g,NULL,NULL,&femmodel,inputs,analysis_type,sub_analysis_type);
    178180       
    179181        _printf_("%s\n","      saving final results...");

  • issm/trunk/src/c/parallel/diagnostic.cpp

    r458 r465  
    5959        _printf_("read and create finite element model:\n");
    6060        _printf_("\n   reading diagnostic horiz model data:\n");
    61         CreateFemModel(&femmodels[0],fid,"diagnostic_horiz");
     61        CreateFemModel(&femmodels[0],fid,"diagnostic","horiz");
    6262        _printf_("\n   reading diagnostic vert model data:\n");
    63         CreateFemModel(&femmodels[1],fid,"diagnostic_vert");
     63        CreateFemModel(&femmodels[1],fid,"diagnostic","vert");
    6464        _printf_("\n   reading diagnostic stokes model data:\n");
    65         CreateFemModel(&femmodels[2],fid,"diagnostic_stokes");
     65        CreateFemModel(&femmodels[2],fid,"diagnostic","stokes");
    6666        _printf_("\n   reading diagnostic hutter model data:\n");
    67         CreateFemModel(&femmodels[3],fid,"diagnostic_hutter");
     67        CreateFemModel(&femmodels[3],fid,"diagnostic","hutter");
    6868        _printf_("\n   reading surface and bed slope computation model data:\n");
    69         CreateFemModel(&femmodels[4],fid,"slope_compute");
     69        CreateFemModel(&femmodels[4],fid,"slope_compute","");
    7070
    7171        _printf_("initialize inputs:\n");

  • issm/trunk/src/c/parallel/diagnostic_core.cpp

    r434 r465  
    7474                       
    7575                if(debug)_printf_("%s\n","computing surface slope (x and y derivatives)...");
    76                 diagnostic_core_linear(&slopex,fem_sl,inputs,SurfaceSlopeComputeXAnalysisEnum());
    77                 diagnostic_core_linear(&slopey,fem_sl,inputs,SurfaceSlopeComputeYAnalysisEnum());
     76                diagnostic_core_linear(&slopex,fem_sl,inputs,SlopeComputeAnalysisEnum(),SurfaceXAnalysisEnum());
     77                diagnostic_core_linear(&slopey,fem_sl,inputs,SlopeComputeAnalysisEnum(),SurfaceYAnalysisEnum());
    7878
    7979                if (dim==3){
    … …  
    8989
    9090                if(debug)_printf_("%s\n"," computing hutter velocities...");
    91                 diagnostic_core_linear(&ug,fem_dhu,inputs,DiagnosticHutterAnalysisEnum());
     91                diagnostic_core_linear(&ug,fem_dhu,inputs,DiagnosticAnalysisEnum(),HutterAnalysisEnum());
    9292
    9393                if(debug)_printf_("%s\n"," computing pressure according to MacAyeal...");
    … …  
    106106               
    107107                if(debug)_printf_("%s\n"," computing horizontal velocities...");
    108                 diagnostic_core_nonlinear(&ug,NULL,NULL,fem_dh,inputs,DiagnosticHorizAnalysisEnum());
     108                diagnostic_core_nonlinear(&ug,NULL,NULL,fem_dh,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
    109109
    110110                if(debug)_printf_("%s\n"," computing pressure according to MacAyeal...");
    … …  
    121121                if(debug)_printf_("%s\n"," computing vertical velocities...");
    122122                inputs->Add("velocity",ug_horiz,numberofdofspernode_dh,numberofnodes);
    123                 diagnostic_core_linear(&ug_vert,fem_dv,inputs,DiagnosticVertAnalysisEnum());
     123                diagnostic_core_linear(&ug_vert,fem_dv,inputs,DiagnosticAnalysisEnum(),VertAnalysisEnum());
    124124
    125125                if(debug)_printf_("%s\n"," combining horizontal and vertical velocities...");
    … …  
    138138
    139139                        if(debug)_printf_("%s\n","computing bed slope (x and y derivatives)...");
    140                         diagnostic_core_linear(&slopex,fem_sl,inputs,BedSlopeComputeXAnalysisEnum());
    141                         diagnostic_core_linear(&slopey,fem_sl,inputs,BedSlopeComputeYAnalysisEnum());
     140                        diagnostic_core_linear(&slopex,fem_sl,inputs,SlopeComputeAnalysisEnum(),BedXAnalysisEnum());
     141                        diagnostic_core_linear(&slopey,fem_sl,inputs,SlopeComputeAnalysisEnum(),BedYAnalysisEnum());
    142142                        SlopeExtrudex( slopex, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials);
    143143                        SlopeExtrudex( slopey, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials);
    … …  
    160160                        if(debug)_printf_("%s\n"," computing stokes velocities and pressure ...");
    161161                        VecFree(&ug);
    162                         diagnostic_core_nonlinear(&ug,NULL,NULL,fem_ds,inputs,DiagnosticStokesAnalysisEnum());
     162                        diagnostic_core_nonlinear(&ug,NULL,NULL,fem_ds,inputs,DiagnosticAnalysisEnum(),StokesAnalysisEnum());
    163163               
    164164                        //decondition" pressure

  • issm/trunk/src/c/parallel/diagnostic_core_linear.cpp

    r434 r465  
    1111#include "../issm.h"
    1212
    13 void diagnostic_core_linear(Vec* pug,FemModel* fem,ParameterInputs* inputs,int analysis_type){
     13void diagnostic_core_linear(Vec* pug,FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
    1414
    1515        /*parameters:*/
    … …  
    4141        //*Generate system matrices
    4242        if (debug) _printf_("   Generating matrices\n");
    43         SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type);
     43        SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type,sub_analysis_type);
    4444
    4545        /*!Reduce matrix from g to f size:*/

  • issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp

    r464 r465  
    1111#include "../issm.h"
    1212
    13 void diagnostic_core_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, FemModel* fem,ParameterInputs* inputs,int analysis_type){
     13void diagnostic_core_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
    1414
    1515
    … …  
    8484                if (debug) _printf_("   Generating matrices\n");
    8585                //*Generate system matrices
    86                 SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type);
     86                SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type,sub_analysis_type);
    8787
    8888                if (debug) _printf_("   Generating penalty matrices\n");
    8989                //*Generate penalty system matrices
    90                 PenaltySystemMatricesx(Kgg, pg,fem->elements,fem->nodes,loads,fem->materials,kflag,pflag,inputs,analysis_type);
     90                PenaltySystemMatricesx(Kgg, pg,fem->elements,fem->nodes,loads,fem->materials,kflag,pflag,inputs,analysis_type,sub_analysis_type);
    9191
    9292                if (debug) _printf_("   reducing matrix from g to f set\n");
    … …  
    122122                inputs->Add("velocity",ug,numberofdofspernode,numberofnodes);
    123123               
    124                 PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,loads,fem->materials,inputs,analysis_type);
     124                PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,loads,fem->materials,inputs,analysis_type,sub_analysis_type);
    125125
    126126                //Figure out if convergence is reached.
    … …  
    167167                kflag=1; pflag=0; //stiffness generation only
    168168       
    169                 SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type);
     169                SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type,sub_analysis_type);
    170170       
    171171                Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets);

  • issm/trunk/src/c/parallel/objectivefunctionC.cpp

    r434 r465  
    3939        double* p_g_copy=NULL;
    4040        int     analysis_type;
     41        int     sub_analysis_type;
    4142        Vec     u_g=NULL;
    4243        double* u_g_double=NULL;
    … …  
    5960        femmodel->parameters->FindParam((void*)&fit,"fit");
    6061        femmodel->parameters->FindParam((void*)&analysis_type,"analysis_type");
     62        femmodel->parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type");
    6163        femmodel->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
    6264
    … …  
    7577
    7678        //Run diagnostic with updated parameters.
    77         diagnostic_core_nonlinear(&u_g,NULL,NULL,femmodel,inputs,analysis_type);
     79        diagnostic_core_nonlinear(&u_g,NULL,NULL,femmodel,inputs,analysis_type,sub_analysis_type);
    7880        VecToMPISerial(&u_g_double,u_g); VecFree(&u_g);
    7981
    … …  
    8183        inputs->Add("fit",fit[n]);
    8284        Misfitx( &J, femmodel->elements,femmodel->nodes, femmodel->loads, femmodel->materials,
    83                 u_g_double,u_g_obs, inputs,analysis_type);
     85                u_g_double,u_g_obs, inputs,analysis_type,sub_analysis_type);
    8486
    8587

  • issm/trunk/src/c/parallel/parallel.h

    r434 r465  
    1212
    1313void diagnostic_core(Vec* pug, Vec* ppg,FemModel* fems, ParameterInputs* inputs);
    14 void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* fem,ParameterInputs* inputs,int analysis_type);
    15 void diagnostic_core_linear(Vec* ppg,FemModel* fem,ParameterInputs* inputs,int  analysis_type);
     14void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
     15void diagnostic_core_linear(Vec* ppg,FemModel* fem,ParameterInputs* inputs,int  analysis_type,int sub_analysis_type);
    1616int cielothermal_core(Vec** pt_g,ParameterInputs* inputs,FemModel* femmodel);
    1717
    … …  
    2929//int ParameterUpdate(double* search_vector,int step, WorkspaceParams* workspaceparams,BatchParams* batchparams);
    3030void OutputDiagnostic(Vec u_g,Vec p_g, FemModel* femmodels,char* filename);
    31 int OutputThermal(Vec* t_g,Vec* tpartition,char* filename,char* analysis_type);
     31int OutputThermal(Vec* t_g,Vec* tpartition,char* filename,char* analysis_type,int sub_analysis_type);
    3232void OutputControl(Vec u_g,double* p_g, double* J, int nsteps, Vec partition,char* filename,NodeSets* nodesets);
    3333void WriteLockFile(char* filename);
    3434
    35 void CreateFemModel(FemModel* femmodel,ConstDataHandle MODEL,char* analysis_type);
     35void CreateFemModel(FemModel* femmodel,ConstDataHandle MODEL,char* analysis_type,char* sub_analysis_type);
    3636//int BatchDebug(Mat* Kgg,Vec* pg,FemModel* femmodel,char* filename);
    3737

  • issm/trunk/src/c/parallel/thermal.cpp

    r386 r465  
    1 /*
    2  * cielothermalsteady.c:
    3  */
     1/*!\file:  thermal.cpp
     2 * \brief: thermal solution
     3 */ 
    44
    5 #include "../include/cielo.h"
    6 #include "../modules.h"
     5#include "../issm.h"
    76#include "./parallel.h"
    87
    98#undef __FUNCT__
    10 #define __FUNCT__ "cielothermalsteady"
     9#define __FUNCT__ "thermal"
     10
     11#ifdef HAVE_CONFIG_H
     12        #include "config.h"
     13#else
     14#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
     15#endif
    1116
    1217int main(int argc,char* *argv){
    13        
     18
    1419        /*I/O: */
    1520        FILE* fid=NULL;
    … …  
    1722        char* outputfilename=NULL;
    1823        char* lockname=NULL;
    19         char* analysis_type_t="thermalsteady";
    20         char* analysis_type_m="melting";
     24        int   numberofnodes;
    2125
    22         /*Intermediary: */
     26        /*Fem models : */
     27        FemModel femmodels[2];
     28
     29        /*initial velocity and pressure: */
     30        Vec u_g=NULL;
     31        Vec p_g=NULL;
     32
    2333       
    24         #if defined(_PARALLEL_) && defined(_HAVE_PETSC_)
    25         FemModel femmodel_t;
    26         FemModel femmodel_m;    Vec* t_g=NULL;
    27         #endif
     34        /*solution vectors: */
     35        Vec t_g=NULL;
     36        Vec m_g=NULL;
     37
    2838        ParameterInputs* inputs=NULL;
     39        int waitonlock=0;
     40       
     41        MODULEBOOT();
    2942
    3043        #if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_))
    31                 _printf_("%s%s\n",__FUNCT__," error message: parallel executable was compiled without support of parallel libraries!");
    32                 return 1;
    33         #else
     44        throw ErrorException(__FUNCT__," parallel executable was compiled without support of parallel libraries!");
     45        #endif
    3446
    35                 /*Initialize MPI environment: */
    36                 PetscInitialize(&argc,&argv,(char *)0,""); 
     47        PetscInitialize(&argc,&argv,(char *)0,""); 
    3748
    38                 /*Size and rank: */
    39                 MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); 
    40                 MPI_Comm_size(MPI_COMM_WORLD,&num_procs);
     49        /*Size and rank: */
     50        MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); 
     51        MPI_Comm_size(MPI_COMM_WORLD,&num_procs);
    4152
    42                 /*Some checks on size of cluster*/
    43                 if (num_procs<=1){
    44                         _printf_("\nSize of MPI COMM WORLD is 1, needs to be at least 2. Include more nodes\n");
    45                         PetscFinalize();
    46                         return 0;
    47                 }
     53        _printf_("recover , input file name and output file name:\n");
     54        inputfilename=argv[2];
     55        outputfilename=argv[3];
     56        lockname=argv[4];
    4857
    49                 /*Recover dbdir, input file name and output file name: */
    50                 dbdir=argv[1];
    51                 inputfilename=argv[2];
    52                 outputfilename=argv[3];
    53                 lockname=argv[4];
     58        /*Open handle to data on disk: */
     59        fid=fopen(inputfilename,"rb");
     60        if(fid==NULL) throw ErrorException(__FUNCT__,exprintf("%s%s%s","could not open file ",inputfilename," for binary reading"));
    5461
     62        _printf_("read and create finite element model:\n");
     63        CreateFemModel(&femmodels[0],fid,"thermal","");
     64        CreateFemModel(&femmodels[1],fid,"melting","");
     65
     66        _printf_("initialize inputs:\n");
     67        femmodels[0].parameters->FindParam((void*)&u_g,"u_g");
     68        femmodels[0].parameters->FindParam((void*)&p_g,"p_g");
     69        femmodels[0].parameters->FindParam((void*)&numberofnodes,"numberofnodes");
     70
     71        inputs=new ParameterInputs;
     72        //inputs->Add("velocity",u_g_initial,3,numberofnodes);
     73
     74        //erase velocities:
     75        //param=(Param*)femmodels[0].parameters->FindParamObject("u_g");
     76        //femmodels[0].parameters->DeleteObject((Object*)param);
     77
     78        _printf_("call computational core:\n");
     79        diagnostic_core(&u_g,&p_g,femmodels,inputs);
     80
     81        _printf_("write results to disk:\n");
     82        OutputDiagnostic(u_g,p_g,&femmodels[0],outputfilename);
     83
     84        _printf_("write lock file:\n");
     85        femmodels[0].parameters->FindParam((void*)&waitonlock,"waitonlock");
     86        if (waitonlock){
     87                WriteLockFile(lockname);
     88        }
    5589               
    56                
    57                 /*Read and create thermal finite element model: */
    58                 fid=fopen(inputfilename,"rb");
    59                 if(fid==NULL){
    60                         _printf_("%s%s\n",__FUNCT__," error message: could not open file ",inputfilename," for binary reading");
    61                         return 0;
    62                 }
    63                 if(!CreateFemModel(&femmodel_t,fid,analysis_type_t)){
    64                         _printf_("%s\n",__FUNCT__," error message: could not read melting finite element model!\n");
    65                         return 0;
    66                 }
    67                 fclose(fid);
     90        _printf_("closing MPI and Petsc\n");
     91        MPI_Barrier(MPI_COMM_WORLD);
    6892
    69                 /*Read and create melting finite element model: */
    70                 fid=fopen(inputfilename,"rb");
    71                 if(!CreateFemModel(&femmodel_m,fid,analysis_type_m)){
    72                         _printf_("%s\n",__FUNCT__," error message: could not read thermal finite element model!\n");
    73                         return 0;
    74                 }
    75                 fclose(fid);
     93        /*Close MPI libraries: */
     94        PetscFinalize();
    7695
    77                 /*Initialize inputs: */
    78                 inputs=NewParameterInputs();
    79                        
    80                 ParameterInputsAddFromMat(inputs,WorkspaceParamsGetVelocity(femmodel_t.workspaceparams),femmodel_t.workspaceparams->gsize,"velocity");
    81                 ParameterInputsAddFromMat(inputs,WorkspaceParamsGetPressure(femmodel_t.workspaceparams),femmodel_t.workspaceparams->gsize,"pressure");
    8296
    83                 ParameterInputsAddFromMat(inputs,&femmodel_t.batchparams->dt,1,"dt");
    84                
    85                 /*Call computational core: */
    86                 if(!cielothermal_core(&t_g,inputs,&femmodel_t)){
    87                         _printf_("%s\n",__FUNCT__," error message: could not run computational core!\n");
    88                         return 0;
    89                 }
    90                
    91                 /*Write results to disk: */
    92                 OutputThermal(t_g,femmodel_t.tpartition,outputfilename,analysis_type_t);
    93 
    94                 /*Write lock file if requested: */
    95                 if (femmodel_t.batchparams->waitonlock){
    96                         WriteLockFile(lockname);
    97                 }
    98                                        
    99                 cleanup_and_return:
    100                 _printf_("done.\n");
    101                
    102                 /*Synchronize everyone before exiting: */
    103                 MPI_Barrier(MPI_COMM_WORLD);
    104 
    105                 /*Close MPI libraries: */
    106                 PetscFinalize();
    107                
    108                 return 0; //unix success return;
    109 
    110         #endif
     97        /*end module: */
     98        MODULEEND();
     99       
     100        return 0; //unix success return;
    111101}
     102       

  • issm/trunk/src/c/shared/Dofs/dofs.h

    r434 r465  
    77
    88double* dofsetgen(int numdofs,int* doflist,int dofspernode,int totaldofs);
     9int DistributeNumDofs(int *pnumdofs,int analysis_type,int sub_analysis_type);
    910
    1011

  • issm/trunk/src/c/shared/Numerics/OptFunc.cpp

    r377 r465  
    2323        double J;
    2424
    25         mxArray*   inputs[8];
     25        mxArray*   inputs[9];
    2626        mxArray*   psearch_scalar=NULL;
    2727        mxArray*   mxJ=NULL;
    … …  
    3030        psearch_scalar=mxCreateDoubleScalar(scalar);
    3131        inputs[0]=psearch_scalar;
    32        
    3332        inputs[1]=optargs->m;
    3433        inputs[2]=optargs->inputs;
    … …  
    3837        inputs[6]=optargs->n;
    3938        inputs[7]=optargs->analysis_type;
     39        inputs[8]=optargs->sub_analysis_type;
    4040
    41         mexCallMATLAB( 1, &mxJ, 8, (mxArray**)inputs, optargs->function_name);
     41        mexCallMATLAB( 1, &mxJ, 9, (mxArray**)inputs, optargs->function_name);
    4242
    4343        /*extract misfit from mxArray*/

  • issm/trunk/src/c/shared/Numerics/numerics.h

    r117 r465  
    1515void cross(double* result,double* vector1,double* vector2);
    1616double norm(double* vector);
    17 int DistributeNumDofs(int *pnumdofs,int analysis_type);
    1817
    1918#endif //ifndef _NUMERICS_H_

  • issm/trunk/src/m/classes/@model/model.m

    r387 r465  
    255255        md.solverstring='';
    256256        md.analysis_type='';
     257        md.sub_analysis_type='';
    257258
    258259        %management of large models

  • issm/trunk/src/m/classes/public/BuildQueueingScript.m

    r1 r465  
    1 function BuildQueueingScript(md,solutiontype,executionpath,codepath)
     1function BuildQueueingScript(md,executionpath,codepath)
    22%BUILDQUEUEINGSCRIPT -
    33%
    44%   Usage:
    5 %      BuildQueueingScript(md,solutiontype,executionpath,codepath)
     5%      BuildQueueingScript(md,executionpath,codepath)
    66
    77disp('building queueing script');
    … …  
    1111if exist(function_name,'file'),
    1212        %Call this function:
    13         eval([function_name '(md,solutiontype,executionpath,codepath);']);
     13        eval([function_name '(md,executionpath,codepath);']);
    1414else
    1515        %Call the generic BuildQueueingScript:
    16         BuildQueueingScriptGeneric(md,solutiontype,executionpath,codepath);
     16        BuildQueueingScriptGeneric(md,executionpath,codepath);
    1717end

  • issm/trunk/src/m/classes/public/BuildQueueingScriptGeneric.m

    r1 r465  
    1 function BuildQueueingScriptGeneric(md,solutiontype,executionpath,codepath)
     1function BuildQueueingScriptGeneric(md,executionpath,codepath)
    22%BUILDQUEUEINGSCRIPTGENERIC - ...
    33%
    … …  
    1717fprintf(fid,'mpirun -np %i ',md.np);
    1818
    19 if strcmpi(solutiontype,'diagnostic_horiz') |  strcmpi(solutiontype,'diagnostic'),
     19if strcmpi(md.analysis_type,'diagnostic'),
    2020        fprintf(fid,'%s/diagnostic.exe',codepath);
    21 elseif strcmpi(solutiontype,'control'),
     21elseif strcmpi(md.analysis_type,'control'),
    2222        fprintf(fid,'%s/control.exe',codepath);
    23 elseif strcmpi(solutiontype,'thermalsteady'),
    24         fprintf(fid,'%s/thermalsteady.exe',codepath);
     23elseif strcmpi(md.analysis_type,'thermal'),
     24        fprintf(fid,'%s/thermal.exe',codepath);
    2525else
    2626        error('BuildQueueingScriptGeneric error message: unsupported solution type!');
    2727end
    2828
    29 
    3029fprintf(fid,' %s %s.bin %s.outbin %s.lock 2> %s.errlog >%s.outlog & ',executionpath,md.name,md.name,md.name,md.name,md.name);
    3130fclose(fid);

  • issm/trunk/src/m/classes/public/BuildQueueingScriptcosmos.m

    r1 r465  
    22%BUILDQUEUEINGSCRIPTCOSMOS - build queueing script for batch system when running parallel
    33%
    4 %   solutiontype is 'diagnostic','prognostic','transient','thermalsteady','thermaltransient','parameters' or 'control'
     4%   solutiontype is 'diagnostic','prognostic','transient','thermal','parameters' or 'control'
    55%   and varargin is an optional package name ('cielo', 'ice' or 'macayeal')
    66%
    … …  
    3838        fprintf(fid,'mpirun.lsf /home/larour/bin/alloc_cleanup.exe\n');
    3939end
    40 if strcmpi(solutiontype,'diagnostic_horiz') |  strcmpi(solutiontype,'diagnostic'),
     40if strcmpi(solutiontype,'diagnostic'),
    4141        fprintf(fid,'mpirun.lsf %s/diagnostic.exe %s %s.bin %s.outbin %s.lock',codepath,executionpath,md.name,md.name,md.name);
    4242elseif strcmpi(solutiontype,'control') ,
    4343        fprintf(fid,'mpirun.lsf %s/control.exe %s %s.bin %s.outbin %s.lock',codepath,executionpath,md.name,md.name,md.name);
    44 elseif strcmpi(solutiontype,'thermalsteady') ,
    45         fprintf(fid,'mpirun.lsf %s/thermalsteady.exe %s %s.bin %s.outbin %s.lock',codepath,executionpath,md.name,md.name,md.name);
     44elseif strcmpi(solutiontype,'thermal') ,
     45        fprintf(fid,'mpirun.lsf %s/thermal.exe %s %s.bin %s.outbin %s.lock',codepath,executionpath,md.name,md.name,md.name);
    4646else
    47         error('BuildQueueingScriptcosmso error message: unsupported solution type!');
     47        error('BuildQueueingScriptcosmos error message: unsupported solution type!');
    4848end
    4949fclose(fid);

  • issm/trunk/src/m/classes/public/BuildQueueingScriptgemini.m

    r1 r465  
    2222fprintf(fid,'rm -rf %s.outlog %s.errlog %s.lock\n',md.name,md.name,md.name);
    2323
    24 if strcmpi(solutiontype,'diagnostic_horiz') |  strcmpi(solutiontype,'diagnostic'),
     24if strcmpi(solutiontype,'diagnostic') ,
    2525        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);
    2626elseif strcmpi(solutiontype,'control') ,
    2727        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);
    28 elseif strcmpi(solutiontype,'thermalsteady') ,
    29         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);
     28elseif strcmpi(solutiontype,'thermal') ,
     29        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);
    3030else
    3131        error('BuildQueueingScriptgemini error message: unsupported solution type!');

  • issm/trunk/src/m/classes/public/ismodelselfconsistent.m

    r440 r465  
    1 function bool=ismodelselfconsistent(md,solutiontype,package),
     1function bool=ismodelselfconsistent(md,package),
    22%ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem.
    33%
    44%   Usage:
    5 %      bool=ismodelselfconsistent(md,solutiontype,package),
     5%      bool=ismodelselfconsistent(md,package),
    66
    77%tolerance we use in litmus tests for the consistency of the model
    88tolerance=10^-12;
    9 if (nargin~=3  )
     9if (nargin~=2  )
    1010        help ismodelselfconsistent
    1111        error('ismodelselfconsistent error message: wrong usage');
    … …  
    115115
    116116%SIZE NUMBEROFGRIDS
    117 fields={'x','y','z','B','drag','gridondirichlet_diag','melting','accumulation','surface','thickness','bed','gridonbed','gridonsurface','deadgrids'};
     117fields={'x','y','z','B','drag','gridondirichlet_diag','melting','accumulation','surface','thickness','bed','gridonbed','gridonsurface'};
    118118for i=1:length(fields),
    119119        if find(length(eval(['md.' char(fields(i))]))~=md.numberofgrids),
    … …  
    159159
    160160%DIAGNOSTIC
    161 if strcmpi(solutiontype,'diagnostic')
     161if strcmpi(md.analysis_type,'diagnostic')
    162162
    163163        %HUTTER ON ICESHELF WARNING
    … …  
    189189
    190190%PROGNOSTIC
    191 if strcmp(solutiontype,'prognostic'),
     191if strcmp(md.analysis_type,'prognostic'),
    192192        %old testing
    193193        %       if isempty(find(md.gridondirichlet_diag)),
    … …  
    208208
    209209%THERMAL TRANSIENT
    210 if strcmp(solutiontype,'thermaltransient')
    211 
    212         %INITIAL TEMPERATURE, MELTING AND ACCUMULATION
    213         if isempty(md.temperature),
    214                 disp(['An  initial temperature is needed for a transient thermal computation'])
    215                 bool=0;return;
    216         end
    217         if isstruct(md.temperature) |  isstruct(md.melting) | isstruct(md.accumulation),
    218                 disp(['The initial temperature, melting or accumulation should be a list and not a structure'])
    219                 bool=0;return;
     210if strcmp(md.analysis_type,'thermal')
     211        if strcmp(md.sub_analysis_type,'transient')
     212
     213                %INITIAL TEMPERATURE, MELTING AND ACCUMULATION
     214                if isempty(md.temperature),
     215                        disp(['An  initial temperature is needed for a transient thermal computation'])
     216                        bool=0;return;
     217                end
     218                if isstruct(md.temperature) |  isstruct(md.melting) | isstruct(md.accumulation),
     219                        disp(['The initial temperature, melting or accumulation should be a list and not a structure'])
     220                        bool=0;return;
     221                end
    220222        end
    221223end
    222224
    223225%THERMAL STEADY AND THERMAL TRANSIENT
    224 if strcmpi(solutiontype,'thermalsteady') |  strcmp(solutiontype,'thermaltransient'),
     226if strcmpi(md.analysis_type,'thermal'),
    225227
    226228        %EXTRUSION
    227229        if strcmp(md.type,'2d'),
    228                 disp(['For a ' solutiontype ' computation, the model must be 3d, extrude it first!'])
     230                disp(['For a ' md.analysis_type ' computation, the model must be 3d, extrude it first!'])
    229231                bool=0;return;
    230232        end
    … …  
    242244
    243245%THERMAL TRANSIENT AND TRANSIENT
    244 if strcmp(solutiontype,'thermaltransient') | strcmp(solutiontype,'transient'),
     246if strcmp(md.analysis_type,'thermal'),
    245247
    246248        %DT and NDT
    … …  
    261263
    262264%PARAMETERS
    263 if strcmp(solutiontype,'parameters')
     265if strcmp(md.analysis_type,'parameters')
    264266
    265267        %OUTPUT
    … …  
    289291
    290292%CONTROL
    291 if strcmpi(solutiontype,'control'),
     293if strcmpi(md.analysis_type,'control'),
    292294
    293295        %CONTROL TYPE
    … …  
    320322
    321323%QMU
    322 if strcmpi(solutiontype,'qmu'),
     324if strcmpi(md.analysis_type,'qmu'),
    323325        if ~strcmpi(md.cluster,'none'),
    324326                if md.waitonlock==0,
    … …  
    330332
    331333%MESH
    332 if strcmpi(solutiontype,'mesh'),
     334if strcmpi(md.analysis_type,'mesh'),
    333335        %this solution is a little special. It should come right after the md=model;  operation. So a lot less checks!
    334336
    … …  
    338340
    339341%MESH2GRID
    340 if strcmpi(solutiontype,'mesh2grid'),
     342if strcmpi(md.analysis_type,'mesh2grid'),
    341343        if ~strcmpi(md.cluster,'none'),
    342344                disp(['model is not correctly configured: mesh2grid not supported in parallel yet!']);

  • issm/trunk/src/m/classes/public/isresultconsistent.m

    r27 r465  
    1 function bool=isresultconsistent(md,solutiontype),
     1function bool=isresultconsistent(md)
    22%ISRESULTCONSISTENT: check that results are consistent
    33%
    … …  
    66%
    77%   Usage:
    8 %      bool=IsModelSelfConsistent(model,solutiontype)
     8%      bool=IsModelSelfConsistent(model)
    99
    1010%tolerance we use in litmus tests for the consistency of the model
    1111tolerance=10^-2;
    1212
    13 if (nargin~=2  )
     13if (nargin~=1  )
    1414        IsResultConsistentUsage();
    1515        error(' ');
    1616end
    1717%Check velocities i(no penalties for 2d meshes)
    18 if (strcmp(solutiontype,'diagnostic') || strcmp(solutiontype,'diagnostic_horiz') ),
     18if strcmp(md.analysis_type,'diagnostic'),
    1919
    2020        if strcmpi(md.type,'3d'),
    … …  
    3737end
    3838
    39 if strcmp(solutiontype,'thermalsteady')
     39if strcmp(md.analysis_type,'thermal')
    4040
    4141        %check melting
    4242        if any(abs(md.melting(md.numberofgrids2d+1:end))>tolerance)
    43                 disp(['''thermalsteady'' result not consistent: increase penalty_melting (negative melting)']);
     43                disp(['''thermal'' result not consistent: increase penalty_melting (negative melting)']);
    4444                bool=0; return;
    4545        end
    … …  
    4747end
    4848
    49 if strcmp(solutiontype,'thermaltransient')
     49if strcmp(md.analysis_type,'thermaltransient')
    5050
    5151        for iter=1:length(md.thermaltransient_results)
    … …  
    6161
    6262
    63 if (strcmp(solutiontype,'transient') & strcmpi(md.type,'3d')),
     63if (strcmp(md.analysis_type,'transient') & strcmpi(md.type,'3d')),
    6464        for iter=1:length(md.transient_results)
    6565
    … …  
    9999function IsResultConsistentUsage()
    100100disp(' ');
    101 disp('   IsResultConsistent usage: flag=IsResultConsistent(model,solutiontype)');
     101disp('   IsResultConsistent usage: flag=IsResultConsistent(model)');
    102102disp('         where model is an instance of the @model class, and flag a boolean');

  • issm/trunk/src/m/classes/public/loadresultsfromcluster.m

    r227 r465  
    1 function md=loadresultsfromcluster(md,solutiontype)
     1function md=loadresultsfromcluster(md)
    22%LOADRESULTSFROMCLUSTER - load results of solution sequence from cluster
    33%
    44%   Usage:
    5 %      md=loadresultsfromcluster(md,solutiontype);
     5%      md=loadresultsfromcluster(md);
    66
    77%Get cluster.rc location
    … …  
    3131
    3232%If we are here, no errors in the solution sequence, call loadresultsfromdisk.
    33 md=loadresultsfromdisk(md,[md.name '.outbin'],solutiontype);
     33md=loadresultsfromdisk(md,[md.name '.outbin'],md.analysis_type);
    3434
    3535%erase the log and output files

  • issm/trunk/src/m/classes/public/loadresultsfromdisk.m

    r460 r465  
    1 function md=loadresultsfromdisk(md,filename,solutiontype)
     1function md=loadresultsfromdisk(md,filename)
    22%LOADRESULTSFROMDISK - load results of solution sequence from disk file "filename"           
    33%
    … …  
    103103disp(sprintf('%s\n','checking result consistency'));
    104104
    105 if ~isresultconsistent(md,solutiontype),
     105if ~isresultconsistent(md),
    106106        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...)
    107107end

  • issm/trunk/src/m/classes/public/marshall.m

    r381 r465  
    1 function marshall(md,solutiontype,package)
     1function marshall(md)
    22%MARSHALL - output Cielo compatible binary file from @model md, for certain solution type.
    33%
    44%   The routine creates a Cielo compatible binary file from @model md
    5 %   solutiontype can be 'diagnostic','prognostic' or 'control'.
    65%   This binary file will be used for parallel runs in cielo.
    76%
    87%   Usage:
    9 %      marshall(md,solutiontype,package)
     8%      marshall(md)
    109
    1110%some checks on list of arguments
    … …  
    2322end
    2423
    25 if strcmp(solutiontype,'diagnostic'),
    26         WriteData(fid,'diagnostic_horiz','String','analysis_type');
    27 else
    28         WriteData(fid,solutiontype,'String','analysis_type');
    29 end
     24WriteData(fid,md.analysis_type','String','analysis_type');
     25WriteData(fid,md.sub_analysis_type','String','sub_analysis_type');
    3026WriteData(fid,md.type,'String','type');
    3127WriteData(fid,md.numberofgrids,'Integer','numberofgrids');

  • issm/trunk/src/m/classes/public/solve.m

    r309 r465  
    1 function md=solve(md,solutiontype,varargin)
     1function md=solve(md,varargin)
    22%SOLVE - apply solution sequence for this model
    33%
    4 %   solutiontype is 'diagnostic','prognostic','transient','thermalsteady','thermaltransient','parameters' or 'control'
    5 %   and varargin is an optional package name ('cielo', 'ice' or 'macayeal')
    6 %
    74%   Usage:
    8 %       md=solve(md,solutiontype,varargin)
     5%       md=solve(md,varargin)
     6%       where varargin is a lit of paired arguments.
     7%       arguments can be: 'analysis_type': 'diagnostic','thermal','prognostic','transient'
     8%       arguments can be: 'sub_analysis_type': 'transient',''
     9%       arguments can be: 'package': 'macayeal','ice','cielo'
    910%
    1011%   Examples:
    11 %      md=solve(md,'diagnostic','macayeal');
    12 %      md=solve(md,'control','cielo');
     12%      md=solve(md,'analysis_type','diagnostic','package','cielo');
     13%      md=solve(md,'analysis_type','control','package','ice');
     14%      md=solve(md,'analysis_type','thermal','sub_analysis_type','transient','package','ice');
     15%      md=solve(md,'analysis_type','thermal','sub_analysis_type','','package','cielo');
    1316
    1417%some checks on list of arguments
    15 
    1618global ISSM_DIR
    1719
    18 solutions={'mesh2grid','mesh','thermaltransient','thermalsteady','qmu','diagnostic','diagnostic_horiz','prognostic','transient','parameters','control'};
    19 found=0;
    20 for i=1:length(solutions),
    21         if strcmpi(solutiontype,solutions{i}),
    22                 found=1;
    23                 break;
    24         end
    25 end
     20%recover options
     21options=recover_solve_options(md,varargin{:});
    2622
    27 if found==0,
    28         error(['solve error message: solution type ' solutiontype ' not supported yet!']);
    29 end
     23%add default options
     24options=process_solve_options(options);
    3025
    31 %we are good with this solutiontype, put in the analysis_type field of md:
    32 md.analysis_type=solutiontype;
    33 
    34 %Recover type of package being used:
    35 if nargin==2,
    36         package='Ice';
    37 else
    38         package=varargin{1};
    39 end
    40 
    41 if ~ischar(package),
    42         error('Package specified in varargin can only be ''ice'', or ''cielo''');
    43 end
    44 
    45 if ~(strcmpi(package,'ice') || strcmpi(package,'cielo') || strcmpi(package,'macayeal'))
    46         error('Package specified in varargin can only be ''ice'', ''macayeal'', or ''cielo''');
    47 end
     26%recover some fields
     27md.analysis_type=options.analysis_type;
     28md.sub_analysis_type=options.sub_analysis_type;
     29package=options.package;
    4830
    4931%Use package to set solution namespace
    5032usenamespace(package);
    5133
    52 if ~ismodelselfconsistent(md,solutiontype,package),
    53                 error(' '); %previous error messages should explain what is going on.
     34if ~ismodelselfconsistent(md,package),
     35        error(' '); %previous error messages should explain what is going on.
    5436end
    5537
    5638displaystring(md.debug,'\n%s\n','launching solution sequence');
    5739
     40
    5841%If running in parallel, we have a different way of launching the solution
    59 %sequences.
    60 if ~strcmpi(solutiontype,'qmu'),
     42%sequences. qmu is the only solution sequence that cannot run in parallel
     43if ~strcmpi(md.analysis_type,'qmu'),
    6144        if ~strcmpi(md.cluster,'none'),
    62                 md=solveparallel(md,solutiontype,package);
     45                md=solveparallel(md)
    6346                return;
    6447        end
    … …  
    6649
    6750%Lauch correct solution sequence
    68 if strcmpi(solutiontype,'diagnostic'),
     51if strcmpi(md.analysis_type,'diagnostic'),
    6952        md=diagnostic(md);
    7053
    71 elseif strcmpi(solutiontype,'mesh'),
     54elseif strcmpi(md.analysis_type,'mesh'),
    7255        md=mesh(md);
    7356
    74 elseif strcmpi(solutiontype,'transient'),
     57elseif strcmpi(md.analysis_type,'transient'),
    7558        md=transient(md);
    7659
    77 elseif strcmpi(solutiontype,'qmu'),
     60elseif strcmpi(md.analysis_type,'qmu'),
    7861        md=qmu(md,package);
    7962
    80 elseif strcmpi(solutiontype,'parameters'),
    81         md=parameters(md);;
    82 
    83 elseif strcmpi(solutiontype,'mesh2grid'),
     63elseif strcmpi(md.analysis_type,'mesh2grid'),
    8464        md=mesh2grid(md);;
    8565
    86 elseif strcmpi(solutiontype,'prognostic'),
     66elseif strcmpi(md.analysis_type,'prognostic'),
    8767        md=prognostic(md);;
    8868
    89 elseif strcmpi(solutiontype,'control'),
     69elseif strcmpi(md.analysis_type,'control'),
    9070        md=control(md);
    9171
    92 elseif strcmpi(solutiontype,'thermalsteady') | strcmpi(solutiontype,'thermaltransient'),
    93         md=thermal(md,solutiontype);
     72elseif strcmpi(md.analysis_type,'thermal'),
     73        md=thermal(md);
    9474else
    9575        error('solution type not supported for this model configuration.');
    … …  
    9878%Check result is consistent
    9979displaystring(md.debug,'%s\n','checking result consistency');
    100 if ~isresultconsistent(md,solutiontype),
     80if ~isresultconsistent(md),
    10181        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...)
    10282end
    … …  
    11090function solveusage();
    11191disp(' ');
    112 disp('   Solve usage: md=solve(md,solutiontype,package)');
    113 disp('         solutiontype can be ''diagnostic'',''transient'', ''thermalsteady'',''thermaltransient'',''parameters'',''mesh2grid''  or ''control''');
     92disp('   Solve usage: md=solve(md,md.analysis_type,package)');
     93disp('         md.analysis_type can be ''diagnostic'',''transient'', ''thermal'',''thermaltransient'',''parameters'',''mesh2grid''  or ''control''');
    11494disp('         package is either ''cielo'' or ''ice''');

  • issm/trunk/src/m/classes/public/solveparallel.m

    r227 r465  
    1 function md=solveparallel(md,solutiontype,varargin)
     1function md=solveparallel(md)
    22%SOLVEPARALLEL - solution sequence using a cluster in parallel mode
    33%
    44%   Usage:
    5 %      md=solveparallel(md,solutiontype,varargin)
    6 
    7 %Recover type of package being used:
    8 if nargin==2,
    9         package='Ice';
    10 else
    11         package=varargin{1};
    12 end
    13 
    14 if ~ischar(package),
    15         error('Package specified in varargin can only be ''ice'', or ''cielo''');
    16 end
    17 
    18 if ~(strcmpi(package,'ice') || strcmpi(package,'cielo') || strcmpi(package,'macayeal'))
    19         error('Package specified in varargin can only be ''ice'', ''macayeal'', or ''cielo''');
    20 end
     5%      md=solveparallel(md);
    216
    227%Get cluster.rc location
    … …  
    2712
    2813%Marshall model data into a binary file.
    29 marshall(md,solutiontype,package);
     14marshall(md);
    3015
    3116%Now, we need to build the queuing script, used by the cluster to launch the job.
    32 BuildQueueingScript(md,solutiontype,executionpath,codepath);
     17BuildQueueingScript(md,executionpath,codepath);
    3318
    3419%Now, launch the queueing script
    … …  
    4025        waitonlock([executionpath '/' md.name '.lock']);
    4126        %load results
    42         md=loadresultsfromcluster(md,solutiontype);
     27        md=loadresultsfromcluster(md);
    4328else
    4429        return;

  • issm/trunk/src/m/solutions/cielo/GradJCompute.m

    r422 r465  
    1 function [u_g grad_g]=GradJCompute(m,inputs, u_g_obs,analysis_type);
     1function [u_g grad_g]=GradJCompute(m,inputs, u_g_obs,analysis_type,sub_analysis_type);
    22
    33%Recover solution for this stiffness and right hand side:
    … …  
    55        disp('         computing velocities...')
    66end
    7 [u_g K_ff0 K_fs0 ]=diagnostic_core_nonlinear(m,inputs,analysis_type);
     7[u_g K_ff0 K_fs0 ]=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type);
    88
    99%Buid Du, difference between observed velocity and model velocity.
    … …  
    1111        disp('         computing Du...')
    1212end
    13 [Du_g]=Du(m.elements,m.nodes,m.loads,m.materials,m.parameters, u_g,u_g_obs,inputs);
     13[Du_g]=Du(m.elements,m.nodes,m.loads,m.materials,m.parameters, u_g,u_g_obs,inputs,analysis_type,sub_analysis_type);
    1414
    1515%Reduce adjoint load from g-set to f-set
    … …  
    2626
    2727%Compute gradJ
    28 grad_g=Gradj(m.elements,m.nodes,m.loads,m.materials,m.parameters, u_g, lambda_g,inputs);
     28grad_g=Gradj(m.elements,m.nodes,m.loads,m.materials,m.parameters, u_g, lambda_g,inputs,analysis_type,sub_analysis_type);

  • issm/trunk/src/m/solutions/cielo/control.m

    r422 r465  
    66        %Build all models requested for control simulation
    77        md.analysis_type='control';
     8        md.sub_analysis_type='';
    89        m=CreateFemModel(md);
    910
    … …  
    4041
    4142                disp('      computing gradJ...');
    42                 [u_g c(n).grad_g]=GradJCompute(m,inputs,u_g_obs,md.analysis_type);
     43                [u_g c(n).grad_g]=GradJCompute(m,inputs,u_g_obs,md.analysis_type,md.sub_analysis_type);
    4344                inputs=add(inputs,'velocity',u_g,'doublevec',2,m.parameters.numberofnodes);
    4445                disp('      done.');
    … …  
    5859               
    5960                disp('      optimizing along gradient direction...');
    60                 [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);
     61                [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);
    6162                disp('      done.');
    6263
    … …  
    7677                %some temporary saving
    7778                if(mod(n,5)==0),
    78                         solution=controlfinalsol(c,m,p_g,inputs,md.analysis_type);
     79                        solution=controlfinalsol(c,m,p_g,inputs,md.analysis_type,md.sub_analysis_type);
    7980                        save temporary_control_results solution
    8081                end
    … …  
    8586        %Create final solution
    8687        disp('      preparing final velocity solution...');
    87         solution=controlfinalsol(c,m,p_g,inputs,md.analysis_type);
     88        solution=controlfinalsol(c,m,p_g,inputs,md.analysis_type,md.sub_analysis_type);
    8889        disp('      done.');
    8990       

  • issm/trunk/src/m/solutions/cielo/controlfinalsol.m

    r377 r465  
    1 function solution=controlfinalsol(c,m,p_g,inputs,analysis_type);
     1function solution=controlfinalsol(c,m,p_g,inputs,analysis_type,sub_analysis_type);
    22
    33%From parameters, build inputs for icediagnostic_core, using the final parameters
    44inputs=add(inputs,m.parameters.control_type,p_g,'doublevec',2,m.parameters.numberofnodes);
    5 u_g=diagnostic_core_nonlinear(m,inputs,analysis_type);
     5u_g=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type);
    66
    77%Build partitioning vectors to recover solution

  • issm/trunk/src/m/solutions/cielo/diagnostic.m

    r358 r465  
    1010        %Build all models requested for diagnostic simulation
    1111        displaystring(md.debug,'%s',['reading diagnostic horiz model data']);
    12         md.analysis_type='diagnostic_horiz'; m_dh=CreateFemModel(md);
     12        md.analysis_type='diagnostic'; md.sub_analysis_type='horiz'; m_dh=CreateFemModel(md);
    1313       
    1414        displaystring(md.debug,'\n%s',['reading diagnostic vert model data']);
    15         md.analysis_type='diagnostic_vert'; m_dv=CreateFemModel(md);
     15        md.analysis_type='diagnostic'; md.sub_analysis_type='vert'; m_dv=CreateFemModel(md);
    1616       
    1717        displaystring(md.debug,'\n%s',['reading diagnostic stokes model data']);
    18         md.analysis_type='diagnostic_stokes'; m_ds=CreateFemModel(md);
     18        md.analysis_type='diagnostic'; md.sub_analysis_type='stokes'; m_ds=CreateFemModel(md);
    1919
    2020        displaystring(md.debug,'\n%s',['reading diagnostic hutter model data']);
    21         md.analysis_type='diagnostic_hutter'; m_dhu=CreateFemModel(md);
     21        md.analysis_type='diagnostic'; md.sub_analysis_type='hutter'; m_dhu=CreateFemModel(md);
    2222       
    2323        displaystring(md.debug,'\n%s',['reading surface and bed slope computation model data']);
    24         md.analysis_type='slope_compute'; m_sl=CreateFemModel(md);
     24        md.analysis_type='slope_compute'; md.sub_analysis_type=''; m_sl=CreateFemModel(md);
    2525       
    2626        % figure out number of dof: just for information purposes.

  • issm/trunk/src/m/solutions/cielo/diagnostic_core.m

    r394 r465  
    1616
    1717        displaystring(debug,'\n%s',['computing surface slope (x and y derivatives)...']);
    18         slopex=diagnostic_core_linear(m_sl,inputs,'surface_slope_compute_x');
    19         slopey=diagnostic_core_linear(m_sl,inputs,'surface_slope_compute_y');
     18        slopex=diagnostic_core_linear(m_sl,inputs,'slope_compute','surfacex');
     19        slopey=diagnostic_core_linear(m_sl,inputs,'slope_compute','surfacey');
    2020
    2121        if dim==3,
    … …  
    3030
    3131        displaystring(debug,'\n%s',['computing hutter velocities...']);
    32         u_g=diagnostic_core_linear(m_dhu,inputs,'diagnostic_hutter');
     32        u_g=diagnostic_core_linear(m_dhu,inputs,'diagnostic','hutter');
    3333
    3434        displaystring(debug,'\n%s',['computing pressure according to MacAyeal...']);
    … …  
    4646
    4747        displaystring(debug,'\n%s',['computing horizontal velocities...']);
    48         u_g=diagnostic_core_nonlinear(m_dh,inputs,'diagnostic_horiz');
     48        u_g=diagnostic_core_nonlinear(m_dh,inputs,'diagnostic','horiz');
    4949
    5050        displaystring(debug,'\n%s',['computing pressure according to MacAyeal...']);
    … …  
    5959        displaystring(debug,'\n%s',['computing vertical velocities...']);
    6060        inputs=add(inputs,'velocity',u_g_horiz,'doublevec',m_dh.parameters.numberofdofspernode,m_dh.parameters.numberofnodes);
    61         u_g_vert=diagnostic_core_linear(m_dv,inputs,'diagnostic_vert');
     61        u_g_vert=diagnostic_core_linear(m_dv,inputs,'diagnostic','vert');
    6262
    6363        displaystring(debug,'\n%s',['combining horizontal and vertical velocities...']);
    … …  
    7575
    7676                displaystring(debug,'\n%s',['computing bed slope (x and y derivatives)...']);
    77                 slopex=diagnostic_core_linear(m_sl,inputs,'bed_slope_compute_x');
    78                 slopey=diagnostic_core_linear(m_sl,inputs,'bed_slope_compute_y');
     77                slopex=diagnostic_core_linear(m_sl,inputs,'slope_compute','bedx');
     78                slopey=diagnostic_core_linear(m_sl,inputs,'slope_compute','bedy');
    7979                slopex=SlopeExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,slopex);
    8080                slopey=SlopeExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,slopey);
    … …  
    9595
    9696                displaystring(debug,'\n%s',['computing stokes velocities and pressure ...']);
    97                 u_g=diagnostic_core_nonlinear(m_ds,inputs,'diagnostic_stokes');
     97                u_g=diagnostic_core_nonlinear(m_ds,inputs,'diagnostic','stokes');
    9898       
    9999                %"decondition" pressure

  • issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m

    r394 r465  
    1 function [u_g varargout]=diagnostic_core_nonlinear(m,inputs,analysis_type)
     1function [u_g varargout]=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type)
    22%INPUT function [ru_g varargout]=cielodiagnostic_core_nonlinear(m,inputs)
    33       
    … …  
    3030               
    3131                %system matrices
    32                 [K_gg_nopenalty , p_g_nopenalty]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type);
     32                [K_gg_nopenalty , p_g_nopenalty]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
    3333               
    3434                %penalties
    35                 [K_gg , p_g]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type);
     35                [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);
    3636
    3737
    … …  
    6060       
    6161                %penalty constraints
    62                 [loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints( m.elements,m.nodes, loads, m.materials,m.parameters,inputs);
     62                [loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints( m.elements,m.nodes, loads, m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
    6363
    6464        %   Figure out if convergence is reached.
    … …  
    107107                inputs=add(inputs,'velocity',soln(count).u_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes);
    108108                m.parameters.kflag=1; m.parameters.pflag=0;
    109                 [K_gg, p_g]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type);
     109                [K_gg, p_g]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
    110110                [K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets);
    111111                varargout(1)={K_ff};

  • issm/trunk/src/m/solutions/cielo/objectivefunctionC.m

    r377 r465  
    1 function J =objectivefunctionC(search_scalar,m,inputs,p_g,u_g_obs,grad_g,n,analysis_type);
     1function J =objectivefunctionC(search_scalar,m,inputs,p_g,u_g_obs,grad_g,n,analysis_type,sub_analysis_type);
    22       
    33%recover some parameters
    … …  
    1313
    1414%Run diagnostic with updated parameters.
    15 u_g=diagnostic_core_nonlinear(m,inputs,analysis_type);
     15u_g=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type);
    1616
    1717%Compute misfit for this velocity field.
    18 J=Misfit(m.elements,m.nodes,m.loads,m.materials,m.parameters, u_g, u_g_obs,inputs);
     18J=Misfit(m.elements,m.nodes,m.loads,m.materials,m.parameters, u_g, u_g_obs,inputs,analysis_type,sub_analysis_type);

  • issm/trunk/src/m/solutions/cielo/thermal.m

    r163 r465  
    77        %timing
    88        t1=clock;
    9        
    10         %md.analysis_type is set to thermalsteady or thermaltransient. save it.
    11         solutiontype=md.analysis_type;
    129
    13         %Build all models requested for thermal and melting  simulation
    14         md.analysis_type=solutiontype; m_t=CreateFemModel(md);
    15         md.analysis_type='melting' m_m=CreateFEMModel(md);
     10        %Build all models requested for diagnostic simulation
     11        displaystring(md.debug,'%s',['reading thermal model data']);
     12        md.analysis_type='thermal'; m_t=CreateFemModel(md);
     13
     14        displaystring(md.debug,'%s',['reading melting model data']);
     15        md.analysis_type='melting'; m_m=CreateFemModel(md);
    1616
    1717        % figure out number of dof: just for information purposes.
    18         md.dof=m_t.nodesets.fsize; %biggest dof number
     18        md.dof=m_t.nodesets.fsize;
    1919
    20 
    21         if strcmpi(solutiontype,'thermalsteady'),
    22 
    23                 %initialize velocity and pressure
    24                 u_g=m_t.parameters.u_g;
    25                 p_g=m_t.parameters.p_g;
    26 
    27                 inputs=struct('pressure',p_g,'velocity',u_g);
    28 
    29                 disp(sprintf('\n%s\n','computing temperature...'));
    30                 [t_g,m_t]=thermal_core(m_t,inputs);
    31 
    32                 disp(sprintf('\n%s\n','computing melting...'));
    33                 inputs=struct('pressure',p_g,'temperature',t_g);
     20        %initialize inputs
     21        displaystring(debug,'\n%s',['setup inputs...']);
     22        inputs=inputlist;
     23        inputs=add(inputs,'velocity',m_t.parameters.u_g,'doublevec',3,m_t.parameters.numberofnodes);
     24        inputs=add(inputs,'pressure',m_t.parameters.p_g,'doublevec',1,m_t.parameters.numberofnodes);
     25        inputs=add(inputs,'dt',m_t.parameters.dt,'double');
     26       
     27        if strcmpi(md.sub_analysis_type,'steady'),
     28       
     29                displaystring(debug,'\n%s',['computing temperatures...']);
     30                [t_g loads_t melting_offset]=thermal_core(m_t,inputs,'thermal','steady');
    3431               
    35                 melting_g=melting_core(m_m,inputs);
    36 
    37                 %load results onto model
     32                displaystring(debug,'\n%s',['computing melting...']);
     33                inputs=add(inputs,'melting_offset',melting_offset,'double');
     34                melting_g=melting_core(m_m,inputs,'melting','steady');
     35               
     36                displaystring(debug,'\n%s',['load results...']);
    3837                md.temperature=t_g;
    3938                md.melting=melting_g*md.yts; %from m/s to m/a
    … …  
    4140        else
    4241
    43                 %initialize velocity, pressure  and temperature
    44                 u_g=m.parameters.u_g;
    45                 p_g=m.parameters.p_g;
    46                 t_g=m.parameters.t_g;
    47                 melting_g=m.parameters.melting_g;
     42                %initialize velocity, pressure, temperature and melting
     43                u_g=m_t.parameters.u_g;
     44                p_g=m_t.parameters.p_g;
     45                t_g=m_t.parameters.t_g;
     46                melting_g=m_t.parameters.melting_g;
    4847
    4948                nsteps=m_t.parameters.ndt/m_t.parameters.dt;
    … …  
    5554                for n=1:nsteps,
    5655
    57                         disp(sprintf('\n%s%i/%i\n','time step: ',n,md.ndt/md.dt));
    58                         disp('   computing temperature...');
    59                
    60                         inputs=struct('pressure',p_g,'temperature',soln(n).t_g,'velocity',u_g,'dt',m_t.parameters.dt);
    61                         [soln(n+1).t_g,m_t]=thermal_core(m_t,inputs);
     56                        displaystring(debug,'\n%s%i/%i\n','time step: ',n,nsteps);
     57                       
     58                        displaystring(debug,'\n%s',['    computing temperatures...']);
     59                        inputs=add(inputs,'temperature',soln(n).t_g,'doublevec',1,m_t.parameters.numberofnodes);
     60                        [soln(n+1).t_g loads_t melting_offset]=thermal_core(m_t,inputs);
    6261                       
    6362                        disp('   computing melting...');
    … …  
    8079                md.melting=solution_melting;
    8180        end
    82        
     81
     82        %stop timing
     83        t2=clock;
     84        displaystring(md.debug,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']);   

  • issm/trunk/src/mex/ControlOptimization/ControlOptimization.cpp

    r377 r465  
    5151        optargs.n=STEP;
    5252        optargs.analysis_type=ANALYSIS;
     53        optargs.sub_analysis_type=SUBANALYSIS;
    5354
    5455        optpars.xmin=xmin;
    … …  
    7374{
    7475        _printf_("\n");
    75         _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__);
     76        _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__);
    7677        _printf_("\n");
    7778}

  • issm/trunk/src/mex/ControlOptimization/ControlOptimization.h

    r377 r465  
    2828#define STEP (mxArray*)prhs[9]
    2929#define ANALYSIS (mxArray*)prhs[10]
     30#define SUBANALYSIS (mxArray*)prhs[11]
    3031
    3132/* serial output macros: */
    … …  
    3738#define NLHS  2
    3839#undef NRHS
    39 #define NRHS  11
     40#define NRHS  12
    4041
    4142

  • issm/trunk/src/mex/Du/Du.cpp

    r246 r465  
    1515        DataSet* loads=NULL;
    1616        DataSet* materials=NULL;
    17         int      analysis_type;
    1817        double*  u_g=NULL;
    1918        double*  u_g_obs=NULL;
    2019        ParameterInputs* inputs=NULL;
     20        char*             analysis_type_string=NULL;
     21        int               analysis_type;
     22        char*             sub_analysis_type_string=NULL;
     23        int               sub_analysis_type;
    2124
    2225        /* output datasets: */
    … …  
    3538        FetchData((void**)&loads,NULL,NULL,LOADS,"DataSet",NULL);
    3639        FetchData((void**)&materials,NULL,NULL,MATERIALS,"DataSet",NULL);
    37         FetchData((void**)&analysis_type,NULL,NULL,mxGetField(PARAMETERS,0,"analysis_type"),"Integer",NULL);
    3840        FetchData((void**)&u_g,NULL,NULL,UG,"Vector","Vec");
    3941        FetchData((void**)&u_g_obs,NULL,NULL,UGOBS,"Vector","Vec");
     42        FetchData((void**)&analysis_type_string,NULL,NULL,ANALYSIS,"String",NULL);
     43        analysis_type=AnalysisTypeAsEnum(analysis_type_string);
     44        FetchData((void**)&sub_analysis_type_string,NULL,NULL,SUBANALYSIS,"String",NULL);
     45        sub_analysis_type=AnalysisTypeAsEnum(sub_analysis_type_string);
    4046
    4147        /*Fetch inputs: */
    … …  
    4450
    4551        /*!Call core code: */
    46         Dux(&du_g, elements,nodes,loads,materials,u_g,u_g_obs,inputs,analysis_type);
     52        Dux(&du_g, elements,nodes,loads,materials,u_g,u_g_obs,inputs,analysis_type,sub_analysis_type);
    4753
    4854        /*write output : */
    … …  
    5864        VecFree(&du_g);
    5965        delete inputs;
     66        xfree((void**)&analysis_type_string);
     67        xfree((void**)&sub_analysis_type_string);
    6068
    6169        /*end module: */

  • issm/trunk/src/mex/Du/Du.h

    r1 r465  
    2525#define UGOBS (mxArray*)prhs[6]
    2626#define INPUTS (mxArray*)prhs[7]
     27#define ANALYSIS (mxArray*)prhs[8]
     28#define SUBANALYSIS (mxArray*)prhs[9]
    2729
    2830/* serial output macros: */
    … …  
    3335#define NLHS  1
    3436#undef NRHS
    35 #define NRHS  8
     37#define NRHS  10
    3638
    3739

  • issm/trunk/src/mex/Gradj/Gradj.cpp

    r246 r465  
    1515        DataSet* loads=NULL;
    1616        DataSet* materials=NULL;
    17         int      analysis_type;
    1817        char*    control_type=NULL;
    1918        double*  u_g=NULL;
    2019        double*  lambda_g=NULL;
    2120        ParameterInputs* inputs=NULL;
     21        char*             analysis_type_string=NULL;
     22        int               analysis_type;
     23        char*             sub_analysis_type_string=NULL;
     24        int               sub_analysis_type;
    2225
    2326        /* output datasets: */
    … …  
    3639        FetchData((void**)&loads,NULL,NULL,LOADS,"DataSet",NULL);
    3740        FetchData((void**)&materials,NULL,NULL,MATERIALS,"DataSet",NULL);
    38         FetchData((void**)&analysis_type,NULL,NULL,mxGetField(PARAMETERS,0,"analysis_type"),"Integer",NULL);
    3941        FetchData((void**)&control_type,NULL,NULL,mxGetField(PARAMETERS,0,"control_type"),"String",NULL);
    4042        FetchData((void**)&u_g,NULL,NULL,UG,"Vector","Vec");
    4143        FetchData((void**)&lambda_g,NULL,NULL,LAMBDAG,"Vector","Vec");
     44        FetchData((void**)&analysis_type_string,NULL,NULL,ANALYSIS,"String",NULL);
     45        analysis_type=AnalysisTypeAsEnum(analysis_type_string);
     46        FetchData((void**)&sub_analysis_type_string,NULL,NULL,SUBANALYSIS,"String",NULL);
     47        sub_analysis_type=AnalysisTypeAsEnum(sub_analysis_type_string);
    4248
    4349        /*Fetch inputs: */
    … …  
    4652
    4753        /*!Call core code: */
    48         Gradjx(&grad_g, elements,nodes,loads,materials,u_g,lambda_g,inputs,analysis_type,control_type);
     54        Gradjx(&grad_g, elements,nodes,loads,materials,u_g,lambda_g,inputs,analysis_type,sub_analysis_type,control_type);
    4955
    5056        /*write output : */
    … …  
    6167        delete inputs;
    6268        VecFree(&grad_g);
     69        xfree((void**)&analysis_type_string);
     70        xfree((void**)&sub_analysis_type_string);
    6371
    6472        /*end module: */

  • issm/trunk/src/mex/Gradj/Gradj.h

    r1 r465  
    2525#define LAMBDAG (mxArray*)prhs[6]
    2626#define INPUTS (mxArray*)prhs[7]
     27#define ANALYSIS (mxArray*)prhs[8]
     28#define SUBANALYSIS (mxArray*)prhs[9]
    2729
    2830/* serial output macros: */
    … …  
    3335#define NLHS  1
    3436#undef NRHS
    35 #define NRHS  8
     37#define NRHS  10
    3638
    3739

  • issm/trunk/src/mex/Misfit/Misfit.cpp

    r246 r465  
    1515        DataSet* loads=NULL;
    1616        DataSet* materials=NULL;
    17         int      analysis_type;
    1817        double*  u_g=NULL;
    1918        double*  u_g_obs=NULL;
    2019        ParameterInputs* inputs=NULL;
     20        char*             analysis_type_string=NULL;
     21        int               analysis_type;
     22        char*             sub_analysis_type_string=NULL;
     23        int               sub_analysis_type;
    2124
    2225        /* output datasets: */
    … …  
    3437        FetchData((void**)&loads,NULL,NULL,LOADS,"DataSet",NULL);
    3538        FetchData((void**)&materials,NULL,NULL,MATERIALS,"DataSet",NULL);
    36         FetchData((void**)&analysis_type,NULL,NULL,mxGetField(PARAMETERS,0,"analysis_type"),"Integer",NULL);
    3739        FetchData((void**)&u_g,NULL,NULL,UG,"Vector","Vec");
    3840        FetchData((void**)&u_g_obs,NULL,NULL,UGOBS,"Vector","Vec");
     41        FetchData((void**)&analysis_type_string,NULL,NULL,ANALYSIS,"String",NULL);
     42        analysis_type=AnalysisTypeAsEnum(analysis_type_string);
     43        FetchData((void**)&sub_analysis_type_string,NULL,NULL,SUBANALYSIS,"String",NULL);
     44        sub_analysis_type=AnalysisTypeAsEnum(sub_analysis_type_string);
    3945
    4046        /*Fetch inputs: */
    … …  
    4349
    4450        /*!Call core code: */
    45         Misfitx(&J, elements,nodes,loads,materials,u_g,u_g_obs,inputs,analysis_type);
     51        Misfitx(&J, elements,nodes,loads,materials,u_g,u_g_obs,inputs,analysis_type,sub_analysis_type);
    4652
    4753        /*write output : */
    … …  
    5561        xfree((void**)&u_g);
    5662        xfree((void**)&u_g_obs);
    57         delete inputs
     63        delete inputs;
     64        xfree((void**)&analysis_type_string);
     65        xfree((void**)&sub_analysis_type_string);
    5866
    5967        /*end module: */

  • issm/trunk/src/mex/Misfit/Misfit.h

    r1 r465  
    2525#define UGOBS (mxArray*)prhs[6]
    2626#define INPUTS (mxArray*)prhs[7]
     27#define ANALYSIS (mxArray*)prhs[8]
     28#define SUBANALYSIS (mxArray*)prhs[9]
    2729
    2830/* serial output macros: */
    … …  
    3335#define NLHS  1
    3436#undef NRHS
    35 #define NRHS  8
     37#define NRHS  10
    3638
    3739

  • issm/trunk/src/mex/ModelProcessor/ModelProcessor.cpp

    r394 r465  
    3232
    3333        /*Create elements, nodes and materials: */
    34         CreateElementsNodesAndMaterials(&elements,&nodes,&materials,model,MODEL);
    35 
    36         /*Create constraints: */
    37         CreateConstraints(&constraints,model,MODEL);
    38 
    39         /*Create loads: */
    40         CreateLoads(&loads,model,MODEL);
    41         /*Create parameters: */
    42         CreateParameters(&parameters,model,MODEL);
     34        CreateDataSets(&elements,&nodes,&materials,&constraints, &loads, &parameters, model,MODEL);
    4335
    4436        /*Write output data: */

  • issm/trunk/src/mex/PenaltyConstraints/PenaltyConstraints.cpp

    r246 r465  
    1616        DataSet* materials=NULL;
    1717        ParameterInputs* inputs=NULL;
     18        char*             analysis_type_string=NULL;
    1819        int               analysis_type;
     20        char*             sub_analysis_type_string=NULL;
     21        int               sub_analysis_type;
    1922
    2023        /*output: */
    … …  
    3538       
    3639        /*parameters: */
    37         FetchData((void**)&analysis_type,NULL,NULL,mxGetField(PARAMETERS,0,"analysis_type"),"Integer",NULL);
     40        FetchData((void**)&analysis_type_string,NULL,NULL,ANALYSIS,"String",NULL);
     41        analysis_type=AnalysisTypeAsEnum(analysis_type_string);
     42
     43        FetchData((void**)&sub_analysis_type_string,NULL,NULL,SUBANALYSIS,"String",NULL);
     44        sub_analysis_type=AnalysisTypeAsEnum(sub_analysis_type_string);
    3845
    3946        /*Fetch inputs: */
    … …  
    4249
    4350        /*!Generate internal degree of freedom numbers: */
    44         PenaltyConstraintsx(&converged, &num_unstable_constraints, elements,nodes,loads,materials,inputs,analysis_type);
     51        PenaltyConstraintsx(&converged, &num_unstable_constraints, elements,nodes,loads,materials,inputs,analysis_type,sub_analysis_type);
    4552
    4653        /*write output datasets: */
    … …  
    5663        delete materials;
    5764        delete inputs;
     65        xfree((void**)&analysis_type_string);
     66        xfree((void**)&sub_analysis_type_string);
    5867
    5968        /*end module: */

  • issm/trunk/src/mex/PenaltyConstraints/PenaltyConstraints.h

    r1 r465  
    2323#define PARAMETERS (mxArray*)prhs[4]
    2424#define INPUTS (mxArray*)prhs[5]
     25#define ANALYSIS (mxArray*)prhs[6]
     26#define SUBANALYSIS (mxArray*)prhs[7]
    2527
    2628/* serial output macros: */
    … …  
    3335#define NLHS  3
    3436#undef NRHS
    35 #define NRHS  6
     37#define NRHS  8
    3638
    3739#endif  /* _PENALTYCONSTRAINTS_H */

  • issm/trunk/src/mex/PenaltySystemMatrices/PenaltySystemMatrices.cpp

    r364 r465  
    2121        int               analysis_type;
    2222        char*             analysis_type_string=NULL;
     23        int               sub_analysis_type;
     24        char*             sub_analysis_type_string=NULL;
    2325       
    2426        /*Boot module: */
    … …  
    4244        analysis_type=AnalysisTypeAsEnum(analysis_type_string);
    4345
     46        FetchData((void**)&sub_analysis_type_string,NULL,NULL,SUBANALYSIS,"String",NULL);
     47        sub_analysis_type=AnalysisTypeAsEnum(sub_analysis_type_string);
     48
    4449
    4550        /*Fetch inputs: */
    … …  
    4853
    4954        /*!Generate stiffnesses from penalties: */
    50         PenaltySystemMatricesx(Kgg, pg,elements,nodes,loads,materials,kflag,pflag,inputs,analysis_type);
     55        PenaltySystemMatricesx(Kgg, pg,elements,nodes,loads,materials,kflag,pflag,inputs,analysis_type,sub_analysis_type);
    5156
    5257        /*write output datasets: */
    … …  
    6368        VecFree(&pg);
    6469        xfree((void**)&analysis_type_string);
     70        xfree((void**)&sub_analysis_type_string);
    6571
    6672        /*end module: */

  • issm/trunk/src/mex/PenaltySystemMatrices/PenaltySystemMatrices.h

    r364 r465  
    2626#define INPUTS (mxArray*)prhs[7]
    2727#define ANALYSIS (mxArray*)prhs[8]
     28#define SUBANALYSIS (mxArray*)prhs[9]
    2829
    2930/* serial output macros: */
    … …  
    3536#define NLHS  2
    3637#undef NRHS
    37 #define NRHS  9
     38#define NRHS  10
    3839
    3940

  • issm/trunk/src/mex/SystemMatrices/SystemMatrices.cpp

    r371 r465  
    2121        char*             analysis_type_string=NULL;
    2222        int               analysis_type;
     23        char*             sub_analysis_type_string=NULL;
     24        int               sub_analysis_type;
    2325       
    2426        /* output datasets: */
    … …  
    4648        analysis_type=AnalysisTypeAsEnum(analysis_type_string);
    4749
     50        FetchData((void**)&sub_analysis_type_string,NULL,NULL,SUBANALYSIS,"String",NULL);
     51        sub_analysis_type=AnalysisTypeAsEnum(sub_analysis_type_string);
     52
    4853        /*Fetch inputs: */
    4954        inputs=new ParameterInputs;
    … …  
    5156
    5257        /*!Generate internal degree of freedom numbers: */
    53         SystemMatricesx(&Kgg, &pg,elements,nodes,loads,materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type);
     58        SystemMatricesx(&Kgg, &pg,elements,nodes,loads,materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type,sub_analysis_type);
    5459
    5560        /*write output datasets: */
    … …  
    6065        /*Free ressources: */
    6166        xfree((void**)&analysis_type_string);
     67        xfree((void**)&sub_analysis_type_string);
    6268        delete elements;
    6369        delete nodes;

  • issm/trunk/src/mex/SystemMatrices/SystemMatrices.h

    r364 r465  
    2424#define INPUTS (mxArray*)prhs[5]
    2525#define ANALYSIS (mxArray*)prhs[6]
     26#define SUBANALYSIS (mxArray*)prhs[7]
    2627
    2728/* serial output macros: */
    … …  
    3334#define NLHS  2
    3435#undef NRHS
    35 #define NRHS  7
     36#define NRHS  8
    3637
    3738
Note: See TracChangeset for help on using the changeset viewer.