F5R.c

Go to the documentation of this file.

#include "F5R.h"
#include "F5B.h"
#include "F5X.h"
#include "F5A.h"
#include "F5private.h"
 
#include <string.h>
#include <assert.h>
 
int     F5Rset_timestep(F5Path*path, long timestep)
{
hid_t   attr_id,
        space_id;
long    have_timestep;
 
        if (!path || path->Grid_hid<0)
        {
                F5printf(1,"F5Rset_timestep(,%ld): got invalid grid", timestep);
                return 0;
        }
        if (F5Rget_timestep(path, &have_timestep) )
        {
                if (have_timestep==timestep)
                        return 1;
 
                F5printf(1,"F5Rset_timestep(%ld): invalid timestep specified, "
                         "grid is already bound to %ld!", timestep, have_timestep);
                return 0;
        }
 
        space_id= H5Screate(H5S_SCALAR);
 
        attr_id = H5Acreate2(path->Grid_hid,
                            FIBER_HDF5_TIMESTEP_ATTRIB, H5T_NATIVE_LONG, space_id,
                            H5P_DEFAULT, H5P_DEFAULT);
 
        H5Awrite(attr_id, H5T_NATIVE_LONG, &timestep);
 
        H5Aclose(attr_id);
        H5Sclose(space_id);
        return 1;
}
 
long*F5Rget_timestep(F5Path*path, long*timestep)
{
hid_t   attr_id,
        space_id;
 
        if (!path || path->Grid_hid<0 || !timestep)
        {
                F5printf(1,"F5Rget_timestep(): got invalid grid or parameter.");
                return NULL;
        }
 
        H5E_BEGIN_TRY
        attr_id = F5Aopen_name(path->Grid_hid, FIBER_HDF5_TIMESTEP_ATTRIB);
        H5E_END_TRY
 
        if (attr_id<0)
                return NULL;
 
        space_id = H5Aget_space(attr_id);
        assert( space_id >0);;
        if (H5Sget_simple_extent_npoints(space_id) != 1)
                return NULL;
 
        F_H5Aread(attr_id, H5T_NATIVE_LONG, timestep, FIBER_HDF5_TIMESTEP_ATTRIB);
 
        H5Aclose(attr_id);
        H5Sclose(space_id);
 
        return timestep;
}
 
 
/**
  A convenience function to create a 3D cartesian chart topology.
  To create topologies in different coordinate systems call F5LTcreate( ... ) directly.
*/
F5Path*F5Rcreate_coordinate_topology(hid_t File_id, const double*time,
                                     const char*gridname,
                                     const char*coordinate_system,
                                     const char*TopologyName,
                                     int IndexDepth, int SkeletonDimensionality,
                                     int Dimensionality,
                                     const hsize_t*refinement)
{
const char*cs = coordinate_system ? coordinate_system : FIBER_HDF5_DEFAULT_CHART;
ChartDomain_IDs*ChartDomain = (strcmp(cs, FIBER_HDF5_DEFAULT_CHART)==0)
                                ? F5B_standard_cartesian_chart3D() : 0;
F5Path*f;
/*
        if (ChartDomain == F5B_standard_cartesian_chart3D() )
        {
                F5printf(45, "F5Rcreate_coordinate_topology(): Using Std chart.");
                //assert(ChartDomain->perm_vector == F5_FORTRAN_ORDER);
        }
        if (!ChartDomain) F5printf(45, "F5Rcreate_coordinate_topology(): No chart domain.");
*/
        f = F5LTcreate(File_id, time,
                       gridname,
                       ChartDomain,
                       F5B_new_global_cartesian_chart3D,
                       cs,
                       TopologyName,
                       IndexDepth,
                       SkeletonDimensionality, Dimensionality,
                       refinement);
 
        if (f)  F5LTmake_enum_type(f);
 
        return f;
}
 
 
F5Path*F5Rcreate_cartesian_nD(hid_t File_id, double time,
                              const char*gridname,
                              int       Dims,
                              const char*coordinate_system)
{
F5Path*myPath;
        F5printf(3, "  ---- F5Rcreate_cartesian(T=%g, grid=`%s', coordinates=`%s')",
                 time, gridname, (coordinate_system?coordinate_system:"Default coordinates"));
 
        F5printf_indent++;
        myPath = F5Rcreate_coordinate_topology(File_id, &time,
                                               gridname,
                                               coordinate_system,
                                               FIBER_HDF5_POINTS,
                                               0, Dims, Dims, 0);
        F5printf_indent--;
        return myPath;
}
 
