Changeset 16996

Timestamp:

12/03/13 11:28:58 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: cleaning up

Location:

issm/trunk-jpl/src/c/classes/Elements

Files:

• Penta.cpp (modified) (1 diff)
• Penta.h (modified) (2 diffs)
• Tria.cpp (modified) (2 diffs)
• Tria.h (modified) (2 diffs)

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

@@ -4520 +4520 @@
 }
 /*}}}*/
-/*FUNCTION Penta::KMatrixGLSstabilization{{{*/
-void Penta::KMatrixGLSstabilization(ElementMatrix* Ke){
-
-        int        numdof  = NUMVERTICES*NDOF4;
-
-        /*Intermediaries */
-        int        i,j;
-        IssmDouble Jdet,viscosity,FSreconditioning,diameter,rigidity;
-        IssmDouble xyz_list[NUMVERTICES][3];
-        IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
-        GaussPenta *gauss=NULL;
-
-        /*Stabilization*/
-        IssmDouble D_scalar;
-        IssmDouble dbasis[3][6];
-        IssmDouble dmu[3];
-        IssmDouble mu;
-        IssmDouble dnodalbasis[6][6][3];
-        IssmDouble SU[6][4][4];
-        IssmDouble SW[6][4][4];
-        int p,q,ii;
-        int c=3; //index of pressure
-
-        /*Retrieve all inputs and parameters*/
-        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-        parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
-        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);
-
-        /*Find minimal length and B*/
-        rigidity=material->GetB();
-        diameter=MinEdgeLength(&xyz_list[0][0]);
-
-                gauss=new GaussPenta();
-                for(int iv=0;iv<6;iv++){
-                        gauss->GaussVertex(iv);
-                        GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
-                        for(i=0;i<6;i++){
-                                for(j=0;j<3;j++){
-                                        dnodalbasis[i][iv][j] = dbasis[j][i];
-                                }
-                        }
-                }
-                delete gauss;
-
-        /* Start  looping on the number of gaussian points: */
-        gauss=new GaussPenta(5,5);
-        for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-                gauss->GaussPoint(ig);
-
-                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-
-                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
-                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
-
-                D_scalar=gauss->weight*Jdet;
-
-                /*Add stabilization*/
-                GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
-                dmu[0]=0.; dmu[1]=0.; dmu[2]=0.;
-                mu = 2.*viscosity;
-                for(p=0;p<6;p++) for(i=0;i<4;i++) for(j=0;j<4;j++){
-                        SU[p][i][j]=0.;
-                        SW[p][i][j]=0.;
-                }
-                for(p=0;p<6;p++){
-                        for(i=0;i<3;i++){
-                                SU[p][i][c] += dbasis[i][p];
-                                SW[p][c][i] += dbasis[i][p];
-                                for(j=0;j<3;j++){
-                                        SU[p][i][i] += -dmu[j]*dbasis[j][p];
-                                        SU[p][i][j] += -dmu[i]*dbasis[i][p];
-                                        for(ii=0;ii<6;ii++){
-                                                SU[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii];
-                                                SU[p][i][j] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii];
-                                        }
-                                        SW[p][i][i] += -dmu[j]*dbasis[j][p];
-                                        SW[p][j][i] += -dmu[j]*dbasis[i][p];
-                                        for(ii=0;ii<6;ii++){
-                                                SW[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii];
-                                                SW[p][j][i] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii];
-                                        }
-                                }
-                        }
-                }
-                IssmDouble tau = 1./3.*pow(diameter,2)/(8.*2.*viscosity);
-                for(p=0;p<6;p++){
-                        for(q=0;q<6;q++){
-                                for(i=0;i<4;i++){
-                                        for(j=0;j<4;j++){
-                                                Ke->values[(p*4+i)*numdof+q*4+j] += gauss->weight*Jdet*tau*SW[p][i][0]*SU[q][0][j];
-                                                Ke->values[(p*4+i)*numdof+q*4+j] += gauss->weight*Jdet*tau*SW[p][i][1]*SU[q][1][j];
-                                                Ke->values[(p*4+i)*numdof+q*4+j] += gauss->weight*Jdet*tau*SW[p][i][2]*SU[q][2][j];
-                                        }
-                                }
-                        }
-                }
-        }
-
-        /*Clean up*/
-        delete gauss;
-}
-/*}}}*/
 #endif
 

issm/trunk-jpl/src/c/classes/Elements/Penta.h

@@ -212 +212 @@
                 void           GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype);
                 Node*          GetNode(int node_number);
-                void           GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss, Input* vx_input, Input* vy_input);
-                void           GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, Gauss* gauss, Input* vx_input, Input* vy_input, Input* vz_input);
                 Penta*         GetUpperElement(void);
                 Penta*         GetLowerElement(void);
@@ -257 +255 @@
                 ElementVector* CreateDVectorStressbalanceHoriz(void);
                 ElementVector* CreateDVectorStressbalanceFS(void);
