Skip to content
Snippets Groups Projects
Commit 36df5498 authored by Stefen Guzik's avatar Stefen Guzik
Browse files

CGNS bugs

parent 96aa0df8
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Geo/CGNSOptions.h
+ 5
− 1
View file @ 36df5498
......@@ -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&")
{ }
};
......
Geo/GModelIO_CGNS.cpp
+ 282
− 160
View file @ 36df5498
......@@ -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();
......
Geo/MZoneBoundary.cpp
+ 106
− 49
View file @ 36df5498
......@@ -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
......
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