v1hpg.c

一个用来实现偏微分方程中网格的计算库

		}
	}
	return xlen;
}


/* Write a NC_attrarray to the header */
static int
v1h_put_NC_attrarray(v1hs *psp, const NC_attrarray *ncap)
{
	int status;

	assert(psp != NULL);

	if(ncap == NULL
#if 1
		/* Backward:
		 * This clause is for 'byte for byte'
		 * backward compatibility.
		 * Strickly speaking, it is 'bug for bug'.
		 */
		|| ncap->nelems == 0
#endif
		)
	{
		/*
		 * Handle empty netcdf
		 */
		const size_t nosz = 0;

		status = v1h_put_NCtype(psp, NC_UNSPECIFIED);
		if(status != ENOERR)
			return status;
		status = v1h_put_size_t(psp, &nosz);
		if(status != ENOERR)
			return status;
		return ENOERR;
	}
	/* else */

	status = v1h_put_NCtype(psp, NC_ATTRIBUTE);
	if(status != ENOERR)
		return status;
	status = v1h_put_size_t(psp, &ncap->nelems);
	if(status != ENOERR)
		return status;

	{
		const NC_attr **app = (const NC_attr **)ncap->value;
		const NC_attr *const *const end = &app[ncap->nelems];
		for( /*NADA*/; app < end; app++)
		{
			status = v1h_put_NC_attr(psp, *app);
			if(status)
				return status;
		}
	}
	return ENOERR;
}


/* Read a NC_attrarray from the header */
static int
v1h_get_NC_attrarray(v1hs *gsp, NC_attrarray *ncap)
{
	int status;
	NCtype type = NC_UNSPECIFIED;

	assert(gsp != NULL && gsp->pos != NULL);
	assert(ncap != NULL);
	assert(ncap->value == NULL);

	status = v1h_get_NCtype(gsp, &type);
	if(status != ENOERR)
		return status;
	status = v1h_get_size_t(gsp, &ncap->nelems);
	if(status != ENOERR)
		return status;
	
	if(ncap->nelems == 0)
		return ENOERR;
	/* else */
	if(type != NC_ATTRIBUTE)
		return EINVAL;

	ncap->value = (NC_attr **) malloc(ncap->nelems * sizeof(NC_attr *));
	if(ncap->value == NULL)
		return NC_ENOMEM;
	ncap->nalloc = ncap->nelems;

	{
		NC_attr **app = ncap->value;
		NC_attr *const *const end = &app[ncap->nelems];
		for( /*NADA*/; app < end; app++)
		{
			status = v1h_get_NC_attr(gsp, app);
			if(status)
			{
				ncap->nelems = (size_t)(app - ncap->value);
				free_NC_attrarrayV(ncap);
				return status;
			}
		}
	}

	return ENOERR;
}

/* End NC_attr */
/* Begin NC_var */

/*
 * How much space will the xdr'd var take.
 * Formerly
NC_xlen_var(vpp)
 */
static size_t
ncx_len_NC_var(const NC_var *varp, size_t sizeof_off_t)
{
	size_t sz;

	assert(varp != NULL);
	assert(sizeof_off_t != 0);

	sz = ncx_len_NC_string(varp->name);
	sz += X_SIZEOF_SIZE_T; /* ndims */
	sz += ncx_len_int(varp->ndims); /* dimids */
	sz += ncx_len_NC_attrarray(&varp->attrs);
	sz += X_SIZEOF_NC_TYPE; /* type */
	sz += X_SIZEOF_SIZE_T; /* len */
	sz += sizeof_off_t; /* begin */

	return(sz);
}


