F5Bchart.c

Go to the documentation of this file.

#include "F5Bchart.h"
#include "F5F.h"
#include "F5X.h"
#include "F5defs.h"
#include "F5types.h"
#include "F5A.h"
 
#include "F5private.h"
 
#include "F5coordinates.h"
 
#include "F5T.h"
 
#include <string.h>
 
#ifndef F5_MEM_H
/* #if defined(__APPLE__) && defined(__MACH__) && defined(__POWERPC__) */
#include <stdlib.h>
/* #else */
/* #include <malloc.h> */
/* #endif */
#include <assert.h>
 
#ifdef  _WIN32
#define strdup _strdup
#endif
 
 
#define new(type, size) (type*)calloc(size, sizeof(type))
#define delete(ptr)     if (ptr) free(ptr)
 
#define true    1
#define false   0
#endif
 
const int*F5B_fortran_permute_vector()
{
static  int vals[FIBER_MAX_RANK];
static  int init = 0;
        if (!init)
        {
        int     i;
                for(i=0; i< FIBER_MAX_RANK; i++)
                {
                        vals[i] = FIBER_MAX_RANK - 1 - i;
                }
 
                init = 1;
        }
 
        return vals;
}
 
static  int is_fortran_permutation(
        int perm_vector[FIBER_MAX_RANK],
        int perm_vector_size)
{
int     i;
        if (!perm_vector)
                return 42;
 
        for(i=0; i<perm_vector_size; i++)
        {
                if (perm_vector[i] != -1)
                        return 0;
        }
 
        return 42;
}
 
 
hsize_t*F5B_permute_dimensions(hsize_t*target_dims,
                               int rank,
                               const hsize_t*source_dims,
                               int perm_vector[FIBER_MAX_RANK],
                               int perm_vector_size)
{
hsize_t tmp[FIBER_MAX_RANK];
hsize_t* dst;
 
        if (!source_dims)
        {
        int     i;
                for(i=0; i< rank; i++)
                        target_dims[i] = 0;
                F5printf(10, "F5B_permute_dimensions(): No source dimensions provided!");
                return NULL;
        }
 
        if (perm_vector == F5_C_ORDER)
        {
        int     i;
                F5printf(50, "F5B_permute_dimensions(): C order permutation vector!");
                if (target_dims == source_dims)
                        return target_dims;
 
                for(i=0; i< rank; i++)
                        target_dims[i] = source_dims[i];
 
                return target_dims;
        }
 
        if (target_dims == source_dims)
        {
                if (rank<2)
                {
                        F5printf(40, "F5B_permute_dimensions(): no permutation for 1D!");
                        /* no permutation possible for 1D data... */
                        return target_dims;
                }
                dst = tmp;
        }
        else
                dst = target_dims;
 
/*      if (perm_vector == F5_FORTRAN_ORDER) */
        if (is_fortran_permutation(perm_vector, perm_vector_size) )
        {
        int     i;
                F5printf(40, "F5B_permute_dimensions(): FORTRAN order permutation vector!");
                for(i=0; i< rank; i++)
                        dst[i] = source_dims[rank-1-i];
        }
        else if (rank>1) /* only then it makes sense to copy stuff */
        {
        int     i;
        int     cdims = perm_vector_size,
                   cd = rank - cdims;
 
                assert(perm_vector);
                F5printf(50, "F5B_permute_dimensions(): generic order permutation vector of size %d, data are size %d", cdims, rank);
/*
                for(i=0; i< cdims; i++)
                {
                int     pv = perm_vector[i];
                        F5printf(10, " perm_vector[%d] = %d", i, pv );
                }
                for(i=0; i< rank; i++)
                {
                int     pv = source_dims[i];
                        F5printf(10, " source_dims[%d] = %d", i, pv );
                }
*/
 
                if (rank > cdims)
                {
                        F5printf(50, "F5B_permute_dimensions(): permutation on %d-higher dimensional subspace!", cd);
                        assert( cd == 1 && "higher dimensional subspaces not supported yet!");
 
                        for(i=0; i<cd; i++)
                        {
/*                      int     pv = perm_vector[i+cd]; */
                                /*
                                  Assume FORTRAN convention on subspace dimensions,
                                  but should actually make use of the perm_vector here
                                  as well within the subspace.
                                  If the dimension exceeds the dimension of the subspace,
                                  then things become interesting what to do then...
                                 */
                        int     pv = rank-1-i;
                                assert(  pv < rank );
                                dst[i] = source_dims[ pv ];
                        }
                        for(; i<rank; i++)
                        {
                        int     pv = perm_vector[i-cd];
                                assert(  pv >= 0 );
                                assert(  pv < rank );
                                dst[i] = source_dims[ pv ];
/*                              F5printf(0, " MAP from %d => %d value %d", pv, i, dst[i] ); */
                        }
/*
                        for(i=0; i<rank; i++)
                        {
                        int     pv = dst[i];
                                F5printf(0, " target_dims[%d] = %d", i, pv );
                        }
*/
                }
                else for(i=0; i< rank; i++)
                {
                int     pv = perm_vector[i-cd];
                /*      F5printf(41, " perm_vector[%d] = %d\n", i-cd, pv );  */
                /*      printf(" perm_vector[%d] = %d (rank=%d, cdims=%d, cd=%d)\n", i-cd, pv, rank, cdims, cd );*/
                        assert(  pv >= 0 );
                        assert(  pv < rank );
                        dst[i] = source_dims[ pv ];
                }
        }
        else if (rank==1)
        {
                target_dims[0] = source_dims[0];
                F5printf(50, "~F5B_permute_dimensions() 1D");
                return target_dims;
        }
 
        if (dst != target_dims)
        {
        int     i;
                F5printf(50, "F5B_permute_dimensions(): Copy back permutation vector!");
                for(i=0; i< rank; i++)
                        target_dims[i] = dst[i];
        }
        F5printf(50, "~F5B_permute_dimensions()");
        return target_dims;
}
 
