source: issm/oecreview/Archive/16554-17801/ISSM-17000-17001.diff

../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h

@@ -61 +61 @@
                 void GetBHOFriction(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
                 void InputUpdateFromSolutionHO(IssmDouble* solution,Element* element);
                 /*FS*/
+                ElementVector* CreateDVectorFS(Element* element);
                 ElementMatrix* CreateJacobianMatrixFS(Element* element);
                 ElementMatrix* CreateKMatrixFS(Element* element);
                 ElementMatrix* CreateKMatrixFSViscous(Element* element);

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

@@ -812 +812 @@
 
 /*Finite Element Analysis*/
 ElementVector* StressbalanceAnalysis::CreateDVector(Element* element){/*{{{*/
-        /*Default, return NULL*/
+
+        int approximation;
+        element->GetInputValue(&approximation,ApproximationEnum);
+        switch(approximation){
+                case FSApproximationEnum: 
+                        return CreateDVectorFS(element);
+                default:
+                        return NULL; //no need for doftypes outside of FS approximation
+        }
         return NULL;
+
 }/*}}}*/
 ElementMatrix* StressbalanceAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
 
@@ -2380 +2389 @@
 }/*}}}*/
 
 /*FS*/
+ElementVector* StressbalanceAnalysis::CreateDVectorFS(Element* element){/*{{{*/
+
+        int         meshtype,dim;
+
+        /*Get problem dimension*/
+        element->FindParam(&meshtype,MeshTypeEnum);
+        switch(meshtype){
+                case Mesh2DverticalEnum: dim = 2; break;
+                case Mesh3DEnum:         dim = 3; break;
+                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+        }
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int vnumnodes = element->GetNumberOfNodesVelocity();
+        int pnumnodes = element->GetNumberOfNodesPressure();
+
+        /*Initialize output vector*/
+        ElementVector* de = element->NewElementVector(FSvelocityEnum);
+
+        for(int i=0;i<vnumnodes;i++){
+                for(int j=0;j<dim;j++) de->values[i*dim+j]=VelocityEnum;
+        }
+        for(int i=0;i<pnumnodes;i++){
+                de->values[vnumnodes*dim+i]=PressureEnum;
+        }
+
+        return de;
+
+}/*}}}*/
 ElementMatrix* StressbalanceAnalysis::CreateJacobianMatrixFS(Element* element){/*{{{*/
 
         /*Intermediaries */

../trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp

@@ -101 +101 @@
         }
 
         /*Create dof vector for stiffness matrix preconditioning*/
-        for (i=0;i<femmodel->elements->Size();i++){
-                element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
-                element->CreateDVector(df);
+        if(pdf){
+                for(i=0;i<femmodel->elements->Size();i++){
+                        element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                        ElementVector* de=analysis->CreateDVector(element);
+                        if(de) de->InsertIntoGlobal(df);
+                        delete de;
+                }
         }
 
         /*Assemble matrices and vector*/

../trunk-jpl/src/c/classes/Elements/Element.h

@@ -103 +103 @@
                 virtual void       Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0;
                 virtual void       SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters)=0;
                 virtual void       SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum)=0;
-                virtual void   CreateDVector(Vector<IssmDouble>* df)=0;
                 virtual void   ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0;
                 virtual void   EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure)=0;
                 virtual IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure)=0;

../trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -190 +190 @@
         return sqrt(2*this->GetArea());
 }
 /*}}}*/
-/*FUNCTION Tria::CreateDVector {{{*/
-void  Tria::CreateDVector(Vector<IssmDouble>* df){
-
-        /*Nothing done yet*/
-}
-/*}}}*/
 /*FUNCTION Tria::ComputeBasalStress {{{*/
 void  Tria::ComputeBasalStress(Vector<IssmDouble>* eps){
         _error_("Not Implemented yet");

../trunk-jpl/src/c/classes/Elements/Tria.h

@@ -66 +66 @@
                 void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
                 void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
                 void        SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
-                void        CreateDVector(Vector<IssmDouble>* df);
                 void        Delta18oParameterization(void);
                 void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");};
                 void        EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");};

../trunk-jpl/src/c/classes/Elements/Penta.cpp

@@ -378 +378 @@
         this->inputs->Configure(parameters);
 }
 /*}}}*/
