From dad0466d6db3b5f3a95ab91d56c996d0998e2dd8 Mon Sep 17 00:00:00 2001
From: Haoyuan Li <hylli@ucdavis.edu>
Date: Mon, 4 Nov 2024 20:16:16 -0800
Subject: [PATCH] follow

---
 include/aspect/material_model/utilities.h |  7 ++++++
 include/aspect/structured_data.h          |  3 +++
 source/material_model/utilities.cc        | 29 +++++++++++++++++++----
 source/structured_data.cc                 |  8 +++++++
 4 files changed, 42 insertions(+), 5 deletions(-)

diff --git a/include/aspect/material_model/utilities.h b/include/aspect/material_model/utilities.h
index aa05a0b7bf7..1d401dc62b7 100644
--- a/include/aspect/material_model/utilities.h
+++ b/include/aspect/material_model/utilities.h
@@ -653,6 +653,13 @@ namespace aspect
           */
           std::string data_directory;
           std::vector<std::string> material_file_names;
+          std::vector<double> minimum_temperature;
+          std::vector<double> maximum_temperature;
+          std::vector<double> interval_temperature;
+          std::vector<double> minimum_pressure;
+          std::vector<double> maximum_pressure;
+          std::vector<double> interval_pressure;
+
 
           /**
            * List of pointers to objects that read and process data we get from
diff --git a/include/aspect/structured_data.h b/include/aspect/structured_data.h
index 121e6d821aa..72a1da884f1 100644
--- a/include/aspect/structured_data.h
+++ b/include/aspect/structured_data.h
@@ -248,6 +248,9 @@ namespace aspect
          */
         double get_maximum_component_value(const unsigned int component) const;
 
+        //todo_PT
+        unsigned get_table_points(const unsigned int dimension) const;
+
       private:
         /**
          * The number of data components read in (=columns in the data file).
diff --git a/source/material_model/utilities.cc b/source/material_model/utilities.cc
index 551c8a6f057..651f86f848f 100644
--- a/source/material_model/utilities.cc
+++ b/source/material_model/utilities.cc
@@ -1167,6 +1167,12 @@ namespace aspect
                      ExcMessage(" The 'discrete phase function' requires that the number of compositional fields "
                                 "matches the number of lookup tables."));
 
+        minimum_temperature = std::vector<double>(material_file_names.size());
+        maximum_temperature = std::vector<double>(material_file_names.size());
+        interval_temperature = std::vector<double>(material_file_names.size());
+        minimum_pressure = std::vector<double>(material_file_names.size());
+        maximum_pressure = std::vector<double>(material_file_names.size());
+        interval_pressure = std::vector<double>(material_file_names.size());
         for (unsigned i = 0; i < material_file_names.size(); ++i)
           {
             material_lookup
@@ -1174,8 +1180,14 @@ namespace aspect
 
             material_lookup[i]->load_file(data_directory+material_file_names[i],
                                           this->get_mpi_communicator());
-          }
 
+            minimum_temperature[i] = material_lookup[i]->get_interpolation_point_coordinates(0).front();
+            maximum_temperature[i] = material_lookup[i]->get_interpolation_point_coordinates(0).back();
+            interval_temperature[i] = (maximum_temperature[i]-minimum_temperature[i])/(material_lookup[i]->get_table_points(0)-1);
+            minimum_pressure[i] = material_lookup[i]->get_interpolation_point_coordinates(1).front();
+            maximum_pressure[i] = material_lookup[i]->get_interpolation_point_coordinates(1).back();
+            interval_pressure[i] = (maximum_pressure[i]-minimum_pressure[i])/(material_lookup[i]->get_table_points(1)-1);
+          }
       }
 
 
@@ -1189,12 +1201,19 @@ namespace aspect
 
         // lookup the most abundant phase
 
-        // Convert pressure from Pa to bar, bar is used in the table.
-        Point<2> temperature_pressure(in.temperature, in.pressure / 1.e5);
+        const unsigned i = 0;
+
+        Assert (in.temperature >= minimum_temperature[i] && in.temperature < maximum_temperature[i], ExcInternalError());
+        Assert (in.pressure/1e5 >= minimum_pressure[i] && in.pressure/1e5 < maximum_pressure[i], ExcInternalError());
+
+        // const unsigned maph = unsigned(material_lookup[0]->get_data(temperature_pressure, 7));
+        const std::vector<double>& maph_points = material_lookup[0]->get_interpolation_point_coordinates(7);
 
-        const unsigned maph = unsigned(material_lookup[0]->get_data(temperature_pressure, 7));
+        const unsigned i_T = static_cast<unsigned>(std::round((in.temperature - minimum_temperature[i]) / interval_temperature[i]));
+        const unsigned i_p = static_cast<unsigned>(std::round((in.pressure/1e5 - minimum_pressure[i]) / interval_pressure[i]));
+        const unsigned i_point = i_p * material_lookup[0]->get_table_points(0) + i_T;
 
-        // Convert pressure from Pa to bar, bar is used in the table.
+        const unsigned maph = unsigned(maph_points[i_point]);
 
         // determine the value
         if (in.phase_index < maph)
diff --git a/source/structured_data.cc b/source/structured_data.cc
index 62da14cd9b8..40b0e290f5d 100644
--- a/source/structured_data.cc
+++ b/source/structured_data.cc
@@ -144,6 +144,14 @@ namespace aspect
     {
       return maximum_component_value[component];
     }
+        
+    //todo_PT
+    template <int dim>
+    unsigned 
+    StructuredDataLookup<dim>::get_table_points(const unsigned int dimension) const
+    {
+      return table_points[dimension];
+    }