From 10351df6eddde592501c041d7b95ec3de7f771ed Mon Sep 17 00:00:00 2001
From: Compte local pour Alexandre Gangloff <Alexandre.GANGLOFF@cea.fr>
Date: Fri, 9 Jun 2023 15:53:00 +0200
Subject: [PATCH] add label for coupling boundary condition

---
 src/language/modules/SchemeModule.cpp         |   6 +-
 .../CouplingBoundaryConditionDescriptor.hpp   |  11 +-
 src/scheme/IBoundaryConditionDescriptor.hpp   |   2 +-
 src/scheme/LocalDtAcousticSolver.cpp          | 119 ++++++++++--------
 4 files changed, 81 insertions(+), 57 deletions(-)

diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index e0bd1af64..696d4c02d 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 3b4b26d32..30ab627d5 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 5092f1973..54e714932 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 e535c6e8d..0e865976f 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}
   {
     ;
   }
-- 
GitLab