Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 8653) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 8654) @@ -3874,4 +3874,579 @@ } /*}}}*/ +/*FUNCTION Penta::GetBasalElement{{{1*/ +Penta* Penta::GetBasalElement(void){ + + /*Output*/ + Penta* penta=NULL; + + /*Go through all elements till the bed is reached*/ + penta=this; + for(;;){ + /*Stop if we have reached the surface, else, take lower penta*/ + if (penta->IsOnBed()) break; + + /* get lower Penta*/ + penta=penta->GetLowerElement(); + _assert_(penta->Id()!=this->id); + } + + /*return output*/ + return penta; +} +/*}}}*/ +/*FUNCTION Penta::GetDofList {{{1*/ +void Penta::GetDofList(int** pdoflist,int approximation_enum,int setenum){ + + int i,j,count=0; + int numberofdofs=0; + int* doflist=NULL; + + /*First, figure out size of doflist: */ + for(i=0;i<6;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); + + /*Allocate: */ + doflist=(int*)xmalloc(numberofdofs*sizeof(int)); + + /*Populate: */ + count=0; + for(i=0;i<6;i++){ + nodes[i]->GetDofList(doflist+count,approximation_enum,setenum); + count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); + } + + /*Assign output pointers:*/ + *pdoflist=doflist; +} +/*}}}*/ +/*FUNCTION Penta::GetDofList1 {{{1*/ +void Penta::GetDofList1(int* doflist){ + + int i; + for(i=0;i<6;i++) doflist[i]=nodes[i]->GetDofList1(); + +} +/*}}}*/ +/*FUNCTION Penta::GetElementType {{{1*/ +int Penta::GetElementType(){ + + /*return PentaRef field*/ + return this->element_type; +} +/*}}}*/ +/*FUNCTION Penta::GetHorizontalNeighboorSids {{{1*/ +int* Penta::GetHorizontalNeighboorSids(){ + + /*return PentaRef field*/ + return &this->horizontalneighborsids[0]; + +} +/*}}}*/ +/*FUNCTION Penta::GetLowerElement{{{1*/ +Penta* Penta::GetLowerElement(void){ + + Penta* upper_penta=NULL; + + upper_penta=(Penta*)verticalneighbors[0]; //first one (0) under, second one (1) above + + return upper_penta; +} +/*}}}*/ +/*FUNCTION Penta::GetNodeIndex {{{1*/ +int Penta::GetNodeIndex(Node* node){ + + _assert_(nodes); + for(int i=0;iGetInput(enumtype); + if (!input) _error_("Input %s not found in element",EnumToStringx(enumtype)); + + /*Checks in debugging mode*/ + _assert_(pvalue); + + /* Start looping on the number of vertices: */ + gauss=new GaussPenta(); + for (int iv=0;ivGaussVertex(iv); + input->GetParameterValue(&pvalue[iv],gauss); + } + + /*clean-up*/ + delete gauss; +} +/*}}}*/ +/*FUNCTION Penta::GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue) {{{1*/ +void Penta::GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue){ + + /*Intermediaries*/ + double value[NUMVERTICES]; + GaussPenta *gauss = NULL; + + /*Recover input*/ + Input* input=inputs->GetInput(enumtype); + + /*Checks in debugging mode*/ + _assert_(pvalue); + + /* Start looping on the number of vertices: */ + if (input){ + gauss=new GaussPenta(); + for (int iv=0;ivGaussVertex(iv); + input->GetParameterValue(&pvalue[iv],gauss); + } + } + else{ + for (int iv=0;ivGetInput(enumtype); + if(!input) _error_("No input of type %s found in tria",EnumToStringx(enumtype)); + + GaussPenta* gauss=new GaussPenta(); + gauss->GaussVertex(this->GetNodeIndex(node)); + + input->GetParameterValue(pvalue,gauss); + delete gauss; +} +/*}}}*/ +/*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 = Tr(sigma * eps) + * = Tr(sigma'* eps) + * = 2 * eps_eff * sigma'_eff + * = 4 * mu * eps_eff ^2*/ + *phi=4*pow(epsilon_eff,2.0)*viscosity; +} +/*}}}*/ +/*FUNCTION Penta::GetSidList{{{1*/ +void Penta::GetSidList(int* sidlist){ + + int i; + for(i=0;iGetSidList(); + +} +/*}}}*/ +/*FUNCTION Penta::GetSolutionFromInputs{{{1*/ +void Penta::GetSolutionFromInputs(Vec solution){ + + int analysis_type; + + /*retrive parameters: */ + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */ + if (analysis_type==DiagnosticHorizAnalysisEnum){ + int approximation; + inputs->GetParameterValue(&approximation,ApproximationEnum); + if(approximation==StokesApproximationEnum || approximation==NoneApproximationEnum){ + GetSolutionFromInputsDiagnosticStokes(solution); + } + else if (approximation==MacAyealApproximationEnum || approximation==PattynApproximationEnum || approximation==HutterApproximationEnum){ + GetSolutionFromInputsDiagnosticHoriz(solution); + } + else if (approximation==MacAyealPattynApproximationEnum || approximation==PattynStokesApproximationEnum || approximation==MacAyealStokesApproximationEnum){ + return; //the elements around will create the solution + } + } + else if(analysis_type==DiagnosticHutterAnalysisEnum){ + GetSolutionFromInputsDiagnosticHutter(solution); + } + else if(analysis_type==DiagnosticVertAnalysisEnum){ + GetSolutionFromInputsDiagnosticVert(solution); + } + else if(analysis_type==ThermalAnalysisEnum){ + GetSolutionFromInputsThermal(solution); + } + else if(analysis_type==EnthalpyAnalysisEnum){ + GetSolutionFromInputsEnthalpy(solution); + } + else{ + _error_("analysis: %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type)); + } +} +/*}}}*/ +/*FUNCTION Penta::GetSolutionFromInputsDiagnosticHoriz{{{1*/ +void Penta::GetSolutionFromInputsDiagnosticHoriz(Vec solution){ + + const int numdof=NDOF2*NUMVERTICES; + + int i; + int approximation; + int* doflist=NULL; + double vx,vy; + double values[numdof]; + GaussPenta* gauss; + + /*Get approximation enum and dof list: */ + inputs->GetParameterValue(&approximation,ApproximationEnum); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + + /*If the element is a coupling, do nothing: every node is also on an other elements + * (as coupling is between MacAyeal and Pattyn) so the other element will take care of it*/ + GetDofList(&doflist,approximation,GsetEnum); + + /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ + /*P1 element only for now*/ + gauss=new GaussPenta(); + for(i=0;iGaussVertex(i); + vx_input->GetParameterValue(&vx,gauss); + vy_input->GetParameterValue(&vy,gauss); + values[i*NDOF2+0]=vx; + values[i*NDOF2+1]=vy; + } + + /*Add value to global vector*/ + VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); + + /*Free ressources:*/ + delete gauss; + xfree((void**)&doflist); +} +/*}}}*/ +/*FUNCTION Penta::GetSolutionFromInputsDiagnosticHutter{{{1*/ +void Penta::GetSolutionFromInputsDiagnosticHutter(Vec solution){ + + const int numdof=NDOF2*NUMVERTICES; + + int i; + int* doflist=NULL; + double vx,vy; + double values[numdof]; + GaussPenta* gauss=NULL; + + /*Get dof list: */ + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + + /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ + /*P1 element only for now*/ + gauss=new GaussPenta(); + for(i=0;iGaussVertex(i); + vx_input->GetParameterValue(&vx,gauss); + vy_input->GetParameterValue(&vy,gauss); + values[i*NDOF2+0]=vx; + values[i*NDOF2+1]=vy; + } + + /*Add value to global vector*/ + VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); + + /*Free ressources:*/ + delete gauss; + xfree((void**)&doflist); +} +/*}}}*/ +/*FUNCTION Penta::GetSolutionFromInputsDiagnosticVert{{{1*/ +void Penta::GetSolutionFromInputsDiagnosticVert(Vec solution){ + + const int numdof=NDOF1*NUMVERTICES; + + int i; + int* doflist=NULL; + double vz; + double values[numdof]; + GaussPenta* gauss=NULL; + + /*Get dof list: */ + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); + + /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ + /*P1 element only for now*/ + gauss=new GaussPenta(); + for(i=0;iGaussVertex(i); + vz_input->GetParameterValue(&vz,gauss); + values[i]=vz; + } + + /*Add value to global vector*/ + VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); + + /*Free ressources:*/ + delete gauss; + xfree((void**)&doflist); +} +/*}}}*/ +/*FUNCTION Penta::GetSolutionFromInputsDiagnosticStokes{{{1*/ +void Penta::GetSolutionFromInputsDiagnosticStokes(Vec solution){ + + const int numdof=NDOF4*NUMVERTICES; + + int i; + int* doflist=NULL; + double vx,vy,vz,p; + double stokesreconditioning; + double values[numdof]; + GaussPenta *gauss; + + /*Get dof list: */ + GetDofList(&doflist,StokesApproximationEnum,GsetEnum); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); + Input* p_input =inputs->GetInput(PressureEnum); _assert_(p_input); + + /*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*/ + gauss=new GaussPenta(); + for(i=0;iGaussVertex(i); + vx_input->GetParameterValue(&vx,gauss); + vy_input->GetParameterValue(&vy,gauss); + vz_input->GetParameterValue(&vz,gauss); + p_input ->GetParameterValue(&p ,gauss); + values[i*NDOF4+0]=vx; + values[i*NDOF4+1]=vy; + values[i*NDOF4+2]=vz; + values[i*NDOF4+3]=p/stokesreconditioning; + } + + /*Add value to global vector*/ + VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); + + /*Free ressources:*/ + delete gauss; + xfree((void**)&doflist); +} +/*}}}*/ +/*FUNCTION Penta::GetSolutionFromInputsThermal{{{1*/ +void Penta::GetSolutionFromInputsThermal(Vec solution){ + + const int numdof=NDOF1*NUMVERTICES; + + int i; + int* doflist=NULL; + double values[numdof]; + double temp; + GaussPenta *gauss=NULL; + + /*Get dof list: */ + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + Input* t_input=inputs->GetInput(TemperatureEnum); _assert_(t_input); + + gauss=new GaussPenta(); + for(i=0;iGaussVertex(i); + t_input->GetParameterValue(&temp,gauss); + values[i]=temp; + } + + /*Add value to global vector*/ + VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); + + /*Free ressources:*/ + delete gauss; + xfree((void**)&doflist); +} +/*}}}*/ +/*FUNCTION Penta::GetSolutionFromInputsEnthalpy{{{1*/ +void Penta::GetSolutionFromInputsEnthalpy(Vec solution){ + + const int numdof=NDOF1*NUMVERTICES; + + int i; + int* doflist=NULL; + double values[numdof]; + double enthalpy; + GaussPenta *gauss=NULL; + + /*Get dof list: */ + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + Input* h_input=inputs->GetInput(EnthalpyEnum); _assert_(h_input); + + gauss=new GaussPenta(); + for(i=0;iGaussVertex(i); + h_input->GetParameterValue(&enthalpy,gauss); + values[i]=enthalpy; + } + + /*Add value to global vector*/ + VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); + + /*Free ressources:*/ + delete gauss; + xfree((void**)&doflist); +} +/*}}}*/ +/*FUNCTION Penta::GetStabilizationParameter {{{1*/ +double Penta::GetStabilizationParameter(double u, double v, double w, double diameter, double rho_ice, double heatcapacity, double thermalconductivity){ + /*Compute stabilization parameter*/ + + double normu; + double tau_parameter; + + normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5); + if(normu*diameter/(3*2*thermalconductivity/(rho_ice*heatcapacity))<1){ + tau_parameter=pow(diameter,2)/(3*2*2*thermalconductivity/(rho_ice*heatcapacity)); + } + else tau_parameter=diameter/(2*normu); + + return tau_parameter; +} +/*}}}*/ +/*FUNCTION Penta::GetStrainRate3dPattyn{{{1*/ +void Penta::GetStrainRate3dPattyn(double* epsilon,double* xyz_list, GaussPenta* 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){ + _error_("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, GaussPenta* 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){ + _error_("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*)verticalneighbors[1]; //first one under, second one above + + return upper_penta; +} +/*}}}*/ +/*FUNCTION Penta::GetVectorFromInputs{{{1*/ +void Penta::GetVectorFromInputs(Vec vector,int input_enum){ + + int doflist1[NUMVERTICES]; + + /*Get out if this is not an element input*/ + if (!IsInput(input_enum)) return; + + /*Prepare index list*/ + this->GetDofList1(&doflist1[0]); + + /*Get input (either in element or material)*/ + Input* input=inputs->GetInput(input_enum); + if(!input) _error_("Input %s not found in element",EnumToStringx(input_enum)); + + /*We found the enum. Use its values to fill into the vector, using the vertices ids: */ + input->GetVectorFromInputs(vector,&doflist1[0]); +} +/*}}}*/ +/*FUNCTION Penta::GetZcoord {{{1*/ +double Penta::GetZcoord(GaussPenta* gauss){ + + int i; + double z; + double xyz_list[NUMVERTICES][3]; + double z_list[NUMVERTICES]; + + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + for(i=0;iGetParameterValue(&approximation,ApproximationEnum); - if(IsOnShelf() || !IsOnBed()) return; - switch(approximation){ case MacAyealApproximationEnum: - tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. - tria->GradjDragMacAyeal(gradient); - delete tria->matice; delete tria; + GradjDragMacAyeal(gradient); break; case PattynApproximationEnum: - tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. - tria->GradjDragMacAyeal(gradient); - delete tria->matice; delete tria; - //GradjDragPattyn(gradient); + GradjDragPattyn(gradient); break; case StokesApproximationEnum: - tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. - tria->GradjDragStokes(gradient); - delete tria->matice; delete tria; - //GradjDragStokes(gradient); + GradjDragStokes(gradient); break; case NoneApproximationEnum: @@ -3921,37 +4485,13 @@ inputs->GetParameterValue(&approximation,ApproximationEnum); switch(approximation){ - case MacAyealApproximationEnum: - - /*This element should be collapsed into a tria element at its base*/ - if (!IsOnBed()) return; - - /*Depth Average B*/ - this->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum); - - tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face). - tria->GradjBMacAyeal(gradient); - delete tria->matice; delete tria; - - /*delete B average*/ - this->matice->inputs->DeleteInput(RheologyBbarEnum); + GradjBbarMacAyeal(gradient); break; - case PattynApproximationEnum: case StokesApproximationEnum: - - /*Gradient is computed on bed only (Bbar)*/ - if (!IsOnBed()) return; - - /*Depth Average B*/ - this->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum); - - /*B is a 2d field, use MacAyeal(2d) gradient even if it is Stokes or Pattyn*/ - tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face). - tria->GradjBMacAyeal(gradient); - delete tria->matice; delete tria; - - /*delete B average*/ - this->matice->inputs->DeleteInput(RheologyBbarEnum); + case PattynApproximationEnum: + GradjBbarPattyn(gradient); break; - + case StokesApproximationEnum: + GradjBbarStokes(gradient); + break; case NoneApproximationEnum: /*Gradient is 0*/ @@ -3999,579 +4539,236 @@ } /*}}}*/ -/*FUNCTION Penta::GetBasalElement{{{1*/ -Penta* Penta::GetBasalElement(void){ - - /*Output*/ - Penta* penta=NULL; - - /*Go through all elements till the bed is reached*/ - penta=this; - for(;;){ - /*Stop if we have reached the surface, else, take lower penta*/ - if (penta->IsOnBed()) break; - - /* get lower Penta*/ - penta=penta->GetLowerElement(); - _assert_(penta->Id()!=this->id); - } - - /*return output*/ - return penta; -} -/*}}}*/ -/*FUNCTION Penta::GetDofList {{{1*/ -void Penta::GetDofList(int** pdoflist,int approximation_enum,int setenum){ - - int i,j,count=0; - int numberofdofs=0; - int* doflist=NULL; - - /*First, figure out size of doflist: */ - for(i=0;i<6;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); - - /*Allocate: */ - doflist=(int*)xmalloc(numberofdofs*sizeof(int)); - - /*Populate: */ - count=0; - for(i=0;i<6;i++){ - nodes[i]->GetDofList(doflist+count,approximation_enum,setenum); - count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); - } - - /*Assign output pointers:*/ - *pdoflist=doflist; -} -/*}}}*/ -/*FUNCTION Penta::GetDofList1 {{{1*/ -void Penta::GetDofList1(int* doflist){ - - int i; - for(i=0;i<6;i++) doflist[i]=nodes[i]->GetDofList1(); - -} -/*}}}*/ -/*FUNCTION Penta::GetElementType {{{1*/ -int Penta::GetElementType(){ - - /*return PentaRef field*/ - return this->element_type; -} -/*}}}*/ -/*FUNCTION Penta::GetHorizontalNeighboorSids {{{1*/ -int* Penta::GetHorizontalNeighboorSids(){ - - /*return PentaRef field*/ - return &this->horizontalneighborsids[0]; - -} -/*}}}*/ -/*FUNCTION Penta::GetLowerElement{{{1*/ -Penta* Penta::GetLowerElement(void){ - - Penta* upper_penta=NULL; - - upper_penta=(Penta*)verticalneighbors[0]; //first one (0) under, second one (1) above - - return upper_penta; -} -/*}}}*/ -/*FUNCTION Penta::GetNodeIndex {{{1*/ -int Penta::GetNodeIndex(Node* node){ - - _assert_(nodes); - for(int i=0;iGetInput(enumtype); - if (!input) _error_("Input %s not found in element",EnumToStringx(enumtype)); - - /*Checks in debugging mode*/ - _assert_(pvalue); - - /* Start looping on the number of vertices: */ - gauss=new GaussPenta(); - for (int iv=0;ivGaussVertex(iv); - input->GetParameterValue(&pvalue[iv],gauss); - } - - /*clean-up*/ - delete gauss; -} -/*}}}*/ -/*FUNCTION Penta::GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue) {{{1*/ -void Penta::GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue){ - - /*Intermediaries*/ - double value[NUMVERTICES]; - GaussPenta *gauss = NULL; - - /*Recover input*/ - Input* input=inputs->GetInput(enumtype); - - /*Checks in debugging mode*/ - _assert_(pvalue); - - /* Start looping on the number of vertices: */ - if (input){ - gauss=new GaussPenta(); - for (int iv=0;ivGaussVertex(iv); - input->GetParameterValue(&pvalue[iv],gauss); - } - } - else{ - for (int iv=0;ivGetInput(enumtype); - if(!input) _error_("No input of type %s found in tria",EnumToStringx(enumtype)); - - GaussPenta* gauss=new GaussPenta(); - gauss->GaussVertex(this->GetNodeIndex(node)); - - input->GetParameterValue(pvalue,gauss); - delete gauss; -} -/*}}}*/ -/*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 = Tr(sigma * eps) - * = Tr(sigma'* eps) - * = 2 * eps_eff * sigma'_eff - * = 4 * mu * eps_eff ^2*/ - *phi=4*pow(epsilon_eff,2.0)*viscosity; -} -/*}}}*/ -/*FUNCTION Penta::GetSidList{{{1*/ -void Penta::GetSidList(int* sidlist){ - - int i; - for(i=0;iGetSidList(); - -} -/*}}}*/ -/*FUNCTION Penta::GetSolutionFromInputs{{{1*/ -void Penta::GetSolutionFromInputs(Vec solution){ - - int analysis_type; - - /*retrive parameters: */ +/*FUNCTION Penta::GradjDragMacAyeal {{{1*/ +void Penta::GradjDragMacAyeal(Vec gradient){ + + /*Gradient is 0 if on shelf or not on bed*/ + if(IsOnShelf() || !IsOnBed()) return; + + /*Spawn tria*/ + Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. + tria->GradjDragMacAyeal(gradient); + delete tria->matice; delete tria; + +} /*}}}*/ +/*FUNCTION Penta::GradjDragPattyn {{{1*/ +void Penta::GradjDragPattyn(Vec gradient){ + + int i,j,ig; + int drag_type,analysis_type; + int doflist1[NUMVERTICES]; + double vx,vy,lambda,mu,alpha_complement,Jdet; + double bed,thickness,Neff,drag; + double xyz_list[NUMVERTICES][3]; + double xyz_list_tria[NUMVERTICES2D][3]={0.0}; + double dk[NDOF2]; + double grade_g[NUMVERTICES]={0.0}; + double grade_g_gaussian[NUMVERTICES]; + double basis[6]; + Friction* friction=NULL; + GaussPenta *gauss=NULL; + + /*Gradient is 0 if on shelf or not on bed*/ + if(IsOnShelf() || !IsOnBed()) return; + + /*Retrieve all inputs and parameters*/ parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */ - if (analysis_type==DiagnosticHorizAnalysisEnum){ - int approximation; - inputs->GetParameterValue(&approximation,ApproximationEnum); - if(approximation==StokesApproximationEnum || approximation==NoneApproximationEnum){ - GetSolutionFromInputsDiagnosticStokes(solution); - } - else if (approximation==MacAyealApproximationEnum || approximation==PattynApproximationEnum || approximation==HutterApproximationEnum){ - GetSolutionFromInputsDiagnosticHoriz(solution); - } - else if (approximation==MacAyealPattynApproximationEnum || approximation==PattynStokesApproximationEnum || approximation==MacAyealStokesApproximationEnum){ - return; //the elements around will create the solution - } - } - else if(analysis_type==DiagnosticHutterAnalysisEnum){ - GetSolutionFromInputsDiagnosticHutter(solution); - } - else if(analysis_type==DiagnosticVertAnalysisEnum){ - GetSolutionFromInputsDiagnosticVert(solution); - } - else if(analysis_type==ThermalAnalysisEnum){ - GetSolutionFromInputsThermal(solution); - } - else if(analysis_type==EnthalpyAnalysisEnum){ - GetSolutionFromInputsEnthalpy(solution); - } - else{ - _error_("analysis: %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type)); - } -} -/*}}}*/ -/*FUNCTION Penta::GetSolutionFromInputsDiagnosticHoriz{{{1*/ -void Penta::GetSolutionFromInputsDiagnosticHoriz(Vec solution){ - - const int numdof=NDOF2*NUMVERTICES; - - int i; - int approximation; - int* doflist=NULL; - double vx,vy; - double values[numdof]; - GaussPenta* gauss; - - /*Get approximation enum and dof list: */ - inputs->GetParameterValue(&approximation,ApproximationEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - - /*If the element is a coupling, do nothing: every node is also on an other elements - * (as coupling is between MacAyeal and Pattyn) so the other element will take care of it*/ - GetDofList(&doflist,approximation,GsetEnum); - - /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ - /*P1 element only for now*/ - gauss=new GaussPenta(); - for(i=0;iGaussVertex(i); + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + for(i=0;iGetInput(AdjointxEnum); _assert_(adjointx_input); + Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* dragcoefficient_input=inputs->GetInput(DragCoefficientEnum); _assert_(dragcoefficient_input); + + /*Build frictoin element, needed later: */ + inputs->GetParameterValue(&drag_type,DragTypeEnum); + friction=new Friction("2d",inputs,matpar,analysis_type); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussPenta(0,1,2,4); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetTriaJacobianDeterminant(&Jdet, &xyz_list_tria[0][0],gauss); + GetNodalFunctionsP1(basis, gauss); + + /*Build alpha_complement_list: */ + if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum); + else alpha_complement=0; + + dragcoefficient_input->GetParameterValue(&drag, gauss); + adjointx_input->GetParameterValue(&lambda, gauss); + adjointy_input->GetParameterValue(&mu, gauss); vx_input->GetParameterValue(&vx,gauss); vy_input->GetParameterValue(&vy,gauss); - values[i*NDOF2+0]=vx; - values[i*NDOF2+1]=vy; - } - - /*Add value to global vector*/ - VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); - - /*Free ressources:*/ + dragcoefficient_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss); + + /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ + for (i=0;iweight*basis[i]; /*basis are 0 for the 3 upper nodes*/ + } + + /*Add gradje_g_gaussian vector to gradje_g: */ + for(i=0;iGetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - - /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ - /*P1 element only for now*/ - gauss=new GaussPenta(); - for(i=0;iGaussVertex(i); - vx_input->GetParameterValue(&vx,gauss); - vy_input->GetParameterValue(&vy,gauss); - values[i*NDOF2+0]=vx; - values[i*NDOF2+1]=vy; - } - - /*Add value to global vector*/ - VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); - - /*Free ressources:*/ + delete friction; +} +/*}}}*/ +/*FUNCTION Penta::GradjDragStokes {{{1*/ +void Penta::GradjDragStokes(Vec gradient){ + + int i,j,ig; + int drag_type,analysis_type; + int doflist1[NUMVERTICES]; + double bed,thickness,Neff; + double lambda,mu,xi,Jdet,vx,vy,vz; + double alpha_complement,drag; + double surface_normal[3],bed_normal[3]; + double xyz_list[NUMVERTICES][3]; + double xyz_list_tria[NUMVERTICES2D][3]={0.0}; + double dk[NDOF2]; + double basis[6]; + double grade_g[NUMVERTICES]={0.0}; + double grade_g_gaussian[NUMVERTICES]; + Friction* friction=NULL; + GaussPenta* gauss=NULL; + + /*Gradient is 0 if on shelf or not on bed*/ + if(IsOnShelf() || !IsOnBed()) return; + + /*Retrieve all inputs and parameters*/ + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + for(i=0;iGetParameterValue(&drag_type,DragTypeEnum); + Input* drag_input =inputs->GetInput(DragCoefficientEnum); _assert_(drag_input); + Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input =inputs->GetInput(VzEnum); _assert_(vz_input); + Input* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input); + Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input); + Input* adjointz_input=inputs->GetInput(AdjointzEnum); _assert_(adjointz_input); + + /*Build frictoin element, needed later: */ + inputs->GetParameterValue(&drag_type,DragTypeEnum); + friction=new Friction("2d",inputs,matpar,analysis_type); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussPenta(0,1,2,4); + for(ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + /*Recover alpha_complement and drag: */ + if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum); + else alpha_complement=0; + drag_input->GetParameterValue(&drag,gauss); + + /*recover lambda mu and xi: */ + adjointx_input->GetParameterValue(&lambda,gauss); + adjointy_input->GetParameterValue(&mu ,gauss); + adjointz_input->GetParameterValue(&xi ,gauss); + + /*recover vx vy and vz: */ + vx_input->GetParameterValue(&vx, gauss); + vy_input->GetParameterValue(&vy, gauss); + vz_input->GetParameterValue(&vz, gauss); + + /*Get normal vector to the bed */ + SurfaceNormal(&surface_normal[0],xyz_list_tria); + + bed_normal[0]=-surface_normal[0]; //Function is for upper 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]; + + /* Get Jacobian determinant: */ + GetTriaJacobianDeterminant(&Jdet,&xyz_list_tria[0][0],gauss); + GetNodalFunctionsP1(basis, gauss); + + /*Get k derivative: dk/dx */ + drag_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss); + + /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ + for (i=0;iweight*basis[i]; + } + + /*Add gradje_g_gaussian vector to gradje_g: */ + for( i=0; iGetInput(VzEnum); _assert_(vz_input); - - /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ - /*P1 element only for now*/ - gauss=new GaussPenta(); - for(i=0;iGaussVertex(i); - vz_input->GetParameterValue(&vz,gauss); - values[i]=vz; - } - - /*Add value to global vector*/ - VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); - - /*Free ressources:*/ - delete gauss; - xfree((void**)&doflist); -} -/*}}}*/ -/*FUNCTION Penta::GetSolutionFromInputsDiagnosticStokes{{{1*/ -void Penta::GetSolutionFromInputsDiagnosticStokes(Vec solution){ - - const int numdof=NDOF4*NUMVERTICES; - - int i; - int* doflist=NULL; - double vx,vy,vz,p; - double stokesreconditioning; - double values[numdof]; - GaussPenta *gauss; - - /*Get dof list: */ - GetDofList(&doflist,StokesApproximationEnum,GsetEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); - Input* p_input =inputs->GetInput(PressureEnum); _assert_(p_input); - - /*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*/ - gauss=new GaussPenta(); - for(i=0;iGaussVertex(i); - vx_input->GetParameterValue(&vx,gauss); - vy_input->GetParameterValue(&vy,gauss); - vz_input->GetParameterValue(&vz,gauss); - p_input ->GetParameterValue(&p ,gauss); - values[i*NDOF4+0]=vx; - values[i*NDOF4+1]=vy; - values[i*NDOF4+2]=vz; - values[i*NDOF4+3]=p/stokesreconditioning; - } - - /*Add value to global vector*/ - VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); - - /*Free ressources:*/ - delete gauss; - xfree((void**)&doflist); -} -/*}}}*/ -/*FUNCTION Penta::GetSolutionFromInputsThermal{{{1*/ -void Penta::GetSolutionFromInputsThermal(Vec solution){ - - const int numdof=NDOF1*NUMVERTICES; - - int i; - int* doflist=NULL; - double values[numdof]; - double temp; - GaussPenta *gauss=NULL; - - /*Get dof list: */ - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - Input* t_input=inputs->GetInput(TemperatureEnum); _assert_(t_input); - - gauss=new GaussPenta(); - for(i=0;iGaussVertex(i); - t_input->GetParameterValue(&temp,gauss); - values[i]=temp; - } - - /*Add value to global vector*/ - VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); - - /*Free ressources:*/ - delete gauss; - xfree((void**)&doflist); -} -/*}}}*/ -/*FUNCTION Penta::GetSolutionFromInputsEnthalpy{{{1*/ -void Penta::GetSolutionFromInputsEnthalpy(Vec solution){ - - const int numdof=NDOF1*NUMVERTICES; - - int i; - int* doflist=NULL; - double values[numdof]; - double enthalpy; - GaussPenta *gauss=NULL; - - /*Get dof list: */ - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - Input* h_input=inputs->GetInput(EnthalpyEnum); _assert_(h_input); - - gauss=new GaussPenta(); - for(i=0;iGaussVertex(i); - h_input->GetParameterValue(&enthalpy,gauss); - values[i]=enthalpy; - } - - /*Add value to global vector*/ - VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); - - /*Free ressources:*/ - delete gauss; - xfree((void**)&doflist); -} -/*}}}*/ -/*FUNCTION Penta::GetStabilizationParameter {{{1*/ -double Penta::GetStabilizationParameter(double u, double v, double w, double diameter, double rho_ice, double heatcapacity, double thermalconductivity){ - /*Compute stabilization parameter*/ - - double normu; - double tau_parameter; - - normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5); - if(normu*diameter/(3*2*thermalconductivity/(rho_ice*heatcapacity))<1){ - tau_parameter=pow(diameter,2)/(3*2*2*thermalconductivity/(rho_ice*heatcapacity)); - } - else tau_parameter=diameter/(2*normu); - - return tau_parameter; -} -/*}}}*/ -/*FUNCTION Penta::GetStrainRate3dPattyn{{{1*/ -void Penta::GetStrainRate3dPattyn(double* epsilon,double* xyz_list, GaussPenta* 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){ - _error_("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, GaussPenta* 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){ - _error_("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*)verticalneighbors[1]; //first one under, second one above - - return upper_penta; -} -/*}}}*/ -/*FUNCTION Penta::GetVectorFromInputs{{{1*/ -void Penta::GetVectorFromInputs(Vec vector,int input_enum){ - - int doflist1[NUMVERTICES]; - - /*Get out if this is not an element input*/ - if (!IsInput(input_enum)) return; - - /*Prepare index list*/ - this->GetDofList1(&doflist1[0]); - - /*Get input (either in element or material)*/ - Input* input=inputs->GetInput(input_enum); - if(!input) _error_("Input %s not found in element",EnumToStringx(input_enum)); - - /*We found the enum. Use its values to fill into the vector, using the vertices ids: */ - input->GetVectorFromInputs(vector,&doflist1[0]); -} -/*}}}*/ -/*FUNCTION Penta::GetZcoord {{{1*/ -double Penta::GetZcoord(GaussPenta* gauss){ - - int i; - double z; - double xyz_list[NUMVERTICES][3]; - double z_list[NUMVERTICES]; - - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - for(i=0;iInputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum); + + /*Collapse element to the base*/ + Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face). + tria->GradjBMacAyeal(gradient); + delete tria->matice; delete tria; + + /*delete Average B*/ + this->matice->inputs->DeleteInput(RheologyBbarEnum); + +} /*}}}*/ +/*FUNCTION Penta::GradjBbarPattyn {{{1*/ +void Penta::GradjBbarPattyn(Vec gradient){ + + /*Gradient is computed on bed only (Bbar)*/ + if (!IsOnBed()) return; + + /*Depth Average B*/ + this->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum); + + /*Collapse element to the base*/ + Tria* tria=(Tria*)SpawnTria(0,1,2); + tria->GradjBMacAyeal(gradient); //We use MacAyeal as an estimate for now + delete tria->matice; delete tria; + + /*delete Average B*/ + this->matice->inputs->DeleteInput(RheologyBbarEnum); +} /*}}}*/ +/*FUNCTION Penta::GradjBbarStokes {{{1*/ +void Penta::GradjBbarStokes(Vec gradient){ + + /*Gradient is computed on bed only (Bbar)*/ + if (!IsOnBed()) return; + + /*Depth Average B*/ + this->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum); + + /*Collapse element to the base*/ + Tria* tria=(Tria*)SpawnTria(0,1,2); + tria->GradjBMacAyeal(gradient); //We use MacAyeal as an estimate for now + delete tria->matice; delete tria; + + /*delete Average B*/ + this->matice->inputs->DeleteInput(RheologyBbarEnum); +} /*}}}*/ /*FUNCTION Penta::Sid {{{1*/ int Penta::Sid(){ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 8653) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 8654) @@ -88,4 +88,10 @@ void GetVectorFromInputs(Vec vector,int NameEnum); void Gradj(Vec gradient,int control_type); + void GradjDragMacAyeal(Vec gradient); + void GradjDragPattyn(Vec gradient); + void GradjDragStokes(Vec gradient); + void GradjBbarMacAyeal(Vec gradient); + void GradjBbarPattyn(Vec gradient); + void GradjBbarStokes(Vec gradient); int Sid(); void InputControlUpdate(double scalar,bool save_parameter); Index: /issm/trunk/src/c/objects/Elements/PentaRef.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/PentaRef.cpp (revision 8653) +++ /issm/trunk/src/c/objects/Elements/PentaRef.cpp (revision 8654) @@ -1044,5 +1044,4 @@ *Jdet=SQRT3/6.0*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2.0)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2.0)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2.0),0.5); if(*Jdet<0) _error_("negative jacobian determinant!"); - } /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 8653) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 8654) @@ -2993,5 +2993,4 @@ double xyz_list[NUMVERTICES][3]; double basis[3],epsilon[3]; - double dbasis[NDOF2][NUMVERTICES]; double grad[NUMVERTICES]={0.0}; GaussTria *gauss = NULL; @@ -3027,5 +3026,4 @@ GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); GetNodalFunctions(basis,gauss); - GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss); /*standard gradient dJ/dki*/ @@ -3050,5 +3048,4 @@ double bed,thickness,Neff,drag; double xyz_list[NUMVERTICES][3]; - double dbasis[NDOF2][NUMVERTICES]; double dk[NDOF2]; double grade_g[NUMVERTICES]={0.0}; @@ -3085,5 +3082,4 @@ GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); GetNodalFunctions(basis, gauss); - GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss); /*Build alpha_complement_list: */ @@ -3109,5 +3105,4 @@ } } - VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES); @@ -3115,108 +3110,4 @@ delete gauss; delete friction; -} -/*}}}*/ -/*FUNCTION Tria::GradjDragStokes {{{1*/ -void Tria::GradjDragStokes(Vec gradient){ - - int i,ig; - int drag_type,analysis_type; - int doflist1[NUMVERTICES]; - double bed,thickness,Neff; - double lambda,mu,xi,Jdet,vx,vy,vz; - double alpha_complement,drag; - double surface_normal[3],bed_normal[3]; - double xyz_list[NUMVERTICES][3]; - double dbasis[NDOF2][NUMVERTICES]; - double dk[NDOF2]; - double basis[3]; - double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ - double grade_g[NUMVERTICES]={0.0}; - double grade_g_gaussian[NUMVERTICES]; - Friction* friction=NULL; - GaussTria* gauss=NULL; - - /*retrive parameters: */ - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*retrieve inputs :*/ - inputs->GetParameterValue(&drag_type,DragTypeEnum); - Input* drag_input =inputs->GetInput(DragCoefficientEnum); _assert_(drag_input); - Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input =inputs->GetInput(VzEnum); _assert_(vz_input); - Input* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input); - Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input); - Input* adjointz_input=inputs->GetInput(AdjointzEnum); _assert_(adjointz_input); - - /*Get out if shelf*/ - if(IsOnShelf())return; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList1(&doflist1[0]); - - /*Build frictoin element, needed later: */ - inputs->GetParameterValue(&drag_type,DragTypeEnum); - friction=new Friction("2d",inputs,matpar,analysis_type); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(4); - for(ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - /*Recover alpha_complement and drag: */ - if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum); - else alpha_complement=0; - drag_input->GetParameterValue(&drag,gauss); - - /*recover lambda mu and xi: */ - adjointx_input->GetParameterValue(&lambda,gauss); - adjointy_input->GetParameterValue(&mu ,gauss); - adjointz_input->GetParameterValue(&xi ,gauss); - - /*recover vx vy and vz: */ - vx_input->GetParameterValue(&vx, gauss); - vy_input->GetParameterValue(&vy, gauss); - vz_input->GetParameterValue(&vz, gauss); - - /*Get normal vector to the bed */ - SurfaceNormal(&surface_normal[0],xyz_list); - - 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]; - - /* Get Jacobian determinant: */ - GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss); - - /* Get nodal functions value at gaussian point:*/ - GetNodalFunctions(basis, gauss); - - /*Get nodal functions derivatives*/ - GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss); - - /*Get k derivative: dk/dx */ - drag_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss); - - /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ - for (i=0;iweight*basis[i]; - } - - /*Add gradje_g_gaussian vector to gradje_g: */ - for( i=0; iGetParameterDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss); @@ -3335,5 +3224,4 @@ double thickness,Jdet; double basis[3]; - double dbasis[NDOF2][NUMVERTICES]; double Dlambda[2],dp[2]; double xyz_list[NUMVERTICES][3]; @@ -3357,5 +3245,4 @@ GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); GetNodalFunctions(basis, gauss); - GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss); adjoint_input->GetParameterDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss); @@ -5371,5 +5258,4 @@ double xyz_list[NUMVERTICES][3]; GaussTria *gauss = NULL; - double dbasis[NDOF2][NUMVERTICES]; double dH[2]; @@ -5391,5 +5277,4 @@ /* Get Jacobian determinant: */ GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss); /*Get parameters at gauss point*/ Index: /issm/trunk/src/c/objects/Loads/Friction.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Friction.cpp (revision 8653) +++ /issm/trunk/src/c/objects/Loads/Friction.cpp (revision 8654) @@ -185,5 +185,5 @@ } /*}}}*/ -/*FUNCTION Friction::GetAlphaComplement {{{1*/ +/*FUNCTION Friction::GetAlphaComplement(double* palpha_complement, GaussTria* gauss,int vxenum,int vyenum) {{{1*/ void Friction::GetAlphaComplement(double* palpha_complement, GaussTria* gauss,int vxenum,int vyenum){ @@ -239,10 +239,68 @@ if(vmag==0 && (s-1)<0) _error_("velocity is 0 ans (s-1)=%g<0, alpha_complement is Inf",s-1); - alpha_complement=pow(Neff,r)*pow(vmag,(s-1)); - _assert_(!isnan(alpha_complement)); + alpha_complement=pow(Neff,r)*pow(vmag,(s-1)); _assert_(!isnan(alpha_complement)); /*Assign output pointers:*/ *palpha_complement=alpha_complement; - +} +/*}}}*/ +/*FUNCTION Friction::GetAlphaComplement(double* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum) {{{1*/ +void Friction::GetAlphaComplement(double* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum){ + + /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. + * FrictionGetAlphaComplement is used in control methods on drag, and it computes: + * alpha_complement= Neff ^r * vel ^s*/ + + /*diverse: */ + int i; + double Neff; + double r,s; + double vx; + double vy; + double vmag; + double drag_p,drag_q; + double drag_coefficient; + double bed,thickness; + double gravity,rho_ice,rho_water; + double alpha_complement; + + /*Recover parameters: */ + inputs->GetParameterValue(&drag_p,DragPEnum); + inputs->GetParameterValue(&drag_q,DragQEnum); + this->GetParameterValue(&thickness, gauss,ThicknessEnum); + this->GetParameterValue(&bed, gauss,BedEnum); + this->GetParameterValue(&drag_coefficient, gauss,DragCoefficientEnum); + + /*Get material parameters: */ + gravity=matpar->GetG(); + rho_ice=matpar->GetRhoIce(); + rho_water=matpar->GetRhoWater(); + + + //compute r and q coefficients: */ + r=drag_q/drag_p; + s=1./drag_p; + + //From bed and thickness, compute effective pressure when drag is viscous: + Neff=gravity*(rho_ice*thickness+rho_water*bed); + + /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because + the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour + for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should + replace it by Neff=0 (ie, equival it to an ice shelf)*/ + if (Neff<0)Neff=0; + + //We need the velocity magnitude to evaluate the basal stress: + this->GetParameterValue(&vx, gauss,vxenum); + this->GetParameterValue(&vy, gauss,vyenum); + vmag=sqrt(pow(vx,2)+pow(vy,2)); + + /*Checks that s-1>0 if v=0*/ + if(vmag==0 && (s-1)<0) _error_("velocity is 0 ans (s-1)=%g<0, alpha_complement is Inf",s-1); + + alpha_complement=pow(Neff,r)*pow(vmag,(s-1)); _assert_(!isnan(alpha_complement)); + + /*Assign output pointers:*/ + *palpha_complement=alpha_complement; } /*}}}*/ Index: /issm/trunk/src/c/objects/Loads/Friction.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Friction.h (revision 8653) +++ /issm/trunk/src/c/objects/Loads/Friction.h (revision 8654) @@ -30,4 +30,5 @@ void GetAlpha2(double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum); void GetAlphaComplement(double* alpha_complement, GaussTria* gauss,int vxenum,int vyenum); + void GetAlphaComplement(double* alpha_complement, GaussPenta* gauss,int vxenum,int vyenum); void GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type); void GetParameterValue(double* pvalue,GaussPenta* gauss,int enum_type);