diff --git a/applications/MappingApplication/custom_mappers/coupling_geometry_mapper.cpp b/applications/MappingApplication/custom_mappers/coupling_geometry_mapper.cpp index 44f6403c8002..7c4fbc72cc2d 100644 --- a/applications/MappingApplication/custom_mappers/coupling_geometry_mapper.cpp +++ b/applications/MappingApplication/custom_mappers/coupling_geometry_mapper.cpp @@ -32,10 +32,13 @@ void CouplingGeometryLocalSystem::CalculateAll(MatrixType& rLocalMappingMatrix, EquationIdVectorType& rDestinationIds, MapperLocalSystem::PairingStatus& rPairingStatus) const { + const IndexType slave_index = (mIsDestinationIsSlave) ? 1 : 0; + const IndexType master_index = 1 - slave_index; + const auto& r_geometry_master = (mIsProjection) - ? mpGeom->GetGeometryPart(0) // set to master - get projected 'mass' matrix - : mpGeom->GetGeometryPart(1); // set to slave - get consistent slave 'mass' matrix - const auto& r_geometry_slave = mpGeom->GetGeometryPart(1); + ? mpGeom->GetGeometryPart(master_index) // set to master - get projected 'mass' matrix + : mpGeom->GetGeometryPart(slave_index); // set to slave - get consistent slave 'mass' matrix + const auto& r_geometry_slave = mpGeom->GetGeometryPart(slave_index); const bool is_dual_mortar = (!mIsProjection && mIsDualMortar) ? true @@ -131,6 +134,7 @@ CouplingGeometryMapper::CouplingGeometryMapper( mMapperSettings(JsonParameters) { JsonParameters.ValidateAndAssignDefaults(GetMapperDefaultSettings()); + const bool destination_is_slave = mMapperSettings["destination_is_slave"].GetBool(); mpModeler = (ModelerFactory::Create( mMapperSettings["modeler_name"].GetString(), @@ -139,17 +143,21 @@ CouplingGeometryMapper::CouplingGeometryMapper( // adds destination model part mpModeler->GenerateNodes(rModelPartDestination); - mpModeler->SetupGeometryModel(); mpModeler->PrepareGeometryModel(); // here use whatever ModelPart(s) was created by the Modeler mpCouplingMP = &(rModelPartOrigin.GetModel().GetModelPart("coupling")); - mpCouplingInterfaceOrigin = mpCouplingMP->pGetSubModelPart("interface_origin"); - mpCouplingInterfaceDestination = mpCouplingMP->pGetSubModelPart("interface_destination"); - mpInterfaceVectorContainerOrigin = Kratos::make_unique(*mpCouplingInterfaceOrigin); - mpInterfaceVectorContainerDestination = Kratos::make_unique(*mpCouplingInterfaceDestination); + mpCouplingInterfaceMaster = (destination_is_slave) + ? mpCouplingMP->pGetSubModelPart("interface_origin") + : mpCouplingMP->pGetSubModelPart("interface_destination"); + mpCouplingInterfaceSlave = (destination_is_slave) + ? mpCouplingMP->pGetSubModelPart("interface_destination") + : mpCouplingMP->pGetSubModelPart("interface_origin"); + + mpInterfaceVectorContainerMaster = Kratos::make_unique(*mpCouplingInterfaceMaster); + mpInterfaceVectorContainerSlave = Kratos::make_unique(*mpCouplingInterfaceSlave); this->CreateLinearSolver(); this->InitializeInterface(); @@ -162,8 +170,9 @@ void CouplingGeometryMapper::InitializeInterface(Krat // compose local element mappings const bool dual_mortar = mMapperSettings["dual_mortar"].GetBool(); const bool precompute_mapping_matrix = mMapperSettings["precompute_mapping_matrix"].GetBool(); - CouplingGeometryLocalSystem ref_projector_local_system(nullptr, true, dual_mortar); - CouplingGeometryLocalSystem ref_slave_local_system(nullptr, false, dual_mortar); + const bool direct_map_to_destination = mMapperSettings["destination_is_slave"].GetBool(); + CouplingGeometryLocalSystem ref_projector_local_system(nullptr, true, dual_mortar, direct_map_to_destination); + CouplingGeometryLocalSystem ref_slave_local_system(nullptr, false, dual_mortar, direct_map_to_destination); MapperUtilities::CreateMapperLocalSystemsFromGeometries(ref_projector_local_system, mpCouplingMP->GetCommunicator(), @@ -176,26 +185,26 @@ void CouplingGeometryMapper::InitializeInterface(Krat AssignInterfaceEquationIds(); // Has to be done every time in case of overlapping interfaces! // assemble projector interface mass matrix - interface_matrix_projector - const std::size_t num_nodes_interface_slave = mpCouplingInterfaceDestination->NumberOfNodes(); - const std::size_t num_nodes_interface_master = mpCouplingInterfaceOrigin->NumberOfNodes(); + const std::size_t num_nodes_interface_slave = mpCouplingInterfaceSlave->NumberOfNodes(); + const std::size_t num_nodes_interface_master = mpCouplingInterfaceMaster->NumberOfNodes(); mpMappingMatrix = Kratos::make_unique(num_nodes_interface_slave, num_nodes_interface_master); // TODO Philipp I am pretty sure we should separate the vector construction from the matrix construction, should be independent otherwise no clue what is happening MappingMatrixUtilities::BuildMappingMatrix( mpMappingMatrixSlave, - mpInterfaceVectorContainerDestination->pGetVector(), - mpInterfaceVectorContainerDestination->pGetVector(), - mpInterfaceVectorContainerDestination->GetModelPart(), - mpInterfaceVectorContainerDestination->GetModelPart(), + mpInterfaceVectorContainerSlave->pGetVector(), + mpInterfaceVectorContainerSlave->pGetVector(), + mpInterfaceVectorContainerSlave->GetModelPart(), + mpInterfaceVectorContainerSlave->GetModelPart(), mMapperLocalSystemsSlave, 0); // The echo-level is no longer neeed here, refactor in separate PR MappingMatrixUtilities::BuildMappingMatrix( mpMappingMatrixProjector, - mpInterfaceVectorContainerOrigin->pGetVector(), - mpInterfaceVectorContainerDestination->pGetVector(), - mpInterfaceVectorContainerOrigin->GetModelPart(), - mpInterfaceVectorContainerDestination->GetModelPart(), + mpInterfaceVectorContainerMaster->pGetVector(), + mpInterfaceVectorContainerSlave->pGetVector(), + mpInterfaceVectorContainerMaster->GetModelPart(), + mpInterfaceVectorContainerSlave->GetModelPart(), mMapperLocalSystemsProjector, 0); // The echo-level is no longer neeed here, refactor in separate PR @@ -225,7 +234,7 @@ void CouplingGeometryMapper::InitializeInterface(Krat SparseMatrixMultiplicationUtility::MatrixMultiplication(*mpMappingMatrixSlave, *mpMappingMatrixProjector, *mpMappingMatrix); } else { - MappingMatrixUtilities::InitializeSystemVector(mpTempVector, mpInterfaceVectorContainerDestination->GetModelPart().NumberOfNodes()); + MappingMatrixUtilities::InitializeSystemVector(mpTempVector, mpInterfaceVectorContainerSlave->GetModelPart().NumberOfNodes()); if (precompute_mapping_matrix) CalculateMappingMatrixWithSolver(*mpMappingMatrixSlave, *mpMappingMatrixProjector); } @@ -246,23 +255,23 @@ void CouplingGeometryMapper::MapInternal( const bool dual_mortar = mMapperSettings["dual_mortar"].GetBool(); const bool precompute_mapping_matrix = mMapperSettings["precompute_mapping_matrix"].GetBool(); - mpInterfaceVectorContainerOrigin->UpdateSystemVectorFromModelPart(rOriginVariable, MappingOptions); + mpInterfaceVectorContainerMaster->UpdateSystemVectorFromModelPart(rOriginVariable, MappingOptions); if (dual_mortar || precompute_mapping_matrix) { TSparseSpace::Mult( *mpMappingMatrix, - mpInterfaceVectorContainerOrigin->GetVector(), - mpInterfaceVectorContainerDestination->GetVector()); // rQd = rMdo * rQo + mpInterfaceVectorContainerMaster->GetVector(), + mpInterfaceVectorContainerSlave->GetVector()); // rQd = rMdo * rQo } else { TSparseSpace::Mult( *mpMappingMatrixProjector, - mpInterfaceVectorContainerOrigin->GetVector(), + mpInterfaceVectorContainerMaster->GetVector(), *mpTempVector); // rQd = rMdo * rQo - mpLinearSolver->Solve(*mpMappingMatrixSlave, mpInterfaceVectorContainerDestination->GetVector(), *mpTempVector); + mpLinearSolver->Solve(*mpMappingMatrixSlave, mpInterfaceVectorContainerSlave->GetVector(), *mpTempVector); } - mpInterfaceVectorContainerDestination->UpdateModelPartFromSystemVector(rDestinationVariable, MappingOptions); + mpInterfaceVectorContainerSlave->UpdateModelPartFromSystemVector(rDestinationVariable, MappingOptions); } template @@ -274,23 +283,23 @@ void CouplingGeometryMapper::MapInternalTranspose( const bool dual_mortar = mMapperSettings["dual_mortar"].GetBool(); const bool precompute_mapping_matrix = mMapperSettings["precompute_mapping_matrix"].GetBool(); - mpInterfaceVectorContainerDestination->UpdateSystemVectorFromModelPart(rDestinationVariable, MappingOptions); + mpInterfaceVectorContainerSlave->UpdateSystemVectorFromModelPart(rDestinationVariable, MappingOptions); if (dual_mortar || precompute_mapping_matrix) { TSparseSpace::TransposeMult( *mpMappingMatrix, - mpInterfaceVectorContainerDestination->GetVector(), - mpInterfaceVectorContainerOrigin->GetVector()); // rQo = rMdo^T * rQd + mpInterfaceVectorContainerSlave->GetVector(), + mpInterfaceVectorContainerMaster->GetVector()); // rQo = rMdo^T * rQd } else { - mpLinearSolver->Solve(*mpMappingMatrixSlave, *mpTempVector, mpInterfaceVectorContainerDestination->GetVector()); + mpLinearSolver->Solve(*mpMappingMatrixSlave, *mpTempVector, mpInterfaceVectorContainerSlave->GetVector()); TSparseSpace::TransposeMult( *mpMappingMatrixProjector, *mpTempVector, - mpInterfaceVectorContainerOrigin->GetVector()); // rQo = rMdo^T * rQd + mpInterfaceVectorContainerMaster->GetVector()); // rQo = rMdo^T * rQd } - mpInterfaceVectorContainerOrigin->UpdateModelPartFromSystemVector(rOriginVariable, MappingOptions); + mpInterfaceVectorContainerMaster->UpdateModelPartFromSystemVector(rOriginVariable, MappingOptions); } template diff --git a/applications/MappingApplication/custom_mappers/coupling_geometry_mapper.h b/applications/MappingApplication/custom_mappers/coupling_geometry_mapper.h index 6f2bfe345c4e..7f330acc76f9 100644 --- a/applications/MappingApplication/custom_mappers/coupling_geometry_mapper.h +++ b/applications/MappingApplication/custom_mappers/coupling_geometry_mapper.h @@ -41,11 +41,13 @@ class CouplingGeometryLocalSystem : public MapperLocalSystem explicit CouplingGeometryLocalSystem(GeometryPointerType pGeom, const bool IsProjection, - const bool IsDualMortar + const bool IsDualMortar, + const bool IsDestinationIsSlave ) : mpGeom(pGeom), mIsProjection(IsProjection), - mIsDualMortar(IsDualMortar) + mIsDualMortar(IsDualMortar), + mIsDestinationIsSlave(IsDestinationIsSlave) {} void CalculateAll(MatrixType& rLocalMappingMatrix, @@ -62,7 +64,7 @@ class CouplingGeometryLocalSystem : public MapperLocalSystem MapperLocalSystemUniquePointer Create(GeometryPointerType pGeometry) const override { - return Kratos::make_unique(pGeometry, mIsProjection, mIsDualMortar); + return Kratos::make_unique(pGeometry, mIsProjection, mIsDualMortar, mIsDestinationIsSlave); } /// Turn back information as a string. @@ -73,9 +75,22 @@ class CouplingGeometryLocalSystem : public MapperLocalSystem bool mIsProjection; // Set to true is we are projecting the master onto the slave. // Set to false if we are projecting the slave onto the slave. bool mIsDualMortar = false; + bool mIsDestinationIsSlave = true; }; +// CouplingGeometryMapper +// +// The mapper always forward maps from the master to the slave. +// Normally: +// master = interface origin +// slave = interface destination +// +// However, this can be reversed by setting 'destination_is_slave' = false. +// This yields: +// master = interface destination +// slave = interface origin + template class CouplingGeometryMapper : public Mapper { @@ -259,8 +274,8 @@ class CouplingGeometryMapper : public Mapper ModelPart& mrModelPartOrigin; ModelPart& mrModelPartDestination; ModelPart* mpCouplingMP = nullptr; - ModelPart* mpCouplingInterfaceOrigin = nullptr; - ModelPart* mpCouplingInterfaceDestination = nullptr; + ModelPart* mpCouplingInterfaceMaster = nullptr; + ModelPart* mpCouplingInterfaceSlave = nullptr; Parameters mMapperSettings; @@ -275,8 +290,8 @@ class CouplingGeometryMapper : public Mapper MapperLocalSystemPointerVector mMapperLocalSystemsProjector; MapperLocalSystemPointerVector mMapperLocalSystemsSlave; - InterfaceVectorContainerPointerType mpInterfaceVectorContainerOrigin; - InterfaceVectorContainerPointerType mpInterfaceVectorContainerDestination; + InterfaceVectorContainerPointerType mpInterfaceVectorContainerMaster; + InterfaceVectorContainerPointerType mpInterfaceVectorContainerSlave; LinearSolverSharedPointerType mpLinearSolver = nullptr; @@ -285,8 +300,8 @@ class CouplingGeometryMapper : public Mapper void AssignInterfaceEquationIds() { - MapperUtilities::AssignInterfaceEquationIds(mpCouplingInterfaceDestination->GetCommunicator()); - MapperUtilities::AssignInterfaceEquationIds(mpCouplingInterfaceOrigin->GetCommunicator()); + MapperUtilities::AssignInterfaceEquationIds(mpCouplingInterfaceSlave->GetCommunicator()); + MapperUtilities::AssignInterfaceEquationIds(mpCouplingInterfaceMaster->GetCommunicator()); } void MapInternal(const Variable& rOriginVariable, @@ -324,6 +339,7 @@ class CouplingGeometryMapper : public Mapper "modeler_parameters" : {}, "consistency_scaling" : true, "row_sum_tolerance" : 1e-12, + "destination_is_slave" : true, "linear_solver_settings" : {} })"); } diff --git a/applications/MappingApplication/tests/test_MappingApplication.py b/applications/MappingApplication/tests/test_MappingApplication.py index bc1a84c70be4..30a0ea1fab5e 100644 --- a/applications/MappingApplication/tests/test_MappingApplication.py +++ b/applications/MappingApplication/tests/test_MappingApplication.py @@ -56,6 +56,7 @@ def AssembleTestSuites(): nightSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_coupling_geometry_mapper.TestCouplingGeometryMapper])) nightSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_coupling_geometry_mapper.TestDualMortarCouplingGeometryMapper])) + nightSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_coupling_geometry_mapper.TestSlaveOriginCouplingGeometryMapper])) nightSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_coupling_geometry_mapper.TestComputeMappingMatrixCouplingGeometryMapper])) # Create a test suit that contains all the tests from every testCase diff --git a/applications/MappingApplication/tests/test_coupling_geometry_mapper.py b/applications/MappingApplication/tests/test_coupling_geometry_mapper.py index cfe331270671..bf5d05eef2cb 100644 --- a/applications/MappingApplication/tests/test_coupling_geometry_mapper.py +++ b/applications/MappingApplication/tests/test_coupling_geometry_mapper.py @@ -27,7 +27,7 @@ def test_inverse_map_forces(self): SetConstantVariable(self.interface_model_part_destination,KM.FORCE,reference_force) self.mapper.InverseMap(KM.FORCE, KM.FORCE,KratosMapping.Mapper.USE_TRANSPOSE) mapped_results = GetInterfaceResult(self.interface_model_part_origin,KM.FORCE) - reference_result = [0.2501913843336516, 0.2501913843336516, 0.2501913843336516, 1.292449164738119, 1.292449164738119, 1.292449164738119, 0.7004426959773076, 0.7004426959773076, 0.7004426959773076, 0.8902247629970814, 0.8902247629970814, 0.8902247629970814, 0.9474688054530651, 0.9474688054530651, 0.9474688054530651, 0.919223186500775, 0.919223186500775, 0.919223186500775] + reference_result = [0.2380991480071958, 0.2380991480071958, 0.2380991480071958, 1.3120351229689677, 1.3120351229689677, 1.3120351229689677, 0.6908309106360845, 0.6908309106360845, 0.6908309106360845, 0.9063686826513201, 0.9063686826513201, 0.9063686826513201, 0.9261336708771284, 0.9261336708771284, 0.9261336708771284, 0.9265324648593039, 0.9265324648593039, 0.9265324648593039] self.assertVectorAlmostEqual(mapped_results,reference_result) @@ -61,6 +61,37 @@ def test_dual_mortar(self): reference_result = [1.0, 1.0, 1.0, 0.9999999999999999, 0.9999999999999999, 0.9999999999999999, 1.0000000000000004, 1.0000000000000004, 1.0000000000000004, 0.9999999999999998, 0.9999999999999998, 0.9999999999999998, 1.0000000000000004, 1.0000000000000004, 1.0000000000000004] self.assertVectorAlmostEqual(mapped_results,reference_result) +class TestSlaveOriginCouplingGeometryMapper(KratosUnittest.TestCase): + @classmethod + def setUpClass(self): + self.mapper_parameters = KM.Parameters("""{ + "mapper_type": "coupling_geometry", + "echo_level" : 0, + "precompute_mapping_matrix" : false, + "dual_mortar": false, + "consistency_scaling" : true, + "destination_is_slave" : false, + "modeler_name" : "MappingGeometriesModeler", + "modeler_parameters":{ + "origin_model_part_name" : "origin", + "destination_model_part_name" : "destination", + "is_interface_sub_model_parts_specified" : true, + "origin_interface_sub_model_part_name" : "origin.line_tri", + "destination_interface_sub_model_part_name" : "destination.line_quad" + } + }""") + + SetupModelParts(self) + CreateMapper(self) + + def test_slave_origin_mortar(self): + reference_displacement = 1.0 + SetConstantVariable(self.interface_model_part_destination,KM.DISPLACEMENT,reference_displacement) + self.mapper.Map(KM.DISPLACEMENT, KM.DISPLACEMENT) + mapped_results = GetInterfaceResult(self.interface_model_part_origin,KM.DISPLACEMENT) + reference_result = [0.9999999999999998, 0.9999999999999998, 0.9999999999999998, 1.0000000000000002, 1.0000000000000002, 1.0000000000000002, 0.9999999999999998, 0.9999999999999998, 0.9999999999999998, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] + self.assertVectorAlmostEqual(mapped_results,reference_result) + class TestComputeMappingMatrixCouplingGeometryMapper(KratosUnittest.TestCase): @classmethod