Changeset 14435

Timestamp:

03/05/13 19:37:53 (12 years ago)

Author:

Mathieu Morlighem

Message:

NEW: CreateDVector is now independent of CreateKMatrix. It is more consistent

Location:

issm/trunk-jpl/src/c/classes

Files:

• FemModel.cpp (modified) (3 diffs)
• matrix/ElementVector.cpp (modified) (1 diff)
• objects/Elements/Element.h (modified) (1 diff)
• objects/Elements/Penta.cpp (modified) (4 diffs)
• objects/Elements/Penta.h (modified) (1 diff)
• objects/Elements/Tria.cpp (modified) (2 diffs)
• objects/Elements/Tria.h (modified) (1 diff)

issm/trunk-jpl/src/c/classes/FemModel.cpp

@@ -826 +826 @@
                 for (i=0;i<this->elements->Size();i++){
                         element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
-                        element->CreateKMatrix(Kff_temp,NULL,NULL);
+                        element->CreateKMatrix(Kff_temp,NULL);
                 }
                 for (i=0;i<this->loads->Size();i++){
@@ -845 +845 @@
         for (i=0;i<this->elements->Size();i++){
                 element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
-                element->CreateKMatrix(Kff,Kfs,df);
+                element->CreateKMatrix(Kff,Kfs);
         }
         for (i=0;i<this->loads->Size();i++){
@@ -872 +872 @@
                         if(load->InAnalysis(configuration_type)) load->PenaltyCreatePVector(pf,kmax);
                 }
+        }
+
+        /*Create dof vector for stiffness matrix preconditioning*/
+        for (i=0;i<this->elements->Size();i++){
+                element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+                element->CreateDVector(df);
         }
 

issm/trunk-jpl/src/c/classes/matrix/ElementVector.cpp

@@ -190 +190 @@
         IssmDouble* localvalues=NULL;
 
-        if(this->fsize && pf){
+        if(this->fsize){
                 /*first, retrieve values that are in the f-set from the g-set values vector: */
                 localvalues=xNew<IssmDouble>(this->fsize);

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

@@ -30 +30 @@
                 virtual void   SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters)=0;
                 virtual void   SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int set1_enum,int set2_enum)=0;
-                virtual void   CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>*  Kfs,Vector<IssmDouble>* df)=0;
+                virtual void   CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>*  Kfs)=0;
+                virtual void   CreateDVector(Vector<IssmDouble>* df)=0;
                 virtual void   CreatePVector(Vector<IssmDouble>* pf)=0;
                 virtual void   CreateJacobianMatrix(Matrix<IssmDouble>* Jff)=0;

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

@@ -405 +405 @@
 /*}}}*/
 /*FUNCTION Penta::CreateKMatrix {{{*/
-void  Penta::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,Vector<IssmDouble>* df){
+void  Penta::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){
 
         /*retrieve parameters: */
         ElementMatrix* Ke=NULL;
-        ElementVector* De=NULL;
         int analysis_type;
         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
@@ -424 +423 @@
                 #ifdef _HAVE_DIAGNOSTIC_
                 case DiagnosticHorizAnalysisEnum:
-                        Ke=CreateKMatrixDiagnosticHoriz(); De=CreateDVectorDiagnosticHoriz();
+                        Ke=CreateKMatrixDiagnosticHoriz();
                         break;
                 case AdjointHorizAnalysisEnum:
@@ -467 +466 @@
                 delete Ke;
         }
+}
+/*}}}*/
+/*FUNCTION Penta::CreateKMatrixPrognostic {{{*/
+ElementMatrix* Penta::CreateKMatrixPrognostic(void){
+
+        if (!IsOnBed()) return NULL;
+
+        /*Depth Averaging Vx and Vy*/
+        this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
+        this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
+
+        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
+        ElementMatrix* Ke=tria->CreateKMatrixPrognostic();
+        delete tria->material; delete tria;
+
+        /*Delete Vx and Vy averaged*/
+        this->inputs->DeleteInput(VxAverageEnum);
+        this->inputs->DeleteInput(VyAverageEnum);
+
+        /*clean up and return*/
+        return Ke;
+}
+/*}}}*/
+/*FUNCTION Penta::CreateBasalMassMatrix{{{*/
+ElementMatrix* Penta::CreateBasalMassMatrix(void){
+
+        if (!IsOnBed()) return NULL;
+
+        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
+        ElementMatrix* Ke=tria->CreateMassMatrix();
+        delete tria->material; delete tria;
+
+        /*clean up and return*/
+        return Ke;
+}
+/*}}}*/
+/*FUNCTION Penta::CreateDVector {{{*/
+void  Penta::CreateDVector(Vector<IssmDouble>* df){
+
+        /*retrieve parameters: */
+        ElementMatrix* Ke=NULL;
+        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_DIAGNOSTIC_
+                case DiagnosticHorizAnalysisEnum:
+                        De=CreateDVectorDiagnosticHoriz();
+                        break;
+                #endif
+        }
+
         /*Add to global Vector*/
         if(De){
@@ -472 +528 @@
                 delete De;
         }
-}
-/*}}}*/
-/*FUNCTION Penta::CreateKMatrixPrognostic {{{*/
-ElementMatrix* Penta::CreateKMatrixPrognostic(void){
-
-        if (!IsOnBed()) return NULL;
-
-        /*Depth Averaging Vx and Vy*/
-        this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
-        this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
-
-        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
-        ElementMatrix* Ke=tria->CreateKMatrixPrognostic();
-        delete tria->material; delete tria;
-
-        /*Delete Vx and Vy averaged*/
-        this->inputs->DeleteInput(VxAverageEnum);
-        this->inputs->DeleteInput(VyAverageEnum);
-
-        /*clean up and return*/
-        return Ke;
-}
-/*}}}*/
-/*FUNCTION Penta::CreateBasalMassMatrix{{{*/
-ElementMatrix* Penta::CreateBasalMassMatrix(void){
-
-        if (!IsOnBed()) return NULL;
-
-        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
-        ElementMatrix* Ke=tria->CreateMassMatrix();
-        delete tria->material; delete tria;
-
-        /*clean up and return*/
-        return Ke;
 }
 /*}}}*/

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

@@ -81 +81 @@
                 void   SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters);
                 void   SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int set1_enum,int set2_enum);
-                void   CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,Vector<IssmDouble>* df);
+                void   CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs);
+                void   CreateDVector(Vector<IssmDouble>* df);
                 void   CreatePVector(Vector<IssmDouble>* pf);
                 void   CreateJacobianMatrix(Matrix<IssmDouble>* Jff);

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

@@ -214 +214 @@
 /*}}}*/
 /*FUNCTION Tria::CreateKMatrix {{{*/
-void  Tria::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,Vector<IssmDouble>* df){
+void  Tria::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){
 
         /*retreive parameters: */
@@ -561 +561 @@
         delete gauss;
         return Ke;
+}
+/*}}}*/
+/*FUNCTION Tria::CreateDVector {{{*/
+void  Tria::CreateDVector(Vector<IssmDouble>* df){
+
+        /*Nothing done yet*/
 }
 /*}}}*/

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

@@ -77 +77 @@
                 void   SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters);
                 void   SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int set1_enum,int set2_enum);
-                void   CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,Vector<IssmDouble>* df);
+                void   CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs);
+                void   CreateDVector(Vector<IssmDouble>* df);
                 void   CreatePVector(Vector<IssmDouble>* pf);
                 void   CreateJacobianMatrix(Matrix<IssmDouble>* Jff);