Changeset 5929

Timestamp:

09/21/10 16:17:43 (14 years ago)

Author:

Mathieu Morlighem

Message:

more ElementVectors

Location:

issm/trunk/src/c/objects/Elements

Files:

• Penta.cpp (modified) (9 diffs)
• Penta.h (modified) (2 diffs)

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

@@ -753 +753 @@
                         pe=CreatePVectorDiagnosticHoriz();
                         if(pe) pe->AddToGlobal(pg,pf);
-                        delete pe;
+                        if(pe) delete pe;
                         break;
                 case AdjointHorizAnalysisEnum:
                         pe=CreatePVectorAdjointHoriz();
                         if(pe) pe->AddToGlobal(pg,pf);
-                        delete pe;
+                        if(pe) delete pe;
                         break;
                 case DiagnosticHutterAnalysisEnum:
-                        CreatePVectorDiagnosticHutter( pg);
+                        pe=CreatePVectorDiagnosticHutter();
+                        if(pe) pe->AddToGlobal(pg,pf);
+                        if(pe) delete pe;
                         break;
                 case DiagnosticVertAnalysisEnum:
@@ -767 +769 @@
                         break;
                 case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
-                        CreatePVectorSlope( pg);
+                        pe=CreatePVectorSlope();
+                        if(pe) pe->AddToGlobal(pg,pf);
+                        if(pe) delete pe;
                         break;
                 case PrognosticAnalysisEnum:
@@ -3271 +3275 @@
 /*}}}*/
 /*FUNCTION Penta::CreatePVectorDiagnosticHutter{{{1*/
-void  Penta::CreatePVectorDiagnosticHutter(Vec pg){
+ElementVector* Penta::CreatePVectorDiagnosticHutter(void){
 
         int i,j,k;
-        
+        int      ig;
         const int numdofs=NDOF2*NUMVERTICES;
-        int*      doflist=NULL;
-        double    pe_g[numdofs]={0.0};
         double    xyz_list[NUMVERTICES][3];
         double    xyz_list_segment[2][3];
@@ -3285 +3287 @@
         int       node0,node1;
         GaussPenta* gauss=NULL;
-        int      ig;
         double   slope[2];
-
         double   z_g;
         double   rho_ice,gravity,n,B;
@@ -3293 +3293 @@
         double   ub,vb;
         double   surface,thickness;
-
-        /*flags: */
         int  connectivity[2];
 
-        /*If on water, skip: */
-        if(IsOnWater())return;
-
-        /*recover doflist: */
-        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-
-        /*recover parameters: */
+        /*Initialize Element vector and return if necessary*/
+        if(IsOnWater()) return NULL;
+        ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters);
+
+        /*Retrieve all inputs and parameters*/
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
         rho_ice=matpar->GetRhoIce();
         gravity=matpar->GetG();
         n=matice->GetN();
         B=matice->GetB();
-
-        //recover extra inputs
         Input* thickness_input=inputs->GetInput(ThicknessEnum);  ISSMASSERT(thickness_input);
         Input* surface_input=inputs->GetInput(SurfaceEnum);      ISSMASSERT(surface_input);
         Input* slopex_input=inputs->GetInput(SurfaceSlopeXEnum); ISSMASSERT(slopex_input);
         Input* slopey_input=inputs->GetInput(SurfaceSlopeYEnum); ISSMASSERT(slopey_input);
-
-        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
         for(i=0;i<NUMVERTICES;i++)z_list[i]=xyz_list[i][2];
 
@@ -3349 +3342 @@
 
                         if (IsOnSurface()){
-                                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];
+                                for(j=0;j<NDOF2;j++) pe->values[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];
+                                for(j=0;j<NDOF2;j++) pe->values[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss->weight*2/(double)connectivity[1];
                         }
                 }
@@ -3363 +3356 @@
                         vb=constant_part*slope[1];
 
-                        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 */
-        xfree((void**)&doflist);
+                        pe->values[2*node0]+=ub/(double)connectivity[0];
+                        pe->values[2*node0+1]+=vb/(double)connectivity[0];
+                }
+        }
+
+        /*Clean up and return*/
+        return pe;
 }
 /*}}}*/
@@ -3773 +3763 @@
 /*}}}*/
 /*FUNCTION Penta::CreatePVectorSlope {{{1*/
-void Penta::CreatePVectorSlope( Vec pg){
-
-        if (!IsOnBed() || IsOnWater()) return;
+ElementVector* Penta::CreatePVectorSlope(void){
+
+        if (!IsOnBed() || IsOnWater()) return NULL;
 
         /*Call Tria function*/
@@ -3781 +3771 @@
         ElementVector* pe=tria->CreatePVectorSlope();
         delete tria->matice; delete tria;
-        pe->AddToGlobal(pg,NULL);
-        delete pe;
 
         /*clean up and return*/
-        return;
+        return pe;
 }
 /*}}}*/

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

@@ -159 +159 @@
                 ElementVector* CreatePVectorCouplingPattynStokesFriction(void);
                 ElementVector* CreatePVectorDiagnosticHoriz(void);
-                void      CreatePVectorDiagnosticHutter( Vec pg);
+                ElementVector* CreatePVectorDiagnosticHutter(void);
                 ElementVector* CreatePVectorDiagnosticMacAyeal(void);
                 ElementVector* CreatePVectorDiagnosticMacAyealPattyn(void);
@@ -171 +171 @@
                 void      CreatePVectorMelting( Vec pg);
                 void      CreatePVectorPrognostic( Vec pg);
-                void      CreatePVectorSlope( Vec pg);
+                ElementVector* CreatePVectorSlope(void);
                 void      CreatePVectorThermal( Vec pg);
                 void      GetDofList(int** pdoflist,int approximation_enum,int setenum);