Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5165) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5166) @@ -744,9 +744,25 @@ /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ if (analysis_type==DiagnosticHorizAnalysisEnum){ - CreateKMatrixDiagnosticHoriz( Kgg); + int approximation; + inputs->GetParameterValue(&approximation,ApproximationEnum); + if(approximation==MacAyealApproximationEnum){ + CreateKMatrixDiagnosticMacAyeal( Kgg); + } + else if(approximation==PattynApproximationEnum){ + CreateKMatrixDiagnosticPattyn( Kgg); + } + else ISSMERROR("Approximation %s not supported yet",EnumToString(approximation)); } else if (analysis_type==AdjointHorizAnalysisEnum){ /*Same as diagnostic*/ - CreateKMatrixDiagnosticHoriz( Kgg); + int approximation; + inputs->GetParameterValue(&approximation,ApproximationEnum); + if(approximation==MacAyealApproximationEnum){ + CreateKMatrixDiagnosticMacAyeal( Kgg); + } + else if(approximation==PattynApproximationEnum){ + CreateKMatrixDiagnosticPattyn( Kgg); + } + else ISSMERROR("Approximation %s not supported yet",EnumToString(approximation)); } else if (analysis_type==DiagnosticHutterAnalysisEnum){ @@ -798,8 +814,24 @@ /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ if (analysis_type==DiagnosticHorizAnalysisEnum){ - CreatePVectorDiagnosticHoriz( pg); + int approximation; + inputs->GetParameterValue(&approximation,ApproximationEnum); + if(approximation==MacAyealApproximationEnum){ + CreatePVectorDiagnosticMacAyeal( pg); + } + else if(approximation==PattynApproximationEnum){ + CreatePVectorDiagnosticPattyn( pg); + } + else ISSMERROR("Approximation %s not supported yet",EnumToString(approximation)); } else if (analysis_type==AdjointHorizAnalysisEnum){ - CreatePVectorAdjointHoriz( pg); + int approximation; + inputs->GetParameterValue(&approximation,ApproximationEnum); + if(approximation==MacAyealApproximationEnum){ + CreatePVectorAdjointMacAyeal( pg); + } + else if(approximation==PattynApproximationEnum){ + CreatePVectorAdjointPattyn( pg); + } + else ISSMERROR("Approximation %s not supported yet",EnumToString(approximation)); } else if (analysis_type==DiagnosticHutterAnalysisEnum){ @@ -1915,5 +1947,7 @@ this->inputs->AddInput(new IntInput(ApproximationEnum,HutterApproximationEnum)); } - else ISSMERROR("Approximation type %i (%s) not supported yet",*(iomodel->elements_type+2*index+0),EnumToString(*(iomodel->elements_type+2*index+0))); + else{ + ISSMERROR("Approximation type %s not supported yet",EnumToString(*(iomodel->elements_type+2*index+0))); + } } } @@ -2041,242 +2075,4 @@ return; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixDiagnosticHoriz {{{1*/ -void Penta::CreateKMatrixDiagnosticHoriz( Mat Kgg){ - - /* local declarations */ - int i,j; - - /* node data: */ - const int numgrids=6; - const int numdof=2*numgrids; - double xyz_list[numgrids][3]; - int* doflist=NULL; - - /* 3d gaussian points: */ - int num_gauss,ig; - double* first_gauss_area_coord = NULL; - double* second_gauss_area_coord = NULL; - double* third_gauss_area_coord = NULL; - double* fourth_gauss_vert_coord = NULL; - double* area_gauss_weights = NULL; - double* vert_gauss_weights = NULL; - int ig1,ig2; - double gauss_weight1,gauss_weight2; - double gauss_coord[4]; - int order_area_gauss; - int num_vert_gauss; - int num_area_gauss; - double gauss_weight; - - /* 2d gaussian point: */ - int num_gauss2d; - double* first_gauss_area_coord2d = NULL; - double* second_gauss_area_coord2d = NULL; - double* third_gauss_area_coord2d = NULL; - double* gauss_weights2d=NULL; - double gauss_l1l2l3[3]; - - /* material data: */ - double viscosity; //viscosity - double oldviscosity; //viscosity - double newviscosity; //viscosity - - /* strain rate: */ - double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - double oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - - /* matrices: */ - double B[5][numdof]; - double Bprime[5][numdof]; - double L[2][numdof]; - double D[5][5]={0.0}; // material matrix, simple scalar matrix. - double D_scalar; - double DL[2][2]={0.0}; //for basal drag - double DL_scalar; - - /* local element matrices: */ - double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix - - double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. - double Ke_gg_drag_gaussian[numdof][numdof]; //stiffness matrix contribution from drag - double Jdet; - - /*slope: */ - double slope[2]={0.0}; - double slope_magnitude; - - /*friction: */ - double alpha2_list[3]; - double alpha2; - - double MAXSLOPE=.06; // 6 % - double MOUNTAINKEXPONENT=10; - - /*parameters: */ - double viscosity_overshoot; - - /*Collapsed formulation: */ - Tria* tria=NULL; - - /*inputs: */ - bool onwater; - bool onbed; - bool shelf; - int approximation; - Input* vx_input=NULL; - Input* vy_input=NULL; - Input* vxold_input=NULL; - Input* vyold_input=NULL; - - /*retrieve inputs :*/ - inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); - inputs->GetParameterValue(&approximation,ApproximationEnum); - - /*retrieve some parameters: */ - this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum); - - /*If on water, skip stiffness: */ - if(onwater)return; - - /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the - bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build - the stiffness matrix. */ - - if ((approximation==MacAyealApproximationEnum) && (onbed==0)){ - /*This element should be collapsed, but this element is not on the bedrock, therefore all its - * dofs have already been frozen! Do nothing: */ - return; - } - else if ((approximation==MacAyealApproximationEnum) && (onbed==1)){ - - /*This element should be collapsed into a tria element at its base. Create this tria element, - *and use its CreateKMatrix functionality to fill the global stiffness matrix: */ - - /*Depth Averaging B*/ - this->InputDepthAverageAtBase(RheologyBEnum,RheologyB2dEnum,MaterialsEnum); - - /*Call Tria function*/ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrix(Kgg); - delete tria->matice; delete tria; - - /*Delete B averaged*/ - this->matice->inputs->DeleteInput(RheologyB2dEnum); - - return; - } - else{ - - /*Implement standard penta element: */ - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist); - - /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore - get tria gaussian points as well as segment gaussian points. For tria gaussian - points, order of integration is 2, because we need to integrate the product tB*D*B' - which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian - points, same deal, which yields 3 gaussian points.*/ - - order_area_gauss=5; - num_vert_gauss=5; - - GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss); - - /*Retrieve all inputs we will be needing: */ - vx_input=inputs->GetInput(VxEnum); - vy_input=inputs->GetInput(VyEnum); - vxold_input=inputs->GetInput(VxOldEnum); - vyold_input=inputs->GetInput(VyOldEnum); - - /* Start looping on the number of gaussian points: */ - for (ig1=0; ig1GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input); - this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss_coord,vxold_input,vyold_input); - - /*Get viscosity: */ - matice->GetViscosity3d(&viscosity, &epsilon[0]); - matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]); - - /*Get B and Bprime matrices: */ - GetBPattyn(&B[0][0], &xyz_list[0][0], gauss_coord); - GetBprimePattyn(&Bprime[0][0], &xyz_list[0][0], gauss_coord); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); - - /*Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant - onto this scalar matrix, so that we win some computational time: */ - - newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); - D_scalar=newviscosity*gauss_weight*Jdet; - for (i=0;i<5;i++){ - D[i][i]=D_scalar; - } - - /* Do the triple product tB*D*Bprime: */ - TripleMultiply( &B[0][0],5,numdof,1, - &D[0][0],5,5,0, - &Bprime[0][0],5,numdof,0, - &Ke_gg_gaussian[0][0],0); - - /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */ - for( i=0; iCreateKMatrixDiagnosticHorizFriction(Kgg); - delete tria->matice; delete tria; - } - - } - - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&fourth_gauss_vert_coord); - xfree((void**)&area_gauss_weights); - xfree((void**)&vert_gauss_weights); - xfree((void**)&first_gauss_area_coord2d); - xfree((void**)&second_gauss_area_coord2d); - xfree((void**)&third_gauss_area_coord2d); - xfree((void**)&gauss_weights2d); - xfree((void**)&doflist); } /*}}}*/ @@ -2364,4 +2160,252 @@ } +/*FUNCTION Penta::CreateKMatrixDiagnosticMacAyeal{{{1*/ +void Penta::CreateKMatrixDiagnosticMacAyeal( Mat Kgg){ + + /*Collapsed formulation: */ + Tria* tria=NULL; + + /*inputs: */ + bool onwater; + bool onbed; + + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + + /*If on water, skip stiffness: */ + if(onwater)return; + + /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the + bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build + the stiffness matrix. */ + + if (onbed==0){ + /*This element should be collapsed, but this element is not on the bedrock, therefore all its + * dofs have already been frozen! Do nothing: */ + return; + } + else if (onbed==1){ + + /*This element should be collapsed into a tria element at its base. Create this tria element, + *and use its CreateKMatrix functionality to fill the global stiffness matrix: */ + + /*Depth Averaging B*/ + this->InputDepthAverageAtBase(RheologyBEnum,RheologyB2dEnum,MaterialsEnum); + + /*Call Tria function*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->CreateKMatrix(Kgg); + delete tria->matice; delete tria; + + /*Delete B averaged*/ + this->matice->inputs->DeleteInput(RheologyB2dEnum); + + return; + } +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixDiagnosticPattyn{{{1*/ +void Penta::CreateKMatrixDiagnosticPattyn( Mat Kgg){ + + /* local declarations */ + int i,j; + + /* node data: */ + const int numgrids=6; + const int numdof=2*numgrids; + double xyz_list[numgrids][3]; + int* doflist=NULL; + + /* 3d gaussian points: */ + int num_gauss,ig; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* fourth_gauss_vert_coord = NULL; + double* area_gauss_weights = NULL; + double* vert_gauss_weights = NULL; + int ig1,ig2; + double gauss_weight1,gauss_weight2; + double gauss_coord[4]; + int order_area_gauss; + int num_vert_gauss; + int num_area_gauss; + double gauss_weight; + + /* 2d gaussian point: */ + int num_gauss2d; + double* first_gauss_area_coord2d = NULL; + double* second_gauss_area_coord2d = NULL; + double* third_gauss_area_coord2d = NULL; + double* gauss_weights2d=NULL; + double gauss_l1l2l3[3]; + + /* material data: */ + double viscosity; //viscosity + double oldviscosity; //viscosity + double newviscosity; //viscosity + + /* strain rate: */ + double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + double oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + + /* matrices: */ + double B[5][numdof]; + double Bprime[5][numdof]; + double L[2][numdof]; + double D[5][5]={0.0}; // material matrix, simple scalar matrix. + double D_scalar; + double DL[2][2]={0.0}; //for basal drag + double DL_scalar; + + /* local element matrices: */ + double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix + double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. + double Ke_gg_drag_gaussian[numdof][numdof]; //stiffness matrix contribution from drag + double Jdet; + + /*slope: */ + double slope[2]={0.0}; + double slope_magnitude; + + /*friction: */ + double alpha2_list[3]; + double alpha2; + + double MAXSLOPE=.06; // 6 % + double MOUNTAINKEXPONENT=10; + + /*parameters: */ + double viscosity_overshoot; + + /*Collapsed formulation: */ + Tria* tria=NULL; + + /*inputs: */ + bool onwater; + bool onbed; + bool shelf; + Input* vx_input=NULL; + Input* vy_input=NULL; + Input* vxold_input=NULL; + Input* vyold_input=NULL; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + + /*retrieve some parameters: */ + this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum); + + /*If on water, skip stiffness: */ + if(onwater)return; + + /*Implement standard penta element: */ + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist); + + /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore + get tria gaussian points as well as segment gaussian points. For tria gaussian + points, order of integration is 2, because we need to integrate the product tB*D*B' + which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian + points, same deal, which yields 3 gaussian points.*/ + + order_area_gauss=5; + num_vert_gauss=5; + + GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss); + + /*Retrieve all inputs we will be needing: */ + vx_input=inputs->GetInput(VxEnum); + vy_input=inputs->GetInput(VyEnum); + vxold_input=inputs->GetInput(VxOldEnum); + vyold_input=inputs->GetInput(VyOldEnum); + + /* Start looping on the number of gaussian points: */ + for (ig1=0; ig1GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input); + this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss_coord,vxold_input,vyold_input); + + /*Get viscosity: */ + matice->GetViscosity3d(&viscosity, &epsilon[0]); + matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]); + + /*Get B and Bprime matrices: */ + GetBPattyn(&B[0][0], &xyz_list[0][0], gauss_coord); + GetBprimePattyn(&Bprime[0][0], &xyz_list[0][0], gauss_coord); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); + + /*Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant + onto this scalar matrix, so that we win some computational time: */ + + newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); + D_scalar=newviscosity*gauss_weight*Jdet; + for (i=0;i<5;i++){ + D[i][i]=D_scalar; + } + + /* Do the triple product tB*D*Bprime: */ + TripleMultiply( &B[0][0],5,numdof,1, + &D[0][0],5,5,0, + &Bprime[0][0],5,numdof,0, + &Ke_gg_gaussian[0][0],0); + + /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */ + for( i=0; iCreateKMatrixDiagnosticHorizFriction(Kgg); + delete tria->matice; delete tria; + } + + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&fourth_gauss_vert_coord); + xfree((void**)&area_gauss_weights); + xfree((void**)&vert_gauss_weights); + xfree((void**)&first_gauss_area_coord2d); + xfree((void**)&second_gauss_area_coord2d); + xfree((void**)&third_gauss_area_coord2d); + xfree((void**)&gauss_weights2d); + xfree((void**)&doflist); +} /*}}}*/ /*FUNCTION Penta::CreateKMatrixDiagnosticStokes {{{1*/ @@ -3108,4 +3152,70 @@ } /*}}}*/ +/*FUNCTION Penta::CreatePVectorAdjointMacAyeal{{{1*/ +void Penta::CreatePVectorAdjointMacAyeal(Vec p_g){ + + int i; + Tria* tria=NULL; + + /*inputs: */ + bool onwater; + bool onbed; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + + /*If on water, skip: */ + if(onwater) return; + + /*Bail out if this element if: + * -> Non MacAyeal and not on the surface + * -> MacAyeal(2d model) and not on bed) */ + if (!onbed){ + return; + } + else{ + + /*This element should be collapsed into a tria element at its base. Create this tria element, + * and compute pe_g*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). + tria->CreatePVectorAdjointHoriz(p_g); + delete tria->matice; delete tria; + return; + } +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorAdjointPattyn{{{1*/ +void Penta::CreatePVectorAdjointPattyn(Vec p_g){ + + int i; + Tria* tria=NULL; + + /*inputs: */ + bool onwater; + bool onsurface; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); + + /*If on water, skip: */ + if(onwater) return; + + /*Bail out if this element if: + * -> Non MacAyeal and not on the surface + * -> MacAyeal(2d model) and not on bed) */ + if (!onsurface){ + return; + } + else{ + + tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). + tria->CreatePVectorAdjointHoriz(p_g); + delete tria->matice; delete tria; + return; + } +} +/*}}}*/ /*FUNCTION Penta::CreatePVectorBalancedthickness {{{1*/ void Penta::CreatePVectorBalancedthickness( Vec pg){ @@ -3180,6 +3290,193 @@ } /*}}}*/ -/*FUNCTION Penta::CreatePVectorDiagnosticHoriz {{{1*/ -void Penta::CreatePVectorDiagnosticHoriz( Vec pg){ +/*FUNCTION Penta::CreatePVectorDiagnosticHutter{{{1*/ +void Penta::CreatePVectorDiagnosticHutter(Vec pg){ + + int i,j,k; + + const int numgrids=6; + const int NDOF2=2; + const int numdofs=NDOF2*numgrids; + int* doflist=NULL; + double pe_g[numdofs]={0.0}; + double xyz_list[numgrids][3]; + double z_list[numgrids]; + double z_segment[2]; + double slope2,constant_part; + double gauss[numdofs][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; + int node0,node1; + BeamRef* beam=NULL; + + /*gaussian points: */ + int num_gauss; + double* segment_gauss_coord=NULL; + double gauss_coord; + double* gauss_weights=NULL; + double gauss_weight; + int ig; + double slope[2]; + + double z_g; + double rho_ice,gravity,n,B; + double Jdet; + double ub,vb; + double surface,thickness; + + /*flags: */ + bool onwater; + bool onbed; + bool onsurface; + int connectivity[2]; + + /*Inputs*/ + Input* thickness_input; + Input* surface_input; + Input* slopex_input; + Input* slopey_input; + + /*recover some inputs: */ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + + /*If on water, skip: */ + if(onwater)return; + + /*recover doflist: */ + GetDofList(&doflist); + + /*recover some inputs: */ + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); + + /*recover parameters: */ + rho_ice=matpar->GetRhoIce(); + gravity=matpar->GetG(); + n=matice->GetN(); + B=matice->GetB(); + + //Initilaize beam + beam=new BeamRef(); + + //recover extra inputs + thickness_input=inputs->GetInput(ThicknessEnum); + surface_input=inputs->GetInput(SurfaceEnum); + slopex_input=inputs->GetInput(SurfaceSlopeXEnum); + slopey_input=inputs->GetInput(SurfaceSlopeYEnum); + + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + for(i=0;iGetParameterValue(&slope[0], &gauss[i][0]); + slopey_input->GetParameterValue(&slope[1], &gauss[i][0]); + surface_input->GetParameterValue(&surface, &gauss[i][0]); + thickness_input->GetParameterValue(&thickness, &gauss[i][0]); + + //compute slope2 slopex and slopey + slope2=pow(slope[0],2)+pow(slope[1],2); + + //%compute constant_part + constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2)); + + z_segment[0]=z_list[i]; + z_segment[1]=z_list[i+3]; + + node0=i; + node1=i+3; + + connectivity[0]=nodes[node0]->GetConnectivity(); + connectivity[1]=nodes[node1]->GetConnectivity(); + + /*Loop on the Gauss points: */ + for(ig=0;igGetParameterValue(&z_g, &z_segment[0],gauss_coord); + + //Get Jacobian determinant: + beam->GetJacobianDeterminant(&Jdet, &z_segment[0],gauss_coord); + + if (onsurface){ + for(j=0;jGetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + + /*If on water, skip load: */ + if(onwater)return; + + /*Figure out if this pentaelem is Macayeal. If so, then bailout, except if it is at the + bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build + the load vector. */ + + if (onbed==0){ + /*This element should be collapsed, but this element is not on the bedrock, therefore all its + * dofs have already been frozen! Do nothing: */ + return; + } + else if (onbed==1){ + + /*This element should be collapsed into a tria element at its base. Create this tria element, + *and use its CreatePVector functionality to return an elementary load vector: */ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->CreatePVector(pg); + delete tria->matice; delete tria; + return; + } +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorDiagnosticPattyn{{{1*/ +void Penta::CreatePVectorDiagnosticPattyn( Vec pg){ int i,j; @@ -3223,11 +3520,6 @@ double pe_g_gaussian[numdof]; - /*Spawning: */ - Tria* tria=NULL; - /*inputs: */ bool onwater; - bool onbed; - int approximation; Input* surface_input=NULL; Input* thickness_input=NULL; @@ -3235,89 +3527,65 @@ /*retrieve inputs :*/ inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - inputs->GetParameterValue(&approximation,ApproximationEnum); /*If on water, skip load: */ if(onwater)return; - /*Figure out if this pentaelem is Macayeal. If so, then bailout, except if it is at the - bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build - the load vector. */ - - if ((approximation==MacAyealApproximationEnum) && (onbed==0)){ - /*This element should be collapsed, but this element is not on the bedrock, therefore all its - * dofs have already been frozen! Do nothing: */ - return; - } - else if ((approximation==MacAyealApproximationEnum) && (onbed==1)){ - - /*This element should be collapsed into a tria element at its base. Create this tria element, - *and use its CreatePVector functionality to return an elementary load vector: */ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVector(pg); - delete tria->matice; delete tria; - return; - } - else{ - - /*Implement standard penta element: */ - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist); - - /*Get gaussian points and weights :*/ - order_area_gauss=2; - num_vert_gauss=3; - - GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss); - - /*Retrieve all inputs we will be needing: */ - thickness_input=inputs->GetInput(ThicknessEnum); - surface_input=inputs->GetInput(SurfaceEnum); - - /* Start looping on the number of gaussian points: */ - for (ig1=0; ig1GetParameterValue(&thickness, gauss_coord); - - /*Compute slope at gaussian point: */ - surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss_coord); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); - - /*Get nodal functions: */ - GetNodalFunctionsP1(l1l6, gauss_coord); - - /*Compute driving stress: */ - driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG(); - - /*Build pe_g_gaussian vector: */ - for (i=0;iGetInput(ThicknessEnum); + surface_input=inputs->GetInput(SurfaceEnum); + + /* Start looping on the number of gaussian points: */ + for (ig1=0; ig1GetParameterValue(&thickness, gauss_coord); + + /*Compute slope at gaussian point: */ + surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss_coord); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); + + /*Get nodal functions: */ + GetNodalFunctionsP1(l1l6, gauss_coord); + + /*Compute driving stress: */ + driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG(); + + /*Build pe_g_gaussian vector: */ + for (i=0;iGetParameterValue(&onwater,ElementOnWaterEnum); - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); - inputs->GetParameterValue(&approximation,ApproximationEnum); - - /*If on water, skip: */ - if(onwater) return; - - /*Bail out if this element if: - * -> Non MacAyeal and not on the surface - * -> MacAyeal(2d model) and not on bed) */ - if ((approximation!=MacAyealApproximationEnum && !onsurface) || (approximation==MacAyealApproximationEnum && !onbed)){ - return; - } - else if (approximation==MacAyealApproximationEnum){ - - /*This element should be collapsed into a tria element at its base. Create this tria element, - * and compute pe_g*/ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - tria->CreatePVectorAdjointHoriz(p_g); - delete tria->matice; delete tria; - return; - } - else{ - - tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). - tria->CreatePVectorAdjointHoriz(p_g); - delete tria->matice; delete tria; - return; - } -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorDiagnosticHutter{{{1*/ -void Penta::CreatePVectorDiagnosticHutter(Vec pg){ - - int i,j,k; - - const int numgrids=6; - const int NDOF2=2; - const int numdofs=NDOF2*numgrids; - int* doflist=NULL; - double pe_g[numdofs]={0.0}; - double xyz_list[numgrids][3]; - double z_list[numgrids]; - double z_segment[2]; - double slope2,constant_part; - double gauss[numdofs][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; - int node0,node1; - BeamRef* beam=NULL; - - /*gaussian points: */ - int num_gauss; - double* segment_gauss_coord=NULL; - double gauss_coord; - double* gauss_weights=NULL; - double gauss_weight; - int ig; - double slope[2]; - - double z_g; - double rho_ice,gravity,n,B; - double Jdet; - double ub,vb; - double surface,thickness; - - /*flags: */ - bool onwater; - bool onbed; - bool onsurface; - int connectivity[2]; - - /*Inputs*/ - Input* thickness_input; - Input* surface_input; - Input* slopex_input; - Input* slopey_input; - - /*recover some inputs: */ - inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - - /*If on water, skip: */ - if(onwater)return; - - /*recover doflist: */ - GetDofList(&doflist); - - /*recover some inputs: */ - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); - - /*recover parameters: */ - rho_ice=matpar->GetRhoIce(); - gravity=matpar->GetG(); - n=matice->GetN(); - B=matice->GetB(); - - //Initilaize beam - beam=new BeamRef(); - - //recover extra inputs - thickness_input=inputs->GetInput(ThicknessEnum); - surface_input=inputs->GetInput(SurfaceEnum); - slopex_input=inputs->GetInput(SurfaceSlopeXEnum); - slopey_input=inputs->GetInput(SurfaceSlopeYEnum); - - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - for(i=0;iGetParameterValue(&slope[0], &gauss[i][0]); - slopey_input->GetParameterValue(&slope[1], &gauss[i][0]); - surface_input->GetParameterValue(&surface, &gauss[i][0]); - thickness_input->GetParameterValue(&thickness, &gauss[i][0]); - - //compute slope2 slopex and slopey - slope2=pow(slope[0],2)+pow(slope[1],2); - - //%compute constant_part - constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2)); - - z_segment[0]=z_list[i]; - z_segment[1]=z_list[i+3]; - - node0=i; - node1=i+3; - - connectivity[0]=nodes[node0]->GetConnectivity(); - connectivity[1]=nodes[node1]->GetConnectivity(); - - /*Loop on the Gauss points: */ - for(ig=0;igGetParameterValue(&z_g, &z_segment[0],gauss_coord); - - //Get Jacobian determinant: - beam->GetJacobianDeterminant(&Jdet, &z_segment[0],gauss_coord); - - if (onsurface){ - for(j=0;j