00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include "AbstractTetrahedralMesh.hpp"
00030
00032
00034
00035
00036 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00037 void AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::SetElementOwnerships()
00038 {
00039 unsigned lo=this->GetDistributedVectorFactory()->GetLow();
00040 unsigned hi=this->GetDistributedVectorFactory()->GetHigh();
00041 for (unsigned element_index=0; element_index<mElements.size(); element_index++)
00042 {
00043 Element<ELEMENT_DIM, SPACE_DIM>* p_element = mElements[element_index];
00044 p_element->SetOwnership(false);
00045 for (unsigned local_node_index=0; local_node_index< p_element->GetNumNodes(); local_node_index++)
00046 {
00047 unsigned global_node_index = p_element->GetNodeGlobalIndex(local_node_index);
00048 if (lo<=global_node_index && global_node_index<hi)
00049 {
00050 p_element->SetOwnership(true);
00051 break;
00052 }
00053 }
00054 }
00055 }
00056
00057 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00058 AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::AbstractTetrahedralMesh()
00059 : mMeshIsLinear(true)
00060 {
00061 }
00062
00063 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00064 AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::~AbstractTetrahedralMesh()
00065 {
00066
00067 for (unsigned i=0; i<mElements.size(); i++)
00068 {
00069 delete mElements[i];
00070 }
00071
00072 for (unsigned i=0; i<mBoundaryElements.size(); i++)
00073 {
00074 delete mBoundaryElements[i];
00075 }
00076 }
00077
00078 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00079 unsigned AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::GetNumElements() const
00080 {
00081 return mElements.size();
00082 }
00083
00084 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00085 unsigned AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::GetNumLocalElements() const
00086 {
00087 return GetNumElements();
00088 }
00089
00090 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00091 unsigned AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::GetNumAllElements() const
00092 {
00093 return mElements.size();
00094 }
00095
00096 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00097 unsigned AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::GetNumAllBoundaryElements() const
00098 {
00099 return mBoundaryElements.size();
00100 }
00101
00102 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00103 unsigned AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::GetNumBoundaryElements() const
00104 {
00105 return mBoundaryElements.size();
00106 }
00107
00108
00109 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00110 unsigned AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::GetNumCableElements() const
00111 {
00112 return 0u;
00113 }
00114
00115 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00116 Element<ELEMENT_DIM, SPACE_DIM>* AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::GetElement(unsigned index) const
00117 {
00118 unsigned local_index = SolveElementMapping(index);
00119 return mElements[local_index];
00120 }
00121
00122 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00123 BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>* AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::GetBoundaryElement(unsigned index) const
00124 {
00125 unsigned local_index = SolveBoundaryElementMapping(index);
00126 return mBoundaryElements[local_index];
00127 }
00128
00129 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00130 typename AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::BoundaryElementIterator AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::GetBoundaryElementIteratorBegin() const
00131 {
00132 return mBoundaryElements.begin();
00133 }
00134
00135 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00136 typename AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::BoundaryElementIterator AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::GetBoundaryElementIteratorEnd() const
00137 {
00138 return mBoundaryElements.end();
00139 }
00140
00141 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00142 void AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::GetInverseJacobianForElement(
00143 unsigned elementIndex,
00144 c_matrix<double, SPACE_DIM, ELEMENT_DIM>& rJacobian,
00145 double& rJacobianDeterminant,
00146 c_matrix<double, ELEMENT_DIM, SPACE_DIM>& rInverseJacobian) const
00147 {
00148 mElements[SolveElementMapping(elementIndex)]->CalculateInverseJacobian(rJacobian, rJacobianDeterminant, rInverseJacobian);
00149 }
00150
00151 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00152 void AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::GetWeightedDirectionForBoundaryElement(
00153 unsigned elementIndex,
00154 c_vector<double, SPACE_DIM>& rWeightedDirection,
00155 double& rJacobianDeterminant) const
00156 {
00157 mBoundaryElements[SolveBoundaryElementMapping(elementIndex)]->CalculateWeightedDirection(rWeightedDirection, rJacobianDeterminant );
00158 }
00159
00160
00161 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00162 void AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ConstructLinearMesh(unsigned width)
00163 {
00164 assert(ELEMENT_DIM == 1);
00165
00166 for (unsigned node_index=0; node_index<=width; node_index++)
00167 {
00168 Node<SPACE_DIM>* p_node = new Node<SPACE_DIM>(node_index, node_index==0 || node_index==width, node_index);
00169 this->mNodes.push_back(p_node);
00170 if (node_index==0)
00171 {
00172 this->mBoundaryNodes.push_back(p_node);
00173 this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(0, p_node) );
00174 }
00175 if (node_index==width)
00176 {
00177 this->mBoundaryNodes.push_back(p_node);
00178 this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(1, p_node) );
00179 }
00180 if (node_index>0)
00181 {
00182 std::vector<Node<SPACE_DIM>*> nodes;
00183 nodes.push_back(this->mNodes[node_index-1]);
00184 nodes.push_back(this->mNodes[node_index]);
00185 this->mElements.push_back(new Element<ELEMENT_DIM,SPACE_DIM>(node_index-1, nodes) );
00186 }
00187 }
00188
00189 this->RefreshMesh();
00190 }
00191
00192
00193 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00194 void AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ConstructRectangularMesh(unsigned width, unsigned height, bool stagger)
00195 {
00196 assert(SPACE_DIM == 2);
00197 assert(ELEMENT_DIM == 2);
00198
00199
00200 unsigned node_index=0;
00201 for (unsigned j=0; j<height+1; j++)
00202 {
00203 for (unsigned i=0; i<width+1; i++)
00204 {
00205 bool is_boundary=false;
00206 if (i==0 || j==0 || i==width || j==height)
00207 {
00208 is_boundary=true;
00209 }
00210
00211 assert(node_index==(width+1)*(j) + i);
00212 Node<SPACE_DIM>* p_node = new Node<SPACE_DIM>(node_index++, is_boundary, i, j);
00213 this->mNodes.push_back(p_node);
00214 if (is_boundary)
00215 {
00216 this->mBoundaryNodes.push_back(p_node);
00217 }
00218 }
00219 }
00220
00221
00222 unsigned belem_index=0;
00223
00224 for (unsigned i=0; i<width; i++)
00225 {
00226 std::vector<Node<SPACE_DIM>*> nodes;
00227 nodes.push_back(this->mNodes[height*(width+1)+i]);
00228 nodes.push_back(this->mNodes[height*(width+1)+i+1]);
00229 assert(belem_index==i);
00230 this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,nodes));
00231 }
00232
00233 for (unsigned j=1; j<=height; j++)
00234 {
00235 std::vector<Node<SPACE_DIM>*> nodes;
00236 nodes.push_back(this->mNodes[(width+1)*(j+1)-1]);
00237 nodes.push_back(this->mNodes[(width+1)*j-1]);
00238 assert(belem_index==width+j-1);
00239 this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,nodes));
00240 }
00241
00242 for (unsigned i=0; i<width; i++)
00243 {
00244 std::vector<Node<SPACE_DIM>*> nodes;
00245 nodes.push_back(this->mNodes[i+1]);
00246 nodes.push_back(this->mNodes[i]);
00247 assert(belem_index==width+height+i);
00248 this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,nodes));
00249 }
00250
00251 for (unsigned j=0; j<height; j++)
00252 {
00253 std::vector<Node<SPACE_DIM>*> nodes;
00254 nodes.push_back(this->mNodes[(width+1)*(j+1)]);
00255 nodes.push_back(this->mNodes[(width+1)*(j)]);
00256 assert(belem_index==2*width+height+j);
00257 this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,nodes));
00258 }
00259
00260
00261 unsigned elem_index = 0;
00262 for (unsigned j=0; j<height; j++)
00263 {
00264 for (unsigned i=0; i<width; i++)
00265 {
00266 unsigned parity=(i+(height-j))%2;
00267 unsigned nw=(j+1)*(width+1)+i;
00268 unsigned sw=(j)*(width+1)+i;
00269 std::vector<Node<SPACE_DIM>*> upper_nodes;
00270 upper_nodes.push_back(this->mNodes[nw]);
00271 upper_nodes.push_back(this->mNodes[nw+1]);
00272 if (stagger==false || parity == 1)
00273 {
00274 upper_nodes.push_back(this->mNodes[sw+1]);
00275 }
00276 else
00277 {
00278 upper_nodes.push_back(this->mNodes[sw]);
00279 }
00280 assert(elem_index==2*(j*width+i));
00281 this->mElements.push_back(new Element<ELEMENT_DIM,SPACE_DIM>(elem_index++,upper_nodes));
00282 std::vector<Node<SPACE_DIM>*> lower_nodes;
00283 lower_nodes.push_back(this->mNodes[sw+1]);
00284 lower_nodes.push_back(this->mNodes[sw]);
00285 if (stagger==false ||parity == 1)
00286 {
00287 lower_nodes.push_back(this->mNodes[nw]);
00288 }
00289 else
00290 {
00291 lower_nodes.push_back(this->mNodes[nw+1]);
00292 }
00293 this->mElements.push_back(new Element<ELEMENT_DIM,SPACE_DIM>(elem_index++,lower_nodes));
00294 }
00295 }
00296
00297 this->RefreshMesh();
00298 }
00299
00300
00301
00302
00303 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00304 void AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ConstructCuboid(unsigned width,
00305 unsigned height,
00306 unsigned depth)
00307 {
00308 assert(SPACE_DIM == 3);
00309 assert(ELEMENT_DIM == 3);
00310
00311
00312 unsigned node_index = 0;
00313 for (unsigned k=0; k<depth+1; k++)
00314 {
00315 for (unsigned j=0; j<height+1; j++)
00316 {
00317 for (unsigned i=0; i<width+1; i++)
00318 {
00319 bool is_boundary = false;
00320 if (i==0 || j==0 || k==0 || i==width || j==height || k==depth)
00321 {
00322 is_boundary = true;
00323 }
00324 assert(node_index == (k*(height+1)+j)*(width+1)+i);
00325 Node<SPACE_DIM>* p_node = new Node<SPACE_DIM>(node_index++, is_boundary, i, j, k);
00326
00327 this->mNodes.push_back(p_node);
00328 if (is_boundary)
00329 {
00330 this->mBoundaryNodes.push_back(p_node);
00331 }
00332 }
00333 }
00334 }
00335
00336
00337
00338 unsigned elem_index = 0;
00339 unsigned belem_index = 0;
00340 unsigned element_nodes[6][4] = {{0, 1, 5, 7}, {0, 1, 3, 7},
00341 {0, 2, 3, 7}, {0, 2, 6, 7},
00342 {0, 4, 6, 7}, {0, 4, 5, 7}};
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355 std::vector<Node<SPACE_DIM>*> tetrahedra_nodes;
00356
00357 for (unsigned k=0; k<depth; k++)
00358 {
00359 if (k!=0)
00360 {
00361
00362 assert(belem_index == 2*(height*width+k*2*(height+width)) );
00363 }
00364 for (unsigned j=0; j<height; j++)
00365 {
00366 for (unsigned i=0; i<width; i++)
00367 {
00368
00369 unsigned global_node_indices[8];
00370 unsigned local_node_index = 0;
00371
00372 for (unsigned z = 0; z < 2; z++)
00373 {
00374 for (unsigned y = 0; y < 2; y++)
00375 {
00376 for (unsigned x = 0; x < 2; x++)
00377 {
00378 global_node_indices[local_node_index] = i+x+(width+1)*(j+y+(height+1)*(k+z));
00379
00380 local_node_index++;
00381 }
00382 }
00383 }
00384
00385 for (unsigned m = 0; m < 6; m++)
00386 {
00387
00388
00389 tetrahedra_nodes.clear();
00390
00391 for (unsigned n = 0; n < 4; n++)
00392 {
00393 tetrahedra_nodes.push_back(this->mNodes[global_node_indices[element_nodes[m][n]]]);
00394 }
00395
00396 assert(elem_index == 6 * ((k*height+j)*width+i)+m );
00397 this->mElements.push_back(new Element<ELEMENT_DIM,SPACE_DIM>(elem_index++, tetrahedra_nodes));
00398 }
00399
00400
00401 std::vector<Node<SPACE_DIM>*> triangle_nodes;
00402
00403 if (i == 0)
00404 {
00405 triangle_nodes.clear();
00406 triangle_nodes.push_back(this->mNodes[global_node_indices[0]]);
00407 triangle_nodes.push_back(this->mNodes[global_node_indices[2]]);
00408 triangle_nodes.push_back(this->mNodes[global_node_indices[6]]);
00409 this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,triangle_nodes));
00410 triangle_nodes.clear();
00411 triangle_nodes.push_back(this->mNodes[global_node_indices[0]]);
00412 triangle_nodes.push_back(this->mNodes[global_node_indices[6]]);
00413 triangle_nodes.push_back(this->mNodes[global_node_indices[4]]);
00414 this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,triangle_nodes));
00415 }
00416 if (i == width-1)
00417 {
00418 triangle_nodes.clear();
00419 triangle_nodes.push_back(this->mNodes[global_node_indices[1]]);
00420 triangle_nodes.push_back(this->mNodes[global_node_indices[5]]);
00421 triangle_nodes.push_back(this->mNodes[global_node_indices[7]]);
00422 this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,triangle_nodes));
00423 triangle_nodes.clear();
00424 triangle_nodes.push_back(this->mNodes[global_node_indices[1]]);
00425 triangle_nodes.push_back(this->mNodes[global_node_indices[7]]);
00426 triangle_nodes.push_back(this->mNodes[global_node_indices[3]]);
00427 this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,triangle_nodes));
00428 }
00429 if (j == 0)
00430 {
00431 triangle_nodes.clear();
00432 triangle_nodes.push_back(this->mNodes[global_node_indices[0]]);
00433 triangle_nodes.push_back(this->mNodes[global_node_indices[5]]);
00434 triangle_nodes.push_back(this->mNodes[global_node_indices[1]]);
00435 this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,triangle_nodes));
00436 triangle_nodes.clear();
00437 triangle_nodes.push_back(this->mNodes[global_node_indices[0]]);
00438 triangle_nodes.push_back(this->mNodes[global_node_indices[4]]);
00439 triangle_nodes.push_back(this->mNodes[global_node_indices[5]]);
00440 this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,triangle_nodes));
00441 }
00442 if (j == height-1)
00443 {
00444 triangle_nodes.clear();
00445 triangle_nodes.push_back(this->mNodes[global_node_indices[2]]);
00446 triangle_nodes.push_back(this->mNodes[global_node_indices[3]]);
00447 triangle_nodes.push_back(this->mNodes[global_node_indices[7]]);
00448 this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,triangle_nodes));
00449 triangle_nodes.clear();
00450 triangle_nodes.push_back(this->mNodes[global_node_indices[2]]);
00451 triangle_nodes.push_back(this->mNodes[global_node_indices[7]]);
00452 triangle_nodes.push_back(this->mNodes[global_node_indices[6]]);
00453 this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,triangle_nodes));
00454 }
00455 if (k == 0)
00456 {
00457 triangle_nodes.clear();
00458 triangle_nodes.push_back(this->mNodes[global_node_indices[0]]);
00459 triangle_nodes.push_back(this->mNodes[global_node_indices[3]]);
00460 triangle_nodes.push_back(this->mNodes[global_node_indices[2]]);
00461 this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,triangle_nodes));
00462 triangle_nodes.clear();
00463 triangle_nodes.push_back(this->mNodes[global_node_indices[0]]);
00464 triangle_nodes.push_back(this->mNodes[global_node_indices[1]]);
00465 triangle_nodes.push_back(this->mNodes[global_node_indices[3]]);
00466 this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,triangle_nodes));
00467 }
00468 if (k == depth-1)
00469 {
00470 triangle_nodes.clear();
00471 triangle_nodes.push_back(this->mNodes[global_node_indices[4]]);
00472 triangle_nodes.push_back(this->mNodes[global_node_indices[7]]);
00473 triangle_nodes.push_back(this->mNodes[global_node_indices[5]]);
00474 this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,triangle_nodes));
00475 triangle_nodes.clear();
00476 triangle_nodes.push_back(this->mNodes[global_node_indices[4]]);
00477 triangle_nodes.push_back(this->mNodes[global_node_indices[6]]);
00478 triangle_nodes.push_back(this->mNodes[global_node_indices[7]]);
00479 this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,triangle_nodes));
00480 }
00481 }
00482 }
00483 }
00484
00485 this->RefreshMesh();
00486 }
00487
00488 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00489 void AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ConstructRegularSlabMesh(double spaceStep, double width, double height, double depth)
00490 {
00491 assert(spaceStep>0.0);
00492 assert(width>0.0);
00493 if (ELEMENT_DIM > 1)
00494 {
00495 assert(height>0.0);
00496 }
00497 if (ELEMENT_DIM > 2)
00498 {
00499 assert(depth>0.0);
00500 }
00501 unsigned num_elem_x=(unsigned)((width+0.5*spaceStep)/spaceStep);
00502 unsigned num_elem_y=(unsigned)((height+0.5*spaceStep)/spaceStep);
00503 unsigned num_elem_z=(unsigned)((depth+0.5*spaceStep)/spaceStep);
00504
00505
00506 {
00507 double actual_width_x=num_elem_x*spaceStep;
00508 double actual_width_y=num_elem_y*spaceStep;
00509 double actual_width_z=num_elem_z*spaceStep;
00510
00511
00512
00513
00514 if ( fabs (actual_width_x - width) > DBL_EPSILON*width
00515 ||( height!= 0.0 && fabs (actual_width_y - height) > DBL_EPSILON*height)
00516 ||( depth != 0.0 && fabs (actual_width_z - depth) > DBL_EPSILON*depth ) )
00517 {
00518 EXCEPTION("Space step does not divide the size of the mesh");
00519 }
00520 }
00521 switch (ELEMENT_DIM)
00522 {
00523 case 1:
00524 this->ConstructLinearMesh(num_elem_x);
00525 break;
00526 case 2:
00527 this->ConstructRectangularMesh(num_elem_x, num_elem_y, true);
00528 break;
00529 default:
00530 case 3:
00531 this->ConstructCuboid(num_elem_x, num_elem_y, num_elem_z);
00532 }
00533 this->Scale(spaceStep, spaceStep, spaceStep);
00534 }
00535
00536 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00537 bool AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::CalculateDesignatedOwnershipOfBoundaryElement( unsigned faceIndex )
00538 {
00539
00540 unsigned tie_break_index = this->GetBoundaryElement(faceIndex)->GetNodeGlobalIndex(0);
00541
00542
00543 if (this->GetDistributedVectorFactory()->IsGlobalIndexLocal(tie_break_index))
00544 {
00545 return true;
00546 }
00547 else
00548 {
00549 return false;
00550 }
00551 }
00552
00553 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00554 bool AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::CalculateDesignatedOwnershipOfElement( unsigned elementIndex )
00555 {
00556
00557 unsigned tie_break_index = this->GetElement(elementIndex)->GetNodeGlobalIndex(0);
00558
00559
00560 if (this->GetDistributedVectorFactory()->IsGlobalIndexLocal(tie_break_index))
00561 {
00562 return true;
00563 }
00564 else
00565 {
00566 return false;
00567 }
00568 }
00569
00570 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00571 unsigned AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::CalculateMaximumNodeConnectivityPerProcess() const
00572 {
00573 unsigned max_num=0u;
00574 unsigned connected_node_index=0u;
00575 for (unsigned local_node_index=0; local_node_index<this->mNodes.size(); local_node_index++)
00576 {
00577 unsigned num=this->mNodes[local_node_index]->GetNumContainingElements();
00578 if (num>max_num)
00579 {
00580 max_num=num;
00581 connected_node_index=local_node_index;
00582 }
00583 }
00584 if (max_num == 0u)
00585 {
00586 #define COVERAGE_IGNORE
00587
00588
00589
00590
00591
00592 assert(this->mNodes.size() == 0u);
00593 return(1u);
00594
00595
00596
00597 #undef COVERAGE_IGNORE
00598 }
00599
00600 std::set<unsigned> forward_star_nodes;
00601 unsigned nodes_per_element = this->mElements[0]->GetNumNodes();
00602 for (typename Node<SPACE_DIM>::ContainingElementIterator it = this->mNodes[connected_node_index]->ContainingElementsBegin();
00603 it != this->mNodes[connected_node_index]->ContainingElementsEnd();
00604 ++it)
00605 {
00606 Element<ELEMENT_DIM, SPACE_DIM>* p_elem=this->GetElement(*it);
00607 for (unsigned i=0; i<nodes_per_element; i++)
00608 {
00609 forward_star_nodes.insert(p_elem->GetNodeGlobalIndex(i));
00610 }
00611 }
00612 return forward_star_nodes.size();
00613 }
00614
00615 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00616 void AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::GetHaloNodeIndices(std::vector<unsigned>& rHaloIndices) const
00617 {
00618
00619 rHaloIndices.clear();
00620 }
00621
00622 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00623 void AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ConstructFromMesh(AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>& rOtherMesh)
00624 {
00625 for (unsigned i=0; i<rOtherMesh.GetNumNodes(); i++)
00626 {
00627 Node<SPACE_DIM>* p_node = rOtherMesh.GetNode(i);
00628 assert(!p_node->IsDeleted());
00629 c_vector<double, SPACE_DIM> location=p_node->rGetLocation();
00630 bool is_boundary=p_node->IsBoundaryNode();
00631
00632 Node<SPACE_DIM>* p_node_copy = new Node<SPACE_DIM>(i, location, is_boundary);
00633 this->mNodes.push_back( p_node_copy );
00634 if (is_boundary)
00635 {
00636 this->mBoundaryNodes.push_back( p_node_copy );
00637 }
00638 }
00639
00640 for (unsigned i=0; i<rOtherMesh.GetNumElements(); i++)
00641 {
00642 Element<ELEMENT_DIM, SPACE_DIM>* p_elem = rOtherMesh.GetElement(i);
00643 assert(!p_elem->IsDeleted());
00644 std::vector<Node<SPACE_DIM>*> nodes_for_element;
00645 for (unsigned j=0; j<p_elem->GetNumNodes(); j++)
00646 {
00647 nodes_for_element.push_back(this->mNodes[ p_elem->GetNodeGlobalIndex(j) ]);
00648 }
00649 Element<ELEMENT_DIM, SPACE_DIM>* p_elem_copy=new Element<ELEMENT_DIM, SPACE_DIM>(i, nodes_for_element);
00650 p_elem_copy->RegisterWithNodes();
00651 this->mElements.push_back(p_elem_copy);
00652 }
00653
00654 for (unsigned i=0; i<rOtherMesh.GetNumBoundaryElements(); i++)
00655 {
00656 BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>* p_b_elem = rOtherMesh.GetBoundaryElement(i);
00657 assert(!p_b_elem->IsDeleted());
00658 std::vector<Node<SPACE_DIM>*> nodes_for_element;
00659 for (unsigned j=0; j<p_b_elem->GetNumNodes(); j++)
00660 {
00661 nodes_for_element.push_back(this->mNodes[ p_b_elem->GetNodeGlobalIndex(j) ]);
00662 }
00663 BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>* p_b_elem_copy=new BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>(i, nodes_for_element);
00664 p_b_elem_copy->RegisterWithNodes();
00665 this->mBoundaryElements.push_back(p_b_elem_copy);
00666 }
00667 this->RefreshMesh();
00668
00669 }
00670
00671 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00672 void AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::CalculateNodeExchange(
00673 std::vector<std::vector<unsigned> >& rNodesToSendPerProcess,
00674 std::vector<std::vector<unsigned> >& rNodesToReceivePerProcess)
00675 {
00676 assert( rNodesToSendPerProcess.empty() );
00677 assert( rNodesToReceivePerProcess.empty() );
00678
00679
00680 std::vector<std::set<unsigned> > node_sets_to_send_per_process;
00681 std::vector<std::set<unsigned> > node_sets_to_receive_per_process;
00682
00683 node_sets_to_send_per_process.resize(PetscTools::GetNumProcs());
00684 node_sets_to_receive_per_process.resize(PetscTools::GetNumProcs());
00685 std::vector<unsigned> global_lows = this->GetDistributedVectorFactory()->rGetGlobalLows();
00686
00687 for (typename AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ElementIterator iter = this->GetElementIteratorBegin();
00688 iter != this->GetElementIteratorEnd();
00689 ++iter)
00690 {
00691 std::vector <unsigned> nodes_on_this_process;
00692 std::vector <unsigned> nodes_not_on_this_process;
00693
00694 for (unsigned i=0; i<ELEMENT_DIM+1; i++)
00695 {
00696 unsigned node_index=iter->GetNodeGlobalIndex(i);
00697 if (this->GetDistributedVectorFactory()->IsGlobalIndexLocal(node_index))
00698 {
00699 nodes_on_this_process.push_back(node_index);
00700 }
00701 else
00702 {
00703 nodes_not_on_this_process.push_back(node_index);
00704 }
00705 }
00706
00707
00708
00709
00710 if (nodes_on_this_process.empty())
00711 {
00712 continue;
00713 }
00714
00715 if(!nodes_not_on_this_process.empty())
00716 {
00717 for (unsigned i=0; i<nodes_not_on_this_process.size(); i++)
00718 {
00719
00720 unsigned remote_process=global_lows.size()-1;
00721 for(; global_lows[remote_process] > nodes_not_on_this_process[i]; remote_process--)
00722 {
00723 }
00724
00725
00726 node_sets_to_receive_per_process[remote_process].insert(nodes_not_on_this_process[i]);
00727
00728
00729 for (unsigned j=0; j<nodes_on_this_process.size(); j++)
00730 {
00731 node_sets_to_send_per_process[remote_process].insert(nodes_on_this_process[j]);
00732 }
00733 }
00734 }
00735 }
00736
00737 for (unsigned process_number = 0; process_number < PetscTools::GetNumProcs(); process_number++)
00738 {
00739 std::vector<unsigned> process_send_vector( node_sets_to_send_per_process[process_number].begin(),
00740 node_sets_to_send_per_process[process_number].end() );
00741 std::sort(process_send_vector.begin(), process_send_vector.end());
00742
00743 rNodesToSendPerProcess.push_back(process_send_vector);
00744
00745 std::vector<unsigned> process_receive_vector( node_sets_to_receive_per_process[process_number].begin(),
00746 node_sets_to_receive_per_process[process_number].end() );
00747 std::sort(process_receive_vector.begin(), process_receive_vector.end());
00748
00749 rNodesToReceivePerProcess.push_back(process_receive_vector);
00750 }
00751
00752 }
00753
00754
00756
00758
00759 template class AbstractTetrahedralMesh<1,1>;
00760 template class AbstractTetrahedralMesh<1,2>;
00761 template class AbstractTetrahedralMesh<1,3>;
00762 template class AbstractTetrahedralMesh<2,2>;
00763 template class AbstractTetrahedralMesh<2,3>;
00764 template class AbstractTetrahedralMesh<3,3>;