Index: /issm/trunk/src/c/EnumDefinitions/EnumAsString.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumAsString.cpp (revision 4884) +++ /issm/trunk/src/c/EnumDefinitions/EnumAsString.cpp (revision 4885) @@ -225,4 +225,5 @@ case P1Enum : return "P1"; case P1DGEnum : return "P1DG"; + case MiniEnum : return "Mini"; case DoubleElementResultEnum : return "DoubleElementResult"; case TriaVertexElementResultEnum : return "TriaVertexElementResult"; Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 4884) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 4885) @@ -254,4 +254,5 @@ P1Enum, P1DGEnum, + MiniEnum, /*}}}*/ /*Results{{{1*/ Index: /issm/trunk/src/c/EnumDefinitions/StringAsEnum.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/StringAsEnum.cpp (revision 4884) +++ /issm/trunk/src/c/EnumDefinitions/StringAsEnum.cpp (revision 4885) @@ -223,4 +223,5 @@ else if (strcmp(name,"P1")==0) return P1Enum; else if (strcmp(name,"P1DG")==0) return P1DGEnum; + else if (strcmp(name,"Mini")==0) return MiniEnum; else if (strcmp(name,"DoubleElementResult")==0) return DoubleElementResultEnum; else if (strcmp(name,"TriaVertexElementResult")==0) return TriaVertexElementResultEnum; Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 4884) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 4885) @@ -88,4 +88,8 @@ penta=new Penta(); + + //deal with PentaRef mother class + penta->element_type_list=(int*)xmalloc(this->numanalyses*sizeof(int)); + for(i=0;inumanalyses;i++) penta->element_type_list[i]=this->element_type_list[i]; //deal with PentaHook mother class @@ -161,7 +165,10 @@ /*this time, no need to get enum type, the pointer directly points to the beginning of the *object data (thanks to DataSet::Demarshall):*/ - memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id); memcpy(&numanalyses,marshalled_dataset,sizeof(numanalyses));marshalled_dataset+=sizeof(numanalyses); + + /*demarshall Ref: */ + this->element_type_list=(int*)xmalloc(this->numanalyses*sizeof(int)); + for(i=0;iMarshallSize() +hmatpar->MarshallSize() @@ -633,4 +646,7 @@ parametersin->FindParam(&analysis_counter,AnalysisCounterEnum); + /*Get Element type*/ + this->element_type=this->element_type_list[analysis_counter]; + /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective * datasets, using internal ids and offsets hidden in hooks: */ @@ -641,23 +657,4 @@ /*Now, go pick up the objects inside the hooks: */ - if (this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); - else this->nodes=NULL; - this->matice=(Matice*)this->hmatice->delivers(); - this->matpar=(Matpar*)this->hmatpar->delivers(); - this->neighbors=(Penta**)this->hneighbors->deliverp(); - - /*point parameters to real dataset: */ - this->parameters=parametersin; -} -/*}}}*/ -/*FUNCTION Penta::SetCurrentConfiguration {{{1*/ -void Penta::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){ - - int analysis_counter; - - /*go into parameters and get the analysis_counter: */ - parametersin->FindParam(&analysis_counter,AnalysisCounterEnum); - - /*Pick up the objects inside the hooks: */ if (this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); else this->nodes=NULL; @@ -958,54 +955,4 @@ } /*}}}*/ -/*FUNCTION Penta::GradjDrag {{{1*/ -void Penta::GradjDrag(Vec gradient){ - - int i; - Tria* tria=NULL; - TriaVertexInput* triavertexinput=NULL; - double temp_gradient[6]={0,0,0,0,0,0}; - - /*inputs: */ - bool onwater; - bool onbed; - bool shelf; - int analysis_type; - - /*retrieve parameters: */ - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*retrieve inputs :*/ - inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); - - /*If on water, skip: */ - if(onwater)return; - - /*If on shelf, skip: */ - if(shelf)return; - - /*Bail out if this element does not touch the bedrock: */ - if (!onbed) return; - - if (analysis_type==AdjointHorizAnalysisEnum){ - - /*MacAyeal or Pattyn*/ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->GradjDrag(gradient); - delete tria; - } - else if (analysis_type==AdjointStokesAnalysisEnum){ - - /*Stokes*/ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->GradjDragStokes(gradient); - delete tria; - } - else ISSMERROR("analysis %s not supported yet",EnumAsString(analysis_type)); - - -} -/*}}}*/ /*FUNCTION Penta::InputAXPY{{{1*/ void Penta::InputAXPY(int YEnum, double scalar, int XEnum){ @@ -1692,4 +1639,26 @@ } /*}}}*/ +/*FUNCTION Penta::SetCurrentConfiguration {{{1*/ +void Penta::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){ + + int analysis_counter; + + /*go into parameters and get the analysis_counter: */ + parametersin->FindParam(&analysis_counter,AnalysisCounterEnum); + + /*Get Element type*/ + this->element_type=this->element_type_list[analysis_counter]; + + /*Pick up the objects inside the hooks: */ + if (this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); + else this->nodes=NULL; + this->matice=(Matice*)this->hmatice->delivers(); + this->matpar=(Matpar*)this->hmatpar->delivers(); + this->neighbors=(Penta**)this->hneighbors->deliverp(); + + /*point parameters to real dataset: */ + this->parameters=parametersin; +} +/*}}}*/ /*FUNCTION Penta::SurfaceArea {{{1*/ double Penta::SurfaceArea(void){ @@ -1756,4 +1725,8 @@ /*P1 Discontinuous Galerkin*/ penta_type=P1DGEnum; + } + else if (analysis_type==StokesAnalysisEnum){ + /*MINI Continuous Galerkin*/ + penta_type=MiniEnum; } else{ @@ -1968,99 +1941,1998 @@ /*Penta specific routines: */ -/*FUNCTION Penta::SpawnTria {{{1*/ -Tria* Penta::SpawnTria(int g0, int g1, int g2){ +/*FUNCTION Penta::CreateKMatrixBalancedthickness {{{1*/ + +void Penta::CreateKMatrixBalancedthickness(Mat Kgg){ + + /*Collapsed formulation: */ + Tria* tria=NULL; + + /*flags: */ + bool onwater; + bool onbed; + + /*recover some inputs: */ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + + /*If on water, skip: */ + if(onwater)return; + + /*Is this element on the bed? :*/ + if(!onbed)return; + + /*Spawn Tria element from the base of the Penta: */ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->CreateKMatrix(Kgg); + delete tria; + return; + +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixBalancedvelocities {{{1*/ + +void Penta::CreateKMatrixBalancedvelocities(Mat Kgg){ + + /*Collapsed formulation: */ + Tria* tria=NULL; + + /*flags: */ + bool onbed; + bool onwater; + + /*recover some inputs: */ + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + + /*If on water, skip: */ + if(onwater)return; + + /*Is this element on the bed? :*/ + if(!onbed)return; + + /*Spawn Tria element from the base of the Penta: */ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->CreateKMatrix(Kgg); + delete tria; + 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[numdof]; + int numberofdofspernode; + + /* 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 collapse; + 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(&collapse,CollapseEnum); + 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; + + /*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 ((collapse==1) && (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 ((collapse==1) && (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: */ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->CreateKMatrix(Kgg); + delete tria; + return; + } + else{ + + /*Implement standard penta element: */ + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /*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; + } + + } + + 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); + +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixDiagnosticHutter{{{1*/ +void Penta::CreateKMatrixDiagnosticHutter(Mat Kgg){ + + /*Collapsed formulation: */ + Beam* beam=NULL; + int i; + + /*flags: */ + bool onwater; + + /*recover some inputs: */ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + + /*If on water, skip: */ + if(onwater)return; + + /*Spawn 3 beam elements: */ + for(i=0;i<3;i++){ + beam=(Beam*)SpawnBeam(i,i+3); //[0 3], [1 4] and [2 5] are the four vertical edges of the Penta + beam->CreateKMatrix(Kgg); + } + + /*clean up*/ + delete beam; + +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixDiagnosticStokes {{{1*/ +void Penta::CreateKMatrixDiagnosticStokes( Mat Kgg){ + + int i,j; + + const int numgrids=6; + const int DOFPERGRID=4; + const int numdof=numgrids*DOFPERGRID; + int doflist[numdof]; + int numberofdofspernode; + + const int numgrids2d=3; + const int numdof2d=numgrids2d*DOFPERGRID; + + int dofs[3]={0,1,2}; + + double K_terms[numdof][numdof]={0.0}; + + /*Material properties: */ + double gravity,rho_ice,rho_water; + + /*Collapsed formulation: */ + Tria* tria=NULL; + + /*Grid data: */ + double xyz_list[numgrids][3]; + + /*parameters: */ + double xyz_list_tria[numgrids2d][3]; + double surface_normal[3]; + double bed_normal[3]; + double thickness; + + /*matrices: */ + double Ke_temp[27][27]={0.0}; //for the six nodes and the bubble + double Ke_reduced[numdof][numdof]; //for the six nodes only + double Ke_gaussian[27][27]; + double Ke_drag_gaussian[numdof2d][numdof2d]; + double B[8][27]; + double B_prime[8][27]; + double LStokes[14][numdof2d]; + double LprimeStokes[14][numdof2d]; + double Jdet; + double Jdet2d; + double D[8][8]={0.0}; + double D_scalar; + double tBD[27][8]; + double DLStokes[14][14]={0.0}; + double tLDStokes[numdof2d][14]; + + /* gaussian points: */ + int num_area_gauss; + int igarea,igvert; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* vert_gauss_coord = NULL; + double* area_gauss_weights = NULL; + double* vert_gauss_weights = NULL; + double gaussgrids[numgrids][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}}; + + /* specific gaussian point: */ + double gauss_weight,area_gauss_weight,vert_gauss_weight; + double gauss_coord[4]; + double gauss_coord_tria[3]; + int area_order=5; + int num_vert_gauss=5; + + double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + double viscosity; + double alpha2_gauss; + Friction* friction=NULL; + + /*parameters: */ + double stokesreconditioning; + int analysis_type; + + /*inputs: */ + bool onwater; + bool onbed; + bool shelf; + bool isstokes; + + /*inputs: */ + Input* vx_input=NULL; + Input* vy_input=NULL; + Input* vz_input=NULL; + + /*retrive parameters: */ + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + inputs->GetParameterValue(&isstokes,IsStokesEnum); + + /*If on water or not Stokes, skip stiffness: */ + if(onwater || !isstokes) return; + + /*recovre material parameters: */ + rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetRhoIce(); + gravity=matpar->GetG(); + + /*retrieve some parameters: */ + this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /* 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.*/ + + area_order=5; + num_vert_gauss=5; + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss); + + /*Retrieve all inputs we will be needing: */ + vx_input=inputs->GetInput(VxEnum); + vy_input=inputs->GetInput(VyEnum); + vz_input=inputs->GetInput(VzEnum); + + /* Start looping on the number of gaussian points: */ + for (igarea=0; igareaGetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input); + + /*Get viscosity: */ + matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + + /*Get B and Bprime matrices: */ + GetBStokes(&B[0][0],&xyz_list[0][0],gauss_coord); + GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0], gauss_coord); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],&gauss_coord[0]); + + /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant + * onto this scalar matrix, so that we win some computational time: */ + D_scalar=gauss_weight*Jdet; + for (i=0;i<6;i++){ + D[i][i]=D_scalar*viscosity; + } + for (i=6;i<8;i++){ + D[i][i]=-D_scalar*stokesreconditioning; + } + + /* Do the triple product tB*D*Bprime: */ + MatrixMultiply(&B[0][0],8,27,1,&D[0][0],8,8,0,&tBD[0][0],0); + MatrixMultiply(&tBD[0][0],27,8,0,&B_prime[0][0],8,27,0,&Ke_gaussian[0][0],0); + + /*Add Ke_gaussian and Ke_gaussian to terms in pKe. Watch out for column orientation from matlab: */ + for(i=0;i<27;i++){ + for(j=0;j<27;j++){ + Ke_temp[i][j]+=Ke_gaussian[i][j]; + } + } + } + } + + if((onbed==1) && (shelf==0)){ + + /*build friction object, used later on: */ + friction=new Friction("3d",inputs,matpar,analysis_type); + + for(i=0;iGetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria); + + /*Get L matrix if viscous basal drag present: */ + GetLStokes(&LStokes[0][0], gauss_coord_tria); + GetLprimeStokes(&LprimeStokes[0][0], &xyz_list[0][0], gauss_coord_tria, gauss_coord); + + /*Compute strain rate: */ + this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input); + + /*Get viscosity at last iteration: */ + matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + + /*Get normal vecyor to the bed */ + SurfaceNormal(&surface_normal[0],xyz_list_tria); + + bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result + bed_normal[1]=-surface_normal[1]; + bed_normal[2]=-surface_normal[2]; + + /*Calculate DL on gauss point */ + friction->GetAlpha2(&alpha2_gauss, gauss_coord,VxEnum,VyEnum,VzEnum); + + DLStokes[0][0]=alpha2_gauss*gauss_weight*Jdet2d; + DLStokes[1][1]=alpha2_gauss*gauss_weight*Jdet2d; + DLStokes[2][2]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[0]*bed_normal[2]; + DLStokes[3][3]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[1]*bed_normal[2]; + DLStokes[4][4]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[0]*bed_normal[2]; + DLStokes[5][5]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[1]*bed_normal[2]; + DLStokes[6][6]=-viscosity*gauss_weight*Jdet2d*bed_normal[0]; + DLStokes[7][7]=-viscosity*gauss_weight*Jdet2d*bed_normal[1]; + DLStokes[8][8]=-viscosity*gauss_weight*Jdet2d*bed_normal[2]; + DLStokes[9][8]=-viscosity*gauss_weight*Jdet2d*bed_normal[0]/2.0; + DLStokes[10][10]=-viscosity*gauss_weight*Jdet2d*bed_normal[1]/2.0; + DLStokes[11][11]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[0]; + DLStokes[12][12]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[1]; + DLStokes[13][13]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[2]; + + /* Do the triple product tL*D*L: */ + MatrixMultiply(&LStokes[0][0],14,numdof2d,1,&DLStokes[0][0],14,14,0,&tLDStokes[0][0],0); + MatrixMultiply(&tLDStokes[0][0],numdof2d,14,0,&LprimeStokes[0][0],14,numdof2d,0,&Ke_drag_gaussian[0][0],0); + + for(i=0;iGetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); + + /*If on water, skip stiffness: */ + if(onwater)return; + + /*If this element is on the surface, we have a dynamic boundary condition that applies, as a stiffness + * matrix: */ + if(onsurface){ + tria=(Tria*)SpawnTria(3,4,5); //nodes 3,4 and 5 are on the surface + tria->CreateKMatrixDiagnosticSurfaceVert(Kgg); + delete tria; + } + + /*Now, onto the formulation for the vertical velocity: */ + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /*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=2; + num_vert_gauss=2; + + 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); + + /* Start looping on the number of gaussian points: */ + for (ig1=0; ig1GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + + /*If on water, skip: */ + if(onwater)return; + + if (!onbed){ + return; + } + else{ + + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->CreateKMatrixMelting(Kgg); + delete tria; + return; + } +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixPrognostic {{{1*/ + +void Penta::CreateKMatrixPrognostic(Mat Kgg){ + + /*Collapsed formulation: */ + 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; + + /*Is this element on the bed? :*/ + if(!onbed)return; + + /*Spawn Tria element from the base of the Penta: */ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->CreateKMatrix(Kgg); + delete tria; + return; + +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixSlope {{{1*/ + +void Penta::CreateKMatrixSlope(Mat Kgg){ + + /*Collapsed formulation: */ + 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; + + /*Is this element on the bed? :*/ + if(!onbed)return; + + /*Spawn Tria element from the base of the Penta: */ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->CreateKMatrix(Kgg); + delete tria; + return; + +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixThermal {{{1*/ +void Penta::CreateKMatrixThermal(Mat Kgg){ + + /* local declarations */ + int i,j; + int found=0; + + /* node data: */ + const int numgrids=6; + const int NDOF1=1; + const int numdof=NDOF1*numgrids; + double xyz_list[numgrids][3]; + int doflist[numdof]; + int numberofdofspernode; + + /* gaussian points: */ + int num_area_gauss,igarea,igvert; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* vert_gauss_coord = NULL; + double* area_gauss_weights = NULL; + double* vert_gauss_weights = NULL; + double gauss_weight,area_gauss_weight,vert_gauss_weight; + double gauss_coord[4]; + double gauss_l1l2l3[3]; + + int area_order=5; + int num_vert_gauss=5; + + int dofs[3]={0,1,2}; + double K[2][2]={0.0}; + + double u,v,w; + + /*matrices: */ + double K_terms[numdof][numdof]={0.0}; + double Ke_gaussian_conduct[numdof][numdof]; + double Ke_gaussian_advec[numdof][numdof]; + double Ke_gaussian_artdiff[numdof][numdof]; + double Ke_gaussian_transient[numdof][numdof]; + double B[3][numdof]; + double Bprime[3][numdof]; + double B_conduct[3][numdof]; + double B_advec[3][numdof]; + double B_artdiff[2][numdof]; + double Bprime_advec[3][numdof]; + double L[numdof]; + double D_scalar; + double D[3][3]; + double l1l2l3[3]; + double tl1l2l3D[3]; + double tBD[3][numdof]; + double tBD_conduct[3][numdof]; + double tBD_advec[3][numdof]; + double tBD_artdiff[3][numdof]; + double tLD[numdof]; + + double Jdet; + + /*Material properties: */ + double gravity,rho_ice,rho_water; + double heatcapacity,thermalconductivity; + double mixed_layer_capacity,thermal_exchange_velocity; + + /*parameters: */ + double dt,epsvel; + bool artdiff; + + /*Collapsed formulation: */ + Tria* tria=NULL; + + + /*inputs: */ + bool onwater; + bool onbed; + bool shelf; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + + /*If on water, skip: */ + if(onwater)return; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + // /*recovre material parameters: */ + rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetRhoIce(); + gravity=matpar->GetG(); + heatcapacity=matpar->GetHeatCapacity(); + thermalconductivity=matpar->GetThermalConductivity(); + + /*retrieve some parameters: */ + this->parameters->FindParam(&dt,DtEnum); + this->parameters->FindParam(&artdiff,ArtDiffEnum); + this->parameters->FindParam(&epsvel,EpsVelEnum); + + /* 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.: */ + + /*Get gaussian points: */ + area_order=2; + num_vert_gauss=2; + + GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss); + + /* Start looping on the number of gaussian points: */ + for (igarea=0; igareaGetParameterValue(&u, gauss_coord,VxEnum); + inputs->GetParameterValue(&v, gauss_coord,VyEnum); + inputs->GetParameterValue(&w, gauss_coord,VzEnum); + + D_scalar=gauss_weight*Jdet; + + if(dt){ + D_scalar=D_scalar*dt; + } + + D[0][0]=D_scalar*u;D[0][1]=0; D[0][2]=0; + D[1][0]=0; D[1][1]=D_scalar*v;D[1][2]=0; + D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar*w; + + /* Do the triple product B'*D*Bprime: */ + MatrixMultiply(&B_advec[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_advec[0][0],0); + MatrixMultiply(&tBD_advec[0][0],numdof,3,0,&Bprime_advec[0][0],3,numdof,0,&Ke_gaussian_advec[0][0],0); + + /*Transient: */ + if(dt){ + GetNodalFunctionsP1(&L[0], gauss_coord); + D_scalar=gauss_weight*Jdet; + D_scalar=D_scalar; + + /* Do the triple product L'*D*L: */ + MatrixMultiply(&L[0],numdof,1,0,&D_scalar,1,1,0,&tLD[0],0); + MatrixMultiply(&tLD[0],numdof,1,0,&L[0],1,numdof,0,&Ke_gaussian_transient[0][0],0); + } + else{ + for(i=0;iCreateKMatrixThermal(Kgg); + delete tria; + } +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorBalancedthickness {{{1*/ +void Penta::CreatePVectorBalancedthickness( Vec pg){ + + /*Collapsed formulation: */ + Tria* tria=NULL; + + /*flags: */ + bool onbed; + bool onwater; + + /*recover some inputs: */ + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + + /*If on water, skip: */ + if(onwater)return; + + /*Is this element on the bed? :*/ + if(!onbed)return; + + /*Spawn Tria element from the base of the Penta: */ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->CreatePVector(pg); + delete tria; + return; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorBalancedvelocities {{{1*/ +void Penta::CreatePVectorBalancedvelocities( Vec pg){ + + /*Collapsed formulation: */ + Tria* tria=NULL; + + /*flags: */ + bool onbed; + bool onwater; + + /*recover some inputs: */ + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + + /*If on water, skip: */ + if(onwater)return; + + /*Is this element on the bed? :*/ + if(!onbed)return; + + /*Spawn Tria element from the base of the Penta: */ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->CreatePVector(pg); + delete tria; + return; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorDiagnosticHoriz {{{1*/ +void Penta::CreatePVectorDiagnosticHoriz( Vec pg){ + + int i,j; + + /* node data: */ + const int numgrids=6; + const int NDOF2=2; + const int numdof=NDOF2*numgrids; + double xyz_list[numgrids][3]; + int doflist[numdof]; + int numberofdofspernode; + + /* parameters: */ + double slope[3]; //do not put 2! this goes into GetParameterDerivativeValue, which addresses slope[3] also! + double driving_stress_baseline; + double thickness; + + /* 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; + double gauss_coord[4]; + int order_area_gauss; + int num_vert_gauss; + int num_area_gauss; + int ig1,ig2; + double gauss_weight1,gauss_weight2; + double gauss_weight; + + /* Jacobian: */ + double Jdet; + + /*nodal functions: */ + double l1l6[6]; + + /*element vector at the gaussian points: */ + double pe_g[numdof]={0.0}; + double pe_g_gaussian[numdof]; + + /*Spawning: */ + Tria* tria=NULL; + + /*inputs: */ + bool onwater; + bool collapse; + bool onbed; + Input* surface_input=NULL; + Input* thickness_input=NULL; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&collapse,CollapseEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + + /*If on water, skip load: */ + 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 load vector. */ + + if ((collapse==1) && (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 ((collapse==1) && (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; + return; + } + else{ + + /*Implement standard penta element: */ + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /*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;iGetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&collapse,CollapseEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); + + /*If on water, skip: */ + if(onwater) return; + + /*Bail out if this element if: + * -> Non collapsed and not on the surface + * -> collapsed (2d model) and not on bed) */ + if ((!collapse && !onsurface) || (collapse && !onbed)){ + return; + } + else if (collapse){ + + /*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; + 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; + return; + } +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorDiagnosticHutter{{{1*/ +void Penta::CreatePVectorDiagnosticHutter(Vec pg){ + + /*Collapsed formulation: */ + Beam* beam=NULL; + int i; + + /*flags: */ + bool onwater; + + /*recover some inputs: */ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + + /*If on water, skip: */ + if(onwater)return; + + /*Spawn 3 beam elements: */ + for(i=0;i<3;i++){ + beam=(Beam*)SpawnBeam(i,i+3); //[0 3], [1 4] and [2 5] are the four vertical edges of the Penta + beam->CreatePVector(pg); + } + + /*clean up*/ + delete beam; + +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorDiagnosticStokes {{{1*/ +void Penta::CreatePVectorDiagnosticStokes( Vec pg){ + + /*indexing: */ + int i,j; + + const int numgrids=6; + const int DOFPERGRID=4; + const int numdof=numgrids*DOFPERGRID; + const int numgrids2d=3; + int numdof2d=numgrids2d*DOFPERGRID; + int doflist[numdof]; + int numberofdofspernode; + + /*Material properties: */ + double gravity,rho_ice,rho_water; + + /*parameters: */ + double xyz_list_tria[numgrids2d][3]; + double xyz_list[numgrids][3]; + double surface_normal[3]; + double bed_normal[3]; + double bed; + + /* gaussian points: */ + int num_area_gauss; + int igarea,igvert; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* vert_gauss_coord = NULL; + double* area_gauss_weights = NULL; + double* vert_gauss_weights = NULL; + + /* specific gaussian point: */ + double gauss_weight,area_gauss_weight,vert_gauss_weight; + double gauss_coord[4]; + double gauss_coord_tria[3]; + + int area_order=5; + int num_vert_gauss=5; + + double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + double viscosity; + double water_pressure; + int dofs[3]={0,1,2}; + + /*matrices: */ + double Pe_temp[27]={0.0}; //for the six nodes and the bubble + double Pe_gaussian[27]={0.0}; //for the six nodes and the bubble + double Ke_temp[27][3]={0.0}; //for the six nodes and the bubble + double Pe_reduced[numdof]; //for the six nodes only + double Ke_gaussian[27][3]; + double L[3]; //for the three nodes of the bed + double l1l7[7]; //for the six nodes and the bubble + double B[8][27]; + double B_prime[8][27]; + double B_prime_bubble[8][3]; + double Jdet; + double Jdet2d; + double D[8][8]={0.0}; + double D_scalar; + double tBD[27][8]; + double P_terms[numdof]={0.0}; + + Tria* tria=NULL; + + /*parameters: */ + double stokesreconditioning; + + /*inputs: */ + bool onwater; + bool onbed; + bool shelf; + bool isstokes; + Input* vx_input=NULL; + Input* vy_input=NULL; + Input* vz_input=NULL; + Input* bed_input=NULL; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + inputs->GetParameterValue(&isstokes,IsStokesEnum); + + /*retrieve some parameters: */ + this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); + + /*If on water or not Stokes, skip load: */ + if(onwater || !isstokes) return; + + /*If on water, skip load: */ + if(onwater)return; + + /*recovre material parameters: */ + rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetRhoIce(); + gravity=matpar->GetG(); + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /* 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.*/ + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss); + + /*Retrieve all inputs we will be needing: */ + vx_input=inputs->GetInput(VxEnum); + vy_input=inputs->GetInput(VyEnum); + vz_input=inputs->GetInput(VzEnum); + bed_input=inputs->GetInput(BedEnum); + + /* Start looping on the number of gaussian points: */ + for (igarea=0; igareaGetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input); + matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); + + /* Get nodal functions */ + GetNodalFunctionsMINI(&l1l7[0], gauss_coord); + + /* Build gaussian vector */ + for(i=0;iGetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria); + + /* Get bed at gaussian point */ + bed_input->GetParameterValue(&bed, gauss_coord); + + /*Get L matrix : */ + tria->GetL(&L[0], &xyz_list[0][0], gauss_coord_tria,1); + + /*Get water_pressure at gaussian point */ + water_pressure=gravity*rho_water*bed; + + /*Get normal vecyor to the bed */ + SurfaceNormal(&surface_normal[0],xyz_list_tria); + + bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result + bed_normal[1]=-surface_normal[1]; + bed_normal[2]=-surface_normal[2]; + + for(i=0;iGetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); + + /*If on water, skip: */ + if(onwater || !onsurface)return; + + /*Call Tria's function*/ + tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). + tria->CreatePVectorAdjointStokes(p_g); + delete tria; + return; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorDiagnosticVert {{{1*/ +void Penta::CreatePVectorDiagnosticVert( Vec pg){ + + int i; + + /* node data: */ + const int numgrids=6; + const int NDOF1=1; + const int numdof=NDOF1*numgrids; + double xyz_list[numgrids][3]; + int doflist[numdof]; + int numberofdofspernode; + + /* 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; + double gauss_coord[4]; + int order_area_gauss; + int num_vert_gauss; + int num_area_gauss; + int ig1,ig2; + double gauss_weight1,gauss_weight2; + double gauss_weight; + + /* Jacobian: */ + double Jdet; + + /*element vector at the gaussian points: */ + double pe_g[numdof]={0.0}; + double pe_g_gaussian[numdof]; + double l1l6[6]; + + /*Spawning: */ + Tria* tria=NULL; + + /*input parameters for structural analysis (diagnostic): */ + double du[3]; + double dv[3]; + double dudx,dvdy; + int dofs1[1]={0}; + int dofs2[1]={1}; + + /*inputs: */ + bool onwater; + bool onbed; + Input* vx_input=NULL; + Input* vy_input=NULL; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + + /*If on water, skip: */ + if(onwater)return; + + /*If we are on the bedrock, spawn a tria on the bedrock, and use it to build the + *diagnostic base vertical stifness: */ + if(onbed){ + tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 are on the bedrock + tria->CreatePVectorDiagnosticBaseVert(pg); + delete tria; + } + + /*Now, handle the standard penta element: */ + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /*Get gaussian points and weights :*/ + order_area_gauss=2; + num_vert_gauss=2; + + 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); + + /* Start looping on the number of gaussian points: */ + for (ig1=0; ig1GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss_coord); + vy_input->GetParameterDerivativeValue(&dv[0],&xyz_list[0][0],gauss_coord); + dudx=du[0]; + dvdy=dv[1]; + + /* Get Jacobian determinant: */ + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); + + /*Get nodal functions: */ + GetNodalFunctionsP1(l1l6, gauss_coord); + + /*Build pe_g_gaussian vector: */ + for (i=0;iGetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); - /*go into parameters and get the analysis_counter: */ - this->parameters->FindParam(&analysis_counter,AnalysisCounterEnum); - - /*out of grids g0,g1 and g2 from Penta, build a tria element: */ - Tria* tria=NULL; - int indices[3]; - int zero=0; - Parameters* tria_parameters=NULL; - Inputs* tria_inputs=NULL; - Results* tria_results=NULL; - - indices[0]=g0; - indices[1]=g1; - indices[2]=g2; - - tria_parameters=this->parameters; - tria_inputs=(Inputs*)this->inputs->SpawnTriaInputs(indices); - tria_results=(Results*)this->results->SpawnTriaResults(indices); - - tria=new Tria(); - tria->id=this->id; - tria->inputs=tria_inputs; - tria->results=tria_results; - tria->parameters=tria_parameters; - tria->element_type=P1Enum; //Only P1 CG for now - this->SpawnTriaHook(dynamic_cast(tria),&indices[0]); - - /*recover nodes, matice and matpar: */ - tria->nodes=(Node**)tria->hnodes[analysis_counter]->deliverp(); - tria->matice=(Matice*)tria->hmatice->delivers(); - tria->matpar=(Matpar*)tria->hmatpar->delivers(); - - return tria; -} -/*}}}*/ -/*FUNCTION Penta::SpawnBeam {{{1*/ -void* Penta::SpawnBeam(int g0, int g1){ - - int i; - - /*out of grids g0,g1 and g2 from Penta, build a beam element: */ - Beam* beam=NULL; - int indices[2]; - int zero=0; - Parameters *beam_parameters = NULL; - Inputs *beam_inputs = NULL; - - indices[0]=g0; - indices[1]=g1; - - beam_parameters=this->parameters; - beam_inputs=(Inputs*)this->inputs->SpawnBeamInputs(indices); - - beam=new Beam(); - beam->id=this->id; - beam->inputs=beam_inputs; - beam->parameters=beam_parameters; - - /*now deal with ndoes,matice and matpar: */ - beam->nodes=(Node**)xmalloc(2*sizeof(Node*)); - for(i=0;i<2;i++)beam->nodes[i]=this->nodes[indices[i]]; - beam->matice=this->matice; - beam->matpar=this->matpar; - - return beam; -} -/*}}}*/ -/*FUNCTION Penta::SpawnSing {{{1*/ -void* Penta::SpawnSing(int index){ - - Sing* sing=NULL; - int zero=0; - Parameters *sing_parameters = NULL; - Inputs *sing_inputs = NULL; - - sing_parameters=this->parameters; - sing_inputs=(Inputs*)this->inputs->SpawnSingInputs(index); - - sing=new Sing(); - sing->id=this->id; - sing->inputs=sing_inputs; - sing->parameters=sing_parameters; - - /*now deal with nodes,matice and matpar: */ - sing->node=this->nodes[index]; - sing->matice=this->matice; - sing->matpar=this->matpar; - - return sing; + /*If on water, skip: */ + if(onwater)return; + + /*Is this element on the bed? :*/ + if(!onbed)return; + + /*Spawn Tria element from the base of the Penta: */ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->CreatePVector(pg); + delete tria; + return; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorSlope {{{1*/ + +void Penta::CreatePVectorSlope( Vec pg){ + + /*Collapsed formulation: */ + 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; + + /*Is this element on the bed? :*/ + if(!onbed)return; + + /*Spawn Tria element from the base of the Penta: */ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->CreatePVector(pg); + delete tria; + return; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorThermal {{{1*/ +void Penta::CreatePVectorThermal( Vec pg){ + + /*indexing: */ + int i,j; + int found=0; + + const int numgrids=6; + const int NDOF1=1; + const int numdof=numgrids*NDOF1; + int doflist[numdof]; + int numberofdofspernode; + + /*Grid data: */ + double xyz_list[numgrids][3]; + + /* gaussian points: */ + int num_area_gauss,igarea,igvert; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* vert_gauss_coord = NULL; + double* area_gauss_weights = NULL; + double* vert_gauss_weights = NULL; + double gauss_weight,area_gauss_weight,vert_gauss_weight; + double gauss_coord[4]; + int area_order=2; + int num_vert_gauss=3; + + double temperature_list[numgrids]; + double temperature; + + /*Material properties: */ + double gravity,rho_ice,rho_water; + double mixed_layer_capacity,heatcapacity; + double beta,meltingpoint,thermal_exchange_velocity; + + /* element parameters: */ + int friction_type; + + int dofs[3]={0,1,2}; + int dofs1[1]={0}; + + /*matrices: */ + double P_terms[numdof]={0.0}; + double L[numdof]; + double l1l2l3[3]; + double alpha2_list[3]; + double basalfriction_list[3]={0.0}; + double basalfriction; + double epsilon[6]; + double epsilon_sqr[3][3]; + double epsilon_matrix[3][3]; + + double Jdet; + double viscosity; + double epsilon_eff; + double phi; + double t_pmp; + double scalar; + double scalar_def; + double scalar_ocean; + double scalar_transient; + + /*Collapsed formulation: */ + Tria* tria=NULL; + + /*parameters: */ + double dt; + + /*inputs: */ + bool onwater; + bool onbed; + bool shelf; + Input* vx_input=NULL; + Input* vy_input=NULL; + Input* vz_input=NULL; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + + /*retrieve some parameters: */ + this->parameters->FindParam(&dt,DtEnum); + + /*If on water, skip: */ + if(onwater)return; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /*recovre material parameters: */ + rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetRhoIce(); + gravity=matpar->GetG(); + heatcapacity=matpar->GetHeatCapacity(); + beta=matpar->GetBeta(); + meltingpoint=matpar->GetMeltingPoint(); + + /* 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.: */ + + /*Get gaussian points: */ + GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss); + + /*Retrieve all inputs we will be needing: */ + vx_input=inputs->GetInput(VxEnum); + vy_input=inputs->GetInput(VyEnum); + vz_input=inputs->GetInput(VzEnum); + + /* Start looping on the number of gaussian points: */ + for (igarea=0; igareaGetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input); + matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); + + /* Get nodal functions */ + GetNodalFunctionsP1(&L[0], gauss_coord); + + /*Build deformational heating: */ + GetPhi(&phi, &epsilon[0], viscosity); + + /*Build pe_gaussian */ + scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss_weight; + if(dt){ + scalar_def=scalar_def*dt; + } + + for(i=0;iGetParameterValue(&temperature, gauss_coord,TemperatureEnum); + scalar_transient=temperature*Jdet*gauss_weight; + for(i=0;iCreatePVectorThermalShelf(pg); + delete tria; + } + + /* Geothermal flux on ice sheet base and basal friction */ + if(onbed && !shelf){ + + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->CreatePVectorThermalSheet(pg); + delete tria; + } + extern int my_rank; + + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&vert_gauss_coord); + xfree((void**)&area_gauss_weights); + xfree((void**)&vert_gauss_weights); + } /*}}}*/ @@ -2104,4 +3976,539 @@ } /*}}}*/ +/*FUNCTION Penta::GetDofList {{{1*/ +void Penta::GetDofList(int* doflist,int* pnumberofdofspernode){ + + int i,j; + int doflist_per_node[MAXDOFSPERNODE]; + int numberofdofspernode; + + for(i=0;i<6;i++){ + nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode); + for(j=0;jGetDofList1(); + } + +} +/*}}}*/ +/*FUNCTION Penta::GetElementType {{{1*/ +int Penta::GetElementType(){ + + /*return PentaRef field*/ + return this->element_type; + +} +/*}}}*/ +/*FUNCTION Penta::GetMatrixInvert {{{1*/ +void Penta::GetMatrixInvert(double* Ke_invert, double* Ke){ + /*Inverse a 3 by 3 matrix only */ + + double a,b,c,d,e,f,g,h,i; + double det; + int calculationdof=3; + + /*Take the matrix components: */ + a=*(Ke+calculationdof*0+0); + b=*(Ke+calculationdof*0+1); + c=*(Ke+calculationdof*0+2); + d=*(Ke+calculationdof*1+0); + e=*(Ke+calculationdof*1+1); + f=*(Ke+calculationdof*1+2); + g=*(Ke+calculationdof*2+0); + h=*(Ke+calculationdof*2+1); + i=*(Ke+calculationdof*2+2); + + det=a*(e*i-f*h)-b*(d*i-f*g)+c*(d*h-e*g); + + *(Ke_invert+calculationdof*0+0)=(e*i-f*h)/det; + *(Ke_invert+calculationdof*0+1)=(c*h-b*i)/det; + *(Ke_invert+calculationdof*0+2)=(b*f-c*e)/det; + *(Ke_invert+calculationdof*1+0)=(f*g-d*i)/det; + *(Ke_invert+calculationdof*1+1)=(a*i-c*g)/det; + *(Ke_invert+calculationdof*1+2)=(c*d-a*f)/det; + *(Ke_invert+calculationdof*2+0)=(d*h-e*g)/det; + *(Ke_invert+calculationdof*2+1)=(b*g-a*h)/det; + *(Ke_invert+calculationdof*2+2)=(a*e-b*d)/det; + +} +/*}}}*/ +/*FUNCTION Penta::GetParameterValue(double* pvalue, double* v_list,double* gauss_coord) {{{1*/ +void Penta::GetParameterValue(double* pvalue, double* v_list,double* gauss_coord){ + + const int numgrids=6; + double l1l6[numgrids]; + + GetNodalFunctionsP1(&l1l6[0], gauss_coord); + + *pvalue=l1l6[0]*v_list[0]+l1l6[1]*v_list[1]+l1l6[2]*v_list[2]+l1l6[3]*v_list[3]+l1l6[4]*v_list[4]+l1l6[5]*v_list[5]; +} +/*}}}*/ +/*FUNCTION Penta::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype) {{{1*/ +void Penta::GetParameterValue(double* pvalue,Node* node,int enumtype){ + + /*Output*/ + double value; + + /*Intermediaries*/ + const int numnodes=6; + int grid=-1; + int i; + double gauss[numnodes][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}}; + + /*go through 3 nodes (all nodes for tria) and identify 1st and 2nd nodes: */ + ISSMASSERT(nodes); + for(i=0;iGetParameterValue(pvalue,&gauss[grid][0],enumtype); + return; + +} +/*}}}*/ +/*FUNCTION Penta::GetPhi {{{1*/ +void Penta::GetPhi(double* phi, double* epsilon, double viscosity){ + /*Compute deformational heating from epsilon and viscosity */ + + double epsilon_matrix[3][3]; + double epsilon_eff; + double epsilon_sqr[3][3]; + + /* Build epsilon matrix */ + epsilon_matrix[0][0]=*(epsilon+0); + epsilon_matrix[1][0]=*(epsilon+3); + epsilon_matrix[2][0]=*(epsilon+4); + epsilon_matrix[0][1]=*(epsilon+3); + epsilon_matrix[1][1]=*(epsilon+1); + epsilon_matrix[2][1]=*(epsilon+5); + epsilon_matrix[0][2]=*(epsilon+4); + epsilon_matrix[1][2]=*(epsilon+5); + epsilon_matrix[2][2]=*(epsilon+2); + + /* Effective value of epsilon_matrix */ + epsilon_sqr[0][0]=pow(epsilon_matrix[0][0],2); + epsilon_sqr[1][0]=pow(epsilon_matrix[1][0],2); + epsilon_sqr[2][0]=pow(epsilon_matrix[2][0],2); + epsilon_sqr[0][1]=pow(epsilon_matrix[0][1],2); + epsilon_sqr[1][1]=pow(epsilon_matrix[1][1],2); + epsilon_sqr[2][1]=pow(epsilon_matrix[2][1],2); + epsilon_sqr[0][2]=pow(epsilon_matrix[0][2],2); + epsilon_sqr[1][2]=pow(epsilon_matrix[1][2],2); + epsilon_sqr[2][2]=pow(epsilon_matrix[2][2],2); + + epsilon_eff=1/pow(2,0.5)*pow((epsilon_sqr[0][0]+epsilon_sqr[0][1]+ epsilon_sqr[0][2]+ epsilon_sqr[1][0]+ epsilon_sqr[1][1]+ epsilon_sqr[1][2]+ epsilon_sqr[2][0]+ epsilon_sqr[2][1]+ epsilon_sqr[2][2]),0.5); + *phi=2*pow(epsilon_eff,2.0)*viscosity; +} +/*}}}*/ +/*FUNCTION Penta::GetSolutionFromInputsDiagnosticHoriz{{{1*/ +void Penta::GetSolutionFromInputsDiagnosticHoriz(Vec solution){ + + int i; + + const int numvertices=6; + const int numdofpervertex=2; + const int numdof=numdofpervertex*numvertices; + double gauss[numvertices][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 doflist[numdof]; + double values[numdof]; + double vx; + double vy; + + int dummy; + + /*Get dof list: */ + GetDofList(&doflist[0],&dummy); + + /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ + /*P1 element only for now*/ + for(i=0;iGetParameterValue(&vx,&gauss[i][0],VxEnum); + inputs->GetParameterValue(&vy,&gauss[i][0],VyEnum); + values[i*numdofpervertex+0]=vx; + values[i*numdofpervertex+1]=vy; + } + + /*Add value to global vector*/ + VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); + +} +/*}}}*/ +/*FUNCTION Penta::GetSolutionFromInputsDiagnosticHutter{{{1*/ +void Penta::GetSolutionFromInputsDiagnosticHutter(Vec solution){ + + int i; + + const int numvertices=6; + const int numdofpervertex=2; + const int numdof=numdofpervertex*numvertices; + double gauss[numvertices][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 doflist[numdof]; + double values[numdof]; + double vx; + double vy; + + int dummy; + + /*Get dof list: */ + GetDofList(&doflist[0],&dummy); + + /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ + /*P1 element only for now*/ + for(i=0;iGetParameterValue(&vx,&gauss[i][0],VxEnum); + inputs->GetParameterValue(&vy,&gauss[i][0],VyEnum); + values[i*numdofpervertex+0]=vx; + values[i*numdofpervertex+1]=vy; + } + + /*Add value to global vector*/ + VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); + +} +/*}}}*/ +/*FUNCTION Penta::GetSolutionFromInputsDiagnosticVert{{{1*/ +void Penta::GetSolutionFromInputsDiagnosticVert(Vec solution){ + + int i; + + const int numvertices=6; + const int numdofpervertex=1; + const int numdof=numdofpervertex*numvertices; + double gauss[numvertices][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 doflist[numdof]; + double values[numdof]; + double vz; + + int dummy; + + /*Get dof list: */ + GetDofList(&doflist[0],&dummy); + + /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ + /*P1 element only for now*/ + for(i=0;iGetParameterValue(&vz,&gauss[i][0],VxEnum); + values[i]=vz; + } + + /*Add value to global vector*/ + VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); +} +/*}}}*/ +/*FUNCTION Penta::GetSolutionFromInputsDiagnosticStokes{{{1*/ +void Penta::GetSolutionFromInputsDiagnosticStokes(Vec solution){ + + int i; + + const int numvertices=6; + const int numdofpervertex=4; + const int numdof=numdofpervertex*numvertices; + double gauss[numvertices][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 doflist[numdof]; + double values[numdof]; + double vx,vy,vz,p; + + int dummy; + double stokesreconditioning; + + /*Get dof list: */ + GetDofList(&doflist[0],&dummy); + + /*Recondition pressure: */ + this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); + + /*Ok, we have vx vy vz and P in values, fill in vx vy vz P arrays: */ + /*P1 element only for now*/ + for(i=0;iGetParameterValue(&vx,&gauss[i][0],VxEnum); + inputs->GetParameterValue(&vy,&gauss[i][0],VyEnum); + inputs->GetParameterValue(&vz,&gauss[i][0],VzEnum); + inputs->GetParameterValue(&p ,&gauss[i][0],PressureEnum); + values[i*numdofpervertex+0]=vx; + values[i*numdofpervertex+1]=vy; + values[i*numdofpervertex+2]=vz; + values[i*numdofpervertex+3]=p/stokesreconditioning; + } + + /*Add value to global vector*/ + VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); + +} +/*}}}*/ +/*FUNCTION Penta::GetSolutionFromInputsThermal{{{1*/ +void Penta::GetSolutionFromInputsThermal(Vec solution){ + + int i; + + const int numvertices=6; + const int numdofpervertex=1; + const int numdof=numdofpervertex*numvertices; + double gauss[numvertices][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 doflist[numdof]; + double values[numdof]; + double vz; + + int dummy; + + /*Get dof list: */ + GetDofList(&doflist[0],&dummy); + + /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ + /*P1 element only for now*/ + for(i=0;iGetParameterValue(&vz,&gauss[i][0],TemperatureEnum); + values[i]=vz; + } + + /*Add value to global vector*/ + VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); +} +/*}}}*/ +/*FUNCTION Penta::GetStrainRate3dPattyn{{{1*/ +void Penta::GetStrainRate3dPattyn(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input){ + /*Compute the 3d Blatter/PattynStrain Rate (5 components): + * + * epsilon=[exx eyy exy exz eyz] + * + * with exz=1/2 du/dz + * eyz=1/2 dv/dz + * + * the contribution of vz is neglected + */ + + int i; + + double epsilonvx[5]; + double epsilonvy[5]; + + /*Check that both inputs have been found*/ + if (!vx_input || !vy_input){ + ISSMERROR("Input missing. Here are the input pointers we have for vx: %p, vy: %p\n",vx_input,vy_input); + } + + /*Get strain rate assuming that epsilon has been allocated*/ + vx_input->GetVxStrainRate3dPattyn(epsilonvx,xyz_list,gauss); + vy_input->GetVyStrainRate3dPattyn(epsilonvy,xyz_list,gauss); + + /*Sum all contributions*/ + for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]; + +} +/*}}}*/ +/*FUNCTION Penta::GetStrainRate3d{{{1*/ +void Penta::GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input, Input* vz_input){ + /*Compute the 3d Strain Rate (6 components): + * + * epsilon=[exx eyy ezz exy exz eyz] + */ + + int i; + + double epsilonvx[6]; + double epsilonvy[6]; + double epsilonvz[6]; + + /*Check that both inputs have been found*/ + if (!vx_input || !vy_input || !vz_input){ + ISSMERROR("Input missing. Here are the input pointers we have for vx: %p, vy: %p, vz: %p\n",vx_input,vy_input,vz_input); + } + + /*Get strain rate assuming that epsilon has been allocated*/ + vx_input->GetVxStrainRate3d(epsilonvx,xyz_list,gauss); + vy_input->GetVyStrainRate3d(epsilonvy,xyz_list,gauss); + vz_input->GetVzStrainRate3d(epsilonvz,xyz_list,gauss); + + /*Sum all contributions*/ + for(i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i]; + +} +/*}}}*/ +/*FUNCTION Penta::GetUpperElement{{{1*/ +Penta* Penta::GetUpperElement(void){ + + Penta* upper_penta=NULL; + upper_penta=(Penta*)neighbors[1]; //first one under, second one above + return upper_penta; + +} +/*}}}*/ +/*FUNCTION Penta::GradjB {{{1*/ +void Penta::GradjB(Vec gradient){ + + int i; + Tria* tria=NULL; + TriaVertexInput* triavertexinput=NULL; + + /*inputs: */ + bool onwater; + bool collapse; + bool onbed; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&collapse,CollapseEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + + /*If on water, skip: */ + if(onwater)return; + + if (collapse){ + /*Bail out element if collapsed (2d) and not on bed*/ + if (!onbed) return; + + /*This element should be collapsed into a tria element at its base. Create this tria element, + * and compute gardj*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). + tria->GradjB(gradient); + delete tria; + } + else{ + /*B is a 2d field, use MacAyeal(2d) gradient even if it is Stokes or Pattyn*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). + tria->GradjB(gradient); + delete tria; + } + + +} +/*}}}*/ +/*FUNCTION Penta::GradjDrag {{{1*/ +void Penta::GradjDrag(Vec gradient){ + + int i; + Tria* tria=NULL; + TriaVertexInput* triavertexinput=NULL; + double temp_gradient[6]={0,0,0,0,0,0}; + + /*inputs: */ + bool onwater; + bool onbed; + bool shelf; + int analysis_type; + + /*retrieve parameters: */ + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + + /*If on water, skip: */ + if(onwater)return; + + /*If on shelf, skip: */ + if(shelf)return; + + /*Bail out if this element does not touch the bedrock: */ + if (!onbed) return; + + if (analysis_type==AdjointHorizAnalysisEnum){ + + /*MacAyeal or Pattyn*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->GradjDrag(gradient); + delete tria; + } + else if (analysis_type==AdjointStokesAnalysisEnum){ + + /*Stokes*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->GradjDragStokes(gradient); + delete tria; + } + else ISSMERROR("analysis %s not supported yet",EnumAsString(analysis_type)); + + +} +/*}}}*/ +/*FUNCTION Penta::InputExtrude {{{1*/ +void Penta::InputExtrude(int enum_type,bool only_if_collapsed){ + + bool onbed,collapse=false; + Penta *penta = NULL; + Input *original_input = NULL; + + /*recover parameters: */ + if (only_if_collapsed) inputs->GetParameterValue(&collapse,CollapseEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + + /*First: if only_if_collapsed, check wether this penta is collapsed*/ + if (only_if_collapsed && !collapse) return; + + /*Are we on the base, not on the surface?:*/ + if(onbed){ + + /*OK, we are on bed. we will follow the steps: + * 1: find input and extrude it. + * 2: follow the upper element until we reach the surface + * 3: for each element, we will add a copy of the extruded input*/ + + /*Step1: Extrude the original input: */ + original_input=(Input*)this->inputs->GetInput(enum_type); + if(!original_input) ISSMERROR("%s%s"," could not find input with enum:",EnumAsString(enum_type)); + original_input->Extrude(); + + /*Stop if there is only one layer of element*/ + if (this->IsOnSurface()) return; + + /*Step 2: this input has been extruded for this element, now follow the upper element*/ + penta=this; + for(;;){ + + /* get upper Penta*/ + penta=penta->GetUpperElement(); + ISSMASSERT(penta->Id()!=this->id); + + /*Add input of the basal element to penta->inputs*/ + Input* copy=NULL; + copy=(Input*)original_input->copy(); + penta->inputs->AddInput((Input*)copy); + + /*Stop if we have reached the surface*/ + if (penta->IsOnSurface()) break; + + } + } + + return; +} +/*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticHoriz {{{1*/ void Penta::InputUpdateFromSolutionDiagnosticHoriz(double* solution){ @@ -2707,183 +5114,4 @@ penta=penta->GetUpperElement(); ISSMASSERT(penta->Id()!=this->id); } -} -/*}}}*/ -/*FUNCTION Penta::GetSolutionFromInputsDiagnosticHoriz{{{1*/ -void Penta::GetSolutionFromInputsDiagnosticHoriz(Vec solution){ - - int i; - - const int numvertices=6; - const int numdofpervertex=2; - const int numdof=numdofpervertex*numvertices; - double gauss[numvertices][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 doflist[numdof]; - double values[numdof]; - double vx; - double vy; - - int dummy; - - /*Get dof list: */ - GetDofList(&doflist[0],&dummy); - - /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ - /*P1 element only for now*/ - for(i=0;iGetParameterValue(&vx,&gauss[i][0],VxEnum); - inputs->GetParameterValue(&vy,&gauss[i][0],VyEnum); - values[i*numdofpervertex+0]=vx; - values[i*numdofpervertex+1]=vy; - } - - /*Add value to global vector*/ - VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); - -} -/*}}}*/ -/*FUNCTION Penta::GetSolutionFromInputsDiagnosticHutter{{{1*/ -void Penta::GetSolutionFromInputsDiagnosticHutter(Vec solution){ - - int i; - - const int numvertices=6; - const int numdofpervertex=2; - const int numdof=numdofpervertex*numvertices; - double gauss[numvertices][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 doflist[numdof]; - double values[numdof]; - double vx; - double vy; - - int dummy; - - /*Get dof list: */ - GetDofList(&doflist[0],&dummy); - - /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ - /*P1 element only for now*/ - for(i=0;iGetParameterValue(&vx,&gauss[i][0],VxEnum); - inputs->GetParameterValue(&vy,&gauss[i][0],VyEnum); - values[i*numdofpervertex+0]=vx; - values[i*numdofpervertex+1]=vy; - } - - /*Add value to global vector*/ - VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); - -} -/*}}}*/ -/*FUNCTION Penta::GetSolutionFromInputsDiagnosticVert{{{1*/ -void Penta::GetSolutionFromInputsDiagnosticVert(Vec solution){ - - int i; - - const int numvertices=6; - const int numdofpervertex=1; - const int numdof=numdofpervertex*numvertices; - double gauss[numvertices][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 doflist[numdof]; - double values[numdof]; - double vz; - - int dummy; - - /*Get dof list: */ - GetDofList(&doflist[0],&dummy); - - /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ - /*P1 element only for now*/ - for(i=0;iGetParameterValue(&vz,&gauss[i][0],VxEnum); - values[i]=vz; - } - - /*Add value to global vector*/ - VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); -} -/*}}}*/ -/*FUNCTION Penta::GetSolutionFromInputsDiagnosticStokes{{{1*/ -void Penta::GetSolutionFromInputsDiagnosticStokes(Vec solution){ - - int i; - - const int numvertices=6; - const int numdofpervertex=4; - const int numdof=numdofpervertex*numvertices; - double gauss[numvertices][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 doflist[numdof]; - double values[numdof]; - double vx,vy,vz,p; - - int dummy; - double stokesreconditioning; - - /*Get dof list: */ - GetDofList(&doflist[0],&dummy); - - /*Recondition pressure: */ - this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); - - /*Ok, we have vx vy vz and P in values, fill in vx vy vz P arrays: */ - /*P1 element only for now*/ - for(i=0;iGetParameterValue(&vx,&gauss[i][0],VxEnum); - inputs->GetParameterValue(&vy,&gauss[i][0],VyEnum); - inputs->GetParameterValue(&vz,&gauss[i][0],VzEnum); - inputs->GetParameterValue(&p ,&gauss[i][0],PressureEnum); - values[i*numdofpervertex+0]=vx; - values[i*numdofpervertex+1]=vy; - values[i*numdofpervertex+2]=vz; - values[i*numdofpervertex+3]=p/stokesreconditioning; - } - - /*Add value to global vector*/ - VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); - -} -/*}}}*/ -/*FUNCTION Penta::GetSolutionFromInputsThermal{{{1*/ -void Penta::GetSolutionFromInputsThermal(Vec solution){ - - int i; - - const int numvertices=6; - const int numdofpervertex=1; - const int numdof=numdofpervertex*numvertices; - double gauss[numvertices][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 doflist[numdof]; - double values[numdof]; - double vz; - - int dummy; - - /*Get dof list: */ - GetDofList(&doflist[0],&dummy); - - /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ - /*P1 element only for now*/ - for(i=0;iGetParameterValue(&vz,&gauss[i][0],TemperatureEnum); - values[i]=vz; - } - - /*Add value to global vector*/ - VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); } /*}}}*/ @@ -2926,3289 +5154,4 @@ } /*}}}*/ -/*FUNCTION Penta::GetUpperElement{{{1*/ -Penta* Penta::GetUpperElement(void){ - - Penta* upper_penta=NULL; - upper_penta=(Penta*)neighbors[1]; //first one under, second one above - return upper_penta; - -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixBalancedthickness {{{1*/ - -void Penta::CreateKMatrixBalancedthickness(Mat Kgg){ - - /*Collapsed formulation: */ - Tria* tria=NULL; - - /*flags: */ - bool onwater; - bool onbed; - - /*recover some inputs: */ - inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - - /*If on water, skip: */ - if(onwater)return; - - /*Is this element on the bed? :*/ - if(!onbed)return; - - /*Spawn Tria element from the base of the Penta: */ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrix(Kgg); - delete tria; - return; - -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixBalancedvelocities {{{1*/ - -void Penta::CreateKMatrixBalancedvelocities(Mat Kgg){ - - /*Collapsed formulation: */ - Tria* tria=NULL; - - /*flags: */ - bool onbed; - bool onwater; - - /*recover some inputs: */ - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - - /*If on water, skip: */ - if(onwater)return; - - /*Is this element on the bed? :*/ - if(!onbed)return; - - /*Spawn Tria element from the base of the Penta: */ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrix(Kgg); - delete tria; - 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[numdof]; - int numberofdofspernode; - - /* 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 collapse; - 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(&collapse,CollapseEnum); - 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; - - /*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 ((collapse==1) && (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 ((collapse==1) && (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: */ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrix(Kgg); - delete tria; - return; - } - else{ - - /*Implement standard penta element: */ - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /*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: */ - GetB(&B[0][0], &xyz_list[0][0], gauss_coord); - GetBPrime(&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; - } - - } - - 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); - -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixDiagnosticHutter{{{1*/ -void Penta::CreateKMatrixDiagnosticHutter(Mat Kgg){ - - /*Collapsed formulation: */ - Beam* beam=NULL; - int i; - - /*flags: */ - bool onwater; - - /*recover some inputs: */ - inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - - /*If on water, skip: */ - if(onwater)return; - - /*Spawn 3 beam elements: */ - for(i=0;i<3;i++){ - beam=(Beam*)SpawnBeam(i,i+3); //[0 3], [1 4] and [2 5] are the four vertical edges of the Penta - beam->CreateKMatrix(Kgg); - } - - /*clean up*/ - delete beam; - -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixDiagnosticStokes {{{1*/ -void Penta::CreateKMatrixDiagnosticStokes( Mat Kgg){ - - int i,j; - - const int numgrids=6; - const int DOFPERGRID=4; - const int numdof=numgrids*DOFPERGRID; - int doflist[numdof]; - int numberofdofspernode; - - const int numgrids2d=3; - const int numdof2d=numgrids2d*DOFPERGRID; - - int dofs[3]={0,1,2}; - - double K_terms[numdof][numdof]={0.0}; - - /*Material properties: */ - double gravity,rho_ice,rho_water; - - /*Collapsed formulation: */ - Tria* tria=NULL; - - /*Grid data: */ - double xyz_list[numgrids][3]; - - /*parameters: */ - double xyz_list_tria[numgrids2d][3]; - double surface_normal[3]; - double bed_normal[3]; - double thickness; - - /*matrices: */ - double Ke_temp[27][27]={0.0}; //for the six nodes and the bubble - double Ke_reduced[numdof][numdof]; //for the six nodes only - double Ke_gaussian[27][27]; - double Ke_drag_gaussian[numdof2d][numdof2d]; - double B[8][27]; - double B_prime[8][27]; - double LStokes[14][numdof2d]; - double LprimeStokes[14][numdof2d]; - double Jdet; - double Jdet2d; - double D[8][8]={0.0}; - double D_scalar; - double tBD[27][8]; - double DLStokes[14][14]={0.0}; - double tLDStokes[numdof2d][14]; - - /* gaussian points: */ - int num_area_gauss; - int igarea,igvert; - double* first_gauss_area_coord = NULL; - double* second_gauss_area_coord = NULL; - double* third_gauss_area_coord = NULL; - double* vert_gauss_coord = NULL; - double* area_gauss_weights = NULL; - double* vert_gauss_weights = NULL; - double gaussgrids[numgrids][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}}; - - /* specific gaussian point: */ - double gauss_weight,area_gauss_weight,vert_gauss_weight; - double gauss_coord[4]; - double gauss_coord_tria[3]; - int area_order=5; - int num_vert_gauss=5; - - double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - double viscosity; - double alpha2_gauss; - Friction* friction=NULL; - - /*parameters: */ - double stokesreconditioning; - int analysis_type; - - /*inputs: */ - bool onwater; - bool onbed; - bool shelf; - bool isstokes; - - /*inputs: */ - Input* vx_input=NULL; - Input* vy_input=NULL; - Input* vz_input=NULL; - - /*retrive parameters: */ - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*retrieve inputs :*/ - inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); - inputs->GetParameterValue(&isstokes,IsStokesEnum); - - /*If on water or not Stokes, skip stiffness: */ - if(onwater || !isstokes) return; - - /*recovre material parameters: */ - rho_water=matpar->GetRhoWater(); - rho_ice=matpar->GetRhoIce(); - gravity=matpar->GetG(); - - /*retrieve some parameters: */ - this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /* 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.*/ - - area_order=5; - num_vert_gauss=5; - - /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ - GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss); - - /*Retrieve all inputs we will be needing: */ - vx_input=inputs->GetInput(VxEnum); - vy_input=inputs->GetInput(VyEnum); - vz_input=inputs->GetInput(VzEnum); - - /* Start looping on the number of gaussian points: */ - for (igarea=0; igareaGetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input); - - /*Get viscosity: */ - matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); - - /*Get B and Bprime matrices: */ - GetBStokes(&B[0][0],&xyz_list[0][0],gauss_coord); - GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0], gauss_coord); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],&gauss_coord[0]); - - /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant - * onto this scalar matrix, so that we win some computational time: */ - D_scalar=gauss_weight*Jdet; - for (i=0;i<6;i++){ - D[i][i]=D_scalar*viscosity; - } - for (i=6;i<8;i++){ - D[i][i]=-D_scalar*stokesreconditioning; - } - - /* Do the triple product tB*D*Bprime: */ - MatrixMultiply(&B[0][0],8,27,1,&D[0][0],8,8,0,&tBD[0][0],0); - MatrixMultiply(&tBD[0][0],27,8,0,&B_prime[0][0],8,27,0,&Ke_gaussian[0][0],0); - - /*Add Ke_gaussian and Ke_gaussian to terms in pKe. Watch out for column orientation from matlab: */ - for(i=0;i<27;i++){ - for(j=0;j<27;j++){ - Ke_temp[i][j]+=Ke_gaussian[i][j]; - } - } - } - } - - if((onbed==1) && (shelf==0)){ - - /*build friction object, used later on: */ - friction=new Friction("3d",inputs,matpar,analysis_type); - - for(i=0;iGetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria); - - /*Get L matrix if viscous basal drag present: */ - GetLStokes(&LStokes[0][0], gauss_coord_tria); - GetLprimeStokes(&LprimeStokes[0][0], &xyz_list[0][0], gauss_coord_tria, gauss_coord); - - /*Compute strain rate: */ - this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input); - - /*Get viscosity at last iteration: */ - matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); - - /*Get normal vecyor to the bed */ - SurfaceNormal(&surface_normal[0],xyz_list_tria); - - bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result - bed_normal[1]=-surface_normal[1]; - bed_normal[2]=-surface_normal[2]; - - /*Calculate DL on gauss point */ - friction->GetAlpha2(&alpha2_gauss, gauss_coord,VxEnum,VyEnum,VzEnum); - - DLStokes[0][0]=alpha2_gauss*gauss_weight*Jdet2d; - DLStokes[1][1]=alpha2_gauss*gauss_weight*Jdet2d; - DLStokes[2][2]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[0]*bed_normal[2]; - DLStokes[3][3]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[1]*bed_normal[2]; - DLStokes[4][4]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[0]*bed_normal[2]; - DLStokes[5][5]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[1]*bed_normal[2]; - DLStokes[6][6]=-viscosity*gauss_weight*Jdet2d*bed_normal[0]; - DLStokes[7][7]=-viscosity*gauss_weight*Jdet2d*bed_normal[1]; - DLStokes[8][8]=-viscosity*gauss_weight*Jdet2d*bed_normal[2]; - DLStokes[9][8]=-viscosity*gauss_weight*Jdet2d*bed_normal[0]/2.0; - DLStokes[10][10]=-viscosity*gauss_weight*Jdet2d*bed_normal[1]/2.0; - DLStokes[11][11]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[0]; - DLStokes[12][12]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[1]; - DLStokes[13][13]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[2]; - - /* Do the triple product tL*D*L: */ - MatrixMultiply(&LStokes[0][0],14,numdof2d,1,&DLStokes[0][0],14,14,0,&tLDStokes[0][0],0); - MatrixMultiply(&tLDStokes[0][0],numdof2d,14,0,&LprimeStokes[0][0],14,numdof2d,0,&Ke_drag_gaussian[0][0],0); - - for(i=0;iGetParameterValue(&onwater,ElementOnWaterEnum); - inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); - - /*If on water, skip stiffness: */ - if(onwater)return; - - /*If this element is on the surface, we have a dynamic boundary condition that applies, as a stiffness - * matrix: */ - if(onsurface){ - tria=(Tria*)SpawnTria(3,4,5); //nodes 3,4 and 5 are on the surface - tria->CreateKMatrixDiagnosticSurfaceVert(Kgg); - delete tria; - } - - /*Now, onto the formulation for the vertical velocity: */ - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /*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=2; - num_vert_gauss=2; - - 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); - - /* Start looping on the number of gaussian points: */ - for (ig1=0; ig1GetParameterValue(&onwater,ElementOnWaterEnum); - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - - /*If on water, skip: */ - if(onwater)return; - - if (!onbed){ - return; - } - else{ - - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrixMelting(Kgg); - delete tria; - return; - } -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixPrognostic {{{1*/ - -void Penta::CreateKMatrixPrognostic(Mat Kgg){ - - /*Collapsed formulation: */ - 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; - - /*Is this element on the bed? :*/ - if(!onbed)return; - - /*Spawn Tria element from the base of the Penta: */ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrix(Kgg); - delete tria; - return; - -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixSlope {{{1*/ - -void Penta::CreateKMatrixSlope(Mat Kgg){ - - /*Collapsed formulation: */ - 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; - - /*Is this element on the bed? :*/ - if(!onbed)return; - - /*Spawn Tria element from the base of the Penta: */ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrix(Kgg); - delete tria; - return; - -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixThermal {{{1*/ -void Penta::CreateKMatrixThermal(Mat Kgg){ - - /* local declarations */ - int i,j; - int found=0; - - /* node data: */ - const int numgrids=6; - const int NDOF1=1; - const int numdof=NDOF1*numgrids; - double xyz_list[numgrids][3]; - int doflist[numdof]; - int numberofdofspernode; - - /* gaussian points: */ - int num_area_gauss,igarea,igvert; - double* first_gauss_area_coord = NULL; - double* second_gauss_area_coord = NULL; - double* third_gauss_area_coord = NULL; - double* vert_gauss_coord = NULL; - double* area_gauss_weights = NULL; - double* vert_gauss_weights = NULL; - double gauss_weight,area_gauss_weight,vert_gauss_weight; - double gauss_coord[4]; - double gauss_l1l2l3[3]; - - int area_order=5; - int num_vert_gauss=5; - - int dofs[3]={0,1,2}; - double K[2][2]={0.0}; - - double u,v,w; - - /*matrices: */ - double K_terms[numdof][numdof]={0.0}; - double Ke_gaussian_conduct[numdof][numdof]; - double Ke_gaussian_advec[numdof][numdof]; - double Ke_gaussian_artdiff[numdof][numdof]; - double Ke_gaussian_transient[numdof][numdof]; - double B[3][numdof]; - double Bprime[3][numdof]; - double B_conduct[3][numdof]; - double B_advec[3][numdof]; - double B_artdiff[2][numdof]; - double Bprime_advec[3][numdof]; - double L[numdof]; - double D_scalar; - double D[3][3]; - double l1l2l3[3]; - double tl1l2l3D[3]; - double tBD[3][numdof]; - double tBD_conduct[3][numdof]; - double tBD_advec[3][numdof]; - double tBD_artdiff[3][numdof]; - double tLD[numdof]; - - double Jdet; - - /*Material properties: */ - double gravity,rho_ice,rho_water; - double heatcapacity,thermalconductivity; - double mixed_layer_capacity,thermal_exchange_velocity; - - /*parameters: */ - double dt,epsvel; - bool artdiff; - - /*Collapsed formulation: */ - Tria* tria=NULL; - - - /*inputs: */ - bool onwater; - bool onbed; - bool shelf; - - /*retrieve inputs :*/ - inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); - - /*If on water, skip: */ - if(onwater)return; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - // /*recovre material parameters: */ - rho_water=matpar->GetRhoWater(); - rho_ice=matpar->GetRhoIce(); - gravity=matpar->GetG(); - heatcapacity=matpar->GetHeatCapacity(); - thermalconductivity=matpar->GetThermalConductivity(); - - /*retrieve some parameters: */ - this->parameters->FindParam(&dt,DtEnum); - this->parameters->FindParam(&artdiff,ArtDiffEnum); - this->parameters->FindParam(&epsvel,EpsVelEnum); - - /* 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.: */ - - /*Get gaussian points: */ - area_order=2; - num_vert_gauss=2; - - GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss); - - /* Start looping on the number of gaussian points: */ - for (igarea=0; igareaGetParameterValue(&u, gauss_coord,VxEnum); - inputs->GetParameterValue(&v, gauss_coord,VyEnum); - inputs->GetParameterValue(&w, gauss_coord,VzEnum); - - D_scalar=gauss_weight*Jdet; - - if(dt){ - D_scalar=D_scalar*dt; - } - - D[0][0]=D_scalar*u;D[0][1]=0; D[0][2]=0; - D[1][0]=0; D[1][1]=D_scalar*v;D[1][2]=0; - D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar*w; - - /* Do the triple product B'*D*Bprime: */ - MatrixMultiply(&B_advec[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_advec[0][0],0); - MatrixMultiply(&tBD_advec[0][0],numdof,3,0,&Bprime_advec[0][0],3,numdof,0,&Ke_gaussian_advec[0][0],0); - - /*Transient: */ - if(dt){ - GetNodalFunctions(&L[0], gauss_coord); - D_scalar=gauss_weight*Jdet; - D_scalar=D_scalar; - - /* Do the triple product L'*D*L: */ - MatrixMultiply(&L[0],numdof,1,0,&D_scalar,1,1,0,&tLD[0],0); - MatrixMultiply(&tLD[0],numdof,1,0,&L[0],1,numdof,0,&Ke_gaussian_transient[0][0],0); - } - else{ - for(i=0;iCreateKMatrixThermal(Kgg); - delete tria; - } -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorBalancedthickness {{{1*/ -void Penta::CreatePVectorBalancedthickness( Vec pg){ - - /*Collapsed formulation: */ - Tria* tria=NULL; - - /*flags: */ - bool onbed; - bool onwater; - - /*recover some inputs: */ - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - - /*If on water, skip: */ - if(onwater)return; - - /*Is this element on the bed? :*/ - if(!onbed)return; - - /*Spawn Tria element from the base of the Penta: */ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVector(pg); - delete tria; - return; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorBalancedvelocities {{{1*/ -void Penta::CreatePVectorBalancedvelocities( Vec pg){ - - /*Collapsed formulation: */ - Tria* tria=NULL; - - /*flags: */ - bool onbed; - bool onwater; - - /*recover some inputs: */ - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - - /*If on water, skip: */ - if(onwater)return; - - /*Is this element on the bed? :*/ - if(!onbed)return; - - /*Spawn Tria element from the base of the Penta: */ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVector(pg); - delete tria; - return; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorDiagnosticHoriz {{{1*/ -void Penta::CreatePVectorDiagnosticHoriz( Vec pg){ - - int i,j; - - /* node data: */ - const int numgrids=6; - const int NDOF2=2; - const int numdof=NDOF2*numgrids; - double xyz_list[numgrids][3]; - int doflist[numdof]; - int numberofdofspernode; - - /* parameters: */ - double slope[3]; //do not put 2! this goes into GetParameterDerivativeValue, which addresses slope[3] also! - double driving_stress_baseline; - double thickness; - - /* 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; - double gauss_coord[4]; - int order_area_gauss; - int num_vert_gauss; - int num_area_gauss; - int ig1,ig2; - double gauss_weight1,gauss_weight2; - double gauss_weight; - - /* Jacobian: */ - double Jdet; - - /*nodal functions: */ - double l1l6[6]; - - /*element vector at the gaussian points: */ - double pe_g[numdof]={0.0}; - double pe_g_gaussian[numdof]; - - /*Spawning: */ - Tria* tria=NULL; - - /*inputs: */ - bool onwater; - bool collapse; - bool onbed; - Input* surface_input=NULL; - Input* thickness_input=NULL; - - /*retrieve inputs :*/ - inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - inputs->GetParameterValue(&collapse,CollapseEnum); - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - - /*If on water, skip load: */ - 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 load vector. */ - - if ((collapse==1) && (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 ((collapse==1) && (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; - return; - } - else{ - - /*Implement standard penta element: */ - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /*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: */ - GetNodalFunctions(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(&collapse,CollapseEnum); - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); - - /*If on water, skip: */ - if(onwater) return; - - /*Bail out if this element if: - * -> Non collapsed and not on the surface - * -> collapsed (2d model) and not on bed) */ - if ((!collapse && !onsurface) || (collapse && !onbed)){ - return; - } - else if (collapse){ - - /*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; - 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; - return; - } -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorDiagnosticHutter{{{1*/ -void Penta::CreatePVectorDiagnosticHutter(Vec pg){ - - /*Collapsed formulation: */ - Beam* beam=NULL; - int i; - - /*flags: */ - bool onwater; - - /*recover some inputs: */ - inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - - /*If on water, skip: */ - if(onwater)return; - - /*Spawn 3 beam elements: */ - for(i=0;i<3;i++){ - beam=(Beam*)SpawnBeam(i,i+3); //[0 3], [1 4] and [2 5] are the four vertical edges of the Penta - beam->CreatePVector(pg); - } - - /*clean up*/ - delete beam; - -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorDiagnosticStokes {{{1*/ -void Penta::CreatePVectorDiagnosticStokes( Vec pg){ - - /*indexing: */ - int i,j; - - const int numgrids=6; - const int DOFPERGRID=4; - const int numdof=numgrids*DOFPERGRID; - const int numgrids2d=3; - int numdof2d=numgrids2d*DOFPERGRID; - int doflist[numdof]; - int numberofdofspernode; - - /*Material properties: */ - double gravity,rho_ice,rho_water; - - /*parameters: */ - double xyz_list_tria[numgrids2d][3]; - double xyz_list[numgrids][3]; - double surface_normal[3]; - double bed_normal[3]; - double bed; - - /* gaussian points: */ - int num_area_gauss; - int igarea,igvert; - double* first_gauss_area_coord = NULL; - double* second_gauss_area_coord = NULL; - double* third_gauss_area_coord = NULL; - double* vert_gauss_coord = NULL; - double* area_gauss_weights = NULL; - double* vert_gauss_weights = NULL; - - /* specific gaussian point: */ - double gauss_weight,area_gauss_weight,vert_gauss_weight; - double gauss_coord[4]; - double gauss_coord_tria[3]; - - int area_order=5; - int num_vert_gauss=5; - - double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - double viscosity; - double water_pressure; - int dofs[3]={0,1,2}; - - /*matrices: */ - double Pe_temp[27]={0.0}; //for the six nodes and the bubble - double Pe_gaussian[27]={0.0}; //for the six nodes and the bubble - double Ke_temp[27][3]={0.0}; //for the six nodes and the bubble - double Pe_reduced[numdof]; //for the six nodes only - double Ke_gaussian[27][3]; - double L[3]; //for the three nodes of the bed - double l1l7[7]; //for the six nodes and the bubble - double B[8][27]; - double B_prime[8][27]; - double B_prime_bubble[8][3]; - double Jdet; - double Jdet2d; - double D[8][8]={0.0}; - double D_scalar; - double tBD[27][8]; - double P_terms[numdof]={0.0}; - - Tria* tria=NULL; - - /*parameters: */ - double stokesreconditioning; - - /*inputs: */ - bool onwater; - bool onbed; - bool shelf; - bool isstokes; - Input* vx_input=NULL; - Input* vy_input=NULL; - Input* vz_input=NULL; - Input* bed_input=NULL; - - /*retrieve inputs :*/ - inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); - inputs->GetParameterValue(&isstokes,IsStokesEnum); - - /*retrieve some parameters: */ - this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); - - /*If on water or not Stokes, skip load: */ - if(onwater || !isstokes) return; - - /*If on water, skip load: */ - if(onwater)return; - - /*recovre material parameters: */ - rho_water=matpar->GetRhoWater(); - rho_ice=matpar->GetRhoIce(); - gravity=matpar->GetG(); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /* 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.*/ - - /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ - GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss); - - /*Retrieve all inputs we will be needing: */ - vx_input=inputs->GetInput(VxEnum); - vy_input=inputs->GetInput(VyEnum); - vz_input=inputs->GetInput(VzEnum); - bed_input=inputs->GetInput(BedEnum); - - /* Start looping on the number of gaussian points: */ - for (igarea=0; igareaGetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input); - matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); - - /* Get nodal functions */ - GetNodalFunctionsStokes(&l1l7[0], gauss_coord); - - /* Build gaussian vector */ - for(i=0;iGetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria); - - /* Get bed at gaussian point */ - bed_input->GetParameterValue(&bed, gauss_coord); - - /*Get L matrix : */ - tria->GetL(&L[0], &xyz_list[0][0], gauss_coord_tria,1); - - /*Get water_pressure at gaussian point */ - water_pressure=gravity*rho_water*bed; - - /*Get normal vecyor to the bed */ - SurfaceNormal(&surface_normal[0],xyz_list_tria); - - bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result - bed_normal[1]=-surface_normal[1]; - bed_normal[2]=-surface_normal[2]; - - for(i=0;iGetParameterValue(&onwater,ElementOnWaterEnum); - inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); - - /*If on water, skip: */ - if(onwater || !onsurface)return; - - /*Call Tria's function*/ - tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). - tria->CreatePVectorAdjointStokes(p_g); - delete tria; - return; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorDiagnosticVert {{{1*/ -void Penta::CreatePVectorDiagnosticVert( Vec pg){ - - int i; - - /* node data: */ - const int numgrids=6; - const int NDOF1=1; - const int numdof=NDOF1*numgrids; - double xyz_list[numgrids][3]; - int doflist[numdof]; - int numberofdofspernode; - - /* 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; - double gauss_coord[4]; - int order_area_gauss; - int num_vert_gauss; - int num_area_gauss; - int ig1,ig2; - double gauss_weight1,gauss_weight2; - double gauss_weight; - - /* Jacobian: */ - double Jdet; - - /*element vector at the gaussian points: */ - double pe_g[numdof]={0.0}; - double pe_g_gaussian[numdof]; - double l1l6[6]; - - /*Spawning: */ - Tria* tria=NULL; - - /*input parameters for structural analysis (diagnostic): */ - double du[3]; - double dv[3]; - double dudx,dvdy; - int dofs1[1]={0}; - int dofs2[1]={1}; - - /*inputs: */ - bool onwater; - bool onbed; - Input* vx_input=NULL; - Input* vy_input=NULL; - - /*retrieve inputs :*/ - inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - - /*If on water, skip: */ - if(onwater)return; - - /*If we are on the bedrock, spawn a tria on the bedrock, and use it to build the - *diagnostic base vertical stifness: */ - if(onbed){ - tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 are on the bedrock - tria->CreatePVectorDiagnosticBaseVert(pg); - delete tria; - } - - /*Now, handle the standard penta element: */ - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /*Get gaussian points and weights :*/ - order_area_gauss=2; - num_vert_gauss=2; - - 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); - - /* Start looping on the number of gaussian points: */ - for (ig1=0; ig1GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss_coord); - vy_input->GetParameterDerivativeValue(&dv[0],&xyz_list[0][0],gauss_coord); - dudx=du[0]; - dvdy=dv[1]; - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); - - /*Get nodal functions: */ - GetNodalFunctions(l1l6, gauss_coord); - - /*Build pe_g_gaussian vector: */ - for (i=0;iGetParameterValue(&onwater,ElementOnWaterEnum); - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - - /*If on water, skip: */ - if(onwater)return; - - /*Is this element on the bed? :*/ - if(!onbed)return; - - /*Spawn Tria element from the base of the Penta: */ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVector(pg); - delete tria; - return; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorSlope {{{1*/ - -void Penta::CreatePVectorSlope( Vec pg){ - - /*Collapsed formulation: */ - 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; - - /*Is this element on the bed? :*/ - if(!onbed)return; - - /*Spawn Tria element from the base of the Penta: */ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVector(pg); - delete tria; - return; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorThermal {{{1*/ -void Penta::CreatePVectorThermal( Vec pg){ - - /*indexing: */ - int i,j; - int found=0; - - const int numgrids=6; - const int NDOF1=1; - const int numdof=numgrids*NDOF1; - int doflist[numdof]; - int numberofdofspernode; - - /*Grid data: */ - double xyz_list[numgrids][3]; - - /* gaussian points: */ - int num_area_gauss,igarea,igvert; - double* first_gauss_area_coord = NULL; - double* second_gauss_area_coord = NULL; - double* third_gauss_area_coord = NULL; - double* vert_gauss_coord = NULL; - double* area_gauss_weights = NULL; - double* vert_gauss_weights = NULL; - double gauss_weight,area_gauss_weight,vert_gauss_weight; - double gauss_coord[4]; - int area_order=2; - int num_vert_gauss=3; - - double temperature_list[numgrids]; - double temperature; - - /*Material properties: */ - double gravity,rho_ice,rho_water; - double mixed_layer_capacity,heatcapacity; - double beta,meltingpoint,thermal_exchange_velocity; - - /* element parameters: */ - int friction_type; - - int dofs[3]={0,1,2}; - int dofs1[1]={0}; - - /*matrices: */ - double P_terms[numdof]={0.0}; - double L[numdof]; - double l1l2l3[3]; - double alpha2_list[3]; - double basalfriction_list[3]={0.0}; - double basalfriction; - double epsilon[6]; - double epsilon_sqr[3][3]; - double epsilon_matrix[3][3]; - - double Jdet; - double viscosity; - double epsilon_eff; - double phi; - double t_pmp; - double scalar; - double scalar_def; - double scalar_ocean; - double scalar_transient; - - /*Collapsed formulation: */ - Tria* tria=NULL; - - /*parameters: */ - double dt; - - /*inputs: */ - bool onwater; - bool onbed; - bool shelf; - Input* vx_input=NULL; - Input* vy_input=NULL; - Input* vz_input=NULL; - - /*retrieve inputs :*/ - inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); - - /*retrieve some parameters: */ - this->parameters->FindParam(&dt,DtEnum); - - /*If on water, skip: */ - if(onwater)return; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /*recovre material parameters: */ - rho_water=matpar->GetRhoWater(); - rho_ice=matpar->GetRhoIce(); - gravity=matpar->GetG(); - heatcapacity=matpar->GetHeatCapacity(); - beta=matpar->GetBeta(); - meltingpoint=matpar->GetMeltingPoint(); - - /* 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.: */ - - /*Get gaussian points: */ - GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss); - - /*Retrieve all inputs we will be needing: */ - vx_input=inputs->GetInput(VxEnum); - vy_input=inputs->GetInput(VyEnum); - vz_input=inputs->GetInput(VzEnum); - - /* Start looping on the number of gaussian points: */ - for (igarea=0; igareaGetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input); - matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); - - /* Get nodal functions */ - GetNodalFunctions(&L[0], gauss_coord); - - /*Build deformational heating: */ - GetPhi(&phi, &epsilon[0], viscosity); - - /*Build pe_gaussian */ - scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss_weight; - if(dt){ - scalar_def=scalar_def*dt; - } - - for(i=0;iGetParameterValue(&temperature, gauss_coord,TemperatureEnum); - scalar_transient=temperature*Jdet*gauss_weight; - for(i=0;iCreatePVectorThermalShelf(pg); - delete tria; - } - - /* Geothermal flux on ice sheet base and basal friction */ - if(onbed && !shelf){ - - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVectorThermalSheet(pg); - delete tria; - } - extern int my_rank; - - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&vert_gauss_coord); - xfree((void**)&area_gauss_weights); - xfree((void**)&vert_gauss_weights); - -} -/*}}}*/ -/*FUNCTION Penta::VecExtrude {{{1*/ -void Penta::VecExtrude(Vec vector,double* vector_serial,int iscollapsed){ - - /* node data: */ - int i; - Node* node=NULL; - int extrude=0; - - /*inputs: */ - bool collapse; - bool onbed; - - /*retrieve inputs :*/ - inputs->GetParameterValue(&collapse,CollapseEnum); - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - - /*Figure out if we should extrude for this element: */ - if (iscollapsed){ - /*From higher level, we are told to extrude only elements that have the collapse flag on: */ - if (collapse)extrude=1; - else extrude=0; - } - else{ - /*From higher level, we are told to extrude all elements: */ - extrude=1; - } - - /*Now, extrusion starts from the bed on, so double check this element is on - * the bedrock: */ - if(onbed==0)extrude=0; - - /*Go on and extrude vector: */ - if (extrude){ - - /* node data: */ - int dof1; - double vectorel; - - /*this penta is a collapsed macayeal. For each node on the base of this penta, - * we grab the vector. Once we know the vector, we follow the upper nodes, - * inserting the same vector value into vector, until we reach the surface: */ - for(i=0;i<3;i++){ - - node=nodes[i]; //base nodes - dof1=node->GetDofList1(); - - /*get vector for this base node: */ - vectorel=vector_serial[dof1]; - - //go throvectorn all nodes which sit on top of this node, until we reach the surface, - //and plvector vector in vector - for(;;){ - - dof1=node->GetDofList1(); - VecSetValues(vector,1,&dof1,&vectorel,INSERT_VALUES); - - if (node->IsOnSurface())break; - /*get next node: */ - node=node->GetUpperNode(); - } - } - } -} -/*}}}*/ -/*FUNCTION Penta::InputExtrude {{{1*/ -void Penta::InputExtrude(int enum_type,bool only_if_collapsed){ - - bool onbed,collapse=false; - Penta *penta = NULL; - Input *original_input = NULL; - - /*recover parameters: */ - if (only_if_collapsed) inputs->GetParameterValue(&collapse,CollapseEnum); - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - - /*First: if only_if_collapsed, check wether this penta is collapsed*/ - if (only_if_collapsed && !collapse) return; - - /*Are we on the base, not on the surface?:*/ - if(onbed){ - - /*OK, we are on bed. we will follow the steps: - * 1: find input and extrude it. - * 2: follow the upper element until we reach the surface - * 3: for each element, we will add a copy of the extruded input*/ - - /*Step1: Extrude the original input: */ - original_input=(Input*)this->inputs->GetInput(enum_type); - if(!original_input) ISSMERROR("%s%s"," could not find input with enum:",EnumAsString(enum_type)); - original_input->Extrude(); - - /*Stop if there is only one layer of element*/ - if (this->IsOnSurface()) return; - - /*Step 2: this input has been extruded for this element, now follow the upper element*/ - penta=this; - for(;;){ - - /* get upper Penta*/ - penta=penta->GetUpperElement(); - ISSMASSERT(penta->Id()!=this->id); - - /*Add input of the basal element to penta->inputs*/ - Input* copy=NULL; - copy=(Input*)original_input->copy(); - penta->inputs->AddInput((Input*)copy); - - /*Stop if we have reached the surface*/ - if (penta->IsOnSurface()) break; - - } - } - - return; -} -/*}}}*/ -/*FUNCTION Penta::GetB {{{1*/ -void Penta::GetB(double* B, double* xyz_list, double* gauss_coord){ - - /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. - * For grid i, Bi can be expressed in the actual coordinate system - * by: - * Bi=[ dh/dx 0 ] - * [ 0 dh/dy ] - * [ 1/2*dh/dy 1/2*dh/dx ] - * [ 1/2*dh/dz 0 ] - * [ 0 1/2*dh/dz ] - * where h is the interpolation function for grid i. - * - * We assume B has been allocated already, of size: 5x(NDOF2*numgrids) - */ - - int i; - const int numgrids=6; - const int NDOF3=3; - const int NDOF2=2; - - double dh1dh6[NDOF3][numgrids]; - - /*Get dh1dh6 in actual coordinate system: */ - GetNodalFunctionsDerivatives(&dh1dh6[0][0],xyz_list, gauss_coord); - - /*Build B: */ - for (i=0;iGetDofList(&doflist_per_node[0],&numberofdofspernode); - for(j=0;jGetDofList1(); - } - -} -/*}}}*/ -/*FUNCTION Penta::GetJacobian {{{1*/ -void Penta::GetJacobian(double* J, double* xyz_list,double* gauss_coord){ - - const int NDOF3=3; - int i,j; - - /*The Jacobian is constant over the element, discard the gaussian points. - * J is assumed to have been allocated of size NDOF2xNDOF2.*/ - - double A1,A2,A3; //area coordinates - double xi,eta,zi; //parametric coordinates - - double x1,x2,x3,x4,x5,x6; - double y1,y2,y3,y4,y5,y6; - double z1,z2,z3,z4,z5,z6; - - /*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */ - A1=gauss_coord[0]; - A2=gauss_coord[1]; - A3=gauss_coord[2]; - - xi=A2-A1; - eta=SQRT3*A3; - zi=gauss_coord[3]; - - x1=*(xyz_list+3*0+0); - x2=*(xyz_list+3*1+0); - x3=*(xyz_list+3*2+0); - x4=*(xyz_list+3*3+0); - x5=*(xyz_list+3*4+0); - x6=*(xyz_list+3*5+0); - - y1=*(xyz_list+3*0+1); - y2=*(xyz_list+3*1+1); - y3=*(xyz_list+3*2+1); - y4=*(xyz_list+3*3+1); - y5=*(xyz_list+3*4+1); - y6=*(xyz_list+3*5+1); - - z1=*(xyz_list+3*0+2); - z2=*(xyz_list+3*1+2); - z3=*(xyz_list+3*2+2); - z4=*(xyz_list+3*3+2); - z5=*(xyz_list+3*4+2); - z6=*(xyz_list+3*5+2); - - - *(J+NDOF3*0+0)=0.25*(x1-x2-x4+x5)*zi+0.25*(-x1+x2-x4+x5); - *(J+NDOF3*1+0)=SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*zi+SQRT3/12.0*(-x1-x2+2*x3-x4-x5+2*x6); - *(J+NDOF3*2+0)=SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*eta+1/4*(x1-x2-x4+x5)*xi +0.25*(-x1+x5-x2+x4); - - *(J+NDOF3*0+1)=0.25*(y1-y2-y4+y5)*zi+0.25*(-y1+y2-y4+y5); - *(J+NDOF3*1+1)=SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*zi+SQRT3/12.0*(-y1-y2+2*y3-y4-y5+2*y6); - *(J+NDOF3*2+1)=SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*eta+0.25*(y1-y2-y4+y5)*xi+0.25*(y4-y1+y5-y2); - - *(J+NDOF3*0+2)=0.25*(z1-z2-z4+z5)*zi+0.25*(-z1+z2-z4+z5); - *(J+NDOF3*1+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*zi+SQRT3/12.0*(-z1-z2+2*z3-z4-z5+2*z6); - *(J+NDOF3*2+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*eta+0.25*(z1-z2-z4+z5)*xi+0.25*(-z1+z5-z2+z4); - -} -/*}}}*/ -/*FUNCTION Penta::GetJacobianDeterminant {{{1*/ -void Penta::GetJacobianDeterminant(double* Jdet, double* xyz_list,double* gauss_coord){ - - /*On a penta, Jacobian varies according to coordinates. We need to get the Jacobian, and take - * the determinant of it: */ - const int NDOF3=3; - - double J[NDOF3][NDOF3]; - - GetJacobian(&J[0][0],xyz_list,gauss_coord); - - *Jdet= J[0][0]*J[1][1]*J[2][2]-J[0][0]*J[1][2]*J[2][1]-J[1][0]*J[0][1]*J[2][2]+J[1][0]*J[0][2]*J[2][1]+J[2][0]*J[0][1]*J[1][2]-J[2][0]*J[0][2]*J[1][1]; - - if(*Jdet<0){ - ISSMERROR("%s%i","negative jacobian determinant on element ",id); - } -} -/*}}}*/ -/*FUNCTION Penta::GetJacobianInvert {{{1*/ -void Penta::GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss_coord){ - - double Jdet; - const int NDOF3=3; - - /*Call Jacobian routine to get the jacobian:*/ - GetJacobian(Jinv, xyz_list, gauss_coord); - - /*Invert Jacobian matrix: */ - MatrixInverse(Jinv,NDOF3,NDOF3,NULL,0,&Jdet); -} -/*}}}*/ -/*FUNCTION Penta::GetLStokes {{{1*/ -void Penta::GetLStokes(double* LStokes, double* gauss_coord_tria){ - - /* - * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. - * For grid i, Li can be expressed in the actual coordinate system - * by: - * Li=[ h 0 0 0] - * [ 0 h 0 0] - * [ 0 0 h 0] - * [ 0 0 h 0] - * [ h 0 0 0] - * [ 0 h 0 0] - * [ h 0 0 0] - * [ 0 h 0 0] - * [ 0 0 h 0] - * [ 0 0 h 0] - * [ 0 0 h 0] - * [ h 0 0 0] - * [ 0 h 0 0] - * [ 0 0 h 0] - * where h is the interpolation function for grid i. - */ - - int i; - const int numgrids2d=3; - int num_dof=4; - - double l1l2l3[numgrids2d]; - - - /*Get l1l2l3 in actual coordinate system: */ - l1l2l3[0]=gauss_coord_tria[0]; - l1l2l3[1]=gauss_coord_tria[1]; - l1l2l3[2]=gauss_coord_tria[2]; - - /*Build LStokes: */ - for (i=0;i<3;i++){ - *(LStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i]; - *(LStokes+num_dof*numgrids2d*0+num_dof*i+1)=0; - *(LStokes+num_dof*numgrids2d*0+num_dof*i+2)=0; - *(LStokes+num_dof*numgrids2d*0+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*1+num_dof*i)=0; - *(LStokes+num_dof*numgrids2d*1+num_dof*i+1)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*1+num_dof*i+2)=0; - *(LStokes+num_dof*numgrids2d*1+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*2+num_dof*i)=0; - *(LStokes+num_dof*numgrids2d*2+num_dof*i+1)=0; - *(LStokes+num_dof*numgrids2d*2+num_dof*i+2)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*2+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*3+num_dof*i)=0; - *(LStokes+num_dof*numgrids2d*3+num_dof*i+1)=0; - *(LStokes+num_dof*numgrids2d*3+num_dof*i+2)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*3+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*4+num_dof*i)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*4+num_dof*i+1)=0; - *(LStokes+num_dof*numgrids2d*4+num_dof*i+2)=0; - *(LStokes+num_dof*numgrids2d*4+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*5+num_dof*i)=0; - *(LStokes+num_dof*numgrids2d*5+num_dof*i+1)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*5+num_dof*i+2)=0; - *(LStokes+num_dof*numgrids2d*5+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*6+num_dof*i)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*6+num_dof*i+1)=0; - *(LStokes+num_dof*numgrids2d*6+num_dof*i+2)=0; - *(LStokes+num_dof*numgrids2d*6+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*7+num_dof*i)=0; - *(LStokes+num_dof*numgrids2d*7+num_dof*i+1)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*7+num_dof*i+2)=0; - *(LStokes+num_dof*numgrids2d*7+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*8+num_dof*i)=0; - *(LStokes+num_dof*numgrids2d*8+num_dof*i+1)=0; - *(LStokes+num_dof*numgrids2d*8+num_dof*i+2)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*8+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*9+num_dof*i)=0; - *(LStokes+num_dof*numgrids2d*9+num_dof*i+1)=0; - *(LStokes+num_dof*numgrids2d*9+num_dof*i+2)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*9+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*10+num_dof*i)=0; - *(LStokes+num_dof*numgrids2d*10+num_dof*i+1)=0; - *(LStokes+num_dof*numgrids2d*10+num_dof*i+2)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*10+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*11+num_dof*i)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*11+num_dof*i+1)=0; - *(LStokes+num_dof*numgrids2d*11+num_dof*i+2)=0; - *(LStokes+num_dof*numgrids2d*11+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*12+num_dof*i)=0; - *(LStokes+num_dof*numgrids2d*12+num_dof*i+1)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*12+num_dof*i+2)=0; - *(LStokes+num_dof*numgrids2d*12+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*13+num_dof*i)=0; - *(LStokes+num_dof*numgrids2d*13+num_dof*i+1)=0; - *(LStokes+num_dof*numgrids2d*13+num_dof*i+2)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*13+num_dof*i+3)=0; - - } -} -/*}}}*/ -/*FUNCTION Penta::GetLprimeStokes {{{1*/ -void Penta::GetLprimeStokes(double* LprimeStokes, double* xyz_list, double* gauss_coord_tria, double* gauss_coord){ - - /* - * Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. - * For grid i, Lpi can be expressed in the actual coordinate system - * by: - * Lpi=[ h 0 0 0] - * [ 0 h 0 0] - * [ h 0 0 0] - * [ 0 h 0 0] - * [ 0 0 h 0] - * [ 0 0 h 0] - * [ 0 0 dh/dz 0] - * [ 0 0 dh/dz 0] - * [ 0 0 dh/dz 0] - * [dh/dz 0 dh/dx 0] - * [ 0 dh/dz dh/dy 0] - * [ 0 0 0 h] - * [ 0 0 0 h] - * [ 0 0 0 h] - * where h is the interpolation function for grid i. - */ - - int i; - const int numgrids2d=3; - int num_dof=4; - - double l1l2l3[numgrids2d]; - double dh1dh6[3][6]; - - - /*Get l1l2l3 in actual coordinate system: */ - l1l2l3[0]=gauss_coord_tria[0]; - l1l2l3[1]=gauss_coord_tria[1]; - l1l2l3[2]=gauss_coord_tria[2]; - - GetNodalFunctionsDerivatives(&dh1dh6[0][0],xyz_list,gauss_coord); - - /*Build LprimeStokes: */ - for (i=0;i<3;i++){ - *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i]; - *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+1)=0; - *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+2)=0; - *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+3)=0; - *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i)=0; - *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+1)=l1l2l3[i]; - *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+2)=0; - *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+3)=0; - *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i)=l1l2l3[i]; - *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+1)=0; - *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+2)=0; - *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+3)=0; - *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i)=0; - *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+1)=l1l2l3[i]; - *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+2)=0; - *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+3)=0; - *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i)=0; - *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+1)=0; - *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+2)=l1l2l3[i]; - *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+3)=0; - *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i)=0; - *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+1)=0; - *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+2)=l1l2l3[i]; - *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+3)=0; - *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i)=0; - *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+1)=0; - *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+2)=dh1dh6[2][i]; - *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+3)=0; - *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i)=0; - *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+1)=0; - *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+2)=dh1dh6[2][i]; - *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+3)=0; - *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i)=0; - *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+1)=0; - *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+2)=dh1dh6[2][i]; - *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+3)=0; - *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i)=dh1dh6[2][i]; - *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+1)=0; - *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+2)=dh1dh6[0][i]; - *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+3)=0; - *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i)=0; - *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+1)=dh1dh6[2][i]; - *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+2)=dh1dh6[1][i]; - *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+3)=0; - *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i)=0; - *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+1)=0; - *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+2)=0; - *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+3)=l1l2l3[i]; - *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i)=0; - *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+1)=0; - *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+2)=0; - *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+3)=l1l2l3[i]; - *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i)=0; - *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+1)=0; - *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+2)=0; - *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+3)=l1l2l3[i]; - - } -} -/*}}}*/ -/*FUNCTION Penta::GetMatrixInvert {{{1*/ -void Penta::GetMatrixInvert(double* Ke_invert, double* Ke){ - /*Inverse a 3 by 3 matrix only */ - - double a,b,c,d,e,f,g,h,i; - double det; - int calculationdof=3; - - /*Take the matrix components: */ - a=*(Ke+calculationdof*0+0); - b=*(Ke+calculationdof*0+1); - c=*(Ke+calculationdof*0+2); - d=*(Ke+calculationdof*1+0); - e=*(Ke+calculationdof*1+1); - f=*(Ke+calculationdof*1+2); - g=*(Ke+calculationdof*2+0); - h=*(Ke+calculationdof*2+1); - i=*(Ke+calculationdof*2+2); - - det=a*(e*i-f*h)-b*(d*i-f*g)+c*(d*h-e*g); - - *(Ke_invert+calculationdof*0+0)=(e*i-f*h)/det; - *(Ke_invert+calculationdof*0+1)=(c*h-b*i)/det; - *(Ke_invert+calculationdof*0+2)=(b*f-c*e)/det; - *(Ke_invert+calculationdof*1+0)=(f*g-d*i)/det; - *(Ke_invert+calculationdof*1+1)=(a*i-c*g)/det; - *(Ke_invert+calculationdof*1+2)=(c*d-a*f)/det; - *(Ke_invert+calculationdof*2+0)=(d*h-e*g)/det; - *(Ke_invert+calculationdof*2+1)=(b*g-a*h)/det; - *(Ke_invert+calculationdof*2+2)=(a*e-b*d)/det; - -} -/*}}}*/ -/*FUNCTION Penta::GetNodalFunctions {{{1*/ -void Penta::GetNodalFunctions(double* l1l6, double* gauss_coord){ - - /*This routine returns the values of the nodal functions at the gaussian point.*/ - - l1l6[0]=gauss_coord[0]*(1-gauss_coord[3])/2.0; - - l1l6[1]=gauss_coord[1]*(1-gauss_coord[3])/2.0; - - l1l6[2]=gauss_coord[2]*(1-gauss_coord[3])/2.0; - - l1l6[3]=gauss_coord[0]*(1+gauss_coord[3])/2.0; - - l1l6[4]=gauss_coord[1]*(1+gauss_coord[3])/2.0; - - l1l6[5]=gauss_coord[2]*(1+gauss_coord[3])/2.0; - -} -/*}}}*/ -/*FUNCTION Penta::GetNodalFunctionsDerivatives {{{1*/ -void Penta::GetNodalFunctionsDerivatives(double* dh1dh6,double* xyz_list, double* gauss_coord){ - - /*This routine returns the values of the nodal functions derivatives (with respect to the actual coordinate system: */ - - - int i; - const int NDOF3=3; - const int numgrids=6; - - double dh1dh6_ref[NDOF3][numgrids]; - double Jinv[NDOF3][NDOF3]; - - /*Get derivative values with respect to parametric coordinate system: */ - GetNodalFunctionsDerivativesReference(&dh1dh6_ref[0][0], gauss_coord); - - /*Get Jacobian invert: */ - GetJacobianInvert(&Jinv[0][0], xyz_list, gauss_coord); - - /*Build dh1dh3: - * - * [dhi/dx]= Jinv*[dhi/dr] - * [dhi/dy] [dhi/ds] - * [dhi/dz] [dhi/dn] - */ - - for (i=0;iGetParameterValue(pvalue,&gauss[grid][0],enumtype); - return; - -} -/*}}}*/ -/*FUNCTION Penta::GetPhi {{{1*/ -void Penta::GetPhi(double* phi, double* epsilon, double viscosity){ - /*Compute deformational heating from epsilon and viscosity */ - - double epsilon_matrix[3][3]; - double epsilon_eff; - double epsilon_sqr[3][3]; - - /* Build epsilon matrix */ - epsilon_matrix[0][0]=*(epsilon+0); - epsilon_matrix[1][0]=*(epsilon+3); - epsilon_matrix[2][0]=*(epsilon+4); - epsilon_matrix[0][1]=*(epsilon+3); - epsilon_matrix[1][1]=*(epsilon+1); - epsilon_matrix[2][1]=*(epsilon+5); - epsilon_matrix[0][2]=*(epsilon+4); - epsilon_matrix[1][2]=*(epsilon+5); - epsilon_matrix[2][2]=*(epsilon+2); - - /* Effective value of epsilon_matrix */ - epsilon_sqr[0][0]=pow(epsilon_matrix[0][0],2); - epsilon_sqr[1][0]=pow(epsilon_matrix[1][0],2); - epsilon_sqr[2][0]=pow(epsilon_matrix[2][0],2); - epsilon_sqr[0][1]=pow(epsilon_matrix[0][1],2); - epsilon_sqr[1][1]=pow(epsilon_matrix[1][1],2); - epsilon_sqr[2][1]=pow(epsilon_matrix[2][1],2); - epsilon_sqr[0][2]=pow(epsilon_matrix[0][2],2); - epsilon_sqr[1][2]=pow(epsilon_matrix[1][2],2); - epsilon_sqr[2][2]=pow(epsilon_matrix[2][2],2); - - epsilon_eff=1/pow(2,0.5)*pow((epsilon_sqr[0][0]+epsilon_sqr[0][1]+ epsilon_sqr[0][2]+ epsilon_sqr[1][0]+ epsilon_sqr[1][1]+ epsilon_sqr[1][2]+ epsilon_sqr[2][0]+ epsilon_sqr[2][1]+ epsilon_sqr[2][2]),0.5); - *phi=2*pow(epsilon_eff,2.0)*viscosity; -} -/*}}}*/ -/*FUNCTION Penta::GetStrainRate3dPattyn{{{1*/ -void Penta::GetStrainRate3dPattyn(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input){ - /*Compute the 3d Blatter/PattynStrain Rate (5 components): - * - * epsilon=[exx eyy exy exz eyz] - * - * with exz=1/2 du/dz - * eyz=1/2 dv/dz - * - * the contribution of vz is neglected - */ - - int i; - - double epsilonvx[5]; - double epsilonvy[5]; - - /*Check that both inputs have been found*/ - if (!vx_input || !vy_input){ - ISSMERROR("Input missing. Here are the input pointers we have for vx: %p, vy: %p\n",vx_input,vy_input); - } - - /*Get strain rate assuming that epsilon has been allocated*/ - vx_input->GetVxStrainRate3dPattyn(epsilonvx,xyz_list,gauss); - vy_input->GetVyStrainRate3dPattyn(epsilonvy,xyz_list,gauss); - - /*Sum all contributions*/ - for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]; - -} -/*}}}*/ -/*FUNCTION Penta::GetStrainRate3d{{{1*/ -void Penta::GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input, Input* vz_input){ - /*Compute the 3d Strain Rate (6 components): - * - * epsilon=[exx eyy ezz exy exz eyz] - */ - - int i; - - double epsilonvx[6]; - double epsilonvy[6]; - double epsilonvz[6]; - - /*Check that both inputs have been found*/ - if (!vx_input || !vy_input || !vz_input){ - ISSMERROR("Input missing. Here are the input pointers we have for vx: %p, vy: %p, vz: %p\n",vx_input,vy_input,vz_input); - } - - /*Get strain rate assuming that epsilon has been allocated*/ - vx_input->GetVxStrainRate3d(epsilonvx,xyz_list,gauss); - vy_input->GetVyStrainRate3d(epsilonvy,xyz_list,gauss); - vz_input->GetVzStrainRate3d(epsilonvz,xyz_list,gauss); - - /*Sum all contributions*/ - for(i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i]; - -} -/*}}}*/ -/*FUNCTION Penta::GradjB {{{1*/ -void Penta::GradjB(Vec gradient){ - - int i; - Tria* tria=NULL; - TriaVertexInput* triavertexinput=NULL; - - /*inputs: */ - bool onwater; - bool collapse; - bool onbed; - - /*retrieve inputs :*/ - inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - inputs->GetParameterValue(&collapse,CollapseEnum); - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - - /*If on water, skip: */ - if(onwater)return; - - if (collapse){ - /*Bail out element if collapsed (2d) and not on bed*/ - if (!onbed) return; - - /*This element should be collapsed into a tria element at its base. Create this tria element, - * and compute gardj*/ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - tria->GradjB(gradient); - delete tria; - } - else{ - /*B is a 2d field, use MacAyeal(2d) gradient even if it is Stokes or Pattyn*/ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - tria->GradjB(gradient); - delete tria; - } - - -} -/*}}}*/ /*FUNCTION Penta::ReduceMatrixStokes {{{1*/ void Penta::ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp){ @@ -6305,4 +5248,101 @@ } /*}}}1*/ +/*FUNCTION Penta::SpawnTria {{{1*/ +Tria* Penta::SpawnTria(int g0, int g1, int g2){ + + int i; + int analysis_counter; + + /*go into parameters and get the analysis_counter: */ + this->parameters->FindParam(&analysis_counter,AnalysisCounterEnum); + + /*out of grids g0,g1 and g2 from Penta, build a tria element: */ + Tria* tria=NULL; + int indices[3]; + int zero=0; + Parameters* tria_parameters=NULL; + Inputs* tria_inputs=NULL; + Results* tria_results=NULL; + + indices[0]=g0; + indices[1]=g1; + indices[2]=g2; + + tria_parameters=this->parameters; + tria_inputs=(Inputs*)this->inputs->SpawnTriaInputs(indices); + tria_results=(Results*)this->results->SpawnTriaResults(indices); + + tria=new Tria(); + tria->id=this->id; + tria->inputs=tria_inputs; + tria->results=tria_results; + tria->parameters=tria_parameters; + tria->element_type=P1Enum; //Only P1 CG for now + this->SpawnTriaHook(dynamic_cast(tria),&indices[0]); + + /*recover nodes, matice and matpar: */ + tria->nodes=(Node**)tria->hnodes[analysis_counter]->deliverp(); + tria->matice=(Matice*)tria->hmatice->delivers(); + tria->matpar=(Matpar*)tria->hmatpar->delivers(); + + return tria; +} +/*}}}*/ +/*FUNCTION Penta::SpawnBeam {{{1*/ +void* Penta::SpawnBeam(int g0, int g1){ + + int i; + + /*out of grids g0,g1 and g2 from Penta, build a beam element: */ + Beam* beam=NULL; + int indices[2]; + int zero=0; + Parameters *beam_parameters = NULL; + Inputs *beam_inputs = NULL; + + indices[0]=g0; + indices[1]=g1; + + beam_parameters=this->parameters; + beam_inputs=(Inputs*)this->inputs->SpawnBeamInputs(indices); + + beam=new Beam(); + beam->id=this->id; + beam->inputs=beam_inputs; + beam->parameters=beam_parameters; + + /*now deal with ndoes,matice and matpar: */ + beam->nodes=(Node**)xmalloc(2*sizeof(Node*)); + for(i=0;i<2;i++)beam->nodes[i]=this->nodes[indices[i]]; + beam->matice=this->matice; + beam->matpar=this->matpar; + + return beam; +} +/*}}}*/ +/*FUNCTION Penta::SpawnSing {{{1*/ +void* Penta::SpawnSing(int index){ + + Sing* sing=NULL; + int zero=0; + Parameters *sing_parameters = NULL; + Inputs *sing_inputs = NULL; + + sing_parameters=this->parameters; + sing_inputs=(Inputs*)this->inputs->SpawnSingInputs(index); + + sing=new Sing(); + sing->id=this->id; + sing->inputs=sing_inputs; + sing->parameters=sing_parameters; + + /*now deal with nodes,matice and matpar: */ + sing->node=this->nodes[index]; + sing->matice=this->matice; + sing->matpar=this->matpar; + + return sing; +} +/*}}}*/ /*FUNCTION Penta::SurfaceNormal {{{1*/ void Penta::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){ @@ -6331,2 +5371,66 @@ } /*}}}*/ +/*FUNCTION Penta::VecExtrude {{{1*/ +void Penta::VecExtrude(Vec vector,double* vector_serial,int iscollapsed){ + + /* node data: */ + int i; + Node* node=NULL; + int extrude=0; + + /*inputs: */ + bool collapse; + bool onbed; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&collapse,CollapseEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + + /*Figure out if we should extrude for this element: */ + if (iscollapsed){ + /*From higher level, we are told to extrude only elements that have the collapse flag on: */ + if (collapse)extrude=1; + else extrude=0; + } + else{ + /*From higher level, we are told to extrude all elements: */ + extrude=1; + } + + /*Now, extrusion starts from the bed on, so double check this element is on + * the bedrock: */ + if(onbed==0)extrude=0; + + /*Go on and extrude vector: */ + if (extrude){ + + /* node data: */ + int dof1; + double vectorel; + + /*this penta is a collapsed macayeal. For each node on the base of this penta, + * we grab the vector. Once we know the vector, we follow the upper nodes, + * inserting the same vector value into vector, until we reach the surface: */ + for(i=0;i<3;i++){ + + node=nodes[i]; //base nodes + dof1=node->GetDofList1(); + + /*get vector for this base node: */ + vectorel=vector_serial[dof1]; + + //go throvectorn all nodes which sit on top of this node, until we reach the surface, + //and plvector vector in vector + for(;;){ + + dof1=node->GetDofList1(); + VecSetValues(vector,1,&dof1,&vectorel,INSERT_VALUES); + + if (node->IsOnSurface())break; + /*get next node: */ + node=node->GetUpperNode(); + } + } + } +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 4884) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 4885) @@ -136,22 +136,8 @@ void CreatePVectorSlope( Vec pg); void CreatePVectorThermal( Vec pg); - double* GaussFromNode(Node* node); - void GetB(double* pB, double* xyz_list, double* gauss_coord); - void GetBPrime(double* B, double* xyz_list, double* gauss_coord); - void GetBPrime_vert(double* B, double* xyz_list, double* gauss_coord); - void GetBStokes(double* B, double* xyz_list, double* gauss_coord); - void GetB_advec(double* B_advec, double* xyz_list, double* gauss_coord); - void GetB_artdiff(double* B_artdiff, double* xyz_list, double* gauss_coord); - void GetB_conduct(double* B_conduct, double* xyz_list, double* gauss_coord); - void GetB_vert(double* B, double* xyz_list, double* gauss_coord); - void GetBprimeStokes(double* B_prime, double* xyz_list, double* gauss_coord); - void GetBprime_advec(double* Bprime_advec, double* xyz_list, double* gauss_coord); + double* GaussFromNode(Node* node); void GetDofList(int* doflist,int* pnumberofdofs); void GetDofList1(int* doflist); - void GetJacobian(double* J, double* xyz_list,double* gauss_coord); - void GetJacobianDeterminant(double* Jdet, double* xyz_list,double* gauss_coord); - void GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss_coord); - void GetLStokes(double* LStokes, double* gauss_coord_tria); - void GetLprimeStokes(double* LprimeStokes, double* xyz_list, double* gauss_coord_tria, double* gauss_coord); + int GetElementType(void); void GetMatrixInvert(double* Ke_invert, double* Ke); void GetNodalFunctions(double* l1l6, double* gauss_coord); @@ -161,5 +147,4 @@ void GetNodalFunctionsDerivativesStokes(double* dh1dh7,double* xyz_list, double* gauss_coord); void GetNodalFunctionsStokes(double* l1l7, double* gauss_coord); - void GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_coord); void GetParameterValue(double* pvalue, double* v_list,double* gauss_coord); void GetParameterValue(double* pvalue,Node* node,int enumtype); Index: /issm/trunk/src/c/objects/Elements/PentaRef.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/PentaRef.cpp (revision 4884) +++ /issm/trunk/src/c/objects/Elements/PentaRef.cpp (revision 4885) @@ -52,2 +52,859 @@ /*Reference Element numerics*/ +/*FUNCTION PentaRef::GetBPattyn {{{1*/ +void PentaRef::GetBPattyn(double* B, double* xyz_list, double* gauss){ + /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. + * For grid i, Bi can be expressed in the actual coordinate system + * by: + * Bi=[ dh/dx 0 ] + * [ 0 dh/dy ] + * [ 1/2*dh/dy 1/2*dh/dx ] + * [ 1/2*dh/dz 0 ] + * [ 0 1/2*dh/dz ] + * where h is the interpolation function for grid i. + * + * We assume B has been allocated already, of size: 5x(NDOF2*numgrids) + */ + + int i; + const int numgrids=6; + const int NDOF3=3; + const int NDOF2=2; + + double dh1dh6[NDOF3][numgrids]; + + /*Get dh1dh6 in actual coordinate system: */ + GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss); + + /*Build B: */ + for (i=0;ielement_type_list=(int*)xmalloc(this->numanalyses*sizeof(int)); - for(i=0;i /* DBL_EPSILON */ #include "../../include/include.h" - +/*}}}*/ + +/*FUNCTION TripleMultiply {{{1*/ int TripleMultiply( double* a, int nrowa, int ncola, int itrna, double* b, int nrowb, int ncolb, int itrnb, double* c, int nrowc, int ncolc, int itrnc, double* d, int iaddd){ + /*TripleMultiply Perform triple matrix product a*b*c+d.*/ int idima,idimb,idimc,idimd; @@ -116,44 +95,8 @@ xfree((void **)&dtemp); } -} - - - - - -/* - - MatrixMultiply Perform matrix multiplication a*b+c. - - - Input: - double* a first matrix - int nrowa number of rows in a - int ncola number of columns in a - int itrna transpose flag for a - double* b second matrix - int nrowb number of rows in b - int ncolb number of columns in b - int itrnb transpose flag for b - double* c addition matrix - int iaddc addition flag for c - - Output: - double* c product matrix - int function return value - -*/ - -int MatrixMultiply( double* a, - int nrowa, - int ncola, - int itrna, - double* b, - int nrowb, - int ncolb, - int itrnb, - double* c, - int iaddc ) -{ +}/*}}}*/ +/*FUNCTION MatrixMuliply {{{1*/ +int MatrixMultiply( double* a, int nrowa, int ncola, int itrna, double* b, int nrowb, int ncolb, int itrnb, double* c, int iaddc ){ + /*MatrixMultiply Perform matrix multiplication a*b+c.*/ int noerr=1; int i,j,k,ipta,iptb,iptc; @@ -161,5 +104,4 @@ /* set up dimensions and increments for matrix a */ - if (!itrna) { nrowc=nrowa; @@ -168,5 +110,4 @@ jinca=1; } - else { nrowc=ncola; @@ -177,5 +118,4 @@ /* set up dimensions and increments for matrix b */ - if (!itrnb) { ncolc=ncolb; @@ -184,5 +124,4 @@ jincb=1; } - else { ncolc=nrowb; @@ -225,24 +164,8 @@ return noerr; -} - - -/* - - MatrixInverse Perform matrix inversion and linear equation solution. - - - Input: - double* a coefficient matrix [ndim][ndim] - int ndim row and column dimension of a - int nrow number of rows in a - double* b right-hand vectors [ndim][nvec] (may be nULL) - int nvec number of vectors in b (0 for NULL b) - - Output: - double* a inverted matrix [ndim][ndim] - double* b solution vectors [ndim][nvec] (may be NULL) - double* pdet determinant (may be NULL) - int error flag +}/*}}}*/ +/*FUNCTION MatrixInverse {{{1*/ +int MatrixInverse( double* a, int ndim, int nrow, double* b, int nvec, double* pdet ){ +/* MatrixInverse Perform matrix inversion and linear equation solution. This function uses Gaussian elimination on the original matrix @@ -253,15 +176,5 @@ Gaussian elimination is inherently inefficient, and so this is - intended for small matrices. - -*/ - -int MatrixInverse( double* a, - int ndim, - int nrow, - double* b, - int nvec, - double* pdet ) -{ + intended for small matrices. */ int noerr=1; int i,j,k,ipt,jpt,irow,icol,ipiv,ncol; @@ -424,3 +337,3 @@ return noerr; -} +}/*}}}*/