write P1 element variables in segregated solve (goma#460)

wortiz · Apr 9, 2024 · 4b42444 · 4b42444
1 parent b4f1081
commit 4b42444
Show file tree

Hide file tree

Showing 4 changed files with 65 additions and 3 deletions.
diff --git a/include/wr_soln.h b/include/wr_soln.h
@@ -53,6 +53,7 @@ extern void write_solution_segregated(char output_file[],
                                       double *gv,
                                       struct Results_Description **rd,
                                       double **gvec,
+                                      double ****gvec_elem,
                                       int *nprint,
                                       dbl delta_t,
                                       dbl theta,

diff --git a/src/rf_solve_segregated.c b/src/rf_solve_segregated.c
@@ -1015,7 +1015,8 @@ void solve_problem_segregated(Exo_DB *exo, /* ptr to the finite element mesh dat
           }
 
           write_solution_segregated(ExoFileOut, resid_vector, x, x_old, xdot, xdot_old, tev_post,
-                                    gv, rd, gvec, &nprint, delta_t, theta, time1, NULL, exo, dpi);
+                                    gv, rd, gvec, gvec_elem, &nprint, delta_t, theta, time1, NULL,
+                                    exo, dpi);
 
           if (ProcID == 0) {
             printf("\n Steady state reached \n");
@@ -1346,7 +1347,8 @@ void solve_problem_segregated(Exo_DB *exo, /* ptr to the finite element mesh dat
 
     if (Write_Initial_Solution) {
       write_solution_segregated(ExoFileOut, resid_vector, x, x_old, xdot, xdot_old, tev_post, gv,
-                                rd, gvec, &nprint, delta_t, theta, time1, NULL, exo, dpi);
+                                rd, gvec, gvec_elem, &nprint, delta_t, theta, time1, NULL, exo,
+                                dpi);
       nprint++;
     }
 
@@ -2159,7 +2161,8 @@ void solve_problem_segregated(Exo_DB *exo, /* ptr to the finite element mesh dat
         if (i_print) {
           if (Write_Intermediate_Solutions == 0) {
             write_solution_segregated(ExoFileOut, resid_vector, x, x_old, xdot, xdot_old, tev_post,
-                                      gv, rd, gvec, &nprint, delta_t, theta, time1, NULL, exo, dpi);
+                                      gv, rd, gvec, gvec_elem, &nprint, delta_t, theta, time1, NULL,
+                                      exo, dpi);
             nprint++;
           }
 

diff --git a/src/wr_exo.c b/src/wr_exo.c
@@ -1257,6 +1257,13 @@ void create_truth_table_segregated(struct Results_Description **rd,
               tev += Num_Var_In_Type[pg->imtrx][j];
               exo->truth_table_existance_key[j - V_FIRST] = 1;
             }
+          } else if (pd_glob[mat_num]->i[pg->imtrx][j] == I_P1) {
+            if (exo->truth_table_existance_key[j - V_FIRST] == 0) {
+              /* We just found a candidate for an element variable */
+              int dofs = getdofs(type2shape(exo->eb_elem_itype[eb_indx]), I_P1);
+              tev += dofs;
+              exo->truth_table_existance_key[j - V_FIRST] = dofs;
+            }
           }
         }
       }
@@ -1359,6 +1366,24 @@ void create_truth_table_segregated(struct Results_Description **rd,
                  this test will prevent that. - RRL */
               asdv(&gvec_elem[pg->imtrx][eb_indx][ev_indx - 1], exo->eb_num_elems[eb_indx]);
             }
+          } else if (pd_glob[mat_num]->i[pg->imtrx][j] == I_P1) {
+            int dof = getdofs(type2shape(exo->eb_elem_itype[eb_indx]), I_P1);
+            /* We just found a candidate for an element variable */
+            for (int k = 0; k < dof; k++) {
+              exo->elem_var_tab[i++] = 1;
+              found_match = TRUE;
+              ev_indx++;
+              /* malloc the entry for this block by number of elems for this block
+                 but - only if the variable exists for this block! (by the truth table) */
+
+              if (has_been_called == 0) {
+                /* NOTE: this final array dim is only to be malloc'd once; when a user
+                   is annealing the mesh, anneal mesh calls wr_result_prelim_exo again,
+                   and hence create_truth_table, which would realloc this dim of gvec_elem.
+                   this test will prevent that. - RRL */
+                asdv(&gvec_elem[pg->imtrx][eb_indx][ev_indx - 1], exo->eb_num_elems[eb_indx]);
+              }
+            }
           }
         }
         if (found_match == FALSE && exo->truth_table_existance_key[j - V_FIRST] == 1) {

diff --git a/src/wr_soln.c b/src/wr_soln.c
@@ -162,6 +162,7 @@ void write_solution_segregated(char output_file[],
                                double *gv,
                                struct Results_Description **rd,
                                double **gvec,
+                               double ****gvec_elem,
                                int *nprint,
                                dbl delta_t,
                                dbl theta,
@@ -238,6 +239,38 @@ void write_solution_segregated(char output_file[],
       }
     }
   }
+
+  /* Now element quantities */
+  for (pg->imtrx = 0; pg->imtrx < upd->Total_Num_Matrices; pg->imtrx++) {
+    for (i = 0; i < rd[pg->imtrx]->nev; i++) {
+      bool is_P1 = FALSE;
+      int dof = 0;
+      for (int eb_index = 0; eb_index < exo->num_elem_blocks; eb_index++) {
+        int mn = Matilda[eb_index];
+        if (exo->eb_num_elems[eb_index] > 0) {
+          if (pd_glob[mn]->i[upd->matrix_index[rd[pg->imtrx]->evtype[i]]]
+                            [rd[pg->imtrx]->evtype[i]] == I_P1) {
+            dof = MAX(getdofs(type2shape(exo->eb_elem_itype[eb_index]), I_P1), dof);
+            is_P1 = TRUE;
+          }
+        }
+      }
+      if (is_P1) {
+        for (int k = 0; k < dof; k++) {
+          extract_elem_vec(x[pg->imtrx], i, rd[pg->imtrx]->evtype[i], gvec_elem[pg->imtrx], exo, k);
+          step = (*nprint) + 1;
+          wr_elem_result_exo(exo, output_file, gvec_elem[pg->imtrx], i, step, time_value,
+                             rd[pg->imtrx]);
+          i++;
+        }
+      } else {
+        extract_elem_vec(x[pg->imtrx], i, rd[pg->imtrx]->evtype[i], gvec_elem[pg->imtrx], exo, 0);
+        step = (*nprint) + 1;
+        wr_elem_result_exo(exo, output_file, gvec_elem[pg->imtrx], i, step, time_value,
+                           rd[pg->imtrx]);
+      }
+    }
+  }
   //
   /* Now element quantities */
   //  for(i = 0; i < tev; i++) {