Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

[Mapping] Switch slave sides of coupling geom mapper #7865

Merged
merged 3 commits into from
Nov 24, 2020
Merged

[Mapping] Switch slave sides of coupling geom mapper #7865

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
b075e32
add option for destination to be slave or master
peterjwilson Nov 24, 2020
549dd52
add slave origin test
peterjwilson Nov 24, 2020
1676df1
correct inverse force mapping test due to updated integration weights
peterjwilson Nov 24, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
75 changes: 42 additions & 33 deletions applications/MappingApplication/custom_mappers/coupling_geometry_mapper.cpp
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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