Changeset 5810

Timestamp:

09/14/10 17:14:45 (15 years ago)

Author:

seroussi

Message:

added PVectorCouplingPattynStokes

Location:

issm/trunk/src/c/objects/Elements

Files:

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

issm/trunk/src/c/objects/Elements/Penta.cpp

@@ -809 +809 @@
                         CreatePVectorDiagnosticPattyn( pg);
                         CreatePVectorDiagnosticStokes( pg);
-                        //CreatePVectorDiagnosticPattynStokes( pg);
+                        CreatePVectorCouplingPattynStokes( pg);
                 }
                 else ISSMERROR("Approximation %s not supported yet",EnumToString(approximation));
@@ -3160 +3160 @@
 }
 /*}}}*/
+/*FUNCTION Penta::CreatePVectorCouplingPattynStokes {{{1*/
+void Penta::CreatePVectorCouplingPattynStokes( Vec pg){
+
+        /*indexing: */
+        int i,j;
+        const int numdof=NUMVERTICES*NDOF4;
+        int*      doflist=NULL;
+
+        /*parameters: */
+        double             xyz_list_tria[NUMVERTICES2D][3];
+        double         xyz_list[NUMVERTICES][3];
+        double             bed_normal[3];
+
+        /* gaussian points: */
+        int     ig;
+        GaussPenta *gauss=NULL;
+
+        double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
+        double  viscosity, w, alpha2_gauss;
+        double  dw[3];
+
+        /*matrices: */
+        double     Pe_gaussian[24]={0.0}; //for the six nodes
+        double     l1l6[6]; //for the six nodes of the penta
+        double     dh1dh6[3][6]; //for the six nodes of the penta
+        double     Jdet;
+        double     Jdet2d;
+
+        Tria*     tria=NULL;
+        Friction* friction=NULL;
+
+        /*parameters: */
+        double stokesreconditioning;
+        int    approximation;
+        int    analysis_type;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&approximation,ApproximationEnum);
+
+        /*retrieve some parameters: */
+        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+        this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
+
+        /*If on water or not Stokes, skip load: */
+        if(IsOnWater() || approximation!=PattynStokesApproximationEnum) return;
+
+        /* Get node coordinates and dof list: */
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+        GetDofList(&doflist,StokesApproximationEnum,GsetEnum);
+
+        /*Retrieve all inputs we will be needing: */
+        Input* vx_input=inputs->GetInput(VxEnum);               ISSMASSERT(vx_input);
+        Input* vy_input=inputs->GetInput(VyEnum);               ISSMASSERT(vy_input);
+        Input* vz_input=inputs->GetInput(VzEnum);               ISSMASSERT(vz_input);
+        Input* vzpattyn_input=inputs->GetInput(VzPattynEnum);   ISSMASSERT(vzpattyn_input);
+
+        /* Start  looping on the number of gaussian points: */
+        gauss=new GaussPenta(5,5);
+        for (ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                /*Compute strain rate and viscosity: */
+                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
+                matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
+
+                /* Get Jacobian determinant: */
+                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+
+                /* Get nodal functions */
+                GetNodalFunctionsP1(&l1l6[0], gauss);
+                GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
+
+                /*Compute the derivative of VzPattyn*/
+                vzpattyn_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
+
+                /* Build gaussian vector */
+                for(i=0;i<NUMVERTICES;i++){
+                        Pe_gaussian[i*NDOF4+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dh1dh6[2][i]/2;
+                        Pe_gaussian[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dh1dh6[2][i]/2;
+                        Pe_gaussian[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dh1dh6[0][i]+dw[1]*dh1dh6[1][i]+dw[2]*dh1dh6[2][i])/2;
+                        Pe_gaussian[i*NDOF4+3]+=Jdet*gauss->weight*stokesreconditioning*dw[2]*l1l6[i];
+                }
+        }
+        delete gauss;
+
+        /*Deal with 2d friction at the bedrock interface: */
+        if(IsOnBed() && !IsOnShelf()){
+
+                for(i=0;i<NUMVERTICES2D;i++){
+                        for(j=0;j<3;j++){
+                                xyz_list_tria[i][j]=xyz_list[i][j];
+                        }
+                }
+
+                /*build friction object, used later on: */
+                friction=new Friction("3d",inputs,matpar,analysis_type);
+
+                /* Start looping on the number of gauss 2d (nodes on the bedrock) */
+                gauss=new GaussPenta(0,1,2,2);
+                for(ig=gauss->begin();ig<gauss->end();ig++){
+
+                        gauss->GaussPoint(ig);
+
+                        /*Get the Jacobian determinant */
+                        GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
+
+                        /*Get L matrix : */
+                        GetNodalFunctionsP1(l1l6, gauss);
+
+                        /*Get normal vecyor to the bed */
+                        BedNormal(&bed_normal[0],xyz_list_tria);
+
+                        /*Get Viscosity*/
+                        this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
+                        matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
+
+                        /*Get friction*/
+                        friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum);
+
+                        /*Get w and its derivatives*/
+                        vzpattyn_input->GetParameterValue(&w, gauss);
+                        vzpattyn_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
+
+                        for(i=0;i<NUMVERTICES2D;i++){
+                                Pe_gaussian[i*NDOF4+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+viscosity*dw[2]*bed_normal[0])*l1l6[i];
+                                Pe_gaussian[i*NDOF4+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+viscosity*dw[2]*bed_normal[1])*l1l6[i];
+                                Pe_gaussian[i*NDOF4+2]+=Jdet2d*gauss->weight*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*l1l6[i];
+                        }
+                }
+                /*Free ressources:*/
+                delete gauss;
+        } //if ( (IsOnBed()==1) && (IsOnShelf()==1))
+
+        /*Add Pe_reduced to global vector pg: */
+        VecSetValues(pg,numdof,doflist,(const double*)Pe_gaussian,ADD_VALUES);
+
+        /*Free ressources:*/
+        xfree((void**)&doflist);
+
+}
+/*}}}*/
 /*FUNCTION Penta::CreatePVectorDiagnosticHutter{{{1*/
 void  Penta::CreatePVectorDiagnosticHutter(Vec pg){

issm/trunk/src/c/objects/Elements/Penta.h

@@ -135 +135 @@
                 void      CreatePVectorAdjointMacAyeal( Vec pg);
                 void      CreatePVectorAdjointPattyn( Vec pg);
+                void      CreatePVectorCouplingPattynStokes( Vec pg);
                 void      CreatePVectorDiagnosticHutter( Vec pg);
                 void      CreatePVectorDiagnosticMacAyeal( Vec pg);