Changeset 5921

Timestamp:

09/21/10 14:02:29 (15 years ago)

Author:

Mathieu Morlighem

Message:

Some ElementVector, more to come

Location:

issm/trunk/src/c/objects/Elements

Files:

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

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

@@ -2984 +2984 @@
 void  Penta::CreatePVectorAdjointMacAyeal(Vec p_g){
 
-        int i;
-        Tria* tria=NULL;
-
-        /*If on water, skip: */
-        if(IsOnWater()) return;
-
-        /*Bail out if this element if:
-         * -> Non MacAyeal and not on the surface
-         * -> MacAyeal(2d model) and not on bed) */
-        if (!IsOnBed()){
-                return;
-        }
-        else{ 
-
-                /*This element should be collapsed into a tria element at its base. Create this tria element, 
-                 * and compute pe_g*/
-                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
-                tria->CreatePVectorAdjointHoriz(p_g);
-                delete tria->matice; delete tria;
-                return;
-        }
+        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->CreatePVectorAdjointHoriz();
+        delete tria->matice; delete tria;
+        pe->AddToGlobal(p_g,NULL);
+        delete pe;
+
+        /*clean up and return*/
+        return;
 }
 /*}}}*/
@@ -3010 +3000 @@
 void  Penta::CreatePVectorAdjointPattyn(Vec p_g){
 
-        int i;
-        Tria* tria=NULL;
-
-        /*If on water, skip: */
-        if(IsOnWater()) return;
-
-        /*Bail out if this element if:
-         * -> Non MacAyeal and not on the surface
-         * -> MacAyeal(2d model) and not on bed) */
-        if (!IsOnSurface()){
-                return;
-        }
-        else{
-
-                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
-                tria->CreatePVectorAdjointHoriz(p_g);
-                delete tria->matice; delete tria;
-                return;
-        }
+        if (!IsOnSurface() || IsOnWater()) return;
+
+        /*Call Tria function*/
+        Tria* tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
+        ElementVector* pe=tria->CreatePVectorAdjointHoriz();
+        delete tria->matice; delete tria;
+        pe->AddToGlobal(p_g,NULL);
+        delete pe;
+
+        /*clean up and return*/
+        return;
 }
 /*}}}*/
@@ -3371 +3353 @@
 void  Penta::CreatePVectorDiagnosticMacAyeal(Vec pg){
 
-        /*Spawning: */
-        Tria* tria=NULL;
-
-        /*If on water, skip load: */
-        if(IsOnWater())return;
-
-        /*Figure out if this pentaelem is Macayeal. If so, then bailout, except if it is at the 
+        /*Figure out if this penta is collapsed. If so, then bailout, except if it is at the 
           bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build 
-          the load vector. */
-
-        if (!IsOnBed()){
-                /*This element should be collapsed, but this element is not on the bedrock, therefore all its 
-                 * dofs have already been frozen! Do nothing: */
-                return;
-        }
-        else{
-
-                /*This element should be collapsed into a tria element at its base. Create this tria element, 
-                 *and use its CreatePVector functionality to return an elementary load vector: */
-                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-                tria->CreatePVector(pg,NULL);
-                delete tria->matice; delete tria;
-                return;
-        }
+          the stiffness matrix. */
+        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->CreatePVectorDiagnosticMacAyeal();
+        delete tria->matice; delete tria;
+        pe->AddToGlobal(pg,NULL);
+        delete pe;
+
+        /*clean up and return*/
+        return;
 }
 /*}}}*/

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

@@ -751 +751 @@
         switch(analysis_type){
                 case DiagnosticHorizAnalysisEnum:
-                        CreatePVectorDiagnosticMacAyeal( pg,pf);
+                        pe=CreatePVectorDiagnosticMacAyeal();
+                        pe->AddToGlobal(pg,pf);
+                        delete pe;
                         break;
                 case AdjointHorizAnalysisEnum:
-                        CreatePVectorAdjointHoriz( pg);
+                        pe=CreatePVectorAdjointHoriz();
+                        pe->AddToGlobal(pg,pf);
+                        delete pe;
                         break;
                 case DiagnosticHutterAnalysisEnum:
@@ -3666 +3670 @@
 /*}}}*/
 /*FUNCTION Tria::CreatePVectorDiagnosticMacAyeal {{{1*/
-void Tria::CreatePVectorDiagnosticMacAyeal( Vec pg, Vec pf){
+ElementVector* Tria::CreatePVectorDiagnosticMacAyeal(){
 
         /*Constants*/
@@ -3678 +3682 @@
         double         slope[2];
         double         l1l2l3[3];
-        double         pe_g[numdof]={0.0};
         double         pe_g_gaussian[numdof];
         GaussTria*     gauss=NULL;
-        ElementVector* pe=NULL;
-
-        /*First, if we are on water, return empty vector: */
-        if(IsOnWater()) return;
-
-        /*Retrieve all Inputs and 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);
         inputs->GetParameterValue(&drag_type,DragTypeEnum);
         Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); 
         Input* surface_input=inputs->GetInput(SurfaceEnum);     ISSMASSERT(surface_input);
         Input* drag_input=inputs->GetInput(DragCoefficientEnum);ISSMASSERT(drag_input);
-
-        /*Get node coordinates and dof list: */
-        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 
         /* Start  looping on the number of gaussian points: */
@@ -3728 +3729 @@
                         }
                 }
-
-                for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i];
-        }
-
-        pe=NewElementVector(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
-        pe->AddValues(&pe_g[0]);
-        pe->AddToGlobal(pg,pf);
-
-        /*Free ressources:*/
-        delete pe;
+                for(i=0;i<numdof;i++) pe->values[i]+=pe_g_gaussian[i];
+        }
+
+        /*Clean up and return*/
         delete gauss;
