Forked from
gmsh / gmsh
16337 commits behind the upstream repository.
-
Christophe Geuzaine authoredChristophe Geuzaine authored
GModelIO_CGNS.cpp 18.98 KiB
// Gmsh - Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
// Copyright (C) 2006 S. Guzik, C. Geuzaine, J.-F. Remacle
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program 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 General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA.
//
// Please report all bugs and problems to <gmsh@geuz.org>.
#include <string.h>
#include <iostream> // DBG
#include <list>
#include <map>
#include <set>
#include <string>
#include <vector>
#include "GModel.h"
#include "Message.h"
#include "MElement.h"
#include "MNeighbour.h"
#include "MVertex.h"
#if defined(HAVE_LIBCGNS)
#include <cgnslib.h>
int cgnsErr(const int cgIndexFile = -1)
{
Msg::Error(cg_get_error());
if(cgIndexFile != -1)
if(cg_close(cgIndexFile)) Msg::Error("Even unable to close CGNS file");
return 0;
}
/*==============================================================================
* File scope structures used for CGNS I/O
*============================================================================*/
//--Structure describing the zones sharing a vertex or element
struct ZoneConn_t {
unsigned int indexZone; // Index of the zone
unsigned int indexConn; // Index of the vertex or element used
// to describe the connectivity
ZoneConn_t(const unsigned int iz, const unsigned int ic) :
indexZone(iz), indexConn(ic) { }
};
//--Class to make C style CGNS name strings act like C++ types
class CgnsNameStr {
private:
char name[33];
public: // Constructors
CgnsNameStr() { name[0] = '\0'; }
CgnsNameStr(const char *const cstr) {
std::strncpy(name, cstr, 32);
name[32] = '\0';
}
CgnsNameStr(const CgnsNameStr& cgs) { std::strcpy(name, cgs.name); }
// Assignment
CgnsNameStr& operator=(const CgnsNameStr& cgs) {
if ( &cgs != this ) {
std::strcpy(name, cgs.name);
}
return *this;
}
CgnsNameStr& operator=(const char *const cstr) {
std::strncpy(name, cstr, 32);
name[32] = '\0';
return *this;
}
// Return the C string
char *c_str() { return name; }
const char *c_str() const { return name; }
};
/*==============================================================================
*
* Routine readCGNS
*
* Purpose
* =======
*
*
*
* I/O
* ===
*
* name - (I) file name
*
*============================================================================*/
int GModel::readCGNS(const std::string &name)
{
return 1;
}
/*==============================================================================
*
* Routine writeCGNS
*
* Purpose
* =======
*
* Writes out mesh in CGNS format
*
* I/O
* ===
*
* name - (I) file name
* scalingFactor - (I) scaling for coordinates
*
*============================================================================*/
int GModel::writeCGNS(const std::string &name, double scalingFactor)
{
//----- Type definitions -----//
typedef std::map<int, std::vector<GEntity*> >::iterator zone_iterator; typedef std::vector<MTriangle*>::const_iterator const_triangle_iterator;
typedef std::vector<MQuadrangle*>::const_iterator const_quadrangle_iterator;
typedef std::list<MVertex*> BoundaryVertList;
typedef std::set<MVertex*> InteriorVertSet;
typedef std::list<MElement*> BoundaryElemList;
typedef std::set<MElement*> InteriorElemSet;
//----- Parameters -----//
enum {
vertex,
edge,
face,
region
};
int meshDim = 0; // Dimension of the mesh
std::map<int, std::vector<GEntity*> > groups[4];
// vector of entities that belong to
// each physical group (in each
// dimension)
std::multimap<MVertex*, ZoneConn_t> sharedVertex;
// map of shared vertices describing the
// zones that share it and the index of
// the vertex in each of the zones.
// std::multimap<MEdge, ZoneConn_t> sharedEdge;
// std::multimap<MFace, ZoneConn_t> sharedFace;
MNeighbour meshNeighbours; // Database of neighbours (computed for
// each zone)
//--Get a list of groups in each dimension and each of the entities in that
//--group. The groups will define the zones. If no 2D or 3D groups exist, just
//--store all mesh elements in one zone.
getPhysicalGroups(groups);
if(groups[face].size() + groups[region].size() == 0) {
// If there are no 2D or 3D physical groups, save all elements in one zone.
// Pretend that there is one physical group ecompassing all the elements.
for(riter it = firstRegion(); it != lastRegion(); ++it)
if((*it)->tetrahedra.size() + (*it)->hexahedra.size() +
(*it)->prisms.size() + (*it)->pyramids.size() > 0)
groups[region][1].push_back(*it);
if(!groups[region].size()) {
// No 3D elements were found. Look for 2D elements.
for(fiter it = firstFace(); it != lastFace(); ++it)
if((*it)->triangles.size() + (*it)->quadrangles.size() > 0)
groups[face][1].push_back(*it);
if(!groups[face].size()) {
Msg::Error("No mesh elements were found");
return 0;
}
}
}
//--Open the file and get ready to write the zones
// Will this be a 2D or 3D mesh?
if(groups[region].size()) meshDim = 3;
else meshDim = 2;
// open the file
int cgIndexFile;
if(cg_open(name.c_str(), MODE_WRITE, &cgIndexFile)) return cgnsErr();
// write the base node
int cgIndexBase;
if(cg_base_write(cgIndexFile, "Base_0000", 2, 2, &cgIndexBase))
return cgnsErr();
// write information about who generated the mesh if(cg_goto(cgIndexFile, cgIndexBase, "end")) return(cgnsErr());
if(cg_descriptor_write("About", "Created by Gmsh")) return cgnsErr();
unsigned int indexZone = 0; // Index of a zone
const unsigned int nZone = groups[face].size() + groups[region].size();
// Total number of zones
std::vector<CgnsNameStr> zoneName(nZone);
// zone names
std::vector<BoundaryVertList> zoneBoundaryVert(nZone);
/*--------------------------------------------------------------------*
* Write the 3D zone information
*--------------------------------------------------------------------*/
// TODO: add
/*--------------------------------------------------------------------*
* Write the 2D zone information
*--------------------------------------------------------------------*/
for(zone_iterator itZone = groups[face].begin(); itZone != groups[face].end();
++itZone) {
const int nFace = itZone->second.size();
std::cout << "Working on zone " << indexZone+1 << " " << itZone->first
<< std::endl;
std::sprintf(zoneName[indexZone].c_str(), "Zone_%d\0", indexZone+1);
//--Compute the mesh neighbours for this group
meshNeighbours.add_elements_in_entities<GFace>(itZone->second.begin(),
itZone->second.end());
//--Load interior lists with all vertices and all elements
unsigned int numTri3 = 0;
unsigned int numTri6 = 0;
unsigned int numQuad4 = 0;
unsigned int numQuad8 = 0;
unsigned int numQuad9 = 0;
InteriorVertSet interiorVert;
BoundaryElemList boundaryElem;
InteriorElemSet interiorElem;
for(unsigned int iFace = 0; iFace != nFace; ++iFace) {
GFace* face = static_cast<GFace*>(itZone->second[iFace]);
const std::vector<MTriangle*>::const_iterator itEndTri =
face->triangles.end();
for(std::vector<MTriangle*>::const_iterator itTri =
face->triangles.begin(); itTri != itEndTri ; ++itTri) {
interiorElem.insert(*itTri);
interiorVert.insert((*itTri)->getVertex(0));
interiorVert.insert((*itTri)->getVertex(1));
interiorVert.insert((*itTri)->getVertex(2));
switch ( (*itTri)->getTypeForMSH() ) {
case MSH_TRI_3:
++numTri3;
break;
case MSH_TRI_6:
++numTri6;
interiorVert.insert((*itTri)->getVertex(3));
interiorVert.insert((*itTri)->getVertex(4));
interiorVert.insert((*itTri)->getVertex(5));
break;
}
}
const std::vector<MQuadrangle*>::const_iterator itEndQuad =
face->quadrangles.end();
for(std::vector<MQuadrangle*>::const_iterator itQuad =
face->quadrangles.begin(); itQuad != itEndQuad ; ++itQuad) { interiorElem.insert(*itQuad);
interiorVert.insert((*itQuad)->getVertex(0));
interiorVert.insert((*itQuad)->getVertex(1));
interiorVert.insert((*itQuad)->getVertex(2));
interiorVert.insert((*itQuad)->getVertex(3));
switch ( (*itQuad)->getTypeForMSH() ) {
case MSH_QUA_4:
++numQuad4;
break;
case MSH_QUA_8:
++numQuad8;
interiorVert.insert((*itQuad)->getVertex(4));
interiorVert.insert((*itQuad)->getVertex(5));
interiorVert.insert((*itQuad)->getVertex(6));
interiorVert.insert((*itQuad)->getVertex(7));
break;
case MSH_QUA_9:
++numQuad9;
interiorVert.insert((*itQuad)->getVertex(4));
interiorVert.insert((*itQuad)->getVertex(5));
interiorVert.insert((*itQuad)->getVertex(6));
interiorVert.insert((*itQuad)->getVertex(7));
interiorVert.insert((*itQuad)->getVertex(8));
break;
}
}
}
//--Check edges for zone boundaries and move vertices and elements from interior
//--sets to boundary lists
const MNeighbour::Edge_const_iterator itEndEdge =
meshNeighbours.edge_end();
for(MNeighbour::Edge_const_iterator itEdge = meshNeighbours.edge_begin();
itEdge != itEndEdge; ++itEdge) {
if(itEdge.num_neighbours() == 1) { // Boundary edge
// Find the boundary element
MElement *element;
meshNeighbours.edge_neighbours(itEdge, &element);
// Move primary vertices from interior to boundary
InteriorVertSet::iterator itVert =
interiorVert.find(itEdge->getVertex(0));
if(itVert != interiorVert.end()) {
zoneBoundaryVert[indexZone].push_back(*itVert);
interiorVert.erase(itVert);
}
itVert = interiorVert.find(itEdge->getVertex(1));
if(itVert != interiorVert.end()) {
zoneBoundaryVert[indexZone].push_back(*itVert);
interiorVert.erase(itVert);
}
// If this is a higher-order element, find the edge on the boundary
// and then the mid-edge vertex
if(element->getNumVertices() != element->getNumPrimaryVertices()) {
// Second-order element
int iEPE = 0;
while(element->getEdge(iEPE) != *itEdge) ++iEPE;
// Add mid-vertex on edge iEPE
itVert = interiorVert.find
(element->getVertex(iEPE + element->getNumPrimaryVertices()));
if(itVert != interiorVert.end()) {
zoneBoundaryVert[indexZone].push_back(*itVert);
interiorVert.erase(itVert);
}
}
}
}
//--Loop through all the boundary vertices
int numVert = 1;
BoundaryVertList::const_iterator itEndBVert =
zoneBoundaryVert[indexZone].end();
{
MElement **boElem = new MElement*[meshNeighbours.max_vertex_neighbours()];
for(BoundaryVertList::iterator itVert =
zoneBoundaryVert[indexZone].begin(); itVert != itEndBVert;
++itVert) {
// Any elements with a vertex in the boundary vertex list should be
// moved to the boundary element list.
const int nElem = meshNeighbours.vertex_neighbours(*itVert, boElem);
for(int iElem = 0; iElem != nElem; ++iElem) {
InteriorElemSet::iterator itElem = interiorElem.find(boElem[iElem]);
if(itElem != interiorElem.end()) {
boundaryElem.push_back(boElem[iElem]);
interiorElem.erase(itElem);
}
}
// Renumber the boundary vertices
(*itVert)->setNum(numVert);
// Add all boundary vertices to the map of shared vertices
sharedVertex.insert(std::make_pair(*itVert, ZoneConn_t(indexZone,
numVert)));
++numVert;
}
delete[] boElem;
}
//--Renumber the interior mesh vertices
InteriorVertSet::const_iterator itEndIVert = interiorVert.end();
for(InteriorVertSet::iterator itVert = interiorVert.begin();
itVert != itEndIVert; ++itVert) (*itVert)->setNum(numVert++);
//--Write the zone node
int cgIndexZone;
int cgZoneSize[3];
cgZoneSize[0] = zoneBoundaryVert[indexZone].size() + interiorVert.size();
cgZoneSize[1] = boundaryElem.size() + interiorElem.size();
cgZoneSize[2] = zoneBoundaryVert[indexZone].size();
if(cg_zone_write(cgIndexFile, cgIndexBase, zoneName[indexZone].c_str(),
cgZoneSize, Unstructured, &cgIndexZone)) return cgnsErr();
//--Write the grid node
int cgIndexGrid;
if(cg_grid_write(cgIndexFile, cgIndexBase, cgIndexZone, "GridCoordinates",
&cgIndexGrid)) return cgnsErr();
//--Write the grid coordinates
{
std::vector<double> coordx(cgZoneSize[0]);
std::vector<double> coordy(cgZoneSize[1]);
int iVert = 0;
for(BoundaryVertList::const_iterator itVert =
zoneBoundaryVert[indexZone].begin(); itVert != itEndBVert;
++itVert) {
coordx[iVert] = (*itVert)->x()*scalingFactor;
coordy[iVert] = (*itVert)->y()*scalingFactor;
++iVert;
}
for(InteriorVertSet::const_iterator itVert = interiorVert.begin();
itVert != itEndIVert; ++itVert) {
coordx[iVert] = (*itVert)->x()*scalingFactor;
coordy[iVert] = (*itVert)->y()*scalingFactor;
++iVert;
}
int cgIndexCoord; if(cg_coord_write(cgIndexFile, cgIndexBase, cgIndexZone, RealDouble,
"CoordinateX", &coordx[0], &cgIndexCoord))
return cgnsErr();
if(cg_coord_write(cgIndexFile, cgIndexBase, cgIndexZone, RealDouble,
"CoordinateY", &coordy[0], &cgIndexCoord))
return cgnsErr();
}
//--Write out the element connectivity
{
int mixedElem = -1;
int numVPE;
int numElemData = 0;
ElementType_t cgElemType;
if(numTri3) {
++mixedElem;
numVPE = 3;
numElemData += numTri3*3;
cgElemType = TRI_3;
}
if(numTri6) {
++mixedElem;
numVPE = 6;
numElemData += numTri6*6;
cgElemType = TRI_6;
}
if(numQuad4) {
++mixedElem;
numVPE = 4;
numElemData += numQuad4*4;
cgElemType = QUAD_4;
}
if(numQuad8) {
++mixedElem;
numVPE = 8;
numElemData += numQuad8*8;
cgElemType = QUAD_8;
}
if(numQuad9) {
++mixedElem;
numVPE = 9;
numElemData += numQuad9*9;
cgElemType = QUAD_9;
}
if(mixedElem) {
numVPE = 0;
numElemData += cgZoneSize[1];
cgElemType = MIXED;
}
std::vector<int> elementConn(numElemData);
int iElemData = 0;
// Load boundary elements
BoundaryElemList::const_iterator itEndBElem = boundaryElem.end();
if(mixedElem) {
for(BoundaryElemList::const_iterator itElem = boundaryElem.begin();
itElem != itEndBElem; ++itElem) {
switch((*itElem)->getTypeForMSH()) {
case MSH_TRI_3:
elementConn[iElemData++] = TRI_3;
break;
case MSH_TRI_6:
elementConn[iElemData++] = TRI_6;
break;
case MSH_QUA_4:
elementConn[iElemData++] = QUAD_4;
break;
case MSH_QUA_8: elementConn[iElemData++] = QUAD_8;
break;
case MSH_QUA_9:
elementConn[iElemData++] = QUAD_9;
break;
}
const int nVPE = (*itElem)->getNumVertices();
for(int iVPE = 0; iVPE != nVPE; ++iVPE)
elementConn[iElemData++] = (*itElem)->getVertex(iVPE)->getNum();
}
}
else {
const int nVPE = numVPE;
for(BoundaryElemList::const_iterator itElem = boundaryElem.begin();
itElem != itEndBElem; ++itElem)
for(int iVPE = 0; iVPE != nVPE; ++iVPE)
elementConn[iElemData++] = (*itElem)->getVertex(iVPE)->getNum();
}
// Load interior elements
InteriorElemSet::const_iterator itEndIElem = interiorElem.end();
if(mixedElem) {
for(InteriorElemSet::const_iterator itElem = interiorElem.begin();
itElem != itEndIElem; ++itElem) {
switch((*itElem)->getTypeForMSH()) {
case MSH_TRI_3:
elementConn[iElemData++] = TRI_3;
break;
case MSH_TRI_6:
elementConn[iElemData++] = TRI_6;
break;
case MSH_QUA_4:
elementConn[iElemData++] = QUAD_4;
break;
case MSH_QUA_8:
elementConn[iElemData++] = QUAD_8;
break;
case MSH_QUA_9:
elementConn[iElemData++] = QUAD_9;
break;
}
const int nVPE = (*itElem)->getNumVertices();
for(int iVPE = 0; iVPE != nVPE; ++iVPE)
elementConn[iElemData++] = (*itElem)->getVertex(iVPE)->getNum();
}
}
else {
const int nVPE = numVPE;
for(InteriorElemSet::const_iterator itElem = interiorElem.begin();
itElem != itEndIElem; ++itElem)
for(int iVPE = 0; iVPE != nVPE; ++iVPE)
elementConn[iElemData++] = (*itElem)->getVertex(iVPE)->getNum();
}
// Write to CGNS file
int cgIndexSection;
if(cg_section_write(cgIndexFile, cgIndexBase, cgIndexZone, "Section_Main",
cgElemType, 1, cgZoneSize[1], boundaryElem.size(),
&elementConn[0], &cgIndexSection)) return cgnsErr();
}
meshNeighbours.clear();
++indexZone;
}
//--Close the file
if(cg_close(cgIndexFile)) return cgnsErr();
std::cout << "Done\n";
return 1;
}
#else
int GModel::readCGNS(const std::string &name)
{
Msg::Error("This version of Gmsh was compiled without CGNS support");
return 0;
}
int GModel::writeCGNS(const std::string &name, double scalingFactor)
{
Msg::Error("This version of Gmsh was compiled without CGNS support");
return 0;
}
#endif