Changeset 4924

Timestamp:

08/02/10 12:58:09 (15 years ago)

Author:

seroussi

Message:

PentaHutter done in Penta (still not working)

File:

• issm/trunk/src/c/objects/Elements/Penta.cpp (modified) (2 diffs)

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

@@ -3318 +3318 @@
 void  Penta::CreatePVectorDiagnosticHutter(Vec pg){
 
-        /*Collapsed formulation: */
-        Beam*  beam=NULL;
-        int    i;
+        int i,j,k;
+        
+        const int numgrids=6;
+        const int NDOF2=2;
+        const int numdofs=NDOF2*numgrids;
+        int       doflist[numdofs];
+        int       numberofdofspernode;
+        double    pe_g[numdofs]={0.0};
+        double    xyz_list[numgrids][3];
+        double    z_list[numgrids];
+        double    z_segment[2];
+        double    slope2,constant_part;
+        double    gauss[numdofs][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+        int       node0,node1;
+        BeamRef*  beam=NULL;
+
+        /*gaussian points: */
+        int      num_gauss;
+        double*  segment_gauss_coord=NULL;
+        double   gauss_coord;
+        double*  gauss_weights=NULL;
+        double   gauss_weight;
+        int      ig;
+        double   slope[2];
+
+        double   z_g;
+        double   rho_ice,gravity,n,B;
+        double   Jdet;
+        double   ub,vb;
+        double   surface,thickness;
+        double   slopex,slopey;
 
         /*flags: */
         bool onwater;
+        bool onbed;
+        bool onsurface;
+        int  connectivity[2];
+
+        /*Inputs*/
+        Input* thickness_input;
+        Input* surface_input;
+        Input* slopex_input;
+        Input* slopey_input;
 
         /*recover some inputs: */
@@ -3331 +3368 @@
         if(onwater)return;
 
-        /*Spawn 3 beam elements: */
+        /*recover doflist: */
+        GetDofList(&doflist[0],&numberofdofspernode);
+
+        /*recover some inputs: */
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
+
+        /*recover parameters: */
+        rho_ice=matpar->GetRhoIce();
+        gravity=matpar->GetG();
+        n=matice->GetN();
+        B=matice->GetB();
+
+        //Initilaize beam
+        beam=new BeamRef();
+
+        //recover extra inputs
+        thickness_input=inputs->GetInput(ThicknessEnum);
+        surface_input=inputs->GetInput(SurfaceEnum);
+        slopex_input=inputs->GetInput(SurfaceSlopeXEnum);
+        slopey_input=inputs->GetInput(SurfaceSlopeYEnum);
+
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+        for(i=0;i<numgrids;i++)z_list[i]=xyz_list[i][2];
+
+        //Get gaussian points and weights. order 2 since we have a product of 2 nodal functions
+        num_gauss=3;
+        GaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);
+
+        /*Loop on the three segments*/
         for(i=0;i<3;i++){
-                beam=(Beam*)SpawnBeam(i,i+3); //[0 3], [1 4] and [2 5] are the four vertical edges of the Penta
-                beam->CreatePVector(pg);
-        }
-
-        /*clean up*/
+                slopex_input->GetParameterValue(&slopex, &gauss[i][0]);
+                slopey_input->GetParameterValue(&slopey, &gauss[i][0]);
+                surface_input->GetParameterValue(&surface, &gauss[i][0]);
+                thickness_input->GetParameterValue(&thickness, &gauss[i][0]);
+
+                //compute slope2 slopex and slopey
+                slope2=pow(slopex,2)+pow(slopey,2);
+
+                //%compute constant_part
+                constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2));
+
+                z_segment[0]=z_list[i];
+                z_segment[1]=z_list[i+3];
+
+                node0=i;
+                node1=i+3;
+
+                connectivity[0]=nodes[node0]->GetConnectivity();
+                connectivity[1]=nodes[node1]->GetConnectivity();
+
+                /*Loop on the Gauss points: */
+                for(ig=0;ig<num_gauss;ig++){
+
+                        //Pick up the gaussian point and its weight:
+                        gauss_weight=gauss_weights[ig];
+                        gauss_coord=segment_gauss_coord[ig];
+
+                        beam->GetParameterValue(&z_g, &z_segment[0],gauss_coord);
+
+                        //Get Jacobian determinant:
+                        beam->GetJacobianDeterminant(&Jdet, &z_segment[0],gauss_coord);
+
+                        if (onsurface){
+                                for(j=0;j<NDOF2;j++){
+                                        pe_g[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss_weight/(double)connectivity[1];
+                                }
+                        }
+                        else{//connectivity is too large, should take only half on it
+                                for(j=0;j<NDOF2;j++){
+                                        pe_g[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss_weight*2/(double)connectivity[1];
+                                }
+                        }
+                }
+
+                //Deal with lower surface
+                if (onbed){
+
+                        //compute ub
+                        constant_part=-1.58*pow((double)10.0,-(double)10.0)*rho_ice*gravity*thickness;
+                        ub=constant_part*slopex;
+                        vb=constant_part*slopey;
+
+                        //Add to pe: 
+                        pe_g[2*node0]+=ub/(double)connectivity[0];
+                        pe_g[2*node0+1]+=vb/(double)connectivity[0];
+                }
+        }
+
+        /*Add pe_g to global vector pg: */
+        VecSetValues(pg,numdofs,doflist,(const double*)pe_g,ADD_VALUES);
+
+        /*Clean up */
         delete beam;
-
 }
 /*}}}*/