/**\ingroup F5B
  Creates a new chart domain object that holds just a "point" type,
  but no tangential types.
 */
ChartDomain_IDs*F5B_new_global_domain(const char*coordnames[],
                                      int dimension,
                                      const char*name,
                                      const int perm_vector[],
                                      int cell_dimensionality,
                                      hid_t SinglePrecisionType,
                                      hid_t DoublePrecisionType)
{
ChartDomain_IDs *IDs = (ChartDomain_IDs*)malloc(sizeof(ChartDomain_IDs));
int     i;
        memset (IDs, 0, sizeof (ChartDomain_IDs));
 
        IDs->domain_name = strdup(name);
 
        IDs->refs = 1;
        IDs->cell_dimensionality = cell_dimensionality;
 
        F5printf(40, "F5B_new_global_domain(,dims=%d,,name=%s,perm_vector=%p,cell_dims=%d)\n",
                 dimension, name, perm_vector, cell_dimensionality);
 
        if (perm_vector == F5_FORTRAN_ORDER)
        {
                for(i = 0; i<cell_dimensionality; i++)
                {
                        IDs->perm_vector[i] = cell_dimensionality - 1 - i;
                }
        }
        else if (perm_vector == F5_C_ORDER)
        {
                for(i = 0; i<cell_dimensionality; i++)
                {
                        IDs->perm_vector[i] = i;
                }
        }
        else
        {
                for(i = 0; i<cell_dimensionality; i++)
                {
                        IDs->perm_vector[i] = perm_vector[i];
                }
        }
 
        /*
          Set unused dimensions to -1
         */
        if (IDs->perm_vector)
        for(i = cell_dimensionality; i<FIBER_MAX_RANK; i++)
        {
                IDs->perm_vector[i] = -1;
        }
 
        F5B_init_tangential_types(&IDs->SinglePrecision,
                                  coordnames,
                                  dimension,
                                  SinglePrecisionType);
 
        F5B_init_tangential_types(&IDs->DoublePrecision,
                                  coordnames,
                                  dimension,
                                  DoublePrecisionType);
 
        IDs->TensorTypes = NULL;
        /*
        IDs->Point_hid_t = H5Tcreate(H5T_COMPOUND, dimension*Tsize );
        for(i=0; i<dimension; i++)
        {
                H5Tinsert(IDs->Point_hid_t, coordnames[i], i*Tsize, coord_type);
        }
 
        IDs->Vector_hid_t       = -1;
        IDs->Covector_hid_t     = -1;
        IDs->Metric_hid_t       = -1;
        IDs->Bilinearform_hid_t = -1;
        IDs->TensorTypes = NULL;
        */
        return IDs;
}
 
void    F5B_init_tangential_types(F5_ChartPrecisionTypes*IDs,
                                  const char*coordnames[],
                                  int dimension,
                                  hid_t coord_type)
{
size_t  Tsize = H5Tget_size(coord_type);
int     i;
 
        IDs->Point_hid_t = H5Tcreate(H5T_COMPOUND, dimension*Tsize );
        // F5printf(40, "F5B_init_tangential_types, IDs->Point_hid_t = %lli\n", IDs->Point_hid_t);
        for(i=0; i<dimension; i++)
        {
                H5Tinsert(IDs->Point_hid_t, coordnames[i], i*Tsize, coord_type);
        }
        // F5printf(40, "F5B_init_tangential_types, IDs->Point_hid_t = %lli         inserted hid_t %lli\n", IDs->Point_hid_t, coord_type);
        IDs->Vector_hid_t       = 0;
        IDs->Covector_hid_t     = 0;
        IDs->Bivector_hid_t     = 0;
        IDs->Bicovector_hid_t   = 0;
 
        IDs->Pseudoscalar_hid_t   = 0;
        IDs->Pseudocoscalar_hid_t = 0;
 
        IDs->Metric_hid_t       = 0;
        IDs->Cometric_hid_t     = 0;
 
        for(i=0; i<F5ChartPrecisionMaxNumberOfFileTypes; i++)
                IDs->AdditionalFileTypes[ i ] = 0;
}
 
/**\ingroup F5B
   Adds tangential types (vector, covector, metric) to a given
   chart domain. These types are derived from the specified coord_type.
   @todo Get dimension and coordnames from the chart domain's points.
   @todo CALL this function for various numerical precisions
 */
