1 // ATC header files 2 #include "FE_Element.h" 3 #include "FE_Mesh.h" 4 #include "LammpsInterface.h" 5 #include "ATC_Error.h" 6 #include "OutputManager.h" 7 #include "Utility.h" 8 #include <sstream> 9 #include <cassert> 10 #include <algorithm> 11 #include <cmath> 12 #include <functional> 13 14 // Other headers 15 #include <iostream> 16 17 18 using namespace std; 19 using ATC_Utility::dbl_geq; 20 using ATC_Utility::parse_min; 21 using ATC_Utility::parse_max; 22 using ATC_Utility::parse_minmax; 23 using ATC_Utility::split_values; 24 using ATC_Utility::plane_coords; 25 using ATC_Utility::norm3; 26 using ATC_Utility::to_string; 27 using ATC_Utility::tolerance; 28 29 30 31 32 33 namespace ATC { 34 35 // constants 36 const static double tangentTolerance = 0.01; 37 const static double tol = 1.0e-10; 38 39 40 // ============================================================= 41 // class FE_Mesh 42 // ============================================================= FE_Mesh()43 FE_Mesh::FE_Mesh() 44 : decomposition_(false), 45 lammpsPartition_(false), 46 partitioned_(false), 47 nNodes_(0), 48 nNodesUnique_(0), 49 feElement_(nullptr), 50 twoDimensional_(false), 51 hasPlanarFaces_(false) 52 53 { 54 } 55 56 // ------------------------------------------------------------- ~FE_Mesh()57 FE_Mesh::~FE_Mesh() 58 { 59 if (feElement_) delete feElement_; 60 } 61 62 // ------------------------------------------------------------- 63 // modify 64 // ------------------------------------------------------------- modify(int narg,char ** arg)65 bool FE_Mesh::modify(int narg, char **arg) 66 { 67 bool match = false; 68 69 if (strcmp(arg[0],"mesh")==0) 70 { 71 /*! \page man_mesh_create_faceset_box fix_modify AtC mesh create_faceset box 72 \section syntax 73 fix_modify AtC mesh create_faceset <id> box 74 <xmin> <xmax> <ymin> <ymax> <zmin> <zmax> <in|out> [units] 75 - <id> = id to assign to the collection of FE faces 76 - <xmin> <xmax> <ymin> <ymax> <zmin> <zmax> = coordinates of 77 the bounding box that is coincident with the desired FE faces 78 - <in|out> = "in" gives inner faces to the box, 79 "out" gives the outer faces to the box 80 - units = option to specify real as opposed to lattice units 81 \section examples 82 <TT> fix_modify AtC mesh create_faceset obndy box -4.0 4.0 -12 12 -12 12 out </TT> 83 \section description 84 Command to assign an id to a set of FE faces. 85 \section restrictions 86 Only viable for rectangular grids. 87 \section related 88 \section default 89 The default options are units = lattice and the use of outer faces 90 */ 91 /*! \page man_mesh_create_faceset_plane fix_modify AtC mesh create_faceset plane 92 \section syntax 93 fix_modify AtC mesh create_faceset <id> plane 94 <x|y|z> <val1> <x|y|z> <lval2> <uval2> [units] 95 - <id> = id to assign to the collection of FE faces 96 - <x|y|z> = coordinate directions that define plane on which faceset lies 97 - <val1>,<lval2>,<uval2> = plane is specified as the x|y|z=val1 plane bounded by 98 the segments x|y|z = [lval2,uval2] 99 - units = option to specify real as opposed to lattice units 100 \section examples 101 <TT> fix_modify AtC mesh create_faceset xyplane plane y 0 x -4 0 </TT> 102 \section description 103 Command to assign an id to a set of FE faces. 104 \section restrictions 105 Only viable for rectangular grids. 106 \section related 107 \section default 108 The default option is units = lattice. 109 */ 110 if (strcmp(arg[1],"create_faceset")==0) 111 { 112 int argIdx = 2; 113 string tag = arg[argIdx++]; 114 if (strcmp(arg[argIdx],"plane")==0) 115 { 116 argIdx++; 117 int ndir, idir[3], isgn; 118 double xlimits[3][2]; 119 parse_plane(argIdx, narg, arg, ndir, idir, isgn, xlimits); 120 if (xlimits[idir[1]][0] == xlimits[idir[1]][1]) 121 split_values(xlimits[idir[1]][0],xlimits[idir[1]][1]); 122 if (xlimits[idir[2]][0] == xlimits[idir[2]][1]) 123 split_values(xlimits[idir[2]][0],xlimits[idir[2]][1]); 124 parse_units(argIdx,narg,arg, 125 xlimits[0][0],xlimits[0][1], 126 xlimits[1][0],xlimits[1][1], 127 xlimits[2][0],xlimits[2][1]); 128 if (ndir > 1) { 129 create_faceset(tag, xlimits[idir[0]][0], idir[0], isgn, 130 idir[1], xlimits[idir[1]][0], xlimits[idir[1]][1], 131 idir[2], xlimits[idir[2]][0], xlimits[idir[2]][1]); 132 } 133 else { 134 create_faceset(tag, xlimits[idir[0]][0], idir[0], isgn); 135 } 136 match = true; 137 } 138 // bounding_box 139 else 140 { 141 if (strcmp(arg[argIdx],"box")==0) argIdx++; 142 double xmin = parse_min(arg[argIdx++]); 143 double xmax = parse_max(arg[argIdx++]); 144 double ymin = parse_min(arg[argIdx++]); 145 double ymax = parse_max(arg[argIdx++]); 146 double zmin = parse_min(arg[argIdx++]); 147 double zmax = parse_max(arg[argIdx++]); 148 bool outward = true; 149 if (narg > argIdx && (strcmp(arg[argIdx++],"in") == 0)) 150 outward = false; 151 parse_units(argIdx,narg,arg, xmin,xmax,ymin,ymax,zmin,zmax); 152 create_faceset(tag, xmin, xmax, ymin, ymax, zmin, zmax, outward); 153 match = true; 154 } 155 } 156 /*! \page man_mesh_create_nodeset fix_modify AtC mesh create_nodeset 157 \section syntax 158 fix_modify AtC mesh create_nodeset <id> 159 <xmin> <xmax> <ymin> <ymax> <zmin> <zmax> 160 - <id> = id to assign to the collection of FE nodes 161 - <xmin> <xmax> <ymin> <ymax> <zmin> <zmax> = coordinates of 162 the bounding box that contains only the desired nodes 163 \section examples 164 <TT> fix_modify AtC mesh create_nodeset lbc -12.1 -11.9 -12 12 -12 12 </TT> 165 \section description 166 Command to assign an id to a set of FE nodes to be used subsequently 167 in defining boundary conditions. 168 \section restrictions 169 \section related 170 \section default 171 Coordinates are assumed to be in lattice units. 172 */ 173 else if (strcmp(arg[1],"create_nodeset")==0) { 174 int argIdx = 2; 175 string tag = arg[argIdx++]; 176 double xmin, xmax, ymin, ymax, zmin, zmax; 177 if (strcmp(arg[argIdx],"plane")==0) { 178 argIdx++; 179 int ndir, idir[3], isgn; 180 double xlimits[3][2]; 181 parse_plane(argIdx, narg, arg, ndir, idir, isgn, xlimits); 182 xmin = xlimits[0][0]; 183 xmax = xlimits[0][1]; 184 ymin = xlimits[1][0]; 185 ymax = xlimits[1][1]; 186 zmin = xlimits[2][0]; 187 zmax = xlimits[2][1]; 188 } 189 else { 190 xmin = parse_min(arg[argIdx++]); 191 xmax = parse_max(arg[argIdx++]); 192 ymin = parse_min(arg[argIdx++]); 193 ymax = parse_max(arg[argIdx++]); 194 zmin = parse_min(arg[argIdx++]); 195 zmax = parse_max(arg[argIdx++]); 196 } 197 if (xmin == xmax) split_values(xmin,xmax); 198 if (ymin == ymax) split_values(ymin,ymax); 199 if (zmin == zmax) split_values(zmin,zmax); 200 parse_units(argIdx,narg,arg, xmin,xmax,ymin,ymax,zmin,zmax); 201 create_nodeset(tag, xmin, xmax, ymin, ymax, zmin, zmax); 202 match = true; 203 } 204 /*! \page man_mesh_add_to_nodeset fix_modify AtC mesh add_to_nodeset 205 \section syntax 206 fix_modify AtC mesh add_to_nodeset <id> 207 <xmin> <xmax> <ymin> <ymax> <zmin> <zmax> 208 - <id> = id of FE nodeset to be added to 209 - <xmin> <xmax> <ymin> <ymax> <zmin> <zmax> = coordinates of 210 the bounding box that contains the desired nodes to be added 211 \section examples 212 <TT> fix_modify AtC mesh add_to_nodeset lbc -11.9 -11 -12 12 -12 12 </TT> 213 \section description 214 Command to add nodes to an already existing FE nodeset. 215 \section restrictions 216 \section related 217 \section default 218 Coordinates are assumed to be in lattice units. 219 */ 220 else if (strcmp(arg[1],"add_to_nodeset")==0) { 221 string tag = arg[2]; 222 double xmin = parse_min(arg[3]); 223 double xmax = parse_max(arg[4]); 224 double ymin = parse_min(arg[5]); 225 double ymax = parse_max(arg[6]); 226 double zmin = parse_min(arg[7]); 227 double zmax = parse_max(arg[8]); 228 add_to_nodeset(tag, xmin, xmax, ymin, ymax, zmin, zmax); 229 match = true; 230 } 231 /*! \page man_mesh_create_elementset fix_modify AtC mesh create_elementset 232 \section syntax 233 fix_modify AtC create_elementset <id> 234 <xmin> <xmax> <ymin> <ymax> <zmin> <zmax> 235 - <id> = id to assign to the collection of FE element 236 - <xmin> <xmax> <ymin> <ymax> <zmin> <zmax> = coordinates of 237 the bounding box that contains only the desired elements 238 \section examples 239 <TT> fix_modify AtC mesh create_elementset middle -4.1 4.1 -100 100 -100 1100 </TT> 240 \section description 241 Command to assign an id to a set of FE elements to be used subsequently 242 in defining material and mesh-based operations. 243 \section restrictions 244 Only viable for rectangular grids. 245 \section related 246 \section default 247 Coordinates are assumed to be in lattice units. 248 */ 249 else if (strcmp(arg[1],"create_elementset")==0) { 250 int argIdx = 2; 251 string tag = arg[argIdx++]; 252 double xmin = 0; 253 double xmax = 0; 254 double ymin = 0; 255 double ymax = 0; 256 double zmin = 0; 257 double zmax = 0; 258 if (narg > 4) { 259 xmin = parse_min(arg[argIdx++]); 260 xmax = parse_max(arg[argIdx++]); 261 ymin = parse_min(arg[argIdx++]); 262 ymax = parse_max(arg[argIdx++]); 263 zmin = parse_min(arg[argIdx++]); 264 zmax = parse_max(arg[argIdx++]); 265 if (xmin == xmax) split_values(xmin,xmax); 266 if (ymin == ymax) split_values(ymin,ymax); 267 if (zmin == zmax) split_values(zmin,zmax); 268 parse_units(argIdx,narg,arg, xmin,xmax,ymin,ymax,zmin,zmax); 269 } 270 else { 271 string regionName = arg[argIdx++]; 272 LammpsInterface::instance()-> 273 region_bounds(regionName.c_str(),xmin,xmax,ymin,ymax,zmin,zmax); 274 } 275 create_elementset(tag, xmin, xmax, ymin, ymax, zmin, zmax); 276 match = true; 277 } 278 /*! \page man_mesh_nodeset_to_elementset fix_modify AtC mesh nodeset_to_elementset 279 \section syntax 280 fix_modify AtC nodeset_to_elementset <nodeset_id> <elementset_id> <max/min> 281 - <nodeset_id> = id of desired nodeset from which to create elementset 282 - <elementset_id> = id to assign to the collection of FE element 283 - <max/min> = flag to choose either the maximal or minimal elementset 284 \section examples 285 <TT> fix_modify AtC mesh nodeset_to_elementset myNodeset myElementset min </TT> 286 \section description 287 Command to create an elementset from an existing nodeset. Either the minimal element set 288 of elements with all nodes in the set, or maximal element set with all elements with at 289 least one node in the set, can be created 290 \section restrictions 291 None. 292 \section related 293 \section default 294 Unless specified, the maximal element set is created 295 */ 296 else if (strcmp(arg[1],"nodeset_to_elementset")==0) { 297 string nodesetId = arg[2]; 298 string elementsetId = arg[3]; 299 set<int> elementSet; 300 if (narg < 5) { 301 this->nodeset_to_maximal_elementset(nodesetId,elementSet); 302 } 303 else { 304 if (strcmp(arg[4],"max")==0) { 305 this->nodeset_to_maximal_elementset(nodesetId,elementSet); 306 } 307 else if (strcmp(arg[4],"min")==0) { 308 this->nodeset_to_minimal_elementset(nodesetId,elementSet); 309 } 310 else { 311 return false; 312 } 313 } 314 elementSetMap_[elementsetId] = elementSet; 315 match = true; 316 } 317 /*! \page man_mesh_output fix_modify AtC mesh output 318 \section syntax 319 fix_modify AtC mesh output <file_prefix> 320 \section examples 321 <TT> fix_modify AtC mesh output meshData </TT> \n 322 \section description 323 Command to output mesh and associated data: nodesets, facesets, and 324 elementsets. This data is only output once upon initialization since 325 currently the mesh is static. Creates (binary, "gold" format) Ensight 326 output of mesh data. 327 \section restrictions 328 none 329 \section related 330 \section default 331 none 332 */ 333 else if (strcmp(arg[1],"output")==0) { 334 string outputPrefix = arg[2]; 335 output(outputPrefix); 336 match = true; 337 } 338 } 339 return match; 340 } 341 // ------------------------------------------------------------- parse_units(int & argIdx,int narg,char ** arg,double & xmin,double & xmax,double & ymin,double & ymax,double & zmin,double & zmax)342 void FE_Mesh::parse_units(int & argIdx, int narg, char ** arg, 343 double & xmin, double & xmax, double & ymin, double & ymax, 344 double & zmin, double & zmax) 345 { 346 if (narg > argIdx && (strcmp(arg[argIdx++],"units") == 0)) {} 347 else 348 { // scale from lattice units to physical units 349 xmin *= xscale_; xmax *= xscale_; 350 ymin *= yscale_; ymax *= yscale_; 351 zmin *= zscale_; zmax *= zscale_; 352 } 353 } 354 // ------------------------------------------------------------- 355 // parse plane 356 // ------------------------------------------------------------- parse_plane(int & argIdx,int narg,char ** arg,int & ndir,int * idir,int & isgn,double xlimits[][2])357 void FE_Mesh::parse_plane(int & argIdx, int narg, char ** arg, 358 int & ndir, int * idir, int & isgn, double xlimits[][2]) 359 { 360 ndir = 0; 361 xlimits[0][0] = parse_min("-INF"); 362 xlimits[0][1] = parse_max("INF"); 363 xlimits[1][0] = parse_min("-INF"); 364 xlimits[1][1] = parse_max("INF"); 365 xlimits[2][0] = parse_min("-INF"); 366 xlimits[2][1] = parse_max("INF"); 367 string n(arg[argIdx++]); 368 int i1dir = -1, i2dir = -1, i3dir = -1; 369 string_to_index(n, i1dir, isgn); 370 idir[ndir++] = i1dir; 371 idir[ndir] = (i1dir+1) % 3; 372 idir[ndir+1] = (i1dir+2) % 3; 373 xlimits[i1dir][0] = parse_minmax(arg[argIdx++]); 374 xlimits[i1dir][1] = xlimits[i1dir][0]; 375 if (narg > argIdx ) { 376 if (string_to_index(arg[argIdx],i2dir)) { 377 argIdx++; 378 xlimits[i2dir][0] = parse_min(arg[argIdx++]); 379 xlimits[i2dir][1] = parse_max(arg[argIdx++]); 380 idir[ndir++] = i2dir; 381 } 382 if (narg > argIdx ) { 383 if (string_to_index(arg[argIdx],i3dir)) { 384 argIdx++; 385 xlimits[i3dir][0] = parse_min(arg[argIdx++]); 386 xlimits[i3dir][1] = parse_max(arg[argIdx++]); 387 } 388 } 389 else { 390 i3dir = 0; 391 if ((i1dir == 0 && i2dir == 1) || (i1dir == 1 && i2dir == 0) ) { 392 i3dir = 2; 393 } 394 else if ((i1dir == 0 && i2dir == 2) || (i1dir == 2 && i2dir == 0) ) { 395 i3dir = 1; 396 } 397 } 398 idir[ndir++] = i3dir; 399 } 400 if ((idir[0]==idir[1]) || (idir[0]==idir[2]) || (idir[1]==idir[2]) ) { 401 throw ATC_Error( "inconsistent directions in plane:"+to_string(idir[0]+1)+" "+to_string(idir[1]+1)+" "+to_string(idir[2]+1)); 402 } 403 } 404 // ------------------------------------------------------------- 405 // initialize 406 // ------------------------------------------------------------- initialize(void)407 void FE_Mesh::initialize(void) 408 { 409 410 bool aligned = is_aligned(); 411 if (!aligned) { 412 feElement_->set_projection_guess(CENTROID_LINEARIZED); 413 ATC::LammpsInterface::instance()->print_msg_once("WARNING: mesh is not aligned with the coordinate directions atom-to-element mapping will be expensive"); 414 // if HEX8 -> orient(); 415 } 416 bool twoD = is_two_dimensional(); 417 if (twoD) { 418 feElement_->set_projection_guess(TWOD_ANALYTIC); 419 if (feElement_->order()< 3) hasPlanarFaces_ = true; 420 ATC::LammpsInterface::instance()->print_msg_once(" mesh is two dimensional"); 421 } 422 } 423 //----------------------------------------------------------------- 424 425 //----------------------------------------------------------------- write_mesh(string meshFile)426 void FE_Mesh::write_mesh(string meshFile) 427 { 428 ofstream out; 429 out.open(meshFile.c_str()); 430 DENS_MAT & x = *(coordinates()); 431 Array2D<int> & conn = *(connectivity()); 432 int nNodes = x.nCols(); // transpose 433 int ndm = x.nRows(); 434 int nElems = conn.nCols(); 435 int nodesPerElem = conn.nRows(); // transpose 436 out << "Coordinates " << nNodes << "\n"; 437 for (int n = 0; n < nNodes; n++) { 438 for (int i = 0; i < ndm; i++) { 439 out << " " << std::setprecision(16) << x(i,n); 440 } 441 out << "\n"; 442 } 443 out << "\n"; 444 string type = element_type(); 445 out << "Elements " << nElems << " " << type << "\n"; 446 for (int n = 0; n < nElems; n++) { 447 for (int i = 0; i < nodesPerElem; i++) { 448 out << 1+conn(i,n) << " "; 449 } 450 out << "\n"; 451 } 452 out << "\n"; 453 if (nodeSetMap_.size()) { 454 out << "Nodesets " << nodeSetMap_.size() <<"\n"; 455 NODE_SET_MAP::const_iterator niter; 456 map<string,DENS_MAT> nodesets; 457 for (niter = nodeSetMap_.begin(); niter != nodeSetMap_.end(); niter++) { 458 string name = niter->first; 459 const set<int> & nset = niter->second; 460 out << name << " " << nset.size() << "\n"; 461 set<int>::const_iterator iter; 462 for (iter = nset.begin(); iter != nset.end(); iter++) { 463 out << *iter << " " ; 464 } 465 out << "\n"; 466 } 467 } 468 } 469 // ------------------------------------------------------------- 470 // test whether almost structured 471 // ------------------------------------------------------------- is_aligned(void) const472 bool FE_Mesh::is_aligned(void) const 473 { 474 vector<bool> foundBestMatch(nSD_,false); 475 vector<DENS_VEC> tangents(nSD_); 476 DENS_VEC xi0(nSD_); 477 xi0 = 0; 478 DENS_MAT eltCoords; 479 for (int ielem = 0; ielem < nElts_; ielem++) { 480 element_coordinates(ielem,eltCoords); 481 feElement_->tangents(eltCoords,xi0,tangents,true); 482 483 for (unsigned i = 0; i < tangents.size(); i++) { 484 // find maximum value for which global axis its closest to 485 int maxIndex = 0; 486 double maxValue = abs(tangents[i](0)); 487 for (int j = 1; j < nSD_; j++) { 488 if (abs(tangents[i](j)) > maxValue) { 489 maxValue = abs(tangents[i](j)); 490 maxIndex = j; 491 } 492 } 493 494 // make sure no other tangent is associated with this direction 495 if (foundBestMatch[maxIndex]) { 496 return false; 497 } 498 else { 499 foundBestMatch[maxIndex] = true; 500 } 501 502 // compute deviation from a perfectly aligned vector 503 double error = 0.; 504 for (int j = 1; j < nSD_; j++) { 505 if (j != maxIndex) { 506 error += abs(tangents[i](j)); 507 } 508 } 509 error /= maxValue; 510 if (error > tangentTolerance) { 511 return false; 512 } 513 } 514 } 515 return true; 516 } 517 518 // ------------------------------------------------------------- 519 // element_type 520 // ------------------------------------------------------------- element_type(void) const521 string FE_Mesh::element_type(void) const { 522 int npe = feElement_->num_elt_nodes(); 523 if (npe == 4) { return "TET4"; } 524 else if (npe == 8) { return "HEX8"; } 525 else if (npe == 20) { return "HEX20"; } 526 else if (npe == 27) { return "HEX27"; } 527 return "UNKNOWN"; 528 } 529 530 // ------------------------------------------------------------- 531 // element_coordinates 532 // ------------------------------------------------------------- element_coordinates(const int eltID,DENS_MAT & xCoords) const533 void FE_Mesh::element_coordinates(const int eltID, 534 DENS_MAT & xCoords) const 535 { 536 const int nne = num_nodes_per_element(); 537 xCoords.reset(nSD_, nne, false); 538 for (int inode=0; inode<nne; inode++) { 539 const int id = element_connectivity_global(eltID, inode); 540 for (int isd=0; isd<nSD_; isd++) { 541 xCoords(isd,inode) = nodalCoords_(isd,id); 542 } 543 } 544 545 } 546 // ------------------------------------------------------------- 547 // position 548 // ------------------------------------------------------------- position(const int eltID,const VECTOR & xi,DENS_VEC & x) const549 void FE_Mesh::position(const int eltID, 550 const VECTOR & xi, 551 DENS_VEC & x) const 552 { 553 const int nne = num_nodes_per_element(); 554 DENS_VEC N; 555 feElement_->shape_function(xi,N); 556 x.reset(nSD_); 557 for (int inode=0; inode<nne; inode++) { 558 const int id = element_connectivity_global(eltID, inode); 559 for (int isd=0; isd<nSD_; isd++) { 560 x(isd) += nodalCoords_(isd,id)*N(inode); 561 } 562 } 563 } 564 565 // ------------------------------------------------------------- 566 // element size in each direction 567 // ------------------------------------------------------------- bounding_box(const int ielem,DENS_VEC & xmin,DENS_VEC & xmax)568 void FE_Mesh::bounding_box(const int ielem, 569 DENS_VEC & xmin, DENS_VEC & xmax) 570 { 571 xmin.reset(nSD_); 572 xmax.reset(nSD_); 573 int nne = num_nodes_per_element(); 574 for (int isd=0; isd<nSD_; isd++) { 575 int id = element_connectivity_global(ielem, 0); 576 double x = nodalCoords_(isd,id); 577 xmin(isd) = x; 578 xmax(isd) = x; 579 for (int inode=1; inode<nne; inode++) { 580 id = element_connectivity_global(ielem, inode); 581 x = nodalCoords_(isd,id); 582 xmin(isd) = min(xmin(isd), x ); 583 xmax(isd) = max(xmax(isd), x ); 584 } 585 } 586 } 587 588 // ------------------------------------------------------------- 589 // element size in each direction 590 // ------------------------------------------------------------- element_size(const int ielem,double & hx,double & hy,double & hz)591 void FE_Mesh::element_size(const int ielem, 592 double & hx, double & hy, double & hz) 593 { 594 DENS_VEC xmin(nSD_), xmax(nSD_); 595 bounding_box(ielem,xmin,xmax); 596 hx = xmax(0)-xmin(0); 597 hy = xmax(1)-xmin(1); 598 hz = xmax(2)-xmin(2); 599 } 600 601 // ------------------------------------------------------------- 602 // face_coordinates 603 // ------------------------------------------------------------- face_coordinates(const PAIR face,DENS_MAT & xCoords) const604 void FE_Mesh::face_coordinates(const PAIR face, DENS_MAT & xCoords) const 605 { 606 const int eltID=face.first, faceID=face.second; 607 const int nnf = num_nodes_per_face(); 608 const Array2D <int> & local_conn = local_face_connectivity(); 609 610 xCoords.reset(nSD_, nnf, false); 611 612 for (int inode=0; inode < nnf; inode++) 613 { 614 int id = element_connectivity_global(eltID, local_conn(faceID,inode)); 615 for (int isd=0; isd<nSD_; isd++) 616 xCoords(isd,inode) = nodalCoords_(isd,id); 617 } 618 } 619 620 // ------------------------------------------------------------- 621 // nodal_coordinates 622 // ------------------------------------------------------------- nodal_coordinates(const int nodeID) const623 DENS_VEC FE_Mesh::nodal_coordinates(const int nodeID) const 624 { 625 DENS_VEC xCoords(nSD_, false); 626 const int id = uniqueToGlobalMap_(nodeID); 627 for (int isd=0; isd<nSD_; isd++) 628 xCoords(isd) = nodalCoords_(isd, id); 629 return xCoords; 630 } global_coordinates(const int nodeID) const631 DENS_VEC FE_Mesh::global_coordinates(const int nodeID) const 632 { 633 DENS_VEC xCoords(nSD_, false); 634 for (int isd=0; isd<nSD_; isd++) 635 xCoords(isd) = nodalCoords_(isd, nodeID); 636 return xCoords; 637 } 638 639 // ------------------------------------------------------------- 640 // query_nodeset 641 // ------------------------------------------------------------- query_nodeset(const string & name) const642 bool FE_Mesh::query_nodeset(const string & name) const 643 { 644 if (name == "all") return true; 645 if (nodeSetMap_.find(name) == nodeSetMap_.end()) return false; 646 return true; 647 } 648 649 // ------------------------------------------------------------- 650 // get_nodeset 651 // ------------------------------------------------------------- nodeset(const string & name) const652 const set<int> & FE_Mesh::nodeset(const string & name) const 653 { 654 NODE_SET_MAP::const_iterator iter = nodeSetMap_.find(name); 655 if (name == "all") return nodeSetAll_; 656 else if (iter == nodeSetMap_.end()) 657 throw ATC_Error( "No nodeset with name " + name + " found."); 658 else return iter->second; 659 } 660 661 // ------------------------------------------------------------- 662 // get_elementset 663 // ------------------------------------------------------------- elementset(const string & name) const664 const set<int> & FE_Mesh::elementset(const string & name) const 665 { 666 NODE_SET_MAP::const_iterator iter = elementSetMap_.find(name); 667 if (name == "all") return elementSetAll_; 668 else if (iter == elementSetMap_.end()) 669 throw ATC_Error( "No elementset with name " + name + " found."); 670 else return iter->second; 671 } 672 673 // ------------------------------------------------------------- 674 // nodeset_to_minimal_elementset 675 // ------------------------------------------------------------- nodeset_to_minimal_elementset(const string & name,set<int> & elemSet) const676 void FE_Mesh::nodeset_to_minimal_elementset 677 (const string & name, set<int> & elemSet) const 678 { 679 if (name == "all") { 680 for (int ielem = 0; ielem < nElts_; ielem++) { 681 elemSet.insert(ielem); 682 } 683 } 684 else { 685 NODE_SET_MAP::const_iterator iter = nodeSetMap_.find(name); 686 if (iter == nodeSetMap_.end()) 687 throw ATC_Error( "No nodeset with name " + name + " found."); 688 nodeset_to_minimal_elementset(iter->second,elemSet); 689 if (elemSet.size()==0) { 690 throw ATC_Error("No elements found in minimal condensation of nodeset " + name); 691 } 692 } 693 } 694 695 // ------------------------------------------------------------- 696 // nodeset_to_minimal_elementset 697 // ------------------------------------------------------------- nodeset_to_minimal_elementset(const set<int> & nodeSet,set<int> & elemSet) const698 void FE_Mesh::nodeset_to_minimal_elementset 699 (const set<int> & nodeSet, set<int> & elemSet) const 700 { 701 int npe = num_nodes_per_element(); 702 for (int ielem=0; ielem < nElts_; ielem++) { 703 int inode = 0; 704 bool in = true; 705 while (in && inode < npe) { 706 int node = element_connectivity_unique(ielem, inode); 707 set<int>::const_iterator iter = nodeSet.find(node); 708 if (iter == nodeSet.end()) { in=false; } 709 inode++; 710 } 711 if (in) elemSet.insert(ielem); 712 } 713 } 714 715 // ------------------------------------------------------------- 716 // nodeset_to_maximal_elementset 717 // ------------------------------------------------------------- nodeset_to_maximal_elementset(const string & name,set<int> & elemSet) const718 void FE_Mesh::nodeset_to_maximal_elementset(const string &name, set<int> &elemSet) const 719 { 720 if (name == "all") { 721 for (int ielem = 0; ielem < nElts_; ielem++) { 722 elemSet.insert(ielem); 723 } 724 } 725 else { 726 NODE_SET_MAP::const_iterator iter = nodeSetMap_.find(name); 727 if (iter == nodeSetMap_.end()) 728 throw ATC_Error( "No nodeset with name " + name + " found."); 729 nodeset_to_maximal_elementset(iter->second,elemSet); 730 if (elemSet.size()==0) { 731 throw ATC_Error("No elements found in maximal condensation of nodeset " + name); 732 } 733 } 734 } 735 736 // ------------------------------------------------------------- 737 // nodeset_to_maximal_elementset 738 // ------------------------------------------------------------- nodeset_to_maximal_elementset(const set<int> & nodeSet,set<int> & elemSet) const739 void FE_Mesh::nodeset_to_maximal_elementset(const set<int> &nodeSet, set<int> &elemSet) const 740 { 741 int npe = num_nodes_per_element(); 742 for (int ielem = 0; ielem < nElts_; ielem++) { 743 int inode = 0; 744 bool in = false; 745 while (!in && inode < npe) { 746 int node = element_connectivity_unique(ielem, inode); 747 set<int>::const_iterator iter = nodeSet.find(node); 748 if (iter != nodeSet.end()) { in = true; } 749 inode++; 750 } 751 if (in) elemSet.insert(ielem); 752 } 753 } 754 755 // ------------------------------------------------------------- 756 // elementset_to_nodeset 757 // ------------------------------------------------------------- elementset_to_nodeset(const set<int> & elemSet,set<int> nodeSet) const758 void FE_Mesh::elementset_to_nodeset 759 (const set<int> & elemSet, set<int> nodeSet) const 760 { 761 int npe = num_nodes_per_element(); 762 set<int>::const_iterator itr; 763 for (itr = elemSet.begin(); itr != elemSet.end(); itr++) 764 { 765 int ielem = *itr; 766 for (int inode=0; inode < npe; inode++) 767 { 768 int node = element_connectivity_global(ielem, inode); 769 nodeSet.insert(node); 770 } 771 } 772 } 773 774 // ------------------------------------------------------------- 775 // elementset_to_nodeset 776 // ------------------------------------------------------------- elementset_to_nodeset(const string & name,set<int> nodeSet) const777 void FE_Mesh::elementset_to_nodeset 778 (const string & name, set<int> nodeSet) const 779 { 780 if (name == "all") 781 for (int ielem = 0; ielem < nElts_; ielem++) 782 nodeSet.insert(ielem); 783 784 else 785 { 786 ELEMENT_SET_MAP::const_iterator iter = elementSetMap_.find(name); 787 if (iter == elementSetMap_.end()) 788 throw ATC_Error( "No elementset with name " + name + " found."); 789 790 int npe = num_nodes_per_element(); 791 const set<int> &elemSet = iter->second; 792 set<int>::const_iterator itr; 793 for (itr = elemSet.begin(); itr != elemSet.end(); itr++) 794 { 795 int ielem = *itr; 796 for (int inode=0; inode < npe; inode++) 797 { 798 int node = element_connectivity_unique(ielem, inode); 799 nodeSet.insert(node); 800 inode++; // XXX: is this correct? 801 } 802 } 803 } 804 } elementset_to_nodeset(const string & name) const805 set<int> FE_Mesh::elementset_to_nodeset 806 (const string & name) const 807 { 808 set<int> nset; 809 elementset_to_nodeset(name,nset); 810 return nset; 811 } 812 813 // ------------------------------------------------------------- 814 // elementset_to_minimal_nodeset 815 // ------------------------------------------------------------- elementset_to_minimal_nodeset(const string & name,set<int> & nodeSet) const816 void FE_Mesh::elementset_to_minimal_nodeset 817 (const string & name, set<int> & nodeSet) const 818 { 819 // return: set - complement_of_set 820 if (name == "all") { return;} 821 else 822 { 823 elementset_to_nodeset(name,nodeSet); 824 set<int> compElemSet; 825 elementset_complement(name,compElemSet); 826 int npe = num_nodes_per_element(); 827 set<int>::const_iterator itr; 828 for (itr = compElemSet.begin(); itr != compElemSet.end(); itr++) 829 { 830 int ielem = *itr; 831 for (int inode=0; inode < npe; inode++) 832 { 833 int node = element_connectivity_unique(ielem, inode); 834 nodeSet.erase(node); 835 inode++; // XXX: is this correct? 836 } 837 } 838 } 839 } 840 841 // ------------------------------------------------------------- 842 // elementset_complement 843 // ------------------------------------------------------------- elementset_complement(const string & name,set<int> & cElemSet) const844 void FE_Mesh::elementset_complement 845 (const string & name, set<int> & cElemSet) const 846 { 847 // return: set - complement_of_set 848 if (name == "all") { return;} 849 else 850 { 851 ELEMENT_SET_MAP::const_iterator iter = elementSetMap_.find(name); 852 if (iter == elementSetMap_.end()) 853 throw ATC_Error( "No elementset with name " + name + " found."); 854 855 const set<int> &elemSet = iter->second; 856 for (int ielem = 0; ielem < nElts_; ielem++) 857 { 858 if(elemSet.find(ielem) == elemSet.end() ) cElemSet.insert(ielem); 859 } 860 } 861 } 862 863 // ------------------------------------------------------------- 864 // elementset_complement 865 // ------------------------------------------------------------- elementset_complement(const set<int> & elemSet,set<int> & cElemSet) const866 void FE_Mesh::elementset_complement 867 (const set<int> & elemSet, set<int> & cElemSet) const 868 { 869 for (int ielem = 0; ielem < nElts_; ielem++) 870 { 871 if(elemSet.find(ielem) == elemSet.end() ) cElemSet.insert(ielem); 872 } 873 } 874 // ------------------------------------------------------------- 875 // faceset_to_nodeset 876 // ------------------------------------------------------------- faceset_to_nodeset(const string & name,set<int> & nodeSet) const877 void FE_Mesh::faceset_to_nodeset(const string &name, set<int> &nodeSet) const 878 { 879 if (name == "all") { 880 for (int inode = 0; inode < nNodesUnique_; inode++) 881 nodeSet.insert(inode); 882 } 883 else 884 { 885 FACE_SET_MAP::const_iterator faceset = faceSetMap_.find(name); 886 if (faceset == faceSetMap_.end()) 887 throw ATC_Error( "No faceset with name " + name + " found."); 888 const set<PAIR> & faceSet = faceset->second; 889 set<PAIR>::const_iterator iter; 890 Array <int> conn; 891 for (iter = faceSet.begin(); iter != faceSet.end(); iter++) 892 { 893 PAIR face = *iter; 894 face_connectivity_unique(face,conn); 895 for (int i = 0; i < conn.size() ; ++i) { 896 int inode = conn(i); 897 nodeSet.insert(inode); 898 } 899 } 900 } 901 } faceset_to_nodeset(const set<PAIR> & faceSet,set<int> & nodeSet) const902 void FE_Mesh::faceset_to_nodeset(const set<PAIR> &faceSet, set<int> &nodeSet) const 903 { 904 set<PAIR>::const_iterator iter; 905 Array <int> conn; 906 for (iter = faceSet.begin(); iter != faceSet.end(); iter++) 907 { 908 PAIR face = *iter; 909 face_connectivity_unique(face,conn); 910 for (int i = 0; i < conn.size() ; ++i) { 911 int inode = conn(i); 912 nodeSet.insert(inode); 913 } 914 } 915 } 916 // ------------------------------------------------------------- 917 // faceset_to_nodeset 918 // ------------------------------------------------------------- faceset_to_nodeset_global(const string & name,set<int> & nodeSet) const919 void FE_Mesh::faceset_to_nodeset_global(const string &name, set<int> &nodeSet) const 920 { 921 if (name == "all") { 922 for (int inode = 0; inode < nNodes_; inode++) 923 nodeSet.insert(inode); 924 } 925 else 926 { 927 FACE_SET_MAP::const_iterator faceset = faceSetMap_.find(name); 928 if (faceset == faceSetMap_.end()) 929 throw ATC_Error( "No faceset with name " + name + " found."); 930 const set<PAIR> & faceSet = faceset->second; 931 set<PAIR>::const_iterator iter; 932 Array <int> conn; 933 for (iter = faceSet.begin(); iter != faceSet.end(); iter++) 934 { 935 PAIR face = *iter; 936 face_connectivity(face,conn); 937 for (int i = 0; i < conn.size() ; ++i) { 938 int inode = conn(i); 939 nodeSet.insert(inode); 940 } 941 } 942 } 943 } faceset_to_nodeset_global(const set<PAIR> & faceSet,set<int> & nodeSet) const944 void FE_Mesh::faceset_to_nodeset_global(const set<PAIR> &faceSet, set<int> &nodeSet) const 945 { 946 set<PAIR>::const_iterator iter; 947 Array <int> conn; 948 for (iter = faceSet.begin(); iter != faceSet.end(); iter++) 949 { 950 PAIR face = *iter; 951 face_connectivity(face,conn); 952 for (int i = 0; i < conn.size() ; ++i) { 953 int inode = conn(i); 954 nodeSet.insert(inode); 955 } 956 } 957 } 958 959 // ------------------------------------------------------------- 960 // get_faceset 961 // ------------------------------------------------------------- faceset(const string & name) const962 const set<PAIR> &FE_Mesh::faceset(const string & name) const 963 { 964 FACE_SET_MAP::const_iterator iter = faceSetMap_.find(name); 965 if (iter == faceSetMap_.end()) 966 { 967 throw ATC_Error( "No faceset with name " + name + " found."); 968 } 969 return iter->second; 970 } 971 972 // ------------------------------------------------------------- 973 // create_nodeset 974 // ------------------------------------------------------------- create_nodeset(const string & name,const set<int> & nodeSet)975 void FE_Mesh::create_nodeset(const string & name, 976 const set<int> & nodeSet) 977 { 978 // Make sure we don't already have a nodeset with this name 979 NODE_SET_MAP::iterator iter = nodeSetMap_.find(name); 980 if (iter != nodeSetMap_.end()) { 981 string message("A nodeset with name " + name + " is already defined."); 982 throw ATC_Error( message); 983 } 984 nodeSetMap_[name] = nodeSet; 985 986 if (ATC::LammpsInterface::instance()->rank_zero()) { 987 stringstream ss; 988 ss << "created nodeset " << name 989 << " with " << nodeSet.size() << " nodes"; 990 ATC::LammpsInterface::instance()->print_msg_once(ss.str()); 991 } 992 } create_nodeset(const string & name,double xmin,double xmax,double ymin,double ymax,double zmin,double zmax)993 void FE_Mesh::create_nodeset(const string & name, 994 double xmin, 995 double xmax, 996 double ymin, 997 double ymax, 998 double zmin, 999 double zmax) 1000 { 1001 // Make sure we don't already have a nodeset with this name 1002 NODE_SET_MAP::iterator iter = nodeSetMap_.find(name); 1003 if (iter != nodeSetMap_.end()) { 1004 string message("A nodeset with name " + name + " is already defined."); 1005 throw ATC_Error( message); 1006 } 1007 1008 set<int> nodeSet; 1009 1010 // Loop over nodes and add their unique id's to the set if they're 1011 // in the correct range 1012 for (int inode = 0; inode < nNodes_; inode++) { 1013 double x = nodalCoords_(0,inode); 1014 double y = nodalCoords_(1,inode); 1015 double z = nodalCoords_(2,inode); 1016 if ( (xmin <= x) && (x <= xmax) && 1017 (ymin <= y) && (y <= ymax) && 1018 (zmin <= z) && (z <= zmax) ) { 1019 int uid = globalToUniqueMap_(inode); 1020 nodeSet.insert(uid); 1021 } 1022 } 1023 if (nodeSet.size() == 0) { 1024 string message("nodeset " + name + " has zero size."); 1025 throw ATC_Error( message); 1026 } 1027 1028 nodeSetMap_[name] = nodeSet; 1029 1030 if (ATC::LammpsInterface::instance()->rank_zero()) { 1031 stringstream ss; 1032 ss << "created nodeset " << name 1033 << " with " << nodeSet.size() << " nodes"; 1034 ATC::LammpsInterface::instance()->print_msg_once(ss.str()); 1035 } 1036 } 1037 1038 // ------------------------------------------------------------- 1039 // add_to_nodeset 1040 // ------------------------------------------------------------- add_to_nodeset(const string & name,double xmin,double xmax,double ymin,double ymax,double zmin,double zmax)1041 void FE_Mesh::add_to_nodeset(const string & name, 1042 double xmin, 1043 double xmax, 1044 double ymin, 1045 double ymax, 1046 double zmin, 1047 double zmax) 1048 { 1049 // Make sure we already have a nodeset with this name 1050 NODE_SET_MAP::iterator iter = nodeSetMap_.find(name); 1051 if (iter == nodeSetMap_.end()) { 1052 string message("A nodeset with name " +name + " is not already defined."); 1053 throw ATC_Error( message); 1054 } 1055 1056 set<int> nodeSet; 1057 1058 // Loop over nodes and add their unique id's to the set if they're 1059 // in the correct range 1060 for (int inode = 0; inode < nNodes_; inode++) { 1061 double x = nodalCoords_(0,inode); 1062 double y = nodalCoords_(1,inode); 1063 double z = nodalCoords_(2,inode); 1064 if ( (xmin <= x) && (x <= xmax) && 1065 (ymin <= y) && (y <= ymax) && 1066 (zmin <= z) && (z <= zmax) ) { 1067 int uid = globalToUniqueMap_(inode); 1068 nodeSet.insert(uid); 1069 } 1070 } 1071 if (nodeSet.size() == 0) { 1072 string message("nodeset " + name + " has zero size."); 1073 throw ATC_Error( message); 1074 } 1075 1076 nodeSetMap_[name].insert(nodeSet.begin(),nodeSet.end()); 1077 1078 if (ATC::LammpsInterface::instance()->rank_zero()) { 1079 stringstream ss; 1080 ss << "added " << nodeSet.size() << " nodes to nodeset " << name ; 1081 ATC::LammpsInterface::instance()->print_msg_once(ss.str()); 1082 } 1083 } 1084 1085 // ------------------------------------------------------------- 1086 // create_faceset 1087 // ------------------------------------------------------------- create_faceset(const string & name,double xmin,double xmax,double ymin,double ymax,double zmin,double zmax,bool outward)1088 void FE_Mesh::create_faceset(const string & name, 1089 double xmin, 1090 double xmax, 1091 double ymin, 1092 double ymax, 1093 double zmin, 1094 double zmax, 1095 bool outward) 1096 { 1097 // Make sure we don't already have a nodeset with this name 1098 FACE_SET_MAP::iterator iter = faceSetMap_.find(name); 1099 if (iter != faceSetMap_.end()) 1100 throw ATC_Error( "A faceset with name " + name + " is already defined."); 1101 1102 set<PAIR> faceSet; 1103 // Loop over face and add their unique id's to the set if they concide 1104 // with region 1105 const int nf = num_faces_per_element(); 1106 const int npf = num_nodes_per_face(); 1107 const Array2D<int> & face_conn = local_face_connectivity(); 1108 for (int ielem = 0; ielem < nElts_; ielem++) 1109 { 1110 for (int iface = 0; iface < nf; iface++) 1111 { 1112 bool in = true; 1113 bool on_xmin = true, on_xmax = true; 1114 bool on_ymin = true, on_ymax = true; 1115 bool on_zmin = true, on_zmax = true; 1116 bool x_neg = false, x_pos = false; 1117 bool y_neg = false, y_pos = false; 1118 bool z_neg = false, z_pos = false; 1119 double x,y,z; 1120 for (int inode = 0; inode < npf; inode++) 1121 { 1122 x = nodalCoords_(0,connectivity_(face_conn(iface,inode),ielem)); 1123 y = nodalCoords_(1,connectivity_(face_conn(iface,inode),ielem)); 1124 z = nodalCoords_(2,connectivity_(face_conn(iface,inode),ielem)); 1125 1126 if ( x + tol < xmin) { in = false; break; } 1127 if ( x - tol > xmax) { in = false; break; } 1128 if ( y + tol < ymin) { in = false; break; } 1129 if ( y - tol > ymax) { in = false; break; } 1130 if ( z + tol < zmin) { in = false; break; } 1131 if ( z - tol > zmax) { in = false; break; } 1132 1133 on_xmin = on_xmin && fabs(x-xmin) <= tol; 1134 on_xmax = on_xmax && fabs(x-xmax) <= tol; 1135 on_ymin = on_ymin && fabs(y-ymin) <= tol; 1136 on_ymax = on_ymax && fabs(y-ymax) <= tol; 1137 on_zmin = on_zmin && fabs(z-zmin) <= tol; 1138 on_zmax = on_zmax && fabs(z-zmax) <= tol; 1139 } 1140 if (in) { 1141 // note based on structured grid 1142 if (outward) 1143 { 1144 if (on_xmin && iface==0) { x_neg = true;} 1145 if (on_xmax && iface==1) { x_pos = true;} 1146 if (on_ymin && iface==2) { y_neg = true;} 1147 if (on_ymax && iface==3) { y_pos = true;} 1148 if (on_zmin && iface==4) { z_neg = true;} 1149 if (on_zmax && iface==5) { z_pos = true;} 1150 } 1151 else 1152 { 1153 if (on_xmin && iface==1) { x_pos = true;} 1154 if (on_xmax && iface==0) { x_neg = true;} 1155 if (on_ymin && iface==3) { y_pos = true;} 1156 if (on_ymax && iface==2) { y_neg = true;} 1157 if (on_zmin && iface==5) { z_pos = true;} 1158 if (on_zmax && iface==4) { z_neg = true;} 1159 } 1160 1161 if ( (x_neg || x_pos) || (y_neg || y_pos) || (z_neg || z_pos) ) { 1162 PAIR face(ielem,iface); 1163 faceSet.insert(face); 1164 } 1165 } 1166 } 1167 } 1168 if (faceSet.empty()) throw ATC_Error( "faceset "+name+" is empty."); 1169 1170 faceSetMap_[name] = faceSet; 1171 if (ATC::LammpsInterface::instance()->comm_rank() == 0) { 1172 stringstream ss; 1173 ss << "created faceset " << name 1174 << " with " << faceSet.size() << " faces"; 1175 ATC::LammpsInterface::instance()->print_msg_once(ss.str()); 1176 } 1177 } 1178 create_faceset(const string & name,double xRef,int nIdx,int nSgn,int nIdx2,double x2lo,double x2hi,int nIdx3,double x3lo,double x3hi)1179 void FE_Mesh::create_faceset(const string & name, 1180 double xRef, 1181 int nIdx, int nSgn, 1182 int nIdx2, double x2lo, double x2hi, 1183 int nIdx3, double x3lo, double x3hi) 1184 { 1185 double xtol = tolerance(xRef); 1186 // Make sure we don't already have a faceset with this name 1187 FACE_SET_MAP::iterator iter = faceSetMap_.find(name); 1188 if (iter != faceSetMap_.end()) 1189 throw ATC_Error( "A faceset with name "+name+" is already defined."); 1190 1191 bool finite2 = (nIdx2 >= 0); 1192 bool finite3 = (nIdx3 >= 0); 1193 1194 set<PAIR> faceSet; 1195 // Loop over faces i& add unique id's to the set if concide w/ plane 1196 int nf = num_faces_per_element(); 1197 int npf = num_nodes_per_face(); 1198 const Array2D<int> & face_conn = local_face_connectivity(); 1199 for (int ielem = 0; ielem < nElts_; ielem++) 1200 { 1201 for (int iface = 0; iface < nf; iface++) 1202 { 1203 bool in = true; 1204 // all nodes must be on the plane 1205 for (int inode = 0; inode < npf; inode++) { 1206 int node = connectivity_(face_conn(iface,inode),ielem); 1207 double x = nodalCoords_(nIdx,node); 1208 if ( fabs(x-xRef) > xtol){ in = false; break;} 1209 if (finite2) { 1210 double y = nodalCoords_(nIdx2,node); 1211 if ( y < x2lo || y > x2hi){ in = false; break;} 1212 } 1213 if (finite3) { 1214 double y = nodalCoords_(nIdx3,node); 1215 if ( y < x3lo || y > x3hi){ in = false; break;} 1216 } 1217 } 1218 // check correct orientation 1219 if (in) 1220 { 1221 if ( (nIdx == 0 && iface==0 && nSgn == -1) 1222 || (nIdx == 0 && iface==1 && nSgn == 1) 1223 || (nIdx == 1 && iface==2 && nSgn == -1) 1224 || (nIdx == 1 && iface==3 && nSgn == 1) 1225 || (nIdx == 2 && iface==4 && nSgn == -1) 1226 || (nIdx == 3 && iface==5 && nSgn == 1) ) 1227 { 1228 PAIR face(ielem,iface); 1229 faceSet.insert(face); 1230 } 1231 } 1232 } 1233 } 1234 1235 if (faceSet.empty()) 1236 throw ATC_Error( "faceset "+name+" is empty."); 1237 1238 faceSetMap_[name] = faceSet; 1239 if (ATC::LammpsInterface::instance()->comm_rank() == 0) { 1240 stringstream ss; 1241 ss << "created faceset " << name 1242 << " with " << faceSet.size() << " faces"; 1243 ATC::LammpsInterface::instance()->print_msg_once(ss.str()); 1244 } 1245 } 1246 1247 // ------------------------------------------------------------- 1248 // create_elementset 1249 // ------------------------------------------------------------- create_elementset(const string & name,double xmin,double xmax,double ymin,double ymax,double zmin,double zmax)1250 void FE_Mesh::create_elementset(const string & name, 1251 double xmin, 1252 double xmax, 1253 double ymin, 1254 double ymax, 1255 double zmin, 1256 double zmax) 1257 { 1258 // Make sure we don't already have a elementset with this name 1259 ELEMENT_SET_MAP::iterator iter = elementSetMap_.find(name); 1260 if (iter != elementSetMap_.end()) { 1261 string message("An elementset with name "+name+" is already defined."); 1262 throw ATC_Error( message); 1263 } 1264 1265 set<int> nodeSet; 1266 1267 // Loop over nodes and add their unique id's to the set if they're 1268 // in the correct range 1269 for (int inode = 0; inode < nNodes_; inode++) { 1270 double x = nodalCoords_(0,inode); 1271 double y = nodalCoords_(1,inode); 1272 double z = nodalCoords_(2,inode); 1273 if ( (xmin <= x) && (x <= xmax) && 1274 (ymin <= y) && (y <= ymax) && 1275 (zmin <= z) && (z <= zmax) ) { 1276 int uid = globalToUniqueMap_(inode); 1277 nodeSet.insert(uid); 1278 } 1279 } 1280 if (nodeSet.size() == 0) { 1281 string message("elementset " + name + " has zero size."); 1282 throw ATC_Error( message); 1283 } 1284 1285 // create a minimal element set from all the nodes included in the region 1286 set<int> elemSet; 1287 int npe = num_nodes_per_element(); 1288 for (int ielem=0; ielem < nElts_; ielem++) 1289 { 1290 int inode = 0; 1291 bool in = true; 1292 while (in && inode < npe) 1293 { 1294 int node = connectivityUnique_(inode, ielem); 1295 set<int>::const_iterator iter = nodeSet.find(node); 1296 if (iter == nodeSet.end()) { in=false; } 1297 inode++; 1298 } 1299 if (in) elemSet.insert(ielem); 1300 } 1301 if (elemSet.size() == 0) { 1302 string message("element set " + name + " has zero size."); 1303 throw ATC_Error( message); 1304 } 1305 elementSetMap_[name] = elemSet; 1306 1307 if (ATC::LammpsInterface::instance()->comm_rank() == 0) { 1308 stringstream ss; 1309 ss << "created elementset " << name 1310 << " with " << elemSet.size() << " elements"; 1311 ATC::LammpsInterface::instance()->print_msg_once(ss.str()); 1312 } 1313 } 1314 // ------------------------------------------------------------- 1315 // get_numIPsPerElement() 1316 // ------------------------------------------------------------- num_ips_per_element() const1317 int FE_Mesh::num_ips_per_element() const 1318 { 1319 return feElement_->num_ips(); 1320 } 1321 // ------------------------------------------------------------- 1322 // get_numNodesPerElement() 1323 // ------------------------------------------------------------- num_nodes_per_element() const1324 int FE_Mesh::num_nodes_per_element() const 1325 { 1326 return feElement_->num_elt_nodes(); 1327 } 1328 // ------------------------------------------------------------- 1329 // get_numFacesPerElement() 1330 // ------------------------------------------------------------- num_faces_per_element() const1331 int FE_Mesh::num_faces_per_element() const 1332 { 1333 return feElement_->num_faces(); 1334 } 1335 // ------------------------------------------------------------- 1336 // get_num_ips_per_face() 1337 // ------------------------------------------------------------- num_ips_per_face() const1338 int FE_Mesh::num_ips_per_face() const 1339 { 1340 return feElement_->num_face_ips(); 1341 } 1342 // ------------------------------------------------------------- 1343 // get_num_nodes_per_face() 1344 // ------------------------------------------------------------- num_nodes_per_face() const1345 int FE_Mesh::num_nodes_per_face() const 1346 { 1347 return feElement_->num_face_nodes(); 1348 } 1349 // ------------------------------------------------------------- 1350 // mappings from element id to associated nodes 1351 // ------------------------------------------------------------- 1352 1353 element_connectivity_global(const int eltID,Array<int> & nodes) const1354 void FE_Mesh::element_connectivity_global(const int eltID, 1355 Array<int> & nodes) const 1356 { 1357 const int npe = num_nodes_per_element(); 1358 nodes.reset(npe); 1359 1360 // use connectivity arrays 1361 if (decomposition_ && partitioned_) { 1362 for (int inode = 0; inode < npe; inode++) { 1363 nodes(inode) = myConnectivity_(inode, map_elem_to_myElem(eltID)); 1364 } 1365 } else { 1366 for (int inode = 0; inode < npe; inode++) { 1367 nodes(inode) = connectivity_(inode, eltID); 1368 } 1369 } 1370 } 1371 // ------------------------------------------------------------- 1372 // 1373 // ------------------------------------------------------------- element_connectivity_unique(const int eltID,Array<int> & nodes) const1374 void FE_Mesh::element_connectivity_unique(const int eltID, 1375 Array<int> & nodes) const 1376 { 1377 const int npe = num_nodes_per_element(); 1378 nodes.reset(npe); 1379 1380 // use connectivity arrays 1381 if (decomposition_ && partitioned_) { 1382 for (int inode = 0; inode < npe; inode++) { 1383 nodes(inode) = myConnectivityUnique_(inode, map_elem_to_myElem(eltID)); 1384 } 1385 } else { 1386 for (int inode = 0; inode < npe; inode++) { 1387 nodes(inode) = connectivityUnique_(inode, eltID); 1388 } 1389 } 1390 } 1391 1392 // ------------------------------------------------------------- 1393 // 1394 // ------------------------------------------------------------- element_connectivity_global(const int eltID,const int inode) const1395 int FE_Mesh::element_connectivity_global(const int eltID, 1396 const int inode) const 1397 { 1398 if (decomposition_ && partitioned_) { 1399 return myConnectivity_(inode, map_elem_to_myElem(eltID)); 1400 } else { 1401 return connectivity_(inode, eltID); 1402 } 1403 } 1404 // ------------------------------------------------------------- 1405 // 1406 // ------------------------------------------------------------- element_connectivity_unique(const int eltID,const int inode) const1407 int FE_Mesh::element_connectivity_unique(const int eltID, 1408 const int inode) const 1409 { 1410 if (decomposition_ && partitioned_) { 1411 return myConnectivityUnique_(inode, map_elem_to_myElem(eltID)); 1412 } else { 1413 return connectivityUnique_(inode, eltID); 1414 } 1415 } 1416 // ------------------------------------------------------------- 1417 // 1418 // ------------------------------------------------------------- element_connectivity_global(const int eltID) const1419 AliasArray<int> FE_Mesh::element_connectivity_global(const int eltID) const 1420 { 1421 if (decomposition_ && partitioned_) { 1422 return myConnectivity_.column(map_elem_to_myElem(eltID)); 1423 } else { 1424 return connectivity_.column(eltID); 1425 } 1426 } 1427 // ------------------------------------------------------------- 1428 // 1429 // ------------------------------------------------------------- element_connectivity_unique(const int eltID) const1430 AliasArray<int> FE_Mesh::element_connectivity_unique(const int eltID) const 1431 { 1432 if (decomposition_ && partitioned_) { 1433 return myConnectivityUnique_.column(map_elem_to_myElem(eltID)); 1434 } else { 1435 return connectivityUnique_.column(eltID); 1436 } 1437 } 1438 // ------------------------------------------------------------- 1439 // local_face_connectivity() 1440 // ------------------------------------------------------------- local_face_connectivity() const1441 const Array2D<int> &FE_Mesh::local_face_connectivity() const 1442 { 1443 return feElement_->local_face_conn(); 1444 } 1445 1446 // ------------------------------------------------------------- 1447 // maps to/from partitioned element data 1448 // ------------------------------------------------------------- 1449 map_elem_to_myElem(int elemID) const1450 int FE_Mesh::map_elem_to_myElem(int elemID) const 1451 { 1452 1453 return elemToMyElemMap_.find(elemID)->second; 1454 } 1455 map_myElem_to_elem(int myElemID) const1456 int FE_Mesh::map_myElem_to_elem(int myElemID) const 1457 { 1458 return myElts_[myElemID]; 1459 } 1460 1461 // ------------------------------------------------------------- 1462 // shape function evaluation 1463 // ------------------------------------------------------------- 1464 1465 // set quadrature scheme pass-through set_quadrature(FeIntQuadrature type)1466 void FE_Mesh::set_quadrature(FeIntQuadrature type) 1467 { 1468 feElement_->set_quadrature(type); 1469 } 1470 1471 // shape function evaluation shape_functions(const VECTOR & x,DENS_VEC & N,Array<int> & nodeList) const1472 void FE_Mesh::shape_functions(const VECTOR &x, 1473 DENS_VEC &N, 1474 Array<int> &nodeList) const 1475 { 1476 // get element id from global coordinates 1477 int eltID = map_to_element(x); 1478 1479 // call appropriate function below, with eltID 1480 shape_functions(x,eltID,N,nodeList); 1481 } 1482 shape_functions(const DENS_VEC & x,DENS_VEC & N,DENS_MAT & dNdx,Array<int> & nodeList) const1483 void FE_Mesh::shape_functions(const DENS_VEC &x, 1484 DENS_VEC &N, 1485 DENS_MAT &dNdx, 1486 Array<int> &nodeList) const 1487 { 1488 // get element id from global coordinates 1489 int eltID = map_to_element(x); 1490 1491 // call appropriate function below, with eltID 1492 shape_functions(x,eltID,N,dNdx,nodeList); 1493 } 1494 shape_functions(const VECTOR & x,const int eltID,DENS_VEC & N,Array<int> & nodeList) const1495 void FE_Mesh::shape_functions(const VECTOR &x, 1496 const int eltID, 1497 DENS_VEC &N, 1498 Array<int> &nodeList) const 1499 { 1500 // Get element node coordinates from mesh 1501 DENS_MAT eltCoords; 1502 element_coordinates(eltID, eltCoords); 1503 1504 // pass through call 1505 feElement_->shape_function(eltCoords,x,N); 1506 1507 // determine nodes which correspond to shape function indices 1508 element_connectivity_unique(eltID,nodeList); 1509 } 1510 shape_functions(const DENS_VEC & x,const int eltID,DENS_VEC & N,DENS_MAT & dNdx,Array<int> & nodeList) const1511 void FE_Mesh::shape_functions(const DENS_VEC &x, 1512 const int eltID, 1513 DENS_VEC &N, 1514 DENS_MAT &dNdx, 1515 Array<int> &nodeList) const 1516 { 1517 // Get element node coordinates from mesh 1518 DENS_MAT eltCoords; 1519 element_coordinates(eltID,eltCoords); 1520 1521 // pass through call 1522 feElement_->shape_function(eltCoords,x,N,dNdx); 1523 1524 // determine nodes which correspond to shp function indices 1525 element_connectivity_unique(eltID,nodeList); 1526 } 1527 shape_function_derivatives(const DENS_VEC & x,const int eltID,DENS_MAT & dNdx,Array<int> & nodeList) const1528 void FE_Mesh::shape_function_derivatives(const DENS_VEC &x, 1529 const int eltID, 1530 DENS_MAT &dNdx, 1531 Array<int> &nodeList) const 1532 { 1533 // Get element node coordinates from mesh 1534 DENS_MAT eltCoords; 1535 element_coordinates(eltID,eltCoords); 1536 1537 // pass through call 1538 feElement_->shape_function_derivatives(eltCoords,x,dNdx); 1539 1540 // determine nodes which correspond to shp function indices 1541 element_connectivity_unique(eltID,nodeList); 1542 } 1543 shape_function(const int eltID,DENS_MAT & N,DIAG_MAT & weights) const1544 void FE_Mesh::shape_function(const int eltID, 1545 DENS_MAT &N, 1546 DIAG_MAT &weights) const 1547 { 1548 // unused data (but required to calc weights) 1549 vector<DENS_MAT> dN; 1550 1551 // call below function with dN to avoid duplicated code 1552 shape_function(eltID,N,dN,weights); 1553 } 1554 shape_function(int eltID,DENS_MAT & N,vector<DENS_MAT> & dN,DIAG_MAT & weights) const1555 void FE_Mesh::shape_function(int eltID, 1556 DENS_MAT &N, 1557 vector<DENS_MAT> &dN, 1558 DIAG_MAT &weights) const 1559 { 1560 // Get element node coordinates from mesh 1561 DENS_MAT eltCoords; 1562 element_coordinates(eltID,eltCoords); 1563 1564 // pass through call 1565 feElement_->shape_function(eltCoords,N,dN,weights); 1566 } 1567 face_shape_function(const PAIR & face,DENS_MAT & N,DENS_MAT & n,DIAG_MAT & weights) const1568 void FE_Mesh::face_shape_function(const PAIR &face, 1569 DENS_MAT &N, 1570 DENS_MAT &n, 1571 DIAG_MAT &weights) const 1572 { 1573 int eltID = face.first; 1574 int faceID = face.second; 1575 1576 // Get element node coordinates from mesh 1577 DENS_MAT eltCoords; 1578 element_coordinates(eltID,eltCoords); 1579 1580 // pass through call 1581 feElement_->face_shape_function(eltCoords,faceID,N,n,weights); 1582 } 1583 face_shape_function(const PAIR & face,DENS_MAT & N,vector<DENS_MAT> & dN,vector<DENS_MAT> & Nn,DIAG_MAT & weights) const1584 void FE_Mesh::face_shape_function(const PAIR &face, 1585 DENS_MAT &N, 1586 vector<DENS_MAT> &dN, 1587 vector<DENS_MAT> &Nn, 1588 DIAG_MAT &weights) const 1589 { 1590 int eltID = face.first; 1591 int faceID = face.second; 1592 1593 // Get element node coordinates from mesh 1594 DENS_MAT eltCoords; 1595 element_coordinates(eltID,eltCoords); 1596 1597 // pass through call 1598 feElement_->face_shape_function(eltCoords,faceID,N,dN,Nn,weights); 1599 } 1600 face_normal(const PAIR & face,int ip,DENS_VEC & normal) const1601 double FE_Mesh::face_normal(const PAIR &face, 1602 int ip, 1603 DENS_VEC &normal) const 1604 { 1605 int eltID = face.first; 1606 int faceID = face.second; 1607 1608 // Get element node coordinates from mesh 1609 DENS_MAT eltCoords; 1610 element_coordinates(eltID,eltCoords); 1611 1612 // pass through call 1613 double J = feElement_->face_normal(eltCoords,faceID,ip,normal); 1614 return J; 1615 } 1616 1617 //----------------------------------------------------------------------- output(string prefix) const1618 void FE_Mesh::output(string prefix) const 1619 { 1620 set<int> otypes; 1621 otypes.insert(ENSIGHT); 1622 otypes.insert(FULL_GNUPLOT); 1623 OutputManager meshSets(prefix,otypes); 1624 meshSets.write_geometry(&nodalCoords_, & connectivity_); 1625 OUTPUT_LIST subsetData; 1626 int size = nNodesUnique_; 1627 // material 1628 // DENS_MAT material(nNodes_,1); 1629 // material = 1; 1630 // subsetData["material"] = &material; 1631 //string name = "global_to_unique_map"; 1632 // nodesets 1633 NODE_SET_MAP::const_iterator niter; 1634 map<string,DENS_MAT> nodesets; 1635 for (niter = nodeSetMap_.begin(); niter != nodeSetMap_.end(); niter++) { 1636 string name = niter->first; 1637 const set<int> & nset = niter->second; 1638 string nodeset = "nodeset_"+name; 1639 nodesets[nodeset].reset(size,1); 1640 set<int>::const_iterator iter; 1641 for (iter = nset.begin(); iter != nset.end(); iter++) { 1642 (nodesets[nodeset])(*iter,0) = 1; 1643 } 1644 subsetData[nodeset] = & nodesets[nodeset]; 1645 } 1646 // facesets 1647 FACE_SET_MAP::const_iterator fiter; 1648 map<string,DENS_MAT> facesets; 1649 for (fiter = faceSetMap_.begin(); fiter != faceSetMap_.end(); fiter++) { 1650 string name = fiter->first; 1651 string faceset = "faceset_"+name; 1652 facesets[faceset].reset(size,1); 1653 set<int> nset; 1654 faceset_to_nodeset(name,nset); 1655 set<int>::const_iterator iter; 1656 for (iter = nset.begin(); iter != nset.end(); iter++) { 1657 (facesets[faceset])(*iter,0) = 1; 1658 } 1659 subsetData[faceset] = & facesets[faceset]; 1660 } 1661 // elementsets 1662 ELEMENT_SET_MAP::const_iterator eiter; 1663 map<string,DENS_MAT> elemsets; 1664 for (eiter = elementSetMap_.begin();eiter != elementSetMap_.end();eiter++) { 1665 string name = eiter->first; 1666 const set<int> & eset = eiter->second; 1667 string elemset = "elementset_"+name; 1668 elemsets[elemset].reset(size,1); 1669 set<int>::const_iterator iter; 1670 for (iter = eset.begin(); iter != eset.end(); iter++) { 1671 Array<int> nodes; 1672 element_connectivity_unique(*iter, nodes); 1673 for(int i = 0; i < nodes.size(); ++i) { 1674 (elemsets[elemset])(nodes(i),0) = 1; 1675 } 1676 } 1677 subsetData[elemset] = & elemsets[elemset]; 1678 } 1679 // output 1680 const int * map = globalToUniqueMap_.data(); 1681 if (subsetData.size() > 0 ) 1682 meshSets.write_data(0.0,&subsetData,map); 1683 else 1684 if (LammpsInterface::instance()->comm_rank() == 0) { 1685 stringstream ss; 1686 ss << "Warning mesh output requested without any mesh entities, output suppressed"; 1687 ATC::LammpsInterface::instance()->print_msg(ss.str()); 1688 } 1689 } 1690 is_owned_elt(int elt) const1691 bool FE_Mesh::is_owned_elt(int elt) const 1692 { 1693 return (find(myElts_.begin(), myElts_.end(), elt) != myElts_.end()); 1694 } 1695 1696 1697 // ------------------------------------------------------------- 1698 // ------------------------------------------------------------- 1699 // class FE_3DMesh 1700 // ------------------------------------------------------------- 1701 // ------------------------------------------------------------- FE_3DMesh(const string elementType,const int nNodes,const int nElements,const Array2D<int> * connectivity,const DENS_MAT * nodalCoordinates,const Array<bool> periodicity,const Array<pair<string,set<int>>> * nodeSets)1702 FE_3DMesh::FE_3DMesh(const string elementType, 1703 const int nNodes, 1704 const int nElements, 1705 const Array2D<int> *connectivity, 1706 const DENS_MAT *nodalCoordinates, 1707 const Array<bool> periodicity, 1708 const Array< pair< string, set<int> > > *nodeSets): 1709 FE_Mesh(), 1710 minEltSize_(0), 1711 tree_(nullptr) 1712 { 1713 // Pick which element class to make 1714 if (elementType == "HEX8") { 1715 feElement_ = new FE_ElementHex(8,4,2); 1716 } else if (elementType == "HEX20") { 1717 feElement_ = new FE_ElementHex(20,8,3); 1718 } else if (elementType == "HEX27") { 1719 feElement_ = new FE_ElementHex(27,9,3); 1720 } else if (elementType == "TET4") { 1721 feElement_ = new FE_ElementTet(4,3,2); 1722 hasPlanarFaces_ = true; 1723 } else { 1724 throw ATC_Error("Unrecognized element type specified."); 1725 } 1726 1727 1728 nSD_ = 3; 1729 nNodes_ = nNodes; 1730 nNodesUnique_ = nNodes; 1731 nElts_ = nElements; 1732 xscale_ = yscale_ = zscale_ = 1; 1733 periodicity_ = periodicity; 1734 1735 // Connectivity and coordinates describe the mesh geometry. 1736 connectivity_.reset(connectivity->nRows(),connectivity->nCols()); 1737 connectivity_ = (*connectivity); 1738 nodalCoords_.reset(nodalCoordinates->nRows(),nodalCoordinates->nCols()); 1739 nodalCoords_ = (*nodalCoordinates); 1740 1741 // set minimum element size 1742 minEltSize_ = 1.e20; 1743 for (int i=0; i< connectivity_.nCols(); ++i) { 1744 int n1 = connectivity_(0,i); 1745 int n2 = connectivity_(1,i); 1746 double dx[3] = {fabs(nodalCoords_(0,n1)-nodalCoords_(0,n2)), 1747 fabs(nodalCoords_(1,n1)-nodalCoords_(1,n2)), 1748 fabs(nodalCoords_(2,n1)-nodalCoords_(2,n2))}; 1749 minEltSize_ = min(minEltSize_,norm3(dx)); 1750 } 1751 1752 // create global-unique maps 1753 coordTol_ = this->coordinate_tolerance(); 1754 setup_periodicity(); 1755 1756 // Create the "all" elementset, the "all" nodeset, and the read-in nodesets. 1757 for (int elem = 0; elem < nElts_; elem++) elementSetAll_.insert(elem); 1758 for (int node = 0; node < nNodesUnique_; node++) nodeSetAll_.insert(node); 1759 const Array<pair<string,set<int> > > & sets = *nodeSets; 1760 for (int nodeSet = 0; nodeSet < sets.size(); ++nodeSet) { 1761 const set<int> & nset = sets(nodeSet).second; 1762 set<int> copy; 1763 if (compactRemap_.size() > 0) { 1764 for (set<int>::iterator itr = nset.begin(); itr != nset.end(); itr++) { 1765 copy.insert(globalToUniqueMap_(compactRemap_(*itr))); 1766 } 1767 } 1768 else { 1769 for (set<int>::iterator itr = nset.begin(); itr != nset.end(); itr++) { 1770 copy.insert(globalToUniqueMap_(*itr)); 1771 } 1772 } 1773 create_nodeset(sets(nodeSet).first, copy); 1774 } 1775 1776 // Insert nodes and elements into KD-tree for PIE search. 1777 if (tree_ == nullptr) { 1778 tree_ = KD_Tree::create_KD_tree(feElement_->num_elt_nodes(), nNodes_, 1779 &nodalCoords_, nElts_, connectivity_); 1780 } 1781 1782 } 1783 ~FE_3DMesh()1784 FE_3DMesh::~FE_3DMesh() { 1785 if (tree_) delete tree_; 1786 } 1787 1788 // ------------------------------------------------------------- 1789 // setup_periodicity 1790 // ------------------------------------------------------------- setup_periodicity(double)1791 void FE_3DMesh::setup_periodicity(double /* tol */) 1792 { 1793 // unique <-> global id maps 1794 globalToUniqueMap_.reset(nNodes_); 1795 uniqueToGlobalMap_.reset(nNodes_); 1796 1797 for (int node = 0; node < nNodes_; node++) { 1798 globalToUniqueMap_(node) = node; 1799 } 1800 1801 // recursive fix and setup global to unique map 1802 if (periodicity_(2)) { fix_periodicity(3); } 1803 if (periodicity_(1)) { fix_periodicity(2); } 1804 if (periodicity_(0)) { fix_periodicity(1); } 1805 1806 1807 // renumber to compact unique numbering 1808 // unique nodes map to the same id with global to unique 1809 if (nNodesUnique_ < nNodes_) { 1810 int n = 0; 1811 int m = nNodesUnique_; 1812 compactRemap_.reset(nNodes_); // old to new global numbering 1813 compactRemap_ = -1; 1814 for (int node = 0; node < nNodes_; node++) { 1815 if (globalToUniqueMap_(node) == node) { // unique nodes 1816 compactRemap_(node) = n++; 1817 } 1818 else { // periodic nodes 1819 compactRemap_(node) = m++; 1820 } 1821 } 1822 if (n != nNodesUnique_) throw ATC_Error("didn't compact numbering"); 1823 int npe = num_nodes_per_element(); 1824 // temporary copy 1825 DENS_MAT coor = DENS_MAT(nodalCoords_); 1826 Array2D<int> conn(connectivity_); 1827 Array<int> oldGlobalToUniqueMap = globalToUniqueMap_; 1828 for (int node = 0; node < nNodes_; node++) { 1829 // unique = remap * global2unique global 1830 globalToUniqueMap_(compactRemap_(node)) = compactRemap_(oldGlobalToUniqueMap(node)); 1831 if (compactRemap_(node) < 0) throw ATC_Error("mis-map to compact numbering"); 1832 // swap coordinates 1833 for (int i = 0; i < nSD_; i++) { 1834 nodalCoords_(i,compactRemap_(node)) = coor(i,node); 1835 } 1836 } 1837 for (int elem = 0; elem < nElts_; elem++) { 1838 for (int i = 0; i < npe; i++) { 1839 connectivity_(i,elem) = compactRemap_(conn(i,elem)); 1840 } 1841 } 1842 } 1843 1844 for (int node = 0; node < nNodes_; node++) { 1845 uniqueToGlobalMap_(globalToUniqueMap_(node)) = node; 1846 } 1847 set_unique_connectivity(); 1848 //amended_ = true; // would always be true since setup_always called 1849 } 1850 fix_periodicity(int idir)1851 void FE_3DMesh::fix_periodicity(int idir) 1852 { 1853 set<int> topNodes,botNodes; 1854 int ntop = find_boundary_nodes( idir,topNodes); 1855 int nbot = find_boundary_nodes(-idir,botNodes); 1856 if (ntop != nbot) 1857 throw ATC_Error("can't fix periodicity, number of top and bottom nodes are different "); 1858 bool match = match_nodes(idir,topNodes,botNodes,globalToUniqueMap_); 1859 if (!match) { 1860 stringstream ss; 1861 ss << "can't match periodic nodes with tolerance " << coordTol_; 1862 throw ATC_Error(ss.str()); 1863 } 1864 } 1865 find_boundary_nodes(int idir,set<int> & nodes)1866 int FE_3DMesh::find_boundary_nodes(int idir, set<int> & nodes) 1867 { 1868 nodes.clear(); 1869 double limit = 0; 1870 int idm = abs(idir)-1; 1871 DENS_MAT & x = nodalCoords_; 1872 if (idir > 0) limit = x.row_max(idm); 1873 else limit = x.row_min(idm); 1874 for (int i=0; i < x.nCols(); ++i) { 1875 double xi = x(idm,i); 1876 if (fabs(xi-limit) < coordTol_) nodes.insert(i); 1877 } 1878 // stringstream ss; 1879 // ss << "found " << nodes.size() << " nodes at x_" << abs(idir) << "= "<< limit ; 1880 // ATC::LammpsInterface::instance()->print_msg_once(ss.str()); 1881 return nodes.size(); 1882 } 1883 match_nodes(int idir,set<int> & nodes1,set<int> & nodes2,Array<int> & map)1884 bool FE_3DMesh::match_nodes(int idir, set<int> & nodes1, set<int> & nodes2, 1885 Array<int> & map) 1886 { 1887 int i1=0,i2=1; 1888 plane_coords(idir-1,i1,i2); 1889 vector<bool> found(nNodes_,false); 1890 DENS_MAT & x = nodalCoords_; 1891 for (set<int>::iterator it1=nodes1.begin(); it1 != nodes1.end(); it1++) { 1892 int n1 = *it1; 1893 double x1 = x(i1,n1); 1894 double x2 = x(i2,n1); 1895 for (set<int>::iterator it2=nodes2.begin(); it2 != nodes2.end(); it2++) { 1896 int n2 = *it2; 1897 if (!found[n2]) { 1898 double y1 = x(i1,n2); 1899 double y2 = x(i2,n2); 1900 if (fabs(x1-y1) < coordTol_ && fabs(x2-y2) < coordTol_) { 1901 map(n1) = n2; 1902 found[n2] = true; 1903 // coincidence 1904 x(i1,n2) = x1; 1905 x(i2,n2) = x2; 1906 } 1907 } 1908 } 1909 if (map(n1) == n1) return false; 1910 } 1911 nNodesUnique_ -= nodes1.size(); 1912 stringstream ss; 1913 ss << "condensed " << nodes1.size() << " periodic nodes in the " << abs(idir) << " direction"; 1914 ATC::LammpsInterface::instance()->print_msg_once(ss.str()); 1915 return true; 1916 } 1917 set_unique_connectivity(void)1918 void FE_3DMesh::set_unique_connectivity(void) 1919 { 1920 int numEltNodes = feElement_->num_elt_nodes(); 1921 connectivityUnique_.reset(numEltNodes, nElts_); 1922 uniqueNodeToElementMap_.reset(nNodes_); 1923 for (int node = 0; node < nNodes_; node++) { 1924 uniqueNodeToElementMap_(node) = vector<int>(); 1925 } 1926 for (int elem = 0; elem < nElts_; elem++) { 1927 for (int node = 0; node < numEltNodes; node++) { 1928 int global_node = connectivity_(node, elem); 1929 int unique_node = globalToUniqueMap_(global_node); 1930 connectivityUnique_(node,elem) = unique_node; 1931 uniqueNodeToElementMap_(unique_node).push_back(elem); 1932 } 1933 } 1934 } 1935 1936 // orient the local coordinate with the global one orient(int idm)1937 bool FE_3DMesh::orient(int idm) 1938 { 1939 int numEltNodes = feElement_->num_elt_nodes(); 1940 if (numEltNodes != 8) throw ATC_Error("can't currently orient non HEX8 elements"); 1941 DENS_MAT x; 1942 for (int elem = 0; elem < nElts_; elem++) { 1943 element_coordinates(elem,x); 1944 double xmax = x.row_max(idm); 1945 double xmin = x.row_min(idm); 1946 set<int> top,bot; 1947 for (int node = 0; node < numEltNodes; node++) { 1948 // find top nodes 1949 if ((x(idm,node) - xmax) < coordTol_) top.insert(node); 1950 // find bottom nodes 1951 else if ((x(idm,node) - xmin) < coordTol_) bot.insert(node); 1952 else return false; 1953 } 1954 // order by rh rule 1955 // start with one 1956 } 1957 throw ATC_Error("not completely implemented function: FE_3DMesh::orient"); 1958 return true; 1959 } 1960 1961 // ------------------------------------------------------------- 1962 // amend mesh for cut at specified faces 1963 // ------------------------------------------------------------- cut_mesh(const set<PAIR> & faceSet,const set<int> & nodeSet)1964 void FE_3DMesh::cut_mesh(const set<PAIR> & faceSet, const set<int> & nodeSet) 1965 { 1966 int nsd = nSD_; 1967 // get all nodes to create an old to new number map 1968 set<int> dupNodes; 1969 map<int,int> oldToNewMap; 1970 faceset_to_nodeset_global(faceSet,dupNodes); 1971 // remove edge nodes 1972 Array<int> & node_map = globalToUniqueMap_; 1973 set<int>::const_iterator itr; 1974 for (itr = dupNodes.begin(); itr != dupNodes.end(); itr++) { 1975 int gnode = *itr; 1976 int unode = node_map(gnode); 1977 if (nodeSet.find(unode) != nodeSet.end()) { 1978 dupNodes.erase(gnode); 1979 } 1980 } 1981 int nNodesAdded = dupNodes.size(); 1982 // copy coordinates and nodemap 1983 DENS_MAT &coordinates = nodalCoords_; 1984 int nNodesNew = coordinates.nCols() + nNodesAdded; 1985 coordinates.resize(nsd,nNodesNew,true); 1986 node_map.resize(nNodesNew,true); 1987 // add duplicate coordinates 1988 int iNew = nNodes_; 1989 int iNewUnique = nNodesUnique_; 1990 for (itr = dupNodes.begin(); itr != dupNodes.end(); 1991 itr++,iNew++) { 1992 int iOld = *itr; 1993 oldToNewMap[iOld] = iNew; // global ids 1994 if (iOld == node_map(iOld)) { // non-image atom 1995 node_map(iNew) = iNewUnique++; 1996 } else { 1997 node_map(iNew) = -1; 1998 } 1999 for(int j = 0; j < nsd; j++) { 2000 coordinates(j,iNew) = coordinates(j,iOld); 2001 } 2002 } 2003 nNodes_ = iNew; 2004 nNodesUnique_ = iNewUnique; 2005 for (itr = dupNodes.begin(); itr != dupNodes.end(); 2006 itr++,iNew++) { 2007 int iOld = *itr; 2008 iNew = oldToNewMap[iOld]; // global ids 2009 int iOldImage = node_map(iOld); 2010 if (iOld != iOldImage) { // image atom 2011 int iNewImage = oldToNewMap[iOldImage]; 2012 node_map(iNew) = node_map(iNewImage); 2013 } 2014 } 2015 // update element connectivities 2016 const int nnf = num_nodes_per_face(); 2017 const Array2D <int> & local_conn = feElement_->local_face_conn(); 2018 set< PAIR >::iterator iter; 2019 for (iter = faceSet.begin(); iter != faceSet.end(); iter++) 2020 { 2021 PAIR face = *iter; 2022 int eltID=face.first, faceID=face.second; 2023 for (int inode=0; inode < nnf; inode++) { 2024 int lid = local_conn(faceID,inode); 2025 int id = connectivity_(lid, eltID); 2026 if (oldToNewMap.find(id) != oldToNewMap.end() ) { 2027 int new_id = (*oldToNewMap.find(id)).second; 2028 connectivity_(lid,eltID) = new_id; 2029 connectivityUnique_(lid, eltID) = node_map(new_id); 2030 } 2031 } 2032 } 2033 } 2034 2035 2036 // ------------------------------------------------------------- 2037 // amend mesh for deleted elements 2038 // ------------------------------------------------------------- delete_elements(const set<int> & elementList)2039 void FE_3DMesh::delete_elements(const set<int> &elementList) 2040 { 2041 int nPE = num_nodes_per_element(); 2042 set<int> elementsNew; 2043 elementset_complement(elementList,elementsNew); 2044 int nElementsNew = elementsNew.size(); 2045 set<int> newToOld; 2046 map<int,int> oldToNewMap; 2047 elementset_to_nodeset(elementsNew,newToOld); 2048 int nNodesNew = newToOld.size(); 2049 set<int>::const_iterator itr; 2050 2051 // coordinates & node map (from nodes to data) 2052 const DENS_MAT &coordinates = nodal_coordinates(); 2053 const Array<int> & node_map = global_to_unique_map(); 2054 2055 DENS_MAT *newCoords = new DENS_MAT(nSD_,nNodesNew); 2056 Array<int> *newNodeMap = new Array<int> (nNodesNew); 2057 Array2D<int> * newConnectivity = new Array2D<int>(nPE,nElementsNew); 2058 int k = 0, i = 0; 2059 for (itr = newToOld.begin(); itr != newToOld.end(); itr++) { 2060 int node = *itr; 2061 oldToNewMap[node] = i++; 2062 (*newNodeMap)(k) = node_map(node); 2063 for(int j = 0; j < nSD_; j++) { 2064 (*newCoords)(j,k) = coordinates(j,node); 2065 } 2066 k++; 2067 } 2068 // nNodes_ = nNodesNew; ??? 2069 // connectivity 2070 k = 0; 2071 for (itr = elementsNew.begin(); itr != elementsNew.end(); itr++) { 2072 int ielem = *itr; 2073 for(int j = 0; j < nPE; j++) { 2074 int old_node = connectivity_(j,ielem); 2075 map<int,int>::iterator map_itr = oldToNewMap.find(old_node); 2076 if (map_itr == oldToNewMap.end()) { 2077 stringstream ss; 2078 ss << "map failure " << old_node << "\n"; 2079 ATC::LammpsInterface::instance()->print_msg(ss.str()); 2080 } 2081 int node = map_itr->second; 2082 (*newConnectivity)(j,k) = node; 2083 } 2084 k++; 2085 } 2086 delete newCoords; 2087 delete newNodeMap; 2088 delete newConnectivity; 2089 } 2090 2091 // ------------------------------------------------------------- 2092 // partition_mesh 2093 // ------------------------------------------------------------- 2094 partition_mesh()2095 void FE_3DMesh::partition_mesh() 2096 { 2097 if (lammpsPartition_) { 2098 lammps_partition_mesh(); 2099 } 2100 2101 if (partitioned_) return; 2102 2103 2104 int nProcs, myrank; 2105 MPI_Comm_size(MPI_COMM_WORLD, &nProcs); // get the number of processors 2106 MPI_Comm_rank(MPI_COMM_WORLD, &myrank); // get the current processor's rank 2107 2108 // use the KD tree for partitioning, getting more blocks than 2109 // processors 2110 if (tree_ == nullptr) { 2111 tree_ = KD_Tree::create_KD_tree(feElement_->num_elt_nodes(), 2112 nNodes_, &nodalCoords_, 2113 nElts_, connectivity_); 2114 } 2115 2116 // Get the divisions of the elements from the KD-tree. 2117 int depth = ceil(log(nProcs)/log(2)); 2118 // In order to make sure there are enough divisions to evenly 2119 // divide between all processors, we get the next-highest 2120 // power of 2. 2121 vector<vector<int> > procEltLists = tree_->getElemIDs(depth); 2122 2123 // Make sure the KD tree is behaving as expected. 2124 assert(procEltLists.size() >= nProcs); 2125 2126 // If the KD-tree was not able to return enough divisions, 2127 // duplicate the largest list. 2128 2129 // elements, then the division would be more even. 2130 vector<vector<int> >::iterator it; 2131 if (numNonempty(procEltLists) < nProcs) { 2132 // Find the list of maximum size and assign it to empty processors 2133 vector<int> maxSizeList (0); 2134 for (it = procEltLists.begin(); it != procEltLists.end(); ++it) { 2135 if (it->size() > maxSizeList.size()) 2136 maxSizeList.assign(it->begin(), it->end()); 2137 } 2138 for (it = procEltLists.begin(); it != procEltLists.end(); ++it) { 2139 if (it->empty()) { 2140 if (numNonempty(procEltLists) >= nProcs) break; 2141 it->assign(maxSizeList.begin(), maxSizeList.end()); 2142 } 2143 } 2144 } 2145 2146 // We will store the owning processor for each element. 2147 int * eltToOwners = new int[nElts_]; 2148 for (int i = 0; i < nElts_; ++i) { 2149 eltToOwners[i] = -1; // -1 here means the element is unowned. 2150 } 2151 2152 // Prune elements that appear on more than one processor. 2153 // Simultaneously fill in the ownership array. 2154 prune_duplicate_elements( procEltLists, eltToOwners); 2155 2156 // If we have more lists than processors, get rid of the 2157 // extras and redistribute their elements. 2158 if (numNonempty(procEltLists) > nProcs) { 2159 redistribute_extra_proclists(procEltLists, eltToOwners, nProcs); 2160 } 2161 2162 // Sort the lists so that the fuller ones are in the front. 2163 sort(procEltLists.begin(), procEltLists.end(), vectorCompSize); 2164 2165 // Assign each processor a list of elements. 2166 myElts_ = procEltLists[myrank]; 2167 2168 //mesh_distribute(eltStartIDs); 2169 delete[] eltToOwners; 2170 2171 // We should do nodes later. 2172 2173 if (decomposition_) distribute_mesh_data(); 2174 partitioned_ = true; 2175 } 2176 departition_mesh()2177 void FE_3DMesh::departition_mesh() 2178 { 2179 if (!partitioned_) return; 2180 partitioned_ = false; 2181 } 2182 prune_duplicate_elements(vector<vector<int>> & procEltLists,int * eltToOwners)2183 void FE_3DMesh::prune_duplicate_elements(vector<vector<int> > & procEltLists, 2184 int * eltToOwners) 2185 { 2186 int procID = 0; 2187 vector<vector<int> >::iterator topIt; 2188 vector<int> * conflictingProc; 2189 vector<int>::iterator it, toErase; 2190 for (topIt = procEltLists.begin(); topIt != procEltLists.end(); ++topIt) { 2191 // Simultaneously fill in eltToOwners and prune, if an element 2192 // appears on multiple processors. 2193 for (it = topIt->begin(); it != topIt->end(); ++it) { 2194 // If the element has no corresponding processor in eltToOwners, 2195 // record it as belonging to processor *it. 2196 if (eltToOwners[*it] == -1) { 2197 eltToOwners[*it] = procID; 2198 } 2199 else { 2200 // If it does have a processor in eltToOwners, then we need 2201 // to remove it from either processor *topIt or from the processor 2202 // listed in eltToOwners. We discriminate based on size. 2203 conflictingProc = &(procEltLists[eltToOwners[*it]]); 2204 if (conflictingProc->size() <= procEltLists[procID].size()) { 2205 // Delete element from processor *topIt, if it has more elements. 2206 it = topIt->erase(it); 2207 --it; 2208 } 2209 else { 2210 // Delete element from conflicting processor otherwise. 2211 toErase = find(conflictingProc->begin(), conflictingProc->end(), *it); 2212 conflictingProc->erase(toErase); 2213 eltToOwners[*it] = procID; 2214 } 2215 } 2216 } 2217 ++procID; 2218 } 2219 } 2220 2221 // ------------------------------------------------------------- 2222 // lammps_partition_mesh 2223 // ------------------------------------------------------------- 2224 lammps_partition_mesh()2225 void FE_3DMesh::lammps_partition_mesh() 2226 { 2227 if (LammpsInterface::instance()->domain_triclinic()) { 2228 LammpsInterface::instance()->print_msg_once("Cannot use lammps partitioning, domain is triclinic"); 2229 return; 2230 } 2231 LammpsInterface::instance()->print_msg_once("Warning: Using native lammps partitioning"); 2232 double xlo, xhi, ylo, yhi, zlo, zhi; 2233 double boxxlo, boxxhi, boxylo, boxyhi, boxzlo, boxzhi; 2234 LammpsInterface::instance()->sub_bounds(xlo, xhi, ylo, yhi, zlo, zhi); 2235 LammpsInterface::instance()->box_bounds(boxxlo, boxxhi, boxylo, boxyhi, boxzlo, boxzhi); 2236 2237 2238 myElts_.clear(); 2239 double xCent, yCent, zCent; 2240 2241 // Assign elements to processors based on the centroid of the element. 2242 int numNodes = num_nodes_per_element(); 2243 for (int i = 0; i < nElts_; ++i) 2244 { 2245 xCent = 0.0; 2246 yCent = 0.0; 2247 zCent = 0.0; 2248 for (int j = 0; j < numNodes; ++ j) 2249 { 2250 xCent += nodalCoords_(0, connectivity_(j,i)); 2251 yCent += nodalCoords_(1, connectivity_(j,i)); 2252 zCent += nodalCoords_(2, connectivity_(j,i)); 2253 } 2254 xCent /= numNodes; 2255 yCent /= numNodes; 2256 zCent /= numNodes; 2257 if (xCent < boxxlo) xCent = boxxlo; 2258 if (xCent < boxxhi) xCent = boxxhi; 2259 if (yCent < boxylo) yCent = boxylo; 2260 if (yCent < boxyhi) yCent = boxyhi; 2261 if (zCent < boxzlo) zCent = boxzlo; 2262 if (zCent < boxzhi) zCent = boxzhi; 2263 if ( dbl_geq(xCent, xlo) && 2264 ((xhi == boxxhi) || !dbl_geq(xCent, xhi)) && 2265 dbl_geq(yCent, ylo) && 2266 ((yhi == boxyhi) || !dbl_geq(yCent, yhi)) && 2267 dbl_geq(zCent, zlo) && 2268 ((zhi == boxzhi) || !dbl_geq(zCent, zhi))) { 2269 myElts_.push_back(i); 2270 } 2271 } 2272 2273 2274 // if decomposing add in myAtomElts list based on nodal locations, i.e., all elements with a local node 2275 // myElts: for FE assembly 2276 // myAndGhost : for atom ops like restrict 2277 if (decomposition_) { 2278 set<int> elms; 2279 vector<int>::const_iterator itr; 2280 for (itr=myElts_.begin(); itr!=myElts_.end(); itr++) {elms.insert(*itr); } 2281 set<int> nodes; 2282 elementset_to_nodeset(elms,nodes); 2283 elms.clear(); 2284 nodeset_to_maximal_elementset(nodes,elms); 2285 myAndGhostElts_.clear(); 2286 set<int>::const_iterator iter; 2287 for (iter=elms.begin(); iter!=elms.end(); iter++) 2288 {myAndGhostElts_.push_back(*iter);} 2289 distribute_mesh_data(); 2290 } 2291 partitioned_ = true; 2292 return; 2293 2294 } 2295 redistribute_extra_proclists(vector<vector<int>> & procEltLists,int * eltToOwners,int nProcs)2296 void FE_3DMesh::redistribute_extra_proclists(vector<vector<int> > &procEltLists, 2297 int *eltToOwners, int nProcs) 2298 { 2299 DENS_MAT faceAdjacencies(nElts_, num_faces_per_element()); 2300 faceAdjacencies = -1; // Set all values to -1, indicating uninitialized/uncalculated 2301 2302 int currentElt, adjacentElt, procID; 2303 2304 // Put all of the hobos onto one master list, allHomelessElts. 2305 list<int> allHomelessElts; 2306 vector<int> oneHomelessList; 2307 vector<vector<int> >::iterator current; 2308 int nHoboLists = numNonempty(procEltLists) - nProcs; 2309 for (int i = 0; i < nHoboLists; ++i) { 2310 current = min_element(procEltLists.begin(), procEltLists.end(), vectorCompSize); 2311 oneHomelessList = *current; 2312 allHomelessElts.insert(allHomelessElts.end(), 2313 oneHomelessList.begin(), oneHomelessList.end()); 2314 current->clear(); 2315 } 2316 2317 // Make sure the hobos lose their association with their old processor. 2318 list<int>::iterator it; 2319 for (it = allHomelessElts.begin(); it != allHomelessElts.end(); it++){ 2320 eltToOwners[*it] = -1; 2321 } 2322 2323 // Figure out which elements the hobos are adjacent to. That way, they 2324 // will know what processors they can be redistributed to. 2325 compute_face_adjacencies(allHomelessElts, faceAdjacencies); 2326 2327 // Place homeless elements onto lists that correspond to actual processors. 2328 while (!allHomelessElts.empty()) { 2329 currentElt = allHomelessElts.back(); 2330 2331 // This will store the ID of the processor with the fewest elements 2332 // so far that has an element adjacent to currentElt. 2333 PAIR smallestProc(-1, INT_MAX); 2334 2335 // Iterate over the faces, check the processors of adjacent elements, 2336 // and slate the element to go on the adjacent processor with the fewest 2337 // elements. 2338 for (int localFaceID = 0; localFaceID < num_faces_per_element(); ++localFaceID) { 2339 adjacentElt = faceAdjacencies(currentElt, localFaceID); 2340 2341 // This means that there is no adjacency through this face. 2342 if (adjacentElt >= nElts_) continue; 2343 2344 procID = eltToOwners[adjacentElt]; 2345 // The procID > -1 check makes sure we're not adjacent to another 2346 // homeless element by this face, in which case it won't have a 2347 // processor to put currentElt onto yet. 2348 if (procID > -1 && ((int) procEltLists[procID].size()) < smallestProc.second) { 2349 smallestProc = PAIR(procID, procEltLists[procID].size()); 2350 } 2351 } 2352 2353 allHomelessElts.pop_back(); 2354 2355 // If we couldn't find an adjacent element that had a processor, 2356 // skip for now and come back to it later. 2357 if (smallestProc.first == -1) { 2358 allHomelessElts.push_front(currentElt); 2359 } 2360 // Otherwise, put it onto the processor with the fewest elements that 2361 // we found. 2362 else { 2363 procEltLists[smallestProc.first].push_back(currentElt); 2364 eltToOwners[currentElt] = smallestProc.first; 2365 } 2366 } 2367 } 2368 2369 compute_face_adjacencies(const list<int> & elts,DENS_MAT & faceAdjacencies)2370 void FE_3DMesh::compute_face_adjacencies(const list<int> &elts, 2371 DENS_MAT &faceAdjacencies) 2372 { 2373 list<int>::const_iterator cit; 2374 for (cit = elts.begin(); cit != elts.end(); ++cit) { 2375 // For each element, look at ever face, get the nodes on that face, 2376 // and find out which elements are in the intersection of elements 2377 // containing that node. Those two elements are adjacent, and we 2378 // mark it as such in the faceAdjacencies array. 2379 for (int localFaceID = 0; localFaceID < num_faces_per_element(); ++localFaceID) { 2380 Array<int> faceNodes; 2381 face_connectivity(PAIR(*(cit), localFaceID), faceNodes); 2382 // Put the first node's elements into the accumulator to start. 2383 vector<int> vIntersect = uniqueNodeToElementMap_(faceNodes(0)); 2384 vector<int> vCurrent; 2385 // set_intersect requires a vector large enough to contain the 2386 // max possible intersect, which cannot be larger than the entirety 2387 // of the first vector involved. 2388 vector<int> vTemp(vIntersect.size(), -1); 2389 // Find the intersection of each of the nodes' element vectors. 2390 for (int ithOnFace = 1; ithOnFace < num_nodes_per_face(); ++ithOnFace) { 2391 vCurrent = uniqueNodeToElementMap_(faceNodes(ithOnFace)); // Vector of elements for this node 2392 set_intersection(vCurrent.begin(), vCurrent.end(), 2393 vIntersect.begin(), vIntersect.end(), vTemp.begin()); 2394 vIntersect = vTemp; 2395 // Because we initialized the vector to be larger than necessary, maybe, 2396 // we remove all of the meaningless values used in initialization. 2397 while (vIntersect.back() == -1) 2398 vIntersect.pop_back(); 2399 vTemp.clear(); 2400 vTemp.resize(vIntersect.size(),-1); 2401 } 2402 // This means there is an adjacent face, since there 2403 // are two elements sharing all of the nodes on the face. 2404 if (vIntersect.size() == 2) { 2405 // We want to choose the element id of NOT the current 2406 // element to be listed as the adjacency. 2407 2408 // well, but that requires more complicated memoization and 2409 // this doesn't take much extra time. 2410 if (*cit == vIntersect[0]) { 2411 faceAdjacencies(*cit, localFaceID) = vIntersect[1]; 2412 } 2413 else { 2414 faceAdjacencies(*cit, localFaceID) = vIntersect[0]; 2415 } 2416 } 2417 // This means the element is on the border. 2418 else if (vIntersect.size() == 1) { 2419 faceAdjacencies(*cit, localFaceID) = INT_MAX; 2420 } 2421 else { 2422 // This should never ever happen! The nodes should at least 2423 // share one element, since they are all on one face! 2424 // There should also never be more than two elements on the 2425 // same face... that would defy mathematics and physics in 2426 // every way. 2427 } 2428 } 2429 } 2430 } 2431 2432 // Sometimes we want to count the number of vectors that actually 2433 // have stuff in them. We use this. numNonempty(vector<vector<int>> & v)2434 int FE_3DMesh::numNonempty(vector<vector<int> > & v) 2435 { 2436 int result = 0; 2437 vector<vector<int> >::iterator it; 2438 for (it = v.begin(); it != v.end(); ++it){ 2439 if (!(it->empty())){ 2440 result++; 2441 } 2442 } 2443 return result; 2444 } 2445 map_to_element(const DENS_VEC & query) const2446 int FE_3DMesh::map_to_element(const DENS_VEC &query) const 2447 { 2448 DENS_MAT eltCoords; 2449 vector<int> candidates = tree_->find_nearest(query); 2450 vector<int> matches = vector<int>(); 2451 2452 // Search through each of the nearest elements 2453 for (vector<int>::iterator elem = candidates.begin(); 2454 elem < candidates.end(); elem++) { 2455 if (contains_point(*elem, query)) { 2456 matches.push_back(*elem); // Keep track of the elements 2457 // which contain it 2458 } 2459 } 2460 2461 // Return the index of the parent element which does contain the point 2462 if (matches.size() == 1) { // x is conclusively in a single element 2463 return matches[0]; 2464 } else if (matches.size() == 0) { // not so much 2465 throw ATC_Error("FE_3DMesh::map_to_element could not find an element"); 2466 } else { // definitely not so much 2467 throw ATC_Error("FE_3DMesh::map_to_element found multiple elements"); 2468 } 2469 } 2470 contains_point(const int eltID,const DENS_VEC & x) const2471 bool FE_3DMesh::contains_point(const int eltID, 2472 const DENS_VEC &x) const 2473 { 2474 DENS_MAT eltCoords; 2475 element_coordinates(eltID, eltCoords); 2476 return feElement_->contains_point(eltCoords, x); 2477 } 2478 2479 //----------------------------------------------------------------------- distribute_mesh_data()2480 void FE_3DMesh::distribute_mesh_data() 2481 { 2482 myNElts_ = myElts_.size(); 2483 2484 //create elemToMyElemMap_ 2485 elemToMyElemMap_.clear(); 2486 for (int myID = 0; myID < myNElts_; ++myID) { 2487 int baseID = myElts_[myID]; 2488 elemToMyElemMap_[baseID] = myID; 2489 } 2490 2491 // create myConnectivity_, myConnectivityUnique_ 2492 int numEltNodes = feElement_->num_elt_nodes(); 2493 myConnectivity_.reset(numEltNodes, myNElts_); 2494 myConnectivityUnique_.reset(numEltNodes, myNElts_); 2495 for (int elem = 0; elem < myNElts_; ++elem) { 2496 for (int node = 0; node < numEltNodes; ++node) { 2497 myConnectivity_(node,elem) = connectivity_(node,map_myElem_to_elem(elem)); 2498 myConnectivityUnique_(node,elem) = connectivityUnique_(node,map_myElem_to_elem(elem)); 2499 } 2500 } 2501 } 2502 2503 // ------------------------------------------------------------- 2504 // ------------------------------------------------------------- 2505 // class FE_Rectangular3DMesh 2506 // ------------------------------------------------------------- 2507 // ------------------------------------------------------------- 2508 FE_Rectangular3DMesh(const Array<double> & hx,const Array<double> & hy,const Array<double> & hz,const double xmin,const double xmax,const double ymin,const double ymax,const double zmin,const double zmax,const Array<bool> periodicity,const double xscale,const double yscale,const double zscale)2509 FE_Rectangular3DMesh::FE_Rectangular3DMesh( 2510 const Array<double> & hx, 2511 const Array<double> & hy, 2512 const Array<double> & hz, 2513 const double xmin, const double xmax, 2514 const double ymin, const double ymax, 2515 const double zmin, const double zmax, 2516 const Array<bool> periodicity, 2517 const double xscale, 2518 const double yscale, 2519 const double zscale) 2520 : hx_(hx), hy_(hy), hz_(hz) 2521 { 2522 tree_ = nullptr; 2523 hasPlanarFaces_ = true; 2524 xscale_ = xscale; 2525 yscale_ = yscale; 2526 zscale_ = zscale; 2527 2528 borders_[0][0] = xmin; 2529 borders_[0][1] = ymin; 2530 borders_[0][2] = zmin; 2531 borders_[1][0] = xmax; 2532 borders_[1][1] = ymax; 2533 borders_[1][2] = zmax; 2534 L_[0] = xmax-xmin; 2535 L_[1] = ymax-ymin; 2536 L_[2] = zmax-zmin; 2537 n_[0] = hx_.size(); 2538 n_[1] = hy_.size(); 2539 n_[2] = hz_.size(); 2540 // Compute region size and element size 2541 double Lx = 0; 2542 for (int i = 0; i < n_[0]; ++i) { Lx += hx_(i); } 2543 double Ly = 0; 2544 for (int i = 0; i < n_[1]; ++i) { Ly += hy_(i); } 2545 double Lz = 0; 2546 for (int i = 0; i < n_[2]; ++i) { Lz += hz_(i); } 2547 // rescale to fit box 2548 double ax = L_[0]/Lx; 2549 for (int i = 0; i < n_[0]; ++i) { hx_(i) *= ax; } 2550 double ay = L_[1]/Ly; 2551 for (int i = 0; i < n_[1]; ++i) { hy_(i) *= ay; } 2552 double az = L_[2]/Lz; 2553 for (int i = 0; i < n_[2]; ++i) { hz_(i) *= az; } 2554 2555 // fill node locations 2556 nSD_ = 3; 2557 x_.reserve(nSD_); 2558 for (int i = 0; i < nSD_; ++i) {x_.push_back(Array<double>(n_[i]+1)); } 2559 Array<double> & xI = x_[0]; 2560 xI(0) = xmin; 2561 for (int i = 0; i < n_[0]; ++i) { xI(i+1) = xI(i)+hx_(i); } 2562 Array<double> & yI = x_[1]; 2563 yI(0) = ymin; 2564 for (int i = 0; i < n_[1]; ++i) { yI(i+1) = yI(i)+hy_(i); } 2565 Array<double> & zI = x_[2]; 2566 zI(0) = zmin; 2567 for (int i = 0; i < n_[2]; ++i) { zI(i+1) = zI(i)+hz_(i); } 2568 2569 // Member data setup 2570 nElts_ = n_[0] * n_[1] * n_[2]; 2571 nNodes_ = (n_[0]+1) * (n_[1]+1) * (n_[2]+1); 2572 2573 periodicity_ = Array<bool>(periodicity); 2574 2575 feElement_ = new FE_ElementRect(); 2576 nodalCoords_.reset(3, nNodes_); 2577 connectivity_.reset(feElement_->num_elt_nodes(), nElts_); 2578 2579 // Fill nodal coordinates 2580 double x[3] = {xmin,ymin,zmin}; 2581 int inode = 0; 2582 for (int k = 0; k <= n_[2]; ++k) { 2583 for (int j = 0; j <= n_[1]; ++j) { 2584 for (int i = 0; i <= n_[0]; ++i) { 2585 for (int m = 0; m < 3; ++m) { 2586 nodalCoords_(m,inode) = x[m]; 2587 } 2588 ++inode; 2589 if (i < n_[0]) x[0] += hx_(i); 2590 } 2591 if (j < n_[1]) x[1] += hy_(j); 2592 x[0] = xmin; 2593 } 2594 if (k < n_[2]) x[2] += hz_(k); 2595 x[1] = ymin; 2596 } 2597 2598 // Compute element connectivities 2599 int ielt = 0; 2600 int noffx = 1; 2601 int noffy = n_[0] + 1; 2602 int noffz = (n_[0]+1) * (n_[1]+1); 2603 for (int k = 0; k < n_[2]; ++k) { 2604 for (int j = 0; j < n_[1]; ++j) { 2605 for (int i = 0; i < n_[0]; ++i) { 2606 int i1 = i + j*noffy + k*noffz; 2607 connectivity_(0,ielt) = i1; 2608 connectivity_(1,ielt) = i1 + noffx; 2609 connectivity_(2,ielt) = i1 + noffx + noffy; 2610 connectivity_(3,ielt) = i1 + noffy; 2611 connectivity_(4,ielt) = i1 + noffz; 2612 connectivity_(5,ielt) = i1 + noffx + noffz; 2613 connectivity_(6,ielt) = i1 + noffx + noffy + noffz; 2614 connectivity_(7,ielt) = i1 + noffy + noffz; 2615 ++ielt; 2616 } 2617 } 2618 } 2619 setup_periodicity(); 2620 } 2621 2622 // ------------------------------------------------------------- 2623 // partition_mesh 2624 // ------------------------------------------------------------- 2625 partition_mesh()2626 void FE_Rectangular3DMesh::partition_mesh() 2627 { 2628 if (lammpsPartition_) { 2629 lammps_partition_mesh(); 2630 } 2631 2632 if (partitioned_) return; 2633 2634 2635 // Currently this code has been rather naively copied from 2636 // FE_Uniform3DMesh::partition_mesh() 2637 2638 // Determine dimensions of mesh in order to partition according to largest dimension. 2639 double xmin = borders_[0][0]; 2640 double xmax = borders_[1][0]; 2641 double ymin = borders_[0][1]; 2642 double ymax = borders_[1][1]; 2643 double zmin = borders_[0][2]; 2644 double zmax = borders_[1][2]; 2645 2646 int numProcs; 2647 MPI_Comm_size(MPI_COMM_WORLD, &numProcs); 2648 2649 int processorRank; 2650 MPI_Comm_rank(MPI_COMM_WORLD, &processorRank); 2651 2652 // Spatially partition along the largest dimension. 2653 procs_.clear(); 2654 if (max(max(L_[0], L_[1]), L_[2]) == L_[0]) { 2655 partitionAxis_ = 0; 2656 for (int i = 0; i < numProcs; ++i) { 2657 procs_.push_back(xmin + (L_[0]*i)/numProcs); 2658 } 2659 procs_.push_back(xmax); 2660 } 2661 else if (max(max(L_[0], L_[1]), L_[2]) == L_[1]) { 2662 partitionAxis_ = 1; 2663 for (int i = 0; i < numProcs; ++i) { 2664 procs_.push_back(ymin + (L_[1]*i)/numProcs); 2665 } 2666 procs_.push_back(ymax); 2667 } 2668 else { 2669 partitionAxis_ = 2; 2670 for (int i = 0; i < numProcs; ++i) { 2671 procs_.push_back(zmin + (L_[2]*i)/numProcs); 2672 } 2673 procs_.push_back(zmax); 2674 } 2675 2676 // Distribute each node to the correct processor 2677 myNodes_.clear(); 2678 for (int i = 0; i < nNodes_; ++i) { 2679 // Allocate node to this processor if it lies between processor's left and right boundaries. 2680 if ( dbl_geq(nodalCoords_(partitionAxis_, i), procs_[processorRank]) && 2681 !dbl_geq(nodalCoords_(partitionAxis_, i), procs_[processorRank + 1])) { 2682 myNodes_.push_back(i); 2683 } 2684 // Also allocate nodes on the right boundary to the last processor. 2685 if ((processorRank == numProcs - 1) && 2686 dbl_geq(nodalCoords_(partitionAxis_, i), procs_[processorRank + 1])) { 2687 myNodes_.push_back(i); 2688 } 2689 } 2690 2691 // Distribute each element to the correct processor - assign it to the processor 2692 // which owns its node of lowest index. (this partitioning scheme is unambiguous) 2693 myElts_.clear(); 2694 for (int i = 0; i < nElts_; ++i) { 2695 int min = INT_MAX; 2696 for (int j = 0; j < connectivity_.nRows(); j++) { 2697 if (connectivity_(j, i) < min) 2698 min = connectivity_(j, i); 2699 } 2700 if (find(myNodes_.begin(), myNodes_.end(), min) != myNodes_.end()) { 2701 myElts_.push_back(i); 2702 } 2703 } 2704 2705 /* Commented out because ghost nodes are never used and dx_ is not a member of 2706 FE_Rectangular3DMesh. 2707 2708 // Compute the facesets that describes the left and right boundaries 2709 // in order to determine ghost nodes. 2710 int leftMult = 0; 2711 while ((leftMult+1)*dx_[partitionAxis_] < procs_[processorRank]) { 2712 ++leftMult; 2713 } 2714 int rightMult = 0; 2715 while ((rightMult)*dx_[partitionAxis_] < procs_[processorRank+1]) { 2716 ++rightMult; 2717 } 2718 // Compute our ghost nodes - nodes that we need that belong to adjacent processors, 2719 // and our shared nodes - our nodes that are ghosted on the adjacent processors. 2720 for (int i = 0; i < nNodes_; ++i) { 2721 if (nodalCoords_(partitionAxis_, i) == leftMult*dx_[partitionAxis_]) 2722 ghostNodesL_.push_back(i); 2723 else if (nodalCoords_(partitionAxis_, i) == rightMult*dx_[partitionAxis_]) 2724 ghostNodesR_.push_back(i); 2725 else if (nodalCoords_(partitionAxis_, i) == (leftMult+1)*dx_[partitionAxis_]) 2726 shareNodesL_.push_back(i); 2727 else if (nodalCoords_(partitionAxis_, i) == (rightMult-1)*dx_[partitionAxis_]) 2728 shareNodesR_.push_back(i); 2729 }*/ 2730 if (decomposition_) distribute_mesh_data(); 2731 partitioned_ = true; 2732 } 2733 departition_mesh()2734 void FE_Rectangular3DMesh::departition_mesh() 2735 { 2736 if (!partitioned_) return; 2737 partitioned_ = false; 2738 } 2739 2740 // ------------------------------------------------------------- 2741 // setup_periodicity 2742 // ------------------------------------------------------------- setup_periodicity()2743 void FE_Rectangular3DMesh::setup_periodicity() 2744 { 2745 nNodesUnique_ = 1; 2746 for (int i = 0; i < 3; i++) { 2747 nNodesUnique_ *= (n_[i] + 1 - periodicity_(i)); 2748 } 2749 2750 // form maximal nodeset 2751 for (int i = 0; i < nNodesUnique_; i++) { 2752 nodeSetAll_.insert(i); 2753 } 2754 2755 // Create global-to-unique map: globalToUniqueMap_(ig) = iu 2756 globalToUniqueMap_.reset(nNodes_); 2757 uniqueToGlobalMap_.reset(nNodesUnique_); 2758 for (int k = 0; k <= n_[2]; ++k) { 2759 int kper = (k == n_[2] && periodicity_(2)) ? 0 : k; 2760 for (int j = 0; j <= n_[1]; ++j) { 2761 int jper = (j == n_[1] && periodicity_(1)) ? 0 : j; 2762 for (int i = 0; i <= n_[0]; ++i) { 2763 int iper = (i == n_[0] && periodicity_(0)) ? 0 : i; 2764 int id = i + j*(n_[0]+1) + k*(n_[0]+1)*(n_[1]+1); 2765 int uid = iper + jper*(n_[0]+1-periodicity_(0)) 2766 + kper*(n_[0]+1-periodicity_(0))*(n_[1]+1-periodicity_(1)); 2767 globalToUniqueMap_(id) = uid; 2768 uniqueToGlobalMap_(uid) = id; 2769 } 2770 } 2771 } 2772 set_unique_connectivity(); 2773 2774 // form maximal elementset 2775 for (int i = 0; i < nElts_; i++) { 2776 elementSetAll_.insert(i); 2777 } 2778 } 2779 map_to_element(const DENS_VEC & x) const2780 int FE_Rectangular3DMesh::map_to_element(const DENS_VEC &x) const 2781 { 2782 int ix[3]; // ix[i] is the element in the ith direction 2783 for (int i = 0; i < 3; i++) { 2784 // map to box 2785 double y = x(i); 2786 if (periodicity_(i)) { 2787 double diff = y-borders_[0][i]; 2788 int shift = int(diff/L_[i]); 2789 if (diff < 0.) shift--; 2790 y -= shift*L_[i]; 2791 } 2792 // project into element 2793 ix[i] = x_[i].index(y); 2794 if (fabs(y-borders_[0][i]) < tol) { ix[i] = 0; } // on the lower boundary 2795 if (ix[i] < 0 || ix[i] >= n_[i]) { 2796 string msg = "FE_Rectangular3DMesh:: point maps outside of mesh, coordinate " 2797 + index_to_string(i) + "=" + to_string(x(i)) + " image=" + to_string(y) 2798 + " not in " + to_string(borders_[0][i]) + ":" + to_string(borders_[1][i]); 2799 throw ATC_Error(msg); 2800 } 2801 } 2802 int elt = ix[2]*(n_[0]*n_[1]) + ix[1]*n_[0] + ix[0]; 2803 return elt; 2804 } 2805 // ------------------------------------------------------------- 2806 // ------------------------------------------------------------- 2807 // class FE_Uniform3DMesh 2808 // ------------------------------------------------------------- 2809 // ------------------------------------------------------------- 2810 FE_Uniform3DMesh(const int nx,const int ny,const int nz,const double xmin,const double xmax,const double ymin,const double ymax,const double zmin,const double zmax,const Array<bool> periodicity,const double xscale,const double yscale,const double zscale)2811 FE_Uniform3DMesh::FE_Uniform3DMesh(const int nx, 2812 const int ny, 2813 const int nz, 2814 const double xmin, const double xmax, 2815 const double ymin, const double ymax, 2816 const double zmin, const double zmax, 2817 const Array<bool> periodicity, 2818 const double xscale, 2819 const double yscale, 2820 const double zscale) 2821 { 2822 hasPlanarFaces_ = true; 2823 tree_ = nullptr; 2824 xscale_ = xscale; 2825 yscale_ = yscale; 2826 zscale_ = zscale; 2827 n_[0] = nx; 2828 n_[1] = ny; 2829 n_[2] = nz; 2830 2831 borders_[0][0] = xmin; 2832 borders_[1][0] = xmax; 2833 borders_[0][1] = ymin; 2834 borders_[1][1] = ymax; 2835 borders_[0][2] = zmin; 2836 borders_[1][2] = zmax; 2837 2838 periodicity_ = Array<bool>(periodicity); 2839 2840 // Compute region size and element size 2841 for (int i = 0; i < 3; i++) { 2842 L_[i] = borders_[1][i] - borders_[0][i]; 2843 dx_[i] = L_[i]/n_[i]; 2844 } 2845 2846 // Member data setup 2847 nSD_ = 3; 2848 nElts_ = n_[0] * n_[1] * n_[2]; 2849 nNodes_ = (n_[0]+1) * (n_[1]+1) * (n_[2]+1); 2850 2851 feElement_ = new FE_ElementRect(); 2852 2853 nodalCoords_.reset(3, nNodes_); 2854 connectivity_.reset(feElement_->num_elt_nodes(), nElts_); 2855 2856 // Fill nodal coordinates 2857 double ix[3]; 2858 int inode = 0; 2859 for (int k = 0; k <= n_[2]; ++k) { 2860 ix[2] = borders_[0][2] + k*dx_[2]; 2861 for (int j = 0; j <= n_[1]; ++j) { 2862 ix[1] = borders_[0][1] + j*dx_[1]; 2863 for (int i = 0; i <= n_[0]; ++i) { 2864 ix[0] = borders_[0][0] + i*dx_[0]; 2865 for (int m = 0; m < 3; ++m) { 2866 nodalCoords_(m,inode) = ix[m]; 2867 } 2868 ++inode; 2869 } 2870 } 2871 } 2872 2873 // Compute element connectivities 2874 int ielt = 0; 2875 int noffx = 1; 2876 int noffy = n_[0] + 1; 2877 int noffz = (n_[0]+1) * (n_[1]+1); 2878 for (int k = 0; k < n_[2]; ++k) { 2879 for (int j = 0; j < n_[1]; ++j) { 2880 for (int i = 0; i < n_[0]; ++i) { 2881 int i1 = i + j*noffy + k*noffz; 2882 connectivity_(0,ielt) = i1; 2883 connectivity_(1,ielt) = i1 + noffx; 2884 connectivity_(2,ielt) = i1 + noffx + noffy; 2885 connectivity_(3,ielt) = i1 + noffy; 2886 connectivity_(4,ielt) = i1 + noffz; 2887 connectivity_(5,ielt) = i1 + noffx + noffz; 2888 connectivity_(6,ielt) = i1 + noffx + noffy + noffz; 2889 connectivity_(7,ielt) = i1 + noffy + noffz; 2890 ++ielt; 2891 } 2892 } 2893 } 2894 2895 setup_periodicity(); 2896 } 2897 2898 // ------------------------------------------------------------- 2899 // destructor 2900 // ------------------------------------------------------------- ~FE_Uniform3DMesh()2901 FE_Uniform3DMesh::~FE_Uniform3DMesh() 2902 { 2903 // Clean up is currently unimplemented 2904 } 2905 2906 // ------------------------------------------------------------- 2907 // partition_mesh 2908 // ------------------------------------------------------------- 2909 partition_mesh()2910 void FE_Uniform3DMesh::partition_mesh() 2911 { 2912 if (lammpsPartition_) { 2913 lammps_partition_mesh(); 2914 } 2915 2916 if (partitioned_) return; 2917 2918 2919 // Determine dimensions of mesh in order to partition according to largest dimension. 2920 double xmin = borders_[0][0]; 2921 double xmax = borders_[1][0]; 2922 double ymin = borders_[0][1]; 2923 double ymax = borders_[1][1]; 2924 double zmin = borders_[0][2]; 2925 double zmax = borders_[1][2]; 2926 2927 int numProcs; 2928 MPI_Comm_size(MPI_COMM_WORLD, &numProcs); 2929 2930 int processorRank; 2931 MPI_Comm_rank(MPI_COMM_WORLD, &processorRank); 2932 2933 // Spatially partition along the largest dimension. 2934 procs_.clear(); 2935 if (max(max(L_[0], L_[1]), L_[2]) == L_[0]) { 2936 partitionAxis_ = 0; 2937 for (int i = 0; i < numProcs; ++i) { 2938 procs_.push_back(xmin + (L_[0]*i)/numProcs); 2939 } 2940 procs_.push_back(xmax); 2941 } 2942 else if (max(max(L_[0], L_[1]), L_[2]) == L_[1]) { 2943 partitionAxis_ = 1; 2944 for (int i = 0; i < numProcs; ++i) { 2945 procs_.push_back(ymin + (L_[1]*i)/numProcs); 2946 } 2947 procs_.push_back(ymax); 2948 } 2949 else { 2950 partitionAxis_ = 2; 2951 for (int i = 0; i < numProcs; ++i) { 2952 procs_.push_back(zmin + (L_[2]*i)/numProcs); 2953 } 2954 procs_.push_back(zmax); 2955 } 2956 2957 // Distribute each node to the correct processor 2958 myNodes_.clear(); 2959 for (int i = 0; i < nNodes_; ++i) { 2960 // Allocate node to this processor if it lies between processor's left and right boundaries. 2961 if ( dbl_geq(nodalCoords_(partitionAxis_, i), procs_[processorRank]) && 2962 !dbl_geq(nodalCoords_(partitionAxis_, i), procs_[processorRank + 1])) { 2963 myNodes_.push_back(i); 2964 } 2965 // Also allocate nodes on the right boundary to the last processor. 2966 if ((processorRank == numProcs - 1) && 2967 dbl_geq(nodalCoords_(partitionAxis_, i), procs_[processorRank + 1])) { 2968 myNodes_.push_back(i); 2969 } 2970 } 2971 2972 // Distribute each element to the correct processor - assign it to the processor 2973 // which owns its node of lowest index. (this partitioning scheme is unambiguous) 2974 myElts_.clear(); 2975 for (int i = 0; i < nElts_; ++i) { 2976 int min = INT_MAX; 2977 for (int j = 0; j < connectivity_.nRows(); j++) { 2978 if (connectivity_(j, i) < min) 2979 min = connectivity_(j, i); 2980 } 2981 if (find(myNodes_.begin(), myNodes_.end(), min) != myNodes_.end()) { 2982 myElts_.push_back(i); 2983 } 2984 } 2985 2986 // Compute the facesets that describes the left and right boundaries 2987 // in order to determine ghost nodes. 2988 int leftMult = 0; 2989 while ((leftMult+1)*dx_[partitionAxis_] < procs_[processorRank]) { 2990 ++leftMult; 2991 } 2992 int rightMult = 0; 2993 while ((rightMult)*dx_[partitionAxis_] < procs_[processorRank+1]) { 2994 ++rightMult; 2995 } 2996 // Compute our ghost nodes - nodes that we need that belong to adjacent processors, 2997 // and our shared nodes - our nodes that are ghosted on the adjacent processors. 2998 for (int i = 0; i < nNodes_; ++i) { 2999 if (nodalCoords_(partitionAxis_, i) == leftMult*dx_[partitionAxis_]) 3000 ghostNodesL_.push_back(i); 3001 else if (nodalCoords_(partitionAxis_, i) == rightMult*dx_[partitionAxis_]) 3002 ghostNodesR_.push_back(i); 3003 else if (nodalCoords_(partitionAxis_, i) == (leftMult+1)*dx_[partitionAxis_]) 3004 shareNodesL_.push_back(i); 3005 else if (nodalCoords_(partitionAxis_, i) == (rightMult-1)*dx_[partitionAxis_]) 3006 shareNodesR_.push_back(i); 3007 } 3008 if (decomposition_) distribute_mesh_data(); 3009 partitioned_ = true; 3010 } 3011 departition_mesh()3012 void FE_Uniform3DMesh::departition_mesh() 3013 { 3014 if (!partitioned_) return; 3015 partitioned_ = false; 3016 } 3017 3018 // ------------------------------------------------------------- 3019 // map_to_element 3020 // ------------------------------------------------------------- map_to_element(const DENS_VEC & x) const3021 int FE_Uniform3DMesh::map_to_element(const DENS_VEC &x) const 3022 { 3023 // countx[i] is the number of the element, where 1 is the 3024 // element adjacent to the lower border, in the ith direction 3025 int countx[3]; 3026 for (int i = 0; i < 3; ++i) { 3027 // handle points on upper boundary; not sure why this is 3028 // hard-coded in, though... 3029 if (fabs(x(i)-borders_[1][i]) < tol) { 3030 countx[i] = n_[i] - 1; 3031 } else { 3032 // find the x, y, and z bins for the point in the mesh 3033 countx[i] = (int)floor((x(i)-borders_[0][i])/dx_[i]); 3034 } 3035 3036 // handle points out-of-range [0:nx-1] w/ periodicity 3037 if (countx[i] < 0 || countx[i] >= n_[i]) { 3038 if (periodicity_(i)) { 3039 countx[i] = countx[i] % n_[i]; 3040 // handles c++ ambiguous mod problems 3041 if (countx[i] < 0) countx[i] += n_[i]; 3042 } else { 3043 string msg = " point maps outside " 3044 "of mesh, coordinate " + 3045 index_to_string(i) + " " + to_string(x(i)) + 3046 " not in " + to_string(borders_[0][i]) + 3047 " : " + to_string(borders_[1][i]); 3048 throw ATC_Error(FILELINE,msg); 3049 } 3050 } 3051 } 3052 3053 int elt = countx[2]*(n_[0]*n_[1]) + 3054 countx[1]*n_[0] + 3055 countx[0]; 3056 return elt; 3057 } 3058 3059 3060 } // namespace ATC 3061