Changeset 5373

Timestamp:

08/18/10 12:12:46 (15 years ago)

Author:

seroussi

Message:

keep working on coupling

Location:

issm/trunk/src/c/objects

Files:

• Elements/Penta.cpp (modified) (12 diffs)
• Loads/Icefront.cpp (modified) (6 diffs)
• Loads/Icefront.h (modified) (1 diff)
• Materials/Matice.cpp (modified) (2 diffs)

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

@@ -773 +773 @@
                 }
                 else if(approximation==MacAyealPattynApproximationEnum){
-                        ISSMERROR("coupling MacAyeal/Pattyn not implemented yet");
                         CreateKMatrixDiagnosticMacAyeal( Kgg);
                         CreateKMatrixDiagnosticPattyn( Kgg);
@@ -2377 +2376 @@
         /* local element matrices: */
         double Ke_gg[numdofp][numdofm]={0.0}; //local element stiffness matrix 
+        double Ke_gg_transp[numdofm][numdofp]={0.0}; //local element stiffness matrix 
         double Ke_gg_gaussian[numdofp][numdofm]; //stiffness matrix evaluated at the gaussian point.
         double Jdet;
@@ -2490 +2490 @@
         } //for (ig1=0; ig1<num_area_gauss; ig1++)
 
-        /*Add Ke_gg to global matrix Kgg: */
-        // one need to be transposed
+        /*Add Ke_gg and its transpose to global matrix Kgg: */
+        MatrixTranspose(&Ke_gg_transp[0][0],&Ke_gg[0][0],12,6);
         MatSetValues(Kgg,numdofp,doflistp,numdofm,doflistm,(const double*)Ke_gg,ADD_VALUES);
-        MatSetValues(Kgg,numdofm,doflistm,numdofp,doflistp,(const double*)Ke_gg,ADD_VALUES);
+        MatSetValues(Kgg,numdofm,doflistm,numdofp,doflistp,(const double*)Ke_gg_transp,ADD_VALUES);
 
         //Deal with 2d friction at the bedrock interface
@@ -2505 +2505 @@
         }
 
-        delete tria;
-        delete pentabase;
+        delete tria->matice; delete tria;
 
         xfree((void**)&first_gauss_area_coord);
@@ -4973 +4972 @@
         double       vx;
         double       vy;
-
-        /*Get dof list: */
-        GetDofList(&doflist);
+        int          approximation;
+
+        /*Get approximation enum and dof list: */
+        inputs->GetParameterValue(&approximation,ApproximationEnum);
+
+        /*If the element is a coupling, do nothing: every grid is also on an other elements 
+         * (as coupling is between MacAyeal and Pattyn) so the other element will take care of it*/
+        if(approximation==MacAyealPattynApproximationEnum){
+         return;
+        }
+        else{
+                GetDofList(&doflist,approximation);
+        }
 
         /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
@@ -5458 +5467 @@
 
         /*Get dof list: */
-        GetDofList(&doflist);
+        GetDofList(&doflist,MacAyealApproximationEnum);
 
         /*Use the dof list to index into the solution vector: */
@@ -5529 +5538 @@
 
         const int    numvertices=6;
-        const int    numdofpervertex=4;
-        const int    solutionnumdof=2;
-        const int    numdof=solutionnumdof*numvertices;
-        int*         doflist=NULL;
-        int*         doflistbase=NULL;
+        const int    numvertices2d=3;
+        const int    numdofpervertex=2;
+        const int    numdof=numdofpervertex*numvertices;
+        const int    numdof2d=numdofpervertex*numvertices2d;
+        int*         doflistp=NULL;
+        int*         doflistm=NULL;
         double       macayeal_values[numdof];
         double       pattyn_values[numdof];
@@ -5556 +5566 @@
 
         /*Get dof listof this element (pattyn dofs) and of the penta at base (macayeal dofs): */
-        GetDofList(&doflist);
-        penta->GetDofList(&doflistbase);
+        GetDofList(&doflistp,PattynApproximationEnum);
+        penta->GetDofList(&doflistm,MacAyealApproximationEnum);
 
         /*Get node data: */
@@ -5563 +5573 @@
 
         /*Use the dof list to index into the solution vector: */