-
-
+        return pe;
 }
 /*}}}*/
@@ -3798 +3792 @@
 /*}}}*/
 /*FUNCTION Tria::CreatePVectorAdjointHoriz{{{1*/
-void Tria::CreatePVectorAdjointHoriz(Vec p_g){
-
-        int i;
-
-        /* node data: */
-        const int    numdof=2*NUMVERTICES;
+ElementVector* Tria::CreatePVectorAdjointHoriz(void){
+
+        /*Constants*/
+        const int    numdof=NDOF2*NUMVERTICES;
         double       xyz_list[NUMVERTICES][3];
-        int*         doflist=NULL;
-
-        /* grid data: */
         double vx_list[NUMVERTICES];
         double vy_list[NUMVERTICES];
@@ -3815 +3804 @@
         double duy_list[NUMVERTICES];
         double weights_list[NUMVERTICES];
-
-        /* gaussian points: */
+        int i;
         int     ig;
         GaussTria* gauss=NULL;
-
-        /* parameters: */
         double  obs_velocity_mag,velocity_mag;
         double  dux,duy;
         double  meanvel, epsvel;
-
-        /*element vector : */
-        double  pe_g[numdof]={0.0};
         double  pe_g_gaussian[numdof];
-
-        /* Jacobian: */
         double Jdet;
-
-        /*nodal functions: */
         double l1l2l3[3];
-
-        /*relative and algorithmic fitting: */
         double scalex=0;
         double scaley=0;
@@ -3842 +3819 @@
         int    response;
 
-        /*retrieve some 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);
         this->parameters->FindParam(&meanvel,MeanVelEnum);
         this->parameters->FindParam(&epsvel,EpsVelEnum);
-
-        /* Get node coordinates and dof list: */
-        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
-        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-
-        /* Recover input data: */
         GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum);
         GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum);
@@ -3856 +3832 @@
         GetParameterListOnVertices(&vy_list[0],VyEnum);
         GetParameterListOnVertices(&weights_list[0],WeightsEnum);
-
         inputs->GetParameterValue(&response,CmResponseEnum);
         if(response==SurfaceAverageVelMisfitEnum){
@@ -3986 +3961 @@
                 gauss->GaussPoint(ig);
 
-                /* Get Jacobian determinant: */
                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
-
-                /* Get nodal functions value at gaussian point:*/
                 GetNodalFunctions(l1l2l3, gauss);
 
-                /*Build due_g_gaussian vector: we have three cases here, according to which type of misfit we are using. */
-
-                /*Compute absolute(x/y) at gaussian point: */
                 TriaRef::GetParameterValue(&dux, &dux_list[0],gauss);
                 TriaRef::GetParameterValue(&duy, &duy_list[0],gauss);
 
-                /*compute Du*/
                 for (i=0;i<NUMVERTICES;i++){
                         pe_g_gaussian[i*NDOF2+0]=dux*Jdet*gauss->weight*l1l2l3[i]; 
@@ -4004 +3972 @@
                 }
 
-                /*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 p_g: */
-        VecSetValues(p_g,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

@@ -143 +143 @@
                 void      CreatePVectorBalancedvelocities(Vec pg);
                 void      CreatePVectorDiagnosticBaseVert(Vec pg);
-                void      CreatePVectorDiagnosticMacAyeal(Vec pg,Vec pf);
-                void      CreatePVectorAdjointHoriz(Vec pg);
+                ElementVector* CreatePVectorDiagnosticMacAyeal(void);
+                ElementVector* CreatePVectorAdjointHoriz(void);
                 void      CreatePVectorAdjointStokes(Vec pg);
                 void      CreatePVectorAdjointBalancedthickness(Vec pg);