FEM: Add faces by face workaround for force constraint

src/Mod/Fem/App/FemMesh.cpp

@@ -126,9 +126,7 @@ void FemMesh::copyMeshData(const FemMesh& mesh)
     SMDS_NodeIteratorPtr aNodeIter = mesh.myMesh->GetMeshDS()->nodesIterator();
     for (;aNodeIter->more();) {
         const SMDS_MeshNode* aNode = aNodeIter->next();
-        int id = aNode->GetID();
         double temp[3];
-        Base::Console().Message("CopyData Mesh ID: %i\n", aNode->getMeshId());
         aNode->GetXYZ(temp);
         meshds->AddNodeWithID(temp[0],temp[1],temp[2], aNode->GetID());
     }
@@ -399,6 +397,7 @@ std::set<long> FemMesh::getSurfaceNodes(long ElemId, short FaceId, float Angle)
  */
 std::list<std::pair<int, int> > FemMesh::getVolumesByFace(const TopoDS_Face &face) const
 {
+    //TODO: This function is broken with SMESH7 as it is impossible to iterate volume faces
     std::list<std::pair<int, int> > result;
     std::set<int> nodes_on_face = getNodesByFace(face);
 
@@ -431,6 +430,38 @@ std::list<std::pair<int, int> > FemMesh::getVolumesByFace(const TopoDS_Face &fac
     return result;
 }
 
+/*! That function returns a list of face IDs.
+ */
+std::list<int> FemMesh::getFacesByFace(const TopoDS_Face &face) const
+{
+    //TODO: This function is broken with SMESH7 as it is impossible to iterate volume faces
+    std::list<int> result;
+    std::set<int> nodes_on_face = getNodesByFace(face);
+
+    SMDS_FaceIteratorPtr face_iter = myMesh->GetMeshDS()->facesIterator();
+    while (face_iter->more()) {
+        const SMDS_MeshFace* face = static_cast<const SMDS_MeshFace*>(face_iter->next());
+        int numNodes = face->NbNodes();
+
+        std::set<int> face_nodes;
+        for (int i=0; i<numNodes; i++) {
+            face_nodes.insert(face->GetNode(i)->GetID());
+        }
+
+        std::vector<int> element_face_nodes;
+        std::set_intersection(nodes_on_face.begin(), nodes_on_face.end(), face_nodes.begin(), face_nodes.end(),
+            std::back_insert_iterator<std::vector<int> >(element_face_nodes));
+
+        // For curved faces it is possible that a volume contributes more than one face
+        if (element_face_nodes.size() == static_cast<std::size_t>(numNodes)) {
+            result.push_back(face->GetID());
+        }        
+    }
+
+    result.sort();
+    return result;
+}
+
 /*! That function returns map containing volume ID and face number
  * as per CalculiX definition for tetrahedral elements. See CalculiX
  * documentation for the details.

src/Mod/Fem/App/FemMesh.h

@@ -98,6 +98,8 @@ public:
     std::set<int> getNodesByVertex(const TopoDS_Vertex &vertex) const;
     /// retrieving node IDs by element ID
     std::list<int> getElementNodes(int id) const;
+    /// retrieving face IDs number by face
+    std::list<int> getFacesByFace(const TopoDS_Face &face) const;
     /// retrieving volume IDs and face IDs number by face
     std::list<std::pair<int, int> > getVolumesByFace(const TopoDS_Face &face) const;
     /// retrieving volume IDs and CalculiX face number by face

src/Mod/Fem/App/FemMeshPy.xml

@@ -84,6 +84,11 @@
                 <UserDocu>Make a copy of this FEM mesh.</UserDocu>
             </Documentation>
         </Methode>
+        <Methode Name="getFacesByFace" Const="true">
+            <Documentation>
+                <UserDocu>Return a list of face IDs which belong to a TopoFace</UserDocu>
+            </Documentation>
+        </Methode>
         <Methode Name="getVolumesByFace" Const="true">
             <Documentation>
                 <UserDocu>Return a dict of volume IDs and face IDs which belong to a TopoFace</UserDocu>

src/Mod/Fem/App/FemMeshPyImp.cpp

@@ -575,6 +575,37 @@ PyObject* FemMeshPy::setTransform(PyObject *args)
     Py_Return;
 }
 
+
+PyObject* FemMeshPy::getFacesByFace(PyObject *args)
+{
+    PyObject *pW;
+    if (!PyArg_ParseTuple(args, "O!", &(Part::TopoShapeFacePy::Type), &pW))
+         return 0;
+
+    try {
+        const TopoDS_Shape& sh = static_cast<Part::TopoShapeFacePy*>(pW)->getTopoShapePtr()->_Shape;
+        if (sh.IsNull()) {
+            PyErr_SetString(Base::BaseExceptionFreeCADError, "Face is empty");
+            return 0;
+        }
+
+        const TopoDS_Face& fc = TopoDS::Face(sh);
+
+        Py::List ret;
+        std::list<int> resultSet = getFemMeshPtr()->getFacesByFace(fc);
+        for (std::list<int>::const_iterator it = resultSet.begin();it!=resultSet.end();++it) {
+            ret.append(Py::Int(*it));
+        }
+
+        return Py::new_reference_to(ret);
+    }
+    catch (Standard_Failure) {
+        Handle_Standard_Failure e = Standard_Failure::Caught();
+        PyErr_SetString(Base::BaseExceptionFreeCADError, e->GetMessageString());
+        return 0;
+    }
+}
+
 PyObject* FemMeshPy::getVolumesByFace(PyObject *args)
 {
     PyObject *pW;

src/Mod/Fem/App/FemMeshShapeNetgenObject.cpp

@@ -76,37 +76,13 @@ App::DocumentObjectExecReturn *FemMeshShapeNetgenObject::execute(void)
 {
 #ifdef FCWithNetgen
 
-    Fem::FemMesh newMesh, mesh2;
+    Fem::FemMesh newMesh;
-    
-    Base::Console().Message("newMesh ID: %i\n", newMesh.getSMesh()->GetMeshDS()->getMeshId());
-    Base::Console().Message("second mesh ID: %i\n", mesh2.getSMesh()->GetMeshDS()->getMeshId());
 
     Part::Feature *feat = Shape.getValue<Part::Feature*>();
     TopoDS_Shape shape = feat->Shape.getValue();
 
-
-//    newMesh.myMesh->ShapeToMesh(shape);
-//    SMESH_Gen *myGen = newMesh.getGenerator();
-//    meshgen->CreateMesh(0, true);
-
-    int hyp=0;
-
     NETGENPlugin_Mesher myNetGenMesher(newMesh.getSMesh(),shape,true);
-/*
+    NETGENPlugin_Hypothesis* tet= new NETGENPlugin_Hypothesis(0,1,newMesh.getGenerator());
-    NETGENPlugin_SimpleHypothesis_2D * tet2 = new NETGENPlugin_SimpleHypothesis_2D(hyp++,1,myGen);
-    static_cast<NETGENPlugin_SimpleHypothesis_2D*>(tet2.get())->SetNumberOfSegments(5);
-    static_cast<NETGENPlugin_SimpleHypothesis_2D*>(tet2.get())->SetLocalLength(0.1);
-    static_cast<NETGENPlugin_SimpleHypothesis_2D*>(tet2.get())->LengthFromEdges();
-    myNetGenMesher.SetParameters(tet2);
-*/
-/*
-    NETGENPlugin_SimpleHypothesis_3D* tet= new NETGENPlugin_SimpleHypothesis_3D(hyp++,1,myGen);
-    static_cast<NETGENPlugin_SimpleHypothesis_3D*>(tet.get())->LengthFromFaces();
-    static_cast<NETGENPlugin_SimpleHypothesis_3D*>(tet.get())->SetMaxElementVolume(0.1);
-    myNetGenMesher.SetParameters( tet);
-*/
-
-    NETGENPlugin_Hypothesis* tet= new NETGENPlugin_Hypothesis(hyp++,1,newMesh.getGenerator());
     tet->SetMaxSize(MaxSize.getValue());
     tet->SetSecondOrder(SecondOrder.getValue());
     tet->SetOptimize(Optimize.getValue());

src/Mod/Fem/FemMeshTools.py

@@ -426,8 +426,8 @@ def get_ref_facenodes_table(femmesh, femelement_table, ref_face):
                         ve_ref_face_nodes.append(nodeID)
                 face_table[veID] = ve_ref_face_nodes  # { volumeID : ( facenodeID, ... , facenodeID ) } only the ref_face nodes
         else:  # the femmesh has face_data
-            volume_faces = femmesh.getVolumesByFace(ref_face)   # (mv, mf)
+            faces = femmesh.getFacesByFace(ref_face)   # (mv, mf)
-            for mv, mf in volume_faces:
+            for mf in faces:
                 face_table[mf] = femmesh.getElementNodes(mf)
     elif is_face_femmesh(femmesh):
         ref_face_nodes = femmesh.getNodesByFace(ref_face)