intersect_swept_face.h
Go to the documentation of this file.
1 /*
2  * This file is part of the Ristra portage project.
3  * Please see the license file at the root of this repository, or at:
4  * https://github.com/laristra/portage/blob/master/LICENSE
5  */
6 
7 #pragma once
8 
11 #include "wonton/support/Point.h"
12 #include "wonton/support/Polytope.h"
13 
14 #ifdef HAVE_TANGRAM
15  #include "tangram/driver/CellMatPoly.h"
16  #include "tangram/driver/driver.h"
17  #include "tangram/support/MatPoly.h"
18  #include "tangram/reconstruct/cutting_distance_solver.h"
19 #endif
20 
21 /* -------------------------------------------------------------------------- */
22 namespace Portage {
23 
36  template<
37  int dim,
38  Entity_kind on_what,
39  class SourceMesh,
40  class SourceState,
41  class TargetMesh,
42  template<class, int, class, class>
43  class InterfaceReconstructor = DummyInterfaceReconstructor,
44  class Matpoly_Splitter = void,
45  class Matpoly_Clipper = void
46  >
48 
49  // useful aliases
50 #ifdef HAVE_TANGRAM
51  using InterfaceReconstructorDriver = Tangram::Driver<
52  InterfaceReconstructor, dim, SourceMesh,
53  Matpoly_Splitter, Matpoly_Clipper>;
54 #endif
55 
56  public:
57 
62  IntersectSweptFace() = delete;
63 
64 #ifdef HAVE_TANGRAM
65 
75  IntersectSweptFace(SourceMesh const &source_mesh,
76  SourceState const &source_state,
77  TargetMesh const &target_mesh,
78  NumericTolerances_t num_tols,
79  std::shared_ptr<InterfaceReconstructorDriver> ir)
80  : source_mesh_(source_mesh),
81  target_mesh_(target_mesh),
82  source_state_(source_state),
83  num_tols_(num_tols),
84  interface_reconstructor(ir) {}
85 
86 #endif
87 
96  IntersectSweptFace(SourceMesh const &source_mesh,
97  SourceState const &source_state,
98  TargetMesh const &target_mesh,
99  NumericTolerances_t num_tols)
100  : source_mesh_(source_mesh),
101  target_mesh_(target_mesh),
102  source_state_(source_state),
103  num_tols_(num_tols) {}
104 
111  IntersectSweptFace& operator=(IntersectSweptFace const& other) = delete;
112 
117  ~IntersectSweptFace() = default;
118 
124  void set_material(int m) { material_id_ = m; }
125 
130  void toggle_displacement_check(bool enable) { displacement_check = enable; }
131 
132 #ifdef HAVE_TANGRAM
133 
144  std::vector<double> compute_face_group_moments(
145  int const cell_id,
146  int const face_group_id,
147  double const swept_volume) const {
148  // see specialization for cells
149  std::cerr << "Sorry: current entity type not supported." << std::endl;
150  return std::vector<double>();
151  }
152 #endif
153 
162  std::vector<Weights_t> operator()(int target_id,
163  std::vector<int> const& stencil) const {
164  // see specialization for cells
165  std::cerr << "Sorry: current entity type not supported." << std::endl;
166  return std::vector<Weights_t>();
167  }
168 
169  private:
170  SourceMesh const &source_mesh_;
171  TargetMesh const &target_mesh_;
172  SourceState const &source_state_;
173  int material_id_ = -1;
174  NumericTolerances_t num_tols_ {};
175  bool displacement_check = false;
176 #ifdef HAVE_TANGRAM
177  std::shared_ptr<InterfaceReconstructorDriver> interface_reconstructor;
178 #endif
179  }; // class IntersectSweptFace
180 
181 
192  template<
193  class SourceMesh, class SourceState, class TargetMesh,
194  template<class, int, class, class>
195  class InterfaceReconstructor,
196  class Matpoly_Splitter, class Matpoly_Clipper
197  >
198  class IntersectSweptFace<2, Entity_kind::CELL,
199  SourceMesh, SourceState,
200  TargetMesh, InterfaceReconstructor,
201  Matpoly_Splitter, Matpoly_Clipper> {
202 
203  // useful aliases
204 #ifdef HAVE_TANGRAM
205  using InterfaceReconstructor2D = Tangram::Driver<
206  InterfaceReconstructor, 2, SourceMesh,
207  Matpoly_Splitter, Matpoly_Clipper>;
208 #endif
209 
210  public:
211 
216  IntersectSweptFace() = delete;
217 
218 #ifdef HAVE_TANGRAM
219 
229  IntersectSweptFace(SourceMesh const &source_mesh,
230  SourceState const &source_state,
231  TargetMesh const &target_mesh,
232  NumericTolerances_t num_tols,
233  std::shared_ptr<InterfaceReconstructor2D> ir)
234  : source_mesh_(source_mesh),
235  target_mesh_(target_mesh),
236  source_state_(source_state),
237  num_tols_(num_tols),
238  interface_reconstructor(ir) {}
239 
240 #endif
241 
250  IntersectSweptFace(SourceMesh const &source_mesh,
251  SourceState const &source_state,
252  TargetMesh const &target_mesh,
253  NumericTolerances_t num_tols)
254  : source_mesh_(source_mesh),
255  target_mesh_(target_mesh),
256  source_state_(source_state),
257  num_tols_(num_tols) {}
258 
265  IntersectSweptFace& operator=(IntersectSweptFace const& other) = delete;
266 
271  ~IntersectSweptFace() = default;
272 
278  void set_material(int m) { material_id_ = m; }
279 
280  private:
281 
296  std::vector<double> compute_moments_divergence_theorem
297  (std::vector<Wonton::Point<2>> const& swept_polygon) const {
298 
299  int const n = swept_polygon.size();
300  constexpr double const factor[] = { 0.5, 1./6., 1./6. }; // evaluated at compile-time
301 
302  std::vector<double> moments(3);
303 
304  for (int i=0; i < n; ++i) {
305  auto const& pi = swept_polygon[i];
306  auto const& pj = swept_polygon[(i + 1) % n];
307  auto const det = (pi[0] * pj[1] - pi[1] * pj[0]);
308  moments[0] += det;
309  moments[1] += (pi[0] + pj[0]) * det; // manually unrolled for perfs
310  moments[2] += (pi[1] + pj[1]) * det;
311  }
312 
313  for (int i=0; i < 3; ++i) { moments[i] *= factor[i]; }
314 
315  return moments;
316  }
317 
324  std::vector<double> compute_source_moments(int source_id) const {
325  double const area = source_mesh_.cell_volume(source_id);
326  Wonton::Point<2> centroid;
327  source_mesh_.cell_centroid(source_id, &centroid);
328  return std::vector<double>{ area, area * centroid[0], area * centroid[1] };
329  }
330 
338  bool centroid_inside_cell(int cell, std::vector<double> const& moment) const {
339 
340  double const cell_area = std::abs(moment[0]);
341  if (cell_area < num_tols_.min_absolute_volume)
342  return false;
343 
344  std::vector<int> edges, dirs, nodes;
345  source_mesh_.cell_get_faces_and_dirs(cell, &edges, &dirs);
346  int const nb_edges = edges.size();
347 
348  Wonton::Point<2> triangle[3];
349  triangle[0] = { moment[1] / cell_area, moment[2] / cell_area };
350 
351  for (int i = 0; i < nb_edges; ++i) {
352  source_mesh_.face_get_nodes(edges[i], &nodes);
353  // set triangle according to edge orientation.
354  if (dirs[i] > 0) {
355  source_mesh_.node_get_coordinates(nodes[0], triangle + 1);
356  source_mesh_.node_get_coordinates(nodes[1], triangle + 2);
357  } else {
358  source_mesh_.node_get_coordinates(nodes[1], triangle + 1);
359  source_mesh_.node_get_coordinates(nodes[0], triangle + 2);
360  }
361 
362  // check determinant sign
363  double const& ax = triangle[0][0];
364  double const& ay = triangle[0][1];
365  double const& bx = triangle[1][0];
366  double const& by = triangle[1][1];
367  double const& cx = triangle[2][0];
368  double const& cy = triangle[2][1];
369  double const det = ax * by - ax * cy
370  - bx * ay + bx * cy
371  + cx * ay - cx * by;
372  if (det < 0.)
373  return false;
374  }
375  return true;
376  }
377  public:
378 
383  void toggle_displacement_check(bool enable) { displacement_check = enable; }
384 
385 #ifdef HAVE_TANGRAM
386 
397  std::vector<double> compute_face_group_moments(
398  int const cell_id,
399  int const face_group_id,
400  double const swept_volume) const {
401 
402  std::vector<int> cfaces, cfdirs;
403  source_mesh_.cell_get_faces_and_dirs(cell_id, &cfaces, &cfdirs);
404 
405  int cface_id = std::distance(
406  cfaces.begin(), std::find(cfaces.begin(), cfaces.end(), face_group_id));
407 #ifdef DEBUG
408  //Face group should be associated with one of the cell's faces
409  int nfaces = cfaces.size();
410  assert(cface_id != nfaces);
411 #endif
412 
413  //Retrieve tolerance used by the interface reconstructor
414  const std::vector<Tangram::IterativeMethodTolerances_t>& ims_tols =
415  interface_reconstructor->iterative_methods_tolerances();
416  double dst_tol = ims_tols[0].arg_eps;
417  double vol_tol = ims_tols[0].fun_eps;
418 
419  //Create a MatPoly for the cell
420  Tangram::MatPoly<2> cell_mp;
421  Tangram::cell_get_matpoly(source_mesh_, cell_id, &cell_mp, dst_tol);
422  cell_mp.set_mat_id(0);
423  //Get the face normal and MatPoly's in the face's group
424  std::vector<Tangram::MatPoly<2>> face_group_polys;
425  Tangram::Plane_t<2> cutting_plane;
426  cutting_plane.normal =
427  cell_mp.face_normal_and_group(cface_id, face_group_id, &face_group_polys);
428 
429  //Find the cutting distance for the given swept volume
430  Tangram::CuttingDistanceSolver<2, Matpoly_Clipper> cds(face_group_polys,
431  cutting_plane.normal, ims_tols[0], true);
432 
433  cds.set_target_volume(std::fabs(swept_volume));
434  std::vector<double> cds_res = cds();
435 
436  //Check if we had enough volume in the face group
437  if (cds_res[1] < std::fabs(swept_volume) - vol_tol) {
438  throw std::runtime_error("Mesh displacement is too big for the implemented swept-face method");
439  }
440 
441  cutting_plane.dist2origin = cds_res[0];
442 
443  //Get the face group MatPoly's with material_id_ from the reconstructor
444  Tangram::CellMatPoly<2> const& cellmatpoly =
445  interface_reconstructor->cell_matpoly_data(cell_id);
446  std::vector<Tangram::MatPoly<2>> group_mat_polys =
447  cellmatpoly.get_face_group_matpolys(face_group_id, material_id_);
448 
449  //Clip the MatPoly's with the plane to get moments
450  Matpoly_Clipper clip_matpolys(vol_tol);
451  clip_matpolys.set_matpolys(group_mat_polys, true);
452  clip_matpolys.set_plane(cutting_plane);
453  std::vector<double> moments = clip_matpolys();
454 
455  //Weights need to be subtracted from a cell for negative swept regions
456  if(swept_volume < 0.0)
457  for(double& moment : moments)
458  moment *= -1;
459 
460  return moments;
461  }
462 #endif
463 
471  std::vector<Weights_t> operator()(int target_id,
472  std::vector<int> const& stencil) const {
473  // convenience function to check if a given cell is within the stencil.
474  auto in_stencil = [&](int cell) -> bool {
475  return std::find(stencil.begin(), stencil.end(), cell) != stencil.end();
476  };
477 
478  // here the source and target cell have the same ID.
479  int const source_id = target_id;
480 
481  std::vector<Weights_t> swept_moments;
482 
483  // Step 1: Add the source cell moments in the first place.
484  // For the single material case, add the moments of the source cell id.
485  // For the multimaterial case, add only the moments for the material
486  // the intersector is working on.
487 #ifdef HAVE_TANGRAM
488  int const nb_mats = source_state_.cell_get_num_mats(source_id);
489  std::vector<int> cellmats;
490  source_state_.cell_get_mats(source_id, &cellmats);
491  // nb_mats == 0 -- no materials ==> single material
492  // material_id_ == -1 -- intersect with mesh not a particular material
493  // nb_mats == 1 && cellmats[0] == material_id_ -- intersection with pure cell
494  // containing material_id_
495  bool const source_cell_mat =
496  (std::find(cellmats.begin(), cellmats.end(), material_id_) != cellmats.end());
497  bool const single_mat_src_cell = !nb_mats || (material_id_ == -1) ||
498  (nb_mats == 1 && source_cell_mat);
499 
500  if (single_mat_src_cell) {
501 #endif
502  // add source cell moments in the first place
503  swept_moments.emplace_back(source_id, compute_source_moments(source_id));
504 
505 #ifdef HAVE_TANGRAM
506  } else if (source_cell_mat) {
507  // mixed cell should contain this material
508  assert(interface_reconstructor != nullptr);
509 
510  // add source cell moments in the first place:
511  // obtain the aggregated moments of all MatPoly's with material_id_
512  Tangram::CellMatPoly<2> const& cellmatpoly =
513  interface_reconstructor->cell_matpoly_data(source_id);
514 
515  swept_moments.emplace_back(source_id, cellmatpoly.material_moments(material_id_));
516  }
517 #endif
518 
519  // Step 2: Obtain all facets and normals of the source cell
520  std::vector<int> edges, dirs, nodes;
521 
522  // retrieve current source cell faces/edges and related directions
523  source_mesh_.cell_get_faces_and_dirs(source_id, &edges, &dirs);
524  int const nb_edges = edges.size();
525 
526 #ifdef DEBUG
527  // ensure that we have the same face/edge index for source and target.
528  std::vector<int> target_edges, target_dirs, target_nodes;
529  target_mesh_.cell_get_faces_and_dirs(target_id, &target_edges, &target_dirs);
530  int const nb_target_edges = target_edges.size();
531 
532  assert(nb_edges == nb_target_edges);
533  for (int j = 0; j < nb_edges; ++j) {
534  assert(edges[j] == target_edges[j]);
535  }
536 #endif
537 
538  // Step 3: Main loop over all facets of the source cell
539  for (int i = 0; i < nb_edges; ++i) {
540 
541  // step 3a: retrieve nodes and reorder them according to edge direction
542  nodes.clear();
543  source_mesh_.face_get_nodes(edges[i], &nodes);
544 
545 #ifdef DEBUG
546  // ensure that we have the same nodal indices for source and target.
547  target_mesh_.face_get_nodes(target_edges[i], &target_nodes);
548  int const nb_source_nodes = nodes.size();
549  int const nb_target_nodes = target_nodes.size();
550 
551  assert(nb_source_nodes == nb_target_nodes);
552  for (int j = 0; j < nb_target_nodes; ++j) {
553  assert(nodes[j] == target_nodes[j]);
554  }
555 #endif
556 
557  // step 3b: construct the swept face polygon
558  std::vector<Wonton::Point<2>> swept_polygon(4);
559 
560  // if the edge has the same orientation as the cell, then reverse
561  // its nodes order such that we have a positive swept volume on
562  // outside and negative swept volume on inside.
563  // otherwise keep the same nodal order.
564  unsigned const j = (dirs[i] > 0 ? 1 : 0);
565  unsigned const k = j ^ 1;
566 
567  source_mesh_.node_get_coordinates(nodes[j], swept_polygon.data());
568  source_mesh_.node_get_coordinates(nodes[k], swept_polygon.data()+1);
569  target_mesh_.node_get_coordinates(nodes[k], swept_polygon.data()+2);
570  target_mesh_.node_get_coordinates(nodes[j], swept_polygon.data()+3);
571 
572  // step 3c: compute swept polygon moment using divergence theorem
573  auto moments = compute_moments_divergence_theorem(swept_polygon);
574 
575  // step 3d: add the source cells and their correct swept-moments to
576  // the weights vector. Approaches to compute the amount
577  // of swept-moment to be added are different for single and multimaterial cases.
578 
579  if (std::fabs(moments[0]) < num_tols_.min_absolute_volume) {
580  // skip if the swept polygon is almost empty.
581  // it may occur when the cell is shifted only in one direction.
582  continue;
583  } else if (moments[0] < 0.0) {
584  // if the computed swept face area is negative then assign its
585  // moments to the source cell: it will be substracted
586  // from the source cell area when performing the interpolation.
587 #ifdef HAVE_TANGRAM
588  if (single_mat_src_cell) {
589 #endif
590  swept_moments.emplace_back(source_id, moments);
591 #ifdef HAVE_TANGRAM
592  } else if (source_cell_mat) {
593  swept_moments.emplace_back(source_id,
594  compute_face_group_moments(source_id, edges[i], moments[0]));
595  }
596 #endif
597  } else {
598  // retrieve the cell incident to the current edge.
599  int const neigh = source_mesh_.cell_get_face_adj_cell(source_id, edges[i]);
600 
601  // just skip in case of a boundary edge
602  if (neigh < 0) {
603  continue;
604  }
605  // sanity check: ensure that incident cell belongs to the stencil.
606  if (!in_stencil(neigh)) {
607  auto id = std::to_string(source_id);
608  throw std::runtime_error("invalid stencil for source cell "+ id);
609  }
610  // sanity check: ensure that swept face centroid remains
611  // inside the neighbor cell and append to list.
612  if (displacement_check && !centroid_inside_cell(neigh, moments)) {
613  throw std::runtime_error("invalid target mesh for swept face");
614  }
615  // append to list as current neighbor moment.
616 #ifdef HAVE_TANGRAM
617  int const adj_cell_nb_mats = source_state_.cell_get_num_mats(neigh);
618  std::vector<int> adj_cellmats;
619  source_state_.cell_get_mats(neigh, &adj_cellmats);
620  // adj_cell_nb_mats == 0 -- no materials ==> single material
621  // material_id_ == -1 -- intersect with mesh not a particular material
622  // adj_cell_nb_mats == 1 && adj_cellmats[0] == material_id_ -- intersection with pure cell
623  // containing material_id_
624  bool const adj_cell_mat =
625  (std::find(adj_cellmats.begin(), adj_cellmats.end(), material_id_) != adj_cellmats.end());
626  bool const single_mat_adj_cell = !adj_cell_nb_mats || (material_id_ == -1) ||
627  (adj_cell_nb_mats == 1 && adj_cell_mat);
628 
629  if (single_mat_adj_cell) {
630 #endif
631  swept_moments.emplace_back(neigh, moments);
632 #ifdef HAVE_TANGRAM
633  } else {
634  //Skip if the neighboring cell doesn't contain material_id_
635  if (!adj_cell_mat)
636  continue;
637 
638  //Compute and append moments for the neighbor
639  swept_moments.emplace_back(neigh,
640  compute_face_group_moments(neigh, edges[i], moments[0]));
641  }
642 #endif
643  }
644  } // end for each edge of current cell
645 
646  return swept_moments;
647  }
648 
649  private:
650  SourceMesh const &source_mesh_;
651  TargetMesh const &target_mesh_;
652  SourceState const &source_state_;
653  int material_id_ = -1;
654  NumericTolerances_t num_tols_ {};
655  bool displacement_check = false;
656 #ifdef HAVE_TANGRAM
657  std::shared_ptr<InterfaceReconstructor2D> interface_reconstructor;
658 #endif
659  }; // class IntersectSweptFace::2D::CELL
660 
671  template<
672  class SourceMesh, class SourceState, class TargetMesh,
673  template<class, int, class, class>
674  class InterfaceReconstructor,
675  class Matpoly_Splitter, class Matpoly_Clipper
676  >
677  class IntersectSweptFace<3, Entity_kind::CELL,
678  SourceMesh, SourceState,
679  TargetMesh, InterfaceReconstructor,
680  Matpoly_Splitter, Matpoly_Clipper> {
681 
682  // useful aliases
683 #ifdef HAVE_TANGRAM
684  using InterfaceReconstructor3D = Tangram::Driver<
685  InterfaceReconstructor, 3, SourceMesh,
686  Matpoly_Splitter, Matpoly_Clipper
687  >;
688 #endif
689 
690  using Polyhedron = Wonton::Polytope<3>;
691 
692  public:
693 
698  IntersectSweptFace() = delete;
699 
700 #ifdef HAVE_TANGRAM
701 
711  IntersectSweptFace(SourceMesh const& source_mesh,
712  SourceState const& source_state,
713  TargetMesh const& target_mesh,
714  NumericTolerances_t num_tols,
715  std::shared_ptr<InterfaceReconstructor3D> ir)
716  : source_mesh_(source_mesh),
717  target_mesh_(target_mesh),
718  source_state_(source_state),
719  num_tols_(num_tols),
720  interface_reconstructor(ir) {}
721 
722 #endif
723 
732  IntersectSweptFace(SourceMesh const& source_mesh,
733  SourceState const& source_state,
734  TargetMesh const& target_mesh,
735  NumericTolerances_t num_tols)
736  : source_mesh_(source_mesh),
737  target_mesh_(target_mesh),
738  source_state_(source_state),
739  num_tols_(num_tols) {}
740 
747  IntersectSweptFace& operator=(IntersectSweptFace const& other) = delete;
748 
753  ~IntersectSweptFace() = default;
754 
760  void set_material(int m) { material_id_ = m; }
761 
762  private:
763 
787  std::vector<double> compute_moments
788  (std::vector<Wonton::Point<3>> const& coords,
789  std::vector<std::vector<int>> const& faces) const {
790 
791  // implementation may change in the future
792  Polyhedron polyhedron(coords, faces);
793  return polyhedron.moments();
794  }
795 
802  std::vector<double> compute_source_moments(int source_id) const {
803 
804  Wonton::Point<3> centroid;
805  double const volume = source_mesh_.cell_volume(source_id);
806  source_mesh_.cell_centroid(source_id, &centroid);
807  return std::vector<double>{volume,
808  centroid[0] * volume,
809  centroid[1] * volume,
810  centroid[2] * volume};
811  }
812 
836  bool valid_displacement(int source_id,
837  std::vector<Weights_t> const& moments) const {
838 
839  int const nb_moments = moments.size();
840  assert(nb_moments > 0);
841 
842  // the first entry of the swept region moment list should be
843  // the source cell moment, that is its own volume and weighted centroid.
844  assert(moments[0].entityID == source_id);
845  double const source_cell_volume = moments[0].weights[0];
846  double source_swept_volume = 0.;
847  assert(source_cell_volume > num_tols_.min_absolute_volume);
848 
849  for (int i = 1; i < nb_moments; ++i) {
850  double const swept_volume = std::abs(moments[i].weights[0]);
851  double const entity_volume = source_mesh_.cell_volume(moments[i].entityID);
852  assert(swept_volume > num_tols_.min_absolute_volume);
853  // each swept region volume should not exceed that of the attached cell.
854  if (swept_volume > entity_volume) {
855  return false;
856  } else if (moments[i].entityID == source_id) {
857  // accumulate the self-contribution for further comparisons.
858  source_swept_volume += swept_volume;
859  }
860  }
861 
862  // Normally, we should have p <= 3 swept regions associated with the source
863  // cell, with p the number of dimensions involved by the displacement.
864  // The volume of each swept region should not normally exceed half the
865  // volume of the source cell itself. Hence to ensure that the displacement
866  // is not too large, we just check that the total self-contribution
867  // is less than twice the source cell volume.
868  // Notice that this criterion could be refined later.
869  return source_swept_volume <= (1.5 * source_cell_volume) + num_tols_.min_absolute_volume;
870  }
871 
872 
873  public:
878  void toggle_displacement_check(bool enable) { displacement_check = enable; }
879 
880 #ifdef HAVE_TANGRAM
881 
892  std::vector<double> compute_face_group_moments(
893  int const cell_id,
894  int const face_group_id,
895  double const swept_volume) const {
896 
897  std::vector<int> cfaces, cfdirs;
898  source_mesh_.cell_get_faces_and_dirs(cell_id, &cfaces, &cfdirs);
899 
900  int cface_id = std::distance(
901  cfaces.begin(), std::find(cfaces.begin(), cfaces.end(), face_group_id));
902 #ifdef DEBUG
903  //Face group should be associated with one of the cell's faces
904  int nfaces = cfaces.size();
905  assert(cface_id != nfaces);
906 #endif
907 
908  //Retrieve tolerance used by the interface reconstructor
909  const std::vector<Tangram::IterativeMethodTolerances_t>& ims_tols =
910  interface_reconstructor->iterative_methods_tolerances();
911  double dst_tol = ims_tols[0].arg_eps;
912  double vol_tol = ims_tols[0].fun_eps;
913 
914  //Create a MatPoly for the cell
915  Tangram::MatPoly<3> cell_mp;
916  Tangram::cell_get_matpoly(source_mesh_, cell_id, &cell_mp, dst_tol);
917  cell_mp.set_mat_id(0);
918  //Get the face normal and MatPoly's in the face's group
919  std::vector<Tangram::MatPoly<3>> face_group_polys;
920  Tangram::Plane_t<3> cutting_plane;
921  cutting_plane.normal =
922  cell_mp.face_normal_and_group(cface_id, face_group_id, &face_group_polys);
923 
924  //Find the cutting distance for the given swept volume
925  Tangram::CuttingDistanceSolver<3, Matpoly_Clipper> cds(face_group_polys,
926  cutting_plane.normal, ims_tols[0], true);
927 
928  cds.set_target_volume(std::fabs(swept_volume));
929  std::vector<double> cds_res = cds();
930 
931  //Check if we had enough volume in the face group
932  if (cds_res[1] < std::fabs(swept_volume) - vol_tol) {
933  throw std::runtime_error("Mesh displacement is too big for the implemented swept-face method");
934  }
935 
936  cutting_plane.dist2origin = cds_res[0];
937 
938  //Get the face group MatPoly's with material_id_ from the reconstructor
939  Tangram::CellMatPoly<3> const& cellmatpoly =
940  interface_reconstructor->cell_matpoly_data(cell_id);
941  std::vector<Tangram::MatPoly<3>> group_mat_polys =
942  cellmatpoly.get_face_group_matpolys(face_group_id, material_id_);
943 
944  //Clip the MatPoly's with the plane to get moments
945  Matpoly_Clipper clip_matpolys(vol_tol);
946  clip_matpolys.set_matpolys(group_mat_polys, true);
947  clip_matpolys.set_plane(cutting_plane);
948  std::vector<double> moments = clip_matpolys();
949 
950  if(swept_volume < 0.0)
951  for(double& moment : moments)
952  moment *= -1;
953 
954  return moments;
955  }
956 #endif
957 
978  std::vector<Weights_t> operator()(int target_id,
979  std::vector<int> const& stencil) const {
980 
981  // convenience lambda to check if some cell is within the stencil.
982  auto in_stencil = [&](int cell) -> bool {
983  return std::find(stencil.begin(), stencil.end(), cell) != stencil.end();
984  };
985 
986  // here the source and target cell have the exact ID.
987  int const source_id = target_id;
988 
989  std::vector<Weights_t> swept_moments;
990 
991 #ifdef HAVE_TANGRAM
992  int const nb_mats = source_state_.cell_get_num_mats(source_id);
993  std::vector<int> cellmats;
994  source_state_.cell_get_mats(source_id, &cellmats);
995  // nb_mats == 0 -- no materials ==> single material
996  // material_id_ == -1 -- intersect with mesh not a particular material
997  // nb_mats == 1 && cellmats[0] == material_id_ -- intersection with pure cell
998  // containing material_id_
999  bool const source_cell_mat =
1000  (std::find(cellmats.begin(), cellmats.end(), material_id_) != cellmats.end());
1001  bool const single_mat_src_cell = !nb_mats || (material_id_ == -1) ||
1002  (nb_mats == 1 && source_cell_mat);
1003 
1004  if (single_mat_src_cell) {
1005 #endif
1006  // add source cell moments in the first place
1007  swept_moments.emplace_back(source_id, compute_source_moments(source_id));
1008 
1009 #ifdef HAVE_TANGRAM
1010  } else if (source_cell_mat) {
1011  // mixed cell should contain this material
1012  assert(interface_reconstructor != nullptr);
1013 
1014  // add source cell moments in the first place:
1015  // obtain the aggregated moments of all MatPoly's with material_id_
1016  Tangram::CellMatPoly<3> const& cellmatpoly =
1017  interface_reconstructor->cell_matpoly_data(source_id);
1018 
1019  swept_moments.emplace_back(source_id, cellmatpoly.material_moments(material_id_));
1020  }
1021 #endif
1022  std::vector<int> faces, dirs, nodes;
1023 
1024  // retrieve current source cell faces/edges and related directions
1025  source_mesh_.cell_get_faces_and_dirs(source_id, &faces, &dirs);
1026  int const nb_faces = faces.size();
1027 
1028 #if DEBUG
1029  // ensure that we have the same face index for source and target.
1030  std::vector<int> target_faces, target_dirs, target_nodes;
1031  target_mesh_.cell_get_faces_and_dirs(target_id, &target_faces, &target_dirs);
1032  int const nb_target_faces = target_faces.size();
1033 
1034  assert(nb_faces == nb_target_faces);
1035  for (int j = 0; j < nb_faces; ++j) {
1036  assert(faces[j] == target_faces[j]);
1037  }
1038 #endif
1039 
1040  for (int i = 0; i < nb_faces; ++i) {
1041  // step 0: retrieve nodes and reorder them according to face orientation
1042  nodes.clear();
1043  source_mesh_.face_get_nodes(faces[i], &nodes);
1044 
1045  int const nb_face_nodes = nodes.size();
1046  int const nb_poly_nodes = 2 * nb_face_nodes;
1047  int const nb_poly_faces = nb_poly_nodes + 2;
1048 
1049 #if DEBUG
1050  // ensure that we have the same nodal indices for source and target.
1051  target_mesh_.face_get_nodes(target_faces[i], &target_nodes);
1052  int const nb_source_nodes = nodes.size();
1053  int const nb_target_nodes = target_nodes.size();
1054 
1055  assert(nb_source_nodes == nb_target_nodes);
1056  for (int j = 0; j < nb_target_nodes; ++j) {
1057  assert(nodes[j] == target_nodes[j]);
1058  }
1059 #endif
1060 
1061  /* step 1: construct the swept volume polyhedron which can be:
1062  * - a prism for a triangular face,
1063  * - a hexahedron for a quadrilateral face,
1064  * - a (n+2)-face polyhedron for an arbitrary n-polygon.
1065  */
1066  std::vector<Wonton::Point<3>> swept_poly_coords(nb_poly_nodes);
1067  std::vector<std::vector<int>> swept_poly_faces(nb_poly_faces);
1068 
1069  /* the swept polyhedron must be formed in a way that the vertices of
1070  * each of its faces are ordered such that their normals outside that
1071  * swept volume - this is necessarily true except for the original face
1072  * inherited from the cell itself.
1073  * for this face, if the ordering of its vertices is such that the normal
1074  * points out of the cell, then the normal points into the swept volume.
1075  * in such a case, the vertex ordering must be reversed.
1076  *
1077  * source hex target hex face swept polyhedron:
1078  * 7'......6'
1079  * .: .: 4'......5'
1080  * 7______6 . : . : /: /:
1081  * /| /| 4'...:..5' : / : / :
1082  * / | / | : :..:..2' / : / :
1083  * 4 __|__5 | : . : . 4____:_5...:
1084  * | 3__|__2 : . : . | /0'| /1'
1085  * | / | / :......: | / | /
1086  * | / | / 0' 1' | / |/
1087  * |/_____|/ |/_____/
1088  * 0 1 0 1
1089  * ∙dirs[f] > 0: [4,5,1,0,4',5',1',0']
1090  * ∙dirs[f] < 0: [0,1,5,4,0',1',5',4']
1091  */
1092  bool const outward_normal = dirs[i] > 0;
1093 
1094  for (int current = 0; current < nb_face_nodes; ++current) {
1095  int const reverse = (nb_face_nodes - 1) - current;
1096  int const index = (outward_normal ? reverse : current);
1097  int const offset = current + nb_face_nodes;
1098  source_mesh_.node_get_coordinates(nodes[index], swept_poly_coords.data() + current);
1099  target_mesh_.node_get_coordinates(nodes[index], swept_poly_coords.data() + offset);
1100  }
1101 
1102  /* now build the swept polyhedron faces, which vertices are indexed
1103  * RELATIVELY to the polyhedron vertices list.
1104  * - first allocate memory for vertices list of each face.
1105  * - then add the original face and its twin induced by sweeping.
1106  * - eventually construct the other faces induced by edge sweeping.
1107  */
1108  for (int current = 0; current < nb_poly_faces; ++current) {
1109  // for each twin face induced by sweeping, its number of vertices is
1110  // exactly that of the current cell face, whereas the number of
1111  // vertices of the other faces is exactly 4.
1112  int const size = (current < 2 ? nb_face_nodes : 4);
1113  swept_poly_faces[current].resize(size);
1114  }
1115 
1116 
1117  /* swept polyhedron face construction rules:
1118  *
1119  * 3'_____2' n_poly_faces: 2 + n_face_edges = 2 + 4 = 6.
1120  * /| /| n_poly_nodes: 2 * n_face_edges = 2 * 4 = 8 = n.
1121  * / | / | ordered vertex list: [3,2,1,0,3',2',1',0']
1122  * / |__/__|
1123  * / / / / 1' absolute relative
1124  * 3___/_2/ / ∙f[0]: (3 |2 |1 |0 ) (0,1,2,3)
1125  * | / | / ∙f[1]: (3'|2'|1'|0') (4,5,6,7)
1126  * | / | / ∙f[2]: (3 |3'|2'|2 ) => (0,4,5,1)
1127  * |/____|/ ∙f[3]: (2 |2'|1'|1 ) (1,5,6,2)
1128  * 0 1 ∙f[4]: (1 |1'|0'|0 ) (2,6,7,3)
1129  * ∙f[5]: (0 |0'|3'|3 ) (3,7,4,0)
1130  *
1131  * let m = n/2 with n the number of polyhedron vertices.
1132  * - twin faces: [0, m-1] and [m-1, n].
1133  * - side faces: [i, i+m, ((i+1) % m)+m, (i+1) % m]
1134  */
1135  for (int current = 0; current < nb_face_nodes; ++current) {
1136  // a) set twin faces vertices
1137  swept_poly_faces[0][current] = current;
1138  swept_poly_faces[1][current] = nb_poly_nodes - current - 1;
1139 
1140  // b) set side faces vertices while keeping them counterclockwise.
1141  int const index = current + 2;
1142  swept_poly_faces[index][0] = current;
1143  swept_poly_faces[index][3] = (current + 1) % nb_face_nodes;
1144  swept_poly_faces[index][1] = swept_poly_faces[index][0] + nb_face_nodes;
1145  swept_poly_faces[index][2] = swept_poly_faces[index][3] + nb_face_nodes;
1146  }
1147 
1148  /* step 2: compute swept polygon moments using divergence theorem */
1149  auto moments = compute_moments(swept_poly_coords, swept_poly_faces);
1150 
1151  /* step 3: assign the computed moments to the source cell or one
1152  * of its neighbors according to the sign of the swept region volume.
1153  */
1154  if (std::abs(moments[0]) < num_tols_.min_absolute_volume) {
1155  // just skip if the swept region is almost flat.
1156  // it may occur when the cell is shifted only in one direction.
1157  continue;
1158  } else if (moments[0] < 0.) {
1159  // if the computed swept region volume is negative then assign its
1160  // moments to the source cell: it will be substracted
1161  // from the source cell area when performing the interpolation.
1162 #ifdef HAVE_TANGRAM
1163  if (single_mat_src_cell) {
1164 #endif
1165  swept_moments.emplace_back(source_id, moments);
1166 #ifdef HAVE_TANGRAM
1167  } else if (source_cell_mat) {
1168  swept_moments.emplace_back(source_id,
1169  compute_face_group_moments(source_id, faces[i], moments[0]));
1170  }
1171 #endif
1172  } else {
1173  // retrieve the cell incident to the current edge.
1174  int const neigh = source_mesh_.cell_get_face_adj_cell(source_id, faces[i]);
1175 
1176  // just skip in case of a boundary edge
1177  if (neigh < 0) {
1178  continue;
1179  }
1180  // sanity check: ensure that incident cell belongs to the stencil.
1181  if (!in_stencil(neigh)) {
1182  auto id = std::to_string(source_id);
1183  throw std::runtime_error("invalid stencil for source cell" + id);
1184  }
1185  // append to list as current neighbor moment.
1186 #ifdef HAVE_TANGRAM
1187  int const adj_cell_nb_mats = source_state_.cell_get_num_mats(neigh);
1188  std::vector<int> adj_cellmats;
1189  source_state_.cell_get_mats(neigh, &adj_cellmats);
1190  // adj_cell_nb_mats == 0 -- no materials ==> single material
1191  // material_id_ == -1 -- intersect with mesh not a particular material
1192  // adj_cell_nb_mats == 1 && adj_cellmats[0] == material_id_ -- intersection with pure cell
1193  // containing material_id_
1194  bool const adj_cell_mat =
1195  (std::find(adj_cellmats.begin(), adj_cellmats.end(), material_id_) != adj_cellmats.end());
1196  bool const single_mat_adj_cell = !adj_cell_nb_mats || (material_id_ == -1) ||
1197  (adj_cell_nb_mats == 1 && adj_cell_mat);
1198 
1199  if (single_mat_adj_cell) {
1200 #endif
1201  swept_moments.emplace_back(neigh, moments);
1202 #ifdef HAVE_TANGRAM
1203  } else {
1204  //Skip if the neighboring cell doesn't contain material_id_
1205  if (!adj_cell_mat)
1206  continue;
1207 
1208  //Compute and append moments for the neighbor
1209  swept_moments.emplace_back(neigh,
1210  compute_face_group_moments(neigh, faces[i], moments[0]));
1211  }
1212 #endif
1213  }
1214  } // end of for each face of current cell
1215 
1216  if (!displacement_check || valid_displacement(source_id, swept_moments)) {
1217  return swept_moments;
1218  } else
1219  throw std::runtime_error("invalid displacement");
1220  }
1221 
1222  private:
1223  SourceMesh const& source_mesh_;
1224  TargetMesh const& target_mesh_;
1225  SourceState const& source_state_;
1226  int material_id_ = -1;
1227  NumericTolerances_t num_tols_ {};
1228  bool displacement_check = false;
1229 #ifdef HAVE_TANGRAM
1230  std::shared_ptr<InterfaceReconstructor3D> interface_reconstructor;
1231 #endif
1232  }; // class IntersectSweptFace::3D::CELL
1233 
1234  /* Define aliases to have the same interface as the other intersectors such
1235  * as R2D and R3D. This simple workaround allow us to discard the dimension
1236  * template parameter when instantiating the kernel.
1237  */
1238  template<
1239  Entity_kind entity_kind,
1240  class SourceMesh, class SourceState, class TargetMesh,
1241  template<class, int, class, class>
1242  class InterfaceReconstructor = DummyInterfaceReconstructor,
1243  class Matpoly_Splitter = void, class Matpoly_Clipper = void
1244  >
1245  using IntersectSweptFace2D = IntersectSweptFace<2, entity_kind,
1246  SourceMesh, SourceState,
1247  TargetMesh,
1248  InterfaceReconstructor,
1249  Matpoly_Splitter,
1250  Matpoly_Clipper>;
1251 
1252  template<
1253  Entity_kind entity_kind,
1254  class SourceMesh, class SourceState, class TargetMesh,
1255  template<class, int, class, class>
1256  class InterfaceReconstructor = DummyInterfaceReconstructor,
1257  class Matpoly_Splitter = void, class Matpoly_Clipper = void
1258  >
1259  using IntersectSweptFace3D = IntersectSweptFace<3, entity_kind,
1260  SourceMesh, SourceState,
1261  TargetMesh,
1262  InterfaceReconstructor,
1263  Matpoly_Splitter,
1264  Matpoly_Clipper>;
1265 
1266  /* ------------------------------------------------------------------------ */
1267 } // namespace Portage
Definition: dummy_interface_reconstructor.h:30
void toggle_displacement_check(bool enable)
Toggle target mesh displacement validity check.
Definition: intersect_swept_face.h:130
std::vector< T > vector
Definition: portage.h:238
IntersectSweptFace & operator=(IntersectSweptFace const &other)=delete
Assignment operator (disabled).
Intersection and other tolerances to handle tiny values.
Definition: portage.h:152
Wonton::Point< OP::dim > centroid(const typename OP::points_t &apts)
Definition: operator.h:616
void set_material(int m)
Set the material we are operating on.
Definition: intersect_swept_face.h:124
~IntersectSweptFace()=default
Destructor.
IntersectSweptFace(SourceMesh const &source_mesh, SourceState const &source_state, TargetMesh const &target_mesh, NumericTolerances_t num_tols)
Constructor for single material case.
Definition: intersect_swept_face.h:250
void toggle_displacement_check(bool enable)
Toggle target mesh displacement validity check.
Definition: intersect_swept_face.h:383
Kernel to compute interpolation weights for swept-face remap.
Definition: intersect_swept_face.h:47
std::vector< Weights_t > operator()(int target_id, std::vector< int > const &stencil) const
Perform the actual swept faces volumes computation.
Definition: intersect_swept_face.h:162
void toggle_displacement_check(bool enable)
Toggle target mesh displacement validity check.
Definition: intersect_swept_face.h:878
Definition: coredriver.h:42
IntersectSweptFace(SourceMesh const &source_mesh, SourceState const &source_state, TargetMesh const &target_mesh, NumericTolerances_t num_tols)
Constructor for single material case.
Definition: intersect_swept_face.h:96
std::vector< Weights_t > operator()(int target_id, std::vector< int > const &stencil) const
Perform the actual swept faces computation.
Definition: intersect_swept_face.h:471
std::vector< Weights_t > operator()(int target_id, std::vector< int > const &stencil) const
Perform the actual swept moments computation for the given cell.
Definition: intersect_swept_face.h:978
IntersectSweptFace(SourceMesh const &source_mesh, SourceState const &source_state, TargetMesh const &target_mesh, NumericTolerances_t num_tols)
Constructor for single material case.
Definition: intersect_swept_face.h:732
std::string to_string(Limiter_type limiter_type)
Definition: portage.h:96
IntersectSweptFace()=delete
Default constructor (disabled).