void    F5Bchart_add_tangential_types(F5_ChartPrecisionTypes*IDs,
                                      const char*coordnames[],
                                      int dimension,
                                      hid_t coord_type)
{
size_t  Tsize = H5Tget_size(coord_type);
        IDs->Vector_hid_t = H5Tcreate(H5T_COMPOUND, dimension*Tsize );
        // F5printf(40, "F5B_add_tangential_types, IDs->Vector_hid_t = %lli\n", IDs->Vector_hid_t);
        {int i;
        for(i=0; i<dimension; i++)
        {
        char    buf[_FIBER_MAX_COORDNAME_LENGTH+2] = "D";
                strcat(buf, coordnames[i]);
                H5Tinsert(IDs->Vector_hid_t, buf, i*Tsize, coord_type);
        }
        }
 
        IDs->Covector_hid_t = H5Tcreate(H5T_COMPOUND, dimension*Tsize );
        {int i;
        for(i=0; i<dimension; i++)
        {
        char    buf[_FIBER_MAX_COORDNAME_LENGTH+2] = "d";
                strcat(buf, coordnames[i]);
                H5Tinsert(IDs->Covector_hid_t, buf, i*Tsize, coord_type);
        }
        }
#if 0
        /* needs revision: do permutation without repetition -> ( 2
                                                                  D )   */
        IDs->Bivector_hid_t = H5Tcreate(H5T_COMPOUND, dimension*dimension*Tsize );
        {int i, j;
        for(i=0; i<dimension; i++)
        for(j=i; j<dimension; j++)
        {
        int     L = i*dimension + j;
        char    buf[_FIBER_MAX_COORDNAME_LENGTH+2] = "D";
                strcat(buf, coordnames[i]);
                strcat(buf, coordnames[j]);
 
                H5Tinsert(IDs->Bivector_hid_t, buf, L*Tsize, coord_type);
        }
        }
#endif
        if (dimension == 3)
        {
        IDs->Bivector_hid_t = H5Tcreate(H5T_COMPOUND, dimension*Tsize );
 
        {
        char    buf[_FIBER_MAX_COORDNAME_LENGTH+2] = "D";
                strcat(buf, coordnames[1]);
                strcat(buf, "^D");
                strcat(buf, coordnames[2]);
 
                H5Tinsert(IDs->Bivector_hid_t, buf, 0*Tsize, coord_type);
        }
 
        {
        char    buf[_FIBER_MAX_COORDNAME_LENGTH+2] = "D";
                strcat(buf, coordnames[2]);
                strcat(buf, "^D");
                strcat(buf, coordnames[0]);
 
                H5Tinsert(IDs->Bivector_hid_t, buf, 1*Tsize, coord_type);
        }
 
        {
        char    buf[_FIBER_MAX_COORDNAME_LENGTH+2] = "D";
                strcat(buf, coordnames[0]);
                strcat(buf, "^D");
                strcat(buf, coordnames[1]);
 
                H5Tinsert(IDs->Bivector_hid_t, buf, 2*Tsize, coord_type);
        }
 
        }
        IDs->Bicovector_hid_t = H5Tcreate(H5T_COMPOUND, dimension*dimension*Tsize );
        {int i, j;
        for(i=0; i<dimension; i++)
        for(j=i; j<dimension; j++)
        {
        int     L = i*dimension + j;
        char    buf[_FIBER_MAX_COORDNAME_LENGTH+2] = "d";
                strcat(buf, coordnames[i]);
                strcat(buf, coordnames[j]);
 
                H5Tinsert(IDs->Bicovector_hid_t, buf, L*Tsize, coord_type);
        }
        }
 
        IDs->Metric_hid_t = H5Tcreate(H5T_COMPOUND, dimension*(dimension+1)/2*Tsize );
        {int i, j, ni, nj;
        int     indices[FIBER_MAX_RANK*FIBER_MAX_RANK];
        int     I = 0; /* the linear index */
 
                for(ni=0, i=0; i<dimension; i++ )
                {
                        for(nj=0, j=0; j<i; j++ )
                        {
                                indices[ni+j] = indices[nj+i] = I++;
                                nj += dimension;
                        }
                        indices[ni+i] = I++;
                        ni += dimension;
                }
 
 
        for(i=0; i<dimension; i++)
        for(j=i; j<dimension; j++)
        {
        int     L = i*dimension + j,
                idx = indices[ L ];
        char    buf[_FIBER_MAX_COORDNAME_LENGTH+2] = "g";
                strcat(buf, coordnames[i]);
                strcat(buf, coordnames[j]);
 
                H5Tinsert(IDs->Metric_hid_t, buf, idx*Tsize, coord_type);
        }
        }
 
        IDs->Cometric_hid_t = H5Tcreate(H5T_COMPOUND, dimension*(dimension+1)/2*Tsize );
        {int i, j, ni, nj;
        int     indices[FIBER_MAX_RANK*FIBER_MAX_RANK];
        int     I = 0; /* the linear index */
 
                for(ni=0, i=0; i<dimension; i++ )
                {
                        for(nj=0, j=0; j<i; j++ )
                        {
                                indices[ni+j] = indices[nj+i] = I++;
                                nj += dimension;
                        }
                        indices[ni+i] = I++;
                        ni += dimension;
                }
 
 
        for(i=0; i<dimension; i++)
        for(j=i; j<dimension; j++)
        {
        int     L = i*dimension + j,
                idx = indices[ L ];
        char    buf[_FIBER_MAX_COORDNAME_LENGTH+2] = "G";
                strcat(buf, coordnames[i]);
                strcat(buf, coordnames[j]);
 
                H5Tinsert(IDs->Cometric_hid_t, buf, idx*Tsize, coord_type);
        }
        }
}
 
 
int     F5B_inject_precision_type(F5_ChartPrecisionTypes*IDs,
                                  hid_t type_id)
{
int     i;
 
        for(i=0; i<F5ChartPrecisionMaxNumberOfFileTypes; i++)
        {
                if (!IDs->AdditionalFileTypes[ i ])
                {
                        IDs->AdditionalFileTypes[ i ] = type_id;
                        return 1;
                }
        }
        return 0;
}
 
 
 
 
ChartDomain_IDs*F5B_new_global_float_chart(const char*coordnames[],
                                           int dimension,
                                           const char*name,
                                           const int perm_vector[FIBER_MAX_RANK])
{
ChartDomain_IDs*IDs = F5B_new_global_domain(coordnames,
                                            dimension,
                                            name,
                                            perm_vector,
                                            dimension,
                                            H5T_NATIVE_FLOAT, H5T_NATIVE_DOUBLE);
 
        F5Bchart_add_tangential_types(&IDs->SinglePrecision, coordnames, dimension, H5T_NATIVE_FLOAT);
        F5Bchart_add_tangential_types(&IDs->DoublePrecision, coordnames, dimension, H5T_NATIVE_DOUBLE);
 
        return IDs;
}
 
