From 086e73ab75b5d2e5f85d8dd906e7ea58fb387ac3 Mon Sep 17 00:00:00 2001
From: skl <skl@opencascade.com>
Date: Thu, 26 Jan 2006 07:14:20 +0000
Subject: [PATCH] Changes for working with quadratic elements

---
 src/Controls/SMESH_Controls.cxx | 252 +++++++++++++++++++++++++-------
 src/SMDS/SMDS_Mesh.cxx          |   7 -
 src/SMDS/SMDS_VolumeTool.cxx    |  54 +++++--
 3 files changed, 235 insertions(+), 78 deletions(-)

diff --git a/src/Controls/SMESH_Controls.cxx b/src/Controls/SMESH_Controls.cxx
index 717f8e66a..153a0baee 100644
--- a/src/Controls/SMESH_Controls.cxx
+++ b/src/Controls/SMESH_Controls.cxx
@@ -46,6 +46,7 @@
 #include "SMDS_MeshElement.hxx"
 #include "SMDS_MeshNode.hxx"
 #include "SMDS_VolumeTool.hxx"
+#include "SMDS_QuadraticFaceOfNodes.hxx"
 
 
 /*
@@ -162,14 +163,30 @@ bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem,
     return false;
 
   // Get nodes of the element
-  SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
-  if ( anIter != 0 )
-  {
-    while( anIter->more() )
-    {
-      const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
-      if ( aNode != 0 ){
-        theRes.push_back( gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) );
+
+  const SMDS_QuadraticFaceOfNodes* F =
+    static_cast<const SMDS_QuadraticFaceOfNodes*>(anElem);
+  if(F) {
+    // use special nodes iterator
+    SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
+    if ( anIter != 0 ) {
+      while( anIter->more() ) {
+        const SMDS_MeshNode* aNode = anIter->next();
+        if ( aNode != 0 ) {
+          theRes.push_back( gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) );
+        }
+      }
+    }
+  }
+
+  if(theRes.size()==0) {
+    SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
+    if ( anIter != 0 ) {
+      while( anIter->more() ) {
+        const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
+        if ( aNode != 0 ) {
+          theRes.push_back( gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) );
+        }
       }
     }
   }
@@ -189,10 +206,12 @@ void  NumericalFunctor::SetPrecision( const long thePrecision )
 
 double NumericalFunctor::GetValue( long theId )
 {
+cout<<"theId="<<theId<<endl;
   TSequenceOfXYZ P;
   if ( GetPoints( theId, P ))
   {
     double aVal = GetValue( P );
+cout<<"aVal="<<aVal<<endl;
     if ( myPrecision >= 0 )
     {
       double prec = pow( 10., (double)( myPrecision ) );
@@ -285,6 +304,7 @@ double AspectRatio::GetValue( const TSequenceOfXYZ& P )
   // According to "Mesh quality control" by Nadir Bouhamau referring to
   // Pascal Jean Frey and Paul-Louis George. Maillages, applications aux elements finis.
   // Hermes Science publications, Paris 1999 ISBN 2-7462-0024-4
+  // PAL10872
 
   int nbNodes = P.size();
 
@@ -840,17 +860,16 @@ SMDSAbs_ElementType Skew::GetType() const
 */
 double Area::GetValue( const TSequenceOfXYZ& P )
 {
-  double aArea = 0;
-  if ( P.size() == 3 )
-    return getArea( P( 1 ), P( 2 ), P( 3 ) );
-  else if (P.size() > 3)
-    aArea = getArea( P( 1 ), P( 2 ), P( 3 ) );
-  else
-    return 0;
-
-  for (int i=4; i<=P.size(); i++)
-    aArea += getArea(P(1),P(i-1),P(i));
-  return aArea;
+  gp_Vec aVec1( P(2) - P(1) );
+  gp_Vec aVec2( P(3) - P(1) );
+  gp_Vec SumVec = aVec1 ^ aVec2;
+  for (int i=4; i<=P.size(); i++) {
+    gp_Vec aVec1( P(i-1) - P(1) );
+    gp_Vec aVec2( P(i) - P(1) );
+    gp_Vec tmp = aVec1 ^ aVec2;
+    SumVec.Add(tmp);
+  }
+  return SumVec.Magnitude() * 0.5;
 }
 
 double Area::GetBadRate( double Value, int /*nbNodes*/ ) const