/* Write a NC_var to the header */
static int
v1h_put_NC_var(v1hs *psp, const NC_var *varp)
{
	int status;

	status = v1h_put_NC_string(psp, varp->name);
	if(status != ENOERR)
		return status;

	status = v1h_put_size_t(psp, &varp->ndims);
	if(status != ENOERR)
		return status;

	status = check_v1hs(psp, ncx_len_int(varp->ndims));
	if(status != ENOERR)
		return status;
	status = ncx_putn_int_int(&psp->pos,
			varp->ndims, varp->dimids);
	if(status != ENOERR)
		return status;

	status = v1h_put_NC_attrarray(psp, &varp->attrs);
	if(status != ENOERR)
		return status;

	status = v1h_put_nc_type(psp, &varp->type);
	if(status != ENOERR)
		return status;

	status = v1h_put_size_t(psp, &varp->len);
	if(status != ENOERR)
		return status;

	status = check_v1hs(psp, psp->version == 1 ? 4 : 8);
	if(status != ENOERR)
		 return status;
	status = ncx_put_off_t(&psp->pos, &varp->begin, psp->version == 1 ? 4 : 8);
	if(status != ENOERR)
		return status;

	return ENOERR;
}


/* Read a NC_var from the header */
static int
v1h_get_NC_var(v1hs *gsp, NC_var **varpp)
{
	NC_string *strp;
	int status;
	size_t ndims;
	NC_var *varp;

	status = v1h_get_NC_string(gsp, &strp);
	if(status != ENOERR)
		return status;

	status = v1h_get_size_t(gsp, &ndims);
	if(status != ENOERR)
		goto unwind_name;

	varp = new_x_NC_var(strp, ndims);
	if(varp == NULL)
	{
		status = NC_ENOMEM;
		goto unwind_name;
	}

	status = check_v1hs(gsp, ncx_len_int(ndims));
	if(status != ENOERR)
		goto unwind_alloc;
	status = ncx_getn_int_int((const void **)(&gsp->pos),
			ndims, varp->dimids);
	if(status != ENOERR)
		goto unwind_alloc;

	status = v1h_get_NC_attrarray(gsp, &varp->attrs);
	if(status != ENOERR)
		goto unwind_alloc;

	status = v1h_get_nc_type(gsp, &varp->type);
	if(status != ENOERR)
		 goto unwind_alloc;

	status = v1h_get_size_t(gsp, &varp->len);
	if(status != ENOERR)
		 goto unwind_alloc;

	status = check_v1hs(gsp, gsp->version == 1 ? 4 : 8);
	if(status != ENOERR)
		 goto unwind_alloc;
	status = ncx_get_off_t((const void **)&gsp->pos,
			       &varp->begin, gsp->version == 1 ? 4 : 8);
	if(status != ENOERR)
		 goto unwind_alloc;
	
	*varpp = varp;
	return ENOERR;

unwind_alloc:
	free_NC_var(varp); /* frees name */
	return status;

unwind_name:
	free_NC_string(strp);
	return status;
}


/* How much space in the header is required for this NC_vararray? */
static size_t
ncx_len_NC_vararray(const NC_vararray *ncap, size_t sizeof_off_t)
{
	size_t xlen = X_SIZEOF_NCTYPE;	/* type */
	xlen += X_SIZEOF_SIZE_T;	/* count */
	if(ncap == NULL)
		return xlen;
	/* else */
	{
		const NC_var **vpp = (const NC_var **)ncap->value;
		const NC_var *const *const end = &vpp[ncap->nelems];
		for( /*NADA*/; vpp < end; vpp++)
		{
			xlen += ncx_len_NC_var(*vpp, sizeof_off_t);
		}
	}
	return xlen;
}


