Changeset 1211

Timestamp:

07/02/09 09:00:31 (16 years ago)

Author:

Mathieu Morlighem

Message:

Added gradjDrag for Stokes CM

Location:

issm/trunk/src/c/objects

Files:

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

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

@@ -1225 +1225 @@
         
         /*Bail out if this element does not touch the bedrock: */
-        if (!onbed){
-                return;
-        }
-        else{
-                
+        if (!onbed) return;
+
+        if (sub_analysis_type==HorizAnalysisEnum()){
+
+                /*MacAyeal or Pattyn*/
                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
                 tria->GradjDrag( grad_g,inputs,analysis_type,sub_analysis_type);
@@ -1235 +1235 @@
                 return;
         }
+        else if (sub_analysis_type==StokesAnalysisEnum()){
+
+                /*Stokes*/
+                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
+                tria->GradjDragStokes( grad_g,inputs,analysis_type,sub_analysis_type);
+                delete tria;
+                return;
+        }
+        else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"));
 }
 

issm/trunk/src/c/objects/Tria.cpp

@@ -2266 +2266 @@
 
         /* parameters: */
-        double  viscosity2;
         double  dvx[NDOF2];
         double  dvy[NDOF2]; 
@@ -2430 +2429 @@
 
 #undef __FUNCT__ 
