Feature puma patches [8/x] Batch processing of tree ref_coords [3/x] Curved geometries

example/geometry/t8_example_geometries.cxx

@@ -99,18 +99,21 @@ struct t8_geometry_sincos: public t8_geometry
    */
   void
   t8_geom_evaluate (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, const size_t num_coords,
-                    double out_coords[3]) const
+                    double *out_coords) const
   {
-    if (num_coords != 1)
-      SC_ABORT ("Error: Batch computation of geometry not yet supported.");
-    double x = ref_coords[0];
-    if (gtreeid == 1) {
-      /* Translate ref coordinates by +1 in x direction for the second tree. */
-      x += 1;
+    for (size_t i_coord = 0; i_coord < num_coords; ++i_coord) {
+      const int offset_2d = 2 * i_coord;
+      const int offset_3d = 3 * i_coord;
+      double x = ref_coords[offset_2d];
+      if (gtreeid == 1) {
+        /* Translate ref coordinates by +1 in x direction for the second tree. */
+        x += 1;
+      }
+      out_coords[offset_3d] = x;
+      out_coords[offset_3d + 1] = ref_coords[offset_2d + 1];
+      out_coords[offset_3d + 2]
+        = 0.2 * sin (ref_coords[offset_2d] * 2 * M_PI) * cos (ref_coords[offset_2d + 1] * 2 * M_PI);
     }
-    out_coords[0] = x;
-    out_coords[1] = ref_coords[1];
-    out_coords[2] = 0.2 * sin (ref_coords[0] * 2 * M_PI) * cos (ref_coords[1] * 2 * M_PI);
   }
 
   /* Jacobian, not implemented. */
@@ -171,26 +174,25 @@ struct t8_geometry_moebius: public t8_geometry_with_vertices
    */
   void
   t8_geom_evaluate (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, const size_t num_coords,
-                    double out_coords[3]) const
+                    double *out_coords) const
   {
-    if (num_coords != 1)
-      SC_ABORT ("Error: Batch computation of geometry not yet supported.");
     double t;
     double phi;
 
     /* Compute the linear coordinates (in [0,1]^2) of the reference vertex and store in out_coords. */
-    /* No idea why, but indent insert a lot of newlines here */
-    t8_geom_compute_linear_geometry (active_tree_class, active_tree_vertices, ref_coords, 1, out_coords);
-
-    /* At first, we map x from [0,1] to [-.5,.5]
-     * and y to [0, 2*PI] */
-    t = out_coords[0] - .5;
-    phi = out_coords[1] * 2 * M_PI;
-
-    /* We now apply the parametrization for the moebius strip. */
-    out_coords[0] = (1 - t * sin (phi / 2)) * cos (phi);
-    out_coords[1] = (1 - t * sin (phi / 2)) * sin (phi);
-    out_coords[2] = t * cos (phi / 2);
+    t8_geom_compute_linear_geometry (active_tree_class, active_tree_vertices, ref_coords, num_coords, out_coords);
+
+    for (size_t i_coord = 0; i_coord < num_coords; ++i_coord) {
+      const int offset_3d = i_coord * 3;
+      /* At first, we map x from [0,1] to [-.5,.5]
+      * and y to [0, 2*PI] */
+      t = out_coords[offset_3d] - .5;
+      phi = out_coords[offset_3d + 1] * 2 * M_PI;
+      /* We now apply the parametrization for the moebius strip. */
+      out_coords[offset_3d] = (1 - t * sin (phi / 2)) * cos (phi);
+      out_coords[offset_3d + 1] = (1 - t * sin (phi / 2)) * sin (phi);
+      out_coords[offset_3d + 2] = t * cos (phi / 2);
+    }
   }
 
   /* Jacobian, not implemented. */
@@ -237,13 +239,15 @@ struct t8_geometry_cylinder: public t8_geometry
    */
   void
   t8_geom_evaluate (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, const size_t num_coords,
-                    double out_coords[3]) const
+                    double *out_coords) const
   {
-    if (num_coords != 1)
-      SC_ABORT ("Error: Batch computation of geometry not yet supported.");
-    out_coords[0] = cos (ref_coords[0] * 2 * M_PI);
-    out_coords[1] = ref_coords[1];
-    out_coords[2] = sin (ref_coords[0] * 2 * M_PI);
+    for (size_t i_coord = 0; i_coord < num_coords; ++i_coord) {
+      const int offset_3d = i_coord * 3;
+      const int offset_2d = i_coord * 2;
+      out_coords[offset_3d] = cos (ref_coords[offset_2d] * 2 * M_PI);
+      out_coords[offset_3d + 1] = ref_coords[offset_2d + 1];
+      out_coords[offset_3d + 2] = sin (ref_coords[offset_2d] * 2 * M_PI);
+    }
   }
 
   /* Jacobian, not implemented. */
