diff --git a/Geo/CGNSOptions.h b/Geo/CGNSOptions.h
index fe32ec018ca6dac92928cef22947eefc69254981..3052f8d295d263aeb8b58c78e2b1359fdfbef7c6 100644
--- a/Geo/CGNSOptions.h
+++ b/Geo/CGNSOptions.h
@@ -26,6 +26,8 @@ struct CGNSOptions
// Type of CGNS connectivity description
std::string baseName;
std::string zoneName;
+ std::string interfaceName;
+ std::string patchName;
// Default values
CGNSOptions() :
@@ -33,7 +35,9 @@ struct CGNSOptions
gridConnectivityLocation(2), // Assumes Vertex = 2 in cgnslib.h
connectivityNodeType(Generalized),
baseName("Base_0"),
- zoneName("Zone_&I0&")
+ zoneName("Zone_&I0&"),
+ interfaceName("Interface_&I0&"),
+ patchName("Patch_&I&")
{ }
};
diff --git a/Geo/GModelIO_CGNS.cpp b/Geo/GModelIO_CGNS.cpp
index 2670f8392a176a019728a0d644906dd34be1b758..cc69b722e969ad99900d870ad1087776c5a55802 100644
--- a/Geo/GModelIO_CGNS.cpp
+++ b/Geo/GModelIO_CGNS.cpp
@@ -39,6 +39,7 @@
int cgnsErr(const int cgIndexFile = -1)
{
+ Message::Error("This is a CGNS error\n");
Message::Error(cg_get_error());
if(cgIndexFile != -1)
if(cg_close(cgIndexFile)) Message::Error("Unable to close CGNS file");
@@ -54,6 +55,32 @@ typedef std::map<int, std::vector<GEntity*> > PhysGroupMap;
// Type providing a vector of entities
// for a physical group index
+//--Class to gather elements that belong to a partition. This class mimics a
+//--GEntity class.
+
+class DummyPartitionEntity : public GEntity
+{
+ public:
+ DummyPartitionEntity()
+ : GEntity(0, 0)
+ { }
+
+ // number of types of elements
+ int getNumElementTypes() const { return 1; }
+ void getNumMeshElements(unsigned *const c) const
+ {
+ c[0] += elements.size();
+ }
+
+ // get the start of the array of a type of element
+ MElement *const *getStartElementType(int type) const
+ {
+ return &elements[0];
+ }
+
+ std::vector<MElement*> elements;
+};
+
//--Class to make C style CGNS name strings act like C++ types
class CGNSNameStr
@@ -169,10 +196,9 @@ struct ElemSortCGNSType
*============================================================================*/
template <unsigned DIM>
-int write_CGNS_zones(GModel &model, const int zoneDefinition,
- const CGNSOptions &options,
- const double scalingFactor, const int vectorDim,
- const PhysGroupMap &group,
+int write_CGNS_zones(GModel &model, const int zoneDefinition, const int numZone,
+ const CGNSOptions &options, const double scalingFactor,
+ const int vectorDim, const PhysGroupMap &group,
const int cgIndexFile, const int cgIndexBase);
@@ -231,69 +257,142 @@ int GModel::writeCGNS(const std::string &name, const int zoneDefinition,
PhysGroupMap groups[4]; // vector of entities that belong to
// each physical group (in each
// dimension)
+ std::vector<DummyPartitionEntity> partitions;
+ // Dummy entities that store the
+ // elements in each partition
int numZone;
- std::vector<std::vector<MElement*> > zoneElements;
+ int meshDim;
+
+ Msg::Warning("CGNS I/O is at an \"alpha\" software stage");
+
switch(zoneDefinition) {
case 1: // By partition
+
+//--Group the elements of each partition into a dummy entity. Pointers to the
+//--entities are then made available in groups[DIM][0].
+
{
numZone = meshPartitions.size();
- zoneElements.resize(numZone);
- for(int i = 0; i != numZone; ++i)
- zoneElements[i].reserve(getMaxPartitionSize());
-// typedef GRegion Ent;
-// Ent *ent;
-// unsigned numElem[4];
-// numElem[0] = 0; numElem[1] = 0; numElem[2] = 0; numElem[3] = 0;
-// ent->getNumMeshElements(numElem);
-// int nType = ent->getNumTypeElements();
-// for(int iType = 0; iType != nType; ++iType) {
-// MElement *const start = ent->getMeshElement
-// ((iType == 0) ? 0 : numElem[iType-1]);
-// // Then can loop
-// const int nElem = numElem[iType];
-// for(int iElem = 0; iElem != nElem; ++iElem) {
-// }
-// }
-
+ if(numZone == 0) zoneDefinition = 0;
+ else {
+ partitions.resize(numZone);
+ for(int i = 0; i != numZone; ++i)
+ partitions[i].elements.reserve(getMaxPartitionSize());
+ unsigned numElem[4];
+ meshDim = getNumMeshElements(numElem);
+ // Load the dummy entities with the elements in each partition
+ switch(meshDim) {
+ case 3:
+ for(riter it = firstRegion(); it != lastRegion(); ++it) {
+ numElem[0] = 0; numElem[1] = 0; numElem[2] = 0; numElem[3] = 0;
+ (*it)->getNumMeshElements(numElem);
+ const int nType = (*it)->getNumElementTypes();
+ for(int iType = 0; iType != nType; ++iType) {
+ MElement *const *element = (*it)->getStartElementType(iType);
+ const int nElem = numElem[iType];
+ for(int iElem = 0; iElem != nElem; ++iElem) {
+ partitions[element[iElem]->getPartition() - 1].elements
+ .push_back(element[iElem]);
+ }
+ }
+ }
+ break;
+ case 2:
+ for(fiter it = firstFace(); it != lastFace(); ++it) {
+ numElem[0] = 0; numElem[1] = 0; numElem[2] = 0; numElem[3] = 0;
+ (*it)->getNumMeshElements(numElem);
+ const int nType = (*it)->getNumElementTypes();
+ for(int iType = 0; iType != nType; ++iType) {
+ MElement *const *element = (*it)->getStartElementType(iType);
+ const int nElem = numElem[iType];
+ for(int iElem = 0; iElem != nElem; ++iElem) {
+ partitions[element[iElem]->getPartition() - 1].elements
+ .push_back(element[iElem]);
+ }
+ }
+ }
+ break;
+ default:
+ Message::Error("No mesh elements were found");
+ return 0;
+ }
+ // Place pointers to the entities in the 'groups' object
+ std::vector<GEntity*> &ents = groups[meshDim][0];
+ ents.resize(numZone);
+ for(int iPart = 0; iPart != numZone; ++iPart)
+ ents[iPart] = &partitions[iPart];
+ }
}
break;
case 2: // By physical
- break;
- }
//--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.
+//--group.
- getPhysicalGroups(groups);
- if(groups[face].size() + groups[region].size() == 0) zoneDefinition = 0;
- {
+ getPhysicalGroups(groups);
+ if(groups[region].size()) {
+ numZone = groups[region].size();
+ meshDim = 3;
+ }
+ else if(groups[face].size()) {
+ numZone = groups[face].size();
+ meshDim = 2;
+ }
+ else {
+ zoneDefinition = 0; // Use single zone
+ }
+ break;
+ }
- // 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()) {
- Message::Error("No mesh elements were found");
- return 0;
+//--For a single zone, put all the entities for a dimension into groups[DIM][0]
+
+ if(zoneDefinition == 0) {
+ numZone = 1;
+ unsigned numElem[4];
+ numElem[0] = 0; numElem[1] = 0; numElem[2] = 0; numElem[3] = 0;
+ meshDim = getNumMeshElements(numElem);
+ switch(meshDim) {
+ case 3:
+ {
+ groups[region].clear();
+ std::vector<GEntity*> &ents = groups[region][0];
+ ents.resize(getNumRegions());
+ int iEnt = 0;
+ for(riter it = firstRegion(); it != lastRegion(); ++it)
+ ents[iEnt++] = *it;
+ }
+ break;
+ case 2:
+ {
+ groups[face].clear();
+ std::vector<GEntity*> &ents = groups[face][0];
+ ents.resize(getNumFaces());
+ int iEnt = 0;
+ for(fiter it = firstFace(); it != lastFace(); ++it)
+ ents[iEnt++] = *it;
}
+ break;
+ default:
+ Message::Error("No mesh elements were found");
+ return 0;
}
}
+/*--------------------------------------------------------------------*
+ * Summary of three possibilities for the 'zoneDefinition':
+ * 0) single zone: all the entities are placed in physical 0. All the
+ * entities form the zone.
+ * 1) defined by partitions: all the entities are place in physical
+ * 0. Each entity is a zone.
+ * 2) defined by physicals: all the entities in a physical form a
+ * zone.
+ *--------------------------------------------------------------------*/
+
//--Open the file and get ready to write the zones
- // Will this be a 2D or 3D mesh?
- int meshDim = 2;
+ // Get the dimension of a vector in the mesh
int vectorDim = 3;
- if(groups[region].size()) meshDim = 3;
- else vectorDim = options.vectorDim;
+ if(meshDim == 2) vectorDim = options.vectorDim;
// open the file
int cgIndexFile;
@@ -301,7 +400,8 @@ int GModel::writeCGNS(const std::string &name, const int zoneDefinition,
// write the base node
int cgIndexBase;
- if(cg_base_write(cgIndexFile, "Base_0000", meshDim, meshDim, &cgIndexBase))
+ if(cg_base_write(cgIndexFile, options.baseName.c_str(), meshDim, meshDim,
+ &cgIndexBase))
return cgnsErr();
// write information about who generated the mesh
@@ -311,13 +411,13 @@ int GModel::writeCGNS(const std::string &name, const int zoneDefinition,
switch(meshDim) {
case 2:
MZone<2>::preInit();
- write_CGNS_zones<2>(*this, zoneDefinition, options, scalingFactor,
+ write_CGNS_zones<2>(*this, zoneDefinition, numZone, options, scalingFactor,
vectorDim, groups[face], cgIndexFile, cgIndexBase);
MZone<2>::postDestroy();
break;
case 3:
MZone<3>::preInit();
- write_CGNS_zones<3>(*this, zoneDefinition, options, scalingFactor,
+ write_CGNS_zones<3>(*this, zoneDefinition, numZone, options, scalingFactor,
vectorDim, groups[region], cgIndexFile, cgIndexBase);
MZone<3>::postDestroy();
break;
@@ -438,6 +538,64 @@ void translateElementMSH2CGNS(ElementConnectivity *const zoneElemConn)
}
}
+/*******************************************************************************
+ *
+ * Routine expand_name
+ *
+ * Purpose
+ * =======
+ *
+ * Expands variables in a string 's' that are supported by the CGNS I/O. 's'
+ * is overwritten with the expanded string.
+ *
+ * - &I[0][%width]% expands into 'index'. Normally 'index' is assumed to have
+ * C numbering and therefore 1 is added to it. Option [0] prevents the
+ * addition of the one. Option [%width] sets the width of the index to
+ * 'width' and pads with leading zeros.
+ * - &N& expands to 'name'.
+ *
+ ******************************************************************************/
+
+void expand_name(std::string &s, const int index, const char *const name)
+{
+ std::string::size_type r1 = s.find('&');
+ // Look for and replace "&-&" sub-strings
+ while(r1 != std::string::npos) {
+ const std::string::size_type r2 = s.find('&', r1+1);
+ if(r2 == std::string::npos) {
+ s.replace(r1, s.length()-r1, "");
+ }
+ else {
+ const int rn = r2 - r1 + 1;
+ switch(s[r1+1]) {
+ case 'I':
+ // &I[0][%width]&
+ {
+ int iplus = 1;
+ if(s[r1+2] == '0') iplus = 0;
+ char fmt[6] = "%d";
+ const std::string::size_type f = s.find('%', r1+1);
+ if(f != std::string::npos && f < r2) {
+ int w = std::atoi(s.substr(f+1, r2-f-1).c_str());
+ if(w > 0 && w < 33) std::sprintf(fmt, "%%0%dd", w);
+ }
+ s.replace(r1, rn, CGNSNameStr(index+iplus, fmt).c_str());
+ break;
+ }
+ case 'N':
+ // &N&
+ s.replace(r1, rn, name);
+ break;
+ default:
+ s.replace(r1, rn, "");
+ }
+ }
+ if(s.length() > 1024) s = ""; // idiot/recursive check
+ r1 = s.find('&');
+ }
+
+}
+
/*******************************************************************************
*
@@ -454,71 +612,58 @@ void translateElementMSH2CGNS(ElementConnectivity *const zoneElemConn)
*
* model - (I) gmsh model
* zoneDefinition - (I) how to define the zone (see enum in code)
+ * numZone - (I) Number of zones in the domain
* options - (I) options for CGNS I/O
- * numPartitions - (I)
* group - (I) the group of physicals used to define the mesh
* globalZoneIndex - (I/O)
- * globalParittion - (I/O)
- * globalItPhysical - (I/O) a global scope iterator to the physicals
+ * globalPhysicalIt - (I/O) a global scope iterator to the physicals
* zoneIndex - (O) index of the zone
* partition - (O) partition of the zone
- * itPhysicalBegin - (O) begin physical for defining the zone
- * itPhysicalEnd - (O) end physical for defining the zone
+ * physicalItBegin - (O) begin physical for defining the zone
+ * physicalItEnd - (O) end physical for defining the zone
* zoneName - (O) name of the zone
*
******************************************************************************/
int get_zone_definition(GModel &model, const int zoneDefinition,
- const CGNSOptions &options,
- const int numPartitions, const PhysGroupMap &group,
- int &globalZoneIndex, int &globalPartition,
- PhysGroupMap::const_iterator &globalItPhysical,
+ const int numZone, const CGNSOptions &options,
+ const PhysGroupMap &group, int &globalZoneIndex,
+ PhysGroupMap::const_iterator &globalPhysicalIt,
int &zoneIndex, int &partition,
- PhysGroupMap::const_iterator &itPhysicalBegin,
- PhysGroupMap::const_iterator &itPhysicalEnd,
+ PhysGroupMap::const_iterator &physicalItBegin,
+ PhysGroupMap::const_iterator &physicalItEnd,
CGNSNameStr &zoneName)
{
- enum {
- DefineZoneSingle = 0,
- DefineZoneByPartition = 1,
- DefineZoneByPhysical = 2
- };
int status = 0;
const char *_zoneName = "Partition";
-//--Get indices for the zone
+//--Get indices for the zonex
#pragma omp critical (get_zone_definition)
{
- switch(zoneDefinition) {
- case DefineZoneSingle:
- if(globalZoneIndex > 0) status = 1; // No more zones
- else {
+ if(globalZoneIndex >= numZone) status = 1;
+ else {
+ switch(zoneDefinition) {
+ case 0: // Single zone
partition = -1;
- itPhysicalBegin = group.begin();
- itPhysicalEnd = group.end();
- }
- break;
- case DefineZoneByPartition:
- if(globalPartition > numPartitions) status = 1; // No more zones
- else {
- partition = globalPartition++;
- itPhysicalBegin = group.begin();
- itPhysicalEnd = group.end();
- }
- break;
- case DefineZoneByPhysical:
- if(globalItPhysical == group.end()) status = 1; // No more zones
- else {
+ physicalItBegin = group.begin();
+ physicalItEnd = group.end();
+ break;
+ case 1: // Zone defined by partition
+ partition = globalZoneIndex;
+ physicalItBegin = group.begin();
+ physicalItEnd = group.end();
+ break;
+ case 2: // Zone defined by physical
partition = -1;
- _zoneName = model.getPhysicalName(globalItPhysical->first).c_str();
- itPhysicalBegin = globalItPhysical++;
- itPhysicalEnd = globalItPhysical;
+ _zoneName = model.getPhysicalName(globalPhysicalIt->first).c_str();
+ physicalItBegin = globalPhysicalIt++;
+ physicalItEnd = globalPhysicalIt;
+ break;
}
- break;
+ zoneIndex = globalZoneIndex++;
}
- zoneIndex = globalZoneIndex++;
}
//omp: end critical
@@ -526,41 +671,7 @@ int get_zone_definition(GModel &model, const int zoneDefinition,
if(status == 0) {
std::string s = options.zoneName;
- std::string::size_type r1 = s.find('&');
- // Look for and replace "&&" sub-strings
- while(r1 != std::string::npos) {
- const std::string::size_type r2 = s.find('&', r1+1);
- if(r2 == std::string::npos) {
- s.replace(r1, s.length()-r1, "");
- }
- else {
- const int rn = r2 - r1 + 1;
- switch(s[r1+1]) {
- case 'I':
- // &I[0][%width]&
- {
- int iplus = 1;
- if(s[r1+2] == '0') iplus = 0;
- char fmt[6] = "%d";
- const std::string::size_type f = s.find('%', r1+1);
- if(f != std::string::npos && f < r2) {
- int w = std::atoi(s.substr(f+1, r2-f-1).c_str());
- if(w > 0 && w < 33) std::sprintf(fmt, "%%0%dd", w);
- }
- s.replace(r1, rn, CGNSNameStr(zoneIndex+iplus, fmt).c_str());
- break;
- }
- case 'N':
- // &N&
- s.replace(r1, rn, _zoneName);
- break;
- default:
- s.replace(r1, rn, "");
- }
- }
- if(s.length() > 1024) s = ""; // idiot/recursive check
- r1 = s.find('&');
- }
+ expand_name(s, zoneIndex, _zoneName);
if(s.length() == 0) { // If empty
s = "Zone_";
s += CGNSNameStr(zoneIndex+1).c_str();
@@ -569,6 +680,7 @@ int get_zone_definition(GModel &model, const int zoneDefinition,
}
return status;
+
}
@@ -633,26 +745,23 @@ struct ZoneInfo
};
template <unsigned DIM>
-int write_CGNS_zones(GModel &model, const int zoneDefinition,
- const CGNSOptions &options,
- const double scalingFactor, const int vectorDim,
- const PhysGroupMap &group,
+int write_CGNS_zones(GModel &model, const int zoneDefinition, const int numZone,
+ const CGNSOptions &options, const double scalingFactor,
+ const int vectorDim, const PhysGroupMap &group,
const int cgIndexFile, const int cgIndexBase)
{
//--Shared data
- const int numPartitions = model.getMeshPartitions().size();
int threadsWorking = omp_get_num_threads();
// Semaphore for working threads
omp_lock_t threadWLock;
std::queue<ZoneTask<DIM>*> zoneQueue; // Queue for zones that have been
// defined and are ready to be written
omp_lock_t queueLock;
- // Next 3 are locked by omp critical
+ // Next two are locked by an omp critical
int globalZoneIndex = 0;
- int globalPartition = 1;
- PhysGroupMap::const_iterator globalItPhysical = group.begin();
+ PhysGroupMap::const_iterator globalPhysicalIt = group.begin();
//--Initialize omp locks
@@ -676,7 +785,7 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition,
std::vector<double> dBuffer; // Buffer for double-precision data
std::vector<int> iBuffer1, iBuffer2;
// Buffers for integer data
- int indexInterface = 0; // Index for interfaces
+ int interfaceIndex = 0; // Index for interfaces
while (threadsWorking || zoneQueue.size()) {
if (zoneQueue.size()) {
@@ -700,8 +809,6 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition,
writeTask->zoneName.c_str(), cgZoneSize,
Unstructured, &cgIndexZone))
return cgnsErr();
-// std::cout << "Zone size: " << cgZoneSize[0] << ' ' // DBG
-// << cgZoneSize[1] << ' ' << cgZoneSize[2] << std::endl;
// 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())
@@ -794,22 +901,26 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition,
iBuffer2[iVert] = vp.vertexIndex2;
}
int cgIndexInterface;
+ std::string interfaceName = options.interfaceName;
+ expand_name(interfaceName, interfaceIndex++, "Interface");
+ if(interfaceName.length() == 0) {
+ interfaceName = "Interface_";
+ interfaceName += CGNSNameStr(interfaceIndex).c_str();
+ }
// In the first zone
if(cg_conn_write
(cgIndexFile, cgIndexBase, zoneInfo[gCIt->first.zone1].cgIndex,
- CGNSNameStr(++indexInterface, "Interface_%d").c_str(), Vertex,
- Abutting1to1, PointList, nVert, &iBuffer1[0],
- zoneInfo[gCIt->first.zone2].name.c_str(), Unstructured,
- PointListDonor, Integer, nVert, &iBuffer2[0],
+ interfaceName.c_str(), Vertex, Abutting1to1, PointList, nVert,
+ &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,
- CGNSNameStr(++indexInterface, "Interface_%d").c_str(), Vertex,
- Abutting1to1, PointList, nVert, &iBuffer2[0],
- zoneInfo[gCIt->first.zone1].name.c_str(), Unstructured,
- PointListDonor, Integer, nVert, &iBuffer1[0],
+ interfaceName.c_str(), Vertex, Abutting1to1, PointList, nVert,
+ &iBuffer2[0], zoneInfo[gCIt->first.zone1].name.c_str(),
+ Unstructured, PointListDonor, Integer, nVert, &iBuffer1[0],
&cgIndexInterface))
return cgnsErr();
}
@@ -830,15 +941,14 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition,
*--------------------------------------------------------------------*/
else {
- PhysGroupMap::const_iterator itPhysicalBegin;
- PhysGroupMap::const_iterator itPhysicalEnd;
+ PhysGroupMap::const_iterator physicalItBegin;
+ PhysGroupMap::const_iterator physicalItEnd;
int zoneIndex;
int partition;
if(get_zone_definition
- (model, zoneDefinition, options, numPartitions, group,
- globalZoneIndex, globalPartition, globalItPhysical,
- zoneTask.zoneIndex, partition, itPhysicalBegin, itPhysicalEnd,
- zoneTask.zoneName)) {
+ (model, zoneDefinition, numZone, options, group, globalZoneIndex,
+ globalPhysicalIt, zoneTask.zoneIndex, partition, physicalItBegin,
+ physicalItEnd, zoneTask.zoneName)) {
omp_set_lock(&threadWLock);
--threadsWorking;
omp_unset_lock(&threadWLock);
@@ -846,14 +956,21 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition,
else {
zoneTask.zone.clear();
// Loop through all physical in the zone definition
- for(PhysGroupMap::const_iterator itPhysical = itPhysicalBegin;
- itPhysical != itPhysicalEnd; ++itPhysical) {
+ for(PhysGroupMap::const_iterator physicalIt = physicalItBegin;
+ physicalIt != physicalItEnd; ++physicalIt) {
// Add elements from all entities in the physical with defined
// partition number
- zoneTask.zone.template add_elements_in_entities
- <typename std::vector<GEntity*>::const_iterator>
- (itPhysical->second.begin(), itPhysical->second.end(),
- partition);
+ if(partition == -1) {
+ zoneTask.zone.template add_elements_in_entities
+ <typename std::vector<GEntity*>::const_iterator>
+ (physicalIt->second.begin(), physicalIt->second.end());
+ }
+ else {
+ zoneTask.zone.template add_elements_in_entities
+ <typename std::vector<GEntity*>::const_iterator>
+ (physicalIt->second.begin() + partition,
+ physicalIt->second.begin() + (partition+1));
+ }
}
// Process the zone
zoneTask.zone.zoneData();
@@ -904,9 +1021,14 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition,
const int numBCVert = iVert - iVertStart;
int cgIndexBoco;
+ std::string patchName = options.patchName;
+ expand_name(patchName, patchIndex-1, "Patch");
+ if(patchName.length() == 0) {
+ patchName = "Patch_";
+ patchName += CGNSNameStr(patchIndex).c_str();
+ }
if(cg_boco_write(cgIndexFile, cgIndexBase,
- zoneInfo[zoneIndex].cgIndex,
- CGNSNameStr(patchIndex, "Patch %d").c_str(),
+ zoneInfo[zoneIndex].cgIndex, patchName.c_str(),
BCTypeNull, PointList, numBCVert, &iBuffer1[0],
&cgIndexBoco))
return cgnsErr();
diff --git a/Geo/MZoneBoundary.cpp b/Geo/MZoneBoundary.cpp
index 648a814b1c8814d22279290752b1dc07643f5d1e..07223536945934f4152bcf4c71ffaf7481dcaf87 100644
--- a/Geo/MZoneBoundary.cpp
+++ b/Geo/MZoneBoundary.cpp
@@ -187,6 +187,7 @@ void updateBoVec
* ===
*
* vertex - (I) The vertex to find the normal at
+ * zoneIndex - (I) Zone being worked on
* edge - (I) The geometry edge upon which the vertex resides
* and from which the normal will be determined
* faces - (I) All faces on the boundary connected to 'vertex'
@@ -197,10 +198,11 @@ void updateBoVec
*============================================================================*/
int edge_normal
-(const MVertex *const vertex, const GEdge *const gEdge,
+(const MVertex *const vertex, const int zoneIndex, const GEdge *const gEdge,
const CCon::FaceVector<MZoneBoundary<2>::GlobalVertexData<MEdge>::FaceDataB>
&faces, SVector3 &boNormal)
{
+
const double par = gEdge->parFromPoint(vertex->point());
if(par == std::numeric_limits<double>::max()) return 1;
@@ -210,17 +212,22 @@ int edge_normal
SVector3 meshPlaneNormal(0.); // This normal is perpendicular to the
// plane of the mesh
- // The interior point and mesh plane normal are computed from all elements.
+ // The interior point and mesh plane normal are computed from all elements in
+ // the zone.
+ int cFace = 0;
const int nFace = faces.size();
for(int iFace = 0; iFace != nFace; ++iFace) {
- interior += faces[iFace].parentElement->barycenter();
- // Make sure all the planes go in the same direction
- //**Required?
- SVector3 mpnt = faces[iFace].parentElement->getFace(0).normal();
- if(dot(mpnt, meshPlaneNormal) < 0.) mpnt.negate();
- meshPlaneNormal += mpnt;
+ if(faces[iFace].zoneIndex == zoneIndex) {
+ ++cFace;
+ interior += faces[iFace].parentElement->barycenter();
+ // Make sure all the planes go in the same direction
+ //**Required?
+ SVector3 mpnt = faces[iFace].parentElement->getFace(0).normal();
+ if(dot(mpnt, meshPlaneNormal) < 0.) mpnt.negate();
+ meshPlaneNormal += mpnt;
+ }
}
- interior /= nFace;
+ interior /= cFace;
// Normal to the boundary edge (but unknown direction)
boNormal = crossprod(tangent, meshPlaneNormal);
boNormal.normalize();
@@ -229,6 +236,7 @@ int edge_normal
if(dot(boNormal, SVector3(vertex->point(), interior)) < 0.)
boNormal.negate();
return 0;
+
}
@@ -244,8 +252,6 @@ void updateBoVec<2, MEdge>
const CCon::FaceVector<MZoneBoundary<2>::GlobalVertexData<MEdge>::FaceDataB>
&faces, ZoneBoVec &zoneBoVec, BCPatchIndex &patch, bool &warnNormFromElem)
{
- SVector3 boNormal; // Normal to the boundary face (edge in
- // 2D)
GEntity *const ent = vertex->onWhat();
if(ent == 0) {
@@ -274,7 +280,7 @@ void updateBoVec<2, MEdge>
for(int iFace = 0; iFace != nFace; ++iFace) {
if(zoneIndex == faces[iFace].zoneIndex) {
MEdge mEdge = faces[iFace].parentElement->getEdge
- (faces[iFace].parentFace);
+ (faces[iFace].parentFace);
// Get the other vertex on the mesh edge.
MVertex *vertex2 = mEdge.getMinVertex();
if(vertex2 == vertex) vertex2 = mEdge.getMaxVertex();
@@ -282,7 +288,9 @@ void updateBoVec<2, MEdge>
// is also connected to vertex. If so, add it to the container of
// edge entities that will be used to determine the normal
GEntity *const ent2 = vertex2->onWhat();
- if(ent2->dim() == 1) useGEdge.push_back(static_cast<GEdge*>(ent2));
+ if(ent2->dim() == 1) {
+ useGEdge.push_back(static_cast<GEdge*>(ent2));
+ }
}
}
@@ -292,14 +300,17 @@ void updateBoVec<2, MEdge>
switch(useGEdge.size()) {
case 0:
- goto getNormalFromElements;
+// goto getNormalFromElements;
+ // We probably don't want BC data if none of the faces attatched to
+ // this vertex and in this zone are on the boundary.
+ break;
case 1:
{
const GEdge *const gEdge =
static_cast<const GEdge*>(useGEdge[0]);
SVector3 boNormal;
- if(edge_normal(vertex, gEdge, faces, boNormal))
+ if(edge_normal(vertex, zoneIndex, gEdge, faces, boNormal))
goto getNormalFromElements;
zoneBoVec.push_back(VertexBoundary(zoneIndex, gEdge->tag(),
boNormal,
@@ -315,14 +326,14 @@ void updateBoVec<2, MEdge>
const GEdge *const gEdge1 =
static_cast<const GEdge*>(useGEdge[0]);
SVector3 boNormal1;
- if(edge_normal(vertex, gEdge1, faces, boNormal1))
+ if(edge_normal(vertex, zoneIndex, gEdge1, faces, boNormal1))
goto getNormalFromElements;
// Get the second normal
const GEdge *const gEdge2 =
static_cast<const GEdge*>(useGEdge[1]);
SVector3 boNormal2;
- if(edge_normal(vertex, gEdge2, faces, boNormal2))
+ if(edge_normal(vertex, zoneIndex, gEdge2, faces, boNormal2))
goto getNormalFromElements;
if(dot(boNormal1, boNormal2) < 0.98) {
@@ -371,7 +382,8 @@ void updateBoVec<2, MEdge>
{
SVector3 boNormal;
- if(edge_normal(vertex, static_cast<const GEdge*>(ent), faces, boNormal))
+ if(edge_normal(vertex, zoneIndex, static_cast<const GEdge*>(ent), faces,
+ boNormal))
goto getNormalFromElements;
zoneBoVec.push_back(VertexBoundary(zoneIndex, ent->tag(), boNormal,
const_cast<MVertex*>(vertex),
@@ -405,11 +417,13 @@ void updateBoVec<2, MEdge>
// plane of the mesh
const int nFace = faces.size();
for(int iFace = 0; iFace != nFace; ++iFace) {
+ if(faces[iFace].zoneIndex == zoneIndex) {
// Make sure all the planes go in the same direction
//**Required?
- SVector3 mpnt = faces[iFace].parentElement->getFace(0).normal();
- if(dot(mpnt, meshPlaneNormal) < 0.) mpnt.negate();
- meshPlaneNormal += mpnt;
+ SVector3 mpnt = faces[iFace].parentElement->getFace(0).normal();
+ if(dot(mpnt, meshPlaneNormal) < 0.) mpnt.negate();
+ meshPlaneNormal += mpnt;
+ }
}
// Sum the normals from each element. The tangent is computed from all
// faces in the zone attached to the vertex and is weighted by the length of
@@ -417,7 +431,7 @@ void updateBoVec<2, MEdge>
// inwards-pointing normal.
SVector3 boNormal(0.);
for(int iFace = 0; iFace != nFace; ++iFace) {
- if(zoneIndex == faces[iFace].zoneIndex) {
+ if(faces[iFace].zoneIndex == zoneIndex) {
const SVector3 tangent = faces[iFace].parentElement->getEdge
(faces[iFace].parentFace).tangent();
// Normal to the boundary (unknown direction)
@@ -599,15 +613,20 @@ void updateBoVec<3, MFace>
for(std::list<const GFace*>::const_iterator gFIt = useGFace.begin();
gFIt != useGFace.end(); ++gFIt) {
const SPoint2 par = (*gFIt)->parFromPoint(vertex->point());
- if(par.x() == std::numeric_limits<double>::max())
+ if(par.x() == std::numeric_limits<double>::max()) //**?
goto getNormalFromElements; // :P After all that!
SVector3 boNormal = (*gFIt)->normal(par);
SPoint3 interior(0., 0., 0.);
+ int cFace = 0;
const int nFace = faces.size();
- for(int iFace = 0; iFace != nFace; ++iFace)
- interior += faces[iFace].parentElement->barycenter();
- interior /= nFace;
+ for(int iFace = 0; iFace != nFace; ++iFace) {
+ if(faces[iFace].zoneIndex == zoneIndex) {
+ ++cFace;
+ interior += faces[iFace].parentElement->barycenter();
+ }
+ }
+ interior /= cFace;
if(dot(boNormal, SVector3(vertex->point(), interior)) < 0.)
boNormal.negate();
@@ -633,10 +652,15 @@ void updateBoVec<3, MFace>
SVector3 boNormal = static_cast<const GFace*>(ent)->normal(par);
SPoint3 interior(0., 0., 0.);
+ int cFace = 0;
const int nFace = faces.size();
- for(int iFace = 0; iFace != nFace; ++iFace)
- interior += faces[iFace].parentElement->barycenter();
- interior /= nFace;
+ for(int iFace = 0; iFace != nFace; ++iFace) {
+ if(faces[iFace].zoneIndex == zoneIndex) {
+ ++cFace;
+ interior += faces[iFace].parentElement->barycenter();
+ }
+ }
+ interior /= cFace;
if(dot(boNormal, SVector3(vertex->point(), interior)) < 0.)
boNormal.negate();
@@ -672,14 +696,16 @@ void updateBoVec<3, MFace>
SVector3 boNormal(0.);
const int nFace = faces.size();
for(int iFace = 0; iFace != nFace; ++iFace) {
- // Normal to the boundary (unknown direction)
- SVector3 bnt = faces[iFace].parentElement->getFace
- (faces[iFace].parentFace).normal();
- // Inwards normal
- const SVector3 inwards(vertex->point(),
- faces[iFace].parentElement->barycenter());
- if(dot(bnt, inwards) < 0.) bnt.negate();
- boNormal += bnt;
+ if(faces[iFace].zoneIndex == zoneIndex) {
+ // Normal to the boundary (unknown direction)
+ SVector3 bnt = faces[iFace].parentElement->getFace
+ (faces[iFace].parentFace).normal();
+ // Inwards normal
+ const SVector3 inwards(vertex->point(),
+ faces[iFace].parentElement->barycenter());
+ if(dot(bnt, inwards) < 0.) bnt.negate();
+ boNormal += bnt;
+ }
}
boNormal.normalize();
zoneBoVec.push_back(VertexBoundary(zoneIndex, 0, boNormal,
@@ -735,6 +761,12 @@ int MZoneBoundary<DIM>::interiorBoundaryVertices
//--Find or insert this vertex into the global map
+// bool debug = false;
+// if(vMapIt->first->x() == 1. && vMapIt->first->y() == 0.) {
+// std::cout << "Working with vertex(1, 0): " << vMapIt->first
+// << " for zone " << newZoneIndex << std::endl;
+// debug = true;
+// }
std::pair<typename GlobalBoVertexMap::iterator, bool> insGblBoVertMap =
globalBoVertMap.insert(std::pair<const MVertex*, GlobalVertexData<FaceT> >
(vMapIt->first, GlobalVertexData<FaceT>()));
@@ -748,9 +780,21 @@ 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();
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";
+// }
globalVertData.faces.push_back
(typename GlobalVertexData<FaceT>::FaceDataB
(newZoneIndex, zoneVertData.faces[iFace]));
@@ -781,26 +825,33 @@ int MZoneBoundary<DIM>::interiorBoundaryVertices
}
// Update the list of faces attached to this vertex
- const int nGFace = globalVertData.faces.size();
- // This is the maximum number of faces
- // searched from 'globalVertData'. This
- // size may increase as faces are added
- // and the order may also change as
- // faces are deleted. However, only the
- // faces before nGFace are important.
+ unsigned nGFace = globalVertData.faces.size();
+ // This is the number of faces searched
+ // from 'globalVertData'. It will only
+ // decrease if the size() is less.
+ // Since a FaceVector swaps the last
+ // element into the erased index, this
+ // implies that if new faces are added,
+ // then old ones deleted, some of the
+ // faces from zoneVertData may also be
+ // searched. This is okay since they
+ // are all unique.
const typename MZone<DIM>::BoFaceMap::const_iterator *zFace =
&zoneVertData.faces[0];
for(int nZFace = zoneVertData.faces.size(); nZFace--;) {
- int iGFace = 0;
- while(iGFace != nGFace) {
+ bool foundMatch = false;
+ for(int iGFace = 0; iGFace != nGFace; ++iGFace) {
if((*zFace)->first == globalVertData.faces[iGFace].face) {
- // Faces match - delete from 'globalVertData'
+ foundMatch = true;
+ // Faces match - delete from 'globalVertData'.
globalVertData.faces.erase(iGFace);
+ // Erasing from the FaceVector swaps the last element into this
+ // index. We only decrease nGFace if the size is less.
+ nGFace = std::min(globalVertData.faces.size(), nGFace);
break;
}
- ++iGFace;
}
- if(iGFace == nGFace) {
+ if(!foundMatch) {
// New face - add to 'globalVertData'
globalVertData.faces.push_back
(typename GlobalVertexData<FaceT>::FaceDataB(newZoneIndex, *zFace));
@@ -812,6 +863,12 @@ int MZoneBoundary<DIM>::interiorBoundaryVertices
// and it may be deleted.
if(globalVertData.faces.size() == 0)
globalBoVertMap.erase(insGblBoVertMap.first);
+ else {
+ // Update the list of zones attached to this vertex
+ globalVertData.zoneData.push_back
+ (typename GlobalVertexData<FaceT>::ZoneData
+ (zoneVertData.index, newZoneIndex));
+ }
}
} // End loop over boundary vertices