diff --git a/Makefile.in b/Makefile.in index 710d9a8588f..9eb32461d94 100644 --- a/Makefile.in +++ b/Makefile.in @@ -337,22 +337,24 @@ am__libmesh_dbg_la_SOURCES_DIST = src/base/dirichlet_boundary.C \ src/fe/inf_fe_lagrange_eval.C src/fe/inf_fe_legendre_eval.C \ src/fe/inf_fe_map.C src/fe/inf_fe_map_eval.C \ src/fe/inf_fe_static.C src/geom/bounding_box.C src/geom/cell.C \ - src/geom/cell_hex.C src/geom/cell_hex20.C \ - src/geom/cell_hex27.C src/geom/cell_hex8.C src/geom/cell_inf.C \ + src/geom/cell_c0polyhedron.C src/geom/cell_hex.C \ + src/geom/cell_hex20.C src/geom/cell_hex27.C \ + src/geom/cell_hex8.C src/geom/cell_inf.C \ src/geom/cell_inf_hex.C src/geom/cell_inf_hex16.C \ src/geom/cell_inf_hex18.C src/geom/cell_inf_hex8.C \ src/geom/cell_inf_prism.C src/geom/cell_inf_prism12.C \ - src/geom/cell_inf_prism6.C src/geom/cell_prism.C \ - src/geom/cell_prism15.C src/geom/cell_prism18.C \ - src/geom/cell_prism20.C src/geom/cell_prism21.C \ - src/geom/cell_prism6.C src/geom/cell_pyramid.C \ - src/geom/cell_pyramid13.C src/geom/cell_pyramid14.C \ - src/geom/cell_pyramid18.C src/geom/cell_pyramid5.C \ - src/geom/cell_tet.C src/geom/cell_tet10.C \ - src/geom/cell_tet14.C src/geom/cell_tet4.C src/geom/edge.C \ - src/geom/edge_edge2.C src/geom/edge_edge3.C \ - src/geom/edge_edge4.C src/geom/edge_inf_edge2.C \ - src/geom/elem.C src/geom/elem_cutter.C src/geom/elem_quality.C \ + src/geom/cell_inf_prism6.C src/geom/cell_polyhedron.C \ + src/geom/cell_prism.C src/geom/cell_prism15.C \ + src/geom/cell_prism18.C src/geom/cell_prism20.C \ + src/geom/cell_prism21.C src/geom/cell_prism6.C \ + src/geom/cell_pyramid.C src/geom/cell_pyramid13.C \ + src/geom/cell_pyramid14.C src/geom/cell_pyramid18.C \ + src/geom/cell_pyramid5.C src/geom/cell_tet.C \ + src/geom/cell_tet10.C src/geom/cell_tet14.C \ + src/geom/cell_tet4.C src/geom/edge.C src/geom/edge_edge2.C \ + src/geom/edge_edge3.C src/geom/edge_edge4.C \ + src/geom/edge_inf_edge2.C src/geom/elem.C \ + src/geom/elem_cutter.C src/geom/elem_quality.C \ src/geom/elem_refinement.C src/geom/elem_side_builder.C \ src/geom/face.C src/geom/face_c0polygon.C \ src/geom/face_inf_quad.C src/geom/face_inf_quad4.C \ @@ -703,6 +705,7 @@ am__objects_1 = src/base/libmesh_dbg_la-dirichlet_boundary.lo \ src/fe/libmesh_dbg_la-inf_fe_static.lo \ src/geom/libmesh_dbg_la-bounding_box.lo \ src/geom/libmesh_dbg_la-cell.lo \ + src/geom/libmesh_dbg_la-cell_c0polyhedron.lo \ src/geom/libmesh_dbg_la-cell_hex.lo \ src/geom/libmesh_dbg_la-cell_hex20.lo \ src/geom/libmesh_dbg_la-cell_hex27.lo \ @@ -715,6 +718,7 @@ am__objects_1 = src/base/libmesh_dbg_la-dirichlet_boundary.lo \ src/geom/libmesh_dbg_la-cell_inf_prism.lo \ src/geom/libmesh_dbg_la-cell_inf_prism12.lo \ src/geom/libmesh_dbg_la-cell_inf_prism6.lo \ + src/geom/libmesh_dbg_la-cell_polyhedron.lo \ src/geom/libmesh_dbg_la-cell_prism.lo \ src/geom/libmesh_dbg_la-cell_prism15.lo \ src/geom/libmesh_dbg_la-cell_prism18.lo \ @@ -1133,22 +1137,24 @@ am__libmesh_devel_la_SOURCES_DIST = src/base/dirichlet_boundary.C \ src/fe/inf_fe_lagrange_eval.C src/fe/inf_fe_legendre_eval.C \ src/fe/inf_fe_map.C src/fe/inf_fe_map_eval.C \ src/fe/inf_fe_static.C src/geom/bounding_box.C src/geom/cell.C \ - src/geom/cell_hex.C src/geom/cell_hex20.C \ - src/geom/cell_hex27.C src/geom/cell_hex8.C src/geom/cell_inf.C \ + src/geom/cell_c0polyhedron.C src/geom/cell_hex.C \ + src/geom/cell_hex20.C src/geom/cell_hex27.C \ + src/geom/cell_hex8.C src/geom/cell_inf.C \ src/geom/cell_inf_hex.C src/geom/cell_inf_hex16.C \ src/geom/cell_inf_hex18.C src/geom/cell_inf_hex8.C \ src/geom/cell_inf_prism.C src/geom/cell_inf_prism12.C \ - src/geom/cell_inf_prism6.C src/geom/cell_prism.C \ - src/geom/cell_prism15.C src/geom/cell_prism18.C \ - src/geom/cell_prism20.C src/geom/cell_prism21.C \ - src/geom/cell_prism6.C src/geom/cell_pyramid.C \ - src/geom/cell_pyramid13.C src/geom/cell_pyramid14.C \ - src/geom/cell_pyramid18.C src/geom/cell_pyramid5.C \ - src/geom/cell_tet.C src/geom/cell_tet10.C \ - src/geom/cell_tet14.C src/geom/cell_tet4.C src/geom/edge.C \ - src/geom/edge_edge2.C src/geom/edge_edge3.C \ - src/geom/edge_edge4.C src/geom/edge_inf_edge2.C \ - src/geom/elem.C src/geom/elem_cutter.C src/geom/elem_quality.C \ + src/geom/cell_inf_prism6.C src/geom/cell_polyhedron.C \ + src/geom/cell_prism.C src/geom/cell_prism15.C \ + src/geom/cell_prism18.C src/geom/cell_prism20.C \ + src/geom/cell_prism21.C src/geom/cell_prism6.C \ + src/geom/cell_pyramid.C src/geom/cell_pyramid13.C \ + src/geom/cell_pyramid14.C src/geom/cell_pyramid18.C \ + src/geom/cell_pyramid5.C src/geom/cell_tet.C \ + src/geom/cell_tet10.C src/geom/cell_tet14.C \ + src/geom/cell_tet4.C src/geom/edge.C src/geom/edge_edge2.C \ + src/geom/edge_edge3.C src/geom/edge_edge4.C \ + src/geom/edge_inf_edge2.C src/geom/elem.C \ + src/geom/elem_cutter.C src/geom/elem_quality.C \ src/geom/elem_refinement.C src/geom/elem_side_builder.C \ src/geom/face.C src/geom/face_c0polygon.C \ src/geom/face_inf_quad.C src/geom/face_inf_quad4.C \ @@ -1498,6 +1504,7 @@ am__objects_2 = src/base/libmesh_devel_la-dirichlet_boundary.lo \ src/fe/libmesh_devel_la-inf_fe_static.lo \ src/geom/libmesh_devel_la-bounding_box.lo \ src/geom/libmesh_devel_la-cell.lo \ + src/geom/libmesh_devel_la-cell_c0polyhedron.lo \ src/geom/libmesh_devel_la-cell_hex.lo \ src/geom/libmesh_devel_la-cell_hex20.lo \ src/geom/libmesh_devel_la-cell_hex27.lo \ @@ -1510,6 +1517,7 @@ am__objects_2 = src/base/libmesh_devel_la-dirichlet_boundary.lo \ src/geom/libmesh_devel_la-cell_inf_prism.lo \ src/geom/libmesh_devel_la-cell_inf_prism12.lo \ src/geom/libmesh_devel_la-cell_inf_prism6.lo \ + src/geom/libmesh_devel_la-cell_polyhedron.lo \ src/geom/libmesh_devel_la-cell_prism.lo \ src/geom/libmesh_devel_la-cell_prism15.lo \ src/geom/libmesh_devel_la-cell_prism18.lo \ @@ -1925,22 +1933,24 @@ am__libmesh_oprof_la_SOURCES_DIST = src/base/dirichlet_boundary.C \ src/fe/inf_fe_lagrange_eval.C src/fe/inf_fe_legendre_eval.C \ src/fe/inf_fe_map.C src/fe/inf_fe_map_eval.C \ src/fe/inf_fe_static.C src/geom/bounding_box.C src/geom/cell.C \ - src/geom/cell_hex.C src/geom/cell_hex20.C \ - src/geom/cell_hex27.C src/geom/cell_hex8.C src/geom/cell_inf.C \ + src/geom/cell_c0polyhedron.C src/geom/cell_hex.C \ + src/geom/cell_hex20.C src/geom/cell_hex27.C \ + src/geom/cell_hex8.C src/geom/cell_inf.C \ src/geom/cell_inf_hex.C src/geom/cell_inf_hex16.C \ src/geom/cell_inf_hex18.C src/geom/cell_inf_hex8.C \ src/geom/cell_inf_prism.C src/geom/cell_inf_prism12.C \ - src/geom/cell_inf_prism6.C src/geom/cell_prism.C \ - src/geom/cell_prism15.C src/geom/cell_prism18.C \ - src/geom/cell_prism20.C src/geom/cell_prism21.C \ - src/geom/cell_prism6.C src/geom/cell_pyramid.C \ - src/geom/cell_pyramid13.C src/geom/cell_pyramid14.C \ - src/geom/cell_pyramid18.C src/geom/cell_pyramid5.C \ - src/geom/cell_tet.C src/geom/cell_tet10.C \ - src/geom/cell_tet14.C src/geom/cell_tet4.C src/geom/edge.C \ - src/geom/edge_edge2.C src/geom/edge_edge3.C \ - src/geom/edge_edge4.C src/geom/edge_inf_edge2.C \ - src/geom/elem.C src/geom/elem_cutter.C src/geom/elem_quality.C \ + src/geom/cell_inf_prism6.C src/geom/cell_polyhedron.C \ + src/geom/cell_prism.C src/geom/cell_prism15.C \ + src/geom/cell_prism18.C src/geom/cell_prism20.C \ + src/geom/cell_prism21.C src/geom/cell_prism6.C \ + src/geom/cell_pyramid.C src/geom/cell_pyramid13.C \ + src/geom/cell_pyramid14.C src/geom/cell_pyramid18.C \ + src/geom/cell_pyramid5.C src/geom/cell_tet.C \ + src/geom/cell_tet10.C src/geom/cell_tet14.C \ + src/geom/cell_tet4.C src/geom/edge.C src/geom/edge_edge2.C \ + src/geom/edge_edge3.C src/geom/edge_edge4.C \ + src/geom/edge_inf_edge2.C src/geom/elem.C \ + src/geom/elem_cutter.C src/geom/elem_quality.C \ src/geom/elem_refinement.C src/geom/elem_side_builder.C \ src/geom/face.C src/geom/face_c0polygon.C \ src/geom/face_inf_quad.C src/geom/face_inf_quad4.C \ @@ -2290,6 +2300,7 @@ am__objects_3 = src/base/libmesh_oprof_la-dirichlet_boundary.lo \ src/fe/libmesh_oprof_la-inf_fe_static.lo \ src/geom/libmesh_oprof_la-bounding_box.lo \ src/geom/libmesh_oprof_la-cell.lo \ + src/geom/libmesh_oprof_la-cell_c0polyhedron.lo \ src/geom/libmesh_oprof_la-cell_hex.lo \ src/geom/libmesh_oprof_la-cell_hex20.lo \ src/geom/libmesh_oprof_la-cell_hex27.lo \ @@ -2302,6 +2313,7 @@ am__objects_3 = src/base/libmesh_oprof_la-dirichlet_boundary.lo \ src/geom/libmesh_oprof_la-cell_inf_prism.lo \ src/geom/libmesh_oprof_la-cell_inf_prism12.lo \ src/geom/libmesh_oprof_la-cell_inf_prism6.lo \ + src/geom/libmesh_oprof_la-cell_polyhedron.lo \ src/geom/libmesh_oprof_la-cell_prism.lo \ src/geom/libmesh_oprof_la-cell_prism15.lo \ src/geom/libmesh_oprof_la-cell_prism18.lo \ @@ -2717,22 +2729,24 @@ am__libmesh_opt_la_SOURCES_DIST = src/base/dirichlet_boundary.C \ src/fe/inf_fe_lagrange_eval.C src/fe/inf_fe_legendre_eval.C \ src/fe/inf_fe_map.C src/fe/inf_fe_map_eval.C \ src/fe/inf_fe_static.C src/geom/bounding_box.C src/geom/cell.C \ - src/geom/cell_hex.C src/geom/cell_hex20.C \ - src/geom/cell_hex27.C src/geom/cell_hex8.C src/geom/cell_inf.C \ + src/geom/cell_c0polyhedron.C src/geom/cell_hex.C \ + src/geom/cell_hex20.C src/geom/cell_hex27.C \ + src/geom/cell_hex8.C src/geom/cell_inf.C \ src/geom/cell_inf_hex.C src/geom/cell_inf_hex16.C \ src/geom/cell_inf_hex18.C src/geom/cell_inf_hex8.C \ src/geom/cell_inf_prism.C src/geom/cell_inf_prism12.C \ - src/geom/cell_inf_prism6.C src/geom/cell_prism.C \ - src/geom/cell_prism15.C src/geom/cell_prism18.C \ - src/geom/cell_prism20.C src/geom/cell_prism21.C \ - src/geom/cell_prism6.C src/geom/cell_pyramid.C \ - src/geom/cell_pyramid13.C src/geom/cell_pyramid14.C \ - src/geom/cell_pyramid18.C src/geom/cell_pyramid5.C \ - src/geom/cell_tet.C src/geom/cell_tet10.C \ - src/geom/cell_tet14.C src/geom/cell_tet4.C src/geom/edge.C \ - src/geom/edge_edge2.C src/geom/edge_edge3.C \ - src/geom/edge_edge4.C src/geom/edge_inf_edge2.C \ - src/geom/elem.C src/geom/elem_cutter.C src/geom/elem_quality.C \ + src/geom/cell_inf_prism6.C src/geom/cell_polyhedron.C \ + src/geom/cell_prism.C src/geom/cell_prism15.C \ + src/geom/cell_prism18.C src/geom/cell_prism20.C \ + src/geom/cell_prism21.C src/geom/cell_prism6.C \ + src/geom/cell_pyramid.C src/geom/cell_pyramid13.C \ + src/geom/cell_pyramid14.C src/geom/cell_pyramid18.C \ + src/geom/cell_pyramid5.C src/geom/cell_tet.C \ + src/geom/cell_tet10.C src/geom/cell_tet14.C \ + src/geom/cell_tet4.C src/geom/edge.C src/geom/edge_edge2.C \ + src/geom/edge_edge3.C src/geom/edge_edge4.C \ + src/geom/edge_inf_edge2.C src/geom/elem.C \ + src/geom/elem_cutter.C src/geom/elem_quality.C \ src/geom/elem_refinement.C src/geom/elem_side_builder.C \ src/geom/face.C src/geom/face_c0polygon.C \ src/geom/face_inf_quad.C src/geom/face_inf_quad4.C \ @@ -3082,6 +3096,7 @@ am__objects_4 = src/base/libmesh_opt_la-dirichlet_boundary.lo \ src/fe/libmesh_opt_la-inf_fe_static.lo \ src/geom/libmesh_opt_la-bounding_box.lo \ src/geom/libmesh_opt_la-cell.lo \ + src/geom/libmesh_opt_la-cell_c0polyhedron.lo \ src/geom/libmesh_opt_la-cell_hex.lo \ src/geom/libmesh_opt_la-cell_hex20.lo \ src/geom/libmesh_opt_la-cell_hex27.lo \ @@ -3094,6 +3109,7 @@ am__objects_4 = src/base/libmesh_opt_la-dirichlet_boundary.lo \ src/geom/libmesh_opt_la-cell_inf_prism.lo \ src/geom/libmesh_opt_la-cell_inf_prism12.lo \ src/geom/libmesh_opt_la-cell_inf_prism6.lo \ + src/geom/libmesh_opt_la-cell_polyhedron.lo \ src/geom/libmesh_opt_la-cell_prism.lo \ src/geom/libmesh_opt_la-cell_prism15.lo \ src/geom/libmesh_opt_la-cell_prism18.lo \ @@ -3508,22 +3524,24 @@ am__libmesh_prof_la_SOURCES_DIST = src/base/dirichlet_boundary.C \ src/fe/inf_fe_lagrange_eval.C src/fe/inf_fe_legendre_eval.C \ src/fe/inf_fe_map.C src/fe/inf_fe_map_eval.C \ src/fe/inf_fe_static.C src/geom/bounding_box.C src/geom/cell.C \ - src/geom/cell_hex.C src/geom/cell_hex20.C \ - src/geom/cell_hex27.C src/geom/cell_hex8.C src/geom/cell_inf.C \ + src/geom/cell_c0polyhedron.C src/geom/cell_hex.C \ + src/geom/cell_hex20.C src/geom/cell_hex27.C \ + src/geom/cell_hex8.C src/geom/cell_inf.C \ src/geom/cell_inf_hex.C src/geom/cell_inf_hex16.C \ src/geom/cell_inf_hex18.C src/geom/cell_inf_hex8.C \ src/geom/cell_inf_prism.C src/geom/cell_inf_prism12.C \ - src/geom/cell_inf_prism6.C src/geom/cell_prism.C \ - src/geom/cell_prism15.C src/geom/cell_prism18.C \ - src/geom/cell_prism20.C src/geom/cell_prism21.C \ - src/geom/cell_prism6.C src/geom/cell_pyramid.C \ - src/geom/cell_pyramid13.C src/geom/cell_pyramid14.C \ - src/geom/cell_pyramid18.C src/geom/cell_pyramid5.C \ - src/geom/cell_tet.C src/geom/cell_tet10.C \ - src/geom/cell_tet14.C src/geom/cell_tet4.C src/geom/edge.C \ - src/geom/edge_edge2.C src/geom/edge_edge3.C \ - src/geom/edge_edge4.C src/geom/edge_inf_edge2.C \ - src/geom/elem.C src/geom/elem_cutter.C src/geom/elem_quality.C \ + src/geom/cell_inf_prism6.C src/geom/cell_polyhedron.C \ + src/geom/cell_prism.C src/geom/cell_prism15.C \ + src/geom/cell_prism18.C src/geom/cell_prism20.C \ + src/geom/cell_prism21.C src/geom/cell_prism6.C \ + src/geom/cell_pyramid.C src/geom/cell_pyramid13.C \ + src/geom/cell_pyramid14.C src/geom/cell_pyramid18.C \ + src/geom/cell_pyramid5.C src/geom/cell_tet.C \ + src/geom/cell_tet10.C src/geom/cell_tet14.C \ + src/geom/cell_tet4.C src/geom/edge.C src/geom/edge_edge2.C \ + src/geom/edge_edge3.C src/geom/edge_edge4.C \ + src/geom/edge_inf_edge2.C src/geom/elem.C \ + src/geom/elem_cutter.C src/geom/elem_quality.C \ src/geom/elem_refinement.C src/geom/elem_side_builder.C \ src/geom/face.C src/geom/face_c0polygon.C \ src/geom/face_inf_quad.C src/geom/face_inf_quad4.C \ @@ -3873,6 +3891,7 @@ am__objects_5 = src/base/libmesh_prof_la-dirichlet_boundary.lo \ src/fe/libmesh_prof_la-inf_fe_static.lo \ src/geom/libmesh_prof_la-bounding_box.lo \ src/geom/libmesh_prof_la-cell.lo \ + src/geom/libmesh_prof_la-cell_c0polyhedron.lo \ src/geom/libmesh_prof_la-cell_hex.lo \ src/geom/libmesh_prof_la-cell_hex20.lo \ src/geom/libmesh_prof_la-cell_hex27.lo \ @@ -3885,6 +3904,7 @@ am__objects_5 = src/base/libmesh_prof_la-dirichlet_boundary.lo \ src/geom/libmesh_prof_la-cell_inf_prism.lo \ src/geom/libmesh_prof_la-cell_inf_prism12.lo \ src/geom/libmesh_prof_la-cell_inf_prism6.lo \ + src/geom/libmesh_prof_la-cell_polyhedron.lo \ src/geom/libmesh_prof_la-cell_prism.lo \ src/geom/libmesh_prof_la-cell_prism15.lo \ src/geom/libmesh_prof_la-cell_prism18.lo \ @@ -5283,6 +5303,7 @@ am__depfiles_remade = src/apps/$(DEPDIR)/amr_dbg-amr.Po \ src/fe/$(DEPDIR)/libmesh_prof_la-inf_fe_static.Plo \ src/geom/$(DEPDIR)/libmesh_dbg_la-bounding_box.Plo \ src/geom/$(DEPDIR)/libmesh_dbg_la-cell.Plo \ + src/geom/$(DEPDIR)/libmesh_dbg_la-cell_c0polyhedron.Plo \ src/geom/$(DEPDIR)/libmesh_dbg_la-cell_hex.Plo \ src/geom/$(DEPDIR)/libmesh_dbg_la-cell_hex20.Plo \ src/geom/$(DEPDIR)/libmesh_dbg_la-cell_hex27.Plo \ @@ -5295,6 +5316,7 @@ am__depfiles_remade = src/apps/$(DEPDIR)/amr_dbg-amr.Po \ src/geom/$(DEPDIR)/libmesh_dbg_la-cell_inf_prism.Plo \ src/geom/$(DEPDIR)/libmesh_dbg_la-cell_inf_prism12.Plo \ src/geom/$(DEPDIR)/libmesh_dbg_la-cell_inf_prism6.Plo \ + src/geom/$(DEPDIR)/libmesh_dbg_la-cell_polyhedron.Plo \ src/geom/$(DEPDIR)/libmesh_dbg_la-cell_prism.Plo \ src/geom/$(DEPDIR)/libmesh_dbg_la-cell_prism15.Plo \ src/geom/$(DEPDIR)/libmesh_dbg_la-cell_prism18.Plo \ @@ -5349,6 +5371,7 @@ am__depfiles_remade = src/apps/$(DEPDIR)/amr_dbg-amr.Po \ src/geom/$(DEPDIR)/libmesh_dbg_la-surface.Plo \ src/geom/$(DEPDIR)/libmesh_devel_la-bounding_box.Plo \ src/geom/$(DEPDIR)/libmesh_devel_la-cell.Plo \ + src/geom/$(DEPDIR)/libmesh_devel_la-cell_c0polyhedron.Plo \ src/geom/$(DEPDIR)/libmesh_devel_la-cell_hex.Plo \ src/geom/$(DEPDIR)/libmesh_devel_la-cell_hex20.Plo \ src/geom/$(DEPDIR)/libmesh_devel_la-cell_hex27.Plo \ @@ -5361,6 +5384,7 @@ am__depfiles_remade = src/apps/$(DEPDIR)/amr_dbg-amr.Po \ src/geom/$(DEPDIR)/libmesh_devel_la-cell_inf_prism.Plo \ src/geom/$(DEPDIR)/libmesh_devel_la-cell_inf_prism12.Plo \ src/geom/$(DEPDIR)/libmesh_devel_la-cell_inf_prism6.Plo \ + src/geom/$(DEPDIR)/libmesh_devel_la-cell_polyhedron.Plo \ src/geom/$(DEPDIR)/libmesh_devel_la-cell_prism.Plo \ src/geom/$(DEPDIR)/libmesh_devel_la-cell_prism15.Plo \ src/geom/$(DEPDIR)/libmesh_devel_la-cell_prism18.Plo \ @@ -5415,6 +5439,7 @@ am__depfiles_remade = src/apps/$(DEPDIR)/amr_dbg-amr.Po \ src/geom/$(DEPDIR)/libmesh_devel_la-surface.Plo \ src/geom/$(DEPDIR)/libmesh_oprof_la-bounding_box.Plo \ src/geom/$(DEPDIR)/libmesh_oprof_la-cell.Plo \ + src/geom/$(DEPDIR)/libmesh_oprof_la-cell_c0polyhedron.Plo \ src/geom/$(DEPDIR)/libmesh_oprof_la-cell_hex.Plo \ src/geom/$(DEPDIR)/libmesh_oprof_la-cell_hex20.Plo \ src/geom/$(DEPDIR)/libmesh_oprof_la-cell_hex27.Plo \ @@ -5427,6 +5452,7 @@ am__depfiles_remade = src/apps/$(DEPDIR)/amr_dbg-amr.Po \ src/geom/$(DEPDIR)/libmesh_oprof_la-cell_inf_prism.Plo \ src/geom/$(DEPDIR)/libmesh_oprof_la-cell_inf_prism12.Plo \ src/geom/$(DEPDIR)/libmesh_oprof_la-cell_inf_prism6.Plo \ + src/geom/$(DEPDIR)/libmesh_oprof_la-cell_polyhedron.Plo \ src/geom/$(DEPDIR)/libmesh_oprof_la-cell_prism.Plo \ src/geom/$(DEPDIR)/libmesh_oprof_la-cell_prism15.Plo \ src/geom/$(DEPDIR)/libmesh_oprof_la-cell_prism18.Plo \ @@ -5481,6 +5507,7 @@ am__depfiles_remade = src/apps/$(DEPDIR)/amr_dbg-amr.Po \ src/geom/$(DEPDIR)/libmesh_oprof_la-surface.Plo \ src/geom/$(DEPDIR)/libmesh_opt_la-bounding_box.Plo \ src/geom/$(DEPDIR)/libmesh_opt_la-cell.Plo \ + src/geom/$(DEPDIR)/libmesh_opt_la-cell_c0polyhedron.Plo \ src/geom/$(DEPDIR)/libmesh_opt_la-cell_hex.Plo \ src/geom/$(DEPDIR)/libmesh_opt_la-cell_hex20.Plo \ src/geom/$(DEPDIR)/libmesh_opt_la-cell_hex27.Plo \ @@ -5493,6 +5520,7 @@ am__depfiles_remade = src/apps/$(DEPDIR)/amr_dbg-amr.Po \ src/geom/$(DEPDIR)/libmesh_opt_la-cell_inf_prism.Plo \ src/geom/$(DEPDIR)/libmesh_opt_la-cell_inf_prism12.Plo \ src/geom/$(DEPDIR)/libmesh_opt_la-cell_inf_prism6.Plo \ + src/geom/$(DEPDIR)/libmesh_opt_la-cell_polyhedron.Plo \ src/geom/$(DEPDIR)/libmesh_opt_la-cell_prism.Plo \ src/geom/$(DEPDIR)/libmesh_opt_la-cell_prism15.Plo \ src/geom/$(DEPDIR)/libmesh_opt_la-cell_prism18.Plo \ @@ -5547,6 +5575,7 @@ am__depfiles_remade = src/apps/$(DEPDIR)/amr_dbg-amr.Po \ src/geom/$(DEPDIR)/libmesh_opt_la-surface.Plo \ src/geom/$(DEPDIR)/libmesh_prof_la-bounding_box.Plo \ src/geom/$(DEPDIR)/libmesh_prof_la-cell.Plo \ + src/geom/$(DEPDIR)/libmesh_prof_la-cell_c0polyhedron.Plo \ src/geom/$(DEPDIR)/libmesh_prof_la-cell_hex.Plo \ src/geom/$(DEPDIR)/libmesh_prof_la-cell_hex20.Plo \ src/geom/$(DEPDIR)/libmesh_prof_la-cell_hex27.Plo \ @@ -5559,6 +5588,7 @@ am__depfiles_remade = src/apps/$(DEPDIR)/amr_dbg-amr.Po \ src/geom/$(DEPDIR)/libmesh_prof_la-cell_inf_prism.Plo \ src/geom/$(DEPDIR)/libmesh_prof_la-cell_inf_prism12.Plo \ src/geom/$(DEPDIR)/libmesh_prof_la-cell_inf_prism6.Plo \ + src/geom/$(DEPDIR)/libmesh_prof_la-cell_polyhedron.Plo \ src/geom/$(DEPDIR)/libmesh_prof_la-cell_prism.Plo \ src/geom/$(DEPDIR)/libmesh_prof_la-cell_prism15.Plo \ src/geom/$(DEPDIR)/libmesh_prof_la-cell_prism18.Plo \ @@ -7797,6 +7827,7 @@ libmesh_SOURCES = \ src/fe/inf_fe_static.C \ src/geom/bounding_box.C \ src/geom/cell.C \ + src/geom/cell_c0polyhedron.C \ src/geom/cell_hex.C \ src/geom/cell_hex20.C \ src/geom/cell_hex27.C \ @@ -7809,6 +7840,7 @@ libmesh_SOURCES = \ src/geom/cell_inf_prism.C \ src/geom/cell_inf_prism12.C \ src/geom/cell_inf_prism6.C \ + src/geom/cell_polyhedron.C \ src/geom/cell_prism.C \ src/geom/cell_prism15.C \ src/geom/cell_prism18.C \ @@ -8882,6 +8914,8 @@ src/geom/libmesh_dbg_la-bounding_box.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_dbg_la-cell.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) +src/geom/libmesh_dbg_la-cell_c0polyhedron.lo: \ + src/geom/$(am__dirstamp) src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_dbg_la-cell_hex.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_dbg_la-cell_hex20.lo: src/geom/$(am__dirstamp) \ @@ -8906,6 +8940,8 @@ src/geom/libmesh_dbg_la-cell_inf_prism12.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_dbg_la-cell_inf_prism6.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) +src/geom/libmesh_dbg_la-cell_polyhedron.lo: src/geom/$(am__dirstamp) \ + src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_dbg_la-cell_prism.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_dbg_la-cell_prism15.lo: src/geom/$(am__dirstamp) \ @@ -10101,6 +10137,8 @@ src/geom/libmesh_devel_la-bounding_box.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_devel_la-cell.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) +src/geom/libmesh_devel_la-cell_c0polyhedron.lo: \ + src/geom/$(am__dirstamp) src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_devel_la-cell_hex.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_devel_la-cell_hex20.lo: src/geom/$(am__dirstamp) \ @@ -10125,6 +10163,8 @@ src/geom/libmesh_devel_la-cell_inf_prism12.lo: \ src/geom/$(am__dirstamp) src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_devel_la-cell_inf_prism6.lo: \ src/geom/$(am__dirstamp) src/geom/$(DEPDIR)/$(am__dirstamp) +src/geom/libmesh_devel_la-cell_polyhedron.lo: \ + src/geom/$(am__dirstamp) src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_devel_la-cell_prism.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_devel_la-cell_prism15.lo: src/geom/$(am__dirstamp) \ @@ -11251,6 +11291,8 @@ src/geom/libmesh_oprof_la-bounding_box.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_oprof_la-cell.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) +src/geom/libmesh_oprof_la-cell_c0polyhedron.lo: \ + src/geom/$(am__dirstamp) src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_oprof_la-cell_hex.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_oprof_la-cell_hex20.lo: src/geom/$(am__dirstamp) \ @@ -11275,6 +11317,8 @@ src/geom/libmesh_oprof_la-cell_inf_prism12.lo: \ src/geom/$(am__dirstamp) src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_oprof_la-cell_inf_prism6.lo: \ src/geom/$(am__dirstamp) src/geom/$(DEPDIR)/$(am__dirstamp) +src/geom/libmesh_oprof_la-cell_polyhedron.lo: \ + src/geom/$(am__dirstamp) src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_oprof_la-cell_prism.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_oprof_la-cell_prism15.lo: src/geom/$(am__dirstamp) \ @@ -12401,6 +12445,8 @@ src/geom/libmesh_opt_la-bounding_box.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_opt_la-cell.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) +src/geom/libmesh_opt_la-cell_c0polyhedron.lo: \ + src/geom/$(am__dirstamp) src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_opt_la-cell_hex.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_opt_la-cell_hex20.lo: src/geom/$(am__dirstamp) \ @@ -12425,6 +12471,8 @@ src/geom/libmesh_opt_la-cell_inf_prism12.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_opt_la-cell_inf_prism6.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) +src/geom/libmesh_opt_la-cell_polyhedron.lo: src/geom/$(am__dirstamp) \ + src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_opt_la-cell_prism.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_opt_la-cell_prism15.lo: src/geom/$(am__dirstamp) \ @@ -13548,6 +13596,8 @@ src/geom/libmesh_prof_la-bounding_box.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_prof_la-cell.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) +src/geom/libmesh_prof_la-cell_c0polyhedron.lo: \ + src/geom/$(am__dirstamp) src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_prof_la-cell_hex.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_prof_la-cell_hex20.lo: src/geom/$(am__dirstamp) \ @@ -13572,6 +13622,8 @@ src/geom/libmesh_prof_la-cell_inf_prism12.lo: \ src/geom/$(am__dirstamp) src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_prof_la-cell_inf_prism6.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) +src/geom/libmesh_prof_la-cell_polyhedron.lo: src/geom/$(am__dirstamp) \ + src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_prof_la-cell_prism.lo: src/geom/$(am__dirstamp) \ src/geom/$(DEPDIR)/$(am__dirstamp) src/geom/libmesh_prof_la-cell_prism15.lo: src/geom/$(am__dirstamp) \ @@ -15574,6 +15626,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@src/fe/$(DEPDIR)/libmesh_prof_la-inf_fe_static.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_dbg_la-bounding_box.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_dbg_la-cell.Plo@am__quote@ # am--include-marker +@AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_dbg_la-cell_c0polyhedron.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_dbg_la-cell_hex.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_dbg_la-cell_hex20.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_dbg_la-cell_hex27.Plo@am__quote@ # am--include-marker @@ -15586,6 +15639,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_dbg_la-cell_inf_prism.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_dbg_la-cell_inf_prism12.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_dbg_la-cell_inf_prism6.Plo@am__quote@ # am--include-marker +@AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_dbg_la-cell_polyhedron.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_dbg_la-cell_prism.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_dbg_la-cell_prism15.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_dbg_la-cell_prism18.Plo@am__quote@ # am--include-marker @@ -15640,6 +15694,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_dbg_la-surface.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_devel_la-bounding_box.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_devel_la-cell.Plo@am__quote@ # am--include-marker +@AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_devel_la-cell_c0polyhedron.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_devel_la-cell_hex.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_devel_la-cell_hex20.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_devel_la-cell_hex27.Plo@am__quote@ # am--include-marker @@ -15652,6 +15707,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_devel_la-cell_inf_prism.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_devel_la-cell_inf_prism12.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_devel_la-cell_inf_prism6.Plo@am__quote@ # am--include-marker +@AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_devel_la-cell_polyhedron.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_devel_la-cell_prism.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_devel_la-cell_prism15.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_devel_la-cell_prism18.Plo@am__quote@ # am--include-marker @@ -15706,6 +15762,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_devel_la-surface.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_oprof_la-bounding_box.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_oprof_la-cell.Plo@am__quote@ # am--include-marker +@AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_oprof_la-cell_c0polyhedron.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_oprof_la-cell_hex.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_oprof_la-cell_hex20.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_oprof_la-cell_hex27.Plo@am__quote@ # am--include-marker @@ -15718,6 +15775,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_oprof_la-cell_inf_prism.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_oprof_la-cell_inf_prism12.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_oprof_la-cell_inf_prism6.Plo@am__quote@ # am--include-marker +@AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_oprof_la-cell_polyhedron.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_oprof_la-cell_prism.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_oprof_la-cell_prism15.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_oprof_la-cell_prism18.Plo@am__quote@ # am--include-marker @@ -15772,6 +15830,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_oprof_la-surface.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_opt_la-bounding_box.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_opt_la-cell.Plo@am__quote@ # am--include-marker +@AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_opt_la-cell_c0polyhedron.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_opt_la-cell_hex.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_opt_la-cell_hex20.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_opt_la-cell_hex27.Plo@am__quote@ # am--include-marker @@ -15784,6 +15843,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_opt_la-cell_inf_prism.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_opt_la-cell_inf_prism12.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_opt_la-cell_inf_prism6.Plo@am__quote@ # am--include-marker +@AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_opt_la-cell_polyhedron.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_opt_la-cell_prism.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_opt_la-cell_prism15.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_opt_la-cell_prism18.Plo@am__quote@ # am--include-marker @@ -15838,6 +15898,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_opt_la-surface.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_prof_la-bounding_box.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_prof_la-cell.Plo@am__quote@ # am--include-marker +@AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_prof_la-cell_c0polyhedron.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_prof_la-cell_hex.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_prof_la-cell_hex20.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_prof_la-cell_hex27.Plo@am__quote@ # am--include-marker @@ -15850,6 +15911,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_prof_la-cell_inf_prism.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_prof_la-cell_inf_prism12.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_prof_la-cell_inf_prism6.Plo@am__quote@ # am--include-marker +@AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_prof_la-cell_polyhedron.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_prof_la-cell_prism.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_prof_la-cell_prism15.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_prof_la-cell_prism18.Plo@am__quote@ # am--include-marker @@ -18185,6 +18247,13 @@ src/geom/libmesh_dbg_la-cell.lo: src/geom/cell.C @AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ @am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_dbg_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_dbg_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_dbg_la-cell.lo `test -f 'src/geom/cell.C' || echo '$(srcdir)/'`src/geom/cell.C +src/geom/libmesh_dbg_la-cell_c0polyhedron.lo: src/geom/cell_c0polyhedron.C +@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_dbg_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_dbg_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_dbg_la-cell_c0polyhedron.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_dbg_la-cell_c0polyhedron.Tpo -c -o src/geom/libmesh_dbg_la-cell_c0polyhedron.lo `test -f 'src/geom/cell_c0polyhedron.C' || echo '$(srcdir)/'`src/geom/cell_c0polyhedron.C +@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_dbg_la-cell_c0polyhedron.Tpo src/geom/$(DEPDIR)/libmesh_dbg_la-cell_c0polyhedron.Plo +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ $(AM_V_CXX)source='src/geom/cell_c0polyhedron.C' object='src/geom/libmesh_dbg_la-cell_c0polyhedron.lo' libtool=yes @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_dbg_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_dbg_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_dbg_la-cell_c0polyhedron.lo `test -f 'src/geom/cell_c0polyhedron.C' || echo '$(srcdir)/'`src/geom/cell_c0polyhedron.C + src/geom/libmesh_dbg_la-cell_hex.lo: src/geom/cell_hex.C @am__fastdepCXX_TRUE@ $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_dbg_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_dbg_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_dbg_la-cell_hex.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_dbg_la-cell_hex.Tpo -c -o src/geom/libmesh_dbg_la-cell_hex.lo `test -f 'src/geom/cell_hex.C' || echo '$(srcdir)/'`src/geom/cell_hex.C @am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_dbg_la-cell_hex.Tpo src/geom/$(DEPDIR)/libmesh_dbg_la-cell_hex.Plo @@ -18269,6 +18338,13 @@ src/geom/libmesh_dbg_la-cell_inf_prism6.lo: src/geom/cell_inf_prism6.C @AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ @am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_dbg_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_dbg_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_dbg_la-cell_inf_prism6.lo `test -f 'src/geom/cell_inf_prism6.C' || echo '$(srcdir)/'`src/geom/cell_inf_prism6.C +src/geom/libmesh_dbg_la-cell_polyhedron.lo: src/geom/cell_polyhedron.C +@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_dbg_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_dbg_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_dbg_la-cell_polyhedron.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_dbg_la-cell_polyhedron.Tpo -c -o src/geom/libmesh_dbg_la-cell_polyhedron.lo `test -f 'src/geom/cell_polyhedron.C' || echo '$(srcdir)/'`src/geom/cell_polyhedron.C +@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_dbg_la-cell_polyhedron.Tpo src/geom/$(DEPDIR)/libmesh_dbg_la-cell_polyhedron.Plo +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ $(AM_V_CXX)source='src/geom/cell_polyhedron.C' object='src/geom/libmesh_dbg_la-cell_polyhedron.lo' libtool=yes @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_dbg_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_dbg_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_dbg_la-cell_polyhedron.lo `test -f 'src/geom/cell_polyhedron.C' || echo '$(srcdir)/'`src/geom/cell_polyhedron.C + src/geom/libmesh_dbg_la-cell_prism.lo: src/geom/cell_prism.C @am__fastdepCXX_TRUE@ $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_dbg_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_dbg_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_dbg_la-cell_prism.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_dbg_la-cell_prism.Tpo -c -o src/geom/libmesh_dbg_la-cell_prism.lo `test -f 'src/geom/cell_prism.C' || echo '$(srcdir)/'`src/geom/cell_prism.C @am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_dbg_la-cell_prism.Tpo src/geom/$(DEPDIR)/libmesh_dbg_la-cell_prism.Plo @@ -21447,6 +21523,13 @@ src/geom/libmesh_devel_la-cell.lo: src/geom/cell.C @AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ @am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_devel_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_devel_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_devel_la-cell.lo `test -f 'src/geom/cell.C' || echo '$(srcdir)/'`src/geom/cell.C +src/geom/libmesh_devel_la-cell_c0polyhedron.lo: src/geom/cell_c0polyhedron.C +@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_devel_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_devel_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_devel_la-cell_c0polyhedron.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_devel_la-cell_c0polyhedron.Tpo -c -o src/geom/libmesh_devel_la-cell_c0polyhedron.lo `test -f 'src/geom/cell_c0polyhedron.C' || echo '$(srcdir)/'`src/geom/cell_c0polyhedron.C +@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_devel_la-cell_c0polyhedron.Tpo src/geom/$(DEPDIR)/libmesh_devel_la-cell_c0polyhedron.Plo +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ $(AM_V_CXX)source='src/geom/cell_c0polyhedron.C' object='src/geom/libmesh_devel_la-cell_c0polyhedron.lo' libtool=yes @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_devel_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_devel_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_devel_la-cell_c0polyhedron.lo `test -f 'src/geom/cell_c0polyhedron.C' || echo '$(srcdir)/'`src/geom/cell_c0polyhedron.C + src/geom/libmesh_devel_la-cell_hex.lo: src/geom/cell_hex.C @am__fastdepCXX_TRUE@ $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_devel_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_devel_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_devel_la-cell_hex.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_devel_la-cell_hex.Tpo -c -o src/geom/libmesh_devel_la-cell_hex.lo `test -f 'src/geom/cell_hex.C' || echo '$(srcdir)/'`src/geom/cell_hex.C @am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_devel_la-cell_hex.Tpo src/geom/$(DEPDIR)/libmesh_devel_la-cell_hex.Plo @@ -21531,6 +21614,13 @@ src/geom/libmesh_devel_la-cell_inf_prism6.lo: src/geom/cell_inf_prism6.C @AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ @am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_devel_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_devel_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_devel_la-cell_inf_prism6.lo `test -f 'src/geom/cell_inf_prism6.C' || echo '$(srcdir)/'`src/geom/cell_inf_prism6.C +src/geom/libmesh_devel_la-cell_polyhedron.lo: src/geom/cell_polyhedron.C +@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_devel_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_devel_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_devel_la-cell_polyhedron.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_devel_la-cell_polyhedron.Tpo -c -o src/geom/libmesh_devel_la-cell_polyhedron.lo `test -f 'src/geom/cell_polyhedron.C' || echo '$(srcdir)/'`src/geom/cell_polyhedron.C +@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_devel_la-cell_polyhedron.Tpo src/geom/$(DEPDIR)/libmesh_devel_la-cell_polyhedron.Plo +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ $(AM_V_CXX)source='src/geom/cell_polyhedron.C' object='src/geom/libmesh_devel_la-cell_polyhedron.lo' libtool=yes @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_devel_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_devel_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_devel_la-cell_polyhedron.lo `test -f 'src/geom/cell_polyhedron.C' || echo '$(srcdir)/'`src/geom/cell_polyhedron.C + src/geom/libmesh_devel_la-cell_prism.lo: src/geom/cell_prism.C @am__fastdepCXX_TRUE@ $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_devel_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_devel_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_devel_la-cell_prism.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_devel_la-cell_prism.Tpo -c -o src/geom/libmesh_devel_la-cell_prism.lo `test -f 'src/geom/cell_prism.C' || echo '$(srcdir)/'`src/geom/cell_prism.C @am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_devel_la-cell_prism.Tpo src/geom/$(DEPDIR)/libmesh_devel_la-cell_prism.Plo @@ -24709,6 +24799,13 @@ src/geom/libmesh_oprof_la-cell.lo: src/geom/cell.C @AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ @am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_oprof_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_oprof_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_oprof_la-cell.lo `test -f 'src/geom/cell.C' || echo '$(srcdir)/'`src/geom/cell.C +src/geom/libmesh_oprof_la-cell_c0polyhedron.lo: src/geom/cell_c0polyhedron.C +@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_oprof_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_oprof_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_oprof_la-cell_c0polyhedron.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_oprof_la-cell_c0polyhedron.Tpo -c -o src/geom/libmesh_oprof_la-cell_c0polyhedron.lo `test -f 'src/geom/cell_c0polyhedron.C' || echo '$(srcdir)/'`src/geom/cell_c0polyhedron.C +@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_oprof_la-cell_c0polyhedron.Tpo src/geom/$(DEPDIR)/libmesh_oprof_la-cell_c0polyhedron.Plo +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ $(AM_V_CXX)source='src/geom/cell_c0polyhedron.C' object='src/geom/libmesh_oprof_la-cell_c0polyhedron.lo' libtool=yes @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_oprof_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_oprof_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_oprof_la-cell_c0polyhedron.lo `test -f 'src/geom/cell_c0polyhedron.C' || echo '$(srcdir)/'`src/geom/cell_c0polyhedron.C + src/geom/libmesh_oprof_la-cell_hex.lo: src/geom/cell_hex.C @am__fastdepCXX_TRUE@ $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_oprof_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_oprof_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_oprof_la-cell_hex.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_oprof_la-cell_hex.Tpo -c -o src/geom/libmesh_oprof_la-cell_hex.lo `test -f 'src/geom/cell_hex.C' || echo '$(srcdir)/'`src/geom/cell_hex.C @am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_oprof_la-cell_hex.Tpo src/geom/$(DEPDIR)/libmesh_oprof_la-cell_hex.Plo @@ -24793,6 +24890,13 @@ src/geom/libmesh_oprof_la-cell_inf_prism6.lo: src/geom/cell_inf_prism6.C @AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ @am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_oprof_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_oprof_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_oprof_la-cell_inf_prism6.lo `test -f 'src/geom/cell_inf_prism6.C' || echo '$(srcdir)/'`src/geom/cell_inf_prism6.C +src/geom/libmesh_oprof_la-cell_polyhedron.lo: src/geom/cell_polyhedron.C +@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_oprof_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_oprof_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_oprof_la-cell_polyhedron.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_oprof_la-cell_polyhedron.Tpo -c -o src/geom/libmesh_oprof_la-cell_polyhedron.lo `test -f 'src/geom/cell_polyhedron.C' || echo '$(srcdir)/'`src/geom/cell_polyhedron.C +@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_oprof_la-cell_polyhedron.Tpo src/geom/$(DEPDIR)/libmesh_oprof_la-cell_polyhedron.Plo +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ $(AM_V_CXX)source='src/geom/cell_polyhedron.C' object='src/geom/libmesh_oprof_la-cell_polyhedron.lo' libtool=yes @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_oprof_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_oprof_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_oprof_la-cell_polyhedron.lo `test -f 'src/geom/cell_polyhedron.C' || echo '$(srcdir)/'`src/geom/cell_polyhedron.C + src/geom/libmesh_oprof_la-cell_prism.lo: src/geom/cell_prism.C @am__fastdepCXX_TRUE@ $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_oprof_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_oprof_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_oprof_la-cell_prism.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_oprof_la-cell_prism.Tpo -c -o src/geom/libmesh_oprof_la-cell_prism.lo `test -f 'src/geom/cell_prism.C' || echo '$(srcdir)/'`src/geom/cell_prism.C @am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_oprof_la-cell_prism.Tpo src/geom/$(DEPDIR)/libmesh_oprof_la-cell_prism.Plo @@ -27971,6 +28075,13 @@ src/geom/libmesh_opt_la-cell.lo: src/geom/cell.C @AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ @am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_opt_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_opt_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_opt_la-cell.lo `test -f 'src/geom/cell.C' || echo '$(srcdir)/'`src/geom/cell.C +src/geom/libmesh_opt_la-cell_c0polyhedron.lo: src/geom/cell_c0polyhedron.C +@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_opt_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_opt_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_opt_la-cell_c0polyhedron.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_opt_la-cell_c0polyhedron.Tpo -c -o src/geom/libmesh_opt_la-cell_c0polyhedron.lo `test -f 'src/geom/cell_c0polyhedron.C' || echo '$(srcdir)/'`src/geom/cell_c0polyhedron.C +@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_opt_la-cell_c0polyhedron.Tpo src/geom/$(DEPDIR)/libmesh_opt_la-cell_c0polyhedron.Plo +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ $(AM_V_CXX)source='src/geom/cell_c0polyhedron.C' object='src/geom/libmesh_opt_la-cell_c0polyhedron.lo' libtool=yes @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_opt_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_opt_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_opt_la-cell_c0polyhedron.lo `test -f 'src/geom/cell_c0polyhedron.C' || echo '$(srcdir)/'`src/geom/cell_c0polyhedron.C + src/geom/libmesh_opt_la-cell_hex.lo: src/geom/cell_hex.C @am__fastdepCXX_TRUE@ $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_opt_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_opt_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_opt_la-cell_hex.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_opt_la-cell_hex.Tpo -c -o src/geom/libmesh_opt_la-cell_hex.lo `test -f 'src/geom/cell_hex.C' || echo '$(srcdir)/'`src/geom/cell_hex.C @am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_opt_la-cell_hex.Tpo src/geom/$(DEPDIR)/libmesh_opt_la-cell_hex.Plo @@ -28055,6 +28166,13 @@ src/geom/libmesh_opt_la-cell_inf_prism6.lo: src/geom/cell_inf_prism6.C @AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ @am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_opt_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_opt_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_opt_la-cell_inf_prism6.lo `test -f 'src/geom/cell_inf_prism6.C' || echo '$(srcdir)/'`src/geom/cell_inf_prism6.C +src/geom/libmesh_opt_la-cell_polyhedron.lo: src/geom/cell_polyhedron.C +@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_opt_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_opt_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_opt_la-cell_polyhedron.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_opt_la-cell_polyhedron.Tpo -c -o src/geom/libmesh_opt_la-cell_polyhedron.lo `test -f 'src/geom/cell_polyhedron.C' || echo '$(srcdir)/'`src/geom/cell_polyhedron.C +@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_opt_la-cell_polyhedron.Tpo src/geom/$(DEPDIR)/libmesh_opt_la-cell_polyhedron.Plo +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ $(AM_V_CXX)source='src/geom/cell_polyhedron.C' object='src/geom/libmesh_opt_la-cell_polyhedron.lo' libtool=yes @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_opt_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_opt_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_opt_la-cell_polyhedron.lo `test -f 'src/geom/cell_polyhedron.C' || echo '$(srcdir)/'`src/geom/cell_polyhedron.C + src/geom/libmesh_opt_la-cell_prism.lo: src/geom/cell_prism.C @am__fastdepCXX_TRUE@ $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_opt_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_opt_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_opt_la-cell_prism.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_opt_la-cell_prism.Tpo -c -o src/geom/libmesh_opt_la-cell_prism.lo `test -f 'src/geom/cell_prism.C' || echo '$(srcdir)/'`src/geom/cell_prism.C @am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_opt_la-cell_prism.Tpo src/geom/$(DEPDIR)/libmesh_opt_la-cell_prism.Plo @@ -31233,6 +31351,13 @@ src/geom/libmesh_prof_la-cell.lo: src/geom/cell.C @AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ @am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_prof_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_prof_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_prof_la-cell.lo `test -f 'src/geom/cell.C' || echo '$(srcdir)/'`src/geom/cell.C +src/geom/libmesh_prof_la-cell_c0polyhedron.lo: src/geom/cell_c0polyhedron.C +@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_prof_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_prof_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_prof_la-cell_c0polyhedron.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_prof_la-cell_c0polyhedron.Tpo -c -o src/geom/libmesh_prof_la-cell_c0polyhedron.lo `test -f 'src/geom/cell_c0polyhedron.C' || echo '$(srcdir)/'`src/geom/cell_c0polyhedron.C +@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_prof_la-cell_c0polyhedron.Tpo src/geom/$(DEPDIR)/libmesh_prof_la-cell_c0polyhedron.Plo +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ $(AM_V_CXX)source='src/geom/cell_c0polyhedron.C' object='src/geom/libmesh_prof_la-cell_c0polyhedron.lo' libtool=yes @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_prof_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_prof_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_prof_la-cell_c0polyhedron.lo `test -f 'src/geom/cell_c0polyhedron.C' || echo '$(srcdir)/'`src/geom/cell_c0polyhedron.C + src/geom/libmesh_prof_la-cell_hex.lo: src/geom/cell_hex.C @am__fastdepCXX_TRUE@ $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_prof_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_prof_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_prof_la-cell_hex.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_prof_la-cell_hex.Tpo -c -o src/geom/libmesh_prof_la-cell_hex.lo `test -f 'src/geom/cell_hex.C' || echo '$(srcdir)/'`src/geom/cell_hex.C @am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_prof_la-cell_hex.Tpo src/geom/$(DEPDIR)/libmesh_prof_la-cell_hex.Plo @@ -31317,6 +31442,13 @@ src/geom/libmesh_prof_la-cell_inf_prism6.lo: src/geom/cell_inf_prism6.C @AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ @am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_prof_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_prof_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_prof_la-cell_inf_prism6.lo `test -f 'src/geom/cell_inf_prism6.C' || echo '$(srcdir)/'`src/geom/cell_inf_prism6.C +src/geom/libmesh_prof_la-cell_polyhedron.lo: src/geom/cell_polyhedron.C +@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_prof_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_prof_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_prof_la-cell_polyhedron.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_prof_la-cell_polyhedron.Tpo -c -o src/geom/libmesh_prof_la-cell_polyhedron.lo `test -f 'src/geom/cell_polyhedron.C' || echo '$(srcdir)/'`src/geom/cell_polyhedron.C +@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_prof_la-cell_polyhedron.Tpo src/geom/$(DEPDIR)/libmesh_prof_la-cell_polyhedron.Plo +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ $(AM_V_CXX)source='src/geom/cell_polyhedron.C' object='src/geom/libmesh_prof_la-cell_polyhedron.lo' libtool=yes @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_prof_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_prof_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_prof_la-cell_polyhedron.lo `test -f 'src/geom/cell_polyhedron.C' || echo '$(srcdir)/'`src/geom/cell_polyhedron.C + src/geom/libmesh_prof_la-cell_prism.lo: src/geom/cell_prism.C @am__fastdepCXX_TRUE@ $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_prof_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_prof_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_prof_la-cell_prism.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_prof_la-cell_prism.Tpo -c -o src/geom/libmesh_prof_la-cell_prism.lo `test -f 'src/geom/cell_prism.C' || echo '$(srcdir)/'`src/geom/cell_prism.C @am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_prof_la-cell_prism.Tpo src/geom/$(DEPDIR)/libmesh_prof_la-cell_prism.Plo @@ -35577,6 +35709,7 @@ distclean: distclean-recursive -rm -f src/fe/$(DEPDIR)/libmesh_prof_la-inf_fe_static.Plo -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-bounding_box.Plo -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell.Plo + -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_c0polyhedron.Plo -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_hex.Plo -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_hex20.Plo -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_hex27.Plo @@ -35589,6 +35722,7 @@ distclean: distclean-recursive -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_inf_prism.Plo -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_inf_prism12.Plo -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_inf_prism6.Plo + -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_polyhedron.Plo -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_prism.Plo -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_prism15.Plo -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_prism18.Plo @@ -35643,6 +35777,7 @@ distclean: distclean-recursive -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-surface.Plo -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-bounding_box.Plo -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell.Plo + -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_c0polyhedron.Plo -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_hex.Plo -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_hex20.Plo -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_hex27.Plo @@ -35655,6 +35790,7 @@ distclean: distclean-recursive -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_inf_prism.Plo -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_inf_prism12.Plo -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_inf_prism6.Plo + -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_polyhedron.Plo -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_prism.Plo -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_prism15.Plo -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_prism18.Plo @@ -35709,6 +35845,7 @@ distclean: distclean-recursive -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-surface.Plo -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-bounding_box.Plo -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell.Plo + -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_c0polyhedron.Plo -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_hex.Plo -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_hex20.Plo -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_hex27.Plo @@ -35721,6 +35858,7 @@ distclean: distclean-recursive -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_inf_prism.Plo -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_inf_prism12.Plo -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_inf_prism6.Plo + -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_polyhedron.Plo -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_prism.Plo -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_prism15.Plo -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_prism18.Plo @@ -35775,6 +35913,7 @@ distclean: distclean-recursive -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-surface.Plo -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-bounding_box.Plo -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell.Plo + -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_c0polyhedron.Plo -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_hex.Plo -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_hex20.Plo -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_hex27.Plo @@ -35787,6 +35926,7 @@ distclean: distclean-recursive -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_inf_prism.Plo -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_inf_prism12.Plo -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_inf_prism6.Plo + -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_polyhedron.Plo -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_prism.Plo -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_prism15.Plo -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_prism18.Plo @@ -35841,6 +35981,7 @@ distclean: distclean-recursive -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-surface.Plo -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-bounding_box.Plo -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell.Plo + -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_c0polyhedron.Plo -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_hex.Plo -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_hex20.Plo -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_hex27.Plo @@ -35853,6 +35994,7 @@ distclean: distclean-recursive -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_inf_prism.Plo -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_inf_prism12.Plo -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_inf_prism6.Plo + -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_polyhedron.Plo -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_prism.Plo -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_prism15.Plo -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_prism18.Plo @@ -38015,6 +38157,7 @@ maintainer-clean: maintainer-clean-recursive -rm -f src/fe/$(DEPDIR)/libmesh_prof_la-inf_fe_static.Plo -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-bounding_box.Plo -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell.Plo + -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_c0polyhedron.Plo -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_hex.Plo -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_hex20.Plo -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_hex27.Plo @@ -38027,6 +38170,7 @@ maintainer-clean: maintainer-clean-recursive -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_inf_prism.Plo -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_inf_prism12.Plo -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_inf_prism6.Plo + -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_polyhedron.Plo -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_prism.Plo -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_prism15.Plo -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_prism18.Plo @@ -38081,6 +38225,7 @@ maintainer-clean: maintainer-clean-recursive -rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-surface.Plo -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-bounding_box.Plo -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell.Plo + -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_c0polyhedron.Plo -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_hex.Plo -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_hex20.Plo -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_hex27.Plo @@ -38093,6 +38238,7 @@ maintainer-clean: maintainer-clean-recursive -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_inf_prism.Plo -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_inf_prism12.Plo -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_inf_prism6.Plo + -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_polyhedron.Plo -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_prism.Plo -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_prism15.Plo -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_prism18.Plo @@ -38147,6 +38293,7 @@ maintainer-clean: maintainer-clean-recursive -rm -f src/geom/$(DEPDIR)/libmesh_devel_la-surface.Plo -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-bounding_box.Plo -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell.Plo + -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_c0polyhedron.Plo -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_hex.Plo -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_hex20.Plo -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_hex27.Plo @@ -38159,6 +38306,7 @@ maintainer-clean: maintainer-clean-recursive -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_inf_prism.Plo -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_inf_prism12.Plo -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_inf_prism6.Plo + -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_polyhedron.Plo -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_prism.Plo -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_prism15.Plo -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_prism18.Plo @@ -38213,6 +38361,7 @@ maintainer-clean: maintainer-clean-recursive -rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-surface.Plo -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-bounding_box.Plo -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell.Plo + -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_c0polyhedron.Plo -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_hex.Plo -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_hex20.Plo -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_hex27.Plo @@ -38225,6 +38374,7 @@ maintainer-clean: maintainer-clean-recursive -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_inf_prism.Plo -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_inf_prism12.Plo -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_inf_prism6.Plo + -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_polyhedron.Plo -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_prism.Plo -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_prism15.Plo -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_prism18.Plo @@ -38279,6 +38429,7 @@ maintainer-clean: maintainer-clean-recursive -rm -f src/geom/$(DEPDIR)/libmesh_opt_la-surface.Plo -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-bounding_box.Plo -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell.Plo + -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_c0polyhedron.Plo -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_hex.Plo -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_hex20.Plo -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_hex27.Plo @@ -38291,6 +38442,7 @@ maintainer-clean: maintainer-clean-recursive -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_inf_prism.Plo -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_inf_prism12.Plo -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_inf_prism6.Plo + -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_polyhedron.Plo -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_prism.Plo -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_prism15.Plo -rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_prism18.Plo diff --git a/include/Makefile.in b/include/Makefile.in index dedc6c46d93..9f1f2ce7d51 100644 --- a/include/Makefile.in +++ b/include/Makefile.in @@ -703,6 +703,7 @@ include_HEADERS = \ fe/inf_fe_map.h \ geom/bounding_box.h \ geom/cell.h \ + geom/cell_c0polyhedron.h \ geom/cell_hex.h \ geom/cell_hex20.h \ geom/cell_hex27.h \ @@ -715,6 +716,7 @@ include_HEADERS = \ geom/cell_inf_prism.h \ geom/cell_inf_prism12.h \ geom/cell_inf_prism6.h \ + geom/cell_polyhedron.h \ geom/cell_prism.h \ geom/cell_prism15.h \ geom/cell_prism18.h \ diff --git a/include/enums/enum_elem_type.h b/include/enums/enum_elem_type.h index e0ffc6855c9..fa755f25ca7 100644 --- a/include/enums/enum_elem_type.h +++ b/include/enums/enum_elem_type.h @@ -81,6 +81,8 @@ enum ElemType : int { QUADSHELL9 = 38, // An arbitrary polygon with N EDGE2 sides C0POLYGON = 39, + // An arbitrary polyhedron with C0POLYGON sides + C0POLYHEDRON = 40, // Invalid INVALID_ELEM}; // should always be last diff --git a/include/fe/fe.h b/include/fe/fe.h index 30d82169be0..82f0feb0884 100644 --- a/include/fe/fe.h +++ b/include/fe/fe.h @@ -1533,6 +1533,19 @@ OutputShape fe_fdm_deriv(const ElemType type, (const ElemType, const Order, const unsigned int, const Point &)); +template +OutputShape fe_fdm_deriv(const ElemType type, + const Order order, + const Elem * elem, + const unsigned int i, + const unsigned int j, + const Point & p, + OutputShape(*shape_func) + (const ElemType type, const Order, + const Elem *, const unsigned int, + const Point &)); + + template OutputShape fe_fdm_second_deriv(const Elem * elem, @@ -1558,6 +1571,20 @@ OutputShape fe_fdm_second_deriv(const ElemType type, const unsigned int, const Point &)); +template +OutputShape fe_fdm_second_deriv(const ElemType type, + const Order order, + const Elem * elem, + const unsigned int i, + const unsigned int j, + const Point & p, + OutputShape(*deriv_func) + (const ElemType, const Order, + const Elem *, + const unsigned int, + const unsigned int, + const Point &)); + /** * Helper functions for Lagrange-based basis functions. */ diff --git a/include/geom/cell_c0polyhedron.h b/include/geom/cell_c0polyhedron.h new file mode 100644 index 00000000000..e80018da1c2 --- /dev/null +++ b/include/geom/cell_c0polyhedron.h @@ -0,0 +1,201 @@ +// The libMesh Finite Element Library. +// Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner + +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. + +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. + +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + + +#ifndef LIBMESH_FACE_C0POLYHEDRON_H +#define LIBMESH_FACE_C0POLYHEDRON_H + +// Local includes +#include "libmesh/libmesh_common.h" +#include "libmesh/cell_polyhedron.h" +#include "libmesh/enum_order.h" + +namespace libMesh +{ + +/** + * The \p C0Polyhedron is an element in 3D with an arbitrary (but fixed) + * number of polygonal first-order (C0Polygon) sides. + * + * In master space ... I give up. We've practically got to have a + * definition distinct for every physical element anyway, and we've + * got to give FE and QBase objects access to the physical element, so + * we'll just do all our FE and Quadrature calculations there. + * + * \author Roy H. Stogner + * \date 2025 + * \brief A 3D element with an arbitrary number of polygonal + * first-order sides. + */ +class C0Polyhedron : public Polyhedron +{ +public: + + /** + * Constructor. Takes the sides as an input, and allocates nodes + * accordingly. + * + * By default this element has no parent. + * + * We'll set up a simple default tetrahedralization here, but if users + * want to adjust node positions later they'll want to retriangulate() + * after. + */ + explicit + C0Polyhedron (const std::vector> & sides, + Elem * p=nullptr); + + C0Polyhedron (C0Polyhedron &&) = delete; + C0Polyhedron (const C0Polyhedron &) = delete; + C0Polyhedron & operator= (const C0Polyhedron &) = delete; + C0Polyhedron & operator= (C0Polyhedron &&) = delete; + virtual ~C0Polyhedron(); + + /** + * \returns An Elem of the specified type, which should be the same + * type as \p this, wrapped in a smart pointer. + * + * This is not a complete clone() method (since e.g. it does not set + * node pointers; the standard use case reassigns node pointers from + * a different mesh), but it allows us to more easily copy polyhedra + * objects that depend on side objects and require existing nodes + * and so cannot be constructed with Elem::build(). + */ + std::unique_ptr disconnected_clone() const override; + + /** + * \returns \p C0POLYHEDRON. + */ + virtual ElemType type () const override final { return C0POLYHEDRON; } + + /** + * \returns the number of vertices. For the simple C0Polyhedron + * every node is a vertex. + */ + virtual unsigned int n_vertices() const override final { return this->_nodelinks_data.size(); } + + /** + * \returns the number of tetrahedra to break this into for + * visualization. + */ + virtual unsigned int n_sub_elem() const override + { return this->n_subelements(); } + + /** + * \returns \p true if the specified (local) node number is a vertex. + */ + virtual bool is_vertex(const unsigned int i) const override; + + /** + * \returns \p true if the specified (local) node number is an edge. + */ + virtual bool is_edge(const unsigned int i) const override; + + /** + * \returns \p true if the specified (local) node number is a face. + */ + virtual bool is_face(const unsigned int i) const override; + + /** + * \returns \p true if the specified (local) node number is on the + * specified side. + */ + virtual bool is_node_on_side(const unsigned int n, + const unsigned int s) const override; + + virtual std::vector nodes_on_side(const unsigned int s) const override; + + virtual std::vector nodes_on_edge(const unsigned int e) const override; + + /** + * \returns \p true if the specified (local) node number is on the + * specified edge. + */ + virtual bool is_node_on_edge(const unsigned int n, + const unsigned int e) const override; + + /** + * \returns \p true if the element map is definitely affine within + * numerical tolerances. We don't even share a master element from + * element to element, so we're going to return false here. + */ + virtual bool has_affine_map () const override { return false; } + + /** + * \returns FIRST + */ + virtual Order default_order() const override { return FIRST; } + + /** + * Don't hide Elem::key() defined in the base class. + */ + using Elem::key; + + virtual void connectivity(const unsigned int sf, + const IOPackage iop, + std::vector & conn) const override; + + /** + * Element refinement is not implemented for polyhedra. + */ + virtual std::pair + second_order_child_vertex (const unsigned int n) const override; + + /** + * An optimized method for calculating the area of a C0Polyhedron. + */ + virtual Real volume () const override; + + /** + * An optimized method for calculating the centroid of a C0Polyhedron. + */ + virtual Point true_centroid () const override; + + ElemType side_type (const unsigned int s) const override final; + + /** + * Create a triangulation (tetrahedralization) based on the current + * sides' triangulations. + */ + virtual void retriangulate() override final; + +protected: + + /** + * Add to our triangulation (tetrahedralization). + */ + void add_tet(int n1, int n2, int n3, int n4); + +#ifdef LIBMESH_ENABLE_AMR + + /** + * Matrix used to create the elements children. + */ + virtual Real embedding_matrix (const unsigned int /*i*/, + const unsigned int /*j*/, + const unsigned int /*k*/) const override + { libmesh_not_implemented(); return 0; } + +#endif // LIBMESH_ENABLE_AMR + +}; + + +} // namespace libMesh + +#endif // LIBMESH_FACE_C0POLYHEDRON_H diff --git a/include/geom/cell_polyhedron.h b/include/geom/cell_polyhedron.h new file mode 100644 index 00000000000..992cd1cbd6e --- /dev/null +++ b/include/geom/cell_polyhedron.h @@ -0,0 +1,380 @@ +// The libMesh Finite Element Library. +// Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner + +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. + +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. + +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + + +#ifndef LIBMESH_CELL_POLYHEDRON_H +#define LIBMESH_CELL_POLYHEDRON_H + + +// Local includes +#include "libmesh/libmesh_common.h" +#include "libmesh/cell.h" + +namespace libMesh +{ + +// Forward declarations +class Polygon; + +/** + * The \p Polyhedron is an element in 3D with an arbitrary + * number of polygonal faces. + * + * \author Roy H. Stogner + * \date 2025 + * \brief A 3D element with an arbitrary number of polygonal faces + */ +class Polyhedron : public Cell +{ +public: + + /** + * Arbitrary polyhedral element, takes a vector of shared pointers + * to sides (which should already be constructed), and a (probably + * unused) parent link. Derived classes implement 'true' elements. + */ + Polyhedron (const std::vector> & sides, + Elem * p); + + Polyhedron (Polyhedron &&) = delete; + Polyhedron (const Polyhedron &) = delete; + Polyhedron & operator= (const Polyhedron &) = delete; + Polyhedron & operator= (Polyhedron &&) = delete; + virtual ~Polyhedron() = default; + + /** + * \returns true - polyhedron subclasses can have numbers of sides + * and nodes which vary at runtime. + */ + virtual bool runtime_topology() const override { return true; } + + /** + * \returns The \p Point associated with local \p Node \p i, + * in master element rather than physical coordinates. + * + * This implementation returns the master vertices; subclasses with + * more points will need to further override. + * + * Trying to come up with a "master" shape for an arbitrary + * polyhedron is *too* arbitrary, so we're just returning physical + * points here. + */ + virtual Point master_point (const unsigned int i) const override; + + static const int num_children = 0; // Refinement not yet supported + + /** + * \returns the number of nodes in the polyhedron. + */ + virtual unsigned int n_nodes() const override { return this->_nodelinks_data.size(); } + + /** + * \returns the number of sides, which is the number of links we had + * to save minus one parent link and minus one interior_parent link. + */ + virtual unsigned int n_sides() const override final { return _elemlinks_data.size()-2; } + + /** + * \returns the number of edges. + */ + virtual unsigned int n_edges() const override final { return _edge_lookup.size(); } + + /** + * \returns the number of faces, which is the number of links we had + * to save minus one parent link and minus one interior_parent link. + */ + virtual unsigned int n_faces() const override final { return _elemlinks_data.size()-2; } + + /** + * \returns num_children. + */ + virtual unsigned int n_children() const override final { return num_children; } + + /** + * \returns \p true if the specified child is on the + * specified side. + * + * Not implemented ... indefinitely? I don't think we'll be doing + * hierarchic refinement for general polyhedra. + */ + virtual bool is_child_on_side(const unsigned int c, + const unsigned int s) const override; + + /** + * Throws an error. I'm not sure how to define the "opposite + * side" of a polyhedron. + */ + virtual unsigned int opposite_side(const unsigned int s) const override final; + + /** + * Throws an error - \p opposite_side(s) is too hard to define in + * general on polyhedra. + */ + virtual unsigned int opposite_node(const unsigned int n, + const unsigned int s) const override final; + + /** + * Don't hide Elem::key() defined in the base class. + */ + using Elem::key; + + /** + * \returns An id associated with the global node ids of this + * element. The id is not necessarily unique, but should be + * close. + */ + virtual dof_id_type key () const override; + + /** + * \returns An id associated with the \p s side of this element. + * The id is not necessarily unique, but should be close. + */ + virtual dof_id_type key (const unsigned int s) const override; + + /** + * \returns An id associated with the \p s side of this element, as + * defined solely by element vertices. The id is not necessarily + * unique, but should be close. This is particularly useful in the + * \p MeshBase::find_neighbors() routine. + */ + virtual dof_id_type low_order_key (const unsigned int s) const override; + + /** + * \returns The local node id for node \p side_node on side \p side of + * this Elem. + */ + virtual unsigned int local_side_node(unsigned int side, + unsigned int side_node) const override; + + /** + * Similar to Elem::local_side_node(), but instead of a side id, takes + * an edge id and a node id on that edge and returns a local node number + * for the Elem. + */ + virtual unsigned int local_edge_node(unsigned int edge, + unsigned int edge_node) const override; + + /** + * \returns A simple Polygon face for side i. + */ + virtual std::unique_ptr side_ptr (const unsigned int i) override final; + + /** + * Rebuilds a simple Polygon coincident with face i. + */ + virtual void side_ptr (std::unique_ptr & elem, + const unsigned int i) override final; + + /** + * Copies the Polygon side coincident with side i. + */ + virtual std::unique_ptr build_side_ptr (const unsigned int i, + bool proxy=false) override; + + /** + * Copies the Polygon side coincident with side i. + */ + virtual void build_side_ptr (std::unique_ptr & elem, + const unsigned int i) override; + + /** + * \returns An element coincident with edge \p i wrapped in a smart pointer. + */ + virtual std::unique_ptr build_edge_ptr (const unsigned int i) override final; + + /** + * Resets the loose element \p edge, which may currently point to a + * different edge than \p i or even a different element than \p + * this, to point to edge \p i on \p this. + */ + virtual void build_edge_ptr (std::unique_ptr & edge, const unsigned int i) override final; + + /** + * \returns The suggested quality bounds for + * the polyhedron based on quality measure q. + */ + virtual std::pair qual_bounds (const ElemQuality q) const override; + + /** + * \returns the (local) side numbers that touch the specified edge. + */ + virtual std::vector sides_on_edge(const unsigned int e) const override final; + + /** + * \returns \p true if the specified edge is on the specified side. + */ + virtual bool is_edge_on_side(const unsigned int e, + const unsigned int s) const override final; + + /** + * Maybe we have non-identity permutations, but trying to figure out + * how many is an exercise in applied group theory, which is a bit + * much for just expanded unit test coverage. + */ + virtual unsigned int n_permutations() const override { return 1; } + + virtual void permute(unsigned int libmesh_dbg_var(perm_num)) override final + { libmesh_assert_equal_to(perm_num, 0); } + + virtual bool is_flipped() const override final; + + /** + * A flip is one of those general non-identity permutations we can't + * handle. + * + * But we'll call it "not implemented" just in cases someone someday + * wants to write an implementation to handle polyhedra with + * topological symmetry planes. + */ + virtual void flip(BoundaryInfo *) override final { libmesh_not_implemented(); }; + + virtual std::vector edges_adjacent_to_node(const unsigned int n) const override; + + /** + * Create a triangulation (tetrahedralization) from the current node + * locations and face triangulations. + * + * If this is not called by the user, it will be called + * automatically when needed. + * + * If the user moves the polyhedron nodes, the triangulation may + * need to be regenerated manually, or the changed mapping may be + * lower quality or even inverted. + * + * Pure virtual because this strongly depends on what mid-edge or + * mid-face nodes the subclass might have. + */ + virtual void retriangulate() = 0; + + /** + * \returns true iff the polyhedron is convex. Some Polyhedron + * methods currently assume (or in debugging modes, test for and + * throw errors if not finding) convexity. + */ + bool convex(); + + /** + * \returns the size of the triangulation of the polyhedron + */ + unsigned int n_subelements() const { return cast_int(this->_triangulation.size()); } + + /** + * \returns the local indices of points on a subelement of the + * polyhedron + * + * Each subelement here is a tetrahedron + */ + virtual std::array subelement (unsigned int i) const + { + libmesh_assert_less(i, this->_triangulation.size()); + return this->_triangulation[i]; + } + + /** + * \returns the master-space points of a subelement of the + * polyhedron + */ + virtual std::array master_subelement (unsigned int i) const; + + /** + * \returns the index of a subelement containing the master-space + * point \p p, along with (the first three) barycentric coordinates + * for \p; or return invalid_uint if no subelement contains p to + * within tolerance \p tol. + */ + std::tuple + subelement_coordinates (const Point & p, + Real tol = TOLERANCE*TOLERANCE) const; + + virtual bool on_reference_element(const Point & p, + const Real eps = TOLERANCE) const override final; + +protected: + + /** + * \returns Clones of the sides of \p this, wrapped in smart + * pointers. + * + * This factors out much of the code required by the + * disconnected_clone() method of Polyhedron subclasses. + */ + std::vector> side_clones() const; + + /** + * Helper method for finding the non-cached side that shares an + * edge, by examining the local node ids there. + */ + bool side_has_edge_nodes(unsigned int side, + unsigned int min_node, + unsigned int max_node) const; + + /** + * Data for links to parent/neighbor/interior_parent elements. + * + * There should be num_sides neighbors, preceded by any parent link + * and followed by any interior_parent link. + * + * We're stuck with a heap allocation here since the number of sides + * in a subclass is chosen at runtime. + */ + std::vector _elemlinks_data; + + /** + * Data for links to nodes. We're stuck with a heap allocation here + * since the number of nodes in a subclass is chosen at runtime. + */ + std::vector _nodelinks_data; + + /** + * Data for links to sides. We use shared_ptr so we can keep sides + * consistent between neighbors. We also have a bool for each + * polygon to indicate whether the zeta (xi cross eta) direction for + * the polygon points into this polyhedron (true) or out of it + * (false), and a map from side-local node number to element-local + * node number. + */ + std::vector, + bool, + std::vector>> _sidelinks_data; + + /** + * One entry for each polyhedron edge, a pair indicating the side + * number and the edge-of-side number which identifies the edge. + * + * Although each edge can be reached from two sides, only the + * lower-numbered side is in this lookup table. For edges accessed + * via the same side, the side edge numbering matches our edge + * numbering. + */ + std::vector> _edge_lookup; + + /** + * Data for a triangulation (tetrahedralization) of the polyhedron. + * + * Positive int values should correspond to local node numbers. + * + * Negative int values may correspond to "special" points, e.g. a + * centroid or skeleton point. + */ + std::vector> _triangulation; +}; + + +} // namespace libMesh + +#endif // LIBMESH_CELL_POLYHEDRON_H diff --git a/include/geom/face.h b/include/geom/face.h index 7a31e07c449..87fa0741d92 100644 --- a/include/geom/face.h +++ b/include/geom/face.h @@ -78,6 +78,11 @@ class Face : public Elem virtual void build_edge_ptr (std::unique_ptr & edge, const unsigned int i) override final { build_side_ptr(edge, i); } + /** + * For the const build_edge_ptr we use the parent method + */ + using Elem::build_edge_ptr; + /** * is_edge_on_side is trivial in 2D. */ diff --git a/include/include_HEADERS b/include/include_HEADERS index 492c84c5fb1..5a31513b7ee 100644 --- a/include/include_HEADERS +++ b/include/include_HEADERS @@ -102,6 +102,7 @@ include_HEADERS = \ fe/inf_fe_map.h \ geom/bounding_box.h \ geom/cell.h \ + geom/cell_c0polyhedron.h \ geom/cell_hex.h \ geom/cell_hex20.h \ geom/cell_hex27.h \ @@ -114,6 +115,7 @@ include_HEADERS = \ geom/cell_inf_prism.h \ geom/cell_inf_prism12.h \ geom/cell_inf_prism6.h \ + geom/cell_polyhedron.h \ geom/cell_prism.h \ geom/cell_prism15.h \ geom/cell_prism18.h \ diff --git a/include/libmesh/Makefile.am b/include/libmesh/Makefile.am index acd27672ca0..fca975b85e8 100644 --- a/include/libmesh/Makefile.am +++ b/include/libmesh/Makefile.am @@ -92,6 +92,7 @@ BUILT_SOURCES = \ inf_fe_map.h \ bounding_box.h \ cell.h \ + cell_c0polyhedron.h \ cell_hex.h \ cell_hex20.h \ cell_hex27.h \ @@ -104,6 +105,7 @@ BUILT_SOURCES = \ cell_inf_prism.h \ cell_inf_prism12.h \ cell_inf_prism6.h \ + cell_polyhedron.h \ cell_prism.h \ cell_prism15.h \ cell_prism18.h \ @@ -869,6 +871,9 @@ bounding_box.h: $(top_srcdir)/include/geom/bounding_box.h cell.h: $(top_srcdir)/include/geom/cell.h $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ +cell_c0polyhedron.h: $(top_srcdir)/include/geom/cell_c0polyhedron.h + $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ + cell_hex.h: $(top_srcdir)/include/geom/cell_hex.h $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ @@ -905,6 +910,9 @@ cell_inf_prism12.h: $(top_srcdir)/include/geom/cell_inf_prism12.h cell_inf_prism6.h: $(top_srcdir)/include/geom/cell_inf_prism6.h $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ +cell_polyhedron.h: $(top_srcdir)/include/geom/cell_polyhedron.h + $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ + cell_prism.h: $(top_srcdir)/include/geom/cell_prism.h $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ diff --git a/include/libmesh/Makefile.in b/include/libmesh/Makefile.in index 6d632031521..5e5e169cb03 100644 --- a/include/libmesh/Makefile.in +++ b/include/libmesh/Makefile.in @@ -552,26 +552,27 @@ BUILT_SOURCES = auto_ptr.h dirichlet_boundaries.h dof_map.h \ h1_fe_transformation.h hcurl_fe_transformation.h \ hdiv_fe_transformation.h inf_fe.h inf_fe_instantiate_1D.h \ inf_fe_instantiate_2D.h inf_fe_instantiate_3D.h inf_fe_macro.h \ - inf_fe_map.h bounding_box.h cell.h cell_hex.h cell_hex20.h \ - cell_hex27.h cell_hex8.h cell_inf.h cell_inf_hex.h \ - cell_inf_hex16.h cell_inf_hex18.h cell_inf_hex8.h \ - cell_inf_prism.h cell_inf_prism12.h cell_inf_prism6.h \ - cell_prism.h cell_prism15.h cell_prism18.h cell_prism20.h \ - cell_prism21.h cell_prism6.h cell_pyramid.h cell_pyramid13.h \ - cell_pyramid14.h cell_pyramid18.h cell_pyramid5.h cell_tet.h \ - cell_tet10.h cell_tet14.h cell_tet4.h compare_elems_by_level.h \ - edge.h edge_edge2.h edge_edge3.h edge_edge4.h edge_inf_edge2.h \ - elem.h elem_cutter.h elem_hash.h elem_internal.h \ - elem_quality.h elem_range.h elem_side_builder.h face.h \ - face_c0polygon.h face_inf_quad.h face_inf_quad4.h \ - face_inf_quad6.h face_polygon.h face_quad.h face_quad4.h \ - face_quad4_shell.h face_quad8.h face_quad8_shell.h \ - face_quad9.h face_quad9_shell.h face_tri.h face_tri3.h \ - face_tri3_shell.h face_tri3_subdivision.h face_tri6.h \ - face_tri7.h node.h node_elem.h node_range.h plane.h point.h \ - reference_elem.h remote_elem.h side.h sphere.h stored_range.h \ - surface.h default_coupling.h ghost_point_neighbors.h \ - ghosting_functor.h non_manifold_coupling.h overlap_coupling.h \ + inf_fe_map.h bounding_box.h cell.h cell_c0polyhedron.h \ + cell_hex.h cell_hex20.h cell_hex27.h cell_hex8.h cell_inf.h \ + cell_inf_hex.h cell_inf_hex16.h cell_inf_hex18.h \ + cell_inf_hex8.h cell_inf_prism.h cell_inf_prism12.h \ + cell_inf_prism6.h cell_polyhedron.h cell_prism.h \ + cell_prism15.h cell_prism18.h cell_prism20.h cell_prism21.h \ + cell_prism6.h cell_pyramid.h cell_pyramid13.h cell_pyramid14.h \ + cell_pyramid18.h cell_pyramid5.h cell_tet.h cell_tet10.h \ + cell_tet14.h cell_tet4.h compare_elems_by_level.h edge.h \ + edge_edge2.h edge_edge3.h edge_edge4.h edge_inf_edge2.h elem.h \ + elem_cutter.h elem_hash.h elem_internal.h elem_quality.h \ + elem_range.h elem_side_builder.h face.h face_c0polygon.h \ + face_inf_quad.h face_inf_quad4.h face_inf_quad6.h \ + face_polygon.h face_quad.h face_quad4.h face_quad4_shell.h \ + face_quad8.h face_quad8_shell.h face_quad9.h \ + face_quad9_shell.h face_tri.h face_tri3.h face_tri3_shell.h \ + face_tri3_subdivision.h face_tri6.h face_tri7.h node.h \ + node_elem.h node_range.h plane.h point.h reference_elem.h \ + remote_elem.h side.h sphere.h stored_range.h surface.h \ + default_coupling.h ghost_point_neighbors.h ghosting_functor.h \ + non_manifold_coupling.h overlap_coupling.h \ point_neighbor_coupling.h sibling_coupling.h abaqus_io.h \ boundary_info.h boundary_mesh.h checkpoint_io.h \ distributed_mesh.h dyna_io.h ensight_io.h exodusII_io.h \ @@ -1200,6 +1201,9 @@ bounding_box.h: $(top_srcdir)/include/geom/bounding_box.h cell.h: $(top_srcdir)/include/geom/cell.h $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ +cell_c0polyhedron.h: $(top_srcdir)/include/geom/cell_c0polyhedron.h + $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ + cell_hex.h: $(top_srcdir)/include/geom/cell_hex.h $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ @@ -1236,6 +1240,9 @@ cell_inf_prism12.h: $(top_srcdir)/include/geom/cell_inf_prism12.h cell_inf_prism6.h: $(top_srcdir)/include/geom/cell_inf_prism6.h $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ +cell_polyhedron.h: $(top_srcdir)/include/geom/cell_polyhedron.h + $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ + cell_prism.h: $(top_srcdir)/include/geom/cell_prism.h $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ diff --git a/src/fe/fe.C b/src/fe/fe.C index 5ca9c884677..333b985ff97 100644 --- a/src/fe/fe.C +++ b/src/fe/fe.C @@ -745,70 +745,12 @@ void FE::init_base_shape_functions(const std::vector & qp, #endif // LIBMESH_ENABLE_INFINITE_ELEMENTS -template -OutputShape fe_fdm_deriv(const Elem * elem, - const Order order, - const unsigned int i, - const unsigned int j, - const Point & p, - const bool add_p_level, - OutputShape(*shape_func) - (const Elem *, const Order, - const unsigned int, const Point &, - const bool)) -{ - libmesh_assert(elem); - - libmesh_assert_less (j, LIBMESH_DIM); - - // cheat by using finite difference approximations: - const Real eps = TOLERANCE; - Point pp = p, pm = p; - - switch (j) - { - // d()/dxi - case 0: - { - pp(0) += eps; - pm(0) -= eps; - break; - } - - // d()/deta - case 1: - { - pp(1) += eps; - pm(1) -= eps; - break; - } - - // d()/dzeta - case 2: - { - pp(2) += eps; - pm(2) -= eps; - break; - } - - default: - libmesh_error_msg("Invalid derivative index j = " << j); - } - - return (shape_func(elem, order, i, pp, add_p_level) - - shape_func(elem, order, i, pm, add_p_level))/2./eps; -} - +// Helper for FDM methods +namespace { +using namespace libMesh; -template -OutputShape fe_fdm_deriv(const ElemType type, - const Order order, - const unsigned int i, - const unsigned int j, - const Point & p, - OutputShape(*shape_func) - (const ElemType, const Order, - const unsigned int, const Point &)) +std::tuple +fdm_points(const unsigned int j, const Point & p) { libmesh_assert_less (j, LIBMESH_DIM); @@ -846,23 +788,11 @@ OutputShape fe_fdm_deriv(const ElemType type, libmesh_error_msg("Invalid derivative index j = " << j); } - return (shape_func(type, order, i, pp) - - shape_func(type, order, i, pm))/2./eps; + return std::make_tuple(pp, pm, eps); } - -template -OutputShape -fe_fdm_second_deriv(const Elem * elem, - const Order order, - const unsigned int i, - const unsigned int j, - const Point & p, - const bool add_p_level, - OutputShape(*deriv_func) - (const Elem *, const Order, - const unsigned int, const unsigned int, - const Point &, const bool)) +std::tuple +fdm_second_points(const unsigned int j, const Point & p) { // cheat by using finite difference approximations: const Real eps = 1.e-5; @@ -929,6 +859,84 @@ fe_fdm_second_deriv(const Elem * elem, libmesh_error_msg("Invalid shape function derivative j = " << j); } + return std::make_tuple(pp, pm, eps, deriv_j); +} + +} + + +template +OutputShape fe_fdm_deriv(const Elem * elem, + const Order order, + const unsigned int i, + const unsigned int j, + const Point & p, + const bool add_p_level, + OutputShape(*shape_func) + (const Elem *, const Order, + const unsigned int, const Point &, + const bool)) +{ + libmesh_assert(elem); + + auto [pp, pm, eps] = fdm_points(j, p); + + return (shape_func(elem, order, i, pp, add_p_level) - + shape_func(elem, order, i, pm, add_p_level))/2./eps; +} + + +template +OutputShape fe_fdm_deriv(const ElemType type, + const Order order, + const unsigned int i, + const unsigned int j, + const Point & p, + OutputShape(*shape_func) + (const ElemType, const Order, + const unsigned int, const Point &)) +{ + auto [pp, pm, eps] = fdm_points(j, p); + + return (shape_func(type, order, i, pp) - + shape_func(type, order, i, pm))/2./eps; +} + + +template +OutputShape fe_fdm_deriv(const ElemType type, + const Order order, + const Elem * elem, + const unsigned int i, + const unsigned int j, + const Point & p, + OutputShape(*shape_func) + (const ElemType, const Order, + const Elem *, + const unsigned int, const Point &)) +{ + auto [pp, pm, eps] = fdm_points(j, p); + + return (shape_func(type, order, elem, i, pp) - + shape_func(type, order, elem, i, pm))/2./eps; +} + + +template +OutputShape +fe_fdm_second_deriv(const Elem * elem, + const Order order, + const unsigned int i, + const unsigned int j, + const Point & p, + const bool add_p_level, + OutputShape(*deriv_func) + (const Elem *, const Order, + const unsigned int, const unsigned int, + const Point &, const bool)) +{ + auto [pp, pm, eps, deriv_j] = fdm_second_points(j, p); + return (deriv_func(elem, order, i, deriv_j, pp, add_p_level) - deriv_func(elem, order, i, deriv_j, pm, add_p_level))/2./eps; } @@ -946,73 +954,31 @@ fe_fdm_second_deriv(const ElemType type, const unsigned int, const unsigned int, const Point &)) { - // cheat by using finite difference approximations: - const Real eps = 1.e-5; - Point pp = p, pm = p; - unsigned int deriv_j = 0; - - switch (j) - { - // d^2() / dxi^2 - case 0: - { - pp(0) += eps; - pm(0) -= eps; - deriv_j = 0; - break; - } - - // d^2() / dxi deta - case 1: - { - pp(1) += eps; - pm(1) -= eps; - deriv_j = 0; - break; - } - - // d^2() / deta^2 - case 2: - { - pp(1) += eps; - pm(1) -= eps; - deriv_j = 1; - break; - } - - // d^2()/dxidzeta - case 3: - { - pp(2) += eps; - pm(2) -= eps; - deriv_j = 0; - break; - } // d^2()/deta^2 + auto [pp, pm, eps, deriv_j] = fdm_second_points(j, p); - // d^2()/detadzeta - case 4: - { - pp(2) += eps; - pm(2) -= eps; - deriv_j = 1; - break; - } + return (deriv_func(type, order, i, deriv_j, pp) - + deriv_func(type, order, i, deriv_j, pm))/2./eps; +} - // d^2()/dzeta^2 - case 5: - { - pp(2) += eps; - pm(2) -= eps; - deriv_j = 2; - break; - } - default: - libmesh_error_msg("Invalid shape function derivative j = " << j); - } +template +OutputShape +fe_fdm_second_deriv(const ElemType type, + const Order order, + const Elem * elem, + const unsigned int i, + const unsigned int j, + const Point & p, + OutputShape(*deriv_func) + (const ElemType, const Order, + const Elem *, + const unsigned int, const unsigned int, + const Point &)) +{ + auto [pp, pm, eps, deriv_j] = fdm_second_points(j, p); - return (deriv_func(type, order, i, deriv_j, pp) - - deriv_func(type, order, i, deriv_j, pm))/2./eps; + return (deriv_func(type, order, elem, i, deriv_j, pp) - + deriv_func(type, order, elem, i, deriv_j, pm))/2./eps; } @@ -1391,6 +1357,13 @@ Real fe_fdm_deriv(const ElemType, const Order, const unsigned int, (const ElemType, const Order, const unsigned int, const Point &)); +template +Real fe_fdm_deriv(const ElemType, const Order, const Elem *, + const unsigned int, const unsigned int, const Point &, + Real(*shape_func) + (const ElemType, const Order, const Elem *, + const unsigned int, const Point &)); + template RealGradient fe_fdm_deriv(const Elem *, const Order, const unsigned int, @@ -1415,6 +1388,14 @@ fe_fdm_second_deriv(const Elem *, const Order, const unsigned int, (const Elem *, const Order, const unsigned int, const unsigned int, const Point &, const bool)); +template +Real +fe_fdm_second_deriv(const ElemType, const Order, const Elem *, + const unsigned int, const unsigned int, const Point &, + Real(*shape_func) + (const ElemType, const Order, const Elem *, + const unsigned int, const unsigned int, const Point &)); + template RealGradient fe_fdm_second_deriv(const Elem *, const Order, const unsigned int, diff --git a/src/fe/fe_interface.C b/src/fe/fe_interface.C index 77bd32c1ffa..25b5c2c0020 100644 --- a/src/fe/fe_interface.C +++ b/src/fe/fe_interface.C @@ -2193,6 +2193,8 @@ unsigned int FEInterface::max_order(const FEType & fe_t, return 2; case PYRAMID18: return 3; + case C0POLYHEDRON: + return 1; default: return unknown; } @@ -2234,6 +2236,7 @@ unsigned int FEInterface::max_order(const FEType & fe_t, case PYRAMID13: case PYRAMID14: case PYRAMID18: + case C0POLYHEDRON: return unlimited; default: return unknown; @@ -2283,6 +2286,7 @@ unsigned int FEInterface::max_order(const FEType & fe_t, case PYRAMID13: case PYRAMID14: case PYRAMID18: + case C0POLYHEDRON: return 0; default: return unknown; @@ -2326,6 +2330,7 @@ unsigned int FEInterface::max_order(const FEType & fe_t, case PYRAMID13: case PYRAMID14: case PYRAMID18: + case C0POLYHEDRON: return 0; default: return unknown; @@ -2364,6 +2369,7 @@ unsigned int FEInterface::max_order(const FEType & fe_t, case PYRAMID13: case PYRAMID14: case PYRAMID18: + case C0POLYHEDRON: return unlimited; default: return unknown; @@ -2404,6 +2410,7 @@ unsigned int FEInterface::max_order(const FEType & fe_t, case PYRAMID13: case PYRAMID14: case PYRAMID18: + case C0POLYHEDRON: return 0; default: return unknown; @@ -2448,6 +2455,7 @@ unsigned int FEInterface::max_order(const FEType & fe_t, case PYRAMID13: case PYRAMID14: case PYRAMID18: + case C0POLYHEDRON: return 0; default: return unknown; @@ -2497,6 +2505,7 @@ unsigned int FEInterface::max_order(const FEType & fe_t, case PYRAMID13: case PYRAMID14: case PYRAMID18: + case C0POLYHEDRON: return 0; default: return unknown; @@ -2545,6 +2554,7 @@ unsigned int FEInterface::max_order(const FEType & fe_t, case PYRAMID13: case PYRAMID14: case PYRAMID18: + case C0POLYHEDRON: return 0; default: return unknown; diff --git a/src/fe/fe_lagrange.C b/src/fe/fe_lagrange.C index 446eb75600f..88cf0fb6793 100644 --- a/src/fe/fe_lagrange.C +++ b/src/fe/fe_lagrange.C @@ -504,10 +504,11 @@ unsigned int lagrange_n_dofs(const ElemType t, const Elem * e, const Order o) return 0; case C0POLYGON: - // C0Polygon requires using newer FE APIs + case C0POLYHEDRON: + // Polygons and polyhedra require using newer FE APIs if (!e) libmesh_error_msg("Code (see stack trace) used an outdated FE function overload.\n" - "n_dofs() on a C0Polygon is not defined by its ElemType alone."); + "n_dofs() on a polygon or polyhedron is not defined by ElemType alone."); return e->n_nodes(); default: @@ -693,6 +694,7 @@ unsigned int lagrange_n_dofs_at_node(const ElemType t, } case C0POLYGON: + case C0POLYHEDRON: return 1; diff --git a/src/fe/fe_lagrange_shape_3D.C b/src/fe/fe_lagrange_shape_3D.C index af32c1f555a..c06eee9322c 100644 --- a/src/fe/fe_lagrange_shape_3D.C +++ b/src/fe/fe_lagrange_shape_3D.C @@ -21,6 +21,7 @@ #include "libmesh/elem.h" #include "libmesh/fe_lagrange_shape_1D.h" #include "libmesh/enum_to_string.h" +#include "libmesh/cell_c0polyhedron.h" // Anonymous namespace for functions shared by LAGRANGE and // L2_LAGRANGE implementations. Implementations appear at the bottom @@ -32,12 +33,14 @@ using namespace libMesh; template Real fe_lagrange_3D_shape(const ElemType, const Order order, + const Elem * elem, const unsigned int i, const Point & p); template Real fe_lagrange_3D_shape_deriv(const ElemType type, const Order order, + const Elem * elem, const unsigned int i, const unsigned int j, const Point & p); @@ -47,6 +50,7 @@ Real fe_lagrange_3D_shape_deriv(const ElemType type, template Real fe_lagrange_3D_shape_second_deriv(const ElemType type, const Order order, + const Elem * elem, const unsigned int i, const unsigned int j, const Point & p); @@ -416,7 +420,7 @@ Real FE<3,LAGRANGE>::shape(const ElemType type, const unsigned int i, const Point & p) { - return fe_lagrange_3D_shape(type, order, i, p); + return fe_lagrange_3D_shape(type, order, nullptr, i, p); } @@ -427,7 +431,7 @@ Real FE<3,L2_LAGRANGE>::shape(const ElemType type, const unsigned int i, const Point & p) { - return fe_lagrange_3D_shape(type, order, i, p); + return fe_lagrange_3D_shape(type, order, nullptr, i, p); } @@ -442,7 +446,7 @@ Real FE<3,LAGRANGE>::shape(const Elem * elem, libmesh_assert(elem); // call the orientation-independent shape functions - return fe_lagrange_3D_shape(elem->type(), order + add_p_level*elem->p_level(), i, p); + return fe_lagrange_3D_shape(elem->type(), order + add_p_level*elem->p_level(), elem, i, p); } @@ -457,7 +461,7 @@ Real FE<3,L2_LAGRANGE>::shape(const Elem * elem, libmesh_assert(elem); // call the orientation-independent shape functions - return fe_lagrange_3D_shape(elem->type(), order + add_p_level*elem->p_level(), i, p); + return fe_lagrange_3D_shape(elem->type(), order + add_p_level*elem->p_level(), elem, i, p); } @@ -470,7 +474,7 @@ Real FE<3,LAGRANGE>::shape(const FEType fet, const bool add_p_level) { libmesh_assert(elem); - return fe_lagrange_3D_shape(elem->type(), fet.order + add_p_level*elem->p_level(), i, p); + return fe_lagrange_3D_shape(elem->type(), fet.order + add_p_level*elem->p_level(), elem, i, p); } @@ -483,7 +487,7 @@ Real FE<3,L2_LAGRANGE>::shape(const FEType fet, const bool add_p_level) { libmesh_assert(elem); - return fe_lagrange_3D_shape(elem->type(), fet.order + add_p_level*elem->p_level(), i, p); + return fe_lagrange_3D_shape(elem->type(), fet.order + add_p_level*elem->p_level(), elem, i, p); } template <> @@ -493,7 +497,7 @@ Real FE<3,LAGRANGE>::shape_deriv(const ElemType type, const unsigned int j, const Point & p) { - return fe_lagrange_3D_shape_deriv(type, order, i, j, p); + return fe_lagrange_3D_shape_deriv(type, order, nullptr, i, j, p); } @@ -505,7 +509,7 @@ Real FE<3,L2_LAGRANGE>::shape_deriv(const ElemType type, const unsigned int j, const Point & p) { - return fe_lagrange_3D_shape_deriv(type, order, i, j, p); + return fe_lagrange_3D_shape_deriv(type, order, nullptr, i, j, p); } @@ -521,7 +525,7 @@ Real FE<3,LAGRANGE>::shape_deriv(const Elem * elem, libmesh_assert(elem); // call the orientation-independent shape function derivatives - return fe_lagrange_3D_shape_deriv(elem->type(), order + add_p_level*elem->p_level(), i, j, p); + return fe_lagrange_3D_shape_deriv(elem->type(), order + add_p_level*elem->p_level(), elem, i, j, p); } @@ -536,7 +540,7 @@ Real FE<3,L2_LAGRANGE>::shape_deriv(const Elem * elem, libmesh_assert(elem); // call the orientation-independent shape function derivatives - return fe_lagrange_3D_shape_deriv(elem->type(), order + add_p_level*elem->p_level(), i, j, p); + return fe_lagrange_3D_shape_deriv(elem->type(), order + add_p_level*elem->p_level(), elem, i, j, p); } @@ -549,7 +553,7 @@ Real FE<3,LAGRANGE>::shape_deriv(const FEType fet, const bool add_p_level) { libmesh_assert(elem); - return fe_lagrange_3D_shape_deriv(elem->type(), fet.order + add_p_level*elem->p_level(), i, j, p); + return fe_lagrange_3D_shape_deriv(elem->type(), fet.order + add_p_level*elem->p_level(), elem, i, j, p); } @@ -562,7 +566,7 @@ Real FE<3,L2_LAGRANGE>::shape_deriv(const FEType fet, const bool add_p_level) { libmesh_assert(elem); - return fe_lagrange_3D_shape_deriv(elem->type(), fet.order + add_p_level*elem->p_level(), i, j, p); + return fe_lagrange_3D_shape_deriv(elem->type(), fet.order + add_p_level*elem->p_level(), elem, i, j, p); } #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES @@ -574,7 +578,7 @@ Real FE<3,LAGRANGE>::shape_second_deriv(const ElemType type, const unsigned int j, const Point & p) { - return fe_lagrange_3D_shape_second_deriv(type, order, i, j, p); + return fe_lagrange_3D_shape_second_deriv(type, order, nullptr, i, j, p); } @@ -586,7 +590,7 @@ Real FE<3,L2_LAGRANGE>::shape_second_deriv(const ElemType type, const unsigned int j, const Point & p) { - return fe_lagrange_3D_shape_second_deriv(type, order, i, j, p); + return fe_lagrange_3D_shape_second_deriv(type, order, nullptr, i, j, p); } @@ -603,7 +607,7 @@ Real FE<3,LAGRANGE>::shape_second_deriv(const Elem * elem, // call the orientation-independent shape function derivatives return fe_lagrange_3D_shape_second_deriv - (elem->type(), order + add_p_level*elem->p_level(), i, j, p); + (elem->type(), order + add_p_level*elem->p_level(), elem, i, j, p); } @@ -620,7 +624,7 @@ Real FE<3,L2_LAGRANGE>::shape_second_deriv(const Elem * elem, // call the orientation-independent shape function derivatives return fe_lagrange_3D_shape_second_deriv - (elem->type(), order + add_p_level*elem->p_level(), i, j, p); + (elem->type(), order + add_p_level*elem->p_level(), elem, i, j, p); } @@ -634,7 +638,7 @@ Real FE<3,LAGRANGE>::shape_second_deriv(const FEType fet, { libmesh_assert(elem); return fe_lagrange_3D_shape_second_deriv - (elem->type(), fet.order + add_p_level*elem->p_level(), i, j, p); + (elem->type(), fet.order + add_p_level*elem->p_level(), elem, i, j, p); } @@ -649,7 +653,7 @@ Real FE<3,L2_LAGRANGE>::shape_second_deriv(const FEType fet, { libmesh_assert(elem); return fe_lagrange_3D_shape_second_deriv - (elem->type(), fet.order + add_p_level*elem->p_level(), i, j, p); + (elem->type(), fet.order + add_p_level*elem->p_level(), elem, i, j, p); } @@ -666,6 +670,7 @@ using namespace libMesh; template Real fe_lagrange_3D_shape(const ElemType type, const Order order, + const Elem * elem, const unsigned int i, const Point & p) { @@ -790,6 +795,41 @@ Real fe_lagrange_3D_shape(const ElemType type, } } + case C0POLYHEDRON: + { + // Polyhedra require using newer FE APIs + if (!elem) + libmesh_error_msg("Code (see stack trace) used an outdated FE function overload.\n" + "Shape functions on a polyhedron are not defined by ElemType alone."); + + libmesh_assert(elem->type() == C0POLYHEDRON); + + const C0Polyhedron & poly = *cast_ptr(elem); + + // We can't use a small tolerance here, because in + // inverse_map() Newton might hand us intermediate + // iterates outside the polyhedron. + const auto [s, a, b, c] = poly.subelement_coordinates(p, 100); + if (s == invalid_uint) + return 0; + libmesh_assert_less(s, poly.n_subelements()); + + const auto subtet = poly.subelement(s); + + // Avoid signed/unsigned comparison warnings + const int nodei = i; + if (nodei == subtet[0]) + return 1-a-b-c; + if (nodei == subtet[1]) + return a; + if (nodei == subtet[2]) + return b; + if (nodei == subtet[3]) + return c; + + // Basis function i is not supported on p's subtet + return 0; + } default: libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(type)); @@ -1434,7 +1474,7 @@ Real fe_lagrange_3D_shape(const ElemType type, } #else // LIBMESH_DIM != 3 - libmesh_ignore(type, order, i, p); + libmesh_ignore(type, order, elem, i, p); libmesh_not_implemented(); #endif } @@ -1444,6 +1484,7 @@ Real fe_lagrange_3D_shape(const ElemType type, template Real fe_lagrange_3D_shape_deriv(const ElemType type, const Order order, + const Elem * elem, const unsigned int i, const unsigned int j, const Point & p) @@ -1729,6 +1770,94 @@ Real fe_lagrange_3D_shape_deriv(const ElemType type, } } + case C0POLYHEDRON: + { + // Polyhedra require using newer FE APIs + if (!elem) + libmesh_error_msg("Code (see stack trace) used an outdated FE function overload.\n" + "Shape functions on a polyhedron are not defined by ElemType alone."); + + libmesh_assert(elem->type() == C0POLYHEDRON); + + const C0Polyhedron & poly = *cast_ptr(elem); + + // We can't use a small tolerance here, because in + // inverse_map() Newton might hand us intermediate + // iterates outside the polyhedron. + const auto [s, a, b, c] = poly.subelement_coordinates(p, 100); + if (s == invalid_uint) + return 0; + libmesh_assert_less(s, poly.n_subelements()); + + const auto subtet = poly.subelement(s); + + // Find derivatives w.r.t. subtriangle barycentric + // coordinates + Real du_da = 0, du_db = 0, du_dc = 0; + + // Avoid signed/unsigned comparison warnings + const int nodei = i; + if (nodei == subtet[0]) + du_da = du_db = du_dc = -1; + else if (nodei == subtet[1]) + du_da = 1; + else if (nodei == subtet[2]) + du_db = 1; + else if (nodei == subtet[3]) + du_dc = 1; + else + // Basis function i is not supported on p's subtet + return 0; + + // We want to return derivatives with respect to + // xi/eta/zeta for the polyhedron, but what we + // calculated above are with respect to xi and eta + // coordinates for a master *triangle*. We need to + // convert from one to the other. + + const auto master_points = poly.master_subelement(s); + + const RealTensor dXi_dA( + master_points[1](0) - master_points[0](0), master_points[2](0) - master_points[0](0), master_points[3](0) - master_points[0](0), + master_points[1](1) - master_points[0](1), master_points[2](1) - master_points[0](1), master_points[3](1) - master_points[0](1), + master_points[1](2) - master_points[0](2), master_points[2](2) - master_points[0](2), master_points[3](2) - master_points[0](2)); + + // When we vectorize this we'll want a full inverse, but + // when we're querying one component at a time it's + // cheaper to manually compute a single column. + // const RealTensor dabc_dxietazeta_dabc = dxietazeta_dabc.inverse(); + const Real jac = dXi_dA.det(); + + switch (j) + { + // d()/dxi + case 0: + { + const Real da_dxi = (dXi_dA(2,2)*dXi_dA(1,1) - dXi_dA(2,1)*dXi_dA(1,2)) / jac; + const Real db_dxi = -(dXi_dA(2,2)*dXi_dA(1,0) - dXi_dA(2,0)*dXi_dA(1,2)) / jac; + const Real dc_dxi = (dXi_dA(2,1)*dXi_dA(1,0) - dXi_dA(2,0)*dXi_dA(1,1)) / jac; + return du_da*da_dxi + du_db*db_dxi + du_dc*dc_dxi; + } + // d()/deta + case 1: + { + const Real da_deta = -(dXi_dA(2,2)*dXi_dA(0,1) - dXi_dA(2,1)*dXi_dA(0,2) ) / jac; + const Real db_deta = (dXi_dA(2,2)*dXi_dA(0,0) - dXi_dA(2,0)*dXi_dA(0,2)) / jac; + const Real dc_deta = -(dXi_dA(2,1)*dXi_dA(0,0) - dXi_dA(2,0)*dXi_dA(0,1)) / jac; + return du_da*da_deta + du_db*db_deta + du_dc*dc_deta; + } + // d()/dzeta + case 2: + { + const Real da_dzeta = (dXi_dA(1,2)*dXi_dA(0,1) - dXi_dA(1,1)*dXi_dA(0,2) ) / jac; + const Real db_dzeta = -(dXi_dA(1,2)*dXi_dA(0,0) - dXi_dA(1,0)*dXi_dA(0,2)) / jac; + const Real dc_dzeta = (dXi_dA(1,1)*dXi_dA(0,0) - dXi_dA(1,0)*dXi_dA(0,1)) / jac; + return du_da*da_dzeta + du_db*db_dzeta + du_dc*dc_dzeta; + } + default: + libmesh_error_msg("ERROR: Invalid derivative index j = " << j); + } + } default: libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(type)); @@ -3086,7 +3215,7 @@ Real fe_lagrange_3D_shape_deriv(const ElemType type, case PYRAMID18: { - return fe_fdm_deriv(type, order, i, j, p, fe_lagrange_3D_shape); + return fe_fdm_deriv(type, order, elem, i, j, p, fe_lagrange_3D_shape); } default: @@ -3100,7 +3229,7 @@ Real fe_lagrange_3D_shape_deriv(const ElemType type, } #else // LIBMESH_DIM != 3 - libmesh_ignore(type, order, i, j, p); + libmesh_ignore(type, order, elem, i, j, p); libmesh_not_implemented(); #endif } @@ -3112,6 +3241,7 @@ Real fe_lagrange_3D_shape_deriv(const ElemType type, template Real fe_lagrange_3D_shape_second_deriv(const ElemType type, const Order order, + const Elem * elem, const unsigned int i, const unsigned int j, const Point & p) @@ -3325,6 +3455,14 @@ Real fe_lagrange_3D_shape_second_deriv(const ElemType type, } } + // All second derivatives for piecewise-linear polyhedra are + // zero or dirac-type distributions, but we can't put the + // latter in a Real, so beware when integrating... + case C0POLYHEDRON: + { + return 0.; + } + default: libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(type)); } @@ -5157,7 +5295,7 @@ Real fe_lagrange_3D_shape_second_deriv(const ElemType type, case PYRAMID18: { - return fe_fdm_second_deriv(type, order, i, j, p, + return fe_fdm_second_deriv(type, order, elem, i, j, p, fe_lagrange_3D_shape_deriv); } @@ -5172,7 +5310,7 @@ Real fe_lagrange_3D_shape_second_deriv(const ElemType type, } #else // LIBMESH_DIM != 3 - libmesh_ignore(type, order, i, j, p); + libmesh_ignore(type, order, elem, i, j, p); libmesh_not_implemented(); #endif } diff --git a/src/fe/fe_monomial.C b/src/fe/fe_monomial.C index e1077e9fd4e..8d5a6b7ef09 100644 --- a/src/fe/fe_monomial.C +++ b/src/fe/fe_monomial.C @@ -87,6 +87,7 @@ unsigned int monomial_n_dofs(const ElemType t, const Order o) case PYRAMID13: case PYRAMID14: case PYRAMID18: + case C0POLYHEDRON: return 4; case INVALID_ELEM: @@ -140,6 +141,7 @@ unsigned int monomial_n_dofs(const ElemType t, const Order o) case PYRAMID13: case PYRAMID14: case PYRAMID18: + case C0POLYHEDRON: return 10; case INVALID_ELEM: @@ -193,6 +195,7 @@ unsigned int monomial_n_dofs(const ElemType t, const Order o) case PYRAMID13: case PYRAMID14: case PYRAMID18: + case C0POLYHEDRON: return 20; case INVALID_ELEM: @@ -244,6 +247,7 @@ unsigned int monomial_n_dofs(const ElemType t, const Order o) case PYRAMID5: case PYRAMID13: case PYRAMID14: + case C0POLYHEDRON: return 35; case INVALID_ELEM: @@ -294,6 +298,7 @@ unsigned int monomial_n_dofs(const ElemType t, const Order o) case PYRAMID5: case PYRAMID13: case PYRAMID14: + case C0POLYHEDRON: return (order+1)*(order+2)*(order+3)/6; case INVALID_ELEM: diff --git a/src/geom/cell_c0polyhedron.C b/src/geom/cell_c0polyhedron.C new file mode 100644 index 00000000000..14ee98877d9 --- /dev/null +++ b/src/geom/cell_c0polyhedron.C @@ -0,0 +1,798 @@ +// The libMesh Finite Element Library. +// Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner + +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. + +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. + +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +// Local includes +#include "libmesh/cell_c0polyhedron.h" + +#include "libmesh/enum_order.h" +#include "libmesh/face_polygon.h" +#include "libmesh/mesh_tools.h" +#include "libmesh/replicated_mesh.h" +#include "libmesh/side.h" +#include "libmesh/tensor_value.h" + +// C++ headers +#include // std::iota + +namespace libMesh +{ + +// ------------------------------------------------------------ +// C0Polyhedron class member functions + +C0Polyhedron::C0Polyhedron + (const std::vector> & sides, Elem * p) : + Polyhedron(sides, p) +{ + this->retriangulate(); +} + + + +C0Polyhedron::~C0Polyhedron() = default; + + + +std::unique_ptr C0Polyhedron::disconnected_clone() const +{ + auto sides = this->side_clones(); + + std::unique_ptr returnval = + std::make_unique(sides); + + returnval->set_id() = this->id(); +#ifdef LIBMESH_ENABLE_UNIQUE_ID + if (this->valid_unique_id()) + returnval->set_unique_id(this->unique_id()); +#endif + returnval->subdomain_id() = this->subdomain_id(); + returnval->processor_id() = this->processor_id(); + + const auto n_elem_ints = this->n_extra_integers(); + returnval->add_extra_integers(n_elem_ints); + for (unsigned int i = 0; i != n_elem_ints; ++i) + returnval->set_extra_integer(i, this->get_extra_integer(i)); + + returnval->set_mapping_type(this->mapping_type()); + returnval->set_mapping_data(this->mapping_data()); + + return returnval; +} + + + +bool C0Polyhedron::is_vertex(const unsigned int libmesh_dbg_var(i)) const +{ + libmesh_assert (i < this->n_nodes()); + + return true; +} + + + +bool C0Polyhedron::is_edge(const unsigned int libmesh_dbg_var(i)) const +{ + libmesh_assert_less (i, this->n_nodes()); + + return false; +} + + + +bool C0Polyhedron::is_face(const unsigned int libmesh_dbg_var(i)) const +{ + libmesh_assert_less (i, this->n_nodes()); + + return false; +} + + + +bool C0Polyhedron::is_node_on_side(const unsigned int n, + const unsigned int s) const +{ + libmesh_assert_less (n, this->n_nodes()); + libmesh_assert_less (s, this->n_sides()); + + const std::vector & node_map = + std::get<2>(_sidelinks_data[s]); + + const auto it = std::find(node_map.begin(), node_map.end(), n); + + return (it != node_map.end()); +} + +std::vector +C0Polyhedron::nodes_on_side(const unsigned int s) const +{ + libmesh_assert_less(s, this->n_sides()); + return std::get<2>(_sidelinks_data[s]); +} + +std::vector +C0Polyhedron::nodes_on_edge(const unsigned int e) const +{ + auto [s, se] = _edge_lookup[e]; + const Polygon & face = *std::get<0>(_sidelinks_data[s]); + const std::vector & node_map = + std::get<2>(_sidelinks_data[s]); + std::vector nodes_on_edge = + face.nodes_on_side(se); + for (auto i : index_range(nodes_on_edge)) + nodes_on_edge[i] = node_map[nodes_on_edge[i]]; + return nodes_on_edge; +} + + + +bool C0Polyhedron::is_node_on_edge(const unsigned int n, + const unsigned int e) const +{ + libmesh_assert_less(e, _edge_lookup.size()); + auto [s, se] = _edge_lookup[e]; + + const Polygon & face = *std::get<0>(_sidelinks_data[s]); + const std::vector & node_map = + std::get<2>(_sidelinks_data[s]); + std::vector nodes_on_edge = + face.nodes_on_side(se); + for (auto noe : nodes_on_edge) + if (node_map[noe] == n) + return true; + + return false; +} + + + +void C0Polyhedron::connectivity(const unsigned int /*sf*/, + const IOPackage /*iop*/, + std::vector & /*conn*/) const +{ + libmesh_not_implemented(); +} + + + +Real C0Polyhedron::volume () const +{ + // This specialization is good for Lagrange mappings only + if (this->mapping_type() != LAGRANGE_MAP) + return this->Elem::volume(); + + // We use a triangulation to calculate here. + + Real six_vol = 0; + for (const auto & subtet : this->_triangulation) + { + const Point p0 = this->point(subtet[0]); + const Point p1 = this->point(subtet[1]); + const Point p2 = this->point(subtet[2]); + const Point p3 = this->point(subtet[3]); + + const Point v01 = p1 - p0; + const Point v02 = p2 - p0; + const Point v03 = p3 - p0; + + six_vol += triple_product(v01, v02, v03); + } + + return six_vol * (1./6.); +} + + + +Point C0Polyhedron::true_centroid () const +{ + // This specialization is good for Lagrange mappings only + if (this->mapping_type() != LAGRANGE_MAP) + return this->Elem::true_centroid(); + + Real six_vol = 0; + Point six_vol_weighted_centroid; + for (const auto & subtet : this->_triangulation) + { + const Point p0 = this->point(subtet[0]); + const Point p1 = this->point(subtet[1]); + const Point p2 = this->point(subtet[2]); + const Point p3 = this->point(subtet[3]); + + const Point v01 = p1 - p0; + const Point v02 = p2 - p0; + const Point v03 = p3 - p0; + + const Real six_tet_vol = triple_product(v01, v02, v03); + + const Point tet_centroid = (p0 + p1 + p2 + p3) * 0.25; + + six_vol += six_tet_vol; + + six_vol_weighted_centroid += six_tet_vol * tet_centroid; + } + + return six_vol_weighted_centroid / six_vol; +} + + + +std::pair +C0Polyhedron::second_order_child_vertex (const unsigned int /*n*/) const +{ + libmesh_not_implemented(); + return std::pair (0,0); +} + + + +ElemType C0Polyhedron::side_type (const unsigned int libmesh_dbg_var(s)) const +{ + libmesh_assert_less (s, this->n_sides()); + return C0POLYGON; +} + + + +void C0Polyhedron::retriangulate() +{ + this->_triangulation.clear(); + + // We should already have a triangulation for each side. We'll turn + // that into a triangulation of the entire polyhedral surface + // (stored as a mesh, so we can use `find_neighbors()`, then turn + // that into a tetrahedralization by shrinking the volume one tet at + // a time, which should work fine for convex polyhedra. + + Parallel::Communicator comm_self; // the default communicator is 1 rank + ReplicatedMesh surface(comm_self); + + // poly_node_to_local_id[poly_node] is the local id of + // poly_node in the Polyhedron, which is also the global id of the + // corresponding node in the surface mesh. + std::unordered_map poly_node_to_id; + + for (unsigned int s : make_range(this->n_sides())) + { + const auto & [side, inward_normal, node_map] = this->_sidelinks_data[s]; + + for (auto t : make_range(side->n_subtriangles())) + { + Elem * tri = surface.add_elem(Elem::build(TRI3)); + + const std::array subtri = side->subtriangle(t); + + for (int i : make_range(3)) + { + const int side_id = subtri[i]; + const Node * poly_node = side->node_ptr(side_id); + + libmesh_assert_less(side_id, node_map.size()); + const unsigned int local_id = node_map[side_id]; + + Node * surf_node = surface.query_node_ptr(local_id); + if (surf_node) + libmesh_assert_equal_to(*(const Point*)poly_node, + *(const Point*)surf_node); + else + surf_node = surface.add_point(*poly_node, local_id); + + // Get a consistent orientation for surface triangles, + // facing their zeta directions outward + const int tri_node = inward_normal ? i : 2-i; + tri->set_node(tri_node, surf_node); + } + } + } + + surface.allow_renumbering(false); + surface.prepare_for_use(); + + // We should have a watertight surface with consistent triangle + // orientations. That's expensive to check. +#ifdef DEBUG + auto verify_surface = [& surface] () + { + for (const Elem * elem : surface.element_ptr_range()) + { + for (auto s : make_range(3)) + { + const Elem * neigh = elem->neighbor_ptr(s); + libmesh_assert(neigh); + libmesh_assert_equal_to(neigh, + surface.elem_ptr(neigh->id())); + const unsigned int ns = neigh->which_neighbor_am_i(elem); + libmesh_assert_less(ns, 3); + libmesh_assert_equal_to(elem->node_ptr(s), + neigh->node_ptr((ns+1)%3)); + libmesh_assert_equal_to(elem->node_ptr((s+1)%3), + neigh->node_ptr(ns)); + } + } + }; + + verify_surface(); +#endif + + // We'll have to edit this as we change the surface elements, but we + // have a method to initialize it easily. + std::vector> nodes_to_elem_vec_map; + MeshTools::build_nodes_to_elem_map(surface, nodes_to_elem_vec_map); + + // We'll be inserting and deleting entries heavily, so we'll use + // sets rather than vectors. We want to get the same results in + // parallel, so we'll use element ids rather than Elem pointers + std::vector> nodes_to_elem_map; + for (auto i : index_range(nodes_to_elem_vec_map)) + nodes_to_elem_map.emplace_back + (nodes_to_elem_vec_map[i].begin(), + nodes_to_elem_vec_map[i].end()); + + // Now start meshing the volume enclosed by surface, one tet at a + // time, with a greedy heuristic: find the vertex node with the most + // acutely convex (solid) angle and strip it out as tetrahedra. + + // We'll want a vector of surrounding nodes for multiple uses, + // sometimes with a similarly-sorted vector of surrounding elements + // to go with it. + auto surroundings_of = + [&nodes_to_elem_map, & surface] + (const Node & node, + std::vector * surrounding_elems) + { + const std::set & elems_by_node = + nodes_to_elem_map[node.id()]; + + const unsigned int n_surrounding = elems_by_node.size(); + libmesh_assert_greater_equal(n_surrounding, 3); + + if (surrounding_elems) + { + libmesh_assert(surrounding_elems->empty()); + surrounding_elems->resize(n_surrounding); + } + + std::vector surrounding_nodes(n_surrounding); + + Elem * elem = surface.elem_ptr(*elems_by_node.begin()); + for (auto i : make_range(n_surrounding)) + { + const unsigned int n = elem->get_node_index(&node); + libmesh_assert_not_equal_to(n, invalid_uint); + Node * next_node = elem->node_ptr((n+1)%3); + surrounding_nodes[i] = next_node; + if (surrounding_elems) + (*surrounding_elems)[i] = elem; + elem = elem->neighbor_ptr((n+2)%3); + libmesh_assert(elem); + libmesh_assert_equal_to(elem, surface.elem_ptr(elem->id())); + + // We should have a manifold here, but verifying that is + // expensive +#ifdef DEBUG + libmesh_assert_equal_to + (std::count(surrounding_nodes.begin(), + surrounding_nodes.end(), next_node), + 1); +#endif + } + + // We should have finished a loop + libmesh_assert_equal_to + (elem, surface.elem_ptr(*elems_by_node.begin())); + + return surrounding_nodes; + }; + + auto geometry_at = [&surroundings_of](const Node & node) + { + const std::vector surrounding_nodes = + surroundings_of(node, nullptr); + + // Now sum up solid angles from tetrahedra created from the + // trivial triangulation of the surrounding nodes loop + Real total_solid_angle = 0; + const int n_surrounding = + cast_int(surrounding_nodes.size()); + + for (auto n : make_range(n_surrounding-2)) + { + const Point + v01 = static_cast(*surrounding_nodes[n]) - node, + v02 = static_cast(*surrounding_nodes[n+1]) - node, + v03 = static_cast(*surrounding_nodes[n+2]) - node; + + total_solid_angle += solid_angle(v01, v02, v03); + } + + return std::make_pair(n_surrounding, total_solid_angle); + }; + + // We'll keep track of solid angles and node valences so we don't + // waste time recomputing them when they haven't changed. We need + // to be able to search efficiently for the smallest angles of each + // valence, but also search efficiently for a particular node to + // remove and reinsert it when its connectivity changes. + // + // Since C++11 multimap has guaranteed that pairs with matching keys + // are kept in insertion order, so we can use Node * for values even + // in parallel. + typedef std::multimap, Node*> node_map_type; + node_map_type nodes_by_geometry; + std::map node_index; + + for (auto node : surface.node_ptr_range()) + node_index[node] = + nodes_by_geometry.emplace(geometry_at(*node), node); + + // In 3D, this will require nested loops: an outer loop to remove + // each vertex, and an inner loop to remove multiple tetrahedra in + // cases where the vertex has more than 3 neighboring triangles. + + // We'll be done when there are only three "unremoved" nodes left, + // so they don't actually enclose any volume. + for (auto i : make_range(nodes_by_geometry.size()-3)) + { + auto geometry_it = nodes_by_geometry.begin(); + auto geometry_key = geometry_it->first; + auto [valence, angle] = geometry_key; + Node * node = geometry_it->second; + libmesh_ignore(i); + + // If our lowest-valence nodes are all points of non-convexity, + // skip to a higher valence. + while (angle > 2*pi-TOLERANCE) + { + geometry_it = + nodes_by_geometry.upper_bound + (std::make_pair(valence, Real(100))); + libmesh_assert(geometry_it != nodes_by_geometry.end()); + + std::tie(geometry_key, node) = *geometry_it; + std::tie(valence, angle) = geometry_key; + } + + std::vector surrounding_elems; + std::vector surrounding_nodes = + surroundings_of(*node, &surrounding_elems); + + const std::size_t n_surrounding = surrounding_nodes.size(); + + // As we separate surrounding nodes from our center node, we'll + // be marking them as nullptr; we still need to be able to find + // predecessor and successor nodes in order afterward. + auto find_valid_nodes_around = + [n_surrounding, & surrounding_nodes] + (unsigned int j) + { + unsigned int jnext = (j+1)%n_surrounding; + while (!surrounding_nodes[jnext]) + jnext = (jnext+1)%n_surrounding; + + unsigned int jprev = (j+n_surrounding-1)%n_surrounding; + while (!surrounding_nodes[jprev]) + jprev = (jprev+n_surrounding-1)%n_surrounding; + + return std::make_pair(jprev, jnext); + }; + + // We may have too many surrounding nodes to handle with + // just one tet. In that case we'll keep a cache of the + // element qualities that we'd get by making a tet with the + // edge from the center node to each surrounding node, so we + // can build the best tets first. + // + // In the case where we just have 3 nodes, we'll just pretend + // they all have the same positive quality, so we can still + // search this vector. + std::vector local_tet_quality(n_surrounding, 1); + + // From our center node with N surrounding nodes we can make N-2 + // tetrahedra. The first N-3 each replace two surface tets with + // two new surface tets. + // + // My first idea was to greedily pick nodes with the smallest + // local (solid) angles to get the best quality. This works in + // 2D, but such a node can give a pancake tet in 3D. + // + // My second idea was to greedily pick nodes with the highest + // prospective tet quality. This works for the first tets, but + // can leave a pancake tet behind. + // + // My third idea is to try to fix the lowest quality tets first, + // by picking cases where they have higher quality neighbors, + // and creating those neighbors so as to change them. + + auto find_new_tet_nodes = + [& local_tet_quality, & find_valid_nodes_around] + () + { + unsigned int jbest = 0; + auto [jminus, jplus] = find_valid_nodes_around(jbest); + Real qneighbest = std::min(local_tet_quality[jminus], + local_tet_quality[jplus]); + for (auto j : make_range(std::size_t(1), + local_tet_quality.size())) + { + // We don't want to build a bad tet + if (local_tet_quality[j] <= 0) + continue; + + std::tie(jminus, jplus) = find_valid_nodes_around(j); + Real qneighj = std::min(local_tet_quality[jminus], + local_tet_quality[jplus]); + + // We don't want to build a tet that can't fix a neighbor + // if we can build one that can. + if (qneighbest <= 0 && + qneighj > 0) + continue; + + // We want to try for the best possible fix. + if ((local_tet_quality[j] - qneighj) > + (local_tet_quality[jbest] - qneighj)) + { + jbest = j; + qneighbest = qneighj; + } + } + + libmesh_error_msg_if + (local_tet_quality[jbest] <= 0, + "Cannot build non-singular non-inverted tet"); + + std::tie(jminus, jplus) = find_valid_nodes_around(jbest); + + return std::make_tuple(jbest, jminus, jplus); + }; + + if (n_surrounding > 3) + { + // We'll be searching local_tet_quality even after tet + // extraction disconnects us from some nodes; when we do we + // don't want to get one. + constexpr Real far_node = -1e6; + + // Vectors from the center node to each of its surrounding + // nodes are helpful for calculating prospective tet + // quality. + std::vector v0s(n_surrounding); + for (auto j : make_range(n_surrounding)) + v0s[j] = *(Point *)surrounding_nodes[j] - *node; + + // Find the tet quality we'd potentially get from each + // possible choice of tet + auto local_tet_quality_of = + [& surrounding_nodes, & v0s, & find_valid_nodes_around] + (unsigned int j) + { + auto [jminus, jplus] = find_valid_nodes_around(j); + + // Anything proportional to the ratio of volume to + // total-edge-length-cubed should peak for perfect tets + // but hit 0 for pancakes and slivers. + + const Real total_len = + v0s[j].norm() + v0s[jminus].norm() + v0s[jplus].norm() + + (*(Point *)surrounding_nodes[jplus] - + *(Point *)surrounding_nodes[j]).norm() + + (*(Point *)surrounding_nodes[j] - + *(Point *)surrounding_nodes[jminus]).norm() + + (*(Point *)surrounding_nodes[jminus] - + *(Point *)surrounding_nodes[jplus]).norm(); + + // Orientation here is tricky. Think of the triple + // product as (v1 cross v2) dot v3, with right hand rule. + const Real six_vol = + triple_product(v0s[jminus], v0s[jplus], v0s[j]); + + return six_vol / (total_len * total_len * total_len); + }; + + for (auto j : make_range(n_surrounding)) + local_tet_quality[j] = local_tet_quality_of(j); + + // If we have N surrounding nodes, we can make N tets and + // that'll bring us back to the 3-surrounding-node case to + // finish. + for (auto t : make_range(n_surrounding-3)) + { + libmesh_ignore(t); + + auto [jbest, jminus, jplus] = find_new_tet_nodes(); + + // Turn these four nodes into a tet + Node * nbest = surrounding_nodes[jbest], + * nminus = surrounding_nodes[jminus], + * nplus = surrounding_nodes[jplus]; + this->add_tet(nminus->id(), nbest->id(), nplus->id(), + node->id()); + + // Replace the old two triangles around these nodes with the + // proper two new triangles. + Elem * oldtri1 = surrounding_elems[jminus], + * oldtri2 = surrounding_elems[jbest], + * newtri1 = surface.add_elem(Elem::build(TRI3)), + * newtri2 = surface.add_elem(Elem::build(TRI3)); + + const unsigned int c1 = oldtri1->get_node_index(node), + c2 = oldtri2->get_node_index(node); + + newtri1->set_node(0, node); + newtri1->set_node(1, nminus); + newtri1->set_node(2, nplus); + + surrounding_elems[jminus] = newtri1; + + Elem * neigh10 = oldtri1->neighbor_ptr(c1); + Elem * neigh12 = oldtri2->neighbor_ptr((c2+2)%3); + newtri1->set_neighbor(0, neigh10); + neigh10->set_neighbor(neigh10->which_neighbor_am_i(oldtri1), newtri1); + newtri1->set_neighbor(1, newtri2); + newtri1->set_neighbor(2, neigh12); + neigh12->set_neighbor(neigh12->which_neighbor_am_i(oldtri2), newtri1); + + newtri2->set_node(0, nplus); + newtri2->set_node(1, nminus); + newtri2->set_node(2, nbest); + + Elem * neigh21 = oldtri1->neighbor_ptr((c1+1)%3); + Elem * neigh22 = oldtri2->neighbor_ptr((c2+1)%3); + newtri2->set_neighbor(0, newtri1); + newtri2->set_neighbor(1, neigh21); + neigh21->set_neighbor(neigh21->which_neighbor_am_i(oldtri1), newtri2); + newtri2->set_neighbor(2, neigh22); + neigh22->set_neighbor(neigh22->which_neighbor_am_i(oldtri2), newtri2); + + for (unsigned int p : make_range(3)) + { + nodes_to_elem_map[oldtri1->node_id(p)].erase(oldtri1->id()); + nodes_to_elem_map[oldtri2->node_id(p)].erase(oldtri2->id()); + nodes_to_elem_map[newtri1->node_id(p)].insert(newtri1->id()); + nodes_to_elem_map[newtri2->node_id(p)].insert(newtri2->id()); + } + + // We've finished replacing the old triangles. + surface.delete_elem(oldtri1); + surface.delete_elem(oldtri2); + + // The solid angle for the far node should now stay + // unchanged until we're out of this inner loop; let's + // recalculate it here, and then we'll be done with it. + Node * & nbestref = surrounding_nodes[jbest]; + nodes_by_geometry.erase(node_index[nbestref]); + node_index[nbestref] = + nodes_by_geometry.emplace(geometry_at(*nbestref), nbestref); + + // The far node is no longer sharing an edge with our center + // node. Make sure we don't use it again with the center + // node. + local_tet_quality[jbest] = far_node; + nbestref = nullptr; + + // The potential tet qualities using the side nodes have + // changed now that they're directly connected to each + // other. + local_tet_quality[jminus] = + local_tet_quality_of(jminus); + + local_tet_quality[jplus] = + local_tet_quality_of(jplus); + } + } + + // Now we should have just 3 surrounding nodes, with which to + // make one tetrahedron. Put them in a counterclockwise + // (looking from outside) orientation, not the "best, clockwise, + // counterclockwise" we get from the lambda. + auto [j2, j1, j3] = find_new_tet_nodes(); + + // Turn these last four nodes into a tet + Node * n1 = surrounding_nodes[j1], + * n2 = surrounding_nodes[j2], + * n3 = surrounding_nodes[j3]; + this->add_tet(n1->id(), n2->id(), n3->id(), node->id()); + + // Replace the three surface triangles of that tet with the new + // fourth triangle. + Elem * oldtri1 = surrounding_elems[j1], + * oldtri2 = surrounding_elems[j2], + * oldtri3 = surrounding_elems[j3], + * newtri = surface.add_elem(Elem::build(TRI3)); + + const unsigned int c1 = oldtri1->get_node_index(node), + c2 = oldtri2->get_node_index(node), + c3 = oldtri3->get_node_index(node); + + newtri->set_node(0, n1); + newtri->set_node(1, n2); + newtri->set_node(2, n3); + Elem * neigh0 = oldtri1->neighbor_ptr((c1+1)%3); + newtri->set_neighbor(0, neigh0); + neigh0->set_neighbor(neigh0->which_neighbor_am_i(oldtri1), newtri); + Elem * neigh1 = oldtri2->neighbor_ptr((c2+1)%3); + newtri->set_neighbor(1, neigh1); + neigh1->set_neighbor(neigh1->which_neighbor_am_i(oldtri2), newtri); + Elem * neigh2 = oldtri3->neighbor_ptr((c3+1)%3); + newtri->set_neighbor(2, neigh2); + neigh2->set_neighbor(neigh2->which_neighbor_am_i(oldtri3), newtri); + + for (unsigned int p : make_range(3)) + { + nodes_to_elem_map[oldtri1->node_id(p)].erase(oldtri1->id()); + nodes_to_elem_map[oldtri2->node_id(p)].erase(oldtri2->id()); + nodes_to_elem_map[oldtri3->node_id(p)].erase(oldtri3->id()); + nodes_to_elem_map[newtri->node_id(p)].insert(newtri->id()); + } + + // We've finished replacing the old triangles. + surface.delete_elem(oldtri1); + surface.delete_elem(oldtri2); + surface.delete_elem(oldtri3); + + // We should have used up all our surrounding nodes now, and we + // shouldn't have messed up our surface in the process, and our + // center node should no longer be part of the surface. +#ifdef DEBUG + surrounding_nodes[j1] = nullptr; + surrounding_nodes[j2] = nullptr; + surrounding_nodes[j3] = nullptr; + + for (auto ltq : surrounding_nodes) + libmesh_assert(!ltq); + + if (surface.n_elem() > 3) + verify_surface(); + + for (const Elem * elem : surface.element_ptr_range()) + for (auto p : make_range(3)) + libmesh_assert_not_equal_to + (elem->node_ptr(p), node); +#endif + + // We've used up our center node, so it's not something we can + // eliminate again. + nodes_by_geometry.erase(geometry_it); + } + + // At this point our surface should just have two triangles left. + libmesh_assert_equal_to(surface.n_elem(), 2); +} + + +void C0Polyhedron::add_tet(int n1, + int n2, + int n3, + int n4) +{ +#ifndef NDEBUG + const auto nn = this->n_nodes(); + libmesh_assert_less(n1, nn); + libmesh_assert_less(n2, nn); + libmesh_assert_less(n3, nn); + libmesh_assert_less(n4, nn); + + const Point v12 = this->point(n2) - this->point(n1); + const Point v13 = this->point(n3) - this->point(n1); + const Point v14 = this->point(n4) - this->point(n1); + const Real six_vol = triple_product(v12, v13, v14); + libmesh_assert_greater(six_vol, Real(0)); +#endif + + this->_triangulation.push_back({n1, n2, n3, n4}); +} + + +} // namespace libMesh diff --git a/src/geom/cell_polyhedron.C b/src/geom/cell_polyhedron.C new file mode 100644 index 00000000000..d63c51a9fa7 --- /dev/null +++ b/src/geom/cell_polyhedron.C @@ -0,0 +1,744 @@ +// The libMesh Finite Element Library. +// Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner + +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. + +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. + +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +#include "libmesh/cell_polyhedron.h" + +// Local includes +#include "libmesh/face_polygon.h" +#include "libmesh/enum_elem_quality.h" +#include "libmesh/hashword.h" + +// C++ includes +#include +#include +#include + + +namespace libMesh +{ + +// ------------------------------------------------------------ +// Polyhedron class static member initializations +const int Polyhedron::num_children; + +// ------------------------------------------------------------ +// Polyhedron class member functions + + +Polyhedron::Polyhedron (const std::vector> & sides, + Elem * p) : + Cell(/* unused here */ 0, sides.size(), p, nullptr, nullptr), + _elemlinks_data(sides.size()+2), // neighbors + parent + interior_parent + _nodelinks_data(0), // We'll have to resize *this* later too! + _sidelinks_data(sides.size()) + { + // Set our sides, and while we're at it figure out our node maps + // and side normal directions and edge lookup table, and count our + // sides' nodes. If we have internal nodes too then the subclass + // will append those afterward. + unsigned int nn = 0; + std::unordered_map local_node_number; + std::unordered_set, + libMesh::hash> edges_seen; + std::unique_ptr edge; + for (unsigned int s : index_range(sides)) + { + libmesh_assert(sides[s].get()); + auto & side_tuple = _sidelinks_data[s]; + std::get<0>(side_tuple) = sides[s]; + + Polygon & side = *sides[s]; // not const, for writeable nodes + for (auto n : make_range(side.n_nodes())) + { + Node * node = side.node_ptr(n); + if (auto it = local_node_number.find(node); + it != local_node_number.end()) + { + std::get<2>(side_tuple).push_back(it->second); + } + else + { + std::get<2>(side_tuple).push_back(nn); + local_node_number[node] = nn++; + _nodelinks_data.push_back(node); + } + } + + for (unsigned int e : make_range(side.n_edges())) + { + side.build_edge_ptr(edge, e); + auto edge_vertices = std::make_pair(edge->node_ptr(0), edge->node_ptr(1)); + if (edge_vertices.first > edge_vertices.second) + std::swap(edge_vertices.first, edge_vertices.second); + + if (!edges_seen.count(edge_vertices)) + { + edges_seen.insert(edge_vertices); + _edge_lookup.emplace_back(s, e); + } + } + } + + // Do the manual initialization that Elem::Elem and Cell::Cell + // couldn't, now that we've resized both our vectors. No need to + // manually set nullptr, though, since std::vector does that. + this->_elemlinks = _elemlinks_data.data(); + this->_nodes = _nodelinks_data.data(); + this->_elemlinks[0] = p; + + libmesh_assert_equal_to(nn, this->n_nodes()); + + // Figure out the orientation of our sides, now that we've got our + // nodes organized enough to find our center. The algorithm below + // only works for convex polyhedra, but that's all we're + // supporting for now. + Point center; + for (auto n : make_range(nn)) + center.add (this->point(n)); + center /= static_cast(nn); + + for (unsigned int s : index_range(sides)) + { + const Polygon & side = *sides[s]; + const Point x_i = side.point(0); + const Point n_i = + (side.point(1) - side.point(0)).cross + (side.point(0) - side.point(side.n_sides()-1)).unit(); + + bool & inward_normal = std::get<1>(_sidelinks_data[s]); + inward_normal = (n_i * (center - x_i) > TOLERANCE); + } + + // We're betting a lot on "our polyhedra are all convex", so let's + // check that if we have time. +#ifdef DEBUG + for (unsigned int s : index_range(sides)) + { + const Polygon & side = *sides[s]; + const Point x_i = side.point(0); + const bool inward_normal = std::get<1>(this->_sidelinks_data[s]); + + const Point n_i = + (side.point(1) - side.point(0)).cross + (side.point(0) - side.point(side.n_sides()-1)).unit() * + (inward_normal ? -1 : 1); + + for (const Point & node : this->node_ref_range()) + { + const Point d_n = node - x_i; + if (d_n * n_i > TOLERANCE * d_n.norm()) + libmesh_not_implemented_msg + ("Cannot create a non-convex polyhedron"); + } + } +#endif + + // Is this likely to ever be used? We may do refinement with + // polyhedra but it's probably not going to have a hierarchy... + if (p) + { + this->subdomain_id() = p->subdomain_id(); + this->processor_id() = p->processor_id(); + _map_type = p->mapping_type(); + _map_data = p->mapping_data(); + +#ifdef LIBMESH_ENABLE_AMR + this->set_p_level(p->p_level()); +#endif + } + + // Make sure the interior parent isn't undefined + this->set_interior_parent(nullptr); + } + + + +Point Polyhedron::master_point (const unsigned int i) const +{ + return this->point(i); +} + + + +bool Polyhedron::convex() +{ + for (unsigned int s : make_range(this->n_sides())) + { + const Polygon & side = *std::get<0>(this->_sidelinks_data[s]); + const Point x_i = side.point(0); + const bool inward_normal = std::get<1>(this->_sidelinks_data[s]); + + const Point n_i = + (side.point(1) - side.point(0)).cross + (side.point(0) - side.point(side.n_sides()-1)).unit() * + (inward_normal ? -1 : 1); + + for (const Point & node : this->node_ref_range()) + { + const Point d_n = node - x_i; + if (d_n * n_i > TOLERANCE * d_n.norm()) + return false; + } + } + return true; +} + + + +bool Polyhedron::on_reference_element(const Point & p, + const Real eps) const +{ + const unsigned int ns = this->n_sides(); + + // Check that the point is on the same side of all the faces by + // testing whether: + // + // n_i.(p - x_i) <= 0 + // + // for each i, where: + // n_i is the outward normal of face i, + // x_i is a point on face i. + + for (auto i : make_range(ns)) + { + const Polygon & face = *std::get<0>(this->_sidelinks_data[i]); + const bool inward_normal = std::get<1>(this->_sidelinks_data[i]); + + const Point x_i = face.point(0); + + const Point n_i = + (face.point(1) - face.point(0)).cross + (face.point(0) - face.point(face.n_sides()-1)).unit() * + (inward_normal ? -1 : 1); + + // This only works for polyhedra with flat sides. +#ifdef DEBUG + for (auto j : make_range(face.n_sides()-1)) + { + const Point x_j = face.point(j+1); + const Point d_j = x_j - x_i; + if (std::abs(d_j * n_i) > eps * d_j.norm()) + libmesh_not_implemented_msg + ("Polyhedra with non-flat sides are not fully supported."); + } +#endif + + if (n_i * (p - x_i) > eps) + return false; + } + + return true; +} + + + +dof_id_type Polyhedron::key (const unsigned int s) const +{ + libmesh_assert_less (s, this->n_sides()); + + const Polygon & face = *std::get<0>(this->_sidelinks_data[s]); + + return face.key(); +} + + + +dof_id_type Polyhedron::low_order_key (const unsigned int s) const +{ + libmesh_assert_less (s, this->n_sides()); + + const Polygon & face = *std::get<0>(this->_sidelinks_data[s]); + + const unsigned int nv = face.n_vertices(); + std::vector vertex_ids(nv); + for (unsigned int v : make_range(nv)) + vertex_ids[v] = face.node_id(v); + + return Utility::hashword(vertex_ids); +} + + + +unsigned int Polyhedron::local_side_node(unsigned int side, + unsigned int side_node) const +{ + libmesh_assert_less (side, this->n_sides()); + + const std::vector & node_map = + std::get<2>(this->_sidelinks_data[side]); + libmesh_assert_less (side_node, node_map.size()); + + return node_map[side_node]; +} + + + +unsigned int Polyhedron::local_edge_node(unsigned int edge, + unsigned int edge_node) const +{ + libmesh_assert_less (edge, this->n_edges()); + libmesh_assert_less (edge, _edge_lookup.size()); + + auto [side, edge_of_side] = _edge_lookup[edge]; + + const Polygon & face = *std::get<0>(this->_sidelinks_data[side]); + + const std::vector & node_map = + std::get<2>(this->_sidelinks_data[side]); + + return node_map[face.local_edge_node(edge_of_side, edge_node)]; +} + + + +dof_id_type Polyhedron::key () const +{ + std::vector node_ids; + for (const auto & n : this->node_ref_range()) + node_ids.push_back(n.id()); + + return Utility::hashword(node_ids); +} + + + +std::unique_ptr Polyhedron::side_ptr (const unsigned int i) +{ + libmesh_assert_less (i, this->n_sides()); + + Polygon & face = *std::get<0>(this->_sidelinks_data[i]); + std::unique_ptr face_copy = face.disconnected_clone(); + for (auto n : face.node_index_range()) + face_copy->set_node(n, face.node_ptr(n)); + + return face_copy; +} + + + +void Polyhedron::side_ptr (std::unique_ptr & side, + const unsigned int i) +{ + libmesh_assert_less (i, this->n_sides()); + + // Polyhedra are irregular enough that we're not even going to try + // and bother optimizing heap access here. + side = this->side_ptr(i); +} + + + +std::unique_ptr Polyhedron::build_side_ptr (const unsigned int i, + bool /*proxy*/) +{ + auto returnval = this->side_ptr(i); + returnval->set_interior_parent(this); + returnval->set_mapping_type(this->mapping_type()); + returnval->subdomain_id() = this->subdomain_id(); +#ifdef LIBMESH_ENABLE_AMR + returnval->set_p_level(this->p_level()); +#endif + return returnval; +} + + + +void Polyhedron::build_side_ptr (std::unique_ptr & side, + const unsigned int i) +{ + this->side_ptr(side, i); + side->set_interior_parent(this); + side->set_mapping_type(this->mapping_type()); + side->subdomain_id() = this->subdomain_id(); +#ifdef LIBMESH_ENABLE_AMR + side->set_p_level(this->p_level()); +#endif +} + + + +std::unique_ptr Polyhedron::build_edge_ptr (const unsigned int i) +{ + auto [s, se] = _edge_lookup[i]; + Polygon & face = *std::get<0>(_sidelinks_data[s]); + return face.build_edge_ptr(se); +} + + + +void Polyhedron::build_edge_ptr (std::unique_ptr & elem, + const unsigned int i) +{ + auto [s, se] = _edge_lookup[i]; + Polygon & face = *std::get<0>(_sidelinks_data[s]); + face.build_edge_ptr(elem, se); +} + + + +bool Polyhedron::is_child_on_side(const unsigned int /*c*/, + const unsigned int /*s*/) const +{ + libmesh_not_implemented(); + return false; +} + + + +unsigned int Polyhedron::opposite_side(const unsigned int /* side_in */) const +{ + // This is too ambiguous in general. + libmesh_not_implemented(); + return libMesh::invalid_uint; +} + + + +unsigned int Polyhedron::opposite_node(const unsigned int /* n */, + const unsigned int /* s */) const +{ + // This is too ambiguous in general. + libmesh_not_implemented(); + return libMesh::invalid_uint; +} + + + +bool Polyhedron::is_flipped() const +{ + if (this->_triangulation.empty()) + return false; + + auto & tet = this->_triangulation[0]; + + const Point v01 = this->point(tet[1]) - this->point(tet[0]); + const Point v02 = this->point(tet[2]) - this->point(tet[0]); + const Point v03 = this->point(tet[3]) - this->point(tet[0]); + + return (triple_product(v01, v02, v03) < 0); +} + + +std::vector +Polyhedron::edges_adjacent_to_node(const unsigned int n) const +{ + libmesh_assert_less(n, this->n_nodes()); + + // For mid-edge or mid-face nodes, the subclass had better have + // overridden this. + libmesh_assert_less(n, this->n_vertices()); + + const unsigned int ne = this->n_edges(); + libmesh_assert_equal_to(ne, _edge_lookup.size()); + + std::vector adjacent_edges; + + unsigned int next_edge = 0; + + // Look for any adjacent edges on each side. Make use of the fact + // that we number our edges in order as we encounter them from each + // side. + for (auto t : index_range(_sidelinks_data)) + { + const Polygon & face = *std::get<0>(_sidelinks_data[t]); + const std::vector & node_map = + std::get<2>(_sidelinks_data[t]); + + while (_edge_lookup[next_edge].first < t) + { + ++next_edge; + if (next_edge == ne) + return adjacent_edges; + } + + // If we haven't seen the next edge on this or an earlier side + // then we might as well go on to the next. + if (_edge_lookup[next_edge].first > t) + continue; + + const unsigned int fnv = face.n_vertices(); + libmesh_assert_equal_to(fnv, face.n_edges()); + libmesh_assert_equal_to(fnv, face.n_sides()); + libmesh_assert_less_equal(fnv, node_map.size()); + + // Polygon faces have one edge per vertex + for (auto v : make_range(fnv)) + { + libmesh_assert_equal_to (_edge_lookup[next_edge].first, t); + + if (_edge_lookup[next_edge].second > v) + continue; + + while (_edge_lookup[next_edge].first == t && + _edge_lookup[next_edge].second < v) + { + ++next_edge; + if (next_edge == ne) + return adjacent_edges; + } + + if (_edge_lookup[next_edge].first > t) + break; + + const unsigned int vn = node_map[v]; + const unsigned int vnp = node_map[(v+1)%fnv]; + libmesh_assert_less(vn, this->n_vertices()); + libmesh_assert_less(vnp, this->n_vertices()); + if (vn == n || vnp == n) + adjacent_edges.push_back(next_edge); + } + } + + return adjacent_edges; +} + + +std::pair Polyhedron::qual_bounds (const ElemQuality q) const +{ + std::pair bounds; + + switch (q) + { + case EDGE_LENGTH_RATIO: + bounds.first = 1.; + bounds.second = 4.; + break; + + case MIN_ANGLE: + bounds.first = 30.; + bounds.second = 180.; + break; + + case MAX_ANGLE: + bounds.first = 60.; + bounds.second = 180.; + break; + + case JACOBIAN: + case SCALED_JACOBIAN: + bounds.first = 0.5; + bounds.second = 1.; + break; + + default: + libMesh::out << "Warning: Invalid quality measure chosen." << std::endl; + bounds.first = -1; + bounds.second = -1; + } + + return bounds; +} + + + +std::vector> +Polyhedron::side_clones() const +{ + const auto ns = this->n_sides(); + + libmesh_assert_equal_to(ns, _sidelinks_data.size()); + + std::vector> cloned_sides(ns); + + for (auto i : make_range(ns)) + { + const Polygon & face = *std::get<0>(this->_sidelinks_data[i]); + + Elem * clone = face.disconnected_clone().release(); + Polygon * polygon_clone = cast_ptr(clone); + cloned_sides[i] = std::shared_ptr(polygon_clone); + + // We can't actually use a *disconnected* clone to reconstruct + // links between sides, so we'll temporarily give the clone our + // own nodes; user code that typically replaces the usual + // nullptr with permanent nodes will then instead place our + // nodes with permanent nodes. + for (auto n : make_range(face.n_nodes())) + cloned_sides[i]->set_node + (n, const_cast(face.node_ptr(n))); + } + + return cloned_sides; +} + + + +bool Polyhedron::side_has_edge_nodes(unsigned int s, + unsigned int min_node, + unsigned int max_node) const +{ + const Polygon & face = *std::get<0>(_sidelinks_data[s]); + const std::vector & node_map = + std::get<2>(this->_sidelinks_data[s]); + + for (unsigned int e : make_range(face.n_sides())) + { + std::vector nodes_on_edge = + face.nodes_on_side(e); + libmesh_assert_equal_to(nodes_on_edge.size(), 2); + nodes_on_edge[0] = node_map[nodes_on_edge[0]]; + nodes_on_edge[1] = node_map[nodes_on_edge[1]]; + if ((nodes_on_edge[0] == min_node) && + (nodes_on_edge[1] == max_node)) + return true; + if ((nodes_on_edge[1] == min_node) && + (nodes_on_edge[0] == max_node)) + return true; + } + + return false; +} + + + +std::vector +Polyhedron::sides_on_edge(const unsigned int e) const +{ + std::vector returnval(2); + auto [s1, s1e] = _edge_lookup[e]; + returnval[0] = s1; + + const Polygon & face1 = *std::get<0>(_sidelinks_data[s1]); + const std::vector & node_map = + std::get<2>(this->_sidelinks_data[s1]); + + std::vector nodes_on_edge = + face1.nodes_on_side(s1e); + libmesh_assert_equal_to(nodes_on_edge.size(), 2); + nodes_on_edge[0] = node_map[nodes_on_edge[0]]; + nodes_on_edge[1] = node_map[nodes_on_edge[1]]; + + if (nodes_on_edge[0] > nodes_on_edge[1]) + std::swap(nodes_on_edge[0], nodes_on_edge[1]); + + for (unsigned int s2 : make_range(this->n_sides())) + { + if (s2 == s1) + continue; + + if (this->side_has_edge_nodes(s2, nodes_on_edge[0], + nodes_on_edge[1])) + { + returnval[1] = s2; + return returnval; + } + } + + libmesh_error(); + + return returnval; +} + + + +bool Polyhedron::is_edge_on_side(const unsigned int e, + const unsigned int s) const +{ + auto [s1, s1e] = _edge_lookup[e]; + + // Did we get lucky with our cache? + if (s1 == s) + return true; + + const Polygon & face1 = *std::get<0>(_sidelinks_data[s1]); + const std::vector & node_map = + std::get<2>(this->_sidelinks_data[s1]); + std::vector nodes_on_edge1 = + face1.nodes_on_side(s1e); + libmesh_assert_equal_to(nodes_on_edge1.size(), 2); + + nodes_on_edge1[0] = node_map[nodes_on_edge1[0]]; + nodes_on_edge1[1] = node_map[nodes_on_edge1[1]]; + if (nodes_on_edge1[0] > nodes_on_edge1[1]) + std::swap(nodes_on_edge1[0], nodes_on_edge1[1]); + + return this->side_has_edge_nodes(s, + nodes_on_edge1[0], + nodes_on_edge1[1]); +} + + + +std::array Polyhedron::master_subelement (unsigned int i) const +{ + libmesh_assert_less(i, this->_triangulation.size()); + + const auto & tet = this->_triangulation[i]; + + return { this->master_point(tet[0]), + this->master_point(tet[1]), + this->master_point(tet[2]), + this->master_point(tet[3]) }; +} + + +std::tuple +Polyhedron::subelement_coordinates (const Point & p, Real tol) const +{ + std::tuple returnval = + {libMesh::invalid_uint, -1, -1, -1}; + + Real best_bad_coord = -1; + + for (auto s : make_range(this->n_subelements())) + { + const std::array subtet = + this->master_subelement(s); + + // Find barycentric coordinates in subelem + const Point v0 = p - subtet[0]; + // const Point v1 = p - subtet[1]; + + const Point v01 = subtet[1] - subtet[0]; + const Point v02 = subtet[2] - subtet[0]; + const Point v03 = subtet[3] - subtet[0]; + + // const Point v12 = subtet[2] - subtet[1]; + // const Point v13 = subtet[3] - subtet[1]; + + // const Real tp0 = triple_product(v1, v13, v12); + const Real tp1 = triple_product(v0, v02, v03); + const Real tp2 = triple_product(v0, v03, v01); + const Real tp3 = triple_product(v0, v01, v02); + + const Real six_vol = triple_product(v01, v02, v03); + + const Real xi = tp1 / six_vol; + const Real eta = tp2 / six_vol; + const Real zeta = tp3 / six_vol; + + if (xi>=0 && eta>=0 && zeta>=0 && xi+eta+zeta<=1) + return { s, xi, eta, zeta }; + + const Real my_best_bad_coord = + std::min(std::min(std::min(xi, eta), zeta), 1-xi-eta-zeta); + + if (my_best_bad_coord > best_bad_coord) + { + best_bad_coord = my_best_bad_coord; + returnval = { s, xi, eta, zeta }; + } + } + + if (best_bad_coord > -tol) + return returnval; + + return {libMesh::invalid_uint, -1, -1, -1}; +} + + +} // namespace libMesh diff --git a/src/geom/elem.C b/src/geom/elem.C index 046428d70f2..9125a4d54fa 100644 --- a/src/geom/elem.C +++ b/src/geom/elem.C @@ -158,6 +158,7 @@ const unsigned int Elem::type_to_dim_map [] = 2, // QUADSHELL9 2, // C0POLYGON + 3, // C0POLYHEDRON }; const unsigned int Elem::max_n_nodes; @@ -220,6 +221,7 @@ const unsigned int Elem::type_to_n_nodes_map [] = 9, // QUADSHELL9 invalid_uint, // C0POLYGON + invalid_uint, // C0POLYHEDRON }; const unsigned int Elem::type_to_n_sides_map [] = @@ -280,6 +282,7 @@ const unsigned int Elem::type_to_n_sides_map [] = 4, // QUADSHELL9 invalid_uint, // C0POLYGON + invalid_uint, // C0POLYHEDRON }; const unsigned int Elem::type_to_n_edges_map [] = @@ -340,6 +343,7 @@ const unsigned int Elem::type_to_n_edges_map [] = 4, // QUADSHELL9 invalid_uint, // C0POLYGON + invalid_uint, // C0POLYHEDRON }; const Order Elem::type_to_default_order_map [] = @@ -400,6 +404,7 @@ const Order Elem::type_to_default_order_map [] = SECOND, // QUADSHELL9 FIRST, // C0POLYGON + FIRST, // C0POLYHEDRON }; // ------------------------------------------------------------ @@ -484,6 +489,12 @@ std::unique_ptr Elem::build(const ElemType type, case C0POLYGON: return std::make_unique(6, p); + // Building a polyhedron can't currently be done without creating + // its nodes first + case C0POLYHEDRON: + libmesh_not_implemented_msg + ("Polyhedra cannot be built via Elem::build()"); + // 3D elements case TET4: return std::make_unique(p); diff --git a/src/geom/reference_elem.C b/src/geom/reference_elem.C index 25350c01b54..5c6abbf88a0 100644 --- a/src/geom/reference_elem.C +++ b/src/geom/reference_elem.C @@ -212,8 +212,8 @@ void init_ref_elem_table() ref_elem_file[PYRAMID14] = ElemDataStrings::one_pyramid14; ref_elem_file[PYRAMID18] = ElemDataStrings::one_pyramid18; - // No entry here for C0POLYGON - the reference element there - // depends on the number of nodes in the polygon and can't be + // No entry here for C0POLYGON or C0POLYHEDRON - the reference + // element there depends on the number of nodes and can't be // precomputed. } diff --git a/src/libmesh_SOURCES b/src/libmesh_SOURCES index 9b6df52cc16..82a031240f2 100644 --- a/src/libmesh_SOURCES +++ b/src/libmesh_SOURCES @@ -121,6 +121,7 @@ libmesh_SOURCES = \ src/fe/inf_fe_static.C \ src/geom/bounding_box.C \ src/geom/cell.C \ + src/geom/cell_c0polyhedron.C \ src/geom/cell_hex.C \ src/geom/cell_hex20.C \ src/geom/cell_hex27.C \ @@ -133,6 +134,7 @@ libmesh_SOURCES = \ src/geom/cell_inf_prism.C \ src/geom/cell_inf_prism12.C \ src/geom/cell_inf_prism6.C \ + src/geom/cell_polyhedron.C \ src/geom/cell_prism.C \ src/geom/cell_prism15.C \ src/geom/cell_prism18.C \ diff --git a/src/mesh/checkpoint_io.C b/src/mesh/checkpoint_io.C index 461673eb125..8534adee78f 100644 --- a/src/mesh/checkpoint_io.C +++ b/src/mesh/checkpoint_io.C @@ -1178,7 +1178,7 @@ void CheckpointIO::read_connectivity (Xdr & io) unsigned int n_nodes = Elem::type_to_n_nodes_map[elem_data[1]]; if (n_nodes == invalid_uint) - libmesh_not_implemented_msg("Support for C0POLYGON not yet implemented"); + libmesh_not_implemented_msg("Support for Polygons/Polyhedra not yet implemented"); // Snag the node ids this element was connected to std::vector conn_data(n_nodes); diff --git a/src/mesh/dyna_io.C b/src/mesh/dyna_io.C index 4992cb4cf97..a93368027f6 100644 --- a/src/mesh/dyna_io.C +++ b/src/mesh/dyna_io.C @@ -108,7 +108,7 @@ DynaIO::ElementDefinition::ElementDefinition { const unsigned int n_nodes = Elem::type_to_n_nodes_map[type_in]; if (n_nodes == invalid_uint) - libmesh_not_implemented_msg("Support for C0POLYGON not yet implemented"); + libmesh_not_implemented_msg("Support for Polygons/Polyhedra not yet implemented"); nodes.resize(n_nodes); std::iota(nodes.begin(), nodes.end(), 0); } diff --git a/src/mesh/exodusII_io_helper.C b/src/mesh/exodusII_io_helper.C index 5de966a5c93..83775c9a366 100644 --- a/src/mesh/exodusII_io_helper.C +++ b/src/mesh/exodusII_io_helper.C @@ -2604,7 +2604,7 @@ void ExodusII_IO_Helper::write_elements(const MeshBase & mesh, bool use_disconti const auto & conv = get_conversion(elem_t); num_nodes_per_elem = Elem::type_to_n_nodes_map[elem_t]; if (Elem::type_to_n_nodes_map[elem_t] == invalid_uint) - libmesh_not_implemented_msg("Support for C0POLYGON not yet implemented"); + libmesh_not_implemented_msg("Support for Polygons/Polyhedra not yet implemented"); elem_blk_id.push_back(subdomain_id); elem_type_table.push_back_entry(conv.exodus_elem_type().c_str()); @@ -2834,7 +2834,7 @@ void ExodusII_IO_Helper::write_elements(const MeshBase & mesh, bool use_disconti const auto & conv = get_conversion(elem_t); num_nodes_per_elem = Elem::type_to_n_nodes_map[elem_t]; if (Elem::type_to_n_nodes_map[elem_t] == invalid_uint) - libmesh_not_implemented_msg("Support for C0POLYGON not yet implemented"); + libmesh_not_implemented_msg("Support for Polygons/Polyhedra not yet implemented"); // If this is a *real* block, we just loop over vectors of // element ids to add. diff --git a/src/parallel/parallel_elem.C b/src/parallel/parallel_elem.C index ca11121ccd7..db1d5d80d10 100644 --- a/src/parallel/parallel_elem.C +++ b/src/parallel/parallel_elem.C @@ -77,7 +77,7 @@ Packing::packed_size (std::vector::const_iterator Elem::type_to_n_nodes_map[type]; if (n_nodes == invalid_uint) - libmesh_not_implemented_msg("Support for C0POLYGON not yet implemented"); + libmesh_not_implemented_msg("Support for Polygons/Polyhedra not yet implemented"); const unsigned int n_sides = Elem::type_to_n_sides_map[type]; diff --git a/src/quadrature/quadrature.C b/src/quadrature/quadrature.C index 9581f4a6ca5..e3cc9286174 100644 --- a/src/quadrature/quadrature.C +++ b/src/quadrature/quadrature.C @@ -124,7 +124,7 @@ void QBase::init(const ElemType t, { // Some element types require data from a specific element, so can // only be used with newer APIs. - if (t == C0POLYGON) + if (t == C0POLYGON || t == C0POLYHEDRON) libmesh_error_msg("Code (see stack trace) used an outdated quadrature function overload.\n" "Quadrature rules on a C0Polygon are not defined by its ElemType alone."); diff --git a/src/quadrature/quadrature_gauss_2D.C b/src/quadrature/quadrature_gauss_2D.C index 0fe6e131333..ef7520fee7e 100644 --- a/src/quadrature/quadrature_gauss_2D.C +++ b/src/quadrature/quadrature_gauss_2D.C @@ -1322,6 +1322,9 @@ void QGauss::init_2D() { auto master_points = poly.master_subtriangle(t); + // The factor of one half from the triweights cancels out + // the factor of two here, so we don't need to do so + // ourselves. const Real twice_master_tri_area = (- master_points[1](1) * master_points[2](0) - master_points[0](1) * master_points[1](0) @@ -1330,20 +1333,15 @@ void QGauss::init_2D() - master_points[0](0) * master_points[2](1) + master_points[1](0) * master_points[2](1)); + const Point v01 = master_points[1] - master_points[0]; + const Point v02 = master_points[2] - master_points[0]; + for (std::size_t i = 0; i != numtripts; ++i) { - _points[numtripts*t+i](0) = - master_points[0](0) + - (master_points[1](0) - - master_points[0](0)) * tripoints[i](0) + - (master_points[2](0) - - master_points[0](0)) * tripoints[i](1); - _points[numtripts*t+i](1) = - master_points[0](1) + - (master_points[1](1) - - master_points[0](1)) * tripoints[i](0) + - (master_points[2](1) - - master_points[0](1)) * tripoints[i](1); + _points[numtripts*t+i] = + master_points[0] + + v01 * tripoints[i](0) + + v02 * tripoints[i](1); _weights[numtripts*t+i] = triweights[i] * twice_master_tri_area; } diff --git a/src/quadrature/quadrature_gauss_3D.C b/src/quadrature/quadrature_gauss_3D.C index 7fed79d85cd..13abb547e19 100644 --- a/src/quadrature/quadrature_gauss_3D.C +++ b/src/quadrature/quadrature_gauss_3D.C @@ -22,6 +22,7 @@ #include "libmesh/quadrature_conical.h" #include "libmesh/quadrature_gm.h" #include "libmesh/enum_to_string.h" +#include "libmesh/cell_c0polyhedron.h" namespace libMesh { @@ -550,6 +551,57 @@ void QGauss::init_3D() } + //--------------------------------------------- + // Arbitrary polyhedron quadrature rules + case C0POLYHEDRON: + { + QGauss tet_rule(3, _order); + tet_rule.init(TET4, _p_level, true); + + std::vector & tetpoints = tet_rule.get_points(); + std::vector & tetweights = tet_rule.get_weights(); + + std::size_t numtetpts = tetpoints.size(); + + // C0Polyhedron requires the newer Quadrature API + if (!_elem) + libmesh_error(); + + libmesh_assert(_elem->type() == C0POLYHEDRON); + + const C0Polyhedron & poly = *cast_ptr(_elem); + + std::size_t numtets = poly.n_subelements(); + _points.resize(numtetpts*numtets); + _weights.resize(numtetpts*numtets); + + for (std::size_t t = 0; t != numtets; ++t) + { + auto master_points = poly.master_subelement(t); + + const Point v01 = master_points[1] - master_points[0]; + const Point v02 = master_points[2] - master_points[0]; + const Point v03 = master_points[3] - master_points[0]; + + // The factor of one sixth from the tetweights cancels out + // the factor of six here, so we don't need to do so + // ourselves. + const Real six_master_tet_vol = + triple_product(v01, v02, v03); + + for (std::size_t i = 0; i != numtetpts; ++i) + { + _points[numtetpts*t+i] = + master_points[0] + + v01 * tetpoints[i](0) + + v02 * tetpoints[i](1) + + v03 * tetpoints[i](2); + _weights[numtetpts*t+i] = tetweights[i] * + six_master_tet_vol; + } + } + return; + } //--------------------------------------------- // Unsupported type diff --git a/src/quadrature/quadrature_nodal_3D.C b/src/quadrature/quadrature_nodal_3D.C index 6294c417d91..3d1e5994dbe 100644 --- a/src/quadrature/quadrature_nodal_3D.C +++ b/src/quadrature/quadrature_nodal_3D.C @@ -22,6 +22,7 @@ #include "libmesh/quadrature_trap.h" #include "libmesh/quadrature_simpson.h" #include "libmesh/enum_to_string.h" +#include "libmesh/cell_c0polyhedron.h" namespace libMesh { @@ -302,6 +303,37 @@ void QNodal::init_3D() return; } + //--------------------------------------------- + // Arbitrary polygon quadrature rules + case C0POLYHEDRON: + { + // Polyhedra require the newer Quadrature API + if (!_elem) + libmesh_error(); + + libmesh_assert(_elem->type() == C0POLYHEDRON); + + const C0Polyhedron & poly = *cast_ptr(_elem); + + const Real master_poly_volume = _elem->volume(); + + // For polyhedra the master element is the physical element. + // + // Using even weights is not the most effective nodal rule for + // that case, but users who want an effective rule probably + // want a non-nodal rule anyway. + const unsigned int nn = poly.n_nodes(); + const Real weight = master_poly_volume / nn; + _weights.resize(nn, weight); + + _points.resize(poly.n_nodes()); + for (auto n : make_range(nn)) + _points[n] = poly.master_point(n); + + return; + } + + default: libmesh_error_msg("ERROR: Unsupported type: " << Utility::enum_to_string(_type)); } diff --git a/src/utils/string_to_enum.C b/src/utils/string_to_enum.C index 060e7c0e05e..4d18ac69793 100644 --- a/src/utils/string_to_enum.C +++ b/src/utils/string_to_enum.C @@ -102,6 +102,8 @@ std::map elem_type_to_enum { {"C0POLYGON" , C0POLYGON}, + {"C0POLYHEDRON" , C0POLYHEDRON}, + {"TET" , TET4}, {"TET4" , TET4}, {"TET10" , TET10}, diff --git a/tests/fe/fe_lagrange_test.C b/tests/fe/fe_lagrange_test.C index 089297c550f..62323d1a91d 100644 --- a/tests/fe/fe_lagrange_test.C +++ b/tests/fe/fe_lagrange_test.C @@ -41,4 +41,6 @@ INSTANTIATE_FETEST(FIRST, LAGRANGE, PYRAMID5); INSTANTIATE_FETEST(SECOND, LAGRANGE, PYRAMID13); INSTANTIATE_FETEST(SECOND, LAGRANGE, PYRAMID14); INSTANTIATE_FETEST(THIRD, LAGRANGE, PYRAMID18); + +INSTANTIATE_FETEST(FIRST, LAGRANGE, C0POLYHEDRON); #endif diff --git a/tests/fe/fe_monomial_test.C b/tests/fe/fe_monomial_test.C index 1b5129b4833..b5edd46c63e 100644 --- a/tests/fe/fe_monomial_test.C +++ b/tests/fe/fe_monomial_test.C @@ -64,4 +64,6 @@ INSTANTIATE_FETEST(FOURTH, MONOMIAL, PRISM21); INSTANTIATE_FETEST(FIFTH, MONOMIAL, PRISM15); INSTANTIATE_FETEST(FIFTH, MONOMIAL, PRISM18); INSTANTIATE_FETEST(FIFTH, MONOMIAL, PRISM21); + +INSTANTIATE_FETEST(FIRST, MONOMIAL, C0POLYHEDRON); #endif diff --git a/tests/fe/fe_test.h b/tests/fe/fe_test.h index 646e3b5351b..bb91eca18cb 100644 --- a/tests/fe/fe_test.h +++ b/tests/fe/fe_test.h @@ -3,6 +3,7 @@ #include "test_comm.h" +#include #include #include #include @@ -394,6 +395,47 @@ class FETestBase : public CppUnit::TestCase { _mesh->add_elem(std::move(polygon)); _mesh->prepare_for_use(); } + else if (elem_type == C0POLYHEDRON) + { + // Build a plain cube by hand for testing purposes. We could + // upgrade this, but it should be to something non-skewed so a + // MONOMIAL basis will be able to exactly represent the + // polynomials we're projecting. + // + // See elem_test.h for the limitations here. + + _mesh->add_point(Point(0, 0, 0), 0); + _mesh->add_point(Point(1, 0, 0), 1); + _mesh->add_point(Point(1, 1, 0), 2); + _mesh->add_point(Point(0, 1, 0), 3); + _mesh->add_point(Point(0, 0, 1), 4); + _mesh->add_point(Point(1, 0, 1), 5); + _mesh->add_point(Point(1, 1, 1), 6); + _mesh->add_point(Point(0, 1, 1), 7); + + const std::vector> nodes_on_side = + { {0, 1, 2, 3}, // min z + {0, 1, 5, 4}, // min y + {2, 6, 5, 1}, // max x + {2, 3, 7, 6}, // max y + {0, 4, 7, 3}, // min x + {5, 6, 7, 4} }; // max z + + // Build all the sides. + std::vector> sides(nodes_on_side.size()); + + for (auto s : index_range(nodes_on_side)) + { + const auto & nodes_on_s = nodes_on_side[s]; + sides[s] = std::make_shared(nodes_on_s.size()); + for (auto i : index_range(nodes_on_s)) + sides[s]->set_node(i, _mesh->node_ptr(nodes_on_s[i])); + } + + std::unique_ptr polyhedron = std::make_unique(sides); + _mesh->add_elem(std::move(polyhedron)); + _mesh->prepare_for_use(); + } else { MeshTools::Generation::build_cube (*_mesh, diff --git a/tests/fe/fe_xyz_test.C b/tests/fe/fe_xyz_test.C index b7aaec9c6bb..d6ebde8fce7 100644 --- a/tests/fe/fe_xyz_test.C +++ b/tests/fe/fe_xyz_test.C @@ -52,4 +52,7 @@ INSTANTIATE_FETEST(THIRD, XYZ, PRISM21); INSTANTIATE_FETEST(FOURTH, XYZ, PRISM15); INSTANTIATE_FETEST(FOURTH, XYZ, PRISM18); INSTANTIATE_FETEST(FOURTH, XYZ, PRISM21); + +INSTANTIATE_FETEST(FIRST, XYZ, C0POLYHEDRON); +INSTANTIATE_FETEST(SECOND, XYZ, C0POLYHEDRON); #endif diff --git a/tests/geom/elem_test.C b/tests/geom/elem_test.C index 03f55460c62..578ff186835 100644 --- a/tests/geom/elem_test.C +++ b/tests/geom/elem_test.C @@ -137,9 +137,9 @@ public: // test are found in tetrahedra, which have a minimum value of // 2.0. However, we return a default value of 1.0 for 0D and // 1D elements here, and we have a custom distorted-pentagon - // C0Polygon with a 0.5 at 2 nodes, so we handle those cases - // too. - if (elem->dim() < 2) + // C0Polygon with a 0.5 at 2 nodes and a custom C0Polyhedron + // with a 1, so we handle those cases too. + if (elem->dim() < 2 || elem->type() == C0POLYHEDRON) CPPUNIT_ASSERT_GREATEREQUAL(1 - TOLERANCE, jac); else if (elem->type() == C0POLYGON) CPPUNIT_ASSERT_GREATEREQUAL(0.5 - TOLERANCE, jac); @@ -298,6 +298,9 @@ public: if (elem->infinite()) continue; + if (elem->type() == C0POLYHEDRON) + continue; + const Point vertex_avg = elem->vertex_average(); const unsigned int n_sides = elem->n_sides(); @@ -401,11 +404,15 @@ public: // Our map should still be affine. // ... except for stupid singular pyramid maps // ... or the polygons we're deliberately testing non-affine + // ... or the polyhedra we deliberately don't define "affine + // map" for. if ((elem->dim() < 3 || elem->n_vertices() != 5) && - elem_type != C0POLYGON) + elem_type != C0POLYGON && + elem_type != C0POLYHEDRON) CPPUNIT_ASSERT(elem->has_affine_map()); - else if (elem_type == C0POLYGON) + else if (elem_type == C0POLYGON && + elem_type != C0POLYHEDRON) CPPUNIT_ASSERT(!elem->has_affine_map()); // The neighbors and bcids should have flipped back to where @@ -595,7 +602,8 @@ public: elem_type == PYRAMID13 || elem_type == PYRAMID14 || elem_type == PYRAMID18 || - elem_type == C0POLYGON) + elem_type == C0POLYGON || + elem_type == C0POLYHEDRON) return; auto refining_mesh = this->_mesh->clone(); @@ -900,6 +908,9 @@ INSTANTIATE_ELEMTEST(QUADSHELL9); // This just tests with one pentagon, but better than nothing INSTANTIATE_ELEMTEST(C0POLYGON); +// And this is just one truncated cube; also better than nothing +INSTANTIATE_ELEMTEST(C0POLYHEDRON); + #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS INSTANTIATE_ELEMTEST(INFQUAD4); INSTANTIATE_ELEMTEST(INFQUAD6); diff --git a/tests/geom/elem_test.h b/tests/geom/elem_test.h index 743a48787a4..c64bb6b71b6 100644 --- a/tests/geom/elem_test.h +++ b/tests/geom/elem_test.h @@ -5,7 +5,8 @@ #include #include -// Avoiding Elem::build() here because we can't pass an n_sides to it +// Avoiding Elem::build() when we can't pass sides or n_sides to it +#include #include #include "libmesh_cppunit.h" @@ -29,11 +30,14 @@ class PerElemTest : public CppUnit::TestCase _mesh = std::make_unique(*TestCommWorld); - std::unique_ptr test_elem = Elem::build(elem_type); + std::unique_ptr test_elem; + + if (elem_type != C0POLYHEDRON) + test_elem = Elem::build(elem_type); #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS #if LIBMESH_DIM > 1 - if (test_elem->infinite()) + if (test_elem.get() && test_elem->infinite()) { Elem * elem = _mesh->add_elem(std::move(test_elem)); @@ -145,6 +149,111 @@ class PerElemTest : public CppUnit::TestCase _mesh->add_elem(std::move(polygon)); _mesh->prepare_for_use(); } + else if (elem_type == C0POLYHEDRON) + { + // There's even less point in having a build_cube for general + // polyhedra, so again we'll hand-make one for testing and + // we'll plan on making a build_dual_mesh for the future. + +#if 0 + // Try a truncated cube for relatively simple verification. + + // We have some differing orientations on the top sides, to + // test handling of that. + // Lower octagon points, counterclockwise + _mesh->add_point(Point(1/Real(3), 0, 0), 0); + _mesh->add_point(Point(2/Real(3), 0, 0), 1); + _mesh->add_point(Point(1, 1/Real(3), 0), 2); + _mesh->add_point(Point(1, 2/Real(3), 0), 3); + _mesh->add_point(Point(2/Real(3), 1, 0), 4); + _mesh->add_point(Point(1/Real(3), 1, 0), 5); + _mesh->add_point(Point(0, 2/Real(3), 0), 6); + _mesh->add_point(Point(0, 1/Real(3), 0), 7); + + // Lower-middle points, counterclockwise + _mesh->add_point(Point(0, 0, 1/Real(3)), 8); + _mesh->add_point(Point(1, 0, 1/Real(3)), 9); + _mesh->add_point(Point(1, 1, 1/Real(3)), 10); + _mesh->add_point(Point(0, 1, 1/Real(3)), 11); + + // Upper-middle points, counterclockwise + _mesh->add_point(Point(0, 0, 2/Real(3)), 12); + _mesh->add_point(Point(1, 0, 2/Real(3)), 13); + _mesh->add_point(Point(1, 1, 2/Real(3)), 14); + _mesh->add_point(Point(0, 1, 2/Real(3)), 15); + + // Upper octagon points, counterclockwise + _mesh->add_point(Point(1/Real(3), 0, 1), 16); + _mesh->add_point(Point(2/Real(3), 0, 1), 17); + _mesh->add_point(Point(1, 1/Real(3), 1), 18); + _mesh->add_point(Point(1, 2/Real(3), 1), 19); + _mesh->add_point(Point(2/Real(3), 1, 1), 20); + _mesh->add_point(Point(1/Real(3), 1, 1), 21); + _mesh->add_point(Point(0, 2/Real(3), 1), 22); + _mesh->add_point(Point(0, 1/Real(3), 1), 23); + + const std::vector> nodes_on_side = + { {0, 1, 2, 3, 4, 5, 6, 7}, // min z + {0, 1, 9, 13, 17, 16, 12, 8}, // min y + {2, 3, 10, 14, 19, 18, 13, 9}, // max x + {4, 5, 11, 15, 21, 20, 14, 10}, // max y + {6, 7, 8, 12, 23, 22, 15, 11}, // min x + {16, 17, 18, 19, 20, 21, 22, 23}, // max z + {7, 0, 8}, // max nothing + {1, 2, 9}, // max x + {3, 4, 10}, // max xy + {5, 6, 11}, // max y + {23, 16, 12}, // max z + {17, 18, 13}, // max xz + {19, 20, 14}, // max xyz + {21, 22, 15} }; // max yz +#endif + +#if 1 + // Or try a plain cube, for simpler debugging + _mesh->add_point(Point(0, 0, 0), 0); + _mesh->add_point(Point(1, 0, 0), 1); + _mesh->add_point(Point(1, 1, 0), 2); + _mesh->add_point(Point(0, 1, 0), 3); + _mesh->add_point(Point(0, 0, 1), 4); + _mesh->add_point(Point(1, 0, 1), 5); + _mesh->add_point(Point(1, 1, 1), 6); + _mesh->add_point(Point(0, 1, 1), 7); + + // With some combinations of face triangulations, even a + // simple cube has no tetrahedralization that doesn't have + // either interior discontinuities or a 0-volume pancake tet! + // + // The initial "natural" way to orient our sides is commented + // out here; permutations that fix diagonals for us are used + // instead. + const std::vector> nodes_on_side = + { {0, 1, 2, 3}, // min z + {0, 1, 5, 4}, // min y + // {1, 2, 6, 5}, // max x - bad + {2, 6, 5, 1}, // max x + {2, 3, 7, 6}, // max y + // {3, 0, 4, 7}, // min x - bad + {0, 4, 7, 3}, // min x + // {4, 5, 6, 7} }; // max z - bad + {5, 6, 7, 4} }; // max z +#endif + + // Build all the sides. + std::vector> sides(nodes_on_side.size()); + + for (auto s : index_range(nodes_on_side)) + { + const auto & nodes_on_s = nodes_on_side[s]; + sides[s] = std::make_shared(nodes_on_s.size()); + for (auto i : index_range(nodes_on_s)) + sides[s]->set_node(i, _mesh->node_ptr(nodes_on_s[i])); + } + + std::unique_ptr polyhedron = std::make_unique(sides); + _mesh->add_elem(std::move(polyhedron)); + _mesh->prepare_for_use(); + } else { const unsigned int dim = test_elem->dim(); diff --git a/tests/geom/volume_test.C b/tests/geom/volume_test.C index a184b95e9b6..e1cf0a28b76 100644 --- a/tests/geom/volume_test.C +++ b/tests/geom/volume_test.C @@ -1,4 +1,5 @@ // libmesh includes +#include #include #include #include @@ -48,6 +49,7 @@ public: CPPUNIT_TEST( testC0PolygonQuad ); CPPUNIT_TEST( testC0PolygonPentagon ); CPPUNIT_TEST( testC0PolygonHexagon ); + CPPUNIT_TEST( testC0PolyhedronCube ); CPPUNIT_TEST_SUITE_END(); public: @@ -1057,6 +1059,133 @@ protected: } + // Helper function to factor out common tests + void testC0PolyhedronMethods(MeshBase & mesh) + { + Elem * elem = mesh.query_elem_ptr(0); + bool found_elem = elem; + mesh.comm().max(found_elem); + CPPUNIT_ASSERT(found_elem); + + if (!elem) + return; + + CPPUNIT_ASSERT_EQUAL(elem->type(), C0POLYHEDRON); + + const int V = elem->n_vertices(); + const int E = elem->n_edges(); + const int F = elem->n_faces(); + + CPPUNIT_ASSERT_EQUAL(int(elem->n_nodes()), V); + CPPUNIT_ASSERT_EQUAL(int(elem->n_sides()), F); + + // Euler characteristic for polygons homeomorphic to spheres is 2 + CPPUNIT_ASSERT_EQUAL(V-E+F, 2); + + int FE = 0; + int FV = 0; + + std::vector sides_on_vertex(V, 0); + for (auto s : make_range(F)) + { + auto side = elem->build_side_ptr(s); + CPPUNIT_ASSERT_EQUAL(side->type(), C0POLYGON); + FE += side->n_edges(); + const int SV = side->n_vertices(); + FV += SV; + + auto nos = elem->nodes_on_side(s); + CPPUNIT_ASSERT_EQUAL(nos.size(), std::size_t(SV)); + + for (auto v : make_range(SV)) + { + Node * sidenode = side->node_ptr(v); + const unsigned int n = elem->get_node_index(sidenode); + CPPUNIT_ASSERT(n != invalid_uint); + ++sides_on_vertex[n]; + + // We're just looking at first order so far + CPPUNIT_ASSERT(side->is_vertex(v)); + CPPUNIT_ASSERT(!side->is_edge(v)); + CPPUNIT_ASSERT(elem->is_vertex(n)); + CPPUNIT_ASSERT(!elem->is_edge(n)); + CPPUNIT_ASSERT(!elem->is_face(n)); + CPPUNIT_ASSERT(elem->is_node_on_side(n, s)); + CPPUNIT_ASSERT(std::find(nos.begin(), nos.end(), n) != nos.end()); + } + } + + // We hit every edge from 2 faces + CPPUNIT_ASSERT_EQUAL(E*2, FE); + + // We hit every vertex from at least 3 faces + CPPUNIT_ASSERT_LESSEQUAL(FV, V*3); // V*3 <= FV + for (auto sov : sides_on_vertex) + CPPUNIT_ASSERT_LESSEQUAL(sov, 3); // 3 <= sov + + CPPUNIT_ASSERT(!elem->is_flipped()); + } + + + + void testC0Polyhedron(const std::vector> & sides, + Real expected_volume) + { + Mesh mesh(*TestCommWorld); + + std::unique_ptr polyhedron = std::make_unique(sides); + + polyhedron->set_id() = 0; + Elem * elem = mesh.add_elem(std::move(polyhedron)); + + const Real derived_volume = elem->volume(); + const Real base_volume = elem->Elem::volume(); + LIBMESH_ASSERT_FP_EQUAL(base_volume, derived_volume, TOLERANCE*TOLERANCE); + LIBMESH_ASSERT_FP_EQUAL(derived_volume, expected_volume, TOLERANCE*TOLERANCE); + + this->testC0PolyhedronMethods(mesh); + } + + + + void testC0PolyhedronCube() + { + ReplicatedMesh mesh(*TestCommWorld); + + mesh.add_point(Point(0, 0, 0), 0); + mesh.add_point(Point(1, 0, 0), 1); + mesh.add_point(Point(1, 1, 0), 2); + mesh.add_point(Point(0, 1, 0), 3); + mesh.add_point(Point(0, 0, 1), 4); + mesh.add_point(Point(1, 0, 1), 5); + mesh.add_point(Point(1, 1, 1), 6); + mesh.add_point(Point(0, 1, 1), 7); + + // See notes in elem_test.h + const std::vector> nodes_on_side = + { {0, 1, 2, 3}, // min z + {0, 1, 5, 4}, // min y + {2, 6, 5, 1}, // max x + {2, 3, 7, 6}, // max y + {0, 4, 7, 3}, // min x + {5, 6, 7, 4} }; // max z + + // Build all the sides. + std::vector> sides(nodes_on_side.size()); + + for (auto s : index_range(nodes_on_side)) + { + const auto & nodes_on_s = nodes_on_side[s]; + sides[s] = std::make_shared(nodes_on_s.size()); + for (auto i : index_range(nodes_on_s)) + sides[s]->set_node(i, mesh.node_ptr(nodes_on_s[i])); + } + + testC0Polyhedron(sides, 1); + } + + + // Helper function that builds the specified type of Elem from a // vector of Points and returns the value of has_invertible_map() // for that Elem.