[cig-commits] [commit] master: Incomplete fix for EQSim file input, max_slip computation nearly implemented (61ab2bd)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/vc

On branch  : master
Link       : https://github.com/geodynamics/vc/compare/e4325192ad1118379f46ba66899cb98143d09e04...627c49b411aa261c79dbf807e71173eb9e5b1317

>---------------------------------------------------------------

commit 61ab2bd74c5369901499565e5241717933d3dceb
Author: kwschultz <kwschultz at ucdavis.edu>
Date:   Mon Oct 13 16:52:16 2014 -0700

    Incomplete fix for EQSim file input, max_slip computation nearly implemented


>---------------------------------------------------------------

61ab2bd74c5369901499565e5241717933d3dceb
 quakelib/src/QuakeLibIO.cpp | 41 +++++++++++++++++++++++++++++++++++++++--
 1 file changed, 39 insertions(+), 2 deletions(-)

diff --git a/quakelib/src/QuakeLibIO.cpp b/quakelib/src/QuakeLibIO.cpp
index 9e8a2d2..1389e12 100644
--- a/quakelib/src/QuakeLibIO.cpp
+++ b/quakelib/src/QuakeLibIO.cpp
@@ -1574,7 +1574,9 @@ int quakelib::ModelWorld::read_files_eqsim(const std::string &geom_file_name, co
     bool                            read_cond_file;
     quakelib::EQSimGeomSectionMap::const_iterator   sit;
     quakelib::EQSimGeomRectangleMap::const_iterator it;
+    quakelib::eiterator eit;
     quakelib::LatLonDepth           base;
+    std::map<UIndex, double> fault_areas;
 
     // Clear the world first to avoid incorrectly mixing indices
     clear();
@@ -1600,6 +1602,9 @@ int quakelib::ModelWorld::read_files_eqsim(const std::string &geom_file_name, co
     // Take the conversion base as the middle of the section map
     base = quakelib::LatLonDepth(geometry_data.lat_lo(), geometry_data.lon_lo());
     
+    // Initiate converter
+    Conversion          conv(base);
+
     // Triangle elements are currently not supported
     if (geometry_data.num_triangles() > 0) {
         std::cerr << "ERROR: Currently cannot handle EQSim triangle elements. These elements will be ignored." << std::endl;
@@ -1611,8 +1616,8 @@ int quakelib::ModelWorld::read_files_eqsim(const std::string &geom_file_name, co
 
         // Assuming aligned rectangular elements
         for (it=sit->second.rectangles.begin(); it!=sit->second.rectangles.end(); ++it) {
-            quakelib::ModelElement  new_element;
-            unsigned int            i;
+            quakelib::ModelElement   new_element;
+            unsigned int             i;
 
             new_element = it->second.create_model_element();
             new_element.set_section_id(sit->second.sid());
@@ -1624,8 +1629,40 @@ int quakelib::ModelWorld::read_files_eqsim(const std::string &geom_file_name, co
             //new_element.set_static_strength(friction_data.get_static_strength(it->first));
             //new_element.set_dynamic_strength(friction_data.get_dynamic_strength(it->first));
             
+            // Insert partially finished element into the eqsim_world
             eqsim_world.insert(new_element);
+            
+            // Compute area of the current element, add it to the total for this section
+            fault_areas[sit->second.sid()] += eqsim_world.create_sim_element(new_element.id()).area();
+            
+            std::cout << "Computed 1 Area " << std::endl;
+        
+        
         }
+        
+        // Pseudocode:
+        // For each element
+        //      find area of element fault
+        //      assign element max slip based on area
+        //      eqsim_world.element(elem_id).set_max_slip();
+        // ********************** KWS 
+        // Go through the created elements and assign maximum slip based on the fault area for
+        // each element's section
+                
+        for (eit=eqsim_world.begin_element();eit!=eqsim_world.end_element();++eit) {
+        
+            // From Table 2A in Wells Coppersmith 1994
+            double moment_magnitude = 4.07+0.98*log10(conv.sqm2sqkm(fault_areas[eit->get_section_id]));
+             
+            double max_slip = pow(10, (3.0/2.0)*(moment_magnitude+10.7))/(1e7*(eit->lame_mu())*fault_areas[eit->get_section_id]);
+            
+            eit->set_max_slip(max_slip);
+            std::cout << "Max slip: " << max_slip << std::endl;
+            
+            }
+        // ********************** KWS 
+        
+        
     }
 
     insert(eqsim_world);