Correction of #23

Created by: kliegeois

This pull request resolves the bug of the issue #23.

The memory issue was either an issue at the construction or an issue at the destruction of the contacts_list class.

The issue at the construction has been solved not using anymore std::vector of std::vector that generates random issues with gcc (I do not know why).

The issue at the destruction has been solved by allocating the correct size to the contacts_normal Kokkos::view (n \times 6 instead of n \times 3 as done previously).

@etphipp this pull request solves the issue related to the email, "Double free memory".

The tests have been tested both with serial ctest and parallel ctest for both release and debug build and several times for each (from 6 to 10 runs to be sure to catch any random behavior) and none of them fails.

This PR includes a commit related to tbox that had an include error that prevents from building with icc.

mrstlnos/src/Contacts_list.cpp

@@ -16,227 +16,251 @@ using namespace discrete_problem;
 
 Contacts_list::Contacts_list()
 {
-    contacts_number =0;
-    contacts_slave_nodes_size = 0;
-    contacts_slave_elems_size = 0;
-    contacts_master_nodes_size = 0;
-    contacts_master_elems_size = 0;
-    contacts_slave_normal = 0;
-    contacts_master_normal = 0;
-    contacts_slave_nodes = 0;
-    contacts_slave_elems = 0;
-    contacts_master_nodes = 0;
-    contacts_master_elems = 0;
-    //contacts_dg = 0;
-
-    contacts = Kokkos::View<int**, Kokkos::LayoutRight> ("R",contacts_number,0);
+  contacts_number =0;
+  contacts_slave_nodes_size = 0;
+  contacts_slave_elems_size = 0;
+  contacts_master_nodes_size = 0;
+  contacts_master_elems_size = 0;
+  contacts_slave_normal = 0;
+  contacts_master_normal = 0;
+  contacts_slave_nodes = 0;
+  contacts_slave_elems = 0;
+  contacts_master_nodes = 0;
+  contacts_master_elems = 0;
 }
 
 Contacts_list::Contacts_list(Problem &pbl, Teuchos::RCP<Map> map, Teuchos::RCP<Elements_list> elements_list, Teuchos::RCP<Nodes_list> nodes_list)
 {
-    typedef Tpetra::Vector<>::local_ordinal_type local_ordinal_type;
-    typedef Tpetra::Vector<>::global_ordinal_type global_ordinal_type;
-
-    contacts_number = pbl.Contacts.size();
-    local_ordinal_type maxNumSlaveNodesPerContactFace = 0;
-    local_ordinal_type maxNumMasterNodesPerContactFace = 0;
-    local_ordinal_type maxNumSlaveElemsPerContactFace = 0;
-    local_ordinal_type maxNumMasterElemsPerContactFace = 0;
-
-    std::vector<std::vector<global_ordinal_type>> contactSlaveNodes ={};
-    std::vector<std::vector<global_ordinal_type>> contactMasterNodes ={};
-    std::vector<std::vector<global_ordinal_type>> contactSlaveElems ={};
-    std::vector<std::vector<global_ordinal_type>> contactMasterElems ={};
+  typedef Tpetra::Vector<>::local_ordinal_type local_ordinal_type;
+  typedef Tpetra::Vector<>::global_ordinal_type global_ordinal_type;
+
+  contacts_number = pbl.Contacts.size();
+  local_ordinal_type maxNumSlaveNodesPerContactFace = 0;
+  local_ordinal_type maxNumMasterNodesPerContactFace = 0;
+  local_ordinal_type maxNumSlaveElemsPerContactFace = 0;
+  local_ordinal_type maxNumMasterElemsPerContactFace = 0;
+  local_ordinal_type maxNumSlaveDirichletNodesPerContactFace = 0;
+
+  for (auto i = 0; i < contacts_number; ++i)
+  {
+    std::vector<global_ordinal_type> myContactSlaveNodes;
+    std::vector<global_ordinal_type> myContactSlaveElems;
+    for (auto j = 0; j < pbl.Contacts[i]->slave_elems.size(); ++j)
+    {
+      local_ordinal_type e = map->mapElems->getLocalElement(pbl.Contacts[i]->slave_elems[j]->no-1);
+      if      (elements_list->Get_Element_Type(e) == static_cast<int>(ELTYPE::QUAD4));
+      else if (elements_list->Get_Element_Type(e) == static_cast<int>(ELTYPE::TRI3) );
+      else
+        continue;
+      myContactSlaveElems.push_back(e);
+      for ( auto n1 = 0; n1 < elements_list->Get_Element_Size(e); ++n1)
+      {
+        global_ordinal_type i_global = map->mapNodesWO->getGlobalElement(elements_list->Get_Element_Node(e,n1));
+
+        if(std::find(myContactSlaveNodes.begin(),myContactSlaveNodes.end(),i_global)==myContactSlaveNodes.end())
+          myContactSlaveNodes.push_back(i_global);
+      }
+    }
 
-    std::vector<std::vector<global_ordinal_type>> contactSlaveDirichletNodes ={};
+    int myNumContactSlaveNodes = myContactSlaveNodes.size();
+    int myNumContactSlaveElems = myContactSlaveElems.size();
 
+    int myNumContactSlaveDirichletNodes = pbl.Contacts[i]->slave_dirichlet_nodes.size();
 
-    for (auto i = 0; i < contacts_number; ++i)
+    std::vector<global_ordinal_type> myContactMasterNodes;
+    std::vector<global_ordinal_type> myContactMasterElems;
+    for (auto j = 0; j < pbl.Contacts[i]->master_elems.size(); ++j)
     {
-        std::vector<global_ordinal_type> myContactSlaveNodes ={};
-        std::vector<global_ordinal_type> myContactSlaveElems ={};
-        for (auto j = 0; j < pbl.Contacts[i]->slave_elems.size(); ++j)
-        {
-            local_ordinal_type e = map->mapElems->getLocalElement(pbl.Contacts[i]->slave_elems[j]->no-1);
-            if      (elements_list->Get_Element_Type(e) == static_cast<int>(ELTYPE::QUAD4));
-            else if (elements_list->Get_Element_Type(e) == static_cast<int>(ELTYPE::TRI3) );
-            else
-                continue;
-            myContactSlaveElems.push_back(e);
-            for ( auto n1 = 0; n1 < elements_list->Get_Element_Size(e); ++n1)
-            {
-                global_ordinal_type i_global = map->mapNodesWO->getGlobalElement(elements_list->Get_Element_Node(e,n1));
-                if(std::find(myContactSlaveNodes.begin(),myContactSlaveNodes.end(),i_global)==myContactSlaveNodes.end())
-                    myContactSlaveNodes.push_back(i_global);
-
-                //std::cout << pbl.Contacts[i]->slave_elems[j]->no << " " << e << " " <<    i_global << " " << map->mapNodes->getGlobalElement(elements_list->Get_Element_Node(e,n1)) << " " << elements_list->Get_Element_Node(e,n1)  << pbl.msh->elems[e]->nodes[n1]->row << std::endl;
-            }
-        }
-        contactSlaveNodes.push_back(myContactSlaveNodes);
-        contactSlaveElems.push_back(myContactSlaveElems);
-
-        std::vector<global_ordinal_type> myContactSlaveDirichletNodes ={};
-        for (auto j = 0; j < pbl.Contacts[i]->slave_dirichlet_nodes.size(); ++j)
-        {
-            myContactSlaveDirichletNodes.push_back(pbl.Contacts[i]->slave_dirichlet_nodes[j]);
-        }
-        //std::cout << i << ": " << myContactSlaveDirichletNodes.size() << std::endl;
-        contactSlaveDirichletNodes.push_back(myContactSlaveDirichletNodes);
+      local_ordinal_type e = map->mapElems->getLocalElement(pbl.Contacts[i]->master_elems[j]->no-1);
+      if      (elements_list->Get_Element_Type(e) == static_cast<int>(ELTYPE::QUAD4));
+      else if (elements_list->Get_Element_Type(e) == static_cast<int>(ELTYPE::TRI3));
+      else if (elements_list->Get_Element_Type(e) == static_cast<int>(ELTYPE::LINE2));
+      else
+        continue;
+      myContactMasterElems.push_back(e);
+      for ( auto n1 = 0; n1 < elements_list->Get_Element_Size(e); ++n1)
+      {
+        global_ordinal_type i_global = map->mapNodesWO->getGlobalElement(elements_list->Get_Element_Node(e,n1));
+        if(std::find(myContactMasterNodes.begin(),myContactMasterNodes.end(),i_global)==myContactMasterNodes.end())
+          myContactMasterNodes.push_back(i_global);
+      }
+    }
 
-        std::vector<global_ordinal_type> myContactMasterNodes ={};
-        std::vector<global_ordinal_type> myContactMasterElems ={};
-        for (auto j = 0; j < pbl.Contacts[i]->master_elems.size(); ++j)
+    int myNumContactMasterNodes = myContactMasterNodes.size();
+    int myNumContactMasterElems = myContactMasterElems.size();
+
+    if (maxNumSlaveNodesPerContactFace < myNumContactSlaveNodes)
+      maxNumSlaveNodesPerContactFace = myNumContactSlaveNodes;
+    if (maxNumSlaveElemsPerContactFace < myNumContactSlaveElems)
+      maxNumSlaveElemsPerContactFace = myNumContactSlaveElems;
+
+    if (maxNumSlaveDirichletNodesPerContactFace < myNumContactSlaveDirichletNodes)
+      maxNumSlaveDirichletNodesPerContactFace = myNumContactSlaveDirichletNodes;
+
+    if (maxNumMasterNodesPerContactFace < myNumContactMasterNodes)
+      maxNumMasterNodesPerContactFace = myNumContactMasterNodes;
+    if (maxNumMasterElemsPerContactFace < myNumContactMasterElems)
+      maxNumMasterElemsPerContactFace = myNumContactMasterElems;
+  }
+
+  contacts_slave_nodes_size = 0;
+  contacts_slave_elems_size = 1;
+  contacts_master_nodes_size = 2;
+  contacts_master_elems_size = 3;
+  contacts_slave_normal = 4;
+  contacts_master_normal = 7;
+  contacts_slave_nodes = 10;
+  contacts_slave_elems = contacts_slave_nodes + maxNumSlaveNodesPerContactFace;
+  contacts_master_nodes = contacts_slave_elems + maxNumSlaveElemsPerContactFace;
+  contacts_master_elems = contacts_master_nodes + maxNumMasterNodesPerContactFace;
+
+  size_t numMaxDirichletSlaveNode = 1;
+  for (auto i = 0; i < contacts_number; ++i)
+    if(numMaxDirichletSlaveNode < pbl.Contacts[i]->slave_dirichlet_nodes.size() )
+      numMaxDirichletSlaveNode = pbl.Contacts[i]->slave_dirichlet_nodes.size();
+
+  size_t numMaxClosedNode = 1;
+  for (auto i = 0; i < contacts_number; ++i)
+    if(numMaxClosedNode < pbl.Contacts[i]->initially_closed_nodes.size() )
+      numMaxClosedNode = pbl.Contacts[i]->initially_closed_nodes.size();
+
+  Kokkos::resize(contacts, contacts_number,contacts_master_elems+maxNumMasterElemsPerContactFace);
+  Kokkos::resize(contacts_normal,contacts_number,6);
+  Kokkos::resize(shape_function,contacts_number);
+  Kokkos::resize(signorini,contacts_number);
+  Kokkos::resize(initially_open,contacts_number);
+  Kokkos::resize(no_update,contacts_number);
+  Kokkos::resize(sticking,contacts_number);
+  Kokkos::resize(mesh_tying,contacts_number);
+  Kokkos::resize(contacts_scale,contacts_number);
+  Kokkos::resize(contacts_dirichlet,contacts_number,numMaxDirichletSlaveNode+1);
+  Kokkos::resize(initially_closed_nodes,contacts_number,numMaxClosedNode+1);
+
+
+  for (auto i = 0; i < contacts_number; ++i)
+  {
+    contacts_dirichlet(i,0) = pbl.Contacts[i]->slave_dirichlet_nodes.size();
+    for (auto j = 0; j < contacts_dirichlet(i,0); ++j)
+      contacts_dirichlet(i,j+1) = pbl.Contacts[i]->slave_dirichlet_nodes[j];
+  }
+
+  for (auto i = 0; i < contacts_number; ++i)
+  {
+    initially_closed_nodes(i,0) = pbl.Contacts[i]->initially_closed_nodes.size();
+    for (auto j = 0; j < initially_closed_nodes(i,0); ++j)
+      initially_closed_nodes(i,j+1) = pbl.Contacts[i]->initially_closed_nodes[j];
+  }
+
+  for (auto i = 0; i < contacts_number; ++i)
+  {
+    int myContactSlaveNodes_id = 0;
+    int myContactSlaveElems_id = 0;
+
+    std::vector<global_ordinal_type> myContactSlaveNodes;
+
+    for (auto j = 0; j < pbl.Contacts[i]->slave_elems.size(); ++j)
+    {
+      local_ordinal_type e = map->mapElems->getLocalElement(pbl.Contacts[i]->slave_elems[j]->no-1);
+      if      (elements_list->Get_Element_Type(e) == static_cast<int>(ELTYPE::QUAD4));
+      else if (elements_list->Get_Element_Type(e) == static_cast<int>(ELTYPE::TRI3) );
+      else
+        continue;
+      contacts(i,contacts_slave_elems+myContactSlaveElems_id) = e;
+      ++myContactSlaveElems_id;
+      for ( auto n1 = 0; n1 < elements_list->Get_Element_Size(e); ++n1)
+      {
+        global_ordinal_type i_global = map->mapNodesWO->getGlobalElement(elements_list->Get_Element_Node(e,n1));
+        if(std::find(myContactSlaveNodes.begin(),myContactSlaveNodes.end(),i_global)==myContactSlaveNodes.end())
         {
-            local_ordinal_type e = map->mapElems->getLocalElement(pbl.Contacts[i]->master_elems[j]->no-1);
-            if      (elements_list->Get_Element_Type(e) == static_cast<int>(ELTYPE::QUAD4));
-            else if (elements_list->Get_Element_Type(e) == static_cast<int>(ELTYPE::TRI3));
-            else if (elements_list->Get_Element_Type(e) == static_cast<int>(ELTYPE::LINE2));
-            else
-                continue;
-            myContactMasterElems.push_back(e);
-            for ( auto n1 = 0; n1 < elements_list->Get_Element_Size(e); ++n1)
-            {
-                global_ordinal_type i_global = map->mapNodesWO->getGlobalElement(elements_list->Get_Element_Node(e,n1));
-                if(std::find(myContactMasterNodes.begin(),myContactMasterNodes.end(),i_global)==myContactMasterNodes.end())
-                    myContactMasterNodes.push_back(i_global);
-            }
+          myContactSlaveNodes.push_back(i_global);
+          ++myContactSlaveNodes_id;
         }
-
-        contactMasterNodes.push_back(myContactMasterNodes);
-        contactMasterElems.push_back(myContactMasterElems);
-
-        if (maxNumSlaveNodesPerContactFace < myContactSlaveNodes.size())
-            maxNumSlaveNodesPerContactFace = myContactSlaveNodes.size();
-        if (maxNumSlaveElemsPerContactFace < myContactSlaveElems.size())
-            maxNumSlaveElemsPerContactFace = myContactSlaveElems.size();
-
-        if (maxNumMasterNodesPerContactFace < myContactMasterNodes.size())
-            maxNumMasterNodesPerContactFace = myContactMasterNodes.size();
-        if (maxNumMasterElemsPerContactFace < myContactMasterElems.size())
-            maxNumMasterElemsPerContactFace = myContactMasterElems.size();
+      }
     }
+    for (auto j = 0; j < myContactSlaveNodes.size(); ++j)
+      contacts(i,contacts_slave_nodes+j) = myContactSlaveNodes[j];
 
-    contacts_slave_nodes_size = 0;
-    contacts_slave_elems_size = 1;
-    contacts_master_nodes_size = 2;
-    contacts_master_elems_size = 3;
-    contacts_slave_normal = 4;
-    contacts_master_normal = 7;
-    contacts_slave_nodes = 10;
-    contacts_slave_elems = contacts_slave_nodes + maxNumSlaveNodesPerContactFace;
-    contacts_master_nodes = contacts_slave_elems + maxNumSlaveElemsPerContactFace;
-    contacts_master_elems = contacts_master_nodes + maxNumMasterNodesPerContactFace;
+    contacts(i,contacts_slave_nodes_size) = myContactSlaveNodes_id;
+    contacts(i,contacts_slave_elems_size) = myContactSlaveElems_id;
 
-    //contacts_dg = contacts_master_elems+maxNumMasterElemsPerContactFace;
+    int myContactMasterNodes_id = 0;
+    int myContactMasterElems_id = 0;
 
-    contacts = Kokkos::View<int**, Kokkos::LayoutRight> ("R",contacts_number,contacts_master_elems+maxNumMasterElemsPerContactFace);
+    std::vector<global_ordinal_type> myContactMasterNodes;
 
-    contacts_normal = Kokkos::View<double**, Kokkos::LayoutRight> ("R",contacts_number,3);
+    for (auto j = 0; j < pbl.Contacts[i]->master_elems.size(); ++j)
+    {
+      local_ordinal_type e = map->mapElems->getLocalElement(pbl.Contacts[i]->master_elems[j]->no-1);
+      if      (elements_list->Get_Element_Type(e) == static_cast<int>(ELTYPE::QUAD4));
+      else if (elements_list->Get_Element_Type(e) == static_cast<int>(ELTYPE::TRI3));
+      else if (elements_list->Get_Element_Type(e) == static_cast<int>(ELTYPE::LINE2));
+      else
+        continue;
+      contacts(i,contacts_master_elems+myContactMasterElems_id) = e;
+      ++myContactMasterElems_id;
+      for ( auto n1 = 0; n1 < elements_list->Get_Element_Size(e); ++n1)
+      {
+        global_ordinal_type i_global = map->mapNodesWO->getGlobalElement(elements_list->Get_Element_Node(e,n1));
+        if(std::find(myContactMasterNodes.begin(),myContactMasterNodes.end(),i_global)==myContactMasterNodes.end())
+        {
+          myContactMasterNodes.push_back(i_global);
+          ++myContactMasterNodes_id;
+        }
+      }
+    }
 
-    shape_function = Kokkos::View<int*, Kokkos::LayoutRight> ("R",contacts_number);
+    for (auto j = 0; j < myContactMasterNodes.size(); ++j)
+      contacts(i,contacts_master_nodes+j) = myContactMasterNodes[j];
 
-    signorini = Kokkos::View<int*, Kokkos::LayoutRight> ("R",contacts_number);
-    initially_open = Kokkos::View<int*, Kokkos::LayoutRight> ("R",contacts_number);
+    contacts(i,contacts_master_nodes_size) = myContactMasterNodes_id;
+    contacts(i,contacts_master_elems_size) = myContactMasterElems_id;
+  }
 
-    no_update = Kokkos::View<int*, Kokkos::LayoutRight> ("R",contacts_number);
-    sticking = Kokkos::View<int*, Kokkos::LayoutRight> ("R",contacts_number);
-    mesh_tying = Kokkos::View<int*, Kokkos::LayoutRight> ("R",contacts_number);
-    contacts_scale = Kokkos::View<double*, Kokkos::LayoutRight> ("R",contacts_number);
+  for (auto i = 0; i < contacts_number; ++i)
+  {
 
-    for (auto i = 0; i < contacts_number; ++i)
+    for ( auto j = 0; j < 3 ; ++j)
     {
+      contacts(i,contacts_slave_normal+j) = pbl.Contacts[i]->slave_normal.x[j];
+      contacts(i,contacts_master_normal+j) = pbl.Contacts[i]->master_normal.x[j];
 
-        contacts(i,contacts_slave_nodes_size) = contactSlaveNodes[i].size();
-        contacts(i,contacts_slave_elems_size) = contactSlaveElems[i].size();
-        contacts(i,contacts_master_nodes_size) = contactMasterNodes[i].size();
-        contacts(i,contacts_master_elems_size) = contactMasterElems[i].size();
+      contacts_normal(i,j) = pbl.Contacts[i]->slave_normal.x[j];
+      contacts_normal(i,j+3) = pbl.Contacts[i]->master_normal.x[j];
+    }
 
-        for ( auto j = 0; j < 3 ; ++j)
-        {
-            contacts(i,contacts_slave_normal+j) = pbl.Contacts[i]->slave_normal.x[j];
-            contacts(i,contacts_master_normal+j) = pbl.Contacts[i]->master_normal.x[j];
+    if(pbl.Contacts[i]->standard_shape_functions)
+      shape_function (i) = 1;
+    else
+      shape_function (i) = 0;
 
-            contacts_normal(i,j) = pbl.Contacts[i]->slave_normal.x[j];
-            contacts_normal(i,j+3) = pbl.Contacts[i]->master_normal.x[j];
-        }
-
-        //contacts(i,contacts_dg) = pbl.Contacts[i]->dg;
-
-        if(pbl.Contacts[i]->standard_shape_functions)
-            shape_function (i) = 1;
-        else
-            shape_function (i) = 0;
-
-        if(pbl.Contacts[i]->is_Signorini)
-            signorini (i) = 1;
-        else
-            signorini (i) = 0;
-
-        if(pbl.Contacts[i]->mesh_tying)
-            mesh_tying (i) = 1;
-        else
-            mesh_tying (i) = 0;
-            
-        if(pbl.Contacts[i]->no_Update)
-            no_update (i) = 1;
-        else
-            no_update (i) = 0;
-
-        if(pbl.Contacts[i]->sticking)
-        {
-            sticking(i) = true;
-            at_least_one_sticking = true;
-        }
-        else
-            sticking(i) = false;
-
-        if(pbl.Contacts[i]->initial_closed)
-            initially_open(i) = 0;
-        else
-            initially_open(i) = 1;
-      
-        contacts_scale(i) = pbl.Contacts[i]->scale;
-
-        for ( auto j = 0; j < contactSlaveNodes[i].size() ; ++j)
-            contacts(i,contacts_slave_nodes+j) = contactSlaveNodes[i][j];
-        for ( auto j = 0; j < contactSlaveElems[i].size() ; ++j)
-            contacts(i,contacts_slave_elems+j) = contactSlaveElems[i][j];
-        for ( auto j = 0; j < contactMasterNodes[i].size() ; ++j)
-            contacts(i,contacts_master_nodes+j) = contactMasterNodes[i][j];
-        for ( auto j = 0; j < contactMasterElems[i].size() ; ++j)
-            contacts(i,contacts_master_elems+j) = contactMasterElems[i][j];
-
-        //
-
-        size_t numMaxDirichletSlaveNode = 1;
-        for (auto i = 0; i < contacts_number; ++i)
-            if(numMaxDirichletSlaveNode < pbl.Contacts[i]->slave_dirichlet_nodes.size() )
-                numMaxDirichletSlaveNode = pbl.Contacts[i]->slave_dirichlet_nodes.size();
-
-        contacts_dirichlet = Kokkos::View<int**, Kokkos::LayoutRight> ("R",contacts_number,numMaxDirichletSlaveNode+1);
-
-        for (auto i = 0; i < contacts_number; ++i)
-        {
-            contacts_dirichlet(i,0) = pbl.Contacts[i]->slave_dirichlet_nodes.size();
-            for (auto j = 0; j < contacts_dirichlet(i,0); ++j)
-                contacts_dirichlet(i,j+1) = pbl.Contacts[i]->slave_dirichlet_nodes[j];
-        }
+    if(pbl.Contacts[i]->is_Signorini)
+      signorini (i) = 1;
+    else
+      signorini (i) = 0;
 
+    if(pbl.Contacts[i]->mesh_tying)
+      mesh_tying (i) = 1;
+    else
+      mesh_tying (i) = 0;
 
-        //
+    if(pbl.Contacts[i]->no_Update)
+      no_update (i) = 1;
+    else
+      no_update (i) = 0;
 
-        size_t numMaxClosedNode = 1;
-        for (auto i = 0; i < contacts_number; ++i)
-            if(numMaxClosedNode < pbl.Contacts[i]->initially_closed_nodes.size() )
-                numMaxClosedNode = pbl.Contacts[i]->initially_closed_nodes.size();
+    if(pbl.Contacts[i]->sticking)
+    {
+      sticking(i) = true;
+      at_least_one_sticking = true;
+    }
+    else
+      sticking(i) = false;
 
-        initially_closed_nodes = Kokkos::View<int**, Kokkos::LayoutRight> ("R",contacts_number,numMaxClosedNode+1);
+    if(pbl.Contacts[i]->initial_closed)
+      initially_open(i) = 0;
+    else
+      initially_open(i) = 1;
 
-        for (auto i = 0; i < contacts_number; ++i)
-        {
-            initially_closed_nodes(i,0) = pbl.Contacts[i]->initially_closed_nodes.size();
-            for (auto j = 0; j < initially_closed_nodes(i,0); ++j)
-                initially_closed_nodes(i,j+1) = pbl.Contacts[i]->initially_closed_nodes[j];
-        }
-    }
+    contacts_scale(i) = pbl.Contacts[i]->scale;
+  }
 }