@@ -306,26 +310,27 @@ struct t8_geometry_circle: public t8_geometry_with_vertices
    */
   void
   t8_geom_evaluate (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, const size_t num_coords,
-                    double out_coords[3]) const
+                    double *out_coords) const
   {
-    if (num_coords != 1)
-      SC_ABORT ("Error: Batch computation of geometry not yet supported.");
     double x;
     double y;
 
     /* Compute the linear coordinates (in [0,1]^2) of the reference vertex and store in out_coords. */
 
     /* No idea why, but indent insert a lot of newlines here */
-    t8_geom_compute_linear_geometry (active_tree_class, active_tree_vertices, ref_coords, 1, out_coords);
-
-    /* We now remap the coords to match the square [-1,1]^2 */
-    x = out_coords[0] * 2 - 1;
-    y = out_coords[1] * 2 - 1;
-
-    /* An now we apply the formula that projects the square to the circle. */
-    out_coords[0] = x * sqrt (1 - y * y / 2);
-    out_coords[1] = y * sqrt (1 - x * x / 2);
-    out_coords[2] = 0;
+    t8_geom_compute_linear_geometry (active_tree_class, active_tree_vertices, ref_coords, num_coords, out_coords);
+
+    for (size_t i_coord = 0; i_coord < num_coords; ++i_coord) {
+      const int offset_3d = i_coord * 3;
+      /* We now remap the coords to match the square [-1,1]^2 */
+      x = out_coords[offset_3d] * 2 - 1;
+      y = out_coords[offset_3d + 1] * 2 - 1;
+
+      /* An now we apply the formula that projects the square to the circle. */
+      out_coords[offset_3d] = x * sqrt (1 - y * y / 2);
+      out_coords[offset_3d + 1] = y * sqrt (1 - x * x / 2);
+      out_coords[offset_3d + 2] = 0;
+    }
   }
 
   /* Jacobian, not implemented. */
@@ -377,29 +382,31 @@ struct t8_geometry_moving: public t8_geometry
    */
   void
   t8_geom_evaluate (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, const size_t num_coords,
-                    double out_coords[3]) const
+                    double *out_coords) const
   {
-    if (num_coords != 1)
-      SC_ABORT ("Error: Batch computation of geometry not yet supported.");
-    double x = ref_coords[0] - .5;
-    double y = ref_coords[1] - .5;
-    const double time = *ptime;
-    double radius_sqr = x * x + y * y;
-    double phi = radius_sqr * (time > 2 ? 4 - time : time);
-
-    /* Change gridlines by applying a 4th order polynomial mapping
-     * [0,1]^2 -> [0,1]^2.
-     * And then map this to [-0.5,-0.5]^2 */
-    int sign = x < 0 ? 1 : -1;
-    double rho = 0.5 - time / 10;
-    x = sign * (1 - exp (-fabs (-x) / rho)) / (2 * (1 - exp (-0.5 / rho)));
-    sign = y < 0 ? 1 : -1;
-    y = sign * (1 - exp (-fabs (-y) / rho)) / (2 * (1 - exp (-0.5 / rho)));
-
-    /* Rotate the x-y axis and add sincos in z axis. */
-    out_coords[0] = x * (cos (phi)) - y * sin (phi);
-    out_coords[1] = y * (cos (phi)) + x * sin (phi);
-    out_coords[2] = 0;
+    for (size_t i_coord = 0; i_coord < num_coords; ++i_coord) {
+      const int offset_3d = i_coord * 3;
+      const int offset_2d = i_coord * 2;
+      double x = ref_coords[offset_2d] - .5;
+      double y = ref_coords[offset_2d + 1] - .5;
+      const double time = *ptime;
+      double radius_sqr = x * x + y * y;
+      double phi = radius_sqr * (time > 2 ? 4 - time : time);
+
+      /* Change gridlines by applying a 4th order polynomial mapping
+      * [0,1]^2 -> [0,1]^2.
+      * And then map this to [-0.5,-0.5]^2 */
+      int sign = x < 0 ? 1 : -1;
+      double rho = 0.5 - time / 10;
+      x = sign * (1 - exp (-fabs (-x) / rho)) / (2 * (1 - exp (-0.5 / rho)));
+      sign = y < 0 ? 1 : -1;
+      y = sign * (1 - exp (-fabs (-y) / rho)) / (2 * (1 - exp (-0.5 / rho)));
+
+      /* Rotate the x-y axis and add sincos in z axis. */
+      out_coords[offset_3d] = x * (cos (phi)) - y * sin (phi);
+      out_coords[offset_3d + 1] = y * (cos (phi)) + x * sin (phi);
+      out_coords[offset_3d + 2] = 0;
+    }
   }
 
   /* Jacobian, not implemented. */