-/*FUNCTION Penta::CreateDVector {{{*/
-void  Penta::CreateDVector(Vector<IssmDouble>* df){
-
-        /*retrieve parameters: */
-        ElementVector* De=NULL;
-        int analysis_type;
-        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-
-        /*Checks in debugging*/
-        _assert_(this->nodes && this->material && this->matpar && this->verticalneighbors && this->parameters && this->inputs);
-
-        switch(analysis_type){
-                #ifdef _HAVE_STRESSBALANCE_
-                case StressbalanceAnalysisEnum:
-                        De=CreateDVectorStressbalanceHoriz();
-                        break;
-                #endif
-        }
-
-        /*Add to global Vector*/
-        if(De){
-                De->InsertIntoGlobal(df);
-                delete De;
-        }
-}
-/*}}}*/
 /*FUNCTION Penta::DeepEcho{{{*/
 void Penta::DeepEcho(void){
 
@@ -4475 +4449 @@
 /*}}}*/
 #endif
 
-#ifdef _HAVE_STRESSBALANCE_
-/*FUNCTION Penta::CreateDVectorStressbalanceHoriz {{{*/
-ElementVector* Penta::CreateDVectorStressbalanceHoriz(void){
-
-        int approximation;
-        inputs->GetInputValue(&approximation,ApproximationEnum);
-
-        switch(approximation){
-                case FSApproximationEnum:
-                        return CreateDVectorStressbalanceFS();
-                default:
-                        return NULL; //no need for doftypes outside of FS approximation
-        }
-}
-/*}}}*/
-/*FUNCTION Penta::CreateDVectorStressbalanceFS{{{*/
-ElementVector* Penta::CreateDVectorStressbalanceFS(void){
-
-        /*output: */
-        ElementVector* De=NULL;
-
-        /*Initialize Element vector and return if necessary*/
-        int approximation;
-        inputs->GetInputValue(&approximation,ApproximationEnum);
-        if(approximation!=FSApproximationEnum) return NULL;
-
-        /*Fetch number of nodes and dof for this finite element*/
-        int vnumnodes = this->NumberofNodesVelocity();
-        int pnumnodes = this->NumberofNodesPressure();
-
-        De=new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
-
-        for(int i=0;i<vnumnodes;i++){
-                De->values[i*3+0]=VelocityEnum;
-                De->values[i*3+1]=VelocityEnum;
-                De->values[i*3+2]=VelocityEnum;
-        }
-        for(int i=0;i<pnumnodes;i++){
-                De->values[vnumnodes*3+i]=PressureEnum;
-        }
-
-        return De;
-}
-/*}}}*/
-#endif
-
 #ifdef _HAVE_HYDROLOGY_
 /*FUNCTION Penta::CreateEPLDomainMassMatrix {{{*/
 ElementMatrix* Penta::CreateEPLDomainMassMatrix(void){

../trunk-jpl/src/c/classes/Elements/Penta.h

@@ -70 +70 @@
                 int    FiniteElement(void);
                 void   SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
                 void   SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
-                void   CreateDVector(Vector<IssmDouble>* df);
                 void   Delta18oParameterization(void);
                 void   EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure);
                 Element* GetBasalElement(void);
@@ -251 +250 @@
                 Tria*            SpawnTria(int location);
                 IssmDouble     StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa);
 
-                #ifdef _HAVE_STRESSBALANCE_
-                ElementVector* CreateDVectorStressbalanceHoriz(void);
-                ElementVector* CreateDVectorStressbalanceFS(void);
-                #endif
-
                 #ifdef _HAVE_HYDROLOGY_
                 ElementMatrix* CreateEPLDomainMassMatrix(void);
                 void           GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);

../trunk-jpl/src/c/classes/Elements/Seg.h

@@ -68 +68 @@
                 void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters){_error_("not implemented yet");};
                 void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters){_error_("not implemented yet");};
                 void        SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){_error_("not implemented yet");};
-                void        CreateDVector(Vector<IssmDouble>* df){_error_("not implemented yet");};
                 void        Delta18oParameterization(void){_error_("not implemented yet");};
                 void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");};
                 void        EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");};