[GeoMechanicsApplication] Phreatic multi line

applications/GeoMechanicsApplication/custom_processes/apply_constant_phreatic_multi_line_pressure_process.cpp

@@ -0,0 +1,284 @@
+// KRATOS___
+// // ) )
+// // ___ ___
+// // ____ //___) ) // ) )
+// // / / // // / /
+// ((____/ / ((____ ((___/ / MECHANICS
+//
+// License: geo_mechanics_application/license.txt
+//
+// Main authors: Vahid Galavi
+// Jonathan Nuttall
+
+#include "apply_constant_phreatic_multi_line_pressure_process.h"
+
+namespace Kratos
+{
+
+ApplyConstantPhreaticMultiLinePressureProcess::ApplyConstantPhreaticMultiLinePressureProcess(ModelPart& model_part,
+ Parameters rParameters)
+ : Process(Flags()) , mrModelPart(model_part)
+{
+ KRATOS_TRY
+
+ InitializeParameters(rParameters);
+
+ mVariableName = rParameters["variable_name"].GetString();
+ mIsFixed = rParameters["is_fixed"].GetBool();
+ mIsFixedProvided = rParameters.Has("is_fixed");
+ mIsSeepage = rParameters["is_seepage"].GetBool();
+ mSpecificWeight = rParameters["specific_weight"].GetDouble();
+ mPressureTensionCutOff = rParameters["pressure_tension_cut_off"].GetDouble();
+ mOutOfPlaneDirection = rParameters["out_of_plane_direction"].GetInt();
+
+ InitializeCoordinates(rParameters);
+ InitializeGravityDirection(rParameters);
+ InitializeHorizontalDirection();
+
+ ValidateCoordinates(rParameters);
+
+ KRATOS_CATCH("")
+}
+
+void ApplyConstantPhreaticMultiLinePressureProcess::InitializeParameters(Parameters& rParameters) const
+{
+ Parameters default_parameters(R"(
+ {
+ "model_part_name":"PLEASE_CHOOSE_MODEL_PART_NAME",
+ "variable_name": "PLEASE_PRESCRIBE_VARIABLE_NAME",
+ "is_fixed": false,
+ "is_seepage": false,
+ "gravity_direction": 1,
+ "out_of_plane_direction": 2,
+ "x_coordinates": [0.0,1.0],
+ "y_coordinates": [1.0,0.5],
+ "z_coordinates": [0.0,0.0],
+ "specific_weight" : 10000.0,
+ "pressure_tension_cut_off" : 0.0,
+ "table" : [0,1]
+ } )" );
+
+ // Some values are mandatory, since no meaningful default value exist. For this reason try accessing to them
+ // So that an error is thrown if they don't exist
+ rParameters["x_coordinates"];
+ rParameters["y_coordinates"];
+ rParameters["z_coordinates"];
+ rParameters["variable_name"];
+ rParameters["model_part_name"];
+
+ // Now validate the defaults again -- this also ensures no type mismatch
+ rParameters.ValidateAndAssignDefaults(default_parameters);
+}
+
+void ApplyConstantPhreaticMultiLinePressureProcess::ValidateCoordinates(const Parameters &rParameters) const
+{
+ if (GravityDirection() == OutOfPlaneDirection())
+ {
+ KRATOS_ERROR << "Gravity direction cannot be the same as Out-of-Plane directions"
+ << rParameters
+ << std::endl;
+ }
+
+ if (!std::is_sorted(HorizontalDirectionCoordinates().begin(), HorizontalDirectionCoordinates().end()))
+ {
+ KRATOS_ERROR << "The Horizontal Elements Coordinates are not ordered."
+ << rParameters
+ << std::endl;
+ }
+}
+
+void ApplyConstantPhreaticMultiLinePressureProcess::InitializeCoordinates(const Parameters &rParameters)
+{
+ mXCoordinates = rParameters["x_coordinates"].GetVector();
+ mYCoordinates = rParameters["y_coordinates"].GetVector();
+ mZCoordinates = rParameters["z_coordinates"].GetVector();
+}
+
+void ApplyConstantPhreaticMultiLinePressureProcess::InitializeGravityDirection(const Parameters &rParameters)
+{
+ mGravityDirection = rParameters["gravity_direction"].GetInt();
+ switch (GravityDirection()) {
+ case 0:
+ mGravityDirectionCoordinates = XCoordinates();
+ break;
+ case 1:
+ mGravityDirectionCoordinates = YCoordinates();
+ break;
+ case 2:
+ mGravityDirectionCoordinates = ZCoordinates();
+ break;
+ default:
+ KRATOS_ERROR << "The Gravity direction is invalid";
+ }
+}
+
+void ApplyConstantPhreaticMultiLinePressureProcess::InitializeHorizontalDirection()
+{
+ mHorizontalDirection = 0;
+ for (unsigned int i=0; i<N_DIM_3D; ++i)
+ if (i != GravityDirection() && i != OutOfPlaneDirection()) mHorizontalDirection = i;
+
+ switch (HorizontalDirection()) {
+ case 0:
+ mHorizontalDirectionCoordinates = XCoordinates();
+ break;
+ case 1:
+ mHorizontalDirectionCoordinates = YCoordinates();
+ break;
+ case 2:
+ mHorizontalDirectionCoordinates = ZCoordinates();
+ break;
+ default:
+ KRATOS_ERROR << "The Horizontal direction is invalid";
+ }
+}
+
+void ApplyConstantPhreaticMultiLinePressureProcess::ExecuteInitialize()
+{
+ KRATOS_TRY
+
+ KRATOS_ERROR_IF(mrModelPart.NumberOfNodes() <= 0) << "Number of nodes is smaller than or equal to zero";
+
+ const Variable<double> &var = KratosComponents< Variable<double> >::Get(VariableName());
+
+ block_for_each(mrModelPart.Nodes(), [&var, this](Node& rNode) {
+ const double pressure = CalculatePressure(rNode);
+ if (IsSeepage()) {
+ if (pressure < PORE_PRESSURE_SIGN_FACTOR * PressureTensionCutOff()) {
+ rNode.FastGetSolutionStepValue(var) = pressure;
+ if (IsFixed()) rNode.Fix(var);
+ } else {
+ if (IsFixedProvided()) rNode.Free(var);
+ }
+ } else {
+ if (IsFixed()) rNode.Fix(var);
+ else if (IsFixedProvided()) rNode.Free(var);
+ rNode.FastGetSolutionStepValue(var) = std::min(pressure, PORE_PRESSURE_SIGN_FACTOR * PressureTensionCutOff());
+ }
+ });
+
+ KRATOS_CATCH("")
+}
+
+std::string ApplyConstantPhreaticMultiLinePressureProcess::Info() const
+{
+ return "ApplyConstantPhreaticMultiLinePressureProcess";
+}
+
+void ApplyConstantPhreaticMultiLinePressureProcess::PrintInfo(std::ostream &rOStream) const
+{
+ rOStream << "ApplyConstantPhreaticMultiLinePressureProcess";
+}
+
+const std::string &ApplyConstantPhreaticMultiLinePressureProcess::VariableName() const
+{
+ return mVariableName;
+}
+
+bool ApplyConstantPhreaticMultiLinePressureProcess::IsFixed() const
+{
+ return mIsFixed;
+}
+
+bool ApplyConstantPhreaticMultiLinePressureProcess::IsFixedProvided() const
+{
+ return mIsFixedProvided;
+}
+
+bool ApplyConstantPhreaticMultiLinePressureProcess::IsSeepage() const
+{
+ return mIsSeepage;
+}
+
+unsigned int ApplyConstantPhreaticMultiLinePressureProcess::GravityDirection() const
+{
+ return mGravityDirection;
+}
+
+double ApplyConstantPhreaticMultiLinePressureProcess::SpecificWeight() const
+{
+ return mSpecificWeight;
+}
+
+unsigned int ApplyConstantPhreaticMultiLinePressureProcess::OutOfPlaneDirection() const
+{
+ return mOutOfPlaneDirection;
+}
+
+unsigned int ApplyConstantPhreaticMultiLinePressureProcess::HorizontalDirection() const
+{
+ return mHorizontalDirection;
+}
+
+const Vector &ApplyConstantPhreaticMultiLinePressureProcess::HorizontalDirectionCoordinates() const
+{
+ return mHorizontalDirectionCoordinates;
+}
+
+const Vector &ApplyConstantPhreaticMultiLinePressureProcess::GravityDirectionCoordinates() const
+{
+ return mGravityDirectionCoordinates;
+}
+
+const Vector &ApplyConstantPhreaticMultiLinePressureProcess::XCoordinates() const
+{
+ return mXCoordinates;
+}
+
+const Vector &ApplyConstantPhreaticMultiLinePressureProcess::YCoordinates() const
+{
+ return mYCoordinates;
+}
+
+const Vector &ApplyConstantPhreaticMultiLinePressureProcess::ZCoordinates() const
+{
+ return mZCoordinates;
+}
+
+double ApplyConstantPhreaticMultiLinePressureProcess::PressureTensionCutOff() const
+{
+ return mPressureTensionCutOff;
+}
+
+int ApplyConstantPhreaticMultiLinePressureProcess::findIndex(const Node &rNode) const
+{
+ const auto& coords = rNode.Coordinates();
+ const auto number_of_coordinates = static_cast<int>(HorizontalDirectionCoordinates().size());
+ for (int index = 0; index < number_of_coordinates; ++index)
+ {
+ if (HorizontalDirectionCoordinates()[index] >= coords[HorizontalDirection()])
+ {
+ return index == 0 ? index : index - 1;
+ }
+ }
+
+ return number_of_coordinates - 1;
+}
+
+double ApplyConstantPhreaticMultiLinePressureProcess::CalculatePressure(const Node &rNode,
+ std::vector<double> deltaH) const
+{
+ // find nodes in horizontalDirectionCoordinates
+ const int firstPointIndex = findIndex(rNode);
+ const int secondPointIndex = firstPointIndex + 1;
+
+ array_1d<double, 2> y;
+ y[0] = GravityDirectionCoordinates()[firstPointIndex];
+ y[1] = GravityDirectionCoordinates()[secondPointIndex];
+
+ if (!deltaH.empty())
+ {
+ y[0] += deltaH[firstPointIndex];
+ y[1] += deltaH[secondPointIndex];
+ }
+
+
+ const double slope = (y[1] - y[0])
+ / (HorizontalDirectionCoordinates()[secondPointIndex] - HorizontalDirectionCoordinates()[firstPointIndex]);
+
+ const double height = slope * (rNode.Coordinates()[HorizontalDirection()] - HorizontalDirectionCoordinates()[firstPointIndex]) + y[0];
+ const double distance = height - rNode.Coordinates()[GravityDirection()];
+ return - PORE_PRESSURE_SIGN_FACTOR * SpecificWeight() * distance;
+}
+
+}

