/*============================================================================
* Main structure for a nodal representation associated with a mesh
*============================================================================*/
/*
This file is part of the "Finite Volume Mesh" library, intended to provide
finite volume mesh and associated fields I/O and manipulation services.
Copyright (C) 2004-2006 EDF
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
/*----------------------------------------------------------------------------
* Standard C library headers
*----------------------------------------------------------------------------*/
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
/*----------------------------------------------------------------------------
* BFT library headers
*----------------------------------------------------------------------------*/
#include <bft_mem.h>
#include <bft_printf.h>
/*----------------------------------------------------------------------------
* Local headers
*----------------------------------------------------------------------------*/
#include <fvm_config_defs.h>
#include <fvm_defs.h>
#include <fvm_io_num.h>
#include <fvm_parall.h>
#include <fvm_tesselation.h>
/*----------------------------------------------------------------------------
* Header for the current file
*----------------------------------------------------------------------------*/
#include <fvm_nodal.h>
#include <fvm_nodal_priv.h>
/*----------------------------------------------------------------------------*/
#ifdef __cplusplus
extern "C" {
#if 0
} /* Fake brace to force back Emacs auto-indentation back to column 0 */
#endif
#endif /* __cplusplus */
/*============================================================================
* Static global variables
*============================================================================*/
/* Number of vertices associated with each "nodal" element type */
const int fvm_nodal_n_vertices_element[] = {2, /* Edge */
3, /* Triangle */
4, /* Quadrangle */
0, /* Simple polygon */
4, /* Tetrahedron */
5, /* Pyramid */
6, /* Prism */
8, /* Hexahedron */
0}; /* Simple polyhedron */
/*============================================================================
* Private function definitions
*============================================================================*/
/*----------------------------------------------------------------------------
* Reduction of a nodal mesh representation section: only the associations
* (numberings) necessary to redistribution of fields for output are
* conserved, the full connectivity being no longer useful once it has been
* output.
*
* parameters:
* this_section <-> pointer to structure that should be reduced
*
* returns:
* true if connectivity has been reduced
*----------------------------------------------------------------------------*/
static _Bool
_fvm_nodal_section_reduce(fvm_nodal_section_t * this_section)
{
_Bool retval = false;
/* If we have a tesselation of polyhedra (face index != NULL),
we may need to keep the connectivity information, to
interpolate nodal values to added vertices */
if ( this_section->tesselation == NULL
|| this_section->_face_index == NULL) {
/* Connectivity */
this_section->connectivity_size = 0;
if (this_section->_face_index != NULL)
BFT_FREE(this_section->_face_index);
this_section->face_index = NULL;
if (this_section->_face_num != NULL)
BFT_FREE(this_section->_face_num);
this_section->face_num = NULL;
if (this_section->_vertex_index != NULL)
BFT_FREE(this_section->_vertex_index);
this_section->vertex_index = NULL;
if (this_section->_vertex_num != NULL)
BFT_FREE(this_section->_vertex_num);
this_section->vertex_num = NULL;
retval = true;
}
if (this_section->tesselation != NULL)
fvm_tesselation_reduce(this_section->tesselation);
return retval;
}
/*----------------------------------------------------------------------------
* Change entity parent numbering; this is useful when entities of the
* parent mesh have been renumbered after a nodal mesh representation
* structure's creation. As the parent_num[] array is defined only when
* non trivial (i.e. not 1, 2, ..., n), it may be allocated or freed
* by this function. The return argument corresponds to the new
* pointer which should replace the parent_num input argument.
*
*
* parameters:
* parent_num_size <-- size of local parent numbering array
* new_parent_num <-- pointer to local parent renumbering array
* ({1, ..., n} <-- {1, ..., n})
* parent_num <-> pointer to local parent numbering array
* _parent_num <-> pointer to local parent numbering array if
* owner, NULL otherwise
*
* returns:
* pointer to resulting parent_num[] array
*----------------------------------------------------------------------------*/
static fvm_lnum_t *
_renumber_parent_num(fvm_lnum_t parent_num_size,
const fvm_lnum_t new_parent_num[],
const fvm_lnum_t parent_num[],
fvm_lnum_t _parent_num[])
{
int i;
fvm_lnum_t old_num_id;
fvm_lnum_t *parent_num_p = _parent_num;
_Bool trivial = true;
if (parent_num_size > 0 && new_parent_num != NULL) {
if (parent_num_p != NULL) {
for (i = 0; i < parent_num_size; i++) {
old_num_id = parent_num_p[i] - 1;
parent_num_p[i] = new_parent_num[old_num_id];
if (parent_num_p[i] != i+1)
trivial = false;
}
}
else {
BFT_MALLOC(parent_num_p, parent_num_size, fvm_lnum_t);
if (parent_num != NULL) {
for (i = 0; i < parent_num_size; i++) {
old_num_id = parent_num[i] - 1;
parent_num_p[i] = new_parent_num[old_num_id];
if (parent_num_p[i] != i+1)
trivial = false;
}
}
else {
for (i = 0; i < parent_num_size; i++) {
parent_num_p[i] = new_parent_num[i];
if (parent_num_p[i] != i+1)
trivial = false;
}
}
}
}
if (trivial == true)
BFT_FREE(parent_num_p);
return parent_num_p;
}
/*----------------------------------------------------------------------------
* Renumber vertices based on those actually referenced, and update
* connectivity arrays and parent numbering in accordance.
*
* The number of vertices assigned to the nodal mesh (this_nodal->n_vertices)
* is computed and set by this function. If this number was previously
* non-zero (i.e. vertices have already been assigned to the structure),
* those vertices are considered as referenced. This is useful if we want
* to avoid discarding a given set of vertices, such as when building a
* nodal mesh representation containing only vertices.
*
* parameters:
* this_nodal <-> nodal mesh structure
*----------------------------------------------------------------------------*/
static void
_renumber_vertices(fvm_nodal_t *this_nodal)
{
size_t i;
int section_id;
fvm_lnum_t j;
fvm_lnum_t vertex_id;
fvm_lnum_t n_vertices;
fvm_nodal_section_t *section;
fvm_lnum_t *loc_vertex_num = NULL;
fvm_lnum_t max_vertex_num = 0;
/* Find maximum vertex reference */
/*-------------------------------*/
/* The mesh may already contain direct vertex references
(as in the case of a "mesh" only containing vertices) */
if (this_nodal->n_vertices > 0) {
if (this_nodal->parent_vertex_num != NULL) {
for (j = 0; j < this_nodal->n_vertices; j++) {
if (this_nodal->parent_vertex_num[j] > max_vertex_num)
max_vertex_num = this_nodal->parent_vertex_num[j];
}
}
else
max_vertex_num = this_nodal->n_vertices;
}
/* In most cases, the mesh will reference vertices through elements */
for (section_id = 0; section_id < this_nodal->n_sections; section_id++) {
section = this_nodal->sections[section_id];
if (this_nodal->parent_vertex_num != NULL) {
for (i = 0; i < section->connectivity_size; i++) {
fvm_lnum_t vertex_num
= this_nodal->parent_vertex_num[section->vertex_num[i] - 1];
if (vertex_num > max_vertex_num)
max_vertex_num = vertex_num;
}
}
else {
for (i = 0; i < section->connectivity_size; i++) {
if (section->vertex_num[i] > max_vertex_num)
max_vertex_num = section->vertex_num[i];
}
}
}
/* Flag referenced vertices and compute size */
/*-------------------------------------------*/
BFT_MALLOC(loc_vertex_num, max_vertex_num, fvm_lnum_t);
for (vertex_id = 0; vertex_id < max_vertex_num; vertex_id++)
loc_vertex_num[vertex_id] = 0;
for (section_id = 0; section_id < this_nodal->n_sections; section_id++) {
section = this_nodal->sections[section_id];
if (this_nodal->parent_vertex_num != NULL) {
for (i = 0; i < section->connectivity_size; i++) {
vertex_id
= this_nodal->parent_vertex_num[section->vertex_num[i] - 1] - 1;
if (loc_vertex_num[vertex_id] == 0)
loc_vertex_num[vertex_id] = 1;
}
}
else {
for (i = 0; i < section->connectivity_size; i++) {
vertex_id = section->vertex_num[i] - 1;
if (loc_vertex_num[vertex_id] == 0)
loc_vertex_num[vertex_id] = 1;
}
}
}
/* Build vertices renumbering */
/*----------------------------*/
n_vertices = 0;
for (vertex_id = 0; vertex_id < max_vertex_num; vertex_id++) {
if (loc_vertex_num[vertex_id] == 1) {
n_vertices += 1;
loc_vertex_num[vertex_id] = n_vertices;
}
}
this_nodal->n_vertices = n_vertices;
/* Update connectivity and vertex parent numbering */
/*-------------------------------------------------*/
/* If all vertices are flagged, no need to renumber */
if (n_vertices == max_vertex_num)
BFT_FREE(loc_vertex_num);
else {
/* Update connectivity */
for (section_id = 0; section_id < this_nodal->n_sections; section_id++) {
section = this_nodal->sections[section_id];
if (section->_vertex_num == NULL)
fvm_nodal_section_copy_on_write(section, false, false, false, true);
if (this_nodal->parent_vertex_num != NULL) {
for (i = 0; i < section->connectivity_size; i++) {
vertex_id
= this_nodal->parent_vertex_num[section->vertex_num[i] - 1] - 1;
section->_vertex_num[i] = loc_vertex_num[vertex_id];
}
}
else {
for (i = 0; i < section->connectivity_size; i++) {
vertex_id = section->vertex_num[i] - 1;
section->_vertex_num[i] = loc_vertex_num[vertex_id];
}
}
}
/* Build or update vertex parent numbering */
this_nodal->parent_vertex_num = NULL;
if (this_nodal->_parent_vertex_num != NULL)
BFT_FREE(this_nodal->_parent_vertex_num);
if (loc_vertex_num != NULL) {
BFT_MALLOC(this_nodal->_parent_vertex_num, n_vertices, fvm_lnum_t);
for (vertex_id = 0; vertex_id < max_vertex_num; vertex_id++) {
if (loc_vertex_num[vertex_id] > 0)
this_nodal->_parent_vertex_num[loc_vertex_num[vertex_id] - 1]
= vertex_id + 1;
}
this_nodal->parent_vertex_num = this_nodal->_parent_vertex_num;
}
}
/* Free renumbering array */
BFT_FREE(loc_vertex_num);
}
/*----------------------------------------------------------------------------
* Dump printout of a nodal representation structure section.
*
* parameters:
* this_section <-- pointer to structure that should be dumped
*----------------------------------------------------------------------------*/
static void
_fvm_nodal_section_dump(const fvm_nodal_section_t *this_section)
{
fvm_lnum_t n_elements, i, j;
const fvm_lnum_t *idx, *num;
/* Global indicators */
/*--------------------*/
bft_printf(_("\n"
"Entity dimension: %d\n"
"Number of elements: %ld\n"
"Element type: %s\n"),
this_section->entity_dim, (long)this_section->n_elements,
fvm_elements_type_name[this_section->type]);
bft_printf(_("\n"
"Connectivity_size: %lu\n"
"Stride: %d\n"
"Number of faces: %d\n"),
(unsigned long)(this_section->connectivity_size),
this_section->stride,
(long)(this_section->n_faces));
bft_printf(_("\n"
"Pointers to shareable arrays:\n"
" face_index: %p\n"
" face_num: %p\n"
" vertex_index: %p\n"
" vertex_num: %p\n"
" parent_element_num: %p\n"),
this_section->face_index, this_section->face_num,
this_section->vertex_index, this_section->vertex_num,
this_section->parent_element_num);
bft_printf(_("\n"
"Pointers to local arrays:\n"
" _face_index: %p\n"
" _face_num: %p\n"
" _vertex_index: %p\n"
" _vertex_num: %p\n"
" _parent_element_num: %p\n"),
this_section->_face_index, this_section->_face_num,
this_section->_vertex_index, this_section->_vertex_num,
this_section->_parent_element_num);
if (this_section->face_index != NULL) {
bft_printf(_("\nPolyhedra -> faces (polygons) connectivity:\n\n"));
n_elements = this_section->n_elements;
idx = this_section->face_index;
num = this_section->face_num;
for (i = 0; i < n_elements; i++) {
bft_printf("%10d (idx = %10d) %10d\n",
i+1, idx[i], num[idx[i]]);
for (j = idx[i] + 1; j < idx[i + 1]; j++)
bft_printf(" %10d\n", num[j]);
}
bft_printf(_(" end (idx = %10d)\n"), idx[n_elements]);
}
if (this_section->vertex_index != NULL) {
fvm_lnum_t n_faces = (this_section->n_faces > 0) ?
this_section->n_faces : this_section->n_elements;
bft_printf(_("\nPolygons -> vertices connectivity:\n\n"));
idx = this_section->vertex_index;
num = this_section->vertex_num;
for (i = 0; i < n_faces; i++) {
bft_printf("%10d (idx = %10d) %10d\n",
i + 1, idx[i], num[idx[i]]);
for (j = idx[i] + 1; j < idx[i + 1]; j++)
bft_printf(" %10d\n", num[j]);
}
bft_printf(_(" end (idx = %10d)\n"), idx[n_faces]);
}
else {
bft_printf(_("\nElements -> vertices connectivity:\n\n"));
n_elements = this_section->n_elements;
num = this_section->vertex_num;
switch(this_section->stride) {
case 2:
for (i = 0; i < n_elements; i++)
bft_printf("%10d : %10d %10d\n",
i+1, num[i*2], num[i*2+1]);
break;
case 3:
for (i = 0; i < n_elements; i++)
bft_printf("%10d : %10d %10d %10d\n",
i+1, num[i*3], num[i*3+1], num[i*3+2]);
break;
case 4:
for (i = 0; i < n_elements; i++)
bft_printf("%10d : %10d %10d %10d %10d\n",
i+1, num[i*4], num[i*4+1], num[i*4+2],
num[i*4+3]);
break;
case 5:
for (i = 0; i < n_elements; i++)
bft_printf("%10d : %10d %10d %10d %10d %10d\n",
i+1, num[i*5], num[i*5+1], num[i*5+2],
num[i*5+3], num[i*5+4]);
break;
case 6:
for (i = 0; i < n_elements; i++)
bft_printf("%10d : %10d %10d %10d %10d %10d %10d\n",
i+1, num[i*6], num[i*6+1], num[i*6+2],
num[i*6+3], num[i*6+4], num[i*6+5]);
break;
case 8:
for (i = 0; i < n_elements; i++)
bft_printf("%10d : %10d %10d %10d %10d %10d %10d %10d %10d\n",
i+1, num[i*8], num[i*8+1], num[i*8+2], num[i*8+3],
num[i*8+4], num[i*8+5], num[i*8+6], num[i*8+7]);
break;
default:
for (i = 0; i < n_elements; i++) {
bft_printf("%10d :", i+1);
for (j = 0; j < this_section->stride; j++)
bft_printf(" %10d", num[i*this_section->stride + j]);
bft_printf("\n");
}
}
}
/* Faces tesselation */
if (this_section->tesselation != NULL)
fvm_tesselation_dump(this_section->tesselation);
/* Numbers of associated elements in the parent mesh */
bft_printf(_("\nLocal element numbers in parent mesh:\n"));
if (this_section->parent_element_num == NULL)
bft_printf(_("\n Nil\n\n"));
else {
for (i = 0; i < this_section->n_elements; i++)
bft_printf(" %10d %10d\n", i+1, this_section->parent_element_num[i]);
}
/* Global element numbers (only for parallel execution) */
if (this_section->global_element_num != NULL) {
bft_printf(_("\nGlobal element numbers:\n"));
fvm_io_num_dump(this_section->global_element_num);
}
}
/*============================================================================
* Semi-private function definitions (prototypes in fvm_nodal_priv.h)
*============================================================================*/
/*----------------------------------------------------------------------------
* Creation of a nodal mesh section representation structure.
*
* parameters:
* type <-- type of element defined by this section
*
* returns:
* pointer to created nodal mesh section representation structure
*----------------------------------------------------------------------------*/
fvm_nodal_section_t *
fvm_nodal_section_create(const fvm_element_t type)
{
fvm_nodal_section_t *this_section;
BFT_MALLOC(this_section, 1, fvm_nodal_section_t);
/* Global information */
if (type == FVM_EDGE)
this_section->entity_dim = 1;
else if (type >= FVM_FACE_TRIA && type <= FVM_FACE_POLY)
this_section->entity_dim = 2;
else
this_section->entity_dim = 3;
this_section->n_elements = 0;
this_section->type = type;
/* Connectivity */
this_section->connectivity_size = 0;
if (type != FVM_FACE_POLY && type != FVM_CELL_POLY)
this_section->stride = fvm_nodal_n_vertices_element[type];
else
this_section->stride = 0;
this_section->n_faces = 0;
this_section->face_index = NULL;
this_section->face_num = NULL;
this_section->vertex_index = NULL;
this_section->vertex_num = NULL;
this_section->_face_index = NULL;
this_section->_face_num = NULL;
this_section->_vertex_index = NULL;
this_section->_vertex_num = NULL;
this_section->tesselation = NULL;
/* Numbering */
/*-----------*/
this_section->parent_element_num = NULL;
this_section->global_element_num = NULL;
return (this_section);
}
/*----------------------------------------------------------------------------
* Destruction of a nodal mesh section representation structure.
*
* parameters:
* this_section <-> pointer to structure that should be destroyed
*
* returns:
* NULL pointer
*----------------------------------------------------------------------------*/
fvm_nodal_section_t *
fvm_nodal_section_destroy(fvm_nodal_section_t * this_section)
{
/* Connectivity */
if (this_section->_face_index != NULL)
BFT_FREE(this_section->_face_index);
if (this_section->_face_num != NULL)
BFT_FREE(this_section->_face_num);
if (this_section->_vertex_index != NULL)
BFT_FREE(this_section->_vertex_index);
if (this_section->_vertex_num != NULL)
BFT_FREE(this_section->_vertex_num);
if (this_section->tesselation != NULL)
fvm_tesselation_destroy(this_section->tesselation);
/* Numbering */
/*-----------*/
if (this_section->parent_element_num != NULL) {
this_section->parent_element_num = NULL;
BFT_FREE(this_section->_parent_element_num);
}
if (this_section->global_element_num != NULL)
fvm_io_num_destroy(this_section->global_element_num);
/* Main structure destroyed and NULL returned */
BFT_FREE(this_section);
return (this_section);
}
/*----------------------------------------------------------------------------
* Copy selected shared connectivity information to private connectivity
* for a nodal mesh section.
*
* parameters:
* this_section <-> pointer to section structure
* copy_face_index <-- copy face index (polyhedra only) ?
* copy_face_num <-- copy face numbers (polyhedra only) ?
* copy_vertex_index <-- copy vertex index (polyhedra/polygons only) ?
* copy_vertex_num <-- copy vertex numbers ?
*----------------------------------------------------------------------------*/
void
fvm_nodal_section_copy_on_write(fvm_nodal_section_t *this_section,
_Bool copy_face_index,
_Bool copy_face_num,
_Bool copy_vertex_index,
_Bool copy_vertex_num)
{
fvm_lnum_t n_faces;
size_t i;
if (copy_face_index == true
&& this_section->face_index != NULL && this_section->_face_index == NULL) {
BFT_MALLOC(this_section->face_num, this_section->n_elements + 1, fvm_lnum_t);
for (i = 0; i < (size_t)(this_section->n_elements + 1); i++) {
this_section->_face_index[i] = this_section->face_index[i];
}
this_section->face_index = this_section->_face_index;
}
if (copy_face_num == true
&& this_section->face_num != NULL && this_section->_face_num == NULL) {
n_faces = this_section->face_index[this_section->n_elements];
BFT_MALLOC(this_section->face_num, n_faces, fvm_lnum_t);
for (i = 0; i < (size_t)n_faces; i++) {
this_section->_face_num[i] = this_section->face_num[i];
}
this_section->face_num = this_section->_face_num;
}
if ( copy_vertex_index == true
&& this_section->vertex_index != NULL
&& this_section->_vertex_index == NULL) {
if (this_section->n_faces != 0)
n_faces = this_section->n_faces;
else
n_faces = this_section->n_elements;
BFT_MALLOC(this_section->vertex_index, n_faces + 1, fvm_lnum_t);
for (i = 0; i < (size_t)n_faces + 1; i++) {
this_section->_vertex_index[i] = this_section->vertex_index[i];
}
this_section->vertex_index = this_section->_vertex_index;
}
if (copy_vertex_num == true && this_section->_vertex_num == NULL) {
BFT_MALLOC(this_section->vertex_num,
this_section->connectivity_size, fvm_lnum_t);
for (i = 0; i < this_section->connectivity_size; i++) {
this_section->_vertex_num[i] = this_section->vertex_num[i];
}
this_section->vertex_num = this_section->_vertex_num;
}
}
/*----------------------------------------------------------------------------
* Return global number of elements associated with section.
*
* parameters:
* this_section <-- pointer to section structure
*
* returns:
* global number of elements associated with section
*----------------------------------------------------------------------------*/
fvm_gnum_t
fvm_nodal_section_n_g_elements(const fvm_nodal_section_t *this_section)
{
if (this_section->global_element_num != NULL)
return fvm_io_num_get_global_count(this_section->global_element_num);
else
return this_section->n_elements;
}
/*----------------------------------------------------------------------------
* Return global number of vertices associated with nodal mesh.
*
* parameters:
* this_nodal <-- pointer to nodal mesh structure
*
* returns:
* global number of vertices associated with nodal mesh
*----------------------------------------------------------------------------*/
fvm_gnum_t
fvm_nodal_n_g_vertices(const fvm_nodal_t *this_nodal)
{
fvm_gnum_t n_g_vertices;
if (this_nodal->global_vertex_num != NULL)
n_g_vertices = fvm_io_num_get_global_count(this_nodal->global_vertex_num);
else
n_g_vertices = this_nodal->n_vertices;
return n_g_vertices;
}
/*============================================================================
* Public function definitions
*============================================================================*/
/*----------------------------------------------------------------------------
* Creation of a nodal mesh representation structure.
*
* parameters:
* name <-> name that should be assigned to the nodal mesh
*
* returns:
* pointer to created nodal mesh representation structure
*----------------------------------------------------------------------------*/
fvm_nodal_t *
fvm_nodal_create(const char *name)
{
fvm_nodal_t *this_nodal;
BFT_MALLOC(this_nodal, 1, fvm_nodal_t);
/* Global indicators */
if (name != NULL) {
BFT_MALLOC(this_nodal->name, strlen(name) + 1, char);
strcpy(this_nodal->name, name);
}
else
this_nodal->name = NULL;
this_nodal->dim = 3;
this_nodal->num_dom = fvm_parall_get_rank() + 1;
this_nodal->n_doms = fvm_parall_get_size();
this_nodal->n_sections = 0;
/* Local dimensions */
this_nodal->n_cells = 0;
this_nodal->n_faces = 0;
this_nodal->n_edges = 0;
this_nodal->n_vertices = 0;
/* Local structures */
this_nodal->vertex_coords = NULL;
this_nodal->_vertex_coords = NULL;
this_nodal->parent_vertex_num = NULL;
this_nodal->_parent_vertex_num = NULL;
this_nodal->global_vertex_num = NULL;
this_nodal->sections = NULL;
return (this_nodal);
}
/*----------------------------------------------------------------------------
* Destruction of a nodal mesh representation structure.
*
* parameters:
* this_nodal <-> pointer to structure that should be destroyed
*
* returns:
* NULL pointer
*----------------------------------------------------------------------------*/
fvm_nodal_t *
fvm_nodal_destroy(fvm_nodal_t * this_nodal)
{
int i;
/* Local structures */
if (this_nodal->name != NULL)
BFT_FREE(this_nodal->name);
if (this_nodal->_vertex_coords != NULL)
BFT_FREE(this_nodal->_vertex_coords);
if (this_nodal->parent_vertex_num != NULL) {
this_nodal->parent_vertex_num = NULL;
BFT_FREE(this_nodal->_parent_vertex_num);
}
if (this_nodal->global_vertex_num != NULL)
fvm_io_num_destroy(this_nodal->global_vertex_num);
for (i = 0; i < this_nodal->n_sections; i++)
fvm_nodal_section_destroy(this_nodal->sections[i]);
if (this_nodal->sections != NULL)
BFT_FREE(this_nodal->sections);
/* Main structure destroyed and NULL returned */
BFT_FREE(this_nodal);
return (this_nodal);
}
/*----------------------------------------------------------------------------
* Reduction of a nodal mesh representation structure: only the associations
* (numberings) necessary to redistribution of fields for output are
* conserved, the full connectivity being in many cases no longer useful
* once it has been output. If the del_vertex_num value is set
* to true, vertex-based values may not be output in parallel mode
* after this function is called.
*
* parameters:
* this_nodal <-> pointer to structure that should be reduced
* del_vertex_num <-- indicates if vertex parent indirection and
* I/O numbering are destroyed (1) or not (0)
*----------------------------------------------------------------------------*/
void
fvm_nodal_reduce(fvm_nodal_t *this_nodal,
int del_vertex_num)
{
int i;
_Bool reduce_vertices = true;
/* Connectivity */
for (i = 0; i < this_nodal->n_sections; i++) {
if (_fvm_nodal_section_reduce(this_nodal->sections[i]) == false)
reduce_vertices = false;
}
/* Vertices */
if (reduce_vertices == true) {
if (this_nodal->_vertex_coords != NULL)
BFT_FREE(this_nodal->_vertex_coords);
this_nodal->vertex_coords = NULL;
}
/* Depending on this option, output on vertices may not remain possible */
if (del_vertex_num > 0) {
if (this_nodal->parent_vertex_num != NULL) {
this_nodal->parent_vertex_num = NULL;
BFT_FREE(this_nodal->_parent_vertex_num);
}
if (this_nodal->global_vertex_num != NULL)
this_nodal->global_vertex_num
= fvm_io_num_destroy(this_nodal->global_vertex_num);
}
}
/*----------------------------------------------------------------------------
* Change entity parent numbering; this is useful when entities of the
* parent mesh have been renumbered after a nodal mesh representation
* structure's creation.
*
* parameters:
* this_nodal <-- nodal mesh structure
* new_parent_num <-- pointer to local parent renumbering array
* ({1, ..., n} <-- {1, ..., n})
* entity_dim <-- 3 for cells, 2 for faces, 1 for edges,
* and 0 for vertices
*----------------------------------------------------------------------------*/
void
fvm_nodal_change_parent_num(fvm_nodal_t *this_nodal,
const fvm_lnum_t new_parent_num[],
int entity_dim)
{
/* Vertices */
if (entity_dim == 0) {
this_nodal->_parent_vertex_num
= _renumber_parent_num(this_nodal->n_vertices,
new_parent_num,
this_nodal->parent_vertex_num,
this_nodal->_parent_vertex_num);
this_nodal->parent_vertex_num = this_nodal->_parent_vertex_num;
}
/* Other elements */
else {
int i = 0;
fvm_nodal_section_t *section = NULL;
for (i = 0; i < this_nodal->n_sections; i++) {
section = this_nodal->sections[i];
if (section->entity_dim == entity_dim) {
section->_parent_element_num
= _renumber_parent_num(section->n_elements,
new_parent_num,
section->parent_element_num,
section->_parent_element_num);
section->parent_element_num = section->_parent_element_num;
}
}
}
}
/*----------------------------------------------------------------------------
* Build external numbering for entities based on global numbers.
*
* parameters:
* this_nodal <-- nodal mesh structure
* parent_global_number <-- pointer to list of global (i.e. domain splitting
* independent) parent entity numbers
* entity_dim <-- 3 for cells, 2 for faces, 1 for edges,
* and 0 for vertices
*----------------------------------------------------------------------------*/
void
fvm_nodal_init_io_num(fvm_nodal_t *this_nodal,
const fvm_gnum_t parent_global_numbers[],
int entity_dim)
{
int i;
fvm_nodal_section_t *section;
if (entity_dim == 0)
this_nodal->global_vertex_num
= fvm_io_num_create(this_nodal->parent_vertex_num,
parent_global_numbers,
this_nodal->n_vertices,
0);
else {
for (i = 0; i < this_nodal->n_sections; i++) {
section = this_nodal->sections[i];
if (section->entity_dim == entity_dim) {
section->global_element_num
= fvm_io_num_create(section->parent_element_num,
parent_global_numbers,
section->n_elements,
0);
}
}
}
}
/*----------------------------------------------------------------------------
* Preset number and list of vertices to assign to a nodal mesh.
*
* If the parent_vertex_num argument is NULL, the list is assumed to
* be {1, 2, ..., n}. If parent_vertex_num is given, it specifies a
* list of n vertices from a larger set (1 to n numbering).
*
* Ownership of the given parent vertex numbering array is
* transferred to the nodal mesh representation structure.
*
* This function should be called before fvm_nodal_set_shared_vertices()
* or fvm_nodal_transfer_vertices() if we want to force certain
* vertices to appear in the mesh (especially if we want to define
* a mesh containing only vertices).
*
* parameters:
* this_nodal <-> nodal mesh structure
* n_vertices <-- number of vertices to assign
* parent_vertex_num <-- parent numbers of vertices to assign
*----------------------------------------------------------------------------*/
void
fvm_nodal_define_vertex_list(fvm_nodal_t *this_nodal,
fvm_lnum_t n_vertices,
fvm_lnum_t parent_vertex_num[])
{
assert(this_nodal != NULL);
this_nodal->n_vertices = n_vertices;
this_nodal->parent_vertex_num = NULL;
if (this_nodal->_parent_vertex_num != NULL)
BFT_FREE(this_nodal->_parent_vertex_num);
if (parent_vertex_num != NULL) {
this_nodal->_parent_vertex_num = parent_vertex_num;
this_nodal->parent_vertex_num = parent_vertex_num;
}
}
/*----------------------------------------------------------------------------
* Assign shared vertex coordinates to an extracted nodal mesh,
* renumbering vertex numbers based on those really referenced,
* and updating connectivity arrays in accordance.
*
* This function should only be called once all element sections
* have been added to a nodal mesh representation.
*
* parameters:
* this_nodal <-> nodal mesh structure
* vertex_coords <-- coordinates of parent vertices (interlaced)
*----------------------------------------------------------------------------*/
void
fvm_nodal_set_shared_vertices(fvm_nodal_t *this_nodal,
const fvm_coord_t vertex_coords[])
{
assert(this_nodal != NULL);
/* Map vertex coordinates to array passed as argument
(this_nodal->_vertex_coords remains NULL, so only
the const pointer may be used for a shared array) */
this_nodal->vertex_coords = vertex_coords;
/* If the mesh contains only vertices, its n_vertices and
parent_vertex_num must already have been set, and do not
require updating */
if (this_nodal->n_sections == 0)
return;
/* Renumber vertices based on those really referenced */
_renumber_vertices(this_nodal);
}
/*----------------------------------------------------------------------------
* Assign private vertex coordinates to a nodal mesh,
* renumbering vertex numbers based on those really referenced,
* and updating connectivity arrays in accordance.
*
* Ownership of the given coordinates array is transferred to
* the nodal mesh representation structure.
*
* This function should only be called once all element sections
* have been added to a nodal mesh representation.
*
* parameters:
* this_nodal <-> nodal mesh structure
* vertex_coords <-- coordinates of parent vertices (interlaced)
*
* returns:
* updated pointer to vertex_coords (may be different from initial
* argument if vertices were renumbered).
*----------------------------------------------------------------------------*/
fvm_coord_t *
fvm_nodal_transfer_vertices(fvm_nodal_t *this_nodal,
fvm_coord_t vertex_coords[])
{
fvm_lnum_t i;
int j;
fvm_coord_t *_vertex_coords = vertex_coords;
assert(this_nodal != NULL);
/* Renumber vertices based on those really referenced, and
update connectivity arrays in accordance. */
_renumber_vertices(this_nodal);
/* If renumbering is necessary, update connectivity */
if (this_nodal->parent_vertex_num != NULL) {
int dim = this_nodal->dim;
const fvm_lnum_t *parent_vertex_num = this_nodal->parent_vertex_num;
BFT_MALLOC(_vertex_coords, this_nodal->n_vertices * dim, fvm_coord_t);
for (i = 0; i < this_nodal->n_vertices; i++) {
for (j = 0; j < dim; j++)
_vertex_coords[j*dim + 1]
= vertex_coords[(parent_vertex_num[i]-1)*dim + j];
}
BFT_FREE(vertex_coords);
this_nodal->parent_vertex_num = NULL;
if (this_nodal->_parent_vertex_num != NULL)
BFT_FREE(this_nodal->_parent_vertex_num);
}
this_nodal->_vertex_coords = _vertex_coords;
this_nodal->vertex_coords = _vertex_coords;
return _vertex_coords;
}
/*----------------------------------------------------------------------------
* Obtain the name of a nodal mesh.
*
* parameters:
* this_nodal <-- pointer to nodal mesh structure
*
* returns:
* pointer to constant string containing the mesh name
*----------------------------------------------------------------------------*/
const char *
fvm_nodal_get_name(const fvm_nodal_t *this_nodal)
{
assert(this_nodal != NULL);
return this_nodal->name;
}
/*----------------------------------------------------------------------------
* Return spatial dimension of the nodal mesh.
*
* parameters:
* this_nodal <-- pointer to nodal mesh structure
*
* returns:
* spatial dimension.
*----------------------------------------------------------------------------*/
int
fvm_nodal_get_dim(const fvm_nodal_t *this_nodal)
{
return this_nodal->dim;
}
/*----------------------------------------------------------------------------
* Return maximum dimension of entities in a nodal mesh.
*
* parameters:
* this_nodal <-- pointer to nodal mesh structure
*
* returns:
* maximum dimension of entities in mesh (0 to 3)
*----------------------------------------------------------------------------*/
int
fvm_nodal_get_max_entity_dim(const fvm_nodal_t *this_nodal)
{
int section_id;
int max_entity_dim = 0;
assert(this_nodal != NULL);
for (section_id = 0; section_id < this_nodal->n_sections; section_id++) {
const fvm_nodal_section_t *section = this_nodal->sections[section_id];
if (section->entity_dim > max_entity_dim)
max_entity_dim = section->entity_dim;
}
return max_entity_dim;
}
/*----------------------------------------------------------------------------
* Return number of entities of a given dimension in a nodal mesh.
*
* parameters:
* this_nodal <-- pointer to nodal mesh structure
* entity_dim <-- dimension of entities we want to count (0 to 3)
*
* returns:
* number of entities of given dimension in mesh
*----------------------------------------------------------------------------*/
fvm_lnum_t
fvm_nodal_get_n_entities(const fvm_nodal_t *this_nodal,
int entity_dim)
{
fvm_lnum_t n_entities;
assert(this_nodal != NULL);
switch(entity_dim) {
case 0:
n_entities = this_nodal->n_vertices;
break;
case 1:
n_entities = this_nodal->n_edges;
break;
case 2:
n_entities = this_nodal->n_faces;
break;
case 3:
n_entities = this_nodal->n_cells;
break;
default:
n_entities = 0;
}
return n_entities;
}
/*----------------------------------------------------------------------------
* Return global number of vertices associated with nodal mesh.
*
* parameters:
* this_nodal <-- pointer to nodal mesh structure
*
* returns:
* global number of vertices associated with nodal mesh
*----------------------------------------------------------------------------*/
fvm_gnum_t
fvm_nodal_get_n_g_vertices(const fvm_nodal_t *this_nodal)
{
return fvm_nodal_n_g_vertices(this_nodal);
}
/*----------------------------------------------------------------------------
* Return global number of elements of a given type associated with nodal mesh.
*
* parameters:
* this_nodal <-- pointer to nodal mesh structure
* element_type <-- type of elements for query
*
* returns:
* global number of elements of the given type associated with nodal mesh
*----------------------------------------------------------------------------*/
fvm_gnum_t
fvm_nodal_get_n_g_elements(const fvm_nodal_t *this_nodal,
fvm_element_t element_type)
{
int i;
fvm_gnum_t n_g_elements = 0;
assert(this_nodal != NULL);
for (i = 0; i < this_nodal->n_sections; i++) {
fvm_nodal_section_t *section = this_nodal->sections[i];
if (section->type == element_type)
n_g_elements += fvm_nodal_section_n_g_elements(section);
}
return n_g_elements;
}
/*----------------------------------------------------------------------------
* Return local number of elements of a given type associated with nodal mesh.
*
* parameters:
* this_nodal <-- pointer to nodal mesh structure
* element_type <-- type of elements for query
*
* returns:
* local number of elements of the given type associated with nodal mesh
*----------------------------------------------------------------------------*/
fvm_lnum_t
fvm_nodal_get_n_elements(const fvm_nodal_t *this_nodal,
fvm_element_t element_type)
{
int i;
fvm_lnum_t n_elements = 0;
assert(this_nodal != NULL);
for (i = 0; i < this_nodal->n_sections; i++) {
fvm_nodal_section_t *section = this_nodal->sections[i];
if (section->type == element_type)
n_elements += section->n_elements;
}
return n_elements;
}
/*----------------------------------------------------------------------------
* Return local parent numbering array for all entities of a given
* dimension in a nodal mesh.
*
* The number of entities of the given dimension may be obtained
* through fvm_nodal_get_n_entities(), the parent_num[] array is populated
* with the parent entity numbers of those entities, in order (i.e. in
* local section order, section by section).
*
* parameters:
* this_nodal <-- pointer to nodal mesh structure
* entity_dim <-- dimension of entities we are interested in (0 to 3)
* parent_num --> entity parent numbering (array must be pre-allocated)
*----------------------------------------------------------------------------*/
void
fvm_nodal_get_parent_num(const fvm_nodal_t *this_nodal,
int entity_dim,
fvm_lnum_t parent_num[])
{
int section_id;
fvm_lnum_t i;
fvm_lnum_t entity_count = 0;
assert(this_nodal != NULL);
/* Entity dimension 0: vertices */
if (entity_dim == 0) {
if (this_nodal->parent_vertex_num != NULL) {
for (i = 0; i < this_nodal->n_vertices; i++)
parent_num[entity_count++] = this_nodal->parent_vertex_num[i];
}
else {
for (i = 0; i < this_nodal->n_vertices; i++)
parent_num[entity_count++] = i + 1;
}
}
/* Entity dimension > 0: edges, faces, or cells */
else {
for (section_id = 0; section_id < this_nodal->n_sections; section_id++) {
const fvm_nodal_section_t *section = this_nodal->sections[section_id];
if (section->entity_dim == entity_dim) {
if (section->parent_element_num != NULL) {
for (i = 0; i < section->n_elements; i++)
parent_num[entity_count++] = section->parent_element_num[i];
}
else {
for (i = 0; i < section->n_elements; i++)
parent_num[entity_count++] = i + 1;
}
}
} /* end loop on sections */
}
}
/*----------------------------------------------------------------------------
* Compute tesselation a a nodal mesh's sections of a given type, and add the
* corresponding structure to the mesh representation.
*
* If global element numbers are used (i.e. in parallel mode), this function
* should be only be used after calling fvm_nodal_init_io_num().
*
* If some mesh sections have already been tesselated, their tesselation
* is unchanged.
*
* parameters:
* this_nodal <-> pointer to nodal mesh structure
* type <-> element type that should be tesselated
* error_count --> number of elements with a tesselation error
* counter (optional)
*----------------------------------------------------------------------------*/
void
fvm_nodal_tesselate(fvm_nodal_t *this_nodal,
fvm_element_t type,
fvm_lnum_t *error_count)
{
int section_id;
fvm_lnum_t section_error_count;
assert(this_nodal != NULL);
if (error_count != NULL)
*error_count = 0;
for (section_id = 0; section_id < this_nodal->n_sections; section_id++) {
fvm_nodal_section_t *section = this_nodal->sections[section_id];
if (section->type == type && section->tesselation == NULL) {
section->tesselation = fvm_tesselation_create(type,
section->n_elements,
section->face_index,
section->face_num,
section->vertex_index,
section->vertex_num,
section->global_element_num);
fvm_tesselation_init(section->tesselation,
this_nodal->dim,
this_nodal->vertex_coords,
this_nodal->parent_vertex_num,
§ion_error_count);
if (error_count != NULL)
*error_count += section_error_count;
}
}
}
/*----------------------------------------------------------------------------
* Dump printout of a nodal representation structure.
*
* parameters:
* this_nodal <-- pointer to structure that should be dumped
*----------------------------------------------------------------------------*/
void
fvm_nodal_dump(const fvm_nodal_t *this_nodal)
{
fvm_lnum_t i;
fvm_lnum_t num_vertex = 1;
const fvm_coord_t *coord = this_nodal->vertex_coords;
/* Global indicators */
/*--------------------*/
bft_printf(_("\n"
"Mesh name:\"%s\"\n"),
this_nodal->name);
bft_printf(_("\n"
"Mesh dimension: %d\n"
"Domain number: %d\n"
"Number of domains: %d\n"
"Number of sections: %d\n"),
this_nodal->dim, this_nodal->num_dom, this_nodal->n_doms,
this_nodal->n_sections);
bft_printf(_("\n"
"Number of cells: %d\n"
"Number of faces: %d\n"
"Number of edges: %d\n"
"Number of vertices: %d\n"),
this_nodal->n_cells,
this_nodal->n_faces,
this_nodal->n_edges,
this_nodal->n_vertices);
if (this_nodal->n_vertices > 0) {
bft_printf(_("\n"
"Pointers to shareable arrays:\n"
" vertex_coords: %p\n"
" parent_vertex_num: %p\n"),
this_nodal->vertex_coords,
this_nodal->parent_vertex_num);
bft_printf(_("\n"
"Pointers to local arrays:\n"
" _vertex_coords: %p\n"
" _parent_vertex_num: %p\n"),
this_nodal->_vertex_coords,
this_nodal->_parent_vertex_num);
/* Output coordinates depending on parent numbering */
if (this_nodal->parent_vertex_num == NULL) {
bft_printf(_("\nVertex coordinates:\n\n"));
switch(this_nodal->dim) {
case 1:
for (i = 0; i < this_nodal->n_vertices; i++)
bft_printf("%10d : %12.5f\n",
num_vertex++, (double)(coord[i]));
break;
case 2:
for (i = 0; i < this_nodal->n_vertices; i++)
bft_printf("%10d : %12.5f %12.5f\n",
num_vertex++, (double)(coord[i*2]),
(double)(coord[i*2+1]));
break;
case 3:
for (i = 0; i < this_nodal->n_vertices; i++)
bft_printf("%10d : %12.5f %12.5f %12.5f\n",
num_vertex++, (double)(coord[i*3]),
(double)(coord[i*3+1]), (double)(coord[i*3+2]));
break;
default:
bft_printf("coordinates not output\n"
"dimension = %d unsupported\n", this_nodal->dim);
}
}
else { /* if (this_nodal->parent_vertex_num != NULL) */
bft_printf(_("\nVertex parent and coordinates:\n\n"));
switch(this_nodal->dim) {
case 1:
for (i = 0; i < this_nodal->n_vertices; i++) {
coord = this_nodal->vertex_coords
+ (this_nodal->parent_vertex_num[i]-1);
bft_printf("%10d : %12.5f\n",
num_vertex++, (double)(coord[0]));
}
break;
case 2:
for (i = 0; i < this_nodal->n_vertices; i++) {
coord = this_nodal->vertex_coords
+ ((this_nodal->parent_vertex_num[i]-1)*2);
bft_printf("%10d : %12.5f %12.5f\n",
num_vertex++, (double)(coord[0]), (double)(coord[1]));
}
break;
case 3:
for (i = 0; i < this_nodal->n_vertices; i++) {
coord = this_nodal->vertex_coords
+ ((this_nodal->parent_vertex_num[i]-1)*3);
bft_printf("%10d : %12.5f %12.5f %12.5f\n",
num_vertex++, (double)(coord[0]), (double)(coord[1]),
(double)(coord[2]));
}
break;
default:
bft_printf("coordinates not output\n"
"dimension = %d unsupported\n", this_nodal->dim);
}
}
}
/* Global vertex numbers (only for parallel execution) */
if (this_nodal->global_vertex_num != NULL) {
bft_printf(_("\nGlobal vertex numbers:\n\n"));
fvm_io_num_dump(this_nodal->global_vertex_num);
}
/* Dump element sections */
/*-----------------------*/
for (i = 0; i < this_nodal->n_sections; i++)
_fvm_nodal_section_dump(this_nodal->sections[i]);
}
/*----------------------------------------------------------------------------*/
#ifdef __cplusplus
}
#endif /* __cplusplus */
syntax highlighted by Code2HTML, v. 0.9.1