FEM: Add getVolumesByFace and write_face_load functions

getVolumesByFace returns std::map<int, int> with ID of volume
and a number of face as per CalculiX definition. The same function is
accessible for python and returns list with the same information, like
this: [[229, 3], [230, 3], [233, 2], [238, 2]]

write_face_load produces something like this in the .inp file:

***********************************************************
** element + CalculiX face + load in [MPa]
** written by write_face_load function
*DLOAD
** Load on face Face2
229,P3,10.0
230,P3,10.0
233,P2,10.0
238,P2,10.0

Optimised by wmayer

Signed-off-by: wmayer
Signed-off-by: Przemo Firszt <przemo@firszt.eu>

src/Mod/Fem/App/FemMesh.cpp

@@ -404,7 +404,88 @@ std::set<long> FemMesh::getSurfaceNodes(long ElemId, short FaceId, float Angle)
     return result;
 }
 
+std::map<int, int> FemMesh::getVolumesByFace(const TopoDS_Face &face) const
 {
+    std::map<int, int> result;
+    std::set<int> nbf = FemMesh::getNodesByFace(face);
+    std::set<int> nodes_on_face;
+    nodes_on_face.insert(nbf.begin(), nbf.end());
+    nbf.clear();
+
+    static std::map<int, std::vector<int> > elem_order;
+    if (elem_order.empty()) {
+        std::vector<int> c3d4  = boost::assign::list_of(1)(0)(2)(3);
+        std::vector<int> c3d10 = boost::assign::list_of(1)(0)(2)(3)(4)(6)(5)(8)(7)(9);
+
+        elem_order.insert(std::make_pair(c3d4.size(), c3d4));
+        elem_order.insert(std::make_pair(c3d10.size(), c3d10));
+    }
+
+    SMDS_VolumeIteratorPtr vol_iter = myMesh->GetMeshDS()->volumesIterator();
+    std::set<int> element_nodes;
+    int num_of_nodes;
+    while (vol_iter->more()) {
+        const SMDS_MeshVolume* vol = vol_iter->next();
+        num_of_nodes = vol->NbNodes();
+        std::pair<int, std::vector<int> > apair;
+        apair.first = vol->GetID();
+
+        std::map<int, std::vector<int> >::iterator it = elem_order.find(num_of_nodes);
+        if (it != elem_order.end()) {
+            const std::vector<int>& order = it->second;
+            for (std::vector<int>::const_iterator jt = order.begin(); jt != order.end(); ++jt) {
+                int vid = vol->GetNode(*jt)->GetID();
+                apair.second.push_back(vid);
+            }
+        }
+
+        // Get volume nodes on face
+        std::vector<int> element_face_nodes;
+        std::set<int> element_nodes;
+        element_nodes.insert(apair.second.begin(), apair.second.end());
+        std::set_intersection(nodes_on_face.begin(), nodes_on_face.end(), element_nodes.begin(), element_nodes.end(),
+        std::back_insert_iterator<std::vector<int> >(element_face_nodes));
+
+        if ((element_face_nodes.size() == 3 && num_of_nodes == 4) ||
+            (element_face_nodes.size() == 6 && num_of_nodes == 10)) {
+            int missing_node = 0;
+            for (int i=0; i<4; i++) {
+                // search for the ID of the volume which is not part of 'element_face_nodes'
+                if (std::find(element_face_nodes.begin(), element_face_nodes.end(), apair.second[i]) == element_face_nodes.end()) {
+                    missing_node = i + 1;
+                    break;
+                }
+            }
+            /* for tetrahedral elements as per CalculiX definition:
+             Face 1: 1-2-3, missing point 4 means it's face P1
+             Face 2: 1-4-2, missing point 3 means it's face P2
+             Face 3: 2-4-3, missing point 1 means it's face P3
+             Face 4: 3-4-1, missing point 2 means it's face P4 */
+            int face_ccx;
+            switch (missing_node) {
+            case 1:
+                face_ccx = 3;
+                break;
+            case 2:
+                face_ccx = 4;
+                break;
+            case 3:
+                face_ccx = 2;
+                break;
+            case 4:
+                face_ccx = 1;
+                break;
+            default:
+                assert(false); // should never happen
+                break;
+            }
+            result[apair.first] = face_ccx;
+        }
+    }
+
+    return result;
+}
+
 std::set<int> FemMesh::getNodesByFace(const TopoDS_Face &face) const
 {
     std::set<int> result;

src/Mod/Fem/App/FemMesh.h

@@ -93,6 +93,8 @@ class AppFemExport FemMesh : public Data::ComplexGeoData
     std::set<long> getNodesByEdge(const TopoDS_Edge &edge) const;
     /// retrieving by vertex
     std::set<long> getNodesByVertex(const TopoDS_Vertex &vertex) const;
+    /// retrieving volume IDs by face
+    std::map<int, int> getVolumesByFace(const TopoDS_Face &face) const;
     //@}
 
     /** @name Placement control */

src/Mod/Fem/App/FemMeshPy.xml

@@ -84,6 +84,11 @@
               <UserDocu>Make a copy of this FEM mesh.</UserDocu>
           </Documentation>
       </Methode>
+    <Methode Name="getVolumesByFace" Const="true">
+      <Documentation>
+        <UserDocu>Return a list of volume IDs which belong to a TopoFace</UserDocu>
+      </Documentation>
+    </Methode>
     <Methode Name="getNodeById" Const="true">
       <Documentation>
         <UserDocu>Get the node position vector by an Node-ID</UserDocu>

src/Mod/Fem/App/FemMeshPyImp.cpp

@@ -514,6 +514,38 @@ PyObject* FemMeshPy::setTransform(PyObject *args)
     Py_Return;
 }
 
