Index: /issm/trunk/src/c/objects/Loads/Riftfront.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Riftfront.cpp (revision 5937) +++ /issm/trunk/src/c/objects/Loads/Riftfront.cpp (revision 5938) @@ -19,4 +19,11 @@ #include "../objects.h" /*}}}*/ + +/*Element macros*/ +#define NUMVERTICES 2 +#define NDOF1 1 +#define NDOF2 2 +#define NDOF3 3 +#define NDOF4 4 /*Riftfront constructors and destructor*/ @@ -375,4 +382,56 @@ } /*}}}*/ +/*FUNCTION Riftfront::PenaltyCreateKMatrix {{{1*/ +void Riftfront::PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs,double kmax){ + + /*Retrieve parameters: */ + ElementMatrix* Ke=NULL; + int analysis_type; + this->parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + switch(analysis_type){ + case DiagnosticHorizAnalysisEnum: + Ke=PenaltyCreateKMatrixDiagnosticHoriz(kmax); + break; + case AdjointHorizAnalysisEnum: + Ke=PenaltyCreateKMatrixDiagnosticHoriz(kmax); + break; + default: + ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); + } + + /*Add to global Vector*/ + if(Ke){ + Ke->AddToGlobal(Kgg,Kff,Kfs); + delete Ke; + } +} +/*}}}1*/ +/*FUNCTION Riftfront::PenaltyCreatePVector {{{1*/ +void Riftfront::PenaltyCreatePVector(Vec pg,Vec pf,double kmax){ + + /*Retrieve parameters: */ + ElementVector* pe=NULL; + int analysis_type; + this->parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + switch(analysis_type){ + case DiagnosticHorizAnalysisEnum: + pe=PenaltyCreatePVectorDiagnosticHoriz(kmax); + break; + case AdjointHorizAnalysisEnum: + /*No penalty applied on load vector*/ + break; + default: + ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); + } + + /*Add to global Vector*/ + if(pe){ + pe->AddToGlobal(pg,pf); + delete pe; + } +} +/*}}}1*/ /*FUNCTION Riftfront::CreateKMatrix {{{1*/ void Riftfront::CreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs){ @@ -384,16 +443,22 @@ void Riftfront::CreatePVector(Vec pg,Vec pf){ /*do nothing: */ -} -/*}}}1*/ -/*FUNCTION Riftfront::PenaltyCreateKMatrix {{{1*/ -void Riftfront::PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs,double kmax){ - - int i; - int j; - const int numgrids = MAX_RIFTFRONT_GRIDS; + return; +} +/*}}}1*/ +/*FUNCTION Riftfront::InAnalysis{{{1*/ +bool Riftfront::InAnalysis(int in_analysis_type){ + if (in_analysis_type==this->analysis_type) return true; + else return false; +} +/*}}}*/ + +/*Riftfront numerics*/ +/*FUNCTION Riftfront::PenaltyCreateKMatrixDiagnosticHoriz {{{1*/ +ElementMatrix* Riftfront::PenaltyCreateKMatrixDiagnosticHoriz(double kmax){ + + const int numdof = NDOF2*NUMVERTICES; + int i,j; int dofs[1] = {0}; double Ke_gg[4][4]; - const int numdof = 2 *numgrids; - int* doflist=NULL; double thickness; double h[2]; @@ -413,99 +478,74 @@ /*enum of element? */ if(elements[0]->Enum()!=TriaEnum)ISSMERROR(" only Tria element allowed for Riftfront load!"); - - /*recover elements on both side of rift: */ tria1=(Tria*)elements[0]; tria2=(Tria*)elements[1]; - - /* Get node coordinates and dof list: */ - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - - /* Set Ke_gg to 0: */ - for(i=0;iactive) return NULL; + ElementMatrix* Ke=NewElementMatrix(nodes,NUMVERTICES,this->parameters); /*Get some parameters: */ this->parameters->FindParam(&penalty_offset,PenaltyOffsetEnum); this->inputs->GetParameterValue(&friction,FrictionEnum); - - if(this->active){ - - /*There is contact, we need to constrain the normal velocities (zero penetration), and the - *contact slip friction. */ - - /*Recover thickness: */ - tria1->GetParameterValue(&h[0],nodes[0],ThicknessEnum); - tria2->GetParameterValue(&h[1],nodes[1],ThicknessEnum); - - if (h[0]!=h[1])ISSMERROR(" different thicknesses not supported for rift fronts"); - thickness=h[0]; - - /*From Peter Wriggers book (Computational Contact Mechanics, p191): */ - //First line: - Ke_gg[0][0]+=pow(normal[0],2)*kmax*pow(10,penalty_offset); - Ke_gg[0][1]+=normal[0]*normal[1]*kmax*pow(10,penalty_offset); - Ke_gg[0][2]+=-pow(normal[0],2)*kmax*pow(10,penalty_offset); - Ke_gg[0][3]+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset); - //Second line: - Ke_gg[1][0]+=normal[0]*normal[1]*kmax*pow(10,penalty_offset); - Ke_gg[1][1]+=pow(normal[1],2)*kmax*pow(10,penalty_offset); - Ke_gg[1][2]+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset); - Ke_gg[1][3]+=-pow(normal[1],2)*kmax*pow(10,penalty_offset); - //Third line: - Ke_gg[2][0]+=-pow(normal[0],2)*kmax*pow(10,penalty_offset); - Ke_gg[2][1]+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset); - Ke_gg[2][2]+=pow(normal[0],2)*kmax*pow(10,penalty_offset); - Ke_gg[2][3]+=normal[0]*normal[1]*kmax*pow(10,penalty_offset); - //Fourth line: - Ke_gg[3][0]+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset); - Ke_gg[3][1]+=-pow(normal[1],2)*kmax*pow(10,penalty_offset); - Ke_gg[3][2]+=normal[0]*normal[1]*kmax*pow(10,penalty_offset); - Ke_gg[3][3]+=pow(normal[1],2)*kmax*pow(10,penalty_offset); - - /*Now take care of the friction: of type sigma=frictiontangent_velocity2-tangent_velocity1)*/ - - //First line: - Ke_gg[0][0]+=pow(normal[1],2)*thickness*length*friction; - Ke_gg[0][1]+=-normal[0]*normal[1]*thickness*length*friction; - Ke_gg[0][2]+=-pow(normal[1],2)*thickness*length*friction; - Ke_gg[0][3]+=normal[0]*normal[1]*thickness*length*friction; - //Second line: - Ke_gg[1][0]+=-normal[0]*normal[1]*thickness*length*friction; - Ke_gg[1][1]+=pow(normal[0],2)*thickness*length*friction; - Ke_gg[1][2]+=normal[0]*normal[1]*thickness*length*friction; - Ke_gg[1][3]+=-pow(normal[0],2)*thickness*length*friction; - //Third line: - Ke_gg[2][0]+=-pow(normal[1],2)*thickness*length*friction; - Ke_gg[2][1]+=normal[0]*normal[1]*thickness*length*friction; - Ke_gg[2][2]+=pow(normal[1],2)*thickness*length*friction; - Ke_gg[2][3]+=-normal[0]*normal[1]*thickness*length*friction; - //Fourth line: - Ke_gg[3][0]+=normal[0]*normal[1]*thickness*length*friction; - Ke_gg[3][1]+=-pow(normal[0],2)*thickness*length*friction; - Ke_gg[3][2]+=-normal[0]*normal[1]*thickness*length*friction; - Ke_gg[3][3]+=pow(normal[0],2)*thickness*length*friction; - - /*Add Ke_gg to global matrix Kgg: */ - MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES); - } - else{ - /*the grids on both sides of the rift do not penetrate. PenaltyCreatePVector will - *take care of adding point loads to simulate pressure on the rift flanks. But as far as stiffness, - there is none (0 stiffness implies decoupling of the flank rifts, which is exactly what we want): */ - } - /*Free ressources:*/ - xfree((void**)&doflist); - -} -/*}}}1*/ -/*FUNCTION Riftfront::PenaltyCreatePVector {{{1*/ -void Riftfront::PenaltyCreatePVector(Vec pg,Vec pf,double kmax){ - - int i ,j; - const int numgrids = MAX_RIFTFRONT_GRIDS; - double pe_g[4]={0.0}; - const int numdof = 2 *numgrids; - int* doflist=NULL; - + tria1->GetParameterValue(&h[0],nodes[0],ThicknessEnum); + tria2->GetParameterValue(&h[1],nodes[1],ThicknessEnum); + if (h[0]!=h[1])ISSMERROR(" different thicknesses not supported for rift fronts"); + thickness=h[0]; + + /*There is contact, we need to constrain the normal velocities (zero penetration), and the + *contact slip friction. */ + + /*From Peter Wriggers book (Computational Contact Mechanics, p191): */ + Ke->values[0*numdof+0]+= +pow(normal[0],2)*kmax*pow(10,penalty_offset); + Ke->values[0*numdof+1]+= +normal[0]*normal[1]*kmax*pow(10,penalty_offset); + Ke->values[0*numdof+2]+= -pow(normal[0],2)*kmax*pow(10,penalty_offset); + Ke->values[0*numdof+3]+= -normal[0]*normal[1]*kmax*pow(10,penalty_offset); + + Ke->values[1*numdof+0]+= +normal[0]*normal[1]*kmax*pow(10,penalty_offset); + Ke->values[1*numdof+1]+= +pow(normal[1],2)*kmax*pow(10,penalty_offset); + Ke->values[1*numdof+2]+= -normal[0]*normal[1]*kmax*pow(10,penalty_offset); + Ke->values[1*numdof+3]+= -pow(normal[1],2)*kmax*pow(10,penalty_offset); + + Ke->values[2*numdof+0]+= -pow(normal[0],2)*kmax*pow(10,penalty_offset); + Ke->values[2*numdof+1]+= -normal[0]*normal[1]*kmax*pow(10,penalty_offset); + Ke->values[2*numdof+2]+= +pow(normal[0],2)*kmax*pow(10,penalty_offset); + Ke->values[2*numdof+3]+= +normal[0]*normal[1]*kmax*pow(10,penalty_offset); + + Ke->values[3*numdof+0]+= -normal[0]*normal[1]*kmax*pow(10,penalty_offset); + Ke->values[3*numdof+1]+= -pow(normal[1],2)*kmax*pow(10,penalty_offset); + Ke->values[3*numdof+2]+= +normal[0]*normal[1]*kmax*pow(10,penalty_offset); + Ke->values[3*numdof+3]+= +pow(normal[1],2)*kmax*pow(10,penalty_offset); + + /*Now take care of the friction: of type sigma=frictiontangent_velocity2-tangent_velocity1)*/ + + Ke->values[0*numdof+0]+= +pow(normal[1],2)*thickness*length*friction; + Ke->values[0*numdof+1]+= -normal[0]*normal[1]*thickness*length*friction; + Ke->values[0*numdof+2]+= -pow(normal[1],2)*thickness*length*friction; + Ke->values[0*numdof+3]+= +normal[0]*normal[1]*thickness*length*friction; + + Ke->values[1*numdof+0]+= -normal[0]*normal[1]*thickness*length*friction; + Ke->values[1*numdof+1]+= +pow(normal[0],2)*thickness*length*friction; + Ke->values[1*numdof+2]+= +normal[0]*normal[1]*thickness*length*friction; + Ke->values[1*numdof+3]+= -pow(normal[0],2)*thickness*length*friction; + + Ke->values[2*numdof+0]+= -pow(normal[1],2)*thickness*length*friction; + Ke->values[2*numdof+1]+= +normal[0]*normal[1]*thickness*length*friction; + Ke->values[2*numdof+2]+= +pow(normal[1],2)*thickness*length*friction; + Ke->values[2*numdof+3]+= -normal[0]*normal[1]*thickness*length*friction; + + Ke->values[3*numdof+0]+= +normal[0]*normal[1]*thickness*length*friction; + Ke->values[3*numdof+1]+= -pow(normal[0],2)*thickness*length*friction; + Ke->values[3*numdof+2]+= -normal[0]*normal[1]*thickness*length*friction; + Ke->values[3*numdof+3]+= +pow(normal[0],2)*thickness*length*friction; + + /*Clean up and return*/ + return Ke; +} +/*}}}1*/ +/*FUNCTION Riftfront::PenaltyCreatePVectorDiagnosticHoriz {{{1*/ +ElementVector* Riftfront::PenaltyCreatePVectorDiagnosticHoriz(double kmax){ + + const int numdof = NDOF2*NUMVERTICES; + int i,j; double rho_ice; double rho_water; @@ -523,5 +563,4 @@ bool shelf; - /*Objects: */ Element **elements = NULL; @@ -531,5 +570,4 @@ Matpar *matpar = NULL; - /*Recover hook objects: */ elements=(Element**)helements->deliverp(); @@ -539,102 +577,78 @@ /*enum of element? */ if(elements[0]->Enum()!=TriaEnum)ISSMERROR(" only Tria element allowed for Riftfront load!"); - - /*recover elements on both side of rift: */ tria1=(Tria*)elements[0]; tria2=(Tria*)elements[1]; - /* Get node coordinates and dof list: */ - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + /*Initialize Element Matrix*/ + if(this->active) return NULL; /*The penalty is active. No loads implied here.*/ + ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters); /*Get some inputs: */ this->inputs->GetParameterValue(&fill,FillEnum); this->inputs->GetParameterValue(&shelf,SegmentOnIceShelfEnum); - - if(!this->active){ - /*Ok, this rift is opening. We should put loads on both sides of the rift flanks. Because we are dealing with contact mechanics, - * and we want to avoid zigzagging of the loads, we want lump the loads onto grids, not onto surfaces between grids.:*/ - - /*Ok, to compute the pressure, we are going to need material properties, thickness and bed for the two grids. We assume those properties to - * be the same across the rift.: */ - - rho_ice=matpar->GetRhoIce(); - rho_water=matpar->GetRhoWater(); - gravity=matpar->GetG(); - - /*get thickness: */ - tria1->GetParameterValue(&h[0],nodes[0],ThicknessEnum); - tria2->GetParameterValue(&h[1],nodes[1],ThicknessEnum); - - if (h[0]!=h[1])ISSMERROR(" different thicknesses not supported for rift fronts"); - thickness=h[0]; - - tria1->GetParameterValue(&b[0],nodes[0],BedEnum); - tria2->GetParameterValue(&b[1],nodes[1],BedEnum); - - if (b[0]!=b[1])ISSMERROR(" different beds not supported for rift fronts"); - bed=b[0]; - - /*Ok, now compute the pressure (in norm) that is being applied to the flanks, depending on the type of fill: */ - if(fill==WaterEnum){ - if(shelf){ - /*We are on an ice shelf, hydrostatic equilibrium is used to determine the pressure for water fill: */ - pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2 - rho_water*gravity*pow(bed,(double)2)/(double)2; - } - else{ - //We are on an icesheet, we assume the water column fills the entire front: */ - pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2 - rho_water*gravity*pow(thickness,(double)2)/(double)2; - } - } - 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==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) ISSMERROR("%s%i%s","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; + rho_ice=matpar->GetRhoIce(); + rho_water=matpar->GetRhoWater(); + gravity=matpar->GetG(); + tria1->GetParameterValue(&h[0],nodes[0],ThicknessEnum); + tria2->GetParameterValue(&h[1],nodes[1],ThicknessEnum); + if (h[0]!=h[1])ISSMERROR(" different thicknesses not supported for rift fronts"); + thickness=h[0]; + tria1->GetParameterValue(&b[0],nodes[0],BedEnum); + tria2->GetParameterValue(&b[1],nodes[1],BedEnum); + if (b[0]!=b[1])ISSMERROR(" different beds not supported for rift fronts"); + bed=b[0]; + + /*Ok, this rift is opening. We should put loads on both sides of the rift flanks. Because we are dealing with contact mechanics, + * and we want to avoid zigzagging of the loads, we want lump the loads onto grids, not onto surfaces between grids.:*/ + + /*Ok, to compute the pressure, we are going to need material properties, thickness and bed for the two grids. We assume those properties to + * be the same across the rift.: */ + + /*Ok, now compute the pressure (in norm) that is being applied to the flanks, depending on the type of fill: */ + if(fill==WaterEnum){ + if(shelf){ + /*We are on an ice shelf, hydrostatic equilibrium is used to determine the pressure for water fill: */ + pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2 - rho_water*gravity*pow(bed,(double)2)/(double)2; } else{ - ISSMERROR("%s%i%s","fill type ",fill," not supported yet."); + //We are on an icesheet, we assume the water column fills the entire front: */ + pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2 - rho_water*gravity*pow(thickness,(double)2)/(double)2; } - - /*Ok, add contribution to first grid, along the normal i==0: */ - for (j=0;j<2;j++){ - pe_g[j]+=pressure*normal[j]*length; - } - - /*Add contribution to second grid, along the opposite normal: i==1 */ - for (j=0;j<2;j++){ - pe_g[2+j]+= -pressure*normal[j]*length; - } - /*Add pe_g to global vector pg; */ - VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); - + } + 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==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) ISSMERROR("%s%i%s","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{ - /*The penalty is active. No loads implied here.*/ - } - - /*Free ressources:*/ - xfree((void**)&doflist); - -} -/*}}}1*/ -/*FUNCTION Riftfront::InAnalysis{{{1*/ -bool Riftfront::InAnalysis(int in_analysis_type){ - if (in_analysis_type==this->analysis_type) return true; - else return false; -} -/*}}}*/ - -/*Riftfront numerics*/ + ISSMERROR("%s%i%s","fill type ",fill," not supported yet."); + } + + /*Ok, add contribution to first grid, along the normal i==0: */ + for (j=0;j<2;j++){ + pe->values[j]+=pressure*normal[j]*length; + } + + /*Add contribution to second grid, along the opposite normal: i==1 */ + for (j=0;j<2;j++){ + pe->values[2+j]+= -pressure*normal[j]*length; + } + + /*Clean up and return*/ + return pe; +} +/*}}}1*/ /*FUNCTION Riftfront::Constrain {{{1*/ #define _ZIGZAGCOUNTER_ @@ -731,42 +745,4 @@ } /*}}}1*/ -/*FUNCTION Riftfront::GetDofList {{{1*/ -void Riftfront::GetDofList(int** pdoflist,int approximation,int setenum){ - - int i,j; - int numberofdofs=0; - int count=0; - - /*output: */ - int* doflist=NULL; - - /*pointers: */ - Node** nodes=NULL; - - /*recover pointers: */ - nodes=(Node**)hnodes->deliverp(); - - /*Some checks for debugging*/ - ISSMASSERT(nodes); - - /*Figure out size of doflist: */ - for(i=0;iGetNumberOfDofs(approximation,setenum); - } - - /*Allocate: */ - doflist=(int*)xmalloc(numberofdofs*sizeof(int)); - - /*Populate: */ - count=0; - for(i=0;iGetDofList(doflist+count,approximation,setenum); - count+=nodes[i]->GetNumberOfDofs(approximation,setenum); - } - - /*Assign output pointers:*/ - *pdoflist=doflist; -} -/*}}}*/ /*FUNCTION Riftfront::IsFrozen{{{1*/ bool Riftfront::IsFrozen(void){ Index: /issm/trunk/src/c/objects/Loads/Riftfront.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Riftfront.h (revision 5937) +++ /issm/trunk/src/c/objects/Loads/Riftfront.h (revision 5938) @@ -13,7 +13,4 @@ class Inputs; class IoModel; - -#define MAX_RIFTFRONT_GRIDS 2 //max number of grids on a rift flank, only 2 because 2d for now. -#define RIFTFRONTSTRING 20 //max string length /*}}}*/ @@ -81,6 +78,7 @@ /*}}}*/ /*Riftfront specific routines: {{{1*/ - void GetDofList(int** pdoflist,int approximation_enum,int setenum); bool PreStable(); + ElementMatrix* PenaltyCreateKMatrixDiagnosticHoriz(double kmax); + ElementVector* PenaltyCreatePVectorDiagnosticHoriz(double kmax); void SetPreStable(); int PreConstrain(int* punstable);