F5Path*F5Rcreate_cartesian_3D(hid_t File_id, double time,
                             const char*gridname,
                             const char*coordinate_system)
{
        return F5Rcreate_cartesian_nD(File_id, time,
                                      gridname, 3, coordinate_system);
}
 
 
 
F5Path*F5Rcreate_static_cartesian_3D(hid_t File_id,
                                     const char*gridname,
                                     const char*coordinate_system)
{
        F5printf(30, "F5Rcreate_static_cartesian(grid=`%s')\n", gridname);
        return F5Rcreate_coordinate_topology(File_id, NULL,
                                             gridname,
                                             coordinate_system,
                                             FIBER_HDF5_POINTS,
                                             0, 3, 3, 0);
}
 
/*******************************************************************/
F5Path*F5Rcreate_triangles_as_vertices_topology(hid_t File_id, double time, const char*gridname)
{
F5Path*f = F5LTcreate(File_id, &time,
                      gridname,
                      F5B_standard_triangle_cells(), /* ChartDomain */
                      F5B_new_global_triangle_cells, /*F5B_new_global_cartesian_chart3D*/
                      FIBER_HDF5_POINTS,
                      FIBER_HDF5_FACES,
                      1,   /* IndexDepth     */
                      2,   /* Dimensionality */
                      2,
                      0 ); /* refinement     */
 
        if (f)  F5LTmake_enum_type(f);
 
        return f;
}
 
 
F5Path*F5Rcreate_edges_as_vertices_topology(hid_t File_id, double time, const char*gridname)
{
F5Path*f = F5LTcreate(File_id, &time,
                      gridname,
                      F5B_standard_edge_cells(),   /* ChartDomain */
                      F5B_new_global_edge_cells,   /* F5B_new_global_cartesian_chart3D*/
                      FIBER_HDF5_POINTS,
                      FIBER_HDF5_EDGES,
                      1,   /* IndexDepth */
                      1,   /* SkeletonDimensionality */
                      1,   /* Dimensionality */
                      0);  /* Refinement */
 
        if (f)  F5LTmake_enum_type(f);
 
        return f;
}
 
#if 0
F5Path*F5Rcreate_faces_as_vertices_topology(hid_t File_id, double time, const char*gridname)
{
F5Path*f = F5LTcreate(File_id, &time,
                      gridname,
                      F5B_standard_faces_cells(),   /* ChartDomain */
                      F5B_new_global_faces_cells,   /* F5B_new_global_cartesian_chart3D*/
                      FIBER_HDF5_POINTS,
                      FIBER_HDF5_FACES,
                      1,   /* IndexDepth */
                      2,   /* Dimensionality */
                      0);  /* Refinement */
 
        if (f)  F5LTmake_enum_type(f);
 
        return f;
}
#endif
 
F5Path*F5Rcreate_tetrahedrons_as_vertices_topology(hid_t File_id, double time, const char*gridname)
{
F5Path*f = F5LTcreate(File_id, &time,
                      gridname,
                      F5B_standard_tetrahedral_cells(), /* ChartDomain */
                      F5B_new_global_tetrahedral_cells, /*F5B_new_global_cartesian_chart3D*/
                      FIBER_HDF5_POINTS,
                      FIBER_HDF5_CELLS,
                      1,  /* IndexDepth*/
                      3, /* int SkeletonDimensionality */
                      3,  /* Dimensionality*/
                      0); /* refinement*/
 
        if (f)  F5LTmake_enum_type(f);
 
        return f;
}
 
 
F5Path*F5Rcreate_hexaedrons_as_vertices_topology(hid_t File_id, double time, const char*gridname)
{
F5Path*f = F5LTcreate(File_id, &time,
                      gridname,
                      F5B_standard_hexahedral_cells(), /* ChartDomain */
                      F5B_new_global_hexahedral_cells, /*F5B_new_global_cartesian_chart3D*/
                      FIBER_HDF5_POINTS,
                      FIBER_HDF5_CELLS,
                      1, /* IndexDepth */
                      3, /* int SkeletonDimensionality */
                      3, /* Dimensionality */
                      0); /* refinement */
 
        if (f)  F5LTmake_enum_type(f);
 
        return f;
}
 
 
 
 
 
