Merge pull request #867 from DLR-AMR/fix-curved_geometry

Bugfixes in curved geometry
DLR-AMR · Dec 1, 2023 · 5b4ee94 · 5b4ee94
2 parents dc48cb5 + 5075d22
commit 5b4ee94
Show file tree

Hide file tree

Showing 3 changed files with 31 additions and 1 deletion.
diff --git a/src/t8_cmesh/t8_cmesh_readmshfile.cxx b/src/t8_cmesh/t8_cmesh_readmshfile.cxx
@@ -1250,6 +1250,35 @@ t8_cmesh_msh_file_4_read_eles (t8_cmesh_t cmesh, FILE *fp, sc_hash_t *vertices,
               continue;
             }
 
+            /* If one vertex lies on a geometry of a higher dim as the other, we have to check,
+             * if the geometry of lower dimension is on that geometry. */
+            {
+              int is_on_geom = 1;
+              for (int i_edge = 0; i_edge < 2; ++i_edge) {
+                if (edge_geometry_dim == 2 && edge_nodes[i_edge].entity_dim == 1) {
+                  if (!occ_geometry->t8_geom_is_edge_on_face (edge_nodes[i_edge].entity_tag, edge_geometry_tag)) {
+                    is_on_geom = 0;
+                    break;
+                  }
+                }
+                else if (edge_geometry_dim == 2 && edge_nodes[i_edge].entity_dim == 0) {
+                  if (!occ_geometry->t8_geom_is_vertex_on_face (edge_nodes[i_edge].entity_tag, edge_geometry_tag)) {
+                    is_on_geom = 0;
+                    break;
+                  }
+                }
+                else if (edge_geometry_dim == 1 && edge_nodes[i_edge].entity_dim == 0) {
+                  if (!occ_geometry->t8_geom_is_vertex_on_edge (edge_nodes[i_edge].entity_tag, edge_geometry_tag)) {
+                    is_on_geom = 0;
+                    break;
+                  }
+                }
+              }
+              if (!is_on_geom) {
+                continue;
+              }
+            }
+
             /* If both nodes are on a vertex we still got no edge. 
              * But we can look if both vertices share an edge and use this edge. 
              * If not we can skip this edge. */

diff --git a/src/t8_geometry/t8_geometry_helpers.c b/src/t8_geometry/t8_geometry_helpers.c
@@ -324,6 +324,7 @@ t8_geom_get_ref_intersection (int edge_index, const double *ref_coords, double r
     else if (ref_coords[1] == ref_opposite_vertex[1]) {
       ref_intersection[0] = 0;
       ref_intersection[1] = 1;
+      break;
     }
     else {
       /* intersectionX = (x1y2-y1x2)(x3-x4)-(x1-x2)(x3y4-y3x4)

diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_occ.cxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_occ.cxx
@@ -302,7 +302,7 @@ t8_geometry_occ::t8_geom_evaluate_occ_triangle (t8_cmesh_t cmesh, t8_gloidx_t gt
             t8_geom_linear_interpolation (&ref_intersection[0], parameters, 1, 1, &interpolated_curve_parameter);
           }
           /* Retrieve curve */
-          T8_ASSERT (edges[i_edge] <= occ_shape_face_map.Size ());
+          T8_ASSERT (edges[i_edge] <= occ_shape_edge_map.Size ());
           curve = BRep_Tool::Curve (TopoDS::Edge (occ_shape_edge_map.FindKey (edges[i_edge])), first, last);
           /* Check if curve is valid */
           T8_ASSERT (!curve.IsNull ());