Index: /issm/trunk/src/c/DataSet/DataSet.cpp =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.cpp (revision 1783) +++ /issm/trunk/src/c/DataSet/DataSet.cpp (revision 1784) @@ -1146,12 +1146,4 @@ #undef __FUNCT__ -#define __FUNCT__ "RiftConstraints" -void DataSet::RiftConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type){ - - throw ErrorException(__FUNCT__," not implemented yet!"); - -} - -#undef __FUNCT__ #define __FUNCT__ "MeltingIsPresent" int DataSet::MeltingIsPresent(){ @@ -1181,5 +1173,5 @@ #undef __FUNCT__ #define __FUNCT__ "MeltingConstraints" -void DataSet::MeltingConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type){ +void DataSet::MeltingConstraints(int* pconverged, int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type){ /* generic object pointer: */ @@ -1189,4 +1181,5 @@ int unstable=0; int num_unstable_constraints=0; + int converged=0; int sum_num_unstable_constraints=0; @@ -1212,5 +1205,10 @@ #endif + /*Have we converged? : */ + if (num_unstable_constraints==0) converged=1; + else converged=0; + /*Assign output pointers: */ + *pconverged=converged; *pnum_unstable_constraints=num_unstable_constraints; } Index: /issm/trunk/src/c/DataSet/DataSet.h =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.h (revision 1783) +++ /issm/trunk/src/c/DataSet/DataSet.h (revision 1784) @@ -70,7 +70,6 @@ void PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type); int RiftIsPresent(); - void RiftConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type); int MeltingIsPresent(); - void MeltingConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type); + void MeltingConstraints(int* pconverged, int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type); DataSet* Copy(void); void Du(Vec du_g,void* inputs,int analysis_type,int sub_analysis_type); Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp (revision 1783) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp (revision 1784) @@ -92,13 +92,12 @@ /*GEOGRAPHY:*/ int GeographyEnum(void){ return 500; } -int WaterEnum(void){ return 501; } int IceSheetEnum(void){ return 502; } int IceShelfEnum(void){ return 502; } /*FILL:*/ -int FillEnum(void){ return 600; } -int WaterFillEnum(void){ return 601; } -int IceFillEnum(void){ return 602; } -int AirFillEnum(void){ return 603; } +int WaterEnum(void){ return 601; } +int IceEnum(void){ return 602; } +int AirEnum(void){ return 603; } +int MelangeEnum(void){ return 604; } /*functions on enums: */ Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 1783) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 1784) @@ -93,13 +93,12 @@ /*GEOGRAPHY:*/ int GeographyEnum(void); -int WaterEnum(void); int IceSheetEnum(void); int IceShelfEnum(void); /*FILL: */ -int FillEnum(void); -int WaterFillEnum(void); -int IceFillEnum(void); -int AirFillEnum(void); +int WaterEnum(void); +int IceEnum(void); +int AirEnum(void); +int MelangeEnum(void); /*Functions on enums: */ Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 1783) +++ /issm/trunk/src/c/Makefile.am (revision 1784) @@ -604,5 +604,5 @@ bin_PROGRAMS = else -bin_PROGRAMS = diagnostic.exe control.exe thermal.exe prognostic.exe transient.exe +dnl bin_PROGRAMS = diagnostic.exe control.exe thermal.exe prognostic.exe transient.exe endif Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 1783) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 1784) @@ -183,6 +183,6 @@ for(i=0;inumrifts;i++){ - el1=(int)*(model->riftinfo+9*i+2)-1; //matlab indexing to c indexing - el2=(int)*(model->riftinfo+9*i+3)-1; //matlab indexing to c indexing + el1=(int)*(model->riftinfo+RIFTINFOSIZE*i+2)-1; //matlab indexing to c indexing + el2=(int)*(model->riftinfo+RIFTINFOSIZE*i+3)-1; //matlab indexing to c indexing epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids; } Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (revision 1783) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (revision 1784) @@ -54,4 +54,5 @@ int riftfront_fill; double riftfront_friction; + double riftfront_fraction; bool riftfront_shelf; double riftfront_penalty_offset; @@ -66,4 +67,5 @@ int fill; double friction; + double fraction; /*Create loads: */ @@ -158,9 +160,9 @@ ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat"); ModelFetchData((void**)&model->gridoniceshelf,NULL,NULL,model_handle,"gridoniceshelf","Matrix","Mat"); - + if(model->numrifts){ for(i=0;inumrifts;i++){ - el1=(int)*(model->riftinfo+9*i+2); + el1=(int)*(model->riftinfo+RIFTINFOSIZE*i+2); #ifdef _PARALLEL_ if (model->epart[el1-1]!=my_rank){ @@ -172,15 +174,16 @@ /*Ok, retrieve all the data needed to add a penalty between the two grids: */ - el2=(int)*(model->riftinfo+9*i+3); - - grid1=(int)*(model->riftinfo+9*i+0); - grid2=(int)*(model->riftinfo+9*i+1); - - normal[0]=*(model->riftinfo+9*i+4); - normal[1]=*(model->riftinfo+9*i+5); - length=*(model->riftinfo+9*i+6); - - fill = (int)*(model->riftinfo+9*i+7); - friction=*(model->riftinfo+9*i+8); + el2=(int)*(model->riftinfo+RIFTINFOSIZE*i+3); + + grid1=(int)*(model->riftinfo+RIFTINFOSIZE*i+0); + grid2=(int)*(model->riftinfo+RIFTINFOSIZE*i+1); + + normal[0]=*(model->riftinfo+RIFTINFOSIZE*i+4); + normal[1]=*(model->riftinfo+RIFTINFOSIZE*i+5); + length=*(model->riftinfo+RIFTINFOSIZE*i+6); + + fill = (int)*(model->riftinfo+RIFTINFOSIZE*i+7); + friction=*(model->riftinfo+RIFTINFOSIZE*i+8); + fraction=*(model->riftinfo+RIFTINFOSIZE*i+9); strcpy(riftfront_type,"2d"); @@ -205,4 +208,5 @@ riftfront_fill=fill; riftfront_friction=friction; + riftfront_fraction=fraction; riftfront_shelf=(bool)model->gridoniceshelf[grid1-1]; @@ -214,9 +218,10 @@ riftfront_prestable=0; - 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); + 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); loads->AddObject(riftfront); } } + /*free ressources: */ @@ -227,4 +232,5 @@ xfree((void**)&model->gridoniceshelf); + cleanup_and_return: Index: /issm/trunk/src/c/ModelProcessorx/Model.h =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Model.h (revision 1783) +++ /issm/trunk/src/c/ModelProcessorx/Model.h (revision 1784) @@ -10,4 +10,6 @@ #include "../DataSet/DataSet.h" #include "../toolkits/toolkits.h" + +#define RIFTINFOSIZE 10 struct Model { Index: /issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.cpp =================================================================== --- /issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.cpp (revision 1783) +++ /issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.cpp (revision 1784) @@ -32,8 +32,8 @@ * management routine, otherwise, skip : */ if (RiftIsPresent(loads)){ - RiftConstraints(&num_unstable_constraints,loads,inputs,analysis_type,sub_analysis_type); + RiftConstraints(&converged,&num_unstable_constraints,loads,inputs,analysis_type,sub_analysis_type); } else if(loads->MeltingIsPresent()){ - loads->MeltingConstraints(&num_unstable_constraints,inputs,analysis_type,sub_analysis_type); + loads->MeltingConstraints(&converged,&num_unstable_constraints,inputs,analysis_type,sub_analysis_type); } else{ @@ -42,11 +42,7 @@ } - /*Have we converged? : */ - if (num_unstable_constraints==0) converged=1; - else converged=0; - + /*Assign output pointers: */ *pconverged=converged; *pnum_unstable_constraints=num_unstable_constraints; - } Index: /issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.cpp =================================================================== --- /issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.cpp (revision 1783) +++ /issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.cpp (revision 1784) @@ -6,41 +6,94 @@ #include "../EnumDefinitions/EnumDefinitions.h" +#define _ZIGZAGCOUNTER_ + #undef __FUNCT__ #define __FUNCT__ "RiftConstrain" -int RiftConstraints(int* pnum_unstable_constraints,DataSet* loads,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ +int RiftConstraints(int* pconverged, int* pnum_unstable_constraints,DataSet* loads,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ int num_unstable_constraints=0; + int converged=0; int potential; - /*First, is set of constraints pre-stabilised? Ie: all unstable rifts are penalised? :*/ - if(!IsPreStable(loads)){ - - PreConstrain(&num_unstable_constraints,loads,inputs,analysis_type); - - if(!num_unstable_constraints){ - /*Ok, preconstraining is done. Figure out size of unstable constraints pool. Then set - * prestable flag to 1 for all rift grids: */ - - potential=PotentialUnstableConstraints(loads,inputs,analysis_type); - printf(" set of constraints is prestabilised; unstable constraints potentialnumber: %i\n",potential); - - SetPreStable(loads); - num_unstable_constraints=-1; //ensure convergence does not happen! - } - + /*First, has the material non-linearity stabilized? : */ + if(!IsMaterialStable(loads,inputs,analysis_type)){ + /*Do nothing, no constraints activated!: */ + printf("converging material\n"); + num_unstable_constraints=0; + converged=0; } else{ - /*Ok, set of constraints is stable. Start increasing number of inactive constraints, one by one: */ - - /*Figure out max penetration: */ - MaxPenetrationInInputs(loads,inputs,analysis_type); + /*Ok, start constraining: */ + printf("converged material: constraining\n"); + + #ifdef _ZIGZAGCOUNTER_ + if(!IsPreStable(loads)){ + + PreConstrain(&num_unstable_constraints,loads,inputs,analysis_type); + + if(!num_unstable_constraints){ + /*Ok, preconstraining is done. Figure out size of unstable constraints pool. Then set + * prestable flag to 1 for all rift grids: */ + + potential=PotentialUnstableConstraints(loads,inputs,analysis_type); + printf(" set of constraints is prestabilised; unstable constraints potentialnumber: %i\n",potential); + + SetPreStable(loads); + num_unstable_constraints=-1; //ensure convergence does not happen! + } + + } + else{ + /*Ok, set of constraints is stable. Start increasing number of inactive constraints, one by one: */ - /*Constrain rifts: */ - Constrain(&num_unstable_constraints,loads,inputs,analysis_type); + /*Figure out max penetration: */ + MaxPenetrationInInputs(loads,inputs,analysis_type); + + /*Constrain rifts: */ + Constrain(&num_unstable_constraints,loads,inputs,analysis_type); + } + #else + Constrain(&num_unstable_constraints,loads,inputs,analysis_type); + #endif + + if(num_unstable_constraints==0)converged=1; } /*Assign output pointers: */ + *pconverged=converged; *pnum_unstable_constraints=num_unstable_constraints; - +} + +#undef __FUNCT__ +#define __FUNCT__ "IsMaterialStable" +int IsMaterialStable(DataSet* loads,ParameterInputs* inputs,int analysis_type){ + + int i; + + Riftfront* riftfront=NULL; + int found=0; + int mpi_found=0; + + /*go though loads, and if non-linearity of the material has converged, let all penalties know: */ + for (i=0;iSize();i++){ + + if(RiftfrontEnum()==loads->GetEnum(i)){ + + riftfront=(Riftfront*)loads->GetObjectByOffset(i); + + if (riftfront->IsMaterialStable(inputs,analysis_type)){ + found=1; + /*do not break! all penalties should get informed the non-linearity converged!*/ + } + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD); + found=mpi_found; + #endif + + return found; } @@ -239,5 +292,5 @@ riftfront=(Riftfront*)loads->GetObjectByOffset(i); - riftfront->Penetration(&penetration,inputs,analysis_type); + riftfront->MaxPenetration(&penetration,inputs,analysis_type); if (penetration>max_penetration)max_penetration=penetration; Index: /issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.h =================================================================== --- /issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.h (revision 1783) +++ /issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.h (revision 1784) @@ -10,5 +10,5 @@ #include "../DataSet/DataSet.h" -int RiftConstraints(int* pnum_unstable_constraints,DataSet* loads,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); +int RiftConstraints(int* pconverged, int* pnum_unstable_constraints,DataSet* loads,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); int RiftIsPresent(DataSet* loads); @@ -25,3 +25,5 @@ int PotentialUnstableConstraints(DataSet* loads,ParameterInputs* inputs,int analysis_type_enum); + +int IsMaterialStable(DataSet* loads,ParameterInputs* inputs,int analysis_type_enum); #endif Index: /issm/trunk/src/c/objects/Icefront.cpp =================================================================== --- /issm/trunk/src/c/objects/Icefront.cpp (revision 1783) +++ /issm/trunk/src/c/objects/Icefront.cpp (revision 1784) @@ -366,6 +366,6 @@ /*Recover material and fill parameters: */ matpar=(Matpar*)element->GetMatPar(); - if (element->GetShelf())fill=WaterFillEnum(); - else fill=AirFillEnum(); + if (element->GetShelf())fill=WaterEnum(); + else fill=AirEnum(); //check that the element is onbed (collapsed formulation) otherwise:pe=0 @@ -494,6 +494,6 @@ /*Recover material and fill parameters: */ matpar=(Matpar*)element->GetMatPar(); - if (element->GetShelf())fill=WaterFillEnum(); - else fill=AirFillEnum(); + if (element->GetShelf())fill=WaterEnum(); + else fill=AirEnum(); /* Set pe_g to 0: */ @@ -657,6 +657,6 @@ /*Recover material and fill parameters: */ matpar=(Matpar*)element->GetMatPar(); - if (element->GetShelf())fill=WaterFillEnum(); - else fill=AirFillEnum(); + if (element->GetShelf())fill=WaterEnum(); + else fill=AirEnum(); /* Set pe_g to 0: */ @@ -868,5 +868,5 @@ bed=b1*(1+segment_gauss_coord[ig])/2+b2*(1-segment_gauss_coord[ig])/2; - if (fill==WaterFillEnum()){ + if (fill==WaterEnum()){ //icefront ends in water: ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2); @@ -878,5 +878,5 @@ water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2)); } - else if (fill==AirFillEnum()){ + else if (fill==AirEnum()){ ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2); air_pressure=0; @@ -1132,9 +1132,9 @@ //Now deal with water pressure: - if(fill==WaterFillEnum()){ //icefront ends in water + if(fill==WaterEnum()){ //icefront ends in water water_level_above_g_tria=min(0,z_g[i]);//0 if the gaussian point is above water level water_pressure_tria=rho_water*gravity*water_level_above_g_tria; } - else if(fill==AirFillEnum()){ + else if(fill==AirEnum()){ water_pressure_tria=0; } @@ -1381,9 +1381,9 @@ //Now deal with water pressure: - if(fill==WaterFillEnum()){ //icefront ends in water + if(fill==WaterEnum()){ //icefront ends in water water_level_above_g_tria=min(0,z_g[i]);//0 if the gaussian point is above water level water_pressure_tria=rho_water*gravity*water_level_above_g_tria; } - else if(fill==AirFillEnum()){ + else if(fill==AirEnum()){ water_pressure_tria=0; } Index: /issm/trunk/src/c/objects/Riftfront.cpp =================================================================== --- /issm/trunk/src/c/objects/Riftfront.cpp (revision 1783) +++ /issm/trunk/src/c/objects/Riftfront.cpp (revision 1784) @@ -20,8 +20,10 @@ Riftfront::Riftfront(){ + /*in case :*/ + material_converged=0; return; } -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){ +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_fraction,double riftfront_penalty_offset, int riftfront_penalty_lock, bool riftfront_active,int riftfront_counter,bool riftfront_prestable,bool riftfront_shelf){ int i; @@ -51,4 +53,5 @@ fill=riftfront_fill; friction=riftfront_friction; + fraction=riftfront_fraction; penalty_offset=riftfront_penalty_offset; penalty_lock=riftfront_penalty_lock; @@ -58,4 +61,7 @@ shelf=riftfront_shelf; + /*not in the constructor, but still needed: */ + material_converged=0; + return; } @@ -79,5 +85,5 @@ printf("normal [%g,%g], length %g\n",normal[0],normal[1],normal[2]); - printf("fill: %i friction %g\n",fill,friction); + printf("fill: %i friction %g fraction %g\n",fill,friction,fraction); printf("penalty_offset %g\n",penalty_offset); printf("penalty_lock %i\n",penalty_lock); @@ -85,4 +91,5 @@ printf("counter %i\n",counter); printf("prestable %i\n",prestable); + printf("material_converged %i\n",material_converged); printf("shelf %i\n",shelf); } @@ -107,5 +114,5 @@ printf("normal [%g,%g], length %g\n",normal[0],normal[1],normal[2]); - printf("fill: %i friction %g\n",fill,friction); + printf("fill: %i friction %g fraction %g\n",fill,friction,fraction); printf("penalty_offset %g\n",penalty_offset); printf("penalty_lock %i\n",penalty_lock); @@ -113,4 +120,5 @@ printf("counter %i\n",counter); printf("prestable %i\n",prestable); + printf("material_converged %i\n",material_converged); printf("shelf %i\n",shelf); @@ -147,4 +155,5 @@ memcpy(marshalled_dataset,&fill,sizeof(fill));marshalled_dataset+=sizeof(fill); memcpy(marshalled_dataset,&friction,sizeof(friction));marshalled_dataset+=sizeof(friction); + memcpy(marshalled_dataset,&fraction,sizeof(fraction));marshalled_dataset+=sizeof(fraction); memcpy(marshalled_dataset,&penalty_offset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset); memcpy(marshalled_dataset,&penalty_lock,sizeof(penalty_lock));marshalled_dataset+=sizeof(penalty_lock); @@ -152,4 +161,5 @@ memcpy(marshalled_dataset,&counter,sizeof(counter));marshalled_dataset+=sizeof(counter); memcpy(marshalled_dataset,&prestable,sizeof(prestable));marshalled_dataset+=sizeof(prestable); + memcpy(marshalled_dataset,&material_converged,sizeof(material_converged));marshalled_dataset+=sizeof(material_converged); memcpy(marshalled_dataset,&shelf,sizeof(shelf));marshalled_dataset+=sizeof(shelf); @@ -173,4 +183,5 @@ sizeof(fill)+ sizeof(friction)+ + sizeof(fraction)+ sizeof(penalty_offset)+ sizeof(penalty_lock)+ @@ -178,4 +189,5 @@ sizeof(counter)+ sizeof(prestable)+ + sizeof(material_converged)+ sizeof(shelf)+ sizeof(int); //sizeof(int) for enum type @@ -212,4 +224,5 @@ memcpy(&fill,marshalled_dataset,sizeof(fill));marshalled_dataset+=sizeof(fill); memcpy(&friction,marshalled_dataset,sizeof(friction));marshalled_dataset+=sizeof(friction); + memcpy(&fraction,marshalled_dataset,sizeof(fraction));marshalled_dataset+=sizeof(fraction); memcpy(&penalty_offset,marshalled_dataset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset); memcpy(&penalty_lock,marshalled_dataset,sizeof(penalty_lock));marshalled_dataset+=sizeof(penalty_lock); @@ -217,4 +230,5 @@ memcpy(&counter,marshalled_dataset,sizeof(counter));marshalled_dataset+=sizeof(counter); memcpy(&prestable,marshalled_dataset,sizeof(prestable));marshalled_dataset+=sizeof(prestable); + memcpy(&material_converged,marshalled_dataset,sizeof(material_converged));marshalled_dataset+=sizeof(material_converged); memcpy(&shelf,marshalled_dataset,sizeof(shelf));marshalled_dataset+=sizeof(shelf); @@ -410,4 +424,8 @@ double bed; double pressure; + double pressure_litho; + double pressure_air; + double pressure_melange; + double pressure_water; /*Some pointer intialization: */ @@ -445,5 +463,5 @@ /*Ok, now compute the pressure (in norm) that is being applied to the flanks, depending on the type of fill: */ - if(fill==WaterFillEnum()){ + if(fill==WaterEnum()){ if(shelf){ /*We are on an ice shelf, hydrostatic equilibrium is used to determine the pressure for water fill: */ @@ -455,9 +473,23 @@ } } - else if(fill==AirFillEnum()){ + else if(fill==AirEnum()){ pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2; //icefront on an ice sheet, pressure imbalance ice vs air. } - else if(fill==IceFillEnum()){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice) + else if(fill==IceEnum()){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice) pressure=0; + } + else if(fill==IceEnum()){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice) + pressure=0; + } + else if(fill==MelangeEnum()){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice) + + if(!shelf) throw ErrorException(__FUNCT__,exprintf("%s%s%i%s",__FUNCT__," error message: fill type ",fill," not supported on ice sheets yet.")); + + pressure_litho=rho_ice*gravity*pow(thickness,(double)2)/(double)2; + pressure_air=0; + pressure_melange=rho_ice*gravity*pow(fraction*thickness,(double)2)/(double)2; + pressure_water=1.0/2.0*rho_water*gravity* ( pow(bed,2.0)-pow(rho_ice/rho_water*fraction*thickness,2.0) ); + + pressure=pressure_litho-pressure_air-pressure_melange-pressure_water; } else{ @@ -559,4 +591,7 @@ *punstable=unstable; } + + +#define _ZIGZAGCOUNTER_ #undef __FUNCT__ @@ -582,9 +617,12 @@ if(!found)throw ErrorException(__FUNCT__," could not find velocity in inputs!"); + + /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */ + penetration=(vxvy_list[1][0]-vxvy_list[0][0])*normal[0]+(vxvy_list[1][1]-vxvy_list[0][1])*normal[1]; + + #ifdef _ZIGZAGCOUNTER_ + found=inputs->Recover("max_penetration",&max_penetration); if(!found)throw ErrorException(__FUNCT__," could not find max_penetration in inputs!"); - - /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */ - penetration=(vxvy_list[1][0]-vxvy_list[0][0])*normal[0]+(vxvy_list[1][1]-vxvy_list[0][1])*normal[1]; /*Activate or deactivate penalties: */ @@ -624,4 +662,8 @@ } } + #else + if(penetration<0)activate=1; + else activate=0; + #endif //Figure out stability of this penalty @@ -636,4 +678,6 @@ //Set penalty flag this->active=activate; + + //if ((penetration>0) & (this->active==1))printf("Riftfront %i wants to be released\n",GetId()); /*assign output pointer: */ @@ -665,4 +709,39 @@ /*Now, we return penetration only if we are active!: */ if(this->active==0)penetration=0; + + /*assign output pointer: */ + *ppenetration=penetration; + +} + +#undef __FUNCT__ +#define __FUNCT__ "Riftfront::MaxPenetration" +int Riftfront::MaxPenetration(double* ppenetration, void* vinputs, int analysis_type){ + + const int numgrids=2; + int dofs[2]={0,1}; + double vxvy_list[2][2]; //velocities for all grids + double max_penetration; + double penetration=0; + int found; + + ParameterInputs* inputs=NULL; + + inputs=(ParameterInputs*)vinputs; + + //initialize: + penetration=-1; + + found=inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); + if(!found)throw ErrorException(__FUNCT__," could not find velocity in inputs!"); + + /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */ + penetration=(vxvy_list[1][0]-vxvy_list[0][0])*normal[0]+(vxvy_list[1][1]-vxvy_list[0][1])*normal[1]; + + /*Now, we return penetration only if we are active!: */ + if(this->active==0)penetration=-1; + + /*If we are zigzag locked, same thing: */ + if(this->counter>this->penalty_lock)penetration=-1; /*assign output pointer: */ @@ -711,2 +790,26 @@ *punstable=unstable; } + + +#undef __FUNCT__ +#define __FUNCT__ "Riftfront::IsMaterialStable" +int Riftfront::IsMaterialStable(void* vinputs, int analysis_type){ + + int found=0; + ParameterInputs* inputs=NULL; + double converged=0; + + inputs=(ParameterInputs*)vinputs; + + found=inputs->Recover("converged",&converged); + if(!found)throw ErrorException(__FUNCT__," could not find converged flag in inputs!"); + + if(converged){ + /*ok, material non-linearity has converged. If that was already the case, we keep + * constraining the rift front. If it was not, and this is the first time the material + * has converged, we start constraining now!: */ + this->material_converged=1; + } + + return this->material_converged; +} Index: /issm/trunk/src/c/objects/Riftfront.h =================================================================== --- /issm/trunk/src/c/objects/Riftfront.h (revision 1783) +++ /issm/trunk/src/c/objects/Riftfront.h (revision 1784) @@ -41,4 +41,5 @@ int fill; double friction; + double fraction; bool shelf; @@ -50,4 +51,5 @@ int counter; bool prestable; + bool material_converged; @@ -55,5 +57,5 @@ Riftfront(); - 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); + 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); ~Riftfront(); @@ -79,5 +81,7 @@ int Constrain(int* punstable, void* inputs, int analysis_type); int Penetration(double* ppenetration, void* inputs, int analysis_type); + int MaxPenetration(double* ppenetration, void* inputs, int analysis_type); int PotentialUnstableConstraint(int* punstable, void* inputs, int analysis_type); + int IsMaterialStable(void* inputs, int analysis_type); Object* copy(); }; Index: /issm/trunk/src/c/parallel/CreateFemModel.cpp =================================================================== --- /issm/trunk/src/c/parallel/CreateFemModel.cpp (revision 1783) +++ /issm/trunk/src/c/parallel/CreateFemModel.cpp (revision 1784) @@ -73,5 +73,4 @@ DeleteModel(&model); - /*Assign output pointers:*/ femmodel->elements=elements; Index: /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp (revision 1783) +++ /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp (revision 1784) @@ -112,7 +112,13 @@ PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,loads,fem->materials,inputs,analysis_type,sub_analysis_type); + //if(debug)_printf_(" number of unstable constraints: %i\n",num_unstable_constraints); + _printf_(" number of unstable constraints: %i\n",num_unstable_constraints); + /*Figure out if convergence is reached.*/ convergence(&converged,Kff,pf,uf,old_uf,fem->parameters); MatFree(&Kff);VecFree(&pf); + + /*add converged to inputs: */ + inputs->Add("converged",converged); //rift convergence Index: /issm/trunk/src/m/classes/public/plot/plot_riftrelvel.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_riftrelvel.m (revision 1783) +++ /issm/trunk/src/m/classes/public/plot/plot_riftrelvel.m (revision 1784) @@ -63,5 +63,5 @@ %plot relative velocities 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))); - set(h,'Color',[1 0 0]); + set(h,'Color',[0 1 0]); end if isp1 & isp2 Index: /issm/trunk/src/m/classes/public/presolve.m =================================================================== --- /issm/trunk/src/m/classes/public/presolve.m (revision 1783) +++ /issm/trunk/src/m/classes/public/presolve.m (revision 1784) @@ -17,5 +17,5 @@ end -md.riftinfo=zeros(numpairs,9); % 2 for grids + 2 for elements+ 2 for normals + 1 for length + 1 for fill + 1 for friction. +md.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. count=1; @@ -25,4 +25,5 @@ md.riftinfo(count:count+numpairsforthisrift-1,8)=md.rifts(i).fill; md.riftinfo(count:count+numpairsforthisrift-1,9)=md.rifts(i).friction; + md.riftinfo(count:count+numpairsforthisrift-1,10)=md.rifts(i).fraction; count=count+numpairsforthisrift; end Index: /issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m (revision 1783) +++ /issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m (revision 1784) @@ -12,4 +12,7 @@ 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 soln(count).u_f= Reducevectorgtof(soln(count).u_g,m.nodesets); + + %add converged=0 flag first in inputs. + inputs=add(inputs,'converged',converged,'double'); displaystring(m.parameters.debug,'\n%s',[' starting direct shooting method']); @@ -53,7 +56,13 @@ %penalty constraints [loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints( m.elements,m.nodes, loads, m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); + + %displaystring(m.parameters.debug,'%s%i',' number of unstable constraints: ',num_unstable_constraints); + displaystring(1,'%s%i',' number of unstable constraints: ',num_unstable_constraints); %Figure out if convergence have been reached converged=convergence(K_ff,p_f,soln(count).u_f,soln(count-1).u_f,m.parameters); + + %add convergence status into inputs + inputs=add(inputs,'converged',converged,'double'); %rift convergence criterion @@ -61,4 +70,5 @@ converged=0; end + end Index: /issm/trunk/src/mex/TriMeshProcessRifts/TriMeshProcessRifts.cpp =================================================================== --- /issm/trunk/src/mex/TriMeshProcessRifts/TriMeshProcessRifts.cpp (revision 1783) +++ /issm/trunk/src/mex/TriMeshProcessRifts/TriMeshProcessRifts.cpp (revision 1784) @@ -27,5 +27,5 @@ /*empty rifts structure: */ double* pNaN=NULL; - const char* fnames[7]; + const char* fnames[8]; mwSize ndim=2; mwSize dimensions[2] = {1,1}; @@ -239,8 +239,9 @@ fnames[5] = "fill"; fnames[6] = "friction"; + fnames[7] = "fraction"; dimensions[0]=out_numrifts; - pmxa_array=mxCreateStructArray( ndim,dimensions,7,fnames); + pmxa_array=mxCreateStructArray( ndim,dimensions,8,fnames); for (i=0;i