-                void           KMatrixGLSstabilization(ElementMatrix* Ke);
                 #endif
 

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -3062 +3062 @@
 #endif
 
-#ifdef _HAVE_STRESSBALANCE_
-/*FUNCTION Tria::GetYcoord {{{*/
-IssmDouble Tria::GetYcoord(GaussTria* gauss){
-
-        IssmDouble y;
-        IssmDouble xyz_list[NUMVERTICES][3];
-        IssmDouble y_list[NUMVERTICES];
-
-        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-        for(int i=0;i<NUMVERTICES;i++) y_list[i]=xyz_list[i][1];
-        TriaRef::GetInputValue(&y,&y_list[0],gauss,P1Enum);
-
-        return y;
-}
-/*}}}*/
-#endif
-
 #ifdef _HAVE_CONTROL_
 /*FUNCTION Tria::BalancethicknessMisfit{{{*/
@@ -4905 +4888 @@
 #endif
 
-#ifdef _HAVE_DAMAGE_
-/*FUNCTION Tria::DamageEvolutionF{{{*/
-void Tria::DamageEvolutionF(IssmDouble* f){
-
-        /*Intermediaries */
-        IssmDouble c1,c2,c3,healing,stress_threshold;
-        IssmDouble s_xx,s_xy,s_yy;
-        IssmDouble J2s;
-        IssmDouble Xis;
-        IssmDouble Psi;
-        IssmDouble PosPsi;
-        IssmDouble NegPsi;
-        Input* damage_input=NULL;
-        Input* sigma_xx_input  = NULL;
-        Input* sigma_xy_input  = NULL;
-        Input* sigma_yy_input  = NULL;
-        GaussTria* gauss=NULL;
-        IssmDouble damage,sigma_xx,sigma_xy,sigma_yy;
-
-        /*retrieve parameters:*/
-        this->parameters->FindParam(&c1,DamageC1Enum);
-        this->parameters->FindParam(&c2,DamageC2Enum);
-        this->parameters->FindParam(&c3,DamageC3Enum);
-        this->parameters->FindParam(&healing,DamageHealingEnum);
-        this->parameters->FindParam(&stress_threshold,DamageStressThresholdEnum);
-
-        /*Compute stress tensor: */
-        this->ComputeStressTensor();
-
-        /*retrieve what we need: */
-        sigma_xx_input  = inputs->GetInput(StressTensorxxEnum);  _assert_(sigma_xx_input);
-        sigma_xy_input  = inputs->GetInput(StressTensorxyEnum);  _assert_(sigma_xy_input);
-        sigma_yy_input  = inputs->GetInput(StressTensoryyEnum);  _assert_(sigma_yy_input);
-        damage_input  = this->material->inputs->GetInput(DamageDbarEnum);        _assert_(damage_input);
-
-        /*Damage evolution z mapping: */
-        gauss=new GaussTria();
-        J2s=0;
-        for (int iv=0;iv<NUMVERTICES;iv++){
-                gauss->GaussVertex(iv);
-                
-                damage_input->GetInputValue(&damage,gauss);
-                sigma_xx_input->GetInputValue(&sigma_xx,gauss);
-                sigma_xy_input->GetInputValue(&sigma_xy,gauss);
-                sigma_yy_input->GetInputValue(&sigma_yy,gauss);
-
-                s_xx=sigma_xx/(1-damage);
-                s_xy=sigma_xy/(1-damage);
-                s_yy=sigma_yy/(1-damage);
-
-                J2s=1.0/sqrt(2.0)*sqrt(pow(s_xx,2)+pow(s_yy,2)+2*pow(s_xy,2));
-                
-                Xis=sqrt(3.0)*J2s;
-
-                Psi=Xis-stress_threshold;
-
-                PosPsi=max(Psi,0.0);
-                NegPsi=max(-Psi,0.0);
-
-                f[iv]= c1* ( pow(PosPsi,c2)   -  healing * pow(NegPsi,c2) )  *   pow((1 - damage),-c3);
-
-        }
-        
-        /*Clean up and return*/
-        delete gauss;
-}
-/*}}}*/
-#endif
-
 #ifdef _HAVE_DAKOTA_
 /*FUNCTION Tria::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type);{{{*/

issm/trunk-jpl/src/c/classes/Elements/Tria.h

@@ -261 +261 @@
                 void           ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");};
 
-                #ifdef _HAVE_STRESSBALANCE_
-                void           PVectorGLSstabilization(ElementVector* pe);
-                IssmDouble     GetYcoord(GaussTria* gauss);
-                #endif
-
                 void UpdateConstraintsExtrudeFromBase(void);
                 void UpdateConstraintsExtrudeFromTop(void);
@@ -283 +278 @@
                 void           HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);
                 void           ComputeEPLThickness(void);
-
-                #endif
-
-                #ifdef _HAVE_DAMAGE_
-                void           DamageEvolutionF(IssmDouble* flist);
-                #endif
-
-
+                #endif
                 /*}}}*/