Changeset 1784


Ignore:
Timestamp:
08/20/09 09:30:06 (15 years ago)
Author:
Eric.Larour
Message:

Rift formulation

Location:
issm/trunk/src
Files:
20 edited

Legend:

Unmodified
Added
Removed

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

    r1628 r1784  
    11461146
    11471147#undef __FUNCT__
    1148 #define __FUNCT__ "RiftConstraints"
    1149 void  DataSet::RiftConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type){
    1150        
    1151         throw ErrorException(__FUNCT__," not implemented yet!");
    1152 
    1153 }
    1154 
    1155 #undef __FUNCT__
    11561148#define __FUNCT__ "MeltingIsPresent"
    11571149int   DataSet::MeltingIsPresent(){
    … …  
    11811173#undef __FUNCT__
    11821174#define __FUNCT__ "MeltingConstraints"
    1183 void  DataSet::MeltingConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type){
     1175void  DataSet::MeltingConstraints(int* pconverged, int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type){
    11841176
    11851177        /* generic object pointer: */
    … …  
    11891181        int unstable=0;
    11901182        int num_unstable_constraints=0;
     1183        int converged=0;
    11911184        int sum_num_unstable_constraints=0;
    11921185
    … …  
    12121205        #endif
    12131206
     1207        /*Have we converged? : */
     1208        if (num_unstable_constraints==0) converged=1;
     1209        else converged=0;
     1210
    12141211        /*Assign output pointers: */
     1212        *pconverged=converged;
    12151213        *pnum_unstable_constraints=num_unstable_constraints;
    12161214}

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

    r1628 r1784  
    7070                void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
    7171                int   RiftIsPresent();
    72                 void  RiftConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type);
    7372                int   MeltingIsPresent();
    74                 void  MeltingConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type);
     73                void  MeltingConstraints(int* pconverged, int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type);
    7574                DataSet* Copy(void);
    7675                void  Du(Vec du_g,void* inputs,int analysis_type,int sub_analysis_type);

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

    r1728 r1784  
    9292/*GEOGRAPHY:*/
    9393int GeographyEnum(void){                return          500; }
    94 int WaterEnum(void){                    return          501; }
    9594int IceSheetEnum(void){                 return          502; }
    9695int IceShelfEnum(void){                 return          502; }
    9796
    9897/*FILL:*/
    99 int FillEnum(void){                     return          600; }
    100 int WaterFillEnum(void){                return          601; }
    101 int IceFillEnum(void){                  return          602; }
    102 int AirFillEnum(void){                  return          603; }
     98int WaterEnum(void){                    return          601; }
     99int IceEnum(void){                      return          602; }
     100int AirEnum(void){                      return          603; }
     101int MelangeEnum(void){                  return          604; }
    103102
    104103/*functions on enums: */

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

    r1741 r1784  
    9393/*GEOGRAPHY:*/
    9494int GeographyEnum(void);
    95 int WaterEnum(void);
    9695int IceSheetEnum(void);
    9796int IceShelfEnum(void);
    9897
    9998/*FILL: */
    100 int FillEnum(void);
    101 int WaterFillEnum(void);
    102 int IceFillEnum(void);
    103 int AirFillEnum(void);
     99int WaterEnum(void);
     100int IceEnum(void);
     101int AirEnum(void);
     102int MelangeEnum(void);
    104103
    105104/*Functions on enums: */

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

    r1671 r1784  
    604604bin_PROGRAMS =
    605605else
    606 bin_PROGRAMS = diagnostic.exe  control.exe thermal.exe prognostic.exe transient.exe
     606dnl bin_PROGRAMS = diagnostic.exe  control.exe thermal.exe prognostic.exe transient.exe
    607607endif
    608608

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

    r1734 r1784  
    183183
    184184        for(i=0;i<model->numrifts;i++){
    185                 el1=(int)*(model->riftinfo+9*i+2)-1; //matlab indexing to c indexing
    186                 el2=(int)*(model->riftinfo+9*i+3)-1; //matlab indexing to c indexing
     185                el1=(int)*(model->riftinfo+RIFTINFOSIZE*i+2)-1; //matlab indexing to c indexing
     186                el2=(int)*(model->riftinfo+RIFTINFOSIZE*i+3)-1; //matlab indexing to c indexing
    187187                epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids;
    188188        }

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

    r1762 r1784  
    5454        int    riftfront_fill;
    5555        double riftfront_friction;
     56        double riftfront_fraction;
    5657        bool   riftfront_shelf;
    5758        double riftfront_penalty_offset;
    … …  
    6667        int    fill;
    6768        double friction;
     69        double fraction;
    6870       
    6971        /*Create loads: */
    … …  
    158160        ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat");
    159161        ModelFetchData((void**)&model->gridoniceshelf,NULL,NULL,model_handle,"gridoniceshelf","Matrix","Mat");
    160 
     162       
    161163        if(model->numrifts){
    162164                for(i=0;i<model->numrifts;i++){
    163165                               
    164                         el1=(int)*(model->riftinfo+9*i+2);
     166                        el1=(int)*(model->riftinfo+RIFTINFOSIZE*i+2);
    165167                        #ifdef _PARALLEL_
    166168                        if (model->epart[el1-1]!=my_rank){
    … …  
    172174
    173175                        /*Ok, retrieve all the data needed to add a penalty between the two grids: */
    174                         el2=(int)*(model->riftinfo+9*i+3);
    175 
    176                         grid1=(int)*(model->riftinfo+9*i+0);
    177                         grid2=(int)*(model->riftinfo+9*i+1);
    178                        
    179                         normal[0]=*(model->riftinfo+9*i+4);
    180                         normal[1]=*(model->riftinfo+9*i+5);
    181                         length=*(model->riftinfo+9*i+6);
    182                        
    183                         fill = (int)*(model->riftinfo+9*i+7);
    184                         friction=*(model->riftinfo+9*i+8);
     176                        el2=(int)*(model->riftinfo+RIFTINFOSIZE*i+3);
     177
     178                        grid1=(int)*(model->riftinfo+RIFTINFOSIZE*i+0);
     179                        grid2=(int)*(model->riftinfo+RIFTINFOSIZE*i+1);
     180                       
     181                        normal[0]=*(model->riftinfo+RIFTINFOSIZE*i+4);
     182                        normal[1]=*(model->riftinfo+RIFTINFOSIZE*i+5);
     183                        length=*(model->riftinfo+RIFTINFOSIZE*i+6);
     184                       
     185                        fill = (int)*(model->riftinfo+RIFTINFOSIZE*i+7);
     186                        friction=*(model->riftinfo+RIFTINFOSIZE*i+8);
     187                        fraction=*(model->riftinfo+RIFTINFOSIZE*i+9);
    185188       
    186189                        strcpy(riftfront_type,"2d");
    … …  
    205208                        riftfront_fill=fill;
    206209                        riftfront_friction=friction;
     210                        riftfront_fraction=fraction;
    207211                        riftfront_shelf=(bool)model->gridoniceshelf[grid1-1];
    208212
    … …  
    214218                        riftfront_prestable=0;
    215219                       
    216                         riftfront=new Riftfront(riftfront_type,riftfront_id, riftfront_node_ids, riftfront_mparid, riftfront_h,riftfront_b,riftfront_s,riftfront_normal,riftfront_length,riftfront_fill,riftfront_friction, riftfront_penalty_offset, riftfront_penalty_lock, riftfront_active,riftfront_counter,riftfront_prestable,riftfront_shelf);
     220                        riftfront=new Riftfront(riftfront_type,riftfront_id, riftfront_node_ids, riftfront_mparid, riftfront_h,riftfront_b,riftfront_s,riftfront_normal,riftfront_length,riftfront_fill,riftfront_friction, riftfront_fraction, riftfront_penalty_offset, riftfront_penalty_lock, riftfront_active,riftfront_counter,riftfront_prestable,riftfront_shelf);
    217221
    218222                        loads->AddObject(riftfront);
    219223                }
    220224        }
     225
    221226
    222227        /*free ressources: */
    … …  
    227232        xfree((void**)&model->gridoniceshelf);
    228233
     234
    229235        cleanup_and_return:
    230236

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

    r1762 r1784  
    1010#include "../DataSet/DataSet.h"
    1111#include "../toolkits/toolkits.h"
     12
     13#define RIFTINFOSIZE 10
    1214
    1315struct Model {

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

    r1628 r1784  
    3232         * management routine, otherwise, skip : */
    3333        if (RiftIsPresent(loads)){
    34                 RiftConstraints(&num_unstable_constraints,loads,inputs,analysis_type,sub_analysis_type);
     34                RiftConstraints(&converged,&num_unstable_constraints,loads,inputs,analysis_type,sub_analysis_type);
    3535        }
    3636        else if(loads->MeltingIsPresent()){
    37                 loads->MeltingConstraints(&num_unstable_constraints,inputs,analysis_type,sub_analysis_type);
     37                loads->MeltingConstraints(&converged,&num_unstable_constraints,inputs,analysis_type,sub_analysis_type);
    3838        }
    3939        else{
    … …  
    4242        }
    4343
    44         /*Have we converged? : */
    45         if (num_unstable_constraints==0) converged=1;
    46         else converged=0;
    47        
     44               
    4845        /*Assign output pointers: */
    4946        *pconverged=converged;
    5047        *pnum_unstable_constraints=num_unstable_constraints;
    51        
    5248}

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

    r1628 r1784  
    66#include "../EnumDefinitions/EnumDefinitions.h"
    77
     8#define _ZIGZAGCOUNTER_
     9
    810#undef __FUNCT__
    911#define __FUNCT__ "RiftConstrain"
    10 int RiftConstraints(int* pnum_unstable_constraints,DataSet* loads,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
     12int RiftConstraints(int* pconverged, int* pnum_unstable_constraints,DataSet* loads,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
    1113
    1214        int num_unstable_constraints=0;
     15        int converged=0;
    1316        int potential;
    1417
    15         /*First, is set of constraints pre-stabilised? Ie: all unstable rifts are penalised? :*/
    16         if(!IsPreStable(loads)){
    17 
    18                 PreConstrain(&num_unstable_constraints,loads,inputs,analysis_type);
    19 
    20                 if(!num_unstable_constraints){
    21                         /*Ok, preconstraining is done. Figure out size of unstable constraints pool. Then set
    22                          * prestable flag to 1 for all rift grids: */
    23 
    24                         potential=PotentialUnstableConstraints(loads,inputs,analysis_type);
    25                         printf("      set of constraints is prestabilised; unstable constraints potentialnumber: %i\n",potential);
    26 
    27                         SetPreStable(loads);
    28                         num_unstable_constraints=-1; //ensure convergence does not happen!
    29                 }
    30 
     18        /*First, has the material non-linearity stabilized? : */
     19        if(!IsMaterialStable(loads,inputs,analysis_type)){
     20                /*Do nothing, no constraints activated!: */
     21                printf("converging material\n");
     22                num_unstable_constraints=0;
     23                converged=0;
    3124        }
    3225        else{
    33                 /*Ok, set of constraints is stable. Start increasing number of inactive constraints, one by one: */
    34        
    35                 /*Figure out max penetration: */
    36                 MaxPenetrationInInputs(loads,inputs,analysis_type);
     26                /*Ok, start constraining: */
     27                printf("converged material: constraining\n");
     28
     29                #ifdef _ZIGZAGCOUNTER_
     30                if(!IsPreStable(loads)){
     31
     32                        PreConstrain(&num_unstable_constraints,loads,inputs,analysis_type);
     33
     34                        if(!num_unstable_constraints){
     35                                /*Ok, preconstraining is done. Figure out size of unstable constraints pool. Then set
     36                                 * prestable flag to 1 for all rift grids: */
     37
     38                                potential=PotentialUnstableConstraints(loads,inputs,analysis_type);
     39                                printf("      set of constraints is prestabilised; unstable constraints potentialnumber: %i\n",potential);
     40
     41                                SetPreStable(loads);
     42                                num_unstable_constraints=-1; //ensure convergence does not happen!
     43                        }
     44
     45                }
     46                else{
     47                        /*Ok, set of constraints is stable. Start increasing number of inactive constraints, one by one: */
    3748               
    38                 /*Constrain rifts: */
    39                 Constrain(&num_unstable_constraints,loads,inputs,analysis_type);
     49                        /*Figure out max penetration: */
     50                        MaxPenetrationInInputs(loads,inputs,analysis_type);
     51                       
     52                        /*Constrain rifts: */
     53                        Constrain(&num_unstable_constraints,loads,inputs,analysis_type);
     54                }
     55                #else
     56                        Constrain(&num_unstable_constraints,loads,inputs,analysis_type);
     57                #endif
     58
     59                if(num_unstable_constraints==0)converged=1;
    4060        }
    4161
    4262        /*Assign output pointers: */
     63        *pconverged=converged;
    4364        *pnum_unstable_constraints=num_unstable_constraints;
    44 
     65}
     66
     67#undef __FUNCT__
     68#define __FUNCT__ "IsMaterialStable"
     69int IsMaterialStable(DataSet* loads,ParameterInputs* inputs,int analysis_type){
     70
     71        int i;
     72       
     73        Riftfront* riftfront=NULL;
     74        int found=0;
     75        int mpi_found=0;
     76
     77        /*go though loads, and if non-linearity of the material has converged, let all penalties know: */
     78        for (i=0;i<loads->Size();i++){
     79
     80                if(RiftfrontEnum()==loads->GetEnum(i)){
     81
     82                        riftfront=(Riftfront*)loads->GetObjectByOffset(i);
     83
     84                        if (riftfront->IsMaterialStable(inputs,analysis_type)){
     85                                found=1;
     86                                /*do not break! all penalties should get informed the non-linearity converged!*/
     87                        }
     88                }
     89        }
     90
     91        #ifdef _PARALLEL_
     92        MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
     93        MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD);               
     94        found=mpi_found;
     95        #endif
     96
     97        return found;
    4598}
    4699
    … …  
    239292                        riftfront=(Riftfront*)loads->GetObjectByOffset(i);
    240293
    241                         riftfront->Penetration(&penetration,inputs,analysis_type);
     294                        riftfront->MaxPenetration(&penetration,inputs,analysis_type);
    242295
    243296                        if (penetration>max_penetration)max_penetration=penetration;

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

    r1628 r1784  
    1010#include "../DataSet/DataSet.h"
    1111
    12 int RiftConstraints(int* pnum_unstable_constraints,DataSet* loads,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
     12int RiftConstraints(int* pconverged, int* pnum_unstable_constraints,DataSet* loads,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
    1313
    1414int RiftIsPresent(DataSet* loads);
    … …  
    2525
    2626int PotentialUnstableConstraints(DataSet* loads,ParameterInputs* inputs,int analysis_type_enum);
     27
     28int IsMaterialStable(DataSet* loads,ParameterInputs* inputs,int analysis_type_enum);
    2729#endif

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

    r1728 r1784  
    366366        /*Recover material and fill parameters: */
    367367        matpar=(Matpar*)element->GetMatPar();
    368         if (element->GetShelf())fill=WaterFillEnum();
    369         else fill=AirFillEnum();
     368        if (element->GetShelf())fill=WaterEnum();
     369        else fill=AirEnum();
    370370
    371371        //check that the element is onbed (collapsed formulation) otherwise:pe=0
    … …  
    494494        /*Recover material and fill parameters: */
    495495        matpar=(Matpar*)element->GetMatPar();
    496         if (element->GetShelf())fill=WaterFillEnum();
    497         else fill=AirFillEnum();
     496        if (element->GetShelf())fill=WaterEnum();
     497        else fill=AirEnum();
    498498
    499499        /* Set pe_g to 0: */
    … …  
    657657        /*Recover material and fill parameters: */
    658658        matpar=(Matpar*)element->GetMatPar();
    659         if (element->GetShelf())fill=WaterFillEnum();
    660         else fill=AirFillEnum();
     659        if (element->GetShelf())fill=WaterEnum();
     660        else fill=AirEnum();
    661661
    662662        /* Set pe_g to 0: */
    … …  
    868868                bed=b1*(1+segment_gauss_coord[ig])/2+b2*(1-segment_gauss_coord[ig])/2;
    869869
    870                 if (fill==WaterFillEnum()){
     870                if (fill==WaterEnum()){
    871871                        //icefront ends in water:
    872872                        ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
    … …  
    878878                        water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2));
    879879                }
    880                 else if (fill==AirFillEnum()){
     880                else if (fill==AirEnum()){
    881881                        ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
    882882                        air_pressure=0;
    … …  
    11321132
    11331133                        //Now deal with water pressure:
    1134                         if(fill==WaterFillEnum()){ //icefront ends in water
     1134                        if(fill==WaterEnum()){ //icefront ends in water
    11351135                                water_level_above_g_tria=min(0,z_g[i]);//0 if the gaussian point is above water level
    11361136                                water_pressure_tria=rho_water*gravity*water_level_above_g_tria;
    11371137                        }
    1138                         else if(fill==AirFillEnum()){
     1138                        else if(fill==AirEnum()){
    11391139                                water_pressure_tria=0;
    11401140                        }
    … …  
    13811381
    13821382                        //Now deal with water pressure:
    1383                         if(fill==WaterFillEnum()){ //icefront ends in water
     1383                        if(fill==WaterEnum()){ //icefront ends in water
    13841384                                water_level_above_g_tria=min(0,z_g[i]);//0 if the gaussian point is above water level
    13851385                                water_pressure_tria=rho_water*gravity*water_level_above_g_tria;
    13861386                        }
    1387                         else if(fill==AirFillEnum()){
     1387                        else if(fill==AirEnum()){
    13881388                                water_pressure_tria=0;
    13891389                        }

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

    r1735 r1784  
    2020               
    2121Riftfront::Riftfront(){
     22        /*in case :*/
     23        material_converged=0;
    2224        return;
    2325}
    2426
    25 Riftfront::Riftfront(char riftfront_type[RIFTFRONTSTRING],int riftfront_id, int riftfront_node_ids[MAX_RIFTFRONT_GRIDS], int riftfront_mparid, double riftfront_h[MAX_RIFTFRONT_GRIDS],double riftfront_b[MAX_RIFTFRONT_GRIDS],double riftfront_s[MAX_RIFTFRONT_GRIDS],double riftfront_normal[2],double riftfront_length,int riftfront_fill,double riftfront_friction, double riftfront_penalty_offset, int riftfront_penalty_lock, bool riftfront_active,int riftfront_counter,bool riftfront_prestable,bool riftfront_shelf){
     27Riftfront::Riftfront(char riftfront_type[RIFTFRONTSTRING],int riftfront_id, int riftfront_node_ids[MAX_RIFTFRONT_GRIDS], int riftfront_mparid, double riftfront_h[MAX_RIFTFRONT_GRIDS],double riftfront_b[MAX_RIFTFRONT_GRIDS],double riftfront_s[MAX_RIFTFRONT_GRIDS],double riftfront_normal[2],double riftfront_length,int riftfront_fill,double riftfront_friction, double riftfront_fraction,double riftfront_penalty_offset, int riftfront_penalty_lock, bool riftfront_active,int riftfront_counter,bool riftfront_prestable,bool riftfront_shelf){
    2628
    2729        int i;
    … …  
    5153        fill=riftfront_fill;
    5254        friction=riftfront_friction;
     55        fraction=riftfront_fraction;
    5356        penalty_offset=riftfront_penalty_offset;
    5457        penalty_lock=riftfront_penalty_lock;
    … …  
    5861        shelf=riftfront_shelf;
    5962
     63        /*not in the constructor, but still needed: */
     64        material_converged=0;
     65
    6066        return;
    6167}
    … …  
    7985
    8086        printf("normal [%g,%g], length %g\n",normal[0],normal[1],normal[2]);
    81         printf("fill: %i friction %g\n",fill,friction);
     87        printf("fill: %i friction %g fraction %g\n",fill,friction,fraction);
    8288        printf("penalty_offset %g\n",penalty_offset);
    8389        printf("penalty_lock %i\n",penalty_lock);
    … …  
    8591        printf("counter %i\n",counter);
    8692        printf("prestable %i\n",prestable);
     93        printf("material_converged %i\n",material_converged);
    8794        printf("shelf %i\n",shelf);
    8895}
    … …  
    107114
    108115        printf("normal [%g,%g], length %g\n",normal[0],normal[1],normal[2]);
    109         printf("fill: %i friction %g\n",fill,friction);
     116        printf("fill: %i friction %g fraction %g\n",fill,friction,fraction);
    110117        printf("penalty_offset %g\n",penalty_offset);
    111118        printf("penalty_lock %i\n",penalty_lock);
    … …  
    113120        printf("counter %i\n",counter);
    114121        printf("prestable %i\n",prestable);
     122        printf("material_converged %i\n",material_converged);
    115123        printf("shelf %i\n",shelf);
    116124       
    … …  
    147155        memcpy(marshalled_dataset,&fill,sizeof(fill));marshalled_dataset+=sizeof(fill);
    148156        memcpy(marshalled_dataset,&friction,sizeof(friction));marshalled_dataset+=sizeof(friction);
     157        memcpy(marshalled_dataset,&fraction,sizeof(fraction));marshalled_dataset+=sizeof(fraction);
    149158        memcpy(marshalled_dataset,&penalty_offset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset);
    150159        memcpy(marshalled_dataset,&penalty_lock,sizeof(penalty_lock));marshalled_dataset+=sizeof(penalty_lock);
    … …  
    152161        memcpy(marshalled_dataset,&counter,sizeof(counter));marshalled_dataset+=sizeof(counter);
    153162        memcpy(marshalled_dataset,&prestable,sizeof(prestable));marshalled_dataset+=sizeof(prestable);
     163        memcpy(marshalled_dataset,&material_converged,sizeof(material_converged));marshalled_dataset+=sizeof(material_converged);
    154164        memcpy(marshalled_dataset,&shelf,sizeof(shelf));marshalled_dataset+=sizeof(shelf);
    155165
    … …  
    173183                sizeof(fill)+
    174184                sizeof(friction)+
     185                sizeof(fraction)+
    175186                sizeof(penalty_offset)+
    176187                sizeof(penalty_lock)+
    … …  
    178189                sizeof(counter)+
    179190                sizeof(prestable)+
     191                sizeof(material_converged)+
    180192                sizeof(shelf)+
    181193                sizeof(int); //sizeof(int) for enum type
    … …  
    212224        memcpy(&fill,marshalled_dataset,sizeof(fill));marshalled_dataset+=sizeof(fill);
    213225        memcpy(&friction,marshalled_dataset,sizeof(friction));marshalled_dataset+=sizeof(friction);
     226        memcpy(&fraction,marshalled_dataset,sizeof(fraction));marshalled_dataset+=sizeof(fraction);
    214227        memcpy(&penalty_offset,marshalled_dataset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset);
    215228        memcpy(&penalty_lock,marshalled_dataset,sizeof(penalty_lock));marshalled_dataset+=sizeof(penalty_lock);
    … …  
    217230        memcpy(&counter,marshalled_dataset,sizeof(counter));marshalled_dataset+=sizeof(counter);
    218231        memcpy(&prestable,marshalled_dataset,sizeof(prestable));marshalled_dataset+=sizeof(prestable);
     232        memcpy(&material_converged,marshalled_dataset,sizeof(material_converged));marshalled_dataset+=sizeof(material_converged);
    219233        memcpy(&shelf,marshalled_dataset,sizeof(shelf));marshalled_dataset+=sizeof(shelf);
    220234
    … …  
    410424        double       bed;
    411425        double       pressure;
     426        double       pressure_litho;
     427        double       pressure_air;
     428        double       pressure_melange;
     429        double       pressure_water;
    412430       
    413431        /*Some pointer intialization: */
    … …  
    445463
    446464                /*Ok, now compute the pressure (in norm) that is being applied to the flanks, depending on the type of fill: */
    447                 if(fill==WaterFillEnum()){
     465                if(fill==WaterEnum()){
    448466                        if(shelf){
    449467                                /*We are on an ice shelf, hydrostatic equilibrium is used to determine the pressure for water fill: */
    … …  
    455473                        }
    456474                }
    457                 else if(fill==AirFillEnum()){
     475                else if(fill==AirEnum()){
    458476                        pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2;   //icefront on an ice sheet, pressure imbalance ice vs air.
    459477                }
    460                 else if(fill==IceFillEnum()){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice)
     478                else if(fill==IceEnum()){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice)
    461479                        pressure=0;
     480                }
     481                else if(fill==IceEnum()){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice)
     482                        pressure=0;
     483                }
     484                else if(fill==MelangeEnum()){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice)
     485                       
     486                        if(!shelf) throw ErrorException(__FUNCT__,exprintf("%s%s%i%s",__FUNCT__," error message: fill type ",fill," not supported on ice sheets yet."));
     487
     488                        pressure_litho=rho_ice*gravity*pow(thickness,(double)2)/(double)2;
     489                        pressure_air=0;
     490                        pressure_melange=rho_ice*gravity*pow(fraction*thickness,(double)2)/(double)2;
     491                        pressure_water=1.0/2.0*rho_water*gravity*  ( pow(bed,2.0)-pow(rho_ice/rho_water*fraction*thickness,2.0) );
     492
     493                        pressure=pressure_litho-pressure_air-pressure_melange-pressure_water;
    462494                }
    463495                else{
    … …  
    559591        *punstable=unstable;
    560592}
     593
     594
     595#define _ZIGZAGCOUNTER_
    561596
    562597#undef __FUNCT__
    … …  
    582617        if(!found)throw ErrorException(__FUNCT__," could not find velocity in inputs!");
    583618
     619
     620        /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
     621        penetration=(vxvy_list[1][0]-vxvy_list[0][0])*normal[0]+(vxvy_list[1][1]-vxvy_list[0][1])*normal[1];
     622
     623        #ifdef _ZIGZAGCOUNTER_
     624
    584625        found=inputs->Recover("max_penetration",&max_penetration);
    585626        if(!found)throw ErrorException(__FUNCT__," could not find max_penetration in inputs!");
    586 
    587         /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
    588         penetration=(vxvy_list[1][0]-vxvy_list[0][0])*normal[0]+(vxvy_list[1][1]-vxvy_list[0][1])*normal[1];
    589627
    590628        /*Activate or deactivate penalties: */
    … …  
    624662                }
    625663        }
     664        #else
     665        if(penetration<0)activate=1;
     666        else  activate=0;
     667        #endif
    626668
    627669        //Figure out stability of this penalty
    … …  
    636678        //Set penalty flag
    637679        this->active=activate;
     680
     681        //if ((penetration>0) & (this->active==1))printf("Riftfront %i wants to be released\n",GetId());
    638682
    639683        /*assign output pointer: */
    … …  
    665709        /*Now, we return penetration only if we are active!: */
    666710        if(this->active==0)penetration=0;
     711       
     712        /*assign output pointer: */
     713        *ppenetration=penetration;
     714
     715}
     716
     717#undef __FUNCT__
     718#define __FUNCT__ "Riftfront::MaxPenetration"
     719int   Riftfront::MaxPenetration(double* ppenetration, void* vinputs, int analysis_type){
     720
     721        const int     numgrids=2;
     722        int           dofs[2]={0,1};
     723        double        vxvy_list[2][2]; //velocities for all grids
     724        double        max_penetration;
     725        double        penetration=0;
     726        int           found;
     727
     728        ParameterInputs* inputs=NULL;
     729
     730        inputs=(ParameterInputs*)vinputs;
     731
     732        //initialize:
     733        penetration=-1;
     734
     735        found=inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
     736        if(!found)throw ErrorException(__FUNCT__," could not find velocity in inputs!");
     737
     738        /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
     739        penetration=(vxvy_list[1][0]-vxvy_list[0][0])*normal[0]+(vxvy_list[1][1]-vxvy_list[0][1])*normal[1];
     740
     741        /*Now, we return penetration only if we are active!: */
     742        if(this->active==0)penetration=-1;
     743
     744        /*If we are zigzag locked, same thing: */
     745        if(this->counter>this->penalty_lock)penetration=-1;
    667746       
    668747        /*assign output pointer: */
    … …  
    711790        *punstable=unstable;
    712791}
     792
     793
     794#undef __FUNCT__
     795#define __FUNCT__ "Riftfront::IsMaterialStable"
     796int   Riftfront::IsMaterialStable(void* vinputs, int analysis_type){
     797
     798        int found=0;
     799        ParameterInputs* inputs=NULL;
     800        double converged=0;
     801
     802        inputs=(ParameterInputs*)vinputs;
     803
     804        found=inputs->Recover("converged",&converged);
     805        if(!found)throw ErrorException(__FUNCT__," could not find converged flag  in inputs!");
     806
     807        if(converged){
     808                /*ok, material non-linearity has converged. If that was already the case, we keep
     809                 * constraining the rift front. If it was not, and this is the first time the material
     810                 * has converged, we start constraining now!: */
     811                this->material_converged=1;
     812        }
     813
     814        return this->material_converged;
     815}

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

    r1717 r1784  
    4141                int         fill;
    4242                double      friction;
     43                double      fraction;
    4344                bool        shelf;
    4445
    … …  
    5051                int         counter;
    5152                bool        prestable;
     53                bool        material_converged;
    5254
    5355
    … …  
    5557
    5658                Riftfront();
    57                 Riftfront(char type[RIFTFRONTSTRING],int id, int node_ids[MAX_RIFTFRONT_GRIDS], int mparid, double h[MAX_RIFTFRONT_GRIDS],double b[MAX_RIFTFRONT_GRIDS],double s[MAX_RIFTFRONT_GRIDS],double normal[2],double length,int fill,double friction, double penalty_offset, int penalty_lock,bool active,int counter,bool prestable,bool shelf);
     59                Riftfront(char type[RIFTFRONTSTRING],int id, int node_ids[MAX_RIFTFRONT_GRIDS], int mparid, double h[MAX_RIFTFRONT_GRIDS],double b[MAX_RIFTFRONT_GRIDS],double s[MAX_RIFTFRONT_GRIDS],double normal[2],double length,int fill,double friction, double fraction, double penalty_offset, int penalty_lock,bool active,int counter,bool prestable,bool shelf);
    5860                ~Riftfront();
    5961
    … …  
    7981                int   Constrain(int* punstable, void* inputs, int analysis_type);
    8082                int   Penetration(double* ppenetration, void* inputs, int analysis_type);
     83                int   MaxPenetration(double* ppenetration, void* inputs, int analysis_type);
    8184                int   PotentialUnstableConstraint(int* punstable, void* inputs, int analysis_type);
     85                int   IsMaterialStable(void* inputs, int analysis_type);
    8286                Object* copy();
    8387};

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

    r1648 r1784  
    7373        DeleteModel(&model);
    7474
    75 
    7675        /*Assign output pointers:*/
    7776        femmodel->elements=elements;

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

    r1669 r1784  
    112112                PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,loads,fem->materials,inputs,analysis_type,sub_analysis_type);
    113113
     114                //if(debug)_printf_("   number of unstable constraints: %i\n",num_unstable_constraints);
     115                _printf_("   number of unstable constraints: %i\n",num_unstable_constraints);
     116
    114117                /*Figure out if convergence is reached.*/
    115118                convergence(&converged,Kff,pf,uf,old_uf,fem->parameters);
    116119                MatFree(&Kff);VecFree(&pf);
     120               
     121                /*add converged to inputs: */
     122                inputs->Add("converged",converged);
    117123
    118124                //rift convergence

  • issm/trunk/src/m/classes/public/plot/plot_riftrelvel.m

    r1 r1784  
    6363                %plot relative velocities
    6464                h=quiver(md.x(penaltypairs(:,1)),md.y(penaltypairs(:,1)),md.vx(penaltypairs(:,2))-md.vx(penaltypairs(:,1)),md.vy(penaltypairs(:,2))-md.vy(penaltypairs(:,1)));
    65                 set(h,'Color',[1 0 0]);
     65                set(h,'Color',[0 1 0]);
    6666        end
    6767        if isp1 & isp2

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

    r1717 r1784  
    1717end
    1818
    19 md.riftinfo=zeros(numpairs,9); % 2 for grids + 2 for elements+ 2 for  normals + 1 for length + 1 for fill + 1 for friction.
     19md.riftinfo=zeros(numpairs,10); % 2 for grids + 2 for elements+ 2 for  normals + 1 for length + 1 for fill + 1 for friction + 1 for fraction.
    2020
    2121count=1;
    … …  
    2525        md.riftinfo(count:count+numpairsforthisrift-1,8)=md.rifts(i).fill;
    2626        md.riftinfo(count:count+numpairsforthisrift-1,9)=md.rifts(i).friction;
     27        md.riftinfo(count:count+numpairsforthisrift-1,10)=md.rifts(i).fraction;
    2728        count=count+numpairsforthisrift;
    2829end

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

    r1667 r1784  
    1212        soln(count).u_g=get(inputs,'velocity',[1 1 (m.parameters.numberofdofspernode>=3) (m.parameters.numberofdofspernode>=3)]); %only keep vz if running with more than 3 dofs per node
    1313        soln(count).u_f= Reducevectorgtof(soln(count).u_g,m.nodesets);
     14               
     15        %add converged=0 flag first in inputs.
     16        inputs=add(inputs,'converged',converged,'double');
    1417
    1518        displaystring(m.parameters.debug,'\n%s',['   starting direct shooting method']);
    … …  
    5356                %penalty constraints
    5457                [loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints( m.elements,m.nodes, loads, m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
     58       
     59                %displaystring(m.parameters.debug,'%s%i','      number of unstable constraints: ',num_unstable_constraints);
     60                displaystring(1,'%s%i','      number of unstable constraints: ',num_unstable_constraints);
    5561
    5662                %Figure out if convergence have been reached
    5763                converged=convergence(K_ff,p_f,soln(count).u_f,soln(count-1).u_f,m.parameters);
     64                       
     65                %add convergence status into inputs
     66                inputs=add(inputs,'converged',converged,'double');
    5867
    5968                %rift convergence criterion
    … …  
    6170                        converged=0;
    6271                end
     72
    6373        end
    6474                       

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

    r1743 r1784  
    2727        /*empty rifts structure: */
    2828        double* pNaN=NULL;
    29         const   char*   fnames[7];
     29        const   char*   fnames[8];
    3030        mwSize     ndim=2;
    3131        mwSize          dimensions[2] = {1,1};
    … …  
    239239                fnames[5] = "fill";
    240240                fnames[6] = "friction";
     241                fnames[7] = "fraction";
    241242
    242243                dimensions[0]=out_numrifts;
    243244
    244                 pmxa_array=mxCreateStructArray( ndim,dimensions,7,fnames);
     245                pmxa_array=mxCreateStructArray( ndim,dimensions,8,fnames);
    245246               
    246247                for (i=0;i<out_numrifts;i++){
    … …  
    283284                        mxSetField(pmxa_array,i,"penaltypairs",pmxa_array3);
    284285
    285                         /*Friction and fill: set to 0 both */
     286                        /*Friction fraction and fill: set to 0 both */
    286287                        mxSetField(pmxa_array,i,"friction",mxCreateDoubleScalar(0));
    287                         mxSetField(pmxa_array,i,"fill",mxCreateDoubleScalar(WaterFillEnum())); //default is water
    288 
    289 
    290 
     288                        mxSetField(pmxa_array,i,"fill",mxCreateDoubleScalar(IceEnum())); //default is ice
     289                        mxSetField(pmxa_array,i,"fraction",mxCreateDoubleScalar(0)); //default is ice
    291290                }
    292291        }
Note: See TracChangeset for help on using the changeset viewer.