/* Write a NC_vararray to the header */
static int
v1h_put_NC_vararray(v1hs *psp, const NC_vararray *ncap)
{
	int status;

	assert(psp != NULL);

	if(ncap == NULL
#if 1
		/* Backward:
		 * This clause is for 'byte for byte'
		 * backward compatibility.
		 * Strickly speaking, it is 'bug for bug'.
		 */
		|| ncap->nelems == 0
#endif
		)
	{
		/*
		 * Handle empty netcdf
		 */
		const size_t nosz = 0;

		status = v1h_put_NCtype(psp, NC_UNSPECIFIED);
		if(status != ENOERR)
			return status;
		status = v1h_put_size_t(psp, &nosz);
		if(status != ENOERR)
			return status;
		return ENOERR;
	}
	/* else */

	status = v1h_put_NCtype(psp, NC_VARIABLE);
	if(status != ENOERR)
		return status;
	status = v1h_put_size_t(psp, &ncap->nelems);
	if(status != ENOERR)
		return status;

	{
		const NC_var **vpp = (const NC_var **)ncap->value;
		const NC_var *const *const end = &vpp[ncap->nelems];
		for( /*NADA*/; vpp < end; vpp++)
		{
			status = v1h_put_NC_var(psp, *vpp);
			if(status)
				return status;
		}
	}
	return ENOERR;
}


/* Read a NC_vararray from the header */
static int
v1h_get_NC_vararray(v1hs *gsp, NC_vararray *ncap)
{
	int status;
	NCtype type = NC_UNSPECIFIED;

	assert(gsp != NULL && gsp->pos != NULL);
	assert(ncap != NULL);
	assert(ncap->value == NULL);

	status = v1h_get_NCtype(gsp, &type);
	if(status != ENOERR)
		return status;
	
	status = v1h_get_size_t(gsp, &ncap->nelems);
	if(status != ENOERR)
		return status;
	
	if(ncap->nelems == 0)
		return ENOERR;
	/* else */
	if(type != NC_VARIABLE)
		return EINVAL;

	ncap->value = (NC_var **) malloc(ncap->nelems * sizeof(NC_var *));
	if(ncap->value == NULL)
		return NC_ENOMEM;
	ncap->nalloc = ncap->nelems;

	{
		NC_var **vpp = ncap->value;
		NC_var *const *const end = &vpp[ncap->nelems];
		for( /*NADA*/; vpp < end; vpp++)
		{
			status = v1h_get_NC_var(gsp, vpp);
			if(status)
			{
				ncap->nelems = (size_t)(vpp - ncap->value);
				free_NC_vararrayV(ncap);
				return status;
			}
		}
	}

	return ENOERR;
}


/* End NC_var */
/* Begin NC */


/*
 * Recompute the shapes of all variables
 * Sets ncp->begin_var to start of first variable.
 * Sets ncp->begin_rec to start of first record variable.
 * Returns -1 on error. The only possible error is a reference
 * to a non existent dimension, which could occur for a corrupted
 * netcdf file.
 */
static int
NC_computeshapes(NC *ncp)
{
	NC_var **vpp = (NC_var **)ncp->vars.value;
	NC_var *const *const end = &vpp[ncp->vars.nelems];
	NC_var *first_var = NULL;	/* first "non-record" var */
	NC_var *first_rec = NULL;	/* first "record" var */
	int status;

	ncp->begin_var = (off_t) ncp->xsz;
	ncp->begin_rec = (off_t) ncp->xsz;
	ncp->recsize = 0;

	if(ncp->vars.nelems == 0)
		return(0);
	
	for( /*NADA*/; vpp < end; vpp++)
	{
		status = NC_var_shape(*vpp, &ncp->dims);
		if(status != ENOERR)
			return(status);

	  	if(IS_RECVAR(*vpp))	
		{
	  		if(first_rec == NULL)	
				first_rec = *vpp;
			ncp->recsize += (*vpp)->len;
		}
		else
		{
		        if(first_var == NULL)
			        first_var = *vpp;
			/*
			 * Overwritten each time thru.
			 * Usually overwritten in first_rec != NULL clause below.
			 */
			ncp->begin_rec = (*vpp)->begin + (off_t)(*vpp)->len;
		}
	}

	if(first_rec != NULL)
	{
		assert(ncp->begin_rec <= first_rec->begin);
		ncp->begin_rec = first_rec->begin;
		/*
	 	 * for special case of exactly one record variable, pack value
	 	 */
		if(ncp->recsize == first_rec->len)
			ncp->recsize = *first_rec->dsizes * first_rec->xsz;
	}

	if(first_var != NULL)
	{
		ncp->begin_var = first_var->begin;
	}
	else
	{
		ncp->begin_var = ncp->begin_rec;
	}

	assert(ncp->begin_var > 0);
	assert(ncp->xsz <= (size_t)ncp->begin_var);
	assert(ncp->begin_rec > 0);
	assert(ncp->begin_var <= ncp->begin_rec);
	
	return(ENOERR);
}


