Changeset 17526

Timestamp:

03/24/14 12:04:51 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: working onXTH

File:

• issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (modified) (5 diffs)

issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

@@ -3387 +3387 @@
 }/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorFSViscousXTH(Element* element){/*{{{*/
-        _error_("STOP");
-
-        int         i,meshtype,dim;
-        IssmDouble  Jdet,forcex,forcey,forcez;
+
+        int         i,tausize,meshtype,dim;
+        IssmDouble  Jdet,r;
         IssmDouble *xyz_list = NULL;
 
@@ -3396 +3395 @@
         element->FindParam(&meshtype,MeshTypeEnum);
         switch(meshtype){
-                case Mesh2DverticalEnum: dim = 2; break;
-                case Mesh3DEnum:         dim = 3; break;
-                case Mesh3DtetrasEnum:   dim = 3; break;
+                case Mesh2DverticalEnum: dim = 2; tausize = 3; break;
+                case Mesh3DEnum:         dim = 3; tausize = 6; break;
+                case Mesh3DtetrasEnum:   dim = 3; tausize = 6; break;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
+
+        /*FIXME*/
+        r = 1.;
 
         /*Fetch number of nodes and dof for this finite element*/
         int vnumnodes = element->NumberofNodesVelocity();
         int pnumnodes = element->NumberofNodesPressure();
+        int tnumnodes = element->GetNumberOfVertices();      //Tensors, P1 DG
 
         /*Prepare coordinate system list*/
@@ -3413 +3416 @@
 
         /*Initialize vectors*/
-        ElementVector* pe     = element->NewElementVector(FSvelocityEnum);
-        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
+        ElementVector* pe    = element->NewElementVector(FSvelocityEnum);
+        IssmDouble*    Dstar = xNew<IssmDouble>(tausize*tnumnodes);
+        IssmDouble*    tau   = xNew<IssmDouble>(tausize);
 
         /*Retrieve all inputs and parameters*/
         element->GetVerticesCoordinates(&xyz_list);
-        IssmDouble  rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
-        IssmDouble  gravity =element->GetMaterialParameter(ConstantsGEnum);
-        Input*      loadingforcex_input=element->GetInput(LoadingforceXEnum);  _assert_(loadingforcex_input);
-        Input*      loadingforcey_input=element->GetInput(LoadingforceYEnum);  _assert_(loadingforcey_input);
-        Input*      loadingforcez_input=NULL;
-        if(dim==3){
-                loadingforcez_input=element->GetInput(LoadingforceZEnum);  _assert_(loadingforcez_input);
-        }
 
         /* Start  looping on the number of gaussian points: */
@@ -3431 +3427 @@
         for(int ig=gauss->begin();ig<gauss->end();ig++){
                 gauss->GaussPoint(ig);
-
                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-                element->NodalFunctionsVelocity(vbasis,gauss);
-
-                loadingforcex_input->GetInputValue(&forcex,gauss);
-                loadingforcey_input->GetInputValue(&forcey,gauss);
-                if(dim==3) loadingforcez_input->GetInputValue(&forcez,gauss);
-
-                for(i=0;i<vnumnodes;i++){
-                        pe->values[i*dim+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i];
-                        pe->values[i*dim+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i];
-                        if(dim==3){
-                                pe->values[i*dim+2] += +rho_ice*forcez*Jdet*gauss->weight*vbasis[i];
-                                pe->values[i*dim+2] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
-                        }
-                        else{
-                                pe->values[i*dim+1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
-                        }
-                }
+                _error_("STOP");
         }
 
@@ -3458 +3437 @@
         delete gauss;
         xDelete<int>(cs_list);
-        xDelete<IssmDouble>(vbasis);
         xDelete<IssmDouble>(xyz_list);
         return pe;