ChartDomain_IDs*F5B_new_global_int_chart(const char*coordnames[],
                                           int dimension,
                                           const char*name,
                                           const int perm_vector[FIBER_MAX_RANK])
{
 
ChartDomain_IDs*IDs = F5B_new_global_domain(coordnames, dimension, name, perm_vector, dimension,
                                     H5T_NATIVE_INT32, H5T_NATIVE_INT64);
 
        F5Bchart_add_tangential_types(&IDs->SinglePrecision, coordnames, dimension, H5T_NATIVE_INT8);
        F5Bchart_add_tangential_types(&IDs->DoublePrecision, coordnames, dimension, H5T_NATIVE_INT16);
 
        return IDs;
}
 
ChartDomain_IDs*F5B_new_global_int_domain(const char*coordnames[],
                                          int dimension,
                                          const char*name,
                                          const int perm_vector[FIBER_MAX_RANK],
                                          int cell_dimensionality)
{
ChartDomain_IDs*IDs = F5B_new_global_domain(coordnames,
                                            dimension,
                                            name,
                                            perm_vector,
                                            cell_dimensionality,
                                            H5T_NATIVE_INT32,
                                            H5T_NATIVE_INT64);
 
        F5Bchart_add_tangential_types(&IDs->OneBytePrecision,coordnames, dimension, H5T_NATIVE_INT8 );
        F5Bchart_add_tangential_types(&IDs->HalfPrecision,   coordnames, dimension, H5T_NATIVE_INT16);
        F5Bchart_add_tangential_types(&IDs->SinglePrecision, coordnames, dimension, H5T_NATIVE_INT32);
        F5Bchart_add_tangential_types(&IDs->DoublePrecision, coordnames, dimension, H5T_NATIVE_INT64);
 
        return IDs;
}
 
ChartDomain_IDs*F5B_new_global_fragment_neighbour_domain(const char*coordnames[],
                                                         int dimension,
                                                         const char*name,
                                                         const int perm_vector[],
                                                         int cell_dimensionality)
{
ChartDomain_IDs*IDs = F5B_new_global_domain(coordnames,
                                            dimension,
                                            name,
                                            perm_vector,
                                            cell_dimensionality,
                                            H5T_NATIVE_CHAR,
                                            H5T_NATIVE_CHAR);
/*
        F5Bchart_add_tangential_types(&IDs->SinglePrecision, coordnames, dimension, H5T_NATIVE_INT32);
        F5Bchart_add_tangential_types(&IDs->DoublePrecision, coordnames, dimension, H5T_NATIVE_INT64);
*/
        return IDs;
}
 
ChartDomain_IDs*F5B_new_global_fractional_domain(const char*coordnames[],
                                                 int dimension,
                                                 const char*name,
                                                 const int perm_vector[],
                                                 int cell_dimensionality)
{
ChartDomain_IDs*IDs = F5B_new_global_domain(coordnames,
                                            dimension,
                                            name,
                                            perm_vector,
                                            cell_dimensionality,
                                            F5T_INT_FRACTION32,
                                            F5T_INT_FRACTION64);
        return IDs;
}
 
 
void    F5Binsert_tensor_type(ChartDomain_IDs*This, F5Ttensor_t*T)
{
        assert(This);
        if (!This->TensorTypes)
        {
                This->TensorTypes = new(F5Ttensor_t*, 2);
                This->TensorTypes[0] = T;
                This->TensorTypes[1] = NULL;
        }
        else
        {
        int     i;
                for(i=0; This->TensorTypes[i]; i++)
                        ;
 
                This->TensorTypes = (F5Ttensor_t**)realloc(This->TensorTypes, i+1*sizeof(F5Ttensor_t*));
                This->TensorTypes[i  ] = T;
                This->TensorTypes[i+1] = NULL;
        }
}
 
 
void    F5B_delete_tangential_types(F5_ChartPrecisionTypes*ID)
{
int     i;
        assert(ID);
 
        F5Tclose(ID->Point_hid_t);
        F5Tclose(ID->Vector_hid_t);
        F5Tclose(ID->Covector_hid_t);
        F5Tclose(ID->Bivector_hid_t);
        F5Tclose(ID->Bicovector_hid_t);
        F5Tclose(ID->Pseudoscalar_hid_t);
        F5Tclose(ID->Pseudocoscalar_hid_t);
 
        F5Tclose(ID->Metric_hid_t);
        F5Tclose(ID->Cometric_hid_t);
 
        for(i=0; i<F5ChartPrecisionMaxNumberOfFileTypes; i++)
        {
                F5Tclose( ID->AdditionalFileTypes[ i ] );
        }
}
 