#if 0
F5Path*F5Rcreate_cell_topology3D(hid_t File_id, const double*time,
                                 const char*gridname,
                                 const char*coordinate_system,
                                 const char*TopologyName,
                                 int IndexDepth)
{
char    timename[128];
char    FullTargetTname[2048];
 
F5Path*retval = F5LTcreate(File_id, time,
                            gridname,
                            0, /* Create new refinement chart object... */
                            F5B_new_integer_regular_domain3D,
                            FIBER_HDF5_POINTS,
                            FIBER_HDF5_CELLS,
                            1,
                           3, 3,
                            0);
        assert(0);
 
        F5I_timegroup(time, timename, sizeof(timename) );
 
#ifdef  _MSC_VER
        _snprintf
#else
        snprintf
#endif
                (FullTargetTname, sizeof(FullTargetTname),
                 "/%s/%s/%s",
                 timename, gridname, FIBER_HDF5_POINTS );
 
/* TODO: Check for existing representer, and if another value exists, create unique one! */
        F5Glink(retval->Representation_hid, H5G_LINK_SOFT, FullTargetTname, FIBER_HDF5_REPRESENTER);
 
        return retval;
}
#endif
 
 
/*******************************************************************/
 
/*
  Default topology name generated for topology groups.
*/
void TopologyName(char*name, int nlength,
                  const hsize_t*level, int Centering, int dims)
{
int     i;
char    buf[1024];
int     L[FIBER_MAX_RANK];
 
        assert(level);
 
        if (!Centering)
        {
                strcpy(buf,"VertexLevel");
        }
        else
        {
                strcpy(buf,"CellCentering");
                for(i=0; i<dims; i++)
                {
                        if (Centering & (1<<i))
                        {
                                if (dims<4)
                                        buf[13+i] = i["XYZ"];
                                else
                                        buf[13+i] = '1';
                        }
                        else
                                buf[13+i] = '_';
                }
                buf[13+dims] = '\0';
        }
 
        assert(dims==3);
 
 
        for(i=0; i<dims; i++)
                L[i] = level[i];
 
#ifdef _MSC_VER
        _snprintf
#else
        snprintf
#endif
                (name,nlength,"%s_%dx%dx%d", buf, L[0], L[1], L[2]);
}
 
 
F5Path*F5Rcreate_vertex_refinement3D(hid_t File_id, double time,
                                     const char*gridname,
                                     const hsize_t refinement[3],
                                     const char*coordinate_system)
{
char    Tname[1024];
        assert(refinement);
        TopologyName(Tname, sizeof(Tname), refinement, 0,3);
 
        return F5Rcreate_coordinate_topology(File_id, &time,
                                             gridname,
                                             coordinate_system,
                                             Tname,
                                             0, 3, 3, refinement);
 
}
 
