From 6c8c78083a45955a945954ac304d435348e1f4aa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Abdoulaye=20Samak=C3=A9?= Date: Wed, 11 Sep 2024 16:54:31 +0000 Subject: [PATCH 1/2] add support for updating ghost elements 1- Add two new functions: * void updateGhostElements(std::vector& mesh_elt_ctr); * void updateGhostElements(ModelVariable& mesh_elt_ctr); 2- some cleaning up --- core/include/gmshmesh.hpp | 69 +++--- core/src/gmshmesh.cpp | 292 +++++++++-------------- model/finiteelement.cpp | 484 +++++++++++++++++++++++++------------- model/finiteelement.hpp | 40 ++-- 4 files changed, 489 insertions(+), 396 deletions(-) diff --git a/core/include/gmshmesh.hpp b/core/include/gmshmesh.hpp index 87003e1fc..7fca96646 100644 --- a/core/include/gmshmesh.hpp +++ b/core/include/gmshmesh.hpp @@ -23,8 +23,6 @@ #include #include -//#include -//#include #include "debug.hpp" extern "C" @@ -50,14 +48,8 @@ class GmshMesh GmshMesh( Communicator const& comm = Environment::comm() ); - GmshMesh( GmshMesh const& mesh ); - // GmshMesh(std::vector const& nodes, - // std::vector const& edges, - // std::vector const& triangles, - // Communicator const& comm = Environment::comm()); - void readFromFile(std::string const& filename, std::string const& format="ascii"); void readFromFileBinary(std::ifstream& ifs); void readFromFileASCII(std::ifstream& ifs); @@ -79,16 +71,18 @@ class GmshMesh int numGlobalNodes() const {return M_global_num_nodes;} int numGlobalNodesFromSarialMesh() const {return M_global_num_nodes_from_serial;} int numNodes() const {return M_num_nodes;} - int numTriangles() const {return M_num_triangles;} - int numEdges() const {return M_num_edges;} - - int numLocalNodesWithoutGhost() const {return M_nldof_without_ghost;} int numLocalNodesWithGhost() const {return M_nldof_with_ghost;} + int numLocalNodesWithoutGhost() const {return M_nldof_without_ghost;} int numLocalGhost() const {return M_nlghost;} int numGlobalElements() const {return M_global_num_elements;} int numGlobalElementsFromSarialMesh() const {return M_global_num_elements_from_serial;} - int numTrianglesWithoutGhost() const {return M_num_triangles_without_ghost;} + int numTriangles() const {return M_num_triangles;} + int numTrianglesWithGhost() const {return M_nlelt_with_ghost;} + int numTrianglesWithoutGhost() const {return M_nlelt_without_ghost;} + int numLocalGhostElt() const {return M_nlghost_elt;} + + int numEdges() const {return M_num_edges;} void setCommunicator(Communicator const& comm) {M_comm=comm;} void setOrdering(std::string const& order) {M_ordering=order;} @@ -103,23 +97,21 @@ class GmshMesh void stereographicProjection(); - std::vector const& localDofWithoutGhost() const {return M_local_dof_without_ghost;} std::vector const& localDofWithGhost() const {return M_local_dof_with_ghost;} + std::vector const& localDofWithGhostDefault() const {return M_local_dof_with_ghost_default;} + std::vector const& localDofWithoutGhost() const {return M_local_dof_without_ghost;} std::vector const& localGhost() const {return M_local_ghost;} - std::vector const& localDofWithGhostInit() const {return M_local_dof_with_ghost_init;} - std::vector const& trianglesIdWithGhost() const {return M_triangles_id_with_ghost;} + std::vector const& localEltIdWithGhost() const {return M_local_eltid_with_ghost;} + std::vector const& localEltIdWithGhostDefault() const {return M_local_eltid_with_ghost_default;} + std::vector const& localEltIdWithoutGhost() const {return M_local_eltid_without_ghost;} + std::vector const& localGhostElt() const {return M_local_ghost_elt;} - bimap_type const& transferMapInit() const {return M_transfer_map;} - bimap_type const& transferMap() const {return M_transfer_map_reordered;} + bimap_type const& transferMapDof() const {return M_transfer_map_dof;} + bimap_type const& transferMapDofDefault() const {return M_transfer_map_dof_default;} bimap_type const& transferMapElt() const {return M_transfer_map_elt;} - //bimap_type const& transferMapReordered() const {return M_transfer_map_reordered;} - - //bimap_type const& mapNodes() const {return M_reorder_map_nodes;} - //std::map const& mapNodes() const {return M_reorder_map_nodes;} - //bimap_type const& mapElements() const {return M_reorder_map_elements;} - //std::map const& mapElements() const {return M_reorder_map_elements;} + bimap_type const& transferMapEltDefault() const {return M_transfer_map_elt_default;} std::vector const& mapNodes() const {return M_map_nodes;} std::vector const& mapElements() const {return M_map_elements;} @@ -172,42 +164,43 @@ class GmshMesh std::vector M_edges; std::vector M_local_dof_with_ghost; - std::vector M_local_dof_with_ghost_init; + std::vector M_local_dof_with_ghost_default; std::vector M_local_dof_without_ghost; std::vector M_local_ghost; - std::vector M_triangles_id_with_ghost; + + std::vector M_local_eltid_with_ghost; + std::vector M_local_eltid_with_ghost_default; + std::vector M_local_eltid_without_ghost; + std::vector M_local_ghost_elt; int M_global_num_nodes; int M_global_num_nodes_from_serial; int M_num_nodes; - int M_num_triangles; - int M_num_edges; - int M_nldof_with_ghost; int M_nldof_without_ghost; int M_nlghost; int M_global_num_elements; int M_global_num_elements_from_serial; - int M_num_triangles_without_ghost; + int M_num_triangles; + int M_nlelt_with_ghost; + int M_nlelt_without_ghost; + int M_nlghost_elt; + int M_num_edges; + bimap_type M_transfer_map_dof; + bimap_type M_transfer_map_dof_default; - bimap_type M_transfer_map; - bimap_type M_transfer_map_reordered; bimap_type M_transfer_map_elt; + bimap_type M_transfer_map_elt_default; // container for storing the mesh marker names std::map > M_marker_names; - //bimap_type M_reorder_map_nodes; - //std::map M_reorder_map_nodes; (not stored in memory) - //bimap_type M_reorder_map_elements; - //std::map M_reorder_map_elements; (not stored in memory) - std::vector M_map_elements; std::vector M_map_nodes; - std::map > timer; + std::map > M_timer; }; } // Nextsim diff --git a/core/src/gmshmesh.cpp b/core/src/gmshmesh.cpp index 359298545..78685b8c3 100644 --- a/core/src/gmshmesh.cpp +++ b/core/src/gmshmesh.cpp @@ -25,28 +25,33 @@ GmshMesh::GmshMesh(Communicator const& comm) M_version("2.2"), M_ordering("gmsh"), M_nodes(), - M_triangles(), - M_edges(), M_nodes_vec(), M_local_dof_with_ghost(), + M_local_dof_with_ghost_default(), M_local_dof_without_ghost(), M_local_ghost(), - M_triangles_id_with_ghost(), + M_triangles(), + M_local_eltid_with_ghost(), + M_local_eltid_with_ghost_default(), + M_local_eltid_without_ghost(), + M_local_ghost_elt(), + M_edges(), M_num_nodes(0), + M_nldof_with_ghost(0), + M_nldof_without_ghost(0), + M_nlghost(0), M_num_triangles(0), + M_nlelt_with_ghost(0), + M_nlelt_without_ghost(0), + M_nlghost_elt(0), M_num_edges(0), - M_nldof_with_ghost(), - M_nldof_without_ghost(), - M_nlghost(), - M_num_triangles_without_ghost(), - M_transfer_map(), - M_transfer_map_reordered(), + M_transfer_map_dof(), + M_transfer_map_dof_default(), M_transfer_map_elt(), - //M_reorder_map_nodes(), - //M_reorder_map_elements(), + M_transfer_map_elt_default(), M_map_nodes(), M_map_elements(), - timer(), + M_timer(), M_mppfile(Environment::nextsimMppfile()), M_log_level(Environment::logLevel()), M_log_all(Environment::logAll()) @@ -57,30 +62,37 @@ GmshMesh::GmshMesh(GmshMesh const& mesh) M_comm(mesh.M_comm), M_version(mesh.M_version), M_ordering(mesh.M_ordering), - M_mppfile(mesh.M_mppfile), - M_log_level(mesh.M_log_level), - M_log_all(mesh.M_log_all), M_nodes(mesh.M_nodes), - M_triangles(mesh.M_triangles), - M_edges(mesh.M_edges), M_nodes_vec(mesh.M_nodes_vec), M_local_dof_with_ghost(mesh.M_local_dof_with_ghost), + M_local_dof_with_ghost_default(mesh.M_local_dof_with_ghost_default), M_local_dof_without_ghost(mesh.M_local_dof_without_ghost), M_local_ghost(mesh.M_local_ghost), - M_triangles_id_with_ghost(mesh.M_triangles_id_with_ghost), + M_triangles(mesh.M_triangles), + M_local_eltid_with_ghost(mesh.M_local_eltid_with_ghost), + M_local_eltid_with_ghost_default(mesh.M_local_eltid_with_ghost_default), + M_local_eltid_without_ghost(mesh.M_local_eltid_without_ghost), + M_local_ghost_elt(mesh.M_local_ghost_elt), + M_edges(mesh.M_edges), M_num_nodes(mesh.M_num_nodes), - M_num_triangles(mesh.M_num_triangles), - M_num_edges(mesh.M_num_edges), M_nldof_with_ghost(mesh.M_nldof_with_ghost), M_nldof_without_ghost(mesh.M_nldof_without_ghost), M_nlghost(mesh.M_nlghost), - M_num_triangles_without_ghost(mesh.M_num_triangles_without_ghost), - M_transfer_map(mesh.M_transfer_map), - M_transfer_map_reordered(mesh.M_transfer_map_reordered), - //M_reorder_map_nodes(mesh.M_reorder_map_nodes), - //M_reorder_map_elements(mesh.M_reorder_map_elements) + M_num_triangles(mesh.M_num_triangles), + M_nlelt_with_ghost(mesh.M_nlelt_with_ghost), + M_nlelt_without_ghost(mesh.M_nlelt_without_ghost), + M_nlghost_elt(mesh.M_nlghost_elt), + M_num_edges(mesh.M_num_edges), + M_transfer_map_dof(mesh.M_transfer_map_dof), + M_transfer_map_dof_default(mesh.M_transfer_map_dof_default), + M_transfer_map_elt(mesh.M_transfer_map_elt), + M_transfer_map_elt_default(mesh.M_transfer_map_elt_default), M_map_nodes(mesh.M_map_nodes), - M_map_elements(mesh.M_map_elements) + M_map_elements(mesh.M_map_elements), + M_timer(mesh.M_timer), + M_mppfile(mesh.M_mppfile), + M_log_level(mesh.M_log_level), + M_log_all(mesh.M_log_all) {} #if 0 @@ -129,11 +141,11 @@ GmshMesh::readFromFile(std::string const& gmshmshfile, std::string const& format } // create nodal partitions - timer["in.nodal"].first.restart(); + M_timer["in.nodal"].first.restart(); if (M_comm.size() > 1) this->nodalGrid(); - LOG(DEBUG)<<"-------------------INSIDE: NODALGRID done in "<< timer["in.nodal"].first.elapsed() <<"s\n"; + LOG(DEBUG)<<"-------------------INSIDE: NODALGRID done in "<< M_timer["in.nodal"].first.elapsed() <<"s\n"; } void @@ -970,9 +982,6 @@ GmshMesh::nodalGrid() std::sort(all_local_nodes.begin(), all_local_nodes.end()); all_local_nodes.erase(std::unique( all_local_nodes.begin(), all_local_nodes.end() ), all_local_nodes.end()); - //std::sort(all_local_nodes.begin(), all_local_nodes.end()); // already sorted - - //std::sort(M_local_dof_without_ghost.begin(), M_local_dof_without_ghost.end()); // already sorted std::set_difference(all_local_nodes.begin(), all_local_nodes.end(), M_local_dof_without_ghost.begin(), M_local_dof_without_ghost.end(), @@ -1092,74 +1101,6 @@ GmshMesh::nodalGrid() M_local_dof_without_ghost = renumbering[M_comm.rank()]; } -#if 0 - else - { - std::sort(M_local_ghost.begin(), M_local_ghost.end()); - std::vector> global_ghosts; - // gather operations for global ghost numbering - boost::mpi::all_gather(M_comm, - M_local_ghost, - global_ghosts); - - - std::vector global_from_mesh(M_num_nodes); - std::iota(global_from_mesh.begin(), global_from_mesh.end(), 1); - - std::vector global_from_code;//(num_nodes); - - for (int ii=0; ii diff_trs; - std::set_difference(global_from_mesh.begin(), global_from_mesh.end(), - global_from_code.begin(), global_from_code.end(), - std::back_inserter(diff_trs)); - -#if 0 - for (int kk=0; kk dom_ids; - - for (int ii=0; ii reorder; - std::map M_reorder_map_nodes; + std::map reorder_map_global_dof_uv; + std::map reorder_map_global_dof; int cpts = 0; int cpts_dom = 0; @@ -1179,29 +1119,17 @@ GmshMesh::nodalGrid() { int sr = renumbering[ii].size(); - //for (int jj=0; jjsecond; - int rdofv = reorder.find(local_dof_with_ghost[k]+M_num_nodes)->second; + int rdof = reorder_map_global_dof_uv.find(M_local_dof_with_ghost_default[k])->second; + int rdofv = reorder_map_global_dof_uv.find(M_local_dof_with_ghost_default[k]+M_num_nodes)->second; // mapping from new global numbering to local numbering - M_transfer_map_reordered.insert(position(rdof,k+1)); + M_transfer_map_dof.insert(position(rdof,k+1)); M_local_dof_with_ghost[k] = rdof; M_local_dof_with_ghost[k+M_nldof_with_ghost] = rdofv; - M_nodes[k+1] = M_nodes_vec[local_dof_with_ghost[k]-1]; + M_nodes[k+1] = M_nodes_vec[M_local_dof_with_ghost_default[k]-1]; if (k < M_nldof_without_ghost) @@ -1264,7 +1190,7 @@ GmshMesh::nodalGrid() M_global_num_nodes = M_num_nodes; M_num_nodes = M_nodes.size(); - std::vector triangles_num_without_ghost; + //std::vector M_local_eltid_without_ghost; std::vector _triangles = M_triangles; M_triangles.resize(0); @@ -1276,7 +1202,7 @@ GmshMesh::nodalGrid() for (int i=0; i<3; ++i) { - if ((M_transfer_map.left.find(it->indices[i]) == M_transfer_map.left.end())) + if ((M_transfer_map_dof_default.left.find(it->indices[i]) == M_transfer_map_dof_default.left.end())) { _test = true; break; @@ -1291,9 +1217,9 @@ GmshMesh::nodalGrid() // new add for (int i=0; i<3; ++i) { - int rdof = reorder.find(it->indices[i])->second; + int rdof = reorder_map_global_dof_uv.find(it->indices[i])->second; - it->indices[i] = M_transfer_map.left.find(it->indices[i])->second; + it->indices[i] = M_transfer_map_dof_default.left.find(it->indices[i])->second; if (std::binary_search(M_local_ghost.begin(),M_local_ghost.end(),rdof)) { @@ -1306,7 +1232,7 @@ GmshMesh::nodalGrid() if ((M_comm.rank() == it->partition)) { - triangles_num_without_ghost.push_back(it->number); + M_local_eltid_without_ghost.push_back(it->number); } } } @@ -1314,8 +1240,8 @@ GmshMesh::nodalGrid() M_num_triangles = M_triangles.size(); // check the nodal partitions - int elt_size = triangles_num_without_ghost.size(); - int num_elements = boost::mpi::all_reduce(M_comm, elt_size, std::plus()); + int elt_lsize = M_local_eltid_without_ghost.size(); + int num_elements = boost::mpi::all_reduce(M_comm, elt_lsize, std::plus()); if (M_global_num_elements_from_serial != num_elements) { @@ -1327,12 +1253,7 @@ GmshMesh::nodalGrid() // --------------------------------BEGINNING------------------------- renumbering.resize(0); int num_trls; - this->allGather(triangles_num_without_ghost, renumbering, num_trls); - - // // gather operations for global element renumbering - // boost::mpi::all_gather(M_comm, - // triangles_num_without_ghost, - // renumbering); + this->allGather(M_local_eltid_without_ghost, renumbering, num_trls); std::vector diff_trs; @@ -1347,7 +1268,6 @@ GmshMesh::nodalGrid() { for (int jj=0; jj global_trs(all_trs.size());//(M_reorder_map_elements.size()); + std::vector global_trs(all_trs.size()); std::iota(global_trs.begin(), global_trs.end(), 1); std::set_difference(global_trs.begin(), global_trs.end(), @@ -1379,7 +1299,7 @@ GmshMesh::nodalGrid() // elements in partition first and ghost at the end _triangles = M_triangles; M_triangles.resize(0); - M_num_triangles_without_ghost = 0; + //M_nlelt_without_ghost = 0; for (auto it=_triangles.begin(), end=_triangles.end(); it!=end; ++it) { @@ -1393,7 +1313,7 @@ GmshMesh::nodalGrid() if ((M_comm.rank() == min_rank)) { - triangles_num_without_ghost.push_back(it->number); + M_local_eltid_without_ghost.push_back(it->number); it->partition = M_comm.rank(); } } @@ -1402,8 +1322,8 @@ GmshMesh::nodalGrid() if (M_comm.rank() == it->partition) { M_triangles.push_back(*it); - M_triangles_id_with_ghost.push_back(it->number); - ++M_num_triangles_without_ghost; + M_local_eltid_with_ghost.push_back(it->number); + //++M_nlelt_without_ghost; } } @@ -1412,16 +1332,10 @@ GmshMesh::nodalGrid() if (M_comm.rank() != it->partition) { M_triangles.push_back(*it); - M_triangles_id_with_ghost.push_back(it->number); + M_local_eltid_with_ghost.push_back(it->number); } } - for (int k=0; kindices[i])->second; + int rdof = reorder_map_global_dof_uv.find(it->indices[i])->second; - it->indices[i] = M_transfer_map.left.find(it->indices[i])->second; + it->indices[i] = M_transfer_map_dof_default.left.find(it->indices[i])->second; if (std::binary_search(M_local_ghost.begin(),M_local_ghost.end(),rdof)) { @@ -1443,40 +1357,65 @@ GmshMesh::nodalGrid() } // --------------------------------BEGINNING------------------------- - std::sort(triangles_num_without_ghost.begin(), triangles_num_without_ghost.end()); + std::sort(M_local_eltid_without_ghost.begin(), M_local_eltid_without_ghost.end()); renumbering.resize(0); - //int num_trls; - allGather(triangles_num_without_ghost, renumbering, num_trls); - - // // gather operations for global element renumbering - // boost::mpi::all_gather(M_comm, - // triangles_num_without_ghost, - // renumbering); + allGather(M_local_eltid_without_ghost, renumbering, num_trls); + // consider cpts = 0; cpts_dom = 0; - std::map M_reorder_map_elements; + std::map reorder_map_global_elt; for (int ii=0; iisecond; + // mapping from new global numbering to local numbering + M_transfer_map_elt.insert(position(reltid,k+1)); + + M_local_eltid_with_ghost[k] = reltid; + + if (k < M_nlelt_without_ghost) + { + M_local_eltid_without_ghost[k] = reltid; + } + else + { + //M_local_ghost_elt[k-M_nlelt_without_ghost] = reltid; + M_local_ghost_elt.push_back(reltid); + } + } + // --------------------------------END------------------------------- - M_global_num_elements = M_reorder_map_elements.size(); + std::sort(M_local_ghost_elt.begin(), M_local_ghost_elt.end()); + M_nlghost_elt = M_local_ghost_elt.size(); + M_global_num_elements = reorder_map_global_elt.size(); // std::cout<<"["<< this->comm().rank() << "] M_global_num_elements= "<< M_global_num_elements <<"\n"; - // std::cout<<"["<< this->comm().rank() << "] M_num_elements = "<< triangles_num_without_ghost.size() <<"\n"; + // std::cout<<"["<< this->comm().rank() << "] M_num_elements = "<< M_local_eltid_without_ghost.size() <<"\n"; // elements: boost::bimap -> std::map (optimization) M_map_elements.resize(M_global_num_elements); cpts = 0; - for (auto it=M_reorder_map_elements.begin(), end=M_reorder_map_elements.end(); it!=end; ++it) + for (auto it=reorder_map_global_elt.begin(), end=reorder_map_global_elt.end(); it!=end; ++it) { // if ((M_comm.rank()==0) && (cpts < 100)) // std::cout<<"--Key= "<< it->first <<" Value= "<< it->second <<"\n"; @@ -1485,16 +1424,13 @@ GmshMesh::nodalGrid() } // nodes: boost::bimap -> std::map (optimization) - M_map_nodes.resize(M_reorder_map_nodes.size()); + M_map_nodes.resize(reorder_map_global_dof.size()); cpts = 0; - for (auto it=M_reorder_map_nodes.begin(), end=M_reorder_map_nodes.end(); it!=end; ++it) + for (auto it=reorder_map_global_dof.begin(), end=reorder_map_global_dof.end(); it!=end; ++it) { - // if ((M_comm.rank()==0) && (cpts < 100)) - // std::cout<<"--Key= "<< it->first <<" Value= "<< it->second <<"\n"; M_map_nodes[cpts] = it->second-1; ++cpts; } - } void diff --git a/model/finiteelement.cpp b/model/finiteelement.cpp index d97ec956c..5edd8bf56 100644 --- a/model/finiteelement.cpp +++ b/model/finiteelement.cpp @@ -85,13 +85,12 @@ FiniteElement::distributedMeshProcessing(bool start) M_nodes = M_mesh.nodes(); M_num_elements = M_mesh.numTriangles(); - M_ndof = M_mesh.numGlobalNodes(); + M_local_nelements = M_mesh.numTrianglesWithoutGhost(); + M_num_nodes = M_mesh.numLocalNodesWithGhost(); M_local_ndof = M_mesh.numLocalNodesWithoutGhost(); - M_local_ndof_ghost = M_mesh.numLocalNodesWithGhost(); - - M_local_nelements = M_mesh.numTrianglesWithoutGhost(); - M_num_nodes = M_local_ndof_ghost; + M_ndof = M_mesh.numGlobalNodes(); + //M_num_nodes = M_local_ndof_ghost; chrono.restart(); this->bcMarkedNodes(); @@ -110,8 +109,71 @@ FiniteElement::distributedMeshProcessing(bool start) LOG(DEBUG)<<"-------------------CONNECTIVITY done in "<< chrono.elapsed() <<"s\n"; chrono.restart(); - this->initUpdateGhosts(); - LOG(DEBUG)<<"-------------------INITUPDATEGHOSTS done in "<< chrono.elapsed() <<"s\n"; + this->initUpdateGhostNodes(); + LOG(DEBUG)<<"-------------------INITUPDATEGHOSTNODES done in "<< chrono.elapsed() <<"s\n"; + + chrono.restart(); + this->initUpdateGhostElements(); + LOG(DEBUG)<<"-------------------INITUPDATEGHOSTELEMENTS done in "<< chrono.elapsed() <<"s\n"; + +#if 1 + //std::vector contnr(M_num_elements,0.); + ModelVariable contnr;//(M_num_elements,0.); + contnr.resize(M_num_elements); + boost::minstd_rand intgen; + boost::uniform_01 gen(intgen); + int ghstsize = M_mesh.localGhostElt().size(); + + for (int i=0; i " << contnr[i + M_local_nelements] + <<"\n"; + } + //fileout.close(); + + this->updateGhostElements(contnr); + + M_comm.barrier(); + if (ghstsize != 0) + fileout<<"----------------------------F2----------------------------\n"; + if (ghstsize != 0) + fileout<<" ["<< M_rank <<"]" <<"\n"; + + for (int i=0; i " << contnr[i + M_local_nelements] + <<"\n"; + } + fileout.close(); +#endif #if 0 // LOG(DEBUG) << NODES = "<< M_mesh.numGlobalNodes() << " --- "<< M_local_ndof <<"\n"; @@ -149,21 +211,6 @@ FiniteElement::distributedMeshProcessing(bool start) void FiniteElement::bcMarkedNodes() { -#if 0 - M_dirichlet_flags_root.resize(5); - M_neumann_flags_root.resize(11); - - M_dirichlet_flags_root = {2,3,11,12,13}; - M_neumann_flags_root = {4,5,6,7,1,8,9,10,14,15,16}; - - // M_dirichlet_flags_root.resize(16); - // std::iota(M_dirichlet_flags_root.begin(), M_dirichlet_flags_root.end(), 1); - // M_neumann_flags_root.resize(0); - -#endif - - //std::cout<<"["<< M_rank << "] NDOFS= "<< M_num_nodes << " --- "<< M_local_ndof <<"\n"; - std::vector flags_size_root(2); if (M_rank == 0) { @@ -180,23 +227,6 @@ FiniteElement::bcMarkedNodes() { std::copy(M_dirichlet_flags_root.begin(), M_dirichlet_flags_root.end(), std::back_inserter(flags_root)); std::copy(M_neumann_flags_root.begin(), M_neumann_flags_root.end(), std::back_inserter(flags_root)); - -#if 0 - for (int i=0; i sizes_elements = M_sizes_elements_with_ghost; @@ -2051,11 +2080,11 @@ FiniteElement::gatherSizes() int out_size_elements = std::accumulate(sizes_elements.begin(),sizes_elements.end(),0); M_id_elements.resize(out_size_elements); - boost::mpi::gatherv(M_comm, M_mesh.trianglesIdWithGhost(), &M_id_elements[0], sizes_elements, 0); + boost::mpi::gatherv(M_comm, M_mesh.localEltIdWithGhostDefault(), &M_id_elements[0], sizes_elements, 0); } else { - boost::mpi::gatherv(M_comm, M_mesh.trianglesIdWithGhost(), 0); + boost::mpi::gatherv(M_comm, M_mesh.localEltIdWithGhostDefault(), 0); } // ------------------------------------------------------------- @@ -2811,7 +2840,7 @@ FiniteElement::diffuse(std::vector& variable_elt, double diffusivity_par void FiniteElement::scatterElementConnectivity() { - auto transfer_map_local = M_mesh.transferMapElt(); + auto transfer_map_default = M_mesh.transferMapEltDefault(); std::vector sizes_elements = M_sizes_elements_with_ghost; int nb_var_element = 3; @@ -2871,9 +2900,9 @@ FiniteElement::scatterElementConnectivity() if (neighbour_id_global_int>0) { - if (transfer_map_local.left.find(neighbour_id_global_int) != transfer_map_local.left.end()) + if (transfer_map_default.left.find(neighbour_id_global_int) != transfer_map_default.left.end()) { - int neighbour_id_local = transfer_map_local.left.find(neighbour_id_global_int)->second; + int neighbour_id_local = transfer_map_default.left.find(neighbour_id_global_int)->second; M_element_connectivity[nb_var_element*cpt+j] = neighbour_id_local; } else @@ -4000,7 +4029,7 @@ FiniteElement::update(std::vector const & UM_P) D_del_ci_ridge_myi[cpt] = 0.; // TO DO // if ( M_conc[cpt] > 1. ) // D_del_ci_ridge_myi[cpt] = 0.5*(1.-M_conc[cpt]); - } + } else { /* Here we conserve the area of multi-year ice: @@ -4119,13 +4148,13 @@ FiniteElement::update(std::vector const & UM_P) M_thick[cpt] = ((M_thick[cpt]>0.)?(M_thick[cpt] ):(0.)) ; M_thick_myi[cpt] = ((M_thick_myi[cpt]>0.)?(M_thick_myi[cpt] ):(0.)) ; M_snow_thick[cpt] = ((M_snow_thick[cpt]>0.)?(M_snow_thick[cpt]):(0.)) ; - /* This del_ci_ridge only works for ridge_myi_and_fyi=false*/ - D_del_ci_ridge_myi[cpt] = -M_conc_myi[cpt]; - if (newice_type == 4 && use_young_ice_in_myi_reset == true) + /* This del_ci_ridge only works for ridge_myi_and_fyi=false*/ + D_del_ci_ridge_myi[cpt] = -M_conc_myi[cpt]; + if (newice_type == 4 && use_young_ice_in_myi_reset == true) M_conc_myi[cpt] = std::max(0.,std::min(M_conc_myi[cpt],M_conc[cpt]+M_conc_young[cpt])); // Ensure M_conc_myi doesn't exceed total ice conc else M_conc_myi[cpt] = std::max(0.,std::min(M_conc_myi[cpt],M_conc[cpt])); // Ensure M_conc_myi doesn't exceed total ice conc - D_del_ci_ridge_myi[cpt]+=M_conc_myi[cpt]; + D_del_ci_ridge_myi[cpt]+=M_conc_myi[cpt]; }//loop over elements }//update @@ -5199,7 +5228,7 @@ FiniteElement::thermo(int dt) double const I_0 = vm["thermo.I_0"].as(); //! \param I_0 (double) Shortwave penetration into ice [fraction of total shortwave] bool use_young_ice_in_myi_reset = vm["age.include_young_ice"].as(); //! \param use_young_ice_in_myi_reset states if young ice should be included in the calculation of multiyear ice when it is reset (only if newice-type = 4) - const std::string date_string_reset_myi_md = vm["age.reset_date"].as(); //! \param date_string_reset_myi_md is the date (mmdd) of each year that the myi concentration should be reset to M_conc or M_conc+M_conc_young (this depends on use_young_ice_in_myi_reset) + const std::string date_string_reset_myi_md = vm["age.reset_date"].as(); //! \param date_string_reset_myi_md is the date (mmdd) of each year that the myi concentration should be reset to M_conc or M_conc+M_conc_young (this depends on use_young_ice_in_myi_reset) bool const reset_by_date = vm["age.reset_by_date"].as(); //! \param reset_by_date determines whether to reset myi on a certain date or by melt days double const freeze_days_threshold = vm["age.reset_freeze_days"].as(); //! \param freeze_days_threshold determines after how many days of freezing to reset myi, if reset_by_date is false bool equal_melting = vm["age.equal_melting"].as(); // decides if melting should affect myi and fyi the same or just fyi @@ -5254,7 +5283,7 @@ FiniteElement::thermo(int dt) std::vector albedo_young(M_num_elements); std::vector I_young(M_num_elements); if ( M_ice_cat_type==setup::IceCategoryType::YOUNG_ICE ) - { + { bulk_for_young = true ; this->IABulkFluxes(M_tsurf_young, M_hs_young, M_conc_young, Qia_young, Qlw_young, Qsw_young, Qlh_young, Qsh_young, @@ -5641,8 +5670,8 @@ FiniteElement::thermo(int dt) throw std::logic_error("Wrong melt_type"); } } - - // ICE AGE + + // ICE AGE if ( reset_by_date == false ) use_young_ice_in_myi_reset = false; @@ -5657,11 +5686,11 @@ FiniteElement::thermo(int dt) // Update M_freeze_days on last time step of the day if (step_in_day == num_steps_in_day) { - if (M_del_vi_tend[i] > 0.) // It's freezing - { + if (M_del_vi_tend[i] > 0.) // It's freezing + { M_freeze_days[i] += 1.; } - else if (M_del_vi_tend[i] < 0.) // It's melting + else if (M_del_vi_tend[i] < 0.) // It's melting { // melting occurring, so need to adjust to new onset M_freeze_days[i] = 0.; @@ -5765,7 +5794,7 @@ FiniteElement::thermo(int dt) #ifdef OASIS /* In case there is an FSD and M_conc>0: */ if ( (M_num_fsd_bins>0) && (melt_type==3) ) - { + { this->redistributeThermoFSD(i,ddt,lat_melt_rate,young_ice_growth,old_conc,old_conc_young); } /* In case there is melt_type!=3 (no FSD dependent lateral melting), FSD shape is unchanged. @@ -6037,7 +6066,7 @@ FiniteElement::thermo(int dt) } } - // Only reset if we have not already reset this season. So need to check if this is the case + // Only reset if we have not already reset this season. So need to check if this is the case if (date_string_md == "0801" && std::fmod(M_current_time, 1.) == 0.) { M_freeze_onset[i] = 0.; @@ -6061,7 +6090,7 @@ FiniteElement::thermo(int dt) // Now ensure that freeze and melt onsets are 0 or 1 M_freeze_onset[i] = std::round(M_freeze_onset[i]); - double old_conc_myi = M_conc_myi[i]; // delta= -old + new + double old_conc_myi = M_conc_myi[i]; // delta= -old + new double old_thick_myi = M_thick_myi[i]; double c_myi_max = M_conc[i]; double v_myi_max = M_thick[i]; @@ -6072,7 +6101,7 @@ FiniteElement::thermo(int dt) v_myi_max += M_h_young[i]; } - if (reset_myi) // + if (reset_myi) // { if (reset_by_date == false) { @@ -6090,7 +6119,7 @@ FiniteElement::thermo(int dt) } M_conc_myi[i] = std::max(0.,std::min(1.,M_conc_myi[i])); //make sure it doesn't exceed 1 (it shouldn't) M_thick_myi[i] = std::max(0.,M_thick_myi[i]); //make sure volume is positive - del_ci_rplnt_myi =M_conc_myi[i] - old_conc_myi; + del_ci_rplnt_myi =M_conc_myi[i] - old_conc_myi; del_vi_rplnt_myi =M_thick_myi[i] - old_thick_myi; } else // on a non-reset day, myi is only modified by melting, not freezing @@ -6099,21 +6128,21 @@ FiniteElement::thermo(int dt) double del_c_ratio = std::min(M_conc[i]/old_conc,1.); double del_v_ratio = std::min(M_thick[i]/old_vol,1.); if (del_v_ratio < 1.) // if there is some melt of old ice - { + { if (equal_melting) - { - del_ci_mlt_myi = std::min(0.,M_conc_myi[i]*(del_c_ratio-1.)); // <0 + { + del_ci_mlt_myi = std::min(0.,M_conc_myi[i]*(del_c_ratio-1.)); // <0 del_vi_mlt_myi = std::min(0.,M_thick_myi[i]*(del_v_ratio-1.)); // <0 } // Check it does not get crazy M_conc_myi[i] = std::max(0.,std::min(c_myi_max, M_conc_myi[i] +del_ci_mlt_myi)); M_thick_myi[i] = std::max(0.,std::min(v_myi_max, M_thick_myi[i]+del_vi_mlt_myi)); //If del_c removes all fyi, take some from myi. Preferentially removes fyi. Include young ice in total ice conc - //M_conc_myi[i] = std::max(0.,std::min(M_conc[i], M_conc_myi[i]*conc_loss_ratio)); - // Same logic for the volume + //M_conc_myi[i] = std::max(0.,std::min(M_conc[i], M_conc_myi[i]*conc_loss_ratio)); + // Same logic for the volume //M_thick_myi[i] = std::max(0.,std::min(M_thick[i], M_thick_myi[i]*thick_loss_ratio)); // // Recompute effective change - del_ci_mlt_myi = M_conc_myi[i] - old_conc_myi; + del_ci_mlt_myi = M_conc_myi[i] - old_conc_myi; del_vi_mlt_myi = M_thick_myi[i]- old_thick_myi; } } @@ -6463,12 +6492,12 @@ FiniteElement::meltPonds(const int cpt, const double dt, const double hi, // Make sure the pond depth isn't too high --> Holland et al. (2012) double pond_depth = std::min(dep2frac*D_pond_fraction[cpt],0.9*hi); - // If pond depth was too high, just drain the excess + // If pond depth was too high, just drain the excess M_pond_volume[cpt] = pond_depth*D_pond_fraction[cpt]; - + // Make sure the pond depth isn't microscopic pond_depth = std::max(0.05, pond_depth); - D_pond_fraction[cpt] = std::min(D_pond_fraction[cpt], + D_pond_fraction[cpt] = std::min(D_pond_fraction[cpt], (M_lid_volume[cpt]+M_pond_volume[cpt])/pond_depth); double delLidVolume = 0; // Volume increase is always positive! @@ -6964,6 +6993,7 @@ FiniteElement::init() //! - 10) Initialise timers M_timer = Timer(); + }//init // ============================================================================== @@ -7550,7 +7580,7 @@ FiniteElement::initOASIS() + std::string("Set setup.ocean-type to coupled or coupler.with_waves to true to activate the coupling.") ); M_cpl_out = GridOutput(bamgmesh, M_local_nelements, grid, nodal_variables, elemental_variables, vectorial_variables, - cpl_time_step*days_in_sec, true, bamgmesh_root, M_mesh.transferMapElt(), M_comm); + cpl_time_step*days_in_sec, true, bamgmesh_root, M_mesh.transferMapEltDefault(), M_comm); if ( M_ocean_type == setup::OceanType::COUPLED ) { @@ -7894,14 +7924,14 @@ FiniteElement::step() M_timer.tock("regrid"); #ifdef OASIS - /* Only M_cpl_out needs to provide M_mesh.transferMapElt and bamgmesh_root because these + /* Only M_cpl_out needs to provide M_mesh.transferMapEltDefault and bamgmesh_root because these * are needed iff we do conservative remapping and this is only supported in the coupled * case (so far). */ M_timer.tick("resetMeshMean_cpl"); if ( M_rank==0 ) - M_cpl_out.resetMeshMean(bamgmesh, M_regrid, M_local_nelements, M_mesh.transferMapElt(), bamgmesh_root); + M_cpl_out.resetMeshMean(bamgmesh, M_regrid, M_local_nelements, M_mesh.transferMapEltDefault(), bamgmesh_root); else - M_cpl_out.resetMeshMean(bamgmesh, M_regrid, M_local_nelements, M_mesh.transferMapElt()); + M_cpl_out.resetMeshMean(bamgmesh, M_regrid, M_local_nelements, M_mesh.transferMapEltDefault()); if ( M_ocean_type == setup::OceanType::COUPLED ) { @@ -7930,7 +7960,7 @@ FiniteElement::step() else #endif if ( vm["moorings.use_conservative_remapping"].as() ) - M_moorings.resetMeshMean(bamgmesh, M_regrid, M_local_nelements, M_mesh.transferMapElt(), bamgmesh_root); + M_moorings.resetMeshMean(bamgmesh, M_regrid, M_local_nelements, M_mesh.transferMapEltDefault(), bamgmesh_root); else M_moorings.resetMeshMean(bamgmesh, M_regrid, M_local_nelements); @@ -8139,7 +8169,7 @@ FiniteElement::step() //====================================================================== ++pcpt; M_current_time = time_init + pcpt*dtime_step/(24*3600.0); - + //====================================================================== //! 8) Does the post-processing, checks the output and updates moorings. //====================================================================== @@ -8375,6 +8405,7 @@ FiniteElement::run() // ********************************************************************** this->finalise(current_time_system); + }//run @@ -9180,7 +9211,7 @@ FiniteElement::initMoorings() "plat", "plon", "ptheta", GridOutput::interpMethod::conservative, false); M_moorings = GridOutput(bamgmesh, M_local_nelements, grid, nodal_variables, elemental_variables, vectorial_variables, - M_moorings_averaging_period, true, bamgmesh_root, M_mesh.transferMapElt(), M_comm); + M_moorings_averaging_period, true, bamgmesh_root, M_mesh.transferMapEltDefault(), M_comm); } else { // don't use conservative remapping GridOutput::Grid grid( Environment::vm()["moorings.grid_file"].as(), @@ -9201,7 +9232,7 @@ FiniteElement::initMoorings() "plat", "plon", "ptheta", GridOutput::interpMethod::conservative, false); M_moorings = GridOutput(bamgmesh, M_local_nelements, grid, nodal_variables, elemental_variables, vectorial_variables, - M_moorings_averaging_period, true, bamgmesh_root, M_mesh.transferMapElt(), M_comm); + M_moorings_averaging_period, true, bamgmesh_root, M_mesh.transferMapEltDefault(), M_comm); } #endif else @@ -10354,10 +10385,10 @@ FiniteElement::explicitSolve() } M_timer.tock("sub-solve"); - M_timer.tick("updateGhosts"); - // Update the ghosts and leave! - this->updateGhosts(M_VT); - M_timer.tock("updateGhosts"); + M_timer.tick("updateGhostNodes"); + // Update the ghost nodes and leave! + this->updateGhostNodes(M_VT); + M_timer.tock("updateGhostNodes"); // Move the mesh and update total displacement // For EVP and BBM we move the mesh every sub-time step @@ -10432,7 +10463,7 @@ FiniteElement::explicitSolve() M_VT[v_indx] /= num_neighbours; } - this->updateGhosts(M_VT); + this->updateGhostNodes(M_VT); } // Move the mesh in the open water part @@ -11342,7 +11373,7 @@ FiniteElement::initIce() boost::mpi::broadcast(M_comm, &random_number_root[0], M_mesh.numGlobalElements(), 0); // - now set on each processor - auto id_elements = M_mesh.trianglesIdWithGhost(); + auto id_elements = M_mesh.localEltIdWithGhostDefault(); for (int i=0; i1e-14) ? tmp_var : 0.; + M_conc[i] = (tmp_var>1e-14) ? tmp_var : 0.; tmp_var=M_init_thick[i]; M_thick[i] = tmp_var ; tmp_var=0.;//This is a test for the southern ocean - M_snow_thick[i] = (tmp_var>1e-14) ? tmp_var : 0.; + M_snow_thick[i] = (tmp_var>1e-14) ? tmp_var : 0.; //if either c or h equal zero, we set the others to zero as well if(M_conc[i]<=0.) @@ -13729,8 +13760,8 @@ FiniteElement::importBamg(BamgMesh const* bamg_mesh) void FiniteElement::createGraph() { - auto M_local_ghost = M_mesh.localGhost(); - auto M_transfer_map = M_mesh.transferMap(); + auto local_ghost = M_mesh.localGhost(); + auto transfer_map_dof = M_mesh.transferMapDof(); int Nd = bamgmesh->NodalConnectivitySize[1]; std::vector dz; @@ -13747,17 +13778,17 @@ FiniteElement::createGraph() int counter_onnz = 0; int Ncc = bamgmesh->NodalConnectivity[Nd*(i+1)-1]; - int gid = M_transfer_map.right.find(i+1)->second; - //if (std::find(M_local_ghost.begin(),M_local_ghost.end(),gid) == M_local_ghost.end()) - if (!std::binary_search(M_local_ghost.begin(),M_local_ghost.end(),gid)) + int gid = transfer_map_dof.right.find(i+1)->second; + //if (std::find(local_ghost.begin(),local_ghost.end(),gid) == local_ghost.end()) + if (!std::binary_search(local_ghost.begin(),local_ghost.end(),gid)) { for (int j=0; jNodalConnectivity[Nd*i+j]; - int gid2 = M_transfer_map.right.find(currentr)->second; - //if (std::find(M_local_ghost.begin(),M_local_ghost.end(),gid2) == M_local_ghost.end()) - if (!std::binary_search(M_local_ghost.begin(),M_local_ghost.end(),gid2)) + int gid2 = transfer_map_dof.right.find(currentr)->second; + //if (std::find(local_ghost.begin(),local_ghost.end(),gid2) == local_ghost.end()) + if (!std::binary_search(local_ghost.begin(),local_ghost.end(),gid2)) { ++counter_dnnz; } @@ -13766,7 +13797,7 @@ FiniteElement::createGraph() ++counter_onnz; } - //std::cout<<"Connect["<< j <<"]= "<< currentr << " or "<< M_transfer_map.right.find(currentr)->second <<"\n"; + //std::cout<<"Connect["<< j <<"]= "<< currentr << " or "<< transfer_map_dof.right.find(currentr)->second <<"\n"; } d_nnz.push_back(2*(counter_dnnz+1)); @@ -13821,7 +13852,7 @@ FiniteElement::createGraph() std::cout<<"WITH GHOST "<< index+1 <<"\n"; std::cout<<"************02************\n"; - for (int const& index : M_local_ghost) + for (int const& index : local_ghost) std::cout<<"GHOST "<< index <<"\n"; } #endif @@ -13832,149 +13863,284 @@ FiniteElement::createGraph() //! Updates the ghost nodes //! Called by the explicit solver void -FiniteElement::updateGhosts(std::vector& mesh_nodal_vec) +FiniteElement::updateGhostNodes(std::vector& mesh_nodal_vec) { std::vector> extract_local_values(M_comm.size()); - for (int i=0; i> ghost_update_values(M_comm.size()); - for (int const& proc : M_recipients_proc_id) + for (int const& proc : M_recipients_node_procid) M_comm.send(proc, M_rank, extract_local_values[proc]); - for (int const& proc : M_local_ghosts_proc_id) + for (int const& proc : M_ghost_node_procid) M_comm.recv(proc, proc, ghost_update_values[proc]); - for (int i=0; i> local_ghosts_global_index(M_comm.size()); + std::vector> ghost_global_indices(M_comm.size()); - for (int i=0; i " << M_transfer_map.left.find(currentid)->second - // <<" true ghost= "<< globalNumToprocId(currentid) <<"\n"; + // << " <-----> " << transfer_map_dof.left.find(currentid)->second + // <<" true ghost= "<< globalDofToProcId(currentid) <<"\n"; - local_ghosts_global_index[globalNumToprocId(currentid)].push_back(currentid); + ghost_global_indices[globalDofToProcId(currentid)].push_back(currentid); } - M_local_ghosts_proc_id.resize(0); - M_local_ghosts_local_index.resize(M_comm.size()); + M_ghost_node_procid.resize(0); + M_ghost_node_local_indices.resize(M_comm.size()); for (int i=0; i "<< i <<" \n"; - - for (int j=0; jsecond-1; - // std::cout<<" "<< M_local_ghosts_local_index[i][j] <<" "; + int currentindex = ghost_global_indices[i][j]; + M_ghost_node_local_indices[i][j] = transfer_map_dof.left.find(currentindex)->second-1; } - // std::cout<<"\n"; } - std::vector> recipients_proc_id_extended; - boost::mpi::all_gather(M_comm, M_local_ghosts_proc_id, recipients_proc_id_extended); + std::vector> recipients_procid_extended; + boost::mpi::all_gather(M_comm, M_ghost_node_procid, recipients_procid_extended); - M_recipients_proc_id.resize(0); + M_recipients_node_procid.resize(0); for (int i=0; i> extract_global_index(M_comm.size()); - for (int const& proc : M_local_ghosts_proc_id) - M_comm.send(proc, M_rank, local_ghosts_global_index[proc]); + for (int const& proc : M_ghost_node_procid) + M_comm.send(proc, M_rank, ghost_global_indices[proc]); - for (int const& proc : M_recipients_proc_id) + for (int const& proc : M_recipients_node_procid) M_comm.recv(proc, proc, extract_global_index[proc]); - M_extract_local_index.resize(M_comm.size()); + M_extract_local_node_indices.resize(M_comm.size()); for (int i=0; isecond-1; + M_extract_local_node_indices[i][j] = transfer_map_dof.left.find(extract_global_index[i][j])->second-1; } } -} //initUpdateGhosts +} //initUpdateGhostNodes // ------------------------------------------------------------------------------------- -//! Called by initUpdateGhosts +//! Called by initUpdateGhostNodes int -FiniteElement::globalNumToprocId(int global_num) +FiniteElement::globalDofToProcId(int global_dof) { int cpt = 0; for (int i=0; i& mesh_elt_ctr) +{ + std::vector> extract_local_values(M_comm.size()); + + for (int i=0; i> ghost_update_values(M_comm.size()); + + for (int const& proc : M_recipients_elt_procid) + M_comm.send(proc, M_rank, extract_local_values[proc]); + + for (int const& proc : M_ghost_elt_procid) + M_comm.recv(proc, proc, ghost_update_values[proc]); + + for (int i=0; i mesh_elt = mesh_elt_ctr; + this->updateGhostElements(mesh_elt); +} + +// ------------------------------------------------------------------------------------- +//! Initialise maps to update ghost elements +//! Called after remeshing and at init by distributedMeshProcessing +void +FiniteElement::initUpdateGhostElements() +{ + auto transfer_map_elt = M_mesh.transferMapElt(); + auto local_ghost_elt = M_mesh.localGhostElt(); + + std::vector> ghost_global_indices(M_comm.size()); + for (int i=0; isecond-1; + } + } + + std::vector> recipients_procid_extended; + boost::mpi::all_gather(M_comm, M_ghost_elt_procid, recipients_procid_extended); + + M_recipients_elt_procid.resize(0); + + for (int i=0; i> extract_global_index(M_comm.size()); + + for (int const& proc : M_ghost_elt_procid) + { + M_comm.send(proc, M_rank, ghost_global_indices[proc]); + } + + for (int const& proc : M_recipients_elt_procid) + { + M_comm.recv(proc, proc, extract_global_index[proc]); + } + + M_extract_local_elt_indices.resize(M_comm.size()); + + for (int i=0; isecond-1; + } + } + +} //initUpdateGhostElements + +// ------------------------------------------------------------------------------------- +//! Called by initUpdateGhostElements +int +FiniteElement::globalEltIdToProcId(int global_eltid) +{ + int cpt = 0; + for (int i=0; i graphmpi_ptrtype; - // typedef ExternalData external_data; - // typedef ExternalData::Dataset Dataset; - // typedef ExternalData::Grid Grid; - // typedef ExternalData::Dimension Dimension; - // typedef ExternalData::Variable Variable; - // typedef ExternalData::Vectorial_Variable Vectorial_Variable; - // typedef boost::ptr_vector externaldata_ptr_vector; - typedef DataSet Dataset; typedef ExternalData external_data; typedef typename std::vector external_data_vec ; @@ -89,8 +81,6 @@ class FiniteElement FiniteElement(Communicator const& comm = Environment::comm()); - // FiniteElement(Communicator const& comm = Environment::comm()); - mesh_type const& mesh() const {return M_mesh;} void initMesh(); @@ -143,7 +133,6 @@ class FiniteElement void gatherFieldsElement(std::vector& interp_in_elements); void scatterFieldsElement(double* interp_elt_out); - //void gatherUM(std::vector& um); void gatherNodalField(std::vector const& field_local, std::vector& field_root); void scatterNodalField(std::vector const& field_root, std::vector& field_local); @@ -281,7 +270,7 @@ class FiniteElement void PwlInterp2D(); void importBamg(BamgMesh const* bamg_mesh); - void createGraph();//(BamgMesh const* bamg_mesh); + void createGraph(); void assignVariables(); void initVariables(); void calcAuxiliaryVariables(); @@ -296,9 +285,14 @@ class FiniteElement void update(std::vector const & UM_P); void updateSigmaDamage(double const dt); - void updateGhosts(std::vector& mesh_nodal_vec); - void initUpdateGhosts(); - int globalNumToprocId(int global_num); + void updateGhostNodes(std::vector& mesh_nodal_vec); + void initUpdateGhostNodes(); + int globalDofToProcId(int global_dof); + + void updateGhostElements(std::vector& mesh_elt_ctr); + void updateGhostElements(ModelVariable& mesh_elt_ctr); + void initUpdateGhostElements(); + int globalEltIdToProcId(int global_eltid); #ifdef OASIS bool M_couple_waves; @@ -388,7 +382,6 @@ class FiniteElement mesh_type M_mesh_previous; std::map M_nodes; - //std::vector M_nodes; std::vector M_edges; std::vector M_elements; @@ -397,7 +390,7 @@ class FiniteElement int M_ndof; int M_local_ndof; - int M_local_ndof_ghost; + //int M_local_ndof_ghost; int M_local_nelements; int M_rank; Communicator M_comm; @@ -610,10 +603,15 @@ class FiniteElement std::string M_export_path; private: // update solution from explicit solver - std::vector> M_extract_local_index; - std::vector M_recipients_proc_id; - std::vector M_local_ghosts_proc_id; - std::vector> M_local_ghosts_local_index; + std::vector> M_extract_local_node_indices; + std::vector M_recipients_node_procid; + std::vector M_ghost_node_procid; + std::vector> M_ghost_node_local_indices; + + std::vector> M_extract_local_elt_indices; + std::vector M_recipients_elt_procid; + std::vector M_ghost_elt_procid; + std::vector> M_ghost_elt_local_indices; private: From 2a6ffda6ed728de1fc4fc98fa81aa3c69a9519c8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Thu, 31 Oct 2024 08:48:05 +0100 Subject: [PATCH 2/2] Make updateGhostElements accept an array pointer This way the overloaded functions which accept an std::vector and a ModelVariable can call the basis function with &array[0] as an argument. --- model/finiteelement.cpp | 12 ++++++++++-- model/finiteelement.hpp | 1 + 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/model/finiteelement.cpp b/model/finiteelement.cpp index 5edd8bf56..e6f510c5f 100644 --- a/model/finiteelement.cpp +++ b/model/finiteelement.cpp @@ -14008,6 +14008,15 @@ FiniteElement::globalDofToProcId(int global_dof) //! Called by the explicit solver void FiniteElement::updateGhostElements(std::vector& mesh_elt_ctr) +{ + this->updateGhostElements(&mesh_elt_ctr[0]); +} + +// ------------------------------------------------------------------------------------- +//! Updates the ghost elements +//! Called by the explicit solver +void +FiniteElement::updateGhostElements(double* mesh_elt_ctr) { std::vector> extract_local_values(M_comm.size()); @@ -14045,8 +14054,7 @@ FiniteElement::updateGhostElements(std::vector& mesh_elt_ctr) void FiniteElement::updateGhostElements(ModelVariable& mesh_elt_ctr) { - std::vector mesh_elt = mesh_elt_ctr; - this->updateGhostElements(mesh_elt); + this->updateGhostElements(&mesh_elt_ctr[0]); } // ------------------------------------------------------------------------------------- diff --git a/model/finiteelement.hpp b/model/finiteelement.hpp index 5e36c0656..aa5f7a147 100644 --- a/model/finiteelement.hpp +++ b/model/finiteelement.hpp @@ -289,6 +289,7 @@ class FiniteElement void initUpdateGhostNodes(); int globalDofToProcId(int global_dof); + void updateGhostElements(double* mesh_elt_ctr); void updateGhostElements(std::vector& mesh_elt_ctr); void updateGhostElements(ModelVariable& mesh_elt_ctr); void initUpdateGhostElements();