-        for(i=0;i<numdof;i++){
-                pattyn_values[i]=solution[doflist[i]];
-                macayeal_values[i]=solution[doflistbase[i]];
+        for(i=0;i<numdof2d;i++){
+                pattyn_values[i]=solution[doflistp[i]];
+                macayeal_values[i]=solution[doflistm[i]];
+        }
+        for(i=numdof2d;i<numdof;i++){
+                pattyn_values[i]=solution[doflistp[i]];
+                macayeal_values[i]=macayeal_values[i-numdof2d];
         }
 
@@ -5588 +5602 @@
 
         /*Now Compute vel*/
-        for(i=0;i<numvertices;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
+        for(i=0;i<numvertices;i++) {
+                vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
+        }
 
         /*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D, 
@@ -5612 +5628 @@
 
         /*Free ressources:*/
-        xfree((void**)&doflist);
-        xfree((void**)&doflistbase);
+        xfree((void**)&doflistp);
+        xfree((void**)&doflistm);
 }
 /*}}}*/
@@ -5641 +5657 @@
 
         /*Get dof list: */
-        GetDofList(&doflist);
+        GetDofList(&doflist,PattynApproximationEnum);
 
         /*Get node data: */

issm/trunk/src/c/objects/Loads/Icefront.cpp

@@ -448 +448 @@
 
         /* Get dof list and node coordinates: */
-        GetDofList(&doflist);
+        GetDofList(&doflist,MacAyealApproximationEnum);
         GetVerticesCoordinates(&xyz_list[0][0],nodes,numgrids);
         
@@ -571 +571 @@
 
         /* Get dof list and node coordinates: */
-        GetDofList(&doflist);
+        GetDofList(&doflist,PattynApproximationEnum);
         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
         
@@ -641 +641 @@
         normal_norm=norm(normal3);normal3[0]=normal3[0]/normal_norm;normal3[1]=normal3[1]/normal_norm;normal3[2]=normal3[2]/normal_norm;
         normal_norm=norm(normal4);normal4[0]=normal4[0]/normal_norm;normal4[1]=normal4[1]/normal_norm;normal4[2]=normal4[2]/normal_norm;
-
 
         //Compute load contribution for this quad:
@@ -796 +795 @@
 /*}}}*/
 /*FUNCTION Icefront::GetDofList {{{1*/
-void  Icefront::GetDofList(int** pdoflist){
+void  Icefront::GetDofList(int** pdoflist,int approximation_enum){
 
         int i,j;
@@ -825 +824 @@
         /*Figure out size of doflist: */
         for(i=0;i<numberofnodes;i++){
-                numberofdofs+=nodes[i]->GetNumberOfDofs();
+                numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum);
         }
 
@@ -834 +833 @@
         count=0;
         for(i=0;i<numberofnodes;i++){
-                nodes[i]->GetDofList(doflist+count);
-                count+=nodes[i]->GetNumberOfDofs();
+                nodes[i]->GetDofList(doflist+count,approximation_enum);
+                count+=nodes[i]->GetNumberOfDofs(approximation_enum);
         }
 

issm/trunk/src/c/objects/Loads/Icefront.h

@@ -75 +75 @@
                 void  CreatePVectorDiagnosticHorizQuad( Vec pg);
                 void  CreatePVectorDiagnosticStokes( Vec pg);
-                void  GetDofList(int** pdoflist);
+                void  GetDofList(int** pdoflist,int approximation_enum=0);
                 void  QuadPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, 
                                               double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list);

issm/trunk/src/c/objects/Materials/Matice.cpp

@@ -378 +378 @@
                 if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && 
                                 (epsilon[3]==0) && (epsilon[4]==0)){
-                        viscosity3d=pow((double)10,(double)14);
+                        viscosity3d=0.5*pow((double)10,(double)14);
                 }
                 else{
@@ -450 +450 @@
                 if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && 
                                 (epsilon[3]==0) && (epsilon[4]==0) && (epsilon[5]==0)){
-                        viscosity3d=pow((double)10,(double)14);
+                        viscosity3d=0.5*pow((double)10,(double)14);
                 }
                 else{