Changeset 5922

Timestamp:

09/21/10 14:11:25 (15 years ago)

Author:

Mathieu Morlighem

Message:

Some ElementVector, more to come

Location:

issm/trunk/src/c/objects/Elements

Files:

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

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

@@ -3771 +3771 @@
 void Penta::CreatePVectorSlope( Vec pg){
 
-        /*Collapsed formulation: */
-        Tria*  tria=NULL;
-
-        /*If on water, skip: */
-        if(IsOnWater())return;
-
-        /*Is this element on the bed? :*/
-        if(!IsOnBed())return;
-
-        /*Spawn Tria element from the base of the Penta: */
-        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-        tria->CreatePVector(pg,NULL);
+        if (!IsOnBed() || IsOnWater()) return;
+
+        /*Call Tria function*/
+        Tria* tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
+        ElementVector* pe=tria->CreatePVectorSlope();
         delete tria->matice; delete tria;
+        pe->AddToGlobal(pg,NULL);
+        delete pe;
+
+        /*clean up and return*/
         return;
 }

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

@@ -761 +761 @@
                         break;
                 case DiagnosticHutterAnalysisEnum:
-                        CreatePVectorDiagnosticHutter( pg);
+                        pe=CreatePVectorDiagnosticHutter();
+                        pe->AddToGlobal(pg,pf);
+                        delete pe;
                         break;
                 case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
-                        CreatePVectorSlope( pg);
+                        pe=CreatePVectorSlope();
+                        pe->AddToGlobal(pg,pf);
+                        delete pe;
                         break;
                 case PrognosticAnalysisEnum:
@@ -3687 +3691 @@
         /*Initialize Element vector and return if necessary*/
         if(IsOnWater()) return NULL;
-        ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters);
+        ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
 
         /*Retrieve all inputs and parameters*/
@@ -4197 +4201 @@
 /*}}}*/
 /*FUNCTION Tria::CreatePVectorDiagnosticHutter{{{1*/
-void  Tria::CreatePVectorDiagnosticHutter(Vec pg){
+ElementVector* Tria::CreatePVectorDiagnosticHutter(void){
 
         /*Collapsed formulation: */
@@ -4205 +4209 @@
         double    constant_part,ub,vb;
         double    rho_ice,gravity,n,B;
-        double    pe_g[numdofs];
         double    slope[2];
         double    thickness;
@@ -4212 +4215 @@
         GaussTria* gauss=NULL;
 
-        /*If on water, skip: */
-        if(IsOnWater())return;
-
-        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-        
-        /* Get parameters */
+        /*Initialize Element vector and return if necessary*/
+        if(IsOnWater()) return NULL;
+        ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters);
+
+        /*Retrieve all inputs and parameters*/
         rho_ice=matpar->GetRhoIce();
         gravity=matpar->GetG();
         n=matice->GetN();
         B=matice->GetBbar();
-
-        /* Get slopes and thickness */
         Input* slopex_input=inputs->GetInput(SurfaceSlopeXEnum); ISSMASSERT(slopex_input);
         Input* slopey_input=inputs->GetInput(SurfaceSlopeYEnum); ISSMASSERT(slopey_input);
@@ -4241 +4241 @@
                 slope2=pow(slope[0],2)+pow(slope[1],2);
 
-                //compute constant_part
                 constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2));
 
-                //compute ub
                 ub=-1.58*pow((double)10.0,(double)-10.0)*rho_ice*gravity*thickness*slope[0];
                 vb=-1.58*pow((double)10.0,(double)-10.0)*rho_ice*gravity*thickness*slope[1];
 
-                pe_g[2*i]=(ub-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0])/(double)connectivity;
-                pe_g[2*i+1]=(vb-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[1])/(double)connectivity;
-
-        }
-
-        VecSetValues(pg,numdofs,doflist,(const double*)pe_g,ADD_VALUES);
+                pe->values[2*i]  =(ub-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0])/(double)connectivity;
+                pe->values[2*i+1]=(vb-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[1])/(double)connectivity;
+        }
 
         /*Clean up and return*/
         delete gauss;
-        xfree((void**)&doflist);
+        return pe;
 }
 /*}}}*/
@@ -4415 +4410 @@
 /*}}}*/
 /*FUNCTION Tria::CreatePVectorSlope {{{1*/
-void Tria::CreatePVectorSlope( Vec pg){
-
-        int             i,j;
+ElementVector* Tria::CreatePVectorSlope(void){
 
         /* node data: */
         const int    numdof=NDOF1*NUMVERTICES;
+        int             i,j;
+        int     ig;
         double       xyz_list[NUMVERTICES][3];
         int*         doflist=NULL;
-        
-        /* gaussian points: */
-        int     ig;
         GaussTria* gauss=NULL;
-
-        /* Jacobian: */
         double Jdet;
-
-        /*nodal functions: */
         double l1l2l3[3];
-
-        /*element vector at the gaussian points: */
-        double  pe_g[numdof]={0.0};
         double  pe_g_gaussian[numdof];
         double  slope[2];
         int     analysis_type;
 
-        /*inputs :*/
+        /*Initialize Element vector and return if necessary*/
+        if(IsOnWater()) return NULL;
+        ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters);
+
+        /*Retrieve all inputs and parameters*/
+        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
         Input* slope_input=NULL;
-
-        /*retrive parameters: */
-        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-
-        /* Get node coordinates and dof list: */
-        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
-        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-
-        /*Retrieve all inputs we will be needing: */
         if ( (analysis_type==SurfaceSlopeXAnalysisEnum) || (analysis_type==SurfaceSlopeYAnalysisEnum)){
                 slope_input=inputs->GetInput(SurfaceEnum); ISSMASSERT(slope_input);
@@ -4481 +4463 @@
 
                 /*Add pe_g_gaussian vector to pe_g: */
-                for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i];
-
-        }
-
-        /*Add pe_g to global vector pg: */
-        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
+                for(i=0;i<numdof;i++) pe->values[i]+=pe_g_gaussian[i];
+
+        }
 
         /*Clean up and return*/
         delete gauss;
-        xfree((void**)&doflist);
+        return pe;
 }
 /*}}}*/

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

@@ -147 +147 @@
                 void      CreatePVectorAdjointStokes(Vec pg);
                 void      CreatePVectorAdjointBalancedthickness(Vec pg);
-                void      CreatePVectorDiagnosticHutter(Vec pg);
+                ElementVector* CreatePVectorDiagnosticHutter(void);
                 void      CreatePVectorPrognostic(Vec pg);
                 void      CreatePVectorPrognostic_CG(Vec pg);
                 void      CreatePVectorPrognostic_DG(Vec pg);
-                void      CreatePVectorSlope( Vec pg);
+                ElementVector* CreatePVectorSlope(void);
                 void      CreatePVectorThermalSheet( Vec pg);
                 void      CreatePVectorThermalShelf( Vec pg);