/**@ingroup F5R
 Define refinement by rational numbers.
 Here, the vertices of the refined objects may cover fractions of the
 coarser cells. Still, vertices have to line up exactly and no refined
 cell may be split by the boundary of a coarser cell.
*/
F5Path*F5Rcreate_relative_vertex_Qrefinement3D(hid_t File_id, double time,
                                              const char*gridname,
                                              const hsize_t current_refinement[3],
                                              double target_time,
                                              const hsize_t target_refinement[3])
{
char    CurrentTname[1024];
char    TargetTname[1024];
F5Path* retval;
char    timename[128];
char    FullTargetTname[2048];
        assert(current_refinement);
        assert(target_refinement);
        TopologyName(CurrentTname, sizeof(CurrentTname), current_refinement, 0,3);
        TopologyName(TargetTname , sizeof(TargetTname) , target_refinement, 0,3);
 
        F5printf(30, "F5Rcreate_relative_vertex_Qrefinement3D()\n");
        retval = F5LTcreate(File_id, &time,
                            gridname,
                            0,
                            F5B_new_rational_regular_domain3D,
                            TargetTname, /* representation_name */
                            CurrentTname, /* TopologyName */
                            0,            /* IndexDepth */
                            3,          /* int SkeletonDimensionality, */
                            3,            /* Dimensionality */
                            current_refinement);  /* const hsize_t*refinement */
 
 
        F5printf(30, "F5Rcreate_relative_vertex_Qrefinement3D() - timegroup\n");
        F5I_timegroup(&target_time, timename, sizeof(timename) );
#ifdef  _MSC_VER
        _snprintf
#else
        snprintf
#endif
                (FullTargetTname, sizeof(FullTargetTname),
                 "/%s/%s/%s",
                 timename, gridname, TargetTname );
 
/* TODO: Check for existing representer, and if another value exists, create unique one! */
        F5Glink(retval->Representation_hid, H5G_LINK_SOFT, FullTargetTname, FIBER_HDF5_REPRESENTER);
 
        F5printf(30, "~F5Rcreate_relative_vertex_Qrefinement3D()\n");
 
        return retval;
}
 
/**@ingroup F5R
 Define refinement by integer numbers.
 Here, the vertices of the refined objects have to match up exactly
 with the vertices of the coarse levels. It is only possible to
 refine entire cells on the coarse grid.
 */
F5Path*F5Rcreate_relative_vertex_Irefinement3D(hid_t File_id, double time,
                                              const char*gridname,
                                              const hsize_t current_refinement[3],
                                              double target_time,
                                              const hsize_t target_refinement[3])
{
char    CurrentTname[1024];
char    TargetTname[1024];
F5Path* retval;
char    timename[128];
char    FullTargetTname[2048];
        assert(current_refinement);
        assert(target_refinement);
        TopologyName(CurrentTname, sizeof(CurrentTname), current_refinement, 0,3);
        TopologyName(TargetTname , sizeof(TargetTname) , target_refinement, 0,3);
 
        retval = F5LTcreate(File_id, &time,
                            gridname,
                            0,
                            F5B_new_integer_regular_domain3D,
                            TargetTname, /* representation_name */
                            CurrentTname, /* TopologyName */
                            0,            /* IndexDepth */
                            3,          /* int SkeletonDimensionality */
                            3,            /* Dimensionality */
                            current_refinement);  /* const hsize_t*refinement */
 
 
        F5I_timegroup(&target_time, timename, sizeof(timename) );
#ifdef  _MSC_VER
        _snprintf
#else
        snprintf
#endif
                (FullTargetTname, sizeof(FullTargetTname),
                 "/%s/%s/%s",
                 timename, gridname, TargetTname );
 
/* TODO: Check for existing representer, and if another value exists, create unique one! */
        F5Glink(retval->Representation_hid, H5G_LINK_SOFT, FullTargetTname, FIBER_HDF5_REPRESENTER);
 
        return retval;
}
 
/**\ingroup F5R
   Treat some refinement level as the default level for refinement-unaware
   readers.
 
   Creates a symbolic link from some refinement level on vertices to appear
   as the default points topology.
   @param grid               A valid F5Path object pointing to some grid.
   @param target_refinement  The refinement level that shall be visible as default.
   @note Currently only 3D refinements are supported.
   @note If a point topology already exists, it will not be touched.
 */
int F5Rlink_default_vertex_topology(F5Path*grid,
                                    const hsize_t target_refinement[3])
{
char    TargetTname[1024];
        if (!grid) return -1;
        if (grid->Grid_hid<1) return -2;
        if (!target_refinement) return -3;
        TopologyName(TargetTname, sizeof(TargetTname),  target_refinement, 0,3);
        {
                F5Glink(grid->Grid_hid, H5G_LINK_SOFT, TargetTname, FIBER_HDF5_POINTS);
        }
        return 0;
}
 
 
 
/*******************************************************************/
 
/*
 Image/Viewer/Positions {xyz}
             /Intensity {I}
 
 Image@2x2x2/Image/Positions {222}
                  /Intensity {I}
 
 Image@4x4x4/Image@2x2x2/Positions {222}
            /Image      /Positions {444}
 
 */