@@ -894,7 +913,10 @@ double Length2D::GetValue( long theElementId)
 {
   TSequenceOfXYZ P;
 
+  //cout<<"Length2D::GetValue"<<endl;
   if (GetPoints(theElementId,P)){
+    //for(int jj=1; jj<=P.size(); jj++)
+    //  cout<<"jj="<<jj<<" P("<<P(jj).X()<<","<<P(jj).Y()<<","<<P(jj).Z()<<")"<<endl;
 
     double  aVal;// = GetValue( P );
     const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId );
@@ -908,7 +930,11 @@ double Length2D::GetValue( long theElementId)
     case SMDSAbs_Edge:
       if (len == 2){
 	aVal = getDistance( P( 1 ), P( 2 ) );
-	break;
+        break;
+      }
+      else if (len == 3){ // quadratic edge
+	aVal = getDistance(P( 1 ),P( 3 )) + getDistance(P( 3 ),P( 2 ));
+        break;
       }
     case SMDSAbs_Face:
       if (len == 3){ // triangles
@@ -926,6 +952,22 @@ double Length2D::GetValue( long theElementId)
 	aVal = Max(Max(L1,L2),Max(L3,L4));
 	break;
       }
+      if (len == 6){ // quadratic triangles
+	double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
+	double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
+	double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 ));
+	aVal = Max(L1,Max(L2,L3));
+        //cout<<"L1="<<L1<<" L2="<<L2<<"L3="<<L3<<" aVal="<<aVal<<endl;
+	break;
+      }
+      else if (len == 8){ // quadratic quadrangles
+	double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
+	double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
+	double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 7 ));
+	double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 ));
+	aVal = Max(Max(L1,L2),Max(L3,L4));
+	break;
+      }
     case SMDSAbs_Volume:
       if (len == 4){ // tetraidrs
 	double L1 = getDistance(P( 1 ),P( 2 ));
@@ -987,6 +1029,63 @@ double Length2D::GetValue( long theElementId)
 
       }
 
+      if (len == 10){ // quadratic tetraidrs
+	double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 ));
+	double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 ));
+	double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
+	double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
+	double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
+	double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
+	aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+	break;
+      }
+      else if (len == 13){ // quadratic piramids
+	double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
+	double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
+	double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 1 ));
+	double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
+	double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
+	double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
+	double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
+	double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
+	aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+	aVal = Max(aVal,Max(L7,L8));
+	break;
+      }
+      else if (len == 15){ // quadratic pentaidres
+	double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
+	double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
+	double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
+	double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
+	double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
+	double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
+	double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
+	double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
+	double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
+	aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+	aVal = Max(aVal,Max(Max(L7,L8),L9));
+	break;
+      }
+      else if (len == 20){ // quadratic hexaider
+	double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 ));
+	double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 ));
+	double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 ));
+	double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
+	double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
+	double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
+	double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
+	double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
+	double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
+	double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
+	double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
+	double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
+	aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+	aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
+	aVal = Max(aVal,Max(L11,L12));
+	break;
+
+      }
+
     default: aVal=-1;
     }
 