@@ -460,13 +467,16 @@ struct t8_geometry_cube_zdistorted: public t8_geometry
    */
   void
   t8_geom_evaluate (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, const size_t num_coords,
-                    double out_coords[3]) const
+                    double *out_coords) const
   {
-    if (num_coords != 1)
-      SC_ABORT ("Error: Batch computation of geometry not yet supported.");
-    out_coords[0] = ref_coords[0];
-    out_coords[1] = ref_coords[1];
-    out_coords[2] = ref_coords[2] * (0.8 + 0.2 * sin (ref_coords[0] * 2 * M_PI) * cos (ref_coords[1] * 2 * M_PI));
+    for (size_t i_coord = 0; i_coord < num_coords; ++i_coord) {
+      const int offset_3d = i_coord * 3;
+      out_coords[offset_3d] = ref_coords[offset_3d];
+      out_coords[offset_3d + 1] = ref_coords[offset_3d + 1];
+      out_coords[offset_3d + 2]
+        = ref_coords[offset_3d + 2]
+          * (0.8 + 0.2 * sin (ref_coords[offset_3d] * 2 * M_PI) * cos (ref_coords[offset_3d + 1] * 2 * M_PI));
+    }
   }
 
   /* Jacobian, not implemented. */
@@ -651,7 +661,7 @@ t8_analytic_geom (int level, t8_example_geom_type geom_type)
 #if T8_WITH_OCC
     t8_global_productionf ("Creating uniform level %i forests with an cad triangle geometry.\n", level);
 
-    /* Constructing a triangle with one curved edge (f2) */
+    /* Constructing a triangle with one curved edge (f1) */
     Handle_Geom_BSplineCurve cad_curve;
     TColgp_Array1OfPnt point_array (1, 3);
     TopoDS_Shape shape;
@@ -668,7 +678,7 @@ t8_analytic_geom (int level, t8_example_geom_type geom_type)
     shape = BRepBuilderAPI_MakeEdge (cad_curve).Edge ();
 
     /* Create a cad geometry. */
-    t8_cmesh_register_geometry<t8_geometry_cad> (cmesh, 3, shape);
+    t8_cmesh_register_geometry<t8_geometry_cad> (cmesh, 2, shape);
 
     /* The arrays indicate which face/edge carries a geometry. 
        * 0 means no geometry and any other number indicates the position of the geometry 

src/t8_forest/t8_forest_cxx.cxx

@@ -396,12 +396,14 @@ t8_forest_element_from_ref_coords_ext (t8_forest_t forest, t8_locidx_t ltreeid,
                                        const double *ref_coords, const size_t num_coords, double *coords_out,
                                        const double *stretch_factors)
 {
-  double tree_ref_coords[3] = { 0 };
   const t8_eclass_t tree_class = t8_forest_get_tree_class (forest, ltreeid);
+  const int tree_dim = t8_eclass_to_dimension[tree_class];
   const t8_eclass_scheme_c *scheme = t8_forest_get_eclass_scheme (forest, tree_class);
   const t8_cmesh_t cmesh = t8_forest_get_cmesh (forest);
   const t8_gloidx_t gtreeid = t8_forest_global_tree_id (forest, ltreeid);
 
+  double *tree_ref_coords = T8_ALLOC (double, (tree_dim == 0 ? 1 : tree_dim) * num_coords);
+
   if (stretch_factors != NULL) {
 #if T8_ENABLE_DEBUG
     const t8_geometry_type_t geom_type = t8_geometry_get_type (cmesh, gtreeid);
@@ -415,13 +417,15 @@ t8_forest_element_from_ref_coords_ext (t8_forest_t forest, t8_locidx_t ltreeid,
           = 0.5 + ((ref_coords[i_coord * tree_dim + dim] - 0.5) * stretch_factors[dim]);
       }
     }
-
     scheme->t8_element_reference_coords (element, stretched_ref_coords, num_coords, tree_ref_coords);
   }
   else {
     scheme->t8_element_reference_coords (element, ref_coords, num_coords, tree_ref_coords);
   }
+
   t8_geometry_evaluate (cmesh, gtreeid, tree_ref_coords, num_coords, coords_out);
+
+  T8_FREE (tree_ref_coords);
 }
 
 void

src/t8_forest/t8_forest_vtk.cxx

@@ -178,13 +178,14 @@ t8_get_number_of_vtk_nodes (const t8_element_shape_t eclass, const int curved_fl
 #if T8_WITH_VTK
 static void
 t8_forest_vtk_get_element_nodes (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *element, const int vertex,
-                                 double *out_coords)
+                                 const int curved_flag, double *out_coords)
 {
   const t8_eclass_t tree_class = t8_forest_get_tree_class (forest, ltreeid);
   const t8_eclass_scheme_c *scheme = t8_forest_get_eclass_scheme (forest, tree_class);
   const t8_element_shape_t element_shape = scheme->t8_element_shape (element);
   const double *ref_coords = t8_forest_vtk_point_to_element_ref_coords[element_shape][vertex];
-  t8_forest_element_from_ref_coords (forest, ltreeid, element, ref_coords, 1, out_coords);
+  const int num_node = t8_get_number_of_vtk_nodes (element_shape, curved_flag);
+  t8_forest_element_from_ref_coords (forest, ltreeid, element, ref_coords, num_node, out_coords);
 }
 
 /**
@@ -265,16 +266,18 @@ t8_forest_element_to_vtk_cell (
       SC_ABORT_NOT_REACHED ();
     }
   }
+  double *coordinates = T8_ALLOC (double, 3 * num_node);
+  /* Compute coordinates for all vertices inside the domain. */
+  t8_forest_vtk_get_element_nodes (forest, itree, element, 0, curved_flag, coordinates);
   /* For each element we iterate over all points */
