Make the Point Coupling dof based not node based

Now each degree of freedom is considered separately

Author: ischeider

src/core/fem/src/dofset/4C_fem_dofset.cpp

@@ -15,6 +15,7 @@
 #include <Epetra_FECrsGraph.h>
 
 #include <algorithm>
+#include <format>
 #include <iostream>
 
 FOUR_C_NAMESPACE_OPEN
@@ -214,87 +215,102 @@ int Core::DOFSets::DofSet::assign_degrees_of_freedom(
     if ((int)couplingconditions.size() > 0) pccdofhandling_ = true;
 
     // do the nodes first
-    Core::LinAlg::Vector<int> numdfrownodes(*dis.node_row_map());
+    Core::LinAlg::Vector<int> num_dof_rownodes(*dis.node_row_map());
     Core::LinAlg::Vector<int> idxrownodes(*dis.node_row_map());
 
     int numrownodes = dis.num_my_row_nodes();
     for (int i = 0; i < numrownodes; ++i)
     {
       Core::Nodes::Node* actnode = dis.l_row_node(i);
-      numdfrownodes[i] = num_dof_per_node(*actnode);
+      num_dof_rownodes[i] = num_dof_per_node(*actnode);
     }
 
     int minnodegid = get_minimal_node_gid_if_relevant(dis);
-    maxnodenumdf = numdfrownodes.max_value();
+    maxnodenumdf = num_dof_rownodes.max_value();
     get_reserved_max_num_dofper_node(maxnodenumdf);  // XFEM::XFEMDofSet set to const number!
 
+    std::vector<std::vector<int>> onoffcond;  // vector of onoff status for each condition
+    std::vector<int> numdofcond;              // vector of NUMDOF for each condition
+    std::vector<std::vector<int>> nodeids;    // vector of node IDs for each condition
+    std::vector<int> mgids;                   // vector of master node IDs for each condition
+
+    for (int k = 0; k < (int)couplingconditions.size(); ++k)
+    {
+      onoffcond.push_back(couplingconditions[k]->parameters().get<std::vector<int>>("ONOFF"));
+      numdofcond.push_back(couplingconditions[k]->parameters().get<int>("NUMDOF"));
+      nodeids.push_back(*couplingconditions[k]->get_nodes());
+      mgids.push_back(nodeids[k][0]);
+      std::cout << "Coupling condition " << k << " contains " << nodeids[k].size()
+                << " nodes, master: " << mgids[k] << "\n";
+
+      // check if all nodes in this condition are on same processor
+      // (otherwise throw a FOUR_C_THROW for now - not yet implemented)
+      bool allononeproc = true;
+      for (auto nd : nodeids[k])
+      {
+        if (!dis.node_row_map()->MyGID(nd)) allononeproc = false;
+      }
+      if (!allononeproc)
+        FOUR_C_THROW(
+            "ERROR: Nodes in point coupling condition must all be on same processor (for now)");
+    }
+
     for (int i = 0; i < numrownodes; ++i)
     {
-      Core::Nodes::Node* actnode = dis.l_row_node(i);
-      const int gid = actnode->id();
+      const int gid = dis.l_row_node(i)->id();
 
       // **********************************************************************
       // **********************************************************************
       // check for DoF coupling conditions                         popp 02/2016
       // **********************************************************************
       // **********************************************************************
-      int relevantcondid = -1;
+      std::vector<int> applied_condition(maxnodenumdf, -1);
       if (dspos_ == 0)
       {
-        for (int k = 0; k < (int)couplingconditions.size(); ++k)
+        for (int icond = 0; icond < (int)couplingconditions.size(); ++icond)
         {
-          if (couplingconditions[k]->contains_node(gid))
+          // check total number of dofs and determine which dofs are to be coupled
+          if (numdofcond[icond] != num_dof_rownodes[i])
+            FOUR_C_THROW(
+                "ERROR: Number of DoFs in coupling condition (%i) does not match node (%i)",
+                numdofcond[icond], num_dof_rownodes[i]);
+
+          // check if the node i is contained in the coupling condition
+          auto found = std::find(nodeids[icond].begin(), nodeids[icond].end(), gid);
+          if (found != nodeids[icond].end())
           {
-            if (relevantcondid != -1) FOUR_C_THROW("ERROR: Two coupling conditions on one node");
-            relevantcondid = k;
+            for (int k_on = 0; k_on < (int)onoffcond[icond].size(); ++k_on)
+            {
+              if (onoffcond[icond][k_on] == 0) continue;
+              FOUR_C_ASSERT(applied_condition[k_on] == -1,
+                  "ERROR: Two coupling conditions on the same degree of freedom");
+              applied_condition[k_on] = icond;
+            }
           }
         }
       }
-
       // check for node coupling condition and slave/master status
-      bool specialtreatment = false;
-      if (relevantcondid >= 0)
-      {
-        const std::vector<int>* nodeids = couplingconditions[relevantcondid]->get_nodes();
-        if (!nodeids) FOUR_C_THROW("ERROR: Condition does not have Node Ids");
+      bool specialtreatment(false);
 
-        // check if all nodes in this condition are on same processor
-        // (otherwise throw a FOUR_C_THROW for now - not yet implemented)
-        bool allononeproc = true;
-        for (int k = 0; k < (int)(nodeids->size()); ++k)
+      // for (int idim = 0; idim < (int)applied_condition.size(); ++idim);
+      if (*std::ranges::max_element(applied_condition.begin(), applied_condition.end()) >= 0)
+      {
+        // do nothing if the node is the master in all of the applied conditions
+        // do something for second, third, ... (slave) node
+        // also if the node is a master in one condition, but a slave in another
+        bool allmaster = true;
+        for (auto condition_id : applied_condition)
         {
-          int checkgid = (*nodeids)[k];
-          if (!dis.node_row_map()->MyGID(checkgid)) allononeproc = false;
+          if (condition_id == -1) continue;
+          if (mgids[condition_id] != gid) allmaster = false;
         }
-        if (!allononeproc)
-          FOUR_C_THROW(
-              "ERROR: Nodes in point coupling condition must all be on same processor (for now)");
-
-        // do nothing for first (master) node in coupling condition
-        // do something for second, third, ... (slave) node
-        if ((*nodeids)[0] != gid)
+        if (!allmaster)
         {
           // critical case
           specialtreatment = true;
 
-          // check total number of dofs and determine which dofs are to be coupled
-          if (couplingconditions[relevantcondid]->parameters().get<int>("NUMDOF") !=
-              numdfrownodes[i])
-            FOUR_C_THROW(
-                "ERROR: Number of DoFs in coupling condition ({}) does not match node ({})",
-                couplingconditions[relevantcondid]->parameters().get<int>("NUMDOF"),
-                numdfrownodes[i]);
-          const std::vector<int>& onoffcond =
-              couplingconditions[relevantcondid]->parameters().get<std::vector<int>>("ONOFF");
-
-          // get master node of this condition
-          int mgid = (*nodeids)[0];
-          std::vector<int>& mdofs = nodedofset[mgid];
-          if ((int)(mdofs.size()) == 0)
-            FOUR_C_THROW("ERROR: Master node has not yet been initialized with DoFs");
-
           // special treatment
-          int numdf = numdfrownodes[i];
+          int numdf = num_dof_rownodes[i];
           int dof = count + (gid - minnodegid) * maxnodenumdf;
           idxrownodes[i] = dof;
           std::vector<int>& dofs = nodedofset[gid];
@@ -303,9 +319,11 @@ int Core::DOFSets::DofSet::assign_degrees_of_freedom(
           duplicatedofs.reserve(numdf);
           for (int j = 0; j < numdf; ++j)
           {
-            // push back master node DoF ID if coupled
-            if (onoffcond[j] == 1)
+            // push back master node DoF ID if the id is coupled
+            // it might be that the node is the master in one direction and a slave in another
+            if (applied_condition[j] >= 0 && mgids[applied_condition[j]] != gid)
             {
+              std::vector<int>& mdofs = nodedofset[mgids[applied_condition[j]]];
               dofs.push_back(mdofs[j]);
               duplicatedofs.push_back(1);
             }
@@ -322,7 +340,8 @@ int Core::DOFSets::DofSet::assign_degrees_of_freedom(
       // standard treatment for non-coupling nodes and master coupling nodes
       if (!specialtreatment)
       {
-        int numdf = numdfrownodes[i];
+        // now treat only the nodes, which are master and thus did not get the special treatment
+        int numdf = num_dof_rownodes[i];
         int dof = count + (gid - minnodegid) * maxnodenumdf;
         idxrownodes[i] = dof;
         std::vector<int>& dofs = nodedofset[gid];
@@ -339,11 +358,28 @@ int Core::DOFSets::DofSet::assign_degrees_of_freedom(
       // **********************************************************************
       // **********************************************************************
       // **********************************************************************
+      // std::cout << "Node ID " << gid << " - 0:" << nodedofset[gid][0] << "/"
+      // << nodeduplicatedofset[gid][0] << ", 2:" << nodedofset[gid][2] << "/"
+      // << nodeduplicatedofset[gid][2]<< "\n";
+    }
+    for (int i = 0; i < numrownodes; ++i)
+    {
+      const int gid = dis.l_row_node(i)->id();
+      std::string s = std::format("Node {}: ", gid + 1);
+      bool printflag = false;
+      for (int j = 0; j < 3; ++j)
+      {
+        const int dupnode = nodeduplicatedofset[gid][j];
+        s += std::format("{}/{}, ", nodedofset[gid][j], dupnode);
+        if (dupnode == 1) printflag = true;
+      }
+      if (printflag) std::cout << s << std::endl;
     }
 
-    Epetra_Import nodeimporter(numdfcolnodes_->get_block_map(), numdfrownodes.get_block_map());
-    int err = numdfcolnodes_->import(numdfrownodes, nodeimporter, Insert);
-    if (err) FOUR_C_THROW("Import using importer returned err={}", err);
+
+    Epetra_Import nodeimporter(numdfcolnodes_->get_map(), num_dof_rownodes.get_map());
+    int err = numdfcolnodes_->import(num_dof_rownodes, nodeimporter, Insert);
+    if (err) FOUR_C_THROW("Import using importer returned err=%d", err);
     err = idxcolnodes_->import(idxrownodes, nodeimporter, Insert);
     if (err) FOUR_C_THROW("Import using importer returned err={}", err);
 
@@ -443,7 +479,7 @@ int Core::DOFSets::DofSet::assign_degrees_of_freedom(
     }
 
     Epetra_Import elementimporter(
-        numdfcolelements_->get_block_map(), numdfrowelements.get_block_map());
+        numdfcolelements_->get_map(), numdfrowelements.get_map());
     err = numdfcolelements_->import(numdfrowelements, elementimporter, Insert);
     if (err) FOUR_C_THROW("Import using importer returned err={}", err);
     err = idxcolelements_->import(idxrowelements, elementimporter, Insert);