/*
 * Return actual unpadded length (in bytes) of a variable, which
 * doesn't include any extra padding used for alignment.  For a record
 * variable, this is the length in bytes of one record's worth of that
 * variable's data.
 */
static off_t
NC_var_unpadded_len(const NC_var *varp, const NC_dimarray *dims)
{
    size_t *shp;
    off_t product = 1;
    
    if(varp->ndims != 0) {
	for(shp = varp->shape + varp->ndims -1; shp >= varp->shape; shp--) {
	    if(!(shp == varp->shape && IS_RECVAR(varp)))
		product *= *shp;
	}
    }
    product = product * varp->xsz;
    return product;
}

/* How much space in the header is required for the NC data structure? */
size_t
ncx_len_NC(const NC *ncp, size_t sizeof_off_t)
{
	size_t xlen = sizeof(ncmagic);

	assert(ncp != NULL);
	
	xlen += X_SIZEOF_SIZE_T; /* numrecs */
	xlen += ncx_len_NC_dimarray(&ncp->dims);
	xlen += ncx_len_NC_attrarray(&ncp->attrs);
	xlen += ncx_len_NC_vararray(&ncp->vars, sizeof_off_t);

	return xlen;
}


/* Write the file header */
int
ncx_put_NC(const NC *ncp, void **xpp, off_t offset, size_t extent)
{
	int status = ENOERR;
	v1hs ps; /* the get stream */

	assert(ncp != NULL);

	/* Initialize stream ps */

	ps.nciop = ncp->nciop;
	ps.flags = RGN_WRITE;

	if (ncp->flags & NC_64BIT_OFFSET)
	  ps.version = 2;
	else 
	  ps.version = 1;

	if(xpp == NULL)
	{
		/*
		 * Come up with a reasonable stream read size.
		 */
		extent = ncp->xsz;
		if(extent <= MIN_NC_XSZ)
		{
			/* first time read */
			extent = ncp->chunk;
			/* Protection for when ncp->chunk is huge;
			 * no need to read hugely. */
	      		if(extent > 4096)
				extent = 4096;
		}
		else if(extent > ncp->chunk)
		{
			extent = ncp->chunk;
		}
		
		ps.offset = 0;
		ps.extent = extent;
		ps.base = NULL;
		ps.pos = ps.base;

		status = fault_v1hs(&ps, extent);
		if(status)
			return status;
	}
	else
	{
		ps.offset = offset;
		ps.extent = extent;
		ps.base = *xpp;
		ps.pos = ps.base;
		ps.end = (char *)ps.base + ps.extent;
	}

	if (ps.version == 2)
	  status = ncx_putn_schar_schar(&ps.pos, sizeof(ncmagic), ncmagic);
	else
	  status = ncx_putn_schar_schar(&ps.pos, sizeof(ncmagic1), ncmagic1);
	if(status != ENOERR)
		goto release;

	{
	const size_t nrecs = NC_get_numrecs(ncp);
	status = ncx_put_size_t(&ps.pos, &nrecs);
	if(status != ENOERR)
		goto release;
	}

	assert((char *)ps.pos < (char *)ps.end);

	status = v1h_put_NC_dimarray(&ps, &ncp->dims);
	if(status != ENOERR)
		goto release;

	status = v1h_put_NC_attrarray(&ps, &ncp->attrs);
	if(status != ENOERR)
		goto release;

	status = v1h_put_NC_vararray(&ps, &ncp->vars);
	if(status != ENOERR)
		goto release;

release:
	(void) rel_v1hs(&ps);

	return status;
}


