Matrix in Linear Reconstruction dynamic array replace by static array (…

…#73)

* Matrix in Linear Reconstruction dynamic array replace by static array

* PR Changes 1

* PR Change 2

src/pmpo_MPMesh.cpp

@@ -34,7 +34,7 @@ void MPMesh::calcBasis() {
             initArray(basisByArea3d,maxVtxsPerElm,0.0);
 
             // calc basis
-            getBasisByAreaGblFormSpherical2(position3d, numVtx, v3d, radius, basisByArea3d);
+            getBasisByAreaGblFormSpherical(position3d, numVtx, v3d, radius, basisByArea3d);
 
             // fill step
             for(int i=0; i<= numVtx; i++){
@@ -339,7 +339,7 @@ void MPMesh::push(){
     else
         p_MPs->rebuild(); //rebuild pumi-pic
     p_MPs->updateMPElmID(); //update mpElm IDs slices
-    reconstructSlices();
+    reconstructSlices(); 
   } 
   while (anyIsMigrating);
 

src/pmpo_MPMesh_assembly.hpp

@@ -147,25 +147,24 @@ void MPMesh::assemblyVtx1() {
   p_MPs->parallel_for(assemble, "assembly");
 
   //Solve Ax=b for each vertex  
-  Kokkos::View<double*[vec4d_nEntries]> VtxCoeffs("VtxMatrices", p_mesh->getNumVertices());
+  Kokkos::View<double*[vec4d_nEntries]> VtxCoeffs("VtxCoeffs", p_mesh->getNumVertices());
   Kokkos::parallel_for("solving Ax=b", numVtx, KOKKOS_LAMBDA(const int vtx){
     Vec4d v0 = {VtxMatrices(vtx,0,0), VtxMatrices(vtx,0,1), VtxMatrices(vtx,0,2), VtxMatrices(vtx,0,3)};
     Vec4d v1 = {VtxMatrices(vtx,1,0), VtxMatrices(vtx,1,1), VtxMatrices(vtx,1,2), VtxMatrices(vtx,1,3)};
     Vec4d v2 = {VtxMatrices(vtx,2,0), VtxMatrices(vtx,2,1), VtxMatrices(vtx,2,2), VtxMatrices(vtx,2,3)};
     Vec4d v3 = {VtxMatrices(vtx,3,0), VtxMatrices(vtx,3,1), VtxMatrices(vtx,3,2), VtxMatrices(vtx,3,3)};
-    Matrix4d A = {v0,v1,v2,v3};
-
-    //double f_norm=A.frobeniusNorm();
+
+    Matrix4d A = {v0,v1,v2,v3}; 
     double A_trace = A.trace();
     Matrix4d A_regularized = {v0, v1, v2, v3};
     A_regularized.addToDiag(A_trace*1e-8);
-
+   
     double coeff[vec4d_nEntries]={0.0, 0.0, 0.0, 0.0};
     CholeskySolve4d_UnitRHS(A_regularized, coeff);
     for (int i=0; i<vec4d_nEntries; i++) 
       VtxCoeffs(vtx,i)=coeff[i];
   });
-
+  
   //Reconstruct
   auto reconstruct = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
     if(mask) { //if material point is 'active'/'enabled'
@@ -189,7 +188,7 @@ void MPMesh::assemblyVtx1() {
     }
   };
   p_MPs->parallel_for(reconstruct, "reconstruct");
-  
+
   Kokkos::parallel_for("assigning", numVtx, KOKKOS_LAMBDA(const int vtx){
     for(int k=0; k<numEntries; k++)
       meshField(vtx, k) = reconVals(vtx,k);

src/pmpo_materialPoints.hpp

@@ -215,7 +215,7 @@ class MaterialPoints {
         auto tgtPosXYZ = MPs->get<MPF_Tgt_Pos_XYZ>();
         auto rotLatLonIncr = MPs->get<MPF_Rot_Lat_Lon_Incr>();
 
-        auto is_rotated = getRotatedFlag();
+        auto is_rotated = getRotatedFlag(); 
         auto updateRotLatLon = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
             if(mask){
                 auto rotLat = curPosRotLatLon(mp,0) + rotLatLonIncr(mp,0); // phi

src/pmpo_utils.hpp

@@ -214,37 +214,21 @@ class Matrix4d {
 
   private:
 
-    double** data_;
-    int rows_=4;
-    int cols_=4;
+    double data_[4][4];
 
   public:
 
     KOKKOS_INLINE_FUNCTION
     Matrix4d(Vec4d v0, Vec4d v1, Vec4d v2, Vec4d v3){
-      data_ = new double*[rows_];
 
-      for (int i=0; i<rows_; i++){
-        data_[i] = new double[cols_];
-      }
-
-      for (int i=0; i<cols_; i++){
+      for (int i=0; i<4; i++){
         data_[0][i] = v0[i];
         data_[1][i] = v1[i];
         data_[2][i] = v2[i];
         data_[3][i] = v3[i];
       }
     }
-
-    //Destructor    
-    KOKKOS_INLINE_FUNCTION
-    ~Matrix4d(){
-      for (int i=0; i<4; i++){
-        delete[] data_[i];
-      }
-      delete[] data_;
-    }
-
+
     //Retrieval
     KOKKOS_INLINE_FUNCTION
     double& operator()(int i, int j) { return data_[i][j];}
@@ -253,8 +237,8 @@ class Matrix4d {
     KOKKOS_INLINE_FUNCTION
     Matrix4d operator/(double scalar) const {
       Matrix4d result(Vec4d(0, 0, 0, 0), Vec4d(0, 0, 0, 0), Vec4d(0, 0, 0, 0), Vec4d(0, 0, 0, 0));
-      for (int i = 0; i < rows_; i++) {
-        for (int j = 0; j < cols_; j++) {
+      for (int i = 0; i < 4; i++) {
+        for (int j = 0; j < 4; j++) {
           result(i, j) = data_[i][j] / scalar;
         }
       }
@@ -265,8 +249,8 @@ class Matrix4d {
     KOKKOS_INLINE_FUNCTION
     double frobeniusNorm() const {
       double sum = 0.0; // Initialize sum of squares
-      for (int i = 0; i < rows_; i++) {
-        for (int j = 0; j < cols_; j++) {
+      for (int i = 0; i < 4; i++) {
+        for (int j = 0; j < 4; j++) {
           sum += data_[i][j]* data_[i][j]; // Sum of squares
         }
       }
@@ -277,7 +261,7 @@ class Matrix4d {
     KOKKOS_INLINE_FUNCTION
     double trace() const {
       double traceSum = 0.0; // Initialize sum of diagonal elements
-      for (int i = 0; i < rows_; i++) {
+      for (int i = 0; i < 4; i++) {
         traceSum += data_[i][i]; // Sum diagonal elements
       }
       return traceSum;
@@ -286,7 +270,7 @@ class Matrix4d {
     //Function to regularize the matrix
     KOKKOS_INLINE_FUNCTION
     void addToDiag(double eps) {
-      for (int i = 1; i < rows_; i++) {
+      for (int i = 1; i < 4; i++) {
         data_[i][i]=data_[i][i]+eps;
       }
     }

src/pmpo_wachspressBasis.hpp

@@ -342,7 +342,7 @@ void sphericalInterpolation(MPMesh& mpMesh){
             initArray(basisByArea3d,maxVtxsPerElm,0.0);
 
             // calc basis
-            getBasisByAreaGblFormSpherical2(position3d, numVtx, v3d, radius, basisByArea3d);
+            getBasisByAreaGblFormSpherical(position3d, numVtx, v3d, radius, basisByArea3d);
 
             // interpolation step
             for(int entry=0; entry<numEntries; entry++){