static
void    F5Rcreate_sampling(F5Path*Topology, const hsize_t*downsampling,
                           int IndexDepthChange,
                           char NameSeparator)
{
char    oldTopologyName[2048], newTopologyName[2048];
char    *sep;
hid_t   Top_id, space_id, Rep_id, F_id, pos_fieldtype;
hsize_t refinement[FIBER_MAX_RANK], extent[FIBER_MAX_RANK],
        total_refinement[FIBER_MAX_RANK];
int     polynom_order[FIBER_MAX_RANK];
const void*(refinement_coeffs_ptr[FIBER_MAX_RANK]);
hsize_t refinement_coeffs[FIBER_MAX_RANK][2];
int     rank, i;
int     component_map[FIBER_MAX_RANK];
int     IndexDepth = F5LTget_index_depth(Topology->Topology_hid);
        IndexDepth += IndexDepthChange;
 
        H5Iget_name(Topology->Topology_hid, oldTopologyName, sizeof(oldTopologyName) );
        strcpy(newTopologyName, oldTopologyName);
        sep = strrchr(newTopologyName, NameSeparator);
        if (sep) *sep = 0;
        sep = newTopologyName + strlen(newTopologyName);
        *sep++ = '@';
 
        space_id = F5Tget_space(Topology);
        if (space_id < 0)
        {
                puts("void F5Rcreate_sampling(): not a regular topology");
                return;
        }
 
        rank = H5Sget_simple_extent_dims(space_id, extent, NULL);
        H5Sclose(space_id);
 
        F5LAget_dimensions(Topology->Topology_hid, FIBER_HDF5_REFINEMENT_INFO, refinement);
 
        pos_fieldtype = H5Tcreate(H5T_COMPOUND, rank*sizeof(int) );
 
        for(i=0; i<rank; i++)
        {
        char    rank_name[2];
 
                component_map[i]    = i;
                polynom_order[i]    = 1;
                total_refinement[i] = refinement[i]*downsampling[i];
                extent[i]           /= downsampling[i];
 
                refinement_coeffs[i][0]  = 0;
                refinement_coeffs[i][1]  = downsampling[i];
                refinement_coeffs_ptr[i] = refinement_coeffs[i];
 
                rank_name[0] = '0' + i;
                rank_name[1] = 0;
                H5Tinsert(pos_fieldtype, rank_name, i*sizeof(int), H5T_NATIVE_INT);
 
                sprintf(sep, "%u", (unsigned)(total_refinement[i]) );
                sep += strlen(sep);
                if (i<rank-1) *sep++ = 'x';
        }
        *sep = 0;
 
        Top_id = F5LTcreate_topology(Topology->Grid_hid, newTopologyName,
                                     IndexDepth, rank, rank, total_refinement);
 
        Rep_id = F5Gappend(Top_id, oldTopologyName);
        F_id   = F5LwriteX(Rep_id, FIBER_HDF5_POSITIONS_STRING, rank, extent,
                           pos_fieldtype, pos_fieldtype,
                           refinement_coeffs_ptr, polynom_order,
                           component_map,
                           F5P_DEFAULT,
                           Topology->GlobalChart_hid);
 
        H5Tclose(pos_fieldtype);
 
        F5Fclose(Topology);
        H5Gclose(Topology->Representation_hid);
        H5Gclose(Topology->Topology_hid);
        Topology->Topology_hid       = Top_id;
        Topology->Representation_hid = Rep_id;
        Topology->Field_hid          = F_id;
}
 
void    F5Rcreate_downsampling(F5Path*Topology, const hsize_t*downsampling)
{
        F5Rcreate_sampling(Topology, downsampling, +1, '@');
}
 
/*
 Image/Viewer/Positions {xyz}
             /Intensity {I}
 
 Image:2x2x2/Image/Positions {222}
                  /Intensity {I}
 
 Image:4x4x4/Image:2x2x2/Positions {222}
            /Image      /Positions {444}
 
 */
 
void F5Rcreate_subsampling(F5Path*Topology, const hsize_t*subsampling)
{
        F5Rcreate_sampling(Topology, subsampling, 0, ':');
}