void F5B_delete_global_chart(ChartDomain_IDs*ID)
{
        if (ID && --ID->refs == 0)
        {
                free( (void*)ID->domain_name);
                F5B_delete_tangential_types( &ID->OneBytePrecision );
                F5B_delete_tangential_types( &ID->HalfPrecision );
                F5B_delete_tangential_types( &ID->SinglePrecision );
                F5B_delete_tangential_types( &ID->DoublePrecision );
                F5B_delete_tangential_types( &ID->ExtendedPrecision );
                free(ID);
        }
}
 
#if 0
hid_t F5Tmake(hid_t location,  const char*name, hid_t type, int*committed)
{
herr_t  (*old_func)(void*);
void    *old_client_data;
hid_t   t;
 
        H5Eget_auto(&old_func, &old_client_data);
        H5Eset_auto(NULL, NULL);
 
        t = H5Topen2(location, name, H5P_DEFAULT);
        H5Eset_auto(old_func, old_client_data);
 
        if (t<0)
        {
                t = F5Tcommit( location,  name, type);
                *committed = 1;
        }
        else
                *committed = 0;
 
        return t;
}
#endif
 
#define SAVE_MEMORDER_WITH_TYPE
 
 
void F5Bchart_read_tangential_types(F5_ChartPrecisionTypes*H5Types,
                                    hid_t location)
{
        assert( H5Types );
 
        H5Types->Point_hid_t            = H5Topen2( location,  FIBER_HDF5_CHART_POINT_TYPE, H5P_DEFAULT);
 
        H5E_BEGIN_TRY
        H5Types->Vector_hid_t           = H5Topen2( location, FIBER_HDF5_CHART_VECTOR_TYPE, H5P_DEFAULT);
        H5Types->Covector_hid_t         = H5Topen2( location, FIBER_HDF5_CHART_COVECTOR_TYPE, H5P_DEFAULT);
 
        H5Types->Bivector_hid_t         = H5Topen2( location, FIBER_HDF5_CHART_BIVECTOR_TYPE, H5P_DEFAULT);
        H5Types->Bicovector_hid_t       = H5Topen2( location, FIBER_HDF5_CHART_BICOVECTOR_TYPE, H5P_DEFAULT);
 
        H5Types->Pseudoscalar_hid_t     = H5Topen2( location, FIBER_HDF5_CHART_PSEUDOSCALAR_TYPE, H5P_DEFAULT);
        H5Types->Pseudocoscalar_hid_t   = H5Topen2( location, FIBER_HDF5_CHART_PSEUDOCOSCALAR_TYPE, H5P_DEFAULT);
 
        H5Types->Metric_hid_t           = H5Topen2( location, FIBER_HDF5_CHART_METRIC_TYPE, H5P_DEFAULT);
        H5Types->Cometric_hid_t         = H5Topen2( location, FIBER_HDF5_CHART_METRIC_TYPE, H5P_DEFAULT);
        H5E_END_TRY
}
 
ChartDomain_IDs*F5B_read_global_chart(hid_t ChartsID, const char*name)
{
hid_t   location;
ChartDomain_IDs*FileID = 0;
 
        F5printf(200, "F5B_read_global_chart(,%s)\n", name);
 
        location = H5Gopen2(ChartsID, name, H5P_DEFAULT);
        if (location<0)
        {
                F5printf(0, "F5B_read_global_chart(,%s) FAILED\n", name);
                return 0;
        }
 
        FileID = F5B_read_chart_domain(location, name);
        H5Gclose(location);
 
        return FileID;
}
 
