Rounded fault

include/world_builder/features/fault_models/composition/smooth.h

@@ -0,0 +1,96 @@
+/*
+  Copyright (C) 2018 - 2020 by the authors of the World Builder code.
+
+  This file is part of the World Builder.
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published
+   by the Free Software Foundation, either version 2 of the License, or
+   (at your option) any later version.
+
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with this program.  If not, see <https://www.gnu.org/licenses/>.
+*/
+
+#ifndef _world_builder_features_fault_composition_smooth_h
+#define _world_builder_features_fault_composition_smooth_h
+
+#include <world_builder/features/fault_models/composition/interface.h>
+#include <world_builder/features/utilities.h>
+#include <world_builder/world.h>
+
+
+namespace WorldBuilder
+{
+  namespace Features
+  {
+    namespace FaultModels
+    {
+      namespace Composition
+      {
+        /**
+         * This class represents a subducting plate and can implement submodules
+         * for temperature and composition. These submodules determine what
+         * the returned temperature or composition of the temperature and composition
+         * functions of this class will be.
+         */
+        class Smooth: public Interface
+        {
+          public:
+            /**
+             * constructor
+             */
+            Smooth(WorldBuilder::World *world);
+
+            /**
+             * Destructor
+             */
+            ~Smooth();
+
+            /**
+             * declare and read in the world builder file into the parameters class
+             */
+            static
+            void declare_entries(Parameters &prm, const std::string &parent_name = "");
+
+            /**
+             * declare and read in the world builder file into the parameters class
+             */
+            void parse_entries(Parameters &prm) override final;
+
+
+            /**
+             * Returns a composition based on the given position, depth in the model,
+             * and current compostion.
+             */
+            double get_composition(const Point<3> &position,
+                                   const double depth,
+                                   const unsigned int composition_number,
+                                   double composition,
+                                   const double feature_min_depth,
+                                   const double feature_max_depth,
+                                   const WorldBuilder::Utilities::PointDistanceFromCurvedPlanes &distance_from_planes,
+                                   const Utilities::AdditionalParameters &additional_paramters) const override final;
+
+          private:
+            // linear temperature submodule parameters
+            double min_distance;
+            double side_distance;
+            double center_composition;
+            // currenly not using the side composition, but maybe usefu if you want another composition towards the end
+            double side_composition;
+            unsigned int compositions;
+            Utilities::Operations operation;
+
+        };
+      }
+    }
+  }
+}
+
+#endif

include/world_builder/features/utilities.h

@@ -33,7 +33,7 @@ namespace WorldBuilder
     {
       enum class Operations
       {
-        REPLACE,ADD,SUBTRACT
+        REPLACE,ADD,SUBTRACT,MAX
       };
 
       /**
@@ -64,6 +64,9 @@ namespace WorldBuilder
             case Utilities::Operations::SUBTRACT:
               return old_value - new_value;
 
+            case Utilities::Operations::MAX:
+              return std::max(old_value, new_value);
+
             default:
               WBAssert(false,"Operation not found.");
           }

include/world_builder/utilities.h

@@ -412,6 +412,16 @@ namespace WorldBuilder
      */
     std::array<std::array<double,3>,3>
     euler_angles_to_rotation_matrix(double phi1, double theta, double phi2);
+
+    /**
+     * Read a file and distribute the content over all MPI processes.
+     * If WB_WITH_MPI is not defined, this function will just read the file.
+     *
+     * @param filename The name of the file to read.
+     * @return The content of the file.
+    */
+    std::string
+    read_and_distribute_file_content(const std::string &filename);
   } // namespace Utilities
 } // namespace WorldBuilder
 

source/features/fault.cc

@@ -612,13 +612,34 @@ namespace WorldBuilder
                                               section_fraction *
                                               (total_fault_length[next_section] - total_fault_length[current_section]);
 