-  double coordinates[3];
   for (int ivertex = 0; ivertex < num_node; ivertex++, (*point_id)++) {
-    /* Compute the vertex coordinates inside the domain. */
-    t8_forest_vtk_get_element_nodes (forest, itree, element, ivertex, coordinates);
+    const size_t offset_3d = 3 * ivertex;
     /* Insert point in the points array */
-    points->InsertNextPoint (coordinates[0], coordinates[1], coordinates[2]);
+    points->InsertNextPoint (coordinates[offset_3d], coordinates[offset_3d + 1], coordinates[offset_3d + 2]);
 
     pvtkCell->GetPointIds ()->SetId (ivertex, *point_id);
   }
+  T8_FREE (coordinates);
   /* We insert the next cell in the cell array */
   cellArray->InsertNextCell (pvtkCell);
   /*

src/t8_geometry/t8_geometry_helpers.c

@@ -82,6 +82,10 @@ void
 t8_geom_compute_linear_geometry (t8_eclass_t tree_class, const double *tree_vertices, const double *ref_coords,
                                  const size_t num_coords, double *out_coords)
 {
+  double tri_vertices[9];
+  double line_vertices[6];
+  double base_coords[2];
+  double vec[3];
   int i_dim;
   size_t i_coord;
   const int dimension = t8_eclass_to_dimension[tree_class];
@@ -106,9 +110,7 @@ t8_geom_compute_linear_geometry (t8_eclass_t tree_class, const double *tree_vert
                                         out_coords + offset_domain_dim);
     }
     break;
-  case T8_ECLASS_PRISM: {
-    double tri_vertices[9];
-    double line_vertices[6];
+  case T8_ECLASS_PRISM:
     for (i_coord = 0; i_coord < num_coords; i_coord++) {
       const size_t offset_tree_dim = i_coord * dimension;
       const size_t offset_domain_dim = i_coord * T8_ECLASS_MAX_DIM;
@@ -125,7 +127,7 @@ t8_geom_compute_linear_geometry (t8_eclass_t tree_class, const double *tree_vert
       t8_geom_triangular_interpolation (ref_coords + offset_tree_dim, tri_vertices, T8_ECLASS_MAX_DIM, 2,
                                         out_coords + offset_domain_dim);
     }
-  } break;
+    break;
   case T8_ECLASS_LINE:
   case T8_ECLASS_QUAD:
   case T8_ECLASS_HEX:
@@ -136,14 +138,13 @@ t8_geom_compute_linear_geometry (t8_eclass_t tree_class, const double *tree_vert
                                     out_coords + offset_domain_dim);
     }
     break;
-  case T8_ECLASS_PYRAMID: {
-    double base_coords[2];
-    double vec[3];
+  case T8_ECLASS_PYRAMID:
     for (i_coord = 0; i_coord < num_coords; i_coord++) {
       const size_t offset_tree_dim = i_coord * dimension;
       const size_t offset_domain_dim = i_coord * T8_ECLASS_MAX_DIM;
-      /* Pyramid interpolation. After projecting the point onto the base, we use a bilinear interpolation to do a quad
-       * interpolation on the base and then we interpolate via the height to the top vertex  */
+      /* Pyramid interpolation. After projecting the point onto the base,
+      * we use a bilinear interpolation to do a quad interpolation on the base
+      * and then we interpolate via the height to the top vertex */
 
       /* Project point on base */
       if (ref_coords[offset_tree_dim + 2] != 1.) {
@@ -167,10 +168,10 @@ t8_geom_compute_linear_geometry (t8_eclass_t tree_class, const double *tree_vert
       }
     }
     break;
-  }
   default:
     SC_ABORT ("Linear geometry coordinate computation is only supported for "
               "vertices/lines/triangles/tets/quads/prisms/hexes/pyramids.");
+    break;
   }
 }