@@ -1038,38 +1137,81 @@ void Length2D::GetValues(TValues& theValues){
   SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
   for(; anIter->more(); ){
     const SMDS_MeshFace* anElem = anIter->next();
-    SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
-    long aNodeId[2];
-    gp_Pnt P[3];
 
-    double aLength;
-    const SMDS_MeshElement* aNode;
-    if(aNodesIter->more()){
-      aNode = aNodesIter->next();
-      const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
-      P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
-      aNodeId[0] = aNodeId[1] = aNode->GetID();
-      aLength = 0;
-    }
-    for(; aNodesIter->more(); ){
-      aNode = aNodesIter->next();
-      const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
-      long anId = aNode->GetID();
-
-      P[2] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
+    if(anElem->IsQuadratic()) {
+      const SMDS_QuadraticFaceOfNodes* F =
+        static_cast<const SMDS_QuadraticFaceOfNodes*>(anElem);
+      // use special nodes iterator
+      SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
+      long aNodeId[4];
+      gp_Pnt P[4];
+
+      double aLength;
+      const SMDS_MeshElement* aNode;
+      if(anIter->more()){
+        aNode = anIter->next();
+        const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
+        P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
+        aNodeId[0] = aNodeId[1] = aNode->GetID();
+        aLength = 0;
+      }
+      for(; anIter->more(); ){
+        const SMDS_MeshNode* N1 = static_cast<const SMDS_MeshNode*> (anIter->next());
+        P[2] = gp_Pnt(N1->X(),N1->Y(),N1->Z());
+        aNodeId[2] = N1->GetID();
+        aLength = P[1].Distance(P[2]);
+        if(!anIter->more()) break;
+        const SMDS_MeshNode* N2 = static_cast<const SMDS_MeshNode*> (anIter->next());
+        P[3] = gp_Pnt(N2->X(),N2->Y(),N2->Z());
+        aNodeId[3] = N2->GetID();
+        aLength += P[2].Distance(P[3]);
+        Value aValue1(aLength,aNodeId[1],aNodeId[2]);
+        Value aValue2(aLength,aNodeId[2],aNodeId[3]);
+        P[1] = P[3];
+        aNodeId[1] = aNodeId[3];
+        theValues.insert(aValue1);
+        theValues.insert(aValue2);
+      }
+      aLength += P[2].Distance(P[0]);
+      Value aValue1(aLength,aNodeId[1],aNodeId[2]);
+      Value aValue2(aLength,aNodeId[2],aNodeId[0]);
+      theValues.insert(aValue1);
+      theValues.insert(aValue2);
+    }
+    else {
+      SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
+      long aNodeId[2];
+      gp_Pnt P[3];
+
+      double aLength;
+      const SMDS_MeshElement* aNode;
+      if(aNodesIter->more()){
+        aNode = aNodesIter->next();
+        const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
+        P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
+        aNodeId[0] = aNodeId[1] = aNode->GetID();
+        aLength = 0;
+      }
+      for(; aNodesIter->more(); ){
+        aNode = aNodesIter->next();
+        const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
+        long anId = aNode->GetID();
+        
+        P[2] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
+        
+        aLength = P[1].Distance(P[2]);
+        
+        Value aValue(aLength,aNodeId[1],anId);
+        aNodeId[1] = anId;
+        P[1] = P[2];
+        theValues.insert(aValue);
+      }
 
-      aLength = P[1].Distance(P[2]);
+      aLength = P[0].Distance(P[1]);
 
-      Value aValue(aLength,aNodeId[1],anId);
-      aNodeId[1] = anId;
-      P[1] = P[2];
+      Value aValue(aLength,aNodeId[0],aNodeId[1]);
       theValues.insert(aValue);
     }
-
-    aLength = P[0].Distance(P[1]);
-
-    Value aValue(aLength,aNodeId[0],aNodeId[1]);
-    theValues.insert(aValue);
   }
 }
 
@@ -1207,7 +1349,7 @@ void MultiConnection2D::GetValues(MValues& theValues){
       const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode1;
       aNodeId[0] = aNodeId[1] = aNodes->GetID();
     }
-    for(; aNodesIter->more(); ){
+    for(; aNodesIter->more(); ) {
       aNode2 = (SMDS_MeshNode*) aNodesIter->next();
       long anId = aNode2->GetID();
       aNodeId[2] = anId;
@@ -1217,7 +1359,8 @@ void MultiConnection2D::GetValues(MValues& theValues){
       if (aItr != theValues.end()){
 	aItr->second += 1;
 	//aNbConnects = nb;
-      } else {
+      }
+      else {
 	theValues[aValue] = 1;
 	//aNbConnects = 1;
       }
@@ -1227,10 +1370,11 @@ void MultiConnection2D::GetValues(MValues& theValues){
     }
     Value aValue(aNodeId[0],aNodeId[2]);
     MValues::iterator aItr = theValues.find(aValue);
-    if (aItr != theValues.end()){
+    if (aItr != theValues.end()) {
       aItr->second += 1;
       //aNbConnects = nb;
-    } else {
+    }
+    else {
       theValues[aValue] = 1;
       //aNbConnects = 1;
     }
@@ -2088,8 +2232,7 @@ static gp_XYZ getNormale( const SMDS_MeshFace* theFace )
   TColgp_Array1OfXYZ anArrOfXYZ(1,4);
   SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator();
   int i = 1;
-  for ( ; aNodeItr->more() && i <= 4; i++ )
-  {
+  for ( ; aNodeItr->more() && i <= 4; i++ ) {
     SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
     anArrOfXYZ.SetValue(i, gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) );
   }
@@ -2097,8 +2240,7 @@ static gp_XYZ getNormale( const SMDS_MeshFace* theFace )
   gp_XYZ q1 = anArrOfXYZ.Value(2) - anArrOfXYZ.Value(1);
   gp_XYZ q2 = anArrOfXYZ.Value(3) - anArrOfXYZ.Value(1);
   n  = q1 ^ q2;
-  if ( aNbNode > 3 )
-  {
+  if ( aNbNode > 3 ) {
     gp_XYZ q3 = anArrOfXYZ.Value(4) - anArrOfXYZ.Value(1);
     n += q2 ^ q3;
   }
diff --git a/src/SMDS/SMDS_Mesh.cxx b/src/SMDS/SMDS_Mesh.cxx
index 5bc8cd544..7f4eb3269 100644
--- a/src/SMDS/SMDS_Mesh.cxx
+++ b/src/SMDS/SMDS_Mesh.cxx
@@ -2248,7 +2248,6 @@ SMDS_MeshEdge* SMDS_Mesh::AddEdgeWithID(const SMDS_MeshNode * n1,
 					const SMDS_MeshNode * n12, 
 					int ID)
 {
-  cout<<"SMDS_Mesh::AddEdgeWithID"<<endl;
   SMDS_QuadraticEdge* edge = new SMDS_QuadraticEdge(n1,n2,n12);
   if(myElementIDFactory->BindID(ID, edge)) {
     SMDS_MeshNode *node1,*node2, *node12;
@@ -2313,7 +2312,6 @@ SMDS_MeshFace* SMDS_Mesh::AddFaceWithID(const SMDS_MeshNode * n1,
                                         const SMDS_MeshNode * n31, 
                                         int ID)
 {
-  cout<<"SMDS_Mesh::AddFaceWithID"<<endl;
   if(hasConstructionEdges()) {
     // creation quadratic edges - not implemented
   }
@@ -2381,7 +2379,6 @@ SMDS_MeshFace* SMDS_Mesh::AddFaceWithID(const SMDS_MeshNode * n1,
                                         const SMDS_MeshNode * n41, 
                                         int ID)
 {
-  cout<<"SMDS_Mesh::AddFaceWithID"<<endl;
   if(hasConstructionEdges()) {
     // creation quadratic edges - not implemented
   }
@@ -2461,7 +2458,6 @@ SMDS_MeshVolume* SMDS_Mesh::AddVolumeWithID(const SMDS_MeshNode * n1,
                                             const SMDS_MeshNode * n34,
                                             int ID)
 {
-  cout<<"SMDS_Mesh::AddVolumeWithID"<<endl;
   if(hasConstructionFaces()) {
     // creation quadratic faces - not implemented
   }
@@ -2554,7 +2550,6 @@ SMDS_MeshVolume* SMDS_Mesh::AddVolumeWithID(const SMDS_MeshNode * n1,
                                             const SMDS_MeshNode * n45,
                                             int ID)
 {
-  cout<<"SMDS_Mesh::AddVolumeWithID"<<endl;
   if(hasConstructionFaces()) {
     // creation quadratic faces - not implemented
   }
@@ -2657,7 +2652,6 @@ SMDS_MeshVolume* SMDS_Mesh::AddVolumeWithID(const SMDS_MeshNode * n1,
                                             const SMDS_MeshNode * n36,
                                             int ID)
 {
-  cout<<"SMDS_Mesh::AddVolumeWithID"<<endl;
   if(hasConstructionFaces()) {
     // creation quadratic faces - not implemented
   }
@@ -2780,7 +2774,6 @@ SMDS_MeshVolume* SMDS_Mesh::AddVolumeWithID(const SMDS_MeshNode * n1,
                                             const SMDS_MeshNode * n48,
                                             int ID)
 {
-  cout<<"SMDS_Mesh::AddVolumeWithID"<<endl;
   if(hasConstructionFaces()) {
     // creation quadratic faces - not implemented
   }
diff --git a/src/SMDS/SMDS_VolumeTool.cxx b/src/SMDS/SMDS_VolumeTool.cxx
index 0fa7d2f0f..817d3ceb4 100644
--- a/src/SMDS/SMDS_VolumeTool.cxx
+++ b/src/SMDS/SMDS_VolumeTool.cxx
@@ -713,7 +713,7 @@ double SMDS_VolumeTool::GetSize() const
   else 
   {
     const static int ind[] = {
-      0, 1, 3, 6, 11, 19, 32, 35, 40};
+      0, 1, 3, 6, 11, 19, 32, 46, 66};
     const static int vtab[][4] = {
       // tetrahedron
       { 0, 1, 2, 3 },
@@ -740,10 +740,6 @@ double SMDS_VolumeTool::GetSize() const
       { 4, 5, 6, 9 },
       { 4, 7, 8, 9 },
       { 4, 5, 9, 8 },
-      //{ 6, 7, 9, 4 },
-      //{ 5, 6, 9, 4 },
-      //{ 7, 8, 9, 4 },
-      //{ 5, 9, 8, 4 },
 
       // quadratic pyramid
       { 0, 5, 8, 9 },
@@ -761,31 +757,57 @@ double SMDS_VolumeTool::GetSize() const
       { 10, 6, 8, 7 },
 
       // quadratic pentahedron
-      { 0, 1, 2, 3 },
-      { 1, 5, 3, 4 },
-      { 1, 5, 2, 3 },
+      { 12, 0, 8, 6 },
+      { 12, 8, 7, 6 },
+      { 12, 8, 2, 7 },
+      { 12, 6, 7, 1 },
+      { 12, 1, 7,13 },
+      { 12, 7, 2,13 },
+      { 12, 2,14,13 },
+
+      { 12, 3, 9,11 },
+      { 12,11, 9,10 },
+      { 12,11,10, 5 },
+      { 12, 9, 4,10 },
+      { 12,14, 5,10 },
+      { 12,14,10, 4 },
+      { 12,14, 4,13 },
 
       // quadratic hexahedron
-      { 1, 4, 3, 0 },
-      { 4, 1, 6, 5 },
-      { 1, 3, 6, 2 },
-      { 4, 6, 3, 7 },
-      { 1, 4, 6, 3 }
+      { 16, 0,11, 8 },
+      { 16,11, 9, 8 },
+      { 16, 8, 9, 1 },
+      { 16,11, 3,10 },
+      { 16,11,10, 9 },
+      { 16,10, 2, 9 },
+      { 16, 3,19, 2 },
+      { 16, 2,19,18 },
+      { 16, 2,18,17 },
+      { 16, 2,17, 1 },
+
+      { 16, 4,12,15 },
+      { 16,12, 5,13 },
+      { 16,12,13,15 },
+      { 16,13, 6,14 },
+      { 16,13,14,15 },
+      { 16,14, 7,15 },
+      { 16, 6, 5,17 },
+      { 16,18, 6,17 },
+      { 16,18, 7, 6 },
+      { 16,18,19, 7 },
+
     };
 
     int type = GetVolumeType();
     int n1 = ind[type];
     int n2 = ind[type+1];
 
-//cout<<"n1="<<n1<<" n2="<<n2<<endl;
     for (int i = n1; i <  n2; i++) {
-//cout<<"V="<<V<<endl;
       V -= getTetraVolume( myVolumeNodes[ vtab[i][0] ],
                            myVolumeNodes[ vtab[i][1] ],
                            myVolumeNodes[ vtab[i][2] ],
                            myVolumeNodes[ vtab[i][3] ]);
     }
-//cout<<"V="<<V<<endl;
   }
   return V;
 }
-- 
2.39.2