+PyObject* FemMeshPy::getVolumesByFace(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;
+        const TopoDS_Face& fc = TopoDS::Face(sh);
+        if (sh.IsNull()) {
+            PyErr_SetString(Base::BaseExceptionFreeCADError, "Face is empty");
+            return 0;
+        }
+        Py::List ret;
+        std::map<int, int> resultSet = getFemMeshPtr()->getVolumesByFace(fc);
+        for (std::map<int, int>::const_iterator it = resultSet.begin();it!=resultSet.end();++it) {
+		Py::List vol_face;
+                vol_face.append(Py::Int(it->first));
+                vol_face.append(Py::Int(it->second));
+                ret.append(vol_face);
+	}
+
+        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::getNodeById(PyObject *args)
 {
     int id;

src/Mod/Fem/ccxInpWriter.py

@@ -29,6 +29,7 @@ def write_calculix_input_file(self):
         self.write_step_begin(inpfile)
         self.write_constraints_fixed(inpfile)
         self.write_constraints_force(inpfile)
+        self.write_face_load(inpfile)
         self.write_outputs_types(inpfile)
         self.write_step_end(inpfile)
         self.write_footer(inpfile)
@@ -89,10 +90,7 @@ def write_load_node_sets(self, f):
             for o, elem in frc_obj.References:
                 fo = o.Shape.getElement(elem)
                 n = []
-                if fo.ShapeType == 'Face':
-                    print '  AreaLoad (face load) on: ', elem
-                    n = self.mesh_object.FemMesh.getNodesByFace(fo)
-                elif fo.ShapeType == 'Edge':
+                if fo.ShapeType == 'Edge':
                     print '  Line Load (edge load) on: ', elem
                     n = self.mesh_object.FemMesh.getNodesByEdge(fo)
                 elif fo.ShapeType == 'Vertex':
@@ -177,6 +175,21 @@ def write_constraints_force(self, f):
                 f.write(frc_obj_name + ',2,' + v2 + '\n')
                 f.write(frc_obj_name + ',3,' + v3 + '\n\n')
 
+    def write_face_load(self, f):
+        f.write('***********************************************************\n')
+        f.write('** element + CalculiX face + load in [MPa]\n')
+        f.write('** written by {} function\n'.format(sys._getframe().f_code.co_name))
+        for fobj in self.force_objects:
+            frc_obj = fobj['Object']
+            f.write('*DLOAD\n')
+            for o, e in frc_obj.References:
+                elem = o.Shape.getElement(e)
+                if elem.ShapeType == 'Face':
+                    v = self.mesh_object.FemMesh.getVolumesByFace(elem)
+                    f.write("** Load on face {}\n".format(e))
+                    for i in v:
+                        f.write("{},P{},{}\n".format(i[0], i[1], -1 * frc_obj.Force))
+
     def write_outputs_types(self, f):
         f.write('\n** outputs --> frd file\n')
         f.write('** written by {} function\n'.format(sys._getframe().f_code.co_name))