source: issm/oecreview/Archive/25834-26739/ISSM-26222-26223.diff@ 26740

Index: ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp     (revision 26222)
+++ ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp     (revision 26223)
@@ -808,6 +808,11 @@
        delete gauss;
        return Me;
 }/*}}}*/
+ElementVector* DamageEvolutionAnalysis::CreateFctPVector(Element* element){/*{{{*/
+
+       return this->CreatePVector(element);
+
+}/*}}}*/
 void           DamageEvolutionAnalysis::FctKMatrix(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,FemModel* femmodel){/*{{{*/
 
        /*Output*/
@@ -836,6 +841,26 @@
                delete Kfs;
        }
 }/*}}}*/
+void           DamageEvolutionAnalysis::FctPVector(Vector<IssmDouble>** ppf,FemModel* femmodel){/*{{{*/
+
+       /*Output*/
+       Vector<IssmDouble>* pf = NULL;
+
+       /*Initialize P vector*/
+       AllocateSystemMatricesx(NULL,NULL,NULL,&pf,femmodel);
+
+       /*Create and assemble matrix*/
+       for(Object* & object : femmodel->elements->objects){
+               Element*       element = xDynamicCast<Element*>(object);
+               ElementVector* pe      = this->CreateFctPVector(element);
+               if(pe) pe->AddToGlobal(pf);
+               delete pe;
+       }
+       pf->Assemble();
+
+       /*Assign output pointer*/
+       *ppf=pf;
+}/*}}}*/
 void           DamageEvolutionAnalysis::MassMatrix(Matrix<IssmDouble>** pMff,FemModel* femmodel){/*{{{*/
 
        /*Initialize Mass matrix*/
Index: ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h       (revision 26222)
+++ ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h       (revision 26223)
@@ -38,7 +38,9 @@
                /*FCT*/
                ElementMatrix* CreateFctKMatrix(Element* element);
                ElementMatrix* CreateMassMatrix(Element* element);
+               ElementVector* CreateFctPVector(Element* element);
                void           FctKMatrix(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,FemModel* femmodel);
                void           MassMatrix(Matrix<IssmDouble>** pMff,FemModel* femmodel);
+               void                            FctPVector(Vector<IssmDouble>** ppf,FemModel* femmodel);
 };
 #endif
Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp
===================================================================
--- ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp       (revision 26222)
+++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp       (revision 26223)
@@ -417,6 +417,7 @@
                        danalysis = new DamageEvolutionAnalysis();
                        danalysis->MassMatrix(&Mc,femmodel);
                        danalysis->FctKMatrix(&K,NULL,femmodel);
+                       danalysis->FctPVector(&p,femmodel);
                        break;
                default: _error_("analysis type " << EnumToStringx(analysis_type) << " not supported for FCT\n");
        }
Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp        (revision 26222)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp        (revision 26223)
@@ -1815,7 +1815,7 @@
                total_weight+=gauss->weight;
        }
        //normalize to phi. 
-       if(total_weight)for(int i=0;i<NUMVERTICES;i++)weights[i]/=total_weight/phi; 
+       if(total_weight>0.) for(int i=0;i<NUMVERTICES;i++)weights[i]/=total_weight/phi; 
        else for(int i=0;i<NUMVERTICES;i++)weights[i]=0;
 
        /*free ressources:*/
Index: ../trunk-jpl/src/c/classes/SealevelGeometry.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/SealevelGeometry.cpp     (revision 26222)
+++ ../trunk-jpl/src/c/classes/SealevelGeometry.cpp     (revision 26223)
@@ -19,7 +19,7 @@
 SealevelGeometry::SealevelGeometry(int localnelin){ /*{{{*/
        localnel=localnelin;
        for(int i=0;i<SLGEOM_NUMLOADS;i++){
-               for (int j=0;j<MAXVERTICES;j++){
+               for (int j=0;j<SLMAXVERTICES;j++){
                        LoadWeigths[i][j]=xNewZeroInit<IssmDouble>(localnel);
                }
                vlatbarycentre[i]=NULL; //we don't know yet 
@@ -42,7 +42,7 @@
 }; /*}}}*/
 SealevelGeometry::~SealevelGeometry(){ /*{{{*/
        for(int i=0;i<SLGEOM_NUMLOADS;i++){
-               for (int j=0;j<MAXVERTICES;j++){
+               for (int j=0;j<SLMAXVERTICES;j++){
                        xDelete<IssmDouble>(LoadWeigths[i][j]);
                }
                xDelete<IssmDouble>(LoadArea[i]);
Index: ../trunk-jpl/src/c/classes/SealevelGeometry.h
===================================================================
--- ../trunk-jpl/src/c/classes/SealevelGeometry.h       (revision 26222)
+++ ../trunk-jpl/src/c/classes/SealevelGeometry.h       (revision 26223)
@@ -10,7 +10,7 @@
 #define SLGEOM_OCEAN 0 
 #define SLGEOM_ICE 1 
 #define SLGEOM_WATER 2
-#define MAXVERTICES 3
+#define SLMAXVERTICES 3
 
 #include "../toolkits/toolkits.h"
 
@@ -19,7 +19,7 @@
        public: 
 
                int         localnel;
-               IssmDouble* LoadWeigths[SLGEOM_NUMLOADS][MAXVERTICES];
+               IssmDouble* LoadWeigths[SLGEOM_NUMLOADS][SLMAXVERTICES];
                IssmDouble* LoadArea[SLGEOM_NUMLOADS];
                Vector<IssmDouble>* vlatbarycentre[SLGEOM_NUMLOADS];
                Vector<IssmDouble>* vlongbarycentre[SLGEOM_NUMLOADS];