From 48efcc449b87236e8250e42341714379db5f68d0 Mon Sep 17 00:00:00 2001 From: Przemo Firszt Date: Tue, 5 May 2015 20:27:20 +0100 Subject: [PATCH] FEM: Add getccxVolumesByFace and write_face_load functions getccxVolumesByFace returns std::map 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 --- src/Mod/Fem/App/FemMesh.cpp | 82 ++++++++++++++++++++++++++++++++ src/Mod/Fem/App/FemMesh.h | 2 + src/Mod/Fem/App/FemMeshPy.xml | 5 ++ src/Mod/Fem/App/FemMeshPyImp.cpp | 32 +++++++++++++ src/Mod/Fem/ccxInpWriter.py | 21 ++++++-- 5 files changed, 138 insertions(+), 4 deletions(-) diff --git a/src/Mod/Fem/App/FemMesh.cpp b/src/Mod/Fem/App/FemMesh.cpp index b94df6c85..515b7a5e7 100755 --- a/src/Mod/Fem/App/FemMesh.cpp +++ b/src/Mod/Fem/App/FemMesh.cpp @@ -404,7 +404,89 @@ std::set FemMesh::getSurfaceNodes(long ElemId, short FaceId, float Angle) 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. + */ +std::map FemMesh::getccxVolumesByFace(const TopoDS_Face &face) const { + std::map result; + std::set nodes_on_face = FemMesh::getNodesByFace(face); + + static std::map > elem_order; + if (elem_order.empty()) { + std::vector c3d4 = boost::assign::list_of(1)(0)(2)(3); + std::vector 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 element_nodes; + int num_of_nodes; + while (vol_iter->more()) { + const SMDS_MeshVolume* vol = vol_iter->next(); + num_of_nodes = vol->NbNodes(); + std::pair > apair; + apair.first = vol->GetID(); + + std::map >::iterator it = elem_order.find(num_of_nodes); + if (it != elem_order.end()) { + const std::vector& order = it->second; + for (std::vector::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 element_face_nodes; + std::set 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 >(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 FemMesh::getNodesByFace(const TopoDS_Face &face) const { std::set result; diff --git a/src/Mod/Fem/App/FemMesh.h b/src/Mod/Fem/App/FemMesh.h index 7bf8a4051..52f7ffd52 100755 --- a/src/Mod/Fem/App/FemMesh.h +++ b/src/Mod/Fem/App/FemMesh.h @@ -93,6 +93,8 @@ public: std::set getNodesByEdge(const TopoDS_Edge &edge) const; /// retrieving by vertex std::set getNodesByVertex(const TopoDS_Vertex &vertex) const; + /// retrieving volume IDs and CalculiX face number by face + std::map getccxVolumesByFace(const TopoDS_Face &face) const; //@} /** @name Placement control */ diff --git a/src/Mod/Fem/App/FemMeshPy.xml b/src/Mod/Fem/App/FemMeshPy.xml index c0bbb5325..1e2133a7c 100755 --- a/src/Mod/Fem/App/FemMeshPy.xml +++ b/src/Mod/Fem/App/FemMeshPy.xml @@ -84,6 +84,11 @@ Make a copy of this FEM mesh. + + + Return a list of volume IDs which belong to a TopoFace + + Get the node position vector by an Node-ID diff --git a/src/Mod/Fem/App/FemMeshPyImp.cpp b/src/Mod/Fem/App/FemMeshPyImp.cpp index 18a0c0f14..df601a858 100755 --- a/src/Mod/Fem/App/FemMeshPyImp.cpp +++ b/src/Mod/Fem/App/FemMeshPyImp.cpp @@ -514,6 +514,38 @@ PyObject* FemMeshPy::setTransform(PyObject *args) Py_Return; } +PyObject* FemMeshPy::getccxVolumesByFace(PyObject *args) +{ + PyObject *pW; + if (!PyArg_ParseTuple(args, "O!", &(Part::TopoShapeFacePy::Type), &pW)) + return 0; + + try { + const TopoDS_Shape& sh = static_cast(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 resultSet = getFemMeshPtr()->getccxVolumesByFace(fc); + for (std::map::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; diff --git a/src/Mod/Fem/ccxInpWriter.py b/src/Mod/Fem/ccxInpWriter.py index e2ae424f6..d3e9bfdfb 100644 --- a/src/Mod/Fem/ccxInpWriter.py +++ b/src/Mod/Fem/ccxInpWriter.py @@ -29,6 +29,7 @@ class inp_writer: 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 @@ class inp_writer: 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 @@ class inp_writer: 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.getccxVolumesByFace(elem) + f.write("** Load on face {}\n".format(e)) + for i in v: + f.write("{},P{},{}\n".format(i[0], i[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))