applications/GeoMechanicsApplication/custom_processes/apply_constant_phreatic_multi_line_pressure_process.h

@@ -0,0 +1,101 @@
+// KRATOS___
+// // ) )
+// // ___ ___
+// // ____ //___) ) // ) )
+// // / / // // / /
+// ((____/ / ((____ ((___/ / MECHANICS
+//
+// License: geo_mechanics_application/license.txt
+//
+// Main authors: Vahid Galavi
+// Jonathan Nuttall
+
+#pragma once
+
+#include <algorithm>
+#include "includes/kratos_flags.h"
+#include "includes/kratos_parameters.h"
+#include "processes/process.h"
+
+#include "geo_mechanics_application_variables.h"
+
+namespace Kratos
+{
+
+class KRATOS_API(GEO_MECHANICS_APPLICATION) ApplyConstantPhreaticMultiLinePressureProcess : public Process
+{
+
+public:
+ KRATOS_CLASS_POINTER_DEFINITION(ApplyConstantPhreaticMultiLinePressureProcess);
+
+ ApplyConstantPhreaticMultiLinePressureProcess(ModelPart& model_part, Parameters rParameters);
+ ~ApplyConstantPhreaticMultiLinePressureProcess() override = default;
+ ApplyConstantPhreaticMultiLinePressureProcess& operator=(ApplyConstantPhreaticMultiLinePressureProcess const&) = delete;
+ ApplyConstantPhreaticMultiLinePressureProcess(ApplyConstantPhreaticMultiLinePressureProcess const&) = delete;
+
+ /// this function is designed for being called at the beginning of the computations
+ /// right after reading the model and the groups
+ void ExecuteInitialize() override;
+
+ std::string Info() const override;
+ void PrintInfo(std::ostream& rOStream) const override;
+ const std::string& VariableName() const;
+ bool IsFixed() const;
+ bool IsSeepage() const;
+ unsigned int GravityDirection() const;
+ double SpecificWeight() const;
+ unsigned int OutOfPlaneDirection() const;
+ unsigned int HorizontalDirection() const;
+ const Vector& HorizontalDirectionCoordinates() const;
+ const Vector& GravityDirectionCoordinates() const;
+ const Vector& XCoordinates() const;
+ const Vector& YCoordinates() const;
+ const Vector& ZCoordinates() const;
+ double PressureTensionCutOff() const;
+ bool IsFixedProvided() const;
+
+protected:
+ ModelPart& mrModelPart;
+ int findIndex(const Node& rNode) const;
+ double CalculatePressure(const Node& rNode, std::vector<double> deltaH = {}) const;
+
+private:
+ std::string mVariableName;
+ bool mIsFixed;
+ bool mIsFixedProvided;
+ bool mIsSeepage;
+ unsigned int mGravityDirection;
+ double mSpecificWeight;
+ unsigned int mOutOfPlaneDirection;
+ unsigned int mHorizontalDirection;
+ Vector mHorizontalDirectionCoordinates;
+ Vector mGravityDirectionCoordinates;
+ Vector mXCoordinates;
+ Vector mYCoordinates;
+ Vector mZCoordinates;
+ double mPressureTensionCutOff;
+
+ void InitializeHorizontalDirection();
+ void InitializeGravityDirection(const Parameters &rParameters);
+ void InitializeCoordinates(const Parameters &rParameters);
+ void ValidateCoordinates(const Parameters &rParameters) const;
+
+ void InitializeParameters(Parameters &rParameters) const;
+}; // Class ApplyConstantPhreaticMultiLinePressureProcess
+
+/// input stream function
+inline std::istream& operator >> (std::istream& rIStream,
+ ApplyConstantPhreaticMultiLinePressureProcess& rThis);
+
+/// output stream function
+inline std::ostream& operator << (std::ostream& rOStream,
+ const ApplyConstantPhreaticMultiLinePressureProcess& rThis)
+{
+ rThis.PrintInfo(rOStream);
+ rOStream << std::endl;
+ rThis.PrintData(rOStream);
+
+ return rOStream;
+}
+
+} // namespace Kratos.