MAST
Multidisciplinary-design Adaptation and Sensitivity Toolkit (MAST)
level_set_intersection.cpp
Go to the documentation of this file.
1 /*
2  * MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
3  * Copyright (C) 2013-2020 Manav Bhatia and MAST authors
4  *
5  * This library is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU Lesser General Public
7  * License as published by the Free Software Foundation; either
8  * version 2.1 of the License, or (at your option) any later version.
9  *
10  * This library is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13  * Lesser General Public License for more details.
14  *
15  * You should have received a copy of the GNU Lesser General Public
16  * License along with this library; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18  */
19 
20 // C++ includes
21 #include <memory>
22 #include <map>
23 #include <numeric>
24 
25 // MAST includes
28 
29 
30 
32 _tol (1.e-8),
33 _max_iters (10),
34 _max_mesh_elem_id (0),
35 _max_mesh_node_id (0),
36 _max_elem_divs (4),
37 _elem (nullptr),
38 _initialized (false),
39 _phi (nullptr),
40 _if_elem_on_positive_phi (false),
41 _if_elem_on_negative_phi (false),
42 _mode (MAST::NO_INTERSECTION),
43 _node_num_on_boundary (0),
44 _edge_num_on_boundary (0){
45 
46 }
47 
48 
50 
51  this->clear();
52 }
53 
54 
55 const libMesh::Elem&
57 
58  libmesh_assert(_initialized);
59 
60  return *_elem;
61 }
62 
63 
66 
67  libmesh_assert(_initialized);
68 
69  return _mode;
70 }
71 
72 
73 
74 unsigned int
76 
77  libmesh_assert(_initialized);
78  libmesh_assert(_mode == MAST::THROUGH_NODE);
79 
80  return _node_num_on_boundary;
81 }
82 
83 
84 
85 unsigned int
87 
88  libmesh_assert(_initialized);
89  libmesh_assert(_mode == MAST::COLINEAR_EDGE);
90 
91  return _edge_num_on_boundary;
92 }
93 
94 
95 
96 
97 bool
99 
100  libmesh_assert(_initialized);
101 
103 }
104 
105 
106 bool
108 
109  libmesh_assert(_initialized);
110 
112 }
113 
114 
115 
116 bool
118 
119  libmesh_assert(_initialized);
120 
121  return !_if_elem_on_positive_phi;
122 }
123 
124 
125 bool
127 
128  libmesh_assert(_initialized);
129 
130  return !_if_elem_on_negative_phi;
131 }
132 
133 
134 
135 bool
137 
138  return (this->if_intersection_through_elem() ||
140 }
141 
142 
143 
144 
145 bool
147 
148  return ((_mode == MAST::OPPOSITE_EDGES) ||
151  (_mode == MAST::NODE_AND_EDGE) ||
154 }
155 
156 
157 
158 Real
160 
161  libmesh_assert(_initialized);
162 
163  std::vector<const libMesh::Elem*>::const_iterator
164  it = _positive_phi_elems.begin(),
165  end = _positive_phi_elems.end();
166 
167  Real
168  v = 0.,
169  v0 = _elem->volume();
170 
171  for ( ; it != end; it++)
172  v += (*it)->volume();
173 
174  return v/v0;
175 }
176 
177 
178 
179 Real
180 MAST::LevelSetIntersection::get_node_phi_value(const libMesh::Node* n) const {
181 
182  libmesh_assert(_initialized);
183  std::map<const libMesh::Node*, std::pair<Real, bool> >::const_iterator
184  it = _node_phi_vals.find(n),
185  end = _node_phi_vals.end();
186 
187  libmesh_assert(it != end);
188  return it->second.first;
189 }
190 
191 
192 
193 bool
194 MAST::LevelSetIntersection::if_hanging_node(const libMesh::Node* n) const {
195 
196  return _hanging_node.count(n);
197 }
198 
199 
200 void
202 
203  _tol = 1.e-8;
204  _max_iters = 10;
205  _elem = nullptr;
206  _initialized = false;
207  _phi = nullptr;
208  _if_elem_on_positive_phi = false;
209  _if_elem_on_negative_phi = false;
211  _max_mesh_elem_id = 0;
212  _max_mesh_node_id = 0;
215  _positive_phi_elems.clear();
216  _negative_phi_elems.clear();
217  _elem_sides_on_interface.clear();
218  _node_local_coords.clear();
219  _hanging_node.clear();
220 
221  std::vector<libMesh::Elem*>::iterator
222  e_it = _new_elems.begin(),
223  e_end = _new_elems.end();
224 
225  for ( ; e_it != e_end; e_it++)
226  delete *e_it;
227 
228 
229  std::vector<libMesh::Node*>::iterator
230  n_it = _new_nodes.begin(),
231  n_end = _new_nodes.end();
232 
233  for ( ; n_it != n_end; n_it++)
234  delete *n_it;
235 
236  _new_nodes.clear();
237  _new_elems.clear();
238 
239  _interior_nodes.clear();
240  _bounding_nodes.clear();
241 }
242 
243 
244 void
246  const libMesh::Elem& e,
247  const Real t,
248  unsigned int max_elem_id,
249  unsigned int max_node_id) {
250 
251  // make sure that this has not been initialized already
252  libmesh_assert(!_initialized);
253  libmesh_assert_equal_to(e.dim(), 2); // this is only for 2D elements
254 
255  _max_mesh_elem_id = max_elem_id;
256  _max_mesh_node_id = max_node_id;
257 
258  _elem = &e;
259  _phi = &phi;
260 
261  switch (e.type()) {
262  case libMesh::QUAD4:
264  break;
265 
266  case libMesh::QUAD9: {
267 
268  std::unique_ptr<libMesh::Elem>
269  quad4(_first_order_elem(e));
270  _init_on_first_order_ref_elem(phi, *quad4, t);
271  }
272  break;
273 
274  default:
275  // currently only QUAD4/9 are handled.
276  libmesh_error();
277  }
278 }
279 
280 
281 
282 std::unique_ptr<libMesh::Elem>
284 
285  std::unique_ptr<libMesh::Elem>
286  first_order_elem;
287 
288  switch (e.type()) {
289  case libMesh::QUAD9: {
290 
291  first_order_elem.reset(libMesh::Elem::build(libMesh::QUAD4).release());
292 
293  for (unsigned int i=0; i<4; i++)
294  first_order_elem->set_node(i) = const_cast<libMesh::Node*>(e.node_ptr(i));
295 
296  first_order_elem->set_id() = e.id();
297  }
298  break;
299 
300  default:
301  // currently only QUAD4/9 are handled.
302  libmesh_error();
303  }
304 
305  return first_order_elem;
306 }
307 
308 
309 
310 
311 void
314  const libMesh::Elem& e,
315  const Real t) {
316 
317  // make sure that this has not been initialized already
318  libmesh_assert(!_initialized);
319  libmesh_assert_equal_to(e.dim(), 2); // this is only for 2D elements
320  libmesh_assert(_elem);
321 
322  // this assumes that the phi=0 interface is not fully contained inside
323  // the element. That is, the interface is assumed to intersect with the
324  // element edges.
325 
326  // check for each mode of intersection. It is assumed that an element only
327  // sees a single mode of intersection, which is hopefully true for a fine
328  // mesh. However, with larger high-order elements this assumption will
329  // not be valid and more accurate implementations will be needed.
330 
331  _node_phi_vals.clear();
332 
333  unsigned int
334  n_node_intersection = 0;
335 
336  Real
337  val = 0.,
338  min_val = 1.e12, // set arbitrary high values
339  max_val = -1.e12;
340 
341  bool
342  on_level_set = false;
343 
344  // the shape function values are checked on all nodes of the reference elem
345  for (unsigned int i=0; i<_elem->n_nodes(); i++) {
346 
347  const libMesh::Node& n = _elem->node_ref(i);
348  phi(n, t, val);
349  on_level_set = fabs(val) <= _tol;
350  _node_phi_vals[&n] = std::pair<Real, bool>(val, on_level_set);
351  // only check the corner nodes for intersection to maintain the
352  // assumption of simple intersections on quad4. This will be
353  // consistent for a quad4 reference elem, but for a quad8/9 elem
354  // this helps in simplifying the process.
355  if (i < e.n_nodes() && on_level_set) {
357  n_node_intersection++;
358  }
359  min_val = (min_val > val)? val : min_val;
360  max_val = (max_val < val)? val : max_val;
361  }
362 
363 
365  // Check to see if the whole element is on either side of the function
367 
368  // if the sign of function on all nodes is the same, then it is assumed
369  // that the element is not intersected
370  if (min_val > _tol) {
371  // element is completely on the positive side with no intersection
372 
374  _positive_phi_elems.push_back(_elem);
376  _if_elem_on_negative_phi = false;
377  _initialized = true;
378  return;
379  }
380  else if (max_val < -_tol) {
381 
382  // element is completely on the negative side, with no intersection
383 
385  _negative_phi_elems.push_back(_elem);
386  _if_elem_on_positive_phi = false;
388  _initialized = true;
389  return;
390  }
391  else {
392  // if it did not get caught in the previous two cases, then there is
393  // an intersection.
394 
395  _if_elem_on_positive_phi = false;
396  _if_elem_on_negative_phi = false;
397  }
398 
400  // Check to see if the function passes through one node
402  // if only one node intersection was found, then the level set
403  // passes through a node. Node intersection happens when the level set on
404  // a node is zero and the level set does not change sign on the element.
405  Real
406  val1 = std::fabs(min_val)<_tol ? 0.: min_val,
407  val2 = std::fabs(max_val)<_tol ? 0.: max_val;
408  bool
409  sign_change = (val1*val2 < 0.);
410  if (n_node_intersection == 1 &&
411  !sign_change) {
412 
414 
415  // this means that the whole element is going to be on either the
416  // positive or the negative side of the level-set function
417  // now figure out which side of the level-set function the
418  // element lies on
419 
420  if (min_val <= _tol &&
421  max_val <= _tol ) { // element is on the negative side
422  _negative_phi_elems.push_back(_elem);
423  }
424  else {
425 
426  libmesh_assert_greater(max_val, _tol);
427  _positive_phi_elems.push_back(_elem);
428  }
429 
430  std::vector<std::pair<libMesh::Point, libMesh::Point> >
431  side_nondim_points;
432  _add_node_local_coords(e, side_nondim_points, _node_local_coords);
433  _initialized = true;
434  return;
435  }
436 
437 
439  // Check to see if the function is conliniear with an edge
441  // Check to see if all nodes on an edge are on the level set
442  for (unsigned int i=0; i<e.n_sides(); i++) {
443 
444  std::unique_ptr<const libMesh::Elem>
445  s(e.side_ptr(i).release());
446 
447  unsigned int
448  n_nodes_on_level_set = 0;
449 
450  for (unsigned int j=0; j<s->n_nodes(); j++)
451  n_nodes_on_level_set += _node_phi_vals[s->node_ptr(j)].second;
452 
453  if (n_nodes_on_level_set == s->n_nodes()) {
454 
455  // side is on the level set
458 
459  // now that we know that one of the edges is on the function,
460  // we need to identify if the element is on the positive or
461  // the negative side of the function. It is on the positive side
462  // if the maximum value of the function on the nodes is positive.
463  // It is on the negative side if the minimum value of the function
464  // is negative.
465  if (max_val > _tol)
466  _positive_phi_elems.push_back(_elem);
467  else if (min_val < _tol)
468  _negative_phi_elems.push_back(_elem);
469 
470  std::vector<std::pair<libMesh::Point, libMesh::Point> >
471  side_nondim_points;
472  _add_node_local_coords(e, side_nondim_points, _node_local_coords);
474  _initialized = true;
475  return;
476  }
477  }
478 
480  // Identify the edges and locations where intersection occurs
482  // now that we are here, this means that none of the previous modes
483  // identify the intersection of the level set with this element.
484  // Therefore, we will now identify the sides of the element where
485  // intersectio occurs and create the sub-elements for integration.
486  // This is done separately for each element type.
487  switch (e.type()) {
488 
489  case libMesh::QUAD4:
491  break;
492 
493  default:
494  libmesh_error_msg("level-set intersections for elem type not handled.");
495  break;
496  }
497 
498  // set the IDs of the new elements. Both the subdomain ID and
499  // element IDs are set here.
500  //
501  // Note that the subdomain IDs are needed
502  // for querying the BCs, Loads and properties from the physics. We
503  // use the same subdomain ID as its parent ID
504  //
505  // The element id is a little more tricky. Strictly speaking, this
506  // should not be necessary, since we are only creating temporary
507  // elements to compute element quantities. However, the current stress
508  // storage implementation is storing this for element ids. So, we provide
509  // a unique element ID for each new element.
510  //
511 
512  libmesh_assert_less_equal(_new_elems.size(), _max_elem_divs);
513 
514  for (unsigned int i=0; i<_new_elems.size(); i++) {
515 
516  _new_elems[i]->subdomain_id() = e.subdomain_id();
517  _new_elems[i]->set_id(_max_mesh_elem_id+e.id()*_max_elem_divs+i);
518  }
519 
520  _initialized = true;
521 }
522 
523 
524 
525 
526 const std::vector<const libMesh::Elem*>&
528 
529  return _positive_phi_elems;
530 }
531 
532 
533 
534 const std::vector<const libMesh::Elem*>&
536 
537  return _negative_phi_elems;
538 }
539 
540 
541 void
543 get_corner_nodes_on_negative_phi(std::set<const libMesh::Node*>& nodes) const {
544 
545  nodes.clear();
546 
547  // iterate over the corner nodes and return the ones that had a negative phi
548  // value
549  unsigned int
550  n_corner_nodes = 0;
551 
552  switch (_elem->type()) {
553  case libMesh::QUAD4:
554  case libMesh::QUAD8:
555  case libMesh::QUAD9:
556  n_corner_nodes = 4;
557  break;
558 
559  default:
560  libmesh_error(); // other cases not yet handled
561  }
562 
563  std::map<const libMesh::Node*, std::pair<Real, bool>>::const_iterator
564  it,
565  end = _node_phi_vals.end();
566 
567  for (unsigned int i=0; i<n_corner_nodes; i++) {
568 
569  const libMesh::Node* nd = _elem->node_ptr(i);
570 
571  it = _node_phi_vals.find(nd);
572 
573  libmesh_assert(it != end);
574 
575  if (it->second.first < 0.)
576  nodes.insert(it->first);
577  }
578 }
579 
580 
581 
582 bool
584 has_side_on_interface(const libMesh::Elem& e) const {
585 
586  std::map<const libMesh::Elem*, int>::const_iterator it;
587  it = _elem_sides_on_interface.find(&e);
588 
589  libmesh_assert(it != _elem_sides_on_interface.end());
590 
591  return it->second >= 0;
592 }
593 
594 
595 
596 unsigned int
598 get_side_on_interface(const libMesh::Elem& e) const {
599 
600  std::map<const libMesh::Elem*, int>::const_iterator it;
601  it = _elem_sides_on_interface.find(&e);
602 
603  libmesh_assert(it != _elem_sides_on_interface.end());
604  libmesh_assert(it->second >= 0);
605 
606  return it->second;
607 }
608 
609 
610 void
612 ( const libMesh::Elem& e,
613  std::vector<std::pair<libMesh::Point, libMesh::Point> >& side_nondim_points,
614  std::map<const libMesh::Node*, libMesh::Point>& node_coord_map) {
615 
616  switch (e.type()) {
617 
618  case libMesh::QUAD4: {
619 
620  // We create a vector storing the nondimensional locations of nodes
621  // on each side of the element. This is used later to identify the
622  // nondimensional location of all new nodes
623  side_nondim_points.resize(4);
624  // side 0: eta =-1, xi=[-1, 1]
625  side_nondim_points[0] =
626  std::pair<libMesh::Point, libMesh::Point>
627  (libMesh::Point(-1, -1, 0.),
628  libMesh::Point(+1, -1, 0.));
629  side_nondim_points[1] =
630  std::pair<libMesh::Point, libMesh::Point>
631  (libMesh::Point(+1, -1, 0.),
632  libMesh::Point(+1, +1, 0.));
633  side_nondim_points[2] =
634  std::pair<libMesh::Point, libMesh::Point>
635  (libMesh::Point(+1, +1, 0.),
636  libMesh::Point(-1, +1, 0.));
637  side_nondim_points[3] =
638  std::pair<libMesh::Point, libMesh::Point>
639  (libMesh::Point(-1, +1, 0.),
640  libMesh::Point(-1, -1, 0.));
641 
642  // populate the node to nondimensional coordinate map with the
643  // original nodes of the element
644  for (unsigned int i=0; i<4; i++)
645  node_coord_map[e.node_ptr(i)] = side_nondim_points[i].first;
646  }
647  break;
648 
649  default:
650  libmesh_assert(false); // not handled.
651  }
652 
653 }
654 
655 
656 void
659  const libMesh::Elem& e,
660  const Real t,
661  const std::map<const libMesh::Node*, std::pair<Real, bool> >&
662  node_phi_vals) {
663 
664  libmesh_assert_equal_to(e.type(), libMesh::QUAD4);
665  libmesh_assert(!_initialized);
666 
667  std::vector<std::pair<libMesh::Point, libMesh::Point> >
668  side_nondim_points;
669  _add_node_local_coords(e, side_nondim_points, _node_local_coords);
670 
671  const unsigned int
672  n_nodes = 4; // this is equal to n_sides for QUAD4
673 
674  std::vector<std::pair<bool, Real> >
675  side_intersection(n_nodes, std::pair<bool, Real>(false, 0.));
676  std::vector<bool>
677  node_intersection(n_nodes, false);
678 
679  std::map<const libMesh::Node*, std::pair<Real, bool> >::const_iterator
680  it0,
681  it1,
682  it_end = node_phi_vals.end();
683 
684  Real
685  v0 = 0.,
686  v1 = 0.;
687 
688 
689  for (unsigned int i=0; i<n_nodes; i++) {
690 
691  it0 = node_phi_vals.find(e.node_ptr(i));
692  libmesh_assert(it0 != it_end);
693 
694  // identify intersection with nodes
695  if (std::fabs(it0->second.first) < _tol)
696  node_intersection[i] = true;
697  else { // if there is node intersection, then there won't be edge interseciton
698 
699  // identify intersection with edges
700  it1 = node_phi_vals.find(e.node_ptr((i+1)%n_nodes));
701  libmesh_assert(it1 != it_end);
702  v0 = it0->second.first;
703  v1 = it1->second.first;
704  if (std::fabs(v0) > _tol &&
705  std::fabs(v1) > _tol &&
706  v0 * v1 < 0.) {
707 
708  // intersection lies somewhere in the interior of the edge
709  side_intersection[i] = std::pair<bool, Real>
710  (true,
712  *it1->first,
713  phi,
714  t));
715  }
716  }
717  }
718 
719 
720 
721  // currently, we only handle two intersections
722  unsigned int
723  n_intersection = 0;
724 
725  // add up intersections on edges
726  {
727  std::vector<std::pair<bool, Real> >::const_iterator
728  it = side_intersection.begin(),
729  end = side_intersection.end();
730 
731  for (; it != end; it++)
732  n_intersection += it->first;
733  }
734 
735  // add up intersections on nodes
736  {
737  std::vector<bool>::const_iterator
738  it = node_intersection.begin(),
739  end = node_intersection.end();
740 
741  for (; it != end; it++)
742  n_intersection += *it;
743  }
744 
745  //
746  // identify the sets of intersections, which will be processed
747  // individually for creation of the subelements.
748  // The current implementation of linear quad4s only permits one intersection
749  // per edge. This will have to be chanegd for high-order elements.
750  //
751  // As a result of this assumption, the following should hold:
752  // -- OPPOSITE_EDGES: only two intersections can exist
753  // -- OPPOSITE_NODES: only two nodes can intersect that ways
754  // -- NODE_AND_EDGE: two such intersections can occur
755  // -- ADJACENT_EDGES: two such intersections can occur
756  //
757 
758  // now identify which approch to use for creation of the sub-elements.
759  // one case is where two adjacent sides of the element have a
760  unsigned int
761  ref_side = 0,
762  ref_node = 0;
763 
764  switch (n_intersection) {
765  case 2: {
766  for (unsigned int i=0; i<n_nodes; i++) {
767 
768  if (side_intersection[i].first &&
769  side_intersection[(i+1)%n_nodes].first) {
770  // found adjacent mode with ref edge i
772  ref_side = i;
773  break;
774  }
775  else if (side_intersection[i].first &&
776  side_intersection[(i+2)%n_nodes].first) {
777  // found adjacent mode with ref edge i
779  ref_side = i;
780  break;
781  }
782  else if (node_intersection[i] &&
783  node_intersection[(i+2)%n_nodes]) {
784  // found intersection on two opposite nodes
786  ref_node = i;
787  break;
788  }
789  else if (node_intersection[i] &&
790  side_intersection[(i+1)%n_nodes].first) {
791  // found intersection on node and edge
793  ref_node = i;
794  ref_side = i+1;
795  }
796  else if (node_intersection[i] &&
797  side_intersection[(i+2)%n_nodes].first) {
798  // found intersection on node and edge
800  ref_node = i;
801  ref_side = i+2;
802  }
803  }
804  }
805  break;
806 
807  case 3: {
808  // this must be a node with two edge intersection
809  // we first identify the node
810  for (unsigned int i=0; i<n_nodes; i++) {
811  if (node_intersection[i] &&
812  side_intersection[(i+1)%n_nodes].first &&
813  side_intersection[(i+2)%n_nodes].first) {
814 
816  ref_node = i;
817  ref_side = i+1;
818  break;
819  }
820  }
821  }
822  break;
823 
824  case 4: {
825  // this must be two adjacent edge intersections, where all four edges
826  // should have intersections
827  if (side_intersection[0].first &&
828  side_intersection[1].first &&
829  side_intersection[2].first &&
830  side_intersection[3].first) {
831 
833 
834  //
835  // this means that all four edges have an intersection, one
836  // per edge. So, now we identify how to connect the edges.
837  // The choices are: intersection in sides 0-1 and 2-3, or
838  // intersction in side 0-3 and 1-2.
839  //
840  // For this, we need to figure out the gradient. We will
841  // do so by checking the value of the level-set function
842  // in the middle of the element. If it is the same sign as
843  // the zero-th node, then the intersection will be 0-1 and 2-3,
844  // if it is the opposite sign, then the intersection is
845  // 0-3 and 1-2.
846  //
847 
848  Real mid_phi = 0.;
849  phi(e.centroid(), t, mid_phi);
850  it0 = node_phi_vals.find(e.node_ptr(0));
851  v0 = it0->second.first;
852 
853  if (mid_phi * v0 > 0.) // same sign
854  ref_side = 0;
855  else
856  ref_side = 1;
857  }
858  }
859  break;
860 
861  default:
862  // if it gets here, then we have a mode that we have not
863  // created
864  assert(false);
865  }
866 
867  // make sure that atleast one of the modes was found
868  libmesh_assert(_mode == MAST::ADJACENT_EDGES ||
874 
876  // now create the sub elements
878  // For intersection with opposite edges, both sides of the element
879  // will be represented by QUAD4 elements.
880  //
881  // When the intersection on adjacent elements happens, then
882  // one side will be represented with triangles and the other
883  // side is represented with a combination of two quads.
884  //
885  // for opposite nodes, one triangle is used for each side of the
886  // intersection
887  //
888  // for node and edge, one side is represented with a triangle,
889  // and the other with a quad.
890  //
891 
892  Real
893  xi_ref = 0.,
894  xi_other = 0.;
895  libMesh::Point
896  p,
897  side_p0,
898  side_p1;
899 
900  // this should pass, but just to be sure, we will check.
901  libmesh_assert_equal_to(_new_nodes.size(), 0);
902 
903  libMesh::Node
904  *nd = nullptr;
905 
906  std::unique_ptr<const libMesh::Elem>
907  s;
908 
909  bool
910  ref_side_node0_positive = false;
911 
912  // used this to set unique node ids, since some of libMesh's operations
913  // depend on valid and unique ids for nodes.
914  unsigned int
915  node_id_incr = 1;
916 
918  // Create a new node at each intersection: first, on ref_side
919  //
920  // New nodes are needed only when intersection does not exist
921  // opposite nodes. In this case, the current nodes are used for
922  // the new elements.
924  if (_mode == MAST::ADJACENT_EDGES ||
930 
931  xi_ref = side_intersection[ref_side%n_nodes].second;
932  s.reset(e.side_ptr(ref_side%n_nodes).release());
933  p = s->point(0) + xi_ref * (s->point(1) - s->point(0));
934  nd = new libMesh::Node(p);
935  nd->set_id(_max_mesh_node_id + node_id_incr);
936  node_id_incr++;
937  _new_nodes.push_back(nd);
938  _bounding_nodes[nd] = std::make_pair(_elem->node_ptr(ref_side%n_nodes),
939  _elem->node_ptr((ref_side+1)%n_nodes));
940  side_p0 = side_nondim_points[ref_side%n_nodes].first;
941  side_p1 = side_nondim_points[ref_side%n_nodes].second;
942  _node_local_coords[nd] = side_p0 + xi_ref * (side_p1 - side_p0);
943 
944 
946  // before we create the new elements, we will figure out
947  // which portion of the side is on the positive or negative side of
948  // the level set function
950 
951  it0 = node_phi_vals.find(s->node_ptr(0));
952  v0 = it0->second.first;
953 
954  if (v0 > 0.)
955  ref_side_node0_positive = true;
956  }
957  else {
958 
959  it0 = node_phi_vals.find(e.node_ptr(ref_node+1));
960  v0 = it0->second.first;
961 
962  if (v0 > 0.)
963  ref_side_node0_positive = true;
964 
965  }
966 
967 
968 
969  switch (_mode) {
970 
971  case MAST::OPPOSITE_EDGES: {
972 
973  // create a new node at each intersection: next, on ref_side + 2
974  xi_other = side_intersection[(ref_side+2)%n_nodes].second;
975  s.reset(e.side_ptr((ref_side+2)%n_nodes).release());
976  p = s->point(0) + xi_other * (s->point(1) - s->point(0));
977  nd = new libMesh::Node(p);
978  nd->set_id(_max_mesh_node_id + node_id_incr);
979  node_id_incr++;
980  _new_nodes.push_back(nd);
981  _bounding_nodes[nd] = std::make_pair(_elem->node_ptr((ref_side+2)%n_nodes),
982  _elem->node_ptr((ref_side+3)%n_nodes));
983  side_p0 = side_nondim_points[(ref_side+2)%n_nodes].first;
984  side_p1 = side_nondim_points[(ref_side+2)%n_nodes].second;
985  _node_local_coords[nd] = side_p0 + xi_other * (side_p1 - side_p0);
986 
987  libMesh::Elem
988  *e_p = const_cast<libMesh::Elem*>(&e),
989  *e1 = libMesh::Elem::build(libMesh::QUAD4, e_p).release(),
990  *e2 = libMesh::Elem::build(libMesh::QUAD4, e_p).release();
991 
992  _new_elems.push_back(e1);
993  _new_elems.push_back(e2);
994 
995  // create the elements. We create the elements such that:
996  // node(0) = new node on ref_side
997  // node(1) = elem node (ref_side+1)%n_sides
998  // node(2) = elem node (ref_side+2)%n_sides
999  // node(3) = new node on ref_side+2
1000  // this way, the element has the same orientation as the original
1001  // element and the sides have the following association:
1002  // side(0) = coincident with ref_side of original element
1003  // side(1) = coincident side with ref_side+1 of original element
1004  // side(2) = coincident side with ref_side+2 of original element
1005  // side(3) = coincident side with level-set interface
1006 
1007  e1->set_node(0) = _new_nodes[0];
1008  e1->set_node(1) = const_cast<libMesh::Node*>(e.node_ptr((ref_side+1)%n_nodes));
1009  e1->set_node(2) = const_cast<libMesh::Node*>(e.node_ptr((ref_side+2)%n_nodes));
1010  e1->set_node(3) = _new_nodes[1];
1011 
1012  _elem_sides_on_interface[e1] = 3;
1013 
1014  // create a second element and set its nodes so that:
1015  // node(0) = new node on ref_side+2
1016  // node(1) = elem node (ref_side+3)%n_nodes
1017  // node(2) = elem node (ref_side+4)%n_nodes
1018  // node(3) = new node on ref_side
1019  // this way, the element has the same orientation as the original
1020  // element and the sides have the following association:
1021  // side(0) = coincident with ref_side+2 of original element
1022  // side(1) = coincident side with (ref_side+3)%n_sides of original element
1023  // side(2) = coincident side with ref_side of original element
1024  // side(3) = coincident side with level-set interface
1025 
1026  e2->set_node(0) = _new_nodes[1];
1027  e2->set_node(1) = const_cast<libMesh::Node*>(e.node_ptr((ref_side+3)%n_nodes));
1028  e2->set_node(2) = const_cast<libMesh::Node*>(e.node_ptr(ref_side));
1029  e2->set_node(3) = _new_nodes[0];
1030 
1031  _elem_sides_on_interface[e2] = 3;
1032 
1033  if (ref_side_node0_positive) {
1034 
1035  _positive_phi_elems.push_back(e2);
1036  _negative_phi_elems.push_back(e1);
1037  }
1038  else {
1039 
1040  _positive_phi_elems.push_back(e1);
1041  _negative_phi_elems.push_back(e2);
1042  }
1043  }
1044  break;
1045 
1046  case MAST::ADJACENT_EDGES: {
1047 
1048 
1049  // create a new node at each intersection: second on ref_side+1
1050  xi_other = side_intersection[(ref_side+1)%n_nodes].second;
1051  s.reset(e.side_ptr((ref_side+1)%n_nodes).release());
1052  p = s->point(0) + xi_other * (s->point(1) - s->point(0));
1053  nd = new libMesh::Node(p);
1054  nd->set_id(_max_mesh_node_id + node_id_incr);
1055  node_id_incr++;
1056  _new_nodes.push_back(nd);
1057  _bounding_nodes[nd] = std::make_pair(_elem->node_ptr((ref_side+1)%n_nodes),
1058  _elem->node_ptr((ref_side+2)%n_nodes));
1059  side_p0 = side_nondim_points[(ref_side+1)%n_nodes].first;
1060  side_p1 = side_nondim_points[(ref_side+1)%n_nodes].second;
1061  _node_local_coords[nd] = side_p0 + xi_other * (side_p1 - side_p0);
1062 
1063 
1064  // create a new node on the opposite side, which will be used to create
1065  // the QUAD4 elements: second on ref_side+2
1066  s.reset(e.side_ptr((ref_side+2)%n_nodes).release());
1067  p = s->point(1) + xi_ref * (s->point(0) - s->point(1));
1068  nd = new libMesh::Node(p);
1069  nd->set_id(_max_mesh_node_id + node_id_incr);
1070  node_id_incr++;
1071  _new_nodes.push_back(nd);
1072  _bounding_nodes[nd] = std::make_pair(_elem->node_ptr((ref_side+2)%n_nodes),
1073  _elem->node_ptr((ref_side+3)%n_nodes));
1074  side_p0 = side_nondim_points[(ref_side+2)%n_nodes].first;
1075  side_p1 = side_nondim_points[(ref_side+2)%n_nodes].second;
1076  _node_local_coords[nd] = side_p1 + xi_ref * (side_p0 - side_p1);
1077 
1078 
1079  // the algorithm of identifying the ref_side always ensures that the
1080  // region with xi > x0 on ref_side will be triangular. The other side
1081  // will be modeled with two QUAD4 elements
1082 
1083  libMesh::Elem
1084  *e_p = const_cast<libMesh::Elem*>(&e),
1085  *e1 = libMesh::Elem::build( libMesh::TRI3, e_p).release(),
1086  *e2 = libMesh::Elem::build(libMesh::QUAD4, e_p).release(),
1087  *e3 = libMesh::Elem::build(libMesh::QUAD4, e_p).release();
1088 
1089  _new_elems.push_back(e1);
1090  _new_elems.push_back(e2);
1091  _new_elems.push_back(e3);
1092 
1093  // create the elements. We create the elements such that:
1094  // node(0) = new node on ref_side
1095  // node(1) = elem node (ref_side+1)%n_sides
1096  // node(2) = new node on (ref_side+1)%n_sides
1097  // this way, the element has the same orientation as the original
1098  // element and the sides have the following association:
1099  // side(0) = coincident with ref_side of original element
1100  // side(1) = coincident side with ref_side+1 of original element
1101  // side(2) = coincident side with level-set interface
1102 
1103  e1->set_node(0) = _new_nodes[0];
1104  e1->set_node(1) = const_cast<libMesh::Node*>(e.node_ptr((ref_side+1)%n_nodes));
1105  e1->set_node(2) = _new_nodes[1];
1106 
1107  _elem_sides_on_interface[e1] = 2;
1108 
1109  // create a second element and set its nodes so that:
1110  // node(0) = new node on ref_side
1111  // node(1) = new node on (ref_side+2)%n_nodes
1112  // node(2) = elem node (ref_side+3)%n_nodes
1113  // node(3) = elem node ref_side
1114  // this way, the element has the same orientation as the original
1115  // element and the sides have the following association:
1116  // side(0) = new edge connecting ref_side to new node on ref_side+2
1117  // side(1) = coincident side with (ref_side+2)%n_sides of original element
1118  // side(2) = coincident side with (ref_side+3)%n_sides of original element
1119  // side(3) = coincident with ref_side of original element
1120 
1121  e2->set_node(0) = _new_nodes[0];
1122  e2->set_node(1) = _new_nodes[2];
1123  _hanging_node.insert(_new_nodes[2]);
1124  e2->set_node(2) = const_cast<libMesh::Node*>(e.node_ptr((ref_side+3)%n_nodes));
1125  e2->set_node(3) = const_cast<libMesh::Node*>(e.node_ptr(ref_side));
1126 
1127  // this element does not have a side on the interface
1128  _elem_sides_on_interface[e2] = -1;
1129 
1130  // create a third element and set its nodes so that:
1131  // node(0) = new node on ref_side+1
1132  // node(1) = elem node (ref_side+2)%n_nodes
1133  // node(2) = new node on (ref_side+2)%n_nodes
1134  // node(3) = new node on ref_side
1135  // this way, the element has the same orientation as the original
1136  // element and the sides have the following association:
1137  // side(0) = coincident side with (ref_side+1)%n_sides of original element
1138  // side(1) = coincident side with (ref_side+2)%n_sides of original element
1139  // side(2) = new edge connecting ref_side to new node on (ref_side+2)%n_sides
1140  // side(3) = coincident side with level-set interface
1141 
1142  e3->set_node(0) = _new_nodes[1];
1143  e3->set_node(1) = const_cast<libMesh::Node*>(e.node_ptr((ref_side+2)%n_nodes));
1144  e3->set_node(2) = _new_nodes[2];
1145  e3->set_node(3) = _new_nodes[0];
1146 
1147  _elem_sides_on_interface[e3] = 3;
1148 
1149  if (ref_side_node0_positive) {
1150 
1151  _positive_phi_elems.push_back(e2);
1152  _positive_phi_elems.push_back(e3);
1153  _negative_phi_elems.push_back(e1);
1154  }
1155  else {
1156 
1157  _positive_phi_elems.push_back(e1);
1158  _negative_phi_elems.push_back(e2);
1159  _negative_phi_elems.push_back(e3);
1160  }
1161  }
1162  break;
1163 
1164 
1165 
1166  case MAST::NODE_AND_EDGE: {
1167 
1168  libMesh::Elem
1169  *e_p = const_cast<libMesh::Elem*>(&e),
1170  *e1 = libMesh::Elem::build(libMesh::TRI3, e_p).release(),
1171  *e2 = libMesh::Elem::build(libMesh::QUAD4, e_p).release();
1172 
1173  _new_elems.push_back(e1);
1174  _new_elems.push_back(e2);
1175 
1176  if (ref_side-ref_node == 1) {
1177 
1178  // create the elements. We create the elements such that:
1179  // node(0) = elem node (ref_node)
1180  // node(1) = elem node (ref_node+1)%n_nodes
1181  // node(2) = new node on ref_side%n_nodes
1182  // this way, the element has the same orientation as the original
1183  // element and the sides have the following association:
1184  // side(0) = coincident with ref_node of original element
1185  // side(1) = coincident side with ref_node+1 of original element
1186  // side(2) = coincident side with level-set interface
1187 
1188  e1->set_node(0) = const_cast<libMesh::Node*>(e.node_ptr((ref_node)%n_nodes));
1189  e1->set_node(1) = const_cast<libMesh::Node*>(e.node_ptr((ref_node+1)%n_nodes));
1190  e1->set_node(2) = _new_nodes[0];
1191 
1192  _elem_sides_on_interface[e1] = 2;
1193 
1194  // create a second element and set its nodes so that:
1195  // node(0) = new node on ref_side%n_nodes
1196  // node(1) = elem node (ref_side+1)%n_nodes
1197  // node(2) = elem node (ref_side+2)%n_nodes
1198  // node(3) = elem node (ref_node)
1199  // this way, the element has the same orientation as the original
1200  // element and the sides have the following association:
1201  // side(0) = coincident with ref_side of original element
1202  // side(1) = coincident side with (ref_side+1)%n_sides of original element
1203  // side(2) = coincident side with (ref_side+2)%n_sides of original element
1204  // side(3) = coincident side with level-set interface
1205 
1206  e2->set_node(0) = _new_nodes[0];
1207  e2->set_node(1) = const_cast<libMesh::Node*>(e.node_ptr((ref_side+1)%n_nodes));
1208  e2->set_node(2) = const_cast<libMesh::Node*>(e.node_ptr((ref_side+2)%n_nodes));
1209  e2->set_node(3) = const_cast<libMesh::Node*>(e.node_ptr((ref_node)%n_nodes));
1210 
1211  _elem_sides_on_interface[e2] = 3;
1212 
1213 
1214  // node0 is the first node of ref_side
1215  if (ref_side_node0_positive) {
1216 
1217  _positive_phi_elems.push_back(e1);
1218  _negative_phi_elems.push_back(e2);
1219  }
1220  else {
1221 
1222  _positive_phi_elems.push_back(e2);
1223  _negative_phi_elems.push_back(e1);
1224  }
1225  }
1226  else if (ref_side-ref_node==2) {
1227 
1228  // create the elements. We create the elements such that:
1229  // node(0) = elem node (ref_node)
1230  // node(1) = new node on ref_side%n_nodes
1231  // node(2) = elem node (ref_node+3)%n_nodes
1232  // this way, the element has the same orientation as the original
1233  // element and the sides have the following association:
1234  // side(0) = coincident side with level-set interface
1235  // side(1) = coincident side with ref_side%n_nodes of original element
1236  // side(2) = coincident side with (ref_side+1)%n_nodes of original element
1237 
1238  e1->set_node(0) = const_cast<libMesh::Node*>(e.node_ptr((ref_node)%n_nodes));
1239  e1->set_node(1) = _new_nodes[0];
1240  e1->set_node(2) = const_cast<libMesh::Node*>(e.node_ptr((ref_node+3)%n_nodes));
1241 
1242  _elem_sides_on_interface[e1] = 0;
1243 
1244  // create a second element and set its nodes so that:
1245  // node(0) = elem node (ref_node)
1246  // node(1) = elem node (ref_node+1)%n_nodes
1247  // node(2) = elem node (ref_node+2)%n_nodes
1248  // node(3) = new node on ref_side%n_nodes
1249  // this way, the element has the same orientation as the original
1250  // element and the sides have the following association:
1251  // side(0) = coincident with ref_node of original element
1252  // side(1) = coincident side with (ref_node+1)%n_sides of original element
1253  // side(2) = coincident side with (ref_node+2)%n_sides of original element
1254  // side(3) = coincident side with level-set interface
1255 
1256  e2->set_node(0) = const_cast<libMesh::Node*>(e.node_ptr((ref_node)%n_nodes));
1257  e2->set_node(1) = const_cast<libMesh::Node*>(e.node_ptr((ref_node+1)%n_nodes));
1258  e2->set_node(2) = const_cast<libMesh::Node*>(e.node_ptr((ref_node+2)%n_nodes));
1259  e2->set_node(3) = _new_nodes[0];
1260 
1261  _elem_sides_on_interface[e2] = 3;
1262 
1263 
1264  if (ref_side_node0_positive) {
1265 
1266  _positive_phi_elems.push_back(e2);
1267  _negative_phi_elems.push_back(e1);
1268  }
1269  else {
1270 
1271  _positive_phi_elems.push_back(e1);
1272  _negative_phi_elems.push_back(e2);
1273  }
1274  }
1275  else
1276  libmesh_error(); // shoudl not get here
1277 
1278  }
1279  break;
1280 
1281 
1282  case MAST::OPPOSITE_NODES: {
1283 
1284  libMesh::Elem
1285  *e_p = const_cast<libMesh::Elem*>(&e),
1286  *e1 = libMesh::Elem::build(libMesh::TRI3, e_p).release(),
1287  *e2 = libMesh::Elem::build(libMesh::TRI3, e_p).release();
1288 
1289  _new_elems.push_back(e1);
1290  _new_elems.push_back(e2);
1291 
1292  // create the elements. We create the elements such that:
1293  // node(0) = elem node (ref_node)
1294  // node(1) = elem node (ref_node+1)%n_nodes
1295  // node(2) = elem node (ref_node+2)%n_nodes
1296  // this way, the element has the same orientation as the original
1297  // element and the sides have the following association:
1298  // side(0) = coincident with ref_node of original element
1299  // side(1) = coincident side with ref_node+1 of original element
1300  // side(2) = coincident side with level-set interface
1301 
1302  e1->set_node(0) = const_cast<libMesh::Node*>(e.node_ptr((ref_node)%n_nodes));
1303  e1->set_node(1) = const_cast<libMesh::Node*>(e.node_ptr((ref_node+1)%n_nodes));
1304  e1->set_node(2) = const_cast<libMesh::Node*>(e.node_ptr((ref_node+2)%n_nodes));
1305 
1306  _elem_sides_on_interface[e1] = 2;
1307 
1308  // create a second element and set its nodes so that:
1309  // node(0) = elem node on (ref_node+2)%n_nodes
1310  // node(1) = elem node on (ref_node+3)%n_nodes
1311  // node(2) = elem node on (ref_node+4)%n_nodes
1312  // this way, the element has the same orientation as the original
1313  // element and the sides have the following association:
1314  // side(0) = coincident with (ref_node+2)%n_nodes of original element
1315  // side(1) = coincident with (ref_node+3)%n_nodes of original element
1316  // side(2) = coincident side with level-set interface
1317 
1318  e2->set_node(0) = const_cast<libMesh::Node*>(e.node_ptr((ref_node+2)%n_nodes));
1319  e2->set_node(1) = const_cast<libMesh::Node*>(e.node_ptr((ref_node+3)%n_nodes));
1320  e2->set_node(2) = const_cast<libMesh::Node*>(e.node_ptr(ref_node));
1321 
1322  _elem_sides_on_interface[e2] = 2;
1323 
1324  if (ref_side_node0_positive) {
1325 
1326  _positive_phi_elems.push_back(e1);
1327  _negative_phi_elems.push_back(e2);
1328  }
1329  else {
1330 
1331  _positive_phi_elems.push_back(e2);
1332  _negative_phi_elems.push_back(e1);
1333  }
1334  }
1335  break;
1336 
1337 
1338  case MAST::NODE_AND_TWO_EDGES: {
1339 
1340  // create a new node at each intersection: second on ref_side+1
1341  xi_other = side_intersection[(ref_side+1)%n_nodes].second;
1342  s.reset(e.side_ptr((ref_side+1)%n_nodes).release());
1343  p = s->point(0) + xi_other * (s->point(1) - s->point(0));
1344  nd = new libMesh::Node(p);
1345  nd->set_id(_max_mesh_node_id + node_id_incr);
1346  node_id_incr++;
1347  _new_nodes.push_back(nd);
1348  _bounding_nodes[nd] = std::make_pair(_elem->node_ptr((ref_side+1)%n_nodes),
1349  _elem->node_ptr((ref_side+2)%n_nodes));
1350  side_p0 = side_nondim_points[(ref_side+1)%n_nodes].first;
1351  side_p1 = side_nondim_points[(ref_side+1)%n_nodes].second;
1352  _node_local_coords[nd] = side_p0 + xi_other * (side_p1 - side_p0);
1353 
1354 
1355  libMesh::Elem
1356  *e_p = const_cast<libMesh::Elem*>(&e),
1357  *e1 = libMesh::Elem::build(libMesh::TRI3, e_p).release(),
1358  *e2 = libMesh::Elem::build(libMesh::TRI3, e_p).release(),
1359  *e3 = libMesh::Elem::build(libMesh::TRI3, e_p).release(),
1360  *e4 = libMesh::Elem::build(libMesh::TRI3, e_p).release();
1361 
1362  _new_elems.push_back(e1);
1363  _new_elems.push_back(e2);
1364  _new_elems.push_back(e3);
1365  _new_elems.push_back(e4);
1366 
1367  // create the elements. We create the elements such that:
1368  // node(0) = elem node (ref_node)
1369  // node(1) = elem node (ref_node+1)%n_nodes
1370  // node(2) = new node on ref_side%n_nodes
1371  // this way, the element has the same orientation as the original
1372  // element and the sides have the following association:
1373  // side(0) = coincident with ref_node of original element
1374  // side(1) = coincident side with ref_node+1 of original element
1375  // side(2) = coincident side with level-set interface
1376 
1377  e1->set_node(0) = const_cast<libMesh::Node*>(e.node_ptr((ref_node)%n_nodes));
1378  e1->set_node(1) = const_cast<libMesh::Node*>(e.node_ptr((ref_node+1)%n_nodes));
1379  e1->set_node(2) = _new_nodes[0];
1380 
1381  _elem_sides_on_interface[e1] = 2;
1382 
1383 
1384  // create the elements. We create the elements such that:
1385  // node(0) = elem node (ref_node)
1386  // node(1) = new node on (ref_side+1)%n_nodes
1387  // node(2) = elem node (ref_node+3)%n_nodes
1388  // this way, the element has the same orientation as the original
1389  // element and the sides have the following association:
1390  // side(0) = coincident side with level-set interface
1391  // side(1) = coincident side with (ref_side+1)%n_nodes of original element
1392  // side(2) = coincident side with (ref_side+2)%n_nodes of original element
1393 
1394  e2->set_node(0) = const_cast<libMesh::Node*>(e.node_ptr((ref_node)%n_nodes));
1395  e2->set_node(1) = _new_nodes[1];
1396  e2->set_node(2) = const_cast<libMesh::Node*>(e.node_ptr((ref_node+3)%n_nodes));
1397 
1398  _elem_sides_on_interface[e2] = 0;
1399 
1400 
1401  // create a third element and set its nodes so that:
1402  // node(0) = elem node ref_node
1403  // node(1) = new node on (ref_side)%n_nodes
1404  // node(2) = elem node (ref_node+2)%n_nodes
1405  // this way, the element has the same orientation as the original
1406  // element and the sides have the following association:
1407  // side(0) = coincident side with level-set interface
1408  // side(1) = coincident side with (ref_side)%n_sides of original element
1409  // side(2) = diagonal between elem nodes (ref_node) and (ref_node+2)%n_sides
1410 
1411  e3->set_node(0) = const_cast<libMesh::Node*>(e.node_ptr((ref_node)%n_nodes));
1412  e3->set_node(1) = _new_nodes[0];
1413  e3->set_node(2) = const_cast<libMesh::Node*>(e.node_ptr((ref_node+2)%n_nodes));
1414 
1415  _elem_sides_on_interface[e3] = 0;
1416 
1417 
1418  // create a fourth element and set its nodes so that:
1419  // node(0) = elem node ref_node
1420  // node(1) = new node on (ref_side)%n_nodes
1421  // node(2) = elem node (ref_node+2)%n_nodes
1422  // this way, the element has the same orientation as the original
1423  // element and the sides have the following association:
1424  // side(0) = coincident side with level-set interface
1425  // side(1) = coincident side with (ref_side)%n_sides of original element
1426  // side(2) = diagonal between elem nodes (ref_node) and (ref_node+2)%n_sides
1427 
1428  e4->set_node(0) = const_cast<libMesh::Node*>(e.node_ptr((ref_node)%n_nodes));
1429  e4->set_node(1) = const_cast<libMesh::Node*>(e.node_ptr((ref_node+2)%n_nodes));
1430  e4->set_node(2) = _new_nodes[1];
1431 
1432  _elem_sides_on_interface[e4] = 2;
1433 
1434 
1435  // node0 is the first node of ref_side
1436  if (ref_side_node0_positive) {
1437 
1438  _positive_phi_elems.push_back(e1);
1439  _positive_phi_elems.push_back(e2);
1440  _negative_phi_elems.push_back(e3);
1441  _negative_phi_elems.push_back(e4);
1442  }
1443  else {
1444 
1445  _positive_phi_elems.push_back(e3);
1446  _positive_phi_elems.push_back(e4);
1447  _negative_phi_elems.push_back(e1);
1448  _negative_phi_elems.push_back(e2);
1449  }
1450  }
1451  break;
1452 
1453  case MAST::TWO_ADJACENT_EDGES: {
1454 
1455  // create a new node at each intersection: second on ref_side+1
1456  for (unsigned int i=1; i<4; i++) {
1457 
1458  xi_other = side_intersection[(ref_side+i)%n_nodes].second;
1459  s.reset(e.side_ptr((ref_side+i)%n_nodes).release());
1460  p = s->point(0) + xi_other * (s->point(1) - s->point(0));
1461  nd = new libMesh::Node(p);
1462  nd->set_id(_max_mesh_node_id + node_id_incr);
1463  node_id_incr++;
1464  _new_nodes.push_back(nd);
1465  _bounding_nodes[nd] = std::make_pair(_elem->node_ptr((ref_side+i)%n_nodes),
1466  _elem->node_ptr((ref_side+i+1)%n_nodes));
1467  side_p0 = side_nondim_points[(ref_side+i)%n_nodes].first;
1468  side_p1 = side_nondim_points[(ref_side+i)%n_nodes].second;
1469  _node_local_coords[nd] = side_p0 + xi_other * (side_p1 - side_p0);
1470  }
1471 
1472  // the algorithm of identifying the ref_side always ensures that the
1473  // region with xi > x0 on ref_side will be triangular. The other side
1474  // will be modeled with two QUAD4 elements
1475 
1476  libMesh::Elem
1477  *e_p = const_cast<libMesh::Elem*>(&e),
1478  *e1 = libMesh::Elem::build( libMesh::TRI3, e_p).release(),
1479  *e2 = libMesh::Elem::build( libMesh::TRI3, e_p).release(),
1480  *e3 = libMesh::Elem::build(libMesh::QUAD4, e_p).release(),
1481  *e4 = libMesh::Elem::build(libMesh::QUAD4, e_p).release();
1482 
1483  _new_elems.push_back(e1);
1484  _new_elems.push_back(e2);
1485  _new_elems.push_back(e3);
1486  _new_elems.push_back(e4);
1487 
1488  // create the elements. We create the elements such that:
1489  // node(0) = new node on ref_side
1490  // node(1) = elem node (ref_side+1)%n_sides
1491  // node(2) = new node on (ref_side+1)%n_sides
1492  // this way, the element has the same orientation as the original
1493  // element and the sides have the following association:
1494  // side(0) = coincident with ref_side of original element
1495  // side(1) = coincident side with ref_side+1 of original element
1496  // side(2) = coincident side with level-set interface
1497 
1498  e1->set_node(0) = _new_nodes[0];
1499  e1->set_node(1) = const_cast<libMesh::Node*>(e.node_ptr((ref_side+1)%n_nodes));
1500  e1->set_node(2) = _new_nodes[1];
1501 
1502  _elem_sides_on_interface[e1] = 2;
1503 
1504 
1505  // create the second element. We create the elements such that:
1506  // node(0) = new node on ref_side+3
1507  // node(1) = new node on (ref_side+2)%n_sides
1508  // node(2) = elem node (ref_side+3)%n_sides
1509  // this way, the element has the same orientation as the original
1510  // element and the sides have the following association:
1511  // side(0) = coincident with ref_side of original element
1512  // side(1) = coincident side with ref_side+1 of original element
1513  // side(2) = coincident side with level-set interface
1514 
1515  e2->set_node(0) = _new_nodes[3];
1516  e2->set_node(1) = _new_nodes[2];
1517  e2->set_node(2) = const_cast<libMesh::Node*>(e.node_ptr((ref_side+3)%n_nodes));
1518 
1519  _elem_sides_on_interface[e2] = 0;
1520 
1521 
1522  // create a third element and set its nodes so that:
1523  // node(0) = elem node ref_side
1524  // node(1) = new node on (ref_side)%n_nodes
1525  // node(2) = new node on (ref_side+2)%n_nodes
1526  // node(3) = new node on (ref_side+3)%n_nodes
1527  // this way, the element has the same orientation as the original
1528  // element and the sides have the following association:
1529  // side(0) = coincident side with ref_side of original element
1530  // side(1) = new edge connecting new node on ref_side to new node on ref_side+2
1531  // side(2) = coincident side with level-set interface
1532  // side(3) = coincident side with (ref_side+3)%n_sides of original element
1533 
1534  e3->set_node(0) = const_cast<libMesh::Node*>(e.node_ptr(ref_side));
1535  e3->set_node(1) = _new_nodes[0];
1536  e3->set_node(2) = _new_nodes[2];
1537  e3->set_node(3) = _new_nodes[3];
1538 
1539  _elem_sides_on_interface[e3] = 2;
1540 
1541  // create a fourth element and set its nodes so that:
1542  // node(0) = new node on ref_side
1543  // node(1) = new node on ref_side+1
1544  // node(2) = elem node ref_side+2
1545  // node(3) = new node on ref_side+2
1546  // this way, the element has the same orientation as the original
1547  // element and the sides have the following association:
1548  // side(0) = coincident side with level-set interface
1549  // side(1) = coincident side with (ref_side+1)%n_sides of original element
1550  // side(2) = coincident side with (ref_side+2)%n_sides of original element
1551  // side(3) = new edge connecting new node on ref_side to new node on ref_side+2
1552 
1553  e4->set_node(0) = _new_nodes[0];
1554  e4->set_node(1) = _new_nodes[1];
1555  e4->set_node(2) = const_cast<libMesh::Node*>(e.node_ptr((ref_side+2)%n_nodes));
1556  e4->set_node(3) = _new_nodes[2];
1557 
1558  _elem_sides_on_interface[e4] = 0;
1559 
1560  if (ref_side_node0_positive) {
1561 
1562  _positive_phi_elems.push_back(e3);
1563  _positive_phi_elems.push_back(e4);
1564  _negative_phi_elems.push_back(e1);
1565  _negative_phi_elems.push_back(e2);
1566  }
1567  else {
1568 
1569  _positive_phi_elems.push_back(e1);
1570  _positive_phi_elems.push_back(e2);
1571  _negative_phi_elems.push_back(e3);
1572  _negative_phi_elems.push_back(e4);
1573  }
1574  }
1575  break;
1576 
1577  default:
1578  // should not get here
1579  libmesh_error();
1580  }
1581 }
1582 
1583 
1584 
1585 Real
1587 (const libMesh::Point& p0,
1588  const libMesh::Point& p1,
1589  const MAST::FieldFunction<Real>& phi,
1590  const Real t) {
1591 
1592  // this uses the bisection search to find the location of intersection
1593  libMesh::Point
1594  pt,
1595  pt0 = p0,
1596  pt1 = p1;
1597 
1598  Real
1599  xi = 0.,
1600  v = 1.e10, // arbitrarily high value to get the algorithm going
1601  v0 = 0.,
1602  v1 = 0.;
1603 
1604  phi(pt0, t, v0);
1605  phi(pt1, t, v1);
1606 
1607  unsigned int
1608  n_iters = 0;
1609 
1610  //
1611  // a0 + a1 xi = 0 (a0 = v0, a1 = (v1-v0))
1612  // or, xi = -a0/a1
1613  //
1614 
1615  while (fabs(v) > _tol &&
1616  n_iters < _max_iters) {
1617 
1618  xi = -v0 / (v1-v0);
1619  pt = pt0 + (pt1 - pt0)*xi;
1620 
1621  phi(pt, t, v);
1622 
1623  if (v*v1 < 0.) {
1624 
1625  v0 = v;
1626  pt0 = pt;
1627  }
1628  else {
1629 
1630  v1 = v;
1631  pt1 = pt;
1632  }
1633 
1634  n_iters++;
1635  }
1636 
1637  // now find the xi location based on the distance of the new
1638  // point from the old points
1639  pt -= p0;
1640  xi = pt.norm();
1641  pt = p1-p0;
1642  xi /= pt.norm();
1643 
1644  return xi;
1645 }
1646 
1647 
1648 const libMesh::Point&
1650 (const libMesh::Node& n) const {
1651 
1652  libmesh_assert(_initialized);
1653 
1654  std::map<const libMesh::Node*, libMesh::Point>::const_iterator
1655  it = _node_local_coords.find(&n);
1656 
1657  libmesh_assert(it != _node_local_coords.end());
1658 
1659  return it->second;
1660 }
1661 
1662 
1663 
1664 bool
1665 MAST::LevelSetIntersection::if_node_is_new(const libMesh::Node& node) const {
1666 
1667  libmesh_assert(_initialized);
1668 
1669  for (unsigned int i=0; i<_new_nodes.size(); i++)
1670  if (&node == _new_nodes[i])
1671  return true;
1672 
1673  // if it gets here then the node was not created by this intersection
1674  // operation
1675  return false;
1676 }
1677 
1678 
1679 
1680 bool
1681 MAST::LevelSetIntersection::if_interior_node(const libMesh::Node& node) const {
1682 
1683  libmesh_assert(_initialized);
1684 
1685  return _interior_nodes.count(&node);
1686 }
1687 
1688 
1689 
1690 std::pair<const libMesh::Node*, const libMesh::Node*>
1692 
1693  libmesh_assert(_initialized);
1694 
1695  libmesh_assert(this->if_node_is_new(node));
1696 
1697  std::map<const libMesh::Node*, std::pair<const libMesh::Node*, const libMesh::Node*>>::const_iterator
1698  it = _bounding_nodes.find(&node);
1699 
1700  libmesh_assert(it != _bounding_nodes.end());
1701 
1702  return it->second;
1703 }
1704 
1705 
1706 void
1708 get_material_sides_without_intersection(std::set<unsigned int>& sides) const {
1709 
1710  libmesh_assert_equal_to(sides.size(), 0);
1711 
1712  std::map<const libMesh::Node*, std::pair<Real, bool> >::const_iterator
1713  it,
1714  end = _node_phi_vals.end();
1715 
1716  for (unsigned int i=0; i<_elem->n_sides(); i++) {
1717 
1718  std::unique_ptr<const libMesh::Elem> side(_elem->side_ptr(i).release());
1719 
1720  // check to see if all nodes of this side are on the material side
1721  bool
1722  on_material = true;
1723  for (unsigned int j=0; j<side->n_nodes(); j++) {
1724 
1725  it = _node_phi_vals.find(side->node_ptr(j));
1726  // the node level set value should be in the code
1727  libmesh_assert(it != end);
1728 
1729  // neither on level set nor on positive side
1730  if (!(it->second.second || it->second.first > _tol)) {
1731 
1732  on_material = false;
1733  break;
1734  }
1735  }
1736 
1737  // if the on_material flag is still true then the side is added to
1738  // the set
1739  if (on_material)
1740  sides.insert(i);
1741  }
1742 }
1743 
1744 
1745 void
1747 get_nearest_intersection_point(const libMesh::Point &p,
1748  libMesh::Point &pt) {
1749 
1750  libmesh_assert(_initialized);
1751 
1752  // check the element nodes for closest node
1753  std::map<const libMesh::Node*, std::pair<Real, bool> >::const_iterator
1754  it = _node_phi_vals.begin(),
1755  end = _node_phi_vals.end();
1756 
1757  libMesh::Point
1758  dp;
1759 
1760  Real
1761  v = 0.,
1762  dist = 1.e12;
1763 
1764  for ( ; it != end; it++) {
1765 
1766  if (it->second.second) { // on level set
1767 
1768  dp = *it->first - p;
1769 
1770  if (dp.norm() < dist) {
1771  pt = *it->first;
1772  dist = dp.norm();
1773  }
1774  }
1775  }
1776 
1777  // check the new nodes for possible candidates
1778  for (unsigned int i=0; i<_new_nodes.size(); i++) {
1779 
1780  (*_phi)(*_new_nodes[i], 0., v);
1781 
1782  if (std::fabs(v) <= _tol) {
1783 
1784  dp = *_new_nodes[i] - p;
1785 
1786  if (dp.norm() < dist) {
1787  pt = *_new_nodes[i];
1788  dist = dp.norm();
1789  }
1790  }
1791  }
1792 }
unsigned int edge_on_boundary() const
MAST::LevelSet2DIntersectionMode _mode
const MAST::FieldFunction< Real > * _phi
const std::vector< const libMesh::Elem * > & get_sub_elems_positive_phi() const
unsigned int get_side_on_interface(const libMesh::Elem &e) const
Real _find_intersection_on_straight_edge(const libMesh::Point &p0, const libMesh::Point &p1, const MAST::FieldFunction< Real > &phi, const Real t)
const libMesh::Elem & elem() const
void get_material_sides_without_intersection(std::set< unsigned int > &sides) const
identifies the sides of the element that are completely on the material side without any intersection...
MAST::LevelSet2DIntersectionMode get_intersection_mode() const
bool _if_elem_on_negative_phi
true if element is completely on the negative side of level set with no intersection ...
Real get_node_phi_value(const libMesh::Node *n) const
std::unique_ptr< libMesh::Elem > _first_order_elem(const libMesh::Elem &e)
creates a first order element from the given high-order element.
bool if_node_is_new(const libMesh::Node &node) const
identifies if the node from the subelements is a new node or an existing node from the parent element...
bool if_hanging_node(const libMesh::Node *n) const
The case of two adjacent edges results in a new node on an edge that is not coincident with the level...
bool _if_elem_on_positive_phi
true if element is completely on the positive side of level set with no intersection ...
std::map< const libMesh::Elem *, int > _elem_sides_on_interface
std::map< const libMesh::Node *, std::pair< const libMesh::Node *, const libMesh::Node * > > _bounding_nodes
std::pair< const libMesh::Node *, const libMesh::Node * > get_bounding_nodes_for_node(const libMesh::Node &node) const
for new nodes required to create the subelements this method returns the nodes on an edge that bound ...
void _add_node_local_coords(const libMesh::Elem &e, std::vector< std::pair< libMesh::Point, libMesh::Point > > &side_nondim_points, std::map< const libMesh::Node *, libMesh::Point > &node_coord_map)
void _find_quad4_intersections(const MAST::FieldFunction< Real > &phi, const libMesh::Elem &e, const Real t, const std::map< const libMesh::Node *, std::pair< Real, bool > > &node_phi_vals)
libMesh::Real Real
void init(const MAST::FieldFunction< Real > &phi, const libMesh::Elem &e, const Real t, unsigned int max_elem_id, unsigned int max_node_id)
std::vector< libMesh::Elem * > _new_elems
std::vector< const libMesh::Elem * > _positive_phi_elems
std::map< const libMesh::Node *, libMesh::Point > _node_local_coords
const Real _tol
std::vector< libMesh::Node * > _new_nodes
bool if_interior_node(const libMesh::Node &node) const
identifies if the new node is on an edge along the level-set method in the interior of the element (a...
bool has_side_on_interface(const libMesh::Elem &e) const
void get_corner_nodes_on_negative_phi(std::set< const libMesh::Node *> &nodes) const
unsigned int node_on_boundary() const
const std::vector< const libMesh::Elem * > & get_sub_elems_negative_phi() const
std::vector< const libMesh::Elem * > _negative_phi_elems
const libMesh::Point & get_nondimensional_coordinate_for_node(const libMesh::Node &n) const
std::set< const libMesh::Node * > _interior_nodes
std::map< const libMesh::Node *, std::pair< Real, bool > > _node_phi_vals
void _init_on_first_order_ref_elem(const MAST::FieldFunction< Real > &phi, const libMesh::Elem &e, const Real t)
initializes on a reference element that is a first-order counterpart of the given high-order element...
void get_nearest_intersection_point(const libMesh::Point &p, libMesh::Point &pt)
void clear()
clears the data structures
std::set< const libMesh::Node * > _hanging_node