+#define __FUNCT__ "Tria::GradjDragStokes"
+void  Tria::GradjDragStokes(Vec grad_g,void* vinputs,int analysis_type,int sub_analysis_type){
+
+        int i;
+
+        /* node data: */
+        const int    numgrids=3;
+        const int    NDOF2=2;
+        double       xyz_list[numgrids][3];
+        int          doflist1[numgrids];
+        int          numberofdofspernode;
+
+        /* grid data: */
+        double vx_list[numgrids];
+        double vy_list[numgrids];
+        double vz_list[numgrids];
+        double vxvyvz_list[numgrids][3];
+        double adjx_list[numgrids];
+        double adjy_list[numgrids];
+        double adjz_list[numgrids];
+        double adjxyz_list[numgrids][3];
+
+        double drag;
+        int    dofs1[1]={0};
+        int    dofs3[3]={0,1,2};
+
+        /* 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* gauss_weights           =  NULL;
+        double  gauss_weight;
+        double  gauss_l1l2l3[3];
+
+        /* parameters: */
+        double  vx,vy,vz;
+        double  lambda,mu,xi;
+        double  bed,thickness,Neff;
+        double  surface_normal[3];
+        double  bed_normal[3];
+
+        /*drag: */
+        double  pcoeff,qcoeff;
+        double alpha_complement_list[numgrids];
+        double alpha_complement;
+
+        /*element vector at the gaussian points: */
+        double  grade_g[numgrids];
+        double  grade_g_gaussian[numgrids];
+
+        /* Jacobian: */
+        double Jdet;
+
+        /*nodal functions: */
+        double l1l2l3[3];
+
+        /* strain rate: */
+        double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
+
+        ParameterInputs* inputs=NULL;
+
+        /*recover pointers: */
+        inputs=(ParameterInputs*)vinputs;
+
+        /* Get node coordinates and dof list: */
+        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+        GetDofList1(&doflist1[0]);
+
+        /* Set grade_g to 0: */
+        for(i=0;i<numgrids;i++) grade_g[i]=0.0;
+
+        /* recover input parameters: */
+        inputs->Recover("drag",&k[0],1,dofs1,numgrids,(void**)nodes);
+        inputs->Recover("bed",&b[0],1,dofs1,numgrids,(void**)nodes);
+        inputs->Recover("thickness",&h[0],1,dofs1,numgrids,(void**)nodes);
+        if(!inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs3,numgrids,(void**)nodes)){
+                throw ErrorException(__FUNCT__,"missing velocity input parameter");
+        }
+        if(!inputs->Recover("adjoint",&adjxyz_list[0][0],3,dofs3,numgrids,(void**)nodes)){
+                throw ErrorException(__FUNCT__,"missing adjoint input parameter");
+        }
+
+        /*Initialize parameter lists: */
+        for(i=0;i<numgrids;i++){
+                vx_list[i]=vxvyvz_list[i][0];
+                vy_list[i]=vxvyvz_list[i][1];
+                vz_list[i]=vxvyvz_list[i][2];
+                adjx_list[i]=adjxyz_list[i][0];
+                adjy_list[i]=adjxyz_list[i][1];
+                adjz_list[i]=adjxyz_list[i][2];
+        }
+
+        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
+        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4);
+
+#ifdef _DEBUGELEMENTS_
+        if(my_rank==RANK && id==ELID){ 
+                printf("   gaussian points: \n");
+                for(i=0;i<num_gauss;i++){
+                        printf("    %g %g %g : %g\n",first_gauss_area_coord[i],second_gauss_area_coord[i],third_gauss_area_coord[i],gauss_weights[i]);
+                }
+        }
+#endif
+
+        /* Start  looping on the number of gaussian points: */
+        for (ig=0; ig<num_gauss; ig++){
+                /*Pick up the gaussian point: */
+                gauss_weight=*(gauss_weights+ig);
+                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
+                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
+                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
+
+                /*Build alpha_complement_list: */
+                if (!shelf && (friction_type==2)){
+
+                        /*Allocate friction object: */
+                        Friction* friction=NewFriction();
+
+                        /*Initialize all fields: */
+                        friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char));
+                        strcpy(friction->element_type,"2d");
+
+                        friction->gravity=matpar->GetG();
+                        friction->rho_ice=matpar->GetRhoIce();
+                        friction->rho_water=matpar->GetRhoWater();
+                        friction->K=&k[0];
+                        friction->bed=&b[0];
+                        friction->thickness=&h[0];
+                        friction->velocities=&vxvyvz_list[0][0];
+                        friction->p=p;
+                        friction->q=q;
+
+                        if(friction->p!=1) throw ErrorException(__FUNCT__,"non-linear friction not supported yet in control methods!");
+
+                        /*Compute alpha complement: */
+                        FrictionGetAlphaComplement(&alpha_complement_list[0],friction);
+
+                        /*Erase friction object: */
+                        DeleteFriction(&friction);
+                }
+                else{
+                        alpha_complement_list[0]=0;
+                        alpha_complement_list[1]=0;
+                        alpha_complement_list[2]=0;
+                }
+
+                /*Recover alpha_complement and k: */
+                GetParameterValue(&alpha_complement, &alpha_complement_list[0],gauss_l1l2l3);
+                GetParameterValue(&drag, &k[0],gauss_l1l2l3);
+#ifdef _DEBUG_ 
+                printf("Drag complement: %20.20lf Drag: %20.20lf\n",alpha_complement,drag);
+#endif
+
+                /*recover lambda mu and xi: */
+                GetParameterValue(&lambda, &adjx_list[0],gauss_l1l2l3);
+                GetParameterValue(&mu, &adjy_list[0],gauss_l1l2l3);
+                GetParameterValue(&xi, &adjz_list[0],gauss_l1l2l3);
+#ifdef _DEBUG_ 
+                printf("Adjoint vector %20.20lf %20.20lf\n",lambda,mu);
+#endif
+
+                /*recover vx vy and vz: */
+                GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3);
+                GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);
+                GetParameterValue(&vz, &vz_list[0],gauss_l1l2l3);
+#ifdef _DEBUG_ 
+                printf("Velocity vector %20.20lf %20.20lf\n",vx,vy);
+
+                /*Get normal vecyor 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];
+#endif
+
+                /* Get Jacobian determinant: */
+                GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
+#ifdef _DEBUG_ 
+                printf("Element id %i Jacobian determinant: %lf\n",TriaElementGetID(this),Jdet);
+#endif
+
+                /* Get nodal functions value at gaussian point:*/
+                GetNodalFunctions(l1l2l3, gauss_l1l2l3);
+
+                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
+                for (i=0;i<numgrids;i++){
+                        grade_g_gaussian[i]=(
+                                                -lambda*(2*drag*alpha_complement*(vx - vz*bed_normal[0]*bed_normal[2]))
+                                                -mu    *(2*drag*alpha_complement*(vy - vz*bed_normal[1]*bed_normal[2]))
+                                                -xi    *(2*drag*alpha_complement*(-vx*bed_normal[0]*bed_normal[2]-vy*bed_normal[1]*bed_normal[2]))
+                                                )*Jdet*gauss_weight*l1l2l3[i]; 
+                }
+
+                /*Add gradje_g_gaussian vector to gradje_g: */
+                for( i=0; i<numgrids; i++)grade_g[i]+=grade_g_gaussian[i];
+        }
+
+        /*Add grade_g to global vector grad_g: */
+        VecSetValues(grad_g,numgrids,doflist1,(const double*)grade_g,ADD_VALUES);
+
+cleanup_and_return: 
+        xfree((void**)&first_gauss_area_coord);
+        xfree((void**)&second_gauss_area_coord);
+        xfree((void**)&third_gauss_area_coord);
+        xfree((void**)&gauss_weights);
+
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::SurfaceNormal"
+
+void Tria::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){
+
+        int i;
+        double v13[3];
+        double v23[3];
+        double normal[3];
+        double normal_norm;
+
+        for (i=0;i<3;i++){
+                v13[i]=xyz_list[0][i]-xyz_list[2][i];
+                v23[i]=xyz_list[1][i]-xyz_list[2][i];
+        }
+
+        normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
+        normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
+        normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
+
+        normal_norm=sqrt( pow(normal[0],2)+pow(normal[1],2)+pow(normal[2],2) );
+
+        *(surface_normal)=normal[0]/normal_norm;
+        *(surface_normal+1)=normal[1]/normal_norm;
+        *(surface_normal+2)=normal[2]/normal_norm;
+
+}
+
+#undef __FUNCT__ 
 #define __FUNCT__ "Tria::GradjB"
 void  Tria::GradjB(Vec grad_g,void* vinputs,int analysis_type,int sub_analysis_type){

issm/trunk/src/c/objects/Tria.h

@@ -96 +96 @@
                 void  Gradj(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type);
                 void  GradjDrag(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type);
+                void  GradjDragStokes(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type);
+                void  SurfaceNormal(double* surface_normal, double xyz_list[3][3]);
                 void  GradjB(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type);
                 double Misfit(void* inputs,int analysis_type,int sub_analysis_type);