36 #include "DistributedTetrahedralMesh.hpp"
43 #include <boost/scoped_array.hpp>
46 #include "Element.hpp"
47 #include "BoundaryElement.hpp"
50 #include "DistributedVectorFactory.hpp"
51 #include "OutputFileHandler.hpp"
52 #include "NodePartitioner.hpp"
54 #include "RandomNumberGenerator.hpp"
57 #include "TetrahedralMesh.hpp"
58 #include "Warnings.hpp"
62 #if (PARMETIS_MAJOR_VERSION >= 4) //ParMETIS 4.x and above
74 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
77 mTotalNumElements(0u),
78 mTotalNumBoundaryElements(0u),
80 mpSpaceRegion(nullptr),
81 mPartitioning(partitioningMethod)
83 if (ELEMENT_DIM == 1 && (partitioningMethod != DistributedTetrahedralMeshPartitionType::GEOMETRIC))
86 mPartitioning = DistributedTetrahedralMeshPartitionType::DUMB;
90 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
93 for (
unsigned i=0; i<this->mHaloNodes.size(); i++)
95 delete this->mHaloNodes[i];
99 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
103 mPartitioning = DistributedTetrahedralMeshPartitionType::DUMB;
106 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
109 std::set<unsigned>& rNodesOwned,
110 std::set<unsigned>& rHaloNodesOwned,
111 std::set<unsigned>& rElementsOwned,
112 std::vector<unsigned>& rProcessorsOffset)
114 if (mPartitioning == DistributedTetrahedralMeshPartitionType::METIS_LIBRARY)
116 WARNING(
"METIS partitioning is deprecated. Switching to parMETIS");
117 mPartitioning = DistributedTetrahedralMeshPartitionType::PARMETIS_LIBRARY;
119 if (mPartitioning == DistributedTetrahedralMeshPartitionType::PETSC_MAT_PARTITION && !
PetscTools::HasParMetis())
123 WARNING(
"PETSc/parMETIS partitioning requires PETSc to be configured with parMETIS as an option. Current install has PETSc and parMETIS installed independently. Switching to parMETIS");
124 mPartitioning = DistributedTetrahedralMeshPartitionType::PARMETIS_LIBRARY;
134 ParMetisLibraryNodeAndElementPartitioning(rMeshReader, rElementsOwned, rNodesOwned, rHaloNodesOwned, rProcessorsOffset);
141 if (mPartitioning==DistributedTetrahedralMeshPartitionType::PETSC_MAT_PARTITION &&
PetscTools::IsParallel())
149 EXCEPTION(
"Using GEOMETRIC partition for DistributedTetrahedralMesh with local regions not set. Call SetProcessRegion(ChasteCuboid)");
162 for (std::set<unsigned>::iterator iter = rNodesOwned.begin();
163 iter != rNodesOwned.end();
167 rElementsOwned.insert( containing_elements.begin(), containing_elements.end() );
172 std::set<unsigned> node_index_set;
174 for (std::set<unsigned>::iterator iter = rElementsOwned.begin();
175 iter != rElementsOwned.end();
184 std::set_difference(node_index_set.begin(), node_index_set.end(),
185 rNodesOwned.begin(), rNodesOwned.end(),
186 std::inserter(rHaloNodesOwned, rHaloNodesOwned.begin()));
190 for (
unsigned element_number = 0; element_number < mTotalNumElements; element_number++)
194 bool element_owned =
false;
195 std::set<unsigned> temp_halo_nodes;
197 for (std::vector<unsigned>::const_iterator it = element_data.
NodeIndices.begin();
201 if (rNodesOwned.find(*it) != rNodesOwned.end())
203 element_owned =
true;
204 rElementsOwned.insert(element_number);
208 temp_halo_nodes.insert(*it);
214 rHaloNodesOwned.insert(temp_halo_nodes.begin(), temp_halo_nodes.end());
219 if (mPartitioning==DistributedTetrahedralMeshPartitionType::PETSC_MAT_PARTITION &&
PetscTools::IsParallel())
231 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
235 std::set<unsigned> nodes_owned;
236 std::set<unsigned> halo_nodes_owned;
237 std::set<unsigned> elements_owned;
238 std::vector<unsigned> proc_offsets;
242 mTotalNumBoundaryElements = rMeshReader.
GetNumFaces();
248 ComputeMeshPartitioning(rMeshReader, nodes_owned, halo_nodes_owned, elements_owned, proc_offsets);
253 this->mElements.reserve(elements_owned.size());
254 this->mNodes.reserve(nodes_owned.size());
260 std::vector<double> coords;
267 unsigned global_node_index = node_it.GetIndex();
269 RegisterNode(global_node_index);
279 this->mNodes.push_back(p_node);
286 unsigned global_node_index = node_it.
GetIndex();
288 RegisterHaloNode(global_node_index);
289 mHaloNodes.push_back(
new Node<SPACE_DIM>(global_node_index, coords,
false));
296 for (
unsigned node_index=0; node_index < mTotalNumNodes; node_index++)
298 std::vector<double> coords;
303 if (nodes_owned.find(node_index) != nodes_owned.end())
305 RegisterNode(node_index);
314 this->mNodes.push_back(p_node);
318 if (halo_nodes_owned.find(node_index) != halo_nodes_owned.end())
320 RegisterHaloNode(node_index);
332 unsigned global_element_index = elem_it.GetIndex();
334 std::vector<Node<SPACE_DIM>*> nodes;
335 for (
unsigned j=0; j<ELEMENT_DIM+1; j++)
338 nodes.push_back(this->GetNodeOrHaloNode(element_data.
NodeIndices[j]));
341 RegisterElement(global_element_index);
343 this->mElements.push_back(p_element);
356 for (
unsigned face_index=0; face_index<mTotalNumBoundaryElements; face_index++)
359 std::vector<unsigned> node_indices = face_data.
NodeIndices;
363 for (
unsigned node_index=0; node_index<node_indices.size(); node_index++)
366 if (mNodesMapping.find(node_indices[node_index]) != mNodesMapping.end())
380 std::vector<Node<SPACE_DIM>*> nodes;
382 for (
unsigned node_index=0; node_index<node_indices.size(); node_index++)
388 nodes.push_back(this->GetNodeOrHaloNode(node_indices[node_index]));
392 EXCEPTION(
"Face does not appear in element file (Face " << face_index <<
" in "<<this->mMeshFileBaseName<<
")");
398 for (
unsigned j=0; j<nodes.size(); j++)
400 if (!nodes[j]->IsBoundaryNode())
402 nodes[j]->SetAsBoundaryNode();
403 this->mBoundaryNodes.push_back(nodes[j]);
406 nodes[j]->AddBoundaryElement(face_index);
409 RegisterBoundaryElement(face_index);
411 this->mBoundaryElements.push_back(p_boundary_element);
429 assert(this->mNodePermutation.size() != 0);
440 num_owned = proc_offsets[rank+1]-proc_offsets[rank];
444 num_owned = mTotalNumNodes - proc_offsets[rank];
447 assert(!this->mpDistributedVectorFactory);
453 assert(this->mpDistributedVectorFactory);
465 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
468 return this->mNodes.size();
471 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
474 return this->mHaloNodes.size();
477 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
480 return this->mElements.size();
483 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
486 return this->mBoundaryElements.size();
489 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
492 return mTotalNumNodes;
495 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
498 return mTotalNumNodes;
501 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
504 return mTotalNumElements;
507 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
510 return mPartitioning;
513 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
516 return mTotalNumBoundaryElements;
519 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
523 rHaloIndices.clear();
524 for (
unsigned i=0; i<mHaloNodes.size(); i++)
526 rHaloIndices.push_back(mHaloNodes[i]->GetIndex());
530 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
533 if (mpSpaceRegion ==
nullptr)
535 EXCEPTION(
"Trying to get unset mpSpaceRegion");
537 return mpSpaceRegion;
540 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
547 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
550 mpSpaceRegion = pRegion;
553 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
566 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
579 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
582 mNodesMapping[index] = this->mNodes.size();
585 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
588 mHaloNodesMapping[index] = mHaloNodes.size();
591 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
594 mElementsMapping[index] = this->mElements.size();
597 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
600 mBoundaryElementsMapping[index] = this->mBoundaryElements.size();
603 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
606 std::map<unsigned, unsigned>::const_iterator node_position = mNodesMapping.find(index);
608 if (node_position == mNodesMapping.end())
612 return node_position->second;
622 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
625 std::map<unsigned, unsigned>::const_iterator element_position = mElementsMapping.find(index);
627 if (element_position == mElementsMapping.end())
632 return element_position->second;
635 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
638 std::map<unsigned, unsigned>::const_iterator boundary_element_position = mBoundaryElementsMapping.find(index);
640 if (boundary_element_position == mBoundaryElementsMapping.end())
645 return boundary_element_position->second;
648 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
651 std::map<unsigned, unsigned>::const_iterator node_position;
653 if ((node_position=mHaloNodesMapping.find(index)) != mHaloNodesMapping.end())
655 return mHaloNodes[node_position->second];
658 if ((node_position=mNodesMapping.find(index)) != mNodesMapping.end())
661 return this->mNodes[node_position->second];
667 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
673 mNodesMapping.clear();
674 mHaloNodesMapping.clear();
677 for (
unsigned index=0; index<this->mNodes.size(); index++)
679 unsigned old_index = this->mNodes[index]->GetIndex();
680 unsigned new_index = this->mNodePermutation[old_index];
682 this->mNodes[index]->SetIndex(new_index);
683 mNodesMapping[new_index] = index;
686 for (
unsigned index=0; index<mHaloNodes.size(); index++)
688 unsigned old_index = mHaloNodes[index]->GetIndex();
689 unsigned new_index = this->mNodePermutation[old_index];
691 mHaloNodes[index]->SetIndex(new_index);
692 mHaloNodesMapping[new_index] = index;
696 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
699 assert(ELEMENT_DIM == 1);
704 EXCEPTION(
"There aren't enough nodes to make parallelisation worthwhile");
708 if (mPartitioning == DistributedTetrahedralMeshPartitionType::GEOMETRIC)
712 EXCEPTION(
"Space region not set for GEOMETRIC partition of DistributedTetrahedralMesh");
730 this->ConstructFromMeshReader(mesh_reader);
735 mPartitioning=DistributedTetrahedralMeshPartitionType::DUMB;
736 mTotalNumNodes=width+1;
737 mTotalNumBoundaryElements=2u;
738 mTotalNumElements=width;
741 assert(!this->mpDistributedVectorFactory);
743 if (this->mpDistributedVectorFactory->GetLocalOwnership() == 0)
754 bool am_top_most = (this->mpDistributedVectorFactory->GetHigh() == mTotalNumNodes);
756 unsigned lo_node=this->mpDistributedVectorFactory->GetLow();
757 unsigned hi_node=this->mpDistributedVectorFactory->GetHigh();
769 for (
unsigned node_index=lo_node; node_index<hi_node; node_index++)
773 if (node_index<this->mpDistributedVectorFactory->GetLow() ||
774 node_index==this->mpDistributedVectorFactory->GetHigh() )
777 RegisterHaloNode(node_index);
778 mHaloNodes.push_back(p_node);
782 RegisterNode(node_index);
783 this->mNodes.push_back(p_node);
789 this->mBoundaryNodes.push_back(p_node);
790 RegisterBoundaryElement(0);
793 if (node_index==width)
795 this->mBoundaryNodes.push_back(p_node);
796 RegisterBoundaryElement(1);
800 if (node_index>lo_node)
802 std::vector<Node<SPACE_DIM>*> nodes;
803 nodes.push_back(p_old_node);
804 nodes.push_back(p_node);
805 RegisterElement(node_index-1);
814 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
817 assert(SPACE_DIM == 2);
818 assert(ELEMENT_DIM == 2);
822 EXCEPTION(
"There aren't enough nodes to make parallelisation worthwhile");
826 if (mPartitioning == DistributedTetrahedralMeshPartitionType::GEOMETRIC)
830 EXCEPTION(
"Space region not set for GEOMETRIC partition of DistributedTetrahedralMesh");
848 this->ConstructFromMeshReader(mesh_reader);
853 mPartitioning=DistributedTetrahedralMeshPartitionType::DUMB;
855 mTotalNumNodes=(width+1)*(height+1);
856 mTotalNumBoundaryElements=(width+height)*2;
857 mTotalNumElements=width*height*2;
861 unsigned lo_y = y_partition.
GetLow();
862 unsigned hi_y = y_partition.
GetHigh();
864 assert(!this->mpDistributedVectorFactory);
866 if (this->mpDistributedVectorFactory->GetLocalOwnership() == 0)
877 bool am_top_most = (this->mpDistributedVectorFactory->GetHigh() == mTotalNumNodes);
892 for (
unsigned j=lo_y; j<hi_y; j++)
894 for (
unsigned i=0; i<width+1; i++)
896 bool is_boundary=
false;
897 if (i==0 || j==0 || i==width || j==height)
901 unsigned global_node_index=((width+1)*(j) + i);
906 RegisterHaloNode(global_node_index);
907 mHaloNodes.push_back(p_node);
911 RegisterNode(global_node_index);
912 this->mNodes.push_back(p_node);
916 this->mBoundaryNodes.push_back(p_node);
922 unsigned belem_index;
926 for (
unsigned i=0; i<width; i++)
928 std::vector<Node<SPACE_DIM>*> nodes;
929 nodes.push_back(GetNodeOrHaloNode( height*(width+1)+i+1 ));
930 nodes.push_back(GetNodeOrHaloNode( height*(width+1)+i ));
932 RegisterBoundaryElement(belem_index);
938 for (
unsigned j=lo_y+1; j<hi_y; j++)
940 std::vector<Node<SPACE_DIM>*> nodes;
941 nodes.push_back(GetNodeOrHaloNode( (width+1)*j-1 ));
942 nodes.push_back(GetNodeOrHaloNode( (width+1)*(j+1)-1 ));
943 belem_index=width+j-1;
944 RegisterBoundaryElement(belem_index);
951 for (
unsigned i=0; i<width; i++)
953 std::vector<Node<SPACE_DIM>*> nodes;
954 nodes.push_back(GetNodeOrHaloNode( i ));
955 nodes.push_back(GetNodeOrHaloNode( i+1 ));
956 belem_index=width+height+i;
957 RegisterBoundaryElement(belem_index);
963 for (
unsigned j=lo_y; j<hi_y-1; j++)
965 std::vector<Node<SPACE_DIM>*> nodes;
966 nodes.push_back(GetNodeOrHaloNode( (width+1)*(j+1) ));
967 nodes.push_back(GetNodeOrHaloNode( (width+1)*(j) ));
968 belem_index=2*width+height+j;
969 RegisterBoundaryElement(belem_index);
976 for (
unsigned j=lo_y; j<hi_y-1; j++)
978 for (
unsigned i=0; i<width; i++)
980 unsigned parity=(i+(height-j))%2;
981 unsigned nw=(j+1)*(width+1)+i;
982 unsigned sw=(j)*(width+1)+i;
983 std::vector<Node<SPACE_DIM>*> upper_nodes;
984 upper_nodes.push_back(GetNodeOrHaloNode( nw ));
985 upper_nodes.push_back(GetNodeOrHaloNode( nw+1 ));
986 if (stagger==
false || parity == 1)
988 upper_nodes.push_back(GetNodeOrHaloNode( sw+1 ));
992 upper_nodes.push_back(GetNodeOrHaloNode( sw ));
994 elem_index=2*(j*width+i);
995 RegisterElement(elem_index);
997 std::vector<Node<SPACE_DIM>*> lower_nodes;
998 lower_nodes.push_back(GetNodeOrHaloNode( sw+1 ));
999 lower_nodes.push_back(GetNodeOrHaloNode( sw ));
1000 if (stagger==
false ||parity == 1)
1002 lower_nodes.push_back(GetNodeOrHaloNode( nw ));
1006 lower_nodes.push_back(GetNodeOrHaloNode( nw+1 ));
1009 RegisterElement(elem_index);
1017 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1022 assert(SPACE_DIM == 3);
1023 assert(ELEMENT_DIM == 3);
1027 EXCEPTION(
"There aren't enough nodes to make parallelisation worthwhile");
1031 if (mPartitioning == DistributedTetrahedralMeshPartitionType::GEOMETRIC)
1035 EXCEPTION(
"Space region not set for GEOMETRIC partition of DistributedTetrahedralMesh");
1053 this->ConstructFromMeshReader(mesh_reader);
1058 mPartitioning=DistributedTetrahedralMeshPartitionType::DUMB;
1060 mTotalNumNodes=(width+1)*(height+1)*(depth+1);
1061 mTotalNumBoundaryElements=((width*height)+(width*depth)+(height*depth))*4;
1062 mTotalNumElements=width*height*depth*6;
1066 unsigned lo_z = z_partition.
GetLow();
1067 unsigned hi_z = z_partition.
GetHigh();
1070 assert(!this->mpDistributedVectorFactory);
1072 if (this->mpDistributedVectorFactory->GetLocalOwnership() == 0)
1083 bool am_top_most = (this->mpDistributedVectorFactory->GetHigh() == mTotalNumNodes);
1097 unsigned global_node_index;
1098 for (
unsigned k=lo_z; k<hi_z; k++)
1100 for (
unsigned j=0; j<height+1; j++)
1102 for (
unsigned i=0; i<width+1; i++)
1104 bool is_boundary =
false;
1105 if (i==0 || j==0 || k==0 || i==width || j==height || k==depth)
1109 global_node_index = (k*(height+1)+j)*(width+1)+i;
1116 RegisterHaloNode(global_node_index);
1117 mHaloNodes.push_back(p_node);
1121 RegisterNode(global_node_index);
1122 this->mNodes.push_back(p_node);
1127 this->mBoundaryNodes.push_back(p_node);
1135 unsigned element_nodes[6][4] = {{0, 1, 5, 7}, {0, 1, 3, 7},
1136 {0, 2, 3, 7}, {0, 2, 6, 7},
1137 {0, 4, 6, 7}, {0, 4, 5, 7}};
1138 std::vector<Node<SPACE_DIM>*> tetrahedra_nodes;
1140 for (
unsigned k=lo_z; k<hi_z-1; k++)
1142 unsigned belem_index = 0;
1146 belem_index = 2*(height*width+k*2*(height+width));
1149 for (
unsigned j=0; j<height; j++)
1151 for (
unsigned i=0; i<width; i++)
1154 unsigned global_node_indices[8];
1155 unsigned local_node_index = 0;
1157 for (
unsigned z = 0; z < 2; z++)
1159 for (
unsigned y = 0; y < 2; y++)
1161 for (
unsigned x = 0; x < 2; x++)
1163 global_node_indices[local_node_index] = i+x+(width+1)*(j+y+(height+1)*(k+z));
1170 for (
unsigned m = 0; m < 6; m++)
1174 tetrahedra_nodes.clear();
1176 for (
unsigned n = 0; n < 4; n++)
1178 tetrahedra_nodes.push_back(GetNodeOrHaloNode( global_node_indices[element_nodes[m][n]] ));
1180 unsigned elem_index = 6 * ((k*height+j)*width+i)+m;
1181 RegisterElement(elem_index);
1186 std::vector<Node<SPACE_DIM>*> triangle_nodes;
1190 triangle_nodes.clear();
1191 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[0] ));
1192 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[2] ));
1193 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[6] ));
1194 RegisterBoundaryElement(belem_index);
1196 triangle_nodes.clear();
1197 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[0] ));
1198 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[6] ));
1199 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[4] ));
1200 RegisterBoundaryElement(belem_index);
1205 triangle_nodes.clear();
1206 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[1] ));
1207 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[5] ));
1208 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[7] ));
1209 RegisterBoundaryElement(belem_index);
1211 triangle_nodes.clear();
1212 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[1] ));
1213 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[7] ));
1214 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[3] ));
1215 RegisterBoundaryElement(belem_index);
1220 triangle_nodes.clear();
1221 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[0] ));
1222 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[5] ));
1223 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[1] ));
1224 RegisterBoundaryElement(belem_index);
1226 triangle_nodes.clear();
1227 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[0] ));
1228 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[4] ));
1229 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[5] ));
1230 RegisterBoundaryElement(belem_index);
1235 triangle_nodes.clear();
1236 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[2] ));
1237 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[3] ));
1238 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[7] ));
1239 RegisterBoundaryElement(belem_index);
1241 triangle_nodes.clear();
1242 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[2] ));
1243 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[7] ));
1244 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[6] ));
1245 RegisterBoundaryElement(belem_index);
1250 triangle_nodes.clear();
1251 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[0] ));
1252 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[3] ));
1253 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[2] ));
1254 RegisterBoundaryElement(belem_index);
1256 triangle_nodes.clear();
1257 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[0] ));
1258 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[1] ));
1259 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[3] ));
1260 RegisterBoundaryElement(belem_index);
1265 triangle_nodes.clear();
1266 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[4] ));
1267 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[7] ));
1268 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[5] ));
1269 RegisterBoundaryElement(belem_index);
1271 triangle_nodes.clear();
1272 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[4] ));
1273 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[6] ));
1274 triangle_nodes.push_back(GetNodeOrHaloNode( global_node_indices[7] ));
1275 RegisterBoundaryElement(belem_index);
1284 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1290 for (
unsigned i=0; i<mHaloNodes.size(); i++)
1292 c_vector<double, SPACE_DIM>& r_location = mHaloNodes[i]->rGetModifiableLocation();
1295 r_location[2] *= zFactor;
1299 r_location[1] *= yFactor;
1301 r_location[0] *= xFactor;
1307 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1310 std::set<unsigned>& rElementsOwned,
1311 std::set<unsigned>& rNodesOwned,
1312 std::set<unsigned>& rHaloNodesOwned,
1313 std::vector<unsigned>& rProcessorsOffset)
1316 assert(ELEMENT_DIM==2 || ELEMENT_DIM==3);
1325 boost::scoped_array<idxtype> element_distribution(
new idxtype[num_procs+1]);
1326 boost::scoped_array<int> element_counts(
new int[num_procs]);
1328 element_distribution[0] = 0;
1330 for (
unsigned proc_index=1; proc_index<num_procs; proc_index++)
1332 element_distribution[proc_index] = element_distribution[proc_index-1] + num_elements/num_procs;
1333 element_counts[proc_index-1] = element_distribution[proc_index] - element_distribution[proc_index-1];
1336 element_distribution[num_procs] = num_elements;
1337 element_counts[num_procs-1] = element_distribution[num_procs] - element_distribution[num_procs-1];
1342 idxtype first_local_element = element_distribution[local_proc_index];
1343 idxtype last_plus_one_element = element_distribution[local_proc_index+1];
1344 idxtype num_local_elements = last_plus_one_element - first_local_element;
1346 boost::scoped_array<idxtype> eind(
new idxtype[num_local_elements*(ELEMENT_DIM+1)]);
1347 boost::scoped_array<idxtype> eptr(
new idxtype[num_local_elements+1]);
1357 for (idxtype element_index = 0; element_index < first_local_element; element_index++)
1363 unsigned counter = 0;
1364 for (idxtype element_index = 0; element_index < num_local_elements; element_index++)
1370 eptr[element_index] = counter;
1371 for (
unsigned i=0; i<ELEMENT_DIM+1; i++)
1376 eptr[num_local_elements] = counter;
1378 rMeshReader.
Reset();
1380 idxtype numflag = 0;
1388 idxtype ncommonnodes = 3;
1389 if (ELEMENT_DIM == 2)
1394 MPI_Comm communicator = PETSC_COMM_WORLD;
1400 ParMETIS_V3_Mesh2Dual(element_distribution.get(), eptr.get(), eind.get(),
1401 &numflag, &ncommonnodes, &xadj, &adjncy, &communicator);
1408 idxtype weight_flag = 0;
1409 idxtype n_constraints = 1;
1414 boost::scoped_array<real_t> tpwgts(
new real_t[n_subdomains]);
1415 real_t ubvec_value = (real_t)1.05;
1418 tpwgts[proc] = ((real_t)1.0)/n_subdomains;
1420 boost::scoped_array<idxtype> local_partition(
new idxtype[num_local_elements]);
1436 ParMETIS_V3_PartKway(element_distribution.get(), xadj, adjncy,
nullptr,
nullptr, &weight_flag, &numflag,
1437 &n_constraints, &n_subdomains, tpwgts.get(), &ubvec_value,
1438 options, &edgecut, local_partition.get(), &communicator);
1442 boost::scoped_array<idxtype> global_element_partition(
new idxtype[num_elements]);
1445 MPI_Datatype mpi_idxtype = MPI_LONG_LONG_INT;
1446 if (
sizeof(idxtype) ==
sizeof(
int))
1448 mpi_idxtype = MPI_INT;
1450 boost::scoped_array<int> int_element_distribution(
new int[num_procs+1]);
1451 for (
unsigned i=0; i<num_procs+1; ++i)
1453 int_element_distribution[i] = element_distribution[i];
1455 MPI_Allgatherv(local_partition.get(), num_local_elements, mpi_idxtype,
1456 global_element_partition.get(), element_counts.get(), int_element_distribution.get(), mpi_idxtype, PETSC_COMM_WORLD);
1458 local_partition.reset();
1460 for (
unsigned elem_index=0; elem_index<num_elements; elem_index++)
1462 if ((
unsigned) global_element_partition[elem_index] == local_proc_index)
1464 rElementsOwned.insert(elem_index);
1468 rMeshReader.
Reset();
1475 std::vector<unsigned> global_node_partition(num_nodes, UNASSIGNED_NODE);
1477 assert(rProcessorsOffset.size() == 0);
1478 rProcessorsOffset.resize(PetscTools::GetNumProcs(), 0);
1503 for (
unsigned element_number = 0; element_number < mTotalNumElements; element_number++)
1505 unsigned element_owner = global_element_partition[element_number];
1511 for (std::vector<unsigned>::const_iterator node_it = element_data.
NodeIndices.begin();
1519 if (global_node_partition[*node_it] == UNASSIGNED_NODE)
1521 if (element_owner == local_proc_index)
1523 rNodesOwned.insert(*node_it);
1526 global_node_partition[*node_it] = element_owner;
1533 rProcessorsOffset[proc]++;
1538 if (element_owner == local_proc_index)
1541 if (global_node_partition[*node_it] != local_proc_index)
1543 rHaloNodesOwned.insert(*node_it);
1563 rElementsOwned.clear();
1564 rMeshReader.
Reset();
1565 for (
unsigned element_number = 0; element_number < mTotalNumElements; element_number++)
1569 bool element_owned =
false;
1570 std::set<unsigned> temp_halo_nodes;
1572 for (std::vector<unsigned>::const_iterator node_it = element_data.
NodeIndices.begin();
1576 if (rNodesOwned.find(*node_it) != rNodesOwned.end())
1578 element_owned =
true;
1579 rElementsOwned.insert(element_number);
1583 temp_halo_nodes.insert(*node_it);
1589 rHaloNodesOwned.insert(temp_halo_nodes.begin(), temp_halo_nodes.end());
1593 rMeshReader.
Reset();
1598 std::vector<unsigned> local_index(PetscTools::GetNumProcs(), 0);
1600 this->mNodePermutation.resize(this->GetNumNodes());
1602 for (
unsigned node_index=0; node_index<this->GetNumNodes(); node_index++)
1604 unsigned partition = global_node_partition[node_index];
1605 assert(partition != UNASSIGNED_NODE);
1607 this->mNodePermutation[node_index] = rProcessorsOffset[partition] + local_index[partition];
1609 local_index[partition]++;
1613 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1636 c_vector<double, SPACE_DIM> global_minimum_point;
1637 c_vector<double, SPACE_DIM> global_maximum_point;
1638 MPI_Allreduce(&my_minimum_point.
rGetLocation()[0], &global_minimum_point[0], SPACE_DIM, MPI_DOUBLE, MPI_MIN, PETSC_COMM_WORLD);
1639 MPI_Allreduce(&my_maximum_point.
rGetLocation()[0], &global_maximum_point[0], SPACE_DIM, MPI_DOUBLE, MPI_MAX, PETSC_COMM_WORLD);
1647 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1654 double best_node_point_distance = DBL_MAX;
1655 if (best_node_index != UINT_MAX)
1657 best_node_point_distance = norm_2(this->GetNode(best_node_index)->rGetLocation() - rTestPoint.
rGetLocation());
1669 value.node_index = best_node_index;
1670 value.distance = best_node_point_distance;
1672 MPI_Allreduce( &value, &minval, 1, MPI_DOUBLE_INT, MPI_MINLOC, MPI_COMM_WORLD );
1674 return minval.node_index;
1677 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1681 c_vector<double, 2> global_min_max;
1683 MPI_Allreduce(&local_min_max[0], &global_min_max[0], 1, MPI_DOUBLE, MPI_MIN, PETSC_COMM_WORLD);
1684 MPI_Allreduce(&local_min_max[1], &global_min_max[1], 1, MPI_DOUBLE, MPI_MAX, PETSC_COMM_WORLD);
1686 return global_min_max;
1689 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1692 return mHaloNodes.begin();
1695 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1699 for (
unsigned i=0; i<this->mHaloNodes.size(); i++)
1701 c_vector<double, SPACE_DIM>& r_location = this->mHaloNodes[i]->rGetModifiableLocation();
1702 r_location = prod(rotationMatrix, r_location);
1705 for (
unsigned i=0; i<this->mNodes.size(); i++)
1707 c_vector<double, SPACE_DIM>& r_location = this->mNodes[i]->rGetModifiableLocation();
1708 r_location = prod(rotationMatrix, r_location);
1712 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1716 for (
unsigned i=0; i<this->mHaloNodes.size(); i++)
1718 c_vector<double, SPACE_DIM>& r_location = this->mHaloNodes[i]->rGetModifiableLocation();
1719 r_location += rDisplacement;
1722 for (
unsigned i=0; i<this->mNodes.size(); i++)
1724 c_vector<double, SPACE_DIM>& r_location = this->mNodes[i]->rGetModifiableLocation();
1725 r_location += rDisplacement;
1729 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1732 return mHaloNodes.end();
ElementIterator GetElementIteratorBegin()
void SetAttribute(double attribute)
void ConstructLinearMesh(unsigned width)
unsigned GetNumLocalNodes() const
unsigned GetNumBoundaryElements() const
virtual unsigned GetNearestNodeIndex(const ChastePoint< SPACE_DIM > &rTestPoint)
virtual ElementData GetNextElementData()=0
virtual ElementData GetElementData(unsigned index)
c_vector< double, DIM > & rGetLocation()
void RegisterNode(unsigned index)
DistributedTetrahedralMeshPartitionType::type GetPartitionType() const
void Rotate(c_matrix< double, SPACE_DIM, SPACE_DIM > rotationMatrix)
void AddNodeAttribute(double attribute)
NodeIterator GetNodeIteratorEnd()
void SetProcessRegion(ChasteCuboid< SPACE_DIM > *pRegion)
NodeIterator GetNodeIteratorBegin()
#define EXCEPTION(message)
std::vector< Node< SPACE_DIM > * >::const_iterator HaloNodeIterator
virtual void ConstructFromMeshReader(AbstractMeshReader< ELEMENT_DIM, SPACE_DIM > &rMeshReader)
DistributedTetrahedralMesh(DistributedTetrahedralMeshPartitionType::type partitioningMethod=DistributedTetrahedralMeshPartitionType::PARMETIS_LIBRARY)
unsigned GetNumNodes() const
void GetHaloNodeIndices(std::vector< unsigned > &rHaloIndices) const
static void DumbPartitioning(AbstractMesh< ELEMENT_DIM, SPACE_DIM > &rMesh, std::set< unsigned > &rNodesOwned)
static void GeometricPartitioning(AbstractMeshReader< ELEMENT_DIM, SPACE_DIM > &rMeshReader, std::vector< unsigned > &rNodePermutation, std::set< unsigned > &rNodesOwned, std::vector< unsigned > &rProcessorsOffset, ChasteCuboid< SPACE_DIM > *pRegion)
virtual ChasteCuboid< SPACE_DIM > CalculateBoundingBox() const
void SetDistributedVectorFactory(DistributedVectorFactory *pFactory)
bool CalculateDesignatedOwnershipOfBoundaryElement(unsigned faceIndex)
static void PrintAndReset(std::string message)
void ComputeMeshPartitioning(AbstractMeshReader< ELEMENT_DIM, SPACE_DIM > &rMeshReader, std::set< unsigned > &rNodesOwned, std::set< unsigned > &rHaloNodesOwned, std::set< unsigned > &rElementsOwned, std::vector< unsigned > &rProcessorsOffset)
std::string GetOutputDirectoryFullPath() const
void RegisterElement(unsigned index)
HaloNodeIterator GetHaloNodeIteratorBegin() const
virtual ElementData GetNextFaceData()=0
std::vector< unsigned > NodeIndices
virtual bool HasNodePermutation()
HaloNodeIterator GetHaloNodeIteratorEnd() const
void RegisterHaloNode(unsigned index)
unsigned GetNumElements() const
const ChastePoint< SPACE_DIM > & rGetUpperCorner() const
ElementIterator GetElementIteratorEnd()
static void PetscMatrixPartitioning(AbstractMeshReader< ELEMENT_DIM, SPACE_DIM > &rMeshReader, std::vector< unsigned > &rNodePermutation, std::set< unsigned > &rNodesOwned, std::vector< unsigned > &rProcessorsOffset)
virtual unsigned GetNumElements() const =0
unsigned SolveBoundaryElementMapping(unsigned index) const
unsigned SolveNodeMapping(unsigned index) const
unsigned GetNumLocalElements() const
virtual bool HasNclFile()
void Translate(const c_vector< double, SPACE_DIM > &rDisplacement)
virtual void ConstructRectangularMesh(unsigned width, unsigned height, bool stagger=true)
unsigned SolveElementMapping(unsigned index) const
#define EXPORT_TEMPLATE_CLASS_ALL_DIMS(CLASS)
virtual unsigned GetNumFaceAttributes() const
virtual unsigned GetNumElementAttributes() const
virtual c_vector< double, 2 > CalculateMinMaxEdgeLengths()
Node< SPACE_DIM > * GetNodeOrHaloNode(unsigned index) const
virtual ChasteCuboid< SPACE_DIM > CalculateBoundingBox() const
virtual void ConstructCuboid(unsigned width, unsigned height, unsigned depth)
virtual void Scale(const double xFactor=1.0, const double yFactor=1.0, const double zFactor=1.0)
unsigned GetLocalOwnership() const
virtual bool IsFileFormatBinary()
void ConstructCuboid(unsigned width, unsigned height, unsigned depth)
void RegisterBoundaryElement(unsigned index)
unsigned GetNumAllNodes() const
ChasteCuboid< SPACE_DIM > * GetProcessRegion()
DistributedTetrahedralMeshPartitionType::type mPartitioning
virtual unsigned GetNearestNodeIndex(const ChastePoint< SPACE_DIM > &rTestPoint)
unsigned GetNumLocalBoundaryElements() const
virtual std::vector< unsigned > GetContainingElementIndices(unsigned index)
unsigned GetNumHaloNodes() const
virtual void Scale(const double xFactor=1.0, const double yFactor=1.0, const double zFactor=1.0)
void ParMetisLibraryNodeAndElementPartitioning(AbstractMeshReader< ELEMENT_DIM, SPACE_DIM > &rMeshReader, std::set< unsigned > &rElementsOwned, std::set< unsigned > &rNodesOwned, std::set< unsigned > &rHaloNodesOwned, std::vector< unsigned > &rProcessorsOffset)
virtual std::string GetMeshFileBaseName()
virtual void SetDistributedVectorFactory(DistributedVectorFactory *pFactory)
virtual void ConstructLinearMesh(unsigned width)
virtual std::vector< double > GetNextNode()=0
unsigned GetIndex() const
bool CalculateDesignatedOwnershipOfElement(unsigned elementIndex)
virtual void WriteFilesUsingMesh(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > &rMesh, bool keepOriginalElementIndexing=true)
void ConstructRectangularMesh(unsigned width, unsigned height, bool stagger=true)
virtual c_vector< double, 2 > CalculateMinMaxEdgeLengths()
void SetElementOwnerships()
virtual ~DistributedTetrahedralMesh()
virtual unsigned GetNumNodes() const =0
virtual std::vector< double > GetNodeAttributes()
virtual unsigned GetNumFaces() const =0
virtual const std::vector< unsigned > & rGetNodePermutation()
const ChastePoint< SPACE_DIM > & rGetLowerCorner() const