Update waves

Created by: acrovato

Context

This pull request brings the last developments from the flow module as well as some upgrade to waves.

Major changes in flow

• flow is now able to compute 3D flows
• The Jacobian (tangent) matrix for the Netwton method is now fully analytical
• The adjacent elements identification is now based on ```std::unordered_set````

Major changes in waves

• tbox and fwk are now included in a namespace
• waves is now fully included in a namespace
• All elements now have a private memory (Mem) which can be publicly accessed
• All elements now have a cache including Gauss points (Gauss) which can be publicly accessed
• MshDeformation now uses a linear elasticity solver
• MshCrack class implemented (creates duplicate internal boundaries in meshes)

Test

The ctest battery has been run on gaston, 100% tests passed under 110s.

Future work

• Use automated indentation rules for the whole code
• Move methods using VTK outside of MshData
• Clean tbox (mainly MshData, Tri3, Tera4 and Hex8)

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;
+  }
 }