diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp index e0bd1af64b6b657686b9168dd9eb05b2469ebbfe..696d4c02d88b2427521be9ce5f87f1dcc2e9e788 100644 --- a/src/language/modules/SchemeModule.cpp +++ b/src/language/modules/SchemeModule.cpp @@ -269,9 +269,11 @@ SchemeModule::SchemeModule() this->_addBuiltinFunction("coupling", std::function( - [](std::shared_ptr<const IBoundaryDescriptor> boundary) + [](std::shared_ptr<const IBoundaryDescriptor> boundary, + const size_t& label) -> std::shared_ptr<const IBoundaryConditionDescriptor> { - return std::make_shared<CouplingBoundaryConditionDescriptor>(boundary); + return std::make_shared<CouplingBoundaryConditionDescriptor>(boundary, + label); } )); diff --git a/src/scheme/CouplingBoundaryConditionDescriptor.hpp b/src/scheme/CouplingBoundaryConditionDescriptor.hpp index 3b4b26d320ab5318efee7ad95b00d34b53e78c67..30ab627d5d3119fbde47c6a6895dd7d8a7a8c4a5 100644 --- a/src/scheme/CouplingBoundaryConditionDescriptor.hpp +++ b/src/scheme/CouplingBoundaryConditionDescriptor.hpp @@ -17,8 +17,10 @@ class CouplingBoundaryConditionDescriptor : public IBoundaryConditionDescriptor } std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor; + size_t m_label; public: + const IBoundaryDescriptor& boundaryDescriptor() const final { @@ -31,8 +33,13 @@ class CouplingBoundaryConditionDescriptor : public IBoundaryConditionDescriptor return Type::coupling; } - CouplingBoundaryConditionDescriptor(std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor) - : m_boundary_descriptor(boundary_descriptor) + size_t + label() const{ + return m_label; + } + + CouplingBoundaryConditionDescriptor(std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor, size_t label) + : m_boundary_descriptor(boundary_descriptor), m_label{label} { ; } diff --git a/src/scheme/IBoundaryConditionDescriptor.hpp b/src/scheme/IBoundaryConditionDescriptor.hpp index 5092f1973a9988be5342c48ebbdf0445ecb9d5b4..54e714932b17ac9c35341f77aa886287905b7d54 100644 --- a/src/scheme/IBoundaryConditionDescriptor.hpp +++ b/src/scheme/IBoundaryConditionDescriptor.hpp @@ -34,7 +34,7 @@ class IBoundaryConditionDescriptor virtual const IBoundaryDescriptor& boundaryDescriptor() const = 0; virtual Type type() const = 0; - + IBoundaryConditionDescriptor() = default; IBoundaryConditionDescriptor(const IBoundaryConditionDescriptor&) = delete; IBoundaryConditionDescriptor(IBoundaryConditionDescriptor&&) = delete; diff --git a/src/scheme/LocalDtAcousticSolver.cpp b/src/scheme/LocalDtAcousticSolver.cpp index e535c6e8dae0eddb2ebbc7d06ebe7d36cb5e4f43..0e865976f56f718d89d9f746db16f71e316e6dc7 100644 --- a/src/scheme/LocalDtAcousticSolver.cpp +++ b/src/scheme/LocalDtAcousticSolver.cpp @@ -272,8 +272,10 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final break; } case IBoundaryConditionDescriptor::Type::coupling: { + const CouplingBoundaryConditionDescriptor& coupling_bc_descriptor = + dynamic_cast<const CouplingBoundaryConditionDescriptor&>(*bc_descriptor); bc_list.emplace_back( - CouplingBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()))); + CouplingBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),coupling_bc_descriptor.label())); break; } default: { @@ -862,62 +864,69 @@ void LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyCoupl std::visit([&](auto&& bc2) { using T2 = std::decay_t<decltype(bc2)>; if constexpr (std::is_same_v<CouplingBoundaryCondition,T2>){ - const auto& node_list2 = bc2.nodeList(); - if constexpr (Dimension == 1){ - for(size_t i = 0; i<node_list1.size(); i++){ - NodeId node_id1 = node_list1[i]; - NodeId node_id2; - for(size_t j = 0; j<node_list2.size(); j++){ - node_id2 = node_list2[j]; - if(abs(xr1[node_id1][0] - xr2[node_id2][0]) < 1e-14){ - Ar1[node_id1] += Ar2[node_id2]; - Ar2[node_id2] = Ar1[node_id1]; - - br1[node_id1] += br2[node_id2]; - br2[node_id2] = br1[node_id1]; - - break; + if(bc1.label() == bc2.label()){ + const auto& node_list2 = bc2.nodeList(); + if constexpr (Dimension == 1){ + for(size_t i = 0; i<node_list1.size(); i++){ + NodeId node_id1 = node_list1[i]; + NodeId node_id2; + for(size_t j = 0; j<node_list2.size(); j++){ + node_id2 = node_list2[j]; + if(abs(xr1[node_id1][0] - xr2[node_id2][0]) < 1e-13){ + Ar1[node_id1] += Ar2[node_id2]; + Ar2[node_id2] = Ar1[node_id1]; + + br1[node_id1] += br2[node_id2]; + br2[node_id2] = br1[node_id1]; + + break; + } } } - } - } else if constexpr (Dimension == 2){ - for(size_t i = 0; i<node_list1.size(); i++){ - NodeId node_id1 = node_list1[i]; - NodeId node_id2; - for(size_t j = 0; j<node_list2.size(); j++){ - node_id2 = node_list2[j]; - if((abs(xr1[node_id1][0] - xr2[node_id2][0]) < 1e-14) and - (abs(xr1[node_id1][1] - xr2[node_id2][1]) < 1e-14)){ - Ar1[node_id1] += Ar2[node_id2]; - Ar2[node_id2] = Ar1[node_id1]; - - br1[node_id1] += br2[node_id2]; - br2[node_id2] = br1[node_id1]; - - break; + } else if constexpr (Dimension == 2){ + for(size_t i = 0; i<node_list1.size(); i++){ + NodeId node_id1 = node_list1[i]; + NodeId node_id2; + for(size_t j = 0; j<node_list2.size(); j++){ + node_id2 = node_list2[j]; + double err = 0; + for(size_t k = 0; k<Dimension; k++){ + err += (xr1[node_id1][k] - xr2[node_id2][k])*(xr1[node_id1][k] - xr2[node_id2][k]); + } + if(sqrt(err) < 1e-14){ + Ar1[node_id1] += Ar2[node_id2]; + Ar2[node_id2] = Ar1[node_id1]; + + br1[node_id1] += br2[node_id2]; + br2[node_id2] = br1[node_id1]; + + break; + } } } - } - } else { - for(size_t i = 0; i<node_list1.size(); i++){ - NodeId node_id1 = node_list1[i]; - NodeId node_id2; - for(size_t j = 0; j<node_list2.size(); j++){ - node_id2 = node_list2[j]; - if((abs(xr1[node_id1][0] - xr2[node_id2][0]) < 1e-14) and - (abs(xr1[node_id1][1] - xr2[node_id2][1]) < 1e-14) and - (abs(xr1[node_id1][2] - xr2[node_id2][2]) < 1e-14)){ - Ar1[node_id1] += Ar2[node_id2]; - Ar2[node_id2] = Ar1[node_id1]; - - br1[node_id1] += br2[node_id2]; - br2[node_id2] = br1[node_id1]; - - break; + } else { + for(size_t i = 0; i<node_list1.size(); i++){ + NodeId node_id1 = node_list1[i]; + NodeId node_id2; + for(size_t j = 0; j<node_list2.size(); j++){ + node_id2 = node_list2[j]; + double err = 0; + for(size_t k = 0; k<Dimension; k++){ + err += (xr1[node_id1][k] - xr2[node_id2][k])*(xr1[node_id1][k] - xr2[node_id2][k]); + } + if(sqrt(err) < 1e-14){ + Ar1[node_id1] += Ar2[node_id2]; + Ar2[node_id2] = Ar1[node_id1]; + + br1[node_id1] += br2[node_id2]; + br2[node_id2] = br1[node_id1]; + + break; + } } } } - } + } } },boundary_condition2); } @@ -1155,6 +1164,7 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::CouplingBo { private: const MeshNodeBoundary<Dimension> m_mesh_node_boundary; + const size_t m_label; public: const Array<const NodeId>& @@ -1163,8 +1173,13 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::CouplingBo return m_mesh_node_boundary.nodeList(); } - CouplingBoundaryCondition(const MeshNodeBoundary<Dimension>& mesh_node_boundary) - : m_mesh_node_boundary{mesh_node_boundary} + size_t + label() const { + return m_label; + } + + CouplingBoundaryCondition(const MeshNodeBoundary<Dimension>& mesh_node_boundary,const size_t label) + : m_mesh_node_boundary{mesh_node_boundary}, m_label{label} { ; }