/**
@todo Make failsafe if a type definition does not exist.
@todo Add support for dynamically defined types.
*/
ChartDomain_IDs*F5B_read_chart_domain(hid_t ChartDomainLocationID, const char*DefaultChartName)
{
int     Tdims;
hid_t   perm_id,  precision_group;
ChartDomain_IDs*FileID = 0;
 
        FileID = (ChartDomain_IDs*)malloc(sizeof(ChartDomain_IDs));
        memset (FileID, 0, sizeof (ChartDomain_IDs));
 
        FileID->refs = 1;
 
        precision_group = F5Gtry_to_open(ChartDomainLocationID, FIBER_HDF5_CHART_SINGLE_PRECISION);
        if (precision_group>0)
        {
                F5Bchart_read_tangential_types( &FileID->SinglePrecision, precision_group);
                H5Gclose(precision_group);
                precision_group = H5Gopen2(ChartDomainLocationID, FIBER_HDF5_CHART_DOUBLE_PRECISION, H5P_DEFAULT);
                F5Bchart_read_tangential_types( &FileID->DoublePrecision, precision_group);
                H5Gclose(precision_group);
        }
        else
        {
                F5Bchart_read_tangential_types( &FileID->SinglePrecision, ChartDomainLocationID );
                F5Bchart_read_tangential_types( &FileID->DoublePrecision, ChartDomainLocationID );
        }
 
        {
        hid_t   Point_hid_t = FileID->SinglePrecision.Point_hid_t;
                if (Point_hid_t<0)
                {
                        F5printf(0, "F5B_read_chart_domain() Invalid single precision point type!\n");
                }
                else
                {
                char*buf = (char*)malloc(256);
                        if (!F5Aget_string( Point_hid_t, FIBER_HDF5_CHART_DOMAIN, buf, 256 ) )
                        {
                                free(buf);
                                FileID->domain_name = strdup(DefaultChartName);
                        }
                        else
                                FileID->domain_name = buf;
                }
 
                Tdims = H5Tget_nmembers( Point_hid_t );
 
 
/*#ifdef        SAVE_MEMORDER_WITH_TYPE*/
                perm_id = F5Atry_to_open(Point_hid_t, FIBER_HDF5_CHART_MEMLAYOUT);
/*
#else
        perm_id = F5Atry_to_open(location, FIBER_HDF5_CHART_MEMLAYOUT);
#endif
*/
        }
 
        if (perm_id>0)
        {
                assert(FileID->perm_vector);
                F_H5Aread(perm_id, H5T_NATIVE_INT, FileID->perm_vector, FIBER_HDF5_CHART_MEMLAYOUT);
                H5Aclose(perm_id);
        }
        else /* NO MEMLAYOUT attribute: Default is a FORTRAN permutation vector */
        {
        int i;
                assert(FileID->perm_vector);
                for(i=0; i<Tdims; i++)
                        FileID->perm_vector[i] = Tdims - 1 - i;
        }
 
        F5printf(200, "~F5B_read_chart_domain(,%s) DONE\n", DefaultChartName);
 
        return FileID;
}
 
 
/**
   Create File ID's from MemoryID's.
 
@todo Possible optimization: Save only if perm_vector is not FORTRAN ordering .
*/
ChartDomain_IDs*F5B_save_global_chart(const ChartDomain_IDs*MemoryIDs, hid_t Domain_hid)
{
int     i;
ChartDomain_IDs*FileID = (ChartDomain_IDs*)malloc(sizeof(ChartDomain_IDs));
        memset (FileID, 0, sizeof (ChartDomain_IDs));
 
        FileID->refs = 1;
        for(i=0; i<FIBER_MAX_RANK; i++)
                FileID->perm_vector[i] = MemoryIDs->perm_vector[i];
 
        FileID->domain_name = strdup(MemoryIDs->domain_name);
        FileID->cell_dimensionality = MemoryIDs->cell_dimensionality;
 
        /* Save the memory layout ordering attribute */
        /* Possible optimization: Save only if perm_vector is not FORTRAN ordering */
 
        {
        hid_t   SinglePrecisionTypes = H5Gcreate2(Domain_hid, FIBER_HDF5_CHART_SINGLE_PRECISION,
                                                  H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT );
 
                FileID->SinglePrecision.Point_hid_t = F5Tsave_tensor(SinglePrecisionTypes, MemoryIDs->domain_name,
                                                                     FIBER_HDF5_CHART_POINT_TYPE, MemoryIDs->SinglePrecision.Point_hid_t, 0,0, 0);
 
                H5Gclose( SinglePrecisionTypes );
        }
        {
        hid_t   DoublePrecisionTypes = H5Gcreate2(Domain_hid, FIBER_HDF5_CHART_DOUBLE_PRECISION,
                                                  H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT );
 
                FileID->DoublePrecision.Point_hid_t = F5Tsave_tensor(DoublePrecisionTypes, MemoryIDs->domain_name,
                                                                     FIBER_HDF5_CHART_POINT_TYPE, MemoryIDs->DoublePrecision.Point_hid_t, 0,0, 0);
 
                H5Gclose( DoublePrecisionTypes );
        }
 
        /* by default, link to single precision type */
        H5Lcreate_soft( FIBER_HDF5_CHART_SINGLE_PRECISION "/" FIBER_HDF5_CHART_POINT_TYPE,
                        Domain_hid, FIBER_HDF5_CHART_POINT_TYPE, H5P_DEFAULT, H5P_DEFAULT );
 
        {
        int     cd = FileID->cell_dimensionality;
                F5Asave_ints(FileID->SinglePrecision.Point_hid_t, FIBER_HDF5_CELLDIMENSIONALITY, &cd, 1);
                F5Asave_ints(FileID->DoublePrecision.Point_hid_t, FIBER_HDF5_CELLDIMENSIONALITY, &cd, 1);
        }
 
        F5T_save_perm_vector(FileID->SinglePrecision.Point_hid_t, H5Tget_nmembers( FileID->SinglePrecision.Point_hid_t ), FileID->perm_vector);
        F5T_save_perm_vector(FileID->DoublePrecision.Point_hid_t, H5Tget_nmembers( FileID->DoublePrecision.Point_hid_t ), FileID->perm_vector);
        return FileID;
}
 
ChartDomain_IDs*F5Bcopy_chart_domain(const ChartDomain_IDs*srcIDs)
{
int     i;
const char*coordnames[FIBER_MAX_RANK];
        assert( srcIDs );
        {
        int     dimension = H5Tget_nmembers( srcIDs->SinglePrecision.Point_hid_t );
/*hid_t coord_type = H5Tget_member_type( srcIDs->SinglePrecision.Point_hid_t, 0 ); */
        ChartDomain_IDs*newID;
        const char**cn = coordnames;
 
                for(i=0; i<dimension; i++)
                        coordnames[i] = H5Tget_member_name( srcIDs->SinglePrecision.Point_hid_t, i );
 
/*      newID = F5B_new_global_chart(cn, dimension, coord_type, srcIDs->domain_name, srcIDs->perm_vector, srcIDs->cell_dimensionality); */
                newID = F5B_new_global_domain(cn, dimension, srcIDs->domain_name, srcIDs->perm_vector, srcIDs->cell_dimensionality,
                                              H5Tget_member_type( srcIDs->SinglePrecision.Point_hid_t, 0 ),
                                              H5Tget_member_type( srcIDs->DoublePrecision.Point_hid_t, 0 )
                        );
 
                for(i=0; i<dimension; i++)
                        H5free_memory( (char*)(coordnames[i]) );
 
                return newID;
        }
}
 
 
void F5T_save_perm_vector(hid_t location, int rank, int perm_vector[])
{
        if (!perm_vector)
        {
                puts("F5T_save_perm_vector()");
                return;
        }
 
        /* Save the memory layout ordering attribute */
        /* Possible optimization: Save only if perm_vector is not FORTRAN ordering */
        {
        hsize_t Tdims = rank;
        hid_t   index_ds = H5Screate_simple(1, &Tdims, &Tdims ) ;
        hid_t   Rank_id = F5Acreate(location, FIBER_HDF5_CHART_MEMLAYOUT, H5T_NATIVE_INT, index_ds, 0);
                H5Awrite(Rank_id, H5T_NATIVE_INT, perm_vector );
                H5Aclose(Rank_id);
                H5Sclose(index_ds);
        }
}
 
