Skip to content
Snippets Groups Projects
Commit 5a7ddd7a authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

pp

parent 2a58c743
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Geo/GModelIO_CGNS.cpp
+ 118
− 79
View file @ 5a7ddd7a
......@@ -399,13 +399,13 @@ int GModel::writeCGNS(const std::string &name, int zoneDefinition,
if(meshDim == 2) vectorDim = options.vectorDim;
// open the file
int cgIndexFile;
int cgIndexFile=0;
if(cg_open(name.c_str(), MODE_WRITE, &cgIndexFile)) return cgnsErr();
// write the base node
int cgIndexBase;
if(cg_base_write(cgIndexFile, options.baseName.c_str(), meshDim, meshDim,
&cgIndexBase))
int cgIndexBase=0;
if(cg_base_write(cgIndexFile, options.baseName.c_str(), meshDim, meshDim,
&cgIndexBase))
return cgnsErr();
// write information about who generated the mesh
......@@ -415,17 +415,21 @@ int GModel::writeCGNS(const std::string &name, int zoneDefinition,
switch(meshDim) {
case 2:
MZone<2>::preInit();
MZoneBoundary<2>::preInit(); // Add?:
write_CGNS_zones<2>(*this, zoneDefinition, ...);
MZoneBoundary<2>::preInit();
write_CGNS_zones<2>(*this, zoneDefinition, numZone, options,
scalingFactor, vectorDim, groups[face],
cgIndexFile, cgIndexBase);
MZone<2>::postDestroy();
MZoneBoundary<2>::postDestroy(); // Add?:
MZoneBoundary<2>::postDestroy();
break;
case 3:
MZone<3>::preInit();
MZoneBoundary<3>::preInit(); // Add?:
write_CGNS_zones<3>(*this, zoneDefinition, ...);
MZoneBoundary<3>::preInit();
write_CGNS_zones<3>(*this, zoneDefinition, numZone, options,
scalingFactor, vectorDim, groups[region],
cgIndexFile, cgIndexBase);
MZone<3>::postDestroy();
MZoneBoundary<3>::postDestroy(); // Add?:
MZoneBoundary<3>::postDestroy();
break;
}
......@@ -645,6 +649,8 @@ int get_zone_definition(GModel &model, const int zoneDefinition,
int status = 0;
const char *_zoneName = "Partition";
// NBN: pass name to expand_name(): Zone_&N& ==> Zone_Fluid
std::string zn;
//--Get indices for the zonex
......@@ -665,8 +671,9 @@ int get_zone_definition(GModel &model, const int zoneDefinition,
break;
case 2: // Zone defined by physical
partition = -1;
_zoneName = model.getPhysicalName(meshDim, globalPhysicalIt->first)
.c_str();
//_zoneName = model.getPhysicalName(meshDim, globalPhysicalIt->first).c_str();
zn = model.getPhysicalName(meshDim, globalPhysicalIt->first); // NBN:
_zoneName = zn.c_str(); // NBN:
physicalItBegin = globalPhysicalIt++;
physicalItEnd = globalPhysicalIt;
break;
......@@ -809,8 +816,8 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition, const int numZone,
//--Write the zone
// Write the zone node
int cgIndexZone;
int cgZoneSize[3];
int cgIndexZone=0;
int cgZoneSize[3]={0};
cgZoneSize[0] = writeZone->zoneVertVec.size(); // Number of vertices
#ifdef CGNS_TEST1
// Count all the sub-elements in a Triangle 10
......@@ -832,40 +839,49 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition, const int numZone,
if(cg_zone_write(cgIndexFile, cgIndexBase,
writeTask->zoneName.c_str(), cgZoneSize,
Unstructured, &cgIndexZone))
{
return cgnsErr();
}
// Manually maintain the size of the 'zoneInfo vector'. 'push_back'
// is not used because the elements are in order of 'zoneIndex'
if(writeTask->zoneIndex >= zoneInfo.size())
zoneInfo.resize(std::max(2*static_cast<int>(zoneInfo.size()),
writeTask->zoneIndex));
if(writeTask->zoneIndex >= zoneInfo.size()) {
zoneInfo.resize(std::max(2*static_cast<int>(zoneInfo.size()),
writeTask->zoneIndex));
}
zoneInfo[writeTask->zoneIndex].name = writeTask->zoneName;
zoneInfo[writeTask->zoneIndex].cgIndex = cgIndexZone;
// Write the grid node
int cgIndexGrid;
int cgIndexGrid=0;
if(cg_grid_write(cgIndexFile, cgIndexBase, cgIndexZone,
"GridCoordinates", &cgIndexGrid))
return cgnsErr();
// Write the grid coordinates
int cgIndexCoord;
int cgIndexCoord=0;
dBuffer.resize(cgZoneSize[0]);
// x
for(int i = 0; i != cgZoneSize[0]; ++i)
dBuffer[i] = writeZone->zoneVertVec[i]->x()*scalingFactor;
// x-coordinates for this zone
for (int i = 0; i != cgZoneSize[0]; ++i) {
dBuffer[i] = writeZone->zoneVertVec[i]->x()*scalingFactor;
}
if(cg_coord_write(cgIndexFile, cgIndexBase, cgIndexZone, RealDouble,
"CoordinateX", &dBuffer[0], &cgIndexCoord))
return cgnsErr();
// y
for(int i = 0; i != cgZoneSize[0]; ++i)
// y-coordinates for this zone
for(int i = 0; i != cgZoneSize[0]; ++i) {
dBuffer[i] = writeZone->zoneVertVec[i]->y()*scalingFactor;
}
if(cg_coord_write(cgIndexFile, cgIndexBase, cgIndexZone, RealDouble,
"CoordinateY", &dBuffer[0], &cgIndexCoord))
return cgnsErr();
// z
// z-coordinates for this zone
if(vectorDim == 3) {
for(int i = 0; i != cgZoneSize[0]; ++i)
for(int i = 0; i != cgZoneSize[0]; ++i) {
dBuffer[i] = writeZone->zoneVertVec[i]->z()*scalingFactor;
}
if(cg_coord_write(cgIndexFile, cgIndexBase, cgIndexZone, RealDouble,
"CoordinateZ", &dBuffer[0], &cgIndexCoord))
return cgnsErr();
......@@ -951,7 +967,9 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition, const int numZone,
writeZone->zoneElemConn[typeMSHm1].numBoElem + iElemSection,
&writeZone->zoneElemConn[typeMSHm1].connectivity[0],
&cgIndexSection))
{
return cgnsErr();
}
++iElemSection;
}
}
......@@ -963,7 +981,8 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition, const int numZone,
// Write the connectivity for each zone pair
const ZoneConnMap::const_iterator gCEnd = zoneConnMap.end();
for(ZoneConnMap::const_iterator gCIt = zoneConnMap.begin();
gCIt != gCEnd; ++gCIt) {
gCIt != gCEnd; ++gCIt)
{
const int nVert = gCIt->second.vertexPairVec.size();
iBuffer1.resize(nVert);
iBuffer2.resize(nVert);
......@@ -987,7 +1006,9 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition, const int numZone,
&iBuffer1[0], zoneInfo[gCIt->first.zone2].name.c_str(),
Unstructured, PointListDonor, Integer, nVert, &iBuffer2[0],
&cgIndexInterface))
{
return cgnsErr();
}
// In the second zone
if(cg_conn_write
(cgIndexFile, cgIndexBase, zoneInfo[gCIt->first.zone2].cgIndex,
......@@ -995,7 +1016,9 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition, const int numZone,
&iBuffer2[0], zoneInfo[gCIt->first.zone1].name.c_str(),
Unstructured, PointListDonor, Integer, nVert, &iBuffer1[0],
&cgIndexInterface))
{
return cgnsErr();
}
}
//--Pop from the queue
......@@ -1021,7 +1044,8 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition, const int numZone,
if(get_zone_definition
(model, zoneDefinition, numZone, options, DIM, group,
globalZoneIndex, globalPhysicalIt, zoneTask.zoneIndex, partition,
physicalItBegin, physicalItEnd, zoneTask.zoneName)) {
physicalItBegin, physicalItEnd, zoneTask.zoneName))
{
omp_set_lock(&threadWLock);
--threadsWorking;
omp_unset_lock(&threadWLock);
......@@ -1030,7 +1054,8 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition, const int numZone,
zoneTask.zone.clear();
// Loop through all physical in the zone definition
for(PhysGroupMap::const_iterator physicalIt = physicalItBegin;
physicalIt != physicalItEnd; ++physicalIt) {
physicalIt != physicalItEnd; ++physicalIt)
{
// Add elements from all entities in the physical with defined
// partition number
if(partition == -1) {
......@@ -1061,67 +1086,80 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition, const int numZone,
* Write the remaining unconnected vertices as boundary conditions
*--------------------------------------------------------------------*/
if(options.writeBC) {
if(options.writeBC)
{
ZoneBoVec zoneBoVec; // from 'MZoneBoundary.h'
if(mZoneBoundary.exteriorBoundaryVertices
(options.normalSource, zoneBoVec) == 0) {
(options.normalSource, zoneBoVec) == 0)
{
// Sort by zone index and then entity index
const int numBoVert = zoneBoVec.size();
std::vector<int> iZBV(numBoVert);
for(int i = 0; i != numBoVert; ++i) iZBV[i] = i;
std::sort<int*, ZoneBoVecSort>(&iZBV[0], &iZBV[numBoVert-1],
ZoneBoVecSort(zoneBoVec));
dBuffer.reserve(1024);
iBuffer1.reserve(1024);
int iVert = 0;
while(iVert != numBoVert) {
dBuffer.clear();
iBuffer1.clear();
const int zoneIndex = zoneBoVec[iZBV[iVert]].zoneIndex;
const int patchIndex = zoneBoVec[iZBV[iVert]].bcPatchIndex;
const int iVertStart = iVert;
while(iVert != numBoVert &&
zoneBoVec[iZBV[iVert]].zoneIndex == zoneIndex &&
zoneBoVec[iZBV[iVert]].bcPatchIndex == patchIndex) {
VertexBoundary &vertBo = zoneBoVec[iZBV[iVert]];
dBuffer.push_back(vertBo.normal[0]);
dBuffer.push_back(vertBo.normal[1]);
if(vectorDim == 3) dBuffer.push_back(vertBo.normal[2]);
iBuffer1.push_back(vertBo.vertexIndex);
++iVert;
}
const int numBCVert = iVert - iVertStart;
int cgIndexBoco;
std::string patchName = options.patchName;
expand_name(patchName, patchIndex, "Patch");
if(patchName.length() == 0) {
patchName = "Patch_";
patchName += CGNSNameStr(patchIndex+1).c_str();
}
if(cg_boco_write(cgIndexFile, cgIndexBase,
zoneInfo[zoneIndex].cgIndex, patchName.c_str(),
BCTypeNull, PointList, numBCVert, &iBuffer1[0],
&cgIndexBoco))
return cgnsErr();
if(options.normalSource) {
int normalIndex;
if(cg_boco_normal_write(cgIndexFile, cgIndexBase,
zoneInfo[zoneIndex].cgIndex, cgIndexBoco,
&normalIndex, 1, RealDouble, &dBuffer[0]))
if (numBoVert>=1)
{
Msg::Info("writing BoVerts...");
std::vector<int> iZBV(numBoVert);
for(int i = 0; i != numBoVert; ++i) iZBV[i] = i;
std::sort<int*, ZoneBoVecSort>(&iZBV[0], &iZBV[numBoVert-1],
ZoneBoVecSort(zoneBoVec));
dBuffer.reserve(1024);
iBuffer1.reserve(1024);
int iVert = 0;
while(iVert != numBoVert) {
dBuffer.clear();
iBuffer1.clear();
const int zoneIndex = zoneBoVec[iZBV[iVert]].zoneIndex;
const int patchIndex = zoneBoVec[iZBV[iVert]].bcPatchIndex;
const int iVertStart = iVert;
while(iVert != numBoVert &&
zoneBoVec[iZBV[iVert]].zoneIndex == zoneIndex &&
zoneBoVec[iZBV[iVert]].bcPatchIndex == patchIndex)
{
VertexBoundary &vertBo = zoneBoVec[iZBV[iVert]];
dBuffer.push_back(vertBo.normal[0]);
dBuffer.push_back(vertBo.normal[1]);
if(vectorDim == 3) dBuffer.push_back(vertBo.normal[2]);
iBuffer1.push_back(vertBo.vertexIndex);
++iVert;
}
const int numBCVert = iVert - iVertStart;
int cgIndexBoco;
std::string patchName = options.patchName;
expand_name(patchName, patchIndex, "Patch");
if(patchName.length() == 0) {
patchName = "Patch_";
patchName += CGNSNameStr(patchIndex+1).c_str();
}
if(cg_boco_write(cgIndexFile, cgIndexBase,
zoneInfo[zoneIndex].cgIndex, patchName.c_str(),
BCTypeNull, PointList, numBCVert, &iBuffer1[0],
&cgIndexBoco))
{
return cgnsErr();
}
if(options.normalSource) {
int normalIndex;
if(cg_boco_normal_write(cgIndexFile, cgIndexBase,
zoneInfo[zoneIndex].cgIndex, cgIndexBoco,
&normalIndex, 1, RealDouble, &dBuffer[0]))
{
return cgnsErr();
}
}
}
}
}
}
// std::cout << "Leaving master thread\n"; // DBG
} // End master thread instructions
// NBN: MZoneBoundary defines no destructor, so force
// the deallocation of the new pointers with clear()
mZoneBoundary.clear();
Msg::Info("Leaving master thread"); // DBG
} // End master thread instructions
} // End omp parallel section
//--Destroy omp locks
......@@ -1129,6 +1167,7 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition, const int numZone,
omp_destroy_lock(&threadWLock);
omp_destroy_lock(&queueLock);
return 0;
}
#else
......
Geo/MZoneBoundary.cpp
+ 66
− 19
View file @ 5a7ddd7a
......@@ -9,6 +9,8 @@
#if defined(HAVE_LIBCGNS)
#define DEBUG_FaceT 0 // NBN: toggle reporting of face insertions
#include <iostream> // DBG
#include <limits> // ?
......@@ -217,7 +219,7 @@ int edge_normal
&faces, SVector3 &boNormal, const int onlyFace = -1)
{
double par;
double par=0.0;
// Note: const_cast used to match MVertex.cpp interface
if(!reparamMeshVertexOnEdge(vertex, gEdge, par)) return 1;
......@@ -580,7 +582,7 @@ void updateBoVec<3, MFace>
// GEdge, and then the MElement is still on the GFace. There is
// also the unlikely case where the two other MVertex are both on
// edges ... and the MElement is still on the GFace.
if(!matchedFace && nVOnF == 3) {
if(!matchedFace && (3 == nVOnF)) {
const MVertex *vertex2;
const MVertex *vertex3;
switch(vertexOnF) {
......@@ -811,25 +813,54 @@ int MZoneBoundary<DIM>::interiorBoundaryVertices
//--A new vertex was inserted
// if(debug) {
// std::cout << "This vertex is new and has bo faces:\n";
// }
// Copy faces
const int nFace = zoneVertData.faces.size();
int iCount=0;
for(int iFace = 0; iFace != nFace; ++iFace) {
// if(debug) {
// std::cout << " ("
// << zoneVertData.faces[iFace]->first.getVertex(0)->x() << ","
// << zoneVertData.faces[iFace]->first.getVertex(0)->y()
// << "), ("
// << zoneVertData.faces[iFace]->first.getVertex(1)->x() << ","
// << zoneVertData.faces[iFace]->first.getVertex(1)->y()
// << ")\n";
// }
#if (0)
// NBN: changing the member object for a pointer means
// we need to do the following in two steps (see below)
globalVertData.faces.push_back
(typename GlobalVertexData<FaceT>::FaceDataB
(newZoneIndex, zoneVertData.faces[iFace]));
#else
// Using FaceAllocator<T> with face2Pool, the std constructors
// are not called, so FaceDataB std members are incomplete.
// By replacing the FaceT member with FaceT* we fix the issue of
// auto-deleting "un-constructed" containers. Here, the simple
// data type (pointer) member is allocated in the ctor, then we
// create and store its FaceT object once the FaceDataB object
// is safely in the container.
//
// Note: we must now make sure to delete these pointers.
// See adjusted version of MZoneBoundary::clear();
// Step 1: append new FaceDataB<> object
typename GlobalVertexData<FaceT>::FaceDataB& tFDB =
globalVertData.faces.push_back
(typename GlobalVertexData<FaceT>::FaceDataB
(newZoneIndex, zoneVertData.faces[iFace]) );
// Step 2: construct its internal face object
tFDB.face = new FaceT(zoneVertData.faces[iFace]->first);
#if (DEBUG_FaceT)
if (! tFDB.face) {
Msg::Info("MZoneBoundary<DIM>::interiorBoundaryVertices, failed to alloc face object");
} else {
int nv = tFDB.face->getNumVertices();
Msg::Info("interiorBoundaryVertices: allocated FaceT object %5d with %d verts", ++iCount, nv);
}
#endif
#endif
}
// Copy information about the vertex in the zone
globalVertData.zoneData.push_back
(typename GlobalVertexData<FaceT>::ZoneData
......@@ -871,8 +902,11 @@ int MZoneBoundary<DIM>::interiorBoundaryVertices
&zoneVertData.faces[0];
for(int nZFace = zoneVertData.faces.size(); nZFace--;) {
bool foundMatch = false;
for(int iGFace = 0; iGFace != nGFace; ++iGFace) {
if((*zFace)->first == globalVertData.faces[iGFace].face) {
for(int iGFace = 0; iGFace != nGFace; ++iGFace)
{
// NBN: face is now a pointer, so need to de-reference
//if((*zFace)->first == globalVertData.faces[iGFace].face )
if((*zFace)->first == *(globalVertData.faces[iGFace].face)) {
foundMatch = true;
// Faces match - delete from 'globalVertData'.
globalVertData.faces.erase(iGFace);
......@@ -892,8 +926,9 @@ int MZoneBoundary<DIM>::interiorBoundaryVertices
// If there are no more faces, connectivity for this vertex is complete
// and it may be deleted.
if(globalVertData.faces.size() == 0)
if(globalVertData.faces.size() == 0) {
globalBoVertMap.erase(insGblBoVertMap.first);
}
else {
// Update the list of zones attached to this vertex
globalVertData.zoneData.push_back
......@@ -1012,13 +1047,25 @@ int MZoneBoundary<DIM>::exteriorBoundaryVertices
template<>
template<>
MZoneBoundary<2>::GlobalVertexData<MEdge>::FaceDataB::FaceDataB()
: face(0, 0)
:
//face(0, 0), // NBN: replaced this MEdge object
face(NULL), // NBN: with a pointer to MEdge
parentElement(0), // NBN: also, init members
parentFace(0),
faceIndex(0),
zoneIndex(0)
{ }
template<>
template<>
MZoneBoundary<3>::GlobalVertexData<MFace>::FaceDataB::FaceDataB()
: face(0, 0, 0)
:
//face(0, 0, 0), // NBN: replaced this MFace object
face(NULL), // NBN: with a pointer to MFace
parentElement(0), // NBN: also, init members
parentFace(0),
faceIndex(0),
zoneIndex(0)
{ }
......
Geo/MZoneBoundary.h
+ 41
− 4
View file @ 5a7ddd7a
......@@ -75,7 +75,7 @@ struct ZoneConnectivity
int vertexIndex1;
int vertexIndex2;
// Constructors
VertexPair()
VertexPair()
: vertex(0), vertexIndex1(0), vertexIndex2(0)
{ }
VertexPair(MVertex *const _vertex,
......@@ -242,14 +242,28 @@ class MZoneBoundary
{
struct FaceDataB
{
FaceT face;
// NBN: cannot use a FaceT object in FaceVector.
// class FaceT has embedded std::vector objects;
// custom allocator for FaceVector<T> does not call ctors,
// but std:: dtors will try to delete _v, _si
//
// Simple fix: use a pointer to FaceT, then build the
// FaceT object once the FaceDataB structure has been
// safely added to the container (push_back)
//FaceT face; // NBN: FaceT contains std:: containers
FaceT* face; // NBN: use FaceT* (then init in two steps)
MElement *parentElement;
int parentFace;
int faceIndex;
int zoneIndex;
FaceDataB(const int _zoneIndex,
const typename MZone<DIM>::BoFaceMap::const_iterator bFMapIt)
: face(bFMapIt->first), parentElement(bFMapIt->second.parentElement),
:
//face(bFMapIt->first),
face(NULL), // NBN: need to load this after insertion into container
parentElement(bFMapIt->second.parentElement),
parentFace(bFMapIt->second.parentFace),
faceIndex(bFMapIt->second.faceIndex), zoneIndex(_zoneIndex)
{ }
......@@ -271,7 +285,9 @@ class MZoneBoundary
// fails on some compilers (earlier versions of g++?)
// The default constructor is required by 'set_offsets()' in
// class 'FaceAllocator'. This is invoked by preInit() below.
ZoneData() { };
ZoneData()
: vertexIndex(0), zoneIndex(0) // NBN: init members
{ }
friend class CCon::FaceAllocator<ZoneData>;
};
CCon::FaceVector<FaceDataB> faces;
......@@ -299,6 +315,27 @@ class MZoneBoundary
void clear()
{
// NBN: using FaceT* so need to dealloc:
int icount = 0;
GlobalBoVertexMap::iterator itEnd = globalBoVertMap.end();
for (GlobalBoVertexMap::iterator itBoV = globalBoVertMap.begin();
itBoV != itEnd; ++itBoV)
{
// ... clear the faces
typename GlobalVertexData<FaceT>& ref = itBoV->second;
size_t nf = ref.faces.size();
for (int i=0; i<nf; ++i) {
++ icount;
FaceT* p = ref.faces[i].face;
if (p) {
delete (p);
p = NULL;
}
}
}
Msg::Info("cleared %d faces.", icount);
// finally, clear the container
globalBoVertMap.clear();
}
......
Plugin/MathEval.h
+ 1
− 1
View file @ 5a7ddd7a
......@@ -7,7 +7,7 @@
#define _MATH_EVAL_H_
#include "Plugin.h"
;
extern "C"
{
GMSH_Plugin *GMSH_RegisterMathEvalPlugin();
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment