1+ // KRATOS___
2+ // // ) )
3+ // // ___ ___
4+ // // ____ //___) ) // ) )
5+ // // / / // // / /
6+ // ((____/ / ((____ ((___/ / MECHANICS
7+ //
8+ // License: geo_mechanics_application/license.txt
9+ //
10+ // Main authors: Gennady Markelov
11+ //
12+
13+ // Project includes
14+ #include " custom_strategies/strategies/geo_mechanics_newton_raphson_erosion_process_strategy.hpp"
15+ #include " geo_mechanics_fast_suite.h"
16+ #include " test_utilities.h"
17+
18+ #include < geo_mechanics_application.h>
19+ #include < linear_solvers/skyline_lu_factorization_solver.h>
20+ #include < solving_strategies/convergencecriterias/mixed_generic_criteria.h>
21+
22+ namespace Kratos
23+ {
24+ template <class TSparseSpace , class TDenseSpace , class TLinearSolver >
25+ class MockGeoMechanicsNewtonRaphsonErosionProcessStrategy
26+ : public GeoMechanicsNewtonRaphsonErosionProcessStrategy<TSparseSpace, TDenseSpace, TLinearSolver>
27+ {
28+ public:
29+ using SolvingStrategyType = ImplicitSolvingStrategy<TSparseSpace, TDenseSpace, TLinearSolver>;
30+ using TConvergenceCriteriaType = ConvergenceCriteria<TSparseSpace, TDenseSpace>;
31+ using TBuilderAndSolverType = typename SolvingStrategyType::TBuilderAndSolverType;
32+ using TSchemeType = typename SolvingStrategyType::TSchemeType;
33+
34+ MockGeoMechanicsNewtonRaphsonErosionProcessStrategy (ModelPart& rModelPart,
35+ typename TSchemeType::Pointer pScheme,
36+ typename TConvergenceCriteriaType::Pointer pNewConvergenceCriteria,
37+ typename TBuilderAndSolverType::Pointer pNewBuilderAndSolver,
38+ Parameters& rParameters)
39+ : GeoMechanicsNewtonRaphsonErosionProcessStrategy<TSparseSpace, TDenseSpace, TLinearSolver>(
40+ rModelPart, pScheme, pNewConvergenceCriteria, pNewBuilderAndSolver, rParameters),
41+ mrModelPart (rModelPart)
42+ {
43+ }
44+
45+ private:
46+ bool Recalculate () override
47+ {
48+ // The function provides a pressure change over iterations that leads to a non-linear growth of an equilibrium pipe height.
49+ // As a result the height curve crosses twice a pipe height growth curve in the search algorithm implemented
50+ // in GeoMechanicsNewtonRaphsonErosionProcessStrategy. The search stops at the second point.
51+ for (auto & node : mrModelPart.Nodes ()) {
52+ node.FastGetSolutionStepValue (WATER_PRESSURE) =
53+ std::pow (node.FastGetSolutionStepValue (WATER_PRESSURE), 0.9725 );
54+ }
55+ return true ;
56+ }
57+
58+ void BaseClassFinalizeSolutionStep () override
59+ {
60+ // to avoid calling the BaseClassFinalizeSolutionStep in tested class.
61+ }
62+
63+ ModelPart& mrModelPart;
64+ };
65+
66+ Node::Pointer CreateNodeWithSolutionStepValues (ModelPart& rModelPart, int Id, double CoordinateX, double CoordinateY, double WaterPressure)
67+ {
68+ constexpr double coordinate_z = 0.0 ;
69+ auto p_node = rModelPart.CreateNewNode (Id, CoordinateX, CoordinateY, coordinate_z);
70+
71+ p_node->FastGetSolutionStepValue (WATER_PRESSURE) = WaterPressure;
72+ const array_1d<double , 3 > GravitationalAcceleration{0 , -9.81 , 0 };
73+ p_node->FastGetSolutionStepValue (VOLUME_ACCELERATION) = GravitationalAcceleration;
74+
75+ return p_node;
76+ }
77+
78+ Geometry<Node>::Pointer CreateQuadrilateral2D4N (ModelPart& rModelPart,
79+ const std::vector<int >& rNodeIds,
80+ double Xmin,
81+ double Xmax,
82+ double WaterPressureLeft,
83+ double WaterPressureRight)
84+ {
85+ constexpr double y_min = -0.1 ;
86+ constexpr double y_max = 0.1 ;
87+
88+ return make_shared<Quadrilateral2D4<Node>>(
89+ CreateNodeWithSolutionStepValues (rModelPart, rNodeIds[0 ], Xmin, y_min, WaterPressureLeft),
90+ CreateNodeWithSolutionStepValues (rModelPart, rNodeIds[1 ], Xmax, y_min, WaterPressureRight),
91+ CreateNodeWithSolutionStepValues (rModelPart, rNodeIds[2 ], Xmax, y_max, WaterPressureRight),
92+ CreateNodeWithSolutionStepValues (rModelPart, rNodeIds[3 ], Xmin, y_max, WaterPressureLeft));
93+ }
94+
95+ auto SetupPipingStrategy (Model& rModel)
96+ {
97+ using SparseSpaceType = UblasSpace<double , CompressedMatrix, Vector>;
98+ using LocalSpaceType = UblasSpace<double , Matrix, Vector>;
99+ using LinearSolverType = SkylineLUFactorizationSolver<SparseSpaceType, LocalSpaceType>;
100+ using ConvergenceCriteriaType = ConvergenceCriteria<SparseSpaceType, LocalSpaceType>;
101+ using MixedGenericCriteriaType = MixedGenericCriteria<SparseSpaceType, LocalSpaceType>;
102+ using ConvergenceVariableListType = MixedGenericCriteriaType::ConvergenceVariableListType;
103+
104+ auto & r_model_part = rModel.CreateModelPart (" ModelPart" 1 );
105+ r_model_part.AddNodalSolutionStepVariable (WATER_PRESSURE);
106+ r_model_part.AddNodalSolutionStepVariable (VOLUME_ACCELERATION);
107+ const auto & r_process_info = r_model_part.GetProcessInfo ();
108+
109+ auto p_element_props = r_model_part.CreateNewProperties (0 );
110+ p_element_props->SetValue (CONSTITUTIVE_LAW, LinearElastic2DInterfaceLaw ().Clone ());
111+
112+ auto p_element = make_intrusive<SteadyStatePwElement<2 , 4 >>(
113+ 0 , CreateQuadrilateral2D4N (r_model_part, std::vector<int >{13 , 14 , 15 , 16 }, 3.0 , 4.0 , 2000.0 , 2000.0 ),
114+ p_element_props, std::make_unique<PlaneStrainStressState>());
115+ p_element->Initialize (r_process_info);
116+ r_model_part.AddElement (p_element);
117+
118+ // Set the pipe element properties
119+ auto p_piping_element_props = r_model_part.CreateNewProperties (1 );
120+ p_piping_element_props->SetValue (PIPE_D_70, 2e-4 );
121+ p_piping_element_props->SetValue (PIPE_ETA, 0.25 );
122+ p_piping_element_props->SetValue (PIPE_THETA, 30 );
123+ p_piping_element_props->SetValue (DENSITY_SOLID, 2650 );
124+ p_piping_element_props->SetValue (DENSITY_WATER, 1000 );
125+ p_piping_element_props->SetValue (PIPE_ELEMENT_LENGTH, 1.0 );
126+ p_piping_element_props->SetValue (PIPE_MODIFIED_D, false );
127+ p_piping_element_props->SetValue (PIPE_MODEL_FACTOR, 1 );
128+ p_piping_element_props->SetValue (PIPE_START_ELEMENT, 1 );
129+ p_piping_element_props->SetValue (CONSTITUTIVE_LAW, LinearElastic2DInterfaceLaw ().Clone ());
130+
131+ // Create the start piping element
132+ auto p_piping_element = make_intrusive<SteadyStatePwPipingElement<2 , 4 >>(
133+ 1 , CreateQuadrilateral2D4N (r_model_part, std::vector<int >{1 , 2 , 3 , 4 }, 0.0 , 1.0 , 0.0 , 500.0 ),
134+ p_piping_element_props, std::make_unique<PlaneStrainStressState>());
135+ p_piping_element->Initialize (r_process_info);
136+ r_model_part.AddElement (p_piping_element);
137+
138+ // Create other piping elements
139+ p_piping_element = make_intrusive<SteadyStatePwPipingElement<2 , 4 >>(
140+ 2 , CreateQuadrilateral2D4N (r_model_part, std::vector<int >{5 , 6 , 7 , 8 }, 2.0 , 3.0 , 500.0 , 1000.0 ),
141+ p_piping_element_props, std::make_unique<PlaneStrainStressState>());
142+ p_piping_element->Initialize (r_process_info);
143+ r_model_part.AddElement (p_piping_element);
144+
145+ p_piping_element = make_intrusive<SteadyStatePwPipingElement<2 , 4 >>(
146+ 3 , CreateQuadrilateral2D4N (r_model_part, std::vector<int >{9 , 10 , 11 , 12 }, 1.0 , 2.0 , 1000.0 , 1500.0 ),
147+ p_piping_element_props, std::make_unique<PlaneStrainStressState>());
148+ p_piping_element->Initialize (r_process_info);
149+ r_model_part.AddElement (p_piping_element);
150+
151+ Parameters p_parameters (R"(
152+ {
153+ "max_piping_iterations": 25
154+ } )" );
155+
156+ const auto p_builder_and_solver =
157+ Kratos::make_shared<ResidualBasedBlockBuilderAndSolver<SparseSpaceType, LocalSpaceType, LinearSolverType>>(nullptr );
158+ const ConvergenceVariableListType convergence_settings{};
159+ const ConvergenceCriteriaType::Pointer p_criteria =
160+ std::make_shared<MixedGenericCriteriaType>(convergence_settings);
161+
162+ return std::make_unique<MockGeoMechanicsNewtonRaphsonErosionProcessStrategy<SparseSpaceType, LocalSpaceType, LinearSolverType>>(
163+ r_model_part, nullptr , p_criteria, p_builder_and_solver, p_parameters);
164+ }
165+ } // namespace Kratos
166+
167+ namespace Kratos ::Testing
168+ {
169+ KRATOS_TEST_CASE_IN_SUITE (CheckGetPipingElements, KratosGeoMechanicsFastSuiteWithoutKernel)
170+ {
171+ // Arrange
172+ Model model;
173+ auto p_solving_strategy = SetupPipingStrategy (model);
174+
175+ // Act
176+ const auto piping_elements = p_solving_strategy->GetPipingElements ();
177+
178+ // Assert
179+ const auto elements = model.GetModelPart (" ModelPart" Elements ();
180+ constexpr int number_of_elements = 4 ;
181+ KRATOS_EXPECT_EQ (elements.size (), number_of_elements);
182+
183+ constexpr int number_of_piping_elements = 3 ;
184+ KRATOS_EXPECT_EQ (piping_elements.size (), number_of_piping_elements);
185+ KRATOS_EXPECT_EQ (piping_elements[0 ]->GetId (), 1 );
186+ KRATOS_EXPECT_EQ (piping_elements[1 ]->GetId (), 3 );
187+ KRATOS_EXPECT_EQ (piping_elements[2 ]->GetId (), 2 );
188+ }
189+
190+ KRATOS_TEST_CASE_IN_SUITE (CheckFinalizeSolutionStep, KratosGeoMechanicsFastSuiteWithoutKernel)
191+ {
192+ // Arrange
193+ Model model;
194+ auto p_solving_strategy = SetupPipingStrategy (model);
195+
196+ // Act
197+ p_solving_strategy->FinalizeSolutionStep ();
198+
199+ // Assert
200+ auto elements = model.GetModelPart (" ModelPart" Elements ();
201+
202+ constexpr int active = 1 ;
203+ KRATOS_EXPECT_EQ (elements[1 ].GetValue (PIPE_ACTIVE), active);
204+
205+ constexpr double expected_active_pipe_height = 0.0183333 ;
206+ KRATOS_EXPECT_NEAR (elements[1 ].GetValue (PIPE_HEIGHT), expected_active_pipe_height, Defaults::relative_tolerance);
207+
208+ constexpr int inactive = 0 ;
209+ KRATOS_EXPECT_EQ (elements[2 ].GetValue (PIPE_ACTIVE), inactive);
210+ KRATOS_EXPECT_EQ (elements[3 ].GetValue (PIPE_ACTIVE), inactive);
211+
212+ constexpr double expected_inactive_pipe_height = 0.0 ;
213+ KRATOS_EXPECT_NEAR (elements[2 ].GetValue (PIPE_HEIGHT), expected_inactive_pipe_height, Defaults::relative_tolerance);
214+ KRATOS_EXPECT_NEAR (elements[3 ].GetValue (PIPE_HEIGHT), expected_inactive_pipe_height, Defaults::relative_tolerance);
215+ }
216+ } // namespace Kratos::Testing
0 commit comments