int     F5Bget_chart_dims(const ChartDomain_IDs*IDs)
{
        if (!IDs)
                return -1;
 
        return H5Tget_nmembers(IDs->SinglePrecision.Point_hid_t);
}
 
 
static  int F5Tequal(hid_t A, hid_t B)
{
        if (A<=0) return 0;
        if (B<=0) return 0;
 
        return H5Tequal(A,B)>0;
}
 
F5_ChartPrecisionTypes*F5Bget_precision_by_size(ChartDomain_IDs*CD, int TypeSize)
{
        if (!CD) return 0;
 
        switch(TypeSize)
        {
        case 1:   return &CD->OneBytePrecision;
        case 2:   return &CD->HalfPrecision;
        case 4:   return &CD->SinglePrecision;
        case 8:   return &CD->DoublePrecision;
        case 16:  return &CD->ExtendedPrecision;
        default:  return 0;
        }
}
 
F5_ChartPrecisionTypes*F5Bget_precision_by_type(ChartDomain_IDs*CD, hid_t type_id)
{
hid_t   FirstCompoundType = H5Tget_member_type( type_id, 0);
int     CompoundSize = H5Tget_size( FirstCompoundType );
        H5Tclose( FirstCompoundType );
 
        return F5Bget_precision_by_size(CD, CompoundSize);
}
 
 
hid_t   F5file_type(F5Path*fpath, hid_t fieldtype)
{
        if (!fpath || !fpath->myChart)
        {
                F5printf(20, "F5file_type: No chart given, returning memory type.");
                return fieldtype;
        }
 
        F5printf(40, "F5file_type(%p, %d): chart=%p", fpath, (int) fieldtype, fpath->myChart );
 
        switch( H5Tget_class( fieldtype) )
        {
        default:
                F5printf(20, "F5file_type query for unknown type (could be a native one).");
                return fieldtype;
 
        case H5T_INTEGER:
        case H5T_FLOAT:
        case H5T_STRING:
                F5printf(20, "F5file_type query for integer/float/string type (could be a native one).");
                return fieldtype;
 
        case H5T_BITFIELD:
        case H5T_OPAQUE:
        case H5T_REFERENCE:
        case H5T_ENUM:
        case H5T_VLEN:
        case H5T_ARRAY:
                F5printf(20, "F5file_type query for unsupported type (could be a native one).");
                return fieldtype;
 
        case H5T_COMPOUND:
                ;
        }
 
/*
  The following code is executed once we know the fieldtype is a compound type.
 */
 
        {
        hid_t   FirstCompoundType = H5Tget_member_type( fieldtype, 0);
        int     CompoundPrecision = H5Tget_size( FirstCompoundType );
 
        F5_ChartPrecisionTypes*IN   = F5Bget_precision_by_size(fpath->myChart, CompoundPrecision),
                              *OUT  = F5Bget_precision_by_size(fpath->FileIDs, CompoundPrecision);
 
                H5Tclose( FirstCompoundType );
 
                F5printf(20, "F5file_type query for compound type, precision %d.\n", CompoundPrecision);
 
                assert(IN);
                assert(OUT);
                assert( fpath->myChart );
                assert( fpath->FileIDs );
 
                if (F5Tequal(fieldtype, IN->Point_hid_t) )
                {
                        F5printf(20, "F5file_type query for point type.\n");
                        if (OUT->Point_hid_t<1)
                        {
                                F5printf(0, "F5file_type(): Possible internal error, point type not found in file!");
 
                                OUT->Point_hid_t = F5Tmake_precise_tensor(fpath, FIBER_HDF5_CHART_POINT_TYPE, CompoundPrecision,
                                                                          IN->Point_hid_t, 0,0, 0);
                                {
                                int cd = fpath->myChart->cell_dimensionality;
                                        printf("CELL Dimensionality: %d\n", cd);
                                        F5Asave_ints(OUT->Point_hid_t,  FIBER_HDF5_CELLDIMENSIONALITY, &cd, 1);
                                }
 
#ifdef  SAVE_MEMORDER_WITH_TYPE
                                F5T_save_perm_vector(OUT->Point_hid_t,
                                                     H5Tget_nmembers( OUT->Point_hid_t ),
                                                     fpath->FileIDs->perm_vector);
#endif
                        }
                        return OUT->Point_hid_t;
                }
                else if (F5Tequal(fieldtype, IN->Vector_hid_t) )
                {
                        F5printf(20, "F5file_type query for vector type %d", (int) IN->Vector_hid_t );
                        if (OUT->Vector_hid_t<1)
                        {
                        int     covariance[] = {1}; /* Tangential vectors have one upper index */
                                OUT->Vector_hid_t =
                                        F5Tmake_precise_tensor(fpath, FIBER_HDF5_CHART_VECTOR_TYPE, CompoundPrecision,
                                                               IN->Vector_hid_t, 1, covariance, 1);
                        }
                        return OUT->Vector_hid_t;
                }
                else if (F5Tequal(fieldtype, IN->Covector_hid_t) )
                {
                        F5printf(20, "F5file_type query for covector type.\n");
                        if (OUT->Covector_hid_t<1)
                        {
                        int     covariance[] = {-1}; /* Co-vectors have one lower index */
                                OUT->Covector_hid_t =
                                        F5Tmake_precise_tensor(fpath, FIBER_HDF5_CHART_COVECTOR_TYPE, CompoundPrecision,
                                                               IN->Covector_hid_t, 1, covariance, 1);
                        }
                        return OUT->Covector_hid_t;
                }
                else if (F5Tequal(fieldtype, IN->Bivector_hid_t) )
                {
                        F5printf(20, "F5file_type query for bivector type.\n");
                        if (OUT->Bivector_hid_t<1)
                        {
                        int     covariance[] = {1,1}; /* bivectors have two upper indices */
                                OUT->Bivector_hid_t =
                                        F5Tmake_precise_tensor(fpath, FIBER_HDF5_CHART_BIVECTOR_TYPE, CompoundPrecision,
                                                               IN->Bivector_hid_t, 2, covariance, 2);
                        }
                        return OUT->Bivector_hid_t;
                }
                else if (F5Tequal(fieldtype, IN->Bicovector_hid_t) )
                {
                        F5printf(20, "F5file_type query for bicovector type.\n");
                        if (OUT->Bivector_hid_t<1)
                        {
                        int     covariance[] = {-1,-1}; /* bicovectors have two lower indices */
                                OUT->Bicovector_hid_t =
                                        F5Tmake_precise_tensor(fpath, FIBER_HDF5_CHART_BICOVECTOR_TYPE, CompoundPrecision,
                                                               IN->Bicovector_hid_t, 2, covariance, 2);
                        }
                        return OUT->Bicovector_hid_t;
                }
 
                else if (F5Tequal(fieldtype, IN->Metric_hid_t) )
                {
                        F5printf(20, "F5file_type query for metric type.\n");
                        if (OUT->Metric_hid_t<1)
                        {
                                /* g_\mu\nu */
                        int     covariance[] = {-1, -1}; /* Metric tensors have two lower indices */
                                OUT->Metric_hid_t =
                                        F5Tmake_precise_tensor(fpath, FIBER_HDF5_CHART_METRIC_TYPE, CompoundPrecision,
                                                               IN->Metric_hid_t, 2, covariance, 0);
                        }
                        return OUT->Metric_hid_t;
                }
                else if (F5Tequal(fieldtype, IN->Cometric_hid_t) )
                {
                        F5printf(20, "F5file_type query for cometric type.\n");
                        if (OUT->Cometric_hid_t<1)
                        {
                                /* g^\mu\nu */
                        int     covariance[] = {1, 1}; /* Cometric tensors have two upper indices */
                                OUT->Cometric_hid_t =
                                        F5Tmake_precise_tensor(fpath, FIBER_HDF5_CHART_COMETRIC_TYPE, CompoundPrecision,
                                                               IN->Cometric_hid_t, 2, covariance, 0);
                        }
                        return OUT->Cometric_hid_t;
                }
 
 
                F5printf(2, __FILE__ ": Warning: F5file_type() query could not identify file type.\n");
                {
                int     TypeMembers = H5Tget_nmembers( fieldtype );
                const   char*MultiplicityCategory;
                hid_t   MultiplicityCategory_id;
                        switch(TypeMembers)
                        {
                        case 1: MultiplicityCategory = "Chart-free compound types with one member"; break;
                        case 2: MultiplicityCategory = "Chart-free compound types with two members"; break;
                        case 3: MultiplicityCategory = "Chart-free compound types with three members"; break;
                        case 4: MultiplicityCategory = "Chart-free compound types with four members"; break;
                        default:MultiplicityCategory = "Chart-free compound types with many members"; break;
                        }
                        MultiplicityCategory_id = F5Gappend(fpath->GlobalChart_hid, MultiplicityCategory);
 
                        {
#define FIBER_HDF5_ANONYMOUS_TYPENAME "F5::TypeName%04u"
                        char    type_name[1024];
                        unsigned nthType = 0;
                                for(;nthType<10000;nthType++)
                                {
                                hid_t   type_id;
                                        sprintf(type_name, FIBER_HDF5_ANONYMOUS_TYPENAME, nthType);
                                        type_id = F5Tappend(MultiplicityCategory_id, type_name, fieldtype);
 
                                        if (type_id>0)
                                        {
                                                /*
                                                  Inject This type to the chart file identifiers
                                                  for the output file, for mere cleanup upon file
                                                  closing. Otherwise This type ID will NOT be associated
                                                  with the respective chart.
                                                 */
                                                F5B_inject_precision_type(OUT, type_id);
                                                H5Gclose(MultiplicityCategory_id);
                                                return type_id;
                                        }
                                }
                                puts("F5file_type: Creating anonymous type name failed!");
                                H5Gclose(MultiplicityCategory_id);
                                return fieldtype;
                        }
                }
        }
}