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

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

@@ -808 +808 @@
         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 +841 @@
                 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*/

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

@@ -38 +38 @@
                 /*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

../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp

@@ -417 +417 @@
                         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");
         }

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

@@ -1815 +1815 @@
                 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:*/

../trunk-jpl/src/c/classes/SealevelGeometry.cpp

@@ -19 +19 @@
 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 +42 @@
 }; /*}}}*/
 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]);

../trunk-jpl/src/c/classes/SealevelGeometry.h

@@ -10 +10 @@
 #define SLGEOM_OCEAN 0 
 #define SLGEOM_ICE 1 
 #define SLGEOM_WATER 2
-#define MAXVERTICES 3
+#define SLMAXVERTICES 3
 
 #include "../toolkits/toolkits.h"
 
@@ -19 +19 @@
         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];