Skip to content

Commit

Permalink
Merge pull request #7865 from KratosMultiphysics/mapping/coupling_geo…
Browse files Browse the repository at this point in the history 
…ms_choose_sides

[Mapping] Switch slave sides of coupling geom mapper
  • Loading branch information
@peterjwilson
peterjwilson authored Nov 24, 2020
2 parents 4ce824d + 1676df1 commit 0fdfdcd
Show file tree
Hide file tree
Showing 4 changed files with 100 additions and 43 deletions.
75 changes: 42 additions & 33 deletions applications/MappingApplication/custom_mappers/coupling_geometry_mapper.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -131,6 +134,7 @@ CouplingGeometryMapper<TSparseSpace, TDenseSpace>::CouplingGeometryMapper(
mMapperSettings(JsonParameters)
{
JsonParameters.ValidateAndAssignDefaults(GetMapperDefaultSettings());
const bool destination_is_slave = mMapperSettings["destination_is_slave"].GetBool();

mpModeler = (ModelerFactory::Create(
mMapperSettings["modeler_name"].GetString(),
Expand All @@ -139,17 +143,21 @@ CouplingGeometryMapper<TSparseSpace, TDenseSpace>::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<InterfaceVectorContainerType>(*mpCouplingInterfaceOrigin);
mpInterfaceVectorContainerDestination = Kratos::make_unique<InterfaceVectorContainerType>(*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<InterfaceVectorContainerType>(*mpCouplingInterfaceMaster);
mpInterfaceVectorContainerSlave = Kratos::make_unique<InterfaceVectorContainerType>(*mpCouplingInterfaceSlave);

this->CreateLinearSolver();
this->InitializeInterface();
Expand All @@ -162,8 +170,9 @@ void CouplingGeometryMapper<TSparseSpace, TDenseSpace>::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(),
Expand All @@ -176,26 +185,26 @@ void CouplingGeometryMapper<TSparseSpace, TDenseSpace>::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<MappingMatrixType>(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<TSparseSpace, TDenseSpace>(
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<TSparseSpace, TDenseSpace>(
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

Expand Down Expand Up @@ -225,7 +234,7 @@ void CouplingGeometryMapper<TSparseSpace, TDenseSpace>::InitializeInterface(Krat
SparseMatrixMultiplicationUtility::MatrixMultiplication(*mpMappingMatrixSlave, *mpMappingMatrixProjector, *mpMappingMatrix);
}
else {
MappingMatrixUtilities::InitializeSystemVector<TSparseSpace, TDenseSpace>(mpTempVector, mpInterfaceVectorContainerDestination->GetModelPart().NumberOfNodes());
MappingMatrixUtilities::InitializeSystemVector<TSparseSpace, TDenseSpace>(mpTempVector, mpInterfaceVectorContainerSlave->GetModelPart().NumberOfNodes());
if (precompute_mapping_matrix) CalculateMappingMatrixWithSolver(*mpMappingMatrixSlave, *mpMappingMatrixProjector);
}

Expand All @@ -246,23 +255,23 @@ void CouplingGeometryMapper<TSparseSpace, TDenseSpace>::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<class TSparseSpace, class TDenseSpace>
Expand All @@ -274,23 +283,23 @@ void CouplingGeometryMapper<TSparseSpace, TDenseSpace>::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<class TSparseSpace, class TDenseSpace>
Expand Down
34 changes: 25 additions & 9 deletions applications/MappingApplication/custom_mappers/coupling_geometry_mapper.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -62,7 +64,7 @@ class CouplingGeometryLocalSystem : public MapperLocalSystem

MapperLocalSystemUniquePointer Create(GeometryPointerType pGeometry) const override
{
return Kratos::make_unique<CouplingGeometryLocalSystem>(pGeometry, mIsProjection, mIsDualMortar);
return Kratos::make_unique<CouplingGeometryLocalSystem>(pGeometry, mIsProjection, mIsDualMortar, mIsDestinationIsSlave);
}

/// Turn back information as a string.
Expand All @@ -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 TSparseSpace, class TDenseSpace>
class CouplingGeometryMapper : public Mapper<TSparseSpace, TDenseSpace>
{
Expand Down Expand Up @@ -259,8 +274,8 @@ class CouplingGeometryMapper : public Mapper<TSparseSpace, TDenseSpace>
ModelPart& mrModelPartOrigin;
ModelPart& mrModelPartDestination;
ModelPart* mpCouplingMP = nullptr;
ModelPart* mpCouplingInterfaceOrigin = nullptr;
ModelPart* mpCouplingInterfaceDestination = nullptr;
ModelPart* mpCouplingInterfaceMaster = nullptr;
ModelPart* mpCouplingInterfaceSlave = nullptr;

Parameters mMapperSettings;

Expand All @@ -275,8 +290,8 @@ class CouplingGeometryMapper : public Mapper<TSparseSpace, TDenseSpace>
MapperLocalSystemPointerVector mMapperLocalSystemsProjector;
MapperLocalSystemPointerVector mMapperLocalSystemsSlave;

InterfaceVectorContainerPointerType mpInterfaceVectorContainerOrigin;
InterfaceVectorContainerPointerType mpInterfaceVectorContainerDestination;
InterfaceVectorContainerPointerType mpInterfaceVectorContainerMaster;
InterfaceVectorContainerPointerType mpInterfaceVectorContainerSlave;

LinearSolverSharedPointerType mpLinearSolver = nullptr;

Expand All @@ -285,8 +300,8 @@ class CouplingGeometryMapper : public Mapper<TSparseSpace, TDenseSpace>

void AssignInterfaceEquationIds()
{
MapperUtilities::AssignInterfaceEquationIds(mpCouplingInterfaceDestination->GetCommunicator());
MapperUtilities::AssignInterfaceEquationIds(mpCouplingInterfaceOrigin->GetCommunicator());
MapperUtilities::AssignInterfaceEquationIds(mpCouplingInterfaceSlave->GetCommunicator());
MapperUtilities::AssignInterfaceEquationIds(mpCouplingInterfaceMaster->GetCommunicator());
}

void MapInternal(const Variable<double>& rOriginVariable,
Expand Down Expand Up @@ -324,6 +339,7 @@ class CouplingGeometryMapper : public Mapper<TSparseSpace, TDenseSpace>
"modeler_parameters" : {},
"consistency_scaling" : true,
"row_sum_tolerance" : 1e-12,
"destination_is_slave" : true,
"linear_solver_settings" : {}
})");
}
Expand Down
1 change: 1 addition & 0 deletions applications/MappingApplication/tests/test_MappingApplication.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
33 changes: 32 additions & 1 deletion applications/MappingApplication/tests/test_coupling_geometry_mapper.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)


Expand Down Expand Up @@ -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
Expand Down

0 comments on commit 0fdfdcd

Please sign in to comment.