+              const double rounded_corner = sqrt(abs(0.25*thickness_local*thickness_local - distance_from_plane*distance_from_plane));
 
-              // Because both sides return positve values, we have to
-              // devide the thickness_local by two
+              // We need 3d points in order to conert them from natural to cartesian coordinates
+              const bool bool_cartesian = (this->world->parameters.coordinate_system->natural_coordinate_system() == cartesian);
+
+              const Point<3> P1 (bool_cartesian ? coordinates[current_section][0] : starting_radius,
+                                 bool_cartesian ? coordinates[current_section][1] : coordinates[current_section][0],
+                                 bool_cartesian ? starting_radius : coordinates[current_section][1],
+                                 this->world->parameters.coordinate_system->natural_coordinate_system());
+
+              const Point<3> P2 (bool_cartesian ? coordinates[next_section][0] : starting_radius,
+                                 bool_cartesian ? coordinates[next_section][1] : coordinates[current_section][0],
+                                 bool_cartesian ? starting_radius : coordinates[next_section][1],
+                                 this->world->parameters.coordinate_system->natural_coordinate_system());
+
+              Point<3> P1_cartesian(this->world->parameters.coordinate_system->natural_to_cartesian_coordinates(P1.get_array()),cartesian);
+              Point<3> P2_cartesian(this->world->parameters.coordinate_system->natural_to_cartesian_coordinates(P2.get_array()),cartesian);
+
+              const double fault_length_at_surface = (P2_cartesian - P1_cartesian).norm();
+
+              // Because both sides return positive values, we have to
+              // divide the thickness_local by two
               if (std::fabs(distance_from_plane) > 0 &&
                   std::fabs(distance_from_plane) <= thickness_local * 0.5 &&
                   distance_along_plane > 0 &&
-                  distance_along_plane <= max_fault_length)
+                  distance_along_plane <= max_fault_length &&
+                  section_fraction * fault_length_at_surface > - rounded_corner &&
+                  (section_fraction - 1) * fault_length_at_surface < rounded_corner)
                 {
                   // Inside the fault!
                   const Features::Utilities::AdditionalParameters additional_parameters = {max_fault_length,thickness_local};

source/features/fault_models/composition/smooth.cc

@@ -0,0 +1,124 @@
+/*
+  Copyright (C) 2018 - 2020 by the authors of the World Builder code.
+
+  This file is part of the World Builder.
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published
+   by the Free Software Foundation, either version 2 of the License, or
+   (at your option) any later version.
+
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with this program.  If not, see <https://www.gnu.org/licenses/>.
+*/
+
+#include <world_builder/utilities.h>
+#include <world_builder/assert.h>
+#include <world_builder/nan.h>
+#include <world_builder/parameters.h>
+
+#include <world_builder/types/double.h>
+#include <world_builder/types/string.h>
+#include <world_builder/types/object.h>
+#include <world_builder/features/fault_models/composition/smooth.h>
+
+
+namespace WorldBuilder
+{
+  using namespace Utilities;
+
+  namespace Features
+  {
+    namespace FaultModels
+    {
+      namespace Composition
+      {
+        Smooth::Smooth(WorldBuilder::World *world_)
+          :
+          min_distance(NaN::DSNAN),
+          side_distance(NaN::DSNAN),
+          center_composition(NaN::DSNAN),
+          side_composition(NaN::DSNAN),
+          operation(Utilities::Operations::MAX)
+        {
+          this->world = world_;
+          this->name = "smooth";
+        }
+
+        Smooth::~Smooth()
+        { }
+
+        void
+        Smooth::declare_entries(Parameters &prm, const std::string &)
+        {
+          // Add compositions to the required parameters.
+          prm.declare_entry("", Types::Object({"compositions"}), "Compositional model object");
+
+          prm.declare_entry("min distance fault center", Types::Double(0),
+                            "The distance in meters from which the composition of this feature is present.");
+
+          prm.declare_entry("side distance fault center", Types::Double(std::numeric_limits<double>::max()),
+                            "The distance over which the composition is reduced from 1 to 0.");
+
+          prm.declare_entry("center composition", Types::Double(1.),
+                            "The composition at the center of the fault.");
+
+          prm.declare_entry("side composition", Types::Double(0),
+                            "The composition at the sides of this feature.");
+
+          prm.declare_entry("operation", Types::String("max", std::vector<std::string> {"max"}),
+                            "Whether the value should replace any value previously defined at this location (replace) or "
+                            "add the value to the previously define value (add, not implemented). Replacing implies that all values not "
+                            "explicitly defined are set to zero.");
+        }
+
+        void
+        Smooth::parse_entries(Parameters &prm)
+        {
+          min_distance = prm.get<double>("min distance fault center");
+          side_distance = prm.get<double>("side distance fault center");
+          WBAssert(side_distance >= min_distance, "distance at the side needs to be larger or equal than the min distance.");
+          operation = Utilities::string_operations_to_enum(prm.get<std::string>("operation"));
+          center_composition = prm.get<double>("center composition");
+          side_composition = prm.get<double>("side composition");
+          compositions = prm.get<unsigned int>("compositions");
+        }
+
+
+        double
+        Smooth::get_composition( const Point<3> &,
+                                 const double ,
+                                 const unsigned int composition_number,
+                                 double composition_,
+                                 const double ,
+                                 const double ,
+                                 const WorldBuilder::Utilities::PointDistanceFromCurvedPlanes &distance_from_planes,
+                                 const Utilities::AdditionalParameters &additional_paramters) const 
+	{
+          double composition = 0;
+
+          if (compositions == composition_number)
+            {
+
+              // Hyperbolic tangent goes from 0 to 1 over approximately x=(0, 2) without any arguements. The function is written
+              // so that the composition returned 1 to 0 over the side_distance on either sides.
+              composition = (center_composition - std::tanh(10 * (distance_from_planes.distance_from_plane - side_distance/2)/side_distance ) )/2 ;
+
+              return Utilities::apply_operation(operation,composition_,composition);
+            }
+
+          return composition;
+
+        }
+
+        WB_REGISTER_FEATURE_FAULT_COMPOSITION_MODEL (Smooth, smooth)
+      }
+    }
+  }
+}
+

source/features/utilities.cc

@@ -31,6 +31,7 @@ namespace WorldBuilder
       {
         if (operation == "add") return Operations::ADD;
         if (operation == "subtract") return Operations::SUBTRACT;
+        if (operation == "max") return Operations::MAX;
 
         WBAssert(operation == "replace", "Could not find operation: " << operation << ".");
         return Operations::REPLACE;

source/parameters.cc

@@ -80,16 +80,9 @@ namespace WorldBuilder
       }
 
     path_level =0;
-    // Now read in the world builder file into a file stream and
-    // put it into a the rapidjason document
-    std::ifstream json_input_stream(filename.c_str());
-
-    // Get world builder file and check whether it exists
-    WBAssertThrow(json_input_stream.good(),
-                  "Could not find the world builder file at the specified location: " + filename);
-
-    WBAssert(json_input_stream, "Could not read the world builder file.");
-
+    // Now read in the world builder file into a stringstream and
+    // put it into a the rapidjson document
+    std::stringstream json_input_stream(WorldBuilder::Utilities::read_and_distribute_file_content(filename));
     rapidjson::IStreamWrapper isw(json_input_stream);
 
     // relaxing sytax by allowing comments () for now, maybe also allow trailing commas and (kParseTrailingCommasFlag) and nan's, inf etc (kParseNanAndInfFlag)?
@@ -120,7 +113,6 @@ namespace WorldBuilder
                                                                             ));
 
     WBAssertThrow(parameters.IsObject(), "World builder file is is not an object.");
-    json_input_stream.close();
 
 
     SchemaDocument schema(declarations);