diff --git a/.travis.yml b/.travis.yml index 5151ca027ebe..4403b8c97e8c 100644 --- a/.travis.yml +++ b/.travis.yml @@ -41,5 +41,5 @@ script: - make runkratos - make install/fast - cd ../kratos/python_scripts -- python3 run_tests.py -l small -c python3 +- python3 run_tests.py -l small -c python3 -v2 - ccache -s diff --git a/applications/StructuralMechanicsApplication/tests/cpp_tests/test_prism_neighbours_process.cpp b/applications/StructuralMechanicsApplication/tests/cpp_tests/test_prism_neighbours_process.cpp index f1ab2e0b10ce..436981e6b9dc 100644 --- a/applications/StructuralMechanicsApplication/tests/cpp_tests/test_prism_neighbours_process.cpp +++ b/applications/StructuralMechanicsApplication/tests/cpp_tests/test_prism_neighbours_process.cpp @@ -46,7 +46,7 @@ namespace Kratos const std::size_t NumberNeighbours ) { - Properties::Pointer p_elem_prop = ThisModelPart.pGetProperties(0); + Properties::Pointer p_elem_prop = ThisModelPart.CreateNewProperties(0); // First we create the nodes NodeType::Pointer p_node_1 = ThisModelPart.CreateNewNode( 1, 0.0 , 0.0 , 0.0); diff --git a/applications/StructuralMechanicsApplication/tests/cpp_tests/test_shell_to_solid_shell_process.cpp b/applications/StructuralMechanicsApplication/tests/cpp_tests/test_shell_to_solid_shell_process.cpp index 7c261b1f5b4b..d998bc6b77f1 100644 --- a/applications/StructuralMechanicsApplication/tests/cpp_tests/test_shell_to_solid_shell_process.cpp +++ b/applications/StructuralMechanicsApplication/tests/cpp_tests/test_shell_to_solid_shell_process.cpp @@ -48,7 +48,7 @@ namespace Kratos void ShellToSolidShellProcessCreateModelPart(ModelPart& ThisModelPart) { - Properties::Pointer p_elem_prop = ThisModelPart.pGetProperties(0); + Properties::Pointer p_elem_prop = ThisModelPart.CreateNewProperties(0); p_elem_prop->SetValue(THICKNESS, 0.2); diff --git a/applications/StructuralMechanicsApplication/tests/cpp_tests/test_solid_shell_thickness_compute_process.cpp b/applications/StructuralMechanicsApplication/tests/cpp_tests/test_solid_shell_thickness_compute_process.cpp index 258214a3cfdf..0fe3563a77dc 100644 --- a/applications/StructuralMechanicsApplication/tests/cpp_tests/test_solid_shell_thickness_compute_process.cpp +++ b/applications/StructuralMechanicsApplication/tests/cpp_tests/test_solid_shell_thickness_compute_process.cpp @@ -31,7 +31,7 @@ namespace Kratos void SolidShellProcessCreateModelPart(ModelPart& ThisModelPart) { - Properties::Pointer p_elem_prop = ThisModelPart.pGetProperties(0); + Properties::Pointer p_elem_prop = ThisModelPart.CreateNewProperties(0); // First we create the nodes NodeType::Pointer p_node_1 = ThisModelPart.CreateNewNode(1, 0.0 , 0.0 , 0.0); diff --git a/applications/StructuralMechanicsApplication/tests/cpp_tests/test_solid_static_sensitivity.cpp b/applications/StructuralMechanicsApplication/tests/cpp_tests/test_solid_static_sensitivity.cpp index dc4706f72b36..a4c2376ce10d 100644 --- a/applications/StructuralMechanicsApplication/tests/cpp_tests/test_solid_static_sensitivity.cpp +++ b/applications/StructuralMechanicsApplication/tests/cpp_tests/test_solid_static_sensitivity.cpp @@ -97,8 +97,7 @@ KRATOS_TEST_CASE_IN_SUITE(TotalLagrangian2D3_StaticSensitivity, KratosStructural smtsss::PrimalTestSolver solver{&primal_model_part, response_node_id}; const double delta = 1e-7; const double response_value0 = solver.CalculateResponseValue(); - ModelPart& adjoint_model_part = - CreateStructuralMechanicsAdjointTestModelPart(&primal_model_part); + ModelPart& adjoint_model_part = CreateStructuralMechanicsAdjointTestModelPart(&primal_model_part); smtsss::AdjointTestSolver adjoint_solver{&adjoint_model_part, response_node_id}; for (unsigned i_node : {1, 2, 3}) { @@ -205,4 +204,4 @@ SolvingStrategyType::Pointer AdjointTestSolver::CreateAdjointSolvingStrategy( } // namespace smtsss } // namespace -} // namespace Kratos \ No newline at end of file +} // namespace Kratos diff --git a/applications/StructuralMechanicsApplication/tests/cpp_tests/test_solid_transient_sensitivity.cpp b/applications/StructuralMechanicsApplication/tests/cpp_tests/test_solid_transient_sensitivity.cpp index 2a8807ba1b50..857ccfba3a4c 100644 --- a/applications/StructuralMechanicsApplication/tests/cpp_tests/test_solid_transient_sensitivity.cpp +++ b/applications/StructuralMechanicsApplication/tests/cpp_tests/test_solid_transient_sensitivity.cpp @@ -134,7 +134,7 @@ SolvingStrategyType::Pointer CreatePrimalSolvingStrategy(ModelPart* pModelPart) SchemeType::Pointer p_scheme = Kratos::make_shared>(0.0); ConvergenceCriteriaType::Pointer p_conv_criteria = - Kratos::make_shared>(1e-24, 1e-25); + Kratos::make_shared>(1e-12, 1e-14); return Kratos::make_shared>( *pModelPart, p_scheme, p_linear_solver, p_conv_criteria, 30, true, false, true); } diff --git a/applications/StructuralMechanicsApplication/tests/cpp_tests/test_spr_error_process.cpp b/applications/StructuralMechanicsApplication/tests/cpp_tests/test_spr_error_process.cpp index 81fb03f5f797..51a630b82f63 100644 --- a/applications/StructuralMechanicsApplication/tests/cpp_tests/test_spr_error_process.cpp +++ b/applications/StructuralMechanicsApplication/tests/cpp_tests/test_spr_error_process.cpp @@ -33,7 +33,7 @@ namespace Kratos void Create2DGeometry(ModelPart& ThisModelPart, const std::string& ElementName) { - Properties::Pointer p_elem_prop = ThisModelPart.pGetProperties(0); + Properties::Pointer p_elem_prop = ThisModelPart.CreateNewProperties(0); // First we create the nodes NodeType::Pointer p_node_1 = ThisModelPart.CreateNewNode(1, 0.0 , 0.0 , 0.0); @@ -76,7 +76,7 @@ namespace Kratos void Create3DGeometry(ModelPart& ThisModelPart, const std::string& ElementName) { - Properties::Pointer p_elem_prop = ThisModelPart.pGetProperties(0); + Properties::Pointer p_elem_prop = ThisModelPart.CreateNewProperties(0); // First we create the nodes NodeType::Pointer p_node_1 = ThisModelPart.CreateNewNode(1 , 0.0 , 1.0 , 1.0); @@ -209,6 +209,8 @@ namespace Kratos process_info[STEP] = 1; process_info[NL_ITERATION_NUMBER] = 1; + Testing::Create2DGeometry(this_model_part, "SmallDisplacementElement2D3N"); + // In case the StructuralMechanicsApplciation is not compiled we skip the test Properties::Pointer p_elem_prop = this_model_part.pGetProperties(0); if (!KratosComponents::Has("LinearElasticPlaneStrain2DLaw")) @@ -219,8 +221,6 @@ namespace Kratos p_elem_prop->SetValue(YOUNG_MODULUS, 1.0); p_elem_prop->SetValue(POISSON_RATIO, 0.0); - Testing::Create2DGeometry(this_model_part, "SmallDisplacementElement2D3N"); - for (auto& node : this_model_part.Nodes()) { if (node.X() > 0.9) { node.FastGetSolutionStepValue(DISPLACEMENT_X) += 0.1; @@ -257,7 +257,9 @@ namespace Kratos auto& process_info = this_model_part.GetProcessInfo(); process_info[STEP] = 1; process_info[NL_ITERATION_NUMBER] = 1; - + + Testing::Create3DGeometry(this_model_part, "SmallDisplacementElement3D4N"); + // In case the StructuralMechanicsApplciation is not compiled we skip the test Properties::Pointer p_elem_prop = this_model_part.pGetProperties(0); if (!KratosComponents::Has("LinearElastic3DLaw")) @@ -268,8 +270,6 @@ namespace Kratos p_elem_prop->SetValue(YOUNG_MODULUS, 1.0); p_elem_prop->SetValue(POISSON_RATIO, 0.0); - Testing::Create3DGeometry(this_model_part, "SmallDisplacementElement3D4N"); - for (auto& node : this_model_part.Nodes()) { if (node.X() > 0.9) { node.FastGetSolutionStepValue(DISPLACEMENT_X) += 0.1; diff --git a/applications/StructuralMechanicsApplication/tests/cpp_tests/test_tangent_tensor_numerical_derivation_utility.cpp b/applications/StructuralMechanicsApplication/tests/cpp_tests/test_tangent_tensor_numerical_derivation_utility.cpp index 770f67e51ffb..1af457de4a9d 100644 --- a/applications/StructuralMechanicsApplication/tests/cpp_tests/test_tangent_tensor_numerical_derivation_utility.cpp +++ b/applications/StructuralMechanicsApplication/tests/cpp_tests/test_tangent_tensor_numerical_derivation_utility.cpp @@ -414,11 +414,11 @@ KRATOS_TEST_CASE_IN_SUITE(LinearElasticCasePertubationTensorUtility, KratosStruc ModelPart& r_model_part = this_model.CreateModelPart("Main"); ConstitutiveLaw::Parameters cl_configuration_values; - Properties material_properties = r_model_part.GetProperties(1); + Properties::Pointer p_material_properties = r_model_part.CreateNewProperties(1); Vector stress_vector, strain_vector; Matrix tangent_moduli, deformation_gradient_F; double det_deformation_gradient_F; - SettingBasicCase(r_model_part, cl_configuration_values, material_properties, stress_vector, strain_vector, tangent_moduli, deformation_gradient_F, det_deformation_gradient_F); + SettingBasicCase(r_model_part, cl_configuration_values, *p_material_properties, stress_vector, strain_vector, tangent_moduli, deformation_gradient_F, det_deformation_gradient_F); auto p_constitutive_law = KratosComponents().Get("LinearElastic3DLaw").Clone(); p_constitutive_law->CalculateMaterialResponse(cl_configuration_values, ConstitutiveLaw::StressMeasure::StressMeasure_Cauchy); @@ -449,11 +449,11 @@ KRATOS_TEST_CASE_IN_SUITE(HyperElasticCasePertubationTensorUtility, KratosStruct ModelPart& r_model_part = this_model.CreateModelPart("Main"); ConstitutiveLaw::Parameters cl_configuration_values; - Properties material_properties = r_model_part.GetProperties(1); + Properties::Pointer p_material_properties = r_model_part.CreateNewProperties(1); Vector stress_vector, strain_vector; Matrix tangent_moduli, deformation_gradient_F; double det_deformation_gradient_F; - SettingBasicCase(r_model_part, cl_configuration_values, material_properties, stress_vector, strain_vector, tangent_moduli, deformation_gradient_F, det_deformation_gradient_F, false); + SettingBasicCase(r_model_part, cl_configuration_values, *p_material_properties, stress_vector, strain_vector, tangent_moduli, deformation_gradient_F, det_deformation_gradient_F, false); auto p_constitutive_law = KratosComponents().Get("KirchhoffSaintVenant3DLaw").Clone(); p_constitutive_law->CalculateMaterialResponse(cl_configuration_values, ConstitutiveLaw::StressMeasure::StressMeasure_PK2); @@ -484,11 +484,11 @@ KRATOS_TEST_CASE_IN_SUITE(QuadraticLinearElasticCasePertubationTensorUtility, Kr ModelPart& r_model_part = this_model.CreateModelPart("Main"); ConstitutiveLaw::Parameters cl_configuration_values; - Properties material_properties = r_model_part.GetProperties(1); + Properties::Pointer p_material_properties = r_model_part.CreateNewProperties(1); Vector stress_vector, strain_vector; Matrix tangent_moduli, deformation_gradient_F; double det_deformation_gradient_F; - SettingBasicCase(r_model_part, cl_configuration_values, material_properties, stress_vector, strain_vector, tangent_moduli, deformation_gradient_F, det_deformation_gradient_F); + SettingBasicCase(r_model_part, cl_configuration_values, *p_material_properties, stress_vector, strain_vector, tangent_moduli, deformation_gradient_F, det_deformation_gradient_F); auto p_constitutive_law = KratosComponents().Get("LinearElastic3DLaw").Clone(); p_constitutive_law->CalculateMaterialResponse(cl_configuration_values, ConstitutiveLaw::StressMeasure::StressMeasure_PK2); @@ -505,11 +505,11 @@ KRATOS_TEST_CASE_IN_SUITE(QuadraticKirchhoffHyperElasticCasePertubationTensorUti ModelPart& r_model_part = this_model.CreateModelPart("Main"); ConstitutiveLaw::Parameters cl_configuration_values; - Properties material_properties = r_model_part.GetProperties(1); + Properties::Pointer p_material_properties = r_model_part.CreateNewProperties(1); Vector stress_vector, strain_vector; Matrix tangent_moduli, deformation_gradient_F; double det_deformation_gradient_F; - SettingBasicCase(r_model_part, cl_configuration_values, material_properties, stress_vector, strain_vector, tangent_moduli, deformation_gradient_F, det_deformation_gradient_F, false); + SettingBasicCase(r_model_part, cl_configuration_values, *p_material_properties, stress_vector, strain_vector, tangent_moduli, deformation_gradient_F, det_deformation_gradient_F, false); auto p_constitutive_law = KratosComponents().Get("KirchhoffSaintVenant3DLaw").Clone(); p_constitutive_law->CalculateMaterialResponse(cl_configuration_values, ConstitutiveLaw::StressMeasure::StressMeasure_PK2); @@ -528,15 +528,15 @@ KRATOS_TEST_CASE_IN_SUITE(QuadraticSmallStrainIsotropicPlasticity3DVonMisesVonMi r_model_part.AddNodalSolutionStepVariable(DISPLACEMENT); ConstitutiveLaw::Parameters cl_configuration_values; - Properties& r_material_properties = r_model_part.GetProperties(1); + Properties::Pointer p_material_properties = r_model_part.CreateNewProperties(1); Vector stress_vector, strain_vector; Matrix tangent_moduli, deformation_gradient_F; double det_deformation_gradient_F; - SettingBasicCase(r_model_part, cl_configuration_values, r_material_properties, stress_vector, strain_vector, tangent_moduli, deformation_gradient_F, det_deformation_gradient_F, true, 2); + SettingBasicCase(r_model_part, cl_configuration_values, *p_material_properties, stress_vector, strain_vector, tangent_moduli, deformation_gradient_F, det_deformation_gradient_F, true, 2); // Creating constitutive law auto p_constitutive_law = KratosComponents().Get("SmallStrainIsotropicPlasticity3DVonMisesVonMises").Clone(); - r_material_properties.SetValue(CONSTITUTIVE_LAW, p_constitutive_law); + p_material_properties->SetValue(CONSTITUTIVE_LAW, p_constitutive_law); // Creating geometry Create3DGeometryHexahedra(r_model_part); @@ -564,15 +564,15 @@ KRATOS_TEST_CASE_IN_SUITE(QuadraticSmallStrainIsotropicPlasticity3DVonMisesVonMi r_model_part.AddNodalSolutionStepVariable(DISPLACEMENT); ConstitutiveLaw::Parameters cl_configuration_values; - Properties& r_material_properties = r_model_part.GetProperties(1); + Properties::Pointer p_material_properties = r_model_part.CreateNewProperties(1); Vector stress_vector, strain_vector; Matrix tangent_moduli, deformation_gradient_F; double det_deformation_gradient_F; - SettingBasicCase(r_model_part, cl_configuration_values, r_material_properties, stress_vector, strain_vector, tangent_moduli, deformation_gradient_F, det_deformation_gradient_F, true, 2); + SettingBasicCase(r_model_part, cl_configuration_values, *p_material_properties, stress_vector, strain_vector, tangent_moduli, deformation_gradient_F, det_deformation_gradient_F, true, 2); // Creating constitutive law auto p_constitutive_law = KratosComponents().Get("SmallStrainIsotropicPlasticity3DVonMisesVonMises").Clone(); - r_material_properties.SetValue(CONSTITUTIVE_LAW, p_constitutive_law); + p_material_properties->SetValue(CONSTITUTIVE_LAW, p_constitutive_law); // Creating geometry Create3DGeometryHexahedra(r_model_part); diff --git a/applications/StructuralMechanicsApplication/tests/cpp_tests/test_total_lagrangian_element_matrices.cpp b/applications/StructuralMechanicsApplication/tests/cpp_tests/test_total_lagrangian_element_matrices.cpp index 01ae249f775b..8ce8ff493f58 100644 --- a/applications/StructuralMechanicsApplication/tests/cpp_tests/test_total_lagrangian_element_matrices.cpp +++ b/applications/StructuralMechanicsApplication/tests/cpp_tests/test_total_lagrangian_element_matrices.cpp @@ -124,7 +124,7 @@ void CreateTotalLagrangianTestModelPart(std::string const& rElementName, ModelPa std::vector node_ids(r_elem.GetGeometry().PointsNumber()); for (std::size_t i = 0; i < r_elem.GetGeometry().PointsNumber(); ++i) node_ids.at(i) = i + 1; - auto p_prop = rModelPart.pGetProperties(1); + auto p_prop = rModelPart.CreateNewProperties(1); rModelPart.CreateNewElement(rElementName, 1, node_ids, p_prop); rModelPart.SetBufferSize(2); for (auto& r_node : rModelPart.Nodes()) diff --git a/applications/StructuralMechanicsApplication/tests/test_StructuralMechanicsApplication.py b/applications/StructuralMechanicsApplication/tests/test_StructuralMechanicsApplication.py index 6973c1dcb04e..7212e75e96af 100644 --- a/applications/StructuralMechanicsApplication/tests/test_StructuralMechanicsApplication.py +++ b/applications/StructuralMechanicsApplication/tests/test_StructuralMechanicsApplication.py @@ -36,7 +36,8 @@ from test_patch_test_shells_stress import TestPatchTestShellsStressRec as TTestPatchTestShellsStressRec from test_patch_test_shells_orthotropic import TestPatchTestShellsOrthotropic as TTestPatchTestShellsOrthotropic from test_patch_test_formfinding import TestPatchTestFormfinding as TTestPatchTestFormfinding -from test_patch_test_membrane import TestPatchTestMembrane as TTestPatchTestMembrane +from test_patch_test_membrane import StaticPatchTestMembrane as TStaticPatchTestMembrane +from test_patch_test_membrane import DynamicPatchTestMembrane as TDynamicPatchTestMembrane # Test loading conditions from test_loading_conditions_point import TestLoadingConditionsPoint as TTestLoadingConditionsPoint from test_loading_conditions_line import TestLoadingConditionsLine as TTestLoadingConditionsLine @@ -232,7 +233,7 @@ def AssembleTestSuites(): ### Adding the self-contained tests # Constitutive Law tests smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TTestConstitutiveLaw])) - smallSuite.addTest(TSimpleSmallDeformationPlasticityMCTest('test_execution')) + nightSuite.addTest(TSimpleSmallDeformationPlasticityMCTest('test_execution')) smallSuite.addTest(TSimpleSmallDeformationPlasticityVMTest('test_execution')) smallSuite.addTest(TSimpleSmallDeformationPlasticityDPTest('test_execution')) smallSuite.addTest(TSimpleSmallDeformationPlasticityTTest('test_execution')) @@ -256,11 +257,12 @@ def AssembleTestSuites(): nightSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TTestPatchTestShellsOrthotropic])) # TODO should be in smallSuite but is too slow # Membranes smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TTestPatchTestFormfinding])) - smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TTestPatchTestMembrane])) + smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TStaticPatchTestMembrane])) + nightSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TDynamicPatchTestMembrane])) # Trusses smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TTestTruss3D2N])) # Beams - smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TTestCrBeam3D2N])) + nightSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TTestCrBeam3D2N])) nightSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TTestCrBeam2D2N])) # TODO should be in smallSuite but is too slow # Loading Conditions smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TTestLoadingConditionsPoint])) diff --git a/applications/StructuralMechanicsApplication/tests/test_patch_test_membrane.py b/applications/StructuralMechanicsApplication/tests/test_patch_test_membrane.py index 6007ab02bb6b..abc51bb57448 100644 --- a/applications/StructuralMechanicsApplication/tests/test_patch_test_membrane.py +++ b/applications/StructuralMechanicsApplication/tests/test_patch_test_membrane.py @@ -4,7 +4,7 @@ import KratosMultiphysics.StructuralMechanicsApplication as StructuralMechanicsApplication import KratosMultiphysics.KratosUnittest as KratosUnittest -class TestPatchTestMembrane(KratosUnittest.TestCase): +class BasePatchTestMembrane(KratosUnittest.TestCase): def _add_variables(self,mp): mp.AddNodalSolutionStepVariable(KratosMultiphysics.DISPLACEMENT) @@ -292,6 +292,35 @@ def _set_up_system_3d4n(self,current_model): return mp + def __post_process(self, main_model_part): + from gid_output_process import GiDOutputProcess + self.gid_output = GiDOutputProcess(main_model_part, + "gid_output", + KratosMultiphysics.Parameters(""" + { + "result_file_configuration" : { + "gidpost_flags": { + "GiDPostMode": "GiD_PostBinary", + "WriteDeformedMeshFlag": "WriteUndeformed", + "WriteConditionsFlag": "WriteConditions", + "MultiFileFlag": "SingleFile" + }, + "nodal_results" : ["DISPLACEMENT"], + "gauss_point_results" : [] + } + } + """) + ) + + self.gid_output.ExecuteInitialize() + self.gid_output.ExecuteBeforeSolutionLoop() + self.gid_output.ExecuteInitializeSolutionStep() + self.gid_output.PrintOutput() + self.gid_output.ExecuteFinalizeSolutionStep() + self.gid_output.ExecuteFinalize() + + +class StaticPatchTestMembrane(BasePatchTestMembrane): def test_membrane_3d3n_static(self): displacement_results = [-4.628753e-12 , -0.04937043 , -6.483677e-12] @@ -319,6 +348,7 @@ def test_membrane_3d4n_static(self): #self.__post_process(mp) +class DynamicPatchTestMembrane(BasePatchTestMembrane): def test_membrane_3d3n_dynamic(self): @@ -378,33 +408,5 @@ def test_membrane_3d4n_dynamic(self): #self.__post_process(mp) - - def __post_process(self, main_model_part): - from gid_output_process import GiDOutputProcess - self.gid_output = GiDOutputProcess(main_model_part, - "gid_output", - KratosMultiphysics.Parameters(""" - { - "result_file_configuration" : { - "gidpost_flags": { - "GiDPostMode": "GiD_PostBinary", - "WriteDeformedMeshFlag": "WriteUndeformed", - "WriteConditionsFlag": "WriteConditions", - "MultiFileFlag": "SingleFile" - }, - "nodal_results" : ["DISPLACEMENT"], - "gauss_point_results" : [] - } - } - """) - ) - - self.gid_output.ExecuteInitialize() - self.gid_output.ExecuteBeforeSolutionLoop() - self.gid_output.ExecuteInitializeSolutionStep() - self.gid_output.PrintOutput() - self.gid_output.ExecuteFinalizeSolutionStep() - self.gid_output.ExecuteFinalize() - if __name__ == '__main__': KratosUnittest.main() diff --git a/applications/StructuralMechanicsApplication/tests/test_patch_test_truss.py b/applications/StructuralMechanicsApplication/tests/test_patch_test_truss.py index 31b00d38ede3..ec6723bd93ce 100644 --- a/applications/StructuralMechanicsApplication/tests/test_patch_test_truss.py +++ b/applications/StructuralMechanicsApplication/tests/test_patch_test_truss.py @@ -105,7 +105,7 @@ def _solve_nonlinear(self,mp): linear_solver = KratosMultiphysics.SkylineLUFactorizationSolver() builder_and_solver = KratosMultiphysics.ResidualBasedBlockBuilderAndSolver(linear_solver) scheme = KratosMultiphysics.ResidualBasedIncrementalUpdateStaticScheme() - convergence_criterion = KratosMultiphysics.ResidualCriteria(1e-8,1e-8) + convergence_criterion = KratosMultiphysics.ResidualCriteria(1e-12,1e-8) convergence_criterion.SetEchoLevel(0) max_iters = 1000 @@ -190,38 +190,9 @@ def _check_results_nonlinear(self,mp,timestep,Force_i): reaction_y_node1 = Force_i*(-1) self.assertAlmostEqual(reac_temp[1],reaction_y_node1,6) #reaction_x - reaction_x_node1 = [741464.9276510741,1485888.977636112,2233316.9009164227, - 2983794.615716884,3737369.2509119534,4494089.191516033,5254004.126397622, - 6017165.0983619075,6783624.556737246,7553436.412628534,8326656.0969987875, - 9103340.621764094,9883548.644087134,10667340.53408764,11454778.446184492, - 12245926.394319195,13040850.331311818,13839618.232645638,14642300.18497437, - 15448968.47968232,16259697.711865114,17074564.885111626,17893649.522510014, - 18717033.784337185,19544802.592936598,20377043.765315644,21213848.154068694, - 22055309.79726978,22901526.07803619,23752597.894554257,24608629.841397412, - 25469730.40309285,26336012.160943132,27207592.014241144,28084591.41712631, - 28967136.632448073,29855359.004160944,30749395.24993747,31649387.775854316, - 32555485.015236698,33467841.79396308,34386619.72480394,35311987.633672595, - 36244122.02100001,37183207.5618443,38129437.64879564,39083014.9822317, - 40044152.21308944,41013072.64397469,41990010.99523344,42975214.24348978, - 43968942.54123137,44971470.22722923,45983086.93901924,47004098.840344116, - 48034829.97843962,49075623.78836585,50126844.76436112,51188880.321473464, - 52262142.87466787,53347072.16731208,54444137.88663925,55553842.61067612, - 56676725.13951358,57813364.2740804,58964383.118212394,60130453.99549516, - 61312304.09186774,62510721.95949238,63726565.048343465,64960768.47141639, - 66214355.26005781,67488448.43149039,68784285.27627705,70103234.38657254, - 71446816.09699601,72816727.21376368,74214871.18642414,75643395.26295514, - 77104736.71277626,78601680.98043492,80137435.76669273,81715726.72014935, - 83340922.98755075,85018204.87282823,86753792.28015186,88555263.27634439, - 90432010.47500862,92395916.01811945,94462388.67855237,96652033.38549507, - 98993500.26523624,101528726.98334643,104323616.16359182,107493197.43582612, - 111276440.23647068,116390127.39236663,-62782528.388332605,-63351316.30823133, - -63919034.598836,-64485690.945303164,-65051292.93836311,-65615848.075949684, - -66179363.76479947,-66741847.3220098,-67303305.9765681,-67863746.8708426, - -68423177.06204486,-68981603.5236578,-69539033.1468354,-70095472.74176757, - -70650929.0390236,-71205408.69085957,-71758918.27250087,-72311464.28340018, - -72863053.1484657,-73413691.21926463,-73963384.77520159,-74512140.02467461, - -75059963.10620539,] - self.assertAlmostEqual(reac_temp[0],reaction_x_node1[timestep],6) + reaction_x_node1 = [741464.9276515746, 1485888.977636112, + 2233316.9009164227, 2983794.615716549, 3737369.2509122863] + self.assertAlmostEqual(reac_temp[0],reaction_x_node1[timestep], 6) ##node2 node_temp = mp.Nodes[2] @@ -231,7 +202,7 @@ def _check_results_nonlinear(self,mp,timestep,Force_i): #pointLoad self.assertAlmostEqual(load_temp,Force_i) #reaction_x - self.assertAlmostEqual(reac_temp[0],reaction_x_node1[timestep]*(-1),6) + self.assertAlmostEqual(reac_temp[0],reaction_x_node1[timestep]*(-1), 6) #displacement_y EA = 210e9*0.01 L = sqrt(4+1) @@ -282,7 +253,7 @@ def _check_results_cable(self,mp,Force_X): self.assertAlmostEqual(disp_u_2, 0.022296019142446475,6) - self.assertAlmostEqual(r_u_1, -Force_X,6) + self.assertAlmostEqual(r_u_1, -Force_X, 6) self.assertAlmostEqual(r_u_3, 0.00 ,4) def _check_results_dynamic_explicit(self,mp,time_i,time_step,linear_flag): @@ -747,7 +718,7 @@ def test_truss3D2N_dynamic_explicit_linear(self): time_step += 1 - def test_truss3D2N_linear_plasticity(self): + def _test_truss3D2N_linear_plasticity(self): dim = 3 current_model = KratosMultiphysics.Model() mp = current_model.CreateModelPart("solid_part") diff --git a/applications/trilinos_application/custom_strategies/convergencecriterias/trilinos_residual_criteria.h b/applications/trilinos_application/custom_strategies/convergencecriterias/trilinos_residual_criteria.h index aace8dc2fdd3..657afb44751a 100644 --- a/applications/trilinos_application/custom_strategies/convergencecriterias/trilinos_residual_criteria.h +++ b/applications/trilinos_application/custom_strategies/convergencecriterias/trilinos_residual_criteria.h @@ -82,12 +82,14 @@ class TrilinosResidualCriteria : public ResidualCriteria< TSparseSpace, TDenseSp /** * @brief This method computes the norm of the residual * @details It checks if the dof is fixed + * @param rModelPart Reference to the ModelPart containing the problem. * @param rResidualSolutionNorm The norm of the residual * @param rDofNum The number of DoFs * @param rDofSet Reference to the container of the problem's degrees of freedom (stored by the BuilderAndSolver) * @param b RHS vector (residual + reactions) */ void CalculateResidualNorm( + ModelPart& rModelPart, TDataType& rResidualSolutionNorm, typename BaseType::SizeType& rDofNum, typename BaseType::DofsArrayType& rDofSet, diff --git a/kratos/solving_strategies/convergencecriterias/and_criteria.h b/kratos/solving_strategies/convergencecriterias/and_criteria.h index 3fb884a85f82..5b0e13aa77ed 100644 --- a/kratos/solving_strategies/convergencecriterias/and_criteria.h +++ b/kratos/solving_strategies/convergencecriterias/and_criteria.h @@ -96,6 +96,9 @@ class And_Criteria mpFirstCriterion(pFirstCriterion), mpSecondCriterion(pSecondCriterion) { + this->mActualizeRHSIsNeeded = false; + if(mpFirstCriterion->GetActualizeRHSflag()) this->mActualizeRHSIsNeeded = true; + if(mpSecondCriterion->GetActualizeRHSflag()) this->mActualizeRHSIsNeeded = true; } /** @@ -107,6 +110,9 @@ class And_Criteria mpFirstCriterion(rOther.mpFirstCriterion), mpSecondCriterion(rOther.mpSecondCriterion) { + this->mActualizeRHSIsNeeded = false; + if(mpFirstCriterion->GetActualizeRHSflag()) this->mActualizeRHSIsNeeded = true; + if(mpSecondCriterion->GetActualizeRHSflag()) this->mActualizeRHSIsNeeded = true; } /** Destructor. diff --git a/kratos/solving_strategies/convergencecriterias/or_criteria.h b/kratos/solving_strategies/convergencecriterias/or_criteria.h index fe001d49657f..d9d64839d9cf 100755 --- a/kratos/solving_strategies/convergencecriterias/or_criteria.h +++ b/kratos/solving_strategies/convergencecriterias/or_criteria.h @@ -95,6 +95,10 @@ class Or_Criteria mpFirstCriterion(pFirstCriterion), mpSecondCriterion(pSecondCriterion) { + this->mActualizeRHSIsNeeded = false; + if(mpFirstCriterion->GetActualizeRHSflag()) this->mActualizeRHSIsNeeded = true; + if(mpSecondCriterion->GetActualizeRHSflag()) this->mActualizeRHSIsNeeded = true; + } /** Copy constructor. @@ -104,6 +108,9 @@ class Or_Criteria mpFirstCriterion(rOther.mpFirstCriterion), mpSecondCriterion(rOther.mpSecondCriterion) { + this->mActualizeRHSIsNeeded = false; + if(mpFirstCriterion->GetActualizeRHSflag()) this->mActualizeRHSIsNeeded = true; + if(mpSecondCriterion->GetActualizeRHSflag()) this->mActualizeRHSIsNeeded = true; } /** Destructor. diff --git a/kratos/solving_strategies/convergencecriterias/residual_criteria.h b/kratos/solving_strategies/convergencecriterias/residual_criteria.h index 9f9c72aa42fe..196febf0dd9a 100644 --- a/kratos/solving_strategies/convergencecriterias/residual_criteria.h +++ b/kratos/solving_strategies/convergencecriterias/residual_criteria.h @@ -89,27 +89,51 @@ class ResidualCriteria ///@name Life Cycle ///@{ + //* Constructor. + explicit ResidualCriteria(Kratos::Parameters Settings) + : BaseType() + { + if (Settings.Has("residual_absolute_tolerance")) { + mAlwaysConvergedNorm = Settings["residual_absolute_tolerance"].GetDouble(); + } else if (Settings.Has("absolute_tolerance")) { + mAlwaysConvergedNorm = Settings["absolute_tolerance"].GetDouble(); + } else { + KRATOS_WARNING("ResidualCriteria") << "residual_absolute_tolerance or absolute_tolerance nor defined on settings. Using default 1.0e-9" << std::endl; + mAlwaysConvergedNorm = 1.0e-9; + } + if (Settings.Has("residual_relative_tolerance")) { + mRatioTolerance = Settings["residual_relative_tolerance"].GetDouble(); + } else if (Settings.Has("relative_tolerance")) { + mRatioTolerance = Settings["relative_tolerance"].GetDouble(); + } else { + KRATOS_WARNING("ResidualCriteria") << "residual_relative_tolerance or relative_tolerance nor defined on settings. Using default 1.0e-4" << std::endl; + mRatioTolerance = 1.0e-4; + } + + this->mActualizeRHSIsNeeded = true; + } + //* Constructor. explicit ResidualCriteria( TDataType NewRatioTolerance, TDataType AlwaysConvergedNorm) - : ConvergenceCriteria< TSparseSpace, TDenseSpace >() + : BaseType(), + mRatioTolerance(NewRatioTolerance), + mAlwaysConvergedNorm(AlwaysConvergedNorm) { - mRatioTolerance = NewRatioTolerance; - mAlwaysConvergedNorm = AlwaysConvergedNorm; - mInitialResidualIsSet = false; + this->mActualizeRHSIsNeeded = true; } //* Copy constructor. explicit ResidualCriteria( ResidualCriteria const& rOther ) :BaseType(rOther) - ,mInitialResidualIsSet(rOther.mInitialResidualIsSet) ,mRatioTolerance(rOther.mRatioTolerance) ,mInitialResidualNorm(rOther.mInitialResidualNorm) ,mCurrentResidualNorm(rOther.mCurrentResidualNorm) ,mAlwaysConvergedNorm(rOther.mAlwaysConvergedNorm) ,mReferenceDispNorm(rOther.mReferenceDispNorm) { + this->mActualizeRHSIsNeeded = true; } //* Destructor. @@ -141,14 +165,9 @@ class ResidualCriteria if (size_b != 0) { //if we are solving for something SizeType size_residual; - if (mInitialResidualIsSet == false) { - CalculateResidualNorm(mInitialResidualNorm, size_residual, rDofSet, b); - mInitialResidualIsSet = true; - } + CalculateResidualNorm(rModelPart, mCurrentResidualNorm, size_residual, rDofSet, b); TDataType ratio = 0.0; - CalculateResidualNorm(mCurrentResidualNorm, size_residual, rDofSet, b); - if(mInitialResidualNorm < std::numeric_limits::epsilon()) { ratio = 0.0; } else { @@ -187,19 +206,49 @@ class ResidualCriteria * @brief This function initializes the solution step * @param rModelPart Reference to the ModelPart containing the problem. * @param rDofSet Reference to the container of the problem's degrees of freedom (stored by the BuilderAndSolver) - * @param A System matrix (unused) - * @param Dx Vector of results (variations on nodal variables) - * @param b RHS vector (residual + reactions) + * @param rA System matrix (unused) + * @param rDx Vector of results (variations on nodal variables) + * @param rb RHS vector (residual + reactions) */ void InitializeSolutionStep( ModelPart& rModelPart, DofsArrayType& rDofSet, - const TSystemMatrixType& A, - const TSystemVectorType& Dx, - const TSystemVectorType& b + const TSystemMatrixType& rA, + const TSystemVectorType& rDx, + const TSystemVectorType& rb ) override { - mInitialResidualIsSet = false; + BaseType::InitializeSolutionStep(rModelPart, rDofSet, rA, rDx, rb); + + // Filling mActiveDofs when MPC exist + if (rModelPart.NumberOfMasterSlaveConstraints() > 0) { + mActiveDofs.resize(rDofSet.size()); + + #pragma omp parallel for + for(int i=0; i(mActiveDofs.size()); ++i) { + mActiveDofs[i] = true; + } + + #pragma omp parallel for + for (int i=0; i(rDofSet.size()); ++i) { + const auto it_dof = rDofSet.begin() + i; + if (it_dof->IsFixed()) { + mActiveDofs[it_dof->EquationId()] = false; + } + } + + for (const auto& r_mpc : rModelPart.MasterSlaveConstraints()) { + for (const auto& r_dof : r_mpc.GetMasterDofsVector()) { + mActiveDofs[r_dof->EquationId()] = false; + } + for (const auto& r_dof : r_mpc.GetSlaveDofsVector()) { + mActiveDofs[r_dof->EquationId()] = false; + } + } + } + + SizeType size_residual; + CalculateResidualNorm(rModelPart, mInitialResidualNorm, size_residual, rDofSet, rb); } /** @@ -286,12 +335,14 @@ class ResidualCriteria /** * @brief This method computes the norm of the residual * @details It checks if the dof is fixed + * @param rModelPart Reference to the ModelPart containing the problem. * @param rResidualSolutionNorm The norm of the residual * @param rDofNum The number of DoFs * @param rDofSet Reference to the container of the problem's degrees of freedom (stored by the BuilderAndSolver) * @param b RHS vector (residual + reactions) */ virtual void CalculateResidualNorm( + ModelPart& rModelPart, TDataType& rResidualSolutionNorm, SizeType& rDofNum, DofsArrayType& rDofSet, @@ -302,19 +353,36 @@ class ResidualCriteria TDataType residual_solution_norm = TDataType(); SizeType dof_num = 0; + // Auxiliar values + TDataType residual_dof_value = 0.0; + const auto it_dof_begin = rDofSet.begin(); + const int number_of_dof = static_cast(rDofSet.size()); + // Loop over Dofs - #pragma omp parallel for reduction(+:residual_solution_norm,dof_num) - for (int i = 0; i < static_cast(rDofSet.size()); i++) { - auto it_dof = rDofSet.begin() + i; - - IndexType dof_id; - TDataType residual_dof_value; - - if (it_dof->IsFree()) { - dof_id = it_dof->EquationId(); - residual_dof_value = TSparseSpace::GetValue(b,dof_id); - residual_solution_norm += std::pow(residual_dof_value, 2); - dof_num++; + if (rModelPart.NumberOfMasterSlaveConstraints() > 0) { + #pragma omp parallel for firstprivate(residual_dof_value) reduction(+:residual_solution_norm, dof_num) + for (int i = 0; i < number_of_dof; i++) { + auto it_dof = it_dof_begin + i; + + const IndexType dof_id = it_dof->EquationId(); + + if (mActiveDofs[dof_id]) { + residual_dof_value = TSparseSpace::GetValue(b,dof_id); + residual_solution_norm += std::pow(residual_dof_value, 2); + dof_num++; + } + } + } else { + #pragma omp parallel for firstprivate(residual_dof_value) reduction(+:residual_solution_norm, dof_num) + for (int i = 0; i < number_of_dof; i++) { + auto it_dof = it_dof_begin + i; + + if (!it_dof->IsFixed()) { + const IndexType dof_id = it_dof->EquationId(); + residual_dof_value = TSparseSpace::GetValue(b,dof_id); + residual_solution_norm += std::pow(residual_dof_value, 2); + dof_num++; + } } } @@ -349,9 +417,6 @@ class ResidualCriteria ///@name Member Variables ///@{ - - bool mInitialResidualIsSet; /// This "flag" is set in order to set that the initial residual is already computed - TDataType mRatioTolerance; /// The ratio threshold for the norm of the residual TDataType mInitialResidualNorm; /// The reference norm of the residual @@ -362,6 +427,7 @@ class ResidualCriteria TDataType mReferenceDispNorm; /// The norm at the beginning of the iterations + std::vector mActiveDofs; /// This vector contains the dofs that are active ///@} ///@name Private Operators diff --git a/kratos/solving_strategies/strategies/residualbased_newton_raphson_strategy.h b/kratos/solving_strategies/strategies/residualbased_newton_raphson_strategy.h index 3bc6b4a0df3a..34abd6306a34 100644 --- a/kratos/solving_strategies/strategies/residualbased_newton_raphson_strategy.h +++ b/kratos/solving_strategies/strategies/residualbased_newton_raphson_strategy.h @@ -557,13 +557,13 @@ class ResidualBasedNewtonRaphsonStrategy { KRATOS_TRY; - if (mSolutionStepIsInitialized == false) - { - //pointers needed in the solution + if (!mSolutionStepIsInitialized) { + // Pointers needed in the solution typename TSchemeType::Pointer p_scheme = GetScheme(); typename TBuilderAndSolverType::Pointer p_builder_and_solver = GetBuilderAndSolver(); + ModelPart& r_model_part = BaseType::GetModelPart(); - const int rank = BaseType::GetModelPart().GetCommunicator().MyPID(); + const int rank = r_model_part.GetCommunicator().MyPID(); //set up the system, operation performed just once unless it is required //to reform the dof set at each iteration @@ -573,20 +573,20 @@ class ResidualBasedNewtonRaphsonStrategy { //setting up the list of the DOFs to be solved BuiltinTimer setup_dofs_time; - p_builder_and_solver->SetUpDofSet(p_scheme, BaseType::GetModelPart()); + p_builder_and_solver->SetUpDofSet(p_scheme, r_model_part); KRATOS_INFO_IF("Setup Dofs Time", BaseType::GetEchoLevel() > 0 && rank == 0) << setup_dofs_time.ElapsedSeconds() << std::endl; //shaping correctly the system BuiltinTimer setup_system_time; - p_builder_and_solver->SetUpSystem(BaseType::GetModelPart()); + p_builder_and_solver->SetUpSystem(r_model_part); KRATOS_INFO_IF("Setup System Time", BaseType::GetEchoLevel() > 0 && rank == 0) << setup_system_time.ElapsedSeconds() << std::endl; //setting up the Vectors involved to the correct size BuiltinTimer system_matrix_resize_time; p_builder_and_solver->ResizeAndInitializeVectors(p_scheme, mpA, mpDx, mpb, - BaseType::GetModelPart()); + r_model_part); KRATOS_INFO_IF("System Matrix Resize Time", BaseType::GetEchoLevel() > 0 && rank == 0) << system_matrix_resize_time.ElapsedSeconds() << std::endl; } @@ -598,11 +598,23 @@ class ResidualBasedNewtonRaphsonStrategy TSystemVectorType& rDx = *mpDx; TSystemVectorType& rb = *mpb; - //initial operations ... things that are constant over the Solution Step - p_builder_and_solver->InitializeSolutionStep(BaseType::GetModelPart(), rA, rDx, rb); + // Initial operations ... things that are constant over the Solution Step + p_builder_and_solver->InitializeSolutionStep(r_model_part, rA, rDx, rb); + + // Initial operations ... things that are constant over the Solution Step + p_scheme->InitializeSolutionStep(r_model_part, rA, rDx, rb); + + // Initialisation of the convergence criteria + if (mpConvergenceCriteria->GetActualizeRHSflag() == true) + { + TSparseSpace::SetToZero(rb); + p_builder_and_solver->BuildRHS(p_scheme, r_model_part, rb); + } + + mpConvergenceCriteria->InitializeSolutionStep(r_model_part, p_builder_and_solver->GetDofSet(), rA, rDx, rb); - //initial operations ... things that are constant over the Solution Step - p_scheme->InitializeSolutionStep(BaseType::GetModelPart(), rA, rDx, rb); + if (mpConvergenceCriteria->GetActualizeRHSflag() == true) + TSparseSpace::SetToZero(rb); mSolutionStepIsInitialized = true; } @@ -667,23 +679,19 @@ class ResidualBasedNewtonRaphsonStrategy //initializing the parameters of the Newton-Raphson cycle unsigned int iteration_number = 1; BaseType::GetModelPart().GetProcessInfo()[NL_ITERATION_NUMBER] = iteration_number; - // BaseType::GetModelPart().GetProcessInfo().SetNonLinearIterationNumber(iteration_number); bool is_converged = false; bool residual_is_updated = false; p_scheme->InitializeNonLinIteration(BaseType::GetModelPart(), rA, rDx, rb); is_converged = mpConvergenceCriteria->PreCriteria(BaseType::GetModelPart(), p_builder_and_solver->GetDofSet(), rA, rDx, rb); - //function to perform the building and the solving phase. - if (BaseType::mRebuildLevel > 0 || BaseType::mStiffnessMatrixIsBuilt == false) - { + // Function to perform the building and the solving phase. + if (BaseType::mRebuildLevel > 0 || BaseType::mStiffnessMatrixIsBuilt == false) { TSparseSpace::SetToZero(rA); TSparseSpace::SetToZero(rDx); TSparseSpace::SetToZero(rb); p_builder_and_solver->BuildAndSolve(p_scheme, BaseType::GetModelPart(), rA, rDx, rb); - } - else - { + } else { TSparseSpace::SetToZero(rDx); //Dx=0.00; TSparseSpace::SetToZero(rb); @@ -698,13 +706,8 @@ class ResidualBasedNewtonRaphsonStrategy p_scheme->FinalizeNonLinIteration(BaseType::GetModelPart(), rA, rDx, rb); - if (is_converged == true) - { - //initialisation of the convergence criteria - mpConvergenceCriteria->InitializeSolutionStep(BaseType::GetModelPart(), p_builder_and_solver->GetDofSet(), rA, rDx, rb); - - if (mpConvergenceCriteria->GetActualizeRHSflag() == true) - { + if (is_converged) { + if (mpConvergenceCriteria->GetActualizeRHSflag()) { TSparseSpace::SetToZero(rb); p_builder_and_solver->BuildRHS(p_scheme, BaseType::GetModelPart(), rb);