+ 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;
+ }