1 #ifndef FE_MESH_H 2 #define FE_MESH_H 3 4 #include "Array.h" 5 #include "Array2D.h" 6 #include "MatrixLibrary.h" 7 #include "ATC_TypeDefs.h" 8 #include "KD_Tree.h" 9 #include <vector> 10 #include <deque> 11 #include <list> 12 #include <map> 13 #include <set> 14 #include <utility> 15 #include <float.h> 16 #include <string> 17 #include <vector> 18 #include "mpi.h" 19 20 namespace ATC { 21 22 class FE_Element; 23 24 /** 25 * @class FE_Mesh 26 * @brief Base class for a finite element mesh 27 */ 28 class FE_Mesh { 29 30 public: 31 32 /** constructor */ 33 FE_Mesh(); 34 35 /** destructor */ 36 virtual ~FE_Mesh(); 37 38 /** parser/modifier */ 39 bool modify(int narg, char **arg); 40 41 /** initialization */ 42 void initialize(void); 43 44 /** write an unstructured mesh */ 45 void write_mesh(std::string meshFile); 46 47 48 49 // SJL why? will they ever be overridden by derived classes? in what 50 // situation would that be required, or desirable? virtual functions 51 // are slightly less efficient because they cannot be hard-linked at 52 // compile time. 53 is_partitioned()54 bool is_partitioned() const { return partitioned_; } 55 virtual void partition_mesh() = 0; 56 virtual void departition_mesh() = 0; 57 58 int map_elem_to_myElem(int elemID) const; 59 int map_myElem_to_elem(int myElemID) const; set_lammps_partition(bool t)60 void set_lammps_partition(bool t) {lammpsPartition_ = t;} set_data_decomposition(bool t)61 void set_data_decomposition(bool t) {decomposition_ = t;} 62 63 /** set quadrature on element from within interpolate class */ 64 void set_quadrature(FeIntQuadrature type); 65 66 /** evaluate shape function at real coordinates */ 67 void position(const int elem, 68 const VECTOR &xi, 69 DENS_VEC &x) const; 70 71 /** evaluate shape function at real coordinates */ 72 void shape_functions(const VECTOR &x, 73 DENS_VEC &N, 74 Array<int> &nodeList) const; 75 76 /** evaluate shape function at real coordinates */ 77 void shape_functions(const VECTOR &x, 78 DENS_VEC &N, 79 Array<int> &nodeList, 80 const Array<bool> &) const; 81 82 /** evaluate shape function at real coordinates */ 83 void shape_functions(const DENS_VEC &x, 84 DENS_VEC &N, 85 DENS_MAT &dNdx, 86 Array<int> &nodeList) const; 87 88 /** evaluate shape function at real coordinates */ 89 void shape_functions(const VECTOR &x, 90 const int eltID, 91 DENS_VEC &N, 92 Array<int> &nodeList) const; 93 94 /** evaluate shape function at real coordinates */ 95 void shape_functions(const DENS_VEC &x, 96 const int eltID, 97 DENS_VEC &N, 98 DENS_MAT &dNdx, 99 Array<int> &nodeList) const; 100 101 /** evaluate shape function at real coordinates */ 102 void shape_function_derivatives(const DENS_VEC &x, 103 const int eltID, 104 DENS_MAT &dNdx, 105 Array<int> &nodeList) const; 106 107 /** evaluate shape functions for all ip's on an element */ 108 // N is numIPsInElement X numNodesInElement 109 void shape_function(const int eltID, 110 DENS_MAT &N, 111 DIAG_MAT &weights) const; 112 113 /** evaluate shape functions for all ip's on an element */ 114 // N is numIPsInElement X numNodesInElement 115 void shape_function(const int eltID, 116 DENS_MAT &N, 117 std::vector<DENS_MAT> &dN, 118 DIAG_MAT &weights) const; 119 120 /** evaluate shape functions for all ip's on a face */ 121 // N is numIPsInFace X numNodesInElement 122 void face_shape_function(const PAIR &face, 123 DENS_MAT &N, 124 DENS_MAT &n, 125 DIAG_MAT &weights) const; 126 127 void face_shape_function(const PAIR &face, 128 DENS_MAT &N, 129 std::vector<DENS_MAT> &dN, 130 std::vector<DENS_MAT> &Nn, 131 DIAG_MAT &weights) const; 132 133 /** compute normal vector from the specified face */ 134 double face_normal(const PAIR &face, 135 const int ip, 136 DENS_VEC &normal) const; 137 138 139 /** return connectivity (global ID numbers) for element eltID */ 140 void element_connectivity_global(const int eltID, 141 Array<int> & nodes) const; 142 143 void element_connectivity_unique(const int eltID, 144 Array<int> & nodes) const; 145 146 int element_connectivity_global(const int eltID, 147 const int inode) const; 148 149 int element_connectivity_unique(const int eltID, 150 const int inode) const; 151 152 AliasArray<int> element_connectivity_global(const int eltID) const; 153 154 AliasArray<int> element_connectivity_unique(const int eltID) const; 155 face_connectivity(const PAIR & faceID,Array<int> & nodes)156 void face_connectivity(const PAIR & faceID, 157 Array<int> & nodes) const 158 { int nNodesPerFace = num_nodes_per_face(); 159 nodes.reset(nNodesPerFace); 160 int eltID = faceID.first; 161 int localFace = faceID.second; 162 const Array2D<int> & localFaceConn = local_face_connectivity(); 163 for(int i = 0; i < nNodesPerFace; ++i) { 164 nodes(i) = element_connectivity_global(eltID, localFaceConn(localFace,i)); 165 } 166 } face_connectivity_unique(const PAIR & faceID,Array<int> & nodes)167 void face_connectivity_unique(const PAIR & faceID, 168 Array<int> & nodes) const 169 { int nNodesPerFace = num_nodes_per_face(); 170 nodes.reset(nNodesPerFace); 171 int eltID = faceID.first; 172 int localFace = faceID.second; 173 const Array2D<int> & localFaceConn = local_face_connectivity(); 174 for(int i = 0; i < nNodesPerFace; ++i) { 175 nodes(i) = element_connectivity_unique(eltID, localFaceConn(localFace,i)); 176 } 177 } 178 179 /** 180 * return spatial coordinates for element nodes on eltID, 181 * indexed xCoords(isd,inode) 182 */ 183 void element_coordinates(const int eltID, 184 DENS_MAT & xCoords) const; 185 186 void face_coordinates(const PAIR face, 187 DENS_MAT & xCoords) const; 188 189 /** access to the nodal coordinate values */ nodal_coordinates(void)190 const DENS_MAT & nodal_coordinates(void) const {return nodalCoords_ ;} 191 192 /** access to nodal coordinates of a unique node */ 193 DENS_VEC nodal_coordinates(const int nodeID) const; 194 195 /** access to nodal coordinates of a node */ 196 DENS_VEC global_coordinates(const int nodeID) const; 197 198 /** map spatial location to element */ 199 virtual int map_to_element(const DENS_VEC &x) const = 0; 200 201 /** map global node numbering to unique node numbering */ map_global_to_unique(const int global_id)202 int map_global_to_unique(const int global_id) const 203 { 204 return globalToUniqueMap_(global_id); 205 } global_to_unique_map(void)206 inline const Array<int>& global_to_unique_map(void) const 207 { 208 return globalToUniqueMap_; 209 } 210 211 /** map unique node numbering a global node numbering */ map_unique_to_global(const int unique_id)212 int map_unique_to_global(const int unique_id) 213 { 214 return uniqueToGlobalMap_(unique_id); 215 } unique_to_global_map(void)216 inline const Array<int>& unique_to_global_map(void) const 217 { 218 return uniqueToGlobalMap_; 219 } 220 221 /** query whether a nodeset with the given name exists */ 222 bool query_nodeset(const std::string & name) const; 223 224 /** get node set (unique ID's) from the string name assigned to the set */ 225 const std::set<int> & nodeset(const std::string & name) const; 226 227 /** create node set with tag "name" from nodes in given spatial range */ 228 void create_nodeset(const std::string & name, const std::set<int> & nodeset); 229 void create_nodeset(const std::string & name, 230 double xmin, double xmax, 231 double ymin, double ymax, 232 double zmin, double zmax); 233 234 /** add to node set with tag "name" from nodes in given spatial range */ 235 void add_to_nodeset(const std::string & name, 236 double xmin, double xmax, 237 double ymin, double ymax, 238 double zmin, double zmax); 239 240 /** get element set from the string name assigned to the set */ 241 const std::set<int> & elementset(const std::string & name) const; 242 243 /** create element set with tag "name" from nodes in given spatial range */ 244 void create_elementset(const std::string & name, 245 double xmin, double xmax, 246 double ymin, double ymax, 247 double zmin, double zmax); 248 249 250 /** get the minimal element set from a nodeset by name */ 251 void nodeset_to_minimal_elementset(const std::string &name, 252 std::set<int> &elemSet) const; 253 /** get the minimal element set from a set of nodes */ 254 void nodeset_to_minimal_elementset(const std::set<int> &nodeSet, 255 std::set<int> &elemSet) const; 256 /** get the maximal element set from a nodeset by name */ 257 void nodeset_to_maximal_elementset(const std::string &name, 258 std::set<int> &elemSet) const; 259 /** get the maximal element set from a set of nodes */ 260 void nodeset_to_maximal_elementset(const std::set<int> &nodeSet, 261 std::set<int> &elemSet) const; 262 263 /** get complement of element set by name */ 264 void elementset_complement(const std::string &name, 265 std::set<int> &elemSet) const; 266 void elementset_complement(const std::set<int> &elemSet, 267 std::set<int> &elemSetComplement) const; 268 269 /** get the node set from an element set by name */ 270 void elementset_to_minimal_nodeset(const std::string &name, 271 std::set<int> &nodeSet) const; 272 273 void elementset_to_nodeset(const std::string &name, 274 std::set<int> &nodeSet) const; 275 void elementset_to_nodeset(const std::set<int> &elemSet, 276 std::set<int> &nodeSet) const; 277 278 /** convert faceset to nodeset in _unique_ node numbering */ 279 void faceset_to_nodeset(const std::string &name, 280 std::set<int> &nodeSet) const; 281 void faceset_to_nodeset(const std::set<PAIR> &faceSet, 282 std::set<int> &nodeSet) const; 283 284 void faceset_to_nodeset_global(const std::string &name, 285 std::set<int> &nodeSet) const; 286 void faceset_to_nodeset_global(const std::set<PAIR> &faceSet, 287 std::set<int> &nodeSet) const; 288 289 /** get face set from the string name assigned to the set */ 290 const std::set< std::pair<int,int> > & faceset(const std::string & name) const; 291 292 /** create face set with tag "name" from faces aligned with box */ 293 void create_faceset(const std::string & name, 294 double xmin, double xmax, 295 double ymin, double ymax, 296 double zmin, double zmax, 297 bool outward); 298 /** create face set with tag "name" from faces aligned with plane */ 299 void create_faceset(const std::string & name, double x, int idir, int isgn, 300 int nIdx2=-1, double x2lo=0.0, double x2hi=0.0, 301 int nIdx3=-1, double x3lo=0.0, double x3hi=0.0); 302 303 /** cut mesh */ 304 virtual void cut_mesh(const std::set<PAIR> & faceSet, const std::set<int> & nodeSet) = 0; 305 306 /** delete elements */ 307 virtual void delete_elements(const std::set<int> & elementList) = 0; 308 309 /** return number of spatial dimensions */ num_spatial_dimensions()310 int num_spatial_dimensions() const { return nSD_; } 311 312 /** return total number of nodes */ num_nodes()313 int num_nodes() const { return nNodes_; } 314 315 /** return number of unique nodes */ num_nodes_unique()316 int num_nodes_unique() const { return nNodesUnique_; } 317 318 /** return number of elements */ num_elements()319 int num_elements() const { return nElts_; } 320 321 /** return number of elements partitioned to my processor */ my_num_elements()322 int my_num_elements() const { return myNElts_; } 323 324 /** return number of integration points per element */ 325 int num_ips_per_element() const; 326 327 /** return number of nodes per element */ 328 int num_nodes_per_element() const; 329 330 /** return number of faces per element */ 331 int num_faces_per_element() const; 332 333 /** return number of nodes per face */ 334 int num_nodes_per_face() const; 335 336 /** return number of integration points per face */ 337 int num_ips_per_face() const; 338 339 /** return a pointer to the connectivity. This function will only work 340 when mesh is not partitioned. */ connectivity(void)341 Array2D<int> * connectivity(void) { return &connectivity_; } 342 /** return a pointer to the connectivity */ coordinates(void)343 DENS_MAT * coordinates(void) { return &nodalCoords_;} 344 /** Engine nodeMap stuff */ node_map(void)345 Array<int> *node_map(void) { return &globalToUniqueMap_;} 346 347 348 /** return scale in x,y,z */ xscale()349 double xscale() const { return xscale_; } yscale()350 double yscale() const { return yscale_; } zscale()351 double zscale() const { return zscale_; } 352 353 /** local face connectivity */ 354 const Array2D<int> & local_face_connectivity() const; 355 356 /** element size in each direction */ 357 virtual void bounding_box(const int ielem, 358 DENS_VEC & xmin, DENS_VEC & xmax); 359 360 /** element size in each direction */ 361 virtual void element_size(const int ielem, 362 double & hx, double & hy, double & hz); 363 364 /** element size in each direction */ min_element_size(void)365 virtual double min_element_size(void) const {return 0.0 ;} 366 367 /** get nodal coordinates for a given element */ element_field(const int eltIdx,const DENS_MAT f,DENS_MAT & local_field)368 void element_field(const int eltIdx, const DENS_MAT f, 369 DENS_MAT &local_field) 370 { 371 int dof = f.nCols(); 372 Array<int> nodes; 373 element_connectivity_unique(eltIdx,nodes); 374 local_field.reset(num_nodes_per_element(), dof); 375 for (int i = 0; i < nodes.size(); i++) { 376 for (int j = 0; j < dof; j++) local_field(i,j) = f(nodes(i),j); 377 } 378 } 379 380 /** almost structured */ 381 bool is_aligned(void) const; 382 383 /** extruded */ is_two_dimensional(void)384 bool is_two_dimensional(void) const {return twoDimensional_;} 385 coordinate_tolerance(void)386 virtual double coordinate_tolerance(void) const {return 1.e-8;} 387 388 /** element type */ 389 std::string element_type(void) const ; 390 391 /** output mesh subsets */ 392 void output(std::string prefix) const; 393 394 /* Parallelization data members */ 395 396 /** return element vector for this processor */ owned_elts()397 const std::vector<int> & owned_elts() const { return myElts_; } owned_and_ghost_elts()398 const std::vector<int> & owned_and_ghost_elts() const { 399 return (decomposition_) ? myAndGhostElts_: myElts_; } 400 bool is_owned_elt(int elt) const; 401 402 protected: 403 404 /** will this mesh use data decomposition? */ 405 bool decomposition_; 406 407 /** should the mesh use the native lammps partitioning? */ 408 bool lammpsPartition_; 409 410 /** is the element/node data currently partitioned among processors? */ 411 bool partitioned_; 412 413 /** number of spatial dimensions */ 414 int nSD_; 415 416 /** number of elements */ 417 int nElts_; 418 /** number of elements partitioned to this processor */ 419 int myNElts_; 420 421 /** number of nodes */ 422 int nNodes_; 423 int nNodesUnique_; 424 425 /** periodicity flags */ 426 Array<bool> periodicity_; 427 428 /** element type for this mesh */ 429 FE_Element *feElement_; 430 431 /** Nodal coordinates: nodalCoords_(nsd, numnode) */ 432 DENS_MAT nodalCoords_; 433 434 /** Element connectivity: connectivity_(neltnode, nelt) */ 435 Array2D<int> connectivity_; 436 Array2D<int> myConnectivity_; 437 Array2D<int> connectivityUnique_; 438 Array2D<int> myConnectivityUnique_; 439 440 /** map from unique node id to associated elements for periodic meshes */ 441 /** this data structure is only ever accessed from an unpartitioned context */ 442 Array<std::vector<int> > uniqueNodeToElementMap_; 443 444 /** map of global to unique node ID's */ 445 Array<int> globalToUniqueMap_; 446 Array<int> compactRemap_; // for condensing unique numbering 447 448 /** map of unique to global node ID's */ 449 Array<int> uniqueToGlobalMap_; 450 451 /** map of string names to node sets */ 452 NODE_SET_MAP nodeSetMap_; 453 454 /** maximal nodeset */ 455 std::set<int> nodeSetAll_; 456 457 /** map of string names to node sets */ 458 FACE_SET_MAP faceSetMap_; 459 460 /** map of string names to element sets */ 461 ELEMENT_SET_MAP elementSetMap_; 462 463 /** maximal elementset */ 464 std::set<int> elementSetAll_; 465 466 /** length scaling used by lammps */ 467 double xscale_, yscale_, zscale_; 468 469 /** Processor demarcations */ 470 std::vector<double> procs_; 471 472 /** Dimension (x=0, y=1, or z=2) */ 473 int partitionAxis_; 474 475 /** List of nodes for this processor */ 476 std::vector<int> myNodes_; 477 478 /** List of elements for this processor */ 479 std::vector<int> myElts_; 480 std::vector<int> myAndGhostElts_; 481 482 /** maps between my IDs and the total IDs */ 483 std::map<int,int> elemToMyElemMap_; 484 485 /** Lists of ghost nodes/neighbor ghost nodes */ 486 std::vector<int> ghostNodesL_; 487 std::vector<int> ghostNodesR_; 488 std::vector<int> shareNodesL_; 489 std::vector<int> shareNodesR_; 490 491 /** extruded */ 492 bool twoDimensional_; 493 bool hasPlanarFaces_; 494 double coordTol_; 495 496 }; 497 498 /** 499 * @class FE_3DMesh 500 * @brief Derived class for an unstructured 3D mesh 501 */ 502 class FE_3DMesh : public FE_Mesh { 503 public: 504 /** constructor */ FE_3DMesh()505 FE_3DMesh(){}; 506 507 /** constructor for read-in mesh **/ 508 // can later be extended to take nodesets, elementsets, etc. 509 FE_3DMesh(const std::string elementType, 510 const int nNodes, 511 const int nElements, 512 const Array2D<int> *connectivity, 513 const DENS_MAT *nodalCoordinates, 514 const Array<bool> periodicity, 515 const Array<std::pair<std::string,std::set<int> > > *nodeSets); 516 517 /** destructor */ 518 virtual ~FE_3DMesh(); 519 520 void partition_mesh(void); 521 522 void departition_mesh(void); 523 524 void lammps_partition_mesh(void); 525 526 /** Removes duplicate elements that appear in more than one vector 527 within procEltLists. **/ 528 void prune_duplicate_elements(std::vector<std::vector<int> > &procEltLists, 529 int *eltToOwners); 530 531 /** Takes procEltLists, and if there are more than nProcs of them 532 it takes the extra elements and distributes them to other vectors 533 in procEltLists. */ 534 535 // processors if during pruning processors end up 536 // elementless. This is slightly complicated because of 537 // ghost nodes. 538 void redistribute_extra_proclists(std::vector<std::vector<int> > &procEltLists, 539 int *eltToOwners, int nProcs); 540 541 /** This takes in a dense matrix and a list of elements 542 and fills in a standard adjacency list (within the matrix) 543 for those elements. **/ 544 545 // the set intersection, which does redundant computations 546 // right now, and filling in the adjacencies for both elements 547 // simultaneously when two elements share a face. 548 void compute_face_adjacencies(const std::list<int> &elts, 549 DENS_MAT &faceAdjacencies); 550 551 /** Counts the number of nonempty vectors in a vector of vectors. **/ 552 int numNonempty(std::vector<std::vector<int> > &v); 553 554 /** In the partitioning, we want to sort vectors of integers by size, 555 and furthermore we want empty vectors to count as the "largest" 556 possible vector because they dont want to count in the minimum. **/ 557 struct vectorComparer { operatorvectorComparer558 bool operator() (std::vector<int> l, std::vector<int> r) { 559 if (l.empty()) 560 return false; 561 if (r.empty()) 562 return true; 563 return (l.size() < r.size()); 564 } 565 } vectorCompSize; 566 min_element_size(void)567 virtual double min_element_size(void) const {return minEltSize_; } coordinate_tolerance(void)568 virtual double coordinate_tolerance(void) const { 569 return 0.25*(this->min_element_size()); // loose 570 //return 0.5; 571 } 572 virtual void cut_mesh(const std::set<PAIR> &faceSet, 573 const std::set<int> &nodeSet); 574 575 virtual void delete_elements(const std::set<int> &elementSet); 576 577 /** map spatial location to element */ 578 virtual int map_to_element(const DENS_VEC &x) const; 579 580 /** sends out data to processors during partitioning */ 581 void distribute_mesh_data(); 582 protected: 583 /** create global-to-unique node mapping */ 584 virtual void setup_periodicity(double tol = 1.e-8); 585 void fix_periodicity (int idim); 586 int find_boundary_nodes(int idim, std::set<int> & nodes); 587 bool match_nodes(int idim, std::set<int> & top, std::set<int> & bot, 588 Array<int> & map); 589 void set_unique_connectivity(void); 590 bool orient(int idir); 591 592 double minEltSize_; 593 KD_Tree *tree_; 594 595 /** test if a specified element actually contains the given point */ 596 bool contains_point(const int eltID, const DENS_VEC & x) const; 597 598 private: 599 Array<std::vector<int> > nodeToParentElements_; 600 601 }; 602 603 /** 604 * @class FE_Rectangular3DMesh 605 * @brief Derived class for a structured mesh with 606 * variable element sizes in x, y, and z directions 607 */ 608 class FE_Rectangular3DMesh : public FE_3DMesh { 609 public: 610 /** constructor */ FE_Rectangular3DMesh()611 FE_Rectangular3DMesh(){}; 612 FE_Rectangular3DMesh( 613 const Array<double> & hx, 614 const Array<double> & hy, 615 const Array<double> & hz, 616 const double xmin, const double xmax, 617 const double ymin, const double ymax, 618 const double zmin, const double zmax, 619 const Array<bool> periodicity, 620 const double xscale=1, 621 const double yscale=1, 622 const double zscale=1); 623 624 /** destructor */ ~FE_Rectangular3DMesh()625 virtual ~FE_Rectangular3DMesh() {}; 626 627 void partition_mesh(void); 628 629 void departition_mesh(void); 630 631 /** map spatial location to element */ 632 virtual int map_to_element(const DENS_VEC &x) const; 633 634 protected: 635 636 /** Number of elements in each spatial direction */ 637 int n_[3]; 638 639 /** Bounds of region on which mesh is defined */ 640 double borders_[2][3]; 641 642 /** Region size in each direction */ 643 double L_[3]; 644 645 646 647 /** create global-to-unique node mapping */ 648 virtual void setup_periodicity(); // note no "tol" 649 650 private: // only used by this class 651 /** partitions in x,y,z */ 652 Array<double> hx_, hy_, hz_; 653 654 /** nodal locations */ 655 std::vector< Array<double> > x_; 656 }; 657 658 /** 659 * @class FE_Uniform3DMesh 660 * @brief Derived class for a uniform structured mesh with 661 * fixed element sizes in x, y, and z directions 662 */ 663 class FE_Uniform3DMesh : public FE_Rectangular3DMesh { 664 665 public: 666 667 /** constructor */ 668 FE_Uniform3DMesh(const int nx, 669 const int ny, 670 const int nz, 671 const double xmin, const double xmax, 672 const double ymin, const double ymax, 673 const double zmin, const double zmax, 674 const Array<bool> periodicity, 675 const double xscale=1, 676 const double yscale=1, 677 const double zscale=1); 678 679 /** destructor */ 680 virtual ~FE_Uniform3DMesh(); 681 682 void partition_mesh(void); 683 684 void departition_mesh(void); 685 element_size(const int ielem,double hx,double hy,double hz)686 virtual void element_size(const int ielem, 687 double hx, double hy, double hz) 688 { hx = L_[0]/n_[0]; hy = L_[1]/n_[1]; hz = L_[2]/n_[2]; } 689 min_element_size(void)690 virtual double min_element_size(void) const 691 { return std::min(L_[0]/n_[0], std::min(L_[1]/n_[1], L_[2]/n_[2])); } 692 693 /** map spatial location to element */ 694 virtual int map_to_element(const DENS_VEC &x) const; 695 696 private: // only used by this class 697 /** Element size in each direction */ 698 double dx_[3]; 699 700 }; 701 702 }; // namespace ATC_Transfer 703 704 #endif // FE_MESH_H 705