/* Make the in-memory NC structure from reading the file header */
int
nc_get_NC(NC *ncp)
{
	int status;
	v1hs gs; /* the get stream */

	assert(ncp != NULL);

	/* Initialize stream gs */

	gs.nciop = ncp->nciop;
	gs.offset = 0; /* beginning of file */
	gs.extent = 0;
	gs.flags = 0;
	gs.version = 0;
	gs.base = NULL;
	gs.pos = gs.base;

	{
		/*
		 * Come up with a reasonable stream read size.
		 */
	        off_t filesize;
		size_t extent = MIN_NC_XSZ;
		
		extent = ncp->xsz;
		if(extent <= MIN_NC_XSZ)
		{
		        status = ncio_filesize(ncp->nciop, &filesize);
			if(status)
			    return status;
			/* first time read */
			extent = ncp->chunk;
			/* Protection for when ncp->chunk is huge;
			 * no need to read hugely. */
	      		if(extent > 4096)
				extent = 4096;
			if(extent > filesize)
			        extent = filesize;
		}
		else if(extent > ncp->chunk)
		{
			extent = ncp->chunk;
		}

		/*
		 * Invalidate the I/O buffers to force a read of the header
		 * region.
		 */
		status = gs.nciop->sync(gs.nciop);
		if(status)
			return status;

		status = fault_v1hs(&gs, extent);
		if(status)
			return status;
	}

	/* get the header from the stream gs */

	{
		/* Get & check magic number */
		schar magic[sizeof(ncmagic)];
		(void) memset(magic, 0, sizeof(magic));

		status = ncx_getn_schar_schar(
			(const void **)(&gs.pos), sizeof(magic), magic);
		if(status != ENOERR)
			goto unwind_get;
	
		if(memcmp(magic, ncmagic, sizeof(ncmagic)-1) != 0)
		{
			status = NC_ENOTNC;
			goto unwind_get;
		}
		/* Check version number in last byte of magic */
		if (magic[sizeof(ncmagic)-1] == 0x1) {
		  gs.version = 1;
		} else if (magic[sizeof(ncmagic)-1] == 0x2) {
		  gs.version = 2;
		  fSet(ncp->flags, NC_64BIT_OFFSET);
		  /* Now we support version 2 file access on non-LFS systems -- rkr */
#if 0
		  if (sizeof(off_t) != 8) {
		    fprintf(stderr, "NETCDF WARNING: Version 2 file on 32-bit system.\n");
		  }
#endif
		} else {
			status = NC_ENOTNC;
			goto unwind_get;
		}
	}
	
	{
	size_t nrecs = 0;
	status = ncx_get_size_t((const void **)(&gs.pos), &nrecs);
	if(status != ENOERR)
		goto unwind_get;
	NC_set_numrecs(ncp, nrecs);
	}

	assert((char *)gs.pos < (char *)gs.end);

	status = v1h_get_NC_dimarray(&gs, &ncp->dims);
	if(status != ENOERR)
		goto unwind_get;

	status = v1h_get_NC_attrarray(&gs, &ncp->attrs);
	if(status != ENOERR)
		goto unwind_get;

	status = v1h_get_NC_vararray(&gs, &ncp->vars);
	if(status != ENOERR)
		goto unwind_get;
		
	ncp->xsz = ncx_len_NC(ncp, (gs.version == 1) ? 4 : 8);

	status = NC_computeshapes(ncp);
	if(status != ENOERR)
		goto unwind_get;

unwind_get:
	(void) rel_v1hs(&gs);
	return status;
}