* @param orderedNodes list of nodes to reorder (in out)
* @return size of the list
*/
-int SMDS_UnstructuredGrid::getOrderedNodesOfFace(int vtkVolId, std::vector<vtkIdType>& orderedNodes)
+int SMDS_UnstructuredGrid::getOrderedNodesOfFace(int vtkVolId, int& dim, std::vector<vtkIdType>& orderedNodes)
{
int vtkType = this->GetCellType(vtkVolId);
- int cellDim = SMDS_Downward::getCellDimension(vtkType);
- if (cellDim != 3)
- return 0;
- SMDS_Down3D *downvol = static_cast<SMDS_Down3D*> (_downArray[vtkType]);
- int downVolId = this->_cellIdToDownId[vtkVolId];
- downvol->getOrderedNodesOfFace(downVolId, orderedNodes);
+ dim = SMDS_Downward::getCellDimension(vtkType);
+ if (dim == 3)
+ {
+ SMDS_Down3D *downvol = static_cast<SMDS_Down3D*> (_downArray[vtkType]);
+ int downVolId = this->_cellIdToDownId[vtkVolId];
+ downvol->getOrderedNodesOfFace(downVolId, orderedNodes);
+ }
+ // else nothing to do;
return orderedNodes.size();
}
* @param nodeDomains: map(original id --> map(domain --> duplicated node id))
* @return ok if success.
*/
-SMDS_MeshVolume* SMDS_UnstructuredGrid::extrudeVolumeFromFace(int vtkVolId,
+SMDS_MeshCell* SMDS_UnstructuredGrid::extrudeVolumeFromFace(int vtkVolId,
int domain1,
int domain2,
std::set<int>& originalNodes,
for (; it != originalNodes.end(); ++it)
orderedOriginals.push_back(*it);
- int nbNodes = this->getOrderedNodesOfFace(vtkVolId, orderedOriginals);
+ int dim = 0;
+ int nbNodes = this->getOrderedNodesOfFace(vtkVolId, dim, orderedOriginals);
vector<vtkIdType> orderedNodes;
+ bool isQuadratic = false;
switch (orderedOriginals.size())
{
case 3:
- case 4:
- for (int i = 0; i < nbNodes; i++)
- orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain1]);
- for (int i = 0; i < nbNodes; i++)
- orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain2]);
+ if (dim == 2)
+ isQuadratic = true;
break;
case 6:
case 8:
- {
- long dom1 = domain1;
- long dom2 = domain2;
- long dom1_2; // for nodeQuadDomains
- if (domain1 < domain2)
- dom1_2 = dom1 + INT_MAX * dom2;
- else
- dom1_2 = dom2 + INT_MAX * dom1;
- //cerr << "dom1=" << dom1 << " dom2=" << dom2 << " dom1_2=" << dom1_2 << endl;
- int ima = orderedOriginals.size();
- int mid = orderedOriginals.size() / 2;
- //cerr << "ima=" << ima << " mid=" << mid << endl;
- for (int i = 0; i < mid; i++)
- orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain1]);
- for (int i = 0; i < mid; i++)
- orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain2]);
- for (int i = mid; i < ima; i++)
- orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain1]);
- for (int i = mid; i < ima; i++)
- orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain2]);
- for (int i = 0; i < mid; i++)
- {
- int oldId = orderedOriginals[i];
- int newId;
- if (nodeQuadDomains.count(oldId) && nodeQuadDomains[oldId].count(dom1_2))
- newId = nodeQuadDomains[oldId][dom1_2];
- else
- {
- double *coords = this->GetPoint(oldId);
- SMDS_MeshNode *newNode = _mesh->AddNode(coords[0], coords[1], coords[2]);
- newId = newNode->getVtkId();
- std::map<long, int> emptyMap;
- nodeQuadDomains[oldId] = emptyMap;
- nodeQuadDomains[oldId][dom1_2] = newId;
- }
- orderedNodes.push_back(newId);
- }
- }
+ isQuadratic = true;
break;
default:
- ASSERT(0);
+ isQuadratic = false;
+ break;
}
- SMDS_MeshVolume *vol = _mesh->AddVolumeFromVtkIds(orderedNodes);
+ if (isQuadratic)
+ {
+ long dom1 = domain1;
+ long dom2 = domain2;
+ long dom1_2; // for nodeQuadDomains
+ if (domain1 < domain2)
+ dom1_2 = dom1 + INT_MAX * dom2;
+ else
+ dom1_2 = dom2 + INT_MAX * dom1;
+ //cerr << "dom1=" << dom1 << " dom2=" << dom2 << " dom1_2=" << dom1_2 << endl;
+ int ima = orderedOriginals.size();
+ int mid = orderedOriginals.size() / 2;
+ //cerr << "ima=" << ima << " mid=" << mid << endl;
+ for (int i = 0; i < mid; i++)
+ orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain1]);
+ for (int i = 0; i < mid; i++)
+ orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain2]);
+ for (int i = mid; i < ima; i++)
+ orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain1]);
+ for (int i = mid; i < ima; i++)
+ orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain2]);
+ for (int i = 0; i < mid; i++)
+ {
+ int oldId = orderedOriginals[i];
+ int newId;
+ if (nodeQuadDomains.count(oldId) && nodeQuadDomains[oldId].count(dom1_2))
+ newId = nodeQuadDomains[oldId][dom1_2];
+ else
+ {
+ double *coords = this->GetPoint(oldId);
+ SMDS_MeshNode *newNode = _mesh->AddNode(coords[0], coords[1], coords[2]);
+ newId = newNode->getVtkId();
+ std::map<long, int> emptyMap;
+ nodeQuadDomains[oldId] = emptyMap;
+ nodeQuadDomains[oldId][dom1_2] = newId;
+ }
+ orderedNodes.push_back(newId);
+ }
+ }
+ else
+ {
+ for (int i = 0; i < nbNodes; i++)
+ orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain1]);
+ if (dim == 3)
+ for (int i = 0; i < nbNodes; i++)
+ orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain2]);
+ else
+ for (int i = nbNodes-1; i >= 0; i--)
+ orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain2]);
+
+ }
+
+ if (dim == 3)
+ {
+ SMDS_MeshVolume *vol = _mesh->AddVolumeFromVtkIds(orderedNodes);
+ return vol;
+ }
+ else if (dim == 2)
+ {
+ SMDS_MeshFace *face = _mesh->AddFaceFromVtkIds(orderedNodes);
+ return face;
+ }
// TODO update sub-shape list of elements and nodes
- return vol;
+ return 0;
}
// and corresponding volume of this domain, for each shared face.
// a volume has a face shared by 2 domains if it has a neighbor which is not in is domain.
+ //MESSAGE("Domain " << idom);
const TIDSortedElemSet& domain = theElems[idom];
TIDSortedElemSet::const_iterator elemItr = domain.begin();
for (; elemItr != domain.end(); ++elemItr)
if (!anElem)
continue;
int vtkId = anElem->getVtkId();
+ //MESSAGE(" vtkId " << vtkId << " smdsId " << anElem->GetID());
int neighborsVtkIds[NBMAXNEIGHBORS];
int downIds[NBMAXNEIGHBORS];
unsigned char downTypes[NBMAXNEIGHBORS];
{
faceDomains[face][idom] = vtkId; // volume associated to face in this domain
celldom[vtkId] = idom;
+ //MESSAGE(" cell with a border " << vtkId << " domain " << idom);
}
}
}
for (int idomain = 0; idomain < theElems.size(); idomain++)
{
+ //MESSAGE("Domain " << idomain);
const TIDSortedElemSet& domain = theElems[idomain];
itface = faceDomains.begin();
for (; itface != faceDomains.end(); ++itface)
continue;
cellDomains[aCell][idomain] = vtkId;
celldom[vtkId] = idomain;
+ //MESSAGE(" cell " << vtkId << " domain " << idomain);
}
}
}
// the value is the ordered domain ids. (more than 4 domains not taken into account)
std::map<std::vector<int>, std::vector<int> > edgesMultiDomains; // nodes of edge --> ordered domains
- std::map<int, std::vector<int> > mutipleNodes; // nodes muti domains with domain order
+ std::map<int, std::vector<int> > mutipleNodes; // nodes multi domains with domain order
+ std::map<int, std::vector<int> > mutipleNodesToFace; // nodes multi domains with domain order to transform in Face (junction between 3 or more 2D domains)
for (int idomain = 0; idomain < theElems.size(); idomain++)
{
//MESSAGE("multiple Nodes detected on a shared face");
int downId = itface->first.cellId;
unsigned char cellType = itface->first.cellType;
- int nbEdges = grid->getDownArray(cellType)->getNumberOfDownCells(downId);
- const int *downEdgeIds = grid->getDownArray(cellType)->getDownCells(downId);
- const unsigned char* edgeType = grid->getDownArray(cellType)->getDownTypes(downId);
- for (int ie =0; ie < nbEdges; ie++)
+ // --- shared edge or shared face ?
+ if ((cellType == VTK_LINE) || (cellType == VTK_QUADRATIC_EDGE)) // shared edge (between two faces)
{
int nodes[3];
- int nbNodes = grid->getDownArray(edgeType[ie])->getNodes(downEdgeIds[ie], nodes);
- if (mutipleNodes.count(nodes[0]) && mutipleNodes.count(nodes[nbNodes-1]))
+ int nbNodes = grid->getDownArray(cellType)->getNodes(downId, nodes);
+ for (int i=0; i< nbNodes; i=i+nbNodes-1) // i=0 , i=nbNodes-1
+ if (mutipleNodes.count(nodes[i]))
+ if (!mutipleNodesToFace.count(nodes[i]))
+ mutipleNodesToFace[nodes[i]] = mutipleNodes[nodes[i]];
+ }
+ else // shared face (between two volumes)
+ {
+ int nbEdges = grid->getDownArray(cellType)->getNumberOfDownCells(downId);
+ const int* downEdgeIds = grid->getDownArray(cellType)->getDownCells(downId);
+ const unsigned char* edgeType = grid->getDownArray(cellType)->getDownTypes(downId);
+ for (int ie =0; ie < nbEdges; ie++)
{
- vector<int> vn0 = mutipleNodes[nodes[0]];
- vector<int> vn1 = mutipleNodes[nodes[nbNodes - 1]];
- sort( vn0.begin(), vn0.end() );
- sort( vn1.begin(), vn1.end() );
- if (vn0 == vn1)
+ int nodes[3];
+ int nbNodes = grid->getDownArray(edgeType[ie])->getNodes(downEdgeIds[ie], nodes);
+ if (mutipleNodes.count(nodes[0]) && mutipleNodes.count(nodes[nbNodes-1]))
{
- //MESSAGE(" detect edgesMultiDomains " << nodes[0] << " " << nodes[nbNodes - 1]);
- double *coords = grid->GetPoint(nodes[0]);
- gp_Pnt p0(coords[0], coords[1], coords[2]);
- coords = grid->GetPoint(nodes[nbNodes - 1]);
- gp_Pnt p1(coords[0], coords[1], coords[2]);
- gp_Pnt gref;
- int vtkVolIds[1000]; // an edge can belong to a lot of volumes
- map<int, SMDS_VtkVolume*> domvol; // domain --> a volume with the edge
- map<int, double> angleDom; // oriented angles between planes defined by edge and volume centers
- int nbvol = grid->GetParentVolumes(vtkVolIds, downEdgeIds[ie], edgeType[ie]);
- for (int id=0; id < vn0.size(); id++)
+ vector<int> vn0 = mutipleNodes[nodes[0]];
+ vector<int> vn1 = mutipleNodes[nodes[nbNodes - 1]];
+ sort( vn0.begin(), vn0.end() );
+ sort( vn1.begin(), vn1.end() );
+ if (vn0 == vn1)
{
- int idom = vn0[id];
- for (int ivol=0; ivol<nbvol; ivol++)
+ //MESSAGE(" detect edgesMultiDomains " << nodes[0] << " " << nodes[nbNodes - 1]);
+ double *coords = grid->GetPoint(nodes[0]);
+ gp_Pnt p0(coords[0], coords[1], coords[2]);
+ coords = grid->GetPoint(nodes[nbNodes - 1]);
+ gp_Pnt p1(coords[0], coords[1], coords[2]);
+ gp_Pnt gref;
+ int vtkVolIds[1000]; // an edge can belong to a lot of volumes
+ map<int, SMDS_VtkVolume*> domvol; // domain --> a volume with the edge
+ map<int, double> angleDom; // oriented angles between planes defined by edge and volume centers
+ int nbvol = grid->GetParentVolumes(vtkVolIds, downEdgeIds[ie], edgeType[ie]);
+ for (int id=0; id < vn0.size(); id++)
{
- int smdsId = meshDS->fromVtkToSmds(vtkVolIds[ivol]);
- SMDS_MeshElement* elem = (SMDS_MeshElement*)meshDS->FindElement(smdsId);
- if (theElems[idom].count(elem))
+ int idom = vn0[id];
+ for (int ivol=0; ivol<nbvol; ivol++)
{
- SMDS_VtkVolume* svol = dynamic_cast<SMDS_VtkVolume*>(elem);
- domvol[idom] = svol;
- //MESSAGE(" domain " << idom << " volume " << elem->GetID());
- double values[3];
- vtkIdType npts = 0;
- vtkIdType* pts = 0;
- grid->GetCellPoints(vtkVolIds[ivol], npts, pts);
- SMDS_VtkVolume::gravityCenter(grid, pts, npts, values);
- if (id ==0)
- {
- gref.SetXYZ(gp_XYZ(values[0], values[1], values[2]));
- angleDom[idom] = 0;
- }
- else
+ int smdsId = meshDS->fromVtkToSmds(vtkVolIds[ivol]);
+ SMDS_MeshElement* elem = (SMDS_MeshElement*)meshDS->FindElement(smdsId);
+ if (theElems[idom].count(elem))
{
- gp_Pnt g(values[0], values[1], values[2]);
- angleDom[idom] = OrientedAngle(p0, p1, gref, g); // -pi<angle<+pi
- //MESSAGE(" angle=" << angleDom[idom]);
+ SMDS_VtkVolume* svol = dynamic_cast<SMDS_VtkVolume*>(elem);
+ domvol[idom] = svol;
+ //MESSAGE(" domain " << idom << " volume " << elem->GetID());
+ double values[3];
+ vtkIdType npts = 0;
+ vtkIdType* pts = 0;
+ grid->GetCellPoints(vtkVolIds[ivol], npts, pts);
+ SMDS_VtkVolume::gravityCenter(grid, pts, npts, values);
+ if (id ==0)
+ {
+ gref.SetXYZ(gp_XYZ(values[0], values[1], values[2]));
+ angleDom[idom] = 0;
+ }
+ else
+ {
+ gp_Pnt g(values[0], values[1], values[2]);
+ angleDom[idom] = OrientedAngle(p0, p1, gref, g); // -pi<angle<+pi
+ //MESSAGE(" angle=" << angleDom[idom]);
+ }
+ break;
}
- break;
}
}
+ map<double, int> sortedDom; // sort domains by angle
+ for (map<int, double>::iterator ia = angleDom.begin(); ia != angleDom.end(); ++ia)
+ sortedDom[ia->second] = ia->first;
+ vector<int> vnodes;
+ vector<int> vdom;
+ for (map<double, int>::iterator ib = sortedDom.begin(); ib != sortedDom.end(); ++ib)
+ {
+ vdom.push_back(ib->second);
+ //MESSAGE(" ordered domain " << ib->second << " angle " << ib->first);
+ }
+ for (int ino = 0; ino < nbNodes; ino++)
+ vnodes.push_back(nodes[ino]);
+ edgesMultiDomains[vnodes] = vdom; // nodes vector --> ordered domains
}
- map<double, int> sortedDom; // sort domains by angle
- for (map<int, double>::iterator ia = angleDom.begin(); ia != angleDom.end(); ++ia)
- sortedDom[ia->second] = ia->first;
- vector<int> vnodes;
- vector<int> vdom;
- for (map<double, int>::iterator ib = sortedDom.begin(); ib != sortedDom.end(); ++ib)
- {
- vdom.push_back(ib->second);
- //MESSAGE(" ordered domain " << ib->second << " angle " << ib->first);
- }
- for (int ino = 0; ino < nbNodes; ino++)
- vnodes.push_back(nodes[ino]);
- edgesMultiDomains[vnodes] = vdom; // nodes vector --> ordered domains
}
}
}
int vtkVolId = itdom->second;
itdom++;
int dom2 = itdom->first;
- SMDS_MeshVolume *vol = grid->extrudeVolumeFromFace(vtkVolId, dom1, dom2, oldNodes, nodeDomains,
+ SMDS_MeshCell *vol = grid->extrudeVolumeFromFace(vtkVolId, dom1, dom2, oldNodes, nodeDomains,
nodeQuadDomains);
stringstream grpname;
grpname << "j_";
int idg;
string namegrp = grpname.str();
if (!mapOfJunctionGroups.count(namegrp))
- mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(SMDSAbs_Volume, namegrp.c_str(), idg);
+ mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(vol->GetType(), namegrp.c_str(), idg);
SMESHDS_Group *sgrp = dynamic_cast<SMESHDS_Group*>(mapOfJunctionGroups[namegrp]->GetGroupDS());
if (sgrp)
sgrp->Add(vol->GetID());
}
// --- create volumes on multiple domain intersection if requested
+ // iterate on mutipleNodesToFace
// iterate on edgesMultiDomains
if (createJointElems)
{
+ // --- iterate on mutipleNodesToFace
+
+ std::map<int, std::vector<int> >::iterator itn = mutipleNodesToFace.begin();
+ for (; itn != mutipleNodesToFace.end(); ++itn)
+ {
+ int node = itn->first;
+ vector<int> orderDom = itn->second;
+ vector<vtkIdType> orderedNodes;
+ for (int idom = 0; idom <orderDom.size(); idom++)
+ orderedNodes.push_back( nodeDomains[node][orderDom[idom]] );
+ SMDS_MeshFace* face = this->GetMeshDS()->AddFaceFromVtkIds(orderedNodes);
+
+ stringstream grpname;
+ grpname << "m2j_";
+ grpname << 0 << "_" << 0;
+ int idg;
+ string namegrp = grpname.str();
+ if (!mapOfJunctionGroups.count(namegrp))
+ mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(SMDSAbs_Face, namegrp.c_str(), idg);
+ SMESHDS_Group *sgrp = dynamic_cast<SMESHDS_Group*>(mapOfJunctionGroups[namegrp]->GetGroupDS());
+ if (sgrp)
+ sgrp->Add(face->GetID());
+ }
+
+ // --- iterate on edgesMultiDomains
+
std::map<std::vector<int>, std::vector<int> >::iterator ite = edgesMultiDomains.begin();
for (; ite != edgesMultiDomains.end(); ++ite)
{
}
else
{
- //MESSAGE("Quadratic multiple joints not implemented");
+ MESSAGE("Quadratic multiple joints not implemented");
// TODO quadratic nodes
}
}
